基于GPS-InSAR加权融合的矿区地表形变监测研究

王志岗，周文韬

（1.四川科技职工大学 应急管理系，四川 成都 610000；2.西南科技大学 环境与资源学院，四川 绵阳 621010；3.国家遥感中心绵阳科技城分部，四川 绵阳 621010）

矿物被开采后，岩层内部应力的平衡状态遭到破坏，导致地表产生形变并出现塌陷坑[1]。地表塌陷对矿区地下安全开采构成威胁，也给矿区周围环境、建筑等造成了破坏[2]。为保护地面建筑并保证地下开采工作正常进行，监测地表形变十分必要。目前，全球导航定位系统（GPS）技术、合成孔径雷达干涉测量（InSAR）技术已被广泛应用于矿区地表形变监测中，GPS技术具有高时间分辨率和高精度[3]，InSAR技术具有监测范围广、不受天气影响等优点[4]；但这两种技术在实际应用中均有不足，即GPS技术仅能反映离散点的监测结果，InSAR技术仅能反映视线向的相对形变量，若仅将视线向形变投影到垂直方向，该形变往往不能客观反映实际形变。王霞迎[5]、曹海坤[6]等基于升降轨InSAR与GPS数据，建立了融合函数模型反演高精度三维形变场；马艳鸽[7]、丁宁[8]和JI P F[9]等基于先验算法改进融合GPS与InSAR的函数模型，以提高形变监测精度。众多学者均是通过建立融合函数模型，利用GPS数据校正InSAR数据来提高InSAR形变精度[10-12]。然而，当地形情况复杂时（如矿区、冰川等），GPS点的形变不能反映其周围面的形变，具有较高的偶然性；InSAR技术的地面分辨率较低，却能很好地反映面的形变。因此，如何通过加权融合，将降低GPS点偶然性和提高InSAR形变监测精度进行统一成为亟待解决的问题。

鉴于此，本文重点研究地表垂直方向的形变，介绍了不同内插GPS点的形变结果与优势，再利用变异系数法（CV）得到组合内插GPS形变结果，最后运用算术平均法将组合内插GPS形变结果与SBAS-InSAR垂直向累积形变结果进行加权融合，充分发挥二者优势，得到GPS-InSAR形变结果。

1 理论基础

1.1 空间插值方法

经验贝叶斯克里金插值法（EBK）是一种地统计插值方法，可通过构造子集和模拟的方式自动计算构建有效克里金模型过程中的参数。与普通克里金插值[13]算法不同，EBK法考虑了半变异函数估计的不确定性，可通过估计半变异函数来说明引入的误差，因此EBK法降低了预测的标准误差[14]。

局部多项式插值法（LPI）是一种局部加权最小二乘拟合法，根据有限的监测数据，采用多个多项式来拟合表面。它引入了“距离权”的概念，对于未知点的计算，考虑在局部范围内所有已知点对其的贡献，距未知点近的点权重大，反之则权重小。每个未知点的预测值都对应一个多项式，每个多项式都处于特定重叠的邻近区域内，通过最小二乘法求解邻域内多项式组成的方程组，从而得到拟合表面[15]。该方法不仅具有趋势面法考虑全部数据点反映趋势性变化的优点，而且具有距离法反映局部特征的优点。

径向基函数插值法（RBF）是一种人工神经网络方法，根据有限的监测数据，选择合适的径向基函数生成一个具有最小曲率、且到各样点的Z值距离最小的曲面。该方法拟合的表面经过所有点数据，并可计算得到高于或低于点Z值的预测值，适用于监测点数据集大、表面变化平缓的情况[16]。

1.2 CV法

CV法是一种客观计算权重的方法，可直接利用各项评价因子包含的信息，通过计算得到各项指标的权重。为消除各指标量纲不同带来的影响，需利用各指标的变异系数来衡量其取值的差异程度[17]。其计算公式为：

式中，Vi为第i项指标的变异系数，亦称标准差系数；Si为第i项指标的标准差；为第i项指标的平均值。

根据式（1）可以得到各项指标的权重，即

1.3 精度评定指标

本文通过比较GPS点预测值与实测值的均方根误差（RMSE）、平均绝对误差（MAE）和拟合优度（R2）3个指标来评定插值精度。其计算公式为：

式中，xi为第i点的实测值；为第i点的预测值；为实测值x的平均值。RMSE和MAE值越小、R2越趋近1，表示预测精度越高。

2 矿区地表形变监测实验分析

本文以金川某矿西二采区为研究区（以下简称西二采），位于我国甘肃省金昌市中部。矿区地势较平坦，平均海拔约为1 500～1 800 m；地表水系不发育，常年干枯，属温带大陆性气候[18]。近年来，由于地下不断开采，导致地表塌陷、裂缝明显，为铜镍矿资源开采留下了极大的安全隐患。

2.1 数据基础

在研究区范围内，布设GPS点进行静态观测，监测时间为2019年4月—2020年6月，监测周期为60 d，共进行8期监测。为与GPS监测时间统一，实验选取覆盖研究区的38景C波段Sentinel-1A升轨影像，时间跨度为2019年3月22日—2020年6月8日，设置临界基线阈值为45%，时间基线阈值为120 d，生成325个干涉对，利用小基线（SBAS）方法提取垂直向累积形变结果。Sentinel-1A升轨影像覆盖范围和GPS点位分布情况如图1所示。

图1 西二采影像覆盖范围和GPS点位分布图

2.2 地表形变监测实验与分析

2.2.1 基于SBAS-InSAR技术的地表形变监测

本文利用SBAS-InSAR技术获取了研究区2019年3月22日—2020年6月8日的地表垂直向累积形变量，如图2所示，可以看出，矿区地表沉降区域主要分布在矿体西南部，该区域受地下开采影响较大；在15个月的时间里，矿区内最大沉降量达到-168.41 mm，地表整体呈下降趋势，沉降面积约为0.267 km2。

图2 SBAS-InSAR垂直向累积形变量

2.2.2 基于CV的三插GPS形变实验与分析

为融合GPS数据与InSAR数据，需先将GPS点元数据内插为面元数据。EBK法克服了传统克里金法需人机交互式建模来寻找准确结果的缺点，提高了预测精度，并可准确预测不够稳定的数据；但处理时间随监测点数、子集大小或重叠系数的增加而快速增加，降低了计算效率，且对于含有异常值的数据，可能会得到大于或小于输入点值若干个数量级的预测结果，从而影响插值结果。因此，本文结合两种精确插值方法来提高运算效率，降低异常数据对插值结果的影响。LPI法突出点数据反映趋势性变化，同时通过距离法强调局部特征。RBF法可预测大于最大测量值和小于最小测量值的值，适用于表面变化平缓的情况。上述3种插值方法的结果如图3～5所示，其中EBK法插值后的最大累积沉降量为-320.88 mm，最大抬升量为123.47 mm；LPI法插值后的最大累积沉降量为-428.49 mm，最大抬升量为565.24 mm；RBF法插值后的最大累积沉降量为-368.34 mm，最大抬升量为264.89 mm；说明3种插值方法内插后的面数据整体形变趋势较一致，3～8行的沉降量大，9～12行出现轻微抬升。本文综合3种插值方法的优点，根据式（2）得到EBK-LPI-RBF法的插值结果，如图6所示。3种插值方法的权重分配如表1所示。

表1 CV法权重分配

图3 EBK法插值结果

图4 LPI法插值结果

图5 RBF法插值结果

图6 EBK-LPI-RBF法插值结果

不同插值方法的精度比较如表2所示，可以看出，LPI法和RBF法的精度略高于EBK法，基于CV法的EBK-LPI-RBF法 的RMSE为8.68 mm，MAE为5.52 mm，R2为0.99，满足精度要求。

表2 不同插值方法的精度比较

2.2.3 加权融合GPS-InSAR的地表形变实验与分析

加权融合的方法有很多，如变权组合法、算术平均法（等权组合法）、非线性组合法、最优加权法等，本文为实现降低GPS点偶然性和提高InSAR形变监测精度的统一，采用算术平均法融合GPS与InSAR数据获取矿区地表形变场。根据EBK-LPI-RBF法和InSAR形变结果图，统一像元分辨率（16.9×16.9）和坐标基准，得到加权融合的GPS-InSAR形变结果，如图7所示，可以看出，沉降区域分布在矿体西南部，与InSAR和GPS监测结果保持一致。

图7 GPS-InSAR累积形变结果

为验证本文提出方法的有效性，分别采用两种方案进行精度验证：①方案1，将GPS监测值作为真值，分析InSAR、EBK-LPI-RBF和GPS-InSAR三种方法在垂直向的形变精度；②方案2，将InSAR监测值作为真值，分析GPS、EBK-LPI-RBF和GPS-InSAR三种方法在垂直向的形变精度，结果如表3所示。

通过分析表3和诸图可以看出：

表3 不同研究方法的精度比较

1）方案1将GPS监测值作为真值，EBK-LPI-RBF法的内插结果精度达到毫米级，RMSE为8.68 mm，MAE为5.52 mm，R2高达0.99；但InSAR形变结果的精度较低，其RMSE为73.28 mm，R2仅为0.45。方案2将InSAR监测值作为真值，EBK-LPI-RBF法的内插结果精度略优于GPS监测结果。

2）对比两个方案发现，虽然GPS-InSAR法精度低于EBK-LPI-RBF法，但却很好地降低了离散点的偶然性，提升了InSAR的形变监测精度。两种方案中GPS-InSAR法的R2分别为0.88和0.77，很好地说明了该方法的可靠性。

3 结语

传统GPS与InSAR数据融合仅依据融合函数模型来提升InSAR形变精度，未考虑两种数据融合时本身性质不同所造成的影响。鉴于此，本文首先基于SBAS-InSAR技术获取研究区地表2019年4月—2020年6月的形变场，提取垂直向累积形变量；再针对GPS离散点形变反映面形变具有偶然性和InSAR监测不能客观反映实际形变的问题，综合了3种插值方法的优点，利用基于CV法的EBK-LPI-RBF法内插GPS点数据，并计算得到加权融合的GPS-InSAR地表形变结果。由精度评定指标可知，该方法有效降低了GPS点数据的偶然性，提高了InSAR监测精度，监测结果满足精度要求，为GPS与InSAR数据融合提供了新方法。