基于SBAS-INSAR的德兴矿区沉降监测与分析

周子琪 周世健 欧阳双艳

(1. 东华理工大学 测绘工程学院, 江西 南昌 330013； 2. 南昌航空大学 环境科学与工程学院, 江西 南昌 330063)

0 引言

改革开放至今,我国国民经济高速发展的同时,对矿产资源的需求量也逐渐增大。各类矿产资源的大力开发对工业发展和社会经济有着良好的促进作用[1],但开采过程可能会对矿区的地质结构带来较大影响,甚至产生一系列地质灾害问题,最易产生塌陷、地表沉降等现象,更为严重会导致地下水资源污染、岩体变形诱发崩塌、滑坡以及泥石流等地质灾害的发生,严重影响当地地质结构和周围居民的人身财产安全[2]。因此,为了保障矿区的可持续开采以及周边居民的正常生活,除了选择合理的开采方式以及减少矿区的开采量外,对矿区进行沉降监测分析也显得尤为重要。如今,通过技术手段获取地面沉降信息及其分布特征可以及时有效采取对应措施,在一定程度上可以降低发生地质灾害的风险并开展防治工作[3-5]。

在近几十年大面积的沉降监测研究当中,沉降监测的方法也发生了变化,用于地表观测的方式、手段也趋于多样化,以往传统观测的方法包含:精密水准测量、三角高程测量、静态全球导航卫星系统(Global Navigation Satellite System，GNSS)测量和动态GNSS监测等[6],这类传统的监测方式虽有较高的测量精度和可靠性,但对于大面积、需要长期动态监测的区域而言,这些散点式的监测方式不能很好地满足实际的工程应用需要。为了满足大面积的沉降监测需要,近年来,合成孔径雷达干涉测量(Interferometric Synthetic Rader,InSAR)的快速发展弥补了传统变形监测对于大面积监测的不足。InSAR作为新兴遥感技术,其中合成孔径雷达差分干涉测量(Diffe-rential Interferometric Synthetic Aperture Rader,D-InSAR)技术具有连续覆盖地表、高精度、自动化程度高的优点[7-8],已成为获取表面变形信息的有效手段,其精度可以达到厘米甚至毫米级别[9-11]。1989年,Gabrie等首次应用D-InSAR技术获取了加州ImperialValley灌溉区的影像数据,首次说明D-InSAR技术可以探测厘米级的形变信息[12]。但由于D-InSAR技术易受时空失相干性以及大气延迟效应的影响,应用的领域受限[13-15]。为了对D-InSAR技术进行优化,2002年BERARDINO等提出了小基线集技术(SBAS-InSAR),该技术能够针对时间跨度长、规模大的区域进行分析[16]。2011年江利明获取了乌达的34景影像,其研究结果显示，利用SBAS-InSAR技术处理出来的结果与GPS(Global Positioning System)结果一致[17],证明广域差分增强系统(Satellite-Based Augmentation System，SBAS)-InSAR技术在地面沉降监测应用的可靠性以及与其他技术监测结果的一致性[17-18]。本文利用SBAS-InSAR技术,对德兴矿区 2018年12月28日 至 2019年 12月24日 期间的30景SAR影像进行处理,得到该区域的沉降速率和沉降分布,对沉降严重的区域进行分析,此研究结果可以为该区域地质灾害的防治问题提供一定帮助。

1 研究区域概况及实验数据

德兴矿区位于江西省上饶市德兴境内,距德兴市22 km,其矿区位置位于117°43′40″E,29°01′26″N处。此次实验研究区的地势呈东南部偏高、西北部偏低,德兴矿区地理位置属于江南丘陵地区,其整体地貌起伏较大,海拔高度约为65～500 m,矿区内地面切割较为强烈,山坡较陡,坡度在30°左右。德兴矿区内矿产资源种类丰富,存矿量较大,其中大型和特大型铜金矿高度集中,产业特色明显,该矿区是我国的矿业经济中最为发达与发展潜力最为雄厚的地区之一,同时也是国家级重要有色贵金属能源基地之一。德兴矿区周围分布有祝家社区、富家坞社区、詹家坞社区以及梨园社区,矿区的开采在一定程度上间接影响了当地居民的生活环境。同时,德兴矿区在长期的挖掘开采过程中,积累了一定程度的矿区地质问题,这些问题会给矿区从业者和周围村民带来潜在的安全隐患。

本文实验数据选取覆盖德兴市矿区的30景升轨Sentinel-1A数据,Sentinel-1A数据,其时间跨度为2018年12月29日至2019年12月24日,影像时间间隔为12 d。本实验使用的Sentinel-1A数据均为C波段,波长为5.66 cm,入射角大约为39.3°。采用由美国航空航天局(National Aeronautics and Space Administration, NASA)提供的30 m分辨率航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission，SRTM)1数据作为外部参考数字高程模型(Digital Elevation Model，DEM)以消除地形误差对实验的影响。Sentinel-1A数据参数见表1。

表1 Sentinel-1A数据参数

2 SBAS-InSAR基本原理

SBAS-InSAR是以传统D-InSAR技术为基础再进行研究,用来获取工作区地表形变的时间序列图。该技术选取N+1幅SAR影像,通过设置合适的时间基线和空间基线生成N个小基线集差分干涉对,去除地形相位后生成差分干涉图并且进行相位解缠,将常规D-InSAR监测的观测结果用作单个观测值,采用奇异值分解(Singular Value Decomposition,SVD)的方法将每一组小基线数据集连接起来,解决时间上采样过于稀疏的问题,又结合稳定散射体的干涉相位信息,得到更高的空间分辨率,然后根据最小二乘法求出高精度的变形时间序列。

本文实验假定在感兴趣区域的有序时间(t1,t2,…tN)内,采集研究区内同轨道上的一组N+1景SAR图像,根据干涉条件,将辅影像配准到主影像上,在配准完成后根据连接图生成的干涉图进行复相位运算。通过设置合理的时间和空间基线生成M个差分干涉图。其中M需要满足以下条件[19]

(1)

在干涉步骤中,在tA和tB时刻(tA>tB)获取的两幅SAR影像作为主辅影像生成的第j幅差分干涉图,第j幅差分干涉图中任何像元的干涉相位可写成

(2)

式中,λ为雷达波长;d(tA,x,r)和d(tB,x,r)是在tA和tB时刻相对于参考时间t0的雷达视线方向的累积形变;Δφjtop是由所参考的DEM不精确而引起的地形相位误差;Δφjatm为大气延迟相位误差;Δφnoise是噪声影响的相位误差。

为去除地形相位误差、大气延迟误差以及噪声相位误差,通过简化式(2),干涉相位可以表示为

(3)

为了获得含有物理意义的沉降序列,将式(3)中相位与时间之比来表示获取影像时间中的平均相位速度

(4)

其中,第j幅干涉图的相位值可以表示为

(5)

用SBAS-InSAR算法计算时间序列的变形可表示为

AV=δφ

(6)

其中,AV为M×N矩阵;δφ是代表干涉相位值的向量。最后使用SVD方法用于获得最终形变速率。

SBAS-InSAR技术处理流程图如图1所示。

图1 SBAS-InSAR技术处理流程图

3 数据处理

3.1 数据处理平台

实验使用GAMMA软件进行SBAS-InSAR处理。GAMMA软件可以完整地处理雷达信号,其中的功能模块包括:从合成孔径雷达(Synthetic Aperture Radar，SAR)原始信号处理到单视复数(Single Look Complex,SLC)成像、单视以及多视处理、基于雷达信号滤波、正射纠正和配准、DEM提取(干涉)、形变分析(差分干涉、点目标干涉分析)、土地利用等,同时也可以处理各种星载、机载及地基雷达的观测数据。

3.2 数据处理流程

3.2.1影像裁剪

影像覆盖范围广,数据量大,为了提高数据的处理效率和监测精度,根据研究区域的经纬度,裁剪出研究区域的影像数据。

3.2.2差分干涉对组合

实验选取2019年6月15日的SAR影像作为公共主影像,将辅影像配准到主影像上,使配准精度达到0.001像元。将多视系数设置为4∶1,对干涉对进行多视处理。本实验的时间基线和空间基线阈值分别为60 d和200 m,最终构成112对干涉对用于时序分析。图2为干涉对时空基线分布图。

图2 干涉对时空基线分布

3.2.3干涉图的生成及相位解缠

根据干涉对时空基线分布图对影像进行干涉处理,利用外部DEM数据去除地形相位的影响,生成差分干涉图,对差分干涉图进行滤波增强处理,采用最小费用流的方法进行相位解缠,得到研究区域相位的真实值。通过查看每组生成的相干系数图、滤波后的干涉图和解缠图,移除干涉质量较差的干涉对。

3.2.4研究区平均沉降速率的获取

在具有高相干性的像素点上建立模型,通过 SVD 法反演估算形变速率,并去除残余地形,计算相位残差,进行残差分离,最终得到研究区域内的地表沉降速率以及地表沉降分布情况。

4 地表形变监测结果与分析

采用SBAS-InSAR技术对30景升轨Sentinel-1A数据进行处理,其结果表明,在该研究区域内有三处明显沉降,分布在研究区西南部、东南部以及北部。

图3为研究区域图,如图所示,将这三个沉降区域分别命名为A、B、C。

图3 研究区域图

图4为A区域沉降速率图,从实验结果可知:研究区域西南面的A区域沉降最为严重,最大年平均沉降速率为-490 mm/a。经核查该区域采矿活动频繁,地表有多处塌陷且出现积水。为了进一步研究监测期间矿区的累积沉降量,在沉降严重的A区域选取三个特征点进行时序分析,分别为点A#1、A#2和A#3,对特征点进行时间序列分析得到时间序列形变图,通过结果可知,在2018年12月29日～2019年12月24日内,随着时间的推移,该区域形变呈现连续下降的趋势。从图5可知,位于该区域沉降中心的采样点A#3,该点的最大累积沉降量高达440 mm,该点存在连续下沉的趋势,且沉降速率较快。A#1和A#2两点选取在沉降区边缘,A#1和A#2两个特征点的沉降量分别为-184 mm、-131 mm。沉降区中心特征点的累积沉降量明显高于沉降区边缘,该区域沉降形成明显的漏斗状。沉降区域主要集中在祝家社区一带,该区域的沉降会对村民的生活造成影响,相关部门应引起高度重视。

图4 A区域沉降速率图

图5 A区域特征点时间序列形变图

沉降区域B位于研究区东南处,从图6可知其最大年平均沉降速率可达-390 mm/a,在B区域选取三个特征点B#1、B#2和B#3,对其进行时间序列分析得到特征点时间序列形变图,从图7可以看到沉降中心特征点B#3的最大累积沉降量达到358 mm,在整个研究的时间跨度之内一直处于持续下沉的状态,并且未来还有进一步下沉的趋势。B#1和B#2的累积沉降量分别为-120和-118 mm。沉降中心的沉降速率向沉降边缘逐渐降低,特征点的累积形变量从沉降中心往沉降边缘逐渐减少。

图6 B区域沉降速率图

图7 B区域特征点时间序列形变图

沉降区域C位于研究区北处,从图8可知,其最大年平均沉降速率可达233 mm/a,在C区域选取两个特征点C#1、C#2和C#3,对其进行时间序列分析得到时间序列形变图(图9),从图9可以看出,特征点C#1沉降量最大,其累积沉降

图8 C区域沉降速率图

图9 C区域特征点时间序列形变图

量达到-148 mm。由于该区域靠近水体,该地的沉降不仅因为频繁的采矿活动,也和水体周围的地质有着密不可分的关系。

5 结束语

本文采用2018—2019年间30景C波段Sentinel-1A升轨雷达数据为实验样本进行监测,结合SBAS-InSAR技术对江西德兴矿区进行地表沉降监测分析,反演出研究区域真实地表形变,进一步分析该区域沉降速率及严重区域特征点的累积形变量,得出以下结论:

(1)在2018年12月29日～ 2019年 12月24日研究期间内,研究区域有多处沉降,分布在研究区西南部、东南部、北部。西北部未见明显沉降。

(2)研究区中最为严重的沉降区域集中在矿区西南部,年平均沉降速率最大达到490 mm/a,在沉降最严重的区域中心选取特征点进行分析,该点的累积沉降量最大可达440 mm,该特征点未来还有进一步下沉的趋势。

(3)祝家社区位于研究区德兴铜矿内,矿区的过度开采严重影响了当地村民的生活,导致该区域地面塌陷、农田减少、地下水资源污染,严重影响当地居民的生命财产安全。今后需把该区域作为重点监测区域并对其进行持续监测,其应有效合理制定开采计划,防止矿区沉降问题失控、周围居民生活得不到保障等情况出现。

综上所述,通过SBAS-InSAR技术可以对矿区开采过程中产生的地表沉降进行有效监测,可以快速、精准地捕捉兴趣区内沉降区域的具体位置和大致范围,结果可知其沉降范围与矿区的开采过程大致相同。对沉降结果分析有助于防治此类灾害的发生,若持续开采可能会造成各类地质灾害问题。在矿山持续开采的过程中,针对地表沉降严重而易引发地质灾害等相关问题,对沉降数据进行准确预测就显得尤为重要,如何对沉降的时间序列数据进行有效精准预测则是下一步研究的重点。