CN106355576A - 基于mrf图像分割算法的sar图像配准方法 - Google Patents

Abstract

本发明提出了一种基于MRF图像分割算法的SAR图像配准方法，用于解决现有基于特征的SAR图像配准方法中存在的配准效率低及稳定性差的技术问题，实现步骤为：利用MRF图像分割算法对SAR图像参考图和待配准图分割；对分割后的参考图和待配准图进行区域截取；将截取后的分割参考图图像块和分割待配准图像块分别对应到SAR图像参考图和待配准图中；构建SAR图像参考图图像块的SAR‑Harris尺度空间和待配准图图像块的SAR‑Harris尺度空间；提取参考图图像块和待配准图像块的SARSIFT特征点，组成匹配点对集合；利用RANSAC算法去除误匹配点对；利用互信息方法，对匹配点对集合进行优化。

Description

基于MRF图像分割算法的SAR图像配准方法

技术领域

本发明属于图像处理技术领域，涉及一种SAR图像配准方法，具体涉及一种基于MRF图像分割算法的SAR图像配准方法，可用于图像融合和多时相图像变化检测等领域。

背景技术

图像配准是对取自不同时间、不同视角或不同传感器的同一景物的两幅或多幅图像进行匹配、叠加的过程。它的主要目的是消除或减少基准图像和待矫正图像之间由于成像条件不同所引起的几何形变，从而获得几何一致性(最佳空间位置匹配)的两幅图像。它在图像融合、多时相图像变化检测等领域都得到了广泛的应用。

合成孔径雷达(Synthetic Aperture Radar，SAR)是一种全天时、全天候的高分辨率微波遥感成像雷达，SAR成像弥补了红外成像、可见光成像的缺陷，是天基侦查监视系统中不可缺少的重要探测技术。虽然SAR图像有诸多好处，但由于SAR是相干系统，所以斑点噪声是SAR图像的固有性质。图像配准是不同的图像之间融合的基础，也是图像变化检测，提取感兴趣目标的基础，具有重要的意义和广泛的应用前景，其配准方法分为三类：基于灰度的配准、基于变换域的配准以及基于特征的配准。基于灰度的配准方法对尺寸较小、灰度变化不大的图像配准精度高，但对于视角、灰度、结构变化较大的遥感图像配准精度较差；基于变换域的配准方法易于硬件实现，而且抗低频噪声强，能做到实时的配准，但无法处理非线性形变，通常碰到高频噪声配准效果也不太好，所以基于灰度配准和基于变换域的配准方法都不适于SAR图像的配准。基于特征的配准方法，直接提取有效的特征进行匹配,抗噪性好，对灰度变化也不敏感。近年来，利用稳定特征描述子匹配成为图像配准领域研究的热点，其中尺度不变特征变换(Scale-invariant Feature Transform，SIFT)算法对图像配准效果较好，该算法包括以下步骤：建立高斯差分尺度空间；提取SIFT特征；筛选特征点；滤除错误匹配点；仿射变换；叠加得到最终的配准结果，该算法已经在SAR图像配准领域取得了不错的效果。例如，中国专利申请，申请公开号为CN103839265A，名称为“基于SIFT和归一化互信息的SAR图像配准方法”中，提出了一种基于SIFT与归一化互信息相结合的SAR图像配准方法，该方法在原基于SIFT方法配准的基础上，给定互信息配准一个初始值，采用归一化互信息方法对参考图和待配准图进行精配准，求得最终的配准参数，并输出配准后的图像。该方法存在不足之处是，在建立尺度空间和特征点提取的过程中计算量较大，导致图像的配准效率较低。

发明内容

本发明的目的在于克服上述现有技术存在的不足，提出了一种基于MRF图像分割算法的SAR图像配准方法，用于解决现有基于特征的SAR图像配准方法中存在的配准效率低的技术问题。

为实现上述目的，本发明采取的技术方案包括如下步骤：

(1)利用MRF图像分割方法，对SAR图像参考图和待配准图分别进行图像分割，实现步骤为：

(1a)设定MRF图像分割方法的分割类别参数，得到确定类别数目的MRF图像分割方法；

(1b)利用确定类别数目的MRF图像分割方法，对SAR图像参考图和待配准图进行图像分割，得到分割参考图和分割待配准图；

(2)对分割后参考图和分割待配准图分别进行区域截取，得到分割参考图的图像块和分割待配准图的图像块；

(3)将分割参考图的图像块对应到SAR图像参考图中，同时将分割待配准图的图像块对应到SAR图像的待配准图中，得到参考图图像块和待配准图图像块；

(4)利用SARSIFT算法，对参考图图像块和待配准图图像块分别进行处理，实现步骤为：

(4a)构建参考图图像块的SAR-Harris尺度空间S1,同时构建待配准图图像块的SAR-Harris尺度空间S2；

(4b)计算SAR-Harris尺度空间S1中各像素点梯度，同时比较每个像素点梯度与该像素点邻域内的像素点梯度的大小，得到SAR-Harris尺度空间S1的多个特征点；

(4c)计算SAR-Harris尺度空间S2中各像素点梯度，同时比较每个像素点梯度与该像素点邻域内的像素点梯度的大小，得到SAR-Harris尺度空间S2的多个特征点；

(4d)利用直方图统计法，逐一统计SAR-Harris尺度空间S1的多个特征点中每个特征点邻域内像素点的梯度幅值和方向，得到SAR-Harris尺度空间S1中每个特征点对应的尺度直方图，并将该直方图中最大的幅值所对应的方向作为SAR-Harris尺度空间S1特征点的主方向；

(4e)利用直方图统计法，逐一统计SAR-Harris尺度空间S2的多个特征点中每个特征点邻域内像素点的梯度幅值和方向，得到SAR-Harris尺度空间S2中每个特征点对应的尺度直方图，并将该直方图中最大的幅值所对应的方向作为SAR-Harris尺度空间S2特征点的主方向；

(4f)将SAR-Harris尺度空间S1和SAR-Harris尺度空间S2中每个特征点的邻域的坐标轴,旋转到该特征点所在尺度空间中的主方向上，得到每个特征点旋转后的邻域，计算该邻域8个方向的梯度方向直方图，并将每个梯度方向所对应的幅值作为该特征点的描述子，得到参考图图像块特征点描述子集合R1和待配准图图像块特征点描述子集合R2，将R1和R2的联合特征点描述子集合记为R＝{R1,R2}；

(5)利用RANSAC算法,去除联合特征点描述子集合R中的错误匹配点对，实现步骤为：

(5a)从联合特征点描述子集合R中随机选取一个包括p个匹配点对的RANSAC样本集，其中p≥4，并对该p个匹配点对进行仿射变换，得到变换矩阵M；

(5b)设定阈值J,计算联合特征点描述子集合R中剩余的点对经过变换矩阵M变换后的欧式距离z，判断欧式距离z是否小于阈值J，若是，将该点对放入RANSAC样本集合中，并将该样本集合作为初步匹配点对集合,否则，舍去该点对；

(6)利用互信息优化方法，对初步匹配点对集合进行优化,实现步骤为：

(6a)对互信息集合进行初始化，得到空的互信息集合；

(6b)从步骤(5b)中得到的初步匹配点集合中随机选取L个匹配点对，利用该L个匹配点对待配准图图像块进行仿射变换，并计算仿射变换后的待配准图图像块与参考图图像块的互信息；

(6c)将互信息添加到空集互信息集合中,得到包含一个互信息集合；

(6d)重复步骤(6b)～(6c),直到选取完所有的匹配点，得到包含多个互信息的互信息集合，并选出该包含多个互信息的互信息集合中最大值对应的点对作为最优匹配点对；

(6e)利用仿射变换函数对最优匹配点对进行仿射变换，得到最优仿射变换矩阵，并通过该最优的仿射变换矩阵对SAR待配准图进行仿射变换，得到仿射变换的SAR待配准图；

(6f)将仿射变换后的SAR图像待配准图与SAR图像参考图进行叠加，得到SAR图像待配准图与SAR图像参考图的配准图像。

本发明与现有技术相比，具有以下优点：

1.本发明由于在SAR-Harris尺度空间建立的过程中，采用了MRF图像分割方法对图像进行分割，该图像分割方法对图像噪声不敏感，有利于提高图像的分割质量，得到的特征明显、尺寸较小的SAR图像块可以更快地提取更有效的特征，与现有的基于特征的SAR图像配准方法相比，有效的提高了图像的配准效率。

2.本发明由于在图像特征匹配的过程中，采用SARSIFT算法对参考图图像块和待配准图图像块分别进行处理，减少了噪声对图像特征提取的影响，提取到了更多更准确的SAR图像特征，与现有的基于特征的SAR图像配准技术相比，有效的提高了图像配准的鲁棒性。

附图说明

图1是本发明的实现流程框图；

图2是本发明SAR图像参考图和得到的分割参考图；

图3是本发明SAR图像待配准图和得到的分割待配准图；

图4是本发明得到的SAR图像参考图图像块图和待配准图图像块图；

图5是本发明得到的最终的图像配准图。

具体实施方式

以下结合附图和实施例，对本发明作进一步详细说明：

参照图1，本发明包括如下步骤：

步骤1.利用MRF图像分割方法，对SAR图像参考图和待配准图分别进行图像分割，实现步骤为：

(1a)设定MRF图像分割方法的分割类别参数，为了使分割后的图像更容易区分，将分割类别设为3类，得到确定类别数目的MRF图像分割方法；

(1b)利用确定类别数目的MRF图像分割方法，对SAR图像参考图和待配准图进行图像分割，得到分割参考图和分割待配准图，实现步骤为：

(1b1)对SAR图像参考图和待配准图的每个像素分别进行类别初始化，得到SAR图像参考图和待配准图的初始化分割类别场；

(1b2)判断SAR图像参考图和待配准图类别场中当前像素点邻域中是否存在占支配地位的标记类x_p，即判断其中是邻域R_i标记为x_p的像素个数之和，如果成立，则令当前像素点的新标记类x_p；否则随机改变当前像素点的标记类，获得新标记类x_k；

(1b3)计算当前SAR图像参考图和待配准图当前标记类像素点和新标记类像素点的能量，判断新标记类像素点的能量是否低于当前标记类像素点的能量，是接受新标记，否则保持当前分割标记类不变，得到SAR图像参考图和待配准图的当前类别场和SAR图像参考图和待配准图的新标记类别场；

(1b4)设阈值G，计算SAR图像参考图的当前类别场的全局能量和新标记类别场的全局能量，同时计算SAR图像待配准图类别场的全局能量和新标记类别场的全局能量，判断SAR图像参考图类别场和SAR图像待配准图类别场的全局能量的变化量是否小于事先确定的阈值G，若是就将SAR图像参考图类别场和SAR图像待配准图类别场作为SAR图像参考图和SAR图像待配准图的分割图，否则迭代寻优过程直至结束。

步骤2.对分割后参考图和分割待配准图分别进行区域截取，截取分割图像中区间连续，特征显著的区域，得到分割参考图的图像块和分割待配准图的图像块；

步骤3.将分割参考图的图像块对应到SAR图像参考图中，同时将分割待配准图的图像块对应到SAR图像的待配准图中，得到参考图图像块和待配准图图像块；

步骤4.利用SARSIFT算法，对参考图图像块和待配准图图像块分别进行处理，实现步骤为：

(4a)构建参考图图像块的SAR-Harris尺度空间S1,同时构建待配准图图像块的SAR-Harris尺度空间S2，实现步骤为：

(4a1)利用下式计算得到SAR-Harris尺度空间的滤波核T：

$T=\underset{x=R,y=R+}{\∫\∫}I(a+x,b+y)\×e^{-\frac{\left|x\right|+\left|y\right|}{\mathrm{\α}}},$

其中是指数平均加权滤波核，α表示SAR-Harris尺度空间的尺度参数，x和y是图像I的像素点的邻域空间，在该像素点水平方向和垂直方向范围均为[-R,R]；

(4a2)利用SAR-Harris尺度空间滤波核T与参考图图像块和待配准图图像块分别进行卷积就得到SAR参考图图像块和待配准图图像块的SAR-Harris尺度空间图像，卷积过程如下所示：

S1＝H_α*I_r，

S2＝H_α*I_s，

其中，S1，S2分别表示参考图图像块和待配准图图像块对应的SAR-Harris尺度空间图像，H_α表示尺度参数为α对应的指数平均加权滤波核，*表示二维卷积操作，I_r,I_s分别表示参考图图像块和待配准图像块；

得到参考图图像块的SAR-Harris尺度空间S1和待配准图图像块的SAR-Harris尺度空间S2。

(4b)计算SAR-Harris尺度空间S1中各像素点梯度，同时比较每个像素点梯度与该像素点邻域内的像素点梯度的大小，得到SAR-Harris尺度空间S1的多个特征点；

(4c)计算SAR-Harris尺度空间S2中各像素点梯度，同时比较每个像素点梯度与该像素点邻域内的像素点梯度的大小，得到SAR-Harris尺度空间S2的多个特征点；

(4d)利用直方图统计法，逐一统计SAR-Harris尺度空间S1的多个特征点中每个特征点邻域内像素点的梯度幅值和方向，得到SAR-Harris尺度空间S1中每个特征点对应的尺度直方图，并将该直方图中最大的幅值所对应的方向作为SAR-Harris尺度空间S1特征点的主方向；

(4e)利用直方图统计法，逐一统计SAR-Harris尺度空间S2的多个特征点中每个特征点邻域内像素点的梯度幅值和方向，得到SAR-Harris尺度空间S2中每个特征点对应的尺度直方图，并将该直方图中最大的幅值所对应的方向作为SAR-Harris尺度空间S2特征点的主方向；

(4f)将SAR-Harris尺度空间S1和SAR-Harris尺度空间S2中每个特征点的邻域的坐标轴,旋转到该特征点所在尺度空间中的主方向上，得到每个特征点旋转后的邻域，计算该邻域8个方向的梯度方向直方图，并将每个梯度方向所对应的幅值作为该特征点的描述子，得到参考图图像块特征点描述子集合R1和待配准图图像块特征点描述子集合R2，将R1和R2的联合特征点描述子集合记为R＝{R1,R2}；

步骤5.利用RANSAC算法,去除联合特征点描述子集合R中的错误匹配点对，实现步骤为：

(5a)从联合特征点描述子集合R中随机选取一个包括p个匹配点对的RANSAC样本集，其中p≥4，并对该p个匹配点对进行仿射变换，得到变换矩阵M；

(5b)设定阈值J,计算联合特征点描述子集合R中剩余的点对经过变换矩阵M变换后的欧式距离z，判断欧式距离z是否小于阈值J，若是，将该点对放入RANSAC样本集合中，并将该样本集合作为初步匹配点对集合,否则，舍去该点对；

步骤6.利用互信息优化方法，对初步匹配点对集合进行优化,实现步骤为：

(6a)对互信息集合进行初始化，得到空的互信息集合；

(6b)从步骤(5b)中得到的初步匹配点集合中随机选取L个匹配点对，利用该L个匹配点对待配准图图像块进行仿射变换，并计算仿射变换后的待配准图图像块与参考图图像块的互信息；

(6c)将互信息添加到空集互信息集合中,得到包含一个互信息集合；

(6d)重复步骤(6b)～(6c),直到选取完所有的匹配点，得到包含多个互信息的互信息集合，并选出该包含多个互信息的互信息集合中最大值对应的点对作为最优匹配点对；

(6e)利用仿射变换函数对最优匹配点对进行仿射变换，得到最优仿射变换矩阵，并通过该最优的仿射变换矩阵对SAR待配准图进行仿射变换，得到仿射变换的SAR待配准图；

(6f)将仿射变换后的SAR图像待配准图与SAR图像参考图进行叠加，得到SAR图像待配准图与SAR图像参考图的配准图像。

下面结合仿真实验，对本发明的技术效果作进一步的描述。

1、仿真实验条件：

采用一幅大小为1000×954的SAR图像参考图,如图2(a)所示和一幅大小为1047×1000的SAR图像待配准图，如图3(a)所示；硬件平台为：Intel(R)Core(TM)i5-4210U、4GBRAM，软件平台：MATLAB R2015a。

2、实验内容与结果：

仿真1，对SAR图像参考图分割过程进行仿真，其结果如图2(b)所示，图2(a)表示SAR图像参考图，图2(b)表示分割参考图；

参照图2可以看出分割参考图与SAR图像参考图相比，分割参考图会形成特征明显的分割图像区域，对分割参考图中的特征明显的连接区域进行截取，可以获得特征显著的分割参考图图像块；

仿真2，对SAR图像待配准图分割过程进行仿真，其结果如图3(b)所示，图3(a)表示SAR图像待配准图，图3(b)表示分割待配准图；

参照图3可以看出分割待配准图与SAR图像参考图相比，分割待配准图会形成特征明显的分割图像区域，对分割待配准图中的特征明显的连接区域进行截取，可以获得特征显著的分割待配准图像块；

仿真3，对区域截取的图像块对应到原图中的过程进行仿真，其结果如图4所示，图4(a)表示SAR图像参考图分割图像块，图4(b)表示SAR图像待配准图分割图像块；

仿真4，对SAR图像参考图和待配准图的叠加过程进行仿真，其结果如图5所示，参照图4可以看出SAR图像参考图图像块和待配准图图像块与原SAR图像参考图和待配准图相比尺寸较小，而且特征更加明显，这样在建立尺度空间和特征提取的过程中，计算量将显著减少，最终提高了图像配准的效率。

Claims (3)

1.基于MRF图像分割算法的SAR图像配准方法：包括以下步骤：

(1)利用MRF图像分割方法，对SAR图像参考图和待配准图分别进行图像分割，实现步骤为：

(1a)设定MRF图像分割方法的分割类别参数，得到确定类别数目的MRF图像分割方法；

(1b)利用确定类别数目的MRF图像分割方法，对SAR图像参考图和待配准图进行图像分割，得到分割参考图和分割待配准图；

(2)对分割后参考图和分割待配准图分别进行区域截取，得到分割参考图的图像块和分割待配准图的图像块；

(3)将分割参考图的图像块对应到SAR图像参考图中，同时将分割待配准图的图像块对应到SAR图像的待配准图中，得到参考图图像块和待配准图图像块；

(4)利用SARSIFT算法，对参考图图像块和待配准图图像块分别进行处理，实现步骤为：

(4a)构建参考图图像块的SAR-Harris尺度空间S1,同时构建待配准图图像块的SAR-Harris尺度空间S2；

(4b)计算SAR-Harris尺度空间S1中各像素点梯度，同时比较每个像素点梯度与该像素点邻域内的像素点梯度的大小，得到SAR-Harris尺度空间S1的多个特征点；

(4c)计算SAR-Harris尺度空间S2中各像素点梯度，同时比较每个像素点梯度与该像素点邻域内的像素点梯度的大小，得到SAR-Harris尺度空间S2的多个特征点；

(4d)利用直方图统计法，逐一统计SAR-Harris尺度空间S1的多个特征点中每个特征点邻域内像素点的梯度幅值和方向，得到SAR-Harris尺度空间S1中每个特征点对应的尺度直方图，并将该直方图中最大的幅值所对应的方向作为SAR-Harris尺度空间S1特征点的主方向；

(4e)利用直方图统计法，逐一统计SAR-Harris尺度空间S2的多个特征点中每个特征点邻域内像素点的梯度幅值和方向，得到SAR-Harris尺度空间S2中每个特征点对应的尺度直方图，并将该直方图中最大的幅值所对应的方向作为SAR-Harris尺度空间S2特征点的主方向；

(4f)将SAR-Harris尺度空间S1和SAR-Harris尺度空间S2中每个特征点的邻域的坐标轴,旋转到该特征点所在尺度空间中的主方向上，得到每个特征点旋转后的邻域，计算该邻域8个方向的梯度方向直方图，并将每个梯度方向所对应的幅值作为该特征点的描述子，得到参考图图像块特征点描述子集合R1和待配准图图像块特征点描述子集合R2，将R1和R2的联合特征点描述子集合记为R＝{R1,R2}；

(5)利用RANSAC算法,去除联合特征点描述子集合R中的错误匹配点对，实现步骤为：

(5a)从联合特征点描述子集合R中随机选取一个包括p个匹配点对的RANSAC样本集，其中p≥4，并对该p个匹配点对进行仿射变换，得到变换矩阵M；

(5b)设定阈值J,计算联合特征点描述子集合R中剩余的点对经过变换矩阵M变换后的欧式距离z，判断欧式距离z是否小于阈值J，若是，将该点对放入RANSAC样本集合中，并将该样本集合作为初步匹配点对集合,否则，舍去该点对；

(6)利用互信息优化方法，对初步匹配点对集合进行优化,实现步骤为：

(6a)对互信息集合进行初始化，得到空的互信息集合；

(6b)从步骤(5b)中得到的初步匹配点集合中随机选取L个匹配点对，利用该L个匹配点对待配准图图像块进行仿射变换，并计算仿射变换后的待配准图图像块与参考图图像块的互信息；

(6c)将互信息添加到空集互信息集合中,得到包含一个互信息集合；

(6d)重复步骤(6b)～(6c),直到选取完所有的匹配点，得到包含多个互信息的互信息集合，并选出该包含多个互信息的互信息集合中最大值对应的点对作为最优匹配点对；

(6e)利用仿射变换函数对最优匹配点对进行仿射变换，得到最优仿射变换矩阵，并通过该最优的仿射变换矩阵对SAR待配准图进行仿射变换，得到仿射变换的SAR待配准图；

(6f)将仿射变换后的SAR图像待配准图与SAR图像参考图进行叠加，得到SAR图像待配准图与SAR图像参考图的配准图像。

2.根据权利要求1所述的基于MRF图像分割算法的SAR图像配准方法，其特征在于：步骤(1b)所述的利用确定类别数目的MRF图像分割方法，对SAR图像参考图和待配准图进行图像分割，实现步骤为：

(1b1)对SAR图像参考图和待配准图的每个像素分别进行类别初始化，得到SAR图像参考图和待配准图的初始化分割类别场；

(1b2)判断SAR图像参考图和待配准图类别场中当前像素点邻域中是否存在占支配地位的标记类x_p，即判断其中是邻域R_i标记为x_p的像素个数之和，如果成立，则令当前像素点的新标记类为x_p；否则随机改变当前像素点的标记类，获得新标记类x_k；

(1b3)计算当前SAR图像参考图和待配准图当前标记类像素点和新标记类像素点的能量，判断新标记类像素点的能量是否低于当前标记类像素点的能量，若是，接受新标记类，否则，保持当前分割标记类不变，得到SAR图像参考图和待配准图的当前类别场和SAR图像参考图和待配准图的新标记类别场；

(1b4)设阈值G，计算SAR图像参考图的当前类别场的全局能量和新标记类别场的全局能量，同时计算SAR图像待配准图类别场的全局能量和新标记类别场的全局能量，判断SAR图像参考图类别场和SAR图像待配准图类别场的全局能量的变化量是否小于事先确定的阈值G，若是，将SAR图像参考图新标记类别场和SAR图像待配准图新标记类别场作为分割参考图和分割待配准图，否则，重复步骤(1b2)～(1b4)。

3.根据权利要求1所述的基于MRF图像分割算法的SAR图像配准方法，其特征在于：步骤(4a)所述的构建SAR图像参考图图像块的SAR-Harris尺度空间S1,同时构建待配准图图像块的SAR-Harris尺度空间S2，实现步骤为：

(4a1)利用指数平均加权滤波算法，计算待构建SAR-Harris尺度空间的滤波核T_α：

$T_{\mathrm{\α}}=\underset{x=R,y=R+}{\∫\∫}T(a+x,b+y)\×e^{-\frac{\left|x\right|+\left|y\right|}{\mathrm{\α}}},$

其中是指数平均加权滤波核，α表示SAR-Harris尺度空间的尺度参数，x和y是图像I的像素点的邻域空间，在该像素点水平方向和垂直方向范围均为[-R,R]；

(4a2)利用待构建SAR-Harris尺度空间的滤波核T_α与参考图图像块和待配准图图像块分别进行卷积，得到SAR参考图图像块的SAR-Harris尺度空间图像s1_iα和待配准图图像块的SAR-Harris尺度空间图像s2_iα，卷积过程：

S1＝T_α*I_r，

S2＝T_α*I_s，

其中，S1，S2分别表示参考图图像块和待配准图图像块对应的SAR-Harris尺度空间图像，*表示二维卷积操作，I_r,I_s分别表示参考图图像块和待配准图像块；

(4a3)改变尺度参数α，并重复执行步骤(4a2)，得到多个不同尺度参数的参考图图像块对应的SAR-Harris尺度空间图像s1_iα和多个待配准图图像块对应的SAR-Harris尺度空间图像s2_iα，并将该两个尺度空间图像进行叠加，得到参考图图像块的SAR-Harris尺度空间S1和待配准图图像块的SAR-Harris尺度空间S2。