金 鑫，李 霞，宋 穎，王春振，趙華榮

(1.桂林理工大學(xué)廣西環(huán)境污染控制理論與技術(shù)重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541004；2.桂林理工大學(xué)巖溶地區(qū)水污染控制與用水安全保障協(xié)同創(chuàng)新中心，廣西 桂林 541004)

0 引 言

土壤侵蝕是多種因素共同作用的復(fù)雜過程，降雨強(qiáng)度、坡度、降雨歷時(shí)、土壤初始含水率、抗剪切力作為其中的重要影響因素而得到研究者的廣泛關(guān)注。已有的研究表明：坡度和降雨強(qiáng)度均為坡面侵蝕的主要影響因子，隨著坡度和雨強(qiáng)增加，徑流率和產(chǎn)沙量也增加，雨強(qiáng)是徑流率變化的主要因素，坡度對(duì)產(chǎn)沙量變化的影響較為顯著，且存在臨界坡度，產(chǎn)流和產(chǎn)沙隨降雨歷時(shí)先迅速增加后趨于穩(wěn)定[1-4]；初始含水率越高，產(chǎn)流越快，達(dá)到穩(wěn)定入滲率的時(shí)間也越短[5]；與坡度、土壤含水率相比，土壤剪切力對(duì)侵蝕量影響最小，但其能夠顯著提高土壤侵蝕量預(yù)測(cè)精度[6]。

在雨強(qiáng)、坡度、降雨歷時(shí)、初始含水率、剪切力等因素對(duì)土壤侵蝕影響方面已經(jīng)不少學(xué)者分別做了研究，但很少細(xì)化這些因素對(duì)侵蝕作用的程度。在實(shí)際侵蝕過程中，這些因素同時(shí)作用于侵蝕過程，其對(duì)侵蝕的貢獻(xiàn)大小以及與侵蝕產(chǎn)沙的關(guān)系需要進(jìn)行量化，以準(zhǔn)確地反映各因素對(duì)侵蝕過程的影響程度。

偏最小二乘回歸(PLS)是一種新的多元統(tǒng)計(jì)數(shù)據(jù)分析方法，將多元線性回歸、典型相關(guān)分析和主成分分析進(jìn)行有機(jī)結(jié)合，可避免自變量、因變量之間存在的多重線性關(guān)系，具有更大的優(yōu)勢(shì)，從而使模型精度、穩(wěn)健性、實(shí)用性都得到提高[7]。并在大尺度的流域徑流和土壤侵蝕預(yù)測(cè)中得到較好的應(yīng)用[8-11]。本文采用偏最小二乘法對(duì)黃土坡面土壤侵蝕進(jìn)行模擬和預(yù)測(cè)，細(xì)化分析雨強(qiáng)、坡度、初始含水率、剪切力、降雨歷時(shí)等因素對(duì)黃土坡面產(chǎn)流產(chǎn)沙過程的影響程度，以期為認(rèn)識(shí)土壤侵蝕過程和水土流失防治提供數(shù)據(jù)支持和基礎(chǔ)見解。

1 材料和方法

1.1 實(shí)驗(yàn)裝置與材料

本實(shí)驗(yàn)在桂林理工大學(xué)降雨實(shí)驗(yàn)大廳進(jìn)行，實(shí)驗(yàn)所用土壤為取自陜西岔巴溝流域的表層土。實(shí)驗(yàn)采用人工模擬坡面降雨，降雨裝置由土槽與人工模擬降雨器2個(gè)部分組成。實(shí)驗(yàn)土槽為鋼質(zhì)，尺寸長(zhǎng)×寬×高為：400 cm×120 cm×80 cm，土槽坡度可以自動(dòng)調(diào)節(jié)，變化范圍為0～30 °，土槽尾部設(shè)有簸箕型集水口，用來收集坡面徑流，底部有集水口用以收集坡面下滲徑流；人工模擬降雨器為頂噴式全自動(dòng)不銹鋼模擬降雨器，配有旋轉(zhuǎn)下噴式噴頭，為帶壓力垂直下噴式降雨，用揚(yáng)程為50 m的潛水泵供水，降雨高度為6 m，使得雨滴散落在坡面的速度接近自然雨滴終速，雨強(qiáng)為10～200 mm/h，雨滴直徑為0.1～6 mm，雨滴的中值粒徑及降雨動(dòng)能與自然降雨條件較為接近[12]。模擬降雨實(shí)驗(yàn)之前測(cè)定降雨均勻性，測(cè)定結(jié)果顯示降雨均勻性大于80%。綜上所述，實(shí)驗(yàn)裝置設(shè)計(jì)合理。

1.2 實(shí)驗(yàn)和測(cè)量

由于黃土高原降雨分布比較集中，多以暴雨形式出現(xiàn)，且根據(jù)調(diào)查研究，5°以上，細(xì)溝侵蝕較強(qiáng)，并發(fā)生淺溝侵蝕，15°以上，細(xì)溝、淺溝侵蝕強(qiáng)烈[13]。故本實(shí)驗(yàn)設(shè)置了50、60、70、80 mm/h四個(gè)降雨強(qiáng)度，5°、10°、15°、20°四個(gè)坡度，共進(jìn)行16場(chǎng)降雨實(shí)驗(yàn)，每場(chǎng)實(shí)驗(yàn)每場(chǎng)降雨時(shí)長(zhǎng)為15或30 min。每次降雨開始前，測(cè)定坡面土壤的初始含水率、土壤剪切力、容重。土壤初始含水率通過TDR水分測(cè)定儀測(cè)定，土壤剪切力通過三頭抗剪儀測(cè)定，容重采用環(huán)刀法測(cè)定，每個(gè)參數(shù)的測(cè)定結(jié)果均為多次取樣的平均值。

降雨開始后，記錄降雨開始時(shí)間。坡面產(chǎn)流后，記錄產(chǎn)流時(shí)間，并根據(jù)實(shí)際產(chǎn)流及產(chǎn)沙情況每間隔0.5～4 min采一次樣，直到降雨結(jié)束。用量筒收集徑流樣品，秒表記錄相應(yīng)的采集時(shí)間，記錄樣品體積并進(jìn)行稱重，通過計(jì)算得到各個(gè)時(shí)段的徑流率。

樣品含沙量的測(cè)量方法采用烘干稱重法，具體步驟為：采集的樣品靜置沉淀24 h以上，抽出上清液，將剩余泥沙樣倒入燒杯中，放入105 ℃的烘箱中24 h，冷卻后取出稱重，通過計(jì)算得到各個(gè)時(shí)段的產(chǎn)沙率。

1.3 PLS建模

PLS的原理可表達(dá)[14]為：設(shè)有q個(gè)因變量{y1，y2，…，yq}和p個(gè)自變量{x1，x2，…，xp}，為了研究因變量與自變量的統(tǒng)計(jì)關(guān)系，觀測(cè)了n個(gè)樣本點(diǎn)，由此構(gòu)成了自變量與因變量的數(shù)據(jù)表X=(x1，x2，…，xp)n×p和Y=(y1，y2…，yq)n×q，偏最小二乘回歸分別在X與Y中提取出成分t1和u1(也就是說，t1是x1，x2，…，xp的線性組合，u1是y1，y2，…，yq的線性組合)，t1和u1應(yīng)盡可能好地代表數(shù)據(jù)表X和Y，同時(shí)，自變量的成分t1對(duì)因變量的成分u1又有很強(qiáng)的解釋能力，在第1個(gè)成分t1和u1被提取后，偏最小二乘回歸分別實(shí)施X對(duì)t1的回歸以及Y對(duì)t1的回歸，如果回歸方程已經(jīng)達(dá)到滿意的精度，則算法終止；否則，將利用X被t1解釋后的殘余信息以及Y被t1解釋后的殘余信息進(jìn)行第2輪的成分提取，如此反復(fù)，直到能達(dá)到一個(gè)較滿意的精度為止，若最終對(duì)X共提取了m個(gè)成分t1，t2，…，tm，偏最小二乘回歸將通過施行yk(k=1，2，…，q)對(duì)t1，t2，…，tm的回歸，然后表達(dá)成yk關(guān)于原變量x1，x2，…，xp的回歸方程。

PLS建模過程依次為[15]：

(1)數(shù)據(jù)標(biāo)準(zhǔn)化處理；目的是使樣本點(diǎn)的集合重心與坐標(biāo)原點(diǎn)重合，先將X經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為E0=(E01，…，E0p)n×p，再將Y經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為F0=(F01，…，F(xiàn)0q)n×p；

(4)第k成分tk提取；采用步驟(2)及步驟(3)的方法，依次進(jìn)行推求并對(duì)第k成分tk提取，k可用交叉有效性原則進(jìn)行識(shí)別；

2 結(jié)果與分析

2.1 參數(shù)選取

通過對(duì)影響因素的分析并結(jié)合相關(guān)研究，將影響坡面徑流率(Y1)和產(chǎn)沙率(Y2)的主要因素確定為：坡度(X1)、雨強(qiáng)(X2)、時(shí)間(X3)、初始含水率(X4)、剪切力(X5)。選用5°、10°、15°、20°四個(gè)坡度，50、60、70、80 mm/h 4個(gè)降雨強(qiáng)度在不同的初始含水率、剪切力的16場(chǎng)降雨，其中14場(chǎng)0～30 min降雨過程中的172個(gè)徑流率、產(chǎn)沙率數(shù)據(jù)作為建模輸入，15°、70 mm/h，20°、80 mm/h的2場(chǎng)降雨在0～15 min的19個(gè)徑流率、產(chǎn)沙率數(shù)據(jù)進(jìn)行檢驗(yàn)分析。

表1 降雨實(shí)驗(yàn)影響因素表

2.2 偏最小二乘模型

SIMCA-P是基于主成分分析、偏最小二乘回歸分析開發(fā)的數(shù)據(jù)分析軟件，因其方便有效的特點(diǎn)，在各個(gè)領(lǐng)域獲得了廣泛應(yīng)用[7]。本研究通過SIMCA軟件構(gòu)建徑流率、產(chǎn)沙率的偏最小二乘模型。

初步假設(shè)坡度(X1)、雨強(qiáng)(X2)、時(shí)間(X3)、初始含水率(X4)、剪切力(X5)和徑流率(Y1)和產(chǎn)沙率(Y2)的關(guān)系為線性關(guān)系，直接導(dǎo)入數(shù)據(jù)構(gòu)建模型。模型的擬合優(yōu)度R2和預(yù)測(cè)優(yōu)度Q2都大于0.6，達(dá)到了可以擬合的標(biāo)準(zhǔn)。同時(shí)假設(shè)坡度(X1)、雨強(qiáng)(X2)、時(shí)間(X3)、初始含水率(X4)、剪切力(X5)、徑流率(Y1)和產(chǎn)沙率(Y2)的關(guān)系為非線性指數(shù)疊加關(guān)系，通過對(duì)數(shù)轉(zhuǎn)化重新構(gòu)建偏最小二乘模型，模型的R2和Q2分別為0.843、0.831，相比較直接線性模擬，指數(shù)疊加關(guān)系經(jīng)對(duì)數(shù)化轉(zhuǎn)變后，模型擬合的效果要更好，故下文采用對(duì)數(shù)轉(zhuǎn)化后的偏最小二乘模型結(jié)果進(jìn)行論述。

表2為自變量、因變量間的相關(guān)系數(shù)矩陣。從表2可以看到，自變量間雨強(qiáng)、初始含水率、剪切力之間存在一定的相關(guān)性，兩個(gè)因變量之間的相關(guān)性也很強(qiáng)，達(dá)到了0.844。偏最小二乘模型可以有效的解決多重共線性問題，同時(shí)計(jì)算交叉有效性，該模型取3個(gè)自變量便達(dá)到了滿意的精度，則提取3個(gè)自變量建模。采用SIMCA軟件，通過分析得到的偏最小二乘回歸標(biāo)準(zhǔn)化方程為：

表2 變量之間的相關(guān)系數(shù)矩陣

lnY1=0.039 493 4×lnX1+0.408 309×lnX2+

0.824 303×lnX3-0.128 515×lnX4+0.145 553×lnX5

lnY2=0.412 347×lnX1+0.341 94×lnX2+

0.706 01×lnX3-0.042 129 9×lnX4-0.023 738 6×lnX5

同時(shí)得到了原始的回歸方程：

lnY1=0.070 650 6×lnX1+2.014 59×lnX2+

0.810 188×lnX3-1.444 52×lnX4+0.596 852×lnX5-1.167 5

lnY2=1.362 91×lnX1+3.117 17×lnX2+1.282 1×

lnX3-0.874 926×lnX4-0.179 851×lnX5-12.47

即：

變量投影重要性指標(biāo)VIP，是描述自變量對(duì)因變量影響程度的指標(biāo)，若某個(gè)變量的VIP>1，則認(rèn)為該變量是顯著影響因素，值越大說明該自變量對(duì)因變量的影響越大[15]。為分析各個(gè)自變量對(duì)因變量的影響，繪制了變量投影圖。如圖1所示，采用VIP值來描述每一個(gè)自變量對(duì)因變量的作用大小。從圖1可以看出，VIP的順序?yàn)椋篨3>X1>X2>X4>X5。根據(jù)VIP>1則X在解釋因變量時(shí)具有重要作用的原則，時(shí)間X3(min)、坡度X1(°)解釋因變量徑流率Y1(mL/s)和產(chǎn)沙率Y2(g/s)時(shí)具有重要作用，其中時(shí)間X3(min)在解釋因變量Y1、Y2時(shí)具有最重要的作用。

為了更直觀的觀測(cè)5個(gè)自變量對(duì)2個(gè)因變量的具體作用，針對(duì)標(biāo)準(zhǔn)化數(shù)據(jù)的回歸方程，繪制了回歸系數(shù)圖(見圖2)。從圖2中可以看出，初始含水率對(duì)徑流率的影響是要高于對(duì)產(chǎn)沙率的影響。但無論是對(duì)于產(chǎn)沙率還是徑流率而言，初始含水率的影響都是負(fù)值，初始含水率越高，坡面的徑流率和產(chǎn)沙率反而減小，這與文獻(xiàn)[5]的研究結(jié)果不相一致，這可能是由于各場(chǎng)實(shí)驗(yàn)初始含水率本身差值不夠明顯，開始降雨后極短時(shí)間內(nèi)土壤含水率就趨于飽和，難以避免的實(shí)驗(yàn)和觀測(cè)誤差導(dǎo)致模型對(duì)影響程度較低的含水率計(jì)算不夠準(zhǔn)確，從而出現(xiàn)隨著初始含水率增加，徑流率和產(chǎn)沙率總體為減小趨勢(shì)這一現(xiàn)象?？辜羟辛?duì)徑流率和產(chǎn)沙率的影響比較有限，對(duì)徑流率的影響為正值，對(duì)產(chǎn)沙率為負(fù)值，這表明土壤抗剪切性越強(qiáng)，土壤越不易被侵蝕，坡面比較平整，徑流受到的阻礙會(huì)越小。

圖2 回歸系數(shù)圖

坡度的變化主要是對(duì)產(chǎn)沙率產(chǎn)生影響，對(duì)于徑流率，坡度變化的影響是微乎其微的，這與ZHAO[16]等人的研究結(jié)果一致，推測(cè)這種情況的出現(xiàn)可能由于在實(shí)驗(yàn)的坡度范圍存在著臨界坡度，使得徑流率在這個(gè)區(qū)間的變化呈現(xiàn)波動(dòng)，導(dǎo)致模型中無法精確體現(xiàn)出坡度對(duì)徑流率的影響。降雨強(qiáng)度是徑流率變化的一個(gè)主要影響因素，這是因?yàn)榻涤晔菑搅鳟a(chǎn)生的來源，徑流率的變化與雨強(qiáng)的變化緊密相關(guān)。在降雨初期，由于土壤入滲未達(dá)到穩(wěn)定，徑流是逐漸增加的，產(chǎn)沙也相伴增加，二者隨降雨歷時(shí)的變化明顯。產(chǎn)沙率的主要影響因素除時(shí)間之外，坡度和降雨強(qiáng)度的影響相差較小，受到二者共同作用比較顯著，坡度變化，土體的勢(shì)能和水流動(dòng)能都發(fā)生改變，而降雨強(qiáng)度的變化會(huì)使得動(dòng)能的改變進(jìn)一步加劇，產(chǎn)生疊加效應(yīng)，二者交互作用共同影響產(chǎn)沙率的變化。時(shí)間對(duì)徑流率和產(chǎn)沙率變化的影響程度相差不大，徑流率和產(chǎn)沙率都隨產(chǎn)流后降雨時(shí)間的延長(zhǎng)而增大。

2.3 精度分析

利用構(gòu)建的偏最小二乘回歸模型，將15°、70 mm/h，20°、80 mm/h的2場(chǎng)降雨在0～15 min的19個(gè)徑流率、產(chǎn)沙率數(shù)據(jù)作為輸入進(jìn)行預(yù)測(cè)檢驗(yàn)。通過構(gòu)建的偏最小二乘方法分別對(duì)2場(chǎng)降雨在0～15 min的徑流率、產(chǎn)沙率進(jìn)行預(yù)測(cè)，分別得到徑流率、產(chǎn)沙率的預(yù)測(cè)值各19個(gè)，其中1～10為15°、70 mm/h場(chǎng)次的數(shù)據(jù)，11～19為20°、80 mm/h場(chǎng)次的數(shù)據(jù)。圖3為實(shí)測(cè)值與預(yù)測(cè)值的對(duì)比圖。

圖3 實(shí)測(cè)值與預(yù)測(cè)值對(duì)比圖

平均相對(duì)誤差絕對(duì)值作為評(píng)價(jià)精度誤差分析的一個(gè)重要指標(biāo)[15]，采用該指標(biāo)進(jìn)行精度分析。根據(jù)平均相對(duì)誤差絕對(duì)值公式：

平均絕對(duì)值公式中n為數(shù)據(jù)個(gè)數(shù)，y為對(duì)應(yīng)的因變量。將上述2場(chǎng)降雨在0～15 min的19個(gè)徑流率、產(chǎn)沙率及其預(yù)測(cè)數(shù)據(jù)利用上述公式進(jìn)行計(jì)算，得到基于偏小二乘法的平均相對(duì)誤差絕對(duì)值，Y1的平均相對(duì)誤差絕對(duì)值為18.313%，Y2的平均相對(duì)誤差絕對(duì)值為26.698%。同時(shí)得到各點(diǎn)的相對(duì)誤差絕對(duì)值，見表3。

表3 徑流率、產(chǎn)沙率誤差分析表

從圖4中也可以看到，徑流率各點(diǎn)實(shí)測(cè)值與預(yù)測(cè)值的相對(duì)誤差絕對(duì)值基本要小于產(chǎn)沙率各點(diǎn)實(shí)測(cè)值與預(yù)測(cè)值的相對(duì)誤差絕對(duì)值。1～10點(diǎn)的15°、70 mm/h場(chǎng)次無精度分析結(jié)果表明偏最小二乘的預(yù)測(cè)精度較高，可解釋性強(qiáng)，能滿足實(shí)際的需要，特別是相對(duì)于產(chǎn)沙率，徑流率的預(yù)測(cè)精度要高得多，模擬預(yù)測(cè)的效果更好。論是徑流率還是產(chǎn)沙率，誤差值都較11～19點(diǎn)的20°、80 mm/h場(chǎng)次要大，推測(cè)可能是由于在15°附近存在臨界坡度，侵蝕變化規(guī)律相對(duì)復(fù)雜，導(dǎo)致預(yù)測(cè)誤差較大。兩場(chǎng)實(shí)驗(yàn)中前2～3點(diǎn)的誤差均要大于后期誤差值，可能是由于實(shí)驗(yàn)中觀測(cè)者對(duì)產(chǎn)流開始時(shí)間判斷不一致，初始產(chǎn)流階段受人為因素影響較大導(dǎo)致。

圖4 實(shí)測(cè)值和預(yù)測(cè)值相對(duì)誤差圖

總體來看，徑流率的預(yù)測(cè)效果要優(yōu)于產(chǎn)沙率預(yù)測(cè)效果，可能是由于與產(chǎn)沙率相比徑流的影響因素相對(duì)簡(jiǎn)單，隨5個(gè)因變量變化的規(guī)律性更強(qiáng)，更適合進(jìn)行模擬預(yù)測(cè)。

3 結(jié) 論

(1)變量對(duì)數(shù)化處理后，構(gòu)建的偏最小二乘模型模擬精度和預(yù)測(cè)優(yōu)度都較高，均達(dá)到了0.8以上，降雨初期坡面徑流率和產(chǎn)沙率的相關(guān)性較好，相關(guān)系數(shù)達(dá)到了0.844。

(2)時(shí)間X3、坡度X1的VIP值結(jié)果表明，二者都是驅(qū)動(dòng)因變量徑流率Y1和產(chǎn)沙率Y2變化的重要因素，其中時(shí)間X3是最重要的影響因素；雨強(qiáng)對(duì)徑流率和產(chǎn)沙率影響均較大，標(biāo)準(zhǔn)回歸系數(shù)分別為0.408、0.341；初始含水率、剪切力的標(biāo)準(zhǔn)回歸系數(shù)最大僅為0.145，表明二者對(duì)徑流率和產(chǎn)沙率的影響都比較有限，坡度主要對(duì)產(chǎn)沙率起作用，標(biāo)準(zhǔn)回歸系數(shù)達(dá)到0.412，對(duì)徑流率的影響在模型中沒有體現(xiàn)，標(biāo)準(zhǔn)回歸系數(shù)僅0.039；

(3)精度分析結(jié)果表明，偏最小二乘模型在黃土坡面土壤侵蝕中適用性較好，對(duì)徑流率和產(chǎn)沙率模擬和預(yù)測(cè)精度都較高，但與產(chǎn)沙率相比，徑流率的預(yù)測(cè)精度更高，平均相對(duì)誤差絕對(duì)值為18.313%，模擬預(yù)測(cè)的效果更好。

□