第39卷2 0 19,第4期 年4月
谱学与光谱分析
Spectroscopy and Spectral Analysis
V o '39,N o.4,pp1118-1127
A pril,2019
一种稀疏约束的图正则化非负矩阵光谱解混方法
甘玉泉1!2,刘伟华\冯向朋\于涛\胡炳樑\汶德胜
中国科学院西安光学精密机械研究所,陕西西安710119 中国科学院大学,北京100049
由于受到高光谱遥感图像传感器平台的限制,图像的空间分辨率受到一定影响,这导致高光谱遥
感图像的像元通常是多种地物的混合,也叫做混合像元。混合像元的存在制约了高光谱遥感图像的准确分 析和应用领域。采用高光谱解混技术可将混合像元分解为纯净的物质光谱(E n d m em b er 端元)和每种物质 光谱所对应的混合比例(Abundance ,丰度&为获取更多更精细的光谱提供了可能。这对高精度的地物分类 识别、目标检测和定量遥感分析等研究领域具有重要的意义。因此,解混技术成为高光谱遥感图像领域的一 个研究热点#基于线性光谱混合模型(linearspectralm ixingm odel ,LM M ),提出了一种端元丰度联合稀疏约 束的图正贝^化非负矩阵分解(endmemberandabundancesparseconstrainedgrapW regularizednonnegativema-
trix factorization ,EAGLNM F )算法# 该算法通过研究基于非负矩阵分解(nonnegative matrix factorization , N M F )的方法,结合图正则化理论来考虑高光谱数据内部的几何结构,将端元光谱稀疏约束和丰度稀疏约束
应用于其中,从而能够对高光谱数据的内部流形结构进行更为有效的表达。首先,构造了 E A G L N M F 算法 的损失函数,采用V C A -F C L S 方法进行初始化,然后,设定相关参数,包括图正则化权重矩阵参数、端元光 谱稀疏约束因子和丰度矩阵稀疏约束因子,最后,通过推导得到了端元矩阵与丰度矩阵的迭代公式,并且设 置了迭代停止条件。该方法不受图像中是否有纯像元的限制。实际上,在现行高光
谱遥感传感器平台情况 下,高光谱遥感图像中几乎不存在纯像元,因此,E A G L N M F 方法为高光谱遥感图像的实际应用提供了一 种思路。采用合成的高光谱数据,构造了 4个实验来分析该方法的可行性和有效性,实验将该算法与VCA-
FCLS ,标准N M F 及G L N M F 等经典的解混算法进行比较,通过光谱角距离(spectral angle distance ,SAD)
和丰度角距离(abundance angle distance ,A A D )这两个度量标准来进行比较#实验1是总体分析实验#在固 定的信噪比和固定端元数目的情况下,用以上三种经典方法与E A G L N M F 同时进行解混#实验2是S N R 影 响分析实验。在固定端元数目和不同信噪比的情况下,用这四种方法进行解混。实验3端元数目分析实验。 在固定信噪比和不同端元数目的情况下,用四种方法进行解混,并且将结果进行对比。实验结果发现提出的
E A G L N M
F 方法在提取端元精度和估计丰度精度上都更为准确#同时,实验4是稀疏因子分析实验#对端
元稀疏约束和丰度稀疏约束之间的影响因子进行分析,实验结果表明引人的端元稀疏约束对于解混结果也 具有较好的影响,并且端元稀疏约束和丰度稀疏约束之间的影响因子也对解混结果具有一定影响。最
后,将 该算法应用于A V IR IS 所采集的真实高光谱图像数据,将其解混结果与美国地质勘探局光谱库中光谱进行 匹配对比,其提取的平均端元精度相比于其他三种方法要稍好。
关键词高光谱图像'图正则化'稀疏约束'非负矩阵分解'光谱解混中图分类号:TP751. 1
文献标识码:A
DOI : 10. 3964/j. issn. 1000-0593(2019)04-1118-10
然而,由于受到图像传感器平台的限制,其空间分辨率受到
高光谱遥感影像以其纳米级的光谱分辨率所获取的光谱一定影响,同时,由于自然界地物的复杂多样性,所获取的
收稿日期:2018-02-09,修订日期:2018-06-28
基金项目:国家重点研发计划项目(2017YFC1403700),国家自然科学基金项目(61501456),中国科学院西部之光青年项目(XAB2016B20)
资助
作者筒介:甘玉泉,女,1986年生,中国科学院西安光学精密机械研究所副研究员
e-mail : ganyuquan@opt. a. cn
引言
信息而广泛应用于伪装目标探测[12]、环境监测[3]、天文探
测[4]、物质分类识别[5]、精细农业[6]及岩矿识别[7]等领域。
第4期光谱学与光谱分析1119
高光谱遥感图像中单像元得到的光谱不一定只是一种物质的 光谱,可能是几种不同物质光谱的组合。这样的像元被称为混合像元。混合像元的存在制约了高光谱遥感影像的准确分析和应用。为了扩大遥感图像的应用领域和提高其应用精度,必须解决混合像元问题。光谱解混技术就是用来解决混合像元问题的一项技术。它将高光谱图像的混合像元分解为端元和丰度的组合,为更精细的光谱应用提供了可能。因此,光谱解混技术是实现高光谱遥感技术定量化研究的重要条件。
光谱解混技术需要首先建立光谱的混合模型。目前,根 据物质的混合状态和物理作用,光谱混合可以分为线性混合模型和非线性混合两种模型线性混合模型仅考虑到达传感器的光子与某一地物发生作用,忽
略了物质相互之间的影响。然而,当地物分布尺度较小时,光子与多种物质发生作用,导致非线性混合。非线性混合模型不仅考虑光子与地物之间的作用,还要考虑光子在不同物质之间的散射作用[89]。在非线性混合模型中,可分为基于辐射度理论的非线性模型和基于计算理论的混合模型。比较典型的采用基于辐射度理论的非线性光谱混合模型有H ap ke模型[1+]和几何光学模型)1]#
基于辐射度理论的模型通常针对特定地物类型,需要大 量先验知识,从混合光谱产生根源上探究光谱解混工作。然 而,很多先验知识很难获取,需要对特定地物进行大量研究。所以,近年来,从计算方法的角度来进行非线性解混的研究也相继出现,主要有基于神经网络的方法、基于核函数的方法、基于流形学习的方法等。其代表方法有多层感知器MLP方法(multilayer p e r Ceptkm)[12],神经网络A R TM A P 方法(adaptive resonance th eoryM A P)[13]等# 总体来说,相 比于非线性混合模型,线性光谱混合模型具有建模简单、物 理意义明确等优点,也是目前最常用的光谱解混模型。
光谱解混技术通常包括三个研究方向,分别是获取端元数目,提取端元光谱及估计端元丰度。目前获取端元数目方 向的研究较少,较为著名的算法有Hysime[14]及VD[15]。大 量的算法集中于研究提取端元光谱及其相应的丰度。
基于线性混合模型的端元提取方法可以分为两类。一类 是基于几何学的方法,另一类是基于统计学的方法。基于几 何学的方法认为高光谱影像的所有像素点都集中于一个单形 体内,其顶点则认为是端元。
具有代表性的方法有纯像元指数(pixel purity index,PPI)[1=17],N-FINDR[18],顶点成分分 析(vertex component analysis,V C A)[19],正交子空间投影方 法(orthogonal subspace projection,OSP)[2°],及体积增长方法(simplex growing algorithm,SG A)[21]等# 基于几何学的方 法可理解为直接从高光谱图像中提取端元光谱,因此需要假 设图像中存在纯像元。但是在实际高光谱遥感图像中,纯像 元几乎不存在,所以其解混精度相对较低。基于统计学的方法则不需要假设图像中存在纯像元,充分利用数据的统计特性来求解端元光谱。基于非负矩阵分解(N M F)[2]的解混方 法也是统计学方法中的一种较为有代表性的方法。
非负矩阵分解是由L e e和Sueng在1999年提出的将一个非负矩阵分解为两个低秩非负矩阵乘积的矩阵分解方法[22]。而对于基于线性模型的高光谱图像而言,高光谱图像 本身,其端元光谱矩阵及丰度矩阵都可以看作是非负矩阵,因而可以将非负矩阵分解用于求解光谱解混问题。非负矩阵 分解可保证非负性且无需指定迭代步长,同时,非负矩阵分解不需要高光谱图像中纯像元的存在。因此,基于非负矩阵分解的方法在高光谱解混领域具有其优越性。
目前,一些学者在基于N M F的解混技术上已经奠定了研究基础。2〇〇7年,M6o和Q i提出了最小体积约束的NMF(minimum volume constrained N M F,MVC-N M F)解混 算法[23]。由于N M F的目标函数非凸,通过加人体积约束,将最小二乘分析和凸面几何结合起来,取得了较好的解混结果。但是,该方法中体积约束采用行列式来计算,导致梯度 计算较为复杂,并且在每次迭代过程中将负值强制置〇来保
证非负性,这对收敛也产生一定影响#2011年,c a等提出 了采用图正则化(graph regularized N M F,GNM F)的概念来 表示数据[4] #随后,R a a b i等于2011年将G N M F方法应用 于高光谱解混[25],其结果相比于传统的N M F方法精确程度 更高,虽然该方法不仅考虑到了高光谱数据的欧式距离内部结构,同时考虑了其内部黎曼几何结构,但是忽略了丰度的稀疏特性#2013 年,Lu 等提出了Graph Regularized L1/2 NM F(GLNM F)方法来进行光谱解混,它在图正则化基础上加人丰度稀疏约束,并且取得了较好的结果[26],但是其算法 的解混精度仍有提高的空间。
本文针对上述研究中存在的问题,提出了在G L N M F算 法的基础上, 对端元光谱矩阵也加人了稀疏约束的方法,称 为端元丰度稀疏约束的图正则化非负矩阵分解(EAGLN M F)算法。通过合成高光谱数据和真实高光谱数据验证,本文方 法取得了较好的实验结果。
1理论介绍
1.1线性混合模型
在线性混合模型中,高光谱影像的像元为各端元按照一定比例的线性组合,各端元在像元中所占比例称为丰度,丰 度需要满足和为一约束和非负约束。假设第z个像元用
(J1,2,…,N)表示,则可以表达为如式(1)所示
= #a j,1 %n%%)
1-1
其中,9为端元光谱,j,i&S0XN为对应端元的丰度值,0为端元数目,L为光谱通道数,iV为像元数目,'为像 元叫中的噪声。其非负约束(abundance non-negativity con­straint,A N C)) 与和为—'约束(abundancesum-to-one con-straint,ASC)和表示如式%)
#- 1,' 0 %)
1-1
1.2非负矩阵分解
给定一个非负矩阵V和一个正整数r C m i n a,V),通过非负矩阵分解N M F可以到两个低秩非负矩阵W 和丑&K0X N,使其满足
,(WH%)
1120光谱学与光谱分析第39卷
定义一个目标函数来描述分解后乘积对原矩阵的逼近程度,它通过最小化欧氏距离的目标函数来表示。非负矩阵分解的 目标函数表示如式%)
J(W,/&=2||,WH) | #(,,, &(W H)2
其中,I I$II_f表示Frobenius范数,简称F范数。为了控制 步长和矩阵的非负性。在文献[22]中采用乘性迭代规则,表 示如式(5)和式%)
W *W. " (YHT)./(W H H T)%)
H *H. " (WT Y)./(WT W H)%)
式(5)和式(6)即为矩阵W和H的迭代公式。采用乘性 迭代规则,可以自动调整步长以进行迭代,对每个矩阵元素施以不同步长,保证了非负性,消除了参数选择带来的影响。然而,一般情况下,目标函数具有明显的非凸性,存在 大量局部极小值,因此会有不唯一解,这也是N M F方法所 存在的缺点。
1.3 GLNMF
N M F通过其非负约束可在欧氏空间中表示数据,但是 并没有考虑数据内部的几何特性,也无法识别数据空间上的结构特征。G N M F则克服了这个缺点。从几何学角度来看,数据通常是从一个嵌人高维空间的低维流形数据所采样得到 的。因此,G N M F方法通过构造最近邻域图来解译数据内在的几何信息[#]。实际应用中,G N M F方法在目标函数中引人权重矩阵来构造一个新的正则项。但是,G N M F方法没有 考虑到丰度矩阵的系数特性。G E N M F方法则在目标函数中引人了丰度矩阵约束,进一步提高了解混精度。其损失函数可以表示为
J0,1) = * ||2—01 ||_| %A ||1 ||1/2 %—*~#T r(H T)
%)其中,A用来调整丰度矩阵稀疏约束|| 1||1/2的影响权重,I I 1||1/2的定义如式%),T r( $ )表示矩阵的迹,并且,)=0—W,-I=#W g,#是图正则化影响因子。
I
;N
II1II1/2 #E S I/2%)
0—1n—1
假设数据2=[X1,心,…,XN]&Kf N为给定数据,那 么,对每个^(=1,…,N)可代表L维空间中的一个数据点,这N个数据点可作为顶点构造最近邻域图。邻域图的权 重矩阵用W表示。W的定义有很多方法,其中0—1权重,热核权重,点积权重为三种常用的权重定义方法[#]。考虑到 高光谱数据的凸面几何特性,即高光谱数据构成的凸面几何体的顶点为端元光谱,因此,选择热核加权方法来定义权重矩阵W更适合。如果点^是点^的最近邻域之一,那么权 重可表示如式%)
其中,C7用来控制两点之间的相似程度。
2EAGLNMF算法描述
2.1损失函数定义
通常情况下,由于真实的端元光谱矩阵本身并不具备稀疏性,因此在每次迭代过程中,端元光谱矩阵是稀疏的,而且多个稀疏矩阵的乘积形成一个非稀疏矩阵。有关此理论的证明仍是一个开放性问题[27_28]。基于上述理论,本文在GL-N M F方法的基础上引人端元光谱的稀疏约束。无论端元光谱矩阵本身是否为稀疏矩阵,引人稀疏约束后对端元光谱矩阵进行稀疏分解,在降低噪声影响的同时,可进一步提高端元光谱的提取精度。本文将算法的损失函数定义为如式(0) J(0,1) —1||2 — 01 ||f %—2#T r (D1T) %
a I I A I I1/2% 沒 I I1 I I1/2(10)其中,A为估计出的端元,1为分解出的丰度矩阵,参数#
用来控制光谱数据的平滑程度。《 = «。0^,t是优化过程中的迭代次数,ao和r是用来调整稀疏程度的常数,并且设置沒=a,为了加强丰度稀疏约束对解混结果的影响,通常将0设置成一个大于1的数。同时,丰度需要满足A N C和ASC
这两个约束条件((1)1n>0,(2)1>1=1N,其中,1n表示维度为;X N的丰度矩阵,1;表示一个值全为1的;维的列向量,1n表示一个值全为1的N维的列向量。
2.2优化算法
根据f范数的性质,有以下等式成立,l|A||f =Tr (AAT),而且,矩阵具有以下性质,Tr(AB) = Tr(BA)
T r U) —T r(0T)(11)式(10)改写如式(2)
J U,1) —1扑[(2 —A S)(2 —A1)T] % 1#T r(1)1T) % a|A|1/2 %^|1|1/2—1T r{'X2T&2S T AT—
012T g 011T0T) %—1#Tr(1L1T ) % a|A|i,/2 %
沒 II1II1/2—1r(22T) & r(2s T AT) %
1T r(A11T A T) % 2T r(S L1T) % a I A |1/2%
&I I1 I I1/2(12)采用拉格朗日乘子法和K K T条件,假设料和%分别
为常数〜'0和^'0的拉格朗日乘子,并且平=[料],®=[你],则拉格朗日L表示如式(3)
C —1T r C X2T) & TrC2ST AT) % 1-T r(A S1T AT) %
2T r(S L1T)% a|A |1/2% & I1 I1/2% T r()A T) %
T r(*1T) (13)对式(13)的C分别求A和1的偏导数
+C =—2S T+1T A+1u  1 %)(14) d+C=—2A T g A1A T%#L1 % 1[S  1 %*(15)
根据K K T条件有,料心=0和你5# = 0,结合式(14)和式(15)则有
& (21T)a#%1T A(%1〇A —0
(16)
第#期光谱学与光谱分析1121
—(XAT)jkS jk% (ASAT)jkS jk% ju(LS') jks j k% *sik#+
(17)
用式(16)及式(17)可以得到迭代规则如式(18)和式(19)
O S1)%8)
9k *9k(18)
(SST A % 1%A-- )n
(ASAT%y S—~ % i3S)}k
将其改写后,可得到乘性迭代规则如式(20)和式O1)
SSTA%1%A-i
s*s xA7±i m%1)
ASAT g1S—I%fj DS
2.3算法说明
O)不同的初始值会得到不同的结果。通常会采用两种初始化方法,分别是随机初始化和VCA-F C L S方法来进行初始化。随机初始化方法采用随机选择0到1之间的值作为 A和S的初始值。但根据以往研究结果,随机初始化方法往往没有VC A-F C L S初始化方法准确程度高,因此,本文采用 VC A-FC L S方法对所提出的算法进行初始化。
%)要保证丰度矩阵满足A S C和A N C约束。根据式(20)和式(21)的迭代规则,S显然不满足这两个约束条件,因此,在对S进行迭代时,采用文献[26]的方法,将矩阵X 和A用矩阵X£和A£代替,其分别定义如式(22)和式(23)
X
Xt j O2)
+1j v
A
At =T O3)
L+i T」
其中,H为像元个数,P为端元个数,+为一个常数,用来调 节A S C的约束效果,使得丰度和趋近于1。实验结果表明,+ 的选择会对解混精度造成影响,随着+的增大,光谱解混的 结果更为准确,但这也导致了收敛速度变慢。为了平衡解混精度和收敛速度,本文实验选择+=20。
%)本文设置两种迭代停止规则。①当迭代次数达到所设置的最大次数时迭代停止。②当目标函数不再收敛时也迭代停止。设置一个差值£作为停止迭代阈值,实验中选择£= 10一#,当连续十次目标函数差值小于£时,则迭代停止。其 目标函数定义如式(24)所示,停止迭代规则如式(25)所示。
〇=11X — AS 1II(24)
I I O n e w,〇l d||— £(25)
O)如果在迭代过程中,A和S中出现0值或者负值,则加上一个极小值来保证迭代可行。
〇)本文的算法流程表示如下,£表示停止迭代阈值,%0和r是用来调整稀疏程度的常数,A是调整端元稀疏与丰度稀疏程度的因子,即有&=%,#是调整参数,P为端元数目,Tm a x为最大迭代次数。
算法流程:
①输人(数据:原始高光谱数据矩阵(X);
参数:£,%0,(9,r,#,_P,T m a x
输出:端元光谱矩阵A及丰度矩阵S
② 具体步骤:
SteP1分别用V C A和F C L S方法来初始化端元A和丰 度矩阵S,令t=1;
Step2构建.,3矩阵;
S0p3采用式(20)和(21)对A和S进行迭代,令=纟十1;
Step#满足停止迭代条件,则迭代终止,输出A和S。
其算法流程图如图1所示。
图1EAGLNMF算法流程图
Fig. 1Flowchart of EAGLNMF
3实验部分
3.1评价标准
采用合成数据和真实高光谱数据对提出方法进行评估,并将其与VCA-FCLS、标准N M F和GLNM F等三种方法进行比较。采用两个定量评估标准,分别是光谱角距离(SAD) 和丰度角距离(AAD),分别对提取的端元及其对应丰度的准 确性进行评价。S A D用来计算估计出的光谱与真实光谱之间 的光谱角,SA D值越小,则估计出的端元光谱与真实光谱的差别越小。A A D是用来评价估计出的丰度与真实丰度之间的差别,A A D值越小,说明估计丰度越接近其真实丰度。其 公式表示如式(26)和式(27)
A A D( #cos—1 (s n s i i) (7)
为了对所有估计出的端元及丰度进行评价,采用SAD 和A A D的均方根(root m ean square,RM S)来进行评价,用 式%.8) 和式%.9) 表示
.0丄
rm sS A D =\—#(SAD()2](28)
0 1=
1
1122光谱学与光谱分析第39卷
*saad #[** #(AAD、)2] (29)
(=1
3.2合成数据
采用文献[23]的方法进行模拟数据合成。模拟数据大小 为64X64,从U S G S光谱数据库[29]中选择一系列反射光谱来合成的。本文合成数据所选择的光谱数据有22#个谱段,覆盖波长范围为+ 38〜2. 5 (m,光谱分辨率为10nm。将整 个图像分成8X8的小块,每个小块内都用六种不同物质光谱进行填充,六种物质光谱均从U S G S光谱数据库中选择,选中的6种物质光谱曲线如图2所示。所产生的图像然后通过
一个空间低通滤波器来进行线性混合。本文采用9X9的滤波器来合成模拟数据。同时,将丰度大于8〇\的像素点的丰度值用1/^来代替,并且加人〇均值高斯噪声来模拟可能出现的误差和传感器噪声等。算法分析通过#个实验来完成。实验1是总体分析实验,对几种算法进行总体比较分析。实验2是S N R影响分析实验,验证算法对信噪比(4gnal % nrnse mt%,SNR)的敏感程度。实验3是端元数目分析实验,验证端元数目对几种方法的影响程度。实验是系数因子分析实验,分析了稀疏系数《,/?之间影响因子'对解混结果的影响。同时,为保证结果准确性,所有的实验结果都是30次运 行结果的平均值。实验1、实验2和实验#都是用选中的6种物质来进行模拟数据的合成,除了实验3中根据所需要的 端元数目来进行模拟数据的合成。
3.3 实验环境
实验中,所有算法验证均采用M A T L A B2016b平台编 程,运行计算机配置如下:C P U型号为Intel(R)CoreCTM) i3-3110M CPU"2. # G H z,内存为 # G B,操作系统为Win7 64bit Service Pack1。
Fig. 2 Six material signatures selected from USGS
4结果与讨论
4.1合成数据
正则化与稀疏
实验1(总体分析实验)实验1的目的是对几种不同算法进行总体比较。在本实验中,加人了 SN R为20 dB的高斯 噪声。采用 VCA-FCLS,NM F,GLNMF 及 EAGLNMF 几种 方法来进行求解。其中,NM F,GLNM F及EAGLNM F均采 用V C A及FCLS的结果作为初始迭代值。终止迭代规则为式(25),当连续10次迭代结果满足式(25)或者当迭代次数达到设定的最大次数了…61时,则终止迭代。E A G L N M F的参 数设置为 e J10—4,的=0. 1,『=25,'=2,"=0. 1O=6, Tm a x=  3 000。为了保证A S C约束,采用式(22)和式(23)中方 法来实现。图3(a)和(b)分别给出几种方法下rm^AD和 ?*M ad及其标准偏差。在图3(a)中,VCA-FCLS,N M F,GL-NMF和 EAGLNMF方法的 rm^SAD分别为 0. 164 0, 0. 193 0, 0.084 0和0.076 7,其对应的标准偏差分别为0.032 2,
0.057 4, 0. 023 1 和 0. 015 9。从图 3(a)中可以看出,VCA-
F C L S比标准N M F方法的rm sSAD及其标准偏差的值都更小,说明VCA-F C L S比N M F的端元提取精度更高,这是由于标 准N M F方法仅仅从非负矩阵分解的角度进行解混,没有考 虑到高光谱数据的几何特性。
G L N M F和E A G L N M F方法 的端元提取效果比VC A-F C L S和N M F方法要好一些,同时,从数据对比可以看出,E A G L N M F比G L N M F的端元提 取结果稍好,这是因为E A G L N M F方法对端元增加了稀疏
约束,在迭代过程中降低了噪声的影响。在图3(b)中:乂匸八- FCLS,N M F,GLNMF和 EAGLNMF方法的 rm_s AAD分别为0.384 1, 0. 545 6, 0. 291 4和0. 275 3,其对应的标准偏差分 别为 0. 077 6, 0. 080 7, 0. 064 1 和 0. 056 0。从图 3(b)中可 以看出,其结果与rm&A D具有一致性。总体来说,E A G L-N M F方法的解混效果要优于其他几种方法。
VCA-FCLS NMF GLNMF EAGLNMF
图3几种方法的(a)rmssA D和(b)rms AAD的对比
Fig. 3 Comparison of rmsS A D and rmsA A D for VCA-FCLS, NMF,GLNMF and EAGLNMF
实验2(SN R影响分析实验):本实验目的是评估噪声对算法的影响程度。实验选择相同的端元数目和不同程度信噪比的情况下,来对算法进行分析。此实验中,端元数目
5,信噪比分别为15 , 20, 25, 30, 35, 40和45d B的高斯噪 声加人了合成数据。其中,合成数据的生成方法与实验1中相同,参数设置也与实验1中相同。从图4(a)中可以看出

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。