北京市門頭溝區(qū)達摩溝泥石流發(fā)育特征與危險性評價

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

(1.國家電網(wǎng)有限公司，北京 100031；2.中國電力科學(xué)研究院有限公司，北京 100192；3.中國地質(zhì)大學(xué)(北京)工程技術(shù)學(xué)院，北京 100083)

0 引 言

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

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

本文通過對北京門頭溝區(qū)達摩溝泥石流開展詳細野外調(diào)查，在摸清泥石流物源類型、分布特征、物源儲量與流域育災(zāi)環(huán)境的基礎(chǔ)上，計算泥石流動力學(xué)參數(shù)，基于FLO-2D對達摩溝在4種極端降雨條件下的泥石流運動學(xué)特征進行研究，依據(jù)泥石流流速、堆積深度、堆積范圍和沖出規(guī)模等模擬結(jié)果進行泥石流危險性分區(qū)，為達摩溝泥石流的預(yù)測預(yù)警與風(fēng)險防控提供科學(xué)依據(jù)。

1 達摩溝流域概況

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

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

2 研究方法

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

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

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

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

(3)不考慮溝道侵蝕；

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

(5)假設(shè)各參數(shù)在每個網(wǎng)格單元中一致；

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

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

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

(1)

(2)

(3)

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

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

(4)

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

(5)

(6)

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

(7)

τy=α1eβ1Cv

(8)

η=α2eβ2Cv

(9)

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

3 達摩溝泥石流數(shù)值模擬

3.1 泥石流動力學(xué)參數(shù)計算

(1)曼寧系數(shù)n。

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

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

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

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

(3)層流阻滯系數(shù)K。

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

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

(4)流變參數(shù)。

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

τy=α2eβ2Cv

(10)

η=α1eβ1Cv

(11)

式中：τy為屈服應(yīng)力；η為賓漢黏滯系數(shù)；α1、α2、β1和β2為經(jīng)驗系數(shù)。

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

(5)體積濃度Cv。

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

(12)

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

表3 物源總儲量與動儲量轉(zhuǎn)換關(guān)系Table 3 Relationship between total source reserve and dynamic reserve

3.2 集水點和匯水面積

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

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

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

QC=(1+φ)Qp×Dc

(13)

Qp=ψFs

(14)

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

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

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

BF=1/(1-Cv)

(15)

計算結(jié)果如表4所示。

4 結(jié)果分析

4.1 達摩溝物源分布與儲量

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

4.2 泥石流流速特征

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

20年一遇(5%)、50年一遇(2%)、100年一遇(1%)三種降雨頻率下泥石流流速特征與10年一遇降雨頻率下大致相似。但隨著降雨量增大，整體流速均逐級增大，0～1 m/s區(qū)域比例逐漸減小，而大于1 m/s區(qū)域占比明顯增大。在不同降雨強度下，流速最大點位置不變，均為下游流通區(qū)溝道與石梯的第一個交點。最大流速依次增大為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%)下，模擬結(jié)果顯示泥石流最大流深為10.82 m，位于下游堆積區(qū)最前端。泥石流發(fā)生時流動深度隨時間逐漸增大，直到到達泥石流溝下游堆積扇區(qū)域時，泥石流流深達到最大值。整體而言，泥石流溝道內(nèi)的流動深度以2～4 m為主，溝道東側(cè)的漫延擴展區(qū)流深相對較小，以0～1 m為主。

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

表5 模擬計算結(jié)果統(tǒng)計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 達摩溝泥石流危險性評價

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

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

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

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

表8 危險性分區(qū)結(jié)果Table 8 Risk zoning result

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

6 結(jié) 論

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

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

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

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

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