Exemple de krigeage ordinaire étape par étape?

8

J'ai suivi des tutoriels en ligne pour le krigeage spatial avec geoRet gstat(et aussi automap). Je peux effectuer le krigeage spatial et je comprends les principaux concepts qui le sous-tendent. Je sais comment construire un semi-variogramme, comment y adapter un modèle et comment effectuer un krigeage ordinaire.

Ce que je ne comprends pas, c'est comment les poids des valeurs mesurées environnantes sont déterminés. Je sais qu'ils dérivent du semi-variogramme et dépendent de la distance de l'emplacement de prédiction et de la disposition spatiale des points mesurés. Mais comment?

Quelqu'un pourrait-il faire un modèle de krigeage ordinaire (non bayésien) avec 3 points mesurés au hasard et 1 emplacement de prédiction? Ce serait instructif.

Pigna
la source
1
juste pour la curiosité, pourquoi ne voulez-vous pas voir la réponse bayésienne? Cela rend les choses beaucoup plus simples lorsque vous traitez avec des processus gaussiens.
DeltaIV
@DeltaIV parce que je veux d'abord apprendre la voie fréquentiste. Les statistiques bayésiennes sont encore nuageuses pour moi
Pigna
1
" Ce que je ne comprends pas, c'est comment les poids des valeurs mesurées environnantes sont déterminés. ". Si quelqu'un est intéressé, j'ai publié une réponse dans GIS SE avec un exemple sur la façon de les calculer ( gis.stackexchange.com/questions/270274/… ). Mais la réponse est déjà excellente ici!
Andre Silva

Réponses:

9

Je vais d'abord décrire mathématiquement le krigeage ordinaire avec trois points. Supposons que nous ayons un champ aléatoire intrinsèquement stationnaire.

Krigeage ordinaire

Nous essayons de prédire la valeur utilisant les valeurs connues La prédiction que nous voulons est de la forme où sont les poids d'interpolation. Nous supposons une valeur moyenne constante . Afin d'obtenir un résultat non biaisé, nous fixons . On obtient alors le problème suivant: En utilisant la méthode du multiplicateur de Lagrange, nous obtenons les équations: Z(x0)Z=(Z(X1),Z(X2),Z(X3))

Z^(X0)=λTZ
λ=(λ1,λ2,λ3)μλ1+λ2+λ3=1
minE(Z(X0)-λTZ)2stλT1=1.
j=13λjγ(Xje-Xj)+m=γ(Xje-X0),je=1,2,3,
j=13λj=1,
où est le multiplicateur de décalage et est le (semi) variogramme. De cela, nous pouvons observer quelques choses:mγ
  • Les poids ne dépendent pas de la valeur moyenne .μ
  • Les poids ne dépendent pas des valeurs de à tous. Uniquement sur les coordonnées (dans le cas isotrope sur la distance uniquement)Z
  • Chaque poids dépend de l'emplacement de tous les autres points.

Le comportement précis des poids est difficile à voir uniquement à partir de l'équation, mais on peut dire très grossièrement :

  • Plus le point est de , plus son poids est faible ("plus loin" par rapport aux autres points).X0
  • Cependant, être proche d'autres points diminue également le poids.
  • Le résultat dépend beaucoup de la forme, de la plage et, en particulier, de l'effet pépite du variogramme. Il serait assez éclairant de considérer le krigeage sur avec seulement deux points et de voir comment le résultat change avec différents paramètres de variogramme.R

Je vais cependant me concentrer sur l'emplacement des points dans un avion. J'ai écrit cette petite fonction R qui prend des points de et trace les poids de krigeage (pour la fonction de covariance exponentielle avec zéro pépite).[0,1]2

library(geoR)

# Plots prediction weights for kriging in the window [0,1]x[0,1] with the prediction point (0.5,0.5)
drawWeights <- function(x,y){
 df <- data.frame(x=x,y=y, values = rep(1,length(x)))
  data <- as.geodata(df, coords.col = 1:2, data.col = 3)

  wls <- variofit(bin1,ini=c(1,0.5),fix.nugget=T)
  weights <- round(as.numeric(krweights(data$coords,c(0.5,0.5),krige.control(obj.mod=wls, type="ok"))),3)

  plot(data$coords, xlim=c(0,1),  ylim=c(0,1))
  segments(rep(0.5,length(x)), rep(0.5,length(x)),x, y, lty=3 )
  text((x+0.5)/2,(y+0.5)/2,labels=weights)
}

Vous pouvez jouer avec en utilisant la clickpppfonction de spatstat :

library(spatstat)
points <- clickppp()
drawWeights(points$x,points$y)

Voici quelques exemples

Points équidistants de et les uns des autresX0

deg <- seq(0,2*pi,length.out=4)
deg <- head(deg,length(deg)-1)
x <- 0.5*as.numeric(lapply(deg, cos)) + 0.5
y <- 0.5*as.numeric(lapply(deg, sin)) + 0.5
drawWeights(x,y)

entrez la description de l'image ici

Les points proches les uns des autres partageront les poids

deg <- c(0,0.1,pi)
x <- 0.5*as.numeric(lapply(deg, cos)) + 0.5
y <- 0.5*as.numeric(lapply(deg, sin)) + 0.5
drawWeights(x,y)

entrez la description de l'image ici

Point à proximité "voler" les poids

deg <- seq(0,2*pi,length.out=4)
deg <- head(deg,length(deg)-1)
x <- c(0.6,0.5*as.numeric(lapply(deg, cos)) + 0.5)
y <- c(0.6,0.5*as.numeric(lapply(deg, sin)) + 0.5)
drawWeights(x,y)

https://i.imgur.com/MeuPvFT.png

Il est possible d'obtenir des poids négatifs

entrez la description de l'image ici

J'espère que cela vous donne une idée du fonctionnement des poids.

Dahn
la source