砒砂巖坡面侵蝕產(chǎn)沙規(guī)律及多元回歸估算模型研究

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

(1.東南大學(xué)土木工程學(xué)院,南京 210096;2.江蘇大學(xué)環(huán)境與安全工程學(xué)院,江蘇鎮(zhèn)江 212013)

黃河流域砒砂巖區(qū)處于黃土高原到西北沙漠類(lèi)型的過(guò)渡區(qū),環(huán)境異質(zhì)性極為突出,生態(tài)退化和水土流失非常嚴(yán)重,侵蝕模數(shù)可達(dá)到30 000～40 000 t/(km2·a),是我國(guó)乃至世界上侵蝕最為劇烈的地區(qū)之一,也是黃河下游“地上懸河”的粗泥沙來(lái)源核心區(qū)[1]。長(zhǎng)期以來(lái),砒砂巖區(qū)被列為國(guó)家生態(tài)環(huán)境建設(shè)和水土流失治理的重點(diǎn)區(qū),水土流失防治取得了較為顯著的治理成效[2-4]。

針對(duì)砒砂巖的水力侵蝕,主要集中在沖刷條件下的坡面侵蝕規(guī)律和水力學(xué)特征。其中:蘇濤等[5-6]利用室內(nèi)和野外徑流沖刷試驗(yàn)研究了砒砂巖坡面的侵蝕規(guī)律和徑流水動(dòng)力學(xué)特性。楊吉山等[7]通過(guò)野外放水沖刷試驗(yàn)對(duì)比研究了白色和紅色原狀砒砂巖坡面的產(chǎn)流、產(chǎn)沙過(guò)程,分析了砒砂巖剝蝕率與水力學(xué)參數(shù)之間的相似關(guān)系。Liang等[8]采用室內(nèi)模擬降雨試驗(yàn)分析了不同降雨強(qiáng)度、坡度、植被覆蓋率條件下的產(chǎn)沙規(guī)律。楊振奇等[9]運(yùn)用系統(tǒng)聚類(lèi)和線(xiàn)性回歸相結(jié)合的方法,劃分裸露砒砂巖區(qū)降雨類(lèi)型,研究了不同降雨類(lèi)型對(duì)裸露坡面產(chǎn)流和產(chǎn)沙的影響,分析了不同降雨條件下坡面的微地形變化規(guī)律。李龍等[10]研究了自然降雨條件下砒砂巖坡面細(xì)溝微形態(tài)及其侵蝕特征,分析了細(xì)溝微形態(tài)變化過(guò)程對(duì)產(chǎn)流產(chǎn)沙的影響。董曉宇等[11]采用野外原位模擬沖刷試驗(yàn)研究了裸露砒砂巖區(qū)坡面侵蝕過(guò)程中地表粗糙度與水力侵蝕特征參數(shù)的關(guān)系。然而針對(duì)砒砂巖坡面的水土流失模擬及侵蝕產(chǎn)沙預(yù)測(cè)模型等研究較少,其中只有葉俊道等[12]對(duì)WEPP模型在砒砂巖地區(qū)不同坡面的土壤侵蝕進(jìn)行了適用性探討。因此,亟需針對(duì)砒砂巖坡面的特點(diǎn)建立適用于該地區(qū)的產(chǎn)沙過(guò)程預(yù)報(bào)模型。

為了更好地研究砒砂巖坡面在降雨侵蝕作用下侵蝕產(chǎn)沙規(guī)律,本文開(kāi)展室內(nèi)模擬降雨試驗(yàn),研究不同砒砂巖邊坡的坡度、降雨強(qiáng)度及植被覆蓋率條件下的侵蝕產(chǎn)沙量,利用多元回歸分析方法,分析不同參數(shù)與侵蝕產(chǎn)沙量之間的相互關(guān)系,初步建立侵蝕產(chǎn)沙估測(cè)模型,以期為砒砂巖水土流失區(qū)治理工作提供一定的科學(xué)依據(jù)。

1 材料及方法

1.1 試驗(yàn)材料

試驗(yàn)用砒砂巖土取自黃河中游上游皇甫川流域的圪秋溝,含水量約為9.0%,土壤顆粒級(jí)配使用馬爾文3000激光粒度儀測(cè)定,＞1.00 mm 的顆粒組成占0.83%,0.50～0.05 mm 的顆粒組成占90.74%,＜0.01 mm 的顆粒組成占8.43%。所選用的植物為狗牙根〔Cynodondactylon(L.)Persoon〕,它與砒砂巖區(qū)廣泛生長(zhǎng)的野牛草極為類(lèi)似,也具有根系較小、葉面窄長(zhǎng)、移植成活率高等特點(diǎn)。

1.2 試驗(yàn)設(shè)計(jì)

采用標(biāo)準(zhǔn)小區(qū)5 m×1 m,進(jìn)行室內(nèi)模擬降雨試驗(yàn),主要考察降雨強(qiáng)度、坡度和植被覆蓋度的影響。根據(jù)全面調(diào)查立體條件和氣候條件的基礎(chǔ)上,砒砂巖區(qū)相對(duì)穩(wěn)定坡度多在30°左右[10]。因此,試驗(yàn)中考察砒砂巖的坡度為20°,30°,40°。根據(jù)當(dāng)?shù)囟嗄杲涤陻?shù)據(jù)資料選擇降雨強(qiáng)度為20,50,80 mm/h,植被覆蓋率為10%,30%,50%。此外,在砒砂巖區(qū)溝坡系統(tǒng)的陡坡面上很少有植被生長(zhǎng),因此,僅在20°的緩坡條件下進(jìn)行植被覆蓋度的影響研究。試驗(yàn)設(shè)計(jì)見(jiàn)表1,砒砂巖裸地坡面依次編號(hào)為A1—A9,有植被覆蓋坡面依次編號(hào)為B1—B9。

表1 試驗(yàn)設(shè)計(jì)方案Table 1 Design of experiments

1.3 試驗(yàn)過(guò)程與測(cè)試方法

模擬降雨試驗(yàn)在黃河水利科學(xué)研究院模型黃河降雨試驗(yàn)大廳進(jìn)行,其中模擬降雨裝置采用下噴式的自動(dòng)模擬降雨系統(tǒng),其降雨的均勻性可達(dá)到90%以上。試驗(yàn)用的土槽為移動(dòng)式可變坡度的鋼槽,其坡度調(diào)節(jié)范圍為5°～45°,土槽的長(zhǎng)為5.0 m,寬為2.0 m,高為0.6 m,寬度方向上分為2個(gè)1.0 m 的土槽。狗牙草則采用試驗(yàn)前一個(gè)月的苗子進(jìn)行移植,并通過(guò)對(duì)種植時(shí)的密度控制及葉面修剪來(lái)控制覆蓋度大小,從而保證試驗(yàn)中的植被覆蓋率要求。

模擬砒砂巖坡面采用控制容重的方法填土,層層壓實(shí),容重均控制在1.4～1.5 g/cm3,使其接近野外實(shí)際砒砂巖容重。為便于滲水,填土之前,在坡面底部預(yù)先鋪上10 cm 厚的沙子。每場(chǎng)降雨試驗(yàn)的前一天,先用非常小的雨量進(jìn)行預(yù)降雨,在基本保證不產(chǎn)生侵蝕的條件下直至形成產(chǎn)流為止,從而可以有效消除不同工況下的邊界差異。同時(shí),在每次降雨前,均需要用塑料薄膜將坡面蓋住,并在坡面的4個(gè)角放置雨量筒,初始降雨15 min,用于校核降雨強(qiáng)度大小及其均勻性。降雨持續(xù)1 h,降雨過(guò)程中每3 min接一個(gè)產(chǎn)流泥沙樣,采用集沙桶收集產(chǎn)流產(chǎn)沙情況,后用烘干法測(cè)得產(chǎn)沙量。

1.4 數(shù)據(jù)分析和建模方法

模擬試驗(yàn)結(jié)果利用SPSS軟件和非線(xiàn)性回歸分析方法進(jìn)行擬合,建立回歸方程。

2 結(jié)果與分析

2.1 砒砂巖裸地坡面產(chǎn)沙特征分析

從圖1中可以看出,不同坡度條件下,降雨初期單位時(shí)間的產(chǎn)沙量均隨降雨歷時(shí)的增加不斷增加,隨后逐漸趨于穩(wěn)定。其中,降雨強(qiáng)度越小的情況下,單位時(shí)間內(nèi)的產(chǎn)沙量最終趨于穩(wěn)定的時(shí)間越早。而在降雨歷時(shí)半小時(shí)以后,所有工況下的單位侵蝕產(chǎn)沙量基本都趨于穩(wěn)定。已有研究表明[13],黃土類(lèi)坡面的侵蝕產(chǎn)沙量隨降雨歷時(shí)的變化曲線(xiàn)會(huì)有3種形式:平緩型、單峰型和多峰型,而砒砂巖坡面侵蝕產(chǎn)沙量隨時(shí)間變化均為平緩型。在本次試驗(yàn)研究的時(shí)間范圍內(nèi)尚未出現(xiàn)峰值,這其中可能的原因是砒砂巖在遇水后短時(shí)間內(nèi)出現(xiàn)侵蝕潰散,并呈現(xiàn)出侵蝕惡化的趨勢(shì),直至達(dá)到水流可攜帶泥沙的極大值。另外,侵蝕產(chǎn)沙量受降雨強(qiáng)度的影響為極其顯著(p＜0.01),即為降雨強(qiáng)度越大,砒砂巖坡面形成的產(chǎn)沙量也越大,曲線(xiàn)波動(dòng)的程度也更加明顯。而在試驗(yàn)中的3種降雨強(qiáng)度條件下,侵蝕產(chǎn)沙量均隨坡度的增加而呈現(xiàn)輕微的增加趨勢(shì),但該趨勢(shì)并不明顯(p＞0.05),說(shuō)明坡度對(duì)侵蝕產(chǎn)沙量的影響小于降雨強(qiáng)度的影響[14]。

圖1 裸坡坡面產(chǎn)沙量隨降雨歷時(shí)的變化曲線(xiàn)Fig.1 Curve of sand yield on the bare slope with the duration of rainfall

2.2 砒砂巖植被覆蓋坡面產(chǎn)沙特征分析

圖2為不同降雨強(qiáng)度和植被覆蓋度下,砒砂巖坡面產(chǎn)沙量隨時(shí)間的變化曲線(xiàn),從圖中可以看出,砒砂巖植被覆蓋坡面產(chǎn)沙量隨時(shí)間變化規(guī)律同裸露坡面的情況相似,其中初期產(chǎn)沙量隨降雨歷時(shí)的增加而不斷增大,之后逐漸趨于穩(wěn)定。砒砂巖坡面植被減蝕效果也明顯受降雨強(qiáng)度影響,不同降雨強(qiáng)度時(shí),植被的減蝕率相差較大。當(dāng)降雨強(qiáng)度為20 mm/h時(shí),不同覆蓋度的植被減沙效果均非常明顯(p＜0.01),所有植被覆蓋度的減沙率均超過(guò)70.00%。但當(dāng)降雨強(qiáng)度為80 mm/h時(shí),只有植被覆蓋度為50%對(duì)應(yīng)的產(chǎn)沙量與裸坡相比變化明顯,減沙率為86.72%,其他覆蓋度條件,產(chǎn)沙量與裸坡相比變化并不明顯(p＞0.05)。

圖2 植被覆蓋坡面產(chǎn)沙量隨降雨歷時(shí)的變化曲線(xiàn)Fig.2 Curve of sediment yield on vegetation cover slope with rainfall duration

3 討論

根據(jù)模擬降雨試驗(yàn)結(jié)果,得到了累計(jì)產(chǎn)沙量隨降雨時(shí)間的關(guān)系(圖3),可以看出,不同坡面的累計(jì)產(chǎn)沙量隨降雨歷時(shí)增加而增加,且有增加變陡的趨勢(shì)(曲率增加),這是由于隨著時(shí)間推移,坡面侵蝕程度增大,徑流中所含的泥沙在輸移過(guò)程中增強(qiáng)了侵蝕動(dòng)力,加劇了侵蝕產(chǎn)沙,該結(jié)果與肖培青等[15]得出累計(jì)產(chǎn)沙量與降雨歷時(shí)成冪函數(shù)的結(jié)論基本一致。

圖3 累計(jì)產(chǎn)沙量隨降雨歷時(shí)的關(guān)系曲線(xiàn)Fig.3 Relationship curve of cumulative sediment production with rainfall duration

3.1 累計(jì)產(chǎn)沙量擬合曲線(xiàn)

由圖3看出,在不同的試驗(yàn)工況下,每分鐘產(chǎn)沙量相差很大,導(dǎo)致這些累計(jì)產(chǎn)沙與降雨歷時(shí)的關(guān)系曲線(xiàn)偏離也很大,難以統(tǒng)一分析。如果采用冪函數(shù)擬合,每一個(gè)試驗(yàn)條件都得到一個(gè)具有兩個(gè)參數(shù)的公式,得不到統(tǒng)一的對(duì)應(yīng)關(guān)系。因此,考慮對(duì)累計(jì)產(chǎn)沙量進(jìn)行無(wú)量綱的歸一化處理,使得數(shù)值在0～1[16],即為:將累計(jì)產(chǎn)沙量值都除以其1 h總產(chǎn)沙量,降雨歷時(shí)以1 h的百分?jǐn)?shù)表示。并得到了歸一化后的關(guān)系曲線(xiàn)(圖4),所有工況條件下得到曲線(xiàn)重合性好。

圖4 累計(jì)產(chǎn)沙量隨降雨歷時(shí)變化的歸一化關(guān)系曲線(xiàn)Fig.4 Normalized relationship curve of cumulative sediment production change with rainfall duration

由圖4可以看出,歸一化曲線(xiàn)基本都位于y=x下方,且非常接近以其為玄的圓弧,因此考慮采用曲線(xiàn)擬合,并假定為一條圓弧,圓心設(shè)為(a,b),在曲線(xiàn)的左上方,半徑為r,其中圓弧同時(shí)經(jīng)過(guò)點(diǎn)(0,0)和(1,1),則該弧線(xiàn)上的點(diǎn)(x,y)滿(mǎn)足以下關(guān)系:

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

從而獲得目標(biāo)的擬合方程為:

對(duì)每條歸一化后的曲線(xiàn)按照公式(2)進(jìn)行非線(xiàn)性擬合,得到a值和擬合相關(guān)系數(shù)R2(表2),其中相關(guān)系數(shù)均在0.95以上,因此擬合的誤差較小。而B(niǎo)1,B2和B9這3個(gè)組別中a值特別大,說(shuō)明該3種工況下,圓弧曲率很大,曲線(xiàn)接近于直線(xiàn),且都屬于具有覆蓋率的兩個(gè)邊緣情況,說(shuō)明該擬合曲線(xiàn)存在一定的臨界條件。在試驗(yàn)組的基礎(chǔ)上,去除該三組數(shù)據(jù)后,對(duì)其他數(shù)據(jù)進(jìn)行求平均值,從而保證其統(tǒng)一擬合的效果,從而得到a值約為-2.480,方差為3.154,各組數(shù)據(jù)基本均勻分布在擬合曲線(xiàn)的附近,且符合歸一化曲線(xiàn)下凹、斜率逐漸增大的特點(diǎn),可以較好地反映試驗(yàn)條件下累計(jì)產(chǎn)沙量隨降雨歷時(shí)的變化規(guī)律,故可以得到砒砂巖坡面侵蝕的累計(jì)產(chǎn)沙量與降雨歷時(shí)(t)的擬合曲線(xiàn)方程為:

表2 擬合參數(shù)a 值及相關(guān)系數(shù)Table 2 Fitting parameter a values and correlation coefficients

式中:t為降雨歷時(shí)(h),取值范圍0～1。

3.2 多元非線(xiàn)性回歸分析模型

多元非線(xiàn)性分析是利用數(shù)理統(tǒng)計(jì)的方法建立多個(gè)自變量與一個(gè)因變量之間非線(xiàn)性的函數(shù)關(guān)系,在諸多領(lǐng)域得到廣泛應(yīng)用[17-19]。

(1)傳統(tǒng)方法是采用曲線(xiàn)化直線(xiàn),利用6種初等函數(shù)對(duì)自變量進(jìn)行變換,然后根據(jù)經(jīng)典的最小二乘法原理,建立預(yù)測(cè)因變量與變換后自變量的線(xiàn)性關(guān)系。這種模型有一定精度,但沒(méi)有考慮因素間的交互影響。其數(shù)學(xué)表達(dá)式一般為:

(2)另一種有效的方法是直接提出非直線(xiàn)形式的期望函數(shù),作為模型框架,采用高斯-牛頓法進(jìn)行參數(shù)估計(jì),反復(fù)迭代直至參數(shù)估計(jì)收斂。該方法的方差分析F值和p值的意義不大。其數(shù)學(xué)表達(dá)式一般為:

影響坡面水力侵蝕產(chǎn)沙的因素,主要包括降雨、坡度、土壤特性(含水量、粒徑大小、黏粒含量、內(nèi)摩擦角等)、植被,以及坡長(zhǎng)和人為作用等因素[20]。試驗(yàn)條件下,降雨強(qiáng)度、坡度和植被覆蓋度是影響產(chǎn)沙量的3個(gè)最主要因素。大量的研究表明,侵蝕產(chǎn)沙量與降雨強(qiáng)度間存在冪函數(shù)關(guān)系[21-22],并在該研究中得到了驗(yàn)證。對(duì)于具有臨界坡度的情況下,在臨界坡度之前,與坡度間可能存在二次函數(shù)關(guān)系[23];而與植被覆蓋度間的關(guān)系比較復(fù)雜,尚未明確,為方便之后建模,暫定為線(xiàn)性關(guān)系。預(yù)測(cè)模型的因變量為總侵蝕產(chǎn)沙量Y(單位kg),自變量分別為:降雨強(qiáng)度I(為方便計(jì)算,單位選擇mm/min);坡面坡度數(shù)值S;植被覆蓋度的百分?jǐn)?shù)C。

采用第一種方法(單因素分析,表明Y與I2.5呈線(xiàn)性關(guān)系),可得到回歸函數(shù)如下:

其相關(guān)系數(shù)R2=0.839,擬合效果一般,且各參數(shù)的標(biāo)準(zhǔn)誤差均很大,用作回歸模型不太理想。

采用第二種方法,需預(yù)先確定期望函數(shù)。根據(jù)前面對(duì)不同條件下產(chǎn)沙量變化過(guò)程的分析,降雨強(qiáng)度是主因,且自變量中降雨強(qiáng)度同坡度、植被覆蓋度間存在較強(qiáng)的響應(yīng)關(guān)系,建立回歸函數(shù)如下:

其相關(guān)系數(shù)R2=0.980,模型擬合效果很好。研究發(fā)現(xiàn)公式(7)中的常數(shù)項(xiàng)同總產(chǎn)沙量的實(shí)測(cè)值比較,可以忽略,因此考慮修改期望函數(shù),重新建立回歸函數(shù),其表達(dá)式如下:

其相關(guān)系數(shù)仍為R2=0.980,但形式更為簡(jiǎn)單。式中:Y為試驗(yàn)小區(qū)砒砂巖坡面的侵蝕產(chǎn)沙量(kg);I為降雨強(qiáng)度(mm/min);S為坡度值,取值范圍10°～40°;C為植被覆蓋度的百分?jǐn)?shù),取值范圍0～50%;a,b,c,d,e為模型擬合參數(shù)。

表2是試驗(yàn)的實(shí)測(cè)值和公式(8)估計(jì)值的比較(圖5)。從圖5可以看到,即使侵蝕產(chǎn)沙量波動(dòng)范圍為4.20～509.78 kg,估測(cè)值與實(shí)際值的誤差大部分都在10%以?xún)?nèi),說(shuō)明公式(8)的估計(jì)效果依然很好,且形式簡(jiǎn)單,能很好地反映變量間關(guān)系。

圖5 總侵蝕產(chǎn)沙量值和估測(cè)值對(duì)比Fig.5 Comparison of total erosion sand production value and estimates

3.3 砒砂巖坡面侵蝕產(chǎn)沙量估測(cè)模型

通過(guò)對(duì)不同條件下砒砂巖坡面侵蝕的累計(jì)產(chǎn)沙量變化規(guī)律和總產(chǎn)沙量的數(shù)值模擬結(jié)果與分析,得到了試驗(yàn)小區(qū)條件下,侵蝕產(chǎn)沙量與降雨強(qiáng)度、坡度、植被覆蓋度和降雨歷時(shí)的相互關(guān)系,形成了如下估測(cè)模型,如公式(9)所示。

式中:Y為試驗(yàn)小區(qū)砒砂巖坡面的侵蝕產(chǎn)沙量(kg);I為降雨強(qiáng)度(mm/min);S為坡度值,取值范圍10°～40°;C為植被覆蓋度的百分?jǐn)?shù),取值范圍0～50%;t為降雨歷時(shí)(h),取值范圍0～1。研究表明,坡面侵蝕產(chǎn)沙量的確與降雨強(qiáng)度間呈冪函數(shù)關(guān)系,而與坡度和植被覆蓋度在一定范圍內(nèi)呈線(xiàn)性關(guān)系,且隨降雨歷時(shí)符合弧線(xiàn)式增長(zhǎng)。

4 結(jié)論

(1)砒砂巖坡面產(chǎn)沙量受降雨強(qiáng)度的影響最為明顯,受坡度的影響相對(duì)較弱。不同條件下,產(chǎn)沙量隨降雨歷時(shí)的變化曲線(xiàn)均為平緩型,初期產(chǎn)沙量隨降雨歷時(shí)的增加而不斷增加,之后趨于穩(wěn)定。

(2)植被的減蝕效果受降雨強(qiáng)度的影響很大,小降雨強(qiáng)度時(shí),植被減蝕效果明顯,然而隨著降雨強(qiáng)度的增加,低植被覆蓋度的效果越來(lái)越差,當(dāng)降雨強(qiáng)度達(dá)到80 mm/h時(shí),只有覆蓋度為50%時(shí),植被減沙效果明顯。

(3)建立了砒砂巖坡面侵蝕產(chǎn)沙量的多元非線(xiàn)性回歸模型。模型結(jié)果表明,產(chǎn)沙量同降雨強(qiáng)度呈冪函數(shù)關(guān)系,而同坡度和植被覆蓋度在一定范圍內(nèi)呈線(xiàn)性關(guān)系,隨降雨歷時(shí)呈現(xiàn)很好的弧線(xiàn)式增長(zhǎng)。然而,該模型目前對(duì)標(biāo)準(zhǔn)小區(qū)的特定條件下的擬合程度和適用性效果較好,針對(duì)野外侵蝕估測(cè)還需要更多的數(shù)據(jù)和模型修正,從而可以得到更加精確的預(yù)測(cè)模型。