砒砂巖坡面侵蝕產沙規律及多元回歸估算模型研究

梁止水,陳煜,孫悅,高海鷹,吳智仁

(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)建立了砒砂巖坡面侵蝕產沙量的多元非線性回歸模型。模型結果表明,產沙量同降雨強度呈冪函數關系,而同坡度和植被覆蓋度在一定范圍內呈線性關系,隨降雨歷時呈現很好的弧線式增長。然而,該模型目前對標準小區的特定條件下的擬合程度和適用性效果較好,針對野外侵蝕估測還需要更多的數據和模型修正,從而可以得到更加精確的預測模型。