gstat | 空间插值（三）——克里金插值之泛克里金和简单克里金

本篇接着上篇继续介绍克里金插值。首先加载相关工具包和上篇使用的示例数据：

library(gstat)
library(sf)
library(tidyverse)
library(readxl)
load("G:/RStudies/空间插值/wh.rdata")
load("G:/RStudies/空间插值/stations.rdata")
data <- read_xlsx("G:/RStudies/空间插值/站点数据.xlsx", 1)
data <- left_join(stations, data, by = c("点位名称" = "监测点位"))
grid <- st_make_grid(wh, n = c(100, 100)) %>%
  st_transform(st_crs(data)) %>% st_sf()

3 泛克里金

在普通克里金的假设下，属性值在各个空间位置的数学期望都是同一个常数，而与空间位置和其他属性无关。泛克里金（Universal Kriging）则允许属性的数学期望与空间位置或其他属性相关，如在实际中对空气污染浓度进行插值时不仅要考虑已知监测点的浓度，还可以考虑海拔、人口密度等因素的影响。

除了其他要素外，属性值还可能具有一定的空间趋势，即空间位置本身也可以作为影响变量。在介绍泛克里金之前，先对趋势面分析做一个简单介绍。

3.1 趋势面分析

计算已知点和待插值点的投影坐标：

coord_data <- st_coordinates(st_geometry(data))
coord_grid <- st_coordinates(st_geometry(grid)) %>%
  data.frame %>% unique() %>%
  group_by(L1, L2) %>%
  summarise(X = mean(X), Y = mean(Y))

data$X = coord_data[,1]
data$Y = coord_data[,2]
grid$X = coord_grid$X
grid$Y = coord_grid$Y

趋势面分析可以使用回归模型来度量空间坐标与属性值之间的关系，这里使用最简单的线性回归（由于本篇使用的已知点过少，这里仅做演示）：

lm.wh <- lm(可吸入颗粒物 ~ X + Y, data)
summary(lm.wh)

## 
## Call:
## lm(formula = 可吸入颗粒物 ~ X + Y, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5020 -0.8074  0.6665  0.8134  2.7700 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  9.805e+02  4.928e+02   1.990   0.0938 .
## X            1.961e-04  9.466e-05   2.072   0.0837 .
## Y           -1.845e-04  1.305e-04  -1.414   0.2071  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.374 on 6 degrees of freedom
## Multiple R-squared:  0.5028, Adjusted R-squared:  0.3371 
## F-statistic: 3.034 on 2 and 6 DF,  p-value: 0.1229

使用回归结果对待插值点的属性值进行预测：

grid$pred <- predict(lm.wh, grid)

对预测结果进行可视化：

lm.wh <- st_interp(grid, st_union(wh))
wh.border <- st_boundary(wh)

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

3.2 泛克里金

趋势面分析没有使用到半变异函数，因此不属于克里金插值。在泛克里金中，可以将趋势面分析融入其中。

泛克里金与普通克里金的主要差别就在于其表达式。在上篇介绍普通克里金的推文中，使用的表达式是可吸入颗粒物 ~ 1，即只有截距而没有解释变量，而在泛克里金中，可以在表达式右侧加入相应的解释变量，本篇使用的表达式是可吸入颗粒物 ~ X + Y（这里仅做技术演示，在实际中可加入地形、人口等因子）。

泛克里金的函数和插值步骤和普通克里金类似：

vg.wh <- variogram(可吸入颗粒物 ~ X + Y, data)
vgm.wh <- vgm(psill = 10, model = "Bes",
              range = 10000, nugget = 0)
fit.wh <- fit.variogram(vg.wh, vgm.wh)
plot(vg.wh, fit.wh)

• variogram函数中的表达式为可吸入颗粒物 ~ X + Y。

universal.wh <- krige(可吸入颗粒物 ~ X + Y,
                      data, grid, fit.wh)
save(universal.wh, file = "31-1.universal.wh.rdata")

## [using universal kriging]

• krige函数中的表达式也为可吸入颗粒物 ~ X + Y。

插值结果可视化：

load("31-1.universal.wh.rdata")
universal.wh <- st_interp(universal.wh, st_union(wh))

ggplot() +
  geom_sf(data = universal.wh, aes(fill = var1.pred), col = NA) +
  geom_sf(data = wh.border) +
  theme_bw() +
  theme(text = element_text(family = "mono")) +
  scale_fill_continuous(low = "white", high = "red", name = "PM10")

4 简单克里金

简单克里金（Simple Kriging）与泛克里金的区别在于，其表达式中的截距和解释变量中的系数是已知的，在krige函数中使用beta参数表示。

这里通过对泛克里金的插值结果进行线性回归作为已知的系数：

lm(var1.pred ~ x + y, universal.wh)

## 
## Call:
## lm(formula = var1.pred ~ x + y, data = universal.wh)
## 
## Coefficients:
## (Intercept)            x            y  
##   1.290e+03    3.482e-04   -2.068e-04

simple.wh <- krige(可吸入颗粒物 ~ X + Y,
                   data, grid, fit.wh, 
                   beta = c(1290, 0.0003482, -0.0002068))
save(simple.wh, file = "31-2.simple.wh.rdata")

## [using simple kriging]

插值结果可视化：

load("31-2.simple.wh.rdata")
simple.wh <- st_interp(simple.wh, st_union(wh))

ggplot() +
  geom_sf(data = simple.wh, aes(fill = var1.pred), col = NA) +
  geom_sf(data = wh.border) +
  theme_bw() +
  theme(text = element_text(family = "mono")) +
  scale_fill_continuous(low = "white", high = "red", name = "PM10")

对于普通克里金而言，当其截距已知时也变成了简单克里金。

这里使用监测点的平均浓度作为已知的截距：

mean(data$可吸入颗粒物)

## [1] 26

vg2.wh <- variogram(可吸入颗粒物 ~ 1, data)
vgm2.wh <- vgm(psill = 10, model = "Bes",
               range = 10000, nugget = 0)
fit2.wh <- fit.variogram(vg2.wh, vgm2.wh)

simple2.wh <- krige(可吸入颗粒物 ~ 1,
                    data, grid, fit2.wh, 
                    beta = 26)
save(simple2.wh, file = "31-3.simple2.wh.rdata")

## [using simple kriging]

插值结果可视化：

load("31-3.simple2.wh.rdata")
simple2.wh <- st_interp(simple2.wh, st_union(wh))

ggplot() +
  geom_sf(data = simple2.wh, aes(fill = var1.pred), col = NA) +
  geom_sf(data = wh.border) +
  theme_bw() +
  theme(text = element_text(family = "mono")) +
  scale_fill_continuous(low = "white", high = "red", name = "PM10")

几种插值方法的结果比较：

想获取示例数据的读者请在公众号后台发送关键词“示例数据”（上篇已获取的读者无需重复获取）。

---