CN103745474A - 基于惯性传感器和摄像机的图像配准方法 - Google Patents

Abstract

本发明提出一种基于惯性传感器和摄像机的图像配准方法。同步控制摄像机和惯性传感器工作，同时获取同一场景的彩色图像、深度图像以及摄像机的三轴运动角速度和三轴运动加速度；利用卡尔曼滤波对三轴角速度和三轴加速度进行滤波，使用欧拉动力学方程式将三轴角速度转换成欧拉角，使用匀加速运动公式将三轴加速度转换成三轴位移，从而获得摄像机外参矩阵；利用彩色图像和深度图像的仿射变换关系，和深度图像灰度值与像素深度信息的比例关系求解出角点在深度图像中的深度信息；使用OFCM模型计算获得图像背景光流；使用RANSAC算法完成图像配准。本发明方法计算量小，实时性和精度均较高，能适用于手持设备和移动设备。

Description

基于惯性传感器和摄像机的图像配准方法

技术领域

本发明属于数字图像处理领域，具体涉及一种基于惯性传感器和摄像机的图像配准方法。

背景技术

图像配准技术的关键问题是对摄像机运动参数的估计，一般通过摄像机标定来完成。现有技术中，摄像机自标定的方法有：

李立春等人提出了基于八点法基础矩阵三维重建算法[李立春,邱志强,王鲲鹏等.基于匹配测度加权求解基础矩阵的三维重建算法[J].计算机应用,2007,27(10):2530-2533]，该方法复杂并且计算量很大，难以实现图像处理的实时性。

张洋等人提出了FOE模型验证了摄像机运动的平移效果[张洋,凯瑟勒威驰,柏森.可见光范围内基于移动相机的运动目标检测算法.计算机视觉和模式识别工作室.IEEE国际会议,2006,130-131]，该方法没有考虑摄像机的旋转运动和摄像机的内参矩阵，精度不高。

郑红等人提出了基于主动视觉的标定方法主要是利用云台等已知运动的硬件平台来提供摄像机运动信息[郑红,刘振强.基于精确模型的云台摄像机自标定[J].机器人，2013.5,35(3):1994-2013]，该方法由于云台等设备较大，虽然能够很好地实现实时性，却很难普及到手持设备和移动设备中，使用面有限。

发明内容

本发明提出一种基于惯性传感器和摄像机的图像配准方法，该方法计算量小，实时性和精度均较高，能适用于手持设备和移动设备。

为了解决上述技术问题，本发明提供一种基于惯性传感器和摄像机的图像配准方法，包括以下步骤：

步骤一、同步控制摄像机和惯性传感器工作，同时获取同一场景的彩色图像I和J、深度图像I’和J’以及摄像机的三轴运动角速度[w_x，w_y，w_z]和三轴运动加速度[a_x，a_y，a_z]；

步骤二、提取彩色图像I的角点m₁(x₁,y₁)；

步骤三、利用卡尔曼滤波对摄像机的三轴角速度[w_x，w_y，w_z]和三轴加速度[a_x，a_y，a_z]进行滤波，获得卡尔曼滤波估计三轴角速度和卡尔曼滤波估计三轴加速度使用欧拉动力学方程式将

转换成欧拉角[α_k，β_k，γ_k]，使用匀加速运动公式将

转换成三轴位移[t_xk,t_yk,t_zk]，从而获得摄像机外参矩阵[R,T]中的旋转矩阵R和平移矩阵T；

步骤四、利用彩色图像和深度图像的仿射变换关系，求解出彩色图像I的角点m₁(x₁,y₁)在深度图像I’和J’中的坐标位置m_1ir(x₁,y₁)和m_2ir(x₁,y₁)；利用深度图像灰度值与像素深度信息的比例关系求解出角点m₁(x₁,y₁)在深度图像I’和J’中的深度信息u₁和u₂；

步骤五、将步骤三获得的旋转矩阵R和平移矩阵T以及步骤四获得的彩色图像的深度信息u₁、u₂带入OFCM模型，计算获得图像J的背景光流m₂(x₂,y₂)，OFCM模型如公式（1）所示，

$m_2=\frac{u_1}{u_2}{\mathrm{KRK}}^{-1}{m}_1+\frac{1}{u_2}\mathrm{KT}---\left(1\right)$

式（1）中，K为摄像机的内参矩阵，K^-1为K的逆矩阵；

步骤六、使用RANSAC算法对彩色图像I的背景光流m₂进行鲁棒性估计，将背景光流m₂中存在的由运动物体产生的光流和误配准光流剔除掉，从而完成彩色图像I和彩色图像J的图像配准，然后利用保留的彩色图像I的正确背景光流对彩色图像J进行初始化，进行后续图像的配准。

主要发明原理：

1、本发明结合欧式空间的刚体变化关系提出了图像背景光流与摄像机运动参数的关系模型(OFCM)，该模型主要推导方式如下：

假设彩色CCD为针孔模型，场景中的背景只存在与彩色CCD的相对运动，摄像机坐标系的建立是以彩色CCD的光心为中心O，光轴所在的直线为Z轴，再根据左手原则建立X，Y轴。

假设摄像机运动前摄像机坐标系与世界坐标系重合，那么空间背景中的任意一点P的坐标满足：

P_W＝P_C         (2)

式（2）中，P_W为空间某点的世界坐标系坐标，P_W为空间某点的摄像机坐标系坐标，空间背景点P与其成像平面上的成像点的齐次坐标m₁(x₁,y₁,1)满足公式(3)：

u₁m₁＝KP_C         (3)

式（3）中，m₁(x₁,y₁)为成像平面上的成像点，利用模型计算时采用Harris角点代替；u₁可近似为m₁(x₁,y₁)对应空间点的深度信息；K为摄像机的内参矩阵，考虑到图像中像素坐标系和物理坐标系的关系，可得出摄像机内参矩阵的五参数模型，如公式（4）所示，

$K=\left[\begin{array}{ccc}{f}_x& s& u_0\\ 0& f_y& v_0\\ 0& 0& 1\end{array}\right]---\left(4\right)$

式（4）中，

其中，dx为像素点在x方向上的物理尺寸，dy为像素点在y方向上的物理尺寸，[u₀,v₀]为主点在像素坐标系中的坐标，s为相机畸变系数。

假设摄像机存在任意方向的平移运动和旋转运动，平移运动可用平移向量T＝[t_x,t_y,t_z]^T来表示，旋转运动可用3*3的旋转矩阵R来表示，此时空间背景点P在摄像机坐标系中的坐标满足：

P'_C＝R×P_W+T＝R×P_C+T         (5)

式（5）中，旋转矩阵R可用欧式空间的欧拉角α,β,γ表示，如式（6）所示，

$R=\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\β}& \cos \mathrm{\γ}\sin \mathrm{\β}\sin \mathrm{\α}-\sin \mathrm{\γ}& \cos \mathrm{\γ}\sin \mathrm{\β}\cos \mathrm{\α}-\sin \mathrm{\α}\sin \mathrm{\γ}\\ \sin \mathrm{\γ}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\γ}\cos \mathrm{\α}& \sin \mathrm{\γ}\sin \mathrm{\β}\cos \mathrm{\α}+\cos \mathrm{\γ}\sin \mathrm{\α}\\ -\sin \mathrm{\β}& \cos \mathrm{\β}\sin \mathrm{\α}& \cos \mathrm{\β}\cos \mathrm{\α}\end{array}\right]---\left(6\right)$

摄像机运动后，空间点P与其成像平面上的成像点的齐次坐标m₂(x₂,y₂,1)满足关系:

u₂m₂＝KP'_C＝K(R×P_C+T)         (7)

根据公式(3)和公式(5)，可以获得摄像机运动前后空间点P在成像平面上的运动光流，即OFCM的最终形式：

$m_2=\frac{u_1}{u_2}{\mathrm{KRK}}^{-1}{m}_1+\frac{1}{u_2}\mathrm{KT}---\left(8\right)$

其中，u₁和u₂可以近似为空间点P在摄像机坐标系中的深度信息，通过摄像机获取的深度图像计算得出。外参矩阵[R,T]通过硬件平台获取并处理得出，m₁可以通过Harris角点检测获取。假定摄像机的深度CCD和彩色CCD均工作在固定焦距f的条件下，此时像平面与物平面是一一对应的。因此，彩色CCD的内参矩阵K在整个图像获取过程中是一个恒定不变的量，可以经过一次标定获取后作为经验矩阵使用。

优势：现今存在的摄像机模型有FOE模型、MOF模型，FOE模型和MOF模型均没有考虑摄像机的内参矩阵，且FOE模型只支持摄像机的平移运动，因此在摄像机复杂运动下无法保持较高的精度。本发明提出的OFCM模型包含摄像机的内参矩阵，支持摄像机的三轴平移运动和三轴旋转运动，能够克服场景亮度突变。因此，OFCM模型在摄像机复杂运动和场景亮度突变情况下，都能保持高精度的图像配准。

2、基于惯性传感器和摄像机的硬件平台：

本发明利用OFCM模型来计算图像的背景光流，由OFCM模型可知，图像背景光流的计算需要获知摄像机的外参矩阵[R,T]和图像像素对应的深度信息u₁、u₂。本发明采用硬件平台来完成摄像机外参矩阵和图像像素深度信息的获取。

本发明所述硬件平台包括三个部分：控制器件FPGA、惯性传感器和摄像机，摄像机可以采用Kinect摄像机，惯性传感器中包括陀螺仪和加速度计，Kinect摄像机包括深度图像获取CCD和彩色图像获取CCD，FPGA和惯性传感器的最小***均固定在深度图像获取CCD上，惯性传感器位于Kinect摄像机的正上方，陀螺仪用于获取摄像机的三轴运动角速度[w_x，w_y，w_z]，加速度计用于获取摄像机的三轴运动加速度[a_x，a_y，a_z]，Kinect摄像机用于拍摄彩色图像I、J和深度图像I’、J’。硬件平台上电后，FPGA通过控制脉冲控制Kinect摄像机的快门工作，从而完成彩色图像和深度图像的获取，FPGA通过SPI协议同步控制Kinect摄像机和惯性传感器工作，以保证本步骤获取的图像序列和采集的摄像机运动参数一一对应。

优势：基于本模型的图像配准算法计算量小，***实时性好；由于惯性传感器具有体积小、精度高等优点，和云台等获取摄像机运动参数的设备相比，能够更好地普及在手持设备和移动设备中。

3、基于惯性传感器的卡尔曼滤波模型的建立：

由于惯性传感器采集到的摄像机运动参数存在零点偏移、刻度因子误差和随机噪声等误差，因此需要对采集到的摄像机运动参数进行滤波去噪处理。

在使用卡尔曼滤波的经典公式时，需要根据惯性传感器的性能与参数建立惯性传感器的误差模型和状态转移关系式：

z′_k=z_k(M+1)+B_f+U         (9)

z_k=Ψ_k,k-1z_k-1+V         (10)

其中，z′_k为k时刻***的实际观测值，z_k为k时刻***的真实值，M为刻度因子误差，B_f为零点漂移，V和U为随机高斯白噪声，Ψ_k,k-1为状态转移矩阵。

根据公式(9)并结合时间序列分析法求解出惯性传感器的刻度因子误差M，零点漂移B_f，和高斯白噪声V和U；根据公式(10)并结合时间序列分析法求解出参数间的状态转移矩阵Ψ_k,k-1。

将公式公式(9)和公式(10)计算的结果带入卡尔曼滤波的经典公式(11)，对k时刻的参数进行滤波，求出估计值

$\left\{\begin{array}{c}{z}_{k,k-1}={\mathrm{\Ψ}}_{k,k-1}{z}_{k-1}\\ \widehat{z_k}=z_{k,k-1}+K_k\left[{z^{\′}}_k(M+1)z_{k,k-1}\right]\\ P_{k,k-1}={\mathrm{\Ψ}}_{k,k-1}{P}_{k-1}{{\mathrm{\Ψ}}_{k,k-1}}^T+Q_{k-1}\\ K_k=P_{k,k-1}{(M+1)}^T{P}_k^{-1}\\ P_k=[I-K_k(M+1)]P_{k,k-1}\end{array}\right.---\left(11\right)$

其中，Q为***的噪声方差，K_k为卡尔曼增益，P_k为状态协方差，w_k,k-1为状态预估值。

本发明与现有技术相比，其显著优点在于：

（1）由于摄像机参数的估计是通过惯性传感器和摄像机实时获取的，其中惯性传感器中的陀螺仪和加速度计分别获取摄像机的旋转运动参数和平移运动参数；摄像机包含深度图像获取CCD(深度CCD)和彩色图像获取CCD(彩色CCD)，二者可以结合获取彩色图像的深度信息。因此本方法虽然基于硬件平台，但由于惯性传感器具有体积小、精度高等优点，与摄像机绑定后，可以很好地普及在手持设备和移动设备中。

（2）从摄影几何角度出发，结合欧式空间的刚体变化关系推导出图像背景光流与摄像机运动参数的关系模型(optical flow and camera’s motion)，记作OFCM模型。此模型包含摄像机的内参矩阵，并且支持摄像机的三轴平移运动和三轴旋转运动；同时，在亮度突变的情况下，仍然有着较高的精度，相较于FOE模型和光流匹配法有了较大的突破。

附图说明

图1为本发明基于惯性传感器和摄像机的图像配准方法流程图。

图2为在场景亮度突变情况下使用本发明方法获得的图像配准效果图。

图3为在摄像机发生平移运动时使用本发明方法获得的图像配准效果图。

图4为在摄像机发生旋转运动时使用本发明方法获得的图像配准效果图。

图5为在摄像机发生平移运动和旋转运动时使用本发明方法获得的图像配准效果图。

具体实施方式

如图1所示，本发明基于惯性传感器和摄像机的图像配准方法，包括以下步骤

步骤一、在硬件平台中，同步控制摄像机和惯性传感器工作，同时获取同一场景的彩色图像I和J、深度图像I’和J’以及摄像机的运动参数，即三轴运动角速度[w_x，w_y，w_z]和三轴运动加速度[a_x，a_y，a_z]，具体实现如下：

本发明所述硬件平台包括三个部分：控制器件FPGA、惯性传感器和摄像机，惯性传感器包括陀螺仪和加速度计，摄像机包括深度图像获取CCD和彩色图像获取CCD，用于拍摄彩色图像I、J和深度图像I’、J’，FPGA和惯性传感器的最小***均固定在深度图像获取CCD上，陀螺仪用于获取摄像机的三轴运动角速度[w_x，w_y，w_z]，加速度计用于获取摄像机的三轴运动加速度[a_x，a_y，a_z]。硬件平台上电后，FPGA通过控制脉冲控制摄像机的快门工作，从而完成彩色图像和深度图像的获取，FPGA通过SPI协议同步控制摄像机和惯性传感器工作，以保证本步骤获取的图像序列和采集的摄像机运动参数一一对应。

步骤二、利用Harris角点检测提取彩色图像I的角点m₁(x₁,y₁)。具体实现如下：

2.1计算彩色图像I中像素点p[x,y]的相关矩阵M，计算方法如公式(1)所示，

$M=\underset{x,y}{\mathrm{\Σ}}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]=w(x,y)\⊗\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]---\left(2\right)$

式(1)中，(x,y)为像素点p的图像坐标，w(x,y)是窗口大小为w*w的高斯权重函数，I_x代表行方向的导数，I_y代表列方向的导数；

2.2计算像素点m₁(x₁,y₁)的角点响应R(x,y)，计算方法如公式(2)所示，

R(x,y)＝detM-k*(traceM)²         (2)

式(2)中，detM代表相关矩阵M行列式的值，traceM代表相关矩阵M的迹，k为常数，一般为0.04-0.06；

2.3判断像素点m₁(x₁,y₁)是否为角点，判断方法为公式(3)所示，若像素点m₁(x₁,y₁)的角点响应R(x,y)满足公式(3)，则像素点m₁(x₁,y₁)为角点；若像素点m₁(x₁,y₁)的角点响应R(x,y)不满足公式(11)，则将m₁(x₁,y₁)舍弃；

[R(x,y)≥R(x′,y′)]∩[R(x,y)≥0.01R_max]         (3)

式(3)中，R(x′,y′)代表窗口w*w中除像素点m₁(x₁,y₁)之外的其它像素点的角点响应，R_max为边缘特征信息图像I_f中的最大角点响应。

步骤三、利用卡尔曼滤波对摄像机的三轴角速度[w_x，w_y，w_z]和三轴加速度[a_x，a_y，a_z]进行滤波，获得卡尔曼滤波估计三轴角速度

和卡尔曼滤波估计三轴加速度使用欧拉动力学方程式将

转换成欧拉角[α_k，β_k，γ_k]，使用匀加速运动公式将

转换成三轴位移[t_xk,t_yk,t_zk]从而获得摄像机外参矩阵[R,T]中的旋转矩阵R和平移矩阵T；

3.1利用卡尔曼滤波对摄像机的运动参数，即三轴角速度[w_x，w_y，w_z]和三轴加速度[a_x，a_y，a_z]进行滤波，获得卡尔曼滤波估计三轴角速度和卡尔曼滤波估计三轴加速度

的具体计算方式如下：

3.1.1根据惯性传感器的性能与运动参数建立惯性传感器的误差模型和状态转移关系式，结合时间序列分析法求解出惯性传感器的刻度因子误差M、零点漂移B_f、高斯白噪声U、运动参数间的状态转移矩阵Ψ_k,k-1以及高斯白噪声V；

误差模型如公式（4）所示，状态转移关系式公式（5）所示，

z′_k=z_k(M+1)+B_f+U         (4)

z_k=Ψ_k,k-1z_k-1+V         (5)

式（4）和式（5）中，z代表三轴角速度[w_x，w_y，w_z]和三轴加速度[a_x，a_y，a_z]中的某一运动参数，z′_k为k时刻惯性传感器测量获得的某一运动参数的测量值，z_k为k时刻某一运动参数的真实值，M为刻度因子误差，B_f为零点漂移，V和U为随机高斯白噪声，Ψ_k,k-1为状态转移矩阵；

根据公式(4)并结合时间序列分析法求解出惯性传感器的刻度因子误差M，零点漂移B_f，和高斯白噪声V和U；根据公式(5)并结合时间序列分析法求解出运动参数间的状态转移矩阵Ψ_k,k-1；

3.1.2将步骤3.1.1所求取的结果带入卡尔曼滤波的经典公式(6)对k时刻的运动参数z_k进行滤波，求出该运动参数的估计值

$\left\{\begin{array}{c}{z}_{k,k-1}={\mathrm{\Ψ}}_{k,k-1}{z}_{k-1}\\ \widehat{z_k}=z_{k,k-1}+K_k\left[{z^{\′}}_k(M+1)z_{k,k-1}\right]\\ P_{k,k-1}={\mathrm{\Ψ}}_{k,k-1}{P}_{k-1}{{\mathrm{\Ψ}}_{k,k-1}}^T+Q_{k-1}\\ K_k=P_{k,k-1}{(M+1)}^T{P}_k^{-1}\\ P_k=[I-K_k(M+1)]P_{k,k-1}\end{array}\right.---\left(6\right)$

式（6）中，Q为***的噪声方差，K_k为卡尔曼增益，P_k为状态协方差，z_k,k-1为状态预估值。

按照上述方法分别对三轴运动角速度[w_x，w_y，w_z]和三轴运动加速度[a_x，a_y，a_z]进行卡尔曼滤波，获得和

3.2使用欧拉动力学方程将k时刻滤波后的三轴运动角速度w_xk，w_yk，w_zk转换成欧式坐标系中的欧拉角α_k，β_k，γ_k，将所得的α_k，β_k，γ_k带入旋转矩阵R中，完成彩色CCD旋转矩阵R的计算；

所述欧拉动力学方程如公式(7)所示，

$\left[\begin{array}{c}{\stackrel{\·}{\mathrm{\α}}}_k\\ {\stackrel{\·}{\mathrm{\β}}}_k\\ {\stackrel{\·}{\mathrm{\γ}}}_k\end{array}\right]=\left[\begin{array}{ccc}-\sin {\mathrm{\α}}_k/\tan {\mathrm{\β}}_k& {\cos \mathrm{\α}}_k/\tan {\mathrm{\β}}_k& 1\\ {\cos \mathrm{\α}}_k& \sin {\mathrm{\α}}_k& 0\\ \sin {\mathrm{\α}}_k/\sin {\mathrm{\β}}_k& -\cos {\mathrm{\α}}_k/\sin {\mathrm{\β}}_k& 0\end{array}\right]*\left[\begin{array}{c}{w}_{\mathrm{xk}}\\ w_{\mathrm{yk}}\\ w_{\mathrm{zk}}\end{array}\right]---\left(7\right)$

式（7）中，为欧拉角α_k，β_k，γ_k的一阶导数，将k-1时刻求得的参量α_k-1，β_k-1，γ_k-1作为公式(7)的初值；

所述旋转矩阵R如公式(8)所示：

$R=\left[\begin{array}{ccc}\cos {\mathrm{\γ}}_k\cos {\mathrm{\β}}_k& \cos {\mathrm{\γ}}_k\sin {\mathrm{\β}}_k\sin {\mathrm{\α}}_k-\sin {\mathrm{\γ}}_k& \cos {\mathrm{\γ}}_k\sin {\mathrm{\β}}_k\cos {\mathrm{\α}}_k-\sin {\mathrm{\α}}_k\sin {\mathrm{\γ}}_k\\ \sin {\mathrm{\γ}}_k\cos {\mathrm{\β}}_k& \sin {\mathrm{\α}}_k\sin {\mathrm{\β}}_k\sin {\mathrm{\γ}}_k+\cos {\mathrm{\γ}}_k\cos {\mathrm{\α}}_k& \sin {\mathrm{\γ}}_k\sin {\mathrm{\β}}_k\cos {\mathrm{\α}}_k+\cos {\mathrm{\γ}}_k\sin {\mathrm{\α}}_k\\ -\sin {\mathrm{\β}}_k& \cos {\mathrm{\β}}_k\sin {\mathrm{\α}}_k& \cos {\mathrm{\β}}_k\cos {\mathrm{\α}}_k\end{array}\right]---\left(8\right)$

3.3使用如公式(9)所示的速度位移公式将三轴加速度a_xk，a_yk，a_zk转换成三轴位移t_xk,t_yk,t_zk。由于摄像机运动参数的采集频率比较高，因此可以将摄像机的平移运动近似为匀加速运动，计算公式如公式(9)所示：

$\left\{\begin{array}{c}{t}_k=v_{k}T+\frac{1}{2}{a}_k{T}^2\\ v_k=v_{k-1}+a_{k-1}\×T\end{array}\right.----\left(9\right)$

式(9)中，t_k代表三轴位移t_xk,t_yk,t_zk中任意一坐标轴方向的位移，T为摄像机的拍摄周期，a_k-1为k-1时刻对应坐标轴方向的加速度，v_k-1和v_k分别为k-1时刻和k时刻的初始速度。本发明给定***初始时刻的初始速度0，即v₀＝0。将三轴加速度a_xk，a_yk，a_zk分别带入公式(9)，可计算获得三轴平移运动的平移向量T＝[t_xk,t_yk,t_zk]。

步骤四、利用彩色图像和深度图像的仿射变换关系，分别求解出彩色图像I的角点m₁(x₁,y₁)在深度图像I’和J’中的坐标位置m_1ir(x₁,y₁)和m_2ir(x₁,y₁)；利用深度图像中灰度值与像素深度信息的比例关系求解出角点m₁(x₁,y₁)在深度图像I’和J’中的深度信息u₁和u₂，具体实现如下：

4.1利用彩色图像和深度图像的仿射变换关系，分别求解出彩色图像I的角点m₁(x₁,y₁)在深度图像I’和J’中的坐标位置m_1ir(x₁,y₁)和m_2ir(x₁,y₁)，具体实现如下：

根据直接线性变换算法计算出彩色图像和深度图像的仿射变换关系，具体如公式(10)所示：

m_ir＝Hm         (10)

式(10)中，m_ir为深度图像中的像素坐标，m为彩色图像中的像素坐标。将彩色图像I的角点m₁(x₁,y₁)带入公式(9)求解出角点m₁(x₁,y₁)在深度图像I’和J’中的坐标位置m_1ir(x₁,y₁)和m_2ir(x₁,y₁)。

4.2利用深度图像中灰度值与像素深度信息的比例关系求解出彩色图像角点m₁(x₁,y₁)在深度图像I’和J’中的深度信息u₁和u₂，具体计算过程如下

灰度值与像素深度信息的比例关系可以通过反向工程提出的经验关系式，如公式(11)所示：

u＝l×m         (11)

式(11)中，m(x,y)为像素坐标(x,y)处的灰度值，u为像素坐标(x,y)处对应的深度值，l为比例系数；将步骤4.2求解出的m_1ir(x₁,y₁)，m_2ir(x₁,y₁)带入公式(11)求解出对应的深度信息u₁和u₂。

步骤五、将步骤三获得的运动外参矩阵R和T以及步骤四获得的彩色图像的深度信息u₁、u₂带入OFCM模型，计算获得图像J的背景光流m₂(x₂,y₂)，具体实现如下：

将步骤三获得的运动外参矩阵R和T以及步骤四获得的彩色图像的深度信息u₁、u₂带入如公式（12）所示的OFCM模型，计算获得第k帧角点m₁在第k+1帧图像中的背景光流m₂，

$m_2=\frac{u_1}{u_2}{\mathrm{KRK}}^{-1}{m}_1+\frac{1}{u_2}\mathrm{KT}---\left(12\right)$

式（12）中，K为摄像机的内参矩阵，K^-1为K的逆矩阵。

步骤六、使用RANSAC算法对彩色图像I的背景光流m₂进行鲁棒性估计，将背景光流m₂中存在的由运动物体产生的光流和误配准光流剔除掉，从而完成彩色图像I和彩色图像J的图像配准，然后利用保留的彩色图像I的正确背景光流对彩色图像J进行初始化，进行后续图像的配准。具体实现如下：

所述彩色图像I的RANSAC计算过程为：

6.1初始化当前所有的角点m₁(x₁,y₁)，m₂(x₂,y₂)为局内点；

6.2利用DLT算法计算出仿射变换矩阵H₁，根据公式（13）结合角点m₁(x₁,_y1)估计出角点估计值m₂′(x₂′,y₂′)：

m₂′=H₁m₁         (13)

6.3根据公式(14)计算角点m₂(x₂,y₂)和角点估计值m₂′(x₂′,y₂′)的欧氏距离d，若d小于阈值d′，则将角点m₁(x₁,y₁)，m₂(x₂,y₂)置为局内点(即背景光流)，若d大于等于阈值d′，则将角点m₁(x₁,y₁)，m₂(x₂,y₂)置为局外点，局外点即为运动物体光流或误配准光流。当完成一次RANSAC计算后，初始化当前局内点为局外点，重复步骤6.2和6.3，当不再产生新的局内点或达到规定重复次数时，结束当前RANSAC算法，

d＝sqrt((x₂-x₂′)²+(y₂-y₂′)²)         (14)

所述对彩色图像J进行初始化的过程为：首先遍历该帧所有的角点，若该角点与彩色图像I保留的正确背景光流重合，则将该角点作为RANSAC计算的初始局内点。重复步骤6.1、6.2、6.3，逐帧完成RANSAC的估计。

根据上述步骤，不断地将图像中的运动物体光流或误配准光流剔除掉；利用RANSAC算法提取的角点光流，完成图像的配准。

本发明的有益效果可以通过以下实验进一步说明：

本实验采用的主控制器件为ALTERA公司的型号为EP2C5Q208C6的FPGA，传感器包含两部分：型号为ADIS16405的惯性传感器和Kinect摄像机，其中惯性传感器中的陀螺仪和加速度计分别获取摄像机的旋转运动参数和平移运动参数；Kinect摄像机包含深度图像获取CCD(深度CCD)和彩色图像获取CCD(彩色CCD)，二者可以结合获取彩色图像的深度信息。由于惯性传感器的尺寸为35mm*30mm，所以硬件平台体积小，可以普及在移动设备中。

利用Kinect相机采集得到大小均为640*480的彩色图像和深度图像。根据对PC端算法的实时性分析，本实施例设置硬件平台的图像采集频率和惯性传感器的参数采集频率均为40帧/s。Harris角点检测的窗口选取为5*5，RANSAC算法的计算次数设置为30次。

为了更好地验证本发明方法提出的OFCM模型的效果，本实验分别验证了Kinect摄像机在亮度突变、随机平移运动、随机旋转运动、随机平移旋转复杂运动条件下，使用本发明方对同一场景进行图像配配准，分别获得如图2、图3、图4、图5所示的效果图。在图2、图3、图4、图5中，每幅图的左半部分为左图、右半部分为右图，左右图间隔均为10帧，其中左图为第k帧图像，右图为第k+10帧图像。从图2、图3、图4、图5看出，本发明方法能够精准的完成图像的配准。

表1为上述不同摄像机运动条件下，PC端图像配准性能统计分析。根据表1，对同种运动方式的各项指标进行纵向比较，可以发现本发明方法在Kinect摄像机发生平移运动和旋转运动时，图像配准都取得了较好的结果，平均处理精度可以达到69.5%；对不同运动方式下图像配准的情况的各项指标进行横向比较，可以发现本发明方法对Kinect摄像机单一运动的配准精度高于混合运动的配准精度，此误差主要是由陀螺仪和加速度计的误差叠加造成的；就处理实时性来看，在PC机上的平均处理时间约22.92ms，时间上的差别主要是Kinect摄像机混合运动产生误匹配点的概率变大，因此RANSAC算法估计相对耗时；根据PC端的时间估计，综合考虑硬件采集的实际难度，设置Kinect摄像机的图像获取帧频和惯性传感器的参数的同步采集频率为30Hz。

表1PC端图像配准性能统计表

Claims (8)

1.基于惯性传感器和摄像机的图像配准方法，其特征在于，包括以下步骤：

步骤一、同步控制摄像机和惯性传感器工作，同时获取同一场景的彩色图像I和J、深度图像I’和J’以及摄像机的三轴运动角速度[w_x，w_y，w_z]和三轴运动加速度[a_x，a_y，a_z]；

步骤二、提取彩色图像I的角点m₁(x₁,y₁)；

步骤三、利用卡尔曼滤波对摄像机的三轴角速度[w_x，w_y，w_z]和三轴加速度[a_x，a_y，a_z]进行滤波，获得卡尔曼滤波估计三轴角速度

和卡尔曼滤波估计三轴加速度；使用欧拉动力学方程式将

转换成欧拉角[α_k，β_k，γ_k]，使用匀加速运动公式将

转换成三轴位移[t_xk,t_yk,t_zk]，从而获得摄像机外参矩阵[R,T]中的旋转矩阵R和平移矩阵T；

步骤四、利用彩色图像和深度图像的仿射变换关系，求解出彩色图像I的角点m₁(x₁,y₁)在深度图像I’和J’中的坐标位置m_1ir(x₁,y₁)和m_2ir(x₁,y₁)；利用深度图像灰度值与像素深度信息的比例关系求解出角点m₁(x₁,y₁)在深度图像I’和J’中的深度信息u₁和u₂；

步骤五、将步骤三获得的旋转矩阵R和平移矩阵T以及步骤四获得的彩色图像的深度信息u₁、u₂带入OFCM模型，计算获得图像J的背景光流m₂(x₂,y₂)，OFCM模型如公式（1）所示，

$m_2=\frac{u_1}{u_2}{\mathrm{KRK}}^{-1}{m}_1+\frac{1}{u_2}\mathrm{KT}---\left(1\right)$

式（1）中，K为摄像机的内参矩阵，K^-1为K的逆矩阵；

步骤六、使用RANSAC算法对彩色图像I的背景光流m₂进行鲁棒性估计，将背景光流m₂中存在的由运动物体产生的光流和误配准光流剔除掉，从而完成彩色图像I和彩色图像J的图像配准，然后利用保留的彩色图像I的正确背景光流对彩色图像J进行初始化，进行后续图像的配准。

2.如权利要求1所述的基于惯性传感器和摄像机的图像配准方法，其特征在于，步骤一中，FPGA通过SPI协议同步控制摄像机和惯性传感器，惯性传感器包括陀螺仪和加速度计，摄像机包括深度图像获取CCD和彩色图像获取CCD，FPGA和惯性传感器的最小***均固定在深度图像获取CCD上。

3.如权利要求1所述的基于惯性传感器和摄像机的图像配准方法，其特征在于，步骤二中使用Harris角点检测提取彩色图像I的角点m₁(x₁,y₁)，计算过程如下：

3.1计算彩色图像I中像素点p[x,y]的相关矩阵M，计算方法如公式(2)所示，

$M=\underset{x,y}{\mathrm{\Σ}}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]=w(x,y)\⊗\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]---\left(2\right)$

式(2)中，(x,y)为像素点p的图像坐标，w(x,y)是窗口大小为w*w的高斯权重函数，I_x代表行方向的导数，I_y代表列方向的导数；

3.2计算像素点m₁(x₁,y₁)的角点响应R(x,y)，计算方法如公式(3)所示，

R(x,y)＝detM-k*(traceM)²           (3)

式(3)中，detM代表相关矩阵M行列式的值，traceM代表相关矩阵M的迹，k为常数，一般为0.04-0.06；

3.3判断像素点m₁(x₁,y₁)是否为角点，判断方法为公式(4)所示，若像素点m₁(x₁,y₁)的角点响应R(x,y)满足公式(4)，则像素点m₁(x₁,y₁)为角点；若像素点m₁(x₁,y₁)的角点响应R(x,y)不满足公式(4)，则将m₁(x₁,y₁)舍弃；

[R(x,y)≥R(x',y')]∩[R(x,y)≥0.01R_max]        (3)

式(4)中，R(x',y')代表窗口w*w中除像素点m₁(x₁,y₁)之外的其它像素点的角点响应，R_max为边缘特征信息图像I_f中的最大角点响应。

4.如权利要求1所述的基于惯性传感器和摄像机的图像配准方法，其特征在于，步骤三中，所述获得卡尔曼滤波估计三轴角速度

和卡尔曼滤波估计三轴加速度的计算过程为：

4.1根据惯性传感器的性能与运动参数建立惯性传感器的误差模型和状态转移关系式，结合时间序列分析法求解出惯性传感器的刻度因子误差M、零点漂移B_f、高斯白噪声U、运动参数间的状态转移矩阵Ψ_k,k-1以及高斯白噪声V；

误差模型如公式（4）所示，状态转移关系式公式（5）所示，

z'_k＝z_k(M+1)+B_f+U       (4)

z_k＝Ψ_k,k-1z_k-1+V              (5)

式（4）和式（5）中，z代表三轴角速度[w_x，w_y，w_z]和三轴加速度[a_x，a_y，a_z]中的某一运动参数，z′_k为k时刻惯性传感器测量获得的某一运动参数的测量值，z_k为k时刻某一运动参数的真实值；

4.2将步骤4.1所求取的结果带入卡尔曼滤波公式对k时刻的运动参数z_k进行滤波，求出该运动参数的估计值

卡尔曼滤波公式如公式(6)所示，

$\left\{\begin{array}{c}{z}_{k,k-1}={\mathrm{\Ψ}}_{k,k-1}{z}_{k-1}\\ \widehat{z_k}=z_{k,k-1}+K_k\left[{z^{\′}}_k(M+1)z_{k,k-1}\right]\\ P_{k,k-1}={\mathrm{\Ψ}}_{k,k-1}{P}_{k-1}{{\mathrm{\Ψ}}_{k,k-1}}^T+Q_{k-1}\\ K_k=P_{k,k-1}{(M+1)}^T{P}_k^{-1}\\ P_k=[I-K_k(M+1)]P_{k,k-1}\end{array}\right.---\left(6\right)$

式（6）中，Q为***的噪声方差，K_k为卡尔曼增益，P_k为状态协方差，z_k,k-1为状态预估值。

5.如权利要求1所述的基于惯性传感器和摄像机的图像配准方法，其特征在于，步骤三中，所述获得摄像机旋转矩阵R和平移矩阵T的计算过程为：

5.1使用欧拉动力学方程将k时刻滤波后的三轴运动角速度w_xk，w_yk，w_zk转换成欧式坐标系中的欧拉角α_k，β_k，γ_k，将所得的α_k，β_k，γ_k带入旋转矩阵R中，完成彩色CCD旋转矩阵R的计算；

所述欧拉动力学方程如公式(7)所示，

$\left[\begin{array}{c}{\stackrel{\·}{\mathrm{\α}}}_k\\ {\stackrel{\·}{\mathrm{\β}}}_k\\ {\stackrel{\·}{\mathrm{\γ}}}_k\end{array}\right]=\left[\begin{array}{ccc}-\sin {\mathrm{\α}}_k/\tan {\mathrm{\β}}_k& {\cos \mathrm{\α}}_k/\tan {\mathrm{\β}}_k& 1\\ {\cos \mathrm{\α}}_k& \sin {\mathrm{\α}}_k& 0\\ \sin {\mathrm{\α}}_k/\sin {\mathrm{\β}}_k& -\cos {\mathrm{\α}}_k/\sin {\mathrm{\β}}_k& 0\end{array}\right]*\left[\begin{array}{c}{w}_{\mathrm{xk}}\\ w_{\mathrm{yk}}\\ w_{\mathrm{zk}}\end{array}\right]---\left(7\right)$

式（7）中，

为欧拉角的一阶导数，将k-1时刻求得的参量α_k-1，β_k-1，γ_k-1作为公式(7)的初值；

所述旋转矩阵R如公式(8)所示：

$R=\left[\begin{array}{ccc}\cos {\mathrm{\γ}}_k\cos {\mathrm{\β}}_k& \cos {\mathrm{\γ}}_k\sin {\mathrm{\β}}_k\sin {\mathrm{\α}}_k-\sin {\mathrm{\γ}}_k& \cos {\mathrm{\γ}}_k\sin {\mathrm{\β}}_k\cos {\mathrm{\α}}_k-\sin {\mathrm{\α}}_k\sin {\mathrm{\γ}}_k\\ \sin {\mathrm{\γ}}_k\cos {\mathrm{\β}}_k& \sin {\mathrm{\α}}_k\sin {\mathrm{\β}}_k\sin {\mathrm{\γ}}_k+\cos {\mathrm{\γ}}_k\cos {\mathrm{\α}}_k& \sin {\mathrm{\γ}}_k\sin {\mathrm{\β}}_k\cos {\mathrm{\α}}_k+\cos {\mathrm{\γ}}_k\sin {\mathrm{\α}}_k\\ -\sin {\mathrm{\β}}_k& \cos {\mathrm{\β}}_k\sin {\mathrm{\α}}_k& \cos {\mathrm{\β}}_k\cos {\mathrm{\α}}_k\end{array}\right]---\left(8\right)$

5.2使用如公式(9)所示的速度位移公式将三轴加速度a_xk，a_yk，a_zk转换成三轴位移t_xk,t_yk,t_zk，

$\left\{\begin{array}{c}{t}_k=v_{k}T+\frac{1}{2}{a}_k{T}^2\\ v_k=v_{k-1}+a_{k-1}\×T\end{array}\right.----\left(9\right)$

式(9)中，t_k代表三轴位移t_xk,t_yk,t_zk中任意一坐标轴方向的位移，T为摄像机的拍摄周期，a_k-1为k-1时刻对应坐标轴方向的加速度，v_k-1和v_k分别为k-1时刻和k时刻的初始速度，给定***初始时刻的初始速度为0。

6.如权利要求1所述的基于惯性传感器和摄像机的图像配准方法，其特征在于，步骤四中根据直接线性变换算法计算出彩色图像和深度图像的仿射变换关系，仿射变换关系如公式(10)所示：

m_ir＝Hm         (10)

式(10)中，m_ir为深度图像中的像素坐标，m为彩色图像中的像素坐标。

7.如权利要求1所述的基于惯性传感器和摄像机的图像配准方法，其特征在于，步骤四中通过反向工程提出的经验关系式获得灰度值与像素深度信息的比例关系，灰度值与像素深度信息的比例关系如公式(11)所示：

u＝l×m         (11)

式(11)中，m(x,y)为像素坐标(x,y)处的灰度值，u为像素坐标(x,y)处对应的深度值，l为比例系数。

8.如权利要求1所述的基于惯性传感器和摄像机的图像配准方法，其特征在于，步骤六中，

所述使用RANSAC算法对彩色图像I的背景光流m₂进行鲁棒性估计的计算过程为：

8.1初始化当前所有的角点m₁(x₁,y₁)，m₂(x₂,y₂)为局内点；

8.2利用DLT算法计算出仿射变换矩阵H₁，根据公式（12）结合角点m₁(x₁,y₁)估计出角点估计值m₂′(x₂′,y₂′)：

m₂′=H₁m₁         (12)

8.3根据公式(13)计算角点m₂(x₂,y₂)和角点估计值m₂′(x₂′,y₂′)的欧氏距离d，若d小于阈值d′，则将角点m₁(x₁,y₁)，m₂(x₂,y₂)置为局内点(即背景光流)，若d大于等于阈值d′，则将角点m₁(x₁,y₁)，m₂(x₂,y₂)置为局外点；当完成一次RANSAC计算后，初始化当前局内点为局外点，重复步骤8.2和8.3，当不再产生新的局内点或达到规定重复次数时，结束当前RANSAC算法，

d＝sqrt((x₂-x₂')²+(y₂-y₂')²)         (13)

所述对彩色图像J进行初始化的过程为：首先遍历该帧所有的角点，若该角点与彩色图像I保留的正确背景光流重合，则将该角点作为RANSAC计算的初始局内点，重复步骤8.1、8.2、8.3所述方法，逐帧完成RANSAC的估计。

Images