-
目前国内外已有不少学者对林分生长可视化模拟进行了相关研究，主要侧重于利用结构功能模型对不同发育阶段的单木细节形态进行绘制[1-3]，或者通过建立描述林分生长的数学函数，而后基于三维引擎或者现有软件进行林分生长三维模拟[4-7]。结构功能模型难以与林分空间结构建立联系，生长函数中仅表达了林分空间结构对生长单方面的影响，未考虑林分空间结构与生长的交互过程。而林分空间结构作为林分中与林木空间位置相关的空间特征[8]，是林木在林地上的空间分布格局及其属性在空间上的排列方式[9]，决定了林木之间的竞争势及空间生态位，影响着林木的生长和林分的稳定性[10]，反过来林分生长又改变着林分空间结构，有学者使用可变生长率法建立单木直径年生长量模型和年存活率模型，考虑了单木在生长期间，林分因子及单木因子的变化引起的单木年存活率及直径连年生长量的变化[11]，其中所蕴含的迭代思想算是交互的一种形式，但对于交互的表达过于抽象。谢小华等[12]基于本体规则推理和案例推理相结合的方法表达林木生长与其周围环境之间的复杂关系，通过建立本体知识库开展林木种植适宜性判断以及林木立地质量评价，进而对林木生长进行模拟，实现了林木种植智能化决策，弥补了数学模型表达抽象这一缺陷。本研究以杉木(Cunninghamia lanceolata (Lamb.) Hook.)林分为研究对象，通过构建杉木生长样本库，研究一种考虑林分空间结构与生长之间交互过程的生长模拟方法，同时结合三维可视化技术，实现林分生长动态可视化模拟。
HTML
-
试验区位于湖南省攸县黄丰桥国有林场，介于113°04′～113°43′ E、26°43′～27°06′ N之间，最高海拔1 270 m，最低海拔115 m，属亚热带季风湿润气候区，年平均气温17.8℃，年降水量1 410.8 mm，森林覆盖率达90.07%，是湖南省重点速生丰产林基地，杉木大径材推广示范基地。
-
利用全站仪测量杉木根部相对坐标(X, Y)，以常规外业方式调查胸径、树高、冠幅、活枝下高、生长状态等，按1年时间间隔采集2012—2017年连续调查数据。从所有样地中选择地位指数一致(18，参考地位指数表见文献[13])的杉木人工林样地6块，其核心区(距边界5 m)属性分布如表 1所示。
样地序号
Sample plots No.面积
Area/(m×m)年龄跨度
Age span/a株数
Stem number胸径分布D
distribution/cm树高分布H
distribution/m活枝下高分布UBH
distribution/m冠幅分布CW
distribution/mA 60×60 17～22 534 7.1～16.2～28.6 4.3～10.6～18.0 1.7～5.0～9.3 0.6～2.7～3.8 B 50×80 23～28 362 13.3～22.6～33.4 8.4～15.4～20.2 1.8～8.21～2.2 1.3～2.9～4.2 C 50×50 11～16 309 5.9～14.2～21.2 5.3～9.5～13.5 1.9～4.3～7.5 0.8～2.7～4.0 D 50×80 10～15 955 3.2～11.6～22.4 3.2～8.4～13.7 1.3～4.4～8.0 0.3～2.2～3.8 E 50×50 16～21 230 10.5～20.0～29.8 6.7～12.8～16.8 2.6～6.6～9.2 1.6～3.0～4.5 F 30×30 13～18 120 6.2～16.5～23.2 5.6～11.0～13.8 1.8～6.3～9.5 1.8～3.4～5.6 注：属性分布中，a～b～c表示最小值～均值～最大值。
Note: In the attribute distribution, a～b～c represents the minimum value～average value～maximum value.Table 1. Summary statistics of sample plots
1.1. 试验区概况
1.2. 数据采集
-
对于林木生长而言，在水平空间层面，林木生长主要受生存空间的影响，生存空间可表达为对象木i与周围竞争邻体j之间的冠幅挤压[14]。据此建立能表达周围竞争邻体与对象木之间挤压程度的水平空间结构参数PHi，计算公式如下：
式中：CWi为对象木i的冠幅；CWj为竞争邻体j的冠幅；dij为i与j的水平距离。
林木之间的竞争是相互的，同样的冠幅挤压对于不同大小的林木所造成的影响也必然有所差异，故利用(Dj-Di)表达二者之间的相对大小差异，对冠幅挤压程度进行加权，得到最终水平空间结构参数PHi如下所示：
式中：Di为对象木i的胸径；Dj为竞争邻体j的胸径。当PHi越大，对象木所受到的冠幅挤压越强，当PHi为负值时，对象木处于竞争优势地位。对于一个空间结构单元而言，对象木受到n株竞争邻体的挤压，n值由象限补树法确定[15]，则平均水平空间结构参数A_PHi为：
-
林木生长在垂直方向上取决于光照的分布，光照主要受竞争邻体上方遮盖的影响。对于对象木i来讲，如果邻体j的高度(Hj)比其高度(Hi)大(如图 1c)，则上方的遮盖由高差(ΔH=Hj-Hi)引起，其影响还与对象木和邻体之间的水平距离dij有关，这种关系可用交角α表示[16]。对象木的高度等于邻体的高度时(如图 1b)，二者在垂直方向上互无遮盖，此时α=0°。对象木的高度高于邻体的高度时(图 1a)，对象木对邻体造成遮盖，此时令α为负。可见α不仅能表达对象木与邻体之间的竞争，而且能对二者的大小差异进行表示。
Figure 1. Vertical space structure parameter diagram
基于上述分析，垂直空间结构参数PVi可表达为：
式中：Hi为对象木i的树高；Hj为竞争邻体j的树高；dij为对象木i与竞争邻体j的水平距离。当PVi越大，对象木所受到的光照遮盖越强，当PVi为负值时，对象木处于竞争优势地位。对于对象木i的n株竞争邻体来讲，则其平均垂直空间结构参数为：
-
林分空间结构与生长交互本质上是一个迭代递归的过程，传统的数学生长模型未考虑二者之间的交互过程，本研究采用样本库作为交互媒介解决这一难题。利用2012—2017年连续调查数据构建杉木生长样本库，通过因子相关性分析从林木生长影响因子中选取样本库指标，以连年生长率作为预测下一年生长属性的依据。
待模拟林分中每株杉木都可以看作是一个包含胸径、树高、冠幅、活枝下高、水平、垂直空间结构参数信息的多维空间中的一点(D, H, CW, UBH, PH, PV)，采用最短马氏距离作为待模拟林木与样本库的相似度匹配准则，得到考虑空间结构与生长交互过程的林分生长模拟算法，具体步骤如下：
(1) 根据初始年龄林分中林木的位置(X, Y)及属性数据(D, H, CW, UBH)，计算得到林木基于空间结构单元的空间结构参数(PH, PV)；
(2) 将待模拟林木的属性(A, D, H, CW, UBH)及空间结构参数信息(PH, PV)与样本库中所有林木的特征指标进行相似度匹配(年龄应一致)，得到与待模拟林木马氏距离最短的最相似木，利用最相似木的连年生长率计算待模拟林木下一年的属性(D1, H1, CW1, UBH1)，实现林木在不同空间结构下的不同生长；
(3) 判断是否达到终止年龄，如果否，则返回步骤(1)，计算新的空间结构参数(PH1, PV1)，实现林木在不同生长情况下产生不同的空间结构，如此迭代，直到达到终止年龄；如果是，则结束循环，得到林分生长过程中每株林木的属性、存活状态及空间结构参数变化等信息。
-
目前大部分的林分生长可视化采用树木模型整体放缩的方法，导致生长形态变化与自然真实不符，而采用结构功能模型进行精细绘制会占用大量计算资源，因此本研究基于树木模型剖分和动态组织对林分生长动态进行可视化模拟，将树干与树冠模型剖分并分别控制，打破林分可视化模拟“千树一面”的局限，实现林分中林木的多样化展现。基于Direct3D所构建的杉木三维模型如图 2所示，模型树高H=13.0 m，胸径D=15.0 cm，冠幅CW=4.0 m，活枝下高UBH=4.0 m。采用Mogre三维渲染引擎支持的左手坐标系，树干模型的坐标原点设于根部，而树冠模型的坐标原点设于树顶处。
Figure 2. Schematic diagram of Chinese fir model split
采用节点动画技术实现林分生长动态可视化模拟，将树冠模型节点作为树干模型的子节点，子节点的位置是相对于父节点而言的，且子节点会继承父节点的移动缩放旋转，反过来子节点的变化无法影响父节点。因此林分生长渲染过程如下所示：
(1) 林分初始状态渲染。设杉木根部坐标位置为(X, Y, Z)，初始时杉木胸径为D1，树高为H1，冠幅为CW1，活枝下高为UBH1，则初始时树干关键帧的位置Trunk_Position1=(X, Y, Z)，X-Z平面缩放系数Trunk_Scale_XZ1=D1/D，Y方向上缩放系数Trunk_Scale_Y1=H1/H。初始时树冠关键帧的位置Crown_Position1=(0, H, 0)，X-Z平面缩放系数Crown_Scale_XZ1=(CW1/CW)/Trunk_Scale_XZ1，Y方向上缩放系数Crown_Scale_Y1=((H1-UBH1)/(H-UBH))/Trunk_Scale_Y1；
(2) 下一时刻生长状态渲染。设此时杉木胸径为D2，树高为H2，冠幅为CW2，活枝下高为UBH2，则此刻树干关键帧的位置Trunk_Position2=(X, Y, Z)，X-Z平面缩放系数Trunk_Scale_XZ2=D2/D1，Y方向上缩放系数Trunk_Scale_Y2=H2/H1。此刻树冠关键帧的位置Crown_Position2=(0, H, 0)，X-Z平面缩放系数Crown_Scale_XZ2=(CW2/CW1)/Trunk_Scale_XZ2，Y方向上缩放系数Crown_Scale_Y2=((H2-UBH2)/(H1-UBH1))/Trunk_Scale_Y2；
(3) 如此循环往复，得到每株杉木任一时刻关键帧的状态，从而形成每一株杉木的动画轨迹，进而组成整个林分的生长动画。
2.1. 水平空间结构参数
2.2. 垂直空间结构参数
2.3. 考虑空间结构与生长交互过程的林分生长模拟
2.4. 生长动态可视化模拟
-
为尽可能减少样本库指标之间的数据冗余，对胸径、树高、冠幅、活枝下高的连年生长率与胸径、树高、冠幅、活枝下高、水平、垂直空间结构参数之间进行相关性分析。由于研究对象为杉木同龄林分，所以应该分年龄进行相关性分析，得到某一年龄时林分生长的影响因素。对10—28年林木数据逐年进行因子相关性分析，总体来说，树高、冠幅、活枝下高的连年生长率与各指标的关系较弱，无明显规律，胸径连年生长率与胸径、树高、冠幅、水平垂直空间结构参数之间关系较为紧密(如图 3即为16年杉木林分的相关系数矩阵)，而林木胸径与树高、胸径与冠幅、树高与活枝下高之间的线性关系明显(如图 4所示，利用所有样地数据建模)。因此本研究采用样本库预测胸径生长，采用图 4中的线性模型预测树高、冠幅、活枝下高的生长。
Figure 3. Correlation coefficient matrix (16a)
Figure 4. The forest attributes linear model
从图 3发现，胸径与树高、冠幅之间相关性较高，已有研究也表明三者之间关系密切，意味着三者之间存在较大冗余。考虑到胸径易于测量，精度较高，故剔除树高、冠幅。水平与垂直空间结构参数之间相关系数也较高，但其所代表的是林木所处空间结构单元的两个不同方面，因此保留二者。年龄作为影响林木生长活力的不可或缺因素，必须加以考虑，所以组成胸径生长样本库的指标包括年龄、胸径、水平与垂直空间结构参数及胸径连年生长率。图 5为16年杉木林分的空间结构参数与胸径连年生长率关系示意图。
Figure 5. Spatial structure parameters and growth relationship diagram (16a)
从图 5可以发现，空间结构参数与胸径连年生长率之间呈负相关，随着空间结构参数的增加，所受到的竞争越强，其胸径连年生长率降低。使用A、B、C、D、E样地中边界木以外(边界木只参与计算，而不加入样本库)的杉木生长数据建立样本库，得到的样本库包含杉木8 934株，其中枯死97株。样本库中存活和枯损杉木的胸径、水平和垂直空间结构参数三者的联合分布如图 6所示。
Figure 6. Sample library parameter distribution
从图 6可以发现样本库中存活木和枯损木在胸径、水平、垂直空间结构参数三者联合分布方面差异明显。当PH、PV越大，D越小，数据点越接近于左下方黑色正方形点处时，林木越偏向于枯损，当数据点越接近于右上方绿色棱形点处时，PH、PV越小，D越大，此时林木越偏向于存活。此图表明所创建的水平与垂直空间结构参数符合林木生长现实状况，利用胸径、水平、垂直空间结构参数创建的样本库合乎常理，因此最终样本库的数据结构如表 2所示。
序号
ID字段名
Field name数据类型
Data type长度
Length允许空
Null主键
Primary key默认值
Default字段说明
Field discription1 ID int 4 否(No) 是 0 ID标识列，自动编号ID 2 Age int 4 否(No) 0 杉木年龄Age/a 3 D double 20 否(No) 0 胸径DBH/cm 4 PH double 20 否(No) 0 水平空间结构参数PH 5 PV double 20 否(No) 0 垂直空间结构参数PV 6 D_Rate double 20 否(No) 0 胸径连年生长率dD/dt 注：林木枯损时，其胸径连年生长率为0。
Note: the annual growth rate of DBH is 0 when tree is dead.Table 2. Sample database structure
-
将F样地13—18年实际观测数据作为检验数据，以13年作为起点，模拟18年的状态，利用Spss24.0对F样地观测值与预测值进行配对样本T检验及线性回归，得到检验结果如表 3、4及图 7、8所示。
属性Attribute 胸径D 树高H 冠幅CW 活枝下高UBH 观测值&预测值(Observed & Predicted) 相关性R 0.962 0.746 0.694 0.677 显著性Sig. 0.000 0.000 0.000 0.000 观测值-预测值(Observed-Predicted) 平均值Mean 0.117 0.031 0.054 -0.123 标准差Std. 0.961 1.151 0.461 1.522 均值标准误差S.E. Mean 0.089 0.107 0.043 0.141 自由度n-1 115 115 115 115 T值T value 0.939 0.285 1.207 -0.867 显著性Sig. 0.350 0.776 0.208 0.388 Table 3. Result of T test
属性
Attribute观测值与预测拟合模型
Fitted model of observed & predicted决定系数
R2均方根误差
RMSE自由度
n-1F值
F value显著性
Sig.D OD=1.101 PD- 1.743 0.925 0.919 115 1 399.223 0.000 H OH=0.786 PH+ 3.063 0.556 1.049 115 142.705 0.000 CW OCW=0.551 PCW+ 1.746 0.482 0.364 115 106.041 0.000 UBH OUBH=0.493 PUBH+ 4.345 0.459 1.110 115 96.696 0.000 Table 4. Result of linear regression
Figure 7. Residual distribution chart(followed by DBH, H, CW, UBH)
Figure 8. Linear fit chart
从表 3可以看出，胸径、树高、冠幅、活枝下高观测值与预测值之间呈极显著相关，相关系数依次降低，胸径观测值与预测值相关最为密切，相关系数达0.962。观测值与预测值之差的显著性都大于0.05，表明胸径、树高、冠幅、活枝下高的观测值与预测值之间并不存在显著差异，生长模拟具有统计学意义。胸径、树高、冠幅观测值与预测值之差的平均值为正，活枝下高观测值与预测值之差的平均值为负，说明就总体状态来说，胸径、树高、冠幅预测偏小，活枝下高预测偏高。树高观测值与预测值之差的平均值最接近于0，总体来看，树高预测值最接近观测值，但其数据波动相对较大。冠幅观测值与预测值之差的平均值、标准差、均值标准误均较小，但预测值与观测值之间相关性较弱，胸径观测值与预测值之差的均值标准误仅高于冠幅，且相关性极高。根据图 7的残差分布图看出，胸径、树高的残差分布较为均匀，拟合效果较好，而冠幅、活枝下高的残差呈对角分布，当冠幅、活枝下高预测值较小时，预测值小于实际值，而当冠幅、活枝下高较大时，预测值偏大。因此总体而言基于样本库的胸径预测效果最佳，树高次之，其次是冠幅，活枝下高的预测相对最差。
根据表 4及图 8可以看出，胸径预测值与观测值拟合效果极好，决定系数达0.925，其次为树高、冠幅、活枝下高，这也与上述配对T检验的分析正好吻合。样本库在实现了生长预测的同时，也对枯损进行了模拟，结果如表 5所示。
年龄(Age) 13 14 15 16 17 18 观测木1 Observed tree 1 Ο Ο Ο × × × 观测木2 Observed tree 2 Ο × × × × × 预测木1 Predicted tree 1 Ο Ο × × × × 预测木2 Predicted tree 2 Ο Ο × × × × 预测木3 Predicted tree 3 Ο Ο Ο × × × 预测木4 Predicted tree 4 Ο Ο × × × × 注：Ο表示存活Survive；×代表枯损Mortality。 Table 5. schematic diagram of tree mortality
由表 5可以发现，实际观测结果中杉木枯损2株(观测木1、观测木2)，而预测结果中除此2株(预测木1、预测木2)外，还额外枯损了2株(预测木3、预测木4)，且预测木1、2的枯损时间与观测木1、2有所差异，说明样本库虽然具有预测枯损的能力，但是一方面预测枯损率偏大，另一方面枯损预测时间和过程与实际观测结果有所不同。
-
本研究以Microsoft Visual Studio 2008为平台，使用C#为开发语言，基于Mogre三维渲染引擎，以节点动画(Node Animation)技术实现了F样地杉木林分(120株)的生长动态可视化模拟，其效果见图 9(左侧为初始状态13年，右侧为18年，上为侧视图、下为俯视图)。
Figure 9. Comparison chart of stand growth
以节点动画技术实现的林分生长动态模拟过程流畅自然，基本能保持50帧/秒的帧速。基于样本库实现的考虑林分空间结构与生长交互的生长模拟方法，实时计算得到各个年龄的胸径、树高、冠幅、活枝下高及枯损信息，保证了林分中林木的生长合乎自然规律，从而为研究林分生长规律提供直观形象的可视化手段，图 10即为F样地中一株杉木的空间结构与生长交互过程。
Figure 10. Spatial structure and growth interaction process
3.1. 杉木生长样本库
3.2. 样本库生长枯损预测精度检验
3.3. 杉木林分生长动态可视化及交互过程
-
本文以黄丰桥国有林场6块长期固定样地为数据源，研建基于空间结构单元的水平、垂直空间结构参数，参数与生长量呈显著负相关。在此基础上构建样本库实现考虑空间结构与生长交互的生长模拟，对胸径的模拟效果最好、其次是树高、冠幅、活枝下高。同时基于节点动画技术对林分生长进行动态可视化模拟，效果直观流畅，能满足林业科研和生产实践要求，但也存在一些不足之处。
(1) 基于样本库的树高、冠幅、活枝下高的生长和林木枯损的模拟相对于胸径生长的模拟精度略差，可能受样本库指标选择、模型误差累积、样本库容量、匹配方法的影响。在今后的研究中，应构建更为合理的空间结构参数，其次样本库应耦合立地质量指标，使得样本库具备更大的适用范围，同时须选取更合适的数据降维及精确匹配方法，保证样本库精度和效率之间的均衡，进而提高样本库的模拟预测水平。
(2) 基于模型剖分、动态组织和节点动画技术实现了样地级别的林分生长动态可视化模拟，能满足胸径、树高、冠幅、活枝下高4个指标的生长模拟，但存在着渲染精细程度与计算机消耗资源这一固有矛盾，今后研究中一方面要在数量和质量之间权衡，根据具体研究要求选择对应的模型剖分精细程度，另一方面要研究其他的更为有效的生长动态可视化模拟形式。总而言之，本研究初步提出了一种基于样本库的林分生长动态可视化模拟方法，但其中存在的问题与不足仍需继续深入研究。
DownLoad: