基于优化K-P-Means解混方法的高光谱图像矿物识别

孙 肖， 徐林林， 王晓阳， 田 野， 王 伟， 张中跃

(1.中国地质调查局廊坊自然资源综合调查中心，廊坊 065000； 2.中国地质大学(北京)土地科学技术学院，北京 100083)

0 引言

高光谱遥感具有“图谱合一”的特点，已广泛应用于矿物识别、固体矿产和油气勘探、环境保护及环境监测、植被分类及检测、月球和行星探测等方面[1-2]。现有高光谱遥感矿物识别技术主要采用基于矿物波谱曲线的光谱匹配技术，如混合像元分解技术、光谱角技术和混合调制匹配滤波技术等[3-5]。混合像元分解技术是一种常用的高光谱矿物识别技术。由于成像光谱仪的硬件问题，高光谱遥感图像普遍光谱分辨率高而空间分辨率低，容易导致影像像元出现同谱异物现象，因此必须研究解混方法予以克服。

目前常用的解混方法包括像元纯净指数、N-Findr、顶点成分分析、最小二乘端元提取算法等[6-9]。在前人工作的基础上，Roberts等[10]提出了经典的迭代光谱解混方法(MESMA)； Asner等[11]引入蒙特卡洛理论，将其与迭代解混理论有效结合起来，提出了AutoMUC方法； Song[12]介绍了一种基于贝叶斯决策的BSMA方法； 吴柯等[13]提出了基于神经网络的端元可变解混方法； 林红磊等[14]提出基于单次散射反照率的矿物高光谱稀疏解混方法。上述高光谱解混方法虽然反演精度有了很大的提高，但是经典的方法大都依赖于纯净像元，由于图像的空间分辨率和地表目标的复杂程度，纯净的像元是很难直接从高光谱图像中获取的； 此外，这些方法在解释混合像元的成像机理方面也还不完善。

Xu等[15]利用线性光谱混合模型来解释混合像元的成像机理，利用净化像元的概念提出K-P-Means算法来进行高光谱混合像元分解。该算法分2步迭代(丰度估计和端元优化)，通过基于高斯混合模型的期望最大(expectation maximization，EM)估计算法直到最终的端元估计收敛。由于K-P-Means算法估计端元和丰度的能力较强，因此，本文将该算法应用于矿物识别。该算法端元的优化借助于计算出的地物类别标签，由于地物类别标签可能不正确，导致计算出的净化像元存在多个类别，这种现象会导致端元优化不理想。Fischler等[16]提出的随机抽样一致性算法(random sample consensus，RANSAC)，可以有效剔除因分类标签不正确导致的净化像元中的异常值的影响。因此本文提出基于RANSAC的稳健的K-P-Means算法(RANSAC based robust K-P-Means，RR-K-P-Means)。在K-P-Means算法的基础上，利用RANSAC算法剔除异常影响，提高端元优化的精度。最终将该算法用于美国内华达州铜矿区的AVIRIS高光谱传感器Cuprite数据的矿物识别，将估计出的光谱曲线与美国地质调查局网站提供的矿物波谱曲线进行匹配，确定矿物种类。

1 高光谱混合像元分解方法

1.1 基本原理

1.1.1 线性光谱混合模型

线性光谱混合模型是一种常用的高光谱遥感图像表达方式[17]。假设高光谱图像像元集X由端元矩阵A和丰度矩阵S以及独立同分布的高斯噪声N组成，即

X=SAT+N，

(1)

(2)

(3)

式中：si(i=1,2,…，m)为E×1的丰度向量，用来表达第j个端元光谱向量aj(j=1,2,…,E)在P×1维高光谱像元xi中的贡献；m为图像像元个数；E为端元个数；P为波段数；n为高斯噪声向量。

1.1.2 K-P-Means模型

由于同一类的混合像元允许多种端元存在，为了更进一步移除丰度值较低的端元的影响，使用丰度值最大的端元(而不是用整个xi)来估计aj是很合理的。Xu等[15]把移除丰度值较低端元贡献后的xi称为净化像元。

K-P-Means模型是传统K-Means模型的衍生，考虑到实际中丰度非负，因此算法可以表述为：

(4)

式中k为地物类别数。

相比K-Means考虑物理过程的整体效果，K-P-Means探究影响观测的物理过程的源头，其目标函数可以表述为:

(5)

(6)

式中：l为类别；yi为第k类的净化像元向量。

基于上面描述的模型，由于未知数远远多于观测方程的个数，这是一个严重奇异的问题。EM算法对于奇异问题的参数求解具有较好的效果，该算法通过迭代的方法寻找统计模型的最大似然估计[18]。K-P-Means模型依据EM算法分2步迭代进行丰度估计和端元优化： 首先,通过非负最小二乘法进行丰度的估计; 然后,通过计算的丰度对端元进行优化。

1.1.3 丰度估计

给定{aj}，假设噪声满足高斯分布模型，则其噪声的概率密度函数p(n)为：

(7)

(8)

估计丰度的目标函数可以表达为：

(9)

由于方差可以通过均匀区域法获取，并作为权重计算丰度，地物类别li通过sik最大确定，即

(10)

因此，式(9)中丰度的估计本质上就是一个加权非负最小二乘(weight nonnegative least square， WNNLS)问题[19]。

1.1.4 端元优化

(11)

1.2 RR-K-P-Means算法

从上述K-P-Means算法的原理很容易发现，端元的优化完全依赖于由丰度确定的标签。类别标签错误或者存在异常值都会导致端元的优化精度降低。实践过程中，从净化像元值的直方图可以明显发现存在多个类别(图1)。

图1 净化像元直方图Fig.1 Histogram of purified pixels

净化像元值存在多个类别的现象降低了利用K-P-Means算法进行混合像元分解的精度。在计算机视觉以及其他很多研究领域，RANSAC算法对错误率超过50%的数据仍然能够得到理想的处理结果，是最有效的稳健估计算法之一[20]。因此，本文利用RANSAC算法对K-P-Means算法混合像元分解过程中产生的净化像元进行优化，剔除错误标签的影响，实现对K-P-Means算法的优化。

基于RANSAC算法的稳健的K-P-Means算法基本过程可以描述为： ①利用式(1)—(10)原理获取带有错误标签的净化像元； ②寻找一个模型(一般为线性模型)适应于假设的正确的净化像元(初始采用随机点)，利用寻找的净化像元估算该模型的参数(图2)； ③用第二步得到的模型去测试所有的其他端元，若某个端元适用于估计的模型，则认为它也是该类别的端元，其他端元满足的条件设置为大于最小端元差值(设定阈值)，如果有足够多的端元被归类为假设的正确类别的端元，那么估计的模型就足够合理，用所有假设的正确的端元去重新估计选取的模型，直到估计的正确的端元数量满足给定条件(设置最低错误率)为止； ④重复第二步和第三步过程来估计更加稳健的参数，选择更加合理的端元，将获取的合理的端元带入式(11)，获取优化后的端元。

图2 基于RANSAC算法的端元优化原理示意图Fig.2 Schematic diagram of the principle of end-memberoptimization based on RANSAC

2 实验与结果

2.1 仿真数据实验

从美国地质调查局公布的地物光谱库中，随机选择4条地物光谱曲线，按照如下程序混合成64×64大小的图像。利用这4个端元，首先把整个图像分成8×8的同质块； 然后通过7×7窗口大小的空间低通滤波对这些同质块降级； 再通过添加零均值高斯噪声对图像再降级，从而获取和真实情况更为接近的仿真数据。

丰度和端元的估计实验中，如果2个端元估计值的光谱角度距离(spectral angle distance，SAD)小于一个给定的数值τ或者预先规定迭代的次数iters，那么2步迭代就停止。本次实验中τ=0.01，iters=50。每次迭代中采用顶点成分分析的结果作为端元的初值[8]。利用RR-K-P-Means算法对端元进行优化时，模型采用线性模型，其他端元距离该直线的最小端元差值设置为5，最低错误率设置为15%。

分别利用K-P-Means算法和RR-K-P-Means算法进行丰度和端元的估计。用广泛使用的SAD和光谱信息散度(spectral information divergence, SID)来评价估计端元和真实端元的一致性[21]。为了便于表示，将SID的数值统一扩大10 000倍。利用结构相似性(structural similarity，SSIM)和峰值信噪比(peak signal-to-noise ratio，PSNR)来衡量利用估计的丰度和端元重新获得的图像与原图像的相似性[22-23]。

SAD和SID值越小，估计的光谱曲线与真实光谱曲线越相似。SSIM和PSNR越大，表明利用估计的丰度和端元重新获得的图像与原始图像越相似。从表1中数据可以看出，利用RR-K-P-Means算法估计的端元和丰度SAD达到0.73，SID达到3.1，估计值和真实值的一致性分别提高8.8%和13.89%。

表1 优化前后图像各指标对比Tab.1 Comparison of image parametersbefore and after optimization

利用估计的端元和丰度重新获取的图像和原始图像相似性SSIM达到0.997，PSNR达到35.67，相似性分别提高10.17%和62.80%。RR-K-P-Means算法明显优于传统的K-P-Means算法。从PSNR来看，RR-K-P-Means算法能更好地减弱噪声对估计的影响。

RR-K-P-Means算法可以剔除个数较少的异常类别的影响，而且选择的数据符合数据的分布形式(图3)。尽管RR-K-P-Means算法不能在大范围上改善分类结果，但是在细节上相比K-P-Means算法有了很大的改善。RR-K-P-Means算法估计的端元和丰度如图4所示。图4(a)—(d)为4种地物端元的估计值和真实值，两者吻合程度比较高； 图4(e)—(h)和(i)—(l)分别为4种地物的真实丰度和估计丰度，说明本文算法可以很好地计算丰度。较高精度的端元和丰度估计结果，为后续的矿物识别做好了准备。

图3 真实标签与优化前后2种算法确定的地物标签

图4 优化后估计的端元、丰度与真实值对比

2.2 真实高光谱数据

将RR-K-P-Means算法用于真实高光谱遥感图像数据的矿物提取。实验中采用美国内华达州铜矿区的AVIRIS高光谱传感器的Cuprite数据集(图5)，数据获取网站地址： http: //www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes。Cuprite数据集除了少数植被覆盖外大部分为含指示矿物的蚀变带，是进行高光谱遥感地质研究的典型样区。该影像数据为1997年6月19日机载AVIRIS获得，包括224个波段，地面瞬时视场约为20 m。剔除水汽影响严重的波段，共选取159个波段进行实验。

图5 内华达州铜矿区的Cuprite数据集Fig.5 Cuprite data sets of Nevadacopper mining area

利用RR-K-P-Means算法对Cuprite数据集进行端元和丰度的估计，通过计算估计出的端元和美国地质调查局矿物波谱库的波谱曲线的相关性确定端元所属的矿物类别。由于实际地物比较复杂，因此实验中的端元数设置为30，多于本地矿物的种类。选择相关系数大于0.75的端元进行研究。

通过相关系数的匹配得出8个类别的主要矿物，分别为： 黄钾铁矾、水铵长石、玉髓、叶蜡石、绿脱石、绿泥石、蒙脱石和白云母。研究区主要矿物绿脱石的端元估计结果如图6所示，波长在400～600 nm的端元估计结果和真实值匹配度一致，其他部分基本相同。

图6 绿脱石的端元估计值与真实值Fig.6 Estimated and true endmember values of nontronite

研究区主要矿物黄钾铁矾的端元估计结果如图7所示，估计端元波谱曲线与真实值基本相同。在实际应用中，利用RR-K-P-Means算法估计的端元可以有效提取矿物的波谱曲线。

图7 黄钾铁矾的端元估计值与真实值Fig.7 Estimated and true endmember values of jarosite

8种主要矿物的丰度图见图8。由图8可见，矿物聚集性比较好，与已知样本分布相比较，利用RR-K-P-Means算法获取的8种主要矿物基本与实际情况一致。

图8-2 识别出的8种矿物的丰度图

3 结论

针对K-P-Means算法中端元优化受异常值影响的问题，提出基于RANSAC的稳健的K-P-Means算法(RR-K-P-Means算法)，并将该算法用于矿物识别，通过仿真和真实数据实验，都取得了理想的效果。得出如下结论：

1)仿真实验证明，利用RR-K-P-Means算法估计的端元和丰度与真实值的一致性提高明显。利用估计的端元和丰度重新获取的图像和原始图像相似性提高较多。从PSNR来看，RR-K-P-Means算法能更好地减弱噪声对估计的影响。本文算法可以剔除个数较少的异常类别的影响，而且选择的数据符合数据的分布形式。

2)通过真实数据验证，利用RR-K-P-Means算法可以较好地识别研究区主要矿物，和美国地质调查局地物波谱库提供的标准地物波谱相比具有较高的一致性。在高光谱遥感图像矿物识别中，RR-K-P-Means算法可以有效提取多种地物光谱，矿物识别效果较好。

本文端元优化采用算数平均值，丰度较小的标签未参与计算，未来可以将加权平均的方法应用到本研究中。除此之外，RR-K-P-Means不仅在高光谱遥感图像矿物识别中可以得到很好的应用，同时可以推广到分类、去噪、超分辨率重建等研究中。