0引言遥感图像配准技术是从不同时间、不同拍摄角度通过不同的传感器捕获多组复杂的遥感图像进行叠加和匹配的研究过程。该项技术在评估气候变化、环境污染监测和城市布局变化研究等方面有十分广泛的应用。为了提高配准细节的准确性和应用类别的多样性，高精度的配准方法和高效率的配准模型是当前遥感图像配准中研究的重点目标。遥感图像的配准方法主要有3种，分别为传统经典算法、半传统半深度学习相结合的算法以及基于端到端深度学习的算法。1) 传统经典算法。主要应用在遥感图像配准的早期阶段，采用基于灰度、变换域和特征的3种传统手工设计的方法进行配准。(1)基于灰度的方法是通过两幅图像之间灰度的相似性进行计算，对灰度变化比较敏感，但是无法对两幅灰度值变化较大的图像进行处理(Ma等，2021)。(2)基于变换域的方法是利用平移、旋转或缩放等变换方式进行频域内的配准，但是容易受到视角变化的干扰，配准精度较低(梁勇等，2010)。(3)基于特征的方法是利用手工设计的特征点进行两幅图像的匹配，计算量明显少于基于灰度的方法和基于变换域的方法。如Yu和Morel(2011)提出的ASIFT(affine-SIFT)方法、Rublee等人(2011)提出的ORB(oriented fast and rotated brief)和刘瑞红(2015)提出的SIFT-Like算法。并且手工设计的局部特征点的提取方法具有一定的尺度不变性、鲁棒性和抗干扰能力，能够克服一定的视角变化和灰度变化的干扰。随着图像配准领域的不断发展，传统配准算法逐渐出现了手工成本过大、算法运行时间过长和缺乏自适应的特征学习等缺点(Wang等，2018)。在遥感图像配准领域中越来越多的学者开始尝试使用深度学习的方法来代替传统方法进行特征提取(许东丽和胡忠正，2019)。2) 半传统半深度学习算法。Quan等人(2016)利用深度神经网络提取图像特征并学习匹配关系，缩短了特征提取网络运行的时间，然后使用改进的RANSAC(random sample consensus)算法进行错误匹配点的剔除；Yang等人(2018)提出利用预训练的VGG(Visual Geometry Group)网络进行特征提取，结合改进的CPD(coherent point drift)算法进行遥感图像配准；叶发茂等人(2019)研究表明经过预训练微调后的卷积神经网络(convolutional neural networks，CNN)得到的特征比传统方法得到的特征要具有更强的鲁棒性，并且提高了遥感图像配准整体的性能。以上方法相比传统方法(Wang等，2015)可以自适应地学习特征，但是对复杂的图像来说配准效果不好，运行时间仍然过长。为了提高模型的运行效率和配准性能，学者们开始尝试端到端的方法，代替传统方法构造一体化的网络模型。3) 端到端深度学习算法。Rocco等人(2017)设计了一种端到端的配准框架，对特征提取、特征匹配和参数回归3个部分整体训练。测试时调用训练好的模型，运行时间比半传统半深度学习的方法提高了很多，但是缺乏考虑来自不同背景的噪音对匹配关系的影响。之后Seo等人(2018)使用偏移感知机制处理图片的噪音，并对匹配部分进行改进，得到了比Rocco方法更好的匹配关系。为了进一步提高配准的精度，Kim等人(2019b)使用预训练的残差网络(He等，2016)进行特征提取，同时考虑到图像受时间和天气等变化因素的影响，利用皮尔逊相关性改进了从源图像到目标图像方向上的匹配关系，实验结果表明此方法比Rocco等人(2017)和Seo等人(2018)提出的方法效果要好。但是Kim等人(2019b)使用预训练网络提取特征时容易忽视合适的目标区域，造成无关点提取过多，重要目标点提取过少。同时在特征匹配过程中忽视了方向的一致性原理(张凤晶等，2016)，仅仅得到从源图像到目标图像单方向的配准关系，回归出单方向的变换参数指导配准，容易导致配准结果不理想。针对上述配准方法存在的问题，本文通过关注重要特征(Veličković等，2017)、研究匹配结构(王丽芳等，2020)，提出基于参数合成空间变换网络的双向一致性遥感图像配准方法。本文的主要贡献有：1) 参考逆向合成空间变换网络(Lin和Lucey，2017)的原理和李红艳等人(2019)使用注意力机制改进骨干网络的方法，对具有空间注意力机制的空间变换网络(Jaderberg等，2015)进行改进，提出参数合成的空间变换网络进行特征提取。2) 参考郑莹和李光耀(2011)提出的信息结合双向配准方法和Kim等人(2019a)提出的循环一致性原理，增加一条匹配和参数回归分支，进行双向一致性的匹配，增强了模型的鲁棒性。3) 参考Ji等人(2015)提出的双向匹配对加权的思想，将匹配后得到的两个参数加权合成，根据不同方向的匹配关系对配准结果的影响，设置不同的权重，提高了配准的精度。1相关工作1.1空间变换网络传统的卷积神经网络由于池化层本身的局限性，缺乏空间不变性，从而影响网络特征提取的准确性。为了提高网络的性能，Jaderberg等人(2015)设计了空间变换模块，并将此模块插入卷积神经网络中，提出了具有空间不变性的空间变换网络。该网络利用空间变换模块聚焦图像中感兴趣的部分，为得到下一层的期望形式进行了铺垫。空间变换网络由定位网络、网格生成器和采样器3部分组成。第1个组成部分是定位网络。输入的特征图通过由卷积层、池化层和全连接层组成的定位网络预测出所需的变换参数$\mathit{\pmb{θ}}$。如果采用仿射变换的方式处理图像，定位网络会得到$θ_{1}$~$θ_{6}$这6个参数，即 1 $\boldsymbol{\theta}=\left[\begin{array}{lll}\theta_{1} & \theta_{2} & \theta_{3} \\\theta_{4} & \theta_{5} & \theta_{6}\end{array}\right]$ 然后该参数进入第2个组成部分，即网格生成器。网格生成器根据输出图像中的坐标$(x_{i}, y_{i})$找到输入图像中对应的坐标$(x′_{i}, y′_{i})$，从而得到输入图像与输出图像之间的映射关系${\mathit{\boldsymbol{T}}}_{θ}$，为进入第3个组成部分采样器做准备。通过网格生成器得到的映射关系${\mathit{\boldsymbol{T}}}_{θ}$为 2 $\left[\begin{array}{c}x_{i}^{\prime} \\y_{i}^{\prime}\end{array}\right]=\boldsymbol{T}_{\theta}=\left[\begin{array}{lll}\theta_{1} & \theta_{2} & \theta_{3} \\\theta_{4} & \theta_{5} & \theta_{6}\end{array}\right]\left[\begin{array}{c}x_{i} \\y_{i} \\1\end{array}\right]$ 采样器利用线性插值的方法，根据网格生成器得到的映射关系，采样输入图像的像素值放入对应输出图像的坐标中，最后得到输出图像。1.2逆向合成空间变换网络空间变换网络最后得到的输出图像实质上是由变换参数$\mathit{\pmb{θ}}$决定的。定位网络通过旋转、平移和裁剪等一系列的仿射变换，自适应地找到显著性区域，得到变换参数，最终得到感兴趣的部分。但是在这个过程中存在一个问题，输出图像由变换参数决定，如果变换参数出现错误，得到的就不是理想的输出图像。所以为了增强网络的准确性，Lin和Lucey(2017)针对变换参数的优化问题，提出了逆向合成的空间变换网络，网络结构如图 1所示。 图1 逆向合成空间变换网络结构图 Structure diagram of inverse compositional spatial transformation networkFig 1由于输入图像每一次经过空间变换模块得到的输出图像都不同，针对感兴趣区域的提取或多或少有些差距。所以输入图像在经过一次空间变换之后逆向返回，再次通过定位网络，得到一个新的变换参数$\mathit{\pmb{θ}}′$。$\mathit{\pmb{θ}}′$将与第1次得到的变换参数$\mathit{\pmb{θ}}$进行合成，合成后得到的参数一方面代替原来的$\mathit{\pmb{θ}}$成为初始变换参数，另一方面经过网格生成器与采样操作，生成新的图像。新的图像再次返回定位网络，将重新得到的变换参数代替原来的$\mathit{\pmb{θ}}′$。以此往复，参数不断合成。最终得到一个最优的参数，网络学习到准确的感兴趣区域，输出精准的变换图像。2本文方法遥感图像配准由于图像信息的复杂性，缺少针对性的特征提取，同时匹配关系很容易出现错误，因此造成配准精度很难提高。本文提出了一种改进的端到端的配准模型。利用参数合成的空间变换网络进行图像的显著性特征提取，同时利用双向一致性的特征匹配和参数估计，提高模型的准确率和鲁棒能力，得到精确的配准结果。本文将从算法整体结构、改进的特征提取网络、改进的匹配回归部分、改进的参数合成方法以及算法的配准流程这5个部分进行介绍。2.1算法整体结构本文算法结构如图 2所示，主要分为特征提取、特征匹配、参数回归、参数加权合成和得到配准结果5个部分。 图2 本文算法框架结构图 The algorithm framework diagram is presented in this paperFig 21) 特征提取。使用参数合成的空间变换网络提取特征，两幅输入图像经网络得到两个变换参数，再经过网格变换和采样生成，可得到代表各自图像中的关键区域。将两个参数合成再进行变换和采样操作，得到了有共同关键区域的特征图，为特征匹配部分进行了铺垫，提高了网络的鲁棒性。2) 特征匹配。利用特征提取得到的显著性特征图进行互相关操作。根据源特征${\mathit{\boldsymbol{S}}}$中的每个特征向量与目标特征${\mathit{\boldsymbol{T}}}$的特征向量的映射关系，得到从源到目标的相关关系$ST$；然后根据目标特征${\mathit{\boldsymbol{T}}}$中的每个特征向量与源特征${\mathit{\boldsymbol{S}}}$的特征向量的映射关系，得到从目标到源的相关关系$TS$。不同的匹配顺序往往对配准结果有很大的影响，单方向的匹配关系会影响网络结构的平衡性，本文构造了两个方向的匹配结构，保证了配准方向的一致性。图 2中特征匹配部分的$i$和$j$代表了特征点的位置，$k$为索引。3) 参数回归。采取与特征匹配相同的思想，根据匹配关系回归得到双向的参数，同样保证了配准的一致性。回归的双向参数分别为图 2中的源图像到目标图像方向上的参数$\mathit{\pmb{θ}}_{\rm {S→T}}$和目标图像到源图像方向上的参数$\mathit{\pmb{θ}}_{\rm {T→S}}$。4) 参数加权合成。经参数回归得到的两个方向的参数进入损失函数向前传播，从特征提取开始迭代训练，训练出最优的两个参数。由于不同方向的回归参数对配准结果的影响不同，本文将得到的参数加权相加，合成最终的变换参数$\mathit{\pmb{θ}}_{\rm {C}}$。双向一致性加权合成的参数不仅会比单向回归的参数准确率更高，而且能够提高网络的配准精度，达到理想的效果。5) 配准结果。对源图像利用合成后的参数进行仿射变换得到最终的配准结果。2.2利用参数合成的空间变换网络特征提取本文提出了一种参数合成的空间变换网络，能够提取更多且具有针对性的特征，为下一步得到正确的匹配关系进行铺垫。基本框架如图 3所示，主要由参数合成空间变换模块和特征提取模块组成。 图3 参数合成的空间变换网络整体结构图 The whole structure diagram of parametric synthesis' spatial transformation networkFig 32.2.1参数合成空间变换模块空间变换网络利用空间变换模块可以自动地将焦点移动到输入图像中更具有辨别力的部分，从而为后续进行特征提取操作提供有利的条件。利用残差学习中的跳跃连接结构增加空间变换模块中定位网络的深度，将不同尺度的特征信息进行融合，提高定位网络的性能。同时将输入的两幅图像经过定位网络得到的变换参数$\mathit{\pmb{θ}}_{1}$和$\mathit{\pmb{θ}}_{2}$进行合成，得到更精确的变换参数$\mathit{\pmb{θ}}$。合成后的参数进入网格生成器和采样器中得到显著性图像，从而减少了其他不相关信息的干扰。本文设计的定位网络结构图如图 4所示。图中参数$k$表示卷积核的尺寸，$n$表示特征通道的数量，$s$表示步长的大小。Pool为池化层，FC为全连接层。为了方便表示，使用Conv$p\_q$的形式表示第$p$组卷积层中的第$q$个子层，如图 4中Conv1_1表示第1组卷积层中的第1个子层。本文增加16个卷积层加深网络结构，由于网络层数变多，网络可能会出现梯度消失、收敛过慢和过拟合这3种情况。所以每一步卷积操作之后都加入随机失活(Dropout)、ReLU激活函数与批量归一化(batch normalization, BN)进行处理。为防止网络训练过程中出现问题，增加的卷积层按照特征通道数的不同划分为4个全卷积模块，每个全卷积模块包含4个卷积层。为了得到多尺度的特征信息，每两个卷积层和每个全卷积模块恒等跳跃连接，使低层特征与高层特征融合，保证后一层的特征信息比前一层的特征信息更多。 图4 定位网络结构 Location network structureFig 4定位网络分为3个部分：第1个部分由两个具有$7×7$和$3×3$大小滤波器的卷积层和两个$3×3$大小滤波器的池化层构成。首先输入两幅遥感图像，然后通过卷积和池化操作提取浅层特征送入网络的第2个部分。第2个部分分别由特征通道数为64、128、256和512的4个全卷积模块组成，且每个模块中的卷积层具有尺寸大小为$3×3$的卷积核。输入的浅层特征${\mathit{\boldsymbol{x}}}_{l}$与其经过两层卷积操作后得到的特征${\mathit{\boldsymbol{x}}}_{l+1}$通过恒等跳跃连接结构进行相加融合，得到信息更全面的特征${\mathit{\boldsymbol{x}}}_{L_{1}}$。然后${\mathit{\boldsymbol{x}}}_{L_{1}}$再与其经过两个卷积层之后得到的特征${\mathit{\boldsymbol{x}}}_{L_{1}+1}$进行相加融合，从而得到更深层的特征${\mathit{\boldsymbol{x}}}_{L_{2}}$。循环此操作，直至通过4个全卷积模块，最终得到所需的多尺度融合的特征。在这个过程中，参数$l、L_{1}、L_{2}$表示的层数由浅入深，每个全卷积模块内部有如下的关系，即 3 $\boldsymbol{x}_{L_{1}}=\boldsymbol{x}_{l}+\boldsymbol{x}_{l+1}$ 4 $\boldsymbol{x}_{L_{2}}=\boldsymbol{x}_{L_{1}}+\boldsymbol{x}_{L_{1}+1}$ 式中，$l+1$表示$l$经过两个卷积层后的所在层，$L_{1}+1$表示$L_{1}$经过两个卷积层后的所在层。同时每个全卷积模块为了保持每个模块的初始特征能够传递到最后，减少模块内信息丢失的情况，使用跳跃连接维系每个全卷积模块的前端与末端(如图 4中橙色的跳跃连接所示)。第3个部分由1层池化层和3层全连接层构成。通过池化层进行降维，然后利用全连接层高度提纯特征。最后一层全连接层为回归层，作用是生成变换参数。输出参数的维度由空间变换的类型决定，本文采用仿射变换的方法进行空间变换，所以回归得到的参数向量$\mathit{\pmb{θ}}$是由6个参数$x_{1}\sim x_{6}$构成的两行三列矩阵，在齐次坐标中，变换矩阵${\mathit{\boldsymbol{M}}}({\mathit{\boldsymbol{x}}})$可写为 5 $\boldsymbol{\theta}=\boldsymbol{M}(\boldsymbol{x})=\left[\begin{array}{ccc}1+x_{1} & x_{2} & x_{3} \\x_{4} & 1+x_{5} & x_{6} \\0 & 0 & 1\end{array}\right]$ 由于后续的网格变换操作和采样操作对变换参数$\mathit{\pmb{θ}}$有很强的依赖性，所以为了提高$\mathit{\pmb{θ}}$的准确率，将输入的两幅遥感图像经过定位网络回归后得到的两个变换参数$\mathit{\pmb{θ}}_{1}$和$\mathit{\pmb{θ}}_{2}$进行合成，用矩阵${\mathit{\boldsymbol{M}}}({\mathit{\boldsymbol{x}}}′)$表示参数向量$\mathit{\pmb{θ}}_{1}$，矩阵${\mathit{\boldsymbol{M}}}({\mathit{\boldsymbol{x}}}″)$表示参数向量$\mathit{\pmb{θ}}_{2}$，结构如图 3中参数合成空间变换模块所示，合成关系表示为 6 $\boldsymbol{M}(\boldsymbol{x})=\boldsymbol{M}\left(\boldsymbol{x}^{\prime}\right) \times \boldsymbol{M}\left(\boldsymbol{x}^{\prime \prime}\right)$ 得到的合成参数经过网格生成操作和采样操作，找到变换后对应的坐标点进行像素值的填充，最后生成显著性图像作为特征提取模块的输入部分。本文利用双线性插值的方法进行填充，具体表示为 7 $\begin{gathered}\boldsymbol{I}_{i}^{t}=\sum\limits_{n} \sum\limits_{m} I_{n m}^{s} \cdot \max \left(0, 1-\left|x_{i}-m\right|\right) \times \\\max \left(0, 1-\left|y_{i}-n\right|\right)\end{gathered}$ 式中，$I^{t}_{i}$表示输出图像对应的像素值，$I^{s}_{nm}$表示输入图像在通道5处坐标为$(n, m)$处对应的像素值，$x_{i}$和$y_{i}$为输出图像中的第$i$个横坐标点和第$j$个纵坐标点对应输入图像中的横纵坐标。2.2.2特征提取模块通过参数合成的空间变换模块对原始图像进行处理后，本文利用修改后的残差网络结构(ResNet-34)作为特征提取模块，针对上一个模块得到的显著性图片进行特征提取。修改后的结构保留其初始的卷积层和池化层，以及前3个通道数为64、128、256的残差层模块，移除后面所有的残差块与全连接层。在本文方法中，提取第3个残差层模块卷积之后输出的特征图，模块结构图如图 5所示。 图5 特征提取模块结构图 Feature extraction module structure diagramFig 52.3一致性特征匹配和参数回归2.3.1皮尔逊相关一致性匹配特征匹配通过找到两幅特征图之间的相似度，经过匹配函数运算，得到匹配后的相关关系。Rocco等人(2017)根据两幅图像间的语义相似性进行建模，建立具有相关关系的相关向量。向量中元素的顺序不考虑源图像中特征(源特征)的位置，而是基于各个目标图像中特征(目标特征)的绝对坐标，单方向地保证相关向量中的每个元素在某个位置的源特征到目标特征存在相应的映射关系。得到的匹配函数为 8 $C_{\rm {src \cdot trg }}(i, j, z)=\boldsymbol{f}_{\text {src }}(i, j)^{\mathrm{T}} \boldsymbol{f}_{\text {trg }}\left(i_{z}, j_{z}\right)$ 式中，$C_{\rm{src·trg}}∈ {\bf{R}}^{HW×W×H}$表示高度为$H$、宽度为$W$、通道数为$HW$的两幅特征图之间的相关关系；${\mathit{\boldsymbol{f}}}_{\rm {src}}(i, j)$表示在$(i, j)$位置处的源特征向量；${\mathit{\boldsymbol{f}}}_{\rm {trg}}(i_{z}, j_{z})$表示空间扁平化后的目标特征图。为了进一步处理由时间和天气等因素对遥感图像产生的非线性因素，本文采用皮尔逊相关性处理非线性变化，得到的匹配函数为 9 $C_{\mathrm{src} \cdot \mathrm{trg}}^{\text {Person }}(i, j, z)=\frac{\left(\boldsymbol{f}_{\mathrm{src}}(i, j)-\mu_{\mathrm{src}}\right)^{\mathrm{T}}\left(\boldsymbol{f}_{\mathrm{trg}}\left(i_{z}, j_{z}\right)-\mu_{\mathrm{trg}}\right)}{\left\|\boldsymbol{f}_{\mathrm{src}}(i, j)-\mu_{\mathrm{src}}\right\|\left\|\boldsymbol{f}_{\mathrm{trg}}\left(i_{z}, j_{z}\right)-\mu_{\mathrm{trg}}\right\|}$ 式中，$C^{\rm Person}_{\rm {src·trg}}∈ {\bf{R}} ^{HW×W×H}$表示高度为$H$、宽度为$W$、通道数为$HW$的两幅特征图之间皮尔逊相关关系；$μ_{\rm {src}}$与$μ_{\rm {trg}}$分别为源特征图${\mathit{\boldsymbol{f}}}_{\rm {src}}$和目标特征图${\mathit{\boldsymbol{f}}}_{\rm {trg}}$的平均值，其他参数的性质与式(8)相同。上述特征匹配中的相关向量仅仅基于目标特征的绝对坐标，得到从源图像到目标图像单方面的相关关系，并由此相关关系决定配准结果的变换参数。为了增加配准的精度，提高变换参数的准确度，本文受循环一致性原理的启发，基于相关关系中元素的顺序，进行方向一致性的特征匹配，交换源特征与目标特征的输入顺序，添加一条相反方向的相关关系分支，同时得到源特征到目标特征的映射与目标特征到源特征的映射。2.3.2一致性参数回归参数回归网络由一个滤波器为$7×7$大小的卷积层、一个滤波器大小为$5×5$的卷积层和一个全连接层组成。将匹配后得到的两幅关系特征图送入参数回归网络中，利用一致性原理进行双向回归后得到了两个方向的参数。然后计算两个参数输入前与输入后的损失，不断迭代优化。采取网格距离函数$L_{\rm {grid}}(\mathit{\pmb{θ}}, \mathit{\pmb{θ}}^{\rm gt})$作为损失函数，即 10 $L_{\text {grid }}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {gt }}\right)=\frac{1}{N} \sum\limits_{i, j=1}^{N}\left[d\left(T_{\theta}\left(x_{i}, y_{j}\right), T_{\theta{\rm{g t}}}\left(x_{i}, y_{j}\right)\right)^{2}\right]$ 式中，距离$d$为变换前图像的真实情况的点$T_{θ^{{\rm gt}}}(x_{i}, y_{i})$与变换后的输出图像上的点$T_{θ}(x_{i}, y_{i})$之间的平方差。$\mathit{\pmb{θ}}^{{\rm gt}}$为真实情况的参数，$\mathit{\pmb{θ}}$为变换后输出的参数。网格点的总数量为$M$，$i, j∈M$, 将网格距离函数看成优化问题，经过训练可以最小化目标函数值，从而减小了真实情况参数与输出参数之间的差异，得到多次训练后效果最好的一组匹配参数值。2.4加权合成匹配参数匹配是有方向性的，双向匹配相较于单向匹配能够降低错误匹配的概率(刘焕敏等，2009)，正确指导配准。宁静静和孔令德(2012)提出了一种双向互匹配的方法，从源图像到目标图像进行匹配，同时从目标图像到源图像进行匹配，有效提高匹配精度。Ji等人(2015)通过实验设置双向匹配对的权重并合成，得到最优的配准结果。本文通过实验对得到的双向参数进行加权。$\mathit{\pmb{θ}}_{\rm {S→T}}$为从源图像到目标图像方向上回归得到的参数，$\mathit{\pmb{θ}}_{\rm {T→S}}$为从目标图像到源图像方向上回归得到的参数。考虑到配准是利用源图像与目标图像之间的最优映射关系，对源图像向目标图像进行变换的过程。所以，参数$\mathit{\pmb{θ}}_{\rm {T→S}}$用来提高匹配的精度，剔除错误匹配。为了保持参数在方向上的一致性，变换匹配方向，取逆$(\mathit{\pmb{θ}}_{\rm {T→S}})^{－1}$，得到从源图像到目标图像方向上的参数$\mathit{\pmb{θ}}_{\rm {S→T}}$。在此处，参数为矩阵参数，每个参数都含有仿射变换所需的6个参数(平移、旋转等)，求逆则是求原矩阵的逆矩阵。此过程可用数学公式表示为 11 $\boldsymbol{A} \times \boldsymbol{X}_{\mathrm{T}}=\boldsymbol{X}_{\mathrm{S}}$ 12 $\boldsymbol{A}^{-1} \boldsymbol{A} \times \boldsymbol{X}_{\mathrm{T}}=\boldsymbol{A}^{-1} \times \boldsymbol{X}_{\mathrm{S}}$ 13 $\boldsymbol{X}_{\mathrm{T}}=\boldsymbol{A}^{-1} \times \boldsymbol{X}_{\mathrm{S}}$ 式(11)表示从目标图像到源图像(target→source)的变换过程。${\mathit{\boldsymbol{X}}}_{\rm {S}}$表示源图像(source)，${\mathit{\boldsymbol{X}}}_{\rm {T}}$表示目标图像(target)，${\mathit{\boldsymbol{A}}}$为变换关系，也就是矩阵参数$\mathit{\pmb{θ}}_{\rm {T→S}}$。式(12)对式(11)左右两边同乘一个${\mathit{\boldsymbol{A}}}^{-1}$，即变换关系矩阵的逆矩阵，得到式(13)，表示若想得到source→target的变换，则求${\mathit{\boldsymbol{A}}}^{-1}$，即$(\mathit{\pmb{θ}}_{\rm {T→S}})^{－1}$。计算时，将其转变成齐次坐标的形式，即 14 $\boldsymbol{\theta}_{\mathrm{T} \rightarrow \mathrm{S}}=\left[\begin{array}{ccc}\beta_{1} & \beta_{2} & \beta_{3} \\\beta_{4} & \beta_{5} & \beta_{6} \\0 & 0 & 1\end{array}\right]$ 15 $\left(\boldsymbol{\theta}_{\mathrm{T} \rightarrow \mathrm{S}}\right)^{-1}=\left[\begin{array}{ccc}\beta_{1} & \beta_{2} & \beta_{3} \\\beta_{4} & \beta_{5} & \beta_{6} \\0 & 0 & 1\end{array}\right]^{-1}$ 式中，$β_{1}$~$β_{6}$为仿射变换的变换参数。对参数$\mathit{\pmb{θ}}_{\rm {S→T}}$与参数$(\mathit{\pmb{θ}}_{\rm {T→S}})^{－1}$加权合成，此过程遵循如下的公式进行合成，即 16 $\boldsymbol{\theta}_{\mathrm{C}}=\mu_{1} \cdot \boldsymbol{\theta}_{\mathrm{S} \rightarrow \mathrm{T}}+\mu_{2} \cdot\left(\boldsymbol{\theta}_{\mathrm{T} \rightarrow \mathrm{S}}\right)^{-1}$ 式中，$\mathit{\pmb{θ}}_{\rm {C}}$是加权合成后的参数，$μ_{1}$为源图像到目标图像方向上的权重，$μ_{2}$为目标图像到源图像方向上的权重。最后根据合成后的参数对源图像进行仿射变换，得到配准结果。2.5算法配准流程本文算法的流程图如图 6所示，主要步骤有：1)利用参数合成的空间变换网络作为骨干网络进行显著特征提取; 2)利用皮尔逊相关性建立匹配关系，进行一致性匹配; 3)根据匹配后的结果进行一致性参数回归; 4)利用网格距离损失函数不断更新网络参数，指导回归准确的变换参数; 5)加权合成一致性回归后的参数，采样得到最终的配准结果。 图6 本文流程图 Flow chart of this paperFig 63实验结果与分析3.1实施细节3.1.1数据集与实验环境Aerial Image Dataset数据集是Park等人(2020)在不同的时间(2015年、2017年、2019年)利用不同的传感器(Landsat-7，Landsat-8，WorldView和QuickBird)从不同的角度拍摄的9 000组韩国区域遥感航空影像对，然后对每组图像的相同区域进行中心裁剪和仿射变换操作生成的复杂数据集。本文使用调整后的Aerial Image Dataset数据集(Park等，2020)训练网络，根据计算机运行速度，依据模拟环境(内存)，将训练集调整为3 000对图像，测试集调整为100对图像，进行定量与定性的测试评估。并根据经验值对网络参数不断调试，将网络的学习率调整为$4×10^{－4}$，批处理大小设定为32，权重衰减为0，动量参数为0.9。测试集同时使用Corpus数据集、VIS-NIR数据集(王爽等，2018)和SUIRD(small UAV image registration dataset)数据集(Gong等，2020)上的图像进行配准测试。其中，Corpus数据集有8种类别的图像，每种类别有两种不同时间段拍摄的卫星图。VIS-NIR数据集包括9组多源图像，为可见光图像与红外图像。SUIRD数据集针对不同类别的场景从不同的角度拍摄，捕获到多种类别的小型无人机图像，每种类别有60对图像，共1 200幅图像。实验使用Pytorch作为深度学习框架，硬件平台具有Intel core i5(2.3 GHz)的处理器和8 GB的内存，在操作系统为Mac OS的PC机上运行实现。3.1.2评估指标本文使用6种图像配准中常用的评估指标作为定量衡量算法性能的依据，分别为关键点正确估计的比例(percentage of correct keypoints，PCK)、平均绝对误差(mean absolute error，MAE)、均方根误差(root mean square error，RMSE)、网格损失(loss)、平均网格损失(average loss)和配准所用的时间(time)。使用棋盘格图和重叠图作为定性衡量配准性能的依据。1) 关键点正确估计的比例。表示为在一定的阈值范围内，计算配准后图像检测到的关键点与其对应的人工标记点(ground truth)间的归一化距离小于设定阈值的比例，此比例即为关键点正确匹配的比率，即 17 $P C K=\frac{1}{N} \sum\limits_{i=1}^{N}\left[d\left(T_{\mathrm{fin}_{i}}\left(p_{s_{i}}\right), p_{t_{i}}\right)\alpha \cdot \max (h, w)\right]$ 式(17)为正确检测到关键点的比率，$N$表示测试图像对的总数，$T{\rm {fin}}_{i}$是模型中的最终变换参数，$p_{{s}_{i}}$是第$i$个图像对经过变换得到的源关键点，$p_{{t}_{i}}$是第$i$个图像对中真实标注的目标关键点，而$α·$max$(h, w)$表示最大阈值范围。若检测到的点与人工标注的点之间的距离$d$小于阈值，则表明是检测正确的点。在高为$h$，宽为$w$的图像里，如果$α$越大(系数$α$不会超过1)，阈值范围越大，则可以更加全局地测量配准情况，一般来说$α$取0.1更为合适。该指标的值越大，配准精度越高。2) 平均绝对误差。用来衡量图像配准后的关键点与人工标记的关键点之间的绝对误差，平均绝对误差越小，准确率越高。3) 均方根误差。用来衡量预测后的关键点与人工标记的真实值之间的偏差，均方根误差值越小，配准效果越好。4) 配准所用时间。输入两幅图像经过训练好的网络处理得到的配准结果所用的时间，所用的时间越少越好。5) 网格损失是每次迭代后的损失。针对的是批处理块的损失，平均网格损失是网格损失里的一个参数，也是计算每次迭代后的损失，但代表的是图片的损失均值；两者的计算方式不一样，但是都能表示网络模型性能的优劣和精度大小。损失越小，配准精度越高。3.1.3网络训练细节网络具体的训练过程如下所示：1) 重构输入的源图像${\mathit{\boldsymbol{I}}}_{\rm {S}}$与目标图像${\mathit{\boldsymbol{I}}}_{\rm {T}}$的大小，通过中心裁剪调整两幅输入图像。输入的图像来自Aerial Image Dataset数据集，由于硬件平台限制原因，选其中3 000对图像进行训练。2) 经过参数合成空间变换网络进行向前传播，训练后得到从源图像到目标图像的变换参数$\mathit{\pmb{θ}}_{\rm {S→T}}$，以及从目标图像到源图像方向的变换参数$\mathit{\pmb{θ}}_{\rm {T→S}}$。3) 利用网格损失函数$loss{\rm {grid}}$计算损失，同时使用Adam优化算法最小化损失后反向传播，不断更新网络的权重，直至epoch 40次后完成训练。训练过程需要花费约4天的时间。3.2本文实验结果本文在Aerial Image Dataset数据集上进行配准测试，如图 7所示。第1组图像来源于数据集中的港口图，经配准后，港口方向和港口形状大致相同。第2组图像来源于数据集中的城市图，图中主要有道路和规划好的住宅区域，配准时以道路为主，住宅区域为辅，得到较好的配准效果。第3组图像来自数据集中的田野图，图中重点区域为耕地，其次是田间的道路，最后是田间的树林区域，经本文方法的配准后发现，图中重点区域匹配效果较好，配准后的图像与目标图像较为一致。第4组是海湾图，主要配准相应的海湾的边界区域。第5组是源图像带雾的乡村图，而目标图像不带雾，经本文方法配准后即使是有白雾干扰，配准结果依然准确。第6组是不同时期拍摄的岛屿图，本文算法对岛屿边界的配准效果较好，整体配准精准。 图7 本文实验结果图 Experimental results in this paperFig 7本文算法验证了Aerial Image Dataset数据集上图片，各类复杂图像的配准结果都较为精准，取得不错的配准效果。3.3算法性能对比将本文方法与SIFT(scale-invariant feature transform)、SURF(speeded up robust features)、CNNGeo(convolutional neural network architecture for geometric matching)、CNN-Registration(multi-temporal remote sensing image registration)和RMNet(robust matching network)这5种方法进行比较。SIFT(Lowe，2004)和SURF(Bay等，2006)这两种方法是经典的传统算法，而CNNGeo(Rocco等，2017)、CNN-Registration(Yang等，2018)和RMNet(Kim等，2019b)这3种方法是新方法。其中CNN-Registration属于半传统和半深度学习结合的算法，CNNGeo和RMNet属于端到端的深度学习的方法。为了确保对比的一致性，使用本文的数据集对上述方法进行训练和测试。同时，对比算法的测试场景涉及多种类型的图像，保证了对比的可靠性。3.3.1定性比较针对多时相、多视角、多源和复杂情况下拍摄的图像(包含多时相、多视角和多源3种情况)进行配准，并将配准结果与经典的传统算法和近几年提出的基于深度学习的方法进行对比，结果如图 8所示。 图8 不同算法在多种情况下的定性比较 Qualitative comparison of different algorithms under various conditionsFig 8((a)source images; (b)target images; (c)SIFT; (d)SURF; (e)CNNGeo; (f)CNN-Registration; (g)RMNet; (h)ours) 图 8中的第1行是针对两个时间段上河流变化的配准应用，图像来自Corpus数据集；第2行是针对从不同的角度拍摄的山体照片进行分析，图像来自SUIRD数据集；第3行是利用不同的传感器对无人机拍摄的光学图与热红外图进行配准，图像来自VIS-NIR数据集；第4行则是对包含了以上3种复杂情况的图像进行处理，拍摄时间不同，拍摄的角度不同，同时拍摄所用传感器也不同，图像来自Aerial Image Dataset数据集。观察实验对比结果，本文方法无论是在单独某一个方面变化(如图 8中前3行)的情况还是针对多种复杂情况(如图 8中第4行所示)，可以看出本文计算图像都可以取得比较好的效果。但是本文方法与RMNet(Kim等，2019b)的配准结果从整体上来看差别不大。为了进一步定性地对比两种方法，对复杂情况下的图片进行配准，同时采用棋盘格图和重叠图对配准细节进行观察，如图 9所示。 图9 RMNet方法与本文方法在复杂情况下的定性比较 A qualitative comparison between RMNet method and the present method in complex casesFig 9((a)source images; (b)target images; (c)RMNet registration results; (d)our registration results; (e)RMNet checkboards; (f)our checkboards; (g)RMNet overly; (h)our overly) 图 9为RMNet方法和本文方法针对复杂图像的配准情况进行分析。图中观察3组来自Aerial Image Dataset数据集的复杂图像，对配准后得到的棋盘格图与重叠图存在偏差的位置用红色方框突出。可以观察到两种配准方法各自存在不同的缺点，但是本文算法(图 9(f)(h))的红色偏差框数量要比RMNet方法少，且RMNet方法重点区域的边缘配准细节的误差较大，所以从定性的角度在视觉上分析本文算法配准效果较好。3.3.2定量比较本文利用上述提到的6个评价指标对算法性能进行定量评估。首先利用PCK评估在Aerial Image Dataset测试集上计算各种方法的正确点匹配比率，根据比率推断各种方法配准的精度。评估结果如表 1所示。 表1 不同算法的PCK评估结果 方法 算法 α为0.1的PCK 传统算法 SIFT 0.548 SURF 0.151 半传统+半深度学习算法 CNN-Registration 0.374 端到端深度学习算法 CNNGeo 0.611 RMNet 0.893 本文 0.910 PCK evaluation results of different algorithmsTable 1 加粗字体为最优结果。表 1中可以看出端到端深度学习的性能优于其他两种方法，且PCK的值要比传统算法和半传统算法的值高30 % 以上。其中CNNGeo方法的PCK值可以达到61.1 %，RMNet方法的PCK精度达到了89.3 %，而本文算法的PCK精度达到了91.0 %，分析可知本文算法的配准精度略优于RMNet方法。接着利用在Aerial Image Dataset数据集上训练40次的网格损失折线图和平均网格损失折线图对比CNNGeo、RMNet和本文算法这3种端到端学习方法的配准精度。观察图 10中每次训练后得到的网格损失变化以及图 11中每次训练后得到的平均网格损失变化，发现RMNet方法和本文方法的网格损失与平均网格损失下降的梯度最大且训练到最后的损失值最小，说明这两种方法的网络性能与算法精度明显比CNNGeo方法要好。 图10 网格损失随迭代次数变化的折线图 A line graph in grid loss with the number of iterationsFig 10 图11 平均网格损失随迭代次数变化的折线图 A line graph in average grid loss with the number of iterationsFig 11由于在图 10和图 11中本文方法与RMNet方法的网格损失与平均网格损失变化迭代图的差距很小，观察不出明显变化，所以将无明显梯度变化的后30次迭代单独制作柱状对比图，观察损失变化。观察图 12和图 13柱状对比图，蓝色代表RMNet方法的损失，红色代表本文方法的损失，发现无论是网格损失还是平均网格损失，本文算法的损失都比RMNet方法的损失值低。为了更容易看出变化，两幅图中的蓝色虚线为RMNet方法的损失变化趋势，红色虚线为本文方法的损失变化趋势，红色趋势线处于蓝色趋势线下方，表明本文方法的损失更低、精度更高。 图12 最后30次迭代的网格损失对比图 Grid loss comparison for the last 30 iterationsFig 12 图13 最后30次迭代的平均网格损失对比图 Average grid loss comparison for the last 30 iterationsFig 13然后利用MAE、RMSE和时间这3个指标，针对图 8中的4组图像在测试集上进行评估。如表 2所示，MAE和RMSE评估最优的结果都是在RMNet方法和本文方法中出现(表 2中加粗字体的结果)，且每组的差距很小。但是对于时间这一指标，每次都是本文方法用时最少，且对比于精度很高的RMNet方法，第1组多时相实验测试时间提高了4.06 s，第2组多视角的时间提高了3.94 s，第3组多源的时间提高了3.69 s，第4组复杂图像的测试时间提高了4.12 s。如果在实际应用中，将大大提高运行效率。 表2 4种情况下不同算法的定量评估 实验 算法 MAE RMSE 时间/s 多时相 SIFT 26.12 30.34 10.21 SURF - - 8.13 CNN-Registration 30.28 37.67 49.26 CNNGeo 42.60 51.34 2.03 RMNet 19.16 22.54 5.03 本文 20.19 23.65 0.97 多视角 SIFT - - 9.25 SURF - - 10.45 CNN-Registration - - 50.12 CNNGeo 27.65 32.49 1.88 RMNet 17.45 19.20 4.98 本文 17.30 19.48 1.04 多源 SIFT - - 11.04 SURF - - 10.25 CNN-Registration - - 52.13 CNNGeo - - 2.16 RMNet 2.895 3.254 4.57 本文 2.936 3.186 0.88 复杂情况 SIFT 35.78 52.34 9.38 SURF - - 9.27 CNN-Registration - - 50.85 CNNGeo 22.38 25.90 1.56 RMNet 19.50 21.56 5.39 本文 21.23 24.87 1.27 Quantitative evaluation of different algorithms under four conditionsTable 2 加粗字体为各实验中最优情况; “-”表示配准结果失真或标记点不在图像内，结果差别过大不纳入计算。3.4创新点有效性分析3.4.1消融实验本文将参数合成的空间变换网络(简写为STN)、一致性匹配回归(简写为consistence)和加权合成匹配参数(简写为concentrate)3个创新点，分别加入基准配准框架CNNGeo方法中。其中参数加权合成考虑以源图像到目标图像的变换参数为主，目标图像到源图像的变换参数为辅，设置不同的系数分别为0.5和0.5、0.6和0.4、0.7和0.3、0.8和0.2以及0.9和0.1(系数1和0代表单向配准为CNNGeo使用的方法，本文使用的系数为0.7和0.3)，同时利用PCK指标在Aerial Image Dataset数据集上定量分析加入每个创新点后网络整体的配准性能和配准精度。考虑到RMNet方法是基于配准框架CNNGeo方法的改进，保留了原有特征提取网络，对特征匹配部分有较大改动。本文利用STN模块代替RMNet网络的特征提取模块，分别与原有RMNet网络和本文方法进行特征提取部分和特征匹配部分的性能对比。基于CNNGeo方法的消融分析结果如表 3所示，基于RMNet方法的消融分析结果如表 4所示。 表3 基于CNNGeo的消融分析表 方法 α为0.1的PCK CNNGeo 0.611 CNNGeo+STN 0.657 CNNGeo+consistence 0.650 CNNGeo+concentrate(0.5+0.5) 0.712 CNNGeo+concentrate(0.6+0.4) 0.682 CNNGeo+concentrate(0.7+0.3) 0.752 CNNGeo+concentrate(0.8+0.2) 0.674 CNNGeo+concentrate(0.9+0.1) 0.624 Ablation analysis table based on CNNGeoTable 3 表4 基于RMNet的消融分析表 方法 α为0.1的PCK RMNet 0.893 RMNet+STN 0.898 本文 0.910 Ablation analysis table based on RMNetTable 4表 3第1行为基准框架CNNGeo在本文测试集中的PCK精准度；第2行表示使用本文提出的参数合成的空间变换网络代替原有算法的特征提取网络，配准后PCK精度提高了4.6 %；第3行为利用一致性匹配和回归代替原有结构的单向匹配和回归，PCK的值提高了3.9 %；从第4行到最后一行对比参数的加权系数，结果表明本文使用的0.7和0.3的系数精度最高，且比基准框架CNNGeo方法提高了14.1 % 的精度。表 4第2行本文提出的参数合成的空间变换网络代替原有算法的特征提取部分，配准后PCK精度提高了0.5 %，表明本文方法的特征提取网络的性能优于RMNet方法；第3行为本文方法，与第2行加入STN模块的RMNet方法比较，PCK精度提高了1.2 %，表明在保持特征提取网络一致的情况下，本文特征匹配部分的性能优于RMNet方法。实验结果证明本文提出的创新点都能在一定程度上提高配准精度，证明了创新点的有效性。3.4.2特征提取分析本文提出的参数合成空间变换网络在经过空间变换模块后会进入特征提取模块，最后得到一组特征提取后的特征图。将本文网络得到的特征图与CNNGeo、CNN-Registration两种算法经特征提取后的特征图进行对比(RMNet与CNNGeo使用同种特征提取网络)，如图 14所示。 图14 不同算法特征图对比 Comparison of feature graphs of different algorithmsFig 14((a)origin; (b)CNN-Registration; (c)CNNGeo; (d)ours) 从图中可以看出，本文方法提取的感兴趣区域的特征较为明显，干扰特征较少，而其他两种方法都存在一些无关特征。3.4.3一致性匹配回归和参数加权合成分析为了观察不同算法间点阵的匹配情况，将传统算法SIFT、基于端到端深度学习算法RMNet与本文方法进行对比。使用匹配线图来观察对应特征点的匹配情况，如图 15—图 17所示。 图15 SIFT算法匹配图 SIFT algorithm matches graphsFig 15 图16 本文方法匹配图 The method in this paper matches the graphFig 16 图17 RMNet方法匹配图 RMNet method matches graphsFig 17首先对比SIFT算法和本文算法的匹配图。观察图 15可知，SIFT方法的匹配结果局限在右边的港口处。而本文算法，由图 16可知，从整体上进行匹配，兼顾到了每一个港口，配准的结果会比SIFT更精准。对比RMNet方法与本文方法的匹配图，观察图 16与图 17发现两种方法从整体上进行了匹配且兼顾到了每一个港口，无法区分出两者的区别。所以通过计算本文算法与RMNet方法的正确匹配点的数目和正确匹配对的数目来对比两种算法的匹配精度，如表 5所示。 表5 3种方法匹配结果 方法 正确匹配对的数目 正确匹配点的数目 SIFT 32 128 RMNet 80 352 本文 86 420 Matching result of the three methodsTable 5RMNet方法依赖于从源图像到目标图像单方向的匹配关系，并通过单方向的回归得到一个用于配准的匹配参数；本文方法增加一条分支，通过双向一致性的匹配关系，回归得到两个方向的参数，然后对两个参数加权合成得到一个用于配准的匹配参数；表中本文方法的正确匹配对数目比RMNet方法多了6对，而正确匹配点的数目多了68个，通过这两种方法的匹配精度对比，本文一致性的加权合成匹配回归方法的精度会略高于RMNet单向匹配回归方法。为了进一步证明本文提出的一致性匹配回归和参数加权合成的有效性，将此方法运用到传统SIFT方法中，对比双向加权的SIFT方法与传统的单向匹配SIFT方法，如图 18所示。 图18 单向SIFT与双向加权SIFT配准结果对比图 Comparison of registration results between unidirectional SIFT and bidirectional weighted SIFTFig 18((a)source images; (b)target images; (c)single SIFT; (d)bidirectional weighted SIFT) 图 18第2组和第3组的图像来自于Aerial Image Dataset数据集，第3组图像为多时相图像对(Yang等，2018)。对以上3组图像对进行配准，观察配准结果可知，双向加权SIFT方法明显比单向SIFT方法的配准效果更好，从而证明了一致性匹配回归和参数加权合成的有效性。4结论针对遥感图像存在特征提取困难的问题和单向匹配关系指导配准引发误差的现象，提出基于参数合成的空间变换网络一致性遥感配准的解决方法。在本文算法中，采用改进的参数合成空间变换网络作为特征提取网络，对显著性区域提取特征，具有较强的针对性和鲁棒性。同时，本文采用双向匹配关系回归双向参数进行一致性的匹配和回归，将双向参数加权合成指导配准，增强了配准的可靠性，提高了配准的精度。实验结果表明，在公开的数据集上，本文使用6种评估指标将本文方法与其他五种具有代表性的遥感配准算法进行对比，不仅在整体的配准效果上优于其他方法，而且在细节的配准效果上也能取得优异的结果。同时可以适用于多种类型的遥感图像配准，实现多场景多模态的配准。本文方法也存在提升的空间，今后准备从以下两方面进行改进：1)利用无监督学习的方法对网络进行改进，减少人工干预；2)将研究领域着眼于3D遥感图像的配准，具有更强实用性和挑战性。