基于光谱信息散度与光谱角匹配的高光谱解混算法

基于光谱信息散度与光谱角匹配的高光谱解混算法
刘万军;杨秀红;曲海成;孟煜
【摘要】针对采用线性逆卷积(LD)算法进行端元初选过程中,端元子集中存在相似端元光谱,影响解混精度的问题,提出了一种基于光谱信息散度(SID)与光谱角匹配(SAM)算法的端元子集优选光谱解混算法.通过在端元进行二次选择时,采用以光谱信息散度和光谱角(SID-SA)混合法准则作为最相似端元选择的判据,去除相似端元,降低相似端元对解混精度的影响.实验结果表明,基于SID与SAM的高光谱解混算法将重构影像的均方根误差(RMSE)降低到0.0104,该方法比传统方法提高了端元的选择精度,减少了丰度估计误差,误差分布更加均匀.
【期刊名称】《计算机应用》
【年(卷),期】2015(035)003
【总页数】5页(P844-848)
【关键词】光谱解混;端元选择;去除端元;解混算法
【作者】刘万军;杨秀红;曲海成;孟煜
【作者单位】辽宁工程技术大学软件学院,辽宁葫芦岛125105;辽宁工程技术大学软件学院,辽宁葫芦岛125105;辽宁工程技术大学软件学院,辽宁葫芦岛125105;哈尔滨工业大学电子与信息工程学院,哈尔滨150006;辽宁工程技术大学软件学院,辽宁葫芦岛125105
【正文语种】中文
【中图分类】TP751.1
0 引言
随着高光谱遥感成像技术的发展，高光谱遥感技术在遥感领域已经成为一个快速发展的分支。

混合像元解混的问题［1］一直存在于高光谱遥感图像处理中。

高光谱解混是指从混合像元中提取不同的基本组成单元(一般称端元)，并求解出这些基本组成单元所占比例(一般称丰度)的过程［2］。

国内外很多知名学者为了解决高光
谱混合像元的问题，提出了多种高光谱解混的方法。

高光谱混合像元分解［3］一般包括两个方面，分别是端元提取和丰度估计。

当前端元提取方法有N-FINDR、顶点成分分析(Vertex Component Analysis，VCA)［4］、迭代误差分析(Iterative Error Analysis，IEA)［5］、纯像元指数(Pure Pixel Index，PPI)等。

传统的丰度估计算法主要有非限制性最小二乘法、非负限制性最小二乘法、和为一限制性最小二乘法和全限制性最小二乘法。

端元组中包含的端元个数一般要大于一个像元包含的端元个数。

在进行端元丰度估计的时候，若采用的端元恰好与像元包含同样的物质，那么得到的估计精度是最高的。

如果选择的端元数量不合理，就会增大高光谱解混的误差。

所以对高光谱图像中的某个像元或者某条混合光谱进行解混时，需要尽可能选择最精确的端元。

本文提出一种基于光谱信息散度与光谱角(Spectral Information Divergence-Spectral Angle，SID-SA)混合法的高光谱解混算法，对进行端元初选后出现的相似端元，采用SID-SA作为最相似端元选择的判据，提取出最准确的端元组，去除相似端元，实现高精度的混合像元分解。

1 光谱区分方法
假设两光谱曲线x和y的长度为n，并且x=(x1，x2，…，xN)，y=(y1，y2，…，yN)，鉴别两种光谱曲线的相似性通常采用光谱角匹配和光谱信息散度这两种光谱匹配方法。

1.1 光谱信息散度
光谱信息散度(Spectral Information Divergence，SID)法是一种基于信息论衡量两条光谱之间差异的波谱分类方法［6］。

x和y两条光谱的概率向量分别是
a=(a1，a2，…，aN)和b=(b1，b2，…，bN)，其中，由信息理论，可以得到x 和y的自信息为:
通过式(1)和式(2)，可得y关于x的相对熵:
同样可以得到x关于y的相对熵:
因而x和y的光谱信息散度为:
光谱信息散度法的优点是可以对两条光谱进行整体上的比较。

1.2 光谱角
光谱角匹配(Spectral Angle Mapping，SAM)法是根据目标光谱和测试光谱之间的夹角数值大小来判断两条光谱曲线的相似性。

如果两个光谱间的夹角数值越小，则它们的匹配程度越高［7］。

即:
其中:SA(x，y)是SAM光谱匹配模型，计算测试光谱和目标光谱的光谱夹角，方向夹角越小，两条光谱的相似程度就越高。

光谱角匹配法虽然能够比较光谱在形状上的相似程度，但是较难区分光谱在局部特征上的差异性。

1.3 SID-SA 混合法
光谱角匹配法是通过计算测试光谱与目标光谱这两个矢量间的广义夹角来确定两者
的相似性。

当两光谱间的夹角小于给定阈值时，则说明两者是相似光谱;当两光谱
间的夹角大于给定阈值时，则说明它们不是相似光谱。

光谱角匹配法的原理是把光谱作为矢量投影到N维空间上，这里的维数是指实验中选取的波段总数。

在N维空间中，各光谱曲线被认为是既有方向又有长度的矢量，各光谱之间形成的夹角就是光谱角，如图1所示。

光谱角匹配法仅能对两条光谱进行形状上的比较，不能
从物理上反映出反射能量值的差异。

而光谱信息散度法可以弥补光谱角匹配法的不足，它是一种考虑光谱概率分布的随机方法，该方法是从光谱曲线的形状出发计算各个信息点包含的信息熵，比较信息熵的数值大小判断两条不同曲线的相似程度。

在分析以上两种方法特点的基础上，将光谱角匹配法和光谱信息散度法结合在一起，从而能够更准确地判断出两条光谱的相似性。

SID-SA混合方法即:
图1 多维光谱角度
2 基于SID-SA的高光谱解混
首先在端元初次提取时先去除部分端元;再进行端元的二次提取迭代循环，去掉相
似端元组，剩下的端元为最佳端元;最后再利用最佳端元进行丰度的估计，最终实
现高光谱的混合像元分解。

算法流程如图2所示。

图2 基于SID-SA的高光谱解混算法流程
在进行端元初次选择过程中，根据信噪比(Signal-to-Noise Ratio，SNR)的不同
采用两种不同的规则进行端元的初次选择。

在低信噪比的情况下，使用全约束最小二乘法(Full Constraint Least Square，FCLS)［8］进行端元的初选。

在数据的信噪比较高的情况下，采用线性逆卷积(Linear Deconvolution，LD)规则［9］，
该规则以比较少的循环次数去掉比较多的无用端元。

当信噪比较高的时候，LD规则选取端元的精度也是比较高的，这样不但可以提高算法的运行速度，也可以获得比较高的精度。

采用全约束最小二乘法每个循环只能去除一个端元，效率没有LD
规则高，但是在数据的SNR较低时端元选取的能力较高，所以当SNR较低时用它进行端元的初选比较好。

根据文献［10］认为选择阈值40比较合适。

端元二次选择的目的是进一步从原始高光谱影像中找到纯像元，从而准确求解各纯物质在混合像元中所占的比例。

以SID-SA作为最相似的端元选择准则，能够选择出最相似的端元组，最相似端元组经加权处理得到新端元，用新端元代替得到的相似端元组，从而降低相似端元对解混精度的影响。

端元二次选择优选算法步骤如下:
步骤1 找出第一组相似端元。

将端元集合利用光谱匹配算法SID-SA进行第一次循环，找出最相似的一组端元，分别是端元A和端元B。

步骤2 找出与端元A相似的样本点。

从图像中随机提取出K个样本点，利用光谱角匹配法将这K个样本点分别和端元A作光谱角匹配。

步骤3 找出与端元B相似的样本点。

即将这K个样本点分别和端元B作光谱角匹配。

步骤4 计算出样本点中分别和端元A、B相似的比例。

即计算出有m%的样本点和端元A相似，有n% 的样本点和端元B相似。

步骤5 得到替换端元A的新端元A1。

利用期望最大化(Expectation Maximization，EM)算法［11］，将 K 个样本点中与端元A相似的样本点作它们的数学期望，形成一个新端元A1，将其替换端元A。

步骤6 得到替换端元B的新端元B1。

利用EM算法，将K个样本点中与端元B 相似的样本点作它们的数学期望，形成一个新端元B1，将其替换端元B。

步骤7 得到替换相似端元组的新端元C。

将新的端元A1和新的端元B1按m%和n%的比例进行加权得到新的端元C。

步骤8 形成新的端元集合E。

将端元集合中新的端元A1和新的端元B1去除，并将得到的新端元C加入到剩余的端元集合中。

步骤9 形成新的端元集合F。

将步骤8中形成的端元集合E，重复步骤1，利用光谱匹配算法SID-SA找出最相似的两个端元，分别是端元G和端元H。

将端元G 和端元H按m%和n%进行加权得到新的端元J。

再将端元集合中的端元G和端元H去除，并将得到的新端元J加入到剩余的端元集合中。

步骤10 形成新的端元集合L。

将步骤9中形成的端元集合F，重复步骤1，利用光谱匹配算法SID-SA找出最相似的两个端元，分别是端元R和端元T。

将端元R 和端元T按m%和n%进行加权得到新的端元Q。

再将端元集合中的端元R和端元T去除，并将得到的新端元Q加入到剩余的端元集合中。

重复此步骤，直至选出均方根误差最小的端元组，最终形成新的端元集合L。

选出优化后的端元组之后，再采用最小二乘法求解高光谱最终的丰度。

3 算法验证方法
3.1 均方根误差
均方根误差(Root Mean Square Error，RMSE)［12］是一个度量准则，用来评价两幅图像的相似性，均方根误差越小，说明二者越相近。

其中:假设端元丰度矩阵为S=，真实丰度矩阵为Z=和zij分别为矩阵S和Z的第i 行第j列元素。

3.2 丰度估计误差
其中:m是图像中混合光谱条数;n是端元组中端元个数;qij是第i条光谱第j个端元的计算丰度值;pij是第i条光谱第j个端元的假定丰度值;fa是丰度估计误差。

3.3 相对光谱识别概率和相对光谱识别熵
本文采用相对光谱识别概率(Relative Spectral Discriminatory ProBability，RSDPB)［13］和相对光谱识别熵(Relative Spectral Discriminatory
Entropy，RSDE)［14］来评价SAM、SID和SID-SA混合法这三种不同光谱匹配方法。

1)相对光谱识别概率。

假设t是一种特定地物类型，{sv}Vv=1是V个特征化光谱信号，则相对光谱识别
概率的公式如式(10)所示:
其中:v的取值范围是［1，V］。

式中m(t，sv)是任一种方法在某一特定地物t中
第v个样本值。

RSDPB值越小，其光谱识别能力越强;反之，光谱识别能力越弱。

2)相对光谱识别熵。

相对光谱识别熵的公式如式(11)所示:
相对光谱识别概率的函数是相对光谱识别熵，当概率值的变化范围是［0，1］时，相对光谱识别熵函数先增大后减小，RSDPB的最大值是0.35。

随着RSDE值的减小，地物光谱的识别能力也随之增强。

4 相关实验及结果分析
4.1 实验数据
实验数据1 数据集为高光谱研究中经常采用的Pavia数据，该研究区域位于意大
利北部的帕维亚，该地区地形复杂，主要有灌木、水体、森林、农田、草地和城镇等地物类型。

实验数据2 数据集为高光谱解混研究中经常采用的美国Nevada州南部沙漠地区
的Cuprite数据。

该数据中共有224个波段，波长范围为0.37μm到2.48μm，
光谱分辨率是10 nm。

4.2 相对光谱识别熵对比
根据式(11)，利用实验数据一作为研究区域，计算出SAM、SID、SID-SA这三种
光谱识别参量的相对光谱识别熵。

SAM的相对光谱识别熵是2，SID的相对光谱识别熵是1.75，SID-SA的相对光谱识别熵是1.25。

从得出的数据可以看出，SAM、SID和SID-SA对地物的识别能力具有明显的差异。

从数据明显看出SID-SA的RSDE的数值是最小的，说明SID-SA对识别地物的能力是最强的。

因此，用SID-SA混合法找出相似端元组的能力强于SAM和SID。

4.3 相对光谱识别概率对比
根据式(10)，利用实验数据一作为研究区域，计算SAM、SID、SID-SA这3种光谱识别参量在不同地表覆盖类型下的相对光谱识别概率，如图3所示。

从图3可以看出，三种光谱识别参量SAM、SID和SID-SA在不同的地表覆盖类型下RSDPB数值变化顺序不一致，说明SAM、SID和SID-SA对不同覆盖物的光谱识别能力有所不同。

在草地、农田、城镇用地和森林区域，SID-SA的RSDPB 值比SAM和SID的RSDPB值小，说明SID-SA对草地、农田、城镇用地和森林的识别能力明显强于SAM和SID。

对于灌木和水体，SID-SA的RSDPB值是最大的，因此SID-SA识别灌木和水体的能力是最弱的。

所以，不难看出，针对本文算法SID-SA混合法寻找相似端元的能力强于SAM和SID。

图3 SAM、SID、SID-SA在不同地表覆盖条件下的RSDPB变化
4.4 均方根误差对比
RMSE是均方根误差，是一个度量准则，用来评价两幅图像的相似性，均方根误差越小，说明二者越相近。

利用实验数据二，根据本文优选算法分别提取出不同的端元数，得到的RMSE随着端元个数的增加的变化情况，如图4所示。

从图4中可以看出，利用本文算法提取出11个端元的均方根误差是最小的，所以得到8组相似的端元，根据8组相似的端元计算出其对应的SID-SA值，见表1。

图4 端元数目和RMSE值的关系
表1 8组相似端元元组号 SID-SA值相似端元组相似端元组号 SID-SA值相似相
似端端元组1 3.90E －05 (6，14) 5 2.49E －04 (14，15)2 1.95E －04 (8，16) 6 2.48E －04 (7，14)3 2.20E －04 (8，13) 7 2.69E －04 (4，6)4 2.35E －04 (9，14) 8 2.88E －04 (7，8)
4.5 优选前后端元及误差分布对比
在实验中，取不同的样本点得到的权值也就不同，实验中取了5组样本点K值，
从而得到了5组m和n的值，然后利用EM算法计算出这5组值的数学期望，从而得到了m和n的取值，m 取值是 38.06，n 的取值是 62.02，见表 2。

表2 加权百分比m和n的取值说明 K的取值 m的取值 n 的取值样本1 2000 38.65 61.35样本 2 2500 38.40 61.60样本 3 3000 37.33 62.67样本 4 3500 37.83 62.48样本 5 4000 38.11 61.98数学期望38.06 62.02
利用实验数据二进行优选前后的端元以及误差分布对比。

在端元初选时，利用LD 规则进行端元的一次提取处理，提取出19个端元的光谱特征曲线，如图5所示。

得到的误差分布，如图6所示。

对19条端元进行端元的二次提取处理后，得到
11个端元的光谱特征曲线，如图7所示。

得到误差分布，如图8所示。

图5 初次提取19种光谱特征曲线
图6 初次提取误差分布
图7 二次提取11种光谱特征曲线
图8 二次提取误差分布
对比图5和图7，可以看出端元的数目减少了，降低了相似端元对解混精度的影响。

对比图6和图8，可以看出优选后误差分布更加均匀，提高了端元选择的精度，降低了误差，分布更加均匀，减少了相似端元对解混精度的影响。

4.6 整体丰度估计误差对比
为了研究噪声对本文算法性能的影响，本文在混合光谱中加入了随机噪声。

加入噪
声的公式为:
其中:r'是加噪声之后的光谱;r是无噪声的光谱;u(0，1)是均值为0，方差为1的随机噪声;SNR是信噪比;X是噪声的假想发射率，假定为 0.5。

本文实验中将本文算法、LD算法以及FCLS算法进行对比，这3种算法的丰度估计误差随信噪比的变化情况，如图9所示。

图9 3种算法的丰度估计误差随信噪比变化曲线
从图9中可以看出，SNR在12～100时，本文算法丰度估计误差fa是最小的，并且在SNR是100时，本文算法的fa小于0.1;在SNR小于25时，LD算法的fa 大于FCLS算法。

因此在数据SNR较低时，比较适合采用FCLS算法;当数据SNR 较高时，比较适合采用LD算法。

5 结语
本文提出了一种基于SID-SA的高光谱解混算法，在进行端元初选后，采用以SID-SA作为最相似端元选择的判据进行端元的二次选择，去除相似端元，降低相似端元对解混精度的影响。

从实验的结果可以看出，提出的基于SID-SA的算法比传统算法提高了端元的选择精度，减少了丰度估计误差，分布也更加均匀;同时也减少了重构影像与原始影像之间的均方根误差(RMSE)。

该方法对高光谱遥感影像进行深度解译具有十分重要的意义。

参考文献:
［1］ MIANJIF A，ZHANG Y.Robust hyperspectral classification using relevance vectormachine ［J］.IEEE Transactionson Geoscience and Remote Sensing，2011，49(6):2100 －2112.
［2］ BIOUCAS-DIAS J M，PLAZA A，DOBIGEON N，et al.Hyperspectral unmixing overview:geometrical，statistical， and sparse regression-based
approaches ［J］.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2012，5(2):354 －379.
［3］ SOMERSB，ASNER G P，TITSL，et al.Endmember variability in spectral mixture analysis:a review ［J］.Remote Sensing of Environment，2011，115(7):1603 －1616.
［4］ KRAFFT C，DIDERHOSHAN M A，RECKNAGEL P，et al.Crisp and s
oft multivariate methods visualize individual cell nuclei in Raman images of liver tissue sections［J］.Vibrational Spectroscopy，2011，55(1):90 －100. ［5］ FANG Y，YAN H.Error analysis for remote nonlinear iterative learning control system with wireless channel noise ［J］.Journal of Shanghai University:English Edition，2011，15(3):7 －11.
［6］ JAROSIK N，BENNETTC L，DUNKLEY J，et al.Seven-year Wilkinson Microwave Anisotropy Probe(WMAP)observations:sky maps，systematic errors， and basic results ［J］.The Astrophysical Journal Supplement Series，2011，192(2):14.
［7］ CUBERO S，ALEIXOS N，MOLTO E，et al.Advances in machine vision applications for automatic inspection and quality evaluation of fruits and vegetables ［J］.Food and Bioprocess Technology，2011，4(4):487－504.
［8］ CHANG C I.An information-theoretic approach to spectral variability，similarity， and discrimination for hyperspectral image analysis［J］.IEEE Transactions on Information Theory，2000，45(5):1927－1932.
［9］ RAMSEY M S，CHRISTENSEN P R.Mineral abundance determination:quantitative deconvolution of thermal emission spectra
［J］.Journal of Geophysical Research:Solid Earth，1998，103(B1):577 －596.
［10］ XU Z，ZHAO H.Advanced spectral unmixing algorithm based on spectral information divergence ［J］.Journal of Beijing University of Aeronautics and Astronautics，2009，35(9):1091 －1094.(徐州，赵慧洁.基于光谱信息散度的光谱解混算法［J］.北京航空航天大学学报，2009，35(9):1091
－1094.)
［11］ MENG W，XIAO W，XIE L.An efficient EM algorithm for energy-based multisource localization in wireless sensor networks ［J］.IEEE Transactions on Instrumentation and Measurement，2011，60(3):1017－1027.
［12］ CHAI T，DRAXLER R R.Root Mean Square Error(RMSE)or Mean Absolute Error(MAE)？［J］.Geoscientific model Development Discussions，2014，7(1):1525 －1534.
［13］ NARESH KUMAR M，SESHASAIM V R，VARA PRASAD K S，et al.A new hybrid spectral similarity measure for discrimination among Vigna species［J］.International Journal of Remote Sensing，2011，32(14):4041 －4053.
［14］ XU M，WANG S，GU X，et al.Density estimation based on RSDE
for domain adaption ［J］.Control and Decision，2013，28(1):125－130.(许敏，王士同，顾鑫，等.基于RSDE的领域自适应概率密度估计方法［J］.控制与决策，2013，28(1):125 －130.)。