CN106485690A - 基于点特征的点云数据与光学影像的自动配准融合方法 - Google Patents

Abstract

本发明公开了一种基于点特征的点云数据与光学影像的自动配准融合方法，包括：对点云数据进行滤波处理；分别采用自适应支撑权重稠密立体算法和Delaunay三角剖分算法确定光学影像和点云数据的深度图；通过尺度不变特征变换算法得到光学影像和点云数据的深度图之间的二维匹配关系；通过两步RANSAC算法剔除错误匹配点对，并得到相机位置参数估计；进行点云数据和光学影像的颜色纹理映射，得到融合后的三维图像。本发明无需GPS/INS等初始位置先验，不依赖于场景中人造建筑物的强特征，自动化程度高，鲁棒性好，同时具有较高的配准正确率。

Description

基于点特征的点云数据与光学影像的自动配准融合方法

技术领域

本发明涉及一种改进的异源图像配准算法，特别是一种基于点特征的点云数据与光学影像的自动配准融合方法。

背景技术

机载激光雷达(Airborne Light Detection and Ranging，LiDAR)是将两种传感器(激光扫描仪与光学投影仪)与其他的测量仪器集于一个航空工作平台，通过主动向地面发射激光，接收反射回来的剩余激光来测量地面点的位置，已经成为航空领域内具有高度应用价值的遥感技术。机载LiDAR作为一种主动式遥感技术，其点云数据直接获得地面的三维坐标数据，缺少光谱信息。数字光学影像作为机载LiDAR的另一种数据来源，具有较高的空间分辨率，地物影像几何特征连续，物理特征明显，尤其地物的颜色信息突出。因此，结合激光雷达点云和影像的方法可综合影像的平面精度和点云的高程精度提高建筑物模型的重建精度，弥补单一数据的缺陷，更大程度上使用户获取测区图像足够的信息量，满足实际应用的需要。

目前国内外学者针对机载激光雷达数据和光学影像数据的配准做了大量的研究工作，主要可以分为两类：第一类是基于特征的二维与三维之间的配准，第二类是基于三维重构的三维与三维之间的配准。针对第一种方法，主要是选取点特征、线特征等实现二维光学影像与激光雷达点云数据的直接配准继而得到相机参数估计；Oberkampf等人提出了通过检测多个共面点对相机参数(位置、方向)进行估计，Lee等人在LiDAR点云数据和光学影像中同时根据直线检测影灭点并结合改进的位置估计方法得到相机的平移和旋转参数，Ding等人利用影灭点提取二维正交特征点(2DOCs)完成相机初始位置估计并利用M估计RANSAC方法得到精确的相机位置估计，随后又提出3CS特征用于匹配，此外Christy等人提出了基于线特征检测的配准方法。Liu等人和Kaminsky等通过运动恢复结构(SFM)算法对序列图像生成离散三维点云，继而实现三维到三维点云的配准，Zhao等人利用立体视觉技术通过视频序列进行三维重构，然后利用ICP算法实现三维点云之间的配准。

以上方法普遍需要依赖GPS/INS等初始位置先验，当复杂多变的场景中缺少建筑物的点、线特征时，算法有效性将会大大降低；对于三维之间的配准，算法复杂度很高， 而且依赖于三维重构的精度。

发明内容

本发明的目的在于提供一种基于点特征的点云数据与光学影像的自动配准融合方法。

实现本发明目的的技术解决方案为：

一种基于点特征的点云数据与光学影像的自动配准融合方法，包括以下步骤：

步骤1，对点云数据进行滤波处理；

步骤2，分别采用自适应支撑权重稠密立体算法和Delaunay三角剖分算法确定光学影像和点云数据的深度图；

步骤3，通过尺度不变特征变换算法得到光学影像和点云数据的深度图之间的二维匹配关系；

步骤4，通过两步RANSAC算法剔除错误匹配点对，得到相机位置参数估计；

步骤5，进行点云数据和光学影像的颜色纹理映射，得到融合后的三维图像。

本发明与现有技术相比，其显著优点为：本发明的基于点特征的点云数据与光学影像的自动配准融合方法无需GPS/INS等先验知识，不依赖于场景中人造建筑物的强特征(线、边缘等)、自动化程度高，鲁棒性好，同时具有较高的配准正确率。

附图说明

图1为本发明基于点特征的点云数据与光学影像的自动配准融合方法流程图。

图2(a)为本发明实施例的点云数据，图2(b)为可见光图像，图2(c)-图2(f)为融合后的三维图像。

具体实施方式

假设摄像机内参(焦距、主点坐标、畸变系数)已知，摄像机模型可表示如下：

x＝PX＝K(R|T)X

其中，x是图像坐标，X是三维空间下的世界坐标，P表示摄像机的投影矩阵，K和[R T]分别代表摄像机的内参和外参；R是摄像机的3×3旋转矩阵，T是摄像机的3×1平移矩阵，α，β和γ分别表示偏航、俯仰和横滚的欧拉角，则旋转矩阵R表示如下：

其中，c和s分别代表cos和sin，摄像机内参矩阵K表示如下：

其中，f_x和f_y分别表示摄像机x和y方向的焦距，(x₀，y₀)表示主点坐标，t代表畸变参数；

对极几何是两幅图像之间的内在的射影几何，它独立于景物结构，只依赖于摄像机内参数和相对位置姿态；本质上，两幅图像之间的对极几何是图像平面与以基线(基线是连接两摄像机中心的直线)为轴的平面束的交的几何。

空间点X在两幅图像中成的像分别为x和x′，则空间点X、像点x、像点x′以及摄像机中心共面，这种包含基线的平面成为对极平面，它是由基线和空间点所确定的，并且是一个单参数簇。所有的对极平面都包含基线，因此它们与图像平面的交点必是基线与像平面的交点，称为对极点。

空间点X、像点x以及光心C形成的直线L在第二幅图像中被投影成直线l′，因此，直线L上的空间点X在第二幅图像上的投影x′必然在直线l′上，这条直线称为点x在第二幅图像上的对极线，它是由像点x和基线所确定的对极平面与图像的交线。一张对极平面与左或右像平面交与对极线，并定义了对极线之间的对应关系，以上称为对极几何约束。

对极几何可用一个3*3基本矩阵描述,所以恢复对极几何的问题就转化为估算基本矩阵，基本矩阵是对极几何的代数表示。当给定一对图像，对于左幅图像上的像点x，在右幅图像中存在一条对应的对极线l′。在第二幅图像上，与像点x匹配的点x′必然在对极线l′上。因此，存在一个从左幅图像上的点到右幅图像上与之对应的对极线之间的映射关系：x→l′；

空间中任一不通过摄像机中心的平面π，经过左摄像机中心点C与像点x的射线与平 面π相交于空间点X，空间点X再投影到右摄像机的像点为x′。因为空间点X在过x与C形成的射线上，所有投影点x′必然对应于这条射线的投影即对极线l′。对极线l′、极点e′以及像点x′的关系为：

l′＝e′×x′＝[e′]_×x′

式中，[e′]_×表示3×3的斜对称矩阵

因为存在一个2D的映射，可得x′＝H_πx：

l′＝[e′]_×H_πx＝Fx

F＝[e′]_×H_π，F为基础矩阵，H_π为平面π关于两幅图像之间的单应矩阵；因为[e′]_x的秩为2，而H_π的秩为3，所以F的秩为2；从几何上来讲，F表示由左幅图像的2维射影平面IP²到通过对极点e′的对极线束的映射。

因此，要实现机载LiDAR点云数据和光学影像的自动配准的目的就是寻求最优的摄像机位置估计(内参、旋转矩阵、平移矩阵)，使其投影误差最小化。

下面结合附图对本明进行详细说明。

结合图1，本发明的基于点特征的点云数据与光学影像的自动配准融合方法，包括以下步骤：

步骤1、采用数学形态学滤波方法对点云数据进行滤波处理：

形态学滤波作为一种自下而上的基于局部高差突变思想的滤波器，目的去剔除噪声点、冗余信息以及高程粗差；腐蚀和膨胀运算是形态学图像处理的基础，用于“减少”(腐蚀)或“增大”(膨胀)图像中特征形状的尺寸；基于灰度图像的腐蚀和膨胀运算是在结构元素定义的邻域内选择图像像素值相作用后的最小或最大像素值。

设LiDAR点云数据观测值序列为p(x，y，z)，则在x和y处对高度z的膨胀操作定义为：

式中，(x_p，y_p，z_p)代表p点的邻域窗口w内的点，w又称为结构元素尺寸；邻域窗口可以是一维的直线，也可以使二维的矩形或其他形状；膨胀运算的结果是邻域窗口中的最大高程值；

腐蚀作为膨胀的对偶算子，定义为：

将膨胀和腐蚀进行组合，就得到可直接用于LiDAR滤波的开运算和闭运算，开运算时对数据线进行腐蚀，再进行膨胀，而闭运算正好相反；

采用线性结构元素对LiDAR点云数据的某一纵剖面进行开运算滤波后，若树木目标尺寸小于结构单元，则经过腐蚀运算后可得以去除，大型建筑物则在膨胀运算中得到重建。

步骤2、分别采用自适应支撑权重稠密立体算法和Delaunay三角剖分算法确定光学影像和点云数据的深度图，具体为：

三维点云数据中的任一离散数据点表示为P′(x，y，z)，首先通过Delaunay三角化算法构造LiDAR点云数据的网格化，然后根据高度值z进行灰度级的颜色映射将网格化的LiDAR数据转换成深度图D_P；

对于光学影像I，采用自适应支撑权重稠密立体匹配算法确定光学影像I与具有重叠区域的相邻影像I′之间最优的立体匹配；

待匹配像素点p和q之间的支撑权重表示为：

w(p，q)＝f(Δc_pq，Δd_pq)

＝f(Δc_pq)·f(Δd_pq)Δc_pq和Δd_pq分别表示待匹配像素点p和q之间的颜色差异和欧氏距离差异，f(·)表示差异强度，分别表示如下：

τ_c，τ_d分别为固定常数，实验中取τ_c＝5，τ_d＝17.5。

由于左右图像在拍摄过程中受光照不同的影响，直接进行匹配代价运算会影响确定好的匹配权重关系，因此对匹配的点先进行rank变换，消除两幅图像因光照不同而产生相应的匹配错误，然后再进行匹配代价计算，具体的rank变换表达式如下：

其中I_p为当前匹配像素点的灰度值，I_q为窗口内像素点的灰度值，τ₁τ₂为rank变换分类条件，实际中为固定常数，本实施方式中取τ₁＝9，τ₂＝2；由于窗口内的匹配点进行rank变换后在不同分类等级上为一固定常数，因此在匹配过程中可以消除光照不同对左右图像匹配结果的影响，提高左右图像匹配精度；

光学影像左视图I的匹配点p和光学影像右视图I′的匹配点在匹配窗口内的对应点为q和其匹配代价值为表示如下：

N表示匹配窗口大小，表示左右图像在q和处的匹配值：

在左右图像进行搜索匹配时，采用极线约束，P点的初始视差为，

对初始视差进行左右一致性检查、中值滤波后，由视差图结合对极几何估计得到光 学影像的深度图D_I。

步骤3、通过尺度不变特征变换算法得到光学影像和点云数据的深度图之间的二维匹配关系，具体步骤如下：

步骤3-1、形成尺度空间图像：采用唯一线性高斯卷积核构造高斯差分方程，与深度图进行卷积运算，进而得到高斯差分图像：

D(x，y，σ)＝[G(x，y，kσ)-G(x，y，σ)]×D_I(x，y)

＝L(x，y，kσ)-L(x，y，σ)

其中，x、y为图像坐标值，σ为尺度坐标，L(x，y，σ)为在原图像尺度空间下的高斯图像，k为乘积因子，G(x，y，σ)为尺度可变的二维高斯函数：

光学影像的深度图D_I与点云数据的深度图D_P分别通过不同尺度的高斯核函数进行连续卷积和下采样，得到多组不同尺寸的高斯模糊图像；由高斯差分图像和高斯模糊图像分别构建高斯差分金字塔和高斯金字塔；

步骤3-2、高斯差分金字塔和高斯金字塔中极值点的检测与定义：检测点与其同尺度的邻域及相邻尺度共26个点相比较，如果此检测点为极小值或者极大值点，将该点选为候补特征点；

步骤3-3、生成特征点的特征描述符：选取0，1/4π，1/2π，3/4π，5/4π，3/2π，7/4π这8个方向，这8个方向由4对相反的方向组成，每对相反的方向梯度分别表示减少量和增加量；对这4对相反的方向进行向量运算得到每对的绝对值，将8方向的描述减少到4方向；每个特征点的特征描述符由4方向梯度组成；

步骤3-4、通过设置描述子向量门(本实施方式中门限值取0.2)以限剔除候补特征点中因为非线性光照、相机局部饱和等情况造成的某些方向的梯度值过大的边缘相应特征点向量，然后，对设置门限的特征向量再进行一次归一化处理，得到最后的特征点所对应的归一化特征向量；

步骤3-5、采用欧式距离实现深度图D_P与深度图D_I之间的特征点匹配。

步骤4、通过两步RANSAC算法剔除错误匹配点对，得到相机位置参数估计，具体步骤如下：

步骤4-1、根据空间邻近度，假定的匹配特征对被分成不同的组，即一幅图像被分割成不同的窗口，窗口大小设置为4M*4M像素，M为正整数；窗口从图像的左上角从左至右从上到下开始滑动，滑动单位为M，从而进行下一步分割，以使每个组里的匹配特征对数量在一个合适的阈值范围内每个组内的匹配点对可能会有交叉；

步骤4-2、对每个组的匹配特征对进行一步RANSAC算法计算，剔除局外点：在迭代过程中，两个匹配特征对被同时采样计算对应的单应矩阵，其他的匹配特征对如果符合这个矩阵则被认为是局内点，否则为局外点，选为局内点的匹配特征对为基本匹配特征；

步骤4-3、一个面向全局的一步RANSAC算法再次对不同组内的基本匹配特征进行计算；每两个组的基本匹配特征被采样计算对应的单应矩阵，其他的基本匹配特征如果符合这个矩阵则被认为是局内点，否则为局外点，每个迭代的分数用下面的公式计算。假设有j组内点，第k组的内点个数为n_k，每次迭代计算的得分为：

表示第k组的权重，当内点集中于少数几组时，权重较大；当内点分布于所有组时，权重将会下降；i＝1～j，表示组数；

步骤4-4、保留所有被步骤4-3计算后称为局内点的匹配特征对，被称为局外点的匹配特征对则重新采用分数最高的单应矩阵进行再次计算，以保留更多的特征对；遍历结束后，所有的局内点即为我们需要的匹配特征对。

由于在特征点匹配中，映射模型即为从一个平面上的特征点到另外一平面上的特征点的射影关系，反应为射影矩阵；在经过两步RANSAC消除了错误的匹配点，即局外点之后，局内点即被用来计算映射矩阵。

步骤5、基于上述两步RANSAC算法，剔除错误的匹配点对，得到高精度的相位位置参数估计，进行点云数据与光学影像之间的重投影，完成点云数据与光学影像的颜色 纹理映射，得到融合后的三维图像。

下面结合具体实施例对本发明作进一步说明。

实施例

本实施例采用机载LiDAR点云数据和光学影像，主要覆盖了部分中心地区，包含主要建筑物、低矮的丛林等，选取了部分机载LiDAR点云数据与光学影像具有重叠区域的数据，采用基于点特征的点云数据与光学影像的自动配准融合方法进行处理。

图2(a)为点云数据，图2(b)为可见光图像，图2(c)～图2(f)为配准的分解图，可以看出，图像配准精度高，实现了对城市的三维重建，使点云图像拥有了可见光图像的色彩信息。

本发明基于点特征的点云数据与光学影像的自动配准融合方法有效地提高了图像配准的精度，自动化程序高，而且无需人工参与以及GPS/INS等先验，不依赖场景的强特征，非常适合于场景特征不丰富的情况，实现点云数据与光学影像的自动配准融合。

Claims (6)

1.一种基于点特征的点云数据与光学影像的自动配准融合方法，其特征在于，包括以下步骤：

步骤1，对点云数据进行滤波处理；

步骤2，分别采用自适应支撑权重稠密立体算法和Delaunay三角剖分算法确定光学影像和点云数据的深度图；

步骤3，通过尺度不变特征变换算法得到光学影像和点云数据的深度图之间的二维匹配关系；

步骤4，通过两步RANSAC算法剔除错误匹配点对，得到相机位置参数估计；

步骤5，进行点云数据和光学影像的颜色纹理映射，得到融合后的三维图像。

2.根据权利要求1所述的基于点特征的点云数据与光学影像的自动配准融合方法，其特征在于，步骤1中采用数学形态学方法对点云数据进行滤波处理。

3.根据权利要求1所述的基于点特征的点云数据与光学影像的自动配准融合方法，其特征在于，步骤2中采用Delaunay三角剖分算法确定点云数据的深度图，具体过程为：

首先通过Delaunay三角剖分算法构造点云数据的网格化，然后根据高度值z进行灰度级的颜色映射，将网格化后的点云数据转换成深度图D_P。

4.根据权利要求1所述的基于点特征的点云数据与光学影像的自动配准融合方法，其特征在于，步骤2中采用自适应支撑权重稠密立体算法确定光学影像的深度图，即采用自适应支撑权重稠密立体匹配算法确定光学影像I与具有重叠区域的相邻影像I′之间最优的立体匹配，具体过程为：

步骤2-1，待匹配像素点p和q之间的支撑权重表示为，

$\begin{array}{c}w(p,q)=f({\mathrm{\Δc}}_{pq},{\mathrm{\Δd}}_{pq})\\ =f\left({\mathrm{\Δc}}_{pq}\right)\·f\left({\mathrm{\Δd}}_{pq}\right)\end{array}$

Δc_pq和Δd_pq分别表示待匹配像素点p和q之间的颜色差异和欧氏距离差异，f(·)表示差异强度，分别表示如下：

$f\left({\mathrm{\Δc}}_{pq}\right)=\exp (-\frac{{\mathrm{\Δc}}_{pq}}{{\mathrm{\τ}}_c})$

$f\left({\mathrm{\Δd}}_{pq}\right)=\exp (-\frac{{\mathrm{\Δd}}_{pq}}{{\mathrm{\τ}}_d})$

其中τ_c，τ_d分别为固定常数；

光学影像左视图I的匹配点p和光学影像右视图I′的匹配点在匹配窗口内的对应点为q和其匹配代价为表示如下：

$E(p,{\stackrel{\‾}{p}}_d)=\frac{{\mathrm{\Σ}}_{q,{\stackrel{\‾}{q}}_d\⋐N}w(p,q)w({\stackrel{\‾}{p}}_d,{\stackrel{\‾}{q}}_d)e(q,{\stackrel{\‾}{q}}_d)}{{\mathrm{\Σ}}_{q,{\stackrel{\‾}{q}}_d\⋐N}w(p,q)w({\stackrel{\‾}{p}}_d,{\stackrel{\‾}{q}}_d)}$

为左右图像在q和处的匹配值，N表示匹配窗口大小；

步骤2-2，采用极线约束，p点的初始视差为：

$d_p=\arg \max E(p,{\stackrel{\‾}{p}}_d)$

步骤2-3，对初始视差进行左右一致性检查和中值滤波，由视差图结合对极几何估计得到光学影像的深度图D_I。

5.根据权利要求1所述的基于点特征的点云数据与光学影像的自动配准融合方法，其特征在于，步骤3具体为：

步骤3-1、采用唯一线性高斯卷积核构造高斯差分方程，与深度图进行卷积运算，得到高斯差分图像：

$\begin{array}{c}D(x,y,\mathrm{\σ})=\[G(x,y,k\mathrm{\σ})-G(x,y,\mathrm{\σ})\]\×D_I(x,y)\\ =L(x,y,k\mathrm{\σ})-L(x,y,\mathrm{\σ})\end{array}$

其中，x，y为图像坐标值，σ为尺度坐标，L(x，y，σ)为在原图像尺度空间下的高斯图像，k为乘积因子，G(x，y，σ)为尺度可变的二维高斯函数：

$G(x,y,\mathrm{\σ})=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{\frac{-(x^2+y^2)}{2{\mathrm{\σ}}^2}}$

光学影像的深度图D_I与点云数据的深度图D_P分别通过不同尺度的高斯核函数进行连续卷积和下采样，得到多组不同尺寸的高斯模糊图像；由高斯差分图像和高斯模糊图像分别构建高斯差分金字塔和高斯金字塔；

步骤3-2、检测点与其同尺度的邻域及相邻尺度共26个点相比较，如果此检测点为极小值或者极大值点，将该点选为候补特征点；

步骤3-3、生成特征点的特征描述符：选取0,1/4π,1/2π,3/4π,5/4π,3/2π,7/4π这8个方向，8个方向由4对相反的方向组成，每对相反的方向梯度分别表示减少量和增加量；对这4对相反的方向进行向量运算得到每对的绝对值，将8方向的描述减少到4方向；每个特征点的特征描述符由4方向梯度组成；

步骤3-4、通过设置描述子向量门限剔除候补特征点中梯度值异常的边缘相应特征点，返回步骤3-3，对设置门限的特征向量再进行一次归一化处理，得到最后的特征点所对应的归一化特征向量；

步骤3-5、采用欧式距离实现深度图D_P与深度图D_I之间的特征点匹配。

6.根据权利要求1所述的基于点特征的点云数据与光学影像的自动配准融合方法，其特征在于，步骤4具体为：

步骤4-1、根据空间邻近度，假定的匹配特征对被分成不同的组，即一幅图像被分割成不同的窗口，窗口大小设置为4M*4M像素，M为正整数；窗口从图像的左上角从左至右从上到下开始滑动，滑动单位为M，从而进行下一步分割；

步骤4-2、对每个组的匹配特征对进行一步RANSAC算法计算，剔除局外点：在迭代过程中，两个匹配特征对被同时采样计算对应的单应矩阵，其他的匹配特征对如果符合这个矩阵则被认为是局内点，否则为局外点；选为局内点的匹配特征对为基本匹配特征；

步骤4-3、采用一步RANSAC算法再次对不同组内的基本匹配特征进行计算；每两个组的基本匹配特征被采样计算对应的单应矩阵，其他的基本匹配特征如果符合这个矩阵则被认为是局内点，否则为局外点；每个迭代的分数通过下面的公式得到：

假设有j组内点，第k组的内点个数为n_k，每次迭代计算的得分为：

$G={\mathrm{\Σ}}_{i=1}^j{n}_k\·\sqrt{n_k}$

为第k组的权重，i＝1～j，表示组数；

步骤4-4、保留步骤4-3中被认为局内点的匹配特征对，被称为局外点的匹配特征对则采用分数最高的单应矩阵进行再次计算；遍历结束后，所有的局内点即为所需的匹配特征对；

在特征点匹配中，映射模型即为从一个平面上的特征点到另外一个平面上的特征点的射影关系，通过局内点确定映射矩阵，得到相位位置参数估计。