刘 帅 李建军 卿东升 朱凯文 马振燕

(中南林业科技大学 长沙 410004)

气候变化是21世纪全球面临的最严峻挑战之一(Parmesanetal.，2003)。据政府间气候变化专门委员会(IPCC，2001)预测，2100年全球平均温度将比2000年升高1.4～5.8 ℃，约是20世纪升温幅度的2～10倍。森林对温度、降水和二氧化碳浓度改变通常较为敏感，未来气候变化将不可避免地影响森林生长和收获、树种组成和分布、森林生产力及其各项功能(Battlesetal.，2008)，且这些影响会随着森林类型、生长地点以及气候场景不同而异(Oboiteetal.，2019)。因此，在气候变化背景下，深入了解气候对森林生长的影响有助于林业工作者应对未来气候变化所带来的挑战(Temperlietal.，2012)。

基于森林生长和收获模型预估未来林木生长动态，一般将竞争、立地条件和林分结构作为协变量(Bigingetal.，1995；Peng，2000；Portéetal.，2002；Bravo-Oviedoetal.，2008)，很少考虑气候因素，即假定森林生长过程中气候条件不发生改变，因此气候与林木生长之间的关系尚不明确；但一旦未来气候发生改变，将给森林经营带来很大的不确定性(Subedietal.，2011；2013)。近年来，气候与林木生长的关系开始得到越来越多关注，一些研究已将气候变量纳入到林分或单木生长建模中(Wangetal.，2007；Cortinietal.，2011；Soniaetal.，2012；Nighetal.，2012；Subedietal.，2011；2013；Sharmaetal.，2015；Aubry-Kientzetal.，2017)。Subedi等(2011；2013)研究发现，林木直径生长与气候变量之间存在复杂的非线性关系，年平均温度和降水显著影响加拿大东部北方森林中短叶松(Pinusbanksiana)和黑云杉(Piceamariana)的胸径生长。Pokharel(2009)研究指出，年平均温度是影响加拿大安大略省4种常见树种断面积生长的主要气候因子。其他类似的研究还有很多，如Huang等(2011)基于广义可加模型研究气候变化对美国北部火炬松(Pinustaeda)的影响、余黎等(2014)探讨气候因子对长白落叶松(Larixolgensis)-云杉(Piceajezoensisvar.microsperma)冷杉(Abiesnephrolepis)林6种主要树种单木胸径生长的影响及不同树种对气候变化的响应等。可见，温度和降水是影响林木直径生长的关键气候因子(Subedietal.，2011；2013；Sharmaetal.，2015；Oboiteetal.，2019)。

林木生长数据通常以单木、林分或区域进行嵌套或分组，一般具有纵向(多层)结构，且由于其获取方式依赖于每隔一段时间对同一目标进行重复测量，数据往往存在较强的相关性，因此采用传统回归分析方法如最小二乘法可能导致模型参数有偏估计。而混合效应模型通过引入随机效应描述组内和组间变异，为重复或嵌套数据提供了解决方案(Bates，1990；Cudecketal.，2002)。在森林生长建模中，混合效应模型相比传统回归模型具有更好的拟合效果和预测精度(Fangetal.，2001；Jordanetal.，2005；Yangetal.，2011；Subedietal.，2011；2013；Sharmaetal.，2015；Saudetal.，2019)。

第八次全国森林资源清查结果显示，我国共有栎类林约1 672万hm2，占全国森林总面积的10.15%(Wangetal.，2019)。青冈栎(Cyclobalanopsisglauca)是构成我国亚热带森林的主要成分之一，因具有良好的材质(如木质坚硬、耐久耐腐等)已成为我国南方地区最重要的用材树种，且在水土保持、维持生物多样性等方面发挥着重要的生态功能(梁贵等，2015)。目前，已有一些关于青冈栎生长的研究，如刘洵等(2019)评价了湖南省天然栎类林的立地质量，编制了天然栎类林胸径指数表；龙时胜等(2018)采用序类聚类方法建立了不同生长阶段的青冈栎生长方程。此外，也有关于青冈栎生物量的研究(丁增发，2014)。但这些研究均未涉及气候变化，关于青冈栎生长过程中是否受气候影响以及受气候影响的程度如何依然不明确。鉴于此，本研究首先构建气候敏感的青冈栎混合效应模型，识别显著影响青冈栎胸径生长的气候变量，在此基础上，量化和预测未来不同气候场景下青冈栎胸径生长动态，以期为未来气候变化条件下的青冈栎林经营决策提供理论依据。

1 数据与方法

1.1 数据来源

青冈栎单木生长数据来源于中南林业科技大学芦头林场(28°31′—28°38′N，113°51′—113°58′E)。该林场地处罗霄山脉北段，属大陆性季风气候，年均气温18.5 ℃，平均降水量1 450.8 mm，平均日照时数1 930 h。全场总面积5 308 hm2，地貌类型以中山、低山为主，平均坡度25°左右，土壤以红壤和山地黄壤为主。

在林场内干扰较少的青冈栎次生林中设置8块固定样地，样地大小均为20 m×20 m。在样地中选择30株33～54年生、长势良好、无病虫害的青冈栎，于2016年7月将其伐倒。采用横剖法(周生祥等，1989)解析树干，测量并记录每个龄级的胸径、树高、材积以及树龄对应的年度等，共获得320组单木生长数据。限于篇幅考虑，本研究仅利用胸径、树龄数据探讨青冈栎生长。为便于后续处理，以树号进行分组，组内按龄级排序。

气候数据来源于ClimateAP气候模型，该模型以无标度数据作为基线，根据纬度、经度和海拔计算指定位置的月、季和年度气候变量，时间跨度为1901—2100年。通过对比亚太地区1 805个气象站的观测数据，验证了气候数据的准确性(Wangetal.，2015；2017)。本研究选择与林木生长相关的常用气候变量(Subedietal.，2011；2013；Sharmaetal.，2015；Saudetal.，2019；Oboiteetal.，2019)，包括年平均温度(mean annual temperature，MAT)、最热月均温(mean warmest month temperature，MWMT)、最冷月均温(mean coldest month temperature，MCMT)、MWMT与MCMT温差(temperature difference between MWMT and MCMT，TD)、年平均降水(mean annual precipitation，MAP)、干燥指数(annual heat-moisture index，AHM)。采用ClimateAP气候模型提取树龄对应年度的气候数据。林木生长和气候数据统计如表1所示。

表1 林木生长和气候数据统计

1.2 研究方法

1.2.1 基础模型和再参数化模型 本研究通过树龄将青冈栎林木生长和气候数据汇总成一个数据集，采用最小二乘方法(least squares method)拟合Mitscherlich、Gompertz、Korf、Richards和Logistic 5个候选指数生长方程，各生长方程的评价结果见附表1。根据赤池信息准则AIC(Akaike’s information criterion)、贝叶斯信息准则BIC(Bayesian information criterion)、对数似然Loglik(Log-likelihood)、残差平方和RSS(residual sum of squares)、平均相对误差MRE(mean relative error)、均方根误差RMSE(root mean squared error)的评价结果，Mitscherlich生长方程性能最好。因此，本研究选择Mitscherlich生长方程作为基础模型。Mitscherlich生长方程及相关评价指标定义如下：

y=β0(1-e-β1t)；

(1)

(2)

(3)

(4)

式中：y为胸径；t为树龄；β0和β1分别为待估计的渐近线(asymptote)参数和生长速率(rate)参数；obji为实际测量值；esti为相应的估计值；n为样本数。

在基础模型中，林木生长完全依赖于树龄(唯一协变量)。为了量化气候变化对林木胸径生长的影响，本研究将基础模型中的参数(β0和β1)表示为气候变量的线性组合，即再参数化(Wangetal.，2007；Bravo-Oviedoetal.，2008；Sharmaetal.，2015；Zangetal.，2016)。同时，为了避免模型过拟合及多重共线性，首先采用向后逐步回归方法(backward stepwise regression)对变量进行初步筛选，然后计算剩余变量的方差膨胀因子(variance inflation factors, VIF)。将VIF最大的变量从模型中剔除，再重新拟合；重复该过程直至所有剩余变量的VIF均小于5(一般VIF低于5，可认为自变量之间无共线性)为止。

1.2.2 非线性混合效应模型 构建混合效应模型是为了处理重复和嵌套的测量数据，其一般形式如下：

yi=f(α,β)+εi。

(5)

其中，

β=λ+γ；

γ～N(0,D)；

εi～N(0,Ri)；

式中：yi为测量值向量；f(.)为协变量向量α的非线性函数；β为待估计的混合效应参数向量，通常由λ和γ2部分组成：λ为固定效应参数向量；γ为随机效应参数向量，假设随机效应γ服从正态分布；D为其方差-协方差矩阵；εi为误差项，假设服从正态分布；Ri为误差项方差-协方差矩阵，用于表示组内误差；σ2为模型误差方差值；Φi为描述组内异方差(heteroscedasticity)的对角矩阵；Γi为描述组内相关性的自相关结构。

以上计算过程均在R软件(版本3.5.2)中完成。青冈栎单木生长基础模型和再参数化模型采用最小二乘法拟合，由R包的“stats”完成；青冈栎胸径生长混合效应模型以单株样木为随机效应，其建模和预测由R包的“NLME”完成。

2 结果与分析

2.1 再参数化模型

经逐步回归和多重共线性筛选，模型保留气候变量包括年均温MAT、最热月均温MWMT和最冷月均温MCMT。因此，原基础模型中的参数β0和β1可以表示为上述气候变量的线性组合：

β0=λ00+λ01MAT+λ02MWMT+λ03MCMT；

β1=λ10+λ11MAT+λ12MWMT+λ13MCMT。

(6)

式中：λ00～λ03和λ10～λ13为待估计参数。

然而，并非每个气候变量都对青冈栎胸径生长具有显著影响，因此需对变量对应的参数进行显著性检验。每次剔除P最大的不显著变量，再重新拟合模型，重复以上过程直至所有变量均显著(P< 0.05)为止。经处理，含气候变量且参数统计显著的再参数化模型如下：

yij=β0(1-e-β1Tij)。

(7)

其中，

β0=λ00+λ01MCMTij；

β1=λ10+λ11MCMTij。

式中：yij、Tij和 MCMTij分别为第i株林木第j次测量时的胸径、树龄和气候变量。

2.2 非线性混合效应模型

yij=β0(1-e-β1Tij)+εij。

(8)

其中，β0=(λ00+γ00)+(λ01+γ01)MCMTij；

β1=(λ10+γ10)+λ11MCMTij；

εij～N(0,Ri)。

如表2所示，渐近线参数项中的MCMT系数λ01为负，表明MCMT负相关于青冈栎胸径生长，MCMT升高可能导致青冈栎胸径生长量减少；而指数项中的MCMT系数λ11为正，将加速上述过程。

表2 非线性混合效应胸径生长模型的参数估计

2.3 模型评价

在逐步完善青冈栎生长建模过程中，本研究依次构建了基础模型【式(1)】、再参数化模型【式(7)】和混合效应模型【式(8)】。为便于比较，表3列出了3个生长模型的性能评价指标。可以看出，再参数化模型由于引入气候变量，除BIC外的各项评价指标均略优于基础模型，表明再参数化模型可用于解释气候变化对林木生长的影响；但基础模型和再参数化模型均采用传统最小二乘法拟合，未考虑林木生长数据的纵向和分组特征，无法修正数据间的相关性和残差异方差，因此二者的性能包括拟合优度和误差水平均远逊于混合效应模型。与再参数化模型相比，混合效应模型的拟合优度指标AIC从1 070.342 0降至582.217 8，BIC从1 087.745 0降至620.504 9，LogLik从-530.171 1升至-280.108 9；其误差指标RSS从1 165.440 9大幅减至57.733 2，MRE由0.263 0降至0.139 1，RMSE从2.208 2显著降至0.491 5。

表3 模型拟合统计比较

经典的回归分析要求误差满足同方差性，否则可能影响模型拟合与模型检验。图1a、b所示为基础模型和再参数化模型的残差分布，可见存在明显异方差性，直观体现在随着胸径增加，残差分布逐渐呈现喇叭状开口；但在混合效应模型中，异方差性被消除(图1c)。由于混合效应模型本身控制异方差效果较好，故本研究没有在误差项中进一步指定方差函数和自相关结构。

图1 3个生长模型的残差

图2所示为部分青冈栎的胸径-树龄生长曲线。每幅子图内均有2条曲线，其中，蓝色曲线由再参数化模型(仅含固定效应)拟合而成，代表所有青冈栎在生长期内的平均胸径；红色曲线由混合效应模型拟合而成，代表单株青冈栎的胸径生长。为了对比，每幅子图中还标记有该株青冈栎胸径的实测值(散点)。可见，蓝色曲线大多偏离于实测值，而红色曲线则与实测值变化相吻合，更能反映林木个体生长情况，也体现了混合效应模型的低误差水平。2条曲线差异是由单木水平(以单木为分组依据)的随机效应决定的，通过引入单木随机效应，可使得混合效应模型无论是误差水平还是拟合优度相比基础模型和再参数化模型均有实质性提升。

图2 由固定效应和混合效应拟合的胸径-树龄生长曲线

2.4 模型预测

为了预测气候变化下的青冈栎胸径生长，本研究采用3种典型浓度路径(representative concentration pathways，RCP)作为未来气候场景(Vuurenetal.，2011)，分别为RCP2.6、RCP4.5和RCP8.5，依次代表未来温室气体低、中、高3种排放标准(2.6、4.5和8.5 W·m-2)。所有气候场景下的气候数据均通过ClimateAP 气候模型预测，由于ClimateAP气候模型以30年为时间间隔批量产生气候数据，因此将2011—2100年分为3个等间隔时间段(2011—2040年、2041—2070年和2071—2100年)，分别预测不同气候背景下青冈栎1～30年的生长过程。表4所示为3个时间段不同气候场景下的气候数据统计情况。为便于比较，当前气候条件也作为一种气候场景。

表4 不同气候场景下2011—2100年气候变量(MCMT)的预测值

采用非线性混合效应模型【式(8)】预测未来青冈栎胸径生长，结果如图3所示。图3分阶段展示了2011—2100年不同气候场景下青冈栎胸径-树龄(1～30年)生长曲线，其中胸径为所有单木预测值的平均值。2011—2040年，青冈栎胸径生长差异非常小(图3a)；2041—2070年，青冈栎胸径生长随着不同气候场景开始出现明显分化(图3b)；2071—2100年，青冈栎胸径生长差异非常大(图3c)。在树龄相同情况下，3种典型浓度路径尤其是RCP8.5场景下青冈栎胸径显著低于气候条件不变时，其生长随着RCP浓度递增依次减少。这表明，随着碳排放量增加，未来气候条件很可能不利于青冈栎胸径生长。气候对青冈栎胸径生长的影响2071—2100年相比2041—2070年更显著，当2100年树龄30年时，RCP2.6、4.5和8.5场景下青冈栎胸径相比气候条件不变时分别下降6.3%、15.6%和53.1%；而在2070年树龄30年时，RCP2.6、4.5和8.5场景下青冈栎胸径相比气候条件不变时仅下降3.7%、6.2%和21%。可见，未来随着时间推移，青冈栎胸径生长受气候因素的负面影响将更加强烈。

图3 基于MCMT平均值的不同气候场景下青冈栎胸径-树龄生长曲线

3 讨论

树木生长是一个复杂的非线性过程，通常采用指数生长方程描述(Budhathokietal.，2008；Subedietal.，2013；Trasobaresetal.，2016；Timilsinaetal.，2013)。本研究中，Mitscherlich生长方程被选为基础模型用于初步描述青冈栎胸径生长，其渐近线和生长速率参数被再参数化为气候变量的线性组合，并采用混合效应建模方法拟合胸径生长模型参数。结果表明，加入气候变量的青冈栎生长模型能够响应气候变化对林木生长的影响(Wangetal.，2007；Sharmaetal.，2015；Zangetal.，2016；Saudetal.，2019)，而随机效应引入则进一步提升了生长模型的拟合优度和预测精度(Doradoetal.，2006；Timilsinaetal.，2013；Yangetal.，2013；Zhuetal.，2019)。

气候被认为是林木生长的潜在驱动，Latta等(2010)指出，森林生长和生产力受环境条件特别是气候条件的影响。一般而言，气温升高通过产生额外热量促进林木生长(Monserudetal.，2008)，且温暖的气候条件还会延长树木生长期，从而有利于林木径向生长(Bronsonetal.，2009；Huangetal.，2010；2013)。但也有研究表明，气温升高引起的蒸散量增加导致土壤水分减少，可能会对树木生长产生负面影响，如果降雨量跟不上气温上升速度，森林生产力反而可能下降(Soniaetal.，2012)。Subedi等(2013)研究发现，气温升高可能是导致北方针叶林中黑云杉胸径生长减少的主要原因。Aubry-Kientz等(2017)指出，如果未来气候变暖持续下去，内华达山脉中一些树种的生长量可能会大幅下降。本研究亦发现，最冷月均温与青冈栎胸径生长呈显著负相关。

本研究选取3种典型浓度路径(RCP2.6、4.5和8.5)作为未来气候场景，预测2011—2100年青冈栎胸径生长。相比气候条件不变时，高排放的RCP8.5对青冈栎胸径生长的不利影响更大，而低排放的RCP2.6对青冈栎胸径生长的负面影响相对较小，且这些影响在2071—2100年尤其显著。结合模型预测结果，随着温室气体排放增加，最冷月均温对青冈栎胸径生长的负面效应还会加剧。类似地，Sharma等(2018)量化分析了气候变化对加拿大安大略地区红松(Pinusresinosa)林生产力的影响，结果发现，在RCP8.5和RCP2.6场景下红松林林分高度均显著下降，且随不同地区和气候差异而有所变化。Battles等(2008)认为，未来所有气候场景都将减缓加利福尼亚内华达山脉的针叶混交林生长，并导致针叶树种生产力显著下降。

本研究中，最冷月均温是唯一显著影响青冈栎胸径生长的气候变量，其他气候变量如年平均温度、最热月均温和年平均降水等在统计上均不显著，因此暂时无法评估这些气候变量对青冈栎胸径生长的影响；特别是降水，通常被认为是树木生长的另一类关键气候因素，当温度适宜时，充沛的降雨可以促进树木和森林生长(Wangetal.，2007；Sharmaetal.，2015；Zangetal.，2016)。本研究数据集中有2个与降水相关的气候变量，即年平均降水和干燥指数，干燥指数结合了年平均温度和年平均降水2个重要指标，非常适合描述一年中气候干燥情况(Zhangetal.，2014)，但这2个与降水相关的气候变量由于统计不显著而被排除在生长模型之外。

林分结构、土壤和地形等有关变量有助于提高气候变化影响林木生长的预测(Rubénetal.，2015；Mansoetal.，2015；Oboiteetal.，2019)。Lo等(2010)研究发现，只有气候变量还不足以解释林木胸径生长变化，且胸径生长与林分密度和相邻木之间的竞争密切相关。然而，由于缺乏相应数据，本研究生长模型中未包含与林分结构、土壤和地形相关的协变量，如果将来收集到相关数据，可以考虑将这些变量纳入到模型中。此外，由于本研究林木数据来源于同一样地，由空间差异引起的气候变化无法在生长模型中得到体现。同时，考虑到气候变化强烈影响着森林生长、树种组成和分布(Chiangetal.，2008；Schelleretal.，2008)，后续研究有必要扩大研究区域，收集更多与气候相关的青冈栎生长数据。

4 结论

全球气候变暖正深刻影响着森林生态系统的各方面。本研究以我国南方典型阔叶树种青冈栎为研究对象，基于芦头林场青冈栎解析木数据和气候数据，构建青冈栎单木混合效应模型，选取代表未来温室气体排放标准的典型浓度路径RCP2.6、4.5和8.5作为气候场景，预测青冈栎在未来不同时期不同气候场景下的生长趋势。本研究构建的青冈栎单木混合效应模型具有气候敏感、统计可靠和预测有效等优点，3种典型浓度路径尤其是RCP8.5场景下青冈栎胸径显著低于气候条件不变时。随着时间推移和碳排放量增加，未来气候条件很可能不利于青冈栎生长。本研究对气候变化影响青冈栎生长作了有益探索，研究结果有助于林业工作者在经营实践中应对未来气候变化所带来的挑战。

附录Appendix

附表1 5种常用生长方程的比较①

附表2 15个非线性混合效应生长模型的比较①