砒砂岩坡面侵蚀产沙规律及多元回归估算模型研究

梁止水,陈煜,孙悦,高海鹰,吴智仁

(1.东南大学土木工程学院,南京 210096;2.江苏大学环境与安全工程学院,江苏镇江 212013)

黄河流域砒砂岩区处于黄土高原到西北沙漠类型的过渡区,环境异质性极为突出,生态退化和水土流失非常严重,侵蚀模数可达到30 000～40 000 t/(km2·a),是我国乃至世界上侵蚀最为剧烈的地区之一,也是黄河下游“地上悬河”的粗泥沙来源核心区[1]。长期以来,砒砂岩区被列为国家生态环境建设和水土流失治理的重点区,水土流失防治取得了较为显著的治理成效[2-4]。

针对砒砂岩的水力侵蚀,主要集中在冲刷条件下的坡面侵蚀规律和水力学特征。其中:苏涛等[5-6]利用室内和野外径流冲刷试验研究了砒砂岩坡面的侵蚀规律和径流水动力学特性。杨吉山等[7]通过野外放水冲刷试验对比研究了白色和红色原状砒砂岩坡面的产流、产沙过程,分析了砒砂岩剥蚀率与水力学参数之间的相似关系。Liang等[8]采用室内模拟降雨试验分析了不同降雨强度、坡度、植被覆盖率条件下的产沙规律。杨振奇等[9]运用系统聚类和线性回归相结合的方法,划分裸露砒砂岩区降雨类型,研究了不同降雨类型对裸露坡面产流和产沙的影响,分析了不同降雨条件下坡面的微地形变化规律。李龙等[10]研究了自然降雨条件下砒砂岩坡面细沟微形态及其侵蚀特征,分析了细沟微形态变化过程对产流产沙的影响。董晓宇等[11]采用野外原位模拟冲刷试验研究了裸露砒砂岩区坡面侵蚀过程中地表粗糙度与水力侵蚀特征参数的关系。然而针对砒砂岩坡面的水土流失模拟及侵蚀产沙预测模型等研究较少,其中只有叶俊道等[12]对WEPP模型在砒砂岩地区不同坡面的土壤侵蚀进行了适用性探讨。因此,亟需针对砒砂岩坡面的特点建立适用于该地区的产沙过程预报模型。

为了更好地研究砒砂岩坡面在降雨侵蚀作用下侵蚀产沙规律,本文开展室内模拟降雨试验,研究不同砒砂岩边坡的坡度、降雨强度及植被覆盖率条件下的侵蚀产沙量,利用多元回归分析方法,分析不同参数与侵蚀产沙量之间的相互关系,初步建立侵蚀产沙估测模型,以期为砒砂岩水土流失区治理工作提供一定的科学依据。

1 材料及方法

1.1 试验材料

试验用砒砂岩土取自黄河中游上游皇甫川流域的圪秋沟,含水量约为9.0%,土壤颗粒级配使用马尔文3000激光粒度仪测定,＞1.00 mm 的颗粒组成占0.83%,0.50～0.05 mm 的颗粒组成占90.74%,＜0.01 mm 的颗粒组成占8.43%。所选用的植物为狗牙根〔Cynodondactylon(L.)Persoon〕,它与砒砂岩区广泛生长的野牛草极为类似,也具有根系较小、叶面窄长、移植成活率高等特点。

1.2 试验设计

采用标准小区5 m×1 m,进行室内模拟降雨试验,主要考察降雨强度、坡度和植被覆盖度的影响。根据全面调查立体条件和气候条件的基础上,砒砂岩区相对稳定坡度多在30°左右[10]。因此,试验中考察砒砂岩的坡度为20°,30°,40°。根据当地多年降雨数据资料选择降雨强度为20,50,80 mm/h,植被覆盖率为10%,30%,50%。此外,在砒砂岩区沟坡系统的陡坡面上很少有植被生长,因此,仅在20°的缓坡条件下进行植被覆盖度的影响研究。试验设计见表1,砒砂岩裸地坡面依次编号为A1—A9,有植被覆盖坡面依次编号为B1—B9。

表1 试验设计方案Table 1 Design of experiments

1.3 试验过程与测试方法

模拟降雨试验在黄河水利科学研究院模型黄河降雨试验大厅进行,其中模拟降雨装置采用下喷式的自动模拟降雨系统,其降雨的均匀性可达到90%以上。试验用的土槽为移动式可变坡度的钢槽,其坡度调节范围为5°～45°,土槽的长为5.0 m,宽为2.0 m,高为0.6 m,宽度方向上分为2个1.0 m 的土槽。狗牙草则采用试验前一个月的苗子进行移植,并通过对种植时的密度控制及叶面修剪来控制覆盖度大小,从而保证试验中的植被覆盖率要求。

模拟砒砂岩坡面采用控制容重的方法填土,层层压实,容重均控制在1.4～1.5 g/cm3,使其接近野外实际砒砂岩容重。为便于渗水,填土之前,在坡面底部预先铺上10 cm 厚的沙子。每场降雨试验的前一天,先用非常小的雨量进行预降雨,在基本保证不产生侵蚀的条件下直至形成产流为止,从而可以有效消除不同工况下的边界差异。同时,在每次降雨前,均需要用塑料薄膜将坡面盖住,并在坡面的4个角放置雨量筒,初始降雨15 min,用于校核降雨强度大小及其均匀性。降雨持续1 h,降雨过程中每3 min接一个产流泥沙样,采用集沙桶收集产流产沙情况,后用烘干法测得产沙量。

1.4 数据分析和建模方法

模拟试验结果利用SPSS软件和非线性回归分析方法进行拟合,建立回归方程。

2 结果与分析

2.1 砒砂岩裸地坡面产沙特征分析

从图1中可以看出,不同坡度条件下,降雨初期单位时间的产沙量均随降雨历时的增加不断增加,随后逐渐趋于稳定。其中,降雨强度越小的情况下,单位时间内的产沙量最终趋于稳定的时间越早。而在降雨历时半小时以后,所有工况下的单位侵蚀产沙量基本都趋于稳定。已有研究表明[13],黄土类坡面的侵蚀产沙量随降雨历时的变化曲线会有3种形式:平缓型、单峰型和多峰型,而砒砂岩坡面侵蚀产沙量随时间变化均为平缓型。在本次试验研究的时间范围内尚未出现峰值,这其中可能的原因是砒砂岩在遇水后短时间内出现侵蚀溃散,并呈现出侵蚀恶化的趋势,直至达到水流可携带泥沙的极大值。另外,侵蚀产沙量受降雨强度的影响为极其显著(p＜0.01),即为降雨强度越大,砒砂岩坡面形成的产沙量也越大,曲线波动的程度也更加明显。而在试验中的3种降雨强度条件下,侵蚀产沙量均随坡度的增加而呈现轻微的增加趋势,但该趋势并不明显(p＞0.05),说明坡度对侵蚀产沙量的影响小于降雨强度的影响[14]。

图1 裸坡坡面产沙量随降雨历时的变化曲线Fig.1 Curve of sand yield on the bare slope with the duration of rainfall

2.2 砒砂岩植被覆盖坡面产沙特征分析

图2为不同降雨强度和植被覆盖度下,砒砂岩坡面产沙量随时间的变化曲线,从图中可以看出,砒砂岩植被覆盖坡面产沙量随时间变化规律同裸露坡面的情况相似,其中初期产沙量随降雨历时的增加而不断增大,之后逐渐趋于稳定。砒砂岩坡面植被减蚀效果也明显受降雨强度影响,不同降雨强度时,植被的减蚀率相差较大。当降雨强度为20 mm/h时,不同覆盖度的植被减沙效果均非常明显(p＜0.01),所有植被覆盖度的减沙率均超过70.00%。但当降雨强度为80 mm/h时,只有植被覆盖度为50%对应的产沙量与裸坡相比变化明显,减沙率为86.72%,其他覆盖度条件,产沙量与裸坡相比变化并不明显(p＞0.05)。

图2 植被覆盖坡面产沙量随降雨历时的变化曲线Fig.2 Curve of sediment yield on vegetation cover slope with rainfall duration

3 讨论

根据模拟降雨试验结果,得到了累计产沙量随降雨时间的关系(图3),可以看出,不同坡面的累计产沙量随降雨历时增加而增加,且有增加变陡的趋势(曲率增加),这是由于随着时间推移,坡面侵蚀程度增大,径流中所含的泥沙在输移过程中增强了侵蚀动力,加剧了侵蚀产沙,该结果与肖培青等[15]得出累计产沙量与降雨历时成幂函数的结论基本一致。

图3 累计产沙量随降雨历时的关系曲线Fig.3 Relationship curve of cumulative sediment production with rainfall duration

3.1 累计产沙量拟合曲线

由图3看出,在不同的试验工况下,每分钟产沙量相差很大,导致这些累计产沙与降雨历时的关系曲线偏离也很大,难以统一分析。如果采用幂函数拟合,每一个试验条件都得到一个具有两个参数的公式,得不到统一的对应关系。因此,考虑对累计产沙量进行无量纲的归一化处理,使得数值在0～1[16],即为:将累计产沙量值都除以其1 h总产沙量,降雨历时以1 h的百分数表示。并得到了归一化后的关系曲线(图4),所有工况条件下得到曲线重合性好。

图4 累计产沙量随降雨历时变化的归一化关系曲线Fig.4 Normalized relationship curve of cumulative sediment production change with rainfall duration

由图4可以看出,归一化曲线基本都位于y=x下方,且非常接近以其为玄的圆弧,因此考虑采用曲线拟合,并假定为一条圆弧,圆心设为(a,b),在曲线的左上方,半径为r,其中圆弧同时经过点(0,0)和(1,1),则该弧线上的点(x,y)满足以下关系:

且有a≤0,b＞0,b＞y。

从而获得目标的拟合方程为:

对每条归一化后的曲线按照公式(2)进行非线性拟合,得到a值和拟合相关系数R2(表2),其中相关系数均在0.95以上,因此拟合的误差较小。而B1,B2和B9这3个组别中a值特别大,说明该3种工况下,圆弧曲率很大,曲线接近于直线,且都属于具有覆盖率的两个边缘情况,说明该拟合曲线存在一定的临界条件。在试验组的基础上,去除该三组数据后,对其他数据进行求平均值,从而保证其统一拟合的效果,从而得到a值约为-2.480,方差为3.154,各组数据基本均匀分布在拟合曲线的附近,且符合归一化曲线下凹、斜率逐渐增大的特点,可以较好地反映试验条件下累计产沙量随降雨历时的变化规律,故可以得到砒砂岩坡面侵蚀的累计产沙量与降雨历时(t)的拟合曲线方程为:

表2 拟合参数a 值及相关系数Table 2 Fitting parameter a values and correlation coefficients

式中:t为降雨历时(h),取值范围0～1。

3.2 多元非线性回归分析模型

多元非线性分析是利用数理统计的方法建立多个自变量与一个因变量之间非线性的函数关系,在诸多领域得到广泛应用[17-19]。

(1)传统方法是采用曲线化直线,利用6种初等函数对自变量进行变换,然后根据经典的最小二乘法原理,建立预测因变量与变换后自变量的线性关系。这种模型有一定精度,但没有考虑因素间的交互影响。其数学表达式一般为:

(2)另一种有效的方法是直接提出非直线形式的期望函数,作为模型框架,采用高斯-牛顿法进行参数估计,反复迭代直至参数估计收敛。该方法的方差分析F值和p值的意义不大。其数学表达式一般为:

影响坡面水力侵蚀产沙的因素,主要包括降雨、坡度、土壤特性(含水量、粒径大小、黏粒含量、内摩擦角等)、植被,以及坡长和人为作用等因素[20]。试验条件下,降雨强度、坡度和植被覆盖度是影响产沙量的3个最主要因素。大量的研究表明,侵蚀产沙量与降雨强度间存在幂函数关系[21-22],并在该研究中得到了验证。对于具有临界坡度的情况下,在临界坡度之前,与坡度间可能存在二次函数关系[23];而与植被覆盖度间的关系比较复杂,尚未明确,为方便之后建模,暂定为线性关系。预测模型的因变量为总侵蚀产沙量Y(单位kg),自变量分别为:降雨强度I(为方便计算,单位选择mm/min);坡面坡度数值S;植被覆盖度的百分数C。

采用第一种方法(单因素分析,表明Y与I2.5呈线性关系),可得到回归函数如下:

其相关系数R2=0.839,拟合效果一般,且各参数的标准误差均很大,用作回归模型不太理想。

采用第二种方法,需预先确定期望函数。根据前面对不同条件下产沙量变化过程的分析,降雨强度是主因,且自变量中降雨强度同坡度、植被覆盖度间存在较强的响应关系,建立回归函数如下:

其相关系数R2=0.980,模型拟合效果很好。研究发现公式(7)中的常数项同总产沙量的实测值比较,可以忽略,因此考虑修改期望函数,重新建立回归函数,其表达式如下:

其相关系数仍为R2=0.980,但形式更为简单。式中:Y为试验小区砒砂岩坡面的侵蚀产沙量(kg);I为降雨强度(mm/min);S为坡度值,取值范围10°～40°;C为植被覆盖度的百分数,取值范围0～50%;a,b,c,d,e为模型拟合参数。

表2是试验的实测值和公式(8)估计值的比较(图5)。从图5可以看到,即使侵蚀产沙量波动范围为4.20～509.78 kg,估测值与实际值的误差大部分都在10%以内,说明公式(8)的估计效果依然很好,且形式简单,能很好地反映变量间关系。

图5 总侵蚀产沙量值和估测值对比Fig.5 Comparison of total erosion sand production value and estimates

3.3 砒砂岩坡面侵蚀产沙量估测模型

通过对不同条件下砒砂岩坡面侵蚀的累计产沙量变化规律和总产沙量的数值模拟结果与分析,得到了试验小区条件下,侵蚀产沙量与降雨强度、坡度、植被覆盖度和降雨历时的相互关系,形成了如下估测模型,如公式(9)所示。

式中:Y为试验小区砒砂岩坡面的侵蚀产沙量(kg);I为降雨强度(mm/min);S为坡度值,取值范围10°～40°;C为植被覆盖度的百分数,取值范围0～50%;t为降雨历时(h),取值范围0～1。研究表明,坡面侵蚀产沙量的确与降雨强度间呈幂函数关系,而与坡度和植被覆盖度在一定范围内呈线性关系,且随降雨历时符合弧线式增长。

4 结论

(1)砒砂岩坡面产沙量受降雨强度的影响最为明显,受坡度的影响相对较弱。不同条件下,产沙量随降雨历时的变化曲线均为平缓型,初期产沙量随降雨历时的增加而不断增加,之后趋于稳定。

(2)植被的减蚀效果受降雨强度的影响很大,小降雨强度时,植被减蚀效果明显,然而随着降雨强度的增加,低植被覆盖度的效果越来越差,当降雨强度达到80 mm/h时,只有覆盖度为50%时,植被减沙效果明显。

(3)建立了砒砂岩坡面侵蚀产沙量的多元非线性回归模型。模型结果表明,产沙量同降雨强度呈幂函数关系,而同坡度和植被覆盖度在一定范围内呈线性关系,随降雨历时呈现很好的弧线式增长。然而,该模型目前对标准小区的特定条件下的拟合程度和适用性效果较好,针对野外侵蚀估测还需要更多的数据和模型修正,从而可以得到更加精确的预测模型。