Abstract

A new nonconvex smooth rank approximation model is proposed to deal with HSI mixed noise in this paper. The low-rank matrix with Laplace function regularization is used to approximate the nuclear norm, and its performance is superior to the nuclear norm regularization. A new phase congruency lp norm model is proposed to constrain the spatial structure information of hyperspectral images, to solve the phenomenon of “artificial artifact” in the process of hyperspectral image denoising. This model not only makes use of the low-rank characteristic of the hyperspectral image accurately, but also combines the structural information of all bands and the local information of the neighborhood, and then based on the Alternating Direction Method of Multipliers (ADMM), an optimization method for solving the model is proposed. The results of simulation and real data experiments show that the proposed method is more effective than the competcing state-of-the-art denoising methods.

1. Introduction

Due to the influence of many factors in the process and transmission of hyperspectral images, the acquired hyperspectral images often contain some complex mixed noises, including gauss noise, salt-and-pepper noise, and dead-line noise. It is very difficult to analyze and apply the hyperspectral image to high-level applications. In the early days, a great number of methods were proposed to remove noise, such as Fourier Transform, Wavelet Transform, nonlocal means (NLM) filter, block-matching, and 3D filtering (BM3D). However, most of the above methods require some prior knowledge of noise and can only deal with one or two types of noise. In addition, HSI data in the real-world are often mixed into the real data by various noise combinations. To solve this problem, some multidimensional methods are proposed to deal with both spectral and spatial information. In [1], a multidimensional wiener filtering (MWF) algorithm was proposed, which represents the histogram of an image as a three-dimensional tensor and uses tensor analysis to remove the noise. In addition, spectral and spatial information [2] 3D wavelets have obtained good performance.

In recent years, low-rank Matrix Recovery (LRMR) [3] has been used in HIS denoising. Different from the traditional method, LRMR can process different types of noise without any prior information of noise; many methods of HSI image denoising [47] based on low rank and spectral correlation have been proposed. Due to the difficulty in solving the low-rank constraint model directly, the kernel norm is used to approximate the low-rank model, and the good results of HSI mixed noise removal are obtained. The representative methods include weighted low-rank model (WLRM) [8], rank minimization (RM) [9], structure tensor total variance weighted nuclear norm minimization (STTV-WNNM) [10], weighted nuclear norm minimization (WNNM) [11], and nonlocal low-rank approximation (NLRA).

These methods mainly have taken advantage of the low-rank feature of low-rank description image with nuclear norm approximation and preserving edge structure of image with total variation regularization. However, the low-rank model based on nuclear norm approximation is not a good approximation rank function but lacks image edge sparsity and structure smoothness regularization. If NSS prior cannot be well used to represent the structure of the image in the low-rank model, the important structure will inevitably be lost, and an “artificial artifact” will appear.

Therefore, it is necessary to approach the rank function better. To overcome the above limitation, the method based on nonconvex smooth rank approximation (SRA) is studied in this paper. The key idea is to use nonconvex smooth functions to approximate rank functions directly and to provide a more rigorous approximation than traditional methods.

In addition, the method based on the fractional band total variation regularization has the following defects: (1) The total variation regularization is essentially a first-order partial differential equation, which is a kind of ill-posed first-order inversion problem; if the original image is disturbed, it is possible to have a large oscillation on the partial derivative. (2) The total variation model is based on the assumption that the image is piecewise linear continuous, and the sharp edges of the image can be effectively preserved. However, only the gradients of adjacent pixels are used in the total variation model discretization; it cannot describe the true edge and structure information of the image accurately and effectively, which may lead to the phenomenon of “artificial artifact.” Total generalized variation (TGV) [12], high-order total variation (HOTV) [13], and Schatten p-norm constraint models [14] are used to deal with the “artificial artifacts” caused by the total variation model; it is more accurate to describe image edge information in a natural image, but this kind of method is not suitable for hyperspectral image processing.

Aiming at the shortcomings of the existing subband TV regularization and low-rank tensor denoising methods, which cannot effectively utilize the local neighborhood information, it is easy to cause the “artificial artifact” phenomenon, especially in the curved edge. The gradient cannot accurately describe the true structure of the image. Therefore, in this paper, a new phase congruency lp norm constraint is proposed to constrain the spatial structure information of hyperspectral images, to solve the phenomenon of “artificial artifact” in the process of hyperspectral image denoising.

2. State of the Art

2.1. Phase Congruency

Phase congruency (PC) [12] is a phase-based frequency-domain feature detector proposed by Morrone and Owens, which can detect a wide range of visual features and is invariant to local smooth illumination variations. Unlike the gradient in the spatial domain, the PC uses the frequency domain to get the spectrum information at the maximum overlap of the phases, for the high-order edge and corner, line, step edge and roof Ridge, and other visual sensitive important image feature capture more complete. In [12], the phase consistency of signal x is defined by Fourier series expansion shown aswhere is the Fourier series of the n term, and is the local phase information of the Fourier series for signal x, and the PC value is the maximum value of equation (1) for the parameter , which is the weighted average of all the amplitudes corresponding to the Fourier terms. Because the PC function defined by equation (1) is greatly influenced by noise, it is easy to produce edge offset, which leads to the loss of part of image structure and contour features.

A phase congruency (MPC) [13] was presented for extracting three orthogonal characteristic components of an image, namely, amplitude, direction, and phase, from a single-pass signal theory and applied to image quality evaluation [14] and image denoising applications, and the satisfactory results of time efficiency and detection effect were obtained.

The phase congruency based on single-pass signal theory is defined as follows:where is the weight function, is an approximate gain factor for sharpening the edge response, and the value range [1, 2], T =  is the compensation noise effect. W (X) is defined bywhere is the gain factor, is the cut-off value of the filter response expansion, and is the fractional scatter measure, whose value is divided by the response amplitude and the highest response value , where is the scale quantity.

Phase consistency lp norm (PCSP) is defined as the Schatten p norm of the PC value on the definition of the singleton phase consistency function in equation (2).where and are the height and width of the image, respectively, . Equation (4) is difficult to optimize the Schatten p norm due to the existence of convolution kernel operation. To solve this problem, equation (4) needs to be transformed into an equivalent form that is easy to be solved.

Given a Matrix , then the singular value of the Matrix can be decomposed into , where U and V are the left and right singular value eigenvectors of , respectively, and the diagonal element of the Matrix is the singular eigenvalue. The Schatten P norm of a Matrix is defined bywhere is the ith singular eigenvalue of the Matrix Y, corresponding to the (i, i) element value. From equation (5), the PCTV equivalent of equation (4) can be defined by

2.2. Nonconvex Low-Rank Approximation Functions

Because the nuclear norm cannot approximate the rank function well, it is unreasonable to replace the rank function with the nuclear norm directly. In the field of image restoration, the method of approximating rank minimization by nonconvex substitution has received extensive attention. The common nonconvex functions and corresponding hypergradients are shown in Table 1.

To provide more rigorous approximations than the nuclear norm, we have conducted numerical experiments on one-dimensional data. We selected five experiments of rank approximation of regularization terms with the best performance, as shown in Figure 1. As can be seen from Figure 2, nonconvex optimization is often superior to convex optimization.

As can be seen from Figure 2, the nuclear norm deviates from the true rank, so all singular values are treated equally. The weighted nuclear norm and the Garman norm are the neutralization between the rank minimization and the kernel norm, which can increase the penalty to the small value and decrease the penalty to the large value. The Logdet norm is poor at small singular values, especially those close to 0. The exponential function is used to deal with different singular values, so that Laplace modules and real rank have obvious consistency. To maintain accuracy and speed, the Laplace norm is the best choice for approximating the actual rank of the NSS matrix composed of a similar patch.

In this section, the Laplace function is chosen to approximate the rank function. Compared with other tensor nuclear norms, the Laplace function proxy norm is a better method to measure the rank of the tensor. We introduce the proposed proxy into a low-rank tensor separation model and solve the model by using the ADMM algorithm, which can effectively complete the missing elements in the tensor, to achieve the goal of hyperspectral image denoising. A large number of experiments show that this method is better than the existing methods.

3. HSI Denoising Based on Nonconvex Low-Rank Tensor Approximation

3.1. The Motivation of the Model

The method based on tensor nuclear norm minimization and TV regularization has achieved some noise reduction results in hyperspectral denoising applications, but due to the limitation of TV total variation and the fact that the kernel norm cannot accurately describe the feature of low tensor rank, as a result, its noise removal performance is limited, and it is difficult to meet the demand of hyperspectral image denoising ability for practical applications. Based on the minimization of nonconvex low-rank approximation, the weight of singular value can be adjusted according to the main features of the image, and the global low-rank characteristic of the HSI image can be used well. In addition, phase congruency can describe and capture image visual features more comprehensively than TV constraints based on gradient information. Therefore, in this subsection, nuclear norms are replaced by tensor-based nonconvex proxy functions. Because of the limitation of low-rank feature detection, phase consistency is used to preserve image edge structure. Therefore, a new HSI denoising method based on nonconvex low-rank tensor approximation and phase consistency constraints is proposed in this section.

3.2. Model Description
3.2.1. Constraint Term of Nonconvex Function and Its Properties

Before elaborating on the denoising model, we first review the definition of the normalized Laplace function, which is shown asand .

Laplace functions have the following useful properties.

Proposition 1.

Proof. for
And
So come to the conclusion.

Proposition 2. ,
For any tensor, and hold.

Proof. and are orthogonal matrices, because the singular value of the Matrix does not change with the multiplication of the orthogonal Matrix, and and have the same individual value; thus, can be obtained.

Proposition 3. holds for any tensor , if and only if .

3.2.2. NLRTAPC Model

To eliminate the mixed noise of HSI, considering the global low-rank and local piecewise smooth characteristics of the band image, and combined with the PCTV regularization term in the nonconvex low-rank tensor approximation model, the mixed noise reduction method based on Nonconvex Low-Rank Tensor Approximation and Phase Consistency for Mixed Denoising (NLRTAPC) is presented bywhere , , are nonnegative parameters, is the restoring low-rank tensor, is the sparse noise component, and is the gauss noise component. The first term of the model is the low-rank tensor approximated by a nonconvex function with formula 2.1. The second term is the Schatten p norm constraint term based on phase congruency, and the last two are for noise error constraints. The NLRTAPC model not only preserves the low-rank features of the image more accurately than the tensor kernel norm model, but also uses the phase consistency to preserve the structural information such as edge and texture in the band restoration component.

3.3. Solution of NLRTAPC Model

NLRTAPC model is a nonconvex optimization problem, which is solved by decomposition subproblem. The ADMM algorithm can be used to solve the corresponding model iteratively. By introducing the secondary variable , the NLRTAPC model can be expressed as

Since the four variables , , and in the model are separable, the objective function can be solved by ADMM with formula 2.2, which is expressed aswhere and are the Lagrange multiplier with and constraints, respectively, and is the penalty parameter. Iterating through the NLRTAPC model according to the ADMM framework can be broken down into the following five steps, each of which solves the corresponding variable, as follows:(1)fix other tensor variables , , , and , and update the estimator . Under the iterative framework of (k + 1), the estimated true image is as follows:where .According to Proposition 2, given any given tensor , it can be decomposed into by T-SVD, and can be solved by a weighted tensor singular threshold method.The tensor rank approximation problem of nonconvex normalized is used. The problem can be solved by the singular threshold operator of tensor t-SVD decomposition. When is decomposed into by T-SVD, the optimal solution is , where is the diagonal tensor of f, and the elements of the diagonal tensor are obtained by Fourier domain calculation; that is, .(2)fix the other tensor variables , , , and , update the lp norm constraint of PC, and save the spatial smoothness and edge information.where . First, we use singular value decomposition tensor with T-SVD, . Then, lp norm approximation is performed for the f diagonal tensor , denoted as . Finally, is reconstructed, and , where dialog (·) means that the vector elements are expressed in a matrix form, is a diagonal matrix convenient operation.(3)To remove impulse noise, other tensor variables , , , and are fixed, and is updated.where , and then the closed-form solution of formula (13) is obtained by using the contraction operator of Matrix elements.(4)To remove gauss noise, fix , , , , and , and update .where , and is the Lagrange multiplier of the k-th iteration. Formula (14) is a standard least square regression problem, which is easy to solve.(5)Update Lagrange multiplier and penalty parameters.where is the contraction parameter; set to accelerate the rate of convergence and to be the maximum value of .

The NLRTAPC model is described by ADMM in Algorithm 1. To obtain good phase consistency, the front slice of the HSI image is a waveband image, and the phase consistency is extracted in Matrix form to get pc ().

input: Observational data , Index set , Parameters
 Initializing: , , k = 0, num = 1000 and ,i = 0.
(1)when or k < num.
(2)k ++;
(3);
(4)Cycle (I = 1) until the maximum band number B of the observed data
(5)The updated estimate is calculated according to formula (11);
(6)The updated estimate is calculated according to formula (12);
(7)The updated estimate is calculated according to formula (13);
(8)The updated estimate is calculated according to formula (14);
(9)According to formula (15), the Lagrange multiplier and the penalty parameters are updated respectively.
(10)End;
(11)Using to restore the tensor structure , , and , ;
(12);
 End;
 output: denoising tensor .

4. Experimental Results and Analysis

In this section, experiments were carried out to demonstrate the mixed noise removal capability of our model. To better illustrate the superiority of the combination of the nonconvex smooth rank approximation model and the lp norm constraint, the validity of the proposed method is verified by simulation and real experiments, and the quantitative and visual performance of the four advanced HSI denoising methods are compared with the denoising results of this method. These methods include block-matched 4D filtering (BM4D), Nonlrma [14], global spatial something spectral total variation (GSSTV) [15], and weighted group sparse regularized low-rank tensor decomposition (WGLRTD). The code of all comparison methods is Matlab Code. All experiments of comparison method are carried out on Intel Core i7-4970 CPU (3.60 GHz) and 16 GB RAM computer using Matlab R2018a. The parameters of these methods in the experiment are set according to the method suggestions to obtain the best performance.

4.1. Simulation Data Experiment

Two data sets, namely, WDC and Indian data sets, are used in the simulation data experiment. The WDC data set was collected by a hyperspectral digital image acquisition experiment (HYDICE) sensor at a mall in Washington, D. C., at 191 wavelengths. The “Indian” is a collection of 145 ∗ 145 pixels and 224 spectral reflectivity bands collected by the AVIRIS sensors over the Indian pine proving ground in northwestern Indiana. Since the data set contains a portion of atmospheric absorption bands that are not useful for subsequent applications, the 200-band noise-free images were selected for the experiment. Experiment parameters are set by .

To simulate the complex noise situation in the real scene, Matlab is used to generate 4 kinds of mixed noise in the two clean HSI data sets. For noise type 1 to noise type 9, the intensity of gauss noise, random noise, and band noise is shown in Table 2.

4.2. Analysis of Simulation Experiment Results

There are four methods in the experiment, BM4D, NonLRMA, GSSTV, and WGLRTD. Three noise types (noise types 1, 4, and 7) are chosen randomly to compare the performance of each algorithm model. Figure 1 shows the result of a 58 band data recovery from a mall in Washington, D. C., under noise Type 4 conditions. For better visual contrast, some areas of the image are deliberately magnified, as shown in Figure 1:

As can be seen from Figure 1, BM4D has a serious mold and loss of image detail due to oversmoothing. NonLRMA and GSSTV still retain a small amount of noise, and the image is darkened due to the offset of the whole pixel value. In the enlarged area, the WGLRTD left a lot of random noise, and the details of the image were not well preserved. On the contrary, the NLRTAPC model can eliminate all the mixed noises and preserve the edges and details of the image effectively, which shows that the NLRTAPC method has better performance in recovering WDC data set than other current typical methods.

Figure 3 shows a comparison of restoration at the 165 bands of the Indian dataset under noise type 1, from the overall view of Figure 3; since the damage to the image structure caused by noise type 1 is not very serious, several methods can ensure that the overall structure recovers well, but it can be found that NLRTAPC, the Algorithm model designed in this chapter, has a clear structure in the edge region, while it still keeps the smooth effect in the smooth region and presents a better visual effect.

In Figure 4, the 165 band data recovery results are for the Indian Dataset. By comparing and analyzing the recovery results of all the comparison methods, a more direct-viewing result is obtained. A more direct result can be obtained from Figure 4. Because the original image has been destroyed seriously, the recovery result of the GSSTV method is ambiguous, the recovery result of the NonLRMA method still has an “artificial gradient” phenomenon, and bM4D and WGLRTD methods have some modulus on the edge of the image structure. In this chapter, HLRTD-SSTV Algorithm is presented to restore the result of the overall structure clear and smooth, and the visual effect is good; overall, it is superior to other methods. To demonstrate the excellent performance of the NLRTAPC algorithm model, two objective quantitative evaluation indexes are introduced to prove the performance of the NLRTAPC algorithm model in all simulation experiments, which are average Peak signal-to-noise ratios (MPSNR) of all bands and average structural similarity (MSSIM) of all bands.

Tables 3 and 4 show the MPSNR and MSSIM values of the Indian and WDC datasets under various mixed noise types, respectively. In the experiment, MPSNR and MSSIM were evaluated as the average of all bands of the data set, and NLRTAPC gets good indicator values. From the evaluation data of nine different noise types, the NLRTAPC model achieves the best performance in most cases. With the increase of noise level, the performance degradation is relatively slow compared with other methods. In particular, the MSSIM Algorithm can improve the performance obviously, and the key is that the phase consistency constraint is used to capture the edge information of the image more comprehensively and more accurately than TV.

4.3. Real Dataset Experiment

In the simulation experiment, the performance of the NLRTAPC under various noise conditions is evaluated by visual, quantitative, and qualitative methods. This section performs a hyperspectral digital image acquisition experiment (HYDICE) on a real urban dataset to verify the effectiveness of NLRTAPC’s real hyperspectral image restoration. The whole dataset is applied in the experiment of denoising algorithm. Some brands of the dataset are heavily polluted by the atmosphere, hygroscopicity, and some thermal noise; in addition, some bands are polluted by striations, gauss noise, and random noise. As a preprocessing, the pixel values for each band are normalized to [0,1].

Figure 5 shows the results of band 76 recoveries in this data. As can be seen from the graph, the original band image is corrupted by various noises, including gauss noise and unknown noise. After using different HSI reconstruction methods, the noise was removed, BM4D mode was burnt, and the WGLRTD method still had some noise in the denoising band image. From the enlarged area and the whole image, it can be seen that the NLRTAPC Algorithm in this chapter is effective to remove the mixed noise while preserving the image structure. Figure 6 shows the results of band 136 recoveries in this data. In contrast, the NLRTAPC approach can completely suppress all kinds of noise, effectively preserving details. Compared with other methods, the image restored by this method is smoother. These visualization results show that the NLRTAPC method can remove more complex noises hidden in HSI than other methods.

To demonstrate the results of image reconstruction, Figure 7 uses different methods to estimate spectral characteristic curves for 76 band images in real-world urban datasets. In Figure 7, the ordinal of bands is shown on the horizontal axis, and the average of the estimated spectral characteristics of each band is shown on the vertical axis. As shown in Figure 5(a), there are many fluctuations in the spectral characteristic curve of the original band image due to striations and other noises. It can be observed that the comparison method suppressed the volatility to some extent after the recovery of the different methods shown in Figures 5(b)5(f). The results of spectral characteristic curves show that the NLRTAPC recovery curve is more smooth than other methods, which indicates that NLRTAPC has a better denoising effect. As shown in Figures 7(b) and 7(c), there are also small fluctuations in the curve, indicating that some of the mixed noise in the image remains in the image.

5. Conclusion

In this paper, a hybrid noise removal method based on nonconvex low-rank tensor approximation is proposed, based on the optimization method ADMM proposed to solve the model by using the structural information of all bands and the local information of the neighborhood. NLRTAPC model can make use of the structure information of all bands and local neighborhood information more effectively, promote the removal of mixed noise, and greatly alleviate the “artificial ladder” phenomenon. The results of simulation and real data experiments show that the proposed method is more effective than the competing state-of-the-art denoising methods. The future improved direction of this approach is to optimize parameters of NLRTAPC to improve model adaptability.

Data Availability

The data used to support the findings of this study were supplied by Li bo under license and so cannot be made freely available. Requests for access to these data should be made to Li Bo ([email protected]).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was jointly funded by the Sichuan Science and Technology Program, Grant/Award Numbers: 2019ZYZF0169; the Science and Technology Support Project of Panzhihua City of Sichuan Province of China under Grant 2021CY-S-6, and Tourism resources big data application technology innovation team project.