黃土高原南北樣帶不同土層土壤容重變異分析與模擬

易小波 邵明安,2 趙春雷 張晨成,2

(1.西北農(nóng)林科技大學資源環(huán)境學院， 陜西楊凌 712100；2.西北農(nóng)林科技大學黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點實驗室， 陜西楊凌 712100)

黃土高原南北樣帶不同土層土壤容重變異分析與模擬

易小波1邵明安1,2趙春雷1張晨成1,2

(1.西北農(nóng)林科技大學資源環(huán)境學院， 陜西楊凌 712100；2.西北農(nóng)林科技大學黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點實驗室， 陜西楊凌 712100)

為探明黃土高原南北樣帶土壤容重空間分布特征，為土壤水文過程模擬與預測提供水力參數(shù)，采用經(jīng)典統(tǒng)計學方法，分析了樣帶不同土層深度(0～10 cm、10～20 cm、20～40 cm)土壤容重的空間變異特征，并用多元逐步回歸、傳遞函數(shù)方程和一階自回歸狀態(tài)空間模型方法分別對土壤容重的空間分布進行了模擬。結(jié)果表明：樣帶0～20 cm深度土壤容重的變異為中等程度變異，20～40 cm為弱變異。狀態(tài)空間方程轉(zhuǎn)換系數(shù)表明，不同土層深度土壤容重的影響因素不同，0～10 cm主要為有機碳含量、黏粒和砂粒體積分數(shù)，10～20 cm為有機碳含量、黏粒和砂粒體積分數(shù)和降水量，20～40 cm為黏粒和砂粒體積分數(shù)、降水量和土地利用。狀態(tài)空間模型的模擬效果均優(yōu)于經(jīng)典統(tǒng)計的多元逐步回歸方程和傳遞函數(shù)方程，基于黏粒和砂粒體積分數(shù)、降水量和土地利用的狀態(tài)空間模型可以解釋樣帶20～40 cm容重92.3%的變異。一階自回歸狀態(tài)空間模型可用于田間條件下土壤容重分布特征的預測。

容重； 狀態(tài)空間； 傳遞函數(shù)； 模擬； 空間變異

引言

容重是土壤重要的基本物理性質(zhì)，對土壤的透氣性、入滲性、持水性、溶質(zhì)遷移以及土壤抗蝕能力等特性都有重要影響，是衡量土壤環(huán)境好壞的重要指標之一[1-2]。土壤容重通過改變土壤硬度影響植物的生長，隨著容重的提高，土壤硬度隨之增大，土壤中根系伸長速度減慢且根變短變粗[3]。容重對于養(yǎng)分在土壤中的擴散也具有十分顯著的影響，由于土壤固體顆粒的排列不同，造成了不同的曲折通路，單位土體內(nèi)的電荷密度隨容重的改變而改變，進而影響擴散系數(shù)[4]。此外，土壤容重也是準確估算土壤碳氮儲量所必需的重要參數(shù)[1]。

土壤容重的傳統(tǒng)測量方法——環(huán)刀法，操作耗時，且在野外大尺度實際應用中存在工作繁重，人力、物力耗費較大等問題[5]。近年來，土壤容重測量的新方法也不斷涌現(xiàn)。研究人員利用土壤質(zhì)量含水率、體積含水率和容重的關(guān)系，分別采用MP水分分析儀和TDR反射水分儀，同時結(jié)合土鉆法間接測量土壤容重。但是，受限于土壤體積含水率獲取的局限性，這種間接測量方法仍不能實現(xiàn)野外大面積的實時連續(xù)測量[6]。土壤容重預測模型作為簡便、快捷的土壤容重獲取方法，越來越多地受到國內(nèi)外學者的關(guān)注和應用[7-10]。已有的土壤容重預測模型中，所需要的輸入?yún)?shù)過多，而且大部分模型模擬精度相對較低，嚴重影響了模型的實用性[8-9]。WANG等[7]利用傳遞函數(shù)方程對黃土高原區(qū)域尺度表層土壤容重進行了模擬，并基于黏粒和粉粒體積分數(shù)、有機碳含量、坡度及其轉(zhuǎn)換組合建立了傳遞函數(shù)預測模型。然而，其研究并未考慮土層深度對土壤容重模擬的影響。

由于成土母質(zhì)、過程、氣候、生物及耕作的影響，自然條件下的土壤容重具有較高的空間變異性[11]。傳統(tǒng)的空間變異理論未考慮土壤參數(shù)與地理位置間的空間關(guān)系，對空間變異的解釋不夠充分，狀態(tài)空間方程對土壤及植被的模擬效果要明顯優(yōu)于線性回歸和多元回歸等傳統(tǒng)方法[12]。賈小旭等[13]對黃土高原北部典型坡面草地表層土壤水分進行了狀態(tài)空間模擬，結(jié)果表明狀態(tài)空間方程可用于環(huán)境因素復雜的黃土高原水蝕風蝕交錯區(qū)，其模擬效果明顯優(yōu)于經(jīng)典統(tǒng)計的線性回歸模型。趙春雷等[14]對黃土高原坡面尺度土壤飽和導水率進行模擬，結(jié)果表明基于容重和土壤顆粒的狀態(tài)空間方程可以解釋坡面飽和導水率90%以上的變異。因此，與經(jīng)典統(tǒng)計的回歸方程相比，一階自回歸狀態(tài)空間模型更能揭示空間上變量與變量之間的關(guān)系，是一種對變量進行模擬和預測的有效工具。

本文選取可能影響土壤容重的環(huán)境因素，包括黏粒和砂粒體積分數(shù)、有機碳含量、降水量、土地利用、海拔、坡度和坡向等作為潛在的預測變量，用一階自回歸狀態(tài)空間模型及經(jīng)典統(tǒng)計的多元逐步回歸和傳遞函數(shù)方程分別對不同深度土層土壤容重的空間分布狀況進行模擬，探索區(qū)域樣帶尺度不同土層深度土壤容重的空間分布特征及主要影響因素，以期為該區(qū)土壤水文過程模擬與預測提供水力參數(shù)。

1 材料與方法

1.1 研究區(qū)概況

黃土高原指日月山、賀蘭山以東，太行山以西，秦嶺以北，陰山以南的廣大地區(qū)，面積64萬km2。黃土高原年均降水量為150～800 mm，由東南向西北方向遞減，多年平均降水量466 mm，多集中于6—9月份，占全年降水量的55%～78%。黃土高原屬典型的溫帶大陸性季風氣候，冬春季寒冷干燥，夏秋季溫暖濕潤，雨熱同步，年平均氣溫3.6～14.3℃。該區(qū)光能資源豐富，年總輻射量為(5.02～6.70)×105J/cm2，呈東南向西北遞增趨勢。黃土高土壤主要為風成黃土，土層深厚，厚度多在50～200 m，部分地區(qū)深達300 m，土質(zhì)疏松，易于侵蝕，且具有明顯的地帶性特征，由東南向西北，依次出現(xiàn)褐色土、塿土、黑壚土、黃綿土、灰鈣土、栗鈣土、棕鈣土、棕漠土、風沙土等。黃土高原地貌類型主要有丘陵、塬、茆、溝、梁、階地、平原、土石山地等[15]。

1.2 樣品采集

2013年4月份在黃土高原南北方向布設一維樣帶，總長860 km。樣帶具有明顯的環(huán)境梯度條件，年均降水量從南到北變化范圍為620～400 mm，年均氣溫12.3～6.8℃，干燥度2.3～5.6，海拔高度383～1 529 m(圖1)；土壤類型從南到北依次為塿土、黃綿土、風沙土；土地利用依次為農(nóng)地(主要指一年生農(nóng)作物地)、林地(包括天然/人工喬木林及灌木林地)、草地(包括人工草地和天然草地)[16]。根據(jù)地表景觀和地形特征，沿樣帶每隔10 km布設一個樣點，共86個樣點，由南至北依次編號為1～86。在樣帶86個樣點中，包含10個農(nóng)地、54個林地和22個草地，為定量刻畫土地利用對土壤容重的影響，將農(nóng)地、林地和草地分別表示為0、1和2[7](圖1)。王云強等[17]基于地統(tǒng)計分析發(fā)現(xiàn)，黃土高原區(qū)域尺度表層(0～5 cm)和底層(20～25 cm)土壤容重變程分別為448 km和300 km。根據(jù)變程可制定合理的采樣策略，F(xiàn)LATMAN等[18]建議最佳的采樣間距為變程的1/4～1/2之間。因此，研究黃土高原區(qū)域尺度土壤容重的空間變異性，表層土壤容重最低的采樣間距應為112～224 km，底層應為75～150 km。本研究以典型黃土區(qū)為對象，結(jié)合黃土高原地區(qū)的數(shù)字地形圖(DEM)和自然地理條件，采用高密度布點(樣點間距約為10 km)策略，采樣密度在其變程范圍內(nèi)，因此，樣點設置合理，所得數(shù)據(jù)足以滿足空間上的自相關(guān)要求。

圖1 降水量、氣溫、相對高程和土地利用沿樣帶分布情況Fig.1 Distributions of precipitation, air temperature, relative elevation and land use types along south to north transect

利用RTK-GPS記錄每個樣點的經(jīng)緯度坐標及海拔高度，利用羅盤儀測定每個樣點的坡度和坡向，同時記錄樣點土地利用方式。在每個樣點挖取40 cm深土壤剖面，用100 cm3的環(huán)刀分3層采集原狀土壤樣品(0～10 cm、10～20 cm、20～40 cm)。將原狀土樣帶回實驗室，利用烘干法測定土壤容重。同時采集不同土層深度的擾動土壤樣品，風干后分別過0.25 mm和1.0 mm篩預處理，用于土壤黏粒、砂粒體積分數(shù)和有機碳含量測定[19]。

利用黃土高原68個國家氣象觀測站1951—2012年月值氣象數(shù)據(jù)，包括月均降水量和氣溫，使用反距離加權(quán)插值方法將氣象要素插值生成空間上連續(xù)分布的氣象數(shù)據(jù)(空間分辨率為100 m×100 m)，然后提取每一樣點的降水量和氣溫數(shù)據(jù)[16]。

1.3 分析方法

1.3.1 經(jīng)典統(tǒng)計學空間變異分析理論

經(jīng)典統(tǒng)計學進行空間變異分析時要求分析的空間變量在研究區(qū)域內(nèi)符合隨機分布特征，對變異特征的分析主要以變異系數(shù)Cv進行，其定義為

(1)

根據(jù)變異程度分級，Cv≤10%屬于弱變異性，10%

1.3.2 狀態(tài)空間方程

狀態(tài)空間方程是基于馬爾可夫(Markovian)過程的特征建立起來的，主要用于描述在同一地點采集到的j組數(shù)據(jù)Yj(xi)的隨機過程。根據(jù)給定現(xiàn)有狀態(tài)的數(shù)據(jù)，可以依據(jù)過去的過程，預測未來的狀況。在這些系統(tǒng)中，過程狀態(tài)包括用來預測將來所需的過去所有信息[20-21]。狀態(tài)空間模型包含2個方程：狀態(tài)方程和觀測方程。系統(tǒng)在位置i處的狀態(tài)通過狀態(tài)空間方程與位置i-h(h=1，2，…，n)的狀態(tài)方程發(fā)生聯(lián)系。當h=1時，狀態(tài)方程可表示為

Xi=φXi-1+wi

(2)

式中Xi——在位置i處的狀態(tài)向量，通常指若干個變量在i處的值

φ——狀態(tài)轉(zhuǎn)化系數(shù)矩陣

wi——模型誤差

φ是一個表示狀態(tài)在位置i與i-1處的空間線性回歸情況的p′×p′矩陣，wi也是一個普通的p′×p′協(xié)方差矩陣，且在空間位置上不相關(guān)、方差恒定呈正態(tài)分布且均值為0。狀態(tài)向量以及轉(zhuǎn)化系數(shù)矩陣中的參數(shù)可通過迭代回歸求得。狀態(tài)向量Xi通過與觀測向量相聯(lián)系，從而產(chǎn)生觀測方程。觀測方程為

Yi=AiXi+vi

(3)

式中vi——在空間位置上不相關(guān)，方差恒定呈正態(tài)分布，且均值為0的協(xié)方差矩陣

在觀測方程中，通過觀測矩陣Ai以及觀測誤差vi將系統(tǒng)的觀測向量和系統(tǒng)的真實狀態(tài)向量聯(lián)系起來。狀態(tài)方程和觀測方程中2個誤差項wi和vi相互獨立。在應用過程中，Xi為系統(tǒng)的真實狀態(tài)，Yi為測量值，但不一定能全部反映系統(tǒng)的真實狀態(tài)，而只是系統(tǒng)狀態(tài)的間接反映，它等于系統(tǒng)的真實狀態(tài)再加上一個不能確定的誤差。狀態(tài)方程和觀測方程的求解同時進行，狀態(tài)向量以及轉(zhuǎn)化系數(shù)矩陣φ中的參數(shù)可通過迭代回歸求得[22]。

為消除不同量綱、變量自身變異和數(shù)值對模擬結(jié)果的影響，狀態(tài)空間模擬要求對數(shù)據(jù)進行標準化，標準化公式為

zi=[Zi-(m-2s)]/(4s)

(4)

式中Zi——標準化前的變量zi——標準化后的變量m——標準化前變量的平均值

標準化后，各變量的均值為0.5，標準差為0.25。標準化后可以根據(jù)轉(zhuǎn)換系數(shù)來反映變量之間的相對貢獻率[23]。本文各種模擬分析中均采用標準化后的值。在評價狀態(tài)空間模型的優(yōu)劣時，統(tǒng)一使用決定系數(shù)(R2)和均方根誤差(RMSE)來評價模型的優(yōu)劣。

決定系數(shù)R2可揭示模型對樣帶土壤容重空間變異的解釋量。R2值越大，表示模擬值所能解釋土壤容重的變異越多。R2值越大，RMSE值越小，表示模擬效果越好。

1.4 數(shù)據(jù)處理

本研究中不同土層深度下土壤容重的統(tǒng)計特征值、Pearson相關(guān)分析、土壤容重的傳遞函數(shù)方程和多元逐步回歸方程的建立均利用SPSS 13.0進行。土壤容重的狀態(tài)空間模型在ASTSA軟件下進行。

2 結(jié)果與分析

2.1 土壤容重、有機碳含量、黏粒和砂粒體積分數(shù)沿樣帶分布情況

圖1顯示了樣點降水量、氣溫、相對高程和土地利用沿樣帶的分布狀況。樣點降水量和氣溫由南向北逐漸降低，這由研究區(qū)氣候特征決定。1～10號樣點相對高程較低，11～86號樣點相對高程較高，這是因為1～10號樣點分布在關(guān)中平原地區(qū)，而11～86號樣點分布于黃土高原丘陵溝壑區(qū)，海拔高度較高，地形地貌復雜[24]。南北樣帶0～10 cm、10～20 cm和20～40 cm深度土壤容重的平均值分別為1.24、1.33、1.37 g/cm3(表1和圖2)，表明隨土層深度的增加，土壤容重逐漸增大。這可能與土壤有機碳含量有關(guān)，本研究中0～10 cm土壤有機碳含量高于10～20 cm土層，20～40 cm土壤有機碳含量最低(圖2)，這與劉志鵬等[12]在區(qū)域尺度下的研究結(jié)果一致，有機碳含量在一定程度上影響著容重沿剖面的分布。WANG等[7]分析研究了黃土高原地區(qū)0～25 cm土壤容重分布狀況，結(jié)果表明容重與有機碳含量顯著負相關(guān)。BENITES等[25]通過研究巴西區(qū)域尺度0～100 cm土壤容重狀況，也發(fā)現(xiàn)有機碳含量與容重存在顯著的負相關(guān)關(guān)系。樣帶0～10 cm和10～20 cm土層深度土壤容重的變異系數(shù)分別為11.01%和11.07%，為中等程度變異；20～40 cm土層容重變異系數(shù)僅為8.82%，為弱變異。這主要是由于上層土壤相較于下層土壤更容易受到環(huán)境因素的影響[7]，因此，空間變異性較強。樣帶砂粒和黏粒體積分數(shù)沿樣帶分別呈現(xiàn)出逐漸增加和減小的趨勢，這主要是由于黃土高原形成過程中風搬運顆粒逐級分選造成的[26]。

表1 不同土層深度下土壤容重統(tǒng)計特征值

注：N為服從正態(tài)分布。

2.2 容重與環(huán)境變量的Pearson相關(guān)分析

表2為所選變量與不同土層深度土壤容重(B)的相關(guān)系數(shù)。

結(jié)果表明，0～10 cm土壤容重與有機碳含量、黏粒體積分數(shù)、降水量和坡度顯著負相關(guān)，與土地利用、海拔高度和坡向無顯著相關(guān)關(guān)系；10～20 cm土壤容重與有機碳含量、黏粒體積分數(shù)、降水量和土地利用顯著負相關(guān)，與海拔高度、坡度和坡向無顯著相關(guān)關(guān)系；20～40 cm土壤容重與黏粒體積分數(shù)、降水量、土地利用和海拔高度顯著負相關(guān)，與有機碳含量、坡度和坡向無顯著相關(guān)關(guān)系；此外，各土層容重均與砂粒體積分數(shù)呈極顯著正相關(guān)關(guān)系，表明砂粒體積分數(shù)越高，土壤容重越大，這與WANG等[7]在黃土高原的研究結(jié)果一致。本研究只選擇與不同土層深度土壤容重顯著相關(guān)的變量分別建立狀態(tài)空間方程、傳遞函數(shù)方程和多元逐步回歸方程，篩選出不同方法最優(yōu)的預測方程，并比較不同方程對樣帶不同深度土壤容重的模擬效果。因此，對于0～10 cm土層，以有機碳含量、黏粒和砂粒體積分數(shù)、降水量和坡度作為預測變量；對于10～20 cm土層，以有機碳含量、黏粒和砂粒體積分數(shù)、降水量和土地利用作為預測變量；對于20～40 cm土層，以黏粒和砂粒體積分數(shù)、降水量、土地利用和海拔高度作為預測變量。

2.3 土壤容重的多元逐步回歸方程和傳遞函數(shù)方程

多元逐步回歸方程作為經(jīng)典的統(tǒng)計預測方法已得到廣泛應用。本文采用經(jīng)典統(tǒng)計的多元逐步回歸方法，將與土壤容重顯著相關(guān)的變量作為自變量對因變量土壤容重進行線性回歸模擬。不同土層深度多元逐步回歸方程均經(jīng)過顯著性檢驗，且均在p<0.01 水平上顯著。在3個多元逐步回歸方程中均包含了黏粒和砂粒體積分數(shù)，表明在研究區(qū)，黏粒和砂粒體積分數(shù)是影響土壤容重的重要因素(表3)。一般來說，隨著砂粒體積分數(shù)的增加和黏粒體積分數(shù)的減小，由于土壤總孔隙度的減小，容重呈現(xiàn)增加的趨勢。 WANG等[7]運用人工神經(jīng)網(wǎng)絡和線性回歸的方法模擬預測了黃土高原區(qū)域尺度表層土壤容重狀況，研究結(jié)果表明，包含黏粒和砂粒體積分數(shù)等因素在內(nèi)的線性回歸模型可以較好地預測區(qū)域尺度下0～25 cm深度土壤容重。本研究中，不同土層深度的逐步回歸方程決定系數(shù)R2為0.444～0.493，均方根誤差RMSE為0.177～0.185 g/cm3，表明在當前研究尺度下，模型變量所能解釋的土壤容重空間變異僅為44.4%～49.3%，即模擬精度較低。

圖2 不同土層深度土壤容重、有機碳含量、黏粒和砂粒體積分數(shù)沿樣帶分布情況Fig.2 Distributions of soil bulk density, soil organic carbon content and clay and sand contents in different soil layers along south to north transect

表2 不同土層深度土壤容重與各變量的相關(guān)系數(shù)

注：*表示相關(guān)關(guān)系在p<0.05水平上具有顯著性(雙尾檢驗)；** 表示相關(guān)關(guān)系在p<0.01水平上具有顯著性(雙尾檢驗)。

為了提高統(tǒng)計方程的模擬精度，進一步將與土壤容重顯著相關(guān)的變量以及其4種常見轉(zhuǎn)換(指數(shù)、對數(shù)、倒數(shù)和平方轉(zhuǎn)換)作為自變量再次進行逐步回歸分析，以期得到對土壤容重空間變異具有更強解釋能力的傳遞函數(shù)方程，并選出各土層模擬精度最優(yōu)的傳遞函數(shù)方程(表3)。不同土層深度傳遞函數(shù)方程均經(jīng)過顯著性檢驗，且均在p<0.01水平上顯著。結(jié)果表明，不同土層深度包含了變量及其常見轉(zhuǎn)換的傳遞函數(shù)方程在模擬精度上得到了不同程度的提高。傳遞函數(shù)方程分別可以解釋0～10 cm、10～20 cm和20～40 cm土層深度土壤容重空間變異的50.4%、58.7%和54.0%，RMSE分別為0.174、0.154、0.165 g/cm3，均低于多元逐步回歸方程的RMSE。與多元逐步回歸方程相比，0～10 cm、10～20 cm和20～40 cm土層深度傳遞函數(shù)方程模擬精度分別提高了6.3%、9.5%和32.2%。這主要是因為傳遞函數(shù)方程包含了不同變量的轉(zhuǎn)換，可以更好地刻畫變量與容重之間的關(guān)系[7]。

表3 不同土層深度下土壤容重模擬的狀態(tài)空間方程、傳遞函數(shù)方程和多元逐步回歸方程

2.4 土壤容重的狀態(tài)空間模擬

將與不同土層深度土壤容重顯著相關(guān)的變量及其不同因素組合分別引入狀態(tài)空間模型中，以期找出能夠模擬樣帶不同土層深度土壤容重的最佳因素組合。所有狀態(tài)空間方程均經(jīng)過顯著性檢驗，且均在p<0.01水平上顯著。表3為不同土層深度土壤容重模擬效果最好的狀態(tài)空間方程。結(jié)果表明，不同土層深度最優(yōu)狀態(tài)空間方程包含變量不同，0～10 cm為有機碳含量、黏粒和砂粒體積分數(shù)；10～20 cm為有機碳含量、黏粒和砂粒體積分數(shù)和降水量；20～40 cm為黏粒和砂粒體積分數(shù)、降水量和土地利用。由決定系數(shù)R2可知，最優(yōu)狀態(tài)空間方程分別可以解釋樣帶0～10 cm、10～20 cm和20～40 cm土層深度土壤容重空間變異信息的69.1%、72.7%和92.3%。因此，基于黏粒和砂粒體積分數(shù)、降水量和土地利用的狀態(tài)空間方程對樣帶20～40 cm容重模擬效果最好，且95%置信區(qū)間最窄。圖3顯示了不同土層深度最優(yōu)狀態(tài)空間方程以及土壤容重實測值和模擬值的分布情況。從圖3可以看出，20～40 cm土層土壤容重的模擬值與實測值重合性很好，且置信區(qū)間相對較窄，置信區(qū)間內(nèi)包括了所有的實測值，而0～10 cm和10～20 cm土層模型的置信區(qū)間相對較寬且有部分實測值在置信區(qū)間以外，其模擬效果遜于20～40 cm，這可能與不同土層深度土壤容重的變異強弱有關(guān)。0～20 cm深度土壤容重為中等程度變異，20～40 cm為弱變異(表 1)，即表層土壤容重比深層土壤容重可能受外界因素影響更加復雜，而依據(jù)有限的變量表層模擬精度要低于深層。

圖3 不同土層深度模擬土壤容重最優(yōu)的狀態(tài)空間模型Fig.3 Optimal state-space models for simulating soil bulk density in different soil layers

圖4 不同土層深度最優(yōu)狀態(tài)空間模型中各變量對土壤容重的相對貢獻率Fig.4 Relative contribution of each variable in optimal state-space model for soil bulk density in different depths

根據(jù)模擬結(jié)果的轉(zhuǎn)換系數(shù)，可計算各狀態(tài)空間方程中i-1處變量對相鄰位置處土壤容重的相對貢獻率(圖4)。在0～10 cm土層，有機碳含量、黏粒和砂粒體積分數(shù)對容重的相對貢獻率分別為16.1%、9.4%和14.8%；在10～20 cm土層，有機碳含量、黏粒和砂粒體積分數(shù)和降水量對容重的相對貢獻率分別為34.2%、5.1%、9.1%和23.7%；在20～40 cm土層，黏粒和砂粒體積分數(shù)、降水量和土地利用對容重的相對貢獻率分別為28.8%、40.4%、12.6%和17.3%。因此，在0～20 cm土層，有機碳含量對土壤容重空間分布的相對貢獻率最大，是影響土壤容重空間分布的主要因素。然而，10～20 cm土層中有機碳含量的相對貢獻率高于0～10 cm，這主要是由于表層土壤容重的影響因素相對復雜，有機碳含量對其變異的刻畫能力相對較弱[25]。在20～40 cm土層，有機碳含量沒有被引入最優(yōu)的狀態(tài)空間方程，這主要是由于深層土壤中有機碳含量較低且變異較小(圖2)，這與KAUR等[27]的結(jié)論一致。值得注意的是，在表層0～10 cm，影響容重空間分布主要是土壤因素，而在10～20 cm和20～40 cm土層，氣候因素和土地利用應當被考慮。

不同土層深度最優(yōu)的狀態(tài)空間方程的RMSE分別為0.118、0.114、0.028 g/cm3，遠低于多元逐步回歸方程和傳遞函數(shù)方程。因此，由R2和RMSE可知，在3個土層中，狀態(tài)空間方程對樣帶土壤容重空間分布的模擬效果均優(yōu)于經(jīng)典統(tǒng)計的多元逐步回歸方程和傳遞函數(shù)方程，可解釋土壤容重空間變異的69.1%～92.3%。這主要是因為經(jīng)典統(tǒng)計的回歸分析不考慮采樣點的空間位置，假設樣點與樣點之間是相互獨立的，即在整個研究區(qū)域上把土壤容重與其他因素的相關(guān)性看作是單一的，而忽略容重與其他因素的關(guān)系在空間位置上的變異性[28]，而狀態(tài)空間模型則考慮了變量之間相關(guān)性的空間變異。此外，狀態(tài)空間方程考慮了模型的誤差和觀測誤差，進而通過迭代算法調(diào)整模型參數(shù)，最終獲取預測精度較高的模型參數(shù)[12]。因此，相對于經(jīng)典統(tǒng)計的回歸分析，狀態(tài)空間模型更能揭示其他因素與土壤容重之間的關(guān)系，是一種對土壤容重進行預測和模擬的有效工具。

3 結(jié)論

(1)黃土高原南北樣帶0～20 cm深度土壤容重變異程度為中等變異，20～40 cm為弱變異。

(2)狀態(tài)空間方程的模擬效果優(yōu)于經(jīng)典統(tǒng)計的多元逐步回歸方程和傳遞函數(shù)方程，基于黏粒和砂粒體積分數(shù)、降水量和土地利用的狀態(tài)空間方程對20～40 cm深度土壤容重的模擬效果最好，決定系數(shù)R2在0.9以上，可用于田間條件下區(qū)域土壤容重分布特征的預測。

(3)樣帶不同土層深度土壤容重的影響因素不同，0～10 cm主要為有機碳含量、黏粒和砂粒體積分數(shù)，10～20 cm為有機碳含量、黏粒和砂粒體積分數(shù)和降水量，20～40 cm為黏粒和砂粒體積分數(shù)、降水量和土地利用。

1 鄭紀勇, 邵明安, 張興昌. 黃土區(qū)坡面表層土壤容重和飽和導水率空間變異特征[J]. 水土保持學報, 2004, 18(3):53-56. ZHENG Jiyong, SHAO Ming’an, ZHANG Xingchang. Spatial variation of surface soils bulk density and saturated hydraulic conductivity on slope in loess region[J]. Journal of Soil and Water Conservation, 2004, 18(3):53-56. (in Chinese)

2 ALEXANDER E B. Bulk densities of California soils in relation to other soil properties[J]. Soil Science Society of America Journal, 1980, 44(4):689-692.

3 劉晚茍, 山侖, 鄧西平. 不同土壤水分條件下土壤容重對玉米根系生長的影響[J]. 西北植物學報, 2002, 22(4):831-838. LIU Wan’gou, SHAN Lun, DENG Xiping. Effects of soil bulk density on the growth of maize root system under different water condition[J]. Acta Botanica Boreali-occidentalia Sinica, 2002, 22(4):831-838. (in Chinese)

4 王輝, 王全九, 邵明安. 表層土壤容重對黃土坡面養(yǎng)分隨徑流遷移的影響[J]. 水土保持學報, 2007, 21(3):10-13. WANG Hui, WANG Quanjiu, SHAO Ming’an. Effect of soil bulk density on soil nutrient in runoff from loess slope[J]. Journal of Soil and Water Conservation, 2007, 21(3):10-13. (in Chinese)

5 韓光中, 王德彩, 謝賢健. 中國主要土壤類型的土壤容重傳遞函數(shù)研究[J]. 土壤學報, 2016, 53(1):93-102. HAN Guangzhong,WANG Decai, XIE Xianjian. Pedotransfer functions for prediction of soil bulk density for major types of soils in China[J]. Acta Pedologica Sinica, 2016, 53(1):93-102. (in Chinese)

6 王巧利, 林劍輝, 許彥峰. 基于BP神經(jīng)網(wǎng)絡的土壤容重預測模型[J]. 中國農(nóng)學通報, 2014, 30(24):237-245. WANG Qiaoli, LIN Jianhui, XU Yanfeng. The model for predicting soil bulk density based on the BP neural network[J]. Chinese Agricultural Science Bulletin, 2014, 30(24):237-245. (in Chinese)

7 WANG Yunqiang, SHAO Ming’an, LIU Zhipeng, et al. Prediction of bulk density of soils in the Loess Plateau region of China[J]. Surveys in Geophysics, 2014, 35(2):395-413.

8 BALLAND V, POLLACCO J A P, ARP P A. Modeling soil hydraulic properties for a wide range of soil conditions[J]. Ecological Modelling, 2008, 219(3-4):300-316.

9 KAUR R, KUMAR S, GURUNG H P. A pedo-transfer function (PTF) for estimating soil bulk density from basic soil data and its comparison with existing PTFs[J]. Australian Journal of Soil Research, 2002, 40(5):847-858.

10 HAN Guangzhong, ZHANG Ganlin, GONG Zitong, et al. Pedotransfer functions for estimating soil bulk density in China[J]. Soil Science, 2012, 177(3):158-164.

11 汪懋華. “精細農(nóng)業(yè)”發(fā)展與工程技術(shù)創(chuàng)新[J]. 農(nóng)業(yè)工程學報, 1999, 15(1):1-8. WANG Maohua. Development of precision agriculture and innovation of engineering technologies[J]. Transactions of the CSAE, 1999, 15(1):1-8. (in Chinese)

12 劉志鵬, 邵明安, 王云強. 區(qū)域尺度下黃土高原土壤全鉀含量的空間模擬[J]. 農(nóng)業(yè)工程學報,2012,28(22): 132-140． LIU Zhipeng, SHAO Ming’an, WANG Yunqiang. Spatial simulation of soil total potassium in regional scale for Loess Plateau region[J]. Transactions of the CSAE, 2012, 28(22):132-140. (in Chinese)

13 賈小旭, 邵明安, 魏孝榮, 等. 黃土高原北部草地表層土壤水分狀態(tài)空間模擬[J]. 農(nóng)業(yè)工程學報,2010,26(10):38-44． JIA Xiaoxu, SHAO Ming’an，WEI Xiaorong, et al. State-space simulation of soil surface water content in grassland of northern Loess Plateau[J]. Transactions of the CSAE, 2010, 26(10):38-44. (in Chinese)

14 趙春雷, 邵明安, 賈小旭. 黃土高原北部坡面尺度土壤飽和導水率分布與模擬[J]. 水科學進展, 2014, 25(6):806-815. ZHAO Chunlei, SHAO Ming’an, JIA Xiaoxu. Distribution and simulation of saturated soil hydraulic conductivity at a slope of northern Loess Plateau[J].Advances in Water Science, 2014, 25(6):806-815. (in Chinese)

15 楊勤業(yè), 張伯平, 鄭度. 關(guān)于黃土高原空間范圍的討論[J]. 自然資源學報, 1988, 3(1):9-15. YANG Qinye, ZHANG Boping, ZHENG Du. On the boundary of the Loess Plateau[J]. Journal of Natural Resources, 1988,3(1):9-15. (in Chinese)

16 JIA X, SHAO M, ZHANG C, et al. Regional temporal persistence of dried soil layer along south-north transect of the Loess Plateau, China[J]. Journal of Hydrology, 2015, 528: 152-160.

17 王云強，邵明安. 黃土高原地區(qū)土壤干層的空間分布與影響因素[D]. 楊凌：中國科學院研究生院, 2010. WANG Yunqiang, SHAO Ming’an.Spatial distribution and influence factors of dried soil layers across the Loess Plateau[D]. Yangling: University of Chinese Academy of Sciences, 2010. (in Chinese)

18 FLATMAN G T, YFANTIS A A. Geostatistical strategies for soil sampling: the survey and the census [J]. Environmental Monitoring and Assessment, 1984, 4(4): 335-349.

19 ZHAO C, SHAO M, JIA X, et al. Using pedotransfer functions to estimate soil hydraulic conductivity in the Loess Plateau of China[J]. CATENA, 2016, 143:1-6.

20 WENDROTH O, REUTER H I, CHRISTIAN K K. Predicting yield of barley across a landscape: a state-space modeling approach[J]. Journal of Hydrology, 2003, 272(1): 250-263.

21 TIMM L C, REICHARDT K, OLIVEIRA J C M, et al. State-space approach for evaluating the soil-plant-atmosphere system[R]. Lectures Given at the College on Soil Physics Trieste, 2003： 387-428.

22 SHUMWAY R H, STOFFER D S. An approach to time series smoothing and forecasting using the EM algorithm[J]. Journal of Time Series Analysis, 1982, 3(4): 253-264.

23 WENDROTH O, JURSCHIK P, KERSEBAUM K C, et al. Identifying, understanding, and describing spatial processes in agricultural landscapes—four case studies[J]. Soil & Tillage Research, 2001, 58(3-4): 113-127.

24 WANG Y, SHAO M, LIU Z, et al. Investigation of factors controlling the regional-scale distribution of dried soil layers under forestland on the Loess Plateau, China[J]. Surveys in Geophysics, 2012, 33(2): 311-330.

25 BENITES V M, MACHADO P L O A, FIDALGO E C C, et al. Pedotransfer functions for estimating soil bulk density from existing soil survey reports in Brazil[J]. Geoderma, 2007, 139(1-2):90-97.

26 ZHAO C, SHAO M, JIA X, et al. Particle size distribution of soils (0～500 cm) in the Loess Plateau, China[J]. Geoderma Regional, 2016, 7(3):251-258.

27 KAUR R, KUMAR S, GURUNG H P. A pedo-transfer function (PTF) for estimating soil bulk density from basic soil data and its comparison with existing PTFs[J]. Australian Journal of Soil Research, 2002, 40(5):847-858.

28 JIA X, SHAO M，WEI X, et al. Estimating total net primary productivity of managed grasslands by a state-space modeling approach in a small catchment on the Loess Plateau，China[J]. Geoderma, 2011, 160(3-4):281-291.

Variation Analysis and Simulation of Soil Bulk Density within Different Soil Depths along South to North Transect of Loess Plateau

YI Xiaobo1SHAO Ming’an1,2ZHAO Chunlei1ZHANG Chencheng1,2

(1.CollegeofResourcesandEnvironment,NorthwestA&FUniversity,Yangling，Shaanxi712100,China2.StateKeyLaboratoryofSoilErosionandDrylandFarmingontheLoessPlateau,NorthwestA&FUniversity,Yangling，Shaanxi712100,China)

To provide soil hydraulic parameters for simulating soil hydrological processes, classical statistics were used to characterize the spatial distribution of soil bulk density (BD) within different soil depths (0～10 cm, 10～20 cm and 20～40 cm) along south to north transect of the Loess Plateau. Furthermore, the multiple stepwise regressions, pedotransfer functions and the first order autoregressive state-space models were applied to simulate the distribution of BD. The magnitude of BD variability in the 0～20 cm soil layer was moderate according to the coefficient of variations, while the BD variability in the 20～40 cm soil layer was weak. The key factors affecting the spatial distribution of BD differed in different soil depths. Soil organic carbon, clay and sand contents were the key factors in the 0～10 cm soil layer; soil organic carbon content, clay and sand contents and precipitation were the key factors in the 10～20 cm soil layer; while in the 20～40 cm soil layer, clay and sand contents, precipitation and land use were the key factors to affect the spatial distribution of BD along south to north transect. State-space models were consistently more effective than multiple stepwise regression functions and pedotransfer functions for estimating spatial distribution of BD. State-space model that included clay and sand contents, precipitation and land use factors showed the best simulation result, and the combination of such variables explained 92.3% of the total variation of BD. State-space models were recommended for studying spatial relations between soil bulk density and other variables on the Loess Plateau.

bulk density; state-space; pedotransfer function; simulation; spatial variation

10.6041/j.issn.1000-1298.2017.04.026

2016-08-14

2016-09-28

國家自然科學基金項目(41530854、41501233)

易小波(1986—)，男，博士生，主要從事環(huán)境生態(tài)研究，E-mail： abobobi@163.com

邵明安(1956—)，男，研究員，博士生導師，主要從事土壤物理研究，E-mail： shaoma@igsnrr.ac.cn

S152.7

A

1000-1298(2017)04-0198-08