杉木三维模型各方向枝下高分布研究

崔泽宇，张怀清*，左袁青，杨廷栋，刘 洋，张 京，王林龙,3

(1.中国林业科学研究院资源信息研究所，北京 100091；2.国家林业和草原局森林经营与生长模拟重点实验室，北京 100091；3.中国林业科学研究院林业科技信息研究所，北京 100091)

在林木三维多态性建模过程中，林木的各个关键参数(树高、胸径、冠幅、枝下高等)是实现多态性的重要因素。其中，枝下高是表达树冠形态特征的重要参数。枝下高是指林木树冠第一活枝距地面的高度[1]，是反映林木生长活力、树冠形态特征的重要指标[2-3]。枝下高与空间竞争强度存在关联，因此也是反映林分内竞争水平的重要指标[4]。在传统林学研究中，枝下高具有唯一性，即每株树只有一个枝下高数值，且不关注枝下高所处方向。但在林分环境中，由于林木受到各个方向的空间竞争强度不同，各方向枝下高存在一定差异，按照传统林学研究采集的数据，不利于林木形态结构的多态性表现；而林木处于生长的过程，已结束的外业调查时间段内林木形态结构参数不可再次获取。因此利用传统林学调查数据进行三维建模，存在模拟效果不真实的问题；而对样地再次测量，则存在无法获得林木之前状态的数据，也无法满足模拟需求，并会加大外业工作强度的问题。

以传统林学研究为出发点的枝下高研究较多，Ritchie等[5]利用Logistic模型以树高、胸径、树冠竞争因子、林分断面积等作为变量拟合了14个树种的枝下高模型；Rijal等[6]以林木尺度与竞争因子为模型变量，利用Logistic模型构建了北美阿卡迪地区13个树种的枝下高模型。李想等[7]以联立方程组的方式将林分断面积与优势高作为拟合变量，加入枝下高模型；而混合效应模型的提出进一步提高了枝下高模型的预测精度[8-11]，但并未将枝下高分布与林木树冠形态相结合考虑。在以林木三维模型构建为目标的研究中，马载阳等[12]分别量测林木东西南北各方向活枝下高，计算对应方向垂直空间结构参数以拟合枝下高模型，考虑到了林木不同方向枝下高分布存在差异的问题，但是该模型需要林木年龄数据；李思佳等[13]通过样本库进行枝下高的预测；朱念福等[14]以结构单元为单位，计算中心木水平空间结构参数与垂直空间结构参数，在不考虑林木年龄的前提下，拟合杉木枝下高。通过对已有研究分析，将空间结构参数引入枝下高模型，使得模型与空间方向间建立关系。

杉木(Cunninghamialanceolata)为单轴分枝针叶树种，一级分枝方位角在东南西北方向内近似均匀分布[15]，各方向枝下高的差异对树冠形态影响较大，因此本研究以杉木为研究对象，以朱念福等[14]的枝下高模型为基础，分析空间结构与枝下高的关系，并基于传统林学调查数据，探究林木三维模型各个方向枝下高位置分枝模型分布差异的问题，以期提高已有数据利用率，降低外业工作强度，增强林木三维模型的多态性表现。

1 材料与方法

1.1 试验区概况

试验区位于江西省新余市分宜县亚热带林业实验中心山下林场的杉木人工林区，该林场位于大岗山林区，海拔85～300 m，属低山丘陵地形；年均气温17.5 ℃左右，年降雨量为1 100～1 700 mm，属亚热带季风湿润性气候。林场面积1 790.4 hm2，包括科研实验林246.2 hm2，杉木、木荷、松类等示范样板林755.3 hm2，活立木蓄积量共11.5万m3，森林覆盖率高达95.6%。

1.2 数据获取

考虑到不同年龄和密度的林分会存在树高差异以及不同的空间竞争强度，因此以试验区内不同林分密度、林龄12～30 a的杉木林分为调查范围，在该调查范围内，随机圈定8块20 m×20 m临时样地。样地分布的地形较为均匀，坡度25°～45°，且所有样地分布未出现跨越坡面情况。对样地内239株杉木，按照传统林学研究调查方式进行每木检尺。使用卷尺测量单木胸径(DBH)，使用激光测高仪测量单木树高(H)、冠高以及每株树第一活枝距离地面高度(Hb)，通过塔尺测量单木东西南北4个方向冠幅(Cw)长度，利用全站仪获取样地内每株林木相对位置。样地基本信息见表1。

表1 样地杉木基本信息

2 杉木枝下高三维模型的构建

2.1 枝下高模型与模型变量选择

已有研究构建了枝下高与诸多因子间的关系，从不同角度模拟不同条件下枝下高的情况，所构建模型的形式也呈现多样化。常见的模型形式有指数模型、Logistic模型[5-6]等。常见的枝下高模型变量有树高[16]、胸径等单木形态结构参数，也有研究将冠幅竞争因子、断面积和、气象因子[17]、地位指数[18-19]、空间结构[13]等加入模型。本研究是为实现基于传统林学研究调查数据，对林木三维模型东西南北各方向枝下高进行差异化表现。因此，对于模型变量的选择，除了树高与胸径等单木测量值，需要加入与方向相关的因子，所以选择空间结构作为模型的剩余变量。

朱念福等[14]构建的枝下高模型中，变量包含树高、胸径以及垂直空间结构参数，根据本研究需求选择该枝下高模型作为本研究的理论模型，模型形式为：

Hb=H/(1+ea+b×DBH+c×PV)。

(1)

式中：Hb为枝下高，H为树高，DBH为胸径，PV为垂直空间结构参数，a、b、c为变量参数。

朱念福等[14]认为，水平空间结构的参数(PH)也可提高模型拟合精度，但是提高程度有限，在最终模型中并未加入水平空间结构参数作为变量。而笔者为了研究枝下高的方向分布关系，将与方向相关的垂直空间结构参数与水平空间结构参数都选入模型变量，因此本研究所用枝下高模型为：

Hb=H/(1+ea+b×DBH+c×PH+d×PV)。

(2)

式中，d为变量参数。

2.2 空间结构因子构建

惠刚盈等[20]提出的最佳空间结构单元一般以最近4株木作为周围木，与中心木构建空间结构单元，但是分析周围木对中心木各个方向的影响时，简单地利用距离判别周围木并不准确。本研究以此方法为基础，引入空间分析中的常见方法，即缓冲区构建，以更加精准地判别对中心木可以造成影响的周围木。

缓冲区是对点、线、面等地理实体建立周围一定宽度范围的区域，本研究以林木水平位置作为点，以样地平均冠幅作为缓冲区半径构建缓冲区，表现样地林木树冠影响的范围。通过缓冲区范围分析，剔除距离较小但是被其他相邻木树冠遮挡，无法对中心木造成直接影响的周围木；添加距离较大，却可以直接影响中心木的林木作为周围木，进一步提高所构建空间结构与林木形态结构的相关性。在此基础上，分别记录林木东西南北4个方向对中心木有影响的林木编号(图1)，构建新的空间结构单元，并以此计算水平空间结构参数与垂直空间结构参数。以图1为例，树1为构建结构单元的中心木，以最近4株木方法构建结构单元，周围木为3、4、5、6号树，而通过结合缓冲区构建的方法判断，树5受到3、4号树遮挡，无法直接影响树1，因此周围木被确定为2、3、4、6号树。

1)水平空间结构参数是表示周围木对中心木树冠挤压强度的因子[20]，通过树冠相互挤压抑制分枝的生长，而不同大小的树受到挤压造成的影响不同，因此通过周围木与中心木的胸径比(Dj/Di)表现相邻木在不同挤压强度下其树冠受到影响的差异。构建水平空间结构参数(图2a)，其表达式为：

(3)

式中：i表示中心木，j表示周围木；CWi与CWj分别为中心木与周围木相邻方向冠幅，Di与Dj分别为周围木与中心木的胸径，dij为i与j的水平距离；n在计算单一方向PH时为单一方向周围木数量，计算空间结构单元内PH时为i的所有周围木数量。

2)垂直空间结构参数是反映相邻木之间相互遮挡强度的因子，由于林木可能生长在坡地，各树垂直位置的高差会影响相邻木之间的遮挡关系，因此将树高与林木位置相对高度之和作为相对树高，计算垂直空间结构参数(图2b)，表达式为：

(4)

式中：Hjr为周围木j的相对树高，Hir为中心木i的相对树高；n在计算单一方向PV时为单一方向周围木数量，计算空间结构单元内PV时为中心木所有周围木数量。

2.3 林木三维模型枝下高分布

利用Unity3D渲染引擎，对林木的主干、分枝模型按照调查数据进行加载，东西南北4个方向的枝下高分布由空间竞争强度决定。构建枝下高测量值与东西南北4个方向的整体空间结构关系，判断枝下高与空间结构的相关性。

由于树高是影响枝下高的重要因子之一[16]，直接构建枝下高与空间结构的关系会受到树高影响，因此在建立关系前应该消除树高影响。对每株树与最高树间的高差与实测枝下高求和，使调整后实测枝下高处于统一树高的条件下。调整公式如下：

(5)

式中：ΔH为样地内每株树与最高树的高差，Hmax为样地最大树高，Hi为每株树高，Hba为调整后实测枝下高，Hb为实测枝下高。

通过调整后实测枝下高与空间结构参数的相关关系，进一步判别实测枝下高在三维模型中的分布方向，并计算所有方向枝下高估计值。计算结果分两种情况:

1)当存在除被判别方向外，还有其他方向估计值低于实测枝下高时，考虑估计误差的存在，将实测值与实测枝下高方向估计值的差值作为修正系数，调整剩余方向枝下高估计值，实现在林木三维模型中反映实测数据，并可以表现不同方向枝下高分布的差异。

(6)

式中：Hbe为实测枝下高方向估计值，Hbm为实测枝下高，Hbp为其他方向枝下高估计值，ΔHb为估计与实测枝下高差，Hbpa为其他方向调整后枝下高。

2)当其他方向估计值都大于实测枝下高时，东西南北4个方向枝下高的估计值不做调整，可直接作为相应方向枝下高参数建模。

2.4 林木三维模型构建

主干、分枝模型加载的方法具有对于分枝分布调控便捷的优点。根据调查数据，求解林木冠形曲线，计算公式为：

(7)

式中：hup、hdown分别为所构建冠形曲线方向冠高上、下部分树冠宽度x对应的高度；Hc为冠高；α1、α2为冠形指数；CW为冠幅。

根据冠形曲线公式所做林木冠形曲线如图3所示。

不同的方向枝下高分布不同，直接影响到冠形曲线形态，进一步影响加载分枝模型的起始位置，最终影响林木三维模型的多态性表现。各分部位模型库由主干与分枝模型构成，将分枝模型按照林木形态结构特征与测量数据加载到主干相应位置，实现林木三维模型的构建。

3 结果与分析

3.1 模型参数拟合结果

对调查的239株杉木数据，按照4∶1的比例随机抽取190条数据，利用SPSS 20.2对以树高、胸径以及每株树总的水平空间结构参数与垂直空间结构参数为变量的选定模型变量参数进行拟合得：

Hb=H/(1+e0.349-0.005 DBH-0.008 PH-0.249 PV)。

(8)

对剩余49条检验组数据，按照模型估算枝下高，进行模型配对样本t检验。检验结果：实测值平均数为8.424，预测值平均数为8.473，t检验值为-0.322，P值为0.749。由检验结果可见，枝下高预测值略大于实测值且P>0.05，表明二者之间无显著性差异，数据沿对角线方向分布，残差满足正态性假设，所选模型可用于本实验区数据。

3.2 杉木三维模型枝下高分布判别结果

将实测枝下高与中心木整体的水平空间结构参数及垂直空间结构参数分别进行线性拟合，获得实测枝下高与两个空间结构参数的相关性，分别为-0.182与-0.030。将调整后实测枝下高与水平空间结构参数及垂直空间结构参数分别进行线性拟合，获得调整后实测枝下高与两个空间结构参数的相关性分别为0.410与0.782。

由结果可见，实测枝下高与空间结构参数相关性不显著，且负相关趋势，在调整之后与水平空间结构参数相关性为0.410，与垂直空间结构参数相关性为0.782，枝下高与二者都呈现较强正相关性。由此可知二者与枝下高存在密切关联。在水平方向，周围木对中心木挤压强度越大，水平空间结构参数越大，中心木分枝枯死掉落情况增强；垂直方向，周围木相对高度越高，中心木受到遮挡越多，垂直空间结构参数越大，分枝自然整枝情况增加。枝下高与空间竞争强度呈正相关，而在调查中，实测枝下高为最低分枝高度，因此实测枝下高应分布在东西南北4个方向中空间竞争强度最小的方向。

将水平空间结构参数与垂直空间结构参数以相关系数作为权重，计算空间竞争强度：

PCi=0.41PHi+0.782PVi。

(9)

式中：PCi为林木i单方向空间竞争强度，PHi为对应方向水平空间结构参数，PVi为对应方向垂直空间结构参数。

对每株杉木东西南北4个方向空间竞争强度分别进行计算，以空间竞争强度最小的方向作为实测枝下高分布方向。利用实测枝下高与空间竞争强度最小方向的空间结构参数再次拟合模型，求解枝下高模型参数，结果为：

Hb=H/(1+e0.047+0.003 DBH-0.044 PH-0.206 PV)

(10)

利用检验数据，进行模型配对样本t检验，则实测值平均数为8.424，预测值平均数为8.248，t检验值为1.814，P值为0.076。由检验结果可见，枝下高预测值略小于实测值且P>0.05，表明二者之间无显著性差异，数据沿对角线方向分布，残差满足正态性假设，模型可用于枝下高预测。

3.3 杉木三维模型验证

利用unity3D引擎，通过实测数据判别实测枝下高在东西南北4个方向中的分布方向，并计算剩余方向枝下高，根据以上数据加载已构建的杉木主干、分枝模型。以样地1中编号为24的杉木为例，通过空间结构构建方法判别，24号杉木的周围木有15、16、23、25、28、29号树，24号树枝下高测量值为8.5 m。通过周围木判别以及空间结构参数计算，中心木与周围木空间结构关系见表2。

表2 24号杉木空间结构关系

通过表2数据分析，24号杉木北侧空间竞争强度0.232为4个方向中最小，故在24号杉木三维模型构建时将实测枝下高分配到模型北侧。4个方向枝下高估计值通过公式(10)将对应方向水平空间结构参数与垂直空间结构参数代入计算，最终计算值为：东侧枝下高8.9 m、西侧8.6 m、南侧8.7 m、北侧8.3 m。可见除北侧外，其余方向枝下高估计值均大于实测枝下高，不需要调整枝下高估计值。因此以东侧8.9 m、西侧8.6 m、南侧8.7 m、北侧8.5 m作为建模数据代入，构建三维模型。

根据以上数据加载各部分模型，所构建以24号杉木为中心木的结构单元内7株杉木的三维模型如图4所示。所构建中心木，即24号杉木三维模型东西方向如图5a所示，南北方向如图5b所示，整体模型如图5c所示，未考虑东西南北各方向枝下高分布差异所构建的24号杉木三维模型如图5d所示。

4 讨 论

对比两种方法所构建的杉木三维模型可见，与目前仅输入一个枝下高数值的参数化建模方法不同，本研究充分考虑了空间竞争影响下杉木不同方向枝下高分布的差异；与需要测量不同方向枝下高的方法相比，本研究简化了杉木各个方向枝下高数值的获取方式；与传统林学研究中有关枝下高模型的研究相比，本研究的目的是对杉木东西南北4个方向最低的活枝进行三维表述，综合考虑了模型的复杂程度，以及模型变量与方向的相关性，构建了以林木属性与单一方向空间结构为变量的枝下高模型。最终凭借分枝、主干模型加载方法便于对不同方向分枝分布进行调节的优势，使得杉木三维模型树冠部分各方向差异能直观表现，更加符合林分环境中树冠形态结构的特征。研究面向林木三维模型构建，此方法在保证实测数据使用的基础上，根据林学规律对各方向的枝下高进行估计，可实现以实测数据为限制的林木三维模型各方向枝下高差异性表现，提高了基于传统林学调查数据构建林木三维模型的多态性表现效果，同时也适用于林分密度较大、枝下高量测困难的情况，还可通过采集相对精确的数据估计量测困难、量测误差较大的数据，减少外业工作量。根据针叶树形态结构特征，除了杉木，本研究方法应适用于其他针叶树种，但需要对不同树种拟合构建新的枝下高模型。

本研究在传统空间结构构建的基础上，引入缓冲区分析方法，相对于按照相对距离判别空间结构单元周围木的方法，考虑相邻木之间的影响作用，对直接影响林木形态结构的相邻木进行筛选，进一步提高所构建空间结构单元中心木与周围木的相关性，构建的空间结构单元也更符合林木实际生长情况。但也有不足之处：①根据模型变量选用已有的枝下高模型，并未重新构建模型以及与其他形式模型进行对比，下一步研究可以参考其他有关枝下高模型的研究，筛选更多环境因子作为变量，对比不同形式模型，构建模拟精度更高的单方向枝下高模型；②在构建空间结构单元时可以考虑结合无人机影像等更加精确判别对中心木形态结构可以直接造成影响的周围木；③可以同时在数据采集阶段测量更多环境条件中林木参数，使得模拟结果适用于更广泛的林分条件。

以江西省分宜县山下林场为例，通过缓冲区构建的方法调整空间结构单元构建过程，分析空间竞争强度与枝下高分布方向的关系。结果表明：①将空间竞争强度最小方向判别为实测枝下高分布方向，利用已有枝下高模型；进一步估计剩余方向枝下高，可较好地模拟杉木不同方向枝下高分布；②结合unity3D渲染引擎加载杉木主干、分枝模型，实现基于实测数据与估计数值的林木三维模型构建，模拟结果可反映不同方向空间结构对枝下高的影响，体现林木形态结构多态性，模拟效果逼真，对林木多态性建模及林业可视化发展有积极作用。