21世纪以来，在全球气候变化的大背景下，火灾、干旱和病虫害等干扰的频率不断增加，造成的破坏逐渐加剧，严重威胁着森林生态系统的持久性^[1]。森林生态系统中的多种不同生态学过程(如种子扩散、种内和种间竞争、干扰、环境异质性等)直接或间接影响着植物的空间格局，因此，目前空间格局研究的主要内容集中在认识格局形成的生态学过程^[2]。以往的研究大多集中在因树木死亡造成的分布格局变化对生态系统稳定性及持久性的影响^[3]、森林群落中不同林分下幼树幼苗更新分布格局或不同树种的空间分布格局变化^[4]以及林窗下小尺度空间的天然更新模式^[5]等，较少关注林分的天然更新对生态系统持久性的影响。作为北方针叶林的主要组成树种，欧洲赤松(Pinus sylvestris L.)的天然更新也受到了一定的关注^[6]。在我国，欧洲赤松变种的樟子松(Pinus sylvestris L. var. mongolica Litv.)主要分布在大兴安岭北部山地及呼伦贝尔草原向大兴安岭过渡带的沙地上^[7]，其中呼伦贝尔沙地是其最为集中的天然分布区，主要林分为1950年以来封山育林及禁伐禁牧后自然恢复形成的中龄林^[8]。同时由于樟子松的抗旱和耐沙等特性，在我国北方干旱半干旱区的造林和植被恢复中得到广泛应用^[9-10]，但近年来大面积的樟子松人工林面临着衰退的危险^[11]，而较差的幼树天然更新状况则进一步削弱了林分的稳定性和持久性^[12]。以往有关幼树更新影响因子的研究大多集中在生物因子，如同种或不同种的大树对幼树生存率和保存率的影响上^[13-14]，而对非生物因子如地表扰动状况特别是对林分中广泛分布的蚁穴对幼树更新影响的研究尚不多见。鉴于此，选择樟子松原产地呼伦贝尔沙地红花尔基镇的天然樟子松纯林为研究对象，对其林分的更新状况及其与林分中大树及地表广泛分布的蚁穴间的相互作用关系进行研究，通过天然林的更新研究结果探求人工林更新的经营管理措施，从而为我国干旱区大规模的樟子松人工林的科学经营管理提供依据。

1 研究方法 1.1 研究区概况

研究区为呼伦贝尔沙地红花尔基镇樟子松天然林分布区，该地位于内蒙古自治区呼伦贝尔沙地南端、大兴安岭西坡中部向内蒙古高原的过渡带上，北纬47°36′~48°35′，东经118°58′~120°32′；该地区为半湿润半干旱区，年平均温度-1.5 ℃，大于10 ℃的年积温为2 000 ℃，年平均降水量344 mm，主要集中在7—8月份，年平均蒸发量为1 174 mm，干燥度为1.33，平均无霜期90 d左右；主要地貌类型为垄状、波状起伏的沙地，由厚度不等的沙层和沙砾层为主的第四纪冰川沉积物组成，海拔高度在700~1 100 m之间^[15]；主要植被类型以温带草甸草原和温带丛生禾草草原为主，东南部分布有兴安落叶松[Larix gmelinii (Rupr.) Kuzen.]和樟子松等温带山地针叶林和山杨(Populus davidiana Dode)、白桦(Betula platyphylla Suk.)、花楸[Sorbus pohuashanensis (Hance) Hedl.]等温带落叶阔叶林^[16]。

1.2 野外调查方法

2015年6—8月，在研究区天然樟子松中龄纯林中，选取一块天然更新状况较好的面积1 hm² (100 m×100 m)的正方形样地作为典型调查样地，该样地处在沙丘间较为广阔的平地上。为便于野外调查，将样地分割成5 m×5 m的方格，以样地一个顶点为原点，调查记录样地中大树和幼树的位置坐标，其中大树树高≥1.5 m，有可以测定的胸径(diameter at breast height，DBH)，幼树树高 < 1.5 m，没有可以测定的DBH，同时对每一株幼树的树龄进行了记录。

1.3 数据分析 1.3.1 幼树存活率及保存率

首先将幼树根据年龄划分为4组，即1年生、2年生、3年生和3年生以上，对其数量及高度进行比较分析。考虑到呼伦贝尔沙地气象站点非常少，为了分析水分因子对幼树存活率和保存率的影响，采用高精度的TerraClimate全球气候水文数据序列(1958—2015年，分辨率4 km)获取样地所在栅格近4 a的降水量及水分亏缺数据^[17]，通过使用Programita软件完成计算分析。

1.3.2 单变量空间格局分析

将所有的树木根据其坐标位置概化为二维平面上的一个点，同时根据蚁穴所在多边形质心的坐标将其概化为一个点，然后采用空间点格局的分析方法对所有点进行单变量分析。采用成对相关函数g(r)对大树、幼树和蚁穴的空间格局进行分析。函数g(r)的单变量格局的计算公式如下：

$ \hat g(r) = \frac{1}{{2{\rm{ \mathsf{ π} }}r}} \cdot \frac{{{A^2}}}{{{n^2}}}\sum\limits_{i = 1}^n {w_{ij}^{ - 1}} {k_h}\left( {r - \left| {{x_i} - {x_j}} \right|} \right) $

(1)

式中：r为点对间的距离(m)；A为样地面积(m²)；n为样方中个体的数量；w_ij为边缘校正的权重因子；k_h为Epanecnikov核函数，这是一个权重函数，当点对间的距离等于r时被赋值为1；h为圆环宽度参数，可采用公式h=c (n/A)^-0.5进行计算，其中0.1≤c≤0.2，x_i和x_j分别代表点i和j所处位置，j≠i。将圆环的宽度确定为1 m。在某一距离r处，g(r)>1表示点对处于聚集状态，而g(r) < 1则表示点对处于规则分布状态，g(r)=1则表示两者间不存在相互作用关系。采用完全空间随机零模型对樟子松大树和蚁穴的空间分布进行分析，考虑到幼树表现出明显的空间聚集状态，采用托马斯过程(Thomas process)零模型对幼树的空间分布进行分析^[18-19]。

1.3.3 双变量格局分析

成对相关函数g(r)也可以进行双变量点格局分析，即计算以物种1的每一个个体为圆心、半径为R、指定宽度圆环区域内物种2的数量，计算公式如下：

$ {{\hat g}_{12}}(r) = \frac{1}{{2{\rm{ \mathsf{ π} }}r}} \cdot \frac{{{A^2}}}{{{n_1}{n_2}}}\sum\limits_{i = 1}^{{n_1}} {\frac{{{n_2}}}{2}} w_{ij}^{ - 1}{k_h}\left( {R - \left| {{x_i} - {y_j}} \right|} \right) $

(2)

式中：n₁和n₂分别代表第1类和第2类的点数；x_i (i=1, 2, …，n₁) x_i和y_j (j=1, 2, …, n₂)分别为两类点所处的位置。g₁₂(r) >1表示两类点间为正相关，g₁₂(r) < 1表示两类点间为负相关，g₁₂(r) =1则表示两类点间不存在相互作用关系。

为了分析樟子松幼树与大树、蚁穴间的关系，选用环形移动零模型^[20-21]对两类样本间格局的独立性进行非参数检验，模型中将研究样地看作环形，然后样本类型1保持不动，而样本类型2围绕样本类型1随机移动，g₁₂(r)值大于或小于模拟包迹线分别表示两类样本间相互吸引或排斥。根据野外调查的结果，蚁穴的短轴半径平均为59 cm (最大154 cm，最小26 cm)，为了更为准确地揭示其与幼树间的关系，本研究将分析的圆环宽度分别设定为0.2 m和1 m。

1.3.4 三变量随机标记分析

为了检验樟子松大树及蚁穴共同存在时对樟子松幼树的影响，分别分析幼树存活的概率与樟子松大树及蚁穴间距离的依赖关系，即检验是否已有的格局会对幼树的存活产生影响。此处采用的分析方法为三变量随机标记法^[22]，在该方法中传统的随机标记被用于樟子松的大树(下标2)与幼树(下标1)间的格局分析，但统计检验被用来估计幼树存活的概率与蚁穴的距离或大树+幼树(下标1+2)间的依赖关系。计算双变量的蚁穴与幼树间的邻体密度和蚁穴与樟子松个体(大树+幼树)间的邻体密度，并采用两者间的比值来反映幼树的存活概率，具体计算公式如下：

$ \rho_{h, 1}(r)=\frac{\lambda_{1}}{\left(\lambda_{1}+\lambda_{2}\right)} \cdot \frac{g_{h, 1}(r)}{g_{h, 1+2}(r)} $

(3)

式中：λ₁和λ₂分别为樟子松幼树和大树的强度系数，由Programita软件计算得出; g_{h, 1}(r)和g_{h, 1+2}(r)分别为蚁穴与幼树间及蚁穴与樟子松个体间的双变量g(r)函数，当ρ_{h, 1}>λ₁/(λ₁+λ₂)时，表示蚁穴的存在对幼树的影响为正向(相互促进)的，反之则是负向(相互抑制)的。

采用相似的检验过程分析了樟子松大树对幼树存活概率的影响，具体计算公式如下：

$ \rho_{1+2, 1}(r)=\frac{\lambda_{1}}{\lambda_{1}+\lambda_{2}} \cdot \frac{g_{1+2, 1}(r)}{g_{1+2, 1+2}(r)} $

(4)

式中：g_{1+2, 1}(r)和g_{1+2, 1+2}(r)分别为樟子松大树与幼树间及所有樟子松个体间的双变量g(r)函数。当ρ_{1+2, 1}>λ₁/(λ₁+λ₂)时，表示大树的存在对幼树的影响为正向的，反之则是负向的。

为了检验上述零模型下函数计算结果的显著性，采用零模型的199次蒙特卡罗模拟产生95%置信区间的模拟包迹线，所有空间格局分析均采用点格局分析软件Programita完成^[23]。柱状图及线形图用Origin 9.0绘制，单因素方差分析(One Way-ANOVA)用SPSS 21.0软件计算。

2 结果与分析 2.1 幼树的成活率及保存率分析

图 1为不同树龄樟子松幼树(1~4 a)高度、株数特征及过去4 a水分亏缺和降水量分布状况。从不同幼树的数量统计结果来看[图 1(a)]，林分中1年生和2年生幼树的数量很大，这表明林分中樟子松的种子萌发和成活率还是很高的，但3 a及以上的幼树的数量则显著下降，显示幼树的保存率并不高，这会严重阻碍林分的更新演替；而从不同幼树的高度的统计结果来看[图 1(b)]，随着樟子松幼树年龄的增加，幼树高度呈显著性递增(P < 0.05)，说明在幼苗的更新过程中，虽然幼树的保存率并不高，但成活的幼树是可以正常成长发育的。考虑到调查区域已经实施禁牧，因此，由放牧采食和牲畜践踏而导致幼树死亡是可以排除的；从过去4 a的降水量及水分亏缺量来看，降水量基本维持在350 mm以下，而水分亏缺量为300~450 mm[图 1(c)]，可见持续干旱是幼树保存率不高的原因之一。

2.2 林分中树木个体及蚁穴空间分布单变量分析

基于托马斯过程零模型的樟子松幼树空间分布的单变量分析结果如图 2所示。总体上看，林分中所有樟子松幼树在40 m的尺度范围内均呈显著的聚集分布状态[图 2(a)]，表明在林隙和微地貌的共同作用下，林分内形成的资源异质性决定了幼树的团块状聚集分布。对不同年龄幼树的分析结果则显示，1年生和2年生的幼树在约15 m的尺度范围内呈显著聚集分布，3年生的幼树成随机分布状态，3年生以上的幼树仍呈现一定范围内的聚集分布[图 2(b~e)]，出现这种现象可能是由于幼树过度聚集，导致对资源的强烈竞争，3年生的幼树开始大量死亡，导致幼树的分布整体上呈稀疏状态，从而表现为一定的随机分布，随着稀疏度的逐步加大，又进一步呈现稀疏的聚集分布，只是尺度有所缩小。同时近几年均存在不同程度的干旱，其累加效果加剧了幼树间对水分的竞争，从而导致幼树大量死亡。

基于完全随机零模型对林分内樟子松大树和蚁丘的空间分布结果如图 3所示。在4 m的尺度范围内，樟子松大树呈显著的聚集分布[图 3(a)]，而蚂蚁所形成的蚁丘在所有尺度上均呈现随机分布[图 3(b)]。

2.3 幼树与大树及蚁穴间关系的双变量分析

幼树与大树及蚁穴间关系的双变量分析结果如图 4、图 5所示。樟子松幼树与大树之间总体上并不存在着显著的作用关系[图 4(a)]，只是随着幼树年龄的增长，在较小尺度上出现了一定程度的竞争，但这种竞争并不显著[图 4(b~e)]；表明在目前的林分结构下，种内竞争并不是导致幼树死亡的主导因子；而蚁穴的存在在较小尺度上(5 m)对樟子松幼树产生了一定的抑制作用[图 5(a)]，特别是对1年生幼树而言，尽管这种作用并不十分显著[图 5(b)]，究其原因，可能是围绕蚁穴周围所形成的蚁丘土壤较为疏松，在高温作用下迅速变干，从而导致根系较浅的幼树无法获取到生长所需要的水分，出现萎蔫乃至死亡的现象。

2.4 幼树存活与大树和蚁穴间协同作用间的依赖关系

三变量随机标记分析的结果(图 6)显示，在樟子松大树存在的情况下蚁穴对幼树的影响与双变量分析的结果相似，在较小尺度上对幼树的存活起到了不显著的抑制作用，由于樟子松大树的存在，这种抑制作用的尺度进一步缩小为0~2 m，而在3~4 m的尺度上则表现为一定程度的促进作用[图 6(a)]；而在蚁穴存在的情况下，樟子松大树对幼树的影响则与双变量分析的结果存在一定的差异，其结果与蚁穴的影响较为接近，即在0~2 m尺度上抑制，3~4 m尺度上促进，但其抑制作用稍弱于蚁穴[图 6(b)]。

3 讨论与结论

樟子松更新苗分布的空间格局与该物种的生长习性和亲代的散布习性密切相关，由物种自身的内在适应性决定，其中单株母树的龄级、生活力、冠幅、树高和所处的林分状况以及种子的重量、生活力都是影响樟子松幼树成活与保存的形态特征因素^[24]；除此之外，有研究表明樟子松天然林的生长近半个世纪以来一直处于暖干化气候下的胁迫之中，并导致出现生长衰退的迹象^[25]，长期、持续的干旱对树木的生理生态过程具有重要的作用，可以影响和改变树木的生活史进程^[26]，同时对樟子松的径向生长有着重要影响^[27]，因此可以从降水在时间和空间上丰亏变化的角度解释樟子松成树衰退及幼树更新的水分驱动机制^[28-29]。

本研究结果表明，林分中樟子松幼树在较大尺度内(40 m)呈显著的聚集分布状态，大树在较小尺度内(4 m)呈显著聚集状态，而蚁穴则在所有尺度内呈随机分布状态。对同种的樟子松幼树和大树而言，两者的关系在不考虑其它影响因子(蚁穴)的情况下总体上不存在着显著的作用关系，而在考虑蚁穴作用的情况下则在0~2 m尺度上表现为抑制作用，3~4 m尺度上表现为促进作用，但两者均不显著；对樟子松幼树和蚁穴而言，两者的关系在不考虑其它影响因子(樟子松大树)的情况下，在较小尺度上(5 m)对樟子松幼树产生了一定的抑制作用，而在考虑樟子松大树作用的情况下，不显著的抑制作用进一步缩小为0~2 m尺度上，而在3~4 m的尺度上则表现为不显著的促进作用，其影响略强于樟子松大树。结合幼树成活率和保存率的分析结果，可以得出天然樟子松林幼树保存率低，在很大程度上取决于气候因子特别是干旱带来的水分亏缺的影响。

在以往的许多研究中，同种或不同种的树木个体间的空间格局分析大多采用单变量和双变量的分析方法^[30]，而没有考虑其它外界因子对两两作用的影响，因此得到的结果可能与真实情况间存在着一定的偏差，本研究引入了三变量随机标记分析方法，可以更为准确地反映樟子松大树和蚁穴对樟子松幼树更新格局的影响。

此外，本研究将林分中广为分布的蚁穴作为幼树存活的主要影响因子进行分析，可以较为全面地考虑了森林生态系统的多样性和复杂性，可以获取到更为科学、可信的研究结果，这对天然及人工樟子松林的可持续经营管理具有重要的理论价值。