甘肃黑方台黄土水分运动参数试验研究

李志萍，王欢，曾磊，赵明星，韩翔宇，周翔

(1．华北水利水电大学，河南 郑州 450045;2．中国地质调查局西安地质调查中心，陕西 西安 710054)

灌溉渗透诱发型黄土崩滑灾害在黄土高原具有普遍性，大量分布于陕、甘、宁、青等省的主要灌区周边，因长期引水灌溉导致地下水位上升，引发台缘周边斜坡地带滑坡频发，规模以大、中型黄土崩滑为主，且多呈带状群发，甚或形成连绵数十公里长的滑坡带．其中以黑方台灌区最具典型性和代表性．

灌溉渗透诱发型黄土崩滑灾害机理具有特殊性，且具有群发性，造成的损失严重，对人民生命财产构成极大威胁，已经成为西北黄土高原乃至全国地质灾害最为频繁和严重的地段之一，在一定程度上制约了区域社会经济的发展．

从非饱和土力学角度看，灌溉导致黄土含水量增加，亦即引起基质吸力的变化，进而引起土体强度发生改变．有研究指出，在相同围压下，饱和黄土在更低的偏应力水平下发生破坏，而偏应力水平相等时，饱和黄土的蠕变量远大于天然含水量黄土，表明黄土在灌溉后更易发生蠕变变形和破坏．这说明，黄土含水量这个土壤水分运动参数会对黄土自身的变形和破坏产生一定的影响，进而造成黄土崩滑灾害发生的可能．

甘肃省临夏回族自治州永靖县盐锅峡镇黑方台为湟水河与黄河交汇处的四级阶地台面，北为湟水河，南为虎狼沟，东西分别被黄河和磨石沟沟谷切割，整体台面为“孤岛”状旱台，台面总体较为平坦、宽广，其南部前缘与黄河阶地相连，高差大于100 m．台面常年灌溉引发台缘周边形成长达十余公里的滑坡群，严重威胁着台缘周边群众的生命财产安全以及黄河八盘峡库区的安全运营．因滑坡形成时代新，诱发因素典型，发生频率高，每年都要发生3 ～5 起．区内交通条件相对较好，有县级公路通过．

此次试验选取黑方台原状黄土与扰动黄土作为试验介质，利用德国公司生产的Ku －pF 非饱和导水率测定系统测定这两种试验介质的土壤水分特征曲线并进行对比分析，采用Van 方程参数拟合方法，结合MATLAB 软件中的非线性最小二乘法函数lsqcurvefit 分别对原状样与扰动样的实测土壤水分特征曲线进行拟合并对比分析，并且得出相关拟合参数，以期为该类型崩滑灾害的综合整治提供土壤水分运动方面的理论依据和技术支撑．

1 试验方法

试验主要利用Ku－pF 非饱和导水率测定系统测定土壤水分特征曲线，主要试验方法与步骤如下．

1．1 取 样

将1 号环刀垂直打压入土壤中，直至环刀完全没入土壤，再用铁锹挖出，回到实验室，用钢刀将1号环刀两侧多余的土抹掉，制备原状样;事先称重小环刀，用小环刀取原状样，取样步骤同1 号环刀;再取一些扰动土带回实验室，以备第2 步制样时使用．

1．2 制 样

1)将取有原状土的小环刀与扰动土一同放入烘箱中(105 ℃)烘干24 h．

2)称重取有原状土的小环刀，计算小环刀中原状土的干容重，然后换算出2 号环刀中需要装入烘干后扰动土的质量(与原状土样质量不能相差太多，应比原状土样重)，再加上事先已称重的底盘、滤纸、环刀的质量，计算出2 号试验样品的总质量．

3)将烘干的扰动土装填入2 号环刀内，再将2号试验样品置于电子天平上称重，以控制质量，如此反复，达到上一步已经计算出的2 号试验样品的总质量为止．

4)将1 号、2 号试验样品用蒸馏水浸泡若干天，直至土样的上部看到水或能够摸到水为止，使土样完全饱和．

1．3 装 样

将1 号、2 号环刀加底上盖，并与张力计联接，放置于Ku－pF 非饱和导水率测定系统的星型吊臂上，通过传感器，将所测得的试验数据(Tension top、Tension bottom、Weight)传输到电脑控制界面．

2 试验数据处理

2．1 土壤水分特征曲线

土壤水的吸力或负压随土壤含水量的变化而变化，其关系曲线称为土壤水分特征曲线或土壤持水曲线．土壤水分特征曲线表示土壤水的能量和数量之间的关系，是反映土壤水分保持和运动的特征曲线，其基本形状如图1 所示． 由图可以看出，负压h与含水量θ 的关系曲线分脱湿曲线和吸湿曲线，脱湿曲线和吸湿曲线是土壤脱湿(由湿变干)过程和吸湿(由干变湿)过程测得的水分特征曲线．

脱湿曲线:θs表示饱和含水量，θr表示残余含水量，当负压逐渐增大时，含水量逐渐减小．h －θ 曲线自θs开始上升，以θ =θr为渐近线．吸湿曲线:当负压达到一定值之后，再慢慢减小负压，含水量开始增加，直至h 近于零．此时的h －θ 曲线并不与脱湿过程的h－θ 曲线重合．

图1 土壤水分特征曲线示意图

此次做脱湿过程的试验．当土壤中的水分处于饱和状态时，含水量为饱和含水量θs，而负压为零．若对土壤施加微小的负压，土壤中尚无水排出，则含水量维持饱和值． 当负压增加至某一临界值hc后，由于土壤中最大孔隙不能抗拒所施加的负压而继续保持水分，于是土壤开始排水，相应地含水量开始减小．当负压进一步提高，次大的孔隙接着排水，土壤含水量随之进一步减小． 如此，随着负压不断增加，土壤中的孔隙由大到小依次不断排水，含水量越来越小．当负压很高时，只在十分狭小的孔隙中才能保持着极为有限的水分［1］．

利用Ku－pF 非饱和导水率测定系统所测得的试验数据，可经简单的计算处理得到一系列hi和θi值．计算公式［1］如下:

式中:hi为负压，cm;θi为含水量，cm3/cm3;i 为时间，min;(h1)i和(h2)i分别为i 时刻上、下端张力计的张力，即(Tension top)i和(Tension bottom)i，cm;W水i为i 时刻土壤水的质量，g;Wi为i 时刻试验样品和仪器的总质量，即(Weight)i，g;W0为烘干后的土样和仪器(包括提篮、上盖、底盖、滤纸、环刀)的总质量，g;V水i为i 时刻土壤水的体积，cm3;V环刀为环刀的体积，cm3．

根据达西定律［1］，可算出i 时刻土壤非饱和渗透系数．

式中:Ki为i 时刻土壤非饱和渗透系数，cm/min;vi为i 时刻渗透流速，cm/min;L 为渗流路径的直线长度，在此次试验中为环刀侧壁纵向排列的两个圆孔之间的距离，cm;ΔHi为i 时刻渗流路径始末断面的总水头差，在此次试验中为i 时刻上、下两端张力计张力总和的一半，cm;ΔHi/L 是i 时刻的水力梯度;Hi为i 时刻的总水头或总水势，cm;T 为时间，min;S环刀为环刀的底面积，cm2．

根据计算得到的hi、θi与Ki，可绘制土壤水分特征曲线，即负压h 与含水量θ 的关系曲线，h －θ曲线;渗透系数K 与含水量θ 的关系曲线，K －θ 曲线;渗透系数K 与负压h 的关系曲线，K－h 曲线．由于此次试验主要研究原状土与扰动土的土壤水分特征曲线，因此，不再绘制K－θ 曲线与K－h 曲线．

2．2 数据拟合

2．2．1 拟合公式

常用的描述土壤特征曲线的经验公式主要有Garder、Brooks－Corey、Campbell、Van Genuchten、Garder－Russo、Hutson－Cass 等6 种类型［2－5］．

在已经建立的众多数学模型中，Van Genuchten模型以其线型与实测数据曲线拟合程度好而得到广泛应用［6－7］．

Van Genuchten 公式［2］:

式中:θ 为土壤体积含水量，cm3/cm3;h 为压力水头，cm，饱和带中为正值，非饱和带中为负值，潜水面上为0;θs、θr分别为饱和含水量和残余含水量，cm3/cm3;α 为进气值的倒数，即α 越小，进气值越大，土壤持水性能越强，cm;n 为经验拟合参数，无量纲;m=1 －1/n，0 ＜m ＜1．

其中，h 与θ 是由所测得的试验数据经简单的运算得到;θs、θr、α、m、n 等基本参数是通过MATLAB 软件中的非线性最小二乘法函数lsqcurvefit 拟合得到．

对于非饱和渗透系数K，拟合公式如下:

式中:s 为饱和度;Ks为饱和渗透系数，cm/min．

式(6)中，s 是通过MATLAB 软件中的非线性最小二乘法函数lsqcurvefit 拟合得到，Ks由MATLAB软件拟合程序的反复调参过程中得到．

3 试验结果分析

3．1 实测负压－含水量分析

图2 为负压－含水量关系曲线．由图可以看出，不同试验介质在试验过程中，负压随着含水量的减小而增大．由于实验仪器测量范围有限，难以测得试验样品的残余含水量θr和饱和含水量θs及其附近的值，所以根据试验数据绘出的水分特征曲线的形状与理论水分特征曲线存在差异，仅为理论水分特征曲线的中间部分．

扰动样的负压－含水量关系曲线在原状样的上方．也就是说，在实验仪器测量范围内，所测得的扰动样的最小含水量与最大含水量的数值均比原状样的大，而且，在相等的含水量条件下，扰动样的负压值也较原状样的要高． 造成这种现象的原因主要是天然状态黄土中存在着大孔隙，而原状土被扰动重塑后，消除了黄土中分布不均匀的大孔隙，土颗粒重新组合形成分散结构，并形成较多细小、连通性不好的小孔隙，失去了天然状态黄土所具有的结构性［8］．当对试验样品施加一定的吸力(或负压)时，原状样中的大孔隙不能抗拒所施加的吸力而继续保持水分，率先开始排水，相应地含水量开始减小． 当吸力进一步提高，扰动样中的小孔隙开始排水，含水量开始减小．

图2 负压－含水量关系曲线

3．2 拟合负压－含水量分析

图3 为拟合负压－含水量关系曲线．由于实验仪器测量范围有限，难以测得试验样品的残余含水量θr和饱和含水量θs及其附近的值．但是，土壤的θr和θs这两个水分参数不仅对于黑方台地区黄土崩滑灾害防治，而且对于科研工作均具有非常重要的意义．尤其是残余含水量θr，它与农田水利专业领域里的植物凋萎系数具有非常密切的关系． 掌握了θr和θs这两个土壤水分参数，对于黑方台地区黄土崩滑灾害防治、农田灌溉、优化配置水资源均具有重要的指导意义［9－11］．

图3 拟合负压－含水量关系曲线

图3 是在实测数据的基础上，运用MATLAB 软件中的非线性最小二乘法函数lsqcurvefit［12］，套用Van Genuchten 公式拟合出来的负压－含水量关系曲线．由图3 可以看出，当拟合负压控制在5 000 cm时，原状样和扰动样的负压－含水量曲线均出现了各自的纵向渐近线，其与横坐标含水量的交点即分别为原状样和扰动样的残余含水量θr．可以看出，扰动样的残余含水量较原状样的小，说明扰动样的失水率较大．当对试验样品施加的负压达到1 000 cm时，原状样几乎不再失水，基本达到残余含水量;随着负压进一步提高，扰动样继续失水，当负压为2 000 cm时，基本达到残余含水量．

3．3 拟合参数分析

表1 为拟合参数，原状样与扰动样的拟合参数和图3 中二者的拟合曲线对比，结果基本吻合．原状样的残余含水量较扰动样的大，而饱和含水量较扰动样的小．原状样的饱和渗透系数较扰动样的小．原状样的拟合参数α 和n 均较扰动样的大．这些拟合参数的变化规律，对实际工程建设具有很重要的参考价值．

表1 拟合参数表

但是，需要指出的是，Van 方程参数拟合方法虽然具有较好的拟合效果，然而，由于所选初值有时对拟合值有较大的影响，因此在这种情况下就需要多次调初值才可获得较为满意的拟合效果，这是该方法的局限所在．不过，该程序所具有的友好人机交互界面使这种调整变得方便而快捷．另外，由于θ 和K的数量级相差较大，而且存在权重问题，利用Van方程参数拟合方法拟合K 值存在较大的误差，需要进一步寻求更好的处理方法，以完善和改进该种求解参数的方法［7］．

4 结 语

1)原状黄土中的大孔隙对土壤水分特征曲线具有显著影响．

2)扰动黄土的残余含水量较原状黄土的小，而饱和含水量较原状黄土的大． 说明原状黄土被扰动重塑后，失水率变大，持水性能变强．

3)原状黄土所具有的特殊的结构性随着水体的侵入而逐渐丧失，颗粒重组． 因此，在湿陷性黄土地区地基处理等工程中，应针对不同工程问题，采取相应措施．

4)此次试验采用张力计法测定土壤水分特征曲线．由于张力计测定的负压范围为0 ～500 cm，因此，它只能用来测定低负压范围的水分特征曲线，在高负压段，根据实测数据，运用MATLAB 软件进行拟合．

5)对于负压－含水量实测曲线，Van 方程参数拟合方法具有较好的拟合效果，但是拟合K 则存在较大的误差． 因此，需要进一步寻求更好的处理方法．

［1］雷志栋，杨诗秀，谢森传． 土壤水动力学［M］． 北京:清华大学出版社，1988．

［2］彭建平，邵爱军．用MatLab 确定土壤水分特征曲线参数［J］．土壤，2007，39(3):433 －438．

［3］来剑斌，王全九．土壤水分特征曲线模型比较分析［J］．水土保持学报，2003，17(1):71 －74．

［4］李云龙，郭春颖，徐敏． 非饱和带水分特征曲线经验公式研究［J］．中国矿业，2010，19(8):105 －109．

［5］魏义长，刘作新，康玲玲，等． 土壤持水曲线Van Genuchten 模型求参的MatLab 实现［J］． 土壤学报，2004，41(3):380 －386．

［6］肖建英，李永涛，王丽．利用Van Genuchten 模型拟合土壤水分特征曲线［J］．地下水，2007，29(5):46 －47．

［7］李春友，任理，李保国．利用优化方法求算Van Genuchten 方程参数［J］．水科学进展，2001，12(4):473 －478．

［8］王辉，岳祖润，叶朝良． 原状黄土及重塑黄土渗透特性的试验研究［J］．石家庄铁道学院学报，2009，22(2):20－22．

［9］戴长雷，邓柏旺，李治军．包气带垂向含水率日变化特征试验研究［J］．黑龙江水专学报，2009，36(2):16－20．

［10］宋孝玉，李亚娟，蒋俊，等．非饱和土壤水分运动参数空间变异性研究进展与展望［J］．地球科学进展，2008，23(6):613 －618．

［11］刘苏峡，毛留喜，莫兴国，等．黄河沿岸陕豫区土壤水分的空间变化特征及其驱动力因子分析［J］．气候与环境研究，2008，13(5):645 －657．

［12］王金生，杨志峰，陈家军，等．包气带土壤水分滞留特征研究［J］．水利学报，2000(2):1 －5．