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

(1.南京農(nóng)業(yè)大學(xué)理學(xué)院，南京 210095；2.南京炮兵學(xué)院，南京 211132；3.南京地質(zhì)礦產(chǎn)研究所，南京 210016)

?

應(yīng)用PCA和ICA法對巖心光譜的定量解譯

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

(1.南京農(nóng)業(yè)大學(xué)理學(xué)院，南京 210095；2.南京炮兵學(xué)院，南京 211132；3.南京地質(zhì)礦產(chǎn)研究所，南京 210016)

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

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

0引言

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

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

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

1巖心光譜解譯程序

1.1光譜預(yù)處理

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

1.2混合像元分離

ICA是由盲源分離(BlindSourceSeparation，BSS)技術(shù)發(fā)展起來的一種多道信號處理方法，它以非高斯源信號為研究對象，在源信號統(tǒng)計獨立的假設(shè)下，分離出隱含在混合信號中的獨立信號。ICA已在語音信號分離方面得到廣泛的應(yīng)用，在核磁共振成像[15]、光譜遙感圖像分類[16]、圖像特征提取[17]等領(lǐng)域近年也得到較好應(yīng)用，但在巖心光譜礦物識別和礦物豐度反演領(lǐng)域還少有應(yīng)用。采用ICA算法分離巖心混合像元，將分離得到的獨立分量作為端元估計光譜進(jìn)行礦物識別和礦物豐度反演，突破了端元未知對巖心混合像元分離的限制。

1.3礦物豐度反演

ICA方法得到的獨立分量具有順序和幅度的不確定性。獨立分量順序的不確定性不影響礦物豐度反演，但獨立分量幅度的不確定性直接影響豐度反演結(jié)果，必須對獨立分量進(jìn)行歸一化處理。首先，計算獨立分量與標(biāo)準(zhǔn)礦物光譜庫中相應(yīng)礦物的光譜角余弦，如果光譜角余弦為負(fù)值，將獨立分量乘以-1；其次，對獨立分量的各數(shù)據(jù)項進(jìn)行處理，保證其值在[0,1]。經(jīng)過歸一化處理后的獨立分量，將作為端元估計光譜用于礦物豐度反演。

2實驗方法

研究中首先采用仿真數(shù)據(jù)(從光譜庫中選取數(shù)種標(biāo)準(zhǔn)礦物光譜，將它們按照不同比例進(jìn)行混合，生成數(shù)個混合像元)，在仿真數(shù)據(jù)無噪聲、有高斯白噪聲的2種情況下，分別采用以下2種方法進(jìn)行處理：直接用ICA分離混合像元并進(jìn)行礦物識別和豐度反演；對混合像元進(jìn)行PCA處理，再用ICA分離混合像元進(jìn)行礦物識別和豐度反演。所得結(jié)果分別與理論值比較。

3實驗結(jié)果

3.1仿真數(shù)據(jù)

從美國地質(zhì)勘查局(USGS)標(biāo)準(zhǔn)礦物光譜庫中選取綠簾石、伊利石和滑石3種礦物，礦物光譜見圖1。將3種礦物光譜按不同比例進(jìn)行混合，生成5個混合像元，礦物光譜有420個波段，混合像元光譜數(shù)據(jù)矩陣的每一列為一個混合像元光譜。

3.1.1無噪聲混合像元分離

直接采用ICA分離A420×5，得到3個獨立分量IC1，IC2和IC3，波形見圖2，它們與綠簾石、伊利石和滑石之間的光譜角余弦見表1。

表1中，IC1與滑石的光譜角余弦絕對值最大，兩者的光譜角最小，IC2與綠簾石光譜角最小，IC3與伊利石光譜角最小。據(jù)圖1、圖2的波形及表1中光譜角余弦得知，IC1為滑石的估計光譜，IC2為綠簾石的估計光譜，IC3為伊利石的估計光譜。

將IC1，IC2，IC3進(jìn)行歸一化處理，再進(jìn)行豐度反演，結(jié)果見表2?；旌舷裨胸S度較大的綠簾石、伊利石理論值和計算值誤差較小，豐度較小的滑石理論值和計算值誤差較大。比較不同混合像元中IC2豐度值，從大到小為4，5，1，3，2(混合像元序號)，與綠簾石豐度理論值大小排列順序完全一致，IC3(對應(yīng)于伊利石)也有同樣的結(jié)果，IC1(對應(yīng)于滑石)有誤差。

表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進(jìn)行PCA處理，計算A420×5協(xié)方差矩陣的特征值λ1，λ2，…，λ420和特征矢量e1，e2，…，e420，最大的3個特征值為λ1(2.412)，λ2(6.964E-004)，λ420(1.547E-013)，相應(yīng)的特征矢量為e1，e2，e3。由于λ1，λ2比其他特征值大得多(兩者之和大于特征值總和的99%)，數(shù)據(jù)矩陣隱含的信息主要集中在第一個主成分e1和第二個主成分e2上。

實驗結(jié)果顯示，當(dāng)保留1個主成分時，混合像元中能量最小的獨立分量(與滑石對應(yīng))被濾除，因此PCA能夠濾除混合像元中能量較小的成分。表3中的IC5，IC6和表5中的IC7，IC8與對應(yīng)礦物的光譜角較小，表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，對應(yīng)于混合像元中的噪聲。

表3　IC4，IC5和IC6中綠簾石、伊利石

表4　IC4，IC5和IC6豐度反演結(jié)果

表5　IC7，IC8中綠簾石、伊利石和滑石的光譜角余弦

表6　IC7，IC8的豐度反演結(jié)果

將前面疊加噪聲者進(jìn)行PCA處理并保留2個主成分，再用ICA分離得到3個獨立分量(IC14，IC15和IC16)，各分量中綠簾石、伊利石和滑石的光譜角余弦見表7，豐度反演結(jié)果見表8；如果PCA處理只保留1個主成分，ICA分離只能得到2個獨立分量(IC17和IC18)，分量中綠簾石、伊利石和滑石的光譜角見表9，豐度反演結(jié)果見表10。

實驗結(jié)果表明，PCA能夠有效濾除混合像元中能量較小的成分和噪聲，表4為未加噪聲混合像元經(jīng)過PCA和ICA處理后得到的端元豐度，表8為疊加噪聲混合像元經(jīng)過PCA和ICA處理后得到的端元豐度，兩表結(jié)果基本相同。經(jīng)過PCA和ICA處理后，混合像元中各端元豐度的計算值與理論值有誤差，但不同混合像元中同一種端元豐度的計算值與理論值變化規(guī)律相同，這對分析鉆孔不同深度巖心中礦物含量變化具有意義。

3.2實際測量數(shù)據(jù)

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

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

表8　IC14，IC15和IC16豐度反演結(jié)果

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

表10　IC17和IC18豐度反演結(jié)果

對D601×5直接進(jìn)行ICA分離，得到5個獨立分量。如果巖心光譜數(shù)量更多，ICA分離得到的獨立分量數(shù)目將更多，后續(xù)的礦物識別及礦物豐度反演比較困難。對D601×5進(jìn)行PCA處理，保留2個最大特征值對應(yīng)的特征矢量，ICA分離得到3個獨立分量(IC19，IC20和IC21)(圖3)。

分別計算IC19，IC20和IC21與USGS光譜庫中標(biāo)準(zhǔn)礦物光譜的光譜角余弦并取最大值，IC19與柯綠泥石0.825，IC20與皂石0.942，IC21與松石0.692。將IC19，IC20和IC21進(jìn)行歸一化處理，再進(jìn)行豐度反演，結(jié)果見表11。用同類分析軟件處理D601×5，礦物識別結(jié)果基本一致。

圖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結(jié)論

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

參考文獻(xiàn)：

[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]. 遙感學(xué)報，2009，13(2)：217-223.

[4]林娜，楊武年，劉漢湖．基于高光譜遙感的巖礦端元識別及信息提取研究[J]．遙感應(yīng)用，2011(5)：114-117,.

[5]劉濤，嚴(yán)海英，于曼竹．高光譜影像混合像元分解比較研究[J]．地理空間信息，2011，9(3)：3-5.

[6]于鉞，孫衛(wèi)東．目標(biāo)光譜指導(dǎo)下的高光譜圖像混合像元分解方法[J]．高技術(shù)通訊，2012，22(3)：240-248.

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

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

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

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

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

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

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

中圖分類號：P627

文獻(xiàn)標(biāo)識碼：A