基于PS-InSAR技术的太原地面沉降监测研究

李 路,洪友堂

(中国地质大学(北京)，土地科学技术学院，北京 100083)

地表形变是一种由人为或自然原因造成地表变形现象，世界上各个国家都受到不同程度的地表形变影响。据统计，自1949～2010年，我国由于地面沉降造成的经济损失累计高达5 000亿元，同时也造成了重大人员伤亡事件发生。目前，传统的地表形变测量仪器主要有水准仪、经纬仪、全站仪、GPS等，它们虽然精度高但需要人员实地监测，且耗费大量人力物力。随着科技发展，合成孔径雷达差分干涉测量(DInSAR)技术凭借全天时不接触、不易受环境影响、监测效率高等优势，极其适用于对大面积地面进行微小形变研究，精度可达毫米级，相比于传统测量技术，具有极大优势。但是由于地球大气层中的大气噪声、形变前后地面物体发生位置变化导致时空失相干等众多因素存在，限制了该技术的发展。2002年，Ferretti等人首次提出PS-InSAR技术并对研究区进行滑坡监测，处理结果与实测数据基本吻合。随后，Berardino等人提出小基线集差分干涉测量(SBAS-InSAR)技术，并证明有效性。经过十余年的发展，PS-InSAR技术与SBAS-InSAR技术得到了显著发展与应用。而PS-DInSAR技术主要依靠散射特性稳定点目标，而城镇区域中的房角、电杆等稳定地物较多，因此该技术极其适用于对城镇区域进行地表沉降监测研究。

2015年，杨海瑞[1]等人研究了太原盆地地下水流和土体变形耦合的机理,并对未来20年太原盆地地下水水位及地面沉降情况进行预测。2016年孙晓涵[2]等人全面分析了太原盆地地裂缝和地下水开采、地面沉降在时间、空间、活动上的相互关系。2017年赵强[3]等人结合太原盆地地面沉降的发育情况,从地下水开采、粘性土分布不均、断层构造、地面载荷等方面,分析了影响地面沉降的因素。2018年周艳萍[4]等人根据五个沉降中心中30个典型的水准观测点的累积沉降量建立了灰色Verhulst预测模型，最后预测了2010年与2015年的地面沉降发展趋势。本文利用分辨率为3 m的COSMO-SkyMed数据，提出基于PS-InSAR技术对2013年6月到2016年6月期间太原区域的地面沉降进行分析，结果显示太原地区共形成以万柏林、南堰、小店地区为中心的3个沉降漏斗，其中沉降最严重的是小店地区，最大累计沉降量达-152.7 mm。

1 PS-DInSAR原理及关键步骤

1.1 PS-DInSAR方法

该方法是意大利学者Ferretti等人于2002年提出，并对研究区进行滑坡监测，处理结果与实测数据基本吻合，证明了方法的有效性。该方法主要以获得的在大时空基线情况下仍保持较高相干性(如房角、电杆等)的永久散射体(PS点)为基础，这些PS点受时空影响较小，更易于分离大气、地形与轨道误差相位，获得精确的地表形变相位。该方法可以有效地消除或削弱大气效应、时空失相关的影响，超越传统时间、空间基线距上的临界基线条件限制，提高了数据利用率。经过10多年的研究，利用该技术已经成功在昆明、北京等地进行地面沉降监测研究，此外一些学者在主影像选择、PS点选取、相位解缠等方面进行了改进，获得了更为精确的地表结果信息。

PS-DInSAR方法的核心思想为：从一系列N景SAR影像中选取其中一景作为公共主影像，其他影像作为辅影像，将所有辅影像与主影像进行配准、重采样、干涉处理与去除地形相位，获得N-1景差分干涉图，选择PS点，建立函数模型，并利用回归分析估计线性形变速率和高程残差，利用时空滤波分离大气相位成分、非线性形变成分与噪声相位成分，最终获得精确的地表形变相位信息[5]。具体流程图如图1所示。

1.2 方法处理关键步骤

1.2.1 PS点组合选取

在PS-DInSAR方法处理流程中，PS点提取是处理的第一步，PS点选择的质量和数量决定数据处理结果的优劣性，这些点分布于裸露岩石、高大建筑物、大金属支架等地物上，具有在大时空基线情况下仍保持较高相干性的特点。常用的点提取方法有相干系数法和振幅离差法。其中相干系数法为主要方法，依据公式(1)，计算目标像元(a，b)在整个时序的相干系数，然后与给定阈值做比较，以确定其是否可选为PS点。通过幅度离差法可以滤除海面上的点。

(1)

式中，M、N为窗口的高度和宽度，φ(a,b)为地形相位项。

1.2.2 相位解缠优化

相位解缠是 InSAR 数据处理流程中重要步骤之一，相位解缠的准确度在很大程度上决定了获取地表形变信息的准确性，算法有路径跟踪算法、最小二乘算法、最小费用流算法。本文为了更准确计算残余相位，首先通过积分求得平均DEM残余和平均线性形变相位，再从差分干涉相位中将其去除，得到残余相位。为了提高残余相位解缠的稳健性，对其进行空间域滤波处理后再解缠。考虑到相位解缠中实际情况以提高解缠能力，引入权值思想，求得最小最优解缠模糊数，见公式(2)。式中，K(a,b)为影像中(a,b)处的相位值，wm为权值，本文将像元振幅值作为权值，在M个时域内进行相位解缠。

(2)

2 研究区概况与数据处理

2.1 研究区概况

研究区域为太原市城区及其部分周边地区，主要位于山西省中北部，太原盆地北部。地理坐标范围为东经112°19′～112°48′，北纬37°40′～38°1′，覆盖面积为36×32 km2。以研究区地形图为底图，叠加影像的强度均值图，如图2所示。据太原市114个二等水准测量点实测资料，1956～2003年地面沉降范围南北长约39 km，东西宽约15 km，至2003年已形成西张、万柏林、下元和吴家堡4个沉降中心[6-7]。最大沉降中心为吴家堡地区， 年均沉降速率为63.0 mm/a。近年来，由于超采地下水引发的地面沉降也愈演愈烈，沉降范围向盆地边缘扩展，沉降漏斗面积扩大，对太原地区进行地表变形监测刻不容缓[8]。

2.2 数据处理

图2 研究区域

本文使用瑞士的GAMMA雷达数据处理软件对37景高分辨率COSMO-SkyMed SAR数据进行实验，时间跨度从2013年6月～2016年6月，X波段(3.1 cm)，地面分辨率为3 m×3 m，极化方式为HH，入射角为24.94°～28.36°，右视，升轨,幅宽为40 km×47 km。采用美国NASA提供的分辨率为30 m的SRTM1 DEM数据。

表1 研究区域数据详细情况

选择时空基线居中的2014年8月9日影像作为单一主影像，对配准好的辅影像进行干涉处理，形成36对干涉对，如图3所示。平均垂直基线385 d，平均时间基线179 d，可以看出基线分布较稳定，相干性较高，DEM误差对形变相位影响较小，结果更精确。

本文采用相干系数法和振幅离差法组合选点方法，相干系数阈值设置为1.5，PS点数为33 735个，振幅离差阈值设置为0.3，PS点数为3 772 458个，共选择初始PS点3 787 263个。本次研究区域为太原城区，因此PS点密度较高，平均每平方公里PS点达两千余个，且广泛分布于研究区域中的房顶、桥梁、火车站、机场等地，更有利于进行区域地表变形分析。基于选取的初始候选点，提取其高程、基线、相位等信息，并利用DEM高程数据去除地形相位，形成基于干涉点的差分干涉相位。对各点的差分干涉相位先进行空间滤波降噪处理，采用的是50×50的均值滤波窗口，之后，对获得的稳健的点差分干涉相位进行多路径回归分析处理，来估计高程改正、线性形变、残余相位等。对于残余相位再次进行缠绕、滤波、解缠处理，以精确估计大气相位。将获得的高程改正、大气相位代入原式进行修正，以进行迭代处理。通过迭代，PS点数量逐步降低，残差标准差降低，点质量得到提升，模型也越来越稳健，以达到优化目的。采用基于加权最小费用流相位解缠算法，并使用最小二乘算法来校正基线信息，以获得可靠的结果。另外，利用SVD算法对每个干涉点目标的时间序列相位进行处理，求算2013年6月～2016年6月期间的年形变速率与累计形变量。

图3 时空基线分布

3 沉降结果分析与验证

3.1 沉降结果

图4 研究区域年沉降速率图

以光学遥感影像作为底图，数据处理获得太原研究区域的年沉降速率图，如图4所示，将区域A、B的年沉降速率情况放大至右侧。A区域为万柏林地区；B地区为南堰地区；C地区为小店地区；D地区为吴家堡地区；E地区为西张地区。可以很明显得看出，万柏林、南堰、小店地区有明显沉降出现，而吴家堡、西张地区沉降平稳。在太原南部的小店地区沉降现象最为严重，最大值累计沉降量为-152.7 mm，平均年沉降速率达50.901 mm/a。

3.2 地表变形分析

查阅太原的沉降历史资料发现，在2000年前市区主要有四个沉降中心，分别为西张、万柏林、下元和吴家堡，而小店地区沉降不明显，在2003年后，这四个地区的沉降现象趋于缓和，局部区域出现反弹，南部区域沉降现象加剧。因此本文主要基于此四个沉降中心与小店地区进行地表变形分析。

A区域-万柏林地区，地质构造为太原西山向斜煤盆地，汾河流域陷落区，属于旧有沉降区，1989年至2000年该地区年沉降速率为46.73 mm/a，随后数年缓和很多，在2010～2013年仅为-14.6 mm/a。根据《太原市万柏林和平老工业区搬迁改造实施方案 (2013年～2020年)》，重点推进下元、南寒两个村的改造，从数据处理结果看，在2013年6月至2016年6月期间，该地区最大累计沉降量达-93.61 mm，年平均沉降速率达-31.20 mm/a，符合区域实际情况。该地区在2016年初有回升趋势，可以初步估计该区沉降速率有所缓和，查阅资料发现该地区基层社区较多，共有7个社区居委会和2个村委会，该时间段内新建高层建筑诸多，正常沉降，暂时没有安全隐患。万柏林地区2013年6月～2016年6月的时间序列沉降量图(如图5所示)。对万柏林地区的37°52′5″纬度线上、东西向约2 km，做时序分析，沉降漏斗剖面图表现年沉降速率变化情况，如图6所示，从图中可以看出该地区很明显为沉降漏斗，最低点的年沉降速率最大可达31 mm/a。

图5 A区域地表监测累计沉降量

图6 A区域沉降漏斗剖面

B地区-南堰地区，与吴家堡地区相近，同属旧有沉降区，在2000年以前为太原市区最大的沉降区，本文数据处理结果发现研究期间该地区累计沉降量达-105.78 mm，平均沉降速率为35.26 mm/a。从遥感影像上发现该地区存在大面积的建筑用地，说明正在进行施工建设或拆迁工程，从提取的PS点数量较少可以说明此地时空失相干较为严重，也可以证明这点。

C地区-小店地区，位于太原市区的东南部，晋中盆地的北端，有明显沉降漏斗出现。在2000年前该地区为非沉降中心，在2003～2010年的年沉降速率为27.9 mm/a，在2010～2013年的年沉降速率为36.6 mm/a，随着太原地区经济南移，使得小店地区经济发展迅猛，小店地区地面沉降现象日趋严重、沉降速率日趋增加、沉降中心逐渐形成。从结果看，在2013年6月～2016年6月期间，小店高新区的沉降中心范围逐渐扩大，沉降量逐渐变大，最大沉降量达152.7 mm，年沉降速率达50.901 mm/a。整个区域沉降趋势较大且有加剧迹象，整个沉降漏斗面积达40 km2。在小店区域的中部、东部部分地区有较严重沉降现象，而东南部较稳定。对纬度为37°44′40″的东西向5 km区域做时序分析，形成沉降漏斗剖面图，如图8所示。从图中可以发现小店地区主要有2个沉降漏斗出现：在a处沉降漏斗区域，东西向约500 m，面积较小，然而仅仅200 m差异却有较大沉降起伏，西向高点年沉降速率为32.825 mm/a，东向高点年沉降速率为32.034 mm/a，而两者间的沉降最低点为47.71 mm/a，有每年15 mm的沉降差异。考察发现该处有数栋高耸建筑，且地处在b处沉降漏斗区域，呈现自距起点2.7 km处(也为a沉降漏斗东向沉降最高点)起，向东侧沉降愈发严重，在4.13 km处达到了小店地区沉降漏斗最大年沉降速率点，年沉降速率达50.901 mm/a。在4.7 km处有明显沉降浮动出现，根据距离、坐标、位置判断为田庄断层，与2014年刘媛媛等人进行太原市地面沉降监测后得出的田庄断层结论一致，之后剖面向东进入农用地范围，沉降稳定。

D地区-吴家堡地区，自1956年沉降观测以来，吴家堡地区一直为沉降中心，在1981～1989年年沉降速率达到了114 mm/a，随后逐步缓和，根据最新的2010～2013年的沉降观测数据显示为6.6 mm/a，几近稳定。在本次监测时间段内，吴家堡地区最大沉降速率为7.929 mm/a，最大沉降量为-23.787 mm。其余多数点年沉降速率在2～6 mm/a之间，整个地区稳定。

图7 C区域地表监测累计沉降量

图8 C区域沉降漏斗剖面

E地区-西张地区，资料查询发现该地区沉降不明显，且在2003～2010年的年沉降速率为2.3 mm/a，而本文结果显示该地区无沉降。

总之，根据太原市区年平均沉降速率图与分析结果显示，在2013～2016年期间，太原地表总计形成万柏林、南堰、小店三个沉降中心，吴家堡、西张地区地表沉降现象较为稳定，而小店地区沉降现象最为严重，应引起注意。

3.3 可靠性分析

2014年，刘媛媛[7]等人使用20景ENVISAT ASAR数据采用PS-InSAR技术对太原市区进行地面沉降监测研究，结果显示2010～2013年太原形成吴家堡、小店两个沉降中心，并且在小店地区的剖线上发现田庄断层，本文也在小店地区纬度为37°44′40″东西向的5 km剖线上发现田庄断层，且剖面趋势与刘媛媛等人研究结果基本一致。2015年，刘瑾[9]等人对太原地表各阶段的各沉降漏斗中心的地下水水位变化与地面沉降速率进行相关联性分析，北部区域基本停止沉降， 中部区域趋缓，南部区域呈恶化趋势。2018年，周艳萍[4]基于灰色Verhulst模型对山西太原地面沉降进行趋势分析，得出西张稳定，万柏林和下元沉降减缓，吴家堡变化不大，小店年均沉降速率为45 mm/a的结果。2018年，许慧鹏[10]等人基于水准测量对太原地面沉降进行时空特征分析与模拟，小店地区沉降加剧。2016年，李军的《太原市现今地面沉降特征分析》[11]中，GPS监测点TY07为观象台，2009年7月～2014年6月的累计沉降量为-115.5 mm，年平均速率为23.1 mm/a。众多论文数据或分析结果与本文分析结论一致，证明本文数据处理成果的可靠性。

4 结 论

采用PS-InSAR技术对覆盖太原市区的37景2013年6月～2016年6月的COSMO-SkyMed数据进行处理，得到该区域的形变速率图和时间序列累计沉降量图。结果显示在该研究期内太原共计形成万柏林、南堰、小店三个沉降中心，尤其注意的是小店地区沉降范围逐渐扩大，沉降量逐渐变大，且最大沉降量达-152.7 mm，年沉降速率达50.901 mm/a，而吴家堡、西张地区沉降现象较为稳定。