CN104574308A - SAR image denoising method based on sampling matrix direction optimization - Google Patents

Abstract

The invention discloses an SAR (synthetic aperture radar) image denoising method based on sampling matrix direction optimization, and belongs to the technical field of synthetic aperture radar image denoising. The directional wavelet transformation of non-downsampling and a Gaussian scale mixture model (GSM) are combined, so that a sampling matrix direction optimization method for a non-downsampling Gaussian scale mixture model (ND-GSM) is provided; binary wavelet transformation is carried out on segmented subgraphs to determine a direction optimization sampling matrix of an SAR image, and the direction optimization sampling matrix is applied to all the subgraphs to form a domain model of a non-downsampling Directionlet domain decomposing coefficient for sampling matrix direction optimization; local denoising on a subgraph transformation domain is carried out by bayesian mean-square estimation, denoised segmented subgraphs are synthesized, and a denoised SAR image is obtained. According to the method, the problem of poor approximation effect of the image is solved; the relevance of inter-domain coefficients can be fully embodied; the method has the advantage in image edge feature maintaining aspect, and the visual effect of the image is improved.

Description

A kind of SAR image denoising method optimized based on sampling matrix direction

Technical field

The invention belongs to diameter radar image noise-removed technology field, be specifically related to a kind of SAR image denoising method of the ND-GSM model based on the optimization of sampling matrix direction.

Background technology

Synthetic-aperture radar (SAR), as the remote sensing sources of a new generation, has round-the-clock, multipolarization, various visual angles and stronger earth's surface loose material penetration capacity, has given play to increasing effect in fields such as military affairs, remote sensing.Coherent speckle noise (speckle) is the intrinsic a kind of determinacy interference of SAR image, therefore, suppresses coherent speckle noise to become an important step of SAR image process.

Speckle noise removal technology looks smoothing technique and the large class of imaging post filtering technology two before can being divided into imaging more, and post-processing technology can be divided into airspace filter and frequency domain filtering again.For the former, as Lee, Gamma-MAP etc. are widely used based on the adaptive filter method of local statistic, and the Speckle reduction ability of this type of wave filter is directly proportional to the size of stationary window, are difficult to the maintenance taking into account the level and smooth of homogeneous area and detailed information.Recent two decades comes, and wavelet analysis achieves and more successfully applies in SAR image process.The people such as nineteen ninety-five Donoho propose to carry out threshold process to reach the object of denoising to wavelet coefficient, and after this wavelet threshold denoising method is widely used in various denoising, and achieves immense success.But the two-dimentional separable wavelets formed by tensor product by one dimension small echo converts (Separable wavelet) only has limited direction, " optimum " can not represents and contain the unusual two dimensional image of line or face.For this vital problem, many scholars propose " multi-scale geometric analysis " (Multiscale Geometric Analysis, MGA) research method.

The Directionalet (direction small echo) proposed at first by people such as Vladan Velisavljevic and Baltasar Beferull Lozano for 2004 is a kind of new multi-scale geometric analysis instrument, based on one dimension small echo, to the expansion of wavelet transformation on multi-direction.It inherits the feature of wavelet transformation separability, according to theory and the integer lattice theory of digital line in computer graphics, constructs multi-direction framework and multi-directional bases.Direction small echo overcomes the shortcoming that standard two-dimensional small echo can not catch image direction characteristic well in image procossing, after being Contourlet (contourlet), based on another new effective image processing tool of discrete domain.

The Directionlet of non-lower sampling converts (Nonsubsampled Directionlet, ND) on the basis of Directionlet conversion, change the filtering mode of coset, no longer carry out down-sampling, thus the redundancy produced makes to have better correlativity between conversion coefficient, and effectively can catch the anisotropic character in image.For each sub-band coefficients after conversion, Gaussian Mixture yardstick (GSM) model can represent its marginal distribution effectively, embodies the strong correlation between neighbour coefficient simultaneously.It is all fixing angle that the direction that non-lower sampling Directionlet converts is chosen, and does not convert according to the feature of image, has cut and distortion phenomenon in graph line edge.Therefore, the sampling matrix direction proposing a kind of non-lower sampling Directionlet territory Gaussian Mixture yardstick (ND-GSM) model is optimized and SAR image denoising, in segmentation subgraph, the non-lower sampling Directionalet of sampling matrix optimization conversion is combined with GSM method, the neighbourhood model of the ND-GSM coefficient of dissociation that structure sampling matrix direction is optimized, utilize Bayes least mean-square estimate antithetical phrase figure transform domain to carry out local denoising, finally synthesize the segmentation subgraph after denoising.Relative experimental results shows, the method effectively can suppress coherent speckle noise, more intactly keeps the detailed information such as the edge in image, embodies the correlativity of coefficient between neighborhood more fully, obtain good visual effect, its integrated performance index is all better than the additive methods such as small echo.

Summary of the invention

The object of the invention is the denoising method that a kind of new SAR image will be provided.This method solve when in the direction and image of the Directionlet basis function of non-lower sampling, anisotropy target is inconsistent, poor to the Approximation effect of image, even deteriorate to the phenomenon of small echo.

First adopt the sampling matrix that the direction optimization method investigation of the sampling matrix of non-lower sampling Directionlet conversion is optimized, its feature comprises following steps:

Step 1: geometrical plane segmentation is carried out to image, require the horizontal resolution of divided image equal with vertical resolution and be all 64 integral multiple, this integral multiple numerical value is the number of times of segmentation, and the subgraph size Pixels size after segmentation should be 64 × 64;

Step 2: carry out dyadic wavelet transform to each subgraph, obtains the detail view v (i, j) of horizontal direction detail view h (i, j) and vertical direction;

Wherein (i, j) represents the position of dyadic wavelet transform coefficient in current sub-band, i, j=1,2 ..., 64;

Step 3: according to h (i, j) and v (i, j), computed segmentation subgraph is in direction θ (i, j) at (i, j) place, and method is as follows:

Step 3.1: if v (i, j) > > h (i, j);

Step 3.2: θ (i, j)=0, if h (i, j) > > v (i, j);

Step 3.3: $\mathrm{\θ}(i,j)=\arctan \frac{v(i,j)}{h(i,j)},$ Other situation;

Step 4: count the value of each segmentation subimage in direction θ (i, j) at (i, j) place by step 3, find out the both direction that occurrence number is maximum, be designated as θ ₁and θ ₂,

If multiple directions are equal during statistics occurrence number, computing method are as follows:

θ _i＝(θ ₁+θ ₂+...+θ _n)/n

Wherein θ ₁, θ ₂... θ _nthe multiple directions that occurrence number is identical, θ _ithen be designated as the direction of this occurrence number.

If there is the Main way of two adjacent sub-images basically identical, computing method are as follows:

Wherein θ ₁, θ ₂the Main way of one of them segmentation subgraph, θ ' ₁with θ ' ₂be the Main way of segmentation subgraph adjacent with it, then these two adjacent subimages carried out merging statistics in the value in direction θ (i, j) at (i, j) place, again find out Main way θ ₁and θ ₂,

Step 5: obtain approximate reasonable slope r by Main way ₁and r ₂, r ₁=arctan θ ₁=b ₁/ a ₁, r ₂=arctan θ ₂=b ₂/ a ₂, according to r ₁and r ₂the sampling matrix of the direction optimization of structure non-lower sampling Directionlet conversion:

Wherein, a ₁, a ₂, b ₁and b ₂all positive integer, wherein along slope r ₁direction be called and change direction, along slope r ₂direction be called queue direction.

Step 6: carrying out sampling matrix with queue direction to each subimage along changing direction is M _∧oblique Anisotropic Wavelet Transform S-AWT (n ₁, n ₂), obtain (ρ=n that entire image anisotropy rate is ρ ₁/ n ₂) non-lower sampling Directionlet convert coefficient;

After the Directionlet decomposition of the non-lower sampling that SAR image is optimized through sampling matrix, the neighborhood y of its observed differential can be as follows with GSM model representation:

$y=x+w=\sqrt{z}u+w---\left(1\right)$

W is zero-mean gaussian vector, and corresponding covariance matrix is C _w, suppose C _wall neighborhoods of same subband are kept constant;

Under condition z, obtaining observed differential neighborhood covariance by formula (1) is C _y|z=zC _u+ C _w, because stochastic variable z, u, w are separate, z is got and expects to substitute into, obtain C _y=E{z}C _u+ C _w, E{z}=1 is set, then:

C _u＝C _y-C (2)

Image obtains the sub-band coefficients of each yardstick through decomposing, suppose coefficient x _cfield coefficient x around meets GSM model, then random vector x can be expressed as zero-mean gaussian vector u and independent positive yardstick random factor product: "=" number expression has identical distribution, and factor z is called weight coefficient, and the probability density of vector x is by the covariance matrix C of u _uwith coefficient probability density p _zz () determined:

N is the dimension (being herein the size of neighborhood) of x and u. establish E{z}=1, then C _x=C _u.

$\begin{array}{c}{p}_x\left(x\right)=\∫p\left(x\right|z)p_z\left(z\right)\mathrm{dz}\\ =\∫\frac{\exp \left(\frac{-x^T{\left(z{C}_u\right)}^{-1}x}{2}\right)}{{\left(2\mathrm{\π}\right)}^{N/2}{\left|{\mathrm{zC}}_u\right|}^{1/2}}{p}_z\left(z\right)\mathrm{dz}\end{array}---\left(3\right)$

The neighborhood of N × N size is used to estimate it, C _ycovariance for observed differential neighborhood in neighborhood:

C _y＝E{(y-μ _u)·(y-μ _u) ^T} (4)

Wherein μ _y=E{y} represents the expectation value of y;

Field of noise covariance C _wby analytic function:

$\mathrm{\σ}\sqrt{N_y{N}_x}\mathrm{\δ}(n,m)---\left(5\right)$

Wherein (N _y, N _x) be image size, this δ signal has identical power spectrum with noise signal;

Use Bayes's least mean-square estimate: $E\left\{x_c\right|y\}=\underset{0}{\overset{\∞}{\∫}}p\left(z\right|y)E\{x_c|y,z\}\mathrm{dz}---\left(6\right)$

The inventive method is used for the center coefficient x under design conditions z _cthe average of the Bayes's least mean-square estimate being weight with posterior density p (z|y), the key property of GSM model is, the neighborhood x of coefficient vector is the gaussian variable under condition z, utilize the character of additive white Gaussian noise, formula (4) expects that being a simple Wei Na estimates:

E{x|y,z}＝zC _u(zC _u+C _w) ^-1y (7)

Use matrix zC _u+ C _wdiagonalization reduce the dependence of above formula to z;

Described denoising method, is characterised in that, the method can adapt to the Main way of anisotropy target in image, and denoising method comprises following steps:

Step 2.1: carry out log-transformation to original SAR image, makes it meet additive noise hypothesis;

Step 2.2: the non-lower sampling Directionlet carrying out the optimization of sampling matrix direction converts, and method is as follows:

Step 2.2.1: compartition is carried out to image, determine the number of times of segmentation according to the size of image, subgraph size should be 64 × 64;

Step 2.2.2: carry out dyadic wavelet transform to each subgraph, finally determines the sampling matrix M optimized in direction _∧;

Step 2.2.3: carrying out sampling matrix to each segmentation subgraph is M _∧sampling, obtains | det (M _∧) | (M _∧the absolute value of determinant) individual coset, each coset corresponds to a displacement vector S _k, and each coset is by displacement vector S _kdetermined;

Step 2.2.4: carry out non-lower sampling Directionlet conversion respectively to each coset of each segmentation subgraph, obtains each segmentation subgraph corresponding high and low frequency coefficient subband;

Step 2.3: utilize GSM model to estimate each each sub-band coefficients noise of segmentation subgraph except low frequency;

Step 2.3.1: according to picture noise standard deviation, calculates neighborhood noise covariance C _w, the covariance C of estimation neighbourhood coefficient _y;

Step 2.3.2: utilize C _yand C _westimate C _u;

Step 2.3.3: simplifying E{x|y, z} is that local Wei Na estimates;

Step 2.3.4: utilize Bayes's least mean-square estimate to calculate its center coefficient x to each neighborhood in subband _c;

Step 2.4: carry out non-lower sampling Directionlet inverse transformation to low frequency sub-band with through the high-frequency sub-band of filtering process;

Step 2.5: be weighted comprehensively according to the direction of displacement vector corresponding to the coset selected, reconstructs each segmentation subgraph;

Step 2.6: by each segmentation subgraph of reconstruct by its position synthesis in original image, obtain the image after denoising.

Accompanying drawing explanation

Fig. 1 sampling matrix direction Optimizing Flow;

Fig. 2 denoising method;

Fig. 3 a-f field SAR image denoising result;

Fig. 4 field SAR image denoising details 1 amplification effect;

Fig. 5 bay SAR image denoising effect;

Fig. 6 base SAR image denoising effect.

Embodiment

In order to verify the validity superiority of the inventive method, the method that contrast experiment chooses comprises: Lee filtering, Gamma-MAP filtering, wavelet soft thresholding, DT-Bayes (± 45 °) threshold method, wherein Lee and Gamma-MAP wave filter all adopts the filter window of 5 × 5 sizes, wavelet method chooses db9-7 as wavelet basis, and the threshold value that wavelet soft-threshold method is chosen is σ represents that noise criteria is poor.

In the inventive method, anisotropy wavelet basis selects AWT (2,1), although Directionlet can construct multiple different changing direction, but correspondingly can increase computation complexity, therefore experiment is chosen 2 and is changed direction, be respectively 45 ° and-45 ° of directions, corresponding sampling matrix is: $M_1=\left[\begin{array}{cc}1& 1\\ -1& 1\end{array}\right],$ $M_2=\left[\begin{array}{cc}-1& 1\\ 1& 1\end{array}\right].$

Fig. 3 gives the inventive method with other control methodss to the result of field image denoising.As can be seen from Fig. 3 (b), Fig. 3 (c), Lee filtering and Gamma-MAP filtering are carried out level and smooth preferably to the speckle noise of image, but obscurity boundary after Lee filtering and Gamma-MAP filtering as can be seen from Figure 4 and Figure 5, image detail information is lost serious.Fig. 3 (d) is the result of wavelet soft thresholding respectively, smoother, but image detail part is seriously lost, and illustrates that noise filtering is limited in one's ability.Fig. 3 (e) is the image after the process of DT-Bayes threshold method, homogeneous area is more level and smooth, but due to Directionlet change direction choose be ± 45 °, do not convert according to the feature of image, have cut and distortion phenomenon in graph line edge.Fig. 3 (f) carries out the result of speckle suppression for ND-GSM method that sampling matrix direction is optimized, effective filtering speckle noise, the texture at edge, field and edge is more clear, other homogeneous area is more level and smooth, distortion is less, and there is higher Y-PSNR, improve overall sharpness, in picture quality and visual effect, all have larger improvement.

Fig. 4 and Fig. 5 is the details enlarged drawing after the image filtering of field, contrast finds that wavelet method exists significantly " ringing effect " in edge, the inventive method then maintains the minutias such as the edge of image preferably, improves the flatness in image uniform region simultaneously.

Table 1 lists said method to the filtered property indices of field SAR image.ENLBLOCK1 and ENLBLOCK2 represents the ENL value of 2 pieces of homogeneous areas selected in Fig. 3 (a) respectively.By contrast images homogeneous area and overall filtering performance, the inventive method is all better than wavelet thresholding method and the GSM method based on small echo.The overall equivalent number of Lee filtering and Gamma-MAP filtering is a little more than the inventive method, this is to caused by speckle noise excess smoothness, observes Fig. 3 (a), (b) can find that it visually reduces. Fig. 5 (b), (c) and Fig. 6 (b), (c) sets forth DT-Bayes threshold method and the inventive method denoising result to bay and base SAR image.The ND-GSM method that sampling matrix direction is optimized all is better than additive method in the raising of ENL and the improvement two of picture quality, again demonstrates its validity.

Table 1 lists various method carries out denoising performance index to field SAR image, can find out, in most of the cases, DT-Bayes denoising result is better than small echo, but also there is the situation of too late small echo, illustrating when setting arbitrarily image conversion direction, Directionlet better performance more excellent in small echo can not be ensured.Under same noise level, the performance index obtained after the Directionlet transform method image denoising that sampling matrix direction of the present invention is optimized, substantially all higher than other method, achieve better visual effect.

Last it is noted that above example only in order to illustrate the present invention and and unrestricted technical scheme described in the invention; Therefore, although this instructions with reference to above-mentioned example to present invention has been detailed description, those of ordinary skill in the art should be appreciated that and still can modify to the present invention or equivalent to replace; And all do not depart from technical scheme and the improvement thereof of the spirit and scope of invention, it all should be encompassed in the middle of right of the present invention.

The each method of table 1 contrasts field SAR image denoising performance.

Claims (1)

1., based on the SAR image denoising method that sampling matrix direction is optimized, adopt the sampling matrix that the direction optimization method investigation of the sampling matrix of non-lower sampling Directionlet conversion is optimized, its feature comprises the following steps:

Step 1.1: geometrical plane segmentation is carried out to image, require the horizontal resolution of divided image equal with vertical resolution and be all 64 integral multiple, this integral multiple numerical value is the number of times of segmentation, and the subgraph size Pixels size after segmentation should be 64 × 64;

Step 1.2: carry out dyadic wavelet transform to each subgraph, obtains the detail view v (i, j) of horizontal direction detail view h (i, j) and vertical direction;

Wherein (i, j) represents the position of dyadic wavelet transform coefficient in current sub-band, i, j=1,2 ..., 64;

Step 1.3: according to h (i, j) and v (i, j), computed segmentation subgraph is in direction θ (i, j) at (i, j) place, and method is as follows:

Step 1.3.1: if v (i, j) > > h (i, j);

Step 1.3.2: θ (i, j)=0, if h (i, j) > > v (i, j);

Step 1.3.3: other situation;

Step 1.4: count the value of each segmentation subimage in direction θ (i, j) at (i, j) place by step 1.3, find out the both direction that occurrence number is maximum, be designated as θ ₁and θ ₂,

If multiple directions are equal during statistics occurrence number, computing method are as follows:

θ _i＝(θ ₁+θ ₂+...+θ _n)/n

Wherein θ ₁, θ ₂θ _nthe multiple directions that occurrence number is identical, θ _ithen be designated as the direction of this occurrence number;

If there is the Main way of two adjacent sub-images basically identical, computing method are as follows:

Wherein θ ₁, θ ₂the Main way of one of them segmentation subgraph, θ ' ₁with θ ' ₂be the Main way of segmentation subgraph adjacent with it, then these two adjacent subimages carried out merging statistics in the value in direction θ (i, j) at (i, j) place, again find out Main way θ ₁and θ ₂,

Step 1.5: obtain approximate reasonable slope r by Main way ₁and r ₂, r ₁=arctan θ ₁=b ₁/ a ₁, r ₂=arctan θ ₂=b ₂/ a ₂, according to r ₁and r ₂the sampling matrix of the direction optimization of structure non-lower sampling Directionlet conversion:

Wherein, a ₁, a ₂, b ₁and b ₂all positive integer, wherein along slope r ₁direction be called and change direction, along slope r ₂direction be called queue direction;

Step 1.6: carrying out sampling matrix with queue direction to each subimage along changing direction is M _∧oblique Anisotropic Wavelet Transform S-AWT (n ₁, n ₂), obtain (ρ=n that entire image anisotropy rate is ρ ₁/ n ₂) non-lower sampling Directionlet convert coefficient;

After the Directionlet decomposition of the non-lower sampling that SAR image is optimized through sampling matrix, the neighborhood y of its observed differential can be as follows with GSM model representation:

W is zero-mean gaussian vector, and corresponding covariance matrix is C _w, suppose C _wall neighborhoods of same subband are kept constant;

Under condition z, show that observed differential neighborhood covariance is C _y|z=zC _u+ C _w, because stochastic variable z, u, w are separate, z is got and expects to substitute into, obtain C _y=E{z}C _u+ C _w, E{z}=1 is set, then:

C _u＝C _y-C _w

Image obtains the sub-band coefficients of each yardstick through decomposing, suppose coefficient x _cfield coefficient x around meets GSM model, then random vector x can be expressed as zero-mean gaussian vector u and independent positive yardstick random factor product: "=" number expression has identical distribution, and factor z is called weight coefficient, and the probability density of vector x is by the covariance matrix C of u _uwith coefficient probability density p _zz () determined:

N is the dimension (being herein the size of neighborhood) of x and u. establish E{z}=1, then C _x=C _u.

The neighborhood of N × N size is used to estimate it, C _ycovariance for observed differential neighborhood in neighborhood:

C _y＝E{(y-μ _u) ·(y-μ _u) ^T}

Wherein μ _y=E{y} represents the expectation value of y;

Field of noise covariance C _wby analytic function:

Wherein (N _y, N _x) be image size, this δ signal has identical power spectrum with noise signal;

Use Bayes's least mean-square estimate:

The inventive method is used for the center coefficient x under design conditions z _cthe average of the Bayes's least mean-square estimate being weight with posterior density p (z|y), the key property of GSM model is, the neighborhood x of coefficient vector is the gaussian variable under condition z, utilize the character of additive white Gaussian noise, expect that being a simple Wei Na estimates:

E{x|y,z}＝zC _u(zC _u+C _w) ^-1y

Use matrix zC _u+ C _wdiagonalization reduce the dependence of above formula to z;

Described denoising method, is characterized in that, the method can adapt to the Main way of anisotropy target in image, and denoising method comprises the steps:

Step 2.1: carry out log-transformation to original SAR image, makes it meet additive noise hypothesis;

Step 2.2: the non-lower sampling Directionlet carrying out the optimization of sampling matrix direction converts, and method is as follows:

Step 2.2.1: compartition is carried out to image, determine the number of times of segmentation according to the size of image, subgraph size should be 64 × 64;

Step 2.2.2: carry out dyadic wavelet transform to each subgraph, finally determines the sampling matrix M optimized in direction _∧;

Step 2.2.3: carrying out sampling matrix to each segmentation subgraph is M _∧sampling, obtains | det (M _∧) | (M _∧the absolute value of determinant) individual coset, each coset corresponds to a displacement vector S _k, and each coset is by displacement vector S _kdetermined;

Step 2.2.4: carry out non-lower sampling Directionlet conversion respectively to each coset of each segmentation subgraph, obtains each segmentation subgraph corresponding high and low frequency coefficient subband;

Step 2.3: utilize GSM model to estimate each each sub-band coefficients noise of segmentation subgraph except low frequency;

Step 2.3.1: according to picture noise standard deviation, calculates neighborhood noise covariance C _w, the covariance C of estimation neighbourhood coefficient _y;

Step 2.3.2: utilize C _yand C _westimate C _u;

Step 2.3.3: simplifying E{x|y, z} is that local Wei Na estimates;

Step 2.3.4: utilize Bayes's least mean-square estimate to calculate its center coefficient x to each neighborhood in subband _c;

Step 2.4: carry out non-lower sampling Directionlet inverse transformation to low frequency sub-band with through the high-frequency sub-band of filtering process;

Step 2.5: be weighted comprehensively according to the direction of displacement vector corresponding to the coset selected, reconstructs each segmentation subgraph;

Step 2.6: by each segmentation subgraph of reconstruct by its position synthesis in original image, obtain the image after denoising.