采空區—滑坡—泥石流鏈式災害隱患綜合遙感識別與評價

焦潤成 郭學飛 南 赟 楊紅磊 曹 穎 王晟宇 閆 馳韓建鋒 馬曉雪 趙丹凝 倪 璇 馬 馳

(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°;滑坡主體均位于龍門組、九龍山組地層單元內,滑坡前緣緊鄰窯坡組。該處龍門組、九龍山組巖性以礫巖、粗砂巖為主,巖石較為堅硬?，F階段,上述滑坡體的變形破壞形式以局部滑塌為主,后緣及坡面上可見拉張裂隙。

本研究以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號滑坡偏小。該處窯坡組巖性以頁巖、粉砂巖及煤層為主,巖石較軟?，F階段,上述滑坡體變形破壞形式以后緣拉裂、前緣鼓脹崩滑為主。

本研究以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,發生泥石流的可能性“中”?？紤]滑坡影響的評價結果更準確,更符合實際。這也證明了通過綜合遙感技術,對于采空區—滑坡—泥石流鏈式災害的識別是必要的,有助于科學、全面地認識采空區地質災害風險。

5 結 論

(1)利用大比例尺Mesh模型、機載LiDAR及時序InSAR開展了礦區采空區—滑坡—泥石流鏈式災害高精度綜合遙感識別是有效且可靠的。

(2)隨著礦區地面沉陷持續發展,有可能進一步形成滑坡隱患。本研究在東江溝采空區內準確識別出10處中型滑坡,測得其在2016—2022年間的LOS方向地表形變速率為8.1～30.5 mm/a。

(3)10處滑坡成為東江溝的主要物源,并強烈影響了多項泥石流易發因子,使得東江溝發生泥石流的可能性由“中”提升為“大”。通過綜合遙感技術開展采空區—滑坡—泥石流鏈式災害識別和評價,有利于科學、全面認識采空區地質災害風險。

(4)北京西山存在諸多地下采空區,部分區域地面沉陷十分發育。有必要在區內開展高精度鏈式災害識別與評價,支撐地質災害防治決策。