李林澤, 常 鳴, 李宏杰, 胡錦鳳, 徐恒志, 劉沛源

(1. 成都理工大學 地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室, 四川 成都 610059;2. 山東正元數(shù)字城市建設有限公司, 山東 煙臺 264670)

0 引言

我國是世界上震后地質(zhì)災害活動最頻繁、損失最嚴重的國家之一,其中泥石流災害是震后山區(qū)地質(zhì)災害中危險性最大的災害之一。強烈的地震作用在震后短期內(nèi)具有誘發(fā)大量山體滑坡并產(chǎn)生松散堆積物分布在溝道內(nèi)及溝道兩岸的能量,對山體造成的內(nèi)部損傷影響亦可達數(shù)百年之久[1-2]。與此同時,地震也大大降低了誘發(fā)泥石流災害的降雨閾值,導致災難性的泥石流事件頻繁發(fā)生,規(guī)模加大,人民的生命財產(chǎn)安全受到極大地危害和威脅[3-4]。如2008年5月12日汶川特大地震后在地震的影響范圍內(nèi)接連暴發(fā)了大規(guī)模泥石流,該地震10年后雖然仍有大規(guī)模泥石流暴發(fā),但在震后短期內(nèi)泥石流啟動降雨閾值普遍降低的情況下雨季暴發(fā)泥石流的頻次更高,規(guī)模更大,受災更嚴重。因此,對震后短期內(nèi)可能因強降雨引發(fā)泥石流的運動過程和沖出規(guī)模需要更為深入的研究。

近二十年來,國內(nèi)外大量學者對泥石流運動過程及沖出規(guī)模開展的研究可分為三類:基于體積方量的預測模型、動態(tài)預測模型和基于地形參數(shù)的預測模型[5-9]。隨著科技水平的提高,AI、數(shù)值計算等方法不斷出現(xiàn)并更新迭代,其中泥石流數(shù)值計算是理論基礎和計算數(shù)學的結合,將現(xiàn)實中的泥石流運動全過程簡化為離散方程后進行求解,加之數(shù)值計算平臺搭建完善后使用成本低且操作簡捷,近年來被大量運用到泥石流研究中[10-13]。但如FLO-2D、FLO-3D、VolcFlow、PFC等數(shù)值計算模型均是以獲得泥石流流深及流速為目的的單一沖出模擬,不能很好地耦合降雨、啟動、侵蝕、運移、堆積整個動態(tài)過程,并且忽略了泥石流運動過程中尤為重要的侵蝕沖刷。

OpenLISEM數(shù)值計算模型是一個集降雨輸入、徑流、啟動、侵蝕、運移、沖出、堆積于一體的綜合水文模型,其原理是耦合分布式流域水文模型與二相動量守恒定律,特點是把流域劃分成若干網(wǎng)格區(qū)域,對每個網(wǎng)格區(qū)域分別輸入不同植被、土壤、巖土參數(shù)和高程等屬性,輸入實時降雨條件后分別計算產(chǎn)流量;通過比較相鄰網(wǎng)格的高程確定各網(wǎng)格的流向,演算降雨過程中泥石流的啟動、運移、堆積整個過程[14]。該模型不但考慮了地形、植被、物源和水文過程,同時也避免了以往常規(guī)模擬軟件需要的泥石流流量過程線計算與輸入、啟動點位置選取等問題的不準確性,能更加真實地重現(xiàn)泥石流暴發(fā)時的場景。

受瀘定6.8級地震影響,四川省瀘定縣落井溝流域新增大量物源,未來雨季,在泥石流激發(fā)降雨閾值降低的條件下研究區(qū)極有可能暴發(fā)泥石流,威脅溝口聯(lián)合村隧道、在建的瀘定—石棉高速與溝口居民生命財產(chǎn)安全。本文以落井溝為研究區(qū),通過分析其物源分布特征,選取可以較好地模擬泥石流運動侵蝕全過程的OpenLISEM數(shù)值計算模型對不同降雨強度下的落井溝泥石流運動全過程進行模擬預測,為該區(qū)域的泥石流危險評價、預警預報及工程治理提供具有一定適用性和參考性的方法。

1 研究區(qū)概況

1.1 區(qū)域環(huán)境地質(zhì)條件

落井溝位于四川省瀘定縣得妥鎮(zhèn)中部,溝口坐標102°10′42″E,29°33′33″N,其形成區(qū)距2022年9月5日瀘定6.8級地震震中約1.2 km,位處Ⅸ級烈度區(qū)(地理位置如圖1)。流域較小,總面積約3.4 km2,地勢復雜,東高西低,地形起伏范圍為1 089～2 572 m,主溝長約2.7 km,流向自東南向西北,平均縱坡降約為390‰,溝谷狹長,兩岸坡度較陡,為明顯的“V”字型,屬典型的高山峽谷區(qū)。流域內(nèi)發(fā)育有兩條小型支溝,于主溝中部交匯并最終流入與溝口正交且流向為東北至西南的大渡河。

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

圖2 瀘定縣2022-09-05—2023-04-15地震概況Fig.2 General situation of earthquakes in Luding County during 2022-09-05—2023-04-15

圖3 落井溝坡度分級與統(tǒng)計Fig.3 Slope grading and statistics in Luojing Gully

流域年平均氣溫15.7 ℃,年最低氣溫-4.5 ℃,年平均降雨量664.4 mm。從山巔至谷底可觀察到明顯的氣候垂直分帶現(xiàn)象,海拔低處氣候較溫暖濕潤,高處較嚴寒干燥,屬典型的亞熱帶季風氣候。研究區(qū)屬華南地層,形成區(qū)主要發(fā)育花崗巖、砂巖;流通區(qū)及堆積區(qū)主要發(fā)育二疊、三疊系變質(zhì)巖,以板巖為主,同時整個流域發(fā)育較多沉積巖,第四系堆積物雜亂且分布不均。得妥鄉(xiāng)位于川滇南北向構造帶北段,為多組構造的交叉部位,構造環(huán)境復雜,新構造運動強烈,主要斷裂為鮮水河斷裂東南段(也稱磨西斷裂)及大渡河斷裂,且大渡河斷裂橫穿落井溝堆積區(qū)[15]。

1.2 地震概述

2022年9月5日12時52分,四川省甘孜州瀘定縣發(fā)生6.8級地震,震源深度16 km,震中位于29.59°N,102.08°E,震源機制解以走滑型為主。對四川省地震局提供的震訊進行統(tǒng)計,發(fā)現(xiàn)從2022年9月5日6.8級地震截止到2023年4月15日,瀘定縣共發(fā)生32次地震,其中包括3次5級及以上地震,即2022年9月5日6.8級、2022年10月22日5.0級及2023年1月26日5.6級地震,為典型的“主震-余震型”地震。

經(jīng)現(xiàn)場調(diào)查,“9·5”瀘定地震直接造成落井溝內(nèi)產(chǎn)生大量滑坡,溝道內(nèi)堆積較多的松散堆積物,后續(xù)接連暴發(fā)的31次地震更使得溝道兩側山體內(nèi)部損傷程度及溝道堆積物松散程度逐步增加。泥石流物源在震后得到大量補給,一旦暴發(fā)將對溝口聯(lián)合村隧道及在建的瀘定—石棉高速的使用和建造造成巨大影響,泥石流若發(fā)生堵江也可能對溝口東南部約11 km處正在運營的大崗山水電站造成難以估量的損失。

2 落井溝泥石流特征

2.1 地形特征

陡峻的溝道地形條件是形成泥石流至關重要的一環(huán)[16]。落井溝流域發(fā)育有兩條小型支溝,主溝頂部至溝口平均縱坡降為390‰,地勢陡峭,為明顯的“V”字型,屬高山峽谷深度切割區(qū)域。其溝谷狹長,兩岸坡度較陡,統(tǒng)計發(fā)現(xiàn)近70%的流域面積分布在坡度為30°～50°范圍內(nèi),滑坡松散物質(zhì)運移到溝道的地形條件極好。這樣的地形條件對溝道與坡面物源穩(wěn)定程度及降水形成徑流的水動力強度均具有正向的交互效應,利于泥石流發(fā)生。

2.2 物源特征

近年來對泥石流物源的識別技術愈加豐富,如通過遙感高精度影像進行滑坡的人工解譯[17]與通過合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar,ISAR)技術進行潛在滑坡的隱患識別[18],這些識別技術無一例外都需通過現(xiàn)場復核以檢驗其適用性及準確性,而一些極為困難且危險的野外調(diào)查及測繪工作難以正常進行,因此無人機低空數(shù)字攝影野外測量技術開始廣泛應用[19]。

瀘定縣6.8級地震后,通過遙感解譯技術獲取落井溝流域新增滑坡37處,面積約0.125 km2。震后新增滑坡位置分布見圖4(b)。

圖4 滑坡震前震后影像對比Fig.4 Comparison between landslide images before and after the earthquake

震后不久對其進行現(xiàn)場調(diào)查并使用搭載長測程激光雷達的無人機進行仿地飛行采集數(shù)據(jù),通過大量點云數(shù)據(jù)處理等方式獲得了分辨率為0.2 m的高精度機載激光雷達(Light Detection and Ranging,LiDAR)數(shù)據(jù),包括數(shù)字高程模型(Digital Elevation Model,DEL)及數(shù)字正射影像(Digital Orthophoto,DO)數(shù)據(jù),用于輔助識別落井溝流域內(nèi)因植被覆蓋未能識別的崩滑體,并作為基礎數(shù)據(jù)進行后期泥石流運動全過程數(shù)值計算。通過圖5中b1與b2對比,發(fā)現(xiàn)去植被后的LiDAR數(shù)字高程模型能識別出影像不易察覺的滑坡堆積物,同時,由于較多滑坡碎屑流僅是淺層物質(zhì)剝落,故通過LiDAR航測處理得到的山體陰影數(shù)據(jù)進行精準識別,剔除未產(chǎn)生與產(chǎn)生極少物源的滑坡,得到實際滑坡堆積面積約7.798 萬m2(初始滑坡分布見圖1,處理后的滑坡分布見圖5)。

圖5 LiDAR數(shù)字正射影像及去植被的數(shù)字高程模型(圖a1與b1分別是落井溝兩處去植被后 的溝道與坡面滑坡LiDAR數(shù)據(jù);圖a2與b2分別是落井溝兩處未去植被的LiDAR數(shù)字 正射影像數(shù)據(jù);c為滑坡取樣點)Fig.5 LiDAR digital orthophoto and digital elevation model for vegetation removal

結合滑坡體積-面積經(jīng)驗公式[20][式(1)]得到震后落井溝溝道及坡面新增的實際滑坡松散堆積物方量為15.52 萬m3,且主要集中在主溝中下段(如圖5)。溝道堆積物易被強降雨作用下形成的地表徑流沖刷侵蝕,溝岸坡腳堆積物被沖刷帶走后形成臨空面又會不斷失穩(wěn)直至運移到溝道,增加溝道物源堆積,且隨海拔下降物源侵蝕會不斷增加,進而暴發(fā)泥石流災害。

v=1.315*A1.201 8

(1)

式中:A為同震滑坡面積(m2);v為同震滑坡體積(m3)。

2.3 水文地質(zhì)特征

震后大量崩滑體在強降雨作用下極易轉化為泥石流[21]。瀘定縣年平均降雨量664.4 mm,雨季降雨充沛,最高月降雨量可達167 mm,受“V”型峽谷地貌影響,該流域匯水形成徑流的能力較強,加之大渡河從溝口流經(jīng),地下水豐富,流域內(nèi)主要為松散堆積層孔隙水和基巖裂隙水,補給方式通常為大氣降水。

泥石流災害的發(fā)生是地形、物源及降雨等多因素共同作用的結果[22]。通過對落井溝地形、物源及降雨三個影響因素分析發(fā)現(xiàn),近70%的流域面積分布在坡度為30°～50°范圍內(nèi),“9·5”瀘定震后溝道及坡面的實際新增滑坡松散堆積物方量為15.52 萬m3,且主要集中在主溝中下段,加之研究區(qū)雨季降雨充沛,地下水豐富,在未來的強降雨條件下具備暴發(fā)泥石流的初始條件,故對落井溝泥石流啟動、物質(zhì)運移及堆積規(guī)模借助數(shù)值計算工具進行深入分析。

3 泥石流運動演化分析

3.1 數(shù)值計算基本原理

OpenLISEM數(shù)值計算模型通過建立柵格單元,融合SWATRE滲透模型、理查茲方程、圣維南二維非恒定流方程、二相動量守恒定律及無限邊坡理論,實現(xiàn)災害啟動、侵蝕堆積和控制流體的微分方程求解,以此模擬地表水沙平衡,其具有很明顯的時空差異,是一個離散的二維數(shù)值模型。該模型可用于模擬流域在降雨期間或降雨后短時間內(nèi)的物源變化、土地利用類型變化、防治工程等制約因素對泥石流的沖出趨勢、物源侵蝕、洪水等動態(tài)流動趨勢的影響。

3.2 數(shù)據(jù)獲取

模擬所需數(shù)據(jù)主要包括物源數(shù)據(jù)、地形數(shù)據(jù)和降雨數(shù)據(jù),其中地形數(shù)據(jù)利用QGIS平臺的轉換功能將0.2 m的機載LiDAR數(shù)字高程模型(DEM)轉換成模擬所需map格式文件,其他地形數(shù)據(jù)利用PCraster 軟件和Arcgis平臺制作成map格式文件,包括流向網(wǎng)、地勢梯度及表面粗糙度等。物源的物理力學特性指標包括滲透系數(shù)、初始含水率、內(nèi)摩擦角、孔隙率、內(nèi)部黏聚力、平均粒徑及特征粒徑D50和D90等,需通過野外調(diào)查及現(xiàn)場取樣進行室內(nèi)試驗。降雨數(shù)據(jù)僅需輸入泥石流暴發(fā)期間的降雨過程,但由于震后落井溝尚未進入雨季且未暴發(fā)泥石流,無法獲取降雨啟動泥石流后運動過程中的實際降雨歷時數(shù)據(jù),故結合5%、2%、1%降雨強度得到落井溝洪水流量過程線,輸入到OpenLISEM數(shù)值計算模型中進行模擬。

(1) 物源數(shù)據(jù)

通過0.2 m分辨率的LiDAR航拍數(shù)據(jù)處理得到的落井溝溝道與坡面滑坡松散堆積物體積方量約15.52 萬m3。此外,還需獲取物源的土體級配特征D50及D90進行建模,通過篩除大于20 mm的現(xiàn)場樣品并進行室內(nèi)顆粒分析試驗最終得到土體級配曲線,獲取了D50及D90分別為5 mm、11 mm(圖6)。取樣點在流域中的位置分布見圖5中c點。

圖6 落井溝滑坡堆積物級配曲線(<20 mm)Fig.6 Grading curve of landslide deposit in Luojing Gully (<20 mm)

通過不同土體的初始含水率、孔隙度及滲透系數(shù)等參考值[23]加上現(xiàn)場調(diào)查與室內(nèi)試驗綜合得到如表1所列的物源指標。

表1 物源特征指標及取值

(2) 洪峰流量

OpenLISEM進行泥石流數(shù)值計算的降雨數(shù)據(jù)不需要通過計算獲取泥石流流量過程,僅需提供實際降雨過程便能結合物源及地形參數(shù)進行自動化模擬,得出最符合實際的泥石流流量過程。但由于落井溝尚未暴發(fā)泥石流,故采用四川省中小流域暴雨洪水計算手冊進行計算,獲得落井溝泥石流20年、50年、100年一遇暴雨重現(xiàn)期的洪水流量水文圖。采用以下公式(2)～(9)進行計算:

(2)

(3)

(4)

(5)

m=0.221θ0.204

(6)

(7)

(8)

(9)

式中:QB為暴雨流量(m3/s);ψ為洪峰徑流系數(shù);i為最大平均暴雨強度(mm/h);F為流域面積(km2);s為暴雨雨力(mm/h);n為暴雨公式指數(shù);τ為流域匯水時間(h)。K6、K24為不同重現(xiàn)期的模比系數(shù),可從皮爾遜Ⅲ型曲線模比系數(shù)表中獲取;θ為流域特征參數(shù);m為匯流參數(shù);τ0為當ψ=1時的流域匯流時間(h);μ為產(chǎn)流參數(shù)(mm/h)。相關取值列于表2,洪水流量過程線見圖7。

表2 清水洪峰流量相關參數(shù)

圖7 落井溝洪水流量過程線Fig.7 Flood flow hydrograph in Luojing Gully

3.3 模擬步驟、結果及分析

模擬過程需提前設置四個基礎全英文命名的文件,包含地形和物源的map文件、模擬運行的命令流run文件、降雨數(shù)據(jù)rain文件與存儲輸出結果的res文件,設置完畢后對落井溝泥石流基于OpenLISEM模型進行運動全過程模擬。

結果顯示:5%、2%、1%降雨強度下,落井溝泥石流運動過程的最大泥深分別為6、8、11 m,最大流速分別為5、7、10 m/s,整個運動過程的泥沙堆積顯示溝道內(nèi)淤積較多泥石流固體物質(zhì),推測受落井溝流域地形影響,于拐彎處發(fā)生彎道超高,降低了泥石流流速,導致泥石流淤積溝道的固體物質(zhì)較多(圖8)。泥沙隨空間的變化關系呈現(xiàn)為逐步從主溝上端運移至溝口(圖9)。

通過分析落井溝泥石流模擬的最終沖出規(guī)模發(fā)現(xiàn):5%、2%降雨頻率下落井溝泥石流雖沖出溝口但對大渡河的堵江程度較小,而1%降雨頻率下落井溝泥石流造成大渡河半堵江狀態(tài)。

4 結語

通過遙感解譯、現(xiàn)場調(diào)查及機載LiDAR航拍,對2022年9月5日瀘定6.8級地震后落井溝泥石流新增滑坡的發(fā)育特征進行分析,結合OpenLISEM數(shù)值計算模型對落井溝泥石流在5%、2%、1%降雨強度下的運動全過程進行模擬,得到運動過程中的最大堆積深度、最大流動速度及最終沖出范圍,并給出1%降雨強度下泥石流的泥沙演化全過程,綜合反映出落井溝泥石流的運動侵蝕堆積情況,并對大渡河的堵塞程度給出建議。主要結論或建議如下:

(1) 受瀘定6.8級地震直接影響,流域新增滑坡37處,滑坡面積約0.125 km2。通過現(xiàn)場調(diào)查及0.2 m高精度分辨率機載LiDAR數(shù)據(jù)對滑坡進行優(yōu)化識別,并結合滑坡體積-面積經(jīng)驗公式得到落井溝溝道及坡面的實際滑坡松散堆積物方量為15.52萬m3,且主要集中在主溝中下段。

(2) 基于OpenLISEM數(shù)值模擬平臺,結合精細化地形、精準物源及5%、2%、1%降雨強度對落井溝泥石流進行運動全過程模擬,獲得運動過程最大泥深分別為6、8、11 m,最大流速分別為5、7、10 m/s。通過對溝口最終的堆積狀況分析,5%、2%降雨頻率下落井溝泥石流雖沖出溝口但堵江程度極小,1%降雨頻率下落井溝泥石流造成大渡河半堵江狀態(tài)。

(3) 落井溝泥石流溝域面積較小,即便在雨季暴發(fā)對溝口聯(lián)合村隧道和在建的瀘定—石棉高速的使用和建造在大渡河下游庫岸相連的大崗山水電站的正常運營及沿線居民的生命財產(chǎn)安全造成的威脅仍然較小。但值得注意的是,自2022年9月5日瀘定地震以來共發(fā)生大小地震32起,包括3次5級以上地震,這種小而多的地震對溝道及坡面滑坡物源的各方面影響(尤其是潛在滑坡臨界失穩(wěn)狀態(tài))有著肉眼不易察覺的特征,雨季豐沛的降雨及形成的徑流對該流域的潛在物源也有著較大影響,對瀘定震區(qū)此類物源豐富、地勢陡峻的溝道仍需密切關注。