CN104517286A - 基于自适应阈值分割和组合优化的sar图像配准 - Google Patents

Abstract

本发明提出了一种基于自适应阈值分割和组合优化的SAR图像配准方法，以便降低灰度变化差异的影响，采用一种由粗到精的优化寻优策略，从而提高图像配准精度，采用图像小波分解策略，降低Powell局部寻优的运算时间复杂度。其实现过程是：(1)首先对同一地区的两幅图像做自适应阈值分割；(2)通过PSO搜索算法搜索配准参数使得两幅图像之间的互信息达到最大；(3)对图像做三层小波分解，取每层低频分量，以PSO搜索得到的配准参数为Powell搜索算法的初始点，从最高层第三层开始逐层搜索直至原图，得到最后的配准参数；(4)通过得到的配准参数将待配准图像通过插值变换回去，并叠加得到最终的配准结果。

Description

基于自适应阈值分割和组合优化的SAR图像配准

技术领域

本发明属于遥感图像处理领域，涉及SAR图像的配准问题，具体涉及一种基于自适应阈值分割和组合优化的SAR图像配准方法，可用于对SAR图像的配准。

背景技术

合成孔径雷达(synthetic aperture radar,SAR)具有分辨率高，全天时、全天候工作的特点，与可见光、红外传感器比较具有独特的优势和无法替代的作用，被广泛应用于工农业生产、科研和军事等领域。图像配准的过程就是将在不同时间、不同传感器、不同视角及不同成像条件下获取的同一场景下的多幅图像进行匹配和叠加的过程。从方法的角度对图像配准进行分类，可以将图像的配准技术分为三类：基于灰度的配准、基于变换域的配准以及基于特征的配准。

基于灰度的配准有一个前提，那就是两幅图像的灰度必须要相近，不能相差太大，所以一般用于同源图像的配准，不太适合于异源图像的配准。这类方法一般是利用参考图和待配图之间的灰度统计信息作为相似度评价标准，然后利用优化算法寻优找到最优的变换参数。它的一个很明显的优势在于不需要复杂的前期预处理过程，也没有复杂的特征提取过程，而且精度高，鲁棒性好。但是它的缺点也很明显，对灰度变化敏感，没有充分利用灰度统计特性，对每一点的灰度依赖大，同时对光照变化也很敏感，对非线性形变不适用。还有一点，由于采用的是优化算法寻优，运算的时间复杂度很高。

基于变换域的配准是基于一些快速算法实现的，比如fft快速算法、小波分解等等。所以它的优势非常明显，那就是高速，易于硬件实现，而且抗低频噪声强，能做到实时的配准。经过一系列的实验验证，这类方法要求两幅图像之间的重叠区域要大，而且它只能处理平移、旋转、尺度变换，无法处理非线性形变，通常碰到高频噪声配准效果也不太好。

基于特征的图像配准实际上只利用少量的特征信息来精确的定位目标信息，实现图像的配准。其实，直接从视觉感官的角度来讲，我们能直接看出两幅是否是同一场景，这是基于我们对目标相似性的一个判断。比如两幅图像里的某些区域是一致的，某些曲线是对应的等等。所以我们从最直观的角度出发，寻找我们需要的就像肉眼看到的一些图像里的显著的特征。由于特征的多样性，基于特征的配准方法也五花八门。由于我们只采用了图像里的部分特征信息，这就减小了我们的计算量。同时特征有很多优点，它抗噪性好，对灰度变化也不敏感，而且它对图像形变以及光照不均都有很好的适应力(以点特征为代表)。缺点也显而易见，特征提取比较难，特别是对于异源图像。

发明内容

本发明的目的是提供一种基于自适应阈值分割和组合优化的SAR图像配准方法。

本发明的技术方案是，基于自适应阈值分割和组合优化的SAR图像配准，包括如下步骤：

步骤101：导入两幅同一地区的待配准的SAR图像；

步骤102：采用自适应阈值分割提取出最佳阈值，将图像中的显著区域与其它区域分离开来；

步骤103：将分割出的目标区域的灰度归一化到同一灰度级，从0～255归一化到0～25；

步骤104：对归一化的区域图像采用PSO(粒子群算法)对两幅图像之间的互信息进行寻优，搜索到使互信息达到最大的配准参数；

步骤105：对归一化的区域图像做三层小波分解，只取低频分量用于Powell寻优，寻优策略还是搜索配准参数使两幅图像之间的互信息达到最大，由于小波分解改变了图像的大小，我们搜索的空间也将变小，这就会降低我们的寻优时间；

步骤106：假定旋转角度为θ，水平平移为Δx，垂直平移为Δy，尺度缩放为α，将PSO搜索得到的配准参数(θ,Δx,Δy,α)变为(θ,Δx/8,Δy/8,α)，作为Powell算法在小波分解第3层搜索的初始点，搜索配准参数(θ₁,Δx₁,Δy₁,α₁)；

步骤107：因为每层之间图像大小都是两倍关系，将得到的(θ₁,Δx₁,Δy₁,α₁)变换为(θ₁,2Δx₁,2Δy₁,α₁)作为下一层搜索的初始点，继续搜索配准参数；

步骤108：重复步骤107，直到搜索到最高层即原图大小，则得到的配准参数就是最终的配准参数；

步骤109：根据得到的配准参数将待配准图像变换回去，然后与参考图叠加得到配准结果图；

所述的步骤102，包括如下步骤：

步骤201：取一幅SAR图像I₁，寻找图像中的最小灰度值Z_min和最大灰度值Z_max，求平均得到我们的初始阈值；

$T_0=\frac{Z_{\max }+Z_{\min }}{2}$

步骤202：T_k将图像分离成两部分，一部分就是我们要的目标区域，另一部分就是我们不需要的背景区域，分别对这两个区域求其平均灰度Z₀和Z_B；

$Z_0=\frac{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

$Z_B=\frac{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

其中，像素点(i,j)的灰度值是Z(i,j)，相应灰度值出现次数为N(i,j)；

步骤203：再对得到的两个区域的灰度值取平均得到新的阈值；

$T_{k+1}=\frac{Z_0+Z_B}{2}$

步骤204：若T_k＝T_k+1结束，否则k+1→k转步骤202；

步骤205：步骤204结束后T_k即为我们想要的阈值。

所述的步骤104，包括如下步骤：

步骤301：初始化旋转、平移、尺度构成的参数解空间，设定PSO算法各项参数，随机分配粒子的位置；

步骤302：将待配准图像根据配准参数采用PV插值的方法进行变换，计算变换后的图像与参考图像之间的互信息，互信息最大的解对应的参数才是最优配准参数；

步骤303：进一步迭代，记录当前迭代下粒子的速度V_i和位置X_i，并记录当前迭代下的最优位置以及目前为止的全局最优位置。粒子的更新公式如下：

V_i＝αV_i-c₁r₁(p_i-X_i)+c₂r₂(p_g-X_i)

X_i＝X_i+βV_i

上述两式中,i＝1,…,n，n为种群大小；V_i表示粒子i的速度，X_i表示粒子i的位置，p_i表示当前这次迭代得到的粒子的最佳位置，p_g表示当前粒子群到目前为止得到的最佳位置为全局最优解；c₁,c₂为学习因子(非负)；r₁,r₂为0～1之间的随机数；α为惯性因子；

步骤304：若p_g构成的全局最优解F小于最小允许误差ε,或者迭代步数超过最大迭代次数，则结束迭代,p_g对应的解即是最终要求的解(旋转角度、水平方向平移分量、垂直方向平移分量),否则转步骤302和步骤303；

所述的步骤302，包括如下步骤：

步骤401：根据配准参数对图像做PV插值，对于待配准图像上任意一点P₀(x,y)，其最近邻四个像素点为n₁、n₂、n₃和n₄，PV插值的过程如图5所示；

步骤402：统计联合灰度直方图。对于两幅图像A和B，a和b是两幅图像中相对应的像素的灰度值，p_A(a)表示图像A边缘概率分布，p_B(b)表示图像B的边缘概率分布，p_AB(a,b)表示图像A和B的联合概率分布，联合概率分布可以通过联合灰度直方图h得到。

$p_A\left(a\right)=\underset{b}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)$

$p_B\left(b\right)=\underset{a}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)$

$p_{\mathrm{AB}}(a,b)=\frac{h(a,b)}{{\mathrm{\&Sigma;}}_{a,b}h(a,b)}$

$h=\left|\begin{array}{cccc}h\left(\mathrm{0,0}\right)& h\left(\mathrm{0,1}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& h(0,N-1)\\ h\left(\mathrm{1,0}\right)& h\left(\mathrm{1,1}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& h(1,N-1)\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ h(M-\mathrm{1,0})& h(M-\mathrm{1,1})& \&CenterDot;\&CenterDot;\&CenterDot;& h(M-1,N-1)\end{array}\right|$

步骤403：根据互信息计算公式计算两幅图像之间的互信息。

$I(A,B)=\underset{a,b}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)\log \frac{p_{\mathrm{AB}}(a,b)}{p_A\left(a\right)p_B\left(b\right)}$

所述的步骤106，包括如下步骤：

步骤501：参数初始化，给定误差δ、初始点X₀和n个线性无关的方向d₁,d₂,d₃,…,d_n；

步骤502：将待配准图像根据配准参数采用PV插值的方法进行变换，计算变换后的图像与参考图像之间的互信息，互信息最大的解对应的参数才是最优配准参数；

步骤503：从X₀出发依次沿方向d₁,d₂,d₃,…,d_n进行一维搜索，得到X₁,X₂,X₃,…,X_n.再从X_n出发沿X_n与X₀连线方向进行一维搜索，得到X_k；

步骤504：如果||X_k-X_k-1||＜δ,停止搜索，得到点X_k；否则置d(k+1,j)＝d(k,j+1)，j＝1,2,…,n，返回步骤503；

步骤505：Powell搜索结束后X即为最终的配准参数。

本发明为了降低灰度变化差异的影响，采用一种由粗到精的优化寻优策略，从而提高图像配准精度；采用图像小波分解策略，降低Powell局部寻优的运算时间复杂度。本发明与现有的技术相比具有如下有益效果：

1、本发明是在区域特征的基础上用互信息进行配准，滤出掉了一些不相关的背景区域的灰度对互信息的影响；同时对区域的灰度信息进行了一个归一化，这也削弱了灰度信息变化对互信息的影响；只使用区域特征在很大程度上也削弱了SAR图像相干斑的影响；

2、本发明从高精度配准的角度出发，采用一种PSO和Powell组合优化的由粗到精的配准策略，很好地利用了Powell局部寻优能力强的优势弥补了PSO容易陷入局部最优的缺点，极大地提高了配准的精度。同时考虑到Powell寻优慢，计算复杂度高的缺点，采用小波分解的策略，对小波分解低频部分逐层搜索，大大的缩短了算法运行的时间；

3、仿真结果表明，本发明采用的基于自适应阈值分割的PSO和Powell组合优化的SAR图像配准方法，配准精度非常高，比单一的搜索策略精度高，而且减少了算法的运行时间。

附图说明

图1是基于自适应阈值分割和组合优化的SAR图像配准主流程图；

图2是自适应阈值分割流程图；

图3是PSO算法搜索配准参数流程图；

图4是基于小波分解和Powell算法搜索配准参数流程图；

图5是PV插值示意图；

图6是加拿大渥太华(Ottawa)地区，分别于1997年5月和1997年8月由加拿大RadarSat SAR卫星拍摄的两幅已经配准好的SAR图像图6(a)和图6(b)；

图7是人工变换的待配准图和配准结果图。图7(a)就是图6(b)人工变换通过水平平移-12.11，垂直平移17.3，旋转8.77°后得到的待配准图像。图7(b)则是配准效果叠加图；

图8是参考图像图8(a)和待配准图像图8(b)。这两幅图像是通过RadarSat-2卫星在2008年和2009年获取的黄河口SAR图像中黄河的部分区域；

图9是两幅图像自适应阈值分割后的结果图；

图10是两种优化算法得到的配准结果图。图10(a)是PSO的配准结果图，图10(b)是Powell的配准结果图；

图11是最终的配准效果叠加图。

具体实施方式

本发明提出了一种基于自适应阈值分割和组合优化的SAR图像配准方法，其特征是：(1)首先对同一地区的两幅图像做自适应阈值分割，分割出目标区域；(2)通过PSO搜索算法搜索配准参数使得两幅图像之间的互信息达到最大；(3)对图像做三层小波分解，取每层低频分量，以PSO搜索得到的配准参数为Powell搜索算法的初始点，从第三层开始逐层搜索最优参数，将上一层的结果作为下一层的初始点，继续搜索直至最顶层，得到最后的配准参数；(4)通过得到的配准参数将待配准图像通过插值变换回去，并叠加得到最终的配准结果。

如图1所示。

主流程图步骤如下：

步骤101：导入两幅同一地区的SAR图像，一幅作为参考图像，另一幅作为待配准图像；

步骤102：对两幅图像分别采用自适应阈值分割提取出最佳阈值，使用最佳阈值分别将两幅图像中的目标区域图像与背景区域图像分离开来；

步骤103：将两幅图像分割出来的目标区域图像的灰度归一化到同一灰度级；

步骤104：采用粒子群算法(PSO)对两幅归一化的目标区域图像之间的互信息进行寻优，搜索到使互信息达到最大的配准参数(θ,Δx,Δy,α)，其中，θ为旋转角度，Δx为水平平移，Δy为垂直平移，α为尺度缩放；

步骤105：对两幅归一化的目标区域图像做三层小波分解，只取低频分量用于Powell算法寻优；

步骤106：将步骤104中通过PSO搜索得到的配准参数(θ,Δx,Δy,α)变为(θ,Δx/8,Δy/8,α)，作为Powell算法在小波分解第3层搜索的初始点，搜索得到配准参数(θ₁,Δx₁,Δy₁,α₁)；

步骤107：因为对两幅归一化的目标区域图像做三层小波分解后的每层之间图像大小都是两倍关系，故将得到的(θ₁,Δx₁,Δy₁,α₁)变换为(θ₁,2Δx₁,2Δy₁,α₁)作为下一层搜索的初始点，继续搜索配准参数；

步骤108：重复步骤107，直到搜索至原图大小(即第三层→第二层→第一层→原图)，得到的配准参数就是最终的配准参数；

步骤109：使用步骤108得到的最终的配准参数将待配准图像变换回去，然后与参考图像叠加得到配准结果图。

如图2所示。

所述的步骤102，包括如下步骤：

步骤201：取上述两幅同一地区的SAR图像中的一幅SAR图像I₁，寻找该图像中的最小灰度值Z_min和最大灰度值Z_max，求平均得到初始阈值；

$T_0=\frac{Z_{\max }+Z_{\min }}{2}$

步骤202：每次更新的阈值T_k将图像分离成两部分，小于等于阈值T_k的部分是目标区域图像，其它的部分是背景区域图像，分别求目标区域的平均灰度Z₀和背景区域的平均灰度Z_B；

$Z_0=\frac{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

$Z_B=\frac{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

其中，像素点(i,j)的灰度值是Z(i,j)，灰度值Z(i,j)出现次数为N(i,j)；

步骤203：再对步骤202得到的目标区域图像的灰度值和背景区域图像的灰度值取平均得到新的阈值；

$T_{k+1}=\frac{Z_0+Z_B}{2}$

步骤204：若T_k和T_k+1相等，结束，得到的T_k即为采用自适应阈值分割的最佳阈值，若Z(i,j)＜＝T_k，则像素点(i,j)为目标区域图像的像素点，若Z(i,j)＞T_k，则像素点(i,j)为背景区域图像的像素点；若T_k和T_k+1不相等，则转步骤202和步骤203，继续更新阈值T_k。

如图3所示。

所述的步骤104，包括如下步骤：

步骤301：初始化由旋转角度、水平平移、垂直平移和尺度缩放四个参数构成的参数解空间，设定PSO算法各项参数，随机分配粒子的位置；

步骤302：将待配准图像的目标区域图像根据参数解空间参数采用部分体积插值(PV插值)的方法进行变换，计算变换后的图像与参考图像的目标区域图像之间的互信息，将互信息最大的解对应的参数解空间参数作为最优配准参数；

步骤303：记录当前迭代下粒子的速度V_i和位置X_i，并记录当前迭代下的最优位置以及目前为止的全局最优位置，进一步迭代，粒子的更新公式如下：

V_i＝αV_i-c₁r₁(p_i-X_i)+c₂r₂(p_g-X_i)

X_i＝X_i+βV_i

上述两式中,i＝1,…,n，其中，n为种群大小，V_i表示粒子i的速度，X_i表示粒子i的位置，p_i表示当前这次迭代得到的粒子的最佳位置，p_g表示当前粒子群到目前为止得到的最佳位置为全局最优解，c₁,c₂为非负的学习因子，r₁,r₂为0～1之间的随机数，α为惯性因子；

步骤304：若p_g构成的全局最优解F小于最小允许误差ε,或者迭代步数超过最大迭代次数，则结束迭代,p_g对应的解即是最终要求的解，该解包括旋转角度、水平方向平移分量、垂直方向平移分量和尺度缩放,否则转步骤302和步骤303。

所述的步骤302，包括如下步骤：

步骤401：根据配准参数对待配准图像的目标区域图像做PV插值，直接更新待配准图像的目标区域图像的联合灰度直方图。对于待配准图像的目标区域图像上任意一点P₀(x,y)，其最近邻四个像素点为n₁、n₂、n₃和n₄，假设参考图像的目标区域图像为R，待配准图像的目标区域图像为F，P₀(x,y)根据配准参数计算出来的新的坐标为(x',y')，dx和dy分别表示x'和y'的小数部分，那么该点与其最近邻四个像素点的坐标差可以表示为dx、1-dx、dy和1-dy，该点与其最近邻四个像素点的联合灰度直方图可以写为下式：

h(R(n_i),F(P₀))＝h(R(n_i),F(P₀))+w_i i＝1,2,3,4

其中四个权值w_i的计算公式如下：

w₁＝(1-dx)×(1-dy)   w₂＝dx×(1-dy)

w₃＝dx×dy            w₄＝(1-dx)×dy

步骤402：对于参考图像的目标区域图像A和PV插值后的目标区域图像B，a和b是参考图像的目标区域图像A和PV插值后的目标区域图像B中相对应的像素的灰度值，p_A(a)表示参考图像的目标区域图像A的边缘概率分布，p_B(b)表示PV插值后的目标区域图像B的边缘概率分布，p_AB(a,b)表示两幅同一地区的SAR图像的目标区域图像A和B的联合概率分布，联合概率分布p_AB(a,b)可以通过联合灰度直方图h得到，边缘概率分布p_A(a)和p_B(b)也可以通过以下公式计算得出。

$h=\left|\begin{array}{cccc}h\left(\mathrm{0,0}\right)& h\left(\mathrm{0,1}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& h(0,N-1)\\ h\left(\mathrm{1,0}\right)& h\left(\mathrm{1,1}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& h(1,N-1)\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ h(M-\mathrm{1,0})& h(M-\mathrm{1,1})& \&CenterDot;\&CenterDot;\&CenterDot;& h(M-1,N-1)\end{array}\right|$

$p_{\mathrm{AB}}(a,b)=\frac{h(a,b)}{{\mathrm{\&Sigma;}}_{a,b}h(a,b)}$

$p_A\left(a\right)=\underset{b}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)$

$p_B\left(b\right)=\underset{a}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)$

步骤403：根据互信息计算公式计算两幅目标区域图像之间的互信息。

$I(A,B)=\underset{a,b}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)\log \frac{p_{\mathrm{AB}}(a,b)}{p_A\left(a\right)p_B\left(b\right)}$

如图4所示。

所述的步骤106，包括如下步骤：

步骤501：取两幅归一化的目标区域图像的小波分解第3层的低频分量用于Powell搜索，将步骤104中通过PSO搜索得到的配准参数(θ,Δx,Δy,α)变为(θ,Δx/8,Δy/8,α)，(θ,Δx/8,Δy/8,α)作为Powell算法搜索的初始点X₀，然后给定误差δ和n个线性无关的方向d₁,d₂,d₃,…,d_n；

步骤502：将待配准图像的目标区域图像根据配准参数采用PV插值的方法进行变换，计算变换后的目标区域图像与参考图像的目标区域图像之间的互信息，互信息最大的解对应的参数才是最优配准参数；

步骤503：从X₀出发依次沿方向d₁,d₂,d₃,…,d_n进行一维搜索，得到X₁,X₂,X₃,…,X_n.再从X_n出发沿X_n与X₀连线方向进行一维搜索，得到X_k；

步骤504：如果||X_k-X_k-1||＜δ,停止搜索，得到点X_k；否则置d(k+1,j)＝d(k,j+1)，j＝1,2,…,n，返回步骤503；

步骤505：Powell搜索结束时得到的X_k就是最终的配准参数(θ₁,Δx₁,Δy₁,α₁)。

本发明的效果可以通过以下仿真进一步说明：

1、仿真内容

本实验都是在CPU为Intel Pentium 2.4GHZ、内存2G、Windows7系统上使用Matlab2009进行仿真。

仿真内容1，精度测试，对变化检测图像图6(a)和图6(b)(已配准好)，自己设定参数，变换后(图7(a))再用这里提出的算法配准得到配准参数与真实设定参数进行比对，配准结果图如图7(b)所示；

仿真内容2，对参考图像图8(a)和待配准图像图8(b)利用自适应阈值分割策略，分割出图像中的目标区域。结果图如图9(a)和图9(b)所示；

仿真内容3，计算提取出的区域图像之间的互信息，利用PSO优化算法搜索配准参数使得图像间的互信息达到最大，PSO的配准结果图如图10(a)所示；

仿真内容4，对图像做三层小波分解，取每层低频分量，以PSO搜索得到的配准参数为Powell搜索算法的初始点，从第三层开始逐层搜索最优参数，使互信息达到最大，将上一层的结果作为下一层的初始点，继续搜索直至最顶层，得到最后的配准参数，Powell的配准结果图如图10(b)所示。

2、仿真实验结果及分析

仿真内容1采用的数据是加拿大渥太华(Ottawa)地区，分别于1997年5月和1997年8月由加拿大RadarSat SAR卫星拍摄。其中，图6(a)和图6(b)是变化检测图像对，图7(a)就是人工变换通过水平平移-12.11，垂直平移17.3，旋转8.77°后得到的待配准图像，配准图像对是图6(a)和图7(a)。图7(b)则是配准效果叠加图。

表1配准参数的计算

参数	水平平移	垂直平移	旋转角度	时间
					设定参数	-12.11	17.30	8.77
PSO配准参数	-12.3128	16.9513	8.7923	31s
					Powell配准参数	-12.2021	17.13	8.7831	18s

从表1可以看出我们的算法配准精度还是非常高的，达到了亚像素级。

仿真内容2、3、4采用的数据是通过RadarSat-2卫星在2008年和2009年获取的黄河口SAR图像中黄河的部分区域。图8(a)和图8(b)分别是参考图像和待配准图像。图9(a)和图9(b)是自适应阈值分割的结果图。图10(a)是PSO的配准结果图，图10(b)是Powell的配准结果图。图11是最终的配准效果叠加图。

表2 Powell+小波分解各层搜索参数

小波分层	水平平移	垂直平移	旋转角度
				PSO初始值	65.2216	-96.3827	16.2113
3层	7.9729	-11.8719	16.4728
				2层	16.2330	-22.6485	16.7712
1层	33.8012	-47.3963	17.0461
				原图	70.0991	-98.5685	17.3134

表2展示了两幅图像在做Powell精配时每层的搜索参数。以PSO得到的参数作为小波分解第3层Powell搜索的初始值，当然要对平移量都要除以8，因为每层之间图像大小都是2倍，相邻层平移参数呈2倍关系。

表3 Powell精配时间分析

图像大小	Powell	Powell+小波分解
			350×290	11.08s	9.31s
400×400	15.68s	13.44s
			700×700	65.73s	47.27s

表3展示了不同图像大小Powell算法采用小波分解和不用的运行时间对比。从表3我们可以看出随着图像大小的增大，我们采用的小波分解策略的优势逐渐体现出来，在图像较大时，能很大程度上减小算法运行时间。

综上来看，本发明采用的基于自适应阈值分割和组合优化的SAR图像配准方法，配准精度非常高，比单一的搜索策略精度高，而且减少了算法的运行时间，同时抑制了噪声和灰度变化对配准的影响。

Claims (6)

1.一种基于自适应阈值分割和组合优化的SAR图像配准方法，包括如下步骤：

步骤101：导入两幅同一地区的SAR图像，一幅作为参考图像，另一幅作为待配准图像；

步骤102：对两幅图像分别采用自适应阈值分割提取出最佳阈值，使用最佳阈值分别将两幅图像中的目标区域图像与背景区域图像分离开来；

步骤103：将两幅图像分割出来的目标区域图像的灰度归一化到同一灰度级；

步骤104：采用粒子群算法(PSO)对两幅归一化的目标区域图像之间的互信息进行寻优，搜索到使互信息达到最大的配准参数(θ,Δx,Δy,α)，其中，θ为旋转角度，Δx为水平平移，Δy为垂直平移，α为尺度缩放；

步骤105：对两幅归一化的目标区域图像做三层小波分解，只取低频分量用于Powell算法寻优；

步骤106：将步骤104中通过PSO搜索得到的配准参数(θ,Δx,Δy,α)变为(θ,Δx/8,Δy/8,α)，作为Powell算法在小波分解第3层搜索的初始点，搜索得到配准参数(θ₁,Δx₁,Δy₁,α₁)；

步骤107：因为对两幅归一化的目标区域图像做三层小波分解后的每层之间图像大小都是两倍关系，故将得到的(θ₁,Δx₁,Δy₁,α₁)变换为(θ₁,2Δx₁,2Δy₁,α₁)作为下一层搜索的初始点，继续搜索配准参数；

步骤108：重复步骤107，直到搜索至原图大小(即第三层→第二层→第一层→原图)，得到的配准参数就是最终的配准参数；

步骤109：使用步骤108得到的最终的配准参数将待配准图像变换回去，然后与参考图像叠加得到配准结果图。

2.根据权利要求1所述的基于自适应阈值分割和组合优化的SAR图像配准方法，其中，所述步骤102包括如下步骤：

步骤201：取上述两幅同一地区的SAR图像中的一幅SAR图像I₁，寻找该图像中的最小灰度值Z_min和最大灰度值Z_max，求平均得到初始阈值；

$T_0=\frac{Z_{\max }+Z_{\min }}{2}$

步骤202：每次更新的阈值T_k将图像分离成两部分，小于等于阈值T_k的部分是目标区域图像，其它的部分是背景区域图像，分别求目标区域的平均灰度Z₀和背景区域的平均灰度Z_B；

$Z_0=\frac{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

$Z_B=\frac{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

其中，像素点(i,j)的灰度值是Z(i,j)，灰度值Z(i,j)出现次数为N(i,j)；

步骤203：再对步骤202得到的目标区域图像的灰度值和背景区域图像的灰度值取平均得到新的阈值；

$T_{k+1}=\frac{Z_0+Z_B}{2}$

步骤204：若T_k和T_k+1相等，结束，得到的T_k即为采用自适应阈值分割的最佳阈值，若Z(i,j)＜＝T_k，则像素点(i,j)为目标区域图像的像素点，若Z(i,j)＞T_k，则像素点(i,j)为背景区域图像的像素点；若T_k和T_k+1不相等，则转步骤202和步骤203，继续更新阈值T_k。

3.根据权利要求1所述的基于自适应阈值分割和组合优化的SAR图像配准方法，其中，所述步骤104包括如下步骤：

步骤301：初始化由旋转角度、水平平移、垂直平移和尺度缩放四个参数构成的参数解空间，设定PSO算法各项参数，随机分配粒子的位置；

步骤302：将待配准图像的目标区域图像根据参数解空间参数采用部分体积插值(PV插值)的方法进行变换，计算变换后的图像与参考图像的目标区域图像之间的互信息，将互信息最大的解对应的参数解空间参数作为最优配准参数；

步骤303：记录当前迭代下粒子的速度V_i和位置X_i，并记录当前迭代下的最优位置以及目前为止的全局最优位置，进一步迭代，粒子的更新公式如下：

V_i＝αV_i-c₁r₁(p_i-X_i)+c₂r₂(p_g-X_i)

X_i＝X_i+βV_i

上述两式中,i＝1,…,n，其中，n为种群大小，V_i表示粒子i的速度，X_i表示粒子i的位置，p_i表示当前这次迭代得到的粒子的最佳位置，p_g表示当前粒子群到目前为止得到的最佳位置为全局最优解，c₁,c₂为非负的学习因子，r₁,r₂为0～1之间的随机数，α为惯性因子；

步骤304：若p_g构成的全局最优解F小于最小允许误差ε,或者迭代步数超过最大迭代次数，则结束迭代,p_g对应的解即是最终要求的解，该解包括旋转角度、水平方向平移分量、垂直方向平移分量和尺度缩放,否则转步骤302和步骤303。

4.根据权利要求3所述的基于自适应阈值分割和组合优化的SAR图像配准方法，其中，所述步骤302包括如下步骤：

步骤401：根据配准参数对待配准图像的目标区域图像做PV插值，直接更新待配准图像的目标区域图像的联合灰度直方图。对于待配准图像的目标区域图像上任意一点P₀(x,y)，其最近邻四个像素点为n₁、n₂、n₃和n₄，假设参考图像的目标区域图像为R，待配准图像的目标区域图像为F，P₀(x,y)根据配准参数计算出来的新的坐标为(x',y')，dx和dy分别表示x'和y'的小数部分，那么该点与其最近邻四个像素点的坐标差可以表示为dx、1-dx、dy和1-dy，该点与其最近邻四个像素点的联合灰度直方图可以写为下式：

h(R(n_i),F(P₀))＝h(R(n_i),F(P₀))+w_i i＝1,2,3,4

其中四个权值w_i的计算公式如下：

w₁＝(1-dx)×(1-dy) w₂＝dx×(1-dy)

w₃＝dx×dy w₄＝(1-dx)×dy

步骤402：对于参考图像的目标区域图像A和PV插值后的目标区域图像B，a和b是参考图像的目标区域图像A和PV插值后的目标区域图像B中相对应的像素的灰度值，p_A(a)表示参考图像的目标区域图像A的边缘概率分布，p_B(b)表示PV插值后的目标区域图像B的边缘概率分布，p_AB(a,b)表示两幅同一地区的SAR图像的目标区域图像A和B的联合概率分布，联合概率分布p_AB(a,b)可以通过联合灰度直方图h得到，边缘概率分布p_A(a)和p_B(b)也可以通过以下公式计算得出。

$h=\left|\begin{array}{cccc}h\left(\mathrm{0,0}\right)& h\left(\mathrm{0,1}\right)& ...& h(0,N-1)\\ h\left(\mathrm{1,0}\right)& h\left(\mathrm{1,1}\right)& ...& h(1,N-1)\\ ...& ...& ...& ...\\ h(M-\mathrm{1,0})& h(M-\mathrm{1,1})& ...& h(M-1,N-1)\end{array}\right|$

$p_{\mathrm{AB}}(a,b)=\frac{h(a,b)}{{\mathrm{\&Sigma;}}_{a,b}h(a,b)}$

$p_A\left(a\right)=\underset{b}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)$

$p_B\left(b\right)=\underset{a}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)$

步骤403：根据互信息计算公式计算两幅目标区域图像之间的互信息。

$I(A,B)=\underset{a,b}{\mathrm{\&Sigma;}}{p}_{\mathrm{AB}}(a,b)\log \frac{p_{\mathrm{AB}}(a,b)}{p_A\left(a\right)p_B\left(b\right)}$

5.根据权利要求1所述的基于自适应阈值分割和组合优化的SAR图像配准方法，其所述的步骤106包括如下步骤：

步骤501：取两幅归一化的目标区域图像的小波分解第3层的低频分量用于Powell搜索，将步骤104中通过PSO搜索得到的配准参数(θ,Δx,Δy,α)变为(θ,Δx/8,Δy/8,α)，(θ,Δx/8,Δy/8,α)作为Powell算法搜索的初始点X₀，然后给定误差δ和n个线性无关的方向d₁,d₂,d₃,…,d_n；

步骤502：将待配准图像的目标区域图像根据配准参数采用PV插值的方法进行变换，计算变换后的目标区域图像与参考图像的目标区域图像之间的互信息，互信息最大的解对应的参数才是最优配准参数；

步骤503：从X₀出发依次沿方向d₁,d₂,d₃,…,d_n进行一维搜索，得到X₁,X₂,X₃,…,X_n.再从X_n出发沿X_n与X₀连线方向进行一维搜索，得到X_k；

步骤504：如果||X_k-X_k-1||＜δ,停止搜索，得到点X_k；否则置d(k+1,j)＝d(k,j+1)，j＝1,2,…,n，返回步骤503；

步骤505：Powell搜索结束时得到的X_k就是最终的配准参数(θ₁,Δx₁,Δy₁,α₁)。

6.根据权利要求1所述的基于自适应阈值分割和组合优化的SAR图像配准方法，其中，步骤103中将两幅图像分割出来的目标区域图像的灰度归一化到同一灰度级包括将两幅图像分割出来的目标区域图像的灰度从0～255归一化到0～25。