Gewöhnliches Kriging-Beispiel Schritt für Schritt?

8

Ich habe Online-Tutorials für räumliches Kriging mit beiden geoRund gstat(und auch automap) verfolgt. Ich kann räumliches Kriging durchführen und verstehe die Hauptkonzepte dahinter. Ich weiß, wie man ein Semivariogramm erstellt, wie man ein Modell daran anpasst und wie man gewöhnliches Kriging durchführt.

Was ich nicht verstehe ist, wie die Gewichte der umgebenden Messwerte bestimmt werden. Ich weiß, dass sie aus dem Semivariogramm stammen und von der Entfernung vom Vorhersageort und von der räumlichen Anordnung der gemessenen Punkte abhängen. Aber wie?

Könnte jemand bitte ein gewöhnliches Kriging-Modell (nicht bayesianisch) mit 3 zufällig gemessenen Punkten und 1 Vorhersageort erstellen? Es wäre aufschlussreich.

Pigna
quelle
1
Warum willst du nicht aus Neugier die Bayes'sche Antwort sehen? Es macht die Dinge viel einfacher, wenn Sie mit Gaußschen Prozessen umgehen.
DeltaIV
@ DeltaIV, weil ich zuerst den frequentistischen Weg lernen möchte. Bayesian Statistiken sind immer noch bewölkt für mich
Pigna
1
" Was ich nicht verstehe, ist, wie die Gewichte der umgebenden Messwerte bestimmt werden. " Falls jemand interessiert ist, habe ich in GIS SE eine Antwort mit einem Beispiel zur Berechnung veröffentlicht ( gis.stackexchange.com/questions/270274/… ). Aber die Antwort hier ist schon großartig!
Andre Silva

Antworten:

9

Zuerst beschreibe ich gewöhnliches Kriging mit drei Punkten mathematisch. Angenommen, wir haben ein intrinsisch stationäres Zufallsfeld.

Gewöhnliches Kriging

Wir versuchen , den Wert vorherzusagen unter Verwendung der bekannten Werte Z = ( Z ( x 1 ) , Z ( x 2 ) , Z ( x 3 ) ) Die Vorhersage wollen wir von der Form Z ( x 0 ) = λ T Z wobei λ = ( λ 1 , λ 2 , λ 3 )Z(x0)Z=(Z(x1),Z(x2),Z(x3))

Z^(x0)=λTZ
λ=(λ1,λ2,λ3)sind die Interpolationsgewichte. Wir nehmen einen konstanten Mittelwert . Um ein unverzerrtes Ergebnis zu erhalten, setzen wir λ 1 + λ 2 + λ 3 = 1 . Wir erhalten dann folgendes Problem: minμλ1+λ2+λ3=1 Unter Verwendung der Lagrange-Multiplikatormethode erhalten wir die Gleichungen: 3 j = 1 λ j γ ( x i - x j ) + m = γ ( x i - x 0 ) ,
minE(Z(X0)λTZ)2s.t.λT1=1.
3 j = 1 λ j = 1 , wobei m der Lagrange-Multiplikator und γ das (Halb-) Variogramm ist. Daraus können wir einige Dinge beobachten:
j=13λjγ(xixj)+m=γ(xix0),i=1,2,3,
j=13λj=1,
mγ
  • Die Gewichte hängen nicht vom Mittelwert .μ
  • Die Gewichte hängen überhaupt nicht von den Werten von ab. Nur an den Koordinaten (im isotropen Fall nur an der Entfernung)Z
  • Jedes Gewicht hängt von der Position aller anderen Punkte ab.

Das genaue Verhalten der Gewichte ist nur aus der Gleichung schwer zu erkennen, aber man kann sehr grob sagen:

  • x0
  • Die Nähe zu anderen Punkten verringert jedoch auch das Gewicht.
  • R

[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)
}

Sie können mit der clickpppFunktion von spatstat damit spielen :

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

Hier einige Beispiele

x0

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)

Geben Sie hier die Bildbeschreibung ein

Punkte nahe beieinander teilen sich die Gewichte

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)

Geben Sie hier die Bildbeschreibung ein

In der Nähe Punkt "stehlen" die Gewichte

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

Es ist möglich, negative Gewichte zu erhalten

Geben Sie hier die Bildbeschreibung ein

Ich hoffe, dies gibt Ihnen ein Gefühl dafür, wie die Gewichte funktionieren.

Dahn
quelle