研究区位于内蒙古牙克石市乌尔旗汗林业局，地处大兴安岭西坡，该区为寒温带大陆性季风半湿润森林气候，森林资源丰富，总面积59.36万hm²，森林覆盖率为78.97%，主要树种有落叶松(Larix gmelinii (Rupr.) Kuzen.)、白桦(Betula platyphylla Suk.)和山杨(Populus davidiana var. davidiana f. pendula (Skv.) C. Wang et tung)等。

2018年7月，采用空间代替时间法，在火烧干扰后恢复12 a、22 a和31 a以及采伐干扰后恢复14 a、20 a和29 a的森林生态系统中随机布设样地，研究区内无人为影响，所有乔木均自然恢复，样地大小为20 m×20 m。火烧迹地按照样地内胸径大于>5 cm且有明显火烧痕迹的枯立木株数来划分火烧程度，不考虑后期自然恢复的植被由于其他原因导致的死亡。采伐迹地的采伐方式为皆伐。干扰类型一共有4种（轻度火烧迹地、中度火烧迹地、重度火烧迹地和采伐），每种类型迹地有3个重复，共计36块样地，恢复时间是干扰发生时间距调查时（2018年）的时间间隔。

调查内容包括各样地的GPS坐标、海拔、坡度和坡向，对所有胸径≥1 cm的木本植物每木检尺，记录物种名、胸径、树高和死木信息，调查数据中共计样木2 692株，其中活立木2 319株，枯立木373株，本研究只针对活立木，数据处理时剔除了每个样地的枯立木。

利用每木检尺数据，得到各样地的平均胸径（d）、平均树高（h）、株数密度（N）和树种组成（优势树种）等，树种组成按样地内某一树种胸高断面积占样地断面积的百分比确定，恢复时间可近似认为是林分的平均年龄（表1）。

非生物因子（环境条件）主要包括地形因子、气候因子以及干扰类型。

地形因子包含3个变量：海拔、坡向和坡度。36个样地的海拔为782~990 m，地势较为平坦；坡度按定性变量，划分标准为：Ⅰ级为平坡（<5°），Ⅱ级为缓坡（5°～14°），Ⅲ级为斜坡（15°～24°），Ⅳ级为陡坡（25°～34°），Ⅴ级为急坡（35°～44°），Ⅵ级为险坡（≥45°），其中仅5个样地为Ⅲ级（斜坡），其余大都为平坡或缓坡；坡向由0~360°转换成0-1之间的值，公式^[17]如下：

式中：P为坡向指数，A为坡向角度，转换后的数值越小表明生境越湿冷，反之越干热。

气候数据从ClimateAP中提取^[18]，气候因子变量有：年平均温度、年平均降水、月平均最高温、月平均最低温、近30年极端最低气温等，计算时取干扰发生时间至2016年期间的平均值作为各样地的气候因子。1987—2016年样地的年平均气温和年总降水量变化趋势并不明显，年平均气温在−2.3 ℃左右波动，年总降水量在430 mm左右波动。

此次调查所涉及的树种有落叶松、白桦、山杨、云杉(Picea asperata Mast.)、柳树(Salix matsudana var. pseudo-matsudana (Y. L. Chou et Skv.) Y. L. Chou)、稠李(Padus racemosa (Lam.) Gilib.)等，使用相关生物量模型结合各组分含碳系数计算碳储量。对于活立木的碳储量计算原理如下：

式中，c为单木碳储量；B_i为第i分项的单木生物量，P_i为第i分项的含碳系数， stem、bark、branch、leaf、below分别表示单木干材、树皮、树枝、树叶和地下等组分。（2）式充分考虑了林木不同器官含碳率的差异，可以推算单株活立木的碳储量。

基于每木检尺的胸径、树高，使用二元生物量模型估算单木各组分生物量，结合各树种不同组分的含碳系数计算单木的碳储量，其中云杉、落叶松、白桦以及山杨采用行业标准《生物量模型及碳计量参数》中东北和内蒙古东部地区的二元生物量模型以及各组分的含碳系数计算；柳树和稠李等树种按软阔类，采用《中国森林植被生物量和碳储量评估》中的全国生物量模型与软阔类全树平均含碳系数计算^[19]。林分碳储量C由样地内单木碳储量的结果累加得到。

使用逐步回归^[20]的方差分析法。具体过程为：首先将所有因子作为自变量，以林分碳储量作为因变量，使用方差分析法得到所有因子的影响，然后按每个因子P值的大小排序，删除P值最大的因子（即最不显著的因子），接下来用剩余因子与林分碳储量继续做方差分析，重复上述步骤，直至所有影响因子显著（P<0.05）为止，方差分析在SPSS软件进行。

结构方程模型（SEM）方法中，首先假设所有因子均对林分碳储量有影响，并考虑变量间可能出现的相互影响来建立初始结构方程模型，然后估计模型中的路径系数，其显著性用t检验判断，若路径系数不显著，则删除该条路径，重新拟合模型，再根据模型修正指数（M.I）优化模型。使用卡方检验、渐进残差均方和平方根(RMSEA)、适配度指数（GFI）、规准适配指数（NFI）、Akaike讯息效标（AIC）等评价指标判断结构方程模型的拟合优度，SEM估计在Amos软件中完成。

表2所示，平均胸径和株数密度对林分碳储量的影响最显著（P<0.001），平均树高对林分碳储量没有显著影响，使用d和N即可很好的描述林分碳储量；月平均最高温和月平均最低温对林分碳储量有影响（P<0.05），表明林分碳储量对极端气候敏感；干扰类型对林分碳储量有影响，比较不同干扰类型下林分碳储量的差异（表3）发现，采伐迹地、中度火烧迹地和重度火烧迹地的林分碳储量没有显著性差异，而轻度火烧迹地的林分碳储量与它们有显著性差异；地形因子（海拔、坡向、坡位）和恢复时间在逐步回归的过程中被剔除，其对林分碳储量没有影响，地形因子无显著影响的主要原因是各样地的地势比较平坦，生境条件较为相似；至于恢复时间，由于不是连续观测的数据以及样本数量的限制，可能导致恢复时间对林分碳储量没有显著性影响。

建立林分碳储量的通用模型时，参考林分碳储量多重比较的结果，将干扰类型（I）设置为2个水平，模型如（3）式所示，确定系数为R²=0.890（调整后的R²=0.872）。林分碳储量随d和N的增加而增加，月平均最高温和月平均最低温均对林分碳储量分别为负相关和正相关的影响，林分碳储量随温度的变化说明在极端气温范围内，有着林分碳储量的适宜生长区间，与轻度火烧迹地相比，采伐迹地、中度和重度火烧迹地的林分碳储量要小，恢复起来更加困难。

上式中，C为林分碳储量，d为林分平均胸径，N为株数密度，Te_max为月平均最高温度，Te_min为月平均最低温度，I为干扰类型。

SEM最终模型的拟合结果很好（表4），模型经卡方检验，显著性概率p>0.05，未达显著水平，接受虚无假设，表明该模型能够契合样本数据的结构，其余适配指标均达到适配标准，AIC的结果表明模型没有出现过拟合现象，各变量间的路径系数通过了t检验，可以使用该模型的结果探究林分碳储量的驱动因素。

SEM的结果反映了林分因子和非生物环境因子对林分碳储量的关系（图1），模型的拟合优度略低于多元线性回归的结果（SEM中林分碳储量部分的模型确定系数为0.757）。SEM的结果可以进一步解析为各因素对林分碳储量的直接影响、间接影响和总影响：d和N对林分碳储量有直接正向影响，N还对林分碳储量有负向间接影响，通过影响d间接影响林分碳储量，而月平均最高温只对d有负向的直接影响，对林分碳储量只有负向的间接影响（−0.52），表5是各因素对林分碳储量的影响汇总，d对林分碳储量的影响最大（0.94），而N对林分碳储量的影响最小（仅0.12），月平均最高温对林分碳储量的影响介于平均胸径和株数密度之间。

图  1  林分碳储量结构方程模型

Figure 1.  Result of structural equation model for stand carbon

本研究分别使用方差分析法和结构方程模型法判断了退化森林生态系统中林分碳储量的驱动因素。从研究方法上来看，方差分析容易实现，很多软件都可以完成，也是大多数学者比较熟悉的方法，而SEM通常需要有理论支撑、经验法则或实验基础，从而得到合理的假设模型，最后才能进行模型参数的估计。方差分析对数据的要求比较高（如独立性、正态性），而SEM不对数据有任何分布假设的限制。关于样本数量，张璇和王嘉宇^[21]认为线性模型中，样本数量取30就可以得到较为可靠的估计值，SEM模型的估计允许样本少于60个^[22]。本研究有36个样本，2种方法的估计结果均良好。其中多元线性回归模型的R²为0.890，SEM方法的R²为0.757（以林分碳储量作为因变量），SEM略低于多元线性回归模型，但是SEM能够考虑变量间的相互作用，并解析为直接影响和间接影响，更能体现林分碳储量的驱动因素，方差分析则得到各变量对林分碳储量产生的独立影响，没有考虑变量间的相互作用。

从研究结果上看，方差分析和SEM的结果均表明林分因子对林分碳储量的影响最大，其中起主要作用的变量为平均胸径和株数密度，气候因子对林分碳储量有影响，但两种方法挑选出的变量有所不同。对于不同的干扰类型，方差分析的结果显示不同干扰类型对林分碳储量有影响，而SEM的结果显示干扰类型对林分碳储量没有影响。两种方法均表明地形因子和恢复时间对林分碳储量没有显著性影响。虽然两种方法的计算原理并不相同，但结果相似。方差分析法操作简单，精度较高，但缺乏对多元驱动机制的深入解释，SEM方法更能揭示多变量间的结构关系，但需要有理论支撑或经验基础才能得到结果。

从林分因子上看，平均胸径是一个重要的林分特征指标，也是野外调查中较为容易且测量准确的因子，许多研究均表明其与林分碳储量间有良好的统计关系^[23-24]，在本研究结果中也得到了进一步印证。干扰主要导致森林郁闭度降低，产生林间空隙或形成林窗，林窗的形成则是因为原林分的大树倒下，极大地降低了生态系统中乔木层的生物量^[25]，进而影响林分碳储量，林窗大小影响了光照条件，使得迹地内林木的萌发数量（即林分密度）有所不同，影响了林分碳储量的积累，李海奎等^[26]的研究结果表明株数密度对碳计量参数有显著影响，进而影响林分碳储量，与本研究2种方法得到的株数密度对林分碳储量有正向影响的结果一致。Fan等^[27]的认为阔叶树的生物量大于针叶树的生物量；由针叶林向针阔混交林转变时，碳储量平均增加16.8%^[28]，这说明不同林分类型的碳储量的累积能力有所不同^[29]，然而本研究的研究结果表明优势树种对林分碳储量没有显著影响，但质量比假说（Mass Ratio Hypothesis）则认为生态系统的功能是由群落中优势物种的特征驱动的^[30]，也就是说优势树种能决定林分碳储恢复量的大小。但在大兴安岭的退化森林生态系统中，白桦是先锋树种，干扰后的萌生种多为丛生白桦，林分属性发生改变，萌发的白桦长得快，老化得也快，萌生白桦林寿命短，形成许多小老树，进而导致林分优势树种以及恢复时间对林分碳储量没有显著影响。

在非生物环境因子中，许多学者研究了气候因子与林分生物量（碳储量）的关系，Condes等^[31]在月平均气温和月平均降水等10个生物气候指标中挑选出了对林分生长有显著影响的变量，最后的结论认为气象因子对林分生长的影响效果不明显，但在本研究中，气候因子对林分碳储量的标准化影响达50%以上，气候对林分碳储量有较大的影响，这与气候数据的来源、挑选的气候变量以及样本数量有关。在地形因子中，谭珊珊等^[8]的研究结果表明海拔对林分地上生物量有负向影响，王轶夫^[32]的研究结果也认为森林生物量跟海拔、坡度和坡向之间存在明显的相关性，这与本研究结果并不相同，其原因是本研究区的地势都比较平坦，地形差异并不大。但值得注意的是，地形与气候之间也有十分紧密的关系，比如随着海拔的升高，气温会逐渐降低，所以气候信息也代表了的部分地形信息，鉴于36块样地的地理位置较为接近，气候差异、空间差异可能未得到充分体现。从干扰类型上看，张玉红等^[33]的研究结果认为林分生物量随着火烧程度的增加而降低，陈永富等^[34]的研究结果表明过度采伐不利于森林碳储量的积累，与本研究结果有所不同，这与萌生白桦的生长习性、干扰前的林分状态以及样本数量有很大关系，同时考虑到退化森林生态系统功能恢复机制的复杂性，最终导致了非生物环境因子对林分碳储量的驱动因素与其他研究结果的不同。

由于不同研究的样地布设方法^[37]、单木生物量模型形式和适用地区以及含碳系数有所差异，导致林分碳储量推算的精度有所不同，单木生物量模型的精度则对林分碳储量的估算有着最直接的影响，而使用0.5作为“生物量—碳储量转化系数”则会进一步影响林分碳储量估算的精度。因为林分碳储量估算方法的差异，也会导致林分碳储量驱动因素的不同，虽然本研究得到了不错的模型拟合效果，但受到数据获取方面的困难，样本数量仍然偏少，结果还需要进行验证，若增加样本数量应该能进一步改善模型拟合效果，深入探究退化森林生态系统中林分碳储量的驱动因素。

方差分析法和SEM法都可以对林分碳储量的驱动因素做出有效判断，2种方法的结果基本一致：退化森林生态系统中，林分因子和气候因子是林分碳储量的主要影响因素，而地形因子则无显著影响，平均胸径、株数密度和月平均最高温是2种方法都挑选出的对林分碳储量有显著影响的变量。