基于光譜相關角和光譜信息散度的高光譜蝕變信息提取

吳 浩,徐 元 進,高 冉

(1.中國地質大學數學地質遙感地質研究所，湖北 武漢 430074；2.杭州科瀾信息技術有限公司，浙江 杭州 310000)

基于光譜相關角和光譜信息散度的高光譜蝕變信息提取

吳 浩1,徐 元 進1,高 冉2

(1.中國地質大學數學地質遙感地質研究所，湖北 武漢 430074；2.杭州科瀾信息技術有限公司，浙江 杭州 310000)

針對高光譜遙感蝕變信息提取過程中，由于混合像元的不可避免，導致蝕變礦物光譜曲線存在不同程度的失真而影響目標礦物識別精度的問題，提出一種基于光譜相關角(Spectral Correlation Angle,SCA)和光譜信息散度(Spectral Information Divergence,SID)的高光譜遙感蝕變信息提取算法(SIDSCAtan)。利用植被覆蓋度表征植被信息在混合像元中所占百分比，劃分出6種植被失真類型。采用不同區分方法分別比較失真光譜與理想光譜的差異，實驗表明，當輸入的光譜信息具有微小差異時，方法SIDSCAtan能夠做出較大的響應，在識別光譜整體形態的前提下，增強了對光譜局部特征差異的區分能力。以云南省個舊西區為研究區，運用該方法提取區內蝕變信息，應用效果較好。

高光譜遙感；蝕變信息提取；光譜相關角；光譜信息散度

0 引言

在礦產勘查中，蝕變巖作為一種重要的找礦標志已有數百年的歷史，其中熱液型礦床常伴有不同程度的近礦圍巖蝕變現象[1,2]。遙感找礦異常提取的目的是探測和圈定與成礦密切相關的蝕變礦物、蝕變巖帶和礦體表生氧化帶[3]。高光譜遙感地表巖性識別的關鍵是目標礦物光譜特征的識別，利用相似度衡量待匹配光譜特征差異以達到識別地物的目的，已得到廣泛的應用[4-7]。

目前計算相似度的方法有很多，傳統的方法主要包括光譜角(SAM)、光譜相關角(SCA)、波譜特征擬合(SFF)、光譜信息散度(SID)、混合調制匹配濾波等[8-15]。由于光譜角是通過計算光譜矢量之間的夾角來衡量其相似性，夾角的大小只與光譜矢量方向有關，與其輻射亮度無關，因此當兩種礦物的光譜矢量方向相似而輻射亮度大小有差別時，區分效果較差，并且光譜角只考慮夾角的絕對值，而不能識別待匹配光譜的正負相關性。光譜相關角則有效避免了光譜負相關，并保持了光譜角最小化陰影影響的優勢。Du等提出了光譜角和光譜信息散度的組合方法(SIDSAMtan)[16-18]，對比于單一匹配算法，提高了礦物識別能力，在高光譜純凈端元提取中應用廣泛[19-20]。

本文結合已有研究成果，提出基于光譜相關角和光譜信息散度(SIDSCAtan)的高光譜蝕變信息提取方法。實驗表明，當輸入的光譜信息具有微小差異時，該方法能夠做出較大的響應，同時通過在云南省個舊西區實際蝕變信息提取中應用該方法，取得較好的效果，驗證了基于光譜相關角和光譜信息散度(SIDSCAtan)的高光譜蝕變信息提取方法的可靠性。

1 原理及方法

兩光譜信號Si、Sj可表示為L維向量(si1,si2,……,siL)T和(sj1,sj2,……,sjL)T的形式。基于SID-SCA的混合方法在考慮光譜整體形態和局部特征差異的基礎上，相對于傳統巖性識別方法，期望可以提高相似光譜的區分能力，從而達到精確提取蝕變信息的目的。

SID-SCA模型計算光譜Si、Sj相似性定義為：

SIDSCAtan(si,sj)=(ΨSID(si,sj))(tan(ΨSCA(si,sj)))

(1)

或SIDSCAsin(si,sj)=(ΨSID(si,sj))(sin(ΨSCA(si,sj)))

(2)

在巖性光譜識別過程中，主要利用光譜Si、Sj的垂向距離而不是其投影距離，因此這里計算光譜相關角的正切(1)和正弦(2)，而不考慮其余弦。

光譜相關角ΨSCA(si,sj)定義如下：

(3)

(4)

式中：L為遙感數據的波段數；i和j分別代表進行匹配的兩個光譜信號；m為光譜信號中的某一波段。

式(4)的結果介于-1和1之間，反映了Si、Sj的線性相關程度。因此光譜相關角Ψ的返回值介于0和1.570796 rad之間，度量兩光譜信號的相似度。

光譜信息散度ΨSID(si,sj)定義如下：

ΨSID(si,sj)=D(si‖sj)+D(sj‖si)

(5)

Si、Sj的概率向量分別是p=(p1,p2,……,pL)和q=(q1,q2,……,qL)，其中:

由信息理論可分別得到光譜Si、Sj中波段k的自信息為Ik(si)=-log2(pk)和Ik(sj)=-log2(qk)。因此可得Sj關于Si的相對熵以及D(si‖sj)和Si關于Sj的相對熵D(sj‖si)，如式(6)、式(7)：

(6)

(7)

2 方法分析與比較

本研究使用的高光譜數據來源于美國E0-1衛星搭載的Hyperion高光譜成像儀，為L1R產品。該產品共有242個波段，經過輻射定標的有198個(VNIR8-57，SWIR77-224)，其中VNIR56、57，SWIR77、78重疊，同時部分波段受水汽影響嚴重，所以只使用155個波段(426.816～2 353.180nm)。在使用之前，對影像進行絕對輻射值轉換、壞線修復、條紋去除、大氣校正、幾何校正處理(圖1)。

植被覆蓋度(Fc)為單位面積內植被地上部分垂直投影面積所占百分比(圖2，見封2)，在高光譜影像中，(1-Fc)可表示為單個植被像元內摻雜其他干擾信息(土壤、巖石、道路等)百分比[21]。為了對其他學者提出的方法與本文基于SID-SCA方法的信息提取精度進行比較，本文從高光譜影像中選擇摻雜不同比例干擾信息的植被像元光譜(圖3,見封2)，利用不同區分方法比較失真光譜與理想光譜(圖4)的差別。

本文采用像元二分模型，利用NDVI估算植被覆蓋度。假設高光譜遙感影像一個像元只由植被和非植被(裸土)組成，所占面積為Fc和(1-Fc)，光譜信息分別為Sveg和Ssoil，則混合像元光譜信息定義為：

圖1 高光譜遙感影像

Fig.1 Hyperspectral RS image

圖4 USGS光譜庫植被光譜曲線

Fig.4 Reflectance spectra of vegetation in USGS library spectra

S=Sveg×Fc+Ssoil×(1-Fc)

(8)

通常每個像元的NDVI可以看成是植被部分NDVIveg和非植被部分NDVIsoil的加權平均，那么混合像元植被指數定義如式(9)：

NDVI=NDVIveg×Fc+NDVIsoil×(1-Fc)

(9)

則植被覆蓋度(Fc)如式(10)：

Fc=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(10)

如圖2，將植被覆蓋度信息由高到低劃分為6個等級，分別表征植被像元中的干擾信息摻雜比由小到大，對應圖3中6條光譜曲線，比較圖4可知，當干擾信息摻雜比逐漸增大，光譜失真愈加嚴重。

相似度的模值體現植被像元間的差異性，模值越大，表示其差異性越大，則更易被區分。由式(11)可知，由于0<θ(si，sj)<1，顯然SIDSCAtan的相似度模值將大于SIDSCAsin。因此，下面只分析比較SIDSCAtan、SIDSAMtan和SIDSGAtan的區分能力，如表1。

(11)

表1 光譜識別方法(SIDSCAtan、SIDSAMtan、SIDSGAtan)相似度計算結果

Table 1 Results of the similarity calculated by different methods

混合植被像元類型方法SIDSCAtanSIDSAMtanSIDSGAtan10.01030.00550.014320.01070.00550.014930.01140.00570.015040.01450.00710.017550.02200.01070.022960.03370.01620.0311

為了更直觀地對比和評價上述方法的區分能力，每種區分方法以最小失真類型1為基準，比較其他類型的相似度對于類型1的增量，如圖5a。由圖5a可知，隨著植被像元信息失真程度的增大，三種方法的區分能力也呈上升趨勢，并且SIDSCAtan方法的相似度增幅明顯高于SIDSAMtan和SIDSGAtan。為了測試三種方法區分能力的增速，分別計算其相似度梯度(圖5b)。由圖5b可知，由類型1至類型2，方法SIDSAMtan的相似度梯度增速較快，SIDSCAtan高于SIDSGAtan，至類型3，SIDSAMtan增速降為最小。隨后，方法SIDSCAtan的區分能力一直保持最高增速。

圖5 不同方法計算相似度曲線

Fig.5 Maps of the similarity calculated by different methods

綜上可知，SIDSCAtan區分方法對植被光譜信息變化的敏感性更強，即對目標地物特征吸收譜帶的微小失真有較大的響應，表現出較強的地物識別能力。由于混合像元的影響，蝕變礦物光譜在高光譜影像中也存在不同程度的失真，特征譜帶發生變化。因此本文試圖將SIDSCAtan區分方法應用于蝕變信息的提取，期望得到理想的結果。

3 實例應用與結果分析

3.1 研究區

本文研究區選在位于中國云南省個舊市西區(102°57′36.51″～103°2′51.78″E，23°28′8.56″～23°33′6.98″N),區內圍巖蝕變類型較多，主要類型有矽卡巖化、硅化、綠泥石化、絹云母化、云英巖化、鈉長石化，其次有鉀長石化、電氣石化、黃鐵礦化、碳酸鹽化、鐵錳礦化、赤鐵礦化、褐鐵礦化等，且不同巖石類型中蝕變類型及強度均有明顯差異。

3.2 蝕變信息提取及結果分析

蝕變礦物在426.816～2 353.180 nm光譜區間內具有明顯的診斷波譜特征(圖6)，是直接識別礦物類型與組分的關鍵。本研究使用圖6中的礦物光譜為參考光譜，利用SIDSCAtan方法提取個舊西區“鐵化”和“泥化”蝕變信息，如圖7、圖8所示。

圖6 USGS光譜庫中礦物光譜曲線

Fig.6 Reflectance spectra of mineral in USGS library spectra

圖7 鐵化礦物蝕變信息

Fig.7 Iron alteration minerals maps

圖8 粘土礦物蝕變信息

Fig.8 Clay alteration minerals maps

個舊礦區圍巖蝕變作用主要是花崗巖侵入及與之有關的成礦作用引起，因此，蝕變信息在空間上依附于花崗巖侵入巖體呈有規律的帶狀分布。由圖7可知，鐵的氧化物信息在位置和輪廓上大致相同，多沿北西向和北東向大斷裂呈線狀分布，野外地質資料證實斷裂帶內均出現帶狀分布的褐鐵礦化石英脈(圖9a)。經風化淋濾作用后地表多呈褐紅色，原蝕變礦物為黃鐵礦，出露地表易氧化為赤鐵礦、針鐵礦、褐鐵礦等鐵的氧化物。如圖7a、圖7d，赤鐵礦、黃鉀鐵礬的一些零星點狀異常，大多數分布于斷裂帶、花崗巖侵入體與地層接觸帶、不同巖性突變邊界，如圖9b，在北西向斷裂帶內發育強烈的鐵化蝕變信息。泥化蝕變礦物易受流水動力作用而沿水系富集，如圖8b、圖8c，絹云母和高嶺土信息在研究區西南角呈條帶狀分布，野外調查顯示此處為干涸山溝，溝內高嶺土、葉臘石、絹云母、蛇紋石較為富集，含少量赤鐵礦及黃鉀鐵礬，可知此處為真實蝕變信息，但屬找礦假異常。如圖8a、圖8d，綠泥石和伊利石信息呈團塊狀，并與鐵化異常輪廓相似，出現于重要的地質體邊界(斷裂、侵入邊界、地層界線)，如圖9c，斷層裂隙中充填綠泥石(圖9d)，花崗巖侵入體內有強烈泥化蝕變信息。綜上可知，鐵化和泥化信息大多數為真實蝕變信息。

圖9 野外驗證照片

Fig.9 The photos of field validation

4 結論

本文提出一種基于光譜相關角和光譜信息散度的高光譜遙感蝕變信息提取算法，在識別光譜整體形態的前提下，增強了對光譜局部特征差異的區分能力，由對混合植被像元的識別實驗可知，SIDSCAtan方法對光譜的細微差異有較大的響應。利用SIDSCAtan方法提取了云南個舊市西區“鐵化”和“泥化”蝕變信息，鐵的氧化物信息分布情況較為相似，野外驗證發現大多數分布于斷裂帶、花崗巖侵入體與地層接觸帶、不同巖性突變邊界，為真實蝕變信息；泥化信息雖受流水動力作用易沿水系富集，為找礦假信息，但由野外驗證可知，此信息為真實蝕變信息。綜上可知，SIDSCAtan方法提取的蝕變信息具有較高的準確性，在熱液型礦床勘查中具有重要的指示意義。

本文在利用SIDSCAtan提取蝕變信息時，相似系數閾值是通過經驗值確定的，具有一定的盲目性。因此，如何更合理地確定SIDSCAtan閾值，是今后有待解決的重點問題。本文結合野外地質調查工作，對SIDSCAtan方法在個舊西區的蝕變信息提取效果做了定性評價，如何定量評價信息提取精度是今后進一步研究的重點內容。

[1] MURPHY R J,MONTEIRO S T.Mapping the distribution of ferric iron minerals on a vertical mine face using derivative analysis of hyperspectral imagery(430-97002nm)[J].Isprs Journal of Photogrammetry & Remote Sensing,2013,75(328):29-39.

[2] 王潤生,熊盛青,聶洪峰，等.遙感地質勘查技術與應用研究[J].地質學報,2011,85(11):1699-1743.

[3] BAISSA R,LABBASSI K,LAUNEAU P,et al.Using HySpex SWIR-320m hyperspectral data for the identification and mapping of minerals in hand specimens of carbonate rocks from the Ankloute formation[J].Journal of African Earth Sciences,2011,61(1):1-9.

[4] POUR A B,HASHIM M,GENDEREN J V.Detection of hydrothermal alteration zones in a tropical region using satellite remote sensing data:Bau goldfield,Sarawak,Malaysia[J].Ore Geology Reviews,2013,54(8):181-196.

[5] ENTON B.Mineral mapping in the Kap Simpson complex,central East Greenland,using HyMap and ASTER remote sensing data[J].Advances in Space Research,2011,47(1):60-73.

[6] 王潤生.遙感地質技術發展的戰略思考[J].國土資源遙感,2008(1):1-12.

[7] 代晶晶,王瑞江,王潤生,等.基于蝕變信息提取的西藏班公湖-怒江成礦帶中段斑巖銅礦找礦預測[J].地球學報,2012,33(5):755-762.

[8] XU Y,MA H,PENG S.Study on identification of altered rock in hyperspectral imagery using spectrum of field object[J].Ore Geology Reviews,2014,56:584-595.

[9] ZHANG X Y,LI P J.Lithological mapping from hyperspectral data by improved use of spectral angle mapper[J].International Journal of Applied Earth Observation & Geoinformation,2014,31(5):95-109.

[10] MAHDIEH H Z,TANGESTANI M H,ROLDAN F V,et al.Sub-pixel mineral mapping of a porphyry copper belt using EO-1 Hyperion data[J].Advances in Space Research,2014,53(3):440-451.

[11] CHEN X,WARNER T A,CAMPAGNA D J.Integrating visible,near-infrared and short-wave infrared hyperspectral and multispectral thermal imagery for geological mapping at Cuprite,Nevada[J].International Journal of Remote Sensing,2007,110(3):344-356.

[12] 張修寶,袁艷,景娟娟,等.信息散度與梯度角正切相結合的光譜區分方法[J].光譜學與光譜分析,2011,31(3):853-857.

[13] DU B,ZHANG Y,ZHANG L,et al.A hypothesis independent subpixel target detector for hyperspectral images[J].Signal Processing,2015,110:244-249.

[14] 閻繼寧,周可法,王金林,等.基于SAM與SVM的高光譜遙感蝕變信息提取[J].計算機工程與應用,2013,49(19):141-146.

[15] 劉萬軍,楊秀紅,曲海成,等.基于光譜信息散度與光譜角匹配的高光譜解混算法[J].計算機應用,2015,35(3):844-848.

[16] DU Y,CHANG C I,REN H,et al.A new hyperspectral discrimination measure for spectral similarity[C].Proceeding of the SPIE,2003.430-439.

[17] DU Y,CHANG C I,REN H,et al.A new hyperspectral discrimination measure for spectral characterization[J].Optical Engineering,2004,43(8):1777-1786.

[18] DU Y,IVES R,ETTER D,et al.A one-dimensional approach for iris identification[C].Proceeding of the SPIE,2004.237-247.

[19] CHANG C I.An information-theoretic approach to spectral variability,similarity,and discrimination for hyperspectral image analysis[J].IEEE Transactions on Information Theory,2000,46(5):1927-1932.

[20] NARESH K M,SESHASAI M,VARA PRASAD V R,et al.A new hybrid spectral similarity measure for discrimination among Vigna species[J].International Journal of Remote Sensing,2011,32(14):4041-4053.

[21] 李曉松,李增元,高志海,等.基于NDVI與偏最小二乘回歸的荒漠化地區植被覆蓋度高光譜遙感估測[J].中國沙漠,2011,31(1):162-167.

Extraction of Alteration Information from Hyperspectral Imagery Based on SCA and SID

WU Hao1,XU Yuan-jin1,GAO Ran2

(1.Institute of Mathematical and Remote Sensing Geology,China University of Geosciences ,Wuhan 430074;2.Hangzhou Kelan Information Technology Co.Ltd,Hangzhou 310000,China)

Hyperspectral remote sensing has been widely used in alteration information extraction.However,the spectra of altered minerals have varying degrees of distortion because of mixed pixel,which affects the accuracy of target minerals identification.A new approach to extract the alteration information from hyperspectral imagery based on SCA and SID has been presented.Using vegetation coverage to characterize the percentage of vegetation information in mixed pixel to divide the vegetation of distortion into six types.The differences between the distortion spectra and the true spectrum is distinguished by different spectral discrimination methods.The experiment shows that when the input spectral information with a slight difference,SIDSCAtancould make a big response,not only identify the whole spectrum shapes,but also enhance the partial differences of the spectral charateristics.The authors selected the west zone in Gejiu County of Yunnan Province as the study area,extracted the information of altered minerals from the hyperspectral image of the study area,the results obtained bySIDSCAtanhad higher reliability.

hyperspectral remote sensing;alteration information extraction;spectral correlation angle;spectral information divergence

2015-07-17；

2015-10-16

國家自然科學基金項目(41072246);中央高校基本科研業務費專項資金資助項目(CUG120116)

吳浩(1989-)，男，碩士研究生，主要研究方向為遙感地質、國土資源遙感。E-mail:wuhaowzd@163.com

10.3969/j.issn.1672-0504.2016.01.009

P627

A

1672-0504(2016)01-0044-05