李海萍，杜佳琪，唐浩竣

（中國人民大學(xué)環(huán)境學(xué)院，北京 100872）

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)

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

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

1.2 數(shù)據(jù)

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

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

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

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

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

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

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

表1 半變異函數(shù)擬合模型參數(shù)

比較發(fā)現(xiàn)，指數(shù)模型殘差最小而決定系數(shù)最大，0.484的塊基比說明土壤有機(jī)碳密度的空間變異是隨機(jī)性因素與結(jié)構(gòu)性因素共同引起的。故選擇指數(shù)模型進(jìn)行插值試驗(yàn)，當(dāng)步長為10 時(shí)擬合效果最好。

根據(jù)插值結(jié)果并結(jié)合土壤類型，計(jì)算土壤有機(jī)碳儲量，公式如下：

因原始數(shù)據(jù)缺乏土壤容重屬性，故依據(jù)Song等［17］的研究結(jié)果進(jìn)行估計(jì)，公式如下：

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

1.4 隨機(jī)森林模型估算土壤有機(jī)碳密度

1.4.1 隨機(jī)森林模型

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

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

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

隨機(jī)森林的構(gòu)建方式使其更加穩(wěn)定，構(gòu)建森林時(shí)使用數(shù)據(jù)集的最佳特征降低每棵樹的強(qiáng)度，對于任意劃分的特征G，對應(yīng)的任意劃分點(diǎn)S 劃分成的數(shù)據(jù)集A1與A2，且A1與A2各自的均方差與均方差之和最小時(shí)所對應(yīng)的特征值為劃分點(diǎn)，隨機(jī)選擇每個(gè)節(jié)點(diǎn)的特征，對比不同情況下的誤差以減少樹間的相關(guān)性泛化誤差，最終的回歸預(yù)測值為所有樹的預(yù)測值均值。

1.4.2 隨機(jī)森林建模

土壤有機(jī)碳生成機(jī)理復(fù)雜且影響因素較多，隨機(jī)森林又具有較強(qiáng)的高維數(shù)據(jù)處理能力，還能對各個(gè)影響因素的重要性進(jìn)行解釋，適于進(jìn)行土壤有機(jī)碳的定量研究。

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

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

2 結(jié)果與分析

2.1 表層土壤有機(jī)碳儲量估算結(jié)果

2.1.1 變量的描述性統(tǒng)計(jì)

對687 個(gè)采樣點(diǎn)0～20 cm 深度的土壤有機(jī)質(zhì)含量、有機(jī)碳含量、有機(jī)碳密度、土壤容重等屬性，以及NDVI、亮度、綠度、濕度等環(huán)境參數(shù)和坡度、坡向、曲率、地形濕度指數(shù)等地形參數(shù)的描述性統(tǒng)計(jì)結(jié)果見表2。

表2 土壤有機(jī)碳及其相關(guān)變量描述性統(tǒng)計(jì)

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

2.1.2 表層土壤有機(jī)碳總儲量

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

2.2 土壤有機(jī)碳密度空間分布

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

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

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

2.2.2 隨機(jī)森林空間分布特征

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

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

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

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

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

2.3 精度驗(yàn)證

2.3.1 普通克里金插值的精度驗(yàn)證

將數(shù)據(jù)分成訓(xùn)練集（2/3）和測試集（1/3），用測試數(shù)據(jù)集驗(yàn)證克里金插值的精度表現(xiàn)，比較不同樣點(diǎn)的測量值與預(yù)測值，計(jì)算均方根誤差、標(biāo)準(zhǔn)平均值、標(biāo)準(zhǔn)均方根誤差、平均標(biāo)準(zhǔn)誤差等精度評價(jià)指標(biāo)，結(jié)果見表4。

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

精度驗(yàn)證的參考標(biāo)準(zhǔn)一般為均方根誤差越小越好、標(biāo)準(zhǔn)平均值接近0，平均標(biāo)準(zhǔn)誤差接近于均方根誤差，平均標(biāo)準(zhǔn)誤差接近于1的估計(jì)為最優(yōu)估計(jì)。本研究的標(biāo)準(zhǔn)平均值近于0，標(biāo)準(zhǔn)平均值很小，平均標(biāo)準(zhǔn)誤差接近于均方根誤差，但均方根誤差達(dá)20.77，說明克里金插值結(jié)果與最優(yōu)估計(jì)還存在一定差距。主要由于實(shí)際采樣困難，難以達(dá)到插值所要求的樣點(diǎn)均勻分布條件。

2.3.2 隨機(jī)森林模型精度驗(yàn)證

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

表5 隨機(jī)森林精度驗(yàn)證

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

3 討論

空間插值主要根據(jù)有機(jī)質(zhì)、容重、采樣深度等計(jì)算出的采樣點(diǎn)土壤有機(jī)碳含量進(jìn)行。不同土壤的有機(jī)質(zhì)轉(zhuǎn)換系數(shù)應(yīng)有所不同，本研究由于數(shù)據(jù)的缺陷而采用了統(tǒng)一的轉(zhuǎn)換系數(shù)，故不能精確解釋每個(gè)采樣點(diǎn)的情況。此外，土壤容重主要通過環(huán)刀法測定，而本研究的采樣數(shù)據(jù)未進(jìn)行容重測定，基于經(jīng)驗(yàn)公式的模擬結(jié)果精度有限。再者，樣點(diǎn)空間分布不均勻也導(dǎo)致了插值結(jié)果的空間分布不平滑，樣點(diǎn)稀疏區(qū)域比較粗糙。

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

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

此外，本文僅對0～20 cm的表層土壤有機(jī)碳密度進(jìn)行預(yù)測，未能對100 cm 內(nèi)的土壤進(jìn)行分層預(yù)測，對缺失的容重參數(shù)進(jìn)行擬合也影響了結(jié)果的精度。

解釋變量的選取是影響隨機(jī)森林預(yù)測結(jié)果的重要因素，若能獲取更深層的土壤數(shù)據(jù)，并輔之以人類活動(dòng)、氣候、土地利用、pH 值等更多解釋變量，就不會(huì)出現(xiàn)本研究因解釋變量值缺失所造成的人類活動(dòng)區(qū)域結(jié)果為0的情況。

4 結(jié)論

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