CN103839265A - 基于sift和归一化互信息的sar图像配准方法 - Google Patents

Abstract

本发明提供了一种基于SIFT和归一化互信息的SAR图像配准方法，实现过程是(1)输入基准图像I₁和待配准图像I₂，分别进行预处理；(2)用MM-SIFT方法对预处理后的I₁和I₂进行特征提取，得到初始特征点对Fc和SIFT特征向量Fv₁、Fv₂；(3)用Fv₁、Fv₂进行初始匹配；(4)用单应矩阵模型的RANSAC策略对Fc进行二次筛选，得到最终正确的匹配点对Fm，采用最小二乘法求出配准参数pr；(5)用仿射变换对I₂进行空间变换，经插值和重采样，得到粗配准后的图像I₃；(6)使用pr作为归一化互信息配准的初始值，采用归一化互信息方法对I₁和I₂进行精配准，求得最终配准参数pr₁，并输出配准后的图像I₄。本发明快速、有效和稳定，提高了SAR图像配准的精度和鲁棒性。

Description

基于SIFT和归一化互信息的SAR图像配准方法

技术领域

本发明属于图像处理技术领域，具体涉及一种基于SIFT和归一化互信息的SAR图像配准方法，应用于图像融合、目标检测、地景与地图的匹配和图像导航等领域。

背景技术

合成孔径雷达（Synthetic Aperture Radar,SAR）是一种远距离、非接触的对地观测工具。作为主动式成像雷达，具有全天候、全天时工作的特点，并且能有效地识别伪装和穿透地表植被，得到高分辨率的SAR图像，其日益成为最具有代表性的对地观测手段之一。配准的方法大致可分为两类：基于灰度相关的图像配准方法和基于特征的图像配准方法。

基于灰度相关的图像配准方法虽然比较直观，容易实现，但容易受到噪声的影响，而且匹配利用的是图像的统计特性，计算量也很大。

基于特征点的图像配准方法对灰度变化有较强适应能力，且计算量小，能够处理图像间大失配的情况。常见的特征点检测算法有：Moravec特征点检测算法、SUSAN特征点检测算法、Harris特征点检测算法和SIFT特征点检测算法。四种算子中，Moravec算子抗噪性比较差，对旋转敏感；Harris算子是最不稳定的一种；SUSAN算子匹配性能较差，点匹配的应用上较少；SIFT特征点的信息丰富，相比其它三种算子，基于SIFT特征向量欧氏距离的匹配策略要更准确和更稳定，这正是SIFT特征点独有的优势。

1999年DavidG.Lowe总结了现有的基于不变量技术的特征检测方法，并正式提出了一种基于尺度空间的、对图像缩放、旋转甚至仿射变换保持不变性的图像局部特征描述算子－SIFT，这种算法在2004年被加以完善；2012年，文章《Modified SIFT for multi-modal remote sensing image registration》在原SIFT方法的基础上进行修改，提出了MM-SIFT方法，以改善遥感图像的配准效果。

SIFT算法在自然图像配准中效果比较好，而应用到SAR图像时，会出现很多问题。首先，受SAR图像相干斑噪声的影响，很可能将噪点检测为特征点，这样就会导致错误的配准结果；其次，机载SAR图像受实际条件的干扰，可能亮度、对比度比较暗，不能直接对SAR图像进行匹配，否则，得不到正确的匹配点对；最后，在初始匹配后，误匹配不可避免，这对配准结果影响很大。如何快速有效地剔除误匹配点对、提高配准的精度是非常重要的研究方向。

发明内容

本发明的目的是针对上述已有技术的不足，提出一种基于SIFT和归一化互信息的SAR图像配准方法，以提高SAR图像配准的精度和鲁棒性，实现SAR图像配准。

本发明的技术方案是提出一种基于SIFT和归一化互信息的SAR图像配准方法，其特征在于：包含以下步骤：

步骤1：输入两幅SAR图像，其中，一幅是基准图像I₁，另一幅是待配准图像I₂，对两幅SAR图像分别进行预处理，首先采用瑞利分布，增强系数为0.2对两幅SAR图像I₁、I₂进行限制对比度直方图均衡(Contrast Limited AdaptiveHistogram Equalization,CLAHE)增强，然后采用db1小波，分解层数为2～3进行小波分解；

步骤2：用MM-SIFT方法对两幅SAR图像分别进行特征提取，得到初始特征点对Fc和SIFT特征向量Fv₁、Fv₂；

步骤3：用SIFT特征向量Fv₁、Fv₂进行初始匹配；

步骤4：用单应矩阵模型的RANSAC策略对初始特征点对Fc进行二次筛选，得到最终正确的匹配点对Fm，然后采用最小二乘法求出配准参数pr；

步骤5：用仿射变换对待配准图像I₂进行空间变换，然后经过插值和重采样，得到粗配准后的图像I₃；

步骤6：使用上述配准参数pr作为归一化互信息配准的初始值，采用归一化互信息方法对两幅SAR图像I₁和I₂进行精配准，求得最终配准参数pr₁，并且输出配准后的图像I₄。

上述步骤2中所述的MM-SIFT方法进行特征提取的具过程包括如下：

2.1高斯模糊和尺度空间生成

对于二维图像I(x,y)，在不同尺度下的尺度空间L(x,y,σ)可由图像I(x,y)与高斯核G(x,y,σ)的卷积得到：

L(x,y,σ)＝G(x,y,σ)*I(x,y)

其中，

(x,y)代表I上的点，σ是尺度因子；高斯金子塔的构建过程可分为两步：(1)对图像做高斯平滑；(2)对图像做降采样；

为了让尺度体现其连续性，在简单下采样的基础上加上了高斯滤波，一幅图像可以产生几组(octave)图像，一组图像包括几层(interval)图像，高斯图像金字塔共O组、S层，则有：

σ(o,s)＝σ₀2^o+s/S，o∈o_min+[0,...,O-1]，s∈[-1,....,S+1]

σ₀-初始尺度；S-每组层数；s-组内某一层数；o-某一组；

高斯金字塔初始尺度σ₀，当图像通过相机拍摄时，相机的镜头已经对图像进行了一次初始的模糊，所以根据高斯模糊的性质：

${\mathrm{\σ}}_0=\sqrt{{\mathrm{\σ}}_{\mathrm{init}}\×{\mathrm{\σ}}_{\mathrm{init}}-{\mathrm{\σ}}_{\mathrm{pre}}\×{\mathrm{\σ}}_{\mathrm{pre}}}$

σ_init-第0层尺度；σ_pre-被相机镜头模糊后的尺度；

高斯金字塔的组数：

O＝[log₂(min(M,N))]-3

其中，M、N分别为图像的行数和列数；

DOG(Difference of Gaussian)即相邻两尺度空间函数之差，用D(x,y,σ)来表示，即：

D(x,y,σ)＝(G(x,y,kσ)-G(x,y,σ))*I(x,y)＝L(x,y,kσ)-L(x,y,σ)

DOG金字塔的第1层的尺度因子与高斯金字塔的第1层的是一致的，其它阶也是一样；

2.2极值点检测和关键点精确定位

为了能够检测到DOG空间的极大值和极小值，DOG尺度空间除去最底层和最顶层的中间层中每个像素点都要与同一层相邻的8个像素点以及它上下两层邻近的9个像素点总共26个相邻像素点进行逐个比较，以确保在二维图像空间和尺度空间都能检测到局部极值；

对尺度空间DOG函数进行曲线拟合，利用DOG所在尺度空间的Taylor展开式：

$D\left(X\right)=D+\frac{{\∂D}^T}{\∂X}X+\frac{1}{2}{X}^T\frac{{\∂}^{2}D}{{\∂X}^2}X$

求导并让方程等于零，可以得到极值点：

$\hat{X}=\frac{{\∂D}^T}{\∂X}{\left(\frac{{\∂}^{2}D}{{\∂X}^2}\right)}^{-1}$

为修正值，对应极值点，方程的值：

$D\left(\hat{X}\right)=D+\frac{1}{2}\frac{{\∂D}^T}{\∂X}\hat{X}$

将

的极值点视为低对比度的不稳定的极值点，进行剔除；

DOG函数的主曲率通过计算在该点位置和尺度的2×2的Hessian矩阵H得到，导数可由相邻采样点的差来估计：

$H=\left[\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{xy}}& D_{\mathrm{yy}}\end{array}\right]$

令α为较大的特征值，β为小的特征值，α＝rβ,则：

Tr(H)＝D_xx+D_yy＝α+βDet(H)＝D_xxD_yy-(D_xy)²＝αβ

$\mathrm{ratio}=\frac{\mathrm{Tr}{\left(H\right)}^2}{\mathrm{Det}\left(H\right)}=\frac{{(\mathrm{\α}+\mathrm{\β})}^2}{\mathrm{\α\β}}=\frac{{(r+1)}^2}{r}$

Lowe论证建议r=10，若

时将关键点保留，反之剔除；

2.3关键点方向分配

使用关键点邻域像素梯度方向的分布，可以指定每个关键点的方向参数，从而使描述子对图像旋转具有不变性；

像素点的梯度表示：

$\mathrm{grad}I(x,y){=}_4(\frac{\∂I}{\∂x},\frac{\∂I}{\∂y})$

梯度幅值：

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(x,y+1)-L(x,y-1))}^2}$

梯度方向：

$\mathrm{\θ}(x,y)={\tan }^{-1}\left[\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\right]$

关键点的方向分配实现步骤：(1)计算关键点邻域直方图的高斯权重参数；(2)生成含有36个柱的梯度直方图，其中范围为0～360度，每10度一个柱，由半径为图像区域生成，柱所代表的方向为像素点的梯度方向，柱的长短表示梯度的幅值。然后对梯度直方图做两次平滑滤波；(3)求取关键点方向；(4)二次拟合方向直方图的Taylor展开式，精确关键点的方向；

关键点主方向：梯度方向直方图的主峰值，也称为是特征点方向；关键点的辅方向：当存在其它的相当于梯度直方图中主峰值80%能量的峰值；

2.4关键点描述子的生成

128维关键点的特征向量生成步骤如下：

(1)确定计算关键点描述子所需要的图像块的大小，图像区域表示关键点所在的尺度空间，其半径通过下式计算：

$\mathrm{radius}=\frac{3{\mathrm{\σ}}_{\mathrm{oct}}\×\sqrt{2}\×(d+1)+1}{2}$

其中σ_oct为关键点的组内尺度，d＝4；

(2)将坐标移到关键点的主方向，旋转后的新坐标表示为：

$\left(\begin{array}{c}\hat{x}\\ \hat{y}\end{array}\right)=\left(\begin{array}{cc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right)\×\left(\begin{array}{c}x\\ y\end{array}\right)$

(3)在图像区域内计算每个像素点的梯度方向和幅值，然后对每个求得的梯度幅值乘以高斯权重参数，生成方向直方图：

$\mathrm{weight}=\left|\mathrm{grad}\left(I_{\mathrm{\σ}}(x,y)\right)\right|\×\exp (-\frac{x_k^2+y_k^2}{2{\mathrm{\σ}}_w})\×(1-d_r)\×(1-d_c)\×(1-d_o)$

x_k：表示该像素点和中心点的列距离；y_k：表示该像素点和中心点的行距

离；σ_w：表示描述子窗口的宽度3σ×直方图列数（取4）的一半；

(4)在每个窗口宽度为2×2的区域内求得8个方向的梯度直方图，计算每个梯度方向的累加值，形成一个种子点。共生成16个种子点，每一个种子点得到一个8位的特征向量，这样就生成4×4×8共128维的特征向量；

(5)描述子向量元素门限化：方向直方图每个方向上梯度幅值限制在一定门限值以下，一般门限值取0.2，描述子向量元素规范化：

$\begin{array}{cc}{l}_j=w_j/\sqrt{\underset{i=1}{\overset{128}{\mathrm{\Σ}}}{w}_i}& j=\mathrm{1,2},...,128\end{array}$

其中，W＝(w₁,w₂,...,w₁₂₈)为得到的128维特征向量，L＝(l₁,l₂,...,l₁₂₈)为规范化的向量，关键点描述子向量的规范化正是可去除满足此模型的光照影响，图像各点的梯度是通过邻域像素相减得到，对于图像灰度值整体漂移，也能去除；

设置窗口大小为标准SIFT窗口的1.67倍，把4×4的增大到6×6，产生6×6×8共288维向量描述子，增加正确匹配的点数，减弱噪声的影响，

上述步骤3所述的SIFT特征向量匹配，通过如下公式进行：

参考图中关键点描述子：R_i＝(r_i1,r_i2,...,r_i128)，实时图中关键点描述子：S_i＝(s_i1,s_i2,...,s_i128)，任意两描述子相似性度量：

$d(R_i,S_i)=\sqrt{\underset{j=1}{\overset{128}{\mathrm{\Σ}}}{(r_{\mathrm{ij}}-s_{\mathrm{ij}})}^2}$

要得到配对的关键点描述子，d(R_i,S_j)需满足：

上述步骤4中所述的单应矩阵模型的RANSAC策略，通过如下公式进行：

4.1给定两幅图像I1，I2的初始匹配点集X1((x_i,y_i),i＝1,2,...,n)，X2((x_j,y_j),j＝1,2,...,n)，将它们分别写成齐次坐标的形式：X1((x_i,y_i,1),i＝1,2,...,n)，X2((x_j,y_j,1),j＝1,2,...,n)；

4.2从匹配点对集合中，任选四个点对，计算单应矩阵H；

4.3在所有匹配点对集合中，以欧氏距离为判据，选取所有符合d||HX1-X2||＜t的点对，将它们作为内点，不符合此条件的点对认为是误匹配点对，予以剔除，并记录满足此约束的内点数量，实验采用阈值t＝36；

4.4重复循环4.2、4.3至N次，记录每一次的内点数量，实验采用N＝2000；

4.5选取对应内点数最大的H，并把此次内点作为最终的内点，即正确匹配点对。

上述步骤4中所述的最小二乘法求解配准参数，通过如下公式进行：

$\left[\begin{array}{cccccc}x& y& 0& 0& 1& 0\\ 0& 0& x& y& 0& 1\\ & & ...& & & \\ & & ...& & & \end{array}\right]\left[\begin{array}{c}{a}_{11}\\ a_{12}\\ a_{21}\\ a_{22}\\ t_x\\ t_y\end{array}\right]=\left[\begin{array}{c}u\\ v\\ .\\ .\\ .\end{array}\right]$

上式可以表示为：Ax＝b，用最小二乘法求得：x＝[A^TA]^-1ATb，进而可以求出，旋转角度：angle＝-(180/π)atan(a₂₁/a₁₁)，水平平移量：Δx＝t_x，垂直平移量：Δy＝t_y，尺度变化：

上述步骤6中所述的基于归一化互信息的SAR图像配准方法，通过如下公式进行：

6.1归一化互信息配准参数的初始值修正为angle、Δx和Δy；

6.2利用Powell算法依据归一化互信息判别所求得的参数是否最优，若不是，则继续优化更新几何变换，直至获得最大归一化互信息值，搜索到满足精度要求的参数，达到实时图的最优变换，求得配准参数和输出配准后的图像；

Powell算法步骤如下：

(1)给定允差ε＞0，初始点X⁽⁰⁾和n个d^(1,1),d^(1,2),...,d^(1,n)线性无关的方向，置k＝1；

其中，ε为允许误差，X⁽⁰⁾为初始点，d^(1,1),d^(1,2),...,d^(1,n)为n个线性无关的方向，k为循环次数；

(2)置x^(k,0)＝x^(k-1)，从x^(k,0)出发，依次沿方向d^(k,1),d^(k,2),...,d^(k,n)做一次一维搜索，得到点x^(k,1),x^(k,2),...,x^(k,n)，求m，使得

$f\left(x^{(k,m-1)}\right)-f\left(x^{(k,m)}\right)=\underset{j=1,...,n}{\max }\{f\left(x^{(k,j-1)}\right)-f\left(x^{(k,j)}\right)\}$

令d^(k,n+1)＝x^(k,n)-x^(k,0)，若||x^(k,n)-x^(k,0)||≤ε，则停止计算；否则，进行步骤(3)；

其中，x^(k,0)为每次搜索的初始位置，d^(k,1),d^(k,2),...,d^(k,n)为n个线性无关的搜索方向，x^(k,1),x^(k,2),...,x^(k,n)为得到的参数值，f(x)为互信息的值，d^(k,n+1)为新的搜索方向；

(3)求λ_n+1，使得 $f(x^{(k,0)}+{\mathrm{\λ}}_{n+1}{d}^{(k,n+1)})=\underset{\mathrm{\λ}}{\min }f(x^{(k,0)}+\mathrm{\λ}{d}^{(k,n+1)}),$ 令x^(k+1,0)＝x^(k)＝x^(k,0)+λ_n+1d^(k,n+1)，若||x^(k)-x^(k-1)||≤ε，则停止计算，得到x^(k)；否则，进行步骤（4）；

其中，λ_n+1，λ为搜索方向的斜率，x^(k)为得到的参数值；

(4)若

则令d^(k+1,j)＝d^(k,j)，j＝1,...,m-1，d^(k+1,j)＝d^(k,j+1)，j＝m,....,n，k＝k+1；否则，令d^(k+1,j)＝d^(k,j)，j＝1,...,n，k＝k+1，转到步骤(2)；

其中，d^(k+1,j)为得到的新的搜索方向，k为循环次数，n为方向个数。

本发明的有益效果是：

第一，本发明采用图像预处理，包括CLAHE增强和小波分解；通过实验验证，如果图像在800×800以上，采用小波分解不但满足一定的去噪效果，而且可以大大缩短匹配的时间；如果图像在800×800以下，时间区别不大，则不采用小波分解，这样能够保证后续配准精度。

第二，本发明对SIFT特征提取，关键点描述子进行修改；相比标准SIFT，采用保留边缘点，对特征向量窗口大小和特征点描述子维数进行了修改，这样能够增加正确的匹配点对，减弱噪声的影响，有利于提高配准精度；采用单应矩阵约束的RANSAC策略能够有效地去除误匹配点对，而且比较稳定。

第三，本发明对归一化互信息配准的初始值进行修正；利用改进SIFT配准得到的参数作为归一化互信息精配准的初始值，然后利用改进Powell算法不断优化更新几何变换，直至搜索到满足精度要求的参数，获得最大归一化互信息值，输出配准后的图像，相比直接用互信息搜索，速度有很大提高，并且配准精度能够精确到亚像素级。

以下将结合附图对本发明做进一步详细说明。

附图说明

图1是本发明的流程图；

图2(a)是大小为400×400的待配准SAR图像；

图2(b)是大小为400×400的基准SAR图像；

图2(c)是基于改进SIFT的SAR图像配准方法配准的结果图；

图2(d)是基于SIFT的自动图像配准方法配准的结果图；

图2(e)是基于本发明的图像配准方法配准的结果图；

图3(a)是大小为600×400的待配准SAR图像；

图3(b)是大小为600×400的基准SAR图像；

图3(c)是基于改进SIFT的SAR图像配准方法配准的结果图；

图3(d)是基于SIFT的自动图像配准方法配准的结果图；

图3(e)是基于本发明的图像配准方法配准的结果图；

图4(a)是大小为700×700的待配准SAR图像；

图4(b)是大小为700×700的基准SAR图像；

图4(c)是基于改进SIFT的SAR图像配准方法配准的结果图；

图4(d)是基于SIFT的自动图像配准方法配准的结果图；

图4(e)是基于本发明的图像配准方法配准的结果图；

图5(a)是大小为400×400的待配准SAR图像；

图5(b)是大小为400×400的基准SAR图像；

图5(c)是基于改进SIFT的SAR图像配准方法配准的结果图；

图5(d)是基于SIFT的自动图像配准方法配准的结果图；

图5(e)是基于本发明的图像配准方法配准的结果图；

具体实施方式

下面结合附图对本发明进一步详细描述。

实施例1：

本发明提出了一种基于SIFT和归一化互信息的SAR图像配准方法，本例的仿真是在主频2.60GHz的CPU intel Pentium Dual-Core E5300、内存2GB的硬件环境和MATLAB R2011a的软件环境的Windows XP SP3系统下进行的。

本发明提出一种基于SIFT和归一化互信息的SAR图像配准方法，以提高SAR图像配准的精度和鲁棒性，实现SAR图像配准，如图1所示，本发明实现SAR图像配准过程包含以下步骤：

步骤1：输入两幅SAR图像，其中，一幅是基准图像I₁，另一幅是待配准图像I₂，对两幅SAR图像分别进行预处理，首先采用瑞利分布，增强系数为0.2对两幅SAR图像I₁、I₂进行限制对比度直方图均衡(Contrast Limited AdaptiveHistogram Equalization,CLAHE)增强，然后采用db1小波，分解层数为2～3进行小波分解；

步骤2：用MM-SIFT方法对两幅SAR图像分别进行特征提取，得到初始特征点对Fc和SIFT特征向量Fv₁、Fv₂；

步骤3：用SIFT特征向量Fv₁、Fv₂进行初始匹配；

步骤4：用单应矩阵模型的RANSAC策略对初始特征点对Fc进行二次筛选，得到最终正确的匹配点对Fm，然后采用最小二乘法求出配准参数pr；

步骤5：用仿射变换对待配准图像I₂进行空间变换，然后经过插值和重采样，得到粗配准后的图像I₃；

步骤6：使用上述配准参数pr作为归一化互信息配准的初始值，采用归一化互信息方法对两幅SAR图像I₁和I₂进行精配准，求得最终配准参数pr₁，并且输出配准后的图像I₄。

实施例2：

基于SIFT和归一化互信息的SAR图像配准方法同实施例1，为了具备可实施性，对本发明进一步的详细描述如下：

步骤2中所述的MM-SIFT方法进行特征提取的具过程包括如下：

2.1高斯模糊和尺度空间生成

对于二维图像I(x,y)，在不同尺度下的尺度空间L(x,y,σ)可由图像I(x,y)与高斯核G(x,y,σ)的卷积得到：

L(x,y,σ)＝G(x,y,σ)*I(x,y)

其中，

(x,y)代表I上的点，σ是尺度因子。高斯金子塔的构建过程可分为两步：(1)对图像做高斯平滑；(2)对图像做降采样；

为了让尺度体现其连续性，在简单下采样的基础上加上了高斯滤波，一幅图像可以产生几组(octave)图像，一组图像包括几层(interval)图像，高斯图像金字塔共O组、S层，则有：

σ(o,s)＝σ₀2^o+s/S，o∈o_min+[0,...,O-1]，s∈[-1,....,S+1]

σ₀-初始尺度；S-每组层数；s-组内某一层数；o-某一组；

高斯金字塔初始尺度σ₀，当图像通过相机拍摄时，相机的镜头已经对图像进行了一次初始的模糊，所以根据高斯模糊的性质：

${\mathrm{\σ}}_0=\sqrt{{\mathrm{\σ}}_{\mathrm{init}}\×{\mathrm{\σ}}_{\mathrm{init}}-{\mathrm{\σ}}_{\mathrm{pre}}\×{\mathrm{\σ}}_{\mathrm{pre}}}$

σ_init-第0层尺度；σ_pre-被相机镜头模糊后的尺度；

高斯金字塔的组数：

O＝[log₂(min(M,N))]-3

其中，M、N分别为图像的行数和列数；

DOG(Difference of Gaussian)即相邻两尺度空间函数之差，用D(x,y,σ)来表示，即：

D(x,y,σ)＝(G(x,y,kσ)-G(x,y,σ))*I(x,y)＝L(x,y,kσ)-L(x,y,σ)

DOG金字塔的第1层的尺度因子与高斯金字塔的第1层的是一致的，其它阶也是一样；

2.2极值点检测和关键点精确定位

为了能够检测到DOG空间的极大值和极小值，DOG尺度空间中间层（除去最底层和最顶层）的每个像素点都要与同一层相邻的8个像素点以及它上下两层邻近的9个像素点总共26个相邻像素点进行逐个比较，以确保在二维图像空间和尺度空间都能检测到局部极值；

为了提高关键点的稳定性，同时去除低对比度的关键点和不稳定的边缘响应点(因为DOG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力，需要对尺度空间DOG函数进行曲线拟合；

利用DOG所在尺度空间的Taylor展开式：

$D\left(X\right)=D+\frac{{\∂D}^T}{\∂X}X+\frac{1}{2}{X}^T\frac{{\∂}^{2}D}{{\∂X}^2}X$

求导并让方程等于零，可以得到极值点：

$\hat{X}=\frac{{\∂D}^T}{\∂X}{\left(\frac{{\∂}^{2}D}{{\∂X}^2}\right)}^{-1}$

为修正值，对应极值点，方程的值：

$D\left(\hat{X}\right)=D+\frac{1}{2}\frac{{\∂D}^T}{\∂X}\hat{X}$

的值对于剔除低对比度的不稳定的极值点十分有用，通常将

的极值点视为低对比度的不稳定的极值点，进行剔除；

仅去除低对比度的极值点对于特征点的稳定是远远不够的，DOG函数在图像边缘有较强的边缘响应，因此我们还需要去除边缘响应；

在DOG函数（欠佳的）的峰值点中，较大的主曲率存在于横跨边缘的方向，而在垂直边缘的方向主曲率较小，主曲率通过计算在该点位置和尺度的2×2的Hessian矩阵H得到，导数可由相邻采样点的差来估计：

$H=\left[\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{xy}}& D_{\mathrm{yy}}\end{array}\right]$

为了避免直接计算这些特征值，而只是考虑它们的之间的比率，令α为较大的特征值，β为小的特征值，α＝rβ,则：

Tr(H)＝D_xx+D_yy＝α+βDet(H)＝D_xxD_yy-(D_xy)²＝αβ

$\mathrm{ratio}=\frac{\mathrm{Tr}{\left(H\right)}^2}{\mathrm{Det}\left(H\right)}=\frac{{(\mathrm{\α}+\mathrm{\β})}^2}{\mathrm{\α\β}}=\frac{{(r+1)}^2}{r}$

(r+1)²/r在两特征值相同时最小，随r的增大，比率值也在逐渐增大。r越大，即在某一方向的梯度值越大，而在另一方向的梯度值越小，而边缘恰好就是这种情况，所以为了去除边缘相应，要让该比值小于一定的阈值；

Lowe论证建议r=10，若

时将关键点保留，反之剔除；

在SIFT特征提取中，低对比度和边缘上的点被去除，由于对SAR图像提取的关键点比较少，我们应该保留低对比度和边缘上的点，这样能够产生更多正确的特征点对，有利于提高配准精度；

2.3关键点方向分配

使用关键点邻域像素梯度方向的分布，可以指定每个关键点的方向参数，从而使描述子对图像旋转具有不变性；

像素点的梯度表示：

$\mathrm{grad}I(x,y){=}_4(\frac{\∂I}{\∂x},\frac{\∂I}{\∂y})$

梯度幅值：

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(x,y+1)-L(x,y-1))}^2}$

梯度方向：

$\mathrm{\θ}(x,y)={\tan }^{-1}\left[\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\right]$

确定关键点的方向采用梯度直方图统计法，统计以关键点为原点，一定区域内的图像像素点对关键点方向生成所作的贡献，方向分配实现步骤：(1)计算关键点邻域直方图的高斯权重参数；(2)生成含有36个柱的梯度直方图，其中范围为0～360度，每10度一个柱，由半径为图像区域生成，柱所代表的方向为像素点的梯度方向，柱的长短表示梯度的幅值。然后对梯度直方图做两次平滑滤波；(3)求取关键点方向（可能是多个方向）；(4)二次拟合方向直方图的Taylor展开式，精确关键点的方向；

关键点主方向：梯度方向直方图的主峰值，也称为是特征点方向；关键点的辅方向：当存在其它的相当于梯度直方图中主峰值80%能量的峰值；

图像中的关键点检测完毕后，每个关键点包含三个信息：位置、尺度和方向，这就使关键点具备平移、旋转和缩放不变性；

2.4关键点描述子的生成

128维关键点的特征向量生成步骤如下：

(1)确定计算关键点描述子所需要的图像块的大小，图像区域表示关键点所在的尺度空间，其半径通过下式计算：

$\mathrm{radius}=\frac{3{\mathrm{\σ}}_{\mathrm{oct}}\×\sqrt{2}\×(d+1)+1}{2}$

其中σ_oct为关键点的组内尺度，d＝4；

(2)将坐标移到关键点的主方向，旋转后的新坐标表示为：

$\left(\begin{array}{c}\hat{x}\\ \hat{y}\end{array}\right)=\left(\begin{array}{cc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right)\×\left(\begin{array}{c}x\\ y\end{array}\right)$

(3)在图像区域内计算每个像素点的梯度方向和幅值，然后对每个求得的梯度幅值乘以高斯权重参数，生成方向直方图：

$\mathrm{weight}=\left|\mathrm{grad}\left(I_{\mathrm{\σ}}(x,y)\right)\right|\×\exp (-\frac{x_k^2+y_k^2}{2{\mathrm{\σ}}_w})\×(1-d_r)\×(1-d_c)\×(1-d_o)$

x_k：表示该像素点和中心点的列距离；y_k：表示该像素点和中心点的行距离；σ_w：表示描述子窗口的宽度3σ×直方图列数（取4）的一半；

(4)在每个窗口宽度为2×2的区域内求得8个方向的梯度直方图，计算每个梯度方向的累加值，形成一个种子点。共生成16个种子点，每一个种子点得到一个8位的特征向量，这样就生成4×4×8共128维的特征向量；

(5)描述子向量元素门限化：方向直方图每个方向上梯度幅值限制在一定门限值以下（门限一般取0.2），描述子向量元素规范化：

$\begin{array}{cc}{l}_j=w_j/\sqrt{\underset{i=1}{\overset{128}{\mathrm{\Σ}}}{w}_i}& j=\mathrm{1,2},...,128\end{array}$

其中，W＝(w₁,w₂,...,w₁₂₈)为得到的128维特征向量，L＝(l₁,l₂,...,l₁₂₈)为规范化的向量，关键点描述子向量的规范化正是可去除满足此模型的光照影响，图像各点的梯度是通过邻域像素相减得到，对于图像灰度值整体漂移，也能去除；

由于待配准的SAR图像有平移、旋转、拍摄时雷达视角不同，配准的图像包含很多相似的区域，这使得用特征向量匹配存在困难。为了减少误匹配的影响，必须增大特征点描述子窗口的大小，通过实验验证，窗口大小是标准SIFT窗口的1.67倍时，效果最好；

由于窗口大小增大1.67倍，在计算描述子时，也要相应增加特征向量的维数，否则，每一个2×2子区域不能提取足够多的特征，为了解决这个问题，我们把4×4的增大到6×6，这将会产生6×6×8共288维向量描述子，这种方法会增加正确匹配的点数，减弱噪声的影响，如果再增大特征向量的维数，这样会导致不稳定的梯度直方图，对图像形变比较敏感，增加匹配的时间。

步骤3所述的SIFT特征向量匹配，通过如下公式进行：

参考图中关键点描述子：R_i＝(r_i1,r_i2,...,r_i128)，实时图中关键点描述子：S_i＝(s_i1,s_i2,...,s_i128)，任意两描述子相似性度量：

$d(R_i,S_i)=\sqrt{\underset{j=1}{\overset{128}{\mathrm{\Σ}}}{(r_{\mathrm{ij}}-s_{\mathrm{ij}})}^2}$

要得到配对的关键点描述子，d(R_i,S_j)需满足：

当distRatio取0.75～0.85时，可以兼顾匹配特征点数目和匹配正确率，匹配效果最佳。

步骤4中所述的单应矩阵模型的RANSAC策略，通过如下公式进行：

4.1给定两幅图像I1，I2的初始匹配点集X1((x_i,y_i),i＝1,2,...,n)，X2((x_j,y_j),j＝1,2,...,n)，将它们分别写成齐次坐标的形式：X1((x_i,y_i,1),i＝1,2,...,n)，X2((x_j,y_j,1),j＝1,2,...,n)；

4.2从匹配点对集合中，任选四个点对，计算单应矩阵H；

4.3在所有匹配点对集合中，以欧氏距离为判据，选取所有符合d||HX1-X2||＜t的点对，将它们作为内点，不符合此条件的点对认为是误匹配点对，予以剔除，并记录满足此约束的内点数量，实验采用阈值t＝36；

4.4重复循环4.2、4.3至N次，记录每一次的内点数量，实验采用N＝2000；

4.5选取对应内点数最大的H，并把此次内点作为最终的内点，即正确匹配点对。

步骤4中所述的最小二乘法求解配准参数，具体到SAR图像配准中，变换公式如下：

$\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{t}_x\\ t_y\end{array}\right]$

P=[a₁₁ a₁₂ a₂₁ a₂₂ t_x t_y]为需要求解的仿射变换模型参数。其中，[t_x t_y]^T表示水平和垂直平移量，a_ii表示尺度和旋转量，[x,y]为待配准图像中特征点的坐标，[u,v]为参考图像中对应的特征点的坐标，一般至少需要三对匹配点来求解空间变换模型参数，仿射模型参数通过如下公式进行：

$\left[\begin{array}{cccccc}x& y& 0& 0& 1& 0\\ 0& 0& x& y& 0& 1\\ & & ...& & & \\ & & ...& & & \end{array}\right]\left[\begin{array}{c}{a}_{11}\\ a_{12}\\ a_{21}\\ a_{22}\\ t_x\\ t_y\end{array}\right]=\left[\begin{array}{c}u\\ v\\ .\\ .\\ .\end{array}\right]$

上式可以表示为：Ax＝b，用最小二乘法求得：x＝[A^TA]^-1A^Tb，进而可以求出，旋转角度：angle＝-(180/π)atan(a₂₁/a₁₁)，水平平移量：Δx＝t_x，垂直平移量：Δy＝t_y，尺度变化：

步骤6中所述的基于归一化互信息的SAR图像配准方法，通过如下公式进行：

6.1归一化互信息配准参数的初始值修正为angle、Δx和Δy；

6.2利用Powell算法依据归一化互信息判别所求得的参数是否最优，若不是，则继续优化更新几何变换，直至获得最大归一化互信息值，搜索到满足精度要求的参数，达到实时图的最优变换，求得配准参数和输出配准后的图像。

Powell算法步骤如下：

(2)给定允差ε＞0，初始点X⁽⁰⁾和n个d^(1,1),d^(1,2),...,d^(1,n)线性无关的方向，置k＝1；

其中，ε为允许误差，X⁽⁰⁾为初始点，d^(1,1),d^(1,2),...,d^(1,n)为n个线性无关的方向，k为循环次数；

(2)置x^(k,0)＝x^(k-1)，从x^(k,0)出发，依次沿方向d^(k,1),d^(k,2),...,d^(k,n)做一次一维搜索，得到点x^(k,1),x^(k,2),...,x^(k,n)，求m，使得

$f\left(x^{(k,m-1)}\right)-f\left(x^{(k,m)}\right)=\underset{j=1,...,n}{\max }\{f\left(x^{(k,j-1)}\right)-f\left(x^{(k,j)}\right)\}$

令d^(k,n+1)＝x^(k,n)-x^(k,0)，若||x^(k,n)-x^(k,0)||≤ε，则停止计算；否则，进行步骤(3)；

其中，x^(k,0)为每次搜索的初始位置，d^(k,1),d^(k,2),...,d^(k,n)为n个线性无关的搜索方向，x^(k,1),x^(k,2),...,x^(k,n)为得到的参数值，f(x)为互信息的值，d^(k,n+1)为新的搜索方向；

(3)求λ_n+1，使得 $f(x^{(k,0)}+{\mathrm{\λ}}_{n+1}{d}^{(k,n+1)})=\underset{\mathrm{\λ}}{\min }f(x^{(k,0)}+\mathrm{\λ}{d}^{(k,n+1)}),$ 令x^(k+1,0)＝x^(k)＝x^(k,0)+λ_n+1d^(k,n+1)，若||x^(k)-x^(k-1)||≤ε，则停止计算，得到x^(k)；否则，进行步骤（4）；

其中，λ_n+1，λ为搜索方向的斜率，x^(k)为得到的参数值；

(4)若则令d^(k+1,j)＝d^(k,j)，j＝1,...,m-1，d^(k+1,j)＝d^(k,j+1)，j＝m,....,n，k＝k+1；否则，令d^(k+1,j)＝d^(k,j)，j＝1,...,n，k＝k+1，转到步骤(2)；

其中，d^(k+1,j)为得到的新的搜索方向，k为循环次数，n为方向个数。

实施例3：

基于SIFT和归一化互信息的SAR图像配准方法同实施例1-2，本发明的SAR图像配准效果可以通过以下实验进一步说明：

仿真实验环境为：MATLAB R2011a，CPU intel Pentium Dual-Core E53002.60GHz，内存2G，Windows XP SP3。

实验用基于改进SIFT的SAR图像配准方法和基于SIFT的自动图像配准方法分别对SAR图像进行匹配，并用本发明和上述两种配准方法进行比较，其中基于改进SIFT的SAR图像配准方法称为方法1；基于SIFT的自动图像配准方法称为方法2。图2(a)(b)，图3(a)(b)，图4(a)(b)和图5(a)(b)大小分别为400×400、600×400、700×700和400×400，是对黄河入海口的两幅图像分别截取的，黄河入海口的两幅图像分别于2008年和2009年均通过Radarsat-2卫星获取。

如图2所示，其中图2(a)为待配准图像，图2(b)为基准图像，图2(c)为方法1配准的结果，图2(d)为方法2配准的结果，图2(e)为本发明配准的结果；如图3所示，其中图3(a)为待配准图像，图3(b)为基准图像，图3(c)为方法1配准的结果，图3(d)为方法2配准的结果，图3(e)为本发明配准的结果；如图4所示，其中图4(a)为待配准图像，图4(b)为基准图像，图4(c)为方法1配准的结果，图4(d)为方法2配准的结果，图4(e)为本发明配准的结果；如图5所示，其中图5(a)为待配准图像，图5(b)为基准图像，图5(c)为方法1配准的结果，图5(d)为方法2配准的结果，图5(e)为本发明配准的结果。

从图2(c)、图3(c)、图4(c)和图5(c)均可以看出，对于SAR1、SAR2、SAR3和SAR4图像对，基于改进SIFT的SAR图像配准方法在图像块之间明显有错接，旋转角度、平移量误差都比较大。这主要因为在利用尺度和旋转准则二次筛选后，还有部分误匹配点对，造成配准结果误差较大。

从图2(d)、图3(d)、图4(d)和图5(d)均可以看出，对于SAR1、SAR2、SAR3和SAR4图像对，基于SIFT的自动图像配准方法在图像块之间也有偏移，相对于基于改进SIFT的SAR图像配准方法，旋转角度，平移量的误差较小。

从图2(e)、图3(e)、图4(e)和图5(e)均可以看出，对于SAR1、SAR2、SAR3和SAR4图像对，本发明得到的棋盘图在图像块之间能更好地衔接起来，效果比方法1和方法2都好。这源于改进SIFT能提取更多的特征点对，而且单应约束的RANSAC策略速度快，基本上去除所有的误匹配点对，效果好又稳定。然后利用改进SIFT粗配准得到的参数作为归一化互信息配准的初始值，配准精度和速度都有很大提高。

采用客观评价指标Δx,Δy、Δθ、RMSE和匹配时间，对本发明与上述两种方法的匹配性能进行比较，其结果如表1和表2所示，其中，理论值为人工选取的特征点对配准的结果，每次选取30对点，对5次结果求平均值。

表1配准结果比较

表2配准结果比较

从表1和表2可以看出，与现有的两种配准方法相比，本方法旋转角度误差不超过0.3度，平移误差不超过0.5个像素，RMSE误差不超过0.9个像素，配准精度达到亚像素级。

综上所述，本发明快速、有效和稳定，SAR图像配准精度能够精确到亚像素级水平，而且对于不同波段，不同视角，不同分辨率的SAR图像，本发明也具有良好的配准效果，为实际中的各种应用提供技术支持和实践价值。

Claims (6)

1.基于SIFT和归一化互信息的SAR图像配准方法，其特征在于：包含以下步骤： 

步骤1：输入两幅SAR图像，其中，一幅是基准图像I₁，另一幅是待配准图像I₂，对两幅SAR图像分别进行预处理，首先采用瑞利分布，增强系数为0.2对两幅SAR图像I₁、I₂进行限制对比度直方图均衡增强，然后采用db1小波，分解层数为2～3进行小波分解； 

步骤2：用MM-SIFT方法对预处理后的两幅SAR图像分别进行特征提取，得到初始特征点对Fc和SIFT特征向量Fv₁、Fv₂； 

步骤3：用SIFT特征向量Fv₁、Fv₂进行初始匹配； 

步骤4：用单应矩阵模型的RANSAC策略对初始特征点对Fc进行二次筛选，得到最终正确的匹配点对Fm，然后采用最小二乘法求出配准参数pr； 

步骤5：用仿射变换对待配准图像I₂进行空间变换，然后经过插值和重采样，得到粗配准后的图像I₃； 

步骤6：使用上述配准参数pr作为归一化互信息配准的初始值，采用归一化互信息方法对两幅SAR图像I₁和I₂进行精配准，求得最终配准参数pr₁，并且输出配准后的图像I₄。 

2.如权利要求1所述的基于SIFT和归一化互信息的SAR图像配准方法，其特征在于：步骤2中所述的MM-SIFT方法进行特征提取的具过程包括如下： 

2.1高斯模糊和尺度空间生成 

对于二维图像I(x,y)，在不同尺度下的尺度空间L(x,y,σ)可由图像I(x,y)与高斯核G(x,y,σ)的卷积得到： 

L(x,y,σ)＝G(x,y,σ)*I(x,y) 

其中，(x,y)代表I上的点，σ是尺度因子；高斯金子塔的构建过程可分为两步：(1)对图像做高斯平滑；(2)对图像做降采样； 

为了让尺度体现其连续性，在简单下采样的基础上加上了高斯滤波，一幅图像可以产生几组(octave)图像，一组图像包括几层(interval)图像，高斯图像金字塔共O组、S层，则有： 

σ(o,s)＝σ₀2^o+s/S，o∈o_min+[0,...,O-1]，s∈[-1,....,S+1] 

σ₀-初始尺度；S-每组层数；s-组内某一层数；o-某一组； 

高斯金字塔初始尺度σ₀，当图像通过相机拍摄时，相机的镜头已经对图像进行了一次初始的模糊，所以根据高斯模糊的性质： 

σ_init-第0层尺度；σ_pre-被相机镜头模糊后的尺度； 

高斯金字塔的组数： 

O＝[log₂(min(M,N))]-3 

其中，M、N分别为图像的行数和列数； 

DOG即相邻两尺度空间函数之差，用D(x,y,σ)来表示，即： 

D(x,y,σ)＝(G(x,y,kσ)-G(x,y,σ))*I(x,y)＝L(x,y,kσ)-L(x,y,σ) 

DOG金字塔的第1层的尺度因子与高斯金字塔的第1层的是一致的，其它阶也是一样； 

2.2极值点检测和关键点精确定位 

为了能够检测到DOG空间的极大值和极小值，DOG尺度空间除去最底层和最顶层的中间层中每个像素点都要与同一层相邻的8个像素点以及它上下两层邻近的9个像素点总共26个相邻像素点进行逐个比较，以确保在二维图像空间和尺度空间都能检测到局部极值； 

对尺度空间DOG函数进行曲线拟合，利用DOG所在尺度空间的Taylor展开式： 

求导并让方程等于零，可以得到极值点： 

为修正值，对应极值点，方程的值： 

将

的极值点视为低对比度的不稳定的极值点，进行剔除； 

DOG函数的主曲率通过计算在该点位置和尺度的2×2的Hessian矩阵H得到，导数可由相邻采样点的差来估计： 

令α为较大的特征值，β为小的特征值，α＝rβ,则： 

Tr(H)＝D_xx+D_yy＝α+βDet(H)＝D_xxD_yy-(D_xy)²＝αβ 

Lowe论证建议r=10，若

时将关键点保留，反之剔除； 

2.3关键点方向分配 

像素点的梯度表示： 

梯度幅值： 

梯度方向： 

关键点的方向分配实现步骤：(1)计算关键点邻域直方图的高斯权重参数；(2)生成含有36个柱的梯度直方图，其中范围为0～360度，每10度一个柱，由半径为图像区域生成，柱所代表的方向为像素点的梯度方向，柱的长短表示梯度的幅值。然后对梯度直方图做两次平滑滤波；(3)求取关键点方向；(4)二次拟合方向 直方图的Taylor展开式，精确关键点的方向； 

关键点主方向：梯度方向直方图的主峰值，即是特征点方向；关键点的辅方向：当存在其它的相当于梯度直方图中主峰值80%能量的峰值； 

2.4关键点描述子的生成 

128维关键点的特征向量生成步骤如下： 

(1)确定计算关键点描述子所需要的图像块的大小，图像区域表示关键点所在的尺度空间，其半径通过下式计算： 

其中σ_oct为关键点的组内尺度，d＝4； 

(2)将坐标移到关键点的主方向，旋转后的新坐标表示为： 

(3)在图像区域内计算每个像素点的梯度方向和幅值，然后对每个求得的梯度幅值乘以高斯权重参数，生成方向直方图： 

x_k：表示该像素点和中心点的列距离；y_k：表示该像素点和中心点的行距离；σ_w：表示描述子窗口的宽度3σ×直方图列数（取4）的一半； 

(4)在每个窗口宽度为2×2的区域内求得8个方向的梯度直方图，计算每个梯度方向的累加值，形成一个种子点。共生成16个种子点，每一个种子点得到一个8位的特征向量，这样就生成4×4×8共128维的特征向量； 

(5)描述子向量元素门限化：方向直方图每个方向上梯度幅值限制在一定门限值以下，一般门限值取0.2，描述子向量元素规范化： 

其中，W＝(w₁,w₂,...,w₁₂₈)为得到的128维特征向量，L＝(l₁,l₂,...,l₁₂₈)为规范化的向量，关键点描述子向量的规范化正是可去除满足此模型的光照影响，图像 各点的梯度是通过邻域像素相减得到，对于图像灰度值整体漂移，也能去除； 

设置窗口大小为标准SIFT窗口的1.67倍，把4×4的增大到6×6，产生6×6×8共288维向量描述子，增加正确匹配的点数，减弱噪声的影响。 

3.如权利要求1所述的基于SIFT和归一化互信息的SAR图像配准方法，其特征在于：步骤3所述的SIFT特征向量匹配，通过如下公式进行： 

参考图中关键点描述子：R_i＝(r_i1,r_i2,...,r_i128)，实时图中关键点描述子：S_i＝(s_i1,s_i2,...,s_i128)，任意两描述子相似性度量： 

要得到配对的关键点描述子，d(R_i,S_j)需满足： 

。

4.如权利要求1所述的基于SIFT和归一化互信息的SAR图像配准方法，其特征在于：步骤4中所述的单应矩阵模型的RANSAC策略，通过如下公式进行： 

4.1给定两幅图像I1，I2的初始匹配点集X1((x_i,y_i),i＝1,2,...,n)，X2((x_j,y_j),j＝1,2,...,n)，将它们分别写成齐次坐标的形式：X1((x_i,y_i,1),i＝1,2,...,n)，X2((x_j,y_j,1),j＝1,2,...,n)； 

4.2从匹配点对集合中，任选四个点对，计算单应矩阵H； 

4.3在所有匹配点对集合中，以欧氏距离为判据，选取所有符合d||HX1-X2||＜t的点对，将它们作为内点，不符合此条件的点对认为是误匹配点对，予以剔除，并记录满足此约束的内点数量，实验采用阈值t＝36； 

4.4重复循环4.2、4.3至N次，记录每一次的内点数量，实验采用N＝2000； 

4.5选取对应内点数最大的H，并把此次内点作为最终的内点，即正确匹配点对。 

5.如权利要求1所述的基于SIFT和归一化互信息的SAR图像配准方法，其特征在于：步骤4中所述的最小二乘法求解配准参数，通过如下公式进行： 

上式可以表示为：Ax＝b，用最小二乘法求得：x＝[A^TA]^-1A^Tb，进而可以求出，旋转角度：angle＝-(180/π)atan(a₂₁/a₁₁)，水平平移量：Δx＝t_x，垂直平移量：Δy＝t_y，尺度变化：

6.如权利要求1所述的基于SIFT和归一化互信息的SAR图像配准方法，其特征在于：步骤6中所述的基于归一化互信息的SAR图像配准方法，通过如下公式进行： 

6.1归一化互信息配准参数的初始值修正为angle、Δx和Δy； 

6.2利用Powell算法依据归一化互信息判别所求得的参数是否最优，若不是，则继续优化更新几何变换，直至获得最大归一化互信息值，搜索到满足精度要求的参数，达到实时图的最优变换，求得配准参数和输出配准后的图像； 

Powell算法步骤如下： 

(1)给定允差ε＞0，初始点X⁽⁰⁾和n个d^(1,1),d^(1,2),...,d^(1,n)线性无关的方向，置k＝1；其中，ε为允许误差，X⁽⁰⁾为初始点，d^(1,1),d^(1,2),...,d^(1,n)为n个线性无关的方向，k为循环次数； 

(2)置x^(k,0)＝x^(k-1)，从x^(k,0)出发，依次沿方向d^(k,1),d^(k,2),...,d^(k,n)做一次一维搜索，得到点x^(k,1),x^(k,2),...,x^(k,n)，求m，使得 

令d^(k,n+1)＝x^(k,n)-x^(k,0)，若||x^(k,n)-x^(k,0)||≤ε，则停止计算；否则，进行步骤(3)； 

其中，x^(k,0)为每次搜索的初始位置，d^(k,1),d^(k,2),...,d^(k,n)为n个线性无关的搜索方 向，x^(k,1),x^(k,2),...,x^(k,n)为得到的参数值，f(x)为互信息的值，d^(k,n+1)为新的搜索方向； 

(3)求λ_n+1，使得令x^(k+1,0)＝x^(k)＝x^(k,0)+λ_n+1d^(k,n+1)，若||x^(k)-x^(k-1)||≤ε，则停止计算，得到x^(k)；否则，进行步骤（4）；其中，λ_n+1，λ为搜索方向的斜率，x^(k)为得到的参数值； 

(4)若则令d^(k+1,j)＝d^(k,j)，j＝1,...,m-1，d^(k+1,j)＝d^(k,j+1)，j＝m,....,n，k＝k+1；否则，令d^(k+1,j)＝d^(k,j)，j＝1,...,n，k＝k+1，转到步骤(2)；其中，d^(k+1,j)为得到的新的搜索方向，k为循环次数，n为方向个。 

Images