CN101697008A - Hyperspectral unmixing method for estimating regularized parameter automatically - Google Patents

Abstract

The invention provides an L-curve-based hyperspectral unmixing method for estimating a regularized parameter automatically. An expected target is reached by using the method in hyperspectral unmixing. The method comprises the following steps: (1) reading hyperspectral remote sensing image data V with a computer; (2) determining the number of end members in the hyperspectral image by a virtual dimensionality (CD) method; (3) initializing an end member matrix W and an abundance matrix H by a vertex component analysis (VCA) algorithm, and performing projection operations on the W and the H; and (4) determining a regularized item, determining an iteration expression of the regulated parameter alpha, the end member matrix W and the abundance matrix H, and iterating the alpha, the W and the H until simultaneous constringency of a target function and the regularized parameter. The algorithm is characterized by overcoming the unreasonable disadvantages caused by selecting the regularized parameter manually in an NMF in the prior art, utilizing a fixed point algorithm to estimate an optimal regularized parameter, and dynamically weighing a fitting degree function and the regularized item automatically so as to achieve a better unmixing effect.

Description

A kind of high spectrum of automatic estimation regularization parameter is separated mixing method

Technical field:

The high spectrum that the invention provides a kind of automatic estimation regularization parameter is separated mixing method, belongs to high-spectrum remote sensing and separates mixed field.

Background technology:

Three during the decade in the past, continuous development along with imaging spectral technology (imaging spectroscopy), the remote sensing images that the imaging spectrometer that carries on aircraft or satellite platform (imaging spectrometer) collects (remote sensing images) have comprised more and more abundanter space, radiation and spectral information, thereby provide strong means for the information extraction and the target detection of terrestrial materials.Spectral resolution is as an important indicator of development of remote sensing, develops into high spectrum, the development of forward superelevation spectrum by multispectral at present.

Because the lower spatial resolution limitations of remote sensing instrument and the complicated diversity of nature atural object, what make the spectrum reflection that obtains at single pixel point place not necessarily is a kind of characteristic of material, and may be the mixing of several different material spectrum, such feature pixel is called as mixed pixel (mixed pixel).In order to improve the precision of remote sensing application, just must solve the resolution problem of mixed pixel, make remote sensing application reach inferior pixel level by the pixel level.Enter pixel inside, mixed pixel is decomposed into different " basic composition unit ", or claim " end member " (endmember), and try to achieve the shared ratio of these solvents, or title " abundance " (abundance fractions), process (spectral unmixing) that Here it is so-called " spectrum is separated mixed " (Tong Qingxi, Zhang Bing, Zheng Lanfen. high-spectrum remote-sensing---principle, technology and application [M]. Higher Education Publishing House 2006).

Usually, spectrum is separated to mix and can be divided into two steps: (1) utilizes end member to select (endmember selection) technology to obtain the composition spectrum that is present in the high spectrum image; (2) determine each end member shared ratio in mixed pixel by mixed pixel (the mixed pixel decomposition) technology of decomposing.Because these class methods are that condition is set up with end member spectrum, therefore must be earlier to supervise or the method for non-supervision finds so-called end member spectrum.In recent years, along with Blind Signal Separation (Blind Source Separation, BSS) and nonnegative matrix decompose that (spectrum of non-supervision is separated the attention that the technology of mixing is subjected to the remote sensing scholars gradually for Nonnegative Matrix Factorization, NMF) rise of technology and development.It is meant under the complete condition of unknown of end member information, directly starts with from remote sensing images itself, according to information such as the spectral model of mixed pixel and constraint conditions, utilizes the signal processing method of non-supervision to obtain end member spectrum and component information thereof.

The spectrum of present non-supervision is separated mixed algorithm and is mainly concentrated on both direction: decompose (being NMF) based on the simplex geometry method of Convex Set Theory with based on the independent component analysis (being ICA) and the nonnegative matrix of Blind Signal Separation (being BSS) technology.

Decompose on (NMF) basis because our invention is based on nonnegative matrix, now it is done and briefly introduce.

Decompose in (NMF) in nonnegative matrix, known nonnegative matrix V seeks suitable nonnegative matrix W and H, makes:

V≈WH

Be the set V ∈ R of given L dimension data vector ^{L * K}, wherein K is the number of data sample in the set, this matrix can be similar to is decomposed into matrix W ∈ R ^{L * P}With matrix H ∈ R ^{P * K}Product.Generally speaking, P should be less than L and K.If suppose v _jAnd h _jBe column vector corresponding among matrix V and the H, then following formula can also be write as the form of column vector: v _j=Wh _jThat is to say each sample v _jCan regard the non-negative wire combination of the column vector of nonnegative matrix W approx as, combination coefficient is h _jComponent.So matrix W can be regarded as the one group of base that the data matrix V is carried out linear proximity, H then is the non-negative projection coefficient of sample set V on basic W.Usually can characterize the lot of data vector with a spot of base vector group, if we seek suitable base vector group, make its can representative data between potential structural relation, then can obtain well to approach and represent effect, this is the marrow place that nonnegative matrix is decomposed, and also is the major reason that it can be used widely.

At present based on the simplex geometry method of Convex Set Theory with based on independent component analysis (the Independent Component Analysis of Blind Signal Separation (BSS) technology, ICA) and nonnegative matrix decompose (NMF) and on mixed pixel decomposes, can obtain comparatively accurate decomposition result, but also face the problem that some need solve simultaneously.For example, simplex method often needs each end member in the tentation data all to have pure pixel, can be subjected to appreciable impact otherwise decompose precision.Decompose in (NMF) in nonnegative matrix, objective function is

$F(W,H)=\frac{1}{2}\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{(V_{\mathrm{ij}}-{\left(\mathrm{WH}\right)}_{\mathrm{ij}})}^2=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2,s.\mathrm{tW}\≥0,H\≥0$

It is simultaneously protruding for matrix W and H right and wrong, has a lot of local minimums, and this can cause the not unique of decomposition, has influenced the performance of algorithm greatly.Effectively address this problem dual mode is roughly arranged: (1) is according to characteristics and some existing prioris of practical problems, in the NMF process, add some the constraint condition of clear and definite physical significance is arranged, can effectively dwindle solution space, make the result of decomposition close to overall minimum solution, here be called the regularization nonnegative matrix decompose (Regularization NonnegativeMatrix Factorization, RNMF); (2) consider how nonnegative matrix to be decomposed to carry out effective initialization, make it near globally optimal solution.In our invention, that initialization is used is vertex component analysis (Vertex Component Analysis, VCA) method (J.M.P Nascimento and J.M.B.Dias.Vertex component analysis:a fast algorithm to uremixhyperspectral data[J] .IEEE Transactions on Geoscience and Remote Sensing, April 2005,43 (4): 898-910.), VCA belongs to the simplex geometry method based on Convex Set Theory, because its physical significance is very clear and definite and the advantage of the rapidity of algorithm, can effectively reduce iterations, make separate to overall minimum solution close.

Our inventive method mainly is to estimate automatically at the regularization parameter that adds regularization term, briefly introduces at this.Add regularization term in objective function after, objective function becomes:

G(W，H)＝||V-WH|| ²+δJ _W(W)+αJ _H(H)

Wherein δ, α are regularization parameter.

Existing document (Sen Jia and Yuntao Qian.Constrained Nonnegative Matrix Factorization forHyperspectral Unmixing[J] .IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, JANUARY 2009, VOL.47, NO.1:161-173.) added a large amount of regularization term at derivation algorithm, it is artificially given that but regularization parameter needs, and this is irrational to a great extent.

Summary of the invention:

Need artificial given irrationality at existing method regularization parameter, the high spectrum that the purpose of this invention is to provide a kind of automatic estimation regularization parameter is separated mixing method, abbreviate autoNMF as, the method can be according to the mutual balance between data fitting degree item and the regularization term, can get a suitable regularization parameter and realize separating mixed, and can obtain result preferably.The present invention is achieved by the following technical solutions:

Computer configuration: Intel (R) Pentium (R) Dual E2180@2.00GHz.

The present invention realizes under MATLAB R2008a language environment, shown in the method flow diagram of Fig. 1, after computing machine reads the high-spectrum remote sensing data, at first utilize virtual dimension (Virtual Dimensionality, VD) method (C.-I Chang and Q.Du.Estimation of number of spectrally distinct signal sources in hyperspectral imagery[J] .IEEETransactions on Geoscience andRemote Sensing, Mar 2004,42 (3): 608-619) determine the end member number, utilize vertex component analysis (VCA) algorithm initialization end member matrix W, we adopt the Moore-Penrose pseudoinverse to calculate abundance matrix H, non-negative in order to guarantee initialization value, opposite end variable matrix W and abundance matrix H carry out projection operation.Set the iteration expression formula of regularization term and regularization parameter then according to the L-curve basic thought, based on the nonnegative matrix decomposition algorithm, determine the update algorithm of end member matrix W and abundance matrix H, constantly iteration α, W, H are up to satisfying objective function and α restrains simultaneously.

The high spectrum of a kind of automatic estimation regularization parameter of the present invention is separated mixing method, abbreviates autoNMF as, and the process flow diagram of this method is seen shown in Figure 1, and it comprises the steps:

1. with calculating machine-readable fetching data.Computing machine reads the high-spectrum remote sensing data of having integrated under MATLAB R2008a environment; Data from the remote sensing images that optical spectrum imagers collects, what obtain is data cube, utilizes computing machine that the data cube data conversion is represented for matrix form, has just obtained the high-spectrum remote sensing data V that is read;

2. determine the end member number.After obtaining remote sensing image data, utilize virtual dimension (VD) method to determine end member number in the scene; Separate the problem of mixing for spectrum, primary problem is to determine end member number in the scene; People such as Chang propose to describe with " virtual dimension " notion (VD) number of end member in the high spectrum image, in our invention, will determine the number of end member in the high spectrum image by means of this method;

3. initialization end member matrix and abundance matrix.After determining the end member number, to the initialization of carrying out of the inventive method.

(1), utilizes vertex component analysis (VCA) initialization end member matrix W.

(2), for abundance matrix H, we adopt the Moore-Penrose pseudoinverse of matrix to calculate: H=pinv (W) * V.

(3), because VCA is based on the simplex geometry method of Convex Set Theory, and the end member spectrum W that extracts might be negative value, and matrix H also has the negative value element, this does not obviously satisfy nonnegative matrix and decomposes requirement to nonnegativity.We use the projection function matrix that initialization obtains to VCA to carry out projection, make it satisfy nonnegativity.

W←P _Ω[W]，H←P _Ω[H]

P _ΩThe formula that embodies of [ξ] is:

$P_{\mathrm{\&Omega;}}\left[\mathrm{\&xi;}\right]=\left\{\begin{array}{cc}\mathrm{\&xi;},& \mathrm{\&xi;}\&GreaterEqual;0\\ 0,& \mathrm{\&xi;}<0\end{array}\right.$

We mention vertex component analysis (VCA) method in background technology, it belongs to the simplex geometry method based on Convex Set Theory, though it is relatively poor that it mixes effect to separating of the high-spectral data that do not have pure pixel, but because itself physical significance is very clear and definite and the advantage of the rapidity of algorithm, we can adopt VCA to come autoNMF is carried out initialization, the calculating that can significantly reduce algorithm like this is consuming time, makes algorithm as far as possible close to the minimum solution of the overall situation;

Go on foot through above-mentioned three: the VCA algorithm extracts end member, pseudoinverse is calculated abundance matrix, projection, just can finish the initialization of autoNMF method.

4. after initialization is finished, determine that regularization term and expression formula thereof, regularization parameter, end member matrix W, abundance matrix H are iterative, bring in constant renewal in iteration α, W, H,, finish to separate mixed work up to satisfying the condition of convergence.

Determine that regularization term and expression formula thereof, regularization parameter, end member matrix W, abundance matrix H are iterative, these expression are shown determines to design before carrying out iteration, does not change in determining iterative process.

(1) determines regularization term and expression formula thereof.

Slickness is the general hypothesis of a big class physical phenomenon.It has described the things spatially consistance and the homogeney at (perhaps on a period of time interval), on image, be embodied in, pixel value between the neighbor shows clearly consistance and slickness except undergoing mutation in edge in a certain zone.Therefore, separate mixed problem for high spectrum, the slickness of abundance is an important priori, can be used for constructing regularization term and come the constrained non-negative matrix decomposition.We adopt and markov random file (Markov Random Fields, MRF) Gibbs that is correlated with of model is smooth comes regularization.

Be that regularization term is:

$J_H\left(H\right)=U\left(H\right)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{i\&Element;N_k}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_{\mathrm{ki}}\mathrm{\&psi;}(h_{\mathrm{pk}}-h_{\left(\mathrm{pk}\right)i},\mathrm{\&gamma;})$

N wherein _kBe the field pixel of pixel k, in order more comprehensively to reflect the relation between the neighborhood, we select eight neighbour structures.ω _KiBe weight factor, the weight that influences between the reflection neighbor.ψ (ξ, γ) be the potential function of certain form of variable ξ, we get document (P.Green.Bayesian reconstructions from emission tomography data using a modifiedEM algorithm[J] .IEEE Transactions on Medical Imaging, 1990, vol.9:84-93) the potential function form in:

ψ(ξ，γ)＝γlog[cosh(ξ/γ)]

(2) determine that regularization parameter, end member matrix W, abundance matrix H are iterative,

1) determines that regularization parameter is iterative.

In the nonnegative matrix (NMF) of classics, estimate through the Euclidean distance of estimating commonly used.

Estimate for Euclidean distance, objective function arranged:

$F(W,H)=\frac{1}{2}\underset{i}{\mathrm{\&Sigma;}}\underset{j}{\mathrm{\&Sigma;}}{(V_{\mathrm{ij}}-{\left(\mathrm{WH}\right)}_{\mathrm{ij}})}^2=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2,s.\mathrm{tW}\&GreaterEqual;0,H\&GreaterEqual;0$

But because objective function is simultaneously protruding for matrix W and H right and wrong, have a lot of local minimums, this can cause the not unique of decomposition, has influenced the performance of algorithm greatly.In order effectively to address this problem, can be according to characteristics and some existing prioris of practical problems, in the NMF process, add some the constraint condition of clear and definite physical significance is arranged, such as sparse property, slickness, orthogonality etc., the constraint condition that adds can effectively be dwindled solution space, make the result of decomposition close, be called regularization NMF (RNMF) here to overall minimum solution;

Add regularization term in objective function after, objective function becomes:

G(W，H)＝||V-WH|| ²+δJ _W(W)+αJ _H(H)

Carried out algorithm design at the regularization term that only comprises H in the objective function, so expression formula is:

$\min G(W,H)=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2+\mathrm{\&alpha;}{J}_H\left(H\right)$

$=F(W,H)+\mathrm{\&alpha;}{J}_H\left(H\right)$

s.f.W≥0，H≥0

Wherein (W is the data fitting degree H) to F, measures the degree of fitting of the data and the raw data of reconstruction; J _H(H) be regularization term, the regularization term utilization is abundance slickness priori in the high spectrum image in our invention.

Based on L-curve (Soontorn Oraintara, William C.Karl, David A.Castanon, Truong Q.Nguyen:A Method for Choosing the Regularization Parameter in GeneralizedTikhonov Regularized Linear Inverse Problems[C] .International Conference of ImageProcessing, 2000) in the method, be about the function of regularization parameter:

Wherein (y is the data fitting degree x) to Φ, measures the degree of fitting of the data and the raw data of reconstruction; Ψ (x) is regularization term, the priori of modeling unknown images, for example slickness, orthogonality, sparse property.According to the basic thought of L-curve, we define the fixed point method of finding the solution optimized parameter α among the autoNMF and are Wherein (W H) is objective function among the classical NMF, J to F _H(H) regularization term for adding, β is and the tangent negative slope that gets straight line of L-curve flex point.

2) determine that W and H are iterative.

In classical nonnegative matrix (NMF), update rule is:

W←W.*(VH ^T)./(WHH ^T+ε) H←H.*(W ^TV)./(W ^TWH+ε)

Add regularization term in objective function, and W and H are carried out projection operation, update rule becomes:

W←W.*P _Ω[VH ^T]./(WHH ^T+ε)

$H\&LeftArrow;H.*P_{\mathrm{\&Omega;}}[W^{T}V-{\mathrm{\&alpha;}}_H{\&dtri;}_H{J}_H\left(H\right)]./(W^T\mathrm{WH}+\mathrm{\&epsiv;})$

W and H just upgrade by above-mentioned update rule in the iterative process.

Close in order further to reduce solution space, to make to separate to the minimum solution of the overall situation, in our invention abundance matrix added and be 1 and the restriction of sparse property.

(3) iteration stopping condition.Stop when variation in the process of value in adjacent twice iteration of objective function is little, the stop condition that uses in our method is:

$\frac{|G(W^{k+1},H^{k+1})-G(W^k,H^k)|}{G(W^k,H^k)}\&le;\mathrm{tol}$

According to update algorithm α, W, H are carried out iteration, restrain simultaneously, stop iteration, finish to separate mixed work up to objective function and α.

Advantage of the present invention is: its adopts vertex component analysis (VCA) algorithm initialization, and the calculating that can significantly reduce algorithm is consuming time, abundance matrix H added and is 1 and the constraint of sparse property separating to mix in the iterative process, makes algorithm close to overall minimum solution as far as possible.Separate mixed effect in order further to improve, we use restraint to matrix W and H, have added regularization term in objective function, and the value of regularization parameter directly has influence on the rule of iteration of two matrixes, can bigger influence be arranged to separating mixed result.Method in the past all is artificial given regularization parameter, this people gets fixed parameter high spectrum to be separated to mix the result be irrational to a great extent, the nonnegative matrix decomposition method that we have proposed a kind of automatic estimation regularization parameter based on L-curve at this situation solves, abbreviate autoNMF as, the method is in the test of high spectrum emulated data and the various algorithms of True Data, have a clear superiority in respect to other algorithm, can make relative better estimation to end member spectrum and abundance distribution simultaneously.Can be generalized to big classes such as machine learning, pattern-recognition, image processing and need estimate among the finding the solution of regularization parameter problem, be with a wide range of applications and be worth.

Description of drawings:

Fig. 1 the method for the invention process flow diagram

Fig. 2 a the present invention is used for emulated data and separates mixed end member spectrum drawing for estimate

Fig. 2 b the present invention is used for emulated data and separates mixed abundance drawing for estimate

Fig. 3 a the present invention is used for Cuprite mineral data and separates mixed end member spectrum drawing for estimate

Fig. 3 b the present invention is used for Cuprite mineral data and separates mixed abundance drawing for estimate

" end member " is different " the basic composition unit " that mixed pixel decomposes in the scene, " abundance " these solvents shared ratio in scene.(Fig. 2 a, Fig. 3 are exactly to utilize the curve that the basic composition unit is drawn at the reflectivity of different wave length in the scene a) to " end member spectrum " figure, horizontal ordinate corresponding end unit wavelength (wavelength) wherein, the corresponding end member reflectivity (reflectance) of ordinate. " abundance estimations " schemed (Fig. 2 b, Fig. 3 b) and utilized solvent shared scale in scene to draw.

Embodiment:

In order to understand technical scheme of the present invention better, below embodiments of the present invention are further described:

The high spectrum of a kind of automatic estimation regularization parameter of the present invention is separated mixing method, abbreviates autoNMF as, and the process flow diagram of this method is seen shown in Figure 1, and it comprises the steps:

Step 1, with calculating machine-readable fetching data.Computing machine reads the high-spectrum remote sensing data under MATLAB R2008a environment.Data from the remote sensing images that imaging spectrometer collects, what obtain is data cube, utilizes computing machine that the data cube data conversion is represented for matrix form, has just obtained the high-spectrum remote sensing data V that is read.

Step 2, determine the end member number.After obtaining remote sensing image data, utilize virtual dimension (VD) method to determine end member number in the scene.Separate the problem of mixing for spectrum, primary problem is to determine end member number in the scene; People such as Chang propose to describe with " virtual dimension " notion (VD) number of end member in the high spectrum image, in our invention, will determine the number of end member in the high spectrum image by means of this method;

Step 3, initialization end member matrix and abundance matrix.After determining the end member number, to the initialization of carrying out of the inventive method.

(1), utilizes vertex component analysis (VCA) initialization end member matrix W.

(2), for abundance matrix H, we adopt the Moore-Penrose pseudoinverse of matrix to calculate: H=pinv (W) * V.

(3), because VCA is based on the simplex geometry method of Convex Set Theory, and the end member spectrum W that extracts might be negative value, and matrix H also has the negative value element, this does not obviously satisfy nonnegative matrix and decomposes requirement to nonnegativity.We use projection function that the matrix that initialization obtains is carried out projection, make it satisfy nonnegativity.

W←P _Ω[W]，H←P _Ω[H]

P _ΩThe formula that embodies of [ξ] is:

$P_{\mathrm{\&Omega;}}\left[\mathrm{\&xi;}\right]=\left\{\begin{array}{cc}\mathrm{\&xi;},& \mathrm{\&xi;}\&GreaterEqual;0\\ 0,& \mathrm{\&xi;}<0\end{array}\right.$

P _Ω[W] and P _ΩW in [H] and H are respectively end member matrix W and the abundance matrix H in the rule of iteration, and expression formula is

W←W.*P _Ω[VH ^T]./(WHH ^T+ε)

$H\&LeftArrow;H.*P_{\mathrm{\&Omega;}}[W^{T}V-{\mathrm{\&alpha;}}_H{\&dtri;}_H{J}_H\left(H\right)]./(W^T\mathrm{WH}+\mathrm{\&epsiv;})$

Concrete derivation is seen step 4.

Go on foot through above-mentioned three: the VCA algorithm extracts end member, pseudoinverse is calculated abundance matrix, projection, just can finish the initialization of autoNMF method.

After step 4, initialization are finished, determine that regularization term and expression formula thereof, regularization parameter, end member matrix W, abundance matrix H are iterative, bring in constant renewal in iteration α, W, H,, finish to separate mixed work up to satisfying the condition of convergence.

Determine that regularization term and expression formula thereof, regularization parameter, end member matrix W, abundance matrix H are iterative, these expression are shown determines to design before carrying out iteration, does not change in determining iterative process.

(1) determines regularization term and expression formula thereof.

Slickness is the general hypothesis of a big class physical phenomenon.It has described the things spatially consistance and the homogeney at (perhaps on a period of time interval), on image, be embodied in, pixel value between the neighbor shows clearly consistance and slickness except undergoing mutation in edge in a certain zone.Therefore, separate mixed problem for high spectrum, the slickness of abundance is an important priori, can be used for constructing regularization term and come the constrained non-negative matrix decomposition.We adopt and markov random file (Markov Random Fields, MRF) Gibbs that is correlated with of model is smooth comes regularization.

Be that regularization term is:

$J_H\left(H\right)=U\left(H\right)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{i\&Element;N_k}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_{\mathrm{ki}}\mathrm{\&psi;}(h_{\mathrm{pk}}-h_{\left(\mathrm{pk}\right)i},\mathrm{\&gamma;})$

N wherein _kBe the field pixel of pixel k, in order more comprehensively to reflect the relation between the neighborhood, we select eight neighbour structures.ω _KiBe weight factor, the weight that reflection influences between neighbor---(when i, j) neighbor up and down, we get ω in pixel _Ki=1, (i during j) diagonal pixels, gets in pixel

γ is a scale factor, gets γ=0.1; ψ (ξ, γ) be the potential function of certain form of variable ξ, we get document (P.Green.Bayesian reconstructions fromemission tomography data using a modified EM algorithm[J] .IEEE Transactions on MedicalImaging, 1990, vol.9:84-93) the potential function form in:

ψ(ξ，γ)＝γlog[cosh(ξ/γ)]

(2) determine that regularization parameter, end member matrix W, abundance matrix H are iterative,

1) determines that regularization parameter is iterative.

In the nonnegative matrix (NMF) of classics, estimate through the Euclidean distance of estimating commonly used.

Estimate for Euclidean distance, objective function arranged:

$F(W,H)=\frac{1}{2}\underset{i}{\mathrm{\&Sigma;}}\underset{j}{\mathrm{\&Sigma;}}{(V_{\mathrm{ij}}-{\left(\mathrm{WH}\right)}_{\mathrm{ij}})}^2=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2,s.\mathrm{tW}\&GreaterEqual;0,H\&GreaterEqual;0$

But because objective function is simultaneously protruding for matrix W and H right and wrong, have a lot of local minimums, this can cause the not unique of decomposition, has influenced the performance of algorithm greatly.In order effectively to address this problem, can be according to characteristics and some existing prioris of practical problems, in the NMF process, add some the constraint condition of clear and definite physical significance is arranged, such as sparse property, slickness, orthogonality etc., the constraint condition that adds can effectively be dwindled solution space, make the result of decomposition close, be called regularization NMF (RNMF) here to overall minimum solution;

Add regularization term in objective function after, objective function becomes:

G(W，H)＝||V-WH|| ²+δJ _W(W)+αJ _H(H)

Carried out algorithm design at the regularization term that only comprises H in the objective function, so expression formula is:

$\min G(W,H)=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2+\mathrm{\&alpha;}{J}_H\left(H\right)$

$=F(W,H)+\mathrm{\&alpha;}{J}_H\left(H\right)$

s.t.W≥0，H≥0

Wherein (W is the data fitting degree H) to F, measures the degree of fitting of the data and the raw data of reconstruction; JH (H) is a regularization term, and the regularization term utilization is abundance slickness priori in the high spectrum image in our invention.

In method, be about the function of regularization parameter based on L-curve:

Wherein (y is the data fitting degree x) to Φ, measures the degree of fitting of the data and the raw data of reconstruction; Ψ (x) is regularization term, the priori of modeling unknown images, for example slickness, orthogonality, sparse property.According to the basic thought of L-curve, we define the fixed point method of finding the solution optimized parameter α among the autoNMF and are

Wherein (W H) is objective function among the classical NMF, J to F _H(H) regularization term for adding, β is and the tangent negative slope that gets straight line of L-curve flex point, rule of thumb is taken as 0.1.

2) determine that W and H are iterative.

In classical nonnegative matrix (NMF), update rule is:

W←W.*(VH ^T)./(WHH ^T+ε) H←H.*(W ^TV)./(W ^TWH+ε)

Add regularization term in objective function, and W and H are carried out projection operation, update rule becomes:

W←W.*P _Ω[VH ^T]./(WHH ^T+ε)

$H\&LeftArrow;H.*P_{\mathrm{\&Omega;}}[W^{T}V-{\mathrm{\&alpha;}}_H{\&dtri;}_H{J}_H\left(H\right)]./(W^T\mathrm{WH}+\mathrm{\&epsiv;})$

W and H just upgrade by above-mentioned update rule in the iterative process.

Close in order further to reduce solution space, to make to separate to the minimum solution of the overall situation, in our invention abundance matrix added and be 1 and the restriction of sparse property.

(3) iteration stopping condition.Stop when variation in the process of value in adjacent twice iteration of objective function is little, the stop condition that uses in our method is:

$\frac{|G(W^{k+1},H^{k+1})-G(W^k,H^k)|}{G(W^k,H^k)}<\mathrm{tol}$

Get tol=10 among the present invention ^-3

According to update algorithm α, W, H are carried out iteration, restrain simultaneously, stop iteration, finish to separate mixed work up to objective function and α.

Beneficial effect:

Experimental result: in order to verify method of the present invention---the validity of autoNMF, we use emulated data and True Data to experimentize respectively, and have obtained and separated mixed effect preferably.

Emulated data is separated mixed figure as a result shown in Fig. 2 a, Fig. 2 b.

Fig. 2 a separates to mix the end member spectrum drawing for estimate that obtains.Relatively separate mixed effect for the ease of observing, we have provided comparison diagram, and wherein thick line is the standard spectrum curve, the curve of spectrum of fine rule for estimating.

Fig. 2 b separates to mix the end member abundance drawing for estimate that obtains.

True Data separate mix as a result figure as Fig. 3 a, Fig. 3 b, shown in.

Fig. 3 a Cuprite mineral data are separated and are mixed the end member spectrum drawing for estimate that obtains, and wherein thick line is the standard spectrum curve, the curve of spectrum of fine rule for estimating.

Fig. 3 b Cuprite mineral data are separated and are mixed the end member abundance drawing for estimate that obtains.

From experimental result as can be seen, we can not only estimate regularization parameter by the algorithm of invention automatically, and can obtain and separate mixed result preferably, the method can be applied to spectrum and separate mixed, pattern-recognition, machine learning or the like a big class and need estimate among the finding the solution of regularization parameter problem, has broad application prospects and is worth.

Claims (1)

1. a high spectrum of estimating regularization parameter is automatically separated mixing method, and it is characterized in that: it has following concrete steps:

Step 1, with calculating machine-readable fetching data, computing machine reads the high-spectrum remote sensing data of having integrated under MATLAB R2008a environment; Data from the remote sensing images that optical spectrum imagers collects, utilize computing machine that the data cube data conversion is represented for matrix form, just obtained the high-spectrum remote sensing data V that is read;

Step 2, determine the end member number, obtain remote sensing image data after, utilizing virtual dimension is that the VD method is determined end member number in the scene; Separate the problem of mixing for spectrum, primary problem is to determine end member number in the scene;

Step 3, initialization end member matrix W and abundance matrix H, determine the end member number after, carry out initialization:

(1), utilizing vertex component analysis is VCA initialization end member matrix W;

(2), for abundance matrix H, we adopt the Moore-Penrose pseudoinverse of matrix to calculate: H=pinv (W) * V;

(3), because VCA is based on the simplex geometry method of Convex Set Theory, and opposite end variable matrix W and abundance matrix H carry out projection, make it satisfy nonnegativity, the projection expression formula is:

W←P _Ω[W]，H←P _Ω[H]

P _ΩThe formula that embodies of [ξ] is:

$P_{\mathrm{\&Omega;}}\left[\mathrm{\&xi;}\right]=\left\{\begin{array}{cc}\mathrm{\&xi;},& \mathrm{\&xi;}\&GreaterEqual;0\\ 0,& \mathrm{\&xi;}<0\end{array}\right.$

Go on foot through above-mentioned three: the VCA algorithm extracts end member, pseudoinverse is calculated abundance matrix, projection, just finishes the initialization of autoNMF method;

After step 4, initialization are finished, determine that regularization term and expression formula thereof, regularization parameter, end member matrix W, abundance matrix H are iterative, bring in constant renewal in iteration α, W, H,, finish to separate mixed work up to satisfying the condition of convergence; The concrete practice is:

(1) determine regularization term and expression formula thereof:

To adopt with the markov random file model be that the relevant Gibbs of MRF model is smooth comes regularization;

Be that regularization term is:

$J_H\left(H\right)=U\left(H\right)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{i\&Element;N_k}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_{\mathrm{ki}}\mathrm{\&psi;}(h_{\mathrm{pk}}-h_{\left(\mathrm{pk}\right)i},\mathrm{\&gamma;})$

N wherein _kBe the field pixel of pixel k, ω _KiBe weight factor, ψ (ξ γ) is the potential function of certain form of variable ξ, its potential function form:

ψ(ξ，γ)＝γlog[cosh(ξ/γ)]

(2) determine that regularization parameter, end member matrix W, abundance matrix H are iterative,

1. determine that regularization parameter is iterative:

Nonnegative matrix at classics is among the NMF, estimates through the Euclidean distance of estimating commonly used,

Estimate for Euclidean distance, objective function arranged:

$F(W,H)=\frac{1}{2}\underset{i}{\mathrm{\&Sigma;}}\underset{j}{\mathrm{\&Sigma;}}{(V_{\mathrm{ij}}-{\left(\mathrm{WH}\right)}_{\mathrm{ij}})}^2=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2$ s.tW≥0，H≥0

Add regularization term in objective function after, objective function becomes:

G(W，H)＝||V-WH|| ²+δJ _W(W)+αJ _H(H)

Carried out algorithm design at the regularization term that only comprises H in the objective function, so expression formula is:

$\min G(W,H)=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2+\mathrm{\&alpha;}{J}_H\left(H\right)$

$=F(W,H)+\mathrm{\&alpha;}{J}_H\left(H\right)$

s.t.W≥0，H≥0

Wherein (W is the data fitting degree H) to F, measures the degree of fitting of the data and the raw data of reconstruction; J _H(H) be regularization term, its regularization term utilization be abundance slickness priori in the high spectrum image;

In method, be about the function of regularization parameter based on L-curve: Wherein (y is the data fitting degree x) to Φ, measures the degree of fitting of the data and the raw data of reconstruction; Ψ (x) is regularization term, the priori of modeling unknown images, for example slickness, orthogonality, sparse property; The fixed point method of finding the solution optimized parameter α among the definition autoNMF is

Wherein (W H) is objective function among the classical NMF, J to F _H(H) regularization term for adding, β is and the tangent negative slope that gets straight line of L-curve flex point;

2. determine that W and H are iterative:

In classical nonnegative matrix is among the NMF, and update rule is:

W←W.*(VH ^T)./(WHH ^T+ε)H←H.*(W ^TV)./(W ^TWH+ε)

Add regularization term in objective function, and W and H are carried out projection operation, update rule becomes:

W←W.*P _Ω[VH ^T]./(WHH ^T+ε)

$H\&LeftArrow;H.*P_{\mathrm{\&Omega;}}[W^{T}V-{\mathrm{\&alpha;}}_H{\&dtri;}_H{J}_H\left(H\right)]./(W^T\mathrm{WH}+\mathrm{\&epsiv;})$

W and H just upgrade by above-mentioned update rule in the iterative process;

The iteration stopping condition: stop when variation in the process of value in adjacent twice iteration of objective function is little, the stop condition that uses in this method is:

$\frac{|G(W^{k+1},H^{k+1})-G(W^k,H^k)|}{G(W^k,H^k)}<\mathrm{tol}$

According to update algorithm α, W, H are carried out iteration, restrain simultaneously, stop iteration, finish to separate mixed work up to objective function and α.

Images