黃土區不同土層土壤容重空間變異與模擬*

于冬雪 賈小旭? 黃來明 邵明安 王 嬌

（1 中國科學院地理科學與資源研究所生態系統網絡觀測與模擬重點實驗室，北京 100101）

（2 中國科學院大學資源與環境學院，北京 100190）

土壤容重（Bulk density, ρb）是土壤基本物理性質之一，是衡量土壤質量和土壤生產力的重要指標［1］，對土壤的透氣性、入滲性能、持水能力、溶質遷移特征以及土壤的抗侵蝕能力均有顯著影響［2］，也是評估土壤有機碳和養分貯量的重要參數［3-4］。土壤容重的大小與空間分布受土壤質地與結構、土地利用方式、地形、氣候等因素的影響［2,5-8］。黃土高原地形復雜多變、水土流失嚴重，了解黃土區土壤容重的空間變異規律及其影響因素，建立適于該區土壤容重的統計預測模型，有助于為該區的水土流失與侵蝕預報提供基本參數。

目前，已有大量學者對不同尺度下土壤容重的空間異質性進行研究，并取得重要進展［2-3,7,9-14］。耿韌等［11］研究了黃土丘陵區淺溝表層土壤容重的空間變異特征，發現土壤容重從淺溝上部到下部總體上呈逐漸減小的趨勢，土壤容重的空間異質性以各向同性為主；傅子洹等［12］結合經典統計學與地統計學方法，研究了黃土高原小流域尺度土壤容重的時空變化規律，發現小流域土壤容重在月際尺度上變化趨勢較為一致，表層土壤容重總體差異較小；易小波等［13］采用經典統計學方法分析了黃土高原南北樣帶不同土層深度土壤容重的變異特征，發現0～20 cm土壤容重為中等程度變異，20～40 cm為弱變異；柴華和何念鵬［14］基于全國3 361個樣地的11 845個土壤容重數據，全面闡釋了中國陸地生態系統土壤容重的基本特征和變異規律，結果表明，土壤容重隨土層深度增加而增大，不同土壤類型間的容重具有較大差異，這為準確評估區域乃至全國土壤碳、氮貯量提供了重要參數。

傳遞函數模型根據容易獲得的土壤理化性質可以快速獲得不同尺度上的水力參數，其構建方法主要有線性回歸、非線性回歸、神經網絡等［15］。由于大面積直接測量獲取土壤容重存在困難，國內外學者利用更易實現的土壤容重傳遞函數模型來預測土壤容重。Wang等［3］利用傳遞函數方程對黃土高原區域尺度表層土壤容重進行模擬，結果表明，基于黏粒和粉粒含量、有機碳含量、坡度建立土壤容重傳遞函數預測模型可以達到合理的精度。然而，該研究只分析了0～25 cm表層土壤容重，未進行土壤容重的分層模擬，而傳遞函數的重要應用之一就是對深層土壤容重進行模擬預測；韓光中等［8］對中國主要土壤類型的土壤容重傳遞函數研究發現，基于土壤系統分類的數據分組后建立的土壤容重傳遞函數能夠明顯提高預測精度，利用土壤容重傳遞函數時需注意研究區及適用范圍。

縱觀已有研究，對黃土高原區域尺度不同土層、不同土地利用方式下的土壤容重空間變異與模擬的研究較少。本研究獲取黃土區0～10、10～20和20～40 cm的土壤容重，采用經典統計學與地統計學方法，分析不同土層深度土壤容重的空間變異規律，并基于土壤、地形、土地利用和氣候因子，利用多元逐步回歸方程和傳遞函數方程對不同土層的土壤容重進行模擬預測，以期為該區的水土資源管理和生態建設提供基本參數和依據。

1 材料與方法

1.1 研究區概況

黃土高原位于中國西北地區（3 3°4 3′～41°16′N，100°54′～114°33′E），是指日月山、賀蘭山以東，太行山以西，秦嶺以北，陰山以南的廣大地區，面積約為64萬km2。該地區氣候呈現地帶性分布，降水量、平均氣溫、干燥度均由東南向西北方向遞減［16］。年降水量為150～800 mm，且集中在6—9月份，占全年降水量的55%～78%，多年平均降水量為466 mm，年平均氣溫為3.6～14.3 ℃，年蒸發量為1 400～2 000 mm。典型黃土區特指我國晉、陜、甘、寧、豫、內蒙諸省間黃土分布最連續、厚度最大、侵蝕地形最典型的地段（圖1），面積約為43萬km2，該區域塬、梁、峁、溝壑等黃土地貌發育最典型，水土流失最嚴重［17］，是黃土高原進行生態建設的重點區域。

圖1 黃土區位置以及采樣點分布Fig. 1 Location of the loess area and distribution of the sampling sites

1.2 數據采集與分析

2012年5—8月，以典型黃土區數字化地圖為底圖，采用間距約為40 km×40 km的網格進行取樣。具體為：在野外采樣之前，使用GIS軟件在典型黃土區的數字地形圖上加載一個40 km×40 km的格網作為參考；此外選擇合適的道路系統覆蓋整個區域構成采樣路徑，沿著選擇的道路方向取樣，采樣點之間間隔約為40 km。實際采樣過程中，具體采樣點的原則既要具有隨機性且能代表網格范圍內的主要土地利用類型、地形特征等。同時，要考慮采樣點的可達性和實際的可操作性。在每個樣點，利用環刀分3層（0～10、10～20和20～40 cm）采集原狀和擾動土壤樣品。將原狀土壤樣品帶回實驗室，采取烘干法測定其容重。擾動土壤樣品自然風干，處理通過1 mm篩孔，經預處理后利用激光粒度儀測定土壤顆粒組成。

為了更好地描述和體現各個網格內地表景觀的代表性，取樣時選擇最能代表該網格土壤和植被類型的樣點，若網格內地表景觀復雜多變，則考慮適當增加代表其他土地利用方式的樣點。經統計，本研究在黃土區共計布設243個樣點，其中農地（主要指1年生農作物地）46個，草地（包括人工草地和天然草地）88個，林地（包括天然和人工喬木林及灌木林地）109個（圖1）。

利用GPS記錄每個采樣點的經緯度坐標和海拔，使用羅盤儀測定樣點的坡度和坡向，同時記錄各樣點的土地利用方式和植被類型。為了分析方便，將土地利用方式分別進行數字編碼，農地=1，草地=2，林地=3［18］。利用黃土高原74個氣象站1951—2012年月尺度氣象數據，包括降水量、氣溫、蒸發量、干燥度（干燥度=蒸發量/降水量），使用反距離加權插值法將各氣象要素插值生成空間上連續分布的氣象數據（空間分辨率為100 m×100 m），然后提取每一樣點的氣象數據。由于土壤實際的蒸發量數據難以獲取，因此，本文利用氣象站常規觀測指標—水面蒸發量來反映樣點所在地區的蒸發能力。

1.3 經典統計學空間變異分析理論

土壤容重的空間變異特征用變異系數（CV）表示：

式中，S代表標準偏差，x 代表平均值。根據變異程度分級，CV≤100%為弱變異性，10%＜CV＜100%為中等變異性，C V≥1 0 0%為強變異性［2, 19］。

1.4 地統計學空間變異分析理論

地統計學進行空間變量變異特征分析是根據半方差變化規律進行的，半方差函數的計算公式如下：

式中，γ(h)代表半方差函數值，N(h)代表相距為的點對數，Z(xi)代表x=xi處變量Z的實測值，Z(xi+h)代表x=xi+h處變量的Z實測值。以滯后距h為橫軸，半方差函數值為縱軸可以繪制半方差圖。根據半方差圖可以得到塊金值C0，基臺值C0+C，變程A三個重要特征值。C0表示隨機部分引起的空間變異性，C0+C表示變量的最大變異程度，A表示變程即變量自相關范圍。通過公式C0/(C0+C)計算可以得到一個描述空間變量依賴性的重要指標即空間異質比，根據Cambardella的劃分標準［9］，C0/(C0+C)≤25%時表示強的空間依賴性，25%＜C0/(C0+C)＜75%時表示中等空間依賴性，C0/(C0+C)≥75%時表示弱的空間依賴性［2,10,20］。基于半方差函數的最佳擬合模型，利用普通克里格插值得到不同土層土壤容重的空間分布圖［10］。

1.5 逐步回歸方程和傳遞函數模型評價指標

本研究采用決定系數（R2）、平均誤差（AE）、平均絕對誤差（MAE）、均方根誤差（RMSE）和模型效率系數（MEC）五個參數來評價逐步回歸方程和傳遞函數模型的精度。其中，決定系數可以表示模型對土壤容重空間變異的解釋量，其值越大，說明模擬值所能解釋的變異越多。模型效率系數最早被用來評價水文模型的預測效果［21］，其值越接近于1表明精度越高。平均誤差、平均絕對誤差和均方根誤差越小，代表傳遞函數方程擬合的效果越好，預測的精度越高。

1.6 數據處理

使用Microsoft Excel2013和SPSS20.0進行土壤容重的描述性統計分析、土壤容重與環境因子的皮爾森相關分析和多元逐步回歸分析以及不同土層土壤容重的傳遞函數方程的模擬。利用GS+9.0進行半方差函數模型的擬合，采用ArcGIS10.2地統計模塊進行土壤容重的克里格插值及制圖。

2 結果與討論

2.1 土壤容重的基本統計特征

表1顯示了不同土層土壤容重的統計特征值。結果表明，3個土層深度中10～20 cm土壤容重的變化范圍最大。黃土區0～10、10～20和20～40 cm土壤容重的平均值分別為1.25、1.32、1.37 g·cm-3，這與易小波等［13］所測黃土高原南北樣帶3個土層土壤容重平均值基本一致，表明黃土區土壤容重隨土層深度的增加而增大，主要原因可能是表層土壤有機質含量高，植物根系也主要分布在淺層土壤，土壤孔隙狀況良好［22］。Wang等［3］和呂殿青等［23］研究表明容重與有機碳含量和孔隙數量呈負相關關系。經典統計結果表明，3個土層土壤容重的變異系數分別為10.0%、11.4%和10.4%，即3個深度的土壤容重均屬于中等程度變異。通過偏度、峰度以及單樣本Kolmogorov-Smirnov（K-S）檢驗來判斷容重數據是否符合正態分布，3個土層土壤容重的K-S值分別為0.697、0.130、0.221，表明不同土層容重均符合正態分布，可進行進一步的統計分析。

表1 不同土層土壤容重的描述性統計特征Table 1 Descriptive statistics characteristics of soil ρb relative to soil layer

2.2 基于地統計學的土壤容重空間變異特征

利用地統計半方差理論對不同土層土壤容重的空間變異性進行分析，結果見表2。0～10和10～20 cm土壤容重的最佳擬合模型均為指數模型，20～40 cm土壤容重最佳擬合模型為球狀模型，最佳擬合模型的決定系數均較高，表明所選模型能夠很好地反映土壤容重的空間變異特征。3個土層的塊金值變化范圍為0.002～0.013，表明隨機因素引起的土壤容重變異性較小。空間異質比分別為13.3%、50.0%和50.0%，表明0～10 cm土壤容重具有強烈的空間依賴性，10～20和20～40 cm土壤容重具有中等程度的空間依賴性［3,13］。因此，底層土壤容重由隨機因素引起的變異性更強。0～10 cm容重變程為22.3 km，10～20和20～40 cm變程分別為283.4 km和781.5 km，表明10～40 cm深度土壤容重空間連續性的尺度范圍更大。根據Flatman 和Yfantis［24］的建議，最佳采樣間距為變程的1/4～1/2，因此，3個土層深度容重最佳采樣間距分別為5.6～11.2 km、70.9～141.7 km和195.4～390.8 km。0～10 cm容重變程大于理論最佳采樣間距，為了達到反映其空間結構的目的，該土層深度需要減小采樣間距，增加采樣點數量，而10～40 cm土層可增大采樣間距，適當減少采樣數量。

表2 不同土層土壤容重的半變異函數模型及其結構參數Table 2 Semivariance model and structural parameters of soil ρb relative to soil layer

2.3 土壤容重的空間分布格局

為了直觀地反映黃土區不同土層土壤容重的空間分布情況，基于半方差函數的建立，采用普通克里格插值法對3個土層深度土壤容重進行空間插值。由圖2可以看出，不同土層土壤容重的連續性較好。總體而言，0～10 cm土層容重表現為黃土區內的甘肅、陜西中部和山西中部較低，20～40 cm土壤容重表現為黃土區內的陜西南部、河南和內蒙古較高。這可能是因為陜西南部和河南主要的土地利用類型為農地，農業活動導致犁底層被壓實，因而土壤容重較大。從局部來看，部分地區會出現土壤容重偏大或者偏小的值，這可能與人類的擾動、牧群的踐踏、鼠群的活動、植物的斑塊化生長有關［25］。

圖2 黃土區不同土層土壤容重的空間分布圖Fig. 2 Spatial distribution of soil ρb relative to soil layer in the loess area

2.4 土壤容重與各變量的Pearson相關性

土壤容重在不同尺度下受土壤類型、氣候、地形、土地利用以及生物擾動等因素的綜合影響［7,12］，因此，本研究選取11個可能影響土壤容重的環境因子（包括黏粒、粉粒、砂粒、土地利用方式、海拔、坡度、坡向、多年平均降水量、多年平均氣溫、蒸發量和干燥度）進行Pearson相關性分析，以期揭示區域尺度與土壤容重顯著相關的因素。對11個因子進行單樣本K-S檢驗，發現砂粒含量、坡度、多年平均氣溫和干燥度4個因子經過對數轉換后符合正態分布，故對此4個因子進行對數轉換后做進一步分析。

表3為不同土層土壤容重與土壤、氣候、地形和土地利用方式的皮爾森相關分析結果。可以看出，3個土層土壤容重與粉粒含量均呈極顯著負相關關系。粉粒含量的增加，可使土壤孔隙結構和透氣性有較大提高，從而導致土壤容重降低［3］。3個土層深度土壤容重均與海拔和坡度呈顯著負相關關系，可能是由于海拔較高、坡度較大的樣點土壤易被侵蝕，土壤較松散，以致土壤容重較小。此外，坡度較大的地方，受動物或人為活動干擾較小，故在一定程度上也可減小對土壤的壓實。此外，3個土層深度土壤容重均與土地利用方式呈顯著負相關關系。

針對土地利用方式對土壤容重的影響進行了單因素方差分析（圖3）。結果表明，對于同一種土地利用方式，土層越深，土壤容重越大；對于同一個土層，土壤容重表現為農地最大，草地與林地相近，這與連綱等［7］得到的灌木林地和荒草地土壤容重較大，林地和坡耕地土壤容重較小，梯田土壤容重居中的結果并不一致。這可能與研究區不同有關。本文以整個黃土區為研究對象，農地主要分布在黃土區南部，土壤黏粒含量高，土質密實，而林草地主要分布在黃土區中部和西北部，土壤顆粒組成以粉粒為主，土質疏松，加之林草地土壤根系含量高，人為活動導致的土壤壓實干擾小，所以林草地土壤容重較農地小。

2.5 土壤容重多元逐步回歸方程和傳遞函數方程

土壤容重的傳統測量方法與近年來出現的間接方法均很難實現大面積的連續測定，故造成土壤容重數據庫缺失。因此，利用土壤其他屬性來預測土壤容重對完善土壤容重數據庫具有重要意義［8,26］。本研究根據皮爾森相關分析結果，利用243個樣點中的80%（195個）作為建模點。選取與每個土層土壤容重顯著相關的因子進行多元逐步回歸分析，得到3個土層的多元逐步回歸方程（表4）；結合與每個土層土壤容重顯著相關的因子及其6種常見的轉換函數（對數、倒數、平方、開方、余弦、顯著相關因子兩兩相乘）作為自變量進行逐步回歸分析，將通過顯著回歸的因子進行線性回歸，得到3個土層土壤容重的傳遞函數方程（表4）。

結果表明，3個多元逐步回歸方程均包含粉粒含量、海拔、土地利用方式。不同土層深度的傳遞函數方程包含的變量不同，0～10 cm主要為粉粒含量、土地利用方式、海拔以及坡度，10～20 cm主要為粉粒含量、海拔、多年平均氣溫、干燥度及土地利用方式，20～40 cm主要為粉粒含量、海拔、土地利用方式、多年平均降水量、坡度及干燥度。結合皮爾森相關分析、多元逐步回歸方程以及傳遞函數方程模擬結果可以看出，黃土區土壤容重受土壤、地形、土地利用方式、氣候等大尺度與小尺度因子的綜合影響。

表3 不同土層土壤容重與各變量的Pearson相關系數Table 3 Pearson's correlation coefficients between soil ρb and other variables relative to soil layer

圖3 不同土地利用方式下土壤容重單因素方差分析（平均值±標準誤差）Fig. 3 One-way analysis of variance of soil ρb relative to land use (mean ± standard error)

利用243個樣點中的20%（48個）進行模型精度驗證，驗證結果見表4。多元逐步回歸方程3個土層的決定系數R2分別為0.39、0.35、0.34，傳遞函數方程3個土層土壤容重的決定系數R2分別為0.40、0.38和0.52。因此，兩種方法對0～10 cm土層土壤容重的模擬效果相近，可以解釋土壤容重空間變異的40%。與0～10 cm模擬結果不同，10～20和20～40 cm深度傳遞函數模型模擬效果明顯優于逐步回歸方程，模擬精度較逐步回歸方程分別提高3%和18%。傳遞函數方程平均誤差、平均絕對誤差均小于0.1，均方根誤差均小于0.2，模型效率系數也均超過0.3，精度驗證參數總體優于多元逐步回歸方程，表明傳遞函數方程在預測底層土壤容重方面效果更佳。這主要是因為傳遞函數方程包含了不同變量的轉換，可以更好地描述變量與容重之間的關系［3］。同時，根據Leij等［27］提出的地形因子可以提高傳遞函數模型的效率，本研究將坡度、海拔等地形因子均考慮在內。此外，Han等［28］研究表明土壤容重與土層深度呈顯著正相關，本研究對不同土層土壤容重分別進行模擬，可剔除土層深度對模擬效果的干擾。與Wang等［3］研究相比，本研究采集的樣點土層更深，因而對深層土壤容重研究更具有指導意義。此外，傳遞函數模型中容重的影響因子也較容易獲取，因此，該傳遞函數模型可用于田間條件下黃土高原區域尺度不同土層土壤容重的模擬與預測。

表4 不同土層土壤容重模擬的多元逐步回歸方程（a）、傳遞函數方程（b）和預測精度Table 4 Simulation of soil ρb with multi stepwise regression equation(a), and pedotransfer function equation(b) and their prediction accuracies

3 結 論

黃土區不同土層、不同土地利用方式下土壤容重的空間分布具有顯著差異，土壤容重隨土層深度的增加而增大，且不同土層中農地的容重最大，其次為林地和草地。0～10 cm容重具有強烈的空間依賴性，而10～40 cm具有中等程度空間依賴性，各土層容重的半變異函數可用指數模型和球狀模型進行較好擬合。0～10、10～20和20～40 cm深度土壤容重最佳采樣間距分別為5.6～11.2、70.9～141.7和195.4～390.8 km。黃土區不同土層土壤容重空間變異受土壤、地形、氣候和土地利用方式的共同影響。0～10 cm主要為粉粒含量、土地利用、海拔以及坡度，10～20 cm主要為粉粒含量、海拔、多年平均氣溫、干燥度及土地利用，20～40 cm主要為粉粒含量、海拔、土地利用、多年平均降水量、坡度及干燥度。傳遞函數模型對底層土壤容重模擬效果優于多元逐步回歸方程，可用于田間復雜環境條件下區域尺度土壤容重空間分布特征的預測。