gstat | 空间插值（四）——克里金插值之协同克里金和交叉验证

本篇是介绍克里金插值的第三篇推文，也是最后一篇。

因为前面两篇使用的数据中已知点的样本太少，本篇使用gstat工具包说明文档中的数据集。该数据集来自sp工具包。

实际上，gstat工具包的方法对sp和sf两种格式的空间矢量对象同样适用，且使用方法也一致。gstat的说明文档中是以sp为例的，而本系列的推文是以sf对象为例。

加载数据集：

library(sp)
data("meuse")
data("meuse.area")

class(meuse)

## [1] "data.frame"

names(meuse)

##  [1] "x"       "y"       "cadmium" "copper"  "lead"    "zinc"    "elev"   
##  [8] "dist"    "om"      "ffreq"   "soil"    "lime"    "landuse" "dist.m"

meuse是数据框结构，记录的是墨兹河附近区域的信息：

• x和y：点的空间位置；

• cadmium-zinc：镉、铜、铅、锌四种金属元素的含量；

• elev：高度；

• dist：距离墨兹河的距离。

前面已有推文介绍过通过数据框创建sf对象的方法了：

• sf | 创建空间矢量对象及其投影设置

library(sf)
data <- st_as_sf(meuse, coords = c(x = "x",
                                   y = "y"))

class(meuse.area)

## [1] "matrix"

colnames(meuse.area)

## [1] "x" "y"

meuse.area是矩阵结构，记录的是区域边界的点坐标。这里将其转换为sf线对象，作为插值边界：

border <- st_multilinestring(list(meuse.area))

创建插值网格、可视化锌含量的空间分布：

grid <- st_make_grid(border, n = c(100, 100))

plot(data["zinc"], reset = F)
plot(border, add = T)

加载相关工具包：

library(gstat)
library(tidyverse)

5 协同克里金

协同克里金（Cokriging）假设多种属性值相互之间存在相关性，因此可以同时对多个属性变量进行插值。协同克里金使用的函数是gstat工具包的同名函数：

gstat(g, id, formula, locations, data, model = NULL, beta,
      nmax = Inf, nmin = 0, omax = 0, maxdist = Inf, force = FALSE,
      dummy = FALSE, set, fill.all = FALSE,
      fill.cross = TRUE, variance = "identity", weights = NULL, merge, 
      degree = 0, vdist = FALSE, lambda = 1.0)

考虑到克里金插值要求变量满足空间平稳性，这里以四种金属元素含量的对数作为插值对象。计算已知点的协同半变异函数值的方法如下：

g <- gstat(NULL, "logCd", log(cadmium) ~ 1, data)
g <- gstat(g, "logCu", log(copper) ~ 1, data)
g <- gstat(g, "logPb", log(lead) ~ 1, data)
g <- gstat(g, "logZn", log(zinc) ~ 1, data)
vg <- variogram(g)

plot(vg$dist, vg$gamma)

拟合协同半变异函数的方法如下：

vgm <- vgm(psill = 1, model = "Sph",
           range = 800, nugget = 1)
fit <- fit.lmc(vg, g, vgm)
plot(vg, fit)

插值过程使用predict函数：

cokrige <- predict(fit, grid)
save(cokrige, file = "32-1.cokrige.rdata")

## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]

插值结果：

load("32-1.cokrige.rdata")
cokrige <- st_interp(cokrige, st_polygonize(border))
names(cokrige)

##  [1] "x"               "y"               "logCd.pred"      "logCd.var"      
##  [5] "logCu.pred"      "logCu.var"       "logPb.pred"      "logPb.var"      
##  [9] "logZn.pred"      "logZn.var"       "cov.logCd.logCu" "cov.logCd.logPb"
## [13] "cov.logCu.logPb" "cov.logCd.logZn" "cov.logCu.logZn" "cov.logPb.logZn"
## [17] "geometry"

• logCd.pred、logCu.pred等即为插值结果，logCd.var、logCu.var等为相应的插值结果的方差。

以锌为例对插值结果进行可视化：

ggplot(data = cokrige) +
  geom_sf(aes(fill = exp(logZn.pred)), col = NA) +
  geom_sf(data = border) +
  theme_bw() +
  theme(text = element_text(family = "mono")) +
  scale_fill_continuous(low = "white", high = "red", name = "Zn")

6 交叉验证

可以使用交叉验证（Cross Validation，CV）对插值结果进行评价。常使用的k折交叉验证（K-fold cross-validation）将样本平均分为k份，其中k-1份作为训练样本，剩余1份作为验证样本，重复k次。

单变量克里金插值（即非协同克里金）的交叉验证函数为krige.cv：

krige.cv(formula, locations, model = NULL, ..., beta = NULL, 
         nmax = Inf, nmin = 0, maxdist = Inf,
         nfold = nrow(locations), 
         verbose = interactive(), debug.level = 0)

• nfold：k参数；

• verbose：控制是否显示进度条的参数。

v.fit <- vgm(0.59, "Sph", 874, 0.04)
cv <- krige.cv(log(zinc) ~ 1, data, v.fit,
               nfold = 5, verbose = T)

cor(cv$var1.pred, cv$observed)

## [1] 0.8510259

协同克里金的交叉验证函数为gstat.cv：

gstat.cv(object, nfold, remove.all = FALSE,
         verbose = interactive(), 
         all.residuals = FALSE, ...)

gstat.cv函数默认只对第一个属性变量进行交叉验证：

• object：gstat函数的输出结果；

• all.residuals：若为TRUE，则生成所有变量的预测残差。

cv2 <- gstat.cv(g, verbose = T)
cor(cv2@data$logCd.pred, cv2@data$observed)

## [1] 0.5563733

---