基于XGBoost模型的安徽省土壤pH空間建模

邱士其，趙明松*，蘆園園，盧宏亮，徐少杰，陳宣強

(1.安徽理工大學空間信息與測繪工程學院，淮南 232001；2.礦山采動災害空天地協(xié)同監(jiān)測與預警安徽省教育廳重點實驗室，淮南 232001；3.礦區(qū)環(huán)境與災害協(xié)同監(jiān)測煤炭行業(yè)工程研究中心，淮南 232001；4.國家環(huán)境保護土壤管理與污染控制重點實驗室，南京 210042；5.生態(tài)環(huán)境部南京環(huán)境科學研究所，南京 210042)

土壤pH不僅對土壤中微生物活性、土壤養(yǎng)分有重要影響，而且影響土壤理化性質(zhì)以及植被對土壤養(yǎng)分的吸收、利用[1]，是土壤中影響范圍極為廣泛的化學指標[2]。對于土壤pH的研究，有利于及時了解土壤狀況，為合理施肥、改良土壤、加強土壤環(huán)境管理起到重要作用。

隨著地信和遙感技術(shù)的發(fā)展，逐漸形成了以定量土壤-景觀模型為理論基礎(chǔ)、以空間分析和數(shù)學方法為技術(shù)手段的土壤調(diào)查與制圖的現(xiàn)代化技術(shù)體系-數(shù)字土壤制圖(digital soil mapping,DSM)[3-4]。數(shù)字土壤制圖起初主要運用判別分析、線性回歸等方法來建立土壤與環(huán)境之間的關(guān)系并繪制成圖。隨著研究人員對數(shù)字土壤制圖技術(shù)的不斷研究與完善，目前數(shù)字土壤制圖方法包括統(tǒng)計學方法、基于專家知識的方法和機器學習等方法[5-6]。統(tǒng)計學方法是根據(jù)土壤與環(huán)境變量之間的統(tǒng)計關(guān)系，推測土壤屬性空間分布的方法，常用方法包括：線性模型、判別分析。克里格[7]和地理加權(quán)回歸[8]是空間統(tǒng)計中較為常用的方法。基于專家知識的方法是指從土壤專家處獲得關(guān)于土壤與環(huán)境變量之間關(guān)系的知識，將獲取的知識與語義模型結(jié)合，并借助地信技術(shù)來完成土壤信息制圖的一種方法。Qin等[9]基于地理信息系統(tǒng)(geographic information system，GIS)、模糊邏輯和專家知識，建立了土壤一環(huán)境推理模型, 并考慮了感興趣區(qū)位置和樣本位置之間的空間距離信息,使得模型精度進一步提高。

基于機器學習的數(shù)字土壤制圖方法在處理多維、非線性海量數(shù)據(jù)、提高模型泛化能力、預測速度、精度等方面具有良好表現(xiàn)，目前已被應用到諸多領(lǐng)域[10-11]。盧宏亮等[12]研究發(fā)現(xiàn)，將Boruta算法與SVR模型結(jié)合可以提高土壤pH的預測制圖精度，且模型的泛化能力較強。王世航等[13]研究表明，使用特征挖掘方法可有效提高廣義提升回歸和隨機森林模型預測精度，且可以起到降維的作用。Liu等[14]基于先進的機器學習方法與高效的并行運算環(huán)境，大幅提高了現(xiàn)有制圖的準確性和精細度，并提供了空間預測的不確定性信息，更好地表征了中國土壤屬性的空間變異特征。

在數(shù)字土壤制圖方法中，常見的機器學習模型包括隨機森林[15-16](random forest，RF)、支持向量機(support vector machine，SVM)、神經(jīng)網(wǎng)絡(neural networks，NNs)等。

極端梯度提升模型[17](extreme gradient boosting,XGBoost)采用Boosting思想并使用多種方法防止模型過擬合。目前，XGBoost算法多應用在醫(yī)學[18-19]、金融[20-21]、工程等領(lǐng)域，在數(shù)字土壤制圖領(lǐng)域的應用包括：結(jié)合SFLA-XGBoost模型與高光譜的土壤有機質(zhì)含量預測[22]、基于極限梯度提升和長短期記憶網(wǎng)絡相融合的土壤溫度預測[23]、土壤含水量估算等，但少有研究者將其用于土壤pH的預測和制圖中。因此，探討XGBoost模型在土壤pH預測及制圖中的可行性以及模型所能達到的精度具有實際的研究意義。

鑒于此，以安徽省為研究區(qū)域，采用數(shù)字高程模型(digital elevation model,DEM)和野外實測土壤樣本數(shù)據(jù)，利用3S技術(shù)提取17個環(huán)境變量，以土壤pH為研究目標。建立XGBoost模型并對土壤pH進行預測及制圖，與隨機森林模型結(jié)果進行對比分析。探討兩種模型下安徽省表層土壤pH的空間分布規(guī)律及差異，定量分析土壤pH的影響因素，以可視化的方式給出兩模型的不穩(wěn)定性表達，以期為XGBoost模型在數(shù)字土壤制圖的應用提供參考并為安徽省加強土壤環(huán)境管理提供數(shù)據(jù)基礎(chǔ)及安徽省土壤pH預測、研究提供理論基礎(chǔ)。

1 材料與方法

1.1 研究區(qū)域

安徽省(114°54′E～119°37′E，29°41′N～34°38′N)，居中靠東，沿江通海，是最具活力的長江三角洲組成部分，總面積達14.01×104km2。地形地貌呈現(xiàn)多樣性，境內(nèi)長江和淮河自西向東橫貫全境，將全省劃分成3個自然區(qū)域。境內(nèi)主要山脈包括大別山、黃山、九華山等，最高山峰是黃山蓮花峰，海拔 1 864 m。安徽省地處暖溫帶與亞熱帶過渡地區(qū)，四季分明、雨量充沛、氣候宜人。安徽省土壤的過渡特征明顯，以林為主的土壤主要有棕壤、草甸土、黃棕壤、黃壤、石灰(巖)土、紫色土、麻骨土、石質(zhì)土等；農(nóng)林兼用的土壤有紅壤、黃褐土等；以農(nóng)為主的土壤有砂姜黑土、潮土、水稻土等。

1.2 數(shù)據(jù)處理

研究所使用的土壤數(shù)據(jù)來源于《中國土系志·安徽卷》[24]，所獲數(shù)據(jù)包含采樣點分布、環(huán)境條件和土壤理化性質(zhì)等。采樣時間為2010年，包含140個樣點的表層土壤pH屬性數(shù)據(jù)。地形數(shù)據(jù)為空間分辨率為90 m 的SRTM數(shù)字高程模型，來源于地理空間數(shù)據(jù)云。坡向(aspect)、坡度(slope)、高程(elevation)、平面曲率(plan curvature)和剖面曲率(profile curvature)通過ArcGis10.6獲得；徑流強度指數(shù)(SPI)、匯聚指數(shù)(CI)、多尺度谷底平坦度(MRVBF)、多尺度脊頂平坦度(MRRTF)、地形濕度指數(shù)(TWI)及地形位置指數(shù)(TPI)在SAGAGIS 6.3.0中得到。歸一化植被指數(shù)(NDVI)和增強植被指數(shù)(EVI)來源于MODIS陸地產(chǎn)品，空間分辨率為250 m。氣候數(shù)據(jù)使用來自中國農(nóng)業(yè)科學院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所中國生態(tài)環(huán)境背景層面建造項目完成的柵格數(shù)據(jù)(1 km分辨率)，為1980—1999年的逐月平均值計算生成年均溫(MAT)和年均降水量(MAP)。此外，將x、y坐標作為環(huán)境變量引入建模。

1.3 研究方法

1.3.1 XGBoost模型

XGBoost是由GBDT(gradient boosting decision tree)模型發(fā)展而來的，在傳統(tǒng)的Boosting基礎(chǔ)上，利用中央處理器(central processing unit,CPU)多線程，引入正則化項，加入剪枝，使用列抽樣來控制模型的復雜度。如何在每一步生成合理的樹是 Boosting分類器的核心。

對于一顆含有t個基礎(chǔ)模型的集成樹來說，關(guān)鍵點就是第t個基礎(chǔ)模型的選擇。對于該問題，需要尋找一個能夠使目標函數(shù)盡可能最大化降低的ft，故構(gòu)造的目標函數(shù)可表示為

(1)

對于一般的損失函數(shù)來說，可以使用泰勒展開對損失函數(shù)做近似估計。移除常數(shù)項(真實值與上一輪的預測值之差)，目標函數(shù)只依賴于每個數(shù)據(jù)點在誤差函數(shù)上的一階導數(shù)和二階導數(shù)。

(2)

式(2)中：gi為數(shù)據(jù)點在損失函數(shù)上的一階導數(shù)；hi為數(shù)據(jù)點在損失函數(shù)上的二階導數(shù)。

XGBoost模型參數(shù)有3種類型，即：通用參數(shù)、輔助參數(shù)和任務參數(shù)[25]。通用參數(shù)包括：booster、silent、nthread等。輔助參數(shù)有：eta、gamma、max_depth、subsample等。任務參數(shù)包含：objective、eval_metric、base_score等。對nrounds、max_depth和eta 3個參數(shù)在R語言中利用caret包進行格網(wǎng)搜索及重復交叉驗證并以均方根誤差值為評價指標，選擇均方根誤差值最小時的參數(shù)組合作為最優(yōu)參數(shù)，利用XGBoost包進行建模并進行預測制圖。

1.3.2 隨機森林模型

與XGBoost模型不同，隨機森林模型在Bagging的基礎(chǔ)上，進一步在訓練過程中引入隨機屬性選擇，且模型對于高維數(shù)據(jù)集的處理能力很好，泛化性能優(yōu)越。

隨機森林模型是n棵決策樹{T1(X),T2(X),…，Tn(X)}的集合，其中X={x1，x2，…，xp}是與目標變量相關(guān)的特征的p維向量，結(jié)果為n棵樹的輸出值Y={Y1=T1(X),Y2=T2(X),…，Yn=Tn(X)}，其中Yn為第n棵樹的預測值。對于分類問題，Y由最多投票數(shù)獲得；對于回歸問題，Y是多棵樹預測的平均值獲得。

隨機森林模型中有兩個重要參數(shù)：分裂次數(shù)(mtry)和決策樹數(shù)量(ntree)。在R語言中利用e1071包采用格網(wǎng)搜索和十折交叉驗證進行參數(shù)調(diào)優(yōu)，并以均方誤差為精度評價指標，并用random Forest包進行建模。

1.3.3 精度評價

模型精度評價選用均方根誤差(root mean squared error,RMSE)、平均絕對誤差(mean absolute deviation,MAE)以及決定系數(shù)(R2)3個指標,其計算公式分別為

(3)

(4)

(5)

模型不確定性是數(shù)字土壤制圖中的重要內(nèi)容，不確定性可表示為90%置信區(qū)間的預測下限和上限。隨機森林模型的不確定性可通過R語言中的quantregForest包對最優(yōu)模型中各個樹的結(jié)果進行統(tǒng)計獲得，XGBoost模型的不確定性則是通過對100次最優(yōu)模型預測結(jié)果統(tǒng)計所得。

2 結(jié)果與分析

2.1 土壤屬性特征統(tǒng)計

如表1所示，安徽省土壤pH分布范圍在4.58～8.67，平均水平為6.37，為中性土壤。從總體上看，總樣本、訓練集和驗證集數(shù)據(jù)分布極為相似，標準差在1.16附近，變異系數(shù)也都分布在18%～19%，屬中等變異。在安徽省土壤樣本中，強酸性樣本36個(26%)，酸性樣本45個(32%)，中性樣本33個(24%)，微堿性樣本20個(14%)，強堿性樣本6個(4%)。就整體而言，安徽省土壤樣本中酸性占比大(58%)，堿性樣本則占18%。安徽省土壤樣點位置分布如圖1所示。

表1 土壤pH特征基本統(tǒng)計

圖1 土壤pH采樣點空間分布

2.2 環(huán)境變量與土壤pH相關(guān)性分析

如表2所示，土壤pH與X3(r=0.63)和X4(r=0.61)呈顯著正相關(guān)，與X1、X2、X6、X10、X13、X14顯著負相關(guān)；與X8在0.05水平上顯著負相關(guān)。此外，X14與土壤pH的相關(guān)性最高，說明隨著年均降雨量的增加，pH呈顯著下降的趨勢。除X1、X3、X4和X14與土壤pH均具有較高的相關(guān)性外，其余環(huán)境變量與土壤pH的相關(guān)系數(shù)并不高。此外，X1與X2、X3、X4、X10、X13、X14、X15顯著相關(guān)，與X6、X8相關(guān)(P<0.05)；X2與X3、X4、X10、X13、X14、X15顯著相關(guān)，其余變量之間也存在類似的情況。X13與X15之間存在著最高的相關(guān)性(r=-0.83)。

表2 土壤pH與環(huán)境變量的相關(guān)性分析

2.3 土壤pH建模及參數(shù)優(yōu)化

如圖2所示，在學習率eta=0.005、0.01、0.02時，無論樹深max_depth為何值，隨著nrounds的值不斷變大時，RMSE逐漸減少并在到達一定的值時趨于穩(wěn)定。eta=0.03時，與其他3種情況稍有不同，在max_depth=2的情況中，迭代次數(shù)nrounds在200后值不斷增加時，RMSE也在不斷變大，但在nrounds=1 000后，值同樣趨于穩(wěn)定。可以看出，XGBoost模型對于eta是較為敏感。隨著eta不斷地變化，RMSE范圍越來越小：0.5～4.0、0.6～2.6、0.75～1.20、0.75～0.85，且RMSE趨于穩(wěn)定時，nrounds的值也在不斷變小：1 000、500、300。max_depth 對于RMSE影響較小，在eta=0.005、0.01時，無明顯的差異；在 eta=0.02、0.03時，才能看出明顯差異。最終的參數(shù)組合為：nrounds=860，eta=0.01，max_depth=2。

圖2 XGBoost模型參數(shù)調(diào)優(yōu)結(jié)果

2.4 環(huán)境因素對土壤pH建模影響

在XGBoost模型中，可用gain來衡量特征的重要性。gain 增益意味著相應的特征對通過模型中的每個樹采取每個特征的貢獻而計算出的模型的相對貢獻。與其他特征相比，較高值意味著它對于生成預測更為重要。

在隨機森林模型中，可用IncNodePurity來判斷特征的重要性。該值通過殘差平方和來度量，代表了每個變量對分類樹每個節(jié)點上觀測值的異質(zhì)性的影響，從而比較變量的重要性。該值越大表示該變量的重要性越大。

如圖3所示，在兩個模型排名前5的環(huán)境變量分別為：緯度(Y)、MRVBF、MAP、MAT、EVI和MAP、Y、MRVBF、EVI、MRRTF。環(huán)境變量的選擇稍有不同，XGBoost模型中包含MAT，而RF模型中則有 MRRTF，其余變量對于土壤pH的影響較低。其中，MAP 在兩個模型中的排名均在前列，說明年均降雨量對于土壤pH有重要影響。在多雨環(huán)境中，土壤及其母質(zhì)浸出強烈導致土壤酸化，而安徽省南部降雨多于北方，這會導致南北土壤酸堿度的差異[26]。而MRVBF和MRRTF可以識別谷底和山坡兩種地形，這兩種地形存在著明顯不同的土壤含水量。一般來說，過量的水分導致土壤酸化，過少的水分導致土壤堿化。

X為經(jīng)度

3 討論

3.1 模型精度及空間制圖

如表3所示，XGBoost模型的精度為：訓練集(RMSE=0.32、MAE=0.25、R2=0.93)、驗證集(RMSE=0.65、MAE=0.50、R2=0.73)；隨機森林模型的精度結(jié)果為：訓練集(RMSE=0.75、MAE=0.60、R2=0.57)，驗證集(RMSE=0.68、MAE=0.53、R2=0.71)。可以看出，XGBoost模型在訓練集的精度是最高的，明顯優(yōu)于驗證集的精度和隨機森林模型的精度，XGBoost模型驗證集和訓練集精度均優(yōu)于隨機森林模型訓練集的精度，較優(yōu)于隨機森林模型驗證集的精度，且隨機森林模型訓練集精度低于驗證集精度。綜上，XGBoost模型對于土壤pH預測結(jié)果精度優(yōu)于隨機森林模型預測結(jié)果。

表3 模型精度

如圖4所示，兩種模型的預測結(jié)果趨勢大體相同，安徽省土壤pH呈“南酸北堿”趨勢。以長江為分界線，長江以南多為酸性或強酸性土壤，長江以北多為堿性土壤或中性土壤。兩模型土壤pH預測的范圍不同：XGBoost模型(4.59～8.96)、隨機森林模型(5.52～8.16)，差別較為明顯。且在部分地區(qū)兩者略有差異。

圖4 模型預測結(jié)果圖

3.2 模型不確定性

如圖5所示，5%的預測下限和95%的預測上限與最優(yōu)預測結(jié)果具有相似的空間模式(安徽省北部土壤pH要高于南部)。隨機森林模型不確定性較大的區(qū)域集中在合肥、阜陽、淮南等長江以北區(qū)域，最大的差值為3.99，最小為1.64；XGBoost模型不確定性較大的區(qū)域則相對集中在靠近長江的地方，包括六安、黃山、宣城等地，最大差值為1.73、最小為0.38。從結(jié)果上看，XGBoost模型的穩(wěn)定性要高于隨機森林模型。

圖5 模型不確定性圖

4 結(jié)論

基于XGBoost、隨機森林模型建立了安徽省表層土壤pH空間預測模型，并對比分析了模型精度和不確定性。得出如下主要結(jié)論。

(1)XGBoost模型的調(diào)參結(jié)果說明eta、max_depth、nrounds對于模型的精度均具有一定的影響。其中，eta的變化對于模型的影響較大。

(2)MAP、Y、MRVBF、EVI對土壤pH建模有重要的影響。

(3)XGBoost模型在訓練集和驗證集預測精度均高于隨機森林模型，訓練集精度明顯高于隨機森林模型訓練集精度，驗證集精度較優(yōu)于隨機森林模型驗證集精度。模型5%的預測下限和95%的預測上限與最優(yōu)預測結(jié)果具有相似的空間模式，從結(jié)果上看，XGBoost模型的穩(wěn)定性要高于隨機森林模型。兩模型預測結(jié)果都表明了安徽省土壤pH呈“南酸北堿”的趨勢，但兩者在部分地區(qū)結(jié)果仍有區(qū)別。