基于SABS-InSAR的保德矿区地表形变监测与分析

樊昱初, 吕义清*, 杨娜

(1.太原理工大学矿业工程学院, 太原 030024; 2.航天宏图信息技术股份有限公司, 北京 100094)

近年来由于煤炭能源的需求量逐渐增大,特别是资源整合后,矿山生产能力得到进一步的提高,保德矿区内形成了大量的采空区。在地下煤炭资源开采的过程中,引起了采场周围岩体的应力重新分布[1],开采工作面及其周边区域失去原始应力平衡。随着采空区范围的扩大,岩层移动逐步向上传递,最终导致地表发生形变,局部出现地面塌陷和地裂缝。破裂煤岩体在应力、地下水等多种因素影响下,其强度、承载力会有一定程度的降低,加大了地质灾害发生的风险。为了保障矿区周边居民的人身财产安全和可持续开采,对于矿区进行沉降监测分析具有重要的研究意义。利用技术手段掌握矿区地表形变信息,对潜在的灾害点进行识别,为后续的防治工作提供科学的指导。

传统的观测方法有水准测量、三角高程测量、全站仪测量、静态全球导航卫星系统(global navigation satellite system,GNSS)和动态全球导航系统等[2],这类型散点式观测手段虽然有较高的精确度,但无法对整个矿区进行大面积,长时间的监测。合成孔径雷达干涉测量(interferometry synthetic aperture radar,InSAR)作为一种新兴的高精度的测量技术,能够很好的解决宏观监测的问题,其中合成孔径雷达差分干涉测量技术(differential interferometry synthetic aperture radar,D-InSAR)可以全天时、全天候、近实时获取大面积地表地形信息,且空间分辨率高,受大气和季节影响较小。王志勇等[3]利用D-InSAR技术对济宁矿区的沉降做了精细化的监测和分析,证明D-InSAR技术非常适合矿区大型变的沉降监测,监测精度可以达到厘米级。但是D-InSAR技术无法获取研究区域长时间的时序变化结果,而且易受到时空失相干和大气延迟效应影响,为了优化这一领域受限的问题,Berardino等[4]在2002年提出了小基线集技术(small baseline subset,SBAS)该技术不仅能够有效的削弱时空失相干和大气延迟效应对监测结果的影响,还能以更好的时空连续性来获取相干点目标信息,使监测结果更加精确和直观。朱建军等[5]详细介绍了InSAR技术目前的发展和主要应用领域,不仅可以对城市的沉降、矿区地表、地震、基础设施等进行形变的监测,还能监测火山,冰川等的活动,同时对滑坡灾害也有一定的识别和预警作用;刘志敏等[6]利用SBAS-InSAR技术对长治矿区的地表形变进行了监测和分析,并揭示了老采空区的残余沉降趋势;麻学飞等[7]先通过D-InSAR技术识别临汾矿区沉陷区位置,后利用SBAS-InSAR技术进行时序监测,共获取105处沉陷,为后续矿区沉陷灾害防治提供了重要依据。

SBAS-InSAR在地表形变监测[8-9]以及地质灾害预警[10]等方面已经有较为成熟的应用体系,但在矿区地表形变监测中,利用InSAR技术对形变驱动性因素的研究尚不充分,因此在此基础上采用SBAS技术对研究区32景Sentinel-1A影像数据进行处理,获取了保德矿区2年时间内的地表形变信息,结合野外勘查与相关资料,选取矿区内地质构造、地层岩性、采矿工程、降雨等因素,定性分析其对研究区内地表形变程度及范围的影响,明确形变诱因,为煤矿安全生产风险性评价提供具体的评价指标,同时识别了潜在的地质灾害风险点,为进一步的开采以及地质灾害的防治提供了科学依据。

1 研究区概况与数据来源

1.1 研究区概况

研究区域所在的保德县位于山西省西北部、地理范围为38°39′00″N～39°6′56″N、110°56′30″E～111°19′40″E(图1)。整体地形东高西低,北高南低,平均海拔为840 m,区域内地表主要被黄土覆盖,水土流失、地形切割较为严重。土壤结构松散,沟壑纵横,梁峁连绵,属于典型的黄土丘陵沟壑地形地貌。大地构造位置处于鄂尔多斯东部大型聚煤盆地东缘河曲—保德断隆以北,山西断隆以西[11]。区内主要发育近南北走向的断裂,整体构造比较简单。

图1 研究区位置Fig.1 Location of the study area

1.2 数据来源

选用覆盖保德区2020-06-24—2022-09-13期间32景Sentinel-1A升轨单视复数(single look complex,SLC)影像作为基础数据,其中path=11,frame=121有24景,path=113,frame=121有8景。卫星对地观测模式为干涉宽幅模式(IW),极化方式为VV,重访周期为12 d。获取影像均为C波段,距离向分辨率为5 m,方位向分辨率为20 m。采用欧空局发布的精密定轨星历数据(precise orbit ephe-merides,POD)对轨道信息进行修正,引入空间分辨率30 m的SRTM数据作为数字高程模型(digital elevation model,DEM)去除地形误差相位。同时获取的还有2020-06—2022-09期间保德县境内五个气象监测点的降水数据,研究区地质数据,以及重点监测区内各煤矿采掘工程信息。

2 数据处理方法

2.1 SBAS-InSAR技术原理

SBAS-InSAR是一种基于多主影像的InSAR时间序列技术,基于短基线原则,将所研究的SAR数据转换组合为有多个主影像的干涉子集,并选择其中一景影像作为公共主影像进行配准,利用时空基线较短的干涉对来提取地表形变信息。将在t1～tm时间内获取的m幅SAR影像进行干涉组合,并对集合间影像进行差分干涉处理,形成N幅干涉条纹图,满足的条件为

(1)

对tA和tB两个时间点的影像生成的第i幅干涉图中,第x个像素的干涉相位表达式为

φint,x,i=φdef,x,i+φtopo,x,i+φaps,x,i+φorb,x,i+φnoise,x,i

(2)

式(2)中:φdef,x,i为斜距向变形即形变相位误差;φtopo,x,i为地形相位误差;φabs,x,i为大气延迟相位引起的误差;φorb,x,i为基线轨道引起的轨道相位误差;φnoise,x,i为噪声误差。各误差相位可在时间序列上采用滤波方法或多项式模型予以去除。

由最小二乘法求解得到监测区地表时间序列形变信息,假设不同干涉图之间的变形速率vk,k+1,在tA～tB时间内的累积形变表达式为

(3)

将得到的N幅干涉图进行解缠处理,可以进一步求出不同SAR影像获取时间的形变速率,以上提及处理方式的具体计算方法可以参考文献[4]。

2.2 数据处理流程

使用ENVI/SARscape软件对数据进行计算处理,利用Arcgis[12]软件对结果进一步分析,具体流程如下。

(1)为了提高影像数据的处理效率,首先根据矿区矢量范围裁剪影像,通过设置时空基线阈值和多视比例值(距离向为5,方位向为1),生成连接图和干涉对组合,选取2021-07-01影像为超级主影像,其他影像依次与主影像配对。

(2)干涉工作流处理,通过导入外部DEM来去除干涉图中的地形相位,选用Goldstein方法对得到的差分干涉图进行滤波处理,采用基于Delaunay三角网的最小费用流法(minimum cost flow,MCF)对差分干涉图进行相位解缠,减小大气延迟误差、基线轨道误差、相干斑噪声误差等对计算结果的影响,提高干涉条纹的清晰度。

(3)轨道精炼和重去平,在地形平坦远离主形变区,或已知形变的区域设置36个控制点,通过这一步估算和去除残余的恒定相位和解缠后还存在得相位坡道。

(4)采用奇异值分解(singular value decomposition, SVD)法估算平均形变速率和残余地形相位,并且对合成的干涉图去平,重新进行相位解缠和精炼,生成更加优化的结果。通过进一步反演,进行定制的大气滤波,估算去除大气相位,得到更精确的时间序列位移结果。

(5)将生成的结果进行地理编码,得到研究区雷达视线向(line of sight,LOS)的平均形变速率和累积形变量。

3 时序监测结果与分析

3.1 矿区地表形变结果分析

3.1.1 形变空间分布特征

在2020-06-24—2022-09-13监测时间段内,通过SBAS-InSAR技术处理得到研究区的年平均形变速率(图2)和累积形变量(图3),其中负值表示下沉,正值表示抬升。

图2 保德矿区地表年平均形变速率Fig.2 Average annual surface deformation rate in Baode mining area

图3 保德矿区地表累积形变Fig.3 Cumulative surface deformation in Baode mining area

结合相关资料分析,研究区主要变形区域与区内煤矿工作面位置空间分布一致,K1～K6重点监测点位于各沉降区域中心位置,监测时间段内中心点K5最大的沉降速率为-43.7 mm/a,最大累积形变量为-84.1 mm。从累积形变图可以看出,研究区整体形变趋势以下沉为主,采矿工程影响范围内区域的平均沉降速率为-20.2～-5.5 mm/a,沉降量为-39～-12 mm;累积形变量为-23.1～-4.7 mm的区域在研究区内分布范围最广,该区域内分布有多个自然村落,地面沉降主要受各种人类活动影响;研究区部分边缘地区在监测期内有较为明显的抬升,平均形变速率为0.4～21.3 mm/a,最大累积抬升为38.1 mm。

3.1.2 矿区形变时空演变特征

以2020-08-23影像为基准,共选取9景影像数据来研究矿区地表形变的空间范围演变特征(图4)。

图4 保德矿区地表形变时间序列图Fig.4 Time series map of surface deformation in Baode mining area

以图1重点研究区位置为基本划分,各重点研究区内的特征点即重点监测点的累积形变曲线(图5)为参考,分析形变时间序列图,随着时间的推移,矿区形成多处由中心向四周发展的漏斗状下沉区域,整体沉降分布面积先增加后减少。在监测时间段内,形变量为2.0～38.1 mm的区域范围有显著增加,主要分布在监测区域边缘;沉降主体外围形变量为-4.7～2.0 mm的区域处于较为稳定的状态,范围基本保持不变;形变量为-23.1～-4.7 mm的区域范围在前期有一定程度的增加,在2021-02-19～2021-09-11期间持续减小后逐渐趋于稳定;形变量为-84.1～-23.1 mm的沉降中心区域在监测时间段内面积有明显的增加;随着采矿工作面的推进,采掘工程对地表的影响范围也进一步扩大,形变量为-39～-23.1 mm的区域面积也在持续增加。通过分析时间序列图可以得出在未来一段时间内,各重点区域内沉降会持续进行,主形变区域的沉降范围发展和沉降速率将会进一步加快。

图5 重点监测点形变曲线图Fig.5 Deformation curves of key monitoring points

从重点监测点K1～K6曲线(图3)整体来看,累积形变量与时间基本呈线性关系,各点形变均以下沉趋势为主,由于所处区域煤矿开采工程工况区别较大,所以各点的累积形变量大小和沉降速率都有明显的差距,其中K1点累积沉降量为-38.36 mm,K4点累积沉降量为-81.13 mm,最大形变量绝对值差为42.77 mm。

3.1.3 重点研究区形变分析

选取煤矿分布集中,形变程度较大的6个重点研究区(图1)进行详细分析,将研究区年平均形变速率图与Google卫星影像相结合(图6)。其中有部分关停煤矿以及个人经营煤矿由于资料缺失,不计入分析。

图6 重点研究区年平均形变速率图Fig.6 The average annual deformation rate in key research areas

研究区内主要涉及12个产能为百万吨级以上的煤矿,其中1号研究区主体为南茆露天煤矿,其年平均形变速率最小,为17.64 mm/a,区域内累积最大形变量为-68.93 mm;5号研究区内年平均形变速率为-27.36 mm/a,沉降程度最为严重,区域内累积最大形变量为-84.87 mm,形变主体为泰安煤矿和孙家沟煤矿。通过分析研究区的地表形变时间序列图和重点区域的形变速率图可以得出,煤矿的采掘工程和相应的人类工程活动是矿区出现局部沉降,形成沉降中心的主要因素。

3.2 研究区地表形变驱动因素分析

为了进一步研究控制和影响保德矿区地表形变的因素,下文将从矿区的地层岩性和地质构造,煤矿的开采工程以及研究区内降水量这4个方面进行具体的分析。

3.2.1 开采工程

以5号研究区中的泰安煤矿生产井田为例(图7),井田内含煤地层主要为石炭系上统太原组(C3t)和二叠系下统山西组(P1s),其中有5层全区可采煤层,自上而下为山西组的8号煤层和太原组的11、12、13、15号煤层,目前煤矿主要开采8、11、12号煤层。

图7 泰安煤矿采掘煤层工作面位置Fig.7 Location of coal seam working face in Tai’an Coal Mine

(1)8号煤层。该煤层井田内西南部大面积采空,煤层厚度为3.31～6.55 m,平均4.25 m,煤层赋存标高为1 100～720 m,其中8101工作面于2013年采空,8102工作面于2014年采空,8103、8104工作面于2015年采空,8105、8106工作面于2016年采空,8107工作面于2017年采空。

(2)11号煤层。该煤层在井田中南部部分采空,上距8号煤层23.9～44.48 m,平均为35.72 m,煤层厚度为1.30～1.95 m,平均为1.53 m,煤层赋存标高为1 060～680 m,其中11101工作面于2017年采空,11102工作面于2018年采空,11103工作面Ⅴ区在2019年回采347 m,Ⅵ区在2018年回采588 m。

(3)12号煤层。该煤层井田内东部大面积采空,上距11号煤层2.43～7.84 m,平均为3.63 m,煤层厚度1.00～6.50 m,平均为3.21 m,煤层赋存标高1 050～690 m,其中东部大面积采空区于2012年采空,12101工作面Ⅶ区在2019年采空,Ⅷ区在2018年采空,12102工作面Ⅴ区在2020年采空,Ⅵ区在2019年采空。

将各煤层采矿工作面与图6中泰安煤矿矿区内地表年平均形变速率图对比可以看出,地表发生变形的区域与采掘工作面空间分布相对应,其中11号煤层中的11104、11105工作面以及12号煤层中12102～12106工作面均在2020年后有采掘工程,2020-06—2022-09,泰安矿区内局部最大累积形变量达为-84.87 mm,12103、12104工作面对应空间的地表平均形变速率可达-33.25 mm/a;8号煤层中的8106及8107工作面分别于2016年和2017年采空,在监测时间内,其累积最大形变量为-30.91 mm,对应区域的平均形变速率为-7.29 mm/a;12号煤层东部采空区于2012年已完成开采,但在监测时间内,仍能观测到较为明显的地面沉降,其中最大累计形变量为-23.11 mm,采空区域内平均形变速率为-4.35 mm/a。根据相关文献资料,采空区在形成十年后,残余沉降基本消失,完成总变形量90%以上[13-14],地表一般已经进入相对稳定的状态,造成东部采空区异常沉降的原因推测是因为采空区与12号煤层活跃开采工作面距离较近,且受到地下水抽取以及地下矿井巷道工作影响,导致其仍存在一定程度的下沉。通过对比正在开采、采空三到四年以及开采完成十年以上不同工况工作面对应区域的形变速率可以得出:人为的采矿工程是造成和加剧矿区地表形变的主要原因。

3.2.2 地层岩性与地质构造

地层岩性与地质构造是地面沉降发生的重要地质背景,研究区内断层构造的分布以及地面可压缩层特性等因素对沉降的发展都有一定程度的控制[15-16]。因此研究区域构造背景是分析矿区地表形变驱动性因素的重要一环。保德地区位于鄂尔多斯盆地东缘的北部,以升降运动为主体,整体表现为面状隆起与局部相对下沉。

为了进一步研究断层构造对矿区地表形变的影响,将构造图(图8)与矿区地表形变速率图(图9)进行叠加分析,可以直观看出矿区整体范围内沉降发展趋势与断层方向一致,断层两侧的地面沉降速率受断层的上、下盘滑动及相对升降运动的影响有较为明显的差异,另一方面断层两侧地层岩性即可压缩层厚度的不同,也是导致沉降速率变化的原因之一。根据形变速率图和监测点数据计算显示,断层两侧年沉降速率的平均差异约为3.1 mm/a;保德县地处晋西北黄土高原,其黄土分布范围较广,第四系堆积物占总面积的69.98%。根据岩层的成因类型、结构和工程地质特征,区内的工程地质岩组划分为3个大类,9个亚类(表1)。

图8 保德矿区地质构造要图Fig.8 Geological structure map of Baode mining area

图9 保德矿区地表形变与断层分布Fig.9 Surface deformation and fault distribution in Baode mining area

研究区内砂土由第四系全新统(Q4)冲、洪积物组成,岩性为砂卵砾石和砂层,属于低压缩性土。黏性土主要为新近系保德组红色黏土(N2),底部为砾岩,压缩系数为0.02～0.082 MPa-1。特殊性岩土则按照形成时间分为新黄土和老黄土,据土工试验,老黄土颗粒成份中黏粒含量高于新黄土,其压缩系数为0.005～0.132 MPa-1,新黄土则为0.002～0.076 MPa-1。相关研究表明,在相同开采工况下,以粉土、黏性土层为主的区域,单位沉降量是以砂层、砂砾石层为主要岩性区域的2～3倍[17],以碳酸盐岩带为基岩的区域的平均形变速率要快于非碳酸盐岩带区域。在研究区内第四系堆积物分布广泛的地质背景下,相关的人类地面活动,地下工程的开发都会导致软土层的压实。此外,人工开采引起的周边区域内地下水水位持续波动,不仅会导致含水砂层压实还会进一步潜蚀碳酸盐岩[18],这也是开采区附近区域沉降速率与其余地区有较大差异的重要原因之一。

3.2.3 降雨影响

据保德县气象局1980—2020年共计40年统计资料显示,区域多年平均降水量为475.9 mm,且年际变化较大,丰水年平均降水量为861.4 m,枯水年平均降水量为194.8 m,其中70%的降水量集中在7、8、9三个月,并且多以暴雨出现。为了研究降雨量对矿区地表形变的影响,向相关气象部门收集了研究区域内自2020年6月起5个地面气象监测点逐月降水资料,涉及义门镇、桥头镇、孙家沟乡、南河沟乡四个区域,由于乡镇之间的降水量有一定的差异,因此分别进行了单独统计,将重点研究区内监测点(图6)的地面累积形变量与监测时间内区域月降水量进行叠加分析(图10)。

图10 保德矿区降水量与累积地表形变Fig.10 Precipitation and cumulative surface deformation in Baode mining area

2～6号重点监测区内,虽然各煤矿采矿的工况有所不同,但从整体形变量走势分析,降水量的增加会减缓地表形变的发展。矿区地表主要被黄土覆盖,黄土由粉土、粉质黏土组成,透水性一般较差,因此在降雨强度较小的4、5、10月时间段内,一定量的雨水渗入土壤,浅层土壤中水分含量会有显著增加[19],此时体现在地表形变图中,监测点下沉趋势减缓,部分测点由于土壤吸水膨胀,地表会发生小幅度抬升。而在此后的11月到来年2月期间,因水分蒸发等原因,表层土壤的含水量会进一步下降,此时地表会发生季节性沉降。而在7、8、9月,由于多暴雨且降雨时间集中,降水会在短时间内汇集,形成的地表水流会在黄土构造节理、风化裂隙、陷穴等发育部位的空隙下渗或灌入,在相对隔水的部位形成上层滞水或饱水带,增大岩土体重力[20],体现在形变图中,强降雨会加速地表下沉进程。图6中C2、B4、A5、C6测点的形变状态对降雨量大小的响应相对及时。以B4测点为例,该点位于王家岭煤矿2019年采空区附近,因此采矿工程扰动对地表形变影响较小,在2021-05—2021-07降雨期间,其地表累积抬升7.32 mm,在2021-11—2022-02期间,地表累积下沉21.07 mm,B4测点在监测时间内地表形变趋势以下降为主,但从折线图可以看出降雨对该点的形变发展起到了一定的控制所用。

1号研究区主体为南茆露天煤矿,A1、B1测点分别在两个黄土开挖边坡工作面附近,C1测点位于煤矿排土场东侧,排土场是露天煤矿开采过程中产生的废弃土石的堆放场所,成份主要由黄土、泥岩、砂质泥岩、粉砂质泥岩等组成,其中还混有前期开采排弃的碎石,砂质土[21-22]。虽然排土场经过一定程度的压实处理,但结构仍然为松散多孔,整体强度较低,在丰沛降雨时间段内,大量雨水会入渗排土场中,导致黄土等岩土体的抗剪强度降低,密度变大,内部的力学性质发生改变,致使其整体有一定量的下沉。而在集中暴雨情况下,降水会在排土场的表面和坡面形成径流,在雨水的侵蚀和冲刷作用下,排土场局部会发生塌陷,在极端降雨条件下,排土场边坡可能会发生失稳,滑坡等灾害。从C1测点的形变折线图分析,在2021-08—2021-09时间内,累积下沉达11.35 mm,2022-06—2022-09期间,累积下沉17.78 mm,总体来看,降雨加速了该点的沉降速率。A1、B1测点整体形变发展受降雨影响较小,但在少量降雨的月份,地表会发生小幅度抬升,无雨的时间段内也发生了一定的季节性沉降。需要说明的是,由于露天开采矿区部分区域有较大程度的形变(图11),因此在数据处理过程中会出现形变失相干的现象,导致这部分区域形变信息的缺失,如大规模岩土体的转移,局部发生严重的塌陷,滑坡等,但通过分析已知测点的形变发展,仍然可以得出降雨量作为诱发地面沉降的重要自然因素之一,与矿区的地表形变有很大程度的相关性。

图11 1号重点监测区卫星影像Fig.11 Satellite image of No. 1 key monitoring area

4 结论

利用SBAS技术对覆盖研究区的32景Sentinel-1A影像数据进行处理,获取保德矿区2020—2022期间形变信息,分析了矿区的形变特征以及引起形变的驱动性因素,得出如下主要结论。

(1)研究区存在6个较大的形变区域,区内地表形变以下沉为主且局部形成多处由中心向四周发展的漏斗状沉降中心,主形变区域监测点累积形变量与时间基本呈线性关系。

(2)监测时间段内,区域的平均地表形变速率分布在-44～24 mm/a,最大累积沉降为85 mm,最大累积抬升为38 mm。

(3)采矿工程是研究区发生地表形变的主要驱动因素,沉降中心空间分布范围与各煤矿工作面位置有良好的一致性,断层分布和地表可压缩层厚度差异对研究区整体的形变起到了一定的控制作用。从监测点累积形变量整体来看,大多数情况下降水量的增加会减缓矿区地表形变的发展,而在极端降雨条件和某些特殊区域中,降雨则会加快地表形变速率,甚至会引发崩塌,滑坡等地质灾害。

(4)通过对形变信息的分析,明确研究区内地表形变诱因,在后续工作中以上述驱动因素为评价指标,结合研究区内人口密度、建筑密度、道路密度、人类工程活动等数据对矿区地质灾害发生的风险性以及易损性做出合理的评价,为煤矿安全生产和形变区的综合治理提供了科学依据,但SBAS-InSAR方法对矿区最大形变量探测能力有限,需要结合其他勘测方法才能更加准确、全面地反映矿区实际情况。