北京市門頭溝區達摩溝泥石流發育特征與危險性評價

聶 瓊，聶治豹，陳 劍，丁士君，吳賽兒，李 多，葛潤澤，陳瑞琛

(1.國家電網有限公司，北京 100031；2.中國電力科學研究院有限公司，北京 100192；3.中國地質大學(北京)工程技術學院，北京 100083)

0 引 言

泥石流是山區溝谷中常見的自然災害，嚴重威脅山區居民的生命財產安全[1]。北京山區的區域性泥石流活動相對活躍，其爆發受到降雨強度影響，且泥石流流域的重現周期尚不明確，對位于老泥石流溝谷內部的居民聚落存在重大安全隱患[2-3]。近年來，受地震頻率加劇與極端天氣的影響，泥石流災害日益加劇[4-5]，北京市門頭溝區清水鎮達摩溝歷史上曾發育多期泥石流災害，目前仍存在大量人工堆積物和崩塌堆積構成泥石流的潛在物源，存在嚴重的泥石流隱患。

數值方法常被用于模擬泥石流的運動與堆積過程，量化分析泥石流的流速與影響范圍等特征[6-9]。美國學者O’Brien等[10]于1988年提出了二維動力模型(FLO-2D)，目前已廣泛應用于泥石流災害的危險性評價研究中：黃勛等[11]利用FLO-2D數值模型，實現泥石流重要承災體的風險量化；Lin等[12]運用FLO-2D數值程序模擬了中國臺灣松河地區泥石流的流動條件，并按不同的危險程度生成風險分布圖，取得較好的結果；劉福臻等[13]通過FLO-2D軟件，以結底崗村泥石流為例，模擬泥石流的沖淤特征，確定泥石流強度并劃定危險性分區；李寶幸等[14]利用FLO-2D軟件模擬麥多溝在4種極端降雨工況下(10%、5%、2%和1%)泥石流爆發過程，對潛在威脅區進行危險性分區。上述研究在泥石流危險性評價方面取得了一定進展，但較少對研究區泥石流的物源分布及儲量開展詳細的調查。

本文通過對北京門頭溝區達摩溝泥石流開展詳細野外調查，在摸清泥石流物源類型、分布特征、物源儲量與流域育災環境的基礎上，計算泥石流動力學參數，基于FLO-2D對達摩溝在4種極端降雨條件下的泥石流運動學特征進行研究，依據泥石流流速、堆積深度、堆積范圍和沖出規模等模擬結果進行泥石流危險性分區，為達摩溝泥石流的預測預警與風險防控提供科學依據。

1 達摩溝流域概況

北京西北部門頭溝區達摩溝位于太行山東麓、北京西山百花山西北側，清水河中游南岸，北部與北京昌平區及河北省懷來縣鄰近，南部連接房山區，西部與河北省琢鹿縣、淶水縣接壤，地理位置介于東經115°37′—115°39′、北緯39°52′—40°56′之間，系清水河南岸的季節性一級支流(圖1)。

達摩溝流域為典型“V”形谷，溝域面積為3.96 km2，地表起伏度(最大相對高差)為647 m。主溝溝向325°，長2680 m，溝寬30～100 m，平均15 m，溝床平均縱坡降23%。溝谷兩岸較陡，坡度30°～60°，平均坡度40°，植被覆蓋率為70%。年平均降水量為472.9 mm，年平均降水天數為79天。年降水量的74.8%集中于夏季；而冬季降水僅占年降水量的2%。達摩溝內出露地層主要為侏羅系砂巖與粉砂巖。同時，溝谷內廣發分布礫石、巖屑等殘積物、坡積物等松散堆積物。

2 研究方法

本研究基于遙感解譯、野外調查、數值模擬對達摩溝流域進行了多維度的調查分析。地面數字高程模型(DEM)被用來定量評估地形特征與數值模擬。研究所采用的DEM和正射影像(DOM)均通過OSGB格式傾斜攝影三維模型轉化而來，處理后的DEM分辨率可達到1 m(圖2)。

數值模擬基于FLO-2D進行。FLO-2D基于非牛頓流體模型及中央有限差分數值方法來計算二維洪水、泥石流運動控制方程的集成模擬軟件[15]。此軟件在滿足以下假設的條件下能夠較為精準地模擬流體的流速、堆積深度與堆積范圍[16-18]：

(1)假設泥石流運動屬于淺水波模式；

(2)假設流體在差分時間間隔內為恒流；

(3)不考慮溝道侵蝕；

(4)假設流體壓力分布為靜水壓力分布；

(5)假設各參數在每個網格單元中一致；

(6)不考慮流動過程中的跳躍及震蕩情況；

(7)不考慮泥石流對工程結構的損毀作用。

連續方程與運動方程作為FLO-2D基本控制方程被使用，其中連續方程為體積質量守恒方程(公式(1))，運動方程為力平衡動量方程，如公式(2)與公式(3)所示：

(1)

(2)

(3)

式中：t為時間；u為x方向的平均流速；h為泥深；v為y方向的平均流速；Sfx、Sfy為摩擦坡降；i為有效降雨強度；Sox、Soy為溝床坡降；g為重力加速度。

流變方程采用了O’Brien等[10]設計的適用于高含砂水流、泥石流和泥流的流變模式：

(4)

C=ρml2+f(ρm,Cm)ds2

(5)

(6)

將公式(4)可改寫成坡度形式：

(7)

τy=α1eβ1Cv

(8)

η=α2eβ2Cv

(9)

式中：ρm為土石流體質量密度；ρs為粒間壓力；ds為沉積物大??；ai為經驗系數；C為慣性剪切應力系數；Cm為含沙量；Cs為泥沙顆粒的最大靜態體積濃度；Sy為屈服坡降；η為黏滯系數；Sv為黏滯坡降；γm為土石流體相對密度；Std為紊流-擴散坡降；K為層流阻滯系數；τy為屈服應力，單位為Pa；e為自然常數，取值約為2.72；n為曼寧系數。

3 達摩溝泥石流數值模擬

3.1 泥石流動力學參數計算

(1)曼寧系數n。

曼寧系數可反映地表粗糙情況。不同的地表特征對應不同的曼寧系數，但在實際模擬中，通常設置相同的曼寧系數n值，以達方便簡潔的目的。一般來說，曼寧系數與地表植被狀況有關，表現為覆蓋率越高，地表粗糙程度越大，曼寧系數越大。本文研究的達摩溝植被覆蓋率為90%，植被以灌木為主，查閱曼寧系數推薦值(表1)后，選定本研究區曼寧系數為0.28。

表1 曼寧系數n取值范圍參考信息Table 1 Reference information on the Manning coefficient (n)range

(2)固體重度r以及泥沙修正系數φ。

固體重度反映泥石流溝道中松散堆積物的多少，并對泥石流運動速度、堆積深度、堆積范圍都有一定影響，表現為在體積濃度一定的情況下，固體重度越大，泥石流運動所需的水動力越大。本次模擬研究的泥石流溝松散堆積物大多來源于“7·21”特大暴雨條件下崩落的巖土體、山體兩側人工采石堆積形成的人工堆積體、以及溝道外坡面上的殘坡積物。參照《災害防治工程勘察規范》中的表G-2(數量化評分與重度、1+φ關系對照表)，取固體重度r=1.495 t/m3，泥沙修正系數φ=0.467。

(3)層流阻滯系數K。

層流阻滯系數K反映泥石流流體在流動狀態下，地表所給予反向運動的阻力。地表的植被種類、高度、覆蓋率等都會對K值產生影響。大多數泥石流在流動過程中受到阻力作用，由于顆粒之間黏性不夠，顆粒流動速度不同而產生分層流動，而層流間的摩擦系數則用K表示。根據表2，將本研究區K取值為2400。

表2 層流阻滯系數取值范圍參考信息Table 2 Reference information on the range of laminar flow arrest coefficient

(4)流變參數。

流變參數的選擇會直接影響屈服應力和黏滯系數的大小，而黏滯系數和屈服應力又會影響到泥石流的受力、變形和流動特征，因此選擇合適的流變參數非常重要。泥石流的體積濃度上升會引起流體中的賓漢屈服應力上升，以及賓漢黏性系數強度的增加，從而導致泥石流流動性變差。它們之間的表達式為：

τy=α2eβ2Cv

(10)

η=α1eβ1Cv

(11)

式中：τy為屈服應力；η為賓漢黏滯系數；α1、α2、β1和β2為經驗系數。

根據FLO-2D使用手冊[16]與詹錢登等[19]的實驗結果，結合野外調查“達摩溝”大部分屬于自然土壤，因此?。害?=0.811，α2=0.00462，β1=13.72，β2=11.24。

(5)體積濃度Cv。

體積濃度反映泥石流中固體物質體積的多少，數值上等于固體物質體積除以泥石流總體積。泥石流是由固體、液體、氣體組成的復雜的多相體。不同形態的物質構成了復雜的動力學特征。因此，體積濃度對泥石流運動堆積特征有著重要的影響。本次模擬體積濃度計算采用以下公式：

(12)

其中，固體物質體積為野外調查得到的人工堆積、沖洪積、崩坡積等松散物物源動儲量。由于沖洪積堆積于溝口，判斷其不參與泥石流流動過程。因此，固體物質體積采用溝域內人工堆積、殘破積、崩坡積三種松散物物源總儲量，再根據喬建平等[20]提出的物源總儲量與物源動儲量轉換，計算參與泥石流運動過程的松散物源動儲量(即泥石流運動過程中的固體物質體積，轉換關系見表3)。計算結果顯示，達摩溝泥石流物源總儲量為44.712×104m3，固體物質體積為6.7068×104m3，體積濃度為0.51。

表3 物源總儲量與動儲量轉換關系Table 3 Relationship between total source reserve and dynamic reserve

3.2 集水點和匯水面積

泥石流的集水點即為泥石流災害發生的啟動點，因此集水點的選取應在松散物源大量聚集的地方[7]，此點的選取還要求水動力條件足夠充足。本次模擬選取達摩溝主溝支溝各一個集水點，集水點位置如圖3所示。匯水面積則選取包括主溝道線在內的高程最高地三角區，沿山脊線繪制。匯水面積為1.566 km2。

圖3 達摩溝集水點位置圖Fig.3 Location map of Damogou Gully water collection spots

本文選取10年一遇(頻率為10%)、20年一遇(5%)、50年一遇(2%)和100年一遇(1%)4種降雨頻率下的降雨強度進行模擬計算。采用雨洪法，計算泥石流暴雨洪峰流量，進而求出泥石流峰值流量，再乘以相對應的放大系數，得出泥石流運動過程中的再生成洪峰流量[21]。最終繪制出模擬所需的泥石流運動過程中通過集水點的流量過程線：

QC=(1+φ)Qp×Dc

(13)

Qp=ψFs

(14)

式中：QC為泥石流洪峰值流量(m3/s)；ψ為暴雨徑流系數；Qp為泥石流暴雨洪峰流量(m3/s)；F為截面匯水面積(km2)；φ為泥石流泥沙修正系數；s為小時雨強(mm/h)；Dc為泥石流堵塞系數。

根據相關規范，選取泥沙修正系數為0.467。根據達摩溝流域面積以及松散物儲量大小，選取阻塞系數為1.2。北京地區泥石流溝均處于節理裂隙發育的山地，因此暴雨徑流系數取0.35。在ArcGIS上測量的匯水面積以及小時雨強如表4所示。由于軟件在模擬泥石流發生的過程中無法模擬出泥石流對于溝道的側蝕作用，因此需要將清水流量線B與流量放大系數F相乘，得到最終的清水流量線，BF計算公式如下：

表4 達摩溝各等級降雨頻率下集水點流量Table 4 Water collection spot flow chart under each grade of rainfall frequency in Damogou Gully

BF=1/(1-Cv)

(15)

計算結果如表4所示。

4 結果分析

4.1 達摩溝物源分布與儲量

通過野外調查發現研究區內泥石流的物源靜儲量為44.7×104m3。其中，溝道松散物分布在溝道中、下游，物質成分及結構主要為黏土與原始棄渣，寬度30～80 m(平均40 m)，厚度3～5 m(平均4 m)，大石塊直徑為0.5～2.0 m，平均粒徑為0.005～0.010 m，方量為2.86 m3；崩滑塌物源分布在溝域下游的溝谷兩側階地上，高于河床15 m，物質成分及結構為松散碎石，大石塊直徑為0.5～0.8 m，平均粒徑為0.2 m，方量為2.1×104m3；人工棄渣以煤矸石為主，分布于溝谷中上游兩側較危險的區域，粒徑為0.05～0.30 m，方量為39.7×104m3。

4.2 泥石流流速特征

在10年一遇降雨頻率(10%)下，泥石流流速整體以0～2 m/s為主，最大流速可達6.08 m/s，位于下游流通區溝道與石梯第一個交點處。溝道內部泥石流流速基本大于2 m/s，而溝道東側的漫延擴展區流速減慢，流通區上半部分流速較下游大，上游以大于1 m/s為主，下游以小于1 m/s為主。在泥石流淤積的前緣部位，受道路影響，泥石流流速較大，而在植被生長區流速較小。

20年一遇(5%)、50年一遇(2%)、100年一遇(1%)三種降雨頻率下泥石流流速特征與10年一遇降雨頻率下大致相似。但隨著降雨量增大，整體流速均逐級增大，0～1 m/s區域比例逐漸減小，而大于1 m/s區域占比明顯增大。在不同降雨強度下，流速最大點位置不變，均為下游流通區溝道與石梯的第一個交點。最大流速依次增大為6.31 m/s、6.60 m/s、6.94 m/s(圖4)。

圖4 達摩溝泥石流隱患點不同降雨頻率下泥石流流速分布圖Fig.4 Debris flow velocity distribution maps of different rainfall frequencies at hidden debris flow spots in Damogou Gully

4.3 泥石流流動深度及堆積面積

在10年一遇降雨頻率(10%)下，模擬結果顯示泥石流最大流深為10.82 m，位于下游堆積區最前端。泥石流發生時流動深度隨時間逐漸增大，直到到達泥石流溝下游堆積扇區域時，泥石流流深達到最大值。整體而言，泥石流溝道內的流動深度以2～4 m為主，溝道東側的漫延擴展區流深相對較小，以0～1 m為主。

20年一遇(5%)、50年一遇(2%)、100年一遇(1%)三種降雨頻率下流深特征與10年一遇降雨頻率下大致相似。但隨著降雨增大，整體流深均逐級增大，0～1 m區域比例逐漸減小，1～4 m區域比例明顯增大。流深最大點位置不變，仍為下游堆積扇最前端，最大流速依次增大為14.21 m、16.25 m、17.20 m。泥石流影響范圍也依次增大(表5和圖5)。

表5 模擬計算結果統計Table 5 Statistics of simulation results

圖5 達摩溝泥石流隱患點不同降雨頻率下泥石流流深分布圖Fig.5 Distribution maps of debris flow depth under different rainfall frequencies of debris flow potential spots in Damogou Gully

5 達摩溝泥石流危險性評價

目前，國內外一直應用很廣泛的危險性分區方法是利用泥石流流速與流深之間的函數關系來判別及評定泥石流溝的影響強度(泥石流危險性判定見表6)。一般將泥石流溝分為高、中、低三個危險性區域，各區域特點見表7。

表6 泥石流危險性判定依據Table 6 Basis for debris flow risk determination

表7 各危險性分區特點Table 7 Characteristics of each hazard zone

根據模擬結果分析，隨著降雨頻率逐漸增大，研究區內低危險區所占比例逐漸減小，高危險區所占比例逐漸增大。低危險區域占該流域泥石流影響范圍的比例為50%左右，主要分布于主溝道與支溝道兩側斜坡處以及堆積區內植被較茂密區域；中危險區域占比在30%左右，主要分布在溝道兩側以及流通區及堆積區植被裸露或道路修建的區域，此類區域離溝道較近，且地表粗糙度低，因此危險性較植被茂密區域高。高危險區域占比為10%～20%，主要分布在溝道內部及下游堆積扇地帶，沿溝道線上游零星分布，主要波及裸地和居民點(表8和圖6)。

表8 危險性分區結果Table 8 Risk zoning result

圖6 達摩溝泥石流隱患點不同降雨強度下危險性區劃圖Fig.6 Debris flow risk zoning map of hidden danger spots under different rainfall intensities in Damogou Gully

6 結 論

本文查明了達摩溝物源儲量與分布情況，并針對10年一遇(降雨頻率為10%)、20年一遇(5%)、50年一遇(2%)、100年一遇(1%)(降雨強度分別為49 mm/h、58 mm/h、68 mm/h、80 mm/h)四種極端降雨頻率工況，基于FLO-2D分析了達摩溝泥石流的運動特征，模擬了泥石流隱患影響范圍并進行危險性分區評價，得出如下結論：

(1)達摩溝泥石流物源靜儲量為44.7×104m3。其中，溝道松散物分布在溝道中、下游，物質成分及結構主要為黏土與原始棄渣，寬度30～80 m(平均40 m)，厚度3～5 m(平均4 m)，平均粒徑為0.005～0.010 m，方量為2.86 m3；崩滑塌物源分布在溝域下游的溝谷兩側階地上，物質成分及結構為松散碎石，平均粒徑為0.2 m，方量為2.1×104m2；人工棄渣以煤矸石為主，分布于溝谷中上游兩側，粒徑為0.05～0.30 m，方量為39.7×104m3。

(2)數值模擬結果顯示，達摩溝泥石流運動流速隨暴雨強度的增大而明顯增加。其中，溝道內速度可達6.08～6.82 m/s。同時，泥石流堆積深度也呈現出隨時間逐漸增大的特征，在到達堆積區最低點時最大深度可達10.82～17.02 m。泥石流整體沖出方量可達2.89×105～3.86×105m3。

(3)達摩溝泥石流在模擬中涉及的危險區總面積可達11.4×104～12.5×104m2，危險區總面積隨雨強增加逐漸增大。其中，高危險區占比在13.1%～20.0%之間，主要分布在溝口堆積區及主、支溝匯流處。中危險區和低危險區分別占比為29.9%～31.2%和48.8%～57.0%。數值模擬結果表明，高危險區和中危險區比例均隨著小時雨量強度增加而增大。

(4)達摩溝泥石流作為泥石流災害隱患點，溝口堆積區的居民點和主溝與支溝匯流處受到威脅程度最大。宜在達摩溝流通區設置上游攔截、下游排導的防治措施，并在暴雨季節加強形成區的雨量監測預警。