粵北山區某泥石流風險的精細化評價研究

何宏偉，溫漢輝，何利平

（廣東省有色金屬地質局九四〇隊，廣東·清遠 511520）

廣東山地面積約占全省土地面積的三分之二，隨著城市建設的不斷推進，在極端天氣、強震和人類活動等多因子耦合作用下，泥石流災害爆發的頻率居高不下。基于此，建立基于動力過程的泥石流風險評估體系，定量識別潛在泥石流風險等級，是山區城市發展規劃的迫切需要。

危險性是指災害發生的可能性、規模和強度。目前泥石流危險性評估以逐漸由經驗公式法、物理實驗法和統計類比法，發展到以數值模擬為主[1]，如：基于FLO-2D對汶川地震后的都江堰龍溪河流域內的12條泥石流溝災害進行了數值反演分析，綜合得出了該地區的泥石流災害危險圖[2]；對不同降雨頻率下的泥石流進行數值模擬，并表明潰決/非潰決型泥石流的沖出量比值能達2.88倍[3]；構建了基于災害動力過程的風險評估與風險制圖方法，進一步提出可用于泥石流減災防災風險調控技術[4]。泥石流的易損性是指承災體可能遭受地質災害破壞的嚴重程度，有學者總結了國外對易損性的定義和評估方法，將泥石流災害的潛在受害者類型分為財產（經濟）和人口（生命）兩大基本類型，并結合中國西南地區的調查和統計數據，建立了擬合資產、人口和易損性的冪函數方程[5]；以舟曲泥石流為背景，揭示了受災體的破壞特征以及破壞模式，并得出泥石流的密度、深度、速度以及與受損建筑物的結構、材料之間的關系[6]；基于非線性回歸方法，以泥石流的流深及建筑物高度等參數推導出了阿爾卑斯山脈泥石流堆積區內的建筑物易損性的經驗性公式[7]；以清平鄉泥石流為研究對象，將泥石流堆積區的土地類型劃分為建筑、道路、農田并建立了其敏感性因子，在此基礎上，以經濟損失為主要的評價指標，建立了該沖積扇上泥石流受災體的易損性分布圖[8]；針對不同計算模型的不確定性，結合現場調查，構建出泥石流強度與房屋建筑的易損度曲線[9]。

本文通過開展不同降雨頻率下的泥石流數值分析，獲取泥石流的動力學參數（流速、流深、沖擊力等），分別分析泥石流易發區的危險性及受災體易損性，最后采用泥石流風險評價方法，對泥石流易發區的風險性進行定量化評價，為泥石流防災減災提供一定的借鑒意義。

1 研究區概況

研究區位于粵北山區某村后山，其地貌類型屬丘陵地貌，地勢整體西高東低，相對高差100 m左右。山體整體坡度較緩，下部20°～30°，山頂附近較陡40°～50°，山間多溝谷，坡腳為現狀村民房屋，地形起伏大，沖溝發育，地形地貌較為復雜，目前主要已形成兩條泥石流支溝N1、N2，總流域面積11810 m2（圖1）。

圖1 研究區平面(a)與剖面(b)圖Fig.1 Plan view (a) and profile section (b) of the study area

2020年6月8日凌晨3時20分，受暴雨影響，該村后山發生泥石流地質災害（圖2）。泥石流裹挾花崗巖塊石及樹木順坡面地勢低凹處奔流而下，堆積于坡腳村莊內，泥石流一次堆積方量約3000 m3，其規模屬于小型。此次泥石流事件沖毀2棟、損毀1棟泥磚房，直接經濟損失約80萬元，潛在威脅村莊26戶157人生命財產安全。

圖2 泥石流裹挾的塊石(a)與沖垮的房屋(b)Fig.2 Boulders carried (a) and collapsed houses (b) by the debris flow

2 基于動力學模型的泥石流強度分析

2.1 模型構建原理

采用深度積分連續介質力學方法分析上述泥石流運動過程[10]。對x、y、z方向上的質量和動量守恒方程可以表示為：

其中：ρ代表流體的密度；h是流體深度；u和v分別代表流體速度在x和y方向的分量；(τzx)b和(τzy)b分別代表zx和zy方向上的切應力；z是地形標高；kap是側向土壓力系數。

Voellmy模型廣泛應用于泥石流數值模擬。基礎摩擦應力τ由Voellmy模型計算如下：

其中：μ是摩擦系數，ξ是湍流系數，V是泥石流的平均速度。

公式（1）～（3）由Massflow軟件進行求解，其是基于深度積分連續介質力學原理和改進的MacCormack-TVD有限差分方法，該方法在保證計算精度的同時簡化了計算過程，極大縮減了計算時間與計算內存需求，所以被廣泛運用于地表流災害數值分析領域。除此之外Massflow采用Fortran語言進行編程，可以進行二次開發。用戶可以結合自身研究對象特點對本構模型進行優化，使模型更貼合研究對象[11]。

本文數值模擬的地形數據由無人機傾斜攝影測量獲取的高精度地形（分辨率0.05 m），數值模擬網格設置為2 m。泥石流危險評價所用的參數值根據不同重現期的泥石流災害進行設定。我們使用兩個步驟來確定泥石流的參數，第一步參考以往相關研究成果[12-13]，基底摩擦系數μ的參考取值范圍為0.1～0.3，湍流系數ζ對應的參考取值范圍為100～300 m/s2。第二步根據第一步所獲取的參數取值范圍，利用反演的方法確定具體的計算參數，見表1所示。

表1 數值模擬參數Table 1 Numerical simulation of parameters

2.2 動力學參數設定

泥石流重度采用現場配漿法測定，將漿體準確配置成泥石流運動時的狀態，再結合漿體體積計算漿體重度。計算公式如下：

式中：γc是泥石流重度（t/m3）；Gc是配置泥石流漿體重量（t）；V為泥石流漿體體積（m3）。根據現場配漿結果得到泥石流重度為1.9 t/m3。

采用雨洪法計算泥石流的峰值流量。該方法假定泥石流與降雨同時發生，且具有相同的重現期。泥石流峰值流量（QC）可用以下公式估算。

其中：D是堵塞系數，與溝道的淤積程度有關；φ是泥石流容重的修正系數，定義為：

其中：γc為泥石流的密度（kN/m3），γw為水的密度（10 kN/m3），γH為泥石流固體物質密度（26.5 kN/m3）。

洪峰流量（QB）根據《泥石流災害防治工程勘察規范》（T/GAGHP 006-2018），計算公式如下：

其中：QB(1%)為百年一遇的洪峰流量。F為流域面積(km2)，對于研究區，峰值流量(QB)和不同重現期的流量可由以下公式確定：QB(2%)=0.8QB(1%)，QB(5%)=0.5QB(1%)，QB(5%)=0.3QB(1%)。

3 泥石流災害的風險評價方法

3.1 危險性評價方法

泥石流對承災體的破壞能力主要體現在兩個方面：淤埋破壞能力以及沖擊破壞能力。基于有關評價技術方案[14]，提出一套結合最大流體深度以及最大沖擊動量的雙因素評價模型來定量化評價泥石流的破壞能力（危險性）。這個模型將泥石流強度劃分為：極高、高、中、輕度危險區（表2）。

表2 泥石流危險性分區標準Table 2 Criteria for zoning of hazard of the debris flow

3.2 易損性評價方法

易損性指所研究的地質災害以一定強度發生時對潛在承災體可能造成的損失程度，常用0-1之間的數字來表征。此外，據相關研究結果[15]，表明承災體易損性可綜合采用社會易損性和物質易損性來表示。最后結合研究區的實際情況，采用社會易損性和物質易損性兩個指標建立人口和建筑易損性評價體系對泥石流易損性開展研究。

（1）易損性指標選取

由于研究區承災體結構較簡單，本文選取研究區內人口損失（社會易損性）和建筑物的經濟價值（物質易損性）作為指標因素進行易損性評價。

（2）易損性評價計算模型

為使兩種指標能置于同一數量級上進行比較，需利用公式對其進行轉換，結合單溝小流域泥石流易損性的計算公式，研究區泥石流的易損性評價計算模型如下[16]：

式中：V為易損性（0-1）；FV1為建筑物的經濟價值指標V1（萬元）的轉換函數（0-1）；FV2為人口損失指標V2（人/km2）的轉換函數（0-1）

式中：V1為建筑物的經濟價值指標；A為建筑物單層的面積（m2）；F為建筑物的層數；I為建筑物的造價（元/m2）

式中：V2為人口損失指標；a為≤13歲的兒童以及≥60歲的老人所占比例（%）；r為研究區5年內的人口自然增長率（‰）；D為每棟建筑物的人口密度（人/km2）。

（3）易損性指標分析

通過逐戶實地調查及參考當地建筑工程造價，研究區建筑物單價取值見表3。

表3 建筑單價取值表Table 3 Values of the construction unit price

利用每棟建筑物的人口數量及年齡結構，計算得到a的值，并結合該地區國家統計數據以及詢問當地村民獲取的數據，研究區村莊的人口自然增長率取值1.29%。人口密度則通過建筑物的總人口除以建筑物的面積得到。

3.3 風險評價方法

泥石流災害風險評價是將危險性評價結果與受災體易損性評價結果進行耦合得到，本文泥石流風險評價方法采用聯合國風險計算公式R（風險性）=H（危險性）×V（易損性）將泥石流的危險性及易損性結果進行耦合，計算出不同重現期泥石流風險度值，再通過分區標準將村莊的承災體劃分為不同風險區，具體明確地展示出風險水平。

參考有關研究[17-19]，首先將危險區分區圖低、中、高、極高危險性分別賦值為1、2、3、4，同時將易損性分區圖低、中、高、極高易損性分別賦值為1、2、3、4；再運用ArcGIS柵格計算器將危險性和易損性圖層進行賦值計算得到泥石流威脅區域內不同承災體的風險度值，最后通過風險計算公式計算風險度，計算標準見圖3，當R＞9時，屬于極高風險；當4＜R 9，屬于高風險；當2＜R 4，屬于中風險；當R＜2，屬于低風險。

圖3 泥石流定量風險評價標準Fig.3 Criteria for quantitative risk assessment of the debris flow

4 結果分析

4.1 模型驗證

村莊后山屬于典型的溝道型泥石流，根據現場調查訪問、泥石流監測及視頻資料分析，將啟動點設置在1#斷面（啟動點1）處和2#斷面（啟動點2）處，在啟動點處設值流量過程線。為簡單起見，將流量過程線簡化為概化五邊形模型進行計算，根據現場調查訪問將流量過程線持續時間T設置為300 s，泥石流最終堆積圖如圖4所示。

圖4 泥石流反演結果對比Fig.4 Comparison between results of back analyses of the debris flow

為驗證數值結果的準確性和合理性，結合相關研究成果[20]，使用Ω開展數值模型的驗證分析，Ω表示模擬結果和實際發生的泥石流堆積體的重疊程度，Ω=2表示完美的重疊，Ω= -2表示沒有重疊，其中Ω的定義如下：

式中：AX表示正確判斷面積，AY表示錯誤判斷面積，AZ表示缺失判斷面積，VX正確判斷的區域，A0表示現場調查面積，V0表示現場調查體積。

經計算Ω=1.44，精度為72%，表明此次泥石流數值計算結果與現場調查結果較為一致，證明該套數值模擬參數和模型能較好應用于泥石流災害數值計算。

4.2 泥石流危險性評價結果

采用已驗證的計算參數和數值模型開展不同重現期泥石流災害數值模擬。數值模擬結果表明：不同重現期泥石流均對村莊造成嚴重威脅，其主要區別在于泥石流運動過程中運動速度v、運動深度h及堆積區面積和堆積深度不同（圖5）。其中10年一遇最大堆積深度1.8 m，20年一遇泥石流堆積區堆積深度一般介于0.1～2 m，最大堆積深度2 m，主要威脅約6棟房屋；50年一遇泥石流堆積區堆積深度一般介于0.1～2.3 m，最大堆積深度2.3 m，主要威脅約19棟房屋；100年一遇泥石流堆積區堆積深度一般介于0.1～2.7 m，主要威脅約26棟房屋。根據數值模擬的結果，在不同重現期下大量泥石流主要沿溝道淤積，相對來說，村莊北部區域則處于相對安全位置。

圖5 不同重現期泥石流堆積深度圖Fig.5 Depth of debris flow accumulation with different recurrence intervals

各重現期泥石流的最大速度云圖如圖6所示，從圖中可以看出10年一遇、20年一遇、50年一遇、100年一遇，最大速度分別為5.5 m/s、5.7 m/s、6.0 m/s、6.7 m/s，隨著降雨量的增大，泥石流最大速度逐漸增大，且最大速度出現在泥石流中部位置。

圖6 不同重現期泥石流速度云圖Fig.6 Cloud map of debris flow velocity with different recurrence intervals

基于3.1節泥石流最大流深h和最大沖擊動量hv的雙因素分區標準，并使用ArcGIS軟件中柵格處理工具將泥石流危險性劃分為極高、高、中和低四種危險區，見圖7。

圖7 不同重現期泥石流危險性分區圖Fig.7 Zoning of debris flow hazards with different recurrence intervals

根據數值計算結果流深和動量綜合考慮泥石流淤埋破壞和沖擊破壞能力，定量評價不同重現期泥石流危險性，見表4所示。模擬結果表明：10年一遇危險性相對較小，20年一遇、50年一遇、100年一遇的泥石流災害均會對村莊建筑物、道路和農田造成較為嚴重的危害，其中20年一遇極高危險區和高危險區面積分別為1300 m2和3384 m2，中危險區和低危險區面積分別為5640 m2和10872 m2；50年一遇極高危險區和高危險區面積分別為2164 m2和4272 m2，中危險區和低危險區面積分別為9120 m2和9856 m2；100年一遇極高危險區和高危險區面積分別為3816 m2和6028 m2，中危險區和低危險區面積分別為6196 m2和10680 m2。

表4 不同重現期泥石流危險區統計表Table 4 Statistics of debris flow caution areas with different recurrence intervals

4.3 易損性評價結果

易損性評價結果計算主要有三步：首先根據現場調查的建筑物、人口等數據分別代入公式（12）和（13）得出建筑物的經濟價值指標V1和人口損失指標V2的值；再將V1和V2分別代入公式（10）和（11），得出建筑易損性指標的函數轉換值FV1和人口易損性指標的函數轉換值FV2，計算結果見表5。從表中可以看出人口易損性相比建筑易損性的數值較低，建筑易損性主要集中在0.4～0.6，人口易損性主要集中在0.1～0.3。

表5 易損性指標分布表Table 5 Distribution of vulnerability indicators

根據現場調查訪問，村莊人口分布結果分布見圖8a，建筑結構結果分布見圖8b；最后將FV1和FV2帶入公式（9）中，得到研究區承災體易損性綜合評價結果8c，通過GIS統計可知，極高易損性建筑面積為721 m2，高易損性建筑面積為3156 m2，中易損性建筑面積為1133 m2，低易損性建筑面積為221 m2。從結果可以看出易損性與建筑物結構類型及人口數量密切相關，結構越差，人口越多則易損性越高。

4.4 風險評價結果

通過上述泥石流風險評價方法對研究區風險性進行劃分，得到了不同重現期泥石流定量風險評價結果，見表6及圖9所示，具體明確地展示了不同承災體的風險水平。模擬結果表明：10年一遇泥石流風險區面積17004 m2，其中高風險區占比5.0%，中風險區占比9.3%，低風險區占比85.7%，主要威脅農田；20年一遇泥石流風險區面積21401 m2，其中極高風險區占比0.6%，分布建筑物1棟，高風險區占比7.1%，分布建筑物2棟，中風險區占比16.0%，分布建筑物2棟，低風險區占比76.3%，分布建筑1棟；50年一遇泥石流風險區面積為26112 m2，其中極高風險區占比4.1%，分布建筑物6棟，高風險區占比11.5%，分布建筑物7棟，中風險區占比14.1%，分布建筑物3棟，低風險區占比70.3%，分布建筑物3棟；100年一遇泥石流風險區面積為27754 m2，其中極高風險區占比5.6 %，分布建筑物9棟，高風險區占比18.1%，分布建筑物10棟，中風險區占比18.7%，分布建筑物5棟，低風險區占比57.6%，分布建筑物2棟。隨著降雨頻率的增加，受威脅的房屋逐漸增加，因此，建議完善溝內攔擋治理措施及加強泥石流監測預警工作。

表6 不同重現期泥石流風險區統計Table 6 Statistics of debris flow risk areas with different recurrence intervals

圖9 不同重現期泥石流風險分區圖Fig.9 Risk of the debris flow with different recurrence intervals

5 結論

本文以深度積分連續介質力學方法作為泥石流動力學過程的定量化獲取手段，建立耦合泥石流淤埋深度和沖擊動量的泥石流強度評價指標，考慮人口和建筑易損性構建泥石流風險評價體系，主要得出以下結論：

（1）基于N-S方程構建泥石流運動演進物理模型，使用該動力學模型對此次泥石流反演結果與現場調查結果較為一致，在此基礎上對不同重現期（10年一遇、20年一遇、50年一遇和100年一遇）泥石流災害進行數值計算。結果表明：隨著降雨強度的增加，危險面積由17004 m2增加至26720 m2，其中百年一遇極高危險區、高危險區、中危險區和低危險區面積分別占比14.3%、22.5%、23.2%和40.0%。

（2）根據統計的建筑物和人口數據，研究區一層砌體結構建筑物占92%，每棟建筑物的人口數量在3～6人左右，以幼老年人為主。易損性主要受建筑物結構類型及人口數量影響，由綜合考慮建筑物結構類型和人口的易損性評價模型得極高易損性房屋5棟，高易損性房屋15棟，中易損性房屋4棟，低易損性房屋2棟。

（3）根據風險評價的結果，10年一遇泥石流主要威脅對象為農田；而20年一遇、50年一遇、100年一遇泥石流則會對建筑物造成威脅，其中20年一遇泥石流位于極高風險區建筑物1棟，高風險區建筑物2棟，中風險區建筑物2棟，低風險區建筑物1棟；50年一遇泥石流位于極高風險區建筑物6棟，高風險區建筑物7棟，中風險區建筑物3棟，低風險區建筑物3棟；100年一遇泥石流位于極高風險區建筑物9棟，高風險區建筑物10棟，中風險區建筑物5棟，低風險區建筑物2棟。建議完善溝內攔擋治理措施及加強泥石流監測預警工作。