CN103020966B - 一种基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法 - Google Patents

Abstract

一种基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法：首先分别从航空、地面LiDAR数据中提取建筑物轮廓，简称航空轮廓、地面轮廓；再从航空轮廓、地面轮廓中提取出建筑物角点，简称航空角点、地面角点；然后以航空轮廓与地面轮廓间的匹配度为约束，计算航空角点与地面角点之间初始转换矩阵，并获取初始匹配角点对；最后使用ICP算法计算初始匹配角点对之间的修正转换矩阵，并用初始转换矩阵和修正转换矩阵依次对待匹配地面点云数据进行转换，实现航空与地面LiDAR数据的自动高精度配准。本发明使用轮廓线做约束，在配准的可靠性与精确性方面都有很大的优势；同时，本发明仅从待匹配LiDAR数据与基准LiDAR数据出发，无需借助其他辅助数据便可实现两者之间的精确配准。

Description

一种基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法

技术领域

本发明涉及一种航空和地面LiDAR数据配准方法，特别是涉及一种基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法。

背景技术

目前，激光雷达技术(LiDAR)正在蓬勃发展，激光雷达大家庭中有航空LiDAR、地面LiDAR、车载LiDAR、室内LiDAR。不同平台的激光雷达性能各不相同，应用范围互相补充。伴随着激光雷达技术的不断进步，多种平台的激光雷达的融合处理已经逐渐成为一种趋势。目前使用最多的激光雷达是航空LiDAR和地面LiDAR。航空LiDAR具有较大的扫描范围，能够获取物体顶部信息，然而点云条带现象明显，地物侧面信息缺失；地面LiDAR能够获取地物详尽的侧面信息，扫描精度也极高，然而扫描范围有限，顶部信息也难以获得。它们之间各有优缺，两者的集成能够全面地反应地物各个尺度，各个方向的信息。这些年来，两者的集成应用也爆炸式地出现了，最典型的应用如下：1)地质勘探，如地形制图，侵蚀量计算，滑坡、滚石监测；2)森林应用，如森林蓄积量计算，冠层结构调查；3)水文研究，如洪水模型，河流环境变化；4)3D场景构建，如表面模型，城市模型(Bremer，Ruiz，Jaboyedoff，Heckmann，Jung，Lovell，Sampson，Hohenthal，Andrews，Jaw，Fruh等)。

尽管航空和地面LiDAR数据集成应用的研究越来越多，目前对两者的配准研究却很少，很多应用都是通过人工选择控制点进行配准，配准精度较低。而高精度的配准是激光雷达技术集成应用的一个先决条件，因此研究航空和地面激光雷达数据自动配准的方法具有重要的意义。由于航空LiDAR数据和地面LiDAR数据本身的差异，两者的配准难度很大：1)不同视角。航空LiDAR以很小的视角俯视获取数据，顶部信息较多，侧面信息较少，而地面LiDAR平视时或仰视获取数据，侧面信息详尽，顶部信息较少，两者的公共信息较少。2)不同平台。航空平台是移动平台，地面平台是静止平台。3)不同分辨率。航空LiDAR距离扫描目标几百米至上千米不等，数据分辨率在米级或分米级，而地面LiDAR距离扫描目标几十米，数据分辨率最高可达毫米级，理论上准确的一对一配准可能变为一对多配准。4)不同范围。航空LiDAR能够获取大范围的数据，提高大尺度的配准基元，而地面LiDAR扫描范围有限，提供的配准基元较少，位置集中。5)点云数据的离散性。点云数据本身具有离散性，从中获取匹配特征比较困难。

当前航空LiDAR与地面LiDAR配准的方法主要可以分为两类：1)借助第三方数据进行辅助配准；2)单纯使用LiDAR数据进行配准。其中第一种方法的思路是借助GPS、航空影像等其他数据，获取地面扫描仪的位置，以此为参照点实现两种数据的配准(Bohm，Hohethal，Bremer，HeckHarm，Fruh，Zakhor等)，然而这种方法的数据可获取性和数据精度得不到保障，因此实现起来有一定难度。第二种方法的思路是提取出航空LiDAR与地面LiDAR数据中的公共配准基元(包括点基元、线基元和面基元)，通过配准基元之间的匹配来实现数据之间的配准。这种方法不依靠外部数据，单纯使用LiDAR数据进行配准，是自动配准方法发展的方向。然而现阶段对这类方法的研究还不够成熟，在稳定性、计算量和自动化程度等方面还存在一定的问题，如何从航空和地面LiDAR数据中获取准确的配准基元，并利用这些基元实现两种数据的高精度配准仍然有待研究。

发明内容

本发明要解决技术问题是：克服现有技术缺点，提出一种基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，提高航空和地面LiDAR数据的配准可靠性，同时可实现航空和地面LiDAR数据的高精度配准。

为了解决上述技术问题，本发明提出的技术方案是：一种基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，包括以下步骤：

第一步、提取建筑物轮廓——从航空LiDAR数据中提取建筑物轮廓，称为航空轮廓；从地面LiDAR数据中提取建筑物轮廓，称为地面轮廓；

第二步、提取建筑物角点——从航空轮廓中提取建筑物角点，称为航空角点；从地面轮廓中提取建筑物角点，称为地面角点；

第三步、寻找轮廓线段约束下的初始转换矩阵——使用航空角点与地面角点迭代计算转换矩阵，用该转换矩阵对地面轮廓进行转换，并使用航空轮廓与转换后地面轮廓的匹配度作为控制约束条件，当航空轮廓与转换后地面轮廓之间成功匹配的线段对数满足给定阈值时停止迭代，相应的转换矩阵即为初始转换矩阵；

第四步、获得初始匹配角点对——使用第三步获得的初始转换矩阵对地面角点进行转换，根据空间距离寻找航空角点中与其配对的角点，得到初始匹配角点对；

第五步、寻找修正转换矩阵——以初始匹配角点对为源数据，寻找两者间的修正转换矩阵，保证经该修正转换矩阵配准后，两者的均方根误差小于预设的极限值ε，极限值ε的取值范围为0.25-0.35；

第六步、LiDAR数据配准——使用初始转换关系与修正转换关系依次对地面LiDAR数据进行转换，得到最终配准结果。

本发明使用轮廓线做约束，在配准的可靠性与精确性方面都有很大的优势；同时，本发明仅从待匹配LiDAR数据与基准LiDAR数据出发，无需借助其他辅助数据便可实现两者之间的精确配准。

本发明基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，还具有如下改进：

1)、本发明第一步从航空LiDAR数据中提取建筑物轮廓的方法如下：构建1m*1m的水平格网，根据点面空间关系计算每个格网中最高点与最低点的高差，保留高差大于实验区最低建筑物高程的格网得到轮廓格网，对轮廓格网使用多尺度的Hough变换，得到航空LiDAR建筑物轮廓线段。

2)、本发明第一步中从地面LiDAR数据中提取建筑物轮廓的方法如下：使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓；在此基础上使用轮廓延伸密度方法对提取的建筑物轮廓进行恢复，形成完整的建筑物轮廓。

3)、本发明第二步中提取建筑物角点的方法如下：将建筑物轮廓投影到三维坐标系的XY平面内寻找二维相交点，如果任两条构成相交点的轮廓的高程差小于1m，则判定两条轮廓在实际的三维空间中相交，两条轮廓的相交点为一个建筑物角点，并将所述两条轮廓的高程均值作为该建筑物角点的高程。

4)、使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓，具体步骤如下：

1a)地面LiDAR点云投影至XY平面——将地面LiDAR点云投影至三维坐标系的XY平面，并保留各个点的X、Y、Z属性；

1b)提取粗略轮廓格网——在所述XY平面内构建1m*1m的粗略格网，计算每个粗略格网中LiDAR投影点的数量，即得到该粗略格网的格网密度，根据建筑物边缘轮廓处的粗略格网密度阈值对所述粗略格网进行筛选，保留格网密度大于所述粗略格网密度阈值的粗略格网，得到粗略轮廓格网；

1c)提取精确轮廓格网——在粗略轮廓格网中构建0.2m*0.2m的精细格网，计算精细格网内LiDAR投影点的数量即得到精细格网的格网密度，根据建筑物边缘轮廓处的精细格网密度阈值对所述精细格网进行筛选，保留格网密度大于所述精细格网密度阈值的精细格网，得到精确轮廓格网；

1d)格网高差筛选——遍历所有精确轮廓格网，如果精确轮廓格网内的最高LiDAR点和最低LiDAR点的高差大于相应实验区建筑最低高程则保留该精确轮廓格网，否则剔除；

1e)获取轮廓线段——对筛选后的精确轮廓格网使用Hough变换得到二维矢量轮廓线段。

5)、上述步骤1b)和1c)中格网密度阈值的确定方法如下：

假设O点为仪器中心点，A点为水平垂直于仪器的墙面点，扫描仪对准A点时的角度为0°，B点为格网靠近仪器一侧，C点为格网远离仪器一侧，D点为B点竖直方向上墙面最高点，设OA＝D_V，CO＝D_M，水平方向格网的边长为D_G，建筑高H_B，仪器高H_L，在A点处水平向相邻两LiDAR点的间距为D_R，则格网密度计算方法如下：

2a)计算水平方向格网内LiDAR点的列数，记2*α为扫描仪每次旋转角度，记格网中水平方向上最靠近于B点的角度为β， $\mathrm{\β}=(2*[\frac{\arctan ((\sqrt{D_M^2{-D}_V^2}-D_G)/D_V)-\mathrm{\α}}{2*\mathrm{\α}}]+3)*\mathrm{\α},$ 则线段BC上LiDAR点数为则水平方向格网内LiDAR点的列数为N_col；

2b)计算每一列LiDAR点的数量，第i列LiDAR点的数量为

$N_{\mathrm{row}}^i=\left[\frac{\arctan \left(\frac{H_B-H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right)-\mathrm{\α}}{2*\mathrm{\α}}\right]+\left[\frac{\arctan \left(\frac{H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right)-\mathrm{\α}}{2*\mathrm{\α}}\right]+3;$

2c)将每列LiDAR点数量相加得到格网处LiDAR点总数网格密度阈值threshod＝rate*N，其中参数rate为描述墙面凹凸、窗户多少的阈值，墙面凹凸越多、窗户越多，则该参数越小，rate的取值范围为0.2-1。

6)、本发明中，使用轮廓延伸密度的方法进行建筑物轮廓的恢复，具体如下：

3a)寻找步骤1e)中获得的二维矢量轮廓线段周边1m范围内格网，将寻找到的所有格网内LiDAR点最大高程的平均值作为二维矢量轮廓线段的高程，将二维矢量轮廓线段变换为三维建筑物轮廓线段；

3b)对三维建筑物轮廓线段构建半径为1m的缓冲区，建缓冲区内LiDAR点数量除以缓冲区体积获得原有轮廓LiDAR点密度；

3c)沿轮廓线段方向以单位距离为延伸步长构建半径为1m的缓冲区，缓冲区内LiDAR点数量除以相应缓冲区体积获得待延伸方向的LiDAR点密度，所述单位距离的取值范围为0.1-0.3m；

3d)若待延伸方向的LiDAR点密度与原有轮廓LiDAR点密度的差异小于20％，则该轮廓沿轮廓线段方向延伸单位距离并重复步骤3c)；否则停止延伸，形成完整的建筑物轮廓。

7)、本发明所述第三步中寻找初始转换矩阵的具体方法如下：

4a)设航空角点的点集分别为A＝{A_i，i＝0，1，2，...，u}；地面角点的点集为B＝{B_i，i＝0，1，2，...，v}，航空轮廓的线段集为LA＝{LA_i，i＝0，1，2，...，m}；地面轮廓的线段集为LB＝{LB_i，i＝0，1，2，...，n}，u为航空角点的数量，v为地面角点的数量，m为航空轮廓的线段数量，n为地面轮廓的线段数量；

4b)从点集A和B中分别选取1个点A_x和B_x，计算由B_x至A_x的平移矩阵，利用该平移矩阵对点集B中的每个点进行平移，得到点集M＝{M_i，i＝0，1，2，...，v}；

4c)从点集A和M中分别选取1个点A_y和M_y，要求A_y≠A_x，M_y≠B_x，以点A_x为原点，计算点M_y旋转至点A_y位置的旋转矩阵，使用该旋转矩阵对点集M中每个点进行旋转，得到点集R＝{R_i，i＝0，1，2，...，v}；

4d)使用步骤4b)中得到的平移矩阵和步骤4c)中得到的旋转矩阵对地面轮廓的线段集LB进行转换，得到转换后的轮廓线段集LC＝{LC_i，i＝0，1，2，...，n}；

4e)遍历线段集LA中的所有线段，寻找线段集LA和线段集LC之间满足给定匹配条件的线段对数量：

从线段集LA中取出一条未检验的线段，计算其与线段集LC中所有尚未匹配轮廓线段的互异度，若最小互异度小于1，则与该最小互异度对应的两条轮廓线为已匹配线段对，重复本过程直至线段集LA中不存在未检验的线段；两条轮廓线段的互异度计算公式如下：

dif＝w₁×lenDif+w₂×lDis+w₃×lAng+w₄×cpDis

其中，dif为两条轮廓线段的互异度；lenDif为两条轮廓线段的长度之差，lDis为两条轮廓线段所在直线的距离，lAng为两条轮廓线段所在直线的夹角， $\mathrm{lAng}=(\frac{\mathrm{\π}}{2}-|\arccos \left(\frac{{\stackrel{\→}{l}}_1\•{\stackrel{\→}{l}}_2}{\left|{\stackrel{\→}{l}}_1\right|\•\left|{\stackrel{\→}{l}}_2\right|}\right)-\frac{\mathrm{\π}}{2}\left|\right)*\frac{180}{\mathrm{\π}};$ cpDis为两轮廓线中点之间的距离， $\mathrm{cpDis}=\sqrt{{(x_{l_{11}}-x_{l_{21}})}^2+{(y_{l_{11}}-y_{l_{21}})}^2+{(z_{l_{11}}-z_{l_{21}})}^2};$ w1、w2、w3、w4分别表示以上4个参量的权重，取值范围分别为[0.4，0.6]，[0.1，0.3]，[0.1，0.2]，[0.1，0.3]，分别代表两条轮廓线段的方向向量：

${\stackrel{\→}{l}}_1=(x_{l_{12}}-x_{l_{11}},y_{l_{12}}-y_{l_{11}},z_{l_{12}}-z_{l_{11}})$

${\stackrel{\→}{l}}_2=(x_{l_{22}}-x_{l_{21}},y_{l_{22}}-y_{l_{21}},z_{l_{22}}-z_{l_{21}})$

其中分别代表LA中参与匹配的轮廓线段的起点和终点x、y、z坐标，分别代表LC中参与匹配的轮廓线段的起点和终点x、y、z坐标；

4f)若线段集LA和线段集LC中匹配的线段对数量不少于6对，则认为该转换矩阵可靠，该转换矩阵为初始转换矩阵；否则转至步骤4b)重新进行转换矩阵计算。

8)、第四步中，使用第三步获得的初始转换矩阵对地面角点的点集B进行转换得到点集F＝{F_i，i＝0，1，2，...，v}，根据空间距离从航空角点的点集A＝{A_i，i＝0，1，2，...，u}和点集F＝{F_i，i＝0，1，2，...，v}选出成功配对的角点，分别记录为点集P＝{P_i，i＝0，1，2，...，m}和点集U＝{U_i，i＝0，1，2，...，n}，其中m＝n，得到初始匹配角点对；第五步寻找修正转换矩阵的方法如下：

5a)在点集P中找出距点集U中每一个点的最近点，组成点集Q＝{Q_i，i＝0，1，2，...n}；

5b)采用最小均方根法计算点集U与点集Q之间的配准关系，得到配准转换矩阵；

5c)对点集U用配准转换矩阵进行坐标转换，得到点集U₁；

5d)计算点集U₁与点集Q之间的均方根误差，如小于预设的极限值ε，极限值ε取0.3，则以该配准转换矩阵作为修正转换矩阵；否则，以点集U₁替换U，转至步骤5a)重新进行配准转换矩阵计算。

9)、所述第三步获得的初始转换矩阵包括旋转矩阵R、平移矩阵T，第五步获得的修正转换矩阵包括旋转矩阵R′、平移矩阵T′，第六步中，用初始转换矩阵对地面LiDAR数据PB进行初始转换，得到初始配准点云数据PB′，PB′＝{PB′_i＝R×PB_i+T，i＝1，2，3，…，C_B}；用修正转换矩阵初始配准点云数据PB′进行修正转换，得到最终配准点云数据PB″，PB″＝{PB″_i＝R′×PB′_i+T′，i＝1，2，3，…，C_B}，其中C_B为地面LiDAR数据中LiDAR点的数量。

本发明的有益成果是：1)、一般情况下从LiDAR数据中提取出的建筑物角点都会比提取出的轮廓线段少，能够用于计算转换矩阵的公共角点则更少，本发明针对这一特点，使用轮廓作为约束条件，可以有效规避由于角点过少引起的错误匹配，提高点云数据配准的精确度和可靠性；2)、本发明可以自动从航空LiDAR数据与地面LiDAR数据中寻找配准基元——轮廓线和角点，计算两者之间的转换关系，无需其他辅助数据就能实现两种数据的精确配准；3)本发明使用一种分层次格网密度方法提取地面LiDAR数据的建筑物轮廓，并使用理论估计的方法对格网密度阈值进行确定，能够从地面LiDAR数据中提取准确的建筑物轮廓线段，从而提取高精度的地面角点；4)使用了轮廓密度延伸的方法对提取的建筑轮廓进行恢复，能够将不完整的轮廓恢复成较为完整的轮廓，提高了提取轮廓的准确性。

附图说明

下面结合附图对本发明作进一步的说明。

图1为本发明实施例的流程图。

图2-a为本发明实施例中航空LiDAR数据示意图。

图2-b为本发明实施例中地面LiDAR数据示意图。

图3-a从图2-a所示点云数据中提取出的航空轮廓线段示意图。

图3-b从图2-b所示点云数据中提取出并经恢复的地面轮廓线段示意图。

图4-a为从图3-a所示航空轮廓线段数据中提取出的航空角点示意图。

图4-b为从图3-b所示地面轮廓线段数据中提取出的地面角点示意图。

图5-a为实施轮廓线作为约束条件下的建筑物角点匹配后的角点示意图。

图5-b为实施轮廓线作为约束条件下的建筑物角点匹配后的轮廓线示意图。

图6为初始匹配角点对示意图。

图7为发明实施例最终配准后两种数据叠加的示意图。

图8-a为验证可靠性实施例的地面轮廓与地面角点示意图。

图8-b为验证可靠性实施例的航空轮廓与航空角点示意图。

图8-c为验证可靠性实施例中利用本发明方法得到的配准结果示意图。

图8-d为验证可靠性实施例中单独使用角点进行配准的配准结果示意图。

图9为格网密度阈值计算示意图。

具体实施方式

本实施例以航空LiDAR数据PA＝{PA_i，i＝0，1，2，...，C_A}为基准，将LeicaScanStation2扫描得到的地面LiDAR数据PB＝{PB_i，i＝0，1，2，...，C_B}配准到航空数据上。航空LiDAR数据平均点间距1m，高程精度15cm，平面精度30cm，点数约1100万个(见图2-a)；地面LiDAR数据由9站数据通过标靶拼接而成，数据精度20cm/100m，点数约3000万个(见图2-b)。

本发明实施例基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法(流程图见图1)，包括以下步骤：

第一步、提取建筑物轮廓——从航空LiDAR数据中提取建筑物轮廓，称为航空轮廓；从地面LiDAR数据中提取建筑物轮廓，称为地面轮廓。

从航空LiDAR数据中提取建筑物轮廓的方法如下：构建1m*1m的水平格网，根据点面空间关系计算每个格网中最高点与最低点的高差，保留高差大于实验区最低建筑物高程的格网得到轮廓格网，本实施例中最低建筑物高程为10m，保留高差大于10m的格网，对轮廓格网使用多尺度的Hough变换，得到航空LiDAR建筑物轮廓线段LA＝{LA_i，i＝0，1，2，...，m}，如图3-a所示。

从地面LiDAR数据中提取建筑物轮廓的方法如下：使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓；在此基础上使用轮廓延伸密度方法对提取的建筑物轮廓进行恢复，形成完整的建筑物轮廓。

其中，使用分层次的格网密度方法从地面Li DAR数据中提取建筑物轮廓，具体步骤如下：

1a)地面LiDAR点云投影至XY平面——将地面LiDAR点云投影至三维坐标系的XY平面，并保留各个点的X、Y、Z属性；

1b)提取粗略轮廓格网——在所述XY平面内构建1m*1m的粗略格网，计算每个粗略格网中LiDAR投影点的数量，即得到该粗略格网的格网密度，根据建筑物边缘轮廓处的粗略格网密度阈值对所述粗略格网进行筛选，保留格网密度大于所述粗略格网密度阈值的粗略格网，得到粗略轮廓格网；

本实施例中，使用理论估计的方法计算墙面筛选密度阈值，本例中，扫描的最小楼高为20m，最小水平距离为8m，仪器高为1.5m，建筑距离测站最远距离约32m，扫描精度为100m远处20cm一个点，考虑到墙面窗户较多，设置参数rate为0.5，对于1m*1m的格网，经过理论估计后得到最终阈值为1100，即格网密度大于1100的所有格网都为粗略轮廓格网。

1c)提取精确轮廓格网——在粗略轮廓格网中构建0.2m*0.2m的精细格网，计算精细格网内LiDAR投影点的数量即得到精细格网的格网密度，根据建筑物边缘轮廓处的精细格网密度阈值对所述精细格网进行筛选，保留格网密度大于所述精细格网密度阈值的精细格网，得到精确轮廓格网；

本实施例在提取得到的1m*1m的轮廓格网内，构建0.2m*0.2m的精细格网，使用理论估计的方法计算得到筛选阈值为550。

1d)格网高差筛选——遍历所有精确轮廓格网，如果精确轮廓格网内的最高LiDAR点和最低LiDAR点的高差大于相应实验区建筑最低高程则保留该精确轮廓格网，否则剔除；本例中实验区最低建筑物高程为10m。

1e)获取轮廓线段——对筛选后的精确轮廓格网使用Hough变换得到二维矢量轮廓线段LB＝{LB_i，i＝0，1，2，...，n}。考虑到大尺度的Hough变换有助于获取比较完整的线段；而小尺度的Hough变换有助于获取比较零碎的线段；因此本实施例分两个尺度对轮廓区域进行Hough变换，首先对完整的精确轮廓格网进行Hough变换，然后将精确轮廓格网分为16个小块分别进行Hough变换，最后将各个结果拼接融合。经过该这样的变换处理后，轮廓提取效果更好。

本实施例在上述步骤1b)和1c)中格网密度阈值的确定方法如下：

如图9所示，假设O点为仪器中心点，A点为水平垂直于仪器的墙面点，扫描仪对准A点时的角度为0°，B点为格网靠近仪器一侧，C点为格网远离仪器一侧，D点为B点竖直方向上墙面最高点，墙面上的圆点为仪器扫描获得的LiDAR点，从图中可见，LiDAR点在墙面上呈现阵列式分布，由于扫描仪每次旋转的角度是固定的，因此离扫描仪越近的墙面上LiDAR点分布越密，相反的，离扫描仪越远的墙面上LiDAR点分布越疏，设OA＝D_V，CO＝D_M，水平方向格网的边长为D_G，建筑高H_B，仪器高H_L，在A点处水平向相邻两LiDAR点的间距为D_R，则格网密度阈值的具体计算方法如下：

2a)计算水平方向格网内LiDAR点的列数，记2*α为扫描仪每次旋转角度，记格网中水平方向上最靠近于B点的角度为β， $\mathrm{\β}=(2*[\frac{\arctan ((\sqrt{D_M^2{-D}_V^2}-D_G)/D_V)-\mathrm{\α}}{2*\mathrm{\α}}]+3)*\mathrm{\α},$ 则线段BC上LiDAR点数为则水平方向格网内LiDAR点的列数为N_col；

2b)计算每一列LiDAR点的数量，第i列LiDAR点的数量为

$N_{\mathrm{row}}^i=\left[\frac{\arctan \left(\frac{H_B-H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right)-\mathrm{\α}}{2*\mathrm{\α}}\right]+\left[\frac{\arctan \left(\frac{H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right)-\mathrm{\α}}{2*\mathrm{\α}}\right]+3;$

2c)将每列LiDAR点数量相加得到格网处LiDAR点总数网格密度阈值threshod＝rate*N，其中参数rate的取值范围为0.2-1，该参数取值与墙面凹凸情况与窗户多少有关，墙面凹凸越多、窗户越多，取值越小，墙面平滑、窗户越少取值越大，当墙面平滑不含有窗户时，参数rate取1，当墙面凹凸起伏含有极大量窗户时rate取0.2，在本例中，rate取0.5。

上述步骤的格网密度阈值推导过程如下：

如图9所示，A点为水平垂直于仪器的墙面点，扫描仪对准A点时的角度为0°；位置E点为格网外水平方向上最靠近B点的扫描点(即E点后面的一个扫描点落入格网范围之内)。

那么，其中 $\∠\mathrm{AOB}=\arctan ((\sqrt{D_M^2-D_V^2}-D_G)/D_V).$

F点为格网中水平方向上最靠近于B点的点，扫描仪从A点扫描至F点所转过的角度∠FOA记为 $\mathrm{\β}=(2*[\frac{\arctan ((\sqrt{D_M^2{-D}_V^2}-D_G)/D_V)-\mathrm{\α}}{2*\mathrm{\α}}]+3)*\mathrm{\α},$ 则BC上点数为其中∠AOC＝arccos(D_M/D_V)。

记当前为N_col列中的第i列，其与水平方向AC的交点为I点，与建筑物顶部交点为M点，与建筑物底部交点为N点，则扫描仪从B点扫描至I点所转过的角度∠BOI＝β+i*2*α，第i列点云的数量包括仪器水平线(0点)以上的点数和仪器水平线(0点)以下的点数，因此第i列点云的数量线段MI上的LiDAR点数为： $N_{\mathrm{Above}}^i=\left[\frac{\∠\mathrm{MOI}-\mathrm{\α}}{2*\mathrm{\α}}\right]+1,$ $\∠\mathrm{MOI}=\arctan \left(\frac{H_B-H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right);$ 线段NI上的LiDAR点数为 $N_{\mathrm{Below}}^i=\left[\frac{\∠\mathrm{NOI}-\mathrm{\α}}{2*\mathrm{\α}}\right]+1,$ $\∠\mathrm{NOI}=\arctan \left(\frac{H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right),$ 故而，得到第i列LiDAR点的数量为：

$N_{\mathrm{row}}^i=\left[\frac{\arctan \left(\frac{H_B-H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right)-\mathrm{\α}}{2*\mathrm{\α}}\right]+\left[\frac{\arctan \left(\frac{H_L}{D_V/\cos (\mathrm{\β}+i*2*\mathrm{\α})}\right)-\mathrm{\α}}{2*\mathrm{\α}}\right]+3,$ 网格内的每列LiDAR点数量相加得到格网处点云总数即为格网点云密度。

由于墙面的凹凸起伏和墙面窗户的反射，同一侧墙面点云投射到XY平面时，其格网密度也会产生较大差异。为了保证所有的墙面格网都有效提取，需要设定一个密度阈值对个网进行删选，网格密度阈值threshod＝rate*N，式中，参数rate取值范围为0.2-1。

本实施例使用轮廓延伸密度的方法进行建筑物轮廓的恢复，具体如下：

3a)寻找步骤1e)中获得的二维矢量轮廓线段周边1m范围内格网，将寻找到的所有格网内LiDAR点最大高程的平均值作为二维矢量轮廓线段的高程，将二维矢量轮廓线段变换为三维建筑物轮廓线段；

3b)对三维建筑物轮廓线段构建半径为1m的缓冲区，建缓冲区内LiDAR点数量除以缓冲区体积获得原有轮廓LiDAR点密度；

3c)沿轮廓线段方向以单位距离为延伸步长构建半径为1m的缓冲区，缓冲区内LiDAR点数量除以相应缓冲区体积获得待延伸方向的LiDAR点密度，本例中单位距离的取0.2m；延伸的单位距离越小，精度越高；

3d)若待延伸方向的LiDAR点密度与原有轮廓LiDAR点密度的差异小于20％，则该轮廓沿轮廓线段方向延伸单位距离并重复步骤3c)；否则停止延伸，形成完整的建筑物轮廓。恢复后的地面轮廓如图3-b所示。

本例中，从航空LiDAR数据中提取到航空轮廓线段103条，地面轮廓线段31条。

第二步、提取建筑物角点——从航空轮廓中提取建筑物角点，称为航空角点；从地面轮廓中提取建筑物角点，称为地面角点。

提取建筑物角点的方法如下：将建筑物轮廓投影到三维坐标系的XY平面内寻找二维相交点，如果任两条构成相交点的轮廓的高程差小于1m，则判定两条轮廓在实际的三维空间中相交，两条轮廓的相交点为一个建筑物角点，并将所述两条轮廓的高程均值作为该建筑物角点的高程。用这种方法分别从地面轮廓与航空轮廓中提取出航空角点A＝{A_i，i＝0，1，2，...，u}如图4-a所示，图中中心带黑点的圆圈代表航空角点，地面角点B＝{B_i，i＝0，1，2，...，v}如图4-b所示，图中三角形代表地面角点。航空角点提取到58个，地面角点提取到15个。

第三步、寻找轮廓线段约束下的初始转换矩阵——使用航空角点与地面角点迭代计算转换矩阵，用该转换矩阵对地面轮廓进行转换，并使用航空轮廓与转换后地面轮廓的匹配度作为控制约束条件，当航空轮廓与转换后地面轮廓之间成功匹配的线段对数满足给定阈值时停止迭代，相应的转换矩阵即为初始转换矩阵。如图5-a所示，为实施轮廓线作为约束条件下的建筑物角点匹配后的角点示意图，图中中心带黑点的圆圈代表航空角点，黑色五角星代表初始配准后的地面角点，图5-b为相应的轮廓线示意图，图中灰色实线为截取的部分航空轮廓，黑色虚线代表初始匹配后的地面轮廓。

本步骤中寻找初始转换矩阵的具体方法如下：

4a)设航空角点的点集分别为A＝{A_i，i＝0，1，2，...，u}；地面角点的点集为B＝{B_i，i＝0，1，2，...，v}，航空轮廓的线段集为LA＝{LA_i，i＝0，1，2，...m}；地面轮廓的线段集为LB＝{LB_i，i＝0，1，2，...，n}，u为航空角点的数量，v为地面角点的数量，m为航空轮廓的线段数量，n为地面轮廓的线段数量；

4b)从点集A和B中分别选取1个点A_x和B_x，计算由B_x至A_x的平移矩阵，利用该平移矩阵对点集B中的每个点进行平移，得到点集M＝{M_i，i＝0，1，2，...，v}；

4c)从点集A和M中分别选取1个点A_y和M_y，要求A_y≠A_x，M_y≠B_x，以点A_x为原点，计算点M_y旋转至点A_y位置的旋转矩阵，使用该旋转矩阵对点集M中每个点进行旋转，得到点集R＝{R_i，i＝0，1，2，...，v}；

4d)使用步骤4b)中得到的平移矩阵和步骤4c)中得到的旋转矩阵对地面轮廓的线段集LB进行转换，得到转换后的轮廓线段集LC＝{LC_i，i＝0，1，2，...，n}；

4e)遍历线段集LA中的所有线段，寻找线段集LA和线段集LC之间满足给定匹配条件的线段对数量：

从线段集LA中取出一条未检验的线段，计算其与线段集LC中所有尚未匹配轮廓线段的互异度，若最小互异度小于1，则与该最小互异度对应的两条轮廓线为已匹配线段对，重复本过程直至线段集LA中不存在未检验的线段；两条轮廓线段的互异度计算公式如下：

dif＝w₁×lenDif+w₂×lDis+w₃×lAng+w₄×cpDis

其中，dif为两条轮廓线段的互异度；lenDif为两条轮廓线段的长度之差，

lDis为两条轮廓线段所在直线的距离，lAng为两条轮廓线段所在直线的夹角， $\mathrm{lAng}=(\frac{\mathrm{\π}}{2}-|\arccos \left(\frac{{\stackrel{\→}{l}}_1\•{\stackrel{\→}{l}}_2}{\left|{\stackrel{\→}{l}}_1\right|\•\left|{\stackrel{\→}{l}}_2\right|}\right)-\frac{\mathrm{\π}}{2}\left|\right)*\frac{180}{\mathrm{\π}};$ cpDis为两轮廓线中点之间的距离， $\mathrm{cpDis}=\sqrt{{(x_{l_{11}}-x_{l_{21}})}^2+{(y_{l_{11}}-y_{l_{21}})}^2+{(z_{l_{11}}-z_{l_{21}})}^2};$ w1、w2、w3、w4分别表示以上4个参量的权重，取值范围分别为[0.4，0.6]，[0.1，0.3]，[0.1，0.2]，[0.1，0.3]，本例中w1、w2、w3、w4取值分别为0.5、0.2、0.1、0.2，分别代表两条轮廓线段的方向向量：

${\stackrel{\→}{l}}_1=(x_{l_{12}}-x_{l_{11}},y_{l_{12}}-y_{l_{11}},z_{l_{12}}-z_{l_{11}})$

${\stackrel{\→}{l}}_2=(x_{l_{22}}-x_{l_{21}},y_{l_{22}}-y_{l_{21}},z_{l_{22}}-z_{l_{21}})$

其中分别代表LA中参与匹配的轮廓线段的起点和终点x、y、z坐标，分别代表LC中参与匹配的轮廓线段的起点和终点x、y、z坐标；

4f)若线段集LA和线段集LC中匹配的线段对数量不少于6对，则认为该转换矩阵可靠，该转换矩阵为初始转换矩阵；否则转至步骤4b)重新进行转换矩阵计算。

经本步骤后即可成功找到初始转换矩阵，包括旋转矩阵R和平移矩阵T。

第四步、获得初始匹配角点对——使用初始转换矩阵(旋转矩阵R和平移矩阵T)对地面角点点集B进行转换得到点集F＝{F_i＝R×B_i+T，i＝0，1，2，...，v}，从航空角点的点集A＝{A_i，i＝0，1，2，...，u}和点集F＝{F_i，i＝0，1，2，...，v}选出成功配对的角点，分别记录为点集P＝{P_i，i＝0，1，2，...，m}和点集U＝{U_i，i＝0，1，2，...，n}，其中m＝n，得到初始匹配角点对。本例中共得到12对匹配成功的角点，如图6所示，灰色圆圈代表匹配成功的航空角点，黑色五角星代表匹配成功的地面角点。

第五步、寻找修正转换矩阵——以初始匹配角点对为源数据，寻找两者间的修正转换矩阵，保证经该修正转换矩阵配准后，两者的均方根误差小于预设的极限值ε；ε，极限值ε的取值范围为0.25-0.35，本实施例中，ε取0.3。

本步骤具体方法如下：

5a)在点集P中找出距点集U中每一个点的最近点，组成点集Q＝{Q_i，i＝0，1，2，...n}；

5b)采用最小均方根法计算点集U与点集Q之间的配准关系，得到配准转换矩阵；

5c)对点集U用配准转换矩阵进行坐标转换，得到点集U₁；

5d)计算点集U₁与点集Q之间的均方根误差，如小于预设的极限值ε，本例中极限值ε取0.3，则以该配准转换矩阵作为修正转换矩阵；否则，以点集U₁替换U，转至步骤5a)重新进行配准转换矩阵计算。

本步骤可以获得修正转换矩阵，包括旋转矩阵R′与平移矩阵T′，R′是3×3的旋转矩阵，T′是3×1的平移矩阵。

第六步、LiDAR数据配准——使用初始转换关系与修正转换关系依次对地面LiDAR数据进行转换，得到最终配准结果：

用初始转换矩阵对地面LiDAR数据PB进行初始转换，得到初始配准点云数据PB′，PB′＝{PB′_i＝R×PB_i+T，i＝1，2，3，…，C_B}；用修正转换矩阵初始配准点云数据PB′进行修正转换，得到最终配准点云数据PB″，PB″＝{PB″_i＝R′×PB′_i+T′，i＝1，2，3，…，C_B}，其中C_B为地面LiDAR数据中LiDAR点的数量。配准结果如图7所示。

验证实施例：

下面为验证本发明方法可靠性，以实例进行说明。

如图8-a所示，为本实例提取得到的地面轮廓与地面角点(以三角形表示)，由图可知，从地面LiDAR数据中提取出的角点较少，而轮廓线较为丰富；如图8-b所示，为从地面LiDAR数据中提取的航空轮廓与航空角点(中心为黑点的圆圈)。根据轮廓间的位置关系可以很明显地看出角点配准后的正确位置应该如图8-c的A区域，利用本发明方法进行自动配准得到的结果符合前述判断(配准的角点落在A区)，而单独使用角点进行配准时，配准的角点(黑色十字符号表示)落在图8-d的B区内，属于明显的错误。出现图8-d的B区域中所示错误的原因在于可供计算转换关系的公共角点数量太少，无法验证结果的正确性，而多数情况下，从LiDAR数据中提取出的建筑物轮廓线会比角点多出不少，因此使用轮廓做约束可以比较稳定地控制配准结果的精度，避免由于角点数据过少引起的配准错误。

可见，本发明方法可以提高航空和地面Li DAR数据的配准可靠性，避免由于可参考的建筑物角点较少而导致的匹配错误，并且实现了自动化提取。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。

Claims (5)

1.一种基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，包括以下步骤：

第一步、提取建筑物轮廓——从航空LiDAR数据中提取建筑物轮廓，称为航空轮廓；从地面LiDAR数据中提取建筑物轮廓，称为地面轮廓；

第二步、提取建筑物角点——从航空轮廓中提取建筑物角点，称为航空角点；从地面轮廓中提取建筑物角点，称为地面角点；

第三步、寻找轮廓线段约束下的初始转换矩阵——使用航空角点与地面角点迭代计算转换矩阵，用该转换矩阵对地面轮廓进行转换，并使用航空轮廓与转换后地面轮廓的匹配度作为控制约束条件，当航空轮廓与转换后地面轮廓之间成功匹配的线段对数满足给定阈值时停止迭代，相应的转换矩阵即为初始转换矩阵；

第四步、获得初始匹配角点对——使用第三步获得的初始转换矩阵对地面角点进行转换，根据空间距离寻找航空角点中与其配对的角点，得到初始匹配角点对；

第五步、寻找修正转换矩阵——以初始匹配角点对为源数据，寻找两者间的修正转换矩阵，保证经该修正转换矩阵配准后，两者的均方根误差小于预设的极限值ε，极限值ε的取值范围为0.25-0.35；

第六步、LiDAR数据配准——使用初始转换矩阵与修正转换矩阵依次对地面LiDAR数据进行转换，得到最终配准结果。

2.根据权利要求1所述的基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，其特征在于：所述第一步从航空LiDAR数据中提取建筑物轮廓的方法如下：构建1m*1m的水平格网，根据点面空间关系计算每个格网中最高点与最低点的高差，保留高差大于实验区最低建筑物高程的格网得到轮廓格网，对轮廓格网使用多尺度的Hough变换，得到航空LiDAR建筑物轮廓线段。

3.根据权利要求1所述的基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，其特征在于：第一步中从地面LiDAR数据中提取建筑物轮廓的方法如下：使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓；在此基础上使用轮廓延伸密度方法对提取的建筑物轮廓进行恢复，形成完整的建筑物轮廓。

4.根据权利要求1所述的基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，其特征在于：第二步中提取建筑物角点的方法如下：将建筑物轮廓投影到三维坐标系的XY平面内寻找二维相交点，如果任两条构成相交点的轮廓的高程差小于1m，则判定两条轮廓在实际的三维空间中相交，两条轮廓的相交点为一个建筑物角点，并将所述两条轮廓的高程均值作为该建筑物角点的高程。

5.根据权利要求3所述的基于建筑物轮廓约束的航空与地面LiDAR数据自动配准方法，其特征在于：使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓，具体步骤如下：

1a)地面LiDAR点云投影至XY平面——将地面LiDAR点云投影至三维坐标系的XY平面，并保留各个点的X、Y、Z属性；

1b)提取粗略轮廓格网——在所述XY平面内构建1m*1m的粗略格网，计算每个粗略格网中LiDAR投影点的数量，即得到该粗略格网的格网密度，根据建筑物边缘轮廓处的粗略格网密度阈值对所述粗略格网进行筛选，保留格网密度大于所述粗略格网密度阈值的粗略格网，得到粗略轮廓格网；

1c)提取精确轮廓格网——在粗略轮廓格网中构建0.2m*0.2m的精细格网，计算精细格网内LiDAR投影点的数量即得到精细格网的格网密度，根据建筑物边缘轮廓处的精细格网密度阈值对所述精细格网进行筛选，保留格网密度大于所述精细格网密度阈值的精细格网，得到精确轮廓格网；

1d)格网高差筛选——遍历所有精确轮廓格网，如果精确轮廓格网内的最高LiDAR点和最低LiDAR点的高差大于相应实验区建筑最低高程则保留该精确轮廓格网，否则剔除；

1e)获取轮廓线段——对筛选后的精确轮廓格网使用Hough变换得到二维矢量轮廓线段。