近年来，随着无人机航拍技术的普及，航拍图像在目标检测^[1-3]、目标跟踪^[4-5]等各个领域得到广泛应用。无人机航拍图像与卫星影像进行匹配可以对无人机进行图像辅助定位或者实现航拍图像与卫星影像融合制作高清地图等，具有一定的应用价值。但是由于成像机制的不同，实现无人机航拍图像与卫星影像的匹配存在着几个难点。

一方面是这2种图像之间的尺度差异极大。卫星影像是一种大尺度、宽视野图像，其单个像素点所覆盖的实际地面区域比较大。相对来说，无人机航拍图像是一种小尺度、窄视野的图像，由于无人机飞行高度限制，航拍图像所覆盖的实际地面区域有限，其单个像素点所覆盖的区域较小。另一方面是2种图像之间的特征差异较大。由于成像机理、成像时间、成像角度以及光照等多方面的差异，航拍图像与卫星影像在地貌、颜色、亮度等多方面存在较大的特征差异，因而实现航拍图像与卫星影像的匹配具有较大的挑战。

目前，图像匹配技术主要可以分为基于特征的图像匹配和基于模板的匹配两类。基于特征的图像匹配技术有SIFT^[6-7]特征、SURF^[8-9]特征以及ORB^[10-11]特征等，它们通过图像的局部特征来实现匹配，对图像的小尺度差异具有一定的鲁棒性。但是无法实现航拍图像与卫星影像这种具有大尺度差异图像的有效匹配。基于模板的匹配是目前应用比较广泛的图像匹配技术，文献[12]提出一种BBS(best-buddies-similarity)匹配算法，该方法通过将模板图像和搜索图像分为若干个小图像块，然后统计模板图像和参考图像中互为最近邻匹配的图像块数量作为相似度衡量标准。这种分块统计的匹配方式在一定程度上减小了图像遮挡对匹配结果的影响，但是对于存在微小扭曲形变的图像匹配没有适应能力。文献[13]提出一种DDIS(deformable diversity similarity)匹配方法，该方法充分考虑了模板图像与目标匹配区域之间的变形问题，对于存在复杂变形和遮挡变化等干扰具有较强的适应性，但是在非均匀光照条件下的匹配鲁棒性不高。

现有的图像匹配方法在常用的匹配场景下具有较高的匹配精度，但是无法对尺度和特征都存在较大差异的航拍图像与卫星影像实现有效匹配。因此，本文提出了一种适应于具有大尺度差异的航拍图像与卫星影像匹配方法。该方法将航拍图像成像时刻的相机位姿参数与图像匹配结合起来，充分利用相机位姿参数对图像进行配准以减小尺度差异对匹配的影响; 再基于CNN特征对配准后的图像进行匹配，从而能够有效地实现航拍图像与卫星影像的匹配。

1 航拍图像与卫星影像匹配概述

本文的目标是实现无人机航拍图像与卫星影像的匹配。由于航拍图像与卫星影像存在大尺度差异和旋转关系难以实现匹配，本文以卫星影像为基准，通过充分利用航拍图像相机成像时刻的位姿参数等信息来对航拍图像进行旋转和放缩操作，实现方向和尺度的配准，再将配准的航拍图像与卫星影像进行匹配。在匹配过程中，为了提高匹配效率，分两步对图像进行匹配。首先，利用航拍图像的相机位姿参数在卫星影像中找到航拍图像可能匹配的区域截取卫星子图得到一个粗匹配结果; 然后，将卫星子图与配准的航拍图像输入到一个预训的CNN特征提取网络提取图像特征向量; 最后基于CNN特征对航拍图像与卫星影像进行精匹配，具体匹配流程见图 1。

在图 1的匹配方法流程中，航拍图像的配准以卫星影像为基准。首先通过相机的方位角确定航拍图像与卫星影像的旋转关系，并通过图像旋转实现两者方向校正。然后根据图像实际地面覆盖区域和像素尺寸关系分别计算出卫星影像和航拍图像单个像素点的地面覆盖区域大小，并通过图像重采样对航拍图像进行放缩，使得航拍图像与卫星影像单个像素点的地面覆盖区域大小一致，从而实现尺度空间的配准。在图像的粗匹配阶段，根据航拍图像成像时刻的相机位置参数和卫星影像的经纬坐标信息在卫星影像上截取包含航拍图像在内的卫星影像子图作为粗匹配结果。在精匹配阶段，首先使用特征提取网络提取到航拍图像特征向量与卫星影像粗匹配结果的卫星子图特征向量；然后做余弦相乘并根据本文计算方法计算它们的特征响应得分；最后寻找特征响应得分最大的区域作为匹配区域，从而实现航拍图像与卫星影像的精匹配。

在本文中，卫星影像为俯拍成像，以图像水平方向为经度方向，垂直方向为纬度方向。航拍图像为无人机机载相机俯拍成像，以正北方向为航拍相机方位角0°方向，取顺时针方向为正值，逆时针方向为负值，航拍图像以卫星影像为基准的示意见图 2。

图 2中实线框表示卫星影像，虚线框表示航拍图像。根据示意图可见，无人机航拍过程中，航拍相机方位角A_z的改变会导致航拍图像与卫星影像之间存在旋转变换关系，旋转关系的存在会影响航拍图像与卫星影像的匹配精度，因此需要对航拍图像进行旋转方向的配准。

2 航拍图像与卫星影像匹配实现方法 2.1 基于相机位姿参数的航拍图像配准

航拍图像以卫星影像为基准，首先进行方向配准，然后进行尺度配准。本文通过解码图像成像时航拍相机的方位角A_z，然后将图像按方位角相反方向进行旋转实现图像的方向配准。图像旋转的本质是像素点坐标位置绕图像中心的旋转变换，这一变化过程可以通过变换矩阵来实现。但是由于图像的像素坐标以图像左上角为原点，因此图像变换前首先需要将图像由像素坐标系变换到图像中心坐标系，具体变换关系见图 3。

根据图 3所示的图像像素坐标系和中心坐标系的关系可得，当假设图像高度为H_u，宽度为W_u时，像素坐标系的一点(X₀, Y₀)变换到中心坐标系下的坐标为(X, Y)，其具体变换见式(1)。

$\left[\begin{array}{l} X \\ Y \\ 1 \end{array}\right]=\left[\begin{array}{ccc} 1 & 0 & -0.5 W_u \\ 0 & -1 & 0.5 H_u \\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{l} X_0 \\ Y_0 \\ 1 \end{array}\right] $

(1)

由于图像的显示是基于像素坐标系的，因此在图像中心坐标系下实现图像旋转后需要将中心坐标系下的坐标进行逆变换到像素坐标系下，由图像中心坐标系变换到像素坐标系的变换见式(2)。

$\left[\begin{array}{l} X_0 \\ Y_0 \\ 1 \end{array}\right]=\left[\begin{array}{rrc} 1 & 0 & 0.5 W_u \\ 0 & -1 & 0.5 H_u \\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{l} X \\ Y \\ 1 \end{array}\right] $

(2)

图像在中心坐标系下实现旋转如图 3(b)所示，假设图中点(x₀, y₀)按顺时针方向旋转θ度后的坐标为(x, y)，根据几何关系可以写出坐标关系如下:

$\left[\begin{array}{l} x \\ y \\ 1 \end{array}\right]=\left[\begin{array}{ccc} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{l} x_0 \\ y_0 \\ 1 \end{array}\right] $

(3)

综合上述推导，图像像素坐标系下的一点(X₀, Y₀)绕图像中心顺时针旋转θ度后的坐标(X′₀, Y′₀)见式(4)。

$ \begin{aligned} {\left[\begin{array}{l} X_0^{\prime} \\ Y_0^{\prime} \\ 1 \end{array}\right]=} & {\left[\begin{array}{rrc} 1 & 0 & 0.5 W_u \\ 0 & -1 & 0.5 H_u \\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{ccc} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array}\right] \cdot } \\ & {\left[\begin{array}{rrr} 1 & 0 & -0.5 W_u \\ 0 & -1 & 0.5 H_u \\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{l} X_0 \\ Y_0 \\ 1 \end{array}\right] } \end{aligned} $

(4)

由于图像的像素坐标为正整数，经过旋转后的图像对像素坐标进行取整时会导致部分像素点没有灰度值而出现像素离散的现象。为了克服这一缺点，本文采用反向映射和双线性插值^[14]的方式对变换后的图像进行处理。反向映射从旋转后的图像出发，计算旋转后的像素坐标(x, y)在原图像中的坐标位置(x′, y′)，然后对其周围4个像素点的灰度值进行插值运算得到灰度值，并将该灰度值赋值给旋转后的图像，具体过程见图 4。

在图 4(b)所示的双线性插值中，点P的灰度值表示为f(P)，首先对点P周围的4个像素点在x轴方向做2次线性插值得到在点R₁和R₂上的值，然后在y轴方向上做线性插值得到点P的灰度值，最终插值计算见式(5)。

$\begin{aligned} f(P)= & \frac{\left(y_2-y^{\prime}\right)\left(x_2-x^{\prime}\right)}{\left(y_2-y_1\right)\left(x_2-x_1\right)} f\left(Q_{11}\right)+\frac{\left(y_2-y^{\prime}\right)\left(x^{\prime}-x_1\right)}{\left(y_2-y_1\right)\left(x_2-x_1\right)} f\left(Q_{21}\right)+ \\ & \frac{\left(y^{\prime}-y_1\right)\left(x_2-x^{\prime}\right)}{\left(y_2-y_1\right)\left(x_2-x_1\right)} f\left(Q_{12}\right)+\frac{\left(y^{\prime}-y_1\right)\left(x^{\prime}-x_1\right)}{\left(y_2-y_1\right)\left(x_2-x_1\right)} f\left(Q_{22}\right) \end{aligned} $

(5)

通过式(5)实现图像方向配准后，再对图像尺度进行配准，本文通过计算图像的实际地面覆盖距离与图像像素尺寸大小来计算图像尺度。首先计算卫星影像的实际地面覆盖距离，在图 5的卫星影像中，影像左上角A的经纬坐标为(l_{A, lon}, l_{A, lat}), 右下角B的经纬坐标为(l_{B, lon}, l_{B, lat})。由于所取的卫星影像地面覆盖区域大小相对于地球来说非常的小，因此可以忽略地球球面的影响，参考文献[15]的计算方式，在水平方向和竖直方向的地面覆盖区域L_ax和L_ay可以根据式(6)、(7)求解。

$ L_{a x}=\left|l_{\mathrm{B}, \mathrm{lon}}-l_{\mathrm{A}, \mathrm{lon}}\right| \cdot \frac{\cos \left(l_{\text {lat }}\right) \times\left(R-l_{\text {lat }} \times(R-r) / 90\right) \times \pi}{180} $

(6)

$ L_{a y}=\left|l_{\mathrm{B}, \text { lat }}-l_{\mathrm{A}, \text { lat }}\right| \cdot \frac{\left(R-l_{\text {lat }} \times(R-r) / 90\right) \times \pi}{180} $

(7)

式中：l_lat为A、B两点纬度的均值，R=6 378 137 m为地球赤道半径大小，r=6 356 752.314 m为地球极轴半径大小。

如图 6所示，无人机航拍图像实际地面覆盖距离可以根据式(8)、(9)求解。

$ L_x=2 L \cdot \tan \frac{\theta_x}{2} $

(8)

$ \begin{aligned} L_y=2 L \cdot \tan \frac{\theta_y}{2} \end{aligned} $

(9)

式中：L为相机到地面的距离，θ_x×θ_y为相机视场角大小。

假设卫星影像图像像素大小为H_a×W_a，H_a为图像行像素点数，W_a为列像素点数。相应地，无人机航拍图像像素大小为h×w。对于卫星影像可以计算出水平和竖直方向上单个像素点实际占据地面范围大小，见式(10)、(11):

$ \begin{aligned} \mathrm{d} U=\frac{L_{a x}}{W_a} \end{aligned} $

(10)

$\mathrm{d} V=\frac{L_{a y}}{H_a} $

(11)

同理，根据无人机航拍图像地面覆盖区域长L_x和宽L_y以及图像的大小h×w可以求得航拍图像在水平和竖直方向上单位像素大小分别为du、dv。因此，航拍图像与卫星影像在水平和竖直方向上的尺度关系式(12)、(13)如下：

$ \sigma_x=\frac{\mathrm{d} u}{\mathrm{~d} U} $

(12)

$ \begin{aligned} \sigma_y=\frac{\mathrm{d} v}{\mathrm{~d} V} \end{aligned} $

(13)

得到航拍图像与卫星影像的尺度关系后，通过对航拍图像进行重采样实现尺度配准。假设无人机航拍图像重采样调整后的图像大小为h′×w′，其数值式(14)、(15)如下:

$ \begin{aligned} w^{\prime} =\left\lfloor\sigma_x \cdot w\right\rfloor \end{aligned} $

(14)

$ h^{\prime} =\left\lfloor\sigma_y \cdot h\right\rfloor $

(15)

式中$\left\lfloor \cdot \right\rfloor $为对数据向下取整。

2.2 基于分步的航拍图像与卫星影像匹配方法

在图像匹配阶段，本文采用粗匹配与精匹配结合的分步匹配方式。粗匹配以无人机航拍图像成像时刻记录的相机经纬位置信息(l_{T, lon}, l_{T, lat})为航拍图像中心点的经纬位置，然后根据相机定位精度设定一个阈值R_t。在卫星影像中以航拍图像经纬度为中心，阈值R_t为半径在全局卫星影像中截取卫星子图作为与无人机航拍图像粗匹配的结果。卫星影像子图的选取示意见图 7，其中虚线框区域为根据地面覆盖区域距离计算以及阈值R_t大小选取的卫星影像子图区域。

经过粗匹配缩小匹配范围后的卫星影像子图与配准后的航拍图像进行精匹配。在精匹配阶段，本文以配准后的航拍图像为模板图像T，卫星影像子图定义为搜索图像S，匹配的目标是在搜索图像中找到模板图像的最佳匹配区域。根据参考文献[16-17] 定义一个特征匹配质量得分函数Q(s, t)，其中t和s分别是模板图像T和搜索图像S中的一小块区域，它们对应的特征向量表示为f_t和f_s。

给定一个搜索图像S中的块区域s，然后定义模板图像T中的块区域t与之相匹配的概率为

$ P(t \mid s)=\frac{\exp \left\{\alpha \cdot \rho\left(f_t, f_s\right)\right\}}{\sum _{t^{\prime} \in T} \exp \left\{\alpha \cdot \rho\left(f_{t^{\prime}}, f_s\right)\right\}} $

(16)

式中：α为正权重，ρ(·)为特征向量余弦相乘运算。P(t|s)为当前块区域t相对于模板图像T中所有其他块区域与搜索图像S中块区域s的匹配程度的排序，其值越大，块区域t和块区域s匹配的可能性越高。

根据P(t|s)的定义，同理可以定义搜索图像S的一个块区域s与模板图像T给定块区域t的匹配概率为P(s|t)，因此特征匹配质量得分函数Q(s, t)如下：

$Q(s, t)=P(t \mid s) \cdot P(s \mid t) $

(17)

得到特征匹配质量得分函数Q(s, t)后，根据整体匹配质量最优化原则，在匹配质量得分函数中找出匹配质量得分响应极大值区域即为最佳匹配区域，寻找最佳匹配区域为

$ R=\underset{S}{\operatorname{argmax}}\left\{\sum\limits_{s \in S} \max \{Q(s, t) \mid t \in T\}\right\} $

(18)

式中R为在搜索图像S中找到的与模板图像T的最佳匹配区域。

3 匹配方法性能分析 3.1 算法复杂度分析

假设配准后的航拍图像大小为$T_w \times T_h \triangleq l $，卫星影像大小为$ S_w \times S_h \triangleq L$，经过本文粗匹配方法后得到的卫星影像子图大小为$ S_{w^{\prime}} \times S_h{ }^{\prime} \triangleq L^{\prime}$，提取的特征维度大小为d。根据匹配过程，首先进行图像特征向量余弦相乘计算，计算结果用一个L′×l维的矩阵D来存储，D_ij表示块区域s_i和t_j特征向量余弦相乘的结果，则其计算时间复杂度为O(L′ld)，内存占用复杂度为O(L′l)。然后计算特征匹配质量得分函数Q(s, t)，根据Q(s, t)的定义同样可以用一个L′×l维的矩阵Q来存储Q(s, t)，其计算时间复杂度为O(L′l)，内存占用复杂度为O(L′l)。最后在Q(s, t)中寻找匹配质量得分最高的匹配区域，其寻找过程为首先在Q的每一行中寻找出最大值得到一个L′×1的向量，然后将向量维度调整为S_w′×S_h′维的得分矩阵，最后以T_w×T_h为滑动窗口在得分矩阵中寻找得分最高的区域作为最佳匹配区域输出，其计算的时间复杂度为O(L′l)，内存占用复杂度为O(L′)。

将上述复杂度进行综合并舍弃数量级较低的变量, 即可得到本文所提匹配算法总的时间复杂度为O(L′ld)，内存占用复杂度为O(L′l)。根据文献[12]，BBS算法的时间复杂度为O(Ll²d)，内存占用复杂度为$O\left(l^2 \sqrt{L}\right) $。根据文献[13]，DDIS算法的时间复杂度为O(Ll+d′llog l)，内存占用复杂度为O(L)，其中d′为特征向量d经过主成分分析算法降维后的特征维度大小。由于经过粗匹配后L′ $ \ll $L，因此相比于BBS算法和DDIS算法，本文所提匹配算法在时间复杂度和内存占用复杂度上都具有较大的优势。

3.2 实验条件设置

航拍图像采集设备为大疆MATRICE 200 V2无人机，搭载大疆禅思X5S云台相机。图像采集地点为哈尔滨工业大学田径场，无人机航拍飞行高度为100 m，航拍图像像素大小为5 280×2 970，对应的实际地面覆盖范围为115.33 m×64.67 m。相应的卫星影像从BIGEMAP软件上下载，其像素大小为1 432×879，对应的地面覆盖经纬坐标范围大小分别为: 经度126.624 550°E~ 126.632 227°E，纬度45.734 169°N~45.738 873°N，具体的航拍图像和卫星影像实例见图 8。

在本文实验中，粗匹配阶段选取的阈值R_t大小为150 m。在精匹配阶段利用在ImageNet上预训练的VGG19网络来提取图像的CNN特征。为了客观地评价本文所提算法的优势，本文将匹配结果与现有的优秀模板匹配算法BBS和DDIS进行比较，并利用匹配结果与真实匹配位置的交并比I_IoU(intersection over union，IoU)来定性评价图像匹配质量的好坏，I_IoU的具体定义公式如下：

$ I_{\mathrm{IoU}}=\frac{S_{A B}}{S_A+S_B-S_{A B}} $

(19)

式中：S_A为真实匹配位置面积，S_B为匹配结果位置面积，S_AB为匹配结果和真实匹配的交叠部分面积。

3.3 结果分析

本次实验共采集了8张无人机航拍图像与卫星影像进行匹配分析，具体匹配结果见图 9。

在图 9中，绿色框为航拍图像与卫星影像匹配的真值框，蓝色框为BBS算法匹配结果，黄色框为DDIS算法匹配结果，红色框为本文所提算法匹配结果。在航拍图像与卫星影像方向和尺度严格配准的情况下，红色匹配框与绿色真值框的形状应该是一种面积等大，边界平行的关系。通过对比匹配结果的红色匹配框和绿色真值框之间的方向和大小关系可见，利用相机位姿参数进行方向和尺度配准的航拍图像与卫星影像之间是存在一定畸变的。这是因为无人机在航拍过程中存在抖动以及风力干扰等因素导致相机位姿参数存在一定的偏差，使得配准结果存在微小畸变，影响图像匹配的精度，从而导致匹配结果红色框和绿色真值框之间存在一定的匹配误差。从匹配结果来看，在存在微小畸变的情况下，本文利用CNN特征进行匹配的结果仍然能够定位到图像的匹配位置，说明所提的匹配算法具有一定的图像畸变鲁棒性。

通过图 9所示3种图像匹配方法比较结果可见，本文所提的匹配方法能够很好的匹配到航拍图像的真实匹配区域，而BBS算法和DDIS算法仅有少数情况下能够实现匹配，其误匹配现象比较严重。3种匹配方法的匹配结果I_IoU评估计算见表 1, 情况1~8为图 9中的8种不同情况。

由表 1可见，虽然在个别情况下BBS和DDIS匹配结果的I_IoU高于本文所提算法，但是本文所提算法匹配结果的I_IoU都在0.6以上，而BBS和DDIS均出现I_IoU为0的不可匹配的情况，并且本文所提算法匹配结果的I_IoU均值远高于BBS算法和DDIS算法。这说明本文算法在匹配大尺度差异的航拍图像与卫星影像上具有较高的鲁棒性，能够很好的适应图像尺度差异变化对图像匹配结果的影响。由此证明了本文所提的图像匹配方法在匹配具有大尺度差异的航拍图像与卫星影像场景下具有较大的优越性。

4 结论

提出了一种实现具有大尺度差异的航拍图像与卫星影像匹配方法。该方法以卫星影像为基准，根据航拍相机的位姿参数对航拍图像进行方向和尺度空间的配准，然后根据航拍位置信息在卫星影像中截取出卫星子图与配准后的航拍图像进行匹配。实验仿真分析结果表明，现有图像匹配算法在匹配尺度差异较大的航拍图像与卫星影像时存在比较严重的误匹配现象，它们匹配结果I_IoU均值都在0.5以下，而本文所提算法匹配结果I_IoU均值为0.72，远高于现有的图像匹配算法。说明本文所提算法具有较高的鲁棒性，能够较好的实现尺度差异较大的航拍图像与卫星影像的匹配。这对于航拍图像与卫星影像融合、无人机基于卫星影像的辅助视觉定位以及目标跟踪等应用领域具有较大的应用价值。