星载InSAR技术支持下的昆明地表沉降监测

麻源源，陈云波，左小清，麻卫峰，吴文豪(. 昆明理工大学国土资源工程学院，云南 昆明 65009； . 昆明市规划编制与信息中心，云南 昆明 650500； . 云南师范大学旅游与地理科学学院，云南 昆明 650000; . 湖南科技大学煤炭资源清洁利用与矿山环境保护湖南省重点实验室，湖南 湘潭 0)

昆明地面沉降与高大建筑群的成片开发及地质因素有关[1]。主要表现为在一定的范围内发生地面形变，形成多个沉降漏斗，影响了昆明城市规划和发展。昆明市是目前为止中国大陆西部高原地区唯一一个发生沉降现象的城市[2]。合成孔径雷达干涉测量(interferometry synthetic aperture rader，InSAR)技术，因大面积、非接触、全天候、高精度等特点，与传统的监测测量(GPS和精密水准)技术相比有着一定的优越性。在近几年发展的时间序列差分干涉测量技术中，小基线集(small baseline subset，SBAS-InSAR)技术和永久散射体雷达干涉技术(persistent scatters InSAR，PS-InSAR)克服常规D-InSAR(differential InSAR)在地面形变监测过程中受到失相干和大气延迟的影响，可以进行长时间监测及提高InSAR技术监测精度。针对昆明沉降现象，国内许多学者通过不同方法和技术进行了大量研究，文献[3]通过二等水准测量对昆明市区在1986—1987年期间的垂直沉降进行了分析。文献[4]从矿物学和土质结构的角度对昆明在1987—1998年期间地面沉降进行了分析。文献[5]利用历史资料对1979—1998年期间昆明地面沉降发展过程及其特征进行了研究。文献[6]利用ALOS PLASAR数据基于SBAS-InSAR技术对昆明市区在2008—2010年的沉降进行了分析。文献[7]基于3S技术获取高精度的昆明地面沉降信息。本文利用欧空局发布的Sentinel-1A数据，提出基于升降轨SBAS-InSAR技术和数据融合技术对昆明在2014年12月至2017年3月期间的地面沉降进行分析，试验结果与上述学者所分析结果在形变趋势上存在一致性。呈现出以小板桥、东菊新村、河尾村、曹家村、罗家村为中心的多个沉降漏斗，而且以河尾村、曹家村、东菊新村和罗家村为中心的沉降漏斗有连成一片的趋势。其中沉降最严重的是以小板桥为中心的沉降漏斗，最大沉降速率达-54.2 mm/a。

1 小基线集(SBAS)技术原理

SBAS-InSAR技术先将时序SAR影像组成若干时空基线较短的干涉对集合，再利用奇异值分解方法将多个小基线集联合起来求地表形变速率最小二乘解，有效地解决了不同InSAR数据集间空间基线过长造成时间不连续的问题。SBAS-InSAR技术处理SAR影像时，除影像相位中形变相位外，还有区域地形误差带来的相位，由获取影像时大气不均质带来的大气相位异常、时空基线失相干现象和热噪声等引起相位误差，在差分相位干涉图中，两个高相干点的观测方程表示为[8]

(1)

2 SAR影像数据与处理

2.1 Sentinel-1数据来源

数据采用欧空局发布Sentinel-1A卫星数据[10]，Sentinel-1A雷达卫星于2014-04-03发射升空，是欧空局哥白尼计划发射的首颗环境卫星，该卫星数据实施免费对外开放政策。Sentinel-1A数据在IW和EW模式下采用TOPS(terrain observation with progressive)模式，该模式控制天线旋转方向与聚束模式正好相反，由后向前以固定速度旋转，使得所有目标在相同方位向天线模式下观测，降低扇贝效应和方位向变化模糊度[11-12]。Sentinel-1A数据具有精密轨道控制技术，采用交规盲区控制技术(across track dead-band control)确保任何重轨时卫星每个轨道点与地固参考轨道的位置误差均方根(RMS)在100 m范围内。试验覆盖区域如图1所示。

图1 试验区域覆盖示意图

利用Sentinel-1A卫星数据通过SBAS-InSAR技术可以对同一地区重访周期内地面沉降快速监测，非常适合昆明城区的地表形变监测，本文收集IW模式下Sentinel-1A数据，数据参数见表1。

表1 Sentinel-1数据参数

除此之外，卫星轨道数据由欧空局提供精密轨道数据和美国宇航局提供的SRTM3 DEM，分辨率为90 m，精度为±30 m，用于去除地形相位。由于升降轨数据开始日期和结束日期均很接近，为升降轨模式下的观测数据融合提供了可能性。

2.2 SBAS-InSAR主要流程

(1) 生成连接图。基于空间基线阈值ΔB=200 m和时间基线阈值ΔT=550 d为基线条件，共选取升轨232对干涉数据和降轨227对干涉数据，其中降轨数据以获取日期2015-10-16为主影像，升轨数据以获取日期2015-10-02为主影像。Sentinel-1A数据升降轨干涉数据时空间基线分布如图2所示。

图2

(2) 影像配准：将所有辅影像进行配准并重采样至主影像。Sentinel-1A卫星采用TOPS模式，配准精度要求达到千分之一像素级，相当于方位向零多谱时间配准真误差小于2 μs，由于地形起伏会对配准造成影响，多项式难以准确地拟合出主副影像几何偏移畸变的映射，因此利用几何配准误差为系统误差的特点采用哨兵卫星精密轨道进行初配准，再利用增强谱分集算法(enhanced spectral diversity，ESD)进行精配准来消除片段间出现的相位跳变现象[13-14]。最后进行方位谱去斜，达到干涉需要的配准精度获取正确的干涉图。

(3) 对干涉对进行集合D-InSAR处理：包括干涉对生成、去地平效应、干涉图滤波、相干性计算和相位解缠。干涉图滤波采用Goldstein方法进行干涉图滤波，相位解缠方法采用Delaunay三角网最小费用流量法(minimum cost flow，MCF)方法进行三维相位解缠时，升降轨γx均选取0.35作为阈值选取PS候选点。

(4) 轨道精炼和重去平：在参考DEM不准确的情况下，需要用地面控制点(GCP)作为稳定参考点，从这些点中计算出误差相位从而去除并且修正SAR数据。使用同一组GCP数据修正升降轨干涉图，目的是减少两个图像间相位误差和使用不同GCP带来的误差。

(5) SBAS-InSAR数据反演。估算形变速率和残余地形，基于线性模型计算所有像对形变。线性模型相比其他模型无需密集连接图和高相干性即可得到可靠结果且稳定性较好。计算时间序列上位移，在得到形变速率基础上，本文通过时域高通滤波和频域低通滤波去除残余DEM误差、大气干扰与轨道误差影响。

(6) SAR数据地理编码和辐射定标：将上述结果编码到GCS-WGS-84坐标系。

2.3 数据融合简化模型

升降轨SBAS-InSAR形变速率可以区分垂直方向和水平方向的形变速率。由于升降轨两组观测数据几乎从相反的两个方向获得。对于地表相同目标而言，东西方向水平形变速率在两组观测数据中的表现是基本相反的变化，但是垂直方向上则表现为基本相同形变速率变化。视线向上的形变速率量(DLos)的构成为

DLos=(UNsinφ-UEcosφ)sinθ+UVcosθ

(2)

式中，UN、UE和UV分别为南北、东西和垂直方向地表形变速率分量；φ为升降轨方位角；θ为入射角。哨兵-1雷达卫星的干涉宽幅模式(IW)，升轨入射角θ约为40°，降轨入射角θ约为39.5°。雷达视线方向(line of sight，LOS)形变对UN方向形变不敏感，视线向上的形变速率量表达式可进一步简化为

DLos≈±UNsinθ+UVcosθ

(3)

根据文献[12]可知，地表变形以垂直方向上形变速率为主，水平方向上的形变速率极为微小，即UNsinθ≪DLos。由于入射角在研究区的变化微弱，因此可用均值来估计地表垂直方向上的形变速率分量。由文献[15]可知入射角变化的简化对方位角的变化影响微弱，可以忽略不计。因而，可以推出升降轨观测模式下垂直方向的形变速率量公式[15-16]为

(4)

(5)

利用式(4)可将升降轨高相干点上雷达视线方向(LOS)形变转换为垂直方向地表形变速率分量UVA(D)，再根据式(5)求得升降轨数据融合后的垂直方向上速率分量UV，其前提是要求升降轨模式下同一观测区域具有相同的观测区域。升降轨下的昆明垂直沉降速率如图3所示。

图3

从图3(a)和图3(b)可知：两种模式下的沉降漏斗中心位置、漏斗形状、沉降区域大小及沉降速率大小几乎一致，均在-55～10 mm/a之间。

3 升降轨数据相互验证和数据融合

3.1 升降轨数据相互验证

由于雷达波入射角和入射方向的影响,同一相干目标在升降轨SAR观测下对象的位置略有差异[17-18]，为了对升降轨观测值进行相互对比验证，本文采用地形校正方法统一升降轨坐标系消除基准偏差。依据单个点目标进行直接比较难以实现，利用统计方法进行对比分析，从整体上验证升降轨在地表形变垂直形变的相关性。从升降轨地表垂直沉降速率图中提取11 072个同名点，其两组数据计算得到的昆明升降轨垂直沉降相关分析如图4所示。

图4 昆明升降轨垂直沉降相关分析图

由图4得知，两组数据线性拟合的方程系数K及两组数据相干性计算得到判定系数R×R、相关系数coherence及均方差根误差(RMES)[19]，见表2。

表2 升降轨垂直沉降速率关系

除此之外，图4中拟合的线性方程Y=0.962 5X-0.799 5与理想状态下的直线y=x方程对比可知两直线非常接近。因此试验结果表明两组数据具有较高的一致性和相关性，当SAR影像为30景以上，升降轨垂直沉降速率监测互检验精度优于2 mm/a。

3.2 升降轨数据数据融合

为了充分利用升降轨SBAS-InSAR技术监测结果和避免因雷达叠掩、透视收缩、阴影等引起失真现象[20]，本文对垂直沉降速率通过式(5)进行数据均值融合[21-22]。目的是提高昆明地表沉降监测精度和可靠性。图5为融合后的2014年12月至2017年3月昆明垂直沉降速率图。

图5 2014年12月至2017年3月期间昆明垂直沉降速率序列图

4 沉降分析

由图5分析可知昆明城区2014—2017年期间地面沉降分布的基本情况，昆明五华区沉降较慢，相对稳定，而昆明西山区和官渡区内的滇池北岸附近区域沉降较快，形成多个沉降漏斗。由于缺少同期的水准测量和GPS测量数据，无法准确地评价本文监测数据的精度，本文将对比不同历史时期[2-5]昆明城区的沉降速率及对昆明城区地表沉降的演变规律进行分析，不同时期昆明城区地表垂直沉降见表3。

表3 不同时期昆明城区地表垂直沉降 mm/a

试验结果分析：SBAS-InSAR技术监测昆明城区地面沉降分布与历史沉降监测数据得到的沉降漏斗分布是基本一致的，但也有新的情况变化。2014—2017年间沉降基本情况是五华区内的严家山一带和吴家营-大塘子一带，该地区较为稳定，与历史时期的沉降速率相比变化不大，沉降速率保持在-13 mm/a以内。其主要原因是该地区沉积物主要以岩石为主，其压缩性很小，地面形变量不大。而昆明地面沉降严重区域主要分布在官渡区和西山区，其中靠近滇池的区域出现大面积沉降，形成多个漏斗中心，形成以河尾村、曹家村、小板桥、罗家村和东菊新村为中心的5个沉降漏斗，5个沉降中心垂直沉降速率均达到-20 mm/a以上且有连成一片的趋势。并且罗家村-曹家村一带在历史时期并没有出现沉降现象，而在2014—2017年间却形成了两个沉降漏斗，沉降速率达到-25 mm/a以上。在几个沉降区域中，小板桥一带的沉降最明显,其最大平均地表沉降速率达-54.2 mm/a。说明目前昆明滇池附近和小板桥沉降情况较为严重，并且严重性逐渐增大，严重区域正在不断扩大。主要原因与该地区的地质因素有很大的关系，昆明地区的地基沉降主要以湖滨地带软土沉降为主，近年来滇池周围大量施工，高层建筑加压会导致软土的水被挤出，土层被压密实。加上工业用水剧增，导致承压水位降低，地层向更为致密状态变化而引起地基下沉。对滇池附近进行实地调查发现，滇池北岸周围存在道路地基错位、房屋倾斜，房体与地面分离、屋内地面下沉等现象，如图6所示。

图6 外业实地核查照片

外业实地调查结果表明，昆明地面沉降已经危害到公共安全和人民的居住安全。该调查也验证了InSAR技术可以快速、直观、有效地发现地面沉降严重区域及沉降中心，说明InSAR技术用于城市地面形变监测是可行的。

5 结 语

本文利用Sentinel-1A数据，基于升降轨SBAS-InSAR技术监测昆明城区地表沉降情况，获取升降轨两种不同模式下的昆明城区2014—2017年间地表沉降速率及城区的沉降漏斗分布情况，通过对两组观测数据进行相互验证表明卫星两种飞行模式在地面监测具有一致性。由于获取升降轨SAR数据的起始和结束时间十分接近，具有相近的观测时间段，因此对两种数据均值融合可以避免卫星单一轨道下观测数据失真。总体而言，SBAS-InSAR技术监测地表信息效果理想，与传统形变监测技术相比可以在短时间内得到大面积、高空间密度的地表形变场，具有成本较低无需布设地面控制点等优点。

参考文献：

[1] 郭乐萍，岳建平，岳顺.SBAS技术在南京河西地表沉降监测中的应用[J].测绘通报，2017(3):26-28.

[2] 姜朝松，邵德晟，樊友心,等.昆明市地面沉降发展过程及其特征[J].地震研究，2001，24(1):55-60.

[3] 黄崐，万建平，吴彬生.昆明地区近期地面垂直升降运动及分析[J].地震研究，1989，12(3):249-253.

[4] 薛传东，刘星，李保珠,等.昆明市区地面沉降的机理分析[J].中国地质灾害与防治学报，2004，15(3):47-54.

[5] 孟国涛.昆明南市区地面沉降研究[D].昆明：昆明理工大学，2003.

[6] 尹振兴，钟丽云，许兵,等.基于SBAS-InSAR的昆明地面沉降监测研究[J].地矿测绘，2016，32(4):1-5.

[7] 高照忠，罗志清，魏海霞.3S技术在昆明市地面沉降监测中的应用[J].地矿测绘，2007，23(1):27-29.

[8] BERARDINO P，FORNARO G，LANARI R，et al.A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，2003，40(11):2375-2383.

[9] 周志伟，鄢子平，刘苏,等.永久散射体与短基线雷达干涉测量在城市地表形变中的应用[J].武汉大学学报(信息科学版)，2011，36(8):928-931.

[10] 张艳梅，王萍，罗想,等.利用Sentinel-1数据和SBAS-InSAR技术监测西安地表沉降[J].测绘通报，2017(4):93-97.

[11] 吴文豪，周志伟，李陶，等.精密轨道支持下的哨兵卫星TOPS模式干涉处理[J].测绘学报，2017，46(9)：1156-1164.

[12] SOWTER A,AMAT MBC,CIGNA F,et al.Mexico City Land Subsidence in 2014—2015 with Sentinel-1 IW TOPS：Results Using the Intermittent SBAS (ISBAS) Technique[J].International Journal of Applied Earth Observations and Geoinformation,2016,52：230-242.

[13] 吴文豪.哨兵雷达卫星TOPS模式干涉处理研究[D].武汉：武汉大学,2016.

[14] 白泽朝，靳国旺，张红敏，等.天津地区Sentinel-1A雷达影像PSInSAR地面沉降监测[J].测绘科学技术学报,2017,34(3):283-288.

[15] 王涛，顾丽娟，詹华明,等.基于D-InSAR技术的天津地区地面沉降监测[J].测绘科学，2013，38(6)：23-28.

[16] 王艳，张玲，葛大庆,等.升降轨PSInSAR观测反演沉降与水平向位移试验[J].国土资源遥感，2014，26(4):97-102.

[17] 刘波，晏王波，范雪婷,等.基于SBAS技术的武进区地面沉降监测应用研究[J].现代测绘，2016，39(4):9-12.

[18] 王艳，葛大庆，张玲,等.升降轨PSInSAR地面沉降监测结果的互检验与时序融合[J].国土资源遥感，2014，26(4):125-130.

[19] 刘一霖.黄河三角洲地面沉降时序InSAR技术监测与地下流体开采相关性分析[D].青岛：中国科学院研究生院(海洋研究所)，2016.

[20] 丁琼，刘国祥，蔡国林,等.InSAR DEM精度与地形特征的关系分析[J].测绘科学，2009，34(1):147-148.

[21] 王琴，陈蜜，刘书军,等.利用升降轨道SAR数据获取DEM的试验研究[J].测绘通报，2015(6):39-43.

[22] 王霞迎，张菊清，张勤,等.升降轨InSAR与GPS数据集成反演西安形变场[J].测绘学报，2016，45(7):810-817.