基于小波包变换和权重光谱角制图的岩心高光谱蚀变信息提取

田青林，潘 蔚，李 瑶，张 川，陈雪娇，余长发

(1.核工业北京地质研究院遥感信息与图像分析技术国家级重点实验室，北京 100029；2.德州农工大学土木工程系，德克萨斯州 77843)

0 引言

高光谱遥感图像能够在获取被测地物空间信息的同时，得到一个光谱区间内(如可见光—短波红外波段)几十甚至几百个连续的窄波段光谱信息，具有“图谱合一”的特点，是当前遥感领域重要的研究方向之一[1-3]，已经在精准农业、目标探测和地物分类等领域的应用中取得了飞速的发展[4-6]。近年来，高光谱技术手段也被广泛应用于区域地质填图和矿产勘查领域，这是高光谱遥感应用最成功的领域之一。利用高光谱数据进行矿物填图，提取矿化蚀变信息，在快速缩小找矿范围、圈定找矿远景区和优选靶区等方面发挥着越来越重要的作用[7]。便携式地物波谱仪和地面成像光谱系统的逐步推广应用，使得岩矿光谱采集与分析工作逐步为遥感地质工作者所熟悉，并与星载、机载和地面遥感平台共同形成了“星空地深”一体化的遥感立体探测体系，大大推动了岩矿光谱分析与矿物填图应用，尤其是利用地面成像光谱测量系统对钻孔岩心样本进行高光谱扫描测量，可以提取并展现蚀变矿物在岩心整体垂向的发育状况，为分析钻孔岩心蚀变矿物空间分带规律和研究深部成矿环境、预测深部矿化前景提供依据[8]。

光谱角制图(spectral angle mapper，SAM)算法是通过计算光谱之间夹角来确定二者之间整体相似性的算法，其突出特点是简单高效，并可以在一定程度上克服照度和地形的影响[9]。但由于SAM算法对光谱曲线的局部微小特征差异刻画不够精细，当出现矿物种类不同而光谱曲线特征相似的“异物同谱”现象时，区分效果较差，一些学者尝试通过设置局部光谱区间权重[10-11]和选取波段组合[12]等方法改进SAM算法。小波包分解是在小波分解的基础上进一步发展起来的一种更加完善的分解方法，能够同时对信号的低频和高频部分进行分解，并且这种分解既无冗余，也无疏漏[13]，因此能够对包含有大量中、高频细节信息的遥感影像、生物医学和语音等信号进行更好的时频局部化分析。

针对上述问题，本文提出了一种基于小波包分解和权重SAM的岩心高光谱蚀变信息提取方法，选取小波包分解系数中信息熵统计量构建特征矢量，通过设置权重得到差异特征区间，并结合SAM算法进行岩心高光谱蚀变信息提取研究，旨在提高岩心高光谱蚀变矿物的识别精度。

1 算法描述

SAM算法是利用图像像元光谱和参考光谱之间矢量的夹角大小来表示光谱曲线间的相似程度，并通过选择阈值对图像未知像元进行归类，以达到识别地物的目的。其计算公式为[14]

(1)

式中：x和y分别代表图像像元光谱和参考光谱；α为像元光谱与参考光谱二者之间的夹角，取值范围为[0，π/2]。α的值越小，说明待分类的未知像元光谱与参考光谱的相似度越高，划归为同类的可能性也就越高[15-16]。但由于SAM算法要求地物的整体光谱均参与计算，致使非特征光谱及噪声的比重会压制诊断性特征光谱对分类的影响，为后续分类中选择相似性阈值带来难度。

权重SAM算法的思路是在图像像元光谱与参考光谱差异较大的特征区间进行权重设置，以突出局部特征差异。光谱曲线包含的波段总数为N，假设差异较大的特征区间包含有N1个波段，对应的像元光谱和参考光谱反射率值分别为xav和yav，v∈[1,N1]，差异较小的普通区间包含有N2个波段，对应的像元光谱和参考光谱的反射率值分别为xbu和ybu，u∈[1,N2]，且N1+N2=N，将差异较大的特征区间添加权重系数γ，γ≥1，当γ=1时即为传统SAM算法，计算公式(1)变化为

(2)

(3)

信号总能量E可以分解为所有不同小波包子分量的能量和，即

(4)

相对小波能量lp反映了每个子分量中小波包能量的分布情况，可由子分量能量与信号总能量之比求得，即

lp=Ep/E。

(5)

信息熵用来衡量各种概率分布所包含信息量的大小，可以作为一个系统复杂程度的度量。每个子分量小波包信息熵统计量定义为[17-18]

Wp=-lplnlp。

(6)

通过计算得到每个像元光谱曲线X所对应的p个子分量小波包信息熵统计量，构造p维信息熵特征矢量T，即

T=(W1,W2,…,Wp)。

(7)

不同类别矿物的诊断性光谱特征不同，反映在各个子分量上具有不同的能量分布特性，使得各个子分量的小波包信息熵统计量也有所差异。因此，可以通过构造信息熵特征矢量来描述不同矿物的光谱曲线特征，并找到图像像元光谱与参考光谱信息熵矢量差异较大的特征区间，对其设置权重，采用权重SAM算法分别计算参考光谱和图像像元光谱小波包信息熵特征矢量的夹角，并根据夹角余弦值的大小进行蚀变矿物匹配，达到信息提取的目的。

2 数据与实验

2.1 数据获取

实验数据是利用HySpex岩心地面成像高光谱测量系统对我国南方某火山岩型铀矿区钻孔岩心样本扫描得到的短波红外(SWIR)波段数据，矿区热液活动期形成了碱性和酸性蚀变相互叠加的矿化和对称中心式近矿围岩蚀变。早期碱交代矿化蚀变主要包括钠长石化、绿泥石化和赤铁矿化等，晚期酸性矿化蚀变主要表现为萤石化和水云母化(伊利石或伊蒙混层黏土矿物)等，钻孔岩心样本岩性总体表现为碎斑流纹岩，斑晶主要为钾长石、斜长石、石英和黑云母。HySpex成像光谱仪包含可见光—近红外(VNIR)和SWIR波段2个传感器，能够同时采集获取0.4～2.5 μm光谱范围内的地面成像高光谱数据，其中SWIR传感器主要技术指标见表1。

表1 HySpex SWIR传感器技术指标Tab.1 Technical specifications of HySpex SWIR sensor

2.2 数据预处理

通过地面扫描得到的钻孔岩心高光谱数据，在采集过程中受仪器本身因素和外部环境影响，原始高光谱图像数据除含有岩心样本自身有效信息以外，还包含了噪声干扰信息[19]。因此，为了减少噪声对蚀变信息提取的干扰，需要对原始高光谱数据进行辐射校正、反射率转换和小波去噪等预处理。

具体操作为：采用系统自带软件，依据辐射定标参数完成辐射校正，得到辐射亮度数据；为了利用反射光谱特征提取矿化蚀变信息，采用经验线性法对辐射亮度值进行反射率转换；干扰信息在反射率光谱曲线上多表现为小的锯齿噪声，因此在不影响光谱曲线整体形态和诊断特征位置的前提下，运用小波变换方法对反射率数据进行光谱去噪，减弱或消除明显的锯齿噪声影响，提高信息提取的精度。

2.3 蚀变矿物端元光谱提取

端元光谱提取是高光谱遥感影像分类和目标识别的基础。以往研究中，矿物端元光谱的确定通常有2种方式：①使用光谱仪测量样品或者直接从矿物标准光谱库中获取；②直接通过数学方法从高光谱影像上选择端元，提取矿物的端元光谱[20]。本文采用第2种方法，采用ENVI软件最小噪声分离(minimum noise fraction，MNF)变换对高光谱数据进行噪声白化和降维处理，利用纯净像元指数(pixel purity index，PPI)方法，提取高光谱图像上的相对纯像元，将纯像元投影到N维可视化空间进行N维散度分析得到各端元光谱曲线，再结合标准矿物光谱库及专家知识人工确定最终的矿物端元光谱曲线。基于HySpex岩心高光谱数据提取并确定的端元矿物主要包括高岭石、绿泥石、地开石和伊利石(图1)。

图1 蚀变矿物端元光谱曲线Fig.1 Endmember spectrum curves of altered minerals

2.4 小波基选取

不同的小波基函数性质各异，对光谱信号的变换结果影响很大。使用不同的小波基函数进行小波包分解会产生不同的分解系数，进而影响所构建的小波包信息熵矢量及最终分类效果。而不同的分解层数同样影响着小波包信息熵特征矢量表征原始光谱曲线特征的好坏。因此，选取适用于钻孔岩心高光谱数据处理的最优小波基和最佳分解层数显得尤为重要。本研究中通过计算不同小波基函数和分解层数与矿物光谱曲线信息熵矢量夹角的关系来确定最佳小波基和分解层数。夹角余弦值越大，不同矿物类别间的可区分性越高，信息熵矢量更有利于刻画光谱曲线特征。通过对比最终确定daubechies4作为最优小波基，最佳分解层数为8层。

2.5 差异特征区间选取

利用权重SAM算法识别某种矿物时，可找到与其他矿物差异较大的特征区间，并对其设置权重系数γ，增加不同矿物之间的可区分性。在信息熵特征空间中，通过比较所构建的4种矿物端元光谱小波包信息熵矢量，采用差异较大的特征区间为小波包分解的前8个子分量信息熵(图2)。

图2 4种矿物部分子分量信息熵Fig.2 Partial subcomponent entropy of four kinds of minerals

从图2可以发现，4种矿物的小波包信息熵在此区间各个子分量上差异较大，表现出不同的分布特性，通过试验权重系数取γ=2。

3 结果与讨论

利用Matlab语言将本文算法编写程序实现，对实验数据进行处理。通过对岩心高光谱影像上提取的高岭石、绿泥石、地开石和伊利石4种矿物端元光谱曲线进行8层小波包分解，根据各自分解系数生成相对应的小波包信息熵特征矢量。依据SAM算法原理分别计算4种矿物端元光谱曲线两两间光谱矢量和小波包信息熵矢量夹角余弦值大小，结果如表2所示。通过对比发现，原始光谱曲线光谱矢量夹角余弦值均大于基于小波包分解系数生成的信息熵矢量夹角余弦值，经过小波包分解后的信息熵矢量增大了不同矿物类别间的夹角，提高了可区分度，更有利于刻画光谱曲线的特征。

表2 光谱矢量和信息熵矢量夹角余弦值Tab.2 Cosine for angles of spectral vectors and information entropy vectors

图3为利用岩心高光谱数据分别进行传统SAM算法(图3(b))以及本文方法(图3(c))的蚀变信息提取结果。

(a)岩心HySpex SWIR B80(R)，B200(G)，B48(B)假彩色合成图像

(b)传统SAM算法提取结果

(c)本文方法提取结果

图3 HySpex钻孔岩心高光谱数据蚀变信息提取结果Fig.3 Extraction results of alteration information using HySpex drill core hyperspectral data

从岩心图像提取结果可以看出，本文方法提取的蚀变矿物比传统SAM算法提取的范围更为精确细致。通过整个钻孔岩心提取结果分析发现，高岭石、地开石和明矾石等酸性蚀变矿物在钻孔上部含量多，下部含量少，而伊利石、绿泥石和碳酸盐岩等中、碱性蚀变矿物上部含量少，下部尤其是深部矿化段含量明显增多，总体分布符合本地区“上酸下碱”的蚀变特征。

为了定量评价2种方法提取结果的有效性，以钻孔柱状图、目视判别及像元光谱曲线为依据，在图像中随机选取5 000个像素点(包括背景值)作为真值进行检验，构建混淆矩阵，选取总体分类精度和Kappa系数指标进行精度评价，统计结果如表3所示。

表3 不同匹配方法分类精度统计Tab.3 Classification accuracy statistics of different matching methods

从表3可以发现，与传统SAM算法相比，本文方法的分类精度更高些，总体分类精度由71.51%提高到75.33%，Kappa系数由0.667 8增加到0.706 3。

通过对本文涉及的蚀变矿物及其光谱曲线分析发现，高岭石和地开石多伴生发育，光谱曲线特征接近，使得依据光谱相似性度量准则的传统SAM算法对其区分效果有限，而经过小波包变换后得到的信息熵矢量夹角被增大，同时通过设置差异特征区间，突出了局部特征信息，使得二者之间可区分性被进一步增强，分类精度得到了提高。

实验结果表明，本文提出的基于小波包变换和权重SAM算法的岩心高光谱蚀变信息提取方法，定义的信息熵特征矢量可以较好地表征原始光谱信息，结合差异特征区间的权重设置，增大了不同矿物类别间的可区分性，分类精度优于传统SAM算法，具有较好的适用性。

4 结论

本文运用小波包变换和权重SAM相结合的方法对HySpex岩心高光谱数据进行蚀变信息提取，并与传统SAM算法进行对比实验。主要结论如下：

1)本文方法提高了蚀变矿物提取的总体精度，由传统SAM算法的71.51%提高到75.33%，具有较好的适用性。

2)依据小波包分解系数计算得到的信息熵特征矢量可有效表征不同矿物光谱差异，结合差异特征区间的权重设置，增大了不同类别矿物间的可区分性。

3)提取得到的钻孔岩心矿化蚀变信息，可展现蚀变矿物在岩心整体垂向的发育状况，为进一步划分酸、碱蚀变分带，分析铀矿富集规律和预测深部矿化前景提供依据。

在后续的研究中拟对距离和相关系数等不同的匹配方法进行对比研究，同时对数据预处理中不同的光谱去噪方法进行试验，进一步提高信息提取精度。