A kind of high spectrum image solution mixing method returned based on weighting joint sparse
Technical field
The invention provides a kind of high spectrum image solution mixing method returned based on weighting joint sparse, belong to high spectrum image analysis field.
Background technology
High-spectrum remote-sensing plays an important role in the observation to earth land, ocean, air, but due to the lower and complicated atural object distribution character of high light spectrum image-forming spectrometer spatial resolution, there is a large amount of mixed pixels in high spectrum image.In order to make full use of the information in hyperspectral image data, need to be by the set of some constituents (also referred to as end member) and their relative scale (also referred to as abundance) by each Decomposition of Mixed Pixels.
Two ten years in the past, have many high spectrum image solutions to mix algorithm in succession to be proposed, they are mostly based on linear mixed model, the spectrum signal of the high spectrum image that this model hypothesis is gathered by imaging spectrometer can be expressed as the weighted linear combination form of end member, and its weighting coefficient is also the abundance of corresponding end member.The early stage high spectrum image solution mixing method based on linear mixed model, as N-FINDR(M. E. Winter, " N-FINDR:An algorithm for fast autonomous spectral end-member determination in hyperspectral data, " in Proc. SPIE Conf. Imag. Spectrom., Pasadena, CA, Oct. 1999, pp. 266 – 275.), convex constituent analysis (VCA) (J. Nascimento, J. Bioucas-Dias, Vertex component analysis:A fast algorithm to unmix hyperspectral data, IEEE Trans. Geosci. Remote Sens., vol. 43, no. 4, pp. 898-910, Apr. 2005.), independent component analysis (ICA) (J. Bayliss, J. Gualtieri, and R. Cromp, " Analyzing hyperspectral data with independent component analysis, " in Proc. SPIE, 1997, vol. 3240, pp. 133-143.) etc., pure pixel is there is in usual supposition in high spectrum image.But because the spatial resolution of current imaging spectrometer is relatively low, and the existence of mixing phenomena under various yardstick, the hypothesis that pure pixel exists is usually and (M. Iordache of being false, J.M. Bioucas-Dias, and A. Plaza, Collaborative Sparse Regression for Hyperspectral Unmixing, IEEE Trans. Geosci. Remote Sens., vol.52, no.1, pp.341-354, Jan. 2014.).In order to solve the non-existent problem of pure pixel, some are suggested in succession based on constrained non-negative matrix decomposition (CNMF) with based on the high spectrum image solution mixing method of constraint sparse regression (CSR).The invention belongs to the high spectrum image solution mixing method based on constraint sparse regression.
A given panel height spectrum picture
y, first the solution mixing method based on CNMF needs to estimate end member number, then under abundance nonnegativity restrictions (ANC) and abundance and a constraint (ASC), and will
ybe decomposed into end member matrix
aand abundance matrix
xproduct.
HU algorithm based on NMF need know the quantity of end member in advance, and needs to solve end member matrix and abundance matrix simultaneously.The another kind of HU algorithm based on LMM corresponding with constraint NMF algorithm mixes algorithm based on the solution of sparse regression.Such algorithm is by compressed sensing (CS), some thoughts in sparse regression (SR), utilize the spectrum storehouse be made up of the spectrum of various material being collected in ground or laboratory, structure spectrum storehouse matrix, and mixed pixel observation is modeled as the linear combination problem of composing in storehouse, the function class of this spectrum storehouse matrix is like the dictionary in CS and SR, then HU problem is converted into constraint sparse regression (CSR) problem (D. Iordache, J. Bioucas-Dias, and A. Plaza. Sparse unmixing of hyperspectral data. IEEE Transactions on Geoscience and Remote Sensing, 49 (6): 2014 – 2039, 2011.).This estimates and this step of Endmember extraction based on the necessary end member number of method of NMF before avoiding.
Still use
a=[
a 1.,
a m ] (
) represent by
mindividually to comprise
lthe spectrum storehouse matrix of the spectrum signal formation of individual spectral coverage, the data matrix of given high spectrum image
y=[y
1.., y
n ] (
), can obtain according to LMM
Y=
AX+
E
Wherein
x=[
x 1..,
x n ] (
) be abundance matrix,
e=[
e 1...,
e n ] be noise matrix.Right
ysolution mixed mean the signal subset that optimum found out by needs from (may be very large) spectrum storehouse, this signal subspace energy collecting represents best
yin each mixed pixel
y i .Ideally, if
yin do not comprise spectrum signal
a i , then
x?
irow should be zero entirely, namely
xin only have
ain row just non-zero corresponding to active member (Active members); Generally speaking,
ycan suppose to be mixed by a small amount of end member signal, and
ain the general signal that comprises more, thus
xin should have be entirely much 0 row.For this reason, document (M. Iordache,, J.M. Bioucas-Dias, and A. Plaza, Collaborative Sparse Regression for Hyperspectral Unmixing, IEEE Trans. Geosci. Remote Sens., vol.52, no.1, pp.341-354, Jan. 2014.) propose a kind of based on
x's
l 2,1mixing norm
(
x(
k:) represent
x?
kjoint sparse OK) returns separates mixed (CLSUnSAL) algorithm, and the optimization problem of its correspondence is
,
Wherein
l r+
(
x) be indicative function, equal 0 when representing the whole non-negative of X, otherwise be infinity.Experiment shows, CLSUnSAL is a kind of more outstanding solution mixing method.
The present invention is on the basis of CLSUnSAL, and propose a kind of high spectrum image solution mixing method returned based on weighting joint sparse, its model is
,
Wherein
for getting the nonnegative constant of 1 or 2,
(wherein
x(
i:) and representing matrix
x?
ioK), namely
xin all row vectors
-norm sum,
1 p with
1 n representation dimension is respectively
pwith
nelement be entirely 1 column vector; And utilizing variable division and alternating direction multiplier method (ADMM) to solve this model, experimental result shows that carried model and algorithm is all better than CLSUnSAL in the mixed precision of solution and operational efficiency.
Summary of the invention
1, object: the object of this invention is to provide a kind of high spectrum image solution mixing method returned based on weighting joint sparse, the method is by introducing weighting
mixing norm, and based on variable division and the heavy weighting algorithm of alternating direction multiplier method structure iteration, quick, the exact solution that realize EO-1 hyperion solution image are mixed.
2, technical scheme: the present invention is achieved by the following technical solutions.
Based on the high spectrum image solution mixing method that weighting joint sparse returns, it is characterized in that, comprise the following steps:
Step one, read data: the remote sensing images collected data from imaging spectrometer, obtain data cube, hyperspectral image data is removed by the wave band of water vapor absorption and the lower wave band of signal to noise ratio (S/N ratio), hyperspectral image data is arranged by pixel, obtains original high spectrum image matrix
y; If high spectrum image has
lindividual wave band, total
nindividual pixel,
y 0=[
y 1,
y 2...,
y n ], wherein
y i (
i=1..
n) be high spectrum image
ithe spectrum column vector of individual pixel is one
ldimensional vector; Read existing high-spectral data database data, select pure substance spectra data construct library of spectra matrix in library of spectra
a; If the quantity of pure substance spectra is
p,
a 0=[
a 1,
a 2...,
a p ],
a j (
j=1..
p) be in library of spectra
jthe spectrum column vector of individual pure material is one equally
ldimensional vector;
Step 2, high spectrum image solution is mixed: establish in high-spectral data storehouse containing enough abundant pure substance spectra data, the end member so contained in high spectrum image only accounts for a few part, that is in high spectrum image, the curve of spectrum of each pixel is made up of the curve of spectrum linear combination of the pure material of minority in high-spectral data storehouse, and this embodies the openness expression of end member; End member is sparse in high-spectral data storehouse, and the abundance corresponding to them also has openness, and namely abundance matrix is sparse, openness as regularization term using abundance, realistic physical significance; Here, adopt linear mixed model, by the curve of spectrum that a detection obtains, be decomposed into the form of pure substance spectra linear combination in library of spectra, its coefficient is corresponding abundance; If abundance matrix is
x, its size is
px
n, its all elements meets non-negativity constraint, and linear mixed model is expressed as
y 0=
a 0 x; Order
w 1,
w 2...,
w p for
pindividual positive weighted number,
wfor with
w 1,
w 2...,
w p for the diagonal matrix of diagonal entry,
x(
i:) and representing matrix
x?
ithe row vector that row is formed,
| represent
x(
i:)
-norm; The optimal model mixed based on the high spectrum image solution of weighting joint sparse recurrence is:
,
Wherein
,
represent
a 0 x-
y 0frobenius norm square, namely matrix (
a 0 x-
y 0)
t(
a 0 x-
y 0) mark,
for regularization parameter,
1 p with
1 n representation dimension is respectively
pwith
nelement be entirely 1 column vector; Order
,
(wherein
for the constant preset, be as usually taken as 15), above-mentioned optimal model being relaxed is
;
Wherein
l r+
(
x) indicative function: if Z=0, then
l r+
(
z)=0, otherwise
l r+
(
z)=
; Introduce variable
z, and requirement
wZ=
x, with season
m=
aW, then have:
Utilize the method that variable divides, introduce free variable V, and above-mentioned optimization problem is converted into following constraint type
The Augmented Lagrangian Functions of above formula is:
Wherein
dfor following constraint
v=
zcorresponding Lagrange multiplier,
for non-negative punishment parameter, represent and ask matrix trace; By alternating direction multiplier method (ADMM), obtain following iterative process:
Wherein
krepresent iterations;
With diag (
a 1,
a 2...,
a n ) represent from the upper left corner with
a 1,
a 2...,
a n for the diagonal matrix of the elements in a main diagonal,
i p represent
prank unit matrix,
x(
r:) and representing matrix
x?
rthe row vector that row is formed; To any real number
xand constant
, definition soft-threshold contracting function
, to arbitrary real matrix
b,
represent matrix
bcarry out shrinking computing by element, namely
, be also
?
irow
jthe element of row is
, wherein
b i,j representing matrix
b?
irow
jthe element of row; Definition matrix soft-threshold contracting function
, wherein
pfor matrix
bline number,
b(
i:) (
i=1...
p) represent
b?
ioK, vect-soft is vectorial soft-threshold contracting function, given vector
b,
vect-
softbe defined as
; The concrete steps mixed based on the high spectrum image solution of weighting joint sparse recurrence are as follows:
(1) suitable parameter is selected
,
,
,
and
, order
,
;
(2) initialization abundance matrix
x 0,
m=
a, and make
w 0=
i p ,
z 0=
x 0,
d 0=0,
k=0;
(3) following steps are repeated:
(a)
V k+1
=
;
(b) Z
k+1=max{(
M T M+
I P )
-1[
M T Y+
(
V k+1
-
D k )], 0};
(c)
X k+1
=
W k Z k+1
;
(d)
W k+1
=
;
(e)
D k+1
=
D k +(
Z k+1
-
V k+1
);
(f)
M←
AW k+1
;
(g)
μ←max(
ρμ, μ max ),
k←
k+1;
(4) until
or k>
n iter;
(5) make
x=
x k , export
x;
Wherein
εfor the limits of error given in advance,
n iterfor greatest iteration step number given in advance;
Step 3, obtains abundance figure and real end member: in acquisition abundance matrix
xafter, suitable threshold value is set,
xin be less than threshold value element be set to zero, the element being not less than threshold value does not process; Find out
xin containing the line label of nonzero element,
xline label and library of spectra matrix
a 0row label corresponding; Take out library of spectra matrix
a 0in the row of corresponding row label, be real endmember spectra; Take out abundance matrix
xrow corresponding to middle line label, is the abundance figure corresponding to end member; The large percentage shared by region representation end member that color is brighter in abundance figure, the ratio corresponding compared with dark region is little, and color is that the region representation of black is not containing this end member.
accompanying drawing explanation:
Fig. 1 is the method for the invention process flow diagram;
Fig. 2 is the gray level image that the 30th wave band of high spectrum image AVIRIS Cuprite is corresponding;
Fig. 3 maps drawing for estimate according to the inventive method to the abundance that high spectrum image AVIRIS Cuprite carries out composing the mixed rear four kinds of mineral obtained of solution, and the abundance of subgraph (a), (b), (c), the corresponding mineral Kaolin/Smect KLF508 of (d) difference, Alunite GDS84, Chalcedony CU91-6A, Jarosite GDS99 K maps drawing for estimate.
embodiment:
Below in conjunction with accompanying drawing, the present invention is described in further detail.
What this example was selected waits that separating mixed high spectrum image is famous AVIRIS Cuprite, and this high spectrum image has 224 spectral coverages, uniform fold 0.2 ~ 2.4
μthe spectral range of m.Due to the cause be only absorbed by the water and signal to noise ratio (S/N ratio) is lower, separate before mixing, spectral coverage 1-2,105-115,150-170 and 223-224 will be removed, and only leave 188 spectral coverages altogether.Fig. 2 corresponding to the 30th wave band of high spectrum image AVIRIS Cuprite gives at gray level image, and image size is 250 x 191, adds up to number of pixels
n=47750.The corresponding length of each pixel removing the AVIRIS Cuprite after contaminated general section is the vector of 188, all pixels of the AVIRIS Cuprite after contaminated general section of removal is sequentially sequentially formed and waits to separate mixed data matrix
y 0all row, then
y 0size be 188 x 47750.
The spectrum storehouse that this example uses comes from United States Geological Survey (USGS) and is published on the spectrum storehouse splib06(in September, 2007 and can be loaded in http://speclab.cr.usgs.gov/spectral.lib06 down), spectrum signal in this spectrum storehouse comprises 224 general section equally, uniform fold 0.2 ~ 2.4
μthe spectral range of m.This example is by a subset of Shi Pu storehouse, the storehouse splib06 of use, and the spectrum signal number wherein comprised is
p=342, require that the angle of any two spectrum signals selected is greater than 2 degree.With to wait to separate mixed high spectrum image the same, contaminated spectral coverage 1-2,105-115,150-170 and 223-224 will be removed, and only leave 188 spectral coverages altogether.With
a 0represent the spectrum storehouse matrix that the spectrum that this example is selected is formed,
a 0corresponding one of each row eliminate the spectrum signal after contaminated spectral coverage, then
a 0size be 188 x 342.
This example is got
.
Concrete implementation step based on the high spectrum image solution mixing method of weighting joint sparse recurrence is as follows:
(1) read in data to wait to separate mixed data matrix
y 0and library of spectra matrix
a 0;
(2) implement spectrum solution to mix, comprise following steps:
(1) suitable parameter is selected
=0.005,
,
,
=10
μ,
=10
-3,
n iter=100, get
, and make
,
;
(2) initialization abundance matrix
x 0for full null matrix, and make
m=
a,
w 0=
i p ,
z 0=
x 0,
d 0=0,
k=0;
(3) following steps are repeated:
(a)
V k+1
=
;
(b) Z
k+1=max{(
M T M+
I P )
-1[
M T Y+
(
V k+1
-
D k )],0};
(c)
X k+1
=
W k Z k+1
;
(d)
W k+1
=
;
(e)
D k+1
=
D k +(
Z k+1
-
V k+1
);
(f)
M←
AW k+1
;
(g)
μ←max(
ρμ, μ max ),
k←
k+1;
(4) until
or k>
n iter;
(5) make
x=
x k , export
x;
(3), abundance figure and real end member is obtained: in acquisition abundance matrix
xafter, suitable threshold value is set,
xin be less than threshold value element be set to zero, the element being not less than threshold value does not process; Find out
xin containing the line label of nonzero element,
xline label and library of spectra matrix
a 0row label corresponding; Take out library of spectra matrix
a 0in the row of corresponding row label, be real endmember spectra; Take out abundance matrix
xrow corresponding to middle line label, is the abundance figure corresponding to end member; The large percentage shared by region representation end member that color is brighter in abundance figure, the ratio corresponding compared with dark region is little, and color is that the region representation of black is not containing this end member.
Four kinds of mineral such as Kaolin/Smect KLF508, Alunite GDS84, Chalcedony CU91-6A, Jarosite GDS99 K are existed
a 0the row sequence number of planting corresponds to
xline number, take out corresponding row data, and to be translated into size be 250 x 191 pseudo color images, result is as shown in Fig. 3 (a)-Fig. 3 (d).Fig. 3 (a)-Fig. 3 (d) clearly demonstrate that four kinds of mineral distribution situation everywhere and enrich degree in high spectrum image AVIRIS Cuprite.The reference mapping graph (referring to http://speclab.cr.usgs.gov/cuprite95.tgif.2.2um map.gif) that this result and USGS provide is basically identical, shows that the high spectrum image solution mixing method returned based on weighting joint sparse of the present invention quite accurately can mix the distribution of the mineral estimated in given scenario according to high spectrum image solution.
The invention provides a kind of high spectrum image solution mixing method returned based on weighting joint sparse; the above is only the preferred embodiment of the present invention; should be understood that; for those skilled in the art; under the premise without departing from the principles of the invention; can also make some improvements and modifications, these improvements and modifications also should be considered as protection scope of the present invention.The all available prior art of each ingredient not clear and definite in the present embodiment is realized.