基于增强Frost局部滤波及单木距离图重构标记的CHM树冠分割

张华聪,谭新建,喻龙华,厉月桥,陈永富,刘 仁,张怀清*

(1.中国林业科学研究院资源信息研究所,北京 100091;2.国家林业和草原局林业遥感与信息技术重点实验室,北京 100091;3. 中国林业科学研究院亚热带林业实验中心,江西 新余 336600)

激光雷达(light detection and ranging,LiDAR)是快速、准确获取林分空间结构的重要信息源[1-2],可对单木乃至林分层次上的森林参数进行精准估测及数字化表达,为森林资源监测[3-5]和三维建模提供数据基础[6-7]。树冠大小及形态是森林参数结构中的基本特征[8],树冠分割精度直接影响LiDAR提取单木树高、胸径以及坐标等其他信息的准确性。树冠分割数据源多为高空间分辨率遥感图像[9-12],但在郁闭度较高及林分背景复杂的林分中,容易受分辨率和树冠同质性影响。LiDAR因其空间属性被广泛应用于建筑、国土、交通、游戏等行业[13-16],也为快速、高效的树冠信息提取提供了更加精细的数据支撑[17-18],是目前国内外林业领域研究热点之一。

LiDAR提取单木树冠信息方法可分为2类,一是基于空间属性的点云分割,二是基于LiDAR生成的树冠高度模型(canopy height model,CHM)分割树冠。其中,基于空间属性的点云分割的原理是利用激光雷达点云间的空间结构关系和属性进行分割或分类。Ayrey等[19]图层堆叠算法,将整个森林点云进行分层切片,合并每个分割层内的同一单木轮廓即为分割树。Hui等[20]通过引入自适应搜索范围的方式改进均值漂移(mean shift)单木分割算法,在不需要复杂参数设置的情况下实现了很好的分割效果。随着人工智能的发展,PointNet[21]、快速回归卷积神经网络(faster R-CNN)[22]等深度学习方法也被应用与树冠分割,取得了较好的分割效果,相较于基于CHM树冠分割算法,避免了点云栅格化过程中造成的信息丢失[23]。但由于点云空间属性的复杂性,这些方法需要大量的先验知识,且深度学习算法的分割精度依赖于大样本和时间,对硬件要求高,多存在运算效率低的问题。CHM树冠分割的核心思想是利用树冠CHM中树顶亮度值大、边缘亮度低的特点进行分割,具有实操性强、数据量小等特点,方法主要有分水岭算法、局部最大值法、区域生长算法等。Vincent等[24]在1991年针对图像分割场景提出了经典的分水岭算法。Chen等[25]提出了一种基于标记控制的分水岭分割方法,通过历遍标记CHM中的树冠顶点,一定程度上避免了单木过分割。王鑫运等[26]利用局部最大值法分别在3种不同分辨率的冠层高度模型中搜索树顶位置,结合分层处理思想通过模拟泛洪算法实现单木分割。黄昕晰等[23]利用局部极大值算法和连通性生长算法优化树冠顶点提取,最后利用标记的分水岭分割方法实现树冠分割。虽然基于高分辨率CHM的方法具有较好的树冠分割效果,但对于高郁闭度林分而言,树冠间粘连重叠现象也使得无法简单的利用树冠边缘灰度值变化进行分割,复杂的树种组成也将导致树冠分割效果存在差异,且高分辨率CHM往往会出现更多的空洞和灰度变化不平滑现象。学者目前多采用Frost、Lee、高斯和中值等图像滤波方法对CHM进行填充与优化,如李岩等[27]在树冠提取前,利用高斯滤波结合半自动的空洞填充算法对CHM进行优化;张海清等[28]利用法中值滤波对CHM进行优化,结合DSM实现了单木树高提取。但上述滤波方法均会影响CHM非噪声点的真值。Lopes等[29]提出了增强Frost滤波方法,比原Frost滤波有着更好的细节保持能力,但依然会引起一些假细节边缘现象[30]。因此,改进CHM滤波优化和树冠CHM分割方法是提高LiDAR分割树冠的有效途径。本研究选取阔叶混交林、针阔混交林、针叶混交林3种不同林分类型为研究对象,利用采集的研究区LiDAR数据生成CHM,通过增强Frost局部滤波优化CHM,然后提出单木距离图标记重构技术对单木树冠进行分割与估算,比较分析不同分辨率和不同方法之间的提取效果,为森林结构参数的调查监测提供参考依据。

1 材料与方法

1.1 研究区概况

研究区位于江西省分宜县中国林业科学研究院亚热带林业实验中心山下实验林场(114°35′～114°40′E,27°40′～27°50′N)。该地区森林资源丰富、森林覆盖度高,地带性植被为常绿阔叶林,人工林树种主要以杉木(Cunninghamialanceolata)、马尾松(Pinusmassoniana)、木荷(Schimasuperba)、香樟(Cinnamomumcamphora)、刨花楠(Machiluspauhoi)等为主,属亚热带湿润性气候,雨量充沛,阳光充足,年降水量1 600 mm,年均气温15～17 ℃,全年无霜期270 d。

1.2 数据获取

1)激光点云数据。研究使用的无人机激光点云数据采集于2020年10月,坐标系为CGCS2000。无人机飞行平台为DJI M600 PRO,设定水平飞行速度为15 km/h,飞行高度为180 m,定位精度为10 cm;搭载的激光扫描仪是RIEGL miniVUX-1UAV激光扫描装置,配备了RIEGL miniVUX 3D激光扫描仪和Trimble APX-15 GNSS惯性集成导航OEM,获取的点云平均密度为 300点/m2。数据采集时天气晴朗,微风无云。在Global mapper 23.0中进行LiDAR点云拼接、裁剪、去噪和分类处理,并利用分类后的树冠点云生成不同分辨率的CHM数据。

2)实测数据。野外调查时间为2020年10月,在无人机所飞区域内随机选取针叶混交林分、阔叶混交林分以及针阔混交林分成熟人工林各1处,其中:针叶混交林林分密度为689株/hm2,树种主要包括杉木、水杉(Metasequoiaglyptostroboides)、马尾松等;针阔混交林林分密度为667株/hm2,树种主要包括杉木、香樟、刨花楠、木荷等;阔叶混交林林分密度为500株/hm2,树种主要包括香樟、枫香(Liquidambarformosana)、木荷。每林分设置30 m×30 m样地3个,通过每木检尺测量林木东西南北冠幅、树高、枝下高,并利用天宝TrimbleS3 RTK卫星定位仪采集每个单木的坐标信息,样地数据作为实测数据验证分割结果。

1.3 基于增强Frost的局部滤波

为改善高分辨率CHM中的孔隙与噪音问题,同时保持树冠间边缘及CHM的非噪声点真值,本研究提出了一种基于增强Frost的局部滤波方法。即先对树冠中间的噪声点进行判断,在对判定点集进行增强Frost滤波。

首先,设定滤波窗口大小均为3×3,以分辨率为0.2 m×0.2 m的CHM为实验对象,利用3×3卷积计算CHM中每个像元值与其周边共9个像元值的标准差S1,以及去除该像元值后周边8个像元值的标准差S2的差值V;通过设定阈值C=10,确定该像元是否属于噪声点,对于V大于阈值C的像元点进行标记为集合U1,V小于阈值C的像元点进行标记为集合U2。

(1)

V=S1-S2;

(2)

(3)

对判断的像元集合U1进行增强Frost滤波,增强Frost滤波是通过计算目标像素值的加权均值和目标像素周围圆形核的方差来确定每个目标像素的灰度值,其中加权值为目标像元到邻域像元绝对距离的倒数。在增强Frost滤波过程中,每个目标像元被分为3类。

(4)

式中:I为滤波卷积核内像元均值;IR为滤波卷积核内的加权均值;Ic为目标像元值原值;L为CHM的图像视数。

1.4 基于树冠标记及单木距离图重构分割算法

1.4.1 “树冠核心区”的提取与定位

利用形态学灰度重建方法[31]对CHM的局部极大值进行查找与二值化,定义查找到的局部极大值区域为树冠核心区,像元值设定为255,非局部极大值区域像元值设定为0,得二值化后“树冠核心”区图层。由于单木树冠顶点一般位于冠层的中心区域,对“树冠核心”在CHM进行定位,得到树冠分层梯度三值图F1,实现的具体路径为:

(5)

式中:Aij为树冠核心二值图像中第i行第j列处像元值;Bij为树冠图二值图像中第i行第j列处像元值;Hij为提取的树冠分层梯度三值图像第i行第j列处像元值(图1)。

1.4.2 编号“树冠核心”区

第1步,将树冠分层梯度三值图作为输入图像F1,该图的“核心区域”灰度值为255,“非核心区域”的灰度值为1,背景区域为0。

第2步,设无符号灰度变量为N,并初始化变量N=2,从上到下,从左至右依次扫描寻找树冠分层梯度三值图中灰度值为255的像素点并赋值为N,判断该像素点的8邻域内是否有灰度值等于255的像素点存在:若存在,则将灰度值为255的像素点的灰度值同样赋值为N,并重新判断被赋值点的8邻域内是否有灰度值等于255的像素点存在,若存在继续赋值,直至8邻域内无灰度值等于255的像元。

第3步,将N=N+1,重复第1步,至图F1中无灰度值等于255,将图中灰度值1值赋值为0,输出为图F2。

1.4.3 单木距离图重构与标记分割

第1步,以图像F1、图像F2为输入图像,F1图的“核心区域”灰度值为255,“非核心区域”的灰度值为1,背景区域为0;图像F2的“核心区域”灰度值为大于2的编号,其他区域为0。

第2步,设一无符号灰度变量为M,并初始化变量M=255,从上到下、从左至右依次扫描寻找图像F1中灰度值为1的像素点F1(i,j),判断的F1(i,j)8邻域内是否有灰度值等于M的像素点存在F1(i+x,j+y),若存在,则将F1(i,j)的灰度值重新赋值为M-1,同时将图像F2中对应点F2(i+x,j+y)中灰度值赋值给F2(i,j);如果不存在,则跳过当前扫描的像素点,继续扫描寻找下一个灰度值为1的像素点,直到全部扫描寻找完图F1中所有灰度值为1的像素点,该过程定义为单木距离图重构与标记。其中:i与j代表图像的行列号;i+x和j+y为8邻域位置信息,x、y取值1或-1。F(i,j)表示图像F中第i行,第j列位置信息。

第3步,将M-1赋值给M,重复第2步,直到图F1中的所有灰度值为1的像素点都被重新赋值,则停止扫描与赋值。最终得到基于距离图重建的单木间距离图像F3和标记分割图像F4。

1.5 精度评价

1.5.1 滤波评价指标

本研究取图像平滑指数(SI,式中记为IS)、信噪比(SNR,式中记为RSN)、等效视数(ENL,式中记为NEL)以及边缘保持指数(EPI,式中记为IEP)4个指标来评估CHM滤波优化效果。SI值越高,则表示平滑作用越强,斑点去除效果越好,滤波后灰度越稳定。SNR是衡量滤波器去噪能力的通用指标,信噪比值越高,则影像的去噪效果越好。ENL用于客观评价匀质区域的平滑程度,也就是对于噪声的抑制程度。等效视数越大,表示斑点抑制能力越强,噪声越弱。EPI的取值范围是[0,1],边缘保持的定义为滤波后影像的边缘与滤波前之比。边缘保持系数越接近1,则边缘保持的效果越好。在每种林分类型中随机选取3个30 m×30 m的区域作为感兴趣区求算评价指标,指标计算公式如下:

(6)

RSN=10log[max(I)2/σ2];

(7)

(8)

(9)

1.5.2 分割评价指标

利用总体精度(OA,式中记为OA)、误判误差 (CE,式中记为CE)、漏判误差(OE,式中记为OE)3个指标评价树冠分割结果的好坏。参考样地内实测位置信息,并利用正射影像目视解译出真实的单木冠幅范围,通过过对比真实树冠投影区域与分割后树冠投影区域来判别分割正确性[32-33]。指标的计算公式分别如下:

(10)

CE=FN/TN×100%;

(11)

OE=FP/TN×100%。

(12)

式中:TN为样地实际树冠个数;TP为分割正确的真实树冠个数,判定标准为真实树冠投影区域与分割后投影区域重叠率 70% 以上;FN为错分的真实树冠个数,判定标准为分割后投影区域与真实树冠投影区域低于重叠率 70%,或判定标准为分割树冠投影区域与多个真实树冠投影区域重叠,且每个真实树冠投影区域均未占到70%以上;FP为漏判的真实树冠个数,判定标准为真实树冠投影没有分割树冠投影与之匹配重叠。

1.5.3 冠幅估算精度评价指标

冠幅估算精度的评价指标为决定系数(coeffcient of determination,R2)及实际观察值与估计值之间的差(REPIdual,δ),计算公式如下:

(13)

(14)

2 结果与分析

2.1 3种林分类型的滤波结果

采用不同滤波方法对阔叶林冠层、针叶林冠层及针阔混交林冠层3种类型的滤波结果规律基本一致(图2、表1)。总体而言,增强Frost、增强Lee 2种滤波方法所得到图像的信噪比(SNR)及等效视数(ENL)均高于局部Sigma和本研究的滤波方法,去噪能力强,对CHM有着很好的平滑作用,但边缘保持指数(EPI)显著低于局部Sigma和本研究的滤波方法,图像细节保留能力差。

表1 不同滤波方法指标比较

A、B、C分别表示阔叶混交林、针叶混交林及针阔混交林冠层。下同。A, B and C respectively represent the canopy layer of mixed broad-leave forest, mixed coniferous forest, and coniferous broad-leaved forest. The same below.冠层类型后1、2、3、4、5、6分别表示原始图像、增强Frost滤波结果、增强Lee滤波结果、局部Sigma滤波结果、中值滤波结果和本研究方法滤波结果。In the figure,1, 2, 3, 4, 5 and 6 respectively represent the original image, the result of the enhanced Frost filtering, the result of the enhanced Lee filtering, the result of the local Sigma filtering, the result of median filtering, and the result of the filtering method used in this study.图2 不同滤波方法的结果比较Fig. 2 Comparison of results of different filtering methods

局部Sigma滤波边缘保持指数(EPI)最高,为0.89,边缘保持能力最好,图像细节保留程度优势明显,但是去噪能力较弱,平滑指数(SI)、信噪比(SNR)及等效视数(ENL)均为最低。本研究使用的增强Frost局部滤波算法的平滑指数(SI)和等效视数(ENL)均值分别为2.85和8.28,优于局部Sigma和中值滤波,斑点去除效果好,滤波后灰度稳定;信噪比(SNR)为14.02;边缘保持指数(EPI)为0.86,仅次于局部Sigma滤波,使得图像细节保留能力均明显优于增强Frost、增强Lee和中值滤波。

2.2 不同分辨率CHM分割结果

单木距离图重构标记分割方法对滤波优化后不同分辨率CHM的分割结果见图3。

从图3可知,可以看到不同分辨率下针叶混交林的分割效果均为最优,其次是针阔混交林,然后是阔叶混交林。在3种林分中,0.1 m分辨率的CHM均会出现过分割现象,0.2 m分辨率的冠层CHM分割效果最好,0.5、1.0 m分辨率的冠层CHM出现了欠分割的情况,且随着分辨率数值的增加,欠分割现象越严重。

2.3 不同分割方法的结果比较

以均值漂移算法类(图4A-3、4B-3、4C-3)、分水岭方法(图4A-4、4B-4、4C-4)以及本研究方法(图4A-5、4B-5、4C-5)的进行分割。

冠层类型后1、2、3、4、5分别表示原始点云数据及样地位置(按高程着色)、滤波处理后的CHM图像、均值漂移方法分割结果、分水岭分割方法和本研究方法分割结果,为样地边界,为单木位置。In the figure, 1, 2, 3, 4 and 5 represent the original point cloud data and plot location (colored by elevation), the CHM image after filtering, the segmentation results of the mean shift method, the watershed segmentation method, and the segmentation results of the method used in this study. denotes the plot boundaries, and indicates the position of individual trees.图4 不同分割方法的树冠分割结果Fig. 4 Tree crown segmentation results of different segmentation methods

对0.2 m分辨率CHM的分割结果进行统计可知(表2),阔叶混交林、针叶混交林、针阔混交林3种不同林分类型的样地平均分割精度OA分别为0.75、0.96和0.84,误判误差 (CE)分别为0.23、0.03和0.12,漏判误差(OE)分别为0.02、0.01 和0.04。分割总体精度均以针叶混交林为最高,针阔混交林次之,分割效果最差的为阔叶混交林,3种不同林分类型的样地分割过程中漏判误差均较低,所有单木树冠均进行了树冠分割。同时,证明本研究提出方法的优势,将本研究方法分割的0.2 m分辨率CHM分割结果与分水岭分割方法(Watershed)、F算法(Mean shift)的分割结果进行了对比,本研究方法的总体精度(OA)、误判误差 (CE)、漏判误差(OE)指标总体均值分别为0.85 、0.13 和0.02,均优于分水岭算法的0.64、0.34、0.02以及均值漂移算法的0.55、0.43和0.02。

表2 不同方法分割精度比较

2.4 树冠估算结果

在本研究方法分割结果的基础上,测量东西、南北长度均值为林分冠幅估测值,选取样地中估算的298棵单木冠幅数据,以样地区域实测值为验证数据进行回归及残差分析,结果如图5所示。本研究方法估测针叶混交林、针阔混交林和阔叶混交林树冠值的R2分别为0.83、0.82及0.71,其中针叶树种的树冠估算精度要优于阔叶混交林,阔叶混交林的估算残差较针叶混交林、针阔混交林更均匀,而针阔混交林出现高估。

图5 基于本研究方法的树冠提取精度Fig. 5 Crown extraction accuracy based on the method in this paper

3 讨 论

对同一株单木而言,树冠孔隙是造成CHM空洞和无效值的主要原因,这种噪音往往随着CHM分辨率的升高而增加,同时CHM质量决定着树冠分割效果。目前,对图像进行滤波是CHM优化的主要方式。李岩等[27]和张海清等[28]在利用CHM进行树冠分割前均进行了滤波处理,分割效果得到了改善。但传统的滤波在抑制噪声与保留细节信息上存在矛盾,树冠边缘需要更多的细节,树冠中间则需要更少的灰度值突变。本研究提出的增强Frost局部滤波方法利用方差变异情况对灰度值突变进行检索,在既减少CHM孔隙和噪声的同时又能保持树冠间边缘等区域CHN非噪声点的真值,不同树种类型基本不影响滤波效果,该滤波方法具有普适性。同时不同分辨率的CHM在分割效果上也存在这显著差异,虽然低分辨率的CHM的空洞和无效值显著降低,但随着分辨率的降低树冠间像元同质性显著提升,导致欠分割现象逐级增,表明高分辨的CHM在树冠分割效果上优于低分辨率CHM。

林分类型也是影响树冠分割的主要因素,在本研究实验的针叶混交林、阔叶混交林以及针阔混交林3种森林类型中,针叶林树冠由于边界清晰且树冠顶点凸出单一,具有最好的分割效果,而在针阔混交林中的阔叶树,由于针叶树种的存在,较在阔叶混交林中更易得到识别,分割效果较好。自然状态下的阔叶树往往存在多个平行顶点和树簇,在林分垂直结构上具有更高的复杂性,无论是基于CHM的分割方法还是基于点云的分割聚类方法,分割效果的都有待提高。孙振峰等[12]利用标记分水岭算法提取树冠,其研究结果表明杉木林单木树冠提取效果结果要优于桉树林。Chen等[21]利用深度学习PointNet算法进行单木分割,结果表明以针叶树种为主的圃地内树冠分割效果最好,其次是混交林,分割效果最差的是落叶阔叶林,这与本研究的结论一致。此外,在不同林龄、树种以及卷积核大小同样可能影响本研究分割方法的效果,可在接下来的研究进行深入探讨。

林分中单木树冠的大小差异及粘连现象使得均值漂移算法和分水岭算法的分割结果并不理想,距离图重构标记分割方法是在提取树冠局部最大值和林冠范围的基础上,基于距离变换的思想,在生成距离矩阵时进行分割标记,有效控制树冠分割范围的同时对形状进行了一定约束,避免了单木距离变换后的二次分割。分割效果虽然较分水岭算法和均值漂移算法有明显的改善,在阔叶林中依然存在一定的过分割,在今后的研究中可对不同林分类型及冠幅大小设定自适应窗口进行探讨与改进。

本研究方法的树冠分割结果可满足单木冠幅的估算需求,分割结果的优势使得针叶混交林树种的树冠估算精度要优于其他2种林分类型;但阔叶混交林的估算残差较针叶混交林、针阔混交林更均匀,针阔混交林更易出现高估。这可能是因为在郁闭度较低的针叶混交林中,地形因子和林下植被均会提高点云分布的复杂性[28,34],影响分类精度,增加CHM的分割范围,导致冠幅相对高估。

本研究以阔叶混交林、针阔混交林、针叶混交林3种不同林分类型为研究对象,通过改进增强Frost局部滤波和单木距离图标记重构技术对CHM进行单木分割与树冠估算,结果表明:

1)研究所使用的增强Frost局部滤波在优化CHM的效果上优于传统的增强Frost、增强Lee、局部Sigma以及中值滤波,在既减少CHM孔隙和噪声的同时又能保持树冠间边缘等区域CHN非噪声点的真值。

2)基于距离图重构标记的树冠分割方法优于传统的分水岭和均值漂移算法,在多种林分类型单木树冠分割上具有较高的精度,方法具有一定的鲁棒性,树冠提取精度满足森林调查和监测的基本需求。