随着CO₂等温室气体浓度不断升高, 全球气候变暖趋势日益加剧, 对当今人类生存环境提出了严重的挑战^[1-2]。作为陆地生态系统最大的碳库^[3], 森林生态系统具有巨大的碳汇能力, 充分发挥森林的碳汇作用是有效缓解全球气候变暖的重要手段之一, 其关键和前提是准确估算和掌握森林生态系统的碳储量、碳密度、净初级生产力(NPP)及其响应机制。目前国内外学者已针对森林碳储量及NPP开展了大量的研究^[4-8], 但由于数据来源、估算方法各异及植被类型的多样性, 致使估算结果差异较大^[9-10], 为更精准的估算森林碳库的大小及其固碳能力, 需要针对特定区域、特定树种进行更小尺度的碳储量及NPP的精准评估。

近年来, 落叶松林生长及分布格局的变化正受到越来越多的关注。有研究表明：气候变暖引起的区域环境暖干化使得限制兴安落叶松径向生长的主导因子正由温度转化为水分, 如果气候继续变暖, 兴安落叶松可能在我国境内消失^[11-13]；干旱程度的加重使得蒙古大部分地区落叶松生长区域缩小^[14]；新疆西伯利亚落叶松出现“高固碳速率北移, 低固碳速率南移”的趋势^[15]。阿尔泰山落叶松林作为林区内分布最为广泛的植被类型^[16]和阿尔泰山泰加林典型代表之一, 具有涵养水源、保持水土、防风固沙等生态功能^[17], 天然更新能力很差, 林区内林木天然更新多呈现云杉取代落叶松的趋势^[18], 且其固碳速率空间分布格局正逐渐由“东南高，西北低”向“西北高，东南低”转变^[15]。阿尔泰山落叶松林生长区域、固碳速率的变化必将对整个阿尔泰山森林生态系统的碳循环及碳收支平衡产生重大影响。当前有关阿尔泰山碳储量及NPP的估算研究多以阿尔泰山森林生态系统为研究对象, 着重研究了整个阿尔泰山的固碳特征^{[8, 19-20]}, 其研究成果无法清晰揭示林区内特定树种尤其是落叶松林固碳特征的时空变化规律, 忽略了不同树种间相互影响、相互制约对阿尔泰山森林生态系统固碳的影响。

现阶段, 针对阿尔泰山落叶松林的研究多集中于树轮、生长量与气候变化等方面^[21-23], 针对落叶松林固碳机制的研究不多, 尤其缺乏针对阿尔泰山落叶松林碳储量和NPP时空分布格局及其气候成因的深入探讨。本文利用2001—2016年的新疆多期森林资源清查资料, 估算阿尔泰山落叶松林碳密度、碳储量及NPP, 探讨其时空分布格局和气候成因, 旨在揭示阿尔泰山落叶松林固碳规律及机制, 进而为有关部门进行科学决策提供有效的理论和决策支持。

1 材料与方法 1.1 研究区概况

阿尔泰山位于新疆北部边界, 全长约2000km, 山体走向为西北-东南走向。地形呈现出“西北高而宽, 东南低而窄”的特征；气候上属温带大陆性气候, 夏季温暖多雨, 冬季寒冷干燥。阿尔泰山林区森林资源丰富, 是我国唯一的南西伯利亚区系原始森林, 多分布于海拔1200—2600m的山地地区, 多年平均气温为-2℃左右, 年降水量为400—600 mm, 随海拔的升高而递增, 并由东南向西北递增^[19]。林区内植被类型主要以西伯利亚落叶松(Larix sibirica)为主, 其次为西伯利亚云杉(Picea obovata), 其他树种还有西伯利亚冷杉(Abies sibirica)、西伯利亚红松(Pinus sibirica)、疣枝桦(Betula pendula)、欧洲山杨(Populus tremula)等。其中, 西伯利亚落叶松是乔木针叶树种, 为阳性树, 具有寿命长、抗寒、抗贫瘠等特征。

1.2 数据获取及预处理

本文主要数据包括森林资源清查数据、高程数据和气象数据。森林资源清查数据包括新疆森林资源第六次(2001年)、第七次(2006年)、第八次(2011年)和第九次(2016年)连续清查数据。该数据由样地数据和样木数据组成。样地全部为地面固定样地, 按系统抽样布设于比例尺1:50000地形图的公里网交点上, 形状为正方形, 面积0.08hm²(样地边长28.28m), 单个样地控制面积约1198.42hm²。样地数据主要包含样地坐标, 高程, 郁闭度, 坡度, 坡位, 坡向, 经纬度, 林地更新情况和病虫害等级等信息。样木数据主要包含样地号, 样木号, 样木胸径。样木检尺起测胸径为5cm, 一律用钢围尺, 读数记到0.1cm, 检尺位置为树干距上坡根颈1.3m高度(长度)处, 胸径测量误差小于1.5%。本文选择阿尔泰山落叶松林分布区内的325个样地数据及各样地内样木数据。高程数据为SRTM_30m分辨率的高程数据, 数据源于地理空间数据云(http://www.gscloud.cn/)。气象数据包括降水和气温数据, 来源于中国气象数据网(http://data.cma.cn/site/index.html)的中国地面气候资料日值数据集, 使用Excel对2001、2006、2011和2016年的新疆各气象站点降水、气温数据进行统计整理后利用ANUSPLIN软件以高程数据作为协变量进行插值, 可得新疆相应年份连续分布的降水和气温栅格数据, 在此基础上利用ArcGIS软件提取各样地点对应年份的降水和气温数值。

1.3 研究方法 1.3.1 碳储量与碳密度估算方法

阿尔泰山落叶松地上生物量和地下生物量估算模型选用高亚琪等^[24]、朱雅丽等^[17]人的研究成果, 计算公式如下：

(1)

(2)

式中, D表示胸径, 单位cm；B_t表示单株地上总生物量, 单位kg。(2)式中, B_r表示根生物量, 单位kg。

(3)

式中, C表示单株含碳量, 其单位均为kg/株, 其地上部分含碳系数为0.4759, 其根部含碳系数为0.4768^[12]。用式(3)计算得到单株样木含碳量, 对各样地所有的样木含碳量分别求和, 得到各样地总碳储量, 除以1000将单位kg转化为t, 再除以样地面积(0.08hm²), 即得到各样地单位面积落叶松碳储量, 即碳密度, 其单位为t/hm²。

1.3.2 生产力(NPP)计算方法

NPP通过森林群落年生长量P和年凋落量L来估算^[25-26], 计算公式如下:

(4)

(5)

(6)

式中: B_t为地上生物量; A为年龄; P为森林群落年生长量；L为年凋落量。

1.3.3 影响因子分析方法

本文采用标准系数法和因子分析法^{[6, 23]}评估阿尔泰山落叶松林碳储量及NPP时空变化特征的气候成因, 包括3个步骤：(1)各要素数据的均值化处理；(2)多元回归建模；(3)影响因子系数计算。

均值化处理的目标是为了消除碳储量、NPP、降水、气温等要素因单位或权重的不同所产生的影响, 进而使得各要素间对比性更强, 其处理过程如下：

(7)

(8)

(9)

(10)

式中, CI_i、NPPI_i、PI_i、TI_i分别为均值化后的第i个样地的碳储量、NPP、年降水和年均温, C_i 、NPP_i、P_i 、T_i分别为第i个样地的实际碳储量(10⁴t)、NPP(t hm^-2 a^-1)、年降水(mm)和年均温(℃)； 分别为各个样地的平均碳储量、平均NPP、平均年降水和平均年均温。

利用SPSS的多元回归分析功能和均值化后的各要素数据, 以碳储量和NPP为因变量, 年均气温和降水为自变量, 进而进行回归建模, 其基本模型为：

(11)

公式(11)中, y为因变量, x₁、x₂为自变量, a、b、c为回归方程系数。根据这些系数的大小, 可初步判断参变量对森林碳密度影响的重要程度, 这种方法叫标准系数法。

由于各系数的标准误差σ不同, 回归系数必须大于其相应的标准误差才有意义。为此, 对各元素的回归系数除以其相应的标准误差, 进而计算影响因子系数, 计算过程如下：

(12)

公式(12)中, E为影响因子系数；ST_i为受气温和降水的影响, 即当降水或气温与碳储量或NPP呈正相关时, ST=1, 反之则ST=-1；b_i为方程中不同要素的回归系数; e_i为各回归系数的标准误差值。只有当E>1时, b_i所对应参变量对因变量的回归才有意义, 且E越大, b_i所对应参变量对因变量贡献越大。因此, 可根据E的大小判断参变量对森林碳储量影响程度, 这种方法叫因子分析法。

1.3.4 方差分析法

方差分析又叫变异数分析或F检验, 反映样本在不同水平下是否具有显著性差异, 如果具有显著差异说明样本在该水平下对因变量存在显著影响。本文利用单因素方差分析和双因素方差分析分析不同因素对阿尔泰山落叶松林固碳的影响, 其显著性P>0.05表示无显著差异, < 0.05表示存在显著差异, < 0.01表示存在极显著差异, 分析工具使用SPSS软件。

2 结果与分析 2.1 阿尔泰山落叶松林固碳变化趋势

2001—2016年间, 阿尔泰山落叶松林碳储量、碳密度及NPP呈增长趋势, 且增长幅度不大(表 1)。落叶松林碳储量由2001年时的2247.45×10⁴t增加到2016年时的2345.20×10⁴t, 净增加97.75×10⁴t, 平均每年增加量为6.52×10⁴t·a^-1, 年均增长率为0.29%, 2001—2006、2006—2011、2011—2016年净增长量分别为5.79×10⁴t、37.43×10⁴t、54.53×10⁴t, 年均增长率分别为0.05%、0.33%和0.48%, 可以看出2001—2006年的碳汇作用弱于其他2个时期；碳密度由2001年的57.70t/hm²增加至2016年的60.21t/hm², 年均增长率为0.29%；2001—2006年阿尔泰山落叶松林生产力略有下降, 年均下降率为0.08%, 下降幅度不明显, 2006—2016年则呈逐年增长趋势, 年均增长率分别为0.16%和0.57%。利用SPSS软件, 不同调查期为分组变量, 以阿尔泰山落叶松林碳储量、碳密度及NPP为因变量分别进行方差分析发现：阿尔泰山落叶松林在不同调查期内的碳储量、碳密度及NPP在4个调查期内均无显著差异(P>0.05), 说明阿尔泰山落叶松林的碳储量、碳密度及NPP在调查期内变化不大。

2.2 不同林区固碳差异分析

阿尔泰山落叶松林碳储量呈“富蕴>阿勒泰>哈巴河、青河>布尔津>福海”的空间分布特征, 碳密度和NPP在青河林区最高, 在布尔津林区最低, 总体表现出“东南高, 西北低”的分布特征(表 2)。富蕴林区是阿勒泰山落叶松林最大的碳库, 在调查期内碳储量占比在38%以上；福海林区固碳量最低, 其碳储量占比小于8%。碳密度和NPP方面, 青河林区固碳能力最强, 2001—2016年平均碳密度和NPP分别为87.92t/hm²和9.47t hm^-2 a^-1；布尔津最低, 分别为44.50t/hm²和5.82t hm^-2 a^-1。利用SPSS软件分别针对4个调查期, 以不同林区为分组因子, 以碳储量、碳密度和NPP为因变量作方差分析, 结果显示2001、2006、2011、2016年阿尔泰山落叶松林碳储量、碳密度和NPP在不同林区分布均表现出极显著差异(P < 0.01), 说明阿尔泰山落叶松林固碳空间差异较大；分别针对富蕴、阿勒泰、哈巴河、青河、布尔津、福海6个林区, 以调查年份为分组因子, 以碳储量、碳密度和NPP为因变量作方差分析, 结果显示阿尔泰山落叶松林各林区碳储量、碳密度和NPP在4个调查期内均无显著差异(P>0.05), 说明阿尔泰山各林区固碳在时间尺度上均无显著变化。

2.3 不同龄组固碳差异分析

对比分析各龄组的固碳差异发现：2001—2016年, 阿尔泰山落叶松林的碳储量及碳密度在各龄组均呈现出“过熟林>成熟林>近熟林>中龄林>幼龄林”的龄组分布特征, NPP则呈现“过熟林>近熟林>成熟林>中龄林>幼龄林”的特征(表 3)。过熟林的固碳量和固碳能力均明显强于其他龄组, 在各龄组落叶松林间占据绝对优势, 表现为调查期内过熟林碳储量百分比达到49—53%, 碳密度和NPP明显高于其他龄组, 分别在78t/hm²和8.3t hm^-2 a^-1以上；幼龄林的固碳最差, 其碳储量、碳密度及NPP均远远低于其他龄组, 在阿尔泰山落叶松林固碳中作用微小；近熟林和成熟林在固碳量和固碳能力方面则表现出一定的差异, 近熟林的碳储量低于成熟林, NPP却高于成熟林, 碳密度在2001年低于成熟林, 至2016年时已超过成熟林。利用SPSS软件分别针对4个调查期, 以龄组为分组因子, 以碳储量、碳密度和NPP为因变量作方差分析, 结果显示2001、2006、2011、2016年阿尔泰山落叶松林碳储量、碳密度和NPP在不同龄组分布均表现出极显著差异(P < 0.01), 说明阿尔泰山落叶松林固碳龄组分布差异较大。时间尺度上, 各龄组碳储量均呈增加趋势, 幼龄林和中龄林碳密度及NPP呈下降趋势, 近熟林、成熟林、过熟林的碳密度及NPP呈增长趋势(表 3)。利用SPSS软件分别针对各龄组, 以调查年份为分组因子, 以碳储量、碳密度和NPP为因变量作方差分析, 结果显示：近熟林在4个调查期内碳储量、碳密度及NPP呈现出显著差异(P < 0.05), 其他各龄组碳储量、碳密度和NPP在4各调查期内均无显著差异(P>0.05), 说明近熟林固碳在2001—2016年有显著变化, 其他各龄组固碳均变化不大。

2.4 不同地形条件固碳差异分析

阿尔泰山落叶松林碳储量、碳密度及NPP随海拔升高呈现先增后降的变化特征, 在北坡和西北坡相对较高(表 4)。在不同海拔梯度上, 阿尔泰山落叶松林碳储量分布在2001—2300m最大, 其次为1701—2000m, 在1401m以下和2600m以上最小；碳密度在2001—2600m最大, 其次为1701—2000m, 在1700m以下和2600m以上相对较小；NPP在2001—2300m最大, 其次为1701—2000m和2301—2600m, 在1400m及以下最小。为探讨海拔、坡向对阿尔泰山落叶松林碳储量、碳密度及NPP的影响及其显著性水平, 分别针对4个调查期, 采用SPSS软件以海拔、坡向作为分组因子, 以碳储量、碳密度及NPP作为因变量进行双因素方差分析, 结果显示：海拔对碳储量、碳密度及NPP的影响极显著(P < 0.01), 坡向对碳储量、碳密度及NPP无显著影响(P>0.05), 海拔与坡向的交互作用对碳储量、碳密度无显著影响, 对NPP的影响在2001年和2006年显著(P < 0.05), 在2011年、2016年极显著(P < 0.01)。

2.5 阿尔泰山落叶松林固碳气候成因分析

利用标准系数法和因子分析法, 采用SPSS软件对均值化的四个调查年份的阿尔泰山落叶松林碳密度、NPP与年降水、气温数据进行回归分析并计算其影响因子系数, 其回归方程均通过了 α=0.01的显著性检验, 结果见图 2。在碳密度和NPP与降水、气温的回归方程中, 气温因子的回归系数均大于零, 降水因子的回归系数均小于零(图 2), 依据标准系数法可认为阿尔泰山落叶松林的碳密度及NPP与气温呈正相关, 与降水呈负相关, 即阿尔泰山落叶松林的碳密度随气温的上升而增加, 但随降水的增加而减少。根据影响因子分析法, 气温对碳密度的影响因子为10.967, 降水对碳密度的影响因子为3.867(图 2a)；气温对NPP的影响因子为11.570, 降水对NPP的影响因子为4.066(图 2b)。由此可见：气温对阿尔泰山落叶松林碳密度及NPP的贡献远远大于降水, 说明气温是影响阿尔泰山碳密度及NPP变化的主导气候因素。本研究结果与丁晓娟等^[22]认为的降水是影响落叶松生长的主要限制因子有所不同, 其原因可能是：(1)丁晓娟等人以1991—2011年间近20年为研究时段, 本文将2001—2016年作为研究周期, 研究时段的不同可能会造成研究结果的差异；(2)丁晓娟等着重研究了阿尔泰山南坡西伯利亚落叶松生长量与气候因子的关系, 本文主要关注于阿尔泰山落叶松林碳密度及NPP与气候的关系, 二者研究区域及研究对象的差异也可能造成研究结果的不同。

2.6 碳密度及NPP估算结果的对比

由于本研究与前人研究在研究区、森林类型、分布面积等方面存在较大的差异, 因此本研究不对碳储量作详细对比, 仅针对更具普遍对比性的碳密度和NPP的研究结果与前人研究成果进行对比(表 5)。从对比结果看, 郑拴丽、方岳等^{[19, 27]}分别以阿尔泰山落叶松和喀纳斯自然保护区落叶松为研究对象, 估算碳密度结果为64.43t/hm²和69.18t/hm²；包艳丽、许文强等^[28-29]分别研究了天山地区的云杉林和雪岭云杉, 其碳密度估算结果分别为59.55t/hm²和92.57t/hm²；王蒙等^[30]、王云霓等^[31]研究了黑龙江省落叶松人工林和华北落叶松的碳密度, 结果为30.38t/hm²和27.95t/hm²。本文估算所得阿尔泰山落叶松林在2000—2016年的平均碳密度为58.64t/hm², 与郑拴丽、方岳等的估算差异保持在±20%以内, 略低于二人的估算结果, 主要是因为郑拴丽、方岳等的碳密度估算结果不仅包括乔木层, 还包括草本层和枯落层, 本文估算所得碳密度仅为乔木层；与包艳丽、许文强在天山地区的估算结果对比具有一定差异, 主要是由于森林类型不同所导致；与王蒙、王云霓等在黑龙江和华北落叶松的估算也有一定差异, 这可能是由于所处地区不同。NPP方面, 王斌、余超等^[25-26]估算1973—2008年的全国落叶松林的平均NPP分别为7.40t hm^-2 a^-1和7.63t hm^-2 a^-1, 本研究在四个调查期的平均NPP为7.41t hm^-2 a^-1, 与王斌、余超等的研究结果差异较小。

3 讨论

阿尔泰山落叶松林固碳的时空变化特征是封山育林政策、气候条件、地形地貌、落叶松生长习性等多重因素共同影响的结果。封山育林能够显著改善林区内包括落叶松林在内的各林分类型立地条件和生态环境, 促进林区植被的生长和恢复, 进而提高林区内林分蓄积量^[32]。近年来, 阿尔泰山林区封山育林工程愈加深入, 落叶松林碳储量、碳密度及NPP呈增长趋势。同时长期的封山育林, 使得各林区内包括落叶松林在内的乔木主要以自然更新为主导, 阿尔泰山落叶松林相较于其他乔木树种处于弱势地位, 天然更新困难^[18], 进而造成成熟林和过熟林分布比重过高, 而幼龄林分布较少的现状, 加之落叶松自身寿命长, 生长缓慢的特点, 因此在调查期内落叶松林碳储量、碳密度及NPP增长幅度相对较小。水热分布状况的变化会引起树木生长的差异^[12]。来自于大西洋、北冰洋的冷湿气流沿西北方向进入阿尔泰山, 由于山体西北部高而宽, 可保留大量的冷湿气流, 进而形成了山体西北部寒冷湿润的气候特征；冷湿气流经西北部山体削弱后到达山体东南部, 由于海拔的降低及西伯利亚蒙古高压干燥反气旋的影响, 形成了山体东南部相对干燥温暖的气候特征。由2.5节知, 气温对碳密度及NPP的影响因子远大于降水, 说明气温是影响阿尔泰山落叶松林固碳的主导气候因素。这主要是由于在水分充足的条件下, 降水对树木生长的限制作用逐渐减弱, 而温度对树木生长的影响逐渐显著, 生长季内温度的升高有利于增加净光合产量^[33]。阿尔泰山整体降水相对充足, 气温升高有利于落叶松的径向生长^[34], 因此气温成为影响落叶松生长的主导因素。由于山体东南部气温较高, 因而更适宜落叶松生长, 故其碳储量、碳密度及NPP较大。阿尔泰山落叶松固碳呈现随海拔升高先升后降的特征, 这主要是由于海拔升高引起落叶松生长水热条件变化所致。阿尔泰山低海拔地区多分布于山体东南部, 该地区气温较高, 降水相对较少, 水分成为限制落叶松生长的主要限制因素。随着海拔的升高, 降水逐渐增多, 可用于落叶松生长的水分不断增加, 因此落叶松林固碳表现出随海拔升高而增加的特征；高海拔地区多位于山体西北部, 该区域降水充足, 气温成为限制落叶松生长的主要限制因素, 随着海拔的继续升高, 气温不断下降, 抑制落叶松生长, 因此落叶松林固碳表现出随海拔升高而降低的特征。在坡向分布上, 阿尔泰山落叶松林碳储量、碳密度及NPP在北坡和西北坡相对较高, 这可能与阿尔泰山“西北-东南”走向有关。本文的研究结果有效揭示了阿尔泰山落叶松林固碳时空变化规律, 为了解和掌握该区域碳汇能力及今后科学合理的森林经营政策的制定提供了有效的数据和决策支持。

4 结论

本文从时间尺度、水平格局和垂直格局出发, 以碳储量、碳密度及NPP为主要指标, 总结了阿尔泰山落叶松林2001—2016年间固碳的时空特征, 得出以下结论：

(1) 从时间尺度变化看, 阿尔泰山落叶松林固碳量和固碳能力均呈增长趋势, 但增长幅度较小。2001—2016年间, 阿尔泰山落叶松林碳储量2247.45×10⁴t增加至2345.20×10⁴t, 年均增长率为0.29%；碳密度由2001年的57.70t/hm²增加至2016年的60.21t/hm², 年均增长率为0.29%；NPP由7.35t hm^-2 a^-1增长至7.41t hm^-2 a^-1, 年均增长率为0.05%。

(2) 从空间水平分布看, 阿尔泰山落叶松林固碳呈“东南高, 西北低”的分布特征。碳储量呈“富蕴>阿勒泰>哈巴河、青河>布尔津>福海”的空间分布特征, 碳密度和NPP在青河林区最高, 在布尔津林区最低, 总体表现出“东南高, 西北低”的分布特征。从垂直分布看, 阿尔泰山落叶松林碳储量、碳密度及NPP随海拔升高呈现先增后降的变化特征, 在北坡和西北坡相对较高。

(3) 利用标准系数法和影响因子分析法对阿尔泰山落叶松林固碳的气候成因分析表明, 气温因子与碳密度及NPP的回归系数均大于零, 其影响因子远大于降水, 说明气温的升高有利于阿尔泰山碳储量、碳密度及NPP的增加, 是影响阿尔泰山落叶松林固碳的主导气候因素；降水因子与碳密度及NPP的回归系数均小于零, 说明降水的增加不利于阿尔泰山碳储量、碳密度及NPP的提高, 是影响阿尔泰山落叶松林固碳的次要气候因素。