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

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

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

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

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

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

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

引言

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 樣品采集

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

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

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

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

1.3 分析方法

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

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

(1)

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

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

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

Xi=φXi-1+wi

(2)

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

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

wi——模型誤差

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

Yi=AiXi+vi

(3)

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

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

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

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

(4)

式中Zi——標(biāo)準(zhǔn)化前的變量zi——標(biāo)準(zhǔn)化后的變量m——標(biāo)準(zhǔn)化前變量的平均值

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

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

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

圖2 不同土層深度土壤容重、有機(jī)碳含量、黏粒和砂粒體積分?jǐn)?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水平上具有顯著性(雙尾檢驗(yàn))；** 表示相關(guān)關(guān)系在p<0.01水平上具有顯著性(雙尾檢驗(yàn))。

為了提高統(tǒng)計(jì)方程的模擬精度，進(jìn)一步將與土壤容重顯著相關(guān)的變量以及其4種常見轉(zhuǎn)換(指數(shù)、對(duì)數(shù)、倒數(shù)和平方轉(zhuǎn)換)作為自變量再次進(jìn)行逐步回歸分析，以期得到對(duì)土壤容重空間變異具有更強(qiáng)解釋能力的傳遞函數(shù)方程，并選出各土層模擬精度最優(yōu)的傳遞函數(shù)方程(表3)。不同土層深度傳遞函數(shù)方程均經(jīng)過顯著性檢驗(yàn)，且均在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%。這主要是因?yàn)閭鬟f函數(shù)方程包含了不同變量的轉(zhuǎn)換，可以更好地刻畫變量與容重之間的關(guān)系[7]。

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

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

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

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

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

3 結(jié)論

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

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

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

1 鄭紀(jì)勇, 邵明安, 張興昌. 黃土區(qū)坡面表層土壤容重和飽和導(dǎo)水率空間變異特征[J]. 水土保持學(xué)報(bào), 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 劉晚茍, 山侖, 鄧西平. 不同土壤水分條件下土壤容重對(duì)玉米根系生長(zhǎng)的影響[J]. 西北植物學(xué)報(bào), 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 王輝, 王全九, 邵明安. 表層土壤容重對(duì)黃土坡面養(yǎng)分隨徑流遷移的影響[J]. 水土保持學(xué)報(bào), 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 韓光中, 王德彩, 謝賢健. 中國(guó)主要土壤類型的土壤容重傳遞函數(shù)研究[J]. 土壤學(xué)報(bào), 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)絡(luò)的土壤容重預(yù)測(cè)模型[J]. 中國(guó)農(nóng)學(xué)通報(bào), 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 汪懋華. “精細(xì)農(nóng)業(yè)”發(fā)展與工程技術(shù)創(chuàng)新[J]. 農(nóng)業(yè)工程學(xué)報(bào), 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 劉志鵬, 邵明安, 王云強(qiáng). 區(qū)域尺度下黃土高原土壤全鉀含量的空間模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào),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è)工程學(xué)報(bào),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 趙春雷, 邵明安, 賈小旭. 黃土高原北部坡面尺度土壤飽和導(dǎo)水率分布與模擬[J]. 水科學(xué)進(jìn)展, 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]. 自然資源學(xué)報(bào), 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 王云強(qiáng)，邵明安. 黃土高原地區(qū)土壤干層的空間分布與影響因素[D]. 楊凌：中國(guó)科學(xué)院研究生院, 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

國(guó)家自然科學(xué)基金項(xiàng)目(41530854、41501233)

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

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

S152.7

A

1000-1298(2017)04-0198-08