克里金插值算法在形变监测中的应用

刘 洪

（1.广西有色勘察设计研究院，广西 南宁 530031）

地表的垂向形变监测在生产生活以及科学研究中均发挥重要作用，如城市沉降监测、地震形变、火山形变等研究均需要高精度的地表垂直形变场［1-5］。高精度的地表垂直形变场在城市沉降监测中能提供直观的空间信息，故可以采用GPS台站来进行监测研究地表形变［6］。其常规做法是，在成本范围内，以一定距离均匀或随机地布设GPS台站，观测并计算各台站位置处的垂直形变后，利用克里金插值算法将其插值为地表垂直形变场［7］。此类方法中，GPS台站的分布等直接决定其布设成本和结果精度，故需要进行比较分析。本文拟构建两种不同的GPS台站分布：均匀布设和随机布设，分析GPS台站在克里金插值算法时的优缺点。

1 克里金插值算法

假设已有N个GPS台站观测得到数据，z=[z1,z2,z3,…,zN]T，待求解位置为待估点z0。根据无偏估计和方差最小，可求得估计值为：

式中，ωi为待估加权系数，满足条件：

根据方差最小的原则，可以得到：

式中，γ(zi,zj)和 γ(zi,z0)分别为是 zi与 zj，zi与 z0之间的变差函数。构造函数：

式中，ϕ为待定的拉格朗日算子，可导出插值的解应满足：

将其写成矩阵形式：

求得ωi值和ϕ值后，便可通过式（1）求出待估位置z0处的最优估值。

2 模型构建

图1 已知地表垂直形变场

图2 GPS台站均匀布设时克里金插值得到的垂直形变场结果（案例I）

假设地表形变如图1所示，该形变范围为100 km× 100 km，有5个主要形变区域：A、B、C、D和E。其中，A、B和D为地表抬升区域，C和E为地表下降区域。假设100个GPS台站，其均匀布设时台站（图2中黑色圆点）间距为10×10 km，较符合GPS台站成本。本文模拟假设随机布设分布为图3、4中黑色圆点所示，首先利用假设的已知形变提取出GPS台站处的形变观测值；本文模拟数据噪声均值为0，标准差为5 mm的高斯白噪声，并加在GPS台站的形变观测值中。

图3 GPS台站随机布设时克里金插值得到的垂直形变场结果（案例II）

图4 GPS台站随机布设时克里金插值得到的垂直形变场结果（案例III）

3 计算结果与分析

图2为均匀布设GPS台站时，利用克里金插值算法得到的形变场结果，及其与已知地表形变场的差异，黑色圆点表示均匀分布的100个GPS台站位置。结果显示，3个形变抬升区域与2个形变下降区域（图2a）均可通过克里金插值算法得到，然而其形状与假设已知形变有一些差异。图3、4为随机布设GPS台站时，利用克里金插值算法得到的形变场及其与已知形变场的差异。

均匀布设时，利用克里金插值算法得到的形变场与已知假设形变场之间的差异（图2b）最大值位于D区域，这是由于D区域为抬升区域，处于下降区域（C和E）中间位置，故其形变梯度较大。此问题在随机布设案例中（图3、4）也有体现。

随机布设案例I中的最大差异（图3b）位于A区域的东侧，该位置形变量和形变梯度均较大，但是GPS台站稀疏，故克里金插值的结果较差。随机布设案例II中的最大差异（图4b）同样位于形变梯度较大且GPS台站稀疏的位置。

均匀和随机布设案例中，在监测区域（100 km×100 km）的边缘位置，利用克里金插值算法得到的形变均较差，这是由于这些位置GPS台站稀疏导致。

对比均匀和随机布设案例发现，均匀布设GPS台站时，使用克里金插值的结果与已知形变之间差异的平均值和标准差较小（表1），其最大差异-1.809 9 cm位于监测区域的边缘位置。故均匀布设相对更优。在形变梯度较大位置，GPS台站需要加密布设［8］。

表1 均匀与随机布设的GPS台站对插值结果的影响

4 结 语

比较了均匀布设和随机布设GPS台站时，利用克里金插值算法得到的地表垂向形变场差异。利用均匀布设的GPS台站数据进行插值的效果整体优于随机布设；而在形变量较大或形变梯度较大区域，随机布设时由于分布有更多GPS台站，效果更优。未来形变监测工作的GPS台站布设可考虑大范围内均匀布设，主要形变区加密均匀布设，在形变量较大或形变梯度较大区域，则需加密GPS台站布设。

[1] 周琦,赵文吉．基于 D-InSAR 三轨法的北京平原区地表形变监测[J]． 地理空间信息,2013(2):146-148

[2] 路旭,匡绍君,贾有良,等． 用 InSAR 作地面沉降监测的试验研究[J]． 大地测量与地球动力学,2002,22(4):66-70

[3] 刘晓,余洁,刘利敏．利用 PSInSAR 技术监测北京地面沉降[J]．地理空间信息,2016,14(1):86-88

[4] 余景波,刘国林,曹振坦,等． D-InSAR 在提取地面形变信息中的应用[J]．地理空间信息,2011(6):4

[5] 胡亚轩,施行觉,王庆良,等．腾冲火山区地表垂直形变分析[J]．大地测量与地球动力学,2003,23(2):37-41

[6] 黄立人,匡绍君,杨贵业．GPS 观测得到的天津地区的现今变形[J] 大地测量与地球动力学,2002,22(4):17-20

[7] Gudmundsson S, Sigmundsson F, Carstensen J M． Three-Dimensional Surface Motion Maps Estimated From Combined Interferometric Synthetic Aperture Radar and GPS Data[J]．Journal of Geophysical Research: Solid Earth, 2002（B10）:1-14

[8] 井研,胡庆武,艾明耀,等． 基于离散沉降观测的城市地面沉降监测[J]． 地理空间信息,2016(11):84-87