Simulation and prediction of tree species composition change of forests in Taizhou City

台州市(28°01′~29°20′N，120°17′~121°56′E)位于浙江省中部沿海，总面积约10 050 km²。该区属中亚热带季风气候，年均日照时数为1 800~2 037 h，全年降水日数为132~171 d，年均降水量为1 632 mm。夏、秋季受台风影响比较大。地带性植被为中亚热带常绿阔叶林，主要的阔叶乔木树种有壳斗科Fagaceae的甜槠Castanopsis eyrei、山茶科Theaceae的木荷Schima superba等，针叶树种主要有马尾松Pinus massoniana、杉木Cunninghamia lanceolata等。主要森林类型有常绿阔叶林、落叶阔叶林、针阔混交林、常绿落叶阔叶混交林、竹林等^[14]，森林覆盖率为67.0%。

选用台州市6期(1994、1999、2004、2009、2014和2019年)森林资源连续清查固定样地和固定样木的资料作为研究数据，各期固定样地总数(含非林地、无林地等非乔木林样地)分别为379、390、388、388、388和388个。1994和1999年2期样地数和后面各期相比略有差异，是因为前3期每期有1/3的样地进行了移动。移动部分样地的目的是为了检验固定样地是否受到特殊对待，即样地内部和外部是否具有相同的经营措施，结果表明没有特殊对待的情况发生，即样地固定后对总体的代表性不会改变。固定样地为800 m²的正方形。森林资源连续清查体系是基于概率抽样理论的机械抽样，每5 a复测1次，是森林资源调查中最具权威性的数据，对于森林的动态变化分析及预测十分有利^[15]，同时也便于更有效地传播和利用现有数据进行可持续的森林生态系统经营和管理^[16]。本研究的对象是乔木树种，所以只选择其中的乔木林样地，根据地类选择样地不改变其随机性和代表性。由于经济林比较特殊，所以本研究不包括乔木经济林。

由于树种众多，对每一树种进行模拟是不现实的，因此需要将具有共同特性的树种聚集为组进行模拟^[17]。将样木分为3个树种组：①松类包括马尾松、黑松P. thunbergii、黄山松P. taiwanensis、湿地松P. elliottii等，其中马尾松占绝对多数；②杉类包括杉木、柳杉Cryptomeria fortunei、水杉Metasequoia glyptostroboides、池杉Taxodium ascendens、红豆杉Taxus wallichiana var. chinensis等，其中杉木占绝对多数；③阔叶类包括栎属Quercus、栲属Castanopsis、枫香Liquidambar、木荷、樟科Lauraceae等植物。

将乔木林样地按优势树种组(样地中蓄积比例最大的树种组)统计，乔木林样木胸径最大值、均值、中位数和标准差等见表1。其中，最小值均为森林资源连续清查规定的5.0 cm，故没列出。

将所有乔木林样地上的样木合在一起，统计3个树种组的总样木数和总材积及其比例。材积数据来自森林资源连续清查数据库，是根据一元材积表计算的。这些数据将用于建立自约束模型。统计各树种组的进界木和采枯木株数(受松材线虫病危害的林木在死亡后及时得到清理，所以调查时根据伐根难以区分采伐木和枯损木，所以本研究将两者合在一起，简称采枯木)，用于解释模型的预测结果。

设y₁、y₂和y₃分别为松类、杉类和阔叶类树种组的株数比例或材积比例，t为年份，可得：

式(1)中：λ₁、λ₂、γ₁和γ₂为模型参数。由于变量和参数自动满足0＜y₁，y₂，y₃＜1，y₁+y₂+y₃=1，λ₁，λ₂＞0，γ₁，γ₂＜0，所以本研究将式(1)命名为非线性方程组自约束树种结构模型。式(1)推导过程繁复，本研究略去推导过程的介绍。年份t的起点可以自由设定，本研究为了方便取起点，即1994年时的t取5，1999年时的t取10等，当然也可以取实际数1994、1999等。y₁，y₂，y₃存在理论极限值y_1∞，y_2∞，y_3∞(极限的点估计)，根据式(1)可得：

极限作为今后的发展趋势，具有重要意义。设$ {y_\infty } = {({y_{1\infty }},{y_{2\infty }},{y_{3\infty }})^{\text{T}}} $为极限的点估计，各极限点估计的方差D为：

式(3)中：$ \sigma _{_{{\lambda _1}}}^2 $、$ \sigma _{_{{\lambda _2}}}^2 $、$ {\sigma _{{\lambda _1}{\lambda _2}}} $分别为参数λ₁、λ₂的方差和它们之间的协方差。

式(1)有4个独立参数，而用于估计参数的数据只有6组，自由度只有2，所以只能采用小样本t分布方法。同时由于比例数不可能小于0、大于1，所以应该采用截尾t分布。设自由度为k的t分布概率密度函数为$ p(t,k) $，则截尾t分布概率密度函数p_a(t, k)为：

用式(4)进行区间估计，具体操作利用Excel的对应函数实现。设α为显著性水平，$ {t_{\rm{l}}}{(k)_{\alpha /2}} $、$ {t_{\rm{r}}}{(k)_{\alpha /2}} $满足以下条件：

则极限真值$ {Y_\infty } $的1−α双侧区间估计为：

若$ {t_{\rm{l}}}{(k)_\alpha } $、$ {t_{\rm{r}}}{(k)_\alpha } $满足以下条件：

则极限真值$ {Y_\infty } ＞ {t_{\rm{l}}}{(k)_\alpha }\sqrt {D({y_\infty })} = {\Delta _{{\rm{l}} \cdot \alpha }} $的概率为1−α，同理，$ {Y_\infty } ＜ 1 - {t_{\rm{r}}}{(k)_\alpha }\sqrt {D({y_\infty })} = 1 - {\Delta _{{\rm{r}} \cdot \alpha }} $的概率为1−α。

从表2和表3可以看出：株数比例和材积比例表现出了相同的规律，松类比例持续下降，阔叶类比例持续上升，杉类比例呈先上升后下降的变化趋势。

从表4可见：松类和杉类的进界木株数在下降，阔叶类的进界木株数在上升，且总量远大于松类和杉类之和。5 a进界率和5 a采枯率分别为监测周期内发生的进界木总株数和采枯木总株数与初期的活立木总株数的比值。从表5可见：松类平均单株材积单调上升，杉类早期有波动，后期也是明显上升，阔叶类则是先下降后上升，树种组的总平均值也整体呈上升趋势。

基于表2和表3数据，利用Stata 15.1拟合自约束模型，结果见表6、图1和图2。从图1和图2可知2019年之后40 a株数和材积的预测情况，株数比例模型和材积比例模型的$ {\sigma _{{\lambda _1}{\lambda _2}}} $分别为0.239 4和0.051 9。

Figure 1.  Graph of predicted percentages of the numbers of trees by species group

Figure 2.  Graph of predicted percentages of the volumes of trees by species group

从表7可见：区间估计的可靠性均为95%，即α统一取0.05。各t值均为自由度为2的标准化截尾t分布取值。双侧分布区间估计给出95%分布的下限和上限，例如松类极限值在1.50%~7.63%的概率为95%。单侧分布区间估计下限表示极限值大于这个下限的概率为95%，例如松类株数比例的极限值大于2.10%的概率为95%。同理单侧分布上限表示极限值小于这个上限的概率为95%，例如松类株数比例的极限值小于6.50%的概率为95%。

本研究表明：松类的株数比例和材积比例均在持续单调下降，杉类均是先上升后下降，阔叶类均在快速单调增长。自约束模型很好地拟合了1994—2019年25 a间台州市乔木林3大树种组的株数比例和材积比例的变化情况。

松类、杉类、阔叶类株数比例的极限点估计分别为4.09%、7.02%、88.89%，材积比例为10.07%、11.94%、77.99%。按95%的可靠性，松类、杉类、阔叶类株数比例极限分布双侧区间估计分别为1.50%~7.63%、2.67%~12.86%、79.83%~95.71%，材积比例分别为3.64%~18.90%、4.35%~22.31%、59.90%~91.76%。单侧区间估计表明：松类、杉类、阔叶类株数比例极限分别大于2.10%、3.71%、82.72%，分别小于6.50%、10.99%、94.07%，材积比例分别大于5.13%、6.12%、65.52%，分别小于16.09%、19.01%、88.51%的可能性为95%。可见，阔叶树将占绝对优势，松类的比例将会是最低的，但永远不会消失，将以低比例一直存在，这与已有研究结论相似^[18]，杉类的比例比松类略高，也将以低比例一直存在。

本研究还表明：松类的单株材积在单调上升，尤其近些年上升速度更快。杉类整体也在上升，且也是近年上升更快。阔叶类是先下降后上升，但上升速度比松类和杉类低，尤其和松类相比低更多。阔叶类与松类单株材积的比值，从1994年的0.901下降到2019年的0.425，阔叶类与杉木的比值从1994年的1.295下降到2019年的0.626，中间有波动。单株材积的变化规律导致了松类和杉类的株数极限比例分别小于各自的材积比例，而阔叶类正好相反，株数比例大于材积比例。

平均材积的变化差异主要原因在于进界木株数的差异。本研究表明：从1999—2004年这个复查期开始，松类的进界木株数持续大幅减少，进界率也在下降，造成松类小径阶林木的占比减少，导致松类的平均单株材积上升，杉类也有类似的情况，而阔叶类进界木持续增加，导致其平均单株材积增长缓慢。

本研究表明：以松类为优势的样地数从2004年的105个下降到2019年的41个，减少了60%；株数占比和材积占比分别从1994年的0.729和0.785单调下降至2019年的0.140和0.246，下降速度很快。松类资源下降的主要原因是松材线虫病的危害^[19]。2019年，浙江省新增松材线虫病发生面积为29.88万hm²，新增病枯死松树达413.74万株^[20]。从本研究可以看出，松类的采枯率高出阔叶类很多。随着疫情的加重和疫区的扩大，正常的商业采伐和抚育采伐比例越来越小，而因疫情死亡的比例越来越高。因为松类树种易受疫情危害，所以松类人工造林急剧下降。由于林分整体密度在增加，整体竞争逐步加剧，马尾松等阳性树种在正常的自然演替中处于相对劣势地位，这也导致了随着演替的进程马尾松枯死现象的增多^[21−24]。

杉类树种的株数比例及材积比例在1994至2004年间增加，其后缓慢减少。杉木优势的样地数由1994年的22个增加至2004年的54个，而同时期的松类样地数变化并不大，说明这个时期台州市造林以杉类为首选。此后样地数略有下降。杉类中杉木占绝大部分，近年由于经济效益的下降，杉木造林不多。但杉木具有萌芽的特性，样地数减少的速度比松类低，进界率近年也比松类略高，但比阔叶类低很多。杉类由于不受松材线虫病的危害，采枯率比松类低很多，但略高于阔叶类。

阔叶类占优势的样地数从1994年的7个增加到2019年的95个，株数比例和材积比例分别从1994年的0.052和0.051上升到2019年的0.612和0.457。进界木株数绝对值从1994年的833株上升到2019年的4 946株，阔叶类进界率是松类的4.8倍以上，对比杉类，除了1994—1999是1.6倍外，其余是3.3倍以上。而最近3个监测期的采枯率低于松类的2/5。阔叶林资源的增加主要有以下几个方面：①为了公益林建设，大面积的封山育林和林种结构调整，使得阔叶幼林和部分灌木林发展成为阔叶林；②随着近自然经营的兴起，出现了大量的“针改阔”，通过针叶树采伐或死亡后补植阔叶树种或促进阔叶树幼树生长的形式调整林分的树种组成^[25]；③随着人们对森林生态系统认识的提高，以及对病虫害的考虑，有意识地增加了造林时对阔叶树种的选择^[26]；④松材线虫病危害导致了松林资源的损失，但同时为阔叶林资源的发展提供了空间；⑤随着演替的不断进行，密度增大，马尾松等阳性先锋树种的竞争力下降，并逐渐被淘汰，阔叶树种的优势逐渐扩大，随之达到树种结构的平衡并趋于稳定^[21−24]。

一般来说常绿阔叶林是中亚热带的顶级森林群落，具有最好的生态功能，但适量的松类和杉类树种有利于生物多样性的保持和生态功能的优化。在控制松材线虫病的前提下，松类是很好的景观树种和用材树种，杉类也是很好的用材树种，尤其是大径材培育的优先选择树种。本研究表明：松类和杉类的比例虽然在下降，但始终会以一定的比例存在而不会消亡。

虽然3个树种组的变化各异，但3个树种组的总株数和总材积均在稳定增长。乔木林的总样地数从2004年开始基本稳定，但总株数和总材积却在持续增长，说明林分的相对密度在持续增大。需引起重视的是阔叶类采枯木在快速增加，这其中采伐木的比例不大，主要是由于密度过大导致枯死增加，这说明需要适当加强阔叶林的抚育。

本研究是一个区域的乔木林树种结构变化规律，需要大量的数据支撑。由于2个原因没有将数据分为建模和检验2个部分：一是乔木林样地不足200个，如果分成2个部分则会造成建模和检验的数据都不够充分，得不到可靠的结果；二是6期数据的乔木林样地并不是固定不变的，前3期每期样地都有1/3移动，后面虽然不再移动，但样地的地类有变化，有的从乔木林变为其他地类，有的从其他地类变成乔木林，乔木林样地不固定，乔木林总样地数也在变化，这给数据建模和检验造成了麻烦。所以本研究采用基于建模结果进行可靠性理论分析的方法。

式(1)具有4个独立参数，而建模数据只有6组，自由度只有2，所以本研究采用自由度为2的小样本t分布方法进行极限分布的区间估计。由于比例数据不可能小于0、大于1，所以又进行了截尾调整，即采用了截尾t分布模型。松类的株数比例极限分布95%的可能性在1.50%~7.63%，大于2.10%的可能性为95%，是具有实际参考价值的。

非线性方程组自约束树种结构模型适用于大区域树种组结构的模拟和预测。它反映区域的整体变化情况，但不能反映局部细节。25 a间台州市乔木林资源总量在持续增加，但呈现“针减阔增”的现象：①受松材线虫病危害以及自然演替中的不利因素影响松类比例持续降低；②受市场因素影响，杉类的比例在1999—2004年间达到最大值后开始下降，但幅度比松类小；③阔叶类的比例呈现持续、稳定增加趋势，将占绝对优势。松类和杉类的比例虽然在下降，但最终会稳定在一定的比例而不会消亡。这表明台州市森林生态系统正朝着以阔叶树种为主、结构更合理的植被群落方向正向演替。