基于數(shù)值模擬的城鎮(zhèn)地質(zhì)災(zāi)害危險性評價方法

龔凌楓， 徐 偉， 鐵永波， 盧佳燕,2， 張 玙， 高延超

(1.中國地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心，四川 成都 610081； 2.中國地質(zhì)大學(xué)(北京)，北京 100083)

0 引言

隨著城鎮(zhèn)化水平的不斷提升，對重點區(qū)城鎮(zhèn)地質(zhì)災(zāi)害風(fēng)險精細化評價需求及評價結(jié)果的準(zhǔn)確性要求進一步提升[1-2]。當(dāng)前危險性評價方法主要包括定性評價方法(專家打分法[3]、模糊綜合評判[4-7])、半定量評價方法(信息量法、層次分析法[7-8]、證據(jù)權(quán)法[9]、頻率比法、邏輯回歸模型[10]、判別分析模型、聚類分析模型[11]、人工智能模型[12-13])及定量評價方法(無限斜坡模型[14]、3D極限平衡方法、斜坡水文模型[15])，上述方法都是地質(zhì)災(zāi)害危險性評價中常用的分析方法。

定性評價方法需豐富的災(zāi)害點數(shù)據(jù)支撐，且對災(zāi)害點數(shù)據(jù)的分布均勻性有較高的要求，否則難以避免因調(diào)查范圍、精度差異造成的評價結(jié)果與災(zāi)害點分布的自然屬性不符，該評價方案常用于小比例尺的地質(zhì)災(zāi)害危險性評價[16]。定量評價方法是以巖土體力學(xué)性質(zhì)、氣象及地表水文、地下水滲流等[6,17-20]與災(zāi)害形成直接相關(guān)的特性為基礎(chǔ)，針對不同區(qū)域地質(zhì)災(zāi)害主控因素輔以對應(yīng)模型進行計算，適用于重點調(diào)查區(qū)精細化地質(zhì)災(zāi)害危險性評價，可實現(xiàn)大比例尺地質(zhì)災(zāi)害危險性評價。

目前，較多的地質(zhì)災(zāi)害單體危險性評價應(yīng)用為重要災(zāi)害點危險性評價提供了較多的案例，但區(qū)域性危險性評價較少采用數(shù)值計算方法，主要是受到群發(fā)性地質(zhì)災(zāi)害評價方法的制約，而單體地質(zhì)災(zāi)害應(yīng)用于區(qū)域性評價不可避免造成工作量劇增，因而難以推廣應(yīng)用。淺層滑坡、群發(fā)性泥石流優(yōu)勢發(fā)育地區(qū)的數(shù)值計算方法較為成熟[15-16]，可以用于精度要求高的重點區(qū)和典型城鎮(zhèn)地質(zhì)災(zāi)害危險性評價。

本文以典型紅層地區(qū)淺層滑坡為研究對象，選用滑坡流體力學(xué)觸發(fā)模型(hydromechanical triggering of landslides,LHT)進行重點研究區(qū)地質(zhì)災(zāi)害危險性評價。2012年8月31日，四川涼山州喜德縣遭遇百年一遇的特大暴雨，孫水河上游米市監(jiān)測當(dāng)日最大降雨量達到149.2 mm； 強降雨發(fā)生后，在孫水河流域兩河口以上流域的山坡多處可見“山剝皮”現(xiàn)象，大量植被被破壞，對生態(tài)環(huán)境產(chǎn)生較大的負(fù)面影響。本文以“8·31”群發(fā)性滑坡泥石流事件中災(zāi)害解譯點為驗證數(shù)據(jù)，探索不同降雨頻率條件下重點區(qū)淺層滑坡危險性評價方法及基于數(shù)值模擬的不同降雨頻率下1∶1萬地質(zhì)災(zāi)害危險性定量計算方法。

1 研究區(qū)概況及數(shù)據(jù)來源

研究重點區(qū)為四川省涼山州喜德縣洛哈鎮(zhèn)的集鎮(zhèn)區(qū)，面積44.5 km2，其地理位置、地質(zhì)概況及地質(zhì)災(zāi)害點分布如圖1所示。研究區(qū)具有冬季干燥無嚴(yán)寒，夏季溫涼多雨，四季不分明，氣溫日差大，年差較小，風(fēng)多，夜雨多，冰雹多的典型亞熱帶季風(fēng)和高原氣候特征。

圖1 研究區(qū)地質(zhì)概況及地質(zhì)災(zāi)害點分布

1.1 工程地質(zhì)概況

區(qū)域內(nèi)主要構(gòu)造為米市向斜，軸線經(jīng)過米市鎮(zhèn)、洛哈鎮(zhèn)，呈NE27°，南起額尼村，向北延伸。向斜軸部及兩翼地層均為上白堊統(tǒng)雷打樹組上段(K2l2)砂質(zhì)泥巖。東翼傾向北南270°～318°，傾角2°～19°； 西翼傾向南東120°～145°，傾角10°～15°。局部地形起伏大，沿路及沿河陡坡區(qū)域卸荷裂隙發(fā)育，岸坡風(fēng)化層厚度大，為小規(guī)模的崩塌及大量的淺層滑坡提供了良好的地形條件，且由于該區(qū)域植被覆蓋度較差，尤其是海拔較高的位置山地裸露，為長歷時及強降雨作用下淺層滑坡的發(fā)生提供了的條件。

研究區(qū)巖性為灰黃色質(zhì)泥質(zhì)砂巖和紫紅色泥巖。河流及筑路切坡處陡坡基巖以泥巖為主，夾薄層砂巖。強風(fēng)化基巖破碎呈碎塊狀、中風(fēng)化基巖呈短柱狀，強度較低。人類居住區(qū)和耕值區(qū)的地形受人類活動影響改變較大，主要表現(xiàn)為農(nóng)作物耕植過程中在坡體上開挖形成梯田狀小土坎，長度一般為8～15 m，坡高一般為0.5～1.5 m?？傮w來說，區(qū)內(nèi)人類工程活動較強烈，對地形的改變?yōu)榛碌幕瑒犹峁┝擞欣淖冃慰臻g。

1.2 數(shù)據(jù)來源

研究區(qū)數(shù)據(jù)主要包括研究區(qū)高分二號衛(wèi)星數(shù)據(jù)(分辨率為0.8 m，拍攝時間為2020年)、1∶10 000地形圖、1∶50 000地質(zhì)圖、1∶50 000地質(zhì)災(zāi)害分布圖、綜合遙感解譯圖及野外地質(zhì)調(diào)查資料等，可通過上述資料獲取該區(qū)域危險性評價基礎(chǔ)圖件(圖2)。通過前期資料收集及野外調(diào)查，區(qū)內(nèi)原有災(zāi)害點3處，初步光學(xué)遙感解譯及InSAR識別獲取滑坡點17處，通過精細光學(xué)遙感解譯共獲取淺層滑坡點409處，滑坡面積共計1.07 km2，占研究區(qū)總面積的2.4%。

圖2 危險性評價基礎(chǔ)數(shù)據(jù)

2 研究方法

本文以淺層滑坡發(fā)育區(qū)地質(zhì)災(zāi)害危險性精細化評價為目標(biāo)，選取數(shù)值模擬方法進行評價。通過對國內(nèi)外已有成果分析可知，危險性精細化評價主要在于巖土體穩(wěn)定性計算方法選取及不同頻率條件下的降雨歷程曲線計算。

2.1 數(shù)值計算方法選取

對重點區(qū)調(diào)查采用基于柵格計算斜坡穩(wěn)定性的LHT模型[21]。LHT是一種基于力學(xué)公式的數(shù)值模型，能完整模擬坡體在降雨條件下土體發(fā)生漸進性破壞，最終整體失穩(wěn)的過程。模型把覆蓋在基巖上的土體離散成一組機械連接的六角形土柱[22]，通過物理模型將土體的水文狀態(tài)和力學(xué)狀態(tài)聯(lián)系起來，以摩擦力、黏聚力、毛細力及植物根強度來保持坡面土體的穩(wěn)定性。

六角土柱的基本受力模型如圖3所示，在降雨條件下，單個六棱柱的力學(xué)參數(shù)包括土柱體積重度W、正應(yīng)力σN和不飽和抗剪強度τS,計算公式為

W=Hsd[θρw+(1-φ)ρr]gcosβsinθ，

(1)

σN=Hs[θρw+(1-φ)ρr]gcos2β，

(2)

τS=csoli+{Hsd[θρw+(1-φ)ρr]gcos2β-
χhρwg}tanγ。

(3)

式中：Hsd為土壤厚度，m；θ為體積含水量,無量綱；h為毛細管頭，m；ρw為水的密度，kg/m3；ρr為土體密度，kg/m3；β為坡角,(°)；φ為孔隙率,%；Csoli為土壤黏聚力,kPa；χ為剪切強度與毛細管壓力的比例因子,無量綱；γ為土體內(nèi)摩擦角,(°)。

圖3 六棱柱觸發(fā)模型基本單元[23]

該模型將滑坡過程建模為一個瞬時事件，沒有去描述滑坡土體運動的具體過程，也沒有跟蹤滑坡的質(zhì)量變化，通過由微觀到宏觀的力學(xué)描述，實現(xiàn)了降雨條件下滑坡的觸發(fā)模擬。本項目采用STEM TRAMM軟件進行計算，需要輸入的數(shù)據(jù)包括柵格數(shù)據(jù)(土層厚度、植被覆蓋率、地形數(shù)據(jù))和土體強度指標(biāo)。

(1)土體厚度。以研究區(qū)地形數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)，分別選取高程、坡度(以坡降計)、地形濕度指數(shù)(TWI)、距水系距離、方位、曲率，綜合計算土體厚度柵格圖。在研究區(qū)，選取因防火通道開挖而形成的開挖邊坡進行線路調(diào)查，共選取了162個點統(tǒng)計其土層深度，其中81個點用于柵格深度校準(zhǔn)，另81個點作為效果分析點位。

(2)植被覆蓋度。植被發(fā)育為基于Landsat 8影像提取的歸一化植被指數(shù)(NDVI)圖，利用公式NDVI=(IR-R)/(IR+R)計算NDVI值，其中IR和R分別為近紅外(NIR)和紅色波段。NDVI數(shù)據(jù)柵格大小為10 m×10 m。

(3)土體強度。參考文獻[24]及區(qū)內(nèi)21份勘察報告實驗結(jié)果，綜合選取用于數(shù)值計算的參數(shù)，主要有土體黏聚力、根系強度、內(nèi)摩擦角及初始含水率，對應(yīng)參數(shù)取值見表1。

表1 土體強度參數(shù)

2.2 降雨特征曲線

根據(jù)前文闡述，降雨條件需考慮10 a、20 a、50 a及100 a一遇降雨工況，但受制于降雨監(jiān)測數(shù)據(jù)精度及歷史的制約，雨量數(shù)據(jù)常獲取困難，無法根據(jù)要求選用不同頻率的日降雨量參與分析評價。因此，對降雨數(shù)據(jù)的分析主要需探索不同頻率下的降雨量及降雨歷程曲線。

關(guān)于降雨量，考慮到四川省內(nèi)實際情況，各研究區(qū)可直接參考《四川省暴雨統(tǒng)計參數(shù)圖集(2010年)》中所附的暴雨量等值線圖，其他地區(qū)可參照地方圖集或全國暴雨參數(shù)圖集。表2給出了研究區(qū)中米市—洛哈鎮(zhèn)1/6 h、1 h、6 h、24 h多年暴雨均值、變異系數(shù)等參數(shù)，并通過查詢皮爾遜Ⅲ型曲線得到不同頻率下的模比系數(shù)(Kp)，求得不同頻率下的雨強值。

表2 研究區(qū)不同頻率下暴雨統(tǒng)計參數(shù)及累積雨量

不同的降雨過程影響地質(zhì)災(zāi)害發(fā)育分布特征及災(zāi)害規(guī)模[25]，且上述數(shù)值模型的運用需將降雨過程納入模型的運算中，因此，評價研究區(qū)降雨的過程特征是進行不同區(qū)域、不同頻率降雨誘發(fā)地質(zhì)災(zāi)害研究的基礎(chǔ)。根據(jù)對地質(zhì)災(zāi)害發(fā)生與降雨過程的響應(yīng)過程分析，在該區(qū)域典型的3次災(zāi)害事件中，中期集中降雨誘發(fā)地質(zhì)災(zāi)害的特征最為明顯，且馬廠溝1985年7·1泥石流、1973年6·30群發(fā)地質(zhì)災(zāi)害及1983年7·1群發(fā)地質(zhì)災(zāi)害均為6 h左右。因此，針對該區(qū)域降雨特征，將降雨概化曲線設(shè)定為降雨歷時為6 h的中期集中降雨型事件。

本次評價主要選取累計降雨量(Ra)—持續(xù)時間(t)的關(guān)系式[26-27](Ra=-1.804 9t3+2.317 3t2+0.471 9t)來代表本區(qū)域誘發(fā)災(zāi)害事件的降雨雨型，累計雨量如表2所示。

此外，也可根據(jù)地區(qū)特征，選取國家雨量站數(shù)據(jù)(https://rp5.ru；http://data.cma.cn/)提取3 h或1 h降雨數(shù)據(jù)產(chǎn)品，經(jīng)過處理獲取月累計平均降雨量、日降雨量等數(shù)據(jù)，進而通過統(tǒng)計學(xué)方法獲取目標(biāo)區(qū)域的降雨歷程曲線。工作區(qū)降雨量在空間上存在顯著差異，需進一步考慮該指標(biāo)隨高程變化的空間差異性。由于本研究區(qū)高差較小，所以研究中未考慮降雨的空間差異性。

3 計算結(jié)果與效果分析

土體深度及降雨過程曲線是數(shù)值計算的基礎(chǔ)，以下分別對土體深度計算、土體深度數(shù)據(jù)分析及降雨過程曲線進行分析，并基于數(shù)值模擬計算淺層滑坡穩(wěn)定性指數(shù)。

3.1 土體深度計算與數(shù)據(jù)校準(zhǔn)

通過提取野外調(diào)查線路上的柵格值與對應(yīng)調(diào)查位置處的土體深度，對比校核獲取研究區(qū)土層深度的柵格圖(圖4(a))，并選取剩余81處調(diào)查點進行精度分析(圖4(c))。根據(jù)圖4(a)中深度分布柵格圖可知，研究區(qū)坡體土層深度范圍為0.05～10.18 m，采取調(diào)查路線上對應(yīng)點的柵格數(shù)據(jù)可獲取圖4(c)中灰色虛線所示的深度值。柵格圖中部分點土層深度波動較為劇烈，通過移動平均法繪制黑色趨勢線。選取圖4(b)中因消防通道開挖揭露良好的土層深度值，結(jié)合實測點深度標(biāo)準(zhǔn)誤差繪制散點圖(圖4(c))，結(jié)果表明研究區(qū)柵格圖能夠較好地反映研究區(qū)坡體處土層的深度分布特征。

(a) 研究區(qū)土層深度柵格圖 (b) 土層深度調(diào)查路線 (c) 土層深度及誤差分析

3.2 降雨歷程曲線計算

結(jié)合安寧河流域降雨誘發(fā)災(zāi)害及喜德縣歷次群發(fā)性地質(zhì)災(zāi)害降雨歷程可知[24, 28]，喜德縣群發(fā)性淺層滑坡均發(fā)生于前期長歷時降雨之后的短時強降雨期間。對比前述雨型分布特征的描述，喜德縣群發(fā)性地質(zhì)災(zāi)害發(fā)生條件與式(1)中雨型分布特征更為一致，因此，不同降雨頻率下的危險性分析均在此雨型條件下設(shè)定。根據(jù)研究區(qū)不同頻率下6 h累計降雨量所計算的降雨分布過程曲線如圖5所示。

圖5 研究區(qū)不同頻率降雨歷程模擬概化曲線

3.3 淺層滑坡穩(wěn)定性指數(shù)及分析方法

根據(jù)前文中所示的降雨歷程輸入數(shù)值計算模型中，結(jié)合易發(fā)性評價中地形(DEM)、覆蓋層厚度、歸一化植被指數(shù)(NDVI)及巖土體力學(xué)參數(shù)(表1)進行降雨條件下的穩(wěn)定性計算，將計算結(jié)果中基于柵格的穩(wěn)定性指數(shù)穩(wěn)定(0)—不穩(wěn)定(1)作為不同降雨頻率下的危險性評價值，評價結(jié)果如圖6所示。

根據(jù)圖7所示的評價結(jié)果統(tǒng)計可知，研究區(qū)在10 a一遇降雨頻率下穩(wěn)定性指數(shù)的平均值為0.42，標(biāo)準(zhǔn)差為0.32； 20 a一遇降雨頻率下穩(wěn)定性指數(shù)的平均值為0.44，標(biāo)準(zhǔn)差為0.33； 50 a一遇降雨頻率下穩(wěn)定性指數(shù)的平均值為0.58，標(biāo)準(zhǔn)差為0.37； 100 a一遇降雨頻率下穩(wěn)定性指數(shù)的平均值為0.58，標(biāo)準(zhǔn)差為0.37。可見，在圖5所示的降雨頻率下，隨著降雨強度及降雨量的增加，不穩(wěn)定指數(shù)的平均值也在增加，表明不穩(wěn)定坡體不斷增加。但隨著降雨量進一步增加到100 a一遇時，穩(wěn)定性指數(shù)平均值及標(biāo)準(zhǔn)偏差均無變化，同時，從圖7中也可以看出評價圖的指數(shù)均無變化。這說明研究區(qū)在50 a一遇降雨條件下即可導(dǎo)致區(qū)域不穩(wěn)定坡體達到較高含水量出現(xiàn)垮塌、溜滑。因此，50 a一遇降雨值可以作為區(qū)域淺層滑坡、泥石流及流域山洪的參考閾值。

圖7 不同降雨條件下穩(wěn)定性指數(shù)柵格值頻率統(tǒng)計直方圖

4 討論

上述基于不同頻率降雨的淺層滑坡危險性評價是針對降雨統(tǒng)計的預(yù)測值，其準(zhǔn)確性可以通過研究區(qū)典型的山洪泥石流事件后淺層滑坡分布位置與對應(yīng)降雨條件下穩(wěn)定性指數(shù)計算進行有效性評價。以下以2012年“8·31”洪水泥石流事件為實例探討上述方法的準(zhǔn)確性。

4.1 “8·31”山洪泥石流事件危險性評價

解明禮等[29]通過收集易發(fā)性相關(guān)文獻中具有代表性的200條成果進行統(tǒng)計，分析不同分級方法的等級劃分準(zhǔn)確性。由圖8可知，基于數(shù)值模擬獲取的不同降雨條件下穩(wěn)定性指數(shù)具有聚類的特征，運用K值聚類分析方法進行分級也取得了較好的效果(表3、圖8)。

圖8 基于K值聚類分析的8·31山洪泥石流事件危險性分級

表3 K值聚類分析

由淺層滑坡解譯的數(shù)據(jù)來源可知，本文淺層滑坡主要是在2013年3月的光學(xué)影像基礎(chǔ)上通過目視解譯獲取，基本反映了“8·31”孫水河流域山洪泥石流事件中淺層滑坡的發(fā)育特征。因此，將2012年8月31日降雨及其形成的山洪泥石流事件作為實例，利用其降雨曲線(圖9(a))進行計算，獲取如圖9(b)所示的此次降雨條件下的地質(zhì)災(zāi)害危險性評價結(jié)果。

(a) 降雨歷時曲線

(b) 地質(zhì)災(zāi)害危險性評價

(c) 評價結(jié)果與淺層滑坡遙感解譯對比

4.2 精度分析

合理的危險性評價結(jié)果應(yīng)實現(xiàn)現(xiàn)有地質(zhì)災(zāi)害點分布于高危險區(qū)、極高危險區(qū)，而非滑坡災(zāi)害點主要位于低易發(fā)區(qū)。量化檢驗評價結(jié)果的方法主要有AUC值、Kappa系數(shù)、經(jīng)驗概率值、總體精度、混淆矩陣及歷史數(shù)據(jù)論證等方法。AUC值評價法在危險性評價中較常用，且依賴于2012年“8·31”山洪泥石流事件后的淺層滑坡解譯點作為樣本，具有較高的可信度。因此，本文選取AUC值探討“8·31”山洪泥石流事件中淺層滑坡點的評價準(zhǔn)確性。

分別提取滑坡解譯面文件對應(yīng)柵格值及非滑坡區(qū)隨機生成的測試點柵格值，圖10給出了滑坡驗證成功點的驗證曲線，基于“8·31”淺層滑坡災(zāi)害點的驗證曲線表明該事件中AUC值為0.65，具有較高的準(zhǔn)確性。由圖9(c)可知，受柵格數(shù)據(jù)的原始精度、光學(xué)影像精度、光學(xué)影像拉伸等因素的影響，較大區(qū)域數(shù)值模擬結(jié)果與實際位置對應(yīng)性受到限制。故通過數(shù)值計算表達了區(qū)域大部分災(zāi)害點發(fā)生位置已具有良好的適用性。

圖10 滑坡區(qū)穩(wěn)定性指數(shù)驗證

此外，研究區(qū)易發(fā)性評價也可采用斜坡單元進行表示，因此需將基于柵格的危險性評價圖運用ArcGIS中的分區(qū)統(tǒng)計進行分區(qū)?？筛鶕?jù)需要選取平均值、眾數(shù)、最大值等指標(biāo)，本文通過對比現(xiàn)有地質(zhì)災(zāi)害分布位置，采用眾數(shù)作為危險性分區(qū)的指標(biāo)，進行斜坡單元的危險性等級劃分。

5 結(jié)論

本文以涼山州喜德縣洛哈鎮(zhèn)砂泥巖地層大量發(fā)育的淺層滑坡為研究對象，探索了紅層地區(qū)淺層滑坡危險性評價方法。

(1)以高程、坡度、地形濕度指數(shù)(TWI)、距水系距離、方位、曲率等因子為基礎(chǔ)，結(jié)合野外調(diào)查取樣核準(zhǔn)計算的方法用于區(qū)域性土體深度計算，有良好的適應(yīng)性。

(2)區(qū)別于區(qū)域降雨統(tǒng)計的頻率計算，根據(jù)有效降雨及雨型統(tǒng)計方法可實現(xiàn)重點區(qū)不同頻率降雨時間曲線的計算或預(yù)測。

(3)基于LHT的土體漸進式破壞模型可用于由降雨入滲而誘發(fā)的松散土體滑坡，可較準(zhǔn)確地預(yù)測類似地區(qū)的滑坡危險性。