张 媛,张杰林,赵学胜,袁 博
(1.中国矿业大学(北京)地球科学与测绘工程学院,北京 100083;2.核工业北京地质研究院遥感重点实验室,北京 100029)
高光谱遥感将传统图像的空间维与光谱维信息融合为一体,其高光谱分辨率的特点使得本来在宽波段多光谱遥感图像中不可探测的地物在高光谱遥感图像中能够被探测和区分[1]。钻探是地质勘查的重要技术手段之一,岩心是钻探获取的重要成果。传统的岩心矿物识别主要是通过人工方法判别,这种方法对人的经验依赖性较高,容易出现错判或漏判等现象,且需耗费大量人力和时间。矿化蚀变信息是找矿的重要标志,钻孔岩心作为铀矿勘探成果的主要载体记录着地质体垂向变化信息[2],通过对其中与铀成矿有关的蚀变信息的提取,可对缩小找矿范围和提高找矿准确度发挥重要作用。利用高光谱遥感技术对岩心进行测量,获取岩心的光谱数据,并利用其光谱连续性强、分辨率高和图谱合一等特点[3]进行岩心光谱数据分析和岩心矿化蚀变信息提取,是目前矿物识别与填图及岩心编录的重要发展方向之一。
光谱匹配是一种高光谱矿化蚀变信息提取的重要方法。光谱角填图(spectral angle mapping,SAM)法[4-5]是光谱匹配中最重要的方法之一,其简单高效、标量乘不变的特点使其得到了广泛的应用,是目前岩心矿化蚀变信息提取中常用的方法。但由于使用SAM方法时岩心的全部光谱都参与运算[6],导致大部分次要光谱的相似性会压制少数特征光谱之间的差异,给后续处理中的阈值设定带来困难[7]。另外,由于岩心光谱的同质多象和类质同象等现象[8](主要表现在吸收谱带漂移、吸收深度以及波形宽度的不同),使得光谱特征匹配存在很多不确定性,导致SAM的实际应用效果有限。对于矿物种类不同而光谱特征相似的“异物同谱”现象,为了加大光谱相似度之间的差异、更好地区分不同矿物,一些学者通过采取在差异性较大的波谱区间设置权重[9]、分组计算和归一化计算[10]等方法改进SAM,提高了矿物识别和区分精度。目前SAM仍是大规模岩心矿化蚀变信息提取的重要方法[11-12]。在采用该方法提取岩心矿化蚀变信息时,一般都希望将属于同一蚀变种类而光谱特征不同的蚀变矿物分为一类,利用其诊断性光谱吸收特征来判断不同类别的蚀变矿物。但实际操作中,有时会将本属同一类的蚀变矿物判断为不同类别的蚀变矿物而造成漏分、错分现象;尽管可以同时取多条波谱曲线作为参考端元来改善这种情况,但端元数量越多,误差累积越大,对阈值设置的要求越高,且在短时间内很难找全代表该类蚀变矿物的全部光谱曲线。
针对上述问题,本文采取在同类蚀变矿物光谱曲线中寻找特征区间、并对诊断性吸收峰位置设置权重的方法,对传统的SAM法进行改进,旨在提高典型矿化蚀变矿物的高光谱异常信息识别精度及提取效率。
传统的SAM法是通过计算光谱矢量之间的夹角来判断参考波谱与像元波谱之间的相似性。设参考波谱为x,图像中的像元波谱为y,则它们之间的光谱角为
夹角α越小,说明越相似。参考光谱可以是实验室光谱或野外测定光谱,也可以是从图像中提取的像元光谱[13]。光谱角的大小只与光谱矢量的方向有关,与辐射亮度无关,从而减弱了照度和地形的影响[6]。然而,SAM是基于图像整体光谱特征的匹配方法,即对图像所有波段的信息同等处理;对于吸收峰相同而整体曲线有较大差异的同质光谱,则会增加光谱角计算的误差,阈值较小时会出现漏判,阈值较大时会将其他类别错分到此类。尽管通过选取同一种蚀变的多条典型光谱曲线作为参考端元可以改善这一缺陷,但这要求对参考端元的提取要全面而准确,需要用足够的时间和精力去寻找更好的方法来提取有“同类异谱”现象的参考端元。通常,从不同深度和不同空间位置取得的岩心中,同一种矿化蚀变的波谱形态都不尽相同,即这一段岩心中的参考端元未必适应于下一段岩心。因此,采用遍历所有岩心图像去全面、准确、快捷地寻找参考端元的方法缺乏可操作性而不可取。同时,由于端元波谱识别的过程是人工判别的过程,其数量的增多也会提升误差出现的概率。
笔者通过对图像像元进行分析发现:同类蚀变矿物的光谱曲线在某一波段区间的波形和吸收峰特征都十分相似;而在其他区间的差别则较大(图1),这是造成SAM在同类异谱现象时产生较大误差的主要原因,故将波形和吸收峰特征都十分相似的区间称为特征区间。
图1 同类异谱蚀变矿物典型波谱曲线Fig.1 Typical spectrum curves of altered minerals with same kind but different spectrum
为了减小非特征区间光谱波形差异对光谱角带来的影响,可对每类蚀变矿物取特征区间[b1,b2](b1,b2分别代表特征区间的起、止波长),仅对像元在该特征区间内的光谱反射率进行处理(即局部光谱角分类)。由于每类蚀变的诊断性光谱吸收峰均包含在该类蚀变矿物的特征区间内,因而特征区间光谱也是区别不同类蚀变的重要局部信息。
因同一蚀变矿物的结构和组分不同,特征区间内吸收峰的形状和宽度亦存在微小差异,且不同岩心图像数据中同类蚀变矿物的诊断性光谱吸收峰的位置是确定的(有的可能会有1个像元内的微弱漂移),而不同物质的峰值位置则明显不同。选取特征区间内的像元光谱进行的局部光谱角分类可以增强该区间波谱范围内的光谱特征,凸显光谱的吸收峰特性,拉大不同矿物在该区间的差异。考虑到峰值可能的漂移,本文在特征区间增加了光谱吸收峰位置附近同一波段范围内光谱反射率的权重ω。对于取自遥感图像的参考端元和图像中的每一个像元,设反射率数组为 r=(x1,x2,…,xn),变换后反射率数组为r'=(y1,y2,…,yn)(其中n为波段总数);将有m个特征峰值的波段位置组成的数组表示为P=(a1,a2,…,am),其中aj表示第j个波峰的波段位置。i为第i个波段,在特征光谱区间[b1,b2]范围内,对于每一个i∈[b1,b2],j∈[1,m]有
若i≠aj且i≠aj-1 且i≠aj+1,
则yi=xi;
否则
式中:ω为权重系数,应根据人工经验进行选择,ω=1时即为传统的 SAM。由此再对特征区间[b1,b2]内的参考端元和图像像元加权后的反射率数值计算光谱角,得到的结果可以拉大不同蚀变矿物在该光谱区间的差距,减小同类蚀变矿物光谱之间的差异,重点突出局部光谱信息;同时兼顾了峰值的位置,在阈值一定的情况下,能够从遥感图像中较好地提取蚀变矿物信息。此外,由于特征区间的缩小和吸收峰位置的强化,仅需1条代表该类型蚀变矿物的光谱曲线即可进行蚀变矿物信息提取,从而省去繁琐的端元提取和人工判别、筛选等过程。
本文采用的岩心高光谱图像是利用高光谱测量系统对研究区内的岩心标本扫描得到的短波红外(SWIR)数据,其技术指标及参数见表1。
表1 高光谱测量系统主要技术指标及参数Tab.1 Main technical performances of hyperspectral measurement system
为了更好地分析岩心高光谱数据的吸收特征和更准确地提取蚀变矿物信息,首先需要对截取的一小段岩心数据进行反射率转换和噪声去除。
本文采用统计学模型中的经验线性法(empirical line,EL)对岩心数据计算反射率。EL定标方法假设图像的DN值与反射率值R之间存在线性关系[14],即
式中:Gain为增益;Bias为偏移。在实际计算过程中,利用已测得的白板数据作为已知点的地面光谱,求得增益和偏移,进而计算岩心高光谱图像反射率。
噪声去除应尽量不改变光谱曲线在小范围内的形态,即不改变吸收峰的位置。对反射率图像中的每一个像元,利用其每奇数个波段求平均得到中间波段的反射率值,即
式中:b'i为处理后图像的反射率值;bi为处理前图像的反射率值;i为波段号;n为图像的波段总数;m可根据去噪的程度选择合适的值,噪声较小时m可选择1或2。“=bi”表示前m个波段和后m个波段的反射率值没有参与计算,仍用原反射率值(如m=2时,表示第1和2波段与倒数第1和2波段的值没有变化)。
该方法最大的优点是简单易行,既有效地保留了波段间的相关性,又消除了噪声对图像光谱在小范围内引起的剧变,提高了后续图像处理和蚀变矿物信息提取的精度。
本文取m=2的方法可能会造成峰值在一个波段范围内位置的改变(此情况出现的概率很小),但由于后续的光谱特征分析和识别、光谱端元提取以及对图像的匹配分类均是在应用该方法去噪的基础上进行的,因而不会影响图像匹配分类的精度。
蚀变矿物中不同离子(或基团)(如Fe2+,Fe3+,C和OH-等)的电子跃升、振动和转动[15]使其在遥感图像中具有诊断性光谱吸收特征,这些特征是识别蚀变矿物类别的基础。本文中的岩心数据所在矿区围岩蚀变发育,种类较多,与铀矿有关的围岩蚀变主要有赤铁矿化、伊利石化、绿泥石化、碳酸盐化、萤石化和碱交代等。综合分析研究区主要热液蚀变与铀矿化的物谱关系,以便为深部铀矿找矿提供基础信息[2]。为便于分析,截取部分岩心数据进行试验。岩矿鉴定结果显示,该岩心段主要的蚀变为伊利石化(Al-OH)、绿泥石化(Mg-OH)和少量的碳酸盐化(C)等3种类型,蚀变矿物几乎覆盖了整个岩心表面。
由于同种岩石或矿物因发育过程和发育状态不同,其成分、结构及光谱特征就会产生一定的差异;同一种矿化蚀变也会因区域地质背景不同而使其波谱曲线呈现一定差异(如诊断性光谱吸收特征峰的漂移)。因此,为提高光谱匹配的精度,笔者选择在图像中确定波谱端元。
本文的岩心高光谱数据中,伊利石化、绿泥石化和碳酸盐化矿物的诊断性光谱吸收特征[16]主要表现在短波红外(SWIR,1 300~2 500 nm)波段:①伊利石化矿物由于铝阳离子的电子跃迁以及Al-OH基团伸缩振动的合频与倍频作用,在短波红外光谱段产生光谱吸收特征[17],特别是在2 210 nm和2 354 nm处形成明显的吸收峰“二元结构”。在该岩心段的B80(R)B200(G)B48(B)假彩色合成图像中显示为亮白色;②绿泥石化矿物诊断性光谱吸收特征主要体现在2 258 nm(Mg-OH基团)和2 360 nm处,整体反射率较低,在该岩心段假彩色合成图像中呈暗绿色点状分布;③碳酸盐化矿物由于受C影响,在2 342 nm处形成强烈的吸收峰(且左宽右窄),据此特征可区别于其他矿物。除上述特征外,3种蚀变矿物均在1 417 nm和1 910 nm处存在水汽吸收峰,且同一种物质的吸收峰均可能存在一个波段范围内漂移。3种蚀变矿物在岩心高光谱图像中的诊断光谱标志如表2所示。
表2 3种蚀变矿物在图像中的诊断光谱特征Tab.2 Diagnostic spectrum features of three kinds of typical alteration minerals
根据伊利石化、绿泥石化及碳酸盐化矿物的光谱特征,从岩心高光谱图像中分别对3种蚀变矿物提取出多条端元光谱曲线(图1)。
通过分析图1中3种蚀变矿物的典型波谱曲线可以看出,在波段[182,256]区间内,伊利石化矿物在2 054~2499 nm波形相似且包括了所有的诊断性光谱吸收特征,吸收峰位置分别为2 210 nm,2 354 nm和2 445 nm;绿泥石化矿物在1 916 nm处的吸收峰之后波形相似特征明显,但由于后几个波段噪声大,故取1 952~2 451 nm,即波段[165,248]为特征区间,主要的吸收峰位置为2 258 nm和2 360 nm;碳酸盐化矿物的主要吸收峰为2 342 nm,所以取2 246~2 499 nm,即特征区间为[214,256]。由于仅需要对添加了峰值权重的特征区间内的曲线部分进行光谱角匹配,便可将同一类蚀变矿物的多条“异谱”端元(图1中仅为一部分)简化成1条端元参与分类计算。
分别取3种蚀变矿物的1条波谱端元作为参考光谱,根据其不同的特征区间和峰值位置,利用IDL编程语言实现上述蚀变信息提取算法。在阈值相同的情况下,权重系数分别取ω=1(普通SAM),ω=2,ω=3和ω=4进行试验。图2示出用传统的SAM算法(图2(a))和基于不同峰值权重的匹配算法(图2(b)-(d))提取蚀变矿物的结果。
图2 不同权重下蚀变矿物提取结果对比Fig.2 Comparison among results of altered mineral classification by setting different weights
从图2的分类图中可以直观地看出,用基于峰值权重的匹配算法提取的蚀变矿物比用传统的SAM方法提取的蚀变矿物更为全面而准确。在该岩心段图像与成矿有关的蚀变中,伊利石化占主导,成面状分布;伴有绿泥石化呈点状分布;而碳酸盐化较少。通过对蚀变矿物提取结果的分析和统计,可以较为精确地对蚀变信息进行定位,并促进定量化遥感在岩心矿化蚀变信息提取中的应用。
为了验证本文矿化蚀变信息提取方法的有效性,分别对传统SAM和在不同权重系数下基于峰值权重匹配方法提取的蚀变矿物分类结果进行总体精度和Kappa系数的计算和对比。总体精度反映了正确分类像元占总分类像元的比例;Kappa系数则从统计意义上表明被评价分类结果在多大程度上优于随机分类结果。利用ERDAS软件中的Accuracy Assessment模块进行了分类精度评估,该方法是将专题分类图像中的特定像元与已知类别的参考像元进行比较(即对分类数据与地面真值进行对比)。由于截取的区域较小,在图像中随机选取200个点(包括背景值),对每个随机点查验该点像元的光谱曲线,判断分类结果是否与地面真值相符,得到分类精度(表3)。
表3 不同权重系数匹配方法的分类精度Tab.3 C lassification accuracy of differentm atching methods by setting different weights
从表3可以看出,由于强化了特征区间内吸收峰的特点,基于峰值权重的方法比传统的SAM方法精度高。在ω=2时,分类总体精度可提高20%以上。但ω的取值并不是越大越好,而是随着ω的增大,提取精度降低。在实际工作中,ω要根据图像信噪比情况和同类物质吸收峰差异的大小来选择。
从岩心高光谱数据中提取矿化蚀变信息常用的SAM方法,由于自身基于整体光谱匹配的局限性,使其在信息提取时造成较大误差。因此,本文通过对局部光谱特征区间中的峰值合理地添加权重,对传统的SAM法进行改进,使其能够克服上述缺陷。实验结果证明:
1)用改进的SAM法提取蚀变矿物的总体精度从58%提高到70%~80%(取决于不同的权重系数)。
2)由于本文方法对端元的数量要求不高,每种蚀变矿物仅需1~2条端元光谱曲线即可参与光谱角计算,因而在提高分类精度的同时,也提高了端元选择和提取的效率。特别是在大数据量、大规模的岩心蚀变矿化信息提取和复杂繁琐的岩心编录工作中,本文方法应用效果显著。
3)通过对与铀成矿关系密切的3种矿化蚀变信息的提取和分析,可为进一步找出铀矿富集规律、缩小普查勘探的范围、提高找矿效率提供依据。
下一步的工作是开展如何根据岩心高光谱图像的信噪比及光谱吸收峰深度进一步确定峰值位置的权重系数,并从中寻找规律的研究。
[1] 张婷婷,唐菊兴,郭 娜,等.高光谱遥感在斑岩矿床蚀变信息提取中的应用[J].矿物学报,2011(s1):922.Zhang T T,Tang JX,Guo N,etal.Application of hyperspectral remote sensing in alteration information extraction of porphyry deposits[J].Acta Minalogica Sinica,2011(s1):922.
[2] 张杰林,黄艳菊,王俊虎,等.铀矿勘查钻孔岩心高光谱编录及三维矿物填图技术研究[J].铀矿地质,2013,29(4):249-255.Zhang JL,Huang Y J,Wang JH,et al.Hyperspectral drilling core logging and 3Dmineralmapping technology for uranium exploration[J].Uranium Geology,2013,29(4):249-255.
[3] 甘甫平,王润生.高光谱遥感技术在地质领域中的应用[J].国土资源遥感,2007,19(4):57-60.doi:10.6046/gtzyyg.2007.04.13.Gan F P,Wang R S.The application of the hyperspectral imaging technique to geological investigation[J].Remote Sensing for Land and Resources,2007,19(4):57-60.doi:10.6046/gtzyyg.2007.04.13.
[4] Kruse F A,Lefkoff A B,Boardman JW,et al.The spectral image processing system(SIPS)-interactive visualization and analysis of imaging spectrometer data[J].Remote Sensing of Environment,1993,44(2/3):145-163.
[5] Robila S.An investigation of spectralmetrics in hyperspectral image preprocessing for classification[C]//Geospatial Goes Global:from Your Neighborhood to theWhole Planet.ASPRSAnnual Conference,Baltimore,Maryland.2005:7-11.
[6] Hecker C,van der Meijde M,van derWerff H,et al.Assessing the influence of reference spectra on synthetic SAM classification results[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(12):4162-4172.
[7] 王 瀛,郭 雷,梁 楠.逐波段修正负相关的光谱角填图算法[J].火力与指挥控制,2013,38(2):22-25.Wang Y,Guo L,Liang N.A spectral anglemapper algorithm modified negative correlation band by band[J].Fire Control& Command Control,2013,38(2):22-25.
[8] 甘甫平,王润生.遥感岩矿信息提取基础与技术方法研究[M].北京:地质出版社,2004.Gan F P,Wang R S.Studies on Basis and Technique Methods of Remote Sensing Minerals Information Extraction[M].Beijing:Geological Publishing House,2004.
[9] 何中海,何彬彬.基于权重光谱角制图的高光谱矿物填图方法[J].光谱学与光谱分析,2011,31(8):2200-2204.He Z H,He B B.Weight spectral angle mapper(WSAM)method for hyperspectral mineralmapping[J].Spectroscopy and Spectral Analysis,2011,31(8):2200-2204.
[10] 唐 宏,杜培军,方 涛,等.光谱角制图模型的误差源分析与改进算法[J].光谱学与光谱分析,2005,25(8):1180-1183.Tang H,Du P J,Fang T,et al.The analysis of error sources for SAM and its improvement algorithms[J].Spectroscopy and Spectral Analysis,2005,25(8):1180-1183.
[11] 别小娟,张廷斌,孙传敏,等.藏东罗布莎蛇绿岩遥感岩矿信息提取方法研究[J].国土资源遥感,2013,25(3):72-78.doi:10.6046/gtzyyg.2013.03.13.Bie X J,Zhang T B,Sun CM,etal.Study ofmethods for extraction of remote sensing information of rocks and altered minerals from Luobusha ophiolite in east Tibet[J].Remote Sensing for Land and Resources,2013,25(3):72- 78.doi:10.6046/gtzyyg.2013.03.13.
[12] 张志军,甘甫平,李贤庆,等.基于ASTER数据的蚀变矿物信息提取——以哈密黄山铜镍矿区为例[J].国土资源遥感,2012,24(2):85-91.doi:10.6046/gtzyyg.2012.02.16.Zhang Z J,Gan F P,Li X Q,etal.The extraction of altered mineral information based on ASTER data:A case study of the Huangshan copper- nickel ore district in Hami[J].Remote Sensing for Land and Resources,2012,24(2):85-91.doi:10.6046/gtzyyg.2012.02.16.
[13] 徐元进.面向找矿的高光谱遥感岩矿信息提取方法研究[D].武汉:中国地质大学,2009.Xu Y J.Research of Prospecting- Oriented Approaches to Information Extraction of Rocks and Minerals Using Hyperspectral Remote Sensing Data[D].Wuhan:China University of Geosciences,2009.
[14] 邓书斌.ENVI遥感图像处理方法[M].北京:科学出版社,2010.Deng S B.Processing Methods of ENVIRemote Sensing Imaging[M].Beijing:Science Press,2010.
[15] 地质部情报研究所.矿物岩石的可见-中红外光谱及其应用遥感专辑(第一辑)[M].北京:地质出版社,1980.Information Research Institute of Geology Ministry.Remote Sensing Album of the Mineral Rock Visible-infrared Spectroscopy and Its Application(the First Series)[M].Beijing:Geological Publishing House,1980.
[16] 傅 锦,裴承凯,韩晓青,等.基于多元统计技术的铀矿蚀变信息高光谱模型[J].世界核地质科学,2007,24(3):166-171.Fu J,Pei C K,Han X Q,etal.Hyperspectralmodel for information of uranium mineral alteration technique based onmathematical statistics analysis[J].World Nuclear Geoscience,2007,24(3):166-171.
[17] 雷天赐,祝明强,周万蓬.高光谱数据挖掘在蚀变矿物识别与提取中的应用[J].资源环境与工程,2005,19(3):213-219.Lei TC,Zhu M Q,Zhou W P.Application on datamining of hyperspectrum to identification and extraction of alterationminerals[J].Resources Environment& Engineering,2005,19(3):213-219.