泥石流启动临界土体含水量及其预警机制分析

　0引言
　　
　　自然界中，泥石流往往由上游沟岸或坡面的固体物质失稳进入沟道，并在沟道水流的动力作用下形成。初始形成的泥石流规模不大，在运动过程中通过侵蚀、裹挟松散物质，规模才逐渐增长［12］。
　　
　　一般而言，坡面土体要达到一定的含水量(体积比，下同)才能启动形成泥石流。坡面土体的含水量是反映泥石流形成的直接参数。目前常用的前期降雨和当期降雨量预警指标，间接体现了降雨使土体含水量增加，抗剪强度和土体稳定性降低，从而使泥石流发生的可能性增大［36］。国内外许多泥石流预警模型多以“前期有效雨量雨强/雨强历时”为基本模式［710］。但是前期有效雨量与流域下垫面特性、岩土体物理特性等有关，其计算过程中涉及的递减系数、前期降水衰减系数等重要常数的确定需要长时间的土体含水量和降雨量观测，且不同学者使用的计算方法以及得到的结果可能不一样［1114］。不同泥石流沟的有效前期雨量衰减规律也不同。因此，前期降雨量预警模式存在局地性强、预警精度不高、参数不易确定等缺点。
　　
　　针对“有效前期雨量雨强/雨强历时”预警模型所存在的不足，笔者通过分析国内外大量泥石流形成过程中的土体含水量变化数据，直接以土体含水量为预警指标，提出泥石流启动临界土体含水量的概念和经验计算公式；经过云南东川蒋家沟泥石流观测数据的验证，所得到的方法和公式可应用于泥石流预警工作中。
　　
　　1泥石流临界土体含水量
　　
　　1.1泥石流形成过程中土体含水量变化
　　
　　Cannon等对美国加利福利亚、科罗拉多火灾地区大量野外土体含水量和泥石流暴发时间监测表明，不同深度坡面土体(一般监测最大深度为60 cm)在泥石流暴发时的含水量都未达到饱和［1517］。以加利福尼亚南部桑加布里埃尔山区2009年发生的泥石流为例，167 mm的累积雨量共造成7条流域产生泥石流，而流域中的土体含水量却始终维持在22%，远小于饱和度40%［18］。这些地区的泥石流在较好的地表植被条件下（比如火灾过后植被逐渐恢复后），可由径流触发转变为滑坡激发，1 h激发雨量也逐渐增大［19］；蒋家沟暴雨泥石流的多年监测也表明：降雨过程中，雨水在坡面松散土体最大可以下渗60~80 cm，斜坡坡脚位置的土体含水量最大，但是都未超过孔隙度 ［2021］；一些模型试验结果也表明，在土体启动时大部分土体处于未饱和状态［2224］。这与传统认识的泥石流形成理论并不一致，尤其从非饱和土力学角度来看，泥石流启动的原因应是土体饱和后短历时雨强造成孔隙水压力剧增、有效应力下降或丧失。事实上，降雨过程中土体含水量空间分布不均匀，垂直深度上土体含水量是非线性关系。对孔隙度大、强度低、颗粒间黏结程度差的坡面松散物源体来说，土体失稳不仅仅是由于含水量增加、抗剪强度减小、下滑力增加导致的，土体失稳在泥石流形成过程中只是初步阶段也是必不可少的阶段，降雨下渗在土体中形成壤中流以及表面形成地表径流也是非常重要的原因。在一些地区甚至径流作用在泥石流形成作用中更是起主导作用。例如，根据蒋家沟流域坡面径流和沟道汇流汇集过程的观测，大多是清水变成泥石流体的过程，在很多情况下，在坡面上就初步形成了小泥石流体或浆体［25］。美国加利福尼亚南部、科罗拉多地区火灾后的泥石流多是流域坡面土体含水量达到一定值后，由峰值降雨时段在引发一定水力条件的坡面径流激发形成的［26］。
　1.2临界土体含水量定义和经验计算公式
　　
　　无论是土力类还是水力类的泥石流，暴发时刻大多在峰值降雨时段附近。随着雨水入渗，土体含水量逐渐增加，当含水量达到某个临界值时，源区坡面土体达到极限平衡状态，或者坡面土体入渗和失水达到动态平衡。后续的降雨强度如果能够达到或者超过渗透率，坡面产生地表径流而且土体失稳下滑，形成最初的泥石流。这一临界含水量称为泥石流启动的临界土体含水量。国内外研究中有关临界土体含水量的概念已有所涉及。例如，在Brocca等提出的一些降雨径流模型中，土体水分平衡主要考虑了降雨入渗率、雨强、蒸发率以及因壤中流和深层土体渗透的土体排水率，临界土体含水量与饱和度并不一致［27］。杨大文等将遂川江流域实测土体含水量与不同时间的雨量结合起来，建立流域内不同点土体饱和度和警戒雨量的关系［28］。该结果间接反映了含水量越大，警戒雨量越小。实测的土体饱和度中也反映了临界土体含水量的概念。
　　
　　通过对国内外大量降雨激发浅表层滑坡、泥石流过程中土体含水量变化与土体物理参数的分析（表1），发现泥石流暴发时的土体含水量与土体渗透系数、孔隙度以及颗粒级配存在正相关关系，用Matlab的Stepwise函数作交互式逐步回归分析，得到多重线性关系式
　　Wa =-1.12K+0.46n+0.12Cc-0.165(1)
　　式中：Wa为临界土体含水量；K为土体渗透系数；n为土体孔隙度；Cc为土体颗粒的曲率系数，Cc=d230/(d10d60)，其中d10、d30和d60分别是土体颗粒质量累计含量（质量分数，下同）为10%、30%和60%时的粒径。
　　
　　式（1）拟合所用数据如表1，判定系数为0.896 5，调整判定系数为0.883，均方根误差为0.055 6。对该方程显著性进行F检验，查表得到F0.05(3,23)值为3028，远小于统计量F值（66.4），因此，回归公式显著。
　　
　　从式(1)可以看出，临界土体含水量与渗透系数成负相关关系，与孔隙度和土体颗粒曲率系数成正相关关系。渗透系数越大，说明土体在一定水力梯度下可以很快下渗并在短时间内转化为壤中流。因内部壤中流流动和深层渗漏的土体排水速度与雨强大小接近时，土体含水量才达到临界值，并短时间维持在相对稳定的水平。孔隙度越大，说明水、气两相占据土体内部空间越大。一般而言，在非饱和阶段，土体内部颗粒之间具有一定的黏聚力和咬合力并使土体具备一定的抗剪强度，从而保持稳定。国内外许多研究证明土体启动时是非饱和的，土体内部需要更多的水才能使抗剪强度下降，下滑力增加。因此，临界土体含水量与孔隙度成正相关关系说明水分在松散碎屑坡面体稳定性中具有重要作用。土体颗粒的曲率系数反映了土颗粒粒径分布曲线形态。从表1中曲率系数可以看出，泥石流源区粒径级配累计曲线斜率比较连续，细颗粒在泥石流源区土体仍 占很大部分。尽管级配连续，但是泥石流源区土体是不均匀的，存在不连续粒径，且均匀系数变化较大。由此说明临界土体含水量与曲率系数成正相关关系，曲率系数越大，土体内部存在不连续粒径，土颗粒间的存水空间越多，临界土体含水量越大。
　　
　　另外，式(1)中的3个变量综合反映了流域内土体的平均最大蓄水量。临界土体含水量实际上与平均最大蓄水量物理涵义一致。超过临界土体含水量或者平均最大蓄水量的降水将从地表流走，形成地表径流，并聚集形成足够水动力条件的沟道水流。坡面汇流而来的初步小规模泥石流浆体与沟道水流
　　
　　表1泥石流形成时的土体含水量、渗透系数、孔隙度以及曲率系数
　　Tab.1Soil Moisture, Permeability Coefficient, Porosity and
　　Coefficient of Curvature when Debris Flow Formed
　　编号WaK/(mm·s-1)nCc
　　1
　　
　　2
　　30.500 00.023 5100.753.333 333
　　0.590 00.014 4930.763.333 333
　　0.660 00.004 6670.693.333 333
　　0.540 00.007 6190.723.333 333
　　0.480 00.003 7040.763.333 333
　　0.540 00.005 0720.763.333 333
　　0.240 00.005 6260.542.000 000
　　0.133 00.005 6260.450.408 333
　　0.120 00.010 2600.580.888 889
　　0.113 00.087 7600.580.888 889
　　0.124 00.010 2600.580.888 889
　　0.103 60.087 7600.580.888 889
　　0.161 80.010 2600.580.888 889
　　0.152 60.005 6260.450.408 333
　　0.117 40.010 2600.580.888 889
　　0.152 90.005 6260.451.125 000
　　0.326 00.003 8190.691.481 481
　　0.310 00.003 8190.691.481 481
　　0.357 00.003 8190.691.481 481
　　0.237 00.003 8190.691.481 481
　　0.340 00.003 8190.691.481 481
　　0.340 00.003 8190.691.481 481
　　0.352 00.003 8190.691.481 481
　　0.378 00.003 8190.691.481 481
　　0.313 00.003 8190.691.481 481
　　0.347 00.003 8190.691.481 481
　　0.347 00.003 8190.691.481 481
　　注：编号1的数据引自文献［29］；编号2的数据引自文献［30］；编号3的数据引自文献［24］。
　　汇集增大了其中的含沙量和携带固体物质能力。沟道水流在流经泥石流动床时强烈侵蚀固体物质并形成泥石流。由于式（1）中的3个参数是泥石流源区土体物理特征值，所以该公式主要针对泥石流形成区的土体，尤其是细颗粒含量多的泥石流物源体。
　　2基于临界土体含水量和实时降雨的预警方法
　　根据临界土体含水量的定义以及泥石流预警模式中有效前期降雨难以精确计算的问题，笔者提出一种基于临界土体含水量和实时降雨的泥石流预警方法(图1)，该方法的具体流程如下。
　　
　　图1基于土体含水量和实时降雨的泥石流预警方法实施过程
　　Fig.1Flow Chart of Debris Flow Forecasting System Based on Critical Soil Moisture and Realtime Rainfall
　　
　　(1)计算泥石流启动的临界土体含水量。根据泥石流形成区的岩土体特征，利用式（1）计算对应临界土体含水量Wa。
　　
　　(2)计算当前土体含水量与临界土体含水量的差值。在得到临界土体含水量的情况下，通过土壤含水量传感器，以太阳能板和蓄电池作为电源，运用GPRS无线网络信息传输和室内数据接收终端，得到降雨过程中t时刻的土体含水量Wt。通过室内数据处理系统反复计算差值含水量ΔW，其表达式为
　　ΔW=Wa-Wt(2)
　　(3)根据当前雨强计算达到临界土体含水量所需要的时间T。具体计算方法为：通过雨量传感器，以无线网络方式发出和接收雨量值，通过室内系统判定t时刻雨强Rt和土体渗透系数K的大小；并选择大于渗透系数的计算公式和小于渗透系数的计算公式对达到临界土体含水量所需要的时间进行计算。
　　T＝(Wa-Wt)/KRt>K
　　(Wa-Wt)/RtRt≤K(4)
　　(4)根据计算的T值进行预警，并每隔一定时间重复第(3)、(4)步直至达到临界土体含水量或降雨结束。
　　如果T>0，则多通道数据反复确认降雨过程中土体含水量是否达到临界土体含水量，并每隔一定时间(比如10 s)重复第(2)、(3)步操作内容。如果降雨一直持续且T≤0，那么安排现场查看或发出泥石流即将发生的警报。如果降雨结束，则停止泥石流预警。
　　
　　
　　3实例分析
　　
　　笔者利用提出的基于临界土体含水量和实时降雨的泥石流预警方法，对云南东川蒋家沟1999年7月16日泥石流进行演算。
　　
　　3.1临界土体含水量确定
　　
　　蒋家沟角砾土密度为1.954 g·cm-3，孔隙度n为0.381 8，渗透系数k为0.008 07 mm·s-1，d10为0.01 mm，d30为0.25 mm，d60为3 mm，角砾土颗粒级配曲线的曲率系数Cc为3.125(图2)。由此得到蒋家沟泥石流暴发的临界土体含水量为401%。陈晓清等在蒋家沟人工降雨激发滑坡失稳形成泥石流的试验中发现，尽管不同深度的土体含水量不同，但每一组试验中土体含水量最大值都介于45%和35%之间，且土体破坏前的土体含水量大致在40%上下剧烈波动［23］。土体在降雨作用下达到破坏前状态，其含水量剧烈变动，但是变动的幅度基本在40%左右。这与通过式(1)计算得到的土体含水量值基本相同。因此，蒋家沟源区土体含水量接近越该值，泥石流暴发的可能性就越大。
　　
　　图2云南东川蒋家沟角砾土的颗粒级配曲线
　　Fig.2Particle Grading Curve of Breccia Soil in Jiangjiagou of Dongchuan, Yunnan
　　
　　3.2演算过程
　　
　　根据东川泥石流观测站提供的降雨数据，激发蒋家沟1999年7月16日泥石流的降雨过程见图3。
　　
　　图31999年7月16日泥石流暴发前的降雨过程和前20 d的雨量过程
　　Fig.3Rainfall Processes Before Occurrence of
　　Debris Flow on 16 July, 1999 and During the 20 Days Before the Debris Flow
　　
　　
　　(1)计算差值含水量。初始时土体含水量为4%，未达到临界土体含水量(40.1%)，因此，在此时没有泥石流发生，差值含水量ΔW为36.1%。
　　
　　(2)比较雨强、渗透系数大小并计算每单位时间的土体含水量。蒋家沟角砾土的渗透系数为0008 07 mm& middot;s-1，相当于每10 min降雨484 mm，与该沟的始发雨强（5 mm）非常接近。渗透系数0008 07 mm·s-1是经过原位渗透试验得到的参数，即土体达到稳定渗透阶段时的渗透系数。渗透系数是随时间和土体含水量变化的，但无论渗透系数在降雨过程中如何变化，稳定渗透阶段的渗透系数在不同测试手段下差异不是很大。比如陈宁生等经人工降雨试验测得降雨开始时的土体初始渗透系数为0009 2 mm·s-1 ［31］，相当于每10 min降雨552 mm的等效雨强。
　　
　　
　　通过实际10 min雨量(R10)与渗透系数的对比，以单位面积和单位垂直深度的土体作为分析对象得到每单位10 min末的含水量
　　Wt=WtΔt+ΔtRAHA(5)
　　式中：WtΔt为tΔt时刻(对应t时刻之初)的土体含水量；t为实际降雨记录时刻；Δt为单位雨量记录持续时间(这里为10 min)；R为雨强；
　　A为土体面积；ΔtRA为时间Δt内进入土体的降水体积；HA为土体的总体积，H为一般含水量传感器的探针长度（60 mm)。
　　
　　(3)根据笔者提出的预警方法演算。由于蒋家沟没有实测的泥石流暴发过程土体含水量变化，这里在得到临界土体含水量、确定基本土体含水量和实际雨量过程后，根据提出的方法具体流程进行演算。演算结果见表2。
　　
　　表2基于临界土体含水量和实时降雨的泥石流预警方法演算过程
　　Tab.2Calculation Process by the Forcasting
　　Method for Debris Flow Based on Critical Soil Moisture and Realtime Rainfall
　　编号时间段R10/mmWtΔtWa-WtΔtT/min
　　121:50~22:000.50.040 00.361 0433.200 0
　　222:00~22:101.70.048 30.352 7124.482 4
　　322:10~22:201.60.076 70.324 3121.612 5
　　422:20~22:301.10.103 30.297 7162.381 8
　　522:30~22:400.80.121 70.279 3209.475 0
　　622:40~22:500.90.135 00.266 0177.333 3
　　722:50~23:001.30.150 00.251 0115.846 2
　　823:00~23:102.00.171 70.229 368.790 0
　　923:10~23:204.40.205 00.196 026.727 3
　　1023:20~23:301.40.278 30.122 752.585 7
　　1123:30~23:401.10.301 70.099 354.163 6
　　1223:40~23:503.00.320 00.081 016.200 0
　　1323:50~次日00:002.30.370 00.031 08.087 0
　　14次日00:00~次日00:101.30.408 3<0.000 0<0.000 0
　　注：次日01:12:34时刻，监测到泥石流；R10值始终小于渗透系数。
　　3.3精度比较
　　
　　将本文提出的预警方法与蒋家沟泥石流预报临界线和暴发线判别式进行结果对比，得到以下结果
　　R10=5.5-0.098(Pa0+h)>0.5 mm(6)
　　R10=6.9-0.123(Pa0+h)>1.0 mm(7)
　　式中：Pa0为泥石流暴发前某一天的指数；h为泥石流暴发前的当日降雨量。
　　
　　临界线式(6)的物理意义是：在10 min雨强大于0.5 mm的降水过程中，某10 min降水量只要等于5.5-0.098(Pa0+h)，则蒋家沟泥石流就可能暴发；暴发线式(7)的物理意义为：在10 min雨强大于1 mm的降水过程中，某10 min降水量只要等于69-0.123(Pa0+h)，则蒋家沟就会暴发泥石流。
　　
　　
　　利用该次泥石流过程之前的降雨量过程和蒋家沟前期雨量计算公式，得到临界线和暴发线10 min雨量分别为2.26、2.86 mm(图3)。从该次降雨过程来看，23:10~23:20时段4.4 mm的降雨是造成该次泥石流的主要原因。但泥石流暴发时并不与该次降雨过程的峰值雨量时段重合，滞后近1 h。表1演算结果表明，该时段土体含水量并未达到临界土体含水量。而泥石流暴发时段雨强仅0.9 mm，大于临界线而小于暴发线。从达到临界土体含水量和泥石流暴发时间上来看，达到临界土体含水量和泥石流暴发之间相差约1 h，而原方法可提前预警17~200 min［8］。因此，从临界土体含水量结合实时降雨过程来判别泥石流的发生更准确。笔者提出的基于临界土体含水量和实时降雨的预警方法比传统利用临界线和暴发线判别泥石流的物理意义更明确，方法更可靠。
　　4结语
　　
　　（1）在泥石流形成过程中，一般理论认为前期降雨使源区土体饱和，短历时雨强造成饱和后的土体产生高孔隙水压力使土体失稳并转化为泥石流。但是国内外大量泥石流形成过程监测表明，源区坡面土体降雨过程和泥石流形成过程中土体都未达到饱和，且土体含水量存在一个临界值。由此，本文提出了临界土体含水量的概念，并通过拟合国内外野外监测数据得到计算临界土体含水量的经验公式。
　　
　　（2）基于临界土体含水量的概念和经验公式，发展了一种基于临界土体含水量和实时降雨的泥石流预警方法。通过云南东川蒋家沟1999年7月16日暴发的泥石流的实际观测资料，对该方法进行了实例演算。演算结果表明，该场泥石流暴发时刻并未与峰值降雨时段重合，而是在达到临界土体含水量后约1 h。
　　
　　（3）由于临界土体含水量计算公式是经验性的，在后续研究中有必要从土体降雨入渗以及激发坡面土体失稳的物理过程并结合泥石流形成区监测进行深入研究，提出更具有物理意义的临界土体含水量概念，建立以水文学、水力学、泥沙运动学等为基础的预警模型。
　　参考文献：
　　References:
　　［1］崔鹏.泥石流起动条件及机理的实验研究［J］.科学通报,1991,36(21):16501652.
　　CUI ment Study on the Mechanism and Condition of Starting Up of Debris Flow［J］.Chinese Science Bulletin,1991,36(21):16501652.
　　［2］崔鹏.中国山地灾害研究进展与未来应关注的科学问题［J］.地理科学进展,2014,33(2):145152.
　　CUI ss and Prospects in Research on Mountain Hazards in China［J］.Progress in Geography,2014,33(2):145152.
　　［3］CAINE Rainfall Intensity:Duration Control of Shallow Landslides and Debris Flows［J］.Geografiska A,Physical Geography,1980,62(1/2):2327.
　　［4］WIECZOREK G F,GLADE icc Factors Influencing Occurrence of Debris Flows［M］∥JAKOB M,HUNGR flow Hazards and Related Pheno:Springer,2005:325362.
　　［5］马超,胡凯衡,宋国虎,等.汶川地震灾区帽壳子滑坡形成泥石流的过程和特征［J］.地球科学与环境学报,2013,35(4):98103.
　　MA Chao,HU Kaiheng,SONG Guohu,et ses and Characteristics of Debris Flows Induced by Maoqiaozi Landslide in Wenchuan Earthquake Stricken Area［J］.Journal of Earth Sciences and Environment,2013,35(4):98103.
　　［6］谢洪,刘维明,赵晋恒,等.四川石棉2012年“7·14”唐家沟泥石流特征［J］.地球科学与环境学报,2013,35(4):9097.
　　XIE Hong,LIU Weiming,ZHAO Jinheng,et cteristics of Tangjiagou Debris Flow in Shimian of Sichuan in July 14,2012［J］.Journal of Earth Sciences and Environment,2013,35(4):9097.
　　［7］谭万沛,王成华,姚令侃,等.暴雨泥石流滑坡的区域预测与预报:以攀西地区为例［M］.成都:四川科学技术出版社,1994.
　　TAN Wanpei,WANG Chenghua,YAO Lingkan,et al Forecasting and Predicting of Rainfall Induced Debris Flows and Landslides:Take the Western Panzhihua as an Example［M］.Chengdu:Sichuan Science and Technology Press,1994.
　　［8］陈景武.降雨预报泥石流的原理及方法［C］∥中国科学院水利部成都山地灾害与环境研究所.第二届全国泥石流学术会议论文集.北京:科学出版社,1989:8490.
　　CHEN Jing Principle and Method of Debris Flow Forecast Based on Rainfall［C］∥Institute of Mountain Hazards and Environment,Chinese Academy of Sciences and Ministry of Water ding of the Second National Debris Flow g:Science Press,1989:8490.