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

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

(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)提取得到的鉆孔巖心礦化蝕變信息，可展現蝕變礦物在巖心整體垂向的發育狀況，為進一步劃分酸、堿蝕變分帶，分析鈾礦富集規律和預測深部礦化前景提供依據。

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