CN103473755B - Based on the sparse denoising method of SAR image that change detects - Google Patents

Abstract

The invention discloses a kind of sparse denoising method of SAR image detected based on change, implementation step is: (1) input picture; (2) region of variation image is obtained; (3) setting change class; (4) two classification binary map are obtained; (5) two area images are obtained; (6) not region of variation denoising; (7) estimating noise standard deviation; (8) composite diagram is obtained; (9) sparse dictionary denoising; (10) Output rusults.The present invention be directed to the noise suppression preprocessing that SAR image change detects, have region of variation holding consistency better, well can suppress the noise of not region of variation, retain SAR image particularly detects impact region of variation greatly texture, profile, edge detail information on later stage change simultaneously, the cut produced after removal of images denoising, reduce the fake information because squelch produces, improve the advantage of the precision that change detects.

Description

SAR image sparse denoising method based on change detection

Technical Field

The invention belongs to the technical field of image processing, and further relates to a Synthetic Aperture Radar (SAR) image sparse denoising method based on change detection in the field of change detection. The method can be used for carrying out denoising pretreatment on SAR image change detection, and can effectively improve the SAR image change detection precision.

Background

The image denoising technology solves the problem of image quality reduction caused by various noise interferences on the image in the processes of acquisition, encoding, transmission and the like, improves the image quality, and is an important link and research content in image processing. The SAR image denoising method can keep the detail information of the image on the basis of effectively inhibiting speckle noise.

Generally, the processing after denoising is segmentation, classification, and the like of the SAR image, and is performed for a single image. However, for the denoising preprocessing of the change detection between the SAR images of two or more time phases, because the speckle noises of the imaging are different, in order to prevent the subsequent change detection from causing false detection of the change area and the non-change area due to the difference of texture traces or information loss caused by the denoising results at different phases, the speckle denoising consistency of the non-change area of the SAR images at different phases needs to be maintained.

Before SAR imaging, multi-view smoothing processing can be carried out to inhibit speckle noise, but the multi-view processing improves the signal-to-noise ratio by N/2 times at the cost of reducing the spatial resolution by N times, so that the prior filtering pretreatment method for SAR image change detection after imaging mainly comprises two methods of spatial filtering and frequency filtering. The spatial domain includes mean filtering, median filtering, morphological filtering, etc. These filtering processes smooth the noise of the homogeneous region to some extent, but have a disadvantage that it is difficult to remove the noise of the homogeneous region while maintaining edge information of the image. The frequency domain speckle suppression method is commonly a wavelet hidden Markov tree based denoising method WHMT and the like.

The iterative weight maximum likelihood denoising PPB method is an SAR image denoising method provided in the document' iterative weighted maximum likelihood denoising patch-based weights [ J ]. IEEETransactionson Imageprocessing,2009,18(12): 2661-once 2672 ], which is characterized in that the noise of the homogeneous region can be well smoothed and the noise of the texture detail region can be well removed, but the homogeneous region can generate scratches and easily generate an over-smoothing phenomenon, and pseudo change information is introduced for subsequent change detection.

A sparse representation-based SAR image speckle suppression method is proposed in a patent technology ' SAR image speckle suppression method based on sparse representation ' (patent application number: 201110346349.7, publication number: CN 102346908B) owned by the university of Western's electronics science and technology. When the technology is used for SAR image change detection denoising preprocessing, sparse approximation of an image on a redundant dictionary is adopted to achieve denoising. Although the method can inhibit noise of homogeneous regions and simultaneously maintain detailed information, the method still has the defects that the error control in dictionary learning is in actual operation, partial texture information of the image is easy to lose, and the false detection rate of later-stage change detection is increased.

Disclosure of Invention

The invention aims to overcome the defects of the prior art and provides a SAR image sparse denoising method based on change detection. The method can effectively remove the noise of the homogeneous region, and can well reserve the detail information such as texture, contour, edge and the like of the change region which has larger influence on the SAR image change detection.

The technical idea for realizing the invention is as follows: firstly, finding out the change area of two time phase images, and secondly, dividing each time phase image into an unchanged area image and a changed area image; then, denoising the invariant region image by using a wavelet hidden Markov WHMT; and finally, denoising the composite image by using a sparse dictionary method to obtain a denoising result image of the two-time phase image.

The specific steps for realizing the purpose of the invention are as follows:

(1) inputting an image:

1a) selecting two single-polarized SAR images which are generated by satellites at the same place and different time and have the same registered size;

1b) reading a time phase image 1 in a first single-polarized SAR image matrix;

1c) a phase image 2 in a second single-polarized SAR image matrix is read in.

(2) Obtaining a change area image:

2a) taking difference between the two time phase SAR images, and taking an absolute value of the difference operation result to obtain a difference image;

2b) carrying out mean filtering on the difference image by using a template with the window size of 3 multiplied by 3 to obtain a filtered difference image;

2c) calculating according to the following formula to obtain the average energy of any pixel point in the filtered difference image,

$E=\frac{1}{9}\underset{q\∈B}{\mathrm{\Σ}}{g}_q^2$

wherein E represents the average energy of any pixel point, q represents the pixel serial number in the neighborhood of 3 × 3 pixels with the pixel point as the center, q =1,2, …, 9, B represents the neighborhood of 3 × 3 pixels with the pixel point as the center, g_qRepresenting the gray value of the q-th pixel in the neighborhood B;

2d) replacing the gray value of the pixel point in each difference image with the average energy of the pixel point in each difference image to obtain an average energy map;

2e) calculating the gray level mean value of the average energy map, setting the mean value as a threshold value, assigning the gray level value which is greater than or equal to the threshold value as 1, and assigning the gray level value which is smaller than the threshold value as 0 to obtain a binary average energy map;

2f) performing morphological erosion operation on the binary average energy map by using 3 x 3 structural elements to obtain a change detection image, wherein a region with a pixel value of 1 in the change detection image represents a change region of the detected two-time phase image;

2g) the gray values of pixels out of the spatial positions corresponding to the change areas of the two time phase images in the time phase image 1 and the time phase image 2 are respectively set to 0, and the initial change area image of the time phase image 1 and the initial change area image of the time phase image 2 are obtained.

(3) Setting a change class:

3a) subtracting the initial change area image of the time phase image 2 from the initial change area image of the time phase image 1 to obtain an initial change area difference image;

3b) setting the pixel value in the difference image of the initial change area to be less than or equal to the set threshold value T₁The pixel points of (1) are set to be in a non-change class; the pixel value in the difference image of the initial change area is larger than the threshold value T₁The pixel points are set as variation classes;

3c) setting the pixel value in the initial change area of the time phase image 1 corresponding to the spatial position of the pixel point set as the change class in the initial change area difference image as 1; and setting the pixel value in the initial change area of the time phase image 2 corresponding to the spatial position of the pixel point set as the change class in the initial change area difference image as 1.

(4) Obtaining two classification binary images:

4a) counting the number of pixels which are larger than 0 in a neighborhood of 3 multiplied by 3 and taking the pixel as the center for each pixel which belongs to a non-change type in the difference image of the initial change area;

4b) comparing the number of pixels greater than 0 in the neighborhood with a threshold T₂、T₃The following operations are performed to obtain two classification binary maps of the time phase image 1 and two classification binary maps of the time phase image 2:

if the number of pixels larger than 0 in the neighborhood is larger than or equal to the threshold value T₂Setting the value of the pixel point of the initial change area image of the time phase image 1 corresponding to the spatial position of the non-change pixel satisfying the condition as 1;

if the number of pixels larger than 0 in the neighborhood is less than the threshold value T₃Setting the value of the pixel point of the initial change area image of the time phase image 2 corresponding to the spatial position of the non-change pixel satisfying the condition as 1;

if the number of pixels larger than 0 in the neighborhood is less than the threshold value T₂And is greater than or equal to threshold T₃Setting the values of the pixel points of the initial change area image of the time phase image 1 and the initial change area image of the time phase image 2 corresponding to the spatial position of the pixel of the non-change class meeting the condition as 1;

the value of the unmarked pixel points in the initial change area image of the time phase image 1 and the initial change area image of the time phase image 2 is set to 0.

(5) Two area images are obtained:

5a) multiplying the two classification binary images of the time phase image 1 by pixel points of pixel values corresponding to the spatial position of the image of the time phase image 1 to obtain a change area image of the time phase image 1; subtracting the change area image of the time phase 1 image from the time phase image 1 to obtain an unchanged area image of the time phase image 1;

5b) multiplying the two classified binary images of the time phase image 2 by pixel points of pixel values corresponding to the space position of the time phase image 2 to obtain a change area image of the time phase image 2; the change area image of the phase image 2 is subtracted from the phase image 2 to obtain the unchanged area image of the phase image 2.

(6) Denoising the unchanged area:

denoising the invariant region images of the time phase image 1 and the time phase image 2 respectively by using a Wavelet Hidden Markov Tree (WHMT) method to obtain denoised images of the invariant regions of the time phase image 1 and the time phase image 2.

(7) Estimating the standard deviation of noise:

7a) obtaining the noise standard deviation of the unchanged area image of the two time phase images according to the following formula

$\mathrm{\σ}=\frac{1}{Z}\underset{Q}{\mathrm{\Σ}}(v_1-v_2)-\frac{1}{Z}\underset{Q}{\mathrm{\Σ}}{(v_1-v_2)}^2$

Wherein σ represents a noise standard deviation of the invariant region image of the time phase image, Z represents a total number of pixel points in the invariant region image of the time phase image, Q represents the invariant region image of the time phase image, v represents a noise standard deviation of the invariant region image of the time phase image, and₁gray value, v, of pixel points of an invariant region image representing a time phase image₂Representing sum v in denoised image of invariant regions of phase image₁And gray values of pixel points corresponding to the spatial positions.

(8) Obtaining a synthetic graph:

and respectively adding the denoised image of the invariant region of each time phase image and the pixel value of the space position corresponding to the image of the variant region of the time phase image to obtain a composite image of the time phase image 1 and a composite image of the time phase image 2.

(9) Denoising the sparse dictionary:

9a) respectively initializing over-complete sparse representation dictionaries of the two time phase images, wherein the over-complete sparse representation dictionaries are selected as Discrete Cosine Transform (DCT) dictionaries;

9b) image blocks of a composite image of two time-phase images are obtained according to the following formula:

I_ij=R_ijI

wherein, I_ijImage blocks representing a composite image, i, j representing the reference numerals of the rows and columns of the image blocks, R_ijRepresenting an image decimation matrix, I representing a composite image;

9c) judging whether the noise standard deviation of the image of the unchanged area is less than or equal to 5, if so, executing the step 9d) for 5 times; otherwise, executing step 9d)10 times;

9d) adopting a dictionary design method KSVD method of sparse representation, respectively updating an overcomplete sparse representation dictionary and an image block I of a synthesized image relative to a time phase image 1 and a time phase image 2_ijSparse representation coefficients of (a);

9e) and respectively estimating the denoising result graph of the time phase image by adopting a reconstruction formula to obtain a denoising result graph of the time phase image 1 and a denoising result of the time phase image 2.

(10) And outputting a result:

and outputting the denoising result graphs of the phase image 1 and the phase image 2.

Compared with the prior art, the invention has the following advantages:

first, the invention detects the change of two time phase images by de-noising pretreatment aiming at SAR image change detection, and overcomes the problem that the prior art introduces pseudo information into a change area because the two time phase images have inconsistent noise, so that the invention has the advantage of better keeping the consistency of the change area.

Secondly, the invention denoises the invariant region by using the wavelet hidden Markov method, overcomes the problem that the edge texture information which has larger influence on the change detection result is lost when the noise is suppressed in the prior art, and has the advantages of well suppressing the noise of the invariant region and simultaneously reserving the texture, contour and edge detail information of the SAR image, particularly the change region which has larger influence on the later change detection.

Thirdly, the invention denoises the synthetic image by using the sparse dictionary, overcomes the problems of scratches generated by image denoising and introduction of pseudo change information in the prior art, and has the advantages of eliminating scratches generated after image denoising, reducing pseudo information generated by noise suppression and improving the precision of change detection.

Drawings

FIG. 1 is a flow chart of the present invention;

FIG. 2 is a two-phase image input by the present invention;

FIG. 3 shows the denoising result of the two-time phase image by the method of the present invention and the PPB method;

FIG. 4 is a logarithmic ratio plot of denoised two-phase images of the method of the present invention and the PPB method;

Detailed Description

The invention is further described below with reference to the accompanying drawings.

Referring to fig. 1, the steps implemented by the present invention are as follows:

step 1, inputting an image.

And selecting two satellites to generate registered single-polarized SAR images with the same size at the same place and different time.

A phase image 1 in a first single-polarized SAR image matrix is read in.

A phase image 2 in a second single-polarized SAR image matrix is read in.

And 2, obtaining a change area image.

And (4) performing difference on the SAR images in the two time phases, and taking an absolute value of the difference operation result to obtain a difference image.

And carrying out average filtering on the difference image by using a template with the window size of 3 multiplied by 3 to obtain the filtered difference image.

Calculating according to the following formula to obtain the average energy of any pixel point in the filtered difference image,

$E=\frac{1}{9}\underset{q\∈B}{\mathrm{\Σ}}{g}_q^2$

wherein E represents the average energy of any pixel point, q represents the pixel serial number in the neighborhood of 3 × 3 pixels with the pixel point as the center, q =1,2, …, 9, B represents the neighborhood of 3 × 3 pixels with the pixel point as the center, g_qRepresenting the gray value of the q-th pixel within the neighborhood B.

And replacing the gray value of the pixel point in each difference image by the average energy of the pixel point in each difference image to obtain an average energy map.

And calculating the gray level mean value of the average energy map, setting the mean value as a threshold value, assigning the gray level values which are greater than or equal to the threshold value as 1, and assigning the gray level values which are smaller than the threshold value as 0 to obtain a binary average energy map.

Morphological erosion operation is performed on the binarized average energy map with 3 × 3 structural elements to obtain a change detection image in which a region having a pixel value of 1 indicates a change region of the detected two-time phase image.

The gray values of pixels out of the spatial positions corresponding to the change areas of the two time phase images in the time phase image 1 and the time phase image 2 are respectively set to 0, and the initial change area image of the time phase image 1 and the initial change area image of the time phase image 2 are obtained.

And 3, setting a variation class.

The initial change area image of the time phase image 2 is subtracted from the initial change area image of the time phase image 1 to obtain an initial change area difference image.

Setting the pixel value in the difference image of the initial change area to be less than or equal to the set threshold value T₁The pixel points of (1) are set to be in a non-change class; the pixel value in the difference image of the initial change area is larger than the threshold value T₁Is set to a variation class, threshold T₁The value range of (A) is 16-24, and the threshold value T of the invention₁Is 20.

Setting the pixel value in the initial change area of the time phase image 1 corresponding to the spatial position of the pixel point set as the change class in the initial change area difference image as 1; and setting the pixel value in the initial change area of the time phase image 2 corresponding to the spatial position of the pixel point set as the change class in the initial change area difference image as 1.

And 4, obtaining two classification binary images.

For each pixel point belonging to the non-change class in the initial change region difference image, the number of pixels larger than 0 is counted in a neighborhood of 3 × 3 with the pixel point as the center.

Comparing the number of pixels greater than 0 in the neighborhood with a threshold T₂、T₃The following operations are performed to obtain two classification binary maps of the time phase image 1 and two classification binary maps of the time phase image 2:

if the number of pixels larger than 0 in the neighborhood is larger than or equal to the threshold value T₂Then, the value of the pixel point of the initial change region image of the time phase image 1 corresponding to the spatial position of the pixel of the non-change class satisfying this condition is set to 1.

If the number of pixels larger than 0 in the neighborhood is less than the threshold value T₃Then, the value of the pixel point of the initial change region image of the time phase image 2 corresponding to the spatial position of the pixel of the non-change class satisfying this condition is set to 1.

If the number of pixels larger than 0 in the neighborhood is less than the threshold value T₂And is greater than or equal to threshold T₃Then, the values of the pixels of the initial change region image of the time phase image 1 and the initial change region image of the time phase image 2 corresponding to the spatial position of the pixel of the non-change class satisfying this condition are both set to 1.

Setting the value of unmarked pixel points in the initial change area image of the time phase image 1 and the initial change area image of the time phase image 2 as 0, and setting the threshold value T₂Is in the range of 4 to 7, and has a threshold value T₃The value range of (A) is more than or equal to T₃≤T₂Of the invention T₂And T₃The values of (a) and (b) are 6 and 4, respectively.

And 5, obtaining two-region images.

Multiplying the two classification binary images of the time phase image 1 by pixel points of pixel values corresponding to the spatial position of the image of the time phase image 1 to obtain a change area image of the time phase image 1; the change area image of the phase 1 image is subtracted from the phase image 1 to obtain the unchanged area image of the phase image 1.

Multiplying the two classified binary images of the time phase image 2 by pixel points of pixel values corresponding to the space position of the time phase image 2 to obtain a change area image of the time phase image 2; the change area image of the phase image 2 is subtracted from the phase image 2 to obtain the unchanged area image of the phase image 2.

And 6, denoising the unchanged area.

Denoising the invariant region images of the time phase image 1 and the time phase image 2 respectively by using a Wavelet Hidden Markov Tree (WHMT) method to obtain denoised images of the invariant regions of the time phase image 1 and the time phase image 2.

The first step is that 3-layer wavelet transformation is respectively carried out on the unchanged area images of the time phase images, and noisy wavelet coefficient sets of the unchanged area images of the time phase images are respectively obtained.

And secondly, adopting a tree-shaped Markov chain to fit the conditional distribution of the wavelet coefficients among the scales to establish a unidirectional transmission hidden Markov tree model for the noisy wavelet coefficient set obtained after wavelet transformation. Obtaining a hidden Markov tree model parameter set according to the following formula:

$V=\{u_{\mathrm{k\ξ}},{\mathrm{\σ}}_{\mathrm{k\ξ}}^2,p\left(s_k\right|s_{\mathrm{\ρ}\left(k\right)}),p\left(s_0\right)\}$

wherein V represents a parameter set of the hidden Markov tree model, u_kξMeans of the node, σ, representing the wavelet coefficient at the state ξ, k position_kξWavelet system representing ξ, k-position nodeNumber standard deviation, p(s)_k|s_ρ(k)) Probability, s, of k-position node at state ξ when parent node representing k-position node is at state ξ_kState variable, s, representing a node at k position_ρ(k)State variables representing the parent of k-position nodes, [ rho (k) represents the parent of k-position nodes, [ p(s) ]₀) Representing the initial state distribution, s, of the root node₀Representing the state variables of the root node.

And thirdly, training the noisy wavelet coefficient by using an expectation maximization algorithm to obtain an optimal estimated hidden Markov tree model parameter set.

Fourthly, estimating the noiseless wavelet coefficient according to the following formula for the obtained optimal estimated hidden Markov tree model parameter set

$y_k=\underset{\mathrm{\ξ}}{\mathrm{\Σ}}p(s_k=\mathrm{\ξ}|W,V){\mathrm{\σ}}_{\mathrm{k\ξ}}^2/({\mathrm{\σ}}_{\mathrm{k\ξ}}^2+{\mathrm{\σ}}_1^2)$

Wherein, y_kRepresenting the estimated noise-free wavelet coefficients, ξ representing the node states, p(s)_k= ξ W, V) represents the state variable s of the k-position node with the knowledge of the optimal estimated set of hidden markov tree model parameters_kProbability of being state ξ, s_kRepresenting state variables of k-position sub-nodes, and W representing a noisy wavelet coefficient setIn sum, V represents the set of parameters, σ, of the hidden Markov tree model for optimal estimation_kξStandard deviation, σ, of node representing state ξ, k position₁Representing the standard deviation contained in the wavelet coefficients.

A fifth step of estimating the wavelet coefficients y_kAnd respectively carrying out wavelet inverse transformation to obtain the denoised image of the invariant region image of the time phase image.

And 7, estimating the noise standard deviation.

Obtaining the standard deviation of the noise of the unchanged area image of the two time phase images according to the following formula:

$\mathrm{\σ}=\frac{1}{Z}\underset{Q}{\mathrm{\Σ}}(v_1-v_2)-\frac{1}{Z}\underset{Q}{\mathrm{\Σ}}{(v_1-v_2)}^2$

wherein σ represents a noise standard deviation of the invariant region image of the time phase image, Z represents a total number of pixel points in the invariant region image of the time phase image, Q represents the invariant region image of the time phase image, v represents a noise standard deviation of the invariant region image of the time phase image, and₁gray value, v, of pixel points of an invariant region image representing a time phase image₂Representing sum v in denoised image of invariant regions of phase image₁And gray values of pixel points corresponding to the spatial positions.

And 8, obtaining a synthetic graph.

And respectively adding the denoised image of the invariant region of each time phase image and the pixel value of the space position corresponding to the image of the variant region of the time phase image to obtain a composite image of the time phase image 1 and a composite image of the time phase image 2.

And 9, denoising the sparse dictionary.

9a) And respectively initializing the over-complete sparse representation dictionaries of the two time phase images, wherein the over-complete sparse representation dictionaries are selected as Discrete Cosine Transform (DCT) dictionaries.

9b) Image blocks of a composite image of two time-phase images are obtained according to the following formula:

I_ij=R_ijI

wherein, I_ijImage blocks representing a composite image, i, j representing the reference numerals of the rows and columns of the image blocks, R_ijDenotes an image decimation matrix, and I denotes a composite image.

9c) Judging whether the noise standard deviation of the image of the unchanged area is less than or equal to 5, if so, executing the step 9d) for 5 times; otherwise, step 9d) is performed 10 times.

9d) Adopting a dictionary design method KSVD method of sparse representation, respectively updating an overcomplete sparse representation dictionary and an image block I of a synthesized image relative to a time phase image 1 and a time phase image 2_ijRepresents coefficients.

First, an image block I of the composite image is calculated according to the following formula_ijSparse representation coefficient of (c):

$\begin{array}{cc}{a}_{\mathrm{ij}}=\arg \underset{a_{\mathrm{ij}}}{\min }\left|\right|a_{\mathrm{ij}}|{|}_0& s.t.{\left|\right|D\×a_{\mathrm{ij}}-I_{\mathrm{ij}}\left|\right|}_2^2\≤{\left(1.15\mathrm{\σ}\right)}^2\end{array}$

wherein, a_ijImage block I representing a composite image_ijI, j represents the reference numbers of the rows and columns of the image blocks of the two-time-phase image composite image, argmin represents the variable value operation that makes the objective function take the minimum value, | |, (|₀Expressing the operation of taking the number of non-zero elements, D expressing a sparse representation dictionary, sigma expressing the noise standard deviation of an invariant region of the two-time phase image, and s.t expressing the constraint condition of a formula;

secondly, calculating the residual error of the image block of the synthetic image according to the following formula:

$e_{\mathrm{ij}}^h=I_{\mathrm{ij}}-\underset{t\≠h}{\mathrm{\Σ}}{d}_t{a}_{\mathrm{ij}}\left(t\right)$

wherein,representing residuals of image blocks of the composite image, I, j representing references of rows and columns of image blocks of the two-phase image composite image, each column of the dictionary being called an atom, h representing the position of one atom in all atoms of the dictionary, h =1,2_ijImage blocks representing a composite image, d_tIs the t-th atom of dictionary D, a_ij(t) is the t-th element in the sparse coefficient;

thirdly, obtaining a set consisting of residuals of image blocks of the synthetic image according to the following formula:

$E=\left\{e_{\mathrm{ij}}^h\right|(i,j)\∈c\}$

wherein E represents a set of residual components of image blocks of the composite image,representing the residual of the image block of the composite image, h representing the position of an atom in all atoms of the dictionary, ∈ representing the symbol "included in" in the logical representation, c representing the use of the atom d_hOf the composite imageA set of coordinate positions of the blocks, i, j representing the labels belonging to the rows and columns respectively in the set c of coordinate positions;

fourthly, performing singular value decomposition on the set E consisting of the residual errors according to the following formula:

E=UΔV^T

wherein U represents a left singular matrix, Δ represents a singular value matrix, and V^TRepresenting the transpose of the right singular matrix;

fifth, replace atom d with the first column in left singular matrix U_h，d_hRepresenting the h-th atom of the dictionary D, transposed V by means of the right singular matrix^TIs multiplied by delta (h, h) instead of a_ij(h) And (6) updating.

9e) And respectively estimating the denoising result graph of the time phase image by adopting a reconstruction formula to obtain a denoising result graph of the time phase image 1 and a denoising result of the time phase image 2.

The reconstruction formula is as follows:

$J={[\mathrm{\λX}+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R}_{\mathrm{ij}}^T{R}_{\mathrm{ij}}]}^{-1}[\mathrm{\λI}+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R}_{\mathrm{ij}}^{T}D{a}_{\mathrm{ij}}]$

wherein, J tableThe denoising result of the time-phase image is shown, wherein lambda represents a Lagrange multiplier, lambda = 30/sigma, sigma represents the noise standard deviation of an unchanged area of the time-phase image, X represents a unit matrix with the same size as the synthesized image, i, j represents the reference number of the row and the column of the image block, Ri_jA matrix representing the decimation of the image is shown,denotes Ri_jIs a transformation of [ ·]^-1Representing the inverse operation of the matrix, I representing the composite image, D representing the updated overcomplete sparse representation dictionary, a_ijThe updated sparse representation coefficients are represented.

And step 10, outputting the result.

And outputting the denoising result graphs of the phase image 1 and the phase image 2.

The effect of the present invention is further explained by combining simulation as follows:

1. description and setting of simulation experiment conditions:

the input images used in the simulation experiment of the present invention are, as shown in fig. 2 (a) and 2 (b), 386 × 321 in size, and are SAR images captured by ALOS satellites in different time regions of the same west ampere region formatted with BMP. In simulation experiments, the method and the comparison method are both realized by programming in MATLAB2009a software.

2. And (3) simulation result analysis:

the present invention uses two comparison methods: the Frost method frequently used in SAR image change detection can compare the advantages of the method relative to the common method of change detection; the existing PPB algorithm with better SAR image denoising effect can show that the SAR image denoising effect is good but the change detection effect is not necessarily good, and the advantages of the invention in change detection are verified again; the advantage of the denoising effect of the invention can be seen, and the later-period change detection is also verified, so that the effect is improved.

Calculating the equivalent index of the quantitative evaluation index of the denoising effect according to the following formula:

$\mathrm{ENL}=A\frac{u^2}{{\mathrm{\σ}}^2}$

where u and σ are the mean and standard deviation, respectively, of a homogeneous region in the image after speckle suppression, for intensity image a =1, and for amplitude image a = 4/pi-1. The equivalent index is an index for measuring the intensity of coherent speckles in an image, and the deeper the degree of suppression of coherent speckles, the larger the equivalent index, the more homogeneous regions of 20 × 30 size are selected for both the phase image 1 and the phase image 2, and the region 1 and the region 2 in fig. 2 (a) represent two regions of the phase image 1 selected to calculate the equivalent index, and the region 1 and the region 2 in fig. 2 (b) represent two regions of the phase image 2 selected to calculate the equivalent index.

The equivalent index ENL is used as a quantitative evaluation index of the denoising effect, a Frost filtering method and a PPB filtering method are adopted, and the denoising results of the method are compared as shown in the following table.

Fig. 3 is a denoising result obtained by denoising fig. 2 by using the method of the present invention and the PPB method, fig. 3 (a) shows the denoising result of the PPB method to fig. 2 (a), fig. 3 (b) shows the denoising result of the PPB method to fig. 2 (b), and the PPB filtering method does not filter the speckle noise well and maintains the edge texture.

Fig. 3 (c) shows the denoising result of the method of the present invention for fig. 2 (a), and fig. 3 (d) shows the denoising result of the method of the present invention for fig. 2 (b), which can effectively denoise, and simultaneously retain the texture, contour, edge detail information of the SAR image, especially the change region having a large influence on the later change detection, and suppress the scratch appearing in the homogeneous region.

Fig. 4 is a ratio graph between two time-phase images which are not denoised and denoised, fig. 4 (a) shows the ratio graph between the two time-phase images which are not denoised, fig. 4 (b) shows the ratio graph between the two time-phase images which are denoised by the PPB method, and fig. 4 (c) shows the ratio graph between the two time-phase images which are denoised by the present invention.

Compared with other denoising methods, the method has absolute advantages in ENL.

The experiments show that compared with the existing method, the method has better denoising effect, simultaneously reserves the texture, contour and edge detail information of the SAR image, particularly the change region with large influence on the later change detection, and inhibits the scratches in the homogeneous region. The invention keeps the consistency of the change area better, reduces the pseudo information generated by noise suppression and improves the precision of change detection.

Claims (7)

1. A SAR image sparse denoising method based on change detection comprises the following steps:

(1) inputting an image:

1a) selecting two single-polarized SAR images which are generated by satellites at the same place and different time and have the same registered size;

1b) reading a time phase image 1 in a first single-polarized SAR image matrix;

1c) reading a time phase image 2 in a second single-polarized SAR image matrix;

(2) obtaining a change area image:

2a) taking difference between the two time phase SAR images, and taking an absolute value of the difference operation result to obtain a difference image;

2b) carrying out mean filtering on the difference image by using a template with the window size of 3 multiplied by 3 to obtain a filtered difference image;

2c) calculating according to the following formula to obtain the average energy of any pixel point in the filtered difference image,

$E=\frac{1}{9}\underset{q\∈B}{\Σ}{g}_q^2$

wherein E represents the average energy of any pixel, q represents the pixel number in the neighborhood of 3 × 3 pixels centered on the pixel, q is 1,2, …, 9, B represents the neighborhood of 3 × 3 pixels centered on the pixel, g_qRepresenting the gray value of the q-th pixel in the neighborhood B;

2d) replacing the gray value of the pixel point in each difference image with the average energy of the pixel point in each difference image to obtain an average energy map;

2e) calculating the gray level mean value of the average energy map, setting the mean value as a threshold value, assigning the gray level value which is greater than or equal to the threshold value as 1, and assigning the gray level value which is smaller than the threshold value as 0 to obtain a binary average energy map;

2f) performing morphological erosion operation on the binary average energy map by using 3 x 3 structural elements to obtain a change detection image, wherein a region with a pixel value of 1 in the change detection image represents a change region of the detected two-time phase image;

2g) setting the gray values of pixels outside the space positions corresponding to the change areas of the two time phase images in the time phase image 1 and the time phase image 2 as 0 respectively to obtain an initial change area image of the time phase image 1 and an initial change area image of the time phase image 2;

(3) setting a change class:

3a) subtracting the initial change area image of the time phase image 2 from the initial change area image of the time phase image 1 to obtain an initial change area difference image;

3b) setting the pixel value in the difference image of the initial change area to be less than or equal to the set threshold value T₁The pixel points of (1) are set to be in a non-change class; the pixel value in the difference image of the initial change area is larger than the threshold value T₁The pixel points are set as variation classes;

3c) setting pixel values in an initial change area of a time phase image 1 and an initial change area of a time phase image 2 corresponding to the spatial positions of pixel points set as change classes in the initial change area difference image as 1;

(4) obtaining two classification binary images:

4a) counting the number of pixels which are larger than 0 in a neighborhood of 3 multiplied by 3 and taking the pixel as the center for each pixel which belongs to a non-change type in the difference image of the initial change area;

4b) comparing the number of pixels greater than 0 in the neighborhood with a threshold T₂、T₃The following operations are performed to obtain two classification binary maps of the time phase image 1 and two classification binary maps of the time phase image 2:

if the number of pixels larger than 0 in the neighborhood is larger than or equal to the threshold value T₂Setting the value of the pixel point of the initial change area image of the time phase image 1 corresponding to the spatial position of the non-change pixel satisfying the condition as 1;

if the number of pixels larger than 0 in the neighborhood is less than the threshold value T₃Setting the value of the pixel point of the initial change area image of the time phase image 2 corresponding to the spatial position of the non-change pixel satisfying the condition as 1;

if the number of pixels larger than 0 in the neighborhood is less than the threshold value T₂And is greater than or equal to threshold T₃Setting the values of the pixel points of the initial change area image of the time phase image 1 and the initial change area image of the time phase image 2 corresponding to the spatial position of the pixel of the non-change class meeting the condition as 1;

setting the value of pixel points which are not marked in the initial change area image of the time phase image 1 and the initial change area image of the time phase image 2 as 0;

(5) two-region images are obtained:

5a) multiplying the two classification binary images of the time phase image 1 by pixel points of pixel values corresponding to the spatial position of the image of the time phase image 1 to obtain a change area image of the time phase image 1; subtracting the change area image of the time phase 1 image from the time phase image 1 to obtain an unchanged area image of the time phase image 1;

5b) multiplying the two classified binary images of the time phase image 2 by pixel points of pixel values corresponding to the space position of the time phase image 2 to obtain a change area image of the time phase image 2; subtracting the changed area image of the time phase image 2 from the time phase image 2 to obtain an unchanged area image of the time phase image 2;

(6) denoising the unchanged area:

denoising the invariant region images of the time phase image 1 and the time phase image 2 respectively by using a Wavelet Hidden Markov Tree (WHMT) method to obtain denoised images of the invariant regions of the time phase image 1 and the time phase image 2;

(7) estimating the standard deviation of noise:

7a) the standard deviation of the noise of the invariant region images of the time phase image 1 and the time phase image 2 is obtained according to the following formula

$\mathrm{\σ}=\frac{1}{Z}\underset{Q}{\Σ}(v_1-v_2)-\frac{1}{Z}\underset{Q}{\Σ}{(v_1-v_2)}^2$

Wherein σ represents a noise standard deviation of the invariant region image of the time phase image, Z represents a total number of pixel points in the invariant region image of the time phase image, Q represents the invariant region image of the time phase image, v represents a noise standard deviation of the invariant region image of the time phase image, and₁gray value, v, of pixel points of an invariant region image representing a time phase image₂Representing sum v in denoised image of invariant regions of phase image₁Gray values of pixel points corresponding to the spatial positions;

(8) obtaining a synthetic graph:

respectively adding the denoised image of the invariant region of each time phase image and the pixel value of the space position corresponding to the image of the variant region of the time phase image to obtain a composite image of the time phase image 1 and a composite image of the time phase image 2;

(9) denoising the sparse dictionary:

9a) respectively initializing over-complete sparse representation dictionaries of the two time phase images, wherein the over-complete sparse representation dictionaries are selected as Discrete Cosine Transform (DCT) dictionaries;

9b) image blocks of a composite image of two time-phase images are obtained according to the following formula:

I_ij＝R_ijI

wherein, I_ijImage blocks representing a two-time-phase image composite image, i, j representing the reference numerals of the rows and columns of the image blocks of the two-time-phase image composite image, R_ijRepresenting an image extraction matrix, and I representing a two-time phase image composite image;

9c) judging whether the noise standard deviation of the image of the unchanged area is less than or equal to 5, if so, executing the step 9d) for 5 times; otherwise, executing step 9d)10 times;

9d) adopting a dictionary design method KSVD method of sparse representation, respectively updating an overcomplete sparse representation dictionary and an image block I of a synthesized image relative to a time phase image 1 and a time phase image 2_ijSparse representation coefficients of (a);

9e) respectively estimating a denoising result graph of the time phase image by adopting a reconstruction formula to obtain a denoising result graph of the time phase image 1 and a denoising result graph of the time phase image 2;

(10) and outputting a result:

and outputting the denoising result graphs of the phase image 1 and the phase image 2.

2. The SAR image sparse denoising method based on change detection as claimed in claim 1, characterized in that: step 3b) said threshold T₁The value range of (1) is 16-24.

3. The SAR image sparse denoising method based on change detection as claimed in claim 1, characterized in that: step 4b) the threshold value T₂The value range of (1) is 4-7.

4. The SAR image sparse denoising method based on change detection as claimed in claim 1, characterized in that: step 4b) the threshold value T₃The value range of (A) is more than or equal to T₃≤T₂。

5. The SAR image sparse denoising method based on change detection as claimed in claim 1, characterized in that: the wavelet hidden markov tree WHMT method for the unchanged area images of the phase image 1 and the phase image 2 in the step (6) is carried out according to the following steps:

the method comprises the steps that firstly, 3-layer wavelet transformation is respectively carried out on an unchanged area image of a time phase image, and a noise-containing wavelet coefficient set of the unchanged area image of the time phase image is respectively obtained;

secondly, adopting a tree-shaped Markov chain to fit the conditional distribution of the wavelet coefficients among scales to the noisy wavelet coefficient set obtained after wavelet transformation, establishing a unidirectional transmission hidden Markov tree model, and obtaining a hidden Markov tree model parameter set according to the following formula:

$V=\{u_{k\mathrm{\ξ}},{\mathrm{\σ}}_{k\mathrm{\ξ}}^2,p\left(s_k|s_{\mathrm{\ρ}\left(k\right)}\right),p\left(S_0\right)\}$

wherein V represents a parameter set of the hidden Markov tree model, u_kξMeans of the node, σ, representing the wavelet coefficient at the state ξ, k position_kξWavelet coefficient standard deviation, p(s) for k-position node representing state ξ_k|s_ρ(k)) Probability, s, of k-position node at state ξ when parent node representing k-position node is at state ξ_kState variable, s, representing a node at k position_ρ(k)State variables representing the parent of k-position nodes, [ rho (k) represents the parent of k-position nodes, [ p(s) ]₀) Representing the initial state distribution, s, of the root node₀A state variable representing a root node;

thirdly, training the noisy wavelet coefficient by using an expectation-maximization algorithm to obtain an optimal estimated hidden Markov tree model parameter set;

fourthly, estimating the noiseless wavelet coefficient of the obtained optimal estimated hidden Markov tree model parameter set according to the following formula:

$y_k=\underset{\mathrm{\ξ}}{\mathrm{\Σ}}p(s_k=\mathrm{\ξ}|W,V){\mathrm{\σ}}_{k\mathrm{\ξ}}^2/({\mathrm{\σ}}_{k\mathrm{\ξ}}^2+{\mathrm{\σ}}_1^2)$

wherein, y_kRepresenting the estimated noise-free wavelet coefficients, ξ representing the node states, p(s)_kξ | W, V) represents the state variable s of the k-position node with the knowledge of the optimal estimated hidden markov tree model parameter set_kProbability of being state ξ, s_kState variables representing k-position sub-nodes, W represents a noisy set of wavelet coefficients, V represents a parameter set of an optimally estimated hidden Markov tree model, σ_kξStandard deviation, σ, of node representing state ξ, k position₁Representing the standard deviation contained in the wavelet coefficients;

a fifth step of estimating the wavelet coefficients y_kAnd respectively carrying out wavelet inverse transformation to obtain the denoised image of the invariant region image of the time phase image.

6. The SAR image sparse denoising method based on change detection as claimed in claim 1, characterized in that: the method for designing the dictionary by using the sparse representation KSVD in the step 9d) is carried out according to the following steps:

first, an image block I of the composite image is calculated according to the following formula_ijSparse representation coefficient of (c):

$\begin{array}{cc}{a}_{ij}=\arg \underset{a_{ij}}{min}\left|\right|a_{ij}|{|}_0& s.t.\left|\right|D\×a_{ij}-I_{ij}|{|}_2^2\≤{\left(1.15\mathrm{\σ}\right)}^2\end{array}$

wherein, a_ijImage block I representing a composite image_ijI, j represents the reference numbers of the rows and columns of the image blocks of the two-time-phase image composite image, argmin represents the variable value operation that makes the objective function take the minimum value, | |, (|₀Expressing the operation of taking the number of non-zero elements, D expressing a sparse representation dictionary, sigma expressing the noise standard deviation of an invariant region of the two-time phase image, and s.t expressing the constraint condition of a formula;

secondly, calculating the residual error of the image block of the synthetic image according to the following formula:

$e_{ij}^h=I_{ij}-\underset{t\≠h}{\Σ}{d}_t{a}_{ij}\left(t\right)$

wherein,representing the residuals of the image blocks of the composite image, I, j representing the indices of the rows and columns of the image blocks of the two-phase image composite image, each column of the dictionary being called an atom, h representing the position of an atom in all atoms of the dictionary, h 1,2_ijImage blocks representing a composite image, d_tIs the t-th atom of dictionary D, a_ij(t) is the t-th element in the sparse coefficient;

thirdly, obtaining a set consisting of residuals of image blocks of the synthetic image according to the following formula:

$E=\left\{e_{ij}^h\right|\left(i,j\right)\∈c\}$

wherein E represents a set of residual components of image blocks of the composite image,representing the residual of the image block of the composite image, h representing the position of an atom in all atoms of the dictionary, ∈ representing the symbol "included in" in the logical representation, c representing the use of the atom d_hI, j respectively represent the labels of the rows and columns belonging to the set c of coordinate positions;

fourthly, performing singular value decomposition on the set E consisting of the residual errors according to the following formula:

E＝UΔV^T

wherein U represents a left singular matrix, Δ represents a singular value matrix, and V^TRepresenting the transpose of the right singular matrix;

fifth, replace atom d with the first column in left singular matrix U_h，d_hRepresenting the h-th atom of the dictionary D, transposed V by means of the right singular matrix^TIs multiplied by delta (h, h) instead of a_ij(h) And (6) updating.

7. The SAR image sparse denoising method based on change detection as claimed in claim 1, characterized in that: the reconstruction formula of step 9e) is as follows:

$J={\[\mathrm{\λ}X+\underset{ij}{\Σ}{R}_{ij}^T{R}_{ij}\]}^{-1}\[\mathrm{\λ}I+\underset{ij}{\Σ}{R}_{ij}^{T}D\×a_{ij}\]$

where J denotes a denoising result of the phase image, λ denotes a lagrange multiplier, λ is 30/σ, σ denotes a noise standard deviation of an invariant region of the phase image, X denotes an identity matrix having the same size as the synthesized image, i, J denotes a reference numeral of a row and a reference numeral of a column of the image block, and R denotes a reference numeral of a column of the image block_ijA matrix representing the decimation of the image is shown,represents R_ijIs a transformation of [ ·]^-1Representing the inverse operation of the matrix, I representing the composite image, D representing the updated overcomplete sparse representation dictionary, a_ijThe updated sparse representation coefficients are represented.