基于林分生长模型的天山云杉碳汇潜力估测

张景路, 张绘芳, 高 健, 朱雅丽, 地力夏提·包尔汉

(新疆林业科学院现代林业研究所,新疆 乌鲁木齐 830000)

森林是全球陆地生态系统中最大也是最重要的植被碳库和碳汇[1-2]。自世界各地频繁出现气候异常现象[3]、林业增汇减排纳入国家应对气候变化规划[4]以来,森林便以其巨大的碳贮能力和碳汇潜力成为被关注的焦点。准确估算森林生态系统固碳潜力,有利于理解森林碳汇的动态过程,也是加强森林生态系统科学管理和制定林业增汇减排具体措施的前提。

目前在大的空间尺度上,碳汇潜力的预测一般采用模型模拟法和遥感法,如汲玉河等[5]采用BIOME4模型法模拟了现存森林分布区的潜在顶级森林生态系统的碳储量,用顶级森林生态系统碳储量与现存森林生态系统碳储量之差作为森林生态系统的固碳潜力;周瑞伍等[6]根据CART模型模拟得到云南省不同森林植被类型的适宜生境面积,以适宜生境的森林植被碳储量与现实森林植被碳储量的差值作为森林植被固碳潜力,并据此测算了云南省森林植被总固碳潜力;李秀娟等[7]基于生态系统观测数据和遥感数据等建立生态系统碳周转模型(TECO-R),并模拟了1982-1999年森林生态系统的碳汇量及其年际变化,估算了森林植被与非植被的固碳潜力。在区域小尺度上,碳汇潜力预测的方法有时间序列法、时空替代法等。时间序列法是利用长期调查或观测资料,得到生态系统固碳量的时间序列,进而推算生态系统的固碳速率和固碳潜力,即通过2个观测点的碳储量的差来推演森林生态系统碳蓄积的动态过程的一种方法,如徐冰等[8]以森林资源清查资料和林业发展规划为依据,通过生物量密度与林龄的关系预测了中国森林2000-2050年的碳汇潜力;李斌等[9]以同龄组林分单位面积最高年净固碳量为标准,估算现有不同龄组林分通过龄组转移达到成熟林时所固定的碳量,即未来通过自然生长增加的部分,并据此计算了湖南省杉木林碳汇潜力;陈科宇等[10]在假定森林面积不变的前提下,根据森林演替理论把各林型成熟林碳密度的平均值作为该林型演替终点的碳密度,并作为该林型理论最大固碳潜力的参照系估算了研究区单位面积固碳潜力。在这些基于时间序列的研究中,尤其是对于慢生树种,往往存在着龄组划分间隔过大的问题,动辄10～20年龄组的龄级期限使得研究结果往往呈现阶梯状分布,预测局限性大。此外,在不同环境条件下的林分往往具有不同的生长发育水平和碳汇潜力,这些不同在以往学者的研究中很少加以区分。

本研究依据林分调查因子计算天山云杉地位级指数、林分密度等多种指标,以4种生长曲线为基础构建了林分蓄积生长模型,并通过生物量蓄积量回归模型和含碳系数最终构建林分碳储量生长模型。设定天山云杉林分极限年龄时的碳储量作为林分最大碳容量,将林分最大碳容量与当前时点的碳储量的差值作为该林分的碳汇潜力。此外,将特定目标年和特定环境条件下天山云杉碳密度的年均增长量作为碳汇潜力的另一指标,通过林分碳储量生长模型计算不同林分条件下天山云杉由当前生长到180 a碳密度年均增长量,以定量评价不同立地条件和林分密度下的研究区天山云杉林分碳汇潜力,并预测30 a、60 a后的天山云杉林分碳储量。

1 材料与方法

1.1 研究区概况

研究区(图1)为天山西部国有林管理局巩留分局东部山区和西天山国家级自然保护区[11],位于新疆伊犁州巩留县,地处伊犁河谷中部,天山支脉那拉提山北麓,地理坐标为81°34′E～83°35′E,42°54′N～43°38′N。气候属温带大陆性半干旱气候,年平均温度7.4 ℃,年日照时数2 731.7 h,年降水量200～780 mm[12];研究区东部、南部和北部有海拔3 500～4 200 m的那拉提山系和伊什格力克山系,阻隔了蒙古高压干燥气流、南下的西伯利亚寒流和新疆南部干热风,地势向西倾斜而开阔,使得来自大西洋、里海的暖湿气流进入且受高山阻挡形成丰富降水。区域内水系丰富,有发源那拉提山系的乌勒肯库尔德宁、沙特布拉克、协天德3条主要河流,形成基本平行且南北走向的3条河谷。林区土壤主要为半淋溶土纲的灰褐土和钙层土纲的黑钙土、栗钙土[13]。研究区景观格局具有中亚山地森林特征[14],林地和草地交错分布,林型简单,树种单纯,保存着完整的原始云杉类型和山地植被自然垂直带谱,是天山林区生长最好、分布最为集中的原始云杉林,也是天山山地森林生态系统的典型代表。

图1 研究区地理位置及样点分布示意

1.2 数据来源与处理

采用研究区森林资源二类调查数据,提取研究区天山云杉相对纯林小班,以小班的平均年龄、平均树高、平均胸径和林分密度为基础数据。以10为步长划定各林阶,以正负2倍标准差为依据剔除各林阶范围内林分平均树高、胸径异常的样本,最终得到133个天山云杉样本数据。各林阶平均年龄、平均树高、平均胸径、林分密度详见表1。

表1 样本描述性统计

根据样本的极值将样本数据分组,以频率和组距的商为高、组距为底绘制频率直方图和正态曲线,计算样本的偏度、峰度并进行柯尔莫戈洛夫-斯米诺正太分布检验(K-S检验)。结果(表2)显示,平均年龄、平均树高与平均胸径偏度小于0,呈负偏态分布,林分密度呈正偏态分布状态;4组数据峰度均小于0,即表征数据分布比正态分布的峰度高且狭窄。样本数据的偏度与峰度绝对值均小于1。4个样本的K-S检验显著性均大于0.05,说明样本服从正态分布,即在0.05水平上符合正态分布。

1.3 研究方法

1.3.1 地位级导向曲线方程 基于各样本平均树高和平均年龄数据,拟合Gompertz、Logistic、Mitscherlich和Schumacher等4个常用树高生长曲线方程(表3);比较各方程决定系数(R2)、残差平方和(Q),选择最符合树木生长规律和拟合效果最佳的方程作为地位级导向曲线。

表3 地位级导向曲线方程

1.3.2 地位级指数(SCI) 采用林分平均树高及林分平均年龄为自变量建立地位级指数方程,拟合得到样本地位级指数作为立地指标。

式中,H为林分平均树高;age为林分平均年龄;age0为林分基准年龄;b、c为参数。

1.3.3 林分密度指数(SDI) 林分密度指数[15]是林分在标准平均胸径时所具有的单位面积的株数,即林木对其所占空间的利用程度,是森林生长的重要影响因子[16]。这一指标能综合反映平均胸径和株数的情况,且与立地条件和林龄不共线[17]。计算公式为:

SDI=N×(Dg/D0)α

式中,N为林分密度(株,1 hm2),Dg为林分现实平均胸径,D0为林分标准平均胸径,α为自然稀疏率。

1.3.4 林分生长模型

1.3.4.1 林分平均直径生长模型 对于生长未受干扰的天然林来说,结合林龄、地位级指数和林分密度指数3个自变量即可得到拟合效果良好的生长模型,实现天然林林分生长的预测[18]。将林龄、地位级指数和林分密度指数指标分别引入Gompertz、Logistic、Mitscherlich和Schumacher等4个常用生长曲线方程(表4),得到林分平均直径生长模型。

表4 林分平均直径生长模型

1.3.4.2 林分蓄积量生长模型 将林分平均年龄、地位级指数和林分平均直径作为基础变量引入Gompertz、Logistic、Mitscherlich和Schumacher等方程中,即可求得林分蓄积量生长模型(表5)。

表5 林分蓄积量生长模型

1.3.4.3 林分碳储量生长模型 根据张煜星等[19]模拟的云杉单位生物量与单位蓄积量回归模型,可将林分蓄积量转换为林分生物量,模型如下:

Bα=0.62VOL0.984 4

式中,Bα为林分地上生物量;VOL为林分蓄积量。

根据张绘芳等[20]、朱雅丽等[21]等研究的云杉含碳系数,可将林分生物量转换为林分碳储量,即可得到林分碳储量生长模型:

CS=0.462 4Bα×R

式中:CS为林分碳储量;R为林分含碳系数。

1.4 精度验证

为了比较不同地位级导向曲线方程所建立的模型的预测精度,选取了5个评价指标(表6),分别为决定系数(R2)、估计值的标准差(SEE)、总系统误差(TRE)、平均系统误差(MSE)、模型精度(P),从不同维度对模型的预测精度进行检验。

表6 模型精度验证方法

1.5 碳汇潜力计算

当前关于中国森林碳汇潜力的研究较少,评估方法也处于探索阶段,不同学者对碳汇潜力的概念有着不同的定义。本研究将碳汇潜力定义为某一森林系统根据森林演替理论演替为顶级自然林时所达到的最大碳储量与现实碳储量之间的差值,计算方法为:

CSC=CSCmax-CSCt

式中,CSC为潜在碳储量;CSCmax为生态系统最大碳储量;CSCt为现实碳储量。

2 结果与分析

2.1 林分因子计算结果

2.1.1 地位级导向曲线 由表7可知,各方程拟合精度差异不大,本研究综合分析各方程拟合曲线发现,Gompertz方程是具有初始值的典型“S”形单调递增函数,拐点约位于最大值的三分之一处,在生物学工程中适用范围较广[22];Logistic方程是生物种群动态模拟中最常用的方程,具有对称型“S”形生长曲线,拐点约位于最大值的二分之一处,适合用于描述慢生树种的生长趋势[23];Mitscherlich方程是单分子方程,曲线无拐点,形态接近肩形,适合初期生长较快、无拐点树种的生长过程[24];Schumcher是Korf方程的特定形式,是基于同龄纯林树木年均增长量的同比增长率与年龄呈反向相关的假设原理得出的方程,其曲线存在1个拐点[25]。综上所述,结合天山云杉的生长特点,选取Logistic方程作为地位级指数导向曲线最终方程,即:

表7 地位级导向曲线拟合参数及精度

Y=28.014/[1+8.774exp(-0.031age)]

2.1.2 地位级指数 基准年龄是指树高生长趋于稳定且能灵敏反映立地差异的年龄[26]。确定基准年龄通常考虑树种生活史平均年龄(主伐年龄)一半左右的年数,或树高、林分蓄积量生长过程中连年生长量或平均生长量最大的时期[27]。取天山云杉主伐年龄160 a的一半(80 a)作为基准年龄,将表7中Logistic方程的参数b和c带入地位级指数方程得到:

SCI=H×{[1-exp(-0.031×80)8.774]/[1-exp(-0.031×age)8.774]}

式中,SCI为地位级指数,H为林分平均树高,age为林分平均年龄。

2.1.3 林分密度指数 选取天山云杉基准年龄的平均胸径作为标准平均胸径,即D0=24 cm,方程斜率(α)受不同树种影响很小,因此借鉴方怀龙[17]的成果取-1.605。最终得到林分密度指数计算公式:

SDI=N×(Dg/24)-1.605

式中,SDI为林分密度指数,N为林分密度,Dg为林分平均胸径。

2.2 林分生长模型

2.2.1 林分平均直径生长模型 采用SPSS非参数回归模块,以地位级指数、林分密度指数和林龄为自变量,林分平均胸径为因变量,计算林分平均直径生长模型参数并验证模型精度。由表8可知,Gompertz方程和Schumacher方程估计值的标准差、总系统误差、平均系统误差的绝对值均显著低于Logistic方程和Mitscherlich方程,且前两者模型精度(P)均高于85.000%,其中,Gompertz方程的总相对误差较低而模型精度较高,因此选取Gompertz方程作为天山云杉平均直径生长模型,表达式为:

表8 林分平均直径生长模型参数拟合及精度验证

DBH=-1.095×SCI0.373×(-5.154×SDI-0.082)×exp(0.004×age)

式中,DBH为林分平均直径,SCI为地位级指数,SDI为林分密度指数,age为林分平均年龄。

2.2.2 林分蓄积量生长模型 采用SPSS非参数回归模块,以地位级指数、林分平均直径和林龄为自变量,以林分每公顷蓄积量为因变量,计算林分蓄积量生长模型参数并验证模型精度。由表9可知,各生长曲线间误差和精度差异均不大,其中Logistic方程模型精度最高,Schumacher方程的误差的绝对值最小。进一步分析4个方程的差异,可知Logistic、Gompertz和Mitscherlich方程和林分年龄相关的参数均过小,使得生长曲线不能反映林分蓄积量随林龄增长的变化趋势,导致曲线存在初始值过大而增长量过小的问题,而Schumacher方程曲线形态更符合林分生长过程,因此选取Schumacher方程作为天山云杉林分蓄积量生长模型,表达式为:

表9 林分蓄积量生长模型参数拟合及精度验证

VOL=0.668×SCI2.101×0.652×DBH-0.018×exp(-3.124/age)

式中,VOL为林分蓄积量,SCI为地位级指数,DBH为林分平均直径,age为林分平均年龄。

2.2.3 林分碳储量生长模型 根据张煜星等[19]、张绘芳等[20]、朱雅丽等[21]的研究成果,得到云杉林分碳储量生长模型:

简化后得到:

式中,CS为林分碳储量,SCI为地位级指数,DBH为林分平均直径,age为林分平均年龄。

2.3 云杉碳汇潜力计算

以研究区森林资源二类调查数据为基础,提取研究区的乔木林地和疏林地小班数据并计算林分平均年龄、平均树高、平均胸径和林分密度等林分因子。由于地位级指数反映的是林木生长环境的综合条件[28],所以在可预见的时间范围内不会产生大的变动。林分密度指数受林分自然更新和枯损影响。提取和分析1996-2016年20年间研究区30块固定复测样地数据,研究发现,林分的更新和枯损受自然和人为的各类因素影响,具有一定程度的随机性,与林分、林阶、树高、胸径等均无显著相关性。综上,保持立地指数和林分密度指数不变,以林分平均年龄为基础变量对研究区天山云杉碳汇潜力进行计算。

2.3.1 林分碳汇增长量 以新疆森林资源补充调查操作细则中规定的天山云杉过熟林下限160 a加天山云杉龄级期限20 a后(180 a)作为林分平均年龄极限。将研究区天山云杉的地位级指数以4为跨度分为8级,再将林分密度指数以400为跨度分为5级,扣除9个无现实分布的空值后得到阶梯式分布的不同地位级指数和林分密度的林分类型共31个。基于林分碳储量生长模型,取每一级别的中值计算不同立地条件和林分密度下0～180 a林分碳密度年均增长量,结果见表10。由表10可知,林分碳密度年均增长量主要受立地条件影响,在相同林分密度条件下,随着林分立地条件的变化,林分碳密度0～180 a年平均增长量为0.020～0.641 t/(hm2·a);林分密度对林分碳密度年均增长量的影响较小,随着林分密度的增加,林分碳密度年均增长量略有提升,在立地条件较差的区域林分密度带来的影响可忽略不计。采用各林分面积结合各林分类型碳汇平均增长量加权平均求得研究区全域天山云杉林分碳密度年平均增长量为0.299 t/(hm2·a),年均增长量拐点位于30 a处。

表10 不同林分条件下碳密度年均增长量

2.3.2 林分未来时刻的碳汇预测 由于研究区天山云杉林分年龄整体较大,林龄超过120 a的占比达到80.00%以上,即60 a后研究区大部分林分平均年龄将达到过熟林上限,因此选择预测30 a后和60 a后2个时间节点的碳储量。以实际调查时间为起始时间节点,采用林分碳储量生长模型对研究区天山云杉林分碳储量当前值、未来30 a和未来60 a的林分碳储量进行预测,得到结果分别为3.439×106t碳、3.447×106t碳、3.450×106t碳,未来30 a、未来30～60 a增长量分别为8×103t碳、3×103t碳,涨幅分别为0.233%和0.087%。

将林分平均年龄达到过熟林上限(180 a)时的森林碳储量作为森林最大碳容量,将森林最大碳容量和林分当前时点森林碳容量之间的差值作为碳汇潜力,根据林分碳储量生长模型求得研究区天山云杉林分最大碳容量为3.451×106t碳,碳汇潜力为1.245×104t碳。

3 结论与讨论

本研究分别基于Gompertz、Logistic、Mitscherlich和Schumacher等4个常用生长曲线方程,采用林龄、平均树高、平均胸径和林分密度等构建林分蓄积量生长模型,通过生物量蓄积量回归模型和含碳系数建立林分碳储量生长模型。其中Schumacher方程曲线形态更贴合天山云杉生长状况,所构建模型的精度较高,误差绝对值最小,模型精度89.082%,估计值的标准差13.006,总系统误差-0.293,平均系统误差-5.943,决定系数0.895,表明所构建的林分碳储量生长模型具有较高的精度和稳定性,可用于研究区林分碳汇潜力的估测。

以林分碳储量生长模型为基础计算出研究区天山云杉林分碳密度0～180 a年平均增长量为0.020～0.641 t/(hm2·a),研究区全域天山云杉林分碳密度年平均增长量为0.299 t/(hm2·a),年均增长量拐点位于30 a处。研究区全域天山云杉碳汇潜力为1.245×104t;当前、未来30 a和未来60 a林分碳储量分别为3.439×106t碳、3.447×106t碳、3.450×106t碳,未来30 a、未来30～60 a增长量分别为8×103t碳、3×103t碳,涨幅分别为0.233%、0.087%。研究区天山云杉成熟林、过熟林占比较高,林分碳汇潜力低,需进行林龄结构优化,以促进天山云杉林的可持续发展。

本研究采用研究区天山云杉的林龄、平均树高、平均胸径和林分密度等数据,建立了天山云杉的林分平均直径生长模型和林分蓄积量生长模型,并通过天山云杉生物量和蓄积量回归模型、含碳系数的转换构建了天山云杉林分碳储量生长模型。相较于需要借助土壤特性、气候因子或生物量的其他碳潜力预测数学模型方法,该方法利用现有林地资源调查数据中的主要调查因子即可实现,在实际的林业工作中模型需要的数据具有易得性和普适性。相较于生物量清单法和样点法中普遍采用的以龄组划分作为时间间隔的模式,林分生长模型法在时间序列上具有更精细的分辨率,可以灵活地预测未来任意时刻的林分因子和碳储量并计算其碳汇潜力。

研究区天山云杉碳储量增长率和单位增长量都低于李奇等[29]和郑拴丽[30]研究中的平均水平。这主要是因为研究区天山云杉林分具有不尽合理的年龄结构,森林资源二类调查数据显示,天山云杉平均林分年龄小于60 a的幼龄林仅占总面积的7.35%,成熟林和过熟林占总面积的80.00%以上,因此,其碳储量增长率和单位增长量均较低。如果没有人工干预(择伐、更新、造林、抚育等),在自然生长的情况下,研究区天山云杉林分碳储量将很快到达最大值。

在对林分碳密度及其变化的分析中,笔者发现林分碳密度年均增长量受立地条件的影响大而受林分密度的影响小。由于生态系统的承载力是有限的,在光照、水文、土壤等条件的约束下,森林系统的生产能力、贮藏能力不可能无限制地增长,在承载力一定的情况下,林分密度如共享潜在生产力的分母,分母越大,每株单木所分享到的资源就越少,其生长态势也随之由快到慢。从本研究中对碳密度年均增长量研究的部分也可以看出,研究区林分密度越大的地方,立地指数的分布区间就越低,在立地指数最大的区域,即同树种同林龄下树高最高的区域,林分密度没有超过1 hm2400株。最终在立地条件和林分密度的双重影响下,出现了林分碳密度增长量随林分密度的增长呈现略微增加这一结果。

森林的可再生性,即森林的更替,是其一切潜力的来源。森林的碳汇潜力不仅与环境条件、森林质量有关,更受到森林面积变化的影响。现有天山云杉林之外的区域是否适宜其更新、生长?该适宜区面积、固碳潜力有多大?这些将是下一步研究的方向。