随着中国城镇化进程的不断推进，为了更有效地保护优质耕地，国土资源部提出了开展低丘缓坡等未利用土地开发试点的工作。残坡积土是低丘缓坡重要的组成部分，高质量的残坡积土层厚度信息在土地资源管理、工程建设、地质灾害预警预报等方面具都有重要意义。获取比较精确的土层厚度数据的方法之一是布设高密度的钻孔, 但由于受到经济水平、技术手段和地形条件的限制, 钻孔点的数量是一定的，这为确定土层厚度带来了较大的误差。随着地统计学理论在工程界不断应用，空间插值的方法也日益成为解决上述问题的经典方法。将利用统计学与地理信息系统相结合的方法, 基于已知钻孔点的观测数据进行空间插值来获得估计点的土层厚度。

空间插值是利用已知的部分空间样本信息, 对未知地理空间的特征属性值进行估计, 是地理信息系统的重要功能模块之一^[1], 并在矿业、水文、气候预报、农业等领域有着广泛的应用。目前已发展了较多的空间插值方法：如泰森多边形法^[2]、克里金(Kriging)插值法^[3]、反距离加权平均法^[4]、趋势面分析法^[5]、多项式回归法^[6]等。但这些方法只是局限于已知点的单一属性值，没有考虑到其他因素对待估计点属性值的影响(如高程对土层厚度的影响等)。

基于上述方法在估计空间属性值时对变量考虑的单一性，学者提出了一种能考虑多个相关变量互相影响的协同克里金法。Yates等^[7]将裸土表面温度和土壤沙粒含量作为协变量，利用协同克里金法(Co-Kriging)估计质量含水率；Ghadermazi等^[8]利用pH值作为辅助变量，估计了饮用水中硝酸盐的含量。胡丹桂等^[9]用考虑降雨量的协克里金法来研究空气的湿度的空间变化；许晶玉等^[10]用考虑土壤粗砂含量和全氮含量的协克里金法来研究山东省种植区地下水硝态氮污染空间变异及分布规律；杜文凤等^[11]将地震数据引入协克里金法中估计煤层厚度的分布规律；黄大年等^[12]首次将重力梯度引入协克里金法中来研究岩脉倾向的问题。

半变异函数是克里金法中反映区域化变量空间变化特征的有效数学模型, 由其确定拟合模型参数直接影响插值精度。严华雯等^[13]通过利用加权最小二乘法优化遗传算法中的适度函数，改进普通基于遗传算法优化的克里金插值方法；张强等^[14]采用基于线性递减权值的粒子群算法估计半变异函数参数的方法，提高了插值的精度；贾雨等^[15]将约束粒子群优化算法用于克里金的插值研究，结果表明能获得较好的插值精度。

本文采用粒子群优化算法与协同克里金法相结合的方法对土层厚度的空间插值进行研究，并对插值结果进行对比，以期提高土层厚度插值分布的精度。

1 方法原理

1.1 普通克里金插值法

克里金插值法又称空间自协方差最佳插值法，是以南非矿业工程师D.G.Krige的名字命名的一种最优内插法。该方法首先考虑的是空间属性在空间位置上的变异分布，确定对一个待插点值有影响的距离范围，然后用此范围内的采样点来估计待插点的属性值。该方法考虑了已知信息与待估计点相互间的空间结构特性。为达到线性、无偏和最小估计方差的估计，而对每一个样品赋予一定的权重，最后进行加权平均来估计待预测点属性值的一种方法。

在克里金法中，用来衡量各个样本点之间空间相关程度的是半变异函数

$ r\left( h \right) = \frac{1}{{2N}} \cdot \sum\limits_{i = 1}^N {{{\left[ {Z\left( {{x_i}} \right) - Z\left( {{x_i} + h} \right)} \right]}^2}} $

(1)

式中：h为两点之间距离；N是由h分开的成对样本点的数量；Z是点的属性值。

一种典型的半变异函数图像如图 1所示，半变异值随距离增大而增大，其中有两个非常重要的点：间距为0时的点以及函数趋于平稳时的拐点，前者表示的是两个非常接近的样本点之间的误差及空间变异，后者表示两样本点超过此间离后将不存在空间相关性。

图 1

图 1 半变异函数图 Fig. 1 Semivariance image

通过式(1)确定半变异函数云图，然后再拟合出相应的模型表达式，应用式(2)确定内插所需要的权重，并通过式(3)进行未知样本点属性值的估计。

$ \left\{ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^N {{\lambda _i} \cdot r\left( {{h_{ij}}} \right) + u} = r\left( {{h_{0j}}} \right)}\\ {\sum\limits_{i = 1}^N {{\lambda _i}} = 1} \end{array}} \right. $

(2)

$ {Z^*}\left( {{x_0}} \right) = \sum\limits_{i = 1}^N {{\lambda _i} \cdots {Z_i}\left( x \right)} $

(3)

式中：λ_i为待定的权重；u为拉格朗日乘子；r(h_ij)为x_i与x_j两点的半变异函数；N为样本点数量；x₀为待插值点，Z为其属性值。

1.2 协同克里金法

协同克里金插值法可以利用同一变量在不同时空或不同变量在同一时空上的协同区域化性质，用易于测定的变量来对那些难以测定的属性或变量进行估值，或者用样品多的变量对样品少的变量进行估值；如果两种以及多种属性具有显著的空间相关性，则可以优先考虑协同克里金插值法。

在协同克里金法中，γ₁₁、γ₂₂分别表示变量1与变量2的半变异模型函数，其计算方法与普通克里金插值法中的半方差计算方法相同；γ₁₂为协半变异函数，是用来衡量两个变量各个样本点之间空间相关程度的表达式

$ \begin{array}{*{20}{c}} {{\gamma _{12}} = \frac{1}{{2N}} \cdot \sum\limits_{i = 1}^N {\left[ {{Z_1}\left( {{x_i}} \right) - {Z_1}\left( {{x_i} + h} \right)} \right] \cdot } }\\ {\left[ {{Z_2}\left( {{x_i}} \right) - {Z_2}\left( {{x_i} + h} \right)} \right]} \end{array} $

(4)

式中：h为两点之间距离；N是由h分开的成对样本点的数量; Z₁是点关于变量1属性值；Z₂是点关于变量2属性值。

拟合后用式(5)即可确定空间插值所需要的权重，并通过式(6)进行未知样本点属性值的估计。

$ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}{\gamma _{11}}\left( {{x_{1i}} - {x_i}} \right)} + \sum\limits_{j = 1}^{{N_2}} {{\lambda _{2j}}{\gamma _{21}}\left( {{x_{2j}} - } \right.} }\\ {\left. {{x_i}} \right) + {u_2} = {\gamma _{21}}\left( {{x_0} - {x_i}} \right)} \end{array}\\ \begin{array}{*{20}{c}} {\sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}{\gamma _{21}}\left( {{x_{1i}} - {x_j}} \right)} + \sum\limits_{j = 1}^{{N_2}} {{\lambda _{2j}}{\gamma _{22}}\left( {{x_{2j}} - } \right.} }\\ {\left. {{x_j}} \right) + {u_2} = {\gamma _{22}}\left( {{x_0} - {x_j}} \right)} \end{array}\\ \sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}} = 0,\sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}} = 0 \end{array} \right. $

(5)

$ Z_2^*\left( {{x_0}} \right) = \sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}} {Z_1}\left( {{x_{1i}}} \right) + \sum\limits_{j = 1}^{{N_2}} {{\lambda _{2j}}{Z_2}\left( {{x_{2j}}} \right)} $

(6)

式中：λ_1i为变量1待定的权重；λ_2j为变量2待定的权重；N₁为变量1样本点数量；N₂为变量2样本点数量；γ₁₁为变量1的半方差函数；γ₂₂为变量2的半方差函数；γ₁₂为变量1与变量2的协半方差函数，γ₁₂=γ₂₁；Z₁(x_1i)为变量1的属性值；Z₂(x_2j)为变量2的属性值；x₀为待插值点，Z₂^*(x₀)为未知点的属性值

1.3 粒子群优化算法

粒子群优化算法(PSO)是一种进化演变技术，最早是由美国的心理学研究者Kennedy博士和从事计算智能研究的Eberhart博士受到人工生命的研究结果启发于1995年提出的一种基于群智能的优化算法^[16-17]。

粒子群算法的实质是一种信息共享，粒子根据自身的个体最优信息及群体的最优信息不断更新自己的速度和位置，最后收敛于全局最优解。每个粒子都有一个有待优化问题所决定的适应度值用来评价该粒子的优良程度，一个粒子还有一个速度用于决定它的运动距离即速度，粒子根据自身的运动经验和群体中其它粒子的运动经验来调整自己的运动。粒子从此时刻运动到下一时刻的速度由式(7)确定，位置由式(8)计算得到。具体流程见图 2。

图 2

图 2 粒子群优化算法流程图 Fig. 2 The flow chart of particle swarm optimization algorithm

$ \left\{ {\begin{array}{*{20}{c}} {{v^{ij}}\left( {t + 1} \right) = \omega \cdot {v^{ij}}\left( t \right) + A + B}\\ {A = {c_1} \cdot r_1^{ij}\left( t \right) \cdot \left( {x_{{\rm{pbest}}}^{ij}\left( t \right) - {x^{ij}}\left( t \right)} \right)}\\ {B = {c_2} \cdot r_2^{ij}\left( t \right) \cdot \left( {x_{{\rm{gbest}}}^j\left( t \right) - {x^{ij}}\left( t \right)} \right)} \end{array}} \right. $

(7)

$ {x^{ij}}\left( {t + 1} \right) = {x^{ij}}\left( t \right) + {v^{ij}}\left( {t + 1} \right) $

(8)

式中：i为第i个粒子，j为第j个维度；t为更新次数；v为粒子速度，x为粒子位置；r₁、r₂为0到1均匀分布的随机数；x_pbest粒子自身最好位置，x_gbest全局最好位置；ω为惯性权重值。

1.4 基于粒子群优化算法的克里金插值法

拟合半变异函数的模型有很多种，但无论是哪种拟合模型，都涉及到3个拟合参数，分别是图 1中的核(nugget)、变程(range)和梁(sill)。在通常情况下，无论是普通克里金法还是协同克里金法，在确定上述3个参数时用的均为最小二乘法，即在满足方差最小的条件下得到的，这就会使得3个参数存在较大的误差。利用粒子群优化算法能快速寻找到全局最优解的特点，将该优化算法应用到确定半变异模型3个参数上。其中该优化算法中的适应度函数为插值结果的均方根误差。

同时，大量研究表明，PSO算法易陷入局部最优和早熟收敛等缺陷。采用一种基于高斯变异的方法来提高种群的多样性^[15]。在算法出现过早收敛时，能够使粒子在其他区域进行搜索，跳出局部最优，寻找更优的解。即在迭代到一半次数后，开始对粒子进行变异，对每个粒子进行高斯变异^[18]。公式为

$ {x_d} = {\rm{gbest}}\left( d \right) \times \left( {0.5 + \sigma } \right) $

(9)

式中：gbest(d)为全局最优在d维的值；σ为高斯白噪声。

1.5 预测精度评价方法

为了能够对上述方法的空间插值精度进行比较, 使各方法间的结果更具有可比性, 可通过计算钻孔点的勘察数据与预测数据的误差来评估各种方法的优劣^[19]。均方根误差(RMSE)可以用来评价预测值和真实值之间的接近程度。利用协同克里金法与普通克里金的均方根误差减少的百分数(RRMSE)来表示预测精度的提高程度。

$ {\rm{RMSE = }}\sqrt {\left\{ {\sum\limits_{i = 1}^N {{{\left[ {Z\left( {{x_i}} \right) - {Z^*}\left( {{x_i}} \right)} \right]}^2}} } \right\}/N} $

(10)

$ {\rm{RRMSE = }}\frac{{{\rm{RMS}}{{\rm{E}}_{{\rm{OK}}}} - {\rm{RMS}}{{\rm{E}}_{{\rm{COK}}}}}}{{{\rm{RMS}}{{\rm{E}}_{{\rm{OK}}}}}} $

(11)

式中：Z(x_i)与Z^*(x_i)分别为x_i处的勘察值与估计值；N为样本数量。

2 工程实例

2.1 研究区域概况与数据处理

2.1.1 研究区域概况

试验区位于重庆万盛经济开发区，地貌单元为丘陵斜坡地貌；总体趋势为西侧高，东侧低；地表总体坡度2°~34°。上图中的红色部分为钻孔点的位置，已知数据有勘察区的平面图、剖面图及钻孔柱状图。试验区最高海拔为418.23 m，最低海拔为334.77 m；共有124个钻孔点。具体情况见图 3。

图 3

图 3 试验区钻孔分布图 Fig. 3 The distribution of borehole

2.1.2 钻孔点数据处理与分析

对试验区124个钻孔的高程数据和土层厚度分析，结果显示：土层厚度和钻孔点的高程有一定的规律(为了便于与土层厚度比较，高程基准点设为320 m)，具体表现为:在一定的高程范围内，钻孔的土层厚度大致和高程值成正相关(见图 4)。

图 4

图 4 土层厚度和高程关系图 Fig. 4 The correlation between soil thickness and altitude

协同克里金法中，作为辅助变量的前提条件是该变量与待估计变量之间存在着相关性。用SPSS统计软件中的数据分析功能得到高程值和土层厚度的相关系数为0.552，且在0.01的显著性水平的条件下，相关性程度为显著性相关(表 1)。因此，可将高程值作为提高土层厚度插值精度的辅助变量。

表 1

表 1 SPSS统计软件计算结果 Table 1 The result calculated by SPSS statistical software 高程 土层厚度 高程 Pearson相关性 1 0.552** 显著性(双侧)   0.000 N 124 124 土层厚度 Pearson相关性   1   0.552** Pearson相关性   0.552**   1 显著性(双侧) 0.000 N 124 124 注：**水平(双侧)上显著相关。

表 1 SPSS统计软件计算结果 Table 1 The result calculated by SPSS statistical software

2.2 空间插值模型试验

采用交叉验证的方法^[20]来估计土层厚度，并用数理统计的方法来比较不同估计模型的估值精度，即首先将待预测点的属性值Z(x_i)暂时剔除，然后将最小二乘法得到的半变异函数表达式用于普通克里金法和考虑高程的协同克里金法来预测Z(x_i)的值。结合相应的公式并利用MATLAB编制程序求出拉克朗日系数和待估计点的影响范围内的各个钻孔点的权重系数，最后利用式(3)与式(6)分别求出待估计点的属性值Z^*(x_i)。重复以上过程直到对所有的124个观测点进行估计。

基于粒子群优化算法的克里金插值法与基于最小二乘法的克里金插值法区别在于：粒子群优化算法中的半变异模型中的3个参数是未知的；随机的给每个粒子赋值，这样每一个粒子就对应了一个半变异函数的模型；然后将每个粒子所对应的半变异函数模型，带入用MATLAB编制的相应插值程序中，求出均方根误差，并且将该均方根误差作为该粒子的适应度函数的适应值。当循环结束后，最小的适应度值的粒子即为所求。进而求出了半变异函数的表达式。

无论是哪一种方法, 每一个钻孔点可以得到真实值和估计值两类数据。将各个插值方法得到的插值结果与真实值做成火柴棍图，见图 5~图 8。

图 5

图 5 最小二乘法的O-Kriging插值结果与真实土层厚度 Fig. 5 O-Kriging based on the least squares method interpolation results and the true soil thickness

图 6

图 6 粒子群优化算法的O-Kriging插值结果与真实土层厚度 Fig. 6 O-Kriging based on PSO method interpolation results and the true soil thickness

图 7

图 7 最小二乘法的CO-Kriging插值结果与真实土层厚度 Fig. 7 CO-Kriging based on the least squares method interpolation results and the true soil thickness

图 8

图 8 粒子群优化算法的Co-Kriging插值结果与真实土层厚度 Fig. 8 Co-Kriging based on PSO method interpolation results and the true soil thickness

从图 5~图 8中可以看出，在上述的4类方法中，基于最小二乘法的普通克里金插值法与勘察值之间的误差是最大的；而考虑高程的协同克里金法的插值的精度要高于基于粒子群优化算法的普通克里金法；基于粒子群优化算法的协同克里金插值法的精度在4类方法中是最高的，除了极个别的土层厚度较厚的钻孔点外，其他钻孔点均与实际情况吻合得较好。

用数理统计的方法来比较4种插值结果的精度的提高程度。其结果见表 2。

表 2

表 2 O-Kriging与Co-Kriging预测精度的比较 Table 2 Error analysis of O-Kriging and Co-Kriging 插值方法 变量 样本数 RRMSE/% 基于最小二乘法的O-Kriging 土层厚度 124 基于粒子群优化法的O-Kriging 土层厚度 124 28.79 基于最小二乘法的CO-Kriging 土层厚度、高程 124 39.32 基于粒子群优化法的CO-Kriging 土层厚度、高程 124 48.45

表 2 O-Kriging与Co-Kriging预测精度的比较 Table 2 Error analysis of O-Kriging and Co-Kriging

从表 2可以看出，当采用普通克里金插值法对边坡的土层厚度进行空间插值时，基于粒子群优化算法的普通克里金插值法较基于最小二乘法的普通克里金插值法的均方根误差减少了28.79%；当均采用最小二乘法去拟合半变异函数时；协同克里金法的均方根误差较普通克里金法减少了39.32%；当都采用PSO法去拟合模型参数时；协同克里金法的均方根误差较普通克里金法减少了19.66%；同时还可以从表中得出结论：当考率单一因素时，考虑高程的协同克里金法较基于粒子群优化算法的普通克里金插值法的插值精度提高了10.53%。

3 结论

以重庆某坡地为研究对象，利用该区域勘察报告中的124个钻孔点的勘察数据为试验样本，采用交叉验证的方法分别用普通克里金插值法和考虑高程的协同克里金插值法对土层厚度进行虚拟差值。还对在拟合半变异函数模型参数时存在的误差进行了研究，用粒子群优化算法拟合半变异函数模型，并对4种不同方法的插值结果进行了对比分析，得到的结论如下：

1) 将最小二乘法拟合的半变异函数用于普通克里金法插值时会产生较大误差。

2) 对土层厚度进行克里金插值时，考虑高程的协同克里金法较普通的克里金法的插值精度高。

3) 通过交叉验证表明，粒子群优化算法能在一定程度上提高普通克里金插值法和协同克里金插值法的插值精度。

参考文献

[1]	BURROUGH P A, MCDONNELL R A. Principle of geographic information systems[J]. Oxford University Press, 1998, 12(1): 102–102.
[2]	闫庆武, 卞正富, 王红. 利用泰森多边形和格网平滑的人口密度空间化研究——以徐州市为例[J]. 武汉大学学报(信息科学版), 2011(8): 987–990+1010-1011. YAN Q W, BIAN Z F, WANG H. The population density by using tyson polygon-taking Xuzhou city as an example[J]. Geomatics and Information Science of Wuhan University, 2011(8): 987–990+1010-1011. (in Chinese)
[3]	OLIVER M A, WEBSTER R. Kriging:A method of interpolation for geographical information systems[J]. International Journal of Geographical Information Systems, 1990, 4(3): 313–332. DOI:10.1080/02693799008941549
[4]	MESNARD L D. Pollution models and inverse distance weighting:Some critical remarks[J]. Computers & Geosciences, 2013, 52(1): 459–469.
[5]	胡绪福, 唐玮, 王博, 等. 趋势面分析法应用效果研究——以YL油田某区块为例[J]. 长江大学学报(自然科学版理工卷), 2010(3): 536–538. HU X F, TANG W, WANG B, et al. The Application of trend surface analysis method-taking a block in YL oilfield as an example[J]. Journal of Yangtze University (Natural Science Edition), 2010(3): 536–538. (in Chinese)
[6]	刘小飞, 张寄阳, 刘祖贵, 等. 喷灌大田土壤水分不同空间插值方法对比分析[J]. 灌溉排水学报, 2008(4): 116–118. LIU X F, ZHANG J Y, LIU Z G, et al. Comparative analysis of different spatial interpolation in sprinkler irrigation field[J]. Journal of Irrigation and Drainage, 2008(4): 116–118. (in Chinese)
[7]	YATES S R, WARRICK A W. Estimating soil water content using cokriging[J]. Soil Science Society of America Journal, 1987, 51(1): 23–30. DOI:10.2136/sssaj1987.03615995005100010005x
[8]	GHADERMAZI J, SAYYAD G, MOHAMMADI J, et al. Spatial prediction of nitrate concentration in drinking water using pH as auxiliary co-kriging variable[J]. Procedia Environmental Sciences, 2011, 3: 130–135. DOI:10.1016/j.proenv.2011.02.023
[9]	胡丹桂, 舒红. 基于协同克里金空气湿度空间插值研究[J]. 湖北农业科学, 2014(9): 2045–2049. HU D G, SHU H. Air humidity based on Co-kriging interpolation[J]. Hubei Agricultural Sciences, 2014(9): 2045–2049. DOI:10.3969/j.issn.0439-8114.2014.09.018 (in Chinese)
[10]	许晶玉.山东省种植区地下水硝态氮污染空间变异及分布规律研究[D].长沙: 中南大学, 2011. XU J Y. Spatial variability and distribution pattern of groundwater nitrate pollution in farming regions of Shandong Province[D].Changsha: Central South University, 2011.(in Chinese)
[11]	杜文凤, 彭苏萍. 利用地质统计学预测煤层厚度[J]. 岩石力学与工程学报, 2010(Sup1): 2762–2767. DU W F, PENG S P. Coal seam thickness prediction with geostatistics[J]. Chinese Journal of Rock Mechanics and Engineering, 2010(Sup1): 2762–2767. (in Chinese)
[12]	黄大年, 高秀鹤, 孙思源, 等. 重力梯度数据协克里金三维反演确定岩脉倾向[J]. 吉林大学学报(地球科学版), 2016: 1–8. HUANG D N, GAO X H, SUN S Y, et al. Identify the dip angle of the dipping dike model based on co-kriging inversion of gravity gradient data[J]. Journal of JiLin University, 2016: 1–8. (in Chinese)
[13]	严华雯, 吴健平. 加权最小二乘法改进遗传克里金插值方法研究[J]. 计算机技术与发展, 2012(3): 92–95. YAN W H, WU J P. Research on genetic algorithm kriging optimized by weight least square[J]. Computer Technology and Development, 2012(3): 92–95. DOI:10.3969/j.issn.1673-629X.2012.03.024 (in Chinese)
[14]	张强, 许少华, 于文涛, 等. 粒子群算法在克里金三维地质建模中的应用[J]. 大庆石油学院学报, 2011(1): 85–89, 120. ZHANG Q, XU S H, YU W T, et al. The particle swarm is applied in kerrin 3D modeling[J]. Journal of Daqing Petroleum Institute, 2011(1): 85–89, 120. DOI:10.3969/j.issn.2095-4107.2011.01.018 (in Chinese)
[15]	贾雨, 邓世武, 姚兴苗, 等. 基于约束粒子群优化的克里金插值算法[J]. 成都理工大学学报(自然科学版), 2015(1): 104–109. JIA Y, DENG S W, YAO X M, et al. Kriging interpolation algorithm based on constraint particle swarm optimization[J]. Journal of Chengdu University of Technology (Science & Technology Edition), 2015(1): 104–109. DOI:10.3969/j.issn.1671-9727.2015.01.13 (in Chinese)
[16]	KENNEDY J, EBERHART R. Particle swarm optimization[C]//Proc IEEE Int Conf on Neural Networks, Perth, 1995: 1942-1948.
[17]	EBERHART R, KENNEDY J. A new optimizer using particle swarm theory[C]//Proc 6th Int Symposium on Micro Machine and Human Science, Nagoya, 1995: 39-43.
[18]	KROHLING R A. Gaussian swarm: a novel particle swarm optimization algorithm[C]//Cybernetics and Intelligent Systems, 2004 IEEE Conference on IEEE, 2005: 372-376.
[19]	史文娇, 岳天祥, 石晓丽, 等. 土壤连续属性空间插值方法及其精度的研究进展[J]. 自然资源学报, 2012(1): 163–175. SHI W J, YUE T X, SHI X L, et al. The accuracy of soil Interpolation[J]. Journal of Natural Resources, 2012(1): 163–175. (in Chinese)
[20]	GOOVAERTS P. Geostatistics for natural resources evaluation[M]. Oxford: Oxford University Press, 1997: 437-438.

X