说明：昨天的推文误把可吸入颗粒物当作PM2.5，实应该为PM10，这里修正后重发。

从本篇开始计划分三篇介绍克里金插值。与反距离权重插值不同，克里金插值是无偏估计，其中也涉及到模型估计。本篇先对普通克里金的原理进行简单介绍，然后重点介绍在R语言中进行克里金插值的步骤。

library(gstat)

1 示例数据

本篇主要对武汉市主城区范围的空气质量监测点数据进行插值。

加载空间矢量对象：

library(sf)
library(tidyverse)
library(readxl)
load("G:/RStudies/空间插值/wh.rdata")
load("G:/RStudies/空间插值/stations.rdata")ggplot() +geom_sf(data = wh) +geom_sf(data = stations, col = "red") +theme_bw() +theme(text = element_text(family = "mono"))

• 红点为空气质量监测点位置。

属性数据是2021年5月15日的监测数据：

library(readxl)data <- read_xlsx("G:/RStudies/空间插值/站点数据.xlsx", 1)
head(data)## # A tibble: 6 x 12
##   日期                监测点位 二氧化硫 二氧化氮 可吸入颗粒物 一氧化碳  臭氧
##   <dttm>              <chr>       <dbl>    <dbl>        <dbl>    <dbl> <dbl>
## 1 2021-05-15 00:00:00 东湖梨园        6       23           23       23    47
## 2 2021-05-15 00:00:00 汉阳月湖        6       33           26       19    42
## 3 2021-05-15 00:00:00 汉口花桥        3       35           26       17    45
## 4 2021-05-15 00:00:00 武昌紫阳        2       39           25       23    43
## 5 2021-05-15 00:00:00 青山钢花        5       35           31       24    30
## 6 2021-05-15 00:00:00 沌口新区        5       39           25       32    45
## # ... with 5 more variables: 细颗粒物 <dbl>, 空气质量指数 <dbl>,
## #   首要污染物 <chr>, AQI指数级别 <chr>, AQI指数类别 <chr>

• 数据中的变量是以中文命名的，但这并不影响后续分析。

进行属性数据连接、在武汉市主城区范围内生成100*100的网格、对监测点可吸入颗粒物即PM10的浓度进行可视化：

data <- left_join(stations, data, by = c("点位名称" = "监测点位"))
grid <- st_make_grid(wh, n = c(100, 100)) %>% st_transform(st_crs(data))ggplot() +geom_sf(data = wh) +geom_sf(data = data, aes(col = 可吸入颗粒物)) +theme_bw() +theme(text = element_text(family = "mono")) +scale_color_continuous(low = "white", high = "red", name = "PM10")

2 普通克里金

2.1 理论基础

与反距离权重插值一样，克里金插值同样是使用已知点的属性加权值来作为未知点属性的预测值：

不同点在于，克里金插值是使用广义最小二乘法对以上表达式进行无偏和最优估计来得到权重值的。

普通克里金（Ordinary Kriging）是克里金插值的一种，它的假设是属性值在空间上是平稳的，即

• 表示属性值，表示空间位置，为常数，表示方差恒定的随机扰动项。

基于上述假设，普通克里金的权重估计需要满足以下约束条件：

无偏估计

因为，则。

最优估计

广义最小二乘法要求属性预测值与真实值之差的方差最小：

由于由加权而得，那么就必然与各点属性值的协方差有关了。[https://xg1990.com/blog/archives/222]详细推导了普通克里金的求解过程，这里直接给出最终结果：

其中是变异函数，也叫半方差函数或半变异函数：

半变异函数可以认为是两个位置属性相关性的度量。普通克里金插值假设，半变异函数只与两个位置的空间距离有关：

因此，求解上述方程的关键就是寻找的拟合函数并用以代替。

2.2 插值函数和步骤

进行克里金插值共需要三个步骤：首先是计算已知点之间的半变异函数值，然后寻找合适的拟合函数，最后使用拟合函数对模型进行求解和预测。这里以可吸入颗粒物（PM10）的浓度为例。

计算已知点之间的半变异函数值

计算函数为variogram：

variogram(object, locations, X, cutoff, width = cutoff/15,alpha = 0, beta = 0, tol.hor = 90/length(alpha), tol.ver = 90/length(beta),cressie = FALSE, dX = numeric(0), boundaries = numeric(0),cloud = FALSE, trend.beta = NULL, debug.level = 1,cross = TRUE, grid, map = FALSE, g = NULL, ..., projected = TRUE, lambda = 1.0, verbose = FALSE, covariogram = FALSE, PR = FALSE, pseudo = -1)

主要参数有：

• object：模型对象；

• locations：已知点所在的空间矢量对象。

vg.wh <- variogram(可吸入颗粒物 ~ 1, data)
vg.wh##   np      dist     gamma dir.hor dir.ver   id
## 1  1  3386.683  0.000000       0       0 var1
## 2  2  5133.450  0.500000       0       0 var1
## 3  3  7166.853  2.166667       0       0 var1
## 4  2  7523.818 20.000000       0       0 var1
## 5  2  8480.413  0.250000       0       0 var1
## 6  1  9661.245 24.500000       0       0 var1
## 7  2 10230.623  6.250000       0       0 var1
## 8  2 11136.533 10.250000       0       0 var1
## 9  1 11936.677  4.500000       0       0 var1

• 输出结果中的gamma变量即计算的半变异函数值。

使用plot函数绘制距离dist与边缘函数gamma之间的散点图：

plot(vg.wh$dist, vg.wh$gamma)

拟合函数

典型的半变异函数如下图所示：

相关术语：

• range：变程；

• nugget：块金；

• partial sill：偏基台值；

• sill：基台，等于块金+偏基台值。

vgm函数用于确定拟合函数的形式：

vgm(psill = NA, model, range = NA, nugget, add.to, anis, kappa = 0.5, ..., covtable,Err = 0)

主要参数有：

• psill：偏基台值的初始值；

• model：拟合函数的模型形式；

• range：变程；

• nugget：块金的初始值。

查看model参数的可选项：

vgm()##    short                                      long
## 1    Nug                              Nug (nugget)
## 2    Exp                         Exp (exponential)
## 3    Sph                           Sph (spherical)
## 4    Gau                            Gau (gaussian)
## 5    Exc        Exclass (Exponential class/stable)
## 6    Mat                              Mat (Matern)
## 7    Ste Mat (Matern, M. Stein's parameterization)
## 8    Cir                            Cir (circular)
## 9    Lin                              Lin (linear)
## 10   Bes                              Bes (bessel)
## 11   Pen                      Pen (pentaspherical)
## 12   Per                            Per (periodic)
## 13   Wav                                Wav (wave)
## 14   Hol                                Hol (hole)
## 15   Log                         Log (logarithmic)
## 16   Pow                               Pow (power)
## 17   Spl                              Spl (spline)
## 18   Leg                            Leg (Legendre)
## 19   Err                   Err (Measurement error)
## 20   Int                           Int (Intercept)

• 常用的模型形式有Exp、Sph、Gau、Mat和Pow。

vgm函数的相关参数主要通过观察已知点之间的变异函数值的散点图来设置，如：

vgm.wh <- vgm(psill = 10, model = "Bes",range = 10000, nugget = 0)
vgm.wh##   model psill range
## 1   Nug     0     0
## 2   Bes    10 10000

fit.variogram用于拟合半变异函数：

fit.variogram(object, model, fit.sills = TRUE, fit.ranges = TRUE,fit.method = 7, debug.level = 1, warn.if.neg = FALSE, fit.kappa = FALSE)

主要参数有：

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

• model：vgm函数的输出结果。

fit.wh <- fit.variogram(vg.wh, vgm.wh)
plot(vg.wh, fit.wh)

属性值预测

克里金插值的预测函数为krige：

krige(formula, locations, newdata, model, ...,beta, nmax = Inf, nmin = 0, omax = 0,maxdist = Inf, block, nsim = 0,indicators = FALSE,na.action = na.pass, debug.level = 1)

• 对于普通克里金而言，大部分主要参数及其含义同前篇介绍的idw函数；

• model：fit.variogram函数的输出结果。

krige.wh <- krige(可吸入颗粒物 ~ 1, data, grid,model = fit.wh)
save(krige.wh, file = "30-1.krige.wh.rdata")## [using ordinary kriging]

load("30-1.krige.wh.rdata")
krige.wh## Simple feature collection with 10000 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -1524661 ymin: 3552782 xmax: -1474924 ymax: 3586638
## Projected CRS: Beijing 1954 / Gauss-Kruger CM 135E
## First 10 features:
##             x       y var1.pred var1.var                       geometry
## ID1  -1524412 3552951  27.72476 17.56884 POLYGON ((-1524661 3552782,...
## ID2  -1523915 3552951  27.82153 17.28171 POLYGON ((-1524164 3552782,...
## ID3  -1523418 3552951  27.91942 17.02860 POLYGON ((-1523666 3552782,...
## ID4  -1522920 3552951  28.01840 16.80727 POLYGON ((-1523169 3552782,...
## ID5  -1522423 3552951  28.11842 16.61532 POLYGON ((-1522672 3552782,...
## ID6  -1521926 3552951  28.21943 16.45019 POLYGON ((-1522174 3552782,...
## ID7  -1521428 3552951  28.32139 16.30922 POLYGON ((-1521677 3552782,...
## ID8  -1520931 3552951  28.42425 16.18964 POLYGON ((-1521180 3552782,...
## ID9  -1520433 3552951  28.52796 16.08860 POLYGON ((-1520682 3552782,...
## ID10 -1519936 3552951  28.63249 16.00325 POLYGON ((-1520185 3552782,...

• 使用克里金插值不仅能预测位置点的属性值var1.pred，还能估计预测值的方差var1.var。

插值结果可视化

krige.wh <- st_interp(krige.wh, st_union(wh))
wh.border <- st_boundary(wh)ggplot() +geom_sf(data = krige.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")

想获取示例数据的读者请在公众号后台发送关键词“示例数据”。