基于小波包變換和權(quán)重光譜角制圖的巖心高光譜蝕變信息提取

田青林，潘 蔚，李 瑤，張 川，陳雪嬌，余長發(fā)

(1.核工業(yè)北京地質(zhì)研究院遙感信息與圖像分析技術(shù)國家級重點(diǎn)實(shí)驗(yàn)室，北京 100029；2.德州農(nóng)工大學(xué)土木工程系，德克薩斯州 77843)

0 引言

高光譜遙感圖像能夠在獲取被測地物空間信息的同時，得到一個光譜區(qū)間內(nèi)(如可見光—短波紅外波段)幾十甚至幾百個連續(xù)的窄波段光譜信息，具有“圖譜合一”的特點(diǎn)，是當(dāng)前遙感領(lǐng)域重要的研究方向之一[1-3]，已經(jīng)在精準(zhǔn)農(nóng)業(yè)、目標(biāo)探測和地物分類等領(lǐng)域的應(yīng)用中取得了飛速的發(fā)展[4-6]。近年來，高光譜技術(shù)手段也被廣泛應(yīng)用于區(qū)域地質(zhì)填圖和礦產(chǎn)勘查領(lǐng)域，這是高光譜遙感應(yīng)用最成功的領(lǐng)域之一。利用高光譜數(shù)據(jù)進(jìn)行礦物填圖，提取礦化蝕變信息，在快速縮小找礦范圍、圈定找礦遠(yuǎn)景區(qū)和優(yōu)選靶區(qū)等方面發(fā)揮著越來越重要的作用[7]。便攜式地物波譜儀和地面成像光譜系統(tǒng)的逐步推廣應(yīng)用，使得巖礦光譜采集與分析工作逐步為遙感地質(zhì)工作者所熟悉，并與星載、機(jī)載和地面遙感平臺共同形成了“星空地深”一體化的遙感立體探測體系，大大推動了巖礦光譜分析與礦物填圖應(yīng)用，尤其是利用地面成像光譜測量系統(tǒng)對鉆孔巖心樣本進(jìn)行高光譜掃描測量，可以提取并展現(xiàn)蝕變礦物在巖心整體垂向的發(fā)育狀況，為分析鉆孔巖心蝕變礦物空間分帶規(guī)律和研究深部成礦環(huán)境、預(yù)測深部礦化前景提供依據(jù)[8]。

光譜角制圖(spectral angle mapper，SAM)算法是通過計(jì)算光譜之間夾角來確定二者之間整體相似性的算法，其突出特點(diǎn)是簡單高效，并可以在一定程度上克服照度和地形的影響[9]。但由于SAM算法對光譜曲線的局部微小特征差異刻畫不夠精細(xì)，當(dāng)出現(xiàn)礦物種類不同而光譜曲線特征相似的“異物同譜”現(xiàn)象時，區(qū)分效果較差，一些學(xué)者嘗試通過設(shè)置局部光譜區(qū)間權(quán)重[10-11]和選取波段組合[12]等方法改進(jìn)SAM算法。小波包分解是在小波分解的基礎(chǔ)上進(jìn)一步發(fā)展起來的一種更加完善的分解方法，能夠同時對信號的低頻和高頻部分進(jìn)行分解，并且這種分解既無冗余，也無疏漏[13]，因此能夠?qū)Π写罅恐?、高頻細(xì)節(jié)信息的遙感影像、生物醫(yī)學(xué)和語音等信號進(jìn)行更好的時頻局部化分析。

針對上述問題，本文提出了一種基于小波包分解和權(quán)重SAM的巖心高光譜蝕變信息提取方法，選取小波包分解系數(shù)中信息熵統(tǒng)計(jì)量構(gòu)建特征矢量，通過設(shè)置權(quán)重得到差異特征區(qū)間，并結(jié)合SAM算法進(jìn)行巖心高光譜蝕變信息提取研究，旨在提高巖心高光譜蝕變礦物的識別精度。

1 算法描述

SAM算法是利用圖像像元光譜和參考光譜之間矢量的夾角大小來表示光譜曲線間的相似程度，并通過選擇閾值對圖像未知像元進(jìn)行歸類，以達(dá)到識別地物的目的。其計(jì)算公式為[14]

(1)

式中：x和y分別代表圖像像元光譜和參考光譜；α為像元光譜與參考光譜二者之間的夾角，取值范圍為[0，π/2]。α的值越小，說明待分類的未知像元光譜與參考光譜的相似度越高，劃歸為同類的可能性也就越高[15-16]。但由于SAM算法要求地物的整體光譜均參與計(jì)算，致使非特征光譜及噪聲的比重會壓制診斷性特征光譜對分類的影響，為后續(xù)分類中選擇相似性閾值帶來難度。

權(quán)重SAM算法的思路是在圖像像元光譜與參考光譜差異較大的特征區(qū)間進(jìn)行權(quán)重設(shè)置，以突出局部特征差異。光譜曲線包含的波段總數(shù)為N，假設(shè)差異較大的特征區(qū)間包含有N1個波段，對應(yīng)的像元光譜和參考光譜反射率值分別為xav和yav，v∈[1,N1]，差異較小的普通區(qū)間包含有N2個波段，對應(yīng)的像元光譜和參考光譜的反射率值分別為xbu和ybu，u∈[1,N2]，且N1+N2=N，將差異較大的特征區(qū)間添加權(quán)重系數(shù)γ，γ≥1，當(dāng)γ=1時即為傳統(tǒng)SAM算法，計(jì)算公式(1)變化為

(2)

(3)

信號總能量E可以分解為所有不同小波包子分量的能量和，即

(4)

相對小波能量lp反映了每個子分量中小波包能量的分布情況，可由子分量能量與信號總能量之比求得，即

lp=Ep/E。

(5)

信息熵用來衡量各種概率分布所包含信息量的大小，可以作為一個系統(tǒng)復(fù)雜程度的度量。每個子分量小波包信息熵統(tǒng)計(jì)量定義為[17-18]

Wp=-lplnlp。

(6)

通過計(jì)算得到每個像元光譜曲線X所對應(yīng)的p個子分量小波包信息熵統(tǒng)計(jì)量，構(gòu)造p維信息熵特征矢量T，即

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

(7)

不同類別礦物的診斷性光譜特征不同，反映在各個子分量上具有不同的能量分布特性，使得各個子分量的小波包信息熵統(tǒng)計(jì)量也有所差異。因此，可以通過構(gòu)造信息熵特征矢量來描述不同礦物的光譜曲線特征，并找到圖像像元光譜與參考光譜信息熵矢量差異較大的特征區(qū)間，對其設(shè)置權(quán)重，采用權(quán)重SAM算法分別計(jì)算參考光譜和圖像像元光譜小波包信息熵特征矢量的夾角，并根據(jù)夾角余弦值的大小進(jìn)行蝕變礦物匹配，達(dá)到信息提取的目的。

2 數(shù)據(jù)與實(shí)驗(yàn)

2.1 數(shù)據(jù)獲取

實(shí)驗(yàn)數(shù)據(jù)是利用HySpex巖心地面成像高光譜測量系統(tǒng)對我國南方某火山巖型鈾礦區(qū)鉆孔巖心樣本掃描得到的短波紅外(SWIR)波段數(shù)據(jù)，礦區(qū)熱液活動期形成了堿性和酸性蝕變相互疊加的礦化和對稱中心式近礦圍巖蝕變。早期堿交代礦化蝕變主要包括鈉長石化、綠泥石化和赤鐵礦化等，晚期酸性礦化蝕變主要表現(xiàn)為螢石化和水云母化(伊利石或伊蒙混層黏土礦物)等，鉆孔巖心樣本巖性總體表現(xiàn)為碎斑流紋巖，斑晶主要為鉀長石、斜長石、石英和黑云母。HySpex成像光譜儀包含可見光—近紅外(VNIR)和SWIR波段2個傳感器，能夠同時采集獲取0.4～2.5 μm光譜范圍內(nèi)的地面成像高光譜數(shù)據(jù)，其中SWIR傳感器主要技術(shù)指標(biāo)見表1。

表1 HySpex SWIR傳感器技術(shù)指標(biāo)Tab.1 Technical specifications of HySpex SWIR sensor

2.2 數(shù)據(jù)預(yù)處理

通過地面掃描得到的鉆孔巖心高光譜數(shù)據(jù)，在采集過程中受儀器本身因素和外部環(huán)境影響，原始高光譜圖像數(shù)據(jù)除含有巖心樣本自身有效信息以外，還包含了噪聲干擾信息[19]。因此，為了減少噪聲對蝕變信息提取的干擾，需要對原始高光譜數(shù)據(jù)進(jìn)行輻射校正、反射率轉(zhuǎn)換和小波去噪等預(yù)處理。

具體操作為：采用系統(tǒng)自帶軟件，依據(jù)輻射定標(biāo)參數(shù)完成輻射校正，得到輻射亮度數(shù)據(jù)；為了利用反射光譜特征提取礦化蝕變信息，采用經(jīng)驗(yàn)線性法對輻射亮度值進(jìn)行反射率轉(zhuǎn)換；干擾信息在反射率光譜曲線上多表現(xiàn)為小的鋸齒噪聲，因此在不影響光譜曲線整體形態(tài)和診斷特征位置的前提下，運(yùn)用小波變換方法對反射率數(shù)據(jù)進(jìn)行光譜去噪，減弱或消除明顯的鋸齒噪聲影響，提高信息提取的精度。

2.3 蝕變礦物端元光譜提取

端元光譜提取是高光譜遙感影像分類和目標(biāo)識別的基礎(chǔ)。以往研究中，礦物端元光譜的確定通常有2種方式：①使用光譜儀測量樣品或者直接從礦物標(biāo)準(zhǔn)光譜庫中獲取；②直接通過數(shù)學(xué)方法從高光譜影像上選擇端元，提取礦物的端元光譜[20]。本文采用第2種方法，采用ENVI軟件最小噪聲分離(minimum noise fraction，MNF)變換對高光譜數(shù)據(jù)進(jìn)行噪聲白化和降維處理，利用純凈像元指數(shù)(pixel purity index，PPI)方法，提取高光譜圖像上的相對純像元，將純像元投影到N維可視化空間進(jìn)行N維散度分析得到各端元光譜曲線，再結(jié)合標(biāo)準(zhǔn)礦物光譜庫及專家知識人工確定最終的礦物端元光譜曲線?；贖ySpex巖心高光譜數(shù)據(jù)提取并確定的端元礦物主要包括高嶺石、綠泥石、地開石和伊利石(圖1)。

圖1 蝕變礦物端元光譜曲線Fig.1 Endmember spectrum curves of altered minerals

2.4 小波基選取

不同的小波基函數(shù)性質(zhì)各異，對光譜信號的變換結(jié)果影響很大。使用不同的小波基函數(shù)進(jìn)行小波包分解會產(chǎn)生不同的分解系數(shù)，進(jìn)而影響所構(gòu)建的小波包信息熵矢量及最終分類效果。而不同的分解層數(shù)同樣影響著小波包信息熵特征矢量表征原始光譜曲線特征的好壞。因此，選取適用于鉆孔巖心高光譜數(shù)據(jù)處理的最優(yōu)小波基和最佳分解層數(shù)顯得尤為重要。本研究中通過計(jì)算不同小波基函數(shù)和分解層數(shù)與礦物光譜曲線信息熵矢量夾角的關(guān)系來確定最佳小波基和分解層數(shù)。夾角余弦值越大，不同礦物類別間的可區(qū)分性越高，信息熵矢量更有利于刻畫光譜曲線特征。通過對比最終確定daubechies4作為最優(yōu)小波基，最佳分解層數(shù)為8層。

2.5 差異特征區(qū)間選取

利用權(quán)重SAM算法識別某種礦物時，可找到與其他礦物差異較大的特征區(qū)間，并對其設(shè)置權(quán)重系數(shù)γ，增加不同礦物之間的可區(qū)分性。在信息熵特征空間中，通過比較所構(gòu)建的4種礦物端元光譜小波包信息熵矢量，采用差異較大的特征區(qū)間為小波包分解的前8個子分量信息熵(圖2)。

圖2 4種礦物部分子分量信息熵Fig.2 Partial subcomponent entropy of four kinds of minerals

從圖2可以發(fā)現(xiàn)，4種礦物的小波包信息熵在此區(qū)間各個子分量上差異較大，表現(xiàn)出不同的分布特性，通過試驗(yàn)權(quán)重系數(shù)取γ=2。

3 結(jié)果與討論

利用Matlab語言將本文算法編寫程序?qū)崿F(xiàn)，對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理。通過對巖心高光譜影像上提取的高嶺石、綠泥石、地開石和伊利石4種礦物端元光譜曲線進(jìn)行8層小波包分解，根據(jù)各自分解系數(shù)生成相對應(yīng)的小波包信息熵特征矢量。依據(jù)SAM算法原理分別計(jì)算4種礦物端元光譜曲線兩兩間光譜矢量和小波包信息熵矢量夾角余弦值大小，結(jié)果如表2所示。通過對比發(fā)現(xiàn)，原始光譜曲線光譜矢量夾角余弦值均大于基于小波包分解系數(shù)生成的信息熵矢量夾角余弦值，經(jīng)過小波包分解后的信息熵矢量增大了不同礦物類別間的夾角，提高了可區(qū)分度，更有利于刻畫光譜曲線的特征。

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

圖3為利用巖心高光譜數(shù)據(jù)分別進(jìn)行傳統(tǒng)SAM算法(圖3(b))以及本文方法(圖3(c))的蝕變信息提取結(jié)果。

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

(b)傳統(tǒng)SAM算法提取結(jié)果

(c)本文方法提取結(jié)果

圖3 HySpex鉆孔巖心高光譜數(shù)據(jù)蝕變信息提取結(jié)果Fig.3 Extraction results of alteration information using HySpex drill core hyperspectral data

從巖心圖像提取結(jié)果可以看出，本文方法提取的蝕變礦物比傳統(tǒng)SAM算法提取的范圍更為精確細(xì)致。通過整個鉆孔巖心提取結(jié)果分析發(fā)現(xiàn)，高嶺石、地開石和明礬石等酸性蝕變礦物在鉆孔上部含量多，下部含量少，而伊利石、綠泥石和碳酸鹽巖等中、堿性蝕變礦物上部含量少，下部尤其是深部礦化段含量明顯增多，總體分布符合本地區(qū)“上酸下堿”的蝕變特征。

為了定量評價(jià)2種方法提取結(jié)果的有效性，以鉆孔柱狀圖、目視判別及像元光譜曲線為依據(jù)，在圖像中隨機(jī)選取5 000個像素點(diǎn)(包括背景值)作為真值進(jìn)行檢驗(yàn)，構(gòu)建混淆矩陣，選取總體分類精度和Kappa系數(shù)指標(biāo)進(jìn)行精度評價(jià)，統(tǒng)計(jì)結(jié)果如表3所示。

表3 不同匹配方法分類精度統(tǒng)計(jì)Tab.3 Classification accuracy statistics of different matching methods

從表3可以發(fā)現(xiàn)，與傳統(tǒng)SAM算法相比，本文方法的分類精度更高些，總體分類精度由71.51%提高到75.33%，Kappa系數(shù)由0.667 8增加到0.706 3。

通過對本文涉及的蝕變礦物及其光譜曲線分析發(fā)現(xiàn)，高嶺石和地開石多伴生發(fā)育，光譜曲線特征接近，使得依據(jù)光譜相似性度量準(zhǔn)則的傳統(tǒng)SAM算法對其區(qū)分效果有限，而經(jīng)過小波包變換后得到的信息熵矢量夾角被增大，同時通過設(shè)置差異特征區(qū)間，突出了局部特征信息，使得二者之間可區(qū)分性被進(jìn)一步增強(qiáng)，分類精度得到了提高。

實(shí)驗(yàn)結(jié)果表明，本文提出的基于小波包變換和權(quán)重SAM算法的巖心高光譜蝕變信息提取方法，定義的信息熵特征矢量可以較好地表征原始光譜信息，結(jié)合差異特征區(qū)間的權(quán)重設(shè)置，增大了不同礦物類別間的可區(qū)分性，分類精度優(yōu)于傳統(tǒng)SAM算法，具有較好的適用性。

4 結(jié)論

本文運(yùn)用小波包變換和權(quán)重SAM相結(jié)合的方法對HySpex巖心高光譜數(shù)據(jù)進(jìn)行蝕變信息提取，并與傳統(tǒng)SAM算法進(jìn)行對比實(shí)驗(yàn)。主要結(jié)論如下：

1)本文方法提高了蝕變礦物提取的總體精度，由傳統(tǒng)SAM算法的71.51%提高到75.33%，具有較好的適用性。

2)依據(jù)小波包分解系數(shù)計(jì)算得到的信息熵特征矢量可有效表征不同礦物光譜差異，結(jié)合差異特征區(qū)間的權(quán)重設(shè)置，增大了不同類別礦物間的可區(qū)分性。

3)提取得到的鉆孔巖心礦化蝕變信息，可展現(xiàn)蝕變礦物在巖心整體垂向的發(fā)育狀況，為進(jìn)一步劃分酸、堿蝕變分帶，分析鈾礦富集規(guī)律和預(yù)測深部礦化前景提供依據(jù)。

在后續(xù)的研究中擬對距離和相關(guān)系數(shù)等不同的匹配方法進(jìn)行對比研究，同時對數(shù)據(jù)預(yù)處理中不同的光譜去噪方法進(jìn)行試驗(yàn)，進(jìn)一步提高信息提取精度。