空间异质性是指系统的某种属性在空间分布上的不均匀性和复杂程度^[1-3]。在生态系统中，土壤pH值、有机质以及氮、磷、钾等属性因子都普遍存在空间异质性^[4-5]，对植物个体或种群的分布以及植被的空间格局特征具有重要影响^[6]。研究土壤性质的空间异质性及其驱动因素对揭示生态系统功能和生物多样性具有重要意义。

在研究土壤的空间异质性时，很多研究表明施肥、土地利用、耕作等人类活动导致了土壤性质各指标的变化，属随机因素引起的变异^[7-9]；然而，土壤属性指标在很大程度上表现为变异的空间自相关性^[10-11]，这说明结构变异是土壤空间异质性的重要来源。在自然因素引起的结构性变异中，地形是一个重要因素。有研究^{[5, 12-13]}表明，在中小尺度上，土壤性质的空间变异较为明显地受海拔、凹凸度、坡度等地形因子的影响。张娜等^[10]发现，酸性降雨对山脊和坡面的直接冲刷明显大于容易覆盖较厚凋落物的沟谷，导致冲刷后山脊和坡面土壤钙、镁等碱性元素的淋溶增加以及对应氢离子浓度的增加，使得山脊和坡面的pH值相对大于沟谷，说明土壤pH值的空间变异主要受地形因子凹凸度的影响较大。彭晚霞等^[11]认为，海拔和坡度对土壤有机质、全氮和有效磷等有极显著的影响，海拔和坡度等地形因子主要通过控制水土流失对土壤的空间变异产生影响。这些研究可在一定程度上揭示土壤和地形因子间的相关关系，有助于我们进一步认识土壤的空间变异机制；但这些研究多局限于地形对土壤空间变异影响的定性分析，或多是针对农耕地、水土流失防治区等受干扰地域上的研究。

云杉天然次生林是我国华北地区亚高山典型的寒温性植被类型之一，对华北地区的生态平衡和水土保持具有重要意义。该类森林资源现存状况好、人为干扰较少，对其林地土壤性质的研究也一直受到重视。对关帝山4 hm²云杉次生林固定监测样地的研究已表明，云杉次生林土壤性质存在空间上的分化，这种分化可分别解释31.12%、44.22%和32.04%的乔木、小乔木和灌木的分布^[14]。以往研究表明，地形因子对不同土壤性质指标的空间变异有不同程度的影响^[15]。为探究云杉次生林土壤存在的变异特征，地形作为结构性因素对土壤性质的空间变异的驱动作用，以及其对植物的分布格局的影响，笔者在上述4 hm²云杉次生林样地研究的基础上，参考CTFS(Center for Tropical Forest Science)技术规范野外取样并进行土壤指标测定，分析土壤pH值、有效氮、有效磷、有效钾和有机质的空间变异和分布格局，并基于RDA排序和方差分解定量评价各地形因子对土壤各指标空间变异的影响程度。本研究旨在揭示森林土壤空间变异特征及其与地形因子的相关性，为深入认识华北地区植物分布和物种共存提供基础依据。

关帝山云杉次生林固定监测样地(简称为GDS样地)，位于山西省西部、吕梁山脉中段关帝山庞泉沟国家自然保护区(E 111°22′~111°33′，N 37°45′~37°55′)，样地东西和南北长分别为200 m，平均海拔2 007 m，最高海拔2 039 m，最低海拔1 976 m，最大高差63 m (图 1)；地形呈北低南高，以东南坡向为主，样地中部南北方向有一沟谷，西北部区域靠近样地边界处地形较为复杂。

图 1(Figure 1)

图 1 GDS样地等高线图 Figure 1 Contourin GDS plot

样地所处区域属暖温带大陆性季风气候，与省内同纬度地区相比，气温偏低，变幅较大，空气湿度偏高，为典型山地气候。年平均气温4.3 ℃，1月平均温度-10.2 ℃，7月平均温度17.5 ℃；年平均降水量822.6 mm，年平均蒸发量1 268 mm。土壤为山地棕壤，平均土层厚度70~80 cm, 其中腐殖质层厚度10 cm^[16]。样地内植物较为丰富，胸径(DBH)≥1 cm的木本植物11科22属30种，10 266株，以寒温性常绿针叶树种青杄(Picea wilsonii)为主，伴生有华北落叶松(Larix principis-rupprechtii)、红桦(Betula albo-sinensis)、白桦(Betula platyphylla)、辽东栎(Quercus wutaishansea)等，主林层林龄约70 a，林分郁闭度为0.70，平均高23.6 m，群落结构稳定，林下植物更新良好。

2016年7月，参照CTFS的土壤采样方案^[6]，结合GDS样地实际情况，将样地分割成49个30 m×30 m的网格，以每个网格结点为基点沿一定方向(可从东、西、南、北、东南、东北、西南、西北随机选取1个方向)向外延伸取土样(处于样地边界上的基点只从5个方向选取，以保证所有采样点都落在样地内)，从距离基点2、5和15 m处随机选择2处作为延伸采样点，由此共计采样点190个，取样深度0~10 cm，采样点空间分布图如图 2所示。为探讨土壤理化性质在垂直方向上的分布特征，将样地分割成9个90 m×90 m的网格，只在每个网格结点处分0~10，10~30和30~50 cm层分别取样。

图 2(Figure 2)

图 2 GDS样地土壤采样点图 Figure 2 Soil sampling points in GDS plot

土样带回实验室后，过2 mm筛后备用。将装有鲜土的铝盒放在105 ℃的烘箱中干燥12 h，测定土壤含水量。本研究测定的土壤性质指标包括pH值、有效氮、有效磷、有效钾和有机质。土壤pH值采用玻璃电极法(土水比1:2.5) 测定，有效氮采用碱解扩散吸收法测定，有效磷采用0.05 mol·L^-1 HCl-0.025H₂SO₄浸提—钼锑抗比色法测定，有效钾采用醋酸铵提取—火焰光度计法测定，有机质采用重铬酸钾外加热氧化法测定^[17]。

选取海拔、坡度和凹凸度3个地形因子作为分析参数。以样地第1次本底调查时测量的每个20 m×20 m样方4个顶点的海拔数据为基础，分别计算各样方平均海拔^[18-19]、坡度^[20]和凹凸度^[21]。

利用经典统计学方法分析土壤pH值、有效氮、有效磷、有效钾和有机质的极值、平均值、标准差和变异系数，以评价各指标的水平和垂直空间分布特征。

利用地统计学的半方差函数对各项指标进行变异计算，拟合最优变异函数理论模型，并根据其主要参数，定量分析各土壤指标的空间异质性特征。其中，模型参数变程(A)表示各变量空间自相关变异的尺度，块金值(C₀)与基台值(C₀+C)可用来表示结构方差比[1-C₀/(C₀+C)]，是衡量空间自相关变异程度的重要指标，当结构方差比>75%时，变异的空间自相关性很强；结构方差比 < 25%时，变异的空间自相关性很弱；介于25%~75%的结构方差比表示变异具有中等程度的空间自相关性^[22]。

为获取土壤性质各指标的空间分布格局，对样点各指标进行Kriging插值，得到样地土壤pH值、有效氮、有效磷、有效钾和有机质的空间分布图。

为探讨地形因子对土壤性质的影响，通过RDA排序的方法直观描述各地形因子对土壤性质解释量的相对大小以及二者相关性。排序图中箭头长度表征各地形因子对土壤性质解释量的相对大小；各地形因子和土壤各指标箭头间的夹角表征二者的相关性，其中夹角角度为0°~90°时表示2个变量间呈正相关，夹角角度为90°~180°时，表示负相关，夹角角度为90°时表示无显著相关^[23]。

RDA排序时，地形因子矩阵由每个20 m×20 m样方的平均海拔、凹凸度和坡度组成；土壤性质矩阵由经Kriging空间插值后每个20 m×20 m样方的平均pH值、有效氮、有效磷、有效钾和有机质组成。进一步用方差分解定量评价地形因子对土壤性质的总解释量和独立解释量。

其中，经典统计分析在R 3.3.1软件中完成。地统计学的半方差函数分析在GS⁺7.0软件中完成。Kriging空间插值在Arc-GIS 9.3软件中完成。RDA排序和方差分解通过R语言Vegan程序包实现^[24]。

由表 1示出，GDS样地土壤pH值、有效氮、有效磷、有效钾和有机质在不同水平空间样点间差异较大。有效氮、有效磷、有效钾和有机质最大值与最小值之间相差倍数分别达到32.80、41.00、8.78和9.21倍，变异系数分别为41.15%、58.36%、46.60%和48.38%。有效磷的平均值为15.88 mg/kg，变异系数为各指标中最高。土壤pH值的平均值为6.06，变异系数达8.78%。K-S检验结果表明，上述土壤性质各指标均符合正态分布。

表 1(Table 1)

表 1 土壤性质各指标水平分布的描述性统计特征 Table 1 Descriptive statistics characteristics of the horizontal distribution of the soil index (n=199) 土壤理化性质 Soil property 平均值 Mean 最小值 Minimum 最大值 Maximum 标准差 Standard derivation 中值 Median 变异系数 Variance coefficient/% P pH 6.06 4.53 7.47 0.53 6.02 8.78 0.044 AN/(mg·kg-1) 30.72 2.14 70.20 12.64 29.47 41.15 0.017 AP/(mg·kg-1) 15.88 1.71 70.11 9.27 13.49 58.36 < 0.001 AK/(mg·kg-1) 159.40 60.30 529.30 74.26 140.70 46.60 < 0.001 OM/(g·kg-1) 146.40 47.80 440.00 70.84 133.20 48.38 < 0.001 注：AN：有效氮；AP：有效磷；AK：有效钾；OM：有机质，下同。Note: AN: Avaivable N; AP: Available P; AK: Available K; OM: Organic matter, the same below.

表 1 土壤性质各指标水平分布的描述性统计特征 Table 1 Descriptive statistics characteristics of the horizontal distribution of the soil index (n=199)

如表 2所示，研究样地土壤中，除pH值外，有效氮、有效磷、有效钾以及有机质含量均以表层土壤0~10 cm层最高，其中，有效氮、有效钾和有机质在表层土壤中的含量显著高于10~30和30~50 cm层土壤(P < 0.05；n=9)。

表 2(Table 2)

表 2 土壤性质各指标的垂直分布特征 Table 2 Characteristics of the vertical distribution of the soil index (n=9) 土层深Soil depth/cm pH AN/(mg·kg-1) AP/(mg·kg-1) AK/(mg·kg-1) OM/(g·kg-1) 0~10 6.50a±0.45 29.18a±20.31 16.64a±11.35 196.18a±156.55 203.52a±108.85 10~30 6.58a±0.62 21.88b±11.83 9.54a±4.92 116.16b±60.53 99.89b±38.66 30~50 6.53a±0.50 18.52b±14.39 14.17a±7.28 125.46b±66.25 112.07b±76.54 注：统计信息为各指标平均值±标准差; 平均同列不同字母表示各土层间差异显著(P < 0.05)。Note: Statistics information means mean value ± standard deviation. Different letters in the same column denoted significant difference between soil layers at 0.05 level.

表 2 土壤性质各指标的垂直分布特征 Table 2 Characteristics of the vertical distribution of the soil index (n=9)

基于表层土壤性质指标显著高于下层土壤(除pH值差异不显著外)的分析结果，利用半方差函数分析0~10 cm层土壤pH值、有效氮、有效磷、有效钾和有机质的空间异质性特征。结果(表 3)表明，本样地的5个土壤性质指标的空间变异特征均最适合用指数模型来描述，模型拟合度均较高。其中，pH值、有效氮、有效磷和有效钾的决定系数为0.539~0.928，有机质的决定系数也达0.368。pH值和有效氮的结构方差比均在25%~75%之间，表现出中等程度的空间自相关性；有效磷、有效钾和有机质的结构方差比>75%，这3个指标的空间分布在研究尺度上存在强烈的空间自相关性。依据半方差函数模型的模拟结果，本研究样地土壤pH值、有效氮、有效磷、有效钾以及有机质分别在58.8、28.0、10.6、8.2和4.9 m的范围内呈现空间自相关性。

表 3(Table 3)

表 3 土壤各指标的变异函数拟合模型及参数 Table 3 Variation function fitting models and parameters for the soil properties 土壤理化性质 Soil property 最适模型 Model type 块金值 Nugget (C0) 基台值 Sill (C+C0) 结构方差比 [1-C0/(C+C0)] 变程 Range (A) R2 pH 指数模型Exponential model 0.155 0.33 0.555 58.8 0.928 AN/(mg·kg-1) 指数模型Exponential model 82.20 164.50 0.500 28.0 0.891 AP/(mg·kg-1) 指数模型Exponential model 1.30 78.48 0.983 10.6 0.559 AK/(mg·kg-1) 指数模型Exponential model 530.00 5 302.00 0.900 8.2 0.539 OM/(g·kg-1) 指数模型Exponential model 580.00 4 978.00 0.883 4.9 0.368

表 3 土壤各指标的变异函数拟合模型及参数 Table 3 Variation function fitting models and parameters for the soil properties

从土壤性质各指标的空间分布格局(图 3)可以看出，样地土壤pH值、有效氮、有效磷、有效钾和有机质均呈现斑块状分布。其中有效氮和有机质的斑块分布在空间上呈现相似的格局，在样地北边界处和西北部土壤有效氮和有机质含量较高，分别达33和166 g/kg以上；但在样地的西北角、西南角以及中部偏东北部区域，这2个指标的含量均较低，分别为27和138 g/kg以下；在样地中部偏东南方向有效氮(27~33 mg/kg)和有机质(138~166 g/kg)均呈现中等含量的斑块状分布。对有效氮和有机质的相关性分析显示(表 4)，二者呈极显著正相关(P < 0.001), 进一步说明二者的空间关联性较好。有效钾和有机质、有效氮和有效钾也存在相似的空间分布格局。相关性分析也进一步表明(表 4)，有效钾和有机质、有效氮和有效钾之间呈极显著的正相关性(P < 0.001)。

图 3(Figure 3)

图中曲线为样地等高线；不同颜色图例代表土壤各指标含量区间。 The curves stand for elevation contours in plot. The legends with different colors stand for the content range of soil properties. 图 3 GDS样地土壤性质各指标空间分布 Figure 3 Spatial distribution of soil properties of GDS plot

表 4(Table 4)

表 4 土壤性质的相关性分析 Table 4 Correlation analysis of soil properties 项目Item pH AN AP AK OM pH 1.00 AN -0.39*** 1.00 AP 0.27** -0.16 1.00 AK -0.06 0.55*** 0.20 1.00 OM -0.09 0.83*** 0.04 0.56*** 1.00 注Note：***：P < 0.001；**：0.001 < P < 0.01；*：0.01 < P < 0.05。

表 4 土壤性质的相关性分析 Table 4 Correlation analysis of soil properties

本研究样地中，土壤pH值与有效氮呈现不同的空间分布格局，如图 3显示，有效氮含量高的斑块(蓝色分布区)土壤pH值(黄色分布区)较低，有效氮含量较低的斑块土壤pH值反而较高，相关性分析显示(表 4)，二者呈极显著负相关性(P < 0.001)。有效磷则与土壤pH值呈相似的空间分布格局(图 3)，相关性分析显示二者具有极显著正相关性(P < 0.01)。有效磷与有效氮、有效钾和有机质并未表现相关性(表 4)。

由RDA排序图(图 4)可得，对土壤各指标空间变异的影响从大到小依次为海拔、凹凸度、坡度。方差分解(表 5)表明，海拔、坡度和凹凸度能分别单独解释土壤性质空间变异的18.96%、6.97%和6.54%，地形因子可总体解释土壤总变异的44.79%。就土壤各指标而言，地形因子可解释pH值、有效氮、有效磷、有效钾和有机质空间变异的68.66%、63.13%、54.46%、29.20%和32.55%。不同地形因子对各土壤指标空间变异的影响也存在差异，其中，海拔对土壤pH值空间变异的影响最大，可独立解释pH值空间变异的40.43%，且为正相关关系；相对于坡度和凹凸度，有效氮和有效磷均受海拔影响较大，不同的是海拔与有效磷呈正相关，能独立解释其空间变异的28.81%，与有效氮则呈负相关，能独立解释其空间变异的19.48%；相对于海拔和坡度，凹凸度对有效钾和有机质变异的影响最大，且与这两个土壤性质指标均呈负相关，可分别独立解释9.37%和6.88%的土壤有效钾和有机质的空间变异。

图 4(Figure 4)

图 4 土壤性质与地形因子的RDA排序图 Figure 4 Ordination diagram of the RDA analysis of the soil properties and terrain factors

表 5(Table 5)

表 5 不同地形因子对土壤性质各指标空间变异影响的方差分解 Table 5 Variance decomposition of the influences of varied terrain factors on spatial heterogeneity of soil property 自变量Variable [a] [b] [c] [d] [e] [f] [g] [h] 总方差Total variation 土壤性质Soil properties 0.189 6 0.069 7 0.065 4 0.000 0 0.015 6 0.107 6 0.000 0 0.552 1 0.447 9 pH 0.404 3 0.194 2 0.038 4 0.000 0 0.000 0 0.000 0 0.049 7 0.313 4 0.686 6 AN 0.194 8 0.053 8 0.049 4 0.000 0 0.028 1 0.305 2 0.000 0 0.368 7 0.631 3 AP 0.288 1 0.072 6 0.076 9 0.066 3 0.040 7 0.000 0 0.000 0 0.455 4 0.544 6 AK 0.022 9 0.000 0 0.093 7 0.000 0 0.013 6 0.161 8 0.000 0 0.708 0 0.292 0 OM 0.037 8 0.030 0 0.068 8 0.000 0 0.027 0 0.161 9 0.000 0 0.674 5 0.325 5 注：[a]单独海拔的解释量；[b]单独坡度的解释量；[c]单独凹凸度的解释量；[d]海拔和坡度2个因子共同的解释量；[e]海拔和凹凸度2个因子共同的解释量；[f]坡度和凹凸度2个因子共同的解释量；[g]海拔、坡度和凹凸度3个因子共同的解释量。Note:[a]=pure elevation；[b]=pure slope；[c]=pure convexity；[d]=shared between elevation and slope but not by convexity；[e]=shared between elevation and convexity but not by slope；[f]=shared between slope and convexity but not by elevation；[g]=shared among elevation, slope and convexity.

表 5 不同地形因子对土壤性质各指标空间变异影响的方差分解 Table 5 Variance decomposition of the influences of varied terrain factors on spatial heterogeneity of soil property

空间异质性是生态系统生物因子和非生物因子共同的重要属性，是导致空间格局的主要原因^[25-26]。森林土壤的空间变异性较其他类型土壤更为复杂且直接影响森林物种的空间分布。笔者研究表明，云杉次生林0~10 cm层土壤的pH值、有效氮、有效磷、有效钾和有机质含量的变异系数分别达8.78%、41.15%、58.36%、46.60%和48.38%，且在小尺度上均表现较强的空间自相关，在整个研究样地上各土壤指标均呈空间斑块状分布格局。这与关帝山华北落叶松林、杨桦阔叶林、云杉、杨桦混交林、辽东栎林、油松林土壤性质空间变异特征的研究^[27-29]相似。说明土壤性质的空间变异是华北亚高山森林土壤共同的显著特征，这可能与华北亚高山山地地形较为复杂有关，也可能与具有不同分布格局的植被对土壤的改善作用有关。

但从斑块的空间分布状况来看，有效氮、有机质和有效钾呈相似的空间斑块状分布格局，而有效磷和pH值则表现出与上述3指标相反的斑块分布格局。这表明土壤性质各指标对地形、植被等因素影响的异质性响应存在差异。

在受人类活动或人为干扰相对较小的森林生态系统中，地形条件作为结构性因素对土壤空间变异的影响不可忽视。相关研究表明，成土母质、土壤类型、植被等自然因素对土壤空间变异也均有影响^{[12-13, 30]}。本次研究表明，地形对土壤性质空间变异的总解释量达44.79%，占所有结构性因素和人为随机因素解释量的近50%。地形因子对土壤性质5个指标分别解释量也均达到29.20%以上，其中针对pH值空间变异的解释量(68.66%)远高于天童山常绿阔叶林^[10]的研究结果(36.52%)。这表明在华北山地寒温性云杉次生林中，地形成为土壤空间变异的重要因素。

方差分解和RDA排序进一步表明，不同地形因子对土壤pH值、有效氮、有效磷、有效钾和有机质各指标空间变异的影响也存在差异。海拔对土壤pH值空间变异的影响最大，独立解释土壤pH值空间变异的40.43%，且二者间正相关。这可能是由于随海拔增加土壤中的Ca、Mg含量降低，土壤中的OH^-增高^[31]。海拔对有效磷和有效氮空间变异的影响也高于其他地形因子，不同的是海拔与有效磷呈正相关，这可能与土壤中磷的可利用情况首先取决于土壤pH值有关^[32]；而海拔与有效氮则呈负相关，这可能是由于海拔越高的地方，气温越低，蒸发量越小，湿度越大，有机质的输入量可能越大，但物质的分解率越小，所以有机质和氮比较富集，但可利用的氮比较缺乏。同时，分解率低也是导致高海拔地方土壤中有效Ca、Mg含量降低的原因^[31]。凹凸度与有机质和有效钾呈负相关，且相对于坡度和海拔其更为显著地影响着有机质和有效钾的空间变异。可能由于凹凸度越大的地方，由于冲刷等原因，越不利于地表物质的积聚，因此降低了有机质和钾的累积。