北京市门头沟区达摩沟泥石流发育特征与危险性评价

聂 琼，聂治豹，陈 剑，丁士君，吴赛儿，李 多，葛润泽，陈瑞琛

(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]的实验结果，结合野外调查“达摩沟”大部分属于自然土壤，因此取：α1=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)达摩沟泥石流作为泥石流灾害隐患点，沟口堆积区的居民点和主沟与支沟汇流处受到威胁程度最大。宜在达摩沟流通区设置上游拦截、下游排导的防治措施，并在暴雨季节加强形成区的雨量监测预警。