Robust non-negative matrix factorization hyperspectral image demixing method

Abstract

The invention discloses a robust non-negative matrix factorization hyperspectral image demixing method. The method comprises the following steps: step 1, calculating a noise estimation value and a noise estimation matrix of a hyperspectral image; step 2, constructing an improved weighted residual non-negative matrix factorization model based on the noise estimation value and the noise estimation matrix; 3, introducing a second-order adjacency relation to obtain a high-order graph, and carrying out Laplacian regularization to obtain a regularization term; 4, introducing a high-order graph Laplacian regular term to obtain a robust non-negative matrix factorization model based on global structure constraint; and step 5, realizing hyperspectral image unmixing based on the robust non-negative matrix factorization model of the global structure constraint. According to the method, a more robust and more accurate unmixing effect can be obtained, the internal structure of the hyperspectral image can be better protected, and the method is more suitable for scenes with more complex ground features, richer details or higher unmixing precision requirements.

Description

Technical Field

[0001] This invention relates to the field of hyperspectral image processing, and in particular to a robust nonnegative matrix factorization method for demixing hyperspectral images. Background Technology

[0002] Hyperspectral images (HSI) are acquired by imaging spectrometers, recording remote sensing data of reflectance across hundreds or even thousands of spectral bands in a single acquisition of a surveyed scene. Due to its rich spectral information, HSI has been applied in numerous fields, including agriculture, military defense, surface mineral exploration, and biomedicine. Unlike traditional RGB images, HSI can detect not only the two-dimensional spatial information of a scene but also its rich one-dimensional spectral information.

[0003] In hyperspectral unmixing research, unmixing methods are mainly divided into supervised and unsupervised methods depending on the application scenario. Supervised unmixing methods achieve land cover identification by utilizing existing spectral information of land cover in spectral libraries and a large number of labeled samples. However, this method requires that the land cover type be known and that there are sufficient labeled samples, and it is usually applied in fields such as agricultural monitoring and urban planning. Due to the complexity of the actual environment, some land cover types may not be included in existing spectral libraries, which limits the applicability of this method and has thus promoted the development of unsupervised unmixing techniques. Unlike supervised methods, unsupervised unmixing methods can perform endmember extraction and abundance inversion simultaneously when the land cover spectrum is unknown. Therefore, this method can save the overhead of labeled samples and plays an important role in the field of hyperspectral unmixing. Among the many unsupervised unmixing methods, the nonnegative matrix factorization (NMF) method based on linear mixture models has been widely used in the field of hyperspectral image unmixing due to its clear physical meaning, strong interpretability, and high computational efficiency. NMF can extract the spectral features and proportions of ground features simultaneously by decomposing the HSI matrix into two non-negative matrices, simplifying the decontamination process. Summary of the Invention

[0004] Purpose of the invention: The technical problem to be solved by the present invention is to provide a robust nonnegative matrix decomposition hyperspectral image unmixing method to address the shortcomings of the existing technology.

[0005] To address the aforementioned technical problems, this invention discloses a robust nonnegative matrix factorization method for hyperspectral image unmixing, comprising the following steps:

[0006] Step 1: Calculate the noise estimate and noise estimation matrix of the hyperspectral image;

[0007] Step 2: Construct an improved weighted residual nonnegative matrix decomposition model based on the noise estimate and the noise estimation matrix;

[0008] Step 3: Introduce second-order adjacency relations to obtain a higher-order graph and perform Laplace regularization to obtain regularization terms;

[0009] Step 4: Introduce the higher-order graph Laplace regularization term to obtain a robust non-negative matrix decomposition model based on global structural constraints;

[0010] Step 5: Implement hyperspectral image demixing using a robust nonnegative matrix factorization model based on global structural constraints.

[0011] Step 1 is as follows:

[0012] Step 1-1, assuming the acquired hyperspectral image data containing L bands and N pixels is as follows:

[0013] Y = [y1,...,y] N ]∈R L×N

[0014] Define matrix Z = Y T =[z1,z2,…,z L ]∈R N×L , where vector y i Let z represent the spectrum of the i-th pixel. i This represents the reflection value of the i-th band in Y;

[0015] Define an N×(L-1) dimensional matrix Z di =[z1,z2,…,z i-1 ,z i+1 ,…z L ], Z di Let Z be the matrix after removing the i-th column, (·) T Represents the transpose of a matrix;

[0016] Steps 1-2: Based on the correlation between hyperspectral image bands, the reflectance value z i The remaining L-1 vectors can be linearly represented as z. i =Z di β i +n i , where β i ∈R L-1 Let n represent the regression vector. i Represents the error vector;

[0017] Steps 1-3, the minimum mean square estimate of the regression vector is Where i = 1, ..., L, the estimated noise value is obtained from this. The noise estimation matrix is ​​then...

[0018] Step 2 is as follows:

[0019] Step 2-1: The correlation matrix of the noise is calculated using the noise estimation matrix.

[0020]

[0021] Where N represents the number of pixels in the hyperspectral image;

[0022] Step 2-2, the noise power of each band is P N =[C 11 C 22 ,…,C LL ], where C ii This represents the noise power of the i-th band;

[0023] Steps 2-3: Due to the randomness and uncertainty of noise, the noise power in different bands may vary significantly. To effectively allocate the weight of each band, the noise power is normalized as follows:

[0024]

[0025] The band weights are then obtained as follows:

[0026]

[0027] Among them, b i Let α represent the weight of the i-th band, i = 1, ..., L, and α be the scaling parameter. Considering the rationality of the weight, α is set to 8. Through the above strategy, the weight of the band can be adjusted in the range of (0,1). The weight of the band with more severe noise interference is closer to 1, so that the model will focus on the band in subsequent algorithms.

[0028] Steps 2-4 involve obtaining an initial estimate of the endmember matrix E of the hyperspectral image through vertex component analysis. After successful extraction of the endmember matrix E, the initial value of the abundance matrix A is obtained by minimizing the error between the hyperspectral image data Y and the spectrum EA reconstructed by a linear mixture model based on endmember spectral and abundance estimates using the fully constrained least squares method. The fully constrained least squares method can simultaneously satisfy the constraints of non-negativity of abundance and the sum of abundance being one, where E∈R. L×n An endmember matrix is ​​a matrix that represents the unique spectral characteristics of a ground feature or substance, where A ∈ R. n×N The abundance matrix is ​​a matrix that represents the spectral band information in an image, where n represents the number of endmembers.

[0029] Steps 2-5 introduce noise into the weighted residual nonnegative matrix factorization model and assign the initial bands of the hyperspectral image by measuring the degree of noise interference.

[0030] The residual weight matrix of the weighted residual nonnegative matrix factorization model improved by noise estimation is:

[0031] W N =diag[w N,11 ,w N,22 ,…,w N,LL ]

[0032] Define its diagonal elements as:

[0033]

[0034] Among them, w N,ii Let R represent the weight of the i-th band, k represent the number of iterations, k = 0 represent the initial state, ||·||² represent the Euclidean norm, and R is the weight of the i-th band. (k) =YE (k) A (k) R represents the residual of the k-th iteration. i,: Let R represent the i-th row, where μ is a constant characterizing the degree of exponential decay of the residual, and μ = 30.

[0035] Steps 2-6 introduce the improved weighted residual nonnegative matrix decomposition model based on noise estimation:

[0036]

[0037] Among them, ||·|| F Describing the Frobenius norm, and Let N and N represent vectors consisting entirely of 1s, respectively.

[0038] The improved weighting mechanism provides initial weights for each band before iteration, and combined with dynamic adjustments during the iteration process, it can effectively reduce the impact of noise and improve the robustness and unmixing efficiency of the model.

[0039] Step 3 specifically involves:

[0040] Step 3-1: Given a hyperspectral image matrix Y, construct adjacency relationships using pixels as vertices. The first-order adjacency weights are expressed as:

[0041]

[0042] Among them, y i Let y represent the i-th pixel. j Let θ represent the j-th pixel, and θ is a parameter that controls the degree of weight change. θ takes a non-zero real number.

[0043] Step 3-2: To more effectively learn the global structure of the hyperspectral image, the weights of the first-order adjacency matrix are improved by combining pixel spatial distance, making the algorithm pay more attention to pixels outside the local neighborhood. The improved first-order adjacency weights are:

[0044]

[0045] Where d(y) i ,y j θ1 and θ2 represent the Euclidean distance between two pixels, respectively, which control pixel saliency and the degree of change in spatial distance. We set θ1 = 1. For parameter θ2, to measure the influence of pixel distance, assuming the spatial size of the hyperspectral image is r × c, we set θ2 = (r × c) / (r × c) × c. 2 +c 2 ) / (e-1), where e is a constant and takes the value of an infinite non-repeating decimal, approximately 2.71828182845904523536. The improved first-order adjacency weight assigns greater weight to similar pixels that are farther apart, thus representing a closer relationship between the two pixels.

[0046] Step 3-3, define the first-order Laplace matrix as:

[0047] H1 = D1 - W1

[0048] Where D1 is the degree matrix, W1 is the improved first-order adjacency weight matrix;

[0049] Steps 3-4 define the second-order adjacency weight as follows:

[0050]

[0051] Among them, W 1,i Represents the i-th row of W1, W 1,j In the j-th row of W1, when two pixels have a large number of interconnected adjacent pixels, W... 2,ij A value >0 indicates a significant second-order adjacency relationship between the two pixels. When there are no mutually connected adjacent pixels, then W... 2,ij A value of zero indicates that the two pixels do not have a second-order adjacency relationship;

[0052] Further derivation of the second-order Laplace matrix is ​​as follows:

[0053] H2 = D2 - W2

[0054] Where D2 is a diagonal matrix. W2 is a second-order adjacency weight matrix;

[0055] Steps 3-5: Combining first-order and second-order adjacency relations (represented by first-order and second-order Laplace matrices respectively), the higher-order Laplace matrix H... h Diagonal matrix D h and higher-order adjacency weight matrix W h They are defined as follows:

[0056] H h =ρH1+(1-ρ)H2

[0057] D h =ρD1+(1-ρ)D2

[0058] W h =ρW1+(1-ρ)W2

[0059] Wherein, ρ is a parameter that balances the impact of first-order and second-order adjacency relations on the overall performance, and its value ranges from [0,1].

[0060] Steps 3-6, based on the higher-order Laplace matrix H h Diagonal matrix D h and higher-order adjacency weight matrix W h The abundance vectors a corresponding to the i-th and j-th pixels :i and a :j The corresponding minimized regularization term is as follows:

[0061]

[0062] Among them, W h,ij W represents the higher-order adjacency weight matrix h The weight value in the i-th row and j-th column, D h,ii Represents a diagonal matrix D h The weight value in the i-th row and i-th column, Tr(·) represents the trace of the matrix.

[0063] Step 4 is as follows:

[0064] Step 4-1: Calculate the estimated endmember matrix E using an unconstrained nonnegative matrix factorization model. est And estimating the spectrum, where, in the traditional nonnegative matrix factorization model, both the endmember matrix and the abundance matrix must satisfy the nonnegativity constraint, and the unconstrained nonnegative matrix factorization model refers to the variant of the nonnegative matrix factorization model that does not impose endmember and abundance nonnegativity constraints;

[0065] Step 4-2: Calculate the endmember smoothing weight matrix Q, and the weight Q of the endmember smoothing level parameter corresponding to the element in the i-th row and j-th column. ij The formula is:

[0066]

[0067] Where η is a parameter controlling the smoothness, and η takes the form of a non-zero real number, E est；i,j Indicates the estimated endmember matrix E est The element in the i-th row and j-th column is the weight Q of the endmember smoothing level parameter. ij E is measured by the difference between the estimated spectrum and the mean of its two preceding and following reflectance values. est；i+1,j The closer to (E) est；i+2,j +E est；i,j The smoother the surface at this point, the higher the weight Q. ij The larger it is;

[0068] Step 4-3: Calculate and store each abundance column vector a of abundance matrix A. :j Corresponding spatial neighborhood coordinate family The formula is:

[0069]

[0070] Each column of the abundance matrix corresponds to a spectral band. It is a set of spatial coordinates (X,Y) that satisfy the above formula conditions, where T∈(0,1) represents the threshold multiple for determining the number of neighborhoods. Let a be the vector in the j-th column of the abundance matrix A. :j The corresponding pixel spatial spectrum similarity matrix H j ∈R r×c The element in the Xth row and Yth column, max(H) j ) is matrix H j The maximum element value;

[0071] in, The pixel spatial spectral similarity S(y) between pixel a and pixel b a ,y b ) is defined as σ1 is a parameter controlling pixel spectral similarity, set to 1; σ2 is a parameter controlling distance significance, set to...

[0072] Step 4-4, calculate the spectral weighting matrix W spe Calculation formula:

[0073]

[0074] in, This represents the i-th row of the abundance matrix in the k-th iteration, ε = 1 × 10 -6 The constant ε is used to avoid numerical instability and takes a value greater than 0;

[0075] Calculate the spatial weighting matrix W spa The calculation formula is:

[0076]

[0077] Where f(·) is a function characterizing spatial correlation. The element a in the i-th row and j-th column of the abundance matrix represents... ij The neighborhood set of (i = 1, ..., n; j = 1, ..., N), in, Indicates a ij and The neighborhood saliency;

[0078] Steps 4-5: Based on the weighted residual mechanism obtained in steps 2-6, a higher-order graph Laplacian regularization term is introduced to obtain a robust non-negative matrix factorization model based on global structural constraints, as follows:

[0079]

[0080] Where λ1, λ2, and λ3 are the first, second, and third regularization parameters, respectively, and all of them take values ​​greater than 0. a j Let E be the j-th column of the abundance matrix A, ⊙ denotes the Hadamard product, F1(E) is the endmember smoothing constraint term, and E i,j Let E represent the element in the i-th row and j-th column of the endmember matrix E. For larger and more detailed hyperspectral images, by adding a higher-order graph Laplacian regularization term, the abundance similarity of potentially related pixels at greater distances can be enhanced. At the same time, combined with the smoothing effect of adaptive spatial-spectral constraints on adjacent abundances, the global structure of the hyperspectral image can be preserved to a great extent.

[0081] Step 5 specifically involves:

[0082] The optimization of the robust nonnegative matrix factorization model based on global structural constraints is divided into two sub-tasks:

[0083]

[0084] The corresponding optimization strategy includes two steps: endmember estimation and abundance estimation. Specifically, the abundance A is fixed first, the endmember E is optimized, and after the endmember E is optimized, the abundance A is optimized again.

[0085] The endmember estimation specifically refers to:

[0086] Step 5-1-1, define the optimization problem for endmember E as:

[0087]

[0088] Where Ψ∈R L×n Represents the Lagrange multipliers;

[0089] Step 5-1-2: Find the partial derivative of J(E) with respect to the endmember E, and combine it with the Karush-Kuhn-Tucker (KKT) conditions to obtain:

[0090] E.*Ψ=0

[0091]

[0092] in, Describes the gradient of matrix E. It is the partial derivative of F1(E) with respect to E;

[0093] Step 5-1-3, the formula described in step 5-1-2 Multiplying both sides of the matrix E by the matrix, and combining this with E.*Ψ=0, we derive the update formula for E as follows:

[0094]

[0095] In this context, .* and . / represent element-wise multiplication and division of the matrix, respectively.

[0096] Step 5-1-4: Iteratively update the endmember matrix E according to the update formula described in step 5-1-3.

[0097] The abundance estimation specifically refers to:

[0098] Step 5-2-1: For the optimization problem of abundance A, the original matrix is ​​replaced by an augmented form of the hyperspectral image matrix and the endmember matrix to satisfy the abundance sum-to-one constraint (ASC), according to the formula... and Augmenting Y, E and W respectively N Obtain the augmented form of Y Augmented form of E and W N augmented forms in, yes The weights in the i-th row and i-th column are δ and β, which are parameters used to control the degree of influence of ASC, where δ>0 and β>0.

[0099] Step 5-2-2, define the optimization problem of A as:

[0100]

[0101] Step 5-2-3: Calculate the partial derivative of the abundance matrix A, and using the KKT conditions, obtain:

[0102] A.*Φ=0

[0103]

[0104] Among them, represents the gradient of the abundance matrix A;

[0105] Step 5-2-4, multiply both sides of the formula described in Step 5-2-3 by the matrix A, and the updated formula for the abundance matrix A is derived as:

[0106]

[0107] Step 5-2-5, iteratively update the abundance A according to the update formula described in Step 5-2-4;

[0108] Step 5-2-6, if the number of iterations reaches the set value or satisfies:

[0109] |F (k+1) (E, A) - F (k) (E, A)| < TOL, then output the endmember E and the abundance A, otherwise update the diagonal weight matrix W according to Step 2-5 N , and then return to Step 4-4, where TOL is the stopping threshold, and the value of the objective function is reduced through iteration to achieve the final optimization effect.

[0110] Beneficial effects:

[0111] The robust non - negative matrix factorization hyperspectral image unmixing method based on global constraints in the present invention can obtain a more robust and accurate unmixing effect compared with the problem that the standard NMF model has non - convexity and is vulnerable to noise when applied to hyperspectral image unmixing, resulting in poor unmixing performance. Moreover, it can better protect the internal structure of hyperspectral images and is more suitable for scenarios where the ground objects are more complex, the details are richer, or the requirement for unmixing accuracy is higher. Brief description of the drawings

[0112] Figure 1 is the algorithm flow chart of the present invention.

[0113] Figure 2 is the true abundance map of trees, water bodies, soil, and roads.

[0114] Figure 3 is the processed abundance map. Detailed implementation manners

[0115] Traditional nonnegative matrix factorization (NMF) models, which use Euclidean distance as the loss function metric, are susceptible to noise and outliers. Furthermore, existing unmixing algorithms do not fully utilize the spatial structure information of hyperspectral images (HSI), often failing to achieve expected performance when handling highly mixed and complex HSI environments. To address these issues, this invention proposes a robust nonnegative matrix factorization hyperspectral image unmixing method based on global structural constraints.

[0116] like Figure 2 As shown, the true abundance maps of trees, water bodies, soil, and roads are presented in sequence. This embodiment provides a robust non-negative matrix factorization hyperspectral image demixing method, including:

[0117] Step 1 is as follows:

[0118] Step 1-1, assuming the acquired hyperspectral image data containing L bands and N pixels is as follows:

[0119] Y = [y1,...,y] N ]∈R L×N

[0120] Define matrix Z = Y T =[z1,z2,…,z L ]∈R N×L , where vector y i Let z represent the spectrum of the i-th pixel. i This represents the reflection value of the i-th band in Y;

[0121] Define an N×(L-1) dimensional matrix Z di =[z1,z2,…,z i-1 ,z i+1 ,…z L ], Z di Let Z be the matrix after removing the i-th column, (·) T Represents the transpose of a matrix;

[0122] Steps 1-2: Based on the correlation between hyperspectral image bands, the reflectance value z i The remaining L-1 vectors can be linearly represented as z. i =Z di β i +n i , where β i ∈R L-1 Let n represent the regression vector. i Represents the error vector;

[0123] Steps 1-3, the minimum mean square estimate of the regression vector is Where i = 1, ..., L, the estimated noise value is obtained from this. The noise estimation matrix is ​​then...

[0124] Step 2 is as follows:

[0125] Step 2-1: The correlation matrix of the noise is calculated using the noise estimation matrix.

[0126]

[0127] Where N represents the number of pixels in the hyperspectral image;

[0128] Step 2-2, the noise power of each band is P N =[C 11 C 22 ,…,C LL ], where C ii This represents the noise power of the i-th band;

[0129] Steps 2-3: Due to the randomness and uncertainty of noise, the noise power in different bands may vary significantly. To effectively allocate the weight of each band, the noise power is normalized as follows:

[0130]

[0131] The band weights are then obtained as follows:

[0132]

[0133] Among them, b i Let α represent the weight of the i-th band, i = 1, ..., L, and α be the scaling parameter. Considering the rationality of the weight, α is set to 8. Through the above strategy, the weight of the band can be adjusted in the range of (0,1). The weight of the band with more severe noise interference is closer to 1, so that the model will focus on the band in subsequent algorithms.

[0134] Steps 2-4 involve obtaining an initial estimate of the endmember matrix E of the hyperspectral image using a vertex component analysis algorithm. After successful extraction of the endmember matrix E, the initial value of the abundance matrix A is obtained by minimizing the error between the hyperspectral image data Y and the spectrum reconstructed by a linear mixture model based on endmember spectra and abundance estimates using the fully constrained least squares method. The fully constrained least squares method can simultaneously satisfy the constraints of non-negativity of abundance and the sum of abundance being one, where E∈R. L×n An endmember matrix is ​​a matrix that represents the unique spectral characteristics of a ground feature or substance, where A ∈ R. n×N The abundance matrix is ​​a matrix that represents the spectral band information in an image, where n represents the number of endmembers.

[0135] Steps 2-5 introduce noise into the weighted residual NMF model to assign the initial bands of the hyperspectral image by measuring the degree of noise interference.

[0136] The residual weight matrix of the weighted residual NMF model improved by noise estimation is as follows:

[0137] W N =diag[w N,11 ,w N,22 ,…,w N,LL ]

[0138] Define its diagonal elements as:

[0139]

[0140] Among them, w N,ii Let R represent the weight of the i-th band, k represent the number of iterations, k = 0 represent the initial state, ||·||² represent the Euclidean norm, and R is the weight of the i-th band. (k) =YE (k) A (k) R represents the residual of the k-th iteration. i,: Let R represent the i-th row, where μ is a constant characterizing the degree of exponential decay of the residual, and μ = 30.

[0141] Steps 2-6 introduce the improved weighted residual NMF model based on noise estimation:

[0142]

[0143] Among them, ||·|| F Describing the Frobenius norm, and Let N and N represent vectors consisting entirely of 1s, respectively.

[0144] The improved weighting mechanism provides initial weights for each band before iteration, and combined with dynamic adjustments during the iteration process, it can effectively reduce the impact of noise and improve the robustness and unmixing efficiency of the model.

[0145] Step 3 specifically involves:

[0146] Step 3-1: Given a hyperspectral image matrix Y, construct adjacency relationships using pixels as vertices. The first-order adjacency weights are expressed as:

[0147]

[0148] Among them, y i Let y represent the i-th pixel. j Let θ represent the j-th pixel, and θ is a parameter that controls the degree of weight variation.

[0149] Step 3-2: To more effectively learn the global structure of the hyperspectral image, the weights of the first-order adjacency matrix are improved by combining pixel spatial distance, making the algorithm pay more attention to pixels outside the local neighborhood. The improved first-order adjacency weights are:

[0150]

[0151] Where d(y) i ,y j θ1 and θ2 represent the spatial Euclidean distance between two pixels. θ1 and θ2 are parameters controlling the pixel saliency and the degree of change in spatial distance, respectively. We set θ1 = 1. For parameter θ2, to measure the influence of pixel distance, assuming the spatial size of the hyperspectral image is r × c, we set θ2 = (r × c) / (r × c) × c. 2 +c 2 The improved first-order adjacency weight assigns greater weight to similar pixels that are farther apart, thus representing a closer relationship between the two pixels.

[0152] Step 3-3, define the first-order Laplace matrix as:

[0153] H1 = D1 - W1

[0154] Where D1 is the degree matrix, W1 is the improved first-order adjacency weight matrix;

[0155] Steps 3-4 define the second-order adjacency weight as follows:

[0156]

[0157] Among them, W 1,i Represents the i-th row of W1, W 1,j In the j-th row of W1, when two pixels have a large number of interconnected adjacent pixels, W... 2,ij A relatively large value for W indicates a significant second-order adjacency relationship between the two pixels. When there are no mutually connected adjacent pixels, then W... 2,ij A value of zero indicates that the two pixels do not have a second-order adjacency relationship;

[0158] Further derivation of the second-order Laplace matrix is ​​as follows:

[0159] H2 = D2 - W2

[0160] Where D2 is a diagonal matrix. W2 is a second-order adjacency weight matrix;

[0161] Steps 3-5: Combining first-order and second-order adjacency relations (represented by first-order and second-order Laplace matrices respectively), the higher-order Laplace matrix H... h Diagonal matrix D hand higher-order adjacency weight matrix W h They are defined as follows:

[0162] H h =ρH1+(1-ρ)H2

[0163] D h =ρD1+(1-ρ)D2

[0164] W h =ρW1+(1-ρ)W2

[0165] Where ρ is a parameter that balances the impact of first-order and second-order adjacency relations on overall performance, and is taken as ρ = 0.4;

[0166] Steps 3-6, based on the higher-order Laplace matrix H h Diagonal matrix D h and higher-order adjacency weight matrix W h The abundance vectors a corresponding to the i-th and j-th pixels :i and a :j The corresponding minimized regularization term is as follows:

[0167]

[0168] Among them, W h,ij W represents the higher-order adjacency weight matrix h The weight value in the i-th row and j-th column, D h,ii Represents a diagonal matrix D h The weight value in the i-th row and i-th column, Tr(·) represents the trace of the matrix.

[0169] Step 4 is as follows:

[0170] Step 4-1: Calculate the estimated endmember matrix E using an unconstrained NMF model. est And estimate the spectrum, where, in the traditional NMF model, both the endmember matrix and the abundance matrix must satisfy non-negativity constraints, and the unconstrained NMF model refers to the NMF model variant that does not impose endmember and abundance non-negativity constraints;

[0171] Step 4-2: Calculate the endmember smoothing weight matrix Q, and the weight Q of the endmember smoothing level parameter corresponding to the element in the i-th row and j-th column. ij The formula is:

[0172]

[0173] Where η is a parameter controlling the smoothness, and η = 1, E est；i,j Indicates the estimated endmember matrix E est The element in the i-th row and j-th column is the weight Q of the endmember smoothing level parameter. ijE is measured by the difference between the estimated spectrum and the mean of its two preceding and following reflectance values. est；i+1,j The closer to (E) est；i+2,j +E est；i,j The smoother the surface at this point, the higher the weight Q. ij The larger it is;

[0174] Step 4-3: Calculate and store each abundance column vector a of abundance matrix A. :j Corresponding spatial neighborhood coordinate family The formula is:

[0175]

[0176] Each column of the abundance matrix corresponds to a spectral band. It is a set of spatial coordinates (X,Y) that satisfy the conditions, with T = 0.7, where T represents the threshold multiple for determining the number of neighbors, used to scale max(H). j The final threshold is determined by T and max(H). j The product of ) is determined, and the coordinates of elements greater than the threshold are denoted as a. :j The neighborhood spatial coordinates, element a in the i-th row and j-th column of the abundance matrix. ij The neighborhood set is represented as Let a be the vector in the j-th column of the abundance matrix A. :j The corresponding pixel spatial spectrum similarity matrix H j ∈R r×c The element in the Xth row and Yth column, max(H) j ) is matrix H j The maximum element value, where, The pixel spatial spectral similarity S(y) between pixel a and pixel b a ,y b ) is defined as σ1 is a parameter controlling the spectral similarity of pixels, set to σ1 = 1; σ2 is a parameter controlling the significance of distance, set to...

[0177] Step 4-4, calculate the spectral weighting matrix W spe Calculation formula:

[0178]

[0179] in, This represents the i-th row of the abundance matrix in the k-th iteration, ε = 1 × 10 -6 This is used to avoid numerical instability, according to the formula. Calculate the spatial weighting matrix W spa ,

[0180] Where f(·) is a function characterizing spatial correlation. in, Indicates a ij and The neighborhood saliency;

[0181] Steps 4-5 improve the original weighted residual mechanism through noise estimation and introduce a higher-order graph Laplacian regularization term to obtain a robust NMF model based on global structural constraints as follows:

[0182]

[0183] Where λ1, λ2, and λ3 are the first, second, and third regularization parameters, respectively, and are set to λ1 = 0.1, λ2 = 0.1, and λ3 = 0.1. Q' represents the mean of the normalized matrix Q' of the endmember smoothing weight matrix Q. ij Defined as a j Let E be the j-th column of the abundance matrix A, ⊙ denotes the Hadamard product, F1(E) is the endmember smoothing constraint term, and E i,j Let E represent the element in the i-th row and j-th column of the endmember matrix E. For larger and more detailed hyperspectral images, by adding a higher-order graph Laplacian regularization term, the abundance similarity of potentially related pixels at greater distances can be enhanced. At the same time, combined with the smoothing effect of adaptive spatial-spectral constraints on adjacent abundances, the global structure of the hyperspectral image can be preserved to a great extent.

[0184] Step 5 specifically involves:

[0185] The optimization of the robust nonnegative matrix factorization model based on global structural constraints is divided into two sub-tasks:

[0186]

[0187] The corresponding optimization strategy includes two steps: endmember estimation and abundance estimation. Specifically, the abundance A is fixed first and the endmember E is optimized. After optimizing the endmember E, the abundance A is optimized again. Through iteration, the value of the objective function is reduced, and the final optimization effect is achieved.

[0188] Endmember estimation is specifically as follows:

[0189] Step 5-1-1, define the optimization problem for endmember E as:

[0190]

[0191] Where Ψ∈R L×n Represents the Lagrange multipliers;

[0192] Step 5-1-2: Find the partial derivative of J(E) with respect to the endmember E, and combine it with the Karush-Kuhn-Tucker (KKT) conditions to obtain:

[0193] E.*Ψ=0

[0194]

[0195] in, Describes the gradient of matrix E. It is the partial derivative of F1(E) with respect to E;

[0196] Step 5-1-3, the formula described in step 5-1-2 Multiplying both sides of the matrix E by the matrix, and combining this with E.*Ψ=0, we derive the update formula for E as follows:

[0197]

[0198] In this context, .* and . / represent element-wise multiplication and division of the matrix, respectively.

[0199] Step 5-1-4: Iteratively update the endmember matrix E according to the update formula described in step 5-1-3.

[0200] The abundance estimation is as follows:

[0201] Step 5-2-1: For the optimization problem of abundance A, the original matrix is ​​replaced by an augmented form of the hyperspectral image matrix and the endmember matrix to satisfy the abundance sum-to-one constraint (ASC), according to the formula... and Augmenting Y, E and W respectively N Obtain the augmented form of Y Augmented form of E and W N augmented forms in, yes The weight of the i-th row and i-th column, δ and β are parameters used to control the degree of influence of ASC, with δ = 15 and β = 0.5.

[0202] Step 5-2-2, define the optimization problem of A as:

[0203]

[0204] Step 5-2-3: Calculate the partial derivative of the abundance matrix A, and using the KKT conditions, obtain:

[0205] A.*Φ=0

[0206]

[0207] wherein, represents the gradient of the abundance matrix A;

[0208] Step 5-2-4, multiply both sides of the formula described in Step 5-2-3 by the matrix A simultaneously, and the updated formula for the abundance matrix A is derived as:

[0209]

[0210] Step 5-2-5, iteratively update the abundance A according to the updated formula described in Step 5-2-4;

[0211] Step 5-2-6, if the number of iterations reaches the set value or satisfies:

[0212] |F (k+1) (E,A)-F (k) (E,A)| < TOL, then output the endmember E and the abundance A, otherwise update the diagonal weight matrix W according to Step 2-5 N , then return to Step 4-4, set the stop threshold TOL = 1×10 -3 , and the maximum number of iterations is 1000. By using a robust non - negative matrix factorization hyperspectral image unmixing method proposed by the present invention to process the hyperspectral image, the processed abundance map is obtained, as Figure 3 shown, indicating that the present invention can effectively improve the unmixing performance while protecting the global structure of the hyperspectral image.

[0213] The present invention provides a robust non - negative matrix factorization hyperspectral image unmixing method. There are many methods and ways to specifically implement this technical solution. The above description is only the preferred embodiment of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention. Each component not clearly defined in this embodiment can be implemented by existing technologies.

Claims

1. A robust nonnegative matrix factorization method for hyperspectral image unmixing, characterized in that, The steps include the following: Step 1: Calculate the noise estimate and noise estimation matrix of the hyperspectral image; Step 2: Construct an improved weighted residual nonnegative matrix decomposition model based on the noise estimate and the noise estimation matrix; Step 3: Introduce second-order adjacency relations to obtain a higher-order graph and perform Laplace regularization to obtain regularization terms; Step 4: Introduce the higher-order graph Laplace regularization term to obtain a robust non-negative matrix decomposition model based on global structural constraints; Step 5: Implement hyperspectral image demixing using a robust nonnegative matrix factorization model based on global structural constraints.

2. The robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 1, characterized in that, Step 1 is as follows: Step 1-1, assuming the acquired hyperspectral image data containing L bands and N pixels is as follows: Y=[y1,…,y N ]∈R L×N Define matrix Z = Y T = [z1, z2, ..., z L ]∈R N×L , where vector y i Let z represent the spectrum of the i-th pixel. i This represents the reflection value of the i-th band in Y; Define an N×(L-1) dimensional matrix Z di = [z1, z2, ..., z i-1 , z i+1 , ...z L ], Z di Let Z be the matrix after removing the i-th column, (·) T Represents the transpose of a matrix; Steps 1-2: Based on the correlation between hyperspectral image bands, the reflectance value z i The remaining L-1 vectors can be linearly represented as z. i =Z di β i +n i , where β i ∈R L-1 Let n represent the regression vector. i Represents the error vector; Steps 1-3, the minimum mean square estimate of the regression vector is Where i = 1, ..., L, the estimated noise value is obtained from this. The noise estimation matrix is ​​then...

3. The robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 1, characterized in that, Step 2 is as follows: Step 2-1: The correlation matrix of the noise is calculated using the noise estimation matrix. Where N represents the number of pixels in the hyperspectral image; Step 2-2, the noise power of each band is P N =[C 11 C 22 C LL ], where C ii This represents the noise power of the i-th band; Steps 2-3, normalize the noise power as follows: The band weights are then obtained as follows: Among them, b i This represents the weight of the i-th band, where i = 1, ..., L, and α is the scaling parameter. Steps 2-4: Initial estimates of the endmember matrix E of the hyperspectral image are obtained by performing vertex component analysis on the hyperspectral image. After the endmember matrix E is successfully extracted, the initial value of the abundance matrix A is obtained by minimizing the error between the hyperspectral image data Y and the spectrum EA reconstructed by the linear mixture model based on the endmember spectrum and abundance estimation using the fully constrained least squares method. Steps 2-5 introduce noise into the weighted residual nonnegative matrix factorization model and assign the initial bands of the hyperspectral image by measuring the degree of noise interference. The residual weight matrix of the weighted residual nonnegative matrix factorization model improved by noise estimation is: IN N =diag[w N,11 ,In N,22 ,…,In N,LL ] Steps 2-6 introduce the improved weighted residual nonnegative matrix decomposition model based on noise estimation: Among them, ||·|| F Describing the Frobenius norm, and Let N and N represent vectors consisting entirely of 1s, respectively.

4. The robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 1, characterized in that, Step 3 specifically involves: Step 3-1: Given a hyperspectral image matrix Y, construct adjacency relationships using pixels as vertices. The first-order adjacency weights are expressed as: Among them, y i Let y represent the i-th pixel. j Let θ represent the j-th pixel, and θ is a parameter that controls the degree of weight change. θ takes a non-zero real number. Step 3-2: Improve the first-order adjacency matrix weights by combining the spatial distance of pixels. The improved first-order adjacency weights are as follows: Where d(y) i y j θ1 and θ2 represent the spatial Euclidean distance between two pixels, respectively, and are parameters controlling pixel saliency and the degree of change in spatial distance. Assuming the spatial size of the hyperspectral image is r×c, then θ2 is set to (r × c). 2 +c 2 ) / (e-1), where e is a constant; Step 3-3, define the first-order Laplace matrix as: H1 = D1 - W1 Where D1 is the degree matrix, W1 is the improved first-order adjacency weight matrix; Steps 3-4 define the second-order adjacency weight as follows: Among them, W 1，i Represents the i-th row of W1, W 1，j In the j-th row of W1, when two pixels have a large number of interconnected adjacent pixels, W... 2，ij A value greater than 0 indicates a significant second-order adjacency relationship between the two pixels. When there are no mutually connected adjacent pixels, then W... 2，ij A value of zero indicates that the two pixels do not have a second-order adjacency relationship; Further derivation of the second-order Laplace matrix is ​​as follows: H2 = D2 - W2 Where D2 is a diagonal matrix. W2 is a second-order adjacency weight matrix; Steps 3-5: Combining first-order and second-order adjacency relations, the higher-order Laplace matrix H is obtained. h Diagonal matrix D h and higher-order adjacency weight matrix W h They are defined as follows: H h =ρH1+(1-ρ)H2 D h =ρD1+(1-ρ)D2 W h =ρW1+(1-ρ)W2 Wherein, ρ is a parameter that balances the influence of first-order and second-order adjacency relations on the overall performance, and its value ranges from [0, 1]. Steps 3-6, based on the higher-order Laplace matrix H h Diagonal matrix D h and higher-order adjacency weight matrix W h The abundance vectors a corresponding to the i-th and j-th pixels ：i and a ：j The corresponding minimized regularization term is as follows: Among them, W h，ij W represents the higher-order adjacency weight matrix h The weight value in the i-th row and j-th column, D h，ii Represents a diagonal matrix D h The weight value in the i-th row and i-th column, Tr(·) represents the trace of the matrix.

5. The robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 1, characterized in that, Step 4 is as follows: Step 4-1: Calculate the estimated endmember matrix E using an unconstrained nonnegative matrix factorization model. est And estimate the spectrum; Step 4-2: Calculate the endmember smoothing weight matrix Q, and the weight Q of the endmember smoothing level parameter corresponding to the element in the i-th row and j-th column. ij The formula is: Where η is a parameter controlling the smoothness, and η takes the form of a non-zero real number, E est；i，j Indicates the estimated endmember matrix E est The element in the i-th row and j-th column is the weight Q of the endmember smoothing level parameter. ij E is measured by the difference between the estimated spectrum and the mean of its two preceding and following reflectance values. est；i+1，j The closer to (E) est；i+2，j +E est；i，j The smoother the surface at this point, the higher the weight Q. ij The larger it is; Step 4-3: Calculate and store each abundance column vector a of abundance matrix A. ：j Corresponding spatial neighborhood coordinate family The formula is: Each column of the abundance matrix corresponds to a spectral band. It is a set of spatial coordinates (X, Y) that satisfy the above formula conditions, where T∈(0,1) represents the threshold multiple for determining the number of neighborhoods. Let a be the vector in the j-th column of the abundance matrix A. ：j The corresponding pixel spatial spectrum similarity matrix H j ∈R r×c The element in the Xth row and Yth column, max(H) j ) is matrix H j The maximum element value; in, The pixel spatial spectral similarity S(y) between pixel a and pixel b a y b ) is defined as σ1 is a parameter controlling pixel spectral similarity, and σ2 is a parameter controlling distance significance. (Settings are missing from the original text.) Step 4-4, calculate the spectral weighting matrix W spe Calculation formula: in, This represents the i-th row of the abundance matrix in the k-th iteration. The constant ε is used to avoid numerical instability and takes a value greater than 0. Calculate the spatial weighting matrix W spa The calculation formula is: Where f(·) is a function characterizing spatial correlation. The element a in the i-th row and j-th column of the abundance matrix represents... ij The neighborhood set of (i = 1, ..., n; j = 1, ..., N). in, Indicates a ij and The neighborhood saliency; Steps 4-5: Based on the weighted residual mechanism obtained in steps 2-6, a higher-order graph Laplacian regularization term is introduced to obtain a robust non-negative matrix factorization model based on global structural constraints, as follows: Where λ1, λ2, and λ3 are the first, second, and third regularization parameters, respectively, and all of them take values ​​greater than 0. a j Let E be the j-th column of the abundance matrix A, ⊙ denotes the Hadamard product, F1(E) is the endmember smoothing constraint term, and E i，j This represents the element in the i-th row and j-th column of the endmember matrix E.

6. The robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 1, characterized in that, Step 5 specifically involves: The optimization of the robust nonnegative matrix factorization model based on global structural constraints is divided into two sub-tasks: The corresponding optimization strategy includes two steps: endmember estimation and abundance estimation. Specifically, the abundance A is fixed first, the endmember E is optimized, and after the endmember E is optimized, the abundance A is optimized again.

7. A robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 6, characterized in that, The endmember estimation specifically refers to: Step 5-1-1, define the optimization problem for endmember E as: Where Ψ∈R L×n Represents the Lagrange multipliers; Step 5-1-2: Find the partial derivative of J(E) with respect to the endmember E, and combine it with the KKT conditions to obtain: E.*Ψ=0 in, Describes the gradient of matrix E. It is the partial derivative of F1(E) with respect to E; Step 5-1-3, the formula described in step 5-1-2 Multiplying both sides of the matrix E by the matrix, and combining this with E.*Ψ=0, we derive the update formula for E as follows: In this context, .* and . / represent element-wise multiplication and division of the matrix, respectively. Step 5-1-4: Iteratively update the endmember matrix E according to the update formula described in step 5-1-3.

8. A robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 7, characterized in that, The abundance estimation specifically refers to: Step 5-2-1, according to the formula and Augmenting Y, E and W respectively N Obtain the augmented form of Y Augmented form of E and W N augmented forms in, yes The weights in the i-th row and i-th column are δ and β, which are parameters used to control the degree of influence of ASC, where δ > 0 and β > 0. Step 5-2-2, define the optimization problem of A as: Step 5-2-3: Calculate the partial derivative of the abundance matrix A, and using the KKT conditions, obtain: A.*Φ=0 in, This represents the gradient of the abundance matrix A; Step 5-2-4, the formula described in step 5-2-3 By multiplying both sides of matrix A, the update formula for the abundance matrix A is derived as follows: Step 5-2-5: Iteratively update abundance A according to the update formula described in step 5-2-4; Step 5-2-6, if the number of iterations reaches the set value or the following conditions are met: |F (k+1) (E, A)-F (k) If (E, A)| < TOL, then output the endmember E and abundance A; otherwise, update the diagonal weight matrix W according to steps 2-5. N Then return to step 4-4, where TOL is the stopping threshold. By iterating, the value of the objective function is reduced to achieve the final optimization effect.

9. A robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 3, characterized in that, In steps 2-4, E∈R L×n An endmember matrix is ​​a matrix that represents the unique spectral characteristics of a ground feature or substance, where A ∈ R. n×N The abundance matrix is ​​a matrix that represents the spectral band information in an image, where n represents the number of endmembers.

10. A robust nonnegative matrix factorization hyperspectral image unmixing method according to claim 3, characterized in that, The residual weight matrix W described in steps 2-5 N Its diagonal elements are defined as: Among them, w N，ii Let R represent the weight of the i-th band, k represent the number of iterations, k = 0 represent the initial state, ||·||² represent the Euclidean norm, and R is the weight of the i-th band. (k) =YE (k) A (k) R represents the residual of the k-th iteration. i，： Let R represent the i-th row, where μ is a constant characterizing the degree of residual exponential decay.

Images