CN104036509B - Method for unmixing hyperspectral mixed pixel based on compressed sensing - Google Patents

Abstract

The invention relates to a method for unmixing a hyperspectral mixed pixel based on compressed sensing, and relates to the compressed sensing field and the hyperspectral remote sensing field. The method is used for solving the problem of low mixed pixel unmixing speed during hyperspectral image data acquisition by use of the traditional Nyquist sampling theorem. The method comprises the steps of firstly, inputting a measurement matrix Phi and a compressed measurement matrix Y and establishing a spectral mixing model Y=PhiXT=Phi(AS)T by use of the compressed sensing theory, secondly, performing iterative processing on the estimated value S^ of an end member abundance matrix S and the estimated value A^ of an end member spectrum matrix A, and if the difference of absolute values of every corresponding element in the estimated values A^ of the end member spectrum matrix A obtained by two adjacent times of iterative processing is smaller than 0.1, stopping iteration and outputting the end member abundance matrix S^ and completing the unmixing of the hyperspectral mixed pixel, otherwise, continuing the iterative processing. The method is mainly applied to unmixing the hyperspectral mixed pixels.

Description

Hyperspectral mixed pixel decomposition method based on compressed sensing

Technical Field

The invention relates to the field of compressed sensing and the field of hyperspectral remote sensing.

Background

The hyperspectral imaging technology is a novel earth observation technology developed in the field of remote sensing, and typical hardware equipment is an imaging spectrometer. The imaging spectrometer decomposes the electromagnetic wave signal into many tiny, adjacent bands by spectroscopic techniques, and the energy on the corresponding band is received by different sensors. Therefore, compared with the traditional remote sensing imaging technology, the hyperspectral imaging technology has the characteristics of integrated maps, multiple spectral bands, high spectral resolution and the like, and has great advantages in the aspects of surface material identification and classification, feature extraction and the like. The spatial resolution of the hyperspectral image is low, so that mixed pixels are inevitably generated, and the existence of the mixed pixels becomes an obstacle to the further development of the hyperspectral remote sensing application. If the mixed pixels of the hyperspectral image are decomposed to obtain the abundance corresponding to each end member, the method has very important significance in performing blind decomposition on the hyperspectral remote sensing data.

The blind source separation technique is a process of recovering a source signal only from an observed mixed signal according to the statistical characteristics of the source signal under the condition that the source signal and a mixing mode are unknown. "Blind" has two explanations: one is that besides the observation signal, other information is unknown, which is called as 'totally blind'; the other is that some a priori information of the signal or hybrid system is known, called "semi-blind". Aiming at the problem of mixed pixel decomposition, the blind source separation technology simultaneously obtains an end member spectrum and end member abundance, namely the blind decomposition of the mixed pixel is realized. Before the problem of mixed pixel decomposition is solved, a spectrum mixing model needs to be established. The spectrum mixing mode mainly comprises linear mixing and nonlinear mixing, the linear spectrum mixing model is a mixing pixel decomposition model which is most widely applied at present, and the expression is as follows:

$X=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{s}_i=\mathrm{AS}---(1-1)$

wherein X is ═ X₁,...x_l,...x_L]^T∈R^L×1Is a mixed pixel spectrum of L wave bands of a pixel, is a known observed quantity A ∈ R^L×PIs an end-member spectral matrix, wherein each column corresponds to a spectral vector of an end-member, P is the number of end-members, S ∈ R^{P ×1}And the abundance vector of the mixed pixel is obtained. The abundance vector needs to satisfy two constraints: a non-negative constraint, and a 1 constraint. Namely, it is

s_i≥0,i∈{1,2,...,P} (2-1)

$\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_i=1,i\∈\{\mathrm{1,2},...,P\}---(3-1)$

If the hyperspectral image has N picture elements, the formula (1-1) can be further extended to formula (4-1).

X＝AS (4-1)

X ∈ R in formula (4-1)^L×NThe column vector of which represents the mixed spectral vector of each pixel in the L bands, S ∈ R^P×NAn end-member abundance matrix is constructed. The end-member abundance matrix meets the following constraint conditions:

s_i,j≥0,i∈{1,2,...,P},j∈{1,2,...,N} (5-1)

$\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_{i,j}=1,j\∈\{\mathrm{1,2},...,N\}---(6-1)$

the method is characterized in that a blind source separation technology-based hyperspectral image mixed pixel decomposition problem is that under the condition that an end member matrix and an abundance matrix are unknown, an end member spectrum and end member abundance are recovered by using an observed value of a mixed pixel spectrum.

The characteristic that the number of wave bands of the hyperspectral image is large enables the data volume of the hyperspectral image to be huge, and if the traditional Nyquist sampling theorem is adopted to collect data, huge pressure is generated on the sampling speed, the transmission speed and the storage capacity of hardware. The problem is solved by the proposal of the compressive sensing theory, which indicates that when the signal satisfies sparsity or is compressible, the signal can be globally observed at a speed far lower than the Nyquist sampling rate, and then the source signal is recovered by adopting a compressive sensing reconstruction algorithm. The hyperspectral image is sparse in the wavelet domain, so that the compressive sensing theory can be applied to sampling of the hyperspectral image. The specific sampling procedure can be expressed as follows:

Y＝ΦX^T＝Φ(AS)^T

in the above formula,. phi. ∈ R^M×NFor the observation matrix, X ∈ R^L×NBeing a mixed-pixel spectral matrix, Y ∈ R^M×LThe observation matrix is compressed.

Disclosure of Invention

The invention provides a hyperspectral mixed pixel decomposition method based on compressed sensing, and aims to solve the problem that when the traditional Nyquist sampling theorem is adopted to acquire hyperspectral image data, the mixed pixel decomposition speed is low.

A method for decomposing a hyperspectral mixed pixel based on compressed sensing comprises the following steps,

inputting an observation matrix phi and a compressed observation matrix Y, and establishing a spectrum mixed model by using a compressed sensing theory:

Y＝ΦX^T＝Φ(AS)^T(1)

Φ∈R^M×Nan observation matrix of M × N, R is a real number,

X∈R^L×Na mixed pixel spectral matrix of L × N,

Y∈R^M×La compressed observation matrix of M × L,

S∈R^P×Nis an end-member abundance matrix of P × N,

A∈R^L×Pan end-member spectral matrix of L × P,

initializing, and randomly selecting an end member spectrum matrix A as an estimated value of the end member spectrum matrix AAnd isIs a matrix of L × P and,

let the estimate of the end-member abundance matrix S:

$\hat{S}={\left[\begin{array}{ccc}0& ...& 0\\ & ...& \\ 0& & 0\end{array}\right]}_{N\×P}---\left(2\right)$

wherein,is a matrix of N × P and,

step three, setting an initial value of an iteration variable t as 1;

step four, estimating the abundance matrix S of the end-to-end elementEstimate of sum-end spectral matrix ACarrying out iterative processing;

step five, multiplying the two sides of the formula (1) simultaneouslyThe pseudo-inverse of (1) is transformed into

Y₁＝ΦS^T(3)，

Step six, if the wavelet basis is theta and sigma is a sparse coefficient of the abundance of the end members, the formula (3) is deformed into

Y₁＝Φ(ΘΣ)^T(4)，

Solving the following formula by using a BP algorithm:

argmin||Σ||₁s.t. Y₁＝Φ(ΘΣ)^T(5)，

obtaining an estimate of an end-member abundance matrix S

$\hat{S}=---\mathrm{\Θ\Σ}\left(6\right),$

Step seven, the probability density function is P (S)_k)∝α exp(α|S_k|) normalized end-member spectral matrix AUpdating the estimated value of the end-member spectral matrix A，

$\hat{A}=A_1+\mathrm{\λ}\{-A_1(B{\hat{S}}^T+I)\}---\left(7\right),$

Wherein A is₁For the estimated value of the end member spectrum matrix obtained in the last iteration, λ is the iteration step length, I is the identity matrix, B is the vector set, and B ═ B₁，B₂，B₃，......B_k}，Is the prior probability distribution of the kth source signal, k being an integer, α being a real number greater than zero,

step eight, if the estimated value of the end member spectrum matrix A obtained by two adjacent timesStopping iteration and executing the step nine if the difference of the absolute values of each corresponding element in the step nine is less than 0.1, otherwise, making t equal to t +1 and returning to the step four;

ninth, output end member abundance matrixAnd (4) decomposing the hyperspectral mixed pixel.

The method has the advantages that the average signal-to-noise ratio SNR is adopted for the evaluation index of the end element abundance matrix S_avgThe realization method is realized in the way that,

${\mathrm{SNR}}_{\mathrm{avg}}=\frac{1}{P}\underset{k=1}{\overset{P}{\mathrm{\Σ}}}20\×{\log }_{10}\left(\frac{\left|\right|S_k{\left|\right|}_2}{{\left|\right|S_k-{\hat{S}}_k\left|\right|}_2}\right)---\left(8\right)$

the number of end members P of the hyperspectral image used in the experiment is 3, the number of wave bands L is 40, and the number of pixels N of the hyperspectral image is 32 × 32. When the compression sensing sampling is carried out, the compression rate is 0.8, the average signal-to-noise ratio of the end member abundance matrix decomposed by the mixed pixels is 17.78dB, the separation can be better realized, and the decomposition speed of the mixed pixels is improved by more than 40%.

Detailed Description

The first embodiment is as follows: the method for decomposing the hyperspectral mixed pixel based on compressed sensing in the embodiment comprises the following steps,

inputting an observation matrix phi and a compressed observation matrix Y, and establishing a spectrum mixed model by using a compressed sensing theory:

Y＝ΦX^T＝Φ(AS)^T(1)

Φ∈R^M×Nan observation matrix of M × N, R is a real number,

X∈R^L×Na mixed pixel spectral matrix of L × N,

Y∈R^M×La compressed observation matrix of M × L,

S∈R^P×Nis an end-member abundance matrix of P × N,

A∈R^L×Pan end-member spectral matrix of L × P,

initializing, and randomly selecting an end member spectrum matrix A as an estimated value of the end member spectrum matrix AAnd isIs a matrix of L × P and,

let the estimate of the end-member abundance matrix S:

$\hat{S}={\left[\begin{array}{ccc}0& ...& 0\\ & ...& \\ 0& & 0\end{array}\right]}_{N\×P}---\left(2\right)$

wherein,is a matrix of N × P and,

step three, setting an initial value of an iteration variable t as 1;

step four, estimating the abundance matrix S of the end-to-end elementEstimate of sum-end spectral matrix ACarrying out iterative processing;

step five, multiplying the two sides of the formula (1) simultaneouslyThe pseudo-inverse of (1) is transformed into

Y₁＝ΦS^T(3)，

Step six, if the wavelet basis is theta and sigma is a sparse coefficient of the abundance of the end members, the formula (3) is deformed into

Y₁＝Φ(ΘΣ)^T(4)，

Solving the following formula by using a BP algorithm:

argmin||Σ||₁s.t. Y₁＝Φ(ΘΣ)^T(5)，

obtaining an estimate of an end-member abundance matrix S

$\hat{S}=---\mathrm{\Θ\Σ}\left(6\right),$

Step seven, the probability density function is P (S)_k)∝α exp(α|S_k|) normalized end-member spectral matrix AUpdating the estimated value of the end-member spectral matrix A

$\hat{A}=A_1+\mathrm{\λ}\{-A_1(B{\hat{S}}^T+I)\}---\left(7\right),$

Wherein A is₁For the estimated value of the end member spectrum matrix obtained in the last iteration, λ is the iteration step length, I is the identity matrix, B is the vector set, and B ═ B₁，B₂，B₃，......B_k}，Is the prior probability distribution of the kth source signal, k being an integer, α being a real number greater than zero,

step eight, if the estimated value of the end member spectrum matrix A obtained by two adjacent timesStopping iteration and executing the step nine if the difference of the absolute values of each corresponding element in the step nine is less than 0.1, otherwise, making t equal to t +1 and returning to the step four;

ninth, output end member abundance matrixAnd (4) decomposing the hyperspectral mixed pixel.

Claims (1)

1. A method for decomposing a hyperspectral mixed pixel based on compressed sensing is characterized in that,

inputting an observation matrix phi and a compressed observation matrix Y, and establishing a spectrum mixed model by using a compressed sensing theory:

Y＝ΦX^T＝Φ(AS)^T(1)

Φ∈R^M×Nan observation matrix of M × N, R is a real number,

X∈R^L×Na mixed pixel spectral matrix of L × N,

Y∈R^M×La compressed observation matrix of M × L,

S∈R^P×Nis an end-member abundance matrix of P × N,

A∈R^L×Pan end-member spectral matrix of L × P,

initializing, and randomly selecting an end member spectrum matrix A as an estimated value of the end member spectrum matrix AAnd isIs a matrix of L × P and,

let the estimate of the end-member abundance matrix S:

$\hat{S}={\left[\begin{array}{ccc}0& \mathrm{...}& 0\\ & \mathrm{...}& \\ 0& & 0\end{array}\right]}_{P\×N}---\left(2\right)$

wherein,is a matrix of P × N and,

step three, setting an initial value of an iteration variable t as 1;

step four, estimating the abundance matrix S of the end-to-end elementEstimate of sum-end spectral matrix ACarrying out iterative processing;

step five, multiplying the two sides of the formula (1) simultaneouslyThe pseudo-inverse of (1) is transformed into

Y₁＝ΦS^T(3)，

Step six, if the wavelet basis is theta and sigma is a sparse coefficient of the abundance of the end members, the formula (3) is deformed into

Y₁＝Φ(ΘΣ)^T(4)，

Solving the following formula by using a BP algorithm:

arg min||Σ||₁s.t. Y₁＝Φ(ΘΣ)^T(5)，

obtaining an estimate of an end-member abundance matrix S

$\hat{S}=\mathrm{\Θ}\Σ---\left(6\right),$

Step seven, the probability density function is P (S)_k)∝αexp(α|S_k|) normalized end-member spectral matrix AUpdating the estimated value of the end-member spectral matrix A

$\hat{A}=A_1+\mathrm{\λ}\{-A_1(B{\hat{S}}^T+I)\}---\left(7\right),$

Wherein A is₁For the estimated value of the end member spectrum matrix obtained in the last iteration, λ is the iteration step length, I is the identity matrix, B is the vector set, and B ═ B₁，B₂，B₃，......B_k}， Is the prior probability distribution of the kth source signal, k being an integer, α being a real number greater than zero, S_kRepresents the k-th end-member abundance matrix S,represents an estimate of the k-th end-member abundance matrix S,

step eight, if the estimated value of the end member spectrum matrix A obtained by two adjacent timesWhen the difference between the absolute values of each corresponding element in (b) is less than 0.1, thenStopping iteration, executing the step nine, otherwise, making t equal to t +1, and returning to the step four;

ninth, output end member abundance matrixAnd (4) decomposing the hyperspectral mixed pixel.