克里格空间插值Kriging及其在R的实现

克里格是一种空间插值的方法. 克里格插值, 是法国地统计学家 Georges Matheron (1930-2000)基于南非工程师Danie G. Krige (1919-2013)的硕士论文发展起来的方法. D. Krige硕士论文研究的问题是如何根据金矿探测点的含金量推断含金量在空间上的变化.

所谓空间插值, 就是通过一些列具有空间关系的点数据, 来推广到面上数据的方法.

除了克里格方法之外, 从点数据推断面的数据变化, 还有线性插值, 距离反平方插值等多种方法. 而克里格其实也是一类方法的通称, 它包括普通克里格, 通用克里格, 通用模块克里格等等. 但是克里格方法可以求得模型的无偏估计,并且方便得获得预测的误差, 这对于其他方法来讲常常难易实现.

这里只简单介绍用gstat程序包实现普通克里格, 以及通过一套模拟数据获得高程图.除了gstat程序包外, geoR程序包以及fields程序包也提供了克里格插值方法.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#### 首先生成一套模拟数据simdat, 通过Volcano数据随机提取1000个点, 用raster的extract函数提取

library(raster)
rrr <- raster(t(volcano), xmn=0, xmx= 870 , ymn=0, ymx=610)
x11()
plot(rrr, col = terrain.colors(100))

set.seed(12345)
loc <- data.frame(x = round(runif(1000, min = 0, max = 870), 2), y = round(runif(1000,min = 0, max = 610), 2))
value <- extract(rrr, loc)
simdat <- data.frame(loc, alt = value) ### The simulated data


#### Krig准备工作开始, 导入程序包
library(sp)
library(gstat)

#### simdat需要转换为sp格式, 需要指定坐标
coordinates(simdat) <- ~x+y
### bubble(simdat, "alt")
spplot(simdat, "alt") ### 查看每个点的值 

#### 创建克里格所用的网格, 20m方格的中心位置
xgrid <- seq(ceiling(min(simdat$x)), floor(max(simdat$x)), by = 20) 
ygrid <- seq(ceiling(min(simdat$y)), floor(max(simdat$y)), by = 20) 
basexy <- expand.grid(xgrid, ygrid)
#### plot(y ~ x, basexy)
colnames(basexy) <- c("x", "y")
coordinates(basexy) <- ~x+y
gridded(basexy) <- TRUE

### 因Kriging需要基于 Variogram以及相应的模型.Variogram横轴是距离, 纵轴是方差. 做模型拟合时, 应先查看 Variogram. Variogram根据点与点之间的距离, 分成若干段, 并计算每一段的方差, 再通过拟合Variogram模型, 获得Krig所需要的参数. 拟合Variogram模型时, 需要提供模型的初始值,再进行迭代优化, 否则很难获得模型的准确参数. 

######### 若所进行的是普通克里格插值, 则需要插值的变量只与自身取样点的距离相关,用公式表示为 alt~1, 这里的~1是variogram的规定 (参见 ?variogram). 
vgm1 <- variogram(alt~1, simdat)
plot(vgm1, plot.numbers = TRUE) 

### variogram的模型常包括
### Circular
### Spherical
### Exponential
### Gaussian
### Linear
### gstat程序包提供的模型可以通过 vgm() 查看

m <- fit.variogram(vgm1,vgm(.01,"Sph",100,300))
plot(vgm1, model=m)

在拟合gm模型时, 需要提供以下几个参数的初始值

  • psill: partial sill
  • model: 模型的类型
  • range: 范围
  • nugget: 块金值

图1 Range, Sill, Partial Sill, and Nugget of the Variogram

各参数的几何意义如图一所示. 初始值要通过在variogram的图中读取, 尽量接近真实值, 这样后续的迭代才会成功.

例如以下例子中:

  • .01 为 psill
  • “Sph” 为球状模型
  • 300 为范围
  • nugget为 块金值. 
1
2
3
4
5
6
7
8
### 设定vgm模型的初始值
m <- vgm(.01,"Sph",300, 0.2)

### 进行克里格
krige_res <- krige(alt~1, simdat, basexy, model = m)

### 查看克里格插值的结果
spplot(krige_res, zcol = "var1.pred", main = "Predictions of altitude based on the randomly sampled data", col.regions = terrain.colors(100))

参考:

http://planet.botany.uwc.ac.za/nisl/GIS/spatial/images/pic024.jpg
http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=How%20Krige%20and%20Variogram%20work

http://geostat-course.org/**system**/files/lewis_tutortPM.pdf
苏姝,林爱文,刘庆华 (2004) 普通 Kriging 法在空间内插中的运用 江南大学学报 (自然科学版) 3(1):18-21