利用短基线集InSAR技术监测抚顺市地面沉降

杨　帆，张　磊，张子文，赵增鹏

(辽宁工程技术大学测绘与地理科学学院，辽宁 阜新 123000)

抚顺煤田煤炭资源开发已有百余年历史，抚顺市经济的发展在很大程度上得益于煤炭资源的开采。但是目前来看，抚顺的两大露天开采区与地下采空区直接或间接的造成了城市内一些地质灾害的发生，如地面沉陷、地裂缝、滑坡、矿震等地质灾害，2014年1月，抚顺西露天矿矿坑南坡形成了一条约3100 m的地裂缝，并以每天8～10 cm的速度向北移动。这些地质灾害危及抚顺市人民的生命安全，制约着抚顺老工业基地的振兴和发展。由于多种因素的制约，抚顺市短期内还不能完全避免地质灾害，只能通过制定地质灾害防灾、减灾措施，降低地质灾害发生的可能性，提高其安全性。

合成孔径雷达差分干涉测量技术(differential interferometric synthetic aperture radar，D-InSAR)是近几十年发展起来的空间大地测量技术，具有全天候、无接触、大面积、高空间分辨率、高精度等优势。短基线集(small baseline subset，SBAS)技术是时间序列D-InSAR技术中一种常用的形变监测分析方法，能够克服时空失相干和大气效应的影响，监测长时间范围内地表形变，获得研究区域的沉降规律和演化特征[1]。与永久散射体InSAR(PS-InSAR)技术相比，SBAS技术对影像数量的要求较低，可以提高影像的利用率，在地面沉降监测方面得到了广泛的应用。李永生等利用SBAS反演结果，了解了北京的地面沉降空间分布特征[2]；刘志敏等采用SBAS技术对长治矿区地表形变进行监测，得到了研究区地面沉降的空间分布及时间序列的相对形变量[3]；邱志伟等基于短基线集技术反演出南京地区的地表形变场及累计形变量[4]；李珊珊等将SBAS技术应用到青藏高原季节性冻土的形变监测[5]；张金芝等将SBAS技术应用到现代黄河三角洲的地面沉降监测[6]；张永红等改进SBAS中的时序形变模型并将其应用到太原市矿区的地面沉降监测[7]。本文利用短基线集(SBAS)技术对覆盖抚顺市部分地区的12景COSMO-SkyMed高分辨率SAR数据进行了处理，获得了该研究区域的地面沉降分布和沉降速率图。

1　研究区概况和数据源

1.1　研究区概况

抚顺市曾经是我国最大的燃料生产基地，素有“煤都”之称，并在煤炭开发的基础上建立起来了相关行业，成为一座具有一定规模的综合型工业城市。其中心地理坐标为东经123°55′，北纬41°52′。抚顺地区以浑河水系为主体，市区位于抚顺西部浑河河谷冲积平原上，平均海拔65～99 m，呈东西走向，南北为山地，浑河由东至西将市区分割成南北两部分。市区的南北宽约6～8 km，东西长约30 km，地势从东向西逐渐降低。抚顺矿区位于抚顺市浑河南岸，长约18 km，宽约2 km，主要有西露天采场和东露天采场等，其中，西露天矿位于抚顺煤田西部，矿坑东西长6.6 km，南北宽2.2 km，总面积为13.2 km2，开采垂直深度388 m。东露天矿位于抚顺市区的南部，抚顺煤田的东部，矿坑东西长5.7 km，南北宽1.9 km，面积9.2 km2。其中抚顺东露天矿在2006—2009年3年的开发建设后，于2009年9月矿山进入正规化、集中化、合理化开采阶段。

1.2　数据源

本试验使用时间跨度为2014年4月至2015年3月、覆盖研究区域的12景降轨COSMO-SkyMed数据。该影像是由意大利航天局和国防部共同研发的卫星星座获取的，成像模式为HIMAGE条带模式，分辨率为3 m×3 m。所用影像数据的基本参数情况见表1。为了提高运行速率，对研究区域加以裁剪，只对研究区域进行数据处理，影像覆盖范围及裁剪研究区域范围如图1所示，外层细线框为影像覆盖范围，内层粗线框为研究区范围。

DEM数据来自美国宇航局(NASA)官方网站，是美国太空总署(NASA)和国防部国家测绘局(NIMA)联合测量的覆盖研究区域的SRTM数据。所使用的DEM数据版本号为4，分辨率为90 m×90 m，用以模拟地形相位，剔除干涉图的地形相位成分[8]。

2　研究方法与数据处理

2.1　短基线集技术

短基线集技术(SBAS-InSAR)的基本原理就是将单次DInSAR得到的形变结果作为观测值，再基于最小二乘法则来获取高精度的形变时间序列[9-10]。为了抑制由于时空基线过长所引起的时空失相关现象，短基线集技术从获取的SAR影像中选择时空基线均小于一定阈值的干涉影像对生成多视差分干涉图，这些干涉影像对会根据基线情况分成若干个集合，每个小集合内用最小二乘法进行解算，求得地表形变时间序列；小集合间利用奇异值分解法(SVD)将多个小基线集联合起来进行求解，得到整个时间段的形变时间序列[11-12]。

表1　CSK影像参数

图1　研究区域与影像覆盖范围

假设同一地域有N+1幅雷达影像，(t0,t1,…,tn代表影像的获取时间)，并假设所有雷达影像均已配准到同一坐标系下，根据干涉条件组合，得到M幅干涉图，其中M满足

在进行去平地效应及DEM相位后，在(x,r)点由tB和tA两个时刻的影像生成的解缠相位，可以得到

(1)

式中，j=1,2,…,M；Δddisp为形变相位差；Δz为DEM误差；Δdatm为大气引起的相位差。

将DEM误差和大气相位误差去除后，可以将式(1)写成如下形式

Bφ=δφ

(2)

式中，B为M×N矩阵；φ表示未知形变相位矢量。同时可以将式(2)的格式进行转换，用两个时间点之间的平均速率来代替其中的未知相位参数，得到

Bυ=δφ

(3)

式中，B表示M×N矩阵。如果SAR数据分割到几个独立SBAS，会造成矩阵B秩亏，导致式(3)有无穷解。然而，在SBAS处理流程中利用SVD方法，对矩阵B进行伪逆运算，就可以得到v的广义逆解。

对式(3)进行SVD运算求出v之后，再对各时间段的速度进行时间域上的积分得到每个时间段上相位的形变值，然后将相位形变值乘以相位形变转换系数(λ/4π)就可以得到整个监测周期内地表在雷达视线方向上的形变时间序列。

2.2　数据处理

本试验使用SARscape雷达影像处理软件，以ENVI5.3为数据处理平台，对覆盖研究区域的12景COSMO影像数据进行预处理。选取比较小的时间和空间基线阈值来抑制时空去相干的影响[13]，设置时间基线阈值为300 d、空间基线阈值为400 m，对12景SAR影像进行自由组合，共产生了26个干涉对。选取成像时间为2014年9月24日的影像为公共主影像，其中空间基线最大约为356 m，最小约12 m，时间间隔最大为272 d，最小为16 d。在干涉图生成阶段对干涉图进行了方位向视数为2、距离向视数为2的多视处理，多视处理后所对应地面分辨率为5 m、利用辅助DEM数据模拟地形相位，从干涉相位中去除地形相位，生成差分干涉图，并使用Goldstein方法对差分干涉图进行滤波[14]，采用3D解缠(Delaunay MCF)算法进行相位解缠[15-16]。处理流程如图2所示。

3　结果与分析

3.1　试验结果

通过对覆盖研究区域的12景影像数据进行处理，得到了研究区域的沉降信息。图3为获得的研究区域的沉降速率场。从图3可以看出，抚顺市在2014年4月—2015年3月监测的时间段内整体存在沉降，大部分沉降速率在-25～-45 mm/a的范围内，沉降严重的区域为新抚区。抚顺西露天矿南侧和东露天矿干涉生成的相干点比较少，现将20140924与20141010两景影像的干涉图与沉降速率图叠加，得到叠加效果图(如图4所示)，图4中显示位置为东露天矿，圆圈内表示在16 d内发生了大的形变的位置，正好与沉降速率图中未得到沉降信息的地方重合，但是从边缘干涉出的点可以看出有比较明显的沉降，这便是造成失相干的原因。

图2　SBAS技术处理流程

图3　抚顺市2014年4月—2015年3月年平均沉降速率分布

为了对研究区进行进一步精细化分析，将栅格格式的沉降速率图转化为矢量格式的沉降速率图，并利用克里金插值算法[17]对其插值，得到面状沉降信息分布图(如图5所示)。从图5中可以看出，抚顺西露天矿南侧、东露天矿内均出现了不同程度的沉降，对沉降最严重的区域进行分析，作A-A′剖面得到其剖面线(如图6所示)，发现剖面线存在两个比较明显的极小值点，分别达到-169和-125 mm/a，形成两个小沉降漏斗，随着时间的推移，两个沉降漏斗有可能合并成一个。此处沉降漏斗的形成与东露天矿采场的生产活动有关。

图4　叠加效果(东露天矿)

图5　研究区域面状沉降信息分布图

图6　A-A′剖面线

3.2　精度分析

利用收集到的4个水准点的观测值对SBAS监测结果进行对比验证，结果见表2。以水准点观测值为形变值的真值，4个水准点的最大误差为6.254 9 mm，最小为2.141 3 mm，误差均小于10 mm。可以认为SBAS-InSAR技术监测地表形变能达到亚厘米级的精度。产生误差的主要原因是水准点与得到的PS点位置存在偏差，不能完全对应，采用了以水准点为中心，将一定范围内的PS点求取的平均值作为SBAS监测值，与水准点进行对比的方法。

表2　水准点与SBAS监测结果对比

4　结　语

本文采用短基线集技术对覆盖抚顺市的12景COSMO-SkyMed影像数据进行干涉处理，得到了抚顺市在研究周期内的地面沉降信息。抚顺市整体存在沉降现状，并且有2个沉降严重的区域，分别位于抚顺西露天矿南侧和东露天矿西北侧，并且有扩大的趋势，相关部门应对这两个地点加强监测，掌握其形变动态，及时采取防治措施。通过试验结果进一步分析可知，SBAS-InSAR技术应用于地面沉降的监测，精度可以达到亚厘米级，但是对于大量级的形变容易造成失相干，从而得不到监测值。

参考文献：

[1] 龙四春.DInSAR改进技术及其在沉降监测中的应用[M].北京:测绘出版社,2012:1-3.

[2] 李永生,张景发,李振洪,等.利用短基线集干涉测量时序分析方法监测北京市地面沉降[J].武汉大学学报(信息科学版),2013,38(11):1374-1377.

[3] 刘志敏,李永生,张景发,等.基于SBAS-InSAR的长治矿区地表形变监测[J].国土资源遥感,2014,26(3):37-42.

[4] 邱志伟,岳建平,汪学琴,等.基于短基线集技术的城市地表沉降监测研究[J].测绘通报,2016(7):25-29.

[5] 李珊珊,李志伟,胡俊,等.SBAS-InSAR技术监测青藏高原季节性冻土形变[J].地球物理学报,2013,56(5):1476-1486.

[6] 张金芝,黄海军,毕海波,等.SBAS时序分析技术监测现代黄河三角洲地面沉降[J].武汉大学学报(信息科学版),2016,41(2):242-248.

[7] 张永红,吴宏安,孙广通.时间序列InSAR技术中的形变模型研究[J].测绘学报,2012,41(6):864-869.

[8] 姜德才,张继贤,张永红,等.百年煤城地表沉降融合PS/SBAS InSAR监测——以徐州市为例[J].测绘通报,2017(1):58-64.

[9] 张子文,杨帆,吴文豪,等.地下水开采与地面沉降关系的短基线集分析[J].测绘科学,2016,41(6):64-69.

[10]SAMSONOV S，D’OREYE N，SMETS B.Ground Deformation Associated with Post-mining Activity at the French-German Border Revealed by Novel InSAR Time Series Method[J].International Journal of Applied Earth Observation & Geoinformation，2013，23(8):142-154.

[11]周吕,郭际明,李昕,等.基于SBAS-InSAR的北京地区地表沉降监测与分析[J].大地测量与地球动力学,2016,36(9):793-797.

[12]孙晓鹏,鲁小丫,文学虎,等.基于SBAS-InSAR的成都平原地面沉降监测[J].国土资源遥感,2016,28(3):123-129.

[13]YANG C S，ZHANG Q，ZHAO C Y，et al.Monitoring Land Subsidence and Fault Deformation Using the Small Baseline Subset InSAR Technique：A Case Study in the Datong Basin，China[J].Journal of Geodynamics，2014，75(4):34-40.

[14]吴文娟.矿区SAR干涉图滤波方法研究[D].徐州:中国矿业大学,2014.

[15]占文俊,李志伟,韦建超,等.一种InSAR大气相位建模与估计方法[J].地球物理学报,2015,58(7):2320-2329.

[16]胡争,谢荣安.SBAS-InSAR技术在地面沉降监测中的应用——以深圳市为例[J].地矿测绘,2013,29(2):10-13.

[17]何亚乐,朱琳,宫辉力,等.TerraSAR的首都机场形变特征分析[J].测绘科学,2016,41(12):14-18.