CN105320959A - End member learning based hyperspectral image sparse unmixing method - Google Patents

Abstract

The invention discloses an end member learning based hyperspectral image sparse unmixing method, which mainly solves the problems of low hyperspectral image unmixing precision, poor reconstruction effect, long consumed time and low efficiency during a low-signal-noise-ratio hyperspectral image unmixing process in the prior art. The method comprises the steps of: inputting hyperspectral data; synthesizing hyperspectral base data; performing end member learning; solving a hyperspectral data abundance matrix; calculating a reconstruction error of the hyperspectral data abundance matrix; and outputting an unmixing result. The method adopts a new solving mode, introduces an end member learning thought, has the advantages of high unmixing precision, good reconstruction effect and high efficiency, is simple in solving step and explicit in principle, and can be used for understanding interpretation of hyperspectral images.

Description

Based on the high spectrum image sparse solution mixing method of end member study

Technical field

The invention belongs to technical field of image processing, further relate to image solution mix in technical field a kind of based on end member study high spectrum image solution mixing method.The present invention, by first carrying out simulation learning to end member, then utilizes the end member learnt to carry out solving of abundance, thus accomplishes the fast understanding decipher to high spectrum image.The present invention can be used for separating mixed process to the high spectrum image of various digital device, effectively can improve the precision that high spectrum image solution is mixed.

Background technology

High-spectrum similarly is be made up of up to a hundred very narrow wave bands, and it not only has the information of spectral domain, further comprises abundant spatial information.But, the deficiency of sensor spatial resolution, make for high spectrum image, pixel is wherein difficult to be pure pixel, but the mixed pixel point merged by many kinds of substance, in order to can be better, effectively utilize hyperspectral image data, just mixed pixel point wherein must be decomposed, be broken down into material collection (being commonly called as end member) existing in image and the product of corresponding proportion set (being commonly called as abundance).

Y.Qian, S.Jia, paper " HyperspectralunmixingviaL0.5sparsity-constrainednonnegat the ivematrixfactorization " (IEEETransactionsonGeoscienceandRemoteSensing that J.ZhouandA.Robles-Kelly delivers at it, vol.49, no11, pp.4282-4297, No.2011). in propose a kind of non-negative matrix factorization method based on L0.5 norm sparse constraint.The method, by carrying out the sparse constraint of L0.5 norm to abundance matrix, uses alternative iteration method, carries out Non-negative Matrix Factorization, thus obtain end member matrix and abundance matrix to high-spectral data matrix.The deficiency that the method exists sneaks out in journey in the high spectrum image solution of low signal-to-noise ratio, and when using alternative iteration method to solve end member matrix and abundance matrix, high spectrum image solution mixes result length consuming time, efficiency is low simultaneously.

The patented technology " a kind of high spectrum image sparse solution mixing method based on accidental projection " (number of patent application: 201110207433.0, Authorization Notice No.: CN102314685A) that Beijing Space aviation university has proposes a kind of high spectrum image sparse solution mixing method based on accidental projection.The method has four large steps: one, by digital independent in software MATLAB; Two, computing machine carries out accidental projection to hyperspectral image data and EO-1 hyperion database data; Three, build the mixed objective function of sparse solution, use division Bregman algorithm optimization objective function to ask extreme value, until reach convergence stop condition.Four, set suitable threshold process abundance matrix, obtain final abundance figure and end member.Present invention utilizes high-spectral data storehouse to select end member, overcome the shortcoming that end member in the past calculated by algorithm and the pure substance spectra in standard high-spectral data storehouse cannot be tightly corresponding; Achieve the qualitative analysis to high spectrum image.The deficiency that the method still exists is, jljl different spectrum phenomenon can cause end member in actual end member and standard high-spectral data storehouse to there is certain difference, directly utilizes end member in high-spectral data storehouse, and high spectrum image solution can be made to mix result precision low, quality reconstruction is poor.

The content of invention

The object of the invention is to the deficiency for above-mentioned prior art, a kind of high spectrum image sparse solution mixing method based on end member study is proposed, mix precision with the solution improving high spectrum image, overcome high spectrum image solution and mix inefficient problem, reduce high spectrum image solution and mix consuming time.

To achieve these goals, step of the present invention comprises as follows:

1., based on a high spectrum image sparse solution mixing method for end member study, comprise the steps:

(1) high-spectral data is inputted;

(2) EO-1 hyperion base data is synthesized:

(2a) from digital spectrum storehouse, select all end members comprised in high-spectral data, obtain alternative end member;

(2b) with alternative end member initialization predeterminable area, alternative area is obtained;

(2c) utilize Di Li Cray method, generate the Abundances of alternative area;

(2d) be multiplied by end member corresponding in alternative area with the Abundances generated, obtain initial base data;

(2e) by low-pass filter, the high frequency signal components in the initial base data of filtering;

(2f) from initial base data, choose the pixel of Abundances lower than predetermined threshold value 0.8, obtain alternate pixel point;

(2g) with end member initialization alternate pixel points all in alternative end member, the Abundances that in the alternate pixel be initialised point, each end member is corresponding is set to the inverse of alternative end member sum, obtains intermediate base data;

(2h) in intermediate base data, add zero mean Gaussian white noise, obtain the EO-1 hyperion base data of synthesis;

(3) end member study:

(3a) to preset size be the screening matrix of L × P and size is the storage matrix of 1 × P, wherein, L represents the wave band number of end member in alternative end member, P to represent in alternative end member end member sum, and the element that to be the screening matrix of L × P and size by size be in the storage matrix of 1 × P is initialized as complete zero;

(3b) according to the following formula, slickness bound term is constructed:

G＝||A|| _F

Wherein, G represents slickness bound term, and A represents end member matrix to be learned, || || _frepresent the operation of getting F norm;

(3c) according to the following formula, structural segmentation slickness bound term:

$R=\underset{l=1,p=1,i=1}{\overset{L,P,2}{\Σ}}(-e^{\left(-{\left(A_{lp}-B_i\right)}^2/\mathrm{\γ}\right)}+1)$

Wherein, R represents sectionally smooth bound term, and ∑ represents gets sum operation, and l represents rower, and the span of l is { 1,2,, L}, L represent the wave band number of end member in alternative end member, and p represents row mark, and the span of p is { 1,2 ..., P}, P represent end member sum in alternative end member, and i represents A _lpthe label of element in the Neighbourhood set of left and right, e ⁽⁾expression take natural number as the index operation at the end, A _lprepresent the capable p column element of l in end member matrix to be learned, B _irepresent A _lpi-th element in the Neighbourhood set of left and right, the span of i is that { 1,2}, γ represent constraining force parameter, and the span of γ is [0,1];

(3d) utilize K averaging method, the EO-1 hyperion base data of synthesis is carried out cluster operation, obtain the EO-1 hyperion base data after cluster;

(3e) from the EO-1 hyperion base data after cluster, random selecting does not also carry out class data of end member study, obtains current subclass EO-1 hyperion base data;

(3f) according to the following formula, end member study is carried out to current subclass EO-1 hyperion base data, obtains the end member matrix of current subclass EO-1 hyperion base data:

$A^{(k+1)}=\mathrm{argmin}\frac{1}{2}\left|\right|Z-A^{\left(k\right)}{Y}^{\left(k\right)}|{|}_F^2+{\mathrm{\λ}}_1{G}^{\left(k\right)}+{\mathrm{\λ}}_2{R}^{\left(k\right)}+\frac{1}{2}{\mathrm{\λ}}_3||Y^{\left(k\right)}-D|{|}_F^2$

Wherein, A ^(k+1)represent the end member matrix of the current subclass EO-1 hyperion base data of kth+1 iteration, k represents iterations used when carrying out end member study to current subclass EO-1 hyperion base data, the span of k is { 1,2, ..., the initial value of 100}, k is set to 1, argmin represents the end member matrix manipulation got when carrying out end member study to current subclass EO-1 hyperion base data and reaching minimum value represent the square operation getting F norm, Z represents current subclass EO-1 hyperion base data, A ^(k)represent the end member matrix of the current subclass EO-1 hyperion base data that kth is secondary, Y ^(k)represent the abundance matrix of the current subclass EO-1 hyperion base data of kth time iteration, λ ₁represent the parameter regulating slickness bound term, λ ₁value be set to 0.9, G ^(k)represent the slickness bound term of kth time iteration, λ ₂represent the parameter regulating sectionally smooth bound term, λ ₂value be set to 0.1, R ^(k)represent the sectionally smooth bound term of kth time iteration, λ ₃represent balance parameters, λ ₃value be set to 1, D and represent the actual value that the abundance matrix of current subclass EO-1 hyperion base data is corresponding;

(3g) utilize digital spectrum storehouse, the end member matrix of the current subclass EO-1 hyperion base data after study is screened, obtains the end member matrix closest to digital spectrum storehouse;

(3h) judge that the whether every class of the EO-1 hyperion base data after cluster has all carried out end member study, if so, obtain the end member matrix after learning and current storage matrix, otherwise, perform step (3e);

(4) high-spectral data abundance matrix is solved:

(4a) the canonical Weighted Constraint item of high-spectral data abundance matrix according to the following formula, is constructed:

$Q=\underset{i=1,j=1}{\overset{N}{\Σ}}\left|\right|y_i-y_j|{|}_2^2{e}^{\left(\frac{\left|\right|y_i-y_j|{|}_2^2}{\mathrm{\ρ}}\right)}$

Wherein, Q represents the canonical Weighted Constraint item of high-spectral data abundance matrix, and ∑ represents sum operation, i represents the numbering of high-spectral data abundance matrix midrange, and the span of i is { 1,2,, N}, j represent the numbering of high-spectral data abundance matrix midrange, the span of j is { 1,2 ..., N}, N represents pixel sum in high-spectral data represent the square operation of amount of orientation 2 norm, y _irepresent the i-th row in high-spectral data abundance matrix, y _jrepresent jth row in high-spectral data abundance matrix, e ⁽⁾expression take natural number as the index operation at the end, and ρ represents constraining force parameter, and the span of ρ is [0,1];

(4b) degree of rarefication of high-spectral data abundance matrix according to the following formula, is calculated:

${\mathrm{\α}}_2=\frac{1}{\sqrt{L}}\underset{l=1}{\overset{L}{\Σ}}\frac{\sqrt{N}-\left|\right|x_l|{|}_1/|\left|x_l\right|{|}_2}{\sqrt{N-1}}$

Wherein, α ₂represent the degree of rarefication of high-spectral data abundance matrix, represent and get radical sign operation, L represents the wave band number of end member in alternative end member, and ∑ represents gets sum operation, and l represents the numbering of high-spectral data line number, the span of l be 1,2 ..., L}, N represent pixel sum in high-spectral data, x _lrepresent that in high-spectral data, l is capable, || || ₁represent the operation of amount of orientation 1 norm, || || ₂represent the operation of amount of orientation 2 norm;

(4c) according to the following formula, high-spectral data abundance matrix is calculated:

$Y^{(k+1)}=\mathrm{argmin}\frac{1}{2}\left|\right|X-{\mathrm{AY}}^{\left(k\right)}|{|}_F^2+\frac{{\mathrm{\α}}_1}{2}{Q}^{\left(k\right)}+{\mathrm{\α}}_2|\left|Y^{\left(k\right)}\right|{|}_{2,1}$

Wherein, Y ^(k+1)represent the high-spectral data abundance matrix of kth+1 iteration, k represents iterations when calculating high-spectral data abundance matrix, and the span of k is { 1,2 ..., 100}, argmin represents the abundance matrix operation of getting when calculating high-spectral data abundance matrix reaches minimum value represent the square operation getting F norm, X represents high-spectral data, and A represents the end member matrix after study, Y ^(k)represent the high-spectral data abundance matrix that kth is secondary, α ₁represent the parameter of the canonical Weighted Constraint item regulating abundance matrix, α ₁value be set to 1, Q ^(k)represent the canonical Weighted Constraint item of the high-spectral data abundance matrix of kth time iteration, || || _2,1represent the sum operation of getting each column vector 2 norm in abundance matrix, α ₂represent the degree of rarefication of high-spectral data abundance matrix;

(5) reconstructed error of high-spectral data abundance matrix according to the following formula, is calculated:

$RMSE={\left(\frac{1}{P\×N}\underset{u=1,t=1}{\overset{PN}{\Σ}}{(Y_{ut}-{\hat{Y}}_{ut})}^2\right)}^{\frac{1}{2}}$

Wherein, RMSE represents the reconstructed error of high-spectral data abundance matrix, and P represents end member sum in alternative end member, N to represent in high-spectral data pixel sum, and u represents rower used when calculating RMSE, and the span of u is { 1,2 ..., P}, t represents row mark used when calculating RMSE, and the span of t is { 1,2, N}, ∑ represents gets sum operation, Y _utrepresent the capable t column element of high-spectral data abundance matrix u, represent the capable t column element of u of high-spectral data abundance matrix actual value;

(6) the mixed result of solution is exported:

Export the reconstructed error of the high-spectral data abundance matrix of separating mixed result.

The present invention compared with prior art, has the following advantages:

First, because the present invention utilizes digital spectrum storehouse in end member study, the end member matrix of the current subclass EO-1 hyperion base data after study is screened, obtain the end member matrix closest to digital spectrum storehouse, overcoming jljl different spectrum phenomenon can cause end member in actual end member and standard high-spectral data storehouse to there is certain difference, directly utilize end member in high-spectral data storehouse, high spectrum image solution can be made to mix result precision low, the problem of quality reconstruction difference, making the present invention have high spectrum image solution, to mix result precision high, the advantage that quality reconstruction is good.

Second, owing to invention introduces the Solution model of end member study, overcoming prior art sneaks out in journey in the high spectrum image solution of low signal-to-noise ratio, when using alternative iteration method to solve end member matrix and abundance matrix simultaneously, high spectrum image solution mixes result length consuming time, efficiency is low, making the present invention have high spectrum image solution, to mix result consuming time short, the advantage that efficiency is high.

Accompanying drawing explanation

Fig. 1 is process flow diagram of the present invention;

Fig. 2 is the oscillogram of the 6 kind end members of synthesis used by simulated data;

Fig. 3 is the present invention and NMF technology, L0.5NMF technology, and GLNMF technology is at the Comparative result figure of the reconstructed error value of the spectral modeling distance value added after white Gaussian noise under 25dB signal to noise ratio (S/N ratio) and high-spectral data abundance matrix;

Fig. 4 is the present invention and NMF technology, L0.5NMF technology, GLNMF technology is adding the Comparative result figure of reconstructed error value of spectral modeling distance value and high-spectral data abundance matrix under 15dB, 20dB, 25dB, 30dB, 35dB, 40dB, 45dB, 100dB signal to noise ratio (S/N ratio) after white Gaussian noise;

Fig. 5 is the line map that the present invention uses in True Data emulation;

Fig. 6 is the abundance figure of 11 kinds of mineral matters that the present invention tries to achieve in True Data emulation;

Embodiment

Below in conjunction with accompanying drawing, the present invention is described in further detail.

With reference to accompanying drawing 1, the step that the present invention realizes is described in further detail.

Step 1, input high-spectral data.

Step 2, synthesis EO-1 hyperion base data:

From digital spectrum storehouse, select all end members comprised in high-spectral data, obtain alternative end member.

With alternative end member initialization predeterminable area, obtain alternative area.

Wherein, the described concrete steps with alternative end member initialization predeterminable area are as follows:

1st step, inputs alternative end member.

2nd step, the initial value of setting number of iterations n, a n is set to 1.

3rd step, presetting a block size is the image-region of 64 × 64.

4th step, is on average divided into the region of 88 × 8, obtains cut zone by default image-region.

5th step, from alternative end member, a random selecting d end member, the region be not also initialised in initialize partition region, the value that number of iterations n adds 1, d produces at random, and the maximal value of d can not exceed the end member sum be selected.

6th step, judges whether number of iterations n is greater than 8, if so, obtains alternative area, otherwise, perform the 5th step.

Utilize Di Li Cray method, generate the Abundances of alternative area.

Be multiplied by end member corresponding in alternative area with the Abundances generated, obtain initial base data.

By low-pass filter, the high frequency signal components in the initial base data of filtering.

From initial base data, choose the pixel of Abundances lower than predetermined threshold value 0.8, obtain alternate pixel point.

With end member initialization alternate pixel points all in alternative end member, the Abundances that in the alternate pixel be initialised point, each end member is corresponding is set to the inverse of alternative end member sum, obtains intermediate base data.

In intermediate base data, add zero mean Gaussian white noise, obtain the EO-1 hyperion base data of synthesis.

Step 3, end member learns:

(3a) to preset size be the screening matrix of L × P and size is the storage matrix of 1 × P, wherein, L represents the wave band number of end member in alternative end member, P to represent in alternative end member end member sum, and the element that to be the screening matrix of L × P and size by size be in the storage matrix of 1 × P is initialized as complete zero.

(3b) according to the following formula, slickness bound term is constructed:

G＝||A|| _F

Wherein, G represents slickness bound term, and A represents end member matrix to be learned, || || _frepresent the operation of getting F norm.

(3c) according to the following formula, structural segmentation slickness bound term:

$R=\underset{l=1,p=1,i=1}{\overset{L,P,2}{\Σ}}(-e^{\left(-{\left(A_{lp}-B_i\right)}^2/\mathrm{\γ}\right)}+1)$

Wherein, R represents sectionally smooth bound term, and ∑ represents gets sum operation, and l represents rower, and the span of l is { 1,2,, L}, L represent the wave band number of end member in alternative end member, and p represents row mark, and the span of p is { 1,2 ..., P}, P represent end member sum in alternative end member, and i represents A _lpthe label of element in the Neighbourhood set of left and right, e ⁽⁾expression take natural number as the index operation at the end, A _lprepresent the capable p column element of l in end member matrix to be learned, B _irepresent A _lpi-th element in the Neighbourhood set of left and right, the span of i is that { 1,2}, γ represent constraining force parameter, and the span of γ is [0,1];

(3d) utilize K averaging method, the EO-1 hyperion base data of synthesis is carried out cluster operation, obtain the EO-1 hyperion base data after cluster.

Wherein, K averaging method concrete steps are as follows:

1st step, the EO-1 hyperion base data of input synthesis.

2nd step, a random selecting K pixel from the EO-1 hyperion base data of synthesis, the span of K be 1,2 ..., 30}, as initial cluster centre.

3rd step, the pixel of a non-cluster is chosen arbitrarily from the EO-1 hyperion base data of synthesis, pixel selected by calculating respectively with the Euclidean distance of a current K cluster centre, find out the cluster centre that Euclidean distance minimum value is corresponding, the cluster centre that the pixel of selected taking-up is corresponding with Euclidean distance minimum value is as same class data.

4th step, judges that in the EO-1 hyperion base data of synthesizing, all whether pixel all completes cluster, if so, performs the 5th step, otherwise, perform the 3rd step.

5th step, the average of each class pixel after calculating cluster, and using the average of each class pixel after cluster as the cluster centre after renewal.

6th step, according to the following formula, calculates the residual error that K cluster centre upgrades front and back respectively:

Cres _i＝||f _i-h _i|| ₂

Wherein, Cres _irepresent the i-th class cluster centre upgrade before and after residual error, the span of i be 1,2 ..., K}, K are the pixel numbers of random selecting, || || ₂represent the operation of amount of orientation 2 norm, f _irepresent the cluster centre after the i-th class renewal, h _irepresent the cluster centre before the i-th class renewal.

7th step, judges that K the cluster centre calculated upgrades maximal value in the residual error of front and back and whether be less than predetermined threshold value 0.2, if, the EO-1 hyperion base data of synthesis is considered as the data of non-cluster, performs the 8th step, otherwise, the EO-1 hyperion base data of synthesis is considered as the data of non-cluster, performs the 3rd step.

8th step, the pixel of a non-cluster is chosen arbitrarily from the EO-1 hyperion base data of synthesis, pixel selected by calculating respectively with the Euclidean distance of a current K cluster centre, find out the cluster centre that Euclidean distance minimum value is corresponding, the cluster centre that the pixel of selected taking-up is corresponding with Euclidean distance minimum value is as same class data.

9th step, judges that in the EO-1 hyperion base data of synthesizing, all whether pixel all completes cluster, if so, obtains the EO-1 hyperion base data after cluster, otherwise, perform the 8th step.

(3e) from the EO-1 hyperion base data after cluster, random selecting does not also carry out class data of end member study, obtains current subclass EO-1 hyperion base data.

(3f) according to the following formula, end member study is carried out to current subclass EO-1 hyperion base data, obtains the end member matrix of current subclass EO-1 hyperion base data:

$A^{(k+1)}=\mathrm{argmin}\frac{1}{2}\left|\right|Z-A^{\left(k\right)}{Y}^{\left(k\right)}|{|}_F^2+{\mathrm{\λ}}_1{G}^{\left(k\right)}+{\mathrm{\λ}}_2{R}^{\left(k\right)}+\frac{1}{2}{\mathrm{\λ}}_3||Y^{\left(k\right)}-D|{|}_F^2$

Wherein, A ^(k+1)represent the end member matrix of the current subclass EO-1 hyperion base data of kth+1 iteration, k represents iterations used when carrying out end member study to current subclass EO-1 hyperion base data, the span of k is { 1,2, ..., the initial value of 100}, k is set to 1, argmin represents the end member matrix manipulation got when carrying out end member study to current subclass EO-1 hyperion base data and reaching minimum value represent the square operation getting F norm, Z represents current subclass EO-1 hyperion base data, A ^(k)represent the end member matrix of the current subclass EO-1 hyperion base data that kth is secondary, Y ^(k)represent the abundance matrix of the current subclass EO-1 hyperion base data of kth time iteration, λ ₁represent the parameter regulating slickness bound term, λ ₁value be set to 0.9, G ^(k)represent the slickness bound term of kth time iteration, λ ₂represent the parameter regulating sectionally smooth bound term, λ ₂value be set to 0.1, R ^(k)represent the sectionally smooth bound term of kth time iteration, λ ₃represent balance parameters, λ ₃value be set to 1, D and represent the actual value that the abundance matrix of current subclass EO-1 hyperion base data is corresponding;

(3g) utilize digital spectrum storehouse, the end member matrix of the current subclass EO-1 hyperion base data after study is screened, obtains the end member matrix closest to digital spectrum storehouse.

Wherein, utilize digital spectrum storehouse, the concrete steps of screening the current subclass EO-1 hyperion base data end member matrix after study are as follows:

1st step, inputs the end member matrix of current subclass EO-1 hyperion base data, and size is the screening matrix of L × P and size is the storage matrix of 1 × P, and L represents the wave band number of end member in alternative end member, and P represents end member sum in alternative end member.

2nd step, to be the span of 1, n be the initial value of setting number of iterations n, a n 1,2 ..., P}, P represent end member sum in alternative end member.

3rd step, utilizes following formula, calculates the spectral modeling distance between the end member matrix of the current subclass EO-1 hyperion base data actual value corresponding with digital spectrum storehouse:

$d=\arccos \left(\frac{m_n^T{a}_n}{\left|\right|m_n\left|{|}_2\right|\left|a_n\right|{|}_2}\right)$

Wherein, d represents the spectral modeling distance between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse, and arccos () represents arc cosine operation, and T represents matrix transpose operation, m _nrepresent the n-th end member in the end member matrix of current subclass EO-1 hyperion base data, a _ncorresponding m in representative digit library of spectra _nactual value, || || ₂represent the operation of amount of orientation 2 norm, n represents current iteration number.

4th step, whether the spectral modeling distance judging between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse meets any one in replacement condition, if so, performs the 5th step, otherwise, perform the 7th step.

Described replacement condition is as follows:

Replacement condition 1: the spectral modeling distance between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse is 0.

Replacement condition 2: the spectral modeling distance between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse is less than corresponding current iteration numerical digit in storage matrix and is set up the value of element.

5th step, by the spectral modeling distance between actual value corresponding with digital spectrum storehouse for the end member matrix of current subclass EO-1 hyperion base data, leaves on position corresponding with current iteration number in storage matrix.

6th step, by the end member of corresponding current iteration number in the end member matrix of current subclass EO-1 hyperion base data, leaves on position corresponding with current iteration number in screening matrix.

7th step, judges whether current iteration number n equals end member sum in alternative end member, if so, obtains the end member matrix closest to digital spectrum storehouse and current storage matrix, otherwise, the value of current iteration number n is added 1, performs the 3rd step.

(3h) judge that the whether every class of the EO-1 hyperion base data after cluster has all carried out end member study, if so, obtain the end member matrix after learning and current storage matrix, otherwise, perform step (3e).

Step 4, according to the following formula, the canonical Weighted Constraint item of structure high-spectral data abundance matrix:

$Q=\underset{i=1,j=1}{\overset{N}{\Σ}}\left|\right|y_i-y_j|{|}_2^2{e}^{\left(\frac{\left|\right|y_i-y_j|{|}_2^2}{\mathrm{\ρ}}\right)}$

Wherein, Q represents the canonical Weighted Constraint item of high-spectral data abundance matrix, and ∑ represents sum operation, i represents the numbering of high-spectral data abundance matrix midrange, and the span of i is { 1,2,, N}, j represent the numbering of high-spectral data abundance matrix midrange, the span of j is { 1,2 ..., N}, N represents pixel sum in high-spectral data represent the square operation of amount of orientation 2 norm, y _irepresent the i-th row in high-spectral data abundance matrix, y _jrepresent jth row in high-spectral data abundance matrix, e ⁽⁾expression take natural number as the index operation at the end, and ρ represents constraining force parameter, and the span of ρ is [0,1];

According to the following formula, the degree of rarefication of high-spectral data abundance matrix is calculated:

${\mathrm{\α}}_2=\frac{1}{\sqrt{L}}\underset{l=1}{\overset{L}{\Σ}}\frac{\sqrt{N}-\left|\right|x_l|{|}_1/|\left|x_l\right|{|}_2}{\sqrt{N-1}}$

Wherein, α ₂represent the degree of rarefication of high-spectral data abundance matrix, represent and get radical sign operation, L represents the wave band number of end member in alternative end member, and ∑ represents gets sum operation, and l represents the numbering of high-spectral data line number, the span of l be 1,2 ..., L}, N represent pixel sum in high-spectral data, x _lrepresent that in high-spectral data, l is capable, || || ₁represent the operation of amount of orientation 1 norm, || || ₂represent the operation of amount of orientation 2 norm;

According to the following formula, high-spectral data abundance matrix is calculated:

$Y^{(k+1)}=\mathrm{argmin}\frac{1}{2}\left|\right|X-{\mathrm{AY}}^{\left(k\right)}|{|}_F^2+\frac{{\mathrm{\α}}_1}{2}{Q}^{\left(k\right)}+{\mathrm{\α}}_2|\left|Y^{\left(k\right)}\right|{|}_{2,1}$

Wherein, Y ^(k+1)represent the high-spectral data abundance matrix of kth+1 iteration, k represents iterations when calculating high-spectral data abundance matrix, and the span of k is { 1,2 ..., 100}, argmin represents the abundance matrix operation of getting when calculating high-spectral data abundance matrix reaches minimum value represent the square operation getting F norm, X represents high-spectral data, and A represents the end member matrix after study, Y ^(k)represent the high-spectral data abundance matrix that kth is secondary, α ₁represent the parameter regulating abundance canonical Weighted Constraint item, α ₁value be set to 1, Q ^(k)represent the canonical Weighted Constraint item of the high-spectral data abundance matrix of kth time iteration, || || _2,1represent the sum operation of getting each column vector 2 norm in abundance matrix, α ₂represent the degree of rarefication of high-spectral data abundance matrix.

Step 5, according to the following formula, calculates the reconstructed error of high-spectral data abundance matrix:

$RMSE={\left(\frac{1}{P\×N}\underset{u=1,t=1}{\overset{PN}{\Σ}}{(Y_{ut}-{\hat{Y}}_{ut})}^2\right)}^{\frac{1}{2}}$

Wherein, RMSE represents the reconstructed error of high-spectral data abundance matrix, and P represents end member sum in alternative end member, N to represent in high-spectral data pixel sum, and u represents rower used when calculating RMSE, and the span of u is { 1,2 ..., P}, t represents row mark used when calculating RMSE, and the span of t is { 1,2, N}, ∑ represents gets sum operation, Y _utrepresent the capable t column element of high-spectral data abundance matrix u, represent the capable t column element of u of high-spectral data abundance matrix actual value;

Step 6, exports and separates mixed result:

Export the reconstructed error of the high-spectral data abundance matrix of separating mixed result.

Below in conjunction with emulation experiment, effect of the present invention is described further.

1. simulated conditions

Be Inter (R) Core (TM) 2DuoT66002.20GHZ at CPU, internal memory 2G, WINDOWS7 system emulates.

2. simulated data emulation

The synthetic method of simulated data is identical with synthesizing EO-1 hyperion base data method in step 2 (Fig. 2 is the oscillogram of the 6 kinds of end members synthesized used by simulated data), this simulated data contains 64 × 64 pixels, 224 wave bands, wavelength coverage is 0.4 to 2.5 microns.Adopt spectral modeling range formula to assess the performance of the present invention (referred to as LENMF) end member study, calculate income value by spectral modeling range formula lower, illustrate that the pure substance spectra degree of correspondence in end member and standard database is higher.Adopt reconstructed error formula assess the present invention try to achieve the performance of abundance, lower by reconstructed error formulae discovery income value, illustrate that high spectrum image sparse solution mixes precision higher, it is better that sparse solution mixes effect.

Fig. 3 is the present invention and NMF technology, L0.5NMF technology, and GLNMF technology is at the fluctuation Comparative result figure of the reconstructed error value of the spectral modeling distance value added after white Gaussian noise under 25dB signal to noise ratio (S/N ratio) and high-spectral data abundance matrix.From experimental result, can find out: the result of NMF technology is the poorest, and fluctuation is also maximum, it is relevant that this and objective function lack bound term, and this point can pass through the experiment show of three kinds of technology below; The result of L0.5NMF technology and the result of GLNMF technology almost, fluctuate also little many; And by acquired results of the present invention, no matter from the angle of spectral modeling distance value, or from the angle of the reconstructed error value of high-spectral data abundance matrix, experimental result is all best, and the fluctuation of result is also minimum.Embody solution of the present invention thus and mix the advantage that precision is high, quality reconstruction is good.

Fig. 4 is the present invention and NMF technology, L0.5NMF technology, GLNMF technology is adding the linear result comparison diagram of reconstructed error value of spectral modeling distance value and high-spectral data abundance matrix under 15dB, 20dB, 25dB, 30dB, 35dB, 40dB, 45dB, 100dB signal to noise ratio (S/N ratio) after white Gaussian noise; From experimental result, can find out: compared to other three kinds of technology, result of the present invention is best, has good robustness.When signal to noise ratio (S/N ratio) is increased to 25dB from 15dB, four kinds of algorithms have obvious change, but results change of the present invention is a little bit smaller comparatively speaking, and the result of four kinds of algorithms is all tending towards smooth afterwards, and this is the cause due to noise decrease.Embody solution of the present invention thus and mix the high advantage of efficiency.

3. True Data emulation

Fig. 5 gives the present invention the line map used in True Data emulation; It is the cuprite photographed at Nevada ,Usa downstate by airborne visual light imaging spectrometer (AVIRIS).The size of image is 250 × 191, wherein containing mineral matters different in 11, they respectively: alunite, andradite, water ammonium feldspar, dumortierite, smalite 1, smalite 2, smectite, white mica, nontronite, vogesite and aspidelite.Cuprite data contain 224 wave bands, and wavelength is from 0.4 micron to 2.5 microns, and in order to better utilize this data in experiment, we will remove low signal-to-noise ratio and water vapor absorption wave band, final remaining 188 wave bands.

Table 1 spectral modeling distance value

Table 1 gives four kinds of technology in True Data emulation, calculates spectral modeling distance value and their mean value of 11 mineral matters of gained.As can be seen from Table 1, the pure substance spectra degree of correspondence in the end member that learns of the present invention and standard database is the highest.Fig. 6 is the abundance figure of 11 kinds of mineral matters that the present invention tries to achieve in True Data emulation.

Claims (4)

1., based on a high spectrum image sparse solution mixing method for end member study, comprise the steps:

(1) high-spectral data is inputted;

(2) EO-1 hyperion base data is synthesized:

(2a) from digital spectrum storehouse, select all end members comprised in high-spectral data, obtain alternative end member;

(2b) with alternative end member initialization predeterminable area, alternative area is obtained;

(2c) utilize Di Li Cray method, generate the Abundances of alternative area;

(2d) be multiplied by end member corresponding in alternative area with the Abundances generated, obtain initial base data;

(2e) by low-pass filter, the high frequency signal components in the initial base data of filtering;

(2f) from initial base data, choose the pixel of Abundances lower than predetermined threshold value 0.8, obtain alternate pixel point;

(2g) with end member initialization alternate pixel points all in alternative end member, the Abundances that in the alternate pixel be initialised point, each end member is corresponding is set to the inverse of alternative end member sum, obtains intermediate base data;

(2h) in intermediate base data, add zero mean Gaussian white noise, obtain the EO-1 hyperion base data of synthesis;

(3) end member study:

(3a) to preset size be the screening matrix of L × P and size is the storage matrix of 1 × P, wherein, L represents the wave band number of end member in alternative end member, P to represent in alternative end member end member sum, and the element that to be the screening matrix of L × P and size by size be in the storage matrix of 1 × P is initialized as complete zero;

(3b) according to the following formula, slickness bound term is constructed:

G＝||A|| _F

Wherein, G represents slickness bound term, and A represents end member matrix to be learned, || || _frepresent the operation of getting F norm;

(3c) according to the following formula, structural segmentation slickness bound term:

$R=\underset{l=1,p=1,i=1}{\overset{L,P,2}{\Σ}}(-e^{(-{\left(A_{lp}-B_i\right)}^2/\mathrm{\γ})}+1)$

Wherein, R represents sectionally smooth bound term, and Σ represents and gets sum operation, and l represents rower, and the span of l is { 1,2, ..., L}, L represent the wave band number of end member in alternative end member, and p represents row mark, and the span of p is { 1,2 ..., P}, P represent end member sum in alternative end member, and i represents A _lpthe label of element in the Neighbourhood set of left and right, e ⁽⁾expression take natural number as the index operation at the end, A _lprepresent the capable p column element of l in end member matrix to be learned, B _irepresent A _lpi-th element in the Neighbourhood set of left and right, the span of i is that { 1,2}, γ represent constraining force parameter, and the span of γ is [0,1];

(3d) utilize K averaging method, the EO-1 hyperion base data of synthesis is carried out cluster operation, obtain the EO-1 hyperion base data after cluster;

(3e) from the EO-1 hyperion base data after cluster, random selecting does not also carry out class data of end member study, obtains current subclass EO-1 hyperion base data;

(3f) according to the following formula, end member study is carried out to current subclass EO-1 hyperion base data, obtains the end member matrix of current subclass EO-1 hyperion base data:

$A^{(k+1)}=\arg \min \frac{1}{2}\left|\right|Z-A^{\left(k\right)}{Y}^{\left(k\right)}|{|}_F^2+{\mathrm{\λ}}_1{G}^{\left(k\right)}+{\mathrm{\λ}}_2{R}^{\left(k\right)}+\frac{1}{2}{\mathrm{\λ}}_3||Y^{\left(k\right)}-D|{|}_F^2$

Wherein, A ^(k+1)represent the end member matrix of the current subclass EO-1 hyperion base data of kth+1 iteration, k represents iterations used when carrying out end member study to current subclass EO-1 hyperion base data, the span of k is { 1,2, ..., the initial value of 100}, k is set to 1, argmin represents the end member matrix manipulation got when carrying out end member study to current subclass EO-1 hyperion base data and reaching minimum value represent the square operation getting F norm, Z represents current subclass EO-1 hyperion base data, A ^(k)represent the end member matrix of the current subclass EO-1 hyperion base data that kth is secondary, Y ^(k)represent the abundance matrix of the current subclass EO-1 hyperion base data of kth time iteration, λ ₁represent the parameter regulating slickness bound term, λ ₁value be set to 0.9, G ^(k)represent the slickness bound term of kth time iteration, λ ₂represent the parameter regulating sectionally smooth bound term, λ ₂value be set to 0.1, R ^(k)represent the sectionally smooth bound term of kth time iteration, λ ₃represent balance parameters, λ ₃value be set to 1, D and represent the actual value that the abundance matrix of current subclass EO-1 hyperion base data is corresponding;

(3g) utilize digital spectrum storehouse, the end member matrix of the current subclass EO-1 hyperion base data after study is screened, obtains the end member matrix closest to digital spectrum storehouse;

(3h) judge that the whether every class of the EO-1 hyperion base data after cluster has all carried out end member study, if so, obtain the end member matrix after learning and current storage matrix, otherwise, perform step (3e);

(4) high-spectral data abundance matrix is solved:

(4a) the canonical Weighted Constraint item of high-spectral data abundance matrix according to the following formula, is constructed:

$Q=\underset{i=1,j=1}{\overset{N}{\Σ}}\left|\right|y_i-y_j|{|}_2^2{e}^{\left(\frac{\left|\right|y_i-y_j|{|}_2^2}{\mathrm{\ρ}}\right)}$

Wherein, Q represents the canonical Weighted Constraint item of high-spectral data abundance matrix, and Σ represents sum operation, i represents the numbering of high-spectral data abundance matrix midrange, and the span of i is { 1,2, ..., N}, j represent the numbering of high-spectral data abundance matrix midrange, the span of j is { 1,2 ..., N}, N represents pixel sum in high-spectral data represent the square operation of amount of orientation 2 norm, y _irepresent the i-th row in high-spectral data abundance matrix, y _jrepresent jth row in high-spectral data abundance matrix, e ⁽⁾expression take natural number as the index operation at the end, and ρ represents constraining force parameter, and the span of ρ is [0,1];

(4b) degree of rarefication of high-spectral data abundance matrix according to the following formula, is calculated:

${\mathrm{\α}}_2=\frac{1}{\sqrt{L}}\underset{l=1}{\overset{L}{\Σ}}\frac{\sqrt{N}-\left|\right|x_l|{|}_1/|\left|x_l\right|{|}_2}{\sqrt{N-1}}$

Wherein, α ₂represent the degree of rarefication of high-spectral data abundance matrix, represent and get radical sign operation, L represents the wave band number of end member in alternative end member, and Σ represents and gets sum operation, and l represents the numbering of high-spectral data line number, the span of l be 1,2 ..., L}, N represent that in high-spectral data, pixel is total, x _lrepresent that in high-spectral data, l is capable, || || ₁represent the operation of amount of orientation 1 norm, || || ₂represent the operation of amount of orientation 2 norm;

(4c) according to the following formula, high-spectral data abundance matrix is calculated:

$Y^{(k+1)}=\mathrm{argmin}\frac{1}{2}\left|\right|X-{\mathrm{AY}}^{\left(k\right)}|{|}_F^2+\frac{{\mathrm{\α}}_1}{2}{Q}^{\left(k\right)}+{\mathrm{\α}}_2|\left|Y^{\left(k\right)}\right|{|}_{2,1}$

Wherein, Y ^(k+1)represent the high-spectral data abundance matrix of kth+1 iteration, k represents iterations when calculating high-spectral data abundance matrix, and the span of k is { 1,2 ..., 100}, argmin represents the abundance matrix operation of getting when calculating high-spectral data abundance matrix reaches minimum value represent the square operation getting F norm, X represents high-spectral data, and A represents the end member matrix after study, Y ^(k)represent the high-spectral data abundance matrix that kth is secondary, α ₁represent the parameter of the canonical Weighted Constraint item regulating abundance matrix, α ₁value be set to 1, Q ^(k)represent the canonical Weighted Constraint item of the high-spectral data abundance matrix of kth time iteration, || || _2,1represent the sum operation of getting each column vector 2 norm in abundance matrix, α ₂represent the degree of rarefication of high-spectral data abundance matrix;

(5) reconstructed error of high-spectral data abundance matrix according to the following formula, is calculated:

$RMSE={\left(\frac{1}{P\×N}\underset{u=1,t=1}{\overset{PN}{\Σ}}\right(Y_{ut}-{\hat{Y}}_{ut}\left)\right)}^{\frac{1}{2}}$

Wherein, RMSE represents the reconstructed error of high-spectral data abundance matrix, and P represents end member sum in alternative end member, N to represent in high-spectral data pixel sum, and u represents rower used when calculating RMSE, and the span of u is { 1,2 ..., P}, t represents row mark used when calculating RMSE, and the span of t is { 1,2, ..., N}, Σ represent and get sum operation, Y _utrepresent the capable t column element of high-spectral data abundance matrix u, represent the capable t column element of u of high-spectral data abundance matrix actual value;

(6) the mixed result of solution is exported:

Export the reconstructed error of the high-spectral data abundance matrix of separating mixed result.

2. the high spectrum image sparse solution mixing method based on end member study according to claim 1, is characterized in that: the described concrete steps with alternative end member initialization predeterminable area of step (2b) are as follows:

1st step, inputs alternative end member;

2nd step, the initial value of setting number of iterations n, a n is set to 1;

3rd step, presetting a block size is the image-region of 64 × 64;

4th step, is on average divided into the same area of 88 × 8, obtains cut zone by default image-region;

5th step, from alternative end member, a random selecting d end member, the region be not also initialised in initialize partition region, the value that number of iterations n adds 1, d produces at random, and the maximal value of d can not exceed the end member sum be selected;

6th step, judges whether number of iterations n is greater than 8, if so, obtains alternative area, otherwise, perform the 5th step.

3. the high spectrum image sparse solution mixing method based on end member study according to claim 1, is characterized in that: the concrete steps of the K averaging method described in step (3d) are as follows:

1st step, the EO-1 hyperion base data of input synthesis;

2nd step, a random selecting K pixel from the EO-1 hyperion base data of synthesis, the span of K be 1,2 ..., 30}, as initial cluster centre;

3rd step, the pixel of a non-cluster is chosen arbitrarily from the EO-1 hyperion base data of synthesis, pixel selected by calculating respectively with the Euclidean distance of a current K cluster centre, find out the cluster centre that Euclidean distance minimum value is corresponding, the cluster centre that the pixel of selected taking-up is corresponding with Euclidean distance minimum value is as same class data;

4th step, judges that in the EO-1 hyperion base data of synthesizing, all whether pixel all completes cluster, if so, performs the 5th step, otherwise, perform the 3rd step;

5th step, the average of each class pixel after calculating cluster, and using the average of each class pixel after cluster as the cluster centre after renewal;

6th step, according to the following formula, calculates the residual error that K cluster centre upgrades front and back respectively:

Cres _i＝||f _i-h _i|| ₂

Wherein, Cres _irepresent the i-th class cluster centre upgrade before and after residual error, the span of i be 1,2 ..., K}, K are the pixel numbers of random selecting, || || ₂represent the operation of amount of orientation 2 norm, f _irepresent the cluster centre after the i-th class renewal, h _irepresent the cluster centre before the i-th class renewal;

7th step, judges that K the cluster centre calculated upgrades maximal value in the residual error of front and back and whether be less than predetermined threshold value 0.2, if, the EO-1 hyperion base data of synthesis is considered as the data of non-cluster, performs the 8th step, otherwise, the EO-1 hyperion base data of synthesis is considered as the data of non-cluster, performs the 3rd step;

8th step, the pixel of a non-cluster is chosen arbitrarily from the EO-1 hyperion base data of synthesis, pixel selected by calculating respectively with the Euclidean distance of a current K cluster centre, find out the cluster centre that Euclidean distance minimum value is corresponding, the cluster centre that the pixel of selected taking-up is corresponding with Euclidean distance minimum value is as same class data;

9th step, judges that in the EO-1 hyperion base data of synthesizing, all whether pixel all completes cluster, if so, obtains the EO-1 hyperion base data after cluster, otherwise, perform the 8th step.

4. the high spectrum image sparse solution mixing method based on end member study according to claim 1, it is characterized in that: step (3g) is described utilizes digital spectrum storehouse, the concrete steps of screening the current subclass EO-1 hyperion base data end member matrix after study are as follows:

1st step, inputs the end member matrix of current subclass EO-1 hyperion base data, and size is the screening matrix of L × P and size is the storage matrix of 1 × P, and L represents the wave band number of end member in alternative end member, and P represents end member sum in alternative end member;

2nd step, to be the span of 1, n be the initial value of setting number of iterations n, a n 1,2 ..., P}, P represent that in alternative end member, end member is total;

3rd step, utilizes following formula, calculates the spectral modeling distance between the end member matrix of the current subclass EO-1 hyperion base data actual value corresponding with digital spectrum storehouse:

$d=arccos\left(\frac{m_n^T{a}_n}{\left|\right|m_n\left|{|}_2\right|\left|a_n\right|{|}_2}\right)$

Wherein, d represents the spectral modeling distance between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse, and arccos () represents arc cosine operation, and T represents matrix transpose operation, m _nrepresent the n-th end member in the end member matrix of current subclass EO-1 hyperion base data, a _ncorresponding m in representative digit library of spectra _nactual value, || || ₂represent the operation of amount of orientation 2 norm, n represents current iteration number;

4th step, whether the spectral modeling distance judging between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse meets any one in replacement condition, if so, performs the 5th step, otherwise, perform the 7th step;

Described replacement condition is as follows:

Replacement condition 1: the spectral modeling distance between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse is 0;

Replacement condition 2: the spectral modeling distance between the actual value that the end member matrix of current subclass EO-1 hyperion base data is corresponding with digital spectrum storehouse is less than corresponding current iteration numerical digit in storage matrix and is set up the value of element;

5th step, by the spectral modeling distance between actual value corresponding with digital spectrum storehouse for the end member matrix of current subclass EO-1 hyperion base data, leaves on position corresponding with current iteration number in storage matrix;

6th step, by the end member of corresponding current iteration number in the end member matrix of current subclass EO-1 hyperion base data, leaves on position corresponding with current iteration number in screening matrix;

7th step, judges whether current iteration number n equals end member sum in alternative end member, if so, obtains the end member matrix closest to digital spectrum storehouse and current storage matrix, otherwise, the value of current iteration number n is added 1, performs the 3rd step.