基于隨機森林的縣域土壤有機碳密度及儲量估算

李海萍，杜佳琪，唐浩竣

（中國人民大學環境學院，北京 100872）

土壤中的有機碳庫相比于無機碳庫更加活躍且總量更大［1］，故大部分氣候變化研究都將有機碳庫作為研究對象，研究表明全球土壤有機碳總量為1.15×103～2.00×103Pg，約占陸地生態系統碳儲量的2/3［2-3］。土壤有機碳在空間上具有高度的變異性，對其進行精確評估對土壤肥力評價及土壤合理利用、全球碳匯及氣候變化研究具有重要意義。

早期的土壤有機碳估算基于少量土壤數據進行，Rubey［4］使用9 個土壤剖面的碳含量，估算出全球土壤有機碳儲量為710 Gt。Bohn［5-6］在Rubey的基礎上根據不同土壤類型和剖面數據，估算出全球土壤有機碳儲量分別為2946 和2200 Gt。Batjes［7］在前人基礎上又增加了厚度、容重、有機碳及礫石含量等參數，估算全球土壤有機碳儲量在1462～1548 Gt 之間。20 世紀80 年代起，有學者利用生命帶類型及其分布面積計算土壤有機碳儲量［2］。這種方法難以統計出全球植被類型及其面積，與土壤類型的對應關系也不夠準確，且未考慮土地利用等其他因素，難以應用于小尺度區域。20 世紀90 年代后，基于GIS的空間差值方法得到應用，Galbraith等［8］通過空間插值估算整個研究區的土壤有機碳儲量，李海萍等［9］將土壤類型法與GIS 方法相結合，在縣域尺度上估算了土壤有機碳儲量，該方法與土壤類型法類似，只是將結果精確直觀地通過GIS 展示出來。

基于土壤類型法和觀測數據的估算，發現該方法成本較高卻精度有限。許多學者采用遙感數據反演土壤有機碳，將高時空分辨率影像、反映生態系統碳循環動態變化的過程模型、實測土壤有機碳結合起來，從而提高了結果的空間分辨率［10］。也有學者通過建立遙感數據與土壤有機碳儲量相關模型進行估算［11-12］，不同學者在不同條件下所使用的模型具有明顯的針對性，鮮有縣域尺度下的土壤有機碳精確估算研究。

模型法能模擬和預測土壤碳儲量在自然界中的動態變化趨勢，也可解決尺度變換問題，因此能全面細致描述土壤有機碳的含量。隨著技術的快速發展，隨機森林模型也出現在土壤有機碳的估算研究中，Yang 等［13］使用增強回歸樹和隨機森林模型繪制了青藏高原東北部的土壤有機碳含量分布圖，Sreenivas 等［14］對比多種土壤有機碳儲量預測方法，通過隨機森林模型給出各個解釋變量重要性的排序。

云南省玉龍納西族自治縣（玉龍縣）有山地、盆地、河谷3 種地貌類型，海拔落差大，山區和半山區占縣域總面積的96%以上，地形復雜。因此以其作為研究對象，采用集合方式確定回歸樹，在縣域尺度上進行土壤有機碳儲量估算不僅能提高小尺度的估算精度及效率，也為大尺度估算提供精確的基礎數據。

1 材料與方法

1.1 研究區

地處云南省西北部的玉龍納西族自治縣，位于青藏高原與云貴高原交界的橫斷山脈中部，金沙江從西北向東北呈近90°的拐彎，形成了萬里長江的第一彎，縣域總面積6198.76 km2（圖1）。

玉龍縣屬典型山地高原，受新構造和河流切割影響，高山峽谷地貌特征明顯，海拔1300～4500 m，土壤垂直地帶性發育完全，依次為亞高山寒漠土（4200～4500 m）、亞高山草甸土（3500～4200 m）、暗針葉林土（3600～3800 m）、暗棕壤（3200～3600 m）、棕壤（2600～3200 m），黃棕壤廣泛分布于2500～2800 m 地勢低凹的濕潤山地，1300～2600 m 為紅壤，包含棕紅壤、黃紅壤、紅壤3 個亞類，非地帶性的紫色土、沼澤土、草甸土及石灰（巖）土等也有分布，見圖2。

1.2 數據

數據主要為土壤采樣數據、遙感影像及其派生數據和基于數字高程模型（DEM）提取的地形參數。

土壤采樣數據主要依據第二次土壤普查技術要求的采樣布點原則，于2009 年5 月至11 月進行實地采樣，使用棋盤法布點10～15 個，取樣后采用四分法將1 kg 土樣分開取樣，經實驗室分析得到的結果，對原始數據進行檢查后共得到686 個有效樣點，空間分布見圖3。

為保持數據的時間一致性，選用2009 年9 月Landsat7-Level2 影像，空間分辨率30 m。對因SLC傳感器故障導致的條帶在圖像處理軟件ENVI 中采用fillgap 插件去除。已有研究表明歸一化植被指數（NDVI）與土壤有機碳存在正相關關系［15］，因此，基于原始遙感影像，計算NDVI；纓帽變換是一種特殊的主成分變換，最早由Kauth 等［16］提出，通過線性變換、多維空間旋轉，不僅可提取植被和土壤信息，還可通過坐標變換將植被與土壤的光譜特征分離，故通過纓帽變換提取地表輻射亮度、植被綠度及土壤濕度等參數。

地形因素可在一定程度上反映土壤有機碳的狀態，不同高度因光照、溫度、水分等差異而使土壤微生物的活動也不相同，間接影響到土壤的腐殖化過程。將高程、坡度、坡向、曲率、地形濕度指數等作為土壤有機碳估算的輔助數據，能有效提高估算結果的精度。地形變量的提取主要基于數字高程模型（DEM），本研究采用SRTM（航天飛機雷達地形測繪使命）90 m DEM，分別提取坡度、坡向、曲率和地形濕度指數，并以45°間距將坡向分為8 個方向。

1.3 基于克里金差值的土壤有機碳儲量估算

土壤在空間上呈連續變化，因此，需要將采樣點數據空間化，采用克里金插值，該方法因考慮樣本的形狀、大小及其與預測點空間位置的關系而成為常用的線性無偏最優估計方法，公式為：

克里金插值考慮空間自相關，因此，先要判定數據的半變異函數并擬合相關參數，在ArcGIS 中將采樣點的經緯度坐標轉換為投影坐標并導出文本文件，采用GS+7.0 地統計分析軟件分別進行線性、球型、指數和高斯4 種常用理論模型的擬合，擬合參數見表1。

表1 半變異函數擬合模型參數

比較發現，指數模型殘差最小而決定系數最大，0.484的塊基比說明土壤有機碳密度的空間變異是隨機性因素與結構性因素共同引起的。故選擇指數模型進行插值試驗，當步長為10 時擬合效果最好。

根據插值結果并結合土壤類型，計算土壤有機碳儲量，公式如下：

因原始數據缺乏土壤容重屬性，故依據Song等［17］的研究結果進行估計，公式如下：

因采樣點均為耕地，各類型土壤的礫石含量較少，故采用耕地容重模型并忽略礫石含量，得到基于柵格的土壤有機碳密度。

1.4 隨機森林模型估算土壤有機碳密度

1.4.1 隨機森林模型

隨機森林是基于分類和回歸樹（CART）延伸出的機器學習新算法，由Breiman 等提出，在運算量沒有顯著提高的前提下提高了預測精度，且對多元共線性不敏感，結果對缺失數據和非平衡數據比較穩健，能很好地預測多達幾千個解釋變量的作用［18］，被譽為當前最好的數據挖掘算法之一［19-20］。

決策樹通過建立對象間屬性與值的映射關系來構建訓練數據的樹模型，既可用于分類，也可進行連續變量的預測。算法將所有數據作為一個根節點，從全部特征中挑選一個分割后生成若干子節點，采用節點-分支的樹結構進行決策，每個非葉子節點是一個判斷條件，每個葉子節點是結論，從根節點開始，判斷每一個子節點，若滿足停止分裂的條件則將該節點設為子節點，多次判斷后輸出節點占比最大的類別。

隨機森林對決策樹進行了優化，在變量（列）和數據（行）的使用上進行隨機化，生成多個分類樹，再匯總分類樹的結果。從解釋變量中隨機抽取n 個子樣本，針對每個樣本建立一棵分類樹，通過對多個樣本數據訓練、分類后進行預測，森林由多個決策樹組成，為避免樹間的相關性，采用套袋法獲取不同的訓練數據以增加其多樣性，再通過放回方式抽取數據集，過程中的數據完全隨機，有可能某一數據被多次使用，也會出現從未使用的情況，結合回歸策略輸出最終結果。

隨機森林的構建方式使其更加穩定，構建森林時使用數據集的最佳特征降低每棵樹的強度，對于任意劃分的特征G，對應的任意劃分點S 劃分成的數據集A1與A2，且A1與A2各自的均方差與均方差之和最小時所對應的特征值為劃分點，隨機選擇每個節點的特征，對比不同情況下的誤差以減少樹間的相關性泛化誤差，最終的回歸預測值為所有樹的預測值均值。

1.4.2 隨機森林建模

土壤有機碳生成機理復雜且影響因素較多，隨機森林又具有較強的高維數據處理能力，還能對各個影響因素的重要性進行解釋，適于進行土壤有機碳的定量研究。

模型最初主要通過R 語言中的隨機森林包、柵格包等實現［21］，2019 年ArcGIS Pro 在其2.3 版本中增加了這一功能模塊，既可對采樣點數據建模，也可對矢量或柵格數據建模，通過加入相應的柵格或矢量變量可直接獲得預測結果的空間分布，并對多個變量的重要性進行排序，比編程更便捷高效。

在ArcGIS Pro 環境下，采用“分析-工具-空間統計工具-空間關系建模-基于森林分類與回歸”功能，將土壤有機碳密度（SOCD）作為預測變量，分別導入坡度、坡向、曲率、DEM、地形濕度指數等地形變量以及NDVI、亮度、綠度、濕度等環境變量，作為模型的解釋變量。

2 結果與分析

2.1 表層土壤有機碳儲量估算結果

2.1.1 變量的描述性統計

對687 個采樣點0～20 cm 深度的土壤有機質含量、有機碳含量、有機碳密度、土壤容重等屬性，以及NDVI、亮度、綠度、濕度等環境參數和坡度、坡向、曲率、地形濕度指數等地形參數的描述性統計結果見表2。

表2 土壤有機碳及其相關變量描述性統計

變異系數顯示，表層土壤的有機質含量、有機碳含量、有機碳密度、綠度、坡度、曲率的離散程度較大，說明研究區植被分布差異大，地形復雜度高，NDVI、亮度、濕度、坡向、地形濕度指數及高程的離散程度較小。

2.1.2 表層土壤有機碳總儲量

分別以克里金插值與隨機森林模型所預測的土壤有機碳密度為屬性值，在ArcGIS 中統計柵格單元的值及其面積，兩者的乘積即為該柵格的土壤有機碳儲量，匯總后得到整個研究區表層土壤的有機碳總儲量，克里金插值結果為2.4×108t，隨機森林結果為1.7×108t。

2.2 土壤有機碳密度空間分布

2.2.1 克里金插值空間分布特征

采用指數模型及其參數進行克里金插值計算，得到0～20 cm 表層土壤的有機碳密度及其空間分布，見圖5。

圖5 顯示，土壤有機碳密度的高值區集中在西北部、中部山地以及東南部金沙江及其支流沿岸的平坦區域，最大值為66.5 kg/m2，東北部的玉龍雪山為低值區，最小值為7.82 kg/m2，中部的水域和城鎮也為低值區，因采樣點的空間集群特征而使插值結果不夠平滑，對空間分布的連續變化特征表達較差。

2.2.2 隨機森林空間分布特征

以土壤有機碳密度為預測變量，以NDVI、亮度、綠度、濕度、曲率、坡度、坡向、地形濕度指數為解釋變量，進行隨機森林回歸預測，對葉子數、樹數多次實驗后發現，當生成1000 棵樹時誤差趨于穩定，故將森林樹量設為1000 棵，葉子節點為5，樹深范圍0～25，平均為7，每棵樹可用的數據百分比為100%，隨機采樣變量數為3，輸出柵格與Landsat 一致，為30 m×30 m，得到土壤有機碳密度的空間分布如圖6。

結果顯示，土壤有機碳密度最高為79.37 kg/m2，最小為13.34 kg/m2，整體上西高東低，西部的高值區主要為河流下游及其沿岸。對照土壤類型分布，紅壤及棕壤的有機碳密度較高，此外，海拔適中地區的有機碳密度也較高，中部地形平坦區及東部邊界區則較低，0 值的空白區為金沙江河谷兩岸及麗江市城區，是河流與人類活動區的多個變量為空值所致。

隨機森林所提供的變量重要性排序能反映其對預測變量的解釋程度，根據去除該變量時模型所受的影響進行估算，結果見表3。

表3 隨機森林模型各變量的重要性排序

NDVI 占比最高，坡向最低，說明NDVI 對土壤有機碳密度的解釋作用最大，坡向則最小，其他變量的解釋程度基本一致。

2.3 精度驗證

2.3.1 普通克里金插值的精度驗證

將數據分成訓練集（2/3）和測試集（1/3），用測試數據集驗證克里金插值的精度表現，比較不同樣點的測量值與預測值，計算均方根誤差、標準平均值、標準均方根誤差、平均標準誤差等精度評價指標，結果見表4。

表4 普通克里金插值模型精度

精度驗證的參考標準一般為均方根誤差越小越好、標準平均值接近0，平均標準誤差接近于均方根誤差，平均標準誤差接近于1的估計為最優估計。本研究的標準平均值近于0，標準平均值很小，平均標準誤差接近于均方根誤差，但均方根誤差達20.77，說明克里金插值結果與最優估計還存在一定差距。主要由于實際采樣困難，難以達到插值所要求的樣點均勻分布條件。

2.3.2 隨機森林模型精度驗證

隨機森林算法每次抽取約2/3的樣本，套袋時第k 棵樹中未被選用的訓練樣本稱為“袋外”（OOB）子集。OOB 子集對于訓練集是不可見的，可用來對模型精度進行驗證，此處采用默認設置，即驗證數據為總數據的10%，驗證精度見表5。

表5 隨機森林精度驗證

決定系數R2是主要的性能指標，本研究的驗證結果大于0.5，說明模型的解釋能力強，且較為穩定，達66.7%。但均方根誤差較大，究其原因，因Landsat7 原始影像存在缺陷，雖進行了條帶修復，但條帶區域的誤差無法避免，因此用其提取其他解釋變量時必然存在誤差傳遞和積累，加之橫斷山區地形復雜，也使驗證結果不甚理想。

3 討論

空間插值主要根據有機質、容重、采樣深度等計算出的采樣點土壤有機碳含量進行。不同土壤的有機質轉換系數應有所不同，本研究由于數據的缺陷而采用了統一的轉換系數，故不能精確解釋每個采樣點的情況。此外，土壤容重主要通過環刀法測定，而本研究的采樣數據未進行容重測定，基于經驗公式的模擬結果精度有限。再者，樣點空間分布不均勻也導致了插值結果的空間分布不平滑，樣點稀疏區域比較粗糙。

隨機森林模型需要較多的解釋變量，本研究基于Landsat7 影像提取其他變量，雖進行了條帶修復，但仍然影響到所提取變量的精度。90 m DEM與30 m Landsat7 空間尺度的不一致也影響了模型結果的精度，基于土壤類型法計算出的有機碳密度本身就存在誤差，再經過解釋變量的回歸模擬必然會出現誤差傳遞和放大效應，增加了結果的不確定性。

兩種方法的差異性表現在隨機森林的高值區分布在西北、中部及東南部，低值區在東北及南部地區，克里金的低值區與隨機森林大體相同，但中部的高值區差別較明顯。

此外，本文僅對0～20 cm的表層土壤有機碳密度進行預測，未能對100 cm 內的土壤進行分層預測，對缺失的容重參數進行擬合也影響了結果的精度。

解釋變量的選取是影響隨機森林預測結果的重要因素，若能獲取更深層的土壤數據，并輔之以人類活動、氣候、土地利用、pH 值等更多解釋變量，就不會出現本研究因解釋變量值缺失所造成的人類活動區域結果為0的情況。

4 結論

基于GIS的普通克里金插值所需變量較少，能夠快速預測土壤有機碳密度及其空間分布，且預測結果與隨機森林結果相差不大，但其結果受采樣點數量及其空間分布的影響較大，用于橫斷山區的玉龍縣這種樣點少且分布不均勻的情況時結果較粗糙。隨機森林模型因所需變量多而效率較低，且算法復雜也使處理速度較慢，但預測結果的精度以及對小尺度區域的細節表現優于普通克里金插值，通過變量的相對重要性評估也可深入了解各變量對模型的影響，因此更適于對高維數據進行分析。