采空区—滑坡—泥石流链式灾害隐患综合遥感识别与评价

2023-02-23焦润成郭学飞杨红磊王晟宇驰韩建锋马晓雪赵丹凝

金属矿山 2023年1期

焦润成郭学飞南赟杨红磊曹颖王晟宇闫驰韩建锋马晓雪赵丹凝倪璇马驰

(1.北京市地质灾害防治研究所,北京 100120;2.中国地质大学(北京)土地科学技术学院,北京 100083)

矿产资源地下开采会形成地裂缝和塌陷坑等地质灾害。开采活动停止后,采空区山体在适当的地质地貌条件下仍会持续变形,形成采空—滑坡—泥石流链式灾害隐患[1],严重影响矿区及周边安全。因此有必要对采空区开展链式灾害隐患识别与评价,对于区域居民防范地质灾害、保障城市地质安全等具有重要意义。

山区地质灾害成因复杂,尤其对于滑坡而言,往往具有高位、隐蔽性强等特点[2-3],传统调查方法周期长、成本高。近年来,大比例尺航空摄影测量、In-SAR、机载LiDAR等技术已经被应用在采空沉陷、滑坡和泥石流等地质灾害隐患识别、监测领域中[4-11]。链式灾害破坏性强、危害大,近年来也有不少学者开展了相关研究。张永双等[12]开展了汶川地区灾害链形成演化过程研究,将地震—滑坡—泥石流灾害链形成、演化过程划分为4个阶段,并提出了高位泥石流的判识指标;种艳等[13]利用“光学遥感回溯分析—地球物理快速探测—数值模型定量评价”方法,在舟曲县立节镇开展了高山峡谷区滑坡—泥石流灾害链成灾模式与危险性评价;曾庆利等[14]采用地面调查、访问和遥感解译方法,分析了新疆叶城“7·6”滑坡泥石流灾害形成条件、灾害链过程与致灾机理;张鹏等[15]采用SINMAP斜坡稳定性分析模型和聚焦模型,开展了甘肃南部典型小流域泥石流灾变链预测研究。上述成果对于灾害识别与防控取得了一定的效果,然而现阶段对于灾害链的研究大都基于已发生的灾害开展回溯分析,对于链式灾害隐患的识别和评估相对欠缺。本研究以北京西山东江沟为例,基于大比例尺航空摄影测量、InSAR、机载LiDAR等综合遥感技术,开展了采空区—滑坡—泥石流链式灾害隐患识别和易发性评价,为采空区灾害防治提供技术依据。

1 研究区概况与数据源

1.1 研究区概况

研究区位于北京市房山区史家营乡东北部的东江沟,在地质构造位置上位于百花山向斜南东翼,地层倾向NW,倾角10°～50°,断裂构造发育,总体呈NE走向,倾角为50°～80°(图1)。区内煤炭资源丰富,主要分布在侏罗系下统窑坡组地层单元内。区内岩石类型以页岩、粉砂岩、煤层为主,岩质软,煤矿开采活动诱发了诸多地质灾害及隐患。根据北京市2022年6月公布的地质灾害隐患点台账,东江沟内发育有采空塌陷1处,小型崩塌1处,小型不稳定斜坡1处(图1),同时东江沟本身也存在泥石流隐患。

图1 研究区地质特征Fig.1 Geological characteristics of the study area

1.2 数据源

本试验选取研究区2016—2022年65景Radar-Sat-2降轨图像(5 m分辨率多视精细模式)、2020年6月倾斜摄影Mesh模型(比例尺1∶500)以及2022年3月研究区局部机载LiDAR数据(比例尺1∶500)作为主要数据源,同时采用2.5 m分辨率DEM辅助开展数据处理。各数据源覆盖范围如图2所示。其中,RadarSat-2数据在北京西山具备良好的干涉效果,时序InSAR结果可达毫米级精度;1∶500比例尺Mesh模型和LiDAR点云可获取研究区的高精度三维数据,并有效识别宽度超过20 cm的塌陷坑、地裂缝。利用上述数据可对采空区分布进行初步判断,对滑坡成因及变形情况进行识别,对泥石流发育风险进行分析,满足采空区—滑坡—泥石流链式灾害隐患识别的精度需求。

图2 研究区数据源覆盖范围Fig.2 Data coverage of the study area

2 研究方法及数据处理

研究区内产生链式灾害的关键因素在于滑坡的发育和演化。面向对象分类、深度学习等方法都已被证明可在滑坡光学遥感识别中发挥作用[16-17],人机交互解译能够融入解译人员的专业知识和经验,也是滑坡光学遥感识别的主要方法。本研究利用大比例尺倾斜摄影测量Mesh模型,寻找滑坡与周边区域色调、纹理的差异,提取拉张裂缝、前缘鼓胀引起的坍塌落石,以及植被、地貌、坡度等地表异常,初步锁定滑坡隐患区域。

为了验证光学遥感提取的隐患区域是否存在变形,本研究采用时序InSAR提取地表形变。相关数据处理流程为:首先选取主影像,对研究区的SAR数据进行配准,之后采用振幅离差指数和相关性系数相结合的方法选取PS点,对PS点构建三角网。确定干涉对组合方式后对PS点进行差分干涉,采用解空间搜索方法得到PS点的线性形变速率。针对大气延迟相位、非线性形变相位和噪声在时间域和空间域的不同频率特性,采用不同的滤波方法得到非线性形变速率,最终得到PS点的形变信息。通过内符合精度检验对InSAR结果进行评估,随机选择稳定区域的形变数据进行统计分析,发现直方图呈正态分布(图3),标准差为2.1 mm,表明本次InSAR处理的结果达到毫米级精度,可作为滑坡变形判定依据。经过地理编码后,可与光学遥感圈定的隐患区域叠加分析,实现滑坡识别。

图3 稳定区域LOS方向累计形变量统计直方图Fig.3 Statistical histogram of the deformation in LOS direction in stable regions

值得注意的是,受植被影响,小规模裂隙及古滑坡体在光学图像上的形态特征不明显,即便在大比例尺Mesh模型上也不易识别;且古滑坡体在无特殊因素扰动情况下一般处于活动衰退期或停止期,地表变形较弱,InSAR监测难以发现。利用机载LiDAR多回波的特性,可对植被进行滤除(图4),从而直观地展示地表破坏情况,实现滑坡的识别。本试验共在研究区识别了10处滑坡(含1处古滑坡),如图5所示。经实地查证,10处滑坡均识别正确。

图4 LiDAR数据滤除植被前后建模效果对比Fig.4 Comparison of the modeling results before and after vegetation removal with LiDAR data

图5 滑坡识别结果Fig.5 Landslide identification results

3 研究区滑坡发育特征

研究区10处滑坡规模均为中型,基本参数取值见表1。其中,5处滑坡发育于流域上游,5处滑坡发育于流域中游;2016—2022年各滑坡LOS方向最大地表形变速率为8.1～30.5 mm/a。

表1 研究区滑坡基本参数Table 1 Basic parameters of landslide in the study area

3.1 流域上游滑坡发育特征

S01～S05号滑坡位于东江沟上游,滑坡体平均坡度为41°～48°;滑坡主体均位于龙门组、九龙山组地层单元内,滑坡前缘紧邻窑坡组。该处龙门组、九龙山组岩性以砾岩、粗砂岩为主,岩石较为坚硬。现阶段,上述滑坡体的变形破坏形式以局部滑塌为主,后缘及坡面上可见拉张裂隙。

本研究以S01号滑坡(图6)为例进行分析。该滑坡主体位于龙门组砾岩层中,上部发育4条明显的拉张裂隙,并形成错台,其中以最下方的裂隙A规模最大,形成的错台高约4 m。滑体中部东侧的块体已经下滑,并形成堆积。西侧块体以裂隙A为上部边界,宽约63 m,高约110 m,整体近似呈锥形,总方量约70 000 m3。目前,西侧块体未发生大规模滑动,2016—2022年时序InSAR结果显示:该块体在LOS方向的形变速率最大值为19.7 mm/a。该块体下边界基本位于龙门组砾岩与窑坡组页岩、粉砂岩的交界部位,龙门组砾岩形成的高陡崖壁为块体下滑提供了天然的临空面,使得该块体成为“高危块体”。S03、S04号滑坡也具备类似特征(图7),2处滑坡共发育a、b、c3条主要拉张裂隙,并形成A、B、C3处“高危块体”,总方量约121 000 m3。2016—2022年时序In-SAR结果显示:3处块体LOS方向的形变速率最大值分别为24.1、22.4、20.7 mm/a。

图6 S01号滑坡Mesh模型及变形特征Fig.6 Mesh model and deformation characteristics of Landslide S01

图7 S03、S04号滑坡Mesh模型及变形特征Fig.7 Mesh models and deformation characteristics of Landslide S03 and S04

相对于S01、S03、S04号滑坡,S02、S05号滑坡地表破坏程度较弱(图8),未发育明显的“高危块体”。其中,S05号滑坡按高程自上而下共发育a、b、c、d4条拉张裂隙,形成多级错台,错台高度0.5～2.0 m不等;坡脚崩滑发育(图8中A区)。2016—2022年时序InSAR结果显示:形变区集中分布在c、d2条拉张裂隙之间及前缘崩滑处(图8中A区),LOS方向形变速率最大为12.8 mm/a。S02号滑坡后缘发育1条总长度近400 m的拉张裂隙,从形态上看,裂隙接近贯通,目前在后缘形成高约1 m的台坎;坡脚崩滑广泛发育(图8中B、C、D、E、F区)。2016—2022年时序InSAR结果显示:形变区集中分布在前缘崩滑处(其中B区最为集中),滑坡体中上部也有形变分布,LOS方向形变速率最大处为11.6 mm/a。S02、S05号滑坡前缘集中发育的崩滑导致山体临空,为滑坡体进一步变形下滑提供了空间。根据Mesh模型和LiDAR结果估算,S01～S05号滑坡体总体方量约1 050 000 m3,为东江沟上游最主要的泥石流松散物质来源。

图8 S02、S05号滑坡Mesh模型及变形特征Fig.8 Mesh models and deformation characteristics of Landslide S02 and S05

3.2 流域中游滑坡发育特征

S06～S10号滑坡位于东江沟中游,均发育于窑坡组地层单元内,滑坡体平均坡度为35°～40°,较S01～S05号滑坡偏小。该处窑坡组岩性以页岩、粉砂岩及煤层为主,岩石较软。现阶段,上述滑坡体变形破坏形式以后缘拉裂、前缘鼓胀崩滑为主。

本研究以S08号滑坡体(图9)为例进行分析。该滑坡体后缘发育圈椅状拉裂,目前基本已贯通,形成高约9 m的台坎。坡体中上部发育多条拉张裂隙,形成高0.5～1.5 m不等的错台。坡体中下部崩、滑较发育,集中分布于南侧(图9中虚线圈定范围)。2016—2022年时序InSAR结果显示:该滑坡体变形较大的区域位于滑坡体中上部,其中裂缝集中分布区尤为突出;中下部崩滑区也有较明显的形变分布。整体而言,S08号滑坡体LOS方向形变速率最大为25.3 mm/a,中上部变形大于下部,表现出推移式滑坡的形变特征。S06、S07、S09号滑坡也具备类似特征(图10),3处滑坡后缘均发育圈椅状拉裂缝,形成高2～5 m不等的滑坡台坎;坡体中下部发育小规模崩滑。2016—2022年时序InSAR结果显示:3处滑坡体上部变形均明显大于下部,LOS方向上变形速率最大分别为30.5、28.2、11.4 mm/a。

图9 S08号滑坡Mesh模型及变形特征Fig.9 Mesh model and deformation characteristics of S08 Landslide

图10 S06、S07、S09号滑坡Mesh模型及变形特征Fig.10 Mesh models and deformation characteristics of landslide S06,S07 and S09

S10号滑坡为一古滑坡体,植被覆盖度高,地表变形弱,大比例尺Mesh模型上滑坡周界可识别程度低,时序InSAR结果也无法有效反映滑坡活动(图11(a))。通过对机载LiDAR数据进行植被滤除和纹理增强后(图11(b)),可清晰识别滑坡周界。由图11可知:S10号滑坡后缘发育高约7 m的台坎,滑坡边界清晰,坡脚普遍发育小型崩滑;现阶段坡体整体变形较小,仅下部存在局部微变形,根据2016—2022年时序InSAR结果,LOS方向上变形速率最大为8.1 mm/a。

图11 S10号滑坡Mesh模型、机载LiDAR数据及变形特征Fig.11 Mesh model,airborne LiDAR data and deformation characteristics of S06,S07 and S09 landslide

根据Mesh模型和LiDAR结果估算,S06～S10号滑坡体总体方量约1 480 000 m3。这些滑坡体有可能成为东江沟中游重要的泥石流沿程物源的补给来源,同时还存在堵塞沟道的风险。

4 东江沟链式灾害隐患特征与易发性评价

根据北京市2022年6月公布的“突发地质灾害隐患点台账”,东江沟为一典型的沟谷型泥石流隐患,发生泥石流的可能性为“中”,威胁对象为下游居民点和道路。2012年7月21日、2016年7月20日,该隐患点发生两次小规模泥石流,累计毁坏房屋4间,农田200 000 m2,道路1 km,无人员伤亡,最大冲淤变幅达2.1 m。

4.1 滑坡成因及驱动力分析

研究区内煤矿开采历史近千年,小煤窑众多。地下开采范围、开采深度、开采厚度等信息的严重缺失一直是区内采空塌陷调查及研究的难点。本研究利用大比例尺Mesh模型,在10处滑坡体及周边区域共提取了地面塌陷坑116个、地裂缝41条、废弃井口44处(图12),说明上述10处滑坡体地下广泛发育采空区。虽然区内煤矿开采活动已于2008年停止,目前广泛的地表沉陷已趋缓,但老采空区山体结构被破坏,具备产生持续变形并形成滑坡的条件。

图12 滑坡周边废弃井口、塌陷坑及地裂缝分布Fig.12 Distribution of abandoned cave mouths,collapse pits and ground cracks around landslides

一般地,地形地貌、地质成分结构等地质因素决定了地质灾害的易发程度;冻融作用、降雨渗流、地震作用和人类工程活动等构成了地质灾害的诱发条件。地质灾害是地质因素、诱发条件耦合作用的结果[18]。

根据北京市2022年6月公布的“突发地质灾害隐患点台账”及实地调查成果,具备上述10处滑坡类似地层岩性及地形地貌条件的其他山区,并未见有典型滑坡发育(采空区除外)。近20 a来,北京西山地区发生多次1～2级地震,但至今未见由于微震导致地质灾害的报道,并且自2008年煤矿停采后,研究区内无其他人类工程活动干扰。根据裴昕等[19]研究,北京西山老采空区内处于等速变形阶段的滑坡隐患,变形速率与降雨量呈正相关[19]。综合分析可知,研究区内滑坡发育并持续变形是地下采空区和降雨共同作用的结果。

4.2 地形条件分析

东江沟下切强烈,为典型的“V”形谷。根据DEM测算,流域面积3.87 km2,山坡平均坡度31.8°,相对高差963 m,主沟纵坡度25.2%。为了对比东江沟与北京地区其他泥石流沟的地形条件,本研究收集了北京市历史上曾经发生过灾害的267条沟谷型泥石流的地形参数数据,如图13所示。根据图13,北京地区最易发生泥石流的沟道条件为:山坡平均坡度25°～32°,流域面积0.2～5.0 km2,相对高差大于500 m,主沟纵坡度大于21.3%。由此可知:东江沟具备北京地区最易形成泥石流的地形条件。

图13 北京地区历史泥石流沟地形参数统计结果Fig.13 Statistical results of terrain parameters of debris flows in Beijing

4.3 泥石流物源条件分析

通过遥感判释,东江沟的泥石流物源可分为煤矿开采堆积物(煤矸石)、自然风化堆积物(残坡积物)及崩滑灾害(新识别滑坡及其他崩滑堆积物)(图14),估算净储量可达2 800 000 m3(约合720 000 m3/km2);其中煤矿开采堆积物净储量151 000 m3,自然风化堆积物净储量164 000 m3,崩滑灾害估算净储量2 485 000 m3(10处滑坡估计方量2 445 000 m3,其他9处小规模崩滑堆积物净储量40 000 m3)。通过对比不难发现,10处滑坡提供的物源量占据了总物源量的87.3%,成为东江沟最主要的泥石流物质来源。

图14显示,S01、S02、S03、S04、S05号滑坡均位于流域上游沟道的汇水点位置(图14(b)),一旦在暴雨情况下产生失稳下滑,堆积物将直接进入沟道,成为泥石流启动的物源。此外,煤矿开采堆积物及自然风化堆积物大量堆积于主沟中下游及各级支沟沟道内,能够为泥石流提供的物源补给长度比约80%;且物源平均厚度较大,尤其是煤矿开采堆积的矸石厚度通常可达10～15 m,个别矸石山堆积高度近30 m。

图14 东江沟物源分布Fig.14 Distribution of debris flow source of Dongjiang Gully

4.4 堵溃风险分析

东江沟地形切割深、沟道狭窄,除了主沟外,还发育东、西两条主要支沟(支沟E、支沟W)。沟道内的部分矸石堆对河道堵塞严重。以支沟W局部段为例,其纵剖面如图15所示。由图15可见A、B2处矸石堆形成了两道高近10 m的天然坝体,不利于沟道排泄。

图15 支沟W局部段纵剖面示意Fig.15 Schematic of the profile of local section of branch W

此外,S06、S07号滑坡分别发育于支沟W与支沟E东侧山坡,S08、S09号滑坡发育于主沟东侧山坡,S10号滑坡发育于主沟西侧山坡(图14(c))。S06、S07、S08、S09号滑坡变形速率较大,在不利因素影响下易失稳下滑。这4处滑坡点位的沟道横截面(其中滑坡体厚度为示意图,不代表其真实厚度)如图16所示。分析可知:沟道底部均十分狭窄(宽20～50 m不等),若暴雨条件下滑坡体突然下滑,可在短时间内形成高数十米的堰塞坝,坝体溃决后极易引发泥石流。

图16 S06、S07、S08、S09号滑坡点沟道横截面示意Fig.16 Schematic of gully cross section at S06,S07,S08 and S09 landslide

4.5 泥石流易发性评价

通过上述分析,东江沟具备形成泥石流的良好地形条件。参考《地质灾害危险性评估技术规范》(DB11/T 893—2021)[20],对东江沟发生泥石流的可能性进行了评价。发现在考虑10处滑坡的情况下,泥石流发生的可能性有显著提升。评价因子及得分见表2。

表2 东江沟泥石流易发性评价结果Table 2 Evaluation results of the susceptibility of debris flow in Dongjiang Gully

泥石流易发程度等级评分按下式计算:

式中,Q为易发程度总得分;Fi为各因子得分。

当Q＞105时,发生泥石流的可能性“大”;当76≤Q≤105时,发生泥石流的可能性“中”;当Q＜76时,发生泥石流的可能性“小”。

在东江沟泥石流发生可能性评价中,10处滑坡对“崩塌、滑坡及水土流失严重程度”“松散物储量”“泥沙沿程补给长度”“河沟堵塞程度”等多个评价因子作出较大贡献。通过计算,当考虑这些泥石流易发因子时,东江沟Q值为112,发生泥石流的可能性“大”;反之,Q值为96,发生泥石流的可能性“中”。考虑滑坡影响的评价结果更准确,更符合实际。这也证明了通过综合遥感技术,对于采空区—滑坡—泥石流链式灾害的识别是必要的,有助于科学、全面地认识采空区地质灾害风险。