臧 颢 刘洪生 黄锦程 张祖栋 欧阳勋志 宁金魁

1.江西农业大学林学院 南昌 330045； 2.崇义县林业局 赣州 341300)

量化树木生长和气候的关系有助于评价气候变化对森林生长的影响(韩士杰等， 2016)。一般认为，气候因子如温度和降水等会影响树木的光合作用、呼吸作用和蒸腾作用，对树木生长具有较强的驱动作用(Toledoetal.， 2011； 余黎等， 2014)，然而树木生长和气候的关系也会受到竞争的调节作用(Ettingeretal.， 2013； Fernndez-de-Uaetal.， 2015； Zhangetal.， 2015； Fordetal.， 2017)。竞争作为评价相邻树木之间拥挤程度的指标，是种群发展的主要因子，对林木个体生存、生长、繁殖及林分的水平分布和垂直分布有着重要影响(Das， 2012； 刁军等， 2013)。竞争反映的是树木之间相互争夺资源而产生的压力(Pretzschetal.， 2010)，而气候因子中如降水等正是树木争夺的资源，因此，竞争可能会改变树木生长对气候的响应(Clarketal.， 2011； Xiangetal.， 2016)，即竞争和气候对树木生长存在交互作用(Clarketal.， 2014； Teietal.， 2014； Fordetal.， 2017)。Magruder等(2013)对美国密歇根州的红松(Pinusresinosa)进行间伐试验时发现，气候对森林生长的影响会受间伐强度的制约。Ford等(2017)在研究美国华盛顿州的山地森林胸径生长时发现，竞争会改变气候对胸径生长的影响。目前，尽管在单木和林分生长模型中同时考虑竞争和气候因素的研究众多(Bravo-Oviedoetal.， 2008； Crookstonetal.， 2010； Nunesetal.， 2011)，但探讨竞争和气候对树木生长是否存在相互作用的报道仍较少。胸径是反映树木生长状况的最重要指标(Sunetal.， 2019； 陈国栋等， 2020)，探讨竞争和气候交互作用对胸径生长的影响有助于研究气候变化背景下单木材积和生物量的估计及立地生产力的评价等。

杉木(Cunninghamialanceolata)是我国南方最重要的速生用材树种之一，据第九次全国森林资源清查结果(国家林业和草原局， 2019)，杉木人工林总面积现已达990.20万hm2、蓄积达7.55亿m3，占全国人工林总面积的17.33%、总蓄积的22.30%。我国是一个木材资源严重缺乏的国家，木材对外依存度高(许传德等， 2015)，在气候变化和木材短缺的压力下，发展人工用材林已成为保障国家生态安全和解决木材供需矛盾的战略选择。鉴于杉木在全国用材林中的重要地位，探索竞争和气候对杉木人工林生长的交互作用，对在气候变化背景下科学经营人工林、实现林业“双增”目标具有重要意义。

鉴于此，本研究基于杉木人工林固定样地数据，采用潜在生长量修正法构建包含竞争指标和气候因子的胸径生长模型，并尝试分析竞争和气候及其交互作用对杉木人工林胸径生长的影响，以期为气候变化背景下模拟抚育间伐、择伐后保留木的生长变化奠定基础，为森林适应性经营中科学合理地对杉木人工林进行间伐、择伐提供依据。

1 研究区概况、数据与研究方法

1.1 研究区概况

研究区位于江西省赣州市的南康区(25°28′—26°14′N，114°29′—114°55′E)、崇义县(25°24′—25°55′N，113°55′—114°38′E)和上犹县(25°42′—26°10′N，114°01′—114°40′E)，罗霄山脉与大庚岭相遇的南端，为诸广山余脉的一部分。属亚热带季风气候，年均气温20.4 ℃，年均降水量2 265.5 mm，以低山、丘陵地貌为主，土壤主要为黄壤、红壤等。研究区种植有大量人工针叶林，主要树种有杉木、马尾松(Pinusmassoniana)、木荷(Schimasuperba)、南酸枣(Choerospondiasaxillaris)、丝栗栲(Castanopsisfargesii)、米槠(Castanopsiscarlesii)、米锥(Castanopsischinensis)等。

1.2 数据资料

1.2.1 固定样地数据 固定样地数据来自南康区、崇义县和上犹县39块杉木人工林样地。样地布设时间为2011—2015年，每块样地进行过1～2次复测，调查间隔期1年。主要调查样地的地理坐标、地形、坡度、海拔、土壤类型和种植时间，并对胸径5 cm以上的树木进行每木检尺，记录树种，测量胸径、树木相对位置等因子，其中， 20块样地没有调查树高。依据胸径测量数据，计算林木2次调查间的胸径生长量，共计6 181个胸径生长量观测值。样地调查因子的详细统计量见表1，样地分布示意见图1。

表1 数据统计①

图1 样地分布示意

1.2.2 气候数据 气候数据依据地理坐标和海拔从Wang等(2017)编写的ClimateAP软件中提取，空间分辨率为4 km×4 km。考虑到固定样地的布设和复测时间多在6月，为模拟气候对胸径生长可能存在的影响，将每年7月至次年6月视为调查间隔期，提取调查间隔期内的月值数据，汇总计算各样地调查间隔期的平均温度、最高温度、最低温度、降水量和大于5 ℃的积温。详细统计量见表1。

1.3 研究方法

潜在生长量修正法是构建单木生长模型的一种常用方法(孟宪宇， 2006)，其基本思路是： 模拟林木的潜在生长量(无竞争压力时的生长量)，利用反映林木所受竞争压力的修正函数对潜在生长量进行调整和修正，从而得到林木实际生长量，用数学形式可以表示为：

id=idpot×M(CI)。

(1)

式中： id为胸径生长量； idpot为胸径潜在生长量；M(CI)为反映林木所受竞争压力的修正函数，且满足0≤M(CI)≤1。

1.3.1 潜在生长量模拟 模拟潜在生长量的步骤如下：

1) 选择基础方程，以气候和地形因子反映立地质量。考虑到林木潜在生长量受立地质量影响(刘洋等， 2012； 王冬至等， 2015； 雷相东等， 2018)，而立地质量可用环境因子评价，如气候(张海平等， 2017； Jiangetal.， 2015； Shenetal.， 2015)和地形因子(杜纪山， 1999)，故本研究尝试以气候和地形因子反映立地质量。综合数据分析和以往研究成果(Pretzschetal.， 2010； Weiskitteletal.， 2011； Burkhartetal.， 2019)，选用柯列尔方程(斯瓦洛夫， 1983； 孟宪宇， 2006)作为潜在生长量的基础方程，采用再参数化方法，将柯列尔方程中的参数表示为关于5个气候因子(调查间隔期的平均温度、最高温度、最低温度、降水量和大于5 ℃的积温)和2个地形因子(海拔和坡度)的函数，通过最小二乘法求解参数，并剔除不显著因子(显著性水平取0.05)，从而得到含环境因子的生长方程。柯列尔方程的数学形式如下：

y=a0xa1e-a2x。

(2)

式中：y和x分别为因变量和自变量；a0、a1和a2为模型参数。

2) 利用分位数回归模拟含环境因子的生长方程。分位数回归模型可拟合因变量任意分位数的预测值，并避免选取样点的主观因素，常应用于最大化和最小化问题(高慧淋等， 2016； Haoetal.， 2007)。潜在生长量，即理论生长最大值，可理解为某一特定树种在特定立地条件下无竞争压力的疏开木生长量(孟宪宇， 2006； Pretzschetal.， 2010)。本研究尝试利用分位数回归模拟含环境因子的生长方程，并将预测结果作为潜在生长量，分位点τ取0.90、0.95和0.99。模型形式如下：

idpot=(a00τ+∑a0jτEj)D(a10τ+∑a1jτEj)e-(a20τ+∑a2jτEj)D。

(3)

式中：D为胸径；Ej为第j个环境因子；a00τ、a0jτ、a10τ、a1jτ、a20τ和a2jτ为第τ分位点模型参数。

第τ分位点模型参数求解可通过下式计算得到(Koenkeretal.， 1978)：

(4)

1.3.2 修正函数构建 为探讨竞争对树木生长的影响，本研究选择6个备用竞争指标用于构建修正函数，以反映胸径生长对竞争的响应： 1) 林分密度指标，包括公顷株数N、公顷断面积BA和林分密度指数SDI； 2) 单木竞争指标，包括相对胸径Rd、大于对象木的断面积之和BAL、简单竞争指数Hegyi(孟宪宇， 2006； Weiskitteletal.， 2011； Zhouetal.， 2019)。部分竞争指标计算公式如下：

(5)

(6)

(7)

(8)

式中：D0为标准平均直径，本研究取10 cm(孟宪宇， 2006)；Dg为林分平均胸径；Di和gi分别为对象木的胸径和断面积； Area为样地面积；gk为样地中大于对象木的第k株样木断面积；Dj为第j株竞争木胸径； Distij为第j株竞争木与对象木的距离；m为竞争木株数。

采用3种方法确定竞争木(汤孟平等， 2007； Kuehneetal.， 2019)： 1) 以对象木为中心，将与对象木最近的1～10 株相邻木作为竞争木； 2) 以对象木为中心，以1～10 m为半径画圆，将各半径内所有非对象木视为竞争木； 3) 选用Voronoi图确定竞争木。考虑到边缘效应的影响，运用8邻体平移法计算Hegyi(申瀚文等， 2012； Ripley， 1979)。

修正函数M(CI)采用指数函数ef(CI)形式。为求得修正函数中的参数，同时考虑到个别为0 cm的胸径生长量无法取对数，参考以往研究成果(李春明， 2012； Pretzschetal.， 2010)，将式(1)整理为ln(id/idpot+1)=f(CI)形式，以避免出现负值，f(CI)为关于竞争指标的线性函数。为反映竞争和气候可能存在的交互作用，在竞争指标构建的修正函数中添加竞争和气候交互作用项，并与不添加竞争和气候交互作用项的函数进行对比。不考虑交互作用和考虑交互作用的修正函数参数如下：

(9)

(10)

式中： CI为竞争指标；CLi为第i个气候因子；b0、b1、ci为模型参数。

考虑到嵌套数据结构可能带来的变量独立性缺失问题(李春明， 2012； 彭娓等， 2018)，基于选定的修正函数形式，采用混合效应模型进行重建。混合效应模型一般形式为：

Yi=f(β,ui,Xi)+εi。

(11)

式中：Yi和Xi分别为因变量向量和自变量向量；εi为误差向量；β和ui分别为固定效应参数向量和随机效应参数向量，且εi～N(0,Ri)、ui～N(0,Ψ)，Ri和Ψ分别为样地内的方差协方差矩阵和随机效应参数的方差协方差矩阵。本研究中，随机效应参数的方差协方差结构采用广义正定矩阵形式。

对混合效应模型进行检验时，未包括在建模样本中的样地须计算其随机效应参数。通过泰勒公式将已有非线性混合效应模型线性化，可以推导出随机效应参数的经验贝叶斯估计值(Davidianetal.， 1995)，公式如下：

(12)

1.3.3 模型评价 采用交叉验证法评价模型，即以全部39块样地所有杉木的单木胸径生长数据为基础，每次将其中38块样地数据作为建模样本，另外1块样地数据作为检验样本进行模型评价，重复39次。使用平均绝对误差(mean absolute error，MAE)、平均相对误差绝对值(relative mean absolute error，RMAE)和平均预估误差(mean predicted error，MPE)评价模型估计精度，计算公式(曾伟生等， 1999； 2011； Ouetal.， 2019)如下：

(13)

(14)

(15)

2 结果与分析

2.1 潜在生长量与环境因子的关系

采用再参数化方法拟合得到含环境因子的潜在生长量方程如下：

idpot=(a00+a01Tmin+a02PPT+

a03Elevation)Da1e-a2D

(16)

式中：Tmin、PPT和Elevation分别为调查间隔期的最低温度、降水量和样地海拔。

剔除不显著的环境因子后(显著性水平=0.05)，最小二乘法的参数估计结果见表2。考虑到共线性可能对参数估计造成的影响，采用方差膨胀因子分析3个环境因子的共线性程度，最低温度、降水量和海拔的方差膨胀因子分别为4.35、4.59和4.84，即可认为最低温度、降水量和海拔之间不存在强共线性。

基于式(16)，利用分位数回归(τ取0.90、0.95和0.99)模拟胸径潜在生长量，参数估计结果见表2。调查间隔期的最低温度、降水量和海拔对杉木人工林胸径潜在生长量存在显著影响，且最低温度和降水量与潜在生长量呈正相关关系，海拔则呈负相关关系。

2.2 竞争对胸径生长的影响

不同竞争指标构建的修正函数对胸径生长量的拟合效果见表3。可以发现，含单木竞争指标的模型精度优于含林分密度指标的模型，与距离无关的单木竞争指标的生长模型估计精度略优于与距离有关的单木竞争指标的生长模型，含BAL的模型评价指标最优，其中，MAE降低0.71%～20.07%，RMAE降低12.52%～33.63%，MPE降低3.21%～25.16%。

表2 式(16)的参数估计

表3 含不同竞争指标的模型评价①

图2 2种方法计算的Hegyi竞争指数对胸径生长量的影响

通过对比不同分位数拟合杉木人工林胸径潜在生长量对模型精度的影响(图2)发现，分位数取0.90时，模型评价指标表现最差，随着分位数增加，MAE、RMAE和MPE逐渐降低，仅分位数取0.95时含BA和Rd模型的MPE大于分位数取0.90时的模型。对比3种分位数下所有竞争指标模型，总体来看，分位数取0.99时含BAL的模型评价指标最优，故将其作为基础模型进一步分析竞争和气候交互作用对杉木人工林胸径潜在生长量的影响。

2.3 竞争和气候交互作用对胸径生长的影响

结合上述分析，以分位数取0.99时含BAL的模型为基础模型，进一步在修正函数中添加BAL与5个气候因子(调查间隔期的平均温度、最高温度、最低温度、降水量和大于5 ℃的积温)的交互作用项，各模型的评价指标见表4。与不含交互作用项的模型(表3)对比，除RMAE外，考虑交互作用项后模型的评价指标更优，其中，MAE降低0.60%～18.69%，MPE降低0.12%～9.72%。考虑BAL与平均温度、最高温度、最低温度交互作用的模型RMAE和MPE最小，即同时考虑BAL与平均温度、最高温度、最低温度交互作用的模型最优。

表4 含不同交互作用项的模型评价

此外，考虑到嵌套数据结构可能带来的变量独立性缺失问题，采用混合效应模型重新模拟最优模型的修正函数部分。由于部分样地仅有2期调查数据，即只有1期生长量数据，故未考虑单木水平的随机效应。最终构建的胸径生长模型形式如下：

(17)

式中： idijk为第i块样地中第j株样木第k次调查间隔期的胸径实测生长量； BALijk和Dijk分别为第i块样地中第j株样木第k次调查间隔的期初BAL和胸径；Tmeanijk、Tmaxijk、Tminijk和PPTijk分别为第i块样地中第j株样木第k次调查间隔期的平均温度、最高温度、最低温度和降水量； Elevationi为第i块样地海拔；a00τ、a10τ、a20τ、a0jτ、a1jτ、a2jτ为第τ分位数回归模型参数(τ取0.99)；b0、b1、b2、b3、b4为固定效应参数；βi为样地水平的随机效应参数；εijk为误差项。

式(17)中潜在生长量参数估计结果见表2，修正因子参数估计结果见表5。

添加随机效应参数后，模型MAE=0.071 1 cm，RMAE=20.37%，MPE=4.90%。相比无随机效应参数模型(表4)，MAE提升35.71%，RMAE降低4.36%，MPE降低37.26%。

基于最终胸径生长模型，各径阶残差分布如图3所示。总体来看，该模型能够较好预测各径阶的胸径生长量，但径阶小于18 cm时，对部分林木实际生长量存在低估情况； 而径阶大于20 cm时，胸径估计结果的精度和稳定性更好。

表5 最终模型修正函数的参数估计①

基于最终模型固定效应参数，不同环境条件下胸径潜在生长量和修正函数的变化如图4、5所示。就潜在生长量而言，胸径20 cm时，潜在生长量达到最大值，且因环境条件改变潜在生长量的最大值也随之变化，其中，调查间隔期的最低温度影响最大，降水量影响最小。就修正因子而言，随着林木所受竞争压力不断增加(BAL增加)，修正函数呈递减趋势，且在BAL>8 m2·hm-2时，减小速度变缓，并逐渐趋于稳定。此外，BAL在0～30 m2·hm-2之间时，平均温度降低会导致修正函数减小，BAL在0～23 m2·hm-2之间时，最高温度增加会导致修正函数减小，BAL在0～15 m2·hm-2之间时，最低温度增加会导致修正函数减小，即平均温度降低、最高温度增加和最低温度增加均会加剧竞争对胸径生长的影响，且该影响随着BAL增加呈先加强后减弱的趋势，BAL≥30 m2·hm-2后不再有明显影响，因此可认为： 对林分中承受竞争压力较小(BAL接近0 m2·hm-2)和竞争压力较大(BAL≥30 m2·hm-2)的林木而言，竞争是对胸径潜在生长量修正的主要因子，气候影响不大，而对林分中BAL在0～30 m2·hm-2之间的林木而言，竞争对胸径潜在生长量的修正效应会受气候影响。

图3 各径阶残差分布

图4 不同环境条件下胸径潜在生长量的预测

图5 不同环境条件下修正函数的平均响应

3 讨论

潜在生长量修正法是构建单木生长模型的一种常用方法，其关键在于构建林木潜在生长量方程及对潜在生长量进行调整的修正函数。理论上，潜在生长量方程应基于始终处于无竞争和外界压力情况下自由生长的疏开木调查数据，但这一点在实际中较难满足。分位数回归模型可拟合因变量任意分位数的预测值，其在国内外林业研究中已有应用，且在最大化问题如计算最大密度限上应用较多(Zhangetal.， 2005； Xueetal.， 2015； 高慧淋等， 2016)，因此，采用分位数回归模型模拟胸径生长量时，通过设定高分位点(如0.99)即可得到该分位点的胸径生长量预测值，并可将其视为胸径潜在生长量。从本研究结果看，最终胸径生长量模型式(17)估计精度良好，因此可认为高分位点对含环境因子的胸径潜在生长量描述较好。Pretzsch等(2010)利用分位数回归模型模拟中欧4个树种的潜在生长量，结果发现，该模型能较好反映出不同立地条件下的胸径潜在生长量。

因部分样地没有调查树高，因此本研究在构建胸径生长模型时，用环境因子描述立地质量对林木潜在生长量的影响，结果显示，调查间隔期的最低温度、降雨量和样地海拔对杉木人工林胸径潜在生长量具有显著影响，且均呈正相关。很多学者也采用气候和地形等环境因子反映立地质量，如杜纪山(1999)在构建落叶松(Larixgmelinii)胸径生长模型时，使用海拔、坡率(坡度的正切值)、坡向及坡率与坡向之间的组合项反映立地质量，并发现坡度的效果最优。Ou等(2019)在模拟落叶松-云杉(Piceaasperata)-冷杉(Abiesfabri)混交林胸径生长量时，考虑坡度、海拔、坡向3个地形因子及44个气候因子，结果发现尽管地形因子能提升模型精度，但效果弱于气候变量。

竞争是树木生长的主要影响因素之一(彭娓等， 2018)，因此竞争指标常作为自变量用于胸径生长模型。依据是否反映出单木的竞争状况差异，可将竞争指标分为反映林分内单木平均竞争状况的林分密度指标和描述各单木不同竞争状况的单木竞争指标(孟宪宇， 2006)。本研究中，通过比较6个竞争指标对胸径生长量的模拟后发现，描述林木竞争状况差异的单木竞争指标明显优于反映林木平均竞争情况的林分密度指标，其中大于对象木的断面积之和拟合效果最好。尽管有研究显示树木胸径生长还与其周围树木的距离具有很高相关性(刘洋等， 2012； 房晓娜等， 2015； Coomesetal.， 2007)，但也有学者发现与距离无关的竞争指标在拟合单木胸径生长时效果更好(Kuehneetal.， 2019)，这可能与研究对象的差异和竞争指标的计算有关。本研究中，样地边缘林木的竞争木运用8邻体平移法确定，这也可能导致计算得到的Hegyi指数与林木的真实竞争状况存在差异，因而降低了对胸径生长量的拟合效果。另外，由于不需要考虑林木间的相对位置，因而在生产实际中，与距离无关的竞争指标应用也更为方便。很多学者在构建胸径生长模型时会选用与距离无关的竞争指标，且估计精度均较高(李春明， 2012； Zhaoetal.， 2013； Fordetal.， 2017)，也有一些学者研究发现，林分密度指标和单木竞争指标均能提高胸径生长模型的估计精度，可同时加入模型中，如吕勇等(1999)在模拟杉木人工林单木胸径生长时，选取1个林分密度指标(相对植距)和1个单木竞争指标(相对优势度)； 马武等(2015)采用逐步回归法，以方差膨胀因子和调整决定系数为依据，筛选用于构建蒙古栎(Quercusmongolica)天然林胸径生长模型的自变量，备选的7个竞争因子(林分断面积、林分密度指数、大于对象木的树木断面积之和、相对胸径、大于对象木的所有林木胸径平方和、对象木直径与林分中最大林木直径之比以及郁闭度)中有1个林分密度指标(林分断面积)和1个单木竞争指标(大于对象木的树木断面积之和)被选入最终模型中。

4 结论

为分析胸径生长受竞争和气候条件的影响，本研究基于江西省赣州市南康区、崇义县和上犹县杉木人工林固定样地数据，采用潜在生长量修正法构建胸径生长模型，结果发现，海拔、调查间隔期的最低温度和降水量对胸径潜在生长量具有显著影响，胸径20 cm时，潜在生长量达到最大值。竞争作为林木个体生存、生长的重要影响因子之一，也对模型拟合好坏有较大影响，且单木竞争指标估计精度优于林分密度指标，与距离无关的单木竞争指标估计精度略优于与距离有关的单木竞争指标。通过分析竞争和气候交互作用对潜在生长量的修正效果发现，竞争和气候交互作用可在一定程度上提高估计精度，但不同竞争指标和气候因子构成的交互作用效果并不一致，其中，同时考虑平均温度、最高温度和最低温度与大于对象木的断面积之和构成的交互作用的模型预估精度最好，且平均温度降低、最高温度和最低温度增加会加剧竞争对胸径生长量的影响。在此基础上，下一步可考虑探讨在气候变化下如何设计间伐、择伐才能最大化经营目的。