杨逍虎，赵维俊，3，4，闫春鸣，任小凤，4，武龙庆，薛春，芦倩，张瑞浩

（1.甘肃省祁连山水源涵养林研究院，甘肃 张掖 734000；2.甘肃张掖生态科学研究院，甘肃 张掖 734000；3.甘肃祁连山森林生态系统国家定位观测研究站，甘肃 张掖 734000；4.甘肃省祁连山生态环境中心，甘肃 兰州 730000；5.甘肃农业大学，甘肃 兰州 730000）

森林空间格局是种群更新或者死亡的动态变化与环境异质性相互作用的结果，这种相互作用的结果作为潜在的驱动因子反过来作用种群的动态变化或可能的变化趋势[1]。在研究森林空间分布格局的过程中，常用林分结构进行其空间分布格局及其变异研究[2]。作为森林结构中的重要生态参数：林分个体的密度和树高[3]，其中林木密度通过种内竞争或种间竞争对其生长所需的水分、养分、光照等环境资源产生作用，进而对林木生长、死亡、产量、干形等产生影响[4]。林木树高在生长的过程中亦影响林木对光照、水分、养分等的利用[5]，特别是对更新苗的建立、存活和生长产生不同程度的影响。而且林木密度和树高也是森林经营管理技术措施中最有可能的调控因子[6]。因此，量化林分的密度和树高空间格局及其异质性，能分析其空间分布格局及其与生态学过程之间的内在联系，也能为森林可持续经营提供最直接的理论支持。

为更好地量化植被与环境的复杂关系，需要进行空间分布格局及其动态分析[7]，地统计学在一定程度上规避了传统分析方法的不足，该方法可以研究空间分布数据的随机性、结构性、空间格局与变异，还可以进行空间数据内插、模拟数据的波动性和离散性，能够很好地分析空间分布上存在某种程度的相关性或者连续性的自然界现象。目前已有专家和学者利用地统计学理论对森林的空间异质性进行了研究，如周隽利用半变异函数、空间相关、分形维数和克里格插值法对幼苗更新格局的空间变异进行了分析[8]，姚丹丹等利用半变异函数和克里格插值法对云冷杉针阔混交林叶面积指数的空间异质性进行了分析[9]。这些研究表明利用地统计学能很好地分析植被结构与环境的相互生态作用。

分布在祁连山山区的青海云杉林，其林分结构对祁连山区乃至整个黑河流域的水文过程具有重要的影响[10]。目前对祁连山森林空间异质性研究多是针对成年树或者胸径比较大的个体，或者是针对更新苗进行空间异质性研究，而将试验样地所有林分个体包括更新苗与成树结合来进行分析的研究较少，二者的结合研究有助于跟踪时间序列上的空间变异，进而连贯地认识与空间变异特征相关的种群相互作用及其动态过程[11]。因此，对青海云杉次生林更新苗和成树空间异质性的研究，有助于为祁连山青海云杉林生态系统的近自然恢复提供理论，有助于评估现有的森林管理措施和区域林水综合管理水平，也有助于预测未来森林发展的时空变化。鉴于此，本研究利用国家级祁连山森林生态站已建立的340 m×300 m 青海云杉林动态监测样地，以样地内分布的青海云杉林更新苗和成树为研究对象，选择表征林分结构特征的两个重要生态参数：密度和高度进行调查，将样地分布的更新苗和成树结合起来对其空间异质性进行经典统计分析和地统计学分析，包括半方差函数、空间自相关、分形维数和Kriging 分析，研究其空间异质性程度与变异尺度，揭示更新苗和成树密度和高度的空间变异格局与生态学过程的内在联系，有望为青海云杉林生态系统的可持续管理提供理论依据，也可为退化的青海云杉林生态系统植被生态恢复提供参考。

1 材料与方法

1.1 研究区概况

研究区位于祁连山北麓中段的大野口流域，该流域面积为68.06 km2，海拔变化范围为2 400～4 000 m，年均降水量为400 mm，年均蒸发量为1 488 mm，年均相对湿度为60%，年均气温为5.4℃，最热月的平均气温为19.6℃，最冷月的平均气温为-12.5℃，气候属于温带大陆性高寒半湿润山地气候[12]。研究区水热同期，冬季阴冷潮湿，在其严酷的气候条件下，植被分布具有明显的梯度性分布，不同植被下分布的土壤类型亦具有明显的梯度性分布，从低海拔到高海拔依次分布山地典型灌丛草原带（栗钙土）、寒温性针叶林带（灰褐土）、高山灌丛草甸带（灌丛草甸土）、高山亚冰雪稀疏植被带（寒漠土）。分布在阴坡、半阴坡的青海云杉和灌木树种与分布在阳坡的草地和稀疏零星分布的祁连圆柏林犬齿交错，形成了明显的林草复合流域。受水热条件的影响，植物多样性种类数量较少，流域分布的乔木优势树种主要是青海云杉Picea crassifolia，灌木优势种主要有吉拉柳Salix gilashanica、金露梅Potentilla fruticosa、银露梅Potentilla glabra、箭叶锦鸡儿Caragana jubata、高山绣线菊Spiraea alpina等，草本优势种主要黑穗薹草Carex atrata、委陵菜Potentilla chinensis、珠芽蓼Polygonum viviparum、马蔺Iris lactea、狼毒Stellera chamaejasme等，还有特别是乔木林下分布的苔藓植物，优势种主要有匐灯藓Plagiomnium cuspidatum、山羽藓Abietinella abietina、欧灰藓Hypnum cupressiforme等。人为活动主要是附近的牧民轮流放牧、采集野蘑菇等。

1.2 试验样地的选取及设置

利用国家级祁连山森林生态站在2010—2011年建立的340 m×300 m 祁连山青海云杉林动态监测样地，该样地是依据中国生物多样性网络关于森林大样地的技术规范和参照国际上CTFS 森林大样地建设标准建立的，该样地仅分布青海云杉树种，地形相对平坦，样地的具体信息详见文献[12]。本研究以“空间代替时间”的方法，并参考相关研究，把树高小于1.5 m 和胸径小于5.0 cm 的幼苗、幼树活体都归为更新苗，把树高大于等于1.5 m 和胸径大于等于5.0 cm 的活体都归为成树[13-15]。根据地统计学理论，利用全站仪将该样地分割为255个20 m×20 m 的样方，然后再将每个样方分割成16 个5 m×5 m 的栅格小样方，考虑到样点分布的均匀性和群落特征梯度性及生态环境的异质性可重复性，同时考虑到可操作性（有些样方的更新苗密集分布，个体数量多，调查困难），本研究对更新苗和成树的调查选择了每个样方左下角的1个5 m×5 m 的小样方为调查基本单元，共255 个调查基本单元，2021年8—9月份调查了每个小样方内更新苗和成树个体，包括测定每个个体的相对坐标（x,y）、株高、冠幅、更新苗的地径和成树的胸径、生长状况（存活、枯死或倒伏），对破损或丢失的林木个体挂牌进行补挂。

1.3 数据分析方法

将样地内采集的小样方更新苗和成树样本数据整理并建立数据集，利用经典统计学和地统计学软件对数据集进行分析处理，在进行数据处理时，将每个调查小样方中心点的坐标作为更新苗和成树密度和高度的采样坐标。经典统计学利用SPSS 19.0 软件进行描述性统计分析和单样本K-S（Kolmogorov-Smirnov）检验，如果K-S 中的Z值大于0.05，则表示该数据集分布符合正态分布，对不符合正态分布的数据集，则利用R 语言包中的forecast 的函数Box-Cox 进行转换，对转换后的数据集再进行K-S 检验，如果样本数据集符合正态分布，则进行下一步分析。地统计学软件利用GS+9.0 软件进行半方差函数、空间自相关、分形维数和Kriging 插值分析，下面简单介绍一下本研究利用的数据分析方法。

1）半方差函数

半方差函数被广泛用来定量描述属性变量的空间结构和空间相关性，包括块金值(CO)、基台值(CO+C)、块基比[CO/(CO+C)]和变程等参数，其具体表征的意义详见文献[16]，根据决定系数选择最佳理论模型[17]，其计算公式为：

式中，γ(h)为半方差函数值，h为两样本间的间隔距离，N(h)为间隔距离为h时的所有观测样点的成对数目，Z(xi)和Z(xi+h)分别为属性变量在位置xi和xi+h位置的实测值。

2）空间自相关

空间自相关用来揭示属性变量的空间自相关性，常用Moran′s I 指数来描述单一值空间变量的相似性，包括全局Moran′s I 指数IN和局部Moran′s I 指数Ii，其具体表征的意义详见文献[16]，其计算公式为：

式中，IN和Ii分别为全局Moran′s I 指数和局部Moran′s I 指数，为属性变量Z 的均值，Zi和Zj分别为属性变量在空间i和j处的数值（j≠i），σ2是属性变量Z的方差，Wij是样本间的距离权重。

3）分形维数

分形维数是在半方差函数分析的基础上得出的参数值，用来定量描述属性变量在各个方向上的复杂程度，其计算公式为：

D=(4-m)/2。

式中，D为分形维数值，m为变异函数曲线图的斜率。D值大小表征的意义详见文献[1-2]。

4）Kriging 插值

Kriging 插值是属性变量获得半变异函数的模拟模型后，对未采样点的属性变量所进行的无偏量的估算，其计算公式为：

式中，z#(x0)为未采样点的属性变量估计值，x0为待插入值点，Z(xi)为变量Z在领域内xi点上的值，λi为权重系数。

2 结果与分析

2.1 更新苗和成树密度、高度的经典统计描述

对青海云杉林动态监测样地的更新苗和成树密度、高度进行了描述性统计分析（表1）。从表1可以看出，样地内更新苗和成树密度均值分别是5 762 株·hm-2和2 016 株·hm-2，对应的 高度均值分别是0.30 m 和8.79 m。从变异系数来看，更新苗密度和高度的变异系数分别为115.62%和120.00%，均属强变异[18]，不同空间样点间更新苗密度和高度的差异很大，在调查的255个小样方中，50%的小样方没有幼苗、幼树，而有的小样方更新株数高达60 000 株·hm-2；而成树密度和高度的变异系数分别为49.70%和50.17%，均属中等变异[18]，不同空间样点之间的变异程度小于更新苗密度和高度的变异程度。更新苗和成树密度、高度的偏度、峰度和K-S 检验表明，除了成树高度符合正态分布外，其他3 个单样本均不符合正态分布，利用函数Box-cox 对3 个单样本进行了转换，转换后的数据参数检验表明更新苗密度和成树密度及更新苗高度数据均符合正态分布。

表1 更新苗和成树密度、高度的描述性统计分析及正态分布检验†Table 1 Descriptive statistical analysis and normal distribution test of the density and height of regenerated and adult trees

2.2 更新苗和成树密度分布的空间格局

2.2.1 半方差函数的分析

对青海云杉林动态监测样地的更新苗密度和成树密度进行了空间结构和变异特征分析。从图1和表2可以看出，更新苗和成树密度的变异函数值分别用球状模型和指数模型描述，能很好地反映更新苗和成树密度的空间结构特征。从表2还可以看出，更新苗密度块基比比值在25%～75%，更新苗密度为中等变异程度的空间相关性[19]；成树密度块基比比值小于25%，成树密度为强烈变异程度的空间相关性[19]。更新苗密度和成树密度的有效变程分别为73.10 和45.60 m，表明更新苗密度的空间自相关范围大于成树空间自相关范围，其生态过程和格局与该尺度范围相关。

表2 更新苗和成树密度变异函数理论模型及其相关参数Table 2 The semi-variogram models and the corresponding parameters of the density of regenerated and adult trees

图1 更新苗和成树密度半方差函数Fig.1 Semi-variogram of the density of regenerated and adult trees

2.2.2 空间自相关分析

对青海云杉林动态监测样地的更新苗密度和成树密度进行了空间自相关系数Moran′s I 分析（图2）。从图2可以看出，更新苗密度和成树密度的Moran′s I 的数值分别是0.269 7 和0.132 6，其变化范围分别是-0.052 4～0.269 7和-0.057 0～0.132 6，更新苗密度的Moran′s I 的数值相对于成树密度Moran′s I 的数值较大，更新苗密度空间依赖性较强，而成树密度空间依赖性较小。更新苗密度和成树密度的Moran′s I 数值均随滞后距的增大而减小，直到下降为负相关并趋于平缓，在这个变化的过程中，出现了拐点，但是其出现的拐点各不相同，更新苗密度在81.51 m 左右出现拐点，而成树密度在134.36 m 左右出现拐点。

图2 更新苗和成树密度的空间相关图Fig.2 Spatial correlation diagram of the density of regenerated and adult trees

2.2.3 分形维数和Kriging 分析

对青海云杉林动态监测样地的更新苗密度和成树密度的各向异性空间变异的分形维数和空间分布格局进行了分析。从表3可以看出，更新苗密度的分形维数是1.936，成树密度的分形维数是1.951。从更新苗密度和成树密度的各向异性即从不同方位来看，更新苗密度分形维数在0°（南北）方向最小，均一性最差，小尺度上的变异更明显；在135°（西北—东南）方向最大，均一性最好，空间自相关弱；从决定系数R2来看，在45°（东北—西南）的值最大，说明该方向是更新苗密度的优势格局。成树密度分形维数在0°（南北）方向最小，均一性最差，空间自相关更强；在90°（东西）方向最大，均一性最好，空间自相关弱；从决定系数R2来看，在全方位的值最大，说明全方位是成树密度的优势格局。

表3 更新苗和成树密度的各向异性空间变异的分形维数Table 3 Fractal dimension of anisotropic spatial variation of the density of regenerated and adult trees

为了更直观地反应青海云杉林动态监测样地更新苗密度和成树密度的空间分布特征，利用半方差函数模型进行了Kriging 插值分析。从图3可以看出（二维结构图），样地更新苗密度和成树密度的空间异质性较强并且变化幅度大，不管是更新苗密度还是成树密度均呈明显的梯度性变化，更新苗密度斑块之间的连接程度低，高密度的更新苗斑块面积所占的比例比较小，破碎度大，呈现明显的聚集性分布格局，而低密度更新苗斑块面积所占的比例比较大，大面积斑块分布的更新苗稀少，其密度甚至为0。成树密度斑块较更新苗密度斑块，其斑块间的连接程度较高，密度低的成树斑块分布相对均匀且占的面积较大，而高密度的成树斑块所占的比例较小，亦呈聚集性分布格局程度。

图3 更新苗和成树密度的等值线图像Fig.3 Contour map of the density of regenerated and adult trees

2.3 更新苗和成树高度分布的空间格局

2.3.1 半方差函数的分析

从图4和表4可以看出，更新苗高度符合指数模型，而非更新高度符合球状模型[16]。从表4还可以看出，更新苗高度和成树高度的块基比比值小于25%，表明样地林木包括更新苗和成树高度为强烈变异程度的空间相关性[19]。2 个变量的有效变程均较小，其中更新苗高度的有效变程最小值为47.40 m，成树高度的有效变程最小值为23.90 m，亦表明更新苗高度的空间自相关范围大于成树空间自相关范围，其生态过程和格局与该尺度范围相关。

图4 更新苗和成树高度的半方差函数Fig.4 Semi-variogram of the height of regenerated and adult trees

表4 更新苗和成树高度变异函数理论模型及其相关参数Table 4 The semi-variogram models and the corresponding parameters of the height of regenerated and adult trees

2.3.2 空间自相关分析

从图5可以看出，更新苗高度和成树的高度Moran′s I 的数值分别是0.178 7 和0.043 7，其变化范围分别是-0.041 9～0.178 7 和-0.036 1～0.043 7，同样是更新苗高度的Moran′s I 的数值相对于成树高度Moran′s I 的数值较大，说明更新苗高度对空间的依赖程度较大，而成树的高度空间依赖性较小。更新苗高度的Moran′s I 的均随滞后距的增大而减小，直到下降为负相关并趋于平缓，在这个变化的过程中，在81.51 m 处出现了拐点，而成树高度随滞后距的增加没有明显的变化规律。

图5 更新苗和成树高度的空间相关分析Fig.5 Spatial correlation diagram of the height of regenerated and adult trees

2.3.3 分形维数和Kriging 分析

从表5可以看出，更新苗高度的分形维数是1.957，成树高度的分形维数是2.000，林木高度在小尺度上的变异更明显。从更新苗高度和成树高度的各向异性即从不同方位来看，更新苗高度分形维数在45°（东北—西南）方向最小，空间异质性最强；在135°（西北—东南）方向最大，空间自相关弱。在45°（东北—西南）的决定系数值最大，东北—西南方向是更新苗高度的优势格局。成树高度分形维数在0°（南北）方向最小，空间异质性最强；在全方位上最大，空间自相关弱。在0°（南北）的决定系数最大，南北方向是成树高度的优势格局。更新苗高度和成树高度的分布优势格局亦不重叠的。

表5 更新苗和成树高度的各向异性分形维数Table 5 Fractal dimension of anisotropic spatial variation of the height of regenerated and adult trees

从图6可以看出（三维结构图），样地更新苗高度和成树高度的空间异质性较强并且变化幅度大，不管是更新苗高度还是成树高度均呈明显的梯度性变化，更新苗高度斑块之间的连接程度低，植株高度较高的更新苗斑块面积所占的比例比较小，破碎度大，呈现明显的聚集性分布格局，而植株高度较低的更新苗斑块面积所占的比例比较大，大面积斑块分布的更新苗稀少。成树高度斑块较更新苗高度斑块，其斑块间的连接程度较高，高度低的成树斑块分布相对均匀且占的面积较大，而植株高度高的成树斑块所占的比例较小，呈现聚集分布格局。更新苗高度和成树高度斑块分布不存在叠加，成树高度的斑块分布格局的这种变化决定了天然更新苗空间分布格局的形成。

图6 更新苗和成树高度的等值线图Fig.6 Contour map of the height of regenerated and adult trees

3 讨 论

1）本研究选择的试验样地是公顷大样地，样地设置范围连续，地形相对平坦，气候条件和成土母质也较为一致，微生境丰富，样本数量也大，采用地统计学能够较好地反应样地更新苗和成树分布及生长的空间格局[16]，利用地统计学测度林木密度和高度2 个参数具有明显的优势，可以表述出这些参数在空间上的变动，同时借助生态理论模型及方法进行尺度推移，更好地反应植被在空间上和时间上的格局变动及其原因。另外，经典描述性统计分析发现，样地更新苗密度和高度属性变量的变异程度（强变异）明显大于成树密度和高度属性变量的变异程度（中等变异），说明更新苗更易受环境因子的影响，如林分结构、微地形、光照、地被物等。

2）从半方差函数模型来看，更新苗和成树密度的理论模型为球状模型或指数模型，反映出更新苗和成树分布呈聚集分布特点。从块基比比值来看，更新苗密度对空间的占据能力中等，说明结构性因素和随机性因素共同起作用，而成树密度对空间的占据能力强，说明引起空间变异的主要因素是结构性因素[20]。这里的结构性因素如微地形、土壤水热、地被物覆盖等自然因素，随机性因素如测定误差、取样设计、干扰等人为因素。从变程上看，更新苗密度空间依赖范围大小为73.10 m，更新苗密度空间上分布存在连续性，表明在此空间范围内任意两个采样点的更新苗密度之间存在一定程度上的相互关系，但是超过这个依赖范围，彼此之间相互作用减弱，并随距离的增加而减小，彼此之间的空间相关系变差。成树密度的空间依赖范围大小为45.60 m，较小于更新苗密度的变程，说明成树密度的空间连续性较更新苗密度强，成树密度决定了更新苗密度的空间分布。更新苗密度的Moran′s I 的数值明显大于成树的Moran′s I 的数值，表明更新苗密度的聚集程度高于成树密度，这可能是青海云杉特殊的生物学特性引起的，非更新密度相对更新苗密度来说具有更大的随机性。正空间相关的距离大致反应了性质相似的斑块平均半径[21]，因此更新苗密度和成树密度存在正相关的尺度为81.51和134.36 m。更新苗和成树密度的分形维数表明，青海云杉林分密度在小尺度的生态学过程作用明显，如青海云杉空间分布除了受自身的生物学特性影响外，还受生长环境的影响。Kriging 分析更新苗密度和成树密度斑块分布不存在叠加，成树密度的斑块分布格局的这种变化决定了天然更新苗空间分布格局的形成。

3）更新苗和成树高度的变异函数的理论模型为指数模型和球状模型，表明林分林木的高度分布呈规律性的空间变异。更新苗和成树高度的块基比比值均小于25%，表明其空间变异均主要由结构性因素引起，青海云杉更新苗在萌生苗、幼苗阶段，可以忍受一定程度的庇荫，或者说在荫蔽的环境中更易存活，但是随着年龄的增长，从幼苗生长到幼树的这个阶段过程中，郁闭的林分不能满足幼树生长对光照为主的因子需要，通过环境的过滤作用，造成了幼树的大量死亡，径级和龄级较大的幼树很少存在，在祁连山青海云杉形成了“只见幼苗不见幼树”的特色现象。更新苗和成树高度的空间结构比大于90%，其格局主要由空间自相关因素引起的，成树高度的空间自相关范围为23.90 m，其值小于更新苗高度的自相关范围47.40 m，同时成树高度Moran′s I 的数值相对更新苗高度的Moran′s I 的数值较低，但其值在±0.043 7 之间，其数值接近0，变异函数的参数值变程和空间自相关拐点值表明成树高度的空间分布和形成过程决定了更新苗的空间分布和形成过程，进而决定了样地林木的空间分布格局[13,22]。分形维数和Kriging 分析表明样地高度和分布的空间格局动态过程明显，更新苗和成树高度在不同方向上的分形维数表现出各向异性，而且二者之间的优势格局不一样，这和林木种子扩散、种内竞争、环境异质性等在各方向上的差异有关[23]。Kriging 分析得出更新苗高度分割的斑块空间面积相对较小，不同大小的斑块在空间上呈现出不连续的分布状态，而成树的高度分割的斑块空间面积相对较大，不同大小的斑块在空间上呈现出连续的分布状态，更新苗和成树高度斑块大小及其空间镶嵌结构的复杂程度表明，在进行植被结构研究时要注意小尺度上的生态学过程[9]，即小尺度的环境空间异质性对更新格局有着重要的影响，如微生境对更新分布和生长的影响。

林木的生长过程是一个动态的变化过程，而且在变化的过程中，林木对环境的选择包括生物和非生物因素存在差异，本研究对更新苗和成树的空间异质性分析只是通过一次调查所进行的分析，需要通过长期动态监测从时间的空间异质性角度进行分析。再者本研究只是基于动态监测样地即大样地进行的青海云杉林更新苗和成树的空间异质性分析，而忽略了海拔梯度上分布的青海云杉林更新苗和成树的空间异质性分析，需要将大样地水平尺度的空间异质性和海拔梯度的空间异质性结合起来进行研究。因此，需要加强时间尺度和空间尺度上的异质性研究，并在此基础上预测青海云杉林生态系统的格局和过程具有重要的意义。

4 结 论

从林木参数密度的分布格局来看，更新苗和成树空间分别表现为中等变异程度和强烈变异程度的空间相关，其有效变程分别为73.10 和45.60 m；更新苗密度空间依赖性较强，而成树密度空间依赖性较小，其Moran′s I 的数值分别是0.269 7 和0.132 6；更新苗密度的分形维数是1.936，成树密度的分形维数是1.951，更新苗和成树密度的优势格局分别是45°和全方位，更新苗的聚集程度比同尺度上成树的聚集程度要高，而且更新苗的空间连续性要比成树的空间连续性低。从林木参数高度的分布格局来看，更新苗和成树空间均表现为强烈变异程度的空间相关性，其有效变程分别为47.40 和23.90 m；更新苗高度空间依赖性较强，而成树高度空间依赖性较小，其Moran′s I 的数值分别是0.178 7 和0.043 7；更新苗高度的分形维数是1.957，成树高度的分形维数是2.000，更新苗和成树高度的优势格局分别是45°和0°，更新苗空间分布格局较成树空间格局具有明显的显斑块性，复杂性程度不高，连接性程度较低。

综上所述，祁连山青海云杉林动态监测样地林木密度和高度均呈现出明显的空间变异格局，并且这种空间变异格局表现出明显的动态变化，林分环境异质性对更新格局有着重要的影响。