基于多變量與RF算法的耕地土壤有機(jī)碳空間預(yù)測研究*——以福建亞熱帶復(fù)雜地貌區(qū)為例

袁玉琦，陳瀚閱，張黎明，任必武，邢世和，童珺玥

（福建農(nóng)林大學(xué)資源與環(huán)境學(xué)院，土壤生態(tài)系統(tǒng)健康與調(diào)控福建省高校重點(diǎn)實(shí)驗(yàn)室，福州 350002）

耕地土壤有機(jī)碳（Soil Organic Carbon，SOC）是衡量土壤質(zhì)量和表征土壤肥力的重要指標(biāo)，也是農(nóng)業(yè)溫室氣體減排潛力的主要來源。掌握土壤有機(jī)碳的含量及空間分布，對提升耕地質(zhì)量、評估土壤健康和緩解全球氣候變化具有重要意義。實(shí)地采樣和分析是獲取土壤屬性空間分布格局的傳統(tǒng)方法，但因費(fèi)時費(fèi)力、成本高、區(qū)域通達(dá)性差限制了其在地區(qū)、國家乃至全球尺度的應(yīng)用[1]。隨著 3S技術(shù)的發(fā)展，越來越多學(xué)者通過建立輔助環(huán)境變量與土壤屬性之間的預(yù)測模型來實(shí)現(xiàn)區(qū)域土壤屬性空間制圖[2]。

目前，用于數(shù)字土壤制圖的預(yù)測模型主要包括逐步線性回歸模型[3]、遞歸偏最小二乘算法[4]及隨機(jī)森林算法（Random Forest，RF）[5]等。其中RF算法是一種機(jī)器學(xué)習(xí)算法，它由分類和回歸樹（Classification and Regression Tree，CART）模型構(gòu)成，因其在處理多元非線性數(shù)據(jù)方面的優(yōu)勢，越來越多地被用于土壤屬性空間預(yù)測[6]。相繼有學(xué)者將RF模型用于大區(qū)域有機(jī)碳含量及儲量[7-9]的空間預(yù)測，并取得了較好的預(yù)測精度（r值在0.7～0.95范圍內(nèi)）。但RF算法在土壤有機(jī)碳空間制圖的應(yīng)用起步較晚，其在不同區(qū)域的預(yù)測精度和可行性方面仍需進(jìn)行深入研究，其中需要關(guān)注的一個研究重點(diǎn)即是環(huán)境變量的篩選[10-11]。

在基于 RF算法的土壤有機(jī)碳空間預(yù)測中，最為常用的輔助環(huán)境因子包括數(shù)字高程模型（Digital Elevation Models，DEM）和遙感植被指數(shù)[7，11]等遙感變量。土壤有機(jī)碳含量高低與地表植被覆蓋狀況和地形地貌密切相關(guān)[12]，而遙感變量因在地形和植被空間信息表達(dá)上具備獨(dú)特優(yōu)勢，且較易獲取，常被用于數(shù)字土壤制圖。但土壤理化性質(zhì)是多方面成土因素綜合作用的結(jié)果，單純利用遙感變量預(yù)測土壤有機(jī)碳有所欠缺。部分學(xué)者[8，13]在此基礎(chǔ)上增加了相對較易獲取的氣候因子（如溫度、降水）和土壤類型、成土母質(zhì)等土壤定性屬性。然而，與土壤有機(jī)碳機(jī)理上互為相關(guān)的土壤定量屬性數(shù)據(jù)因獲取難度較大，以往較少被用于土壤有機(jī)碳空間預(yù)測，但隨著測定技術(shù)的發(fā)展，該類變量獲取難度有所降低，日漸被部分學(xué)者關(guān)注。Liu等[14]在植被要素和氣候因子等較易獲取的環(huán)境變量基礎(chǔ)上增加了容重、土層厚度等土壤屬性，基于RF算法模擬中國北方草地土壤有機(jī)碳密度（R2達(dá)到0.73）。Were等[9]在對比三種機(jī)器學(xué)習(xí)算法（包括RF算法）用于南非地區(qū)SOC儲量預(yù)測的研究中，創(chuàng)新性地選擇Mg、Ca、P、全氮（Total Nitrogen，TN）和pH等土壤定量屬性聯(lián)合遙感數(shù)據(jù)和地形及衍生因子作為模型驅(qū)動因子，結(jié)果顯示土壤屬性中的 TN變量在所有三種模型的變量貢獻(xiàn)率中占絕對主導(dǎo)地位。

由此可見，基于RF模型預(yù)測SOC所用的環(huán)境變量隨科學(xué)測定技術(shù)的發(fā)展有所創(chuàng)新，但國內(nèi)外研究少有基于不同組合的環(huán)境變量對預(yù)測精度影響的評估，尤其是獲取難度較大的土壤定量屬性，作為變量加入對模型精度的貢獻(xiàn)如何尚不十分明確。此外，針對國內(nèi)基于RF算法的SOC預(yù)測多應(yīng)用于地貌類型較單一的小尺度區(qū)域[13，15]，而大范圍復(fù)雜地貌區(qū)因地形、氣候、植被、土壤肥力狀況等因素在水平和垂直方向上存在顯著分異[16]，導(dǎo)致SOC的空間分布規(guī)律及其與環(huán)境變量間的關(guān)系更為復(fù)雜，不同變量對SOC預(yù)測精度的貢獻(xiàn)率仍需進(jìn)一步評價。RF模型正因?yàn)樵谀M變量間復(fù)雜非線性關(guān)系和高階相互作用方面的突出優(yōu)勢[17]，更適用于本文亞熱帶復(fù)雜地貌區(qū)土壤有機(jī)碳含量的空間預(yù)測。

目前，我國已經(jīng)完成兩次全國性的土壤普查、農(nóng)業(yè)農(nóng)村部測土配方施肥重大國家計(jì)劃，也正在進(jìn)行農(nóng)業(yè)長期定位實(shí)驗(yàn)和第三次土地調(diào)查等，積累了大量的土壤樣點(diǎn)及屬性數(shù)據(jù)，這些較難獲取的土壤定量屬性數(shù)據(jù)可為提高模型預(yù)測精度提供極為便利的研究條件?；诖?，本文以典型的亞熱帶復(fù)雜地貌區(qū)——閩東南地區(qū)為研究區(qū)域，采用隨機(jī)森林模型，選取易獲取的遙感變量和氣候因子及不易獲取的土壤屬性作為模型輸入，基于大量土壤有機(jī)碳實(shí)測樣點(diǎn)分別訓(xùn)練兩種不同變量組合數(shù)據(jù)集驅(qū)動的RF模型，并進(jìn)行有機(jī)碳含量空間分布格局的預(yù)測和對比，以期為區(qū)域地形地貌復(fù)雜區(qū)有機(jī)碳儲量的準(zhǔn)確估算提供理論基礎(chǔ)。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于福建省東南部，簡稱閩東南地區(qū)（23°32′～26°04′ N，116°53′～119°91′ E），由漳州市、廈門市、莆田市、泉州市（除永春縣、德化縣和金門縣）和福州市（僅長樂區(qū)、福清市和平潭綜合試驗(yàn)區(qū)）組成（圖1），區(qū)域總面積為29 373 km2，其中耕地面積約4 384 km2。該區(qū)位于閩中大山帶戴云山-博平嶺段東南側(cè)，地貌復(fù)雜、地形起伏，海拔最高1 152 m，最低0 m，從內(nèi)陸向沿海遞減，其中山地丘陵（海拔>200 m）面積高達(dá)44.33%。受亞熱帶季風(fēng)氣候與地形因素綜合作用，年平均氣溫為 20.0℃，年平均降雨量為 1 526 mm，主要集中于 3—8月，占全年降雨量80%以上。從土地利用方式來看，可以用作農(nóng)業(yè)用途的土地占比?。ǜ貎H為14.93%），且受海拔、地形等要素影響，80.15%的耕地位于福建省最主要的農(nóng)業(yè)高產(chǎn)區(qū)——沿海平原臺地（海拔≤200 m）[18]。該區(qū)土層深厚、肥力高，農(nóng)業(yè)生態(tài)氣候條件好，耕地利用強(qiáng)度大，土壤屬性空間變異強(qiáng)烈。全區(qū)耕作田塊分布較為零散，耕地面積小而破碎，土壤類型以水稻土為主，占耕地總面積的 70.97%，其次為赤紅壤，占耕地總面積的19.91%，而紅壤、濱海鹽土、潮土、風(fēng)砂土、黃壤、石灰土和紫色土的分布面積較小，合計(jì)面積僅占總面積10%左右?？傊孛驳奶厥庑耘c復(fù)雜性導(dǎo)致氣候、植被、土壤等空間分異顯著，造成不同地區(qū)、不同海拔高程的耕地自然條件差異明顯[16]。

1.2 數(shù)據(jù)來源

本研究采用的土壤有機(jī)碳實(shí)測數(shù)據(jù)來源于國家農(nóng)業(yè)農(nóng)村部2017年末測土配方施肥調(diào)查樣點(diǎn)數(shù)據(jù)，共計(jì)1 257個，包括調(diào)查樣點(diǎn)地理坐標(biāo)及土壤有機(jī)質(zhì)（Soil Organic Matter，SOM）、有效磷（AP）、速效鉀（AK）、水解性氮（HN）、交換性鎂（Ex-Mg）和 pH等土壤理化屬性（圖1）。每個樣點(diǎn)均按照密度控制、代表性、均勻性、優(yōu)先性及適當(dāng)性原則設(shè)置，采樣深度為 0～20 cm。在所選地塊均勻隨機(jī)采集10～15個耕層土樣，充分混合后采用四分法留取1.5 kg土樣裝袋，自然風(fēng)干過篩后備用。土壤有機(jī)質(zhì)含量采用重鉻酸鉀氧化—外加熱測定，有效磷采用碳酸氫鈉浸提—鉬銻抗比色法測定，速效鉀采用乙酸銨提取—火焰光度計(jì)法測定，水解性氮采用堿解擴(kuò)散法測定，交換性鎂采用乙酸銨浸提—原子吸收分光光度法測定，pH采用酸度計(jì)法測定。

選取樣點(diǎn)中測定的土壤理化屬性，聯(lián)合遙感變量和氣候變量作為SOC空間預(yù)測的環(huán)境變量，具體構(gòu)成見表1。遙感變量中的植被指數(shù)和地形因子分別基于Landsat8 OLI影像和ASTER GDEM高程影像提取得到，數(shù)據(jù)來源于美國地質(zhì)勘探局（United States Geological Survey，https://www.usgs.gov/）和地理空間數(shù)據(jù)云網(wǎng)站（Geospatial Data Cloud，http://www.gscloud.cn/）。影像空間分辨率均為30 m，投影坐標(biāo)系為WGS_1984_UTM_Zone_50N。為使植被指數(shù)能反映采樣時地表真實(shí)狀況，影像選取月份與實(shí)際調(diào)查樣點(diǎn)獲取時間一致，且天氣晴朗，基本無云層覆蓋。氣象因子基于2017年福建全省22個標(biāo)準(zhǔn)氣象站點(diǎn)的地面氣候資料日值數(shù)據(jù)集插值而來（具體站點(diǎn)位置見圖1），數(shù)據(jù)來源于國家氣象中心網(wǎng)站（National Meteorological Information Center，http://data.cma.cn/）。

1.3 環(huán)境變量的獲取、組合與篩選

為了探索多源協(xié)同變量解釋 SOC空間變異的可能性及土壤屬性的加入對提高模型預(yù)測精度的貢獻(xiàn)程度，本研究將表1中16個環(huán)境變量根據(jù)獲取的難易程度分成兩種不同的組合：（1）僅基于易獲取的遙感變量（Remote Sensing Variables，RS）和氣象因子（Climate Factor，CF）訓(xùn)練的模型Training model based on simple-to-obtain variables（RF-S）；（2）基于遙感變量、氣象因子和土壤屬性（Soil Attribute，SA）所有變量訓(xùn)練的預(yù)測模型 Training model based on all variables（RF-A）。

表1 環(huán)境變量的構(gòu)成Table 1 The composition of environmental variables

影響土壤有機(jī)碳含量的環(huán)境變量眾多，模型訓(xùn)練前需利用 RF算法預(yù)測所產(chǎn)生的袋外誤差的大小對部分變量進(jìn)行剔除[10]，即依據(jù)逐次剔除某一變量后RF模型袋外得分（Out-of-bag Score，OOB Score）的增減判斷該變量是否保留，OOB Score值增加則變量剔除，反之保留[11]。

1.4 RF模型構(gòu)建和驗(yàn)證

以各模型篩選后的變量數(shù)據(jù)集為輸入，利用RF模型進(jìn)行有機(jī)碳含量的回歸預(yù)測。RF模型是建立在決策樹基礎(chǔ)上的一種集成學(xué)習(xí)方法，通過多次bootstrap抽樣獲取多個隨機(jī)樣本，并通過這些樣本子集分別構(gòu)建相應(yīng)決策樹，從而構(gòu)建隨機(jī)森林[10]。當(dāng)模型用于回歸預(yù)測時，取所有決策樹預(yù)測結(jié)果的均值作為最終的預(yù)測結(jié)果[19]。模型的運(yùn)算過程中需設(shè)定兩個關(guān)鍵參數(shù)：n_estimators和max_depth，其中 n_estimators為決策樹的數(shù)量，即使用 bootstrap重抽樣的次數(shù)，max_depth為決策樹的最大深度。依據(jù)預(yù)測過程中產(chǎn)生的OOB Score的大小，RF-S和RF-A模型設(shè)定的（n_estimators，max_depth）分別為（1 400，9）和（1 400，12）。

為衡量兩種不同變量組合下模型的表現(xiàn)，將所有訓(xùn)練樣本點(diǎn)作為模型驗(yàn)證點(diǎn)進(jìn)行驗(yàn)證。此外，使用不同抽樣百分比（20%、30%、40%、50%、60%和70%）的驗(yàn)證數(shù)據(jù)集獨(dú)立評估模型性能（例如，80%樣點(diǎn)訓(xùn)練，余下20%樣點(diǎn)進(jìn)行驗(yàn)證，以此類推），并與常用的地統(tǒng)計(jì)插值方法——普通克里格插值模型（Ordinary Kriging，OK）結(jié)果進(jìn)行對比。使用平均絕對誤差（MAE）、均方根誤差（RMSE）、相關(guān)系數(shù)（r）和變異系數(shù)（CV）評估模型預(yù)測的絕對誤差表現(xiàn)，使用相對誤差（Relative Error，RE）和相對均方根誤差（Relative RMSE，RRMSE）定量化模型的準(zhǔn)確性程度，數(shù)值越小，模型準(zhǔn)確性越佳。

1.5 數(shù)據(jù)處理方法

數(shù)據(jù)處理主要為環(huán)境變量提取、RF模型構(gòu)建和驗(yàn)證，以及耕地土壤有機(jī)碳空間分布圖的生成3個方面，具體如下：

環(huán)境變量提取。土壤定量屬性采用ArcGIS10.2地統(tǒng)計(jì)插值模塊中普通克里格法分別進(jìn)行插值以獲取 30 m× 30 m柵格數(shù)據(jù)。遙感變量中的 NDVI和TVI植被指數(shù)通過Landsat8 OLI影像波段運(yùn)算得到，計(jì)算公式分別如下：

式中，近紅外（Near Infrared，NIR）和紅光（Red）波段反射率為Landsat8 OLI影像利用ENVI5.3軟件進(jìn)行大氣校正、鑲嵌、裁剪等預(yù)處理后獲取的研究區(qū)地表反射率。地形因子中的DEM由福建省23幅ASTER GDEM影像數(shù)據(jù)經(jīng)坐標(biāo)轉(zhuǎn)換、鑲嵌及研究區(qū)腌膜運(yùn)算提取。而其他4個地形因子（坡度、坡向等）是利用ArcGIS10.2軟件Spatial Analyst模塊，基于 DEM 計(jì)算得到。氣象因子?xùn)鸥駭?shù)據(jù)由全省氣候資料數(shù)據(jù)集使用ArcGIS10.2地統(tǒng)計(jì)插值模塊中反距離權(quán)重法（Inverse Distance Weighted，IDW）插值后，再通過研究區(qū)腌膜提取得到，空間分辨率與遙感變量一致，為30 m。將上述環(huán)境變量形成的柵格數(shù)據(jù)集，利用 ArcMap提取到與土壤有機(jī)碳實(shí)測樣點(diǎn)空間匹配的訓(xùn)練數(shù)據(jù)集中（如圖1所示），用于RF模型的構(gòu)建和驗(yàn)證。

RF模型構(gòu)建和預(yù)測的實(shí)現(xiàn)均通過 Python scikit-learn庫中RandomForestRegressor包實(shí)現(xiàn)。變量相對重要性排序可直接調(diào)用工具包中 feature_importances屬性實(shí)現(xiàn)。

基于 RF模型和 OK模型生成耕地土壤有機(jī)碳空間分布圖，用于評價不同模型在SOC空間異質(zhì)性表達(dá)上的優(yōu)劣。針對RF模型，將空間分辨率為30 m的遙感變量、氣候因子?xùn)鸥窈屯寥缹傩詵鸥褚罁?jù)對應(yīng)的變量組合分別輸入 RF-S與 RF-A模型，得到SOC空間分布格局。OK模型則是基于所有SOC樣本點(diǎn)，使用普通克里格方法插值后重采樣為30 m獲取。利用ArcGIS10.2制圖模塊完成SOC空間分布專題制圖。

2 結(jié) 果

2.1 同變量組合下RF模型預(yù)測精度對比

如表2所示，基于全部樣本點(diǎn)，RF-A和RF-S模型精度較好，均呈現(xiàn)高度相關(guān)（r>0.8）和中等變異水平（CV>23%），相對誤差RE小于10%。但與RF-S模型相比，加入了土壤屬性變量的RF-A模型預(yù)測誤差有顯著下降，RMSE和 MAE分別下降45.13%和42.68%，表明加入土壤屬性變量（N、P）有利于提升模型擬合度及預(yù)測精度。

表2 兩種不同變量組合下RF模型的預(yù)測精度Table 2 Prediction accuracies of RF models under two different combinations of variables

如表3所示，兩種RF模型預(yù)測值的均值與SOC實(shí)測值的均值非常接近，約等于 14 g·kg–1，但預(yù)測結(jié)果范圍明顯被壓縮，變異系數(shù)和標(biāo)準(zhǔn)差值均變小。由圖2也可發(fā)現(xiàn)，兩種模型在SOC低值區(qū)（累積百分比40%以下）略高估實(shí)測值，而在SOC高值區(qū)（累積百分比75%以上）明顯低估實(shí)測值，但高值樣點(diǎn)數(shù)（SOC>26 g·kg–1）所占比重較小，所以兩種 RF模型的預(yù)測結(jié)果能夠解釋SOC大部分空間變異?？傮w而言，加入了土壤屬性變量的 RF-A模型預(yù)測值的累積分布圖更接近實(shí)際結(jié)果，能更好地表征區(qū)域SOC的動態(tài)變化范圍。

表3 兩種不同變量組合下RF模型的預(yù)測結(jié)果與SOC實(shí)測值對比Table 3 Comparison SOC measured value with prediction results of RF models under two different combinations of variables

2.2 RF模型環(huán)境變量重要性

各模型最終篩選出的用于土壤有機(jī)碳預(yù)測的變量如表4所示，各類別均有環(huán)境變量被保留參與模型構(gòu)建。由表4可知，在RF-S模型中，重要性最高的環(huán)境變量為DEM（23.22%），氣候因子分列重要性排序第二到五位，雖然單個變量重要性稍弱，但累積貢獻(xiàn)率達(dá)58.20%，依舊占據(jù)主導(dǎo)地位。在加入土壤屬性變量的RF-A模型中，N的重要性超過DEM，位列第一，說明這兩個變量是影響閩東南地區(qū)耕地SOC空間變異的主要協(xié)同因子，且N和SOC的關(guān)系更為密切。值得注意的是，遙感變量中的植被指數(shù)NDVI和TVI，全部通過變量篩選被保留參與模型預(yù)測，盡管單因素貢獻(xiàn)較弱，對空間結(jié)果的預(yù)測也是不可或缺的。各類環(huán)境變量中貢獻(xiàn)率最高的因子分別為地形因子DEM、植被指數(shù)NDVI、氣象因子Mint和土壤屬性N。

表4 RF模型特征變量重要性排序Table 4 Ranking of relative importance of environmental variables of RF model

2.3 基于不同抽樣百分比的精度檢驗(yàn)

在不同梯度抽樣百分比下，對兩種 RF模型以及OK模型的訓(xùn)練和預(yù)測精度進(jìn)行檢驗(yàn)，結(jié)果如表5所示。在訓(xùn)練數(shù)據(jù)集中，不同抽樣百分比下，RF-S與RF-A模型預(yù)測的SOC值與實(shí)測值的相關(guān)系數(shù)r分別在0.75～0.85和0.89～0.92之間，而 OK模型r處于0.60左右，個別情況下出現(xiàn)低值0.55，整體精度明顯低于RF模型。由誤差MAE、RMSE、RE和變異系數(shù)CV值亦可發(fā)現(xiàn)，RF-A模型在各梯度下計(jì)算誤差最小，變異程度最大，OK模型表現(xiàn)最弱。顯然，RF-A模型擬合精度最高，RF-S次之，而OK模型最低。三種模型的擬合能力隨訓(xùn)練樣本量呈現(xiàn)不同的變化趨勢：RF-A模型受樣本數(shù)量影響小，整體表現(xiàn)穩(wěn)定；RF-S模型小幅波動，規(guī)律性不強(qiáng)，在訓(xùn)練樣本為70%時精度最低，其余情況較穩(wěn)定；而OK模型波動較大，在訓(xùn)練樣本為 80%時表現(xiàn)最佳（r=0.65），30%時表現(xiàn)最差（r=0.55）。由于多種環(huán)境變量的協(xié)同作用，RF模型魯棒性較好，而OK模型對采樣點(diǎn)要求較高，且插值結(jié)果隨取樣空間尺度增大會產(chǎn)生明顯的平滑效應(yīng)，這與陳飛香等[20]使用克里格法對土壤原始樣點(diǎn)在不同采樣密度下的插值結(jié)果相一致。所以，當(dāng)采樣點(diǎn)較少時，應(yīng)選擇 RF模型，并優(yōu)先選含有土壤屬性變量的 RF-A模型。在驗(yàn)證數(shù)據(jù)集中，相較于其他兩個模型，RF-S整體精度不高，RF-A模型計(jì)算的r除在20%和40%驗(yàn)證數(shù)據(jù)集較OK模型低0.1，其他情況下均略高于OK模型。總體而言，RF-A模型預(yù)測精度普遍優(yōu)于OK模型，表現(xiàn)最好。

2.4 耕地土壤有機(jī)碳含量空間分布

如圖3a、圖3b、圖3c所示，RF-S、RF-A和OK模型預(yù)測結(jié)果的空間分布均呈現(xiàn)為東部沿海較低、西部內(nèi)陸較高。雖然三種模型的總體趨勢比較相似，但在研究區(qū)南部和中部偏北地區(qū)，OK模型與RF-S模型預(yù)測的SOC高值區(qū)明顯小于RF-A模型。在上述預(yù)測差異較大區(qū)域隨機(jī)選取子區(qū)域（圖3a～圖3c中a、b框）放大顯示（圖3a～圖3c右邊的子圖a、b），可發(fā)現(xiàn)RF-A模型的SOC含量分級區(qū)間數(shù)明顯多于其他兩個模型，且空間變異更強(qiáng)。總體來看，RF-A模型無論在模型精度或空間異質(zhì)性表達(dá)上均為最優(yōu)模型，以下僅對最優(yōu)模型預(yù)測的SOC空間分布格局進(jìn)行分析。

RF-A模型反演得到閩東南區(qū)SOC均值為14.70±2.95 g·kg–1，范圍為 3.63～25.51 g·kg–1，其中 13～19 g·kg–1區(qū)間的面積占比最高，超過研究區(qū)耕地總面積的65%，主要分布在西部內(nèi)陸閩中大山帶戴云山-博平嶺段東南側(cè)；小于 10 g·kg–1和大于 19 g·kg–1的面積占比較低，不足10%，分別分布在閩東南地區(qū)三大平原（漳州平原、泉州平原、莆仙平原）和西部海拔最高地；10～13 g·kg–1區(qū)間所占面積在19%左右，位于高低值過渡帶。

通過數(shù)值分析發(fā)現(xiàn)，SOC空間特征與各類環(huán)境變量中貢獻(xiàn)率最高的四個因子（DEM、NDVI、Mint、N）呈現(xiàn)明顯的相關(guān)性（如圖4所示）。針對 DEM和Mint兩個因子，一般海拔越高，相應(yīng)溫度越低，而SOC含量正是與海拔高度呈正相關(guān)（圖4a），而與年最低溫度Mint呈負(fù)相關(guān)（圖4c）。這與楊順華等[21]的觀點(diǎn)一致，認(rèn)為平原丘陵過渡帶土壤有機(jī)碳與高程等穩(wěn)定因素呈極顯著正相關(guān)。NDVI值越大，植被覆蓋度越高，SOC含量越高（圖4b），這是由于生物量也是土壤有機(jī)碳最重要的來源。值得注意的是，SOC含量在一定范圍內(nèi)隨著N含量的升高而增加，但當(dāng) SOC含量均值>17 g·kg–1時，對應(yīng)N含量呈下降趨勢（圖4d），可能原因是土壤碳氮比（C/N）是一個相對固定的數(shù)值，土壤N素含量極大地影響SOC含量，然而在高海拔地區(qū)，雖氮素平均含量有所減少，但鑒于高海拔地區(qū)低溫降水少的條件下有機(jī)質(zhì)分解速率較低，SOC含量并無下降。

3 討 論

3.1 閩東南地區(qū)土壤有機(jī)碳空間預(yù)測及主要環(huán)境變量影響

基于 RF-A模型預(yù)測得到的閩東南地區(qū)耕地SOC含量東部沿海較低、西部內(nèi)陸較高，這與劉素真[22]利用全國第二次土壤普查數(shù)據(jù)和近期野外采樣數(shù)據(jù)開展的福建省有機(jī)碳含量模擬結(jié)果基本一致。西部地區(qū)因海拔高度較高伴隨溫度較低，用作農(nóng)業(yè)生產(chǎn)的可能性降低，植被逐漸轉(zhuǎn)化為枯落物較多的自然植被，同時土壤微生物分解有機(jī)質(zhì)的速度減慢，礦化作用減弱，從而導(dǎo)致有機(jī)質(zhì)的積累量逐漸增加[10]，而東部低海拔區(qū)域，溫度較高，農(nóng)業(yè)生產(chǎn)較便利，頻繁的耕作加劇土壤擾動，促進(jìn)有機(jī)碳分解，使土壤有機(jī)碳周轉(zhuǎn)速率加快，積累量減少，導(dǎo)致SOC含量較低，這與上述模型計(jì)算結(jié)果相一致。由此可見，RF-A模型預(yù)測的研究區(qū)SOC含量空間分布格局是合理的，從這方面證明了該制圖方法的可行性。

特征變量貢獻(xiàn)率分析顯示，水解性氮（N）是預(yù)測土壤有機(jī)碳含量最重要的環(huán)境變量，這與預(yù)期一致。大量研究已證實(shí)陸地生態(tài)系統(tǒng)碳氮循環(huán)存在緊密耦合[23]，堿解氮與土壤有機(jī)碳存在顯著正相關(guān)關(guān)系[24]，但因數(shù)據(jù)獲取不易，以往較少用于土壤碳空間預(yù)測。Were等[9]和謝恩澤等[25]在RF模型研究中也證明全氮TN是解釋SOC空間變異最重要的變量，這與本文研究結(jié)論基本一致。

遙感變量中的DEM的重要性僅次于N。但與N不同，DEM是通過地形地貌差異間接影響SOC的空間分布，屬于外部因素。與一般區(qū)域耕地分布于相對平緩且交通便利地帶不同，本研究區(qū)位于亞熱帶復(fù)雜地貌區(qū)，耕地分布的海拔差異大，不同高度主要分布的土壤類型也有所不同，這就強(qiáng)調(diào)了DEM對SOC預(yù)測結(jié)果的影響。齊雁冰等[10]在利用RF模型反演陜西省土壤有機(jī)質(zhì)的研究中得出了相似的結(jié)果，即在地形地貌復(fù)雜且耕地面積較小的區(qū)域內(nèi)，DEM對SOC空間預(yù)測的貢獻(xiàn)率相對較高。此外，高程和地形可通過影響降水情況間接影響 SOC含量的空間分布。當(dāng)來自太平洋的暖濕氣流進(jìn)入福建沿海并向西北運(yùn)行過程中首先遇到閩中大山帶，由于地形對氣流的抬升作用，在東南坡產(chǎn)生較多的降水，這與SOC的空間分布高度一致。由此可見，高程和降水、溫度對SOC空間分布的影響存在重疊，并且高程的影響更大，所以降低了氣候因子在模型中的重要性。

在基于遙感影像獲取土壤有機(jī)碳的研究中，植被指數(shù)是最常用的變量之一，而本研究中所選取的兩個植被指數(shù)在 RF模型中的重要性并不突出，但二者均通過 OOB Score的篩選被保留參與模型預(yù)測。從數(shù)值分析上看，SOC含量大致與NDVI呈正相關(guān)關(guān)系（圖4b），這與以往研究相一致[26]。

3.2 RF-A模型精度

研究結(jié)果可得，RF-A模型無論在模型精度或空間異質(zhì)性表達(dá)均為閩東南區(qū) SOC空間預(yù)測的最佳選擇。基于全部1 257個樣點(diǎn)訓(xùn)練的RF-A模型的r=0.95，RRMSE值與RE值分別為13.89%和–5.91%（表2）。Hengl等[27]提出，RRMSE≤40%以內(nèi)為可以接受的模型準(zhǔn)確性，當(dāng) RE≤±10%時，模型模擬結(jié)果處于可接受范圍?？梢?，RF-A模型精確度高，且模型的準(zhǔn)確性處于可接受范圍。相較于 RF-S模型，加入土壤屬性因子的RF-A模型RMSE和MAE分別下降45.13%和42.68%，可見，RF模型的預(yù)測精度受目標(biāo)變量與輔助環(huán)境變量之間的相關(guān)性強(qiáng)弱控制，當(dāng)加入與SOC緊密耦合的內(nèi)部影響因子（如N），RF模型的精度顯著提升。謝恩澤等[25]在對比保留或移除輔助因子 TN對蘇南農(nóng)田土壤有機(jī)質(zhì)空間分布預(yù)測精度的影響時，得到了一致的結(jié)論，進(jìn)一步證實(shí)了土壤因子用于亞熱帶復(fù)雜地貌區(qū)有機(jī)碳空間預(yù)測的可行性。

但RF-A模型在驗(yàn)證數(shù)據(jù)集中的r值（表5）與以往文獻(xiàn)[10-11]相比，沒有明顯的優(yōu)勢?？赡茉蛟谟谘芯繀^(qū)為復(fù)雜地貌區(qū)，耕地面積小而破碎，SOC預(yù)測難度較一般地勢平坦、地形地貌單一的小區(qū)域大，且本文采用的大樣本數(shù)據(jù)也會顯著增加驗(yàn)證結(jié)果的誤差?？傮w而言，RF-A模型在全部樣點(diǎn)和不同抽樣百分比的預(yù)測精度均處于準(zhǔn)確度范圍內(nèi)（不同梯度下RRMSE值均小于34%，見表5），且預(yù)測結(jié)果與實(shí)際相符，可用于研究區(qū)SOC空間分布格局的獲取。

4 結(jié) 論

本文采用隨機(jī)森林模型，選擇易獲取的遙感變量和氣候因子及不易獲取的土壤屬性作為模型輸入，基于大量土壤有機(jī)碳實(shí)測樣點(diǎn)分別訓(xùn)練兩種不同變量組合數(shù)據(jù)集驅(qū)動的 RF模型，對比驗(yàn)證加入與未加入土壤屬性變量構(gòu)建的 RF模型精度，并與普通克里格插值模型進(jìn)行對比。結(jié)果顯示，聯(lián)合遙感變量、土壤屬性和氣候因子共同構(gòu)建 RF-A模型精度最高，可作為預(yù)測該研究區(qū)SOC含量的高效方法，在模型驅(qū)動因子中加入N、P等土壤屬性變量，能顯著改進(jìn)模型擬合度與精度，更好地提升區(qū)域SOC空間異質(zhì)性的預(yù)測能力。相較于普通克里格插值模型，RF-A綜合模型在不同抽樣百分比的預(yù)測精度普遍更優(yōu)，且穩(wěn)定性更好，因此在采樣點(diǎn)較少的情況下，應(yīng)優(yōu)先選擇 RF-A模型。變量貢獻(xiàn)率結(jié)果表明，影響研究區(qū)耕地SOC預(yù)測結(jié)果空間分布最重要的環(huán)境變量是與SOC直接相關(guān)的土壤屬性——水解性氮，其次是DEM。