應用提升回歸樹研究碳鋼的土壤腐蝕規律

魯慶，穆志純

(北京科技大學 自動化學院，北京，100083)

埋在土壤中的材料因電化學效應、微生物侵蝕等而發生的腐蝕現象稱作材料的土壤腐蝕。通過研究材料的土壤腐蝕規律，可以為合理選擇工程材料、評估設施狀態和采取有效的防護措施等提供指導原則。近年來，國內外研究學者利用多元回歸[1]、時間序列分析[2]、神經網絡[3-4]、灰色系統[5]等方法建立各種腐蝕規律模型。但由于自然環境腐蝕的實驗數據往往存在高維度、小樣本、高噪聲、樣本間相關性等問題，這些方法雖然各具優勢，也存在一些瑕疵，例如多元線性回歸、時間序列等方法需要復雜的數據預處理過程且模型精確度較低；神經網絡等方法可以較好地擬合實驗數據，但其所建立的模型在解釋腐蝕規律影響因素方面卻存在局限性等等。提升回歸樹是集成算法的一種，由Friedman[6]最早提出，其原理是將大量簡單回歸分類樹(CART)在提升過程中進行集成，以提高樹模型的預測能力。該算法在搜索排名、生物研究等多個領域[7-9]得到廣泛應用，在材料腐蝕研究領域應用還很少。本文作者利用提升回歸樹對以碳鋼為例的材料土壤腐蝕規律進行研究，并提出使用ε 不敏感損失函數、動態收縮系數對原算法進行改進，仿真數據和實驗數據驗證表明，改進后的算法更適合小樣本數據的處理。

1 提升回歸樹

提升是近年來最有效的機器學習思想之一[10]，其代表AdaBoost 及衍生算法在分類問題上得到成功應用，提升回歸樹算法將其擴展以解決回歸問題。

1.1 提升回歸樹算法

不失一般性，設有一個數據總體，輸入變量X={x1，x2，x3，…，xp}，共有p 個因子，X 包含N 組觀測數據。輸出變量為Y。總體的樣本{Yi，Xi}1N為已知的實驗觀測數據。

提升算法本質上是一種利用M 個基函數的加法展開式對目標函數進行逼近的方法[10]，可以表示為

其中：βm為展開式系數；b(X，γm)為基函數。

建立模型的過程就是求取一個函數F*(X)，令

其中：L(Y，F(X))為損失函數。

提升回歸樹算法使用分類回歸樹(CART)作為基函數。為便于描述，一棵包含J 個端節點的樹可以表示為

提升回歸樹模型表示為M 棵樹累加的形式：

要使FM(X)逼近F*(X)，對于每一棵加入模型的基樹，要求

即利用當前的損失函數梯度值和輸入變量X 為新的樣本，來訓練下一棵基樹，進而更新提升樹模型：

其中：η 為收縮系數，用于正則化處理[10]。

1.2 模型的解釋性

在很多應用中，模型的可解釋性十分重要，例如本文所研究的碳鋼土壤腐蝕模型，除了要求準確地預測腐蝕速率，還希望通過模型了解哪些變量是造成碳鋼在土壤中腐蝕的主要因素。

其中：T 為具有L 個終端節點的樹；υt為在非終端節點t 上的被選中的屬性測試條件； It2為對應的平方誤差風險改進。對于如式(4)所示的提升回歸樹模型，通過對式(9)代表的每一棵樹求和然后平均，即可得到輸入變量j 與響應變量的相關性度量[6]。

其中：Xc的值取自訓練樣本，如X1c為樣本中Xc子集的第1 條數據的取值，X2c為對應第2 條數據的取值。盡管Xc的取值會影響F(Xs)，但在一般情況下，式(11)可以提供F(X)關于子集Xs的偏依賴性信息。特別對于一維或二維子集，還可以繪制對應的偏依賴圖來直觀了解數據特性。

2 算法改進

2.1 損失函數選擇

Friedman[6]采用平方誤差、絕對值誤差等作為損失函數訓練提升回歸樹，后者對長尾誤差分布、離群點有很強的健壯性，在大多數實際數據集中表現良好。

具體到本文的碳鋼土壤腐蝕數據，如前所述，由于自然環境的實驗數據獲取困難，只有70 個樣本，對于這樣的小樣本數據，除了考慮健壯性，還要考慮模型的過擬合問題。為此，本文嘗試采用ε 不敏感損失函數改進原算法。

ε 不敏感損失函數[12]的定義為

2.2 動態收縮系數

由于提升回歸樹的每次提升過程均會減小訓練損失L(Y, F(X))，為了保持模型的泛化能力，原算法一方面使用獨立檢驗集或交叉驗證等方式控制提升次數M；另一方面，如式(8)所示，在每棵樹添加到當前模型時，均乘以收縮系數η(0＜η＜1)，用于限制單棵基樹的貢獻度。經驗表明，較小的η 有助于改善模型性能。

Friedman 指出通過增加隨機性，可提高模型的性能，原算法的實現方式是對輸入樣本X 進行隨機抽樣[13]。作者在實驗中發現，對于小樣本數據，抽樣后結果并不穩定，有時會引起提升次數M 大幅度增加，模型的預測準確度也較低。考慮其原因是由于小樣本背景下，部分抽樣樣本訓練的基樹性能較差和回歸樹算法本身的不穩定[14]引發了這一問題，也可能存在更深層次的原因，需要進一步探討。

基于上述考慮，作者提出采用動態的收縮系數ηm對原算法進行改進。

如前所述，基樹的解為

定義動態的收縮系數ηm為

其中：η 為人工設置的收縮系數， 0＜η＜1。根據ηm的定義，τm為第m 棵基樹的平方誤差與前m-1 棵基樹平方誤差均值的比，ηm為指數函數。τm越小，這棵樹在模型中將具有較高的權重，但不會超過上限收縮系數η。當τm較大時，ηm變小以懲罰劣化的基樹，當τm等于0 時，ηm=η。當a 取不同的值時，將影響ηm關于τm的變化速度。在本文的研究中，令a 等于e-1，圖1所示為動態收縮系數。

2.3 改進后的算法檢驗

綜上所述，改進后的提升回歸樹算法如下：

圖1 動態收縮系數Fig.1 Dynamic shrinkage coefficient

給定：回歸樹端節點數L，收縮系數上限η，樣本抽樣率fra初始化單端節點樹 ∑=F 0(X)=iNLYi argmin 1(,γ)對于m=1 到M：對i=1 到N，計算梯度r■?L(=■■Y,F(x))i sign)i ■[Y-F(x i i)],otherwise=imYFx)?F(x i■■F=F 0,if|-(|ε＜m-1■■■i i隨機抽取fra*N 個輸入變量，產生基樹訓練樣本Wm以(Wm，rim)為樣本，訓練新的回歸樹，樹的端節點數為L，進而求出各端節點區域Rj對于每個端節點，計算賦予該節點的常量∑∈-+γ=j arg1γ)min L(Y,F(X)γ i m X R j J更新 ∑=F(x)=F m m-1(x)+η m γI(X∈R)j j j1*XFx輸出)()F=m(

改進后的提升回歸樹算法為了驗證改進算法的有效性，使用UCI 數據集中Friedman#1 仿真數據對算法進行檢驗。Friedman#1 數據 的 原 函 數 為Y=10sin(πX1X2)+20(X3-0.5)2+10X4+5X5+ε，各變量均為N(0，SD)分布。令SD=1，生成樣本數量為50，100 和500 的3 組數據，分別隨機抽取90%數據作為訓練集，另外10%數據作為驗證集。使用原算法和改進算法分別對訓練集建立模型，并使用驗證集檢驗模型。在模型建立過程中，利用10折交叉驗證以避免過擬合現象，結果如表1 所示。

中藥產業具有衛生資源、經濟資源、科技資源、文化資源和生態資源5大資源優勢，同時具有傳統和現代產業的特點。從傳統產業看，中藥產業是中醫藥體系中不可分割的一部分，其產業特征與化學制藥具有明顯的不同，中藥產業鏈長，依賴自然資源、農業生產和氣候條件，且地域性強[7]。從現代產業看，中藥產業是以中藥農業為基礎、中藥工業為核心、中藥商業為紐帶、中藥知識業為動力的完整產業體系，目前各環節日益規模化、標準化，逐步脫離自然經濟的特征。

從表1 可以看出：改進后的算法對于小樣本數據的處理能力較原算法有所增強，特別對比以平方誤差作為損失函數的模型，預測能力有了顯著提高。無論對于提升回歸樹的原算法還是改進算法，建立在較大樣本數基礎上的模型性能更好。

表1 改進算法效果檢驗Table 1 Results verification of improved algorithm

3 結果與討論

3.1 研究樣本和資料庫

以國家材料腐蝕站網中碳鋼在土壤中的腐蝕數據進行分析和建模，數據示例如表2 所示，包含了實驗站點、埋樣時間、土壤pH、有機質含量、全氮含量等理化參數及實驗材料的腐蝕率，其中type 變量為名義型數據，指示了埋件的種類，其余為連續型數據。本研究的目標是訓練模型來估計碳鋼的腐蝕率，并分析腐蝕率的影響因素。剔除目標數據缺失的記錄，可用于建立模型的共有70 條數據。部分作為輸入變量的數據包含缺失值。

3.2 建立模型

3.2.1 參數選擇

如2.3 節算法所示，建立提升回歸樹模型需要確定基樹端節點數L、收縮系數上限η 和樣本抽樣率fra3 個參數。以下將簡述參數的選擇過程。

表2 碳鋼的有關參數Table 2 Relative parameters of carbon steel

(2) 收縮系數上限η 和提升迭代次數M。為確定收縮系數上限η，以50 為步長，增加提升回歸樹模型中樹的數量，對不同η 分別建立模型。隨機抽取的10%數據對模型進行驗證，預測偏差均值如圖2 所示。從圖2 可以看出：當η 為1 時，隨著樹的增加，模型很快出現了過擬合情況，降低收縮系數上限η 值則會明顯提升模型的性能，但這將造成提升迭代次數M 的增加。通過實驗，并綜合考慮計算機的計算性能，確定η=0.005。為避免模型過擬合，采取10 折交叉驗證的方法在建立模型過程中確定最優的提升迭代次數M。

圖2 不同η 的模型預測偏差Fig.2 Predictive deviance of models with different η

(3) 樣本抽樣系數fra。在2.2 節中提到，為提高模型性能，每棵基樹的訓練不是使用全部的輸入樣本，而是對輸入樣本進行隨機抽樣，通常抽樣系數選擇0.5＜fra＜1。對于小樣本數據，此系數的選擇對模型性能影響較大，令fra={0.55,0.65,0.75,0.85,0.95}，對于每個取值分別進行5 次建模，對其10 折交叉驗證偏差進行平均，結果如表3 所示。經過實驗，最終確定fra=0.65。

表3 不同抽樣系數的交叉驗證偏差均值Table 3 Cross-validation deviation mean of different sampling coefficients

3.2.2 模型分析

綜上所述，確定了建立模型所需的參數{η，fra，L}={0.005，0.65，5}。由于觀測數據很少，為提高模型的性能，使用全部數據對模型進行訓練，使用10折交叉驗證保證模型的泛化能力。模型訓練結果的10折交叉驗證偏差最小值為0.580，提升回歸樹模型累加的基樹數量為15 950 棵，訓練時間為2.75 min。

使用1.2 節中的算法對模型進行分析，各因素的重要性如表4 所示，可以得出埋件的埋存時間、土壤有機物含量和土壤全氮含量是碳鋼在土壤環境中的腐蝕主要影響因素，另外SO42-，H2O，K+和土壤電阻R等對碳鋼土壤腐蝕也有著較大的影響。

表4 腐蝕影響因素的重要性Table 4 Importance of corrosion factors

為考察各影響因素與腐蝕率的關系，可以計算一維偏依賴函數，以重要性系數排前六的因素為例，計算偏依賴函數關系曲線如圖3 所示。

圖3 腐蝕因素與腐蝕率的一維偏依賴函數的關系Fig.3 One-dimensional partial dependence function of corrosion factors about corrosion rate

從圖3 可以看出：各因子與腐蝕率的依賴關系，如隨著試件埋存時間的增加，腐蝕率逐漸下降；在有機質含量和氮含量較高的地區，腐蝕率較低；而在K含量較高的土壤中，腐蝕率會較高。將圖3 與文獻[16]關于土壤腐蝕因素分析部分的內容對比，發現了一些對應內容，如文獻[16]中指出“含水率在10%甚至5%以下，腐蝕速率低；含水率在10%～25%時，腐蝕速率最高；當達35%左右飽和含水率時，腐蝕速率降至最低”，從圖 3(e)中可以看到相似的變化趨勢：在15%～25%時，確實存在1 個尖峰。文獻[16]還提到“硫酸根離子和腐蝕速率之間沒有明確的相關關系”，也可以在圖3(d)中得出相似的結論。

通過計算高維偏相關函數，可以進一步探索各輸入變量對于腐蝕率的交互作用，在文獻[16]中提到“鷹潭、廣州、深圳等試驗站點，屬于酸性土壤，電阻率常常高至幾百至幾千歐姆，但腐蝕速率反而比一般的中堿性土壤高得多，對于這種強腐蝕性，目前還無法解釋”。為考察這一問題，計算土壤電阻R 與pH 的偏依賴相關函數，結果如圖4 所示。

由圖4 可見：碳鋼在酸性土壤的腐蝕率高于堿性土壤的腐蝕率，而對于固定的pH，具有高電阻土壤腐蝕率大于低電阻土壤腐蝕率的趨勢；當電阻較小時，對于任意酸堿度土壤，電阻越小，腐蝕率越高，但當電阻大于30 ? 之后，pH 對腐蝕率的影響變得顯著，高電阻土壤腐蝕率變為高于低電阻土壤的腐蝕率，但變化較平緩。通過計算各輸入變量組的偏依賴函數，就可以對變量交互作用進行探索。

圖4 土壤pH、電阻率因素與腐蝕率二維偏依賴函數Fig.4 Two-dimensional partial dependence function of soil pH and resistivity about corrosion rate

根據模型分析和與土壤腐蝕專業文獻的對比，可以認為提升回歸樹模型可以較好地描述碳鋼的土壤腐蝕規律，盡管不能揭示腐蝕發生的深層原因，但可以為腐蝕規律的研究提供有益的線索。

3.2.3 算法性能檢驗

由于模型是使用全部樣本建立的，為了檢驗改進算法的擬合和泛化性能，作者在原數據集中隨機抽取10 個樣本作為檢驗集，使用其余60 個樣本建立模型，并與原算法、SVM 算法、神經網絡算法、分類回歸樹(CART)算法和線性回歸算法進行比對。

表5 模型性能對比Table 5 Comparison of models performance

其中改進算法的參數選擇為{η，fra，L}={0.005，065，5}，原算法的參數根據模型訓練情況，選擇為{η，fra，L}={0.005，0.75，5}。各自的提升次數M 通過10 折交叉檢驗在建模過程中確定。SVR 算法使用R 軟件包E1071 計算，令ε 等于0.01，采用10折交叉檢驗控制迭代次數。神經網絡算法使用R 軟件包NNET 計算，設置為3 層網絡，限制迭代次數為500次。分類回歸樹(CART)算法使用R 軟件包rpart 計算。線性回歸采用逐步回歸算法。對于作為對比的SVR 算法、神經網絡算法、分類回歸樹算法和逐步回歸算法，沒有進行進一步優化，基本采用默認值。

關于數據中的缺失值，由于數據集樣本較小，數據又采集自不同地區的臺站，缺乏必要的信息以完成缺失數據填補，而且提升回歸樹算法數據對缺失值有很好的魯棒性，因此，在原算法、改進算法和分類回歸樹算法中保持了數據缺失狀態。對于其他算法，如刪除含缺失數據的樣本，需要刪去約11%的數據，因此，采用了同地區均值填充的方式進行缺失數據處理。實驗結果如表5 所示。

圖5 不同算法的預測性能對比Fig.5 Comparison of different algorithms predictable performance

各算法建立的模型在測試集上的結果如圖5 所示，其中Actual Value 為實際值，SVR，Regression，Nnet，BRT，CART 和I-BRT 分別為支撐向量回歸、逐步回歸、神經網絡、提升回歸樹算法、分類回歸樹和本文的改進算法的模型預測值。從表5 和圖5 可以看出：改進的提升回歸樹算法可以對碳鋼土壤腐蝕率進行準確預測，性能好于原算法。雖然SVR 和神經網絡算法沒有經過系統優化，至少可以認為提升回歸樹算法的性能可以與前兩者相當。

4 結論

(1) 提出一種基于提升回歸樹算法的新方法，針對實驗數據小樣本情況下的參數選取問題，采用ε 不敏感損失函數、動態收縮系數對原算法進行改進。實驗證明，改進算法適合小樣本數據的處理。

(2) 建立碳鋼在土壤中的腐蝕率模型。模型較好地描述和預測了碳鋼在土壤中的腐蝕規律，并為探索腐蝕影響因素及其因素間交互作用提供了線索。

(3) 提升回歸算法需要人工確定樹端節點數L、收縮系數上限η 和樣本抽樣率fra等參數，由于各參數相互影響，通過多次實驗分別確定參數的方式并不精確(盡管模型結果對一定區間內的參數是不敏感的)，如何優化模型建立的過程將有待下一步研究。

[1] Yahaya N, Noor NM, Othman RS, et al. New technique for studying soil corrosion of underground pipeline[J]. Journal of Applied Sciences, 2011, 11(9): 1510-1521.

[2] 楊瑞成, 王彬, 張天云. 基于MATLAB 的混沌時間序列算法對材料腐蝕行為的預測[J]. 蘭州理工大學學報, 2009, 35(5):5-8.YANG Ruicheng, WANG Bin, ZHANG Tianyun. Prediction of the corrosion behavior of chaotic time series algorithm based on MATLAB[J]. Journal of Lanzhou University of Technology,2009, 35(5): 5-8.

[3] 賴延清, 陳湘濤, 秦慶偉, 等. NiFe2O4基金屬陶瓷耐腐蝕因素分析及腐蝕率預測[J]. 中南大學學報(自然科學版), 2004,35(6): 896-901.LAI Yanqing, CHEN Xiangtao, QIN Qingwei, et al. Corrosion analysis and corrosion rates prediction of NiFe2O4cermet inert anodes[J]. Journal of Central South University (Science and Technology), 2004, 35(6): 896-901.

[4] Sadowski L. Non-destructive investigation of corrosion current density in steel reinforced concrete by artificial neural networks[J]. Archives of Civil and Mechanical Engineering,2013, 13(1): 104-111.

[5] LIANG Ping, DU Cuiwei, LI Xiaogang. Grey relational space analysis of effect of environmental factors on corrosion resistance of X70 pipeline steel in ying tan soil simulated solution[J]. Corrosion & Protection, 2009(4): 23-27.

[6] Friedman J H. Greedy function approximation: A gradient boosting machine[J]. Annals of Statistics, 2001, 29(5):1189-1232.

[7] Buhrmann P, Hothorn T. Boosting algorithms: Regularization,prediction and model fitting[J]. Statistical Science, 2007, 22(4):477-505.

[8] Mohan A, CHEN Zheng, Weinberger K. Web-search ranking with initialized gradient boosted regression trees[J]. Journal of Machine Learning Research, Workshop and Conference Proceedings, 2011(14): 77-89.

[9] Kneib T, Hothorn T, Tutz G. Variable selection and model choice in geoadditive regression models[J]. Biometrics, 2009, 65(2):626-634.

[10] Hastie T, Tibshirani R, Friedman J H. 統計學習基礎—數據挖掘、推理與預測[M]. 范明, 等譯. 北京: 電子工業出版社,2004: 337-384.Hastie T, Tibshirani R, Friedman J H. The elements of statistical learning: data mining, inference, and prediction[M]. FAN Ming,et al, trans. Beijing: Electronic Industry Press, 2004: 337-384.

[11] Breiman L, Ihaka R. Nonlinear discriminant analysis via scaling and ACE[R]. Berkeley: University of California, 1984: 17-38.

[12] Dekel O, Shalev-Shwartz S, Singer Y. Smooth epsiloninsensitive regression by loss summarization[J]. Journal of Machine Learning Research, 2006, 6(1): 711-741.

[13] Friedman J H. Stochastic gradient boosting[J]. Computational Statistics & Data Analysis, 2002, 38(4): 367-378.

[14] 梁茵. 分類回歸樹算法的探討[J]. 廣東技術師范學院學報,2008(6): 29-32.LIANG Yin. A discussion of classification and regression tree[J].Journal of Guang dong Polytechnic Normal University, 2008(6):29-32.

[15] Elith J, Leathwick J R, Hastie T. A working guide to boosted regression trees[J]. Journal of Animal Ecology, 2008, 77(4):802-813.

[16] 曹楚南. 中國材料的自然環境腐蝕[M]. 北京: 化學工業出版社, 2005: 375-383.CAO Chunan. The environment corrosion of materials in China[M]. Beijing: Chemical Industry Press, 2005: 375-383.