CN103065295B - 一种基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法 - Google Patents

Abstract

一种基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，首先从航空LiDAR数据提取建筑物角点（称为航空角点）；再从地面LiDAR数据提取建筑物角点（称为地面角点）；然后对航空角点与地面角点进行初始匹配，分别从航空角点和地面角点中选取任意两个点计算用于坐标变换的转换矩阵，对所有转换矩阵根据最大匹配对数和最小匹配误差确定最优转换矩阵；最后在确定初始匹配角点对基础上，以地面角点为参考，对航空LiDAR角点进行位置修正，最终实现航空LiDAR数据与地面LiDAR数据的自动配准。本发明可修正航空LiDAR角点中误差较大的点，较大地提高了航空LiDAR数据和地面LiDAR数据的配准精度。

Description

一种基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法

技术领域

本发明涉及一种航空和地面LiDAR数据配准方法，特别是涉及一种基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法。

背景技术

激光雷达技术，作为一种大范围的面状三维坐标测量工具，能够提供地物表面详尽而准确的不规则LiDAR点，数据质量优于其他一些测量技术，如数字摄影测量、雷达干涉测量等。目前，使用最多的激光雷达是航空激光雷达和地面激光雷达。航空激光雷达俯视地获取大范围的数据，数据中包含大量地物的顶部信息，然而侧面信息获取困难，数据分辨率也有限；地面激光雷达平视或仰视地获取区域数据，地物侧面信息详尽，数据分辨率较高，然而扫描范围有限，地物顶部信息难以获得。航空与地面激光雷达各有优缺点，为了能够获取地物各个尺度、各个方向的详尽数据，两种平台数据的融合已经成为一种必然趋势。事实上，这些年来越来越多的学者集成两种数据进行应用研究，相关研究涉及地形制图、地质勘探、森林研究、水文研究、以及虚拟现实等等。

为了集成航空和地面LiDAR数据进行研究，首先需要对两种数据进行配准。事实上，航空和地面LiDAR数据存在着较大差异，从两种数据中提取的公共角点数量有限，精度差异也较大：1)不同视角。航空LiDAR俯视获取数据，顶部信息较多，侧面信息较少，而地面LiDAR平视、仰视获取数据，侧面信息详尽，顶部信息较少，两者公共角点较少。2)不同范围。航空LiDAR能够获取大范围的数据，提高大量的配准角点，而地面LiDAR扫描范围有限，提供的配准角点较少，位置集中。3)不同分辨率。航空LiDAR距离扫描目标几百米至上千米不等，数据分辨率在米级或分米级，提取的角点精度较低，而地面LiDAR距离扫描目标几十米，数据分辨率最高可达毫米级，提取的角点精度较高。如何充分利用少量的精度差异较大的公共角点，实现航空和地面LiDAR更为准确的配准还有待于研究。

发明内容

本发明要解决技术问题是：克服现有技术缺点，提出基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，该方法对误差较大的航空角点进行自动移动修正，从而提高航空角点的精度，以便实现航空和地面LiDAR数据的高精度配准。

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

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

第二步、分别从航空角点和地面角点中选取出任意两个点计算用于坐标变换的转换矩阵，对所有可能的转换矩阵根据最大匹配对数和最小匹配误差确定最优转换矩阵，并使用最优转换矩阵对地面角点进行初始变换，对应的成功匹配的航空角点和地面角点即为初始匹配角点对；

第三步、计算每对航空和地面角点之间的距离，将具有最大距离值的航空角点移动至对应的地面角点位置，完成一次误差较大航空角点的修正；

第四步、循环执行第二步和第三步，在循环过程中，计算航空角点与地面角点的总误差值，如果本次循环的总误差值大于前面一次循环总误差值的20％，则停止循环，完成所有误差较大航空角点移动修正；

第五步、采用倒数第二次循环中获得的最优转换矩阵对航空LiDAR数据和地面LiDAR数据进行配准。

本发明可修正航空角点中误差较大的点，较大地提高了航空LiDAR数据和地面LiDAR数据的配准精度，利用本发明可实现航空与地面LiDAR数据的高精度自动配准。

本发明第一步中，从航空LiDAR数据中提取航空角点的方法如下：使用反向迭代数学形态学方法从航空LiDAR数据中提取建筑物区域；对不规则建筑物区域使用规则化方法得到规则的建筑物轮廓线段，两两相交的轮廓线段形成航空角点，并将所述航空角点周围1m范围内的LiDAR点的高程最高值作为该航空角点的高程。

本发明第一步中从地面LiDAR数据中提取地面角点的方法如下：使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓；在此基础上使用轮廓延伸密度方法对提取的建筑物轮廓进行恢复，形成完整的建筑物轮廓；然后将完整的建筑物轮廓投影到三维坐标系的XY平面内寻找二维相交点，如果任两条构成相交点的轮廓的高程差小于1m，则判定两条轮廓在实际的三维空间中相交，两条轮廓的相交点为一个地面角点，并将所述两条轮廓的高程均值作为该地面角点的高程。

其中，使用分层次的格网密度方法从地面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变换得到二维矢量轮廓线段。。

上述步骤1a)和1b)中，格网密度阈值的确定方法，具体步骤如下：

假设0点为仪器中心点，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点的列数，记α为扫描仪每次旋转角度的一半，记格网中水平方向上最靠近于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。

使用轮廓延伸密度的方法进行建筑物轮廓恢复的具体步骤如下：

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

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

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

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

本发明在第二步中根据最大匹配对数和最小匹配误差方法确定初始变换，假设两个点集分别为A＝{A_i，i＝0，1，2，...，u}和B＝{B_i，i＝0，1，2，...，v}，本发明中初始变换的具体方法如下：

4a)从点集A、B中分别选取一个点A₁、B₁，计算平移矩阵T，利用平移矩阵T对点集B中的每个点进行平移，得到点集M＝{M_i，i＝0，1，2，...，v}；

4b)从点集A、M中分别选取一个点A₂、M₂，要求A₂≠A₁，M₂≠B₁，以点A₁为原点，计算点M₂旋转至点A₂位置的旋转矩阵R，使用旋转矩阵R对点集M中每个点进行旋转，得到点集R＝{R_i，i＝0，1，2，...，v}；

4c)使用平移矩阵T和旋转矩阵R对点集B中所有点转换，得到点集C＝{C_i，i＝0，1，2，...，v}；对于点集A中每一点A_i在点集C中寻找与其最近的点C_closet，如果点A_i和C_closet距离小于X米，X取值范围为1-5，则两点匹配；如果点集A中存在两个点在点集C中的最近点为同一个点，则两者距离最近的两点匹配；记录成功匹配的点对MA＝{MA_i，i＝1，2，...，n}和MB＝{MB_i，i＝1，2，...，n}；

4d)重复执行步骤4b)和4c)，计算所有可能的转换矩阵，选取匹配对数最多的转换矩阵作为待选可靠转换矩阵；

4e)对于每组待选可靠转换矩阵，计算匹配点对MA和MB总误差，其中误差最小的待选可靠转换矩阵为最优转换矩阵；

4f)使用最优转换矩阵对点集B中所有点进行转换，完成初始变换，对应的成功匹配的航空和地面角点，即为初始匹配角点对。

本发明还提供了误差较大航空角点的修正方法，设第一步中提取的航空角点集为P＝{P_i，i＝0，1，2，...，n}，地面角点集分别为U＝{U_i，i＝0，1，2，...，n}，第二步中获得的最优转换矩阵包括旋转矩阵R和平移矩阵T，则第三步中，对误差较大航空角点的修正方法具体如下：

5a)利用旋转矩阵R和平移矩阵T对地面角点集U进行转换得到转换后的地面角点，记为Q＝{Q_i，i＝0，1，2，...，n}；

5b)计算航空角点集P和转换后的地面角点集Q对应角点的三维空间距离，得到一维距离矩阵D＝{D(P_i，Q_i)，i＝0，1，2，...，n}，计算角点总匹配误差

5c)寻找一维距离矩阵D中距离的最大值D_max，并找到航空角点集P和转换后的地面角点集Q中对应的角点P_max和Q_max，将航空角点P_max的坐标移动至转换后的地面角点Q_max的坐标，完成一次误差较大航空角点的移动修正。

本发明的有益效果是：1)本发明提供了一种建筑物角点自修正的方法，通过对航空和地面角点进行ICP迭代，修正航空角点中误差较大的点，有效提高了航空角点的精度，从而为实现高精度的航空和地面LiDAR数据的配准提供了条件；2)本发明使用一种分层次格网密度方法提取地面LiDAR数据的建筑物轮廓，并使用理论估计的方法对格网密度阈值进行确定，能够从地面LiDAR数据中提取准确的建筑物轮廓线段，从而提取高精度的地面角点；3)本发明对从地面LiDAR数据中提取的建筑物轮廓进行了恢复，得到完整的建筑物轮廓，从而可以获得准确的地面角点，进一步提高了数据配准的准确性和精度。

附图说明

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

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

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

图3为从图1中提取的建筑物示意图。

图4为图3中提取的航空LiDAR数据轮廓及航空角点示意图。

图5为图2中提取的建筑物粗略轮廓格网示意图。

图6为经高程筛选后得到的建筑物精细轮廓格网示意图。

图7为格网密度阈值理论估算示意图。

图8为图6中提取的地面LiDAR数据轮廓示意图。

图9为图8恢复后的地面LiDAR数据轮廓及地面角点示意图。

图10为航空角点与地面角点的初始匹配示意图。

图11为修正后的航空角点示意图。

图12为使用修正后的航空角点对航空和地面LiDAR数据配准结果示意图。

具体实施方式

本实施例的航空和地面LiDAR数据分别如图1和图2所示，航空LiDAR数据平均点间距为1m，高程精度为15cm，平面精度为30cm，LiDAR点总数约46万；地面LiDAR数据使用LeicaScanStation2分9站扫描获得，LiDAR点总数约3000万，地面LiDAR点分辨率为100m远处20cm一个点。

本实施例的基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，包括以下步骤：

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

本步骤中从航空LiDAR数据中提取航空角点方法是：使用反向迭代数学形态学方法从航空LiDAR数据中提取建筑物区域；对不规则建筑物区域进行规则化得到规则建筑物边缘轮廓，轮廓两两相交得到航空角点，航空角点周围1m范围内的LiDAR点的高程最高值作为该航空角点的高程。

本实例使用反向迭代数学形态学方法的打开操作提取建筑物，对于LiDAR数据来说，打开操作可以消除点集中小于指定窗口大小高值区域。这里设定反向迭代数学形态学方法的初始形态学窗口大小为106，窗口减少的步长为10，高差筛选阈值为15m，粗糙度阈值为1.6，得到的建筑物区域如图3所示，共10个建筑物。

使用规则化方法对10个不规则的面状的建筑物区域进行规则化，得到了如图4所示的规则化结果以及航空角点，共72个角点，相对于实际航空角点平均误差0.91m。

本步骤中从地面LiDAR数据中提取地面角点方法是：使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓；在此基础上使用轮廓延伸密度方法对提取的建筑物轮廓进行恢复，形成完整的建筑物轮廓；然后将完整的建筑物轮廓投影到三维坐标系的XY平面内寻找二维相交点，如果任两条构成相交点的轮廓的高程差小于1m，则判定两条轮廓相交，两条轮廓的相交点为一个地面角点，并将所述两条轮廓的高程均值作为该地面角点的高程。

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

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的所有格网都为粗略轮廓格网。本例提取的粗略轮廓格网见图5。

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

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

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

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

本实施例在上述步骤1a)和1b)中，使用理论估计的方法计算格网密度阈值，如图7所示，假设0点为仪器中心点，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点的列数，记α为扫描仪每次旋转角度的一半，记格网中水平方向上最靠近于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。

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

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

那么， $\∠\mathrm{AOE}=\left[\frac{\∠\mathrm{AOB}-\mathrm{\α}}{2*\mathrm{\α}}\right]+\mathrm{\α},$ 其中 $\∠\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.1-0.3m；

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

延伸的单位距离越小，精度越高。本实例中，单位距离取值0.2m。

本实施例恢复后的轮廓(已转化为二维)如图9所示，共有30条轮廓，相交为16个角点。

第二步、对航空角点与地面角点进行初始匹配——分别从航空角点和地面角点中选取出任意两个点计算用于坐标变换的转换矩阵，对所有可能的转换矩阵根据最大匹配对数和最小匹配误差确定最优转换矩阵，并使用最优转换矩阵对地面角点进行初始变换，对应的成功匹配的航空角点和地面角点即为初始匹配角点对。

本实施例根据最大匹配对数和最小匹配误差方法确定初始变换，由于航空角点共72个，地面角点共16个，假设两个点集分别为A＝{A_i，i＝0，1，2，...，72}和B＝{B_i，i＝0，1，2，...，16}，则初始变换的具体方法如下：

4a)从点集A、B中分别选取一个点A₁、B₁，计算平移矩阵T，利用平移矩阵T对点集B中的每个点进行平移，得到点集M＝{M_i，i＝0，1，2，...，16}；

4b)从点集A、M中分别选取一个点A₂、M₂，要求A₂≠A₁，M₂≠B₁，以点A₁为原点，计算点M₂旋转至点A₂位置的旋转矩阵R，使用旋转矩阵R对点集M中每个点进行旋转，得到点集R＝{R_i，i＝0，1，2，...，16}；

4c)使用平移矩阵T和旋转矩阵R对点集B中所有点转换，得到点集C＝{C_i，i＝0，1，2，...，16}；对于点集A中每一点A_i在点集C中寻找与其最近的点C_closet，如果点A_i和C_closet距离小于3米，则两点匹配；如果点集A中存在两个点在点集C中的最近点为同一个点，则两者距离最近的两点匹配；记录成功匹配的点对MA＝{MA_i，i＝1，2，...，n}和MB＝{MB_i，i＝1，2，...，n}；

4d)重复执行步骤4b)和4c)，对于所有的可能的转换矩阵，实例中共306720个转换矩阵，选取匹配对数最多的转换矩阵作为待选可靠转换矩阵，实例中匹配对数最多的转换矩阵有6084个，匹配对数为13对；

4e)对于这6084个待选可靠转换矩阵，计算其匹配点对MA和MB的距离，其中距离最小的待选可靠转换矩阵为最优转换矩阵，得到一个最优转换矩阵，其中按Z轴旋转149.3°，按Y轴旋转3.9°，在XZ平面旋转28.1°，平移矩阵为[-128569.2，148434.5，21.9]；

4f)使用最优转换矩阵对地面角点集B中所有点进行转换，完成初始变换。

本实施例对航空角点和地面角点进行初始变换的结果如图10所示，其中三角形为航空角点，圆为地面角点。

第三步、对航空角点进行移动修正——计算每对航空和地面角点之间的距离，将具有最大距离值的航空角点移动至对应的地面角点位置，完成一次误差较大航空角点的修正。

假设提取的航空角点集和地面角点集分别为P＝{P_i，i＝0，1，2，...，n}和U＝{U_i，i＝0，1，2，...，n}，则对误差较大航空角点的修正方法如下：

5a)利用第二步中获得的最优转换矩阵对地面角点集U进行转换得到转换后的地面角点，记为Q＝{Q_i，i＝0，1，2，...，n}；

5b)计算航空角点集P和转换后的地面角点集Q对应角点的三维空间距离，得到一维距离矩阵D＝{D(P_i，Q_i)，i＝0，1，2，...，n}，计算角点总匹配误差

5c)寻找一维距离矩阵D中距离的最大值D_max，并找到航空角点集P和转换后的地面角点集Q中对应的角点P_max和Q_max，将航空角点P_max的坐标移动至转换后的地面角点Q_max的坐标，完成一次误差较大航空角点的移动修正。

第四步、循环执行第二步和第三步，在循环过程中，计算航空角点与地面角点的总误差值，如果本次循环的总误差值大于前面一次循环总误差值的20％，则停止循环，完成所有误差较大航空角点移动修正。

在本实例中，共循环了10次，共对航空角点进行了9次修正，修正的航空角点编号依次为：2、1、7、4、5、12、8、2、9、8，经过修正后的角点如图11所示，其中五边形为实际航空角点位置，三角形为修正后航空角点位置，经过修正后航空角点平均误差为0.32m，比原始的航空角点精度提高0.6m。

第五步、采用倒数第二次循环中获得的最优转换矩阵对航空LiDAR数据和地面LiDAR数据进行配准。

在本实例中，使用修正后的航空角点和地面角点进行配准，最终得到的配准结果如图12所示，经过修正后的角点进行配准，角点绝对位置误差仅0.29m，配准精度较高。

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

Claims (7)

1.一种基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，包括以下步骤：

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

第二步、分别从航空角点和地面角点中选取出任意两个点计算用于坐标变换的转换矩阵，对所有可能的转换矩阵根据最大匹配对数和最小匹配误差确定最优转换矩阵，并使用最优转换矩阵对地面角点进行初始变换，对应的成功匹配的航空角点和地面角点即为初始匹配角点对；

第三步、计算每对航空和地面角点之间的距离，将具有最大距离值的航空角点移动至对应的地面角点位置，完成一次误差较大航空角点的修正；

第四步、循环执行第二步和第三步，在循环过程中，计算航空角点与地面角点的总误差值，如果本次循环的总误差值大于前面一次循环总误差值的20％，则停止循环，完成所有误差较大航空角点移动修正；

第五步、采用倒数第二次循环中获得的最优转换矩阵对航空LiDAR数据和地面LiDAR数据进行配准。

2.根据权利要求1所述的基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，其特征在于：所述第一步从航空LiDAR数据中提取航空角点的方法如下：使用反向迭代数学形态学方法从航空LiDAR数据中提取建筑物区域；对不规则建筑物区域使用规则化方法得到规则的建筑物轮廓线段，两两相交的轮廓线段形成航空角点，并将所述航空角点周围1m范围内的LiDAR点的高程最高值作为该航空角点的高程。

3.根据权利要求1所述的基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，其特征在于：第一步中从地面LiDAR数据中提取地面角点的方法如下：使用分层次的格网密度方法从地面LiDAR数据中提取建筑物轮廓；在此基础上使用轮廓延伸密度方法对提取的建筑物轮廓进行恢复，形成完整的建筑物轮廓；然后将完整的建筑物轮廓投影到三维坐标系的XY平面内寻找二维相交点，如果任两条构成相交点的轮廓的高程差小于1m，则判定两条轮廓在实际的三维空间中相交，两条轮廓的相交点为一个地面角点，并将所述两条轮廓的高程均值作为该地面角点的高程。

4.根据权利要求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变换得到二维矢量轮廓线段。

5.根据权利要求4所述的基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，其特征在于，步骤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点的列数，记α为扫描仪每次旋转角度的一半，记格网中水平方向上最靠近于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.根据权利要求1所述的基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，其特征在于，第二步中根据最大匹配对数和最小匹配误差方法确定初始变换；设两个点集分别为A＝{A_i,i＝0,1,2,...,u}和B＝{B_i,i＝0,1,2,...,v}，则初始变换的具体方法如下：

4a)从点集A、B中分别选取一个点A₁、B₁，计算平移矩阵T，利用平移矩阵T对点集B中的每个点进行平移，得到点集M＝{M_i,i＝0,1,2,...,v}；

4b)从点集A、M中分别选取一个点A₂、M₂，要求A₂≠A₁,M₂≠B₁，以点A₁为原点，计算点M₂旋转至点A₂位置的旋转矩阵R，使用旋转矩阵R对点集M中每个点进行旋转，得到点集R＝{R_i,i＝0,1,2,...,v}；

4c)使用平移矩阵T和旋转矩阵R对点集B中所有点转换，得到点集C＝{C_i,i＝0,1,2,...,v}；对于点集A中每一点A_i在点集C中寻找与其最近的点C_closet，如果点A_i和C_closet距离小于X米，X取值范围为1-5，则两点匹配；如果点集A中存在两个点在点集C中的最近点为同一个点，则两者距离最近的两点匹配；记录成功匹配的点对MA＝{MA_i,i＝1,2,...,n}和MB＝{MB_i,i＝1,2,...,n}；

4d)重复执行步骤4b)和4c)，计算所有可能的转换矩阵，选取匹配对数最多的转换矩阵作为待选可靠转换矩阵；

4e)对于每组待选可靠转换矩阵，计算匹配点对MA和MB总误差，其中误差最小的待选可靠转换矩阵为最优转换矩阵；

4f)使用最优转换矩阵对点集B中所有点进行转换，完成初始变换，对应的成功匹配的航空和地面角点，即为初始匹配角点对。

7.根据权利要求1所述的基于建筑物角点自修正的航空与地面LiDAR数据高精度自动配准方法，其特征在于：设第一步中提取的航空角点集为P＝{P_i,i＝0,1,2,...,n}，地面角点集为U＝{U_i,i＝0,1,2,...,n}，则第三步中，对误差较大航空角点的修正方法具体如下：

5a)利用第二步中获得的最优转换矩阵对地面角点集U进行转换得到转换后的地面角点，记为Q＝{Q_i,i＝0,1,2,...,n}；

5b)计算航空角点集P和转换后的地面角点集Q对应角点的三维空间距离，得到一维距离矩阵D＝{D(P_i,Q_i),i＝0,1,2,...,n}，计算角点总匹配误差

5c)寻找一维距离矩阵D中距离的最大值D_max，并找到航空角点集P和转换后的地面角点集Q中对应的角点P_max和Q_max，将航空角点P_max的坐标移动至转换后的地面角点Q_max的坐标，完成一次误差较大航空角点的移动修正。