CN103679210A - Ground object recognition method based on hyperspectral image unmixing - Google Patents

Abstract

The invention discloses a ground object recognition method based on hyperspectral image unmixing. The method mainly solves the problem that according to an existing method, the type of a ground object which mixed pixels belong to cannot be accurately judged. The method includes the steps of inputting a hyperspectral image, and arranging the mixed pixels in the hyperspectral image into a matrix to form a data matrix; adding a constraint term composed of a manifold constraint of a data matrix, a sparsity constraint of an abundance matrix and a smoothness constraint of an end member matrix into a target function of an NMF algorithm to form a new target function; carrying out optimal unmixing on the new target function to obtain an end member matrix and an abundance matrix of the hyperspectral image after unmixing; and judging the type of the ground object of all mixed pixels in the hyperspectral image according to the end member matrix and the abundance matrix after unmixing. The method can improve precision of an end member value and an abundance value obtained through unmixing, and accordingly precision of hyperspectral image ground object recognition is improved and the method can be used for target tracking.

Description

Based on the mixed atural object recognition methods of high spectrum image solution

Technical field

The invention belongs to technical field of remote sensing image processing, is that a kind of based on the mixed atural object recognition methods of high spectrum image solution, the method can be used for the analysis of high spectrum image, and a mixed pixel point is decomposed into end member and corresponding Abundances.

Background technology

High-spectrum similarly is to utilize imaging spectrometer to tens of same ground table section and even hundreds of wave band the imaging simultaneously and 3-D view that obtains is comprised of two-dimensional space information and one dimension spectral information.Utilize these abundant spectral informations that atural object is segmented and differentiated, multi-field, be widely applied.Because the spectral resolution of the spectrum sensor of high spectrum image is high, so the spectral coverage forming is many, but the energy that each spectral coverage is accepted is little, so can only improve the floor area of accepting spectrum, reduces spatial resolution.Thereby so because the spatial resolution of high spectrum image is not high, also have the complicacy of nature atural object to form mixed pixel point.The ubiquity of mixed pixel point not only affects identification and the nicety of grading of atural object, and is that remote sensing technology is to the significant obstacle of quantification development.Therefore how effectively carrying out the decomposition of mixed pixel point is one of key issue of high-spectrum remote sensing application.

In high spectrum image, the model of mixed pixel point generally adopts line style mixture model, and its advantage is that algorithm is simple, explicit physical meaning.The mathematical procedure of this model is briefly described below: a pixel with L spectral coverage is expressed as X _ij∈ R ^{l * 1}, the end member matrix representation with P end member is M ∈ R ^{l * P}, the abundance matrix that M is corresponding is expressed as S _ij∈ R ^{p * 1}.Have: X _ij=MS _ij+ n, wherein n is noise.In actual environment, this model is subject to the restriction of two conditions: 1. M _up>=0, (1≤u≤L, 1≤p≤P) 2.

two formulas above represent that respectively it is certain that the energy of spectrum does not exist the size of the energy of negative value and mixing, can not be infinitely great.Above-mentioned model and restrictive condition all meet the mathematical model of Non-negative Matrix Factorization (Nonnegative Matrix Factorization, NMF), so can separate with NMF algorithm mixed.

The solution mixing method of the high spectrum image based on Non-negative Matrix Factorization (Nonnegative Matrix Factorization, NMF) proposing is at present all to add regular terms on the objective function of NMF algorithm.Because the method does not take into full account the characteristic of high spectrum image, make to separate mixed weak effect, thereby cause the precision of atural object identification not high.

Summary of the invention

The object of the invention is to for existing methodical deficiency, propose a kind of based on the mixed atural object recognition methods of high spectrum image solution, this method is by the flatness of the end member matrix of high spectrum image, the sparse property of abundance matrix and the stream shape of data matrix hypothesis structure combine, make high spectrum image solution mix better effects if, thereby make the precision of atural object identification higher.

The technical scheme that realizes this method is: input one panel height spectrum picture, mixed pixel in this high spectrum image is pressed to row and line up a matrix, form data matrix, stream shape constraint with data matrix, the bound term that the smoothness constraint of the sparse constraint of abundance matrix and end member matrix forms, join in the objective function of NMF algorithm, form new objective function, then this new objective function is optimized to separate and mixes, obtain end member matrix and the abundance matrix of this high spectrum image, then the end member matrix mixing out according to this solution and abundance matrix judge the atural object classification of all mixed pixel points in this high spectrum image.Concrete steps comprise as follows:

(1) input one panel height spectrum picture X ∈ R ^{m * N * L}, and by mixed pixel point X in this high spectrum image _ij∈ R ^{1 * L}by row, arrange composition data matrix Z ∈ R ^{l * B}, the row and column that wherein M and N are two dimensional image, i and be j be two dimensional image horizontal ordinate and ordinate, L is spectral coverage number, B is mixed pixel point sum in high spectrum image, B=M * N, R represents real number set;

(2) theoretical according to stream shape hypothesis, the stream shape bound term of construction data matrix Z:

$\frac{1}{2}\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2{W}_{\mathrm{ij}}=\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_i{D}_{\mathrm{ii}}-\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_j{W}_{\mathrm{ij}}$

$=\mathrm{Tr}\left({\mathrm{SYS}}^T\right)$

Z wherein _ithe i row of Z, z _jthe j row of Z, and z _iz _jk neighbour in one, S is abundance matrix, s _ibe the i row of S, W is the weight matrix of Z, W _ijan element of W,

for z _iand z _jweights, the mark of Tr () representing matrix, the transposition of T representing matrix, D is the diagonal line weight matrix of Z, D _iian element on D diagonal line, D _ii=Σ _jw _ij, Y is flow shape factor matrix, Y=D-W;

(3), according to high spectrum image imaging theory, in abundance matrix S, add L _1/2norm, obtains sparse constraint expression formula ‖ S ‖ _1/2, the sparse constraint item of usining as abundance matrix S;

(4) according to high spectrum image imaging theory, in end member matrix M, add Frobenius norm, obtain smoothness constraint expression formula

the smoothness constraint term of usining as end member matrix M;

(5) three bound terms that step (2)-(4) obtained are added the objective function of NMF algorithm to

in, to form new objective function:

$f\text{'}(M,S)=\frac{1}{2}{\left|\right|Z-\mathrm{MS}\left|\right|}_F^2+\mathrm{\α}{\left|\right|S\left|\right|}_{1/2}+\mathrm{\β}{\left|\right|M\left|\right|}_F^2+\mathrm{\γTr}\left({\mathrm{SYS}}^T\right)$

Wherein, α is the sparse constraint regular parameter of abundance matrix S, and β is the smoothness constraint regular parameter of end member matrix M, and γ is the stream shape constraint regular parameter of data matrix Z;

(6) the objective function f ' (M, S) step (5) being obtained is optimized and solves with iteration multiplication, obtains high spectrum image X ∈ R ^{m * N * L}end member matrix M and abundance matrix S;

(7) by above-mentioned high spectrum image X ∈ R ^{m * N * L}middle mixed pixel point X _ijby step (6), solve end member matrix M and the abundance vector s obtaining _irepresent, i.e. mixed pixel point X _ij=Ms _i;

(8) theoretical according to high spectrum image statistical distribution, the abundance vector s in step (7) _ito mixed pixel point X _ijcarry out the judgement of atural object classification, as max (s _i)=s _aitime, sentence mixed pixel point X _ijbelong to a class, the class label that obtains this mixed pixel point is v _ij=a, wherein max () represents the maximal value in amount of orientation, a=1,2 ..., P represents corresponding atural object classification numbering in this high spectrum image, P represents atural object classification sum in this high spectrum image, s _ais _ia element;

(9) to above-mentioned high spectrum image X ∈ R ^{m * N * L}in all for mixed pixel point the operation of step (8) carry out the judgement of atural object classification, obtain this high spectrum image X ∈ R ^{m * N * L}atural object classification matrix V ∈ R ^{m * N}.

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

1. the present invention compares the architectural characteristic that existing method has been considered high spectrum image more fully, the stream shape constraint of data matrix is joined to separate and mix, with the NMF algorithm that does not add the shape constraint that becomes a mandarin, CNMF algorithm is compared with sectionally smooth NMF algorithm, the precision of the abundance matrix obtaining is higher, under making to judge mixed pixel point, atural object classification is more accurate, thereby reaches higher atural object accuracy of identification.

2. the present invention joins the sparse constraint of the smoothness constraint of end member matrix and abundance matrix to separate simultaneously and mixes, compare with existing GLNMF algorithm, more taken into full account high spectrum image characteristic, greatly improved the precision of end member matrix, make the atural object classification calculating more approach true atural object classification, the requirement of more realistic application.

Accompanying drawing explanation

Fig. 1 is realization flow figure of the present invention;

Fig. 2 is the high spectrum image schematic diagram that the present invention uses;

Fig. 3 is the mixed curve of spectrum of end member and the contrast schematic diagram of actual spectrum curve out of the present invention and existing several algorithm solution;

Fig. 4 is average 30 error bar comparison diagrams that the spectrum angular distance SAD of the end member that mixes out of the present invention and existing several algorithm solution conciliates the root-mean-square error RMSE of the abundance of mixing out.

Embodiment

With reference to Fig. 1, specific embodiment of the invention step is as follows:

Step 1: input high spectrum image, build data matrix, obtain the real atural object classification matrix of this high spectrum image, real end member matrix and real abundance matrix.

1.1) high spectrum image as shown in Figure 2 of input, this image size is 145 * 145, has 16 class atural objects, each the mixed pixel point in this image can regard that the spectrum that the spectral information by 200 spectral coverages forms is vectorial as;

1.2) by this high spectrum image X ∈ R ^{m * N * L}middle mixed pixel point X _ij∈ R ^{1 * L}by row, arrange composition data matrix Z ∈ R ^{l * B}wherein, the row and column that M and N are two dimensional image, i and j are respectively horizontal ordinate and the ordinate of two dimensional image, L is spectral coverage number, P is atural object classification number, and B is mixed pixel point sum in high spectrum image, B=M * N, R represents real number set, L value is that 200, B value is that 21025, P value is 16 in this example;

1.3) the true classification matrix that obtains this high spectrum image is V ∈ R ^{m * N}, true end member is M ∈ R ^{l * P}, true abundance is S ∈ R ^{p * B}.

Step 2: structure constraint item:

1.1) theoretical according to stream shape hypothesis, the stream shape bound term of construction data matrix Z:

1.1a) use the i row z of heat kernel function computational data matrix Z _ij row z with data matrix Z _jweights W _ij:

$W_{\mathrm{ij}}=e^{-\frac{{\left|\right|z_i-z_j\left|\right|}^2}{\mathrm{\σ}}},$

Wherein, σ is thermonuclear parameter, and value is 1;

1.1b) use L ₂the i row s of norm calculation abundance matrix S _ij row s with abundance matrix S _jdistance:

‖s _i-s _j‖ ²；

1.1c) theoretical according to stream shape hypothesis, by step 1.1a) formula and step 1.1b) formula combine and form stream shape bound term:

$\frac{1}{2}\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2{W}_{\mathrm{ij}}=\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_i{D}_{\mathrm{ii}}-\underset{i,j=1}{\overset{B}{\mathrm{\Σ}}}{s}_i^T{s}_j{W}_{\mathrm{ij}}$

$=\mathrm{Tr}\left({\mathrm{SYS}}^T\right)$

Z wherein _iz _jk neighbour in one, W is the weight matrix of Z, the mark of Tr () representing matrix, the transposition of T representing matrix, D is the diagonal line weight matrix of Z, D _iian element on D diagonal line, D _ii=Σ _jw _ij, Y is flow shape factor matrix, Y=D-W;

1.2), according to high spectrum image imaging theory, in abundance matrix S, add L _1/2norm, obtains sparse constraint expression formula ‖ S ‖ _1/2, the sparse constraint item of usining as abundance matrix S;

1.3) according to high spectrum image imaging theory, in end member matrix M, add Frobenius norm, obtain smoothness constraint expression formula

the smoothness constraint term of usining as end member matrix M.

Step 3: construct new objective function

1.1) three bound terms that step 2 obtained are added into the objective function of NMF algorithm

in, to form new objective function:

$f\text{'}(M,S)=\frac{1}{2}{\left|\right|Z-\mathrm{MS}\left|\right|}_F^2+\mathrm{\α}{\left|\right|S\left|\right|}_{1/2}+\mathrm{\β}{\left|\right|M\left|\right|}_F^2+\mathrm{\γTr}\left({\mathrm{SYS}}^T\right),$

Wherein, α is the sparse constraint regular parameter of abundance matrix S, and value is that 2, β is the smoothness constraint regular parameter of end member matrix M, and value is that 1, γ is the stream shape constraint regular parameter of data matrix Z, and value is 0.6.

Step 4: the new objective function f ' (M, S) that step 3 is obtained is optimized and solves with iteration multiplication, obtains end member matrix M and abundance matrix S.

1.1) with random number initialization end member matrix M and abundance matrix S between [0,1], every row of S are normalized;

1.2) input end variable matrix M and abundance matrix S;

1.3) pass through the end member matrix M in new objective function f ' (M, S) and abundance matrix S differentiate, obtain respectively the computing formula of the iteration multiplication of end member matrix M and abundance matrix S:

M’=M.*(ZS ^T-βM)./MSS ^T，

$S.*(M^{T}Z+2\mathrm{\γSW})./(M^T\mathrm{MS}+\frac{1}{2}\mathrm{\α}{S}^{-\frac{1}{2}}+2\mathrm{\γSD})$

Wherein, (.*) the element multiplication of representing matrix, the element division of (. /) representing matrix;

1.4) by step 1.3) the new abundance matrix S ' that calculates of formula and new end member matrix M ', using new abundance matrix S ' and new end member matrix M ' as step 1.2) input;

1.5) repeated execution of steps 1.2)～1.4) n time altogether, end member matrix M and abundance matrix S that output finally need to be obtained, calculate and finish, and wherein, n is for carrying out number of times, and value is 3000.

Step 5: the high spectrum image X ∈ R of determining step 1 input ^{m * N * L}in the atural object classification of all mixed pixel points.

1.1) the high spectrum image X ∈ R to step 1 input ^{m * N * L}in mixed pixel point X _ijby step 3, separate the mixed end member matrix M obtaining and abundance vector s _irepresent, i.e. X _ij=Ms _i;

1.2) theoretical according to high spectrum image statistical distribution, by step 1.1) in the abundance vector s that obtains _ito mixed pixel point X _ijcarry out the judgement of atural object classification, as max (s _i)=s _aitime, sentence mixed pixel point X _ijbelong to a class, the class label that obtains this mixed pixel point is v _ij=a, wherein max () represents the maximal value in amount of orientation, a=1,2 ..., P represents corresponding atural object classification numbering in this high spectrum image, P represents atural object classification sum in this high spectrum image, s _ais _ia element;

1.3) to this high spectrum image X ∈ R ^{m * N * L}in all steps 1.2 for mixed pixel point) operation carry out the judgement of atural object classification, obtain this high spectrum image X ∈ R ^{m * N * L}the atural object classification matrix V ∈ R of solution after mixed ^{m * N}.

Effect of the present invention further illustrates by following emulation experiment:

(1) experiment simulation condition:

The high-spectrum that this experiment is used similarly is typical AVIRIS high spectrum image: take from the Indian remote sensing test site, the Indiana, USA northwestward of taking in June, 1992, atural object classification amounts to 16 classes, and the size of image is 145 * 145.Raw data has 220 spectral coverages, removes by 20 spectral coverages of noise pollution and water pollution, only retains 200 remaining spectral coverages.This experiment is Intel (R) Core (TM) i5-2450, dominant frequency 2.5GHz at CPU, inside saves as in the WINDOWS7 system of 4G and adopts software MATLAB2009a to carry out emulation.

(2) evaluation index:

1.1), for end member, with light spectral corner, apart from SAD, carry out truer end member M _twith estimation end member M _tsimilarity.The curve of spectrum of less two end members of value of spectrum angular distance SAD is more approaching.For abundance, with root-mean-square error RMSE, carry out truer abundance S _twith estimation abundance S _tbetween difference.Less two Abundances of value of root-mean-square error RMSE are more approaching.The formula of above-mentioned two evaluation criterions is respectively:

${\mathrm{SAD}}_t={\cos }^{-1}\left(\frac{M_t^T{M}_t}{\left|\right|M_t^T\left|\right|\left|\right|M_t\left|\right|}\right),$

${\mathrm{RMSE}}_t={\left(\frac{1}{B}{|S_t-S_t|}^2\right)}^{\frac{1}{2}},$

Wherein, M _tfor the t row of true end member matrix M, M _tfor estimating the t row of end member matrix M, S _tfor the t row of true abundance matrix S, S _tfor estimating the t row of abundance matrix S, t represents the classification of atural object, t=1, and 2 ..., P, P is atural object classification sum;

(3) experiment simulation content:

Experiment one

Utilizing the present invention to be optimized to separate to high spectrum image described in step (1) mixes, obtain end member matrix M and abundance matrix S, use respectively again NMF algorithm, CNMF algorithm, GLNMF algorithm and sectionally smooth NMF algorithm are optimized to separate to above-mentioned high spectrum image and mix, contrast with the mixed result obtaining of solution of the present invention, experimental result as shown in Figure 3, wherein:

Fig. 3 (a) is a kind of atural object zunyite curve of spectrum of mixing out by NMF algorithm solution and the contrast schematic diagram of the actual atural object zunyite curve of spectrum;

Fig. 3 (b) is a kind of atural object zunyite curve of spectrum of mixing out by CNMF algorithm solution and the contrast schematic diagram of the actual atural object zunyite curve of spectrum;

Fig. 3 (c) is a kind of atural object zunyite curve of spectrum of mixing out by GLNMF algorithm solution and the contrast schematic diagram of the actual atural object zunyite curve of spectrum;

Fig. 3 (d) is a kind of atural object zunyite curve of spectrum of mixing out by sectionally smooth NMF algorithm solution and the contrast schematic diagram of the actual atural object zunyite curve of spectrum;

Fig. 3 (e) is a kind of atural object zunyite curve of spectrum of mixing out by algorithm solution of the present invention and the contrast schematic diagram of the actual atural object zunyite curve of spectrum;

As seen from Figure 3, the present invention is than other existing methods, and the curve of spectrum of the atural object obtaining more approaches the curve of spectrum of true atural object atural object.

Experiment two

Utilize step (2) to testing an end member matrix M obtaining and abundance matrix S, to calculate respectively the root-mean-square error RMSE value of spectrum angular distance sad value and the abundance matrix S of end member matrix M, experimental result as shown in Figure 4, wherein:

Fig. 4 (a) is average 30 resultant errors rod comparison diagram of the spectrum angular distance sad value of the end member that mixes out of the present invention and above-mentioned 4 kinds of existing algorithm solutions;

Fig. 4 (b) is average 30 resultant errors rod comparison diagram of the root-mean-square error RMSE value of the present invention and the above-mentioned 4 kinds of algorithm solutions abundance of mixing out;

As seen from Figure 4, the present invention is than other existing methods, and the precision of separating the mixed end member value obtaining and Abundances is higher.

Experiment three

Utilize the experiment one end member matrix M obtaining and abundance matrix S to carry out classification judgement to all mixed pixel points in this high spectrum image, obtain separating mixed classification matrix V afterwards, recycling supporting vector machine SVM carries out the measurement of atural object accuracy of identification to separating the mixed classification matrix V obtaining afterwards with true classification matrix V, the present invention is measured to the atural object accuracy of identification of gained and atural object accuracy of identification that above-mentioned 4 kinds of existing methods are measured gained contrasts, as shown in table 1.

The contrast of table 1 atural object accuracy of identification numerical value index

As seen from Table 1, utilize the present invention to carry out atural object and identify the atural object accuracy of identification that the atural object accuracy of identification obtaining has significantly been better than above-mentioned 4 kinds of existing algorithms.

In sum, the present invention can increase substantially the precision of end member value and Abundances in high spectrum image, thereby improved better the precision of atural object identification, so based on high spectrum image solution, mixed atural object recognition methods has broad application prospects in high spectrum image identification field as a kind of in the present invention.

Claims (3)

1. based on the mixed atural object recognition methods of high spectrum image solution, comprise the steps:

(1) input one panel height spectrum picture X ∈ R ^{m * N * L}, and by mixed pixel point X in this high spectrum image _ij∈ R ^{1 * L}by row, arrange composition data matrix Z ∈ R ^{l * B}, the row and column that wherein M and N are two dimensional image, i and be j be two dimensional image horizontal ordinate and ordinate, L is spectral coverage number, B is mixed pixel point sum in high spectrum image, B=M * N, R represents real number set;

(2) theoretical according to stream shape hypothesis, the stream shape bound term of construction data matrix Z:

Z wherein _ithe i row of Z, z _jthe j row of Z, and z _iz _jk neighbour in one, S is abundance matrix, s _ithe i row of S, s _jthe j row of S,, W is the weight matrix of Z, W _ijan element of W,

for z _iand z _jweights, the mark of Tr () representing matrix, the transposition of T representing matrix, D is the diagonal line weight matrix of Z, D _iian element on D diagonal line, D _ii=Σ _jw _ij, Y is flow shape factor matrix, Y=D-W;

(3), according to high spectrum image imaging theory, in abundance matrix S, add L _1/2norm, obtains sparse constraint expression formula ‖ S ‖ _1/2, the sparse constraint item of usining as abundance matrix S;

(4) according to high spectrum image imaging theory, in end member matrix M, add Frobenius norm, obtain smoothness constraint expression formula

the smoothness constraint term of usining as end member matrix M;

(5) three bound terms that step (2)-(4) obtained are added the objective function of NMF algorithm to

in, to form new objective function:

Wherein, α is the sparse constraint regular parameter of abundance matrix S, and β is the smoothness constraint regular parameter of end member matrix M, and γ is the stream shape constraint regular parameter of data matrix Z;

(6) the objective function f ' (M, S) step (5) being obtained is optimized and solves with iteration multiplication, obtains high spectrum image X ∈ R ^{m * N * L}end member matrix M and abundance matrix S;

(7) by above-mentioned high spectrum image X ∈ R ^{m * N * L}middle mixed pixel point X _ijby step (6), solve end member matrix M and the abundance vector s obtaining _irepresent, i.e. mixed pixel point X _ij=Ms _i;

(8) theoretical according to high spectrum image statistical distribution, the abundance vector s in step (7) _ito mixed pixel point X _ijcarry out the judgement of atural object classification, as max (s _i)=s _aitime, sentence mixed pixel point X _ijbelong to a class, the class label that obtains this mixed pixel point is v _ij=a, wherein max () represents the maximal value in amount of orientation, a=1,2 ..., P represents corresponding atural object classification numbering in this high spectrum image, P represents atural object classification sum in this high spectrum image, s _ais _ia element;

(9) to above-mentioned high spectrum image X ∈ R ^{m * N * L}in all for mixed pixel point the operation of step (8) carry out the judgement of atural object classification, obtain this high spectrum image X ∈ R ^{m * N * L}atural object classification matrix V ∈ R ^{m * N}.

According to described in claims 1 based on the mixed atural object recognition methods of high spectrum image solution, described theoretical according to stream shape hypothesis of step (2) wherein, the stream shape bound term of construction data matrix Z, carry out as follows:

2a) use the i row z of heat kernel function computational data matrix Z _ij row z with data matrix Z _jweights W _ij:

Wherein, σ is thermonuclear parameter, and value is 1;

2b) use L ₂the i row s of norm calculation abundance matrix S _ij row s with abundance matrix S _jdistance:

‖s _i-s _j‖ ²；

2c) theoretical according to stream shape hypothesis, by step 2a) formula and step 2b) formula combine and form stream shape bound term:

wherein, B is the sum of mixed pixel point in high spectrum image.

According to described in claims 1 based on the mixed atural object recognition methods of high spectrum image solution, with iteration multiplication Optimization Solution objective function f ' (M, S), carry out as follows in wherein said step (6):

3a), with random number initialization end member matrix M and abundance matrix S between (0,1), every row of abundance matrix S are normalized; Input end variable matrix M and abundance matrix S;

3b) the abundance matrix S in objective function f ' (M, S)=0 is carried out to differentiate, obtains the iteration multiplication computing formula of abundance matrix S:

The element multiplication of (.*) representing matrix wherein, the element division of (. /) representing matrix, W is the weight matrix of Z, D is the diagonal line weight matrix of Z, Z is data matrix, the transposition of T representing matrix, α is the sparse constraint regular parameter of abundance matrix S, γ is the stream shape constraint regular parameter of data matrix Z;

3c) the end member matrix M in objective function f ' (M, S)=0 is carried out to differentiate, obtains the iteration multiplication computing formula of end member matrix M:

M.*(ZS ^T-βM)./MSS ^T

Wherein, β is the smoothness constraint regular parameter of end member matrix M;

3d) use step 3b) and step 3c) the new abundance matrix S ' that calculates of formula and new end member matrix M ', using abundance matrix S ' and new end member matrix M ' as step 3a) input;

3e) repeated execution of steps 3a)～3d) n time, end member matrix M and abundance matrix S that output finally need to be obtained, calculate and finish, and wherein, n is for carrying out number of times, and value is 3000.

Images