CN103886569A - 连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法 - Google Patents

Abstract

本发明公开了一种用于连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法，属于无人机图像处理技术领域。该发明深入分析图像拼接的四个阶段（特征提取、特征匹配、变换矩阵、图像合成）的数据量的大小以及阶段间的串行或并行的依存关系。依据数据的类型及其关系，设计符合无人机大数据量图像处理特点和要求的数据处理的模式。同时将计算机的多核资源、多线程计算方法应用到数据处理过程中以提高计算资源的利用率，加快运算的速度。该方法针对性强，能够显著提高无人机大数据量遥感图像的拼接处理能力，保障无人机图像拼接情报的实时生成。

Description

连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法

技术领域

本发明属于无人机图像处理技术领域，具体涉及一种用连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法。

背景技术

无人机侦察技术属于遥感技术的一种，相对卫星侦察具有其自身优点，如成本低、侦察地域控制灵活、不存在访问时间和周期限制、地面目标分辨率高等；相对有人侦察机来说，有昼夜可持续工作、不考虑飞行员疲劳和伤亡等方面的优越性。近年来，无人机因其低空遥感数据具有高分辨率、高灵活性、高效率和低成本的优势而被广泛应用于自然灾害区域评估、战场侦察、环境监测等领域。例如汶川地震、玉树地震发生时，利用无人机拍摄的遥感图像对灾区的灾害损失快速评估及制定救灾决策提供科学依据；美军大战伊拉克战场上使用了数十架无人机对战场区域进行侦查、监视，为美军制定实时的作战计划提供了有利的信息。

无人机侦察技术虽然在自然灾害区域评估、战场侦察等领域起到了重要的作用，但是由于无人机获得的侦察图像具有数据量大、重叠率高的特点，很难用一帧图像将感兴趣区域的信息全部展现出来，使得感兴趣区域目标的获取和精确定位存在巨大困难，从而给地面站指挥人员带来繁重的工作。基于以上问题，为了从无人机局部区域侦察图像中获得更加有效、全面、精准的信息，需要对侦察图像进行一系列合理的处理，实现无人机侦察图像的重建工作。通过无人机侦察图像的重建，能够得到高分辨率、大视角、宽视野的全景图像，从而为地面站指挥人员提供一种有效且完整的场景表示方法，使其能够更好地统一处理、解译、分析和研究图像信息。

目前，国内外在图像重建领域已经有了相当成熟的理论基础，图像拼接作为图像重建的重要研究方向，已广泛应用于实际生活当中。2004年，Lowe提出一种尺度不变特征提取SIFT（Scale Invariant Feature Transform）算法，是目前来说应用最多的图像配准算法。该算法虽然给图像配准带来良好的配准效果，但是在提高图像拼接速度方面并没有太大改善，大多数学者都是从SIFT算法的参数配置方面来改善速度指标，这样并不能完全达到实时处理的效果。对于大数据量的无人机侦察图像来说，想要采用改变SIFT算法参数的方法来达到图像实时处理的效果难度更大。

发明内容

本发明的目的是为了提高无人机大数据量侦察图像的处理速度，根据图像拼接的处理流程（特征提取、特征匹配、变换模型、图像合成等方面）分析拼接过程中各个阶段待处理数据量的大小以及各阶段的前后依附关系，利用计算机多核资源，创建多个线程，设计图像拼接的加速方法，实现无人机大数据遥感图像的快速拼接。

本发明提供的连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法，具体步骤如下：

第一步，分配任务块，创建线程，提取图像SIFT特征点。

第二步，对第一步提取的图像特征点进行匹配，找出相邻两幅图像之间的同名点。

第三步，根据第二步匹配好的同名特征点计算相邻两幅图像之间的变换矩阵。

第四步，利用第三步获得的变换矩阵H将待拼接的图像合成为高分辨率的全景图像。

本发明的优点在于：

（1）基于拼接流程设计并行处理方法，处理速度快；

（2）采用三步匹配策略，使得特征点配准精度高；

（3）充分利用计算机多核资源，提高了计算机CPU的使用效率；

因此，使用本发明方法能充分利用计算机自身的资源，提高大数据量图像拼接的处理速度，实现图像特征点的精确配准，从而得到拼接效果优良的高分辨率全景图像。

附图说明

图1是本发明用于无人机大数据量侦察图像的并行拼接方法的整体步骤图；

图2是本发明步骤一中任务分块、线程创建并行方法设计示意图；

图3是本发明步骤二中图像特征点三步配准法实现流程图；

图4a、图4b为单幅图像提取的特征点示意图；

图5a为图4a、图4b两幅图像的粗匹配结果；

图5b为图4a、图4b两幅图像中存在的误匹配；

图6为图4a、图4b采用斜率约束法去除误匹配的匹配结果；

图7为实例中20幅图像合成的全景图像。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法，流程如图1所示，包括以下几个步骤：

第一步，分配任务块，创建线程，提取图像SIFT特征点。

在图像拼接的处理流程中，图像的特征提取是图像拼接过程中最耗时的一步，并且相邻两幅图像的特征点之间不存在依附关系，满足并行化处理条件。因此，利用计算机处理器的多核资源，实现图像特征提取的并行化处理，如图2所示。

包括以下几个步骤：

首先，读入待拼接图像的存储路径，统计出待拼接图像的数目n；

然后，根据图像的数目n和处理器的核心数k将待拼接的图像进行任务分块，即将待处理的图像平均分为k组，理想情况下每组的图像数目应该相等，为n/k幅图像，形成k个任务块。

最后，根据分配的任务块数目创建线程，所创建的k个线程同时并行工作，分别用来提取相应任务块中图像的SIFT特征点，其中，利用尺度不变特征提取算法（SIFT）来提取待拼接图像的特征点。

假设每个线程完成任务块的时间为T_i，i＝1,2,3……k，则采用多线程并行计算n幅图像特征点的时间应该为T_并＝max（T₁,T₂,……T_k）。相比于串行处理时所耗时间T_串＝T₁+T₂+……+T_k要缩短很多，同时，计算机处理器核心数目越多，可创建的线程数目越多，处理时间就越短。

第二步，对第一步提取的图像特征点进行匹配，找出相邻两幅图像之间的同名点。

确定稳定的特征点是图像配准和拼接的重要环节，为了使相邻图像之间的特征点达到理想的配准精度，本发明采用三歩配准法逐步消除相邻图像之间的误匹配，如图3所示。具体实现方案如下所述：

首先，通过相邻两幅图像特征点的描述符（128维向量）计算欧氏距离实现特征点的粗匹配。

将待匹配的第2幅图像提取的所有特征点（记为P2_i,i=1,2,3,……N2,N2表示第二幅图像的特征点的数目)构建成k-d树，逐个选定待匹配的第1幅图像的特征点（记为P1_j，j从1到N1，N1表示第一幅图像的特征点数目）基于图像2构建的k-d树进行搜寻，找到与P1_j对应的128维特征点描述符之间欧式距离最近的特征点和次近的特征点，并计算它们之间的欧氏距离分别标记为d₀和d₁，如果最近欧氏距离d₀与次近邻欧氏距离d₁的比值小于某个阈值（根据经验值可以设置为0.8），则认为能够计算出最近邻欧氏距离的两个特征点匹配成功，否则，遍历寻找图像1中下一个特征点的同名点，直至将图像1中的特征点遍历完毕，获得N对相匹配的同名点。

然后，利用斜率约束删除斜率相差较大的匹配点对。

将0-180度的角度范围均分为36等分，每5度为一个级别；分别计算每对同名点对之间的斜率，并将其对应到每个角度级别，对处于（-90°，0°)的负角度值，则将对应的角度值加上180°后将其变换到90°-180°的范围之内；统计所有角度级别的值，找出处于角度级别最多的对应的角度，并将该角度对应的斜率作为剔除阈值k_m，删除斜率对应角度与k_m对应角度之间的差值相差大于5°的那些特征点对。

最后，结合第三步，通过RANSAC算法获取的变换矩阵，将相邻两幅图像中一幅图像的同名特征点通过变换矩阵变换到另一幅图像坐标系下，若变换后的点与其对应的同名点的欧式距离大于某个阈值，则认为这对同名点为外点，将该对同名点剔除，来实现特征点精确匹配。

第三步，根据第二步匹配好的同名特征点计算相邻两幅图像之间的变换矩阵。

首先，采用RANSAC算法确定初始变换矩阵。

令相邻两幅图像上的同名点坐标分别为I_i和I_i+1，I_i＝(x_i,y_i,1)^T，I_i+1＝(x_i+1,y_i+1,1)^T，其中I_i,I_i+1为相邻两幅图像上对应的同名点坐标。两幅图像同名点之间的对应关系可以表示如下：

I_i＝HI_i+1            （1）

其中：

$H=\left(\begin{array}{ccc}{h}_0& h_1& h_2\\ h_3& h_4& h_5\\ h_6& h_7& 1\end{array}\right)---\left(2\right)$

其中H为变换矩阵的表达式，h₀,h₁,h₃,h₄为图像缩放和图像旋转量，h₂,h₅为水平方向和垂直方向的位移量，h₆,h₇为水平方向和垂直方向的变形量。

变换矩阵H需要求解8个参数，所以只要4个同名点对即可。RANSAC算法是采用迭代的方法求解初始变换矩阵，求解之前设置好判定为内点的误差阈值、迭代次数、内点数目阈值，每次迭代都随机抽取4个同名点计算出一个变换矩阵H，利用式1将第二幅图像的特征点变换到第一幅图像坐标系下，计算变换后同名点之间的变化误差，如果在设定的内点误差阈值内，则认为该对同名点为内点，统计所有内点的数目，看是否满足设定的内点数目阈值，满足的话将此次求得的H作为初始变换矩阵，否则，继续迭代。

然后，采用非线性最小二乘法L-M算法，在初始变换矩阵的基础上求解高精度的变换矩阵，以便减小累积误差和噪声的影响，提高拼接质量。根据最小二乘法L-M算法求解变换矩阵H的最优参数，使得以变换矩阵参数为变量构造的函数F(h)取最小值，函数F(h)表示所有匹配同名特征点之间距离和，如式3：

$F\left(h\right)=\underset{i}{\mathrm{\Σ}}{[{(\frac{h_0{x}_i+h_1{y}_i+h_2}{h_6{x}_i+h_7{y}_i+1}-x_i^{\′})}^2+{(\frac{h_3{x}_i+h_4{y}_i+h_5}{h_6{x}_i+h_7{y}_i+1}-y_i^{\′})}^2]}^{\frac{1}{2}}---\left(3\right)$

其中x_i,y_i表示相邻两幅图像中其中一幅图像特征点的图像坐标，x′_i,y′_i表示另一幅与之相匹配的特征点的图像坐标，令

$f_j\left(h\right)={[{(\frac{h_0{x}_i+h_1{y}_i+h_2}{h_6{x}_i+h_7{y}_i+1}-x_i^{\′})}^2+{(\frac{h_3{x}_i+h_4{y}_i+h_5}{h_6{x}_i+h_7{y}_i+1}-y_i^{\′})}^2]}^{\frac{1}{2}}---\left(4\right)$

其中f_i(h)表示公式（3）后半部分的求和分量，即

在给定精度参数ε条件下，采用迭代的方法最小化变换矩阵H，迭代过程如式5

h_k+1＝h_k+h_LM             （5）

其中：h_k+1,h_k表示变换矩阵，k从0开始，h₀表示初始的变换矩阵，h_LM是一表达式如式6所示：

h_LM＝-(J(h)^TJ(h)+μI)^-1J(h)^Tf(h)            （6）

其中：μ为一大于0的常数，一般设置成μ＝0.01，式中的h表示变量集，代表变换矩阵h₀,h₁......h₇8个参数，J(h)为f对矩阵参数h₀,h₁......h₇偏微分的雅克比矩阵，如式7：

$J\left(h\right)=\left|\begin{array}{cccc}\frac{\∂f_1}{{\∂h}_0}& \frac{{\∂f}_1}{{\∂h}_2}& ...& \frac{{\∂f}_1}{{\∂h}_7}\\ \frac{{\∂f}_2}{{\∂h}_0}& \frac{{\∂f}_2}{{\∂h}_1}& ...& \frac{{\∂f}_2}{{\∂h}_7}\\ ...& ...& ...& ...\\ \frac{{\∂f}_n}{{\∂h}_0}& \frac{{\∂f}_n}{{\∂h}_1}& ...& \frac{{\∂f}_n}{{\∂h}_7}\end{array}\right|---\left(7\right)$

其中：表示多元函数f_n(h)对各个自变量h₀,h₁......h₇的偏微分，n从0到K，K表示迭代次数。

经过K次（K大于等于3）迭代运算，可以获取精度相对较高的变换矩阵H＝h_k+1，采用这个矩阵进行图像的拼接，减少了图像拼接过程中累积误差和噪声的影响，可以提高图像拼接质量。

第四步，利用第三步获得的变换矩阵H将待拼接的图像合成为高分辨率的全景图像。

利用已获得的变换矩阵H，将单帧图像按照原有的成像顺序逐一变换到同一坐标系下，然后采用最近邻插值法或双线性插值法将图像重新赋值合成在全景图像上。为了使最终的全景图像具有较为平滑的视觉效果，重叠区域的图像融合过程采用基于小波变换的图像融合方法来实现图像的融合，对待融合的两幅图像完成小波变换的多级分解，不断使用低通滤波器和高通滤波器进行滤波，并进行逐级降2采样，最后获得多级高频分量和低频分量。然后对两幅图像每层的系数图，采用相应的融合策略完成融合操作，最后得到融合后的小波塔形结构，对它进行小波重构变换产生融合后的图像，从而有效地消除相邻图像之间的拼接缝。

实施例：

本实例采用4核计算机对88帧无人机遥感图像进行测试，测试过程及结果如下所述：

第一步：分配任务块，创建线程，提取图像SIFT特征点。

1)读入待拼接图像的存储路径，统计出待拼接图像的数目n=88；

2）根据图像的数目n和处理器的核心数k=4将待拼接的图像进行任务分块，即将待处理的图像平均分为k组，理想情况下每组的图像数目应该相等，为n/k幅图像。

3）根据分配的任务块数目创建线程，所创建的k个线程同时并行工作，分别用来提取相应任务块中图像的特征点，两幅单幅图像的特征点如图4a、4b所示。

测试过程中，根据计算机的内核数目分别开辟1、2、4、8个线程，在线程数目小于等于计算机内核数目的情况下，处理速度跟线程数目成正比，测试结果如表1所示：

表1计算时间对比

线程数目	组数	层数	计算时间（ms）
				1	2	3	66019
2	2	3	33650
				4	2	3	18018
8	2	3	18423

第二步：对第一步提取的图像特征点采用三步匹配法进行匹配，找出相邻两幅图像之间的同名点，该步骤只对88帧图像中的相邻两幅进行说明，以图4a、图4b为例：

（1）通过相邻两幅图像特征点的描述符（128维向量）计算欧氏距离实现特征点的粗匹配，匹配结果如图5所示，图5a为图4a、图4b两幅图像的粗匹配结果，图5b为图4a、图4b两幅图像中存在的误匹配。

（2）利用斜率约束法删除斜率相差较大的匹配点对，匹配结果如图6所示，图6为图4a、图4b采用斜率约束法去除误匹配的匹配结果。

（3）通过RANSAC算法剔除外点来实现特征点精确匹配。

第三步：根据第二步匹配好的同名特征点计算相邻两幅图像之间的变换矩阵，此处的变换矩阵为透视变换矩阵

$H=\left(\begin{array}{ccc}0.998569& 3.93978*{10}^{-5}& -71.290049\\ 5.070968*{10}^{-4}& 0.999078& 23.675688\\ -2.8121*{10}^{-6}& -2.2483*{10}^{-7}& 1\end{array}\right)$

第四步：利用第三步获得的变换矩阵H将待拼接的图像合成为高分辨率的全景图像，由于全景图像有一定的尺度限制，并不能无限制的合成，测试结果给出88帧中的20帧图像的拼接结果，如图7所示。

Claims (3)

1.一种连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法，包括以下几个步骤：

第一步，分配任务块，创建线程，提取图像SIFT特征点；

读入待拼接图像，将待拼接图像分为k组，形成k个任务块，创建k个线程，分别提取相应任务块中图像的SIFT特征点，k为处理器核心数；

第二步，对第一步提取的图像特征点进行匹配，找出相邻两幅图像之间的同名点；

具体包括以下几个步骤：

首先，通过相邻两幅图像特征点的描述符计算欧氏距离，实现特征点的粗匹配，其中，描述符为128维向量；

将待匹配的第2幅图像提取的所有特征点构建成k-d树，第2幅图像特征点记为P2_i,i=1,2,3,……N2，N2表示第2幅图像的特征点的数目，逐个选定待匹配的第1幅图像的特征点，第1幅图像的特征点记为P1_j，j从1到N1，N1表示第一幅图像的特征点数目，基于图像2构建的k-d树进行搜寻，找到与P1_j对应的128维特征点描述符之间欧式距离最近的特征点和次近的特征点，并计算它们之间的欧氏距离，分别标记为d₀和d₁，如果最近欧氏距离d₀与次近邻欧氏距离d₁的比值小于预设阈值，则第1幅图像的特征点与第2幅图像中欧式距离最近的特征点匹配成功，得到第1幅图像中待匹配特征点的同名点，形成一对同名点否则，遍历寻找第1幅图像中下一个特征点的同名点，直至将第1幅图像中的特征点遍历完毕，获得N对相匹配的同名点；

然后，利用斜率约束删除斜率相差大于5°匹配点对；

将0-180度的角度范围均分为36等分，每5度为一个级别；

分别计算每对同名点对之间的斜率，并将其对应到每个角度级别，对处于（-90°，0°)的负角度值，则将对应的角度值加上180°后将其变换到90°-180°的范围之内；统计所有角度级别的值，找出处于角度级别最多的对应的角度，并将该角度对应的斜率作为剔除阈值k_m，删除斜率对应角度与k_m对应角度之间的差值相差大于5°的特征点对；

最后，通过RANSAC算法剔除获取的变换矩阵，将相邻两幅图像中一幅图像的同名特征点通过变换矩阵变换到另一幅图像坐标系下，若变换后的点与其对应的同名点的欧式距离大于阈值，则认为这对同名点为外点，将该对同名点剔除，实现特征点精确匹配；

第三步，根据第二步匹配好的同名特征点计算相邻两幅图像之间的变换矩阵；

首先，采用RANSAC算法确定初始变换矩阵；

令相邻两幅图像上的同名点坐标分别为I_i和I_i+1，I_i＝(x_i,y_i,1)^T，I_i+1＝(x_i+1,y_i+1,1)^T，其中I_i,I_i+1为相邻两幅图像上对应的同名点坐标；两幅图像同名点之间的对应关系可以表示如下：

I_i＝HI_i+1               （1）

其中：

$H=\left(\begin{array}{ccc}{h}_0& h_1& h_2\\ h_3& h_4& h_5\\ h_6& h_7& 1\end{array}\right)---\left(2\right)$

其中H为变换矩阵的表达式，h₀,h₁,h₃,h₄为图像缩放和图像旋转量，h₂,h₅为水平方向和垂直方向的位移量，h₆,h₇为水平方向和垂直方向的变形量；

设置RANSAC算法中判定为内点的误差阈值、迭代次数、内点数目阈值，每次迭代时随机抽取4个同名点，计算出一个变换矩阵H，利用式1将第二幅图像的特征点变换到第一幅图像坐标系下，计算变换后同名点之间的变化误差，如果在设定的内点误差阈值内，则认为该对同名点为内点，统计所有内点的数目，看是否满足设定的内点数目阈值，满足的话将此次求得的H作为初始变换矩阵，否则，继续迭代；

然后，采用非线性最小二乘法L-M算法，在初始变换矩阵的基础上求解高精度的变换矩阵，根据最小二乘法L-M算法求解变换矩阵H的最优参数，使得以变换矩阵参数为变量构造的函数F(h)取最小值，函数F(h)表示所有匹配同名特征点之间距离和，如式3：

$F\left(h\right)=\underset{i}{\mathrm{\Σ}}{[{(\frac{h_0{x}_i+h_1{y}_i+h_2}{h_6{x}_i+h_7{y}_i+1}-x_i^{\′})}^2+{(\frac{h_3{x}_i+h_4{y}_i+h_5}{h_6{x}_i+h_7{y}_i+1}-y_i^{\′})}^2]}^{\frac{1}{2}}---\left(3\right)$

其中：x_i,y_i表示相邻两幅图像中其中一幅图像特征点的图像坐标，x′_i,y′_i表示另一幅与之相匹配的特征点的图像坐标，令

$f_j\left(h\right)={[{(\frac{h_0{x}_i+h_1{y}_i+h_2}{h_6{x}_i+h_7{y}_i+1}-x_i^{\′})}^2+{(\frac{h_3{x}_i+h_4{y}_i+h_5}{h_6{x}_i+h_7{y}_i+1}-y_i^{\′})}^2]}^{\frac{1}{2}}---\left(4\right)$

其中：f_i(h)表示公式（3）后半部分的求和分量，即

在给定精度参数ε条件下，采用迭代的方法最小化变换矩阵H，迭代过程如式5

h_k+1＝h_k+h_LM             （5）

其中h_k+1,h_k表示变换矩阵，k从0开始，h₀表示初始的变换矩阵，h_LM表示如式6：

h_LM＝-(J(h)^TJ(h)+μI)^-1J(h)^Tf(h)              （6）

其中：μ为常数，h表示变量集，代表变换矩阵h₀,h₁......h₇，J(h)为f对矩阵参数h₀,h₁......h₇偏微分的雅克比矩阵，如式7：

$J\left(h\right)=\left|\begin{array}{cccc}\frac{\∂f_1}{{\∂h}_0}& \frac{{\∂f}_1}{{\∂h}_2}& ...& \frac{{\∂f}_1}{{\∂h}_7}\\ \frac{{\∂f}_2}{{\∂h}_0}& \frac{{\∂f}_2}{{\∂h}_1}& ...& \frac{{\∂f}_2}{{\∂h}_7}\\ ...& ...& ...& ...\\ \frac{{\∂f}_n}{{\∂h}_0}& \frac{{\∂f}_n}{{\∂h}_1}& ...& \frac{{\∂f}_n}{{\∂h}_7}\end{array}\right|---\left(7\right)$

其中：表示多元函数f_n(h)对各个自变量h₀,h₁......h₇的偏微分，n从0到K，K表示迭代次数，经过K次迭代运算，获取变换矩阵H＝h_k+1；

第四步，利用第三步获得的变换矩阵H将待拼接的图像合成为高分辨率的全景图像；

根据变换矩阵H，将单帧图像按照原有的成像顺序逐一变换到同一坐标系下，然后将图像重新赋值合成在全景图像上。

2.根据权利要求1所述的一种连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法，所述的第一步包括以下几个步骤：

首先，读入待拼接图像的存储路径，待拼接图的像数目为n；

然后，根据图像的数目n和处理器的核心数k将待拼接的图像进行任务分块，即将待处理的图像平均分为k组，每组的图像数目为n/k幅图像，形成k个任务块；

最后，根据分配的任务块数目创建线程，所创建的k个线程同时并行工作，分别用来提取相应任务块中图像的SIFT特征点，其中，利用尺度不变特征提取算法来提取待拼接图像的SIFT特征点。

3.根据权利要求1所述的一种连续帧多特征点无人机侦察图像并行与匹配精度约束的拼接方法，所述的第四步中，重叠区域的图像融合过程采用基于小波变换的图像融合方法实现图像的融合，对待融合的两幅图像完成小波变换的多级分解，不断使用低通滤波器和高通滤波器进行滤波，并进行逐级降2采样，最后获得多级高频分量和低频分量；然后对两幅图像每层的系数图，采用相应的融合策略完成融合操作，最后得到融合后的小波塔形结构，对它进行小波重构变换产生融合后的图像。

Images