基于SBAS-InSAR时序分析技术的甘肃省红会矿区地面沉降监测及其灾害发育特征研究*

王 颂 张路青 周 剑 张奋翔 韩振华 尹凡龙

(①中国科学院地质与地球物理研究所,中国科学院页岩气与地质工程重点实验室,北京 100029,中国)(②中国科学院大学,地球与行星科学学院,北京 100049,中国)(③中国科学院地球科学研究院,北京 100029,中国)(④北京工业大学,城市与工程安全减灾教育部重点实验室,北京 100124,中国)(⑤北京市勘察设计研究院有限公司,北京 100038,中国)(⑥北京中科地环工程技术有限公司,北京 101499,中国)

0 引 言

人类矿山工程与地质环境之间矛盾尖锐,大规模的采矿活动改变了区域能量交换和物质转移过程,破坏了局部地质体的稳定性,继而容易诱发地质灾害(胡静云等,2014; 周超等,2020)。矿山地质灾害类型众多,既包括发生于井下的冒顶、片帮、岩爆、突水、突泥等,也包括发生于地表的地面沉降、地面塌陷、地裂缝、崩塌、滑坡、泥石流等(宁宇,2009; 郭维君等,2010)。其中:地面沉降、地面塌陷和地裂缝是地下矿体开采后上覆岩层冒落,围岩因内部原有应力平衡状态被改变而不断移动变形,最终造成地表下沉破坏的灾害现象(何芳等,2008,2010; 黄平路等,2010)。高强度矿山开采所诱发的地面沉陷不仅导致了地表水和地下水循环紊乱、植被损伤、矿区沙漠化等一系列生态环境问题,还毁坏了采空区地表的房屋和农田,易使社会矛盾激化(王卓理等,2011)。

甘肃省红会矿区是受高强度地下开采影响而形成的典型地表沉陷区,地表塌陷坑和地裂缝错综复杂,危害显著。目前已有相关文献针对矿区地面沉降进行了初步探索。李彦军(2008)重点论述了矿区地裂缝和塌陷坑的发育特征、成因及其所造成的危害。然而,受当时科学技术条件所限,他所获取的灾害数量和空间分布数据十分有限。随着遥感技术的飞速发展,大范围内的地貌和地物得以快速便捷观测和精确解译。刘欢等(2010)、丁丽等(2010)和李莹莹等(2014)使用高分辨率的IKONOS卫星影像识别并统计了红会矿区内的塌陷群和煤堆、煤矸石堆等土地占用情况,但没有涉及到地裂缝的研究。王娟等(2012)和尹展等(2017)基于RS和GIS对红会矿区地质环境进行了评价和分区,但评价过程中没有考虑体现地面沉降动态过程的相关指标。因此,目前针对红会矿区地面沉陷的研究还很不充分,特别是缺乏对灾害发育特征的定量描述和矿区地表动态变形的监测。

广泛使用的地表变形监测技术有精密水准测量、数字摄影测量、GPS实时监测、分布式光纤监测等(卫建东,2007; 董文文等,2016)。合成孔径雷达干涉测量(InSAR)是一种基于主动式微波遥感区分地表属性和记录地表三维形变信息的新型监测手段。和传统的变形监测技术相比,InSAR具有全天候、大面积、高精度、高分辨率、低成本等突出的优点,适用于火山活动、滑坡活动、冰川运动和地面沉降等诸多研究领域(李珊珊等,2013; 杨泽发等,2015; Rosi et al.,2018; 兰恒星等,2019; 姚佳明等,2020)。目前,应用最多的InSAR变形监测方法包括D-InSAR、PS-InSAR、SBAS-InSAR、IPTA-InSAR、DS-InSAR和MAI技术等(Tao et al.,2014; 朱建军等,2017;Xie et al.,2020; 张严等,2021)。例如,Chen et al.(2012)基于IPTA-InSAR技术研究了2006～2010年中国珠江三角洲地区地质、人类活动与地面沉降的相互作用,发现和第四纪沉积活动相比,人类城市化活动更容易诱发地面沉降。Zhang et al.(2018)使用D-InSAR技术提取了特拉华盆地钾盐矿地面沉降信息,发现地面沉降和钾盐矿的开采强度有密切的相关性,并且开采停止后地面仍然在持续下沉。

上述各种方法中,SBAS-InSAR选取空间、时间基线均较小的干涉像对,能有效缓解数据处理过程中的失相干情况,提高了形变监测精度,更适合应用于城市及矿区地面沉降监测(李珊珊等,2013)。张金芝(2015)基于SBAS时序分析技术提取了1992～2000年黄河三角洲地面沉降速度,并重点分析了石油开采、地下水抽取和地表加载等典型影响因素对地面形变的作用。Zhao et al.(2011)监测到山西大同新城建设区、地下水开采区和采煤塌陷区的地面沉降速度分别是1cm/年、2cm/年和4cm/年。Zhang et al.(2019)将SBAS-InSAR形变分析结果与武汉市56个特定水准点测量结果进行对比,结果表明,SBAS-InSAR技术测量精度可达毫米级; 武汉市分布有4个主要沉降区,地下水开采、地铁开挖和城市建设是造成城市沉降的主要因素。

综上所述,本文以甘肃省白银市平川区红会矿区为研究对象,结合遥感解译和现场调查定量分析地面变形灾害的发育特征,并利用SBAS-InSAR技术系统研究该区域从2017年3月到2019年10月期间地面变形程度和随时间的演化规律,为地面灾害的预警和治理提供依据。

1 地质概况

1.1 区域地质背景

白银市平川区位于甘肃省的中北部,区域海拔为1347～2858m,整体地形呈阶梯状,东南高、西北低,由东南向西北倾斜,多侵蚀-构造中山、低中山。平川区的东北和东南分别有六盘山和秦岭作为屏障阻隔暖湿气流,故区内气候干燥,年平均降雨量不足200mm。北近腾格里沙漠,西伯利亚寒流易入侵,因此冬季寒冷,年平均气温约8.0℃(尹展,2012)。

大地构造上,平川区地处北祁连构造带东段,秦祁褶皱强烈,经历了加里东期、燕山期和喜山期多次大型构造运动,表现为强烈的沉降和局部的隆升。区内NNW向的断裂构造发育,NNE向构造次之,主要有车轮口逆断裂、响泉口逆断裂、大水头逆断裂等(刘欢,2011)。前寒武系至第四系地层在区内均有出露,其中:第四系覆盖面积最大,岩性主要为砂土和亚黏土; 白垩系地层岩性为紫红色砂岩; 侏罗纪陆相含煤沉积形成该区域的含煤地层,可分为延安组和富县组上、下两部分,成矿条件好,煤炭资源总量位居甘肃省第二; 区内还出露小面积岩浆岩,主要是志留纪花岗闪长岩和斜长花岗岩(李莹莹等,2014)。

1.2 研究区工程地质条件

本文研究区主要为白银市平川区中北部的红会矿区及其周边地区(图1)。矿区中心经纬度为36°40′32″ N,105°04′47″ E,高程为1902m。煤层赋存以厚煤层为主,埋藏深度200～500m。研究区内有红会一矿、红会二矿、红会三矿、红会四矿、会通煤矿、长征煤矿、神峰煤矿等多家煤矿企业,均为地下开采。其中:红会煤矿隶属于靖远煤业集团有限责任公司,含煤面积约18km2,是甘肃省重要的煤炭基地之一。矿区内主要涉及的行政村有西合村、牛拜村、玉湾村等,南侧有多条省道和县道穿过,兼有红会铁路直达矿山,交通便利。

图1 研究区概况

研究区地处西格拉滩盆地的西部,地形相对平坦,附近的地层岩性和地质构造如图2示。该区域属于宝积山向斜的北翼部分,北西向分布的压性逆断层发育(刘欢等,2010)。整个研究区内出露的岩性主要有第四系全新统洪积层(Qhpl)的砂、砾及砂土; 中侏罗统龙凤山组(J2l)的含砾砂岩、细砾岩; 上三叠统西大沟组(T3x)的砂岩夹粉砂岩、页岩; 志留纪花岗闪长岩(Sγδ)等。由于坐落在祁连山断裂带和西海固断裂带的交会处,研究区地震活动频繁,根据中国地震动参数区划图(中华人民共和国国家质量监督检验检疫总局等,2015),区内地震烈度为Ⅷ 度,地震动加速度反应谱特征周期为0.45s,地震动峰值加速度为0.3g。

图2 红会矿区周边地质图

矿区地表仅有数条因沟谷汇聚而形成的沙河,在暴雨时可形成洪流并最终汇入黄河。矿区地下水主要为赋存于白垩系、侏罗系砂砾岩中的碎屑岩类孔隙裂隙水和赋存于山前冲洪积物中的松散岩类孔隙水。研究区内主要人类活动是矿产资源开采,农耕等。受地形条件限制,区内崩塌、滑坡、泥石流等山地灾害不发育。但由于大规模高强度的采煤活动,地面塌陷和地裂缝极其发育。

2 灾害发育特征分析

红会矿区历史开采形势复杂,尤其是20世纪中后期有近百家地方及私营企业同时进行近地表煤炭开采,不仅破坏了煤炭赋存的整体性,造成大量资源的浪费,而且使用落后无序的开采方式采煤后在井下形成了大范围不规则采空区,诱发了地面沉降、塌陷坑和地裂缝等地质灾害。尽管进入21世纪后,开采技术大幅提升,规模化、机械化开采迅速推广,但能源需求的增加迫使煤矿开采强度不断加大,使得矿山地面地质灾害加剧。

结合遥感影像,解译和统计了截至2019年红会矿区地面塌陷坑和地裂缝的空间分布,尺寸和数量,如图3和图4所示。在红会矿区地表共发现地裂缝735条,塌陷坑559个,地裂缝平均长度约200m,塌陷坑平均面积约360m2。其中:地裂缝长度主要分布在100～300m范围内,塌陷坑面积主要分布在<400m2的范围内。从图3a中还可以看出,地裂缝和塌陷坑密集分布在各个煤矿周围,呈现出明显的方向性,大致沿北东向分布,与该区域煤层的分布基本一致。

图3 红会矿区塌陷坑和地裂缝空间分布

图4 红会矿区塌陷坑和地裂缝数量及尺寸分布特征

此外,由于红会一矿长年运用走向长壁综采放顶煤采煤方法进行地下采煤活动,较薄的碳酸岩顶板大面积垮落,致使平川区黄峤乡牛拜村附近成为受地面塌陷影响最为严重的区域。调查发现,牛拜村地面沉陷区长约2.6km,宽约1.5km,影响面积约3.9km2,且随着该区域地下开采的推进,矿山地面灾害数量呈逐年递增的趋势。沉陷区内的部分塌陷坑深4m,半径约15m,面积可达700m2; 地裂缝纵横交错,并穿过塌陷坑(图5)。整个牛拜村周围地面沉降大、塌陷群分布广,造成周围村镇的房屋、耕地、道路破坏严重,牛拜村村民的生命和财产安全无法得到保障(图6)。

图5 牛拜村地面沉陷区内的塌陷坑与地裂缝照片

图6 牛拜村地面塌陷产生的主要危害

针对红会矿区一系列的矿山地质环境问题,白银市常年开展环境恢复治理工程,其中:2012年和2013年开展的两期重点治理工程(图3b),主要包括红会一矿周围区域煤矸石堆、废煤堆平整,挖土填方,土地复垦,种植果树等,并且对周围房屋损害严重的部分居民进行搬迁。通过对比治理工程前后的遥感影像并统计相关数据,发现红会矿区内的地裂缝减少183个,塌陷坑减少108个,工程效果显著(图3b、图3c和图4)。然而,治理工程的范围很小,约为整个沉陷区的十分之一。尽管平川区人民政府在2019年的《白银市平川区采煤沉陷区综合治理行动实施方案》中指出,计划在未来5年实施兄弟村、中和村、西合村群众1017户4058人和牛拜村群众95户528人的搬迁安置任务,并继续推进采煤沉陷区的生态修复和重大地质灾害隐患治理(平川区人民政府,2019); 但随着采矿活动的继续进行,矿山地面地质灾害仍有进一步发展的趋势,因此,十分有必要对该区域的地表形变进行跟踪监测,以支撑灾害的预警和治理。

3 小基线集(SBAS)技术

3.1 基本原理

在SBAS技术中通过使用每次D-InSAR计算的形变结果即可获取高精度的形变时间序列(Berardino et al.,2002; 李凌婧等,2014)。首先,选取N+1景覆盖同一区域、时间分布为t0,t1,…,tN的SAR数据影像图,进行基线估算后,选取其中1景作为超级主影像,其余N景为辅影像,将时空基线均小于某一阈值的影像进行配对,共形成M对干涉像对,其中M满足:

(1)

将M对干涉像对采用双通道D-InSAR进行差分干涉处理,去除平地相位和地形相位后,假设在tA和tB时刻获取的两景SAR影像生成第i幅解缠的多视差分干涉图,那么SAR坐标系下任意坐标(x,r)处相干点的相位值可以表示为(Guo et al.,2017):

Δφi(x,r)=φ(tB,x,r)-φ(tA,x,r)

≈Δφi,def(x,r)+Δφi,topo(x,r)+

Δφi,atm(x,r)+Δφi,θ(x,r)+

Δφi,nosie(x,r)

datm(tA,x,r)]+Δn

(2)

式中:λ为雷达波长; Δφi,def(x,r)是累积形变引起的相位; d(tA,x,r)和d(tB,x,r)分别是tA和tB时刻相对于参考时间t0时刻在雷达视线方向(LOS)上的累积形变相位量; Δφi,topo(x,r)是由地形误差所导致的相位,与垂直基线B⊥和高程误差Δz成正比,与入射角θ的正弦成反比; Δφi,atm(x,r)是两幅影像中大气延迟引起的相位datm(tA,x,r)和datm(tB,x,r)所造成的相位误差;Δφi,noise(x,r)为噪声残余相位; Δn表示失相关引起的噪声相位。其中:大气相位误差和噪声误差可以通过时域高通滤波和空域低通滤波与干涉相位分离。

在去除上述各种相位误差分量后,由式(2)可以得到M个方程,方程中含有N个未知量,用矩阵表示为:

Aφ=Δφ

(3)

式中:A对应一个M×N的系数矩阵,每一行和每一列分别对应一个干涉图和一景SAR影像,因此M和N分别为干涉图和SAR影像的数量。

直接求解相位将违反形变的物理规律,故将两景SAR影像获取时间之间的平均相位速度作为未知量来进行求解(李珊珊等,2013),

(4)

式中:φ(t0,x,r)=0。

结合式(3),(4),可以得到一个新的矩阵方程:

Bv=Δφ

(5)

式中:B为M×N的系数矩阵。采用奇异值分解法(SVD)和最小二乘法进行求解方程(4),即可得到形变速率,进一步在时间域上积分得到整个时间序列的形变相位值。

3.2 数据介绍

本文地表形变监测选用由欧洲空间局(ESA)发射的Sentinel-1A卫星获取的覆盖研究区域的32景单视复数影像(SLC)作为实验数据,数据范围如图1。表1显示的是所用SAR影像的基本参数情况,其中时间跨度为2017年3月27日到2019年10月19日,影像数据的极化方式为VV,轨道方向为升轨,分辨率为5m×20m。为了提高影像的轨道数据精度,本研究还采用了欧洲空间局发布的AUX＿POEORB精密轨道参数文件进行轨道校准。此外,所用数字高程模型数据为美国宇航局(NASA)公开的SRTM数据(图7),分辨率为30m×30m,用于消除地形相位的影响。

表1 Sentinel-1A 影像技术参数

3.3 数据处理流程

利用SBAS-InSAR技术进行红会矿区地表形变监测的数据处理流程如图8所示。32景数据裁剪后依次以每一景数据作为超级主影像,进行基线估算,其中最大的绝对空间基线为254.409m,最小的绝对空间基线为1.665m,最大的绝对时间间隔为936d,最小的绝对时间间隔为24d。设置最大空间基线为临界基线的45%,最大时间基线为550d,计算各组像对的空间基线和时间基线如图9。

图8 SBAS-InSAR数据处理流程

图9 时空基线分布图

根据时空基线图设置2018-05-21的影像为超级主影像,其他影像与其配准,经差分干涉生成、轨道精炼和重去平、形变速率和残余地形估算、地形残余相位和大气相位去除、形变序列估算、地理编码等多项处理步骤后获得最终的形变结果。

4 红会矿区地表形变分析

红会矿区各时间段地表变形的累积量及变形速率如图10a。从图中可以看出,计算区域内累积沉降量较大的位置集中在划定的研究区内,2017年3月27日～2019年10月19日近两年半的时间里,研究区内累积最大沉降量可达177mm,最大平均年沉降速度达70mm·a-1。快速沉降区主要集中在红会一矿、红会二矿和红会四矿周围,沉降量从沉降中心向周围递减,共和镇、牛拜村、西合村和玉湾村均受地面沉降灾害影响。将计算范围的结果图进行裁剪得到研究区内的结果图(图10b),并与遥感解译的地裂缝和塌陷坑分布图相叠加,可见,红会矿区现有的地裂缝几乎全部分布在快速沉降区,随着煤矿采空区的进一步增大,极有可能继续扩展。多数塌陷坑则分布在缓慢沉降区,变形相对稳定,进一步灾变的可能性很小。尽管环境恢复治理区没有新的地裂缝和塌陷坑产生,但依旧在遭受采矿活动的剧烈扰动,沉降较为明显。

图10 SBAS-InSAR时序分析获取的2017/03/27～2019/10/19期间的累积沉降量和沉降速度(雷达视线方向)

以2个月为时差展示地面沉降过程,图11所示的研究区形变时间序列。以2017年3月27日为起点,假设此时地面沉降为0,则从图中可以看出,随着矿山开采活动的进行,地面沉降逐渐严重,最先呈现红色的即为沉降中心,并随着时间开始向周围发展。沉降最为严重的为红会一矿和红会二矿周围,这也符合矿山生产的实际情况。

图11 研究区形变时间序列图

为了进一步研究地面沉降速度及其对周围村镇、公路的影响,如图3a所示,选取典型的监测点(P1～P22),给出其整个时间范围内的沉降曲线和距离研究区约30km的靖远县雨量站日降雨量数据(图12)。其中:P1和P2位于2012～2013年的环境恢复治理区,解译结果表明,该区域仍以每年40～50mm的速度沉降,治理效果难以维持。P3和P9代表红会一矿和二矿矿区内沉降比较严重的位置点,在两年半的时间内,几乎以线性速度沉降,最终累积沉降量分别达130mm和120mm。P4～P7所在区域为历史塌陷严重区,该区域还在不断沉降,累积沉降量达80mm,但是其不是目前沉降的中心区域,沉降速度呈降低趋势。P8位于快速沉降区的边缘地带,沉降速度大幅降低,年均沉降速度小于10mm。结合P10～P12空间位置和时序变形数据,表明了距离开采中心越远,地面形变强度逐渐减小。P1～P12所在区域原本为农田,由于地裂缝和塌陷坑导致无法灌溉,已基本荒废。

图12 红会矿区典型监测点时序变形(雷达视线方向)

P13为红会三矿周围的位置点,沉降量很小,不到10mm。P14点位于红会四矿矿区,累积沉降量达70mm,但在2018年7月之前沉降速度较大,之后沉降缓慢,这主要因为四矿资源枯竭,生产减慢,即将面临关停。P15神峰煤矿沉降量约20mm,P16西合村距离红会四矿较近,近两年沉降量很小,不到10mm。P17、P18和P19分别为牛拜村、玉湾村和共和镇位置点,累积沉降量分别为30mm、15mm和25mm。P20、P21和P22为省道308上面的3个位置点,沉降量相差不大,为20～25mm,调查时路面未见明显裂缝,但已立警示牌告知注意路面塌陷。

观察P17～P22监测点时序变形曲线,可见6条曲线的形态大致相同,几乎都在每年7月累积沉降量达到峰值,之后沉降变形减弱,地面甚至开始抬升。考虑到降雨主要集中在每年的6、7、8月份,这说明这些人口聚集区的地面沉降变形,一定程度上还与地下水过量开采有关,雨季时雨水补给地下水,使得原本沉降的地面又有所回升。但从整体沉降趋势来看,这些区域仍处于缓慢下沉阶段。值得注意的是,本文所使用的SAR影像数据仅包含了近两年半的地表变形信息,结合遥感历史数据和现场调查发现,牛拜村、西合村、兄弟村等村庄及农田区域地面沉降严重,伴有大量的塌陷坑和地裂缝,这是煤矿数十年高强度地下开采导致的结果。因此,十分有必要尽快推进沉陷区群众搬迁安置以及生态修复、环境整治和重大地质灾害隐患治理等工作。

图A1 SBAS-InSAR时序分析获取的2017/03/15～2019/11/12期间的累积沉降量和沉降速度(雷达视线方向)

5 结 论

(1)依托遥感解译和现场调查,获取了甘肃省红会矿区地面塌陷坑和地裂缝的空间分布,尺寸和数量,结果表明,红会矿区数十年高强度的地下开采导致至少700条地裂缝和500个塌陷坑沿北东向密集分布于矿区地表,造成牛拜村、西合村、兄弟村等村庄的房屋严重受损,农田无法灌溉。

(2)基于2017年3月至2019年10月期间的32景Sentinel-1A数据,利用SBAS-InSAR技术获取了红会矿区的地面平均形变速率与累积沉降量,并针对典型位置点进行了时序变形分析。结果发现,快速沉降区集中在红会一、二、四矿周围,与地裂缝空间分布位置相吻合,累积最大沉降量达170mm; 环境恢复治理区仍受采煤活动的剧烈扰动,以年均40～50mm的速度沉降。

(3)牛拜村、西合村、兄弟村等人类聚集区的地表变形整体呈缓慢下沉趋势,累积沉降量约15～30mm; 这些区域的地面沉降除受地下采煤影响之外,还可能与地下水过量开采有关,雨季时原本沉降的地面有明显回升。

(4)整个红会矿区地面变形灾害严重且处于发展态势,应规范矿业行为,采用合理的采矿技术方法,进一步控制灾害的恶性发生,同时尽快推进沉陷区群众搬迁安置、生态修复和地质灾害隐患治理等工作。

附 录A

在未能获取同时间阶段内现场实测地表形变数据的条件下,使用SBAS-InSAR技术和位于不同轨道时的数据再次提取了研究区域内的地表变形信息,以对前文所述的监测结果进行交叉验证。其中:所采用的数据为2017年3月15日到2019年11月12日期间33景降轨的Sentinel-1A数据,影像数据的极化方式为VV,其他参数与3.2节中所述的一致。最终,获得的红会矿区各时间段地表变形的累积量及变形速率如图A1所示。

对比图A1和图10a,在两种不同类型观测数据模式下,两者的累积最大沉降量和年均沉降速度有小幅差距,这可能是因为最终获取结果为雷达视线方向形变而不是真实三维形变,故而在一定程度上受雷达观测方向的影响。但整体上,研究区域内两者的地表形变基本吻合,快速沉降区的空间分布保持一致,集中在红会一、二、四矿周围,证明了监测结果可靠性较高。