偏最小二乘模型在黃土坡面土壤侵蝕中的應用研究

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

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

0 引 言

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

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

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

1 材料和方法

1.1 實驗裝置與材料

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

1.2 實驗和測量

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

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

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

1.3 PLS建模

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

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

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

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

2 結果與分析

2.1 參數(shù)選取

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

表1 降雨實驗影響因素表

2.2 偏最小二乘模型

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

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

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

表2 變量之間的相關系數(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

同時得到了原始的回歸方程：

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

即：

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

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

圖2 回歸系數(shù)圖

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

2.3 精度分析

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

圖3 實測值與預測值對比圖

平均相對誤差絕對值作為評價精度誤差分析的一個重要指標[15]，采用該指標進行精度分析。根據(jù)平均相對誤差絕對值公式：

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

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

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

圖4 實測值和預測值相對誤差圖

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

3 結 論

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

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

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

□