CN115880487A - 基于深度学习方法的林木激光点云枝叶分离方法 - Google Patents

基于深度学习方法的林木激光点云枝叶分离方法 Download PDF

Info

Publication number
CN115880487A
CN115880487A CN202111137023.3A CN202111137023A CN115880487A CN 115880487 A CN115880487 A CN 115880487A CN 202111137023 A CN202111137023 A CN 202111137023A CN 115880487 A CN115880487 A CN 115880487A
Authority
CN
China
Prior art keywords
points
point cloud
point
feature
branch
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111137023.3A
Other languages
English (en)
Inventor
薛联凤
云挺
丁竹娴
曹林
白颢
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Forestry University
Original Assignee
Nanjing Forestry University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing Forestry University filed Critical Nanjing Forestry University
Priority to CN202111137023.3A priority Critical patent/CN115880487A/zh
Publication of CN115880487A publication Critical patent/CN115880487A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Image Processing (AREA)

Abstract

本发明公开了基于深度学习方法的林木激光点云枝叶分离方法,包括获取林木激光点云数据;将去燥后的激光点云数据分为地上点和地面点;采用机器学习算法以及人工修正方式对地上点进行枝叶分离操作,对已进行枝叶分离操作的地上点进行体素化剖分;将已经进行枝叶分离操作和体素化剖分的点云数据作为训练样本数据集,并采用数据增广的方法对训练样本数据集进行扩增,采用扩增后的训练样本数据集对深度学习网络开展训练,得到训练好的深度学习网络模型;通过训练好的深度学习网络模型中,实现体素内的点云数据的枝叶分离。本发明通过构造的深度学习网络模型可更有效地提取林木点云的全局与局部特征信息,实现点云枝叶分离,准确率高,稳定性好。

Description

基于深度学习方法的林木激光点云枝叶分离方法
技术领域
本发明属于林业技术领域,具体涉及一种基于深度学习方法的林木激光点云枝叶分离方法。
背景技术
树木骨架空间模型的精准重建在研究林木信息化资源调查和树木表型结构特征反演中发挥着至关重要的作用。近些年,针对树木模型重建的方法大致可以分为两类,分别为基于图像与激光点云的树木重建方法。这两种方法在数据获取上存在差异性,但提取树木空间特征的方式上又具有一定的相似性。
基于图像的树木枝干重建方法通常采用双目、多目视觉及深度相机的方法捕获树木的纹理与景深信息,结合曲率约束实现树木枝干重建,或利用多目视觉的多张图像通过交互式编辑重建树木枝干,并通过叶片分割重建叶子模型。基于图像的树干重建方法虽然简单易行,但提取的树木枝干处于二维层面,很难直观反映树木的三维形态,且提取结果会被多样的环境背景、复杂的树干空间拓扑结构、立体匹配误差、视角遮挡等多方面因素干扰。
随着激光测量技术的进步,树木的数字化建模已逐渐从基于二维图像的方法转变为从三维激光点云中生成。激光扫描具有精度高、密度高等特性,同时,可结合计算机图形学、机器视觉等理论算法开展林木三维点云的处理分析,如:利用机器学习的方法对林木点云进行枝叶分类,利用图论将采集到的点云数据存储为八叉树结构,采用结合点法线的PROSAC算法建模提取林分树干;或使用空间填补的三角面片构造每片树叶来建立树木冠层模型等。
随着新一代人工智能技术不断取得应用性突破,以及激光扫描技术的落地应用,使深度学习技术在三维点云处理中崭露头角,一些知名的深度网络模型应运而生。如2016年提出的PoineNet网络设计了基于对称函数的多层感知机与基于空间转换不变性思想的T-Net架构,解决了输入点云的无序性和部分的旋转不变性。2017年出现的三个面向激光点云处理的深度学习网络如:在体素尺度下对三维点云进行划分并开展分类检测与位置回归的VoxelNet,借鉴卷积神经网络(CNN)的多层感受野的思想并结合球邻域提取点云局部空间特征的PointNet++,以及设计了χ变换与连续的权重函数卷积操作来实现空间目标点云卷积、局部特征抽取与分类的PointCNN。除此之外,基于图论思想的点云采样、分组和池化的深度网络模块DPAM和结合分层数据结构(K-d树)的网络3DContextNet也被相继提出。
虽然深度学***面上的特征开展识别或直接分析处理空间林木点云数据,如:构造卷积神经网络和深度置信网络处理从激光点云生成的肯塔基大学罗宾逊森林的数字表面模型(DSM)和厦门市树木的侧面投影图来开展针叶与不同阔叶树种的分类;或采用Superpoint graph、PointNet、U-Net、Kd-Net和FasterR-CNN等深度网络模型在体素尺度下开展欧洲山毛榉的枝叶分离、安徽池州亚热带森林的单株树冠分割、加拿大芬兰的部分树种识别和海南橡胶林参数提取等任务。
尽管深度学习在激光点云分析上取得了一些进展,但依然存在如下问题:1)树冠中枝干骨架具有拓扑结构各异、形式多样且受到互相遮挡等因素,其对网络的鲁棒性与普适性都有很高的要求,如何设计一个更精准的枝叶点云分割深层网络还处在研究阶段;2)深度学习能够把待分析的二维数据通过卷积操作映射到高维空间中开展特征维度的增广与压缩,但对于三维空间点云的无序性和旋转变化性,很难设置在多尺度空间下开展卷积操作并提取点云的全局与局部有效特征。
发明内容
本发明所要解决的技术问题是针对上述现有技术的不足提供一种基于深度学习方法的林木激光点云枝叶分离方法,本基于深度学习方法的林木激光点云枝叶分离方法构造具有多个特征编码层与特征解码层的枝叶点云分类深度学习网络,可更有效地提取林木点云的全局与局部特征信息,通过该深度学习网络模型实现点云枝叶分离,准确率高,稳定性好。
为实现上述技术目的,本发明采取的技术方案为:
基于深度学习方法的林木激光点云枝叶分离方法,包括:
步骤1:获取林木激光点云数据;
步骤2:对采集的激光点云数据进行去燥处理,将去燥后的激光点云数据分为地上点和地面点;
步骤3:采用机器学习算法以及人工修正方式对地上点进行枝叶分离操作,对已进行枝叶分离操作的地上点进行体素化剖分;
步骤4:将已经进行枝叶分离操作和体素化剖分的点云数据作为训练样本数据集,并采用数据增广的方法对训练样本数据集进行扩增,得到新的训练样本数据集;
步骤5:构建深度学习网络,采用新的训练样本数据集对深度学习网络开展训练,得到训练好的深度学习网络模型;
步骤6:采集待测林地的激光点云数据,对采集的激光点云数据进行去燥处理,将去燥后的激光点云数据分为地上点和地面点;对地上点进行体素化剖分,将单个体素内的点云数据输入到训练好的深度学习网络模型中,以实现体素内的点云数据的枝叶分离。
作为本发明进一步改进的技术方案,所述深度学习网络包括特征编码层和特征解码层;所述特征编码层用于采用下采样方法提取点云特征信息,所述特征解码层用于采用上采样方法传递点云特征信息;
所述特征编码层包括采样模块、第一分组模块和第一PointConv模块;所述特征解码层包括特征插值模块、第二分组模块和第二PointConv模块。
作为本发明进一步改进的技术方案,所述特征编码层共四个,所述特征解码层共四个。
作为本发明进一步改进的技术方案,第j个特征编码层中采样模块的计算过程为:
输入大小为(d+Cj)×Nj的点云矩阵,通过最远点采样法从该点云矩阵中选取Nj+1个子采样点,得到大小为(d+Cj+1)×Nj+1的点云矩阵;
其中:(d+Cj)×Nj是指具有d维坐标和Cj维特征的Nj个点,(d+Cj+1)×Nj+1是指具有d维坐标和Cj+1高维特征的Nj+1个子采样点。
作为本发明进一步改进的技术方案,第j个特征编码层中第一分组模块的计算过程为:
输入采样模块输出的d×Nj+1的子采样点的坐标和最初输入的(d+Cj)×Nj的特征点集,采用最近邻规则分类法,在采样模块提取出的Nj+1个采样点中每个点的邻域内寻找最近的K个近邻点组成一个小组,得到d×K×Nj+1近邻点相对于每个采样中心点的坐标索引集合,并根据坐标索引得到(d+Cj)×K×Nj+1近邻点的特征集合,其中K是Nj+1个采样中心点中每个点的邻域内近邻点的个数。
作为本发明进一步改进的技术方案,第j个特征编码层中第一PointConv模块的输入是三部分数据;其中第一部分数据为:Nj+1个中心点附近K个近邻点的特征,大小为(d+Cj)×K×Nj+1;第二部分数据为:Nj+1个中心点附近K个近邻点相对于中心点的局部坐标,大小为d×K×Nj+1;第三部分数据为:Nj+1个中心点附近K个近邻点的密度,大小为1×K×Nj+1;第一PointConv模块输出大小为(d+Cj+1)×Nj+1的新的局部区域。
作为本发明进一步改进的技术方案,所述第j个特征编码层中第一PointConv模块的计算过程为:
步骤(1)、利用核密度估计方法计算中心点
Figure BDA00032824094700000415
附近K个近邻点中每个输入点的核密度估计,再通过核密度估计计算K个近邻点的逆密度变换,再逆密度变换的结果输入多层感知机MLP2进行激活函数的非线性变换,得到大小为1×K的逆密度系数,对逆密度系数进行Cin次的复制平铺,输出大小为Cin×K的逆密度张量;
步骤(2)、将K个近邻点相对于中心点的局部坐标输入多层感知机MLP1,输出大小为(Cin×Cout)×K的权重张量W,其中Cin为输入特征数,Cout为输出特征数;
步骤(3)、把步骤(1)输出的大小为Cin×K的逆密度张量和K个近邻点的特征Fin进行点乘运算,输出大小为Cin×K的矩阵
Figure BDA0003282409470000041
再对/>
Figure BDA0003282409470000042
进行Cout次的复制平铺,输出大小为(Cin×Cout)×K的中间张量;
步骤(4)、把步骤(2)输出的大小为(Cin×Cout)×K的权重张量W和步骤(3)输出的大小为(Cin×Cout)×K的中间张量进行点乘运算,输出大小为(Cin×Cout)×K的矩阵
Figure BDA0003282409470000043
步骤(5)、对步骤(4)得到的矩阵
Figure BDA0003282409470000044
中的第一维和最后一维进行求和运算,得到中心点/>
Figure BDA0003282409470000045
的输出特征Fout
其中中心点
Figure BDA0003282409470000046
为Nj+1个中心点中的某一点。
作为本发明进一步改进的技术方案,第j个特征解码层中特征插值模块的计算过程为:
以第j-1个特征编码层中采样模块得到的Nj个采样点中每个点
Figure BDA0003282409470000047
为中心,在第j个特征编码层中采样模块得到的Nj+1个采样点中,选取距离Nj个采样点中每个点/>
Figure BDA0003282409470000048
最近的三个邻域点/>
Figure BDA0003282409470000049
得到三个邻域点/>
Figure BDA00032824094700000410
在Nj+1个采样点中的索引,其中t=1,2,3;
按照距离倒数计算出三个邻域点
Figure BDA00032824094700000411
相对于Nj个采样点中对应中心点/>
Figure BDA00032824094700000412
的权重值/>
Figure BDA00032824094700000413
计算公式为:
Figure BDA00032824094700000414
其中
Figure BDA0003282409470000051
表示点/>
Figure BDA0003282409470000052
与点/>
Figure BDA0003282409470000053
之间的欧氏距离,μ表示加权幂指数;
通过这三个邻域点的索引在第j个特征解码层最初输入的(d+C′j)×Nj+1特征点集中得到这三个邻域点对应的特征(d+C′j)×3,再结合对应的权重值
Figure BDA0003282409470000054
对这三个邻域点的特征进行求和平均,得到/>
Figure BDA0003282409470000055
对应的(d+C′j)×1的插值特征/>
Figure BDA0003282409470000056
其中插值特征/>
Figure BDA0003282409470000057
的计算公式为:
Figure BDA0003282409470000058
其中c′t,j∈C′j (2);
因为第j-1个特征编码层共有Nj个采样点,所以得到了(d+C′j)×Nj的插值特征;
得到(d+C′j)×Nj的插值特征后,运用跳跃连接方法将(d+C′j)×Nj插值特征与第j-1个特征编码层的输出的(d+Cj)×Nj的特征进行连接操作,得到特征插值模块的最终输出结果:(2d+Cj+C′j)×Nj
其中(d+C′j)×Nj表示具有d维坐标和C′j维特征的Nj个点。
作为本发明进一步改进的技术方案,第j个特征解码层中第二分组模块的计算过程为:
输入特征插值模块的最终输出结果(2d+Cj+C′j)×Nj和第j-1个特征编码层中采样模块得到的d×Nj的子采样点的坐标,采用最近邻规则分类法,在Nj个采样点中每个点的邻域内寻找最近的K′个近邻点组成一个小组,得到d×K′×Nj近邻点相对于每个采样中心点的坐标索引集合,并根据坐标索引得到(2d+Cj+C′j)×K′×Nj近邻点的特征集合,其中K′是Nj个采样中心点中每个点的邻域内近邻点的个数。
作为本发明进一步改进的技术方案,第j个特征解码层中第二PointConv模块的输入是三部分数据;其中第一部分数据为:Nj个中心点附近K′个近邻点的特征,大小为(2d+Cj+C′j)×K′×Nj;第二部分数据为:Nj个中心点附近K′个近邻点相对于中心点的局部坐标,大小为d×K′×Nj;第三部分数据为:Nj个中心点附近K′个近邻点的密度,大小为1×K′×Nj;第二PointConv模块输出大小为(d+C′j-1)×Nj新的局部区域。
本发明的有益效果为:
本发明把带有枝叶标签的林木激光点云数据体素化剖分作为训练集;构造具有多个特征编码层与特征解码层的枝叶点云分类深度学习网络,该深度学习网络模型包括采用核密度估计计算局部点云特征的改进PointConv模块和通过多尺度特征融合的点云特征插值模块,可更有效地提取林木点云的全局与局部特征信息,通过该深度学习网络模型实现点云枝叶分离,准确率高,稳定性好,无需人工参数调节与结果干预修正。
附图说明
图1中(a)为研究区域海南省儋州市的位置示意图。
图1中(b)为海南岛儋州市橡胶树种植园内的三块橡胶林地示意图。
图2中(a)为橡胶树品种PR107的部分样地运用机器学习方法加人工修正得到的枝叶分割与体素剖分结果示意图。
图2中(b)为橡胶树品种CATAS 7-20-59的部分样地运用机器学习方法加人工修正得到的枝叶分割与体素剖分结果示意图。
图2中(c)为橡胶树品种CATAS 8-79的部分样地运用机器学习方法加人工修正得到的枝叶分割与体素剖分结果示意图。
图3为面向点云枝叶分离的深度学习网络架构图。
图4中(a)为中心点pi联合其近邻点开展局部特征计算示意图。
图4中(b)为改进的PointConv模块框架图。
图5中(a)为深度学习网络分类精度曲线示意图。
图5中(b)为深度学习网络损失函数曲线示意图。
图6中(a1)为橡胶树品种PR107测试样地的激光点云数据以及枝叶分离结果示意图。
图6中(a2)为橡胶树品种CATAS 7-20-59测试样地的激光点云数据以及枝叶分离结果示意图。
图6中(a3)为橡胶树品种CATAS 8-79测试样地的激光点云数据以及枝叶分离结果图。
图6中(b1)为橡胶树品种PR107测试样地枝干圆柱体拟合结果图。
图6中(b2)为橡胶树品种CATAS 7-20-59测试样地枝干圆柱体拟合结果图。
图6中(b3)为橡胶树品种CATAS 8-79测试样地枝干圆柱体拟合结果图。
图6中(c1)为橡胶树品种PR107测试样地主枝干和一级枝干分类结果图。
图6中(c2)为橡胶树品种CATAS 7-20-59测试样地主枝干和一级枝干分类结果图。
图6中(c3)为橡胶树品种CATAS 8-79测试样地主枝干和一级枝干分类结果图。
图6中(d1)为橡胶树品种PR107测试样地基于空间分水岭和一级分枝的树叶聚类结果与叶团簇的凸包显示图。
图6中(d2)为橡胶树品种CATAS 7-20-59测试样地基于空间分水岭和一级分枝的树叶聚类结果与叶团簇的凸包显示图。
图6中(d3)为橡胶树品种CATAS 8-79测试样地基于空间分水岭和一级分枝的树叶聚类结果与叶团簇的凸包显示图。
图7中(a)为三块橡胶树样地树木胸径实测值与计算值的比较结果图。
图7中(b)为三块橡胶树样地计算的树木主枝干与一级枝干夹角与实测值的比对结果图。
图7中(c)为三块橡胶树样地树木一级枝干宽度和对应的叶团簇体积的相关性图。
具体实施方式
下面根据附图对本发明的具体实施方式作出进一步说明:
本实施例提出了深度学习与计算机图形学相融合的树木骨架重建与林木参数反演方法。首先,把带有枝叶标签的橡胶林激光点云数据体素化剖分作为训练集;接着,构造具有多个特征编码层与解码层的枝叶点云分类深度学习网络,该网络模型包括采用核密度估计计算局部点云特征的改进PointConv模块和通过多尺度特征融合的点云插值模块,可更有效地提取林木点云的全局与局部特征信息;再次,面向分割后的林木枝干点云,运用计算机图形学方法建立骨架中心点的连通性,并用圆柱体开展分段枝干拟合,进而从枝干点云中重建整片林地的树木骨架三维模型;最后,与人工实测数据相比对,在叶团簇尺度下精细刻画样地内每株树木,并对不同品种橡胶树拓扑结构与表型参数进行分析。
1.1研究区与数据采集:
研究区位于海南省儋州市,中国热带农业科学院橡胶研究所试验田(北纬19°32′47.89″,东经109°28′29.33″),如图1所示。作为中国第二大橡胶生产基地,海南岛属于热带季风气候,每年5月到10月为雨季,11月到4月为旱季。年平均降水量1815毫米,年平均气温约23.2℃。这一气候有利于橡胶产业的发展,且该地区橡胶树的种植量也在不断增加。经过数十年来的引进、试验,RRIM660、PR107、PB86、海垦2号、CATAS 7-20-59、CATAS 8-79、热垦126等品种的橡胶树逐渐发展成为适宜海南地域规模种植的且具有抗不同自然扰动、复生力强、产量高等特点的橡胶树优良品种。其中,PR107生长速度中等、干胶含量高、耐刺激、抗风性好,是一个抗风高产的橡胶树优良品种;CATAS 7-20-59作为PR107和RRIM600的杂交后代,速生高产而相对晚熟,是抗风性能好、目前普遍推广的橡胶树新品种;CATAS 8-79再生能力强、早熟、产量高且产量平稳,是一个综合性状优良的高产品系。因此,如图1所示,选择了PR107品种橡胶树(1号样地)、CATAS 7-20-59品种橡胶树(2号样地)和CATAS 8-79品种橡胶树(3号样地)作为试验的典型树种。
激光雷达数据于2019年10月15日获得,使用可在背包模式下操作的VelodyneHDL-32E高清激光雷达传感器。该传感器有32对激光/探测器,用于测量具有以下设置参数的橡胶树图:+10.67°至-30.67°垂直视野(FOV),角度分辨率为1.33°,360°水平视野,角度分辨率为0.16°,帧速率为10HZ,测量范围为70m。Velodyne HDL-32E扫描***由一名实验者携带,扫描仪被设置为“连续拍摄模式”,以每秒10转的速度采集数据。实验者背负Velodyne激光扫描***,并根据预定义的测量路线在三个橡胶树地块内行走。调查路线编程为预定义的“来回平行”路线(如图1虚线),两条平行测绘路线之间的间隔为6m,旨在尽可能覆盖三个研究样地的所有林木。同时,由于橡胶树种植园的复杂地形和沉重的扫描仪器,实验者以0.5m/s的速度沿着调查路线移动。Velodyne激光雷达***结合同步定位和映射(SLAM)技术,以快速完成每次扫描的配准,并为每棵目标橡胶树生成高密度点云。三个橡胶树地块获得的激光雷达数据的平均分辨率约为2cm。
图1为研究区域概况,图1中(a)显示了研究区域海南省儋州市的位置,图1中(b)显示了海南岛儋州市橡胶树种植园内的三块橡胶林地;图1中(b)为谷歌地图显示的三块橡胶林地的遥感图像,用不同的矩形标记了不同橡胶林地块的边缘,虚线表示机载激光雷达探测路线。
1.2训练样本与测试样本:
通过激光扫描仪获取树木的点云数据后,采用高斯滤波方法对采集的点云进行了去噪处理,并采用点云地面点滤波(CSF)的方法将去噪后的点云分为地上点和地面点。接着用传统的机器学习算法加人工目视检验与结果修正对试验区地上点进行了枝叶分离操作。从三块橡胶林样地中创建了三个子集,作为后续深度网络的训练样本。每个子集由代表相应橡胶林样地的约200×100m的区域组成,图2中(a)、(b)、(c)的左侧分别显示了三块橡胶林地部分子集数据的枝叶分离结果。同时由于三块样地上不同品种橡胶树(PR107、CATAS7-20-59和CATAS 8-79)具有相同的东西(3m)与南北(7m)种植间距,且林下灌木较少,因此对三块训练样地设置了相同的体素大小(长3m,宽7m,高7m)来对地上林木点云进行体素化剖分,保证每棵树位于每个体素的中心。剖分后,三个训练集的所有激光点云被划分给相应体素,并保证单棵树点云能被上下2个连续体素剖分。接着,采用抽稀化处理,每个体素内的点云随机采样为4096个点,图2中(a)、(b)、(c)的右侧显示了提取的以上下体素为边界的单棵树的枝干与叶子点云。经过以上操作后,共收集2743棵(每个品种900棵左右)枝叶分离后的橡胶树点云数据作为深度学习的训练集。
测试样本由三个橡胶树品种的三个样地子集组成,并保证与训练样本不相交。每个测试样地为35×35m的区域。通过背包激光雷达在相同的扫描分辨率下获取样地数据。整个研究区橡胶树具体参数如表1所示,其中冠幅与冠积根据单株分离算法计算获取。
表1为研究样地中橡胶树参数:
Figure BDA0003282409470000091
大量充足的训练样本是进行深度学习网络参数优化的基础,为了避免过拟合并提升网络的泛化能力,有必要收集尽可能多的训练样本集来参与深度学习网络的训练过程以优化神经元权值。本研究中,使用数据增广技术获取更多的训练样本,在不收集和处理新林木点云数据的情况下,基于现有体素内的树木点云,提出了一种沿随机向量以较小偏移量移动体素中的每个点
Figure BDA0003282409470000092
的策略,即高斯噪声抖动,具体过程如式(1)和(2)所示;
Figure BDA0003282409470000093
其中,a和b是两个0到1之间的随机数,λ为一个功率放大系数,r为产生的扰动,接着把每个点的扰动
Figure BDA0003282409470000095
带入式(2)中,将原来的树木点云从pi变为pi′,式(2)中μnoise和σnoise为高斯噪声的均值与方差。
Figure BDA0003282409470000094
或通过让体素内点云按照空间旋转方程沿Z轴(即树冠中轴)旋转较小的角度(如1~2°),通过调节方程参数,保证变换后出现一个相较原始树木点云的较小的标准偏差值(范围为1~5cm)。综上,橡胶树每个品种的样本增加1100个左右,总共训练样本被扩增至6000个。
图2为运用机器学习加人工修正的方法得到的三块橡胶林训练样地子集点云数据枝叶分割与体素剖分结果;右侧显示了以体素为边界的单棵树木枝干和叶子点云放大图。其中分类的叶子点云显示为灰色、分类的枝干点云显示为黑色;这些数据将作为训练样本集代入深度学习网络开展网络参数获取。图2中(a)为橡胶树品种PR107的部分样地运用机器学习方法加人工修正得到的枝叶分割与体素剖分结果示意图。图2中(b)为橡胶树品种CATAS7-20-59的部分样地运用机器学习方法加人工修正得到的枝叶分割与体素剖分结果示意图。图2中(c)为橡胶树品种CATAS 8-79的部分样地运用机器学习方法加人工修正得到的枝叶分割与体素剖分结果示意图。
1.3基于深度学习网络的林木枝叶分离:
以单个体素内带枝和叶标签的点云作为输入,训练本研究构造的深度学习网络,优化深度学习网络神经元的权重。具体设计的网络架构如图3所示。这里网络主要由四层特征编码层(Feature encoding layers)和四层特征解码层(Feature decoding layers)组成。
1.3.1深度学习网络中的特征编码层:
深度学习网络中特征编码层共四层,在图3中按网络运行次序从左至右分别为第1层、第2层、第3层、第4层,分别记为特征编码层1、特征编码层2、特征编码层3、特征编码层4。特征编码层运用下采样方法提取点云特征,每个特征编码层都包括采样模块、第一分组模块(简称分组模块)、第一PointConv模块(即为优化的PointConv模块,可简称PointConv模块)。特征编码第j层以一个(d+Cj)×Nj的点云矩阵作为输入,即输入具有d维坐标和Cj维特征的Nj个点。输出是一个(d+Cj+1)×Nj+1的矩阵,即输出具有d维坐标和Cj+1高维特征的总结局部上下文的Nj+1个子采样点所组成的二维矩阵。
特征编码层的采样模块:输入当前体素内的点集{p1,p2,…,pn},用最远点采样法从中选取一个子集
Figure BDA0003282409470000101
这样从中抽取出一些较为重要的中心点,进而减少点云样本数量并极大保留点云的原始的空间结构。采样模块的输入是一组大小为3×Nj的点的坐标,输出是大小为3×Nj+1的子采样点的坐标。
特征编码层的第一分组模块:采用最近邻规则分类法,在采样模块提取出的Nj+1个点中每个点的邻域内寻找最近的K个近邻点组成一个小组,结合这些近邻点对中心的Nj+1个点进行特征计算。分组模块的输入是采样模块输出的d×Nj+1子采样点的坐标和最初输入的大小为(d+Cj)×Nj的特征点集,输出两个结果,分别为d×K×Nj+1近邻点相对于每个采样中心点的坐标索引集合和(d+Cj)×K×Nj+1近邻点的特征集合,其中K是采样中心点邻域内近邻点的个数,在编码层中都设置为K=32。
特征编码层的第一PointConv模块,即优化的PointConv模块,解决了点云空间排列无序导致的卷积困难问题:利用密度函数对第一分组模块提取出的Nj+1个中心点附近K个近邻点的特征进行加权卷积,得到第j个特征编码层Nj+1个中心点的局部相关特征,并将该层的特征值送入第j+1个特征编码层继续。第一PointConv模块的输入是三部分数据,分别是Nj+1个中心点附近K个近邻点的特征grouped_feature,大小为(d+Cj)×K×Nj+1、Nj+1个中心点附近K个近邻点相对于中心点的局部坐标grouped_xyz,大小为d×K×Nj+1、Nj+1个中心点附近K个近邻点的密度density,大小为1×K×Nj+1。经过第一PointConv模块输出数据大小为(d+Cj+1)×Nj+1新的局部区域。
举例说明,图3中第1个特征编码层(图第一行左起第一个虚线框)所示,输入是(3)×4096的点云矩阵,其中(3)是点云的坐标维度d=3和特征维度C=0的和,4096是点云的总数。该输入首先进入采样模块,通过最远点采样法寻找N2=2048个点,得到3×2048采样中心点的坐标(sampled_xyz)。该采样中心点的坐标和原始输入的(3)×4096的特征点集再一同进入第一分组模块,第一分组模块在采样模块提取出的N2个中心点中每个点的邻域内寻找最近的K=32个近邻点组成一个小组,得到近邻点相对于每个采样中心点的坐标索引集合3×32×2048(grouped_xyz),并根据坐标索引得到近邻点的特征集合(3)×32×2048(grouped_feature),第1个特征解码层中坐标维度即为特征数。同时,第一分组模块得到的(3)×32×2048近邻点的坐标索引集合(grouped_xyz)通过核密度估计算法得到1×32×2048的近邻点密度(density),该近邻点密度再和近邻点的坐标索引集合(grouped_xyz)、近邻点的特征集合(grouped_feature)一起作为第一PointConv模块的输入,根据第一PointConv模块中多层感知机(MLP)等操作得到该特征编码层的最终输出,即第1个特征解码层2048个点的输出特征集合,为(64)×2048。
本研究构造的深度学习网络中特征编码层第二、三、四层的工作原理类似。通过编码层,不断提取体素内的中心点,使得其个数从4096减少到72,同时增加中心点包含的信息(中心点的特征维度从3增加到512)。经过四个特征编码层,实现点云的高维特征提取。
1.3.2深度学习网络中的特征解码层:
本文构造的深度学习网络中特征解码层的数量与特征编码层相同,也为4层,在图3中按网络运行次序从右至左分别为第4层、第3层、第2层、第1层,分别记为特征解码层4、特征解码层3、特征解码层2、特征解码层1。特征解码层运用上采样方法传递点云特征信息,每个特征解码层都包括特征插值模块、第二分组模块(简称分组模块)、第二PointConv模块(即优化的PointConv模块,可简称PointConv模块)。在特征解码第j层上,将经过特征解码后扩充的点云特征从(d+C′j)×Nj+1降低到(d+C′j-1)×Nj,而点云数量经过解码后逐步上升,其中Nj+1和Nj是特征解码第j层的输入和输出的点集大小。
特征插值模块:联合第j-1个特征编码层中采样模块得到的Nj个采样点
Figure BDA0003282409470000121
(lj-1_xyz)和第j个特征编码层中采样模块得到的Nj+1个采样点/>
Figure BDA0003282409470000122
(lj_xyz)开展特征关联操作(Feature association),然后再采用反距离加权插值的方法(见公式4)来实现特征传播。如以第j-1个特征编码层采样后得到的Nj个采样点中每个点为中心,在第j个特征编码层采样后得到的Nj+1个采样点中中选取距离Nj内每个点/>
Figure BDA0003282409470000123
最近的三个邻域点/>
Figure BDA0003282409470000124
得到三个邻域点/>
Figure BDA0003282409470000125
在Nj+1个采样点中的索引,并按照距离倒数计算出这三个邻域点/>
Figure BDA0003282409470000126
相对于Nj个点中对应点/>
Figure BDA0003282409470000127
的权重值/>
Figure BDA0003282409470000128
计算如式(3)所示,其中L2范数|| ||2表示两点之间的欧氏距离,μ为加权幂指数,默认值为2。
Figure BDA0003282409470000129
接着,通过这三个点的索引在第j个特征解码层最初输入的(d+C′j)×Nj+1特征点集中得到这三个邻域点对应的特征(d+C′j)×3,再结合对应的权重值
Figure BDA00032824094700001210
对这三个点的特征进行求和平均,得到/>
Figure BDA00032824094700001211
对应的(d+C′j)×1的插值特征/>
Figure BDA00032824094700001212
因为第j-1个特征编码层共有Nj个采样点,所以共得到(d+C′j)×Nj的插值特征(图3中下标为“插值结果”的小框所示),其中每个插值特征/>
Figure BDA0003282409470000131
具体计算如下式:
Figure BDA0003282409470000132
其中c′t,j∈C′j (4);
得到(d+C′j)×Nj的插值特征后,运用跳跃连接(Skip links)将该插值特征与第j-1个特征编码层的输出特征(d+Cj)×Nj(lj-1_points)进行连接操作,得到最终特征插值模块的输出结果(2d+Cj+C′j)×Nj(图3下方底色为灰色的小框所示)。
特征解码层的第二分组模块和特征编码层中的第一分组模块类似。输入是特征插值模块的最终结果(2d+Cj+C′j)×Nj和第j-1个特征编码层中采样模块得到的d×Nj采样中心点坐标(lj-1_xyz)。以这Nj个点中每个点为中心,在包括其自身的邻域内选取K′个近邻点,共(2d+Cj+C′j)个特征,扩充结果分别为(2d+Cj+C′j)×K′×Nj和d×K′×Nj,即作为第二分组模块的输出。其中K′是中心点邻域内近邻点的个数,在解码层中都设置为K′=16。
特征解码层的第二PointConv模块和特征编码层中的第一PointConv模块类似。第二PointConv模块的输入是三部分数据,分别是Nj个中心点附近K′个近邻点的特征grouped_feature,大小为(2d+Cj+C′j)×K′×Nj、Nj个中心点附近K′个近邻点相对于中心点的局部坐标grouped_xyz,大小为d×K′×Nj、Nj个中心点附近K′个近邻点的密度density,大小为1×K′×Nj。经过第二PointConv模块输出数据大小为(d+C′j-1)×Nj新的局部区域。
举例说明,如图3中第4个特征解码层所示,输入是(512)×72的特征点集,以第3个特征编码层得到的3×128的采样点(l3_xyz)中每个点为中心,在第4个特征编码层得到的3×72的采样点(l4_xyz)中选取3个距离其最近且不包括自身的点,获得这3个点在这72个采样点(l4_xyz)中的索引,并按照距离倒数计算出各自的权重值。然后通过这3个点的索引在第4个特征解码层最初输入的(512)×72特征点集(l4_input)中找到这3个点对应的特征(512)×3,再结合距离倒数计算得到的权重值对这3个点的特征进行求和平均,得到(512)×1的平均特征,因为第3个特征编码层共有128个采样点,所以得到(512)×128的插值特征。该插值特征通过跳跃连接(Skip links)与第3个特征编码层输出的(256)×128特征点集(l3_points)进行连接,得到(768)×128的特征点集。该特征点集再和第3个特征编码层得到的3×128的采样点坐标(l3_xyz)一同进入第二分组模块,第二分组模块在128个采样点中每个点的邻域内寻找最近的K′=16个近邻点组成一个小组,得到3×16×128近邻点相对于每个采样点的坐标索引集合(grouped_xyz),并根据坐标索引得到(768)×16×128近邻点的特征集合(grouped_feature)。同时,第二分组模块得到的3×16×128近邻点的坐标索引集合(grouped_xyz)通过密度估计算法得到1×16×128的点云密度(density),该点云密度再和点的坐标索引集合3×16×128(grouped_xyz)、点的特征集合(768)×16×128(grouped_feature)一起作为第二PointConv模块的输入,根据第二PointConv模块中多层感知机(MLP)等操作得到该特征解码层的最终输出,即(512)×128的第4个特征解码层128个点的输出结果(图3右边特征解码层4所示)。
深度网络中特征解码层第二、三、四层的工作原理类似。由此可见,四个特征解码层将特征从子采样点云传播到更密集的点云,即随着层数的加深,中心点的个数越来多(中心点的个数从72增加到4096),但每个中心点包含的特征信息越来越少(中心点的特征维度从512减少到128)。经过四个特征解码层,获得所有输入点的特征。
最后,全连接层通过1×1的卷积操将体素内128维的点云特征压缩至2维,即属于枝干和叶子类别的置信度,根据置信度大小实现林木逐点的枝叶分类。
1.3.3改进的PointConv模块(即第一PointConv模块和第二PointConv模块):
第一PointConv模块和第二PointConv模块计算原理相同。
传统的卷积是作用在二维图像数据上,图像数据通常可以表示成密集的有排列规律的网格形式,每个像素之间的相对位置总是有序的。带入进卷积神经网络中,每个滤波器被限制在一个小的局部区域内,按次序滑动窗口开展卷积操作。而三维点云数据是一组空间中无序的三维点,一组点云集合{pi|i=1,2,...,n},其中每个点包含一组位置向量(x,y,z)或其它信息如颜色、表面切线、表面切矢量等。不同于图像,点云的表示形式更灵活。因此,一个点云在不同的局部区域内的顺序和相对位置都是不同的。面向图像的传统离散卷积滤波器不能直接应用于点云,所以,本文提出了改进的PointConv网络开展局部与全局点云特征的提取。
图4显示了以采样点pi为中心点的K个近邻点所组成的局部区域上的操作,并展示了以逐点输入计算的结果。图4为改进的PointConv模块框架图。图4中(a)表示中心点及其近邻点的一个局部区域,其中pi为中心点,p1,p2,p3,p4为pi的近邻点,对应的f1,f2,f3,f4为点的特征;图4中(b)表示在以点pi为中心的局部区域上进行PointConv操作的过程。输入特征来自于以pi为中心的K个最近邻点的特征集合,输出为pi对应的特征信息。
改进的PointConv模块的输入分为三部分:1)grouped_feature是pi的K个近邻点的特征,大小为Cin×K,Cin为输入特征数。2)grouped_xyz是K个近邻点相对于中心点pi的局部坐标,大小为3×K;3)density是K个近邻点的密度,大小为1×K。PointConv模块的输出是中心点pi的特征,大小为Cout×1,Cout为输出特征数。
优化的PointConv模块具体运算步骤如下:
第一步:计算逆密度:利用核密度估计(KDE)方法开展K个近邻点grouped_xyz中每个输入点的核密度估计(即K个近邻点的密度,见公式5),再计算K个近邻点的逆密度变换(见公式8),再把结果输入MLP2进行激活函数的非线性变换,得到1×K大小的逆密度系数,接着对逆密度系数进行Cin次的复制平铺
Figure BDA0003282409470000151
输出大小为Cin×K的逆密度张量。
核密度估计(KDE)公式:
Figure BDA0003282409470000152
其中,
Figure BDA0003282409470000153
为密度函数,h为窗宽(Bandwidth),/>
Figure BDA0003282409470000154
表示对第k个点pk选取邻域点pm的总个数,/>
Figure BDA0003282409470000155
个点从输入集合中获取,选取的/>
Figure BDA0003282409470000156
Figure BDA0003282409470000157
为核函数,使用的核函数/>
Figure BDA0003282409470000158
主要有:
Figure BDA0003282409470000159
/>
其中B(a,b)=Γ(a)Γ(b)/Γ(a+b),伽玛函数为
Figure BDA00032824094700001510
当t=0时,为均匀核;当τ=1时,为Epanechnikov核;当τ=2时,为Biweight核。经过反复试验,这里采用了Biweight核。
最佳的窗宽h通过拇指法则来计算,其中的计算公式如下:
Figure BDA00032824094700001511
其中
Figure BDA00032824094700001512
是输入样本点与邻域点之间距离的标准差,/>
Figure BDA00032824094700001513
是核密度估计对象的维度,这里为1;逆密度变换S为:
Figure BDA00032824094700001514
第二步:计算权重:将K个近邻点的局部坐标grouped_xyz输入MLP1,输出大小为W=(Cin×Cout)×K的权重张量。其中1×1卷积的作用相当于多层感知机,在1×1卷积之后使用一个Relu的非线性层开展非线性变换。
第三步:把第一步的输出结果Cin×K的逆密度张量和K个近邻点的特征Fin即grouped_feature进行点乘运算(两个大小相同的矩阵对应元素相乘),输出大小为
Figure BDA0003282409470000161
的矩阵,再对/>
Figure BDA0003282409470000162
进行Cout次的复制平铺/>
Figure BDA0003282409470000163
输出大小为(Cin×Cout)×K的中间张量。
第四步:把第二步的输出结果(Cin×Cout)×K的权重张量和第三步的输出结果(Cin×Cout)×K的中间张量进行点乘运算,输出大小为
Figure BDA0003282409470000164
的矩阵。整体公式如式(9),其中·代表点乘。
Figure BDA0003282409470000165
第五步:对第四步的输出结果
Figure BDA0003282409470000166
的第一维和最后一维进行求和运算,得到中心点Pi的输出特征Fout=Cout×1。
1.3.4损失函数:
深度学习网络以负对数似然函数(The negative log likelihood lossfunction)作为网络的损失函数,如公式(10)、(11)所示。在训练过程中,损失函数定义如下:
Figure BDA0003282409470000167
Figure BDA0003282409470000168
公式(10)中l={1,2}表示枝叶两个类别。indici,l表示与分类类别有关的指标,如果计算出的类别与体素中第i个点pi的当前类别l相同,则将该指标赋值为1,否则赋值为0。
Figure BDA0003282409470000169
表示每个点pi通过深度学习网络反演得到的当前类别的概率。/>
Figure BDA00032824094700001610
表示通过Softmax函数得到的体素中pi点属于第l个类别的置信度。
1.4树木骨架重建:
面向枝叶分类后的枝干点云数据,运用先前提出的方法,对三块试验林地的橡胶树进行了枝干骨架重建,步骤如下:
首先,进行树干高度分层及每层中心点求取。具体为:根据单株树木的垂直高度,按一定的高度间隔将树干点云自下而上分为若干层。每棵树的分层树取决于树的高度和设定的高度间隔,层数
Figure BDA0003282409470000171
Δh统一设置为0.3m。在获得每一层的枝干点云后,按照DBSCAN密度聚类算法从每一层的枝干点云中提取聚类中心点
Figure BDA0003282409470000172
DBSCAN算法的距离阈值设置为0.10m,根据三个样地的林木平均点云分辨率0.02m以及多次实验判别得到。layer表示第layer层,来标记不同层次的枝干分布,以此获得整棵树木的枝干连接点雏形。
其次,对提取的中心点分类,分别为根结点、分枝结点和边缘结点。从根结点开始,根据最短距离与上层结点连接,将树木骨架转化为从根结点到每层中心点直至边缘结点结束的连接链。属于相邻层的每个中心点之间的连接链由带有空间方向性与自适应半径的圆柱体
Figure BDA0003282409470000173
进行拟合,r为圆柱体的半径,圆柱体的中心轴线根据对应的该段枝干扫描点云Pb拟合而成。具体步骤如下:/>
Figure BDA0003282409470000174
其中
Figure BDA0003282409470000175
代表组成该段枝干的/>
Figure BDA0003282409470000176
个点云,/>
Figure BDA0003282409470000177
为该段枝干所有点云的中心坐标,即点云坐标的平均值。SVD表示奇异值分解,分解后U为/>
Figure BDA0003282409470000178
的矩阵,∑为/>
Figure BDA0003282409470000179
的矩阵,V为3×3的矩阵。取V的第一列V(1)=[V(1,1),V(1,2),V(1,3)],即为所拟合直线的矢量,结合该段枝干点云的中心点/>
Figure BDA00032824094700001710
构造出拟合的空间直线L1方程,其点法式为:
Figure BDA00032824094700001711
接下来,根据公式(13)计算该段枝干所有点云
Figure BDA00032824094700001712
到拟合直线L1的最短距离di。其中,q1,q2分别为拟合直线L1上任意2个点,|| ||2代表2范数。
Figure BDA00032824094700001713
则圆柱体的半径设置为该段枝干所有点云到拟合直线的平均距离,即
Figure BDA00032824094700001714
再次,实现树木主枝干和一级枝干分类。具体为:将树木骨架划分为主枝干和许多一级分枝。其中,主枝干为从根结点开始,基于最短的Dijkstra距离算法,并根据树木生长角的最小变化规律寻找树干上的连接链。对于一级分枝是以主枝干链上的分枝结点一直搜索到另一端边缘结点所组成的链。
最后,根据植物生理学理论,将一个叶团簇定义为附着在同一主枝干或一级分枝上的叶子集合。利用空间三维分水岭分割的方法,以提取的主枝干和一级分枝为聚类中心对叶团簇进行分割。叶团簇分割之后,便可运用凸包进行单株树木叶团簇体积的计算。因此,树木的三维重建可以简化为不同枝干骨架和若干叶团簇组合的模型。
2.1深度学习网络的训练和测试结果:
深度学习网络的训练和测试均是在搭载了Intel(R)Core(TM)i7-10750H [email protected]处理器(Intel Inc.,Santa Clara,CA,USA)和16GB-RAM的Windows10-64位的服务器上进行的。由于深度学习需要大量的训练数据来训练网络,为了提高算力,使用NVIDIAGTX 1650Ti GPU(NVIDIA Inc.,SantaClara,CA,USA)代替CPU减少算法的训练时间。在本研究构建的深度学习网络中,学习率为0.0001,批处理规模Batch size为16,Epoch为200。训练精度和训练损失如图5所示,总的训练时间约为48小时。
图5为本文深度学习网络在不同轮数下的训练精度和损失函数曲线,图5中(a)为深度学习网络分类精度曲线示意图。图5中(b)为深度学习网络损失函数曲线示意图。
随着学习过程的不断进行,训练样本(每个体素中的点云)的训练精度呈上升趋势,训练损失呈下降趋势,说明本研究的深度学习网络是一个全局优化过程。训练准确率和训练损失在前25个Epoch分别有显著的增加和减少,表明深度网络中神经元的权值在快速达到枝叶分类任务的要求。在训练过程中,神经元网络在一个批次中可能会遇到一些复杂的样本,如包含多棵树的部分体素、遮挡产生的数据缺失或叶子和枝干形状模糊等,这些削弱了模型的学习效率,并导致回归损失函数值的剧烈波动。但曲线的整体上升趋势和下降趋势表明本文的深度网络在训练过程中有较好的收敛性。经过200个Epoch后,训练样本的精度和损失分别收敛到97.33%和0.05,说明本研究构建的深度学习网络具备点云枝叶分类的能力。
深度学习方法与传统的机器视觉算法检测的对比结果如表2所示,分别在分类精度、IOU与分类时间三个参量上阐述算法的性能。其中IOU定义如下:
Figure BDA0003282409470000181
上式中
Figure BDA0003282409470000182
和/>
Figure BDA0003282409470000183
分别代表每个点云pi所对应的真实标签/>
Figure BDA0003282409470000184
与算法识别的类别标号/>
Figure BDA0003282409470000185
==代表如果相等,则返回“True”,如果不相等,即≠,则记为“False”。通过式(14)给出了所有点云中分类正确和分类错误的数量比值。
测试样本并代入权值训练好的深度学习网络中开展枝叶分类,不同方法的分类结果如表2和图6中(a1)、(a2)、(a3)所示。传统的机器学习是基于特征描述和检测的方法,对于目标特征的提取依赖于人工设计与获取,整个算法无特征与权重的正反传播,算法的学习表达能力有限,仅获得了75.24%的分割准确率。对于树冠形状不规则且内部结构复杂的森林区域,森林的枝叶分布密度高和冠层之间的重叠遮蔽可会导致分割精度下降。而基于深度学习的三维点云识别通常采用多视角投影、体素剖分、原始点云直接输入等方法,可以依靠池化、卷积、激活与特征压缩和增广来自动提取三维点云有效特征,相较于传统的机器学习方法,深度学习方法直接使用三维无序点云作为输入,通过局部特征和全局特征进行串联的方式来聚合更多的特征信息,根据类别置信度实现点云分割。然而,PointNet通过多层感知机与最大池化操作只能够保留全局最优的特征信息,而丢失了局部基团所包含的点云局部特征,因此分割的准确度受到影响,分割准确率为82.75%。而本研究采用的PointConv方法,既考虑了全局特征的获取,又通过核密度函数对局部点云特征进行加权卷积与变换,并联合插值模块实现点云多尺度特征融合,从不同尺度上增强了对树体点云语义特征的刻画,因此获得更高的枝叶分割精度(90.32%)。表2列出了三种方法对树木点云进行枝叶分割的比较结果,从表中可以看出深度学习方法虽然花费的分割时间较长,但无论是在IoU数值还是整体分割准确率方面无疑都取得了更好的量化结果,可见,本研究的深度学习框架结合了树体点云的局部特征进一步提升了分割的准确性能。
表2为不同方法的点云枝叶分类检测结果:
方法 重叠度(枝/叶)/% 总体精度/% 分类时间(训练/测试)/h
传统机器学习方法 60.52/54.28 75.24 25.45/0.04
PointNet 70.16/66.53 82.75 43.94/0.06
本研究方法 78.24/72.74 90.32 48.14/0.06
2.2橡胶树木骨架重建结果:
针对深度学习提取的骨架数据,并根据1.4树木骨架重建的方法,对三块不同类型的橡胶树测试样地开展树木骨架重建。图6展现了三块橡胶树样地枝叶分离和树木骨架重建的结果。图6中(a1)为橡胶树品种PR107测试样地的激光点云数据以及枝叶分离结果示意图。图6中(a2)为橡胶树品种CATAS 7-20-59测试样地的激光点云数据以及枝叶分离结果示意图。图6中(a3)为橡胶树品种CATAS 8-79测试样地的激光点云数据以及枝叶分离结果图。图6中(b1)为橡胶树品种PR107测试样地枝干圆柱体拟合结果图。图6中(b2)为橡胶树品种CATAS 7-20-59测试样地枝干圆柱体拟合结果图。图6中(b3)为橡胶树品种CATAS 8-79测试样地枝干圆柱体拟合结果图。图6中(c1)为橡胶树品种PR107测试样地主枝干和一级枝干分类结果图。图6中(c2)为橡胶树品种CATAS 7-20-59测试样地主枝干和一级枝干分类结果图。图6中(c3)为橡胶树品种CATAS 8-79测试样地主枝干和一级枝干分类结果图。图6中(d1)为橡胶树品种PR107测试样地基于空间分水岭和一级分枝的树叶聚类结果与叶团簇的凸包显示图。图6中(d2)为橡胶树品种CATAS 7-20-59测试样地基于空间分水岭和一级分枝的树叶聚类结果与叶团簇的凸包显示图。图6中(d3)为橡胶树品种CATAS 8-79测试样地基于空间分水岭和一级分枝的树叶聚类结果与叶团簇的凸包显示图。
图6中(b1)、(b2)、(b3)显示了三块样地的圆柱体拟合结果。在拟合过程中,首先按照一定的高度间隔(0.3m)将枝干点云自下而上划分为多层,对每层点云进行聚类提取中心点并根据最短的Dijkstra距离建立连通链;接着采用带方向性与自适应半径的圆柱体开展枝干的分段拟合,拟合后平均每棵橡胶树枝干约由180个圆柱体组合而成,且每个圆柱体被标注为不同的颜色。
其次,根据生长角度变化与中心点的连通链将树木骨架划分为主枝干和许多一级分枝,如图6中(c1)、(c2)、(c3)所示,提取的主枝干,从主枝干延伸出的各个一级分枝通过不同深度的颜色表示。可以看出,PR107品种橡胶树枝干生长角度的变化较大,枝干向外扩散,树木骨架呈发散状。相较于PR107品种橡胶树,CATAS 7-20-59品种和CATAS 8-79品种橡胶树枝干生长角度的变化较小,枝干多向上方延伸,树木骨架呈现花瓶状。
之后,根据植物生理学理论,将叶团簇定义为在同一主枝干或一级分枝上生长的叶子集合。利用空间三维分水岭分割的方法,以提取的主枝干和一级分枝对应的点云为聚类中心,对每棵树的叶团簇进行分割,并将分割后的叶团簇用不同颜色表示,分割结果如图6(d1)、(d2)、(d3)所示。三块橡胶树样地中,PR107品种橡胶树一级分枝较多,因此每棵树分割的叶团簇为3个左右;而CATAS 7-20-59和CATAS 8-79品种橡胶树一级分枝较少,得到平均每棵树2.5个左右的叶团簇。最后,运用凸包算法对每株树木叶团簇进行标记并开展体积的计算。
2.3橡胶树参数反演结果分析:
表3给出了PR107、CATAS 7-20-59和CATAS 8-79三块橡胶树测试样地树木的生长参数实际测量值和本研究计算值,其中包括胸径、一级枝干直径和主枝干与一级枝干之间夹角,同时,通过三个比对指标R2,均方根误差(RMSE)和相对均方根误差(rRMSE)定量化来分析本研究方法的有效性。图7给出了具体的参数对比结果。
表3为本研究方法获得的林分参数与实地测量值对比:
Figure BDA0003282409470000211
注:(F)代表实际测量值,(O)代表本研究方法的计算值。
图7中(a)显示了本研究方法获取的三块橡胶树样地树木胸径实测值与计算值的比较结果。其中PR107品种橡胶树和CATAS 7-20-59品种橡胶树与实地观测的比对结果分别为(R2=0.95,RMSE=3.01cm,rRMSE=12.22%)和(R2=0.94,RMSE=2.78cm,rRMSE=12.53%)。相对于前两块样地而言,算法在CATAS 8-79品种橡胶树枝干参数提取上取得了更好的效果(R2=0.97,RMSE=2.63cm,rRMSE=10.55%)。这主要因为PR107和CATAS 7-20-59品种橡胶树是在有叶时期采集的点云数据,树木不同器官间相互遮挡会导致数据缺失,且扫描过程中风力扰动使树体抖动产生噪点,进而影响了最终的参数反演结果。而CATAS8-79品种橡胶树不耐寒害,数据采集时处于落叶状态,因此获取的枝干数据较为完整,且CATAS 8-79品种橡胶树枝干较粗,风力扰动带来的树体的震动影响较小,使得获取的点云质量较高、噪声少,参数反演的精度更高。
图7中(b)显示了三块橡胶树样地计算的树木主枝干与一级枝干夹角与实测值的比对结果,可以看出对不同品种的橡胶树都获得了较好的精度。其中PR107品种(R2=0.94,RMSE=4.94°,rRMSE=11.23%)和CATAS 8-79品种(R2=0.92,RMSE=3.48°,rRMSE=10.86%)的相关系数相对于CATAS 7-20-59品种(R2=0.91,RMSE=3.81°,rRMSE=11.45%)略高一些,可以解释为PR107品种橡胶树具有较为清晰的枝干特征和发散型的骨架结构,一级枝干和主枝干之间的夹角较大,在实地测量和算法角度计算中更易获得准确的结果;而CATAS 8-79品种橡胶树林内叶子稀疏导致受遮挡影响小,使得分枝角估算较准确。
图7中(c)显示了三块橡胶树样地树木一级枝干宽度和对应的叶团簇体积的相关性。其中,PR107品种橡胶树由于其较大的分枝夹角为叶片生长提供了更多空间,使得橡胶树具有较大的冠积和一级枝干对应的叶团簇的体积。CATAS 7-20-59品种橡胶树的分枝角度范围小,形成了花瓶状的树冠,因此叶团簇的大小中等且分布均匀。而CATAS 8-79品种橡胶树因为出现了寒害落叶现象,一级枝干对应的叶团簇体积较小且在图中分布较集中。图7中(c)显示一级枝干宽度与对应的叶团簇体积之间总体满足正相关分布,即较粗的枝干支撑起较大的叶团簇,这与植物枝干负责营养传输和叶花果等器官的重力承载原理相符,所产生的局部偏差可解释为人为修剪或遭受的自然环境扰动而带来的局部变化。
3结论:
树木骨架的精准模型重建对于分析树木的表型结构、树体自身的趋光性、林木竞争因子的影响都起到了重要的数据支撑作用。首先,本研究提出了改进的PointConv深度学习算法,面向橡胶林的背包移动激光雷达数据,通过对点云开展体素剖分以及设计四层的编码与解码层,和改进的PointConv模块与插值特征传播模块,在不同尺度下自动获取林木点云的全局最优与局部邻域的核密度特征,旨在精准实现枝叶点云的分类。通过对三个不同品种的橡胶林样地测试数据分析表明,本研究的深度学习方法分类准确率达到90.32%;其次,面向提取出的枝干点云,运用计算机图形学算法,自动提取主枝干胸径、一级分枝直径、分枝夹角、叶团簇体积等林学参数;最后,通过与样地实测值比对,不同品种橡胶树主枝干、一级分枝的直径计算值与真实值的相关性R2均高于0.94,RMSE低于3.01cm;主枝干与一级枝干夹角估测值与真实值的相关性R2高于0.91,RMSE低于4.94°,且橡胶树的叶团簇与对应的一级分枝直径呈现近似的正相关分布。在对林木的生态研究中,同科属不同品种的树体间的表观差异性、所表现的空间异速生长规律和与自然环境间产生的生产驱动力分析等都可从树体的表观形态结构上找寻到规律。本文提出的深度学习网络与计算机图形学算法相结合,能够自动从林木扫描数据中重构树木的骨架空间模型,并适用于对不同品种的橡胶树进行参数反演,是现今人工智能技术在林业领域落地发展的应用点之一。
本发明的保护范围包括但不限于以上实施方式,本发明的保护范围以权利要求书为准,任何对本技术做出的本领域的技术人员容易想到的替换、变形、改进均落入本发明的保护范围。

Claims (10)

1.基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:包括:
步骤1:获取林木激光点云数据;
步骤2:对采集的激光点云数据进行去燥处理,将去燥后的激光点云数据分为地上点和地面点;
步骤3:采用机器学习算法以及人工修正方式对地上点进行枝叶分离操作,对已进行枝叶分离操作的地上点进行体素化剖分;
步骤4:将已经进行枝叶分离操作和体素化剖分的点云数据作为训练样本数据集,并采用数据增广的方法对训练样本数据集进行扩增,得到新的训练样本数据集;
步骤5:构建深度学习网络,采用新的训练样本数据集对深度学习网络开展训练,得到训练好的深度学习网络模型;
步骤6:采集待测林地的激光点云数据,对采集的激光点云数据进行去燥处理,将去燥后的激光点云数据分为地上点和地面点;对地上点进行体素化剖分,将单个体素内的点云数据输入到训练好的深度学习网络模型中,以实现体素内的点云数据的枝叶分离。
2.根据权利要求1所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:所述深度学习网络包括特征编码层和特征解码层;所述特征编码层用于采用下采样方法提取点云特征信息,所述特征解码层用于采用上采样方法传递点云特征信息;
所述特征编码层包括采样模块、第一分组模块和第一PointConv模块;所述特征解码层包括特征插值模块、第二分组模块和第二PointConv模块。
3.根据权利要求2所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:所述特征编码层共四个,所述特征解码层共四个。
4.根据权利要求3所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:第j个特征编码层中采样模块的计算过程为:
输入大小为(d+Cj)×Nj的点云矩阵,通过最远点采样法从该点云矩阵中选取Nj+1个子采样点,得到大小为(d+Cj+1)×Nj+1的点云矩阵;
其中:(d+Cj)×Nj是指具有d维坐标和Cj维特征的Nj个点,(d+Cj+1)×Nj+1是指具有d维坐标和Cj+1高维特征的Nj+1个子采样点。
5.根据权利要求4所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:第j个特征编码层中第一分组模块的计算过程为:
输入采样模块输出的d×Nj+1的子采样点的坐标和最初输入的(d+Cj)×Nj的特征点集,采用最近邻规则分类法,在采样模块提取出的Nj+1个采样点中每个点的邻域内寻找最近的K个近邻点组成一个小组,得到d×K×Nj+1近邻点相对于每个采样中心点的坐标索引集合,并根据坐标索引得到(d+Cj)×K×Nj+1近邻点的特征集合,其中K是Nj+1个采样中心点中每个点的邻域内近邻点的个数。
6.根据权利要求5所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:第j个特征编码层中第一PointConv模块的输入是三部分数据;其中第一部分数据为:Nj+1个中心点附近K个近邻点的特征,大小为(d+Cj)×K×Nj+1;第二部分数据为:Nj+1个中心点附近K个近邻点相对于中心点的局部坐标,大小为d×K×Nj+1;第三部分数据为:Nj+1个中心点附近K个近邻点的密度,大小为1×K×Nj+1;第一PointConv模块输出大小为(d+Cj+1)×Nj+1的新的局部区域。
7.根据权利要求6所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:所述第j个特征编码层中第一PointConv模块的计算过程为:
步骤(1)、利用核密度估计方法计算中心点
Figure FDA0003282409460000021
附近K个近邻点中每个输入点的核密度估计,再通过核密度估计计算K个近邻点的逆密度变换,再逆密度变换的结果输入多层感知机MLP2进行激活函数的非线性变换,得到大小为1×K的逆密度系数,对逆密度系数进行Cin次的复制平铺,输出大小为Cin×K的逆密度张量;
步骤(2)、将K个近邻点相对于中心点的局部坐标输入多层感知机MLP1,输出大小为(Cin×Cout)×K的权重张量W,其中Cin为输入特征数,Cout为输出特征数;
步骤(3)、把步骤(1)输出的大小为Cin×K的逆密度张量和K个近邻点的特征Fin进行点乘运算,输出大小为Cin×K的矩阵
Figure FDA0003282409460000022
再对/>
Figure FDA0003282409460000023
进行Cout次的复制平铺,输出大小为(Cin×Cout)×K的中间张量;
步骤(4)、把步骤(2)输出的大小为(Cin×Cout)×K的权重张量W和步骤(3)输出的大小为(Cin×Cout)×K的中间张量进行点乘运算,输出大小为(Cin×Cout)×K的矩阵
Figure FDA0003282409460000024
步骤(5)、对步骤(4)得到的矩阵
Figure FDA0003282409460000031
中的第一维和最后一维进行求和运算,得到中心点
Figure FDA0003282409460000032
的输出特征Fout
其中中心点
Figure FDA0003282409460000033
为Nj+1个中心点中的某一点。
8.根据权利要求7所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:
第j个特征解码层中特征插值模块的计算过程为:
以第j-1个特征编码层中采样模块得到的Nj个采样点中每个点
Figure FDA0003282409460000034
为中心,在第j个特征编码层中采样模块得到的Nj+1个采样点中,选取距离Nj个采样点中每个点/>
Figure FDA0003282409460000035
最近的三个邻域点/>
Figure FDA0003282409460000036
得到三个邻域点/>
Figure FDA0003282409460000037
在Nj+1个采样点中的索引,其中t=1,2,3;
按照距离倒数计算出三个邻域点
Figure FDA0003282409460000038
相对于Nj个采样点中对应中心点/>
Figure FDA0003282409460000039
的权重值
Figure FDA00032824094600000310
计算公式为:
Figure FDA00032824094600000311
其中
Figure FDA00032824094600000312
表示点/>
Figure FDA00032824094600000313
与点/>
Figure FDA00032824094600000314
之间的欧氏距离,μ表示加权幂指数;
通过这三个邻域点的索引在第j个特征解码层最初输入的(d+C′j)×Nj+1特征点集中得到这三个邻域点对应的特征(d+C′j)×3,再结合对应的权重值
Figure FDA00032824094600000315
对这三个邻域点的特征进行求和平均,得到/>
Figure FDA00032824094600000316
对应的(d+C′j)×1的插值特征/>
Figure FDA00032824094600000317
其中插值特征/>
Figure FDA00032824094600000318
的计算公式为:
Figure FDA00032824094600000319
/>
因为第j-1个特征编码层共有Nj个采样点,所以得到了(d+C′j)×Nj的插值特征;
得到(d+C′j)×Nj的插值特征后,运用跳跃连接方法将(d+C′j)×Nj插值特征与第j-1个特征编码层的输出的(d+Cj)×Nj的特征进行连接操作,得到特征插值模块的最终输出结果:(2d+Cj+C′j)×Nj
其中(d+C′j)×Nj表示具有d维坐标和C′j维特征的Nj个点。
9.根据权利要求8所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:
第j个特征解码层中第二分组模块的计算过程为:
输入特征插值模块的最终输出结果(2d+Cj+C′j)×Nj和第j-1个特征编码层中采样模块得到的d×Nj的子采样点的坐标,采用最近邻规则分类法,在Nj个采样点中每个点的邻域内寻找最近的K′个近邻点组成一个小组,得到d×K′×Nj近邻点相对于每个采样中心点的坐标索引集合,并根据坐标索引得到(2d+Cj+C′j)×K′×Nj近邻点的特征集合,其中K′是Nj个采样中心点中每个点的邻域内近邻点的个数。
10.根据权利要求9所述的基于深度学习方法的林木激光点云枝叶分离方法,其特征在于:
第j个特征解码层中第二PointConv模块的输入是三部分数据;其中第一部分数据为:Nj个中心点附近K′个近邻点的特征,大小为(2d+Cj+C′j)×K′×Nj;第二部分数据为:Nj个中心点附近K′个近邻点相对于中心点的局部坐标,大小为d×K′×Nj;第三部分数据为:Nj个中心点附近K′个近邻点的密度,大小为1×K′×Nj;第二PointConv模块输出大小为(d+C′j-1)×Nj新的局部区域。
CN202111137023.3A 2021-09-27 2021-09-27 基于深度学习方法的林木激光点云枝叶分离方法 Pending CN115880487A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111137023.3A CN115880487A (zh) 2021-09-27 2021-09-27 基于深度学习方法的林木激光点云枝叶分离方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111137023.3A CN115880487A (zh) 2021-09-27 2021-09-27 基于深度学习方法的林木激光点云枝叶分离方法

Publications (1)

Publication Number Publication Date
CN115880487A true CN115880487A (zh) 2023-03-31

Family

ID=85763021

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111137023.3A Pending CN115880487A (zh) 2021-09-27 2021-09-27 基于深度学习方法的林木激光点云枝叶分离方法

Country Status (1)

Country Link
CN (1) CN115880487A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116167668A (zh) * 2023-04-26 2023-05-26 山东金至尊装饰工程有限公司 基于bim的绿色节能建筑施工质量评价方法及***
CN116486261A (zh) * 2023-04-14 2023-07-25 中山大学 一种树木点云木质成分与叶片成分的分离方法及***
CN116893428A (zh) * 2023-09-11 2023-10-17 山东省地质测绘院 基于激光点云的森林资源调查与监测方法及***
CN117710601A (zh) * 2023-12-27 2024-03-15 南京林业大学 一种基于激光点云和图像信息的单木骨架提取方法及***

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116486261A (zh) * 2023-04-14 2023-07-25 中山大学 一种树木点云木质成分与叶片成分的分离方法及***
CN116167668A (zh) * 2023-04-26 2023-05-26 山东金至尊装饰工程有限公司 基于bim的绿色节能建筑施工质量评价方法及***
CN116893428A (zh) * 2023-09-11 2023-10-17 山东省地质测绘院 基于激光点云的森林资源调查与监测方法及***
CN116893428B (zh) * 2023-09-11 2023-12-08 山东省地质测绘院 基于激光点云的森林资源调查与监测方法及***
CN117710601A (zh) * 2023-12-27 2024-03-15 南京林业大学 一种基于激光点云和图像信息的单木骨架提取方法及***
CN117710601B (zh) * 2023-12-27 2024-05-24 南京林业大学 一种基于激光点云和图像信息的单木骨架提取方法及***

Similar Documents

Publication Publication Date Title
Li et al. Automatic organ-level point cloud segmentation of maize shoots by integrating high-throughput data acquisition and deep learning
Li et al. PlantNet: A dual-function point cloud segmentation network for multiple plant species
CN112381861B (zh) 一种基于地基激光雷达的林地点云数据配准和分割方法
CN111898688B (zh) 一种基于三维深度学习的机载LiDAR数据树种分类方法
Koukoulas et al. Mapping individual tree location, height and species in broadleaved deciduous forest using airborne LIDAR and multi‐spectral remotely sensed data
CN115880487A (zh) 基于深度学习方法的林木激光点云枝叶分离方法
CN111340826B (zh) 基于超像素与拓扑特征的航拍图像单株树冠分割算法
CN111598045B (zh) 一种基于对象图谱和混合光谱的遥感耕地变化检测方法
CN113591766B (zh) 一种无人机多源遥感的树种识别方法
CN109146889A (zh) 一种基于高分辨率遥感图像的农田边界提取方法
CN112819830A (zh) 基于深度学习与机载激光点云的单株树冠分割方法
CN110309780A (zh) 基于bfd-iga-svm模型的高分辨率影像房屋信息快速监督识别
CN111652193A (zh) 基于多源影像的湿地分类方法
Chen et al. Predicting individual apple tree yield using UAV multi-source remote sensing data and ensemble learning
CN113033453A (zh) 一种适用于景观破碎区作物类型遥感识别的方法及***
CN113269825B (zh) 基于地基激光雷达技术林木胸径值提取的方法
CN105447274A (zh) 一种利用面向对象分类技术对中等分辨率遥感图像进行滨海湿地制图的方法
CN117409339A (zh) 一种用于空地协同的无人机作物状态视觉识别方法
Fan et al. UAV image crop classification based on deep learning with spatial and spectral features
CN117114147A (zh) 基于雷达和卫星遥感的森林植被碳储量估算的方法和装置
CN107705344A (zh) 激光扫描环境点云数据中树冠模型提取方法
Chen et al. 3D model construction and ecological environment investigation on a regional scale using UAV remote sensing
Aviña-Hernández et al. Predictive performance of random forest on the identification of mangrove species in arid environments
CN114494586B (zh) 晶格投影的深度学习网络阔叶树枝叶分离与骨架重建方法
Zhu et al. Research on deep learning individual tree segmentation method coupling RetinaNet and point cloud clustering

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination