不同方法預測蘇南農田土壤有機質空間分布對比研究*

謝恩澤 趙永存? 陸訪儀 史學正 于東升

（1 土壤與農業可持續發展國家重點實驗室（中國科學院南京土壤研究所），南京 210008）

（2 中國科學院大學，北京 100049）

土壤有機質（Soil organic matter, SOM）是土壤的重要組成部分，盡管只占土壤總質量的很小一部分，但其作為土壤肥力形成的基礎，不但影響土壤質量及關鍵功能，而且在陸地生態系統物質循環中扮演重要角色［1-3］。因此，利用不同空間預測方法對比預測農田SOM，對于準確把握SOM的空間分布規律和土壤資源管理的科學決策，進而保障糧食安全、氣候安全及土壤安全具有重要意義。

地統計學方法廣泛用于土壤屬性空間預測［4］。而克里格法是其中最準確，強大的插值方法之一［5］。徐尚平等［6］基于半方差函數和普通克里格方法分析了內蒙古地區SOM含量的空間結構特征；Mushtaq等［7］采用普通克里格方法描繪了克什米爾農業區土壤微量元素的空間分布；Heba等［8］則利用普通克里格方法對尼羅河三角洲地區的土壤碳庫和氮庫的空間分布特征進行了深入分析。這些結果表明，克里格方法在不同區域、不同土壤性質的空間變化中均能得到較好的預測結果。遙感和數字制圖技術的快速發展，為空間預測提供了大量廉價且與土壤形成過程直接或間接相關的輔助數據，比如，DEM、植被指數及生物量、光譜反射率、數字化的土壤圖、土地利用/覆蓋圖等。這些輔助數據提供了關于土壤性質變異性的先驗信息，從而使得整合輔助數據先驗信息的回歸克里格等混合建模方法得到了快速發展。Hengl等［9］以土壤圖和DEM獲取的地形因子為輔助數據，采用回歸克里格方法預測了克羅地亞SOM的空間分布；Andrew等［10］則采用回歸克里格方法整合紅外光譜數據，大大提高了土壤全碳、氮等土壤屬性的預測精度。回歸克里格方法不但反映了輔助因子對土壤屬性空間分異的潛在影響，同時，趨勢項的分離和殘差的再預測也使得其空間預測精度通常要優于普通克里格方法。近年來，隨機森林作為經典的機器學習模型之一，由于具有處理高維度數據、防止過度擬合、提供輔助因子的重要性度量等優勢，在土壤空間預測，特別是地形分異明顯地區的土壤空間預測中也逐漸得到廣泛應用。比如，Mareike等［11］在厄瓜多爾南部高程起伏1.7～3.2 km的安第斯山脈典型區中，利用隨機森林方法成功地預測了土壤質地的空間分布；Dharumarajan等［12］以印度南部安得拉邦阿嫩達布爾區的Bukkarayasamudrum mandal為研究區，在高程（295～595 m）較復雜多變的情況下利用地形、植被等輔助數據，采用隨機森林方法預測了土壤有機碳、pH和電導率的空間分布；而Zhi等［13］則通過隨機森林方法并整合DEM獲取的地形因子、遙感影像獲取的光譜反射率及植被指數等輔助數據預測了青藏高原地區高山草氈層的空間分布。隨機森林方法在地形復雜且人為活動影響相對較弱的區域土壤空間預測中得到了較為成功的應用，而在地形平坦、人為影響強烈的地區應用效果則依然不清楚。

上述研究可以看出，不同預測方法利用輔助因子在平原地區的應用仍需進一步探討。所以，選擇地勢平坦的江蘇南部（蘇南）為研究區，以影響土壤生態系統服務功能的關鍵因子SOM為研究對象，以輔助因子與SOM的相關性強弱及輔助因子的可獲取性為切入點，通過回歸克里格和隨機森林方法與傳統的普通克里格方法對比分析來研究不同空間預測方法對該區SOM空間分布預測的影響，以便為準確預測該區SOM空間分布規律，進而為通過科學調控不同區域SOM來實現關鍵土壤生態系統服務功能提供科學依據。

1 材料與方法

1.1 研究區概況

研究區為江蘇省南部地區（30°44′～32°16′N，118°29′～121°19′E）。該區屬北亞熱帶季風氣候，四季分明，雨量充沛，年均溫16℃，年均降雨量1 000 mm。全區地勢西高東低，太湖平原地區地勢平坦，河流湖泊較多，主要土壤類型為水稻土，沿江地區則以灰潮土為主。研究區西部及西南部多為丘陵、低山崗地，土壤類型主要為黃棕壤。研究區農耕條件優越，工業經濟發展迅速，城鎮化水平較高，是長三角地區經濟最為發達的地區之一。

1.2 數據來源與處理

土壤采樣點位置根據該區20世紀80年代第二次土壤普查結果進行布設，共413個（圖1），采樣時間為2000年，采樣深度為0～20 cm。采樣同時記錄了樣點的GPS位置信息以及環境條件和主要農業管理措施信息。土壤有機質（SOM）含量采用重鉻酸鉀氧化-外加熱法測定，全氮采用半微量開氏法測定，有效磷采用碳酸氫鈉浸提-鉬銻抗比色法測定，速效鉀采用乙酸銨浸提-火焰光度法測定，pH采用電位法（土︰水=1︰2.5）測定，機械組成采用吸管法測定［14］。其中，土壤全氮、有效磷、速效鉀、pH及黏粒含量主要用于SOM空間預測的輔助土壤屬性因子。

圖1 蘇南地區地形及土壤采樣點分布圖Fig. 1 Topographic map of South Jiangsu and distribution of soil sampling sites

SOM空間預測的其他輔助因子數據包括：20世紀80年代初始SOM含量（Initial SOM,SOM1980），來源于研究區各縣/區的土壤普查報告；江蘇省1︰20萬土壤圖；江蘇省30 m分辨率D E M計算的主要地形因子，包括高程（Elevation, Elev）、坡度（Slope）、平面曲率（Plane curvature, Planc）、剖面曲率（Profile curvature, Profc）和復合地形指數（Compound terrain index, CTI）；1 km×1 km的年均氣溫（Mean annual temperature, Mat）和降雨量（Mean annual precipitation, Map）數據（1976—2005年平均，來源于碳專項）。除了上述常用于SOM空間輔助預測的因子外，本研究也考慮了一些容易獲取且免費的數據產品，以檢驗其在提高SOM空間預測精度中的可行性，這些因子包括：土壤溫度（Soil temperature, ST）和土壤濕度（Soil moisture,SM）1 km×1 km柵格數據（2000年月平均，分別來源于歐洲中期天氣預報中心（ECMWF）和Earthdata）；遙感反演的5 km×5 km凈初級生產力（Net Primary Production, NPP）數據（1980—2000年平均，來源于北京師范大學）；5 km×5 km年均氮肥用量（Nitrogen fertilizer, N-fert）柵格數據（1994—2001年平均，來源于美國國家航空航天局（National Aeronautics and Space Administration,NASA）社會經濟數據應用中心）。其中遙感NPP主要用于估算農田作物還田碳輸入。農作物碳投入主要包括：秸稈還田、作物根茬殘留和作物地下生物量。根據周睿等［15］區域農田有機物質的估算方法利用NPP數據進行農田碳投入的碳估算：

Btotal=α×NPP式中，Btotal為作物整株生物量（包括地上和地下），單位為g·m-2（生物量按C計算）；NPP為凈初級生產力；α為碳-干物質系數。由作物的根冠比、收獲指數，秸稈還田比例計算農田碳投入。

對于土壤SOM1980、TN、Ap、Ak、pH和Clay等點位測定數據，采用普通克里格方法預測提取300 m×300 m柵格數據，插值后的柵格數據作為SOM空間預測的輔助因子。對于30 m分辨率DEM計算得到的Elev、Slope、Planc、Profc和CTI等輔助預測因子，采用鄰域平均法將尺度降至300 m×300 m柵格，以消減土壤樣點可能存在的GPS定位誤差對空間預測的影響。而對于空間分辨率相對較低的ST、SM、和N-fert等輔助預測數據則采用普通克里格方法插值將尺度升至300 m×300 m柵格作為SOM空間預測的輔助因子。

1.3 SOM空間預測輔助因子的篩選

首先，采用皮爾森相關分析對SOM空間預測的輔助因子進行初步篩選。表1相關分析的結果表明，Clay、SOM1980、TN、Ak、pH、SD、CT、Elev、Profc、Mat、Map、InputC、SM、ST、N-fert共15個輔助預測因子與SOM含量顯著相關，可潛在地用于SOM含量空間分布預測。

表1 土壤有機質含量與輔助預測因子的相關系數Table 1 Pearson correlation coefficients between soil organic matter content and auxiliary factors

其次，對于皮爾森相關分析初步篩選的15個輔助預測因子，進一步采用逐步回歸分析進行篩選，以便避免因子間的共線性效應。同時，為探討輔助預測因子與SOM相關性強弱對不同空間預測方法預測精度的潛在影響，在逐步回歸分析過程中，考慮了保留相關性最強的輔助預測因子TN和移除TN兩種情況以便進行對比分析。逐步回歸的篩選主要是以最小信息準則（AIC）為衡量標準，通過選擇最小的AIC信息統計量，來達到刪除變量的目的，直到當去掉某一因子時AIC均增加，逐步回歸分析終止，達到最優模型。由此得到，保留相關性最強的輔助預測因子TN的條件下，逐步回歸篩選得到的SOM輔助預測因子分別為：TN、SOM1980、SD、Mat、SM、Clay；移除與SOM相關性最強的TN預測因子的條件下，逐步回歸保留的輔助預測因子包括SD、Elev、Mat、InputC、Map、ST和SOM1980共七個因子（表2）。

表2 土壤有機質（SOM）預測輔助因子的逐步回歸分析Table 2 Stepwise regression analysis of auxiliary factors for prediction of soil organic matter （SOM）

1.4 SOM的空間預測方法

普通克里格（Ordinary kriging，OK）：基于土壤屬性的空間自相關性，通過鄰近的相關觀測點的權重均值來預測采樣點位置的土壤屬性［16］。本文通過SOM含量半方差函數及擬合參數，采用OK方法預測SOM空間分布。

回歸克里格（Regression kriging，RK）：首先基于逐步回歸分析保留的輔助預測因子，采用多元線性回歸（MLR）對SOM進行預測，然后對多元回歸預測的SOM殘差進行OK預測，最后將回歸預測值與殘差的OK預測值相加得到SOM的預測值［8］。

隨機森林（Random forest，RF）：首先基于逐步回歸篩選的預測因子和SOM觀測數據，對隨機森林模型中的關鍵參數每次分裂向量數（mtry）進行參數優化，然后采用經參數優化后的RF模型對SOM空間分布進行預測。參數優化過程中，分裂節點隨機向量數量ntree設定為1 000，以預測誤差最小化為目標函數，通過遍歷比較確定mtry的最優值［17］（根據預測因子的個數設置mtry參數為1～7，并進行建模），以運行100次的平均RMSE結果尋找最優參數。結果表明，當保留與SOM相關性最強的預測因子TN時，mtry最優值為4，而移除相關性最強的輔助因子TN時,mtry最優值為1。

1.5 空間預測精度分析

SOM的空間預測精度采用獨立驗證進行評價。為了使模型每次驗證的結果覆蓋到整個樣點數據集，將樣點數據集分成預測集和驗證集進行三折獨立驗證（即將樣點數據隨機分成三份，分三次取其中兩份為預測數據集，一份為驗證數據集），循環執行三折獨立驗證100次，通過比較100次得到的驗證數據點位置上SOM均方根誤差進行空間預測精度分析。其中，均方根誤差（RMSE）的計算公式為：

式中，zi為第i個樣點的預測值；zi為第i個樣點的實際值；n為樣點個數。RMSE的結果越小表明預測精度越高。

本研究中輔助數據的篩選、尺度轉換等數據處理均采用R軟件進行分析，其中OK、RK和RF及空間預測精度分析采用R語言編程計算，SOM空間分布采用Arcmap10.2軟件進行制圖。

2 結 果

2.1 SOM含量的描述性統計特征

研究區413個樣點的描述統計分析（圖2）結果表明，農田土壤耕層SOM平均含量為29.29 g·kg-1，中值與均值較為接近，為28.79 g·kg-1，高于江蘇北部的黃淮平原和東部濱海平原，略高于江淮平原地區［18］；標準差為7.18 g·kg-1，變異系數（CV）為25%，變異程度中等；SOM含量介于7.75～58.31 g·kg-1之間，極差為50 g·kg-1。SOM含量的頻率分布直方圖表現為稍向均值右側偏離，偏度系數為0.35，峰度值為0.75，表現為近似正態分布。

圖2 土壤有機質含量的描述統計及頻數分布直方圖Fig. 2 Descriptive statistics and histogram of soil organic matter（SOM） contents

2.2 不同空間預測方法的預測精度

從不同空間預測方法100次三折獨立驗證的結果（圖3）來看。OK方法的預測誤差最大，空間預測的均方根誤差（RMSE）均值為6.97 g·kg-1，僅能解釋SOM方差的9.7%；當整合與SOM相關性最強的輔助預測因子TN后，RK和RF預測方法的RMSE均值分別降低至5.25 g·kg-1和4.97 g·kg-1，分別降低了25%和29%，同時，RK和RF方法預測結果對SOM方差的解釋程度分別提高到了51%和55%，RF方法均方根誤差的四分位距亦較RK方法小，表明整合與SOM相關性最強的輔助預測因子TN時，RF方法預測的穩健性更好。然而，當不整合與SOM相關性最強的輔助預測因子TN時，進入空間預測過程的其余輔助預測因子與SOM的相關性減弱（表1、表2），RK和RF方法預測的RMSE均值也分別增加至6.21 g·kg-1和6.29 g·kg-1，相應的SOM方差解釋程度亦分別降低至29%和28%，同時，RK方法均方根誤差的四分位距也較RF方法稍小，表明在不整合與SOM相關性強的輔助因子TN時，RK方法預測的精度及穩健性較RF方法稍好。此外，盡管不整合與SOM相關性最強的輔助預測因子TN時，SOM的預測精度要較整合輔助預測因子TN時低，但其預測精度依然較OK方法高。

圖3 不同空間預測方法的預測誤差對比Fig. 3 Comparison between the spatial prediction methods in deviation

2.3 SOM的空間分布特征

從不同方法預測的SOM含量空間分布圖來看（圖5），SOM含量空間分布的總體趨勢大體一致，主要表現為研究區西部的丹陽-金壇-溧陽一線含量相對較低，SOM含量變化范圍集中在15～35 g·kg-1，而研究區東部到東南部的常熟-昆山-吳江一線SOM含量相對較高,大多數為25～45 g·kg-1。OK方法預測的SOM含量空間分布模式表現為具有較多的、局部分布的島狀高值區和低值區，且島狀區域分布較為零散。與OK方法相比，RK和RF兩種方法預測的SOM空間分布的低值區連續性稍好，同時，昆山南部和吳江南部的SOM高值區也得到了更好的表達。然而，考慮到是否整合與SOM相關性最強的輔助預測因子時，RK和RF兩種方法預測的SOM空間分布也存在一定的差異性。當在空間預測過程中整合了與SOM相關性最強的輔助預測因子（比如TN）時，除了RF方法能夠進一步表達太倉東部地區的SOM低值區外，RK和RF兩種方法預測的SOM空間分布基本完全一致；當移除與SOM相關性最強的輔助預測因子TN后，其余輔助預測因子與SOM的相關性較弱，此時，與RK方法預測的SOM空間分布模式相比，RF預測結果中，研究區西部丹陽-金壇-溧陽一線的SOM含量低值區范圍變大，同時研究區東部的SOM高值區的范圍亦明顯增加。

3 討 論

3.1 不同方法的預測精度

在樣點數量及樣點空間位置不變的條件下，OK方法的空間預測精度在很大程度上取決于待預測目標土壤屬性的空間變異結構［19］。從SOM含量半方差函數的擬合結果來看（表3），研究區SOM含量的空間相關性較弱，其中塊金值與基臺值的比值高達83%，表明耕作、施肥等隨機因素的強烈影響削弱了SOM的空間自相關性［20］，可能是導致OK方法預測精度低的主要原因之一。

整合與SOM具有一定相關性的輔助預測因子后，RK和RF方法均能在一定程度上提高SOM的空間預測精度。Knotters［21］和Ahmed［22］研究結果表明，在插值方法中使用輔助因子有兩個重要特性需要考慮，首先是輔助因子和目標變量之間的相關性強弱，其次是預測的殘差是否具有一定的空間結構性。當整合與SOM相關性較強的輔助預測因子（比如TN）后，RF的預測精度顯著提高，表明RF方法預測精度主要受目標變量與輔助預測因子間的相關性強弱所控制。而對于RK方法而言，整合高相關性的輔助預測因子的預測精度也較不整合時要高，但由于整合高相關性的輔助預測因子后分離出的殘差缺乏空間相關性而導致RK方法的預測精度低于RF方法。

圖4 不同方法預測SOM空間分布Fig. 4 SOM spatial distribution maps based on prediction using different prediction methods

表3 土壤有機質含量的半方差函數模型及擬合參數Table 3 Semi-variogram model and fitting parameters of soil organic matter content

在不整合高度相關的輔助預測因子TN的條件下，RK方法預測精度仍然較RF方法稍高，表明RK方法的預測精度除了與目標變量（比如本研究中的SOM）和輔助因子間的相關性強弱有關外［23］，還在很大程度上受目標變量多元回歸（MLR）預測殘差的空間變異結構影響［24］。比如，從整合與不整合輔助預測因子TN的MLR殘差的半方差函數對比來看（表3），在輔助數據與SOM相關性較強的條件下（整合TN后），由于SOM的大部分趨勢組分被MLR所分離，導致SOM的殘差半方差函數基本為純塊金效應，即趨勢性組分預測精度高而殘差預測精度降低；而在輔助數據與SOM的相關性較弱的條件下（不整合TN），盡管輔助數據MLR分離的SOM趨勢組分相對較少，但SOM殘差的半方差函數相比于整合TN條件下依然具有一定的結構性，從而殘差的預測精度能夠有所提升。從數據來源分析，由于土壤全氮含量獲取成本較高，而移除全氮因子后回歸模型中保留的氣象、DEM、土壤圖及初始SOM含量等輔助預測數據（表1）則相對較為容易獲取，因此，從這一點來看，在輔助預測數據與SOM有相關性，但相關性相對較弱的條件下，傳統的RK方法在本研究區SOM空間預測上仍然具有較高的應用價值。

3.2 不同方法預測結果局限性

本研究區地勢平坦，DEM獲取的地形因子在SOM預測應用中還有一定的局限性，同時，本研究區土壤受人為活動影響強烈，目前容易獲取、相對廉價且可能用于輔助SOM空間預測的數據在提高本研究區SOM空間預測精度方面還面臨諸多挑戰。比如，相對于本研究區的面積（2.7萬km2）而言，盡管土壤濕度、溫度、氮肥施用量等輔助數據比較容易獲取，但依然存在空間分辨率相對較低，難以和SOM觀測點數據進行尺度匹配等問題。Li［25］利用RK方法預測SOM空間分布的研究結果表明，目標變量和輔助因子的空間尺度匹配性直接影響SOM的空間預測精度。因此，如何對現有容易獲取且廉價但相對粗糙的輔助數據進行尺度轉換以便更好地識別輔助數據與SOM間潛在的相關性是未來需要進一步研究的問題。此外，對本研究區而言，SOM除了受土壤類型等環境因子影響外，SOM含量特別是表層SOM含量還受施肥、耕作、秸稈管理等農業管理措施的強烈影響［26］。作物殘茬的C輸入是影響農田土壤SOM含量的重要因素之一［27］，整合土壤C輸入來預測SOM的空間分布不但能在一定程度上間接地反映農業管理措施等人為活動對SOM空間分布區域差異性的潛在影響，也可能在一定程度上提高SOM空間預測精度。然而，從本研究區C輸入與SOM間的相關性分析來看，其直接的相關性還相對較弱。這不但與NOAA AVHRR遙感數據估算的NPP空間分辨率相對較低且存在混合像元等問題有關，也與研究區的農田土壤C輸入來源相對復雜（比如，除了根茬C輸入外還有有機肥等C輸入等）有關。因此，在平原地區人為活動因子影響強烈的條件下，尋找與SOM相關性強的、容易獲取且廉價的其他輔助預測替代因子，或者開發適宜的輔助數據降尺度方以實現輔助數據與SOM實測點數據的尺度匹配，將是實現平原地區強烈人為作用影響下SOM空間有效預測的重要途徑之一。

4 結 論

OK、RK和RF三種方法對蘇南地區農田表層SOM含量空間預測結果的對比分析表明：（1）三種方法預測的SOM空間分布總體趨勢大體一致，總體表現為東高西低，但SOM空間分布的局部細節存在一定的差異性。從預測精度來看，OK方法預測的平均RMSE為6.97 g·kg-1，RK和RF在保留相關性最強的輔助因子TN時平均RMSE分別為5.25 g·kg-1和4.97 g·kg-1，而在移除相關性最強的輔助因子TN后，平均RMSE分別為增加至6.21 g·kg-1和6.29 g·kg-1，因此，RK和RF預測精度均較OK要高，且RK和RF均隨輔助因子與SOM的相關性變強而預測精度增加。（2）RK和RF的預測精度存在一定的差異，表現為保留相關性最強的TN時，RF的預測精度較RK要高，但移除TN后，RK較RF精度稍高。因此，考慮到TN的獲取成本及移除TN后進行預測模型的輔助因子可獲取性問題，RK方法在輔助因子與SOM相關性較弱的條件下對本研究區依然有較高的應用價值。（3）RK和RF的預測精度均與輔助因子和SOM間的相關性強弱有關，但目前能夠通過公開途徑免費獲取的相關輔助數據對于本研究區而言還存在分辨率相對較低等問題，因此，如何尋找與SOM相關性強的，容易獲取且廉價的其他輔助因子預測替代因子，或者開發適宜的輔助因子降尺度方法以實現輔助因子與SOM實測數據的尺度匹配，將是實現平原地區強烈人為作用影響下SOM空間預測的重要挑戰。