利用Voronoi图评价油松人工林空间结构

孙宇晗，王士博，王润涵，郑小雨，闫 飞

（北京林业大学 林学院 精准林业北京市重点实验室，北京100083）

林分空间结构反映了全林分在空间上的位置特征，决定了林分间各部分相互影响和制约的关系［1－2］。林分空间结构影响林木的竞争程度、空间生态位和全林分生长的稳定性［3－4］，研究林分空间布局能为森林经营决策提供重要的理论支撑［5］。林分空间结构单元作为林分空间上的基本组建单位，是由一株中心木以及其四周邻近木构成的区域［3］。在森林经理学中，通常选取某些参数分析森林空间结构特征，这些参数包括聚集指数、混交度、大小比数和开敞度等，它们反映了林分在竞争、生长空间、树种隔离等不同方面的水平［6－7］。由于空间结构参数的计算依赖于基本单元的划分，因此如何构建最佳空间结构单元一直是研究热点。传统构建森林空间布局单元的方法是固定边长个数为4，固定边长会造成空间结构参数的有偏估量［8］。由荷兰气象学家提出的Voronoi图能够隐含表达林分间林木的所有侧向邻近信息［9－10］。Voronoi又称作泰森多边形［11］，在计算几何、城市规划、气象、地质、地理信息系统、图像处理和机器人路径规划等方面的应用非常广泛［12］。本研究引入Voronoi图评价林分空间结构，Voronoi的特性对于研究林分空间结构单元有3项优势：一是每个泰森多边形是一个小型空间竞争单位，林木在此单位中相互争夺自然资源；二是影响范围。每个结构单元的大小范围不一样，代表每个中心木的影响范围有所区分；三是具有动态针对性。邻近木株数n与Voronoi多边形边长数相同，而每个Voronoi多边形的边数随实际环境的变化发生改变，因此，每株中心木的影响范围都是独立的。此方法更能灵活准确地确定空间结构单元，比传统固定空间结构单元边数的方法更具有针对性，一定程度上克服了固定n＝4时导致的偏差。当前，国内外部分学者已经开展了基于Voronoi图的不同群落空间结构特征基础性研究［13］，但针对油松Pinus tabulaeformis人工林开展的空间结构研究尚未见报道。因此，本研究以北京八达岭油松人工林为研究对象，利用Voronoi图进行森林空间结构单元划分并进行结构参数提取，其结果将为油松人工林生态系统结构优化提供理论依据，对森林群落演替和重构研究［14］、森林空间结构调查信息化精准化具有重要的理论价值。

1 实验区及调查数据

1.1 实验区概况

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

1.2 设置研究区与调查

本项目在具有代表性地区设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。

1.3 坐标转换

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

2 研究方法

2.1 森林空间结构参数选取

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

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

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

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

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

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

2.1.3 大小比数 大小比数为胸径大于中心木的邻近木株数占全部株数的比例［19－20］。计算公式为：

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

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

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

2.2 Delaunay三角网及Voronoi图构建

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

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

图1 步骤二示意图Figure 1 The plot of step 2

图2 森林空间结构制图V1.0软件生成Voronoi图Figure 2 Voronoi diagrams generated by the forest space structure mapping V1.0 software

2.3 研究方法可行性分析及空间结构参数提取

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

3 结果与分析

3.1 空间结构参数的正态检验

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

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

3.2 空间结构参数相关性分析

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

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

3.3 空间结构参数方差分析

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

表1 各指数S—W检验结果Table 1 S—W test results of the structure indexes

表2 各空间结构参数相关性分析Table 2 Spatial structure index correlation analysis results

3.4 森林空间结构分析

3.4.1 林分基本特征 该人工林区经过封山育林和人工栽植油松，其基本情况是森林覆盖率达57%，林木分布不均匀。乔木层共有15种树种，分别为山杏，油松，华山松Pinus armandii，榆树Ulmus pumila，暴马丁香Syringa reticulatavar.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%。

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

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

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

图3 邻近木株数频率图Figure 3 Frequency distribution of adjacent trees

表3 林间树种混交度分布频率Table 3 Frequency distribution of mingling mixing degree of stand and tree species

图4 全林分大小比数分布频率Figure4 Distribution frequencyofneighborhood comparison of whole stand

图5 各树种的平均大小比数Figure 5 Distribution frequency of average neighborhood comparison of each tree species

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

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

图6 不同树种的开敞度Figure 6 Open degree of each tree species

4 结论与讨论

4.1 结论

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

4.2 讨论

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

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

5 参考文献

［1］ 惠刚盈，克劳斯·冯佳多.森林空间结构量化分析方法［M］.北京：中国科学技术出版社，2003：26－28.

［2］ POMMERENING A.Evaluating structural indices by reversing forest structural analysis ［J］.For Ecol Manage,2006,224（3）:266 － 277.

［3］ 赵春燕，李际平，李建军.基于Voronoi图和Delaunay三角网的林分空间结构量化分析［J］.林业科学，2010， 46（6）： 78 － 84.ZHAO Chunyan,LI Jiping,LI Jianjun.Quantitative analysis of forest stand spatial structure based on Voronoi diagram＆ Delaunay triangulated network ［J］.Sci Silv Sin,2010,46（6）:78 － 84.

［4］ 刘帅，吴舒辞，王红，等.基于Voronoi图的林分空间模型及分布格局研究［J］.生态学报，2014，34（6）：1436－1443.LIU Shuai,WU Shuci,WANG Hong,et al.The stand spatial model and pattern based on voronoi diagram ［J］.Acta Ecol Sin,2014,34（6）:1436 － 1443.

［5］ MASON W L,CONNOLLY T,POMMERENING A,et al.Spatial structure of semi-natural and plantation stands of Scots pine （Pinus sylvestrisL.） in northern Scotland ［J］.Forest,2007,80（5）:567 － 586.

［6］ 封尧.基于Voronoi图的杉木生态公益林空间结构研究［D］.长沙：中南林业科技大学，2014.FENG Yao.Study on the Spatial Structure of Cunninghamia Lanceolata Non-commercial Forest based on Voronoi Diagram［D］.Changsha:Central South University of Forestry and Technology,2014.

［7］ 汤孟平.森林空间结构分析与优化经营模型研究［D］.北京：北京林业大学，2003.TANG Mengping.Study on Forest Spatial Structure Analysis and Optimal Management Model［D］.Beijing:Beijing Forestry University,2003.

［8］ 汤孟平，周国模，陈永刚，等.基于Voronoi图的天目山常绿阔叶林混交度［J］.林业科学，2009，45（6）：1－5.TANG Mengping， ZHOU Guomo， CHEN Yonggang,et al.Mingling of evergreen broad-leaved forests in Tianmu mountain based on Voronoi diagram ［J］.Sci Silv Sin,2009,45（6）:1 － 5.

［9］ 陈军，赵仁亮，乔朝飞.基于Voronoi图的GIS空间分析研究［J］.武汉大学学报（信息科学版），2003，28（特刊）：32－37.CHEN Jun,ZHAO Renliang,QIAO Chaofei.Voronoi diagram-based GIS spatial analysis ［J］.Geom Inf Sci Wuhan U-niv,2003,28（spec）:32 － 37.

［10］ 彭仪普，刘文熙.Delaunay三角网与Voronoi图在GIS中的应用研究［J］.测绘工程，2002，11（3）：39－41.PENG Yipu， LIU Wenxi.Study on delaunay triangulation and Voronoi diagram application in GIS ［J］.Eng Survey Map,2002,11（3）:39 － 41.

［11］ 王新生，李全，郭庆胜，等.Voronoi图的扩展、生成及其应用于界定城市空间影响范围［J］.华中师范大学学报（自然科学版）， 2002， 36（1）： 107 － 111.WANG Xinsheng,LI Quan,GUO Qingsheng,et al.The generalization and construction of Voronoi diagram and its application on delimitating city’s affected coverage ［J］.J Cent China Norm Univ Nat Sci,2002,36（1）:107 － 111.

［12］ 秦喜文，张树清，李晓峰，等.基于Voronoi图的丹顶鹤巢址空间格局分析［J］.生态学杂志，2008，27（12）：2118－2122.QIN Xiwen,ZHANG Shuqing,LI Xiaofeng,et al.Spatial pattern of red-crowned crane nest-sites based on Voronoi diagram ［J］.Chin J Ecol,2008,27（12）:2118 － 2122.

［13］ 汤孟平，陈永刚，施拥军，等.基于Voronoi图的群落优势树种种内种间竞争［J］.生态学报，2007，27（11）：4707－4716.TANG Mengping,CHEN Yonggang,SHI Yongjun,et al.Instraspecific and interspecific competition analysis of community dominant plant populations based on Voronoi diagram ［J］.Acta Ecol Sin,2007,27（11）:4707 － 4714.

［14］ 孜来比·买木提名.八达岭地区油松人工林单木竞争模型研究［D］.北京：北京林业大学，2012.Zilaibi Maimutiming.Study on Modeling of Competition Indices of Pinus tabulaeformis Individual Treesat Badaling Area［D］.Beijing:Beijing Forestry University,2012.

［15］ 吴刚，冯宗炜.中国油松林群落特征及生物量的研究［J］.生态学报，1994，14（4）：415－422.WU Gang,FENG Zongwei.Study on the social characteristics and biomass of thePinus tabulaeformisforest systems in China ［J］.Acta Ecol Sin,1994,14（4）:415 － 422.

［16］ 宁杨翠，郑小贤，梁雨，等.北京八达岭油松人工林结构分析［J］.林业科技开发，2008，22（1）：47－49.NING Yangcui,ZHENG Xiaoxian,LIANG Yu,et al.Study on the structure of artificialPinus tabulaeformisstand［J］.China For Sci Technol,2008,22（1）:47 － 49.

［17］ 郝月兰.基于林分空间结构优化的采伐木确定方法研究［D］.北京：中国林业科学研究院，2012.HAO Yuelan.Study on Cutting Tree Determining Method based on Forest Stand Spatial Structure Optimization［D］.Beijing:Chinese Academy of Forestry,2012.

［18］ PIELOU E C.Segregation and symmetry in two-species populations as studied by nearest neighbour relations ［J］.J Ecol,1961,49（2）:255 － 269.

［19］ 惠刚盈，von GADOW K,ALBERT M.一个新的林分空间结构参数：大小比数［J］.林业科学研究，1999，12（1）： 1 － 6.HUI Gangying,von GADOW K,ALBERT M.A new parameter for stand spatial structure:neighbourhood comparison［J］.For Res,1999,12（1）:1 － 6.

［20］ SHI Haijin,ZHANG Lianjun.Local analysis of tree competition and growth ［J］.For Sci,2003,49（6）:938 － 955.

［21］ 林存学.开敞度对次生林林冠下人工栽植红松生长发育的影响研究［D］.哈尔滨：东北林业大学，2013.LIN Cunxue.Effects of Opening Degree on Growth of Korean pine Trees under Canopy of Secondary Forest［D］.Harbin:Northeast Forestry University,2013.

［22］ HESSE M， RODRIGUE J P.The transport geography of logistics and freight distribution ［J］.J Transp Gergr,2004,12（3）:171 － 184.doi:10.1016/j.jtrangco.2003.12.004.

［23］ 胡艳波，惠刚盈.优化林分空间结构的森林经营方法探讨［J］.林业科学研究，2006，19（1）：1－8.HU Yanbo,HUI Gangying.A discussion on forest management method optimizing forest spatial structure ［J］.For Res,2006,19（1）:1 － 8.