CN117765006A - 基于无人机影像与激光点云的多层次密集树冠分割方法 - Google Patents
基于无人机影像与激光点云的多层次密集树冠分割方法 Download PDFInfo
- Publication number
- CN117765006A CN117765006A CN202311791699.3A CN202311791699A CN117765006A CN 117765006 A CN117765006 A CN 117765006A CN 202311791699 A CN202311791699 A CN 202311791699A CN 117765006 A CN117765006 A CN 117765006A
- Authority
- CN
- China
- Prior art keywords
- region
- segmentation
- point
- block
- cluster
- 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
Links
- 230000011218 segmentation Effects 0.000 title claims abstract description 162
- 238000000034 method Methods 0.000 title claims abstract description 86
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 111
- 239000000203 mixture Substances 0.000 claims abstract description 13
- 238000000926 separation method Methods 0.000 claims abstract description 5
- 238000009826 distribution Methods 0.000 claims description 15
- 238000001914 filtration Methods 0.000 claims description 11
- 238000009499 grossing Methods 0.000 claims description 4
- 230000002146 bilateral effect Effects 0.000 claims description 3
- 238000011160 research Methods 0.000 abstract description 16
- 230000006870 function Effects 0.000 description 42
- 238000010586 diagram Methods 0.000 description 30
- 238000000111 isothermal titration calorimetry Methods 0.000 description 14
- 238000001514 detection method Methods 0.000 description 12
- 230000000694 effects Effects 0.000 description 10
- 238000012545 processing Methods 0.000 description 10
- 238000005259 measurement Methods 0.000 description 9
- 230000012010 growth Effects 0.000 description 6
- 241000124033 Salix Species 0.000 description 5
- 238000005520 cutting process Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 239000002023 wood Substances 0.000 description 5
- 238000011161 development Methods 0.000 description 4
- 230000018109 developmental process Effects 0.000 description 4
- 238000000605 extraction Methods 0.000 description 4
- 238000003709 image segmentation Methods 0.000 description 4
- 239000002028 Biomass Substances 0.000 description 3
- 238000013459 approach Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 238000013135 deep learning Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 230000001788 irregular Effects 0.000 description 3
- 238000005192 partition Methods 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- 241000218631 Coniferophyta Species 0.000 description 2
- 241000218998 Salicaceae Species 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 239000003086 colorant Substances 0.000 description 2
- 238000013527 convolutional neural network Methods 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000012417 linear regression Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000010606 normalization Methods 0.000 description 2
- 238000012706 support-vector machine Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000012876 topography Methods 0.000 description 2
- RZVHIXYEVGDQDX-UHFFFAOYSA-N 9,10-anthraquinone Chemical compound C1=CC=C2C(=O)C3=CC=CC=C3C(=O)C2=C1 RZVHIXYEVGDQDX-UHFFFAOYSA-N 0.000 description 1
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 description 1
- 244000025254 Cannabis sativa Species 0.000 description 1
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 241000112525 Salix psammophila Species 0.000 description 1
- 241000422846 Sequoiadendron giganteum Species 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 238000002790 cross-validation Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000003475 lamination Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000011430 maximum method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000008635 plant growth Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 238000011524 similarity measure Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- Y—GENERAL 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
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Image Analysis (AREA)
Abstract
本发明公开基于无人机影像与激光点云的多层次密集树冠分割方法，包括：采集林分点云数据和林分正射影像数据；对林分正射影像进行滤波，对林分点云数据进行滤波：对滤波后的点云归一化处理，用平面栅格法去除地面点云；使用改进超像素算法和能量函数对滤波后的正射影像进行分割；提取分割后的掩膜坐标，将二维图像的掩膜坐标与归一化后的点云坐标对齐进而完成三维分割，采用基于长轴的核密度估计算法确定每个初始分割区域内树冠顶点的个数；采用高斯混合模型进行点云分割，获得树冠分离结果。本发明通过无人机点云和影像数据相结合，能够精确检测密集林分中单个树冠，是评估森林生态系统的先决条件的科学研究方法，为精准林业提供了重要的推动力。
Description
技术领域
本发明涉及林业技术领域，具体是一种基于无人机影像与激光点云的多层次密集树冠分割方法。
背景技术
树木是森林的基本单位，其空间结构、生物特性和组成成分是森林生态环境建模、生物量估计、碳储量估计和资源整合的关键因素。单个树木的检测是集成遥感数据的森林库存的理想产品，单棵树冠分割(ITCs)是精准林业的重要组成部分，一旦准确地分割了一棵树，就可以得到树的高度、种类、冠大小、木材体积和生物量等属性。随着工业的发展和森林的砍伐，森林保护和可持续发展已成为当今社会的重要课题，与天然林相比，人工林具有生长快、生长量高、开发方便等特点，因此以人工林代替自然森林资源是保护森林的重要方式。传统的人工测量在林业领域虽然应用很广，但是需要消耗大量的时间周期和人力资源，且实时性较差，也无法获得大面积区域的林地数据。
遥感技术的快速发展为森林资源的监测创造了有利的条件，由于卫星图像分辨率低，导致树冠在分割正确率上有很大的限制。而无人机遥感因其能在低海拔飞行并获得的高分辨率数据已经成为林业数据获取的主流工具，通过携带RGB传感器获得的高分辨率二维影像具有丰富的形状和纹理信息，利用精细的空间分辨率图像对高郁闭度林分分割具有挑战性，因为单个像素覆盖了像素的不同组成成分，尤其对单个树冠的分割的复杂性而言，边缘分割通常为像素级别。传统的单个树冠分割对树冠大小和形状的均匀性进行假设，但模糊了树冠之间的边界，在多层次林地中，传统深度学习方法如图割法通过调整参数来分割树冠，识别在不同冠层大小的不同地点工作的设置是一个迭代的、耗时的过程，且结果需要视觉评估和微调，其泛化性也不足，近年来，CNNs的深度学习方法被用于植被遥感，该方法不需要使用图像分割和机器学习分类器的特征选择，但是需要标注大量的数据，且泛化性能差，在换一个林分结构的数据后模型性能会下降，此外对林地树冠分割其精度较低。
随着精准的林业的需求，二维图像获取林分的冠层信息的方法已经无法满足要求。虽然2.5维深度相机RGBD在室内建模和小面积信息获取方面得到广泛应用，但林地资源整合需要大面积数据，因此，人们把对林地信息获取转向了更高维度即3D无人机雷达，虽然该点云数据具有较好的空间信息和三维结构，且数据获取不受环境影响，但对于林地的机载激光雷达生成的点云进行分割是一项很大的挑战，现有的方法大致可以分为两类：
第一种是基于点云数据的直接分割，基于卷积神经网络分割算法例如使用PointNet等层出不穷，这些算法在在稀疏林地单棵树冠的检测和分割上表现出了很好的效果。在密集林地和复杂林分中传统深度学习算法更有效，结合树树冠表型和结构几何特点，例如通过随机样本共识(RANSAC)算法，先在点云数据的一个子集拟合一个模型，然后用该模型对剩余的点进行分类。这个过程不断重复，直到达到令人满意的分割效果；有基于轮廓模型的分割算法，该方法通过在点云数据上最小化能量函数实现分割。能量函数是基于相邻点之间的相似性而定义的，能够将点云数据分割成不同的类别；有通过分水岭算法，将点云数据作为地形表面，并进行洪水模拟，然后根据淹没区域的边界对点云数据进行分割；还有一些其他方法，包括区域增长算法、体素的点云分析算法、K-means算法和支持向量机(SVM)算法等。这些方法根据简单的接近规则和可信的树形状来来识别点并将点归类为单颗树，重点在于点之间的关系而不是形状，而且在处理三维点云时，效率低，与基于图像处理的方法相比需要更多的计算能力。
另一种方法是基于两步分割的策略，首先，将点云转换为数字表面模型(DSM)或冠层盖度模型(CHM)映射到平面光栅中，构成二维数字形态来描述森林，然后结合计算机视觉和图像处理实现单棵树冠分割，该方法的分割部分主要有基于梯度的方法，如Li等人在利用LiDAR数据得到的DSM上使用了一种分水岭算法来分割针叶林中的单个树冠，有学者结合了区域生长和流域算法，利用分水岭分割和基于形状的建模的组方法，从高密度LiDAR数据中分割单个树冠。当数据较为复杂时，譬如被分割区域边界重叠，则需要需先确定树的位置再对边界进行分割，一种常见的方法是局部搜索策略的区域生长算法通过寻找种子点确定树的位置，还有固定核带宽平均位移算法设置内核带宽聚合像素的颜色信息，从复杂的数据中分割单个树冠。这些算法的核心理念多基于计算计视觉，通过有效处理点云数据完成更好的分割。
尽管已经提出了许多方法对林地的激光点云数据进行树顶检测和树冠分割，但是针对不同类型林地这些方法都存在一定缺陷。大多数人工林的冠层分割都有一定的困难，其难点主要在于：(1)树冠之间交叉重叠，使得边界点难以界定，而且树冠形状不规则增加了分割的复杂性；(2)无人机雷达提供了自上而下的扫描视角，其数据可能受到场景中其他物体遮挡的影响，从而降低分割的准确性。(3)可变的树木密度导致分割复杂度提升，树木的密度在整个场景中可能会存在差异，从而影响分割算法的性能。因此，分割算法需要具备适应性，以处理不同密度的树木。为了满足精准林业的需求，解决密集林分条件下的ITCs分割困难、总体分割精度较低等问题，需要提出一种更精确的方法，可以实现单个树冠的精确分割，能够检测和测量复杂森林中单个树冠。
发明内容
本发明所要解决的技术问题是针对密集林分条件下的ITCs分割困难、总体分割精度较低等问题提供一种基于无人机影像与激光点云的多层次密集树冠分割方法，本基于无人机影像与激光点云的多层次密集树冠分割方法通过无人机点云和影像数据相结合，能够精确检测密集林分中单个树冠，是评估森林生态系统的先决条件的科学研究方法，并为精准林业提供了重要的推动力。
为实现上述技术目的，本发明采取的技术方案为：
一种基于无人机影像与激光点云的多层次密集树冠分割方法，包括：
步骤1：采集林分点云数据和林分正射影像数据；
步骤2：通过Voronoi划分对林分正射影像进行滤波，通过双边滤波对林分点云数据进行滤波；
步骤3：对滤波后的点云归一化处理，并用平面栅格法去除地面点云；
步骤4：使用改进超像素算法和能量函数对滤波后的正射影像进行分割；
步骤5：提取分割后的掩膜坐标，将二维图像的掩膜坐标与归一化后的点云坐标对齐进而完成三维分割，得到多个初始分割段；
步骤6：采用核密度估计算法确定每个初始分割段内树冠顶点的个数；
步骤7：基于树冠顶点的个数并采用高斯混合模型进行点云分割，获得树冠分离结果。
作为本发明进一步改进的技术方案，所述的步骤4具体为：
步骤4.1、从滤波后的正射影像中提取颜色距离分量、空间距离分量以及纹理相似度分量；
步骤4.2、基于颜色距离分量、空间距离分量以及纹理相似度分量并采用改进超像素算法将正射影像分割成多个区域块；
步骤4.3、基于能量函数合并区域块；
步骤4.3.1、计算每个区域块的能量值，其中，第j个区域块的区域代表点Pj的局部平均能量值简称为区域块j的能量值，区域块j的能量值为：
其中，表示区域块j内像素点r″的能量值；
步骤4.3.2、选择聚类中心：
其中，是以区域块j为聚类中心的区域块j的能量值，其中以区域块j为聚类中心的区域代表点Pj记为聚类中心区域代表点Pj，是所有区域块的能量值的均值，是同一个簇中区域块的区域代表点距聚类中心区域代表点Pj的距离，是所有的期望，表示同一簇中区域块的区域代表点距聚类中心区域代表点Pj的最远距离；
满足公式(2)和公式(3)的区域块j为聚类中心；
初始判定为同一簇的条件：
以最高能量值的区域块的区域代表点Pj为中心，八个方位进行搜索，在八方位中搜索距离中心最远的方位为基准，当位于最远方位的区域块ni的区域代表点的局部平均能量值下降到一定值ε或当位于最远方位的区域块ni的区域代表点的局部平均能量值不再继续下降时，则以该区域块ni的区域代表点和最高能量值的区域块的区域代表点Pj的距离为半径，区域代表点Pj为圆心画圆，该圆范围内涉及到的区域块归属于同一个簇；
优化簇：
以区域块j为聚类中心，当区域块bi与区域块j的能量相似度值大于等于同一簇中位于最远方位的区域块ni与区域块j的能量相似度值时，区域块bi归属于以区域块j为聚类中心的簇；
以区域块j为聚类中心，当区域块bi与区域块j的能量相似度值小于同一簇中位于最远方位的区域块ni与区域块j的能量相似度值时，区域块bi不属于以区域块j为聚类中心的簇；
其中区域块j与区域块bi的能量相似度值为： 表示区域块bi的能量值；区域块ni与区域块j的能量相似度值为： 表示区域块ni的能量值，K1为与聚类中心区域块j同属于同一个簇内的所有区域块的集合；
步骤4.3.3、计算未归属于簇的剩余区域块的能量值，按照步骤4.3.2的方法在剩余区域块中再选取聚类中心，判定归属于同一个簇的区域块，并优化该簇；
步骤4.3.4、将未归属于簇的剩余区域块按照步骤4.3.3的方法再选取聚类中心、判定归属于同一个簇的区域块，并优化该簇，直到所有区域块均归属于对应的簇为止；
步骤4.3.5、若一个区域块共同归属于两个或者两个以上不同的簇，则计算该区域块分别与共同归属的不同簇内聚类中心区域块的能量相似度值并比较，该区域块归属于其中能量相似度值最高所对应的聚类中心区域块所归属的簇；
步骤4.3.6、合并归属于同一个簇的区域块。
作为本发明进一步改进的技术方案，所述步骤6中，核密度估计算法的公式为：
其中，n′为每个初始分割中的点数，h>0为一个平滑参数，被称作带宽，K(x′)为核函数，x′表示样本点，x′γ表示观测点。
作为本发明进一步改进的技术方案，所述的步骤7中，采用高斯混合模型进行点云分割，其公式为：
其中，N′表示核密度估计的树冠点数，X表示投影在剖面上的点，λk′为比例系数，表示各混合分量的先验概率，μk′、δk′表示高斯分布初始参数，分别表示高斯函数的平均值和方差。
本发明的有益效果为：
本发明针对密集林分条件下的ITCs分割困难、总体分割精度较低等问题，结合无人机正射影像和点云数据的特点开发了一种高精度的分割算法。首先，利用无人机正射影像数据开发了一种改进超像素算法和能量函数用于像素块分割和合并，实现分割二维树冠，其中，所提出基于能量函数的图像块合并的算法解决了现有普遍存在的过分割问题，且有更好的分割效果。其次，将二维掩码素点坐标与无人机点云数据坐标对齐完成三维的降维分割。然后，通过核密度函数估计确定每个初始分割段内聚类分量的个数，并采用高斯混合模型分离方法对每棵树的冠层进行分离。最后，通过实验对比验证了本研究分割方法的准确性和泛化性，本算法在不仅对稀疏林地有较高的分割准确率，在密集的林分中仍有较高的分割结果，分割精度高。综上所述，本发明方法通过使用改进超像素算法和能量函数提高了树冠边缘的分割精度，并通过使用核密度函数估计和高斯混合模型细化了点云分割以提高精度，通过对每个掩膜的检测推动了冠层分割朝着更好的解决方向发展。
附图说明
图1为研究地点的位置图。
图2为基于改进超像素算法和能量函数对正射影像进行初始分割流程图。
图3为基于改进超像素算法(加入边缘概率函数)的分割图。
图4为现有超像素分割图。
图5为归一化函数图。
图6为图像块(区域块)合并示意图。
图6中的(a)为两个不同簇不存在共同区域块的示意图。
图6中的(b)为两个不同簇存在共同区域块的示意图。
图7为点云精细分割图。
图7中的(a)为法向单冠点云示意图。
图7中的(b)为异常单冠点云示意图。
图7中的(c)为树顶检测点云示意图。
图7中的(d)为树顶检测直方示意图。
图7中的(e)为单个树冠点云分割结果图。
图8为点密度分布直方图和高斯核密度估计曲线图。
图8中的(a)为统计区间为0.3对应的点密度分布直方图。
图8中的(b)为统计区间为0.5对应的点密度分布直方图。
图8中的(c)为统计区间为0.7对应的点密度分布直方图。
图8中的(d)为带宽为0.3对应的高斯核密度估计曲线图。
图8中的(e)为带宽为0.5对应的高斯核密度估计曲线图。
图8中的(f)为带宽为0.7对应的高斯核密度估计曲线图。
图9为实验结果示意图。
图9中的(a)为低密度样地采用归一化无人机点云方法进行树冠分割结果示意图。
图9中的(b)为低密度样地采用局部极大值算法进行树冠分割结果示意图。
图9中的(c)为低密度样地采用点云分割算法进行树冠分割结果示意图。
图9中的(d)为低密度样地采用种子点堆叠算法进行树冠分割结果示意图。
图9中的(e)为低密度样地采用本文本研究的算法进行树冠分割结果示意图。图9中的(f)为中密度样地采用归一化无人机点云方法进行树冠分割结果示意图。
图9中的(g)为中密度样地采用局部极大值算法进行树冠分割结果示意图。
图9中的(h)为中密度样地采用点云分割算法进行树冠分割结果示意图。
图9中的(i)为中密度样地采用种子点堆叠算法进行树冠分割结果示意图。
图9中的(j)为中密度样地采用本文本研究的算法进行树冠分割结果示意图。图9中的(k)为高密度样地采用归一化无人机点云方法进行树冠分割结果示意图。
图9中的(l)为高密度样地采用局部极大值算法进行树冠分割结果示意图。
图9中的(m)为高密度样地采用点云分割算法进行树冠分割结果示意图。。
图9中的(n)为高密度样地采用种子点堆叠算法进行树冠分割结果示意图。
图9中的(o)为高密度样地采用本文本研究的算法进行树冠分割结果示意图。
图10为四种不同方法得到的实测冠宽与提取冠宽散点图。
图10中的(a)为低密度样地采用局部最大值算法得到的实测冠宽与提取冠宽散点图。
图10中的(b)为低密度样地采用点云分割算法得到的实测冠宽与提取冠宽散点图。
图10中的(c)为低密度样地采用种子点堆叠算法得到的实测冠宽与提取冠宽散点图。
图10中的(d)为低密度样地采用本文本研究的算法得到的实测冠宽与提取冠宽散点图。
图10中的(e)为中密度样地采用局部最大值算法得到的实测冠宽与提取冠宽散点图。
图10中的(f)为中密度样地采用点云分割算法得到的实测冠宽与提取冠宽散点图。
图10中的(g)为中密度样地采用种子点堆叠算法得到的实测冠宽与提取冠宽散点图。
图10中的(h)为中密度样地采用本文本研究的算法得到的实测冠宽与提取冠宽散点图。
图10中的(i)为高密度样地采用局部最大值算法得到的实测冠宽与提取冠宽散点图。
图10中的(j)为高密度样地采用点云分割算法得到的实测冠宽与提取冠宽散点图。
图10中的(k)为高密度样地采用种子点堆叠算法得到的实测冠宽与提取冠宽散点图。
图10中的(l)为高密度样地采用本文本研究的算法得到的实测冠宽与提取冠宽散点图。
具体实施方式
下面根据附图对本发明的具体实施方式作出进一步说明：
针对密集林分条件下的ITCs分割困难、总体分割精度较低等问题，本文结合无人机正射影像和点云数据的特点开发了一种高精度的分割算法。
(1)首先，利用无人机正射影像数据开发了一种改进超像素算法和能量函数用于像素块分割和合并而分割二维树冠。
(2)其次，将二维掩码素点坐标与无人机点云数据坐标对齐完成三维的降维分割。
(3)然后，通过核密度函数估计确定每个初始分割段内聚类分量的个数，并采用高斯混合模型分离方法对每棵树的冠层进行分离。
(4)最后，通过实验对比验证了本研究分割方法的准确性和泛化性，本算法在不仅对稀疏林地有较高的分割准确率，在密集的林分中仍有较高的分割结果。
1、材料与方法：
1.1、研究区域：
本研究使用的数据来自洪泽湖(31°32″N，118°89″E)，数据采集于2023年8月。该区域属于亚热带季风气候区，雨量充沛，光照充足，年平均气温15.5℃，年平均降水量1019.5mm，适宜植物生长。该区域以高冠层密度和树种组成复杂的人工林为主，有两个非常不同的植被层：树层，主要以柳树为主(平均高度6.72m)，草层，有各种各样的一年生物。主要经过几十年的森林管理、伐木和再造林，原来种植的树木原分布规律、林分年龄和垂直结构相对均匀，已演变成多层、复杂的森林结构。不同林分结构给ITCs的分割带来了巨大的挑战。为确保实验的有效性，本研究将种植时间较短、林分结构较单一的数据作为对比，在两片区域内选取了三种密度的柳树作为试验对象。选择区域如图1所示。
1.1.1、遥感数据与现场数据：
本研究选用大疆m600Pro采集影像和点云数据，其配备了三余度A3 Pro飞控、Lightbridge2高清数字图传、智能飞行电池组和电池管理系统，提供厘米级的精准定位。数据均来源于型号为gAirHawk GS-100C采集系统，该系统集成了Livox激光扫描仪、航空定位定向(POS)系统和可见光相机等，设备重量小于1050g。遥感影像数据通过圳市大疆创新科技有限公司生产的M3M无人机系统内置的相机SONY A6000获取，相机的主要参数为：2000万有效像素，视场角83°；最高分辨率5280*3956；图像大小:每图6-24MB；焦距:25mm；沿S型轨迹采集以便后续配准矫正。无人机点云数据同时段通过GS-100C激光传感器系统获取，该传感器测距精度为±20mm，激光脉冲频率为240kHz。波长和光束散度分别约为1550nm和0.5mrad。三次回波72万点/秒的频率确保良好的植被穿透能力。导航定位定向系统(POS)是通过全球定位导航卫星系统(GNSS)获取位置数据作为初始值，配合惯导系统((InertialMeasurement Unit,IMU)获取姿态变化的增量，通过卡尔曼滤波器以及反馈误差控制迭代运算，实时生成导航数据，定位精度在0.02-0.03m可以获取摄影相机的外方位元素和飞机的绝对位置，实现定点摄影成像和无地面控制的高精度对地直接定位，便于后续对数据的进一步处理。研究区域地面实测数据是由本团队共同完成，调查时间为2023年8月，测量的因子包括：树种、胸径、树高、东西和南北方向的冠幅。
无人机数据获取通过将照片和POS数据导入Pix4Dmapper软件中，通过影像配准、空间计算、平差处理和几何校正，最后进行影像镶嵌与正射校正，就完成了正射影像的制作，为保证高分辨率数字正射影像的结果的有效性，本研究设计的无人机飞行轨迹航向重叠度为60％，旁向重叠度为20％。
机载激光扫描系统在完成户外扫描工程作业后，获取了机载与基站的GNSS数据和IMU惯导数据，获取每个时刻所对应的空间位置和姿态对应的POS数据，然后将POS数据与原始波形数据进行融合解算，生成点云las文件，坐标系统一为WGS84。处理过程通过点云计算机软件gAirHawk完成。
1.2、方法论框架：
本研究的总体工作流程包括三个主要步骤：首先，三种密度的无人机影像数据、激光雷达数据和现场数据的采集和预处理。其次，通过改进的超像素分割方法将图片分割成小像素块，并使用能量函数对像素块合并，接着将分割掩膜按照像素坐标与同时空采集的三维点云数据坐标对齐完成三维分割。然后通过核密度估计算法检测树冠位置，并结合高斯混合模型完成三维树冠细分割。最后对三维分割结果检测，通过三维点云的空间优势分割出林下树。将该结果与LM、PCS和LS算法比较验证了该模型的有效性，以及在不同密度下比较验证了本研究模型的泛化能力。
1.2.1、数据处理：
图像预处理：
光学遥感影像中最常用分水岭算法对针叶林或稀疏林地单木分割，是一种基于地形的图像分割算法，能够识别图像中的区域和边界，该算法避免了噪声影响，一般先采用局部最大值滤波法提取树冠区域的特征点，或者采用高斯滤波消除树冠内部纹理和噪声，然后通过多尺度分水岭算法实现对混合落叶林的分割。以上方法对郁闭度较低的林地有较好的分割效果，然而局部最大值滤波常局限于滤波窗口的大小，窗口过大会导致小树遗漏，窗口过小会导致单个树冠被分为多个树冠，应用此方法往往需要对原始影像进行平滑处理消除冠层内部结构和噪声的影响，然而这种处理方式同直接高斯滤波一样会模糊树冠边界，易造成过分割和欠分割误差。目前对树轮廓清晰程度高的人工林树冠分割的研究尚少，密集林冠层具有冠层纹理复杂、排布不规则、局部制高点不突出且冠幅相似等特点，加之分割高密林存在地面与树冠轮廓混淆，树冠重叠严重，难以清晰定义个体。
为解决以上问题，本文参考Voronoi图邻近特性划分空间的原理，将像素点根据颜色阈值限制划分区域进行滤波。首先将图像转换到Lab空间，生成每个素点坐标，然后定义区域内颜色范围，遍历像素并根据颜色范围约束确定所属区域，其具体方法如下：a)根据图像特性设置颜色阈值的上下限，b)检查图像的像素点点颜色阈值是否在设置阈值范围内。同时满足以上两个条件的将划分为同一区。该方法能够快速将草地与林地区分开,同时也能使边界特征更明确。
ALS预处理：
在对点云的细分割之前，需要将点云归一化，因此需要先生成冠层高程模型(CHM)，CHM是一个二维表面模型，表示上冠层表面超过地面的高度。它消除了地形起伏对冠层高度的影响，从而反映了森林冠层高度的波动变化，在森林参数或森林生物量的反演中也起着重要作用。本文借助Arcgis软件，对导入的点云进行双边滤波去噪并采用平面栅格分割地面，然后生成数字高程模型(DEM)和数字表面模型(DSM)，从相应的地面高程值中减去DSM中每个点的高程值生成CHM模型，最终得到归一化的点云。
1.2.2、基于边界能量函数的区域中心点提取：
本文提出一种新的无人机影像分割算法。无人机影像的每个子区域内，各个样本都具有某些共同属性,如颜色、亮度和纹理等特征值。根据采用的方法不同，图像分割算法主要有四类：(a)基于聚类的算法；(b)基于区域合并的算法；(c)基于图论的算法；(d)基于分类的算法。本文采用了一种基于区域的方法。该算法首先进行超像素预分割，然后引入能量函数进行块合并，可以有效解决传统算法对边缘特征提取模糊的问题。其边缘能量函数表示为式：
EE＝BEnergy{jD(Edge)} (1)；
其中，jD(Edge)表示改进超像素算法分割后区域块j，然后基于能量函数将多个区域块合并成最终结果。
无人正射影像图初始分割流程图如图2所示。
图像块预分割：
提取Lab颜色空间和空间位置特征构建5维特征向量度量相似性，完成基于像素局部聚类的分割边界模糊问题。本文利用图像的颜色、坐标及纹理信息构造7维向量。
首先假设输入原图像I有N个像素，预期分割为K个像素块，块尺寸为N/K，每个聚类中心之间的距离定义为在初始化聚类过程中使用尺寸为3×3像素的窗口范围内计算像素点的梯度值，在该窗口范围内找到梯度值最小的像素点，并把聚类中心置于该位置。
然后将每个超像素看成图像处理的基本单元，计算两两超像素之间的欧氏距离，dc(x,k)表示像素点k与聚类中心点x在Lab色彩空间上的欧氏距离，两者距离越近则C(x,k)值越小；dl(x,k)表示像素点k到聚类中心点x的空间距离。
纹理的相似度为其纹元直方图的x2(cx,ck)距离，定义纹理的相似度，也即纹元直方图距离为：
其中cx(i)＝∑k∈w(x)I[T(k)＝i]，cx(i)表示在某个窗口w(x)内，某种纹元出现的频率(纹元直方图)，w(x)表示一个窗口，T(k)表示一个指针函数，cx、ck分别表示2个纹元直方图，式中x2(cx,ck)值越小，表示cx、ck越接近。
定义dw(x,k)为纹理的相似性度量：
dw(x,k)＝x2(cx,ck) (3)；
cx表示聚类中心点的纹元直方图，ck表示以聚类中心点为中心的一定范围内像素点k的纹元直方图。
本文将纹理特征分量融入到总度量计算中，使用边界概率函数g(x)，表示像素点x处于图像边界上的概率。所构建的总目标函数公式如：
其中，D是衡量像素点与聚类中心之间点的相似度(也即距离)，dw(x,k)为像素点k与聚类中心点x之间的纹理相似度(即纹理距离)；dc(x,k)表示像素点k与聚类中心点x在Lab色彩空间上的颜色向量欧氏距离；dl(x,k)表示像素点k到聚类中心点x的空间向量欧式距离，其中dc(x,k)与dl(x,k)采用现有公式计算，dw(x,k)采用公式(2)和公式(3)计算；λ1为常数量。g(x)为边界概率函数，表示像素点x处于图像边界上的概率。该项针对不同的像素点为颜色距离分量C(x,k)、空间距离分量L(x,k)和纹理相似度分量W(x,k)分配不同的权重，权重由g(x)确定。当像素点处于边界上的概率较大时，对这类像素在聚类划分过程就应更多地考虑像素和种子点的颜色距离，这样才能保证分割产生的超像素更能贴合图像边缘。所以此时应增加纹理分量W(x,k)权重，减小空间距离分量L(x,k)权重和颜色相似度分量C(x,k)权重。
其中，边缘概率g(x)是由图像梯度幅值f(x)经边缘概率函数映射得到，像素的梯度幅值在一定程度上可以反映像素落在图像边界上的概率。一般情况下在图像平坦区域梯度较小，图像边缘处梯度较大。虽然梯度幅值可以表征像素落在图像边界上的概率，但不同图像之间梯度信息千差万别，梯度幅值的值域范围较大，导致很难直接利用梯度幅值表征像素点处于边缘上的概率。因此本研究需要归一化函数对梯度幅值进行合理映射，得到像素点的边缘概率，使用如图3所示的归一化函数将梯度幅值映射到(0，1)范围内，即 该函数关于(m，0.5)中心对称；当梯度幅值在(0，t)范围内时，说明该像素点的梯度幅值过小，可认为该像素不是图像边缘处像素，此时边界概率映射函数将该梯度幅值映射为接近0的值；当像素点梯度幅值在(T，+∞)范围内时，说明此时像素点的梯度幅值足够大，可以被认为是图像边缘处像素，边界概率函数会将其映射为接近1的值。其中，a和m为超参数，本研究可通过调整a和b的值来调整两个阈值t和T的值。
利用式子(4)；计算聚类中心点与其2S×2S窗口范围内像素点的相似度，并将最相似聚类中心的标签赋予最小D值所对应的像素点(将每个像素划分到相似度最高的区域)。更新迭代聚类中心直到收敛，为了增加连贯性，将面积较小的超像素块划入相邻像素块中。
本文的改进超像素分割的主要步骤如下：
步骤1：初始化聚类中心。输入图像I，设置预期分割的超像素数目N，每个超像素块尺寸为N/K，聚类中心距离为
步骤2：调整聚类中心。为了防止预期分割的聚类中心处于图像边缘的情况，采取在初始化聚类中心点3×3像素的窗口范围内计算像素点的梯度值，在该窗口范围内找到梯度值最小的像素点，并把聚类中心置于该位置。
步骤3：计算相似度。利用式子(4)计算聚类中心点与其2S×2S窗口范围内每一个像素点的相似度(也即距离)，并将最相似聚类中心的标签赋予最小D值所对应的像素点。
步骤4：更新迭代聚类中心。重复步骤(3)，直到收敛。
步骤5：增加连通性。对于面积较小的超像素块，采取就近划入与其相邻的超像素块的原则。
其分割结果如图4所示。图4中a处的边界点分割效果明显优于图5的a′，体现了加入边缘概率函数后的有效性，图4中b处的像素分割效果明显优于图5的b′，当像素位于边界点时，根据纹理分量的差距将树冠重叠边缘划分开。同时也保证了超像素边界与图像边界的贴合。
基于能量函数的图像块合并：
本研究发现树冠数据呈现两个特征：1)虽然树冠本身比较复杂，但是局部信息比较丰富，即树木冠层比较容易识别，2)局部信息丰富的区域之间距离较远。3)树木冠层呈拟圆形，本研究可以对合并形状进行定义。因此为了更加清楚的识别局部值，本研究对图像进行增强并转换成归一化的能量图。
已知在预分割中将图像分割成K块形状各异的区域块，假设区域块j的聚类中心为Pj，此图像可以表示为：P＝(p1,p2,…,pj,…,pM)，区域j中的点的集合为 其中表示在n在区域块j上的第n个点。则该区域的代表点为Pj，本研究定义各个点的能量值为则区域块j的聚类中心为区域块j的区域代表点Pj，将区域块j的能量值定义为区域代表点Pj的局部平均能量值，区域块j的能量值为：
表示区域块j内像素点r″的能量值。也就是说区域块j的能量值为区域块j中所有像素点的能量值求和并平均。
基于能量函数的图像块合并的步骤包括：
步骤a：按照公式(5)计算每个区域块的能量值；其中某区域块的能量值即为某区域块中所有像素点的能量值求和并平均；
步骤b、步骤b1、选择聚类中心：
其中，是以区域块j为聚类中心的区域块j的能量值，其中以区域块j为聚类中心的区域代表点Pj记为聚类中心区域代表点Pj，是所有区域块的能量值的均值，是同一个簇中区域块的区域代表点距聚类中心区域代表点Pj的距离，是所有的期望，表示同一簇中区域块的区域代表点距聚类中心区域代表点Pj的最远距离；
满足公式(6)和公式(7)的区域块j即为聚类中心；
初始判定为同一簇的条件：
从所有区域块中，选取最高能量值的区域块的区域代表点Pj为中心，八个方位进行搜索，在八方位中搜索距离中心最远的方位为基准，当位于最远方位的区域块ni的区域代表点的局部平均能量值下降到一定值ε或当位于最远方位的区域块ni的区域代表点的局部平均能量值不再继续下降时，则以该区域块ni的区域代表点和最高能量值的区域块的区域代表点Pj的距离为半径，区域代表点Pj为圆心画圆，该圆范围内涉及到的区域块归属于同一个簇；
步骤b2、优化簇：
Bj表示以区域j为中心多个区域块与之为同一簇；
以区域块j为聚类中心，当通过步骤b1判断与区域块j归属于同一簇内的某区域块bi与区域块j的能量相似度值大于等于同一簇中位于最远方位的区域块ni与区域块j的能量相似度值时，值为1，表示该区域块bi归属于以区域块j为聚类中心的簇；
以区域块j为聚类中心，当通过步骤b1判断与区域块j归属于同一簇内的某区域块bi与区域块j的能量相似度值小于同一簇中位于最远方位的区域块ni与区域块j的能量相似度值时，值为0，表示该区域块bi不属于以区域块j为聚类中心的簇，从簇中删除该区域块bi，如图6中，某区域块上写叉号的，即为所删除的区域块；
其中区域块j与区域块bi的能量相似度值为： 表示区域块bi的能量值，表示区域块j的能量值；区域块ni与区域块j的能量相似度值为： 表示区域块ni的能量值，中分子表示区域块ni与区域块j的能量差度值，K1表示通过步骤b1判断与聚类中心区域块j同属于同一个簇内的所有区域块的集合,表示区域块f1的能量值，其中区域块f1属于K1；中分母的含义为将属于K1的所有区域块与区域块j的能量差度值进行求和。
步骤c、计算未归属于簇的剩余区域块的能量值；按照步骤b的方法在剩余区域块中再选取聚类中心，判定归属于同一个簇的区域块，并优化该簇；
步骤d、将未归属于簇的剩余区域块按照步骤c的方法再选取聚类中心、判定归属于同一个簇的区域块，并优化该簇，直到所有区域块均归属于对应的簇为止；
步骤e、若一个区域块共同归属于两个或者两个以上不同的簇，如图6(b)所示，两个圆有共同重叠的区域块，则计算该区域块分别与共同归属的不同簇内聚类中心区域块的能量相似度值并比较，该区域块归属于其中能量相似度值最高所对应的聚类中心区域块所归属的簇；
步骤f、合并归属于同一个簇的多个区域块，多个区域块合并后形成一个大的区域块，因此，步骤f合并后得到多个区域块，但是区域块的数量与基于改进超像素分割后得到的区域块的数量相比大大减少，解决了改进超像素算法过分割的情况。
上述过程中，是以高能量值区域代表点Pj为中心，八个方位进行搜索，并选择最远距离进行合并，在四方位中搜索最远的方位为基准，当某区域块的能量值下降到一定值ε或当某区域块的能量值不在继续下降时，则以该区域块的区域代表点和高能量之中心点Pj的距离为半径，预合并圆内的区域块，在合并范围内有可能存在其他方位能量值不符合问题，所以按照上述条件再筛选，如图6中的(a)所示；同时也会存在一个区域同时在两个树冠预合并范围内的情况，则根据该区域块分别与共同归属的不同簇(不同圆)内聚类中心区域块(为聚类中心的区域块)的能量相似度值进行划分，如图6中的(b)所示。图6中的(a)表示距离较远的合并方式，图6中的(b)表示存在两个树冠同时合并范围内的区域。
所提出基于能量函数的图像块合并的算法解决了现有普遍存在的过分割问题，且有更好的分割效果。
1.2.3、雷达点云分割：
通过Arcgis软件对二维影树冠轮掩膜提取，根据同时空拍摄下点云与影像坐标点相对应的特点，将掩膜的栅格数据与无人机雷达点云的CHM相对应对CHM分割，以分割后的矢量边界为参考拾取CHM中属于边界范围内的点，并随机赋予不同的颜色，这样以二维分割促使三维分割从而达到降维分割的效果。
核密度估计函数对树冠个数检测：
在实地调查中，每一棵树的位置在其茎的底部被正确地测量出来。由于机载激光雷达点云的点密度通常从上到下下降，这导致树干及以下结构细节的点数量相对于树的其他部分较少即林下信息缺乏，特别是在重叠的森林中，机载雷达点云难以准确定位单个树干的位置。与树干类似，树的顶点也被用作单棵树的标记，它代表了点云中单棵树的特定数量。如图7中(c)所示，每棵树中心点的密度通常都很大，从中心到两侧的点密度呈下降趋势，单个树冠在水平面上的投影轮廓(即正射影像分割轮廓)一般近似为圆形，而多个树点云的投影轮廓和空间分布近似为椭圆，所以按照长轴分割点簇能更好的检测出多棵树。因此，本文提出了长轴方向核密度估计来定位潜在的树的顶端，即找出沿树冠投影下椭圆的长轴，以长轴为x轴进行核密度估计来寻找树冠位置。为了准确地检测点密度的局部最大值，采用核密度估计计算各初始段的概率密度函数分布，核密度估计定义为：
其中，n′为每个初始分割中的点数，h>0为一个平滑参数，被称作带宽，K(x′)为核函数，x′表示样本点，x′γ表示观测点，本文选用高斯核函数进行密度估计。
高斯混合模型树冠分割：
在粗分割后，对每个掩膜下的聚类进行精细分割，K-means算法常用于三维点云分割，但对于非凸形状数据分割具有很大的局限性，因此本文采用适合复杂数和非凸形状数据的GMM对树冠点云进行建模，该方法通过判别数据点属性该列将其划分到不同类别完成三维点云的精细分割。在投影到三维点云中时，本研究为每个掩膜下的点簇构建了局部坐标参考框架，构建了当前点聚类的局部坐标参考框架后，以逆时针方向每间隔60°将点簇从0°到180°分割成三组进行进一步分析，以x轴所在方向即0°为起始，以xoy投影面沿x轴最远点的连线记为d1，60°方向时最远点距离记为d2，120°方向时最远点距离记为d3，这意味着本研究在实现中将对0°、60°和120°方向上的每个点簇距离进行分析。本研究选择三个方向距离最大值的剖面投影水平方向作为x轴，即max{d1,d2,d3}，图7中(b)是沿逆时针方向旋转0°时多向三维空间结构图，沿x轴方向进行核密度估计。
然后，基于核密度估计的树冠顶点数(简称树冠点数)并结合高斯混合模型进行点云分割，其公式可表示为：
其中，N′表示上述核密度估计的树冠点数，X表示投影在剖面上的点(xp,zp)，λk′为比例系数，表示各混合分量的先验概率，μk′、δk′表示高斯分布初始参数，分别表示高斯函数的平均值和方差。该算法通过E-step计算分量概率，通过M-step更新高斯混合参数μk′、δk′，不断迭代计算混合分布参数，当参数值小于阈值或迭代次数达到最大EM步骤收敛，从而将点云分类，图7中(e)为优化为两个树冠的结果。
图7为点云精细分割图，其中图7中的(a)表示法向单冠，图7中的(b)表示异常单冠，即包括两个树冠，图7中的(c)和(d)表示树顶检测，(e)表示单个木材分割的结果。
1.2.4、验证指标：
为了避免地形复杂、冠层密度高、树冠不规则生长导致树冠人工测量误差，本研究以数字正射影像数据作为基础图，手工将树冠勾画为真实树冠边界，以进行精度评价。将ITCs分割得到的单个树的位置与真冠重叠率，以确定真阳性(TP)、假阳性(FP)和假阴性(FN)的数量。当参考树冠与分割树冠重叠率≥70％时，认为分割正确，当参考树冠与分割树冠重叠率≤70％时，认为分割错误。本研究计算了3个样本方块的ITS结果，得到了TP、FN和FP。本研究还计算了F_score(f,overall accuracy)分数，以评估考虑错误和遗漏分数的总体准确性。
式中，r表示单个树的检出率；p表示单个树检测的正确率；f为ITCs分割的总体精度。
式中，g为正确分割的单个树的数量；Xτ表示正确分割单个树的冠宽；表示单个树的平均宽度，xτ表示单个分割树的测量冠宽，表示单个分割树的测量平均冠宽，将从本研究的方法中分割提取的冠宽与测量的冠宽进行了比较。采用相对系数R2、均方根误差(RMSE)来评估提取的冠的宽度。
2、结果：
2.1、核密度估计参数选择：
聚在一起的不同树的核密度分布曲线根据核带宽的大小而发生变化，因此可以通过交叉验证或根据领域知识来选择合适的核带宽参数，本文通过对13组树冠组进行试验，通过分类准确度来选择最佳的核带宽参数。
表1：
本文以图7中的两个树冠分割为例，图8表示统计区间分别为0.3、0.5、0.7的密度直方图，以及该带宽情况下的核密度估计曲线，可以看出在带宽为0.5时能正确找到树冠点，即最佳带宽为0.5。
图8表示点密度分布直方图和高斯核密度估计曲线，图8中(a)中的统计区间为0.3，图8中(b)中的统计区间为0.5，图8中(c)中的统计区间为0.7，图8中(d)中的带宽为0.3，图8中(e)中的带宽为0.5，图8中(d)中的带宽为0.7。
2.2、不同点云密度下itcs分割的比较：
将三种不同密度的样本划分，分别为低密度样本、中密度样本和高密度样本。在图9中，用不同四种分割方法分割三种密度的树。局部极大值算法(LM算法)和种子点堆叠算法(LS算法)也能准确地识别树顶的位置。本实验则是通过准确识别边界并根据其垂直特性进行分割，并使用不同颜色可视化结果。图9中低密度样地为稀疏且种植均匀的柳树林，该片林地中有部分树为幼树，单株树冠面积较小；中密度样地为密集但种植均匀的柳树林，该地的柳树较为成熟；高密度样地为密集且种植不均匀的柳树林，其中部分树的周围衍生出了新的小柳树。在四种单独的树冠分割方法中，局部极大值(LM算法)和种子点堆叠方法(LS算法)对树的位置得到了满意的检测结果，其中在有大冠层的林地中，种子点堆叠
(LS算法)检测到多个有小冠层的树，但对一些相邻树的分割不足，用点云分割方法(PCS算法)能检测树冠且能很好地单个树冠的实际大小，但该方法分割单棵树的冠层大小相对均匀，意味着该方法在实际分布不均匀的林地不能很好的分割，即在高密度样地中分割效果较差。而本研究的方法既能有较好的检测率同时能很好地描绘出树冠的边界。
图9为实验结果示意图。图9中(a)-(e)为低密度样地，(f)-(j)为中密度样地，(k)-(o)为高密度样地，(a)、(f)、(k)为归一化无人机点云，(b)、(g)、(l)为局部最大值算法，(c)、(h)、(m)为点云分割算法，(d)、(i)、(n)为种子点堆叠算法，(d)、(j)、(n)为本研究的算法。
通过将3个样地中四种分割方法获得的单个树的分割结果与样地中每棵树的测量位置数据进行匹配，评估其准确性，如表2所示。不同方法对ITCs的分割精度存在明显差异，3种不同密树的r、p和f的浮动范围分别为0.81～0.98、0.71～0.95和0.76～0.97。实验表明本研究的方法最好，本研究的算法分割的树冠相对完整，边界清晰，破碎的树冠较少。PCS算法的分割结果有许多欠分割，LM算法能检测到小的树冠，但是分割准确率不高，在高密度数据中，本研究的与LS算法相比，本研究的r、p和f分别提高了0.02、0.04、0.03。不同密度的样地对ITCs的分割精度也有较大的差异，在低密度区域，四种方法的r、p、f的差异不大，但PCS和本研究的算法在分割精度上略高于LM和LS，在高密度区域，由于存在部分衍生树苗紧挨着大树，所以较难分割，尤其是PCS方法，在密集林中比表现较差。由于本研究的算法在二维分割时对边缘进行细化，所以对于密集树冠来说，树冠分割精度表现最佳。
表2：
由于本算法涉及衍生树苗的优化分割，精度反映了树的分割准确率，因此很明显，本研究的分割的精度提高最大，其次是同时考虑r和精度的f分数。召回提高最小，因为它与过分割样本的数量有关。
2.3、树木参数的准确性：
根据地块类型，建立了提取的树高与实测树高、提取树冠与实测树冠之间的线性关系。计算出了提取的树高与实测树高、提取树冠与实测树冠之间的R2和RMSE，计算结果见表3；
表3：
LM：局部最大值算法，PCS：点云分割算法；LS：种子点堆叠算法；OURS：本研究的算法。
树高参数的准确性：
实验结果表明，不同的分割方法在三种密度的地块之间的分割精度上存在显著差异，但分割结果的精度与实际匹配后计算出的树高提取精度之间的差异很小。R2的最大值和最小值之间的差值在0.06以内，RMSE的差值小于0.12m。这些结果表明，四种方法的树高提取结果是稳定的，四种方法可以提取涉及无人机-激光雷达点云的个别树高参数。
冠宽参数的准确性：
图10显示了不同方法下三个样本的柳树冠宽估计与测量值的线性回归，根据R2和RMSE验证分割质量，虽然由于实际限制，不可能进行真实位置测量，但实验样地的冠层大多是带有明显的顶部，但仍然可以清楚地看到下冠层林分，通过可视化和人工测量的点云获得的单个树的信息，虽然精度有限，但仍可作为ITCs结果的参考值。三种密度林地冠宽的线性回归结果如图10所示。从R2值来看，由分割算法计算出的估计值之间有很强的相关性。
图10为四种不同方法得到的实测冠宽与提取冠宽散点图。图10中的(a)、(b)、(c)、(d)为低密度样地，(e)、(f)、(g)、(h)为中密度样地，(i)、(j)、(k)、(l)为高密度样地，(a)、(e)、(i)为局部最大值算法，(b)、(f)、(j)为点云分割算法，(c)、(g)、(k)为种子点堆叠算法，(d)、(h)、(l)为本研究的算法。
低密度样地中树冠分布较均匀，且树龄相对较小长势较一致，易于分割和测量，四种方法分割下树冠的R2和RMSE相差不大，在中密度样地中，树木分布较均匀，但由于该样地种植较久且密集，出现了部分树生长空间被周围树占领，导致树冠大小不均衡，而高密度样地中，在不均匀的密集种植条件下，样地中衍生出小树，使得树冠更加紧凑且导致树冠更加分割困难，LM、PCS和LM算法的ITCs分割精度均受到环境的不同程度的影响，其中PCS表现最差，本研究的算法表现最好，R2为0.93，RMSE为0.42m，这主要是由于考虑到柳树冠层属性，一定程度上降低了衍生树对ITCs分割的影响。
3、结论：
近几年无人机激光点云的单木分割方法层出不穷，但高郁闭度的单木分割仍然具有挑战性，本研究结合无人机正射影像和点云数据的特点提出了一种高精度的ITCs分割算法。对树冠特征的分析，该方法通过使用改进超像素算法和能量函数提高了树冠边缘的分割精度，通过使用局部最大值法寻找剖面种子点细化了点云分割以提高精度，通过对每个掩膜的检测推动了冠层分割朝着更好的解决方向发展。根据人工测量对3种密度类型的林地进行验证，冠层的检出率R≥0.81，同时树冠宽度估计的决定系数R2≥0.93，均值误差RMSE≤0.42m，高度估计的决定系数R2≥0.86，均值误差RMSE≤0.36m。在低密度样地中，四种方法表现的效果差异较小，在高郁闭度的高密度样地中，LM和LS算法的检测出更多的树冠，但树冠分割准确率不佳，PCS方法适用于树冠大小均匀的样地，因此在高郁闭度样地中分割效果最差，本研究的方法检出率高且树冠边界的分割上具有优越性，而且二维分割投影到三维一定程度上降低了点云分割的时间和复杂度。
本发明的保护范围包括但不限于以上实施方式，本发明的保护范围以权利要求书为准，任何对本技术做出的本领域的技术人员容易想到的替换、变形、改进均落入本发明的保护范围。
Claims (4)
1.一种基于无人机影像与激光点云的多层次密集树冠分割方法，其特征在于：包括：
步骤1：采集林分点云数据和林分正射影像数据；
步骤2：通过Voronoi划分对林分正射影像进行滤波，通过双边滤波对林分点云数据进行滤波；
步骤3：对滤波后的点云归一化处理，并用平面栅格法去除地面点云；
步骤4：使用改进超像素算法和能量函数对滤波后的正射影像进行分割；
步骤5：提取分割后的掩膜坐标，将二维图像的掩膜坐标与归一化后的点云坐标对齐进而完成三维分割，得到多个初始分割区域；
步骤6：采用基于长轴的核密度估计算法确定每个初始分割区域内树冠顶点的个数；
步骤7：基于树冠顶点的个数并采用高斯混合模型进行点云分割，获得树冠分离结果。
2.根据权利要求1所述的基于无人机影像与激光点云的多层次密集树冠分割方法，其特征在于：所述的步骤4具体为：
步骤4.1、从滤波后的正射影像中提取颜色距离分量、空间距离分量以及纹理相似度分量；
步骤4.2、基于颜色距离分量、空间距离分量以及纹理相似度分量并采用改进超像素算法将正射影像分割成多个区域块；
步骤4.3、基于能量函数合并区域块；
步骤4.3.1、计算每个区域块的能量值，其中，第j个区域块的区域代表点Pj的局部平均能量值简称为区域块j的能量值，区域块j的能量值为：
其中，表示区域块j内像素点r′′的能量值；
步骤4.3.2、选择聚类中心：
其中，是以区域块j为聚类中心的区域块j的能量值，其中以区域块j为聚类中心的区域代表点Pj记为聚类中心区域代表点Pj，是所有区域块的能量值的均值，是同一个簇中区域块的区域代表点距聚类中心区域代表点Pj的距离，是所有的期望，表示同一簇中区域块的区域代表点距聚类中心区域代表点Pj的最远距离；
满足公式(2)和公式(3)的区域块j为聚类中心；
初始判定为同一簇的条件：
以最高能量值的区域块的区域代表点Pj为中心，八个方位进行搜索，在八方位中搜索距离中心最远的方位为基准，当位于最远方位的区域块ni的区域代表点的局部平均能量值下降到一定值ε或当位于最远方位的区域块ni的区域代表点的局部平均能量值不再继续下降时，则以该区域块ni的区域代表点和最高能量值的区域块的区域代表点Pj的距离为半径，区域代表点Pj为圆心画圆，该圆范围内涉及到的区域块归属于同一个簇；
优化簇：
以区域块j为聚类中心，当区域块bi与区域块j的能量相似度值大于等于同一簇中位于最远方位的区域块ni与区域块j的能量相似度值时，区域块bi归属于以区域块j为聚类中心的簇；
以区域块j为聚类中心，当区域块bi与区域块j的能量相似度值小于同一簇中位于最远方位的区域块ni与区域块j的能量相似度值时，区域块bi不属于以区域块j为聚类中心的簇；
其中区域块j与区域块bi的能量相似度值为： 表示区域块bi的能量值；区域块ni与区域块j的能量相似度值为： 表示区域块ni的能量值，K1为与聚类中心区域块j同属于同一个簇内的所有区域块的集合；
步骤4.3.3、计算未归属于簇的剩余区域块的能量值，按照步骤4.3.2的方法在剩余区域块中再选取聚类中心，判定归属于同一个簇的区域块，并优化该簇；
步骤4.3.4、将未归属于簇的剩余区域块按照步骤4.3.3的方法再选取聚类中心、判定归属于同一个簇的区域块，并优化该簇，直到所有区域块均归属于对应的簇为止；
步骤4.3.5、若一个区域块共同归属于两个或者两个以上不同的簇，则计算该区域块分别与共同归属的不同簇内聚类中心区域块的能量相似度值并比较，该区域块归属于其中能量相似度值最高所对应的聚类中心区域块所归属的簇；
步骤4.3.6、合并归属于同一个簇的区域块。
3.根据权利要求1所述的基于无人机影像与激光点云的多层次密集树冠分割方法，其特征在于：
所述步骤6中，核密度估计算法的公式为：
公式(4)中，n′每个初始分割中的点数，h>0为一个平滑参数，被称作带宽，K(x′)为核函数，x′表示样本点，x′γ表示观测点。
4.根据权利要求1所述的基于无人机影像与激光点云的多层次密集树冠分割方法，其特征在于：所述的步骤7中，采用高斯混合模型进行点云分割，其公式为：
公式(6)中，N′表示核密度估计的树冠点数，X表示投影在剖面上的点，λk′为比例系数，表示各混合分量的先验概率，μk′、δk′表示高斯分布初始参数，分别表示高斯函数的平均值和方差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311791699.3A CN117765006A (zh) | 2023-12-25 | 2023-12-25 | 基于无人机影像与激光点云的多层次密集树冠分割方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311791699.3A CN117765006A (zh) | 2023-12-25 | 2023-12-25 | 基于无人机影像与激光点云的多层次密集树冠分割方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117765006A true CN117765006A (zh) | 2024-03-26 |
Family
ID=90325217
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311791699.3A Pending CN117765006A (zh) | 2023-12-25 | 2023-12-25 | 基于无人机影像与激光点云的多层次密集树冠分割方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117765006A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117970330A (zh) * | 2024-03-29 | 2024-05-03 | 月明星(北京)科技有限公司 | 基于聚合像素的时序干涉sar滤波方法和装置 |
-
2023
- 2023-12-25 CN CN202311791699.3A patent/CN117765006A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117970330A (zh) * | 2024-03-29 | 2024-05-03 | 月明星(北京)科技有限公司 | 基于聚合像素的时序干涉sar滤波方法和装置 |
| CN117970330B (zh) * | 2024-03-29 | 2024-06-11 | 月明星(北京)科技有限公司 | 基于聚合像素的时序干涉sar滤波方法和装置 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yang et al. | An individual tree segmentation method based on watershed algorithm and three-dimensional spatial distribution analysis from airborne LiDAR point clouds | |
| CN110570428B (zh) | 一种从大规模影像密集匹配点云分割建筑物屋顶面片的方法及系统 | |
| CN108109139B (zh) | 基于灰度体元模型的机载lidar三维建筑物检测方法 | |
| Shendryk et al. | Bottom-up delineation of individual trees from full-waveform airborne laser scans in a structurally complex eucalypt forest | |
| Williams et al. | 3D segmentation of trees through a flexible multiclass graph cut algorithm | |
| CN109146889A (zh) | 一种基于高分辨率遥感图像的农田边界提取方法 | |
| Harikumar et al. | A local projection-based approach to individual tree detection and 3-D crown delineation in multistoried coniferous forests using high-density airborne LiDAR data | |
| CN112819830A (zh) | 基于深度学习与机载激光点云的单株树冠分割方法 | |
| CN112907520B (zh) | 基于端到端深度学习方法的单株树冠检测方法 | |
| CN109961470B (zh) | 一种基于激光点云的活立木叶属性精准估测方法 | |
| CN108074232B (zh) | 一种基于体元分割的机载lidar建筑物检测方法 | |
| CN103218787A (zh) | 多源异构遥感影像控制点自动采集方法 | |
| CN112099046B (zh) | 基于多值体素模型的机载lidar三维平面检测方法 | |
| CN110047036B (zh) | 基于极坐标格网的地面激光扫描数据建筑物立面提取方法 | |
| CN109766824B (zh) | 基于模糊证据理论的主被动遥感数据融合分类方法 | |
| Hao et al. | A hierarchical region-merging algorithm for 3-D segmentation of individual trees using UAV-LiDAR point clouds | |
| CN113781431A (zh) | 基于城市点云数据的绿视率计算方法 | |
| CN117765006A (zh) | 基于无人机影像与激光点云的多层次密集树冠分割方法 | |
| Li et al. | Measuring detailed urban vegetation with multisource high-resolution remote sensing imagery for environmental design and planning | |
| CN116994029A (zh) | 一种用于多源数据的融合分类方法及系统 | |
| CN117036971A (zh) | 自适应局部空谱一致下的机载LiDAR数据建筑物提取方法 | |
| Tóvári | Segmentation based classification of airborne laser scanner data | |
| Huang et al. | Ground filtering algorithm for mobile LIDAR using order and neighborhood point information | |
| Berhe | Towards a UAV based standalone system for estimating and mapping aboveground biomass/carbon stock in Berkelah Tropical Rain Forest, Malaysia | |
| Liu et al. | A Comprehensive Comparison of Individual Tree Crown Delineation of Plantations using UAV-LiDAR Data: A Case Study for Larch (Larix olgensis) Forests in Northeast China |
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 |