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

Abstract

The invention discloses a level set synthetic aperture radar (SAR) image segmentation method based on edge and regional information, and mainly solves the problems that an SAR image with a vague edge is difficult to segment and the real edge of the SAR image cannot be accurately positioned by using an existing level set method. The method has the steps that: first, the edge of the SAR image is detected by an index weighted average ratio operator, the strength modulus absolute value of Rmax of the edge is worked out, and the energy item of the edge is constructed; second, the function phi of the level set is initialized, the SAR image is segmented into two areas omega 1 and omega 2, and the strength means c1 and c2 of the two areas are calculated; the estimated probability density p1 and p2 of the areas omega 1 and omega 2 are calculated according to c1 and c2, and a regional energy item is constructed; finally, a correction energy item for preventing re-initialization is added, an aggregate energy function electronically steerable array radar (ESAR) is constructed, a gradient descending flow equation is worked out through a variational method, the level set phi is upgraded, and new segmented areas are obtained. An experiment result shows that the segmentation method realized according to the invention can achieve a more ideal segmentation effect, and can be used for the edge detection and the target identification of the SAR image.

Description

Level set SAR image partition method based on edge and area information

Technical field

The invention belongs to technical field of image processing, relate to the SAR image and cut apart, a kind of level set SAR image partition method that combines based on edge and area information specifically can be used for the cutting apart of SAR image, rim detection and target identification.

Background technology

Synthetic-aperture radar SAR is the active radar of a kind of high resolving power, but has advantages such as round-the-clock, round-the-clock, resolution height side-looking imaging, can be applicable to numerous areas such as military affairs, agricultural, navigation, geographical supervision.The SAR image is widely used in object detection field, and it then is the important step of handling graphical analysis from image that the SAR image is cut apart, and is target classification and base of recognition.Because SAR is a kind of coherence imaging system, the SAR image is to the electromagnetic scattering characteristic of target and the reaction of architectural characteristic in itself, and its imaging effect depends on radar parameter and region electromagnetic parameter to a great extent.Because the singularity of SAR imaging system, make information representation mode and the optical imagery of SAR image that very big-difference be arranged, and can be subjected to the influence of many geometric properties such as coherent speckle noise and shade, just because of these factors make the dividing method that is applicable to optical imagery no longer suitable to the SAR image.

Conventional on the whole SAR dividing method great majority are still cut apart the SAR image with traditional parted pattern, as some based on the cutting apart of threshold value, based on the edge cut apart and based on dividing method of region growing etc., but a lot of dividing methods wherein are subjected to the influence of coherent speckle noise in the SAR image easily, so these class methods all need to suppress coherent speckle noise by pre-service earlier.These preprocess methods have lost a lot of edges and target information inevitably when suppressing coherent speckle noise, influenced segmentation effect.

Level set image segmentation method is one of important images disposal route.The advantage of Level Set Method is to adapt to changes of topology structure, and algorithm stability is higher.Utilize Level Set Method research SAR image segmentation problem, by setting up suitable energy functional, can in the definition of energy functional, introduce image area information, do not needing to obtain segmentation result more accurately for the SAR image that influenced by coherent speckle noise under the pretreated situation of coherent spot.

Classical Level Set Method is divided based on marginal information with based on two kinds of methods of area information.Geodesic line movable contour model GAC model is based in the Level Set Method of marginal information most widely used a kind of, it can be under not additional any extraneous controlled condition, topologies change when freely handling curvilinear motion, but this model only utilizes the gradient information of image, because the edge in the image motion is not all to be desirable step edges, be difficult to be partitioned into the homogeneous region in the ill-defined image, in addition when the object in the image has darker depression border, the GAC model may make the evolution curve stop at a certain local minimum state of value, and is not consistent with the borderline phase of object.The CV model that Chan and Vese propose is a kind of level set dividing method based on the zone, and as a kind of zone level collection parted pattern that can effectively improve curve evolvement topology adaptive ability, it has utilized the global area information of image.But equally also there is weak point in it, the internal energy term of this model just guarantees zero level collection line smoothing, and do not consider the inwardness that level set function itself is intrinsic, also need in some applications level set function is reinitialized, so that it is near the symbolic distance function, guarantee the stability of numerical solution, in addition, only utilized the area information of image in the model, and the marginal information of another important images is not utilized well, cuts apart at some especially may produce the image border to the inhomogeneous image of gray scale in the application and locate inaccurate defective.

Though the level set dividing method has obtained great success in optical imagery and medical image segmentation, cut apart the field at the SAR image, research also fewer.Present two mechanisms leading of the people such as Chen of people such as the well-known Refregier that France is arranged and the U.S. in the world, and two mechanisms respectively lay particular stress on, the former lay particular stress on research geodesic line active contour GAC model the SAR image cut apart and rim detection aspect application, and this model depends on image gradient definition edge detection operator, because the gradient information itself in the GAC model just is not suitable for the SAR image, so produce many erroneous segmentation, the latter mainly studies geometric active contour CV model cutting apart multizone SAR image, since the CV model be one based on the parted pattern of additive Gaussian noise, if directly this model is applied to the SAR image, because the existence of the property taken advantage of coherent speckle noise, the testing result that obtains is very undesirable, and it has only utilized the area information of SAR image, and the very important marginal information of SAR image is not fully utilized, and can't correctly cut apart the SAR image.

Summary of the invention

The objective of the invention is to the deficiency at above-mentioned prior art, provide a kind of based on edge and the regional level set SAR image partition method that combines, accurately the cutting apart of ill-defined SAR image-region improved the positioning performance at edge, thereby improve the quality of cutting apart of SAR image realizing.

Realize that the object of the invention technical thought is: according to the statistical property of SAR image, using Gamma distributes as the probability density function of SAR image, and calculate the estimated value of being divided the probability density in two zones that obtain by level set, with the logarithm of each estimated value as area information, using simultaneously has the testing result of edge detection operator of the exponential weighting average ratio of fine edge positioning performance to replace Grad as marginal information to the SAR image, add at last and avoid the correction term that reinitializes, the modeling of energy function is carried out in this three partial informations combination, thereby reached desirable segmentation result.Its concrete implementation procedure comprises as follows:

(1) edge detection operator to SAR image I exponential weighted mean ratio to be split carries out rim detection, obtains edge strength mould value | R _Max|;

(2) level set function φ is initialized to the symbolic distance functional form, positive and negative according to the level set function value is divided into two regional Ω with the SAR image ₁And Ω ₂

(3) according to two regional Ω ₁And Ω ₂, calculate its corresponding estimated probability density p ₁And p ₂:

3a) zoning Ω ₁And Ω ₂Strength mean value c _i:

I=1 wherein, 2, (x y) is image coordinate;

3b) according to strength mean value c _i, zoning Ω ₁And Ω ₂Estimated probability density p ₁And p ₂:

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

Wherein L is the number of looking of SAR image, and Γ () is the Gamma function, i=1,2;

(4) integrating step (1)-step (3) is constructed total energy function model E of cutting apart _SAR:

$E_{\mathrm{SAR}}=-\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\∫\∫}_{{\mathrm{\Ω}}_i}{\mathrm{\λ}}_i\log \left(p_i\right)+\mathrm{\μ}{\∫\∫}_{\mathrm{\Ω}}g\left(\right|R_{\max }\left|\right)\left|\mathrm{\φH}\left(\mathrm{\φ}\right)\right|+v{\∫\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿\mathrm{\φ}|-1)}^2,$

Wherein Ω is entire image zone, i.e. Ω=Ω ₁+ Ω ₂,

Be the region energy item,

Be the edge energy item,

Be the correction term energy term, λ _iBe the weights of region energy item, i=1,2, μ are the weights of edge energy item, and v is the weights of correction term, λ _i＞0, μ＞0, v＞0, Be that expression is asked gradient to φ, H (φ) is the Heaviside function,

Expression is asked gradient to H (φ), g (| R _Max|) be to be defined in edge strength mould value | R _Max| on indicator function, expression formula is as follows:

$g\left(\right|R_{\max }\left|\right)=\frac{1}{1+{\left|R_{\max }\right|}^2/k^2},$ Wherein k is positive proportionality constant;

(5) according to the energy function model of cutting apart of step (4) structure the SAR image I is cut apart:

5a) total energy function model of cutting apart is used the variational method, obtain gradient katabatic drainage equation

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\δ}\left(\mathrm{\φ}\right)(-{\mathrm{\λ}}_1\log \left(p_1\right)+{\mathrm{\λ}}_2\log \left(p_2\right))$

$+\mathrm{\μ\δ}\left(\mathrm{\φ}\right)\mathrm{div}(g\left(\right|R_{\max }\left|\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+v(\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)),$

Wherein δ (φ) is the Dirac function, and Δ is Laplace operator;

5b) to gradient katabatic drainage equation

Discretize obtains following expression:

$\frac{{\mathrm{\φ}}^{n+1}-{\mathrm{\φ}}^n}{\mathrm{\Δt}}=\mathrm{\δ}\left(\mathrm{\φ}\right)(-{\mathrm{\λ}}_1\log \left(p_1\right)+{\mathrm{\λ}}_2\log \left(p_2\right))$

$+\mathrm{\μ\δ}\left(\mathrm{\φ}\right)\mathrm{div}(g\left(\right|R_{\max }\left|\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+v(\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)),$

φ wherein ^N+1Represent the level set function after the iteration the n+1 time, φ ⁿRepresent the level set function after the iteration the n time, Δ t is iteration step length;

5c) according to step 5b) try to achieve new level set function φ ^N+1, by φ ^N+1Positive and negative value obtain new cut zone

With

5d) whether the determined level set function is restrained and is reached maximum iterations 100 times, does not then forward step (3) to if do not satisfy, and uses

With

Substitute Ω ₁And Ω ₂Continue iteration, otherwise stop iteration, obtain With It namely is final segmentation result.

The present invention has the following advantages compared with prior art:

1, the present invention is incorporated into the level set model with edge and area information, the area information that had both comprised the SAR image, comprised important marginal information again, not only the image border is had and well cut apart the localization effect, improve the accuracy to object edge location again, also improved segmentation effect and robustness to the inhomogeneous SAR image of edge fog and gray scale simultaneously.

2, the present invention combines the Gamma statistical model that is fit to SAR image coherent speckle noise, directly constructed the area item energy model of level set energy function, parted pattern is carried out under the pretreated situation coherent speckle noise not needing, can obtain relatively accurate segmentation result to the SAR image that influenced by coherent speckle noise.

3, the present invention has used the gradient information that the SAR image is had the edge strength mould value alternative image that the edge detection operator of the exponential weighting average ratio of good edge positioning performance obtains, and cuts apart thereby make algorithm better be used in the SAR image.

Through with the emulation of existing level set dividing method contrast, verified that the present invention is cut apart the SAR image can obtain better segmentation effect.

Description of drawings

Fig. 1 is process flow diagram of the present invention;

Fig. 2 is first width of cloth test pattern;

Fig. 3 is with the present invention and the existing method segmentation effect comparison diagram to Fig. 2;

Fig. 4 is second width of cloth test pattern;

Fig. 5 is with the present invention and the existing method segmentation effect comparison diagram to Fig. 4;

Embodiment

With reference to Fig. 1, specific implementation step of the present invention is as follows:

Step 1 is carried out rim detection with the edge detection operator of SAR image I exponential weighted mean ratio to be split, obtains edge strength mould value | R _Max|.

The detailed process that realizes this step is as follows:

1a) construct causal filter f respectively ₁With non-causal filter f ₂Function expression:

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

Wherein z is argument of function, and z=1,2 ... N, N are the positive integer of N＞1; A and b are constant, and satisfy 0＜b＜e ^-a＜1, H () is the Heaviside function;

1b) according to wave filter f ₁And f ₂The function expression of structural index smoothing filter f:

$f\left(z\right)=\frac{1}{1+b}{f}_1\left(z\right)+\frac{b}{1+b}{f}_2(z-1),$

1c) according to the wave filter f that constructs, f ₁, f ₂, calculating filter f ₁Exponential weighting average μ in the horizontal direction _X1, wave filter f ₂Exponential weighting average μ in the horizontal direction _X2, wave filter f ₁Exponential weighting average μ in vertical direction _Y1, wave filter f ₂Exponential weighting average μ in vertical direction _Y2, each exponential weighting average μ _X1, μ _X2, μ _Y1, μ _Y2:

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

Wherein: x is the coordinate variable of horizontal direction, and y is the coordinate variable of vertical direction, and * represents the convolution of horizontal direction, represents the convolution of vertical direction;

1d) use 1c) the result, ask the intensity mould value R of horizontal direction _Xmax(x, y) and the intensity mould value R of vertical direction _Ymax(x, y):

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

1e) according to R _Xmax(x, y) and R _Ymax(x, y) ask edge strength mould value | R _Max|:

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

Step 2 is initialized to the symbolic distance functional form with level set function φ, and positive and negative according to the level set function value is divided into two regional Ω with the SAR image ₁And Ω ₂

The detailed process that realizes this step is expressed as follows: do a rectangle with PC Tools in image I to be split, this rectangle equation is as initialization level set function φ, when φ＞0, and expression rectangle perimeter Ω ₁, when φ＜0, expression rectangle interior zone Ω ₂, therefore, positive and negative according to the level set function value is divided into two regional Ω with the SAR image ₁And Ω ₂

Step 3 is according to the existing formula of asking SAR image intensity average, zoning Ω ₁Strength mean value c ₁With regional Ω ₂Strength mean value c ₂, expression formula is:

$\left\{\begin{array}{c}{c}_1=\frac{{\∫\∫}_{{\mathrm{\Ω}}_1}I(x,y)\mathrm{dxdy}}{{\∫\∫}_{{\mathrm{\Ω}}_1}\mathrm{dxdy}}\\ c_2=\frac{{\∫\∫}_{{\mathrm{\Ω}}_2}I(x,y)\mathrm{dxdy}}{{\∫\∫}_{{\mathrm{\Ω}}_2}\mathrm{dxdy}}\end{array}\right.,$ Wherein (x y) is image coordinate,

Step 4, the strength mean value c that tries to achieve according to step 3 _iWith the formula of existing SAR image probability density function, zoning Ω ₁And Ω ₂Estimated probability density p ₁And p ₂:

$\left\{\begin{array}{c}{p}_1\left(I\right)=\frac{L^L}{c_1\mathrm{\Γ}\left(L\right)}{\left(\frac{I}{c_1}\right)}^{L-1}{e}^{-\mathrm{LI}/c_1}\\ p_2\left(I\right)=\frac{L^L}{c_2\mathrm{\Γ}\left(L\right)}{\left(\frac{I}{c_2}\right)}^{L-1}{e}^{-\mathrm{LI}/c_2}\end{array}\right.,$

Wherein L is the number of looking of SAR image, and I is SAR image to be split, and Γ () is the Gamma function.

Step 5 is in conjunction with four the total energy function E of step structure in front _SAR

The detailed process that realizes this step is as follows:

5a) edge detection operator according to step 1 exponential weighted mean ratio carries out rim detection, obtains edge strength mould value | R _Max|, calculate in edge strength mould value | R _Max| on edge indicator function g (| R _Max|):

$g\left(\right|R_{\max }\left|\right)=\frac{1}{1+{\left|R_{\max }\right|}^2/k^2},$ Wherein k is positive proportionality constant;

5b) the edge indicator function that calculates according to step (5a), structure edge energy item, its expression formula is:

$\mathrm{\μ}{\∫\∫}_{\mathrm{\Ω}}g\left(\right|R_{\max }\left|\right)|\▿H\left(\mathrm{\φ}\right)|,$

5c) calculate the regional Ω that tries to achieve according to step 4 ₁And Ω ₂Estimated probability density p ₁And p ₂, the structure realm energy term, its expression formula is: $-{\∫\∫}_{{\mathrm{\Ω}}_1}{\mathrm{\λ}}_1\log \left(p_1\right)-{\∫\∫}_{{\mathrm{\Ω}}_2}{\mathrm{\λ}}_2\log \left(p_2\right),$

5d) integrating step (5b) and the edge energy item of (5c) trying to achieve and region energy item add and avoid the punishment energy term that reinitializes, construct total energy function E _SAR, its expression formula is as follows:

$E_{\mathrm{SAR}}=-{\∫\∫}_{{\mathrm{\Ω}}_1}{\mathrm{\λ}}_1\log \left(p_1\right)-{\∫\∫}_{{\mathrm{\Ω}}_2}{\mathrm{\λ}}_2\log \left(p_2\right)+\mathrm{\μ}{\∫\∫}_{\mathrm{\Ω}}g\left(\right|R_{\max }\left|\right)|\▿H\left(\mathrm{\φ}\right)|+v{\∫\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿\mathrm{\φ}|-1)}^2,$

Wherein Ω is entire image zone, i.e. Ω=Ω ₁+ Ω ₂, λ ₁And λ ₂Be respectively interior zone Ω ₁Energy term and perimeter Ω ₂The weights of energy term, μ are the weights of edge energy item, and v is the weights of correction term, λ ₁＞0, λ ₂＞0, μ＞0, v＞0, Be that expression is asked gradient to φ, H (φ) is the Heaviside function,

Expression is asked gradient to H (φ).

Step 6 is to energy function E _SARUse the variational method, obtain gradient katabatic drainage equation And iteration renewal level set function φ, obtain final cut zone

With

The detailed process that realizes this step is as follows:

6a) total energy function model of cutting apart is used the variational method, obtain gradient katabatic drainage equation

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\δ}\left(\mathrm{\φ}\right)(-{\mathrm{\λ}}_1\log \left(p_1\right)+{\mathrm{\λ}}_2\log \left(p_2\right))$

$+\mathrm{\μ\δ}\left(\mathrm{\φ}\right)\mathrm{div}(g\left(\right|R_{\max }\left|\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+v(\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)),$

Wherein δ (φ) is the Dirac function, and Δ is Laplace operator;

6b) to gradient katabatic drainage equation

Discretize obtains following expression:

$\frac{{\mathrm{\φ}}^{n+1}-{\mathrm{\φ}}^n}{\mathrm{\Δt}}=\mathrm{\δ}\left(\mathrm{\φ}\right)(-{\mathrm{\λ}}_1\log \left(p_1\right)+{\mathrm{\λ}}_2\log \left(p_2\right))$

$+\mathrm{\μ\δ}\left(\mathrm{\φ}\right)\mathrm{div}(g\left(\right|R_{\max }\left|\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+v(\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)),$

φ wherein ^N+1Represent the level set function after the iteration the n+1 time, φ ⁿRepresent the level set function after the iteration the n time, Δ t is iteration step length;

6c) according to step 6b) try to achieve new level set function φ ^N+1, by φ ^N+1Positive and negative value obtain new cut zone With

Step 7 judges whether convergence.

The fundamental purpose of this step is the iteration situation according to level set, and whether the determined level set function is restrained and reached maximum iterations, if gradient katabatic drainage equation does not reach steady state (SS) and do not reach maximum iterations, then forwards step (3) to, uses

With

Substitute Ω ₁And Ω ₂Continue iteration, otherwise stop iteration, obtain final zone With It namely is final segmentation result.

Effect of the present invention can further specify by following emulation:

1, simulated conditions

Hardware platform is: Intel Core2 Duo CPU E6550@2.33GHZ, 2GB RAM.

Software platform is: MATLAB 7.9.

2, emulation content and result

Emulation one: use the present invention and existing side ground active contour (GAC) method and existing C V method and respectively the test pattern of Fig. 2 carried out split-run test, its segmentation result as shown in Figure 3, wherein: Fig. 3 (a) is that existing GAC method is to the segmentation result of Fig. 2; Fig. 3 (b) is that existing C V method is to the segmentation result of Fig. 2; Fig. 3 (c) is that the inventive method is to the segmentation result of Fig. 2.

Can be found out that by Fig. 3 (a) the GAC method is relatively more responsive to initialization, the SAR image relatively abundanter to the zone can only detect partial contour, and can not obtain complete segmentation result.Can find out that by Fig. 3 (b) existing C V method has detected many false edges, responsive to noise ratio.Leaking can not appear in comparison diagram 3 (a) and Fig. 3 (c), the segmentation result that method of the present invention obtains cuts apart phenomenon; Comparison diagram 3 (b) and Fig. 3 (c), the segmentation result that method of the present invention obtains the over-segmentation phenomenon can not occur.

Emulation two: use the present invention and existing side ground active contour (GAC) method and existing C V method and respectively the test pattern of Fig. 4 carried out split-run test, its segmentation result as shown in Figure 5, wherein: Fig. 5 (a) is that existing GAC method is to the segmentation result of Fig. 4; Fig. 5 (b) is that existing C V method is to the segmentation result of Fig. 2; Fig. 5 (c) is that the inventive method is to the segmentation result of Fig. 4.

Can be found out that by Fig. 5 (a) the GAC method is relatively more responsive to initialization, the SAR image relatively abundanter to the zone can only detect partial contour, and can not obtain complete segmentation result.Can find out that by Fig. 5 (b) existing C V method has detected many false edges, responsive to noise ratio.Leaking can not appear in comparison diagram 5 (a) and Fig. 5 (c), the segmentation result that method of the present invention obtains cuts apart phenomenon; Comparison diagram 5 (b) and Fig. 5 (c), the segmentation result that method of the present invention obtains the over-segmentation phenomenon can not occur, and because the present invention combines the Gamma statistical model that is fit to SAR image coherent speckle noise, the inventive method is carried out under the pretreated situation coherent speckle noise not needing, can suppress the influence of SAR coherent speckle noise, and strengthened the positioning performance of true edge, thereby obtain relatively accurate segmentation result.

Claims (2)

1. level set SAR image partition method based on edge and area information may further comprise the steps:

(1) edge detection operator to SAR image I exponential weighted mean ratio to be split carries out rim detection, obtains edge strength mould value | R _Max|;

(2) level set function φ is initialized to the symbolic distance functional form, positive and negative according to the level set function value is divided into two regional Ω with the SAR image ₁And Ω ₂

(3) according to two regional Ω ₁And Ω ₂, calculate its corresponding estimated probability density p ₁And p ₂:

3a) zoning Ω ₁And Ω ₂Strength mean value c _i:

I=1 wherein, 2, (x y) is image coordinate;

3b) according to strength mean value c _i, zoning Ω ₁And Ω ₂Estimated probability density p ₁And p ₂:

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

Wherein L is the number of looking of SAR image, I (x, y) be in the image coordinate for (x y) locates the pixel value of pixel, and Γ () is the Gamma function, i=1,2;

(4) integrating step (1)-step (3) is constructed total energy function model E of cutting apart _SAR:

$E_{\mathrm{SAR}}=-\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\∫\∫}_{{\mathrm{\Ω}}_i}{\mathrm{\λ}}_i\log \left(p_i\right)+\mathrm{\μ}{\∫\∫}_{\mathrm{\Ω}}g\left(\right|R_{\max }\left|\right)|\▿H\left(\mathrm{\φ}\right)|+v{\∫\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿\mathrm{\φ}|-1)}^2,$

Wherein Ω is entire image zone, i.e. Ω=Ω ₁+ Ω ₂,

λ _iLog(p _i) be the region energy item,

Be the edge energy item,

Be the correction term energy term, λ _iBe the weights of region energy item, i=1,2, μ are the weights of edge energy item, and ν is the weights of correction term, λ _i＞0, μ＞0, ν＞0, ▽ φ is that expression is asked gradient to φ, and H (φ) is the Heaviside function, and ▽ H (φ) expression is asked gradient to H (φ), g (| R _Max|) be to be defined in edge strength mould value | R _Max| on indicator function, expression formula is as follows:

$g\left(\right|R_{\max }\left|\right)=\frac{1}{1+{\left|R_{\max }\right|}^2/k^2},$ Wherein k is positive proportionality constant;

(5) according to the energy function model of cutting apart of step (4) structure the SAR image I is cut apart:

5a) total energy function model of cutting apart is used the variational method, obtain gradient katabatic drainage equation

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\δ}\left(\mathrm{\φ}\right)(-{\mathrm{\λ}}_1\log \left(p_1\right)+{\mathrm{\λ}}_2\log \left(p_2\right))$

$+\mathrm{\μ\δ}\left(\mathrm{\φ}\right)\mathrm{div}{(g\left(\right|R_{\max }\left|\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+v(\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right))}^{,}$

Wherein δ (φ) is the Dirac function, and Δ is Laplace operator;

5b) to gradient katabatic drainage equation

Discretize obtains following expression:

$\frac{{\mathrm{\φ}}^{n+1}-{\mathrm{\φ}}^n}{\mathrm{\Δt}}=\mathrm{\δ}\left(\mathrm{\φ}\right)(-{\mathrm{\λ}}_1\log \left(P_1\right)+{\mathrm{\λ}}_2\log \left(p_2\right))$

$+\mathrm{\μ\δ}\left(\mathrm{\φ}\right)\mathrm{div}{(g\left(\right|R_{\max }\left|\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+v(\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right))}^{,}$

φ wherein ^N+1Represent the level set function after the iteration the n+1 time, φ ⁿRepresent the level set function after the iteration the n time, Δ t is iteration step length;

5c) according to step 5b) try to achieve new level set function φ ^N+1, by φ ^N+1Positive and negative value obtain new cut zone

With

5d) whether the determined level set function is restrained and is reached maximum iterations 100 times, does not then forward step (3) to if do not satisfy, and uses

With

Substitute Ω ₁And Ω ₂Continue iteration, otherwise stop iteration, obtain

With It namely is final segmentation result.

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

1a) construct causal filter f respectively ₁With non-causal filter f ₂Function expression:

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

Z=1 wherein, 2 ... N, N＞1, z is argument of function, and z=1,2 ... N, N are the positive integer of N＞1; A and b are constant, and satisfy 0＜b＜e ^-a＜1, H () is the Heaviside function;

1b) according to wave filter f ₁And f ₂The function expression of structural index smoothing filter f:

$f\left(z\right)=\frac{1}{1+b}{f}_1\left(z\right)+\frac{b}{1+b}{f}_2(z-1),$

1c) according to the wave filter f that constructs, f ₁, f ₂, calculating filter f ₁Exponential weighting average μ in the horizontal direction _X1, wave filter f ₂Exponential weighting average μ in the horizontal direction _X2, wave filter f ₁Exponential weighting average μ in vertical direction _Y1, wave filter f ₂Exponential weighting average μ in vertical direction _Y2, each exponential weighting average μ _X1, μ _X2, μ _Y1, μ _Y2:

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

Wherein: x is the coordinate variable of horizontal direction, and y is the coordinate variable of vertical direction, and * represents the convolution of horizontal direction, represents the convolution of vertical direction;

1d) use 1c) the result, ask the intensity mould value R of horizontal direction _Xmax(x, y) and the intensity mould value R of vertical direction _Ymax(x, y):

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

1e) according to R _Xmax(x, y) and R _Ymax(x, y) ask edge strength mould value | R _Max|:

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

Images