基于灰度积分投影互相关函数特征的运动参数估计方法

技术领域

本发明属于图像分析技术领域，特别是一种基于灰度积分投影互相关函数特征的运 动参数估计方法。

背景技术

数字图像稳定系统基于数字图像处理技术，检测图像抖动向量并进行实时补偿，它 不需要依托陀螺或棱镜等专用设备，检测精度和补偿精度只与算法实现有关，因此设置 灵活，且能处理大范围抖动，可直接应用于现已普及的视频监视系统后端，与现有系统 兼容，因此得到大范围的研究和应用。数字稳像系统一般包括两个处理单元：运动估计 单元和运动补偿单元。运动估计单元获取抖动图像序列的全局运动参数，运动补偿单元 根据运动估计单元的结果，生成补偿运动向量，获得连续的、平滑的图像序列。运动估 计单元是数字稳像系统的核心和基础。常用的运动估计算法包括图像块匹配法(Blocks  Matching Algorithm，BMA)、灰度积分投影法(Intensity Integral Projections Algorithm， IIPA)、特征跟踪法(Features Tracking Algorithm，FTA)等。

海域安全监控应用领域中的视频图像多数具有低对比度、高分辨率等特点，且风致 振动导致图像抖动严重。去抖操作的全局运动参数估计采用特征跟踪法(FTA)时，发现 低对比度、高分辨率图像的灰度变化不剧烈，角点等特征很不明显，角点的搜索阈值非 常小，且经常搜索失败。采用对图像小邻域特征不敏感的块匹配法(BMA)时需要事先设 定好大小、位置合适的匹配块，这在海域安全监控应用中很难满足，而传统的灰度积分 投影法(IIPA)虽不需要人为过多干预，但在直接应用时的实时性和准确性无法达到应用 要求。

发明内容

本发明的目的在于提供一种精度高、效果好的基于灰度积分投影互相关函数特征的 运动参数估计方法，有效估计低对比度、高分辨率图像的全局运动参数。

实现本发明目的的技术解决方案为：一种基于灰度积分投影互相关函数特征的运动 参数估计方法，包括以下步骤：

步骤1，将目标图像等分为N×N个矩形分块，并对每个分块进行灰度积分投影， 确定每个分块灰度积分投影的方差，遴选方差较大的R个区域作为备选区域，R为整数 且2≤R≤N；

步骤2，对连续两帧目标图像中备选区域的灰度积分投影进行互相关运算，确定互 相关函数极值点的二阶导数，对灰度积分投影的方差和互相关函数极值点的二阶导数进 行加权求和，获得各个备选区域的优劣因子，将优劣因子降序排列，选择排列靠前的Q 个区域作为优质区域，其中Q≤R；

步骤3，分别对所得的优质区域进行局部运动参数估计，获得每个区域的X轴或Y 轴平移运动参数，再对这些局部运动参数进行加权求和获得全局运动参数。

本发明与现有技术相比，其显著优点为：(1)不需要事先了解图像的特点，通过图 像分块计算各个分块的灰度投影以及互相关函数，利用灰度积分投影的方差和互相关函 数在极值点的二阶导数定量描述各个分块区域的优劣，筛选优质区域进行后续计算，避 免了对匹配块大小和位置、以及对角点检测阈值的依赖；(2)运算速度快，由于只对灰 度积分投影方差较大的若干分块区域进行互相关计算，而不进行全幅面的灰度积分投影 互相关运算，因此计算量大大降低；(3)算法精度高，避免了以整幅面低对比度图像的 灰度积分投影为匹配对象，而是基于区域优劣因子选择若干优质区域进行匹配和参数估 计，因此得到的全局运动参数相对准确，精度比IIPA提高了约30％；(4)鲁棒性强， 由于基于分块区域的灰度积分投影，因此随机噪声的影响很小。

附图说明

图1是本发明基于灰度积分投影互相关函数特征的运动参数估计方法的流程图。

图2是本发明第(p,q)分块子图像的灰度积分投影曲线图，其中(a)为行灰度积分投影 曲线图，(b)为列灰度积分投影曲线图。

图3是本发明连续两帧目标图像中备选区域的灰度积分投影互相关曲线图，其中(a) 为行灰度积分投影的互相关曲线图，(b)为列灰度积分投影的互相关曲线图。

图4是本发明连续两帧目标图像中不同备选区域的行灰度积分投影互相关曲线图。

图5是本发明实施例中几种不同方法进行数字稳像后的图像灰度差平方和棒状图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。

如图1所示，本发明基于灰度积分投影互相关函数特征的运动参数估计方法，包括 以下步骤：

步骤1，将目标图像划分为N×N个同等大小矩形分块，并对每个分块进行灰度积 分投影，确定每个分块灰度积分投影的方差，遴选方差较大的R个区域作为备选区域， R为整数且2≤R≤N，具体如下：

(1.1)将目标图像划分为N×N个大小相等的分块，其中N为整数且N≥2；N值根 据图像的分辨率和抖动最大幅值确定。如果图像抖动的横向最大平移为40像素，纵向 25像素，则分块大小设定为400×250左右是合适的。对于高清网络型摄像机，1920×1080 分辨率的图像可采用4×4分区，每个区域大小为480×270，可满足应用需求。

(1.2)对每个分块进行灰度积分投影，如式(1)所示：

$\left(\begin{array}{c}{\mathrm{Col}}_{k,(p,q)}\left(j\right)=\underset{i}{\Sigma }{G}_{k,(p,q)}(i,j)\\ {\mathrm{Row}}_{k,(p,q)}\left(i\right)=\underset{j}{\Sigma }{G}_{k,(p,q)}(i,j)\end{array}\right)---\left(1\right)$

式中，G_k,(p,q)(i,j)为第k帧目标图像第(p,q)分块子图像的灰度矩阵，(p,q)表示第p 行第q列的子区域，i为矩阵的横坐标、j为矩阵的纵坐标，Col_k,(p,q)(j)为第k帧目标图 像第(p,q)分块子图像第j列的灰度积分投影、Row_k,(p,q)(i)为第k帧目标图像第(p,q)分块 子图像第i行的灰度积分投影；Col_k,(p,q)(j)和Row_k,(p,q)(i)组成了两个一维数组，数组大 小分别等于图像矩阵的列数和行数。分别以行序号和列序号为横坐标，以灰度投影值为 纵坐标，可绘制出该图像的两幅灰度投影曲线，如图2所示，图2(a)为第(p,q)分块子图 像的行灰度积分投影曲线，(b)为第(p,q)分块子图像的列灰度积分投影曲线。

确定每个分块灰度积分投影的方差，如式(2)所示：

$\left(\begin{array}{c}{Z}_M=\frac{1}{M}\underset{t=1}{\overset{M}{\Sigma }}{Z}_{k,(p,q)}\left(t\right)\\ {\sigma }^2=\frac{1}{M-1}\underset{t=1}{\overset{M}{\Sigma }}{(Z_{k,(p,q)}\left(t\right)-Z_M)}^2\end{array}\right),Z=x,y---\left(2\right)$

式中，Z_k,(p,q)(t)为第k帧目标图像第(p,q)分块子图像的灰度积分投影，Z_M为第k 帧目标图像第(p,q)分块子图像的灰度积分投影均值，M为第k帧目标图像第(p,q)分块子 图像的灰度积分投影元素个数，σ²为第k帧目标图像第(p,q)分块子图像的灰度积分投影 的方差，Z＝x表示行灰度积分投影、Z＝y表示列灰度积分投影；

(1.3)以灰度积分投影的方差大小对这些分块排序，遴选方差较大的R个区域作 为备选区域，R为整数且2≤R≤N。R值一般根据经验确定，当图像16等分(即N＝4) 时，R最大取值4可满足应用需求。R值还需要利用方差的梯度值进行适合性验证，即 利用排序靠前的R+1个方差进行梯度计算，如果梯度大于设定阈值D_max，则对应该梯度 的序号减1作为新的R值，其后的方差对应的区域舍弃。D_max可根据经验取值4～6。

步骤2，对连续两帧目标图像中备选区域的灰度积分投影进行互相关运算，确定互 相关函数极值点的二阶导数，对灰度积分投影的方差和互相关函数极值点的二阶导数进 行加权求和，获得各个备选区域的优劣因子，将优劣因子降序排列，选择排列靠前的Q 个区域作为优质区域，其中Q≤R；具体如下：

(2.1)采用最小差平方和(SSD)的方法对连续两帧目标图像中备选区域的灰度积分 投影进行互相关运算，如式(3)所示：

$\left(\begin{array}{c}\left(\begin{array}{c}{C}_{c,(p,q)}\left(u\right)=\underset{j=1}{\overset{H}{\Sigma }}{({\mathrm{Col}}_{k,(p,q)}(u+j-1)-{\mathrm{Col}}_{l,(p,q)}(m+j))}^2\\ 1\le u\le 2m+1\end{array}\right)\\ \left(\begin{array}{c}{C}_{r,(p,q)}\left(v\right)=\underset{i=1}{\overset{W}{\Sigma }}{({\mathrm{Row}}_{k,(p,q)}(v+i-1)-{\mathrm{Row}}_{l,(p,q)}(n+i))}^2\\ 1\le v\le 2n+1\end{array}\right)\end{array}\right)---\left(3\right)$

式中，k、l为连续的两帧目标图像编号，C_c,(p,q)(u)表示列投影的互相关值、C_r,(p,q)(v) 表示行投影的互相关值，m为水平方向的最大抖动位移，n为垂直方向的最大抖动位移， H为第k帧目标图像第(p,q)分块子图像灰度矩阵的总列数，W为第k帧目标图像第(p,q) 分块子图像灰度矩阵的总行数。图3为连续两帧目标图像中备选区域的灰度积分投影互 相关曲线，其中图3(a)为行灰度积分投影的互相关曲线图，图3(b)为列灰度积分投影的 互相关曲线图。

(2.2)确定灰度积分投影互相关函数在极值点的二阶导数，连续两帧目标图像第 (p,q)分块子图像的行或列灰度积分投影互相关函数在极值点的二阶导数Ds_z,(p,q)为：

${\mathrm{Ds}}_{z,(p,q)}=\left|\frac{d^2{C}_{z,(p,q)}\left(u\right)}{d{u}^2}{|}_{u=u_{\min }}\right|,Z=x,y---\left(4\right)$

式中，C_z,(p,q)(u_min)为连续两帧目标图像第(p,q)分块子图像的灰度积分投影互相关曲 线的极值点的灰度投影值，Z＝x表示行灰度积分投影，Z＝y表示列灰度积分投影；

传统灰度积分投影算法只注重互相关函数的极值点位置，以判断运动参数，而实际 上互相关函数的极值并不是越小越好，还要看极值点的变化率特征。极值点附近的梯度 变化如果很小，则极值点不容易判断。从物理意义上解释就是连续两帧抖动图像的灰度 投影曲线的相关性不强。因此可以对低对比度高分辨率图像细分获得的不同区域，进行 灰度积分投影的互相关函数特性分析，重点是极值点位置的梯度变化率(用二阶导数描 述)。图4为连续两帧图像三个不同区域的灰度投影相关性曲线，横坐标为前一帧相对 于后续帧的抖动偏移量，纵坐标为互相关的计算结果。可以看出，同一图像不同区域的 极值点位置并不一致，极值数值、极值点附近的曲线变化剧烈程度，以及对应的平移量 都不同。实线代表第(2,1)分块子图像的极值点横坐标大约为2，点划线代表第(3,2)分块 子图像的极值点横坐标大约为5，而虚线代表第(4,4)分块子图像的极值点横坐标大约为 3。极值点附近的变化幅度越大，极值点的判断越明显，因此应优先选择极值点变化幅 度大的黑色虚线曲线对应的区域作为优选区域，以其偏移量作为整幅图像偏移量的主要 参考。

由于图像为离散数据集，将导数改为差分；为简化计算，略去像素间隔的除数，使 用式(5)中的Ds′_x,(p,q)代替Ds_z,(p,q)，即：

Ds′_z,(p,q)＝C_z,(p,q)(u_min-1)+C_z,(p,q)(u_min+1)-2C_z,(p,q)(u_min)，Z＝x、y    (5)

其中，C_z,(p,q)(u_min-1)为连续两帧目标图像第(p,q)分块子图像的灰度积分投影互相 关曲线的极值点的前一点的灰度投影值、C_z,(p,q)(u_min+1)为连续两帧目标图像第(p,q)分 块子图像的灰度积分投影互相关曲线的极值点的后一点的灰度投影值，Z＝x表示行灰度 积分投影，Z＝y表示列灰度积分投影；

对灰度积分投影的方差和互相关函数极值点的二阶导数进行加权求和，获得各个备 选区域的优劣因子，如下式所示：

s_z,(p,q)＝α·σ²+β·Ds′_z,(p,q)，Z＝x、y         (6)

式(6)中，s_z,(p,q)为第(p,q)分块子图像的行或列区域优劣因子，Z＝x时表示第(p,q)分块 子图像的行区域优劣因子、Z＝y时表示第(p,q)分块子图像的列区域优劣因子，α为灰度 积分投影的方差因子，β为连续两帧目标图像第(p,q)分块子图像的灰度积分投影互相关 函数的极值特征因子，为简化计算，α取值0，β取值1，则式(6)简化为：

s_z,(p,q)＝Ds′_z,(p,q)，Z＝x、y          (7)

即直接采用第(p,q)分块子图像的行灰度积分投影互相关函数在极值点的二阶导数作 为该分块的行区域优劣因子，采用第(p,q)分块子图像的列灰度积分投影互相关函数在极 值点的二阶导数作为该分块的列区域优劣因子；

(2.3)将备选区域按照区域优劣因子的降序排列，选择排序靠前的Q个区域为优 质区域，进行后续的局部运动参数估计。Q取值可直接等于R，但Q值也需要利用优劣 因子的梯度值进行适合性验证：设第Q+1位的区域优劣因子梯度变化为τ_Q+1，第Q位的 区域优劣因子梯度变化为τ_Q，如果则删除排在第Q位后面的区域，D_max可 根据经验取值4～6。

步骤3，分别对所得的优质区域进行局部运动参数估计，获得各个优质区域X轴或 Y轴的平移运动参数，再对这些局部运动参数进行加权求和获得全局运动参数。具体步 骤如下：

(3.1)分别对所得的优质区域进行局部运动参数估计，采用式(8)分别计算行、列 优质区域的平移运动参数：

$\left(\begin{array}{c}{d}_{x,(p,q)}=m+1-u_{\min }\\ d_{y,(p,q)}=n+1-v_{\min }\end{array}\right)---\left(8\right)$

式(8)中，u_min为C_c,(p,q)(u)取最小值时的u值，v_min为C_r,(p,q)(v)取最小值时的v值， d_x,(p,q)和d_y,(p,q)分别表示水平和垂直方向的偏移，d_x,(p,q)为正时表示第l帧目标图像相对 于第k帧目标图像左移了|d_x,(p,q)|个像素，d_x,(p,q)为负时表示第l帧目标图像相对于第k 帧目标图像右移了|d_x,(p,q)|个像素；d_y,(p,q)为正时表示第l帧目标图像相对于第k帧目标 图像上移了|d_y,(p,q)|个像素，d_y,(p,q)为负时表示第l帧目标图像相对于第k帧目标图像下 移了|d_y,(p,q)|像素；在全幅图像运动一致性的假设下，可通过求解互相关曲线的最小值， 判断图像的偏移值。

(3.2)对各个优质区域的局部运动参数进行加权求和获得全局运动参数，如下式：

$\left(\begin{array}{c}{f}_{z,(p,q)}=\frac{s_{z,(p,q)}^{\prime }}{\underset{\lambda =1}{\overset{Q}{\Sigma }}{s}_{z,\lambda }^{\prime }}\\ d_z=\underset{\lambda =1}{\overset{Q}{\Sigma }}{f}_{z,\lambda }{d}_{z,\lambda }\end{array}\right),Z=x,y---\left(9\right)$

式(9)中，f_z,(p,q)为第(p,q)分块子图像局部运动参数的加权因子，f_z,(p,q)即本区域的 优劣因子占所有优质区域的优劣因子之和的比率，d_z为最终的全局运动参数，即 Z取x时d_x为行偏移量，Z取y时d_y为列偏移量。

实施例1

以下对照数字稳像常用的块匹配法BMA、特征跟踪法FTA和传统的灰度积分投影 法IIPA，衡量本发明基于灰度积分投影互相关函数特征的运动参数估计方法IIPA-FCF 在进行低对比度高分辨率图像稳像时的有效性，从计算偏移量的快速性、准确性和鲁棒 性等几个方面进行分析。块匹配法假设已选择好待匹配的块大小和位置，图像特征法则 使用图像中的基础特征—角点。

(1)快速性比较

对工程现场实际获取的图像进行分析，图像分辨率为1920×1080，帧率为30帧/秒。 采用普通个人计算机，处理器为Intel Core i5-3210M CPU@2.50GHz，内存为4.00GB， 操作系统为Windows7，64位，显卡为GeForce GT650M。采用Matlab R2012b(8.0.0.783) 作为计算和仿真平台。

表1几种方法运行速度的对比列表

表1为几种不同方法运行在同一平台下并对实验视频进行去抖参数的计算速度对照 表。BMA方法中使用了多种设置：A-F中的目标块匹配区域为[10090](宽100像素， 高90像素)，使用全路径搜索法；A-T的目标区域不变，使用三步法搜索；B-F则调整 目标块匹配区域为[200 180]，使用全路径搜索法；B-T目标区域不变而使用三步法搜索。 由于图像分辨率很高而对比度较差，因此匹配区域事先设定好了大小和位置，这里选择 两种大小合适的区域[100 90]和[200 180]进行实验对比。

FTA方法则使用了两种不同的角点判别方法：H使用经典的HARRIS方法；L使用 区域灰度比较(LIC Local intensity comparison)快速算法。由于图像整体对比度较低， HARRIS算法中角点阈值设置非常小才能获取足够多数量的角点。IIPA法使用了全幅面 灰度积分投影。

表1可以看出，使用BMA方法且采用全路径搜索时运算非常耗时，必须采取三步 法等快速算法才能接近IIPA-FCF方法的运算速度。另外，采用BMA方法时需事先设定 好某个特征明显的区域作为待匹配目标区域。表1中使用的两个目标区域大小不同，其 耗时明显不同，匹配区域面积增大4倍，三步法匹配耗时也基本增加4倍，而全路径搜 索法耗时增加了近10倍。FTA算法的执行效率较低，使用了LIC快速算法也不能达到 预期效果。使用传统的灰度积分投影法IIPA的耗时也较长，本文提出的IIPA-FCF方法 计算效率是最优的。

(2)准确性分析

数字稳像中通常采用稳像前后连续两帧图像的峰值信噪比作为评价稳像准确度的品 质因子。它实际上是一种基于均方差值MSE(Mean square error)的比较方法，但在图像 采用单精度数、8位无符号数、16位无符号数等不同数制表示时，其计算过程略有差别。 又由于这里的图像文件和比较帧序相同，仅考虑不同算法的准确性，因此直接采用均方 差MSE作为比较参数，如式(10)所示。

$R_{\mathrm{MSE}}(I_k,I_l)=\frac{Q(I_k,I_l)}{M\times N}=\frac{1}{M\times N}\underset{\alpha =1}{\overset{M}{\Sigma }}\underset{\beta =1}{\overset{N}{\Sigma }}{[G_k(\alpha ,\beta )-G_l(\alpha ,\beta )]}^2---\left(10\right)$

式(10)中，I_k和I_l分别表示参考帧和当前帧图像，R_MSE(I_k,I_l)表示连续两帧的均方 差，Q(I_k,I_l)表示连续两帧灰度差的平方和。G_k(α,β)和G_l(α,β)的分别表示参考帧和当 前帧图像的第(α,β)点的灰度值，M和N分别表示图像的行列最大值。均方差值R_MSE反 映了图像序列的变化快慢和变化幅度，对于稳像后的相邻两幅图像，R_MSE值越小，说明 两幅图像越吻合，R_MSE为0时表示图像完全重合。实际应用中由于噪声干扰，R_MSE为0 的情况很难出现。由于均方差数值很小，一般为10^-5量级，因此这里用灰度差平方和 Q(I_k,I_l)作为对比参数。表2为上述几种方法进行同一个视频分析时的灰度差平方和列 表，图5为灰度差平方和的棒状图，Unchanged表示未进行图像去抖处理。

表2几种不同方法进行数字稳像后的图像灰度差平方和列表

表2和图5中可以看出，BMA方法和IIPA方法的准确性相当，而FTA方法的准确 性最高。BMA方法应用时，准确性不仅与匹配区域大小有关，还与匹配区域的位置有 关，位置不合适时大区域的匹配准确性反而不高，如图5所示，BMA-B-F匹配区域最 大4倍，但准确性反而降低。另外，BMA方法在计算偏移量时使用全路径搜索法和三 步法具有同等的精度。FTA方法与其它方法相比优势明显。实验中使用的两种特征点— —Harris角点和LIC角点的准确性相当。IIPA计算方法的准确性与BMA方法相当，而 本文提出的基于灰度积分投影互相关函数特征的运动参数估计方法IIPA-FCF与IIPA相 比，精确性提高约1.4倍，也明显优于BMA方法，但离精度较高的FTA相比还有明显 的不足。

(3)鲁棒性分析

鲁棒性也是数字稳像算法中需着重考虑的性能指标。由于本文的主要应用对象是海 域低对比度图像，因此鲁棒性分析主要考虑图像噪声对稳像精度的影响。在实验中人为 对视频录像帧中添加椒盐噪声，观察不同算法下的精度影响。

表3不同稳像算法中加入噪声的匹配精度

表3为不同稳像算法中加入椒盐噪声后进行稳像输出的图像灰度差平方和，与无噪 声的数字稳像相比精度有所下降，尤其是FTA方法，已无法实现正常稳像。图像受噪 声污染对BMA方法的精确性影响较小，甚至在三步法求解时精度略有提高。本文提出 的基于灰度积分投影互相关函数特征的运动参数估计方法IIPA-FCF与传统的灰度积分 投影法IIPA的精度受影响情况基本一致。