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

doi:10.6046/zrzyyg.2021215

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

孙肖^,¹, 徐林林^,², 王晓阳¹, 田野¹, 王伟¹, 张中跃¹

1.中国地质调查局廊坊自然资源综合调查中心,廊坊 065000

2.中国地质大学(北京)土地科学技术学院,北京 100083

Mineral identification from hyperspectral images based on the optimized K-P-Means unmixing method

SUN Xiao^,¹, XU Linlin^,², WANG Xiaoyang¹, TIAN Ye¹, WANG Wei¹, ZHANG Zhongyue¹

1. Langfang Natural Resources Comprehensive Survey Center, China Geological Survey, Langfang 065000, China

2. School of Land Science and Technology, China University of Geosciences(Beijing), Beijing 100083, China

通讯作者: 徐林林(1984-),男,博士,副教授,主要从事高光谱遥感理论研究。Email:linlinxu618@gmail.com。

责任编辑: 陈理

收稿日期: 2021-07-14 修回日期: 2021-08-19

基金资助:

中国地质调查局项目“辽阳市多要素城市地质调查”(DD20191025)

Received: 2021-07-14 Revised: 2021-08-19

作者简介 About authors

孙肖(1988-),男,硕士,助理工程师,主要从事高光谱遥感解译研究。Email: sunxiao@mail.cgs.gov.cn。

摘要

有效的高光谱混合像元分解方法可以提高矿物信息提取的精度。为进一步研究高光谱混合像元分解方法,采用线性光谱混合模型解释高光谱图像的成像机制,用不同矿物端元的线性组合表达混合像元。在最大似然估计的框架下,利用期望最大算法对混合像元的端元和丰度进行估计。针对端元提取易受异常值影响的问题,提出基于随机抽样一致算法的稳健的K-P-Means算法,优化端元提取过程。利用光谱角度距离和光谱信息散度评价估计端元和真实端元的一致性; 利用结构相似性和峰值信噪比衡量估计的丰度和端元重新获得的图像与原图像的相似性。多种仿真数据指标显示,该优化模型可以获得精度较高的端元和丰度估计结果。通过将提取的端元与美国地质调查局光谱库提供的矿物波谱曲线进行匹配来确定矿物类型,真实数据采用内华达州铜矿区的AVIRIS高光谱传感器Cuprite数据集。矿物提取结果显示,该模型对绿泥石等8种类别的主要矿物识别结果较好,矿物聚集性明显,与实际情况一致。该方法不仅可以较好地提取矿物信息,而且可以有效抵制噪声的影响。

关键词： 矿物; 高光谱图像; 端元; 解混; 稳健

Abstract

An effective unmixing method of hyperspectral mixed pixels can improve the precision of mineral information extraction. To further study such unmixing methods, this study explained the imaging mechanism of hyperspectral images using a linear spectral mixing model. The linear combinations of different mineral endmembers were used to express mixed pixels. The expected maximum (EM) algorithm was used to estimate the endmembers and abundance of mixed pixels under the framework of maximum likelihood estimation. A robust K-P-Means algorithm based on a random sampling consensus algorithm was proposed to improve the endmember optimization process, aiming to resist the impacts of anomalies on endmember extraction. The spectral angular distance and the spectral information divergence were used to assess the consistency between the estimated endmembers and the real endmembers. To obtain the similarity between the image and the original image, the structural similarity and the peak signal-to-noise ratio were used to measure the estimated abundance and endmembers. Various simulation data indicators show that the optimized model can obtain more precise estimations of endmembers and abundance. The mineral types were determined by matching the extracted endmembers with the mineral spectrum curves provided by the USGS spectral library. The actual data originated from the Cuprite data set of the AVIRIS hyperspectral sensor for the Nevada copper mining area. The results of mineral extraction showed that the model proposed in this study yielded satisfactory recognition results for eight types of main minerals including chlorite, which showed significant mineral aggregation and were consistent with the actual situation. Therefore, the method proposed in this study can extract precise mineral information while effectively resisting the impacts of noise.

Keywords： mineral; hyperspectral image; endmember; unmixing; robust

PDF (3634KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

孙肖, 徐林林, 王晓阳, 田野, 王伟, 张中跃. 基于优化K-P-Means解混方法的高光谱图像矿物识别[J]. 自然资源遥感, 2022, 34(3): 43-49 doi:10.6046/zrzyyg.2021215

SUN Xiao, XU Linlin, WANG Xiaoyang, TIAN Ye, WANG Wei, ZHANG Zhongyue. Mineral identification from hyperspectral images based on the optimized K-P-Means unmixing method[J]. Remote Sensing for Land & Resources, 2022, 34(3): 43-49 doi:10.6046/zrzyyg.2021215

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$ 组成,即

(1)

X = S A^{T} + N

(2)

(\begin{array}{l} x_{1}^{T} \\ x_{2}^{T} \\ ︙ \\ x_{m}^{T} \end{array}) = (\begin{array}{l} s_{1}^{T} \\ s_{2}^{T} \\ ︙ \\ s_{m}^{T} \end{array}) (a_{1}, a_{2}, \dots, a_{E})^{T} + (\begin{array}{l} n^{T} \\ n^{T} \\ ︙ \\ n^{T} \end{array})

(3)

x_{i} = \overset{E}{\sum_{j = 1}} s_{i j} a_{j} + n

式中: $s_{i} (i = 1,2, \dots, m)$ 为 $E \times 1$ 的丰度向量,用来表达第j个端元光谱向量 $a_{j} (j = 1,2, \dots, E)$ 在 $P \times 1$ 维高光谱像元 $x_{i}$ 中的贡献; m为图像像元个数; E为端元个数; P为波段数; $n$ 为高斯噪声向量。

1.1.2 K-P-Means模型

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

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

(4)

x_{i}^{k} = \overset{E}{\sum_{j = 1}} s_{i j} a_{j} + n, s_{i k} > 0

式中k为地物类别数。

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

(5)

{a_{k}, l} = \underset{l, m}{m i n} \overset{E}{\sum_{k = 1}} \sum_{l_{i} = k} {(y_{i} - a_{k})}_{2}

(6)

y_{i} = (x_{i} - \overset{E}{\sum_{j \neq k}} s_{i j} a_{j}) / s_{i k}

式中 $: l$ 为类别; $y_{i}$ 为第k类的净化像元向量。

K-P-Means排除低控制端元的贡献来估计 $a_{k}$ , $a_{k}$ 在提出的算法中被当作纯净高光谱像元 ${y_{i}^{k}}$ 的均值向量,因此该算法叫做K-P-Means。

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

1.1.3 丰度估计

给定 ${a_{j}}$ ,假设噪声满足高斯分布模型,则其噪声的概率密度函数 $p (n)$ 为:

(7)

p (n) = \frac{1}{\sqrt[]{{(2 π)}^{P} |Λ|}} e x p (- \frac{1}{2} n^{T} Λ^{- 1} n)

(8)

Λ = (\begin{array}{l} σ_{1}^{2} \end{array} \begin{array}{l} σ_{2}^{2} \end{array} \begin{array}{l} ⋱ \end{array} \begin{array}{l} σ_{P}^{2} \end{array})

式中: $Λ$ 为噪声协方差矩阵; $σ_{i}^{2} (i = 1,2, \dots, P)$ 为第 $i$ 波段的噪声方差。

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

(9)

\overset{︿}{S} = a r g \underset{A, S, Λ}{m i n} \sum_{i} (x_{i} - A s_{i})^{T} Λ^{- 1} (x_{i} - A s_{i}) \forall s_{i k} \geq 0

。

由于方差可以通过均匀区域法获取,并作为权重计算丰度,地物类别l_i通过 $s_{i k}$ 最大确定,即

(10)

l_{i} = a r g \underset{k}{m a x} {s_{i k}}

。

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

1.1.4 端元优化

给定{ $y_{i}^{k}$ },K-P-Means通过模型优化为:

(11)

y_{i}^{k} = a_{k} + n

式中 $a_{k}$ 为 $y_{i}^{k}$ 的算术平均值。

1.2 RR-K-P-Means算法

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

图1

新窗口打开| 下载原图ZIP| 生成PPT

图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

新窗口打开| 下载原图ZIP| 生成PPT

图2 基于RANSAC算法的端元优化原理示意图

Fig.2 Schematic diagram of the principle of end-member optimization based on RANSAC

2 实验与结果

2.1 仿真数据实验

从美国地质调查局公布的地物光谱库中,随机选择4条地物光谱曲线,按照如下程序混合成 $64 \times 64$ 大小的图像。利用这4个端元,首先把整个图像分成 $8 \times 8$ 的同质块; 然后通过 $7 \times 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 parameters before and after optimization

算法	SAD	SID	SSIM	PSNR
K-P-Means	0.80	3.6	0.905	21.91
RR-K-P-Means	0.73	3.1	0.997	35.67

新窗口打开| 下载CSV

利用估计的端元和丰度重新获取的图像和原始图像相似性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

新窗口打开| 下载原图ZIP| 生成PPT

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

Fig.3 Real label and feature label determined by two algorithms before and after optimization

图4

新窗口打开| 下载原图ZIP| 生成PPT

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

Fig.4 Comparison of estimated endmembers, abundance and true value after optimization

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

新窗口打开| 下载原图ZIP| 生成PPT

图5 内华达州铜矿区的Cuprite数据集

Fig.5 Cuprite data sets of Nevada copper mining area

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

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

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 绿脱石的端元估计值与真实值

Fig.6 Estimated and true endmember values of nontronite

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

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 黄钾铁矾的端元估计值与真实值

Fig.7 Estimated and true endmember values of jarosite

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

图8

新窗口打开| 下载原图ZIP| 生成PPT

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

Fig.8 Abundance map of the identified 8 minerals

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不仅在高光谱遥感图像矿物识别中可以得到很好的应用,同时可以推广到分类、去噪、超分辨率重建等研究中。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

李根军, 杨雪松, 张兴, 等.

ZY1-02D高光谱数据在地质矿产调查中的应用与分析

[J]. 国土资源遥感, 2021, 33(2):134-140.doi: 10.6046/gtzyyg.2020190.