应用PCA和ICA法对岩心光谱的定量解译

蒋夕平，于瀚文，吴　芳，黄俊杰

(1.南京农业大学理学院，南京 210095；2.南京炮兵学院，南京 211132；3.南京地质矿产研究所，南京 210016)



应用PCA和ICA法对岩心光谱的定量解译

蒋夕平1，于瀚文2，吴芳1，黄俊杰3

(1.南京农业大学理学院，南京 210095；2.南京炮兵学院，南京 211132；3.南京地质矿产研究所，南京 210016)

摘要：岩心光谱属于混合像元光谱，维度多、数据量大，在端元数量、端元光谱及混合矩阵未知的情况下，定量解译岩心光谱以提取岩心所含矿物信息难度大，噪声的存在使问题更加复杂。文章应用PCA和ICA法定量解译岩心光谱主要有三步：采用PCA法预处理混合像元光谱矩阵，在新的特征空间中保留特征值较大的少量特征矢量，有效滤除能量较小的成分，只保留主要成分信息，同时滤除数据中的噪声；采用ICA法分离岩心混合像元，得到混合像元中的端元光谱集，通过矿物识别获取岩心矿物成分；对端元光谱进行归一化处理，基于线性光谱混合模型进行丰度反演，得到岩心混合像元中各端元丰度。通过研究仿真数据获得光谱定量解译的一般规律，创建有效的算法模型，再处理实际测量的岩心光谱，获得了比较理想的结果。

关键词：PCA；ICA；岩心光谱；混合像元分离；丰度反演

0引言

岩心是地质找矿和地质研究工作中常用的重要实物，通过分析钻孔岩心中矿物类型及含量的变化，能够为地质找矿提供信息，但岩心实物保存、查询及调阅存在诸多不便。高光谱具有很高的光谱分辨率，能在电磁波的可见光、近红外、中红外和热红外波段获取大量非常窄的光谱波段信息，被广泛应用于地学领域。利用光谱学原理，由岩心的光谱特性可获取岩心结构、成分及含量，以及岩心形成时及后期变化的物化环境信息。采用高性能的岩心光谱扫描仪，快速获取岩心图像和岩心光谱，解决了岩心实物保存和分析过程中的诸多不便，为地质研究和地质找矿提供了便捷的信息获取方法。

由于受空间分辨率的限制，混合像元问题广泛出现在各类遥感高光谱图像[1-2]和合成孔径雷达图像[3]中。常用的混合像元分离方法首先需要提取端元，再进行混合像元分离及端元丰度计算，端元提取方法如极小体积变换(MinimumVolumeTransform，MVT)[4]、像元纯度指数(PixelPurityIndex，PPI)[5-6]等。岩心光谱扫描仪获得的岩心光谱具有光谱分辨率高、信息量丰富等特点，岩心由多种矿物组成，岩心光谱同样存在混合像元问题，由于矿物种类繁多，岩心光谱中端元类型较多，单个混合像元中端元数量较多。由于测量设备和测量环境影响，岩心光谱中存在随机噪声及背景信号[7-8]，并且岩心光谱缺乏统一标准[9]，使得不同测量条件下获得的同类光谱难以比较分析。解译岩心光谱目的是获取组成岩心的端元及端元丰度，前者展示岩心中矿物的类型，后者表示矿物在岩心中所占比重[10-11]。由于岩心一般不存在纯像元(也称端元)，遥感高光谱图像研究中的端元提取方法对提取岩心端元有较大的局限性。有的研究人员直接采用光谱库中的矿物光谱作为端元光谱进行混合像元分离[12]，并在混合像元分离时加入稀疏性约束[13]，但由于光谱库中的矿物光谱与岩心端元光谱之间存在差异，影响混合像元分离精度。独立成分分析(IndependentComponentAnalysis，ICA)可以在端元数目和端元光谱未知的情况下，通过混合像元盲分离获取端元信息，具有较高的实用价值。

岩心中一般有2～3种矿物含量较高，相应矿物光谱在岩心光谱中的能量较大，因此在用ICA方法获取岩心端元信息前，先对岩心光谱预处理，滤除含量较少的成分以及噪声(它们在岩心光谱中能量较小)。主成分分析(PrincipalComponentsAnalysis，PCA)可保留特征值较大的成分，有效滤除混合信号中能量较小的成分。本研究采用PCA方法对岩心光谱预处理，采用ICA方法获取岩心端元光谱并进行矿物识别，进而实现矿物丰度反演。

1岩心光谱解译程序

1.1光谱预处理

PCA是一种常用的多元统计分析技术，该方法能够去除岩心光谱波段间的二阶相关性，通过保留少量特征值较大的主成分实现岩心光谱降噪减维[14]。将经过降噪减维处理的岩心光谱变换回原来的特征空间，即可滤除岩心光谱中能量较小的成分和噪声。

1.2混合像元分离

ICA是由盲源分离(BlindSourceSeparation，BSS)技术发展起来的一种多道信号处理方法，它以非高斯源信号为研究对象，在源信号统计独立的假设下，分离出隐含在混合信号中的独立信号。ICA已在语音信号分离方面得到广泛的应用，在核磁共振成像[15]、光谱遥感图像分类[16]、图像特征提取[17]等领域近年也得到较好应用，但在岩心光谱矿物识别和矿物丰度反演领域还少有应用。采用ICA算法分离岩心混合像元，将分离得到的独立分量作为端元估计光谱进行矿物识别和矿物丰度反演，突破了端元未知对岩心混合像元分离的限制。

1.3矿物丰度反演

ICA方法得到的独立分量具有顺序和幅度的不确定性。独立分量顺序的不确定性不影响矿物丰度反演，但独立分量幅度的不确定性直接影响丰度反演结果，必须对独立分量进行归一化处理。首先，计算独立分量与标准矿物光谱库中相应矿物的光谱角余弦，如果光谱角余弦为负值，将独立分量乘以-1；其次，对独立分量的各数据项进行处理，保证其值在[0,1]。经过归一化处理后的独立分量，将作为端元估计光谱用于矿物丰度反演。

2实验方法

研究中首先采用仿真数据(从光谱库中选取数种标准矿物光谱，将它们按照不同比例进行混合，生成数个混合像元)，在仿真数据无噪声、有高斯白噪声的2种情况下，分别采用以下2种方法进行处理：直接用ICA分离混合像元并进行矿物识别和丰度反演；对混合像元进行PCA处理，再用ICA分离混合像元进行矿物识别和丰度反演。所得结果分别与理论值比较。

3实验结果

3.1仿真数据

从美国地质勘查局(USGS)标准矿物光谱库中选取绿帘石、伊利石和滑石3种矿物，矿物光谱见图1。将3种矿物光谱按不同比例进行混合，生成5个混合像元，矿物光谱有420个波段，混合像元光谱数据矩阵的每一列为一个混合像元光谱。

3.1.1无噪声混合像元分离

直接采用ICA分离A420×5，得到3个独立分量IC1，IC2和IC3，波形见图2，它们与绿帘石、伊利石和滑石之间的光谱角余弦见表1。

表1中，IC1与滑石的光谱角余弦绝对值最大，两者的光谱角最小，IC2与绿帘石光谱角最小，IC3与伊利石光谱角最小。据图1、图2的波形及表1中光谱角余弦得知，IC1为滑石的估计光谱，IC2为绿帘石的估计光谱，IC3为伊利石的估计光谱。

将IC1，IC2，IC3进行归一化处理，再进行丰度反演，结果见表2。混合像元中丰度较大的绿帘石、伊利石理论值和计算值误差较小，丰度较小的滑石理论值和计算值误差较大。比较不同混合像元中IC2丰度值，从大到小为4，5，1，3，2(混合像元序号)，与绿帘石丰度理论值大小排列顺序完全一致，IC3(对应于伊利石)也有同样的结果，IC1(对应于滑石)有误差。

表1IC1，IC2和IC3与绿帘石、伊利石和滑石的光谱角余弦

Table1Spectralanglecosinebetween

IC1, IC2, IC3 and Epidote, Illite, Talc

图1　绿帘石、伊利石和滑石的光谱Fig.1　Spectrum of epidote, illite, talc

图2　IC1，IC2和IC3的光谱Fig.2　Spectrum of IC1,IC2,IC3

理论值绿帘石伊利石滑石计算值IC2IC3IC1混合像元10.50.420.080.570.30.13混合像元20.30.650.050.150.610.24混合像元30.340.60.060.240.540.21混合像元40.650.330.020.790.160.06混合像元50.570.40.030.650.250.1

对A420×5进行PCA处理，计算A420×5协方差矩阵的特征值λ1，λ2，…，λ420和特征矢量e1，e2，…，e420，最大的3个特征值为λ1(2.412)，λ2(6.964E-004)，λ420(1.547E-013)，相应的特征矢量为e1，e2，e3。由于λ1，λ2比其他特征值大得多(两者之和大于特征值总和的99%)，数据矩阵隐含的信息主要集中在第一个主成分e1和第二个主成分e2上。

实验结果显示，当保留1个主成分时，混合像元中能量最小的独立分量(与滑石对应)被滤除，因此PCA能够滤除混合像元中能量较小的成分。表3中的IC5，IC6和表5中的IC7，IC8与对应矿物的光谱角较小，表4中的IC5，IC6及表6中的IC7，IC8在不同混合像元中丰度计算值变化规律，与表2中绿帘石、伊利石丰度理论值变化规律完全一致。

3.1.2有噪声混合像元分离

在A420×5中叠加高斯白噪声，噪声能量小于原信号能量。直接采用ICA分离，得到5个独立分量(IC9，IC10，IC11，IC12和IC13)。计算各独立分量中绿帘石、伊利石和滑石的光谱角余弦，IC9中滑石光谱角余弦最大(0.8146)，IC10中绿帘石光谱角余弦最大(0.9729)；IC11中伊利石光谱角余弦最大(0.8261)；IC12，IC13中绿帘石、伊利石和滑石的光谱角余弦都小于0.4，对应于混合像元中的噪声。

表3　IC4，IC5和IC6中绿帘石、伊利石

表4　IC4，IC5和IC6丰度反演结果

表5　IC7，IC8中绿帘石、伊利石和滑石的光谱角余弦

表6　IC7，IC8的丰度反演结果

将前面叠加噪声者进行PCA处理并保留2个主成分，再用ICA分离得到3个独立分量(IC14，IC15和IC16)，各分量中绿帘石、伊利石和滑石的光谱角余弦见表7，丰度反演结果见表8；如果PCA处理只保留1个主成分，ICA分离只能得到2个独立分量(IC17和IC18)，分量中绿帘石、伊利石和滑石的光谱角见表9，丰度反演结果见表10。

实验结果表明，PCA能够有效滤除混合像元中能量较小的成分和噪声，表4为未加噪声混合像元经过PCA和ICA处理后得到的端元丰度，表8为叠加噪声混合像元经过PCA和ICA处理后得到的端元丰度，两表结果基本相同。经过PCA和ICA处理后，混合像元中各端元丰度的计算值与理论值有误差，但不同混合像元中同一种端元丰度的计算值与理论值变化规律相同，这对分析钻孔不同深度岩心中矿物含量变化具有意义。

3.2实际测量数据

采集岩心样品(钻孔位于南京江宁)，通过岩心光谱扫描仪获取岩心光谱，波长1 300～2 500nm，波长间隔2nm，测点直径10mm，相邻测点距离50mm。由于测量设备和测量环境的影响，岩心光谱中存在噪声，取相邻5个岩心光谱组成光谱数据矩阵D601×5。

表7　IC14，IC15和IC16中绿帘石、伊利石和滑石的光谱角余弦

表8　IC14，IC15和IC16丰度反演结果

表9　IC17和IC18中绿帘石、伊利石和滑石的光谱角余弦

表10　IC17和IC18丰度反演结果

对D601×5直接进行ICA分离，得到5个独立分量。如果岩心光谱数量更多，ICA分离得到的独立分量数目将更多，后续的矿物识别及矿物丰度反演比较困难。对D601×5进行PCA处理，保留2个最大特征值对应的特征矢量，ICA分离得到3个独立分量(IC19，IC20和IC21)(图3)。

分别计算IC19，IC20和IC21与USGS光谱库中标准矿物光谱的光谱角余弦并取最大值，IC19与柯绿泥石0.825，IC20与皂石0.942，IC21与松石0.692。将IC19，IC20和IC21进行归一化处理，再进行丰度反演，结果见表11。用同类分析软件处理D601×5，矿物识别结果基本一致。

图3　IC19，IC20和IC21的光谱Fig.3　Spectrum of IC19,IC20,IC21

IC19IC20IC21岩心10.410.360.23岩心20.450.220.33岩心30.450.210.34岩心40.440.230.33岩心50.450.250.30

4结论

PCA方法基于光谱概率分布函数的二阶统计信息对混合像元进行处理，通过在新的特征空间中保留特征值较大的主成分实现光谱降噪减维。PCA处理时保留1～2个主成分，其能量占总能量的99%以上，从而滤除能量较小的成分和噪声，保留混合像元的主要信息。在没有端元光谱先验知识的情况下，只要端元光谱是非高斯性信号且满足统计独立性，ICA能够有效分离混合像元得到独立分量，根据独立分量与光谱库中矿物光谱相似性进行矿物识别，基于光谱线性混合模型对岩心混合像元进行矿物丰度反演。运用PCA和ICA解译岩心光谱，能够快速获取岩心的矿物组分及含量信息，定量展示矿物随岩心深度变化的规律，对于充分利用宝贵的岩心资源具有重要意义。

参考文献：

[1]DundarMM,LandgrebeDA.Towardanoptimalsupervisedclassifierfortheanalysisofhyperspectraldata[J].IEEETransactionsonGeoscienceandRemoteSensing. 2004, 42(1): 271-277.

[2]ChangCI,RenH,ChangCC,etal.Estimationofsubpixeltargetsizeforremotelysensedimagery[J].IEEETransactionsonGeoscienceandRemoteSensing,2004,42(6):1309-l320.

[3]曹恒智，余先川，张立保. 基于SL-ICA算法的SAR图像混合像元分解[J]. 遥感学报，2009，13(2)：217-223.

[4]林娜，杨武年，刘汉湖．基于高光谱遥感的岩矿端元识别及信息提取研究[J]．遥感应用，2011(5)：114-117,.

[5]刘涛，严海英，于曼竹．高光谱影像混合像元分解比较研究[J]．地理空间信息，2011，9(3)：3-5.

[6]于钺，孙卫东．目标光谱指导下的高光谱图像混合像元分解方法[J]．高技术通讯，2012，22(3)：240-248.

[7]孙蕾，罗建书. 高光谱遥感图像微分域三维混合去噪方法[J]. 光谱学与光谱分析，2009，29(10)：2717-2720.

[8]刘庆杰，蔺启忠，王钦军，等. 基于连续统快速傅里叶变换的红外光谱处理技术[J]. 光谱学与光谱分析，2009，29(12)：3279-3282.

[9]张雪红，田庆久. 基于连续统去除法的冬小麦叶片氮积累量的高光谱评价[J]. 生态学杂志，2010，29(1)：181-186.

[10]CraigMD.MinimumVolumeTransformsforRemotelySensedData[J].IEEETransactionsonGeoscienceandRemoteSensing,1994,32(3):542-552.

[11]BoardmanJW,KruseFA,GreenRO.MappingtargetsignaturesviapartialunmixingofAVIRISdatainsummaries[C]∥FifthJPLAirbomeEarthScienceWorkshop,JPLPublication1995:23-26.

[12]IordacheMD,PlazaA,DiasJ．Ontheuseofspectrallibrariestoperformsparseunmixingofhyporspectraldata[C]．IEEEGRSSWorkshoponHyperspectralImageandSignalProcessing：EvolutioninRemoteSensing．Reykjavik，Iceland：IEEE，2010：1-4.

[13]BobinJ,MouddenY,StarckJL,etal．Sparsityconstraintsforhyperspectraldataanalysis：Linearmixturemodelandbeyond[M]．ProceedingsofSPIE：WaveletsⅫ．Bellingham，USA：SPIEPress,2009.

[14]JacksonGM,MasonIM,GreenhalghSA.Principalcomponenttransformsoftriaxialrecordingsbysingularvaluedecomposition[J].Geophysies,1991,56(4):528-533.

[15]HyvärinenA,OjaE.Independentcomponentanalysis:Algorithmsandapplications[J].NeuralNetworks,2000,13(4/5): 411-430.

[16]ZhangN,YuXC,DingGS.Anoptimalindependentcomponentanalysisapproachforfunctionalmagneticresonanceimagingdata[J].IEEEIIHMSP,2006:163-166.

[17]HeH,YuXC,PengWL.AComparisonbetweenFastICAandKernelICAinremotesensingimageryclassification[J].GeoscienceandRemoteSensingSymposium,2005, 3:1658-1661.

TheapplicationofPCAandICAmethodsinquantitativeinterpretationofmineralhyper-spectrum

JIANGXiping1,YUHanwen2,WUFang1,HUANGJunjie3

(1.College of Science，Nanjing Agricultural University，Nanjing 210095，China;2.Nanjing Artillery College，Nanjing 211132，China;3.Nanjing Institute of Geology and Mineral Resource，Nanjing 210016，China)

Abstract：The core hyper-spectrum is of mixed pixels with multi-dimensions and massive data. The quantitative interpretation of hyper-spectrum is extremely difficult when the end number of pixel, end spectrum and the mixing matrix are all unknown. If noise appears it is much complicated. with the application of PCA and ICA methods a three-step process is developed for the quantitative interpretation. Firstly PCA method is used to preprocess the pixel spectrum matrix and in new characteristic space several characteristic vectors with larger eigenvalues are kept and simultaneously noise filtered, Secondly, ICA method is used to un-mix the pixels, obtaining collection of end spectrum and identifying mineral components in the core. Finally, the end spectrum is normalized and abundance inversion carried out and each end abundance obtained. The quantitative interpretation law of hyper-spectrum is resulted from study on the simulation data. Algorithm model is built based on the law and the measured spectral of core processed by the model with ideal results.

Key Words：PCA; ICA; core spectra; separation of mixed pixel; abundance inversion

收稿日期：2014-12-12；责任编辑：赵庆

基金项目：国家重大科学仪器设备开发专项“岩心光谱扫描仪研发与产业化”(编号：2012YQ050250)和南京农业大学基本科研业务费专项资金(编号：KYZ201425)资助。

作者简介：蒋夕平(1963—)，男，讲师，主要从事数据建模与信号分析的教学、研究工作。

通信地址：江苏省南京市卫岗1号，南京农业大学理学院；邮政编码：210095；E-mail:jiangxp@njau.edu.cn 通信作者：吴芳(1974—)，女，讲师，硕士，主要从事数值计算、计算机仿真与数据建模研究。 江苏省南京市卫岗1号，南京农业大学理学院；邮政编码：210095；E-mail:wufang318@njau.edu.cn

doi：10. 6053/j. issn.1001-1412. 2016. 01. 016

中图分类号：P627

文献标识码：A