Multi-objective management model and harvesting strategy for Phyllostachys edulis timber forests

研究区位于浙江省湖州市安吉县杭垓镇和孝丰镇(30°30′22.61″~30°38′31.28″N，119°17′35.7″~119°32′24.76″E)。该区域属于亚热带季风气候，年平均气温为17.9 ℃，年平均降水量为1 529.5 mm，年平均日照时数为1 651.7 h。植被区划属亚热带东部常绿阔叶林亚区，森林覆盖率达70.3%，林地面积为13.5万hm²，其中竹林占6.7万hm²，毛竹林占竹林总面积的86.7%，立竹量为1.8亿株^[14−15]。

2022年8月建立野外调查样地30个，其中有25个在2022年处于毛竹林出笋大年，每个监测样地大小为20 m×20 m (400 m²)。

从2022年8月、2023年8月和2024年8月完成3次样地数据调查。每次样地调查的主要内容包括：立竹年龄、胸径、立竹度。立竹年龄用“度”进行测量，定义Ⅰ度竹为当年新竹，Ⅱ度竹为2~3年生，Ⅲ度竹为4~5年生，实际调查根据竹秆皮色法推算年龄^[16]。立竹胸径采用测量离地面1.3 m处毛竹直径获得，单位为cm。样地立竹度为样地内立竹总株数除以样地面积，单位为单位面积上的株数。

本研究中样地总地上碳储量由样地毛竹地上生物量累加再乘以毛竹含碳系数(0.5042)得到，其中单株毛竹生物量采用毛竹单株二元异速生物量生长模型^[17]来估算。计算公式如下：

式(1)中：B为单株毛竹生物量干质量(kg·株⁻¹)；D为胸径(cm)；A为年龄(度)。

对以上调查数据经最小显著差异法进行事后多重比较检验，显著性水平设定为0.05。由显著性检验发现样地各年龄、立竹度及碳储量在3期调查中没有显著差异，可用作毛竹生长更新模型的构建。

应用三参数Weibull 函数^[18]模拟新增幼竹的径阶分布，构建毛竹林更新生长模型。

式(2)中：D 为毛竹胸径(cm)；a为位置参数；b 为尺度参数；c 为形状参数。D、a、b、c＞0。对2022—2024年30个监测样地数据进行Weibull 参数拟合，通过SPSS非线性回归构建毛竹林Weibull参数a、b、c的预测模型。

①采伐策略设计。应用Matlab软件编写代码，构建不同年龄毛竹林采伐策略矩阵，主要包括采伐强度和采伐方式，其中采伐强度用0~100%表示，采伐方式分为均匀伐、两端伐、上层伐和下层伐等4种。具体策略设计如表1所示。②经营期策略设计。由汤孟平等^[19]采用时间序列谱分析法研究可知：毛竹林生长过程中的重要时间分界点存在2和6 a周期，因此本研究将毛竹材用林经营策略分为短经营期(2~6 a)、中经营期(6~12 a)、长经营期(12~18 a)等3种假定情况进行讨论。具体的经营期策略设计如下：短经营期以快速提升经济效益、同时保持碳汇潜力为目标，优先采伐大径阶Ⅲ度竹，获取高价值竹材；Ⅱ度竹生物量积累显著，选择性采伐成熟度高的个体，可优化林分空间结构；Ⅰ度竹生长速率快且固碳潜力大，保留大部分以维持竹林更新能力，疏除密度过高的幼竹，避免过度拥挤，间接提升未来竹材质量。中经营期以平衡经济收益与碳汇能力、优化资源可持续性为目标，采伐林冠层高大、完全木质化的Ⅲ度竹，经济回报率高；淘汰低效Ⅱ度竹个体，提升保留竹的质量与碳汇效率；Ⅰ度竹正处于固碳速率峰值期，采取零采伐原则。长经营期以最大化碳汇功能、兼顾稳定经济收益为目标，尽量保持40%Ⅰ度竹、30%Ⅱ度竹、30%Ⅲ度竹的合理结构，优化林分地上空间分布。由此依据不同年龄的生物学特性和经营期目标，设计采伐策略如表2所示。由于3种经营期所对应的采伐策略不同，各经营期之间可能年份重合，但不影响结果产生。③采伐过程模拟。初始经营模拟综合考虑毛竹年龄、立竹度和胸径，选取了位于杭垓镇唐舍村的9号样地。样地数据具有典型性和普遍性，其初始状态如表3所示。

根据采伐策略设计，生成保留立竹和采伐原竹2种概率矩阵，并计算伐前、伐后林分结构指标。同时更新当前林分状态，依据经营期策略设计进行多次模拟，得到相应采伐结果。

①固碳效益计算。将采伐后保留立竹和更新的立竹生物量相加，结合含碳系数(Ⅰ度竹为0.508 9 g·g⁻¹，Ⅱ度竹为0.514 0 g·g⁻¹，Ⅲ度竹为0.487 5 g·g^−1[20])，累计得到总地上碳储量，用于固碳效益评价。②经济效益计算。根据调查数据收集整理得到浙江省安吉县材用毛竹价格为12.75元·株⁻¹，以本研究30个毛竹林样地的单株毛竹平均生物量，结合竹秆的干质量比(0.44)^[21]，得到单株毛竹平均鲜质量，获得竹材单价为0.32元·kg⁻¹，贴现率设置为0.03。本研究不考虑采伐成本，根据采伐周期折现未来收益，用于经济效益评价。毛竹林净现值(V_NP)计算公式如下：

式(3)中： V_NP为净现值；t为采伐周期；I_t为采伐周期t内的收入；C_t为采伐周期t内的支出；i为贴现率；R为经营期。

①根据经营周期策略设计和采伐策略设计，将传统的“一次采伐循环”部分扩展为多周期迭代，模拟多目标情境下毛竹林长期经营效果，实现经济收益与固碳效益的总体最优。②对输出的碳储量与毛竹林净现值进行制图分析，引入帕累托前沿(Pareto Front)分析其最优解。

根据一般Weibull分布参数值的性质和可能变动范围，选取初始值进行多次迭代，得到3个参数a、b、c分别为5.312、2.895和1.101，此时决定系数(R²)为0.996，同时引入K-S拟合优度检验，得到P＞0.05，满足精度要求。

由图1可知：短经营期内2 a采伐周期碳储量随Ⅱ和Ⅲ度竹的采伐比例增加明显下降。当Ⅰ、Ⅱ、Ⅲ度竹的采伐比例分别为10%、20%、40%时，碳储量达到最大值，为17.61 t·hm⁻²；采伐比例分别为20%、50%、70%时，碳储量最小，为16.12 t·hm⁻²。4 a采伐周期碳储量较2 a采伐周期时下降约20%，当Ⅰ、Ⅱ、Ⅲ度竹的采伐比例分别为20%、20%、40%时，碳储量最大，为15.21 t·hm⁻²；当Ⅰ、Ⅱ、Ⅲ度竹的采伐比例分别为10%、50%、60%时碳储量最小，为12.51 t·hm⁻²。6 a采伐周期碳储量较4 a采伐周期下降约18%，当Ⅰ、Ⅱ、Ⅲ度竹的采伐比例分别为20%、20%、40%时，碳储量最大，为12.44 t·hm⁻²；当Ⅰ、Ⅱ、Ⅲ度竹的采伐比例分别为10%、40%、70%时，碳储量最小，为11.02 t·hm⁻²。

Figure 1.  Three-dimensional relationship diagram between short-term harvesting intensity and carbon storage

由图2可得：短经营期内2 a采伐周期毛竹林净现值随Ⅱ度竹和Ⅲ度竹的采伐比例增加明显上升，当Ⅰ、Ⅱ、Ⅲ度竹的采伐比例分别为20%、50%、70%时，毛竹林净现值达到最大值，为5 532.71元·hm⁻²；当Ⅰ、Ⅱ、Ⅲ度竹的采伐比例分别为10%、20%、40%时，毛竹林净现值最小，为3 064.68元·hm⁻²。4 a采伐周期经济收益随采伐周期延长大幅提升，在相同采伐比例时毛竹林净现值最大为10 477.26元·hm⁻²，最小为6 351.32元·hm⁻²。6 a采伐周期经济收益达到2 a周期的4倍以上，在相同采伐比例时毛竹林净现值最大为14 374.15元·hm⁻²，最小为9 996.68元·hm⁻²。

Figure 2.  Three-dimensional relationship diagram between short-term harvesting intensity and net present value

将短经营期32种采伐策略下毛竹材用林的固碳效益和经济效益综合分析，绘制帕累托前沿(图3)。图3中灰色散点表示在约束条件下的所有采伐策略解集，红色折线为帕累托前沿，指在碳储量和毛竹林净现值中至少在1个目标上优于其他解且在另一目标上不劣于其他解的解。对红色折线包含的多个散点进行对比分析可知：短经营期毛竹林净现值增速放缓，碳储量降幅提高，毛竹林生态环境功能显著降低。在Ⅱ和Ⅲ度竹的采伐比例较低时能够维持高碳储量，但经济收益较低，而连续的高强度采伐则可能会使6 a周期的碳储量跌破阈值，触发生态系统退化。同时考虑到Ⅰ度竹应尽量维持采伐比例为10%，以保障幼竹更新能力，由此选取达到短经营期固碳效益和经济效益协同最优的采伐策略为：Ⅰ度竹采伐10%，Ⅱ度竹采伐40%，Ⅲ度竹采伐60%。这时碳储量为11.14 t·hm⁻²，毛竹林净现值为12 853.18元·hm⁻²。

Figure 3.  Pareto frontiers under different short-term harvesting strategies

由图4可知：中经营期内6 a采伐周期在Ⅱ和Ⅲ度竹的采伐比例分别为30%和20%时碳储量最大，达15.37 t·hm⁻²，采伐比例分别为10%和70%时碳储量最小，为11.68 t·hm⁻²。8 a采伐周期在相同采伐比例下碳储量最大为14.23 t·hm⁻²，最小为11.27 t·hm⁻²。10 a采伐周期碳储量则在采伐比例为30%和20%时达到最大，为13.28 t·hm⁻²，在10%和50%时最小，为10.92 t·hm⁻²。12 a采伐周期碳储量同样在采伐比例为30%和20%时达到最大，为12.36 t·hm⁻²，在10%和40%时最小，为10.54 t·hm⁻²。

Figure 4.  Scatter plot of mid-term harvesting intensity and carbon storage

由图5可知：各个采伐周期的毛竹林净现值分布规律都呈现采伐强度越大毛竹林净现值越高的趋势，最高值均出现在Ⅱ度竹采伐30%、Ⅲ度竹采伐70%时，但毛竹林净现值增长速率随采伐周期增加而减缓，中经营期内8 a采伐周期增长23.1%，10 a周期增长18.3%，12 a周期则为13.5%，毛竹林净现值最高值由11 840.25元·hm⁻²提升至18 835.72元·hm⁻²。

Figure 5.  Scatter plot of mid-term harvesting intensity and net present value

对中经营期的18种采伐策略进行帕累托前沿分析(图6)，可知各周期帕累托前沿呈明显负斜率。6 a采伐周期斜率为−2.8 t·千元⁻¹，8 a周期降至−1.5 t·千元⁻¹，10 a周期毛竹林净现值＜14 000元·hm⁻²时斜率为−1.2 t·千元⁻¹，毛竹林净现值＞14 000元·hm⁻²时斜率为−2.1 t·千元⁻¹，12 a周期前沿退化为近似垂直线，斜率为−0.3 t·千元⁻¹。可以看出经济和生态的矛盾强度随时间增加而加剧，8~10 a是调整采伐策略的关键时期。由此选取中经营期固碳效益和经济效益协同最优的采伐策略为对Ⅱ和Ⅲ度竹均采伐30%，用Ⅱ度竹采伐稳定林分结构，Ⅲ度竹采伐保持经济收益。

Figure 6.  Pareto frontiers under different mid-term harvesting strategies

由图7可知：长经营期碳储量随Ⅱ和Ⅲ度竹的采伐比例增加而上升，各个采伐周期碳储量最大值都位于Ⅱ度竹采伐40%、Ⅲ度竹采伐90%时。Ⅲ度竹占主导地位，是提升碳储量的关键，当Ⅱ度竹的采伐比例≥30%时，碳储量增速放缓。

Figure 7.  Scatter plot of long-term harvesting intensity and carbon storage

由图8可知：长经营期毛竹林净现值随2种年龄采伐强度增加而上升，各个采伐周期毛竹林净现值最大值依然位于Ⅱ度竹采伐40%、Ⅲ度竹采伐90%时。Ⅲ度竹采伐比例每提升10%，毛竹林净现值平均增长12.3%，Ⅱ度竹采伐比例＞30%后边际收益递减，毛竹林净现值增幅仅3.2%。

Figure 8.  Scatter plot of long-term harvesting intensity and net present value

对长经营期的12种采伐策略进行帕累托前沿分析(图9)。可知长经营期帕累托前沿与短经营期和中经营期不同，呈现绝对最优解，即Ⅱ度竹采伐40%、Ⅲ度竹采伐90%。长经营期内12 a采伐周期时最大碳储量为11.94 t·hm⁻²，毛竹林净现值达20 583.35元·hm⁻²。14 a周期时最大碳储量为11.95 t·hm⁻²，毛竹林净现值达22 904.94元·hm⁻²，毛竹林净现值的提升主要得益于Ⅲ度竹的采伐。16 a周期时最大碳储量回到11.94 t·hm⁻²，毛竹林净现值达25 158.92元·hm⁻²，Ⅲ度竹采伐比例每增加10%，毛竹林净现值提升8.2%，碳储量同步增长1.3%。18 a周期时最大碳储量依然保持11.94 t·hm⁻²，毛竹林净现值达27 456.99元·hm⁻²，线性回归分析表明毛竹林净现值与碳储量呈强正相关，决定系数(R²)达0.98 (P＜0.001)。

Figure 9.  Pareto frontiers under different long-term harvesting strategies

毛竹林科学的经营措施包含施肥、挖笋、钩梢和择伐等^[22−23]，能够提高竹林碳汇能力和经济收益，尤其是择伐，对林分结构、植被碳储量、养分分配等具有重要影响^[24]。因此，本研究主要分析采伐年龄、采伐强度和经营期对毛竹材用林固碳效益和经济效益的影响，以探寻协同优化经营方案。本研究在设定毛竹材用林采伐经营周期时将采伐1次的最小周期设为2 a，与毛竹林大小年生长特性相符，也符合诸多学者关于采用该择伐年限来提高毛竹林立地生产力的建议^[25]。本研究建立的3种经营期的采伐状态，通过帕累托前沿量化固碳效益和经济效益替代率，对竹林经营中生态与经济协同优化进行了定量研究^[26]，证实了毛竹林经营管理中2种目标相对达到最优的可行性。

在采伐策略设计上，本研究发现：短经营期时，为了快速提升经济效益，同时保持碳汇潜力，采用Ⅰ度竹留养90%、Ⅱ和Ⅲ度竹分别实施40%和60%大径阶的采伐收获策略，可以实现固碳效益和经济效益的短期协同优化，主要原因在于对Ⅰ度竹的留养能够较好地选留壮笋、伐去弱笋，对新竹进行密度调控，释放竹林生长空间^[27]。对Ⅱ和Ⅲ度竹实施的大径阶采伐，在获得竹材收益的同时可调控竹林地上空间结构，减少聚集竞争，促进毛竹林碳汇能力提升。

在中经营期时，采用Ⅱ度竹30%小径阶采伐、Ⅲ度竹30%大径阶采伐策略，能够达到固碳效益与经济收益协同优化的目标，主要原因在于对Ⅱ度竹的小径阶采伐有选择性地移除部分弱势个体，优化林分结构，降低密度压力，为保留的Ⅱ度竹和后续新竹(Ⅰ度竹)提供更优的生长环境和光照资源，促进了碳积累效率的提升。对Ⅲ度竹实施的大径阶采伐，在获取经济价值较高的成熟竹材的同时，适度疏伐了上层林冠，有效调控了林分地上空间结构，既避免了过度采伐对碳储量的显著削减，又促进了保留竹和新竹的固碳效率^[28]，维持了竹林生产力的稳定和碳汇功能的持续增强。

在长经营期时，采用Ⅱ度竹40%小径阶采伐、Ⅲ度竹90%大径阶采伐策略，能够最大化固碳效益同时兼顾经济收益，主要原因在于对Ⅱ度竹实施的小径阶采伐强度略高于中经营期，更主动地筛选淘汰竞争力较弱的个体，优化林分年龄结构和径阶分布，为长期碳汇能力提供保障。对Ⅲ度竹实施的大径阶采伐，虽然移除了大量现存碳储量，但最大程度地释放了林分空间资源，高效地促进了新竹更新与生长，提升了林分系统的生产力，从而实现竹林碳汇功能的最大化^[29]。同时，采伐获得的大量高价值成熟竹材也带来了长期的经济收益，支撑了经营的可持续性^[30−31]。

本研究对于林分辅助采伐方式未进行深入研究，在实际操作过程中应该持续监测土壤和种群^[32]，设置林分密度、竹笋产量等多种约束指标，配合进行翻耕、施肥、劈山抚育等多种培育措施，改善竹林生长环境^[33]，为林分可持续发展提供更合理精确的经营管理优化方案。

本研究构建了毛竹材用林多目标经营模型，提出了固碳效益和经济效益协同实现的采伐与管理经营策略。针对研究区典型样地可得，短经营期时最优采伐策略为：Ⅰ度竹10%小径阶密度调控，Ⅱ和Ⅲ度竹分别实施40%和60%的大径阶采伐。中经营期时最优采伐策略为：Ⅱ度竹30%小径阶密度调控，Ⅲ度竹30%大径阶采伐。长经营期时最优采伐策略为：Ⅱ度竹40%小径阶采伐，Ⅲ度竹90%的大径阶采伐。