非匀质泥石流运动特征与峰值流量关系研究*

王高峰 田运涛 高幼龙 乐琪浪 杨 强 李瑞冬 叶振南 邓 兵 郭 宁 孙秀娟

(①中国地质调查局水文地质环境地质调查中心，保定 071051，中国)(②自然资源部三峡库区地质灾害监测预警野外科学观测基地，重庆 404100，中国)(③甘肃省地质环境监测院，兰州 730050，中国)

0 引 言

泥石流是山区广泛发生的自然灾害，尤其在我国西部山区泥石流灾害极为严重，典型的2010年8月7日舟曲三眼峪特大泥石流灾害、汶川震区多次强降雨群发性泥石流灾害及2018年10月17日西藏雅鲁藏布江色东普沟堵江事件。为防灾减灾，开展非匀质泥石流运动特征模型研究，分析其与泥石流流量的关系具有十分重要的现实意义。

泥位和流速作为泥石流三大运动要素(泥位、流速、冲击力)的一部分，对泥石流防治工程规划设计、监测预警预报和灾害评估尤为关键，随着泥石流监测预警手段的智能化、集成化，对可监测的泥石流能直接获取其基础数据。如我国科研人员通过观测分别获取了陇南白龙江流域泥湾沟、火烧沟等和小江流域蒋家沟泥石流泥位深度和流速的原始资料。当泥石流运动特征值超过一定值后，泥石流爆发过程无法监测，其特征值只能靠一些公式来计算获取。例如，可通过流体模型FLO-2D得出不同规模的泥石流流速(王高峰等，2019)；采用测量野外沟道泥痕数据推算泥石流流速(Crosta et al.，2001)；在建立科学分界粒径模型的基础上，构建了我国非匀质泥石流流速计算公式(舒安平等，2012)。长期以来，我国无论在泥石流防灾减灾工程设计还是理论研究方面，泥石流流速都是基于均匀条件下的曼宁公式来确定，并在此基础上建立了许多泥石流流速计算模式(表1)，得到较好的应用效果。但对于真实地形上发生的泥石流，均匀流假设很难成立，然而此公式是针对某一泥石流沟道所得出的统计经验模型，但大多数的沟道都缺乏系列观察数据，限制了曼宁公式的普适性(崔鹏等，2016)。此外，该公式针对宽浅式沟道泥石流流速计算效果较为可靠，对窄深式沟道的情形还有待进一步研究。当泥石流经过弯道时表现出极为强烈的爬高、冲淤等特征，基于此现象提出了弯道超高法，上述公式主要适合于泥石流弯道流体角速度相等的情况，但现实中弯道不同位置的角速度很难确定。

表1 国内常用的黏性泥石流运动特征值计算公式

通常泥位主要与泥石流流量、沟床断面形态等有关，具有“规模-频率关系”的分布特征(李泳等，2004)。目前有关泥位的研究多集中在监测预警方面，但针对泥石流泥位特征值定量评价的研究成果却鲜有报道。一些专家学者认为高黏性和黏性泥石流的堆积泥位深度与泥石流切应力、沟床坡度及容重有关，当沟床表面经铺床“润滑”后，该公式不再适用(吴积善等，1990)；当假设非匀质混合相泥石流整体平均流速与某一断面流速相等时，在断面宽度系数确定后，可计算得到一定流量条件下的泥位特征值(舒安平等，2012)。国外学者认为针对大尺度泥石流的研究，试验得出的公式不能完全得到应用(Prochaska et al.，2008)。一些研究者从泥石流流经弯道受力平衡、流体静力学原理等建立了弯道超高公式(周必凡等，1991；陈宁生等，2009)；也有通过修正前人公式系数及考虑爬高效应，并通过实例验证效果较为合理(Hungr et al.，1984；赵晋恒等，2015)。由于泥石流输移、冲刷或堵塞等运动特征的不均匀性，导致现有的弯道超高模型计算结果仍存误差，但一些弯道超高模型与泥石流残留泥痕有很大的相关性，可经过修正换算确定泥石流泥位特征值(表1)。上述问题的存在促使泥石流运动特征值新方法、新模型的进一步探索和试验。

非匀质泥石流是一种典型的固液两相流，通常具有容重大、颗粒级配宽等特征，其表面泥位深度和流速沿沟道纵横剖面是不连续分布的(舒安平等，2013)。实际上泥石流经弯道流动时角速度时变性强，同时由于泥石流运动阻力的存在，泥石流沟床粗糙率难以合理地确定。此外，泥石流流量在充分利用现有小流域暴雨雨洪计算模型的基础上，针对不同地区通过合理修正流量参数，可得到较为可靠的不同频率下泥石流流量预测值。因此本文利用陇南泥石流野外观测点武都区泥湾沟泥石流运动特征值观测资料进行整理分析，基于泥石流流量与泥位深度和流速具有较好的幂相关性，研究非匀质泥石流运动特征与峰值流量关系，建立了泥石流运动特征值修正预测模型, 并以陇南市礼县下胡杨沟为例进行模拟分析，计算不同频率下泥石流流量、流速及泥位的特征值，分析泥石流运动变化规律，为有效防治泥石流灾害及监测预警体系建设提供技术支撑。

1 非匀质泥石流运动特征值预测模型构建及验证

1.1 模型构建

泥石流监测断面为“V”型，位于流通区下游海拔高程1155im处，距离出山口920im，断面宽4.6im，高3im(图1)，该断面处沟道顺直，沟岸稳定。布设超声波泥位和流速监测仪，采用长期观测、定期观测和山洪泥石流应急观测相结合，重点监测泥石流泥位和流速等要素。

图1 泥湾沟流域概况及观测断面位置关系图

泥石流流量资料的获取采用监测断面处实测泥位深度和流速，利用泥石流流量数学模型计算QC=vA(QC为泥石流流量(m3·s-1)；v为泥石流断面平均流速(m·s-1)；A为泥石流过流面积(m2)获得。由于部分场次泥石流泥位或流速过小，未到达观测断面，因而没有泥石流泥位或流速记录，另外，当泥石流洪峰流量达到一定值时，由于泥石流瞬时流速过大，监测仪器存在采集时间间隔问题，导致泥位或流速获取滞后。

因此，筛选两者数据准确且相对应的1998～2008年间爆发的174场有完整泥位和流速观测记录资料的泥石流进行分析，在此基础上采用统计方法构建了泥石流流量(QC)与泥位(h)、流速(v)的数学表达式：

(1)

(2)

式中：h为泥石流沟道泥位深度(m)；v为泥石流流速(m·s-1)；QC为泥石流流量(m3·s-1)。泥石流泥位深度预测模型的复相关系数R=0.9132，泥石流流速预测模型的复相关系数R=0.8247，对应的幂函数相关性高(图2)。上述构建的泥石流运动要素特征值预测的统计模型经过自检验误差分析，均满足精度要求。

图2 泥石流流量与泥位、流速的关系

1.2 模型验证与应用

为检验上述构建的非匀质泥石流运动特征值预测模型的合理性，本文选取了我国西部白龙江流域、汶川震区、藏东南山区、小江流域等泥石流高发区域，具有典型性的14条泥石流沟共20次泥石流事件作为预测模型验证应用的研究对象(表2)。其中泥石流运动特征值以及流量值是通过实地调查和研究文献资料获取，对比分析泥石流泥位深度、流速计算值与实际值，可得到不同地区泥石流运动特征值的修正系数，并据此构建更为可靠的预测模型。

表2 验证区14条泥石流沟基本参数及预测模型误差统计表

根据表2所示验证结果表明，白龙江流域和汶川地震区泥石流泥位深度计算值与实际值的绝对误差在0.49～16.90之间；泥石流流速计算值与实际值的绝对误差在1.09～24.86之间。模型误差范围在该区较为合理，而且该模型的相对误差正负均有，说明该模型对白龙江流域和汶川地震区泥石流运动特征值的预测是均衡的，满足预测要求。另外，小江流域泥石流流速计算值与实际值的绝对误差在1.72～21.93之间，说明模型满足预测要求。但在藏东南地区和小江流域泥石流泥位深度计算值与实际值的绝对误差均大于40，而藏东南地区泥石流流速计算值与实际值的绝对误差在28.27～58.49之间，说明预测模型在该区不合理。

由图3和表2可知原因为：第1，小江流域泥石流属于阵性运动的高黏性泥石流，在泥石流运动时，龙头的速度显著增强，随着固体物质的不断参与，泥石流浓度和流量处于逐渐增加的趋势，则会出现连续淤积的现象。当流量增大至某一临界值时或黏性泥石流向稀性泥石流过渡以及阵性流转变为连续运动时，沟道的淤积层则会遭到破坏，其运动速度开始逐渐降低，但总体平均速度变化太大。再者，所选取蒋家沟泥位检验数据，由于测量断面宽度较大，泥石流运动存在向两侧扩散的二维流，且在沟床表面流动具有时变性即表现为床面流(吴积善等，1990)，致使观测值与计算值泥位特征值缩减比例在0.5～0.7之间。第2，藏东南地区泥石流多属冰川诱发型，具有瞬间暴发、运动速度快、运动距离远及冲出规模大等特点。泥石流形成的动力过程为：在流域上游降雪、降水产生大规模冰崩、滑坡，堵断沟道形成堆积坝，随着汇流过程变形与集中，堆积坝溃决，伴随着铲刮、搬运沟床冰碛物，使流量瞬间增大；之后泥石流经流域中游峡谷段阻挡，长距离输移到达下游宽谷段的逐渐淤积展平，伴随着沿程沟道输移和支沟洪峰汇入，使泥石流浓度发生改变。但是在同等水力条件下，此类泥石流的运动速度一般小于水流，但超常的运动速度是由于泥石流发生堵塞溃决作用，在泥石流活动中，随着固体物质的不断参与，改变了沟床的水动力条件，致使流速增大，而实际中泥石流泥位特征值比计算值缩减比例在0.5左右、流速特征值缩减比例在0.6～0.8之间。

图3 泥石流运动特征值计算值与实际值的关系

为使预测结果更可靠，对藏东南地区和小江流域泥石流运动特征值预测模型进行了修正。修正后的泥石流运动特征参数与峰值流量计算公式见表3。由修正后的预测模型计算结果可知(表4)，泥位深度计算值与实际值的绝对误差在0.49～22.49之间；流速计算值与实际值的绝对误差在1.09～26.22之间，上述验证结果表明修正后的泥石流运动特征值与峰值流量计算模型能够满足预测精度要求。

表3 不同地区泥石流运动特征值修正公式

表4 验证区14条泥石流沟修正后预测模型误差统计表

2 数值模拟计算实例——以下胡杨沟为例

2.1 流域概况

为了定量评价不同流量条件下泥石流运动变化规律，本文以陇南礼县下胡杨沟为例进行模拟分析研究。下胡杨沟位于礼县县城西南侧，西汉水左岸，其流域面积0.61ikm2，区内相对高差422im，主沟道长1.1ikm，沟床平均纵比降为383.6‰，陡峻的地形条件为泥石流输移提供了势能保障。该区域属于典型的温带大陆性季风气候区，多年平均降雨量为489.9imm，历史年平均最大降雨量为965.9imm，日最大降雨量达116.3imm，每年7～9月份暴雨频繁，强度大且历时短，为泥石流的形成与运动提供了充足的降水条件。流域出露地层为泥盆系白云质灰岩夹炭质板岩、千枚岩，第四系崩坡积物及冲洪积堆积物。EW向礼县—罗家堡断裂与NWW向彭家崖—风台山弧形断裂在此交汇穿过，且局部次级断裂和小型褶皱集中发育，及地震活动频繁，导致流域松散固体物源总量达99.92×104im3，大量松散物质堆积于沟道，为泥石流发生提供了丰富的物源条件，同时也增强了泥石流沟床冲刷能力。在流通区上游伴有沟道堵塞体和卡口的存在，增加了沟道堵塞体瞬时溃决发生的可能性。

下胡杨沟是一条高频泥石流沟，根据历史泥石流资料整理分析，下胡杨沟于1964年、1984年、1998年及2013年均暴发过泥石流。结合流通区不同断面处堆积物颗粒级配分析，下胡杨沟泥石流是陇南山地典型的非匀质泥石流，其容重在1.83～2.15it·m-3之间，属于高黏性泥石流。现场调查表明，该泥石流动力过程为：软弱千枚岩风化形成的浅层松散堆积物启动→沿程沟道侵蚀→堵塞体瞬时溃决流量放大→更加强烈的弯道冲淤→停淤堆积或堵塞河道，属典型的沟道侵蚀搬运-溃决混合型泥石流。

2.2 数值模拟验证

FLO-2D模型主要参数获取根据0.5im分辨率DEM数据、5组不同断面处颗粒级配组成及其他流域特征参数。结合泥石流动力特征和甘肃省小流域暴雨洪峰模型计算得出的不同容重1.83it·m-3、1.97it·m-3条件下泥石流流量。实际调查中下胡杨沟道严重堵塞，浅层滑坡堰塞体溃决后形成泥石流会产生一定的放大效应，故输入FLO-2D模型的泥石流流量为计算得出的泥石流流量乘以体积膨胀系数。最后将泥石流参数和泥石流流量(表5)和(表6)输入FLO-2D模型，整个过程未有人为干预，计算结果较为可靠真实。

表5 FLO-2D数值模拟参数及精度分析

表6 下胡杨沟泥石流FLO-2D数值模拟运动特征值对比分析

根据访问和实测断面泥位深度可按照表3公式反推下胡杨沟不同降水条件下泥石流流量，经计算得出下胡杨沟2%频率和5%频率降水条件下泥石流流量分别为103.60im3·s-1、51.63im3·s-1。在50年一遇情景下泥石流流量差异较大，采用本文预测值是小流域暴雨洪峰模型计算结果的1.23倍，绝对误差为23.16%，但是两种方法求得20年一遇情景下泥石流流量基本一致。从整体验证结果看本文通过修正后的泥石流运动特征值预测模型具有一定的可靠性和应用性。

2.3 不同规模泥石流运动特征差异性分析

取下胡杨沟分别暴发50年一遇(2%频率)和20年一遇(5%频率)泥石流泥位和流速动态过程数据，分析泥石流运动规律。根据数值模拟结果绘制50年一遇(2%频率)泥石流运动特征曲线、20年一遇(5%频率)运动特征曲线(图4)，可见下胡杨沟泥石流不同频率条件的规模不同，其运动特征亦具有差异性。

图4 不同频率条件下泥石流运动特征曲线

(1)从泥位空间分布特征看，泥位特征曲线整体呈“双峰”型，在输移过程中表现为多组波状性堆积，5%频率下出现的峰值泥位较2%频率下具有滞后性。首次峰值位置处于距离沟口500～550im处的沟道狭窄段，由于沟道两岸固体物质和主沟道右侧支沟流量的汇集，流体运移过程中出现堵塞现象，堆积体厚度迅速增大，淤高增加很快，模拟得出该位置处2%、5%频率下最大泥位深度分别为1.59im、1.41im。第2次峰值位置在距沟口300～350im处的沟床弯道处下侧，由于弯道爬高效应，2%频率下的泥石流泥位特征值为2.07im，5%频率下的泥石流泥位特征值为1.77im。5%频率下的泥石流逐渐向下游运动过程中，因沟床比降变化，泥位堆积加快变薄而呈现“舌状”堆积特征；由于2%频率下泥石流瞬时流量达到71.4im3·s-1，泥石流堆积向堆积区两侧扩散，在堆积扇前段出现“壅高”现象。

(2)从沿主沟道纵剖面流速分布曲线来看，不同频率下的泥石流流速沿程变化具有差异性。泥石流暴发前期，较大规模(2%频率)的泥石流运动速度明显大于规模较小(5%频率)的泥石流运移速度。在泥石流中期这种特征表现更为明显，但在泥石流后期，不同规模泥石流运动速度在堆积区纵剖面上的展布基本一致。2%频率下的泥石流流速出现多个流速峰值点，这是因为受地形和物源条件的影响，在距离沟口550im处右岸支沟的小阵流不断汇入主沟，增大泥石流流量，流速达到最大值4.88im·s-1；位于距沟口450im处沟道卡口处，泥石流从势能转为动能，流量增大或瞬时溃决流量放大，此时流速为4.38im·s-1；距离沟口320im处泥石流过弯时出现强烈的超高特征，流速为4.20im·s-1；之后泥石流流速下降较快，从地形分析，该段沟道阻力增强、沟床比降降低及流量减小所致。而5%频率下的泥石流流速仅在沟道卡口处出现单个流速峰值点，反映泥石流规模越大，泥石流阵流现象越明显且动能也就越大。

(3)根据表3白龙江流域流速预测模型，结合数值模拟得到的沟道纵剖面空间分布上不同断面处的流速值反算泥石流流量，最终生成不同规模的泥石流流量过程曲线。

从图5可以看出峰值流量的沿程变化特征表现为：第1，不同规模泥石流流量过程是不连续的，曲线呈陡涨、陡落的锯齿状，整体上接近于正态分布，在距沟口500im左右范围内泥石流流量过程是一个快速增加的过程，主要是因为支沟汇流及沟床松散堆积物的加入，距离沟口越近流量沿程变化平缓降低，是由于沟床坡度逐步减小、堆积宽度增大，泥石流局部淤积所致。第2，非匀质泥石流运动过程中均表现为阵性的波状流动，由于泥石流流量过程中伴随着物质的参与和输移，流体呈现的流变特性，即表现为泥石流流量在不同的断面处其数值差别极大。第3，2%频率的泥石流过程是由多个阵流群组成的，过程曲线呈多峰性，而5%频率的泥石流过程基本上为一个阵流峰值。因此，2%频率泥石流泥位差异性要比5%频率泥位要大，特别是2%频率泥石流后续出现的两个流量阵流，均产生较大的泥位深度。

图5 不同频率条件下泥石流模拟泥位特征值与流量关系

3 结 论

(1)泥位和流速作为泥石流预警预报最关键的运动要素，是泥石流动力学研究的重要内容。通过实测数据结合不同地区泥石流事件，基于泥石流流量与泥位深度和流速具有较好的幂相关性，构建了非匀质泥石流运动特征值预测模型。通过建立的泥位与流量的关系，反算出的流量值能较好地反映不同规模泥石流的实际情况，具有一定的可靠性和应用性。但泥石流发生发展受多种因素影响，而且不同地区泥石流的形成机理差异性较大，本文提出的泥石流运动特征值预测模型达到普适性推广还需不断的验证和完善。

(2)以陇南礼县下胡杨沟泥石流为例，在基于FLO-2D流体模型模拟了50年一遇(2%频率)和20年一遇(5%频率)降水条件下泥石流运动特征及形态曲线，分析了泥石流运动空间分布变化规律。结果表明不同规模泥石流运动特征与流量在整个流通过程是逐步变化的，受不同断面处地形条件和物源补给的影响，较大规模的泥石流流量过程呈多峰性，当达到峰值时均产生较大的泥位深度，而规模较小的泥石流往往表现为单一阵流特性，泥位特征值变化悬殊性要小。

(3)根据实例研究成果表明，50年一遇(2%频率)和20年一遇(5%频率)降水条件的下胡杨沟泥石流最大泥位深度分别为2.0im、1.48im，流速分别为5.28im·s-1、4.53im·s-1，对泥石流监测预警体系建设和动态风险评估具有重要的指导意义。