A SAR image filtering method based on non-subsampled shearlet transform

Abstract

The application discloses a SAR image filtering method based on non-subsampled shearlet transform, and belongs to the technical field of image processing. The application comprises the following steps: performing logarithmic transform on a SAR image, and then performing non-subsampled shearlet (NSST) transform; modifying NSST transform coefficients of each high-frequency subband image based on a set threshold function; performing median filtering on each NSST transform coefficient of a low-frequency subband image; performing exponential transform on the processed low-frequency subband image and high-frequency subband image after NSST inverse transform, and obtaining a filtered SAR image. The threshold function adopted by the application is based on a hard threshold and a soft threshold, and the threshold of different layers of high-frequency subband coefficients is calculated respectively, considering that different subbands contain different noise components. The application can effectively filter out coherent speckle noise in the SAR image, and ensure the reservation of the detail information of the SAR image.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of image processing, and particularly relates to a synthetic aperture radar (SAR) image filtering method based on nonsubsampled shearlet transform (NSST). BACKGROUND

[0002] Synthetic aperture radar (SAR) has all-weather and all-day observation capability on the ground. The SAR imaging mechanism is different from the optical image. The electromagnetic wave emitted by SAR is scattered by randomly distributed scatterers on the ground, resulting in different travel distances of the echoes through different scatterers, and the speckle noise is generated due to coherent superposition. The coherent speckle noise is formed by the reflection waves of numerous scatterers in a resolution unit. When the image pixel spacing is comparable to the radar resolution, the noise power is non-correlated. In this case, it can be assumed that the coherent speckle is a kind of incoherent multiplicative noise. The coherent speckle noise seriously affects the signal-to-noise ratio of the SAR image, and causes difficulties for subsequent target recognition and feature extraction. In order to be able to effectively extract information from the SAR image, the coherent speckle noise filtering work is usually carried out before the SAR image processing.

[0003] The coherent speckle suppression technology for SAR image is generally divided into two categories: multi-view smoothing processing technology before imaging and filtering technology after imaging. The filtering technology after imaging can be divided into spatial domain filtering technology, diffusion domain filtering technology, transform domain filtering technology and the like. The spatial domain filtering algorithm is to obtain the local statistical characteristics of the image through a sliding window and then realize filtering processing. The algorithm is easy to implement and has good real-time performance. At present, the Lee filtering algorithm is most widely used. However, the method is prone to produce scallop effect and false fine line phenomenon. The diffusion domain filtering algorithm uses local structure information to constrain the diffusion intensity and diffusion direction to realize the purpose of denoising. The algorithm has good effects in noise suppression and edge preservation. However, the method is prone to image detail blurring and insufficient noise suppression phenomenon. The transform domain processing mainly includes Contourlet transform, non-subsampled Contourlet transform, shearlet transform, NSST transform and the like. In early stage, scholars use the characteristics of discrete wavelet transform (DWT) such as multi-resolution and decorrelation to make the image signal realize good signal-noise separation in the wavelet domain. In 1995, scholars of Stanford University such as Donoho and Johnstone proposed a wavelet threshold denoising method. A suitable threshold is selected. The wavelet coefficient greater than the threshold should be retained because it is considered to be generated by the signal. The coefficient less than the threshold is set to zero. The noise is removed by using this criterion. The method is most studied at present. The two key places are the selection of the threshold and the threshold function. However, the threshold function based on hard threshold denoising can better maintain the image edge characteristics. However, the threshold function has discontinuity, which leads to ringing phenomenon. The threshold function based on soft threshold denoising can completely suppress the noise concentrated coefficient. The threshold function has continuity. However, high-frequency information is lost, which leads to the loss of image edge and texture information. Therefore, a suitable threshold function needs to be designed to achieve good denoising effect. SUMMARY

[0004] The present application is directed to the coherent speckle noise in the SAR image. A SAR image filtering method based on non-subsampled shearlet transform is proposed to ensure that the coherent speckle noise in the SAR image is effectively filtered out while the detail information of the SAR image is preserved.

[0005] To achieve the above object, the technical scheme adopted by the present application is as follows:

[0006] A SAR image filtering method based on non-subsampled shearlet transform, the method comprising the following steps:

[0007] Step 1, logarithmic transformation is performed on the SAR image to make the multiplicative noise (coherent speckle noise) in the SAR image become additive noise;

[0008] Step 2, performing a non-subsampled shearlet (NSST) transform on the logarithmically transformed SAR image to obtain a low-frequency subband coefficient matrix and a plurality of high-frequency subband coefficient matrices;

[0009] Step 3, modifying the transform coefficients of each high-frequency subband coefficient matrix:

[0010] For each high-frequency subband coefficient matrix, estimating the noise variance of each image based on a median estimation method where k represents a high-frequency subband coefficient matrix identifier;

[0011] For each high-frequency subband coefficient matrix, estimating the noise variance of the current high-frequency subband coefficient based on the square average of different subband coefficients

[0012] According to the formula obtain the threshold T of the kth high-frequency subband coefficient matrix k ;

[0013] According to the threshold function, modify the NSST transform coefficients of each high-frequency subband coefficient matrix;

[0014]

[0015] wherein, represents an element in the kth high-frequency subband coefficient matrix, and the subscripts i,j represent the matrix coordinates, represents the modified NSST transform coefficient, and e represents the natural base;

[0016] Step 4, performing median filtering on each NSST transform coefficient of the low-frequency subband coefficient matrix;

[0017] Step 5, performing NSST inverse transform on the low-frequency subband coefficient matrix and the high-frequency subband coefficient matrix processed in steps 3 and 4;

[0018] Step 6, performing exponential transform on the NSST inverse transformed image to obtain the filtered SAR image.

[0019] In summary, due to the adoption of the above technical solutions, the present application has the following advantages:

[0020] In the present application, a denoising method based on NSST transform is realized. In view of the characteristic that the SAR image speckle noise is multiplicative noise, the noise model is converted through logarithmic transformation, and the NSST transform direction sensitivity, sparsity and other characteristics are used to suppress noise in the transform domain through threshold estimation and threshold function design. In the present application, the threshold function adopted is based on the hard threshold and the soft threshold, and considering that different subbands contain different noise components, the threshold is calculated for different layers of high-frequency subband coefficients. BRIEF DESCRIPTION OF DRAWINGS

[0021] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed to be used in the embodiments will be briefly introduced. Obviously, the drawings in the following description only constitute some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor based on these drawings.

[0022] Figure 1 A processing flow chart of a SAR image filtering method based on non-subsampled shear wave transform provided by the embodiment of the present application.

[0023] Figure 2 A non-subsampled shear wave transform decomposition schematic diagram;

[0024] Figure 3 In the embodiment of the present application, a real SAR image is used, which includes two images, SAR1 and SAR2.

[0025] Figure 4 A comparison chart before and after denoising of SAR1, wherein (4-a) is a SAR1 noise chart, and (4-b) is a denoised image by the method of the embodiment of the present application.

[0026] Figure 5 A comparison chart before and after denoising of SAR1, wherein (5-a) is a SAR1 noise chart, and (5-b) is a denoised chart by combining a hard threshold function with NSST transform.

[0027] Figure 6 A comparison chart before and after denoising of SAR1, wherein (6-a) is a SAR1 noise chart, and (6-b) is a denoised chart by combining a soft threshold function with NNSST transform.

[0028] Figure 7 A comparison chart before and after denoising of SAR2, wherein (7-a) is a SAR2 noise chart, and (7-b) is a denoised image by the method of the embodiment of the present application.

[0029] Figure 8 A comparison chart before and after denoising of SAR2, wherein (8-a) is a SAR2 noise chart, and (8-b) is a denoised chart by combining a hard threshold function with NSST transform.

[0030] Figure 9 A comparison chart before and after denoising of SAR2, wherein (9-a) is a SAR2 noise chart, and (9-b) is a denoised chart by combining a soft threshold function with NNSST transform. DETAILED DESCRIPTION

[0031] In order to make the objects, technical solutions and advantages of the present application clearer, the following will further describe the embodiments of the present application in detail with reference to the drawings.

[0032] As a possible implementation, as shown in Figure 1 the embodiment of the present application provides a SAR image filtering method based on non-subsampled shearlet transform, which specifically comprises the following steps.

[0033] Step S1, logarithmic transformation is performed on the SAR image to convert the multiplicative noise contained in the SAR image into additive noise, and the following formula is used to perform logarithmic transformation to obtain the transformed image G:

[0034] G(x, y) = ln(1 + I(x, y))

[0035] Where (x, y) represents the pixel point coordinates, G(x, y) represents the pixel value of the image G, and I(x, y) represents the pixel value of the image I (SAR image).

[0036] Step S2, NSST transform processing is performed on the image G, and the NSST transform processing uses a local shearlet filter to realize reverse localization, which can omit the step of downsampling processing. The NSST transform completes sampling in the decomposition process, and then uses a local shearlet filter to filter the multi-scale coefficient matrix to realize directional localization.

[0037] Where the decomposition process of the NSST transform is as shown in Figure 2 The implementation process mainly consists of the following two parts:

[0038] 1) Multi-scale partitioning.

[0039] The non-subsampled pyramid (NSP) decomposition uses a two-channel non-subsampled filter bank to make the NSST have multi-scale property. After one layer of NSP decomposition of the image G, the low-frequency coefficient and the high-frequency coefficient of the image are obtained. Then, the NSP decomposition of each layer is iterated on the low-frequency component obtained by the decomposition of the upper layer to obtain the singular points of the image. Since there is no downsampling in the NSST process, after m (a preset value) layers of NSP decomposition of the image G, m+1 sub-band images (i.e. sub-band coefficient matrices) with the same size as G can be finally obtained. The m+1 images include 1 low-frequency component and m high-frequency components.

[0040] 2) Directional localization.

[0041] NSST utilizes local shearlet filter to convolve with high frequency subband components to achieve directional localization. NSST transform reflects the local shearlet filter from the pseudo-polar grid coefficients to the Cartesian coordinate system, and calculates its two-dimensional inverse discrete Fourier transform to obtain the non-subsampled shearlet transform coefficients Y.

[0042] Step S3, estimate the high frequency subband noise variance, and calculate the threshold value.

[0043] After NSST transform of the SAR image, the noise information will be distributed in the high frequency subband, and the noise is mainly contained in the small coefficients after the transform. The amount of noise contained in different subbands is different, so the relationship between the noise energy and the signal energy of each subband is estimated respectively. The noise is randomly distributed on the subband coefficients, and the coefficient value is small. By setting the threshold value, the transform coefficients are compared with the threshold value, and the transform coefficients are modified according to the proposed threshold function. For the threshold value calculation in the present application, the MapShrink threshold calculation method is adopted. After NSST transform of the image, the image information X is estimated through the shearlet transform coefficients Y, and according to the Bayesian maximum posterior:

[0044]

[0045] According to the definition of conditional probability defined by Bayes:

[0046]

[0047] Assuming that the noise obeys Rayleigh distribution, the estimated conditional probability can be expressed as:

[0048]

[0049] In the formula, f(X) = log(p X (X)). If p X (X) obeys Gaussian distribution with mean 0 and variance σ 2 , then:

[0050]

[0051] If p X (X) obeys Gaussian distribution with variance σ 2 , then:

[0052]

[0053] The above formula can be regarded as a soft threshold function, and the threshold value is:

[0054]

[0055] In the above formula, σ represents the noise variance, and σ represents the noise variance of the subband coefficient. In the above formula, σ represents the noise variance, and σ represents the noise variance of the subband coefficient.

[0056] The noise variance is very important for describing the noise statistical characteristics, so it is necessary to estimate the noise variance as accurately as possible during denoising. In the embodiment of the present application, The specific calculation method of σ is as follows:

[0057]

[0058] In the formula, y(i,j) represents the high-frequency subband coefficient, and Median represents the method of using median estimation, which can eliminate the influence of useful signals in the subband on the noise variance estimation.

[0059] The noise variance σ of the subband coefficient can be estimated according to the subband coefficient:

[0060]

[0061]

[0062] In the formula, N is the number of subband coefficients, and y(m) represents the subband coefficient.

[0063] In step S4, the NSST denoising is performed by using a threshold method to process the shear wave transform coefficient. Considering that the hard threshold function has good edge performance when processing signals, but due to the discontinuity of the function itself, unnecessary oscillation may occur in the reconstructed signal; the soft threshold function has good continuity, and the denoising result is smoother, but when the wavelet coefficient is large, there is a fixed deviation between the estimated wavelet coefficient and the original coefficient, which causes the high-frequency part of the signal to be lost, resulting in poor approximation of the reconstructed signal compared with the original signal, and distortion is easy to occur. The threshold function used in the embodiment of the present application is as follows:

[0064]

[0065] In the formula, y ij represents the subband coefficient, that is, the element in the kth high-frequency subband coefficient matrix, T k represents the threshold value calculated from the kth high-frequency subband coefficient, represents the subband coefficient after threshold processing.

[0066] It can be found from the above formula that, It can be seen that the function values on the left and right sides of the above formula at y ij = T k are equal, which indicates that the function is continuous at the threshold T k . In addition, when y ij tends to infinity,

[0067]

[0068] From the above formula, when y ij tends to infinity, That is, the modified shear wave transform coefficient is consistent with the true transform coefficient, thereby reducing the deviation between the denoised image and the true noise-free image.

[0069] Step S5, considering that the low-frequency sub-band may contain a small amount of noise components, the low-frequency sub-band coefficient is subjected to median filtering processing. Considering that the low-frequency sub-band component after shear wave transform mainly reflects the contour information of the image, the high-frequency sub-band component mainly contains the edge, texture and other detail information of the image, and the low-frequency component contains less information, so the filtering is not considered to cause loss of detail information. In the present application, the filtering processing of the low-frequency sub-band component adopts a 3x3 median filter.

[0070] Step S6, the NSST transform coefficient after threshold processing is subjected to inverse transform to obtain a denoised logarithmic image G_L.

[0071] Step S7, the image G_L is subjected to exponential transform G_F = e G_L -1, and after the image G_F is subjected to normalization processing, it is multiplied by the highest gray level (usually 255), and finally the denoised image is obtained.

[0072] The effect of the SAR image filtering method based on non-subsampled shear wave transform provided by the embodiment of the present application can be further illustrated by the following simulation experiment:

[0073] 1) Experimental data, in the present embodiment, two real SAR images with image sizes of 276x276 and 277x277 are subjected to filtering experiment, as shown in the following table. Figure 3

[0074] 2) Experimental content and result analysis.

[0075] In order to verify the effectiveness of the method of the embodiment of the present application, three objective evaluations are adopted to measure the effect of the SAR image filtering algorithm: mean square error, peak signal-to-noise ratio and equivalent number.

[0076] Mean square error (MSE): MSE represents the difference between the denoised SAR image and the ideal SAR image, and the smaller the MSE, the closer the denoised SAR image is to the ideal SAR image, in other words, the better the noise suppression effect. The calculation formula is shown in the following table.

[0077]

[0078] In the formula, f i ′ ​f represents the image after denoising i N represents an ideal image, and the size of the image.

[0079] Peak signal-to-noise ratio (PSNR): PSNR represents the ratio of the maximum power of the SAR image to the noise power. The larger the PSNR, the smaller the proportion of noise, in other words, the better the noise suppression effect. The calculation formula is as follows.

[0080]

[0081] Equivalent number of looks (ENL): The calculation of the equivalent number of looks is irrelevant to the original noise-free ideal image. For a uniform area, the larger the ENL value, the better the denoising effect. In general, several uniform areas in the denoised image are selected, and the average value of the ENL of the uniform areas is calculated as an evaluation index.

[0082]

[0083] In the above formula, μ and σ 2 respectively represent the mean and variance of a uniform area in the denoised SAR image.

[0084] 3) SAR image filtering experiment.

[0085] Considering that the peak signal-to-noise ratio and the mean square error require a noise-free ideal image, the SAR image to which multiplicative noise is added is filtered by using the method of the embodiment of the application. Because the ENL value needs to be calculated for a uniform area of the image, the ENL values of the white box part in Figure 3 are calculated in this embodiment.

[0086] The method of the embodiment of the application is compared with the NSST denoising algorithm combined with a hard threshold and the denoising algorithm combined with a soft threshold function. Two real SAR images in Figure 3 to which multiplicative Rayleigh noise with a variance of 0.10 is added are denoised, and the denoising result images of the methods are as shown in Figure 4 to Figure 9 The changes of the three evaluation indexes before and after denoising are compared, and the results are as shown in the following table.

[0087] Table 1 Comparison of indexes before and after denoising of SAR1

[0088] PSNR (dB) ENL MSE Before denoising by the method of the present embodiment 20.71 5.53 551.99 After denoising by the method of the present embodiment 25.18 17.81 197.07 Before denoising by NSST+hard threshold function 20.67 5.60 557.64 After denoising by NSST+hard threshold function 21.56 6.36 453.70 Before denoising by NSST+soft threshold function 20.75 5.64 546.73 After denoising by NSST+soft threshold function 24.49 11.96 231.16

[0089] Table 2 Comparison of indexes before and after denoising of SAR2

[0090] PSNR (dB) ENL MSE Before denoising by the method of the present embodiment 19.79 8.69 682.77 After denoising by the method of the present embodiment 24.26 100.37 243.84 Before denoising by NSST+hard threshold function 19.76 8.70 686.51 After denoising by NSST+hard threshold function 20.84 10.70 535.40 Before denoising by NSST+soft threshold function 19.80 8.63 681.15 After denoising by NSST+soft threshold function 23.04 31.84 323.16

[0091] From Table 1, the peak signal-to-noise ratio of the image after denoising by the method of the embodiment of the application is improved by 4.47, the equivalent visual number is increased by 12.28, and the mean square error is reduced by 354.92; the peak signal-to-noise ratio of the image after denoising by NSST combined with a hard threshold function algorithm is improved by 0.89, the equivalent visual number is increased by 0.76, and the mean square error is reduced by 103.94; the peak signal-to-noise ratio of the image after denoising by NSST combined with a soft threshold function algorithm is improved by 3.74, the equivalent visual number is increased by 6.32, and the mean square error is reduced by 315.57. From Table 2, the peak signal-to-noise ratio of the image after denoising by the method of the embodiment of the application is improved by 4.47, the equivalent visual number is increased by 91.68, and the mean square error is reduced by 438.93; the peak signal-to-noise ratio of the image after denoising by NSST combined with a hard threshold function algorithm is improved by 1.08, the equivalent visual number is increased by 2, and the mean square error is reduced by 151.11; the peak signal-to-noise ratio of the image after denoising by NSST combined with a soft threshold function algorithm is improved by 3.24, the equivalent visual number is increased by 23.21, and the mean square error is reduced by 357.99.

[0092] In summary, the method of the application has good denoising effect, and the method of the application can exhibit good denoising performance in the peak signal-to-noise ratio, the equivalent visual number and the mean square error.

[0093] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the application, but not to limit them; although the application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for some technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the application.

[0094] The above only describes some embodiments of the application. Those skilled in the art can make several modifications and improvements without departing from the inventive concept, and these all belong to the protection scope of the application.

Claims

1. A non-subsampled shearlet transform based SAR image filtering method, characterized in that, The method comprises the following steps: Step 1, logarithmic transformation is performed on the SAR image I so that the multiplicative noise in the SAR image becomes additive noise, and the logarithmic transformation is performed by using the following formula to obtain a transformed image G: G(x, y) = ln(1 + I(x, y)) Where (x, y) represents the pixel point coordinates, G(x, y) represents the pixel value of the image G, and I(x, y) represents the pixel value of the image I; Step 2, a non-subsampled shearlet (NSST) transform is performed on the SAR image after logarithmic transformation to obtain a low-frequency sub-band coefficient matrix and a plurality of high-frequency sub-band coefficient matrices; Wherein the decomposition process of the NSST transform comprises: 1) Multi-scale decomposition: the non-subsampled pyramid (NSP) decomposition uses a two-channel non-subsampled filter bank to make the NSST multi-scale, and the image G is decomposed by one layer of NSP to obtain the low-frequency coefficient and the high-frequency coefficient of the image. Then, the NSP decomposition of each layer is iterated on the low-frequency component obtained by the decomposition of the upper layer to obtain the singular points of the image. After m layers of NSP decomposition of the image G, m+1 sub-band images with the same size as G are finally obtained, including 1 low-frequency component and m high-frequency components; 2) Direction localization: the NSST utilizes a local shearlet filter to perform convolution operation with a high-frequency sub-band component to achieve direction localization, the NSST transform reflects the local shearlet filter from the coefficient in the pseudo-polar grid coordinate system to the Cartesian coordinate system, and calculates the two-dimensional inverse discrete Fourier transform to obtain a non-subsampled shearlet transform coefficient Y; Step 3, the transform coefficients of each high-frequency sub-band coefficient matrix are modified: For each high frequency subband coefficient matrix, estimate the noise variance of the current high frequency subband coefficient matrix based on a median estimation method where k denotes a high frequency subband coefficient matrix distinguisher; denotes an element in the kth high frequency subband coefficient matrix, the subscript i,j denotes a matrix coordinate, and Median {} denotes a median estimation method. For each high frequency subband subband coefficient matrix, estimate the noise variance of the current high frequency subband coefficient based on the square average of different subband coefficients According to the formula The threshold T of the kth high-frequency subband coefficient matrix is obtained k ; The NSST transform coefficients of each high-frequency sub-band coefficient matrix are modified according to a threshold function; wherein denotes an element in the k-th high-frequency subband coefficient matrix, the subscripts i,j denote the matrix coordinates, denotes the modified NSST transform coefficients, e denotes the natural base; Step 4, each NSST transform coefficient of the low-frequency sub-band coefficient matrix is median filtered; Step 5, the low-frequency sub-band coefficient matrix and the high-frequency sub-band coefficient matrix processed in steps 3 and 4 are subjected to NSST inverse transformation; Step 6, exponential transformation is performed on the image after NSST inverse transformation to obtain a filtered SAR image.

2. The method of claim 1, wherein, In step 4, a 3*3 median filter is used to perform median filtering on the low-frequency sub-band coefficients.

Images