基于隨機(jī)森林的農(nóng)耕區(qū)土壤有機(jī)質(zhì)空間分布預(yù)測

楊煜岑，楊聯(lián)安，*，任 麗，李聰莉，朱群娥，王天泰，李新堯

(1.西北大學(xué) 陜西省地表系統(tǒng)與環(huán)境承載力重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710127； 2.西北大學(xué) 城市與環(huán)境學(xué)院，陜西 西安 710127； 3.陜西省周至縣土壤肥料工作站，陜西 西安 710400)

土壤有機(jī)質(zhì)(soil organic matter， SOM)是土壤最具代表性的屬性之一，是評價土壤肥力和土壤質(zhì)量的重要指標(biāo)[1]。土壤有機(jī)質(zhì)的精準(zhǔn)預(yù)測及空間制圖是了解區(qū)域土壤有機(jī)質(zhì)分布狀況、掌握土壤質(zhì)量信息的重要途徑，能夠?yàn)橛行Ч芾硗寥馈⒖茖W(xué)規(guī)劃土地、實(shí)現(xiàn)精準(zhǔn)施肥等提供信息支撐與調(diào)控建議。

在數(shù)字土壤制圖(digital soil mapping， DSM)領(lǐng)域，多種方法和模型都被用來預(yù)測土壤屬性的空間分布[2]。地統(tǒng)計(jì)學(xué)的克里格插值法因其原理簡單、易操作而得到廣泛應(yīng)用。其中：普通克里格(ordinary kriging， OK)插值法對采樣點(diǎn)數(shù)量及密度要求極高，但不考慮環(huán)境變量的影響，限制了其預(yù)測精度[3]；協(xié)同克里格(cokriging，COK)插值法可以利用與目標(biāo)變量具有較好相關(guān)性的輔助變量，獲取較高的插值精度[4]。多元線性回歸模型結(jié)構(gòu)簡單，易于計(jì)算和解釋[5-7]。然而，土壤有機(jī)質(zhì)與輔助變量間并不是簡單的線性關(guān)系，因而傳統(tǒng)預(yù)測難以達(dá)到理想效果。準(zhǔn)確表達(dá)土壤與多種環(huán)境因素間復(fù)雜的非線性關(guān)系已成為土壤養(yǎng)分空間預(yù)測的關(guān)鍵問題。

近年來，神經(jīng)網(wǎng)絡(luò)模型[8]、決策樹[9]等機(jī)器學(xué)習(xí)方法被引入土壤養(yǎng)分預(yù)測領(lǐng)域，雖然提高了預(yù)測精度，但也容易產(chǎn)生過度擬合[10-12]。隨機(jī)森林(random forest， RF)模型是一種自然的非線性建模工具，由決策樹發(fā)展而來，計(jì)算效率高，克服了其他人工智能模型易過度擬合和不穩(wěn)定的缺陷，能夠準(zhǔn)確刻畫變量間的非線性關(guān)系[13-15]，RF算法還能評價指標(biāo)的重要性[16-17]。王茵茵[18]、齊雁冰等[19]運(yùn)用RF模型對陜西省土壤有機(jī)質(zhì)進(jìn)行預(yù)測，并對環(huán)境變量進(jìn)行重要性分析，結(jié)果表明，RF模型準(zhǔn)確可行。鑒于此，以國家測土配方施肥示范縣及耕地地力提升項(xiàng)目示范縣——陜西省周至縣為研究區(qū)域，探索隨機(jī)森林方法在土壤養(yǎng)分空間預(yù)測中的應(yīng)用，以期為土壤養(yǎng)分精準(zhǔn)管理與測土配方施肥提供理論依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

周至縣地處西安市西南部，地處107°39′～108°37′E、33°42′～34°14′N之間，總面積2 947 km2。該縣為暖溫帶大陸性季風(fēng)氣候，年均降水量674.3 mm，年均氣溫15.2 ℃，多年平均無霜期225 d。研究區(qū)位于周至縣農(nóng)耕區(qū)，從南向北地貌類型分屬為黃土臺塬、渭河平原，土壤類型主要有塿土、潮土、黃土和水稻土，主要種植獼猴桃、小麥、玉米。

1.2 數(shù)據(jù)來源及處理

1.2.1 土壤樣品采集與分析

2014年9月，利用GPS技術(shù)，遵循土壤采樣的全面性、均衡性、客觀性、可比性原則，根據(jù)NY/T 1118—2006《測土配方施肥技術(shù)規(guī)范》，分別在小麥-玉米種植區(qū)和獼猴桃果園選取0～20 cm和0～40 cm的耕層采樣，根據(jù)土壤類型及地貌類型等實(shí)際情況調(diào)整采樣點(diǎn)位置。選取每一土樣具有代表性的田塊，采用“S”法均勻隨機(jī)采取8點(diǎn)混勻，用四分法留取1 kg土樣裝袋以備分析，采用重鉻酸鉀加熱法測定土壤有機(jī)質(zhì)[6]。

將測得的土壤養(yǎng)分?jǐn)?shù)據(jù)通過編號匹配地理位置與其他屬性，利用ArcGIS 10.3數(shù)字化后保存到樣點(diǎn)Geodatabase中。為保證研究結(jié)果的精確性和可靠性，采用域值識別法剔除離群值，最終取得研究區(qū)192個有效土壤采樣點(diǎn)數(shù)據(jù)(圖1)。

1.2.2 地形因子

在地理空間數(shù)據(jù)云平臺(http://www.gscloud.cn)獲取30 m分辨率的ASTER GDEM數(shù)據(jù)?；跀?shù)字高程模型(digital elevation model，DEM)，在ArcGIS 10.3軟件中計(jì)算獲取研究區(qū)7個地形因子：坡度(slope)， 坡向(aspect)，地形起伏度(relief)，坡面曲率(slope curvature， SLOC)，剖面曲率(profile curvature， PROC)，徑流強(qiáng)度指數(shù)(stream power index， SPI)，地形濕度指數(shù)(topographic wetness index， TWI)。

1.2.3 氣候因子

從中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn)下載中國地面氣象資料月值數(shù)據(jù)集，獲取西安市地面氣象觀測站數(shù)據(jù)，在ArcGIS 10.3平臺上進(jìn)行線性插值，得到研究區(qū)月平均降水量和月平均氣溫數(shù)據(jù)。

圖1 研究區(qū)采樣點(diǎn)分布Fig.1 Distribution of sampling points in study area

1.3 變量篩選

影響土壤有機(jī)質(zhì)的因素眾多，在RF預(yù)測過程中并非參與預(yù)測的輔助變量越多、結(jié)果越精確，而是與輔助變量對土壤有機(jī)質(zhì)的重要性和影響程度有關(guān)。選取地形、氣候、土壤類型等因子，通過袋外誤差(out of bag error，OOB error)對參與預(yù)測的因子進(jìn)行篩選。RF模型的袋外誤差是對預(yù)測誤差的無偏估計(jì)，是殘差的平方[20]，依據(jù)逐次剔除掉因子后袋外誤差的增減判斷該因子是否保留，袋外誤差增加則刪掉該因子，反之則保留[19]，從而實(shí)現(xiàn)RF模型的因子篩選。通過試驗(yàn)，剔除土壤類型和坡向因子，最后參與RF預(yù)測的輔助變量包括：年均氣溫、年均降水量、DEM、坡度、地形起伏度、SLOC、PROC、SPI、TWI等9個因子。

1.4 研究方法

RF模型是由樹結(jié)構(gòu)的集合組成的分類器，由多棵決策樹構(gòu)成，是用隨機(jī)的方式建立一個森林，每棵決策樹之間沒有關(guān)聯(lián)?；舅枷胧峭ㄟ^bootstrap重采樣方法，從原始訓(xùn)練樣本集中有放回地重復(fù)隨機(jī)抽取多個樣本生成新的訓(xùn)練樣本集合，對每一個bootstrap進(jìn)行決策樹建模，然后采用bagging(bootstrap aggregating)方法組成隨機(jī)森林進(jìn)行預(yù)測。該算法調(diào)解參量少、運(yùn)算效率高，而且能夠處理高維(特征變量多)數(shù)據(jù)，訓(xùn)練速度快，且不會過度擬合。RF算法對特征選取具有較好的穩(wěn)健性，對異常值和噪聲具有較高的容忍度，且可對數(shù)據(jù)進(jìn)行推廣和范化。

1.5 精度檢驗(yàn)

為了驗(yàn)證RF算法的預(yù)測精度，從原始192個土壤樣點(diǎn)中隨機(jī)選取29個采樣點(diǎn)(15%)作為抽樣驗(yàn)證數(shù)據(jù)集。求取驗(yàn)證集土壤有機(jī)質(zhì)的實(shí)測值和預(yù)測值的平均誤差 (ME)、絕對平均誤差 (MAE)、均方根誤差(RMSE)和相關(guān)系數(shù)(r)，作為判斷依據(jù)。

2 結(jié)果與分析

2.1 土壤有機(jī)質(zhì)描述性統(tǒng)計(jì)分析

研究區(qū)土壤有機(jī)質(zhì)值域范圍11.16～30.10 g·kg-1，訓(xùn)練集平均值為19.47 g·kg-1，驗(yàn)證集平均值19.54 g·kg-1。根據(jù)《陜西省第二次土壤普查養(yǎng)分分級標(biāo)準(zhǔn)》，研究區(qū)土壤有機(jī)質(zhì)平均含量處于第四等級(15～20 g·kg-1)，屬于中等水平。從變異程度來看，訓(xùn)練集和驗(yàn)證集的變異系數(shù)均屬于[10%，100%]，表現(xiàn)為中等程度變異性。數(shù)據(jù)正態(tài)性檢驗(yàn)結(jié)果表明，土壤有機(jī)質(zhì)為偏態(tài)分布，通過對數(shù)變換后可轉(zhuǎn)換為近似正態(tài)分布(表1)。

2.2 模型構(gòu)建

2.2.1 隨機(jī)森林預(yù)測

利用R 3.3.2軟件中的random forest軟件包SD， Standard deviation; CV， Coefficient of variation.

對RF模型進(jìn)行擬合。在對RF模型擬合前，需要對模型中的2個參數(shù)進(jìn)行設(shè)定，分別是ntree和mtry。ntree為決策樹的數(shù)量，即用boot-strap重抽樣的次數(shù)；mtry為每次樹模型重建時節(jié)點(diǎn)分裂的次數(shù)[18, 21]。這2個參數(shù)的設(shè)定依據(jù)RF模型運(yùn)算時產(chǎn)生的袋外誤差的大小具體調(diào)整，可通過plot( )函數(shù)查看參數(shù)設(shè)定對誤差的影響。在試驗(yàn)中(圖2)發(fā)現(xiàn)，當(dāng)ntree達(dá)到600之后，誤差趨于穩(wěn)定，經(jīng)過調(diào)整，本研究設(shè)定ntree為800、mtry為2。

從原始樣點(diǎn)中隨機(jī)抽取163個(85%)采樣點(diǎn)作為訓(xùn)練集樣點(diǎn)，利用篩選出的9個影響因子，在R 3.3.2軟件中用訓(xùn)練集樣點(diǎn)對模型進(jìn)行訓(xùn)練。將擬合的RF模型以30 m分辨率預(yù)測農(nóng)耕區(qū)土壤有機(jī)質(zhì)含量空間分布，并將預(yù)測結(jié)果以文本格式保存，在ArcGIS 10.3中轉(zhuǎn)換為柵格格式文件(圖3)，分析研究區(qū)土壤有機(jī)質(zhì)的空間分布特征。

2.2.2 kriging預(yù)測

為了驗(yàn)證RF模型的預(yù)測與制圖效果，基于相同訓(xùn)練集利用OK和COK插值方法對采樣點(diǎn)有機(jī)質(zhì)進(jìn)行空間插值預(yù)測，與RF土壤有機(jī)質(zhì)空間預(yù)測結(jié)果進(jìn)行對比分析。在GS+7.0得出插值時的最優(yōu)半方差模型是球狀模型。進(jìn)行COK插值時，選取2種土壤養(yǎng)分指標(biāo)作為協(xié)變量，根據(jù)相關(guān)性控制協(xié)變量進(jìn)入模型的順序，在SPSS 22.0軟件進(jìn)行Pearson相關(guān)分析，有機(jī)質(zhì)與堿解氮、有效磷、速效鉀的相關(guān)系數(shù)分別是0.324、0.331、0.210，因此依次進(jìn)入的協(xié)變量分別是堿解氮和有效磷。插值結(jié)果如圖3所示。

圖2 決策樹數(shù)量與袋外誤差關(guān)系圖Fig.2 Relation between number of decision tree and OOB error

a，隨機(jī)森林；b，普通克里格；c，協(xié)同克里格。a， Random forest; b， OK; c， COK.圖3 不同模型的插值結(jié)果Fig.3 Spatial distribution of different models

2.2.3 土壤有機(jī)質(zhì)含量空間分布預(yù)測

由圖3可知，RF、OK、COK模型預(yù)測的土壤有機(jī)質(zhì)含量空間分布趨勢基本相似，但RF預(yù)測結(jié)果對局部的刻畫更細(xì)致，克服了克里格插值結(jié)果的不合理突變問題。研究區(qū)中、南部的黑河?xùn)|岸區(qū)，即獼猴桃主產(chǎn)區(qū)土壤有機(jī)質(zhì)含量相對較高，東北部渭河沿岸含量較低，大致以樓觀鎮(zhèn)為中心，由內(nèi)向外呈2個半環(huán)狀分布，在空間上南北差異明顯，整體呈現(xiàn)出由南向北逐漸遞減的趨勢。高值區(qū)域主要分布在樓觀鎮(zhèn)、馬召鎮(zhèn)、廣濟(jì)鎮(zhèn)，以及終南鎮(zhèn)、司竹鎮(zhèn)的部分地區(qū)，可能與當(dāng)?shù)刈匀粭l件、作物類型、施肥習(xí)慣有關(guān)。該區(qū)域位于河流沿岸，土壤含水量高，更易于有機(jī)質(zhì)累積和土壤固碳，且主要種植作物為獼猴桃，與小麥、玉米等大田作物相比需肥量較大，在果園管理中也比較注重增加有機(jī)肥投入，加之果園種草，有利于土壤保墑，故有機(jī)質(zhì)含量偏高。低值區(qū)域主要分布在北部渭河沿岸地區(qū)，密集種植農(nóng)作物，農(nóng)戶施肥行為不合理，偏重氮磷鉀等化肥投入，少施農(nóng)家肥、有機(jī)肥，易造成土壤板結(jié)，僅靠作物根茬、根系分泌物等作物歸還難以維持土壤原有的有機(jī)質(zhì)水平，造成土壤養(yǎng)分比例失衡，有機(jī)質(zhì)含量偏低。

根據(jù)趙業(yè)婷等[22]對周至縣農(nóng)耕區(qū)的現(xiàn)狀分析及實(shí)地調(diào)研，研究區(qū)主要種植物獼猴桃是多年生果樹，枝梢年生長量遠(yuǎn)大于一般果樹，需要從土壤中吸收大量營養(yǎng)，而土壤中可供養(yǎng)分有限，要通過施肥補(bǔ)充樹體所需營養(yǎng)[23]，且我國無公害果樹技術(shù)規(guī)程要求土壤有機(jī)質(zhì)含量達(dá)15 g·kg-1以上，最好能達(dá)到30 g·kg-1[24]，故在今后的施肥管理中，應(yīng)大力增加有機(jī)肥投入量和施用比例，以滿足其養(yǎng)分需求?？偟膩碚f，周至縣北部農(nóng)耕區(qū)的土壤有機(jī)質(zhì)平均含量為19.47 g·kg-1，屬于中等及偏下水平，有機(jī)肥施用不足。研究結(jié)果與雷寶佳等[25]和張彬等[26]的研究結(jié)果相吻合，說明了RF對農(nóng)耕區(qū)土壤有機(jī)質(zhì)預(yù)測的準(zhǔn)確性與適應(yīng)性。

2.3 預(yù)測模型精度驗(yàn)證

通過29個(15%)獨(dú)立驗(yàn)證點(diǎn)對預(yù)測結(jié)果進(jìn)行精度驗(yàn)證。結(jié)果表明，引入輔助變量進(jìn)行土壤有機(jī)質(zhì)預(yù)測的方法(RF和COK)明顯優(yōu)于僅基于樣點(diǎn)值進(jìn)行空間內(nèi)插的預(yù)測方法(OK)，其中又以引入環(huán)境變量的隨機(jī)森林模型(RF)對土壤有機(jī)質(zhì)的預(yù)測精度最高(表2)。RF預(yù)測結(jié)果的ME、MAE、RMSE較OK分別降低了26.6%、27.1%和10.3%，與COK相比，分別降低了7.4%、26.5%和6.9%，預(yù)測精度顯著提高，具有較好的回歸預(yù)測效果。

2.4 變量重要性評價

變量重要性評價是RF算法的一個突出優(yōu)點(diǎn)，基本依據(jù)是當(dāng)某解釋變量加入噪聲后，RF預(yù)測準(zhǔn)確率降低的程度。選擇importance( )函數(shù)提供的節(jié)點(diǎn)純度(IncNodePurity)指數(shù)來評價各解釋變量的重要性，它通過計(jì)算所有樹變量分割的節(jié)點(diǎn)不純性減少值來比較變量的重要性，該值越大表示該變量越重要[27]。計(jì)算可知，影響研究區(qū)土壤有機(jī)質(zhì)的變量按重要性大小排序依次為DEM>年均降水量>年均氣溫>TWI>坡度>地形起伏度>SPI>SLOC>PROC。周至縣北部農(nóng)耕區(qū)對土壤有機(jī)質(zhì)含量影響程度最大的因素是DEM，這與雷寶佳等[25]研究結(jié)果相吻合。

同列數(shù)據(jù)后無相同小寫字母的表示差異顯著(P<0.05)。

Data followed by no same letters within the same column indicated significant difference atP<0.05.

3 結(jié)論與討論

本研究表明，周至縣農(nóng)耕區(qū)土壤有機(jī)質(zhì)含量處于中等偏低水平，中、南部的黑河?xùn)|岸區(qū)土壤有機(jī)質(zhì)含量相對較高，北部渭河沿岸含量較低。高值區(qū)出現(xiàn)在樓觀鎮(zhèn)、馬召鎮(zhèn)、廣濟(jì)鎮(zhèn)，以及終南鎮(zhèn)、司竹鎮(zhèn)的部分地區(qū)；而低值區(qū)主要出現(xiàn)在尚村鎮(zhèn)、竹峪鎮(zhèn)等的局部地區(qū)。根據(jù)實(shí)際情況，農(nóng)戶應(yīng)更加重視耕地培肥和改良，尤其是要增施有機(jī)肥。

DEM和年均降水量是影響該區(qū)域農(nóng)耕區(qū)土壤有機(jī)質(zhì)含量空間分布的主要環(huán)境因素。DEM主要通過不同高程土壤受到的不同外力侵蝕影響土壤有機(jī)質(zhì)的積累，降水量主要通過影響土壤含水率來影響土壤有機(jī)質(zhì)的空間分布。RF方法準(zhǔn)確地揭示了土壤養(yǎng)分與多元環(huán)境因子間復(fù)雜的非線性關(guān)系，與實(shí)際情況吻合最佳，且不需要顧慮一般回歸分析面臨的多元共線性的問題，因此，可以應(yīng)用在農(nóng)業(yè)用地土壤養(yǎng)分回歸預(yù)測上，實(shí)現(xiàn)土壤養(yǎng)分的空間預(yù)測、綜合肥力評價及施肥建議。

采用RF模型對周至縣農(nóng)耕區(qū)土壤有機(jī)質(zhì)含量進(jìn)行空間預(yù)測，分析結(jié)果與前人研究一致且與實(shí)際情況相符合，表明RF模型在土壤養(yǎng)分空間預(yù)測上的應(yīng)用具有一定的科學(xué)性和可靠性。從預(yù)測精度來看，RF模型的ME、MAE、RMSE均小于OK、COK，預(yù)測精度顯著優(yōu)于OK模型和COK模型，說明RF模型可以較好地預(yù)測農(nóng)耕區(qū)土壤有機(jī)質(zhì)含量的空間分布，可為農(nóng)耕區(qū)的土壤有機(jī)質(zhì)預(yù)測提供一種新思路和新途徑。