CN104392462A - 一种基于显著分割子区域对的sar图像配准方法 - Google Patents

Abstract

本发明提供了一种基于显著分割子区域对的SAR图像配准方法，步骤是：1)输入同一地区两个时相的SAR图像；2)对两幅SAR图像进行均值滤波和直方图均衡化抑制相干斑；3)对两幅滤波均衡化后图像进行基于频谱残差的显著区域分割；4)对两幅图像的显著分割区域利用直方图众数排序后的灰度级序列计算最小欧氏距离获得匹配区域；5)利用改进后的综合学习粒子群算法对匹配区域计算配准参数；6)将配准参数利用像素点在子区域和原图像的坐标关系转化为原图像对的配准参数，得到最终配准结果。本发明解决现有方法在区域分割一致性不理想情况下的配准问题和基于灰度信息的配准计算量大的问题，用于变化检测、地物分类、地图修正、数据融合等。

Description

一种基于显著分割子区域对的SAR图像配准方法

技术领域

本发明属于数字图像处理技术领域，利用两幅图像中相似区域信息进行配准，是变化检测、地物分类、地图修正、数据融合等应用的基础；具体是一种基于显著分割子区域对的SAR图像配准方法。

背景技术

合成孔径雷达(SAR)具有全天候、全天时、对云层和地表的穿透等特点被广泛应用于水文、地矿、环境监测、军事等领域。但是由于原始数据拍摄条件的不同，导致图像与图像之间存在相对的几何差异，校正图像差异的图像配准技术成为进一步处理图像的前提条件。

基于区域进行配准的思路，大多首先对两个图像进行分割，然后对两幅图像的分割区域进行匹配，然后对匹配区域提取控制点，或者无需进行区域匹配，直接利用闭合区域提取控制点，之后对控制点进行匹配，进而利用控制点进行配准。

2007年Wang Shuang等人利用Hu不变矩和边缘链码作为区域的特征，利用上述区域特征进行区域匹配。将匹配区域的重心以及边缘角点看做控制点，将利用最小化距离分类的准则建立控制点对，在找到足够正确的点后利用最小二乘法对几何变换模型参数进行估计并完成整个变换。

2010年Wang Zhenhua等人首先利用几何约束限制和Hu不变矩匹配区域，用匹配区域的质心作为控制点实现粗配准，第二步以Harris角点和交叉点作为特征点，在粗配准的基础上每个特征点在相应邻域内以互信息和相关系数作为相似性度量参数采用多重约束来搜索匹配点。

2012年Xiong Boli等人首先利用主动轮廓模型(GAC)提取图像中闭合区域。之后利用多边形拟合定位区域边界上的特征点，然后用几何哈希方法计算匹配点对，进而计算配准参数。

上述方法都存在一个隐形前提，即两幅图像的分割区域一致性很高。当分割区域一致性程序不高时，将会给描述区域形状信息的Hu不变矩、边缘链码和区域面积、区域周长等几何信息进行区域匹配时带来干扰，而基于区域结构或者边界的质心以及多边形拟合的控制点也将无法很好的对应，因而会影响上述方法的配准效果。而对于同一SAR系统在不同时间拍摄的SAR图像来说，由于时间点的不同，同一种分割方法在相同参数设置的情况下，是非常有可能出现分割结果不一致性的情况的。

发明内容

本发明的目的在于克服上述现有技术配准方法的不足，提出如下配准方法：首先利用基于频谱残差的显著区域分割，之后利用直方图众数排序的方法匹配分割区域，最后对匹配显著区域采用以比值熵为相似性度量函数，以改进后的综合学习粒子群算法为搜索算法的方法进行快速配准，之后将匹配显著区域的配准参数转化为原图像对的配准参数，从而有效提高配准的效率、精度和稳定性。

为实现上述目的，本发明提供了一种基于显著分割子区域对的SAR图像配准方法，其技术方案是：一种基于显著分割子区域对的SAR图像配准方法，包括如下步骤：

步骤1，输入两幅在不同时相获取的同一地区的大小为N行N列的SAR图像，分别记为I和J，图像的灰度级范围为0～255；

步骤2，对时相1图像I和时相2图像J分别进行均值滤波，然后再进行直方图均衡化，得到两幅滤波均衡化后图像，分别记为E和F，其中，均值滤波的邻域大小为m×m个像素，m的取值范围为3～8；

步骤3，对滤波均衡化后图像E进行基于频谱残差的显著区域分割，将得到的连通区域中面积大于A的区域分别标记为EA₁、EA₂、…、EA_r、…、EA_re，其中EA_r表示第r个面积大于A的区域，re表示对滤波均衡化后图像E进行基于频谱残差的显著区域分割后得到的连通区域中面积大于面积阈值A的区域个数，1≤r≤re；对滤波均衡化后图像F进行基于频谱残差的显著区域分割，将得到的连通区域中面积大于A的区域分别标记为FA₁、FA₂、…、FA_t、…、FA_tf，其中FA_t表示连通区域中第t个面积大于A的区域，tf表示对滤波均衡化后图像F进行基于频谱残差的显著区域分割后得到的连通区域中面积大于面积阈值A的区域个数，1≤t≤tf；

步骤4，利用直方图众数排序的方法计算时相1的所有分割区域和时相2的所有分割区域中欧氏距离最小的一对匹配区域，在这对匹配区域中，时相1的分割区域记为EA_mj，对应的图像块记为IA_mj，时相2的分割区域记为FA_mj，对应的图像块记为JA_mj；

步骤5，利用改进后的综合学习粒子群算法，计算时相2图像块JA_mj相对于时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy；

步骤6，根据时相2图像块JA_mj相对于时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx、Y轴方向平移参数Δy计算时相2图像J相对于时相1图像I的几何变换矩阵TF；

步骤7，根据几何变换矩阵TF，对时相2图像进行几何变换，其中的插值方法采用双线性插值方法。

上述步骤3具体按如下步骤进行:

3a)对滤波均衡化后图像E根据基于频谱残差的显著图计算方法计算其显著图G；

3b)对显著图G进行阈值分割，得到二值图像Gb，将二值图像Gb中大于阈值的像素赋值为1，小于阈值的像素赋值为0，若大于阈值的像素所形成区域的总面积比小于阈值的像素所形成区域的总面积大，则将大于阈值的像素的值赋为0，同时将小于阈值的像素的值赋为1，否则不做处理；得到一幅阈值分割图像Ew，其中，阈值的取值范围为0.02～0.05；

3c)将阈值分割图像Ew中连通区域中面积大于面积阈值A的区域分别标记为EA₁、EA₂、…、EA_r、…、EA_re，其中面积阈值A≥η×N²，η的取值范围为0.01～0.03，N²为图像的像素个数，EA_r表示第r个面积大于A的区域，re为连通区域中面积大于面积阈值A的区域个数，1≤r≤re；

3d)对滤波均衡化后图像F进行与步骤3a)到步骤3b)同样地处理，将得到的阈值分割图像的连通区域中面积大于A的区域标记为FA₁、FA₂、…、FA_t、…、FA_tf，其中FA_t表示连通区域中第t个面积大于A的区域，tf表示连通区域中面积大于面积阈值A的区域个数，1≤t≤tf。

上述步骤4具体按如下步骤进行：

4a)计算分割区域EA₁的灰度直方图，得到其灰度级像素数目矩阵AN_EA1＝[AN₀，AN₁，AN₂，…，AN_j，…AN₂₅₅]，其中AN_j表示该分割区域的第j个灰度级对应的像素数目，0≤j≤255，该分割区域的灰度级像素数目矩阵所对应的灰度级矩阵为L_EA1＝[0，1，2，…，j，…255]；将灰度级像素数目矩阵AN_EA1进行降序排列，得到灰度级像素数目降序矩阵AN_EA1s，将灰度级矩阵L_EA1按照AN_EA1s中的顺序重新排序，得到排序后灰度级矩阵L_EA1s＝[a₀，a₁，a₂，…，a_j，…，a₂₅₅]，其中a_j表示排序后灰度级矩阵L_EA1s的第j个矩阵元素，为灰度级像素数目降序矩阵AN_EA1s中的第j个矩阵元素所对应的灰度级值，0≤a_j≤255；

4b)对时相1的滤波均衡化后图像E中其余的分割区域EA₂、EA₃、…、EA_r、…、EA_re进行与步骤4a)同样地处理，分别得到对应的排序后灰度级矩阵L_EA2s、L_EA3s、…、L_EArs、…、L_EAres；同样地，对时相2的滤波均衡化后图像F中的所有的分割区域FA₁、FA₂、…、FA_t、…、FA_tf，进行与步骤4a)同样地处理，分别得到对应的排序后灰度级矩阵L_FA1s、L_FA2s、…、L_FAts、…、L_FAtfs；

4c)计算时相1中第j个排序后灰度级矩阵L_EAjs与时相2中各个排序后灰度级矩阵L_FA1s、L_FA2s、…、L_FAts、…、L_FAtfs的欧氏距离，将这些欧氏距离中的最小值所对应的时相2的排序后灰度级矩阵记为L_FAsj，其对应的时相2分割区域记为FA_sj，该最小欧氏距离记为D_sj，其中sj的取值范围为[1，re]；

4d)重复步骤4a)到4c)计算完时相1中所有的排序后灰度级矩阵所对应的最小欧氏距离，得到最小欧氏距离矩阵DS＝[D_s1D_s2…D_sr…D_sre]，最小欧氏距离矩阵DS中的最小值记为Ds_mj，其对应的时相1和时相2的分割区域分别为EA_mj和FA_mj，即为一对匹配区域对，它们对应的图像块分别记为IA_mj和JA_mj。

上述步骤5具体按如下步骤进行：

5a)初始化粒子群的代数Gen和粒子数目Ps，其中，Gen的取值范围为Gen≥10，粒子数目Ps的取值范围为Ps≥150；每个粒子在逐代进化过程中所产生的适应度值中的最大值称为该粒子的历史最大适应度值，该粒子的历史最大适应度值对应的位置称为该粒子的历史最优粒子位置，所有粒子在进化过程所能产生的适应度值中的最大值称为最大适应度值，最大适应度值对应的粒子位置称为最优粒子位置；初始化每个粒子历史最大适应度值保持不变的代数阈值为M，M的取值范围为1～3；

5b)一个粒子的粒子位置为一个三维向量，其第一维参数取值范围为[-90,90]，第二维参数取值范围为[-row,row]，第三维参数取值范围为[-col,col]，row和col分别为时相2图像块JA_mj的行数和列数；同样地，粒子速度也是一个由三个参数构成的三维向量，每一维参数的取值范围均为[0,1]；在粒子位置和粒子速度的每一维参数取值范围内随机初始化第一代种群的每个粒子位置和粒子速度中的每一维参数；

5c)将第i个粒子的历史最大适应度值保持不变的代数fg_i初始化为0；由第i个粒子的位置pos_i(θ_i,Δx_i,Δy_i)，可建立几何变换矩阵 ${\mathrm{tf}}_i=\left[\begin{array}{ccc}\cos \left({\mathrm{\θ}}_i\right)& \sin \left({\mathrm{\θ}}_i\right)& {\mathrm{\Δx}}_i\\ -\sin \left({\mathrm{\θ}}_i\right)& \cos \left({\mathrm{\θ}}_i\right)& {\mathrm{\Δy}}_i\\ 0& 0& 1\end{array}\right];$ 对时相2图像块JA_mj根据几何变换矩阵tf_i进行几何变换，得到时相2图像块的变换后图像块JT_i，进行几何变换所用的插值方法为双线性插值方法；

5d)计算第i个粒子的适应度值公式如下：

f(θ_i,Δx_i,Δy_i)＝-H(RI)

其中，图像RI＝{RI(x,y)|1≤x≤col,1≤y≤row}，其中，像素(x,y)的灰度值RI(x,y)计算公式如下：

其中，IA_mj(x,y)和JT_i(x,y)分别表示图像块IA_mj和JT_i中坐标位于(x,y)处的像素值；max(·)是取最大值函数；rand_(0-255)表示取值范围为[0,255]的随机数；H(RI)为计算图像RI的熵；

5e)将第一代种群粒子中适应度值最大的粒子的位置记为最优粒子位置Gb，其适应度值记为最大适应度值FG；将第i个粒子的历史最优位置Pb_i初始化为其位置,每个粒子的历史最大适应度值FP_i初始化为其适应度值；

5f)对于第一代种群的第i个粒子，i的取值范围为1,2,…Ps，如果其历史最大适应度值保持不变的代数fg_i大于每个粒子历史最大适应度值保持不变的代数阈值M，则按普通粒子群学习策略更新第i个粒子的位置；如果第i个粒子的历史最大适应度值保持不变的代数fg_i大于Mn，则除做上述位置更新之外，执行步骤5g)的克隆选择操作，其中Mn≥M；

5g)将第i个粒子的粒子位置X_i克隆GN_i次，形成GN_i个和粒子位置X_i相同的克隆个体；其中，GN_i≥50，设置变异概率Pm＝0.1，克隆个体的序号数ii的范围1≤ii≤GN_i；

5h)对于任意的第ii个克隆个体Gb_ii按如下公式进行变异：

${\mathrm{Gb}}_{\mathrm{ii}}\←\left\{\begin{array}{cc}{\mathrm{Gb}}_{\mathrm{iin}}& {\mathrm{rand}}_{(0-1)}>\mathrm{Pm}\\ {\mathrm{Gb}}_{\mathrm{ii}}& {\mathrm{rand}}_{(0-1)}\≤\mathrm{Pm}\end{array}\right.$

其中，Gb_iin是一个1×3的矩阵，其第d个矩阵元素取为随机数rand_(m)，rand_(m)的取值范围为[min^d,max^d]，min^d和max^d分别表示第d维参数取值范围的下界和上界，rand_(0-1)表示取随机数，取值范围为[0,1]；

5i)将每一个变异后的克隆个体视为一个粒子，按照步骤5d)中的粒子适应度值公式计算其适应度值，将所有克隆个体适应度值中的最大值所对应的克隆个体记为Gb_ic；如果克隆个体Gb_ic的适应度值大于第i个粒子的适应度值，则将第i粒子的粒子位置X_i更新为克隆个体Gb_ic的粒子位置；

对更新后的第i个粒子按照步骤5d)中的粒子适应度值公式计算其适应度值，如果该适应度值大于最大适应度值FG，则将最大适应度值FG更新为该适应度值，最优粒子位置Gb更新为该粒子的粒子位置；如果该适应度值大于第i个粒子的历史最大适应度值，则将第i个粒子的历史最大适应度值FP_i更新为该适应度值，第i个粒子的历史最优粒子位置Pb_i更新为该粒子的粒子位置；

5j)运用综合学习策略更新第i个粒子的粒子位置；

5k)重复上述步骤5f)到5j)，将计算完所有粒子后得到的最优粒子位置Gb克隆GN次，形成GN个和最优粒子位置Gb相同的克隆个体，GN≥200，设置变异概率Pm＝0.1，克隆个体的序号数ii的范围1≤ii≤GN；执行步骤5h)；然后，将每一个变异后的克隆个体视为一个粒子，按照步骤5d)中的粒子适应度值公式计算其适应度值，所有克隆个体适应度值中的最大值所对应的克隆个体为Gbc，如果该克隆个体Gbc的适应度值大于最大适应度值FG，则将最大适应度值FG更新为该适应度值，最优粒子位置Gb更新为克隆个体Gbc；

5l)重复上述5f)到5k)的过程，依次计算第2、3、…直至第Gen代种群；最优粒子位置Gb的第一维参数、第二维参数、第三维参数分别对应时相2图像块JA_mj相对时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy。

上述步骤6具体按如下步骤进行：

6a)根据时相2图像块JA_mj相对时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy，构建时相2图像块相对时相1图像块的变换矩阵tf：

$\mathrm{tf}=\left[\begin{array}{ccc}\cos \left(\mathrm{\θ}\right)& \sin \left(\mathrm{\θ}\right)& \mathrm{\Δx}\\ -\sin \left(\mathrm{\θ}\right)& \cos \left(\mathrm{\θ}\right)& \mathrm{\Δy}\\ 0& 0& 1\end{array}\right];$

6b)从时相2图像块JA_mj选出3个点JA_p1(JA_p1x,JA_p1y)、JA_p2(JA_p2x,JA_p2y)、JA_p3(JA_p3x,JA_p3y)，其坐标位置按如下方式计算：

JA_p1x＝0.25×JAc，JA_p1y＝0.5×JAr

JA_p2x＝0.5×JAc，JA_p2y＝0.25×JAr

JA_p3x＝0.75×JAc，JA_p3y＝0.5×JAr

其中JAr和JAc分别为时相2图像块JA_mj的行数和列数；

根据变换矩阵tf可得到时相1图像块IA_mj中对应的3个点IA_p1(IA_p1x,IA_p1y)、IA_p2(IA_p2x,,IA_p2y)、IA_p3(IA_p3x,IA_p3y)，计算方式如下：

$\left[\begin{array}{ccc}{\mathrm{IA}}_{p1x}& {\mathrm{IA}}_{p2x}& {\mathrm{IA}}_{p3x}\\ {\mathrm{IA}}_{p1y}& {\mathrm{IA}}_{p2y}& {\mathrm{IA}}_{p3y}\\ 1& 1& 1\end{array}\right]=\mathrm{tf}\×\left[\begin{array}{ccc}{\mathrm{JA}}_{p1x}& {\mathrm{JA}}_{p2x}& {\mathrm{JA}}_{p3x}\\ {\mathrm{JA}}_{p1y}& {\mathrm{JA}}_{p2y}& {\mathrm{JA}}_{p3y}\\ 1& 1& 1\end{array}\right];$

6c)计算时相1图像块IA_mj中点IA_p1(IA_p1x,IA_p1y)、IA_p2(IA_p2x,,IA_p2y)、IA_p3(IA_p3x,IA_p3y)在时相1图像I中对应的点I_p1(I_p1x,I_p1y)、I_p2(I_p2x,,I_p2y)、I_p3(I_p3x,I_p3y)，计算时相2图像块JA_mj中点JA_p1(JA_p1x,JA_p1y)、JA_p2(JA_p2x,,JA_p2y)、JA_p3(JA_p3x,JA_p3y)在时相2图像J中对应的点J_p1(J_p1x,J_p1y)、J_p2(J_p2x,J_p2y)、J_p3(J_p3x,J_p3y)，计算公式如下：

IA_p1x＝IA_p1x+Ix0-1,IA_p1y＝IA_p1y+Iy0-1

IA_p2x＝IA_p2x+Ix0-1,IA_p2y＝IA_p2y+Iy0-1

IA_p3x＝IA_p3x+Ix0-1,IA_p3y＝IA_p3y+Iy0-1

JA_p1x＝JA_p1x+Jx0-1,JA_p1y＝JA_p1y+Jy0-1

JA_p2x＝JA_p2x+Jx0-1,JA_p2y＝JA_p2y+Jy0-1

JA_p3x＝JA_p3x+Jx0-1,JA_p3y＝JA_p3y+Jy0-1

其中(Ix0，Iy0)是时相1图像块IA_mj中坐标位置(1,1)的像素点在时相1图像I中的坐标；(Jx0，Jy0)是时相2图像块JA_mj中坐标位置(1,1)的像素点在时相2图像J中的坐标；

6d)根据时相1图像I中的点I_p1(I_p1x,I_p1y)、I_p2(I_p2x,,I_p2y)、I_p3(I_p3x,I_p3y)和时相2图像J中的点J_p1(J_p1x,J_p1y)、J_p2(J_p2x,J_p2y)、J_p3(J_p3x,J_p3y)计算时相2图像J相对时相1图像I的几何变换矩阵TF，计算公式如下：

$\mathrm{TF}=\frac{\left[\begin{array}{ccc}{I}_{p1x}& I_{p2x}& I_{p3x}\\ I_{p1y}& I_{p2y}& I_{p3y}\\ 1& 1& 1\end{array}\right]}{\left[\begin{array}{ccc}{J}_{p1x}& J_{p2x}& J_{p3x}\\ J_{p1y}& J_{p2y}& J_{p3y}\\ 1& 1& 1\end{array}\right]}.$

本发明的有益效果：与现有的技术相比本发明具有以下优点：

a)本发明从显著区域的角度对图像进行有效分割，有效降低区域分割的复杂性，进而利用两幅图像共有的显著区域或非显著区域对图像进行有效配准。

b)本发明提出了基于灰度直方图的旋转不变性度量方式，可以在区域分割一致性并不是很好的情况下实现分割区域的匹配。

c)本发明提出了通过衡量两幅图像的比值图像的熵来度量两幅图像相似性的方式，可以在优化时作为有效的目标函数使用。并将综合学习粒子群算法(CLPSO)和克隆选择算法(CSA)相结合，实现了更为有效的优化。

d)本发明提出一种利用像素点在子区域和原图像的坐标，有效的将子区域对变换矩阵转化为整幅图像对变换矩阵的方法，避免了直接进行参数转化时几何关系方面的分析和计算。

附图说明

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

图2是本发明仿真使用的不同时相的SAR图像。

图3是本发明实验中基于频谱残差的区域分割结果。

图4是本发明实验中的匹配区域对。

图5是本发明的匹配子区域对的配准结果。

图6是本发明实验中不同时相SAR图像的最终的配准结果。

图7是本发明实验中不同时相SAR图像配准后的图像拼接结果。

具体实施方式

以下参照附图对本发明的具体实现及效果做进一步的详细说明。

参照图1，本发明的实现步骤如下:

步骤1，输入两幅在不同时相获取的同一地区的大小为N行N列的SAR图像，分别记为I和J，图像的灰度级范围为0～255。

步骤2，对时相1图像I和时相2图像J分别进行均值滤波，然后再进行直方图均衡化，得到两幅滤波均衡化后图像，分别记为E和F。其中，均值滤波的邻域大小为m×m个像素，m的取值范围为3～8，本发明实例中取值为7。

步骤3，对滤波均衡化后图像E进行基于频谱残差的显著区域分割，将得到的连通区域中面积大于A的区域分别标记为EA₁、EA₂、…、EA_r、…、EA_re，其中EA_r表示第r个面积大于A的区域，re表示对滤波均衡化后图像E进行基于频谱残差的显著区域分割后得到的连通区域中面积大于面积阈值A的区域个数，1≤r≤re。对滤波均衡化后图像F进行基于频谱残差的显著区域分割，将得到的连通区域中面积大于A的区域分别标记为FA₁、FA₂、…、FA_t、…、FA_tf，其中FA_t表示连通区域中第t个面积大于A的区域，tf表示对滤波均衡化后图像F进行基于频谱残差的显著区域分割后得到的连通区域中面积大于面积阈值A的区域个数，1≤t≤tf。具体步骤如下：

3a)对滤波均衡化后图像E进行傅里叶变换得到变换后矩阵Eft，对傅里叶变换后的矩阵Eft的幅度图像取以10为底的对数，得到幅度对数图像Et。

3b)将幅度对数图像Et减去该图像的均值滤波后图像，得到滤波差值图像，将滤波差值图像作为合成复数矩阵的实部，将傅里叶变换后的矩阵Eft的虚部作为合成复数矩阵的虚部，对合成复数矩阵进行傅里叶反变换并进行高斯低通滤波，得到显著图G。均值滤波的邻域大小为n×n个像素，n的取值范围为3～8，本发明实例中取值为3，高斯低通滤波器的标准差设为2.5。

3c)对显著图G进行阈值分割，得到二值图像Gb，将二值图像Gb中大于阈值的像素赋值为1，小于阈值的像素赋值为0，若大于阈值的像素所形成区域的总面积比小于阈值的像素所形成区域的总面积大，则将大于阈值的像素的值赋为0，同时将小于阈值的像素的值赋为1，否则不做处理。得到一幅阈值分割图像Ew，其中，阈值的取值范围为0.02～0.05。本发明实例中取值为0.035。

3d)将阈值分割图像Ew中连通区域中面积大于面积阈值A的区域分别标记为EA₁、EA₂、…、EA_r、…、EA_re，其中面积阈值A≥η×N²，η的取值范围为0.01～0.03，本发明实例中取值为0.0125，N²为图像的像素个数，EA_r表示第r个面积大于A的区域，re为连通区域中面积大于面积阈值A的区域个数，1≤r≤re。

3e)对滤波均衡化后图像F进行与步骤3a)到步骤3b)同样地处理，将得到的阈值分割图像的连通区域中面积大于A的区域标记为FA₁、FA₂、…、FA_t、…、FA_tf，其中FA_t表示连通区域中第t个面积大于A的区域，tf表示连通区域中面积大于面积阈值A的区域个数，1≤t≤tf。

步骤4，利用直方图众数排序的方法计算时相1的所有分割区域和时相2的所有分割区域中欧氏距离最小的一对匹配区域，在这对匹配区域中，时相1的分割区域记为EA_mj，对应的图像块记为IA_mj，时相2的分割区域记为FA_mj，对应的图像块记为JA_mj。具体步骤如下：

4a)计算分割区域EA₁的灰度直方图，得到其灰度级像素数目矩阵AN_EA1＝[AN₀，AN₁，AN₂，…，AN_j，…AN₂₅₅]，其中AN_j表示该分割区域的第j个灰度级对应的像素数目，0≤j≤255，该分割区域的灰度级像素数目矩阵所对应的灰度级矩阵为L_EA1＝[0，1，2，…，j，…255]。将灰度级像素数目矩阵AN_EA1进行降序排列，得到灰度级像素数目降序矩阵AN_EA1s，将灰度级矩阵L_EA1按照AN_EA1s中的顺序重新排序，得到排序后灰度级矩阵L_EA1s＝[a₀，a₁，a₂，…，a_j，…，a₂₅₅]，其中a_j表示排序后灰度级矩阵L_EA1s的第j个矩阵元素，为灰度级像素数目降序矩阵AN_EA1s中的第j个矩阵元素所对应的灰度级值，0≤a_j≤255。

4b)对时相1的滤波均衡化后图像E中其余的分割区域EA₂、EA₃、…、EA_r、…、EA_re进行与步骤4a)同样地处理，分别得到对应的排序后灰度级矩阵L_EA2s、L_EA3s、…、L_EArs、…、L_EAres。同样地，对时相2的滤波均衡化后图像F中的所有的分割区域FA₁、FA₂、…、FA_t、…、FA_tf，进行与步骤4a)同样地处理，分别得到对应的排序后灰度级矩阵L_FA1s、L_FA2s、…、L_FAts、…、L_FAtfs。

4c)计算时相1中第j个排序后灰度级矩阵L_EAjs与时相2中各个排序后灰度级矩阵L_FA1s、L_FA2s、…、L_FAts、…、L_FAtfs的欧氏距离，将这些欧氏距离中的最小值所对应的时相2的排序后灰度级矩阵记为L_FAsj，其对应的时相2分割区域记为FA_sj，该最小欧氏距离记为D_sj，其中sj的取值范围为[1，re]。

4d)重复步骤4a)到4c)计算完时相1中所有的排序后灰度级矩阵所对应的最小欧氏距离，得到最小欧氏距离矩阵DS＝[D_s1D_s2…D_sr…D_sre]，最小欧氏距离矩阵DS中的最小值记为Ds_mj，其对应的时相1和时相2的分割区域分别为EA_mj和FA_mj，即为一对匹配区域对，它们对应的图像块分别记为为IA_mj和JA_mj。

步骤5，利用改进后的综合学习粒子群算法，计算时相2图像块JA_mj相对于时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy。具体步骤如下：

5a)初始化粒子群的代数Gen和粒子数目Ps，其中，Gen的取值范围为Gen≥10，本发明实例中取值为10，粒子数目Ps的取值范围为Ps≥150，本发明实例中取值为150。每个粒子在逐代进化过程中所产生的适应度值中的最大值称为该粒子的历史最大适应度值，该粒子的历史最大适应度值对应的位置称为该粒子的历史最优粒子位置，所有粒子在进化过程所能产生的适应度值中的最大值称为最大适应度值，最大适应度值对应的粒子位置称为最优粒子位置。初始化每个粒子历史最大适应度值保持不变的代数阈值为M，M的取值范围为1～3，本发明实例中取值为2。

5b)一个粒子的粒子位置对应为一个三维向量，其第一维参数取值范围为[-90,90]，第二维参数取值范围为[-row,row]，第三维参数取值范围为[-col,col]，row和col分别为时相2图像块JA_mj的行数和列数。同样地，粒子速度也是一个由三个参数构成的三维向量，每一维参数的取值范围为[0,1]。在粒子位置和粒子速度的每一维参数取值范围内随机初始化第一代种群的每个粒子位置和粒子速度的每一维参数。

5c)将第i个粒子的历史最大适应度值保持不变的代数fg_i初始化为0。由第i个粒子的位置pos_i(θ_i,Δx_i,Δy_i)，可建立几何变换矩阵 ${\mathrm{tf}}_i=\left[\begin{array}{ccc}\cos \left({\mathrm{\θ}}_i\right)& \sin \left({\mathrm{\θ}}_i\right)& {\mathrm{\Δx}}_i\\ -\sin \left({\mathrm{\θ}}_i\right)& \cos \left({\mathrm{\θ}}_i\right)& {\mathrm{\Δy}}_i\\ 0& 0& 1\end{array}\right].$ 对时相2图像块JA_mj根据几何变换矩阵tf_i进行几何变换，得到时相2图像块的变换后图像块JT_i，进行几何变换时所用到的插值方法为双线性插值方法。

5d)计算第i个粒子的适应度值公式如下：

f(θ_i,Δx_i,Δy_i)＝-H(RI)

其中，图像RI＝{RI(x,y)|1≤x≤col,1≤y≤row}，其中，像素(x,y)的灰度值RI(x,y)计算公式如下：

其中，IA_mj(x,y)和JT_i(x,y)分别表示图像块IA_mj和JT_i中坐标位于(x,y)处的像素值。max(·)是取最大值函数。rand_(0-255)表示取值范围为[0,255]的随机数。H(RI)为计算图像RI的熵。

5e)将第一代种群粒子中适应度值最大的粒子的位置记为最优粒子位置Gb，其适应度值记为最大适应度值FG。将第i个粒子的历史最优位置Pb_i初始化为其位置,每个粒子的历史最大适应度值FP_i初始化为其适应度值。

5f)对于第一代种群的第i个粒子，i的取值范围为1,2,…Ps，如果其历史最大适应度值保持不变的代数fg_i大于每个粒子历史最大适应度值保持不变的代数阈值M，则按普通粒子群学习策略更新第i个粒子的位置，具体公式如下：

$X_i^d\←X_i^d+V_i^d+c_1\×\mathrm{rand}{1}_i^d\×({\mathrm{Pb}}_i^d-X_i^d)+c_2\×{\mathrm{rand}2}_i^d\×({\mathrm{Gb}}_i^d-X_i^d)$

其中，符号←左边的X_i ^d表示为更新后的第i个粒子位置的第d维参数，符号←右边的X_i ^d表示为更新前的第i个粒子位置的第d维参数，符号←表示替代的更新操作，V_i ^d为第i个粒子速度的第d维参数。Pb_i ^d为第i个粒子历史最优位置Pb_i的第d维参数。Gb_i ^d为最优粒子位置Gb的第d维参数。c₁和c₂为学习因子，取值为2。rand1_i ^d和rand2_i ^d为随机数，取值范围都为[01]。

如果第i个粒子的历史最大适应度值保持不变的代数fg_i大于Mn，则除做上述位置更新之外，执行步骤5g)的克隆选择操作，其中Mn≥M，本发明实例中Mn＝M+2。

5g)将第i个粒子的粒子位置X_i克隆GN_i次，形成GN_i个和粒子位置X_i相同的克隆个体。其中，GN_i≥50，本发明实例中GN_i取值为50，设置变异概率Pm＝0.1，克隆个体的序号数ii的范围1≤ii≤GN_i。

5h)对于任意的第ii个克隆个体Gb_ii按如下公式进行变异：

${\mathrm{Gb}}_{\mathrm{ii}}\←\left\{\begin{array}{cc}{\mathrm{Gb}}_{\mathrm{iin}}& {\mathrm{rand}}_{(0-1)}>\mathrm{Pm}\\ {\mathrm{Gb}}_{\mathrm{ii}}& {\mathrm{rand}}_{(0-1)}\≤\mathrm{Pm}\end{array}\right.$

其中，Gb_iin是一个1×3的矩阵，其第d个矩阵元素取为随机数rand_(m)，rand_(m)的取值范围为[min^d,max^d]，min^d和max^d分别表示第d维参数取值范围的下界和上界，rand_(0-1)表示取随机数，取值范围为[0,1]。

5i)将每一个变异后的克隆个体视为一个粒子，按照步骤5d)中的粒子适应度值公式计算其适应度值，将所有克隆个体适应度值中的最大值所对应的克隆个体记为Gb_ic。如果克隆个体Gb_ic的适应度值大于第i个粒子的适应度值，则将第i粒子的粒子位置X_i更新为克隆个体Gb_ic的粒子位置。

对更新后的第i个粒子按照步骤5d)中的粒子适应度值公式计算其适应度值，如果该适应度值大于最大适应度值FG，则将最大适应度值FG更新为该适应度值，最优粒子位置Gb更新为该粒子的粒子位置。如果该适应度值大于第i个粒子的历史最大适应度值，则将第i个粒子的历史最大适应度值FP_i更新为该适应度值，第i个粒子的历史最优粒子位置Pb_i更新为该粒子的粒子位置。

5j)运用综合学习策略更新第i个粒子的粒子位置，具体公式如下：

$X_i^d\←X_i^d+V_i^d$

其中，符号←左边的X_i ^d表示为更新后的第i个粒子位置的第d维参数，符号←右边的X_i ^d表示为更新前的第i个粒子位置的第d维参数，符号←表示替代的更新操作，V_i ^d表示为第i个粒子速度的第d维参数，V_i ^d的更新公式如下：

${V_i}^d\&LeftArrow;\left\{\begin{array}{cc}{W}_k\&times;{V_i}^d+c\&times;{\mathrm{rand}}_i^d\&times;({\mathrm{Pb}}_{f_i}^d-X_i^d)& \mathrm{ran}<{\mathrm{Pc}}_i\\ W_k\&times;{V_i}^d+c\&times;{\mathrm{rand}}_i^d\&times;({\mathrm{Pb}}_i^d-X_i^d)& \mathrm{ran}\&GreaterEqual;{\mathrm{Pc}}_i\end{array}\right.$

其中，ran为随机数，取值范围为[0,l]。符号←左边的V_i ^d表示为更新后的第i个粒子速度的第d维参数，符号←右边的V_i ^d表示为更新前的第i个粒子速度的第d维参数，惯性权重w₀＝0.9，w₁＝0.4，k为当前的种群代数，学习因子c取值为1.49445，rand_i ^d为随机数，取值范围为[0 1]，为第f_i个粒子历史最优位置的第d维参数。

决策概率Pc_i的计算公式如下:

${\mathrm{Pc}}_i=0.05+0.45\&times;\frac{(e^{\frac{10(i-1)}{\mathrm{ps}-1}}-1)}{(e^{10}-1)}.$

f_i的计算公式如下:

$f_i=\left\{\begin{array}{cc}{f}_{1i}& {\mathrm{FP}}_{f_{1i}}>{\mathrm{FP}}_{f_{2i}}\\ f_{2i}& {\mathrm{FP}}_{f_{1i}}\&le;{\mathrm{FP}}_{f_{2i}}\end{array}\right.$

其中，f_1i和f_2i为速度未更新的所有粒子中随机选取的两个粒子。

对更新后的第i个粒子按照步骤5d)中的粒子适应度值公式计算其适应度值，如果该适应度值大于最大适应度值FG，则将最大适应度值FG更新为该适应度值，最优粒子位置Gb更新为该粒子的粒子位置。如果该适应度值大于第i个粒子的历史最大适应度值，则将第i个粒子的历史最大适应度值FP_i更新为该适应度值，第i个粒子的历史最优粒子位置Pb_i更新为该粒子的粒子位置，并将该粒子的历史最大适应度值保持不变的代数fg_i置0，否则将其加1。

5k)重复上述步骤5f)到5j)，将计算完所有粒子后得到的最优粒子位置Gb克隆GN次，形成GN个和最优粒子位置Gb相同的克隆个体，GN≥200，本发明实例中取值为200。设置变异概率Pm＝0.1。克隆个体的序号数ii的范围1≤ii≤GN。执行步骤5h)。然后，将每一个变异后的克隆个体视为一个粒子，按照步骤5d)中的粒子适应度值公式计算其适应度值，所有克隆个体适应度值中的最大值所对应的克隆个体为Gbc，如果该克隆个体Gbc的适应度值大于最大适应度值FG，则将最大适应度值FG更新为该适应度值，最优粒子位置Gb更新为克隆个体Gbc。

5l)重复上述5f)到5k)的过程，依次计算第2、3、…直至第Gen代种群。最优粒子位置Gb的第一维参数、第二维参数、第三维参数分别对应时相2图像块JA_mj相对时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy。

步骤6，根据时相2图像块JA_mj相对于时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx、Y轴方向平移参数Δy计算时相2图像J相对于时相1图像I的几何变换矩阵TF。具体步骤如下：

6a)根据时相2图像块JA_mj相对时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy，构建时相2图像块相对时相1图像块的变换矩阵tf：

$\mathrm{tf}=\left[\begin{array}{ccc}\cos \left(\mathrm{\&theta;}\right)& \sin \left(\mathrm{\&theta;}\right)& \mathrm{\&Delta;x}\\ -\sin \left(\mathrm{\&theta;}\right)& \cos \left(\mathrm{\&theta;}\right)& \mathrm{\&Delta;y}\\ 0& 0& 1\end{array}\right]$

6b)从时相2图像块JA_mj选出3个点JA_p1(JA_p1x,JA_p1y)、JA_p2(JA_p2x,JA_p2y)、JA_p3(JA_p3x,JA_p3y)，其坐标位置按如下方式计算：

JA_p1x＝0.25×JAc，JA_p1y＝0.5×JAr

JA_p2x＝0.5×JAc，JA_p2y＝0.25×JAr

JA_p3x＝0.75×JAc，JA_p3y＝0.5×JAr

其中JAr和JAc分别为时相2图像块JA_mj的行数和列数。

根据变换矩阵tf可得到时相1图像块IA_mj中对应的3个点IA_p1(IA_p1x,IA_p1y)、IA_p2(IA_p2x,,IA_p2y)、IA_p3(IA_p3x,IA_p3y)，计算方式如下：

$\left[\begin{array}{ccc}{\mathrm{IA}}_{p1x}& {\mathrm{IA}}_{p2x}& {\mathrm{IA}}_{p3x}\\ {\mathrm{IA}}_{p1y}& {\mathrm{IA}}_{p2y}& {\mathrm{IA}}_{p3y}\\ 1& 1& 1\end{array}\right]=\mathrm{tf}\&times;\left[\begin{array}{ccc}{\mathrm{JA}}_{p1x}& {\mathrm{JA}}_{p2x}& {\mathrm{JA}}_{p3x}\\ {\mathrm{JA}}_{p1y}& {\mathrm{JA}}_{p2y}& {\mathrm{JA}}_{p3y}\\ 1& 1& 1\end{array}\right]$

6c)计算时相1图像块IA_mj中点IA_p1(IA_p1x,IA_p1y)、IA_p2(IA_p2x,,IA_p2y)、IA_p3(IA_p3x,IA_p3y)在时相1图像I中对应的点I_p1(I_p1x,I_p1y)、I_p2(I_p2x,,I_p2y)、I_p3(I_p3x,I_p3y)，计算时相2图像块JA_mj中点JA_p1(JA_p1x,JA_p1y)、JA_p2(JA_p2x,,JA_p2y)、JA_p3(JA_p3x,JA_p3y)在时相2图像J中对应的点J_p1(J_p1x,J_p1y)、J_p2(J_p2x,J_p2y)、J_p3(J_p3x,J_p3y)，计算公式如下：

IA_p1x＝IA_p1x+Ix0-1,IA_p1y＝IA_p1y+Iy0-1

IA_p2x＝IA_p2x+Ix0-1,IA_p2y＝IA_p2y+Iy0-1

IA_p3x＝IA_p3x+Ix0-1,IA_p3y＝IA_p3y+Iy0-1

JA_p1x＝JA_p1x+Jx0-1,JA_p1y＝JA_p1y+Jy0-1

JA_p2x＝JA_p2x+Jx0-1,JA_p2y＝JA_p2y+Jy0-1

JA_p3x＝JA_p3x+Jx0-1,JA_p3y＝JA_p3y+Jy0-1

其中(Ix0，Iy0)是时相1图像块IA_mj中坐标位置(1,1)的像素点在时相1图像I中的坐标。(Jx0，Jy0)是时相2图像块JA_mj中坐标位置(1,1)的像素点在时相2图像J中的坐标。

6d)根据时相1图像I中的点I_p1(I_p1x,I_p1y)、I_p2(I_p2x,,I_p2y)、I_p3(I_p3x,I_p3y)和时相2图像J中的点J_p1(J_p1x,J_p1y)、J_p2(J_p2x,J_p2y)、J_p3(J_p3x,J_p3y)计算时相2图像J相对时相1图像I的几何变换矩阵TF，计算公式如下：

$\mathrm{TF}=\frac{\left[\begin{array}{ccc}{I}_{p1x}& I_{p2x}& I_{p3x}\\ I_{p1y}& I_{p2y}& I_{p3y}\\ 1& 1& 1\end{array}\right]}{\left[\begin{array}{ccc}{J}_{p1x}& J_{p2x}& J_{p3x}\\ J_{p1y}& J_{p2y}& J_{p3y}\\ 1& 1& 1\end{array}\right]}.$

步骤7，根据几何变换矩阵TF，对时相2图像进行几何变换，其中的插值方法采用双线性插值方法。

本发明的效果可以进一步通过以下内容进行说明：

1.图1为本发明的框架流程图。

2.图2(a1)、图2(b1)是2008年RadarSat-2卫星获取的黄河口图像的部分图像。图2(a2)、图2(b2)是2009年RadarSat-2卫星获取的黄河口图像的部分图像。

3.改进后的综合学习粒子群算法CSA-CLPSO和原综合学习粒子群算法CLPSO进行对比。其中，种群规模都为100，进化代数设置为10，运行次数都为100次，统计100次运行结果的最优适应度值和平均适应度值，对比结果如下：

a)F₁＝21.5+x₁×sin(4π×x₁)+x₂×sin(20π×x₂)

表1F₁的进化结果对比

	CLPSO	CSA-CLPSO
			最优适应度值对应的x1	11.621978	11.625059
最优适应度值对应的x2	5.725319	5.725235
			最优适应度值	38.837766	38.849668
平均适应度值	38.140903	38.662368

b)F₂＝(x₂+x₁)⁴+(x₂+x₁)²+x₂+x₁

表2F₂的进化结果对比

	CLPSO	CSA-CLPSO
			最优适应度值对应的x1	12.098836	12.099771
最优适应度值对应的x2	5.796861	5.798752
			最优适应度值	102902.048796	102966.937106
平均适应度值	102041.619751	102573.613205

c)F₃＝(x₁-x₂)³+(x₁-x₂)²+x₁-x₂

表1F₃的进化结果对比

	CLPSO	CSA-CLPSO
			最优适应度值对应的x1	12.097875	12.099767
最优适应度值对应的x2	4.100414	4.100057
			最优适应度值	583.469628	583.939270
平均适应度值	574.245806	579.535028

4.图3(a1)、图3(a2)、图3(b1)、图3(b2)分别是图2(a1)、图2(a2)、图2(b1)、图2(b2)基于频谱残差的区域分割结果。

5.图4(a1)、图4(a2)是通过步骤4得到的图3(a1)、图3(a2)中最相似的一对匹配区域。图4(b1)、图4(b2)是通过步骤4得到的图3(b1)、图3(b2)中最相似的一对匹配区域。

6.图5(a1)、图5(a3)即为图4(a1)、图4(a2)，图5(b1)、图5(b3)即为图4(b1)、图4(b2)。图5(a2)为经步骤5计算后，图5(a3)几何变换后的结果。图5(b2)为经步骤5计算后，图5(b3)几何变换后的结果。

7.图6(a1)、图6(a3)即为图2(a1)、图2(a2)，图6(b1)、图6(b3)即为图2(b1)、图2(b2)。图6(a2)为经步骤6、步骤7计算后，图6(a3)几何变换后的结果。图6(b2)为经步骤6、步骤7计算后，图6(b3)几何变换后的结果。

8.图7(a)是图2(a1)、图2(a2)根据配准参数的拼接结果，图7(b)是图2(b1)、图2(b2)根据配准参数的拼接结果。

总结分析：

从表1、表2、表3中的数据可以看出，改进后的综合学习粒子群算法在相同参数设置的情况下表现出更好的寻优效果。对图2(a1)和图2(a2)的配准中间结果和最终结果可以看出，本方法有效克服了变化区域对配准的干扰。从图3(b1)、图3(b2)，图4(b1)、图4(b2)，图5(b1)、图5(b2)、图5(b3)可以看出本方法在区域分割一致性不是很好的情况下实现了区域匹配，并对匹配区域实现良好的配准。图2到图7显示的配准中间结果可验证本方法流程的可行性，从最终的配准结果图6以及图像拼接结果图7可以看出本方法的有效性，实现对两幅SAR图像的良好配准。

本实施方式中没有详细叙述的部分属本行业的公知的常用手段，这里不一一叙述。以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (5)

1.一种基于显著分割子区域对的SAR图像配准方法，其特征在于，包括如下步骤：

步骤1，输入两幅在不同时相获取的同一地区的大小为N行N列的SAR图像，分别记为I和J，图像的灰度级范围为0～255；

步骤2，对时相1图像I和时相2图像J分别进行均值滤波，然后再进行直方图均衡化，得到两幅滤波均衡化后图像，分别记为E和F，其中，均值滤波的邻域大小为m×m个像素，m的取值范围为3～8；

步骤3，对滤波均衡化后图像E进行基于频谱残差的显著区域分割，将得到的连通区域中面积大于A的区域分别标记为EA₁、EA₂、…、EA_r、…、EA_re，其中EA_r表示第r个面积大于A的区域，re表示对滤波均衡化后图像E进行基于频谱残差的显著区域分割后得到的连通区域中面积大于面积阈值A的区域个数，1≤r≤re；对滤波均衡化后图像F进行基于频谱残差的显著区域分割，将得到的连通区域中面积大于A的区域分别标记为FA₁、FA₂、…、FA_t、…、FA_tf，其中FA_t表示连通区域中第t个面积大于A的区域，tf表示对滤波均衡化后图像F进行基于频谱残差的显著区域分割后得到的连通区域中面积大于面积阈值A的区域个数，1≤t≤tf；

步骤4，利用直方图众数排序的方法计算时相1的所有分割区域和时相2的所有分割区域中欧氏距离最小的一对匹配区域，在这对匹配区域中，时相1的分割区域记为EA_mj，对应的图像块记为IA_mj，时相2的分割区域记为FA_mj，对应的图像块记为JA_mj；

步骤5，利用改进后的综合学习粒子群算法，计算时相2图像块JA_mj相对于时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy；

步骤6，根据时相2图像块JA_mj相对于时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx、Y轴方向平移参数Δy计算时相2图像J相对于时相1图像I的几何变换矩阵TF；

步骤7，根据几何变换矩阵TF，对时相2图像进行几何变换，其中的插值方法采用双线性插值方法。

2.根据权利要求1所述的一种基于显著分割子区域对的SAR图像配准方法，其特征在于，步骤3具体按如下步骤进行:

3a)对滤波均衡化后图像E根据基于频谱残差的显著图计算方法计算其显著图G；

3b)对显著图G进行阈值分割，得到二值图像Gb，将二值图像Gb中大于阈值的像素赋值为1，小于阈值的像素赋值为0，若大于阈值的像素所形成区域的总面积比小于阈值的像素所形成区域的总面积大，则将大于阈值的像素的值赋为0，同时将小于阈值的像素的值赋为1，否则不做处理；得到一幅阈值分割图像Ew，其中，阈值的取值范围为0.02～0.05；

3c)将阈值分割图像Ew中连通区域中面积大于面积阈值A的区域分别标记为EA₁、EA₂、…、EA_r、…、EA_re，其中面积阈值A≥η×N²，η的取值范围为0.01～0.03，N²为图像的像素个数，EA_r表示第r个面积大于A的区域，re为连通区域中面积大于面积阈值A的区域个数，1≤r≤re；

3d)对滤波均衡化后图像F进行与步骤3a)到步骤3b)同样地处理，将得到的阈值分割图像的连通区域中面积大于A的区域标记为FA₁、FA₂、…、FA_t、…、FA_tf，其中FA_t表示连通区域中第t个面积大于A的区域，tf表示连通区域中面积大于面积阈值A的区域个数，1≤t≤tf。

3.根据权利要求1所述的一种基于显著分割子区域对的SAR图像配准方法，其特征在于，其中步骤4具体按如下步骤进行：

4a)计算分割区域EA₁的灰度直方图，得到其灰度级像素数目矩阵AN_EA1＝[AN₀，AN₁，AN₂，…，AN_j，…AN₂₅₅]，其中AN_j表示该分割区域的第j个灰度级对应的像素数目，0≤j≤255，该分割区域的灰度级像素数目矩阵所对应的灰度级矩阵为L_EA1＝[0，1，2，…，j，…255]；将灰度级像素数目矩阵AN_EA1进行降序排列，得到灰度级像素数目降序矩阵AN_EA1s，将灰度级矩阵L_EA1按照AN_EA1s中的顺序重新排序，得到排序后灰度级矩阵L_EA1s＝[a₀，a₁，a₂，…，a_j，…，a₂₅₅]，其中a_j表示排序后灰度级矩阵L_EA1s的第j个矩阵元素，为灰度级像素数目降序矩阵AN_EA1s中的第j个矩阵元素所对应的灰度级值，0≤a_j≤255；

4b)对时相1的滤波均衡化后图像E中其余的分割区域EA₂、EA₃、…、EA_r、…、EA_re进行与步骤4a)同样地处理，分别得到对应的排序后灰度级矩阵L_EA2s、L_EA3s、…、L_EArs、…、L_EAres；同样地，对时相2的滤波均衡化后图像F中的所有的分割区域FA₁、FA₂、…、FA_t、…、FA_tf，进行与步骤4a)同样地处理，分别得到对应的排序后灰度级矩阵L_FA1s、L_FA2s、…、L_FAts、…、L_FAtfs；

4c)计算时相1中第j个排序后灰度级矩阵L_EAjs与时相2中各个排序后灰度级矩阵L_FA1s、L_FA2s、…、L_FAts、…、L_FAtfs的欧氏距离，将这些欧氏距离中的最小值所对应的时相2的排序后灰度级矩阵记为L_FAsj，其对应的时相2分割区域记为FA_sj，该最小欧氏距离记为D_sj，其中sj的取值范围为[1，re]；

4d)重复步骤4a)到4c)计算完时相1中所有的排序后灰度级矩阵所对应的最小欧氏距离，得到最小欧氏距离矩阵DS＝[D_s1D_s2…D_sr…D_sre]，最小欧氏距离矩阵DS中的最小值记为Ds_mj，其对应的时相1和时相2的分割区域分别为EA_mj和FA_mj，即为一对匹配区域对，它们对应的图像块分别记为IA_mj和JA_mj。

4.根据权利要求1所述的一种基于显著分割子区域对的SAR图像配准方法，其特征在于，其中步骤5具体按如下步骤进行：

5a)初始化粒子群的代数Gen和粒子数目Ps，其中，Gen的取值范围为Gen≥10，粒子数目Ps的取值范围为Ps≥150；每个粒子在逐代进化过程中所产生的适应度值中的最大值称为该粒子的历史最大适应度值，该粒子的历史最大适应度值对应的位置称为该粒子的历史最优粒子位置，所有粒子在进化过程所能产生的适应度值中的最大值称为最大适应度值，最大适应度值对应的粒子位置称为最优粒子位置；初始化每个粒子历史最大适应度值保持不变的代数阈值为M，M的取值范围为1～3；

5b)一个粒子的粒子位置为一个三维向量，其第一维参数取值范围为[-90,90]，第二维参数取值范围为[-row,row]，第三维参数取值范围为[-col,col]，row和col分别为时相2图像块JA_mj的行数和列数；同样地，粒子速度也是一个由三个参数构成的三维向量，每一维参数的取值范围均为[0,1]；在粒子位置和粒子速度的每一维参数取值范围内随机初始化第一代种群的每个粒子位置和粒子速度中的每一维参数；

5c)将第i个粒子的历史最大适应度值保持不变的代数fg_i初始化为0；由第i个粒子的位置pos_i(θ_i,Δx_i,Δy_i)，建立几何变换矩阵

${\mathrm{tf}}_i=\left[\begin{array}{ccc}\cos \left({\mathrm{\&theta;}}_i\right)& \sin \left({\mathrm{\&theta;}}_i\right)& {\mathrm{\&Delta;x}}_i\\ -\sin \left({\mathrm{\&theta;}}_i\right)& \cos \left({\mathrm{\&theta;}}_i\right)& {\mathrm{\&Delta;y}}_i\\ 0& 0& 1\end{array}\right];$ 对时相2图像块JA_mj根据几何变换矩阵tf_i进行几何变换，得到时相2图像块的变换后图像块JT_i，进行几何变换所用的插值方法为双线性插值方法；

5d)计算第i个粒子的适应度值公式如下：

f(θ_i,Δx_i,Δy_i)＝-H(RI)

其中，图像RI＝{RI(x,y)|1≤x≤col,1≤y≤row}，其中，像素(x,y)的灰度值RI(x,y)计算公式如下：

其中，IA_mj(x,y)和JT_i(x,y)分别表示图像块IA_mj和JT_i中坐标位于(x,y)处的像素值；max(·)是取最大值函数；rand_(0-255)表示取值范围为[0,255]的随机数；H(RI)为计算图像RI的熵；

5e)将第一代种群粒子中适应度值最大的粒子的位置记为最优粒子位置Gb，其适应度值记为最大适应度值FG；将第i个粒子的历史最优位置Pb_i初始化为其位置,每个粒子的历史最大适应度值FP_i初始化为其适应度值；

5f)对于第一代种群的第i个粒子，i的取值范围为1,2,…Ps，如果其历史最大适应度值保持不变的代数fg_i大于每个粒子历史最大适应度值保持不变的代数阈值M，则按普通粒子群学习策略更新第i个粒子的位置；如果第i个粒子的历史最大适应度值保持不变的代数fg_i大于Mn，则除做上述位置更新之外，执行步骤5g)的克隆选择操作，其中Mn≥M；

5g)将第i个粒子的粒子位置X_i克隆GN_i次，形成GN_i个和粒子位置X_i相同的克隆个体；其中，GN_i≥50，设置变异概率Pm＝0.1，克隆个体的序号数ii的范围1≤ii≤GN_i；

5h)对于任意的第ii个克隆个体Gb_ii按如下公式进行变异：

${\mathrm{Gb}}_{\mathrm{ii}}\&LeftArrow;\left\{\begin{array}{cc}{\mathrm{Gb}}_{\mathrm{iin}}& {\mathrm{rand}}_{(0-1)}>\mathrm{Pm}\\ {\mathrm{Gb}}_{\mathrm{ii}}& {\mathrm{rand}}_{(0-1)}\&le;\mathrm{Pm}\end{array}\right.$

其中，Gb_iin是一个1×3的矩阵，其第d个矩阵元素取为随机数rand_(m)，rand_(m)的取值范围为[min^d,max^d]，min^d和max^d分别表示第d维参数取值范围的下界和上界，rand_(0-1)表示取随机数，取值范围为[0,1]；

5i)将每一个变异后的克隆个体视为一个粒子，按照步骤5d)中的粒子适应度值公式计算其适应度值，将所有克隆个体适应度值中的最大值所对应的克隆个体记为Gb_ic；如果克隆个体Gb_ic的适应度值大于第i个粒子的适应度值，则将第i粒子的粒子位置X_i更新为克隆个体Gb_ic的粒子位置；

对更新后的第i个粒子按照步骤5d)中的粒子适应度值公式计算其适应度值，如果该适应度值大于最大适应度值FG，则将最大适应度值FG更新为该适应度值，最优粒子位置Gb更新为该粒子的粒子位置；如果该适应度值大于第i个粒子的历史最大适应度值，则将第i个粒子的历史最大适应度值FP_i更新为该适应度值，第i个粒子的历史最优粒子位置Pb_i更新为该粒子的粒子位置；

5j)运用综合学习策略更新第i个粒子的粒子位置；

5k)重复上述步骤5f)到5j)，将计算完所有粒子后得到的最优粒子位置Gb克隆GN次，形成GN个和最优粒子位置Gb相同的克隆个体，GN≥200，设置变异概率Pm＝0.1，克隆个体的序号数ii的范围1≤ii≤GN；执行步骤5h)；然后，将每一个变异后的克隆个体视为一个粒子，按照步骤5d)中的粒子适应度值公式计算其适应度值，所有克隆个体适应度值中的最大值所对应的克隆个体为Gbc，如果该克隆个体Gbc的适应度值大于最大适应度值FG，则将最大适应度值FG更新为该适应度值，最优粒子位置Gb更新为克隆个体Gbc；

5l)重复上述5f)到5k)的过程，依次计算第2、3、…直至第Gen代种群；最优粒子位置Gb的第一维参数、第二维参数、第三维参数分别对应时相2图像块JA_mj相对时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy。

5.根据权利要求1所述的一种基于显著分割子区域对的SAR图像配准方法，其特征在于，其中步骤6具体按如下步骤进行：

6a)根据时相2图像块JA_mj相对时相1图像块IA_mj的旋转角度θ、X轴方向平移参数Δx和Y轴方向平移参数Δy，构建时相2图像块相对时相1图像块的变换矩阵tf：

$\mathrm{tf}=\left[\begin{array}{ccc}\cos \left(\mathrm{\&theta;}\right)& \sin \left(\mathrm{\&theta;}\right)& \mathrm{\&Delta;x}\\ -\sin \left(\mathrm{\&theta;}\right)& \cos \left(\mathrm{\&theta;}\right)& \mathrm{\&Delta;y}\\ 0& 0& 1\end{array}\right];$

6b)从时相2图像块JA_mj选出3个点JA_p1(JA_p1x,JA_p1y)、JA_p2(JA_p2x,JA_p2y)、JA_p3(JA_p3x,JA_p3y)，其坐标位置按如下方式计算：

JA_p1x＝0.25×JAc，JA_p1y＝0.5×JAr

JA_p2x＝0.5×JAc，JA_p2y＝0.25×JAr

JA_p3x＝0.75×JAc，JA_p3y＝0.5×JAr

其中JAr和JAc分别为时相2图像块JA_mj的行数和列数；

根据变换矩阵tf可得到时相1图像块IA_mj中对应的3个点IA_p1(IA_p1x,IA_p1y)、IA_p2(IA_p2x,,IA_p2y)、IA_p3(IA_p3x,IA_p3y)，计算方式如下：

$\left[\begin{array}{ccc}{\mathrm{IA}}_{p1x}& {\mathrm{IA}}_{p2x}& {\mathrm{IA}}_{p3x}\\ {\mathrm{IA}}_{p1y}& {\mathrm{IA}}_{p2y}& {\mathrm{IA}}_{p3y}\\ 1& 1& 1\end{array}\right]=\mathrm{tf}\&times;\left[\begin{array}{ccc}{\mathrm{JA}}_{p1x}& {\mathrm{JA}}_{p2x}& {\mathrm{JA}}_{p3x}\\ {\mathrm{JA}}_{p1y}& {\mathrm{JA}}_{p2y}& {\mathrm{JA}}_{p3y}\\ 1& 1& 1\end{array}\right];$

6c)计算时相1图像块IA_mj中点IA_p1(IA_p1x,IA_p1y)、IA_p2(IA_p2x,,IA_p2y)、IA_p3(IA_p3x,IA_p3y)在时相1图像I中对应的点I_p1(I_p1x,I_p1y)、I_p2(I_p2x,,I_p2y)、I_p3(I_p3x,I_p3y)，计算时相2图像块JA_mj中点JA_p1(JA_p1x,JA_p1y)、JA_p2(JA_p2x,,JA_p2y)、JA_p3(JA_p3x,JA_p3y)在时相2图像J中对应的点J_p1(J_p1x,J_p1y)、J_p2(J_p2x,J_p2y)、J_p3(J_p3x,J_p3y)，计算公式如下：

IA_p1x＝IA_p1x+Ix0-1,IA_p1y＝IA_p1y+Iy0-1

IA_p2x＝IA_p2x+Ix0-1,IA_p2y＝IA_p2y+Iy0-1

IA_p3x＝IA_p3x+Ix0-1,IA_p3y＝IA_p3y+Iy0-1

JA_p1x＝JA_p1x+Jx0-1,JA_p1y＝JA_p1y+Jy0-1

JA_p2x＝JA_p2x+Jx0-1,JA_p2y＝JA_p2y+Jy0-1

JA_p3x＝JA_p3x+Jx0-1,JA_p3y＝JA_p3y+Jy0-1

其中(Ix0，Iy0)是时相1图像块IA_mj中坐标位置(1,1)的像素点在时相1图像I中的坐标；(Jx0，Jy0)是时相2图像块JA_mj中坐标位置(1,1)的像素点在时相2图像J中的坐标；

6d)根据时相1图像I中的点I_p1(I_p1x,I_p1y)、I_p2(I_p2x,,I_p2y)、I_p3(I_p3x,I_p3y)和时相2图像J中的点J_p1(J_p1x,J_p1y)、J_p2(J_p2x,J_p2y)、J_p3(J_p3x,J_p3y)计算时相2图像J相对时相1图像I的几何变换矩阵TF，计算公式如下：

$\mathrm{TF}=\frac{\left[\begin{array}{ccc}{I}_{p1x}& I_{p2x}& I_{p3x}\\ I_{p1y}& I_{p2y}& I_{p3y}\\ 1& 1& 1\end{array}\right]}{\left[\begin{array}{ccc}{J}_{p1x}& J_{p2x}& J_{p3x}\\ J_{p1y}& J_{p2y}& J_{p3y}\\ 1& 1& 1\end{array}\right]}.$