基于3種空間預測方法的黃土區土壤顆粒組成空間分布研究—以寧夏海原縣為例

申哲，張認連，龍懷玉，王轉，朱國龍，石乾雄，喻科凡，徐愛國

（1中國農業科學院農業資源與農業區劃研究所，北京 100081；2北京大學地球與空間科學學院，北京100871）

0 引言

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

1 材料與方法

1.1 研究區概況

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

表1 各土壤類型面積占比及樣點統計Table 1 The area proportion and sample points of each soil type

圖1 研究區位置及樣點分布圖Fig. 1 The map of the study position and sampling sites

1.2 樣品采集與分析

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

1.3 輔助變量來源及處理

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

1.3.1 輔助變量來源

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

1.3.2 輔助變量的提取及處理

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

1.4 預測方法

1.4.1 數據轉換方法 為了使數據滿足正態分布，且使預測結果滿足非負、定和為1的要求，本文首先對土壤顆粒組成進行對稱對數比（symmetry log-ratio，SLR）轉換[32]，轉換公式為：

轉回公式為：

式中，Zij代表第i個樣點第j種粒級的相對含量（%）；代表第i個樣點第j種粒級含量的轉換值；D表示成分數據的維度，本文中D=3；δj為常數，取研究區第j種顆粒除0外最小含量的一半。

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

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

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

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

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

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

1.5 精度驗證

從124個樣點中隨機選取80%的（100個）樣點用于空間預測，20%的（24個）樣點用于精度驗證[39-40]，訓練樣點和驗證樣點的分布如圖1所示。本研究采用驗證集的平均絕對誤差（MAE）、均方根預測誤差（RMSE）、平均Aitchison距離（MAD）[41]進行模型精度評價。MAE能更好地反映預測值誤差的實際情況；RMSE對系統誤差和隨機誤差都很敏感，常用來評估不同模型的預測精度；MAD可以反映成分數據預測值和觀測值之間整體的相似性。MAE、 RMSE、MAD越小，模型精度越高。

Zij代表第i個樣點第j種粒級含量的觀測值；代表第i個樣點第j種粒級含量的預測值；n表示驗證點的數量；D表示成分數據的維度，本文中D=3。

2 結果

2.1 土壤顆粒組成的描述性統計特征

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

2.2 多元線性回歸和RF模型的構建

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

表3 研究區采樣點土壤顆粒組成基本統計特征Table 3 Statistical characters of soil particle size distribution in study area

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

表5 SLR-RF參數擬合結果Table 5 SLR-RF parameter fitting results

2.3 土壤顆粒組成的空間分布預測結果

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

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

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

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

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

2.4 不同預測方法的精度對比

表6列出了3種方法在驗證集的預測精度，可以看出SLR-RF法對3個粒級含量預測的MAE和RMSE均低于其他兩種方法，且SLR-RF的MAD（0.208）最低，表明該方法的預測精度更高。因此，研究區土壤顆粒組成的最優預測模型為SLR-RF。

表6 不同空間預測方法精度對比Table 6 Accuracy comparison of different spatial prediction methods

3 討論

3.1 3種預測方法的精度

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

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

3.2 土壤顆粒組成空間分布與輔助變量的關系

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

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

本研究區地形復雜，地形因子是重要的輔助變量；區內土地利用類型多為天然草地和旱地，且單季種植，土壤物理性狀受人為影響不強烈。因此要對平原地區、受人類活動影響強烈的地區進行顆粒組成的空間預測，仍需探索選取適當的輔助變量和預測方法。

4 結論

3種方法（SLR-EBK、SLR-RK和SLR-RF）預測的海原縣土壤各粒級含量空間分布的趨勢基本一致，表現為砂粒含量西南部低，東北部高，粉粒、黏粒則相反； SLR-RF、SLR-RK預測精度均高于SLR-EBK，其中以SLR-RF預測精度最高，為預測海原縣土壤顆粒組成空間分布的最優方法；在黃土母質區，高程、土壤類型是與土壤顆粒組成的空間變異相關性較強的輔助變量。