CN103839240A - SAR image despeckle method based on probability model and NSDT - Google Patents

Abstract

The invention discloses an SAR image despeckle method based on a probability model and NSDT. The SAR image despeckle method solves the problems that in the prior art, a margin and detail information of an SAR image after despeckle are not insufficiently maintained, and pseudo lines appear in a homogeneous area of the SAR image after despeckle. The SAR image despeckle method specifically comprises the following steps that (1) an image is input, (2) non-subsample direction wave transformation is conducted, (3) tag masks of high-frequency subband coefficients are obtained, (4) a likelihood ratio is obtained, (5) a priori ratio is obtained, (6) the high-frequency subband coefficients after despeckle are obtained, (7) non-subsample direction wave inverse transformation is conducted, and (8) the image is output. According to the SAR image despeckle method based on the probability model and the NSDT, statistical properties of an SAR image in a non-subsample direction wave transformation domain can be effectively captured, a margin and detail information of the SAR image after despeckle can be clearly maintained, and speckle noise in a homogeneous area of the SAR image to be subjected to despeckle can be effectively filtered out.

Description

SAR image speckle method based on probability model and NSDT

Technical field

The invention belongs to technical field of image processing, further relate to synthetic-aperture radar (Synthetic Aperture Radar, SAR) one of image speckle technical field is based on probability model and non-lower sampling direction wave conversion (Nonsubsampled Directionlet, NSDT) SAR image speckle method, the present invention can be used for removing the speckle noise in SAR image.

Background technology

Synthetic-aperture radar (Synthetic Aperture Radar, SAR) is the relevant imaging radar of a kind of high-resolution microwave, has imaging characteristics round-the-clock, round-the-clock, contains a large amount of signal characteristics.Along with the development of SAR image processing techniques, be widely used in military field and civil area.But the coherent imaging mechanism of SAR image causes image to be generally subject to the pollution of speckle noise, and SAR image post-processed is brought to puzzlement, therefore, SAR image speckle method is a key technical problem of SAR image interpretation and analysis.

When the key of SAR image speckle is filtering speckle noise, keep as much as possible edge and the detailed information of image.Along with multi-scale geometric analysis (Multiscale Geometric Analysis, MGA) development, the method of frequency domain filtering has been widely applied to going in the middle of spot of SAR image, conventional frequency domain filtering method has small echo (Wavelet) to remove spot, non-down sampling contourlet (Nonsubsampled Contourlet, NSCT) go spot, non-lower sampling direction wave (Nonsubsampled Directionlet, NSDT) to go the methods such as spot.

Xian Electronics Science and Technology University applies for a patent in " the SAR image method for reducing speckle based on neighborhood directivity information " (number of patent application: 200910022866.1, publication number: CN101566688A) and has proposed a kind of SAR image speckle method based on neighborhood directivity information.This patented claim utilizes non-down sampling contourlet transform (Nonsubsampled Contourlet, NSCT) to picture breakdown, obtain conditional likelihood ratio according to logarithm Gaussian distribution and blended index distribution, obtain priori ratio according to neighborhood directivity model, finally reduce than obtaining the high-frequency sub-band that reduction factor pair decomposites according to likelihood ratio and priori, the coefficient after reduction is changed is reconstructed and obtains spot image.Although the spot image that goes that this patented claim obtains can keep edge and the grain details of SAR image, but the deficiency still existing is, this patented claim goes to homogeneous region in the SAR image after spot to still have pseudo-line, and goes still can lose in spot image some marginal informations and grain details.

Xian Electronics Science and Technology University apply for a patent in " based on non-lower sampling direction wave conversion and the SAR image de-noising method merging " (number of patent application: 201310092935.2, publication number: CN103177428A), proposed a kind of based on the conversion of non-lower sampling direction wave and the SAR image de-noising method merging.SAR image is divided into homogeneous region and target area by this patented claim, image is carried out to the conversion of non-lower sampling direction wave and obtain three prescriptions to wave system number, then three prescriptions are distinguished to denoising to wave system number, coefficient after denoising is reconstructed and obtains three corresponding width denoising images, finally three width denoising images are merged to the SAR image that obtains final denoising.Though this patented claim has advantages of the SAR image target area detailed information after denoising and keeps better, homogeneous region more level and smooth, but the deficiency still existing is, this patented claim does not take into full account the geometric properties of SAR image and the SAR image statistical property at non-lower sampling direction wave transform domain, so go details and marginal information in the SAR image after spot to have partial loss.

Summary of the invention

The object of the invention is to overcome above-mentioned the deficiencies in the prior art, proposed a kind of SAR image speckle method based on probability model and NSDT.The present invention can effectively remove the speckle noise of homogeneous area pixel point in SAR image, has kept edge and the detailed information of SAR image.

Comprise as follows for realizing the concrete steps of the object of the invention:

(1) the SAR image of the width spot to be gone of input option.

(2) non-lower sampling direction wave conversion:

The SAR image of the spot to be gone to input carries out the non-lower sampling direction wave conversion that yardstick is L, obtains high-frequency sub-band coefficient and low frequency sub-band coefficient under different scale.

(3) high-frequency sub-band coefficient is asked to label mask:

According to the following formula, obtain the label mask of high-frequency sub-band coefficient under different scale:

$X=\left\{\begin{array}{cc}0,& |D\&CenterDot;y|<{\mathrm{\&sigma;}}^2\\ 1,& |D\&CenterDot;y|\&GreaterEqual;{\mathrm{\&sigma;}}^2\end{array}\right.$

Wherein, X represents the label mask of high-frequency sub-band coefficient under different scale, 0 represents the label mask of spot SAR noise in image to be gone, 1 represents the label mask at edge and details in spot SAR image to be gone, D represents that the SAR image for the treatment of spot carries out the high-frequency sub-band coefficient under the rear yardstick l of non-lower sampling direction wave conversion, the SAR image that y represents to treat spot carries out the estimated value of the high-frequency sub-band coefficient under yardstick l+1 after the conversion of non-lower sampling direction wave, l ∈ L, σ represents the standard deviation of high-frequency sub-band coefficient under different scale.

(4) ask likelihood ratio:

According to the following formula, ask the likelihood ratio of high-frequency sub-band coefficient amplitude information under different scale according to probability distribution:

ξ＝p(r|1)/p(r|0)

Wherein, ξ represents the likelihood ratio of high-frequency sub-band coefficient amplitude information under different scale, p (r|1) represents the probability density function of the high-frequency sub-band coefficient amplitude information that under different scale, label mask is 1, p (r|0) represents the probability density function of the high-frequency sub-band coefficient amplitude information that under different scale, label mask is 0, and r represents amplitude information corresponding to high-frequency sub-band coefficient under different scale.

(5) ask priori ratio:

From neighbourhood model storehouse, choose 8 different neighbourhood models, according to the following formula, ask the priori ratio of high-frequency sub-band coefficient under different scale:

$\mathrm{\&eta;}=\left\{\begin{array}{cc}\max \left[{\mathrm{\&eta;}}_{\mathrm{\&alpha;}}\right],& M\&GreaterEqual;T,\mathrm{\&alpha;}\&Element;\{\mathrm{0,1},...7\}\\ {\mathrm{\&eta;}}_8,& \mathrm{others}\end{array}\right.$

Wherein, η represents the priori ratio of high-frequency sub-band coefficient under different scale, and max represents to get maxima operation, η _αrepresent the priori ratio of lower 8 the different neighbourhood model medium-high frequency sub-band coefficients of different scale, α represents the label of 8 different neighbourhood models, M represents the number of the high-frequency sub-band coefficient that under different scale, label mask is 1, T represents threshold value, its value is 2/3 of selected neighbourhood model medium-high frequency sub-band coefficients, α ∈ { 0,1, ... 7} represents that 8 different neighbourhood model labels are from 0 to 7, η ₈the priori ratio of 3 × 3 square neighbourhood model medium-high frequency sub-band coefficients under expression different scale.

(6) ask the high-frequency sub-band coefficient after spot:

(6a) ask according to the following formula, the collapse threshold of high-frequency sub-band coefficient:

ρ＝ηξ/(1+ηξ)

Wherein, ρ represents the collapse threshold of high-frequency sub-band coefficient, and η represents the priori ratio of high-frequency sub-band coefficient, and ξ represents the likelihood ratio of high-frequency sub-band coefficient.

(6b) be multiplied by high-frequency sub-band coefficient by the collapse threshold of high-frequency sub-band coefficient, obtain the high-frequency sub-band coefficient going after spot.

(7) non-lower sampling direction wave inverse transformation:

To low frequency sub-band coefficient and remove the high-frequency sub-band coefficient after spot, carrying out yardstick is the non-lower sampling direction wave inverse transformation of L, obtains the SAR image after spot.

(8) the SAR image after spot is removed in output.

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

First, because the present invention carries out the conversion of non-lower sampling direction wave to the spot SAR image to be gone of input, the conversion of non-lower sampling direction wave can effectively catch the detailed information such as profile and edge in SAR image, overcome prior art edge and detailed information and kept inadequate problem, made the present invention there is the speckle noise filtering of homogeneous region fully and edge keeps advantage clearly.

Second, the SAR figure that uses collapse threshold processing to treat spot due to the present invention carries out the high frequency coefficient after the conversion of non-lower sampling direction wave, collapse threshold is to try to achieve than with likelihood ratio according to priori, overcome prior art and can not effectively catch the problem of SAR image in the statistical property of non-lower sampling direction wave transform domain, made the present invention there is the homogeneous region of abundant filtering speckle noise.

Brief description of the drawings

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

Fig. 2 is that the SAR image field figure that the present invention and prior art are treated spot goes spot effect contrast figure;

Fig. 3 is that the SAR image airport figure that the present invention and prior art are treated spot goes spot effect contrast figure.

Embodiment

Below in conjunction with accompanying drawing, the present invention will be further described.

With reference to accompanying drawing 1, concrete steps of the present invention are as follows:

Step 1: input picture.

The SAR image of the one width spot to be gone of input option.The SAR image of the spot to be gone using in the embodiment of the present invention is as shown in Fig. 2 (a) and Fig. 3 (a), Fig. 2 (a) is the SAR image of spot the to be gone field figure that chooses from SAR image library, its size is 256 × 256, Fig. 3 (a) is the SAR image of spot the to be gone airport figure that chooses from SAR image library, and its size is 256 × 256.

Step 2: non-lower sampling direction wave conversion.

The SAR image of the spot to be gone to input carries out the non-lower sampling direction wave conversion that yardstick is L, obtains high-frequency sub-band coefficient and low frequency sub-band coefficient under different scale.

The concrete steps of non-lower sampling direction wave conversion are as follows:

According to the following formula, tectonic transition matrix:

$M_{\mathrm{\&Lambda;}}=\left[\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right]$

Wherein, M _Λrepresent transformation matrix, Λ represents one group of point that the linear combination of the integer vectors of two linear independences forms; a ₁, a ₂, b ₁, b ₂represent respectively transform matrix M _Λfour integer element values.

Choose 45 degree, 90 degree, 135 degree change direction, tectonic transition matrix, the transformation matrix of acquisition is respectively:

$\begin{array}{ccc}\left[\begin{array}{cc}1& 0\\ -1& 1\end{array}\right]& \left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]& \left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\end{array}$

Be multiplied by transformation matrix with the SAR image of spot to be gone, obtain three cosets.

Each coset along continuous straight runs is carried out to one-dimensional wavelet transform twice, vertically carry out one-dimensional wavelet transform one time, obtain seven high-frequency sub-band coefficients and a low frequency sub-band coefficient.Therefore, the SAR image wait removing spot of input is carried out can obtaining 7 × L high-frequency sub-band coefficient and L low frequency sub-band coefficient after non-lower sampling direction wave conversion that yardstick is L.

Step 3: high-frequency sub-band coefficient is asked to label mask.

According to the following formula, ask the label mask of high-frequency sub-band coefficient under different scale:

$X=\left\{\begin{array}{cc}0,& |D\&CenterDot;y|<{\mathrm{\&sigma;}}^2\\ 1,& |D\&CenterDot;y|\&GreaterEqual;{\mathrm{\&sigma;}}^2\end{array}\right.$

Wherein, X represents the label mask of high-frequency sub-band coefficient under different scale, 0 represents the label mask of spot SAR noise in image to be gone, 1 represents the label mask at edge and details in spot SAR image to be gone, D represents that the SAR image for the treatment of spot carries out the high-frequency sub-band coefficient under the rear yardstick l of non-lower sampling direction wave conversion, the SAR image that y represents to treat spot carries out the estimated value of the high-frequency sub-band coefficient under yardstick l+1 after the conversion of non-lower sampling direction wave, l ∈ L, σ represents the standard deviation of high-frequency sub-band coefficient under different scale.

According to the following formula, calculate the standard deviation of high-frequency sub-band coefficient under different scale:

$\mathrm{\&sigma;}=\frac{\mathrm{median}\left(\right|D-\mathrm{median}\left(D\right)\left|\right)}{0.6745}$

Wherein, σ represents the standard deviation of high-frequency sub-band coefficient under different scale, the median operation that represents to average, || represent to take absolute value operation, the SAR image that D represents to treat spot carries out the non-lower sampling direction wave conversion high-frequency sub-band coefficient under different scale afterwards.

Step 4: ask likelihood ratio.

According to the following formula, ask the likelihood ratio of high-frequency sub-band coefficient amplitude information under different scale according to probability distribution:

ξ＝p(r|1)/p(r|0)

Wherein, ξ represents the likelihood ratio of high-frequency sub-band coefficient amplitude information under different scale, p (r|1) represents the probability density function of the high-frequency sub-band coefficient amplitude information that under different scale, label mask is 1, p (r|0) represents the probability density function of the high-frequency sub-band coefficient amplitude information that under different scale, label mask is 0, and r represents amplitude information corresponding to high-frequency sub-band coefficient under different scale.

Label mask is that the probability density function that mixes that 1 high-frequency sub-band coefficient amplitude information is made up of with Gamma probability density function cuclear density function represents:

p(r|1)＝(1-τ)p ₁(r)+τp ₂(r|1)

Wherein, p (r|1) represents the probability density function of the high-frequency sub-band coefficient amplitude information that label mask is 1, and r represents the amplitude information that high-frequency sub-band coefficient is corresponding, and τ represents the mutual restricting relation between two, and its value is τ ∈ [0,1], p ₁(r) the cuclear density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 1, p ₂(r|1) the Gamma probability density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 1.

According to the following formula, asking label mask is the cuclear density function of 1 high-frequency sub-band coefficient amplitude information:

$p_1\left(r\right)=\frac{1}{\mathrm{N\&sigma;}\sqrt{2\mathrm{\&pi;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\exp [-{(r_s-r_j)}^2/{2\mathrm{\&sigma;}}^2]\&times;B[(u_s,v_s),(u_j,v_j)]$

Wherein, p ₁(r) the cuclear density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 1, be made up of gaussian kernel function and box kernel function, r represents the amplitude information that high-frequency sub-band coefficient is corresponding, and N represents the number of high-frequency sub-band coefficient, σ represents the standard deviation of high-frequency sub-band coefficient

represent from 1 to N summation, (r _s-r _j) represent the difference of two high-frequency sub-band coefficient amplitude information, s, j represents the position of high-frequency sub-band coefficient, B[(u _s, v _s), (u _j, v _j)] represent box kernel function, can calculate according to the following formula box kernel function.

$B[(u_s,v_s),(u_j,v_j)]=\left\{\begin{array}{cc}1& |u_s-u_j|\&le;L,|v_s-v_j|\&le;L\\ 0& \mathrm{otherwise}\end{array}\right.$

Wherein, (u _s, v _s) and (u _j, v _j) be illustrated respectively in the position s of high-frequency sub-band coefficient and the coordinate at j place, | u _s-u _j| represent the difference of horizontal ordinate, | v _s-v _j| represent the difference of ordinate.

According to the following formula, asking label mask is the Gamma probability density function of 1 high-frequency sub-band coefficient amplitude information:

$p_2\left(r\right|1)=\frac{{(r/b)}^2}{2b}{e}^{-r/b}$

Wherein, p ₂(r|1) the Gamma probability density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 1, r represents the amplitude information that high-frequency sub-band coefficient is corresponding, b represents parameter.

Label mask is that the probability density function that mixes that 0 high-frequency sub-band coefficient amplitude information is made up of with exponential probability density function cuclear density function represents:

p(r|0)＝(1-τ)p ₀(r)+τp ₃(r|0)

Wherein p (r|0) represents the probability density function of the high-frequency sub-band coefficient amplitude information that label mask is 0, and r represents the amplitude information that high-frequency sub-band coefficient is corresponding, and τ represents the mutual restricting relation between two, and its value is τ ∈ [0,1], p ₀(r) the cuclear density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 0, p ₃(r|0) exponential probability density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 0.

According to the following formula, asking label mask is the cuclear density function of 0 high-frequency sub-band coefficient amplitude information:

$p_0\left(r\right)=\frac{1}{\mathrm{N\&sigma;}\sqrt{2\mathrm{\&pi;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\exp [-{(r_s-r_j)}^2/{2\mathrm{\&sigma;}}^2]\&times;B[(u_s,v_s),(u_j,v_j)]$

Wherein, p ₀(r) the cuclear density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 0, be made up of gaussian kernel function and box kernel function, r represents the amplitude information that high-frequency sub-band coefficient is corresponding, and N represents the number of high-frequency sub-band coefficient, σ represents the standard deviation of high-frequency sub-band coefficient represent from 1 to N summation, (r _s-r _j) represent the difference of two high-frequency sub-band coefficient amplitude information, s, j represents the position of high-frequency sub-band coefficient, B[(u _s, v _s), (u _j, v _j)] represent box kernel function, can calculate according to the following formula box kernel function.

$B[(u_s,v_s),(u_j,v_j)]=\left\{\begin{array}{cc}1& |u_s-u_j|\&le;L,|v_s-v_j|\&le;L\\ 0& \mathrm{otherwise}\end{array}\right.$

Wherein, (u _s, v _s) and (u _j, v _j) be illustrated respectively in the position s of high-frequency sub-band coefficient and the coordinate at j place, | u _s-u _j| represent the difference of horizontal ordinate, | v _s-v _j| represent the difference of ordinate.

According to the following formula, asking label mask is the exponential probability density function of 0 high-frequency sub-band coefficient amplitude information:

$p_3\left(r\right|0)=\frac{1}{a}{e}^{-r/b}$

Wherein, p ₃(r|0) exponential probability density function of the high-frequency sub-band coefficient amplitude information that expression label mask is 0, r represents the amplitude information that high-frequency sub-band coefficient is corresponding, a, b all represents parameter.

Step 5: ask priori ratio.

From neighbourhood model storehouse, choose 8 different neighbourhood models, than formula, ask respectively the priori ratio of high-frequency sub-band coefficient under different scale according to priori, priori is as follows than formula:

$\mathrm{\&eta;}=\left\{\begin{array}{cc}\max \left[{\mathrm{\&eta;}}_{\mathrm{\&alpha;}}\right],& M\&GreaterEqual;T,\mathrm{\&alpha;}\&Element;\{\mathrm{0,1},...7\}\\ {\mathrm{\&eta;}}_8,& \mathrm{others}\end{array}\right.$

Wherein, η represents the priori ratio of high-frequency sub-band coefficient under different scale, and max represents to get maxima operation, η _αrepresent the priori ratio of lower 8 the different neighbourhood model medium-high frequency sub-band coefficients of different scale, α represents the label of 8 different neighbourhood models, M represents the number of the high-frequency sub-band coefficient that under different scale, label mask is 1, T represents threshold value, its value is 2/3 of selected neighbourhood model medium-high frequency sub-band coefficients, α ∈ { 0,1, ... 7} represents that 8 different neighbourhood model labels are from 0 to 7, η ₈the priori ratio of 3 × 3 square neighbourhood model medium-high frequency sub-band coefficients under expression different scale.

Neighbourhood model is chosen according to following rule:

The first step, from 8 different neighbourhood models, chooses label mask and is a maximum neighbourhood model α of high-frequency sub-band coefficient number of 1, is multiplied by 2/3 by the high-frequency sub-band coefficient number in selected neighbourhood model α, and this product is made as to threshold value;

Second step, with high-frequency sub-band coefficient number and threshold value comparison that in the neighbourhood model α choosing, label mask is 1, if this high-frequency sub-band coefficient number is greater than threshold value, choose this model and ask the priori ratio of high-frequency sub-band coefficient according to priori than formula, if this high-frequency sub-band coefficient number is less than threshold value, choose 3 × 3 square neighbourhood model and ask than formula according to priori the priori ratio of high-frequency sub-band coefficient.

According to the following formula, calculate the priori ratio of 8 different neighbourhood model medium-high frequency sub-band coefficients:

η _α＝exp[κ·T ^α]

Wherein, η _αthe priori ratio that represents 8 different neighbourhood model medium-high frequency sub-band coefficients, κ represents smoothing parameter, T ^αrepresent the priori factor of 8 different neighbourhood models.

According to the following formula, calculate priori factor:

$T^{\mathrm{\&alpha;}}=\left\{\begin{array}{cc}\mathrm{\&Sigma;X}-{\mathrm{\&gamma;}}_1,& X\&Element;\left\{\mathrm{0,1}\right\},\mathrm{\&alpha;}\&Element;\{\mathrm{0,1},...7\}\\ \mathrm{\&Sigma;}{(X-{\mathrm{\&gamma;}}_2)}^3,& X\&Element;\left\{\mathrm{0,1}\right\},\mathrm{\&alpha;}=8\end{array}\right.$

Wherein, T ^αthe priori factor that represents 8 different neighbourhood models, α represents the label of 8 different neighbourhood models, and Σ represents sum operation, and X represents the label mask of high-frequency sub-band coefficient, γ ₁represent parameter, X ∈ 0,1} represent the label mask of high-frequency sub-band coefficient can be 0 or 1, α ∈ 0,1 ..., 7} represents that 8 different neighbourhood model labels are from 0 to 7, γ ₂represent parameter, the label of the square neighbourhood model of α=8 expression 3 × 3.

Step 6: ask the high-frequency sub-band coefficient after spot.

According to the following formula, ask the collapse threshold of high-frequency sub-band coefficient under yardstick l:

ρ＝ηξ/(1+ηξ)

Wherein, ρ represents the collapse threshold of high-frequency sub-band coefficient under yardstick l, and η represents the priori ratio of high-frequency sub-band coefficient under yardstick l, and ξ represents the likelihood ratio of high-frequency sub-band coefficient under yardstick l, l ∈ L.

By high-frequency sub-band multiplication under the collapse threshold of high-frequency sub-band coefficient under yardstick l and yardstick l, obtain the high-frequency sub-band coefficient that yardstick l goes down after spot, l ∈ L.

Step 7: non-lower sampling direction wave inverse transformation.

To low frequency sub-band coefficient and remove the high-frequency sub-band coefficient after spot, carrying out yardstick is the non-lower sampling direction wave inverse transformation of L, obtains the SAR image after spot.

Step 8: output image.

The SAR image after spot is removed in output.

Below in conjunction with the simulated effect figure of accompanying drawing 2 and accompanying drawing 3, the present invention will be further described.

1. emulation experiment condition:

Hardware test platform of the present invention is: processor is Inter Core i3 350M, and dominant frequency is 2.27GHz, internal memory 2GB, and software platform is: Windows7 Ultimate 32-bit operating system and Matlab R2010b.Input picture of the present invention is respectively the SAR image field figure of spot to be gone and the SAR image airport figure of spot to be gone, and size is all 256 × 256, and form is all BMP.

2. emulation content:

Two prior aries as a comparison that the present invention uses are as follows respectively:

The spot method that goes to SAR image that the people such as Lee J S propose in document " Speckle analysis and smoothing of synthetic radar images.Computer Graphics and Image Processing; 1981; 17:24-32. ", is called for short enhanced Lee filtering method.

The spot method that goes to SAR image that the people such as Xiaolin Tian propose in document " Despeckling SAR images based on a new probabilistic model in nonsubsampled contourlet transform domain.Signal; Image and Video Processing; 1-16; 2012 ", is called for short H-NSCT method.

Fig. 2 is that the SAR image field figure that in emulation experiment, the present invention and prior art are treated spot goes spot effect contrast figure.Wherein, Fig. 2 (a) is the spot SAR image field figure to be gone of input, and its size is 256 × 256.Fig. 2 (b) removes the design sketch of spot for the SAR image field figure that adopts enhanced Lee filtering method to treat spot, Fig. 2 (c) removes the design sketch of spot for the SAR image field figure that adopts H-NSCT method to treat spot, and the SAR image field figure that Fig. 2 (d) treats spot for the inventive method removes the design sketch of spot.

Fig. 3 is that the SAR image airport figure that in emulation experiment, the present invention and prior art are treated spot removes the comparison diagram of spot effect.Wherein, Fig. 3 (a) is the spot SAR image airport figure to be gone of input, its size is 256 × 256, Fig. 3 (b) removes the design sketch of spot for the SAR image airport figure that adopts enhanced Lee filtering method to treat spot, Fig. 3 (c) removes the design sketch of spot for the SAR image airport figure that adopts H-NSCT method to treat spot, and the SAR image airport figure that Fig. 3 (d) treats spot for the inventive method removes the design sketch of spot.

3. analysis of simulation result:

Fig. 2 and Fig. 3 be the present invention and prior art treat the SAR image field figure of spot and airport figure go spot effect contrast figure.Can find out from Fig. 2 (b) and Fig. 3 (b), enhanced Lee filtering method of the prior art treat spot SAR image go in the homogeneous region of spot result figure, obviously there are a lot of pseudo-lines, go the SAR image ratio after spot fuzzyyer.Can find out from Fig. 2 (c) and Fig. 3 (c), H-NSCT method of the prior art treat spot SAR image go in spot result figure, in the maintenance of edge and grain details, increase, but homogeneous region still has pseudo-line.Can find out from Fig. 2 (d) and Fig. 3 (d), the inventive method can keep edge and the detailed information of SAR after spot better, and compared to existing technology, the flatness in homogeneous region is better.

In sum, the inventive method can effectively be removed the speckle noise of spot SAR to be gone, can also effectively keep edge and the detailed information of spot SAR image to be gone simultaneously.

To Fig. 2, the spot effect of going that in Fig. 3, each method is treated spot SAR image is carried out objective evaluation, and result is respectively as table 1, shown in table 2.

In general, equivalent number (ENL) is that equivalent number is larger for the module of SAR graphical design, goes that spot effect is unreasonable to be thought; Average (mean) is used for weighing the hold facility of gradation of image value, requires to go after spot the average of image more to approach original image better; Standard deviation (std) is used for weighing the smoothing capability that removes spot method, and standard deviation is less, and its smoothing capability is stronger.

Adopt the present invention and enhanced Lee filtering method of the prior art and H-NSCT method to treat the performance index of going spot effect of spot SAR image field figure as shown in table 1, A, B represents the different spot effect of going twice.

Table 1 goes the performance index of spot effect to field figure

Adopt the present invention and enhanced Lee filtering method of the prior art and H-NSCT method to treat the performance index of going spot effect of spot SAR image airport figure as shown in table 2, A, B represents the different spot effect of going twice.

Table 2 goes the performance index of spot effect to airport figure

From table 1, in table 2, can find out, use the average of removing SAR image after spot that the inventive method obtains close to input picture, show that the inventive method treats the radiation characteristic of spot SAR image and keep best.The standard deviation minimum of removing SAR image after spot that uses the inventive method to obtain, shows that the smoothing capability of the inventive method is best.What use that the inventive method obtains goes after spot in SAR image the equivalent number ENL in homogeneous region the highest, and what show the inventive method goes spot effect best.

Claims (3)

1. the SAR image speckle method based on probability model and NSDT, comprises the steps:

(1) the SAR image of the width spot to be gone of input option;

(2) non-lower sampling direction wave conversion:

The SAR image of the spot to be gone to input carries out the non-lower sampling direction wave conversion that yardstick is L, obtains high-frequency sub-band coefficient and low frequency sub-band coefficient under different scale;

(3) high-frequency sub-band coefficient is asked to label mask:

According to the following formula, obtain the label mask of high-frequency sub-band coefficient under different scale:

$X=\left\{\begin{array}{cc}0,& |D\&CenterDot;y|<{\mathrm{\&sigma;}}^2\\ 1,& |D\&CenterDot;y|\&GreaterEqual;{\mathrm{\&sigma;}}^2\end{array}\right.$

Wherein, X represents the label mask of high-frequency sub-band coefficient under different scale, 0 represents the label mask of spot SAR noise in image to be gone, 1 represents the label mask at edge and details in spot SAR image to be gone, D represents that the SAR image for the treatment of spot carries out the high-frequency sub-band coefficient under the rear yardstick l of non-lower sampling direction wave conversion, the SAR image that y represents to treat spot carries out the estimated value of the high-frequency sub-band coefficient under yardstick l+1 after the conversion of non-lower sampling direction wave, l ∈ L, σ represents the standard deviation of high-frequency sub-band coefficient under different scale;

(4) ask likelihood ratio:

According to the following formula, ask the likelihood ratio of high-frequency sub-band coefficient amplitude information under different scale according to probability distribution:

ξ＝p(r|1)/p(r|0)

Wherein, ξ represents the likelihood ratio of high-frequency sub-band coefficient amplitude information under different scale, p (r|1) represents the probability density function of the high-frequency sub-band coefficient amplitude information that under different scale, label mask is 1, p (r|0) represents the probability density function of the high-frequency sub-band coefficient amplitude information that under different scale, label mask is 0, and r represents amplitude information corresponding to high-frequency sub-band coefficient under different scale;

(5) ask priori ratio:

From neighbourhood model storehouse, choose 8 different neighbourhood models, than formula, ask the priori ratio of high-frequency sub-band coefficient under different scale according to priori, priori is as follows than formula:

$\mathrm{\&eta;}=\left\{\begin{array}{cc}\max \left[{\mathrm{\&eta;}}_{\mathrm{\&alpha;}}\right],& M\&GreaterEqual;T,\mathrm{\&alpha;}\&Element;\{\mathrm{0,1},...7\}\\ {\mathrm{\&eta;}}_8,& \mathrm{others}\end{array}\right.$

Wherein, η represents the priori ratio of high-frequency sub-band coefficient under different scale, and max represents to get maxima operation, η _αrepresent the priori ratio of lower 8 the different neighbourhood model medium-high frequency sub-band coefficients of different scale, α represents the label of 8 different neighbourhood models, M represents the number of the high-frequency sub-band coefficient that under different scale, label mask is 1, T represents threshold value, its value is 2/3 of selected neighbourhood model medium-high frequency sub-band coefficients, α ∈ { 0,1, ... 7} represents that 8 different neighbourhood model labels are from 0 to 7, η ₈the priori ratio of 3 × 3 square neighbourhood model medium-high frequency sub-band coefficients under expression different scale;

(6) ask the high-frequency sub-band coefficient after spot:

(6a) ask according to the following formula, the collapse threshold of high-frequency sub-band coefficient:

ρ＝ηξ/(1+ηξ)

Wherein, ρ represents the collapse threshold of high-frequency sub-band coefficient, and η represents the priori ratio of high-frequency sub-band coefficient, and ξ represents the likelihood ratio of high-frequency sub-band coefficient;

(6b) be multiplied by high-frequency sub-band coefficient by the collapse threshold of high-frequency sub-band coefficient, obtain the high-frequency sub-band coefficient going after spot;

(7) non-lower sampling direction wave inverse transformation:

To low frequency sub-band coefficient and remove the high-frequency sub-band coefficient after spot, carrying out yardstick is the non-lower sampling direction wave inverse transformation of L, obtains the SAR image after spot;

(8) the SAR image after spot is removed in output.

2. the SAR image speckle method based on probability model and NSDT according to claim 1, is characterized in that, the described non-lower sampling direction wave conversion of step (2) is carried out in accordance with the following steps:

The first step, according to the following formula, tectonic transition matrix:

$M_{\mathrm{\&Lambda;}}=\left[\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right]$

Wherein, M _Λrepresent transformation matrix, Λ represents one group of point that the linear combination of the integer vectors of two linear independences forms; a ₁, a ₂, b ₁, b ₂represent respectively transform matrix M _Λfour integer element values;

Second step, chooses respectively 45 degree, 90 degree, 135 three of degree and changes direction, and constructs three transformation matrixs as follows:

$\begin{array}{ccc}\left[\begin{array}{cc}1& 0\\ -1& 1\end{array}\right]& \left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]& \left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\end{array}$

The 3rd step, uses the SAR image of spot to be gone to be multiplied by respectively three transformation matrixs, obtains three cosets;

The 4th step, to each coset, along continuous straight runs carries out one-dimensional wavelet transform twice, vertically carries out one-dimensional wavelet transform one time, obtains seven high-frequency sub-band coefficients and a low frequency sub-band coefficient.

3. the SAR image speckle method based on probability model and NSDT according to claim 1, is characterized in that, the standard deviation of high-frequency sub-band coefficient under the described different scale of step (3), is calculated as follows:

$\mathrm{\&sigma;}=\frac{\mathrm{median}\left(\right|D-\mathrm{median}\left(D\right)\left|\right)}{0.6745}$

Wherein, σ represents the standard deviation of high-frequency sub-band coefficient under different scale, the median operation that represents to average, || represent to take absolute value operation, the SAR image that D represents to treat spot carries out the non-lower sampling direction wave conversion high-frequency sub-band coefficient under different scale afterwards.

Images