城市地表沉降作为一种常见、缓慢的地面沉降现象，一直是我国城市最严重、最普遍的地质灾害之一^［1］。昆明地处中国西南部的云贵高原地区，该地区本身存在着继承性下沉倾向，加之近几年城市化进展加快和地下水过度开采，截止到2004年就已形成面积超300 km²的地面沉降区域^［2］。姜朝松等^［3-4］通过早期水准数据分析显示，在1979~1998年间昆明市周围地区存在较为明显且有规律的地面升降运动，升降分界区域在滇池西岸附近，在滇池西侧呈现上升趋势，东侧则呈现下降趋势。东侧存在两个明显的沉降漏斗区域，其中一个位于官渡区的小板桥附近，自1985~1998年间累计沉降量达到227.5 mm，年均沉降速率达到30 mm/a。另一个沉降区的沉降中心位于渔户村西北方向的河尾村，自1993~1997年间沉降量达100 mm，年均沉降速率达到25 mm/a。

传统地面沉降监测主要依托于水准测量和全球定位系统技术，该方法可获得空间内特定点的沉降信息，但难以准确提取完整空间面域沉降特征。InSAR技术是进行地表沉降监测的基础方法，其最初利用干涉相位来测量传感器与目标物之间的距离^［5］。通过引入外部DEM数据或多轨InSAR结果进行差分（D-InSAR）来实现地面沉降的监测。然而，由于受到空间失相关、时间失相关和大气延迟的影响，D-InSAR技术监测地表沉降结果的精度受到明显影响。为了突破D-InSAR技术的这些限制，基于多时相InSAR（Multi-Temporal InSAR，MT-InSAR）技术的PS-InSAR技术、SBAS-InSAR技术和DS-InSAR技术被先后提出来^［6-8］。

PS-InSAR技术的优点在于选取的PS点具备稳定的高散射、高空间分辨率特性，能够有效克服时间、空间失相干和大气延迟效应等的影响，不足之处在于对影像数量的要求过分苛刻，一般需要25景影像以上才能保持较高的形变解算精度^［9-10］。SBAS-InSAR技术的优点在于通过小基线组合增加时间采样，能够获得比PS-InSAR技术更多、时间和空间基线更短的干涉对^［11］。由此我们发现，PS-In SAR和SBAS-InSAR两种技术存在优势互补。基于此，刘利敏等^［12］基于PS-InSAR 和SBAS-InSAR各自的优势，采用融合PS和SBAS的算法获取了北京地面沉降时空分布特征，测量的结果与水准结果相比，误差在±2 mm/a内。

虽然在SBAS基础上引入PS点能够有效改善形变解算质量，但是对于一些植被覆盖茂密或者城市建筑等高相干性点稀疏的区域，该方法仍然无法克服PS点稀疏使得这类算法应用受限的问题^［13-14］。Ferretti等^［15］提出了采用相干性较高的分布式目标（DS）来补充PS点不足的问题，亦可以称为DS-InSAR技术。该技术通过同质点选取算法增强InSAR时序协方差矩阵的估测精度，采用统计分布函数选取具有相同统计分布的同名像素，可以在保留空间分辨率的同时保证监测的精度。Hooper等^［16］研究表明，联合PS-InSAR和DS-InSAR或许可以提高沉降监测产品的精度；此外，DS-InSAR和SBAS-InSAR联合还可以克服SBAS算法中采用常规多视处理或空间滤波降噪方法带来的低分辨率影响，克服复杂沉降特征环境中非沉降特征的误读及细节信息的损失等问题。

近些年，MT-InSAR技术在诸多与昆明地面活动有关的形变监测和反演中体现出比常规测量手段更大的优势，也有部分学者^［17-20］将PS-InSAR、SBAS-InSAR等技术实践于昆明地表沉降变化研究。研究结果显示，自2014年以来，昆明城区地表整体存在形变，平均速率在±10mm/a左右，小板桥和河尾村依然是两个最大的沉降漏斗中心。然而，针对昆明区域海拔高、地形复杂、地表植被和水体覆盖造成的失相干、大气延迟效应和PS点选择困难的情况，采用现有时序InSAR数据及技术获得高精度的沉降信息具有极大挑战^［21］。

基于此，本研究拟获取欧空局（European Space Agency，ESA）2017~2020年Sentinel-1A VV极化数据，将Stacking思想引入SBAS中，通过对多个短时空干涉对进行叠加，在叠加结果中联合强度和相干阈值同时选取PS点和DS点，并对昆明市区地面沉降进行监测，针对典型沉降漏斗区域将该方法与已有的PS+SBAS-InSAR技术监测形变结果进行交叉检验，探索提出的方法在高原城市地表沉降监测中的有效性。

研究区位于昆明市区（24°50′44.16″~25°09′13.08″N，102°25′16.16″~102°49′38.75″E）， 横跨盘龙区、西山区、官渡区和晋宁区，具体研究区位于图1中红色矩形覆盖的区域。研究区在南北方位上整体地势呈北高南低的态势，并且由北向南逐渐降低。在东西方位上则表现为中部隆起，东西两侧低的态势。研究区处于云贵高原中部，平均海拔约1900 m。从地质构造上来看，研究区位于滇中高原地区，处于云南高原规模最大的第四纪继承性断陷盆地，以南北向构造为主，周围断裂带分布密集，包括小江断裂带、黑龙潭—官渡断裂带和南北向的普渡河断裂带。研究区以湖泊地貌为主，土壤类型以湖沼相粉砂土和软粘土为主，易受到流水侵蚀作用。近年来，随着城市化的快速推进，高层建筑迅速增加，工业用水猛增导致承压水位迅速降低，这些因素进而引起研究区内地面沉降现象频发。

本研究共获取了36景C波段Sentinel-1A卫星数据。时序范围覆盖为2017年3月15日~2020年2月16日，共计1 068 d，约3 a时间。C波段Sentinel-1A卫星数据由欧洲航天局（European Space Agency）发射的双星地球观测卫星星座提供，数据具有轨道精度高、重访周期短、影像分辨率高、数据获取便捷等特点，目前被广泛用于基于时序InSAR技术地表形变、滑坡、地震等灾害探测与分析中^［22-28］。C波段Sentinel-1A卫星数据的重访周期为12 d，像元方位向和距离向的空间分辨率分别为2.3 m和13.9 m。为了降低SAR成像参数变化对监测结果的影响，本研究中获取数据的极化方式均为VV极化，入射角均为39.3°。本研究中所有数据获取的具体时间、时间基线和空间基线等详细信息如表1所示。

在干涉处理中，为使时间失相干对影像配准的影像影响最小，选取时序上居中的2018年10月12日获取的影像作为主影像^［13］，其余35景影像作为辅影像。然后基于强度互相关和增强谱分集算法^［10］（Enhanced Spectral Diversity，ESD）来消除Burst间出现的相位跳变，最终实现精确配准，本研究中的配准精度达到了亚像元级，满足本研究中干涉测量的要求。尽管在选择主影像时考虑了尽量缩短时间基线，但是由于C波段数据穿透性弱，在植被、水体等地表受到时间失相干的影响仍然较大，影像中这些区域的整体相干性较低，为了获得较好的干涉相干性，保证足够的干涉对组合及较好的干涉效果，在经过多次实验尝试的基础上，将空间基线的阈值设置为300 m，时间基线的阈值设置为60天。文中最终确定了46个小基线干涉对，干涉对的时空基线信息如图2所示。

融合PS、SBAS、DS的InSAR方法的基本思想是将三种方法的优势互补，较好地解决了采用InSAR技术进行地面沉降监测过程中干涉时空失相干、大气延迟、SAR数据不足、PS点选取困难等问题。算法流程包括：首先采用SBAS方法形成干涉对，然后将Stacking思想引入SBAS算法中，通过对多个干涉对进行叠加，在叠加结果中以一定的干涉阈值同时选取PS点和DS点；以选取的PS点和DS点为基础，借助DEM数据进行差分干涉处理并建立Delaunay三角网，采用最小二乘网络求得各点的沉降速率、DEM误差相位和残余相位；对各残余相位进行SVD分解并滤波进而获得各点的时序形变量。算法的基本原理如下：

设获取了覆盖同一研究区的N+1幅SAR影像，为保证各个干涉对之间有较好的相干性，通过设定时间和空间基线阈值对所有影像进行干涉对组合并进行差分干涉，可得到M幅差分干涉图，其中M满足式（1）。

利用外部DEM数据，根据成像几何关系，去除干涉图中平地相位和地形相位。研究中PS、DS点的提取是基于全局高质量的差分干涉相位来提取，假设全局差分干涉相位均被一个稳定的相位点精确校正和解缠，时序内的任意时刻tii=1,⋯,N相对于参考时刻的差分干涉相位φti和数据处理获取的差分干涉相位δφtk可分别表示为：

得到一组时序差分干涉图之后，基于解缠相干性阈值提取时序内保持稳定相位的PS和DS点。

假定t0时刻为初始形变为0，在干涉基线不为0的情况下，则在第i幅第j点的差分干涉相位Δφjdef可以表示为：

其中：∆φj为第i幅第j点的干涉相位；λ为雷达波长；dtA、dtB分别表示tA和tB相对于t0时刻的累积形变量；∆φjatm为大气扰动引起的延迟相位；Δφjtopo为地形起伏引起的相位；Δφjflat为平地相位；Δφjnoise为噪声相位；φtA为tA时刻的相位值；φtB为tB时刻的相位值。

在经过差分干涉处理之后，假设在tA~tB时刻地表形变速率满足线性变化，则tA~tB之间的累积形变量可表示为

据式（4）可知：对于研究区N个时刻M幅干涉对的形变量可以线性表示为：

其中：系数矩阵A［M×N］的每一行对应于1幅干涉图，每一列对应于1景影像每个元素代表一个干涉相位，∆φ为M幅差分干涉图中高相干点解缠相位组成的矩阵，φ为高相干点在N个时刻的未知形变相位组成的矩阵。

由于组成的干涉图集均限定在唯一的短基线算法内，因此可以利用最小二乘法求得形变相位，即：

进一步将相位信息的求解转化为对相位变化速度的求解，可以得到如下矩阵方程：

其中：B为M（N-1）矩阵，通过式（7）就可以得到速度矢量v的最小范数解。算法实现的具体流程图如图3。

相较于传统PS-InSAR的单一主影像法选取干涉对而言，融合PS、SBAS、DS的方法引入小基线集思想，即通过设定严格的时间和空间基线阈值将各个SAR影像进行自由组合，不仅增加了干涉对的数目，而且还确保了干涉对中PS点的高的相干性。由于Sentinel-1A卫星提供精密的轨道数据，且空间基线一般不超过500 m，因此影响干涉相干性的主要因素为时间失相干。基于此，本研究选用控制空间基线基本一致的情况下，将PS-InSAR选择的干涉相干性与本文方法选择的干涉相干性进行定量分析。两组干涉对时、空基线信息如表2所示。

图4采用直方图对比了两种方法获得的相干性信息。采用PS-InSAR方法获得的整体相干性分别为0.227 1、0.228 4，低于本文方法所获得的0.255 5、0.254 7。在相干性0~0.4间，采用PS-InSAR技术获得的相干像元更多，而在相干性0.4~1区间内，本文提出的方法获得的相干性明显优于PS-InSAR技术。这是因为本文方法融合了小基线集思想，干涉影像的时间基线更短，最长时间间隔只有60 d，而PS-InSAR采用单主影像干涉方式，无法优选干涉影像间的时间和空间基线，在对于监测长周期地表形变背景下，容易造成时间失相干现象。

常规D-InSAR处理中， 理想状态下，在完全消除地形、大气延迟等相位后，差分干涉相位只包含地表形变信息。但是由于没有实时的大气数据进行改正，因此大气延迟相位会严重影响形变信息。本文方法继承了PS+SBAS-InSAR方法，将地表形变相位、DEM误差相位作为未知参数，利用多幅PS干涉图联合求解并将其从差分干涉图中分离出来，然后根据时序差分干涉图中组成残余相位成分的不同时空分布特性，利用高通和低通滤波的方法逐步进行分离，从而将大气延迟相位进行分离^［10］。

针对本研究所用的36景Sentinel-1A数据而言，PS-InSAR技术采用时间序列内一景影像

为主影像，将剩余35影像与主影像进行干涉，最多得到35个干涉对；本研究所用的方法融合了SBAS技术，若不设置严格的时间和空间基线阈值，理论上能够生成630个干涉对，而本次在设置严格的时空基线阈值后，仍然得到46个干涉对，因此在SAR数据量不足的情况下，仍然能够获得更多的干涉对。

融合DS点的PS+SBAS-InSAR技术对部分PS点进行了相位滤波处理，一定程度上降低了PS点的空间分辨率，为了定性和定量分析本文方法反演昆明地区形变速率的有效性，本文提取PS点的sigma（相位标准差）、有效PS点数量和反演的高程误差3个量化因子，综合对比本文方法与PS+SBAS-InSAR方法提取的PS点质量。两种方法求取的高相干目标点的sigma值见图5。

由图5可知，两种时序InSAR方法得到的sigma值均在0~1.3内。在PS+SBAS-InSAR法中，PS点经过3次迭代后，sigma相位标准差和sigma平均值均得到优化，其中simga值平均降低17.03%，sigma值标准差降低14.88%。在本文所提出的融合PS、SBAS、DS-InSAR技术中，经过3次迭代后，sigma值显著降低，其中simga值平均降低31.51%，sigma值标准差降低28.05%，PS点质量优化较PS+SBAS-InSAR方法更加明显。表3记录了PS点的详细质量信息，由表3可知，本研究提出的方法获得的PS点质量较优，这是由于PS+SBAS-InSAR方法采用振幅离差指数和光谱相干性进行提取PS点，同时考虑了时间域像素强度稳定性和点目标的光谱特性信息，因此提取的PS点质量略高，但从PS点数量分析，本研究所提出的方法在PS+SBAS基础之上融合了部分DS点，使得提取的PS点数量相比于PS-InSAR方法大大增加，增加数量达648 075，因此在反演形变速率更有优势。

此外，研究还对两种方法提出的高相干点高程值均进行了改正。由图6直方图统计结果可知，采用PS+SBAS-InSAR方法的所有高程改正量的标准差为4.39，96.79%的点的高程误差改正在 ±10 m以内，而采用PS+SBAS+DS-InSAR方法的所有高程改正量的标准差为4.91，92.59%的点的高程误差改正在±10 m以内，本研究提出的方法结果略优。

通过采用上述两种时序InSAR技术提取昆明市的地表位移时间序列，并将反演的LOS（Line Of Sight，LOS）方向形变速率采用UVA(D)≈DLOScosθ公式转换为垂直方向上形变速率，我们得到了研究区域的沉降灾害结果。图7显示了该结果，图中红色区域代表该地区发生了地面沉降。

由图7可知，两种选点方法反演的结果在空间尺度上有高度一致性，形变量也大致吻合。其中PS+SBAS方法监测的最大沉降速率达到 -99.6 mm/a，本研究方法监测的最大沉降速率达 -81.5 mm/a，PS点的选取方法不同是造成沉降速率存在略微差异的主要原因。从整体分布来看，昆明市地表沉降呈多中心聚集态势，滇池东岸存在5个明显的沉降漏斗区域，且有连成一片的趋势（图7中A、B、C、D、E）。

A区以小板桥至广卫村为沉降漏斗中心，沉降面呈南北走向，形成椭圆形的槽形沉降面，周围包括子君村、晓东村、小喜村、雨龙村和东菊新村等地都存在一定程度的沉降现象，该区域沉降速率多集中在-15~0 mm/a，最大下沉速率达-59.9 mm/a，3 a最大累积沉降量达-175.5 mm。B区位于昆明官渡区的西南方位，该区域沉降分布较为分散，沉降漏斗分布最广且不均匀，沉降区域覆盖多个小型沉降中心，沉降面呈西南-东北走向，包含福保村、陶李村、曹家桥村、陈家营、严家村、叶家村、新螺蛳湾、新二村和海埂村等多个沉降明显区域，沉降速率多集中在-17~0 mm/a，最大下沉速率达-78 mm/a，3 a最大累积沉降量达-220.4 mm。C区包含2个沉降漏斗，分别是以河尾村—渔户村、小渔村的两个沉降中心，涉及新河村、河尾村、大坝村、福海村、渔户村、邬大村、万达广场、小渔村等，沉降趋势呈东北部发展态势， 沉降速率多集中在-12~0 mm/a，最大下沉速率达-86 mm/a，3 a累积最大沉降量达 -242 mm。此外，本次研究还发现了棕树营—下栗树村存在小型聚集型沉降区（D区），且与大塘子沉降中心相连。值得一提的是，本次在滇池西岸及北岸也探测出新型小型沉降带（E区），沉降中心位于碧鸡村附近，并与小渔村沉降中心相连，该地区沿滨河流域沉降速率明显高于其他区域。此外，蒋家营也首次探测出具有明显的沉降现象。

融合PS、SBAS和DS的方法不仅继承了多基线下干涉影像对的强相干特性，还能提取部分相干目标点，极大地增强了PS点选取的密度。本次研究发现，采用经典的PS+SBAS技术共选取出749 327个PS点，融合PS、SBAS和DS的技术共提取出1 191 616个PS点。我们采用图8~图9首先定性化地对比分析了本文方法与经典PS+SBAS方法选取的PS点效果及反演的沉降监测结果。由图8~图9可知，4个PS+DS点位于小板桥至广卫村沉降漏斗中心，另外4个点分布于河尾村沉降漏斗中心。可以看出，融合PS、SBAS和DS的方法与经典PS+SBAS方法反演得到的形变速率存在很好的一致性，沉降速率以线性形变为主。选取出8个典型沉降区的目标点采用两种方法得到的形变速率相关性分别为0.963 6、0.980 2、0.954 7、0.964 5、0.984 3、0.991 2、0.979 6和0.981 1，两者技术形变速率平均误差为6.34 mm/a，融合PS、SBAS和DS结果表现出高精度探测地表形变的潜力。从两个典型沉降区域的年沉降变化和PS点的分布来看（图8、图10），两种结果一致性较好，但也存在微小差异。在2017年3月~2018年2月期间，针对经典PS+SBAS技术采用光谱相干值和振幅偏离差两种选点方法反演出两区域沉降速率基本维持在-20~-10 mm/a，而在2017年3月~2020年2月，小板桥沉降速率增加到-40~-20 mm/a，河尾村增加到-40~-30 mm/a，但融合PS、SBAS和DS的方法计算出的两区域沉降速率增加不明显，大部分还是-30~-20 mm/a，这是因为本文提出的选取方法是全部解缠干涉影像中均保持相位稳定的点，而PS+SBAS选点考虑了影像的散射特性，PS点分布基本都是具有形变特征的城镇区域，所以反演沉降速率略大，但对于沉降变化量大的非城市区域，本研究方法在获取PS点基础上，融合了部分分布式点，因此选取的PS点更密集，即使在非城镇等区域反演的形变结果也更明显。

从图8可以看出，河尾村3 a累积沉降量达到242 mm，而小板桥沉降区3 a累积沉降量为182.8 mm。此外根据图9显示，从2017年3月至9月期间，小板桥和河尾村沉降速率变化最明显，而在2017年10月~2019年9月期间，沉降速率变化较慢，监测曲线存在轻微的阶梯状变化特征。

由于研究区缺少同期水准数据，无法进行多种数据源监测结果的交叉验证。因此，本研究通过获取历史沉降资料^［3-8］，选出位于研究区的8个村落5个时间段内的沉降速率与本研究所用的融合PS、SBAS与DS-InSAR方法监测的沉降速率进行对比，验证其技术识别形变区的准确性。不同时间段内沉降速率对比结果如表4所示。

表4显示了本文提出方法得到的昆明城区地面沉降速率与历史沉降监测数据结果的对比。由表4可知，自1994年以来，各个区域的沉降速率变化不大，基本都在维持在-30~-20 mm/a。但小板桥自2014年起，沉降速率明显加快，据历史资料显示，该区域历史地面沉降诱导因素主要是地下水的过度开采和粘性土层透水性弱，其中地下水抽汲权重高达72%^［8］。小板桥沉降漏斗在2014年开始，新螺蛳湾片区作为子君村的城改项目建设，使得常住人口比之前增加12 631人，建筑用地的不断扩大，加之地下水开采显著增加，当水位降低后，越来越多的土层参与固结压缩，极大的增加了地面荷载，使得该区域沉降速率显著增加^［2］。此外，位于小板桥沉降漏斗中心北部的东菊新村也探测出了轻微的沉降现象，因此对该片区亟需采取相应的措施来缓解大面积的地表形变。结合图8中的各村沉降速率分布图来看，位于古滇池沉积区内的河尾村沉降漏斗区已与小渔村区域连成一片，沉降严重区域主要发生在靠滇池水域边缘部分，这是因为昆明湖滨土质以软土、细粒土为主，具有高含水量、高压缩、低承载力、强固结变形等特点，加上近几年建筑用地开发量大，因此这些地区成了形变较为严重的区域，这与之前学者研究的结论相一致^［5］。另外位于D区的棕树营区在2015~2017年期间首次探测出沉降现象^［5］，而在2017~2020期间也逐步形成小型沉降漏斗，年沉降速率在-23~0 mm/a，且有与小渔村沉降漏斗连成一片的趋势（图11）。

地质构造运动也与昆明地区沉降分布有着较高的相关性，本次形变最为严重的区域分别为小板桥—广卫村中心和河尾村-福海村为沉降中心的区域，其中小板桥沉降中心位于黑龙潭地块南端，南市区位于小江断裂带与普渡河—西山断裂带的南中部^［4］。构造运动的继承性、多期性、间歇性造成了当地地貌形态具备成层分布的特征^［5］。

本研究使用C波段Sentinel-1A雷达数据，融合PS、SBAS、DS 的InSAR技术对昆明市2017~2020年沉降速率分布进行探测，并结合历史水准数据和地质资料对形变结果进行了验证，并综合分析地面沉降变化原因。研究结果表明，昆明整体存在-22~8 mm/a的形变速率，地面沉降区域较以往有所扩张，其中小板桥、河尾村、小渔村、新二村等沉降中心已经连为一体，最大年沉降速率达到81.5 mm/a。昆明沉降分布具有明显的空间分布特性，滇池东南方向主要以沉降为主，西北向则表现为抬升。普吉和棕树营区域较以往研究具有明显沉降变化，东北方向的蒋家营则是本次研究中发现的新沉降点，说明了沉降区域有向东发展的趋势，后期应当加强该区域的治理和维护。此外，通过对比小板桥和河尾村两个历史典型沉降漏斗上选取的8个形变特征点，验证了PS+SBAS+DS-InSAR技术作用于昆明这类复杂地形城市的地表形变监测的有效性。本实验提出的方法可提高非相干目标点的选取数量，且获得了比常规PS+SBAS方法更直观的形变反演结果，但同时也因为相位滤波而牺牲空间分辨率，降低了部分选点质量。因此，未来应加强对于SAR影像的相位质量恢复，并辅以同期地面水准数据进行验证，从而达到更高精度的形变监测。