基于浅水流模型的广州市大坪泥石流预测与防治

郑 伟，尹 培 杰

(1.深圳高速工程顾问有限公司，广东 深圳 518000； 2.长安大学 公路学院，陕西 西安 710000)

0 引 言

中国国土面积广阔，地貌类型多种多样，但同时也使得各种地质灾害频发。泥石流作为一种山区常见的自然地质灾害，对道路交通、村庄和其他基础设施会造成极大危害[1-4]。夏季随着雨水增多，使得原本松散的物源在水动力的驱使下快速运移，沿途侵蚀铲刮沟道岸坡，形成泥石流并最终堆积在下游，对沿途造成不同程度的影响。及时准确预测泥石流的危害范围和危害程度，为相关部门提供及时的灾害预警具有重要的意义。

通过经验公式计算法或现场物理尺度模型等传统方法来定量评价泥石流的影响相对困难，且精确程度受人为判断影响较大[1，5-6]。上述方法在预测泥石流的过程中往往缺少动态流动过程，模型的尺寸效应使得结果和现场吻合度较低。随着计算机技术的高速发展，科研人员对复杂过程的求解能力显著提高。数值模拟方法因成本低廉且结论可靠，成为泥石流精确计算的一种重要方法[7]。多种数值模拟方法，例如浅水流模型[4]、平滑粒子流体动力学模型[8]、离散元模型[9]、两相流甚至多相流模拟[10]，都已应用于泥石流模拟。

浅水流模型因其高效性成为科研和工程领域应用最广泛的一类计算模型。刘铁冀等[1]利用Massflow模拟拦挡坝对泥石流沟的治理情况，认为其能较好地控制沟内不稳定物源的活动情况。潘青等[2]利用离散单元的方法，通过数值模拟手段分析了泥石流运动过程对拦挡结构的冲击破坏作用，取得了良好的效果。张奋翔等[3]利用FlO-2D软件进行数值模拟，发现在允许条件下，将泥石流流速和流深与其发生频率相结合，得到危险评价图可为灾害预测提供良好的建议。Han 等[4]基于浅水流计算开发了SFLOW软件对泥石流和滑坡进行数值模拟，对比FlO-2D其对冲击范围模拟有着较好的结果。丹羽论等[5]依据二维的浅水流方程，研究了泥石流沟谷中设置的拦挡坝(谷坊坝)对泥石流的控制效果，并开展了数值模拟研究。

不同降雨频率下泥石流的冲出范围是泥石流预测的重点，而拦挡坝不同高度和位置对泥石流的减灾效应仍待研究。广东省部分山区长期受到泥石流灾害困扰[11-14]，由此造成的经济损失达数十亿元。本文选取广州市中部山区大坪沟为对象，采用浅水流模型对流域内泥石流的冲出范围展开数值模拟，分析其运动过程并提出相应的治理方案。

1 研究区概况

大坪沟位于广州市从化区，属于广东省中部地区，沟口坐标为东经113°46′23.5″，北纬23°40′44.14″(见图1)。研究区植被较为丰富，多灌木。根据走访当地居民可知，沟内多次发生洪水和泥石流，最近一次为2014年5月份，从化区普降大雨，大坪村数座房屋倒塌，农田被冲毁，公路多处中断且沿途多处塌方。据统计(http：∥roll.sohu.com/20140524/n399975660.shtml)，此次从化区温泉镇和良口镇是暴雨的重灾区，死亡4人，失踪4人，约1.2万人受灾。派潭镇8 h内录得雨量291.6 mm，流溪河洪水已经超过50 a一遇标准。

图1 广州市大坪沟整体景貌及沟谷的坡度和坡向分布范围Fig.1 The overall landscape of Daping gully in Guangzhou City and the distribution of slope and aspect

大坪沟沟谷流域面积约6.27 km2，主沟长度为4.85 km。左侧发育两条主要的支沟，主沟外侧可见鸭洞水河道由此蜿蜒而过。整体上呈现西南高东北低的地势形态，最高的海拔高程为1 168 m，最低海拔为98 m，整体高差为1 050 m，属于低山山地地形。大坪沟的主沟整体走向为西北向，沟道陡缓相接，沟道平均坡度为10°左右，局部坡度可达20°。除沟道坡度外，大坪沟两侧岸坡的坡度多在10°～30°之间，两侧山地稍显陡峭，坡度高于30°，局部山坡坡度甚至可达60°，这种岸坡坡度为泥石流汇聚水能提供了有利条件。此外，现场调查结果显示大坪沟内松散物源主要包括坡面崩滑体、沟道往期泥石流堆积物和人工堆积体，这些物源呈松散堆积，密实度较低，在雨水的冲刷之下可能进入沟道参与泥石流活动。

研究区地处中国大陆的最南端，属于亚热带季风气候。降雨量较高且降水相对集中，每年4～9月的降水量占全年的85%以上，此时是泥石流的暴发高频期。研究区附近观测的多年的降水量平均值为1 696.5 mm，日最大降水量为423.9 mm (24 h)[12]。

2 浅水流模型

2.1 理论基础

二维浅水流方程被应用于描述水平尺度远大于垂直尺度的流体，例如对泥石流、洪水等现象的描述。由于三维不可压Navier-Stokes方程是非线性方程，求解过程比较复杂，可以通过假设所计算的流体在其垂直方向上满足静水压力并且速度统一，将上述三维方程降为二维深度上平均、质量及动量上平衡的方程形式，即二维浅水流方程。其运算过程相对简单且仍保持着良好的准确性，因此在本研究中，将二维浅水流方程作为泥石流运动模型的控制方程。其方程的整体形式为[5-6]

(1)

其中：

q=(h，qx，qy)T= (h，uh，vh)T

(2)

f=(uh，u2h+1/2gh2，uvh)T

(3)

g=(vh，uvh，v2h+1/2gh2)T

(4)

(5)

式中：q为变量失量；f为x方向上的通量矢量；g为y方向上的通量矢量；S为源相矢量；q为平面x、y方向的单宽流量；u，v为平面x，y方向相应的流速；h为流深；zb为沟道海拔高程。

泥石流流体运动特征和普通水体相差较大。运动过程中与边界的摩擦行为较为特殊。基于不同的侧重点提出了多种泥石流的摩擦模型，比如Coulomb摩擦模型、Voellmy摩擦模型、O’Brien 等提出的流变摩擦模型，所得结果也存在一定差异。综合比较几种模型，流变摩擦模型包含了泥石流流体的摩擦性质、黏度性质及固体颗粒接触能量损失等，已经被广泛用于处理泥石流流体的摩阻行为[6]，其表达式为

(6)

(7)

式中：τ为屈服应力；ρm为泥石流中的固体物质密度；K为阻力系数；β为泥石流的黏度；ntd为等效曼宁阻力系数。

2.2 数据获取和参数设置

以网络下载的SRTM DEM为地形数据，通过重采样将其精度转化为8.75 m。利用ArcGIS软件栅格裁剪工具获得研究区范围的栅格数据文件，并将其转换为ASCII码文件带入程序。泥石流入流点对于数值模拟较为关键，通过实地野外调查和遥感解译可以确定研究区泥石流启动点[3，6]。本次模拟选择起点为泥石流沟形成区的下游点。

在模拟过程中设置流量过程曲线对数值模拟进行流量和时间控制。现实中泥石流的过程曲线是复杂的且多具有脉动性，与泥石流过程呈阵发性有关。泥石流按动力来源分为暴雨型泥石流和冰川融雪型泥石流，后者多分布于高山冰川地带。本文研究区处于热带地区，属于暴雨型泥石流，通常认为此类泥石流的暴发频率和降水频率是相互对应吻合的，可将洪峰流量简化为三角形或五边形进行计算，认为是单峰值过程[14]。此次研究将其概化为五边形模型，通过计算一次泥石流暴发的几个关键点，以一次泥石流时间的1/3段作为分界点，以峰值流量的1/4和1/3作为在这两个时间节点下的泥石流流量[15]，绘制出泥石流暴发的流量过程线(见图2)。

图2 五边形法的泥石流过程曲线[14]Fig.2 Debris flow process curve by Pentagon method[14]

在泥石流计算过程中，体积浓度决定着泥石流的基本性质和参数，也影响着泥石流运动状态。体积浓度相对大小直接决定了泥石流的密度，从而影响了流体的黏性，控制着屈服应力和流变参数等水文、力学参数。对于泥石流的运动状态，一方面，泥石流流体中含有着大颗粒的固体，彼此之间的镶嵌作用形成骨架结构对泥石流运动速度有着重大影响，而且在运移过程中会呈现出不均匀的分布状态；另一方面，泥石流的破坏力、侵蚀性等也受到固体物质含量的影响。体积浓度是指固体物质体积与泥石流总体积的比值。O′Brien等[16]通过总结不同条件下的泥石流体积浓度，给出一个经验值为0.55～0.65。根据参考文献[3-4]，本文选取Cv=0.55为泥石流的体积浓度。

(8)

式中：Vs是所含固体体积，Vw是所含液体体积。

2.3 泥石流洪峰计算

数值模拟过程中，泥石流洪峰流量的确定将会影响最终结果。不同的流域特征触发泥石流的降雨阈值呈现出较大差异性，导致泥石流洪峰计算结果出现差异。因此，为了研究不同降雨条件下大坪沟泥石流的冲出规模及危害范围，本文设计了6种降水频率，分别为P=20%，10%，5%，2%，1%，0.5%。根据水文手册计算洪峰流量，其计算公式[5，16]为

(9)

本文依据暴雨洪峰流量，利用雨洪法求取泥石流相应的峰值流量，并考虑泥石流携带一定比例的泥沙物质，对暴雨洪峰进行修正。针对沟道形状对其修正，不同沟道形状对应泥石流的洪峰相差甚大，对其运移的速度和侵蚀速率都有影响。引入清水洪峰修正和沟道堵塞修正系数，最终的公式为[5，16]

Qc=(1+φ)QpDc

(10)

式中：Qc为在降水频率为P时，相应的泥石流峰值流量，m3/s；(1+φ)为清水洪峰修正系数，可以由体积浓度Cv计算得到；Dc为沟道的堵塞修正系数。

最终设计的6种不同降水频率下的泥石流水文曲线如图3所示。

图3 在不同降雨频率条件下的水文曲线Fig.3 Hydrological curves under different rainfall frequency conditions

3 计算结果

3.1 数值模拟结果

对泥石流的数值模拟结果用ArcGIS软件进行处理，可以更加直观地看到泥石流冲出范围和危害情况。图4为大坪沟泥石流流深的数值模拟最终结果。在降雨频率为20%，10%，5%时，泥石流运动路径主要为沿大坪沟沟道运移，堆积厚度变化不大，最大厚度分别为2.1，2.2，2.7 m，一般沿沟道的泥石流堆积仅在1.5 m左右。而在降雨频率小于5%时，也就是降雨重现期超过50 a一遇时，泥石流运移至外侧鸭洞水河道，在沟口形成堆积扇，堆积厚度明显增加，最大厚度增至5.3 m。当降雨频率降低至0.5%时，泥石流在沟口堆积长度近800 m，而且向下游方向堆积范围更广，符合泥石流的运动特征。考虑到大坪沟内暴发特大型泥石流时可能会堵塞鸭洞水河道，导致河水暂时性断流，该情况下下游沿途村庄和道路都受到巨大的威胁，建议相关部门建立泥石流的监测和预警预报系统。但由于泥石流的堆积厚度不大，数值模拟结果为5 m左右，在河水作用下应该可以自然下泄。

在流速分布图5中，可见泥石流沿沟道保持高速运移，最大流速为3.3～6.1 m/s，说明泥石流在沟道内动能较大，表现为侵蚀搬运，尤其是对沟岸两侧和沟道的侵蚀。但是运移至外侧鸭洞水河道流速明显降低，泥石流在此则以堆积作用为主。在降水频率为20%和10%时，泥石流的流速均为较低速，最大不超过4 m/s，与中国著名的泥石流案例相比，例如汶川震后泥石流[17]和蒋家沟泥石流[18]分别为11.9 m/s和8 m/s，相差较大，这主要和泥石流的性质和洪峰流量相关。该泥石流最大流速可达6.1 m/s，虽仍低于上述两值，但此时泥石流已经具备较高的破坏力和冲击力[19-21]。

泥石流的危险对象评价是泥石流灾害模拟的最终目的[22-24]。在降雨频率大于5%时，泥石流对大坪村的影响较小，泥石流主要沿沟道运移，运移2 km左右。由于大坪村房屋建设主要位于主沟右岸地势较高之处，此时泥石流的危害性不大，在适当防护的条件下，对经济建设和生命财产安全影响不大。当降雨频率取值降低，也就是高于50 a一遇的暴雨时，泥石流的危险性明显增加。泥石流沿途可能会冲毁居民的房屋，最终堆积在沟口对岸的公路上，导致当地的道路瘫痪交通受阻。而且泥石流堆积范围也明显增大，仅沟口的堆积范围就达到0.3 km2，直接影响到下游村庄。

在数值模拟过程中，由于设置的集水点位于大坪沟左侧小支沟上游，泥石流运移在支沟沟口出现明显的堆积现象，也就是图4和图5中，物质刚运移时形成的“鼓包状”堆积。由于支沟沟口地形开阔，适宜泥石流堆积。在实际泥石流发生过程中，降水在主沟形成区汇聚开始携带沟道内的物质进行运移。当其运动到下游沟道时，来自侧面支沟的汇流会增大泥石流动力，在其沟口堆积现象不明显，或呈小型的舌状堆积。

图5 6种不同降水频率下泥石流的流速分布Fig.5 Velocity distribution of debris flow under six different precipitation frequencies

3.2 危险性区划

以往的研究表明可以按照泥石流流深对评价结果进行危险级别划分[3]，本文按照0～1.0 m划分为低危险区，1.0～2.0 m为中危险区，2.0～4.0 m为高危险区，大于4.0 m为极高危险区。对处理的结果在ArcGIS中进行栅格统计分析，可以得出条形统计图(见图6)。

图6 泥石流危险区评价和危害指数(HI)Fig.6 Debris flow risk area assessment and hazardous index (HI)

由图6可知，在评价结果中，泥石流的大多数影响范围均为低危险区。当设计的降雨频率高于5%时，低危险区的占比都接近或超过80%，其余的中危险和高危险区均位于支沟沟口和部分沟道。这3种情况的模拟结果均未出现极高危险区，说明此时泥石流危险性较小。但当降雨频率低于5%时，即降雨重复期分别为50 a，100 a和200 a一遇时，评价结果中出现极高危险区，所占比例分别为2.54%，7.59%，14.81%，对应的低危险区的面积明显下降，由原来的80%左右下降至77.65%，50.63%，47.23%，并且中危险区和高危险区也是呈现增大的趋势。在本文中，将高危险区和极高危险区占总体泥石流危害范围的比例定义为泥石流危害指数(HI)，该指数大小可以用来表征泥石流对沿岸道路和居民危害的程度，即

(11)

对于模拟得到的泥石流影响范围，在设计的6个降水频率情况下，泥石流危害指数HI分别为3.71，3.10，3.94，6.57，26.38，30.68，对应图6中的黑色折线显示，呈逐渐递增趋势，说明泥石流的危害性越来越大。泥石流最终的分布范围主要为外侧的鸭洞水河道，尤其是降雨频率为0.5%时。

3.3 拦挡坝治理

通过上述数值模拟的结果可见，泥石流对沟口村庄和道路均有着不同程度的影响。在降雨频率为0.5%时，泥石流运移过程出现范围较广的极高危险区。以往研究表明，拦挡坝对泥石流的治理效果是明显的，因此，建议在大坪沟沟内修建拦挡坝以降低泥石流对沿岸道路及居民的威胁。

本文采用地形重构技术以研究拦挡坝对泥石流的阻拦作用。将大坪沟沟道地形进行修改，具体操作是将地形栅格高程数据转换为ASCII码，以文本文件形式即可在沟道内提高相应的高程数据来模拟修建的拦挡坝。当拦挡坝坝高设置为10 m，即将原沟道高程最小ASCII码增加10，其他值相应增大以保持“拦挡坝”顶部水平。然后进行数值模拟，所用的水文曲线完全根据降雨概率P=0.5%时所设置。拦挡坝设计位置和高度的不同会影响到最终的治理效果，本文在大坪沟内设置两个拦挡坝的位置(A点和B点)，此外设计3种拦挡坝高度(5，10 m和15 m)，以探究其对泥石流的治理效果。

如图7所示，最终的数值模拟结果表明，设置拦挡坝之后，泥石流最大流深主要分布在拦挡坝之上，形成一个小型的“水库”。由图7可见：A，B两点设置5 m的拦挡坝治理效果均不明显(图中红色曲线代表200 a一遇的泥石流范围)，基本和原始情况一致，如图4(f)和7(a)、7(d)所示。在A点设置15 m高的拦挡坝最终的治理效果最为明显，泥石流在沟口外侧的主要影响范围被减小至不足原先的50%，且拦挡坝以下的泥石流流深全部为1.3 m以下，属于中低危险范围内，对沿岸的基础设施和居民影响基本不大，达到预期治理效果。泥石流的最大流深受到拦挡坝高度影响，这主要是最深处集中于拦挡坝上游的位置。

利用3.2节提出的泥石流危害指数HI评价拦挡坝治理效果。在A点设置拦挡坝为5 m时，HI指数为27.72，设置10 m和15 m拦挡坝时，HI指数分别为15.81和11.39，较原本200 a一遇的HI指数(30.68)有明显下降；B点设置拦挡坝时，HI指数对应为24.72，14.25和8.01，因此，据HI指数可见B拦挡坝优于A拦挡坝效果。但是由图7可见，设置A拦挡坝后，其上游出现了大范围的极高危险区和高危险区，导致HI指数明显增大。在实际情况中，依据拦挡坝下游的危险性区划求得HI指数才是正确评价的指标，因此，提出“有效”危险区，公式(11)可改写为以下形式：

(12)

由上式计算得到在A点分别设置坝高为5，10 m和15 m拦挡坝时的有效HI指数为7.33，0，0。

对于大坪沟的设置拦挡坝的位置，本文在主沟A，B两点进行研究发现，3种不同高度的拦挡坝在A点的治理效果均优于B点的治理效果，推断这主要是地形因素所导致的。A点位于一支沟沟口，适合拦挡坝截流由主沟而来的大部分洪水和泥石流，如图7(a)～(c)所示。但B点地形开阔，修建拦挡坝的治理效果较差。因此，增加设计一组数值模拟以验证上述结论，如图8所示，在A、B两点设置两级拦挡坝对沟内泥石流进行治理，拦挡坝高度分别设置为10 m和15 m。最终结果显示泥石流经A阻拦之后，B拦挡坝能发挥的作用较为有限，基本与单独设置A拦挡坝相近，说明A处设置拦挡坝能够有效治理大坪沟的泥石流。

图7 对大坪沟进行不同方案拦挡坝治理效果对比Fig.7 Comparison of dam treatment effects of different schemes for Daping gully

图8 对大坪沟进行两级拦挡坝治理效果对比Fig.8 Comparison of the effect of two-stage barrage treatment on Daping gully

综上所述，根据泥石流冲出危害范围、HI指数和拦挡坝合理选址，对于大坪沟泥石流可选择在A点出建立15 m高的拦挡坝进行治理。

4 结论与分析

本文采用浅水流的数值计算方法对广州市大坪沟泥石流进行数值模拟计算，基于6个不同的降水频率特征，计算了相应的水文过程曲线，并提出8种拦挡坝方案进行灾害治理的对比，分析表明：

(1) 泥石流的流深分布图显示，泥石流对大坪沟沟口居民和道路设施的影响最为显著。设计降雨频率为0.5%时，最大流深为5.3 m，可能沿外侧鸭洞水河道堆积长约800 m的范围。但在降雨频率高于5%时，泥石流主要沿沟道运移堆积，对沿途影响甚微。

(2) 泥石流的最大流速为3.3～6.1 m/s，主要集中在沟道，此时泥石流表现为侵蚀作用。当泥石流运移至外侧河道时，流速明显降低，此时，泥石流表现为堆积作用。

(3) 按照最终的泥石流堆积厚度进行危险区划分，表明当降水频率设计值大于5%时，泥石流影响区域80%以上均为低危险区。但是当设计值高于5%时，出现极高危险区。

(4) 对泥石流沟设置8个拦挡坝的方案进行数值模拟，结果显示，当在A点设置15 m拦挡坝时，主要的泥石流堆积量被控制在拦挡坝之上，此时下游均为中低危险区，有效HI指数为0。