CN104463223A - Hyperspectral image group sparse demixing method based on empty spectral information abundance restraint - Google Patents

Abstract

The invention belongs to the technical field of image processing, and particularly discloses a hyperspectral image group sparse demixing method based on empty spectral information abundance restraint. The method comprises the steps of inputting a hyperspectral image data set and a standard spectrum bank, conducting adaptive grouping on hyperspectral image data with the mean-shift algorithm, conducting group sparse demixing on each group of hyperspectral image data, trimming the spectrum bank with the abundance matrix of each group of hyperspectral image data as the criteria, and outputting a hyperspectral image data sparse demixing result. According to the method, the structural features of hyperspectral data and spectrum bank data are considered, group sparse demixing and spectrum bank trimming are adopted for hyperspectral image demixing, and hyperspectral data sparse demixing precision is improved. The semi-supervision-based hyperspectral image demixing method has the advantages of sparse demixing and is high in demixing precision and effective.

Description

Based on the high spectrum image group sparse solution mixing method of empty spectrum information abundance constraint

Technical field

The invention belongs to technical field of image processing, further relate to remote sensing images technical field, specifically a kind of high spectrum image group sparse solution mixing method based on empty spectrum information abundance constraint.

Background technology

High-spectrum remote sensing technology has very fast development in recent years, and its research is mainly devoted to find learn and identify the true atural object class of high spectrum image technical method with making computer intelligence.High spectrum image has huge application prospect in all many-sides such as city planning, environment measuring, vegetative breakdown, military target detection and the identifications of mineral geology.Due to the impact of hyperspectral imager resolution and atural object terrain complexity, single pixel in image is made to there is the mixing of many kinds of substance, thus formation mixed pixel, have impact on the further decipher to high spectrum image and application, therefore, the high spectrum image solution technology of mixing becomes the subject topic that current remote sensing fields has Research Significance most.It is linear mixed model that high spectrum image solution the most general mixes model, this model mixed pixel described in high spectrum image is the linear combination of one group of a small amount of atural object and end member, and by noise, because this model is simple and practical, explicit physical meaning, often being solved EO-1 hyperion solution as basic model mixes problem.The process that EO-1 hyperion solution is mixed is: first, extracts the endmember spectra information existed in high spectrum image; Secondly, percent information and the abundance messages of different end member in mixed pixel is obtained.In addition, because high spectrum image has, data volume is large, redundant information is many, containing unfavorable factors such as noises, therefore require that mixing technical method in EO-1 hyperion solution has certain anti-interference and high-precision ability.

At present, EO-1 hyperion solution mixing method mainly comprises the three major types method based on geometry, statistics and sparse regression.More and more obtain accreditation and the research of a lot of people as a kind of semi-supervised method based on the EO-1 hyperion solution mixing method of sparse regression, it avoid the estimation of end member number and the extraction of end member spectrum information, utilize library of spectra can realize obtaining end member information and abundance estimation simultaneously.The end member number comprised due to pixel each in high spectrum image is far smaller than the dimension in library of spectra, the coefficient of the rarefaction representation obtained is sparse, meet the characteristic that Sparse represents, therefore, utilize the EO-1 hyperion solution mixing method based on sparse regression can obtain effective result.

The people such as Marian-Daniel Iordache are at paper " Sparse Unmixing of HyperspectralData " (GRS, 2011) propose basic sparse solution in mix, it separates mixed to the carrying out independently of each mixed pixel in high spectrum image, obtain end member information and abundance messages simultaneously, solve sparse solution and mix the method (SUnSAL) that model make use of variable division and augmentation Lagrange, this method is the theory based on alternating direction multiplier method.

Subsequently on this basis, the realization proposing a kind of improvement applies the solution mixing method of sparse constraint simultaneously to all pixels of high spectrum image, because all end member numbers comprised in high spectrum image are much smaller than library of spectra dimension, the abundance matrix obtained is that row is sparse, that is all pixels that a small amount of spectrum signature is used in matching high spectrum image are only had in library of spectra, by solving the problem of collaborative sparse regression, the solution improving high spectrum image mixes effect, what adopt is collaborative variable division and the algorithm (CLSUnSAL) of augmentation Lagrange, it is also the specific implementation of alternating direction multiplier method theory.But the weak point that above-mentioned two kinds of sparse solution mixing methods still exist is: do not consider the regional structure information of hyperspectral image data and weaken the impact that in library of spectra, between object spectrum, high correlation is mixed on sparse solution further.

Summary of the invention

The object of the invention is to overcome above-mentioned the deficiencies in the prior art, a kind of high spectrum image group sparse solution mixing method based on empty spectrum information abundance constraint is proposed, utilize the regional structure information of high spectrum image, weaken the impact of object spectrum height coherency in library of spectra, make high spectrum image sparse solution mix effect and be further enhanced.

Technical scheme of the present invention is:

(1) hyperspectral image data is inputted with standard spectrum storehouse wherein for the spectral signature of the pixel of i-th in hyperspectral image data, for the object spectrum of the jth in standard spectrum storehouse, L is wave band number, and n is pixel number, and m is the object spectrum number comprised in library of spectra, represent real number field;

(2) utilize mean shift algorithm to carry out adaptivenon-uniform sampling to hyperspectral image data Y, be divided into k region, each region is one group, then Y={G ₁... G _i..., G _k, wherein i=1 ..., k is that hyperspectral image data splits the i-th group of pixel collection obtained, | G _i| be the pixel number comprised in i-th group of pixel collection;

(3) by organizing the mixed method of sparse solution to the often group pixel collection G obtained in step (2) _i, i=1 ..., k carries out sparse solution respectively to be mixed:

3a) input splits by hyperspectral image data the one group of pixel collection G obtained _iwith library of spectra A, primary iteration number of times iter=0 is set, object spectrum number t=m in library of spectra A;

The model solution pixel collection G of 3b) utilization group sparse regression _icorresponding abundance matrix X _i, mathematical model is s.t X _i>=0, wherein Section 1 error term, represent the F norm of Arbitrary Matrix H, trace{} is the Section 2 of matrix trace, model the l of matrix _2,1norm, it constrains row matrix sparse, x ^jrepresent abundance matrix X _ijth row, m is the line number of abundance matrix, and λ is regular terms parameter;

3c) make iterations iter=iter+1, utilize step 3b) obtain pixel collection G _icorresponding abundance matrix X _iprune library of spectra A, upgrade library of spectra A, the concrete steps that library of spectra is pruned are:

1st step, input abundance matrix X _i=[x ¹... x ^j..., x ^m] ^twith library of spectra A=[a ₁... a _j..., a _m], x ^jfor the jth row vector of abundance matrix, a _jfor the jth row spectral signature in library of spectra;

2nd step, judges row vector x ^jin all elements whether be all greater than threshold tau, initial threshold τ size generally gets 2 × e ^-3if satisfy condition, then retain a corresponding in library of spectra _j, otherwise, from library of spectra A, reject a _j;

3rd step, obtains new library of spectra A, upgrade threshold tau=iter* τ and parametric t=| A|, | A| is the object spectrum number comprised in new library of spectra A;

3d) repeat step 3b) – 3c), until meet end condition, end condition is remain object spectrum number t≤T in maximum iteration time iter=20 or library of spectra, the span of threshold value T is p < T < 2p, p is the end member number comprised in hyperspectral image data, can be estimated by the EO-1 hyperion signal Recognition Algorithm based on least error to obtain;

(4) k group pixel collection { G is obtained ₁... G _i..., G _kcorresponding abundance matrix { X ₁... X _i..., X _k, wherein X _ibe i-th group of pixel collection G _icorresponding abundance matrix, exports hyperspectral image data solution and mixes result.

Above-mentioned steps 3b) in the concrete solution procedure of mathematical model be:

1st step, for mathematical model s.t X _ivariable G known in>=0 _i, A, λ do pre-service, order g _i=G _i/ const, A=A/const, λ=λ/const, const represents an intermediate variable, | G _i| represent pixel collection G _ithe pixel number inside comprised;

2nd step, introduces auxiliary variable U, makes U=X _i, then mathematical model is equivalent to following form:

$\underset{U,V_1,V_2,V_3}{\min }\frac{1}{2}{\left|\right|V_1-G\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|V_2\left|\right|}_{\mathrm{2,1}}+{\mathrm{\&iota;}}_{R+}\left(V_3\right)$

s.t V ₁＝AU (1)

V ₂＝U

V ₃＝U

Wherein for indicator function, | V ₃| be matrix V ₃columns, V _3ifor matrix V ₃i-th row, ι _r+(V _3i) mathematic(al) representation as follows:

${\mathrm{\&iota;}}_{R+}\left(V_{3i}\right)=\left\{\begin{array}{cc}0& V_{3i}\&GreaterEqual;0\\ \&infin;& V_{3i}<0\end{array}\right.---\left(2\right)$

3rd step, according to alternating direction multiplier method, introduces augmentation Lagrange multiplier D ₁/ μ, D ₂/ μ, D ₃/ μ, μ are constant, and formula (1) equivalence is converted to following form:

$\begin{array}{c}\mathrm{\&Gamma;}(U,V_1,V_2,V_3,D_1,D_2,D_3)\\ =\frac{1}{2}{\left|\right|V_1-G\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|V_2\left|\right|}_{\mathrm{2,1}}+{\mathrm{\&iota;}}_{R+}\left(V_3\right)+\\ \frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{AU}-V_1-D_1\left|\right|}_F^2+\frac{\mathrm{\&mu;}}{2}{\left|\right|U-V_2-D_2\left|\right|}_F^2+\\ \frac{\mathrm{\&mu;}}{2}{\left|\right|U-V_3-D_3\left|\right|}_F^2\end{array}---\left(3\right)$

4th step, if iterations η=0, constant μ>=0, initialization U ⁽⁰⁾, order initial value is zero, fixes its dependent variable, respectively changes persuing amount U, V ₁, V ₂, V ₃, D ₁, D ₂, D ₃value;

1) fixed variable V ₁, V ₂, V ₃, D ₁, D ₂, D ₃, the value of changes persuing amount U

$U^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_U\mathrm{\&Gamma;}(U,V_1^{\left(\mathrm{\&eta;}\right)},V_2^{\left(\mathrm{\&eta;}\right)},V_3^{\left(\mathrm{\&eta;}\right)},D_1^{\left(\mathrm{\&eta;}\right)},D_2^{\left(\mathrm{\&eta;}\right)},D_3^{\left(\mathrm{\&eta;}\right)})$

: U ^(η+1)← (A ^ta+2I) ^-1(A ^tξ ₁+ ξ ₂+ ξ ₃)

Wherein: ${\mathrm{\&xi;}}_1=V_1^{\left(\mathrm{\&eta;}\right)}+D_1^{\left(\mathrm{\&eta;}\right)},{\mathrm{\&xi;}}_2=V_2^{\left(\mathrm{\&eta;}\right)}+D_2^{\left(\mathrm{\&eta;}\right)},{\mathrm{\&xi;}}_3=V_3^{\left(\mathrm{\&eta;}\right)}+D_3^{\left(\mathrm{\&eta;}\right)}$

2) fixed variable U, V ₂, V ₃, D ₁, D ₂, D ₃, ask variable V ₁value

$V_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_{V_1}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1,V_2^{\left(\mathrm{\&eta;}\right)},V_3^{\left(\mathrm{\&eta;}\right)})$

: $V_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;\frac{1}{1+\mathrm{\&mu;}}[G+\mathrm{\&mu;}({\mathrm{AU}}^{\left(\mathrm{\&eta;}\right)}-D_1^{\left(\mathrm{\&eta;}\right)})]$

3) fixed variable U, V ₁, V ₃, D ₁, D ₂, D ₃, ask variable V ₂value

$V_2^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_{V_2}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1^{\left(\mathrm{\&eta;}\right)},V_2,V_3^{\left(\mathrm{\&eta;}\right)})$

Solve herein and need distinguish compute matrix often row wherein, for the r of matrix is capable, function vect_soft () is row vector soft-threshold function, and computing formula is:

4) fixed variable U, V ₁, V ₂, D ₁, D ₂, D ₃, ask variable V ₃value

$V_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_{V_3}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1^{\left(\mathrm{\&eta;}\right)},V_2^{\left(\mathrm{\&eta;}\right)},V_3)$

: $V_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;\max (U^{\left(\mathrm{\&eta;}\right)}-D_3^{\left(\mathrm{\&eta;}\right)},0)$

5) variables D is upgraded ₁, D ₂, D ₃

$D_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_1^{\left(\mathrm{\&eta;}\right)}-{\mathrm{AU}}^{(\mathrm{\&eta;}+1)}+V_1^{(\mathrm{\&eta;}+1)}$

$D_2^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_2^{\left(\mathrm{\&eta;}\right)}-U^{(\mathrm{\&eta;}+1)}+V_2^{(\mathrm{\&eta;}+1)}$

$D_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_3^{\left(\mathrm{\&eta;}\right)}-U^{(\mathrm{\&eta;}+1)}+V_3^{(\mathrm{\&eta;}+1)}$

5th step, makes iterations η=η+1, repeats above-mentioned to variable U, V ₁, V ₂, V ₃, D ₁, D ₂, D ₃solution procedure, until meet end condition, end condition is the error of before and after maximum iteration time η=200 or variable U twice the general value size of threshold epsilon is 1 × e ^-5, the value of output variable U, i.e. abundance matrix X _i.

Beneficial effect of the present invention: the present invention can by the adaptive grouping to high spectrum image, reasonably utilize the regional structure information of high spectrum image, on this basis, the advantage that utilization group sparse solution mixes and the advantage that library of spectra is pruned, realize being separated mixed pixel in esse in high spectrum image, reach the further decipher to image and application.Compared with prior art, the present invention has the following advantages:

First, the present invention mixes field mean shift segmentation algorithm application in high spectrum image solution, it is mainly as the preprocessing means that EO-1 hyperion solution is mixed, the combination that existing high spectrum image sparse solution mixes empty spectrum information will contribute to improving the mixed effect of solution, utilize the sky spectrum prior imformation of high spectrum image, think that the EO-1 hyperion pixel in neighborhood has identical end member and close Abundances, and actual high-spectrum similarly is be made up of different regions, high spectrum image is carried out the group sparse solution that adaptivenon-uniform sampling contributes to carrying out to mix, make the present invention's closing to reality application more.

Second, utilization group sparse solution mixing method of the present invention carries out high spectrum image solution to be mixed, using to the zones of different of high spectrum image adaptivenon-uniform sampling as different pixels group, in utilization group, EO-1 hyperion pixel has the information of identical constituent and close ratio, the further advantage having played group sparse solution and mixed, improve the precision that high spectrum image solution is mixed, the present invention is had and better separates mixed effect.

3rd, present invention utilizes the library of spectra pruning method of abundance constraint, reject partial spectrum feature in library of spectra, reduce the impact that library of spectra spectrum signature high correlation is mixed on sparse solution, the solution making the present invention further increase high spectrum image mixes precision.

Below with reference to accompanying drawing, the present invention is described in further details.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the inventive method;

Fig. 2 is high spectrum image in emulation experiment of the present invention;

Fig. 3 is one piece of region in emulation experiment of the present invention in high spectrum image;

Fig. 4 is the contrast that the inventive method True Data solution mixes result and additive method;

Fig. 4 (a) is algorithm NCLS abundance drawing for estimate;

Fig. 4 (b) is algorithm SUnSAL abundance drawing for estimate;

Fig. 4 (c) is algorithm CLSUnSAL abundance drawing for estimate;

Fig. 4 (d) is the inventive method abundance drawing for estimate.

Concrete implementing measure

Embodiment 1:

1 concrete steps of the present invention are described below by reference to the accompanying drawings:

(1) hyperspectral image data is inputted with standard spectrum storehouse wherein for the spectral signature of the pixel of i-th in hyperspectral image data, for the object spectrum of the jth in standard spectrum storehouse, L is wave band number, and n is pixel number, and m is the object spectrum number comprised in library of spectra, represent real number field;

Simulated data sets size is 224 × 30 × 30, arrange onesize region by nine piece of three row three to form, each area size is 10 × 10, the end member kind comprised in each region is different with number, simulated data concentrates the end member number comprised to be 9, end member is random selecting from library of spectra, and abundance obeys Dirichlet distribute.

True Data collection: Nevada, USA area mineral data set (as shown in Figure 2) obtained by airborne visible ray/Infrared Imaging Spectrometer AVIRIS, wherein in actual emulation process, we utilize wherein one piece of region, size is 250 × 191, and each pixel has 188 spectrums wave band (as shown in Figure 3).

Input standard spectrum database data A, stem from the splib06 library of spectra that US Geological Survey USGS issued in 2007, in library of spectra, the spectral band number of all object spectrums is 224.

(2) utilize mean shift algorithm to carry out adaptivenon-uniform sampling to hyperspectral image data Y, be divided into k region, each region is one group, then Y={G ₁... G _i..., G _k, wherein i=1 ..., k is that hyperspectral image data splits the i-th group of pixel collection obtained, | G _i| be the pixel number comprised in i-th group of pixel collection;

(3) by organizing the mixed method of sparse solution to the often group pixel collection G obtained in step (2) _i, i=1 ..., k carries out sparse solution respectively to be mixed:

3a) input splits by hyperspectral image data the one group of pixel collection G obtained _iwith library of spectra A, primary iteration number of times iter=0 is set, object spectrum number t=m in library of spectra A;

The model solution pixel collection G of 3b) utilization group sparse regression _icorresponding abundance matrix X _i, mathematical model is s.t X _i>=0, wherein Section 1 error term, represent the F norm of Arbitrary Matrix H, trace{} is the Section 2 of matrix trace, model the l of matrix _2,1norm, it constrains row matrix sparse, makes non-zero row as far as possible few, x ^jrepresent abundance matrix X _ijth row, m is the line number of abundance matrix, and λ is regular terms parameter, and the mode generally manually regulated makes final solution mix result to reach necessary requirement;

3c) make iterations iter=iter+1, utilize step 3b) obtain pixel collection G _icorresponding abundance matrix X _iprune library of spectra A, upgrade library of spectra A, the concrete steps that library of spectra is pruned are:

1st step, input abundance matrix X _i=[x ¹... x ^j..., x ^m] ^twith library of spectra A=[a ₁... a _j..., a _m], x ^jfor the jth row vector of abundance matrix, a _jfor the jth row spectral signature in library of spectra;

2nd step, judges row vector x ^jin all elements whether be all greater than threshold tau, initial threshold τ size generally gets 2 × e ^-3if satisfy condition, then retain a corresponding in library of spectra _j, otherwise, from library of spectra A, reject a _j;

3rd step, obtains new library of spectra A, upgrade threshold tau=iter* τ and parametric t=| A|, | A| is the object spectrum number comprised in new library of spectra A;

3d) repeat step 3b) – 3c), until meet end condition, end condition is remain object spectrum number t≤T in maximum iteration time iter=20 or library of spectra, the span of threshold value T is p < T < 2p, p is the end member number comprised in hyperspectral image data, can be estimated by the EO-1 hyperion signal Recognition Algorithm based on least error to obtain;

(4) k group pixel collection { G is obtained ₁... G _i..., G _kcorresponding abundance matrix { X ₁... X _i..., X _k, wherein X _ibe i-th group of pixel collection G _icorresponding abundance matrix, exports hyperspectral image data solution and mixes result.

According to the above-mentioned high spectrum image group sparse solution mixing method based on the constraint of empty spectrum information abundance, step 3b) in the concrete solution procedure of mathematical model be:

1st step, for mathematical model s.t X _ivariable G known in>=0 _i, A, λ do pre-service, order g _i=G _i/ const, A=A/const, λ=λ/const, const represents an intermediate variable, | G _i| represent pixel collection G _ithe pixel number inside comprised;

2nd step, introduces auxiliary variable U, makes U=X _i, then mathematical model is equivalent to following form:

$\underset{U,V_1,V_2,V_3}{\min }\frac{1}{2}{\left|\right|V_1-G\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|V_2\left|\right|}_{\mathrm{2,1}}+{\mathrm{\&iota;}}_{R+}\left(V_3\right)$

s.t V ₁＝AU (1)

V ₂＝U

V ₃＝U

Wherein for indicator function, | V ₃| be matrix V ₃columns, V _3ifor matrix V ₃i-th row, ι _r+(V _3i) mathematic(al) representation as follows:

${\mathrm{\&iota;}}_{R+}\left(V_{3i}\right)=\left\{\begin{array}{cc}0& V_{3i}\&GreaterEqual;0\\ \&infin;& V_{3i}<0\end{array}\right.---\left(2\right)$

3rd step, according to alternating direction multiplier method, introduces augmentation Lagrange multiplier D ₁/ μ, D ₂/ μ, D ₃/ μ, μ are constant, and formula (1) equivalence is converted to following form:

$\begin{array}{c}\mathrm{\&Gamma;}(U,V_1,V_2,V_3,D_1,D_2,D_3)\\ =\frac{1}{2}{\left|\right|V_1-G\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|V_2\left|\right|}_{\mathrm{2,1}}+{\mathrm{\&iota;}}_{R+}\left(V_3\right)+\\ \frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{AU}-V_1-D_1\left|\right|}_F^2+\frac{\mathrm{\&mu;}}{2}{\left|\right|U-V_2-D_2\left|\right|}_F^2+\\ \frac{\mathrm{\&mu;}}{2}{\left|\right|U-V_3-D_3\left|\right|}_F^2\end{array}---\left(3\right)$

4th step, if iterations η=0, constant μ>=0, initialization U ⁽⁰⁾, order initial value is zero, fixes its dependent variable, respectively changes persuing amount U, V ₁, V ₂, V ₃, D ₁, D ₂, D ₃value;

1) fixed variable V ₁, V ₂, V ₃, D ₁, D ₂, D ₃, the value of changes persuing amount U

$U^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_U\mathrm{\&Gamma;}(U,V_1^{\left(\mathrm{\&eta;}\right)},V_2^{\left(\mathrm{\&eta;}\right)},V_3^{\left(\mathrm{\&eta;}\right)},D_1^{\left(\mathrm{\&eta;}\right)},D_2^{\left(\mathrm{\&eta;}\right)},D_3^{\left(\mathrm{\&eta;}\right)})$

: U ^(η+1)← (A ^ta+2I) ^-1(A ^tξ ₁+ ξ ₂+ ξ ₃)

Wherein: ${\mathrm{\&xi;}}_1=V_1^{\left(\mathrm{\&eta;}\right)}+D_1^{\left(\mathrm{\&eta;}\right)},{\mathrm{\&xi;}}_2=V_2^{\left(\mathrm{\&eta;}\right)}+D_2^{\left(\mathrm{\&eta;}\right)},{\mathrm{\&xi;}}_3=V_3^{\left(\mathrm{\&eta;}\right)}+D_3^{\left(\mathrm{\&eta;}\right)}$

2) fixed variable U, V ₂, V ₃, D ₁, D ₂, D ₃, ask variable V ₁value

$V_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_{V_1}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1,V_2^{\left(\mathrm{\&eta;}\right)},V_3^{\left(\mathrm{\&eta;}\right)})$

: $V_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;\frac{1}{1+\mathrm{\&mu;}}[G+\mathrm{\&mu;}({\mathrm{AU}}^{\left(\mathrm{\&eta;}\right)}-D_1^{\left(\mathrm{\&eta;}\right)})]$

3) fixed variable U, V ₁, V ₃, D ₁, D ₂, D ₃, ask variable V ₂value

$V_2^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_{V_2}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1^{\left(\mathrm{\&eta;}\right)},V_2,V_3^{\left(\mathrm{\&eta;}\right)})$

Solve herein and need distinguish compute matrix often row wherein, for the r of matrix is capable, function vectsoft () is row vector soft-threshold function, and computing formula is:

4) fixed variable U, V ₁, V ₂, D ₁, D ₂, D ₃, ask variable V ₃value

$V_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;{\arg \min }_{V_3}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1^{\left(\mathrm{\&eta;}\right)},V_2^{\left(\mathrm{\&eta;}\right)},V_3)$

: $V_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;\max (U^{\left(\mathrm{\&eta;}\right)}-D_3^{\left(\mathrm{\&eta;}\right)},0)$

5) variables D is upgraded ₁, D ₂, D ₃

$D_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_1^{\left(\mathrm{\&eta;}\right)}-{\mathrm{AU}}^{(\mathrm{\&eta;}+1)}+V_1^{(\mathrm{\&eta;}+1)}$

$D_2^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_2^{\left(\mathrm{\&eta;}\right)}-U^{(\mathrm{\&eta;}+1)}+V_2^{(\mathrm{\&eta;}+1)}$

$D_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_3^{\left(\mathrm{\&eta;}\right)}-U^{(\mathrm{\&eta;}+1)}+V_3^{(\mathrm{\&eta;}+1)}$

5th step, makes iterations η=η+1, repeats above-mentioned to variable U, V ₁, V ₂, V ₃, D ₁, D ₂, D ₃solution procedure, until meet end condition, end condition is the error of before and after maximum iteration time η=200 or variable U twice the general value size of threshold epsilon is 1 × e ^-5, the value of output variable U, i.e. abundance matrix X _i.

Embodiment 2:

Below in conjunction with accompanying drawing 3 and accompanying drawing 4, effect of the present invention is described further.

Emulation experiment of the present invention is at Intel Core (TM) 2Duo CPU, dominant frequency 2.00GHz, and the MATLAB R2011b of internal memory 2G, Windows 7 on platform realizes.

Emulation of the present invention is the experiment simulation done on simulated data sets and True Data collection, simulated data is made up of nine fritters of 10 × 10, the end member number that every block comprises is different, end member is random selecting from nine spectrum signatures, the abundance of all fritters obeys Dirichlet distribute, simulated data sets, and size is 224 × 900, data disturb by different stage white Gaussian noise, signal to noise ratio snr (dB)=Ε || Ax|| ²/ Ε || n|| ²be respectively: 20dB, 30dB and 40dB.

True Data is the mining area data set (as shown in Figure 2) in Nevada, USA area, wherein in actual emulation process, we utilize wherein one piece of region, size is 250 × 191, to remove in wave band some by the wave band of water vapor and noise severe jamming, residue has 188 spectrums wave band (as shown in Figure 3)

It is signal and reconstruction error signal ratio (SRE) that simulated data collected explanations or commentaries mixes precision evaluation index, and mathematic(al) representation is as follows:

$\mathrm{SRE}\left(\mathrm{dB}\right)=20\log \left(E\right[{\left|\right|X\left|\right|}_F^2]/E[{\left|\right|X-\stackrel{\&OverBar;}{X}\left|\right|}_F^2\left]\right)$

Wherein: X is true abundance matrix, for estimating abundance matrix.SRE value is larger, represents that the mixed effect of solution is better.

The EO-1 hyperion solution that the present invention does on simulated data sets is mixed the result that effect and nonnegativity restrictions least square (NCLS), SUnSAL, CLSUnSAL obtain and is contrasted, and it is as shown in table 1 that the EO-1 hyperion solution obtained mixes Contrast on effect.

The solution of table 1 distinct methods on simulated data sets mixes precision

SNR	NCLS	SUnSAL	CLSUnSAL	The inventive method
					20dB	4.3704dB	5.5499dB	6.2458dB	7.0340dB
30dB	9.6123dB	11.3869dB	12.5738dB	14.9755dB
					40dB	16.0141dB	19.3443dB	21.8499dB	27.3691dB

As can be seen from Table 1, the simulation result on simulated data sets, under different noise, to mix ratio of precision existing sparse solution mixing method high for solution of the present invention.

The present invention on the mining area data set in Nevada, USA area and algorithms of different NCLS, SUnSAL, CLSUnSAL emulation experiment, result (as shown in Figure 4) can find out that the present invention can obtain higher solution and mix effect by experiment, mineral (mineral matter is followed successively by alunit-alunite, buddingtonite-water ammonium feldspar, chalcedony-calcedony from left to right) the abundance figure obtained is more clear, and often kind of mineral matter distribution is more concentrated.

In sum, the present invention can by the adaptive grouping to high spectrum image, reasonably utilize the regional structure information of high spectrum image, on this basis, the advantage that utilization group sparse solution mixes and the advantage that library of spectra is pruned, realize being separated mixed pixel in esse in high spectrum image, reach the further decipher to image and application.Compared with prior art, the present invention has the following advantages:

First, the present invention mixes field mean shift segmentation algorithm application in high spectrum image solution, it is mainly as the preprocessing means that EO-1 hyperion solution is mixed, the combination that existing high spectrum image sparse solution mixes empty spectrum information will contribute to improving the mixed effect of solution, utilize the sky spectrum prior imformation of high spectrum image, think that the EO-1 hyperion pixel in neighborhood has identical end member and close Abundances, and actual high-spectrum similarly is be made up of different regions, high spectrum image is carried out the group sparse solution that adaptivenon-uniform sampling contributes to carrying out to mix, make the present invention's closing to reality application more.

Second, utilization group sparse solution mixing method of the present invention carries out high spectrum image solution to be mixed, using to the zones of different of high spectrum image adaptivenon-uniform sampling as different pixels group, in utilization group, EO-1 hyperion pixel has the information of identical constituent and close ratio, the further advantage having played group sparse solution and mixed, improve the precision that high spectrum image solution is mixed, the present invention is had and better separates mixed effect.

3rd, present invention utilizes the library of spectra pruning method of abundance constraint, reject partial spectrum feature in library of spectra, reduce the impact that library of spectra spectrum signature high correlation is mixed on sparse solution, the solution making the present invention further increase high spectrum image mixes precision.

Therefore, the present invention is one more effective high spectrum image sparse solution mixing method.

The part do not described in detail in present embodiment belongs to the known conventional means of the industry, does not describe one by one here.More than exemplifying is only illustrate of the present invention, does not form the restriction to protection scope of the present invention, everyly all belongs within protection scope of the present invention with the same or analogous design of the present invention.

Claims (2)

1., based on the high spectrum image group sparse solution mixing method of empty spectrum information abundance constraint, it is characterized in that: comprise the following steps:

(1) hyperspectral image data is inputted with standard spectrum storehouse wherein for the spectral signature of the pixel of i-th in hyperspectral image data, for the object spectrum of the jth in standard spectrum storehouse, L is wave band number, and n is pixel number, and m is the object spectrum number comprised in library of spectra, represent real number field;

(2) utilize mean shift algorithm to carry out adaptivenon-uniform sampling to hyperspectral image data Y, be divided into k region, each region is one group, then Y={G ₁... G _i..., G _k, wherein i=1 ..., k is that hyperspectral image data splits the i-th group of pixel collection obtained, | G _i| be the pixel number comprised in i-th group of pixel collection;

(3) by organizing the mixed method of sparse solution to the often group pixel collection G obtained in step (2) _i, i=1 ..., k carries out sparse solution respectively to be mixed:

3a) input splits by hyperspectral image data the one group of pixel collection G obtained _iwith library of spectra A, primary iteration number of times iter=0 is set, object spectrum number t=m in library of spectra A;

The model solution pixel collection G of 3b) utilization group sparse regression _icorresponding abundance matrix X _i, mathematical model is $\underset{X_i}{\min }1/2{\left|\right|{\mathrm{AX}}_i-G_i\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|X_i\left|\right|}_{\mathrm{2,1}}s.t{X}_i\&GreaterEqual;0,$ Wherein Section 1 error term, represent the F norm of Arbitrary Matrix H, trace{} is the Section 2 of matrix trace, model the l of matrix _2,1norm, it constrains row matrix sparse, x ^jrepresent abundance matrix X _ijth row, m is the line number of abundance matrix, and λ is regular terms parameter;

3c) make iterations iter=iter+1, utilize step 3b) obtain pixel collection G _icorresponding abundance matrix X _iprune library of spectra A, upgrade library of spectra A, the concrete steps that library of spectra is pruned are:

1st step, input abundance matrix X _i=[x ¹... x ^j..., x ^m] ^twith library of spectra A=[a ₁... a _j..., a _m], x ^jfor the jth row vector of abundance matrix, a _jfor the jth row spectral signature in library of spectra;

2nd step, judges row vector x ^jin all elements whether be all greater than threshold tau, initial threshold τ size generally gets 2 × e ^-3if satisfy condition, then retain a corresponding in library of spectra _j, otherwise, from library of spectra A, reject a _j;

3rd step, obtains new library of spectra A, upgrade threshold tau=iter* τ and parametric t=| A|, | A| is the object spectrum number comprised in new library of spectra A;

3d) repeat step 3b) – 3c), until meet end condition, end condition is remain object spectrum number t≤T in maximum iteration time iter=20 or library of spectra, the span of threshold value T is p < T < 2p, p is the end member number comprised in hyperspectral image data, can be estimated by the EO-1 hyperion signal Recognition Algorithm based on least error to obtain;

(4) k group pixel collection { G is obtained ₁... G _i..., G _kcorresponding abundance matrix { X ₁... X _i..., X _k, wherein X _ibe i-th group of pixel collection G _icorresponding abundance matrix, exports hyperspectral image data solution and mixes result.

2. the high spectrum image group sparse solution mixing method based on the constraint of empty spectrum information abundance according to claim 1, described step 3b) in the concrete solution procedure of mathematical model be:

1st step, for mathematical model $\underset{X_i}{\min }1/2{\left|\right|{\mathrm{AX}}_i-G_i\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|X_i\left|\right|}_{\mathrm{2,1}}s.t{X}_i\&GreaterEqual;0$ In known variable G _i, A, λ do pre-service, order g _i=G _i/ const, A=A/const, λ=λ/const, const represents an intermediate variable, | G _i| represent pixel collection G _ithe pixel number inside comprised;

2nd step, introduces auxiliary variable U, makes U=X _i, then mathematical model is equivalent to following form:

$\underset{U,V_1,V_2,V_3}{\min }\frac{1}{2}{\left|\right|V_1-G\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|V_2\left|\right|}_{\mathrm{2,1}}+{\mathrm{\&iota;}}_{R+}\left(V_3\right)$

s.t V ₁＝AU (1)

V ₂＝U

V ₃＝U

Wherein for indicator function, | V ₃| be matrix V ₃columns, V _3ifor matrix V ₃i-th row, ι _r+(V _3i) mathematic(al) representation as follows:

${\mathrm{\&iota;}}_{R+}\left(V_{3i}\right)=\left\{\begin{array}{cc}0& V_{3i}\&GreaterEqual;0\\ \&infin;& V_{3i}<0\end{array}\right.---\left(2\right)$

3rd step, according to alternating direction multiplier method, introduces augmentation Lagrange multiplier D ₁/ μ, D ₂/ μ, D ₃/ μ, μ are constant, and formula (1) equivalence is converted to following form:

$\begin{array}{c}\mathrm{\&Gamma;}(U,V_1,V_2,V_3,D_1,D_2,D_3)\\ =\frac{1}{2}{\left|\right|V_1-G\left|\right|}_F^2+\mathrm{\&lambda;}{\left|\right|V_2\left|\right|}_{\mathrm{2,1}}+{\mathrm{\&iota;}}_{R+}\left(V_3\right)+\\ \frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{AU}-V_1-D_1\left|\right|}_F^2+\frac{\mathrm{\&mu;}}{2}{\left|\right|U-V_2-D_2\left|\right|}_F^2+\\ \frac{\mathrm{\&mu;}}{2}{\left|\right|U-V_3-D_3\left|\right|}_F^2\end{array}---\left(3\right)$

4th step, if iterations η=0, constant μ>=0, initialization U ⁽⁰⁾, order initial value is zero, fixes its dependent variable, respectively changes persuing amount U, V ₁, V ₂, V ₃, D ₁, D ₂, D ₃value;

1) fixed variable V ₁, V ₂, V ₃, D ₁, D ₂, D ₃, the value of changes persuing amount U

$U^{(\mathrm{\&eta;}+1)}\&LeftArrow;\arg {\min }_U\mathrm{\&Gamma;}(U,V_1^{\left(\mathrm{\&eta;}\right)},V_2^{\left(\mathrm{\&eta;}\right)},V_3^{\left(\mathrm{\&eta;}\right)},D_1^{\left(\mathrm{\&eta;}\right)},D_2^{\left(\mathrm{\&eta;}\right)},D_3^{\left(\mathrm{\&eta;}\right)})$

: U ^(η+1)← (A ^ta+2I) ^-1(A ^tξ ₁+ ξ ₂+ ξ ₃)

Wherein: ${\mathrm{\&xi;}}_1=V_1^{\left(\mathrm{\&eta;}\right)}+D_1^{\left(\mathrm{\&eta;}\right)},{\mathrm{\&xi;}}_2=V_2^{\left(\mathrm{\&eta;}\right)}+D_2^{\left(\mathrm{\&eta;}\right)},{\mathrm{\&xi;}}_3=V_3^{\left(\mathrm{\&eta;}\right)}+D_3^{\left(\mathrm{\&eta;}\right)}$

2) fixed variable U, V ₂, V ₃, D ₁, D ₂, D ₃, ask variable V ₁value

$V_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;\arg {\min }_{V_1}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1,V_2^{\left(\mathrm{\&eta;}\right)},V_3^{\left(\mathrm{\&eta;}\right)})$

: $V_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;\frac{1}{1+\mathrm{\&mu;}}[G+\mathrm{\&mu;}({\mathrm{AU}}^{\left(\mathrm{\&eta;}\right)}-D_1^{\left(\mathrm{\&eta;}\right)})]$

3) fixed variable U, V ₁, V ₃, D ₁, D ₂, D ₃, ask variable V ₂value

$V_2^{(\mathrm{\&eta;}+1)}\&LeftArrow;\arg {\min }_{V_2}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1^{\left(\mathrm{\&eta;}\right)},V_2,V_3^{\left(\mathrm{\&eta;}\right)})$

Solve herein and need distinguish compute matrix often row wherein, for the r of matrix is capable, function vectsoft () is row vector soft-threshold function, and computing formula is:

4) fixed variable U, V ₁, V ₂, D ₁, D ₂, D ₃, ask variable V ₃value

$V_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;\arg {\min }_{V_3}\mathrm{\&Gamma;}(U^{\left(\mathrm{\&eta;}\right)},V_1^{\left(\mathrm{\&eta;}\right)},V_2^{\left(\mathrm{\&eta;}\right)},V_3)$

: $V_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;\max (U^{\left(\mathrm{\&eta;}\right)}-D_3^{\left(\mathrm{\&eta;}\right)},0)$

5) variables D is upgraded ₁, D ₂, D ₃

$D_1^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_1^{\left(\mathrm{\&eta;}\right)}-{\mathrm{AU}}^{(\mathrm{\&eta;}+1)}+V_1^{(\mathrm{\&eta;}+1)}$

$D_2^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_2^{\left(\mathrm{\&eta;}\right)}-U^{(\mathrm{\&eta;}+1)}+V_2^{(\mathrm{\&eta;}+1)}$

$D_3^{(\mathrm{\&eta;}+1)}\&LeftArrow;D_3^{\left(\mathrm{\&eta;}\right)}-U^{(\mathrm{\&eta;}+1)}+V_3^{(\mathrm{\&eta;}+1)}$

5th step, makes iterations η=η+1, repeats above-mentioned to variable U, V ₁, V ₂, V ₃, D ₁, D ₂, D ₃solution procedure, until meet end condition, end condition is the error of before and after maximum iteration time η=200 or variable U twice the general value size of threshold epsilon is 1 × e ^-5, the value of output variable U, i.e. abundance matrix X _i.