CN103745487B - The mixed compression sensing method of Bayes's EO-1 hyperion solution based on structural sparse priori - Google Patents

Abstract

The invention discloses the mixed compression sensing method of a kind of Bayes's EO-1 hyperion solution based on structural sparse priori, for solving the technical problem of the mixed Compression of hyperspectral images cognitive method low precision of existing combined spectral solution.Technical scheme is that compression process mixes the low-rank character excavating high-spectral data inherence by linear solution, uses wavelet transformation that Abundances converts to structural sparse signal, uses compressed sensing to obtain compression data afterwards.Process of reconstruction, suitable end member matrix is selected from library of spectra, introduce the structural sparse priori of Abundances wavelet coefficient, then use the Bayesian inference method accurate reconstruction Abundances matrix based on gibbs sampler, finally use linear mixed model to rebuild initial data.About 10% is promoted relative to background technology compressed sensing class arithmetic accuracy.

Description

Bayes hyperspectral unmixing compressive sensing method based on structured sparse prior

Technical Field

The invention relates to a hyperspectral unmixing compression sensing method, in particular to a Bayes hyperspectral unmixing compression sensing method based on structured sparse prior.

Background

Due to the huge data volume of the hyperspectral images and the compression algorithm with high compression ratio, the hyperspectral image compression method has profound significance for acquisition, transmission and processing of the hyperspectral images. The existing hyperspectral image compression algorithms are mainly divided into two types, one type is a traditional compression method based on information coding, the method is expanded from a compression method of a common image, and focuses on removing redundancy inside each wave band and between wave bands in the hyperspectral image for compression, wherein the method comprises clustered pulse-based differential-coding (CDPCM), three-dimensional wavelet transform (3D-DWT), three-dimensional discrete cosine transform (3D-DCT) and the like, but the compression algorithm can only be compressed after image acquisition, still needs a large amount of hardware for acquiring and temporarily storing data, and the compression ratio is usually low; another method is a compression method based on compressive sensing (compressive sensing), which can accurately reconstruct the original sparse signal by collecting a small number of sample points, and provides a theoretical basis for reducing the hardware consumption required for collecting data, and the method focuses more on the reconstruction process.

The literature, "Acompressivegiven and unknown matching spectral image compression sensing, IEEETransactionson Imageprocessing,2012,21(3):1200 1210" discloses a hyperspectral image compression sensing method combining spectral unmixing. Decomposing hyperspectral data into an abundance value matrix and an end member matrix through linear unmixing; the gradient of the abundance value matrix in the spatial dimension has sparsity, and the compressed reconstruction can be carried out by adopting a compressed sensing mode. Firstly, obtaining compressed data by using a compressed sensing matrix; selecting proper end members by combining the spectrum library; then, reconstructing an abundance value matrix containing sparse gradients by using a compressed sensing method; and finally, reconstructing hyperspectral data to be acquired by using the linear mixed model. However, the structural sparsity of a deeper layer of an abundance value matrix cannot be excavated through the gradient sparsity, and the method for reconstructing the abundance value by using the augmented Lagrange method cannot adaptively select the proportion factor of the gradient sparsity regular term. These two factors directly affect the reconstruction accuracy of the compression algorithm.

Disclosure of Invention

In order to overcome the defect that the existing hyperspectral image compressed sensing method combining spectrum unmixing is poor in precision, the invention provides a Bayes hyperspectral unmixing compressed sensing method based on structured sparse prior. According to the method, the inherent low-rank property of the hyperspectral data is explored through linear unmixing in the compression process, the abundance value is converted into a structured sparse signal through wavelet transformation, and then compressed data are obtained through compressed sensing. In the reconstruction process, a proper end member matrix is selected from a spectrum library, structured sparse prior of an abundance value wavelet coefficient is introduced, then a Bayesian inference method based on Gibbs sampling is used for accurately reconstructing the abundance value matrix, and finally a linear mixed model is used for reconstructing original data. The data of the invention is on the urban data shot by a HYDICE satellite, when the compression ratio is 1:16, the normalized mean square error is less than 0.4, and on the Indiana data shot by AVIRIS, when the compression ratio is 1:16, the normalized mean square error is also less than 0.4, and the precision is improved by about 10 percent compared with the precision of a compressed sensing algorithm in the background technology.

The technical scheme adopted by the invention for solving the technical problems is as follows: a Bayes hyperspectral unmixing compressive sensing method based on structured sparse prior is characterized by comprising the following steps:

step one, aiming at high spectrum datan_pRepresenting the total number of pixels, n, contained in the hyperspectral data space_bRepresenting the number of bands, x, contained in the hyperspectral data_i,i＝1,...,n_bThe image representing each band is spread by rows to form a column vector. Each pixel is at n_bThe different reflection values over the individual wavelength bands constitute the spectrum of the pixel. This spectrum is described using a linear mixing model. The linear mixture model considers any one mixtureThe resultant spectrum is a linear combination of all end members, and the proportion of the end members in the mixed spectrum is called abundance value. X represents a linear combination of the abundance value matrix and the end-member matrix:

X＝HW（1）

wherein,is a matrix of abundance values and is,represents the proportion of the ith end member in the mixed spectrum of each pixel, n_eThe number of end-members is indicated,is an end-member matrix in whichRepresenting each end member.

Using a discrete wavelet transform for the abundance value:

Θ＝Ψ^TH（2）

Y＝Ψ^TX＝Ψ^THW＝ΘW

wherein,is a wavelet basis and is a composite of a plurality of wavelets,wavelet coefficients representing abundance values, and Y represents wavelet coefficients of the entire data.Represents h_iHas structural sparsity.

Secondly, reconstructing an abundance value matrix by utilizing the wavelet coefficient, and finally realizing the reconstruction of the whole data:

1. generation of m × n using random variables satisfying a Gaussian distribution_pA matrix of sizes, and each column is normalized to obtain a random perceptual matrixSampling the wavelet transformed original data by phi to obtain compressed data

F＝ΦY＝ΦΨ^TX（3）

Wherein m represents a pair length of n_pAfter compression of the signal, typically m < n_p。

2. For a scale-limited scenario, the influence of environmental factors on the spectrum is ignored, assuming a limited number of end-members. Using the ASTER spectral library, portions of the spectrum are selectively extracted from the spectral library as end-members for a particular scene. Finally determining that n is contained_eA matrix W of individual end-members.

3. And estimating an abundance value matrix from the compressed data by adopting a Bayesian inference mode.

(1) Obtaining data obtained by compressing wavelet coefficients of an abundance value matrix independently by using a least square method through formula (2) and formula (3)

G＝FW^T(WW^T)^-1Phi theta + N (4) where,representing the noise term introduced using the least squares method.

(2) Obtaining a vector form g of the formula (4)_i＝Φθ_i+n_iAssuming noise n_iObedience mean 0 and covarianceGaussian distribution ofAnd the noise vectors are independently and identically distributed, then theta is obtained_i,i＝1,2,...,n_eLikelihood function of (d):

$p\left(g_i\right|{\mathrm{\&theta;}}_i,{\mathrm{\&beta;}}_i^{-1})={\left(2\mathrm{\&pi;}{\mathrm{\&beta;}}_i^{-1}\right)}^{-\frac{m}{2}}\exp (-\frac{1}{2{\mathrm{\&beta;}}_i^{-1}}{\left|\right|g_i-{\mathrm{\&Phi;\&theta;}}_i\left|\right|}_2^2)---\left(5\right)$

β_ifor the accuracy of the noise vector, β is specified_iIs a gamma distribution as follows:

p(β_i|κ₀,υ₀)＝Gamma(κ₀,υ₀)（6）

κ₀＝10^-6,υ₀＝10^-6the shape parameter and the scale parameter of the gamma distribution.

(3) Establishing theta_iThe structured sparsity of each element in (a). The structured sparse prior is mainly reflected in a tree structure of wavelet coefficients, when a parent node is 0, a child node thereof is highly likely to be 0, and vice versa:

$p\left({\mathrm{\&theta;}}_{s,\mathrm{ji}}\right|{\mathrm{\&beta;}}_{s,i}^{-1})=(1-{\mathrm{\&pi;}}_{s,\mathrm{ji}}){\mathrm{\&delta;}}_0+{\mathrm{\&pi;}}_{s,\mathrm{ji}}N(0,{\mathrm{\&beta;}}_{s,i}^{-1})---\left(7\right)$

wherein, s represents the total number of scales of the wavelet transform, and θ represents the total number of scales of the wavelet transform, wherein s is 1_s,jiDenotes theta_iThe jth element at scale s.₀Is a distribution that is completely centered at 0. Pi_s,jiIs a proportionality coefficient when_s,jiTowards 0, theta_s,jiIs 0, otherwise, pi_s,jiWhen going to 1, theta_s,jiFrom a gaussian distribution. Pi_s,jiDetermined by the tree structure of the wavelet transform:

${\mathrm{\&pi;}}_{s,\mathrm{ji}}=\left\{\begin{array}{cc}{\mathrm{\&pi;}}_r& \mathrm{ifs}=1\\ {\mathrm{\&pi;}}_s^0& \mathrm{if}2\&le;s\&le;L,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}=0\\ {\mathrm{\&pi;}}_s^1& \mathrm{if}2\&le;s\&le;L,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}\&NotEqual;0\end{array}\right.---\left(8\right)$

θ_pa(s,ji)is theta_s,jiParent node of when theta_pa(s,ji)When equal to 0, pi_s,jiSet to a high probability of 0When theta is_pa(s,ji)Not equal to 0, pi_s,jiSet to a high probability of 1β_s,iIs theta_s,jiThe accuracy of the gaussian distribution in the prior. Pi_r，And β_s,iThe prior distribution of (a) is as follows:

$\begin{array}{cc}{\mathrm{\&pi;}}_r~\mathrm{Beta}({\mathrm{\&epsiv;}}_0^r,{\mathrm{\&gamma;}}_0^r),& {\mathrm{\&pi;}}_s^0~\mathrm{Beta}({\mathrm{\&epsiv;}}_0^s,{\mathrm{\&gamma;}}_0^s),\\ {\mathrm{\&pi;}}_s^1~\mathrm{Beta}({\mathrm{\&epsiv;}}_1^s,{\mathrm{\&gamma;}}_1^s),& {\mathrm{\&beta;}}_{s,i}~\mathrm{Gamma}({\mathrm{\&kappa;}}_1,{\mathrm{\&upsi;}}_1)\end{array}---\left(9\right)$

wherein, ${\mathrm{\&kappa;}}_1={10}^{-6},{\mathrm{\&upsi;}}_1={10}^{-6},[{\mathrm{\&epsiv;}}_0^r,{\mathrm{\&gamma;}}_0^r]=\left[\mathrm{0.9,0.1}\right]\&times;N_1,[{\mathrm{\&epsiv;}}_0^s,{\mathrm{\&gamma;}}_0^s,]=[1/n_p,1/n_p]\&times;N_s,$ N_sthe number of wavelet coefficients of L layers is expressed as the s-th 1.

(4) Sampling wavelet coefficient theta_s,ji. The wavelet coefficient theta is obtained through the steps (2) and (3)_s,jiLikelihood function of theta_s,jiAnd other prior distributions of unknown parameters. The wavelet coefficient theta to be estimated is obtained through Bayesian inference as follows:

p(Θ|G)∝∫p(G|Θ,β)p(Θ|Ω)p(Ω)p(β)dΩ（10）

wherein, $\mathrm{\&Omega;}=\{{\left\{{\mathrm{\&beta;}}_{s,i}\right\}}_{i=1:n_e}^{s=2:L},{\mathrm{\&pi;}}_r,{\{{\mathrm{\&pi;}}_0^s,{\mathrm{\&pi;}}_1^s\}}_{s=2:L}\},\mathrm{\&beta;}={[{\mathrm{\&beta;}}_1,{\mathrm{\&beta;}}_2,...,{{\mathrm{\&beta;}}_n}_e]}^T$ assuming that the elements within β and Ω are independently identically distributed for all noise accuracies, p (β) and p (Ω) are as follows:

$p\left(\mathrm{\&beta;}\right)=\underset{i=1}{\overset{n_e}{\mathrm{\&Pi;}}}\mathrm{Gamma}({\mathrm{\&kappa;}}_0,{\mathrm{\&upsi;}}_0),p\left(\mathrm{\&Omega;}\right)=\underset{i=1}{\overset{n_e}{\mathrm{\&Pi;}}}\left\{\begin{array}{c}\underset{s=1}{\overset{L}{\mathrm{\&Pi;}}}\mathrm{Gamma}({\mathrm{\&kappa;}}_1,{\mathrm{\&upsi;}}_1)\&times;\mathrm{Beta}({\mathrm{\&epsiv;}}_0^r,{\mathrm{\&gamma;}}_0^r)\\ \&times;\underset{s=1}{\overset{L}{\mathrm{\&Pi;}}}\mathrm{Beta}({\mathrm{\&epsiv;}}_0^s,{\mathrm{\&gamma;}}_0^s)\mathrm{Beta}({\mathrm{\&epsiv;}}_1^s,{\mathrm{\&gamma;}}_1^s)\end{array}\right\}---\left(11\right)$

a markov monte carlo method based on gibbs sampling is employed to approximate the a posteriori distribution of the variables to be estimated by extracting samples from the corresponding conditional distribution. Theta_s,jiThe gibbs sampler as follows:

$\begin{array}{c}p\left({\mathrm{\&theta;}}_{s,\mathrm{ji}}\right|~)=(1-{\mathrm{\&pi;}}_{s,\mathrm{ji}}){\mathrm{\&delta;}}_0+{\mathrm{\&pi;}}_{s,\mathrm{ji}}N({\mathrm{\&mu;}}_{s,\mathrm{ji}},{\mathrm{\&beta;}}_{s,\mathrm{ji}}^{-1})\\ {\mathrm{\&beta;}}_{s,\mathrm{ji}}={\mathrm{\&beta;}}_{s,i}+{\mathrm{\&beta;}}_i{\mathrm{\&Phi;}}_j^T{\mathrm{\&Phi;}}_j,{\mathrm{\&mu;}}_{\mathrm{sji}}={\mathrm{\&beta;}}_{s,\mathrm{ji}}^{-1}{\mathrm{\&Phi;}}_j^T[g_i-\underset{k\&NotEqual;j}{\overset{n_p}{\mathrm{\&Sigma;}}}{\mathrm{\&Phi;}}_k^T{\mathrm{\&Phi;}}_k]\end{array}---\left(12\right)$

wherein-represents a further parameter, Φ, required in the conditional distribution_jRepresents phi_jColumn j in (d).

(5) Sampling noise accuracy β_i. The sampling distribution satisfies the following gamma distribution:

$p\left({\mathrm{\&beta;}}_i\right|~)=\mathrm{Gamma}({{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_i,{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_i),{{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_i={\mathrm{\&kappa;}}_0+\frac{m}{2},{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_i={\mathrm{\&upsi;}}_0+\frac{1}{2}{(g_i-{\mathrm{\&Phi;\&theta;}}_i)}^T(g_i-{\mathrm{\&Phi;\&theta;}}_i)---\left(13\right)$

(6) sampling hyper-parametersThe specific sampling distribution is as follows:

$p\left({\mathrm{\&beta;}}_{s,i}\right|~)=\mathrm{Gamma}({{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_{s,i},{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_{s,i}),{{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_{s,i}={\mathrm{\&kappa;}}_1+\frac{1}{2}\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0),{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_{s,j}={\mathrm{\&upsi;}}_1+\frac{1}{2}\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_{s,\mathrm{ji}}^2---\left(14\right)$

$p\left({\mathrm{\&pi;}}_r\right|~)=\mathrm{Gamma}({{\mathrm{\&epsiv;}}_{\mathrm{\&pi;}}}_r,{{\mathrm{\&gamma;}}_{\mathrm{\&pi;}}}_r),{{\mathrm{\&epsiv;}}_{\mathrm{\&pi;}}}_r={\mathrm{\&epsiv;}}_0^r+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0),{{\mathrm{\&gamma;}}_{\mathrm{\&pi;}}}_r={\mathrm{\&gamma;}}_0^r+\frac{1}{2}\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}=0)---\left(15\right)$

$\begin{array}{c}p\left({\mathrm{\&pi;}}_{s,i}^0\right|~)=\mathrm{Gamma}({\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^0},{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^0})\\ {\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^0}={\mathrm{\&epsiv;}}_0^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}=0),{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^0}={\mathrm{\&gamma;}}_0^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}=0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}=0)\end{array}---\left(16\right)$

$\begin{array}{c}p\left({\mathrm{\&pi;}}_{s,i}^1\right|~)=\mathrm{Gamma}({\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^1},{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^1})\\ {\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^1}={\mathrm{\&epsiv;}}_1^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}\&NotEqual;0),{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^1}={\mathrm{\&gamma;}}_1^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}=0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}\&NotEqual;0)\end{array}---\left(17\right)$

wherein 1(x) is 1 when x is true, otherwise it is 0.

(7) And (4), (5) and (6) are circularly and sequentially executed for each wavelet coefficient in the theta to carry out sampling until a convergence condition or a termination condition is met, and finally estimated wavelet coefficients are obtainedM_sIn order to count the total number of samples,for the result of each sampling. The termination condition used is the number of iterations, where the sampler stabilization process is 400 times and the sampling process is 200 times.

4. Obtaining reconstructed hyperspectral data from the linear mixture model (1)

The invention has the beneficial effects that: according to the method, the inherent low-rank property of the hyperspectral data is explored through linear unmixing in the compression process, the abundance value is converted into a structured sparse signal through wavelet transformation, and then compressed data are obtained through compressed sensing. In the reconstruction process, a proper end member matrix is selected from a spectrum library, structured sparse prior of an abundance value wavelet coefficient is introduced, then a Bayesian inference method based on Gibbs sampling is used for accurately reconstructing the abundance value matrix, and finally a linear mixed model is used for reconstructing original data. The data of the invention is on the urban data shot by a HYDICE satellite, when the compression ratio is 1:16, the normalized mean square error is less than 0.4, and on the Indiana data shot by AVIRIS, when the compression ratio is 1:16, the normalized mean square error is also less than 0.4, and the precision is improved by about 10 percent compared with the precision of a compressed sensing algorithm in the background technology.

The present invention will be described in detail with reference to the following embodiments.

Detailed Description

The Bayes hyperspectral unmixing compressive sensing method based on the structured sparse prior specifically comprises the following steps:

for one high spectrum datan_pRepresenting the total number of pixels, n, contained in the hyperspectral data space_bRepresenting the number of bands, x, contained in the hyperspectral data_i,i＝1,...,n_bThe image representing each band is spread by rows to form a column vector. Each pixel is at n_bThe different reflection values over the individual wavelength bands constitute the spectrum of the pixel. Typically, the only spectrum that a pure substance has is called an end member. Due to factors such as low spatial resolution, the spectrum of the pixel is often a mixture of the spectra of various pure terrestrial objects. This spectral mixture phenomenon can be described using a linear mixture model or a non-linear mixture model. Only linear hybrid models are of interest herein. The model considers that any one mixed spectrum is a linear combination of all end members, and the proportion of the end members in the mixed spectrum is called abundance value. X can be expressed as a linear combination of an abundance value matrix and an end-member matrix, as follows:

X＝HW（1）

wherein,is a matrix of abundance values and is,represents the proportion of the ith end member in the mixed spectrum of each pixel, n_eThe number of end-members is indicated,is an end-member matrix in whichRepresenting each end member. Since adjacent pixels have similar abundance values for each end-member, the abundance value matrix inherits the similarity of spatially adjacent pixels. To tap signals with sparse structures from the abundance value matrix, a Discrete Wavelet Transform (DWT) is used for abundance values:

Θ＝Ψ^TH（2）

Y＝Ψ^TX＝Ψ^THW＝ΘW

wherein,is a wavelet basis and is a composite of a plurality of wavelets,wavelet coefficients representing abundance values, and Y represents wavelet coefficients of the entire data.Represents h_iHas structural sparsity. The method specifically includes the steps that sparse wavelet coefficients need to be reconstructed by using a compressed sensing technology for compressed data, then a fullness value matrix is reconstructed by using the wavelet coefficients, and finally the reconstruction of the whole data is achieved, wherein the method specifically includes the following steps:

1. and obtaining compressed data.

Generation of m × n using random variables satisfying a Gaussian distribution_pA matrix of sizes, and each column is normalized to obtain a random perceptual matrixThe wavelet transformed raw data is sampled using phi,obtaining compressed dataThe following were used:

F＝ΦY＝ΦΨ^TX（3）

wherein m represents a pair length of n_pAfter compression of the signal, typically m < n_p。

2. And selecting end members.

For a scale-limited scenario, the influence of environmental factors on the spectrum is ignored, assuming a limited number of end-members. The ASTER spectral library is used, and partial spectrums are selectively extracted from the spectral library as end members according to specific scenes. Finally determining that n is contained_eA matrix W of individual end-members. Usually, only a small number of end members are included in the fixed scene, so that the low rank characteristic of the hyperspectral data can be mined by linear unmixing.

3. And reconstructing an abundance value matrix through Bayesian inference.

For accurate reconstruction of the raw data, bayesian inference is employed herein to estimate a matrix of abundance values from the compressed data. The specific implementation process comprises the following steps:

(1) obtaining data in which wavelet coefficients of an abundance value matrix are individually compressed by using the least square method (Leastquares) by equations (2) and (3)

G＝FW^T(WW^T)^-1＝ΦΘ+N（4）

Wherein,representing the noise term introduced using the least squares method.

(2) A gaussian likelihood function of the wavelet coefficients is constructed. Obtaining a vector form g of the formula (4)_i＝Φθ_i+n_iAssuming noise n_iObedience mean 0 and covarianceGaussian distribution ofAnd the noise vectors are independently and identically distributed, theta can be obtained_i,i＝1,2,...,n_eLikelihood function of (d):

$p\left(g_i\right|{\mathrm{\&theta;}}_i,{\mathrm{\&beta;}}_i^{-1})={\left(2\mathrm{\&pi;}{\mathrm{\&beta;}}_i^{-1}\right)}^{-\frac{m}{2}}\exp (-\frac{1}{2{\mathrm{\&beta;}}_i^{-1}}{\left|\right|g_i-{\mathrm{\&Phi;\&theta;}}_i\left|\right|}_2^2)---\left(5\right)$

β_ifor the accuracy of the noise vector, β is specified herein_iIs a gamma distribution as follows:

p(β_i|κ₀,υ₀)＝Gamma(κ₀,υ₀)（6）

κ₀＝10^-6,υ₀＝10^-6the shape parameter and the scale parameter of the gamma distribution.

(3) Establishing theta_iThe structured sparsity of each element in (a). The structured sparse prior is mainly reflected in a tree structure of wavelet coefficients, when a parent node is 0, a child node thereof is highly likely to be 0, and vice versa:

$p\left({\mathrm{\&theta;}}_{s,\mathrm{ji}}\right|{\mathrm{\&beta;}}_{s,i}^{-1})=(1-{\mathrm{\&pi;}}_{s,\mathrm{ji}}){\mathrm{\&delta;}}_0+{\mathrm{\&pi;}}_{s,\mathrm{ji}}N(0,{\mathrm{\&beta;}}_{s,i}^{-1})---\left(7\right)$

wherein, s represents the total number of scales of the wavelet transform, and θ represents the total number of scales of the wavelet transform, wherein s is 1_s,jiDenotes theta_iThe jth element at scale s.₀Is a distribution that is completely centered at 0. Pi_s,jiIs a proportionality coefficient when_s,jiTowards 0, theta_s,jiIs 0, otherwise, pi_s,jiWhen going to 1, theta_s,jiFrom a gaussian distribution. Pi_s,jiDetermined by the tree structure of the wavelet transform:

${\mathrm{\&pi;}}_{s,\mathrm{ji}}=\left\{\begin{array}{cc}{\mathrm{\&pi;}}_r& \mathrm{ifs}=1\\ {\mathrm{\&pi;}}_s^0& \mathrm{if}2\&le;s\&le;L,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}=0\\ {\mathrm{\&pi;}}_s^1& \mathrm{if}2\&le;s\&le;L,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}\&NotEqual;0\end{array}\right.---\left(8\right)$

θ_pa(s,ji)is theta_s,jiParent node of when theta_pa(s,ji)When equal to 0, pi_s,jiSet to a high probability of 0When theta is_pa(s,ji)Not equal to 0, pi_s,jiSet to a high probability of 1β_s,iIs theta_s,jiEssence of a prior medium Gaussian distributionAnd (4) degree. Pi_r，And β_s,iThe prior distribution of (a) is as follows:

$\begin{array}{cc}{\mathrm{\&pi;}}_r~\mathrm{Beta}({\mathrm{\&epsiv;}}_0^r,{\mathrm{\&gamma;}}_0^r),& {\mathrm{\&pi;}}_s^0~\mathrm{Beta}({\mathrm{\&epsiv;}}_0^s,{\mathrm{\&gamma;}}_0^s),\\ {\mathrm{\&pi;}}_s^1~\mathrm{Beta}({\mathrm{\&epsiv;}}_1^s,{\mathrm{\&gamma;}}_1^s),& {\mathrm{\&beta;}}_{s,i}~\mathrm{Gamma}({\mathrm{\&kappa;}}_1,{\mathrm{\&upsi;}}_1)\end{array}---\left(9\right)$

wherein, ${\mathrm{\&kappa;}}_1={10}^{-6},{\mathrm{\&upsi;}}_1={10}^{-6},[{\mathrm{\&epsiv;}}_0^r,{\mathrm{\&gamma;}}_0^r]=\left[\mathrm{0.9,0.1}\right]\&times;N_1,[{\mathrm{\&epsiv;}}_0^s,{\mathrm{\&gamma;}}_0^s,]=[1/n_p,1/n_p]\&times;N_s,$ N_sthe number of wavelet coefficients of L layers is expressed as the s-th 1.

(4) Sampling wavelet coefficient theta_s,ji. The wavelet coefficient theta is obtained through the steps (2) and (3)_s,jiLikelihood function of theta_s,jiAnd other prior distributions of unknown parameters. The wavelet coefficient Θ to be estimated can be obtained by bayesian inference as follows:

p(Θ|G)∝∫p(G|Θ,β)p(Θ|Ω)p(Ω)p(β)dΩ（10）

wherein, $\mathrm{\&Omega;}=\{{\left\{{\mathrm{\&beta;}}_{s,i}\right\}}_{i=1:n_e}^{s=2:L},{\mathrm{\&pi;}}_r,{\{{\mathrm{\&pi;}}_0^s,{\mathrm{\&pi;}}_1^s\}}_{s=2:L}\},\mathrm{\&beta;}={[{\mathrm{\&beta;}}_1,{\mathrm{\&beta;}}_2,...,{{\mathrm{\&beta;}}_n}_e]}^T$ assuming that the elements within β and Ω are independently identically distributed for all noise accuracies, p (β) and p (Ω) are as follows:

$p\left(\mathrm{\&beta;}\right)=\underset{i=1}{\overset{n_e}{\mathrm{\&Pi;}}}\mathrm{Gamma}({\mathrm{\&kappa;}}_0,{\mathrm{\&upsi;}}_0),p\left(\mathrm{\&Omega;}\right)=\underset{i=1}{\overset{n_e}{\mathrm{\&Pi;}}}\left\{\begin{array}{c}\underset{s=1}{\overset{L}{\mathrm{\&Pi;}}}\mathrm{Gamma}({\mathrm{\&kappa;}}_1,{\mathrm{\&upsi;}}_1)\&times;\mathrm{Beta}({\mathrm{\&epsiv;}}_0^r,{\mathrm{\&gamma;}}_0^r)\\ \&times;\underset{s=1}{\overset{L}{\mathrm{\&Pi;}}}\mathrm{Beta}({\mathrm{\&epsiv;}}_0^s,{\mathrm{\&gamma;}}_0^s)\mathrm{Beta}({\mathrm{\&epsiv;}}_1^s,{\mathrm{\&gamma;}}_1^s)\end{array}\right\}---\left(11\right)$

since equation (10) is difficult to solve, the markov monte carlo method based on gibbs sampling is employed herein. The a posteriori distribution of the desired estimated variable is approximated by taking samples from the corresponding conditional distribution. Theta_s,jiThe gibbs sampler as follows:

$\begin{array}{c}p\left({\mathrm{\&theta;}}_{s,\mathrm{ji}}\right|~)=(1-{\mathrm{\&pi;}}_{s,\mathrm{ji}}){\mathrm{\&delta;}}_0+{\mathrm{\&pi;}}_{s,\mathrm{ji}}N({\mathrm{\&mu;}}_{s,\mathrm{ji}},{\mathrm{\&beta;}}_{s,\mathrm{ji}}^{-1})\\ {\mathrm{\&beta;}}_{s,\mathrm{ji}}={\mathrm{\&beta;}}_{s,i}+{\mathrm{\&beta;}}_i{\mathrm{\&Phi;}}_j^T{\mathrm{\&Phi;}}_j,{\mathrm{\&mu;}}_{\mathrm{sji}}={\mathrm{\&beta;}}_{s,\mathrm{ji}}^{-1}{\mathrm{\&Phi;}}_j^T[g_i-\underset{k\&NotEqual;j}{\overset{n_p}{\mathrm{\&Sigma;}}}{\mathrm{\&Phi;}}_k^T{\mathrm{\&Phi;}}_k]\end{array}---\left(12\right)$

wherein-represents a further parameter, Φ, required in the conditional distribution_jRepresents phi_jColumn j in (d).

(5) Sampling noise accuracy β_i. The sampling distribution satisfies the following gamma distribution:

$p\left({\mathrm{\&beta;}}_i\right|~)=\mathrm{Gamma}({{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_i,{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_i),{{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_i={\mathrm{\&kappa;}}_0+\frac{m}{2},{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_i={\mathrm{\&upsi;}}_0+\frac{1}{2}{(g_i-{\mathrm{\&Phi;\&theta;}}_i)}^T(g_i-{\mathrm{\&Phi;\&theta;}}_i)---\left(13\right)$

(6) sampling hyper-parametersThe specific sampling distribution is as follows:

$p\left({\mathrm{\&beta;}}_{s,i}\right|~)=\mathrm{Gamma}({{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_{s,i},{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_{s,i}),{{\mathrm{\&kappa;}}_{\mathrm{\&beta;}}}_{s,i}={\mathrm{\&kappa;}}_1+\frac{1}{2}\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0),{{\mathrm{\&upsi;}}_{\mathrm{\&beta;}}}_{s,j}={\mathrm{\&upsi;}}_1+\frac{1}{2}\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_{s,\mathrm{ji}}^2---\left(14\right)$

$p\left({\mathrm{\&pi;}}_r\right|~)=\mathrm{Gamma}({{\mathrm{\&epsiv;}}_{\mathrm{\&pi;}}}_r,{{\mathrm{\&gamma;}}_{\mathrm{\&pi;}}}_r),{{\mathrm{\&epsiv;}}_{\mathrm{\&pi;}}}_r={\mathrm{\&epsiv;}}_0^r+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0),{{\mathrm{\&gamma;}}_{\mathrm{\&pi;}}}_r={\mathrm{\&gamma;}}_0^r+\frac{1}{2}\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}=0)---\left(15\right)$

$\begin{array}{c}p\left({\mathrm{\&pi;}}_{s,i}^0\right|~)=\mathrm{Gamma}({\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^0},{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^0})\\ {\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^0}={\mathrm{\&epsiv;}}_0^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}=0),{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^0}={\mathrm{\&gamma;}}_0^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}=0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}=0)\end{array}---\left(16\right)$

$\begin{array}{c}p\left({\mathrm{\&pi;}}_{s,i}^1\right|~)=\mathrm{Gamma}({\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^1},{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^1})\\ {\mathrm{\&epsiv;}}_{{\mathrm{\&pi;}}_{s,i}^1}={\mathrm{\&epsiv;}}_1^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}\&NotEqual;0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}\&NotEqual;0),{\mathrm{\&gamma;}}_{{\mathrm{\&pi;}}_{s,i}^1}={\mathrm{\&gamma;}}_1^s+\underset{j=1}{\overset{N_s}{\mathrm{\&Sigma;}}}1({\mathrm{\&theta;}}_{s,\mathrm{ji}}=0,{\mathrm{\&theta;}}_{\mathrm{pa}(s,\mathrm{ji})}\&NotEqual;0)\end{array}---\left(17\right)$

wherein 1(x) is 1 when x is true, otherwise it is 0.

(7) And (4), (5) and (6) are circularly and sequentially executed for each wavelet coefficient in the theta to carry out sampling until a convergence condition or a termination condition is met, and finally estimated wavelet coefficients are obtainedM_sIn order to count the total number of samples,for the result of each sampling. The termination condition used herein is the number of iterations, where the sampler stabilizes 400 times and the sampling 200 times.

4. And (5) reconstructing hyperspectral data.

Obtaining reconstructed hyperspectral data from the linear mixture model (1)

Claims (1)

1. A Bayes hyperspectral unmixing compressive sensing method based on structured sparse prior is characterized by comprising the following steps:

step one, aiming at high spectrum datan_pRepresenting the total number of pixels, n, contained in the hyperspectral data space_bRepresenting the number of bands, x, contained in the hyperspectral data_i,i＝1,...,n_bThe image representing each band is line-wiseUnfolding the formed column vector; each pixel is at n_bThe different reflection values on each band constitute the spectrum of the pixel; the spectrum is described by a linear mixed model; the linear mixed model considers that any mixed spectrum is a linear combination of all end members, and the proportion of the end members in the mixed spectrum is called a abundance value; x represents a linear combination of the abundance value matrix and the end-member matrix:

X＝HW(1)

wherein,is a matrix of abundance values and is,represents the proportion of the ith end member in the mixed spectrum of each pixel, n_eThe number of end-members is indicated,is an end-member matrix in whichRepresenting each end member;

using a discrete wavelet transform for the abundance value:

wherein,is a wavelet basis and is a composite of a plurality of wavelets,wavelet coefficients representing abundance values, Y represents wavelet coefficients of the entire data;represents h_iHas structural sparsity;

secondly, reconstructing an abundance value matrix by utilizing the wavelet coefficient, and finally realizing the reconstruction of the whole data:

step 1, generating m × n by using random variables satisfying Gaussian distribution_pA matrix of sizes, and each column is normalized to obtain a random perceptual matrixSampling the wavelet transformed original data by phi to obtain compressed data

F＝ΦY＝ΦΨ^TX(3)

Wherein m represents a pair length of n_pAfter compression of the signal, typically m < n_p；

Step 2, for a scene with limited scale, neglecting the influence of environmental factors on the spectrum, and assuming that the number of end members is limited; using an ASTER spectral library, and selectively extracting partial spectrums from the spectral library as end members aiming at a specific scene; finally determining that n is contained_eA matrix W of individual end members;

step 3, estimating an abundance value matrix from the compressed data by adopting a Bayesian inference mode;

(1) obtaining data obtained by compressing wavelet coefficients of an abundance value matrix independently by using a least square method through formula (2) and formula (3)

G＝FW^T(WW^T)^-1＝ΦΘ+N(4)

Wherein,representing a noise term introduced using a least squares method;

(2) obtaining a vector form g of the formula (4)_i＝Φθ_i+n_iAssuming noise n_iObedience mean 0 and covarianceGaussian distribution ofAnd the noise vectors are independently and identically distributed, then theta is obtained_i,i＝1,2,...,n_eLikelihood function of (d):

β_ifor the accuracy of the noise vector, β is specified_iIs a gamma distribution as follows:

p(β_i|κ₀,υ₀)＝Gamma(κ₀,υ₀)(6)

κ₀＝10^-6,υ₀＝10^-6is the shape parameter and the scale parameter of the gamma distribution;

(3) establishing theta_iStructured sparse prior of each element in (a); the structured sparse prior is mainly reflected in a tree structure of wavelet coefficients, when a parent node is 0, a child node thereof is highly likely to be 0, and vice versa:

wherein, s represents the total number of scales of the wavelet transform, and θ represents the total number of scales of the wavelet transform, wherein s is 1_s,jiDenotes theta_iThe jth element at scale s;₀is a distribution that is completely centered at 0; pi_s,jiIs a proportionality coefficient when_s,jiTowards 0, theta_s,jiIs 0, otherwise, pi_s,jiWhen going to 1, theta_s,jiFrom a gaussian distribution; pi_s,jiDetermined by the tree structure of the wavelet transform:

θ_pa(s,ji)is theta_s,jiParent node of when theta_pa(s,ji)When equal to 0, pi_s,jiSet to a high probability of 0When theta is_pa(s,ji)Not equal to 0, pi_s,jiSet to a high probability of 1β_s,iIs theta_s,jiThe accuracy of the gaussian distribution in the prior; pi_r， And β_s,iThe prior distribution of (a) is as follows:

wherein, κ₁＝10-⁶,υ₁＝10^-6, N_sThe number of wavelet coefficients of the L layer is expressed as the s-th 1.;

(4) sampling wavelet coefficient theta_s,ji(ii) a The wavelet coefficient theta is obtained through the steps (2) and (3)_s,jiLikelihood function of theta_s,jiThe structured sparse prior of (a), and prior distribution of other unknown parameters; the wavelet coefficient theta to be estimated is obtained by Bayesian extrapolation as followsObtaining:

p(Θ|G)∝∫p(G|Θ,β)p(Θ|Ω)p(Ω)p(β)dΩ(10)

wherein,assuming that the elements within β and Ω are independently identically distributed for all noise accuracies, p (β) and p (Ω) are as follows:

adopting a Markov Monte Carlo method based on Gibbs sampling, and approximating the posterior distribution of the variables to be estimated by extracting samples from corresponding condition distribution; theta_s,jiThe gibbs sampler as follows:

wherein-represents a further parameter, Φ, required in the conditional distribution_jRepresents phi_jJ-th column in (1);

(5) sampling noise accuracy β_i(ii) a The sampling distribution satisfies the following gamma distribution:

(6) sampling hyper-parametersThe specific sampling distribution is as follows:

wherein 1(x) is 1 when x is true, otherwise 0;

(7) and (4), (5) and (6) are circularly and sequentially executed for each wavelet coefficient in the theta to carry out sampling until a convergence condition or a termination condition is met, and finally estimated wavelet coefficients are obtainedM_sIn order to count the total number of samples,for each sampled result; the used termination condition is iteration times, wherein the sampler is stabilized for 400 times, and the sampling process is 200 times;

step 4, obtaining the reconstructed hyperspectral data according to the linear mixed model (1)