典型山西黃土高原區土壤有機質的空間異質性及空間插值預測

張小美 ，高春瑞 ，閆曉斌 ，楊 莎 ，喬星星 ，王 超 ，楊武德 ，Fahad Shafiq ，馮美臣 ，李廣信

（1.山西農業大學 農學院，山西 太谷 030801；2.拉合爾政府學院大學植物學系，巴基斯坦 旁遮普省 54000）

土壤作為地球陸地生態系統中最大的碳庫，是陸地生態系統的核心，是碳封存的關鍵，也是全球碳循環的重要通道[1]。據估計，全球1 m 深度的土壤有機碳儲量約為1550 Gt，大約是大氣碳庫的2 倍和陸地植物碳庫的3 倍[2]。土壤有機碳（Soil Organic Carbon，SOC）是土壤肥力形成的基礎，同時也是表征土壤質量和功能、減緩全球氣候變化的重要因素[3-5]。因此，明確區域SOC 含量的空間分布及其影響因素，對于研究SOC 的空間差異性分布及變化具有重要意義[6]。

獲取SOC 空間分布信息的傳統方法主要通過實地采樣、預處理和化學化驗等過程，存在耗時、費力、低效的缺點，難以滿足大樣本、大尺度、快速高效的現實需求。空間插值方法可以將SOC 含量從離散的點狀信息轉化為面狀連續信息，用來表征SOC 空間分布特征以及預測采樣點分布區域的土壤屬性空間信息，是較為常用的SOC 含量估算方法[7-8]。但是應用不同的插值方法得到的土壤有機碳的空間分布也會有所不同，王翔[9]研究發現，普通克里格插值可以很好地利用采樣點模擬空間聚集度；謝夢姣等[10]研究發現，徑向基函數插值神經網絡法可以避免普通克里格插值的平滑效應現象，提高土壤有機質空間預測精度；王海菁等[11]在宜黃縣城的洪水淹沒分析中表示，普通克里格插值相較于反距離權重法而言能使誤差的方差更小，更好地預測洪水的淹沒范圍。因此，探索預測精度更高的插值方法對實現土壤有機碳的高效預測具有重要意義。

SOC 的含量因易受到土壤形成因素（如環境因子、植被因子、水文因子、氣候因子和地形因子）的影響而具有很強的空間異質性。因此，當研究SOC 空間分布的區域尺度較大時，可能會存在影響因素多且不容易明確哪些因素影響比較大的問題，尤其是當氣象站點與采樣點距離較遠，極易造成較大的誤差，但研究小區域則能盡可能的減少這種誤差[12-13]。同時，基于單一因素的單一空間插值方法也會忽略部分環境因素帶來的影響，導致制圖精度降低[14-15]。為了能夠更好地捕捉環境變量在空間上對SOC 的影響，文雯等[16]、張宏帥等[17]將地形因子和氣象因子等輔助變量用于構建區域或國家尺度的SOC 空間插值模型，取得了較好的效果；楊順華等[18]研究發現，結合回歸克里格和空間位置的地理加權回歸克里格能更好地模擬SOC 空間分布格局；QIN 等[19]結合環境推理模型與反距離權重模型也提高了SOC 的作圖精度；江葉楓等[20]研究也指出，結合多個環境因素的回歸克里格模型精度要優于單一的普通克里格模型。因此，利用多模型結合對SOC 的時空變異性研究具有極大的應用潛力。

本研究以山西省忻州市忻府區農田SOC 為對象，通過引入溫度、降水、坡度、高程為環境輔助變量，在利用反距離權重法、徑向基函數插值、普通克里格插值、多元線性回歸、回歸克里格、回歸反距離權重法和回歸徑向基函數插值多種模型插值方法的基礎上，探討小尺度區域條件下SOC 與環境變量的關系，定量評價不同插值方法的表現，尋找最優的預測方法并進行空間制圖，以實現小區域尺度條件下農田SOC 的精確預測制圖，旨在為未來大尺度農田土壤碳庫估算提供一定的理論和技術參考。

1 材料和方法

1.1 研究區概況

忻州市忻府區位于山西省中北部（北緯38°6'～39°40'、東經110°53'～113°58'），屬季風型大陸性氣候，占地面積為1986.53 km2，整體地貌是被黃土廣泛覆蓋的山地型高原，忻府區地形西高東低，逐步傾斜，北、西、南三面環山，東部開闊平坦，是忻定盆地的主體部分。降水集中于7～9月，年均降雨量在462.5 mm 左右，無霜期平均167.1 d，年均溫在5.7 ℃，主要種植玉米、谷子、糜子、向日葵等作物。

1.2 數據來源

1.2.1 SOC 數據及預處理 本研究的SOC 數據來源于山西自然科技資源共享平臺（http://www.sxzrzy.ac.cn/），共獲取907 個土壤樣點數據。該數據是通過GPS 定點并實地考察，通過土種圖與土地利用現狀圖疊加進行評價，根據評價單元的個數及相應面積，在樣點總數的控制范圍內進行合理布局，并根據圖斑大小、種植制度、作物種類、產量水平等因素確定評價單元布點數量和點位，同時根據土壤類型、生態條件、種植面積等因素，選擇有代表性的土壤，采取定點取樣和隨機取樣相結合的方法進行布局及采樣，采樣深度為20 cm，利用外加熱-重鉻酸鉀氧化法測定獲取土壤有機質含量。為檢驗不同插值方法在研究區域的預測精度，運用Arc-GIS10.2 地統計工具下的Subset 將907 個實測樣本按照3∶1 的比例分為預測集和驗證集。

1.2.2 環境變量數據及預處理 選擇氣象因子（降水量和溫度）和地形因子（高程和坡度）作為輔助SOC 預測的環境變量（圖1），其中，氣象數據來自中國氣象數據網（https://data.cma.cn/），高程及坡度數據來源于山西自然科技資源共享平臺（http://www.sxzrzy.ac.cn/），將氣象數據進行整理和預處理，用Excel 計算出年均溫和年均降水量，利用Arc-GIS 10.2 提取出各采樣點的溫度和降水量。

圖1 環境變量空間分布Fig.1 Spatial distribution map of environmental variables

1.3 研究方法

本研究用到了非地統計、地統計、互結合三大類方法，其中，非地統計主要包括徑向基函數、反距離權重插值、多元線性回歸；地統計為普通克里格插值；互結合是基于非地統計和地統計的拓展方法，主要包括回歸克里格法、回歸反距離權重法和回歸徑向基函數插值[19]。本研究中用到的7 種插值方法都可以通過ArcGIS 10.2 實現。

1.3.1 反距離權重法（Inverse Distance Weighting Method，IDW） 其是對插值點與樣本點之間的距離進行加權平均，基于相近相似的原理距離越短，權重越大[16]。

1.3.2 徑向基函數插值（Radial Basis Function Method，RBF） 其是某種沿徑向對稱的標量函數，通常定義為樣本到數據中心之間徑向距離（通常是歐氏距離）的單調函數[21]。

1.3.3 普通克里格插值（Ordinary Kriging，OK）其是依靠所選采樣點數據在空間上加以插值，其插值效果較依賴于采樣點的個數、密度，以及其數據的準確性[22]。

1.3.4 多元線性回歸（Multiple Linear Regression，MLR） 其是對2 個或者多個變量進行回歸分析，本研究用預測集在SPSS 中進行逐多元逐步線性回歸，可以定量描述環境變量對SOC 的相關關系，并且避免了多重共線性[23]。

1.3.5 回歸克里格（Regression Kriging，RK） 其是MLR 與OK 的結合，首先需要通過MLR 得到殘差值，再用OK 對殘差進行插值，并最終將MLR 預測的趨勢項與OK 的插值結果相加，得到RK結果[22]。

1.3.6 回歸反距離權重法（Regression Inverse Distance Weighting Method，MIDW） 其是MLR 與IDW 的結合，通過MLR 模型得到殘差值，再用IDW 對殘差進行插值，最后將MLR 預測的趨勢項與IDW 的插值結果相加，從而得到MIDW 結果。

1.3.7 回歸徑向基函數插值（Regression Radial Basis Function Method，MRBF） 其是MLR 與RBF的結合，通過MLR 模型得到殘差值，再用RBF 對殘差進行插值，最后將MLR 預測的趨勢項與RBF的插值結果相加，從而得到MRBF 結果。

1.4 檢驗方法

為檢驗7 種插值方法的預測精度，本研究通過預測集中的樣點進行插值預測得到結果圖，再根據驗證集的位置提取出驗證集的預測值，對驗證集的實測值和預測值進行相關性分析，通過計算得到預測值和實測值的Pearson 相關系數（r）、平均絕對誤差（MAE）、均方根誤差（RMSE）作為檢驗插值精度的標準。MRE 和RMSE 值越小，說明結果越精確。

2 結果與分析

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

分別對總樣本（n=907）、預測集（n=680）和驗證集（n=227）的SOC 含量進行描述性統計分析，結果如圖2 所示。

圖2 土壤有機碳含量基本描述特征Fig.2 Basic description characteristics of soil organic carbon content

從圖2 可以看出，總樣本SOC 的平均值為15.05 g/kg，最大值為43.4 g/kg，最小值為1.4 g/kg；預測集的平均值為14.66 g/kg，最大值為43.4 g/kg，最小值為1.4 g/kg；驗證集的平均值為14.66 g/kg，最大值為36.2 g/kg，最小值為3.9 g/kg。總樣本變異系數為30%，即SOC 變化幅度較大，有較強的空間異質性。預測集的全距為42，驗證集的全距為32.3，預測集的全距涵蓋了驗證集的全距，從統計學的角度可知，該數據劃分合理。由圖2 可知，該數據呈近似正偏態分布特征，表明數據符合進一步分析的需要。

2.2 SOC 與環境變量的相關關系

利用所提取的坡度（S）、高程（E）、溫度（T）和降水量（P）與SOC 含量進行相關性分析，結果如圖3 所示，在地形因子中，高程（-0.255**）和坡度（-0.214**）與SOC 含量呈極顯著負相關（P＜0.01）；在氣候因子中，降水（-0.085*）和溫度（-0.246**）與SOC 呈顯著（P＜0.05）或極顯著負相關（P＜0.01）。結果表明，溫度與SOC 相關性最高，降水與SOC 相關性最弱。

圖3 環境變量與土壤有機碳的相關關系Fig.3 Correlation between environmental variables and SOC

2.3 7 種插值模型的構建及預測

將環境因素與SOC 數據進行多元線性回歸分析，得到多元線性回歸模型系數的基本統計特征，結果如表1 所示，MLR 模型的非標準化系數為高程、年均降水量、年均溫和坡度，4 個環境變量的方差膨脹因子均小于7.5，則說明這4 個環境因子相互不存在共線性。此外，通過t值和P值可以看出，環境變量與SOC 之間存在顯著相關關系（P＜0.05）；從標準系數看，溫度對SOC 影響最大。MLR 模型的R2=0.484，即環境變量可以解釋SOC 48.4%的變化。以此為基礎，構建RK、MIDW 和MRBF 模型，利用MLR 模型得到殘差值，分別將殘差值進行OK、IDW、RBF 插值，利用柵格計算器加上MLR圖層的有機碳含量值，最終得到RK、MIDW、MRBF的插值模型。

表1 多元線性回歸模型系數的基本統計特征Tab.1 Basic statistical characteristics of the coefficients of the multiple linear regression models

圖4 分別為MLR、RK、IDW、RBF、MIDW、MRBF 得到的SOC 空間插值預測情況，從整體來看，7 種插值模型所反映出SOC 含量的空間分布沒有明顯差異，但在采樣點所在區域附近還是有所不同，其中，IDW、RBF、MIDW、MRBF 相對來說所得到的預測圖不夠平滑，呈現“牛眼”現象，能夠很明顯地看出采樣點的位置；對比這4 種模型的SOC 空間預測圖可以看出，結合MLR 模型的MIDW、MRBF 模型相對于另外2 種更加平滑。

圖4 基于7 種模型的土壤有機碳空間分布預測Fig.4 Soil organic carbon spatial distribution prediction map based on 7 prediction models

2.4 插值結果的驗證與分析

從表2 可以看出，除OK 空間插值方法外，其余6 種方法的預測值與驗證值之間具有極顯著相關性（P＜0.01），相關系數大小為r（MRBF）＞r（RBF）＞r（IDW）＞r（RK）＞r（MIDW）＞r（MLR）；再對比7 種模型的MAE 和RSME 值可知，其中，MLR 的RMSE 最小，小于其他6 種模型，MRBF 的MAE 最小，因此，MRBF 模型預測精度最高，在對土壤有機碳進行空間插值時可以將MLR 模型結合地統計模型來預測環境變量下土壤有機碳的變化。

表2 7 種模型的精度評價指標比較Tab.2 Comparison of precision evaluation of seven models

3 討論

3.1 環境變量與SOC 之間的關系

本研究選用氣候、地形作為環境變量應用于模型，其中，氣候變化影響是SOC 濃度的重要因素[24]；溫度會影響植物的光合作用速率、凋落物分解速率及土壤有機質礦化速率，對土壤碳循環有著顯著影響[25]。本研究通過線性分析發現，溫度與SOC 呈顯著負相關，這可能是因為較低的溫度更有助于SOC的積累。羅梅等[13]研究基于環境變量的SOC 空間分布特征時通過相關性分析發現，溫度與SOC 呈負相關，與本研究結果一致。周莉等[26]研究指出，氣候條件會通過制約植被類型、影響植被生產力進而影響輸入土壤有機碳的輸入量；同時還能通過影響土壤水分和溫度等條件，進而影響微生物對有機碳的分解和轉化。降水通過影響陸地土壤碳密度和土壤的通氣性影響土壤持有的有機碳量[27]，當土壤水分充足時，則透氣性差，原有機碳不易礦化，外源有機殘體在水分作用下更易腐爛降解為有機質，從而提高SOC 含量；土壤水分缺失時，孔隙大，則會加速有機碳分解，不利于SOC 的積累。山西省屬于黃土高原，地質干旱，降雨量小，因此，降雨量與土壤有機碳的含量呈負相關[28]，與本研究結果一致。

此外，地形因子能調節熱量和物質的再分配，對SOC 分布起重要作用。其中，坡度是影響土壤侵蝕的重要因素，受到土壤侵蝕的影響，不同坡度下SOC 含量存在較大差異[27]。本研究中，坡度與土壤有機碳呈顯著負相關。土壤屬性的空間變異受多種因子的控制，而地形作為五大成土因素中的重要因素，能調節熱量和物質的再分配，對SOC 分布起重要作用[28]。李龍等[27]研究指出，坡度是影響農地固碳潛力的主要因素，對于不同的坡度單元，坡度大于25°的耕地退耕后的固碳效應最高，而坡度小于15°的效應最低。本研究中，坡度最大值為50°，大于25°的約占樣本數的1/30，因此，就坡度而言該地區的固碳潛力有待大幅度提升。本研究區域的高程與SOC 含量呈顯著負相關，高程自西向東呈現降低趨勢。高海拔地區通過影響溫度和降雨量制約土壤微生物的活性，進而影響土壤有機碳的分解速率[14]。高程條件通過影響區域的水熱條件來影響SOC 含量的累積[29]。但在本研究中，SOC含量自西向東逐漸升高，主要原因是因為該地區屬于黃土高原地區，較為干旱，溫度較低。在小尺度上氣候、植被以及地形因素在時間和空間上基本一致，沒有太大變化，樣點雖然很多，但是人為因素方面的影響仍然是不可準確估量的。因此，在大尺度上SOC 的空間分布仍然存在不確定性。

3.2 模型預測精度差異性對比

本研究使用MLR、RK、MRBF 和MIDW 模型結合地形和氣候因子預測SOC 含量的空間分布，對比沒有考慮地形和氣候影響的OK、RBF 和IDW模型可知，MRBF 模型精度最高。RK、MRBF 和MIDW 模型包含權重函數的MLR 變體，考慮了SOC 與環境變量之間的關系，MLR 是通過對預測集進行統計分析建模，預測精度主要受到數據量和數據精度的影響，但是OK、RBF、IDW 直接利用預測集進行插值，只考慮了樣本點自身的結構特征，忽略了環境變量的影響，因此，為了得到更精確的小尺度土壤有機碳的空間預測，應充分考慮土壤形成因素帶來的各種影響，MLR 可以得出多個環境變量對SOC 的影響，將MLR 與OK、RBF、IDW 結合得到的RK、MRBF、MIDW 模型即能更綜合地考慮到SOC 與環境變量的關系[17,30]，并且RK、MRBF、MIDW 模型精度相較于OK、RBF、IDW 也有所提升，江葉楓等[19]研究也指出，考慮輔助變量的模型會提高模型精度，更加證實了互結合模型在土壤有機碳空間預測方面的優勢，但是RK、MRBF、MIDW插值相對OK、RBF、IDW 等插值方法，步驟繁多，參數設置要求較高，耗時耗力，在研究目的許可的情況下，當SOC 數據樣點密度足夠大且不需要結合環境變量進行預測時，可以直接考慮簡單靈活的OK、RBF 或者IDW 插值，預測精度相對較高[31-32]。本研究表明，互結合模型在土壤有機碳空間預測方面有很好的應用潛力；隨著信息技術的不斷發展，這種互相結合的模型也將成為一種新的趨勢，本研究采用的是MLR 與其他模型的相互結合，未來可以考慮更多算法模型的相互結合，更好地達到預測效果；此外，本研究只考慮了氣候和地形因子對土壤有機碳的影響，考慮到土壤形成過程受到多種因素的影響，未來可以考慮加入更多的影響因子，更好地完善土壤有機碳的空間預測；如果是需要構建生態模型或是圖形疊加等空間分析，則需要考慮更多因素在空間上的互相影響，選擇插值精度更高、在空間上表現最好的空間插值方法進行預測[20]。

4 結論

本研究以高程、坡度、降水量、溫度為輔助變量，比較OK、RBF、IDW、MLR、RK、MRBF、MIDW等7 種模型的建模結果，研究結果表明，SOC 與高程、坡度、降水量以及溫度呈極顯著負相關關系；其中，溫度和高程與SOC 相關性最高，降水與SOC 相關性最弱；利用插值模型進行SOC 空間分布預測制圖，發現MRBF 模型具有較好的擬合精度，且預測結果優于其余6 種模型。該結果為促進土壤固碳提供了科學依據，并對理解小尺度SOC 的空間異質性及精確制圖提供了一定的理論和實踐參考。