基于最大似然法与矩法的黄土高原小流域土壤碳氮空间变异分析

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

土壤由于受到气候、地形、植被、水文、土地利用方式等因子的综合作用，在空间分布上具有高度异质性［1］。土壤特性的空间变异已成为当前土壤科学研究的热点之一［2］。黄土高原是我国乃至全球土壤侵蚀最为严重的地区之一，长期以来不合理的土地利用，加剧了土壤侵蚀和生态环境恶化［3－4］。退耕还林生态恢复实施以来，该地区土地利用与土地覆盖发生重大变化，土壤特性受到影响。土壤碳氮是土壤特性的重要指标，土壤碳氮循环对全球变化有重要作用［5］。开展退耕后黄土高原土壤碳氮空间分异研究，对揭示土壤碳氮分布与土地利用方式等生态因子之间的关系具有重要意义。地统计方法与经典统计学相比，能够更好地将空间格局与生态过程联系起来［6］，有效地揭示空间格局对生态过程的影响［7］，已被广泛应用于土 壤 碳 氮 空 间 分 布 的 研 究［1，6，8－10］。 估算变异函数是采用地统计方法进行空间变异分析的关键。国内在研究土壤碳氮空间变异中基本上都采用 Matheron［11］提 出 的 矩 法 （method of moments，MOM）求变异函数，但该方法对研究对象的样本数量要求较高。Kerry［12］指出，MOM法要求样本数量在采样间距合理的前提下大于100，才能得到充分描述空间异质性的变异函数。Webster和Oliver［13］研究指出，当样本数量无法达到100时，受限最大似然法（residual maximum likelihood，REML）是可选的变异函数估算方法。Kerry［12］采用矩法和受限最大似然法对不同样本数量下精准农业的土壤性质变异函数进行估算，并将两种方法结果进行比较研究，结果表明当样本数量小于100时，采用受限最大似然法估计的变异函数进行空间插值预测更加准确。可见，REML法是土壤空间变异研究可以采用的方法，而该方法在我国土壤空间变异的研究中还鲜有应用。本研究以退耕后的黄土高原寺底沟小流域为对象，采用受限最大似然法和矩法求变异函数，通过对比选择最优的变异函数模型进行地统计插值，以此分析小流域尺度下土壤有机碳氮的空间变异特征，旨在为进一步研究土壤碳氮变异的生态过程以及改善土壤碳氮水平提供科学依据。

1 研究区概况

研究区寺底沟小流域，位于宁夏自治区南部固原市原州区河川乡内，地理范围为东经106°25′13″—106°26′17″，北纬 35°57′34″—35°59′17″，总面积为5.16km2。该小流域地处黄土高原丘陵沟壑区，海拔1 590～1 815m之间，沟坡地占的88%，坡度小于8°的平缓地仅占12%；属温带半干旱大陆性气候，年均气温与降水量分别是6.9°C和419mm（1982—2002年）；土壤类型主要为黄土母质上发育的黄绵土；本土植被是天然草地，由于过度开垦而一度减少，直到2002年退耕还林后才有所恢复。退耕前植被类型以耕地和天然草地占主导地位，经过退耕和农业结构的调整，植被覆盖发生重大变化，植被类型包括天然草地、人工灌木林〔柠条（C.korshinskii）和沙棘（H.rhamnoides）〕和人工草地〔苜蓿（M.sativa）〕、耕地、果园和弃耕地等。其中，苜蓿与柠条为退耕植被。

2 材料与方法

2.1 土样采集与分析

小流域尺度下土壤类型和气候水文条件可视为基本上同质，考虑土地利用方式和地形因素，按照不同土地利用方式和不同坡位（坡上、坡中、坡下），于2010年7月对寺底沟小流域进行采样，共选取46个代表性的土壤样点（图1）。各样点分0—10cm，10—30cm，30—60cm共3层采集土壤样品，每个样点设5次重复，5次样本充分混合后取1kg带回实验室。土壤样品采回后经自然风干，去除根系石块，研磨过筛后分析土壤有机质和全氮，有机质采用重铬酸钾外加热法测定，全氮采用凯氏定氮法测定。

2.2 数据准备

样点选择所参考的基本图件为2010年寺底沟的土地利用图和1982年1∶1万地形图，其中土地利用图（图1）是由2008年SPOT5高分辨率卫星影像人机交互解译后，再经2010年7月实地调查对每个地块逐个更新获得。采用Kappa指数对该影像分类精度进行定量评估，Kappa指数为0.89，超过了最低允许判别精度 Kappa指数≥0.70的要求［14－15］。土地利用图和地形图投影坐标均统一为阿尔伯斯（Albers）投影和Krasovsky－1940坐标系统。采样点坐标由GPS确定导入ArcGIS软件，原经纬度坐标经投影转换为与上述图件一致的投影坐标系统，形成样点分布图。将46个采样点的各土层有机碳、全氮数据添加到采样点属性表中，以分析土壤碳氮的空间变异特征。

2.3 地统计方法

本研究采用地统计方法分析土壤碳氮的空间变异特征，变异函数用来分析土壤碳氮的空间格局，普通克里金插值方法用来进行土壤碳氮空间分布的预测模拟。估计变异函数参数最常用的方法是矩法，即常用公式（1）计算得到实验半方差，然后采用加权最小二乘法对这些实验半方差进行拟合得到理论变异函数的各参数值［12］。变异函数的理论模型一般包括球面模型、指数模型、高斯模型和线性模型，最优理论模型由反映变异函数拟合程度的决定系数大小来选择。

式中：N（h）——距离等于h时的点对数是；Z（xi）——在位置xi处的数值；Z（xi＋h）——在xi＋h处的数值。

图1 寺底沟土地利用与样点分布

由于本研究样本量小于100，所以同时采用受限最大似然法来估计土壤变量的变异函数。受限最大似然法是Patterson和Thompson提出的［16］，该方法估计变异函数参数是无偏估计［17］，比最小二乘法的结果更加可信［18］。与矩法估计变异函数相比，可直接利用原始数据进行估计，不需要搜索不同空间距离的数据对，避免因局部空间内数据点对数太少而使所求变异函数不可靠［19］。受限最大似然法估计变异函数的原理和公式在文献中［12，19－20］均有叙述。本研究采用Pardo—Iguzquiza［20］提出的 MLREML程序来实现受限最大似然法（REML）方法的变异函数计算。MLREML程序提供了球面模型、指数模型、高斯模型3种理论模型，并且选择其中负对数似然函数最小的作为变异函数最优理论模型。

采用矩法和受限最大似然法两种方法估计变异函数参数，然后通过交叉检验来选择其中最优的变异函数［12］。交叉检验中采用平均误差（ME）和标准化克里金方差（MSDR）来衡量两种方法估计的变异函数进行克里金预测的效果。平均误差是指每个点的实际观察值与普通克里金预测值之间差值的平均值，标准化克里金方差是指每个点的误差平方与克里金方差比值的平均值，分别采用公式（2）—（3）求得：

式中：Zxi——克里金预测值：z（xi）——实际观察值。ME值越接近于0，MSDR值越接近于1，表明变异函数模型越准确。本研究采用GS＋和Excel软件进行模拟变异函数模型，利用ArcGIS 9.2进行普通克里金插值。

3 结果分析与讨论

3.1 土壤有机碳、全氮和碳氮比的统计分析

研究区土壤各指标的描述性统计与正态分布检验结果详见表1。土壤表层（0—10cm）平均有机碳含量和平均全氮含量分别为9.16和0.8g／kg，且都随着土层的加深而减少。各指标的变异系数均在0.1～0.9之间，表明有机碳含量和全氮含量都属于中等变异。从不同土层土壤性质的变异程度来看，有机碳含量的变异程度随着土层加深略有下降，全氮含量的变异程度略有增加。单样本柯尔莫哥洛夫—斯米尔诺夫检验（K—Stest）结果表明，各土层土壤有机碳和全氮均呈正态分布（K—S检验值均大于0.05），其原始数据均可直接用于地统计分析。

表1 研究区土壤有机碳含量和全氮的描述性统计与正态分布检验

3.2 土壤有机碳和全氮的空间异质性

由REML和MOM两种方法估计土壤碳氮变异函数的交叉检验对比结果详见表2。除了10—30cm土壤全氮，其它变量均是采用受限最大似然法求得变异函数的标准克里金方差（MSDR）更接近于1，平均误差（MS）更接近于0。此结果表明大多数情况下，采用受限最大似然法估计的变异函数进行克里金预测，效果优于钜法估计的变异函数。因此，除了10—30 cm土层的土壤全氮克里金插值采用MOM法变异函数外，其它土壤有机碳和全氮变量均采用REML法估计的变异函数参数进行克里金插值。

由土壤有机碳和全氮变异函数模型及相关参数可以看出（表2），块金方差表示实验误差和小于最小取样尺度引起的随机变异，结构方差代表由空间自相关引起的结构性变异，块金效应指块金方差占空间总变异的比例，可以反映空间变异程度的大小［21］。块金效应小于0.25表示变量具有强烈的空间相关性，在0.25～0.75之间表示具有中等空间相关性，大于0.75表示空间相关弱。由表2可以看出，土壤有机碳含量和全氮含量的空间相关性基本上都属于中等。与全氮含量相比，土壤有机碳含量的随机变异比例较大，结构性变异比例较小，空间相关性相对较弱。不同土层深度相比较，无论有机碳还是全氮，都表现为表层土壤比下层土壤随机变异比例较大，结构性变异比例较小，空间相关性较弱。

表2 研究区土壤有机碳（SOC）和全氮（TN）变异函数模型的相关参数及交叉检验结果

变程大小代表空间相关性范围的大小，即变程之内变量具有空间相关性，反之不存在空间相关性。表2显示，在0—10cm，10—30cm，30—60cm这3个不同土层深度下，土壤有机碳变程分别为952.33，605.24和576.82m，而全氮变程分别为799.93，291和486.01m，说明土壤有机碳空间相关性范围大于全氮。土壤全氮空间相关尺度较小与氮肥使用和土壤氮素流失有关。一方面，研究区地形破碎、土壤疏松，土壤侵蚀现象多发，土壤养分流失严重，而氮素在土坡中具有易移动特性，而使得土壤流失及淋溶造成的土壤氮素损失较有机质多，从而有机碳空间相关性范围比全氮大［1，22］。另一方面，氮肥是研究区农户种植作物和果树使用较多的肥料，氮肥施用差异及产生这种差异的随机性可能造成土壤全氮空间相关尺度的减少。比较不同土层下土壤变量的变程大小可以看出，表层土壤（0—10cm）有机碳和全氮空间相关性范围比下层土壤大，其中土壤有机碳变异范围随着土层加深而减小。

3.3 土壤有机碳和全氮空间分布特征

土壤有机碳和全氮的地统计插值结果如图2所示。结合图1可以看出，不同土地利用方式下的土壤有机碳和全氮含量总体表现为灌木林最大、天然草地和弃耕地居中、梯田和人工草地最小。其中，有机碳含量在0—10cm土层的最大值出现在沙棘林地，20—30cm和30—60cm土层的有机碳含量最大值出现在柠条林地；土壤全氮含量在0—10cm土层的最大值出现在沙棘林地，10—30cm土层的最大值出现在弃耕地，30—60cm土层的最大值出现在柠条林地和天然草地。分析结果表明，灌木林、天然草地和弃耕地与其它土地利用方式相比，表现出更好的土壤有机碳氮蓄存能力，证实了退耕还林政策从总体上对增加土壤碳氮储量是有积极作用的。一方面，土地利用方式影响土壤有机碳氮水平，耕地转化为多年生植物，能够增加土壤有机碳含量；相反，当耕地转化为多年生植物时，土壤有机碳会以CO2的形式流失进入大气圈［23］。另一方面，耕地施肥没能提高土壤全氮水平，这与该小流域水土流失严重密切相关。另外，人工草地作为退耕还林中重要的退耕植被，对土壤碳氮储量没有明显贡献，这是由于人工草地（苜蓿）作为当地农民畜牧业和家庭养殖业的天然饲料而被频繁收割，导致有机碳和氮无法在土壤中积累。尽管如此，考虑到人工草地防治水土流失的积极作用以及其作为农户畜牧养殖的主要饲料来源，建议人工草地在控制收割次数的情况下，继续保留作为退耕还林的重要植被类型。

不同土壤深度对土壤有机碳氮空间分布的影响表现为前者随着土壤深度变化不明显，而后者明显表现为表层（0—10cm）土壤有机碳含量及变异程度大于下层土壤。表明表层土壤对于有机碳积累很重要，且土地利用方式对土壤有机碳汇的作用主要发生在土壤表层，可能与表层土壤受到植被覆盖、枯落物和根系的植物量以及人类干扰的影响较大有关。

图2 寺底沟土壤有机碳和全氮空间分布

4 结论

（1）采用REML法与MOM法分别对不同深度下土壤有机碳和全氮的变异函数进行估算，并通过交叉检验结果的比较来选择克里金预测更为精确的变异函数。结果表明，与MOM法相比，在多数情况下REML法求得的变异函数进行克里金插值更精确。该结论是否适用于小流域以外的其它研究尺度，还需要进一步探讨。

（2）土地利用方式对土壤有机碳和全氮的空间分布有重要影响，灌木林和天然草地土壤有机碳和全氮水平最高，弃耕地其次，梯田、果园、人工草地最低。

（3）土层深度对土壤全氮空间变异影响较小，对土壤有机碳影响较大，表层土壤有机碳含量及变异程度明显高于下层土壤。

（4）退耕还林等生态恢复措施，对于黄土高原小流域增加土壤碳氮储量和提高土壤质量是有效可行的。应坚持退耕还林等生态恢复的实施，合理优化土地利用方式，优先选择种植人工灌木林和恢复天然草地，其次是耕地撂荒，最后是种植人工草被。

［1］ 王军，傅伯杰，邱扬，等.黄土高原小流域土壤养分的空间异质性［J］.生态学报，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］ 程励励，文启孝，林心.内蒙古自治区土壤中有机碳、全氮和固定态铵的贮量［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］ 李元寿，张人禾，王根绪，等.青藏高原典型高寒草甸区土壤有机碳氮的变异特征［J］.环境科学，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，Frans 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］ 牛文杰，朱大培，陈其明.滑动邻域克里金插值法的改进［J］.计算机辅助设计与图形学学报，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］ 陈皓，章申.黄土地区氮磷流失的模拟研究［J］.地理科学，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.