一种林区点云快速配准方法
阅读说明：本技术 一种林区点云快速配准方法 (Forest region point cloud rapid registration method ) 是由 王迪 全英汇 徐楷杰 别博文 肖国尧 于 2021-07-26 设计创作，主要内容包括：本发明涉及一种林区点云快速配准方法,包括：步骤1：获取多视角下的每个单站的原始点云数据；步骤2：根据原始点云数据,提取对应的特征点；步骤3：根据特征点,构建Delaunay三角网,计算得到特征点的特征描述向量；步骤4：根据特征描述向量,计算得到不同单站之间对应的相似度矩阵；步骤5：根据相似度矩阵,采用线性分配算法,构建分配总成本TC,使分配总成本TC最小以获得不同单站之间的最终匹配特征点对,根据最终匹配特征点对得到粗配准结果；步骤6：对粗配准结果进行精配准,得到精配准结果。本发明的林区点云快速配准方法,避免了穷竭搜索,极大地提升了计算效率和林区点云配准的效率。(The invention relates to a forest region point cloud rapid registration method, which comprises the following steps: step 1: acquiring original point cloud data of each single station under multiple viewing angles; step 2: extracting corresponding characteristic points according to the original point cloud data; and step 3: constructing a Delaunay triangulation network according to the feature points, and calculating to obtain feature description vectors of the feature points; and 4, step 4: calculating to obtain corresponding similarity matrixes among different single stations according to the feature description vectors; and 5: according to the similarity matrix, a linear distribution algorithm is adopted to construct a distribution total cost TC, the distribution total cost TC is minimized to obtain a final matching characteristic point pair among different single stations, and a coarse registration result is obtained according to the final matching characteristic point pair; step 6: and carrying out fine registration on the coarse registration result to obtain a fine registration result. The forest region point cloud rapid registration method avoids exhaustive search, and greatly improves the calculation efficiency and the forest region point cloud registration efficiency.)
技术领域
本发明属于森林结构测量领域，具体涉及一种林区点云快速配准方法。
背景技术
激光雷达(Light Detection And Ranging，LiDAR)是一种主动遥感技术，通过测定激光信号的传播时间、能量大小、波谱幅度、相位信息变化等获取目标物体的三维几何信息以及辐射特性，该激光信号由传感器发出经目标物反射后再由传感器接收。激光雷达获取的数据是三维空间离散的点集，称为点云，其在刻画研究对象的三维结构信息方面具有独特的优势。
森林在陆地生态系统和人类生产生活中扮演着重要的角色。精准森林结构测量对森林经营管理，森林碳汇估测及其动态变化监测，乃至全球气候变化研究具有重要意义。用来定量描述森林结构的常用参数包括树木位置，胸径，高度，树种等。激光雷达点云因其独特的三维结构测量能力，能够快速、精准、非破坏性地对以上参数进行估测，因此近年来在森林结构测量领域有着广泛应用。由于激光对物体的不可穿透性以及森林环境复杂的遮挡情况，从单个视角获取的激光雷达点云数据无法满足精准森林结构参数测量的要求，从多个视角获取的激光雷达的点云数据需要进行配准从而获得同一坐标框架下的完整信息。
点云配准一般分为粗配准和精配准两个部分，粗配准是将距离远旋转大的点云进行初步的对齐，精配准是对初步对齐的点云进行精确对齐。目前林区点云的配准研究主要集中在粗配准部分。森林点云具有结构复杂度高，树木自相似性，缺乏明显几何特征点等特点，常用的点云特征点如轮廓点、直线点、高曲率点等无法直接应用。目前最常见的林区点云粗配准方法是利用树木的位置信息作为特征点，然后利用随机抽样一致性(RandomConsensus Sampling RANSAC)等计算不同点云的树木位置之间的变换矩阵。但是此方法存在计算复杂度高的问题，而且，三维空间变换的求解需要至少三个同名点，该方法的计算复杂度随树木数量三次方增长。
发明内容
为了解决现有技术中存在的上述问题，本发明提供了一种林区点云快速配准方法。本发明要解决的技术问题通过以下技术方案实现：
本发明提供了一种林区点云快速配准方法，包括：
步骤1：获取多视角下的每个单站的原始点云数据；
步骤2：根据所述原始点云数据，提取对应的特征点；
步骤3：根据所述特征点，构建Delaunay三角网，计算得到所述特征点的特征描述向量；
步骤4：根据所述特征描述向量，计算得到不同单站之间对应的相似度矩阵；
步骤5：根据所述相似度矩阵，采用线性分配算法，构建分配总成本TC，使分配总成本TC最小以获得不同单站之间的最终匹配特征点对，根据所述最终匹配特征点对得到粗配准结果；
步骤6：对所述粗配准结果进行精配准，得到精配准结果。
在本发明的一个实施例中，在所述步骤2中，所述特征点为树木位置。
在本发明的一个实施例中，所述步骤3包括：
步骤31：根据所述特征点，构建Delaunay三角网；
步骤32：根据所述Delaunay三角网，计算得到每个所述特征点与其一级相邻特征点之间的距离以及与其二级相邻特征点之间的距离；
步骤33：所述特征点与其一级相邻特征点之间的距离以及与其二级相邻特征点之间的距离，构成所述特征点对应的特征描述向量F。
在本发明的一个实施例中，所述一级相邻特征点为所述Delaunay三角网中与所述特征点直接连接的特征点；
所述二级相邻特征点为所述Delaunay三角网中与所述特征点的一级相邻特征点直接连接的特征点。
在本发明的一个实施例中，所述步骤4包括：
根据所述特征描述向量，计算不同两个单站中每个特征点的相似度SMi,j，根据所述相似度SMi,j，构建得到所述相似度矩阵，
SMi,j＝Fi∩Fj/min(‖Fi‖,‖Fj‖)，
其中，Fi表示第一单站的原始点云数据S1中的第i个特征点对应的特征描述向量F，Fj表示第二单站的原始点云数据S2中的第j个特征点对应的特征描述向量F，‖Fi‖表示特征描述向量Fi的长度，‖Fj‖表示特征描述向量Fj的长度。
在本发明的一个实施例中，所述步骤5包括：
步骤51：根据所述相似度矩阵，得到成本矩阵Cost，
其中，所述成本矩阵Cost中第i行第j列的值Costi,j表示将第i行分配给第j列的成本，Costi,j＝1-SMi,j；
步骤52：根据所述成本矩阵Cost，构建分配总成本TC，
其中，m表示成本矩阵的总行数，n表示成本矩阵的总列数，p表示已经匹配的行和列的总数，m+n-2p表示未能匹配的行和列的总数，M(i,1)和M(i,2)分别表示一组匹配的行和列；
步骤5.3：求解使得所述分配总成本TC最小，以获得对应的若干匹配特征点对；
步骤5.4：利用随机抽样最大似然估计算法对若干所述匹配特征点对中的错误匹配点对进行去除，得到最终匹配特征点对；
步骤5.5：根据所述最终匹配特征点对，得到空间变换矩阵Q，根据所述空间变换矩阵Q，进行第一次配准，得到粗配准结果。
在本发明的一个实施例中，所述步骤6包括：利用最近点迭代算法对所述粗配准结果进行第二次配准，得到精配准结果。
与现有技术相比，本发明的有益效果在于：
本发明的林区点云快速配准方法，通过计算得到不同单站之间对应的相似度矩阵，根据该相似度矩阵，采用线性分配算法，构建分配总成本TC，利用求解线性分配问题获得极大子集匹配，避免了穷举法测试所有可能的匹配点，将计算复杂度降为线性，避免了穷竭搜索，极大地提升了计算效率和林区点云配准的效率。
上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，而可依照说明书的内容予以实施，并且为了让本发明的上述和其他目的、特征和优点能够更明显易懂，以下特举较佳实施例，并配合附图，详细说明如下。
附图说明
图1是本发明实施例提供的一种林区点云快速配准方法的示意图；
图2是本发明实施例提供的一种林区点云快速配准方法的流程框图；
图3是本发明实施例提供的一种Delaunay三角网以及一级、二级相邻特征点的示意图；
图4是本发明实施例提供的利用线性分配算法求解得到的匹配特征点对；
图5是本发明实施例提供的利用随机抽样最大似然估计算法优化后的最终匹配特征点对；
图6是本发明实施例提供的粗配准结果示意图。
具体实施方式
为了进一步阐述本发明为达成预定发明目的所采取的技术手段及功效，以下结合附图及具体实施方式，对依据本发明提出的一种林区点云快速配准方法进行详细说明。
有关本发明的前述及其他技术内容、特点及功效，在以下配合附图的具体实施方式详细说明中即可清楚地呈现。通过具体实施方式的说明，可对本发明为达成预定目的所采取的技术手段及功效进行更加深入且具体地了解，然而所附附图仅是提供参考与说明之用，并非用来对本发明的技术方案加以限制。
实施例一
请结合参见图1和图2，图1是本发明实施例提供的一种林区点云快速配准方法的示意图；图2是本发明实施例提供的一种林区点云快速配准方法的流程框图。如图所示，本实施例的林区点云快速配准方法，包括：
步骤1：获取多视角下的每个单站的原始点云数据；
由于激光对物体的不可穿透性以及森林环境复杂的遮挡情况，从单个视角获取的激光雷达的点云数据无法满足精准森林结构参数测量的要求，因此，需要从多个视角获取激光雷达的点云数据，然后对其进行配准从而获得同一坐标框架下的完整信息。
在本实施例中，每个单站的原始点云数据也就是从单个视角获取的激光雷达的点云数据。
步骤2：根据原始点云数据，提取对应的特征点；
在本实施例中，特征点为树木位置，也就是将林区中的每个单木位置作为原始点云数据中的特征点。
需要说明的是，可选地，树木位置可由聚类、圆检测等方法从原始点云数据提取得到，在本实施例中，对提取树木位置的具体算法不做任何限制，任何单木位置感知算法均可以。
步骤3：根据特征点，构建Delaunay三角网，计算得到特征点的特征描述向量；
在本实施例中，为避免穷竭搜索，首先计算特征点的特征，描述特征点的空间属性，采用特征点之间的相对空间分布作为特征描述。
具体地，步骤3包括：
步骤31：根据特征点，构建Delaunay三角网；
步骤32：根据Delaunay三角网，计算得到每个特征点与其一级相邻特征点之间的距离以及与其二级相邻特征点之间的距离；
请参见图3，图3是本发明实施例提供的一种Delaunay三角网以及一级、二级相邻特征点的示意图。如图所示，在本实施例中，特征点的一级相邻特征点为Delaunay三角网中与该特征点直接连接的特征点；特征点的二级相邻特征点为Delaunay三角网中与该特征点的一级相邻特征点直接连接的特征点。
步骤33：特征点与其一级相邻特征点之间的距离以及与其二级相邻特征点之间的距离，构成特征点对应的特征描述向量F。
在本实施例中，特征描述向量F为一维长度为n的向量，该特征描述向量F表示某个树木位置(即特征点)距离其一级近邻树木和二级近邻树木的平面位置，其中，该树木位置的一级近邻树木和二级近邻树木的总数为n个。
步骤4：根据特征描述向量，计算得到不同单站之间对应的相似度矩阵；
以第一单站的原始点云数据S1中的每个特征点与第二单站的原始点云数据S2中每个特征点的相似度计算为例，对其进行具体说明：
具体地，步骤4包括：
根据特征描述向量，计算不同两个单站中每个特征点的相似度SMi,j，根据相似度SMi,j，构建得到相似度矩阵，其中，
SMi,j＝Fi∩Fj/min(‖Fi‖,‖Fj‖) (1)，
其中，Fi表示第一单站的原始点云数据S1中的第i个特征点对应的特征描述向量F，Fj表示第二单站的原始点云数据S2中的第j个特征点对应的特征描述向量F，‖Fi‖表示特征描述向量Fi的长度，‖Fj‖表示特征描述向量Fj的长度。
在本实施例中，根据相似度SMi,j，构建得到相似度矩阵，该相似度矩阵的大小为Ns1×Ns2，其中，其中Ns1表示第一单站的原始点云数据S1中特征点的个数，Ns2表示第二单站的原始点云数据S2中特征点的个数。
需要说明的是，在实际计算中，由于特征点(即树木位置)的计算存在误差，因此，在计算特征描述向量相交(即，计算Fi∩Fj)时，允许存在偏差eF，在本实施例中，偏差eF设为0.2米，即假设树木位置估测的误差大小为0.2米。
步骤5：根据相似度矩阵，采用线性分配算法，构建分配总成本TC，使分配总成本TC最小以获得不同单站之间的最终匹配特征点对，根据最终匹配特征点对得到粗配准结果；
其中，线性分配算法是一种将矩阵行分配给列的方法，要求每行都分配给一列，并且分配的总成本最小化或最大化。
具体地，步骤5包括：
步骤51：根据相似度矩阵，得到成本矩阵Cost，
其中，成本矩阵中第i行第j列的值Costi,j表示将第i行分配给第j列的成本，Costi,j＝1-SMi,j；
在本实施例中，成本矩阵Cost的代价即为非相似度，对于每个未匹配的行或列分配一个未分配成本，例如，某个行和列未匹配，则会使总成本增加2CostUnmatched。
步骤52：根据成本矩阵Cost，构建分配总成本TC，
其中，m表示成本矩阵的总行数，n表示成本矩阵的总列数，p表示已经匹配的行和列的总数，m+n-2p表示未能匹配的行和列的总数，M(i,1)和M(i,2)分别表示一组匹配的行和列。
在本实施例中，M的大小为p×2。
在本实施例中，成本矩阵Cost为特征点之间的非相似性，将CostUnmatched设为0.25，表示不同两个单站中的两个特征点之间的非相似度如果大于50％，则不应该匹配，因为非匹配时系统只增加2CostUnmatched的成本，低于二者匹配之后的成本。
步骤5.3：求解使得分配总成本TC最小，以获得对应的若干匹配特征点对；
在本实施例中，通过线性分配算法求解得到一系列的匹配特征点对Mi,j，避免了穷举法测试所有可能的匹配点，将计算复杂度降为线性，极大地提升了计算效率。
步骤5.4：利用随机抽样最大似然估计算法对若干匹配特征点对中的错误匹配点对进行去除，得到最终匹配特征点对；
进一步地，基于线性分配算法求解得到的匹配特征点对Mi,j，利用随机抽样最大似然估计(M-estimator SAmple Consensus，MSAC)算法进一步去除一系列匹配特征点对中错误的匹配特征点对，得到最终匹配特征点对。具体算法与现有MSAC算法类似，在此不再赘述。
随机抽样最大似然估计算法的具体步骤在此不再赘述。
步骤5.5：根据最终匹配特征点对，得到空间变换矩阵Q，根据空间变换矩阵Q，进行第一次配准，得到粗配准结果。
在本实施例中，根据最终匹配特征点对，求取空间变换矩阵Q，以实现平面配准的目的。
其中，空间变换矩阵Q由平移变换矩阵T和旋转矩阵R组成，其矩阵形式为：
其中，tx和tx分别表示X和Y轴方向的平移，θ表示围绕原点的旋转角度。
求解该变换矩阵Q需要至少两对匹配特征点对，MSAC算法每次迭代随机选取两对匹配特征点，求解变换矩阵，并测试变换后两个单站的原始点云数据的匹配程度，即匹配的特征点数量，从而经过多次迭代选取出最优解。
然后，针对Z轴平移变换，分别计算匹配特征点对的各自坐标系的Z值中值，使得中值对齐，得到粗配准结果。
在本实施例中，步骤5在于如何寻找第一单站的原始点云数据S1和第二单站的原始点云数据S2之间的对应点对，从而计算变换关系。求解平面变换需要两对对应点，而三维变换需要三对对应点。传统的方法穷竭测试所有可能的点对，计算量庞大。本实施例在相似矩阵的基础上采用线性分配算法快速求解二维平面变换，再单独求解Z轴平移变换，从而达到粗配准两站三维点云的目的，得到粗配准结果。
步骤6：对粗配准结果进行精配准，得到精配准结果。
具体地，利用最近点迭代算法对粗配准结果进行第二次配准，得到精配准结果。
需要说的是，粗配准后的多站点云仅为初步对齐，需要进一步采取最近点迭代(Iterative Closest Point,ICP)算法进行精配准。
ICP是点云精配准最经典的算法，其基本思想是每一次迭代过程中采用两站数据间的最邻近点对求取空间变换矩阵，经过迭代后两站数据逐步靠近，最终达到收敛。与前述步骤5中的特征点对不同的是，此时，点对仅为空间最近点，因此ICP算法要求两站点云数据已经初步配准，多用于精配准。
与前述步骤5中二维变换矩阵类似，ICP求解三维变换矩阵，包括平移变换：
绕X轴旋转变换：
绕Y轴旋转变换：
绕Z轴旋转变换：
具体算法与现有ICP算法类似，在此不再赘述。
本实施例的林区点云快速配准方法，通过计算得到不同单站之间对应的相似度矩阵，根据该相似度矩阵，采用线性分配算法，构建分配总成本TC，利用求解线性分配问题获得极大子集匹配，避免了穷举法测试所有可能的匹配点，将计算复杂度降为线性，避免了穷竭搜索，极大地提升了计算效率和林区点云配准的效率。
实施例二
本实施例对实施例一的林区点云快速配准方法的配准效果进行了实验验证。
具体地，采用Faro地基激光雷达扫描仪获取的森林点云数据进行测试。请结合参见图4-图6图4是本发明实施例提供的利用线性分配算法求解得到的匹配特征点对；图5是本发明实施例提供的利用随机抽样最大似然估计算法优化后的最终匹配特征点对；图6是本发明实施例提供的粗配准结果示意图(图中实心点表示第一单站的原始点云数据S1中的特征点，空心圈表示第二单站的原始点云数据S2中的特征点)。根据定量测试结果表明，粗配准后旋转误差为0.03°，X轴平移误差为0.34米，Y轴平移误差为0.07米，Z轴平移误差为0.06米；经过ICP精配准之后，旋转误差为0.12°，X轴平移误差为0.02米，Y轴平移误差为0.05米，Z轴平移误差为0.04米。
应当说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的物品或者设备中还存在另外的相同要素。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。