基于PS-InSAR和SBAS技术监测昆明市主城区地面沉降

2020-09-22姜乃齐左小清王志红赵永恒

贵州大学学报(自然科学版) 2020年4期

姜乃齐，左小清*，王志红，游洪，赵永恒

(1.昆明理工大学国土资源工程学院，云南昆明 650093；2.贵州工程应用技术学院，贵州毕节 551700；3.云南省测绘工程院，云南昆明 650033)

地面沉降是一种由于人类工程经济活动影响和地下松散固结压缩而导致的地壳表面标高下降的地质现象，又称地面下沉或低陷。它一般有区域性下沉和局部性下沉两种表现形式，引起城市地面沉降的原因与大量开采地下水、工程活动等有密切关系，可引起建筑物倾斜，地基的稳定性被破坏等，给生产和生活带来很大的影响。昆明市作为云南省的省会城市，也是目前为止我国西南地区唯一发生沉降的城市[1]，其主城区地面沉降现状备受关注。近几年，昆明市正处于城市发展的重要时期，例如地铁建设、建筑施工等一系列城市建设活动都对地面沉降产生很大的影响。传统的地面沉降监测方式有水准测量、全球定位系统测量等，相对于这些传统的变形监测手段，合成孔径雷达差分干涉测量(differential interferometry synthetic aperture radar,D-InSAR)能更加有效地对目标进行全天候、全天时、大范围、高精度的监测。2001年，FERRETTI等提出了永久散射体雷达干涉测量(permanent scatterers interferometry synthetic aperture radar,PS-InSAR)技术[2]，并且在2003年的Fringe大会上被确认为InSAR领域最具有应用前景的技术之一[3]。2002年，由BERARDINO和FORNARO等研究人员提出了短基线集干涉测量(small baseline subset znterterometry synthetic aperture radar, SBAS-InSAR)技术与PS-InSAR技术不同策略的时间序列分析方法[4]。PS-InSAR技术和SBAS-InSAR技术这两种方法的提出有效地解决了传统D-InSAR易受时空失相关与大气效应的影响，能充分利用多景SAR影像和时序分析技术提高InSAR监测地面的精度。

2003年，孟国涛[5]基于1987—1998年四期水准数据分析，研究了昆明南市区地面沉降情况，表明以小板桥为中心的沉降区累计沉降量最大且有明显扩张的趋势；2007年，高照忠等[6]采用以“3S”技术为基础，建立GIS自动监测系统监测了昆明市的地面沉降情况；2016年，尹振兴等[7]利用SBAS-InSAR技术基于2007—2010年SAR影像对昆明市地面沉降进行监测，发现了新的沉降漏斗；2017年，邵九明等[1]采用PS-InSAR技术研究昆明市地面沉降，昆明市地面沉降区域主要集中在南市区及滇池沿岸且空间分布不均，地面沉降总体呈现放缓趋势；2018年，谢文斌等[8]基于SBAS-InSAR技术对昆明市地表沉降进行监测研究，得出了传统沉降中心沉降趋于缓和并发现了滇池会展中心片区形成了新的沉降带。本文选取2018—2019年30景降轨的Sentinel-1A数据，分别利用PS-InSAR和SBAS-InSAR两种技术监测昆明市的地面沉降并作对比分析。

1 基本原理

1.1 PS-InSAR原理

PS-InSAR技术的原理是利用覆盖同一地区的多景(一般要求大于25景)SAR影像，通过统计分析影像的幅度信息，生成时间序列干涉像对，提取出不受时空基线失相关和大气效应影响的永久散射体点，最后分离出这些目标点上的地形相位以此来监测地面沉降。在城市区域可以获得大量的稳定点作为PS候选点，例如人工建筑物等。

假设获得了同一研究区域M+1景不同时间段的单视复数(single look complex, SLC)影像，以多普勒质心频率以及时间基线为依据选出主影像，剩余的M景影像则为辅影像，主影像分别与辅影像进行配准、干涉处理，然后引入外部数字高程模型(digital eleavation model, DEM)，通过二次差分法去除地形相位，生成去平地效应和地形效应后干涉图。其中，每个像素对应的差分相位由式(1)表示。

φ=φdef+φtopo+φatm+φorb+φnoi。

(1)

式中：φ为每个像素的差分相位；φdef为形变相位；φtopo为地形相位；φatm为大气相位；φorb为轨道误差相位；φnoi为噪声相位。其中φtopo和φatm是由于DEM精度问题和大气延迟导致的。本文通过振幅离差指数阈值法来筛选稳定的PS点，振幅离差指数阈值法[9]是由FERRETTI提出，表示为式(2)。

(2)

1.2 SBAS-InSAR原理

SBAS-InSAR技术的原理是将时间序列雷达数据整合成多个小基线集，基于形变速率最小范数准则，利用矩阵奇异值分解(singwlar valne decomposition, SVD)方法获取目标的形变速率及其时间序列[10]。覆盖同一区域范围时间序列为t1,t2,t3，……，tN的N+1幅SLC影像，生成的差分干涉图数量M满足不等式(3)：

(3)

假设由t1时刻和t2时刻的两幅影像干涉获得了第i幅去除地形相位后的干涉图，且t2时刻大于t1时刻，则第i幅干涉图的干涉相位可表示为式(4)：

(4)

(5)

则进一步可得第j幅相位值可得：

(6)

把上式写成矩阵形式可得：

BV=δφj。

(7)

式中：B为M×N的矩阵，当M

2 实例分析

2.1 研究区域及试验数据

昆明市位于云贵高原的中部，总体地势为北高南低，并呈阶梯状逐渐递减[12]，其地区面积约为2 000 km2，其中沉降严重的区域大约为300 km2。由于近几年昆明市轨道交通的飞速发展、城市内高楼的兴建和地下水的抽取等原因，市区内出现大规模的地面沉降现象且有新的沉降中心产生，沉降面积不断扩大。因此研究区域范围选取如图1所示，包括西山区、盘龙区、五华区、官渡区和呈贡区，中心经纬度为25°00′北，102°43′东。本文从欧空局(european space agency，ESA)下载了30景IW模式下Sentinel-1A降轨影像，时间跨度为2018年7月20日至2019年7月27日，极化方式为VV。Sentinel-1A卫星是欧空局发射的C波段雷达遥感卫星。此卫星的默认成像模式为渐进式地形观测模式(terrain observation with progressive，TOPS)此模式提高了影像的质量和干涉性能[13]。数据参数如表1所示。

另外，为了提高影像轨道精度，本文从欧空局获取了星载多谱勒无线电定轨定位系统(delft object-oriented radar interferometric software,DORIS)精密轨道数据，美航局提供的DEM，用于去除地形相位，拼接后的DEM如图2所示。

图1 试验区域范围图2 昆明地区的DEM图 Fig.1 Location of the study area Fig.2 Digital elevation model map of Kunming area

表1 Sentinel-1A数据参数表Tab.1 Table of sentinel-1A data parameter

2.2 数据处理

该试验主要基于SARscape软件和30景Sentinel-1A数据，具体试验处理流程包括SBAS干涉处理和PS干涉处理。

2.2.1SBAS干涉处理

(1)数据预处理：读取30景原始数据，利用矢量文件裁剪出研究的区域。

(2)连接图和连接相对的生成：选取日期为2019-02-29的影像作为超级主影像，设置最大临界基线为2%，最大时间基线阈值设置为250 d，最终生成308个像对，像对时空基线连接图如图3所示，最大时间基线为240 d，最小时间基线为12 d。最大空间基线为100.307 m,最小空间基线为3.173 m。

图3 SBAS干涉处理时空基线分布Fig.3 SBAS spatiotemporal baseline distribution

(3)干涉工作流：将308个像对进行干涉处理，相干性生成，去平、滤波和相位解缠，把所有的数据对都配准到超级主影像上为轨道精炼和重去平以及SBAS反演做好准备。本文选用的解缠方法为Delaunay MCF,该方法能获得全局最优处理结果，滤波方法选择Goldstein。

(4)轨道精炼和重去平：用于估算并去除残余相位，由于卫星轨道不准确或DEM地理定位不准确，用地面控制点(ground control point, GCP)来修正SAR数据，其中，选择GCP需要参考的标准有：1)选择没有残余地形条纹的区域；2)选择没有形变的区域，除非知道该点的形变速率；3)选择没有相位跳跃的点。本文一共选择了50个GCP。

(5)SBAS一次反演：是反演的核心，估算形变速率和残余地形相位，并通过二次解缠来对输入的干涉图进行优化。

(6)SBAS二次反演：SBAS二次解缠的干涉相位图残余相位分离为低通和高通两部分相位，分别用低通滤波和高通滤波对残余相位滤波，获得非线性形变，把线性形变相位和非形变相加，即所得地表精确的形变信息。

(7)地理编码：将地表形变结果投影到地理坐标系上，本文将结果投影到视线方向(LOS)上。

2.2.2PS干涉处理

(1)数据预处理：与SBAS处理一致。

(2)连接图生成：同样选取日期为2019-02-29的影像为主影像，生成连接图如图4所示,最大时间基线为204 d,最小时间基线为12 d,最大空间基线为98.950 m,最小空间基线为5.537 m。

图4 PS干涉处理时空基线分布Fig.4 PS spatiotemporal baseline distribution

(3)干涉工作流：对30景SAR影像进行配准去平差分干涉处理生成29个干涉对，通过振幅离差指数筛选出稳定的PS点，一共2 230 055个有效的PS点，且这些PS点大多分布在建筑物、桥梁等，滇池、绿化区域，基本没有PS点。

(4)PS第一次反演：获得位移速率和残余地形来对合成的干涉图进行去平，主要方法是基于识别一定数量的“相干雷达信号反射”(永久散射体)，要求永久散射体必须是稳定的(变动小于1 mm)且能被SAR天线探测到。

(5)PS第二次反演：这一步去除大气相位得到最终的形变速率。

(6)地理编码。

3 结果分析

3.1 沉降结果分析

经过上述处理过程，得到昆明市主城区地表的年沉降速率如图5所示，可以看出两种方法的沉降漏斗中心位置、沉降速率及沉降区域大小几乎一致，其中A、B、C三个区域为本次监测中沉降相对严重的区域，如图6所示。从图6中可以看出昆明的地表沉降在空间上呈不均匀分布，其中昆明市的地表沉降区主要分布在滇池东北部、东部地区以及昆明中部地区。

图5 PS与SBAS年均沉降速率Fig.5 PS and SBAS average annual settlement rat

图6 沉降严重的区域Fig.6 Area with serious settlemen

A区位于昆明官渡区，出现了多处沉降漏斗：首先是螺蛳湾-广卫村-义路村-羊甫村-枫林盛景小区-星宇园-子君欣景等片区，PS和SBAS两种方法的最大年沉降速率分别为-36.822 mm/a和-33.686 mm/a，沉降面积为5.63 km2；其次是昆磨高速公路和东环城高速公路交界附近，最大年沉降速率分别为-35.275 mm/a和-32.685 mm/a，这两处沉降由于居民地和交通网络密集，土层以软黏土土为主，且附近有昆明地铁四号线在建，地铁隧道的施工、地下水抽取以及施工降水措施以及交通流量和建筑用地施压等造成了土层水分渗出、土体变密、地面下陷等，一系列原因共同造成了这一片区的沉降[14]。小板桥位于官渡区的东南部鸣泉村立交桥西南处附近，最大年沉降速率分别为-23.874 mm/a和-21.814 mm/a，另外从世纪金源购物中心至望春苑一带靠近地铁1号线饵季路站已连成一片，年沉降速率分别为-17.137 mm/a和-14.726 mm/a，根据文献[15]和[16]可知在2014—2017年晓东村沉降严重，本次监测可以得出晓东村一带沉降趋势有明显的缓和，最大年沉降速率分别为-13.171 mm/a和-11.303 mm/a。

B区位于官渡老城区，临近滇池，这一片区出现了多个沉降漏斗并有向市中心扩张的趋势，靠近滇池由北向南分别为小河嘴村、苏王村、滇池国际会展中心一带、龚家村、丁姚村、六甲村、罗家村、五甲河、龙马村、罗衙村等沉降区域，其中沉降较为严重的几块区域：六甲村-罗家村-六甲河一带，最大年沉降速率分别达到-38.679 mm/a和34.501 mm/a，小河嘴村-滇池国际会展中心-龚家村-韩家村一带，最大年沉降速率分别为-38.716 mm/a和-35.482 mm/a，这一沉降区域靠近滇池，滇池附近的土层松软且饱含水分，并夹有松软沉积物，如淤泥等，这些地层孔隙大且含有大量水分使得周边地区的地基松软，并且有大量的地下水向滇池方向径流和排泄，周围的软土下沉导致了这一片区的沉降现象[15]。此外，靠近清水河湖畔之梦小区的南绕城高速公路一带，仟村佳宇小区周边以及润城附近都有明显的沉降现象，这一区域最大年沉降速率分别为-32.070 mm/a和-29.265 mm/a。

C区位于昆明西山区，可以看出棕树营一带出现沉降，最大沉降速率分别为-21.258 mm/a和-18.801 mm/a。

3.2 结果对比分析

3.2.1相关性分析

由于缺少同期水准数据，为了进一步对两种方法进行精度分析，分别在PS和SBAS两种方法的沉降速率图上提取同名点，这些同名点具有相同的地理坐标，在研究区域均匀筛选1 010个相干性较高的同名点，对比PS和SBAS两种方法所得结果的年平均沉降速率，得到的相关系数图如图7所示，两者的线性相关性R达到了0.997 9，且PS点和SBAS点的沉降速率分布均匀，可以得出两种方法得到的年平均沉降速率具有一致性，说明这两种方法对地面沉降监测的可行性。

图7 PS与SBAS年均沉降速率相关系数图Fig.7 Correlation coefficient diagram of annual average settlement rate of PS and SBAS

3.2.2时序分析

分别从A、B、C三个区域选择在沉降漏斗区域选出了三对典型点进行分析，如图8所示。A区典型点位置选择在昆磨高速公路旁义路村附近(a)，B区典型点选择在滇池国际会展中心一带靠近苏王村(b)，C区典型点选择在棕树营(c)，绘制出这三个典型点在2018年7月到2019年7月时间跨度内的形变曲线。

由图8可知PS和SBAS两种方法下的典型点沉降情况和趋势基本一致，但是得到的结果还是有一些细微的差别：图8(a)中，2019年6月左右，SBAS点的沉降呈现放缓趋势，而PS点的沉降有明显的加剧趋势；图8(b)中，SBAS点总体呈持续沉降加剧的趋势，而PS点的沉降情况比SBAS点更加剧烈，且在某些时间段有沉降放缓的情况，例如2018年9月18日至10月、11月17日至12月等；图8(c)中，在2018年10月至11月、2019年5月28日至6月初等时间段，SBAS点呈现出加剧的沉降现象，PS点的沉降却明显放缓，在2018年7月至8月、2019年4月至5月中旬等期间，SBAS点的沉降情况有缓和的趋势，而PS点呈现加剧的趋势。出现这些细微的差别是由于这两种方法选择备选点的算法不同导致的[17]。