如何由离散的单一沉降值推演得到区域沉降趋势面，是进行建筑物荷载作用下地面沉降数据建模与可视化空间分析的基础^[1]。地面沉降趋势面空间插值，就是选择具有代表性的沉降观测点，通过一定的插值算法来模拟整个沉降区域的地表形状。国内外学者对此进行了研究，例如，文献[2-4]对已有空间插值模型进行比较，并在此基础上修订相关的模型和方法，提高区域沉降的插值精度；文献[5-8]将样条函数、平行Kriging等插值方法应用于地学空间模拟的计算和分析中。这些方法主要的研究对象是传统测量方法获得的地面沉降点数据，而对于建筑物荷载引发地面沉降的空间模拟则相对较少。本文根据研究区不同建筑物的荷载情况以及地下岩土特征，对普通Kriging插值方法进行优化，并进行城市地面沉降空间趋势面模拟和计算。

1 理论基础 1.1 单体建筑物沉降量计算方法

建筑物自身荷载产生的地面沉降为^[9]：

(1)

式中，s为单体建筑物沉降量，ψ_s为经验系数，s_i为各土层沉降量，p₀为桩端附加压力，n为土层数，E_si为土层压缩模量，z_i为第i层土厚度。

建筑物荷载引起的地下土体和水体的应力变化在地下土体应力场和地下水体渗流场所构成的地下三维空间中是连续的^[10]，而通过式(1)计算的每栋建筑物荷载所产生的沉降值是离散数据。需要将离散点的沉降量通过空间插值拟合成连续变化的空间趋势面，在地下空间进行组织和重构，从而进行地面沉降的可视化分析。

1.2 应用普通Kriging插值模型模拟地面沉降

普通Kriging插值采用半方差函数作为确定权重参数λ_i的依据^[11]。设x₁, x₂, …, x_n为某区域若干个地面沉降观测点，z(x₁), z(x₂), …, z(x_n)为相应的地面沉降量。根据普通Kriging插值方法计算得到的区域化沉降变量在某一点x₀处的沉降值z^*(x₀)采用一个线性组合来估计：

(2)

进一步推导，可得到普通Kriging插值的公式：

(3)

式中，λ_i为权重系数，μ为拉格朗日系数，γ(x_i-x_j)为变异函数，j=1, …, n。

1.3 改进Kriging插值的构建

普通Kriging插值没有将建筑物荷载考虑进去，更没有解决地下土体受力变形各向异性的问题。我们将普通Kriging插值方法中的z(x_i)用式(1)中的单体建筑物沉降量s替代，再叠加上地下水抽取引发的沉降量s′，从而得到区域内的综合沉降数值：

(4)

式中，z^*(x₀)为x₀处的沉降量，λ_i为权重系数，p₀为桩端附加压力，m为土层数，E_si为土层压缩模量，z_i为第i层土的厚度，s′为地下水抽取引发的沉降量。

2 改进Kriging插值的应用与验证 2.1 研究区与数据源

本文的研究区选取天津市东南部沉降区(117°30′12″~118°54′30″E，38°12′18″~40°42′24″N)，属于冲积海积平原。该区域地下土压缩层厚度在60~80 m，远离天津市地下水开采的第四系承压含水层中埋深170~190 m的第2含水组^[12]，地下水开采引发的沉降在这一区间相对较弱，因此研究重点是建筑物荷载所产生的地面沉降。

计算使用ArcGIS10.0软件。首先根据研究区域内均匀分布的48个水文地质钻孔点的勘探数据得到地下土壤的分层性质(包括分层厚度及各层的压缩模量)以及抽水带来的沉降数值，用普通Kriging方法插值到整个研究区范围；然后根据矢量化得到的研究区内2001~2015年共3 580栋建筑物的平面形状，提取其质心作为每栋建筑物的空间分布位置，并采用掩膜方法得到建筑物下面的土层厚度和压缩模量；再分别计算每栋建筑物的荷载大小，得到其桩端附加压力，以此构建沉降分析的空间数据库，并将上述空间变量代入改进的Kriging插值模型；最后再通过研究区内InSAR数据反映的地面沉降趋势对插值结果进行比较和验证。

2.2 改进Kriging插值模型的应用 2.2.1 选择变异函数模型

由式(1)可知，每栋建筑物之间的距离和沉降量之间的函数关系都与其桩端附加压力、地基土的土层厚度和压缩模量有关，而这3个参数在每栋建筑物不同的空间位置上都是变化的。所以，改进Kriging插值的变异函数也发生改变。

改进Kriging插值变异函数中距离与半方差之间的关系可以是线性关系，也可以是二次函数、三次函数、指数或对数关系等非线性关系。为了确认这种关系，首先计算空间任意两点(x_i, y_i)与(x_j, y_j)之间的距离：

(5)

再根据上述变化的3个参数计算半方差：

(6)

从而得到n²个[D_{(i, j)}，γ(x)]数据对。再将所有沉降点的D_{(i, j)}和γ(x)绘制成散点图，导入统计分析软件SPSS，寻找一个最优函数曲线来拟合这种关系，从而得到函数表达式：

(7)

为了简化计算，先选取研究区2011~2015最近5 a新建建筑物数量(1 460栋)的70%，即1 022栋建筑物的年均沉降量作为样本点进行拟合，余下的30%作为验证点。最终确定最优函数为三次项拟合函数(图 1)，也就是通常所说的球状模型，其R²等于0.816，大于普通Kriging插值的指数函数和其他函数(表 1)。

2.2.2 计算权重系数

将上面的球状模型代入变异函数模型，得：

(8)

式中，c为基台值，α为变程，h=x_i-x_j为滞后距。

在接近原点处，变异函数呈线性形状，在变程α处达到基台值c。合并式(3)和(8)，得：

(9)

由式(9)可知，如果知道γ(x_i-x_j)=γ_{(i, n)}的值，由变异函数与权重系数的函数关系式:

(10)

对式中的γ_in求逆，再与γ_n0相乘，就可以计算出改进Kriging插值模型的权重系数λ_i^[11]。

2.2.3 区域地面沉降的空间插值模拟

在已知每栋建筑物桩端附加压力p₀、建筑物地基中所划分的土层数m、各层土的厚度z_i和压缩模量E_si的情况下，根据式(10)计算出来的权重λ_i，代入式(4)即可采用改进Kriging插值模型模拟区域的地面沉降，从而得到研究区建筑物荷载作用下的地面沉降趋势面。本文使用MATLAB语言实现改进后的Kriging插值模型，对研究区2011~2015年1 022栋建筑物的年均沉降量进行实验，其计算结果如图 2所示。

2.3 插值结果检验

使用研究区2011~2015新建建筑物的30%，即438栋建筑物作为验证点进行检验。首先，将改进Kriging插值模型运算后的插值结果直接从MATLAB拟合的空间曲面中提取出来；然后，使用交叉验证的方法来验证插值方法的准确程度，以此来评判其精度。交叉验证的主要精度指标包括平均绝对误差MAE、平均相对误差MRE和均方根误差RMSE。其中，MAE反映插值结果可能的误差范围，MRE反映采样点数据的估值灵敏度和极值效应，RMSE反映插值结果的离散程度^[14]。

分别对反距离加权插值(IDW插值)、普通Kriging插值、改进Kriging插值方法计算其MAE、MRE和RMSE，并对改进Kriging插值的效果进行评价(表 2)。

可以看出3种插值方法的平均绝对误差，IDW插值和普通Kriging插值都大于2 mm，而改进Kriging插值只有1.425 mm；对于平均相对误差，由于IDW插值属于确定性的插值法，插值的表面都会通过采样点，采样点的取值变化对插值结果会有一定影响，取值变化相对比较均衡的区域，插值结果会比较准确，而沉降变化比较剧烈和频繁的区域，其插值结果差异就较大。因此，如果要想利用IDW插值法模拟沉降趋势，就必须要求研究区域样本点的取值是均匀变化、不剧烈的^[14]，这显然与研究区地下土体受力变形各向异性造成沉降数值大小不一的实际情况不符，故其插值精度比较低。而普通Kriging插值和改进Kriging插值因为是基于统计学的插值法，都是根据采样点在空间位置上的分布进行结构分析，并确定对待计算点有影响的距离范围，再根据与之相关的其他观测点的沉降值进行插值，因此Kriging插值得到的权重较符合实际地形变化。在本例中，普通Kriging插值和改进Kriging插值的MRE均在0.2左右，改进Kriging插值的精度要更高些，为0.197。最后，对比3者的均方根误差可见，改进Kriging插值的总体离散程度明显小于IDW插值和普通Kriging插值，后二者的RMSE数值分别为2.398 mm和2.167 mm，而改进Kriging插值只有1.352 mm。

综合3种方法来看，IDW插值方法的精度最低，而普通Kriging插值和改进Kriging插值相对较好，特别是改进Kriging插值方法，3项评价精度指标均为最小，对区域地面沉降的趋势面拟合较好，是一种比较理想的插值方法。根据上述改进的Kriging插值模型，分别对研究区2001~2005年、2006~2010年、2011~2015年3个时间段的全部建筑物年均沉降量进行插值运算，并转入ArcGIS中叠加显示(图 3)。由图可见，改进Kriging插值较好地反映了沉降量受不同建筑物荷载产生附加应力以及由地基土压缩性质不同而表现出的空间变异性，不但有效地消除了IDW插值的“牛眼”效应和普通Kriging插值的平滑效应，而且其沉降趋势面也最接近于同一时期InSAR数据所反映的年均沉降量，更真实地反映了地面变形趋势。

以2001~2005年、2006~2010年、2011~2015年这3个时间段InSAR数据所反映的年均沉降量作为该时期研究区沉降的真实值，将不同插值方法的结果与之作差，并计算平均绝对误差MAE来进行量化对比，发现改进Kriging插值的MAE值为1.428 mm，而普通Kriging插值的MAE值为1.743 mm，IDW插值的MAE值为3.269 mm。说明改进Kriging插值的精度要高于普通Kriging插值和IDW插值，这与图 3反映的沉降趋势是吻合的。

3 结语

1) 相对于IDW插值方法，Kriging插值模型属于地统计学的范畴，可以解决估算点空间各向异性的问题，得到的结果也较符合实际地形变化。

2) 改进Kriging插值模型进一步优化了普通Kriging插值对观测点权重影响的处理方法，既考虑了其他观测点方向、距离的影响，又具有对建筑物荷载产生的附加应力以及土壤压缩性质的相关性响应机制，从而能更好地反映出建筑物荷载作用下地面沉降的实际空间分布特征，消除估算点与周围区域地基土中附加应力以及压缩模量差异而引起的插值异常。

3) 实际应用表明，改进Kriging插值的精度高于普通Kriging插值方法，计算得到的沉降趋势面也最接近于真实沉降变形趋势面，适合于地面沉降的空间模拟与分析。

需要指出的是，本文是在相对较小的区域内，对建筑物荷载引发的地面沉降进行分析。如果区域进一步扩大，则应将地壳形变、地质构造变异等其他因素考虑进来，进一步优化推演模型，从而更加准确地预测地面沉降的变形趋势。