基于三種空間預(yù)測模型的海南島土壤有機(jī)質(zhì)空間分布研究*

姜賽平 張懷志 張認(rèn)連 李兆君 謝良商 徐愛國?

（1 中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所，北京 100081）

（2 海南省農(nóng)業(yè)科學(xué)院農(nóng)業(yè)環(huán)境與土壤研究所，?？?571100）

海南島處于熱帶地區(qū)，島內(nèi)地形復(fù)雜，土地利用多樣，具有豐富的熱帶作物資源，是我國重要的熱帶作物生產(chǎn)基地。研究該區(qū)土壤有機(jī)質(zhì)（SOM）含量的空間變異規(guī)律，對于了解熱帶氣候條件下、不同熱帶作物種植方式下土壤肥力的空間分布狀況具有重要意義。

近些年來，隨著“3S”等信息技術(shù)的快速發(fā)展，數(shù)字土壤屬性制圖方法被廣泛應(yīng)用到土壤學(xué)領(lǐng)域中，用以描述土壤屬性的空間變異規(guī)律。常見的數(shù)字土壤屬性制圖方法有地統(tǒng)計學(xué)制圖方法、機(jī)器學(xué)習(xí)方法和混合模型方法等。地統(tǒng)計學(xué)制圖方法是應(yīng)用最廣泛也是最為成熟的方法。其中，普通克里格法（Ordinary kriging，OK）和回歸克里格法(Regression-kriging，RK)最具代表性。OK因方法簡單、容易操作而被普遍接受。但土壤屬性的空間變異同時受地形因子、土地利用類型、土壤類型等眾多因素的影響［1-3］。由于未考慮到輔助變量對土壤屬性的影響，因而成圖只能描述土壤屬性的整體空間分布規(guī)律，對局部信息描述不夠詳細(xì)，且平滑效應(yīng)嚴(yán)重［4］。為了更好地揭示土壤屬性空間變異規(guī)律，結(jié)合輔助變量的RK被應(yīng)用于土壤屬性制圖當(dāng)中。趙永存等［5］采用多元線性回歸，泛克里格法和RK分別對河北省土壤有機(jī)碳的空間分布規(guī)律進(jìn)行預(yù)測，結(jié)果表明，結(jié)合地形屬性的RK預(yù)測精度最高，且對土壤有機(jī)碳局部變異信息描述地更加詳細(xì)。連綱等［3］采用該法結(jié)合地形因子與遙感指數(shù)對黃土高原丘陵溝壑區(qū)土壤屬性的空間分布規(guī)律進(jìn)行預(yù)測，結(jié)果表明，該法不但能夠提高制圖精度，且能消除部分平滑效應(yīng)。楊順華等［6］選取相對高程和匯流動力指數(shù)為輔助變量采用該法對湖北省枝城鎮(zhèn)SOM空間分布規(guī)律進(jìn)行預(yù)測，亦取得較高的精度。RK能夠較好地描述地形復(fù)雜區(qū)土壤屬性的空間分布規(guī)律，但該法要求數(shù)據(jù)滿足正態(tài)分布和內(nèi)蘊假設(shè)［7］，給實際問題的處理帶來一定的困難。支持向量機(jī)是一種基于統(tǒng)計學(xué)習(xí)理論的機(jī)器學(xué)習(xí)方法，相對于地統(tǒng)計學(xué)方法，它不需要任何前提假設(shè)，在解決小樣本、高維數(shù)、非線性以及局部極小點等問題方面有著較高的性能，但其要求輔助變量為連續(xù)型變量，且需對輔助變量進(jìn)行標(biāo)準(zhǔn)化處理，減小了輔助變量空間分布的異質(zhì)性，影響輔助變量與土壤屬性間的協(xié)同變化關(guān)系［8］。目前該法在數(shù)字土壤屬性制圖中應(yīng)用相對較少。神經(jīng)網(wǎng)絡(luò)模型、分類與回歸樹模型是近些年來發(fā)展較為迅速的機(jī)器學(xué)習(xí)方法，對于連續(xù)型變量和類別型變量均適用，且不需對數(shù)據(jù)進(jìn)行預(yù)處理，能夠很好地描述土壤屬性與環(huán)境變量間的非線性關(guān)系，雖有研究表明二者能夠提高制圖精度［9-10］，但同時也存在著一些缺陷，例如：神經(jīng)網(wǎng)絡(luò)模型的參數(shù)較多，不易確定［11］，模型的預(yù)測結(jié)果也不易解釋［12］。分類與回歸樹模型的結(jié)構(gòu)取決于樣本數(shù)據(jù)，對于樣本量較小的數(shù)據(jù)，訓(xùn)練樣本的微小變化可能導(dǎo)致分類節(jié)點的不同［13］。并且二者均易過度擬合［14-15］。因此并未廣泛應(yīng)用至土壤屬性制圖當(dāng)中。隨機(jī)森林 (Random Forest，RF)模型在分類與回歸樹模型基礎(chǔ)上發(fā)展而來，模型參數(shù)相對較少，計算簡單，且不易過度擬合［16］，被引入數(shù)字土壤屬性制圖當(dāng)中，在地形相對復(fù)雜地區(qū)，能夠取得較高的制圖精度［17-18］。

目前涉及海南全島的SOM含量空間變異特征研究較少，結(jié)合土地利用類型、土壤類型等信息對不同模型的模擬結(jié)果進(jìn)行對比，從而選出適合該區(qū)SOM含量空間變異模型的研究更是鮮見［19-20］。本文根據(jù)地形復(fù)雜區(qū)模型適用情況，結(jié)合土地利用類型、土壤類型和地形因子等輔助變量，選用OK、RK和RF三種方法預(yù)測該區(qū)SOM含量的空間分布特征，通過精度驗證確定出最優(yōu)模型，從而為地處熱帶、較大區(qū)域尺度地形復(fù)雜區(qū)SOM含量空間變異制圖方法的選取提供依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

海南島位于1 8°1 0′～2 0°1 0′N，108°37′～111°3′E之間，是我國第二大島，面積32 900 km2。該島屬于熱帶季風(fēng)氣候，年均氣溫23～25℃，≥10℃年積溫8 200～9 200℃，年平均降水量為1 720mm。島內(nèi)地形復(fù)雜，由山地、丘陵、臺地、平原組成中間高聳、四周低平的環(huán)形層狀地貌。成土母質(zhì)有花崗巖、砂頁巖、淺海沉積物、玄武巖等10種。土壤類型包括磚紅壤、水稻土、赤紅壤、燥紅土、風(fēng)沙土、紫色土等15個土類。土地利用主要為水田、旱地、林地，島內(nèi)以種植熱帶作物為主，種類繁多。

1.2 樣品采集與分析

土壤樣品的采集時間為2012年11月—2012年12月，采用近似網(wǎng)格法，綜合考慮當(dāng)?shù)氐牡匦蔚孛?、土壤類型和土地利用類型等因素，以樣點具有代表性、均勻分布性為原則，在全島19個縣市布點，采樣間距平均為10 km，共采集樣點163個，每個樣點按0～5、0～20、20～40、40～60 cm分層采樣，分別研究表層、耕層、中層和中下層SOM含量空間分布狀況。其中，0～5和0～20 cm分別取樣，0～5 cm土鉆取樣點緊鄰0～20 cm土鉆點，記錄樣點經(jīng)緯度、地形和土地利用類型等信息。土壤類型包括磚紅壤、水稻土、赤紅壤、燥紅土、火山灰土、紫色土、風(fēng)沙土、石灰（巖）土等8個土類。土地利用類型包括水田、旱地、園地。

土壤樣品SOM含量的測定采用重鉻酸鉀-硫酸消化法，具體步驟參見文獻(xiàn)［21］。

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

輔助變量：RK和RF均在6個地形因子：高程（x1）、坡度（x2）、坡向（x3）、平面曲率（x4）、剖面曲率（x5）、地形濕度指數(shù)（x6），3個環(huán)境因子：歸一化植被指數(shù)（NDVI，x7），土地利用類型（x8）、土壤類型（x9）的基礎(chǔ)上進(jìn)行剔除。

輔助變量數(shù)據(jù)來源：（1）海南島1︰50 000等高線矢量圖；（2）遙感影像：2013年10月26日和2013年12月6日的海南島Landsat8 OLI_TIRS四景衛(wèi)星影像，數(shù)據(jù)來源于中國科學(xué)院計算機(jī)網(wǎng)絡(luò)信息中心國際科學(xué)數(shù)據(jù)鏡像網(wǎng)站（http://www.gscloud.cn）；（3）海南島1︰200 000土壤圖；（4）2010年1︰100 000遙感解譯土地利用圖，數(shù)據(jù)源自中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心。

1.4 模型方法

本研究采用OK、RK和RF對SOM含量的空間分布特征進(jìn)行預(yù)測。其中，RF的基本原理參見文獻(xiàn)［22］。樹節(jié)點預(yù)選的變量個數(shù)（mtry）和隨機(jī)森林中樹的個數(shù)（ntree）是RF的兩個重要參數(shù)。本文通過對不同個數(shù)的輔助變量逐次計算確定出最優(yōu)的mtry值［12］。設(shè)定ntree為500、1 000、1 500、2 000，結(jié)合mtry值，選擇使模型袋外誤差（OOB）最小的參數(shù)組合用于最終預(yù)測。

1.5 數(shù)據(jù)處理與分析

（1）數(shù)字高程模型的生成。采用1︰50 000等高線矢量圖在ArcGIS9.3中生成10 m分辨率的數(shù)字高程模型柵格影像。

（2）啞變量賦值。由于土壤類型和土地利用類型為類別變量，不能直接用于回歸分析，因此采用啞變量賦值方法［23］對其賦值，賦值結(jié)果如下：（x91=1、x92=0、x93=0、x94=0、x95=0），（x91=0、x92=1、x93=0、x94=0、x95=0），（x91=0、x92=0、x93=1、x94=0、x95=0），（x91=0、x92=0、x93=0、x94=1、x95=0），（x91=0、x92=0、x93=0、x94=0、x95=1），（x91=0、x92=0、x93=0、x94=0、x95=0）分別代表磚紅壤、水稻土、赤紅壤、燥紅土、火山灰土和紫色土。（x81=1、x82=0），（x81=0、x82=1），（x81=0、x82=0）分別代表水田、旱地、園地。

（3）異常值處理與驗證集的選取。采用閾值法將區(qū)間（平均值±3×標(biāo)準(zhǔn)差）外的數(shù)值視作異常值并剔除，并參見文獻(xiàn)［2］，將風(fēng)沙土、石灰（巖）土樣點剔除，剩余有效樣點160個，從有效樣點中隨機(jī)選取128個樣點作為訓(xùn)練集，用于空間數(shù)據(jù)分析，其余32個樣點用于驗證。全部保留樣點土壤類型覆蓋海南島91.24%的陸域面積。樣點分布見圖1。

（4）模型精度驗證。采用驗證集的平均預(yù)測誤差（ME）、均方根預(yù)測誤差（RMSE）和決定系數(shù)（R2）來評價預(yù)測的準(zhǔn)確性。

（5）分析工具。地形因子的提取以及SOM含量空間制圖在ArcGIS9.3軟件中完成，遙感影像鑲嵌和NDVI的提取在ENVI5.1軟件中完成，地統(tǒng)計學(xué)分析和隨機(jī)森林規(guī)則分別在R3.3.1的gstat包和randomForest包中完成，統(tǒng)計分析使用SAS9.2軟件和Microsoft Excel 2010工具。

圖1 研究區(qū)采樣點分布圖Fig. 1 Distribution of sampling sites in the study area

2 結(jié)果與討論

2.1 不同土層SOM含量的描述性統(tǒng)計特征

由表1可知，隨著土層深度的增加，SOM含量的最小值、最大值逐漸減小，最小值由4.52減小至1.73 g kg-1，最大值由45.83減小至21.53 g kg-1。從表層（0～5 cm）到深層（40～60 cm），SOM含量的均值分別為19.67、15.89、10.30、8.07 g kg-1。按全國第二次土壤普查養(yǎng)分分級標(biāo)準(zhǔn)，除40～60 cm土層SOM含量的均值處于五級水平外，其他三個土層均處于四級水平。四個土層SOM含量的標(biāo)準(zhǔn)誤差值均較小，表明各個土層數(shù)據(jù)波動不是很大，結(jié)果較可靠。從變異程度看，四個土層SOM含量的變異系數(shù)在47.16%～56.56%，均屬于中等變異。

表1 SOM含量基本統(tǒng)計特征Table 1 Statistical characters of SOM contents

2.2 基于OK的SOM含量空間分布

對訓(xùn)練集數(shù)據(jù)進(jìn)行Kolmogorov-Smirnov檢驗發(fā)現(xiàn)，四個土層原始SOM含量數(shù)據(jù)p均小于0.01，均不符合正態(tài)分布。對其取算數(shù)平方根后，0～5、0～20、20～40和40～60 cm，p＞0.05，滿足正態(tài)分布。

由圖2可知，0～5 cm土層SOM含量共分為6個等級，按占全島面積的比例進(jìn)行排序如下：20～25 （32.63%）＞15～20 （25.83%）＞10～15 （19.36%）＞25～30 （12.99%）＞6～10（8.74%）＞30～40 g kg-1（0.45%），SOM含量＞10 g kg-1的等級占全島面積的91.26%。其中高值（25～30、30～40 g kg-1）主要分布在島內(nèi)東北的瓊山南部、文昌西南部、定安東部、瓊海北部地區(qū)和島中部的瓊中地區(qū)，低值（6～10 g kg-1）分布在島西部和西南部的昌江、東方沿海一帶。整體而言，該島SOM含量呈現(xiàn)出西部地區(qū)低于東部地區(qū)的空間分布趨勢，這與前人［19-20］的研究結(jié)果相一致。0～20 cm土層SOM含量的空間變異規(guī)律與表層相似，但SOM含量等級數(shù)減少，＞10 g kg-1等級的土壤面積占全島的84.53%，其中10～15 g kg-1和15～20 g kg-1等級所占比例相當(dāng)，二者共占全島面積的70.98%。20～40 cm土層＞10 g kg-1等級的土壤面積占全島的44.54%。而40～60 cm土層SOM含量＞10 g kg-1等級的土壤面積僅占全島的2.85%，該土層主要分布在6～10 g kg-1范圍內(nèi)，占全島面積的72.12%。

2.3 基于RK的SOM含量空間分布

本文采用逐步線性回歸方程擬合SOM含量與輔助變量之間的趨勢，然后用簡單克里格對剔除趨勢后的殘差進(jìn)行估計，將二者結(jié)果疊加（即RK）作為最終預(yù)測結(jié)果。

由表2可知，四個土層擬合方程對總方差的解釋率分別為：26.83%、17.12%、18.68%、18.04%，其回歸系數(shù)相對較小，這是因為本研究區(qū)面積較大，SOM與輔助變量之間的關(guān)系較復(fù)雜，而實際進(jìn)入方程的變量較少，故解釋率相對較低［17］。如果考慮氣溫、降水等因素，可能會獲得更高的解釋率。根據(jù)李燕麗等［24］的統(tǒng)計，逐步線性回歸模型對土壤屬性變異的解釋率在15%～82%，本研究結(jié)果在該范圍內(nèi)。且從概率來看，四個土層擬合模型均是極顯著的，表明方程能夠很好地描述SOM含量的變異規(guī)律。NDVI和土壤類型變量在四個土層中均進(jìn)入方程，表明這兩個變量是影響SOM含量的重要因素。有研究表明［1］，NDVI和SOM含量呈極顯著正相關(guān)關(guān)系，即該指數(shù)越大，SOM含量越高。海南島西部NDVI值低于東部，與西部較東部干旱，植被生長茂密程度低于東部相一致。海南島土壤類型繁多，研究表明［25］，在熱帶土壤中，不同的土壤類型，SOM含量差別較大。曾迪等［26］的研究也表明SOM含量的空間變異受土壤類型的影響。從表2中還可以看出，土地利用類型變量僅在0～5 cm和0～20 cm土層進(jìn)入方程，表明該變量對耕層的影響較大，而對深層影響較小。這是因為該區(qū)土地利用類型以水田、旱地和園地為主，受人為影響較大，隨著土層加深，施肥、耕作等人為擾動減小，所以對SOM含量變異的影響也減小。

圖2 OK預(yù)測SOM含量空間分布圖Fig. 2 OK-based SOM content spatial distribution map

表2 各土層逐步線性回歸方程擬合結(jié)果Table 2 Fitting of the soil layers with the stepwise linear regression equation

由圖3可知，應(yīng)用RK預(yù)測SOM含量的空間變異規(guī)律與OK預(yù)測所得結(jié)果略有不同，呈現(xiàn)出西南、北部、東北高，西部、東南沿海地區(qū)低的空間分布趨勢。RK對細(xì)節(jié)刻畫地更加清晰，這是引入了輔助變量的緣故。該法對SOM含量劃分等級相對于OK增多，各個等級所占的面積有所變化，高低值分布略有不同。0～5 cm土層SOM含量各等級占全島面積的比例排序如下：15～20 (30.28%)＞10～15 (22.40%)＞20～25 (21.83%)＞25～30(10.18%)＞6～10 (6.75%)＞30～40 (6.46%)＞高于40 (1.18%)＞低于6 g kg-1(0.84%)，與OK相比，該法增加了高于40 g kg-1和低于6 g kg-1兩個等級，對SOM含量的預(yù)測范圍更廣。高值（25～30、30～40 g kg-1）較OK預(yù)測區(qū)域有所增加，主要分布在瓊山、儋州北部、定安南部、瓊中西南部、五指山東北部、保亭西部、三亞東北部、樂東西北部地區(qū)，低值（6～10 g kg-1）與OK法所得結(jié)果相似，主要分布在島西部地區(qū)。這主要是引入輔助變量土壤類型所致。以下三個土層與0～5 cm分布規(guī)律相似，但等級數(shù)減少，根據(jù)表1統(tǒng)計，隨著土層深度的增加，SOM含量減少。0～20 cm土層SOM含量＞10 g kg-1等級的土壤面積占全島的86.26%，與OK該等級占全島面積的比例相當(dāng)，而10～15 g kg-1和15～20 g kg-1較OK有所減少，分別減少5.56%和3.02%，兩個等級共占全島面積的62.40%；20～40 cm土層＞10 g kg-1等級的土壤面積占全島比例較OK增加10.55%，與OK相比增加了15～20、20～25、＞25 g kg-1三個等級；40～60 cm土層主要分布在6～10 g kg-1（54.35%），與OK法所得結(jié)果相似。

圖3 RK預(yù)測SOM含量空間分布圖Fig. 3 RK-based SOM content spatial distribution map

2.4 基于 RF的SOM含量空間分布

由圖4可知，RF預(yù)測的SOM含量空間分布圖與RK預(yù)測圖相似，呈現(xiàn)出西南、東北高，西部、東南沿海地區(qū)低的空間分布趨勢，因為二者均引入了輔助變量，能夠更好地描述SOM含量空間變異的細(xì)節(jié)信息。對SOM含量等級的劃分，該法在0～5 cm和0～20 cm土層與OK相同，40～60 cm土層與RK相同，20～40 cm土層介于兩者之間。RF四個土層＞10 g kg-1等級的土壤面積占全島比例在三種方法中均最高，分別為99.79%、98.08%、60.32%、34.80%。0～5 cm土層SOM含量主要在15～25 g kg-1，占全島面積的74.6%；0～20 cm土層SOM含量主要在10～20 g kg-1（81.72%），與前兩種方法所得結(jié)果一致；20～40 cm土層SOM含量主要在6～15 g kg-1（84.2%），與RK法所得結(jié)果相同；40～60 cm土層SOM含量主要在6～15 g kg-1（89.15%），與前兩種方法相比，10～15 g kg-1等級所占比例明顯增加。

2.5 不同預(yù)測模型的精度評價

由表3可知，0～5 cm土層，三種預(yù)測方法中，從R2來看，RF ＞ RK ＞ OK，表明RF對SOM含量空間變異的解釋能力最強(qiáng)，其次為RK法，OK法最弱；從RMSE來看，RF ＜ RK ＜ OK,表明RF法預(yù)測精度最高，RK法居中，而OK法預(yù)測精度最低；OK法的ME更接近于0，其次為RK法，而RF法的無偏估計相對較差；從以上三個參數(shù)綜合來看，該土層的最優(yōu)擬合模型為RF。

圖4 RF預(yù)測SOM含量空間分布圖Fig. 4 RF-based SOM content spatial distribution map

以此方法分析其他土層，確定0～20、20～40和40～60 cm三個土層的最優(yōu)擬合模型，分別為RF法、RF法和OK法。

表3 各土層SOM含量預(yù)測模型精度Table 3 Accuracy of the prediction of SOM contents relative to soil layer

雖然RF對該區(qū)SOM含量空間變異的描述在0～5、0～20和20～40 cm土層表現(xiàn)較好。但由表3可知，驗證集中決定系數(shù)（R2）在0.19～0.37之間，較前人的研究結(jié)果［17-18,27］偏低。這是因為以往的研究采樣密度在0.03～2.52個 km-2，是本研究區(qū)采樣密度（0.004個km-2）的8倍～610倍，采樣密度小可能是造成本研究區(qū)RF預(yù)測精度偏低的原因之一。此外，本研究選取的輔助因子相對較少，有研究表明［19,28］，氣溫、降水和人為因素（如施肥等）是影響SOM含量空間變異的重要因素，如果增加上述因子可能會提高模型的解釋率。

OK在本研究中僅適合40～60 cm土層，且從預(yù)測圖（圖2）來看，該法圖斑較大，僅能夠預(yù)測出SOM含量的整體空間分布規(guī)律，對SOM含量空間變異的細(xì)節(jié)信息描述地不夠詳細(xì)。這是因為本研究區(qū)地形復(fù)雜，而有研究表明，OK通常適用于土壤屬性變化較為均勻的區(qū)域，對地形復(fù)雜、土壤屬性變化較為強(qiáng)烈的區(qū)域，其制圖精度不太理想［29］，40～60 cm土層受地形因子的影響相對表層要小，故OK法的制圖精度相對表層和其他方法均有所提高；此外，本研究采樣密度較小，對OK的制圖精度也會有一定的影響。

RK在本研究區(qū)的應(yīng)用效果不是很理想，決定系數(shù)R2與前人［30］的研究相比較小，這可能是因為剔除趨勢的模型選擇不合適，在本研究中采用的是逐步線性回歸方程來剔除趨勢，而土壤屬性與輔助變量間的關(guān)系常常是復(fù)雜的非線性關(guān)系，且因子之間可能存在交互作用，逐步線性回歸模型難以描述上述信息。以后的研究中可以考慮采用其他回歸模型（如回歸樹或RF等）來剔除土壤屬性與輔助變量間的趨勢。

從研究結(jié)果中還可以看出，不同土層所選擇的最優(yōu)模型有所不同，這可能是不同土層的影響因子略有不同，不同的模型根據(jù)自身機(jī)理，對于土壤屬性與輔助變量間的關(guān)系的解釋程度不同，因而導(dǎo)致最終的預(yù)測結(jié)果不同。

3 結(jié) 論

本研究在較大區(qū)域尺度地形復(fù)雜地區(qū)，對海南島128個樣點數(shù)據(jù)，采用OK、RK和RF結(jié)合高程、歸一化植被指數(shù)、土地利用類型和土壤類型等輔助變量對該區(qū)0～5、0～20、20～40、40～60 cm四個土層SOM含量的空間分布規(guī)律進(jìn)行預(yù)測，并以32個驗證點進(jìn)行驗證。主要結(jié)論如下：0～5、0～20和20～40 cm土層的最優(yōu)擬合模型均為RF，而40～60 cm土層的最優(yōu)擬合模型為OK，RK和RF相對于OK對SOM含量的局部變異信息描述地更加詳細(xì)。四個土層SOM含量的均值分別為19.67、15.89、10.30、8.07 g kg-1，隨著土層深度的增加，SOM含量逐漸減小。從空間預(yù)測圖來看，四個土層SOM含量均呈現(xiàn)出西南、東北高，西部、東南沿海地區(qū)低的空間分布趨勢。