Spatial structure of a Pinus tabulaeformis plantation based on the Voronoi diagram

本项目以北京西北延庆县内八达岭森林公园中油松人工林为研究对象, 八达岭林场建场于1958年，位于40°20′89″N，115°59′56″E。最高峰海拔为1 238 m，林木绿化率达96%。八达岭油松人工林的平均年龄为43 a，栽植密度为800~1 200株·hm^-2。林场由于油松林树种单一，生物多样性低，森林结构不合理造成生态功能脆弱。森林抚育管理滞后，造成林木生长缓慢，长势较弱。本研究选取的样地树木位于海拔690~725 m地段，全部处在一个小坡上，林分健康状况一般。

本项目在具有代表性地区设90 m × 90 m的研究区，同时详细分成4个45 m × 45 m大小的调查单元。在单元内对林木挂牌，利用胸径尺和南方NTS-372R型测树全站仪进行每木测量，详实记下各林木编号、树种、树高、胸径、冠幅、调查单元内的相对位置坐标等基本数据。研究区内主要林木类型为针阔混交林，共563株树，主要树种有油松，山杏Armeniaca sibirica，家榆Ulmus pumila，杜梨Pyrus betulifolia和小叶朴Celtis bungeana等^[14-15]。样地平均树高为8.69 m，平均胸径为14.6 cm。

将研究区的基本信息与测得的每木数据按调查单元进行整理、汇总。将调查单元内林木的相对位置坐标统一转化为北京1954坐标系下的真实坐标，再将整理后的数据导入软件中，进行进一步的数据处理与分析。

本研究采用描述空间分布格局的聚集指数（R），表明树种间混交水平的混交度（M），表现林木个体在相互竞争中优势水平的大小比数（U），体现林木个体生长空间大小的开敞度（K）。

聚集指数为中心木与周围最近邻近木的水平距离的均值与期望的平均距离的比例，通常这种计算方法被称作最近邻体分析方法（nearest neighbor analysis，NNA）^[6]。计算公式为：

式（1）中：R为聚集指数；$\overline {{r_A}} $为相邻最近单株距离的平均值；$\overline {{r_E}} $为随机分布下期望的相邻最近单株之间的平均距离；r_i为中心木i与其最近邻近木的距离；N为样地内的林木株数；A为样地面积（m²）；P为样地周长（m）。

实测与预测的偏离程度可利用正态分布进行检验：

式（2）中：$\overline {{r_A}} $，$\overline {{r_E}} $，N，A意义同式（1）；${\sigma _{{r_E}}}$是一个密度为p=$\frac{N}{A}$泊松分布的r_E标准差。

混交度为中心木周围邻近木与中心木是不同树种的个体占总体的比例^[16-18]。计算公式：

式（3）中：M_i为中心木i的混交度；n为邻近木株数；V_ij为离散型变量，若中心木i与第j株邻近木属相同树种，V_ij=0，否则，V_ij=1。

大小比数为胸径大于中心木的邻近木株数占全部株数的比例^[19-20]。计算公式为：

式（4）中：U_i为中心木i的大小比数；n为邻近木株数；若邻近木j的胸径比中心木i小，K_ij＝0，否则，K_ij=1。

开敞度为中心木周围邻近木和中心木的水平间距，与邻近木树高比值的均值^[21]。计算公式为：

式（5）中：K_i为中心木i点的开敞度；n为邻近木株数；D_ij为中心木i与其第j株邻近木的水平间距；H_ij为邻近木j的树高。

建立Voronoi图前首先需要构建Delaunay三角网。本研究使用Bowyer-Watson算法，是一种先粗后细的逐点插入算法。第1步假设一个超级三角形，将其放到链表中。该三角形包含所有的样本点；第2步将离散点逐次插入，并查询出外接圆中包含插入点的三角形，消除受影响三角形的公共边，再把插入点同受影响三角形的顶点相连，完成一个点的链表插入（图 1）；第3步对新三角形LOP算法优化；第4步，循环往复直到所有点插完，完成Delaunay三角网的建立^[22]。

Figure 1.  The plot of step 2

基于Voronoi图与Delaunay三角网的对偶关系，采用间接法生成Voronoi图。利用2.2.1生成的Delaunay三角网，分别对各个三角形3边作垂直平分线，所有垂直平分线相交形成Voronoi图。本研究利用Arc GIS Engine开发平台和C#语言进行二次开发，形成森林空间结构制图V1.0软件。软件生成的Voronoi图如图 2所示。

Figure 2.  Voronoi diagrams generated by the forest space structure mapping V1.0 software

本研究采取频率分布直方图、S—W检验2种方法，分别对基于传统n＝4和Voronoi图方法确定的林分空间结构参数（混交度、大小比数、开敞度）进行正态检验，并通过对2种方法所得数据进行相关性分析和方差分析，判断基于Voronoi图确定空间结构单元的方法是否可行，同时利用Voronoi图进行森林空间结构特征参数的提取和分析。

通过生成各个指数的频率分布直方图可以发现，对于混交度、开敞度2个结构指数来说，基于2种方法得到的数据均为正偏态分布（表 1）。而对于大小比数，基于n＝4方法确定的数据是低于标准正态的正偏态分布，基于Voronoi图的方法获得的数据符合低于标准正态分布的右偏态分布。

本研究样本数量为563株，属小样本量。为进一步验证上述结论正确性，对数据进行S—W检验，通过检验发现所有指数显著性均为P=0.000，所以拒绝原假设，即数据均不符合正态分布，结果如表 1所示。为进一步探证两者关系，进行相关性分析和单因素方差分析（ANOVA），分析2种方法之间是否相关，是否存在显著性差异。

鉴于3种空间结构指数数据均不符合正态分布，因此，本研究选择斯皮尔曼法而非皮尔逊法对n＝4和Voronoi图确定的林分的3个空间结构参数——混交度、大小比数、开敞度两两进行相关性分析。

由表 2可知：基于2种方式得到的各参数之间相关性较高，斯皮尔曼相关性指数均大于0.80，除开敞度外其余参数均高于0.90，高度显著。表明基于n＝4和基于Voronoi图计算的空间结构指数具有高相关性，均为有效参数。以上分析可知基于Voronoi图确定空间结构的方法是正确有效的。

进行单因素ANOVA分析，原始假设利用2种方法确定的各个指数之间无显著差异。分析结果显示：基于2种不同方法计算各空间结构参数的显著性均为P＝0.000＜0.05。所以原假设不成立，2种方法确定的空间结构参数存在显著性差异，两者相互独立。因此，本研究认为Voronoi图明确森林空间布局单元是一种异于基于固定边长n＝4的新方法，且具有显著作用。

该人工林区经过封山育林和人工栽植油松，其基本情况是森林覆盖率达57%，林木分布不均匀。乔木层共有15种树种，分别为山杏，油松，华山松Pinus armandii，榆树Ulmus pumila，暴马丁香Syringa reticulata var. amurensis，小叶朴，刺槐Robinia pseudoacacia，山荆Malus baccata，家榆，桑树Morus alba，多花栒子Cotoneaster multiflorus，杜梨，春榆Ulmus davidiana，蒙古栎Quercus mongolica和椿树Ailanthus altissima。边沿校正后研究区内乔木总数为563株，其中油松占61.3%，山杏18.6%，华山松2.0%，家榆4.0%，暴马丁香2.0%，小叶朴3.0%，杜梨3.0%。

根据Voronoi图确认中心木邻近木株数。如图 3所示：研究区内中心木邻近木株数为3~11株，且大多会集于4~8株，以株数为6最为常见。邻近木株数为3株的情况主要聚集在样地的边沿地带，位于边沿区内的林木不参与空间结构参数计算，并不影响对研究区空间结构的分析。

Figure 3.  Frequency distribution of adjacent trees

通过平均混交度可知：研究区混交程度为中度，油松混交程度为弱度。当以油松为中心木时，图表显示零度混交（M_i=0.00）和弱度混交（M_i∈（0.00，0.25]）的比例较大（表 3）。油松在该林分内占有绝对个体数量优势，以其为中心木聚集生长的空间结构单元数量最多。其他树种，如刺槐、多花栒子、华山松、小叶朴、暴马丁香等的混交程度大多为强度混交（M_i∈（0.50，0.75]）和极强度混交（M_i∈（0.75，1.00]）。说明研究区稳定性较弱，树种间隔离度较低，在今后的森林经营中可以适当补植非油松树种，增大栽植密度，提高树种间隔离程度。

林分之间的竞争度通过林分大小比数来确定。全林大小比数为0.53，表明全林树种优势不明显（图 4和图 5）。油松个体中处于优势、亚优势和中庸的超过该树种整体数目的一半，在研究区中享有优势。优势树种还有小叶朴、桑树、大叶白蜡Fraxinus rhynchophylla，但其株数较少，对全林分结构影响不大。处于劣势或绝对劣势的有椿树、多花栒子、山荆子、春榆、杜梨等，其在空间结构单元中竞争能力较差。通过适当的抚育间伐和补植来降低油松的竞争强度，伐去长势较差而且影响邻近木生长的林木，从而提高其他树种的竞争指数^[23]。

Figure 4.  Distribution frequency of neighborhood comparison of whole stand

Figure 5.  Distribution frequency of average neighborhood comparison of each tree species

将各树种开敞度取算术平均值，得出平均开敞度为0.52，说明研究区内林分生长空间较为充足（图 6）。从树种角度来看，桑树开敞度最高，树高较高，但其株数过少，数据不具有代表性。油松的开敞度为0.49，属于充足状态，其他树种的开敞度均大于0.40。以上数据均表明研究区林分开阔程度高，林木生长空间较为充足。

Figure 6.  Open degree of each tree species

本研究通过聚集指数对研究区的林分空间分布格局进行分析，利用正态检验计算u，分析实测与预测的偏离水平。当|u|＜1.96，林木为随机分布；1.96＜|u|＜2.58，林木为团状或均匀分布；|u|＞2.58且R＜1.00时，林木呈团状分布，当R＞1.00时，林木呈均匀分布。样地校正后的u为0.342，且全林分聚集指数为1.00，表明研究区的林木为随机分布。今后可利用如抚育间伐和补植等森林经营手段，维持研究区内林木随机分布状态。

本研究利用Arc GIS Engine二次开发得到森林空间结构制图软件，通过使用软件生成Delaunay三角网进而得到Voronoi多边形。选取北京市延庆县八达岭国家森林公园的油松人工林为研究对象，将林分中林木个体作为空间平面上的离散点，构建空间结构单元。计算混交度、大小比数、开敞度、聚集指数4个空间结构参数，量化分析人工林的空间结构特性。主要研究结论如下：①本研究利用频率直方图、单样本S—W检验对混交度、大小比数、开敞度等3个结构参数进行正态检验，得到这3种参数均不符合正态分布，其原因与所确定的空间结构单元有关。空间布局单元是会集在空间上一定范围的，且各参数取值离散，导致各项数据不切合正态分布。2种确定林分空间结构单元的方法存在显著性差异且相关性很强，表明基于Voronoi图构建林分空间结构单元是一种新的合理的方法，且相较于常规确定邻近木为4的方法，考虑到了每株中心木的特异性，更加合理全面。②本研究的研究区域为八达岭油松人工林，林分内林木随机分布，林木生长空间较为充足。全林分为中度混交，其中优势树种油松为弱度混交，其他树种则以强度混交和极强度混交为主，表明研究区内种间隔离程度较小，林分稳定性较弱。需要后期的抚育间伐，增大种间隔离度，提高劣势树种的竞争强度；通过建立林分空间结构优化模型，制定相应的林分补植方案，科学地进行森林经营；通过人工更新，保护山杏等树种的林下更新，改善林分内生长态势弱的林木，使其向结构稳定、质量良好的针阔混交林发展，提高林分的综合效益，使林分能得以近自然经营。

本研究只是单纯基于普通Voronoi图确定林分空间结构单元，仍存在一定的偏差。今后将结合林木胸径、树高等因子综合分析，生成加权Voronoi图构建分析林分空间结构，更符合林木实际生长情况。

分析空间结构方面只是选择了混交度、大小比数、开敞度和聚集指数等4个指标，不能完全反映林分空间结构特征。今后研究中需进一步完善。