方 兵,吳思聰,陳弘揚(yáng),宋 強(qiáng),莊紅娟,周鵬飛,楊 斌,張世文*

(1.安徽理工大學(xué)地球與環(huán)境學(xué)院，安徽 淮南 232001；2.中國科學(xué)院南京土壤研究所，江蘇 南京 210000；3.安徽理工大學(xué)空間信息與測繪工程學(xué)院，安徽 淮南 232001)

【研究意義】土壤容重(bulk density，BD)作為土壤物理性質(zhì)之一，對于評估土壤養(yǎng)分、區(qū)域碳儲(chǔ)量、土壤緊實(shí)度和土壤水分運(yùn)移等至關(guān)重要[1]，被廣泛作為各種預(yù)測和描述性土壤模型的輸入?yún)?shù)。盡管土壤容重?cái)?shù)據(jù)十分重要，但在實(shí)際的獲取工作中，常常由于土壤中植物根系和礫石較多等原因，很難或無法通過環(huán)刀法采樣來測定[2]；此外，由于環(huán)刀法操作繁瑣，在野外大量獲取容重?cái)?shù)據(jù)是一件費(fèi)時(shí)費(fèi)力、甚至不切實(shí)際的工作[3]。因此在土壤科學(xué)領(lǐng)域，一些國內(nèi)外學(xué)者開始利用日常易測定的土壤理化性質(zhì)對土壤容重等不易測定的理化性質(zhì)進(jìn)行模擬，這種更易獲取和實(shí)現(xiàn)的數(shù)學(xué)模擬方法被統(tǒng)稱為土壤轉(zhuǎn)換函數(shù)(pedotransfer functions，PTFs)。【前人研究進(jìn)展】國內(nèi)外用于構(gòu)建土壤轉(zhuǎn)換函數(shù)模型的方法眾多，從基本統(tǒng)計(jì)方法到多變量空間統(tǒng)計(jì)均有涉及。采用逐步回歸確定轉(zhuǎn)換函數(shù)模型，并使用驗(yàn)證集法(以均方根誤差RMSE作為誤差度量指標(biāo))對模型測試誤差進(jìn)行估計(jì)是最常被使用的建模流程[3-5]。如WANG等[5]在中國黃土高原地區(qū)，從1254個(gè)土壤樣品中隨機(jī)抽取1003個(gè)土樣作為模型訓(xùn)練集使用逐步多元回歸得到轉(zhuǎn)換函數(shù)模型，再將余下的251個(gè)土樣作為測試集，以RMSE作為誤差度量指標(biāo)來評價(jià)模型的優(yōu)劣。但逐步回歸因追尋高效的運(yùn)算效率從而無法保證在全部的變量組合模型中得到最優(yōu)變量組合；同時(shí)，用以評判模型優(yōu)劣的驗(yàn)證集法，對數(shù)據(jù)訓(xùn)練集和測試集的分割存在隨機(jī)性，使模型測試誤差的估計(jì)波動(dòng)較大，往往高估了模型的預(yù)測能力，進(jìn)而降低了模型的普適性與預(yù)測精度。近年來許多學(xué)者引入了神經(jīng)網(wǎng)絡(luò)的方法來提升模型輸出的穩(wěn)健性與預(yù)測精度，如高如泰等[6]基于BP-神經(jīng)網(wǎng)絡(luò)(BP-ANN)以土壤粒徑分布、土壤容重、土壤有機(jī)質(zhì)等土壤基本理化性質(zhì)來對土壤水力學(xué)參數(shù)進(jìn)行預(yù)測并得到較好的結(jié)果。廖凱華等[7]利用主成分分析(PCA)與神經(jīng)網(wǎng)絡(luò)相結(jié)合的方法(PANN)，利用PANN消除了ANN輸入層參數(shù)的相關(guān)性，降低了網(wǎng)絡(luò)的拓?fù)鋸?fù)雜度，從而使模型的預(yù)測能力得到提升。神經(jīng)網(wǎng)絡(luò)雖在土壤變量預(yù)測上提高了模型輸出的穩(wěn)健性與預(yù)測精度，但該方法不能有效刻畫出土壤轉(zhuǎn)換函數(shù)的具體形式，進(jìn)而影響了其在實(shí)踐中的應(yīng)用與推廣?！颈狙芯壳腥朦c(diǎn)】針對當(dāng)前各種建模方法的不足，本文采用融合十折交叉檢驗(yàn)的最優(yōu)子集和lasso壓縮估計(jì)的方法，選取影響土壤容重(BD,mg/m3)變化的環(huán)境因子，包括土壤有機(jī)質(zhì)(SOM,g/kg)、砂粒(Sand,%)、黏粒(Clay,%)和粉粒(Silt,%)的體積分?jǐn)?shù)、土壤含水量(WC,%)與土壤采樣深度(Depth,cm)作為預(yù)測變量，分不同垂直尺度(0～10、10～20、20～40、0～40 cm)對皖北平原上的土壤容重進(jìn)行模擬，并與現(xiàn)有土壤容重轉(zhuǎn)換函數(shù)進(jìn)行精度比較，同時(shí)探討了最優(yōu)子集和lasso壓縮估計(jì)在不同垂直尺度下的適用范圍。【擬解決的關(guān)鍵問題】本研究通過探討不同建模方法的優(yōu)劣，選取最優(yōu)建模方法，以期為類似區(qū)域的土壤容重研究提供方法支撐。

1 材料與方法

1.1 研究區(qū)概括與數(shù)據(jù)處理

研究區(qū)位于安徽省北部(圖1)，屬于皖北平原區(qū)域，全境屬暖溫帶半濕潤季風(fēng)氣候，并具有以暖溫帶向北亞熱帶漸變的過渡帶氣候特征。研究區(qū)北部與黃河決口扇形地相連，南部與江淮丘崗區(qū)隔淮河相望，大部屬平原地帶，地勢平坦，總面積達(dá)3.9萬km2。皖北平原年降雨量820～950 mm，境內(nèi)河流均屬淮河水系。土壤類型以砂姜黑土為主，并包含少量潮土，土層厚度在100～120 cm。

為建立區(qū)域容重轉(zhuǎn)換模型，對研究區(qū)利用網(wǎng)格布點(diǎn)與分層抽樣相結(jié)合的方法進(jìn)行土壤樣品的分層(0～10、10～20、20～30和30～40 cm)采集工作。采樣時(shí)每個(gè)樣點(diǎn)隨機(jī)設(shè)置3個(gè)重復(fù)，共設(shè)采集樣點(diǎn)34個(gè)，土樣816個(gè)，其中環(huán)刀土樣和混合土樣各408個(gè)。土壤容重與土壤含水量采用環(huán)刀法測定，105 ℃下烘至恒重后，稱重計(jì)算；土壤有機(jī)質(zhì)采用重鉻酸鉀—外加熱法測定；土壤機(jī)械組成采用比重計(jì)法，依據(jù)土壤礦質(zhì)顆粒粒徑分級標(biāo)準(zhǔn)(美國制)：砂粒(2.000～0.050 mm)，粉粒(0.050～0.002 mm)，黏粒(<0.002 mm)[8]。

1.2 研究方法

1.2.1 十折交叉檢驗(yàn) 十折交叉驗(yàn)證(Ten-fold cross-validation)是對一種統(tǒng)計(jì)學(xué)方法測試誤差的估計(jì)。本文中，交叉驗(yàn)證法被用來評價(jià)最優(yōu)子集法中不同變量規(guī)模下的模型表現(xiàn)能力，以及為lasso壓縮估計(jì)的調(diào)節(jié)參數(shù)Lambda選擇一個(gè)合適的值，用來均衡模型的方差—偏差，提升模型的穩(wěn)健性。

十折交叉方法的具體實(shí)施過程，將土壤樣品的數(shù)據(jù)集隨機(jī)等分成10份，輪流將其中的9份用于訓(xùn)練數(shù)據(jù)集，剩下的1份作為檢驗(yàn)數(shù)據(jù)集，并用殘差平方和(residual sum of squares,RSS)與均方根誤差(root mean squared error,RMSE)作為不同模型下的評價(jià)指標(biāo)。

(1)

(2)

1.2.2 最優(yōu)子集法 最優(yōu)子集法(optimal subset method)對研究區(qū)p個(gè)土壤屬性：土壤有機(jī)碳、土壤砂粒、黏粒和粉粒的體積分?jǐn)?shù)、土壤含水量和土壤采樣深度作為預(yù)測變量的所有可能組合，進(jìn)行擬合回歸。該算法對含有一個(gè)預(yù)測變量的模型(p|1)，擬合p個(gè)模型，并以RSS最小作為(p|1)的最優(yōu)模型；對含有兩個(gè)預(yù)測變量的模型(p|2)，擬合p(p-1)/2個(gè)模型，同時(shí)也是以RSS最小為依據(jù)，選取(p|2)的最優(yōu)模型，依次類推，最終在這p個(gè)模型中，以十折交叉驗(yàn)證法(RSS作為度量指標(biāo))，選取一個(gè)測試誤差最低的模型作為最優(yōu)模型。

圖1 研究區(qū)位置與采樣點(diǎn)分布Fig.1 Location of the study area and distribution of sampling points

為了提高土壤容重轉(zhuǎn)換函數(shù)模型的精度，本文引入最優(yōu)子集法得到最優(yōu)的p個(gè)模型，并采用十折交叉驗(yàn)證規(guī)避了訓(xùn)練集和驗(yàn)證集在分割上存在隨機(jī)性，以殘差平方和RSS作為度量指標(biāo)，評價(jià)這p個(gè)模型的表現(xiàn)能力，并找出預(yù)測精度較高，模型變量組合最優(yōu)的土壤容重轉(zhuǎn)換函數(shù)模型。

1.2.3 lasso壓縮估計(jì)法 在土壤容重轉(zhuǎn)換函數(shù)模型的研究中，簡單線性擬合、多元逐步回歸法都是基于最小二乘法得以實(shí)現(xiàn)，該方法通常保證了所選數(shù)據(jù)訓(xùn)練集的偏差值較小，但無法降低模型方差的大小，即模型對于給定數(shù)值的輸出穩(wěn)定性較差。為了克服最小二乘法的缺陷，國內(nèi)外學(xué)者引入了BP神經(jīng)網(wǎng)絡(luò)的方法來提升模型的穩(wěn)健性[3,9-10]，雖在土壤變量預(yù)測時(shí)提高了模型的穩(wěn)健性與預(yù)測精度，但該方法并不能有效刻畫出土壤轉(zhuǎn)換函數(shù)的具體形式，進(jìn)而影響了其在實(shí)踐中的應(yīng)用與推廣。綜上，本文引入lasso回歸，以犧牲小部分模型偏差為代價(jià)，從而提高模型的穩(wěn)健性，并給出具體的模型刻畫形式。

與最小二乘法不同，lasso法并不具有尺度不變性，所以需要事先對訓(xùn)練集數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，如下式：

(3)

lasso回歸系數(shù)估計(jì)值通過最小化下式得到：

(4)

1.2.4 數(shù)據(jù)處理 土壤理化性質(zhì)的描述性統(tǒng)計(jì)、Pearson相關(guān)分析、LSD多重比較以及基于融合十折交叉檢驗(yàn)的最優(yōu)子集法、lasso壓縮估計(jì)法的模型構(gòu)建及其相關(guān)圖片的制作均使用R-3.5.3平臺(tái)進(jìn)行。

2 結(jié)果與分析

2.1 土壤理化性質(zhì)的統(tǒng)計(jì)特征

2.1.1 描述性統(tǒng)計(jì) 土壤容重含量介于1.33～1.72 mg/m3(表1)；土壤有機(jī)質(zhì)介于1.28～44.62 g/kg，土壤含水量介于9.12%～35.10%，砂粒含量介于18.92%～84.83%，粉粒含量介于13.42%～74.27%，黏粒含量介于0.67%～12.59%。在土層0～40 cm下，土壤容重與其它土壤理化性質(zhì)相比(如：土壤有機(jī)質(zhì)、土壤黏粒)，其變異系數(shù)最小，處于低變異水平(CV=0.08%)；土壤含水量、粉粒和黏粒所選數(shù)據(jù)達(dá)到中等變異水平(CV=19.10%、24.13%和28.43%)，土壤有機(jī)質(zhì)與黏粒數(shù)據(jù)達(dá)到高度變異水平(CV=37.90%、52.96%)。土壤有機(jī)質(zhì)、黏粒含量處于高度變異水平可能與研究區(qū)內(nèi)各采樣小區(qū)上所從事的不同農(nóng)業(yè)活動(dòng)有關(guān)[11-12]。

表1 土壤理化性質(zhì)的描述性統(tǒng)計(jì)(n=134)

集中趨勢和離散趨勢是數(shù)據(jù)分布的2個(gè)重要特征，峰度(skewness)與偏態(tài)(kurtosis)可以通過描述數(shù)據(jù)分布形狀的對稱性、偏斜程度以及扁平程度反應(yīng)這些分布特征，并度量所選數(shù)據(jù)分布與正態(tài)分布的差異。所選土壤容重?cái)?shù)據(jù)服從正態(tài)分布，其頻數(shù)分布圖顯示出輕微左偏和平頂(Skew.=-0.04，Kur.=-0.48)；其它土壤理化性質(zhì)數(shù)據(jù)也都較好地服從正態(tài)分布(表2)。

2.1.2 垂直變化規(guī)律 土壤有機(jī)質(zhì)的含量在表層達(dá)到最大值(圖2)，并隨土層厚度的增加而降低，其含量為20.86、15.63、11.82和9.90 g/kg，土壤容重呈現(xiàn)出先上升后降低的趨勢，并在土層深度20 cm處達(dá)到最大值，為1.48、1.56、1.54和1.52 mg/m3。有研究表明，土壤容重與土壤有機(jī)質(zhì)在土壤表層呈極顯著負(fù)相關(guān)關(guān)系，即隨著土層厚度的增加，土壤有機(jī)質(zhì)含量的逐漸減少，土壤容重則呈現(xiàn)出增大的趨勢[13]。土壤含水量與容重的變化趨勢呈現(xiàn)出相反的變化，即先降低后升高，22.81%、20.19%、21.09%與22.60%；所選土壤黏粒、粉粒與砂粒的含量的變化趨勢較為平穩(wěn)，均在10%的幅度內(nèi)波動(dòng)(5.42%～6.29%、47.34%～48.92%和44.98%～47.28%)。綜上，土壤有機(jī)質(zhì)、土壤含水量、土壤容重的含量在表層0～20 cm變化較大，在土層20～40 cm變化較小，因此在后續(xù)研究中將垂直尺度劃分為0～10、10～20、20～40與0～40 cm。

圖2 土壤理化性質(zhì)的垂直變化規(guī)律Fig.2 Vertical variation of soil physical and chemical properties

BD.土壤容重(mg/m3);WC.土壤含水量(%);SOM.土壤有機(jī)質(zhì)(g/kg);Silt.粉粒(%);Sand.砂粒(%);Clay.黏粒(%) BD.Soil bulk density(mg/m3);WC.Soil water content(%);SOM.Soil organic matter(g/kg);Silt.Silt(%);Sand.Sand(%);Clay.Clay(%)圖3 土壤屬性與所選容重之間的相關(guān)系數(shù)Fig.3 Correlation coefficient between soil properties and selected bulk density

2.1.3 土壤容重的相關(guān)性分析 Curtis、Aleaxander與Benmoux等[3,14-16]認(rèn)為土壤有機(jī)質(zhì)、土壤質(zhì)地是影響土壤容重變化的重要環(huán)境因子；土壤含水量在Suuster等[17]開發(fā)的預(yù)測模型中具有顯著的統(tǒng)計(jì)學(xué)意義，并與容重呈反比；Qiao等[18]在構(gòu)建轉(zhuǎn)換函數(shù)模型時(shí)證明了土壤采樣深度在對容重變異性的解釋占據(jù)重要作用。根據(jù)現(xiàn)有研究，本文選取了上述與土壤容重相關(guān)的土壤屬性，對所獲得的土壤樣本分不同垂直尺度進(jìn)行了Pearson相關(guān)性分析(土壤采樣深度為離散性變量，不參與相關(guān)性分析，圖3)。

從圖3-a的下三角矩陣中看出，土壤容重與其它土壤各屬性的散點(diǎn)分布呈現(xiàn)出較高的非線性關(guān)系；上三角矩陣顯示各土壤屬性之間的Pearson線性相關(guān)系數(shù)，其中土壤容重與土壤有機(jī)質(zhì)呈現(xiàn)出顯著的負(fù)相關(guān)性水平(corr.=-0.17*)，這與王霖嬌等[16]的研究結(jié)果一致，Yun qiang W[23]指出土壤有機(jī)質(zhì)含量可以解釋土壤容重81.00%的變異性，是改變土壤空隙結(jié)構(gòu)的重要聚合體，同時(shí)具有較強(qiáng)的吸水性，有機(jī)質(zhì)含量的增加直接影響了土壤容重的下降[27]；容重與其它土壤屬性的相關(guān)系數(shù)介于-0.07～0.20，均處于較低的線性相關(guān)性，與下三角矩陣所得結(jié)論相一致。

從圖3-b可知，0～10 cm土層中，土壤容重與土壤有機(jī)質(zhì)呈現(xiàn)出顯著正相關(guān)線性關(guān)系(corr.=0.33*)，其系數(shù)的絕對值較其它土壤屬性相比達(dá)到最高；20～40 cm土層中，土壤容重與黏粒呈現(xiàn)出顯著的弱線性關(guān)系(corr.=-0.25*)；0～40 cm土層中，土壤容重與土壤有機(jī)質(zhì)和土壤含水量呈現(xiàn)顯著的弱線性關(guān)系(corr.=-0.17*,-0.20*)。

Yang[13]指出土壤容重與土壤其他屬性間可能存在顯著的非線性關(guān)系。為使研究區(qū)土壤容重的預(yù)測精度得到提高，在后續(xù)模型的構(gòu)建中，將會(huì)對現(xiàn)有數(shù)據(jù)中的一些土壤預(yù)測變量進(jìn)行4種常見的轉(zhuǎn)換，如：lnSOM、SOM2、(lnSOM)2和1/SOM。

2.2 模型的建立及其精度比較

2.2.1 土壤容重的最優(yōu)子集回歸模型 由圖4可知,土層厚度為0～10和10～20 cm時(shí)，當(dāng)變量數(shù)個(gè)數(shù)為1時(shí)(分別為：1/SOM、ln(SOM))，模型殘差平方和值(RSS=0.005、0.002)達(dá)到最低。土層厚度為20～40 cm時(shí)，當(dāng)所選3個(gè)變量(SOM、WC和Sand)時(shí)，殘差平方和(RSS=0.005)達(dá)到最低。土層厚度為0～40 cm時(shí)，其所選2個(gè)變量[ln(SOM)、lnWC]使模型普適性達(dá)到最佳(RSS=0.006)。

4種函數(shù)模型精度均較高(RMSE接近0，圖5)，從散點(diǎn)分布趨勢來看，土層0～10、10～20和20～40 cm，容重預(yù)測值與實(shí)際值較為接近1∶1軸，即預(yù)測結(jié)果較好，但土層0～40 cm，容重預(yù)測值與實(shí)際值偏離1∶1軸較大；從土壤容重轉(zhuǎn)換函數(shù)的變量形式來看，易小波等[28]利用狀態(tài)空間方程估算土壤容重值時(shí)提出，不同土層深度下土壤容重的影響因素不同，在本文得以印證。土層0～20 cm時(shí)土壤有機(jī)質(zhì)是其主要的影響因素；而土層20～40 cm時(shí)，土壤有機(jī)質(zhì)、含水量和砂粒是其共同影響因素。土壤有機(jī)碳含量隨土層厚度的加深而逐漸減小，也是影響函數(shù)模型參數(shù)的重要因素，在土壤有機(jī)質(zhì)含量較高的土壤中，容重主要受有機(jī)質(zhì)含量影響較大，而在有機(jī)質(zhì)含量較低時(shí)，土壤質(zhì)地等其它理化性質(zhì)才對容重產(chǎn)生影響[4,19-20]。

土層厚度:(a).0～10 cm;(b).10～20 cm;(c).0～40 cm;(d).20～40 cm Soil thickness:(a).0-10 cm;(b).10-20 cm;(c).0-40 cm;(d).20-40 cm圖4 不同土層深度下的十折交叉驗(yàn)證結(jié)果Fig.4 Ten-fold cross-validation results for different soil depths

圖5 基于最優(yōu)子集法的土壤容重預(yù)測值與實(shí)測值比較Fig.5 Comparison of predicted values of soil bulk density and measured values based on optimal subset method

土層厚度:(a).0～10 cm;(b).10～20cm;(c).0～40 cm;(d).20～40 cm Soil thickness:(a).0-10 cm;(b).10～20cm;(c).0～40 cm;(d).20～40 cm圖6 基于十折交叉驗(yàn)證的Lambda最優(yōu)選擇Fig.6 Lambda optimal selection based on ten-fold cross-validation

2.2.2 土壤容重的lasso壓縮估計(jì)模型 Lambda值選擇在-102～1010的范圍進(jìn)行土壤轉(zhuǎn)換函數(shù)形式的壓縮估計(jì)，該范圍包含了只含截距項(xiàng)的空模型到最小二乘估計(jì)的擬合模型的所有情況，并利用十折交叉檢驗(yàn)法計(jì)算每個(gè)值的交叉驗(yàn)證誤差，以RMSE作為評價(jià)指標(biāo)，選擇誤差最小的值。在土層厚度0～10、10～20、20～40和0～40 cm下最優(yōu)的Lambda=33.442、0.010、0.011和0.006(圖6)；除土層0～10 cm的Lambda對回歸系數(shù)估計(jì)的影響程度較大，對其余值的影響都較小。

圖7 基于lasso壓縮估計(jì)法的土壤容重預(yù)測值與實(shí)測值比較Fig.7 Comparison of soil bulk density predictions and measured values based on lasso compression estimates

將上述所得Lambda的最優(yōu)值帶入lasso模型進(jìn)行各土層下的轉(zhuǎn)換函數(shù)模型的建立(圖7)。lasso模型在土層0～10與10～20 cm下的土壤容重的預(yù)測精度較低，RMSE分別為0.216、0.091 mg/m3，較最優(yōu)子集法相比分別下降了67%和40%；從散點(diǎn)分布趨勢可以看出，容重預(yù)測值的變化范圍較小(1.485～1.487、1.553～1.572 mg/m3)，與土層的容重實(shí)際值的分布偏離較大。土層20～40與0～40 cm下的土壤容重預(yù)測精度(RMSE=0.007、0.029 mg/m3)較最優(yōu)子集法相比得到了較大的提升，分別為90%和53%；從散點(diǎn)分布趨勢來看，散點(diǎn)整體趨勢較最優(yōu)子集法相比，更加傾向于1∶1軸，即土壤容重預(yù)測值與實(shí)際值偏離較小，特別是在土層20～40 cm偏離達(dá)到最低。

2.2.3 兩種土壤容重轉(zhuǎn)換函數(shù)預(yù)測精度的比較分析 對于土壤容重的實(shí)際值、最優(yōu)子集法與lasso法所得土壤容重預(yù)測值進(jìn)行描述性統(tǒng)計(jì)與LSD多重比較(表2)。

由土壤容重的最大值、最小值以及變異系數(shù)可知，在土層0～10、10～20 cm，lasso法預(yù)測所得最大值偏低，最小值偏高，變異系數(shù)較低，由此說明lasso法在模擬土壤容重時(shí)表現(xiàn)出一定的平滑效應(yīng)，而在土層20～40、0～40 cm下，最優(yōu)子集法較lasso法相比平滑效應(yīng)更為明顯；LSD多重比較顯示兩種方法所得土壤容重預(yù)測值之間不存在顯著性的差異，其與土壤容重的實(shí)際值也不存在差異(P>0.05)，從統(tǒng)計(jì)學(xué)角度來看，上述2種方法對于土壤容重的預(yù)測值與實(shí)際值在數(shù)值上相差不大，即均為正確的土壤容重預(yù)測值。但綜合平均值和標(biāo)準(zhǔn)差，土層0～10、10～20 cm，最優(yōu)子集法較lasso法預(yù)測效果更佳，土層20～40、0～40 cm，lasso法預(yù)測值則更為接近實(shí)際土壤容重的分布趨勢。綜上，土層0～10，10～20 cm，使用最優(yōu)子集法較lasso法效果更好；土層20～40、0～40 cm，使用lasso法模擬土壤容重更為精確。

2.2.4 不同土壤容重轉(zhuǎn)換函數(shù)在皖北平原表現(xiàn)能力比較 由于成土過程、氣候、生物、母質(zhì)以及人類活動(dòng)等影響，使土壤自身存在廣泛的空間異質(zhì)性的特點(diǎn)，所以土壤轉(zhuǎn)換函數(shù)具有一定的研究區(qū)和適用限制范圍；韓光中等[3]曾對中國的主要土壤類型分別建立了最適宜的土壤轉(zhuǎn)換函數(shù)，結(jié)果顯示不同土壤類型所具有的轉(zhuǎn)換函數(shù)在形式上與變量上均存在較大差異。本文從現(xiàn)已發(fā)表的國內(nèi)外土壤容重轉(zhuǎn)換函數(shù)模型中，以存在土壤有機(jī)質(zhì)、土壤質(zhì)地、土壤含水量和土壤采樣深度為參數(shù)的選取依據(jù)，結(jié)合本文已有土壤數(shù)據(jù)集，按照相應(yīng)條件代入各傳遞函數(shù)模型中，求得土壤容重的預(yù)測值，并以均方根誤差(RMSE)與平均誤差(ME)為評價(jià)指標(biāo)來衡量在皖北平原區(qū)域上土壤容重預(yù)測的適用性潛力。

表2 2種方法下的土壤容重?cái)?shù)據(jù)的統(tǒng)計(jì)分析

由表3可知，各土壤轉(zhuǎn)換函數(shù)在研究區(qū)的適用性表現(xiàn)不一。從平均誤差值來看，利用Curtis R O(1964)、Adams W A(1973)和Benites V M(2007)模擬得到的ME值為正，說明這些容重傳遞函數(shù)模型高估了容重值；其余函數(shù)模型的ME值均為負(fù)值，說明了函數(shù)模型低估了容重值。通過RMSE作為另一個(gè)評價(jià)指標(biāo)，進(jìn)一步比較現(xiàn)有土壤容重函數(shù)的適用性潛力。所選函數(shù)模型的RMSE值在0.121～24.602 mg/m3，波動(dòng)幅度較大，其中Adams W A(1973)、Yang Y Q (2005)、HAN G Z(2016)-A與HAN G Z(2016)-C的RMSE值最小，在0.012～0.034 mg/m3。綜合RMSE與ME值來看，HAN G Z(2016)-A模型模擬結(jié)果相對最好，該模型的RMSE值接近于0，ME值為-1.197 mg/m3，研究區(qū)土壤有機(jī)質(zhì)含量較高，可能是其模擬結(jié)果較好的原因。Huntington(1989)與Benites V M(2007)模型模擬結(jié)果最差，其原因可能在于，研究區(qū)土壤生成環(huán)境的不同，Huntington的研究區(qū)位于新罕布什爾州一個(gè)23 hm2流域上的森林土壤區(qū)域，皖北平原區(qū)域特殊的地形地貌(淮河平原區(qū))和氣候背景(暖溫帶半濕潤季風(fēng)氣候)與其差異較大。Benite等[4]曾就不同土壤母質(zhì)生成環(huán)境對于現(xiàn)已發(fā)表的土壤容重函數(shù)的使用發(fā)出提醒。

在土層0～40 cm，與現(xiàn)存該地區(qū)最優(yōu)模型HAN G Z(2016)-A相比，使用融合十折交叉檢驗(yàn)的最優(yōu)子集、lasso壓縮估計(jì)所開發(fā)的函數(shù)模型精度從RMSE上來看提高了一個(gè)小數(shù)位(RMSE=0.121至0.063 mg/m3)，這說明了本文所構(gòu)建的模型相比于現(xiàn)存模型，其更好的適用于皖北平原地區(qū)。

3 討 論

最優(yōu)子集法和lasso法分別在不同垂直尺度下對土壤容重的集中與分散趨勢進(jìn)行了較好的預(yù)測。但從散點(diǎn)分布趨勢中可以看出，模型中依然存在一些并未解釋的變異性。除了一些數(shù)據(jù)本身無法解釋的變異外，在模型中引入對時(shí)空要素的刻畫可能會(huì)對模型的解釋能力有顯著的改善。Lark和Cullis等[21-23]利用采樣點(diǎn)的坐標(biāo)信息與克里金法對模型施加了空間相關(guān)結(jié)構(gòu)，使模型的解釋能力得到提升。段良霞[24]則以時(shí)間序列數(shù)據(jù)相互依賴為基礎(chǔ)，通過空間狀態(tài)模型和經(jīng)典線性回歸模型對該區(qū)不同尺度下的土層深度水含量進(jìn)行預(yù)測，并得到了較好的預(yù)測效果。土壤各屬性間的相關(guān)性大小，通常表現(xiàn)為單一剖面下采樣數(shù)據(jù)高于在不同剖面下所獲得的數(shù)據(jù)[25]，這種存在于單一剖面下的特有的相關(guān)性也許可以反映出該剖面對于其內(nèi)部各土壤屬性所產(chǎn)生的各種微小、未知的影響。同時(shí)，這種特有的相關(guān)性也可能廣泛存在于單一的采樣區(qū)、研究區(qū)，甚至包括特定的植被類型、氣候與年份下[16]。這種特有的相關(guān)性是很難被實(shí)際捕獲或量化。若只特定研究某一剖面下土壤容重的轉(zhuǎn)換函數(shù)，這種難以量化的特有的相關(guān)性則可以作為噪音進(jìn)行處理，也即在一般回歸模型下包含固定效應(yīng)(各土壤屬性對于容重的影響)和噪音(特定剖面對土壤各屬性的影響)。但現(xiàn)有國內(nèi)外研究[18,26-27]中，往往是基于大量不同土壤剖面下所采集的數(shù)據(jù)，此時(shí)的數(shù)據(jù)集中存在兩種隨機(jī)因素會(huì)影響模型的預(yù)測能力，一種是某個(gè)特定土壤剖面所具有的隨機(jī)噪聲，另一種則是因?yàn)槠拭媾c剖面在空間或時(shí)間上的不同而形成的隨機(jī)效應(yīng)，這種存在于數(shù)據(jù)集中嵌套的隨機(jī)因素結(jié)構(gòu)通常會(huì)導(dǎo)致對數(shù)據(jù)中實(shí)際可用信息量的誤判，也可能會(huì)對參數(shù)估計(jì)和假設(shè)檢驗(yàn)的結(jié)果產(chǎn)生影響。

表3 國內(nèi)外土壤容重轉(zhuǎn)換函數(shù)

在土層0～40 cm下的土壤容重預(yù)測模型中，土壤采樣深度在最優(yōu)子集法與lasso法中均未被選入。這是因?yàn)楸韺油寥烙袡C(jī)碳含量與土壤采樣深度間的高度的相關(guān)性(corr.=-0.67**)所造成的模型過擬合而被剔除。同時(shí)，在很多土壤容重的預(yù)測模型的研究中，采樣深度均不被認(rèn)作是用于預(yù)測容重的可能性變量。采樣深度可以用作容重的預(yù)測因子，但它只能描述了其變異性中較少的部分[7-8]。但上述這些研究多集中在土壤表層(0～30 cm)，Qiao等[28]在中國黃土高原地區(qū)進(jìn)行了土層50～200 m下的土壤容重模型的構(gòu)建，模型中采樣深度作為的容重變異性的重要解釋因子，但這可能也與土壤有機(jī)碳在土層深處含量較低有關(guān)，進(jìn)而提高了采樣深度的重要性。在之前的研究中，一些學(xué)者將采樣深度作為連續(xù)型變量代入模型中，并得出容重與采樣深度呈現(xiàn)出顯著的非線性關(guān)系[19,29-30]。由于大多數(shù)的土壤樣本是在固定深度下進(jìn)行采集，因此無需將此結(jié)論擴(kuò)展到任意模型中。盡管已有研究說明，表層采樣深度對于容重的預(yù)測并不起決定性作用，但本文還是建議將采樣深度作為分類變量引入任何的模型策略中，因?yàn)殡S著土層深度的增加，土壤中各屬性均發(fā)生變化。

4 結(jié) 論

對于當(dāng)前研究現(xiàn)狀的不足，本文選取了最優(yōu)子集、lasso壓縮估計(jì)與十折交叉檢驗(yàn)法構(gòu)建土壤容重轉(zhuǎn)換函數(shù)，并探究其在研究區(qū)不同垂直尺度下對土壤容重預(yù)測的效果，獲得以下結(jié)論。

(1)土壤容重?cái)?shù)據(jù)服從正態(tài)分布(Skew.=-0.04，Kur.=-0.48)；土層0～40 cm，土壤容重變異系數(shù)為5.07%，處于低變異水平。Pearson相關(guān)性表明，在不同垂直尺度下，土壤容重與土壤有機(jī)質(zhì)、含水量、黏粒、粉粒、砂粒呈現(xiàn)出較低的線性相關(guān)性(corr.<0.5)。

(2)在土壤有機(jī)質(zhì)含量較高土層中，土壤容重的系統(tǒng)變異性主要受土壤有機(jī)質(zhì)影響，當(dāng)有機(jī)質(zhì)含量較低時(shí)，土壤質(zhì)地、土壤含水量對土壤容重的影響程度得到提升；土壤采樣深度因與有機(jī)質(zhì)間高度的相關(guān)性(corr.=-0.67**)所造成的模型過擬合，固在最優(yōu)子集和lasso壓縮估計(jì)法均被剔除。

(3)基于融合十折交叉檢驗(yàn)的最優(yōu)子集與lasso壓縮估計(jì)所得土壤容重預(yù)測值與實(shí)際值沒有顯著差異(P>0.05)，但在土層0～10、10～20 cm，最優(yōu)子集法較lasso法預(yù)測效果更佳(RMSE=0.070、0.054)，土層20～40、0～40 cm，lasso法預(yù)測值則更為接近實(shí)際的土壤容重分布趨勢(RMSE=0.070、0.029)。

上述兩種方法構(gòu)建的模型相比于現(xiàn)有土壤轉(zhuǎn)換函數(shù)模型，更好的適用于皖北平原地區(qū)，可為類似區(qū)域的土壤容重研究提供方法支撐。