1 引 言

航空影像语义信息丰富，是地物提取最大的数据源之一，由于缺乏直接的三维信息，单独利用影像进行数据处理有一定的局限性。机载激光雷达(airborne light detection and ranging，LiDAR)是一种新型的航空对地观测传感器，具有主动观测、受天气影响小、不受阴影影响、对地物间缝隙有一定穿透性并能直接获取地表三维几何信息等优点^[1]。将航空影像和LiDAR点云配准、融合，可同时获取地物的空间信息及语义信息，有效弥补了单一数据源的不足^[2]，在灾害评估、真正射影像制作、道路、建筑物自动建模等领域有很大的应用前景^{[3, 4, 5, 6]}。

与传统配准算法不同，航空影像与机载LiDAR点云的配准，是二维与三维数据、连续影像与离散点云之间的配准问题，传统影像之间的配准方法收效甚微^[7]。近10年来，国内外学者对其进行了一定的研究，大致可以归结为3类算法。①通过航空影像密集匹配获取三维点集，然后通过ICP(iterative closest point)或其改进方法^{[9, 10]}建立其与机载LiDAR点云的对应关系，这本质上是一种3D-3D配准变换。这类算法无法进行单张影像与点云的配准，同时密集匹配生成点云的过程会降低配准精度。②将LiDAR点云按照高程或强度转换为二维影像，再依据传统的图像配准方法实现两种数据的配准。如基于互信息的全局配准方法^[11]、基于SIFT(scale invariant feature transform)算子的配准方法^[12]、基于FFT变换的算法^[13]、6元组松弛法^[14]等。这类算法充分利用了原有的成熟图像配准算法，自动化程度高，然而规则化过程中的误差、点位选取不准的误差，都会降低配准精度。③部分学者尝试建立遥感影像与机载LiDAR点云之间的直接配准模型。文献[15—16]用直线特征或平面特征建立遥感影像与激光点云的直接映射关系。文献^[17]将空间三维点采用参数方程和辅助变量联合表示，通过扩展的共线方程，建立起二维影像点(2D)与的点云空间三维直线(3D)的直接映射关系，避免了LiDAR空间点位选取不准的误差、点云内插误差以及密集匹配的误差，配准模型严密、简单，具有一定的推广价值。该算法的自动化程度取决于配准基元的选取与匹配的自动化程度，需要进一步的研究。

基于上述分析，本文综合第2类和第3类算法的优点，提出一种基于点、线相似不变性的自动配准算法。首先，将点云内插成强度影像(和取反的强度影像)，通过SIFT算法提取足够多的匹配点对，并实现点云与影像的粗配准；同时分别基于点云和航空影像提取边缘特征，并通过一定规则将一条较长的边缘自动的分割为若干较短的边缘直线；然后利用局部区域中点特征与线特征的相似不变性，通过匹配点对搜索匹配直线对；最后采用基于扩展共线方程的2D-3D严密配准模型实现航空影像与LiDAR点云的精配准。技术流程如图 1所示。

图 1 基于点、线相似不变性的城区航空影像与机载LiDAR点云的自动配准技术流程 Fig. 1 The flow chart of automatic registration of urban aerial images with airborne LiDAR points based on line-point similarity invariants

图选项

SIFT算法通过在DOG尺度空间寻找极值点作为特征点，并在邻域计算梯度直方图作为该点的局部特征描述符，它能够较好地克服两幅影像因旋转、比例尺不一以及光源不同导致的差异^[12]。本文采用文献^[12]的建议，分别利用LiDAR点云内插距离影像、强度影像和取反后的强度影像与遥感图像进行匹配，提取足够多的匹配点，并利用RANSAC算法消除错配，最后对所有匹配点按可靠性进行排序，选择可靠性较高的匹配点，基于DLT模型进行粗配准。

本文采用从粗到细的策略从LiDAR点云中提取线特征(主要指建筑物轮廓线)，首先利用文献^[18]提出的技术流程从LiDAR点云中快速提取建筑的大致轮廓线并进行分割，然后根据建筑物的扫描特征，提出一种精确提取建筑物特征线的“微元特征分析法”。详细步骤如下：

(1) 建筑物点检测。将点云栅格化成距离影像DSM，再采用形态学滤波方法^[19]获取DEM数据，进一步求得规则化的nDSM，最后利用高差阈值对nDSM进行分割，求出建筑物区域及其对应的建筑物点。

(2) 建筑物轮廓生成与规则化。采用Alpha-Shapes^[20]算法从检测出的建筑物点中提取建筑物初始轮廓，采用袖子算法^[21]提取原始轮廓的关键点，进而得到轮廓线的骨架，再对骨架线进行规则化处理，最后得到规则化后的建筑物轮廓。

(3) 特征直线提取。获取规则化的建筑物轮廓后，依据线段的长度与角度，将轮廓拆分为多条直线。

(4) 利用“微元特征分析法”提取精确直线特征。

由于目前主流的机载LiDAR硬件所采集的点云空间分辨率不高，并非所用边缘点都正好落在屋顶边缘，同时考虑到建筑物点识别(特别是边缘点)可能存在误差等原因，初步提取的建筑物轮廓通常不是严格分布在实际建筑物边界，规则化后的边缘线可能会向建筑物内部收缩不超过点间距0.5倍的距离^[22]。由于LiDAR的扫描特性，航线两侧特征线偏移的方向不一致，内缩误差在平差运算中可能得到一定的消减，文献^[22]提出为控制点设定相应的精度(权值)，计点云高程精度为控制点高程精度，计0.5个点间距与点云平面精度之和为平面精度，然后代入平差方程求最优解。该方法难以完全消除建筑物边缘的内缩误差。

本文利用LiDAR的扫描特性，通过“微元特征分析法”对精化一部分边缘线。其理论为：当扫描角大于一定的阈值时(如图 2(a)所示)，激光光束打到建筑物墙立面形成点云，利用概略边缘线及其在地面上的投影组成一个特征立面V，该特征平面沿其法线方向的微小移动即为V在该特征上的微小变量(本文dv取点间距的0.6倍，如图 2(b))，定义为“微元特征”。统计“微元特征”内点的个数是否达到给定的阈值，若满足，则利用最小二乘原理拟合墙立面，然后以墙立面的水平位置为约束获取精化的边缘直线。具体步骤如下：

图 2 基于LiDAR点云的直线特征提取 Fig. 2 Linear features extraction based on LiDAR point clouds

图选项

(1) 遍历初始边缘线，搜索每条直线的“微元特征”。

(2) 统计微小立面的点云个数，并计算出微元内的点密度。

(3) 若点密度大于阈值，根据最小二乘原理拟合墙立面，以墙立面与建筑物顶面的交线为约束条件精化边缘线；反之，转到(1)。

(4) 遍历所有的边缘线，提取精化后的直线为特征线。

建筑物是否存在立面点、立面点数量的多少与建筑物的高度、朝向、飞机的高度、飞行方向等众多因素相关。LiDAR工程中，视场角一般在50°左右，即扫描角在-25°和+25°之间。试验发现，当扫描角绝对值大于5°时，除特别小、密集的建筑(如棚屋区)外，大部分建筑物至少有一个墙立面存在点云。因此，在城市地区利用“微元特征分析法”得到的精确特征线，可以覆盖大部分区域，满足精配准的要求。当精确提取的特征线分布不足时，可参考文献^[21]的方案，将其他特征线赋以较小的权重，代入平差方程获得最优解。

基于点云提取的线特征精度稍差但完整性好，基于影像提取的边缘能准确覆盖建筑物内外边缘像素，但是比较破碎。本文采用LSD法提取特征直线^[23]，与传统的Canny算子结合Hough变换提取直线的方法相比，漏检现象少，完整性也较好。同时，为了提高自动匹配算法的运行效率，通过计算线段的长度、长度比、角度、端点间距离、线段间法向距离等参数，对基于影像提取的特征线进行优化。具体可分下述几种情况：

(1) 如图 3(a)所示的轮廓ABCD，通过角度阈值和长度阈值，切割成AB、BC、CD 3条线段。

图 3 基于影像提取的直线特征 Fig. 3 The linear features extracted by images

图选项

(2) 如图 3(b)的轮廓ABCDE，简化成选段AE。

(3) 若一条明显的线段被分割成几条短小线段，如图 3(c)所示，线段AB、CD、EF可以合并为线段AF。

(4) 若多条线段平行且线段间距离非常相近，如图 3(d)所示，选段AB、CD可以合并为线段A′B′。

为了满足共面条件，将基于LiDAR点云提取的特征线沿着Z轴投影到点云强度影像上，利用L₁={q₁,q₂,…,q_n}表示其集合，通过L₂={p₁,p₂,…,p_m}表示基于粗匹配后的航空影像提取的直线段集合，C={c₁,c₂,…,c_k}表示匹配点对的集合，其中c_m={(a_i,b_j),0≤m≤k}。本文中的直线通过系数矩阵表示，如直线Ax+By+C=0，表示为[ABC]^T。

如图 4所示，特征线q位于点云空间，特征线p是q在航空影像上的匹配直线，由于粗匹配后的航空影像已被近似地投影到点云空间，大部分误差已消除，同时一组匹配特征线对所涉及的范围较小，可以假设它们之间只存在微小的角度差异和尺度差异，线对之间的关系可以通过式(1)表示

图 4 局部范围点线之间的对应关系 Fig. 4 The corresponding relation between points and lines in local area

图选项

点击浏览公式

如图 4所示，特征点a₁、a₂位于q附近，a₁、a₂在航空影像上的匹配点b₁、b₂位于p周围。可知匹配点对之间的关系也满足用同一变换矩阵，用表示a₁、a₂、b₁、b₂的齐次坐标，可得

点击浏览公式

假设

点击浏览公式

点击浏览公式

点击浏览公式

式(5)的物理意义是：当两个特征点距离特征直线足够近时，特征点与特征直线之间的距离比相等，即满足相似不变性。利用该特性，即可通过匹配点对搜索匹配的直线对。本文以点云空间的特征线q为例，介绍自动匹配的步骤：

(1) 以特征线q的法线方向α·L(q)和与直线平行方向的β·L(q)为界设置搜索区域(如图 5所示，本文α、β各取2)，找出搜索区域内所有的匹配点对，以及影像空间中与该区域有交集的全部“待匹配特征线”ΔP={Δp₁,Δp₂,…,Δp_n}。

图 5 搜索范围 Fig. 5 Search range

图选项

(2) 通过式(6)计算q与每条“待匹配直线”的相似性Sim(q,Δp)，由式(6)可知，当Δp与q相匹配且参与计算的两个匹配点对也没有错配时，Sim(q,Δp)≈1。

点击浏览公式

(3) 为了尽可能地消除粗匹配中错匹配点的影响，本文针对每条“疑似匹配直线”Δp，与搜索区域内所有的匹配点对一起，计算所有的相似性测度值，然后取最大值max Sim(q)作为相似性测度依据

点击浏览公式

(4) 重复(1)—(3)，直到搜索出所有的匹配直线对。

(5) 本文提取建筑物边缘直线为线配准基元，因此有大量的平行线段存在，假设c₁、c₂为平行的特征线附近的匹配点对，当c₁、c₂组成的“线段”与线特征平行时，仅仅通过max Sim(q)进行判断逻辑不够严密，容易造成错判。为了增强算法的抗差性，引入直线之间的距离作为辅助的相似性测度参数。直线距离定义如图 6所示：q为点云空间的特征线；Δp为影像空间的特征线在点云空间的投影直线；O为Δp的中点，取O一定距离(本文取2.5m)的端点A、B，A到q的距离为d₁，B到q的距离为d₂，则q与Δp的距离为

点击浏览公式

图 6 直线段之间的距离 Fig. 6 Distance between two straight line-segments

图选项

(6) 引入长度比和角度差作为测度依据。理论上，同名特征线长度基本一致，考虑到基于影像提取的直线完整性不如基于点云提取的直线，本文限定点云特征线与影像特征线长度比是1～3倍；特征线之间的夹角也可以起到有效的滤除效果，通过设定夹角阈值15°，可快速去除大部分不匹配的直线段。在(3)之前，通过(5)、(6)进行明显的非匹配直线滤除，可有效提高匹配效率。

共线方程是摄影测量中的基本方程，其实质是二维影像坐标和三维物方坐标之间的转化关系。如果在影像空间特征线p上任取一点α，可以确定α在LiDAR点云空间的同名点β在p的匹配特征线q上，通过特征线端点A、B的坐标并引入辅助参数λ可精确表示β的坐标，如式(9)所示

点击浏览公式

将上述关系代入原始的共线方程，得到扩展的共线方程模型(10)，通过该模型实现航空影像与LiDAR点云的精配准

点击浏览公式

通过最小二乘平差迭代使得点云空间每条直线表示的空间点不断接近真值，最后同时解算出航空影像外方位元素以及辅助参数λ，实现精确配准，当初始外方位元素未知时，可利用该思路扩展多项式模型等其他配准模型实现精配准。

为了验证本文算法的实际效果，采用长春某地区的航空影像与LiDAR点云试验(如图 7所示)，测区人工地物丰富，分布合理，典型性较高。其中，LiDAR点云平均点间距约0.7m，高程精度约0.15m，水平精度约0.4m。相机的CCD像元大小为0.0068mm，焦距60mm，航高约1500m，地面分辨率约0.2m。其中，粗配准的结果如表 1所示，基于点云提取的匹配点和特征线(精化后)的分布如图 8所示，粗配准后的点云与影像叠加效果如图 9所示。

图 7 试验数据 Fig. 7 Experimental data

图选项

图 8 基于点云提取的匹配点以及特征线的分布 Fig. 8 The distribution of Match points and feature lines extracted by LiDAR point clouds

图选项

图 9 粗配准后的点云与航空影像叠加(局部) Fig. 9 The superposition of aerial images with LiDAR point clouds in local area

图选项

由表 1可知：①由于距离影像与航空影像纹理差异太大，错匹配率过高，匹配结果不具备实用性；②点云强度影像能取得较好的匹配效果，可能是机载LiDAR所采用的激光属于红外波段而光学影像属于可见光波段的缘故，将强度值取反再与光学影像进行匹配，获取的匹配点数量更多、可靠性更高；③尽管试验中采取RANSAC算法进行了去噪，但无论采用哪种数据，都难以避免较大的错配率。

由图 8、图 9可知：①特征线分布比较均匀，分布密度与人工建筑的分布呈近似的正相关关系；②每条特征线周围有一定数量的匹配点集，基本保证相似不变性的判断条件；③粗匹配后的影像已经近似投影到点云空间，但是通过图 9发现，屋顶面没有完全套合，尚存在一定的偏差。

为了对比精配准和粗配准之间的差别，分别截取几个不同角落的点云，与通过精配准模型转化后的正射影像叠加，如图 10所示。上边是基于SIFT算法的粗配准结果，下面为精细配准后的结果，从图中可以看出，上面的点云与影像位置存在着较明显的错位，精配准后影像与点云的套合度均得到一定程度的改善。

图 10 精配准前后的效果对比 Fig. 10 Comparison before and after precise registration

图选项

由于点云的特殊性，传统光学影像之间的配准精度评价方法(相关系数、互信息等)，不适用于影像与LiDAR点云的配准精度的精确评价。本文利用线特征作为配准精度的评价依据。试验中，分别在点云和影像的9个角落人工勾画9条同名直线作为“检查线”，然后将影像空间的“检查线”投影到点云空间，利用式(8)计算两者的距离作为残差值。表 2列出了粗配准与精配准各个检查线的残差以及残差平均值与中误差。

由表 2可得，精配准后平均误差、中误差分别由粗配准后的0.92m、0.41m、0.23m减少到0.42m、0.23m、0.11m，配准精度得到了较明显的提高。鉴于试验LiDAR数据的点间距约为0.7m，高程精度约0.15m，影像的地面分辨率约为0.2m，可知本文提出的配准算法已达到较理想的精度。

本文提出的自动配准算法的创新点有：① 提出“微元特征分析法”，利用LiDAR的扫描特征，解决基于LiDAR点云提取建筑物顶边缘精确度不足的问题；② 通过局部区域中点特征与线特征的相似不变性原理，利用粗匹配的点对搜索匹配直线对，有效结合基于点特征的粗匹配和基于线特征的精匹配各自的优点，保证配准精度的同时提升算法自动化程度。

试验证明，该方法在城市地区能够以较高的自动化程度取得较理想的配准结果。在直线特征稀少的山区、郊区，如何结合点特征和曲线特征实现航空影像与LiDAR点云的自动配准，需要进一步研究。另外，不同的影像分辨率和点云密度会影响特征提取的精度和密度，进而影响自动匹配的效果，两者的尺度关系对配准精度的影响，需要更丰富的试验去验证。