CN113723348A - Hyperspectral mixed pixel decomposition method based on abundance sparsity constraint - Google Patents

Abstract

In order to improve the precision of hyperspectral mixed pixel decomposition and solve the problem that abundance sparsity is neglected in a PCHA (principal convex hull prototype analysis) algorithm, a hyperspectral mixed pixel decomposition method based on abundance sparsity constraint is provided. The method combines the matrix decomposition form of the original algorithm, and adds l of abundance on the basis of the original objective function₁-norm sparsity-constraining regularization term; then, the coefficient matrix and the abundance matrix are optimized and solved by using a projection gradient and an alternating direction multiplier method respectively; and finally, obtaining end members and abundance according to the solving result. The method can simultaneously realize the extraction of the end members and the estimation of the abundance, has lower data fitting error, and the solved abundance is closer to the actual situation in the physical sense.

Description

Hyperspectral mixed pixel decomposition method based on abundance sparsity constraint

Technical Field

The invention belongs to the technical field of hyperspectral mixed pixel decomposition, and particularly relates to a hyperspectral mixed pixel decomposition algorithm based on abundance sparsity constraint.

Background

Due to the limited spatial resolution of the hyperspectral image, complex and various natural features and the mixed effect in the image acquisition process, mixed pixels widely exist in the hyperspectral remote sensing image, the hyperspectral classification precision is low, the effect is poor, and the development and the application of the hyperspectral remote sensing technology are restricted. End member extraction and abundance estimation in the spectrum unmixing technology are key technologies for solving the problem of mixed pixels. Wherein the end-member extraction can identify potential characteristics in the image, extract the ground object classes forming the image, and the abundance estimation can obtain the components of each ground object in each pixel.

The traditional mixed pixel decomposition algorithm extracts end members firstly and then estimates the abundance, so that the problem of accumulation of a large amount of errors exists, and the accuracy of abundance estimation depends heavily on the accuracy of end member extraction. The blind decomposition method utilizes the inherent physical characteristics of the end members and the abundance, such as nonnegativity of the end members, nonnegativity of the abundance, 1 sum of the abundance, sparsity of the abundance and the like, can simultaneously obtain the spectral information of the end members and the component information of the end members in each pixel under the condition of not giving the information of the end members, and is more effective compared with the traditional method. The prototype analysis algorithm is one of the blind decomposition algorithms.

Aiming at the problem that abundance sparsity is neglected in a PCHA algorithm of main convex hull prototype analysis, a hyperspectral mixed pixel decomposition method based on abundance sparsity constraint is provided. The method adds l of abundance on the basis of the target function of the original algorithm₁-norm sparsity-constraining regularization term; then, carrying out optimization solution on the objective function by using a projection gradient and alternating direction multiplier method; and finally, obtaining end members and abundance according to the solving result. The method considers the sparse physical significance of abundance and improves the solving precision of the original algorithm while ensuring that the data fitting has the minimum error.

Disclosure of Invention

Technical problem to be solved

The invention aims to further optimize a hyperspectral mixed pixel decomposition method based on a main convex hull prototype analysis algorithm, and the main problem to be solved is to construct an objective function based on abundance sparsity constraint and optimize and solve the objective function.

(II) technical scheme

In order to achieve the purpose, the invention adopts the following technical scheme: a hyperspectral mixed pixel decomposition method based on abundance sparsity constraint. The invention comprises the following steps:

(1) inputting a hyperspectral image R with a size of LxMxN_mixedConverting it into a two-dimensional Lxn shadowAn image matrix R, where L is the total number of bands, n is the total number of pixels, and its size is:

n＝M×N

wherein M is an image R_mixedN is the number of rows of the image;

(2) constructing an objective function based on abundance sparsity constraint main convex hull prototype analysis:

s.t.C≥0,S≥0，1-δ≤α_i≤1+δ,|c_i|₁＝1,|s_j|₁＝1

wherein | · | purple_FFrobenius norm representing the solution matrix, δ being the relaxation factor, α_iIs a scaling factor, alpha is a scaling factor alpha_iThe n-dimensional vector is formed, s.t. represents a constraint condition, λ is an abundance sparse regular term parameter, i is 1,2 …, n, j is 1,2, … p, p is the number of end members, diag is a diagonal function, and is used for extracting diagonal elements of the matrix and creating a diagonal matrix:

(3) initializing a scaling factor alpha₀Initializing coefficient matrix C₀Sum abundance matrix S₀Setting an initial scaling factor alpha₀0, auxiliary abundance sparsity term regularization term parameter λ₀0.1, scaling factor regularization term parameter μ_αCoefficient matrix regularization term parameter μ 1_C1 and the abundance sparsity term assists the regularization parameter mu_SA value of 10 λ, where λ is expressed as an abundance sparsity term canonical term parameter, defined as:

obtaining initialized coefficient matrix C by utilizing FURTHESTSUM initialized coefficient matrix C₀According to abundanceObtaining an initialized abundance matrix S by iterative solution expression of the degree matrix S₀：

Wherein Supdate represents that the iteration update rule of the abundance matrix S is that the iteration times are 10 times:

wherein k | · 1,2, … n, | | | | ceiling₁And |. non chlorine₁All represent solving matrix 1 norm, u, d are auxiliary sequences in the optimization process, and initial u₀，d₀Obtained by an initialization method in SUnSAL, and an end member matrix A is RC₀diag(α₀)，B＝A^TA+μ_SI, I is an identity matrix;

(4) optimizing an objective function, and solving a scaling factor vector alpha, a coefficient matrix C and an abundance matrix S; the objective function constructed above is decomposed into three subproblems:

in the formula, when argmin is expressed as a function and the minimum value is taken, an equation is established, k is the current iteration number, k is 1-10, a matrix of a subscript k represents a value of a matrix after k iterations, a matrix of a subscript k-1 represents a value of the matrix after k-1 iterations, the iteration number is 10, and a solution scaling factor vector alpha is obtained₁₀Coefficient matrix C ═ C₁₀And abundance matrix S ═ S₁₀；

Wherein, the iterative solving expression of the scaling factor vector alpha is as follows:

in the formula of alpha_k,iDenotes alpha_kValue of (i), mu_αIs a regular parameter to solve for the scaling factor,

is the coefficient after the normalization and is,

wherein, the iterative solving expression of the coefficient matrix C is as follows:

in the formula of_CIs the regular parameter for solving the coefficient matrix;

the iterative solution expression of the abundance matrix S is:

(5) substituting the initialized value in the step (3) into the iterative solution expression in the step (4) for iterative calculation, judging whether a cycle termination condition is reached, if not, continuing to execute the step (4), if so, terminating the cycle, and outputting a final coefficient matrix C and a dilute abundance matrix S; the condition of cycle termination is that the cycle frequency is more than 500 or the error R-RCdiag (alpha) S Y₂Less than 10^-6；

(6) Decomposing a hyperspectral image mixed pixel, and obtaining a final end member matrix A and an abundance matrix S according to a coefficient matrix C obtained by solving, wherein the end member matrix A has the calculation formula:

A＝RCdiag(α)。

(III) advantageous effects

The invention has the following advantages and positive effects: aiming at the defect that the original algorithm lacks consideration of the actual physical significance of abundance, the invention provides the hyperspectral mixed pixel decomposition method based on abundance sparsity.

Drawings

FIG. 1 is a schematic flow diagram of the process of the present invention.

Fig. 2 is a photographic image of the hyperspectral dataset Cuprite.

FIG. 3 is a comparison of the accuracy of each comparison algorithm with the method of the present invention.

Detailed Description

The key invention is a hyperspectral mixed pixel decomposition method based on abundance sparse prototype analysis, which mainly aims at improving the problem that the abundance sparse characteristic of a main convex hull prototype analysis algorithm is neglected, and improves the precision of hyperspectral mixed pixel decomposition. The present invention will be described in further detail with reference to fig. 1.

1. A hyperspectral mixed pixel decomposition method based on abundance sparsity constraint is characterized by comprising the following steps:

(1) inputting a hyperspectral image R with a size of LxMxN_mixedConverting the image into a two-dimensional L multiplied by n image matrix R, wherein L is the total wave band number, n is the total pixel number, and the size is as follows:

n＝M×N

wherein M is an image R_mixedN is the number of rows of the image;

(2) constructing an objective function based on abundance sparsity constraint main convex hull prototype analysis:

s.t.C≥0,S≥0，1-δ≤α_i≤1+δ,|c_i|₁＝1,|s_j|₁＝1

wherein | · | purple_FFrobenius norm representing the solution matrix, δ being the relaxation factor, α_iIs a scaling factor, alpha is a scaling factor alpha_iForming an n-dimensional vector, wherein s.t. represents a constraint condition, lambda is an abundance sparse regular term parameter, and i is equal to1,2 …, n, j ═ 1,2, … p, p is the number of end members, diag is a diagonal function, used for extraction of diagonal elements of the matrix and creation of diagonal matrix:

(3) initializing a scaling factor alpha₀Initializing coefficient matrix C₀Sum abundance matrix S₀Setting an initial scaling factor alpha₀0, auxiliary abundance sparsity term regularization term parameter λ₀0.1, scaling factor regularization term parameter μ_αCoefficient matrix regularization term parameter μ 1_C1 and the abundance sparsity term assists the regularization parameter mu_SA value of 10 λ, where λ is expressed as an abundance sparsity term canonical term parameter, defined as:

obtaining initialized coefficient matrix C by utilizing FURTHESTSUM initialized coefficient matrix C₀Obtaining an initialized abundance matrix S according to an iterative solution expression of the abundance matrix S₀：

Wherein Supdate represents that the iteration update rule of the abundance matrix S is that the iteration times are 10 times:

wherein k | · 1,2, … n, | | | | ceiling₁And |. non chlorine₁All represent solving matrix 1 norm, u, d are auxiliary sequences in the optimization process, and initial u₀，d₀Obtained by an initialization method in SUnSAL, and an end member matrix A is RC₀diag(α₀)，B＝A^TA+μ_SI, I is an identity matrix;

(4) optimizing an objective function, and solving a scaling factor vector alpha, a coefficient matrix C and an abundance matrix S; the objective function constructed above is decomposed into three subproblems:

in the formula, when argmin is expressed as a function and the minimum value is taken, an equation is established, k is the current iteration number, k is 1-10, a matrix of a subscript k represents a value of a matrix after k iterations, a matrix of a subscript k-1 represents a value of the matrix after k-1 iterations, the iteration number is 10, and a solution scaling factor vector alpha is obtained₁₀Coefficient matrix C ═ C₁₀And abundance matrix S ═ S₁₀；

Wherein, the iterative solving expression of the scaling factor vector alpha is as follows:

in the formula of alpha_k,iDenotes alpha_kValue of (i), mu_αIs a regular parameter to solve for the scaling factor,

is the coefficient after the normalization and is,

wherein, the iterative solving expression of the coefficient matrix C is as follows:

in the formula of_CIs the regular parameter for solving the coefficient matrix;

the iterative solution expression of the abundance matrix S is:

(5) substituting the initialized value in the step (3) into the iterative solution expression in the step (4) for iterative calculation, judging whether a cycle termination condition is reached, if not, continuing to execute the step (4), if so, terminating the cycle, and outputting a final coefficient matrix C and a dilute abundance matrix S; the condition of cycle termination is that the cycle frequency is more than 500 or the error R-RCdiag (alpha) S Y₂Less than 10^-6；

(6) Decomposing a hyperspectral image mixed pixel, and obtaining a final end member matrix A and an abundance matrix S according to a coefficient matrix C obtained by solving, wherein the end member matrix A has the calculation formula:

A＝RCdiag(α)

(7) and (3) precision evaluation, namely judging the extraction precision of the end member according to the spectral angular distance SAD of the extraction end member and the real end member, and judging the overall precision of image unmixing according to the image reconstruction error RMSE, wherein the smaller the spectral angular distance SAD and the image reconstruction error RMSE is, the better the extraction effect is.

The beneficial effects of the present invention are verified by comparative experiments as follows.

The data adopted in the experiment is a real data set, which is a Cuprice image with the size of 250 multiplied by 190 pixels and the wavelength of 188 wave bands obtained by AVIRIS, the spatial resolution is 20m, and the spectral resolution is 10 nm. The data contains a total of 12 surface feature types. The 12 ground objects are sequentially: limonite, ferromanganese, manganese iron ore, hard montmorillonite, kaolinite _1, kaolinite _2, muscovite, montmorillonite, nontronite, magnalium garnet, sphene and chalcedony. For data normalization consistent with the true reflectance values, the first band was removed, its wavelength was 390.09nm, and the total band number was 187. The red, green and blue wave bands in the true color image respectively correspond to the 28 th, 19 th and 10 th wave bands of the original image, the central wavelengths of the three wave bands are 635.72nm, 547.32nm and 458.89 nm respectively, and the data image is shown in figure 2. The peak component analysis and the full-constraint least square method VCA-FCLS, the minimum volume-constrained non-negative matrix decomposition algorithm MVC-NMF, the near convex prototype analysis algorithm NCAA, the main convex hull prototype analysis algorithm PCHA and the method are respectively adopted for spectrum unmixing, and the extracted results are quantitatively evaluated and compared, as shown in figure 3.

The spectral angular distance SAD of the invention is used for carrying out quantitative evaluation and analysis of end member extraction on the method and other spectral decomposition methods, the SAD value range is 0-pi, and the closer the value is to 0, the better the performance of the unmixing method is. The results are shown in Table 1.

TABLE 1 comparison of the accuracy of the process of the invention with other processes

From the experimental result, the end member extraction effect of the method is obviously superior to that of other comparison algorithms. The result of RMSE is not optimal, but the method realizes sparse extraction of abundance, and the estimation result of the abundance is more in line with practical significance, so the method still has promotion effect.

Claims (1)

1. A hyperspectral mixed pixel decomposition method based on abundance sparsity constraint is characterized by comprising the following steps:

(1) inputting a hyperspectral image R with a size of LxMxN_mixedConverting the image into a two-dimensional L multiplied by n image matrix R, wherein L is the total wave band number, n is the total pixel number, and the size is as follows:

n＝M×N

wherein M is an image R_mixedN is the number of rows of the image;

(2) constructing an objective function based on abundance sparsity constraint main convex hull prototype analysis:

s.t.C≥0,S≥0，1-δ≤α_i≤1+δ,|c_i|₁＝1,|s_j|₁＝1

wherein | · | purple_FFrobenius norm representing the solution matrixDelta is a relaxation factor, alpha_iIs a scaling factor, alpha is a scaling factor alpha_iThe n-dimensional vector is formed, s.t. represents a constraint condition, λ is an abundance sparse regular term parameter, i is 1,2 …, n, j is 1,2, … p, p is the number of end members, diag is a diagonal function, and is used for extracting diagonal elements of the matrix and creating a diagonal matrix:

(3) initializing a scaling factor alpha₀Initializing coefficient matrix C₀Sum abundance matrix S₀Setting an initial scaling factor alpha₀0, auxiliary abundance sparsity term regularization term parameter λ₀0.1, scaling factor regularization term parameter μ_αCoefficient matrix regularization term parameter μ 1_C1 and the abundance sparsity term assists the regularization parameter mu_SA value of 10 λ, where λ is expressed as an abundance sparsity term canonical term parameter, defined as:

obtaining initialized coefficient matrix C by utilizing FURTHESTSUM initialized coefficient matrix C₀Obtaining an initialized abundance matrix S according to an iterative solution expression of the abundance matrix S₀：

Wherein Supdate represents that the iteration update rule of the abundance matrix S is that the iteration times are 10 times:

wherein k | · 1,2, … n, | | | | ceiling₁And |. non chlorine₁All represent the solution matrix 1 normNumber u, d is an auxiliary sequence in the optimization process, initial u₀，d₀Obtained by an initialization method in SUnSAL, and an end member matrix A is RC₀diag(α₀)，B＝A^TA+μ_SI, I is an identity matrix;

(4) optimizing an objective function, and solving a scaling factor vector alpha, a coefficient matrix C and an abundance matrix S; the objective function constructed above is decomposed into three subproblems:

in the formula, when argmin is expressed as a function and the minimum value is taken, an equation is established, k is the current iteration number, k is 1-10, a matrix of a subscript k represents a value of a matrix after k iterations, a matrix of a subscript k-1 represents a value of the matrix after k-1 iterations, the iteration number is 10, and a solution scaling factor vector alpha is obtained₁₀Coefficient matrix C ═ C₁₀And abundance matrix S ═ S₁₀；

Wherein, the iterative solving expression of the scaling factor vector alpha is as follows:

in the formula of alpha_k,iDenotes alpha_kValue of (i), mu_αIs a regular parameter to solve for the scaling factor,

is the coefficient after the normalization and is,

wherein, the iterative solving expression of the coefficient matrix C is as follows:

in the formula of_CIs the regular parameter for solving the coefficient matrix;

the iterative solution expression of the abundance matrix S is:

(5) substituting the initialized value in the step (3) into the iterative solution expression in the step (4) for iterative calculation, judging whether a cycle termination condition is reached, if not, continuing to execute the step (4), if so, terminating the cycle, and outputting a final coefficient matrix C and a dilute abundance matrix S; the condition of cycle termination is that the cycle frequency is more than 500 or the error R-RCdiag (alpha) S Y₂Less than 10^-6；

(6) Decomposing a hyperspectral image mixed pixel, and obtaining a final end member matrix A and an abundance matrix S according to a coefficient matrix C obtained by solving, wherein the end member matrix A has the calculation formula:

A＝RCdiag(α)。

Images