ABSTRACT

Spatial kriging is one of the most important interpolation methods in geostatistics, which models the spatial autocorrelation of the interested variable in the region and provides the unbiased optimized estimation. This article studies the soil moisture in the upstream of the Heihe River, and analyzes the spatial variability of the observations using the R software packages. The spatial kriging model was constructed to calculate the estimation of the soil moisture in the whole region and provide the necessary index for further study of hydrological process in the region.

Keywords:Geostatistics, Regional Variable, Spatial Kriging, Soil Moisture

基于空间克里金插值的黑河流域土壤水分研究

张晞曈

长安大学，陕西 西安

收稿日期：2018年11月20日；录用日期：2018年12月6日；发布日期：2018年12月13日

摘 要

空间克里金是地统计学中重要的一类插值方法，它通过对研究变量在研究区域内的空间自相关性进行建模，对待估计的点给出无偏最优估计。本文以黑河上游的典型灌溉农田土壤水分为研究对象，利用R软件中的相关扩展包，分析观测数据的空间变异性，构建研究区域内土壤水分含量的空间克里金插值模型，得到该区域内土壤水分含量的空间估计，为进一步研究该流域的生态水文过程提供了必要的指标。

关键词 :地统计学，区域化变量，空间克里金，土壤水分

Copyright © 2018 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

地统计学(又称地质统计学)是在传统的统计学基础上，以区域化变量为基础，采用变异函数为主要工具，研究在空间分布上具有自相关性的自然现象的科学 [1] 。克里金(Kriging)插值法，又称空间局部估计法，是地统计学的主要内容之一。它不同于传统基于距离的空间插值法，比如反距离加权法，只考虑了待测点周围的数据，直接通过插值公式得到结果 [2] 。克里金法是从研究变量本身出发，根据已有样本点和待测点之间的空间位置关系来对自相关性进行建模，从而对待测点进行线性无偏的估计，同时还能给出估计的方差，比传统方法能更好地描述所研究的地理现象 [3] 。

土壤是联系大气水、地表水、地下水和植物水并进行水分交换的核心地带，揭示土壤水文性质的空间分布特征是研究流域生态–水文过程的重要一环 [4] ，土壤水分含量是水文过程分析中重要的一个指标。本文以流经青海省的黑河上游八宝河流域的土壤水分为研究对象，根据已有的无线传感器网络点收集的土壤水分数据，利用空间普通克里金法对研究区域内的土壤水分进行插值，得到整个地区土壤水分的分布估计，并给出估计方差，为进一步的水文环境分析提供了必要的指标。

2. 空间克里金方法简介

2.1. 区域化变量

区域化变量是指能用空间分布来表征一个自然现象的变量，也就是将空间位置作为随机函数的自变量，它是地统计学中重要的基础理论 [5] 。设研究区域的位置集合为D ( $D\subseteq R^2$ 或者 $D\subseteq R^3$ )，空间位置表示为 $x\in D$ ，所研究的变量在该处的属性值为 $Z\left(x\right)$ 。也就是说，将所研究的地质现象看作研究区域D内以空间位置 为参数的随机函数 $Z\left(x\right)$ 。

区域化变量在描述研究对象的空间结构性时，是从空间点对的角度出发的。设研究区域内任意两点为 $x_1\in D,x_2\in D$ ，两点之间的距离向量为 $h\in D$ ，则这两点之间属性的空间相关性通常用 $Z\left(x_1\right),Z\left(x_2\right)$ 之间的一阶矩和二阶矩来刻画。

区域化变量 $Z\left(x\right)$ 的自协方差函数定义为空间两点x和 $x+h$ 处的随机变量 $Z\left(x\right)$ 和 $Z\left(x+h\right)$ 的二阶混合中心矩，即

$C\left(x,h\right)=\text{Cov}\left[Z\left(x\right),Z\left(x+h\right)\right]=E\left[Z\left(x\right)Z\left(x+h\right)\right]-E\left[Z\left(x\right)\right]\cdot E\left[Z\left(x+h\right)\right]$ (1)

通常出于研究需要，对于区域化变量的性质需要做出一些假设，最常见的就是如下的平稳假设 [5] 。当区域化变量满足下列两个条件时，称该区域化变量满足二阶平稳假设：

1).在研究区域内，区域化变量 $Z\left(x\right)$ 的数学期望存在且为常数：

$E\left[Z\left(x\right)\right]=E\left[Z\left(x+h\right)\right]=m$ (2)

2). 在研究区域内，区域化变量 的协方差存在且只依赖于两点之间空间距离h：

$\text{Cov}\left[Z\left(x\right),Z\left(x+h\right)\right]=E\left[Z\left(x\right)Z\left(x+h\right)\right]-m^2$ (3)

在二阶平稳假设下，协方差函数 $C\left(x,h\right)$ 可以简记为 $C\left(h\right)$ ，也就是只和点对之间的距离向量有关。

2.2. 变异函数

变异函数是地统计学特有的工具，是许多地统计学计算的基础，它定义为区域化变量 $Z\left(x\right)$ 在空间点x和 $x+h$ 处属性值之差 $Z\left(x\right)-Z\left(x+h\right)$ 的方差的一半，即

$\begin{array}{c}\gamma \left(x,h\right)=\frac{1}{2}Var\left[Z\left(x\right)-Z\left(x+h\right)\right]\\ =\frac{1}{2}E\left\{{\left[Z\left(x\right)-Z\left(x+h\right)\right]}^2\right\}-\frac{1}{2}{\left\{E\left[Z\left(x\right)\right]-E\left[Z\left(x\right)\right]\right\}}^2\end{array}$ (4)

严格来说，变异函数应当叫做半变异函数，因为前面有系数，但是这只是为了公式推导方便，所以通常不做区分。

在二阶平稳的假设下，变异函数和协方差函数之间通常有如下关系 [6] ：

$\gamma \left(x,h\right)=\gamma \left(h\right)=C\left(0\right)-C\left(h\right)$ (5)

变异函数是地统计学研究的出发点，对于空间变量之间的自相关性的刻画都是从变异函数出发而不是从协方差函数出发。一般来说，两点之间距离越远，也就是h越大，变差函数的值越大，此时两点之间的相关性越小。

在二阶平稳假设下变异函数是关于距离h的函数，该函数有指数模型、球状模型和高斯模型等 [3] 。实际应用中，变异函数的具体表达式需要通过拟合得到。首先根据已有样本点数据计算每两个样本点之间的距离h和相应的实验变异函数值 ${\left[Z\left(x\right)-Z\left(x+h\right)\right]}^2$ ，n个样本点两两之间总共可以得到 $n\left(n-1\right)/2$ 对值，以h为横坐标、 ${\left[Z\left(x\right)-Z\left(x+h\right)\right]}^2$ 为纵坐标对这些值绘制散点图可以得到变异函数云图 [7] ，变异函数云图通常用于分析原始数据之间是否存在自相关性。然后对上述 $n\left(n-1\right)/2$ 对值进行聚合，根据如下(6)式计算实验变异函数值 [7] ，

${\gamma }^{*}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{{\displaystyle \sum }}}{\left[Z\left(x_i\right)-Z\left(x_i+h\right)\right]}^2$ (6)

这里 $N\left(h\right)$ 表示距离之差为 的点对数量。最后选定变异函数的模型，使用上述得到的实验变异函数值对其中的参数进行拟合得到理论变异函数的表达式 [6] 。

2.3. 空间克里金插值

不同的克里金插值法在估计的条件和细节上略有不同，本文以普通克里金(Ordinary Kriging)为例，阐述克里金插值法的基本原理 [3] 。克里金插值的基本思想是利用空间上所有已知点的属性值，进行加权求和来估计未知点处的值，本质上也是一个线性估计量。设已知观测点集合为 $\left\{x_i|i=1,2,\cdots ,n\right\}$ ，空间点 $x_i$ 处的属性值为 $Z\left(x_i\right)$ ，简记为 $Z_i$ 。待估计的点为 $x_0$ ，待估计点的实际值为 $Z\left(x_0\right)$ ，简记为 $Z_0$ ，估计值为 $Z{\left(x_0\right)}^{\text{*}}$ ，简记为 $Z_0^{\text{*}}$ 。

待测点 $x_0$ 的估计值 $Z_0^{\text{*}}$ 由如下的线性估计式给出

$Z_0^{*}={\displaystyle {\sum }_{i=1}^n{\lambda }_i{Z}_i}$ (7)

关键在于如何确定这些权重 ${\lambda }_i\left(i=1,2,\cdots ,n\right)$ 的值。

克里金方法是一种在满足无偏性的条件下寻求估计方差最小的权重系数，所以克里金插值问题可以写成一个有约束条件的优化问题 [3] ：

$\{\begin{array}{ll}\min \hfill & \text{Var}\left[Z_0^{*}-Z_0\right]\hfill \\ \text{s}\text{.t}\hfill & E\left[Z_0^{*}\right]=Z_0\hfill \end{array}$ (8)

上述优化问题可以展开成如下的优化问题模型：

$\{\begin{array}{ll}\min \hfill & C\left(Z_0,Z_0\right)+\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\underset{j=1}{\overset{n}{{\displaystyle \sum }}}{\lambda }_i{\lambda }_{j}C\left(Z_i,Z_j\right)-2\underset{i=1}{\overset{n}{{\displaystyle \sum }}}{\lambda }_{i}C\left(Z_0,Z_i\right)\hfill \\ \text{s}\text{.t}\hfill & \underset{i=1}{\overset{n}{{\displaystyle \sum }}}{\lambda }_i=1\hfill \end{array}$ (9)

这里 $C\left(Z_i,Z_j\right)$ 表示 $C\left(\left|x_i-x_j\right|\right)$ ，也就是根据点对 $x_i$ 和 $x_j$ 之间的距离计算得到的协方差值，这个值通常是利用(5)式由已经拟合的理论变异函数计算得到。

(9)式的优化问题可以采用拉格朗日法转化为无约束优化问题求解出一系列的权重值，从而得到待估计点 $x_0$ 的最小方差无偏估计。

3. 黑河土壤水分的空间克里金插值

本文采用R软件中的gstat包来对黑河流域的土壤水分含量进行空间克里金插值。空间克里金插值一般分为三个步骤：

1) 根据已有样本点数据绘制变异函数云图和实验变异函数图，进行初步分析；

2) 选择变异函数模型，拟合实验变异函数；

3) 对待估计点进行空间克里金插值，给出估计方差。

3.1. 数据来源

本文使用了来自于黑河数据管理计划寒区旱区科学数据中心(http://westdc.westgis.ac.cn/)的黑河生态水文遥感试验数据集 [8] ，该数据集由黑河流域上游生态水文无线传感器网络WATERNET于2013年收集，这里使用了当年8月24号的观测数据。该数据集包括黑河上游八宝河流域内40个WATERNET节点的观测数据，这40个节点的数采仪编号为10~49，但是其中第46和48号数采仪在8月24号这天的数据由于外界因素有缺失，所以实际使用的是38个节点提供的观测数据，这些节点的分布如图1所示 [8] 。

图1. 黑河上游区域WSN节点分布

每个节点都含有4 cm，10 cm，20 cm三个深度的土壤水分数据，这里以4 cm深度的土壤水分含量为对象。数采仪每5分钟采集一次数据，这样的观测频率过于密集，所以需要对原始数据进行聚合，这里以天为单位，每个节点取8月24号观测数据的平均值。处理后的各节点分布以及土壤水分含量分布见图2 (单位：cm³/cm³)。

图2. 土壤水分含量观测值分布(单位：cm³/cm³)

3.2. 绘制变异函数云图和实验变异函数图

首先根据已有38个节点的土壤水分含量数据绘制变异函数云图和实验变异函数图。由于不同方向上土壤水分的变化趋势可能不相同，表现为各向异性 [3] ，这里选择了正北(0˚)、东北(45˚)、正东(90˚)和东南(135˚)四个方向分别绘制图像，以判断变化趋势较为明显的方向。绘制的变异函数云图见图3，实验变异函数图见图4。

图3. 四个方向的变异函数云图

由变异函数云图和实验变异函数图可以看出，正北和东北两个方向的点对数量较少，变异函数云图上表现为图3中点的数量较少，不利于变异函数的计算，同时由图4中也可以看出，这两个方向上的实验变异函数散点图过于稀疏，无法反映自相关性。正东和东南两个方向的点对数量较多，但是在图4的实验变异函数图像中，正东方向的变化程度没有东南方向的显著，东南方向的实验变异函数散点图分布趋势较为陡峭，而正东方向的实验变异函数散点图分布较为平缓。所以最终选择东南方向的实验变异函数来进行拟合。

图4. 四个方向的实验变异函数图

同时，应当注意到，在滞后距离(lag distance)超过60 km的时候，图3中东南方向的变异函数云图中的点明显变少，这是导致图4中东南方向的实验变异函数图像忽然下降的原因，点的数目变少后，计算得到的方差也就降低了，也就是说，过少的点数计算得到的方差不具有代表性，所以这里只选择滞后距离在65 km以内的实验变异函数点进行拟合。

3.3. 拟合变异函数模型

得到上述的实验变异函数图像后，需要选择变异函数模型来进行拟合。常用的变异函数模型有球状模型和指数模型 [3] 。球状模型的表达式如(10)式所示：

$\gamma \left(h\right)=\{\begin{array}{ll}0,\hfill & h=0\hfill \\ C_0+C\left(\frac{3}{2}\cdot \frac{h}{a}-\frac{1}{2}\cdot \frac{h^3}{a^3}\right),\hfill & 0<h\le a\hfill \\ C_0+C,\hfill & h>a\hfill \end{array}$ (10)

指数模型的表达式如(11)式所示：

$\gamma \left(h\right)=\{\begin{array}{ll}0,\hfill & h=0\hfill \\ C_0+C\left(1-{\text{e}}^{-\frac{h}{a}}\right),\hfill & h>0\hfill \end{array}$ (11)

上述两式中， $C_0$ 为块金常数，C为偏基台值，a为变程。

拟合后得到的两个模型参数见表1：

表1. 拟合的模型参数

拟合后两个模型的图形见图5：

图5. 变异函数拟合模型，左为指数模型，右为球状模型

由表1可知，指数模型的拟合误差SSE相比球状模型的较小，所以这里最终采用的是指数模型的变异函数。

3.4. 空间克里金插值

得到了变异函数模型之后，选择需要插值的空间点，然后根据(9)式求解相应的优化问题，即可得到空间克里金插值。本文对于整个黑河上游八宝河区域都进行了插值，即先对八宝河流域网格化，然后在对每个网格中的点进行空间克里金插值，得到整个八宝河流域的土壤水分含量估计分布图，同时给出估计的方差。空间克里金插值的结果和方差分布如图6所示。

图6左边是空间克里金在整个区域的插值分布图，右侧是相应的方差分布图。左图可以看出八宝河流域土壤水分含量的大致分布，右图中标星号的点是已有的观测点，可以看出，在观测点附件的估计方差较小，远离观测点的地方估计方差较大，这与经验是相符合的。

4. 结论

本文利用地统计学中的空间克里金模型，根据已有观测点收集的数据，对黑河上游八宝河流域的土壤

图6. 空间克里金插值结果，左为插值图，右为方差图

水分进行了插值估计，并给出了估计方差，得到该流域土壤水分分布的一个估计，为进一步的水文环境分析提供了必要的指标。空间克里金模型相对于反距离加权法等传统的空间插值方法，不仅考虑了空间距离，还考虑了空间中点与点之间的相对位置以及自相关性，并且能给出估计精度，具有较高的实用性和可靠性。

文章引用

张晞曈. 基于空间克里金插值的黑河流域土壤水分研究
Study on Soil Moisture of Heihe River Based on Spatial Kriging[J]. 统计学与应用, 2018, 07(06): 622-629. https://doi.org/10.12677/SA.2018.76071

参考文献