克里金插值的基本介紹可以參考ARCGIS的幫助文檔[1]. 其本質就是根據已知點的數值，確定其周圍點（預測點）的數值。最直觀的方法就是找到已知點和預測點數值之間的關系，從而預測出預測點的數值。比如IDW插值方法，就是假設已知點和預測點的值跟它們相對距離成反比。克里金插值的精妙之處在於它不僅考慮了已知點和預測點的距離關系，還考慮了這些已知點之間的自相關關系。
如何衡量已知點之間的自相關關系呢？通常使用的就是半變異函數，其公式如下[1]：
Semivariogram(distance h) = 0.5 * average((value i – value j)2)
這就是克里金插值不同於其他插值方法的核心所在，通過計算每個距離范圍內所有配對點的離差平方和，就可以繪制出不同距離范圍下的變差值圖，如下示意圖[1]。我們知道，離差平方和是衡量一組數據變化程度的量，而這里通過計算所有已知點與其不同距離范圍內鄰居點集的距離離差平方和，就可以大致衡量出已知點與不同距離范圍鄰居點的變化程度關系，通常的做法是進行曲線擬合，常用的包括指數擬合、球面擬合、高斯擬合等。而這個擬合的結果就是我們插值所需要的塊金值、基台值等。
注意，目前為止這里只涉及到克里金插值是如何分析已知點信息，以及如何構造這些已知點之間的關系，整個插值（預測）過程就是假設這種已知點之間的關系：距離關系和自相關關系，對於預測點同樣適用！
然而，在沒有更多信息的前提下，我們如何知道這種插值是否可信？目前較為合理的方法就是交叉驗證Cross validation。其本質是拿出一些已知點作為預測點，這些被拿出的點不參與上述已知點關系的探索過程，而是作為驗證數據來衡量我們預測是否合理。比如，我們每次拿出一個已知點作為驗證數據，來驗證這個點的預測值，我們就可以得到所有已知點與其預測值之間的偏差，這個所有點的偏差從某種程度上講就為我們提供了整個預測方法是否合理的依據。
一個基於R gstat空間插值包的示例：
如下圖，假設我們已知的數據點的柵格圖，需要插值出其它預測點的值（白色透明區域）：
首先，我們需要將這個柵格數據讀入到R中(需要安裝raster包)
install.packages("raster")
library(raster)
data.observed <- raster("C:/Users/workspace/allmax.img")
因為gstat都是以data frame的格式進行數據處理，首先我們編寫一個函數，將這個柵格數據轉化為data frame：
raster_to_dataframe<- function(inRaster, name) #
{
rp<- rasterToPoints(inRaster)
df = data.frame(rp)
colnames(df) = c("X", "Y", "VALUE")
df$NAME<- name
return(df)
}
name參數可以是隨意的字符串，方便與其他data frame合並的時候辨別數據，這時候我們調用此函數，如下：
data.observed <- raster_to_dataframe (data.observed, "Observed")
查看轉換為data frame的結果，可以使用：
head(data.observed)
X Y VALUE NAME 1 -318466.5 3794841 12 Observed 2 -304466.5 3794841 10 Observed 3 -303466.5 3794841 9 Observed 4 -302466.5 3794841 11 Observed 5 -301466.5 3794841 7 Observed 6 -297466.5 3794841 14 Observed
好了，接下來就是要分析這寫已知點，gstat提供了半變異圖繪制函數variogram:
v<- variogram(object = VALUE~1,data = df.wg,locations =~X+Y, width= 200) plot(v)
通過觀察，已知點的自相關關系semivariance隨着距離distance的增加，呈現出指數形式的減弱（還記得嗎？離差平方和？值越小，代表關系越強），所以通過觀察半變異方差的分布，確定使用指數模型來擬合：
v.fit = fit.variogram(v, vgm(model = "Exp", psill= 45, range = 20000, kappa = 10),fit.sills = 50) plot(v, model=v.fit)
v.fit
如何確定pstill, range等初始值呢，同樣通過觀察大致給出即可，fit.variogram()函數會自動計算出正確的結果，這個初始值存在的作用只是輔助計算，不用十分正確，其結果如下圖：
model psill range 1 Exp 61.39472 5326.663
可以看出，R計算出了新的psill和range值。
接下來，我們就可以利用上述分析的信息，進行克里金插值。這里使用的是krige函數：
kingnewdata<- raster_to_dataframe2(max.allgf, "Interpolate location") aa<- krige(formula = VALUE~1,locations =~X+Y , model = vgm(61.39472, "Exp", 5326.663 , 10) , data = data.observed, newdata= kingnewdata, nmax=12, nmin=10)
參數formula指出已知點和其他信息的關系（假如本例已知點代表土壤含水量，我們也知道這些點對應的降雨信息，那么這里的公式就是土壤含水量和降雨的簡單線性關系），這里我們沒有其他信息，就是普通克里金插值。參數location是已知點的坐標，這里X代表經度，Y值代表緯度。參數model就是我們上邊分析的指數模型，使用擬合后的參數即可。data表示要使用的已知點的數據框。newdata表示要插值的點的位置，注意要包含和data參數所使用的dataframe一樣的變量名稱（本例中位置是大寫X， Y），nmax和nmin是最多和最少搜索的點數，其他參數大家可以參考幫助文檔，插值的結果如下：
好了，我們已經完成了克里金插值所有任務並且得到了我們所需要插值圖，但是怎么才能知道我們的插值結果可信與否呢？當然，如果我們有預測點的實際數據，我們可以評價插值結果的精度，但是很多情況下，我們沒有這些數據，cross validaion應運而生。對於克里金插值，gstat提供的cross validation函數是krige.cv()，對於本文，只需要如下代碼即可完成cross validation：
kriging<- krige.cv(formula = VALUE~1, locations = ~X+Y, model = vgm(61.39472, "Exp", 5326.663 , 10), data = data.observed ,nfold= nrow(data.observed))
head(kriging)
注意，krige.cv()函數本質上是進行很多次克里金插值，然后我們就可以分析被拿出的已知點的值和預測值，估計克里金插值的可信性。這里我們每次拿出一個點，所以nfold的值設置為和data.observed的行數一樣，可以看到結果：
var1.pred var1.var observed residual zscore fold X Y 1 19.020820 31.29921 12 -7.020820 -1.2549348 1 -318466.5 3794841 2 11.220626 27.62193 10 -1.220626 -0.2322499 2 -304466.5 3794841 3 10.758057 24.82611 9 -1.758057 -0.3528407 3 -303466.5 3794841 4 9.200462 24.85426 11 1.799538 0.3609613 4 -302466.5 3794841 5 10.840395 27.07824 7 -3.840395 -0.7380158 5 -301466.5 3794841 6 12.446044 29.40208 14 1.553956 0.2865826 6 -297466.5 3794841
結果包含了預測值，預測值的方差，已知值，殘差，z統計量值等，我們可以繪制出z統計量圖：
ggplot(data = kriging, aes(x=X,y=Y))+
geom_raster( aes(fill= zscore)) +
scale_fill_gradient2( name="zscore",low = "green",
mid = "grey", high = "red", midpoint = 0,
space = "rgb", na.value = "grey50")
結果如下：
或者其直方圖來評價插值方案：
par(mar=c(2,2,2,2)) hist(kriging$zscore)
歡迎大家留言討論，轉載請注明出處 http://www.cnblogs.com/ABMRG/p/5186727.html 。
參考文獻：
[1] http://desktop.arcgis.com/en/arcmap/10.3/tools/3d-analyst-toolbox/how-kriging-works.htm
關於R基本柵格數據處理我推薦一篇博文：http://rstudio-pubs-static.s3.amazonaws.com/1057_1f7e9ac569644689b7e4de78c1fece90.html