基于三種空間預測方法的安慶市耕地土壤速效鉀空間分布預測

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

（1.安徽農業大學資源與環境學院，安徽 合肥 230036；2.安慶市種植業管理局，安徽 安慶 246000）

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

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

1 材料與方法

1.1 研究區概況

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

1.2 采樣點數據來源

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

1.3 環境變量因子來源與選取

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

1.4 空間預測方法

1.4.1 隨機森林

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

本研究采用MATLAB2015 R2015b軟件編程，選取地形變量因子（高程、坡度、坡向、剖面曲率、平面曲率和地形起伏度）、土壤變量因子（成土母質、土壤類型）、氣候變量因子（年平均降水量、年平均溫度）以及空間位置變量因子（經度、緯度），共計12個變量因子作為預測變量。一般來說，mtry取值為變量因子總數三分之一［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時，研究區RF預測精度最高。

表1 RF模型中不同ntree和mtry的驗證樣點R2統計

1.4.2 普通克里金

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

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

1.4.3 反距離權重

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

1.4.4 研究區模型精度檢驗

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

2 結果與分析

2.1 耕地土壤速效鉀描述性統計

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

2.2 研究區模型精度評價

通過對比研究區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最小。可以得出結論：研究區內速效鉀空間分布的3種方法的預測精度高低順序分別為RF＞OK＞IDW。

表3 研究區3種預測方法的耕地土壤速效鉀預測精度統計

通過篩選速效鉀極大值（＞200 mg·kg-1），發現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 研究區3種預測方法的耕地土壤速效鉀極大值（＞200 mg·kg-1）預測精度統計

2.3 隨機森林環境變量因子權重

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

表5 研究區隨機森林環境變量因子權重統計

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

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

3 討論

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

通過3種空間預測方法的極大值和整體精度比較，可以得知，RF的預測結果能夠精確地反映研究區耕地土壤速效鉀分布情況，其精度最高。其本質在于RF屬于多變量因素分析，OK和IDW屬于單變量因素分析，而耕地土壤速效鉀含量是受多種環境變量因素綜合影響［25］。因此，OK和IDW僅僅從空間位置關系上預測速效鉀含量的算法限制了其預測精度。RF的整體分布趨勢與OK、IDW接近，但RF通過大量的決策樹和樣本預測器構造模型，在空間位置關系的基礎上，充分運用環境變量因子，很好地刻畫了研究區耕地土壤速效鉀含量與環境變量因子的非線性關系，所以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在研究區范圍內的預測精度最高，OK次之，IDW預測精度最低，OK和IDW只能在一定趨勢上反映出速效鉀空間分布的特點。其原因在于OK和IDW僅僅以空間位置作為空間預測的變量因子，RF則充分考慮了環境變量因子，故精度最高。

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