李夢佳，王 磊，劉洪斌*，武 偉

(1.西南大學資源環(huán)境學院，重慶 400716；2.西南大學計算機與信息科學學院，重慶 400715；3.重慶市數字農業(yè)重點實驗室，重慶 400716)

【研究意義】土壤有機質(Soil Organic Matter，SOM)是土壤的重要組成部分，與土壤肥力存在正相關關系，是反映土壤肥力的重要特征，與環(huán)境、大氣圈、生物圈等的可持續(xù)發(fā)展存在聯系[1-2]，受到不同因素的綜合影響，具有高度的空間差異性。利用不同模型預測有機質的空間分布，分析不同預測模型的差異，對土壤肥力研究、農業(yè)可持續(xù)發(fā)展有重要的意義[3]。【前人研究進展】近幾十年來，地統(tǒng)計學方法被普遍用來描述和預測土壤數據的空間分布，已有的研究表明地統(tǒng)計學是研究土壤特征及其不一致性空間分配、減少評估誤差的有效方法[4-6]。普通克里格因其操作簡便而成為最廣泛使用的地統(tǒng)計學方法之一，國內外研究者運用普通克里格預測了土壤有機碳、pH、有機質、堿解氮、有效磷、速效鉀及微量元素的空間分布，并分析了與其他預測模型預測結果的差異，認為普通克里格可以用于預測不同土壤數據的空間分布[7-11]。隨著信息技術的發(fā)展，分類與回歸樹、隨機森林等機器學習算法被應用到土壤數據空間分布預測中，國內外研究者運用機器學習模型預測了土壤水分、有機碳、pH、有機質及全氮的空間分布[12-16]，研究表明隨機森林等能較好地模擬環(huán)境變量與土壤數據變量的關系，能夠有效地解決土壤數據變量與環(huán)境變量之間的非線性問題，是現在被廣泛應用的方法。【本研究切入點】普通克里格之所以具有較好的預測效果，是因為它有能力使估計誤差變化最小化[17]。普通克里格方法基于變量的空間自相關性，依賴于所有樣點的空間關系確定，一般情況下，獲得有較好代表性的空間關系需要數量足夠大和空間覆蓋比較好的樣本集，因此普通克里格要求樣本數量較多、分布均勻及代表性好[12]?？死锝鸱ǖ闹饕秉c是當樣本量少時預測將變得不確定，并且無法模擬因海拔高度、氣候和植被等差異引起的土壤特性局部變化，這個缺點可以通過大量采樣和使用與土壤特性高度相關的環(huán)境變量來解決[17]。隨著信息技術的發(fā)展，機器學習預測模型被廣泛應用，但隨機森林等屬于黑箱技術，依據它們所獲得的結果難以直觀地了解土壤與環(huán)境變量之間的定量關系[18]，而且尋找與研究變量相關性強的、容易獲取且廉價的環(huán)境變量具有一定挑戰(zhàn)性[19]，此外研究變量可能受到人為活動的影響，但是人為因子影響不穩(wěn)定且不易獲得定量數據。由此可知，不同的預測模型各有優(yōu)缺點，本研究是要利用充分的采樣點比較不同模型的預測能力?！緮M解決的關鍵問題】本文結合地形、氣候、植被和成土母質等9個環(huán)境變量，利用分類與回歸樹、隨機森林、隨機森林殘差克里格和普通克里格4種預測模型對研究區(qū)土壤有機質含量進行空間分布預測制圖，得到研究區(qū)有機質含量空間分布特征，為土壤屬性空間變異研究提供依據，并比較不同預測模型的預測精度，分析不同預測模型繪制空間分布圖的差異，為高密度采樣區(qū)制圖方法的選取提供參考。

1 材料與方法

1.1 研究區(qū)概況

長壽區(qū)位于重慶市中部、長江三峽庫區(qū)上游，地理位置為106°49′～107°27′E，29°43′～30°12′N，總面積為1423.62 km2，研究區(qū)DEM及樣點分布見圖1。該地區(qū)屬四川盆地東部平行嶺谷褶皺低山丘陵區(qū)，為典型低山丘陵地帶，低山丘陵區(qū)的坡度較緩、起伏不大。亞熱帶濕潤季風氣候，常年平均氣溫17.68 ℃，常年降水量1162.7 mm，常年平均無霜期360 d。研究區(qū)土類包括水稻土、紫色土、潮土、黃壤和石灰土，其中水稻土占總耕地面積61.69 %，紫色土占總耕地面積35.05 %，潮土占總耕地面積0.25 %，黃壤和石灰土占總耕地面積3.02 %[20]。研究區(qū)成土母質包括二疊系和三疊系灰?guī)r，三疊系須家河組砂巖，下侏羅統(tǒng)自流井組粉砂巖，中侏羅統(tǒng)沙溪廟組粉砂巖和上侏羅統(tǒng)遂寧組砂泥巖(圖2)。

1.2 數據獲取

1.2.1 土壤有機質數據獲取 土壤有機質樣點數據從長壽區(qū)測土配方施肥項目中獲取，共有5162個樣本點，采樣時間集中在2009年，采樣深度為20 cm，采樣點空間分布情況如圖1所示。采集土壤樣本的同時，記錄采樣點經緯度、高程、土壤類型、地形部位等信息。采樣點土壤類型包括：紫色土(50.4 %)、水稻土(45.4 %)、黃壤(2.6 %)、潮土(1.3 %)、石灰土(0.3 %)。土壤樣品帶回實驗室后在室內自然風干，過篩后采用重鉻酸鉀測法測定土壤有機質含量。

1.2.2 環(huán)境變量數據獲取 土壤有機質的空間分布受到地形、氣候、植被和成土母質等因素的影響，本文選取9個環(huán)境變量：高程(Elevation,ELE)、相對坡度位置(Relative Slope Position,RSP)、溝谷深度(Valley Depth,VD)、地形濕度指數(Topographic Wetness Index,TWI)、垂直到溝谷距離(Vertical Distance to Channel Network,VDCN))、歸一化植被指數(Normalized Difference Vegetation Index,NDVI)、年均溫(Annual Average Temperature,ANTP)、年降水量(Annual Precipitation,ANPR)和成土母質(Strata,ST)。

(1)數字高程模型(DEM)數據來源于中國科學院計算機網絡信息中心地理空間數據云平臺(http://www.gscloud.cn)，空間分辨率為30 m。DEM數據通過SAGA GIS v.6.4提取4個地形因子：相對坡度位置(RSP)，溝谷深度(VD)，地形濕度指數(TWI)和垂直到溝谷距離(VDCN)。

(2)歸一化植被指數(NDVI)數據從http://ladsweb.nascom.nasa.gov下載獲得，空間分辨率為250 m。下載的遙感影像經過大氣校正、輻射校正和幾何校正，通過Savitzky-Golay濾波去除時間序列中的噪聲，并通過最鄰近法重采樣獲得分辨率30 m的NDVI數據。

(3)年均溫(ANTP)和年降水量(ANPR)數據來源于世界氣象數據庫(http://www.worldclim)，空間分辨率為1000 m。

(4)成土母質(ST)數據來源于全國第二次土壤普查，比例尺為1∶50 000。

1.3 預測模型

1.3.1 分類與回歸樹 分類與回歸樹(Classification and Regression Tree,CART)是一種監(jiān)督分類方法，它假設決策樹為二叉樹，利用訓練樣本來構造二叉樹并進行決策分類。CART可以用于分類或連續(xù)變量預測，可以處理非線性、非正態(tài)分布的數據，包容數據的缺失和錯誤，具有結構清楚、計算簡單和適用性廣等優(yōu)點[21]。構造CART樹分為樹生長和樹剪枝兩個步驟：首先基于總樣本生成一個多層次、多葉節(jié)點的決策樹，當所有葉節(jié)點中的樣本數為1或決策樹高度到達用戶設置的閾值時停止建樹，生成的決策樹足夠大，能充分反映數據之間的聯系；然后對決策樹進行剪枝，從產生的一系列子樹中選擇適當大小的樹，對未知數據進行準確預測[22-23]。

1.3.2 隨機森林 隨機森林(Random Forest,RF)是由多個決策樹形成的組合預測模型，是CART模型的擴展。RF克服了CART模型過度擬合的問題,能夠估計影響變量的重要性[24]，但對樣本數據集的大小較為敏感。RF的計算步驟為：從總樣本中有放回的隨機選出N個樣本構建回歸樹，不進行剪枝操作，余下的樣本作為袋外數據，袋外數據用來評價模型性能和自變量的重要性，通過投票或取平均值的方法，得到回歸樹最終預測結果。RF方法中的兩個重要參數是樹的數量和每個非葉子節(jié)點上可供選擇的變量數。

1.3.3 隨機森林殘差克里格 隨機森林殘差克里格(Random Forest with Residual Kriging,RFRK)是一種結合了RF與OK的混合地統(tǒng)計方法，它考慮了變量間的非線性關系及目標變量的空間自相關性，提高了預測精度。通過RF算法得到的預測值與觀測值之間存在殘差，如果這些殘差存在空間自相關性，可以將殘差當作隨機變量，通過克里格提高對未知點的預測精度。RFRK的計算步驟為：首先使用RF得到研究區(qū)有機質的預測值，計算RF預測殘差，然后對預測殘差進行OK插值，最后將RF預測值和OK插值誤差相加得到預測結果。

1.3.4 普通克里格 普通克里格(Ordinary Kriging,OK)是一種隨機性局部插值法，用已知點的樣本估計未知點的值，考慮變量的隨機性。OK著重于空間自相關因素，用擬合的半變異直接進行插值，通用方程式為：

(1)

式中，z為未知點的估計值，zx為x點的已知值，wx為x點的權重，n為用于估算的已知點數。權重不僅與估算點和已知點之間的半變異有關，還與已知點之間的半變異有關。

1.4 預測模型精度評價

根據模型構建的原則和相關研究，隨機選取總樣點的80 %作為訓練集，余下20 %作為驗證集，訓練集用于構建預測模型，驗證集用于評估模型的預測精度。統(tǒng)計指標包括平均絕對誤差(Mean Absolute Error,MAE)、均方根誤差(Root Mean Square Error,RMSE)和決定系數(Coefficient of Determination,R2)，當MAE和RMSE值越小、R2值越接近1時，模型精度越高。使用多個統(tǒng)計指標很難對預測模型精度進行排序，因此，本文使用能夠代表所有選定統(tǒng)計指標(MAE,RMSE,R2)的整體評價指標(Global Performance Indicator,GPI)對模型精度進行排序，結果顯示GPI值越大，預測模型總體效果越好[25-26]。計算公式分別為：

(2)

(3)

(4)

(5)

2 結果與分析

2.1 土壤有機質含量描述性統(tǒng)計

土壤有機質含量的描述性統(tǒng)計列于表1，可見總樣本均值為17.46 g·kg-1，取值范圍為5.10～78.80 g·kg-1，土壤有機質含量整體上屬于中等水平。土壤有機質含量的變異系數為33.73 %，屬于中等變異程度，表明研究區(qū)內土壤有機質具有一定的變異性，適合空間局部估計[27]。

表1 研究區(qū)采樣點土壤有機質含量描述性統(tǒng)計Table 1 Descriptive statistics of soil organic matter content in the study area

從5162個總樣本點中隨機選取4100個樣本點構建預測模型，剩余1062個樣本點用于評價不同預測模型的精度。普通克里格法要求輸入數據呈正態(tài)分布，但研究區(qū)土壤有機質含量的偏度大于1，不符合正態(tài)分布，因此對研究區(qū)SOM數據進行對數轉換，數據經轉換后符合要求。

2.2 土壤有機質含量與環(huán)境變量的相關性分析

2.2.1 土壤有機質含量與氣候、地形和植被的相關性分析 土壤有機質含量與年均溫(ANTP)、歸一化植被指數(NDVI)、相對坡度位置(RSP)、年降水量(ANPR)和垂直到溝谷距離(VDCN)呈極顯著負相關，相關系數分別為-0.25、-0.13、-0.10、-0.06、-0.06。土壤有機質含量與高程(ELE)、溝谷深度(VD)、地形濕度指數(TWI)呈極顯著正相關，相關系數分別為0.18、0.18、0.11。

通過SOM與ELE、VD、TWI、RSP、NDVI、ANTP、ANPR和VDCN的相關性分析，證實所選的環(huán)境變量對研究區(qū)SOM有顯著影響，應參與CART、RF和RFRK的模型構建(表2)。

表2 土壤有機質含量與環(huán)境變量的相關性Table 2 Pearson’s correlations between soil organic matter content and environmental variable

2.2.2 成土母質對土壤有機質含量的影響 土壤有機質含量在不同成土母質中存在顯著差異(P<0.05)，二疊系和三疊系灰?guī)r發(fā)育的土壤有機質含量最高，為27.34 g·kg-1，而下侏羅統(tǒng)自流井組粉砂巖發(fā)育的土壤有機質含量最低，為15.17 g·kg-1。說明成土母質對土壤有機質含量的空間分布有顯著影響，應參與CART、RF和RFRK的模型構建，這與Guo等學者的研究結果一致[28](表3)。

表3 成土母質對土壤有機質含量的影響Table 3 The influence of strata on soil organic matter content

2.3 土壤有機質含量的半方差分析

通過半方差分析可以得到變量的空間自相關性和最優(yōu)擬合模型，用于OK預測模型構建。研究中除了土壤有機質直接使用OK預測模型進行插值外，在RFRK預測模型插值過程中也用到OK預測模型，因此SOM(對數轉換)和RF預測殘差都需要進行半方差分析。

采用GS+軟件計算半方差，結果(表4)表明，SOM(對數轉換)和RF預測殘差的最優(yōu)擬合模型為指數模型，塊金效應分別為12 %和9 %。塊金效應表示空間相關性強弱，如果塊金效應小于25 %，表示強空間自相關性，說明變量具有很好的空間結構性[29-30]，研究中的塊金效應均小于25 %，屬于強空間自相關性，空間變異主要受結構性因素(例如成土母質、地形和氣候)的影響。

OK預測模型著重考慮空間自相關因素，用擬合的半變異模型直接進行插值，研究中的SOM(對數轉換)和RF預測殘差屬于強空間自相關性，滿足OK預測模型的使用要求，有條件取到較好的插值結果。

2.4 土壤有機質含量預測模型構建及模型精度

構建CART、RF預測模型時，通過網格搜索與交叉驗證確定預測模型的最優(yōu)參數組，最后確定CART預測模型中最小子節(jié)點、最小父節(jié)點和最大樹深度分別為10、20、20，RF預測模型中樹的數量和每個非葉子節(jié)點上可供選擇的變量數量分別為500個和3個，CART、RF預測模型在MatLab軟件中構建。OK預測模型在ArcGIS軟件中構建，需要的擬合模型及參數通過半方差分析得到(表4)。基于RF預測模型得到變量的預測值，計算預測殘差，然后對預測殘差進行OK插值，最后將RF預測值和OK插值誤差相加構建RFRK預測模型。

表4 土壤有機質含量的最優(yōu)半方差函數模型及參數Table 4 The semi-variance model and parameters for soil organic matter content

不同預測模型的精度見表5和圖3。CART、RF、RFRK、OK的GPI值分別為-2.25、-0.80、0、0.71，預測模型整體精度由低到高排序為CART

表5 不同土壤有機質含量預測模型精度Table 5 The accuracy of different prediction models on soil organic matter content

2.5 使用不同預測模型繪制空間分布圖

在使用不同預測模型繪制的空間分布圖(圖4)中，土壤有機質分布的總體趨勢一致，土壤有機質含量整體上屬于中等水平，高值主要集中在西部的明月山地區(qū)，而東部地區(qū)的有機質含量偏低，中部地區(qū)的有機質含量分布不均、高低錯落。模型精度最低的CART預測模型繪制的空間分布圖，出現邊界明顯的塊狀區(qū)域。OK預測模型與RF、RFRK預測模型相比，OK預測模型繪制的空間分布圖更平滑，RF和RFRK預測模型繪制的空間分布圖對局部細節(jié)刻畫更精細。

3 討 論

本研究中的地統(tǒng)計學預測模型OK整體精度高于機器學習預測模型CART、RF 和RFRK。研究區(qū)樣點密度較高、分布較均勻，半方差模型的隨機性會隨著采樣間隔的減小而減小，變量的空間相關性增強，在半方差分析中研究區(qū)土壤有機質含量呈現強烈的空間自相關性，空間自相關性極大地提高了模型的性能[31]，因此依賴于空間自相關性的OK預測模型精度高。Pouladi等[12]的研究結果同樣表明對于采樣密度較大的區(qū)域，克立格法可以不需要輔助變量而直接預測土壤有機質含量；而在觀測點數量有限的地區(qū)，結合輔助因子的方法表現更好[32]。機器學習模型通過模擬環(huán)境變量與土壤有機質含量的非線性關系，預測未知位置的土壤有機質含量，在土壤屬性受多種環(huán)境因素共同作用時，參與構建模型的指標越多，機器學習的預測能力可以得到越好的挖掘，而對于成土環(huán)境相對單一的區(qū)域，使用易于操作的地統(tǒng)計方法可以獲得理想的制圖效果[33]。謝恩澤等[19]在研究中提出：尋找與土壤有機質含量相關性強的、容易獲取且廉價的輔助因子具有一定挑戰(zhàn)性，這增加了機器學習預測模型的復雜性。分析表明采樣密度、空間自相關性和環(huán)境變量對研究結果有影響，所以當研究區(qū)已有土壤數據庫，且數據庫中土壤采樣點密度較高、分布較均勻，同時研究區(qū)位于坡度較緩、起伏不大的低山丘陵區(qū)，可以考慮采用操作簡單的OK預測模型繪制土壤有機質含量空間分布圖，有效降低工作量和工作難度。

機器學習預測模型整體精度由低到高排序為CART

通過對比分析4種不同預測模型繪制的空間分布圖，發(fā)現不同預測模型的土壤有機質空間分布總體趨勢一致，但局部細節(jié)存在差異，這與姜賽平等[35]研究結果一致。RF和RFRK預測模型繪制的空間分布圖不是光滑的連續(xù)曲面，而且CART預測模型繪制的空間分布圖出現邊界明顯的塊狀區(qū)域，這是因為在采用回歸樹方法繪制的土壤有機質空間分布圖過程中，每一個節(jié)點處土壤有機質含量值都是不連續(xù)的，因此導致生成的土壤有機質空間分布圖不是一個平滑的連續(xù)面，并且在節(jié)點較少的情況下，會導致與現實不符的土壤突變現象[36]?？死锔癫逯凳腔谝阎目臻g數據對未知空間數據進行估算，利用擬合模型將離散數據插值為連續(xù)的曲面數據，因此OK預測模型比機器學習模型繪制的空間分布圖更平滑。

4 結 論

在不同預測模型繪制的空間分布圖中，土壤有機質分布的總體趨勢一致，高值主要集中在西部的明月山地區(qū)，而東部地區(qū)的有機質含量偏低，中部地區(qū)的有機質含量分布不均、高低錯落。不同空間分布圖在局部細節(jié)中存在差異，RF和RFRK預測模型繪制的分布圖對局部細節(jié)刻畫更精細，OK預測模型繪制的空間分布圖更平滑。

預測模型整體精度由低到高排序為CART