基于光谱相关角和光谱信息散度的高光谱蚀变信息提取

吴 浩,徐 元 进,高 冉

(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