郭治興袁宇志郭 穎，2孫 慧，2柴 敏

陳澤鵬3Mogens H. Greve4

（1 廣東省生態(tài)環(huán)境技術(shù)研究所廣東省農(nóng)業(yè)環(huán)境綜合治理重點(diǎn)實(shí)驗(yàn)室，廣州 510650）

（2 山西農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院，山西晉中 030800）

（3 廣東省煙草公司，廣州 510610）

（4 Department of Agro-ecology，F(xiàn)aculty of Science and Technology，Aarhus University，8830 Tjele，Denmark）

基于地形因子的土壤有機(jī)碳最優(yōu)估算模型*

郭治興1袁宇志1郭 穎1，2孫 慧1，2柴 敏1

陳澤鵬3Mogens H. Greve4

（1 廣東省生態(tài)環(huán)境技術(shù)研究所廣東省農(nóng)業(yè)環(huán)境綜合治理重點(diǎn)實(shí)驗(yàn)室，廣州 510650）

（2 山西農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院，山西晉中 030800）

（3 廣東省煙草公司，廣州 510610）

（4 Department of Agro-ecology，F(xiàn)aculty of Science and Technology，Aarhus University，8830 Tjele，Denmark）

基于數(shù)字地面模型（Digital Terrain Model，DTM），同時(shí)考慮因子組合和分辨率構(gòu)建土壤有機(jī)碳（SOC）最優(yōu)估算模型。在7 100 km2范圍內(nèi)，選取了71個(gè)分辨率和22個(gè)地形因子中不多于5個(gè)因子的所有可能組合，構(gòu)造了2 514 820個(gè)模型。采樣點(diǎn)隨機(jī)分為兩組，6 362個(gè)訓(xùn)練樣點(diǎn)構(gòu)造數(shù)據(jù)挖掘模型，其他2 208個(gè)為驗(yàn)證樣點(diǎn)。根據(jù)模型相關(guān)系數(shù)r值大小從中選取了不同個(gè)數(shù)因子組合以及相應(yīng)分辨率的最優(yōu)模型，并根據(jù)這些模型生成對(duì)應(yīng)的土壤有機(jī)碳圖。結(jié)果表明：?jiǎn)蝹€(gè)地形因子模型和柵格大小之間的關(guān)系表現(xiàn)出多樣化，并不是分辨率越高模型結(jié)果越好。單因子模型r值的大小并不能決定其在因子組合模型中的重要性。不同的因子及其組合有其特定的最適分辨率，最佳分辨率范圍約為60～150 m。綜合數(shù)據(jù)的存儲(chǔ)空間和計(jì)算量、模型復(fù)雜度、預(yù)測(cè)精度以及空間表達(dá)能力，該地區(qū)最優(yōu)模型由相對(duì)坡位、高程、歸一化高程及多尺度山谷平坦指數(shù)等4個(gè)變量組成，對(duì)應(yīng)分辨率為121.6 m。同時(shí)與多種克里格空間插值方法生成的土壤有機(jī)碳空間分布圖進(jìn)行了對(duì)比分析，發(fā)現(xiàn)無(wú)論幾個(gè)變量的組合，其空間預(yù)測(cè)能力均較克里格空間插值方法更能表達(dá)SOC的空間變化，預(yù)測(cè)精度也較高。

土壤有機(jī)碳；數(shù)字土壤制圖；數(shù)據(jù)挖掘；最優(yōu)估算模型；數(shù)字地面模型；地形參數(shù)

土壤有機(jī)碳（SOC）是全球碳庫(kù)的重要組成部分，是陸地生態(tài)系統(tǒng)中最大的有機(jī)碳庫(kù)［1］，在全球碳循環(huán)和氣候變化過(guò)程中起著極其重要的作用［2-3］。同時(shí)也是人類可持續(xù)利用土地資源進(jìn)行生產(chǎn)，保證全球糧食安全必須要掌握的基礎(chǔ)信息［4］。

土壤有機(jī)碳儲(chǔ)量受自然和人為兩個(gè)方面因素的影響［5］。地球環(huán)境的進(jìn)化決定了一個(gè)特定地區(qū)的氣候、植物、地形、地貌以及土壤等自然屬性，這些屬性之間必定有較高的相關(guān)性。一個(gè)地區(qū)的土壤有機(jī)碳和其所處的地質(zhì)、地貌以及地形等因子是自然形成和演化的共同結(jié)果［6］。在此過(guò)程中，特定的地質(zhì)、地貌以及地形等因子均會(huì)影響土壤有機(jī)碳的形成與演化。因此地形因子的組合與土壤有機(jī)碳必然有其內(nèi)在的聯(lián)系，能夠在一定程度上反映所在地土壤有機(jī)碳的狀態(tài)［7-8］。

土壤有機(jī)碳還受人類活動(dòng)的影響［9-10］。人類利用土壤的方式除受氣候、傳統(tǒng)等因素影響［11］，特定的地形因子也會(huì)直接影響人類利用土地的方式，例如坡度、坡向等。因此，從人類活動(dòng)的角度看，地形因子與土壤有機(jī)碳儲(chǔ)量之間也有其必然的聯(lián)系。

土壤有機(jī)碳的計(jì)算模型較多，大體分為經(jīng)驗(yàn)?zāi)Ｐ停?2-13］、統(tǒng)計(jì)模型［14-15］和機(jī)理模型［16-17］?；揪情g接的方法，目前還沒(méi)有一個(gè)方法可以快速測(cè)量或估算一個(gè)區(qū)域內(nèi)的土壤有機(jī)碳儲(chǔ)量。采用地形因子等相對(duì)較易獲取的參數(shù)建模就是一個(gè)經(jīng)常被采用的方法。在數(shù)字土壤制圖領(lǐng)域，土壤屬性的預(yù)測(cè)和計(jì)算，包括土壤有機(jī)碳的估算，地形因子經(jīng)常會(huì)被納入模型［18-20］。遙感技術(shù)為快速獲取數(shù)字地面模型（DTM），提供了一個(gè)高效的技術(shù)手段［21］。

一個(gè)能夠精確預(yù)測(cè)土壤屬性的模型，還依賴于找到一個(gè)最適柵格大小和由其推導(dǎo)出的地形因子［22-23］。眾所周知遙感數(shù)據(jù)是柵格形式的，而地形因子的取值與柵格的大小息息相關(guān)，模型的結(jié)果也會(huì)隨之變化［24-27］。如何確定這些因子的最適柵格大小的問(wèn)題仍然未解決，只有一些針對(duì)特定區(qū)域的研究和原則性的指南［28］。

通常DTM數(shù)據(jù)是最容易獲取的數(shù)據(jù)源，常規(guī)的做法是直接采用原始分辨率的數(shù)據(jù)，一般認(rèn)為這樣會(huì)提高模型的準(zhǔn)確度和精度［29］。但是不恰當(dāng)?shù)母駯糯笮⊥ǔ?huì)帶來(lái)一些問(wèn)題［30］。如果柵格過(guò)大，地形因子的屬性會(huì)丟失太多的細(xì)節(jié)，進(jìn)而失去預(yù)測(cè)能力。反之，如果柵格過(guò)小，由此導(dǎo)出的地形因子會(huì)帶有太多的細(xì)節(jié)，甚至是“噪聲”，這樣勢(shì)必會(huì)影響最終模型的精度［31］。同時(shí)，高分辨率數(shù)據(jù)由于擁有較高的空間精度，柵格非常小，數(shù)據(jù)量非常大，需要更大的存儲(chǔ)空間；模型會(huì)需要更多的計(jì)算時(shí)間，有的時(shí)候這種增長(zhǎng)是指數(shù)級(jí)別的。一個(gè)合適的分辨率，應(yīng)該是柵格尺寸足夠大，但又可以充分表達(dá)模型所需要的信息［22］。

預(yù)測(cè)土壤屬性采用哪些地形因子也是一個(gè)需要考慮的問(wèn)題［32］。由DTM可以導(dǎo)出很多地形因子，一個(gè)模型采用哪些因子，通常均具有較大的隨意性和主觀性［33-34］。同一個(gè)地形因子對(duì)不同的土壤屬性有著不同的影響。反之，不同土壤屬性有屬于其特有的決定性的地形因子［35］。地形因子之間存在復(fù)雜的相互關(guān)系，單個(gè)因子在土壤模型中的重要性在多個(gè)地形因子組合模型中也會(huì)發(fā)生變化［36-37］。

利用DTM及相關(guān)因子預(yù)測(cè)SOC前人做出了大量工作，但大都基于一個(gè)固定的柵格大小或人為選取的因子或組合［38-39］。本研究試圖通過(guò)DTM及其推導(dǎo)的相關(guān)地形因子，綜合柵格大小、模型精度以及數(shù)據(jù)的冗余度等因素，通過(guò)數(shù)據(jù)挖掘技術(shù)，尋找最優(yōu)因子組合和相匹配的最佳分辨率的土壤有機(jī)碳預(yù)測(cè)模型。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

研究區(qū)位于56°9′～56°33′N及8°6′～10°57′E之間。是一條從西到東貫穿整個(gè)丹麥日德蘭半島（Jutland）的寬約45 km的條狀帶，總面積約7 100 km2（圖1）。整個(gè)區(qū)域地勢(shì)平緩，平均海拔為34 m，最高峰僅為152 m。該區(qū)域冬季和夏季平均溫度分別為0°C和16°C，年平均降雨量為500～800 mm。主要種植冬小麥、春大麥、油菜、牧草以及玉米。

由于歷史上冰川侵蝕活動(dòng)，造成了研究區(qū)內(nèi)復(fù)雜多樣的地貌類型［40］。據(jù)丹麥土壤分類，研究區(qū)范圍內(nèi)的土壤質(zhì)地從東部的重黏土和黏土到壤質(zhì)土到西部的粗砂質(zhì)土。東部地區(qū)的土壤由魏克塞冰期冰磧地貌發(fā)育而成，以壤質(zhì)土為主。而西部主要由嚴(yán)重侵蝕和歷史相對(duì)較長(zhǎng)的撒利冰磧地貌，以及一個(gè)相對(duì)較大的冰水平原組成，以砂質(zhì)土為主。

圖1 研究區(qū)位置Fig. 1 Location of the study area

該區(qū)域包含有丹麥全部8種土壤類型，全部12種景觀類型中10種，全部12種地質(zhì)類型中的11種，具有典型的代表性。

1.2 土壤采樣點(diǎn)和數(shù)據(jù)庫(kù)

研究中用到的土壤樣點(diǎn)主要來(lái)源于全丹麥土壤分類數(shù)據(jù)庫(kù)［41］。該數(shù)據(jù)庫(kù)建成于20世紀(jì)70年代中期，當(dāng)時(shí)建庫(kù)的主要目的是為了土壤制圖和保護(hù)城市化進(jìn)程中日益受到威脅的土地資源，由于該項(xiàng)目要求三年內(nèi)完成，因此其目標(biāo)主要集中在采樣和分析幾個(gè)重要和易分析的指標(biāo)，如土壤質(zhì)地和耕層（0～20 cm）土壤有機(jī)質(zhì)。全國(guó)共采集36 000余個(gè)土壤樣品。在研究區(qū)域內(nèi)共計(jì)有8 570個(gè)樣點(diǎn)數(shù)據(jù)，其中約75%，計(jì)6 362個(gè)樣點(diǎn)用來(lái)生成和優(yōu)選模型，其余25%，計(jì)2 208個(gè)的樣點(diǎn)用來(lái)驗(yàn)證模型。

圖2 采樣點(diǎn)分布圖Fig. 2 Distribution map of sampling sites

1.3 DTM數(shù)據(jù)

丹麥全國(guó)的DTM數(shù)據(jù)由國(guó)家環(huán)境部基于光探測(cè)與測(cè)量（Light Detection and Ranging，LiDAR）數(shù)據(jù)在2005—2007年間生成。采用不規(guī)則三角網(wǎng)方法從LiDAR點(diǎn)云數(shù)據(jù)中提取全國(guó)的地表高程，然后將TIN數(shù)據(jù)轉(zhuǎn)換為1.6m分辨率的柵格數(shù)據(jù)。

采用簡(jiǎn)單求平均值的方法生成其他較低分辨率DTM柵格數(shù)據(jù)集，共采樣計(jì)71個(gè)，小于150 m的設(shè)計(jì)了41組，200～1 000 m之間每50 m取一組，1 000 m以上大約每100 m取一組，最小分組的分辨率為2 304 m?；贒TM計(jì)算其他地形參數(shù)之前，先對(duì)每個(gè)不同分辨率的柵格DTM數(shù)據(jù)進(jìn)行填洼處理，以減少或避免在進(jìn)行地形運(yùn)算，特別是推導(dǎo)與水系網(wǎng)有關(guān)的地形參數(shù)時(shí)，由于細(xì)小洼和峰導(dǎo)致的錯(cuò)誤或出現(xiàn)問(wèn)題［42］。

1.4 導(dǎo)出變量

基于不同分辨率的DTM數(shù)據(jù)集，通過(guò)SAGAGIS計(jì)算其他地形相關(guān)因子。共導(dǎo)出21個(gè)地形因子參數(shù)，每個(gè)因子的說(shuō)明見表1。

1.5 基于規(guī)則的數(shù)據(jù)挖掘模型

本研究使用澳大利亞Rulequest公司的數(shù)據(jù)挖掘工具Cubist 2.08生成SOC預(yù)測(cè)模型，這種模型利用線性的方法模擬非線性的現(xiàn)實(shí)。Cubist設(shè)計(jì)用來(lái)處理海量數(shù)據(jù)，它的數(shù)據(jù)挖掘模型不僅注重模型的預(yù)測(cè)精度，同時(shí)也重視模型的可解釋性。其生成的模型以一組規(guī)則集表達(dá)，每個(gè)規(guī)則采用如果滿足條件則使用一個(gè)線性方程的模式。當(dāng)自變量滿足某一規(guī)則給定的條件時(shí)，它就使用該規(guī)則對(duì)應(yīng)的多變量線性模型進(jìn)行預(yù)測(cè)［43］。

Cubist模型利用相對(duì)坡位（RSP）、多尺度山谷指數(shù)（MRVBF）、高程（Elevation）和歸一化高度（Normalh）四個(gè)變量預(yù)測(cè)SOC含量。該模型由29個(gè)規(guī)則組成，每個(gè)規(guī)則給出一個(gè)條件和一個(gè)對(duì)應(yīng)的多元線性模型。例如，第一個(gè)規(guī)則預(yù)測(cè)1758個(gè)樣點(diǎn)的SOC，這些樣點(diǎn)的取值范圍為0.08～540.1，平均值和估計(jì)誤差分別為20.315和7.662。該規(guī)則列出條件當(dāng)MRVBF＞8.810 939，并且RSP＞0.227 911時(shí)，其對(duì)應(yīng)的用來(lái)預(yù)測(cè)SOC的多變量線性回歸方程。從這條規(guī)則可以看出，這些樣點(diǎn)的SOC隨著RSP的增加而減少呈現(xiàn)負(fù)相關(guān)，隨著MRVBF、Elevation和Normalh的增加而增加。該模型的總體相對(duì)誤差和模型相關(guān)系數(shù)分別為0.54和

表1 模型地形因子列表Table 1 List of topographic factors derived by DTM

2 結(jié) 果

0.71。同時(shí)，模型還給出了變量在模型中貢獻(xiàn)率。貢獻(xiàn)率是指相對(duì)重要性和相對(duì)使用率兩個(gè)參數(shù)，相對(duì)重要性指變量在模型全部規(guī)則所有條件中出現(xiàn)的頻率，相對(duì)使用率指變量在多元線性方程中參與計(jì)算的頻率。在這個(gè)模型中RSP的相對(duì)重要性和相對(duì)使用率分別為89%和81%，較Normalh的9%和52%要高很多。

所有模型均是簡(jiǎn)單性和精度之間的平衡。22個(gè)變量的所有可能的組合數(shù)高達(dá)4 194 303個(gè)，變量過(guò)多會(huì)引入太多的噪聲，增加模型的復(fù)雜程度。綜合模型復(fù)雜度和數(shù)據(jù)的計(jì)算量及存儲(chǔ)空間，選取不超過(guò)5個(gè)變量的組合是比較合適的。這種組合共計(jì)有35 442個(gè)，每個(gè)組合有71個(gè)分辨率，共計(jì)生成2 514 820個(gè)模型并進(jìn)行計(jì)算。

以上所有組合通過(guò)Cubist運(yùn)算獲得模型結(jié)果，每個(gè)模型均有對(duì)應(yīng)的相關(guān)系數(shù)r值。對(duì)比所有模型r值，最大r值所在的因子組合與分辨率即為該因子組合和最佳分辨率。

2.1 數(shù)據(jù)基本統(tǒng)計(jì)特征

整體樣本的平均值和標(biāo)準(zhǔn)差分別為36.0 g kg-1和69.7 g kg-1。最小值、中位數(shù)和最大值分別為0.08 g kg-1、15.58 g kg-1和540.1 g kg-1。變異系數(shù)高達(dá)193.6%，說(shuō)明數(shù)據(jù)整體變化幅度較大，整體為帶偏度的正態(tài)分布。盡管兩組數(shù)據(jù)的很多統(tǒng)計(jì)量有所不同，但通過(guò)柯爾莫諾夫—斯米爾諾夫檢驗(yàn)（Kolmogorov-Smirnov test），訓(xùn)練數(shù)據(jù)組和驗(yàn)證數(shù)據(jù)組之間并不存在顯著性差異（99.9%）［44］，說(shuō)明兩組數(shù)據(jù)來(lái)自同一母體，驗(yàn)證數(shù)據(jù)可以用來(lái)對(duì)模型進(jìn)行驗(yàn)證。

2.2 地形因子與分辨率對(duì)模型的影響

表2中列出了單因子在所有分辨率中的模型的r平均值、最大值以及最大值對(duì)應(yīng)的分辨率。單個(gè)地形變量模型和柵格大小之間的關(guān)系表現(xiàn)出多樣化，并不總是負(fù)相關(guān)。大部分變量與格柵大小呈偏態(tài)的正態(tài)分布，如相對(duì)坡位（RSP）和水系網(wǎng)高程（Chnl_alti）。個(gè)別變量呈正相關(guān)，如匯指數(shù)（Convergence），隨著柵格的增大，模型的r值增大。也有部分變量與柵格的大小無(wú)關(guān)，呈隨機(jī)狀態(tài)，如水系網(wǎng)基準(zhǔn)面（Chnl_base），還有部分變量與柵格沒(méi)有任何關(guān)系，其模型的r值在任何分辨率下均為0，如坡向。最好的地形因子為相對(duì)坡位（RSP），r平均值為0.272，最大值達(dá)到0.5。最差的為坡向，無(wú)論是何種分辨率，坡向的r值均為0，即坡向與SOC沒(méi)有任何關(guān)系。

表2 單因子模型r值及最佳分辨率Table 2 R and optimal resolution of single factor models

對(duì)單個(gè)地形因子模型，無(wú)論其與柵格大小的關(guān)系是正相關(guān)還是負(fù)相關(guān)，r最大值時(shí)的分辨率通常并不是最大或最小柵格。不同的變量有一個(gè)或多個(gè)最佳的柵格大小，例如，相對(duì)坡位（RSP）的最佳分辨為92.8m，歸一化高度（Normalh）則有2個(gè)分辨率，分別為12.8和16m。

任意多個(gè)因子的組合與單個(gè)因子的表現(xiàn)不同，多個(gè)因子的組合表現(xiàn)出了一致的規(guī)律，均為呈偏態(tài)的正態(tài)分布。從最高分辨率開始隨著柵格的增大，模型的r值增加至一個(gè)最大值，此處的分辨率范圍大約為 60～150 m。此后r與分辨率呈負(fù)相關(guān)，在1 500m附近r值基本達(dá)到穩(wěn)定，之后隨著分辨率的改變，r值呈隨機(jī)變化狀態(tài)。因此無(wú)論使用哪些因子組合建模，分辨率均應(yīng)在60～150m之間選取。

圖3 變量組合模型相關(guān)系數(shù)r值與分辨率的關(guān)系Fig. 3 Relationship between r of combination of variables with resolution

單因子模型r值的大小并不能決定其在因子組合模型中的重要性。即2個(gè)r值最大的單因子模型，它們組合后的模型的r值并不是最大的。單個(gè)因子r值最大的6個(gè)因子分別為相對(duì)坡位（RSP）、水系網(wǎng)高程（Chnl_alti）、山谷深度（Vall_depth）、歸一化高度（Standh）、地表曲率（Convexity）和高程（Elevation）（表2）。但是在篩選出的最佳因子組合模型中出現(xiàn)次數(shù)最多的6個(gè)因子是：水系網(wǎng)基準(zhǔn)面（Chnl_base）、高程（Elevation）、水系網(wǎng)高程（Chnl_alti）、多尺度山脊指數(shù)（MRRTF）、歸一化高度（Normalh）和相對(duì)坡位（RSP）。事實(shí)上在最佳模型中出現(xiàn)的12個(gè)因子，并不是最好的12個(gè)單因子。最明顯的例子是多尺度山脊指數(shù)（MRRTF）因子，該因子與分辨率關(guān)系呈隨機(jī)分布狀態(tài)，r值總體上也偏低，但是在因子組合模型中該因子共出現(xiàn)了4次。

無(wú)論何種分辨率，因子的多少?zèng)Q定了r的總體趨勢(shì)，引入的因子越多，r值越大。當(dāng)引用全部22個(gè)因子時(shí)所構(gòu)建的模型，r值達(dá)到最大0.78（圖3），但是模型的預(yù)測(cè)能力并非如此。

2.3 SOC最優(yōu)估算模型

通過(guò)對(duì)比r值，篩選出不同數(shù)量因子組合的最優(yōu)模型和對(duì)應(yīng)的分辨率。不同個(gè)數(shù)因子的最優(yōu)模型列于表3，這些模型的r值變化見圖4。單個(gè)因子的模型為相對(duì)坡位（RSP）分辨率為 92.8 m。2個(gè)因子由相對(duì)坡位（RSP）和水系網(wǎng)基準(zhǔn)面（Chnl_ base）組成，最適分辨率為 60.8 m。3個(gè)因子的模型為水系網(wǎng)高程（Chnl_alti），水系網(wǎng)基準(zhǔn)面（Chnl_base）和多尺度山谷指數(shù)（MRVBF），最佳分辨率為64 m。有6組4個(gè)因子的組合模型其r均為0.71，分辨率從64 m至136 m不等。5個(gè)因子共有2個(gè)組合模型，r值為0.78，分辨率為152 m。每組均包含水系網(wǎng)高程（Chnl_alti）、高程（Elevation）、水系網(wǎng)基準(zhǔn)面（Chnl_base）、多尺度山脊指數(shù)（MRRTF）四個(gè)因子，此外歸一化高度（Normalh）和中坡位（Midslppst）2個(gè)因子分在不同組中。同時(shí)還計(jì)算了全部22個(gè)因子組合模型的r值為0.78，最佳分辨率為40 m或64 m。

表3 不同數(shù)量因子最佳組合及分辨率Table 3 Optimal combination of different numbers of factors and resolution

為了確定最終采用的模型，使用預(yù)留的25% 計(jì)2 208個(gè)樣點(diǎn)對(duì)這些預(yù)選的模型的預(yù)測(cè)能力和有效性進(jìn)行了檢驗(yàn)。r值有著明顯的差異，22個(gè)因子模型的r值最差，甚至小于單個(gè)因子的r值（表3）。這可能是由于因子太多，引入的噪聲也隨之加大，從而無(wú)法真正總結(jié)內(nèi)在的規(guī)律。驗(yàn)證點(diǎn)的r值也是隨著因子數(shù)量的增多而增加，但r值最大的是4個(gè)因子的組合，共有2組，分辨率分別為64 m 和121.6 m?？紤]到數(shù)據(jù)的存儲(chǔ)空間前者是后者的4倍，模型運(yùn)算時(shí)間指數(shù)級(jí)增長(zhǎng)，因此本文選擇相對(duì)坡位（RSP）、高程（Elevation）、歸一化高度（Normalh）和多尺度山谷指數(shù)（MRVBF）組成，柵格大小為121.6 m的模型為最優(yōu)模型。

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

為了表達(dá)不同模型對(duì)原始樣點(diǎn)SOC的再現(xiàn)能力和相互間的對(duì)比，所有圖均采用相同色系和分級(jí)值。分級(jí)值為11、14、18和28，分別為全部原始8 570個(gè)樣點(diǎn)SOC值按從小到大排列后，20%、40%、60%和80%的分段值。這些模型生成的SOC總體空間分布格局是一致的。SOC的高值基本發(fā)生在濕地地區(qū)。除濕地外，SOC總體上西部較東部高。從空間上對(duì)比這些模型生成的SOC空間分布圖（圖5），1個(gè)和2個(gè)因子生成的圖形明顯綜合過(guò)度，粒度偏大，22個(gè)變量的圖形中噪聲非常明顯，4個(gè)和5個(gè)變量的綜合效果最好。

圖4 最優(yōu)模型r曲線Fig. 4 Curve of r with resolutions of optimae model

圖5 不同模型參數(shù)生成的土壤有機(jī)碳圖Fig.6 Predicted maps of SOC

表4 克里格方法與模型模擬結(jié)果對(duì)比Table 4 Results comparison of fitting of the models with Kriging

2.5 與克里格預(yù)測(cè)模型對(duì)比

為了驗(yàn)證模型的預(yù)測(cè)能力和空間表達(dá)能力，與通常采用的克里格空間插值方法生成的空間分布圖進(jìn)行了對(duì)比。采用全局泛克里格、全局普通克里格、全局克里格以及普通克里格等多種空間插值方法生成空間分布圖（圖6）。利用2208個(gè)驗(yàn)證點(diǎn)驗(yàn)證其預(yù)測(cè)精度，r和均方根誤差（RMSE）范圍分別為0.34～0.47，64.6～67.3，與模型的結(jié)果相比差距較大。從空間分布圖也可以明顯看出，空間異質(zhì)性未得到合理的表達(dá)，克里格方法生成的空間數(shù)據(jù)趨向于將數(shù)據(jù)進(jìn)行均一化。從表4中的最大、最小以及方差也可以看出克里格生成的數(shù)據(jù)非常集中。

圖6 不同克里格方法生成的土壤有機(jī)碳空間分布圖Fig. 6 Predicted SOC distribution map relative to Kriging method

3 結(jié) 論

研究通過(guò)數(shù)據(jù)挖掘的方法，尋找最優(yōu)的因子組合及最適柵格大小，構(gòu)建丹麥典型區(qū)域土壤有機(jī)碳估算模型。結(jié)果表明不同的單個(gè)地形因子模型與柵格大小有不同的規(guī)律。絕大部分變量預(yù)測(cè)的r

值與格柵大小呈負(fù)相關(guān)，只有個(gè)別的呈正相關(guān)。

r最大值通常并未像預(yù)料的那樣出現(xiàn)在最大分辨率處。多個(gè)變量組合的模型與柵格大小呈偏態(tài)的正態(tài)分布，不同的變量組合最佳分辨率的位置也不同。最大的r值并不是出現(xiàn)在最高分辨率處，而是出現(xiàn)在60～150 m附近，因此無(wú)論何種因子的組合，選擇的分辨率不能低于200 m?？傮w上，變量越多，模型的r值越大，4到5個(gè)變量的組合的效果均能較好地預(yù)測(cè)SOC。但是由于噪聲的影響，超過(guò)4個(gè)變量的組合后，變量越多對(duì)未知數(shù)據(jù)的預(yù)測(cè)能力越差。綜合數(shù)據(jù)的存貯量和計(jì)算量、模型復(fù)雜度、預(yù)測(cè)精度以及空間表達(dá)能力，由4個(gè)因子相對(duì)坡位（RSP）、高程（Elevation）、歸一化高度（Normalh）和多尺度山谷指數(shù)（MRVBF）組合的，對(duì)應(yīng)的分辨率121.6 m為SOC最優(yōu)估算模型。

［1］Eswaran H，Vandenberg E，Reich P. Organic-carbon in soils of the world. Soil Science Society of America Journal，1993，57（1）：192—194

［2］Lal R. Soil carbon sequestration to mitigate climate change. Geoderma，2004，123（1/2）：1—22

［3］Post W M，Emanuel W R，Zinke P J，et al. Soil carbon pools and world life zones. Nature，1982，298 （5870）：156—159

［4］Lal R. Soil carbon sequestration impacts on global climate change and food security. Science，2004，304 （5677）：1623—1627

［5］Melillo J，Steudler P，Aber J，et al. Soil warming and carbon-cycle feedbacks to the climate system. Science，2002，298（5601）：2173—2176

［6］Nichols J. Relation of organic carbon to soil properties and climate in the southern great plains. Soil Science Society of America Journal，1984，48（6）：1382— 1384

［7］Hontoria C，Saa A，Rodríguez-Murillo J. Relationships between soil organic carbon and site characteristics in peninsular spain. Soil Science Society of America Journal，1999，63（3）：614—621

［8］Percival H J，Parfitt R L，Scott N A. Factors controlling soil carbon levels in new zealand grasslands is clay content important? Soil Science Society of America Journal，2000，64（5）：1623—1630

［9］Stevens A，van Wesemael B. Soil organic carbon stock in the Belgian Ardennes as affected by afforestation and deforestation from 1868 to 2005. Forest Ecology and Management，2008，256（8）：1527—1539

［10］李凌浩. 土地利用變化對(duì)草原生態(tài)系統(tǒng)土壤碳儲(chǔ)量的影響. 植物生態(tài)學(xué)報(bào)，1998，22（4）：300—302

Li L H. Effects of land-use change on soil carbon storage in grassland ecosystems（In Chinese）. Acta Phytoecologica Sinica，1998，22（4）:300—302

［11］Fisher M，Rao I，Ayarza M，et al. Carbon storage by introduced deep-rooted grasses in the south American Savannas. Nature，1994，371（6494）：236—238

［12］Burke I C，Yonker C，Parton W，et al. Texture，climate，and cultivation effects on soil organic matter content in US grassland soils. Soil Science Society of America Journal，1989，53（3）：800—805

［13］于東升，史學(xué)正，孫維俠，等. 基于1∶100萬(wàn)土壤數(shù)據(jù)庫(kù)的中國(guó)土壤有機(jī)碳密度及儲(chǔ)量研究. 應(yīng)用生態(tài)學(xué)報(bào)，2005，16（12）：2279—2283

Yu D S，Shi X Z，Sun W X，et al. Estimation of China soil organic carbon storage and density based on 1∶1 000 000 soil database（In Chinese）. Chinese Journal of Applied Ecology，2005，16（12）：2279—2283

［14］Sims Z，Nielsen G. Organic carbon in montana soils as related to clay content and climate. Soil Science Society of America Journal，1986，50（5）：1269—1271

［15］解憲麗，孫波，周慧珍，等. 中國(guó)土壤有機(jī)碳密度和儲(chǔ)量的估算與空間分布分析. 土壤學(xué)報(bào)，2004，41 （1）：35—43

Xie X L，Sun B，Zhou H Z，et al. Organic carbon density and storage in soils of China and spatial analysis （In Chinese）. Acta Pedologica Sinica，2004，41 （1）：35—43

［16］Franko U，Oelschl?gel B，Schenk S. Simulation of temperature，waterand nitrogen dynamics using the model candy. Ecological Modelling，1995，81（1）：213—222

［17］Li C，F(xiàn)rolking S，Crocker G J，et al. Simulating trends in soil organic carbon in long-term experiments using the dndc model. Geoderma，1997，81（1）：45—60

［18］Hao Y，Lal R，Owens L，et al. Effect of cropland management and slope position on soil organic carbon pool at the north appalachian experimental watersheds. Soil and Tillage Research，2002，68（2）：133—142

［19］McBratney A B，Mendon?a Santos M L，Minasny B. On digital soil mapping. Geoderma，2003，117（1/2）：3—52

［20］陸曉宇，張洪江，程金花，等. 地形因子與人工林土壤有機(jī)碳密度的關(guān)系——以晉西黃土丘陵區(qū)為例. 東北林業(yè)大學(xué)學(xué)報(bào)，2013，41（1）：48—50，58

Lu X Y，Zhang H J，Cheng J H，et al. Ralationship of topographic factors and soil organic carbon density under plantations：A case in loess hilly region in western Shanxi ProvinceChina（In Chinese）. Journal of Northeast Forestry University，2013，41（1）：48—50，58

［21］Mulder V L，de Bruin S，Schaepman M E，et al. The use of remote sensing in soil and terrain mapping — A review. Geoderma，2011，162（1/2）：1—19

［22］Cavazzi S，Corstanje R，Mayr T，et al. Are fine resolution digital elevation models always the best choice in digital soil mapping? Geoderma，2013，195/196：111—121

［23］Gessler P，Chadwick O，Chamran F，et al. Modeling soil-landscape and ecosystem properties using terrain attributes. Soil Science Society of America Journal，2000，64（6）：2046—2056

［24］Behrens T，Zhu A X，Schmidt K，et al. Multi-scale digital terrain analysis and feature selection for digital soil mapping. Geoderma，2010，155（3/4）：175—185

［25］Greve M H，Greve M B，B?cher P K，et al. Generating a danish raster-based topsoil property map combining choropleth maps and point information. Geografisk Tidsskrift-Danish Journal of Geography，2007，107 （2）：1—12

［26］Kim D，Zheng Y. Scale-dependent predictability of dem-based landform attributes for soil spatial variability in a coastal dune system. Geoderma，2011，164 （3/4）：181—194

［27］Smith M P，Zhu A X，Burt J E，et al. The effects of DEM resolution and neighborhood size on digital soil survey. Geoderma，2006，137（1/2）：58—69

［28］Shi X，Girod L，Long R，et al. A comparison of lidarbased dems and usgs-sourced DEMs in terrain analysis for knowledge-based digital soil mapping. Geoderma，2012，170：217—226

［29］Murphy P N C，Ogilvie J，Meng F-R，et al. Modelling and mapping topographic variations in forest soils at high resolution：A case study. Ecological Modelling，2011，222（14）：2314—2332

［30］Thompson J A，Bell J C，Butler C A. Digital elevation model resolution：Effects on terrain attribute calculation and quantitative soil-landscape modeling. Geoderma，2001，100（1/2）：67—89

［31］Sorensen R，Seibert J. Effects of DEM resolution on the calculation of topographical indices：TWI and its components. Journal of Hydrology，2007，347 （1/2）：79—89

［32］Grunwald S. Multi-criteria characterization of recent digital soil mapping and modeling approaches. Geoderma，2009，152（3/4）：195—207

［33］Florinsky I V，Eilers R G，Manning G R，et al. Prediction of soil properties by digital terrain modelling. Environmental Modelling & Software，2002，17 （3）：295—311

［34］Simbahan G C，Dobermann A，Goovaerts P，et al. Fine-resolution mapping of soil organic carbon based on multivariate secondary data. Geoderma，2006，132 （3/4）：471—489

［35］Peng Y，Knadel M，Gislum R，et al. Predicting soil organic carbon at field scale using a national soil spectral library. Journal of Near Infrared Spectroscopy，2013，21（3）：213—222

［36］Adhikari K，Kheir R B，Greve M B，et al. Comparing kriging and regression approaches for mapping soil clay content in a diverse danish landscape. Soil Science，2013，178（9）：505—517

［37］Greve M H，Kheir R B，Greve M B，et al. Quantifying the ability of environmental parameters to predict soil texture fractions using regression-tree model with GIS and lidar data：The case study of Denmark. Ecological Indicators，2012，18（7）：1—10

［38］Kheir R B，Greve M H，Bocher P K，et al. Predictive mapping of soil organic carbon in wet cultivated lands using classification-tree based models：The case study of Denmark. Journal of Environmental Management，2010，91（5）：1150—1160

［39］劉京，朱阿興，張淑杰，等. 基于樣點(diǎn)個(gè)體代表性的大尺度土壤屬性制圖方法. 土壤學(xué)報(bào)，2013，50（1）：12—20

Liu J，Zhu A X，Zhang S J，et al. Large-scaled soil attribute mapping method based on individual representativeness of sample sites（In Chinese）.Acta Pedologica Sinica，2013，50（1）：12—20

［40］Schou A，Nielsen N. Atlas of Denmark：Text and photographs//Hagerup I H.The landscapes. 1949

［41］Madsen H B，Norr A H，Holst K A. The Danish soilclassification atlas of Denmark. Copenhagen：Rosendahl Press，1992

［42］Wang L，Liu H. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. International Journal of Geographical Information Science，2006，20 （2）：193—213

［43］Minasny B，McBratney A B. Regression rules as a tool for predicting soil properties from infrared reflectance spectroscopy. Chemometrics and Intelligent Laboratory Systems，2008，94（1）：72—79

［44］Holliday I E. Kolmogorov-Smirnov test（v1.0.3）in free statistics software（v1.1.23-r7）. Office for Research Development and Education，2012

Optimal Estimation Model of Soil Organic Carbon Based on the Terrain Factor

GUO Zhixing1YUAN Yuzhi1GUO Yin1，2SUN Hui1，2CHAI Min1CHEN Zepeng3Mogens H. Greve4
（1 Guangdong Key Laboratory of Agricultural EnvironmentPollution Integrated Control，Guangdong Institute of Eco-environmental and Soil Sciences，Guangzhou 510650，China）
（2 College of Resources and Environment，Shanxi Agriculture University，Jinzhong，Shanxi 030800，China）
（3 Guangdong Tobacco Company，Guangzhou 510610，China）
（4 Department of Agro-ecology，F(xiàn)aculty of Science and Technology，Aarhus University，8830 Tjele，Denmark）

【Objective】As an important component of the global carbon pool，soil organic carbon （SOC）is the largest organic carbon pool in the terrestrial ecosystem and plays an extremely important role in the global carbon cycle and global warming. The SOC pool is subject to the impacts of both natural and human activities and sure closely related to terrain attributes or factors. There are a number of methodsfor calculation of SOC，which can roughly be sorted into three types，that is，empirical，statistical and mechanismones. But none of them can be used to predict or calculate reapidly soil organic carbon pool of a region rapidly. Remote sensing is an efficient technical means for fast acquisition of DTM，from which numerous information can be derived with the aid of GIS，thus making it possible to constitute a model for rapid calculation of SOC. 【Method】Based on the Digital Terrain Model（DTM）and the topographic attributesderived thereof，an optimal SOC prediction model was built up，taking into account factor combination and resolution with Cubist，a powerful data mining tool for generating rule-based models. This tool works on conditionspecific rules where the output is a set of rules and each rule has a specific multivariate linear model attached. Whenever a situation matches the condition of a rule，the associated model is used to calculate or predictevalues. A total of 8 570 soil samplescollected from the 7 100 km2study area were divided into two groups randomly，6 362 for training and the other 2 208 for model validation，a total of 2 514 820 models were constructed based on 71 selected resolutions and all possible combinations of no more than 5of the 22 terrains attributes. According to the correlation coefficient（R），terrain factors，varying in number，were selected，to form optimal models with their corresponding resolutions，Based on these models，SOC maps were plotted.【Result】Results show that the relationsships between resolution and single-factor models are diversified，it is not true that the higher the resolution，the better the model. The R value of a single-factor model is not necessarily the factor that determines its importance in a multi-factor model. All the multi-factor modelsexhibit a similar rule of skewed normal distribution. Each factor and its combination has a factorspecific optimal resolution，varying in the range of 60～150 m. For models composed of whatever factors，the resolution t be selected should not be lower than 200 m. The variable of the optimal single-factor model is RSP，with resolution being 92.8 m，the variables of the optimal two-factor model are RSP and Chnl_base with resolution being 60.8 m，while the variables of the three-factor model are Chnl_alti，Chnl_base andMRVBF，with resolution being 64 m. There are 6 four-factor models，with R being 0.71 and resolution varying in the range of 64 ～ 136 m，and 2 five-factor models with R being 0.78，and resolution being 152meters. Every model consists at least of Chnl_alti，elevation，Chnl_base and MRRTF，and in the 2 five-factor modelsNormalh or Midslppst is added. The R of all the models consisting of any four ro five of the 22 factors was calculated to be 0.78 with two optimal resolutions，i.e. 40 and 64 m. In general，the more the variables，the higher the R of the models. But owing to impact of the noise，models with more than four factors decine in predictive ability.Four to five is the appropriate number of factors in combination，making the models more capable of predicting SOC. Comparative analysis of the SOC maps plotted with the aid of global fan Kriging，global ordinary Kriging，global kriging and ordinary Kriging shows that regardless of the number of factors in the model，this method is better than all the four Kriging interpolation methods in prediction of spatial variation of SOC and prediction accuracy.【Conclusion】Takinginto comprehensive account storage space，amount of calculationm，sophistication of the model，accuracy of prediction and ability of spatial expression，the optimal model for the study region should be the four-factor model，consisting of relative slope position，elevation，normalheight and MRVBF，with resolution being 121.6 m.

Soil organic carbon（SOC）；Digital soil mapping；Data mining；Digital Terrain Model （DTM）；Terrain attribute；Optimal model

S158

A

10.11766/trxb201608150111

（責(zé)任編輯：檀滿枝）

* 廣東省科技計(jì)劃項(xiàng)目（2015B070701017，2014A040401059，2015A030401068）、國(guó)家自然科學(xué)基金項(xiàng)目（41601558）、廣東省科學(xué)院創(chuàng)新平臺(tái)建設(shè)專項(xiàng)資助 Supported by the Project of the Science and Technology of Guangdong Province（Nos. 2015B070701017，2014A040401059，2015A030401068），National Natural Science Fund（No. 41601558）and SPICC Program（The Scientific Platform and Innovation Capability Construction Program of GDAS）

郭治興（1966—），男，山西太谷人，博士，研究員，主要從事GIS應(yīng)用研究。E-mail：zxguo@soil.cn

2016-08-15；

2016-11-09；優(yōu)先數(shù)字出版日期（www.cnki.net）：2017-01-03