CN102426699B

CN102426699B - Level set synthetic aperture radar (SAR) image segmentation method based on edge and regional information

Info

Publication number: CN102426699B
Application number: CN 201110346314
Authority: CN
Inventors: 侯彪; 焦李成; 刘娜娜; 王爽; 刘芳; 尚荣华
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2011-11-04
Filing date: 2011-11-04
Publication date: 2013-08-14
Anticipated expiration: 2031-11-04
Also published as: CN102426699A

Abstract

The invention discloses a level set SAR image segmentation method based on the combination of edge and region information, which mainly solves the problems that the existing level set method is difficult to segment the SAR image with blurred edges and the real edge of the SAR image is inaccurately located. The implementation steps include: firstly, the exponential weighted mean ratio operator is used to detect the edge of the SAR image, and the edge intensity modulus |R _max | is obtained, and the edge energy item is constructed; secondly, the level set function φ is initialized, and the SAR image is divided into inner and outer regions Ω ₁ and Ω ₂ , and calculate the intensity mean values c ₁ and c ₂ of these two regions; calculate the estimated probability densities p ₁ and p ₂ of the regions Ω ₁ and Ω ₂ according to c ₁ and c ₂ again, and construct the regional energy term ;Finally, add the corrected energy item to avoid reinitialization, construct the total energy function E _SAR , and apply the variational method to find the gradient descent flow equation, and update the level set φ to obtain a new segmentation area

and

Experimental results show that the segmentation method realized by the present invention can obtain relatively ideal segmentation effect, and can be used for edge detection and target recognition of SAR images.

Description

Level Set SAR Image Segmentation Method Based on Edge and Region Information

技术领域 technical field

本发明属于图像处理技术领域，涉及SAR图像分割，具体地说是一种基于边缘和区域信息相结合的水平集SAR图像分割方法，可用于SAR图像的分割、边缘检测和目标识别。The invention belongs to the technical field of image processing and relates to SAR image segmentation, in particular to a level set SAR image segmentation method based on the combination of edge and region information, which can be used for SAR image segmentation, edge detection and target recognition.

背景技术 Background technique

合成孔径雷达SAR是一种高分辨率主动式雷达，具有全天候、全天时、分辨率高、可侧视成像等优点，可应用于军事、农业、导航、地理监视等诸多领域。SAR图像被广泛的应用在目标检测领域，而SAR图像分割则是从图像处理到图像分析的重要步骤，是目标分类和识别的基础。由于SAR是一种相干成像系统，SAR图像在本质上是对目标的电磁散射特性和结构特性的反应，其成像效果很大程度上依赖于雷达参数和地域电磁参数。由于SAR成像系统的特殊性，使SAR图像的信息表达方式和光学图像有很大差异，并且会受到相干斑噪声及阴影等许多几何特征的影响，正因为这些因素使适用于光学图像的分割方法对SAR图像不再适用。Synthetic Aperture Radar (SAR) is a high-resolution active radar with the advantages of all-weather, all-time, high resolution, and side-view imaging. It can be used in many fields such as military, agriculture, navigation, and geographical surveillance. SAR images are widely used in the field of target detection, and SAR image segmentation is an important step from image processing to image analysis, and is the basis of target classification and recognition. Since SAR is a coherent imaging system, the SAR image is essentially a response to the electromagnetic scattering characteristics and structural characteristics of the target, and its imaging effect depends largely on radar parameters and regional electromagnetic parameters. Due to the particularity of the SAR imaging system, the information expression of SAR images is very different from that of optical images, and will be affected by many geometric features such as coherent speckle noise and shadows. It is precisely because of these factors that the segmentation method suitable for optical images No longer applicable for SAR images.

总体上说常规的SAR分割方法大多数还是用传统的分割模型对SAR图像进行分割，如一些基于阈值的分割、基于边缘的分割和基于区域生长的分割方法等，但是其中的很多分割方法容易受到SAR图像中相干斑噪声的影响，所以此类方法都需要先通过预处理抑制相干斑噪声。这些预处理方法在抑制相干斑噪声的同时，不可避免地损失了很多边缘与目标信息，影响了分割效果。Generally speaking, most conventional SAR segmentation methods still use traditional segmentation models to segment SAR images, such as some threshold-based segmentation, edge-based segmentation, and region-growing segmentation methods, etc., but many of them are vulnerable to The influence of coherent speckle noise in SAR images, so such methods need to suppress coherent speckle noise through preprocessing. These preprocessing methods inevitably lose a lot of edge and target information while suppressing speckle noise, which affects the segmentation effect.

水平集图像分割方法是重要的图像处理方法之一。水平集方法的优点是可以适应拓扑结构的变化，而且算法稳定性较高。利用水平集方法研究SAR图像分割问题，通过建立合适的能量泛函，可以在能量泛函的定义中引入图像区域信息，在不需要相干斑预处理的情况下，对于受相干斑噪声影响的SAR图像可以获得比较准确的分割结果。Level set image segmentation method is one of the important image processing methods. The advantage of the level set method is that it can adapt to the change of topology, and the algorithm has high stability. The level set method is used to study the SAR image segmentation problem. By establishing a suitable energy functional, the image region information can be introduced into the definition of the energy functional. Without the need for speckle preprocessing, for SAR affected by speckle noise The image can obtain more accurate segmentation results.

经典的水平集方法分基于边缘信息和基于区域信息的两种方法。测地线活动轮廓模型GAC模型是基于边缘信息的水平集方法中应用最广的一种，它能在不附加任何外界控制条件下，自由处理曲线运动时的拓扑结构变化，但是该模型仅利用图像的梯度信息，由于图像运动中的边缘并非都是理想的阶梯边缘，难以分割出边缘模糊的图像中的同质区域，另外当图像中的对象有较深的凹陷边界时，GAC模型可能使演化曲线停止在某一局部极小值状态，并不与对象的边界相一致。Chan和Vese提出的CV模型是一种基于区域的水平集分割方法，作为一种能有效提高曲线演化拓扑自适应能力的区域水平集分割模型，它利用了图像的全局区域信息。但是其同样也存在不足之处，该模型的内部能量项只是保证零水平集曲线光滑，而没有考虑水平集函数本身所固有的内在性质，在某些应用中还需要对水平集函数进行重新初始化，以使它接近符号距离函数，保证数值解法的稳定性，另外，模型中仅仅利用了图像的区域信息，而另一个重要的图像的边缘信息没有很好地利用，在一些分割应用中尤其是对灰度不均匀的图像可能产生图像边缘定位不准确的缺陷。The classic level set method is divided into two methods based on edge information and based on region information. The geodesic active contour model (GAC) model is the most widely used level set method based on edge information. It can freely deal with the topological structure changes during the curve motion without any external control conditions. However, the model only uses The gradient information of the image, because the edges in the image motion are not all ideal step edges, it is difficult to segment the homogeneous area in the image with blurred edges. In addition, when the object in the image has a deep concave boundary, the GAC model may make the The evolution curve stops at a certain local minimum state, which does not coincide with the boundary of the object. The CV model proposed by Chan and Vese is a region-based level set segmentation method. As a regional level set segmentation model that can effectively improve the adaptive ability of curve evolution topology, it utilizes the global region information of the image. However, it also has shortcomings. The internal energy term of the model only ensures the smoothness of the zero level set curve, without considering the intrinsic properties of the level set function itself. In some applications, the level set function needs to be reinitialized. , so that it is close to the signed distance function and ensures the stability of the numerical solution. In addition, only the area information of the image is used in the model, and the edge information of another important image is not well utilized. In some segmentation applications, especially For images with uneven gray levels, the defect of inaccurate image edge positioning may occur.

虽然水平集分割方法已经在光学图像和医学图像分割中获得了巨大的成功，但是在SAR图像分割领域，研究的还比较少。目前国际上知名的有法国的Refregier等人和美国的Chen等人所领导的两个机构，并且两个机构各有偏重，前者偏重研究测地线活动轮廓GAC模型在SAR图像分割和边缘检测方面的应用，且此模型依赖于图像梯度定义边缘检测算子，由于GAC模型中的梯度信息本身就不适用于SAR图像，所以产生许多错误分割，后者主要研究几何活动轮廓CV模型对多区域SAR图像的分割，由于CV模型是一个基于加性高斯噪声的分割模型，若直接将该模型应用于SAR图像，由于乘性相干斑噪声的存在，得到的检测结果十分不理想，且它仅利用了SAR图像的区域信息，而SAR图像很重要的边缘信息未被充分利用，无法正确地分割SAR图像。Although the level set segmentation method has achieved great success in optical image and medical image segmentation, but in the field of SAR image segmentation, there is still relatively little research. At present, there are two well-known institutions in the world, led by Refregier et al. in France and Chen et al. in the United States, and the two institutions have their own emphasis. The former focuses on the study of the geodesic active contour GAC model in SAR image segmentation and edge detection. The application of this model relies on the image gradient to define the edge detection operator. Since the gradient information in the GAC model itself is not suitable for SAR images, many error segmentations are generated. The latter mainly studies the geometric active contour CV model for multi-region SAR. Image segmentation, since the CV model is a segmentation model based on additive Gaussian noise, if the model is directly applied to SAR images, due to the existence of multiplicative coherent speckle noise, the detection results obtained are very unsatisfactory, and it only uses The area information of the SAR image, and the important edge information of the SAR image are not fully utilized, and the SAR image cannot be segmented correctly.

发明内容 Contents of the invention

本发明的目的在于针对上述已有技术的不足，提供一种基于边缘和区域相结合的水平集SAR图像分割方法，以实现对边缘模糊的SAR图像区域的准确分割，提高边缘的定位性能，从而提高SAR图像的分割的质量。The purpose of the present invention is to address the deficiencies of the above-mentioned prior art, to provide a level set SAR image segmentation method based on the combination of edges and regions, so as to realize accurate segmentation of SAR image regions with blurred edges and improve the positioning performance of edges, thereby Improve the quality of segmentation of SAR images.

实现本发明目的技术思路是：根据SAR图像的统计特性，应用Gamma分布作为SAR图像的概率密度函数，并且计算由水平集划分得到的两个区域的概率密度的估计值，以各估计值的对数作为区域信息，同时应用对SAR图像有很好边缘定位性能的指数加权均值比率的边缘检测算子的检测结果代替梯度值作为边缘信息，最后加入避免重新初始化的修正项，将这三部分信息结合进行能量函数的建模，从而达到理想的分割结果。其具体的实现过程包括如下：The technical thought of realizing the object of the present invention is: according to the statistical characteristic of SAR image, apply Gamma distribution as the probability density function of SAR image, and calculate the estimated value of the probability density of two regions obtained by level set division, with the pair of each estimated value At the same time, the detection result of the exponentially weighted mean ratio edge detection operator that has good edge positioning performance for SAR images is used as the edge information instead of the gradient value as the edge information. Finally, a correction item to avoid re-initialization is added, and the three parts of information Combined with the modeling of the energy function, the ideal segmentation result can be achieved. Its specific implementation process includes the following:

(1)对待分割SAR图像I应用指数加权均值比率的边缘检测算子进行边缘检测，得到边缘强度模值|R_max|；(1) Apply an edge detection operator with an exponentially weighted mean ratio to the SAR image I to be segmented to perform edge detection, and obtain the edge intensity modulus |R _max |;

(2)将水平集函数φ初始化成符号距离函数形式，根据水平集函数值的正负，将SAR图像分割成两个区域Ω₁和Ω₂；(2) The level set function φ is initialized into a signed distance function form, and the SAR image is divided into two regions Ω ₁ and Ω ₂ according to the positive or negative of the level set function value;

(3)根据两个区域Ω₁和Ω₂，计算其对应的估计概率密度p₁和p₂：(3) According to the two regions Ω ₁ and Ω ₂ , calculate their corresponding estimated probability densities p ₁ and p ₂ :

3a)计算区域Ω₁和Ω₂的强度均值c_i：3a) Compute the intensity mean c _i of regions Ω ₁ and Ω ₂ :

其中i＝1，2，(x，y)是图像坐标；

where i=1, 2, (x, y) are image coordinates;

3b)根据强度均值c_i，计算区域Ω₁和Ω₂的估计概率密度p₁和p₂：3b) Calculate the estimated probability densities p ₁ and p ₂ for regions Ω ₁ and Ω ₂ from the intensity mean c _i :

${p p}_{i i} = = \frac{{L L}^{L L}}{{c c}_{i i} Γ Γ ((L L))} {((\frac{I I}{{c c}_{i i}}))}^{L L - - 11} {e e}^{- - LI LI / / {c c}_{i i}},,$

其中L是SAR图像的视数，Γ(·)是Gamma函数，i＝1，2；Wherein L is the view number of the SAR image, Γ(·) is the Gamma function, i=1,2;

(4)结合步骤(1)-步骤(3)，构造总的分割能量函数模型E_SAR：(4) In conjunction with step (1)-step (3), construct the total segmentation energy function model E _SAR :

${E E.}_{SAR SAR} = = - - {Σ Σ}_{i i = = 11}^{22} {&Integral; &Integral; &Integral; &Integral;}_{{Ω Ω}_{i i}} {λ λ}_{i i} log log (({p p}_{i i})) + + μ μ {&Integral; &Integral; &Integral; &Integral;}_{Ω Ω} g g ((| | {R R}_{max max} | |)) | | φH φH ((φ φ)) | | + + v v {&Integral; &Integral; &Integral; &Integral;}_{Ω Ω} \frac{11}{22} {((| | &dtri; &dtri; φ φ | | - - 11))}^{22},,$

其中Ω是整个图像区域，即Ω＝Ω₁+Ω₂，

是区域能量项，

是边缘能量项，

是修正项能量项，λ_i是区域能量项的权值，i＝1，2，μ是边缘能量项的权值，v是修正项的权值，λ_i＞0，μ＞0，v＞0，是表示对φ求梯度，H(φ)是Heaviside函数，

表示对H(φ)求梯度，g(|R_max|)是定义在边缘强度模值|R_max|上的指示函数，表达式如下：Where Ω is the entire image area, ie Ω=Ω ₁ +Ω ₂ ,

is the area energy term,

is the marginal energy term,

is the energy item of the correction item, λ _i is the weight of the area energy item, i=1, 2, μ is the weight of the edge energy item, v is the weight of the correction item, λ _i >0, μ>0, v> 0, It means to find the gradient for φ, H(φ) is the Heaviside function,

Represents the gradient of H(φ), g(|R _max |) is an indicator function defined on the edge strength modulus |R _max |, the expression is as follows:

$g (| R_{\max} |) = \frac{1}{1 + {| R_{\max} |}^{2} / k^{2}},$ 其中k是正的比例常数； $g (| R_{\max} |) = \frac{1}{1 + {| R_{\max} |}^{2} / k^{2}},$ where k is a positive proportionality constant;

(5)根据步骤(4)构造的分割能量函数模型对SAR图像I进行分割：(5) segment the SAR image I according to the segmentation energy function model constructed in step (4):

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

5a) Apply the variational method to the total segmentation energy function model to obtain the gradient descent flow equation

$\frac{&PartialD; &PartialD; φ φ}{&PartialD; &PartialD; t t} = = δ δ ((φ φ)) ((- - {λ λ}_{11} log log (({p p}_{11})) + + {λ λ}_{22} log log (({p p}_{22}))))$

$+ + μδ μδ ((φ φ)) div div ((g g ((| | {R R}_{max max} | | \frac{&dtri; &dtri; φ φ}{| | &dtri; &dtri; φ φ | |})) + + v v ((Δφ Δφ - - div div ((\frac{&dtri; &dtri; φ φ}{| | &dtri; &dtri; φ φ | |})))),,$

其中δ(φ)是Dirac函数，Δ是拉普拉斯算子；Where δ(φ) is the Dirac function, Δ is the Laplacian operator;

5b)对梯度下降流方程

离散化，得到如下表达式：5b) For the gradient descent flow equation

Discretization, the following expression is obtained:

$\frac{{φ φ}^{n no + + 11} - - {φ φ}^{n no}}{Δt Δt} = = δ δ ((φ φ)) ((- - {λ λ}_{11} log log (({p p}_{11})) + + {λ λ}_{22} log log (({p p}_{22}))))$

其中φⁿ⁺¹代表第n+1次迭代后的水平集函数，φⁿ代表第n次迭代后的水平集函数，Δt是迭代步长；Where φ ⁿ⁺¹ represents the level set function after the n+1th iteration, φ ⁿ represents the level set function after the nth iteration, and Δt is the iteration step size;

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

和

5c) Obtain a new level set function φ ⁿ⁺¹ according to step 5b), and obtain a new segmented region from the positive and negative values of φ ⁿ⁺¹

and

5d)判断水平集函数是否收敛且达到最大的迭代次数100次，若不满足则转到步骤(3)，用

和

替代Ω₁和Ω₂继续迭代，否则停止迭代，得到的和即是最终的分割结果。5d) Judging whether the level set function converges and reaches the maximum number of iterations of 100 times, if not, go to step (3), use

and

Substitute Ω ₁ and Ω ₂ to continue iterating, otherwise stop iterating, and get and That is the final segmentation result.

本发明与现有的技术相比具有以下优点：Compared with the prior art, the present invention has the following advantages:

1、本发明将边缘和区域信息结合进水平集模型，既包含了SAR图像的区域信息，又包含了重要的边缘信息，不仅对图像边缘具有很好的分割局域化效果，又提高了对物体边缘定位的精确性，同时还提高了对边缘模糊和灰度不均匀的SAR图像的分割效果和鲁棒性。1. The present invention combines the edge and region information into the level set model, which not only includes the region information of the SAR image, but also includes important edge information, which not only has a good segmentation and localization effect on the edge of the image, but also improves the accuracy of the image. The accuracy of object edge positioning is also improved, and the segmentation effect and robustness of SAR images with blurred edges and uneven gray levels are also improved.

2、本发明结合了适合SAR图像相干斑噪声的Gamma统计模型，直接构造了水平集能量函数的区域项能量模型，使分割模型在不需要对相干斑噪声进行预处理的情况下，能够对受相干斑噪声影响的SAR图像得到相对准确的分割结果。2. The present invention combines the Gamma statistical model suitable for SAR image coherent speckle noise, and directly constructs the regional item energy model of the level set energy function, so that the segmentation model can be used without preprocessing the coherent speckle noise. SAR images affected by coherent speckle noise get relatively accurate segmentation results.

3、本发明应用了对SAR图像具有很好的边缘定位性能的指数加权均值比率的边缘检测算子得到的边缘强度模值代替图像的梯度信息，从而使算法更好的使用于SAR图像分割。3. The present invention uses the edge intensity modulus obtained by the edge detection operator of the exponentially weighted mean ratio that has good edge location performance for SAR images instead of the gradient information of the image, so that the algorithm can be better used in SAR image segmentation.

经过与现有的水平集分割方法的仿真对比，验证了本发明对SAR图像分割能得到更好的分割效果。Through simulation comparison with the existing level set segmentation method, it is verified that the present invention can obtain better segmentation effect for SAR image segmentation.

附图说明 Description of drawings

图1是本发明的流程图；Fig. 1 is a flow chart of the present invention;

图2是第一幅测试图像；Figure 2 is the first test image;

图3是用本发明与已有方法对图2的分割效果对比图；Fig. 3 is a comparison diagram of the segmentation effect of Fig. 2 with the present invention and existing method;

图4是第二幅测试图像；Figure 4 is the second test image;

图5是用本发明与已有方法对图4的分割效果对比图；Fig. 5 is a comparison diagram of the segmentation effect of Fig. 4 with the present invention and existing methods;

具体实施方式Detailed ways

参照图1，本发明的具体实现步骤如下：With reference to Fig. 1, the concrete realization steps of the present invention are as follows:

步骤一，将待分割SAR图像I应用指数加权均值比率的边缘检测算子进行边缘检测，得到边缘强度模值|R_max|。Step 1: Apply an edge detection operator with an exponentially weighted mean ratio to the SAR image I to be segmented for edge detection, and obtain the edge intensity modulus |R _max |.

实现该步骤的具体过程如下：The specific process to realize this step is as follows:

1a)分别构造因果滤波器f₁和非因果滤波器f₂的函数表达式：1a) Construct the functional expressions of the causal filter f ₁ and the non-causal filter f ₂ respectively:

f₁(z)＝ab^zH(z)，f₂(z)＝ab^-zH(-z)，f ₁ (z) = ab ^z H(z), f ₂ (z) = ab ^-z H(-z),

其中z是函数的自变量，且z＝1，2，......N，N为N＞1的正整数；a和b均为常量，且满足0＜b＜e^-a＜1，H(·)是Heaviside函数；Where z is the independent variable of the function, and z=1, 2, ... N, N is a positive integer of N>1; both a and b are constants, and satisfy 0<b<e ^-a <1 , H(·) is the Heaviside function;

1b)根据滤波器f₁和f₂构造指数平滑滤波器f的函数表达式：1b) Construct the functional expression of the exponential smoothing filter f according to the filters f ₁ and f ₂ :

$f f ((z z)) = = \frac{11}{11 + + b b} {f f}_{11} ((z z)) + + \frac{b b}{11 + + b b} {f f}_{22} ((z z - - 11)),,$

1c)根据构造的滤波器f，f₁，f₂，计算滤波器f₁在水平方向上的指数加权均值μ_x1，滤波器f₂在水平方向上的指数加权均值μ_x2，滤波器f₁在垂直方向上的指数加权均值μ_y1，滤波器f₂在垂直方向上的指数加权均值μ_y2，各指数加权均值μ_x1，μ_x2，μ_y1，μ_y2：1c) According to the constructed filters f, f ₁ , f ₂ , calculate the exponentially weighted mean μ _x1 of the filter f ₁ in the horizontal direction, the exponentially weighted mean μ _x2 of the filter f ₂ in the horizontal direction, filter f ₁ The exponentially weighted mean value μ _y1 in the vertical direction, the exponentially weighted mean value μ _y2 of the filter f ₂ in the vertical direction, each exponentially weighted mean value μ _x1 , μ _x2 , μ _y1 , μ _y2 :

$\{\begin{matrix} {μ μ}_{x x 11} = = {f f}_{11} ((x x)) * * ((f f ((y the y)) \cdot &Center Dot; I I ((x x,, y the y)))) \\ {μ μ}_{x x 22} = = {f f}_{22} ((x x)) * * ((f f ((y the y)) \cdot &Center Dot; I I ((x x,, y the y)))) \\ {μ μ}_{y the y 11} = = {f f}_{11} ((y the y)) \cdot &Center Dot; ((f f ((x x)) * * I I ((x x,, y the y)))) \\ {μ μ}_{y the y 22} = = {f f}_{22} ((y the y)) \cdot \cdot ((f f ((x x)) * * I I ((x x,, y the y)))) \end{matrix},,$

其中：x是水平方向的坐标变量，y是垂直方向的坐标变量，*代表水平方向的卷积，·代表垂直方向的卷积；Among them: x is the coordinate variable in the horizontal direction, y is the coordinate variable in the vertical direction, * represents the convolution in the horizontal direction, · represents the convolution in the vertical direction;

1d)应用1c)的结果，求水平方向的强度模值R_xmax(x，y)和垂直方向的强度模值R_ymax(x，y)：1d) Apply the results of 1c) to find the intensity modulus R _xmax (x, y) in the horizontal direction and the intensity modulus R _ymax (x, y) in the vertical direction:

$\{\begin{matrix} {R R}_{x x max max} ((x x,, y the y)) = = max max {{\frac{{μ μ}_{x x 11} ((x x - - 11,, y the y))}{{μ μ}_{x x 22} ((x x + + 11,, y the y))},, \frac{{μ μ}_{x x 22} ((x x + + 11,, y the y))}{{μ μ}_{x x 11} ((x x - - 11,, y the y))}}} \\ {R R}_{y the y max max} ((x x,, y the y)) = = max max {{\frac{{μ μ}_{y the y 11} ((x x,, y the y - - 11))}{{μ μ}_{y the y 22} ((x x,, y the y + + 11))},, \frac{{μ μ}_{y the y 22} ((x x,, y the y + + 11))}{{μ μ}_{y the y 11} ((x x,, y the y - - 11))}}} \end{matrix},,$

1e)根据R_xmax(x，y)和R_ymax(x，y)求边缘强度模值|R_max|：1e) According to R _xmax (x, y) and R _ymax (x, y), calculate the edge strength modulus |R _max |:

$| | {R R}_{max max} ((x x,, y the y)) | | = = \sqrt{{R R}_{x x max max}^{22} ((x x,, y the y)) + + {R R}_{y the y max max}^{22} ((x x,, y the y))} . .$

步骤二，将水平集函数φ初始化成符号距离函数形式，根据水平集函数值的正负，将SAR图像分成两个区域Ω₁和Ω₂。Step 2: Initialize the level set function φ into a signed distance function form, and divide the SAR image into two regions Ω ₁ and Ω ₂ according to the positive or negative value of the level set function.

实现该步骤的具体过程表述如下：用计算机工具在待分割图像I上做一个矩形，此矩形方程作为初始化水平集函数φ，当φ＞0时，表示矩形外部区域Ω₁，当φ＜0时，表示矩形内部区域Ω₂，因此，根据水平集函数值的正负，将SAR图像分成两个区域Ω₁和Ω₂。The specific process of realizing this step is expressed as follows: use computer tools to make a rectangle on the image I to be segmented, and the equation of this rectangle is used as the initial level set function φ. When φ>0, it represents the outer area of the rectangle Ω ₁ . When φ<0 , represents the rectangular inner area Ω ₂ , therefore, according to the positive or negative value of the level set function, the SAR image is divided into two areas Ω ₁ and Ω ₂ .

步骤三，根据现有的求SAR图像强度均值的公式，计算区域Ω₁的强度均值c₁和区域Ω₂的强度均值c₂，表达式为：Step 3, according to the existing formula for calculating the mean value of SAR image intensity, calculate the mean intensity value c ₁ of the region Ω ₁ and the mean value c ₂ of the region Ω ₂ , the expression is:

$\{\begin{matrix} c_{1} = \frac{{&Integral; &Integral;}_{Ω_{1}} I (x, y) dxdy}{{&Integral; &Integral;}_{Ω_{1}} dxdy} \\ c_{2} = \frac{{&Integral; &Integral;}_{Ω_{2}} I (x, y) dxdy}{{&Integral; &Integral;}_{Ω_{2}} dxdy} \end{matrix},$ 其中(x，y)是图像坐标， $\{\begin{matrix} c_{1} = \frac{{&Integral; &Integral;}_{Ω_{1}} I (x, the y) dxdy}{{&Integral; &Integral;}_{Ω_{1}} dxdy} \\ c_{2} = \frac{{&Integral; &Integral;}_{Ω_{2}} I (x, the y) dxdy}{{&Integral; &Integral;}_{Ω_{2}} dxdy} \end{matrix},$ where (x, y) are the image coordinates,

步骤四，根据步骤三求得的强度均值c_i和现有的SAR图像概率密度函数的公式，计算区域Ω₁和Ω₂的估计概率密度p₁和p₂：Step 4: Calculate the estimated probability densities p 1 and _p 2 of regions Ω ₁ and Ω ₂ according to the intensity mean c _i obtained in step ₃ and the formula of the existing SAR image probability density function:

$\{\begin{matrix} {p p}_{11} ((I I)) = = \frac{{L L}^{L L}}{{c c}_{11} Γ Γ ((L L))} {((\frac{I I}{{c c}_{11}}))}^{L L - - 11} {e e}^{- - LI LI / / {c c}_{11}} \\ {p p}_{22} ((I I)) = = \frac{{L L}^{L L}}{{c c}_{22} Γ Γ ((L L))} {((\frac{I I}{{c c}_{22}}))}^{L L - - 11} {e e}^{- - LI LI / / {c c}_{22}} \end{matrix},,$

其中L是SAR图像的视数，I为待分割SAR图像，Γ(·)是Gamma函数。Where L is the view number of the SAR image, I is the SAR image to be segmented, and Γ(·) is the Gamma function.

步骤五，结合前面四个步骤构造总的能量函数E_SAR Step five, combine the previous four steps to construct the total energy function E _SAR

5b)根据步骤(5a)计算的边缘指示函数，构造边缘能量项，其表达式为：5b) Construct the edge energy term according to the edge indicator function calculated in step (5a), and its expression is:

$μ μ {&Integral; &Integral; &Integral; &Integral;}_{Ω Ω} g g ((| | {R R}_{max max} | |)) | | &dtri; &dtri; H h ((φ φ)) | |,,$

5c)根据步骤四计算求得的区域Ω₁和Ω₂的估计概率密度p₁和p₂，构造区域能量项，其表达式为： $- {&Integral; &Integral;}_{Ω_{1}} λ_{1} \log (p_{1}) - {&Integral; &Integral;}_{Ω_{2}} λ_{2} \log (p_{2}),$ 5c) According to the estimated probability densities p ₁ and p ₂ of the regions Ω ₁ and Ω ₂ calculated and obtained in step 4, the regional energy item is constructed, and its expression is: $- {&Integral; &Integral;}_{Ω_{1}} λ_{1} \log (p_{1}) - {&Integral; &Integral;}_{Ω_{2}} λ_{2} \log (p_{2}),$

5d)结合步骤(5b)和(5c)求得的边缘能量项和区域能量项，再加入避免重新初始化的惩罚能量项，构造总的能量函数E_SAR，其表达式如下：5d) Combining the edge energy term and area energy term obtained in steps (5b) and (5c), and then adding the penalty energy term to avoid reinitialization, constructing the total energy function E _SAR , its expression is as follows:

${E E.}_{SAR SAR} = = - - {&Integral; &Integral; &Integral; &Integral;}_{{Ω Ω}_{11}} {λ λ}_{11} log log (({p p}_{11})) - - {&Integral; &Integral; &Integral; &Integral;}_{{Ω Ω}_{22}} {λ λ}_{22} log log (({p p}_{22})) + + μ μ {&Integral; &Integral; &Integral; &Integral;}_{Ω Ω} g g ((| | {R R}_{max max} | |)) | | &dtri; &dtri; H h ((φ φ)) | | + + v v {&Integral; &Integral; &Integral; &Integral;}_{Ω Ω} \frac{11}{22} {((| | &dtri; &dtri; φ φ | | - - 11))}^{22},,$

其中Ω是整个图像区域，即Ω＝Ω₁+Ω₂，λ₁和λ₂分别为内部区域Ω₁能量项和外部区域Ω₂能量项的权值，μ是边缘能量项的权值，v是修正项的权值，λ₁＞0，λ₂＞0，μ＞0，v＞0，是表示对φ求梯度，H(φ)是Heaviside函数，

表示对H(φ)求梯度。Where Ω is the entire image area, that is, Ω=Ω ₁ +Ω ₂ , λ ₁ and λ ₂ are the weights of the Ω ₁ energy item in the inner area and Ω ₂ energy item in the outer area respectively, μ is the weight of the edge energy item, v is the weight of the correction item, λ ₁ >0, λ ₂ >0, μ>0, v>0, It means to find the gradient for φ, H(φ) is the Heaviside function,

Represents the gradient of H(φ).

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

和

Step six, apply the variational method to the energy function E _SAR to obtain the gradient descent flow equation And update the level set function φ iteratively to get the final segmentation area

and

6a)对总的分割能量函数模型应用变分法，得到梯度下降流方程 6a) Apply the variational method to the total segmentation energy function model to obtain the gradient descent flow equation

6b)对梯度下降流方程

离散化，得到如下表达式：6b) For the gradient descent flow equation

Discretization, the following expression is obtained:

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

6c) Obtain a new level set function φ ⁿ⁺¹ according to step 6b), and obtain a new segmented region from the positive and negative values of φ ⁿ⁺¹ and

步骤七，判断是否收敛。Step 7, judging whether to converge.

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

和

替代Ω₁和Ω₂继续迭代，否则停止迭代，得到最终的区域和即是最终的分割结果。The main purpose of this step is to judge whether the level set function converges and reaches the maximum number of iterations according to the iteration of the level set. If the gradient descent flow equation does not reach a steady state and the maximum number of iterations is not reached, then go to step (3) ,use

and

Replace Ω ₁ and Ω ₂ to continue iterating, otherwise stop iterating and get the final area and That is the final segmentation result.

本发明的效果可以通过下面的仿真进一步说明：Effect of the present invention can be further illustrated by following simulation:

1、仿真条件1. Simulation conditions

硬件平台为：Intel Core2 Duo CPU E65502.33GHZ、2GB RAM。The hardware platform is: Intel Core2 Duo CPU E65502.33GHZ, 2GB RAM.

软件平台为：MATLAB 7.9。The software platform is: MATLAB 7.9.

2、仿真内容与结果2. Simulation content and results

仿真一：应用本发明和现有的侧地活动轮廓(GAC)方法和现有的CV方法分别对图2的测试图像进行分割实验，其分割结果如图3所示，其中：图3(a)是现有的GAC方法对图2的分割结果；图3(b)是现有的CV方法对图2的分割结果；图3(c)是本发明方法对图2的分割结果。Simulation one: apply the present invention and existing lateral ground active contour (GAC) method and existing CV method to carry out segmentation experiment respectively to the test image of Fig. 2, its segmentation result is as shown in Fig. 3, wherein: Fig. 3 (a ) is the segmentation result of the existing GAC method to Fig. 2; Fig. 3 (b) is the segmentation result of the existing CV method to Fig. 2; Fig. 3 (c) is the segmentation result of the inventive method to Fig. 2.

由图3(a)可看出，GAC方法对初始化比较敏感，对区域比较丰富的SAR图像只能检测出部分轮廓，而不能得到完整的分割结果。由图3(b)可看出，现有的CV方法检测出了比较多的虚假边缘，对噪声比较敏感。对比图3(a)和图3(c)，本发明的方法得到的分割结果不会出现漏分割现象；对比图3(b)和图3(c)，本发明的方法得到的分割结果不会出现过分割现象。It can be seen from Figure 3(a) that the GAC method is sensitive to initialization, and can only detect part of the contours of the SAR image with rich regions, but cannot obtain a complete segmentation result. It can be seen from Figure 3(b) that the existing CV method detects more false edges and is more sensitive to noise. Comparing Fig. 3 (a) and Fig. 3 (c), the segmentation result that the method of the present invention obtains can not occur omission segmentation phenomenon; Comparing Fig. 3 (b) and Fig. 3 (c), the segmentation result that the method of the present invention obtains does not Oversegmentation occurs.

仿真二：应用本发明和现有的侧地活动轮廓(GAC)方法和现有的CV方法分别对图4的测试图像进行分割实验，其分割结果如图5所示，其中：图5(a)是现有的GAC方法对图4的分割结果；图5(b)是现有的CV方法对图2的分割结果；图5(c)是本发明方法对图4的分割结果。Simulation two: apply the present invention and existing lateral ground active contour (GAC) method and existing CV method to carry out segmentation experiment respectively to the test image of Fig. 4, its segmentation result is as shown in Fig. 5, wherein: Fig. 5 (a ) is the segmentation result of the existing GAC method to Fig. 4; Fig. 5 (b) is the segmentation result of the existing CV method to Fig. 2; Fig. 5 (c) is the segmentation result of the inventive method to Fig. 4.

由图5(a)可看出，GAC方法对初始化比较敏感，对区域比较丰富的SAR图像只能检测出部分轮廓，而不能得到完整的分割结果。由图5(b)可看出，现有的CV方法检测出了比较多的虚假边缘，对噪声比较敏感。对比图5(a)和图5(c)，本发明的方法得到的分割结果不会出现漏分割现象；对比图5(b)和图5(c)，本发明的方法得到的分割结果不会出现过分割现象，且由于本发明结合了适合SAR图像相干斑噪声的Gamma统计模型，使本发明方法在不需要对相干斑噪声进行预处理的情况下，能够抑制SAR相干斑噪声的影响，并增强了真实边缘的定位性能，从而得到相对准确的分割结果。It can be seen from Fig. 5(a) that the GAC method is sensitive to initialization, and can only detect part of the contours of SAR images with rich regions, but cannot obtain complete segmentation results. It can be seen from Figure 5(b) that the existing CV method detects more false edges and is more sensitive to noise. Comparing Fig. 5 (a) and Fig. 5 (c), the segmentation result that the method of the present invention obtains can not occur omission segmentation phenomenon; Comparing Fig. 5 (b) and Fig. 5 (c), the segmentation result that the method of the present invention obtains does not appear There will be over-segmentation, and because the present invention combines the Gamma statistical model suitable for SAR image coherent speckle noise, the method of the present invention can suppress the influence of SAR coherent speckle noise without preprocessing the coherent speckle noise, And enhance the positioning performance of the real edge, so as to obtain relatively accurate segmentation results.

Claims

1. A level set SAR image segmentation method based on edge and regional information, comprising the following steps:

(1) Apply an edge detection operator with an exponentially weighted mean ratio to the SAR image I to be segmented to perform edge detection, and obtain the edge intensity modulus |R _max |;

(2) The level set function φ is initialized into a signed distance function form, and the SAR image is divided into two regions Ω ₁ and Ω ₂ according to the positive or negative of the level set function value;

(3) According to the two regions Ω ₁ and Ω ₂ , calculate their corresponding estimated probability densities p ₁ and p ₂ :

3a) Compute the intensity mean c _i of regions Ω ₁ and Ω ₂ :

where i=1,2, (x,y) are image coordinates;

3b) Calculate the estimated probability densities p ₁ and p ₂ for regions Ω ₁ and Ω ₂ from the intensity mean c _i :

{p p}_{i i} ((I I ((x x,, y the y)))) = = {\frac{{L L}^{L L}}{{c c}_{i i} Γ Γ ((L L))} {((\frac{I I ((x x,, y the y))}{{c c}_{i i}}))}^{L L - - 11} e e}^{- - LI LI ((x x,, y the y)) / / {c c}_{i i}}

Wherein L is the visual number of the SAR image, I (x, y) is the pixel value of the pixel point at the coordinates (x, y) in the image, Γ ( ) is a Gamma function, i=1,2;

(4) in conjunction with step (1)-step (3), construct the total segmentation energy function model E _SAR :

{E E.}_{SAR SAR} = = - - {Σ Σ}_{i i = = 11}^{22} {&Integral; &Integral; &Integral; &Integral;}_{{Ω Ω}_{i i}} {λ λ}_{i i} log log (({p p}_{i i})) + + μ μ {&Integral; &Integral; &Integral; &Integral;}_{Ω Ω} g g ((| | {R R}_{max max} | |)) | | &dtri; &dtri; H h ((φ φ)) | | + + v v {&Integral; &Integral; &Integral; &Integral;}_{Ω Ω} \frac{11}{22} {((| | &dtri; &dtri; φ φ | | - - 11))}^{22},,

Where Ω is the entire image area, ie Ω=Ω ₁ +Ω ₂ ,

λ _i log(p _i ) is the area energy term,

is the marginal energy term,

is the energy item of the correction item, λ _i is the weight of the area energy item, i=1,2, μ is the weight of the edge energy item, ν is the weight of the correction item, λ _i ＞0, μ＞0, ν＞ 0, ▽φ means to calculate the gradient for φ, H(φ) is the Heaviside function, ▽H(φ) means to calculate the gradient for H(φ), g(|R _max |) is defined in the edge strength modulus |R _max The indicator function on |, the expression is as follows:

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

where k is a positive proportionality constant;

(5) segment the SAR image I according to the segmentation energy function model constructed in step (4):

\frac{&PartialD; &PartialD; φ φ}{&PartialD; &PartialD; t t} = = δ δ ((φ φ)) ((- - {λ λ}_{11} log log (({p p}_{11})) + + {λ λ}_{22} log log (({p p}_{22}))))

+ + μδ μδ ((φ φ)) div div {((g g ((| | {R R}_{max max} | | \frac{&dtri; &dtri; φ φ}{| | &dtri; &dtri; φ φ | |})) + + v v ((Δφ Δφ - - div div ((\frac{&dtri; &dtri; φ φ}{| | &dtri; &dtri; φ φ | |}))))}^{,,}

Where δ(φ) is the Dirac function, Δ is the Laplacian operator;

5b) For the gradient descent flow equation

Discretization, the following expression is obtained:

\frac{{φ φ}^{n no + + 11} - - {φ φ}^{n no}}{Δt Δt} = = δ δ ((φ φ)) ((- - {λ λ}_{11} log log (({P P}_{11})) + + {λ λ}_{22} log log (({p p}_{22}))))

+ + μδ μδ ((φ φ)) div div {((g g ((| | {R R}_{max max} | | \frac{&dtri; &dtri; φ φ}{| | &dtri; &dtri; φ φ | |})) + + v v ((Δφ Δφ - - div div ((\frac{&dtri; &dtri; φ φ}{| | &dtri; &dtri; φ φ | |}))))}^{,,}

Where φ ⁿ⁺¹ represents the level set function after the n+1th iteration, φ ⁿ represents the level set function after the nth iteration, and Δt is the iteration step size;

and

5d) Judging whether the level set function converges and reaches the maximum number of iterations of 100 times, if not, go to step (3), use

and

Substitute Ω ₁ and Ω ₂ to continue iterating, otherwise stop iterating, and get

and That is the final segmentation result.

2. The level set SAR image segmentation method according to claim 1, wherein the edge detection operator of the SAR image I to be segmented using the exponential weighted mean ratio of the step (1) described in the step (1) carries out edge detection, as follows:

1a) Construct the function expressions of the causal filter f ₁ and the non-causal filter f ₂ respectively:

f ₁ (z) = ab ^z H(z), f ₂ (z) = ab ^-z H(-z),

Where z=1,2,...N,N>1, z is the independent variable of the function, and z=1,2,...N,N is a positive integer of N>1; Both a and b are constant, and satisfy 0<b<e ^-a <1, H( ) is a Heaviside function;

1b) Construct the functional expression of the exponential smoothing filter f according to the filters f ₁ and f ₂ :

f f ((z z)) = = \frac{11}{11 + + b b} {f f}_{11} ((z z)) + + \frac{b b}{11 + + b b} {f f}_{22} ((z z - - 11)),,

1c) According to the constructed filters f, f ₁ , f ₂ , calculate the exponentially weighted mean value μ _x1 of filter f ₁ in the horizontal direction, the exponentially weighted mean value μ _x2 of filter f ₂ in the horizontal direction, filter f ₁ The exponentially weighted mean value μ _y1 in the vertical direction, the exponentially weighted mean value μ _y2 of the filter f ₂ in the vertical direction, each exponentially weighted mean value μ _x1 , μ _x2 , μ _y1 , μ _y2 :

\{\begin{matrix} {μ μ}_{x x 11} = = {f f}_{11} ((x x)) * * ((f f ((y the y)) \cdot &Center Dot; I I ((x x,, y the y)))) \\ {μ μ}_{x x 22} = = {f f}_{22} ((x x)) * * ((f f ((y the y)) \cdot \cdot I I ((x x,, y the y)))) \\ {μ μ}_{y the y 11} = = {f f}_{11} ((y the y)) \cdot &Center Dot; ((f f ((x x)) * * I I ((x x,, y the y)))) \\ {μ μ}_{y the y 22} = = {f f}_{22} ((y the y)) \cdot \cdot ((f f ((x x)) * * I I ((x x,, y the y)))) \end{matrix},,,

Among them: x is the coordinate variable in the horizontal direction, y is the coordinate variable in the vertical direction, * represents the convolution in the horizontal direction, · represents the convolution in the vertical direction;

1d) Apply the results of 1c) to find the intensity modulus R _xmax (x,y) in the horizontal direction and the intensity modulus R _ymax (x,y) in the vertical direction:

\{\begin{matrix} {R R}_{x x max max} ((x x,, y the y)) = = max max {{\frac{{μ μ}_{x x 11} ((x x - - 11,, y the y))}{{μ μ}_{x x 22} ((x x + + 11,, y the y))},, \frac{{μ μ}_{x x 22} ((x x + + 11,, y the y))}{{μ μ}_{x x 11} ((x x - - 11,, y the y))}}} \\ {R R}_{y the y max max} ((x x,, y the y)) = = max max {{\frac{{μ μ}_{y the y 11} ((x x,, y the y - - 11))}{{μ μ}_{y the y 22} ((x x,, y the y + + 11))},, \frac{{μ μ}_{y the y 22} ((x x,, y the y + + 11))}{{μ μ}_{y the y 11} ((x x,, y the y - - 11))}}} \end{matrix},,

1e) According to R _xmax (x, y) and R _ymax (x, y), calculate the edge strength modulus |R _max |:

| | {R R}_{max max} ((x x,, y the y)) | | = = \sqrt{{R R}_{x x max max}^{22} ((x x,, y the y)) + + {R R}_{y the y max max}^{22} ((x x,, y the y))} . .