UAV (Unmanned Aerial Vehicle, 无人机)具有运行简单、反应迅速、飞行灵活、成本低等特点, 已被广泛应用于抗灾救灾、军事侦察、海上检测、环境保护和土地动态监测等区域^[1]. 而由于无人机航拍序列图片拍摄的视角有限, 受到无人机飞行的高度及相机的参数的限制. 为了对拍摄的区域进行全局的把握与分析, 得到更多的目标区域信息, 图像的快速拼接显得十分迫切与需要.

目前无人机的航拍影像拼接方法主要有两种: 一种是基于图像特征的拼接方法, 另一种是基于POS信息的拼接方法.

基于特征的图像拼接主要包括两个主要的步骤, 分别为图像配准和图像融合. 其中图像配准为核心步骤. 在基于特征检测的图像配准算法中, 1988年CHris Harris提出了Harris角点检测算子, 该算法对角点响应函数R进行阈值处理并且提取R的局部最大值来进行角点的提取, 缺点是该算法不具有尺度不变性并且特征检测时间不是很令人满意. SIFT特征匹配算法(Scale Invariant Feature Transform)^[2]是一种经典的算法, 它通过在空间尺度中寻找极值点, 并提取出其位置、尺度、旋转不变量, 并且提出来的特征比较稳定, 但是缺点为提取出的特征点数量较多, 计算量大, 耗时长. Bay^[3]于2006年提出了SURF算法, 意为加速鲁棒性特征, 是对SIFT算法进行的一种改进, SURF算法利用Hessian矩阵确定候选点然后进行非极大抑制, 降低了计算复杂度, 一般来说, 标准的SURF算子比SIFT算子快好几倍, 并且在多副图片下具有更好的鲁棒性.

基于POS的信息拼接方法主要是利用航拍影像的地理坐标关系来进行配准, 文献[4–7]就是利用特征点检测算法和POS信息结合的方法来进行拼接, 解决了水域拼接问题, 但是拼接速度上并没有得到有效提高, 大部分都不能达到实时拼接. 文献[8–12]都是基于视频拼接, 提高了拼接的速度, 但是融合效果并不是十分理想, 文献[13]对视频进行感兴趣区域提取, 裁剪了非重叠区域, 提高了检测的速度, 但是视频并非无人机航拍视频, 并且重叠区域的计算方法并未给出, 文献[14]利用FAST快速的提取到了图像的灰度特征, 速度上基本满足实时的要求, 但是也并未给出最后的融合效果, 文献[15]在特征定位时采用SURF, 在特征描述上采用局部差分二进制算法, 提高了特征的提取速度和匹配速度, 在特征的稳定性和去分析上也有所提高.

本文提出的方法将特征提取算法与影像的地理坐标相结合, 在拼接效率及视觉效果上比传统的特征提取拼接算法均有较大的提高, 且拼接后的影像有地理信息, 具有一定的实用价值.

1 设计方案概述

本文提出了一种基于POS系统拼接图像的方法, 首先对图像进行几何校正, 然后计算出影像四角点的地理坐标, 以第一张影像地理坐标为基础, 通过提取相邻影像重叠区域的SURF特征去得到相同匹配目标点的位置关系, 从而修正后面影像的地理坐标, 最后采用自适应渐入渐出的融合算法, 得到了一副视觉效果良好的全景图像, 完成图像的良好拼接. 设计方案流程图如图1所示.

2 航拍图像的预处理

无人机在执行航拍任务中, 由于受到飞行器的姿态、高度、速度以及地球自转等因素的影响, 造成图像相对于地面目标实际位置发生挤压、扭曲、拉伸和偏移等, 因此需要对航拍的图像进行几何矫正, 以得到基于同一基准投影面的遥感图像.

2.1 几何矫正

几何矫正处理过程主要包括坐标系统变换, 构建共线方程, 图像重采样等过程.

在几何矫正的过程中, 需要依次建立的坐标系为: 地球坐标系, 地理坐标系, 机体坐标系, 光电平台坐标系, 相机坐标系. 传统校正的具体步骤如下: 设光电平台坐标系到数码相机坐标系的变换矩阵为R₁, 由机体坐标系到光电平台坐标系的变换矩阵为R₂, 由地理坐标系到机体坐标系的变换矩阵为R₃, 由大地坐标系到空间直角坐标系的变换矩阵为R₄, 则有:

$\left\{ \begin{array}{l}R_1 = R\left( {{\beta _P}} \right)R\left( {{\alpha _P}} \right)\\R_2 = R\left( \chi \right)R\left( \omega \right)\\R_3 = {R_x}(\beta ){R_z}(\theta ){R_y}(90 - \lambda ){R_x}\left( {90} \right)\\R_4 = {R_x}^{\rm T}\left( { - \left( {90 - B} \right)} \right){R_z}^{\rm T}\left( {90 + L} \right)\end{array} \right.$

(1)

其中 ${R_x}(\lambda ), {R_{\rm{y}}}(\lambda ), {R_{\rm{z}}}(\lambda )$ 分别为绕x轴、y轴、z轴的旋转的旋转矩阵, 其中α_P为像元扫描角, β_P为像元偏角, χ为方位角, ω为平台滚动角, β为飞机横滚角, θ为航向角, B为经度, L为纬度, H为高度.

因此从相机坐标系到大地坐标系之间的转换可以用如下关系表示:

$\left[ \begin{array}{l}{\rm{Xe}}\\{\rm{Ye}}\\{\rm{Ze}}\end{array} \right]{\rm{ = }}{R_1}^{\rm T}{R_2}^{\rm T}{R_3}^{\rm T}{R_4}^{\rm T}\left[ \begin{array}{l}{\rm{Xc}}\\{\rm{Yc}}\\{\rm{Zc}}\end{array} \right]$

(2)

常规几何矫正由式(2)逐像素点进行计算和采样. 考虑到本课题主要对单一航带进行拼接, 且无人机平台横滚角和俯仰角比较小时(近似垂直向下拍摄时)可以认为图像分辨率率是一致的, 因此为了减少运算量, 提高速度, 几何矫正只对水平方向上的飞机航向角和平台方位角进行校正. 矫正模型如下所示:

$\left[ \begin{array}{l}{\rm{Xe}}\\{\rm{Ye}}\\{\rm{Ze}}\end{array} \right]{\rm{ = }}{R^{\rm T}}\left( {\chi {\rm{ + }}\theta } \right)\left[ \begin{array}{l}{\rm{Xc}}\\{\rm{Yc}}\\{\rm{Zc}}\end{array} \right]$

(3)

其中χ为方位角, θ为航向角, ${R^{\rm T}}\left( {\chi {\rm{ + }}\theta } \right)$ 为变换矩阵, 完成模型转换后再通过双线性差值重采样方法计算出像元在新的坐标系下的灰度值, 从而生成新的图像矩阵, 完成几何矫正.

2.2 重叠区域计算

无人机的航行速度, 航行高度, 当前曝光点的飞机经纬度信息, 都可以通过无人机携带的POS系统(图2)获取到. 这些信息记录了飞机飞行时在每个曝光时刻飞行时的姿态, 可根据此来计算相邻影像的重叠区域, 具体步骤如下.

1) 假设无人机航线为偏离正北θ, 其航速为V, 将速度分解为正北方向上和正东方向上两个方向上的分量. 分别表示为V₁和V₂, 如图3所示.

2) 连续拍摄相邻两张曝光点的影像分别编号为Pic1和Pic2, 间隔为t, 中心点经纬度为LatA、LonA、LatB、LonB, 则有:

$\left\{\begin{array}{*{20}{l}}\begin{aligned}C & = {\rm{sin}}( {LatA} )*{\rm{sin}}( {LatB} )*{\rm{cos}}( {LonA - LonB}) \\ &+ {\rm{cos}}( {LatA} )*{\rm{cos}}( {LatB} )\end{aligned}\\{{L}} = R*{\rm{Arccos}}( C )* {\text{π}} /180\end{array}\right.$

(4)

其中R=6371.004千米, π取3.1415926.

3) 此时以Pic1作为基准, Pic2与Pic1之间的重叠区域可表示为如图4所示的矩形区域S, 将其重叠区域扩大定义到规则区域, 将重叠区域顶点投影到x和y方向上则表示为:

$S = (W - P_1)(H - P_2)$

(5)

通过计算重叠区域减少了拼接的计算量, 并且随着特征点的寻找面积越小, 可能出现的误匹配概率就越小, 提高了检测的效率.

3 基于POS参数的拼接方法 3.1 计算影像地理坐标

对于几何纠正过的航拍影像, 影像的每点像素坐标都有对应的地理坐标, 影像的地理坐标的计算过程可以分为以下几步.

1) 计算地面分辨率, 计算公式如下:

${GSD = }\frac{{{H*P}}}{{f}}$

(6)

其中GSD表示地面分辨率(m), f为镜头焦距(mm), P为成像传感器的像元尺寸(mm), H为无人机对应的航高(m).

2) 计算影像对角线实际地面距离, 根据影像的宽度w和高度h得到图像对角线线之间的地面距离:

${L = }\sqrt {{{w}^2}{\rm{ + }}{{h}^2}} {*GSD}$

(7)

3) 计算影像四角点地理坐标, 根据影像中心点经纬度及另一点相对中心点的距离及方向角, 可以求得以影像中点为圆心, 半径为L/2的对应四角点地理坐标. 具体计算公式为:

${{Lo}}{{{n}}_i} = \frac{{\left( {\frac{{L\sin ({\theta _i})}}{{2({R_j} + ({R_c} - {R_j})*(90 - La{t_a})/90)*\cos (La{t_a})}} + Lo{n_a}} \right)*180}}{\text{π}}$

(8)

${{La}}{{{t}}_i} = \frac{{\left( {\frac{{L\cos ({\theta _i})}}{{2({R_j} + ({R_c} - {R_j})*(90 - La{t_a})/90)}} + La{t_a}} \right)*180}}{\text{π}}$

(9)

其中 ${\theta _i} \in (0,2*{\text{π}})$ , Lon_a, Lat_a为影像中心点的经纬度, R_i为赤道半径取6378 137 m, R_j为极半径取 6356 725 m.

4) 地理坐标转换到空间之间坐标系之间的转换, 转换公式如下:

$\left[ \begin{align}{X_s}\\{Y_s}\\{Z_s}\end{align} \right] =\left[ \begin{aligned}&(N + H){\rm{cos}}(Lon){\rm{cos}}(Lat)\\&(N + H){\rm{cos}}(Lon)\sin (Lat)\\&[N({1} {-} {e^2}) + H]\sin(Lon)\end{aligned} \right]$

(10)

其中N为曲率半径, Lon、Lat、H分别为经度、纬度和高程, 至此, 得到完整的影像坐标.

3.2 利用特征匹配修正地理坐标

由于POS精度低以及几何校正中存在一定的误差, 计算出的坐标映射关系有一定的误差, 此时利用特征匹配算法去修正影像的地理坐标, 具体过程如下:

假设图像1的地理坐标为P1(x1, y1), 图像2的地理坐标为P2(x2, y2), 提取影像特征匹配对后则可以得到同一个点在不同两张影像间的像素坐标位置, 由此可以得到两张影像中同一个目标点的经纬度坐标(Lon1, Lat1)(Lon2, Lat2), 最后以第一张影像地理坐标为基准, 求得第二张影像中目标点与第一张影像之间的偏移量, 公式如下:

$\left\{ \begin{array}{l}x = Lon2 - Lon1\\y = Lat2 - Lat1\end{array} \right.$

(11)

然后用求得的偏移量去修正图像2的地理坐标P2:

$\left\{ \begin{array}{l}x2 = x2 + x\\y2 = y2 + y\end{array} \right.$

(12)

然后将地理坐标投影到空间直角坐标系下, 从而完成图像精确的配准.

文章通过计算重叠区域的SURF特征, 完成特征匹配点对的提取. SURF采用近似Hessian矩阵检测特征点, 并使用积分图像进行卷积运算, 大幅减少了运算, 从而提高了特征提取速度. SURF描述子包含两个主要的部分: 检测特征点和计算特征; 具体实现分为以下几步.

1) 构建Hessian矩阵, 构造尺度空间.

假设图像上某点为X(x, y), 在σ尺度下的矩阵H定义为:

${\rm{H}}\left( {x,\sigma } \right){\rm{ = }}\left[ \begin{array}{l}{L_{xx}}\left( {x,\sigma } \right){L_{x{\rm{y}}}}\left( {x,\sigma } \right)\\{L_{x{\rm{y}}}}\left( {x,\sigma } \right){L_{{\rm{yy}}}}\left( {x,\sigma } \right)\end{array} \right]$

(13)

当Hessian矩阵的的判别式取得局部极大值时, 认为定位关键点的位置.

2) 检测特征点

在上述得到的尺度空间中, 将经过Hessian矩阵处理的每个像素点与二维图像空间和尺度空间邻域内的26个点进行比较, 初步定位出关键点, 再经过滤去除能量比较弱的关键点以及错误定位的关键点, 筛选出最终的稳定的特征点, 检测过程如图5.

3) 确定特征点主方向

以特征点为圆心, 以6σ为半径, 求XY方向上的Haar小波响应, 统计60度扇形内所有点的水平Haar小波特征和垂直Haar小波特征总和, 并设置Haar小波的尺寸变长为4 s, 使得每个扇形都得到了对应值. 然后将60度扇形以一定间隔进行旋转, 最后将锁定最大值扇形的方向作为该特征点的主方向. 该过程的示意图如图6.

4) 计算特征描述符

在特征点周围取一个正方形框, 然后把该框分为16个子区域, 每个子区域统计25个像素的水平方向和垂直方向的Haar小波特征, 这里的水平和垂直方向都是相对主方向而言的, 所以每个特征点就是16*4=64维的向量, 这在特征匹配过程中会大大加快匹配速度.

4 图像融合

对于修正坐标后的图像, 由于每张影像之间有一定的痕迹, 因此需要采用一些策略来解决拼接缝隙之间的颜色过渡差异较大的问题, 使拼接图像更加平滑自然, 渐入渐出自适应权重融合过程如下所示.

假设I₁、I₂、I分别为融合前图像1、图像2、和融合后的图像3, 则可以通过公式(14), 完成图像融合:

$\begin{array}{l}I\left( {x,y} \right) = \left\{ \begin{array}{l}{I_1}\left( {x,y} \right)\\\varphi {I_1}\left( {x,y} \right){\rm{ + }}\left( {1{\rm{ + }}\varphi } \right){I_2}\left( {x,y} \right)\\{I_2}\left( {x,y} \right)\end{array} \right.\\\varphi = {w_i}/W;\varphi \in \left( {0,1} \right)\end{array}$

(14)

式中, W是两幅不同画面重复的总宽度; w是重合区域左边缘与当前像素点之间的横向距离.

5 实验与分析

本文选取我国某处海岸线的两条航带影像作为数据源, 测试平台为CPU Intel(R)Core(TM)i7-4970 CPU@3.60 GHz, 内存16 GB, 软件开发环境为VS2010+OpenCV2.4.1, 图片分辨率为1820*1080, 航拍高度为150 m左右.

5.1 重叠区域特征检测

图7(a)为用传统SURF算法检测到的特征点, 并进行匹配连线, 图7(b)为根据飞行参数计算得到的重叠区域, 图7(c)为通过影像重叠区域计算出特征点之后, 再根据RANSAC算法进一步消除错误特征点, 从而得到最终的匹配连线对.

可以看出, surf算子在图7(a)中非重叠区域之外的特征也进行了检测, 所以存在很多的误匹配对. 图7(c)只在重叠区域进行特征检测, 减少了误匹配率. 从而提高了拼接的效率. 本文通过测试561张图片并对得到的数据进行了平均, 实验中全局区域和重叠区域检测到的特征点个数及检测时间对比如表1所示.

5.2 地理坐标拼接效果

按照本文提出的计算方法计算地理坐标并进行投影后及修正后的对比结果如图8所示: 其中图8(a)为未修正前的地理坐标投影结果, 可以看出修正前的拼接效果存在明显的误差, 图8(b)在图8(a)的基础上, 按照本文提出的方法, 进行坐标的调整后, 图像拼接精度得到了明显的提高, 图8(c)为按照本文方法, 进行多张图片的修正地理坐标拼接效果.

5.3 图像融合结果

根据每张影像的地理坐标可以方便提取出图像的边界, 从而得到图像之间的拼接缝隙的精确位置, 这样重叠区域的目标在融合过程中都是以地理坐标位置进行投影, 消除了传统图像处理算法配准时候误差的积累, 从而避免目标重影的产生, 最后采用渐入渐出的自适应权重融合算法对边界线进行融合, 对比拼接效果如图9所示: 其中图9(a)明显可以看出拼接后的边界, 图9(b)为采用融合策略后生成的图像, 可以看出消除了图像边界之间的差异.

6 总结

本文针对无人机航拍影像的特点, 提出了一张基于无人机POS信息的拼接方法, 首先根据POS信息计算影像的四角点的地理坐标, 然后提取影像的重叠区域的SURF特征, 利用相邻影像的相同特征点的位置关系去修正地理坐标, 从而完成配准, 最后采用自适应渐入渐出融合算法对图像进行平滑过渡, 得到了一副视觉效果良好的全景图像.