基于最大似然法與矩法的黃土高原小流域土壤碳氮空間變異分析

方 炫，安韶山，薛志婧，李壁成

土壤由于受到氣候、地形、植被、水文、土地利用方式等因子的綜合作用，在空間分布上具有高度異質(zhì)性［1］。土壤特性的空間變異已成為當(dāng)前土壤科學(xué)研究的熱點(diǎn)之一［2］。黃土高原是我國乃至全球土壤侵蝕最為嚴(yán)重的地區(qū)之一，長期以來不合理的土地利用，加劇了土壤侵蝕和生態(tài)環(huán)境惡化［3－4］。退耕還林生態(tài)恢復(fù)實(shí)施以來，該地區(qū)土地利用與土地覆蓋發(fā)生重大變化，土壤特性受到影響。土壤碳氮是土壤特性的重要指標(biāo)，土壤碳氮循環(huán)對全球變化有重要作用［5］。開展退耕后黃土高原土壤碳氮空間分異研究，對揭示土壤碳氮分布與土地利用方式等生態(tài)因子之間的關(guān)系具有重要意義。地統(tǒng)計(jì)方法與經(jīng)典統(tǒng)計(jì)學(xué)相比，能夠更好地將空間格局與生態(tài)過程聯(lián)系起來［6］，有效地揭示空間格局對生態(tài)過程的影響［7］，已被廣泛應(yīng)用于土 壤 碳 氮 空 間 分 布 的 研 究［1，6，8－10］。 估算變異函數(shù)是采用地統(tǒng)計(jì)方法進(jìn)行空間變異分析的關(guān)鍵。國內(nèi)在研究土壤碳氮空間變異中基本上都采用 Matheron［11］提 出 的 矩 法 （method of moments，MOM）求變異函數(shù)，但該方法對研究對象的樣本數(shù)量要求較高。Kerry［12］指出，MOM法要求樣本數(shù)量在采樣間距合理的前提下大于100，才能得到充分描述空間異質(zhì)性的變異函數(shù)。Webster和Oliver［13］研究指出，當(dāng)樣本數(shù)量無法達(dá)到100時(shí)，受限最大似然法（residual maximum likelihood，REML）是可選的變異函數(shù)估算方法。Kerry［12］采用矩法和受限最大似然法對不同樣本數(shù)量下精準(zhǔn)農(nóng)業(yè)的土壤性質(zhì)變異函數(shù)進(jìn)行估算，并將兩種方法結(jié)果進(jìn)行比較研究，結(jié)果表明當(dāng)樣本數(shù)量小于100時(shí)，采用受限最大似然法估計(jì)的變異函數(shù)進(jìn)行空間插值預(yù)測更加準(zhǔn)確?？梢?，REML法是土壤空間變異研究可以采用的方法，而該方法在我國土壤空間變異的研究中還鮮有應(yīng)用。本研究以退耕后的黃土高原寺底溝小流域?yàn)閷ο?，采用受限最大似然法和矩法求變異函?shù)，通過對比選擇最優(yōu)的變異函數(shù)模型進(jìn)行地統(tǒng)計(jì)插值，以此分析小流域尺度下土壤有機(jī)碳氮的空間變異特征，旨在為進(jìn)一步研究土壤碳氮變異的生態(tài)過程以及改善土壤碳氮水平提供科學(xué)依據(jù)。

1 研究區(qū)概況

研究區(qū)寺底溝小流域，位于寧夏自治區(qū)南部固原市原州區(qū)河川鄉(xiāng)內(nèi)，地理范圍為東經(jīng)106°25′13″—106°26′17″，北緯 35°57′34″—35°59′17″，總面積為5.16km2。該小流域地處黃土高原丘陵溝壑區(qū)，海拔1 590～1 815m之間，溝坡地占的88%，坡度小于8°的平緩地僅占12%；屬溫帶半干旱大陸性氣候，年均氣溫與降水量分別是6.9°C和419mm（1982—2002年）；土壤類型主要為黃土母質(zhì)上發(fā)育的黃綿土；本土植被是天然草地，由于過度開墾而一度減少，直到2002年退耕還林后才有所恢復(fù)。退耕前植被類型以耕地和天然草地占主導(dǎo)地位，經(jīng)過退耕和農(nóng)業(yè)結(jié)構(gòu)的調(diào)整，植被覆蓋發(fā)生重大變化，植被類型包括天然草地、人工灌木林〔檸條（C.korshinskii）和沙棘（H.rhamnoides）〕和人工草地〔苜蓿（M.sativa）〕、耕地、果園和棄耕地等。其中，苜蓿與檸條為退耕植被。

2 材料與方法

2.1 土樣采集與分析

小流域尺度下土壤類型和氣候水文條件可視為基本上同質(zhì)，考慮土地利用方式和地形因素，按照不同土地利用方式和不同坡位（坡上、坡中、坡下），于2010年7月對寺底溝小流域進(jìn)行采樣，共選取46個(gè)代表性的土壤樣點(diǎn)（圖1）。各樣點(diǎn)分0—10cm，10—30cm，30—60cm共3層采集土壤樣品，每個(gè)樣點(diǎn)設(shè)5次重復(fù)，5次樣本充分混合后取1kg帶回實(shí)驗(yàn)室。土壤樣品采回后經(jīng)自然風(fēng)干，去除根系石塊，研磨過篩后分析土壤有機(jī)質(zhì)和全氮，有機(jī)質(zhì)采用重鉻酸鉀外加熱法測定，全氮采用凱氏定氮法測定。

2.2 數(shù)據(jù)準(zhǔn)備

樣點(diǎn)選擇所參考的基本圖件為2010年寺底溝的土地利用圖和1982年1∶1萬地形圖，其中土地利用圖（圖1）是由2008年SPOT5高分辨率衛(wèi)星影像人機(jī)交互解譯后，再經(jīng)2010年7月實(shí)地調(diào)查對每個(gè)地塊逐個(gè)更新獲得。采用Kappa指數(shù)對該影像分類精度進(jìn)行定量評估，Kappa指數(shù)為0.89，超過了最低允許判別精度 Kappa指數(shù)≥0.70的要求［14－15］。土地利用圖和地形圖投影坐標(biāo)均統(tǒng)一為阿爾伯斯（Albers）投影和Krasovsky－1940坐標(biāo)系統(tǒng)。采樣點(diǎn)坐標(biāo)由GPS確定導(dǎo)入ArcGIS軟件，原經(jīng)緯度坐標(biāo)經(jīng)投影轉(zhuǎn)換為與上述圖件一致的投影坐標(biāo)系統(tǒng)，形成樣點(diǎn)分布圖。將46個(gè)采樣點(diǎn)的各土層有機(jī)碳、全氮數(shù)據(jù)添加到采樣點(diǎn)屬性表中，以分析土壤碳氮的空間變異特征。

2.3 地統(tǒng)計(jì)方法

本研究采用地統(tǒng)計(jì)方法分析土壤碳氮的空間變異特征，變異函數(shù)用來分析土壤碳氮的空間格局，普通克里金插值方法用來進(jìn)行土壤碳氮空間分布的預(yù)測模擬。估計(jì)變異函數(shù)參數(shù)最常用的方法是矩法，即常用公式（1）計(jì)算得到實(shí)驗(yàn)半方差，然后采用加權(quán)最小二乘法對這些實(shí)驗(yàn)半方差進(jìn)行擬合得到理論變異函數(shù)的各參數(shù)值［12］。變異函數(shù)的理論模型一般包括球面模型、指數(shù)模型、高斯模型和線性模型，最優(yōu)理論模型由反映變異函數(shù)擬合程度的決定系數(shù)大小來選擇。

式中：N（h）——距離等于h時(shí)的點(diǎn)對數(shù)是；Z（xi）——在位置xi處的數(shù)值；Z（xi＋h）——在xi＋h處的數(shù)值。

圖1 寺底溝土地利用與樣點(diǎn)分布

由于本研究樣本量小于100，所以同時(shí)采用受限最大似然法來估計(jì)土壤變量的變異函數(shù)。受限最大似然法是Patterson和Thompson提出的［16］，該方法估計(jì)變異函數(shù)參數(shù)是無偏估計(jì)［17］，比最小二乘法的結(jié)果更加可信［18］。與矩法估計(jì)變異函數(shù)相比，可直接利用原始數(shù)據(jù)進(jìn)行估計(jì)，不需要搜索不同空間距離的數(shù)據(jù)對，避免因局部空間內(nèi)數(shù)據(jù)點(diǎn)對數(shù)太少而使所求變異函數(shù)不可靠［19］。受限最大似然法估計(jì)變異函數(shù)的原理和公式在文獻(xiàn)中［12，19－20］均有敘述。本研究采用Pardo—Iguzquiza［20］提出的 MLREML程序來實(shí)現(xiàn)受限最大似然法（REML）方法的變異函數(shù)計(jì)算。MLREML程序提供了球面模型、指數(shù)模型、高斯模型3種理論模型，并且選擇其中負(fù)對數(shù)似然函數(shù)最小的作為變異函數(shù)最優(yōu)理論模型。

采用矩法和受限最大似然法兩種方法估計(jì)變異函數(shù)參數(shù)，然后通過交叉檢驗(yàn)來選擇其中最優(yōu)的變異函數(shù)［12］。交叉檢驗(yàn)中采用平均誤差（ME）和標(biāo)準(zhǔn)化克里金方差（MSDR）來衡量兩種方法估計(jì)的變異函數(shù)進(jìn)行克里金預(yù)測的效果。平均誤差是指每個(gè)點(diǎn)的實(shí)際觀察值與普通克里金預(yù)測值之間差值的平均值，標(biāo)準(zhǔn)化克里金方差是指每個(gè)點(diǎn)的誤差平方與克里金方差比值的平均值，分別采用公式（2）—（3）求得：

式中：Zxi——克里金預(yù)測值：z（xi）——實(shí)際觀察值。ME值越接近于0，MSDR值越接近于1，表明變異函數(shù)模型越準(zhǔn)確。本研究采用GS＋和Excel軟件進(jìn)行模擬變異函數(shù)模型，利用ArcGIS 9.2進(jìn)行普通克里金插值。

3 結(jié)果分析與討論

3.1 土壤有機(jī)碳、全氮和碳氮比的統(tǒng)計(jì)分析

研究區(qū)土壤各指標(biāo)的描述性統(tǒng)計(jì)與正態(tài)分布檢驗(yàn)結(jié)果詳見表1。土壤表層（0—10cm）平均有機(jī)碳含量和平均全氮含量分別為9.16和0.8g／kg，且都隨著土層的加深而減少。各指標(biāo)的變異系數(shù)均在0.1～0.9之間，表明有機(jī)碳含量和全氮含量都屬于中等變異。從不同土層土壤性質(zhì)的變異程度來看，有機(jī)碳含量的變異程度隨著土層加深略有下降，全氮含量的變異程度略有增加。單樣本柯爾莫哥洛夫—斯米爾諾夫檢驗(yàn)（K—Stest）結(jié)果表明，各土層土壤有機(jī)碳和全氮均呈正態(tài)分布（K—S檢驗(yàn)值均大于0.05），其原始數(shù)據(jù)均可直接用于地統(tǒng)計(jì)分析。

表1 研究區(qū)土壤有機(jī)碳含量和全氮的描述性統(tǒng)計(jì)與正態(tài)分布檢驗(yàn)

3.2 土壤有機(jī)碳和全氮的空間異質(zhì)性

由REML和MOM兩種方法估計(jì)土壤碳氮變異函數(shù)的交叉檢驗(yàn)對比結(jié)果詳見表2。除了10—30cm土壤全氮，其它變量均是采用受限最大似然法求得變異函數(shù)的標(biāo)準(zhǔn)克里金方差（MSDR）更接近于1，平均誤差（MS）更接近于0。此結(jié)果表明大多數(shù)情況下，采用受限最大似然法估計(jì)的變異函數(shù)進(jìn)行克里金預(yù)測，效果優(yōu)于鉅法估計(jì)的變異函數(shù)。因此，除了10—30 cm土層的土壤全氮克里金插值采用MOM法變異函數(shù)外，其它土壤有機(jī)碳和全氮變量均采用REML法估計(jì)的變異函數(shù)參數(shù)進(jìn)行克里金插值。

由土壤有機(jī)碳和全氮變異函數(shù)模型及相關(guān)參數(shù)可以看出（表2），塊金方差表示實(shí)驗(yàn)誤差和小于最小取樣尺度引起的隨機(jī)變異，結(jié)構(gòu)方差代表由空間自相關(guān)引起的結(jié)構(gòu)性變異，塊金效應(yīng)指塊金方差占空間總變異的比例，可以反映空間變異程度的大?。?1］。塊金效應(yīng)小于0.25表示變量具有強(qiáng)烈的空間相關(guān)性，在0.25～0.75之間表示具有中等空間相關(guān)性，大于0.75表示空間相關(guān)弱。由表2可以看出，土壤有機(jī)碳含量和全氮含量的空間相關(guān)性基本上都屬于中等。與全氮含量相比，土壤有機(jī)碳含量的隨機(jī)變異比例較大，結(jié)構(gòu)性變異比例較小，空間相關(guān)性相對較弱。不同土層深度相比較，無論有機(jī)碳還是全氮，都表現(xiàn)為表層土壤比下層土壤隨機(jī)變異比例較大，結(jié)構(gòu)性變異比例較小，空間相關(guān)性較弱。

表2 研究區(qū)土壤有機(jī)碳（SOC）和全氮（TN）變異函數(shù)模型的相關(guān)參數(shù)及交叉檢驗(yàn)結(jié)果

變程大小代表空間相關(guān)性范圍的大小，即變程之內(nèi)變量具有空間相關(guān)性，反之不存在空間相關(guān)性。表2顯示，在0—10cm，10—30cm，30—60cm這3個(gè)不同土層深度下，土壤有機(jī)碳變程分別為952.33，605.24和576.82m，而全氮變程分別為799.93，291和486.01m，說明土壤有機(jī)碳空間相關(guān)性范圍大于全氮。土壤全氮空間相關(guān)尺度較小與氮肥使用和土壤氮素流失有關(guān)。一方面，研究區(qū)地形破碎、土壤疏松，土壤侵蝕現(xiàn)象多發(fā)，土壤養(yǎng)分流失嚴(yán)重，而氮素在土坡中具有易移動(dòng)特性，而使得土壤流失及淋溶造成的土壤氮素?fù)p失較有機(jī)質(zhì)多，從而有機(jī)碳空間相關(guān)性范圍比全氮大［1，22］。另一方面，氮肥是研究區(qū)農(nóng)戶種植作物和果樹使用較多的肥料，氮肥施用差異及產(chǎn)生這種差異的隨機(jī)性可能造成土壤全氮空間相關(guān)尺度的減少。比較不同土層下土壤變量的變程大小可以看出，表層土壤（0—10cm）有機(jī)碳和全氮空間相關(guān)性范圍比下層土壤大，其中土壤有機(jī)碳變異范圍隨著土層加深而減小。

3.3 土壤有機(jī)碳和全氮空間分布特征

土壤有機(jī)碳和全氮的地統(tǒng)計(jì)插值結(jié)果如圖2所示。結(jié)合圖1可以看出，不同土地利用方式下的土壤有機(jī)碳和全氮含量總體表現(xiàn)為灌木林最大、天然草地和棄耕地居中、梯田和人工草地最小。其中，有機(jī)碳含量在0—10cm土層的最大值出現(xiàn)在沙棘林地，20—30cm和30—60cm土層的有機(jī)碳含量最大值出現(xiàn)在檸條林地；土壤全氮含量在0—10cm土層的最大值出現(xiàn)在沙棘林地，10—30cm土層的最大值出現(xiàn)在棄耕地，30—60cm土層的最大值出現(xiàn)在檸條林地和天然草地。分析結(jié)果表明，灌木林、天然草地和棄耕地與其它土地利用方式相比，表現(xiàn)出更好的土壤有機(jī)碳氮蓄存能力，證實(shí)了退耕還林政策從總體上對增加土壤碳氮儲(chǔ)量是有積極作用的。一方面，土地利用方式影響土壤有機(jī)碳氮水平，耕地轉(zhuǎn)化為多年生植物，能夠增加土壤有機(jī)碳含量；相反，當(dāng)耕地轉(zhuǎn)化為多年生植物時(shí)，土壤有機(jī)碳會(huì)以CO2的形式流失進(jìn)入大氣圈［23］。另一方面，耕地施肥沒能提高土壤全氮水平，這與該小流域水土流失嚴(yán)重密切相關(guān)。另外，人工草地作為退耕還林中重要的退耕植被，對土壤碳氮儲(chǔ)量沒有明顯貢獻(xiàn)，這是由于人工草地（苜蓿）作為當(dāng)?shù)剞r(nóng)民畜牧業(yè)和家庭養(yǎng)殖業(yè)的天然飼料而被頻繁收割，導(dǎo)致有機(jī)碳和氮無法在土壤中積累。盡管如此，考慮到人工草地防治水土流失的積極作用以及其作為農(nóng)戶畜牧養(yǎng)殖的主要飼料來源，建議人工草地在控制收割次數(shù)的情況下，繼續(xù)保留作為退耕還林的重要植被類型。

不同土壤深度對土壤有機(jī)碳氮空間分布的影響表現(xiàn)為前者隨著土壤深度變化不明顯，而后者明顯表現(xiàn)為表層（0—10cm）土壤有機(jī)碳含量及變異程度大于下層土壤。表明表層土壤對于有機(jī)碳積累很重要，且土地利用方式對土壤有機(jī)碳匯的作用主要發(fā)生在土壤表層，可能與表層土壤受到植被覆蓋、枯落物和根系的植物量以及人類干擾的影響較大有關(guān)。

圖2 寺底溝土壤有機(jī)碳和全氮空間分布

4 結(jié)論

（1）采用REML法與MOM法分別對不同深度下土壤有機(jī)碳和全氮的變異函數(shù)進(jìn)行估算，并通過交叉檢驗(yàn)結(jié)果的比較來選擇克里金預(yù)測更為精確的變異函數(shù)。結(jié)果表明，與MOM法相比，在多數(shù)情況下REML法求得的變異函數(shù)進(jìn)行克里金插值更精確。該結(jié)論是否適用于小流域以外的其它研究尺度，還需要進(jìn)一步探討。

（2）土地利用方式對土壤有機(jī)碳和全氮的空間分布有重要影響，灌木林和天然草地土壤有機(jī)碳和全氮水平最高，棄耕地其次，梯田、果園、人工草地最低。

（3）土層深度對土壤全氮空間變異影響較小，對土壤有機(jī)碳影響較大，表層土壤有機(jī)碳含量及變異程度明顯高于下層土壤。

（4）退耕還林等生態(tài)恢復(fù)措施，對于黃土高原小流域增加土壤碳氮儲(chǔ)量和提高土壤質(zhì)量是有效可行的。應(yīng)堅(jiān)持退耕還林等生態(tài)恢復(fù)的實(shí)施，合理優(yōu)化土地利用方式，優(yōu)先選擇種植人工灌木林和恢復(fù)天然草地，其次是耕地撂荒，最后是種植人工草被。

［1］ 王軍，傅伯杰，邱揚(yáng)，等.黃土高原小流域土壤養(yǎng)分的空間異質(zhì)性［J］.生態(tài)學(xué)報(bào)，2002，22（8）：1173－1178.

［2］ Shen Runping，Sun Bo，Zhao Qiguo.Spatial and temporal variability of N，P and K balances for agroecosystems in China［J］.Pedoshere，2005，15（3）：347－355.

［3］ Fu Bojie，Wang Yafeng，Lu Yihe，et al.The effects of land－use combinations on soil erosion：A case study in the Loess Plateau of China［J］.Progress in Physical Geography，2009，33（6）：793－804.

［4］ Zheng Fenli.Effect of vegetation changes on soil erosion on the Loess Plateau［J］.Pedosphere，2006，16（4）：420－427.

［5］ 程勵(lì)勵(lì)，文啟孝，林心.內(nèi)蒙古自治區(qū)土壤中有機(jī)碳、全氮和固定態(tài)銨的貯量［J］.土壤.1994，26（5）：248－252.

［6］ Wei Jianbing，Xiao Duning，Zhang Xingyi，et al.Spatial variability of soil organic carbon in relation to environmental factors of a typical small watershed in the black soil region，Northeast China［J］.Environmental Monitoring and Assessment，2006，121（1／3）：597－613.

［7］ Rossi R E，Mulls D J，Journel A G，et al.Geostatistical tools for modeling and interpreting ecological spatial dependence［J］.Ecological Monographs，1992，62（2）：277－314.

［8］ Van M M，Pannier J，Hofman G，et al.Spatial and temporal changes of soil C after establishment of a pasture on a long－term cultivated vertisol（Martinique）［J］.Geoderma，1996，94（1）：43－58.

［9］ Wang Zongming，Zhang Bai，Song Kaisan，et al.Spatial variability of soil organic carbon under maize monoculture in the Songnen Plain，Northeast China［J］.Pedosphere，2010，20（1）：80－89.

［10］ 李元壽，張人禾，王根緒，等.青藏高原典型高寒草甸區(qū)土壤有機(jī)碳氮的變異特征［J］.環(huán)境科學(xué)，2009，30（6）：248－253.

［11］ Matheron G L.Variables Régionalisées et leur Estimation：Une Application de la Théorie de Fonctions Aléatoires aux Sciences de la Nature［M］.Paris：Masson et Cie，1965.

［12］ Kerry R，Oliver M A.Comparing sampling needs for variograms of soil properties computed by the method of moments and residual maximum likelihood［J］.Geoderma，Pedometrics，2007，140（4）：383－396.

［13］ Webster R，Oliver M A.Sample adequately to estimate variograms of soil properties［J］.Journal of Soil Sci－ence，2007，43（1）：177－192.

［14］ Lucas I F J，F(xiàn)rans J M，Wel V D.Accuracy assessment of satellite derived land cover data：A review［J］.Potogammetric Engineering ＆ Remote Sensing，1994，60（4）：410－432.

［15］ Fang Xuan，Xue Zhijing，Li Bicheng，et al.Soil organic carbon distribution in relation to land use and its storage in a small watershed of the Loess Plateau，China［J］.Catena，2012，88（1）：6－13.

［16］ Patterson H D，Thompson R.Recovery of interblock information when block sizes are unequal［J］.Biometrika，1971，58（2）：545－554.

［17］ Chai Xurong，Shen Chongyang，Yuan Xiaoyong，et al.Spatial prediction of soil organic matter in the presence of different external trends with REML—EBLUP［J］.Geoderma，2008，148（2）：159－166.

［18］ Lark R M，Cullis B R.Model－based analysis using REML for inference from systematically sampled data on soil［J］.European Journal of Soil Science，2004，55（4）：799－813.

［19］ 牛文杰，朱大培，陳其明.滑動(dòng)鄰域克里金插值法的改進(jìn)［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2001，13（8）：752－756.

［20］ Pardo I E.Inference of spatial indicator covariance parameters by maximum likelihood using MLREML［J］.Computers and Geosciences，1998，24（5）：453－464.

［21］ Cambardella C A，Moorman T B，Novak J M，et al.Field－scale variability of soil properties in Central Iowa soils［J］.Soil Sci.Soc.Am.J.，1994，58（5）：1501－1511.

［22］ 陳皓，章申.黃土地區(qū)氮磷流失的模擬研究［J］.地理科學(xué)，1991，11（2）：142－148.

［23］ Groenendijk F M，Condron L M，Rijkse W C.Effect of afforestation on organic carbon，nitrogen，and sulfur concentration in New Zealand hill country soils［J］.Geoderma，2002，108（1）：91－100.