Exemplo de krigagem comum passo a passo?

8

Segui tutoriais on-line para krigagem espacial com ambos geoRe gstat(e também automap). Posso executar krigagem espacial e compreendo os principais conceitos por trás disso. Eu sei como construir um semivariograma, como ajustar um modelo a ele e como executar krigagem comum.

O que não entendo é como são determinados os pesos dos valores medidos circundantes. Eu sei que eles derivam do semivariograma e dependem da distância do local da previsão e do arranjo espacial dos pontos medidos. Mas como?

Alguém poderia, por favor, criar um modelo comum de krigagem (não bayesiana) com 3 pontos aleatórios medidos e 1 local de previsão? Seria esclarecedor.

Pigna
fonte
1
só por curiosidade, por que você não quer ver a resposta bayesiana? Isso torna as coisas muito mais simples quando você lida com processos gaussianos.
DeltaIV
@ DeltaIV, porque primeiro eu quero aprender da maneira frequentista. Estatística Bayesiana ainda estão nublado para mim
Pigna
1
" O que não entendo é como são determinados os pesos dos valores medidos circundantes. ". Caso alguém esteja interessado, postei uma resposta no GIS SE com um exemplo de como calculá-las ( gis.stackexchange.com/questions/270274/… ). Mas a resposta aqui já é ótima!
Andre Silva

Respostas:

9

Primeiro, descreverei a krigagem comum com três pontos matematicamente. Suponha que temos um campo aleatório intrinsecamente estacionário.

Kriging comum

Estamos tentando prever o valor usando os valores conhecidos Z = ( Z ( x 1 ) , Z ( x 2 ) , Z ( x 3 ) ) A previsão queremos é da forma Z ( x 0 ) = λ T Z onde λ = ( λ 1 , λ 2 , λ 3 )Z(x0 0)Z=(Z(x1),Z(x2),Z(x3))

Z^(x0 0)=λTZ
λ=(λ1,λ2,λ3)são os pesos de interpolação. Assumimos um valor médio constante . Para obter um resultado imparcial, fixamos λ 1 + λ 2 + λ 3 = 1 . Em seguida, obtemos o seguinte problema: minμλ1+λ2+λ3=1 Utilizando o método multiplicador de Lagrange, obtemos as equações: 3 j = 1 λ j γ ( x i - x j ) + m = γ ( x i - x 0 ) ,
minE(Z(X0 0)-λTZ)2stλT1=1
3 j = 1 λ j = 1 , onde m é o multiplicador de intervalo e γ é o (semi) variograma. A partir disso, podemos observar algumas coisas:
j=13λjγ(xEu-xj)+m=γ(xEu-x0 0),Eu=1,2,3,
j=13λj=1,
mγ
  • Os pesos não dependem do valor médio .μ
  • Os pesos não depende dos valores de em tudo. Somente nas coordenadas (apenas no caso isotrópico à distância)Z
  • Cada peso depende da localização de todos os outros pontos.

O comportamento preciso dos pesos é difícil de ver apenas a partir da equação, mas pode -se dizer muito grosso modo :

  • x0 0
  • No entanto, estar perto de outros pontos também reduz o peso.
  • R

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

Você pode brincar usando a clickpppfunção spatstat :

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

Aqui estão alguns exemplos

x0 0

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)

insira a descrição da imagem aqui

Pontos próximos entre si compartilharão os pesos

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)

insira a descrição da imagem aqui

Ponto próximo "roubando" os pesos

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

É possível obter pesos negativos

insira a descrição da imagem aqui

Espero que isso lhe dê uma idéia de como os pesos funcionam.

Dahn
fonte