CN104331880A - Hyper-spectral mixed pixel decomposition method based on geometric spatial spectral structure information - Google Patents

Abstract

The invention belongs to the technical field of remote sensing data processing, and particularly discloses a mixed pixel decomposition method based on geometric spatial spectral structure information, so as to solve the problems that classifications of hyper-spectral image mixed pixel point surface features are not distinct and distribution is not accurate. The method comprises steps: 1) hyper-spectral data are inputted, and the data are arranged in a matrix after pretreatment; 2) a VD method is used for estimating the number of pure end members; 3) an edge contour of the image is extracted; 4) a formula for computing a spatial distance is brought forward according to the edge and the position; 5) a formula for computing a spectral distance is brought forward according to spectral statistic information; 6) a geometric spatial spectral binding term is built according to the spatial distance and the spectral distance, and the binding term is added to an NMF model; and 7) an output end member matrix and an abundance matrix are unmixed in new NMF algorithm, and the scene surface feature classifications and the distribution ratio are judged. The method is well applicable to different hyper-spectral data, and compared with the prior method, the precision of mixed pixel decomposition is improved, and great value is provided for target detection and recognition.

Description

Based on the EO-1 hyperion mixed pixel decomposition method of geometry sky spectrum structural information

Technical field

The invention belongs to high-spectrum remote sensing data processing technology field, relate to high-spectral data solution and mix, specifically a kind of EO-1 hyperion mixed pixel decomposition method based on geometry sky spectrum structural information, can be used for the research that atural object detection and target abnormality detection and sub-pix are differentiated.

Background technology

High spectrum resolution remote sensing technique is the emerging remote sensing technology of one grown up in recent ten years.High-spectrum remote sensing data has more abundant spectral information relative to traditional multispectral data, and its widespread use shows the larger potentiality of high spectrum resolution remote sensing technique.For it in military affairs, the application on industrial and civilian, corresponding treatment technology has: dimensionality reduction, target detection, and change detects, Endmember extraction, Decomposition of Mixed Pixels and classification.Although high spectrum image has higher spectral resolution, its spatial resolution is limited, thus in an instantaneous field of view, often comprises several different atural object, i.e. spectral mixing phenomenon, almost cannot identification, reduces the performance of algorithm of target detection.In order to address this problem, Decomposition of Mixed Pixels technology can be used to the spectral pattern of extraction atural object and predicts the proportional components comprising target in pixel, significant to subsequent treatment research.Usually, under the linear mixed model of hypothesis, unsupervised solution mixing method can first carry out Endmember extraction and carry out abundance estimation again in step, also can both calculate simultaneously; Also can be divided into according to the existence hypothesis of pure end member: pure end member deposits the method that there is not pure end member in method in the data and data.

Wherein, Non-negative Matrix Factorization (NMF) is that a kind of effective EO-1 hyperion solution mixes technology.Be not present in the hypothesis in data at pure end member under, this kind of method that simultaneously can calculate pure end member family curve and abundance distribution ratio is also feasible.Domestic and international researcher proposed the problem of the Non-negative Matrix Factorization solution mixed pixel of various improvement in recent years, comprise the Non-negative Matrix Factorization (MVCNMF) of minimum volume constraint, the Non-negative Matrix Factorization (SNMF) of sparse constraint, Non-negative Matrix Factorization (PNMFSC) of smoothness constraint etc.But these methods exist identification end member to the remotely-sensed data that there is not pure end member is forbidden, distribution proportion calculates coarse problem, and do not consider similar feature statistically between the contiguous similarity of the space structure of data self and spectrum, make existing method there is limitation.

Summary of the invention

The object of the invention is to the deficiency overcoming above-mentioned prior art, consider and a kind of EO-1 hyperion mixed pixel decomposition method based on geometry sky spectrum structural information is proposed to the utilization of the empty spectrum information of data, with the distribution character of clear and definite end member and spatial positional information on the impact of result, improve and separate mixed effect.

Technical scheme of the present invention is: a kind of EO-1 hyperion mixed pixel decomposition method based on geometry sky spectrum structural information, and its step comprises as follows:

(1) input high-spectral data, image size is a × b pixel, and pre-service obtains high-spectral data l is the wave band number of process, and N is pixel number, N=a × b;

(2) use the VD estimation technique, calculate pure end member number p;

(3) extract the profile information of high spectrum image, extract first principal component amount by PCA dimensionality reduction propose Y Robert operator extraction scene profile diagram n=a × b;

(4) according to the geometric space distance metric between marginal information, positional information calculation high-spectral data the formula of this space length is proposed;

If the interior central point of local space window LW that size is m × m pixel is (m _i, n _i), in the window LW of local, judge central point (m by outline line _i, n _i) and surrounding any point (m _j, n _j) whether be in the same area, any point (m in local space window LW _j, n _j) with the geometric space distance metric of central point be:

Herein be the difference between two spectral vector;

(5) spectrum spacing tolerance is calculated according to data spectrum statistical information the formula of this spectrum spacing is proposed;

In spectrum dimension, arest neighbors figure is set up to data, then two spectral vector x _i, x _jbetween geometry spectrum spacing tolerance be defined as:

$d_{i,j}^e=\left\{\begin{array}{c}{e}^{-{\left|\right|x_i-x_j\left|\right|}_2^2},\mathrm{if}{x}_i\∈\mathrm{NB}\left(x_j\right)\\ e^{-\mathrm{\γ}{\left|\right|x_i-x_j\left|\right|}_2^2},\mathrm{if}{x}_i\∉\mathrm{NB}\left(x_j\right)\end{array}\right.$

γ herein=[(m _i-m _j)+(n _i-n _j)] ²be two spectral vector positional distances, NB (x _j) be vector x _jk nearest neighbor point set;

(6) new NMF algorithm frame is constructed

By the geometry sky spectrum distance metric in step (4) and step (5) with if according to data x former in graph theory _i, x _jbetween similar then its decompose after abundance distribution vector between must be similar principle structure empty spectrum Laplce flow shape canonical:

${\left|\right|f\left|\right|}_M^a=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2/d_{i,j}^a=S^T{L}^{a}S$ With ${\left|\right|f\left|\right|}_M^e=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2/d_{i,j}^e=S^T{L}^{e}S$

L herein ^a, L ^ethe Laplacian Matrix of adjacent map, with the diagonal matrix that each node is relevant, ${\tilde{d}}_{\mathrm{ii}}^e={\mathrm{\Σ}}_{j=1}^{n}1/d_{\mathrm{ij}}^e,{\tilde{d}}_{\mathrm{ii}}^a={\mathrm{\Σ}}_{j=1}^{n}1/d_{\mathrm{ij}}^a;$

Construct new Non-negative Matrix Factorization frame-type thus, for separate mixed after the pure end member eigenmatrix that draws, then corresponding abundance distribution scaling matrices is

$f(M,S)=f_e(M,S)+{\left|\right|f\left|\right|}_M^a+{\left|\right|f\left|\right|}_M^e=\frac{1}{2}{\left|\right|X-\mathrm{MS}\left|\right|}_2^2+\mathrm{tr}\left({\mathrm{SL}}^a{S}^T\right)+\mathrm{tr}\left({\mathrm{SL}}^e{S}^T\right)$

(7) end member curve of spectrum matrix and corresponding each end member distribution abundance figure matrix is exported by new Non-negative Matrix Factorization model.

Input high-spectral data described in above-mentioned steps (1), and pre-service obtains high-spectral data carry out in accordance with the following steps:

1a) input original hyperspectral image data s is original wave band number;

1b) three-dimensional data I is arranged in 2-D data

1c) large to noise effect in remotely-sensed data wave band is won, and obtains the spectroscopic data after band selection

1d) wave band is pressed to data matrix Z, namely by row, carry out norm normalization, obtain data

The each end member distribution abundance figure matrix being exported end member curve of spectrum matrix and correspondence by new Non-negative Matrix Factorization model described in above-mentioned steps (7), is carried out in accordance with the following steps:

2a) by EO-1 hyperion sample as input data, p is the pure end member number needing to decompose;

2b) the pure end member eigenmatrix of initialization with abundance scaling matrices

2c) upgrade M and S according to the formula alternating iteration of following multiplicative rule:

M←M.*RS ^T./MSS ^T

$S\←S.*(M^{T}R+S{D}^a+S{D}^e)./(M^T\mathrm{MS}+S{\tilde{D}}^e+S{\tilde{D}}^e)$

2d) detect whether meet the condition of convergence or halt condition, namely whether iterations has reached the upper limit or has decomposed Output rusults whether reach convergence, does not meet and returns 2c).

Beneficial effect of the present invention: the present invention is by similarity measurement between research geometry sky spectrum structural information structure spatial simlanty tolerance and spectrum, and setting new Non-negative Matrix Factorization model, to carry out to data solutions mixed.Have the following advantages compared with prior art:

1) the present invention effectively make use of high spectrum image scene hollow spectrum information, designs the tolerance of space length and spectrum spacing between two kinds of spectral vector;

2) empty spectrum information incorporates in the model framework of Non-negative Matrix Factorization by the present invention, and the iteration made new advances of deriving upgrades multiplication formula and approaches objective function.

The simulation experiment result shows, the EO-1 hyperion mixed pixel decomposition method based on geometry sky spectrum structural information that the present invention proposes can be effectively applied to separate and mix, and is applied to target detection and sub-pixed mapping analysis further.

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

Accompanying drawing explanation

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

Fig. 2 is the process flow diagram of the Non-negative Matrix Factorization of the geometry sky spectrum structural information proposed in the present invention;

Fig. 3 is the gray-scale map that the present invention emulates the real EO-1 hyperion scene of employing;

Fig. 4 is that the present invention emulates the first principal component amount after to True Data dimensionality reduction;

Fig. 5 is the scene profile diagram that simulated extraction of the present invention arrives;

Fig. 6 is the abundance distribution figure that the present invention emulates each atural object obtained by new Algorithms of Non-Negative Matrix Factorization.

Embodiment

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

Step 1. pair EO-1 hyperion original mixed data carry out pre-service.

Choose the high-spectral data I that scene size is a × b pixel, wave band number is S, and pixel number is N, and N=a × b.Because real high-spectral data numerical value is very large, then need to be normalized in spectrum dimension, be beneficial to follow-up calculating.Preprocessing process not only comprises normalized, also needs noisy wave band to filter out, and obtains new high-spectral data

Concrete steps are as follows:

1a) input original hyperspectral image data s is original wave band number;

1b) three-dimensional data I is arranged in 2-D data

1c) large to noise effect in remotely-sensed data wave band is won, and obtains the spectroscopic data after band selection

1d) wave band is pressed to data matrix Z, namely by row, carry out norm normalization, obtain data

The estimation of the pure end member number p of step 2.;

" virtual dimension " (VD) estimates, the characteristic analysis method HFC based on Neyman-Pearson detection theory effectively can estimate VD.

A) calculate correlation matrix R and the covariance matrix K of high-spectral data, and make Eigenvalues Decomposition.The eigenwert of R is designated as the eigenwert of K is designated as { λ ₁>=λ ₂>=...>=λ _l.

B) putative signal source is unknown normal amount, and noise is the Gaussian noise of zero-mean.VD question variation is dualism hypothesis problem

$\begin{array}{c}{H}_0:z_l={\widehat{\mathrm{\λ}}}_l-{\mathrm{\λ}}_l=0\\ H_1:z_l={\widehat{\mathrm{\λ}}}_l-{\mathrm{\λ}}_l>0\end{array},l=\mathrm{1,2},\·\·\·,L$

If c) H ₁be true, then think that existence signal energy contributes to associated eigenvalue, otherwise then think only there is noise contribution.Under this dualism hypothesis, can every a pair eigenwert be regarded as obey following conditional probability density stochastic variable

$p_0\left(z_l\right):p_0\left(z_l\right|H_0)\≅N(0,{\mathrm{\σ}}_{z_l}^2)$

$p_1\left(z_l\right):p_1\left(z_l\right|H_1)\≅N({\mathrm{\μ}}_l,{\mathrm{\σ}}_{z_l}^2)$

Wherein: ${\mathrm{\σ}}_{z_l}^2=\mathrm{var}\left[{\widehat{\mathrm{\λ}}}_l\right]+\mathrm{var}\left[{\mathrm{\λ}}_l\right]-2\mathrm{cov}({\widehat{\mathrm{\λ}}}_l,{\mathrm{\λ}}_l)$

D) according to equation in c), definition detection probability and false-alarm probability

$p_D=\underset{t_l}{\overset{\∞}{\∫}}{p}_1\left(z\right)\mathrm{dz}$

$p_F=\underset{t_l}{\overset{\∞}{\∫}}{p}_0\left(z\right)\mathrm{dz}$

A given false-alarm probability p _f, so t _ldetermined.If then be determined with a signal energy to make contributions to this eigenwert.Respectively t is determined to each wave band _l, make H ₁the wave band number set up is virtual dimension.

Step 3. extracts the profile information of high spectrum image.

High-spectral data extracts first principal component amount Y by principal component analysis (PCA) (the princomp () function in MATLAB) dimensionality reduction, to Y Robert operator extraction scene profile

Step 4. is according to the geometric space distance metric between marginal information, positional information calculation high-spectral data

The profile diagram obtained in step 3 in, if the interior central point of local space window LW that size is m × m pixel is (m _i, n _i), outline line can mark off in forms other points whether with central point in the same area, then any point (m in window LW _j, n _j) with the geometric space distance metric of central point be

Herein be the difference between two spectral vector.

Step 5. computational geometry spectrum spacing tolerance

In spectrum dimension, arest neighbors figure is set up to data, then two spectral vector x _i, x _jbetween geometry spectrum spacing tolerance definable:

$d_{i,j}^e=\left\{\begin{array}{c}{e}^{-{\left|\right|x_i-x_j\left|\right|}_2^2},\mathrm{if}{x}_i\∈\mathrm{NB}\left(x_j\right)\\ e^{-\mathrm{\γ}{\left|\right|x_i-x_j\left|\right|}_2^2},\mathrm{if}{x}_i\∉\mathrm{NB}\left(x_j\right)\end{array}\right.$

γ herein=[(m _i-m _j)+(n _i-n _j)] ²be two spectral vector positional distances, NB (x _j) be vector x _jk nearest neighbor point set.

Step 6. constructs new Non-negative Matrix Factorization framework.

By the geometry sky spectrum distance metric in step 4 and step 5 with according to " if former data x _i, x _jbetween similar then its decompose after abundance distribution vector between must be similar " empty spectrum Laplce (Laplacian) can be constructed flow shape canonical:

${\left|\right|f\left|\right|}_M^a=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2/d_{i,j}^a=S^T{L}^{a}S$ With ${\left|\right|f\left|\right|}_M^e=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2/d_{i,j}^e=S^T{L}^{e}S$

L herein ^a, L ^ethe Laplacian Matrix of adjacent map, with the diagonal matrix that each node is relevant, ${\tilde{d}}_{\mathrm{ii}}^e={\mathrm{\Σ}}_{j=1}^{n}1/d_{\mathrm{ij}}^e,{\tilde{d}}_{\mathrm{ii}}^a={\mathrm{\Σ}}_{j=1}^{n}1/d_{\mathrm{ij}}^a;$

Construct new Non-negative Matrix Factorization frame-type thus, for separate mixed after the pure end member eigenmatrix that draws, then corresponding abundance distribution scaling matrices is

$f(M,S)=f_e(M,S)+{\left|\right|f\left|\right|}_M^a+{\left|\right|f\left|\right|}_M^e=\frac{1}{2}{\left|\right|X-\mathrm{MS}\left|\right|}_2^2+\mathrm{tr}\left({\mathrm{SL}}^a{S}^T\right)+\mathrm{tr}\left({\mathrm{SL}}^e{S}^T\right)$

Step 7. inputs data, goes out endmember spectra curve and corresponding each end member distribution abundance figure, with reference to Fig. 2 by new Non-negative Matrix Factorization model decomposition.

A) by EO-1 hyperion sample as input data, p is the pure end member number needing to decompose;

B) the pure end member eigenmatrix of initialization with abundance scaling matrices

C) M and S is upgraded according to the formula alternating iteration of following multiplicative rule:

M←M.*RS ^T./MSS ^T

$S\←S.\·*(M^{T}R+S{D}^a+S{D}^e)./(M^T\mathrm{MS}+S{\tilde{D}}^e+S{\tilde{D}}^e)$

D) whether meet the condition of convergence or halt condition, do not meet and return c)

1. experiment condition

Experiment simulation is intel pentium 4 at processor, and the MATLAB2011b software of the computing machine of internal memory 3G carries out.As shown in Figure 3, this figure is the L=220 wave band of HYDICE sensor in Washington Square region, New York to the high spectrum image chosen in experiment, and size is the pcolor of the high spectrum image of 70 × 70.In this high spectrum image, the increased surface covering in region is trees, lawn, highway, roof, pond and path.

2. experiment content

Obtain through step 1 pre-service and need to separate mixed experimental data through the dimensionality reduction of PCA, choose first principal component amount as the scene graph extracting profile, as shown in Figure 4.Fig. 5 illustrates the profile diagram that first principal component amount obtains after Robert operator extraction profile.The mixed process of solution of the Non-negative Matrix Factorization of geometry sky spectrum structural information, to obtain the pure end member characteristic curve of blended data and corresponding abundance distribution figure, Fig. 6 respectively show abundance figure (a) Chi Shui (b) trees (c) lawn (d) roof (e) street (f) path of six kinds of atural objects.

Experiment proves that the EO-1 hyperion mixed pixel decomposition method that the present invention is based on geometry sky spectrum structural information obtains good decomposition result, and abundance distribution situation and the fact are more fitted.

Therefore, have the following advantages compared with prior art:

1) the present invention effectively make use of high spectrum image scene hollow spectrum information, designs the tolerance of space length and spectrum spacing between two kinds of spectral vector;

2) empty spectrum information incorporates in the model framework of Non-negative Matrix Factorization by the present invention, and the iteration made new advances of deriving upgrades multiplication formula and approaches objective function.

The simulation experiment result shows, the EO-1 hyperion mixed pixel decomposition method based on geometry sky spectrum structural information that the present invention proposes can be effectively applied to separate and mix, and is applied to target detection and sub-pixed mapping analysis further.

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

Claims (3)

1., based on the EO-1 hyperion mixed pixel decomposition method of geometry sky spectrum structural information, it is characterized in that: comprise the following steps:

(1) input high-spectral data, image size is a × b pixel, and pre-service obtains high-spectral data l is the wave band number of process, and N is pixel number, N=a × b;

(2) use the VD estimation technique, calculate pure end member number p;

(3) extract the profile information of high spectrum image, extract first principal component amount by PCA dimensionality reduction propose Y Robert operator extraction scene profile diagram n=a × b;

(4) according to the geometric space distance metric between marginal information, positional information calculation high-spectral data the formula of this space length is proposed;

If the interior central point of local space window LW that size is m × m pixel is (m _i, n _i), in the window LW of local, judge central point (m by outline line _i, n _i) and surrounding any point (m _j, n _j) whether be in the same area, any point (m in local space window LW _j, n _j) with the geometric space distance metric of central point for:

Herein be the difference between two spectral vector;

(5) spectrum spacing tolerance is calculated according to data spectrum statistical information the formula of this spectrum spacing is proposed;

In spectrum dimension, arest neighbors figure is set up to data, then two spectral vector x _i, x _jbetween geometry spectrum spacing tolerance be defined as:

$d_{i,j}^e=\left\{\begin{array}{c}{e}^{-{\left|\right|x_i-x_j\left|\right|}_2^2},\mathrm{if}{x}_i\∈\mathrm{NB}\left(x_j\right)\\ e^{-\mathrm{\γ}{\left|\right|x_i-x_j\left|\right|}_2^2},\mathrm{if}{x}_i\∉\mathrm{NB}\left(x_j\right)\end{array}\right.$

γ herein=[(m _i-m _j)+(n _i-n _j)] ²be two spectral vector positional distances, NB (x _j) be vector x _jk nearest neighbor point set;

(6) new NMF algorithm frame is constructed

By the geometry sky spectrum distance metric in step (4) and step (5) with if according to data x former in graph theory _i, x _jbetween similar then its decompose after abundance distribution vector between must be similar principle structure empty spectrum Laplce flow shape canonical:

${\left|\right|f\left|\right|}_M^a=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2/d_{i,j}^a=S^T{L}^{a}S$ With ${\left|\right|f\left|\right|}_M^e=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|s_i-s_j\left|\right|}^2/d_{i,j}^e=S^T{L}^{e}S$

L herein ^a, L ^ethe Laplacian Matrix of adjacent map, with the diagonal matrix that each node is relevant, ${\tilde{d}}_{\mathrm{ii}}^e={\mathrm{\Σ}}_{j=1}^{n}1/d_{\mathrm{ij}}^e,{\tilde{d}}_{\mathrm{ii}}^a={\mathrm{\Σ}}_{j=1}^{n}1/d_{\mathrm{ij}}^a;$

Construct new Non-negative Matrix Factorization frame-type thus, for separate mixed after the pure end member eigenmatrix that draws, then corresponding abundance distribution scaling matrices is

$f(M,S)=f_e(M,S)+{\left|\right|f\left|\right|}_M^a+{\left|\right|f\left|\right|}_M^e=\frac{1}{2}{\left|\right|X-\mathrm{MS}\left|\right|}_2^2+\mathrm{tr}\left({\mathrm{SL}}^a{S}^T\right)+\mathrm{tr}\left({\mathrm{SL}}^e{S}^T\right)$

(7) end member curve of spectrum matrix and corresponding each end member distribution abundance figure matrix is exported by new Non-negative Matrix Factorization model.

2. the EO-1 hyperion mixed pixel decomposition method based on geometry sky spectrum structural information according to claim 1, it is characterized in that: the input high-spectral data wherein described in step (1), and pre-service obtains high-spectral data carry out in accordance with the following steps:

1a) input the hyperspectral image data that original size is a × b s is original wave band number;

1b) three-dimensional data I is arranged in 2-D data n=a × b;

1c) large to noise effect in remotely-sensed data wave band is won, and obtains the spectroscopic data after band selection l < S;

1d) wave band is pressed to data matrix Z, namely by row, carry out norm normalization, obtain data

3. the EO-1 hyperion mixed pixel decomposition method based on geometry sky spectrum structural information according to claim 1, it is characterized in that: each end member distribution abundance figure matrix being exported end member curve of spectrum matrix and correspondence by new Non-negative Matrix Factorization model wherein described in step (7), carry out in accordance with the following steps:

2a) by EO-1 hyperion sample as input data, p is the pure end member number needing to decompose;

2b) the pure end member eigenmatrix of initialization with abundance scaling matrices

2c) upgrade M and S according to the formula alternating iteration of following multiplicative rule:

M←M.*RS ^T./MSS ^T

$S\&LeftArrow;S.*(M^{T}R+{\mathrm{SD}}^a+{\mathrm{SD}}^e)./(M^T\mathrm{MS}+S{\tilde{D}}^e+S{\tilde{D}}^e)$

2d) detect whether meet the condition of convergence or halt condition, namely whether iterations has reached the upper limit or has decomposed Output rusults whether reach convergence, does not meet and returns 2c).