CN113094645A - Hyperspectral data unmixing method based on morphological component constraint optimization - Google Patents

Abstract

The invention discloses a morphological component constraint optimization hyperspectral data unmixing method, which mainly utilizes a batch of hyperspectral data and a pure substance spectrum library to construct a spectrum sample expansion matrix; by estimating the band noise standard deviation and the sparsity structure noise decomposition, use l₁Norm to constrain the sparse noise,/_2,0The norm constrains the global row sparsity of abundance coefficients of different pure substances; and establishing a morphological component constraint optimized hyperspectral data unmixing model, and alternately iterating to realize linear unmixing of the mixed spectrum. The invention comprehensively considers Gaussian random noise and sparsityThe structural noise influences the linear unmixing precision, the mixed noise is robust, the waveform morphological structure difference between the same batch of spectral data can be effectively overcome, the rapid and high-precision unmixing is realized through optimization iteration, and the root mean square error of the unmixing is less than 0.0025; the method has wide application prospect for rock mineral identification, hyperspectral remote sensing ground object fine identification and the like.

Description

A hyperspectral data unmixing method for morphological component constrained optimization

technical field

The invention belongs to the technical field of hyperspectral remote sensing image processing, and in particular relates to a hyperspectral data unmixing method for morphological component constraint optimization.

Background technique

Remote sensing technology refers to the technology of obtaining electromagnetic signals reflected by ground objects from long-distance space (aerospace) or outer space (aviation) through optical sensors, obtaining the required target information and performing data analysis and processing. Since the 1980s, with the gradual improvement of the spectral resolution of optical sensors, hyperspectral remote sensing has gradually become a research hotspot in the field. Hyperspectral remote sensing refers to an optical remote sensing technology that uses sensors to obtain electromagnetic wave signals reflected from the surface of materials, and has a spectral resolution of 10 nm near the visible light to near-infrared bands. Due to the limitation of the spatial resolution of the sensor and the complex diversity of the surface, a pixel often contains a variety of ground objects, which is called a mixed pixel. The widespread mixed pixels in hyperspectral data restrict hyperspectral remote sensing. Therefore, how to decompose mixed pixels has become an important research topic in the field of hyperspectral remote sensing. Hyperspectral unmixing, or mixed pixel decomposition, is a method of decomposing mixed pixels into pure substances (called endmembers) and corresponding component ratios (called abundance coefficients). Mineral exploration and other fields have been widely used.

At present, hyperspectral unmixing models can be mainly divided into two categories: linear mixture models and nonlinear mixture models. The nonlinear mixed model considers the multiple scattering effect between endmembers, while the linear mixed model ignores the effect of multiple scattering. It is considered that the mixed pixel spectrum is a linear mixture of the endmember spectrum, so the structure is simple, the physical meaning is clear, and the application is easy. the most extensive.

According to whether the spectral library needs to be provided in the unmixing process, the current hyperspectral data unmixing algorithms can be divided into blind source unmixing algorithms and semi-blind source unmixing algorithms. The blind source unmixing algorithm needs to extract endmembers from the current hyperspectral data, and then perform abundance inversion to obtain abundance coefficients. The semi-blind source unmixing algorithm uses a spectral library composed of a large number of pure endmember spectra as the endmember dictionary, and then selects the appropriate endmembers from the endmember dictionary and calculates the corresponding abundance coefficient. Since the number of endmembers in the spectral library is much larger than the number of endmembers contained in the hyperspectral data, the abundance coefficients will exhibit sparsity, so such algorithms usually add sparsity constraints to the abundance coefficients, for example, Iordache M et al. Considering that the mixed pixels in the local neighborhood are generally composed of the same ground objects, they impose l _2,1 norm constraints on the abundance coefficients, and finally use the Lagrangian method and the alternating direction multiplier method to carry out the model. Solving [Iordache MD, Bioucas-Dias JM, PlazaA. Collaborative sparse regression for hyperspectral unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 52(1): 341-354.]; Rui Wang et al. are using l ₁ On the basis of norm-constrained abundance coefficient, it is proposed to perform double weighting in spectral dimension and spatial dimension, while enhancing the sparsity of endmember matrix and abundance matrix [Wang R, Li HC, Liao W, et al.Double reweightedsparse regression for hyperspectral unmixing[C]//2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). IEEE, 2016:6986-6989.]. However, due to the influence of factors such as the mechanical failure of the sensor, the hyperspectral data will inevitably be polluted by various noises such as Gaussian noise, impulse noise, stripe noise, etc. during the acquisition and transmission process. These methods all default to hyperspectral data. The intensity of Gaussian noise in different bands is the same, and the existence of sparse noise such as impulse noise is ignored, which is often not in line with the actual situation, and the l ₁ norm cannot well express the structural sparsity of abundance coefficients, so that the unmixing algorithm Higher accuracy cannot be achieved.

SUMMARY OF THE INVENTION

The purpose of the present invention is to propose a hyperspectral data unmixing method with robustness to noise and strong sparsity constrained optimization of morphological components for the hyperspectral unmixing problems of mineral exploration, precision agriculture, environmental monitoring and disaster assessment.

The technical solution for realizing the purpose of the present invention is: a hyperspectral data unmixing method for morphological component constraint optimization, the steps are as follows:

Step 1, input hyperspectral data and endmember spectral library;

Step 2, construct the spectrum sample expansion matrix;

Step 3, establish a hyperspectral data unmixing model optimized by morphological component constraints; based on the assumption that the intensity of Gaussian noise in different bands of the hyperspectral data is different and has other sparse noises, by estimating the standard deviation of the band noise and the sparse structure noise decomposition, use l ₁ norm constrains sparse noise, l _2,0 norm constrains the global row sparsity of different pure substance abundance coefficients, and establishes a hyperspectral data unmixing model optimized by morphological component constraints;

Step 4: Calculate the noise standard deviation of each band of the hyperspectral data: first calculate the noise correlation matrix, then calculate the regression vector by band, estimate the noise, and finally calculate the standard deviation of the noise;

Step 5, alternate iterative solution: convert the model into an equivalent augmented Lagrangian form, and solve the optimization problem alternately for each variable according to the alternate direction multiplier method, using the soft threshold limit operator and the hard threshold respectively. The operator solves the subproblems;

Step 6, output the abundance coefficient.

Compared with the prior art, the present invention has the following significant advantages: (1) The influence of Gaussian random noise and sparse structure noise on the linear unmixing accuracy is comprehensively considered, and the _l1 norm is used to constrain the sparse noise, and the _l2,0 norm The global row sparsity of the abundance coefficients of different pure substances is constrained by the number, and the hyperspectral data unmixing model optimized by the morphological component constraints is established; (2) The _l2,0 norm problem is solved by the row hard threshold method to further characterize the abundance coefficients. (3) The simulation experiment results show that the invention has good robustness to mixed noise, and can effectively overcome the difference of waveform morphology and structure between the same batch of spectral data. Fast and high-precision unmixing is realized, and the root mean square error of unmixing is less than 0.0025. The method of the invention has wide application prospects for rock mineral identification and hyperspectral remote sensing ground feature fine identification.

Description of drawings

Fig. 1 is a flow chart of a hyperspectral data unmixing method for morphological component constraint optimization of the present invention.

Figure 2(a) is a true abundance map generated from simulated hyperspectral data.

Figure 2(b) is the abundance map obtained from simulated hyperspectral data using the SUnSAL algorithm.

Figure 2(c) is the abundance map obtained from simulated hyperspectral data using the CLSUnSAL algorithm.

Figure 2(d) is the abundance map obtained from simulated hyperspectral data using CSUnL0 algorithm.

Fig. 2(e) is an abundance map obtained by simulating hyperspectral data using the hyperspectral data unmixing method of morphological component constraint optimization of the present invention.

Figure 3(a) is the original single-strip simulated hyperspectral data.

Figure 3(b) is a single piece of hyperspectral data polluted by sparse noise such as Gaussian noise and impulse noise.

Figure 3(c) is the reconstruction result of a single hyperspectral data by different unmixing algorithms.

Figure 4 is the distribution map of the three minerals and the corresponding abundance maps obtained by different unmixing algorithms on the Cuprite dataset.

Detailed ways

The invention discloses a hyperspectral data unmixing method for morphological component constraint optimization. The unmixing model is established by estimating band noise standard deviation and sparse structure noise decomposition, adding sparsity constraints to sparse components and abundance coefficients, and iterating alternately. To realize the linear unmixing of the mixed spectrum, as shown in Figure 1, the specific steps of the method are:

y_i∈R^B和端元光谱库E∈R^B×M，B是波段数，P是像元数，M是端元数。Step 1, input hyperspectral data and endmember spectral library. Input a batch of hyperspectral data to be unmixed

y _i ∈ R ^B and endmember spectral library E∈R ^B×M , where B is the number of bands, P is the number of pixels, and M is the number of endmembers.

Step 2, construct the spectral sample expansion matrix. The input hyperspectral data is arranged according to the pixel-by-pixel spectral vector to form a spectrum-pixel matrix Y=[y ₁ , y ₂ ,...,y _N ]∈R ^B×P .

Step 3, establish a hyperspectral data unmixing model optimized by morphological component constraints. Based on the assumption that the intensity of Gaussian noise in different bands of hyperspectral data is different and has other sparse noise, a hyperspectral data unmixing model is established:

Y=EA+S+N

Among them, A∈R ^M×P is the abundance coefficient, and the elements in each column represent the proportion of the corresponding endmember spectrum in a single mixed pixel; S∈R ^B×P is the sparse component, indicating that only the hyperspectral is polluted Mixed noises such as impulse noise and stripe noise in a few bands of data; N∈R ^B×P is Gaussian noise. Introducing the sparse constraint, the optimal solution function of the above model is obtained:

s.t.A≥0

stA≥0

对角上的元素为各波段高斯噪声的标准差的倒数；‖·‖_F表示矩阵的F范数；

表示矩阵的l_2,0范数；

表示矩阵的l₁范数。通过增广拉格朗日乘子法将上式转化为：Among them, λ and α are regularization parameters, λ>0, α>0; W is a diagonal matrix,

The elements on the diagonal are the reciprocal of the standard deviation of the Gaussian noise in each band; ‖·‖ _F represents the F-norm of the matrix;

represents the l _2,0 norm of the matrix;

represents the l ₁ norm of the matrix. The above equation is transformed into:

X_i,j表示矩阵X第i行第j列的元素，当X_i,j为非负值时，

等于0，否则等于正无穷；Among them, V ₁ , V ₂ , V ₃ are auxiliary variables, D ₁ , D ₂ , D ₃ are Lagrange multipliers, μ>0 is the penalty factor,

X _i,j represents the element in the i-th row and the j-th column of the matrix X. When X _i,j is a non-negative value,

equal to 0, otherwise equal to positive infinity;

Step 4: Calculate the noise standard deviation of each band of the hyperspectral data to obtain a diagonal matrix W. Proceed as follows:

Step 4.1, calculate the regression vector

表示矩阵

的第i行，Z＝Y^T，

表示相关矩阵，

表示从矩阵

中删除第

行和第

列后获得的矩阵，

表示矩阵

的第i行，

i＝1,…,L；in,

representation matrix

The ith row of , Z=Y ^T ,

represents the correlation matrix,

represents the matrix from

delete the

row and

The matrix obtained after columns,

representation matrix

The i-th row of ,

i=1,...,L;

Step 4.2, Calculate the noise vector

表示估计噪声

的第i行，z_i表示矩阵Z的第i行，

表示从矩阵Z中删除第

列后获得的矩阵；in,

represents the estimated noise

The ith row of , _zi represents the ith row of matrix Z,

means to remove the first from the matrix Z

The matrix obtained after the columns;

Step 4.3, calculate the reciprocal W _ii of the standard deviation of the Gaussian noise of each band, and get the matrix W:

的均值，P表示光谱采样点个数。Among them, β represents

The mean of , and P represents the number of spectral sampling points.

Step 5, alternate iterative solution. Proceed as follows:

Step 5.1, fix other variables, update the abundance coefficient A,

Step 5.2, fix other variables, update the sparse component S,

Among them, S _τ (·) is the soft threshold limit operator, and S _τ (x)=sgn(x)max(|x|-τ,0), where τ is a non-negative parameter, and sgn(·) represents the sign function ;

Step 5.3, fix other variables, update auxiliary variables V ₁ , V ₂ , V ₃ ,

表示行硬阈值算子，且

其中τ为非负参数，B(i,:)表示矩阵B的第i行；in,

represents the row hard threshold operator, and

where τ is a non-negative parameter, and B(i,:) represents the i-th row of matrix B;

Step 5.4, fix other variables, update the Lagrange multipliers D ₁ , D ₂ , D ₃ ,

Step 6, output the abundance coefficient A.

The effect of the present invention can be further illustrated by the following simulation experiments:

Simulation conditions

Simulation experiments use simulated hyperspectral data and real hyperspectral data (Cuprite dataset). Spectral Library The Spectral Library comes from the ground object spectral database released by the United States Geological Survey (USGS), including 224 bands and 240 endmembers of pure substances. For the simulated hyperspectral data, each band contains 3136 pixels. In order to better simulate the actual situation, all the bands are polluted by Gaussian noise of different intensities (signal-to-noise ratio between 20dB and 40dB), 11 bands (60～40dB). 70) Contaminated by 30% impulse noise, 11 bands (120～130) are polluted by stripe noise, that is, the input hyperspectral data is Y=[y ₁ ,y ₂ ,...,y ₃₁₃₆ ]∈R ^{224 ×3136} , and the spectral library is E∈R ^224×240 . The real hyperspectral data comes from the AVIRIS (Airborne Visible Infrared Imaging Spectrometer) Cuprite dataset, after removing 36 water vapor absorption and low signal-to-noise ratio bands (band numbers are 1-2, 105-115, 150-170, 223-224 ), select the remaining 188 bands as the research object, and correspondingly, remove the bands corresponding to the input spectral library, that is, the input hyperspectral data is Y=[y ₁ , y ₂ ,...,y ₄₇₇₅₀ ]∈ R ^188×47750 , the spectral library is E∈R ^{188 ×240} .

The evaluation index used in this simulation experiment is Root Mean Square Error (RMSE, Root Mean Square Error). The smaller the RMSE, the higher the unmixing accuracy.

Simulation content

The present invention uses simulated hyperspectral data and real hyperspectral data (Cuprite data set) to test the unmixing performance of the algorithm. In order to test the performance of the algorithm of the present invention, a hyperspectral data unmixing method (SUnMC) with morphological component constraint optimization proposed is compared with the currently popular unmixing algorithm in the world. The comparison methods include: Sparse Unmixing Algorithm Based on Variable Splitting and Augmented Lagrangian (SUnSAL), Cooperative Sparse Unmixing Algorithm Based on Variable Splitting and Augmented Lagrangian (CLSUnSAL), Synergistic Based on l ₀ Norm Sparse unmixing algorithm (CSUnL0).

Analysis of simulation results

Table 1 shows the comparison results of the root mean square error of simulated hyperspectral data under different unmixing algorithms:

Table 1 RMSE of different algorithms on simulated hyperspectral data

It can be seen from Table 1 that under the condition of simulating hyperspectral data, since SUnMC takes into account the different intensity of Gaussian noise in different bands and the existence of sparse noise, its root mean square error is significantly lower than that of other algorithms. Mixing works best.

Figures 2(a) to 2(e) are abundance maps obtained on simulated hyperspectral data using different unmixing algorithms. It can be seen that under the influence of sparse noise such as Gaussian noise and impulse noise, the abundance maps obtained by SUnSAL, CLSUnSAL and CSUnL0 contain a lot of noise, and the unmixing accuracy is not high. With better robustness, results closer to the true abundance map can be obtained.

In order to verify the effectiveness of SUnMC for single spectral unmixing, we take all the bands on a single pixel in the simulated hyperspectral data as input, that is, the input single hyperspectral data is y ∈ R ^224×1 and the spectral library is E ∈ R ^{224 ×240} . Figures 3(a) to 3(c) show the original spectrum, the spectrum polluted by sparse noise such as Gaussian noise and impulse noise, and the reconstruction results of different unmixing algorithms, respectively. It can be seen from Fig. 3(c) that, compared with other algorithms, the reconstructed spectrum of SUnMC is closer to the original spectrum, indicating that SUnMC has good anti-noise performance and can obtain better results.

Figure 4 shows the abundance maps of alumite, feldspar, and chalcedony obtained by the SUnSAL, CLSUnSAL, CSUnL0 and SUnMC algorithms on the Cuprite dataset. In order to facilitate the comparative analysis of abundance maps, we used Tetracorder, a hyperspectral mineral mapping software developed by the United States Geological Survey (USGS), to identify the mineral information contained in the spectrum and generate a mineral distribution map, as shown in the left column of Figure 4 . Each pixel in the figure indicates whether the point belongs to a certain mineral, and does not show the content of a certain mineral at that point, so we only use it as a qualitative comparative analysis here, not as a real abundance map. The rightmost column in Figure 4 is the unmixing result of the SUnMC algorithm, with hotter hues indicating higher levels of the mineral. For alumite minerals, the abundance maps obtained by the SUnSAL, CLSUnSAL, CSUnL0 and SUnMC algorithms have little difference, and are very similar to the results obtained by Tetracorder. For the two minerals, feldspar and chalcedony, from the visual effect, the abundance map obtained by the SUnMC algorithm has more details, which proves that our algorithm is more effective than other algorithms.

Claims (6)

1. a hyperspectral data unmixing method of morphological component constraint optimization, is characterized in that, comprises the following steps: Step 1, input hyperspectral data and endmember spectral library; Step 2, construct the spectrum sample expansion matrix; Step 3, establish a hyperspectral data unmixing model optimized by morphological component constraints: based on the assumption that the intensity of Gaussian noise in different bands of hyperspectral data is different and has other sparse noises, by estimating the standard deviation of band noise and sparse structure noise decomposition, use l ₁ norm constrains sparse noise, l _2,0 norm constrains the global row sparsity of different pure substance abundance coefficients, and establishes a hyperspectral data unmixing model optimized by morphological component constraints; Step 4: Calculate the noise standard deviation of each band of the hyperspectral data: first calculate the noise correlation matrix, then calculate the regression vector by band, estimate the noise, and finally calculate the standard deviation of the noise; Step 5, alternate iterative solution: transform the model into an equivalent augmented Lagrangian form, and solve the optimization problem alternately for each variable according to the alternate direction multiplier method; Step 6, output the abundance coefficient. 2. The hyperspectral data unmixing method of morphological component constraint optimization according to claim 1, wherein the input hyperspectral data and endmember spectral library described in step 1 are as follows: Input a batch of hyperspectral data to be unmixed

and endmember spectral library E∈R ^B×M , where B is the number of bands, P is the number of spectral sampling points, and M is the number of endmembers. 3. The hyperspectral data unmixing method of morphological component constraint optimization according to claim 1, characterized in that, the construction of the spectral sample expansion matrix described in step 2 is as follows: The input hyperspectral data is arranged according to the pixel-by-pixel spectral vector to form a spectrum-pixel matrix Y=[y ₁ , y ₂ ,...,y _N ]∈R ^B×P . 4. The hyperspectral data unmixing method of morphological component constraint optimization according to claim 1, is characterized in that, establishing the hyperspectral data unmixing model of morphological component constraint optimization described in step 3 is as follows: Based on the assumption that the intensity of Gaussian noise in different bands of hyperspectral data is different and has other sparse noise, a hyperspectral data unmixing model is established: Y=EA+S+N (1) In formula (1), A∈R ^M×P is the abundance coefficient, and the elements in each column represent the proportion of the corresponding endmember spectrum in a single mixed pixel; S∈R ^B×P is the sparse component, N∈ R ^B×P is Gaussian noise; Introducing the sparse constraint, the optimal solution function of the model described in equation (1) is obtained:

In formula (2), λ and α are regularization parameters, λ>0, α>0; W is a diagonal matrix,

The elements on the diagonal are the reciprocal of the standard deviation of the Gaussian noise in each band; ‖·‖ _F represents the F-norm of the matrix;

represents the L _2,0 norm of the matrix;

Represents the l ₁ norm of the matrix; Equation (2) is transformed into:

In formula (3), V ₁ , V ₂ , V ₃ are auxiliary variables, D ₁ , D ₂ , D ₃ are Lagrange multipliers, μ>0 is a penalty factor,

X _i,j represents the element in the i-th row and the j-th column of the matrix X. When X _i,j is a non-negative value,

equal to 0, otherwise equal to positive infinity. 5. The hyperspectral data unmixing method of morphological component constraint optimization according to claim 1, characterized in that, calculating the noise standard deviation of each band of hyperspectral data described in step 4 to obtain a diagonal matrix, comprising the following steps: Calculate the regression vector

In formula (4),

representation matrix

The ith row of , Z=Y ^T ,

represents the correlation matrix,

represents the matrix from

delete the

row and

The matrix obtained after columns,

representation matrix

The i-th row of ,

Calculate the noise vector

In formula (5),

represents the estimated noise

The ith row of , _zi represents the ith row of matrix Z,

means to remove the first from the matrix Z

The matrix obtained after the columns; Compute the reciprocal Wi _ii of the standard deviation of the Gaussian noise for each band to get the matrix W:

In formula (6), β represents

The mean of , and P represents the number of spectral sampling points. 6. The hyperspectral data unmixing method of morphological component constraint optimization according to claim 1, wherein in step 5, the linear unmixing of the mixed spectrum is realized alternately and iteratively, comprising the following steps: Fix other variables and update the abundance coefficient A:

Fix other variables, update the sparse component S:

In formula (8), S _τ (·) is the soft threshold limit operator, and S _τ (x)=sgn(x)max(|x|-τ,0), where τ is a non-negative parameter, sgn(· ) represents a symbolic function; Fix other variables and update auxiliary variables V ₁ , V ₂ , V ₃ :

In formula (10),

represents the row hard threshold operator, and

where τ is a non-negative parameter, and B(i,:) represents the i-th row of matrix B; Fix other variables and update the Lagrange multipliers D ₁ , D ₂ , D ₃ :

Images