CN107705309A - 激光点云中林木参数评估方法 - Google Patents

激光点云中林木参数评估方法 Download PDF

Info

Publication number
CN107705309A
CN107705309A CN201710955957.5A CN201710955957A CN107705309A CN 107705309 A CN107705309 A CN 107705309A CN 201710955957 A CN201710955957 A CN 201710955957A CN 107705309 A CN107705309 A CN 107705309A
Authority
CN
China
Prior art keywords
mrow
msub
msup
point
cloud
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.)
Granted
Application number
CN201710955957.5A
Other languages
English (en)
Other versions
CN107705309B (zh
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.)
Zhejiang Zhongnan Green Carbon Technology Co.,Ltd.
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 CN201710955957.5A priority Critical patent/CN107705309B/zh
Publication of CN107705309A publication Critical patent/CN107705309A/zh
Application granted granted Critical
Publication of CN107705309B publication Critical patent/CN107705309B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20112—Image segmentation details
    • G06T2207/20152—Watershed segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation
    • G06T2207/30188—Vegetation; Agriculture

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Image Processing (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

本发明提供一种测定数据方便，不受环境影响，评估结果准确的激光点云中林木参数评估方法。该方法是对点云数据进行窗口平滑操作，以消除地形高度差；对点云特征进行枝叶分离操作，并获取主枝干信息；计算主枝干点云的法矢量，求取瑞利商获取主枝干倾斜角；结合投影分水岭算法，实现对橡胶林段的单株提取；株株分离后进行反演冠幅与叶片分布，分别计算了冠幅、胸径、冠积、倾斜度、树叶分布与叶面积指数等参数。

Description

激光点云中林木参数评估方法
技术领域
本发明涉及对树林中树木的有关参数如棵数、树高、胸径、冠体积、冠幅等进行测定评估的方法。
背景技术
中国的亚热带森林，是以常绿阔叶林为主要乔木的生物群落，覆盖着全国大约25％的面积。森林对二氧化碳下降、动物群落、调节水文流动和巩固土壤起着重要的作用。此外，森林生物量是评估森林碳贮量的重要参数。若森林遭受破坏，则会产生长期和不可逆的影响。因此，如何对森林进行实时监测，准确而又高效地获取森林相关资源信息，研究森林的生长状况及其动态变化具有非常重要的意义。
正是由于森林的重要性，近年来，国内外许多的学者做了大量的面向森林监测的相关工作。按研究区域的尺度可分为以下三类：
一、在宏观尺度下，运用星载数据的大区域影像，针对林木参数反演开展了较多的工作。其中：基于QuickBird的高分辨率VHR星载图像对城区地区进行定位和树冠描绘；面向WorldView-2卫星数据设计半自动化方法实现对树冠自动定位；面向Ikonos卫星数据对人工针叶林实现树冠分割以及单株树的定位与提取。基于Landsat 8卫星数据针对赞比亚米隆博林进行树冠描绘，进而实现对树冠的定位；面向ALOS PALSAR和Landsat ETM+卫星图像数据智力中部森林中的樟子松和桉树进行树冠分割。虽然宏观尺度下获取到了树木的参数信息，但由于大尺度下精度低和高光谱变异性问题，以及树种、树龄、环境等因素的影响，都会造成基于宏观尺度下的测量方法对不同种类的树冠边缘识别与冠中心定位不准确等误差。
二、在中观尺度下，如机载Lidar、小尺度激光Lidar扫描仪、航空影像、移动式激光扫描和高光谱小型机载光谱成像仪(CASI)等。中观尺度下的研究工作，在树顶遮挡、参数扫描和树木差异性上都存在着问题，这对结果的获取和对单株树木提取都有着一定的影响。
三、在微观尺度下，目前的方法主要有：落叶法、角规测树法、半球摄影法，叶面积测定法和LAI-2200植物冠层分析仪等方法。微观尺度下所采用的方法都较为简单，不需要特殊的仪器，较为方便快捷，但多是人力测量，其耗力大，而且测量工作容易受外界因素影响，如季节、天气、光照和色差等，这些因素都会导致对林分参数的测量照成很大的困难，测定工作历经时间相对较长。
本发明的目的是提供一种测定数据方便，不受环境影响，评估结果准确的激光点云中林木参数评估方法。
本发明所述的激光点云中林木参数评估方法，包括以下步骤：
步骤A:获取点云数据：以激光雷达对树林自下而上扫描，获取树林的点云数据pi(xi，yi，zi)，xi、yi、zi分别表示点pi在x、y、z坐标轴方向上的坐标值；
步骤B:对树木枝叶分离：
对于点云数据中的一点pi＝(xi,yi,zi)T,其中对于点pi，半径为rm内的领域定义为：pj＝(xj,yj,zi)T，且满足条件：||pj-pi||≤rm；
点pi邻域的协方差矩阵定义为Cp：
当时，其中μ表示pi邻域所有点在空间中位置的均值；
重新定义一个新的坐标系，特征向量代表轴方向，相对应的特征值代表沿轴的点方差；设ek,i是矩阵Cp的特征向量，λk,i是相应的特征值，k＝0,1,2且λ0,i≤λ1,i≤λ2,i；λk,i定量的显示沿轴ek,i方向的数据方差；
用e0,i＝{ei,x,ei,y,ei,z}表示最小特征值λ0,i所对应的特征向量，e0,i也代表着点pi处的法向量；
将对应的每个点云pi的方差矩阵的特征值归一化，即：
计算结构张量特征，包括结构张量的平面性特征(c0,i)、熵特征(c1,i)与线性特征(c2,i)，分别计算为
使用以下等式计算每个点云pi的法向量分布：
其中e0,j是点云pi的法向量，是点云pi邻域基团的平均法向量分布，且
Vp三个特征值{l0,l1,l2}的相对幅度满足条件：l0≥l1≥l2，这反映出点云pi邻域的几何特征；
基于上述分析，一系列关于点云pi的特征被计算出来，记为:其中{eix,eiy,eiz}为点pi的法向量，{c0,i,c1,i,c2,i}为点pi的结构张量特征，而{l0,i,l1,i,l2,i}为点pi的形状特征值；
当特征值符合l0≈l1≈l2时，基团点云呈现球形状，对应为果实；当特征值分布为l0≥l1≈l2时，基团点云呈现线形状，对应为枝干；而当特征值分布为l0≈l1≥l2时，基团点云呈现面状，对应为叶片；也就是说，根据获得的点pi的特征fpi，对树木进行枝叶分离操作；
基于枝叶分离后的枝干数据，按照树木自然生长方向，以树总高度的1/2区域作为主枝干数据，其余部分为分支干数据；
步骤C:树冠中心定位，确定树冠的中心点坐标；
步骤D:倾斜度分析，计算出则主枝干倾斜度α、主枝干与分支干之间的夹角β；
步骤E:基于meanshift和分水岭算法进行株株分离，得到树木总棵数和每株树冠的点云，从这些点云中计算得出树高、胸径、冠体积、冠幅。
树高计算：
对于林分中的每株树木的点云数据的z轴数据求最大值，即作为林分树木的树高。胸径及平均胸径计算：
(1)针对枝叶分离后获得的树干数据，从中选取20棵树，计算每棵树木树干左端点坐标(x1,y1)和最右端点的坐标(x2,y2)，并计算出□x＝x2-x1和□y＝y2-y1。根据计算出的□x和□y，由公式计算出树干的胸径，R的值即为胸径。
如，任意一棵树的左右端点坐标为(17.02999964.399994)、(16.63999964.229987)，其□x＝-0.3900(m)，□y＝-0.1700(m)，
同理，计算出其他19棵树的胸径。选取的20棵树的□x、□y和R值为：
(2)R值一列中，去除最大值0.6419(m)和最小值0.3397(m)，其余18组数据求其平均值为0.4726，即把均值0.4726(m)作为选取的20棵树的平均胸径。
冠体积及平均冠体积计算：
(1)基于计算机图形学中的最大凸包算法(convexhull)，对每株树的叶子点云数据求取冠体积。
如，任选15棵树，基于这15棵树的叶子点云数据和最大凸包算法，求得的冠体积为：
(2)对这15组数据，从中去掉一组最大值564.5472(m3)和一组最小值45.8267(m3)后求平均数，求得平均树冠体积为141.7978(m3)，即为选取的15棵树的平均冠体积。冠幅及平均冠幅计算：
以东西方向上的冠幅为例：
(1)从株株分离后的林分树木的点云数据的俯视图中，计算每排树木最左端和最右端点的坐标。东西方向上选取7排树木，其中一排的最左端点坐标为(-23.47008.6700)，最右端点坐标为(39.1700102.8400)，则其□x＝62.6400(m)，□y＝94.1700(m)，则这7排数据的□x、□y和R为：
(2)从俯视图中数出每排树木的数目n，7排树木的n为：37，34，35，32，15，23，30
(3)冠幅＝R/n，则这7排树木的冠幅为：3.0370，3.3737，3.2314，3.3882，6.5646，4.9032，3.7606
(4)去掉最大值6.5646(m)和最小值3.0370(m)，剩余5组数据求平均值(3.3737+3.2314+3.3882+4.9032+3.7606)/5＝3.7314(m)，则3.7314(m)即为该林分树木东西方向上的平均冠幅。其南北方向上的冠幅计算过程同理，求得其平均冠幅为6.5055(m)。
作为对上述的激光点云中林木参数评估方法的进一步改进，步骤C中树冠中心定位方法如下：
以根部往上取相对于整棵树高度的1/20区域的点云数据(p1,...,ph(1/20))为基点，拟合出h(1/20)条沿树干方向的最佳拟合直线L(v(vx,vy,vz),p0(pj,x,pj,y,pj,z))，其中v为直线的方向向量，p0为点云(p1,...,ph(1/20))中的一点，j＝1,...,h(1/20)；最佳拟合直线即为每株树的树干的所有点云数据到直线L的正交距离的平方和最小；根据
令wi＝p0-pi，ci 2＝|wi|2，pi为整棵橡胶树点云数据中的任意一点；式(3)可以写为：
其中：
设矩阵B为点集(p1,...,ph(1/20))的协方差矩阵乘N：
B＝UTU，其中
则α和β可以表示为：
β(v)＝vTBv，α(v)＝vT(kl)v，其中
式(4)则可以写为：
其中M＝kl-B；
表达式是瑞利商，当v等于矩阵M的最小特征值对应的特征向量时，其表达式的值达到最小值；且k为矩阵M的任意特征值，ξ为对应的特征向量，则：
Mξ＝(kl-B)ξ＝klξ-Bξ＝kξ-λξ＝(k-λ)ξ (8)
由式(8)可知，矩阵B的特征值都与矩阵M的特征值相关联，则当取到矩阵M的最小特征值对应的特征向量时，矩阵B取到最大特征值对应的特征向量；因此，选择矩阵B的最大特征值相对应的特征向量作为使点到直线正交距离和最小的拟合直线的方向向量v，则最佳拟合直线为：
对求解出的h(1/20)条直线求得平均值，得到平均直线作为该单株树的最终最佳拟合直线，其中， 且平均直线L1为：
平均直线L1在高度hheig下的交点即为冠中心点，其中hheig为主枝干高度。
作为对上述的激光点云中林木参数评估方法的进一步改进，步骤E中所述的基于meanshift和分水岭算法进行株株分离的具体步骤：基于树冠中心点方法，即根据树冠的中心点坐标，将该三维坐标投影到二维灰度图上，得到点云数据的聚类图，再结合分水岭算法对该灰度图进行分割，以获取到各分割区域，即树木的总数目，进而实现了对单树木的提取。
作为对上述的激光点云中林木参数评估方法的进一步改进，步骤E中所述的基于meanshift和分水岭算法进行株株分离的具体步骤：基于meanshift方法，根据消除地面高度不平后的点云数据，从根部开始，往上取树总高度的1/7，将在这个阈值范围内的点云数据再次筛选，即在半径r1为0.35m的圆柱体内的点云数据选取出来，并将这些点云数据投影至二维的灰度图上，白点标记的就是每棵树的位置信息；再结合meanshift算法，对该灰度图进行处理，根据树的位置信息对林木进行分类；分类好的结果采用分水岭算法根据分类号进行分割，将点云数据分割为多个分水岭区域，从而实现对单株树的分离。
本发明的有益效果：本发明提出了的研究方法——利用激光雷达采用自下而上的扫描方式，获取树林的点云数据，结合算法实现对橡胶林的株株分离。其具体的过程为：首先对数据做消除地面高度不平操作；其次对点云数据进行枝叶分离，获取树干信息；再次基于枝叶分离的结果获取主枝干倾斜角与树冠中心点；最后基于树冠中心点和meanshift算法的点聚类结果，再结合分水岭算法，对林木进行分割，进而实现对单株树木的提取。整个项目的实验流程如图1所示。
本发明利用背负式激光雷达采用自下而上的扫描方式在一定程度上减少了树冠的遮挡影响；另外，由于树木在在长期风力的胁迫下，树木可能出现倾斜现象，所以本专利采用主枝干作为主要分离特征，测试数据更加准确。本方法基本不受环境影响，计算得到的树高、冠幅、胸径、冠积、倾斜度、树叶分布与叶面积指数参数等树木参数更加符合实际。
附图说明
图1是树木参数评估方法流程图。
图2是实验地——儋州橡胶林的扫描图。
图3是为消除地面高度差设置的平滑窗口示意图。
图4是树冠中心点获取过程示意图，其中(a)为cloudcompare中树的3D图，(b)为冠中点示意图。
图5是林段2未消除地面高度差的图。
图6是林段2消除地面高度差后的图。
图7是林段1几株树的枝叶分离结果与枝干提取图。
图8林段2几株树的枝叶分离结果与枝干提取图。
图9是林段1获取的树冠中心点示意图。
图10是林段2获取的树冠中心点示意图。
图11是林段2基于树冠中心点的点聚类图。
图12是林段2基于树冠中心点的分水岭分割图。
图13是林段2基于meanshift算法的点聚类图。
图14是林段2基于meanshift算法的分水岭分割图。
图15是林段2基于分水岭算法的分聚图(俯视图)。
图16是林段2基于分水岭算法的分聚图(侧视图)。
图17是林段1基于下的中心点比较结果图。
图18是林段1和林段2的东西和南北方向上的冠幅箱体图。
图19是林段1和林段2的主枝干倾斜角α和分支与主枝干间的夹角β的箱体图。
图20是林段1叶面分布雷达图。
图21是林段2叶面分布雷达图。
图22是林段1抽取的三幅鱼目图。
图23是林段2抽取的三幅鱼目图。
图24是林段1和2不同天顶角下计算获取的叶面积指数。
图25是林段1叶面积密度结果图。
图26是林段2叶面积密度结果图。
具体实施方式
1研究区概况与数据获取
1.1研究区概况
研究区选在海南儋州的橡胶实验基地，地处海南岛的西部工业走廊北段(19°32'47.89″N,109°28'29.33″E)。儋州市是海南岛最大的橡胶生产基地，橡胶林众多，橡胶资源及其丰富。儋州市属于热带季风海洋性气候区，夏无酷暑，冬无严寒。阳光充足，气候温和，年均气温23.2℃，极端最高气温为40.0℃。雨量充沛但分布不匀，冬春雨量较稀少，夏秋雨水较充足，年均降雨量1815毫米，5月至10月份为雨季，11月至翌年4月份为旱季。强台风平均每4年一遇，受热带风暴及台风影响，部分树木已出现弯折，且橡胶树的树冠顶并未出现形态结构明显的冠形特征。因此如何精确有效地实现强风力载荷下橡胶树林分中单株提取较为困难。研究区概况如图2所示。
1.2激光雷达数据的获取
本研究使用的原始点云数据的获取装置为SCANLOOK V系列肩背激光雷达，ScanLook是高效便捷的移动空间地理信息获取系统(Mobile Acquisition GeospatialInformation System,MAGIS)，相关参数为：重量为2.5KG,扫描速度为70万点/秒，精度为2公分(30米内)，扫描频率为5Hz～20Hz，最大测程为120m。该系统具有重量轻，体积小，便于携带，高分辨率、高精确，可以在不同移动平台间轻松切换，无论无人机、汽车、船、电动车还是肩背均可轻松作业的优势。数据获取时间为2017年7月30日。作业的温度为-20℃～+45℃，激光扫描仪的波长为905nm，扫描视场为40°×360°，回波次数为2次，输入电压为9VDC～32VDC。获取到两个林段点云数据，分别共有1,159,870和1,326,925个点。林段是对实验样地人为划分出几个区域，每个区域下的林木即为一个林段。
激光雷达数据：如：(339728.6900，2162912.2800，152.0600，1，44，194533.5000，8，1)，其中，前三列表示点云数据在x、y、z坐标轴方向上的数值，第六列表示回波数据，其他列是根据GPS定位而得到的一些标定信息。
2研究方法
2.1消除地面高度不平
获取到的激光雷达数据中，其地面是起伏不平，地面的高度并不是处于水平方向。由于地面高度差的影响，导致一些与水平方向成高度对比的地面，在研究中被检测为树干部分，这对研究的结果会造成很大的误差，也对我们的提取单株树木研究带来一定的难度。为消除这一误差，减少造成单株树木提取工作困难的因素，本文设置了平滑窗口，对原数据进行平滑滤波处理，将所有的值放在同一个高度上，以消除实验数据中的地面高度差。设一窗口下的点云数据为(p0(x0,y0,z0),p1(x1,y1,z1),…,pn(xn,yn,zn))，令zmin＝(z0,z1,…,zn)，另一窗口下点云数据中的任一点集为pi(xi,yi,zi)，要想这两个窗口中的点云数据处于同一地面水平上，令height＝min(zi)-zmin，则zi＝zi-height。随着滑动窗口连续平移，不断消除地形高度差异性，进而使得整片林分处于相同的地形高度上。
以图3为例，黑色框1和灰色框2为两个平滑窗口。其中，黑窗口中高度z最小值为zmin＝143.039993，当窗口滑动到灰色窗口位置时，高度min(zi)＝143.149994。为使地面处于同一高度上，令height＝min(zi)-zmin＝0.110001，则zi＝zi-height＝143.039993，使得两个窗口下的树木便处于相同的地形高度上。窗口按照左右上下的方向滑动，进而整个林分处于同一地面水平线上。
2.2林分枝叶分离及主枝干提取
目前，我们所获取到的仅是点云数据，并没有获取到橡胶树木的具体相关信息，如单株树干的位置、树冠的中心点等。所以，为从大量稠密的点云中提取有效的橡胶树信息，本文首先对其进行主枝干提取操作。主枝干提取，即枝叶分离，这一过程对我们的研究有着很大的作用，更为之后的获取单株橡胶树的信息奠定基础。该操作的具体步骤如下：对于点云数据P中的一点pi＝(xi,yi,zi)T,其中对于点pi，半径为rm内的领域定义为：pj＝(xj,yj,zi)T，且满足条件：||pj-pi||≤rm。点pi邻域的协方差矩阵定义为Cp：
当时，其中μ表示pi邻域所有点在空间中位置的均值。在这里，我们重新定义一个新的坐标系，特征向量代表轴方向，相对应的特征值代表沿轴的点方差。设ek,i是矩阵Cp的特征向量，λk,i是相应的特征值，k＝0,1,2且λ0,i≤λ1,i≤λ2,i。λk,i定量的显示沿轴ek,i方向的数据方差。用e0,i＝{ei,x,ei,y,ei,z}表示最小特征值λ0,i所对应的特征向量，经验证，e0,i也代表着点pi处的法向量。接下来，我们将对应的每个点云pi的方差矩阵的特征值归一化，即：结构张量特征，包括结构张量的平面性特征、熵特征与线性特征，分别计算为
另外，我们使用以下等式计算每个点云pi的法向量分布：
其中e0,j是点云pi的法向量，是点云pi邻域基团的平均法向量分布，且值得注意的是，Vp三个特征值{l0,l1,l2}的相对幅度满足条件：l0≥l1≥l2，这反映出点云pi邻域的几何特征；我们认为当特征值符合l0≈l1≈l2时，基团点云呈现球形状，对应为果实；当特征值分布为l0≥l1≈l2时，基团点云呈现线形状，对应为枝干；而当特征值分布为l0≈l1≥l2时，基团点云呈现面状，对应为叶片。
基于上述分析，一系列关于激光点pi的特征被计算出来，记为:fpi＝{eix,eiy,eiz,c0,i,c1,i,c2,i,l0,i,l1,i,l2,i},其中{eix,eiy,eiz}为点pi的法向量，{c0,i,c1,i,c2,i}为点pi的结构张量特征，而{l0,i,l1,i,l2,i}为点pi的形状特征值。
根据获得的点pi的特征fpi，对橡胶林段进行枝叶分离操作。基于枝叶分离后的枝干数据，按照橡胶树自然生长方向，以树总高度的1/2区域作为主枝干数据(其余部分为分支干数据)，从而获取到两个林段的主枝干信息。
(1)取一棵树中的一组机载激光点云数据的前三列(每个点的x、y、z坐标)，如
(2)对于取得的点云数据pi＝(xi,yi,zi)T，在小于半径rm＝0.585m求得pi的领域：pj＝(xj,yj,zi)T，即则根据2.2节公式(1)可得点pi邻域的协方差矩阵Cp；
(3)对Cp进行特征值分解，eig(Cp)求得Cp的特征值和特征向量；并求得其最小特征值和对应的特征向量；
(4)基于以上结果，根据2.2节公式(2)，求得每个点云pi的法向量分布Vp。并对Vp进行特征值分解，获得Vp的特征值和特征向量。
如，若pi＝(339739.7100,2162925.1900,147.5400)T，求得pi的领域内的点有37个，根据公式(1)求得协方差矩阵
对Cp特征值分解，计算得到
特征向量为
则其最小特征值为0.3123，对应的特征向量为(0.2330,0.8951,-0.3802)T。带入公式(2)求得
对Vp特征值分解，计算得到
特征向量为
同理，计算其他点云法向量分布Vp的特征值为：
(5)基于对点云特征的分析和求得的特征值数据，结合高斯模型聚类分析，对点云进行分类，预先定义为两类(枝和叶)。分类的结果如下：
1，1，1，1，1，1，1，1，1，2
其中，为1的定义为叶子，而为2的定义为枝干。
2.3冠中心定位与主干倾斜度分析
2.3.1冠中心定位
由2.2节可知，本文对原数据进行了枝叶分离操作。现我们已经获取到单株橡胶树的主枝干信息，但我们又发现，大部分的橡胶树的树干是倾斜的。由于主枝干的倾斜，造成冠中心与主枝干有偏差，即使已经分离出单株橡胶树的树干的信息，也无法直接获取树冠的中心点。
为解决此问题，我们根据2.2节获取的枝干数据，以根部往上取相对于整棵树高度的1/20区域的点云数据(p1,...,ph(1/20))为基点，拟合出h(1/20)条沿树干方向的最佳拟合直线L(v(vx,vy,vz),p0(pj,x,pj,y,pj,z))，其中v为直线的方向向量，p0为点云(p1,...,ph(1/20))中的一点，j＝1,...,h(1/20)。要求最佳拟合直线，但拟合直线的方向向量v未知，因此只要确定v即可得拟合直线。为求解v，本文采用瑞利商的方法，对拟合直线的法矢量进行求解。可知，最佳拟合直线即为每株树的树干的所有点云数据到直线L的正交距离的平方和最小：
令wi＝p0-pi，ci 2＝|wi|2，pi为整棵橡胶树点云数据中的任意一点。所以式(3)可以写为：
其中：
设矩阵B为点集(p1,...,ph(1/20))的协方差矩阵乘N：
B＝UTU，其中
则α和β可以表示为：
β(v)＝vTBv，α(v)＝vT(kl)v，其中
式(4)则可以写为：
其中M＝kl-B。(7)
表达式是瑞利商，当v等于矩阵M的最小特征值对应的特征向量时，其表达式的
值达到最小值。且k为矩阵M的任意特征值，ξ为对应的特征向量。则：
Mξ＝(kl-B)ξ＝klξ-Bξ＝kξ-λξ＝(k-λ)ξ。 (8)
由式(8)我们可知，矩阵B的特征值都与矩阵M的特征值相关联，则当取到矩阵M的最小特征值对应的特征向量时，矩阵B取到最大特征值对应的特征向量。因此，我们选择矩阵B的最大特征值相对应的特征向量作为使点到直线正交距离和最小的拟合直线的方向向量v。则最佳拟合直线为：
对求解出的h(1/20)条直线求得平均值，得到平均直线作为该单株树的最终最佳拟合直线，其中， 且平均直线L1为：
另外，平均直线L1在高度hheig下的交点即为冠中心点，其中hheig为主枝干高度。示意图如图4所示。
2.3.2主干倾斜度分析
基于2.3.1求出的拟合直线的平均法向量垂直方向上一方向向量设为则主枝干倾斜度
另外，由2.2节树叶分离后，可知主枝干与分支干之间的夹角为β。
如：以林段2的一棵橡胶树为例，
(1)针对株株分离操作后的树干点云数据(为区分株株分离后的每株不同的树木，根据其分类号，对其点云进行RGB值标记。同一类别号的树木，RGB值是相同的。株株分离的具体过程参见2.4节。)，自树干底部向上取树高1/20区域内点云，共13个点，这13个点的坐标为及RGB值如下：
其中，分类号是指，分水岭操作后，对获得的各分割区域进行标号，此标号即为分类号。(2)以点(26.8026139.9091144.9100)为列，由2.3.1节中公式(3)-(7)知，矩阵
其中，wi＝p0-pi；
其中ci 2＝|wi|2。则矩阵
其中，l为单位矩阵。
(3)基于瑞利商的方法，知最小特征值对应的特征向量即为拟合直线的方向向量。对矩阵M进行特征值分解，根据函数eig(M)，求得矩阵M的特征值，其特征值为10397.0674、64057.4834、69112.9510，取最小特征值为10397.0674，则最小特征值对应的特征向量为(-0.2458-0.45410.8564)。同理，求取其他12个点的最小特征值和特征向量。
则这13个点云的最小特征值及对应的特征向量分别为：
(3)基于13个点云的坐标数据及方向向量，求出13条拟合直线的直线方程。由2.3.1节知，拟合直线的方程为其中(pj,x,pj,y,pj,z)为点坐标，(vx,vy,vz)为特征向量。以第一个点的坐标为(26.8879 139.7503 145.2400)和方向向量为(-0.2820-0.4579 0.8431)为例，则直线方程为：同理，其他12条拟合直线的方程为：
(4)
基于公式求出13个点的平均点，其平均点为(27.0943 139.8683 145.1046)；基于公式求出平均向量，其平均向量为(-0.3038-0.4578 0.8340)。由(3)知，这13条拟合直线的平均直线方程为
(5)在树高的5/3的高度下，得到树冠中心点为(26.3020 138.6744 147.2797)，其中树高取1.5647m。即令求出到的(xy z)即为冠中点。
2.4基于meanshift和分水岭算法的株株分离方法
根据已经获取到的单株橡胶树的树干信息和树冠的中心点信息，我们结合分水岭算法，对橡胶林段进行分类，每棵橡胶树即归为一类，进而实现对橡胶林的株株分离。为更好的分析林段1和林段2受风力侵害程度的差异，在这一过程中，我们采用了两种不同的方法：
一、基于树冠中心点。由2.3节知，已提取出树冠的中心点坐标，将该三维坐标投影到二维灰度图上，得到点云数据的聚类图，再结合分水岭算法对该灰度图进行分割，以获取到各林段的分割区域，即树木的总数目，进也实现了对单株橡胶林的提取。
具体过程如下：
(1)基于整个林段的枝干信息，确定图像大小，设置行和列各为500像素；
(2)选取10棵树的冠中心点：
(3)由冠中点的x和y轴坐标，求得各冠中心点的所对应的像素(x、y轴坐标各比上行和列区间值的最小值，即为所对应的像素)
(4)基于imshow函数，将冠中心点坐标投影至与设置的图像等大小的灰度图上。得到矩阵pic，其中，冠中点不为空时，对应的位置pic矩阵记为1，否则为0；
(5)在灰度图的基础上，求得pic的极大值，及个元素间的距离；再结合空间分水岭算法，即watershed函数，对点云进行等区域分割，共分割为10个区域，即共有10棵树。
二、基于meanshift方法。根据消除地面高度不平后的点云数据，从根部开始，往上取树总高度的1/7，将在这个阈值范围内的点云数据再次筛选，即在半径r1为0.35的圆柱体内的点云数据选取出来，并将这些点云数据投影至二维的灰度图上，白点标记的就是每棵树的位置信息。再结合meanshift算法，对该灰度图进行处理，根据树的位置信息对林段进行分类。分类好的结果采用分水岭算法进行分割，将点云数据分割为多个分水岭区域，并对每个分割区域进行标号(即分类号)，从而实现对单株树的分离。其中林段1共155棵橡胶树；林段2共142棵橡胶树。
具体过程如下：
(1)取树高1/7区域，半径r1为0.35的圆柱体内的点云数据：
(2)求各点所对应的像素：
基于imshow函数绘制灰度图，得到矩阵img；
(3)在灰度图的基础上，使用meanshift的方法对img进行处理。使用meanshiftseg对img进行分割，得到初步的分类。找到每一类的中心位置,并确定中心点间的距离。设置阈值为12，将距离小于此阈值的中心点合为一类，得到矩阵pic；
(4)求得pic的极大值，及个元素间的距离，再结合空间分水岭算法，即watershed函数，对点云进行等区域分割，实现对橡胶林段的株株分离。
3结果与分析
基于第2节的方法，本文已实现对单株橡胶树的分离和提取，各研究阶段所获取到的结果如图5-17所示。依据获取的结果，可知林段1的冠幅要大于林段2。另外，由图15可知，林段1基于下的中心点均沿着树干倾斜方向分布，由图17中*、□的分布，可以明确林段1中的各橡胶树的倾斜程度，即林段1在台风载荷下的侵害程度。
本文从林段1和林段2中分别选取十五棵树，基于convexhull算法求得凸包体积，将所得的凸包体积的平均值，作为两个林段的平均冠体积。另外，本文对橡胶树的冠幅、倾斜度、平均胸径和平均冠体积等进行了参数提取，两个林段的各参数准确度的定量评估见表1。
表1不同林段分类的总体准确性评估
由表1中数据，可以得到两林段的东西和南北方向上的冠幅及主枝干倾斜角α和分支与主枝干间的夹角β的箱体图，如图18、19所示。
4讨论
4.1叶面分布雷达图分析
为研究两林段的叶面分布，本文从林段1和林段2中各选取30棵橡胶树，分别对其叶子数据求得各自分布的雷达图，结果如图20、21所示。由图可知，林段1叶面主要分布在东西方向上，而林段2在四个方向上均匀分布。当风力方向如图所示时，由于林段1的叶面主要分布在东西方向上，而南北方向上的叶面几乎没有，这就造成了林段1的叶子易受风力侵害；对于林段2，无论台风是哪一方向，由于林段2的叶子在四个方向上均匀分布，因此，林段2的叶子不受风力侵害。由此表明，林段1的橡胶树受风力侵害程度要大于林段2的橡胶树。
4.2基于鱼目图像验证
为进一步验证林段1受风力侵害程度大于林段2结论的正确性，本文基于随机抽取的三株橡胶树的鱼目图像，在不同天顶角下计算并获取叶面积指数，其结果如图22、23所示。图24表明，林段1的叶面积指数总是高于林段2的叶面积指数。由此再结合两林段的叶面分布，可以验证林段1受风力侵害程度比林段2要大，即林段1的橡胶树比林段2的橡胶树易倒伏。
4.3叶面积密度分析
为进一步研究两橡胶林段的叶面分布，本文基于每株树的叶子点云数据与冠体积，求得相应的叶面积密度。在此基础上，设置一阈值，将低于该阈值的点云筛选出来，获得其林分中叶面积密度较低的林木。其中，虚线圈出来的为低于该阈值的点云数据。结果如图25、26所示。
5总结
在预处理人背负激光雷达Lidar数据的基础上，本研究对两个橡胶林段进行单株提取，分别进行枝叶分离操作和主枝干提取，基于主枝干的信息，进行拟合最佳直线获取单株树冠的冠中心点信息，根据冠中心点和meanshit算法再结合分水岭算法，分别对单株橡胶树进行定位和提取，进而实现株株分离。首先，本文对两种方法下的树冠定位结果进行比较。接着，本文从两个林段中分别随机选取十五棵树，求得它们的叶面分布雷达图并进行比较分析，我们发现林段1叶面主要分布在东西方向上，而林段2在四个方向上均匀分布。由此我们得到以下结论：林段1受风力侵害程度比林段2要大，且林段1的橡胶树比林段2的橡胶树易倒伏。本文又基于三棵橡胶树的鱼目图像获取两林段的叶面积指数，对以上结论进行验证。验证结果表明林段1的叶面积指数总是要高于林段2的叶面积指数，即结论得以验证。
单株提取的研究工作一般采用的研究方法是自上而下的扫描方式，但这仅对树冠顶部凸起的情况有效。基于这方面的考虑，本文选用自下而上的移动式扫描方式，针对Lidar数据进行研究，最终实现对橡胶林段的单株提取。此方法虽可以有效的解决树冠特这不明显的问题，但还是存在一定的遮挡。后期研究内容将侧重于这方面的改进，期望在精度上有一个较大的提升，从而为实际应用提供坚实的理论基础。

Claims (4)

1.激光点云中林木参数评估方法，其特征是：包括以下步骤：
步骤A:获取点云数据：以激光雷达对树林自下而上扫描，获取树林的点云数据pi(xi，yi，zi)，xi、yi、zi分别表示点pi在x、y、z坐标轴方向上的坐标值；
步骤B:对树木株叶分离：
对于点云数据中的一点pi＝(xi,yi,zi)T,其中对于点pi，半径为rm内的领域定义为：pj＝(xj,yj,zi)T，且满足条件：||pj-pi||≤rm；
点pi邻域的协方差矩阵定义为Cp：
当时，其中μ表示pi邻域所有点在空间中位置的均值；
重新定义一个新的坐标系，特征向量代表轴方向，相对应的特征值代表沿轴的点方差；设ek,i是矩阵Cp的特征向量，λk,i是相应的特征值，k＝0,1,2且λ0,i≤λ1,i≤λ2,i；λk,i定量的显示沿轴ek,i方向的数据方差；
用e0,i＝{ei,x,ei,y,ei,z}表示最小特征值λ0,i所对应的特征向量，e0,i也代表着点pi处的法向量；
将对应的每个点云pi的方差矩阵的特征值归一化，即：
计算结构张量特征，包括结构张量的平面性特征(c0,i)、熵特征(c1,i)与线性特征(c2,i)，分别计算为
使用以下等式计算每个点云pi的法向量分布：
其中e0,j是点云pi的法向量，是点云pi邻域基团的平均法向量分布，且
Vp三个特征值{l0,l1,l2}的相对幅度满足条件：l0≥l1≥l2，这反映出点云pi邻域的几何特征；
基于上述分析，一系列关于点云pi的特征被计算出来，记为:其中{eix,eiy,eiz}为点pi的法向量，{c0,i,c1,i,c2,i}为点pi的结构张量特征，而{l0,i,l1,i,l2,i}为点pi的形状特征向量；
当特征值符合l0≈l1≈l2时，基团点云呈现球形状，对应为果实；当特征值分布为l0≥l1≈l2时，基团点云呈现线形状，对应为枝干；而当特征值分布为l0≈l1≥l2时，基团点云呈现面状，对应为叶片；也就是说，根据获得的点pi的特征对树木进行枝叶分离操作；
基于枝叶分离后的枝干数据，按照树木自然生长方向，以树总高度的1/2区域作为主枝干数据，其余部分为分支干数据；
步骤C:树冠中心定位，确定树冠的中心点坐标；
步骤D:倾斜度分析，计算出则主枝干倾斜度α、主枝干与分支干之间的夹角β；
步骤E:基于meanshift和分水岭算法进行株株分离，得到树木总棵数，并计算得出树高、胸径、冠体积、冠幅。
2.如权利要求1所述的激光点云中林木参数评估方法，其特征是：
步骤C中树冠中心定位方法如下：
以根部往上取相对于整棵树高度的1/20区域的点云数据(p1,...,ph(1/20))为基点，拟合出h(1/20)条沿树干方向的最佳拟合直线L(v(vx,vy,vz),p0(pj,x,pj,y,pj,z))，其中v为直线的方向向量，p0为点云(p1,...,ph(1/20))中的一点，j＝1,...,h(1/20)；最佳拟合直线即为每株树的树干的所有点云数据到直线L的正交距离的平方和最小；根据
<mrow> <mi>arg</mi> <mi> </mi> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>v</mi> <mo>(</mo> <mrow> <msub> <mi>v</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mfrac> <mrow> <msup> <mrow> <mo>|</mo> <mrow> <msub> <mi>p</mi> <mn>0</mn> </msub> <mo>-</mo> <msub> <mi>p</mi> <mi>i</mi> </msub> </mrow> <mo>|</mo> </mrow> <mn>2</mn> </msup> <msup> <mrow> <mo>|</mo> <mi>v</mi> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <msup> <mrow> <mo>(</mo> <mo>(</mo> <mrow> <msub> <mi>p</mi> <mn>0</mn> </msub> <mo>-</mo> <msub> <mi>p</mi> <mi>i</mi> </msub> </mrow> <mo>)</mo> <mo>&amp;CenterDot;</mo> <mi>v</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <msup> <mrow> <mo>|</mo> <mi>v</mi> <mo>|</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
令wi＝p0-pi，ci 2＝|wi|2，pi为整棵橡胶树点云数据中的任意一点；式(3)可以写为：
<mrow> <mi>arg</mi> <mi> </mi> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>v</mi> <mo>(</mo> <mrow> <msub> <mi>v</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </msubsup> <msubsup> <mi>c</mi> <mi>i</mi> <mn>2</mn> </msubsup> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>v</mi> <mo>-</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>x</mi> </msub> <msub> <mi>w</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>x</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>v</mi> <mi>y</mi> </msub> <msub> <mi>w</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>y</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> <msub> <mi>w</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>z</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>v</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mi>&amp;alpha;</mi> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>&amp;beta;</mi> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>v</mi> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
其中：
设矩阵B为点集(p1,...,ph(1/20))的协方差矩阵乘N：
B＝UTU，其中
则α和β可以表示为：
β(v)＝vTBv，α(v)＝vT(kl)v，其中
式(4)则可以写为：
<mrow> <mi>arg</mi> <mi> </mi> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>v</mi> <mo>(</mo> <mrow> <msub> <mi>v</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mi>l</mi> <mo>)</mo> </mrow> <mi>v</mi> <mo>-</mo> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>B</mi> <mi>v</mi> </mrow> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>v</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mi>l</mi> <mo>-</mo> <mi>B</mi> <mo>)</mo> </mrow> <mi>v</mi> </mrow> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>v</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>M</mi> <mi>v</mi> </mrow> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>v</mi> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
其中M＝kl-B；
表达式是瑞利商，当v等于矩阵M的最小特征值对应的特征向量时，其表达式的值达到最小值；且k为矩阵M的任意特征值，ξ为对应的特征向量，则：
Mξ＝(kl-B)ξ＝klξ-Bξ＝kξ-λξ＝(k-λ)ξ (8)
由式(8)可知，矩阵B的特征值都与矩阵M的特征值相关联，则当取到矩阵M的最小特征值对应的特征向量时，矩阵B取到最大特征值对应的特征向量；因此，选择矩阵B的最大特征值相对应的特征向量作为使点到直线正交距离和最小的拟合直线的方向向量v，则最佳拟合直线为：
<mrow> <mfrac> <mrow> <mi>x</mi> <mo>-</mo> <msub> <mi>p</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>x</mi> </mrow> </msub> </mrow> <msub> <mi>v</mi> <mi>x</mi> </msub> </mfrac> <mo>=</mo> <mfrac> <mrow> <mi>y</mi> <mo>-</mo> <msub> <mi>p</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>y</mi> </mrow> </msub> </mrow> <msub> <mi>v</mi> <mi>y</mi> </msub> </mfrac> <mo>=</mo> <mfrac> <mrow> <mi>z</mi> <mo>-</mo> <msub> <mi>p</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>z</mi> </mrow> </msub> </mrow> <msub> <mi>v</mi> <mi>z</mi> </msub> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
对求解出的h(1/20)条直线求得平均值，得到平均直线作为该单株树的最终最佳拟合直线，其中， 且平均直线L1为：
<mrow> <mfrac> <mrow> <mi>x</mi> <mo>-</mo> <mover> <msub> <mi>x</mi> <mrow> <mi>a</mi> <mi>v</mi> <mi>g</mi> </mrow> </msub> <mo>&amp;OverBar;</mo> </mover> </mrow> <mover> <mi>x</mi> <mo>&amp;OverBar;</mo> </mover> </mfrac> <mo>=</mo> <mfrac> <mrow> <mi>y</mi> <mo>-</mo> <mover> <msub> <mi>y</mi> <mrow> <mi>a</mi> <mi>v</mi> <mi>g</mi> </mrow> </msub> <mo>&amp;OverBar;</mo> </mover> </mrow> <mover> <mi>y</mi> <mo>&amp;OverBar;</mo> </mover> </mfrac> <mo>=</mo> <mfrac> <mrow> <mi>z</mi> <mo>-</mo> <mover> <msub> <mi>z</mi> <mrow> <mi>a</mi> <mi>v</mi> <mi>g</mi> </mrow> </msub> <mo>&amp;OverBar;</mo> </mover> </mrow> <mover> <mi>z</mi> <mo>&amp;OverBar;</mo> </mover> </mfrac> <mo>=</mo> <mfrac> <mn>5</mn> <mn>3</mn> </mfrac> <msub> <mi>h</mi> <mrow> <mi>h</mi> <mi>e</mi> <mi>i</mi> <mi>g</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
平均直线L1在高度hheig下的交点即为冠中心点，其中hheig为主枝干高度。
3.如权利要求1所述的激光点云中林木参数评估方法，其特征是：步骤E中所述的基于meanshift和分水岭算法进行株株分离的具体步骤：基于树冠中心点方法，即根据树冠的中心点坐标，将该三维坐标投影到二维灰度图上，得到点云数据的聚类图，再结合分水岭算法对该灰度图进行分割，以获取到各分割区域，即树木的总数目，进而实现了对单树木的提取。
4.如权利要求1所述的激光点云中林木参数评估方法，其特征是：步骤E中所述的基于meanshift和分水岭算法进行株株分离的具体步骤：基于meanshift方法，根据消除地面高度不平后的点云数据，从根部开始，往上取树总高度的1/7，将在这个阈值范围内的点云数据再次筛选，即在半径r1为0.35m的圆柱体内的点云数据选取出来，并将这些点云数据投影至二维的灰度图上，白点标记的就是每棵树的位置信息；再结合meanshift算法，对该灰度图进行处理，根据树的位置信息对林木进行分类；分类好的结果采用分水岭算法根据分类号进行分割，将点云数据分割为多个分水岭区域，从而实现对单株树的分离。
CN201710955957.5A 2017-10-15 2017-10-15 激光点云中林木参数评估方法 Active CN107705309B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710955957.5A CN107705309B (zh) 2017-10-15 2017-10-15 激光点云中林木参数评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710955957.5A CN107705309B (zh) 2017-10-15 2017-10-15 激光点云中林木参数评估方法

Publications (2)

Publication Number Publication Date
CN107705309A true CN107705309A (zh) 2018-02-16
CN107705309B CN107705309B (zh) 2020-12-04

Family

ID=61183780

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710955957.5A Active CN107705309B (zh) 2017-10-15 2017-10-15 激光点云中林木参数评估方法

Country Status (1)

Country Link
CN (1) CN107705309B (zh)

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108804854A (zh) * 2018-06-29 2018-11-13 河北农业大学 一种果树机械修枝的树体结构参数估计方法及其系统
CN108921885A (zh) * 2018-08-03 2018-11-30 南京林业大学 一种综合三类数据源联合反演森林地上生物量的方法
CN109002418A (zh) * 2018-06-20 2018-12-14 厦门大学 基于体素生长和地面激光点云的树木胸径自动计算方法
CN109035286A (zh) * 2018-06-29 2018-12-18 北京农业信息技术研究中心 植物三维点云数据高通量获取方法及装置
CN109299739A (zh) * 2018-09-26 2019-02-01 速度时空信息科技股份有限公司 基于法向量的曲面拟合对车载激光点云滤波的方法
CN109446691A (zh) * 2018-11-08 2019-03-08 南京林业大学 基于激光点云与空气动力学的活立木抗风性能分析方法
CN109856139A (zh) * 2019-01-07 2019-06-07 北京林业大学 树木叶片数量估测方法及系统
CN110378909A (zh) * 2019-06-24 2019-10-25 南京林业大学 基于Faster R-CNN的面向激光点云的单木分割方法
CN111220967A (zh) * 2020-01-02 2020-06-02 小狗电器互联网科技（北京）股份有限公司 一种激光雷达数据有效性的检测方法与装置
CN111709986A (zh) * 2020-05-18 2020-09-25 中国能源建设集团江苏省电力设计院有限公司 一种基于激光点云的输电线下林木统计方法
WO2020192586A1 (zh) * 2019-03-22 2020-10-01 南京农业大学 一种田间小麦茎蘖数提取方法
CN111966952A (zh) * 2020-07-21 2020-11-20 电子科技大学 一种叶面积指数数据过滤方法
CN112837309A (zh) * 2021-03-02 2021-05-25 华南农业大学 果树冠层靶标识别装置、方法、计算设备及存储介质
CN112946669A (zh) * 2021-02-05 2021-06-11 国际竹藤中心 基于地基激光雷达的毛竹林数量识别方法及装置
CN113269825A (zh) * 2021-04-06 2021-08-17 云南师范大学 基于地基激光雷达技术林木胸径值提取的方法
CN114396975A (zh) * 2022-01-05 2022-04-26 深圳市金阳盛城市服务集团有限公司 一种树木抗倒伏监控预警方法、装置、电子设备及介质
CN114862872A (zh) * 2022-05-10 2022-08-05 桂林理工大学 一种基于Faster R-CNN的红树林单木分割方法
CN117475002A (zh) * 2023-12-27 2024-01-30 青岛亿联建设集团股份有限公司 基于激光扫描技术的建筑倾斜度测量方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1300999A1 (en) * 2001-10-05 2003-04-09 Siemens Information and Communication Networks S.p.A. Channel estimation
CN102667835A (zh) * 2009-09-30 2012-09-12 公共风险管理有限责任合作公司 用于定位目标树的树干的方法和系统
CN104156964A (zh) * 2014-08-14 2014-11-19 陈荣元 一种综合mrf和贝叶斯网络的遥感影像区域分割方法
CN104504709A (zh) * 2014-12-28 2015-04-08 大连理工大学 一种基于特征球的室外场景三维点云数据的分类方法
CN105373814A (zh) * 2014-08-26 2016-03-02 云挺 一种基于激光点云数据的真实阔叶树器官分类识别方法
CN105957076A (zh) * 2016-04-27 2016-09-21 武汉大学 一种基于聚类的点云分割方法及系统

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1300999A1 (en) * 2001-10-05 2003-04-09 Siemens Information and Communication Networks S.p.A. Channel estimation
CN102667835A (zh) * 2009-09-30 2012-09-12 公共风险管理有限责任合作公司 用于定位目标树的树干的方法和系统
CN104156964A (zh) * 2014-08-14 2014-11-19 陈荣元 一种综合mrf和贝叶斯网络的遥感影像区域分割方法
CN105373814A (zh) * 2014-08-26 2016-03-02 云挺 一种基于激光点云数据的真实阔叶树器官分类识别方法
CN104504709A (zh) * 2014-12-28 2015-04-08 大连理工大学 一种基于特征球的室外场景三维点云数据的分类方法
CN105957076A (zh) * 2016-04-27 2016-09-21 武汉大学 一种基于聚类的点云分割方法及系统

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
ALFRED O. HERO等: "Cutoff Rate and Signal Design for the Quasi-Static Rayleigh-Fading Space–Time Channel", 《IEEE TRANSACTIONS ON INFORMATION THEORY》 *
JING LIU等: "An individual tree crown delineation method based on multi-scale segmentation of imagery", 《ISPRS JOURNAL OF PHOTOGRAMMETRY AND REMOTE SENSING》 *
TING YUN等: "A Novel Approach for Retrieving Tree Leaf Area from Ground-Based LiDAR", 《REMOTE SENSING》 *
李增元等: "激光雷达森林参数反演研究进展", 《遥感学报》 *
王学春等: "基于激光点云数据的活立木林学参数挖掘", 《科学技术与工程》 *

Cited By (28)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109002418A (zh) * 2018-06-20 2018-12-14 厦门大学 基于体素生长和地面激光点云的树木胸径自动计算方法
CN109035286A (zh) * 2018-06-29 2018-12-18 北京农业信息技术研究中心 植物三维点云数据高通量获取方法及装置
CN108804854B (zh) * 2018-06-29 2022-11-22 河北农业大学 一种果树机械修枝的树体结构参数估计方法及其系统
CN108804854A (zh) * 2018-06-29 2018-11-13 河北农业大学 一种果树机械修枝的树体结构参数估计方法及其系统
CN108921885B (zh) * 2018-08-03 2020-05-12 南京林业大学 一种综合三类数据源联合反演森林地上生物量的方法
CN108921885A (zh) * 2018-08-03 2018-11-30 南京林业大学 一种综合三类数据源联合反演森林地上生物量的方法
CN109299739A (zh) * 2018-09-26 2019-02-01 速度时空信息科技股份有限公司 基于法向量的曲面拟合对车载激光点云滤波的方法
CN109446691A (zh) * 2018-11-08 2019-03-08 南京林业大学 基于激光点云与空气动力学的活立木抗风性能分析方法
CN109856139A (zh) * 2019-01-07 2019-06-07 北京林业大学 树木叶片数量估测方法及系统
CN109856139B (zh) * 2019-01-07 2021-08-06 北京林业大学 树木叶片数量估测方法及系统
WO2020192586A1 (zh) * 2019-03-22 2020-10-01 南京农业大学 一种田间小麦茎蘖数提取方法
US20220189053A1 (en) * 2019-03-22 2022-06-16 Nanjing Agricultural University Method of extracting number of stem and tillering for wheat under field condition
CN110378909A (zh) * 2019-06-24 2019-10-25 南京林业大学 基于Faster R-CNN的面向激光点云的单木分割方法
CN110378909B (zh) * 2019-06-24 2023-06-06 南京林业大学 基于Faster R-CNN的面向激光点云的单木分割方法
CN111220967A (zh) * 2020-01-02 2020-06-02 小狗电器互联网科技（北京）股份有限公司 一种激光雷达数据有效性的检测方法与装置
CN111709986A (zh) * 2020-05-18 2020-09-25 中国能源建设集团江苏省电力设计院有限公司 一种基于激光点云的输电线下林木统计方法
CN111709986B (zh) * 2020-05-18 2023-11-03 中国能源建设集团江苏省电力设计院有限公司 一种基于激光点云的输电线下林木统计方法
CN111966952B (zh) * 2020-07-21 2023-04-18 电子科技大学 一种叶面积指数数据过滤方法
CN111966952A (zh) * 2020-07-21 2020-11-20 电子科技大学 一种叶面积指数数据过滤方法
CN112946669A (zh) * 2021-02-05 2021-06-11 国际竹藤中心 基于地基激光雷达的毛竹林数量识别方法及装置
CN112837309A (zh) * 2021-03-02 2021-05-25 华南农业大学 果树冠层靶标识别装置、方法、计算设备及存储介质
CN112837309B (zh) * 2021-03-02 2023-10-20 华南农业大学 果树冠层靶标识别装置、方法、计算设备及存储介质
CN113269825A (zh) * 2021-04-06 2021-08-17 云南师范大学 基于地基激光雷达技术林木胸径值提取的方法
CN114396975A (zh) * 2022-01-05 2022-04-26 深圳市金阳盛城市服务集团有限公司 一种树木抗倒伏监控预警方法、装置、电子设备及介质
CN114396975B (zh) * 2022-01-05 2023-05-23 深圳市金阳盛城市服务集团有限公司 一种树木抗倒伏监控预警方法、装置、电子设备及介质
CN114862872A (zh) * 2022-05-10 2022-08-05 桂林理工大学 一种基于Faster R-CNN的红树林单木分割方法
CN114862872B (zh) * 2022-05-10 2024-05-07 桂林理工大学 一种基于Faster R-CNN的红树林单木分割方法
CN117475002A (zh) * 2023-12-27 2024-01-30 青岛亿联建设集团股份有限公司 基于激光扫描技术的建筑倾斜度测量方法

Also Published As

Publication number Publication date
CN107705309B (zh) 2020-12-04

Similar Documents

Publication Publication Date Title
CN107705309A (zh) 激光点云中林木参数评估方法
CN112381861B (zh) 一种基于地基激光雷达的林地点云数据配准和分割方法
Hosoi et al. 3-D voxel-based solid modeling of a broad-leaved tree for accurate volume estimation using portable scanning lidar
Guan et al. Integration of orthoimagery and lidar data for object-based urban thematic mapping using random forests
CN103364781B (zh) 基于遥感数据与地理信息系统的粮田地面参照点筛选方法
Zhou et al. Individual tree parameters estimation for plantation forests based on UAV oblique photography
Korpela et al. Variation and directional anisotropy of reflectance at the crown scale—Implications for tree species classification in digital aerial images
Dinuls et al. Tree species identification in mixed Baltic forest using LiDAR and multispectral data
González-Yebra et al. Methodological proposal to assess plastic greenhouses land cover change from the combination of archival aerial orthoimages and Landsat data
CN108932521B (zh) 一种基于深度学习的农作物分类方法及系统
CN111767865A (zh) 一种利用航拍影像和激光数据反演红树林生物量的方法
CN105335966A (zh) 基于局域同质性指标的多尺度遥感影像分割方法
CN104881868A (zh) 植物群落空间结构提取方法
CN109447089A (zh) 基于超分技术的高分辨率北极海冰类型提取方法
Rodriguez-Galiano et al. Incorporating the downscaled Landsat TM thermal band in land-cover classification using random forest
CN109409429A (zh) 一种基于激光雷达点云数据的树种分类方法
Li et al. An iterative-mode scan design of terrestrial laser scanning in forests for minimizing occlusion effects
CN115880487A (zh) 基于深度学习方法的林木激光点云枝叶分离方法
CN106780586A (zh) 一种基于地面激光点云的太阳能潜力评估方法
CN112669363A (zh) 城市绿地三维绿量测算方法
CN107705344A (zh) 激光扫描环境点云数据中树冠模型提取方法
Chen et al. 3D model construction and ecological environment investigation on a regional scale using UAV remote sensing
Xu et al. Vegetation information extraction in karst area based on UAV remote sensing in visible light band
CN109635828A (zh) 一种生态保护红线区典型地理国情要素提取系统及方法
CN117114147A (zh) 基于雷达和卫星遥感的森林植被碳储量估算的方法和装置

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
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20210630

Address after: 210000 No.8 Huayuan Road, Xuanwu District, Nanjing, Jiangsu Province

Patentee after: Nanjing Maoting Information Technology Co.,Ltd.

Address before: Nanjing City, Jiangsu province 210037 Longpan Road No. 159

Patentee before: NANJING FORESTRY University

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220826

Address after: Room 808, Building 4, No. 245, Binkang Road, Changhe Street, Binjiang District, Hangzhou City, Zhejiang Province, 310000

Patentee after: Zhejiang Zhongnan Green Carbon Technology Co.,Ltd.

Address before: 210000 No.8 Huayuan Road, Xuanwu District, Nanjing, Jiangsu Province

Patentee before: Nanjing Maoting Information Technology Co.,Ltd.

TR01 Transfer of patent right