CN103942787A - Spectral unmixing method based on core prototype sample analysis - Google Patents

Abstract

The invention relates to a spectral unmixing method based on core prototype sample analysis. The method includes the steps of collecting hyperspectral data to be processed, determining the parameters of the whole process, preprocessing input image data and achieving spectral unmixing on the preprocessed data through a core prototype sample analysis method. The spectral unmixing method is easy to implement, the spectral unmixing process does not need to be independently factorized into an end member extraction process and an unmixing process, the unmixing problem of the inexistence of pure end members can be solved, and the optimal end member selecting and unmixing problem of data at different mixing degrees can be solved as well. In addition, the physical meanings of an ultimate extraction result are definite, and data deciphering capacity is higher. Meanwhile, the result obtained in the method is more stable compared with a non-negative matrix factorization spectral unmixing result, and the precision is better.

Description

A kind of spectrum solution mixing method based on the sample analysis of core prototype

Technical field

The present invention relates to a kind of spectrum solution mixing method based on the sample analysis of core prototype.

Background technology

High-spectrum remote-sensing imager can obtain in the ultraviolet of electromagnetic wave spectrum, visible ray, near infrared and mid infrared region tens of to hundreds of very narrow and continuous spectral bands, and spectral resolution advantage based on high is obtained meticulous terrestrial object information.But due to the restriction of imaging space resolution and the complicated variety on earth's surface, in some pixel of the remote sensing images of imaging gained, often comprise multiple type of ground objects, therefore the ubiquitous fact of mixed pixel makes mixed pixel treatment technology be subject in recent years Chinese scholars to pay close attention to widely and study.How from high spectral mixing data, accurately to extract typical feature (end member) spectrum, and effectively carry out the ratio (abundance) that Decomposition of Mixed Pixels obtains mixing between them and become the important research direction that quantitative remote sensing is analyzed.

Line style mixture model can perform well in the ordinary course of things describing mixed pixel and form mechanism.Decomposition of Mixed Pixels algorithm based on linear model can be divided into two large classes: based on geometric method with based on statistical method.The former supposes that all data of high spectrum image are arranged in a monomorphous, and this monomorphous summit corresponding each end member distributing position respectively.But geometry method is to concentrate and find summit from data with existing, is not suitable for the data set that there is no pure pixel.Overcome this shortcoming based on statistical method, made full use of data statistics characteristic and solve end member spectrum.Non-negative Matrix Factorization receives more concern in recent years as the one in statistical method, is especially suitable for and processes the higher data of mixability.This model is comparatively flexible, can effectively extract the information such as variance, but its expression-form complexity, thus some visual informations obtain very difficult.It adopts random initializtion in addition, is easily absorbed in local extremum, unstable result.No matter be that method of geometry or Non-negative Matrix Factorization are all applicable to the mixed problem of the very applicable high mixability data solution of situation, especially Non-negative Matrix Factorization that a large amount of blended datas exist.But in actual conditions, the situation of low mixability also can exist, in nearly one or two years, there is foreign study scholar to prove, under low mixability, based on without supervision clustering, utilize cluster centre to extract end member and complete the mixed view data being more suitable in low mixability of spectrum solution.

Analyze based on appeal, it is mixed that the present invention proposes to adopt prototype method of sample analysis to realize spectrum solution, obtains having the decomposition result of physical meaning directly perceived.And use the kernel method of its expansion, dynamic adjustments nuclear parameter, changes the mapping ability of similarity between sample in feature space, finally can flexible adaptation in the mixed demand of different mixability high-spectral data spectrum solution.

Summary of the invention

The object of the present invention is to provide the mixed result of a kind of spectrum solution more stable, precision is the spectrum solution mixing method based on the sample analysis of core prototype better.

The object of the present invention is achieved like this:

(1) gather pending high-spectral data X, X ∈ R ^{m × N}, wherein M is the dimension of spectrum vector, N is the number of all pixels of data;

(2) determine the parameter of overall flow, comprise the end member number D that view data will be extracted, nuclear parameter σ is set, relaxation factor δ;

(3) to input image data pre-service: D-1 principal component, i.e. X' ∈ R before utilizing PCA dimension-reduction algorithm to extract ^{(D-1) × N};

(4) in pretreated data, adopt that to realize spectrum solution based on core prototype method of sample analysis mixed:

At data-oriented collection X' ∈ R ^{(D-1) × N}, the number that D is prototype vector, finds the main convex closure that comprises data set, and the D-1 that comprises data set dimension convex closure is:

$\arg {\min }_{C,S}D\left(X^{\′}\right|\mathrm{SCS})=\left|\right|X^{\′}-X^{\′}\mathrm{CS}{\left|\right|}_F^2$

s.t.|c _d| ₁＝1,|s _n| ₁＝1

C≥0,S≥0

Wherein d, n is respectively D, the row sequence number of N representative, C ∈ R ^{n × D}with S ∈ R ^{d × N}will obtain main convex closure;

X'C ∈ R ^{(D-1) × D}for the prototype vector matrix of decomposition estimation gained, S is abundance matrix, obtains end member battle array X'C and abundance matrix S.

When not comprising pure end member vector in mixed land cover, the D-1 that comprises data set dimension convex closure is:

arg min _C,SD(X'|SCS)

s.t.1-δ≤|c _d| ₁≤1+δ,|s _n|=1

C≥0,S≥0。

Introduce yardstick variable α _d, make | c _d| ₁=1 and 1-δ≤| α _d|≤1+ δ, the D-1 that comprises data set dimension convex closure is:

arg min _α,C,SD(X'|SCdiag(α)S)

s.t.1-δ≤|α _d| ₁≤1+δ,|c _d| ₁=1,|s _n|=1

C≥0,S≥0。

The mixed concrete steps that obtain end member battle array X'C and abundance matrix S of spectrum solution comprise:

(4.1) initialization prototype sample analysis:

(4.1.1), to search for apart from current selected some distance and highest distance position, remove first random point of selecting:

$j^{\mathrm{new}}=\arg {\max }_i\left\{\underset{j}{\mathrm{\Σ}}\right||x_i-x_j||,j\∈C\},$

Wherein, x _i∈ X', j ^newfor the position number of new selected element, C refers to current all position numbers of reconnaissance;

(4.1.2) random initializtion abundance matrix S;

(4.1.3) α=1 is set;

(4.2) with more new variables α of projection gradient method iteration, C and S;

(4.2.1) upgrade α;

Compute gradient $g_d^{\mathrm{\α}}={\mathrm{\Σ}}_{n^{\′}}[X^{\′T}{X}^{\′}{\tilde{C}}^T\mathrm{diag}\left(\mathrm{\α}\right)\tilde{S}{\tilde{S}}^T-X^{\′T}{X}^{\′}{\tilde{S}}^T{]}_{n^{\′},d}{\tilde{c}}_{n^{\′}d};$

Upgrade α ← α-u at gradient direction _αg ^α;

α is projected to suitable region alpha _d=P _α(α _d);

U _αfor linear regulation step factor;

(4.2.2) upgrade S;

Compute gradient $G^{\tilde{S}}=\mathrm{diag}\left(\mathrm{\α}\right)({\tilde{C}}^T{T}^{\′T}{X}^{\′}\tilde{C}\mathrm{diag}\left(\mathrm{\α}\right)\tilde{S}-{\tilde{C}}^T{X}^{\′T}{X}^{\′});$

Upgrade at gradient direction $s_{d,n}\←\max \{{\tilde{s}}_{d,n}-u_{\tilde{s}}(g_{d,n}^{\tilde{s}}-\underset{d^{\′}}{\mathrm{\Σ}}{g}_{d^{\′},n}^{\tilde{s}}{\tilde{s}}_{d^{\′}n}),0\};$

S is projected on single file body

(4.2.3) upgrade S;

Compute gradient $G^{\tilde{C}}=(X^{\′T}{X}^{\′}\tilde{C}\mathrm{diag}\left(\mathrm{\α}\right)\tilde{S}{\tilde{S}}^T-X^{\′T}{X}^{\′}{\tilde{S}}^T)\mathrm{diag}\left(\mathrm{\α}\right);$

Upgrade at gradient direction $c_{n,d}\←\max \{{\tilde{c}}_{n,d}-u_{\tilde{c}}(g_{n,d}^{\tilde{c}}-\underset{n^{\′}}{\mathrm{\Σ}}{g}_{n^{\′},d}^{\tilde{c}}{\tilde{c}}_{n^{\′}d}),0\};$

S is projected on single file body

Linear regulation step factor

$P_{\mathrm{\&alpha;}}\left({\mathrm{\&alpha;}}_k\right)=\left\{\begin{array}{ccc}1-\mathrm{\&delta;}& \mathrm{if}& {\mathrm{\&alpha;}}_k<1-\mathrm{\&delta;}\\ 1+\mathrm{\&delta;}& \mathrm{if}& {\mathrm{\&alpha;}}_k>1+\mathrm{\&delta;}\\ {\mathrm{\&alpha;}}_k& & \mathrm{otherwise}\end{array}\right.;$

Gradient core form is:

$g_d^{\mathrm{\&alpha;}}={\mathrm{\&Sigma;}}_{n^{\&prime;}}[K{\tilde{C}}^T\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-K{\tilde{S}}^T{]}_{n^{\&prime;},d}{\tilde{c}}_{n^{\&prime;}d},$

$G^{\tilde{S}}=\mathrm{diag}\left(\mathrm{\&alpha;}\right)({\tilde{C}}^T{T}^{\&prime;T}{X}^{\&prime;}\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}-{\tilde{C}}^{T}K),$

$G^{\tilde{C}}=(K\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-K{\tilde{S}}^T)\mathrm{diag}\left(\mathrm{\&alpha;}\right),$

K(x _i,x _j)＝exp(-(1/2σ ²)||x _i-x _j|| ²)，

(4.3) step-up error threshold value is for termination of iterations process;

(4.4) calculate output and separate mixed result:

End member battle array is E=X'C, and abundance matrix is S.

Beneficial effect of the present invention is:

The invention process is simple, and the mixed process of spectrum solution need not independently be disassembled as end member extracts and conciliate mixed two processes, can process the mixed situation of the solution existing without pure end member, and the different optimum end member of level data that mixes selects to conciliate mixed problem.In addition finally extract result physical meaning clear and definite, stronger to the decipher of data.Simultaneously the mixed result of the relative Non-negative Matrix Factorization spectrum of the method result solution is more stable, and precision is better.

Brief description of the drawings

Fig. 1 is schematic flow sheet of the present invention;

Fig. 2 is the affect schematic diagram of relaxation factor on prototype sample extraction;

Fig. 3 is generated data 1 schematic diagram;

Fig. 4 is generated data 2 schematic diagram;

Fig. 5 is the each method of generated data 1 and extracts end member result schematic diagram while getting different IPs parameter, and wherein red point represents the end member vector finally obtaining, and blue dot represents all pixel vectors;

Fig. 6 is the each method of generated data 1 and extracts the end member curve of spectrum while getting different IPs parameter;

Fig. 7 is the each method of generated data 2 and extracts end member result schematic diagram while getting different IPs parameter, and wherein red point represents the end member vector finally obtaining, and blue dot represents all pixel vectors;

Fig. 8 is the each method of generated data 2 and extracts the end member curve of spectrum while getting different IPs parameter;

Fig. 9 is the each method of True Data and extracts end member result schematic diagram while getting different IPs parameter, (a) in three kinds of colors represent three kinds, (b-d) in red point represent the end member vector finally obtaining, blue dot represents all pixel vectors;

Figure 10 is the mixed abundance figure of True Data solution.

Embodiment

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

The invention provides a kind of spectrum solution mixing method based on the sample analysis of core prototype.Specifically refer to utilize prototype method of sample analysis (Archetypal Analysis, AA), can obtain the end member matrix that reconstruct blended data battle array is approximate (the not necessarily composition of sample in data) and abundance matrix by interactive iteration and realize spectrum solution mixed (Spectral Unmixing, SU) more accurately.Meanwhile data vector inner product form in iterative process is expanded to gaussian kernel matrix and participate in computing, and by regulating gaussian kernel parameter, with the core prototype method of sample analysis (Kernel-AA) of this new expansion, realize the high spectral mixing data of different mixabilities are extracted to the desirable end member and the abundance that are more suitable for characterization data.The inventive method is defined as: core prototype sample analysis spectrum solution mixing method (KAASU).

The present invention adopts a kind of prototype method of sample analysis of similar Non-negative Matrix Factorization, utilizes its kernel form to solve the mixed demand of the high spectral mixing data of different mixabilities spectrum solution.Prototype method of sample analysis itself can estimate approximate end member matrix and abundance matrix simultaneously.In the present invention, randomization is avoided in the system of selection of initial end metaset, ensures that decomposable process has efficient constringency performance.Relaxation factor makes model can process the blended data resolution problem without pure pixel.In addition because data vector in decomposable process presents in computing with the form of inner product, after introducing kernel function, the adjusting of nuclear parameter makes the mapping ability difference to sample inner product, thereby last prototype vector is had to control action in the position of blended data space distribution.Therefore the vertex position at the corresponding monomorphous of data space according to pure end member vector distribution in height blended data, and in low mixability data, pure end member vector is positioned at the assumed condition of this class atural object inside, for the blended data of high mixability, adopt prototype method of sample analysis can obtain desirable result; For the situation of low mixability, the mixing situation that original prototype method of sample analysis can not the actual atural object of actual response, and the present invention utilizes the sample analysis of core prototype, can solve this class situation by less nuclear parameter value; Further expand, if regulate nuclear parameter to get suitable value, core prototype method of sample analysis can be processed the mixed problem of high spectral mixing data spectrum solution of different mixed types flexibly.

The present invention is the same with the Non-negative Matrix Factorization spectrum solution mixing method of widespread use, integrates end member and extracts and two processes of abundance estimation.Its complexity is consistent with Non-negative Matrix Factorization, but initial method of the present invention has been avoided the deficiency of random initializtion, has accelerating convergence and reduces the risk that searches low performance extremal vector.The utilization of kernel method simultaneously makes its actual being widely used flexibly, the most important thing is that result physical meaning is clear and definite, and decipher is strong, has great using value in mixed pixel process field.

For realizing above-mentioned goal of the invention, this patent according to flow process as shown in Figure 1 complete separate mixed.The specific embodiment of the invention is as follows:

1. read in pending high-spectral data X, X ∈ R ^{m × N}, wherein M is the dimension of spectrum vector, N is the number of all pixels of data.

2. the parameter relating in overall flow is set.Provide the end member number D that view data will be extracted, nuclear parameter σ is set, and relaxation factor δ.

3. pair input image data pre-service obtains X'.D-1 principal component, i.e. X' ∈ R before utilizing PCA dimension-reduction algorithm to extract ^{(D-1) × N}.

4. in pretreated data, adopt that to realize spectrum solution based on core prototype method of sample analysis mixed.

At data-oriented collection X' ∈ R ^{(D-1) × N}d is the prototype vector number of (also claiming the end member that extremal vector or the present invention mention), prototype method of sample analysis is used for finding the main convex closure (principal convex hull, PCH) that comprises data set, and (D-1) dimension convex closure Optimized model that it comprises data set is:

$\begin{array}{c}\arg {\min }_{C,S}D\left(X^{\&prime;}\right|\mathrm{SCS})=\left|\right|X^{\&prime;}-X^{\&prime;}\mathrm{CS}{\left|\right|}_F^2\\ s.t|c_d{|}_1=1,|s_n{|}_1=1\\ C\&GreaterEqual;0,S\&GreaterEqual;0\end{array}---\left(1\right)$

Wherein d, n is respectively D, the row sequence number of N representative, optimum C ∈ R ^{n × D}with S ∈ R ^{d × N}will obtain main convex closure.X'C ∈ R ^{(D-1) × D}be the prototype vector matrix of decomposition estimation gained, i.e. the end member matrix of spectrum solution in mixed, S is abundance matrix.| c _d| ₁=1 and C>=0 ensure together the weighted sum that the prototype vector of last gained is sample data, and | s _n| ₁=1 and S>=0 determined that n pixel is by required prototype vector weighted sum gained.

Prototype sample X'C in actual mixed land cover, may not comprise pure end member vector, so may always can not be expressed as the convex combination of available data sample.So Optimized model expands to:

arg min _C,SD(X'|SCS)

s.t.1-δ≤|c _d| ₁≤1+δ，|s _n|＝1

C≥0，S≥0 (2)

Further introduce yardstick variable α _d, make | c _d| ₁=1 and 1-δ≤| α _d|≤1+ δ, therefore objective function is rewritten as:

arg min _α,C,SD(X'|SCdiag(α)S)

s.t.1-δ≤|α _d| ₁≤1+δ，|c _d| ₁＝1，|s _n|＝1

C≥0，S≥0 (3)

Above formula is the model of this paper final optimization pass, finally obtains end member battle array X'C and abundance matrix S.The concrete operations step of optimizing is as follows:

(1) initialization prototype sample analysis

A. for initialization vector can differ greatly as far as possible, apart from each other, using apart from current selected some distance and highest distance position as search principle:

$j^{\mathrm{new}}=\arg {\max }_i\left\{\underset{j}{\mathrm{\&Sigma;}}\right||x_i-x_j||,j\&Element;C\}---\left(4\right)$

Wherein, x _i∈ X', j ^newfor the position number of new selected element, C refers to current all position numbers of reconnaissance.First random point of selecting is removed, and selects in addition a point to replace this point again.

B. random initializtion abundance matrix S.

C., α=1 is set.

(2) with more new variables α, C and S of projection gradient method iteration.

Concrete operations are as follows:

A. upgrade α

Compute gradient $g_d^{\mathrm{\&alpha;}}={\mathrm{\&Sigma;}}_{n^{\&prime;}}[X^{\&prime;T}{X}^{\&prime;}{\tilde{C}}^T\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-X^{\&prime;T}{X}^{\&prime;}{\tilde{S}}^T{]}_{n^{\&prime;},d}{\tilde{c}}_{n^{\&prime;}d};$

Upgrade α ← α-u at gradient direction _αg ^α;

α is projected to suitable region alpha _d=P _α(α _d);

Linear regulation step factor u _α.

B. upgrade S

Compute gradient $G^{\tilde{S}}=\mathrm{diag}\left(\mathrm{\&alpha;}\right)({\tilde{C}}^T{T}^{\&prime;T}{X}^{\&prime;}\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}-{\tilde{C}}^T{X}^{\&prime;T}{X}^{\&prime;});$

Upgrade at gradient direction $s_{d,n}\&LeftArrow;\max \{{\tilde{s}}_{d,n}-u_{\tilde{s}}(g_{d,n}^{\tilde{s}}-\underset{d^{\&prime;}}{\mathrm{\&Sigma;}}{g}_{d^{\&prime;},n}^{\tilde{s}}{\tilde{s}}_{d^{\&prime;}n}),0\};$

S is projected on single file body

Linear regulation step factor

C. upgrade S

Compute gradient $G^{\tilde{C}}=(X^{\&prime;T}{X}^{\&prime;}\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-X^{\&prime;T}{X}^{\&prime;}{\tilde{S}}^T)\mathrm{diag}\left(\mathrm{\&alpha;}\right);$

Upgrade at gradient direction $c_{n,d}\&LeftArrow;\max \{{\tilde{c}}_{n,d}-u_{\tilde{c}}(g_{n,d}^{\tilde{c}}-\underset{n^{\&prime;}}{\mathrm{\&Sigma;}}{g}_{n^{\&prime;},d}^{\tilde{c}}{\tilde{c}}_{n^{\&prime;}d}),0\};$

S is projected on single file body

Linear regulation step factor

It may be noted that in above-mentioned expression formula $P_{\mathrm{\&alpha;}}\left({\mathrm{\&alpha;}}_k\right)=\left\{\begin{array}{ccc}1-\mathrm{\&delta;}& \mathrm{if}& {\mathrm{\&alpha;}}_k<1-\mathrm{\&delta;}\\ 1+\mathrm{\&delta;}& \mathrm{if}& {\mathrm{\&alpha;}}_k>1+\mathrm{\&delta;}\\ {\mathrm{\&alpha;}}_k& & \mathrm{otherwise}\end{array}\right..$ δ needs careful selection, and the characteristic component extracting when too large can not be understood as prototype feature again.Accompanying drawing 2 shows that slack variable δ value can help us to extract potential pure characteristic component, although these features are not the convex combination stack of available sample.

Found out by above-mentioned iteration of variables process, renewal process depends on the inner product relation between sample, and for example form is K=X' ^tthe nuclear matrix of X'.Therefore the iterative process of prototype sample analysis can be converted into the interative computation of the nuclear matrix based on similarity relation between sample.Can be understood as in a potential unlimited Hilbert space and extract main convex closure PCH.Therefore the solution procedure based on above-mentioned prototype sample analysis model, the Kernel-AA method that this patent emphasis provides, its implementation procedure is by α in above-mentioned iterative process, in each renewal process of C and S, the replacement of gradient calculation formula is that kernel form is as follows:

$g_d^{\mathrm{\&alpha;}}={\mathrm{\&Sigma;}}_{n^{\&prime;}}[K{\tilde{C}}^T\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-K{\tilde{S}}^T{]}_{n^{\&prime;},d}{\tilde{c}}_{n^{\&prime;}d}---\left(5\right)$

$G^{\tilde{S}}=\mathrm{diag}\left(\mathrm{\&alpha;}\right)({\tilde{C}}^T{T}^{\&prime;T}{X}^{\&prime;}\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}-{\tilde{C}}^{T}K)---\left(6\right)$

$G^{\tilde{C}}=(K\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-K{\tilde{S}}^T)\mathrm{diag}\left(\mathrm{\&alpha;}\right)---\left(7\right)$

In the present invention, choosing gaussian kernel function is:

K(x _i,x _j)＝exp(-(1/2σ ²)||x _i-x _j|| ²) (8)

By regulating nuclear parameter σ, K (x-y)=exp ((1/2 σ ²) || x-y|| ²) spatial relation between data sample is portrayed to difference.Therefore be mapped to the prototype vector difference of extracting in unlimited Hilbert space, shine upon back the prototype sample position difference that original data space is extracted afterwards.Corresponding in high spectral mixing data, the position difference of the representative end member that can extract.Therefore the Decomposition of Mixed Pixels problem that Kernel-AA can corresponding different mixabilities, the selection of nuclear parameter can realize the extraction that diverse location place is represented to atural object.

(3) step-up error threshold value is for termination of iterations process.

(4) calculate output and separate mixed result.

End member battle array is E=X'C, and abundance matrix is S.

The present invention is corresponding with optimization and the constraint condition of the mixed model of spectrum solution by prototype method of sample analysis model, realizes the demand of high spectral mixing pixel processing.Adopt apart from having selected initial vector distance and principle initialization end member collection farthest, guarantee the validity of decomposable process; And utilize relaxation factor to realize the effective decomposition without pure end member atural object data, make method practical; Application gaussian kernel function obtains data core matrix, can change parameter value, realizes the choose reasonable of end member under the different mixabilities of atural object, guarantees to separate validity and the accuracy of mixed physical meaning.

The present invention selects a sample vector first at random, then selects sample distance and MAXIMUM SELECTION to go out remaining initial end member according to distance.The most random first sample vector of selecting removes, selects in addition a sample to replace this initial end member again.Guarantee that so selected initial end member vector is stable, and the apart from each other that distributes as far as possible each other, approach the summit of blended data convex closure.

The present invention utilizes relaxation factor to realize the processing without pure end member atural object data, loosen the constraint condition of original main convex closure Optimized model the inside C matrix, make finally for the not necessarily convex combination of available data of the end member collection to mixed pixel reconstruct, can effectively process actual blended data.

The present invention adopts gaussian kernel function to replace model iterative process data vector inner product form, by regulating nuclear parameter to get less numerical value, can realize the effective extraction to being positioned at class end member in low mixability data; And the mapping ability of kernel function to spatial relation between sample under bigger numerical is applicable to process the effective extraction that is positioned at convex closure summit place end member in high blended data.

In order to verify performance of the present invention, below in conjunction with Figure of description, the present invention is described in more detail:

Use respectively two width generated datas and a width True Data to carry out the emulation experiment of the mixed content of correlation spectrum solution in this patent.Experimental result, according to data prior imformation, adopts respectively the index form such as spectrum angular distance (Spectral Angle Distance, SAD), root-mean-square error (Root Mean Square Error, RMSE) and abundance figure to evaluate and separates mixed effect.

In generated data emulation experiment, in two-dimensional space Graphics Processing effect directly perceived, three end member generated datas for selecting, and with PCA dimensionality reduction to two-dimensional space processes and displays.The first width data are from U.S.Geological Survey(USGS) choose Camallite, Almandine, tri-kinds of pure end member spectrum (224 wave bands of Biotite digital spectrum storehouse, its curve of spectrum is as Fig. 3 (a)) to synthesize mixed pixel number be 10000 blended data, Fig. 3 (b) is shown as the space distribution of 100*100 size.Utilization is uniformly distributed the random number between generation (0,1), constructs different big or small stochastic matrixes and every row of matrix are normalized, and the matrix multiple forming with end member spectrum obtains mixed pixel spectrum.The manually synthetic three kinds of atural object virtual spaces of the second width data distribute as Fig. 4 (a), and taking off decimation factor is 3, and select the pure end member filling that the first width data acquisition is used to obtain blended data as Fig. 4 (b).

Can find out that by the demonstration directly perceived of the first width generated data this data mixing degree religion is high.Prototype sample analysis solution mixing method (being defined as AASU) in this patent can extract desirable end member as shown in Fig. 5 (a).Its kernel method (being defined as KAASU) also can obtain same result as Fig. 5 (b) in the time getting larger nuclear parameter value simultaneously.On the contrary, when nuclear parameter value hour, the prototype sample position of extracting changes as shown in Fig. 5 (c-d).Fig. 6 is shown as the corresponding light spectral curve that extracts end member.Table 1 has provided the corresponding end member spectrum angular distance that extracts result of said method, and separates the root-mean-square error after mixing.Shown in graph results, for height blended data, the solution of AASU is sneaked out the end member extracting in journey and is more approached true atural object end member.KAASU also can extract same accurate prototype vector in the time getting larger numerical value, and root-mean-square error can be less simultaneously.When nuclear parameter value hour, the prototype sample searching out is positioned at mixing sample inside, now prototype sample is the weighting gained of local mixing sample, can not represent desirable end member, therefore the root-mean-square error of reconstruction result is larger.Given the second width generated data mixability of showing lower by Fig. 4.Fig. 7 provides the position of extracting end member, and Fig. 8 is the corresponding curve of spectrum.Shown in Fig. 8 (a), AASU can not extract the representative atural object of the third atural object, only extracts two kinds of atural object reconstruct blended datas that purity is higher.Fig. 8 (b-c) corresponding diagram 7 (b-c) is extracted the curve of spectrum of end member, can find out and adopt KAASU can extract three kinds of object spectrums.Although there is larger error in the third object spectrum and its real spectrum curve, and the root-mean-square error of now separating mixed result also increases, but on sacrificing the basis of Image Reconstruction accuracy, its decipher to view data is correct, this view data is really synthetic by three kinds of atural objects, but not two kinds.

In true high spectrum image data experiment, data are taken from the ROSIS imaging data of the university of Pavia of North of Italy, and spectrum dimension is 103 wave bands, has 9 kinds of atural objects, in this experiment, get that wherein three kinds of atural objects are as experimental data, its true atural object distributes as shown in Figure 10 (a).Data are three kinds of atural object data, and therefore its end member extraction result is shown as shown in Fig. 9 at two-dimensional space.Fig. 9 (a) is that three classifications truly distribute, and can find out, this blended data mixability is lower, and pure atural object should be distributed in the inside that this classification atural object distributes.Based on this hypothesis, KAASU shown in Fig. 9 (d) separates mixed end member and approaches assumed condition, and the root-mean-square error numerical value providing from table 3 also can obtain consistent conclusion.Further, because the prior imformation of the pure end member of True Data Plays lacks, Figure 10 is further from separating the validity of the abundance figure aid illustration this patent algorithm that mixes result.It is pointed out that in abundance figure, darker red position represents that abundance ratio is larger, and this component represents that atural object is more at this substep; On the contrary, bluer abundance is less.The atural object that contrast Figure 10 (a) provides truly distributes, and the decipher that the abundance result that Figure 10 (d) obtains distributes to three kinds of atural objects is more accurate, and it is particularly evident that the mixed result of the second atural object and the third atural object solution embodies.And Figure 10 (b-c) is consistent with the conclusion in generated data, in the time that parameter is larger in KAASU, can obtain the mixed effect of the solution consistent with AASU method.

Above-mentioned for the present invention especially exemplified by embodiment, not in order to limit the present invention.Provided by the inventionly be equally applicable to the processing to other high spectral mixing images based on core prototype sample analysis spectrum solution mixing method.Not departing from the spirit and scope of the invention, can do a little adjustment and optimization, be as the criterion with claim with protection scope of the present invention.

The invention provides a kind of spectrum solution mixing method based on the sample analysis of core prototype.Compare existing mixed pixel spectrum solution mixing method, the present invention is directed to high spectral mixing data exists without the actual conditions consideration of difference to some extent of pure end member, mixability, and consider solution and sneak out the requirement of journey unitarity and result stability, adopt the extend type of archetypal analysis method effectively to solve appeal correspondence problem.By adopting suitable initial method, improve the impact of random initializtion on algorithm validity and stability; Adopt in addition the algorithm introduced after relaxation factor can homographic solution never to have the situation of pure end member; The employing that the most important thing is kernel method strengthens algorithm versatility in the time of solving practical problems.Finally can comprehensively, effectively and more intuitively realize accurately spectrum solution mixed.

The present invention can be widely used in high spectral mixing pixel processing aspect, a kind of method that not only can extract as end member separately, also can be used as the mixed technology of overall spectrum solution, its result superior performance, for the mixed problem of spectrum solution provides an effective common template.There is very strong applicability, meet the demand of mixed pixel processing.

The accuracy evaluation of table 1 generated data 1

Method	e 1(SAD)	e 2(SAD)	e 3(SAD)	RMSE
					AASU	0.0050	0.0018	0.0280	0.2284
KAASU(σ＝1200)	0.0050	0.0018	0.0280	0.1060
					KAASU(σ＝100)	0.0821	0.0699	0.3091	0.1420
KAASU(σ＝10)	0.1207	0.1239	0.3390	0.2076

The accuracy evaluation of table 2 generated data 2

Method	e 1(SAD)	e 2(SAD)	e 3(SAD)	RMSE
					AASU	0	0.0361	Infinitely great	0.0754
KAASU(σ＝100)	0.0007	0.0361	0.5833	0.1400
					KAASU(σ＝10)	0.0047	0.0341	0.5801	0.3170

The accuracy evaluation of table 3 True Data 1

Method	RMSE
		AASU	0.0499
KAASU(σ＝100)	0.0268
		KAASU(σ＝0.5)	0.0193

Claims (4)

1. the spectrum solution mixing method based on the sample analysis of core prototype, is characterized in that:

(1) gather pending high-spectral data X, X ∈ R ^{m × N}, wherein M is the dimension of spectrum vector, N is the number of all pixels of data;

(2) determine the parameter of overall flow, comprise the end member number D that view data will be extracted, nuclear parameter σ is set, relaxation factor δ;

(3) to input image data pre-service: D-1 principal component, i.e. X' ∈ R before utilizing PCA dimension-reduction algorithm to extract ^{(D-1) × N};

(4) in pretreated data, adopt that to realize spectrum solution based on core prototype method of sample analysis mixed:

At data-oriented collection X' ∈ R ^{(D-1) × N}, the number that D is prototype vector, finds the main convex closure that comprises data set, and the D-1 that comprises data set dimension convex closure is:

$\arg {\min }_{C,S}D\left(X^{\&prime;}\right|\mathrm{SCS})=\left|\right|X^{\&prime;}-X^{\&prime;}\mathrm{CS}{\left|\right|}_F^2$

s.t.|c _d| ₁=1,|s _n| ₁=1

C≥0,S≥0

Wherein d, n is respectively D, the row sequence number of N representative, C ∈ R ^{n × D}with S ∈ R ^{d × N}will obtain main convex closure;

X'C ∈ R ^{(D-1) × D}for the prototype vector matrix of decomposition estimation gained, S is abundance matrix, obtains end member battle array X'C and abundance matrix S.

2. a kind of spectrum solution mixing method based on the sample analysis of core prototype according to claim 1, is characterized in that: in described step (4), when not comprising pure end member vector in mixed land cover, the D-1 that comprises data set dimension convex closure is:

arg min _C,SD(X'|SCS)

s.t.1-δ≤|c _d| ₁≤1+δ,|s _n|=1

C≥0,S≥0。

3. a kind of spectrum solution mixing method based on the sample analysis of core prototype according to claim 1, is characterized in that: in described step (4), introduce yardstick variable α _d, make | c _d| ₁=1 and 1-δ≤| α _d|≤1+ δ, the D-1 that comprises data set dimension convex closure is:

arg min _α,C,SD(X'|SCdiag(α)S)

s.t.1-δ≤|α _d| ₁≤1+δ,|c _d| ₁=1,|s _n|=1

C≥0,S≥0。

4. according to a kind of spectrum solution mixing method based on the sample analysis of core prototype described in any one in claims 1 to 3, it is characterized in that: the mixed concrete steps that obtain end member battle array X'C and abundance matrix S of described spectrum solution comprise:

(4.1) initialization prototype sample analysis:

(4.1.1), to search for apart from current selected some distance and highest distance position, remove first random point of selecting:

$j^{\mathrm{new}}=\arg {\max }_i\left\{\underset{j}{\mathrm{\&Sigma;}}\right||x_i-x_j||,j\&Element;C\},$

Wherein, x _i∈ X', j ^newfor the position number of new selected element, C refers to current all position numbers of reconnaissance;

(4.1.2) random initializtion abundance matrix S;

(4.1.3) α=1 is set;

(4.2) with more new variables α of projection gradient method iteration, C and S;

(4.2.1) upgrade α;

Compute gradient $g_d^{\mathrm{\&alpha;}}={\mathrm{\&Sigma;}}_{n^{\&prime;}}[X^{\&prime;T}{X}^{\&prime;}{\tilde{C}}^T\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-X^{\&prime;T}{X}^{\&prime;}{\tilde{S}}^T{]}_{n^{\&prime;},d}{\tilde{c}}_{n^{\&prime;}d};$

Upgrade α ← α-u at gradient direction _αg ^α;

α is projected to suitable region alpha _d=P _α(α _d);

U _αfor linear regulation step factor;

(4.2.2) upgrade S;

Compute gradient $G^{\tilde{S}}=\mathrm{diag}\left(\mathrm{\&alpha;}\right)({\tilde{C}}^T{T}^{\&prime;T}{X}^{\&prime;}\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}-{\tilde{C}}^T{X}^{\&prime;T}{X}^{\&prime;});$

Upgrade at gradient direction $s_{d,n}\&LeftArrow;\max \{{\tilde{s}}_{d,n}-u_{\tilde{s}}(g_{d,n}^{\tilde{s}}-\underset{d^{\&prime;}}{\mathrm{\&Sigma;}}{g}_{d^{\&prime;},n}^{\tilde{s}}{\tilde{s}}_{d^{\&prime;}n}),0\};$

S is projected on single file body

(4.2.3) upgrade S;

Compute gradient $G^{\tilde{C}}=(X^{\&prime;T}{X}^{\&prime;}\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-X^{\&prime;T}{X}^{\&prime;}{\tilde{S}}^T)\mathrm{diag}\left(\mathrm{\&alpha;}\right);$

Upgrade at gradient direction $c_{n,d}\&LeftArrow;\max \{{\tilde{c}}_{n,d}-u_{\tilde{c}}(g_{n,d}^{\tilde{c}}-\underset{n^{\&prime;}}{\mathrm{\&Sigma;}}{g}_{n^{\&prime;},d}^{\tilde{c}}{\tilde{c}}_{n^{\&prime;}d}),0\};$

S is projected on single file body

Linear regulation step factor

$P_{\mathrm{\&alpha;}}\left({\mathrm{\&alpha;}}_k\right)=\left\{\begin{array}{ccc}1-\mathrm{\&delta;}& \mathrm{if}& {\mathrm{\&alpha;}}_k<1-\mathrm{\&delta;}\\ 1+\mathrm{\&delta;}& \mathrm{if}& {\mathrm{\&alpha;}}_k>1+\mathrm{\&delta;}\\ {\mathrm{\&alpha;}}_k& & \mathrm{otherwise}\end{array}\right.;$

Gradient core form is:

$g_d^{\mathrm{\&alpha;}}={\mathrm{\&Sigma;}}_{n^{\&prime;}}[K{\tilde{C}}^T\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-K{\tilde{S}}^T{]}_{n^{\&prime;},d}{\tilde{c}}_{n^{\&prime;}d},$

$G^{\tilde{S}}=\mathrm{diag}\left(\mathrm{\&alpha;}\right)({\tilde{C}}^T{T}^{\&prime;T}{X}^{\&prime;}\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}-{\tilde{C}}^{T}K),$

$G^{\tilde{C}}=(K\tilde{C}\mathrm{diag}\left(\mathrm{\&alpha;}\right)\tilde{S}{\tilde{S}}^T-K{\tilde{S}}^T)\mathrm{diag}\left(\mathrm{\&alpha;}\right),$

K(x _i,x _j)=exp(-(1/2σ ²)||x _i-x _j|| ²)，

(4.3) step-up error threshold value is for termination of iterations process;

(4.4) calculate output and separate mixed result:

End member battle array is E=X'C, and abundance matrix is S.