申哲，張認(rèn)連，龍懷玉，王轉(zhuǎn)，朱國(guó)龍，石乾雄，喻科凡，徐愛(ài)國(guó)

（1中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所，北京 100081；2北京大學(xué)地球與空間科學(xué)學(xué)院，北京100871）

0 引言

【研究意義】土壤質(zhì)地是十分穩(wěn)定的土壤自然屬性，也是最基礎(chǔ)的土壤物理性質(zhì)之一，它影響土壤的通氣透水性、保水保肥性[1]，決定著諸多其他土壤物理化學(xué)行為。黃土高原地區(qū)是我國(guó)主要的農(nóng)業(yè)生產(chǎn)區(qū)，也是水土流失嚴(yán)重、生態(tài)脆弱區(qū)[2]，因此，利用空間預(yù)測(cè)方法實(shí)現(xiàn)對(duì)黃土區(qū)土壤顆粒組成空間分布的精準(zhǔn)預(yù)測(cè)，了解其空間分布規(guī)律，有助于研究黃土區(qū)域土壤侵蝕和土壤化學(xué)元素行為特征，為區(qū)域農(nóng)田施肥、灌水的精準(zhǔn)化管理和土壤質(zhì)量評(píng)價(jià)提供科學(xué)依據(jù)?！厩叭搜芯窟M(jìn)展】地統(tǒng)計(jì)學(xué)方法是預(yù)測(cè)土壤屬性空間分布最常見(jiàn)、最成熟的方法，以克里格插值最具代表性。傳統(tǒng)的普通克里格（ordinary Kriging，OK）適用于較均一、土壤屬性變化不強(qiáng)烈的環(huán)境[3-4]，在小尺度和均質(zhì)景觀區(qū)域應(yīng)用較多[5-7]。黃土高原地形復(fù)雜，土壤屬性空間變異程度大，其空間分布受地形等因素影響顯著[8-9]，結(jié)合輔助變量的混合插值模型和常規(guī)統(tǒng)計(jì)模型在黃土區(qū)取得了較好的空間預(yù)測(cè)效果。文雯等[10]利用基于土地利用類(lèi)型修正的OK對(duì)黃土丘陵區(qū)土壤有機(jī)碳進(jìn)行了空間插值。連綱等[11]以地形因子和土地利用方式為輔助變量，采用多元線性回歸（multiple linear regression，MLR）和回歸克里格（regression Kriging，RK）對(duì)黃土高原縣域土壤養(yǎng)分進(jìn)行空間預(yù)測(cè)，結(jié)果表明RK有效地減小了預(yù)測(cè)殘差，精度高于MLR。LIU等[12]利用空間狀態(tài)方程和 MLR結(jié)合地形、氣候等要素對(duì)整個(gè)黃土高原的土壤有機(jī)碳進(jìn)行空間模擬，結(jié)果表明空間狀態(tài)方程取得了更好的預(yù)測(cè)效果。除加入輔助變量提高預(yù)測(cè)精度以外，經(jīng)驗(yàn)貝葉斯克里格（empirical Bayesian Kriging，EBK）也被證明是一種可對(duì)非均質(zhì)景觀區(qū)域進(jìn)行空間預(yù)測(cè)的插值方法[13]，但目前該方法在黃土區(qū)應(yīng)用較少。地統(tǒng)計(jì)和傳統(tǒng)統(tǒng)計(jì)模型均難以擬合土壤屬性和環(huán)境變量之間的非線性關(guān)系，且面臨著處理類(lèi)別型變量的困難。為克服這些局限，分類(lèi)與回歸樹(shù)（classification and regression tree，CART）、神經(jīng)網(wǎng)絡(luò)（artificial neural network，ANN）等機(jī)器學(xué)習(xí)方法被引入到土壤屬性空間預(yù)測(cè)中，研究表明這些模型能成功預(yù)測(cè)類(lèi)別型變量、提高土壤屬性預(yù)測(cè)的精度[14-15]，但兩者均易過(guò)度擬合[16-17]，且神經(jīng)網(wǎng)絡(luò)過(guò)程繁瑣。隨機(jī)森林（random forest，RF）模型在分類(lèi)與回歸樹(shù)模型基礎(chǔ)上發(fā)展而來(lái)，具有計(jì)算簡(jiǎn)單、不易過(guò)度擬合、反映輔助變量重要程度等特點(diǎn)[18]，目前在平原和山區(qū)、流域等母質(zhì)多樣化地區(qū)應(yīng)用較多。LIU等[19]基于MODIS獲取的地表溫度資料，用RF法完成了江蘇平原土壤顆粒組成、有機(jī)質(zhì)和PH的空間制圖。郭澎濤等[20]利用RF結(jié)合多源環(huán)境變量（成土母質(zhì)、平均降雨量、平均氣溫和歸一化植被指數(shù)）對(duì)橡膠園土壤全氮含量進(jìn)行空間預(yù)測(cè)，并與MLR、CART比較，結(jié)果表明RF模型表現(xiàn)最優(yōu)。MAREIKE等[21]在厄爾多瓜安第斯山脈對(duì)土壤顆粒組成空間分布的研究和GERALD 等[22]在布基納法索西南部小流域?qū)ν寥李w粒組成、有機(jī)碳空間預(yù)測(cè)的研究均表明，RF法優(yōu)于其他機(jī)器學(xué)習(xí)方法?！颈狙芯壳腥朦c(diǎn)】就研究對(duì)象而言，前人對(duì)黃土區(qū)土壤屬性空間預(yù)測(cè)的研究集中在土壤養(yǎng)分指標(biāo)，對(duì)土壤顆粒組成等物理屬性研究相對(duì)較少。與土壤養(yǎng)分指標(biāo)相比，土壤顆粒組成受人類(lèi)活動(dòng)等隨機(jī)因素的影響較小，受母質(zhì)、地形、氣候等結(jié)構(gòu)性因素影響較大，空間分布規(guī)律性更強(qiáng)。且土壤顆粒組成屬于成分?jǐn)?shù)據(jù)，預(yù)測(cè)結(jié)果要滿(mǎn)足非負(fù)、定和為1的要求。就空間預(yù)測(cè)方法而言，目前黃土區(qū)已有的研究中對(duì)較新近的EBK使用較少，且缺乏地統(tǒng)計(jì)學(xué)方法與機(jī)器學(xué)習(xí)方法的對(duì)比研究。就輔助變量的選擇而言，平原地區(qū)多選擇高光譜數(shù)據(jù)[23-24]、地表溫度[19]等作為輔助變量進(jìn)行預(yù)測(cè)；山區(qū)、流域多為母質(zhì)多樣化地區(qū)，以地形因子結(jié)合成土母質(zhì)[20,25]為輔助變量能取得理想的預(yù)測(cè)結(jié)果。但在地形復(fù)雜、母質(zhì)較為單一的區(qū)域，成土母質(zhì)無(wú)法作為輔助變量參與空間預(yù)測(cè)，僅使用地形因子預(yù)測(cè)可能精度并不高，需考慮加入其他環(huán)境變量提高預(yù)測(cè)精度。因此，對(duì)地形復(fù)雜、單一母質(zhì)地區(qū)的土壤屬性進(jìn)行空間預(yù)測(cè)時(shí)，在方法的選擇、輔助變量的選擇方面仍需進(jìn)一步研究?！緮M解決的關(guān)鍵問(wèn)題】本文選擇黃土、黃土狀母質(zhì)的寧夏自治區(qū)海原縣為研究區(qū)，選用3種方法（SLR-EBK、SLR-RK、SLR-RF）結(jié)合地形、土壤類(lèi)型、歸一化植被指數(shù)等輔助變量進(jìn)行表層土壤顆粒組成的空間制圖，并進(jìn)行精度對(duì)比選出最優(yōu)模型，以期為地形復(fù)雜的黃土母質(zhì)區(qū)土壤顆粒組成的空間預(yù)測(cè)提供依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

海原縣地處黃土高原西北部（36°6′—37°4′N(xiāo)、105°9′—106°10′E），面積 689 900 hm2。該區(qū)屬大陸季風(fēng)氣候，年均氣溫7℃，≥10℃年積溫2 398℃，年平均降水量 363 mm，降水量由南向北遞減?？h內(nèi)地勢(shì)由西南向東北傾斜（圖1），地貌以黃土丘陵為主，南部月亮山、中部南華山、西部西華山、北部油井山呈孤立零星分布，被黃土丘陵包圍，東北角與同心縣交界處為狹長(zhǎng)的河谷沖積平原，面積較小，為提灌區(qū)。研究區(qū)內(nèi)絕大部分成土母質(zhì)為黃土、黃土狀母質(zhì)，僅在南華山、西華山零星分布有巖石風(fēng)化物。土壤類(lèi)型包括灰褐土、黑壚土、黃綿土、灰鈣土、新積土、紅黏土和粗骨土7種，各土壤類(lèi)型面積占比見(jiàn)表1。其中黃綿土、黑壚土、灰鈣土 3個(gè)土類(lèi)面積最大，共占全縣面積的81.51%。以鹽池鄉(xiāng)-李旺一線為界，此線以北主要為灰鈣土，以南黃綿土和黑壚土插花分布，南華山、西華山發(fā)育有灰褐土，丘陵間低地及沖積平原有新積土分布。海原縣以雨養(yǎng)農(nóng)業(yè)為主，土地利用類(lèi)型主要為旱地、草地，種植作物主要為玉米、蕎麥、馬鈴薯。

表1 各土壤類(lèi)型面積占比及樣點(diǎn)統(tǒng)計(jì)Table 1 The area proportion and sample points of each soil type

圖1 研究區(qū)位置及樣點(diǎn)分布圖Fig. 1 The map of the study position and sampling sites

1.2 樣品采集與分析

采用網(wǎng)格法布點(diǎn)，在此基礎(chǔ)上綜合考慮地形地貌、土壤類(lèi)型等因素，在南華山、西華山等地形復(fù)雜地區(qū)加密樣點(diǎn)，共選取具有代表性的樣點(diǎn)124個(gè)（圖1）。以鹽池鄉(xiāng)-南華山-九彩鄉(xiāng)一線為界，南部面積占32%，采樣點(diǎn)共53個(gè)，占全部采樣點(diǎn)的43%，平均每4 300 hm2一個(gè)樣點(diǎn)；北部面積占68%，采樣點(diǎn)共71個(gè)，平均每6 600 hm2一個(gè)樣點(diǎn)。研究區(qū)全部樣點(diǎn)平均間隔7—8 km。東南部邊緣因地勢(shì)陡峭無(wú)路，取樣較困難，故3個(gè)樣點(diǎn)取在海原縣外3—5 km范圍內(nèi)。各土壤類(lèi)型的采樣點(diǎn)個(gè)數(shù)見(jiàn)表1。土壤樣品采集時(shí)間為2018年9月、2019年4月，采樣深度為0—20 cm，記錄樣點(diǎn)經(jīng)緯度、地形地貌、土壤類(lèi)型、土地利用類(lèi)型等信息。土壤顆粒組成采用濕篩-吸管法測(cè)定[26]，粒級(jí)劃分采用國(guó)際制（砂粒2—0.02 mm，粉粒0.02—0.002 mm，黏粒＜0.002 mm）。

1.3 輔助變量來(lái)源及處理

海原縣內(nèi)地形復(fù)雜，地勢(shì)起伏較大，海拔落差大于1 600 m，地貌類(lèi)型主要包括山地、黃土丘陵、平原。因地形地貌不同，各區(qū)域溫度、降水量、植被差異較大，直接影響成土過(guò)程，進(jìn)而影響土壤顆粒組成的空間分布。土壤類(lèi)型是五大成土因素綜合作用的產(chǎn)物，一定程度上包含土壤質(zhì)地的信息；土壤類(lèi)型也含有母質(zhì)信息，而本研究區(qū)母質(zhì)相對(duì)單一，因此，其他因素對(duì)質(zhì)地的影響會(huì)更大；另外，許多研究表明，土壤類(lèi)型與土壤有機(jī)質(zhì)、容重有一定相關(guān)性[27-28]。溫度、降水量對(duì)土壤顆粒組成有影響，由于缺乏研究區(qū)較精細(xì)的溫度、降水?dāng)?shù)據(jù)，而植被覆蓋情況直接受溫度、降水影響，且研究區(qū)基本為雨養(yǎng)農(nóng)業(yè)，故將反映植被覆蓋程度的歸一化植被指數(shù)作為輔助變量反映氣候因素。同時(shí)地形、土壤類(lèi)型在一定程度上也包含氣象信息。與土壤養(yǎng)分相比，土壤顆粒組成穩(wěn)定，受人為耕作影響較小，研究區(qū)土地利用類(lèi)型主要為旱地和草地，受灌水耕作影響很小，因此本研究區(qū)黏土礦物分解更加緩慢。綜上所述，本研究選取的輔助變量為地形因子、土壤類(lèi)型、歸一化植被指數(shù)。

1.3.1 輔助變量來(lái)源

（1）寧夏自治區(qū)1∶50 000等高線矢量圖；（2）寧夏自治區(qū)1∶200 000土壤圖；（3）本文使用處于植物生長(zhǎng)季（8—9 月）、且云量覆蓋度少（小于5%）的衛(wèi)星影像計(jì)算植被指數(shù)。2017年8月28日和9月6日的Landsat8 OLI_TIRS 四景衛(wèi)星影像覆蓋海原縣全境，影像均來(lái)源于地理空間數(shù)據(jù)云網(wǎng)站（http://www.gscloud.cn）。由于缺乏研究區(qū)2018年5月以后的衛(wèi)星影像，故采用前一年相應(yīng)月份的影像代替。

1.3.2 輔助變量的提取及處理

（1）用 1∶50 000等高線矢量圖在 ArcGIS10.2中生成10 m分辨率的數(shù)字高程模型（digital elevation model，DEM）柵格影像，從中提取 5個(gè)地形因子：高程（elevation，Ele）、坡度（slope，Slo）、平面曲率（horizontal curvature，HC）、剖面曲率（profile curvature，PC）、地形濕度指數(shù)（topographic wetness index，TWI）。使用SAGA GIS軟件基于DEM提取風(fēng)力作用指數(shù)（wind exposition index，WEI）[29]，該指數(shù)是各個(gè)風(fēng)向?qū)^(qū)域的風(fēng)影響（wind effect）[29]的平均值。WEI越大，表明該區(qū)域越容易受風(fēng)力影響。（2）土壤類(lèi)型（soil types，ST）為類(lèi)別變量，本文采用“算術(shù)平均值變換”[30]，以不同土壤類(lèi)型下定量因變量的算術(shù)平均值代替該土壤類(lèi)型，即將土壤類(lèi)型變?yōu)橐粋€(gè)具有7個(gè)數(shù)值的定量變量。由于有砂粒、粉粒、黏粒3個(gè)因變量，則土壤類(lèi)型共轉(zhuǎn)化為3組定量變量，分別用于構(gòu)建砂粒、粉粒、黏粒的預(yù)測(cè)模型。算術(shù)平均值變換用動(dòng)態(tài)思維（類(lèi)別自變量與定量因變量的關(guān)系）建立起自變量各水平與定量結(jié)果變量之間數(shù)量關(guān)系，比傳統(tǒng)的啞變量變換更具合理性[30]；此外，啞變量變換將具有n個(gè)水平的類(lèi)別變量變?yōu)椋╪-1）個(gè)二值變量，這些二值變量非相互獨(dú)立[31]，需全部進(jìn)入回歸模型中，造成模型中無(wú)關(guān)變量太多，擬合精度不高，而算術(shù)均值變換將類(lèi)別變量變?yōu)榫哂衝個(gè)數(shù)值的定量變量，提高了回歸擬合精度[30]。（3）在ENVI5.1中對(duì)遙感影像進(jìn)行鑲嵌、計(jì)算，提取歸一化植被指數(shù)（normalized difference vegetation index，NDVI）。（4）為消除變量之間的量綱影響，對(duì)變量進(jìn)行歸一化處理。

1.4 預(yù)測(cè)方法

1.4.1 數(shù)據(jù)轉(zhuǎn)換方法 為了使數(shù)據(jù)滿(mǎn)足正態(tài)分布，且使預(yù)測(cè)結(jié)果滿(mǎn)足非負(fù)、定和為1的要求，本文首先對(duì)土壤顆粒組成進(jìn)行對(duì)稱(chēng)對(duì)數(shù)比（symmetry log-ratio，SLR）轉(zhuǎn)換[32]，轉(zhuǎn)換公式為：

轉(zhuǎn)回公式為：

式中，Zij代表第i個(gè)樣點(diǎn)第j種粒級(jí)的相對(duì)含量（%）；代表第i個(gè)樣點(diǎn)第j種粒級(jí)含量的轉(zhuǎn)換值；D表示成分?jǐn)?shù)據(jù)的維度，本文中D=3；δj為常數(shù)，取研究區(qū)第j種顆粒除0外最小含量的一半。

1.4.2 經(jīng)驗(yàn)貝葉斯克里格（EBK） EBK在OK的基礎(chǔ)上發(fā)展而來(lái)，可以通過(guò)構(gòu)造子集、建立局部模型對(duì)非平穩(wěn)數(shù)據(jù)進(jìn)行空間插值[33]。它克服了OK的兩大局限：第一，EBK則通過(guò)構(gòu)造子集，通過(guò)評(píng)估較小子集上的半方差函數(shù)，在研究區(qū)的不同區(qū)域中得到不同的半方差函數(shù)模型，克服了OK假設(shè)單一半方差函數(shù)可準(zhǔn)確表示所有位置數(shù)據(jù)的空間結(jié)構(gòu)的局限。第二，在半方差函數(shù)的參數(shù)評(píng)估方面，OK只依賴(lài)于單一半方差函數(shù)，忽略了模型的不確定性[34]，因此在構(gòu)建半方差函數(shù)的過(guò)程中，需對(duì)每個(gè)參數(shù)進(jìn)行手動(dòng)反復(fù)調(diào)試，以使該模型更接近于真實(shí)的半方差函數(shù)。而EBK在每個(gè)子集中，自動(dòng)評(píng)估半方差函數(shù)，并以此半方差函數(shù)來(lái)模擬子集中的新數(shù)據(jù)值。然后使用這些模擬的數(shù)據(jù)值來(lái)評(píng)估子集的新半方差函數(shù)。此模擬和評(píng)估過(guò)程重復(fù)多次，并且在每個(gè)子集中生成多個(gè)模擬的半方差函數(shù)，最后將這些模擬混合在一起以生成最終的預(yù)測(cè)圖。本研究中，訓(xùn)練集砂粒、粉粒、黏粒含量經(jīng)SLR轉(zhuǎn)換后的全局Moran’s I（莫蘭指數(shù)）分別為0.293、0.260、0.192，且P值均小于0.05，表明經(jīng)SLR轉(zhuǎn)換的土壤顆粒組成在研究區(qū)域內(nèi)有一定的空間自相關(guān)性。ArcGIS趨勢(shì)分析發(fā)現(xiàn)，轉(zhuǎn)換后的砂粒、粉粒含量均存在二階趨勢(shì)，不滿(mǎn)足空間平穩(wěn)假設(shè)。Voronoi圖也顯示兩者的局部空間規(guī)律變化較大，因此使用EBK進(jìn)行空間插值。具體步驟為：（1）將數(shù)據(jù)分為多個(gè)特定大小的重疊子集（通過(guò)逐次嘗試，最終本文每個(gè)子集設(shè)置為70個(gè)點(diǎn)，子集重疊系數(shù)設(shè)置為3，即每個(gè)點(diǎn)可落入3個(gè)子集中）；（2）對(duì)子集中的數(shù)據(jù)進(jìn)行半方差函數(shù)擬合；（3）使用該半方差函數(shù)，在每個(gè)輸入位置生成新值；（4）使用新的模擬數(shù)據(jù)生成新的半方差函數(shù)；（5）將步驟（3）和（4）重復(fù)100次。此過(guò)程將為每個(gè)子集生成100個(gè)半方差函數(shù)，每個(gè)都是子集真實(shí)半方差函數(shù)的估計(jì)。將其繪制在一起生成按密度著色的半方差函數(shù)分布；（6）對(duì)于空間內(nèi)每個(gè)位置，都使用唯一的半方差函數(shù)分布生成預(yù)測(cè)，該分布是通過(guò)周?chē)蛹陌敕讲詈瘮?shù)分布加權(quán)綜合計(jì)算得出的；子集距離預(yù)測(cè)位置越近，給定的權(quán)重就越高。這就保證了每個(gè)預(yù)測(cè)值使用的是臨近數(shù)據(jù)所定義的模型，而不是遠(yuǎn)處數(shù)據(jù)所定義的模型[13]。最后，將局部預(yù)測(cè)結(jié)果混合生成全局預(yù)測(cè)結(jié)果。

1.4.3 回歸克里格（RK） RK是應(yīng)用多元線性回歸擬合土壤屬性與輔助變量間的關(guān)系，然后對(duì)回歸殘差應(yīng)用克里格法插值，最后將回歸預(yù)測(cè)結(jié)果與殘差插值結(jié)果結(jié)合起來(lái)的一種空間預(yù)測(cè)方法。首先用皮爾森相關(guān)分析對(duì)輔助變量進(jìn)行初步篩選，表2相關(guān)分析的結(jié)果表明，高程、歸一化植被指數(shù)和土壤類(lèi)型與土壤顆粒組成顯著相關(guān)。其中，高程和NDVI均與土壤砂粒含量呈負(fù)相關(guān)，與粉粒和黏粒含量呈正相關(guān)。為避免輔助變量間的共線性問(wèn)題，基于篩選出的變量，采用逐步線性回歸的方法對(duì)經(jīng) SLR轉(zhuǎn)換后的土壤顆粒組成進(jìn)行預(yù)測(cè)（F值的概率Fp＜0.05為變量進(jìn)入方程的標(biāo)準(zhǔn)，F(xiàn)p＞0.10為剔除變量的標(biāo)準(zhǔn)），并對(duì)預(yù)測(cè)殘差進(jìn)行簡(jiǎn)單克里格插值，最后將回歸預(yù)測(cè)值與殘差相加并轉(zhuǎn)回得到預(yù)測(cè)結(jié)果。

表2 土壤各粒級(jí)含量與輔助變量的相關(guān)性分析Table 2 Pearson correlation coefficients between soil particle size distribution and auxiliary variables

1.4.4 隨機(jī)森林（RF） 與其他模型一樣，隨機(jī)森林（RF）可以反映若干自變量對(duì)因變量的作用[35]。其主要過(guò)程包括：從原始訓(xùn)練集中進(jìn)行bootstrap重新抽樣（有放回的隨機(jī)抽?。?gòu)成新的訓(xùn)練集生成決策樹(shù)；同時(shí)，隨機(jī)地選擇部分變量進(jìn)行決策樹(shù)節(jié)點(diǎn)的確定；隨機(jī)地生成幾百至幾千個(gè)決策樹(shù)構(gòu)成隨機(jī)森林。對(duì)于分類(lèi)問(wèn)題，所有決策樹(shù)預(yù)測(cè)結(jié)果中得票最高的分類(lèi)結(jié)果為最終結(jié)果；對(duì)于回歸問(wèn)題，所有決策樹(shù)預(yù)測(cè)結(jié)果的平均值為最終的預(yù)測(cè)結(jié)果[18]。每次未被抽到的樣本組成袋外樣本，故RF可根據(jù)袋外誤差估計(jì)模型誤差，對(duì)于分類(lèi)問(wèn)題，袋外誤差是分類(lèi)錯(cuò)誤率；對(duì)于回歸問(wèn)題，袋外誤差通過(guò)回歸殘差計(jì)算，本文采用均方誤差（mean square eroor，MSE）。

RF不要求原始數(shù)據(jù)符合正態(tài)分布，但為了更好地同其他兩種方法做對(duì)比，且使預(yù)測(cè)結(jié)果滿(mǎn)足定和為 1的要求，本研究使用經(jīng)SLR轉(zhuǎn)換后的數(shù)據(jù)進(jìn)行建模。由于 RF不是簡(jiǎn)單的線性回歸，皮爾森相關(guān)分析篩選變量的結(jié)果不一定適合 RF模型，故用兩步法對(duì) RF輔助變量進(jìn)行篩選：對(duì)變量進(jìn)行重要性排序，初步剔除不重要變量[36]；對(duì)剩下的變量，基于隨機(jī)森林逐步選擇[37]。本文采用Breiman經(jīng)典隨機(jī)森林版本中的評(píng)價(jià)方法，通過(guò)變量值的置換計(jì)算變量重要性得分，剔除得分小于0的變量。結(jié)果表明平面曲率、剖面曲率兩個(gè)輔助變量重要性得分均小于 0，故剔除。剩下的變量按照重要性評(píng)分降序排列，以不同的變量子集所構(gòu)建森林運(yùn)行100次的MSE均值作為逐步篩選變量的標(biāo)準(zhǔn)，從最重要變量開(kāi)始，遍歷所有變量，MSE最小結(jié)束（變量篩選過(guò)程中，RF參數(shù)均為默認(rèn)參數(shù)）?；诤Y選的輔助變量對(duì) RF模型進(jìn)行參數(shù)優(yōu)化：首先以 MSE最小化為目標(biāo)，通過(guò)逐次計(jì)算確定決策樹(shù)節(jié)點(diǎn)分裂時(shí)所用變量個(gè)數(shù)（mtry）；通過(guò)觀察決策樹(shù)數(shù)量（ntree）與誤差關(guān)系圖，在保證誤差穩(wěn)定的前提下，選擇較小的ntree[38]。

本研究中隨機(jī)森林在R3.5.2中完成，相關(guān)性分析和逐步回歸使用SPSS 22軟件，經(jīng)驗(yàn)貝葉斯克里格插值、回歸殘差的簡(jiǎn)單克里格插值和空間制圖在ArcGIS 10.2中完成。

1.5 精度驗(yàn)證

從124個(gè)樣點(diǎn)中隨機(jī)選取80%的（100個(gè)）樣點(diǎn)用于空間預(yù)測(cè)，20%的（24個(gè)）樣點(diǎn)用于精度驗(yàn)證[39-40]，訓(xùn)練樣點(diǎn)和驗(yàn)證樣點(diǎn)的分布如圖1所示。本研究采用驗(yàn)證集的平均絕對(duì)誤差（MAE）、均方根預(yù)測(cè)誤差（RMSE）、平均Aitchison距離（MAD）[41]進(jìn)行模型精度評(píng)價(jià)。MAE能更好地反映預(yù)測(cè)值誤差的實(shí)際情況；RMSE對(duì)系統(tǒng)誤差和隨機(jī)誤差都很敏感，常用來(lái)評(píng)估不同模型的預(yù)測(cè)精度；MAD可以反映成分?jǐn)?shù)據(jù)預(yù)測(cè)值和觀測(cè)值之間整體的相似性。MAE、 RMSE、MAD越小，模型精度越高。

Zij代表第i個(gè)樣點(diǎn)第j種粒級(jí)含量的觀測(cè)值；代表第i個(gè)樣點(diǎn)第j種粒級(jí)含量的預(yù)測(cè)值；n表示驗(yàn)證點(diǎn)的數(shù)量；D表示成分?jǐn)?shù)據(jù)的維度，本文中D=3。

2 結(jié)果

2.1 土壤顆粒組成的描述性統(tǒng)計(jì)特征

表3為研究區(qū)采樣點(diǎn)土壤顆粒組成的描述性統(tǒng)計(jì)結(jié)果，可以看出研究區(qū)土壤砂粒含量最高（訓(xùn)練集和驗(yàn)證集砂粒平均含量分別為61.82%、63.96%），黏粒含量最低（訓(xùn)練集和驗(yàn)證集的平均黏粒含量分別為14.56%、13.36%），粉粒含量居中（訓(xùn)練集和驗(yàn)證集的平均粉粒含量分別為23.62%、22.68%）。訓(xùn)練集中3個(gè)粒級(jí)含量的范圍和標(biāo)準(zhǔn)差都比驗(yàn)證集中的大，相對(duì)應(yīng)地，訓(xùn)練集中砂粒、粉粒、黏粒含量的變異系數(shù)（14.02%、23.96%、27.41%）大于驗(yàn)證集（11.30%、21.91%、25.45%）。方差分析結(jié)果表明訓(xùn)練集與驗(yàn)證集土壤3個(gè)粒級(jí)含量差異均不顯著，表明訓(xùn)練數(shù)據(jù)和驗(yàn)證數(shù)據(jù)都具有較好的代表性。訓(xùn)練集中3個(gè)粒級(jí)含量偏度均接近于0，但經(jīng)Kolmogorov-Smirnov檢驗(yàn)，只有砂粒含量P值大于0.05，滿(mǎn)足正態(tài)分布，粉粒和黏粒含量均不符合正態(tài)分布。經(jīng)SLR轉(zhuǎn)換后，砂粒、粉粒、黏粒含量均符合正態(tài)分布，可用于空間插值和線性回歸。

2.2 多元線性回歸和RF模型的構(gòu)建

由表4可以看出，最終進(jìn)入土壤顆粒組成線性回歸預(yù)測(cè)模型的輔助變量包括高程 Ele和土壤類(lèi)型，NDVI由于與高程、土壤類(lèi)型存在多重共線性，在逐步線性回歸的過(guò)程中被剔除。從調(diào)整決定系數(shù)（AdjustedR2）來(lái)看，砂粒含量的線性回歸模型方差解釋率最高（0.339），擬合度較好；粉粒和黏粒含量的線性回歸方程方差解釋率較低（分別為 0.208，0.205），擬合度較差。表5為RF模型的擬合結(jié)果，可以看出進(jìn)入模型預(yù)測(cè)的輔助變量除了高程 Ele、土壤類(lèi)型外，還包括坡度Slo和風(fēng)力作用指數(shù)WEI，這主要是因?yàn)?RF可以擬合土壤顆粒組成和輔助變量之間的非線性關(guān)系。其中，高程對(duì)3個(gè)粒級(jí)含量的空間預(yù)測(cè)來(lái)說(shuō)，均是最重要的輔助變量，其次是土壤類(lèi)型，坡度、風(fēng)力作用指數(shù)的重要性相對(duì)較低。雖然 RF對(duì)輔助變量間的多重共線性不敏感，但考慮到NDVI對(duì)預(yù)測(cè)精度無(wú)顯著影響，為了簡(jiǎn)化模型，提高計(jì)算速度，本研究未將NDVI作為RF的輔助變量。

表3 研究區(qū)采樣點(diǎn)土壤顆粒組成基本統(tǒng)計(jì)特征Table 3 Statistical characters of soil particle size distribution in study area

表4 土壤顆粒組成逐步線性回歸方程擬合結(jié)果Table 4 Fitting of soil particle size distribution with the stepwise linear regression equation

表5 SLR-RF參數(shù)擬合結(jié)果Table 5 SLR-RF parameter fitting results

2.3 土壤顆粒組成的空間分布預(yù)測(cè)結(jié)果

利用擬合的SLR-EBK、SLR-RK和RF模型分別對(duì)研究區(qū)土壤砂粒、粉粒、黏粒含量的空間分布進(jìn)行預(yù)測(cè)，得到其空間分布圖（圖2—4）。可以看出，3種方法預(yù)測(cè)的土壤各粒級(jí)含量空間分布的趨勢(shì)大體一致。研究區(qū)砂粒含量基本呈現(xiàn)西南部低，東北部高的趨勢(shì)，其中西華山-南華山-九彩鄉(xiāng)海拔較高一帶含量較低，集中在60%以下，北部和東部邊緣含量較高，大多數(shù)為65%—70%，部分達(dá)到70%以上。粉粒和黏粒含量的空間分布趨勢(shì)與砂粒相反，呈現(xiàn)西南部高，東北部低的趨勢(shì)，西華山-南華山-九彩鄉(xiāng)一帶粉粒含量變化范圍集中在 25%—30%，黏粒含量變化范圍集中在 15%—20%，北部和東部邊緣粉粒含量低至 20%以下，黏粒含量大部分為10%—15%。土壤顆粒組成的空間分布趨勢(shì)明顯受到海拔的影響。

雖然3種方法預(yù)測(cè)的土壤顆粒組成空間分布總體趨勢(shì)比較相似，但在局部區(qū)域上不盡相同。首先，SLR-RK和SLR-RF的預(yù)測(cè)結(jié)果對(duì)細(xì)節(jié)刻畫(huà)比較清晰，而 SLR-EBK預(yù)測(cè)結(jié)果只能反映海原縣土壤顆粒組成的宏觀分布趨勢(shì)，難以體現(xiàn)局部變異信息。其次，雖然3種方法預(yù)測(cè)的砂粒、粉粒、黏粒含量的平均值相近（SLR-EBK、SLR-RK、SLR-RF 3種方法預(yù)測(cè)砂粒含量的均值分別為62.32%、62.83%、62.53%，粉粒含量的均值分別為23.34%、23.25%、23.38%，黏粒含量的均值分別為14.34%、13.91%、14.09%），但由圖2—4可以看出，SLR-EBK存在平滑效應(yīng)，與土壤采樣點(diǎn)原始數(shù)據(jù)相比，該方法預(yù)測(cè)的極大值和極小值均向均值方向靠攏，SLR-RK與SLR-RF的預(yù)測(cè)范圍更接近實(shí)測(cè)值。

圖2 SLR-EBK預(yù)測(cè)土壤顆粒組成空間分布圖Fig. 2 Spatial distribution of soil particle size distribution predicted by SLR-EBK

圖3 SLR-RK預(yù)測(cè)土壤顆粒組成空間分布圖Fig. 3 Spatial distribution of soil particle size distribution predicted by SLR-RK

圖4 SLR-RF預(yù)測(cè)土壤顆粒組成空間分布圖Fig. 4 Spatial distribution of soil particle size distribution predicted by SLR-RF

2.4 不同預(yù)測(cè)方法的精度對(duì)比

表6列出了3種方法在驗(yàn)證集的預(yù)測(cè)精度，可以看出SLR-RF法對(duì)3個(gè)粒級(jí)含量預(yù)測(cè)的MAE和RMSE均低于其他兩種方法，且SLR-RF的MAD（0.208）最低，表明該方法的預(yù)測(cè)精度更高。因此，研究區(qū)土壤顆粒組成的最優(yōu)預(yù)測(cè)模型為SLR-RF。

表6 不同空間預(yù)測(cè)方法精度對(duì)比Table 6 Accuracy comparison of different spatial prediction methods

3 討論

3.1 3種預(yù)測(cè)方法的精度

本文利用3種方法在寧夏南部黃土丘陵區(qū)對(duì)土壤顆粒組成的空間分布進(jìn)行預(yù)測(cè)，結(jié)果顯示，SLR-EBK方法預(yù)測(cè)精度較低，原因在于空間插值方法僅利用了土壤顆粒組成的空間自相關(guān)性，忽略了環(huán)境變量的影響。雖然EBK在一定程度上考慮了土壤屬性變化的不均勻性，但由于本研究區(qū)地形破碎，且南部高山區(qū)與北部黃土丘陵區(qū)氣候條件差異大，故 SLR-EBK法預(yù)測(cè)精度并不高。

SLR-RK、SLR-RF法由于引入輔助變量參與預(yù)測(cè)，提高了土壤顆粒組成的空間預(yù)測(cè)精度，相比SLR-RK，SLR-RF的預(yù)測(cè)精度更高。可以看出，雖然RK通過(guò)對(duì)回歸殘差進(jìn)行簡(jiǎn)單克里格插值進(jìn)一步提高了多元線性回歸的預(yù)測(cè)精度，但 RF通過(guò)精確捕捉土壤屬性同輔助變量之間的非線性關(guān)系，依然取得了比RK更好的預(yù)測(cè)效果。在黃土母質(zhì)區(qū)，地形是影響土壤屬性的主要因素，精確擬合地形因子與土壤屬性之間的關(guān)系是進(jìn)行土壤屬性空間預(yù)測(cè)的關(guān)鍵所在。小尺度下，氣候、植被等環(huán)境條件相對(duì)一致，地形因子與土壤屬性之間的關(guān)系較為簡(jiǎn)單，可直接利用線性回歸擬合[20]。王庫(kù)[42]以高程、坡度等為輔助變量，使用RK對(duì)福建省旗山北麓一塊面積為29 000 hm2的區(qū)域進(jìn)行土壤全氮的空間預(yù)測(cè)，結(jié)果表明，RK取得了較高的預(yù)測(cè)精度，驗(yàn)證集的決定系數(shù)達(dá)到0.682。MOORE等[43]在科羅拉多州一塊 5.4 hm2的坡面上，利用多元線性回歸結(jié)合地形因子（坡度、坡向、平面曲率、剖面曲率、地形濕度指數(shù)等）對(duì)土壤屬性（有機(jī)質(zhì)、pH、砂粒含量、粉粒含量等）進(jìn)行空間預(yù)測(cè)，取得了理想的預(yù)測(cè)結(jié)果。但在大尺度區(qū)域，隨著地形的變化，氣候、植被等也出現(xiàn)明顯變化，土壤屬性與地形因子之間往往表現(xiàn)出非線性的特征，這時(shí)運(yùn)用 RF進(jìn)行空間預(yù)測(cè)就有明顯優(yōu)勢(shì)。姜賽平[28]等以整個(gè)海南島（3.29×106hm2）為研究區(qū)域，利用RK、RF結(jié)合輔助變量預(yù)測(cè)全島的有機(jī)質(zhì)空間分布，結(jié)果表明 RF預(yù)測(cè)精度高于 RK，更適合地形復(fù)雜的大尺度地區(qū)土壤屬性的空間預(yù)測(cè)。本研究區(qū)面積達(dá)689 900 hm2，南北部地形、氣候、植被差異較大，部分地形因子（坡度、風(fēng)力作用指數(shù)）與土壤顆粒組成無(wú)顯著的線性關(guān)系，并未進(jìn)入回歸方程，但卻進(jìn)入 RF模型參與空間預(yù)測(cè)，表明RF模型可以擬合土壤屬性同環(huán)境變量之間的非線性關(guān)系，進(jìn)一步提高預(yù)測(cè)精度。除地形因子外，本研究中土壤類(lèi)型也是預(yù)測(cè)土壤顆粒組成的重要變量。土壤類(lèi)型屬于類(lèi)別變量，與多元線性回歸模型相比，RF對(duì)類(lèi)別變量具有更好的包容性[44]，可以更多地捕捉土壤顆粒組成與土壤類(lèi)型之間的關(guān)系。綜上所述，RF更適合用于預(yù)測(cè)研究區(qū)土壤顆粒組成的空間分布。

3.2 土壤顆粒組成空間分布與輔助變量的關(guān)系

地形作為五大成土因素之一，通過(guò)影響物質(zhì)與能量的再分配，引起土壤的理化性質(zhì)變化[45]。本研究中，高程是預(yù)測(cè)土壤顆粒組成最重要的地形因子。皮爾森相關(guān)分析和空間預(yù)測(cè)結(jié)果均表明，海拔高度越高，砂粒含量越低，粉粒和黏粒含量越高。這與部分學(xué)者[46-48]低海拔地區(qū)土壤黏粒含量更高的研究結(jié)果不符。但也有學(xué)者與本文研究結(jié)果相似，BACIS等[49]在巴西杰尼羅州的研究表明高程與土壤質(zhì)地的變異密切相關(guān)，海拔越高，黏粒含量越高。這是因?yàn)榍罢遊46-48]的研究中，海拔主要影響土壤顆粒的運(yùn)移速度和方向，粒徑較小的顆粒在重力和地表徑流的作用下運(yùn)移的距離更遠(yuǎn)。BACIS等[49]的研究中，與海拔密切相關(guān)的是土壤母質(zhì)，該地區(qū)高海拔區(qū)域土壤主要由前寒武紀(jì)片麻巖原位風(fēng)化形成，質(zhì)地較細(xì)；低海拔區(qū)域以砂質(zhì)坡積物為主，砂粒含量較高。本研究中，海拔主要影響區(qū)域降水、植被等，研究區(qū)內(nèi)區(qū)域降水量差異大，以鹽池鄉(xiāng)-李旺一線為界，南部高海拔地區(qū)降水在400 mm以上，植被覆蓋度高，濕潤(rùn)的條件有利于黏土礦物的形成；北部低海拔地區(qū)降水為200—350 mm，植被覆蓋度低，地表裸露多，蒸發(fā)強(qiáng)烈，土壤顆粒較粗。這也與NDVI和砂粒含量呈顯著負(fù)相關(guān)，和粉粒和黏粒含量呈正相關(guān)的結(jié)果相一致。除高程外，坡度、風(fēng)力作用指數(shù)兩個(gè)地形因子也進(jìn)入RF模型（表5），但重要性相對(duì)較低，表明這兩個(gè)輔助變量雖然也會(huì)影響土壤顆粒組成的空間分布，但影響程度較小。坡度通過(guò)影響表層土壤顆粒運(yùn)動(dòng)、徑流挾沙能力和侵蝕方式[50-51]，從而影響土壤顆粒組成。鄒心雨等[52]等研究表明，坡度對(duì)小尺度下土壤顆粒組成的空間變異影響較大，本文研究尺度較大，坡度的影響并不顯著。

除地形因子外，土壤類(lèi)型也是預(yù)測(cè)土壤顆粒組成的重要變量。土壤類(lèi)型在一定程度上可以反映母質(zhì)、氣候等信息，楊艷麗等[53]研究表明，土壤類(lèi)型與黏粒含量的空間分布顯著相關(guān)。海原縣絕大部分區(qū)域?yàn)辄S土、黃土狀母質(zhì)，在此母質(zhì)上發(fā)育有黃綿土、灰鈣土和黑壚土，與黃綿土和灰鈣土相比，黑壚土形成的環(huán)境條件更為濕潤(rùn)，黏粒含量更高。海拔2 100 m以上的南華山、西華山中上部分布有頁(yè)巖、泥巖和片巖的風(fēng)化物，這些巖石風(fēng)化物是形成山地灰褐土的主要母質(zhì)，其中頁(yè)巖、泥巖屬于沉積巖中的泥質(zhì)巖類(lèi)，發(fā)育形成的土壤顆粒較細(xì)。新積土、粗骨土和紅黏土的面積較小，其中新積土是在水力或重力遷移堆積的物質(zhì)上形成，主要分布于東北角的沖積平原，顆粒較細(xì)；粗骨土僅分布于西華山等山地的陡坡，含較多的砂巖風(fēng)化碎屑，植被覆蓋度很低，土壤顆粒較粗。紅黏土主要分布于研究區(qū)東南角，為裸露的第三紀(jì)紅土層，質(zhì)地較黏重。

本研究區(qū)地形復(fù)雜，地形因子是重要的輔助變量；區(qū)內(nèi)土地利用類(lèi)型多為天然草地和旱地，且單季種植，土壤物理性狀受人為影響不強(qiáng)烈。因此要對(duì)平原地區(qū)、受人類(lèi)活動(dòng)影響強(qiáng)烈的地區(qū)進(jìn)行顆粒組成的空間預(yù)測(cè)，仍需探索選取適當(dāng)?shù)妮o助變量和預(yù)測(cè)方法。

4 結(jié)論

3種方法（SLR-EBK、SLR-RK和SLR-RF）預(yù)測(cè)的海原縣土壤各粒級(jí)含量空間分布的趨勢(shì)基本一致，表現(xiàn)為砂粒含量西南部低，東北部高，粉粒、黏粒則相反； SLR-RF、SLR-RK預(yù)測(cè)精度均高于SLR-EBK，其中以SLR-RF預(yù)測(cè)精度最高，為預(yù)測(cè)海原縣土壤顆粒組成空間分布的最優(yōu)方法；在黃土母質(zhì)區(qū)，高程、土壤類(lèi)型是與土壤顆粒組成的空間變異相關(guān)性較強(qiáng)的輔助變量。