Remote sensing images denoising based on multi-linear weighted nuclear norm minimization

Xiang-yang KONG, Bao-gen XU, Jie ZHOU

(1Ministry of Basic Education, Sichuan Engineering Technical College, Deyang 618000, China)(2School of Automation, Northwestern Polytechnical University, Xi’An 710072, China)(3School of Science, East China Jiaotong University, Nanchang 330013, China)

Abstract: Classical vector-based or matrix-based denoising methods for remote sensing images (RSIs) may cause distortion both in spatial domain and spectral domain. To improve denoising performance and reduce the distortion, a denoising method based on multi-linear weighted nuclear norm is proposed. Firstly, by considering spectral continuity and inter-dependency of three unfolding matrices along three modes of RSI, a multi-linear rank is constructed to model the spatial and spectral nonlocal similarity. Then, to make the proposed method more tractable, a variable splitting-based technique is used to solve the optimization problem. Experiment results reveal that the proposed method can provide substantial improvements over the current state-of-the-art methods in terms of both objective metric and subjective visual quality.

Key words: Denoising, Weighted nuclear norm, Alternating direction method of minimization, Peak signal-to-noise ratio, Structural similarity index measurement

Remote sensing images (RSIs) carry rich spectral and spatial informations, and they have been wildly used in many fields [1-5]. For the non-ideal optical and electronic device, acquired RSI is always corrupted with noise, the existing of noise can not only badly degrades the visual appearance, but also has bad influence on the subsequent applications, including unmixing[1-2], super-resolution[3], classification[4], compressive sensing reconstruction [5] and so on. Therefore, denoising is a critical pre-processing step for the applications of RSIs.

The noisy RSIs can be denoised by two typical ways in general. One way is that the same scene can be captured several times and then the acquired data is averaged. However, this procession costs high imaging time and always is limited by physical conditions. The post-processing methods are the other most common way. Please refer to reference [6] to survey these denoising methods. Here, the patch-based RSI denoising methods, which are most related to this work, are reviewed.

The patch-based methods, which exploit the non-local similarity, have shown popularity and effectiveness over denoising. In reference [7], Dabov K et al. first proposed the Block-Matching and 3D filtering (BM3D) method, by grouping similar blocks into 3D data arrays and then shrinking/filtering in the 3D transform domain, the patches were then aggregated at each location to obtain the denoised image. Dong W et al. [8] introduced nonlocally centralized sparse representation (NCSR) method to make the local sparsity and the nonlocal sparsity constraints into a variational framework. With the low rank prior and the non-local similarity, the nuclear norm minimum (NNM) denoising algorithm shows good performance in dealing with 2-D images. To pursue the convexity, NNM treats every singular value equally, so the denoising result was not ideal. In order to improve the flexibility of the nuclear norm, a weighted nuclear norm minimum (WNNM) was proposed in [9]:

(1)

wherew=[w1,w2,…,wn],wirepresented non-negative weight ofσi(X),σi(X) was thei-th largest singular value ofX.

The RSI can be molded as a tensor of order 3, with one-dimensional spectral information and two-dimensional spatial information [10]. The methods BM3D and WNNM usually regard each spectral band or the unfolding matrix of RSI as a gray-level image, the difference is that the former method processes RSI bandwise, while the later method processes RSI globally. The NCSR makes the local sparsity and the nonlocal sparsity constraints into a variational framework, but it does not take spectral correlations and the inter-dependency of three modes into consideration, and often resulted in a less competitive performance.

Motivated by above methods, a method based on multi-linear weighted nuclear norm is proposed, which can fully exploit the inter-dependency of three unfolding modes and non-local similarity inside each mode. Firstly, the low-rankness of RSI’s unfolding matrices along 3 modes is considered. Secondly, considering the structural similarity inside each mode, a multi-linear weighted nuclear norm is proposed to take advantage of the nonlocal spatial and spectral information jointly. Thirdly, an optimization algorithm is developed to solve the proposed mathematical model. The comparison experiments of RSI denoising which are carried out both on simulated and real datasets show that the proposed method substantially outperforms the current state-of-the-art RSI denoising methods.

1 Notation of tensor algebra

As the extension of matrix and vector, a tensor of orderNis regard as an N-dimensional array, it can be represented asA∈RI1×I2×…×IN. The elementai1,i2,…,iNinAis a scalar, where 1≤in≤In,1≤n≤N. The matrix, vector and scalar are called 2-order, 1-order and 0-order tensor, respectively. Each dimension of tensor is called a mode. Then-th dimension of a tensor is also called mode-n.

The fiber (or mode-nvector) of anN-th order tensorAis aIndimensional vector, which is extracted fromAby varying indexinwhile keeping the other indices fixed. By rearranging all the mode-ifibers ofAto be the columns of a matrix, the tensor can be converted to a matrixA(i)∈RIi×(I1I2…Ii-1Ii+1…IN). This is the process of matrization (or unfolding) of tensor.A(n)is also represented byunfoldn(A). For other tensor algebra operations, please refer to [11].

2 Model proposal and optimization

From a mathematical point of view, there are two kinds of noise to result in the degraded image, one is additive noise, and the other is multiplicative noise. The multiplicative noise can be transformed to additive noise by logarithmic transformation. So, the noise considered here is additive noise. Suppose the observed degraded RSI is denoted asY, the aim of denoising is to find a good estimate of original clear RSIXfromY. The degraded process of the RSI can be modeled as：

Y=X+N

(2)

WhereNis the additive white gaussian noise, and the size ofY,X,NisI1×I2×I3, in whichI1andI2stands for the numbers of row and column in spatial domain, andI3is the dimension of spectral domain.

For tensorA, its multi-linear rank rankmulti(A) is defined as follows[13]:

(3)

The clean RSI’s unfolding matrix along mode-3 can be regarded as low-rankness [12], but it is observed that the other two matrices along mode-1 and mode-2 have the same property. Considering the interdependency of different modes, the RSI can be denoised by 3 modes unfolding matrices separately. Besides, there are similar patches in each matrix, the yellow rectangles in Fig.1 are examples of similar patches. Considering of these facts, the denoising problem can be modeled as weighted multi-linear rank:

Fig.1 Examples of similar patches in each mode, as shown in the yellow rectangles (From top to bottom are matrices unfolding along 1, 2 and 3 mode, respectively)

(4)

The multi-linear rank defined by Eq. (4) can be regarded as a generalization of nuclear norm of the matrix. The denoising model in [13] penalizes the singular values of each unfolding matrix equal, it means that the same soft-threshold will be applied to all singular values. This is not reasonable since different singular values may have different importance and hence they should be treated differently. Then the multi-linear weighted nuclear norm can be used as relaxation of multi-linear weighted rank:

(5)

The problem (4) is an ill-posed inverse problem, based on the definition of multi-linear weighted nuclear norm in Eq. (5), Eq.(4) can be rewritten as:

(6)

To solve the proposed denoising model in Eq. (6), the variables are separated and assign a separable variable to each unfolding mode ofX. That isUi=XorUi,(i)=X(i)(i=1,2,3). Thus Eq. (6) can be rewritten as:

(7)

To obtain different unfoldings of tensor, the problem of Eq. (7) can be divided into three individual problems. For each problem, by initializingY=X, andXi=X, the augmented Lagrange function of Eq. (7) is obtained by introducing the Lagrange multiplierΛi(i=1,2,3),which is formulated as follows：

(8)

whereβ>0 is a penalty parameter.

The alternating direction method of minimization (ADMM) used to solve Eq. (8). It is divided into the following subproblems:

Subproblem 1: updateUiwhile fix other variables

(9)

The problem of Eq. (9) is a weighted nuclear norm minimization problem, while the weighted nuclear norm problem is generally non-convex, and its convexity depends on the order of its weights. For Eq. (9), the following theorem is established [9]:

Theorem 1: For arbitraryτ>0,X,Y∈Rm×n, andws≥ws-1≥…≥w1≥0,(s=min(m,n)), the optimal solution of the following problem

(10)

isX=Dτw(Y)=P∑τwQT, whereY=P∑QT, ∑τw=diag{(σi(Y)-τw)+}.

Therefore, the closed-form solution of Eq. (9) is

(11)

Here fold (·) is the inverse operation of unfolding, which means to reshape a matrix to a tensor.

Subproblem 2: update Lagrange multipliers

Λi=Λi+β(Λi-Xi)

(12)

Subproblem 3: updateXi

(13)

3 Experiment and results

To compare the performance of denoising methods on both synthetic and real RSIs datasets, the RSIs datasets of Washington DC Mall (WDC for short, including 191 bands), Indian Pines (including 220 bands) and Urban (including 162 bands) are used. In synthetic dataset experiment, the quantitative indexes are the peak signal-to-noise ratio (PSNR) [14] and structural similarity index measurement (SSIM) [15]. The unit of PSNR is dB, while the SSIM is a no unit index. The higher the PSNR and SSIM are, the closer the restored image to the original is.

The state-of-the-art RSI denoising methods are used to compare with the proposed method, including BM3D [7] and NCSR [8], WNNM [9]. To keep consistence with the reference papers, the parameters are used as described in these papers. All the experiments are implemented on the Windows7, the Core i5-7300HQCPU@2.5GHz and the 8G memory platform by MATLAB R2014a. In the following experiment, the parameters are set toγ1=0.6，γ2=γ3=0.2.

3.1 Denoising on synthetic datasets

To quantitatively evaluate the effectiveness of the proposed algorithm, the above three methods under varying levels’ Gaussian noise are compared on synthetic datasets.

In Fig.2, the values of PSNR and SSIM on WDC and Urban are plotted under different noisy levels. The proposed method significantly outperforms other three methods in terms of PSNR and SSIM values.

Fig.2 PSNR/SSIM values on WDC ((a) and (b)) and Urban((c)and(d)) datasets with noise variances varying from 5 to 50

To have a clearer view at the exact values of PSNR and SSIM, the PSNR and SSIM values of WDC and Urban under different noise levels are listed in Table 1 and Table 2. It can be seen from these two tables that the proposed method significantly outperforms other three methods under different noisy levels. Compared with WNNM and BM3D, the PSNR values of proposed method are improved from 0.9 dB to 2.6 dB, from 0.7 dB to 2.9 dB and from 3.3 dB to 8 dB respectively. The SSIM values of proposed method are improved from 0.01 to 0.03, from 0.02 to 0.03, and 0.03 to 0.11 respectively.

Table 1 Comparison of PSNR/SSIM of four methods on WDC

Table 2. Comparison of PSNR and SSIM of four methods on Urban

To visually show the denoising performance of the proposed method, the denoised results on WDC dataset and Indian Pines dataset under Gaussian noise variance of 40 are listed in Fig.3 and Fig.4, respectively. Obvious noise is remaining in the image of BM3D. The denoising result of NCSR is blurred and loses more details compared against other three methods. WNNM and the proposed methods both recover many details as well as effectively remove noise. But the edge in WNNM result is oversmoothed. There are additional burrs in the restored image in both WNNM and BM3D. However, these burrs are not in the original image.

Fig.3 Comparisons of visual effects with noise variance 40 on WDC

Fig.4 Comparisons of visual effects with noise variance 40 on Urban

To evaluate the spectral distortion, the spectral curves of original image and the denoised images with noise variance 30 on WDC are listed in Fig.5. Fig.5(a) and Fig.5(b) show the spectral values in the spatial position of (56,56) and (90,90), respectively. It is observed that both WNNM and BM3D show almost identical results, but the proposed method is the closest to the original data. This also proves that the spectral distortion of the proposed algorithm is lower.

Fig.5 Comparisons in terms of spectral values at different position

3.2 Comparison of visual effect of real RSI

In the experiment with real RSI, for lack of reference images, the denoising performance cannot be measured by PSNR and SSIM. Therefore, the non-reference remote sensing images quality assessment [16](NHQA for short) method based on quality-sensitive features extraction is used.

The original Indian Pines data contains noise, so one representative band 110 is used to compare the denoising effect, as shown in Fig.6.

It can be seen from Fig.6 that there is still existing much noise in the NCSR method, while BM3D method keeps nice details but over-smooth in some regions. The proposed method obtains better visual result.

Fig.6 Comparison of denoising visual effects of real data of band 110

Scores of NHQA for the four algorithms are shown in Table 3, the score of the proposed algorithm is smaller than the other three algorithms, which shows that the denoising result of the proposed algorithm is better.

Table.3 Comparison of scores of index in [16] on Indian Pines

4 Conclusion

The methods which reshaped RSI to matrix usually ignored the intrinsic structure and the strong correlation between different modes, to address this problem, the proposed algorithm considers the 3-mode unfolding of RSI, and different weights are also assigned to them depending on the contribution. Combined with the spatial and spectral nonlocal similarity, the denoising algorithm based on multi-linear weighted nuclear norm is proposed. Results from simulated data and real data show that the proposed algorithm improves the visual effect while preserving details. The PSNR and SSIM values verify that the proposed algorithm is superior to the existing advanced algorithms.