CN103345757A - 多层次多特征约束下的光学和sar影像自动配准方法 - Google Patents

Abstract

一种多层次多特征约束下的光学和SAR影像自动配准方法，包括对光学影像和SAR影像预处理，然后进行多尺度水平集分割，得到面分割结果；如果有存在相似的面目标，则计算疑似同名区域的质心点坐标集，否则利用小波变换对图像进行多尺度分析，在最粗尺度影像上进行低层直线特征提取及点集匹配，得到提取低层次上的配准变换参数，并进行高层直线特征提取，同时利用低层次上的变换参数，定义控制点匹配度函数，进行高层点集匹配；最后利用KNN图从结构上精确判断出匹配的点对，同时剔除误匹配对，对匹配点对采用多项式变换模型求解变换参数，得到最终配准结果。

Description

多层次多特征约束下的光学和SAR影像自动配准方法

技术领域

本发明涉及遥感影像处理技术领域，尤其是涉及一种光学影像和SAR影像的特征配准方法。

背景技术

随着成像技术的发展，合成孔径雷达(SAR)和光学传感器系统大量涌现，为同一地区的对地观测提供了多种影像数据。SAR影像反映结构信息较好，而且具有全天时、全天候、穿透性等优点，较光学影像具有更大的侦察范围，可以发现不易被光学传感器发现的目标，但影像不含光谱信息，且目标的微波反射特性导致相同物体的SAR影像可能呈现不同的表现形式。光学影像直观地对目标进行成像，含有丰富的光谱信息，但受大气衰减、天气状况的影响较大。因此，为了获得地物的多层次特性，提高对目标描述、解读的准确性，将两种图像进行融合成为遥感图像处理的热点。图像配准是图像融合的前提，配准精度直接影响图像融合的效果，因此研究多源传感器图像的高精度配准十分必要。

同时，在应用方面，SAR影像和光学影像的配准对灾害监测、军事目标监测和卫星/无人机3-4级产品自动生产起着非常重要的作用。例如，灾害发生时，天气情况往往很恶劣，此时光学传感器因受到天气影响，无法及时进行图像拍摄。全天候、全天时的SAR系统获取的遥感影像可有规律地获取用于监测的实时数据并及时确定灾害发生的区域面积等。对比灾害发生时获取的同一地区的SAR影像和光学影像，光学影像受到云层遮挡，城区目标无法解译。

如上所述，SAR影像和光学影像的配准对遥感应用具有重要的作用，但是SAR影像斑点噪声的影响以及两种影像完全不同的光谱特性增大了影像特征提取及匹配的难度，使SAR影像与光学影像的配准成为两种影像信息有效集成的“瓶颈”。因此，SAR影像与光学影像的配准成为目前多传感器卫星图像配准的研究热点及难题，主要体现如下：

●SAR影像特有的斑点噪声模糊了图像结构，使得同名特征的精确获取十分困难。SAR成像属于主动微波遥感器，它是以斜距成像，并且成像的信息是地物目标对电磁波作用的反映。由于SAR依靠目标地物对电磁波的后向散射系数成像，因此后向散射的能量(回波强度)决定SAR影像的灰度。SAR接受到的信号是散射单元回波的矢量和，每个散射单元回波的相位和它们与传感器的距离及散射物质的特性相关，因此导致接受信号的强度并不完全由地物目标的散射系数决定，而是围绕着散射系数的值有很大的随机起伏，称之为衰落。这使得具有均匀散射系数的区域它的SAR影像中并不具有均匀的灰度，仍然呈现出很强的噪声表现，这种效应称为斑点效应。斑点噪声的存在使得SAR影像不能正确反映地物目标的散射特性，严重影响图像的质量，降低了对图像中目标的信息提取能力。

●SAR影像和光学影像地物的辐射特征存在差异：例如，建筑物密集的城市能有效地把雷达回波反射回雷达方向，也就是后向散射能量大，因而产生强的回波信号。在雷达图像上，城市呈现为具有各种特征性形状的一些亮信号的集合体。另外，SAR对强点目标成像的时候，其旁瓣电平也能在图像上显现出来，并沿着方位向和距离向成十字型分布，这在光学影像上是不会有的。

●SAR影像和光学影像地物的几何特征存在差异：二者成像过程不在同一个平面内，SAR图像成像是在距离向和方位向构成的平面内，而光学影像一般的是垂直照射地面，其一般在水平面上成像。同时，SAR影像中还有一些其它的固有特征，如透视收缩（foreshortening）、后坡延长（Elongation）、顶底位移（layover）和阴影（shadow）等几何形变现象，这些现象增加了图像配准工作的难度。

正是由于上述SAR影像的这种复杂灰度分布特性以及与光学影像目标特性的差异，使得二者的配准难度大大增加，传统的影像配准方法并不适合于光学与SAR影像的配准问题，需要对其进行改进。目前SAR影像和光学影像配准方法大致可分为两类：基于特征的配准方法和基于区域的配准方法。

（1）基于区域的配准方法

基于区域的配准方法又称为直接法或者模板匹配法，该类方法将模板上的灰度信息或者其某种变换作为基准进行配准。它将模板作为单位，并根据某种相似性度量准则，计算图像中各个位置上的当前窗口与模板之间的相似程度。由于该类算法中图像特征不需要被检测，因此相似性度量准则的设计成为算法的重点。经典的相似性度量准则主要有三类（B.Zitova等，2003）：互信息方法（mutual information methods）、相关法（correlation-like methods）、傅立叶方法（Fourier methods）。由于光学与SAR影像不同的成像机理，使得同一地物灰度可能完全不同，所以给基于区域的配准方法带来了困难。基于区域的配准方法大多都利用影像的灰度信息。但是由于SAR影像和光学不同的成像机理，导致两种影像的灰度之间存在复杂的关系，故大多基于区域的配准方法很难得到令人满意的配准结果。

（2）基于特征的配准方法

基于特征的配准算法根据提取的图像特征的不同，可分为基于点特征的图像配准算法、基于边缘（线）特征的图像配准算法、基于区域（面）特征的图像配准算法等。这类配准算法中图像同名特征的提取和匹配是两个关键的步骤。

点特征一般包括角点、高曲率点、区域的质心点、直线的交点等。常用角点提取算法有Harris角点检测算法（HARRIS等，1988），SUSAN角点检测法（Smith等，1997），Lowe提出并总结完善了SIFT算法(LOWE,1999;LOWE,2004)。SIFT算法（David G.Lowe,2004）具有旋转、缩放及仿射不变性对噪声及其视角变换不敏感、稳定性强，已经被广泛应用于目标识别和图像配准研究领域中。杨雪梅（杨雪梅等，2010）及龚俊斌等（龚俊斌等，2011）经过预处理减少光学图像与SAR图像的辐射差异，并对SIFT描述子进行了改进，对光学图像与SAR图像进行配准，但结果仍然不太稳定。

图像中大部分本质结构均可以用边缘特征代表，边缘特征能较好地剔除由姿态变化和传感器类型造成的几何畸变、成像畸变的影响。由于边缘检测计算的快捷性，目前成为基于特征的配准方法中一个较好的选择。刘佳敏等（刘佳敏等，2003）首先对SAR图像及光学图像作小波变换，然后在各自的极值点处提取边缘特征；在得到的边缘上利用均值归一化灰度相关算法进行多分辨率分级匹配，低分辨率匹配结果作为较高分辨层上搜索的依据，最终得到较好的匹配结果。Zhaojun Hu（Zhaojun Hu，2010）利用hough变换提取光学图像与SAR图像的直线段特征，并用直线段端点到直线段的距离和两条直线段的中心距离作为相似度，并使用遗传算法作为搜索策略来获取配准参数，从而获得较好的匹配效果。基于边缘匹配的方法都是以提取较好的边缘特征为条件的，然而由于SAR图像和光学图像的辐射差异，造成边缘提取上差异，从而影响匹配结果。

面特征是图像中最为稳定的特征，也是包含信息最多的特征。王子路（王子路，2007）提出利用基于自适应非线性扩散滤波和C均值分类或者联合MAP准则和规则化SRAD方法的方法作为SAR图像分割方法，而光学图像则利用四阶偏微分方程的滤波和C均值分类，分别提取闭合区域边缘特征，使用改进的Hausdorff距离作为两组特征间的相似性度量，利用遗传算法作为搜索策略估算出参数，从而得到较理想的匹配效果。张雍吉等（张雍吉等，2008）提出利用光学影像与SAR影像共有区域特征进行影像配准的方法，该方法首先提取影像的边缘特征并采用数学形态学方法得到封闭轮廓进而得到闭合区域，并把区域重心作为区域控制点，进而得到配准变换参数，从而实现SAR图像和光学图像的配准。

综上所述，基于特征的配准方法更适合异源影像的配准，然而现有的特征配准算法大多利用相似性测度来对每个目标特征是否相似进行判断，仅考虑了目标本身的特性，对于光学与SAR影像来说，由于存在上述辐射和几何差异问题，且由于SAR影像自身特征提取不精确，仅仅依靠目标本身特征进行判断是不够的。

发明内容

针对上述问题，本发明提供一种多层次多特征约束下的光学和SAR影像自动配准方法，

本发明提供的一种多层次多特征约束下的光学和SAR影像自动配准方法，包括如下步骤：

步骤1，对初始的光学影像和SAR影像分别预处理；

步骤2，对步骤1所得预处理后的光学影像和SAR影像，分别进行多尺度水平集分割，得到面分割结果；

步骤3，利用形状相似曲线相似性测度算法对步骤2获得的光学与SAR影像面分割结果进行判断，如果有存在相似的面目标，则转到步骤4，否则，执行步骤5；

步骤4，步骤3所得所有相似的面目标构成的集合形成疑似同名区域Ω₀，得到粗配准结果，对疑似同名区域Ω₀中的面分割对象计算得到疑似匹配面目标质心点集，转步骤7；

步骤5，利用小波变换对初始的光学影像和SAR影像进行多尺度分析，在最粗尺度影像上进行低层直线特征提取及低层点集匹配，得到低层次上的变换参数；

步骤6，对步骤1所得预处理后的光学和SAR影像，进行高层直线特征提取，同时利用低层次上的变换参数，定义控制点匹配度函数，进行高层点集初匹配，得到初匹配后的候选控制点对；

步骤7，基于步骤4所得疑似匹配面目标质心点集或步骤6所得初匹配后的候选控制点对，利用无向图从结构上精确判断出匹配的点对，同时剔除误匹配对，对匹配点对采用多项式变换模型求解变换参数，得到最终配准结果。

而且，步骤5包括以下子步骤，

a)利用小波变换分别对初始的光学影像和SAR影像进行多尺度分解，得到低分辨率影像；

b)低层特征提取，包括在光学影像的低分辨率影像和SAR影像的低分辨率影像上分别提取直线；

c)低层直线预处理，包括对提取的所有直线按长度从大到小排序，从最长的直线起统计斜率相近的线段，若线段间的距离值小于预设阈值则认为是属于同一条直线段，然后保留同一方向上最长的直线；

d)低层交点求取，包括对预处理后的直线求交点，将求取的交点集作为控制点集；

e)低层点集匹配，包括利用高斯混合模型进行点匹配，获取低层次上的变换参数。

而且，步骤6包括以下子步骤，

a)高层特征提取，包括在初始的光学影像和SAR影像上分别提取直线；

b)高层直线预处理，包括对提取的所有直线按长度从大到小排序，从最长的直线起统计斜率相近的线段，若线段间的距离值小于预设阈值则认为是属于同一条直线段，然后保留同一方向上最长的直线；

c)高层交点求取，包括对预处理后的直线求交点；

d)高层点集初匹配，包括利用步骤5中子步骤e)所获取的低层次上的变换参数，定义控制点匹配度函数，若控制点集中某一点对的匹配度大于预设阈值，则保留作为候选控制点对，

控制点匹配度函数定义为：

$f(P_m,P_n)=\frac{|k_1-k_1^{\′}|}{\left|k_1\right|+\left|k_1^{\′}\right|+1}+\frac{|k_2-k_2^{\′}|}{\left|k_2\right|+\left|k_2^{\′}\right|+1}+\frac{|d_1-d_2|}{d_1+d_2+1}+\frac{d_1+d_2}{2T}$

其中，k₁和k₂是光学影像上两条直线M₁、M₂的斜率，P_m为M₁与M₂的交点，k′₁和k′₂是SAR影像上两条直线N₁、N₂的斜率，P_n为N₁与N₂的交点，d₁是线段M₁的中点到N₁所在直线的距离，d₂是线段M₂的中点到N₂所在直线的距离，T为粗配准精度。

而且，步骤7包括以下子步骤，

a)设疑似匹配面目标质心点集中的质心点对或初匹配后的候选控制点对包括光学影像上的点集 $M_o=\{p_{o_1},p_{o_2},\·\·\·,p_{o_k}\}$ 和SAR影像上的点集 $M_o=\{p_{o_1},p_{o_2},\·\·\·,p_{o_k}\},$ $M_s=\{p_{s_1},p_{s_2},\·\·\·,p_{s_k}\},$ 其中

是初始匹配点对，k为初始匹配点对数量；

b)对点集M_o按照如下方式构建无向图G_o(V_o,E_o)，

点集M_o中的每一点

为图中的一个节点，构成节点集合V_o；边集E_o的构成方式为，如果

和

满足以下两个条件则它们之间存在一条边E_o(i,j)，

1）

是

的K最近邻域中的点，K为预设近邻数；

2）

其中μ为V_o中节点之间所有距离的中值，i,j取值为1，2，…,k，i≠j，对无向图G_o(V_o,E_o)构建邻接矩阵A_O，其元素如下，

$A_o(i,j)=\left\{\begin{array}{cc}1,& (i,j)\∈E_o\\ 0,& (i,j){\∉E}_o\end{array}\right.$

采用同样方式，对点集M_s构造无向图G_s(V_s,E_s),对应的邻接矩阵为A_s；

c)两个邻接矩阵的相似性定义为A＝|A_o-A_s|，其中元素记为A（i,j），两个无向图G_o(V_o,E_o)和G_s(V_s,E_s)中边差最大对应的列如下，

$l=\arg \underset{j=1,\·\·\·,k}{\max }\underset{i=1}{\overset{k}{\mathrm{\Σ}}}A(i,j)$

与所得列相应点对

被认为是误匹配对，从初始匹配点对中剔除误匹配对后返回步骤b）利用剩下的初始匹配点对重新建立无向图并计算邻接矩阵，进一步寻找误匹配点，直到 $\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}A(i,j)<{\mathrm{\&epsiv;}}_2,$ ε₂为预设阈值。

而且，步骤2包括以下子步骤，

a)对步骤1所得预处理后的光学影像和SAR影像，分别利用OTSU算法进行初始分割，获得相应分割阈值t_O*和t_S*，通过下式来分别初始化光学影像和SAR影像的零水平集函数φ_{0_O}(x,y)和φ_{0_S}(x,y)，

φ_{0_O}(x,y)＝I_{S_O}(x,y)-t_O*    (1)

φ_{0_S}(x,y)＝I_{S_S}(x,y)-t_S*    (2)

b)根据相应初始化零水平集函数，对预处理后的光学影像和SAR影像分别进行多尺度水平集分割，

对预处理后的光学影像进行多尺度水平集分割如下，

基于如下水平集能量函数获得面分割结果，

$E_{\mathrm{\&epsiv;}}^{\mathrm{CV}}(c_1,c_2,\mathrm{\&phi;})=\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)|\&dtri;\mathrm{\&phi;}(x,y)|\mathrm{dxdy}$

$+v{\&Integral;}_{\mathrm{\&Omega;}}{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$

$+{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u_0(x,y)-c_1|}^2{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$           (3)

$+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u_0(x,y)-c_2|}^2\left({1-H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\right)\mathrm{dxdy}$

其中，

(c₁,c₂,φ)表示水平集能量函数，φ为水平集函数φ(x,y)的简写，c₁,c₂分别为曲线内部和外部的灰度均值，μ表示长度项，λ₁和λ₂表示正则项，v表示面积项，Ω代表图像域，(x,y)为像素坐标，dxdy表示对x坐标方向和y坐标方向求偏导，表示水平集函数φ(x,y)的梯度的模值，u₀(x,y)表示待分割影像，H_ε(z)和δ_ε(z)分别是海氏函数H(z)和狄拉克函数δ(z)的正则化形式，z为函数自变量，ε为正参数；

采用如下的水平集演化方程，

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}(x,y)}{|\&dtri;\mathrm{\&phi;}(x,y)|}\right)-v-{\mathrm{\&lambda;}}_1(u_0(x,y)-c_1{)}^2+{\mathrm{\&lambda;}}_2(u_0(x,y)-c_2{)}^2]---\left(5a\right)$

φ(0,x,y)＝φ_{0_O}(x,y)inΩ    (5b)

其中，

表示水平集函数φ(x,y)对时间t求偏导数，

表示φ(x,y)的曲率，(5b)为初始条件，其中φ(0,x,y)表示0时刻的水平集函数，即零水平集函数φ_{0_O}(x,y)；式（5a）中φ的解为对预处理后的光学影像进行多尺度水平集分割的结果；

对预处理后的SAR影像进行多尺度水平集分割如下，

将SAR影像分割成区域Ω_c，其中c取值为{1,2}，每个区域内的概率密度函数表示为均值为u_c，有效视数为L的Gamma分布模型如下，

$P_1=\frac{L^L}{u_1(L-1)!}{\left(\frac{u_0(x,y)}{u_1}\right)}^{L-1}{e}^{-\frac{L\&CenterDot;u_0(x,y)}{u_1}}---(6-1)$

$P_2=\frac{L^L}{u_2(L-1)!}{\left(\frac{u_0(x,y)}{u_2}\right)}^{L-1}{e}^{-\frac{L\&CenterDot;u_0(x,y)}{u_2}}---(6-2)$

对于SAR影像，相应的水平集能量泛函转化为：

$F(\mathrm{\&phi;}(x,y),P_1,P_2)=\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}\&dtri;H\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}+v{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$

$-{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}(x,y)\right)\log P_1\mathrm{dxdy}$     (7)

$-{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}(1-H\left(\mathrm{\&phi;}(x,y)\right))\log P_2\mathrm{dxdy}$

对式(5a)求一阶变分，导出如下Euler-Lagrange方程:

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}(x,y)}{|\&dtri;\mathrm{\&phi;}(x,y)|}\right)-v-{\mathrm{\&lambda;}}_1\log P_1+{\mathrm{\&lambda;}}_2\log P_2]---\left(8a\right)$

φ(0,x,y)＝φ_{0_S}(x,y)inΩ    (8b)

式（8a）中φ的解为对预处理后的SAR影像进行多尺度水平集分割的结果。

因此，本发明提出的一种多层次多特征约束下的光学和SAR影像自动配准方法，同时考虑影像的线特征和面特征，可以解决光学与SAR影像自动配准问题。

附图说明

图1为本发明实施例的流程图。

具体实施方式

本发明技术方案可采用计算机软件技术实现自动运行流程。以下结合附图和实施例详细说明本发明技术方案。如图1，实施例的技术方案的流程包括以下步骤：

步骤1，对初始的光学影像和SAR影像预处理：

实施例分别对光学影像和SAR影像进行滤波处理，包括对光学影像进行高斯滤波，可采用高斯滤波器进行滤波处理；对SAR影像进行Frost滤波，可采用Frost滤波器进行滤波处理。在输入光学影像和SAR影像时，需注意光学影像和SAR影像能进行正确配准的前提是必须存在重叠区域。

步骤2，对步骤1所得预处理后的光学影像和SAR影像，分别进行多尺度水平集分割，得到面分割结果，具体包括以下步骤：

a)对步骤1所得预处理后的光学影像和SAR影像，利用OTSU算法进行初始分割，获得分割阈值，利用此阈值来初始化零水平集函数：

首先，将预处理后的光学和SAR影像通过金字塔分解为S个尺度影像，通常取S=2；接着，利用OTSU算法，即现有技术中最大类间差算法分别对光学最粗尺度影像I_{S_O}(x,y)和SAR最粗尺度影像I_{S_S}(x,y)进行计算，分别得到阈值t_O*和t_S*，则可通过下式来分别初始化光学影像和SAR影像的零水平集函数φ_{0_O}(x,y)和φ_{0_S}(x,y)，

φ_{0_O}(x,y)＝I_{S_O}(x,y)-t_O*    (1)

φ_{0_S}(x,y)＝I_{S_S}(x,y)-t_S*    (2)

b)对预处理后的光学影像进行多尺度水平集分割包括如下步骤，基于如下水平集能量函数获得面分割结果，

$E_{\mathrm{\&epsiv;}}^{\mathrm{CV}}(c_1,c_2,\mathrm{\&phi;})=\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)|\&dtri;\mathrm{\&phi;}(x,y)|\mathrm{dxdy}$

$+v{\&Integral;}_{\mathrm{\&Omega;}}{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$

$+{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u_0(x,y)-c_1|}^2{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$     (3)

$+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u_0(x,y)-c_2|}^2\left({1-H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\right)\mathrm{dxdy}$

其中，

(c₁,c₂,φ)表示水平集能量函数，φ为水平集函数φ(x,y)的简写，c₁,c₂分别为曲线内部和外部的灰度均值，μ表示长度项，λ₁和λ₂表示正则项，均为非负常数，通常取经验值，v表示面积项，通常取值为0，Ω代表图像域，(x,y)为像素坐标，dxdy表示对x坐标方向和y坐标方向求偏导，

表示水平集函数φ(x,y)的梯度的模值，u₀(x,y)表示待分割影像，H_ε(z)和δ_ε(z)分别是海氏(Heaviside)函数H(z)和狄拉克(Dirac)函数δ(z)的正则化形式，其中z为函数自变量，ε为一个正参数，可取较小值。H(z)和δ(z)函数的公式如下所示：

$H\left(z\right)=\left\{\begin{array}{cc}1,& \mathrm{ifz}\&GreaterEqual;0\\ 0& \mathrm{ifz}<0\end{array}\right.,\mathrm{\&delta;}\left(z\right)=\frac{d}{\mathrm{dz}}H\left(z\right)---(4-1)$

$H_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\&pi;}}\arctan \left(\frac{z}{\mathrm{\&epsiv;}}\right)),{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(z\right)=\frac{d}{\mathrm{dz}}{H}_{\mathrm{\&epsiv;}}\left(z\right)---(4-2)$

上述最小化问题可以通过求解能量泛函对应的Euler-Lagrange方程来实现，最后可以得到如下的水平集演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}(x,y)}{|\&dtri;\mathrm{\&phi;}(x,y)|}\right)-v-{\mathrm{\&lambda;}}_1(u_0(x,y)-c_1{)}^2+{\mathrm{\&lambda;}}_2(u_0(x,y)-c_2{)}^2]---\left(5a\right)$

φ(0,x,y)＝φ_{0_O}(x,y)inΩ    (5b)

其中，

表示水平集函数φ(x,y)对时间t求偏导数，

表示φ(x,y)的曲率，(5b)为初始条件，其中φ(0,x,y)表示0时刻的水平集函数，即零水平集函数φ_{0_O}(x,y)。

对预处理后的SAR影像进行多尺度水平集分割包括如下步骤，考虑到相干斑噪声影响，将SAR影像分割成区域Ω_c，其中c取值为{1,2}，以区分两类，每个区域内的概率密度函数可以表示为均值为u_c，有效视数为L的Gamma分布模型：

$P_1=\frac{L^L}{u_1(L-1)!}{\left(\frac{u_0(x,y)}{u_1}\right)}^{L-1}{e}^{-\frac{L\&CenterDot;u_0(x,y)}{u_1}}---(6-1)$

$P_2=\frac{L^L}{u_2(L-1)!}{\left(\frac{u_0(x,y)}{u_2}\right)}^{L-1}{e}^{-\frac{L\&CenterDot;u_0(x,y)}{u_2}}---(6-2)$

对于SAR影像，相应的水平集能量泛函可以转化为：

$F(\mathrm{\&phi;}(x,y),P_1,P_2)=\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}\&dtri;H\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}+v{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$

$-{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}(x,y)\right)\log P_1\mathrm{dxdy}$     (7)

$-{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}(1-H\left(\mathrm{\&phi;}(x,y)\right))\log P_2\mathrm{dxdy}$

对能量函数(5a)求一阶变分，可以导出如下Euler-Lagrange方程:

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}(x,y)}{|\&dtri;\mathrm{\&phi;}(x,y)|}\right)-v-{\mathrm{\&lambda;}}_1\log P_1+{\mathrm{\&lambda;}}_2\log P_2]---\left(8a\right)$

φ(0,x,y)=φ_{0_S}(x,y)inΩ    (8b)

式（5a）中φ的解即为光学图像分割结果，式（8a）中φ的解即为SAR影像分割结果。

步骤3，利用形状曲线相似性测度算法对获得的光学与SAR影像面分割结果进行判断，如果有存在相似的面目标，则转到步骤4，否则，执行步骤5；形状曲线相似性测度算法为现有技术，本发明不予赘述；

步骤4，步骤3所得所有相似的面目标构成的集合图中的疑似匹配面目标集合形成疑似同名区域Ω₀，得到粗配准结果，对疑似同名区域Ω₀中的面分割对象计算其质心点坐标集，计算方式为现有技术，图中记为疑似匹配面目标质心点集，光学与SAR影像中对应面分割对象的质心点构成初始匹配点对，转步骤7；

步骤5，利用小波变换对初始的光学和SAR影像进行多尺度分析，在最粗尺度影像上进行低层直线特征提取及点集匹配，得到低层次上的配准变换参数，具体包括以下步骤：

a)利用小波变换分别对光学和SAR影像进行多尺度分解，分别获得光学最粗尺度影像Y_opt(0)和SAR最粗尺度影像Y_sar(0)，使得多尺度分解后的最粗尺度影像大小小于256×256像素，从而获得低分辨率影像；

b)低层特征提取:在a）所得低分辨率影像上提取线特征，其中光学影像的低分辨率影像上采用LSD算法提取直线，SAR影像的低分辨率影像上先采用D1算法提取边缘特征再采用Hough变换提取直线特征，具体提取方式为现有技术；

c)低层直线预处理：对b）所得所有直线按长度从大到小排序，从最长的直线起统计斜率相近的线段(可通过斜率差值是否满足预设斜率差值阈值来确定相近的线段)，计算线段间的距离，若距离值小于预设距离阈值（一般设为3），则认为是属于同一条直线段，删除长度较小的直线，只保留同一方向上最长的直线；

d)低层交点求取：对c）预处理后的直线求交点，为了防止斜率相近的线相交得到极点，仅对线集中夹角为30°到150°的线段求取交点，将求取的交点集作为控制点集；

e)低层点集匹配：基于d）所得控制点集，利用高斯混合模型（Gaussian Mixture Model,GMM）进行点匹配，获取低层次上的变换参数，高斯混合模型为现有技术，为便于参考起见，提供具体实现说明如下。

从d）所得控制点集中，取光学影像的低分辨率影像上的交点集作为基准点集，取SAR影像的低分辨率影像上的交点集作为待匹配点集。假设基准点集为Y_M＝[y₁,y₂,...,y_M]^T，其中任一基准点记为y_m，m的取值为1,2,…,M；待匹配点集为X_N＝[x₁,x₂,...,x_N]^T，其中任一待配点记为x_n，n的取值为1,2,…,N；将基准点集中的各点视为高斯混合模型中各个高斯分量的质心，另外考虑到待配点集中存在噪声等相对于基准点中集无匹配点对（缺失点），因而在基准点集中增加一个虚拟基准点及相应高斯分量，用P(y_M+1)表示缺失点的比例。定义该GMM的概率密度函数为：

$p\left(x_n\right)=\underset{m=1}{\overset{M+1}{\mathrm{\&Sigma;}}}P\left(y_m\right)p\left(x_n\right|y_m)$

其中，m的取值为1,2,…,M+1，n的取值为1,2,…,N，p(x_n)表示第n个点对应的高斯混合模型的概率密度函数，为第m个高斯分量的概率密度，x_n表示第n个待匹配点y_m表示第m个基准点，σ表示每个高斯分量的协方差。

且设为待匹配点集相对于基准点集中缺失点的均匀分布。各个分量的混合概率P(y_m)均相同： $P\left(y_m\right)=\frac{1}{M}.$

任意两点y_m和x_n的匹配概率由后验概率定义：

P(y_m|x_n)＝P(y_m)p(x_n|y_m)/p(x_n)

其中，m的取值为1,2,…,M+1，n的取值为1,2,…,N，假设GMM的各个高斯分量质心按变换参数集

进行几何变换，同时其协方差σ²也随之变化。利用EM算法通过对数似然函数来估计变换参数集

和匹配关系P(y_m|x_n)：

$Q=-\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{m=1}{\overset{M+1}{\mathrm{\&Sigma;}}}{P}^{\mathrm{old}}\left(y_m\right|x_n)\log \left(P^{\mathrm{new}}\left(y_m\right)p^{\mathrm{new}}\left(x_n\right|y_m)\right)$

其中，m的取值为1,2,…,M+1，n的取值为1,2,…,N，Q为目标观测方程，P^old(y_m|x_n)表示根据当前变换参数利用Bayesian原理计算出的后验概率，p^new(y_m)表示更新后的混合概率，p^new(x_n|y_m)表示通过最小化目标函数得到的更新高斯分量概率密度。

变换参数集

包括仿射变换参数B，t，低分辨率上变换模型采用仿射变换模型，定义为T(y_m;B,t)＝By_m+t，从而有：

$Q(B,t,{\mathrm{\&sigma;}}^2)=\frac{1}{{2\mathrm{\&sigma;}}^2}\underset{m,n=1}{\overset{M,N}{\mathrm{\&Sigma;}}}{P}^{\mathrm{old}}\left(y_m\right|x_n){\left|\right|x_n-({\mathrm{By}}_m+t)\left|\right|}^2+\log {\mathrm{\&sigma;}}^2\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{p}^{\mathrm{old}}\left(y_m\right|x_n)$

其中，m的取值为1,2,…,M，n的取值为1,2,…,N。

采用EM算法求解目标方程，通过E-步和M-步交替进行计算，E-步利用更新后的点位置矢量重新计算后验概率，M-步估计点集之间的匹配关系和点集之间的几何变换矩阵，直至收敛。收敛条件为第N次与第N+1次中经变换后的点集与基准点集中有匹配关系的点的平均误差距离的差值Dist_err小于0.01或迭代次数Iter_max大于100。具体求解过程：

（1）初始化参数：B=I，t=0,0≤w≤1，

其中，m的取值为1,2,…,M+1，n的取值为1,2,…,N；

其中，I为单位矩阵，即2×2矩阵[10;01]。

（2）E-步计算匹配概率p_mn：

$p_{\mathrm{mn}}=\frac{\exp (-\frac{1}{{2\mathrm{\&sigma;}}^2}{\left|\right|x_n-({\mathrm{By}}_m+t\left|\right|}^2}{\underset{r=1}{\overset{M}{\mathrm{\&Sigma;}}}\exp (-\frac{1}{{2\mathrm{\&sigma;}}^2}{\left|\right|x_n-({\mathrm{By}}_r+t\left|\right|}^2+\left({2\mathrm{\&pi;\&sigma;}}^2\right)\frac{w}{1-w}\frac{M}{N}}$

其中，变量r的取值为1,2,…,M。

（3）M-步计算变换参数B,t,σ²:

N_P＝1^TP1

${\mathrm{\&mu;}}_x=E\left(X\right)=\frac{1}{N_P}{X}^T{P}^{T}1$

${\mathrm{\&mu;}}_y=E\left(Y\right)=\frac{1}{N_P}{Y}^T{P}^{T}1$

$X^{\&prime;}=X-1{\mathrm{\&mu;}}_x^T$

$Y^{\&prime;}=Y-1{\mathrm{\&mu;}}_y^T$

B＝(X^'TP^TY')(Y^'Td(P1)Y')^-1

t＝μ_x-Bμ_y

${\mathrm{\&sigma;}}^2=\frac{1}{2N_P}(\mathrm{tr}\left(X^{\&prime;T}d\left(P^{T}1\right)X^{\&prime;}\right)-\mathrm{tr}\left(X^{\&prime;T}{P}^T{Y}^{\&prime;}{B}^T\right))$

其中，1可视为元素全为数字1的列矢量，N_P为匹配概率之和，μ_x、μ_y分别为在匹配关系P下待配点集和基准点集取均值E(·)，X'、Y'分别为待配点集与基准点集消除平移影响后的新坐标值；d(P1)表示以取对角矩阵，tr(·)表示矩阵的迹。

从而变换后的点集为：Γ(Y)＝YB^T+1t^T。

步骤6，对步骤1所得预处理后的光学和SAR影像，进行高层直线特征提取，同时利用低层次上的变换参数，定义控制点匹配度函数，进行高层点集初匹配，包括以下步骤：

a)高层特征提取:在原始分辨率影像（即初始的光学和SAR影像）上提取线特征，其中光学影像采用LSD算法，SAR影像先采用D1算法提取边缘特征再采用Hough变换提取直线特征；

b)高层直线预处理：对所有直线按长度从大到小排序，从最长的直线起统计斜率相近的线段，计算线段间的距离，若距离值小于阈值3，则认为是属于同一条直线段，删除长度较小的直线，只保留同一方向上最长的直线；

c)高层交点求取：对预处理后的直线求交点，为了防止斜率相近的线相交得到极点，可以仅对线集中夹角为30°到150°的线段求取交点，将求取的交点集作为控制点集；

d)高层点集初匹配：利用步骤5中子步骤e)所获取的低层次上的变换参数，定义控制点匹配度函数，若控制点集中某一点对的匹配度大于阈值，则认为是匹配的点对：

控制点匹配度函数定义为：

$f(P_m,P_n)=\frac{|k_1-k_1^{\&prime;}|}{\left|k_1\right|+\left|k_1^{\&prime;}\right|+1}+\frac{|k_2+k_2^{\&prime;}|}{\left|k_2\right|+\left|k_2^{\&prime;}\right|+1}+\frac{|d_1-d_2|}{d_1+d_2+1}+\frac{d_1+d_2}{2T}$

其中，k₁和k₂是光学影像上两条直线M₁、M₂的斜率，P_m为M₁与M₂的交点，k'₁和k'₂是SAR影像上两条直线N₁、N₂的斜率，P_n为N₁与N₂的交点，d₁是线段M₁的中点到N₁所在直线的距离，d₂是线段M₂的中点到N₂所在直线的距离。匹配度函数的阈值可由三个归一化差值系数、d₁和d₂的最大值及粗配准精度T决定。例如，可预设归一化斜率差值系数（即上式的第（1）项和第（2）项）的阈值为0.1，归一化距离差值系数（即上式的第（3）项）的阈值为0.1,d₁和d₂的最大值均为5个像素，T取25个像素。

初匹配中若某一点对的匹配度函数值大于阈值0.5，则认为它们不是正确的匹配点对并将其剔除，否则保留点对作为候选控制点对。

步骤7，基于步骤4所得疑似匹配面目标质心点集中的质心点对或步骤6所得初匹配后的候选控制点对，利用KNN图从结构上精确判断出匹配的点对，同时剔除误匹配对，对匹配点对采用多项式变换模型求解变换参数，得到最终配准结果。具体包括以下步骤：

a)设疑似匹配面目标质心点集中的质心点对或初匹配后的候选控制点对包括：光学影像上的点集 $M_o=\{p_{o_1},p_{o_2},\&CenterDot;\&CenterDot;\&CenterDot;,p_{o_k}\},$ SAR影像上的点集 $M_s=\{p_{s_1},p_{s_2},\&CenterDot;\&CenterDot;\&CenterDot;,p_{s_k}\},$ 其中

是初始匹配点对；

b)对点集M_o按照如下方式构建KNN图G_o(V_o,E_o)：点集中的每一点

为图中的一个节点v_oi,即节点集合 $V_o=\{v_{o1},v_{o2},\&CenterDot;\&CenterDot;\&CenterDot;,v_{\mathrm{ok}}\}=\{p_{o_1},p_{o_2},\&CenterDot;\&CenterDot;\&CenterDot;,p_{o_k}\}.$ 边集E_o的构成方式为，如果

和

满足以下两个条件则它们之间存在一条边E_o(i,j)：1）是

的K最近邻域中的点；2） $\left|\right|p_{o_j}-p_{o_i}\left|\right|\&le;\mathrm{\&mu;},$ 其中μ为V_o中节点之间所有距离的中值，即 $\mathrm{\&mu;}={\mathrm{med}}_{(h,l)\&Element;V_O\&times;V_O}\left|\right|p_{o_h}-p_{o_l}\left|\right|.$ 预设近邻数K可取经验值，一般为8。i,j分别为光学和SAR初始匹配点对索引，取值为1，2，…,k，i≠j。

对无向图G_o(V_o,E_o)构建邻接矩阵A_o，其元素如下：

$A_o(i,j)=\left\{\begin{array}{cc}1,& (i,j)\&Element;E_o\\ 0,& (i,j){\&NotElement;E}_o\end{array}\right.$

同理，对点集M_s构造无向图G_s(V_s,E_s),对应的邻接矩阵为A_s：

点集中的每一点M_s为图中的一个节点,即节点集合

边集E_s的构成方式为，如果和

满足以下两个条件则它们之间存在一条边E_s(i,j)：1）是

的K最近邻域中的点；2）

其中μ为V_s中节点之间所有距离的中值。

对无向图G_s(V_s,E_s)构建邻接矩阵A_s，其元素为

$A_s(i,j)=\left\{\begin{array}{c}1,(i,j)\&Element;E_s\\ 0,(i,j){\&NotElement;E}_s\end{array}\right.$

c)两个邻接矩阵的相似性定义为：A＝|A_o-A_s|，其中元素记为A（i,j）。两个无向图中边差最大对应的列为：

$l=\arg \underset{j=1,\&CenterDot;\&CenterDot;\&CenterDot;,k}{\max }\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}A(i,j)$

与上述公式所得列对应的点对

被认为是误匹配对，在从初始匹配点对中剔除误匹配对后返回步骤b）利用剩下的初始匹配点对重新建立KNN图并计算KNN图的邻接矩阵，进一步寻找误匹配点，直到

理论上ε₂应为0，实际上由于提取的线特征可能存在误差等原因，阈值ε₂可取为

，k为步骤2所得疑似匹配面目标质心点集中的质心点对或步骤6后得到的初匹配后的候选控制点对数量，

为向下取整。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (5)

1.一种多层次多特征约束下的光学和SAR影像自动配准方法，其特征在于，包括以下步骤：

步骤1，对初始的光学影像和SAR影像分别预处理；

步骤2，对步骤1所得预处理后的光学影像和SAR影像，分别进行多尺度水平集分割，得到面分割结果；

步骤3，利用形状相似曲线相似性测度算法对步骤2获得的光学与SAR影像面分割结果进行判断，如果有存在相似的面目标，则转到步骤4，否则，执行步骤5；

步骤4，步骤3所得所有相似的面目标构成的集合形成疑似同名区域Ω₀，得到粗配准结果，对疑似同名区域Ω₀中的面分割对象计算得到疑似匹配面目标质心点集，转步骤7；

步骤5，利用小波变换对初始的光学影像和SAR影像进行多尺度分析，在最粗尺度影像上进行低层直线特征提取及低层点集匹配，得到低层次上的变换参数；

步骤6，对步骤1所得预处理后的光学和SAR影像，进行高层直线特征提取，同时利用低层次上的变换参数，定义控制点匹配度函数，进行高层点集初匹配，得到初匹配后的候选控制点对；

步骤7，基于步骤4所得疑似匹配面目标质心点集或步骤6所得初匹配后的候选控制点对，利用无向图从结构上精确判断出匹配的点对，同时剔除误匹配对，对匹配点对采用多项式变换模型求解变换参数，得到最终配准结果。

2.根据权利要求1所述多层次多特征约束下的光学和SAR影像自动配准方法，其特征在于：步骤5包括以下子步骤，

a)利用小波变换分别对初始的光学影像和SAR影像进行多尺度分解，得到低分辨率影像；

b)低层特征提取，包括在光学影像的低分辨率影像和SAR影像的低分辨率影像上分别提取直线；

c)低层直线预处理，包括对提取的所有直线按长度从大到小排序，从最长的直线起统计斜率相近的线段，若线段间的距离值小于预设阈值则认为是属于同一条直线段，然后保留同一方向上最长的直线；

d)低层交点求取，包括对预处理后的直线求交点，将求取的交点集作为控制点集；

e)低层点集匹配，包括利用高斯混合模型进行点匹配，获取低层次上的变换参数。

3.根据权利要求1所述多层次多特征约束下的光学和SAR影像自动配准方法，其特征在于：步骤6包括以下子步骤，

a)高层特征提取，包括在初始的光学影像和SAR影像上分别提取直线；

b)高层直线预处理，包括对提取的所有直线按长度从大到小排序，从最长的直线起统计斜率相近的线段，若线段间的距离值小于预设阈值则认为是属于同一条直线段，然后保留同一方向上最长的直线；

c)高层交点求取，包括对预处理后的直线求交点；

d)高层点集初匹配，包括利用步骤5中子步骤e)所获取的低层次上的变换参数，定义控制点匹配度函数，若控制点集中某一点对的匹配度大于预设阈值，则保留作为候选控制点对，

控制点匹配度函数定义为：

$f(P_m,P_n)=\frac{|k_1-k_1^{\&prime;}|}{\left|k_1\right|+\left|k_1^{\&prime;}\right|+1}+\frac{|k_2-k_2^{\&prime;}|}{\left|k_2\right|+\left|k_2^{\&prime;}\right|+1}+\frac{|d_1-d_2|}{d_1+d_2+1}+\frac{d_1+d_2}{2T}$

其中，k₁和k₂是光学影像上两条直线M₁、M₂的斜率，P_m为M₁与M₂的交点，k'₁和k'₂是SAR影像上两条直线N₁、N₂的斜率，P_n为N₁与N₂的交点，d₁是线段M₁的中点到N₁所在直线的距离，d₂是线段M₂的中点到N₂所在直线的距离，T为粗配准精度。

4.根据权利要求1所述多层次多特征约束下的光学和SAR影像自动配准方法，其特征在于：步骤7包括以下子步骤，

a)设疑似匹配面目标质心点集中的质心点对或初匹配后的候选控制点对包括光学影像上的点集 $M_o=\{p_{o_1},p_{o_2},\&CenterDot;\&CenterDot;\&CenterDot;,p_{o_k}\}$ 和SAR影像上的点集 $M_o=\{p_{o_1},p_{o_2},\&CenterDot;\&CenterDot;\&CenterDot;,p_{o_k}\},$ $M_s=\{p_{s_1},p_{s_2},\&CenterDot;\&CenterDot;\&CenterDot;,p_{s_k}\},$ 其中

是初始匹配点对，k为初始匹配点对数量；

b)对点集M_o按照如下方式构建无向图G_o(V_o,E_o)，

点集M_o中的每一点

为图中的一个节点，构成节点集合V_o；边集E_o的构成方式为，如果

和

满足以下两个条件则它们之间存在一条边E_o(i,j)，

1）

是的K最近邻域中的点，K为预设近邻数；

2）

其中μ为V_o中节点之间所有距离的中值，i,j取值为1，2，…,k，i≠j，对无向图G_o(V_o,E_o)构建邻接矩阵A_O，其元素如下，

$A_o(i,j)=\left\{\begin{array}{cc}1,& (i,j)\&Element;E_o\\ 0,& (i,j){\&NotElement;E}_o\end{array}\right.$

采用同样方式，对点集M_s构造无向图G_s(V_s,E_s),对应的邻接矩阵为A_s；

c)两个邻接矩阵的相似性定义为A＝|A_o-A_s|，其中元素记为A（i,j），两个无向图G_o(V_o,E_o)和G_s(V_s,E_s)中边差最大对应的列如下，

$l=\arg \underset{j=1,\&CenterDot;\&CenterDot;\&CenterDot;,k}{\max }\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}A(i,j)$

与所得列相应点对

被认为是误匹配对，从初始匹配点对中剔除误匹配对后返回步骤b）利用剩下的初始匹配点对重新建立无向图并计算邻接矩阵，进一步寻找误匹配点，直到 $\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}A(i,j)<{\mathrm{\&epsiv;}}_2,$ ε₂为预设阈值。

5.根据权利要求1或2或3或4所述多层次多特征约束下的光学和SAR影像自动配准方法，其特征在于：步骤2包括以下子步骤，

a)对步骤1所得预处理后的光学影像和SAR影像，分别利用OTSU算法进行初始分割，获得相应分割阈值t_O*和t_S*，通过下式来分别初始化光学影像和SAR影像的零水平集函数φ_{0_O}(x,y)和φ_{0_S}(x,y)，

φ_{0_O}(x,y)＝I_{S_O}(x,y)-t_O*    (1)

φ_{0_S}(x,y)＝I_{S_S}(x,y)-t_S*    (2)

b)根据相应初始化零水平集函数，对预处理后的光学影像和SAR影像分别进行多尺度水平集分割，

对预处理后的光学影像进行多尺度水平集分割如下，

基于如下水平集能量函数获得面分割结果，

$E_{\mathrm{\&epsiv;}}^{\mathrm{CV}}(c_1,c_2,\mathrm{\&phi;})=\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)|\&dtri;\mathrm{\&phi;}(x,y)|\mathrm{dxdy}$

$+v{\&Integral;}_{\mathrm{\&Omega;}}{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$            (3)

$+{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u_0(x,y)-c_1|}^2{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$

$+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u_0(x,y)-c_2|}^2\left({1-H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)\right)\mathrm{dxdy}$

其中，

(c₁,c₂,φ)表示水平集能量函数，φ为水平集函数φ(x,y)的简写，c₁,c₂分别为曲线内部和外部的灰度均值，μ表示长度项，λ₁和λ₂表示正则项，v表示面积项，Ω代表图像域，(x,y)为像素坐标，dxdy表示对x坐标方向和y坐标方向求偏导，

表示水平集函数φ(x,y)的梯度的模值，u₀(x,y)表示待分割影像，H_ε(z)和δ_ε(z)分别是海氏函数H(z)和狄拉克函数δ(z)的正则化形式，z为函数自变量，ε为正参数；

采用如下的水平集演化方程，

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}(x,y)}{|\&dtri;\mathrm{\&phi;}(x,y)|}\right)-v-{\mathrm{\&lambda;}}_1(u_0(x,y)-c_1{)}^2+{\mathrm{\&lambda;}}_2(u_0(x,y)-c_2{)}^2]---\left(5a\right)$

φ(0,x,y)＝φ_{0_O}(x,y)inΩ     (5b)

其中，

表示水平集函数φ(x,y)对时间t求偏导数，表示φ(x,y)的曲率，(5b)为初始条件，其中φ(0,x,y)表示0时刻的水平集函数，即零水平集函数φ_{0_O}(x,y)；式（5a）中φ的解为对预处理后的光学影像进行多尺度水平集分割的结果；

对预处理后的SAR影像进行多尺度水平集分割如下，

将SAR影像分割成区域Ω_c，其中c取值为{1,2}，每个区域内的概率密度函数表示为均值为u_c，有效视数为L的Gamma分布模型如下，

$P_1=\frac{L^L}{u_1(L-1)!}{\left(\frac{u_0(x,y)}{u_1}\right)}^{L-1}{e}^{-\frac{L\&CenterDot;u_0(x,y)}{u_1}}---(6-1)$

$P_2=\frac{L^L}{u_2(L-1)!}{\left(\frac{u_0(x,y)}{u_2}\right)}^{L-1}{e}^{-\frac{L\&CenterDot;u_0(x,y)}{u_2}}---(6-2)$

对于SAR影像，相应的水平集能量泛函转化为：

$F(\mathrm{\&phi;}(x,y),P_1,P_2)=\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}\&dtri;H\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}+v{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}(x,y)\right)\mathrm{dxdy}$

$-{\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}(x,y)\right)\log P_1\mathrm{dxdy}$     (7)

$-{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}(1-H\left(\mathrm{\&phi;}(x,y)\right))\log P_2\mathrm{dxdy}$

对式(5a)求一阶变分，导出如下Euler-Lagrange方程:

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}(x,y)\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}(x,y)}{|\&dtri;\mathrm{\&phi;}(x,y)|}\right)-v-{\mathrm{\&lambda;}}_1\log P_1+{\mathrm{\&lambda;}}_2\log P_2]---\left(8a\right)$

φ(0,x,y)＝φ_{0_S}(x,y)inΩ    (8b)

式（8a）中φ的解为对预处理后的SAR影像进行多尺度水平集分割的结果。

Images