CN104144283A - 一种基于改进的卡尔曼滤波的实时数字视频稳像方法 - Google Patents

Abstract

本发明涉及一种数字视频稳像方法，一种基于改进的卡尔曼滤波的实时数字视频稳像方法，包括以下步骤：步骤1、建立帧间运动模型，步骤2、对图像进行分块，步骤3、改进的菱形快速块匹配，步骤4、基于最小二乘的全局运动估计，步骤5、利用改进的卡尔曼滤波器对图像序列进行补偿，步骤6、对含抖动的图像序列进行补偿；本发明对传统的菱形搜索进行了优化，在达到同样或是较好的估计效果的条件下，能减少块匹配算法的运行时间；在运动补偿方面，也能很好的去除抖动分量、很好地保留追拍分量。

Description

一种基于改进的卡尔曼滤波的实时数字视频稳像方法

技术领域

本发明涉及一种数字视频稳像方法，更具体地说，涉及一种基于改进的卡尔曼滤波的实时数字视频稳像方法。

背景技术

数字稳像技术(Digital Image Stabilization,DIS)是图像处理领域中一项新技术，它能将抖动的视频序列中的抖动分量很好地滤除，同时能将有意的追拍分量进行保留。它的这些特点使其在视频监控、手持和机载摄像、遥感领域有着重要的应用价值。在当前的研究背景下，提出一种有效、鲁棒性强的算法，也是本专利的主旨所在。

典型的数字稳像系统包括运动估计单元和运动补偿单元。运动估计单元中，传统的运动估计算法包括块匹配法、位平面匹配法、灰度投影法、特征点匹配等算法。其中，菱形快速搜索算法因具有快速性和较优异的性能而得到了广泛的应用。但传统意义上的菱形搜索算法的性能还有待进一步优化；运动补偿单元中，常用的补偿算法包括卡尔曼滤波器、比例积分差分控制器(Proportional IntegralDifferential,PID)、运动矢量累加(Motion Vector Integration,MVI)、低通滤波器和帧间位置平滑(Frame Position Smooth,FPS)等方法。然而，FPS方法需要消耗大量的时间来重构相机长时间的运动情况，因此不满足实时稳像的要求；MVI等方法可以很好地对帧间抖动进行平滑，但会对追拍分量造成损失，同时有较明显的时滞现象；传统的卡尔曼滤波器方法和PID方法会对追拍分量有一定的保留，但仍然有明显的丢边现象存在，同时对抖动的平滑效果不如MVI性能好。因此，传统卡尔曼滤波器很难做到在平滑抖动分量的同时，对主观的追拍分量进行很好地保留。

发明内容

为了克服现有技术中存在的不足，本发明目的是基于Similarity运动模型，在传统的菱形搜索块匹配算法的基础上，提供一种基于改进的卡尔曼滤波的实时数字视频稳像方法。该方法对传统的菱形搜索算法和卡尔曼滤波器进行了改进，即便是有高速的主观追拍运动存在时，该方法依然能很好地保留追拍分量，同时取得较好的平滑效果，对丢边现象的抑制明显优于其他补偿算法。

为了实现上述发明目的，解决现有技术中所存在的问题，本发明采取的技术方案是：一种基于改进的卡尔曼滤波的实时数字视频稳像方法，包括以下步骤：

步骤1：建立帧间运动模型：选用Similarity运动模型对图像序列的运动进行描述，以反映图像序列帧间的平移、绕光轴旋转和变焦运动，其运动模型为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=s\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin & \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d}_x\\ d_y\end{array}\right]---\left(1\right)$

其中，(x,y)和(x',y')分别表示每个像素点在参考帧和当前帧的坐标，s代表摄像机的变焦系数，θ代表帧间图像的旋转角度，d_x、d_y分别代表水平和垂直方向的偏移量；对于安防监控或手持成像小幅抖动的情况下，图像序列帧间的偏转角θ很小，认为cosθ≈1，sinθ≈θ；同时，变焦系数s保持不变，设为1，这样，(1)式的简化模型为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=s\left[\begin{array}{cc}1& \mathrm{\θ}\\ -\mathrm{\θ}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d}_x\\ d_y\end{array}\right]---\left(2\right)$

步骤2：对图像进行分块：将图像分成一系列8×8的图像块，图像块应选在具有明显图像特征且不为运动前景目标的区域；

步骤3：改进的菱形快速块匹配：按照绝对平均误差函数(Mean AbsoluteDifference,MAD)准则进行匹配；选取以当前帧选定子块为中心的8×8个像素点，进行对应上一帧的相同位置和周围8个方向的搜索，使用菱形搜索算法大模板，搜索的单位步长为2，具体包括以下子步骤(a)、(b)及(c)；

步骤(a)：若匹配后，最小位置为9个位置中的中心位置，则执行步骤(c)，否则，最小位置为外围的8个位置，则执行步骤(b)；

步骤(b)：记录当前最小位置的MAD值，记为MAD₀，并沿着该最小位置的方向再进行一次匹配，得到新的MAD值，记为MAD₁，若MAD₁小于MAD₀，则循环执行步骤(b)，若MAD₁大于或等于MAD₀，则跳转到步骤3；

步骤(c)：匹配后，最小位置为9个位置中的中心位置，则以步骤(a)得到的最小MAD值的位置为初始位置，进行对应上一帧的相同位置和周围4个方向的搜索，使用菱形搜索算法小模板，得到当前的最小MAD值即为最佳匹配位置：

$\mathrm{MAD}(i,j)=\frac{1}{64}\underset{m=1}{\overset{8}{\mathrm{\Σ}}}\underset{n=1}{\overset{8}{\mathrm{\Σ}}}|f_k(m,n)-f_{k-1}(m+i,n+j)|---\left(3\right)$

步骤4：基于最小二乘的全局运动估计：对式(1)进行整理，得到的矩阵方程：

$\left[\begin{array}{cccc}x& -y& 1& 0\\ y& x& 0& 1\end{array}\right]\left[\begin{array}{c}s\cos \mathrm{\θ}\\ s\sin \mathrm{\θ}\\ d_x\\ d_y\end{array}\right]=\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]---\left(4\right)$

根据每个8×8的图像块得到的运动矢量及匹配点对，利用随机采样一致性(Random Sample Consensus,RANSAC)剔除错误的匹配点对，并利用保留下来的正确值以及式(4)对参数s,θ,d_x,d_y进行最小二乘估计；利用最小二乘估计思想估计出的全局运动参数，即得到了当前帧相对于上一帧的全局运动向量，即GMV＝[GMV_s,GMV_θ,GMV_x,GMV_y]＝[s,θ,d_x,d_y]；

步骤5：利用改进的卡尔曼滤波器对图像序列进行补偿，具体包括子步骤(a)、(b)、(c)、(d)、(e)及(f)；

步骤(a)：根据状态空间理论对稳像系统进行建模，系统的状态方程和观测方程分别为：

$\left\{\begin{array}{c}x(k+1)=\mathrm{\Φx}\left(k\right)+\mathrm{\Γu}\left(k\right)+v\left(k\right)\\ y(k+1)=\mathrm{Hx}\left(k\right)+w\left(k\right)\end{array}\right.---\left(5\right)$

其中，Φ、Γ和H分别为传递矩阵、输入矩阵和观测矩阵；x(k)、y(k)和u(k)分别为状态向量、观测向量和控制向量，而v(k)和w(k)代表噪声向量；为了简化起见，采用具有两个状态和单观测的线性系统进行建模，其中Φ，Γ和H的取值如下：

$\mathrm{\Φ}=\left[\begin{array}{cc}1& 1\\ 0& 1\end{array}\right],\mathrm{\Γ}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],H=\left[\begin{array}{cc}1& 0\end{array}\right]---\left(6\right)$

步骤(b)：用PGMV表示全局运动矢量GMV的累加量，PGMV(k-1)表示GMV从第1帧到第k帧的累加值，将PGMV的各个分量分别送入一维的卡尔曼滤波器中进行滤波，则得到：

$\left\{\begin{array}{c}{\mathrm{JMV}}_s(k-1)={\mathrm{CMV}}_s(k-1)={\mathrm{PGMV}}_s(k-1)-{\mathrm{Ko}}_s(k-1)\\ {\mathrm{JMV}}_{\mathrm{\θ}}(k-1)={\mathrm{CMV}}_{\mathrm{\θ}}(k-1)={\mathrm{PGMV}}_{\mathrm{\θ}}(k-1)-{\mathrm{Ko}}_{\mathrm{\θ}}(k-1)\\ {\mathrm{JMV}}_x(k-1)={\mathrm{CMV}}_x(k-1)={\mathrm{PGMV}}_x(k-1)-{\mathrm{Ko}}_x(k-1)\\ {\mathrm{JMV}}_y(k-1)={\mathrm{CMV}}_y(k-1)={\mathrm{PGMV}}_y(k-1)-{\mathrm{Ko}}_y(k-1)\end{array}\right.---\left(7\right)$

其中，Ko(k-1)＝[Ko_s(k-1),Ko_θ(k-1),Ko_x(k-1),Ko_y(k-1)]表示四个卡尔曼滤波器在第k-1时刻的输出，PGMV(k-1)和Ko(k-1)两个矢量之间的差表示第k-1帧的运动补偿矢量，CMV(k-1)＝[CMV_s(k-1),CMV_θ(k-1),CMV_x(k-1),CMV_y(k-1)]第k-1帧的抖动运动矢量用JMV(k-1)＝[JMV_s(k-1),JMV_θ(k-1),JMV_x(k-1),JMV_y(k-1)]表示，其在数值上与第k-1帧的运动补偿矢量CMV(k-1)相等；

步骤(c)：根据卡尔曼滤波器理论对状态进行预测；

$\left\{\begin{array}{c}\hat{x}(k+1|k)=\mathrm{\Φ}\hat{x}\left(k\right|k)+\mathrm{\Γu}\left(k\right)\\ P(k+1|k)=\mathrm{\ΦP}\left(k\right){\mathrm{\Φ}}^T+R_1\end{array}\right.---\left(8\right)$

其中，和是根据k时刻对k+1时刻的估值，而P(k)为协方差矩阵，P(k+1|k)是其根据当前时刻对下一时刻的估值，另外，噪声信号v(k)和w(k)认为是不相关的，R₁为噪声信号v(k)的方差；

步骤(d)：引入场景识别参数，λ＝[λ_s,λ_θ,λ_x,λ_y]，来识别当前场景为追拍运动或是抖动运动，λ值由下式确定：

$\left\{\begin{array}{c}{\mathrm{\λ}}_x=\frac{|\max \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_x\left(i\right)\right)+\min \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_x\left(i\right)\right)|}{2}\\ {\mathrm{\λ}}_y=\frac{|\max \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_y\left(i\right)\right)+\min \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_y\left(i\right)\right)|}{2}\\ {\mathrm{\λ}}_{\mathrm{\θ}}=\frac{|\max \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_{\mathrm{\θ}}\left(i\right)\right)+\min \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_{\mathrm{\θ}}\left(i\right)\right)|}{2}\end{array}\right.---\left(9\right)$

其中，L表示选取一定时间段内JMV的窗长，当某一个参量方向上存在追拍运动时，参数λ在这个方向上的数值大，反之，λ将是一个接近于0的很小值；

步骤(e)：根据λ的值确定第k帧观测噪声方差R(k)的值：

$\left\{\begin{array}{c}{R}_s\left(k\right)=f\left({\mathrm{\λ}}_s\right)\\ R_{\mathrm{\θ}}\left(k\right)=f\left({\mathrm{\λ}}_{\mathrm{\θ}}\right)\\ R_x\left(k\right)=f\left({\mathrm{\λ}}_x\right)\\ R_y\left(k\right)=f\left({\mathrm{\λ}}_y\right)\end{array}\right.---\left(10\right)$

其中f(·)为非线性函数，而且随着λ的增大而递减，这里给出一个参考的非线性函数，参考函数中，参量k_p的经验值取100；

$\left\{\begin{array}{c}{R}_s\left(k\right)=f\left({\mathrm{\λ}}_s\right)=\frac{k_p}{{\mathrm{\λ}}_s^2}\\ R_{\mathrm{\θ}}\left(k\right)=f\left({\mathrm{\λ}}_{\mathrm{\θ}}\right)=\frac{k_p}{{\mathrm{\λ}}_{\mathrm{\θ}}^2}\\ R_x\left(k\right)=f\left({\mathrm{\λ}}_x\right)=\frac{k_p}{{\mathrm{\λ}}_x^2}\\ R_y\left(k\right)=f\left({\mathrm{\λ}}_y\right)=\frac{k_p}{{\mathrm{\λ}}_y^2}\end{array}\right.---\left(11\right)$

当某一个参数方向存在追拍时，由参量确定的R(k)会很小，因而能很好的发挥卡尔曼滤波器的跟踪特性；反之，当前只存在随机抖动，R(k)的值很大，能很好的发挥卡尔曼滤波器的平滑作用；

步骤(f)：利用改进的卡尔曼滤波器和步骤(e)确定的观测噪声方差R(k)，对状态进行更新：

$\left\{\begin{array}{c}K\left(k\right)=P\left(k\right|k-1)C^T{(H\left(k\right)P\left(k\right|k-1)H{\left(k\right)}^T+R\left(k\right))}^{-1}\\ \hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+K\left(k\right)(y\left(k\right)-C\hat{x}\left(k\right|k-1))\\ P(k+1)=(1-K\left(k\right)C\left(k\right))P(k+1|k)\end{array}\right.---\left(12\right)$

这里，R(k)为标量而非矢量，这是由于采用系统模型是单观测的，若采用多观测的系统模型，则应用R(k)来表示观测噪声的协方差矩阵；

步骤6：对含抖动的图像序列进行补偿：将运动补偿得到的运动补偿向量的各个参量，CMV＝[s',θ',d'_x,d'_y]以及原图像的像素位置代入下式，进行补偿：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=s^{\′}\left[\begin{array}{cc}\cos {\mathrm{\θ}}^{\′}& -\sin {\mathrm{\θ}}^{\′}\\ \sin {\mathrm{\θ}}^{\′}& \cos {\mathrm{\θ}}^{\′}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d_x}^{\′}\\ {d_y}^{\′}\end{array}\right]---\left(13\right)$

并将像素赋给新的坐标位置(x',y')，即完成了稳像；如果采取简化后的模型，应用下式进行补偿：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}1& -{\mathrm{\θ}}^{\′}\\ {\mathrm{\θ}}^{\′}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d_x}^{\′}\\ {d_y}^{\′}\end{array}\right]---\left(14\right)$

本发明有益效果是：一种基于改进的卡尔曼滤波的实时数字视频稳像方法，包括以下步骤：步骤1、建立帧间运动模型，步骤2、对图像进行分块，步骤3、改进的菱形快速块匹配，步骤4、基于最小二乘的全局运动估计，步骤5、利用改进的卡尔曼滤波器对图像序列进行补偿，步骤6、对含抖动的图像序列进行补偿；与已有技术相比，本发明对传统的菱形搜索进行了优化，在达到同样或是较好的估计效果的条件下，能减少块匹配算法的运行时间；在运动补偿方面，也能很好的去除抖动分量、很好地保留追拍分量。

附图说明

图1是改进的菱形搜索(大模板)的示意图。

图2是菱形搜索(小模板)的示意图。

图3是基于优化的菱形搜索和改进的卡尔曼滤波的实时数字视频稳像方法的系统框图。

图4是原始视频(A)参考帧图像图。

图5是原始视频(A)当前帧图像图。

图6是用传统菱形搜索算法处理后图像图。

图7是用本发明方法改进的菱形搜索算法处理后图像图。

图8是原始视频(B)参考帧图像图。

图9是原始视频(B)当前帧图像。

图10是用本发明方法处理后图像图。

图11是原始视频(C)的连续图像序列图。

图12是用本发明方法改进的卡尔曼滤波器处理后连续图像序列图。

图13是用传统的卡尔曼滤波器处理后的连续图像序列图。

图14是利用PID、MVI、传统的卡尔曼滤波器和本方法存在追拍时的补偿效果对比图。

具体实施方式

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

一种基于改进的卡尔曼滤波的实时数字视频稳像方法，包括以下步骤：

步骤1：建立帧间运动模型：选用Similarity运动模型对图像序列的运动进行描述，以反映图像序列帧间的平移、绕光轴旋转和变焦运动，其运动模型为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=s\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d}_x\\ d_y\end{array}\right]---\left(15\right)$

其中，(x,y)和(x',y')分别表示每个像素点在参考帧和当前帧的坐标，s代表摄像机的变焦系数，θ代表帧间图像的旋转角度，d_x、d_y分别代表水平和垂直方向的偏移量；对于安防监控或手持成像小幅抖动的情况下，图像序列帧间的偏转角θ很小，认为cosθ≈1，sinθ≈θ；同时，变焦系数s保持不变，设为1，这样，(1)式的简化模型为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=s\left[\begin{array}{cc}1& \mathrm{\θ}\\ -\mathrm{\θ}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d}_x\\ d_y\end{array}\right]---\left(16\right)$

步骤2：对图像进行分块：将图像分成一系列8×8的图像块，图像块应选在具有明显图像特征且不为运动前景目标的区域；

步骤3：改进的菱形快速块匹配：按照绝对平均误差函数(Mean AbsoluteDifference,MAD)准则进行匹配；选取以当前帧选定子块为中心的8×8个像素点，进行对应上一帧的相同位置和周围8个方向的搜索，使用菱形搜索算法大模板，搜索的单位步长为2，具体包括以下子步骤(a)、(b)及(c)；

步骤(a)：若匹配后，最小位置为9个位置中的中心位置，则执行步骤(c)，否则，最小位置为外围的8个位置，则执行步骤(b)；

步骤(b)：记录当前最小位置的MAD值，记为MAD₀，并沿着该最小位置的方向再进行一次匹配，得到新的MAD值，记为MAD₁，若MAD₁小于MAD₀，则循环执行步骤(b)，若MAD₁大于或等于MAD₀，则跳转到步骤3；

步骤(c)：匹配后，最小位置为9个位置中的中心位置，则以步骤(a)得到的最小MAD值的位置为初始位置，进行对应上一帧的相同位置和周围4个方向的搜索，使用菱形搜索算法小模板，得到当前的最小MAD值即为最佳匹配位置：

$\mathrm{MAD}(i,j)=\frac{1}{64}\underset{m=1}{\overset{8}{\mathrm{\Σ}}}\underset{n=1}{\overset{8}{\mathrm{\Σ}}}|f_k(m,n)-f_{k-1}(m+i,n+j)|---\left(17\right)$

步骤4：基于最小二乘的全局运动估计：对式(1)进行整理，得到的矩阵方程：

$\left[\begin{array}{cccc}x& -y& 1& 0\\ y& x& 0& 1\end{array}\right]\left[\begin{array}{c}s\cos \mathrm{\θ}\\ s\sin \mathrm{\θ}\\ d_x\\ d_y\end{array}\right]=\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]---\left(18\right)$

根据每个8×8的图像块得到的运动矢量及匹配点对，利用随机采样一致性(Random Sample Consensus,RANSAC)剔除错误的匹配点对，并利用保留下来的正确值以及式(4)对参数s,θ,d_x,d_y进行最小二乘估计；利用最小二乘估计思想估计出的全局运动参数，即得到了当前帧相对于上一帧的全局运动向量，即GMV＝[GMV_s,GMV_θ,GMV_x,GMV_y]＝[s,θ,d_x,d_y]；

步骤5：利用改进的卡尔曼滤波器对图像序列进行补偿，具体包括子步骤(a)、(b)、(c)、(d)、(e)及(f)；

步骤(a)：根据状态空间理论对稳像系统进行建模，系统的状态方程和观测方程分别为：

$\left\{\begin{array}{c}x(k+1)=\mathrm{\Φx}\left(k\right)+\mathrm{\Γu}\left(k\right)+v\left(k\right)\\ y(k+1)=\mathrm{Hx}\left(k\right)+w\left(k\right)\end{array}\right.---\left(19\right)$

其中，Φ、Γ和H分别为传递矩阵、输入矩阵和观测矩阵；x(k)、y(k)和u(k)分别为状态向量、观测向量和控制向量，而v(k)和w(k)代表噪声向量；为了简化起见，采用具有两个状态和单观测的线性系统进行建模，其中Φ，Γ和H的取值如下：

$\mathrm{\Φ}=\left[\begin{array}{cc}1& 1\\ 0& 1\end{array}\right],\mathrm{\Γ}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],H=\left[\begin{array}{cc}1& 0\end{array}\right]---\left(20\right)$

步骤(b)：用PGMV表示全局运动矢量GMV的累加量，PGMV(k-1)表示GMV从第1帧到第k帧的累加值，将PGMV的各个分量分别送入一维的卡尔曼滤波器中进行滤波，则得到：

$\left\{\begin{array}{c}{\mathrm{JMV}}_s(k-1)={\mathrm{CMV}}_s(k-1)={\mathrm{PGMV}}_s(k-1)-{\mathrm{Ko}}_s(k-1)\\ {\mathrm{JMV}}_{\mathrm{\θ}}(k-1)={\mathrm{CMV}}_{\mathrm{\θ}}(k-1)={\mathrm{PGMV}}_{\mathrm{\θ}}(k-1)-{\mathrm{Ko}}_{\mathrm{\θ}}(k-1)\\ {\mathrm{JMV}}_x(k-1)={\mathrm{CMV}}_x(k-1)={\mathrm{PGMV}}_x(k-1)-{\mathrm{Ko}}_x(k-1)\\ {\mathrm{JMV}}_y(k-1)={\mathrm{CMV}}_y(k-1)={\mathrm{PGMV}}_y(k-1)-{\mathrm{Ko}}_y(k-1)\end{array}\right.---\left(21\right)$

其中，Ko(k-1)＝[Ko_s(k-1),Ko_θ(k-1),Ko_x(k-1),Ko_y(k-1)]表示四个卡尔曼滤波器在第k-1时刻的输出，PGMV(k-1)和Ko(k-1)两个矢量之间的差表示第k-1帧的运动补偿矢量，CMV(k-1)＝[CMV_s(k-1),CMV_θ(k-1),CMV_x(k-1),CMV_y(k-1)]第k-1帧的抖动运动矢量用JMV(k-1)＝[JMV_s(k-1),JMV_θ(k-1),JMV_x(k-1),JMV_y(k-1)]表示，其在数值上与第k-1帧的运动补偿矢量CMV(k-1)相等；

步骤(c)：根据卡尔曼滤波器理论对状态进行预测；

$\left\{\begin{array}{c}\hat{x}(k+1|k)=\mathrm{\Φ}\hat{x}\left(k\right|k)+\mathrm{\Γu}\left(k\right)\\ P(k+1|k)=\mathrm{\ΦP}\left(k\right){\mathrm{\Φ}}^T+R_1\end{array}\right.---\left(22\right)$

其中，和是根据k时刻对k+1时刻的估值，而P(k)为协方差矩阵，P(k+1|k)是其根据当前时刻对下一时刻的估值，另外，噪声信号v(k)和w(k)认为是不相关的，R₁为噪声信号v(k)的方差；

步骤(d)：引入场景识别参数，λ＝[λ_s,λ_θ,λ_x,λ_y]，来识别当前场景为追拍运动或是抖动运动，λ值由下式确定：

$\left\{\begin{array}{c}{\mathrm{\λ}}_x=\frac{|\max \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_x\left(i\right)\right)+\min \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_x\left(i\right)\right)|}{2}\\ {\mathrm{\λ}}_y=\frac{|\max \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_y\left(i\right)\right)+\min \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_y\left(i\right)\right)|}{2}\\ {\mathrm{\λ}}_{\mathrm{\θ}}=\frac{|\max \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_{\mathrm{\θ}}\left(i\right)\right)+\min \left(\underset{i=k-L}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{JMV}}_{\mathrm{\θ}}\left(i\right)\right)|}{2}\end{array}\right.---\left(23\right)$

其中，L表示选取一定时间段内JMV的窗长，当某一个参量方向上存在追拍运动时，参数λ在这个方向上的数值大，反之，λ将是一个接近于0的很小值；

步骤(e)：根据λ的值确定第k帧观测噪声方差R(k)的值：

$\left\{\begin{array}{c}{R}_s\left(k\right)=f\left({\mathrm{\λ}}_s\right)\\ R_{\mathrm{\θ}}\left(k\right)=f\left({\mathrm{\λ}}_{\mathrm{\θ}}\right)\\ R_x\left(k\right)=f\left({\mathrm{\λ}}_x\right)\\ R_y\left(k\right)=f\left({\mathrm{\λ}}_y\right)\end{array}\right.---\left(24\right)$

其中f(·)为非线性函数，而且随着λ的增大而递减，这里给出一个参考的非线性函数，参考函数中，参量k_p的经验值取100；

$\left\{\begin{array}{c}{R}_s\left(k\right)=f\left({\mathrm{\λ}}_s\right)=\frac{k_p}{{\mathrm{\λ}}_s^2}\\ R_{\mathrm{\θ}}\left(k\right)=f\left({\mathrm{\λ}}_{\mathrm{\θ}}\right)=\frac{k_p}{{\mathrm{\λ}}_{\mathrm{\θ}}^2}\\ R_x\left(k\right)=f\left({\mathrm{\λ}}_x\right)=\frac{k_p}{{\mathrm{\λ}}_x^2}\\ R_y\left(k\right)=f\left({\mathrm{\λ}}_y\right)=\frac{k_p}{{\mathrm{\λ}}_y^2}\end{array}\right.---\left(25\right)$

当某一个参数方向存在追拍时，由参量确定的R(k)会很小，因而能很好的发挥卡尔曼滤波器的跟踪特性；反之，当前只存在随机抖动，R(k)的值很大，能很好的发挥卡尔曼滤波器的平滑作用；

步骤(f)：利用改进的卡尔曼滤波器和步骤(e)确定的观测噪声方差R(k)，对状态进行更新：

$\left\{\begin{array}{c}K\left(k\right)=P\left(k\right|k-1)C^T{(H\left(k\right)P\left(k\right|k-1)H{\left(k\right)}^T+R\left(k\right))}^{-1}\\ \hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+K\left(k\right)(y\left(k\right)-C\hat{x}\left(k\right|k-1))\\ P(k+1)=(1-K\left(k\right)C\left(k\right))P(k+1|k)\end{array}\right.---\left(26\right)$

这里，R(k)为标量而非矢量，这是由于采用系统模型是单观测的，若采用多观测的系统模型，则应用R(k)来表示观测噪声的协方差矩阵；

步骤6：对含抖动的图像序列进行补偿：将运动补偿得到的运动补偿向量的各个参量，CMV＝[s',θ',d'_x,d'_y]以及原图像的像素位置代入下式，进行补偿：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=s^{\′}\left[\begin{array}{cc}\cos {\mathrm{\θ}}^{\′}& -\sin {\mathrm{\θ}}^{\′}\\ \sin {\mathrm{\θ}}^{\′}& \cos {\mathrm{\θ}}^{\′}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d_x}^{\′}\\ {d_y}^{\′}\end{array}\right]---\left(27\right)$

并将像素赋给新的坐标位置(x',y')，即完成了稳像；如果采取简化后的模型，应用下式进行补偿：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}1& -{\mathrm{\θ}}^{\′}\\ {\mathrm{\θ}}^{\′}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{d_x}^{\′}\\ {d_y}^{\′}\end{array}\right]---\left(28\right)$

表1是本发明提出的实时数字视频稳像方法处理不同分辨率视频的耗时(每帧取36个图像块)。

本发明优点在于：即便是有高速的主观追拍运动存在时，该方法依然能很好地保留追拍分量，同时取得较好的平滑效果，对丢边现象的抑制明显优于其他补偿算法。

Claims (1)

1.一种基于改进的卡尔曼滤波的实时数字视频稳像方法，其特征在于包括以下步骤： 

步骤1：建立帧间运动模型：选用Similarity运动模型对图像序列的运动进行描述，以反映图像序列帧间的平移、绕光轴旋转和变焦运动，其运动模型为 

其中，(x,y)和(x',y')分别表示每个像素点在参考帧和当前帧的坐标，s代表摄像机的变焦系数，θ代表帧间图像的旋转角度，d_x、d_y分别代表水平和垂直方向的偏移量；对于安防监控或手持成像小幅抖动的情况下，图像序列帧间的偏转角θ很小，认为cosθ≈1，sinθ≈θ；同时，变焦系数s保持不变，设为1，这样，(1)式的简化模型为： 

步骤2：对图像进行分块：将图像分成一系列8×8的图像块，图像块应选在具有明显图像特征且不为运动前景目标的区域； 

步骤3：改进的菱形快速块匹配：按照绝对平均误差函数(Mean Absolute Difference,MAD)准则进行匹配；选取以当前帧选定子块为中心的8×8个像素点，进行对应上一帧的相同位置和周围8个方向的搜索，使用菱形搜索算法大模板，搜索的单位步长为2，具体包括以下子步骤(a)、(b)及(c)； 

步骤(a)：若匹配后，最小位置为9个位置中的中心位置，则执行步骤(c)，否则，最小位置为外围的8个位置，则执行步骤(b)； 

步骤(b)：记录当前最小位置的MAD值，记为MAD₀，并沿着该最小位置的方向再进行一次匹配，得到新的MAD值，记为MAD₁，若MAD₁小于MAD₀，则循环执行步骤(b)，若MAD₁大于或等于MAD₀，则跳转到步骤3； 

步骤(c)：匹配后，最小位置为9个位置中的中心位置，则以步骤(a)得到的最小MAD值的位置为初始位置，进行对应上一帧的相同位置和周围4个方向的搜索，使用菱形搜索算法小模板，得到当前的最小MAD值即为最佳匹配位置： 

步骤4：基于最小二乘的全局运动估计：对式(1)进行整理，得到的矩阵方程： 

根据每个8×8的图像块得到的运动矢量及匹配点对，利用随机采样一致性(Random Sample Consensus,RANSAC)剔除错误的匹配点对，并利用保留下来的正确值以及式(4)对参数s,θ,d_x,d_y进行最小二乘估计；利用最小二乘估计思想估计出的全局运动参数，即得到了当前帧相对于上一帧的全局运动向量，即GMV＝[GMV_s,GMV_θ,GMV_x,GMV_y]＝[s,θ,d_x,d_y]； 

步骤5：利用改进的卡尔曼滤波器对图像序列进行补偿，具体包括子步骤(a)、(b)、(c)、(d)、(e)及(f)； 

步骤(a)：根据状态空间理论对稳像系统进行建模，系统的状态方程和观测方程分别为： 

其中，Φ、Γ和H分别为传递矩阵、输入矩阵和观测矩阵；x(k)、y(k)和u(k)分别为状态向量、观测向量和控制向量，而v(k)和w(k)代表噪声向量；为了简化起见，采用具有两个状态和单观测的线性系统进行建模，其中Φ，Γ和H的取值如下： 

步骤(b)：用PGMV表示全局运动矢量GMV的累加量，PGMV(k-1)表示GMV从第1帧到第k帧的累加值，将PGMV的各个分量分别送入一维的卡尔曼滤波器中进行滤波，则得到： 

其中，Ko(k-1)＝[Ko_s(k-1),Ko_θ(k-1),Ko_x(k-1),Ko_y(k-1)]表示四个卡尔曼滤波器在第k-1时刻的输出，PGMV(k-1)和Ko(k-1)两个矢量之间的差表示第k-1帧的运动补偿矢量，CMV(k-1)＝[CMV_s(k-1),CMV_θ(k-1),CMV_x(k-1),CMV_y(k-1)]第k-1帧的抖动运动矢量用JMV(k-1)＝[JMV_s(k-1),JMV_θ(k-1),JMV_x(k-1),JMV_y(k-1)]表示，其在数值上与第k-1帧的运动补偿矢量CMV(k-1)相等； 

步骤(c)：根据卡尔曼滤波器理论对状态进行预测； 

其中，和是根据k时刻对k+1时刻的估值，而P(k)为协方差矩阵，P(k+1|k)是其根据当前时刻对下一时刻的估值，另外，噪声信号v(k)和w(k)认为是不相关的，R₁为噪声信号v(k)的方差； 

步骤(d)：引入场景识别参数，λ＝[λ_s,λ_θ,λ_x,λ_y]，来识别当前场景为追拍运动或是抖动运动，λ值由下式确定： 

其中，L表示选取一定时间段内JMV的窗长，当某一个参量方向上存在追拍运动时，参数λ在这个方向上的数值大，反之，λ将是一个接近于0的很小值； 

步骤(e)：根据λ的值确定第k帧观测噪声方差R(k)的值： 

其中f(·)为非线性函数，而且随着λ的增大而递减，这里给出一个参考的非线性函数，参考函数中，参量k_p的经验值取100； 

当某一个参数方向存在追拍时，由参量确定的R(k)会很小，因而能很好的发挥卡尔曼滤波器的跟踪特性；反之，当前只存在随机抖动，R(k)的值很大，能很好的发挥卡尔曼滤波器的平滑作用； 

步骤(f)：利用改进的卡尔曼滤波器和步骤(e)确定的观测噪声方差R(k)， 对状态进行更新： 

这里，R(k)为标量而非矢量，这是由于采用系统模型是单观测的，若采用多观测的系统模型，则应用R(k)来表示观测噪声的协方差矩阵； 

步骤6：对含抖动的图像序列进行补偿：将运动补偿得到的运动补偿向量的各个参量，CMV＝[s',θ',d'_x,d'_y]以及原图像的像素位置代入下式，进行补偿： 

并将像素赋给新的坐标位置(x',y')，即完成了稳像；如果采取简化后的模型，则应用下式进行补偿： 

。