Dynamic Visual Simulation of Stand Multi-objective Management Based on Simulated Annealing Algorithm

实验区设置在湖南省攸县黄丰桥国有林场，林场地处113º04´~113º43´E、26º43´~27º06´N之间，最低海拔115 m，最高海拔1 270 m。隶属于亚热带季风气候，平均气温17.8 ℃，年降水量1 410.8 mm，森林覆盖率为90.07%，是湖南省重点杉木大径材推广示范基地。

用常规的测树学方法对选择的样地进行每木检尺和单木定位，测量其胸径、树高、冠高、冠幅、活枝下高、生长情况以及树木的相对位置$ ({\rm x},{\rm y},{\rm z})$，按照1年时间为间隔期，采集了2012—2017年的连年调查数据，总结其中地位指数相同的5块样地数据如表1所示：

根据惠刚盈等人提出的结构化森林的经营思想，健康的林分结构可以通过角尺度、大小比数、混交度、拥挤度4个结构参数加以描述，分别代表了林分的分布状况、林分的竞争情况、林分种间的隔离程度、林木的营养空间大小，根据前人的研究成果，健康的林分结构在这4个参数中均可以用林分均值定量化来表示^[11]。本研究以杉木人工纯林为研究对象，考虑到人工纯林混交度为零，所以剔除混交度参数，建立空间结构函数，公式如下，当空间结构函数Q(g)越小，林分结构越接近健康水平。

其中g为样地林木向量；$ \bar W$为林分平均角尺度；$ {\bar U}_a$为调整后的林木大小比数；$ {\bar U}_{a-0}$为调整前的大小比数；$ {\bar C}$为林分平均拥挤度。

通过空间结构函数来进行间伐木的判断，按照惠刚盈等人提出的方法需要对每一株木的结构参数进行计算，并且需要对模拟出的多种预间伐后的方案进行比较，虽然这种穷举法可以得出经营的最优解，但是，一旦树木的数量较大，计算起来可能会出现难以忍受的等待时间。考虑节省人力物力的条件下，满足森林经营整体需求的次优解也是被允许的，因此本研究设计了一种基于模拟退火算法(SAA)的最优间伐经营方案选择，模拟退火算法是求解组合问题的主要方法之一，它主要是基于Monte-Carlo迭代求解策略的一种随机寻优算法，对于求解最优组合快速、稳定、有效^[19-20]。基于SAA的最优间伐经营方案选择过程如下：

初始化温度T，设置温度下降的速率r和临界值，计算出未间伐经营条件下的空间结构函数的值为初始解ans。根据竞争指数的比较，建立间伐经营方案链表，在某一温度下随机选取一个索引的间伐经营方案的空间结构函数作为新解tmp，当∆=tmp-ans<0时，接受tmp；当∆=tmp-ans>0时，根据数据量设置该温度阶段的循环次数K=10。当循环终止时，按照exp((tmp-ans)/T)的概率接受次优解tmp并且将温度下降为T=T*r，当温度降到临界值时，得到间伐经营方案的解趋向于全局最优解。

采用了张雄清等人^[17]的基于Hegyi竞争指数的可变生长率的单木生长方程作为林分的生长驱动方程，结合基于模拟退火算法的最优间伐经营方案选择实现林分的生长与经营的动态交互，并在动态的变化中满足经营者的需求。经营动态可视化模拟的过程如下：

①初始化林分场景。通过读取Excel数据表，设杉木根部坐标为$ {\rm{(x,y,z)}}$，初始单木胸径为D₀，树高为H₀，冠幅为CW₀，活枝下高为UBH₀，则初始时林木的位置Tree_Position=${\rm{(x,y,z)}}$，X-Z平面缩放系数Tree_XZscale=D₀/D(D为原始模型的胸径大小，如图1所示)，Y方向上的缩放系数为Tree_Yscale=H₀/H(H为原始模型的树高大小)。

Figure 1.  Initial Chinese fir model &stump model

②林分经营判断与最优经营方案选择。基于空间结构参数分析森林经营紧迫性，根据分析结果和经营者的实际考虑确定是否进行间伐。一旦确定进行间伐，首先通过竞争指数的比较初步筛选出预间伐木，建立间伐经营方案链表，然后基于模拟退火算法确定最优间伐经营方案，被确定为间伐木对象其间伐属性(Cut)设置为1，Tree_Position=$ {\rm{(x,y,z)}}, {\rm{x,y,z}}$设置为无限远，并用Unity3D中Asset Store下载的树桩模型（如图1）替换树木原来的位置。

③林分生长动态变化。基于Unity3D关键帧技术，将树木的生长以一定的间隔化分为若干阶段，每一阶段用一个关键帧表示。在基于竞争指数的可变生长率的生长方程的模拟下，计算出第i关键帧的杉木胸径设为D_i，树高设为H_i，冠幅设为CW_i，活枝下高设为UBH_i，则此刻林木的位置Tree_Position=$ {\rm{(x,y,z)}}$，X-Z平面缩放系数Tree_XZscale=D_i/D_i-1(i>1)，Y方向上的缩放系数为Tree_Yscale=H_i/H_i-1(i>1)。

④以健康林分结构特征和目标林分胸径为经营目标作为停止条件，将步骤②和步骤③进行循环往复，直到满足经营目标为止。

根据惠刚盈等人提出的健康林分的判断标准^[11]，将纯林林分是否间伐划分为4种情况：不经营、适当经营、需要经营、必须经营。将攸县黄丰桥国有林场5号样地数据导入，初始化渲染场景，分析样地的空间结构特征如表2，经营紧迫性为适当经营，表明有其中一个指标不满足条件，由表2可以看出，平均角尺度不在[0.475~0.517]的健康林分区间。

通过改进张雄清等人的基于Hegyi竞争指数的可变生长率的单木生长方程，进行胸径生长拟合，参数方程如公式（2）：

$ {\rm{D}}\left( {i,t + 1} \right)$：表示第i棵树在t+1年的胸径；$ D\left( {i,t} \right)$：表示第i棵树在t年的胸径；At：表示林分的年龄；CIi：表示第i棵树的Hegyi-Ci的竞争指数；B1、B2、B3为参数。

对前4块样地的连年调查数据进行处理，剔除波动较大的异常值数据，然后通过ForStat2.2进行非线性拟合，拟合胸径的生长方程结果如表3所示，参数B1、B2、B3取值的95%的置信区间上下限均没有跨越0，表明参数值可用；表4数据显示B1、B2、B3三个参数的相关性较低，说明函数构造合理，故拟合出的胸径生长的方程。在得出胸径的生长方程后，通过线性回归分析可以得出胸径与树高、冠幅以及树高与枝下高之间存在明显的线性关系，如图2所示，从图中可以看出D-H、D-CW、H-UBH的R²分别是0.790 8、0.541 8、0.578 6，表明胸径树高相关性较高，胸径冠幅和树高枝下高相关性较差一些。

Figure 2.  Linear regression results of three tree factors (DBH-HDBH-CW, H-UBH)

选择5号样地的连年数据对上述单木生长方程结果进行验证，由图3可知胸径与树高估计值与实测值偏差较小，枝下高与冠幅的估计值与实测值偏差较大，这一趋势与上面拟合的单木生长方程R²的大小比较结果一致，虽然枝下高与冠幅的拟合结果与真实值有一定的差距，但总体上具有一定的相关性，因此基本满足生长拟合的要求。

Figure 3.  Comparison of fitting value and measured value of four tree factor

根据间伐量不得大于生长量的原则，在竞争指数的比较中对采伐木进行初步的筛选，本研究根据样地实际情况，筛选出竞争指数较大的前10%作为预间伐木，然后对初步筛选后的间伐木进行组合排列，制定出所有经营方案的集合，通过模拟退火算法对所有间伐经营方案进行最优解的计算，得出最优采伐方案，同时对最优解与未间伐下的空间函数进行比较，以确定是否需要实施经营，由于初始林分（16 a）时林分经营紧迫性分析为适当经营（如图4左），所以对该样地进行连续5 a的生长模拟，再对生长后（21 a）的林分重新进行经营紧迫性分析，分析结果为需要经营（如图4右），因此对其进行经营模拟。图5左为通过竞争指数初步筛选的间伐木结果，从图中可以看出初步筛选了23株预采伐木，图5右为基于模拟退火算法判断的最优间伐经营方案选择结果，实际间伐了4株树木。

Figure 4.  Analysis of management urgency before (left) and after (right) growth

Figure 5.  Judgment of thinning (left)and results of thinning execution(right)

经营动态可视化模拟过程可以抽象成林分状态场景、林分结构与经营紧迫性分析、林分生长、间伐经营方案和经营目标5个模块，它们之间相互连接与迭代构成了林分动态生长与经营过程。如图6所示，林分状态场景用于模拟林分不同时期的三维渲染；林分结构与经营紧迫性分析用于分析不同时期林分结构的状况和经营紧迫性的等级，主要侧重于数字化的表达，从而判断林分是否生长还是间伐；林分生长通过生长率的变化模拟出下一年的生长量；间伐经营方案则是在分析结构与胸径的基础上以间伐树木控制林分的状态；在实施生长或者间伐后，判断经营目标，重新回到林分状态，如果满足经营要求停止活动，否则继续以上步骤。

Figure 6.  Dynamic process of management

本研究以林分平均胸径达到30 cm和健康的林分结构为经营目标，对5号样地进行经营动态可视化模拟，模拟结果如图7所示。5号样地初始林分场景时（如表5所示），林分平均年龄为16 a，平均胸径为20.21 cm，平均角尺度为0.64，即分布状态为团状分布，经营紧迫性为适当经营，所以对5号样地进行采取生长模拟。本研究将间伐间隔期设定为5 a，模拟出5号样地21年生的林分状况，此时林分平均胸径为23.79 cm，平均角尺度没有变化，因为没有对林分进行间伐模拟，林分树木相对位置没有发生变化，林分拥挤度为0.81，表明林分平均冠幅变大，林分逐渐拥挤，经营紧迫性为需要经营，对林分进行间伐模拟，模拟结果为图7中上所示，间伐了其中的4株树木，间伐后重新对林分进行分析，结果表明林分经营紧迫性仍为需要经营，但间伐后如图8所示密集度分布向右偏移，表明拥挤程度有所下降，大小比数向左偏移，表明垂直竞争下降，所以经营仍然有效。由于平均胸径未达到目标，因此继续对林分进行生长模拟。在连续的生长模拟与间伐模拟后，当林分平均年龄模拟到31 a时，林分平均胸径达到30.02 cm，林分经营紧迫性为必须经营，表明仍需要对林分进行经营模拟，对林分进行间伐模拟后，经营紧迫性等级下降一级，间伐后平均胸径为30.10 cm，故结束此次林分模拟。

Figure 7.  Management dynamic visualization simulation process

Figure 8.  Distribution of structural parameter grades before and after management(left 21a, right26a)

本研究以湖南省攸县黄丰桥国有林场5块长期固定样地为数据源，通过林分结构参数，建立了空间结构函数，并采用模拟退火算法进行间伐经营方案的最优选择，用以指导经营采伐；同时，研建了基于竞争指数的可变生长率的生长方程，得到了胸径与树高、胸径与冠幅、树高与枝下高之间具有相关性性，且胸径与树高之间相关性最强；并基于Unity3D三维渲染引擎，模拟了林分动态经营采伐过程，效果形象逼真，运行流畅，但也存在一些不足之处。

1）本研究模拟的为同龄纯林，未涉及混交异龄的经营，因此空间参数选择较少，对于混交度、开敞度等其他空间结构参数等没有引入，因此经营紧迫性划分等级较为粗放，另外没有考虑补植经营措施的影响，因此需要在以后的研究中进行深入。

2）基于模拟退火算法的最优间伐经营方案选择与模拟退火温度的上下限以及下降速率有关，本研究考虑了计算机的性能，适当的降低了温度的上下限和提高了下降速率。另外，考虑到样地林木数量较少，故在竞争指数初步筛选间伐木时，仅选择了竞争指数较大的前10%的树木作为预间伐木，因此可能会出现间伐后的效果不是太明显，间伐结果不一定是最优经营解的情况。

3）对于可视化模拟大面积动态经营的能力还有所欠缺，主要包括两个方面的原因：一方面大面积树木的渲染需要消耗大量GPU性能，虽然可以通过层次细节（LOD）和射线碰撞检测进行缓解，但是问题依旧存在；另一方面前期的空间结构参数的计算以及组合经营方案集的遍历也会随着树木数量的上升，出现等待的情况，在以后的研究中仍需要研究有效的算法来进行组织。

模拟退火算法是解决大量组合排列问题的有效方法之一，通过合理的构建目标函数，对于解决森林经营中存在的大量组合排列问题也同样具有很好的效果。本研究通过可视化的手段形象地描述了模拟退火算法在森林经营中应用的可行性，生动地展示了在多目标需求下的森林生长与经营的整个过程，对于进一步提高经营水平和经营的精细化程度具有一定的作用。