應用PCA和ICA法對巖心光譜的定量解譯

蔣夕平，于瀚文，吳　芳，黃俊杰

(1.南京農業(yè)大學理學院，南京 210095；2.南京炮兵學院，南京 211132；3.南京地質礦產研究所，南京 210016)

?

應用PCA和ICA法對巖心光譜的定量解譯

蔣夕平1，于瀚文2，吳芳1，黃俊杰3

(1.南京農業(yè)大學理學院，南京 210095；2.南京炮兵學院，南京 211132；3.南京地質礦產研究所，南京 210016)

摘要：巖心光譜屬于混合像元光譜，維度多、數據量大，在端元數量、端元光譜及混合矩陣未知的情況下，定量解譯巖心光譜以提取巖心所含礦物信息難度大，噪聲的存在使問題更加復雜。文章應用PCA和ICA法定量解譯巖心光譜主要有三步：采用PCA法預處理混合像元光譜矩陣，在新的特征空間中保留特征值較大的少量特征矢量，有效濾除能量較小的成分，只保留主要成分信息，同時濾除數據中的噪聲；采用ICA法分離巖心混合像元，得到混合像元中的端元光譜集，通過礦物識別獲取巖心礦物成分；對端元光譜進行歸一化處理，基于線性光譜混合模型進行豐度反演，得到巖心混合像元中各端元豐度。通過研究仿真數據獲得光譜定量解譯的一般規(guī)律，創(chuàng)建有效的算法模型，再處理實際測量的巖心光譜，獲得了比較理想的結果。

關鍵詞：PCA；ICA；巖心光譜；混合像元分離；豐度反演

0引言

巖心是地質找礦和地質研究工作中常用的重要實物，通過分析鉆孔巖心中礦物類型及含量的變化，能夠為地質找礦提供信息，但巖心實物保存、查詢及調閱存在諸多不便。高光譜具有很高的光譜分辨率，能在電磁波的可見光、近紅外、中紅外和熱紅外波段獲取大量非常窄的光譜波段信息，被廣泛應用于地學領域。利用光譜學原理，由巖心的光譜特性可獲取巖心結構、成分及含量，以及巖心形成時及后期變化的物化環(huán)境信息。采用高性能的巖心光譜掃描儀，快速獲取巖心圖像和巖心光譜，解決了巖心實物保存和分析過程中的諸多不便，為地質研究和地質找礦提供了便捷的信息獲取方法。

由于受空間分辨率的限制，混合像元問題廣泛出現在各類遙感高光譜圖像[1-2]和合成孔徑雷達圖像[3]中。常用的混合像元分離方法首先需要提取端元，再進行混合像元分離及端元豐度計算，端元提取方法如極小體積變換(MinimumVolumeTransform，MVT)[4]、像元純度指數(PixelPurityIndex，PPI)[5-6]等。巖心光譜掃描儀獲得的巖心光譜具有光譜分辨率高、信息量豐富等特點，巖心由多種礦物組成，巖心光譜同樣存在混合像元問題，由于礦物種類繁多，巖心光譜中端元類型較多，單個混合像元中端元數量較多。由于測量設備和測量環(huán)境影響，巖心光譜中存在隨機噪聲及背景信號[7-8]，并且?guī)r心光譜缺乏統(tǒng)一標準[9]，使得不同測量條件下獲得的同類光譜難以比較分析。解譯巖心光譜目的是獲取組成巖心的端元及端元豐度，前者展示巖心中礦物的類型，后者表示礦物在巖心中所占比重[10-11]。由于巖心一般不存在純像元(也稱端元)，遙感高光譜圖像研究中的端元提取方法對提取巖心端元有較大的局限性。有的研究人員直接采用光譜庫中的礦物光譜作為端元光譜進行混合像元分離[12]，并在混合像元分離時加入稀疏性約束[13]，但由于光譜庫中的礦物光譜與巖心端元光譜之間存在差異，影響混合像元分離精度。獨立成分分析(IndependentComponentAnalysis，ICA)可以在端元數目和端元光譜未知的情況下，通過混合像元盲分離獲取端元信息，具有較高的實用價值。

巖心中一般有2～3種礦物含量較高，相應礦物光譜在巖心光譜中的能量較大，因此在用ICA方法獲取巖心端元信息前，先對巖心光譜預處理，濾除含量較少的成分以及噪聲(它們在巖心光譜中能量較小)。主成分分析(PrincipalComponentsAnalysis，PCA)可保留特征值較大的成分，有效濾除混合信號中能量較小的成分。本研究采用PCA方法對巖心光譜預處理，采用ICA方法獲取巖心端元光譜并進行礦物識別，進而實現礦物豐度反演。

1巖心光譜解譯程序

1.1光譜預處理

PCA是一種常用的多元統(tǒng)計分析技術，該方法能夠去除巖心光譜波段間的二階相關性，通過保留少量特征值較大的主成分實現巖心光譜降噪減維[14]。將經過降噪減維處理的巖心光譜變換回原來的特征空間，即可濾除巖心光譜中能量較小的成分和噪聲。

1.2混合像元分離

ICA是由盲源分離(BlindSourceSeparation，BSS)技術發(fā)展起來的一種多道信號處理方法，它以非高斯源信號為研究對象，在源信號統(tǒng)計獨立的假設下，分離出隱含在混合信號中的獨立信號。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?；旌舷裨胸S度較大的綠簾石、伊利石理論值和計算值誤差較小，豐度較小的滑石理論值和計算值誤差較大。比較不同混合像元中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協(xié)方差矩陣的特征值λ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在不同混合像元中豐度計算值變化規(guī)律，與表2中綠簾石、伊利石豐度理論值變化規(guī)律完全一致。

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處理后，混合像元中各端元豐度的計算值與理論值有誤差，但不同混合像元中同一種端元豐度的計算值與理論值變化規(guī)律相同，這對分析鉆孔不同深度巖心中礦物含量變化具有意義。

3.2實際測量數據

采集巖心樣品(鉆孔位于南京江寧)，通過巖心光譜掃描儀獲取巖心光譜，波長1 300～2 500nm，波長間隔2nm，測點直徑10mm，相鄰測點距離50mm。由于測量設備和測量環(huán)境的影響，巖心光譜中存在噪聲，取相鄰5個巖心光譜組成光譜數據矩陣D601×5。

表7　IC14，IC15和IC16中綠簾石、伊利石和滑石的光譜角余弦

表8　IC14，IC15和IC16豐度反演結果

表9　IC17和IC18中綠簾石、伊利石和滑石的光譜角余弦

表10　IC17和IC18豐度反演結果

對D601×5直接進行ICA分離，得到5個獨立分量。如果巖心光譜數量更多，ICA分離得到的獨立分量數目將更多，后續(xù)的礦物識別及礦物豐度反演比較困難。對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方法基于光譜概率分布函數的二階統(tǒng)計信息對混合像元進行處理，通過在新的特征空間中保留特征值較大的主成分實現光譜降噪減維。PCA處理時保留1～2個主成分，其能量占總能量的99%以上，從而濾除能量較小的成分和噪聲，保留混合像元的主要信息。在沒有端元光譜先驗知識的情況下，只要端元光譜是非高斯性信號且滿足統(tǒng)計獨立性，ICA能夠有效分離混合像元得到獨立分量，根據獨立分量與光譜庫中礦物光譜相似性進行礦物識別，基于光譜線性混合模型對巖心混合像元進行礦物豐度反演。運用PCA和ICA解譯巖心光譜，能夠快速獲取巖心的礦物組分及含量信息，定量展示礦物隨巖心深度變化的規(guī)律，對于充分利用寶貴的巖心資源具有重要意義。

參考文獻：

[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]于鉞，孫衛(wèi)東．目標光譜指導下的高光譜圖像混合像元分解方法[J]．高技術通訊，2012，22(3)：240-248.

[7]孫蕾，羅建書. 高光譜遙感圖像微分域三維混合去噪方法[J]. 光譜學與光譜分析，2009，29(10)：2717-2720.

[8]劉慶杰，藺啟忠，王欽軍，等. 基于連續(xù)統(tǒng)快速傅里葉變換的紅外光譜處理技術[J]. 光譜學與光譜分析，2009，29(12)：3279-3282.

[9]張雪紅，田慶久. 基于連續(xù)統(tǒng)去除法的冬小麥葉片氮積累量的高光譜評價[J]. 生態(tài)學雜志，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；責任編輯：趙慶

基金項目：國家重大科學儀器設備開發(fā)專項“巖心光譜掃描儀研發(fā)與產業(yè)化”(編號：2012YQ050250)和南京農業(yè)大學基本科研業(yè)務費專項資金(編號：KYZ201425)資助。

作者簡介：蔣夕平(1963—)，男，講師，主要從事數據建模與信號分析的教學、研究工作。

通信地址：江蘇省南京市衛(wèi)崗1號，南京農業(yè)大學理學院；郵政編碼：210095；E-mail:jiangxp@njau.edu.cn 通信作者：吳芳(1974—)，女，講師，碩士，主要從事數值計算、計算機仿真與數據建模研究。 江蘇省南京市衛(wèi)崗1號，南京農業(yè)大學理學院；郵政編碼：210095；E-mail:wufang318@njau.edu.cn

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

中圖分類號：P627

文獻標識碼：A