CN110033484B - 一种结合uav影像和tls点云的高郁闭森林样地树高提取方法 - Google Patents

一种结合uav影像和tls点云的高郁闭森林样地树高提取方法 Download PDF

Info

Publication number
CN110033484B
CN110033484B CN201910321231.5A CN201910321231A CN110033484B CN 110033484 B CN110033484 B CN 110033484B CN 201910321231 A CN201910321231 A CN 201910321231A CN 110033484 B CN110033484 B CN 110033484B
Authority
CN
China
Prior art keywords
point
point cloud
sample plot
key
points
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.)
Expired - Fee Related
Application number
CN201910321231.5A
Other languages
English (en)
Other versions
CN110033484A (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.)
Fuzhou University
Original Assignee
Fuzhou 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 Fuzhou University filed Critical Fuzhou University
Priority to CN201910321231.5A priority Critical patent/CN110033484B/zh
Publication of CN110033484A publication Critical patent/CN110033484A/zh
Application granted granted Critical
Publication of CN110033484B publication Critical patent/CN110033484B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/251—Fusion techniques of input or preprocessed data
    • 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
    • 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/64—Analysis of geometric attributes of convexity or concavity
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
    • 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/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10044—Radar image
    • 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/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging
    • 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)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geometry (AREA)
  • Multimedia (AREA)
  • Image Analysis (AREA)

Abstract

本发明涉及一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，结合无人机影像和地面激光雷达点云的森林样地树高提取方法：同时获取高郁闭度森林样地无人机（UAV：Unmanned aerial vehicle）影像点云与该样地的地面激光雷达（TLS：Terrestrial Laser Scanning）点云数据。通过尺度不变特征变换（SIFT）特征点匹配算法实现两种数据匹配融合。基于融合的点云数据，采用改进的局部极大值算法提取样地的单木树高。本发明可实现同时从林下和冠层两个方向来测量森林样地参数，为更全面、更精确的估算亚热带高郁闭度森林样地的树高提供了一种新的途径。

Description

一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法
技术领域
本发明涉及森林树高研究技术领域，特别是一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法。
背景技术
森林是陆地上面积分布最广、组成结构最复杂且物质资源最丰富的生态系统，能够及时清晰的了解森林信息变化情况，高效地进行森林资源管理是非常重要的，所以国内外开展了大量的森林资源调查的工作。树高是森林结构的重要参数，是估算森林蓄积量等其它参数的基础。
在普通的森林区域，TLS可以提供精细的森林直立面结构数据并可精确估算树高等参数。随着航空遥感技术的发展，低成本无人机(Unmanned Aerial Vehicle： UAV)由于其灵活和机动性，其配备的传感器可以实时获取高空间、时间分辨率的像片数据，通过运动结构多视图立体摄影(SfM-MVS)的技术方法可将获取的具有一定重叠度的序列像片数据生成三维点云数据(Ullman，1979；Lowe，2004； Dandois，2013)。从三维立体视角观察整个森林的概况，通过获取精细的地面点数据，为树高等参数的估算提供前提条件(Kattenborn等，2014)。
然而在高郁闭度的森林区域，UAV影像点云无法穿透冠层到达地面，对于无法获取地面点的区域，其树高等参数的估算受到限制。而在同一区域进行TLS点云数据获取时，由于复杂的森林结构或扫描仪自身条件限制等原因，其区域内树木信息不能被完整扫描，部分单木冠层信息会严重缺失。Jung(2016)将在不同角度方向存在数据缺失的TLS点云和ALS点云进行结合，提高了森林单木树高、冠幅等参数的估算精度。但是ALS点云获取成本较高，且实时性较差。
发明内容
有鉴于此，本发明的目的是提出一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，可实现同时从林下和冠层两个方向来测量森林样地参数，为更全面、更精确的估算亚热带高郁闭度森林样地的树高提供了一种新的途径。
本发明采用以下方案实现：一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，具体包括以下步骤：
步骤S1：同时获取高郁闭度森林样地无人机(UAV)影像点云与该样地的地面激光雷达(TLS)点云数据并生成数字表面模型(DSM)；
步骤S2：估算高郁闭度森林样地无人机(UAV)影像点云与该样地的地面激光雷达(TLS)点云数据的数字表面模型DSM上所有点的曲率，存储DSM所有点云尺度不变的局部邻近点，将曲率值高于其邻近点曲率值的点作为关键点；
步骤S3：建立关键点描述符区域，生成无人机影像点云和地面激光雷达样地点云的DSM描述符SIFT；
步骤S4：匹配最相似的关键点描述符，进行关键点坐标转换，实现无人机影像点云和地面激光雷达样地点云数据匹配融合；
步骤S5：基于无人机影像点云和地面激光雷达样地点云的融合数据，进行样地地面点分割，建立样地冠层高度模型CHM；
步骤S6：基于冠层高度模型CHM，采用改进的局部最大值算法提取样地的单木树高。
本发明结合无人机影像和地面激光雷达点云的森林样地树高提取方法：同时获取高郁闭度森林样地无人机(UAV：Unmanned aerial vehicle)影像点云与该样地的地面激光雷达(TLS：Terrestrial Laser Scanning)点云数据。通过尺度不变特征变换(SIFT)特征点匹配算法实现两种数据匹配融合。基于融合的点云数据，采用改进的局部极大值算法提取样地的单木树高。本发明可实现同时从林下和冠层两个方向来测量森林样地参数，为更全面、更精确的估算亚热带高郁闭度森林样地的树高提供了一种新的途径。
进一步地，步骤S1中，设置森林样地的郁闭度值大于0.7；进行无人机影像数据采集时，影像航向和旁向重叠度分别设置为大于等于30％和大于等于70％，影像分辨率为设置为小于等于3cm；对得到的UAV像对，利用运动结构多视图立体摄影(SfM-MVS)技术和算法生成三维点云数据。
进一步地，步骤S2具体包括以下步骤：
步骤S21：估算两种点云数据的数字表面模型DSM上所有点的表面曲率，利用协方差的特征分析有3个特征向量eμ和其对应的特征值λμ；其中eμ定义一个局部三维正交参考坐标系，而λμ定义每个eμ的大小；根据每点的特征值λ0、λ1、λ2，将局部DSM表面量化为一个近似曲面的曲率，计算公式如下：
Figure BDA0002034760890000031
Figure BDA0002034760890000032
式中，λ0表示候选关键点的表面法线，即是在三维表面上局部切平面的最小变化方向；
步骤S22：步骤S21中各个点的邻近点云半径范围是不确定的，在0-1之间的归一化点云坐标系统中，半径尺度范围是0.1-0.5，增量是0.05，将最大曲率值作为关键点邻近点云的半径值；如果候选点的曲率值大于其相邻点云范围内所有点的曲率值，那么将该候选点作为一个关键点。
进一步地，步骤S3具体包括以下步骤：
步骤S31：计算每个关键点的方向，依照这个方向再做进一步的计算；利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性；其中，点(x，y)处梯度的模值和方向计算公式分别如下：
Figure BDA0002034760890000041
Figure BDA0002034760890000042
式中，L所用的尺度为每个关键点各自所在的尺度；
步骤S32：每个关键点有三个信息：位置、所处尺度以及方向，以关键点为中心，确定一个SIFT特征区域；
步骤S33：在关键点邻域内，对其中点的梯度做高斯加权；其中邻域分为4×4 个子区域，每个子区域生成一个描述子，一个描述子中涉及8个方向，最终每个关键点形成128维的SIFT特征向量。
进一步地，步骤S4具体包括以下步骤：
步骤S41：采用关键点特征向量的欧式距离作为两幅图像中关键点的相似性判定度量，两个样地的数据分别生成相应的描述子，将两图中所有尺度的描述子进行分配，在描述符特征空间中，利用关键点与其第一和第二最邻近点的欧氏距离来匹配每个关键点；
步骤S42：在描述符空间使用最邻近点特征匹配时，如果两个图像关键点之间的相对距离dist1小于用户设置的阈值，即在边界内搜索不到另一幅图像的最邻近关键点，则关键点的描述符D1不能与另外一幅图像关键点描述符Dtrue匹配；如果边界内只有一个可匹配的点，该点即可作为最佳匹配的点；
步骤S43：在描述符空间使用最邻近点特征匹配时，如果在一定的边界范围内，
Figure BDA0002034760890000051
比较
Figure BDA0002034760890000052
更接近D1，使用最邻近距离
Figure BDA0002034760890000053
的
Figure BDA0002034760890000054
作为最佳匹配，
Figure BDA0002034760890000055
是第二邻近匹配；
步骤S44：在描述符空间使用最邻近点特征匹配时，如果在一定的边界范围内，
Figure BDA0002034760890000056
和
Figure BDA0002034760890000057
的接近到使用步骤S43无法有效区分时，则使用最邻近距离比法，只有最邻近距离比NNDR小于用户定义的值τ，匹配点才会被接受，其中NNDR的计算公式如下所示：
Figure BDA0002034760890000058
步骤S45：在确认匹配点以后，使用RANSAC方法剔除错误匹配的点；
步骤S46：估算无人机影像点云和地面激光雷达样地点云表面对应关键点的尺度因子、旋转角和平移位移，进行两种数据的坐标转换。
进一步地，步骤S44中，所述用户定义的值τ取0.3。
进一步地，步骤S5具体包括以下步骤：
步骤S51：利用无人机影像点云和地面激光雷达点云的融合数据生成的DSM，与样地的数字高程模型DEM做差值运算获得冠层高度模型CHM；
步骤S52：基于融合数据的树冠高度模型CHM，将获取的初始样地单木树顶点结果作为候选树顶点；
步骤S53：采用自适应的动态窗口对获取的树顶点进行判断，若当前顶点为对应窗口区域的最大值则保存，否则删除。
进一步地，步骤S52中，所述动态窗口的大小通过计算潜在顶点八个剖面方向半方差值的变化值来自适应的确定，其公式计算如下：
Figure BDA0002034760890000061
式中，γ(h)为经验半方差值，xi为影像的像元位置，h为两个像元的空间分割距离，Z(xi)为对应影像xi处的像元值，N为在一定分隔距离下像元对的对数。
与现有技术相比，本发明有以下有益效果：本发明有效解决了单一数据在高郁闭度森林样地参数提取的弊端，可实现同时从林下和冠层两个方向来测量森林样地参数，为更全面、更精确的估算亚热带高郁闭度森林样地的树高提供了一种新的途径。
附图说明
图1为本发明实施例的方法流程示意图。
图2为本发明实施例的研究样地，编号为0310。
图3为本发明实施例的激光扫描仪及其在样地中的位置图。其中A为RIEGL VZ-400地面激光扫描仪及其扫描视场的示意图，B为各扫描站扫描仪在样地中的位置分布图。
图4为本发明实施例的点云数据。其中A为地面激光雷达(TLS)点云数据，B为无人机(UAV)影像点云数据。
图5为本发明实施例的编号0310样地TLS点云、UAV影像点云生成的数字地面模型(DSM)、两种数据检测到的匹配点对示意图。其中(a)为编号0310样地TLS 点云，(b)为UAV影像点云生成的数字地面模型(DSM)，(c)为两种数据检测到的匹配点对示意图。
图6为本发明实施例的单木点云数据及融合数据。
图7为本发明实施例的四个模型示意图。其中A为基于单一TLS点云数据生成的数字表面模型(DSM)，B为基于单一TLS点云数据生成的数字高程模型(DEM)，C 为基于融合数据获取的数字表面模型(DSM)，D为基于融合数据利用DSM与DEM 的差值生成样地树冠高度模型(CHM)。
图8为本发明实施例的检测样地树顶点及单木树高值，绿色点表示单木树顶点的位置及相应单木的树高。
图9为本发明实施例的估算结果对比图。其中A为基于单一TLS点云估算的树高与野外实测单木树高的对比图，B为基于融合数据估算的树高与野外实测单木树高的对比图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
应该指出，以下详细说明都是示例性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和 /或它们的组合。
如图1所示，本实施例提供了一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，具体包括以下步骤：
步骤S1：同时获取高郁闭度森林样地无人机(UAV)影像点云与该样地的地面激光雷达(TLS)点云数据并生成数字表面模型(DSM)；
步骤S2：估算高郁闭度森林样地无人机(UAV)影像点云与该样地的地面激光雷达(TLS)点云数据的数字表面模型DSM上所有点的曲率，存储DSM所有点云尺度不变的局部邻近点，将曲率值高于其邻近点曲率值的点作为关键点；
步骤S3：建立关键点描述符区域，生成无人机影像点云和地面激光雷达样地点云的DSM描述符SIFT；
步骤S4：匹配最相似的关键点描述符，进行关键点坐标转换，实现无人机影像点云和地面激光雷达样地点云数据匹配融合；
步骤S5：基于无人机影像点云和地面激光雷达样地点云的融合数据，进行样地地面点分割，建立样地冠层高度模型CHM；
步骤S6：基于冠层高度模型CHM，采用改进的局部最大值算法提取样地的单木树高。
在本实施例中，步骤S1中，设置森林样地的郁闭度值大于0.7；进行无人机影像数据采集时，影像航向和旁向重叠度分别设置为大于等于30％和大于等于70％，影像分辨率为设置为小于等于3cm；对得到的UAV像对，利用运动结构多视图立体摄影(SfM-MVS)技术和算法生成三维点云数据。
在本实施例中，步骤S2具体包括以下步骤：
步骤S21：估算两种点云数据的数字表面模型DSM上所有点的表面曲率，利用协方差的特征分析有3个特征向量eμ和其对应的特征值λμ；其中eμ定义一个局部三维正交参考坐标系，而λμ定义每个eμ的大小；根据每点的特征值λ0、λ1、λ2，将局部DSM表面量化为一个近似曲面的曲率，计算公式如下：
Figure BDA0002034760890000091
Figure BDA0002034760890000092
式中，λ0表示候选关键点的表面法线，即是在三维表面上局部切平面的最小变化方向；
步骤S22：步骤S21中各个点的邻近点云半径范围是不确定的，在0-1之间的归一化点云坐标系统中，半径尺度范围是0.1-0.5，增量是0.05，将最大曲率值作为关键点邻近点云的半径值；如果候选点的曲率值大于其相邻点云范围内所有点的曲率值，那么将该候选点作为一个关键点。
在本实施例中，步骤S3具体包括以下步骤：
步骤S31：计算每个关键点的方向，依照这个方向再做进一步的计算；利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性；其中，点(x，y)处梯度的模值和方向计算公式分别如下：
Figure BDA0002034760890000101
Figure BDA0002034760890000102
式中，L所用的尺度为每个关键点各自所在的尺度；
步骤S32：每个关键点有三个信息：位置、所处尺度以及方向，以关键点为中心，确定一个SIFT特征区域；
步骤S33：在关键点邻域内，对其中点的梯度做高斯加权；其中邻域分为4×4 个子区域，每个子区域生成一个描述子，一个描述子中涉及8个方向，最终每个关键点形成128维的SIFT特征向量。
在本实施例中，步骤S4具体包括以下步骤：
步骤S41：采用关键点特征向量的欧式距离作为两幅图像中关键点的相似性判定度量，两个样地的数据分别生成相应的描述子，将两图中所有尺度的描述子进行分配，在描述符特征空间中，利用关键点与其第一和第二最邻近点的欧氏距离来匹配每个关键点；
步骤S42：在描述符空间使用最邻近点特征匹配时，如果两个图像关键点之间的相对距离dist1小于用户设置的阈值，即在边界内搜索不到另一幅图像的最邻近关键点，则关键点的描述符D1不能与另外一幅图像关键点描述符Dtrue匹配；如果边界内只有一个可匹配的点，该点即可作为最佳匹配的点；
步骤S43：在描述符空间使用最邻近点特征匹配时，如果在一定的边界范围内，
Figure BDA0002034760890000111
比较
Figure BDA0002034760890000112
更接近D1，使用最邻近距离
Figure BDA0002034760890000113
的
Figure BDA0002034760890000114
作为最佳匹配，
Figure BDA0002034760890000115
是第二邻近匹配；
步骤S44：在描述符空间使用最邻近点特征匹配时，如果在一定的边界范围内，
Figure BDA0002034760890000116
和
Figure BDA0002034760890000117
的值接近到使用步骤S43无法有效区分时，则使用最邻近距离比法，只有最邻近距离比NNDR小于用户定义的值τ，匹配点才会被接受，其中NNDR的计算公式如下所示：
Figure BDA0002034760890000118
步骤S45：在确认匹配点以后，使用RANSAC方法剔除错误匹配的点；
步骤S46：估算无人机影像点云和地面激光雷达样地点云表面对应关键点的尺度因子、旋转角和平移位移，进行两种数据的坐标转换。
在本实施例中，步骤S44中，所述用户定义的值τ取0.3。
在本实施例中，步骤S5具体包括以下步骤：
步骤S51：利用无人机影像点云和地面激光雷达点云的融合数据生成的DSM，与样地的数字高程模型DEM做差值运算获得冠层高度模型CHM；
步骤S52：基于融合数据的树冠高度模型CHM，将获取的初始样地单木树顶点结果作为候选树顶点；
步骤S53：采用自适应的动态窗口对获取的树顶点进行判断，若当前顶点为对应窗口区域的最大值则保存，否则删除。
在本实施例中，步骤S52中，所述动态窗口的大小通过计算潜在顶点八个剖面方向半方差值的变化值来自适应的确定，其公式计算如下：
Figure BDA0002034760890000121
式中，γ(h)为经验半方差值，xi为影像的像元位置，h为两个像元的空间分割距离，Z(xi)为对应影像xi处的像元值，N为在一定分隔距离下像元对的对数。
特别的，本实施例选择的高郁闭度森林样地位于福建省将乐县，样地为1亩、边长为25.82×25.82m的正方形固定森林样地，样地编号为0310。在固定样地中分别获取UAV影像和TLS点云数据。图2为研究区域及0310样地侧视图和俯视图的实景照片。由图可看出，该样地样木较为密集，郁闭度较高。
在2016年11月，运用RIEGL VZ-400地面激光扫描仪采集样地的TLS点云数据，其扫描仪的扫描视场是100°×360°(垂直×水平)，图3中的A是RIEGL VZ-400地面激光扫描仪及其扫描视场的示意图。图3中的B中的三角形为各扫描站扫描仪在样地中的位置，共8站数据，均匀分布在样地中心和四周。在2016年 11月，完成了0310样地的UAV影像数据的采集，影像航向和旁向重叠度分别设置为30％和70％，影像分辨率为2.5cm。
图4中的A和B分别是编号0310样地获取的TLS点云和UAV影像点云数据。由图可看出，该样地的两种点云数据都在不同的角度存在缺失，TLS点云其地面数据可以精细获取，但是部分单木树冠点云未扫描完整，UAV影像点云是缺失林冠下层信息，为获取精确的样地树高，通过以上步骤，将两种点云数据进行匹配(图 5)，得到融合后的点云数据，获取更加完整的描述该森林样地(图6)。
图7中的A、B是样地基于单一TLS点云数据生成的数字表面模型(DSM)和数字高程模型(DEM)，图7中的C基于融合数据获取的数字表面模型(DSM)，由图可看出，融合前后生成的DSM差别明显，弥补了单一TLS点云数据树木冠层点云数据的缺失。图7中的D基于融合数据利用DSM与DEM的差值生成样地树冠高度模型(CHM)。
通过以上步骤检测样地树顶点及单木树高值如图8所示，绿色点表示单木树顶点的位置，并标注了其高度值来展示相应单木的树高。
图9中的A为本实施案例中基于单一TLS点云估算的树高与野外实测单木树高的对比图，图9中的B为本实施案例中基于融合数据估算的树高与野外实测单木树高的对比图，对比发现基于融合数据估算的树高相关性明显增强。
下表是基于融合数据估算样地的平均树高、单一TLS数据估算的样地平均树高与野外实测平均树高的比较。
Figure BDA0002034760890000131
从上表可以得出，基于融合数据估算的样地平均树高精度明显提升，其估算误差可以降到10％，满足了林业参数调查的精度要求。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/ 或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容,依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。

Claims (8)

1.一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：包括以下步骤：
步骤S1：同时获取高郁闭度森林样地无人机影像点云与该样地的地面激光雷达点云数据并生成数字表面模型；
步骤S2：估算高郁闭度森林样地无人机影像点云与该样地的地面激光雷达点云数据的数字表面模型DSM上所有点的曲率，存储DSM所有点云尺度不变的局部邻近点，将曲率值高于其邻近点曲率值的点作为关键点；
步骤S3：建立关键点描述符区域，生成无人机影像点云和地面激光雷达样地点云的DSM描述符SIFT；
步骤S4：匹配最相似的关键点描述符，进行关键点坐标转换，实现无人机影像点云和地面激光雷达样地点云数据匹配融合；
步骤S5：基于无人机影像点云和地面激光雷达样地点云的融合数据，进行样地地面点分割，建立样地冠层高度模型CHM；
步骤S6：基于冠层高度模型CHM，采用改进的局部最大值算法提取样地的单木树高。
2.根据权利要求1所述的一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：步骤S1中，设置森林样地的郁闭度值大于0.7；进行无人机影像数据采集时，影像航向和旁向重叠度分别设置为大于等于30％和大于等于70％，影像分辨率设置为小于等于3cm；对得到的UAV像对，利用运动结构多视图立体摄影技术和算法生成三维点云数据。
3.根据权利要求1所述的一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：步骤S2具体包括以下步骤：
步骤S21：估算两种点云数据的数字表面模型DSM上所有点的表面曲率，利用协方差的特征分析有3个特征向量eμ和其对应的特征值λμ；其中eμ定义一个局部三维正交参考坐标系，而λμ定义每个eμ的大小；根据每点的特征值λ0、λ1、λ2，将局部DSM表面量化为一个近似曲面的曲率，计算公式如下：
Figure FDA0003548053080000021
Figure FDA0003548053080000022
式中，λ0表示候选关键点的表面法线，即是在三维表面上局部切平面的最小变化方向；
步骤S22：步骤S21中各个点的邻近点云半径范围是不确定的，在0-1之间的归一化点云坐标系统中，半径尺度范围是0.1-0.5，增量是0.05，将最大曲率值作为关键点邻近点云的半径值；如果候选点的曲率值大于其相邻点云范围内所有点的曲率值，那么将该候选点作为一个关键点。
4.根据权利要求1所述的一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：步骤S3具体包括以下步骤：
步骤S31：计算每个关键点的方向，依照这个方向再做进一步的计算；利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性；其中，点(x，y)处梯度的模值和方向计算公式分别如下：
Figure FDA0003548053080000031
Figure FDA0003548053080000032
式中，L所用的尺度为每个关键点各自所在的尺度；
步骤S32：每个关键点有三个信息：位置、所处尺度以及方向，以关键点为中心，确定一个SIFT特征区域；
步骤S33：在关键点邻域内，对其中点的梯度做高斯加权；其中邻域分为4×4个子区域，每个子区域生成一个描述子，一个描述子中涉及8个方向，最终每个关键点形成128维的SIFT特征向量。
5.根据权利要求1所述的一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：步骤S4具体包括以下步骤：
步骤S41：采用关键点特征向量的欧式距离作为两幅图像中关键点的相似性判定度量，两个样地的数据分别生成相应的描述子，将两图中所有尺度的描述子进行分配，在描述符特征空间中，利用关键点与其第一和第二最邻近点的欧氏距离来匹配每个关键点；
步骤S42：在描述符空间使用最邻近点特征匹配时，如果两个图像关键点之间的相对距离dist1小于用户设置的阈值，即在边界内搜索不到另一幅图像的最邻近关键点，则关键点的描述符D1不能与另外一幅图像关键点描述符Dtrue匹配；如果边界内只有一个可匹配的点，该点即可作为最佳匹配的点；
步骤S43：在描述符空间使用最邻近点特征匹配时，如果在一定的边界范围内，
Figure FDA0003548053080000041
比较
Figure FDA0003548053080000042
更接近D1，使用最邻近距离
Figure FDA0003548053080000043
的
Figure FDA0003548053080000044
作为最佳匹配，
Figure FDA0003548053080000045
是第二邻近匹配；
步骤S44：在描述符空间使用最邻近点特征匹配时，如果在一定的边界范围内，
Figure FDA0003548053080000046
和
Figure FDA0003548053080000047
的值接近到使用步骤S43无法有效区分时，则使用最邻近距离比法，只有最邻近距离比NNDR小于用户定义的值τ，匹配点才会被接受，其中NNDR的计算公式如下所示：
Figure FDA0003548053080000048
步骤S45：在确认匹配点以后，使用RANSAC方法剔除错误匹配的点；
步骤S46：估算无人机影像点云和地面激光雷达样地点云表面对应关键点的尺度因子、旋转角和平移位移，进行两种数据的坐标转换。
6.根据权利要求5所述的一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：步骤S44中，所述用户定义的值τ取0.3。
7.根据权利要求1所述的一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：步骤S5具体包括以下步骤：
步骤S51：利用无人机影像点云和地面激光雷达点云的融合数据生成的DSM，与样地的数字高程模型DEM做差值运算获得冠层高度模型CHM；
步骤S52：基于融合数据的冠层高度模型CHM，将获取的初始样地单木树顶点结果作为候选树顶点；
步骤S53：采用自适应的动态窗口对获取的树顶点进行判断，若当前顶点为对应窗口区域的最大值则保存，否则删除。
8.根据权利要求7所述的一种结合UAV影像和TLS点云的高郁闭森林样地树高提取方法，其特征在于：步骤S52中，所述动态窗口的大小通过计算潜在顶点八个剖面方向半方差值的变化值来自适应的确定，其公式计算如下：
Figure FDA0003548053080000051
式中，γ(h)为经验半方差值，xi为影像的像元位置，h为两个像元的空间分割距离，Z(xi)为对应影像xi处的像元值，N为在一定分隔距离下像元对的对数。
CN201910321231.5A 2019-04-22 2019-04-22 一种结合uav影像和tls点云的高郁闭森林样地树高提取方法 Expired - Fee Related CN110033484B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910321231.5A CN110033484B (zh) 2019-04-22 2019-04-22 一种结合uav影像和tls点云的高郁闭森林样地树高提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910321231.5A CN110033484B (zh) 2019-04-22 2019-04-22 一种结合uav影像和tls点云的高郁闭森林样地树高提取方法

Publications (2)

Publication Number Publication Date
CN110033484A CN110033484A (zh) 2019-07-19
CN110033484B true CN110033484B (zh) 2022-05-24

Family

ID=67239552

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910321231.5A Expired - Fee Related CN110033484B (zh) 2019-04-22 2019-04-22 一种结合uav影像和tls点云的高郁闭森林样地树高提取方法

Country Status (1)

Country Link
CN (1) CN110033484B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110580468B (zh) * 2019-09-10 2023-07-18 南京林业大学 一种基于影像匹配点云的单木结构参数提取的方法
CN113205548A (zh) * 2021-04-01 2021-08-03 广西壮族自治区自然资源遥感院 一种林区无人机和地基点云的自动化配准方法及系统
CN113793367B (zh) * 2021-11-17 2022-02-15 西南交通大学 一种工程结构转角位移视觉测量及动力识别系统及方法
CN114387288A (zh) * 2022-01-17 2022-04-22 滁州学院 基于车载激光雷达点云数据的单立木三维信息提取方法
CN114969166B (zh) * 2022-03-01 2024-04-26 国家林业和草原局调查规划设计院 对航空激光雷达估算山地林区chm树高精度的分组评价方法
CN114281871B (zh) * 2022-03-01 2022-10-18 国家林业和草原局调查规划设计院 航空激光雷达估算山地林区chm树高的精度验证方法
CN114777703B (zh) * 2022-04-25 2024-04-16 贵州省第三测绘院(贵州省国土资源遥感监测中心) 基于距离匹配的林业样地精准定位方法及装置
CN115100250A (zh) * 2022-06-13 2022-09-23 北京林业大学 一种林木树高参数提取方法
CN116665081B (zh) * 2023-07-31 2023-10-13 自然资源部第一海洋研究所 一种滨海植被地上生物量的估算方法、计算机设备及介质
CN118015685B (zh) * 2024-04-09 2024-07-02 湖北楚天龙实业有限公司 一种一卡通的识别方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7212670B1 (en) * 2002-05-03 2007-05-01 Imagetree Corp. Method of feature identification and analysis
CN106815850A (zh) * 2017-01-22 2017-06-09 武汉地普三维科技有限公司 基于激光雷达技术获取甚高郁闭度森林蓄积量的方法
CN108594262A (zh) * 2018-07-05 2018-09-28 南京林业大学 一种结合无人机和地基激光雷达提取树高的方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7212670B1 (en) * 2002-05-03 2007-05-01 Imagetree Corp. Method of feature identification and analysis
CN106815850A (zh) * 2017-01-22 2017-06-09 武汉地普三维科技有限公司 基于激光雷达技术获取甚高郁闭度森林蓄积量的方法
CN108594262A (zh) * 2018-07-05 2018-09-28 南京林业大学 一种结合无人机和地基激光雷达提取树高的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于无人机影像匹配点云的苗圃单木冠层三维分割;陈崇成等;《农业机械学报》;20180225(第02期);全文 *

Also Published As

Publication number Publication date
CN110033484A (zh) 2019-07-19

Similar Documents

Publication Publication Date Title
CN110033484B (zh) 一种结合uav影像和tls点云的高郁闭森林样地树高提取方法
US11244197B2 (en) Fast and robust multimodal remote sensing image matching method and system
Fan et al. Pothole detection based on disparity transformation and road surface modeling
CN109903313B (zh) 一种基于目标三维模型的实时位姿跟踪方法
CN103218787B (zh) 多源异构遥感影像控制点自动采集方法
JP6899189B2 (ja) ビジョンシステムで画像内のプローブを効率的に採点するためのシステム及び方法
Cheng et al. Building boundary extraction from high resolution imagery and lidar data
JP2007090448A (ja) 二次元コード検出装置及びそのプログラム、並びに、ロボット制御情報生成装置及びロボット
CN103606170B (zh) 基于彩色尺度不变的街景影像特征检测与匹配方法
CN104574401A (zh) 一种基于平行线匹配的图像配准方法
CN107909018B (zh) 一种稳健的多模态遥感影像匹配方法和系统
CN112164117A (zh) 一种基于Kinect相机的V-SLAM位姿估算方法
Yuan et al. Combining maps and street level images for building height and facade estimation
KR101460313B1 (ko) 시각 특징과 기하 정보를 이용한 로봇의 위치 추정 장치 및 방법
CN109671109B (zh) 密集点云生成方法及系统
CN112102375A (zh) 一种点云配准可靠性检测的方法、装置、移动智慧设备
Wang Automatic extraction of building outline from high resolution aerial imagery
CN113989604A (zh) 基于端到端深度学习的轮胎dot信息识别方法
Xia et al. A table method for coded target decoding with application to 3-D reconstruction of soil specimens during triaxial testing
Junior et al. A new variant of the ICP algorithm for pairwise 3D point cloud registration
Gao et al. Synergizing low rank representation and deep learning for automatic pavement crack detection
TW201601119A (zh) 物件辨識與定位方法
Polat et al. An investigation of tree extraction from UAV-based photogrammetric dense point cloud
CN113724277A (zh) 一种基于Radon变换的电力线检测方法
CN110738098A (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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20220524