一种搜跟系统中航拍视频的快速自动拼接技术

摘要

本发明提出一种搜跟系统中航拍视频的快速自动拼接技术。首先将航拍视频序列中连续两帧彩色图像转化成灰度图像；利用图像三值化方法对灰度图像进行三值化处理、利用边缘检测sobel算子提取边缘特征信息图像、使用Harris算法提取特征点；利用光流金字塔算法完成特征点匹配；结合DLT算法和RANSAC算法计算彩色图像的最优单应性矩阵；利用修正几何转换矩阵对彩色图像进行图像修正后完成图像全景的拼接；在图像全景拼接过程中动态调整图像三值化的上下限阈值以及sobel算子的边缘检测的阈值。本发明提高了图像拼接的准确度和稳健性。

说明书

技术领域

本发明属于图像自动拼接技术领域，具体涉及一种搜跟系统中航拍视频的快速自动拼接技术。 

背景技术

图像拼接中最关键的步骤是图像特征点的寻找及匹配。杨涛等人（基于场景复杂度与不变特征的航拍视频实时配准算法.电子学报，2010.5，5期38卷）提出了利用Harris角点和SIFT描述子相结合的方式进行图像匹配，保证了算法具有良好的旋转、尺度、亮度不变性和配准精度使算法。但由于前期特征点提取的数量较多，导致后期计算量过大，无法真正意义上实现图像拼接的实时性。针对这个问题，需要从简化特征点的提取与匹配运算量入手优化算法。 

很多人提出了多种图像自动拼接技术的实现方法，主要有基于特征点的方法、基于光流的方法，基于图像融合的方法等等。这些方法主要是利用图像的几何变换矩阵完成对图像的修正。这种方法对于帧数较少的图像序列来说，拼接效果很好，但随着视频序列帧数的增加，图像拼接的几何变换矩阵就会因为前期不断地叠加而产生误差的累积与放大；同时，由于飞机的抖动和场景变换明显，也会导致图像拼接的准确率降低，同时还会导致图像拼接过程的中断 

发明内容

本发明公开了一种搜跟系统中航拍视频的快速自动拼接技术，提高了图像拼接的准确度和稳健性。 

为了解决上述技术问题，本发明提出一种搜跟系统中航拍视频的快速自动拼接技术，包括以下步骤： 

步骤一：将航拍视频序列中连续两帧彩色图像I和J转化成灰度图像I_g和J_g，并记录帧序号Z和Z+1； 

步骤二：利用图像三值化方法对灰度图像I_g进行三值化处理获得三值化图像，利用边缘检测sobel算子提取三值化图像的边缘特征信息获得边缘特征信息图像I_f，使用 Harris算法提取边缘特征信息图像I_f的特征点，在提取三值化图像边缘特征信息的过程中，设定sobel算子的边缘检测阈值为S₁和S₂，S₁代表水平方向的边缘计算参数阈值，S₂代表垂直方向的边缘计算参数阈值； 

步骤三：利用光流金字塔（KLT）算法计算光流，完成灰度图像I_g和J_g的特征点匹配； 

步骤四：结合DLT算法和RANSAC算法计算得出彩色图像I和J的最优单应性矩阵H_I,J和图像拼接的修正矩阵H'_z，并求解出图像拼接的修正几何转换矩阵P； 

步骤五：利用修正几何转换矩阵P对彩色图像J进行图像修正得到图像J'，设置一个长宽为彩色图像J长宽两倍的黑色背景，将修正后的图像J'的像素信息J'(x',y')贴至黑色背景中，完成图像全景的拼接； 

步骤六：将当前帧图像的单应性矩阵H_I,J与已保存的前帧图像的单应性矩阵H_I-1,I比较，将准确匹配的特征点对数目N与特征点对数目阈值N_f进行比较，根据比较结果动态调整图像三值化的上下限阈值T₁和T₂以及sobel算子的边缘检测的阈值S₁和S₂。 

本发明与现有技术相比，其显著优点在于，采用图像三值化、边缘提取和Harris角点检测三者结合的方法来降低特征点的数量，并筛选出图像中最为稳定的、特征最为突出的特征点；采用了RANSAC算法与DLT算法结合的方法，求解出带有修正矩阵的图像几何转换矩阵矩阵，提高图像拼接的准确度；阈值动态调整有效解决了因图像场景变换明显或飞机抖动等原因造成的程序中断问题。 

附图说明

图1是本发明简要流程图。 

图2是本发明详细流程图。 

图3（a）是仅利用Harris算法提取特征点效果图；图3（b）是本发明步骤二中经过图像三值化和sobel算子边缘检测后，再利用Harris算法提取特征点的效果图。 

图4（a）和图4（b）分别为本发明实施例自动拼接70帧和150帧的总场景图。 

具体实施方式

如图1和图2所示，本发明一种搜跟系统中航拍视频的快速自动拼接技术，包括以 下步骤： 

步骤一：将航拍视频序列中连续两帧彩色图像I和J转化成灰度图像I_g和J_g，并记录帧序号Z和Z+1。 

图像灰度化可以将R,G,B任何一个分量作为图像的灰度值，也可以用R,G,B三者的平均值作为灰度值，也可以利用亮度公式将R,G,B分量转换成灰度值。本发明采用亮度公式的方法，将彩色图像I和J转换成只包含亮度信息的灰度图像I_g和J_g，亮度公式为公式（1）： 

Y＝0.3R+0.59G+0.11B    （1） 

步骤二：利用图像三值化方法对灰度图像I_g进行三值化处理获得三值化图像，利用边缘检测sobel算子提取三值化图像的边缘特征信息获得边缘特征信息图像I_f，使用Harris算法提取边缘特征信息图像I_f的特征点，在提取三值化图像边缘特征信息的过程中，设定sobel算子的边缘检测阈值为S₁和S₂，S₁代表水平方向的边缘计算参数阈值，S₂代表垂直方向的边缘计算参数阈值； 

所述获得三值化图像的计算方法如公式（2）所示， 

$>\left(\begin{array}{c}p[x,y]=0,p[x,y]<T_1\\ p[x,y],T_1<p[x,y]<T_2\\ p[x,y]=1,p[x,y]>T_2\end{array}\right)---\left(2\right)$>

式（2）中，p[x,y]是灰度图像I_g中坐标为(x,y)处的像素值，T₁和T₂为图像三值化的上下限阈值，T₁和T₂的具体取值根据航拍视频的场景进行设置，例如公路、河流和岩石等的灰度图像素值一般大于0.6，树木、植被等的灰度图像素值一般小于0.4； 

所述使用Harris算子提取边缘特征图像I_f的特征点的计算过程如下： 

1.1计算像素点p[x,y]的相关矩阵M，计算方法如公式（3）所示， 

$>M=\underset{x,y}{\Sigma }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)\otimes \left(\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right)---\left(3\right)$>

式（3）中，w(x,y)是窗口大小为w*w的高斯权重函数，I_x代表行方向的导数，I_y代表列方向的导数； 

1.2计算像素点p[x,y]的特征点响应R(x,y)，计算方法如公式（4）所示， 

R(x,y)＝detM-k′*(traceM)²    （4） 

式（4）中，detM代表相关矩阵M行列式的值，traceM代表相关矩阵M的迹，k′为常数，一般为0.04-0.06； 

1.3根据公式（5）判断像素点p[x,y]是否为特征点，如果像素点p[x,y]的特征点响应R(x,y)满足公式（5），则像素点p[x,y]为特征点，将像素点p[x,y]放入特征点集合D中， 

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

式（5）中，R(x',y')代表窗口w*w中除像素点p[x,y]之外的其它像素点的特征点响应，R_max为边缘特征信息图像I_f中的最大特征点响应； 

为了实现算法的实时性，需要降低程序的运算量，最根本的方法是在保证特征点质量的基础上降低特征点的数目，根据航拍视频背景范围大、空旷、颜色单一，并以高山、河流和公路、树木和房屋等轮廓性较强的事物居多的特点。本发明采取的图像三值化、图像边缘提取和Harris角点提取三者结合的方式可以去除图像冗余信息，在突出图像特征的基础上进行特征点的有效选取，有利于降低特征点数量，从而降低后续计算量。 

步骤三：利用光流金字塔（KLT）算法计算光流，完成灰度图像I_g和J_g的特征点匹配，计算过程如下： 

3.1利用下采样公式分别对灰度图I_g和J_g进行下采样，获得灰度图I_g和J_g的图像金字塔，计算方法如公式（6）所示， 

$>\left(\begin{array}{c}{I_g}^L(x,y)=\frac{1}{4}{I_g}^{L-1}(2x,2y)+\\ \frac{1}{8}({I_g}^{L-1}(2x-\mathrm{1,2}y)+{I_g}^{L-1}(2x+\mathrm{1,2}y)+{I_g}^{L-1}(2x,2y-1)+{I_g}^{L-1}(2x,2y+1))+\\ \frac{1}{16}({I_g}^{L-1}(2x-\mathrm{1,2}y-1)+{I_g}^{L-1}(2x+\mathrm{1,2}y+1)+{I_g}^{L-1}(2x-\mathrm{1,2}y+1)+{I_g}^{L-1}(2x+\mathrm{1,2}y+1))\end{array}\right)---\left(6\right)$>

式（6）中，(x,y)表示灰度图I_g中第L层图像中像素坐标为（x，y）的像素值，设原灰度图I_g和J_g为图像金字塔的第一层图像，根据下采样公式可以得到图像金字塔的第二层图像和同理，对第二层图像和进行下采样，得到图像金字塔的第三层图像和一般情况下，图像金字塔的层数取3到5比较合适。 

3.2计算图像中的特征点p[x,y]相对于图像的光流量，计算方法如公式（7）： 

$>\left(\begin{array}{cc}\underset{w*w}{\Sigma }{I}_x{I}_x& \underset{w*w}{\Sigma }{I}_x{I}_y\\ \underset{w*w}{\Sigma }{I}_x{I}_y& \underset{w*w}{\Sigma }{I}_y{I}_y\end{array}\right)d^n=-\left(\begin{array}{c}\underset{w*w}{\Sigma }{I}_x{I}_t\\ \underset{w*w}{\Sigma }{I}_y{I}_t\end{array}\right)---\left(7\right)$>

式（7）中，I_x表示w*w窗口内每个像素点的行方向的导数，I_y表示w*w的窗口内每个像素点的列方向的导数，I_t表示图像与的时间导数，dⁿ＝[uⁿ,vⁿ]表示图像中特征点p[x,y]相对于图像的光流量，其中uⁿ为图像金字塔第n层的行方向的光流量，vⁿ为图像金字塔第n层的列方向的光流量; 

光流计算一般从顶层开始，如果图像金字塔的层数为3，这里便从第三层开始; 

3.3将光流量dⁿ作为初始光流量传递至光流金字塔的第n-1层，即在图像中的像素点q[x,y]上叠加光流量dⁿ，使每一个像素点的坐标变q[x+2*uⁿ,y+2*vⁿ]，然后在图像和叠加了光流量dⁿ的图像之间，重复步骤3.2计算出图像金字塔第n-1层的光流量d^n-1＝[u^n-1,v^n-1]; 

3.4重复步骤3.2和3.3，直至计算获得图像金字塔第一层的光流量d¹＝[u¹,v¹]，即为灰度图像I_g和J_g的总光流量； 

3.5根据第一层光流量d¹＝[u¹,v¹]和灰度图像I_g的特征点p[x,y]坐标，根据公式（8）计算获得灰度图像J_g中的匹配特征点的坐标(x',y')，若求得的x'和y'超出了图像的像素范围，即特征点跟踪匹配失败，则将匹配特征点q[x',y']舍掉，并将对应的特征点p[x,y]从特征点集合D中剔除出去；若求得的x'和y'在图像像素范围内，则将匹配特征点q[x',y']存入特征点集D'中， 

x'＝x-u¹    y'＝y-v¹   （8） 

步骤四：结合DLT算法和RANSAC算法计算得出彩色图像I和J的最优单应性矩阵H_I,J和图像拼接的修正矩阵H'_z，并求解出图像拼接的修正几何转换矩阵P， 

所述最优单应性矩阵H_I,J的计算过程如下： 

4.1利用公式（9）判断每一对特征点对是否为准确匹配的特征点对，将误差较大的特征点对剔除，将准确匹配的特征点对存入特征点对集合S中，并统计出准确匹配的特征点对数N； 

Dpq≥|p[x_i,y_i]-q[x'_i,y'_i]|      （9） 

式（9）中，Dpq为特征点对的平均距离，计算方式如公式（10）所示， 

$>\mathrm{Dpq}=\underset{i=1}{\overset{N}{\Sigma }}\left|p\right[x_i,y_i]-q[{x^{\prime }}_i,{y^{\prime }}_i\left]\right|/N---\left(10\right)$>

4.2利用DLT算法计算获得图像I和图像J的初始单应性矩阵H'_I,J，具体实现方法如下所述： 

在特征点对集合S中任意提取四对特征点对p[x_i,y_i]、_q[x'_i,y'_i]，i=1,2,3,4，构造DLT计算矩阵A，如公式（11）所示， 

$>A=\left(\begin{array}{ccccccccc}p\left[x_1\right]& 0& -p\left[x_1\right]*q\left[{x^{\prime }}_1\right]& p\left[y_1\right]& 0& -p\left[y_1\right]*q\left[{x^{\prime }}_1\right]& 1& 0& -q\left[{x^{\prime }}_1\right]\\ 0& p\left[x_1\right]& -p\left[x_1\right]*q\left[{y^{\prime }}_1\right]& 0& p\left[y_1\right]& -p\left[y_1\right]*q\left[{y^{\prime }}_1\right]& 0& 1& -q\left[{y^{\prime }}_1\right]\\ p\left[x_2\right]& 0& -p\left[x_2\right]*q\left[{x^{\prime }}_2\right]& p\left[y_2\right]& 0& -p\left[y_2\right]*q\left[{x^{\prime }}_2\right]& 1& 0& -q\left[{x^{\prime }}_2\right]\\ 0& p\left[x_2\right]& -p\left[x_2\right]*q\left[{y^{\prime }}_2\right]& 0& p\left[y_2\right]& -p\left[y_2\right]*q\left[{y^{\prime }}_2\right]& 0& 1& -q\left[{y^{\prime }}_2\right]\\ p\left[x_3\right]& 0& -p\left[x_3\right]*q\left[{x^{\prime }}_3\right]& p\left[y_3\right]& 0& -p\left[y_3\right]*q\left[{x^{\prime }}_3\right]& 1& 0& -q\left[{x^{\prime }}_3\right]\\ 0& p\left[x_3\right]& -p\left[x_3\right]*q\left[{y^{\prime }}_3\right]& 0& p\left[y_3\right]& -p\left[y_3\right]*q\left[{y^{\prime }}_3\right]& 0& 1& -q\left[{y^{\prime }}_3\right]\\ p\left[x_4\right]& 0& -p\left[x_4\right]*q\left[{x^{\prime }}_4\right]& p\left[y_4\right]& 0& -p\left[y_4\right]*q\left[{x^{\prime }}_4\right]& 1& 0& -q\left[{x^{\prime }}_4\right]\\ 0& p\left[x_4\right]& -p\left[x_4\right]*q\left[{y^{\prime }}_4\right]& 0& p\left[y_4\right]& -p\left[y_4\right]*q\left[{y^{\prime }}_4\right]& 0& 1& -q\left[{y^{\prime }}_4\right]\end{array}\right)---\left(11\right)$>

对矩阵A进行奇异值分解，得到复杂形式的单应性矩阵H_(I,J)″,将H_(I,J)″进行简化获得简化后的初始单应性矩阵H_(I,J)'，奇异值分解计算过方式如公式（12）所示， 

[U,S,V]＝SVD(A),H_(I,J)″＝V(:,end)    （12） 

简化后的初始单应性矩阵H'_I,J为公式（13）： 

$>H_{{I,J}^{\prime }}=\left(\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\end{array}\right),$>其中h₃₃＝1    （13） 

4.3利用RANSAC算法优化初始单应性矩阵H_(I,J)'，求得最优单应性矩阵H_I,J，具体实现方法如下所述： 

将图像I的特征点集合D中的特征点p[x,y]带入公式（14）中，求得中间量q'[x'_i,y'_i]： 

q'[x'_i,y'_i]＝(H_(I,J)')^-1*p[x_i,y_i]  （14） 

若求得的中间量_q'[x'_i,y'_i]满足公式（15），则将该特征点对保留在集S中；否则，将 该特征点对放入余集C中，最终计算出保留在集S中的匹配点对的数量Num。 

|q[x'_i,y'_i]-q'[x'_i,y'_i]|≤0.1   （15） 

若Num≥M，则结束RANSAC算法，并将当前初始单应性矩阵H_(I,J)'作为最优单应性矩阵H_I,J，其中M为特征点对集S的样本个数阈值； 

若Num＜M，则进一步判断RANSAC算法的循环次数是否在预先设定的阈值N_RAN内，若RANSAC算法的循环次数在阈值N_RAN内，则重复步骤4.2求得新的初始单应性矩H_(I,J)'；若RANSAC算法的循环次数超出阈值N_RAN，则结束RANSAC算法，并将当前初始单应性矩阵H_(I,J)'作为最优单应性矩阵H_I,J，并保存H_I,J。 

所述图像拼接的修正矩阵H'_Z采取隔帧计算的方法，即对当前帧序号进行如式（16）的判断，若图像帧序号不满足公式（16），则直接返回步骤一进行下帧图像的图像拼接操作；若图像帧序号满公式(16)，则进行图像拼接修正矩阵H'_Z的计算，其中Z代表当前帧图像的帧序号： 

Z＝k*n+1，(k＝1,2,3...,n＝1,2,3...)   （16） 

式（16）中，n是自然数，k代表隔帧的数目，具体设置需要考虑航拍视频序列的图像重叠情况，若图像间图像重叠区域较大，k的值可以适当取大一些，若图像间图像重叠区域相对较小，k的值可以适当取小一些，一般为3到7比较合适；图像拼接修正矩阵H'_Z的计算过程为：取出预先存储的第Z-k帧彩色图像J_save和当前帧彩色图像J,重复最优单应性矩阵的计算过程计算获得最优单应性矩阵H_Z-k,Z’将最优单应性矩阵H_Z-k,Z记为修正矩阵H'_Z,计算结束后保存当前帧彩色图像J，将其作为下次修正矩阵计算的J_save。 

所述图像拼接的修正几何转换矩阵P的计算方法如公式（17）所示， 

P＝H_1,2*H_2,3*…^H'_k+1*…H_J-1,J*^H_I,J*H'_Z   （17） 

步骤五、根据公式（18）利用修正几何转换矩阵P对彩色图像J进行图像修正得到图像J'，设置一个长宽为彩色图像J长宽两倍的黑色背景，将修正后的图像J'的像素信息J'(x',y')贴至黑色背景中，完成图像全景的拼接， 

$>\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array}\right)=P*\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)---\left(18\right)$>

式（18）中，(x,y)为原彩色图像J的像素坐标，(x',y')为转换后图像J'的像素坐标。 

步骤六：按照公式（19）或公式（20），将当前帧图像的单应性矩阵H_I,J与已保存的前帧图像的单应性矩阵H_I-1,I,比较，将准确匹配的特征点对数目N与特征点对数目阈值N_f进行比较， 

若公式（19）和公式（20）有任意一个公式成立的，则将图像三值化的上下限阈值T₁和T₂调整为T₁-0.1和T₂+0.1，将sobel算子的边缘检测的阈值S₁和S₂调整为S₁-0.1和S₂-0.1，用于下帧图像的处理， 

若公式（19）和公式（20）都不成立的，则保持图像三值化的上下限阈值T₁和T₂以及sobel算子的边缘检测的阈值S₁和S₂不变， 

|n_ij-n'_ij|≥1.0，i＝1,2,3，j＝1,2,3     （19） 

N≤N_f       （20） 

式（19）中，n_ij为矩阵H_I,J中第i行第j列的参数，n'_ij为矩阵H_I-1,I中第i行第j列的参数；式（20）中N_f为准确匹配的特征点对数目阈值。 

实施例 

本实施例利用basler相机采集的24位真彩色图像，图像大小为640*480。为了使获取的真彩色图像序列相互之间有足够的重叠信息，采集帧频设置为10帧/s。 

本实施例的视频序列背景多为植被和公路且轮廓行非常强，根据视频特点，设置图像三值化的阈值为[0.40.65]，sobel算子的边缘检测的阈值为[0.20.3]，在此基础上进行Harris角点检测。为了在保证算法准确率的前提下降低算法的计算量，Harris角点检测的窗口选取为5*5，,图像金字塔的层数选取为3，KLT光流计算的窗口大小设置为7*7，图像拼接的修正矩阵隔帧数目为k=4。 

图3（b）是使用本发明步骤二利用图像三值化方法、边缘检测算子和Harris算法三者结合的方法获得的特征点。图3（a）为不进行图像三值化和边缘检测而直接进行Harris 角点检测的结果。将图3（b）与图3（a）相比较可以明显的看出，图3（b）中特征点数量明显下降，且大部分集中在背景边缘部分，因此特征点提取的稳定性和效率均得到了提高。 

图4（a）和图4（b）是使用本发明方法获得的最终效果图，图4（a）是连续处理70帧图象所得的效果图，图4（b）是连续处理150帧图象所得的效果图。从图4（a）和图4（b）图中可以看出场景旋转明显，在图像三值化和边缘检测阈值的自动控制下，图像很好地完成了自动拼接，有效地避免了飞机气流和场景噪声的影响。同时，根据对本发明方法算法量的统计，本发明方法的实时性可以支持15帧以上的相机实时拍摄帧频，很好地实现了航拍视频的图像快速自动拼接。