朱福斌，丁世偉，甘曉玉，黃 海，吳錦松，馬友華*

（1.安徽農(nóng)業(yè)大學資源與環(huán)境學院，安徽 合肥 230036；2.安慶市種植業(yè)管理局，安徽 安慶 246000）

優(yōu)越的土壤肥力是保證糧食作物產(chǎn)量的基礎，在現(xiàn)代農(nóng)業(yè)生產(chǎn)中有著決定性作用［1-2］。而土壤速效鉀是土壤肥力中重要的影響因素，是作物生長發(fā)育必不可少的養(yǎng)分元素之一，對作物產(chǎn)量有著重要影響［3-4］。耕地土壤速效鉀空間分布主要受空間位置、地形部位、成土母質、氣候環(huán)境、人為施肥等多種因素影響［5-7］。相對而言，安慶市耕地土壤速效鉀含量整體偏低［8］，因此準確預測其含量及分布情況對指導農(nóng)業(yè)施肥、促進作物增產(chǎn)、提升耕地質量等級尤為重要。

傳統(tǒng)數(shù)字土壤制圖的空間預測方法主要包括地統(tǒng)計插值法和確定性插值法［9］。地統(tǒng)計插值法的代表性方法是普通克里金插值法［10］，確定性插值法的代表性方法是反距離權重插值法［11］。近年來，機器學習方法，尤其是隨機森林模型逐步成為數(shù)字土壤制圖方向的國內外研究熱點。Kalambukattu等［12］在遙感和地形數(shù)據(jù)的基礎上，利用人工神經(jīng)網(wǎng)絡模型繪制了喜馬拉雅流域土壤有機碳和氮元素空間分布圖；高鳳杰等［13］運用貝葉斯最大熵模型預測了黑土區(qū)小流域土壤有機質空間分布，通過與協(xié)同克呂格法比較，結果表明貝葉斯最大熵模型顯著優(yōu)于協(xié)同克呂格法。然而人工神經(jīng)網(wǎng)絡模型需要大量參數(shù)，學習時間過長，不能觀察之間的學習過程；貝葉斯模型對數(shù)據(jù)的表達形式過于敏感；決策樹［14］則有過擬合的缺點。不同學者各自通過隨機森林模型對不同研究區(qū)的土壤養(yǎng)分含量進行預測，結果均表明隨機森林模型的預測精度十分可靠，是數(shù)字土壤制圖的理想方式［15-16］。相對于其他機器學習方法，隨機森林可以輸入大量變量，且學習過程較快，同時當部分數(shù)據(jù)缺失時，仍然可以維持相當高的精準度。因此，本研究采用了隨機森林來訓練模型，探索隨機森林在數(shù)字土壤制圖中的應用效果，同時與傳統(tǒng)插值模型中最為常見的普通克里金和反距離權重方法作比較，旨在準確預測安慶市耕地土壤速效鉀空間分布模式。

1 材料與方法

1.1 研究區(qū)概況

安慶，又稱宜城，轄二市五縣三區(qū)，位于皖西南部，是皖西南區(qū)域中心城市。安慶市屬于亞熱帶季風性濕潤氣候，年平均降水量為1368 mm，年平均溫度為16.5℃。地勢整體西北高，東南低，呈山地——丘陵崗地——平原階梯狀分布。土壤母質類型豐富，主要成土母質為酸性結晶巖類風化物、山河沖積物、泥頁巖類風化物；土壤類型多樣，以水稻土、黃棕壤、黃褐土、潮土為主。

1.2 采樣點數(shù)據(jù)來源

本研究采樣點數(shù)據(jù)來源于2019年安慶市耕地質量等級評價取樣結果，布點要求平原地區(qū)耕地每667 hm2不少于一個點位，山地地區(qū)耕地每133～150 hm2一個點位，并盡可能分布至全市各鄉(xiāng)鎮(zhèn)、土類和成土母質。采樣按照五點取樣法，采取0～20 cm耕地表層土壤化驗，全市共獲取682個采樣點位，隨機選擇15%（103個點位）的樣點作為驗證樣點，85%（579個點位）的樣點作為測試點位。具體分布情況如圖1所示。

1.3 環(huán)境變量因子來源與選取

研究區(qū)地形復雜多樣，地勢起伏巨大，海拔落差峰值達到1700余m，地貌類型主要包括西北大別山區(qū)、中部丘陵崗地區(qū)、東南長江沿江平原區(qū)；成土母質類型豐富，土壤類型多樣；氣候為亞熱帶季風性濕潤氣候，氣候溫和、雨量充沛，但因地形地貌不同，研究區(qū)各區(qū)域氣溫和降水量差異較大，直接影響巖石風化、耕地土壤礦物質的移動和分布，進而直接或間接影響研究區(qū)耕地土壤速效鉀的分布情況。本研究選取的環(huán)境變量因子主要包括空間位置變量因子、地形變量因子、土壤變量因子和氣候變量因子。其中，空間位置變量因子包括經(jīng)度和緯度；地形變量因子包括高程（圖2）、坡度（Slope）、坡向（Aspect）、剖面曲率（Profile curvature）、平面曲率（Plane curvature）和地形起伏度，數(shù)據(jù)來源于地理國監(jiān)測云平臺（http：//www.dsac.cn）獲取的30 m分辨率研究區(qū)數(shù)字高程模型（Digital elevation model），坡度、坡向、剖面曲率、平面曲率，地形起伏度數(shù)據(jù)運用了ArcGIS 10.2的空間分析工具（Spatial Analyst），根據(jù)數(shù)字高程模型計算得來；土壤變量因子包括成土母質和土壤類型，成土母質（圖3）與土壤類型圖，數(shù)據(jù)來源于安慶市第二次全國土壤普查紙質圖件，經(jīng)ArcGIS 10.2矢量化形成；氣候變量因子包括年平均降水量（圖4）和年平均溫度（圖5），數(shù)據(jù)來源于全球氣候數(shù)據(jù)網(wǎng)（http：//www.worldclim.org），獲取了分辨率為1000 m的1970～2000年的年平均降水量及年平均溫度數(shù)據(jù)。

1.4 空間預測方法

1.4.1 隨機森林

隨機森林（Random Forest，RF）是機器學習中常用的一種方法，最早由Breiman于20世紀80年代提出［17］。模型運用bootstrap隨機重復選取樣得到預測集，并對每個預測集組合決策樹，得到?jīng)Q策樹數(shù)目（ntree）。在各個決策樹節(jié)點，選擇樣本預測器的數(shù)目（mtry）作為分割變量，輸入矩陣值，通過匯總分類樹的結果組合成隨機森林，并以此預測變量。公式如下：

本研究采用MATLAB2015 R2015b軟件編程，選取地形變量因子（高程、坡度、坡向、剖面曲率、平面曲率和地形起伏度）、土壤變量因子（成土母質、土壤類型）、氣候變量因子（年平均降水量、年平均溫度）以及空間位置變量因子（經(jīng)度、緯度），共計12個變量因子作為預測變量。一般來說，mtry取值為變量因子總數(shù)三分之一［18］，ntree一般選取默認值為500［19］，但當ntree的增加不能顯著提升預測精度時，ntree需要盡可能小［20］。本研究通過3組12次試驗（ntree=200，300，400，500；mtry=2，3，4），對比各試驗組中驗證樣點的R2值。其中，R2越接近1，擬合精度越高。結果表明，在試驗組中，當ntree=400，mtry=3時，研究區(qū)RF預測精度最高。

表1 RF模型中不同ntree和mtry的驗證樣點R2統(tǒng)計

1.4.2 普通克里金

克里金插值法（Kriging）起源悠久，南非礦產(chǎn)工程師Danie G. Krige在20世紀50年代初開創(chuàng)性地采用此方法找尋金礦［21］，法國統(tǒng)計學家Georges Matheron為此方法奠定了基礎，并命名為Kriging［22］。普通克里金（Ordinary Kriging，OK）插值法作為地統(tǒng)計學的重要內容之一，在空間自相關基礎上，利用半變異函數(shù)對區(qū)域內樣點進行預測分析，OK要求原始數(shù)據(jù)或經(jīng)過函數(shù)變換后的數(shù)據(jù)正態(tài)分布。其公式為：

式中，s表示不同空間位置的坐標點位，可默認為經(jīng)緯度坐標；是s處的變化量；表示趨勢值；表示自相關誤差。

1.4.3 反距離權重

反距離權重（Inverse Distance Weighting，IDW）插值法是基于地理學第一定律的原理：即萬事萬物都是有聯(lián)系的，兩個物體距離越近，其性質越接近；兩個物體距離越遠，其性質差異性越明顯［23］。IDW插值法作為運用最為廣泛和普遍的空間插值方法，通過已知點和預測點之間的距離作為權重并加權平均計算，離已知點越近的預測點權重越大。IDW插值法要求插值對象分布均勻適中，并且插值面有局部因變量。其公式為：

1.4.4 研究區(qū)模型精度檢驗

研究區(qū)樣點由85%預測樣點（579個）和15%驗證樣點（103個）組成，模型精度檢驗是對研究區(qū)103個驗證樣點的誤差估計，通過計算驗證樣點的平均絕對誤差（MAE）、均方根誤差（RMSE）和決定系數(shù)（R2）確定各模型預測效果，其中，MAE、RMSE越小，R2越接近1，模型預測精度越高，預測效果越好。具體計算公式如下：

2 結果與分析

2.1 耕地土壤速效鉀描述性統(tǒng)計

研究區(qū)采樣點隨機劃分為預測樣點組和驗證樣點組。從表1可知，預測樣點速效鉀值預測范圍為4 ～ 691 mg·kg-1，驗證樣點速效鉀值范圍是23 ～349 mg·kg-1，平均值、中值和標準差預測樣點均略大于驗證樣點，預測樣點存在少量極端極大值。變異系數(shù)分別為74.1%和69.8%，均屬于中等變異水平。預測樣點偏度和峰度分別為2.58和11.76，經(jīng)過Box-Cox變換，同時經(jīng)過K-S檢驗后符合正態(tài)分布，變換參數(shù)為0.1。

2.2 研究區(qū)模型精度評價

通過對比研究區(qū)103個驗證樣點與3種模型預測結果，MAE、RMSE和R2，確定各模型精度。其中，RF的MAE、RMSE、R2分別為30.93 mg·kg-1、41.31 mg·kg-1和0.58，OK的MAE、RMSE、R2分別 為31.97 mg·kg-1、44.08 mg·kg-1和0.49，IDW的MAE、RMSE、R2分別為32.77 mg·kg-1、46.21 mg·kg-1和0.47。從表3可以看出，在3種方法中，RF的MAE和RMSE均為最小，R2最大，RF的MAE較其余兩種方法分別減少了1.04和1.84 mg·kg-1，RMSE較其余兩種方法分別減少了2.77和4.9 mg·kg-1，R2較其余兩種方法分別增加了0.09、0.11；OK的MAE和RMSE、R2均 為 第 二，OK的MAE較IDW減少了0.8 mg·kg-1，RMSE較IDW減少 了2.13 mg·kg-1，R2較IDW增 加 了0.02；IDW的MAE和RMSE均為最大，R2最小?？梢缘贸鼋Y論：研究區(qū)內速效鉀空間分布的3種方法的預測精度高低順序分別為RF＞OK＞IDW。

表3 研究區(qū)3種預測方法的耕地土壤速效鉀預測精度統(tǒng)計

通過篩選速效鉀極大值（＞200 mg·kg-1），發(fā)現(xiàn)OK的預測結果中的極大值最低，由表4可得其R2為0.13，MAE為93.92 mg·kg-1，RMSE為105.13 mg·kg-1，IDW雖然因“牛眼”效應（某些極值點在插值結果上所形成的次插值點為圓心的圓斑，且圓斑的值與該極值點較為接近）使得極大值較大，但預測精度在3種方法中仍然最低，其R2僅為0.0002，MAE為107.64 mg·kg-1，RMSE為111.37 mg·kg-1。RF的極大值較大，其R2為0.30，MAE為89.90 mg·kg-1，RMSE為101.90 mg·kg-1，可 以 看出，RF對極大值的預測精度較另外2種有較大提高，但極大值部分的預測精度較整體的預測精度低。

表4 研究區(qū)3種預測方法的耕地土壤速效鉀極大值（＞200 mg·kg-1）預測精度統(tǒng)計

2.3 隨機森林環(huán)境變量因子權重

為確定各個環(huán)境因子在RF模型中的權重，通過移除某個環(huán)境變量因子，核查RF模型RMSE增加的比例，從而判定各個環(huán)境因子的權重及其排序。從表5中可以看出，緯度、年平均溫度、成土母質、高程、經(jīng)度、年平均降水量是影響研究區(qū)RF模型下速效鉀含量預測的重要因素，平面曲率、剖面曲率、地形起伏度、坡度、坡向與土壤類型權重較低。前5個地形變量因子與高程有相關性，土壤類型與成土母質有相關性，可以得知，此類協(xié)同變量因子在RF模型中的權重會因共線性而導致其權重降低，這與任麗等［24］的研究結果一致。經(jīng)度、緯度通過空間位置的變化對速效鉀含量產(chǎn)生影響。成土母質是土壤形成的基礎，是速效鉀含量的初始影響因素。年平均溫度和年平均降水量則通過風化作用，崩解、遷移成土母質中的化學物質，從而影響速效鉀含量。

表5 研究區(qū)隨機森林環(huán)境變量因子權重統(tǒng)計

2.4 耕地土壤速效鉀空間分布預測結果

對比3種方法下研究區(qū)耕地土壤速效鉀空間分布預測結果（圖6），可以看到整體分布趨勢接近，研究區(qū)耕地土壤速效鉀含量在東南沿江處較高，同時研究區(qū)中南部出現(xiàn)小塊高值區(qū)，其原因在于兩部分區(qū)域地貌類型為平原，高程較低，成土母質以長江沖積物為主，易于速效鉀累積。OK和IDW的低值區(qū)較為接近，主要分布在研究區(qū)的北部及中部區(qū)域，RF的低值區(qū)則主要分布在研究區(qū)中部區(qū)域。研究區(qū)中部區(qū)域在3種空間預測方法中均為低值區(qū)，其主要原因是受成土母質影響，成土母質主要是紫色巖類風化物和山河沖積物。這部分區(qū)域母質以風化作用為主，風化物松散，膠結性低，易遭侵蝕，風化淋溶較弱，發(fā)育的土壤主要屬于幼年土壤，剖面發(fā)育差，剖面構型以耕作層（A）-底土層（C）為主，土體淺薄，僅10 cm左右，不利于速效鉀積累。研究區(qū)北部區(qū)域主要受成土母質和地貌影響。成土母質基本為酸性結晶巖類風化物，這類母質發(fā)育的土壤質地多為壤土、土體厚度相對較厚，保水保肥能力較研究區(qū)中部區(qū)域強。地貌類型為山地，山區(qū)坡度較大，地形變化復雜，易發(fā)生水土流失。綜合成土母質與地貌影響，在RF空間分布預測圖中，研究區(qū)北部區(qū)域耕地土壤速效鉀含量不高，但并不處于低值區(qū)。

3 討論

通過篩選速效鉀極大值（＞200 mg·kg-1），發(fā)現(xiàn)其極大值部分的R2仍低于整體值的R2，極大值部分的MAE和RMSE則大于整體值的MAE和RMSE，說明極大值部分的預測精度較整體的預測精度低。推斷可能在于局部區(qū)域受人為因素影響導致局部速效鉀含量增加，諸如人為施肥、農(nóng)田管理等。而本研究只選取了環(huán)境變量因素，并未考慮選取反映局部變量的人為因素作為輔助變量因子，后期可以通過增加一些能夠反映局部信息的非環(huán)境變量因子來提升RF的擬合效果。

通過3種空間預測方法的極大值和整體精度比較，可以得知，RF的預測結果能夠精確地反映研究區(qū)耕地土壤速效鉀分布情況，其精度最高。其本質在于RF屬于多變量因素分析，OK和IDW屬于單變量因素分析，而耕地土壤速效鉀含量是受多種環(huán)境變量因素綜合影響［25］。因此，OK和IDW僅僅從空間位置關系上預測速效鉀含量的算法限制了其預測精度。RF的整體分布趨勢與OK、IDW接近，但RF通過大量的決策樹和樣本預測器構造模型，在空間位置關系的基礎上，充分運用環(huán)境變量因子，很好地刻畫了研究區(qū)耕地土壤速效鉀含量與環(huán)境變量因子的非線性關系，所以RF的預測精度在3種空間預測方法中最高。

4 結論

對比RF、OK、IDW 3種方法對安慶市耕地土壤速效鉀含量空間分布預測結果，3種方法整體趨勢一致。

從預測精度上看，RF的MAE和RMSE均小于OK和IDW，RF的R2大于OK和IDW；OK的MAE和RMSE、R2均為第二；IDW的MAE和RMSE均為最大，R2最小。說明RF在研究區(qū)范圍內的預測精度最高，OK次之，IDW預測精度最低，OK和IDW只能在一定趨勢上反映出速效鉀空間分布的特點。其原因在于OK和IDW僅僅以空間位置作為空間預測的變量因子，RF則充分考慮了環(huán)境變量因子，故精度最高。

通過移除某個環(huán)境變量因子，驗證RF的RMSE變化比例，可以得知12個變量因子中，緯度、年平均溫度、成土母質、高程、經(jīng)度、年平均降水量是影響RF精度的主要因素，剩余變量因子在RF模型中的權重會因協(xié)同重復而導致其權重降低。