毛竹向杉木林扩张不同程度林分空间结构遥感量化分析

陆雪婷,曹碧凤,杨樟平,严夏帆,宋贤芬,余坤勇*,刘 健

(1.福建农林大学 林学院,福建 福州 350002;2.3S技术与资源优化利用福建省高等学校重点实验室,福建 福州 350002;3.永安市林业局,福建 永安 366038)

毛竹(Phyllostachysedulis)和杉木(Cunninghamialanceolata)是我国亚热带地区分布最为广泛的树种,两者具有相似的生长环境和立地条件要求,因此常相邻生长[1],并伴有大量的混交界面,这种混交状态主要是由于毛竹自身具有适应能力强、繁衍速度快的生物学特性[2-3]。毛竹扩张主要通过其地下的鞭根系统,不断向周围林分拓展,引起周边林木的枯死或抑制林下幼苗的生长,最终演化为毛竹纯林[4]。林农借助毛竹这种特性,增扩林地边界,引发林权边界纠纷,一定程度上影响了林权制度改革的有效实施和管理。而且,毛竹扩张也会改变林木个体的空间分布,使林分内温度、水分、光照、养分等微环境的空间分布发生变化[5],进而对森林群落结构、生态系统功能等产生负面影响[6-7]。目前关于毛竹扩张的研究,主要集中在土壤理化性质、微生物群落结构和土壤碳储量等方面,马鑫茹等[8]通过测定常绿阔叶林、竹-阔混交林和毛竹林3种林型的土壤理化性质和微生物群落结构,发现土壤微生群落主要的影响因子为有机碳和铵态氮;佟富春等[9]分析了毛竹扩张过程中土壤甲螨的群落结构及其变化规律;池鑫晨[10]的研究表明毛竹扩张会减少原有森林的固碳潜力,但目前对于毛竹扩张过程中林分空间结构的研究较少。

林分空间结构体现了立木及其属性在空间上的分布特征[11],是影响林分的质量和生产力的重要指标之一[12],决定了林木的生长和动态演替的过程[13],对林业的可持续发展具有不可替代的重要意义。林分空间结构一般包含树种混交、林木空间分布格局和大小分化等[14]。毛竹向杉木林扩张过程中会发生时间和空间上的变化,这些变化会造成林分立木空间分布状态的改变。为了全面反映林木与周边木的关系,利用 Voronoi 图确定林木空间结构单元的方法得到了广泛的应用,汤孟平等[15-16]基于Voronoi图计算、V_Hegyi竞争指数和混交度(Mv),并与传统计算方法进行对比分析,得出V_Hegyi竞争指数和混交度(Mv),可更加准确地描述群落优势树种的种内种间和种间相互隔离关系的结论;李际平等[17]基于Voronoi图确定林分空间结构单元,获取混交度、大小比数等林分空间结构参数,分析了杉木生态公益林空间结构特征;刘帅等[18]通过建立Voronoi图模型表示林分空间结构特征,对林木空间信息进行量化分析;宋语涵等[19]通过Voronoi图计算得到空间结构参数,分析了红树林的空间结构特征。但前人研究过程中创建Voronoi图时,林木位置大多是实测数据,难免要消耗大量人力、物力。近年来,无人机遥感凭借高时效、低成本、高分辨率等特性,被应用到各个领域中[20]。如何借助无人机遥感技术,有效、快速地实现毛竹不同扩张程度林分空间结构基础数据的获取,是揭示毛竹扩张演化机制的重要基础,可为毛竹向杉木林扩张不同程度林分空间结构的演变提供一个动态反演的过程,对提升毛竹资源管理水平和林权制度改革成效具有重要的意义。

本研究以毛竹-杉木混交界面为对象,利用无人机遥感影像,进行单木位置获取及毛竹和杉木识别,结合Voronoi图构建林分空间结构单元,量化混交度、角尺度、大小比数等林分空间结构参数,构建林分空间结构评价指数,分析毛竹向杉木林扩张不同程度林木空间结构的变化规律,为揭示毛竹向杉木林扩张演化机制提供理论指导和数据支撑。

1 研究区概况

研究区位为福建省三明市永安市天宝岩国家级自然保护区,位于福建省中部(117°28′3″-117°35′28″E,25°50′51″-26°1′20″N)(图1)。与上坪乡、西洋镇和青水乡相邻,属中、低山地貌,最低海拔580 m,最高海拔1 604.8 m,气候温暖湿润,四季分明,年平均气温15 ℃,土壤类型多为红壤,主要植被类型为毛竹林、针叶林、常绿阔叶林、针阔混交林等,其中毛竹面积为5.85万hm2。

图1 研究区位置Fig.1 Location of the study area

2 试验设计及数据获取

2.1 试验方案

2019年12月在福建省永安市天宝岩自然保护区设置试验样地(图2),从左到右依次是1～4号试验区,4个试验区都沿毛竹向杉木林扩张方向设置条带,共计15条。各条带均由10 m×10 m的小样地构成,总计85个小样地。以每木检尺调查结果为依据,将毛竹扩张不同程度以毛竹占比0%～20%、20%～40%、40%～60%、60%～80%、80%～100%分为5个阶段(图3),便于后续毛竹向杉木林扩张不同程度空间结构变化规律分析。

图2 研究区样地分布Fig.2 Plot distribution map of the study area

图3 毛竹向杉木林扩张程度划分Fig.3 Division of P.edulis expansion to Cunninghamia lanceolata forest

2.2 数据获取及预处理

2.2.1 基础数据采集 2020年1月,利用测高仪、胸径尺、皮尺等测量工具结合目视判定,对样地内所有毛竹的树高、胸径、竹龄、位置等进行测量。位置的测定以调查单元的东北角顶点为坐标原点,以边界接近正东方向的边为x轴,接近正北方向的边界为y轴,来记录每株毛竹的位置。

2.2.2 无人机数据获取及预处理 与基础数据采集时间同步,利用大疆精灵4多光谱版无人机在福建省永安市天宝岩自然保护区进行拍摄。考虑阳光直射效应的影响,选择无风、无云且光照稳定(11:00-14:00)的时间飞行。无人机飞行参数设置为飞行高度120 m,航向重叠为80%,旁向重叠率为75%,飞行速度8 m·s-1。由于无人机拍摄时会受到大气、风力等客观因素影响,存在一定的畸变,因此采用Pix4D Mapper软件对原始无人机影像分别进行信息检测、拼接及几何校正等预处理。通过导入无人机拍摄的单张影像及其自带的位置数据,然后选择白板对每个波段的影像进行辐射校正,将地面目标测得的影像数值转换为图像反射率,从而得到比较真实反映地表反射率的数值,最后利用动态结构算法进行地面特征点的匹配,得到研究区RGB、DSM 以及多光谱影像。

3 研究方法

3.1 单木位置及树种遥感获取

3.1.1 单木位置遥感提取 借助Python3.9版本的Python-OpenCV以及Python-Skimage模块,对DSM进行直方图均衡化,增强灰度值对比,通过最大滤波消除伪树顶,再结合局部最大值法,得到树顶点分布,最后提取树顶点位置及高程信息并统计,实现单木提取。

3.1.2 树种识别 由于毛竹与杉木在无人机可见光谱和植被指数中存在一定的差异性,该差异可为毛竹和杉木的分类识别提供支撑。基于无人机多光谱影像的5个波段,可获取多种波段组合,在描述植被情况上具有很大的优势。采用单波段、双波段模型以及常用的一些植被指数,作为毛竹和杉木识别的指标,结合随机森林机器学习分类器,实现毛竹和杉木的有效识别。

3.2 林分空间结构单元的确定

林分空间结构单元是分析林分空间结构特征的基础,通常由林分内任意一株林木(中心木)和与之最近的n株相邻木所构成[21]。Voronoi图是以对象集合元素的最近属性为依据,将空间划分为许多单元区域,图内任一凸多边形所包围区域的点到该凸多边形对象点的距离比该点到任意其他对象点的距离都小[22]。因此利用Voronoi图所得到的相邻木,可以有效避免2种现象:其一是因所选取的相邻木数量过小而引起的无法兼顾周围所有相邻木的问题;其二是因选取的相邻木数量过大而将实际的非相邻木纳入考虑范围[19]。利用Voronoi图的最近、邻接的特性,进行林分空间结构的计算。

3.3 林分空间结构评价指数的构建

3.3.1 林分空间结构参数的选取 采用角尺度[23]、混交度[24]和大小比数[25]3个林分空间结构指数来分析毛竹-杉木混交林的林分空间结构,同时为避免处于林分边缘树木对林分整体空间结构的影响,结合样地大小,在原样地四周设置1 m宽的缓冲区。对于85个样地,均采用以下的公式来计算林分的各项空间结构指数。

(1)

式中:n为结构单元i中相邻木株数,若相邻木j比中心木i的树高小,则Kij=1,否则Kij=0。大小比数常用来描述个体大小的分化程度。

(2)

式中:Wi表示中心木的角尺度值;n为相邻木的株数,Zij为角尺度计数器;当中心木i与相邻木j所构成的夹角比标准角小时,Zij=1,否则Zij=0。角尺度的标准角大小和相邻木的数量有关,取值为360°/(n+1)。

(3)

式中:Mi表示中心木i的混交度;n为中心木i和相邻木的株数;Vij表示中心木i与相邻木j的树种差异比较结果,当中心木i的第j株相邻木不属于同种时,Vij=1,否则Vij=0。Mi的取值从0～1,代表从弱到强的混交程度,具有明显的生物学意义。

3.3.2 林分空间结构评价指数的构建 参考刘玉平等[26]的方法,选取经济学中应用最广泛、最经典的柯布-道格拉斯生产函数(C-D生产函数)作为函数模型。以林分空间结构综合指数(FSSCI)作为“产出”,混交度(M)、角尺度、大小比作为“投入”生产要素。假设3个参数在演替过程中影响比重一致,鉴于各取值可能为0,参考张君钰等[27]变异系数法的思想确定各个指标的权重。

(4)

式中:WW为角尺度的权重;σW、σM、σU分别为角尺度、大小比数、混交度的标准差。最终建立林分空间结构评价指数的公式如下。

(5)

式中:FSSCIi代表空间结构的评价指数;Mci为混交度;WM为混交度的权重;WU为大小比的权重。当FSSCIi值越大时,说明该单元的空间结构越优。

利用归一化处理,使各单木林分空间结构评价指数的取值范围在[0,1]内。

(6)

式中:xmax、xmin分别表示归一化前所有单木林分空间结构评价指数的最大值和最小值;x'i、xi分别表示归一化后和归一化前第i株中心木的林分空间结构评价指数。

4 结果与分析

4.1 单木及树种遥感获取

4.1.1 单木提取 利用Python编程,采用局部最大值法得到树顶点,实现1～4号研究区的单木提取(图4)。参考前人单木提取精度评价方法[28],采用总体精度、错分误差和漏分误差3个精度评价指标对研究区所提取的单木进行精度评价(表1)。评价结果显示,1～4号样地的总体精度分别为94.41%、91.22%、92.43%、91.01%,总体精度均>90%,满足研究需要。

表1 单木提取精度评价Table 1 Single wood extraction accuracy evaluation

4.1.2 树种识别 以样地实地调查的毛竹和杉木在遥感影像上的定位为圆心,以0.5 m为半径作圆形缓冲区,通过缓冲区来提取无人机影像中各个缓冲区对应的多光谱影像,作为树种识别的分类样本,并以此形成49个波段组合和15个常用的植被指数。该操作的目的是减少因提取光谱影像过小,使冠层光谱异常而导致提取光谱值不具有代表性。鉴于波段组合与植被指数数量过多,计算量太大,通过随机森林重要性排序进行筛选,选择重要性排行前2位的波段组合B1-B2,B1/B2作为主要指标,用SPSS 26.0对毛竹和杉木冠层光谱表达的差异性作独立样本检验。检验结果表明,基于毛竹与杉木树顶点提取的B1-B2、B1/B2组合波段模型中存在着极显著差异。因此,依据毛竹树顶点提取这2类波段组合,可以实现样地林分中毛竹与杉木的分类,得到研究区树种分类结果(图5),总体分类精度为87%(表2),满足研究需要。

图5 树种识别结果Fig.5 Tree species identification results

4.2 林分空间结构单元的构建

基于单木提取和树种识别结果,得到每个样地基本情况以及单木位置数据,将这些数据导入到ArcMap 10.2软件中,生成样地内的林木空间位置散点图,然后创建Voronoi图,将与中心木直接相邻的Voronoi多边形内的单木界定为相邻木,以此确定林分空间结构单元。由于位于样地边缘的参照木最近相邻木有可能在样地之外,因此会产生边缘效应。为了减小这种误差,采用距离缓冲法进行边缘矫正[18],将缓冲区内的样木仅仅作为相邻木,而不作中心木。鉴于样地较小,且参考样地株行距很小,将样地缓冲距离设置为1 m。最终所构建的林分空间结构单元见图6。

图6 空间结构单元Fig.6 Spatial structural element

4.3 毛竹扩张杉木林分空间结构参数变化及分析

毛竹扩张不同程度林分混交度、角尺度、大小比数的变化情况见图7。从图7可以看出,林分不同扩张程度的平均混交度分别为0.25、0.46、0.48、0.45和0.40,整体平均混交度在0.20～0.50,属于弱度混交,但接近中度混交状态。在扩张第1和第5阶段,林分平均混交度分别为0.25和0.40,处于弱混交状态,树种混交程度不高;在扩张第2和第4阶段,平均混交度分别为0.46和0.45,仍属于弱混交情况;在扩张第3阶段,平均混交度最高,为0.48,接近中度混交。林分不同扩张程度的平均角尺度分别为0.53、0.54、0.54、0.53和0.51,整体平均角尺度在0.50～0.55,处于随机分布的状态。林分不同扩张程度的平均大小比数分别为0.49、0.51、0.51、0.49和0.49,整体平均大小比数在0.48～0.51,处于亚优势和中庸状态。在扩张第1、第4和第5阶段,平均大小比数都为0.49,处于亚优势状态;在扩张第2和第3阶段,林分平均大小比数都为0.51,属于中庸状态。综上可知,毛竹扩张改变了原始林分林木生长态势。

图7 毛竹向杉木林扩张不同阶段空间结构指标变化Fig.7 Changes of spatial structure index of P.edulis at different stages of expansion to C.lanceolata forest

毛竹向杉木林扩张不同程度林分混交度、角尺度、大小比数整体变化趋势一致,都是先增加后减小,出现这种变化趋势的原因可能是,在扩张第1阶段毛竹占比极小,杉木占比极大,所以林分中大部分目标树是杉木,而杉木的相邻木大部分也还是杉木,导致林分平均混交度、大小比数较低;随着毛竹的扩张,毛竹数量逐渐增加,到了扩张第2阶段,混交度和大小比数有所增长;而到了扩张第3阶段,毛竹和杉木数量占比接近于1∶1,这时混交度和大小比数达到最大,而随着毛竹的继续扩张,毛竹在林分中的占比逐渐增大,导致混交度和大小比数逐渐减小。林分平均角尺度虽然随着扩张程度的加深呈现先增后减的趋势,但整体变化非常小,均处于随机分布范围,这可能是由于研究区位于自然保护区,人为干扰较少,林木生长水平分布格局总体趋近于随机分布的状态。

4.4 毛竹扩张杉木林分空间结构评价指数变化及分析

根据变异系数法的原理,利用式(4)确定角尺度、混角度、大小比数3个空间结构的权重,3个指数的标准差及权重见表3。其中角尺度的权重最大,混交度次之,大小比数最小,由此可知,角尺度对评价指数的影响最大。根据式(5)得到各个样地的林分空间结构评价指数,利用式(6)对其进行归一化,使其取值范围为[0,1]。利用Arc Map10.2插值工具,将各个样地林木空间结构评价指数进行插值,得到各个样地林分空间结构评价指数分布图(图8)。空间结构评价指数越高,林分空间结构越接近于理想状态。可以看出,4块研究区的空间结构评价指数在0.5～0.8,林分空间结构整体较好。

表3 角尺度、大小比数及混交度的标准差及权重Table 3 Standard deviation and weight of angle scale,size ratio and mixing degree

图8 林分空间结构评价指数分布Fig.8 Index distribution of stand spatial structure evaluation

依据毛竹向杉木林扩张不同阶段的划分,计算不同扩张程度林分空间结构评价指数(FSSCI)的均值,得到空间结构评价指数随扩张程度的变化情况(图9)。由图9可知,随着毛竹占比的增加,即扩张程度的加深,林分空间结构评价指数依次为0.25、0.50、0.51、0.47和0.40,整体呈现先增后减的趋势,这与所选取的3个林分空间结构指标(角尺度、混交度、大小比数)的变化规律一致。在扩张第1阶段,杉木几乎处于纯林阶段,此时林分物种单一,物种丰富度不足,林分空间结构整体较低,导致生态稳定性极其不理想,极易受到邻近毛竹的扩张。随着毛竹的扩张,空间结构评价指数迅速上升,到了毛竹扩张的第3阶段,林分空间结构评价指数达到最高,代表此时林分的生态稳定性处于最优。而到了扩张第5阶段,林分空间结构评价指数有所下降,这表明混交林的空间结构要优于毛竹或杉木纯林。

图9 空间结构评价指数随扩张程度变化Fig.9 The evaluation index of spatial structure changes with the degree of expansion

5 结论与讨论

5.1 结论

本研究单木提取总体精度在90%以上,毛竹和杉木的识别精度达到87%,满足试验精度需求。

随着毛竹扩张程度的加深,林分平均角尺度、混交度、大小比数和林分空间结构评价指数均呈现先增后减的趋势,林分整体平均混交度在0.20～0.50,属于弱度混交,但接近中度混交状态,平均角尺度在0.50～0.55,处于随机分布的状态,平均大小比数在0.48～0.51,处于亚优势和中庸状态,空间结构评价指数在0.5～0.8,在毛竹和杉木数量占比接近1∶1时,林分空间结构评价指数达到最大值,此时林分空间结构达到最优状态,林分空间结构评价指数能有效地反映毛竹扩张过程中林分空间结构的变化趋势。

5.2 讨论

常用的林分空间结构获取方法是通过将实地调查的单木位置信息,导入ArcGIS软件或Winkelmass软件中,计算得到各个空间结构参数,但这种方法仍需要消耗大量人力、物力进行单木位置数据的获取。为解决这一问题,研究利用无人机遥感,实现了单木位置的快速获取,为林业工作者高效获取森林资源信息提供便利。

研究所得毛竹向杉木林扩张过程中林分空间结构的演化规律与王秀云等[5]得到的毛竹扩张阔叶林变化规律整体一致,但在数值上又有些差异,王秀云等[5]得到的平均混交度在毛竹扩张前期达到0.663,扩张中期为0.54,扩张后期降为0,本研究的结果是在扩张前期为0.25、扩张中期达到0.48、扩张后期又降为0.40,这可能是因为所选取的研究区类型存在差异,王秀云等[5]选取的是常绿阔叶林、竹阔混交林和毛竹纯林,在扩张前期,虽然毛竹干扰不大,但常绿阔叶林分中树种较多,混交度较高,而到了扩张后期,演变为毛竹纯林,树种单一,混交度降为0,而本研究是按毛竹占比0%～20%、20%～40%、40%～60%、60%～80%、80%～100%来选取的,整个扩张过程中树种最多的是毛竹和杉木,因此混交度整体变化差异不大。

研究所构建的林分空间结构评价指数可为毛竹扩张机制研究提供基础数据,进而对毛竹扩张防控以及森林精准经营具有重要的指导意义。