边缘和区域概率密度差相结合的水平集SAR图像分割方法

摘要

本发明公开了一种基于边缘和区域概率密度差相结合的水平集SAR图像分割方法，它属于图像处理技术领域，主要解决现有水平集方法难以分割边缘模糊的SAR图像和对SAR图像真实边缘定位不准的问题。其实现步骤包括：首先应用ROEWA算子检测SAR图像边缘强度模值|R

说明书

技术领域

本发明属于图像处理技术领域，涉及SAR图像分割领域的应用，具体地说是一种基于边缘和区域相结合的水平集SAR图像分割方法，可用于SAR图像的分割、边缘检测和目标识别。

背景技术

合成孔径雷达SAR是一种高分辨率主动式雷达，具有全天候、全天时、分辨率高、可侧视成像等优点，可应用于军事、农业、导航、地理监视等诸多领域。SAR图像被广泛的应用在目标检测领域，而SAR图像分割则是从图像处理到图像分析的重要步骤，是目标分类和识别的基础。由于SAR是一种相干成像系统，SAR图像在本质上是对目标的电磁散射特性和结构特性的反应，其成像效果很大程度上依赖于雷达参数和地域电磁参数。由于SAR成像系统的特殊性，使SAR图像的信息表达方式和光学图像有很大差异，并且会受到相干斑噪声及阴影等许多几何特征的影响，正因为这些因素使适用于光学图像的分割方法对SAR图像不再适用。

总体上说常规的SAR分割方法差不多还是用传统的分割模型对SAR图像进行分割，如一些基于阈值的分割、基于边缘的分割和基于区域生长的分割方法等，但是其中的很多分割方法容易受到SAR图像中相干斑噪声的影响，所以此类方法都需要先通过预处理抑制相干斑噪声。这些预处理方法在抑制相干斑噪声的同时，不可避免地损失了很多边缘与目标信息，影响了分割效果。

水平集图像分割方法是重要的图像处理方法之一。水平集方法的优点是可以适应拓扑结构的变化，而且算法稳定性较高。利用水平集方法研究SAR图像分割问题，通过建立合适的能量泛函，可以在能量泛函的定义中引入图像区域信息，在不需要相干斑预处理的情况下，对于受相干斑噪声影响的SAR图像可以获得比较准确的分割结果。

传统的水平集模型分基于边缘和基于区域两种方法。测地线活动轮廓模型GAC模型能在不附加任何外界控制条件下，自由处理曲线运动时的拓扑结构变化，是基于边缘的模型中应用最广的一种。但是该模型仅利用图像的梯度信息，由于图像运动中的边缘并非都是理想的阶梯边缘，难以分割出边缘模糊的图像中的同质区域，另外当图像中的对象有较深的凹陷边界时，GAC模型可能使演化曲线停止在某一局部极小值状态，并不与对象的边界相一致，且梯度信息本身就不适用于SAR图像。Chan和Vese提出的CV模型作为一种能有效提高曲线演化拓扑自适应能力的区域水平集分割模型，利用了图像匀质区域的全局信息。但是其同样也存在不足之处。该模型的内部能量项只是保证零水平集曲线光滑，而没有考虑水平集函数本身所固有的内在性质，在某些应用中还需要对水平集函数进行重新初始化，以使它接近符号距离函数，保证数值解法的稳定性，另外，模型中仅仅利用了图像的区域信息，而另一个重要的图像的边缘信息没有很好地利用，在一些分割应用中尤其是对灰度不均匀的图像可能产生图像边缘定位不准确的缺陷。

虽然水平集分割方法已经在光学和医学图像分割中获得了巨大的成功，但是在SAR图像分割领域，研究的还比较少，国际上知名的有法国的Refregier等人和加拿大的Ayed等人所领导的两个机构，并且两个机构各有偏重，前者偏重研究参数活动轮廓模型在SAR图像分割和边缘检测方面的应用，后者主要研究几何活动轮廓模型对多区域SAR图像的分割。而国内在水平集SAR图像分割方面研究还很少。

发明内容

本发明的目的在于克服上述已有技术和方法中的不足，提供一种基于边缘和区域相结合的水平集SAR图像分割方法，以实现对边缘模糊的SAR图像区域的准确分割，提高边缘的定位性能，从而提高SAR图像的分割的质量。

实现本发明目的技术思路是：根据SAR图像的统计特性，应用Gamma分布作为SAR图像的概率密度函数，并且分别计算由水平集划分得到的两个区域的概率密度的真实值和估计值，以这两个值差的平方作为区域信息，同时应用对SAR图像有很好边缘定位性能的指数加权均值比率(ROEWA)边缘检测算子的检测结果代替梯度值作为边缘信息，将这两部分信息结合进行能量函数的建模，从而达到理想的分割结果。其具体的实现过程如下：

(1)对待分割SAR图像I应用ROEWA边缘检测算子进行边缘检测，得到边缘强度模值|R_max|；

(2)将水平集函数φ初始化成圆形的符号距离函数形式，根据水平集函数值的正负，将SAR图像分割成两个区域Ω₁和Ω₂；

(3)根据所述的两个区域Ω₁和Ω₂，利用下式分别计算对应的估计概率密度和以及真实估计概率密度p₁和p₂：

3a)利用下式计算区域Ω₁和Ω₂的强度均值：

$c_i=\frac{{\int \int }_{{\Omega }_i}I(x,y)\mathrm{dxdy}}{{\int \int }_{{\Omega }_i}\mathrm{dxdy}},$其中i＝1，2，(x，y)是图像坐标；

3b)结合步骤3a)得到的强度均值，利用下式计算区域Ω₁和Ω₂的估计概率密度和

$p_i^E\left(I\right)=\frac{L^L}{c_i\Gamma \left(L\right)}{\left(\frac{I}{c_i}\right)}^{L-1}{e}^{-\mathrm{LI}/c_i},$

其中L是SAR图像的视数，Γ(·)是Gamma函数，i＝1，2；

3c)利用下式计算区域Ω₁和Ω₂的实际概率密度p₁和p₂：

$p_i\left(I\right)=\frac{{\int \int }_{{\Omega }_i}\delta (I(x_1,y_1)-I){\mathrm{dx}}_1{\mathrm{dy}}_1}{{\int \int }_{{\Omega }_i}{\mathrm{dx}}_1{\mathrm{dy}}_1},$

其中I(x₁，y₁)代表坐标为(x₁，y₁)的图像像素值，i＝1，2；

(4)结合步骤(1)、步骤(2)、步骤(3)构造总的分割能量函数模型E_SAR：

$E_{\mathrm{SAR}}=\underset{i=1}{\overset{2}{\Sigma }}{\int \int }_{{\Omega }_i}{\lambda }_i{(p_i^E-p_i)}^2+\mu {\int \int }_{\Omega }g\left(\right|R_{\max }\left|\right)|▿H\left(\phi \right)|+v{\int \int }_{\Omega }\frac{1}{2}{\left(\right|▿\phi |-1)}^2,$

其中Ω是整个图像区域，即Ω＝Ω₁+Ω₂，λ_i，i＝1，2是第一部分即区域能量项的权值，μ是第二部分边缘能量项的权重，v是第三部分修正项的权值，λ_i＞0，μ＞0，v＞0，是梯度算子，H(·)是Heaviside函数，g(|R_max|)是定义在边缘强度模值|R_max|上的指示函数，表达式如下：

$g\left(\right|R_{\max }\left|\right)=\frac{1}{1+{\left|R_{\max }\right|}^2/k^2},$

其中k是正的比例常数；

(5)根据步骤(4)构造的分割能量函数模型对SAR图像I进行下一步分割：

5a)对总的分割能量函数模型应用最速下降法，得到梯度下降流方程

$\frac{\partial \phi }{\partial t}=\delta \left(\phi \right)(-{\lambda }_1{(p_1-p_1^E)}^2+{\lambda }_2{(p_2-p_2^E)}^2$

$+\mathrm{\mu \delta }\left(\phi \right)\mathrm{div}(g\left(\right|R_{\max }\left|\frac{▿\phi }{|▿\phi |}\right)+v(\mathrm{\Delta \phi }-\mathrm{div}\left(\frac{▿\phi }{|▿\phi |}\right))$

其中δ(·)是Dirac函数，Δ是拉普拉子算子；

5b)对梯度下降流方程离散化，得到如下表达式：

$\frac{{\phi }^{n+1}-{\phi }^n}{\mathrm{\Delta t}}=\delta \left(\phi \right){(-{\lambda }_1(p_1-p_1^E)}^2+{\lambda }_2{(p_2-p_2^E)}^2$

$+\mathrm{\mu \delta }\left(\phi \right)\mathrm{div}(g\left(\right|R_{\max }\left|\frac{▿\phi }{|▿\phi |}\right)+v(\mathrm{\Delta \phi }-\mathrm{div}\left(\frac{▿\phi }{|▿\phi |}\right))$

其中φⁿ⁺¹代表第n+1次迭代后的水平集函数，φⁿ代表第n次迭代后的水平集函数，初始化的水平集函数φ是n＝0时的φⁿ，是迭代步长；

5c)根据步骤5b)求得新的水平集函数φⁿ⁺¹，由φⁿ⁺¹的正负值即得到新的分割区域和

5d)判断水平集函数是否收敛且达到最大的迭代次数，若不满足则转到步骤(3)，用和替代Ω₁和Ω₂继续迭代，否则停止迭代，得到的和即是最终的分割结果。

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

1、本发明将边缘和区域结合进水平集模型，既包含了图像的区域信息，又包含了重要的边缘信息。使分割方法不仅对图像边缘具有很好的局域化效果，又提高了对物体边缘定位的精确性，同时还提高了对边缘模糊和灰度不均匀的图像的分割效果和鲁棒性。

2、本发明结合了SAR图像相干斑噪声的统计模型，直接构造了水平集函数的能量模型，使分割模型在不需要对相干斑进行预处理的情况下，能够对受相干斑影响的SAR图像得到相对准确的分割结果。

3、本发明应用概率密度差异作为SAR图像的全局的区域信息，使总的能量函数模型的区域能量项构造更加合理。

4、本发明应用了对SAR图像具有很好的边缘定位性能的ROEWA边缘检测算子得到的边缘模值代替图像的梯度信息，从而使算法更好的使用于SAR图像分割。

经过与其它一些经典的水平集分割方法的对比，验证了本发明对SAR图像分割可以得到更好的分割效果。

附图说明

图1是本发明的流程图；

图2是本发明中待分割的四副SAR图像；

图3是采用经典的GAC方法对图2的SAR图像的分割结果；

图4是采用经典的CV方法对图2的SAR图像的分割结果；

图5是本发明对图2的SAR图像的分割结果。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤一，对待分割SAR图像I应用ROEWA边缘检测算子进行边缘检测，得到边缘强度模值|R_max|。

实现该步骤的具体过程如下：

(1a)分别构造因果滤波器f₁和非因果滤波器f₂如下：

f₁(z)＝ab^zH(z)，f₂(z)＝ab^-zH(-z)；

其中z是函数的自变量，a和b满足a＝1-b，0＜b＜e^-a＜1，H(·)是Heaviside函数。

(1b)构造指数平滑滤波器f如下：

$f\left(z\right)=\frac{1}{1+b}{f}_1\left(z\right)+\frac{b}{1+b}{f}_2(z-1),$其中z＝1，2，......N，N＞1。

(1c)应用构造的滤波器f，f₁，f₂对SAR图像I卷积求指数加权均值μ_x1，μ_x2，μ_y1，μ_y2，从而求得梯度模值，方法如下：

$\left(\begin{array}{c}{\mu }_{x1}=f_1\left(x\right)*(f\left(y\right)·I(x,y))\\ {\mu }_{x2}=f_2\left(x\right)*(f\left(y\right)·I(x,y))\\ {\mu }_{y1}=f_1\left(y\right)·(f\left(x\right)*I(x,y))\\ {\mu }_{y2}=f_2\left(y\right)·(f\left(x\right)*I(x,y))\end{array}\right)$

其中：x是水平方向的坐标变量，y是垂直方向的坐标变量。

定义：$\left(\begin{array}{c}{R}_{x\max }(x,y)=\max \{\frac{{\mu }_{x1}(x-1,y)}{{\mu }_{x2}(x+1,y)},\frac{{\mu }_{x2}(x+1,y)}{{\mu }_{x1}(x-1,y)}\}\\ R_{y\max }(x,y)=\max \{\frac{{\mu }_{y1}(x,y-1)}{{\mu }_{y2}(x,y+1)},\frac{{\mu }_{y2}(x,.y+1)}{{\mu }_{y1}(x,y-1)}\}\end{array}\right).$

类似梯度模的定义，定义ROEWA边缘检测算子边缘强度模值如下：

$\left|R_{\max }(x,y)\right|=\sqrt{R_{x\max }^2(x,y)+R_{y\max }^2(x,y)}.$

步骤二，将水平集φ(x，y)函数初始化成圆形的符号距离函数形式，根据水平集函数值的正负，将SAR图像分成两个区域Ω₁和Ω₂。

步骤三，计算区域Ω₁和Ω₂的强度均值c₁和c₂。

区域Ω₁和Ω₂的强度均值c₁和c₂表达式如下：

$c_i=\frac{{\int \int }_{{\Omega }_i}I(x,y)\mathrm{dxdy}}{{\int \int }_{{\Omega }_i}\mathrm{dxdy}},$

其中i＝1，2。

步骤四，根据c₁和c₂结合Gamma概率密度函数分别求出两个区域Ω₁和Ω₂的估计概率密度和

区域Ω₁和Ω₂的估计概率密度和表达式如下：

$p_i^E\left(I\right)=\frac{L^L}{c_i\Gamma \left(L\right)}{\left(\frac{I}{c_i}\right)}^{L-1}{e}^{-\mathrm{LI}/c_i},$

其中i＝1，2，L是SAR图像的视数。

步骤五，计算水平集分成的两个区域Ω₁和Ω₂的实际概率密度p₁和p₂。

区域Ω₁和Ω₂的实际概率密度p₁和p₂表达式如下：

$p_i\left(I\right)=\frac{{\int \int }_{{\Omega }_i}\delta (I(x_1,y_1)-I){\mathrm{dx}}_1{\mathrm{dy}}_1}{{\int \int }_{{\Omega }_i}{\mathrm{dx}}_1{\mathrm{dy}}_1},$

其中，i＝1，2，(x₁，y₁)是图像中的一个像素坐标，像素坐标δ(·)是Dirac函数。

步骤六，结合前面五个步骤建立总的能量函数E_SAR。

结合基于边缘和区域的水平集模型，结合避免重新初始化的惩罚项，构造总的能量函数总的能量函数E_SAR表达式如下：

$E_{\mathrm{SAR}}=\underset{i=1}{\overset{2}{\Sigma }}{\int \int }_{{\Omega }_i}{\lambda }_i{(p_i^E-p_i)}^2+\mu {\int \int }_{\Omega }g\left(\right|R_{\max }\left|\right)|▿H\left(\phi \right)|+v{\int \int }_{\Omega }\frac{1}{2}{\left(\right|▿\phi |-1)}^2$

其中Ω是整个图像区域，即Ω＝Ω₁+Ω₂，λ₁，λ₂，μ，v是各部分能量项的权重，λ₁＞0，λ₂＞0，μ＞0，v＞0；

g(x，y)是定义在边缘强度模值|R_max|上的指示函数，表达式是：

$g\left(\right|R_{\max }\left|\right)=\frac{1}{1+{\left|R_{\max }\right|}^2/k^2},$其中k是以比例常数。

步骤七，对能量函数E_SAR应用变分法，得到梯度下降流方程并迭代更新水平集函数φ，得到新的分割区域Ω₁和Ω₂。

实现该步骤的具体过程如下：

(7a)对上一步骤构造的总的能量函数应用最速下降法，得到梯度下降流方程其表达式如下：

$\frac{\partial \phi }{\partial t}=\delta \left(\phi \right)(-{\lambda }_1{(p_1-p_1^E)}^2+{\lambda }_2{(p_2-p_2^E)}^2$

$+\mu *\delta \left(\phi \right)\mathrm{div}(g\left(\frac{▿\phi }{|▿\phi |}\right)+v(\mathrm{\Delta \phi }-\mathrm{div}\left(\frac{▿\phi }{|▿\phi |}\right)).$

(7b)对上一步骤求的梯度下降流方程离散化，得到如下表达式：

$\frac{{\phi }^{n+1}-{\phi }^n}{\mathrm{\Delta t}}=\delta \left(\phi \right){(-{\lambda }_1(p_1-p_1^E)}^2+{\lambda }_2{(p_2-p_2^E)}^2$

$+\mu *\delta \left(\phi \right)\mathrm{div}(g\left(\frac{▿\phi }{|▿\phi |}\right)+v(\mathrm{\Delta \phi }-\mathrm{div}\left(\frac{▿\phi }{|▿\phi |}\right)),$

其中φⁿ⁺¹代表第n+1次迭代后的水平集函数，φⁿ代表第n次迭代后的水平集函数，初始化的水平集函数可看作φ⁰，是迭代步长。

根据离散化后的可得到新的水平集函数φⁿ⁺¹，并得到新的区域Ω₁和Ω₂。

步骤八，判断是否收敛。

本步骤的主要目的是根据水平集的迭代情况，判断是否收敛且是否达到最大迭代次数。如果梯度下降流方程没有达到稳定状态且没有达到最大的迭代次数，则转到步骤(3)继续迭代，否则停止迭代得到最终的分割结果。

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

1、仿真条件

硬件平台为：Intel Core2 Duo CPU E6550@2.33GHZ、2GB RAM。

软件平台为：MATLAB 7.5。

2、仿真内容与结果

本发明将图2所示的(a)、(b)、(c)、(d)四副图像作为测试图像，用本发明所提出的SAR图像分割方法与现有的GAC方法和CV方法进行比较。

图3是应用现有的GAC方法对图2的分割结果图，其中图3(a)式对图2(a)的分割结果，图3(b)是对图2(b)的分割结果，图3(c)是对图2(c)的分割结果，图3(d)是对图2(d)的分割结果，从如图3(a)可以看出图3(a)有一部分虽然梯度比较明显但出现过分割，并且由图3(b)、3(c)、3(d)可看到GAC方法对初始化比较敏感，对区域比较丰富的SAR图像只能检测出部分轮廓，而不能得到完整的分割结果。

图4是应用现有的CV方法对图2的分割结果图，其中图4(a)是对图2(a)的分割结果，图4(b)是对图2(b)的分割结果，图4(c)是对图2(c)的分割结果，图4(d)是对图2(d)的分割结果，从图4可以看出现有的CV方法检测出了比较多的虚假边缘，对噪声比较敏感。

图5是应用本发明的方法对图2的分割结果图，其中图5(a)式对图2(a)的分割结果，图5(b)式对图2(b)的分割结果，图5(c)式对图2(c)的分割结果，图5(d)式对图2(d)的分割结果。对比图3(a)和5(a)可以看出，本发明不会出现过分割现象；对比3(b)、5(b)，3(c)、5(c)和3(d)、5(d)可以看出，本发明对初始化不再敏感；对比4(a)、5(a)，4(b)、5(b)，4(c)、5(c)和4(d)、5(d)可以看出，本发明的方法对噪声具有抑制作用，增强了真实的边缘的定位性能。