一种新的高光谱遥感蚀变矿物提取算法——MSSSt相似光谱匹配组合法

韩海辉，王艺霖，张转，任广利，杨敏

(1.中国地质调查局 西安地质调查中心，西安 710054；2.西北地质科技创新中心，西安 710054；3.长安大学，西安 710054)

0 引言

根据岩石和矿物的光谱特征，利用高光谱图像提取蚀变矿物已成为一个遥感地质研究中的一个热点[1-2]，这是因为应用高光谱成像数据有可能直接识别与成矿作用密切相关的蚀变矿物以及围岩蚀变的范围和强度，而这类信息常常能为矿产勘查提供示矿依据[3-6]。

高光谱遥感时代源于2000年后Hyperion数据的大量获取[1-2]。这一时期，遥感蚀变矿物识别处理方法有了新突破，出现了光谱匹配识别技术与亚像素技术[1]。光谱匹配技术中使用最广的是由Kruse提出的光谱角填图法(SAM)[1，7]，而Chang认为光谱信息散度法(SID)要优于光谱角法[8]，Xu则通过对光谱吸收特征参量的条件约束，提高了光谱特征拟合(SFF)法对矿物的区分能力[9]。亚像素识别方法则是利用已知端元光谱，从一个混合像元中分解出所含矿物的相对丰度信息再进行光谱识别分类。近年在端元光谱识别中使用最广的方法是混合调制匹配滤波(MTMF)，它是线性波谱混合理论与目标信号匹配滤波处理法相结合的产物，该方法综合了匹配滤波不需要其他背景端元光谱和混合调制技术可以减少出现虚假信号概率的优点[1-2，7]。近年，学者又陆续提出了光谱信息散度-光谱角(SID-SAMtan)、光谱信息散度-光谱相关角(SID-SCAtan)、光谱信息散度-光谱梯度角(SID-SGAtan)等光谱匹配组合方法，研究结果显示这类组合法进一步提高了矿物的区分能力[10-12]。

这些高光谱方法取得了良好的应用效果，如Mars等利用光谱特征拟合法，在美国内华达和加州帕斯山两个赤铜矿蚀变区域内提取出泥化蚀变矿物(高岭石、明矾石、地开石)、青盘岩化矿物(方解石、绿泥石、绿帘石)和绢英岩化矿物(绢云母)[13]；Govil利用光谱角填图法提取了印度喜马拉雅山脉Kuta地区的绿泥石、白云母、伊利石、针铁矿，为复杂地貌区的矿产勘查提供了蚀变矿物信息[14]；Carrion等利用混合调制匹配滤波法提取了秘鲁南部Chapi Chiara金矿区的泥化(伊利石、绢云母等)和青磐岩化(绿泥石等)蚀变矿物，为识别岩性和热液蚀变区提供了帮助[15]；甘甫平等通过对岩矿谱形特征组合波段的主成分变换，识别了河北后沟金矿区的蚀变分布，并利用光谱角识别技术从机载高光谱图像上提取了高岭石、伊利石、绢云母等蚀变矿物[16]。

但基于光谱夹角的方法比较了光谱在整体形状上的相似性，忽略了局部特征的变化，容易产生“同物异谱”和“异物同谱”的现象；而光谱信息熵算法是基于概率统计理论，受光谱背景信息的影响比较大，复杂背景下往往难以反映细节特征从而会造成光谱区分度的降低。为了改进上述不足，本文提出了MTMFSFF-SIDSAMtan(MSSSt)组合方法。该算法将波谱特征拟合(SFF)和混合调制匹配滤波(MTMF)加入SID-SAMtan算法中，计算时既考虑了光谱整体的形状和信息熵特征，又通过包络线处理有效突出光谱曲线的局部吸收和反射特征，同时使用背景的协方差对光谱变化进行线性组合建模，从而最大化地放大与目标波谱匹配的像元信号，并抑制未知背景以提高输出的信号信噪比，以达到分离感兴趣目标和背景及噪声信息的目的。

1 MTMFSFF-SIDSAMtan(MSSSt)组合模型原理

MTMFSFF-SIDSAMtan(MSSSt)组合模型原理如下：假设存在2组光谱信号X，Y，其可分别表示为n维向量(X1，X2，Xi……，Xn)T和(Y1，Y2，Yi……，Yn)T，则MSSSt模型计算X和Y相似性可定义为：

(1)

计算结果是每个端元波谱对比每个像元的灰度图像，值越大表示越接近目标，两光谱的相似程度越好。理论上MSSSt组合模型可以将与目标波谱最为相似的像元信号突出增强显示，同时抑制背景及噪声的干扰，其中各因子的模型介绍如下。

1)光谱角SAM(X，Y)是把光谱看作多维矢量，通过计算两光谱之间的广义夹角确定光谱的相似程度，夹角越小，光谱越相似，其最大的特点是只使用了光谱的方向，而不使用其长度，因此夹角值与光谱向量的模无关，只比较光谱在形状上的相似性，故SAM分类对亮度影响不是很敏感[17-18]。计算公式为：

(2)

2)光谱角信息散度 SID(X，Y)是基于信息理论衡量光谱之间的差异，其定义如下[17-18]：

(3)

式中：D(X‖Y)为X关于Y的相对熵，D(Y‖X)为Y关于X的相对熵；p和q分别为X，Y光谱的概率向量。光谱角信息散度值越小，说明两光谱越相似。

3)波谱特征拟合(SFF)是一种基于吸收特征的方法，它使用最小二乘法对经过连续统去除处理的像元光谱曲线与参考光谱曲线进行匹配，即用各个离差的平方和最小来保证每个离差的绝对值都很小。当波谱曲线有明显的吸收特征时该方法的识别效果很好，这也弥补了SAM不能有效体现曲线局部变化特征的不足。计算结果的值越大表示越接近目标，两光谱的相似程度越好[17-18]。

该方法计算时，对于每个像元会输出一个拟合值和一个均方差误差RMS值。本研究中，采用拟合值和与RMS的比值作为拟合度图像，这样既保证了高度拟合，又保证了误差较小。计算公式为：

(4)

式中：a，b可通过下述方程求解

(∑xi2)a+(∑xi)b=∑xiyi

(∑xi)a+nb=∑yi

(5)

4)混合调制匹配滤波 MTMF(X，Y)是在匹配滤波的基础上发展的，它将匹配滤波的优点(并不需要已知所有端元的信息)同混合理论中的物理条件限制(给定像素的信号是包含在该像素中的单一物质成分的线性组合)结合起来，使用线性波谱混合理论，来限制可行性混合的结果，相对单一匹配滤可以减小虚假信号出现的概率[17-18]。计算结果的值越大表示越接近目标，两光谱的相似程度越好。

该方法计算时，对于每个像元会输出一个MF像元匹配值和一个不可行性噪声RMS值。本研究中，采用MF匹配值与不可行性噪声的比值作为拟合度图像，这样既保证了高度匹配，又保证了噪声较小。

2 实验区概况与数据预处理

2.1 实验区地质概况

实验区选择在甘肃北山方山口地区。该区域自20世纪80年代以来相继发现了方山口、明金沟、新老金厂、明水河、拾金坡、花牛山等多个金、银、铜、铅锌的矿床(点)，使该区域呈现出金矿床(点)数量多、成矿类型多样、蚀变矿物类型丰富、矿化蚀变特征明显的特点[19-22]。同时，该区域植被覆盖稀疏，基岩出露良好，地表碎石残留原地，山体相对高差小，是开展遥感地质研究的理想场地。

研究区地质演化历史复杂，矿产资源丰富。地层从太古宇到中生界出露较全，岩性包括变质岩、碎屑岩、碳酸盐等；区内岩浆活动强烈，侵入岩类型多样，从超基性、基性到中酸性岩浆岩均有产出，侵入时代集中为华力西期；区内断裂构造也较为发育，空间展布规律为近东西向和北北东向，另外区内还发育有弧形断裂，已知矿化点的分布与断裂构造关系较紧密；区内存在的蚀变类型有硅化、黄铁矿化(褐铁矿化)、绢云母化、铁碳酸盐岩化和青磐岩化等[19，23-25]。

2.2 遥感数据预处理

CASI/SASI是加拿大ITRES公司研制的两个机载高光谱成像仪。本次使用的CASI/SASI数据来源于核工业北京地质研究院，飞行获取时间为2011年7月6日，影像质量较好，无云及阴影遮盖。

CASI和 SASI 数据共有 388 个波段，1～288(0.380～1.050 μm)为CASI 数据范围；289～388(0.950～2.450 μm)为SASI 数据范围。CASI/SASI数据在0.950～1.050 μm处为重叠波段，由于 SWIR 噪声大于 VNIR，因此通常使用 CASI 数据替代 SASI 数据重叠部分。波段数据替代中会同时去除受到水气吸收(1.400 μm和1.900 μm附近)影响和航拍中低性噪比情况下形成的坏波段。本文使用的CASI/SASI经处理后最终剩余116个波段(CASI：35个波段；SASI：81个波段)，光谱范围为0.404～2.435 5 μm。

CASI/SASI数据的预处理与光谱重建包括辐射定标、几何校正、正射校正、空间重采样、CASI/SASI波段配准、大气校正。其中辐射定标、几何校正和正射校正是利用系统自带的软件和飞行测量时获取的相关数据进行，大气校正采用Flaash大气校正模块进行，利用双线性插值法进行数据空间重采样。

3 MTMFSFF-SIDSAMtan(MSSSt)法与已知方法的比较

为了将MSSSt方法与已知方法进行比较，本文利用美国内华达州Cuprite矿区的AVIRIS(airborne visible/infrared imaging spectrometer)高光谱成像数据进行实验，利用相似度增量和相关光谱区分熵指标，定量比较MSSSt、MTMF、SFF、SIDSAMtan、SID、SAM 6种方法的光谱区分能力。选择Cuprite矿区AVIRIS数据的原因是该矿区自 1987 年以来开展了多次矿物填图研究，有较为精确的蚀变矿物分布资料可参考，而我国没有蚀变矿物专题调查图集。示例的AVIRIS数据共 50个 波段，光谱区域为1.99～2.5 um，平均光谱分辨率10 nm。

比较时首先参考已知蚀变矿物填图结果[18]，在AVIRIS图像上选取一个1×1的相对纯净明矾石像元(根据USGS标准光谱特征对比和已知蚀变矿物填图结果确定)，取其光谱曲线作为参考光谱，另外再选取一个 1×1的明矾石混合像元和干扰像元(高岭土组成)(根据USGS标准光谱特征对比和已知蚀变矿物填图结果确定)，取两者光谱曲线，然后按照0%、15%、30%、45%、60%、75%的比例依次将干扰像元光谱掺杂到明矾石混合像元光谱中，再利用不同区分方法计算掺杂后的各混合光谱与明矾石参考光谱的相似性，对比各方法对光谱变化的区分能力。图1显示了所取的纯净明矾石像元的光谱曲线，以及按不同比例掺杂干扰信息后对应的混合像元的6条光谱曲线，可以看出，随着干扰信息掺杂比例逐渐增大，混合光谱与纯净明矾石光谱的差异越来越大。

图1 掺杂不同比例干扰信息的明矾石混合像元光谱曲线

表1为不同方法计算的明矾石参考光谱与含不同比例干扰信息的混合像元光谱的相似度。可以看出，随着干扰信息掺杂比例逐渐增大，6种方法计算结果都在不同程度的发生变化，即MSSSt、MTMF、SFF计算的光谱相似度不断减小，而SIDSAMtan、SID、SAM方法计算的光谱相似度不断增大，这表明随着光谱失真程度的增大，6种方法对干扰信息的的区分能力均呈上升趋势。

为了更直观对比和评价各方法的区分能力，本文对每种方法以最小失真类型1为标准，比较其他类型的相似度对于类型1的增量，相似度增量越大表明该方法区分度越好。由于各方法计算结果的量纲不同，无法直接利用其相似度增量结果进行比较，因此本文先将各方法计算的增量结果进行标准化处理，标准化方法选用min-max标准化(min-max normalization)，即离差标准化，它是对原始数据进行线性变换，使所有增量结果值都映射到[0，1]之间，这样获取的数据才更具可比性。

图2为不同方法计算的相似度增量曲线，其中相似度增量越大表明该方法区分度越好。可以看出，MSSSt方法的相似度增幅在类型1至类型6中都明显高于其他方法，尤其是由类型1至类型2再至类型3时，尽管曲线仅仅发生了很细微变化，但MSSSt方法的反应还是很灵敏，表明该算法对相似光谱曲线的区分能力最强。此外，该结果中MTMF方法和SFF方法对相似光谱曲线的区分能力要相对强于SAM、SID和SIDSAMtan方法。

表1 不同方法的相似度计算结果

图2 不同方法计算的相似度增量曲线

为了进一步定量比较几种方法对相似光谱的区分能力，本文采用了相关光谱区分熵的评价标准[11，26]，利用其中的相关光谱区分概率和相关光谱区分熵两个标准进行评价。

1)相关光谱区分概率(RSDPB)。假设存在目标光谱t和参考光谱集S={s1 ，s2 ， …，sn}，则目标光谱关于参考光谱集的RSDPB 为

(6)

式中：m(t，si)表示对目标光谱和参考光谱采用某种区分方法进行区分比较的结果。RSDPB是进一步计算RSDE的基础数据，不同方法间该数值本身没有太多对比意义。

2)相关光谱区分熵(RSDE)。相关光谱区分熵定义为：

(7)

该熵值用于评价目标光谱和参考光谱相似性的不确定度， 其值越大， 不确定度越高；反之，则目标光谱被正确识别的机率越大。

表2为计算的不同方法的相关光谱区分概率和熵值。从评价结果来看，MSSSt方法和SIDSAMtan方法的RSDE值显著低于其他几种方法的计算结果，表明这2种方法相对更有效。其中，MSSSt方法的RSDE值仅为1.485 1，为几种方法中的最小值，显示了其光谱区分能力相对最强。

表2 不同方法的相关光谱区分概率和熵

综上，本文提出的MSSSt区分方法对混合像元光谱信息变化的敏感性更强，即对目标矿物光谱曲线的微小失真有较大的响应，且相关光谱区分熵也表明MSSSt法的区分能力是几种方法中最强的，可以增大目标光谱被正确识别的机率。就遥感蚀变矿物信息识别应用来说，尽管目前的高光谱影像拥有较高的光谱分辨率和空间分辨率，但成像形成的像元往往都是混合像元，而混合像元的光谱曲线就存在不同程度的失真现象，因此将本文提出的MSSSt方法应用于高光谱遥感蚀变矿物信息的识别和提取具有较高的可行性。

4 MTMFSFF-SIDSAMtan(MSSSt)法的实例验证

为验证MSSSt组合方法的区分能力，本文又采用MSSSt法对北山方山口试验区CASI-SASI数据进行蚀变矿物提取，然后根据野外调查验证中获取的蚀变矿物和蚀变岩石的地面实测数据(包括岩矿光谱数据、岩矿测试分析数据、岩石岩性野外验证数据)对提取结果的精度进行评价。图3为研究区高光谱蚀变矿物提取结果图(将各类蚀变信息叠加到SASI单波段图像上)，包括了褐铁矿、白云母(绢云母)、白云石、绿泥石、黄钾铁矾等9种蚀变矿物，其中各类蚀变矿物的相似系数阈值是通过经验值确定的。

图3 北山试验区CASI-SASI高光谱遥感蚀变矿物异常图

本文采用点样本检验法[27]对提取结果进行精度评价，从统计结果(表3)可以看出，9种蚀变矿物提取结果的总体精度较为理想，基本都能保持在90%左右，表明本次使用的MSSSt方法较为有效。其中，褐铁矿、黄钾铁矾和绢云母的精度相对较高，而白云石与方解石、绿泥石与绿帘石两组矿物的光谱吸收峰位十分接近，因此利用高光谱数据识别时仍然出现较多的误提结果。

表3 用MSSSt方法提取北山试验区蚀变矿物的混淆矩阵对比验证结果

注：1.褐铁矿(Lim)；2.黄钾铁矾(Jar)；3.绿泥石(Chl)；4.绿帘石(Epi)；5.白云石(Dol)；6.方解石(Cal)；7.高铝绢云母(H-Ser)；8.中铝绢云母(M-Ser)；9.低铝绢云母(L-Ser)

从以上分析结果可以看出，在北山这种岩石出露良好的戈壁荒漠区，采用MSSSt组合法提取蚀变矿物信息具有很好的使用效果，那么在其他区域该方法是否依据适用呢？为此，本文又在青藏高原的东昆仑纳赤台地区开展了尝试。纳赤台地区地处青海格尔木市辖范围，气候属大陆高原气候，地势起伏大，区内植被不发育，基岩裸露一般，地表普遍被第四系浅覆盖，成矿区带划分属于东昆仑南铜、钴、金、钨成矿亚带，是东昆仑地区重要的矿产资源基地。测试中，利用该区域2012年获取的机载高光谱CASI-SASI数据(经过了各项数据预处理)，采用MSSSt组合法提取了褐铁矿、绿泥石、绢云母、方解石等8种蚀变矿物信息(图4)，分布图显示蚀变矿物信息整体分布不均匀，其中以白云岩、方解石、绢云母、绿泥石、绿帘石较为发育，表现出与北东向构造线方向相一致的特点，而褐铁矿发育较弱，仅在局部地段呈带状、点状产出。

根据野外调查工作中获取的查证分析资料对该提取结果的精度进行了统计，由于该区域平均海拔在4 000 m以上，因此野外工作中获取的验证点数据相比北山区域要少。通过混淆矩阵统计结果可以看出(表4)，采用MSSSt组合法提取的8种蚀变矿物的正确率相比北山地区有一定降低，但是各蚀变矿物的总体精度还是超过了84%，显示出较好的提取效果。野外验证中也发现，只要地表覆盖少，基岩出露相对较好，那么该地段提取的蚀变矿物信息基本都能与地面实际状况吻合。这也表明，MSSSt组合法的提取效果主要取决于地表岩石和矿物的出露状况以及遥感图像上地物的光谱准确性，它的计算原理(主要依据光谱相似度计算)决定了它的适用性不会受到遥感数据源和试验区地理位置的约束，因此理论上来说该方法可适用于任何区域和任何高光谱数据。

图4 东昆仑试验区高光谱蚀变矿物信息分布图

表4 用MSSSt方法提取东昆仑试验区蚀变矿物的混淆矩阵对比验证结果

注：1.褐铁矿(Lim)；2.绿泥石(Chl)；3.绿帘石(Epi)；4.白云石(Dol)；5.方解石(Cal)；6.高铝绢云母(H-Ser)；7.中铝绢云母(M-Ser)；8.低铝绢云母(L-Ser)

5 结束语

在SIDSAMtan方法基础上，本文将混合调制匹配滤波和波谱特征拟合算法加入到其中，提出了MTMFSFF-SIDSAMtan(MSSSt)组合法。新方法既考虑了两光谱曲线的形状特征(夹角)和幅值特征(信息熵)的相似性，又通过局部光谱吸收特征和丰度特征的匹配突出对光谱细节特征的反应，从而能够最大化地放大与目标波谱匹配的像元信号，提高对相似光谱的区分能力。

利用MSSSt、MTMF、SFF、SIDSAMtan、SID、SAM 6种光谱区分方法对人工掺杂不同比例干扰信息的明矾石混合像元光谱曲线进行相似度增量指标的计算，结果发现MSSSt方法的相似度增幅明显高于其他方法，尤其是对曲线细微变化的反应最为灵敏，而6种方法的相关光谱区分熵指标对比也显示MSSSt方法的RSDE值显著低于其他几种方法的计算结果，表明MSSSt方法的光谱区分能力相对最强。

根据野外调查验证中获取的蚀变矿物和蚀变岩石的地面实测数据(包括岩矿光谱数据、岩矿测试分析数据、岩石岩性野外验证数据)，采用点样本检验法对MSSSt方法提取的北山高光谱蚀变矿物结果进行精度评价，结果表明9种蚀变矿物的总体精度基本都在90%左右，其中褐铁矿、黄钾铁矾和绢云母的精度相对较高；而东昆仑地区8种蚀变矿物的总体精度超过了84%，显示出较好的提取效果。野外验证中也发现，只要地表覆盖少，基岩出露相对较好，那么该地段提取的蚀变矿物信息基本都能与地面实际状况吻合，这表明本文提出的新方法行之有效。

MSSSt方法提取效果主要取决于地表岩石和矿物的出露状况以及遥感图像上地物的光谱准确性，它的计算原理(主要依据光谱相似度计算)决定了它的适用性不会受到遥感数据源和试验区地理位置的约束，因此理论上来说该方法可适用于任何区域和任何高光谱数据。