天山云杉林土壤养分空间异质性及环境解释

席丽，李思瑶，夏晓莹，王 杰，马小龙，米尔扎提·柯尼加里木，阿丽耶·麦麦提，王卫霞

（新疆农业大学 a.林学与风景园林学院；b.干旱区林业生态与产业技术重点实验室，新疆 乌鲁木齐 830052）

土壤养分作为土壤最基本的属性之一，其丰缺程度及存在形态直接决定着土壤的肥力状况，影响着植被的生长繁衍[1]，同时土壤养分在气候、植被、母质、地形、人类活动等诸多因素的共同作用下往往呈现出高度的空间异质性[2-3]。有研究表明即使在土壤类型、土壤质地、成土母质相对一致区域内，同一时刻土壤养分在空间上也表现出明显的差异，这种差异在小尺度空间上观察依然存在[4-5]。森林土壤是森林环境下形成的一种土壤资源，其理化特性在小尺度上受林分年龄[6]、林分密度[7]、地形因子[8]（海拔、坡度、坡向等）及土壤相关因子[9]等多因素的影响，森林土壤养分空间分布特征将直接影响林地生产力的高低及生态系统恢复与重建的方向，因此，研究森林土壤养分分布特征具有其自身独特性，是充分认识其发育格局与合理经营管理的基础。

地统计学是分析土壤空间特性和研究其变异规律最有效的方法之一[10-12]，很多学者利用地统计学方法结合GIS 技术针对不同区域的土壤养分空间分布特征做了许多工作。王云全等[13]借助GIS 技术和地统计学方法对云南滇西地区耕地土壤养分空间的变化情况及其影响因素进行了分析和研究；张志坚等[5]基于地统计学和GIS 对江西省森林土壤养分空间分布特征进行了研究；梁贵等[14]以湘中丘陵区石栎-青冈栎常绿阔叶林为研究对象，利用地统计学模型对长沙大山冲国有林场土壤钾含量的空间分布规律和影响因素进行了研究。基于GIS 技术和地统计学等方法研究土壤养分空间异质性已日趋成熟，但在新疆地区天山云杉林土壤养分空间异质性及影响因素的研究鲜有报道。天山云杉Picea schrenkiana是新疆山地森林中分布最广、蓄积量最大的森林生态树种，以天山云杉为主体的天山山地森林生态系统对天山的水源涵养、水土保持和林区生态系统的形成和维护起着重要作用。为此，本研究以天山云杉林为研究对象，结合野外调查和室内分析测定，借助地统计学方法和GIS 技术研究天山云杉林土壤养分的空间变异性及其影响因素，以期为新疆山地森林生态系统的生态恢复及天山云杉林的可持续生态经营提供参考依据。

1 研究区概况

本研究以新疆农业大学实习林场为研究区。新疆农业大学实习林场（43°16′～43°26′N，86°46′～87°00′E）位于新疆天山北麓中部的头屯河上游。林场平均海拔约2 200 m，年平均温度3℃，年降水量500～600 mm，年降水分布不均，春夏偏多，约占全年总降水量的70%左右。地势南高北低，坡度多在10～40°。土壤母质主要为黄土状物质和钙质风化物，其土壤类型为普通灰褐色森林土，土层50～60 cm 处开始出现钙积层，质地较黏重，为中壤至重壤。森林主要分布在阴坡，一般生长较好，基本上呈连续的带状分布，以天山云杉纯林为主。

2 研究方法

2.1 数据的采集

2020 年6—8 月，在研究区内采用机械抽样方法布设样地，样地规格为30 m×30 m，共布设195 个样地，样地分布如图1 所示。用GPS 定位每个样地坐标，并记录每个样地的海拔、坡度、坡向、坡位等地形因子，通过林分调查获取郁闭度、林分年龄等林分因子（表1），并在每个样地采集表层0～20 cm 土壤样品，采集的土样带回实验室自然风干，粉碎过筛后进行有机质（OM）、全氮（TN）、碱解氮（AN）、速效磷（AP）、速效钾（AK）、pH 值等指标的测定，同时采集环刀土样，105℃下烘干至恒质量，测定土壤含水量、容重。土壤样品有机质的测定采用重铬酸钾容量法、全氮采用凯氏定氮法、碱解氮采用碱解扩散法、速效磷采用钼锑抗分光光度法、速效钾采用醋酸铵浸提，火焰光度计法、土壤pH 值采用酸度计测定法、土壤含水量采用标准烘干法、土壤容重采用环刀法测定[15]。

图1 研究区及采样点概况图Fig.1 Survey of the study area and sampling plots

表1 研究区采样点属性分布表Table 1 The distribution of sample properties in study area

2.2 数据处理及分析方法

采用SPSS 26.0 软件对土壤养分指标进行描述性分析；应用Canoco 5.0 软件对土壤养分做环境解释RDA 分析，排序轴特征值采用蒙特卡罗检验显著性，前向选择法检验9个影响因子的边际作用，依影响因子特征值进行排序。利用GS+9.0 软件对土壤各养分指标进行半方差分析与理论模型拟合；应用ArcGIS 10.4 软件中的地统计学分析模块，采用普通Kriging 插值法预测未知点的估计值生成土壤养分空间分布图，其中半方差函数分析是用到的核心方法之一，是研究土壤空间异质性有力工具，是进行空间插值的基础,能把土壤养分异质性和生态学过程有效地联系起来[16]，其公式[17]如下：

式中：γ(h)为样点空间距离为h时的半方差函数；N(h)是分割距离为h时的样本点对总数；Z(xi)为Z(x)在空间位置Xi处的实测值；Z(xi+h)是Z(x)在Xi处距离偏离h的实测值。在通常情况下，半方差函数值随着分割距离h的增加而增大，并在一定的间距（变程）升大到一个基本稳定的常数（C+CO），当h=0 时半方差函数值为块金值（CO）。半方差函数模型有球状（Spherical）、高斯（Gaussian）、指数（Exponential）和线性（Linear）等模型[18]。

3 结果与分析

3.1 土壤养分统计特征

土壤养分的描述性统计结果如表2 所示。研究区内土壤有机质、全氮、碱解氮、速效磷、速效钾的平均值分别为139.80 g/kg、3.85 g/kg、113.79 mg/kg、25.95 mg/kg、179.20 mg/kg，参 照全国第二次土壤普查养分分级标准[19],样地中土壤有机质、全氮、速效磷、速效钾含量丰富，尤其是有机质与全氮。根据划分标准，当变异系数CV ≤0.1 时，为弱变异性；当 0.1 ＜CV ＜1.0 时，为中等变异；当CV ≥1.0 时，为强变异性[20]。本研究区内土壤各养分指标的变异系数范围为45.27%～66.04%，表现为中等程度变异，各项土壤养分指标的变异程度表现为有机质＞全氮＞速效磷＞速效钾＞碱解氮。研究区内土壤各养分的偏态数均大于零，说明均为正偏态分布，各指标经过对数转换后均服从正态分布。

表2 土壤养分描述性统计特征Table 2 Descriptive statistical characteristics of soil nutrients

3.2 土壤养分空间异质性

描述性统计分析可以反映土壤养分的基本信息，但无法定量反映土壤养分的结构性、随机性、独立性、相关性及空间变异性[21-22]，而地统计学可以有效地反映这些特性。利用地统计学软件GS+9.0 对研究区各土壤养分指标进行半方差函数拟合，以残差平方和（RSS）最小，决定系数（R2）最大为原则[11,23]，得到各指标的最佳拟合模型及相关参数见表3 和图2。拟合结果显示，有机质、速效磷、速效钾、全氮的最优拟合模型为指数模型，碱解氮的最优拟合模型为高斯模型。各养分指标R2在0.494～0.830，RSS 均接近0，说明各养分变量最佳模型的拟合效果较好，能很好地反映各土壤养分的空间结构特征。

图2 土壤养分的半方差分析Fig.2 Semi-variance analysis of soil nutrients

表3 土壤养分半方差函数和相应参数Table 3 Semi-variance function and corresponding parameters of soil nutrients

块金值（CO）代表随机变异的量，而基台值（C+CO）则代表某变量由随机因素和结构因素共同引起的变异总和。块金值与基台值的比称为空间结构比，CO/(C+CO)可以反映土壤养分变量的空间相关程度。按照Cambardella等[24]的划分标准，若比值＜25%，说明系统具有强烈空间相关性；若比值在25%～75%，说明系统具有中等空间自相关性；比值＞75%则说明系统空间自相关性很弱。由表3 可知，研究区内土壤有机质、全氮、碱解氮、速效钾具有强烈的空间自相关性，其空间变异主要来源于结构因素。速效磷属于中等程度空间自相关性，其变异受随机因素与结构因素共同作用影响。变程反映土壤属性有效自相关距离，研究区土壤养分除速效磷变程为1 722 m，有机质、碱解氮、速效钾、全氮变程波动范围较小，均在97～159 m。

3.3 土壤养分空间分布

为了进一步了解天山云杉林土壤养分空间分布特征，结合半方差函数最优拟合模型及相关参数，采用普通Kriging 插值法[25]进行空间插值，并用交叉验证得知各养分插值过程平均误差（ME）、均方根误差（RMSE）参数的绝对值均接近0（表3），表明普通克里格方法能够提供较为精确的土壤养分空间插值，最后绘制空间分布图（图3），并参照全国第二次土壤普查养分分级标准，将插值结果重新分类，并统计各级面积所占的比例，结果见表4。由图3 可知，研究区内土壤有机质、碱解氮和全氮的变化趋势比较相似，高值区与低值区均出现在相同的位置，且变化趋势相同，均从东北方向的中高值区向西南方向的中低值区过渡。速效磷空间上连续性较好，自北到南呈先增加后降低再增加的趋势，自西到东呈减少的趋势。速效钾在空间上的分布格局较复杂，没有明显的变化规律。总体来看，研究区北边的土壤养分比南边要丰富。由表4 可知，有机质与全氮的含量处于1 级的面积比例分别为100%、92.6%，可见这两指标的含量非常丰富。碱解氮处于3 级的面积比例较大。速效磷的分布主要以2 级、3 级为主，速效钾的分布主要有1、2、3 级。

图3 土壤养分空间分布Fig.3 Spatial distribution of soil nutrients

表4 土壤养分插值后的分级面积比例统计Table 4 Statistics of grading area proportions of the soil nutrients after interpolation

3.4 土壤养分空间变异的影响因素

首先进行DCA 排序，发现第一排序轴的长度LGA ＜3，应选择线性模型，即冗余分析RDA 方法，研究土壤养分空间分布与环境因子的关系。分析结果包括土壤养分与环境因子RDA 相关系数及环境因子边际作用检验结果（表5）、RDA 二维排序图（图4）。根据表5 中环境因子与排序轴的相关性大小可知，土壤含水量与第一排序轴呈正相关关系，相关系数为0.654 9，而土壤pH 值、土壤容重、海拔与第一排序轴呈负相关关系，相关性依次为土壤pH 值＞土壤容重＞海拔，第一排序轴主要由土壤含水量、土壤pH 值、林分年龄等环境因子控制。第二排序轴主要由土壤pH 值控制，相关系数为0.257 4。由边际作用检验结果可知，各环境因子对土壤养分空间分布的解释程度大小依次为：土壤含水量＞土壤pH 值＞林分年龄＞土壤容重＞郁闭度＞海拔＞坡向＞坡度＞坡位。前两排序轴的方差解释累计百分比达98.2%，可以很好地解释所选环境因子与土壤养分的关系。

表5 土壤养分与环境因子RDA 相关系数及环境因子边际作用检验结果Table 5 Correlation coefficient between soil nutrients and environmental factors RDA and the marginal effect test results of environmental factors

图4 土壤养分与环境因子RDA 排序Fig.4 RDA ordination of soil nutrients and environmental factors

土壤养分与环境因子RDA 排序图中只包含方差解释百分比较大的2 条排序轴，各环境因子用红色带箭头的连线表示，红色矢量箭头长度与方向代表土壤与此环境因子的相关关系，土壤含水量的矢量箭头在有机质、碱解氮、全氮上的余弦投影较长，表明土壤含水量对三者积累具有重要意义；相反，土壤容重与有机质、碱解氮、全氮的夹角接近180°，说明对其有负相关影响。速效磷与pH 值、海拔夹角接近180°，说明pH 值、海拔对速效磷空间变异的影响较强，即随着土壤pH值、海拔的上升，土壤速效磷含量逐渐减小。速效钾的空间变异主要受郁闭度影响。本研究环境因子总计解释了土壤养分 50.7%的变异信息，其中第一排序轴、第二排序轴分别解释了土壤养分47.9%、1.8%的变异信息。

4 讨论

土壤养分特性的空间变异无论在区域尺度还是中小尺度上都普遍存在[26]。描述性统计分析结果表明研究区内天山云杉林土壤有机质、全氮、碱解氮、速效磷和速效钾均呈中等程度变异，这与王千军等[27]得出的结果较为一致。半方差拟合模型能够较好地反映土壤养分的空间属性特征，研究区内土壤有机质、全氮、碱解氮和速效钾均表现出强烈的空间自相关特性，说明随机干扰因素对以上这些土壤养分指标的影响相对较小，其空间变异主要受结构性因素的影响。其CO/(C+CO)值均小于天山南坡农区[28]，说明本研究区土壤养分的空间自相关性更强，这是因为天山云杉林受人为干扰较少，而天山南坡农区受人为施肥等随机因素影响较大，随机性因素促使其空间相关性减弱[10]。岩石风化是土壤中磷的主要来源，导致磷积累缓慢，分布相对均匀[29]，本研究结果显示速效磷其空间连续性好，表现为中等程度空间自相关性，即成土母质、气候条件较为一致情况下，其空间分布既受林分结构、地形、土壤微生境等自然因素的影响，又受取样设计、试验误差等随机因素的影响。本研究与祁连山青海云杉林[30]的研究相比既有相同点又有不同点，相同点在于土壤有机质、全氮、碱解氮、速效钾均表现为强烈空间自相关性，同时这些养分指标的空间自相关距离均与本研究结果（97～159 m）相近。不同点在于本研究中土壤速效磷的空间自相关距离最大，且表现为中等空间自相关性。而祁连山青海云杉林中土壤速效磷较其他因子的空间自相关距离相对较小，表现为强烈空间自相关性，这可能与天山北坡天山云杉林森林生态系统具有独特复杂的微生境有关，生态系统的结构与功能对土壤特性空间分布格局产生重要影响。天山云杉林5种土壤养分（除速效磷外）空间自相关范围差异不大，生态过程在97～159 m 尺度上起作用，均大于最小采样距离（50 m），变程对土壤养分的取样设计有效性及区域变量的空间内插和制图具有指导作用，若取样为土壤养分空间插值制图服务，取样网格应小于有效变程，有机质、全氮应取150 m，碱解氮、速效钾分别取90 m、130 m，而速效磷应取1 700 m 左右为宜，一般来说，变程以外的取样距离对任何内插和制图均无效[31-32]。因此，同时采集多个土壤养分指标时，为满足土壤取样设计的要求，须采用尺度嵌套的网格模式[33]。天山云杉林土壤养分存在不同程度的空间异质性与有效变程，且这一现象与马桦薇等[34]研究结果一致，说明不同土壤养分自身形成的主导因素不同。土壤养分的空间分布图可以准确地反映整个研究区域土壤养分的空间分布状况。本研究中，林下土壤有机质、全氮、碱解氮的空间分布格局相似性较高，从东北方向的中高值区向西南方向的中低值区过渡，说明三者形成的主导因素可能是一致的；速效磷空间上连续性较好，自北到南呈先增加再降低再增加的趋势，自西到东呈减少的趋势；速效钾的空间分布没有明显的规律，研究结果表明不同养分的空间分布格局不同。

成土母质相对一致的条件下生物富集过程（植物新陈代谢和凋落物的分解与合成）对土壤养分空间分布的影响最大[9]。其中，土壤含水量反映土壤水分与群落内部的湿润情况，能直接影响凋落物与土壤表层的物质能量交换及土壤盐基养分的淋溶程度[35]。本研究结果表明，土壤含水量是影响天山北坡天山云杉林土壤养分的主要驱动因子，同时RDA 排序图也显示，土壤含水量的矢量箭头在有机质、碱解氮、全氮上的余弦投影较长，表明土壤含水量对这3 种养分指标的积累具有重要意义；相反，土壤容重与有机质、碱解氮、全氮的夹角接近180°，说明土壤容重与土壤有机质、碱解氮、全氮呈负相关关系，这一研究结果与前人研究结果一致[34-36]。土壤容重作为表示土壤紧实度的指标，会对植物根系发育及凋落物分解产生影响，进而影响养分元素在土壤中的转换和积累。土壤容重越大，土壤越紧实，这不仅不利于植物根系的生长，同时也会降低凋落物的分解，最终影响植物与土壤间养分元素流通，使得土壤养分含量较低。

本研究借助地统计学和GIS 技术探究土壤养分的空间变异特征，并通过RDA 分析进行解释，研究结果可为天山云杉天然林可持续经营提供理论依据。但本研究仅结合普通克里格方法进行空间插值，没有结合其他空间插值方法研究天山云杉土壤养分空间特征，今后可以结合其他空间插值方法，更加科学准确了解天山云杉林土壤属性空间特征。另外，本研究所选环境因子对土壤养分空间变异的解释程度为 50.7%，因为天然林生境复杂，影响土壤养分的因子较多，后续应综合多方面因素加以考虑，进一步深入研究。

5 结论

1）研究区内土壤有机质、全氮、碱解氮、速效磷、速效钾的平均值分别为139.80 g/kg、3.85 g/kg、113.79 mg/kg、25.95 mg/kg、179.20 mg/kg。各养分指标的变异系数范围为45.27%～66.04%，属于中等程度变异。

2）林下土壤有机质、全氮、速效磷、速效钾的最适模型为指数模型，碱解氮的最适模型为高斯模型。有机质、全氮、碱解氮、速效钾的块金系数分别为11%、13%、14%和10%，表现出强烈的空间自相关性，表明其主要受结构因素影响；速效磷的块金系数为41%，属于中等程度空间自相关性，其变异受随机因素与结构因素共同作用影响。

3）普通克里格插值表明，天山云杉林土壤养分丰富，尤其是碳、氮含量。

4）林下土壤养分与环境因子RDA 排序表明，土壤含水量、pH 值、林分年龄影响着研究区土壤的养分含量水平，所选环境因子总计解释了土壤养分50.7%的变异信息。