基于特征约束与光流场方法的遥感图像配准

刘 璐，陈 甫，刘建波＋

（1.中国科学院 遥感与数字地球研究所，北京100094；2.中国科学院大学，北京100049）

0 引 言

大部分在轨运行的对地观测光学卫星同时提供高分辨率 全色图像（high spatial resolution panchromatic，PAN）和低分辨率多光谱图像 （low spatial resolution multispectral，LRMS），如LandSat系列、Spot系列及GeoEye系列等。其中，一些卫星的PAN 和LRMS由于成像模式不同，存在复杂的非刚性变形。例如资源三号卫星，采用全色和多光谱2个相机获取PAN 和LRMS，成像时间有差异、成像角度不同，因此获取的PAN 和LRMS存在复杂的非刚性变形，即使经过系统级几何纠正后依然存在显著位移。为了使用PAN 的几何信息和LRMS的光谱信息，进行目视解译、分类等图像处理操作，需要对上述图像进行高精度的配准。

图像配准是对取自不同时间、不同传感器或不同视角的同一景物的2 幅或多幅图像进行匹配和叠加的过程［1］，是遥感图像应用的基础。图像配准的方法有使用遥感图像处理软件 （如ERDAS等）的人工手动配准；基于特征提取的配准，如SIFT［2］；基于区域的配准，如互信息法［3，4］；基于相位的配准，如基于傅里叶的方法［5］。这些配准方法手动或自动地获取待配准图像上的同名特征，并通过假设图像满足一种空间变换模型，对整幅图像进行空间变换。因此，上述的方法也可以归类为基于空间变换的图像配准。另一类配准方法为基于物理模型的配准，如光流场方法［6］，它不假设空间变换模型，以物理学中的热传导等理论为依托构建配准方程，该类方法在遥感领域的研究较少。

光流场方法起源于计算机视觉领域，主要用于运动物体检测。由于运动物体的光流场与图像配准的位移场具有一定的相似性，因此光流场的基本原理也可以用于图像配准。主流的光流场方法是基于Horn等［6］提出的变分法，该方法基于像素灰度值的局部梯度模型和整体平滑性假设。针对上述方法的不足，一系列的改进和延伸模型也被提出，如文献 ［7］利用梯度恒定解决光照问题等，为了解决大位移问题，Lucas等［8］提出了金字塔方法，Brox等［9］提出了描述符匹配等。

针对同一卫星的PAN 和MS的变形复杂、难以用简单的空间变换模型表示这一特点，本文使用光流场方法对上述图像进行配准。为了提高光流场配准的鲁棒性和精度，本文在Brox等［9］提出的方法基础上，选择广泛适用于遥感图像的SIFT［2］特征描述符作为约束，并采用格网分割分配SIFT 光流，对PAN 和LRMS进行高精度亚像素的配准，实验结果表明本文算法可以处理大位移图像、减少误匹配光流、提高光流场的精度。

1 算法原理

光流场方法最先起源于计算机视觉领域，将该方法应用到遥感图像配准中，需要解决一系列遥感图像特有的问题以满足模型的适用条件。首先，由于光流场模型一般用泰勒级数展开并省略高阶项，因此不适用于大位移图像配准，一般位移不超过2个像素，目前大部分的解决方案是选用由粗到细的策略，即影像金字塔方法，本文在采用金字塔方法的基础上同时采用SIFT 特征作为辅助信息，解决大位移问题。其次，为了满足亮度一致性假设，LRMS要通过融合仿真PAN 图像。最后，偏微分方法的一个基本条件是图像连续，因此待配准图像需要平滑处理。

1.1 图像预处理

首先，由于PAN 和LRMS 的分辨率不同，需要对PAN 图像进行降采样。

其次，光流场方法是基于图像灰度的方法，亮度不变假设 （假设同一物体的灰度不随时间变化而改变，在运动中保持辐射亮度一致）是光流场方法的基本假设，为了配准LRMS 和PAN 图像，需要将LRMS 多波段融合模拟PAN 图像，融合后的LRP （low resolution panchromatic）作为计算光流场的源数据。本文采用无常量的线性回归波段拟合法求解LRP，线性回归的公式如下

式中：n——多光谱图像融合的波段数，CK——线性回归系数，CK的求解方法可以基于SRF［10］（spectral response function），也可以由像元灰度值直接求解，本文采用后者求解。

1.2 基于SIFT约束的光流场方法

1.2.1 光流场变分模型

光流场方法是基于图像灰度的算法，它的基本假设的亮度不变假设，由此可得数据项公式

其中，：＝ （x，y）T定义为二维微分算符。仅有上述公式不足以求解光流场需要添加正则项约束。本文采用全局平滑约束，公式如下

全局平滑约束假设整幅图像光流场平滑变化，对光流进行基于全局的平滑。

1.2.2 基于SIFT 约束的光流场模型

上述变分法光流场模型，是基于图像灰度的方法，以每个像素为基准进行配准处理，这种处理方法虽然可以解决复杂变换模型的问题，然而当位移较大或者灰度变化较大时，易造成误匹配。由粗到细的策略诚然可以在一定程度上解决大位移问题，提高算法的稳定性，但是该方法在高层次金字塔上的光流场依赖于底层金字塔，而图像在底层金字塔分辨率降低、图像信息减少，底层的光流场精度降低。同时，考虑到PAN 和LRMS 图像的特殊性，即图像存在云雾等遮挡、PAN 和LRP图像灰度不一致等现象，本文受Brox［9］等人的启发，采用基于SIFT 约束的光流场方法。

SIFT 特征是遥感图像领域广泛采用的特征描述子，具有旋转和位移不变性［2］。可以作为一种附加图像信息添加到光流场方法中。然而，特征匹配有其弊端，它是一种离散的方法，并不能够提供亚像素的精度，因此在光流场模型中并不能作为最高层金字塔的约束条件。综上所述，可以采用SIFT 作为约束项，提高光流场在金字塔底层的精度与稳定性。SIFT 作为约束的原理是，用光流场方法求得的位移应与SIFT 求得的位移相等，因此公式如下所示

式中：uSIFT，vSIFT——由SIFT 约束所得的x，y方向上的位移。τ——正则项的权重，Q——权值，Layer——金字塔的层数，1——最高层，2——次高层，以此类推。GridMax将在下面阐述。

SIFT 特征匹配得到的是离散的同名点，在基于特征的图像配准中，可以利用空间变换模型对待配准图像纠正。本文采用一种网格分割的方法分配SIFT 光流信息。该方法将图像分为GridSize像素大小的块，一般情况下GridSize为奇数，可取 ［1，3，5，9…2n－1］，n是大于1的整数。在每个网格内，由于网格面积较小可以假设网格内的像素点满足某一种空间变换模型，结合遥感图像的特点，仿射变换较为适宜，仿射变换的公式如下

式中：X2，Y2——待求点，x1，y1——已知点，（a，b，c，e，f，g）——仿射变换系数。GridMax是落在每个网格内的同名特征点数量，当采用仿射变换模型时，GridMax最小为3，当GridMax大于3时，可以使用最小二乘法计算，公式如下

式中：A——仿射系数矩阵，X2，X1——已知点矩阵，，——转置矩阵和逆矩阵。计算出仿射变换系数矩阵A后，所有在网格内的点可以近似求出位移量。遥感图像分辨率较小，且图像配准的目的是最大限度的匹配两幅图像上的一致地物，如遥感图像中有运动物体，则原则上较好的匹配算法可以自动甄别运动物体，避免其对其它地区造成配准干扰。本文采用格网分割的方法，正是基于上述考虑，光流场算法是一种像素级的配准方法，容易出现误匹配点，且由于全局约束正则项的存在，误匹配点会影响到其它像素的配准结果。SIFT 特征提取到的同名特征点，可能会存在两幅图像上运动的点或者误匹配的点，如果直接使用在光流场方法中，会影响配准精度和鲁棒性。经过格网分割后，每个格网内至少要存在3个以上的同名特征点才会添加SIFT 约束，因此减少了SIFT 误匹配或运动点对本文配准算法的影响，同时最大限度的利用SIFT 信息，提高配准的鲁棒性。

经过网格分割与特征点计算，可以在特征点满足求解条件的区域得到uSIFT，vSIFT约束，当网格较小时，可以将Q 设置较大值；反之则为Q 设置较小值。综合上述分析，则总能量函数为式 （10），其中λ，α表示数据项和正则项的权重

式 （10）的求解过程可以利用SOR 法迭代求解增量duq＋1与dvq＋1，结果如下所示

其中ω ∈ （0 ，2） 为SOR 的松 弛 因 子，s为SOR 的 迭 代 次数，k表示求解u，v的迭代次数，称之为外层迭代uk＋1＝uk＋duk，vk＋1＝vk＋dvk，q表示求解duk，dvk迭代次数，称之为内层迭代，具体推导过程感兴趣的读者可以参考文献 ［7］。

2 算法实现与实验分析

2.1 算法流程

算法流程如图1所示。首先，对2幅图像进行预处理，生成PAN 和LRP。其次，构造影像金字塔：当Layer＞1时，对每层金字塔影像分别提取SIFT 特征，并进行匹配；构造图像网格，统计每个小网格图像中的SIFT 特征个数，并计算每个小块的仿射变换系数，以此求得小块内每个点的偏移量。第三，计算每层的光流场，其中高层光流场的初值由低层光流场初始化。第四，输出配准后的图像。

图1 基于特征约束的光流场算法流程

2.2 实验分析

本文采用国产高分辨率卫星资源三号 （ZY－3）作为实验数据。ZY－3卫星是中国在2012年发射的民用高分辨率光学传输型立体测图卫星，搭载的前、后、正视和多光谱相机，正视相机空间分辨率2.1 m，多光谱相机空间分辨率6m，本文采用正视影像和多光谱影像进行配准。精度评价采用主观与客观相结合的方法：主观上，通过人工目视判读、差值图像法等观察图像配准效果；客观上，选择相关系数 （correlation coefficient，CC）、均方根误差 （rootmean－square error，RMSE）等统计特性评价配准精度。采用ENVI手动配准作为对比实验。

本文选择大连地区影像为实验数据，数据获取时间为2012／06／26。以LRP作为参考图像，待配准图像为PAN，两者之间的位移横轴约5 个像素，纵轴约175 个像素。2幅图像均为L1级产品，经过辐射纠正，没有几何纠正。

对图2中的图像采用本文方法进行配准。实验中参数根据经验值分别设置为：λ＝1，α＝50，β＝0.001，ω＝1.9，β＝0.8，k＝15，s＝300，q＝1，Q＝1，GridMax＝3，金字塔缩放比例为0.5，平滑方法为高斯平滑，高斯平滑因子为0.8，插值方法为双线性内插。SIFT 特征提取采用文献 ［2］的方法，参数与文献 ［2］相同。

利用ENVI软件人工选择的控制点时，尽量使控制点分布均匀，均方误差控制在1个像素以内，变换模型选择一次多项式，插值方法为双线性内插具体控制点分布和误差如图2所示。

为了便于评价配准精度，选择上述图像位于中心的3000＊3000像素块进行精度评价。如表1所示，RMSE－1、CC－1分别表示配准后均方根误差与相关系数，RMSE－0、CC－0表示配准前的均方根误差与相关系数。由表1 可知，采用本文配准后的相关系数和均方根误差均较ENVI手动配准有较明显的提高

表1 配准精度

图3 山区配准图像细节

由于拍摄角度不同，高程较高地区 （如山区）的变形更为复杂，难以用简单的仿射模型纠正。如图3 （a）和图3 （b）所示，左图是山区的LRMS图像，中间是用ENVI配准后的图像与参考图像的差值图像，右图是本文算法配准后的图像与参考图像的差值图像。从图像可以看出，右图差值图像灰度较为平滑，无较多明显误匹配点，中间ENVI配准图像有较多灰度凸起，配准效果较差。表2列出上述两组图像的均方根误差和相关系数，由表2可知，本文方法在山区图像配准上精度优于ENVI手动配准，但是有些区域精度依然有待提高，这也是后续研究的重点。

表2 配准精度

3 结束语

遥感图像配准是遥感图像处理领域的一项基本研究内容，是诸多遥感图像应用的基础步骤。同一卫星的PAN 和LRMS存在复杂的非刚性位移，为了得到高精度的配准结果，本文提出一种基于特征约束与光流场模型的配准方法，并详细阐述将该方法运用到全色与多光谱图像配准过程中的具体实施步骤。本文实验证明，基于特征约束与光流场模型的配准方法精度可以达到亚像素级，目视效果较为理想。

光流场方法在遥感图像配准中的应用还有很多问题值得探讨与研究：首先该方法需要合理选择参数，这些参数目前没有较客观的计算方法，多数是通过实验的经验值设定；其次是计算速度，由于光流场方法对每一像素进行光流计算，计算量较大，对于实时性要求较高的系统，该方法还需在计算速度上有所改进；最后，可以进一步探讨该方法在多源遥感图像配准中的应用。

［1］Yetik I S，Nehorai A.Performance bounds on image registration ［J］.IEEE Transactions on Signal Processing，2006，54（5）：1737－1749.

［2］Lowe DG.Distinctive image features from scale－invariant keypoints［J］.International Journal of Computer Vision，2004，60 （2）：91－110.

［3］Kern J P，Pattichis M S.Robust multispectral image registration using mutual－information models［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（5）：1494－1505.

［4］Wong A，Clausi D A.ARRSI：Automatic registration of remote－sensing images ［J］.Geoscience and Remote Sensing，2007，45 （5）：1483－1493.

［5］Xu M，Varshney P K，Niu R.A subspace method for Fourier based image registration ［C］／／Signals，Systems and Computers，2006：425－429.

［6］LU Ziyun.The optical flow field calculation and its some optimization technologies ［D］.Hefei：Hefei University of Technology，2012 （in Chinese）.［路子赟.光流场计算及其若干优化技术研究 ［D］.合肥：合肥工业大学，2012.］

［7］Brox T，Bregler C，Malik J.Large displacement optical flow［C］／／Computer Vision and Pattern Recognition，2009：41－48.

［8］Bruhn A，Weickert J，Schnrr C.Lucas／Kanade meets Horn／Schunck：Combining local and global optic flow methods［J］.International Journal of Computer Vision，2005，61 （3）：211－231.

［9］B ro x T，Malik J.Large displacement optical flow：Descriptor matching in variational motion estimation ［J］.Pattern Analysis and Machine Intelligence，2011，33 （3）：500－513.

［10］DOU Wen.Comparison among remotely sensed image fusion method based on spectral response function ［J］.Spectroscopy and Spectral Analysis，2011，31 （3）：746－752 （in Chinese）.［窦闻.基于光谱响应函数的遥感图像融合对比研究 ［J］.光谱学与光谱分析，2011，31 （3）：746－752.］