CN103745487A - Bayes high-spectral unmixing compressive sensing method based on structured sparsity prior - Google Patents

Abstract

The invention discloses a Bayes high-spectral unmixing compressive sensing method based on a structured sparsity prior and aims at solving a technical problem that prior high-spectral image compressive sensing method which unites with spectral unmixing is poor in precision. The technical scheme is that in a compressive process, inherent low-rank property of high-spectral data is explored through linear unmixing and wavelet transformation is used to convert an abundance value into a structured sparse signal and then compressive sensing is used to obtain compressive data; and in a reconstruction process, an appropriate end-member matrix is selected from a spectral library and the structured sparse prior of an abundance-value wavelet factor is introduced and then a Bayes inference method based on Gibbs sampling is used to establish an abundance-value matrix precisely and at last a linear mixing model is used to reconstruct original data. Compared with a background technology compressive-sensing-type algorithm, the precision of the Bayes high-spectral unmixing compressive sensing method based on the structured sparsity prior is improved by about 10%.

Description

The mixed compression sensing method of the high spectrum solution of Bayes based on structural sparse priori

Technical field

The present invention relates to the mixed compression sensing method of a kind of high spectrum solution, particularly relate to the mixed compression sensing method of the high spectrum solution of a kind of Bayes based on structural sparse priori.

Background technology

Because high spectrum image data volume is huge, the compression algorithm of high compression ratio, has very profound significance for collection, transmission and the processing of high spectrum image.Current existing Hyperspectral image compression algorithm is mainly divided into two classes, one class is traditional compression method based on information coding, these class methods are to expand from the compression method of normal image, laying particular emphasis on the redundancy of removing in high spectrum image between each wave band inside and wave band compresses, comprising Impulsive Difference coding (the Cluster-based differential pulse-code modulation of cluster, CDPCM), 3 D wavelet transformation (3dimension discrete wavelet transformation, 3D-DWT), Three-dimensional DCT (3dimension discrete cosine transformation, 3D-DCT) etc., but this class compression algorithm can only be compressed after Image Acquisition, still need a large amount of hardware to gather and temporarily store data, and ratio of compression is low conventionally, another kind of method is the compression method based on compressed sensing (Compressive Sensing), these class methods can be by gathering the original sparse signal of a small amount of sample point Exact Reconstruction, the hardware consumption needing for minimizing image data provides theoretical foundation, and the method more lays particular emphasis on process of reconstruction.

Document " A compressive sensing and unmixing scheme for hyperspectral data processing; IEEE Transactions on Image Processing; 2012,21 (3): 1200 – 1210 " discloses the mixed Compression of hyperspectral images cognitive method of a kind of combined spectral solution.The method is mixed high-spectral data is decomposed into Abundances matrix and end member matrix by linear solution; The gradient of Abundances matrix on Spatial Dimension has sparse property, can adopt the mode of compressed sensing to compress reconstruction.Specific practice is first to use compressed sensing matrix to obtain packed data; In conjunction with library of spectra, select suitable end member; Use afterwards the method for compressed sensing to rebuild the Abundances matrix that comprises sparse gradient; Finally, use linear mixed model to rebuild the high-spectral data that needs collection.But the sparse property of gradient can not be excavated the sparse property of the deeper structure of Abundances matrix, and this article uses the method for augmentation Lagrange to rebuild Abundances, the proportion coefficients that the method can not the sparse property of adaptively selected gradient regular terms.These two factors have directly affected the reconstruction precision of this compression algorithm.

Summary of the invention

In order to overcome the deficiency of the mixed Compression of hyperspectral images cognitive method low precision of existing combined spectral solution, the invention provides the mixed compression sensing method of the high spectrum solution of a kind of Bayes based on structural sparse priori.The method compression process is mixed the low-rank character of excavating high-spectral data inherence by linear solution, use wavelet transformation to convert Abundances to structural sparse signal, uses afterwards compressed sensing to obtain packed data.Process of reconstruction, from library of spectra, select suitable end member matrix, introduce the structural sparse priori of Abundances wavelet coefficient, then use the Bayesian inference method Exact Reconstruction Abundances matrix based on gibbs sampler, finally use linear mixed model to rebuild raw data.In the urban data that data of the present invention are taken at HYDICE satellite, when ratio of compression is 1:16, normalized square error is less than 0.4, in the Indiana data of taking at AVIRIS, when ratio of compression is 1:16, normalized mean squared error is less than 0.4 equally, with respect to background technology compressed sensing class arithmetic accuracy, promotes 10% left and right.

The technical solution adopted for the present invention to solve the technical problems is: the high spectrum solution of a kind of Bayes based on structural sparse priori is mixed compression sensing method, is characterized in comprising the following steps:

Step 1, for a panel height spectroscopic data

n _prepresent the total number of pixels comprising on high-spectral data space, n _brepresent the wave band number that high-spectral data comprises, x _i, i=1 ..., n _bthe image that represents each wave band launches the column vector forming by row.Each pixel is at n _bdifferent reflected values on individual wave band have formed the spectrum of this pixel.This spectrum adopts linear mixed model to describe.Linear mixed model thinks that any one mixed spectra is the linear combination of all end members, and end member shared ratio in mixed spectra is called Abundances.X represents the linear combination of Abundances matrix and end member matrix:

X＝HW （1）

Wherein,

for Abundances matrix,

represent i the ratio that end member is shared in the mixed spectra of each pixel, n _erepresent the number of end member,

for end member matrix, wherein

represent each end member.

For Abundances, use wavelet transform:

Θ＝Ψ ^TH （2）

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

Wherein, for wavelet basis,

the wavelet coefficient that represents Abundances, Y represents the wavelet coefficient of whole data.

represent h _iwavelet coefficient, there is the sparse property of structure.

Step 2, utilize wavelet coefficient to rebuild Abundances matrix, finally realize the reconstruction of whole data:

1, use the stochastic variable that meets Gaussian distribution to produce m × n _pthe matrix of size, and each row is normalized and obtains random perception matrix

raw data after using Φ to wavelet transformation is sampled, and obtains packed data

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

Wherein, m represents that to length be n _psignal compression after data length, general m < n _p.

2, for the limited scene of yardstick, ignore the impact of environmental factor on spectrum, suppose that end member is Limited Number.Use ASTER library of spectra, for specific scene selectively from library of spectra extracting part divide spectrum as end member.Final definite n that comprises _ethe matrix W of individual end member.

3, the data estimation Abundances matrix of the mode of employing Bayesian inference from compression.

(1), by formula (2) and formula (3), utilize the wavelet coefficient of least square method acquisition Abundances matrix by the data after compressing separately

G=FW ^t(WW ^t) ^-1=Φ Θ+N (4) wherein,

the noise item that represents to utilize least square method and introduce.

(2) obtain the vector form g of formula (4) _i=Φ θ _i+ n _i, suppose noise n _iobeying average is 0, and covariance is

gaussian distribution

and independent same distribution between each noise vector, obtain θ _i, i=1,2 ..., n _elikelihood function:

$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 precision of noise vector, specify β _ipriori be that following gamma distributes:

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

κ ₀=10 ^-6, υ ₀=10 ^-6for form parameter and the scale parameter of gamma distribution.

(3) set up θ _iin the structural sparse priori of each element.Structural sparse priori is to be mainly reflected in the tree construction of wavelet coefficient, and when parent node is 0, its child node is also probably 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 yardstick s=1 in wavelet transformation ..., L, L is the total number of yardstick, θ _{s, ji}represent θ _ij the element at yardstick s place.δ ₀it is a distribution that concentrates on 0 place completely.π _{s, ji}be a scale-up factor, work as π _{s, ji}be tending towards at 0 o'clock, θ _{s, ji}0, otherwise, π _{s, ji}be tending towards at 1 o'clock, θ _{s, ji}come from a Gaussian distribution.π _{s, ji}tree construction by wavelet transformation is determined:

${\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)}θ _{s, ji}parent node, work as θ _{pa (s, ji)}=0 o'clock, π _{s, ji}be made as high probability and be 0

work as θ _{pa (s, ji)}≠ 0 o'clock, π _{s, ji}be made as high probability and be 1

β _s,iθ _{s, ji}the precision of Gaussian distribution in priori.π _r, and β _s,iprior distribution 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 _srepresent wavelet transformation s=1 ..., the number of the wavelet coefficient of L layer.

(4) sampling wavelet coefficient θ _{s, ji}.By above-mentioned (2), (3) step, obtained wavelet coefficient θ _{s, ji}likelihood function, θ _{s, ji}structural sparse priori, and the prior distribution of other unknown parameters.Need the wavelet coefficient Θ estimating, by following Bayesian inference, obtain:

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$ For all noise precision, suppose the element independent same distribution in β and Ω, p (β) and p (Ω) are as follows respectively:

$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)$

The Markov monte carlo method of employing based on gibbs sampler, needs the posteriority of predictor to distribute by being similar to from corresponding condition distribution sample drawn.θ _{s, ji}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, other parameters that～expression condition needs in distributing, Φ _jrepresent Φ _jin j row.

(5) sampling noiset precision β _i.Sample distribution meets following gamma and distributes:

$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) the super parameter of sampling

concrete sample 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 at x when being genuine, otherwise is 0.

(7) for each wavelet coefficient circulation in Θ successively execution (4), (5), (6) three steps, sample, until meet the condition of convergence or end condition, obtain the final wavelet coefficient of estimating

m _sfor total number of sampling,

for the result of sampling each time.The end condition using is iterations, wherein, and sampling thief stabilization process 400 times, sampling process 200 times.

4, according to linear mixed model (1), obtain the high-spectral data of rebuilding

The invention has the beneficial effects as follows: the method compression process is mixed the low-rank character of excavating high-spectral data inherence by linear solution, use wavelet transformation to convert Abundances to structural sparse signal, use afterwards compressed sensing to obtain packed data.Process of reconstruction, from library of spectra, select suitable end member matrix, introduce the structural sparse priori of Abundances wavelet coefficient, then use the Bayesian inference method Exact Reconstruction Abundances matrix based on gibbs sampler, finally use linear mixed model to rebuild raw data.In the urban data that data of the present invention are taken at HYDICE satellite, when ratio of compression is 1:16, normalized square error is less than 0.4, in the Indiana data of taking at AVIRIS, when ratio of compression is 1:16, normalized mean squared error is less than 0.4 equally, with respect to background technology compressed sensing class arithmetic accuracy, promotes 10% left and right.

Below in conjunction with embodiment, the present invention is elaborated.

Embodiment

The mixed compression sensing method concrete steps of the high spectrum solution of Bayes that the present invention is based on structural sparse priori are as follows:

For a panel height spectroscopic data n _prepresent the total number of pixels comprising on high-spectral data space, n _brepresent the wave band number that high-spectral data comprises, x _i, i=1 ..., n _bthe image that represents each wave band launches the column vector forming by row.Each pixel is at n _bdifferent reflected values on individual wave band have formed the spectrum of this pixel.Conventionally, unique spectrum that pure material has is known as end member.Due to factor impacts such as low spatial resolutions, the spectrum of pixel is the mixing of multiple pure object spectrum often.This spectral mixing phenomenon can be described with linear mixed model or non-linear mixture model.Only pay close attention to linear mixed model herein.This model thinks that any one mixed spectra is the linear combination of all end members, and end member shared ratio in mixed spectra is called Abundances.X can be expressed as the linear combination of Abundances matrix and end member matrix, as follows:

X=HW (1) wherein, for Abundances matrix,

represent i the ratio that end member is shared in the mixed spectra of each pixel, n _erepresent the number of end member,

for end member matrix, wherein

represent each end member.Because adjacent pixel has approximate Abundances to each end member, therefore, Abundances matrix has been inherited the similarity of neighbor on space.In order excavating, to there is the sparse signal of structure from abundance value matrix, for Abundances, to use wavelet transform (DWT):

Θ＝Ψ ^TH （2）

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

Wherein, for wavelet basis,

the wavelet coefficient that represents Abundances, Y represents the wavelet coefficient of whole data.

represent h _iwavelet coefficient, there is the sparse property of structure.Need to use the wavelet coefficient that compressed sensing reconstruction is sparse for packed data exactly herein, utilize afterwards wavelet coefficient to rebuild Abundances matrix, finally realize the reconstruction of whole data, concrete steps are as follows:

1, packed data obtains.

The stochastic variable that use meets Gaussian distribution produces m × n _pthe matrix of size, and each row is normalized and obtains random perception matrix

raw data after using Φ to wavelet transformation is sampled, and obtains packed data as follows:

F=Φ Y=Φ Ψ ^twherein, m represents that to length be n to X (3) _psignal compression after data length, general m < n _p.

2, end member is selected.

For the limited scene of yardstick, ignore the impact of environmental factor on spectrum, suppose that end member is Limited Number.Use ASTER library of spectra herein, for specific scene selectively from library of spectra extracting part divide spectrum as end member.Final definite n that comprises _ethe matrix W of individual end member.Conventionally in fixed scene, only comprise a small amount of end member, therefore by linear solution, mix the low-rank characteristic that can excavate high-spectral data itself.

3, Bayesian inference is rebuild Abundances matrix.

For the Exact Reconstruction of raw data, adopt the mode of Bayesian inference to carry out the data estimation Abundances matrix from compression herein.Specific implementation process is divided into following step:

(1), by formula (2) and (3), utilize the wavelet coefficient of least square method (Least squares) acquisition Abundances matrix by the data after compressing separately

G=FW ^t(WW ^t) ^-1=Φ Θ+N (4) wherein, the noise item that represents to utilize least square method and introduce.

(2) Gauss's likelihood function of structure wavelet coefficient.Obtain the vector form g of formula (4) _i=Φ θ _i+ n _i, suppose noise n _iobeying average is 0, and covariance is

gaussian distribution

and independent same distribution between each noise vector, can obtain θ _i, i=1,2 ..., n _elikelihood function:

$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 precision of noise vector, specify β herein _ipriori be that following gamma distributes:

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

κ ₀=10 ^-6, υ ₀=10 ^-6for form parameter and the scale parameter of gamma distribution.

(3) set up θ _iin the structural sparse priori of each element.Structural sparse priori is to be mainly reflected in the tree construction of wavelet coefficient, and when parent node is 0, its child node is also probably 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 yardstick s=1 in wavelet transformation ..., L, L is the total number of yardstick, θ _{s, ji}represent θ _ij the element at yardstick s place.δ ₀it is a distribution that concentrates on 0 place completely.π _{s, ji}be a scale-up factor, work as π _{s, ji}be tending towards at 0 o'clock, θ _{s, ji}0, otherwise, π _{s, ji}be tending towards at 1 o'clock, θ _{s, ji}come from a Gaussian distribution.π _{s, ji}tree construction by wavelet transformation is determined:

${\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)}θ _{s, ji}parent node, work as θ _{pa (s, ji)}=0 o'clock, π _{s, ji}be made as high probability and be 0

work as θ _{pa (s, ji)}≠ 0 o'clock, π _{s, ji}be made as high probability and be 1

β _s,iθ _{s, ji}the precision of Gaussian distribution in priori.π _r,

and β _s,iprior distribution 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 _srepresent wavelet transformation s=1 ..., the number of the wavelet coefficient of L layer.

(4) sampling wavelet coefficient θ _{s, ji}.By above-mentioned (2), (3) step, obtained wavelet coefficient θ _{s, ji}likelihood function, θ _{s, ji}structural sparse priori, and the prior distribution of other unknown parameters.Need the wavelet coefficient Θ estimating, can obtain by following Bayesian inference:

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$ For all noise precision, suppose the element independent same distribution in β and Ω, p (β) and p (Ω) are as follows respectively:

$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)$

Because formula (10) is difficult to solve, so adopt the Markov monte carlo method based on gibbs sampler herein.By being similar to from corresponding condition distribution sample drawn, need the posteriority of predictor to distribute.θ _{s, ji}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, other parameters that～expression condition needs in distributing, Φ _jrepresent Φ _jin j row.

(5) sampling noiset precision β _i.Sample distribution meets following gamma and distributes:

$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) the super parameter of sampling

concrete sample 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 at x when being genuine, otherwise is 0.

(7) for each wavelet coefficient circulation in Θ successively execution (4), (5), (6) three steps, sample, until meet the condition of convergence or end condition, obtain the final wavelet coefficient of estimating

m _sfor total number of sampling,

for the result of sampling each time.End condition used herein is iterations, wherein, and sampling thief stabilization process 400 times, sampling process 200 times.

4, rebuild high-spectral data.

According to linear mixed model (1), obtain the high-spectral data of rebuilding

Claims (1)

1. the mixed compression sensing method of the high spectrum solution of the Bayes based on structural sparse priori, is characterized in that comprising the following steps:

Step 1, for a panel height spectroscopic data

n _prepresent the total number of pixels comprising on high-spectral data space, n _brepresent the wave band number that high-spectral data comprises, x _i, i=1 ..., n _bthe image that represents each wave band launches the column vector forming by row; Each pixel is at n _bdifferent reflected values on individual wave band have formed the spectrum of this pixel; This spectrum adopts linear mixed model to describe; Linear mixed model thinks that any one mixed spectra is the linear combination of all end members, and end member shared ratio in mixed spectra is called Abundances; X represents the linear combination of Abundances matrix and end member matrix:

X＝HW （1）

Wherein,

for Abundances matrix,

represent i the ratio that end member is shared in the mixed spectra of each pixel, n _erepresent the number of end member,

for end member matrix, wherein

represent each end member;

For Abundances, use wavelet transform:

Θ＝Ψ ^TH （2）

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

Wherein,

for wavelet basis,

the wavelet coefficient that represents Abundances, Y represents the wavelet coefficient of whole data; represent h _iwavelet coefficient, there is the sparse property of structure;

Step 2, utilize wavelet coefficient to rebuild Abundances matrix, finally realize the reconstruction of whole data:

The stochastic variable that step 1, use meet Gaussian distribution produces m × n _pthe matrix of size, and each row is normalized and obtains random perception matrix

raw data after using Φ to wavelet transformation is sampled, and obtains packed data

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

Wherein, m represents that to length be n _psignal compression after data length, general m < n _p;

Step 2, for the limited scene of yardstick, ignore the impact of environmental factor on spectrum, suppose that end member is Limited Number; Use ASTER library of spectra, for specific scene selectively from library of spectra extracting part divide spectrum as end member; Final definite n that comprises _ethe matrix W of individual end member;

The data estimation Abundances matrix of the mode of step 3, employing Bayesian inference from compression;

(1), by formula (2) and formula (3), utilize the wavelet coefficient of least square method acquisition Abundances matrix by the data after compressing separately

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

Wherein, the noise item that represents to utilize least square method and introduce;

(2) obtain the vector form g of formula (4) _i=Φ θ _i+ n _i, suppose noise n _iobeying average is 0, and covariance is

gaussian distribution

and independent same distribution between each noise vector, obtain θ _i, i=1,2 ..., n _elikelihood function:

β _ifor the precision of noise vector, specify β _ipriori be that following gamma distributes:

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

κ ₀=10 ^-6, υ ₀=10 ^-6for form parameter and the scale parameter of gamma distribution;

(3) set up θ _iin the structural sparse priori of each element; Structural sparse priori is to be mainly reflected in the tree construction of wavelet coefficient, and when parent node is 0, its child node is also probably 0, and vice versa:

Wherein, s represents the yardstick s=1 in wavelet transformation ..., L, L is the total number of yardstick, θ _{s, ji}represent θ _ij the element at yardstick s place; δ ₀it is a distribution that concentrates on 0 place completely; π _{s, ji}be a scale-up factor, work as π _{s, ji}be tending towards at 0 o'clock, θ _{s, ji}0, otherwise, π _{s, ji}be tending towards at 1 o'clock, θ _{s, ji}come from a Gaussian distribution; π _{s, ji}tree construction by wavelet transformation is determined:

θ _{pa (s, ji)}θ _{s, ji}parent node, work as θ _{pa (s, ji)}=0 o'clock, π _{s, ji}be made as high probability and be 0

work as θ _{pa (s, ji)}≠ 0 o'clock, π _{s, ji}be made as high probability and be 1

β _s,iθ _{s, ji}the precision of Gaussian distribution in priori; π _r,

and β _{s, i}prior distribution as follows:

Wherein,

n _srepresent wavelet transformation s=1 ..., the number of the wavelet coefficient of L layer;

(4) sampling wavelet coefficient θ _{s, ji}; By above-mentioned (2), (3) step, obtained wavelet coefficient θ _{s, ji}likelihood function, θ _{s, ji}structural sparse priori, and the prior distribution of other unknown parameters; Need the wavelet coefficient Θ estimating, by following Bayesian inference, obtain:

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

Wherein,

for all noise precision, suppose the element independent same distribution in β and Ω, p (β) and p (Ω) are as follows respectively:

The Markov monte carlo method of employing based on gibbs sampler, needs the posteriority of predictor to distribute by being similar to from corresponding condition distribution sample drawn; θ _{s, ji}gibbs sampler as follows:

Wherein, other parameters that～expression condition needs in distributing, Φ _jrepresent Φ _jin j row;

(5) sampling noiset precision β _i; Sample distribution meets following gamma and distributes:

(6) the super parameter of sampling

concrete sample distribution is as follows:

Wherein, 1 (x) is 1 at x when being genuine, otherwise is 0;

(7) for each wavelet coefficient circulation in Θ successively execution (4), (5), (6) three steps, sample, until meet the condition of convergence or end condition, obtain the final wavelet coefficient of estimating m _sfor total number of sampling,

for the result of sampling each time; The end condition using is iterations, wherein, and sampling thief stabilization process 400 times, sampling process 200 times;

Step 4, according to linear mixed model (1) obtain rebuild high-spectral data