1 引言

高光谱遥感是20世纪80年代发展起来的新兴遥感成像技术。高光谱图像能够对地物的空间信息和光谱信息同时成像,各波段的空间图像叠加在一起形成一个数据立方体。由于高光谱图像往往具有数十甚至数百个波段,每个像元在光谱维度上提取的辐射值可以形成一条近似平滑的光谱曲线,因此可以利用高光谱图像丰富的光谱信息对地物进行精细分类和识别。然而,受图像传感器空间分辨率的限制及地物空间分布复杂多样性的影响,高光谱遥感图像中存在大量的混合像元。这些像元包含多种地表覆盖类型,其光谱特性由多种地物的光谱混合而成。因此,将混合像元分类到任何典型地物都是不准确的。这使得传统的像元级分类方法不再适用。混合像元的存在已成为制约高光谱遥感技术定量化应用的主要障碍。

混合像元分解技术能够获知一个混合像元中所包含的地物种类(端元)及其所占比例(丰度),可以大大改善高光谱图像分类的精度^[1-3]。目前,研究最多、应用最广泛的端元提取方法是基于凸面几何学的方法。自1993年Boardman首次将凸面单体理论用于端元自动提取并获得成功后,国内外学者对基于凸面单体的端元提取方法进行了广泛的研究。凸面单形体方法认为高光谱图像的每个像元在光谱特征空间中都对应一个样本点,这些样本点在特征空间中都分布在一个凸面单形体内部,其中纯像元位于单形体的顶点位置,混合像元位于单形体的内部。在这类算法中,端元提取的过程往往需要所有像元都参与运算,导致运算量较大,算法运行效率较低。为了简化参与运算的数据量,提高端元提取算法的运行效率,有学者做了一些初步的研究。杨可明等^[4]于2014年提出了一种结合光谱信息熵的内部最大体积(N-FINDR)改进算法,实验结果表明,改进后的N-FINDR算法在确保较高端元提取精度的同时,参与运算的数据量大大减小。徐君等^[5]2016年提出了一种利用图像的空间信息在特征空间中设计超球简化高光谱数据的方法,提高了端元提取算法的运行效率。本文在总结现有研究成果的基础上,提出了一种高光谱数据简化改进的非负矩阵分解端元提取方法。通过设置光谱信息熵的阈值,筛选出图像中可能存在纯净像元的图像同质区,只选择同质区中最具代表性的像元参与非负矩阵分解运算,从而减少了端元提取算法的数据量,提高了算法的运行效率。

2 光谱信息熵

1948年,Shannon借用热力学中热熵的概念建立了“信息熵”的概念,解决了信息大小的量化度量问题。信息熵是对信息源不确定性平均值的描述,可以用信息的数学期望来计算。

设信息源X={x₁,x₂,…,x_n}为离散的随机变量,各信息源发生的概率为P={p₁,p₂,…,p_n},且 ∑i=1np_i=1,这里n表示信息源的总数目,p_i表示第i个信息源发生的概率,那么信息熵可表示为

H(X)=Elogb1pi=-∑i=1npilogbpi,(1)

当(1)式中b=2,e,10时,信息熵的单位分别为bit,nat和dit。

光谱信息熵^[4,6]为信息熵扩展应用到高光谱图像中形成的概念。光谱信息熵的计算式为

H(s,t)=-∑j=1mp(xj)log2p(xj),(2)

式中:x_j是高光谱图像上位于(s,t)处的像元第j波段的灰度值;p(x_j)是第j波段上x_j出现的概率,p(x_j)=N_j(s,t)/N_j;N_j(s,t)是第j波段上灰度值与x_j相等的像元个数;N_j是第j波段上的像元总数;m为影像波段数。

在高光谱图像中,纯净像元和地物混合程度均匀的像元在每一波段上的灰度值基本相等,称为同质区;而其他混合程度差别较大的混合像元由于其组成成分及其所占比例各不相同,导致其在不同波段上的灰度值无法保持一直相等。

对于同质区的像元,(2)式中的p(x_j)会比较大,经过各波段的累加后,光谱信息熵H(s,t)也会较大;而对于其他混合程度差别较大的像元,p(x_j)会比较小,而且在不同波段上灰度值也不能一直保持相等,导致其光谱信息熵H(s,t)会较小。因此,通过计算光谱信息熵,将光谱信息熵值较大的区域划分出来,就可以划分出图像的同质区。在地物区块分布比较明显的高光谱图像中,这些同质区大多是纯像元分布区或者是地物混合程度均匀的区域。

3 非负矩阵分解端元提取

非负矩阵分解(NMF)是Lee等^[7]提出的一种矩阵分解算法,它将一个非负矩阵分解为两个非负矩阵的乘积。非负矩阵分解的数学模型与线性光谱混合数学模型十分相似,这为NMF算法应用于线性混合模型下的混合像元分解提供了可能。在线性光谱混合模型中,如果用X、A、S分别表示混合像元的观测值、端元矩阵和丰度矩阵,那么混合像元的非负矩阵分解就是要寻找元素均为非负的l行p列的端元矩阵A∈R^l×p,和p行n列的丰度矩阵S∈R^p×n,R表示矩阵的集合。也可以这样表述:混合像元的非负矩阵分解就是要寻找元素均为非负的l行p列的端元矩阵A和p行n列的丰度矩阵S,满足X≈AS。

非负矩阵分解算法是一种盲源分离方法,不需要假定图像中存在纯像元,能够同时提取端元和丰度,比较适用于光谱混合程度较高的图像^[8]。

NMF的分解过程是一个优化求解的过程,需要构建一个合理的目标函数,然后选择合适的优化法则来求得最优解。Lee等^[7]构造了两种目标函数,当信息中的噪声服从泊松分布时,以K-L(Kullback-Leibler)散度作为目标函数f

minM,Sf(M,S)=∑e,fXe,flgXe,f(MS)e,f-Xe,f+(MS)e,fs.t.Me,f>0,Se,f>0。(3)

式中:e和f分别为矩阵的行数和列数。

当信息中的噪声服从高斯分布时,以矩阵X和估计矩阵M,S之间的欧氏距离构建目标函数。

minM,Sf(M,S)=12‖X-MS‖2s.t.Me,f≥0,Se,f≥0。(4)

欧氏距离比K-L散度更被人熟知,因此,多数NMF分解算法都将其作为目标函数。NMF分解算法包括梯度下降算法、乘法更新规则算法和交替最小二乘算法。NMF混合像元分解算法等价于约束最优问题:

minM,Sf(M,S)=12‖X-MS‖2s.t.Me,f≥0,Se,f≥0,∑e=1pSe,f=1。(5)

与(4)式相比,(5)式多了“和为1”的约束,因此,在每一步迭代后还需要对丰度矩阵S进行归一化处理。

由于NMF目标函数具有明显的非凸性,因此只能求得局部最优解。在应用NMF进行混合像元分解时,一般在算法中增加一定的约束条件来避免算法陷入局部极小值,如最小体积约束^[9]、平滑性约束^[10]、分段平滑和丰度稀疏约束^[11]。

4 采用光谱信息熵简化数据改进NMF混合像元分解算法

计算高光谱影像的光谱信息熵,将具有相同光谱信息熵的像元划分为图像的同质区。在同一个同质区的像元具有相似的光谱特性,这些像元在光谱特征空间的散点图必然聚集成团。根据凸面几何学理论,其中最纯的像元是距离凸面单形体中心最远的像元。在每一个同质区内,只提取一个最纯的像元作为典型像元参与非负矩阵分解端元提取运算,剔除其余像元,简化参与运算的高光谱数据。

改进的NMF混合像元分解算法的具体步骤如下:1) 对图像进行去噪处理,统计每一波段上像元的灰度值,将灰度值彼此相等的像元找出来,并根据(2)式计算高光谱图像中每个像元的光谱信息熵;2) 将光谱信息熵值相等的像元划入一个同质区;3) 计算整幅图像所有像元的均值光谱作为凸面单形体的中心,在每个同质区中选取距离中心点最远的像元作为该同质区的代表性像元,参与NMF混合像元分解运算,剔除同质区中的剩余像元;4) 保留没有划入同质区的单个孤立像元,参与NMF混合像元分解运算。

算法的流程图如图1所示。

图 1. 数据简化改进的NMF混合像元分解算法流程

Fig. 1. Flow chart of NMF mixed pixel decomposition algorithm improved by data simplification

下载图片 查看所有图片

同质区中距离凸面单形体中心相对最远的样本点,距离凸面单形体的顶点也相对最近,因此,可以认为是同质区内最纯的像元。对于NMF算法,选择合适的初始值可以避免陷入局部最小,提高收敛速度和算法精度。本文首先利用Hysime算法^[12]确定初始端元的数目k,然后根据同质区的大小,从大到小,在每个同质区中选取一个最纯的像元,一共选取k个端元组成初始端元矩阵,然后用全限制最小二乘方法(FCLS)进行混合像元分解,将获得的丰度矩阵作为算法的初始丰度矩阵。

5 实验分析

实验采用的数据是1997年由机载可见光/红外成像光谱仪(AVIRIS)获取的美国内华达州Cuprite地区的高光谱图像,截取大小为300 pixel×300 pixel的区域。这个区域中有多种矿物质,并以纯像元的形式存在,这些矿物质的典型光谱都存在美国地质勘探局(USGS)光谱库中,因此,该数据在混合像元分解研究中的应用比较广泛。AVIRIS高光谱图像包含224个波段,波长范围为400~2500 nm,去除受水蒸气吸收影响较大和信噪比较低的波段,剩余188个波段^[13],做去噪处理。选择13、68、158波段作为红、绿、蓝通道,合成假彩色图像,如图2所示。

图 2. Cuprite地区的AVIRIS数据假彩色合成图

Fig. 2. AVIRIS data false-color image of Cuprite region

下载图片 查看所有图片

分别在图像中的每一波段上统计所有像元的灰度值,将灰度值彼此相等的像元找出来并做分组处理,统计p(x_j)值。根据(2)式计算每个像元的光谱信息熵,将光谱信息熵值相等的像元划入同质区。根据凸面几何学的原理,距离凸面单形体中心最远的像元纯度最高。计算整幅图像所有像元的均值光谱作为凸面单形体的中心点,分别在每个同质区中选取距离凸面单形体中心点最远的少量像元作为该同质区的代表性像元,参与NMF混合像元分解运算,剔除同质区中的剩余像元。同时,保留那些没有划入同质区的单个孤立像元,让其也参与到NMF混合像元分解运算中。原始影像有90000个像元,经过上述简化后实际参与端元提取计算的像元数仅为18678,数据量减少了约4/5。

利用Hysime算法确定初始端元的数目为9。通过计算模拟图像上所有像元的光谱信息熵,将像元划入不同的同质区。按照同质区的大小,由大到小选取其中的9个最纯像元作为初始端元组成初始端元矩阵,然后以FCLS分解获得的丰度矩阵作为算法的初始丰度矩阵。

采用最小体积约束(MVC)-NMF算法对数据简化前后的高光谱数据进行混合像元分解,得到9种地物覆盖类型——明矾石、铵长石、玉髓、沙漠地表、高岭石、白云母、蒙脱石、皂石、榍石的端元光谱如图3所示,其中,图3(a)、(b)分别是数据简化前、后MVC-NMF分解得到的端元光谱。

图 3. MVC-NMF算法端元提取结果。(a)数据简化前;(b)数据简化后

Fig. 3. Endmember extraction results of MVC-NMF algorithm. (a) Before data simplification; (b) after data simplification

下载图片 查看所有图片

将数据简化前、后MVC-NMF算法提取的端元光谱与USGS光谱库中对应的光谱进行对比分析^[14],并用光谱夹角距离(SAD)做定量化比较,如表1所示。可以看出,数据简化前提取的端元和数据简化后提取的端元与USGS光谱库中的标准光谱都非常接近。数据简化前9种矿物的光谱角均值为0.1785 rad,数据简化后9种矿物的光谱角均值为0.1789 rad,也就是说数据简化前、后,端元提取的精度基本不变。比较数据简化前、后MVC-NMF混合像元分解算法的运行效率,数据简化前需要用时2686.3 s,而数据简化后仅需539.8 s,算法运行时间减少了约4/5,因此,算法运行效率得到很大提高。

6 结论

高光谱数据是一个庞大的数据立方体,一般的混合像元分解算法需要所有数据都参与运算,运算量很大,运行速度很慢,导致算法运行效率较低,在很多场合下无法满足实时性的要求。利用光谱信息熵理论对高光谱图像进行划分同质区,在每个同质区中只选取最纯的像元参与混合像元分解运算,简化了参与运算的数据量。用MVC-NMF混合像元分解算法进行验证和分析。实验证明,数据简化后MVC-NMF算法的端元提取精度依然能够得到保证,而算法的执行效率却得到大大提高。这说明,在混合像元分解前对高光谱数据进行简化数据预处理,可以减小混合像元分解算法的运算量。这对于其他传统的混合像元分解算法也有借鉴和启发意义。