地面沉降给人类生产生活带来巨大危害^[1]，降雨和海平面上升等是影响地面沉降的自然因素，地下工程施工、地面荷载增加以及地下水抽取等人为因素引发的地质灾害也会影响区域经济和社会可持续发展^[2-4]。地表形变监测的主要手段为水准测量和GNSS监测，水准测量受人力、物力和财力限制，布点少、路线稀疏、监测周期长、时空分辨率较低，难以对地表形变进行高精度大范围监测；而GNSS技术虽然可以获取连续的地表形变监测结果，但受限于稀疏的监测点密度和昂贵的地面设备。合成孔径雷达干涉测量技术(InSAR)是一种先进的影像大地测量技术，拥有全天时、全天候、大范围、高精度等诸多优点，在地面沉降监测中得到广泛应用^[5-6]。多时域永久散射体雷达干涉测量技术(PS-InSAR)是基于传统InSAR技术和D-InSAR技术发展出的一种新时序形变分析技术，通过对多景SAR影像进行联合分析，能够较好地抑制干涉图中时空失相干和大气噪声影响^[7-8]，非常适合城市大面积地表形变监测和灾害普查。

截至2019年，华北平原部分地区的沉降速率大于100 mm/a，其中90%的沉降区位于北京、天津、河北地区^[9]。为探究天津区域新近沉降状况，本文以天津市及其周围区域为研究对象，首先利用2019-01~2022-03共90景Sentinel-1A卫星影像和PS-InSAR方法获取研究区的地表形变空间分布与显著沉降区的形变时间变化趋势；然后分析地表形变的地理分布特征，讨论时序形变与降水量、地下水等因素的相关性。本研究对未来天津市地表沉降的监测治理和水资源的合理规划利用具有重要意义。

1 研究区概况和实验数据

研究区域主要包括天津市及其周边地区(38°34′~40°15′N、116°43′~ 118°04′E)，具体见图 1。实验所用数据为90景Sentinel-1A卫星影像升轨数据，空间分辨率为5 m×20 m，时间间隔为12 d，时间跨度为3 a(2019-01~2022-03)，入射角为39°。为提高信噪比和处理效率，对原始SAR影像进行20 ∶5多视处理，通过裁剪原始影像获取环绕天津市区及其周边地区的SAR影像。图 2为PS-InSAR的时间基线和空间基线图。

2 方法 2.1 基于分块误差校正的PS-InSAR技术

根据SAR干涉测量原理，PS-InSAR干涉相位φ可以分解成若干分量：

$ \varphi=\varphi_{\mathrm{flat}}+\varphi_{\mathrm{topo}}+\varphi_{\mathrm{def}}+\varphi_{\mathrm{atm}}+\varphi_{\text {noise }} $

(1)

式中，φ_flat为参考椭球(又称平地)相位，φ_topo为地形相位，φ_def为沉降引发的相位分量，φ_atm为大气相位。在PS弧段建立包含形变量和高程误差的相位模型，并利用整数最小二乘法估计弧段参数。

在大气影响严重的沿海区域，时空滤波和多项式拟合是常用的误差去除方法，但这2种方法难以去除复杂分布的误差。经典的六次多项式公式如下：

$ z=a_1+a_2 x+a_3 y+a_4 x y+a_5 x^2+a_6 y^2 $

(2)

式中，a_i为多项式待估参数，x、y为像素坐标，z为像素值。本文采用分块误差去除方法，将SAR影像分为相等大小的矩形块，小矩形块内的大气和趋势误差可以利用三次多项式拟合去除：

$ z=a_1+a_2 x+a_3 y $

(3)

具体实现过程如图 3所示。主要步骤如下：

1) 通过SAR影像配准和干涉，在差分相位图基础上，通过建立相位和参数之间的函数模型，利用整数最小二乘估计器来解缠相位并通过空间平差得到PS点的形变速率和高程误差。

2) 对含有误差的速率进行分块去误差处理，矩形块大小设置为200×200，块与块之间重叠区域设置为50。在每个影像块内，利用相干点计算3个参数。

3) 在移除线性形变基础上，对残余相位进行时空滤波，去除大气误差。最后结合线性形变和非线性形变得到最终形变时间序列结果。

2.2 随机森林分类

随机森林分类器可实现Breiman^[10]提出的随机森林分类算法，即使用Bootstrap取样和决策树进行预测。分类实验具体实现步骤如下：

1) 选择强度影像进行样本标注。本文选取2021-03-07的SAR影像，可以有效划分水域边界。人为经验性标注样本，将SAR影像标注为建筑区、水域和农田、草地、裸地区3类。

2) 利用SAR影像强度特征进行模型训练，随机选择70%的样本作为训练集，剩余30%的标注样本作为验证集，用来验证模型精度。

3) 利用训练好的模型对整个研究区进行预测和分类。

3 实验结果与分析 3.1 PS-InSAR形变结果

PS-InSAR形变监测结果如图 4所示，其中蓝色和紫色为沉降区域，绿色为稳定区域，黄色和红色为抬升区域。由图可见，基于分块误差的校正方法可有效去除局部大气和整体趋势，获得更加准确的形变估计结果。由图 4(c)可见，2019-01~2022-03期间研究区发生明显地面沉降，主要分布在天津市外围，位于河北省霸州市胜芳镇。其中a、b两个区域沉降最为严重，最大沉降速率出现在区域a中，沉降速率约为80 mm/a。为进一步验证PS-InSAR监测结果的准确性，获取研究区内相同监测时段的10个二等水准点数据，水准点分布如图 4(c)所示。将水准测量数据与InSAR监测结果进行对比，结果如图 4(d)所示。由图可见，原始形变速率、基于六次多项式误差校正的形变速率以及基于分块误差校正的形变速率与水准数据的均值偏差分别为10.3 mm、31.4 mm和8.1 mm，与水准数据的标准差分别为6.7 mm、9.5 mm和6.8 mm。由此可知，基于分块误差校正结果具有更小的均值偏差和标准差，可用于PS-InSAR结果校正。

3.2 成因分析

由图 4可见，2019-01~2022-03天津区域的最大沉降漏斗发生在河北省廊坊市霸州市胜芳镇。图 5为天津区域土地利用状况和地表沉降的空间分布关系。其中，图 5(a)为SAR强度影像，图 5(b)为利用随机森林(RF)算法和强度影像得到的天津区域地表分类图，包括建筑区、水域和农田、草地、裸地区3类。实验分类精度达到89.1%，说明随机森林分类效果良好。由图 5(c)、(d)可见，天津区域沉降主要发生在建筑区，这些区域的地表沉降与地下水的过度开采以及地面荷载增加有关。此外，农业种植区也发生一定量沉降，这些区域的沉降与农业灌溉有关。由图 5(e)、(f)可见，建筑区沉降主要发生在胜芳镇和左各庄镇中心。

为进一步分析研究区的沉降变化规律和形变成因，以2019-01作为起始时间，假设此时的形变量为0，选取目标区域内a、b和c三个特征点进行形变时序分析(图 6)。由图 6(a)可见，沉降最大点位于胜芳镇中心，沉降点a在3 a内沉降约250 mm，且降水带来的地下水增加并未减缓该建筑区的沉降趋势；由图 6(b)可见，沉降点b位于农业种植区，该监测点的形变量与降雨量密切相关，呈季节性变化趋势，降水量增加会导致深层地下水量增加，在一定程度上会抵消地下水抽取带来的影响，从而减缓地表沉降趋势，使得该监测点在3 a内仅沉降140 mm；由图 6(c)可见，天津市无明显的线性沉降趋势，降水量引起地表发生周期性变化，变化幅度小于10 mm。图 6(c)为大气趋势误差校正前后的时序结果对比，由图可见，相较于误差校正前的明显抬升趋势，误差校正后的周期性形变趋势更加符合实际情况，再次验证了基于分块误差校正的PS-InSAR方法的有效性。

4 结语

本文利用PS-InSAR技术对天津市地表形变情况进行监测和灾害普查。结果表明，该地区形变的主要原因包括地下水过度抽取、地面荷载增加、降雨量增加等。分析天津区域土地利用状况和地表沉降速率发现，研究区形变主要发生在建筑区，尤其是胜芳镇和左各庄镇。此外，部分种植区也出现明显的地表沉降，这些区域的农业灌溉需注意节制抽取地下水。本文可为天津市及其周边区域的沉降治理和灾害预警提供参考。