基于端元子集优选的高光谱解混算法研究

刘万军　杨秀红　曲海成,2　孟　煜　姜庆玲

1(辽宁工程技术大学软件学院　辽宁 葫芦岛 125105)2(哈尔滨工业大学电子与信息工程学院　黑龙江 哈尔滨 150006)3(铁岭师范高等专科学校　辽宁 铁岭 112000)



基于端元子集优选的高光谱解混算法研究

刘万军1杨秀红1曲海成1,2孟煜1姜庆玲3

1(辽宁工程技术大学软件学院辽宁 葫芦岛 125105)2(哈尔滨工业大学电子与信息工程学院黑龙江 哈尔滨 150006)3(铁岭师范高等专科学校辽宁 铁岭 112000)

摘要针对采用最大体积单体MVS(Maximization Volume Simplex)端元提取算法进行端元初选时存在相似端元光谱问题，提出一种光谱信息散度SID(Spectral Information Divergence)和光谱梯度角SGA(Spectral Gradient Angle)相结合以区分两个相似端元光谱的方法。该方法对经过端元初选之后的端元子集进行端元的二次选择，采用以SID_SG作为最相似端元选择的判据，除去相似端元，降低相似端元对解混精度的影响，利用全约束最小二乘法进行丰度估计。实验结果表明，提出的优化方法与传统方法相比，提高了端元的选择精度，重构影像与原始影像之间的均方根误差RMSE(Root Mean Square Error)也有所降低，分布更加均匀。该方法对高光谱遥感影像进行深度解译具有十分重要的意义。

关键词光谱解混端元初选二次提取除去端元解混算法

0引言

高光谱遥感技术已经广泛应用于军事和民用等各个领域。然而，阻碍高光谱遥感技术向定量化方向深入发展的主要因素是混合像元的普遍存在[1]。因为高光谱传感器有限的空间分辨率，每个像元对应着较大的区域，所以每个像元一般包含多种物质，多种物质光谱混合的混合光谱就是传感器探测的光谱[2]。解决混合像元问题最有效的方法就是混合像元分解[3]。高光谱混合像元分解一般包括两个方面，分别是端元提取[4]和丰度估计[5]。当前端元提取方法有N-FINDR[6]、顶点成分分析[7]VCA(Vertex Component Analysis)、迭代误差分析IEA(Iterative Error Analysis)[8]和纯像元指数PPI(Pure Pixel Index)等。传统的丰度估计算法主要有非限制性最小二乘法，非负限制性最小二乘法，和为一限制性最小二乘法和全限制性最小二乘法[9]。在进行高光谱解混时需要尽可能地选择出最精确的端元，降低解混时出现的误差。

本文结合两种常用的光谱区分方法，提出了一种将光谱信息散度和光谱梯度角相结合的方法用以区分两个相似端元光谱。该方法不仅能够从形状上反映两条光谱曲线的匹配程度，而且能够对两条有较小差异的光谱曲线做出较大的响应。首先利用最大体积单体端元提取算法进行端元的初次提取，然后再利用光谱信息散度和光谱梯度角SID_SG(Spectral Information Divergence_ Spectral Gradient)相结合的算法找出端元子集中存在的相似端元，去除相似端元，最后利用全约束最小二乘法FCLS(Full Constraint Least Square)进行端元的丰度估计，从而实现高光谱的混合像元分解。

1高光谱解混算法

1.1最大体积单体MVS

一般而言，高光谱图像中的每个像元都能看作是图像中各端元的线性混合的像元:

(1)

(2)

0≤ki≤1

(3)

式中，x是高光谱图像中任一P维的光谱向量(P为图像的波段数目)，E=[e1,e2,…,en]是P×N的矩阵，N是端元的数目，E中的每一列都是端元向量。c=(c1,c2,…,cn)T是系数向量，ci是像元x中端元ei所占比例，n表示误差项。

图1　两个波段三个端元

在误差项n非常小时，只有满足上面3个式子的点的集合，才能够构成高维空间的凸集，这些端元就坐落在凸面的单形体的各个顶点[10]。这里采用两个波段、三个端元作为例子，从而清晰地表明它们之间存在的几何关系，如图1所示。从图1中看出,三角形体的各个顶点分别是端元a、b和c，那么三角形内部的点就对应着图像中混合像元。所以提取高光谱图像的端元问题就转化成求单形体的顶点问题[11]。

从图1可以看出，图1中的端元a，端元b以及端元c所构成的三角形面积是图像中的任意三个像元所构成的面积集合中的最大的一个，表示如下:

S(a,b,c)=max{S(i,j,k)}

(4)

式中，S(i,j,k)表示的是由像元i、j和k所构成的三角形面积。所以寻找单形体的顶点的问题就转化为求最大体积单体的问题，体积公式如下:

(5)

(6)

式中，ei表示第i个端元的列向量，V表示的是由e1,e2,…,en这N个端元所构成的单形体体积，|·|表示行列式的运算符。因为这里用到了求行列式的运算，所以要求E必须为方阵。

1.2光谱信息散度SID

Ii(x)=-lgai

(7)

Ii(y)=-lgbi

(8)

通过式(7)和式(8)，可得y关于x的相对熵：

(9)

同样可以得到x关于y的相对熵：

(10)

因而x和y的光谱信息散度为：

SID(x，y)=D(x‖y)+D(y‖x)

(11)

光谱信息散度法的优点是可以对两条光谱进行整体上的比较。

1.3光谱梯度角SGA

光谱梯度角是在对两条光谱向量一阶求导的基础上求解它们的梯度向量，然后再得出两个梯度向量间的广义夹角大小[13]。

光谱曲线x的梯度向量为：

sGA(x)=[x2-x1,x3-x2,…,xn-xn-1]

(12)

光谱曲线y的梯度向量为：

sGA(y)=[y2-y1,y3-y2,…,yn-yn-1]

(13)

则梯度向量的广义夹角大小为：

(14)

该方法的优势在于能够反映光谱曲线的斜率变化，能反映出光谱曲线的局部特征的变化。

1.4信息散度与梯度角相结合

两条光谱曲线的锐角夹角范围是(0,π/2)，并且梯度向量的锐角夹角范围同样是(0,π/2)。因为正切函数的特性是当夹角是π/4时值为1；当夹角大于π/4时，其值随着夹角的增大迅速变大。根据正切的函数特征，将光谱梯度夹角控制在[π/4，π/2]范围内，当夹角值改变较小时，其结果就能有明显变化。因此为了不仅可以在整体上比较两条光谱曲线的匹配程度，而且可以从局部上区分两条光谱曲线的差异，提出将SID与SGA相结合的公式如下：

(15)

1.5均方根误差RMSE

RMSE[14]是均方根误差，是一个度量准则，用来评价两幅图像的相似性，均方根误差越小，说明二者越相近。

(16)

式中，假设端元丰度矩阵为S=(s11,s12,…,s1n;…;sp1,sp2,…,spn),真实丰度矩阵为Z=(z11,z12,…,z1n;…;zp1,zp2,…,zpn),sij和zij分别为矩阵S和Z的第i行第j列元素。

1.6丰度估计误差

(17)

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

2基于端元子集优选的高光谱解混

利用最大体积单体端元提取算法进行端元的初次提取，先去除部分端元，再进行端元的二次提取迭代循环，去掉相似的端元，剩下的端元为最佳端元，最后再利用全约束最小二乘法进行丰度的估计，最终实现高光谱的混合像元分解。端元初次提取的目的是为了去除部分多余的端元，端元二次选择的目的是进一步从原始高光谱影像中找到纯像元，从而准确求解各纯物质在混合像元中所占的比例。以SID_SG作为最相似的端元选择准则，能够选择出最相似的端元组，去除均方根误差较大的端元，保留均方根误差较小的端元，从而降低相似端元对解混精度的影响。

算法流程如图2所示。

图2　基于端元子集优选的高光谱解混算法流程图

端元二次选择优选算法步骤：

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

步骤2计算出去除端元A后的均方根误差。去除端元集合中的端元A，利用剩下的端元计算出均方根误差。

步骤3计算出去除端元B后的均方根误差。去除端元集合中的端元B，利用剩下的端元计算出均方根误差。

步骤4比较步骤2和步骤3中的均方根误差。将步骤2和步骤3中计算出的均方根误差进行对比，留下均方根误差比较小的端元。

步骤5重复步骤1～4，直至选出均方根误差最小并且求解丰度所运行的时间是最小的端元组，最终形成新的端元集合。

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

3实验结果分析

3.1端元数目不同时的对比实验

该实验采用的模拟图像大小是64×64，端元数目是3至8个，并且添加的信噪比SNR(Signal to Noise Ratio)为30 dB。在端元数目不同的情况下对比N-FINDR、单体增长算法SGA(Simplex Growing Algorithm)和最大体积单体(MVS)这三种端元提取算法的优劣，见表1所示。

表1　端元数不同时的实验对比

从表1中可以看出，N-FINDR、单体增长算法(SGA)、最大体积单体(MVS)这三种算法在端元数目不同的情况下对端元提取的结果进行相似性度量的对比。随着端元数目的增多， MVS算法的效果是最好的，SGA算法要优于N-FINDR算法，所以采用MVS算法进行端元初次选择的效果要优于SGA算法和N-FINDR算法。

3.2信噪比不同时的对比实验

该实验采用的模拟图像是64×64，端元数目是4条，并且添加的信噪比为10至50 dB。在不同的信噪比下对比N-FINDR，单体增长算法(SGA)，最大体积单体(MVS)这三种端元提取算法的鲁棒性，见表2所示。

表2　信噪比不同时的实验对比

从表2中可以看出，N-FINDR、单体增长算法(SGA)、最大体积单体(MVS)这三种算法在信噪比不同的情况下对端元提取的结果进行相似性度量的对比。对比表2中数据发现N-FINDR、SGA和MVS这三种算法的鲁棒性从高到低依次是MVS算法、SGA算法、N-FINDR算法，所以采用MVS算法进行端元的初次选择的稳定性要优于SGA算法和N-FINDR算法。

3.3分类精度对比实验

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

为了更好地证明SID_SG方法的优越性，利用实验数据将SID_SG方法与光谱角匹配SAM(Spectral Angle Mapping)法，光谱信息散度(SID)法，光谱梯度角(SGA)法这三种方法做对比实验，将它们在分类精度上进行对比验证，见表3所示。

表3　不同分类方法的精度分析

从表3中可以看出，分类精度最高的是SID_SG，达到78.96%，它对所调查的样本点正确分类的个数是30个。分类精度最低的是SAM，它的精度是60.51%，它对所调查的样本点正确分类的个数是24个。SID的分类精度低于SID_SG的分类精度，高于SAM的分类精度，它对所调查的样本点正确分类的个数是27个，其精度是69.75%。SGA的分类精度也低于SID_SG的分类精度，高于SAM的分类精度，它对所调查的样本点正确分类的个数是25个，其精度是62.81%。从这些数据可以看出，SID_SG的分类精度要高于SAM、SGA和SID，所以采用SID_SG算法提取相似端元的效果要优于SAM、SGA和SID。

3.4均方根误差及丰度估计时间对比

实验数据集为高光谱解混研究中经常采用的美国 Nevada 州南部沙漠地区的Cuprite数据[16]。该数据中共有224个波段，波长范围为0.37到2.48 μm，光谱分辨率是10 nm。

RMSE是均方根误差，是一个度量准则，用来评价两幅图像的相似性，均方根误差越小，说明二者越相近。

利用实验数据，根据本文优选算法分别提取出不同的端元数，得到的RMSE随着端元个数增加的变化情况，如图3所示。利用全约束最小二乘法解混求丰度运行的时间，见表4所示。根据图3中不同端元数目得出的均方根误差值以及表4中给出的在不同端元数目下利用全约束最小二乘法解混求解丰度所运行的时间值可以看出去除十个相似的端元是最适宜的，所以在端元二次提取时提取出12个端元是最合适的。

图3　端元数目和RMSE值的关系

端元数目求解丰度所运行时间(s)221156211097208681969318668176341660315578145421352112503

3.5整体丰度估计误差对比

为了研究噪声对本文算法性能的影响, 本文在混合光谱中加入了随机噪声。加入噪声的公式为：

(18)

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

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

图4　三种算法的丰度估计误差随信噪比变化情况

从图4中可以看出, 当信噪比SNR的范围是[12，100]时, 本文算法丰度估计误差fa是最小的，并且当信噪比达到100时，本文算法的丰度估计误差小于0.1，综上所述，本文算法要优于最大体积单体算法和全约束最小二乘算法。

3.6优选前后端元及误差分布对比

利用Cuprite数据进行优选前后的误差分布图对比。在端元初选时，利用最大体积单体端元提取算法进行端元的一次提取处理，提取出22个端元，得到其误差分布图，如图5所示。对22个端元进行端元的二次提取处理后，得到12个端元的误差分布图，如图6所示。

图5　初次提取误差分布图　　图6　二次提取误差分布图

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

4结语

本文提出了一种将光谱信息散度和光谱梯度角相结合的方法用以区分两个相似端元光谱。在利用最大体积单体端元提取算法进行端元初选之后，再采用以SID_SG作为最相似端元选择的判据进行端元的二次选择，除去相似的端元，降低相似端元对解混精度的影响。从实验的结果可以看出, 提出的基于端元子集优选的高光谱解混算法比传统算法提高了端元的选择精度，分布也更加均匀，同时也减少了重构影像与原始影像之间的均方根误差(RMSE)。该方法对高光谱遥感影像进行深度解译具有十分重要的意义。

参考文献

[1] Chen G,Song F,Wang J,et al.FRET spectral unmixing:a ratiometric fluorescent nanoprobe for hypochlorite[J].Chemical Communications,2011,24(4):2949-2951.

[2] Roth K L,Dennison P E,Roberts D A.Comparing endmember selection techniques for accurate mapping of plant species and land cover using imaging spectrometer data[J].Remote Sensing of Environment,2012,32(8):139-152.

[3] Gu Y,Wang S,Jia X.Spectral unmixing in multiple-kernel hilbert space for hyperspectral imagery[J].Gon and Rmo Nng Ranaon on,2013,51(10):3968-3981.

[4] Martin G,Plaza A.Spatial-spectral preprocessing prior to endmember identification and unmixing of remotely sensed hyperspectral data[J].Ld O N Ald Arh Obrvaon and Rmo Nng Jornal of,2012,5(12):380-395.

[5] Gonzalez C,Resano J,Plaza A,et al.FPGA implementation of abundance estimation for spectral unmixing of hyperspectral data using the image space reconstruction algorithm[J].Ld O N Ald Arh Obrvaon and Rmo Nng Jornal of,2012,5(3):248-261.

[6] Xiong W,Chang C,Wu C,et al.Fast algorithms to implement N-FINDR for hyperspectral endmember extraction[J].Ld O N Ald Arh Obrvaon and Rmo Nng Jornal of,2011,4(3):545-564.

[7] Krafft C,Diderhoshan M A,Recknagel P,et al.Crisp and soft multivariate methods visualize individual cell nuclei in Raman images of liver tissue sections[J].Vibrational Spectroscopy,2011,55(1):90-100.

[8] 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.

[9] Chang C I.An information-theoretic approach to spectral variability,similarity,and discrimination for hyperspectral image analysis[J].IEEE Transactionson InformationTheory,2000,45(5):1927-1932.

[10] 耿修瑞,张兵,张霞,等.一种基于高维空间凸面单行体体积的高光谱图像解混算法[J].自然科学进展,2004,7(14):810-814.

[11] 耿修瑞,赵永超,周冠华.一种利用单行体体积自动提取高光谱图像端元的算法[J].自然科学进展,2006,9(16):1196-1200.

[12] Jarosik N,Bennett C 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.

[13] Liu J,Zhang J.A new maximum simplex volume method based on householder transformation for endmember extraction[J].Ranaon on Gon and Rmo Nng,2012,50(6):104-118.

[14] 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.

[15] Brasco L,Pratelli A.Sharp stability of some spectral inequalities[J].Geometric and Functional Analysis,2012,22(1):107-135.

[16] Chen B,Vodacek A,Cahill N D.A novel adaptive scheme for evaluating spectral similarity in high-resolution urban scenes[J].Selected Topics in Applied Earth Observations and Remote Sensing,IEEE Journal of,2013,9(3):1376-1385.

收稿日期：2014-09-28。国家高技术研究发展计划项目(2012AA 12A405)；国家自然科学基金项目(61172144)。刘万军，教授，主研领域：数字图像处理，运动目标检测与跟踪。杨秀红，硕士。曲海成，博士。孟煜，硕士。姜庆玲，博士。

中图分类号TP75

文献标识码A

DOI:10.3969/j.issn.1000-386x.2016.07.058

HYPERSPECTRAL UNMIXING ALGORITHM BASED ON ENDMEMBER SUBSET OPTIMISATION

Liu Wanjun1Yang Xiuhong1Qu Haicheng1,2Meng Yu1Jiang Qingling3

1(SchoolofSoftware,LiaoningTechnicalUniversity,Huludao125105,Liaoning,China)2(SchoolofElectronicsandInformationEngineering,HarbinInstituteofTechnology,Harbin150006,Heilongjiang,China)3(TielingNormalCollege,Tieling112000,Liaoning,China)

AbstractFor the problem of maximisation volume simplex (MVS)endmember extraction algorithm that when applying in primary endmembers selection it would have similar endmembers spectra, we proposed an algorithm which distinguishes two similar endmember spectra by combining the spectral information divergence (SID) and the spectral gradient angle (SGA). This algorithm carries out secondary selection on the endmember subset derived from primary endmembesr selection, and adopts SID_SG rule as the criteria for selecting the most similar endmembers to remove the similar endmembers and to reduce the influence of similar endmembers on unmixing accuracy, and uses the full constraint least square for abundance estimation. Experimental results showed that proposed optimisation algorithm improved the accuracy of endmember selection than the traditional method, the root mean square error (RMSE) between reconstruction images and original image was rather reduced, and the distribution was more evenly as well. The algorithm is of very great significance to the deep interpretation on hyperspectral remote sensing image.

KeywordsSpectral unmixingPrimary endmembers selectionSecondary extractionRemoval of endmemberUnmixing algorithm