一步一步的普通克里金例子？
我已經按照在線教程學習了空間克里金法
geoR和gstat（以及automap）。我可以執行空間克里金法，並且了解其背後的主要概念。我知道如何構建半變異函數、如何擬合模型以及如何執行普通克里金法。我不明白的是如何確定周圍測量值的權重。我知道它們來自半變異函數，並且取決於與預測位置的距離和測量點的空間排列。但是怎麼做？
任何人都可以製作一個具有 3 個隨機測量點和 1 個預測位置的普通克里金（非貝葉斯）模型嗎？這會很有啟發性。
除了這個答案之外，對於gis.stackexchange.com上的類似問題，還有一些不錯的附加答案
首先，我將在數學上用三點描述普通克里金法。假設我們有一個本質上靜止的隨機場。
普通克里金法
我們試圖預測價值 $ Z(x_0) $ 使用已知值 $ Z=(Z(x_1),Z(x_2),Z(x_3)) $ 我們想要的預測是形式 $$ \hat Z(x_0) = \lambda^T Z $$ 在哪裡 $ \lambda = (\lambda_1,\lambda_2,\lambda_3) $ 是插值權重。我們假設一個恆定的平均值 $ \mu $ . 為了獲得無偏的結果，我們修正 $ \lambda_1 + \lambda_2 + \lambda_3 = 1 $ . 然後我們得到以下問題： $$ \text{min} ; E(Z(X_0) - \lambda^T Z)^2 \quad \text{s.t.};; \lambda^T \mathbf{1} = 1. $$ 使用拉格朗日乘子法，我們得到方程： $$ \sum^3_{j=1} \lambda_j \gamma(x_i - x_j) + m = \gamma(x_i - x_0),;; i=1,2,3, $$ $$ \sum^3_{j=1} \lambda_j =1 , $$ 在哪裡 $ m $ 是拉格朗日乘數和 $ \gamma $ 是（半）變異函數。由此，我們可以觀察到幾件事：
- 權重不依賴於平均值 $ \mu $ .
- 權重不取決於 $ Z $ 一點也不。僅在坐標上（在各向同性情況下僅在距離上）
- 每個權重取決於所有其他點的位置。
僅從方程式很難看出權重的精確行為，但可以非常粗略地說：
- 點離得越遠 $ x_0 $ ，其權重越低（相對於其他點“更遠”）。
- 然而，靠近其他點也會降低重量。
- 結果很大程度上取決於變差函數的形狀、範圍，尤其是塊金效應。考慮克里金法會很有啟發性 $ \mathbb R $ 只有兩個點，看看結果如何隨著不同的變異函數設置而變化。
然而，我將專注於平面中點的位置。我寫了這個小 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) }您可以使用 spatstat 的功能來玩它
clickppp：library(spatstat) points <- clickppp() drawWeights(points$x,points$y)這裡有幾個例子
等距的點 $ x_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)
彼此靠近的點將共享權重
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)
附近點“竊取”權重
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)
有可能得到負權重
希望這能讓您對權重的工作方式有所了解。
