A plug-and-play hyperspectral image unmixing method based on non-negative matrix factorization

Abstract

The present invention discloses a plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition. The present invention jointly considers manually selected regularization terms and regularization terms obtained through learning, and then utilizes prior information in hyperspectral images. The method uses a denoiser to learn spectral and spatial information, avoiding manually designing regularizers and solving complex optimization problems. The present invention also uses manual regularization terms to introduce some physically clear priors in the image. Taking sparse regularization terms as an example, the present invention adds prior knowledge of abundance sparsity, illustrating that the combined use of designed regularization terms and learned regularization terms can effectively improve the unmixing accuracy. The denoiser in the present invention can not only introduce prior information, but also the denoising characteristics of the denoiser improve the robustness of the algorithm to noise.

Description

Technical Field

[0001] The invention belongs to the field of image processing, and in particular relates to a plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition. Background Art

[0002] The development of hyperspectral remote sensing imaging technology has greatly improved the spectral resolution of images. Image data can be collected and formed in a narrow band of tens of nanometers to obtain richer spatial spectrum information. The extremely high spectral resolution of hyperspectral image data gives it the ability to diagnose fine spectral features and identify and analyze the types, materials and material components of objects. It has good application prospects in mineral exploration, crop assessment, environmental monitoring and urban planning. However, due to the limitation of spatial resolution, imaging spectrometers often cannot separate all different objects, resulting in the phenomenon that multiple objects share one pixel. This phenomenon brings difficulties to the interpretation and analysis of hyperspectral remote sensing data, especially for hyperspectral applications such as fine object classification and target detection. Therefore, mixed pixel unmixing is one of the important problems that need to be solved in hyperspectral image information processing.

[0003] At present, a lot of work has been done in the field of hyperspectral unmixing, which can be divided into two categories: the first is supervised unmixing technology, which is divided into two steps: endmember extraction and abundance inversion; the second is unsupervised unmixing technology, which extracts endmembers and abundance information simultaneously through algorithms. The non-negative matrix factorization (NMF) method is a widely used unsupervised unmixing method. However, NMF is a non-convex problem and has no unique solution. To overcome this shortcoming, regularization terms can be added to the objective function to constrain the solution space and utilize the spatial and spectral characteristics of hyperspectral images. One of the typical algorithms is based on the total variation (TV) regularization term to enhance the spatial consistency of the estimated abundance. The non-local TV regularization term takes into account non-local spatial information and can make full use of similar structures in the entire image. Since mixed pixels usually contain a subset of endmembers, regularization methods to improve sparsity have also been widely used. Researchers usually use L1 or L0 norms to constrain the spatial domain and spectral domain to enhance the sparsity of the solution. It can be seen that the appropriate introduction of constraints and the clever design of regularization terms play an important role in improving the unmixing performance. However, there are still the following disadvantages:

[0004] First, it lacks flexibility and requires manual selection of the type of regularization term;

[0005] Second, different optimization problem solving methods need to be designed for different regularization terms, which is complex to solve and cannot achieve fast unmixing of hyperspectral images;

[0006] Third, the unmixing effect is poor for images with high noise content, and the calculation results are of low accuracy. To address such limitations, a plug-and-play method can be used to introduce priors into the process of establishing mathematical models through regularization and structural replacement, which can avoid artificial construction and design of image priors and does not increase the difficulty of solving mathematical optimization models. Summary of the invention

[0007] The purpose of the present invention is to overcome the above-mentioned shortcomings and provide a plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition to solve the technical problems of existing unmixing methods, lack of flexibility in regularization term selection, complex optimization problem solving process, poor image unmixing effect, and low calculation result accuracy.

[0008] In order to achieve the above object, the present invention comprises the following steps:

[0009] S1, obtaining a three-dimensional hyperspectral image to be unmixed, and reconstructing the three-dimensional hyperspectral image to be unmixed to obtain a two-dimensional hyperspectral image;

[0010] S2, performing a minimization projection on the two-dimensional hyperspectral image, and determining the number of end members of the two-dimensional hyperspectral image in the orthogonal subspace by using the sum of the noise energy after minimization projection and the projection signal error energy;

[0011] S3, construct a loss function containing a priori regularization terms;

[0012] S4, introduce auxiliary variables, transform the loss function into an equivalent form, and obtain the transformed loss function and the corresponding augmented Lagrangian function;

[0013] S5, use the VCA endmember extraction results to initialize the endmember matrix, and use the FCLS abundance estimation results to initialize the abundance matrix and auxiliary variables;

[0014] S6, estimate the endmembers according to the minimized loss function of the endmembers;

[0015] S7, abundance estimation based on the abundance minimization loss function;

[0016] S8, estimating auxiliary variables according to the minimization loss function of the auxiliary variables;

[0017] S9, loop S6 to S8 until the number of iterative cycles reaches the maximum number of cycles, and output the unmixing result of the hyperspectral image.

[0018] The three-dimensional hyperspectral image is The two-dimensional hyperspectral image is Where L represents the number of spectral channels of the three-dimensional hyperspectral image to be unmixed, H is the spatial dimension length of the three-dimensional hyperspectral image to be unmixed; W is the spatial dimension width of the three-dimensional hyperspectral image to be unmixed, and N = H × W is the total number of pixels of the two-dimensional hyperspectral image.

[0019] In S3, the loss function including the prior regularization term is:

[0020]

[0021] stE≥0,A≥0,1 T A=1

[0022] in, represents the endmember matrix, represents the abundance matrix,

[0023] In S4, the auxiliary variables are constraint The transformed loss function is:

[0024]

[0025] stE≥0,A≥0,1 T A=1,

[0026] The corresponding augmented Lagrangian function is:

[0027]

[0028] stE≥0,A≥0,1 T A=1,

[0029] In S6, the minimization loss function of the end member is:

[0030]

[0031] Where Λ is the strength of the Lagrange multiplier controlling the end member E to be "non-negative". The above equation is differentiated and the derivative is set to 0, as follows:

[0032] EAA T -RA T +Λ=0

[0033] Multiply both sides of the above equation by E. According to the KKT condition E⊙Λ=0, the update solution of the end element E is:

[0034] E←E⊙(RA T ) / (EAA T )

[0035] In S7, when estimating abundance, two augmented matrices are constructed:

[0036]

[0037] The expression of the minimization loss function for estimating abundance is:

[0038]

[0039] Where δ is the strength of the Lagrange multiplier controlling the abundance A to be “non-negative”, and Γ is the strength of the Lagrange multiplier controlling the abundance A to be “total equal to one”. The above equation is differentiated and the derivative is set to 0, as follows:

[0040]

[0041] in,

[0042]

[0043] According to the KKT condition A⊙Γ=0, the update solution of end element A is:

[0044]

[0045] In S8, it is estimated The expression of the minimization loss function is:

[0046]

[0047] Convert it to equivalent:

[0048]

[0049] Solve this problem using the denoiser to obtain the auxiliary variables

[0050] The denoisers include two-dimensional image denoiser NLM, two-dimensional image denoiser BM3D, three-dimensional image denoiser BM4D and three-dimensional image denoiser LRTDTV.

[0051] Compared with the prior art, the present invention jointly considers manually selected regularization terms and regularization terms obtained through learning, and then utilizes prior information in hyperspectral images. This method uses a denoiser to learn spectral and spatial information, avoiding manually designing regularization terms and solving complex optimization problems. The present invention also uses manual regularization terms to introduce some physically clear priors in the image. Taking sparse regularization terms as an example, the present invention adds prior knowledge of abundance sparsity to illustrate that the combined use of designed regularization terms and learned regularization terms can effectively improve the unmixing accuracy. The denoiser in the present invention can not only introduce prior information, but also the denoising characteristics of the denoiser improve the robustness of the algorithm to noise. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 It is a schematic diagram of the process of the present invention;

[0053] Figure 2This is the abundance map of simulated data at 10dB;

[0054] Figure 3 The endmembers estimated by the NML denoiser for the simulated data at 20 dB;

[0055] Figure 4 Abundance estimation results for the Jasper Ridge dataset. DETAILED DESCRIPTION

[0056] In order to make the technical problems, technical solutions and beneficial effects solved by the present invention more clearly understood, the present invention is further described in detail in the following specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0057] Step 1: Obtain the 3D hyperspectral image to be unmixed Treat unmixed 3D hyperspectral images Reconstruct the image to obtain a two-dimensional hyperspectral image. Where L represents the number of spectral channels of the three-dimensional hyperspectral image to be unmixed, H is the spatial dimension length of the three-dimensional hyperspectral image to be unmixed, W is the spatial dimension width of the three-dimensional hyperspectral image to be unmixed, and N = H × W is the total number of pixels of the two-dimensional hyperspectral image;

[0058] Step 2: Use the minimum error identification method to minimize the projection of the two-dimensional hyperspectral image R; use the sum of the noise energy after minimizing the projection and the projection signal error energy to determine the number of end members P of the two-dimensional hyperspectral image in the orthogonal subspace.

[0059] Step 3: Construct a loss function containing a priori regularization terms:

[0060]

[0061] stE≥0,A≥0,1 T A=1

[0062] in, represents the endmember matrix, represents the abundance matrix, is a sparse constraint, α represents the sparse constraint weight, and μ represents the learning prior constraint weight.

[0063] Step 4: Introduce auxiliary variables Transform the loss function into an equivalent form and constrain Get the converted loss function:

[0064]

[0065] stE≥0,A≥0,1T A=1,

[0066] The corresponding augmented Lagrangian function is obtained:

[0067]

[0068] stE≥0,A≥0,1 T A=1,

[0069] Where λ is the penalty factor.

[0070] Step 5: Use the VCA endmember extraction results to initialize the endmember matrix E, and use the FCLS abundance estimation results to initialize the abundance matrix A and auxiliary variables

[0071] Step 6: Estimate the end member E. The expression of the minimization loss function of the estimated end member is:

[0072]

[0073] Where Λ is the strength of the Lagrange multiplier controlling the "non-negative" end member E. Further, the above equation is differentiated and the derivative is set to 0, as follows:

[0074] EAA T -RA T +Λ=0

[0075] Furthermore, multiply both sides of the above equation by E, and according to the KKT condition, E⊙Λ=0. Furthermore, the update solution of the end member E is:

[0076] E←E⊙(RA T ) / (EAA T )

[0077] Step 7: To estimate the abundance A, first construct two augmented matrices:

[0078]

[0079] The expression of the minimization loss function for estimating abundance is:

[0080]

[0081] Where δ is the strength of the Lagrange multiplier controlling the "non-negative" abundance A, and Γ is the strength of the Lagrange multiplier controlling the "sum equals one" constraint of the abundance A. Further, the above equation is differentiated and the derivative is set to 0, as follows:

[0082]

[0083] in,

[0084]

[0085] According to the KKT condition A⊙Γ=0, the update solution of the end element A is further obtained as:

[0086]

[0087] Step 8: Estimate auxiliary variables estimate The expression of the minimization loss function is:

[0088]

[0089] Furthermore, it is converted into an equivalent form:

[0090]

[0091] Solve this problem using the denoiser to obtain the auxiliary variables

[0092] Step 9: Iterate steps 6, 7, and 8 until the number of iterations reaches the maximum number of iterations K, and output the unmixing result of the hyperspectral image.

[0093] Preferably, α=0.1;

[0094] Preferably, the λ=3×10 4 ;

[0095] Preferably, said σ=10;

[0096] Preferably, μ=500;

[0097] Preferably, the noise reducer is selected from a two-dimensional image noise reducer NLM, a two-dimensional image noise reducer BM3D, a three-dimensional image noise reducer BM4D, or a three-dimensional image noise reducer LRTDTV.

[0098] Specific experiments:

[0099] The demixing method and system described in the present invention are used to test a set of simulation data and Jasper Ridge data. In this experiment, a two-dimensional image denoiser NLM, a two-dimensional image denoiser BM3D, a three-dimensional image denoiser BM4D, and a three-dimensional image denoiser LRTDTV are used to verify the flexibility and effectiveness of the present invention.

[0100] The end members used to generate simulated data are extracted from the USGS spectral library. The spectra of these materials consist of 224 bands. The experiment uses the spectral curves of four pure substances as end members (P = 4), and uses the HYDRA toolkit to generate abundances containing spatial information. A total of 256 × 256 pixels are generated to evaluate the unmixing performance of the unmixing algorithm. Zero-mean Gaussian noise with signal-to-noise ratios of 5dB, 10dB, 20dB, and 30dB is added to the data.

[0101] The Jasper Ridge data consists of 224 spectral bands, covering a spectral range of 380nm to 2500nm, with a spectral resolution of up to 9.46nm. After removing channels affected by dense water vapor and atmospheric environment, 198 spectral channel information is retained. In this experiment, the number of end members is set to 4, including water, trees, rocks, and roads.

[0102] In order to prove the effectiveness of the method, four methods are used for comparison: VCA-SUnSAL-TV, CoNMF, TV-RSNMF and NMF-QMV. The method proposed in the present invention is PNMF-NLM, PNMF-BM3D, PNMF-BM4D, PNMF-LRTDTV.

[0103] In this experiment, the root mean square error (RMSE) was used to evaluate the performance of hyperspectral image abundance estimation, and the root mean square error comparison results of the estimated simulated abundance were obtained. The spectral angular distance (SAD) was used to evaluate the hyperspectral image endmember extraction results, and the spectral angular distance comparison results of the extracted endmembers were obtained, as shown in Table 1.

[0104] Table 1 Comparison of RMSE and SAD of simulated data

[0105]

[0106]

[0107] The results show that the proposed method achieves the best RMSE and PSNR results and the lowest SAD value compared to other methods. This highlights the effect of sparse regularization and the superiority of the prior learned by the denoiser. In addition, thanks to the use of the denoiser, the unmixing effect of the proposed method is more significantly enhanced when the noise level is high. This shows that the proposed method is robust to noise.

[0108] As attached Figure 2 As shown, attached Figure 2 The abundance mapping diagram of the simulated data under the 10 dB condition is given, and it can be seen that the abundance estimation result of the method of the present invention has less noise and is closer to the ground truth.

[0109] As attached Figure 3 As shown, attached Figure 3 The endmembers estimated by using the NML denoiser for the simulated data at 20 dB are given. It can be seen that the method proposed in the present invention can achieve good endmember estimation results.

[0110] As attached Figure 4 As shown, attached Figure 4 The abundance estimation results of the Jasper Ridge dataset are given. The results show that the proposed method has lower noise and better smoothness.

[0111] Compared with the existing methods, the unmixing method system described in the present invention jointly considers manually selected regularization terms and regularization terms obtained through learning, and then utilizes the prior information in the hyperspectral image. This method uses a denoiser to learn spectral and spatial information, avoiding the manual design of regularization terms and solving complex optimization problems. The present invention also uses manual regularization terms to introduce some physically clear priors in the image. Taking the sparse regularization term as an example, the addition of the prior knowledge of abundance sparsity shows that the combined use of the designed regularization term and the learned regularization term can effectively improve the unmixing accuracy. The denoiser in the present invention can not only introduce prior information, but also the denoising characteristics of the denoiser improve the robustness of the algorithm to noise.

[0112] The above embodiment is only one of the implementation methods that can realize the technical solution of the present invention. The scope of protection claimed by the present invention is not limited only to this embodiment, but also includes changes, replacements and other implementation methods that can be easily thought of by any technician familiar with the technical field within the technical scope disclosed by the present invention.

Claims

1. A plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition, characterized in that: The following steps are involved: S1, obtaining a three-dimensional hyperspectral image to be unmixed, and reconstructing the three-dimensional hyperspectral image to be unmixed to obtain a two-dimensional hyperspectral image; S2, performing a minimization projection on the two-dimensional hyperspectral image, and determining the number of end members of the two-dimensional hyperspectral image in the orthogonal subspace by using the sum of the noise energy after minimization projection and the projection signal error energy; S3, construct a loss function containing a priori regularization terms; S4, introduce auxiliary variables, transform the loss function into an equivalent form, and obtain the transformed loss function and the corresponding augmented Lagrangian function; S5, use the VCA endmember extraction results to initialize the endmember matrix, and use the FCLS abundance estimation results to initialize the abundance matrix and auxiliary variables; S6, estimate the endmembers according to the minimized loss function of the endmembers; S7, abundance estimation based on the abundance minimization loss function; S8, estimating auxiliary variables according to the minimization loss function of the auxiliary variables; S9, loop S6 to S8 until the number of iterative cycles reaches the maximum number of cycles, and output the unmixing result of the hyperspectral image.

2. The plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition according to claim 1, characterized in that: The three-dimensional hyperspectral image is The two-dimensional hyperspectral image is Where L represents the number of spectral channels of the three-dimensional hyperspectral image to be unmixed, H is the spatial dimension length of the three-dimensional hyperspectral image to be unmixed; W is the spatial dimension width of the three-dimensional hyperspectral image to be unmixed, and N = H × W is the total number of pixels of the two-dimensional hyperspectral image.

3. The plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition according to claim 1, characterized in that: In S3, the loss function including the prior regularization term is: s.t.E≥0,A≥0,1 T A=1 in, represents the endmember matrix, represents the abundance matrix, 4. The plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition according to claim 1, characterized in that: In S4, the auxiliary variables are constraint The transformed loss function is: The corresponding augmented Lagrangian function is:

5. The plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition according to claim 1, characterized in that: In S6, the minimization loss function of the end member is: Where Λ is the strength of the Lagrange multiplier controlling the end member E to be non-negative. The above equation is differentiated and the derivative is set to 0, as follows: EAA T - RA T +Λ=0 Multiply both sides of the above equation by E, according to the KKT condition The update solution of end member E is:

6. The plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition according to claim 1, characterized in that: In S7, when estimating abundance, two augmented matrices are constructed: The expression of the minimization loss function for estimating abundance is: Where δ is the strength of the Lagrange multiplier controlling the abundance A to be non-negative, Γ is the strength of the Lagrange multiplier controlling the abundance A and the constraint of unity. The above equation is differentiated and the derivative is set to 0, as follows: in, According to KKT conditions The update solution of end member A is:

7. The plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition according to claim 1, characterized in that: In S8, it is estimated The expression of the minimization loss function is: Convert it to equivalent: Solve this problem using the denoiser to obtain the auxiliary variables 8. The plug-and-play hyperspectral image unmixing method based on non-negative matrix decomposition according to claim 7, characterized in that: The denoisers include two-dimensional image denoiser NLM, two-dimensional image denoiser BM3D, three-dimensional image denoiser BM4D and three-dimensional image denoiser LRTDTV.

Images