张玉君曾朝铭陈薇
(中国国土资源航空物探遥感中心,北京)
摘要:从原理、模型试验及实际数据处理等方面分析对比了比值法、光谱角法及主分量分析法的优劣;选择了以主分量分析为主、光谱角法为辅进行蚀变遥感异常的提取;引用了误差理论某些基本概念,以标准离差σ作为遥感异常切割的尺度;建立了“去干扰异常主分量门限化技术流程”。以西藏驱龙—甲马蚀变遥感异常群为例,展示了此技术的效果,并与高光谱研究结果进行了比较。
关键词:主分量分析;光谱角法;比值法;异常切割量化尺度;门限化
中图分类号:TP753 文献标识码:A文章编号:1001-070X(2003)02-0044-06
本文系中国地质调查局“西部主要成矿带遥感找矿异常提取与应用研究”项目(200215000008)资助
引言
提取蚀变遥感异常信息的方法有多种,常用或较常用的有比值法、主分量分析法、光谱角法以及它们的混合法[1~7]。近期新出现但较少被使用的有小波分析、神经网络、分数维几何(分形或分几)等方法。为了选择最佳方法,我们以东天山139/31景、东昆仑135/35景和冈底斯137/39景ETM+或TM数据为例,对常用的3种方法异常信息提取效果进行了对比试验,其结果表明,3种方法所获互洽性异常虽有,但占比例不大。因为除139/31景外,其余两景中已知矿床、矿点分布较少,缺少评判3种方法优劣的足够已知条件,故而萌生进行模型试验的想法,即通过简单的模型试验来直观地展示3种方法(比值法、主分量分析法及光谱角法)的特点,对比其优劣,并以此为基础确定技术流程的主体。
1模型试验
模型一
模型设计
取一组二维数据点,它们在平面上呈椭圆形分布,而各点的数值按级差递增,且与坐标相关。第一维数值从左向右递增(彩版附图14(6)上左),第二维数值从下向上递增(彩版附图14(6)上中)。此设计是为了保证试验结果的直观性。这一模型的原始数据(彩版附图14(6)上右)可比拟为一个2波段的TM图像,如TM5和TM7。待提取异常应具备的条件是TM5为高亮度值,TM7为低亮度值。
方法试验结果
对模型—分别进行主分量分析、比值法处理及光谱角填图等方法试验,对所获第二主分量进行高端切割,3种方法所获异常如彩版附图14(6)左下图所示,其中,光谱角法异常提取以第二主分量异常区高位中心点的数值作为参考谱;比值法则为第一维/第二维,并对其进行异常切割。
从彩版附图14(6)可以直观地看出:①主分量分析所获异常位于第一维的高值区;②光谱角法所获异常区各点有相似的特征向量角,异常区内既包括高值点(与主分量异常互洽)也包括低值点;③比值法比较的是比值的大小,故该方法所获异常与原始值的高低无直接关系,和主分量及光谱角法结果的互洽区不大。这就不难理解3种方法所获异常不可能全部互洽的原因。如果提取遥感异常的原则是以大型、特大型矿(应有高异常)为主要目标,那么,提取遥感异常的首选方法应是主分量分析法。
如果在比值计算之前进行数据截取,可以减少低值区异常(彩版附图14(6)的下右),但若先比值后截取,则效果较差(彩版附图14(6)的下中)。当然,对光谱角法也可以通过数据掩模处理来减少低值区异常,且无论掩模在光谱角处理之前或之后做,其结果相同(彩版附图14(6)的下中和下右)。
模型二
模型设计
取模型二图像(540×360×4),分割为180×180×4六个子图,每一子图对角线皆为255个像素。对每一子图的4个图层分别送入某种岩性的TM1、TM4、TM5、TM7波段亮度值,且令其按对角线方向从0渐变至相应波段的TM亮度值。为达此目的,首先取一辅助图像(180×180×1),用试验图像工具(Test Image)在对角线方向形成从0~255的渐变灰阶,并借助此辅助图像将各岩性的各波段亮度值在对角线方向做线性拉伸,然后逐一镶嵌入模型二各子图各波段的相应位置,形成岩性的渐变图像组合。6个子图分别对应于金矿蚀变岩、植被、白色大理岩、辉长岩、冰和盐碱地(数据取自文献[6])。这样,便形成一个有6种已知岩性的TM1、TM4、TM5、TM7四个波段的数据模型。
方法试验结果
对模型二的4个波段数据进行主分量分析,从本征向量可知,第2主分量表征冰,第3主分量表征植被,第4主分量表征蚀变岩。对它们分别进行异常切割,以绿、蓝、红三色示于彩版附图14(7)中。可以看到,白色大理岩形成对蚀变岩的干扰异常,经用光谱角法提取蚀变异常(彩版附图14(7)中的黄色区),可以消除此干扰。比值法所得结果未在彩片中示出,它在除辉长岩外的5种岩性中均形成异常,这是它们在TM5波段上的亮度值均高于TM7的必然结果。
2三种方法的特点讨论
比值法(RM)
(1)原始值的高低不直接左右比值的大小,也就是说,高低值均可有同样的比值,低值区异常不可避免(图1),这就使得比值结果中包含了大量假异常(如阴影区、水域等)。
(2)差和比值法是一种非线性处理方法,该方法高值区压缩,对低值区拉伸,这对于提取高值区异常不利; 的实际值均在0~1/3之间,动态范围较小,也是不利因素(图2)。
图1波段比值原理
图2差和比值原理
(据丰茂森,1992)
(3)数据先截取再比值较先比值再截取更有益于减少低值区异常(图3,彩版附图14(6)下中和下右)。
虽然经数据截取所做比值对低值区异常有所抑制,但并不能从根本上改变比值法的前两个弱点,故在选择遥感异常提取方法时,不得不将比值法结果仅作为参考。
图3数据截取比值原理
图4光谱角法原理
光谱角填图法(SAM)
光谱角法把每一个多维空间点以其空间特征向量来表征,并以空间向量角的相似性作为判据。它是一种监督分类,要求每一类别有一个已知参考谱。此参考谱可以是地面实测入库光谱,也可以是已知条件的图面单元的统计入库结果(又称图像采样)。
为了直观,设三维空间点P在彩色坐标系中的特征向量为OP,以此向量为轴作小角锥(图4),凡位于此小角锥内的空间点都视为相似的。
根据线性代数理论[9],向量α,β间夹角θ为
张玉君地质勘查新方法研究论文集
式中,(α,β)为n维向量α,β的内积,| a |、|β|分别为向量a,β的长度。
当存在已知矿点或矿床时,可以利用光谱角法圈定与其有相似谱特征的靶区,以减少主分量分析所获异常中的非矿异常(如在139/31景中所做);当存在两种以上已知矿点时,可以用光谱角法对主分量分析异常进行类别区分(如在137/39景中所做)。这两点是光谱角法在异常信息提取中对主分量分析法可以起的重要辅助作用。
主分量分析(PCA)
原理
遥感图像各波段间存在一定的相关性,为了减少相关性对分类的影响,常使用主分量分析法去相关。主分量分析基于变量之间的相互关系,在信息总量守恒的前提下,利用线性变换的方法来实现去相关性。由于所获各主分量之间不相关,故各主分量之间信息没有重复或冗余。
主分量分析的这一基本性质在蚀变异常信息提取中被充分利用。TM多波段数据通过PCA所获每一主分量常常代表一定的地质意义,且互不重复,即各主分量的地质意义有其独特性。用TM1、TM4、TM5、TM7等4个波段进行PCA,对代表羟基矿物主分量的判断准则是:构成该主分量的本征向量,其TM5系数应与TM7及TM4的系数符号相反,TM1一般与TM5系数符号相同。依文献[7]中有关地物的波谱特征,羟基信息包含于符合这一判断准则的主分量内,故此主分量可称为羟基异常主分量。
用TM1、TM3、TM4、TM5等4个波段进行PCA,对代表铁染物主分量的判断准则是:构成该主分量的本征向量,其TM3系数应与TM1及TM4的系数符号相反,TM3一般与TM5系数符号相同。同理,可将该主分量称之为铁染异常主分量。
由于异物同谱现象的存在,这些异常主分量还包含着非蚀变的地质因素,有待用其他计算方法乃至借助地质知识经目视解译加以区分。
借助模型试验探讨PCA的灵敏度极限
当用PCA提取遥感异常信息时,无论是羟基异常(用TM1、TM4、TM5、TM7)还是铁染异常(用TM1、TM3、TM4、TM5),都常常出现在第4主分量[3~7],其本征值都是最小的,也就是说,蚀变遥感信息占全景信息量的很小份额,那么,蚀变信息量少到何种程度将被PCA丢失,即PCA的检测灵敏度下限是多少?
为此,首先我们设计图像模型(三)(1020×255×4),并将其分割为255×255×4四个子图,依次由左向右送入金矿蚀变岩、植被、辉长岩及冰的TM1、TM4、TM5、TM7波段数据,并令其按行或列方向递变。经PCA提取切割,获蚀变异常(红)、植被异常(蓝)和冰异常(绿),见彩版附图14(8)。
然后,对蚀变岩区逐步减半,每次都以辉长岩数据取代,再做PCA提取,直至蚀变岩区仅剩一行(255×1×4),PCA仍有相应的结果,但此时相对信息量已小于,程序给出报告是(表1),可视为已达检出限,这时,蚀变岩原始数据在全图的份额(1/4×1/255=1/1020)小于千分之一。因此,可初步认为,用PCA提取蚀变信息的检出限优于千分之一。
表1模型三试验结果统计
将表1第二行本征向量分别除以TM5的系数,得到一组新数,把它们做为投影系数(表2),将TM1、TM4、TM7投影在TM5上,得到图5。这有助于理解PC4的构成。
表2PC4系数换算
图5PC4与TM1、TM4、TM5及TM7的关系
3误差理论某些基本概念的引用
在图像处理中,经常使用概率密度分布曲线(简称直方图),于是便产生两个问题:
(1)如何理解TM数据直方图常接近正态分布?
(2)是否可以使用标准离差σ作为遥感异常切割的尺度?
正态分布
我们知道,观测值的偶然误差成正态分布,而且早在1795年德国数学家高斯就推导出偶然误差或然率曲线的函数表达式,即高斯误差分布定律[8]
张玉君地质勘查新方法研究论文集
式中,σ称为标准误差。如果取k倍的标准误差,那么,任一观测值的误差介于±ka之间的或然率p如表3所示。
表3或然率与误差的关系
结合具体情况,当所面对的地区足够大时,可以近似地认为TM各波段数据以及其线性处理结果的统计性规律具有类似(在不严格的概念下)正态分布的特点:
(l)常常只有一个中心,即众值(当出现双峰时,往往是存在两类主宰性地质因素);
(2)小偏离比大偏离出现的机会多;
(3)大小相等、符号相反的正负偏离机率接近,直方图近似对称于y轴;
(4)极大的正偏离和极小的负偏离机率都很小,直方图向两端迅速衰减。
究其原因,当地区较大时,不同岩性的出现具有随机特性。故而,描述偶然误差的高斯误差分布方程及相关概念可以被借用于TM数据处理中。
σ的借用
TM数据处理以多元分析为基础,在多元分析中,对应于误差理论中称之为标准误差的σ,是标准离差(或标准偏差),其定义为:
张玉君地质勘查新方法研究论文集
既然TM数据及其线性处理结果一般均有近似正态分布的直方图(如图6所示,为137/39景羟基主分量直方图),那么,在做异常切割或数据切割时,便可借用σ这个表征正态分布曲线的尺度。例如,主分量分析结果可以把均值 理解为代表区域背景,利用 确定异常下限和划分异常强度等级。异常总面积可用(1-p)/2近似计算,3σ级异常总面积仅约在×10-3之内。有了这一尺度,切割异常时可以减少主观任意性,并使操作较为规范化。
4遥感异常提取流程
图6137/39景羟基主分量直方图
在以上方法选择及多景图像蚀变异常提取研究的基础上,形成了独具特色的方法技术,即“去干扰异常主分量门限化技术流程”(图7),并以此作为蚀变异常提取快速遥感扫面的基本方法。
图7去干扰异常主分量门限化技术流程
流程图包括预处理和信息提取两大部分,预处理的作用十分关键,它的许多重要问题(如多种干扰去除、大气影响粗校正、太阳照射条件归一化等)尚须另文阐明,信息提取和后处理的技术细节也待补充。
5实例
现以甲马—驱龙火山沉积盆地铜多金属矿田为例,展示应用该流程进行异常提取的效果。
研究区位于冈底斯火山岩浆弧铜金多金属成矿带的冲江—甲马亚带中,区内广泛发育燕山—喜山早期岛弧型火山岩、碎屑岩和碳酸盐岩,重熔、同熔型钙碱性岩浆杂岩,以及新生代碰撞后期生成的花岗斑岩等,具有火山弧→弧内盆地-岩浆弧等构造演化序次和时空叠加特点和与之相关但性质各异的热液流体蚀变和成矿作用,生成了相互关联且各具特色的矿化,积聚了巨大的矿产资源潜力。本区中-晚侏罗世火山弧阶段,以钙碱性火山岩系列组合为主,具多个火山-沉积韵律旋回,火山热液蚀变和矿化较广泛,如驱龙南火山喷流-热液型铜多金属矿,矿体呈似层状产于火山岩中,长70~1000m,厚10~50m,脉石矿物以石榴石、透辉石、符山石等为主,矿石矿物有黄铁矿、黄铜矿、斑铜矿、方铅矿和闪锌矿等;晚侏罗—白垩世弧内盆地阶段,堆积了巨厚的碎屑岩-碳酸盐岩层,成矿作用有海底热水喷流特点,如甲马铜多金属矿床,矿体呈似层状产于碳酸盐岩与上覆砂板岩之间的透辉石、斜长石矽卡岩中,矿体长,平均厚,潜在铜资源量达200万t,成矿作用与海底热水喷流、交代作用有关;白垩—新近纪岩浆弧阶段,大量中酸性岩基、岩株侵位于盆地中,与岩浆活动相关的热液流体蚀变、矿化作用广泛,如驱龙斑岩铜矿,矿化主要见于花岗闪长斑岩中,岩体呈小岩株状,出露面积约,岩体蚀变强烈,分带明显,从中心向外依次为石英-绢云母化带、高岭土-绢云母化带及青盘岩化带等,主要矿物有黄铜矿、斑铜矿和辉钼矿等。
遥感异常信息提取用137/39景数据(彩版附图14(9)),同时处理了两个时相的数据。效果评估准则是:①从多时相数据中提取异常的稳定性;②异常表征已知矿床、矿田的能力;③异常空间分布的规律性。3种方法提取的异常对比结果表明,主分量分析法提取的异常无论在异常的稳定性、分布规律性及与已知矿床、矿田对比性方面,均极大地优于比值法提取结果;光谱角填图法提取的异常,其“三性”稳定程度与参考光谱的选取有密切关系,科学地选取有代表性(波谱涵盖范围)的参考光谱是取得与其他方法之间可比性的关键。此外,甘甫平等[10]在驱龙地区,利用EO-1高光谱成像遥感数据开展的蚀变矿物提取,获得了高铝绢云母、低铝绢云母、高岭石、绿泥石+孔雀石化等4种矿物异常,这为我们提供了对比不同类型数据、不同方法提取效果的机会。二者相比,除大片出现在阴影区的绿泥石+孔雀石化在我们的图中没有对应的异常显示外,前3种蚀变矿物异常区均与我们利用主分量分析和光谱角法从ETM+数据中提取的异常区基本相似,其差异主要表现为,前者提取的是单矿物分布,后者提取的是蚀变矿物的集合体异常。
本景遥感异常提取的效果仅从驱龙和甲马两个异常群的找矿前景就可见一斑。驱龙异常群呈东西走向的带状分布,长约15km宽约5km,全区铁染异常极为发育,羟基异常多呈斑块状叠置于其上。在已圈定的6个异常中,多数具有从中心向四周、从羟基到铁染、场强从强到弱逐步过渡的特点。其中,4、5号异常是目前正在勘探的斑岩铜矿区:1号异常呈北西走向,中部已证实有岩石蚀变,推测两端可能有斑岩体隐伏;2号异常西部及3号、6号异常,深部均可能有斑岩体。全区经光谱角法检测可知,不仅甲马型异常特征显著,而且1、3、5、6号异常均有驱龙斑岩铜矿区(4号异常)蚀变波谱特征。根据本区异常带规模、结构特征和铁染异常强大等推测,该异常群具有寻找大型斑岩型和热液型矿的前景。甲马异常群有5个异常组成,其中1号异常位于甲马铜多金属矿体上方,异常辅呈南北向,长,宽,铁染、羟基异常紧密相伴,前者异常强大,连续完整,后者的高强度异常呈斑块状叠置于铁染之上,形成4个强羟基蚀变中心。根据已知矿体及含矿矽卡岩带上遥感异常显示仅为不多的散点,异常的轴向与已知矿体倾向之间呈450夹角,以及异常规模巨大,具有斑岩型和热液型双重结构特征等推断,该异常形成时间应晚于甲马铜多金属矿,异常区深部可能有呈南北向带状侵位的斑岩体和历经数次的热液活动,甲马已知铜多金属矿深部矿体可能已被岩体侵入和热液活动强烈改造,据此推测该区可能是一个寻找大型富矿的远景区。
6结论和建议
(1)模型一、二试验及137/39景应用实例证明,以寻找大型、特大型矿为目的的遥感异常“遥感扫面”工作方法,应以主分量分析法为主,以光谱角法为辅;
(2)经模型三试验可初步认为,PCA的遥感异常检出限(蚀变岩所占像素比例)优于千分之一,但这一问题尚待进一步研究;
(3)ETM+与高光谱分析结果的高度吻合,增强了利用ETM+数据进行遥感异常快速扫面的信心;
(4)为了便于交流,建议对有关方法和异常的代号作如下规定:
遥感异常RSA(Remote Sensing Anomaly)
羟基异常OHA(Hydroxylate Anomaly)
铁染异常FCA(Ferric Contamination Anomaly)
主分量分析PCA(Principal Component Analysis)
光谱角制图SAM(Spectral Angle Mapper)
比值法RM(Ratio Method)
背景BG(Background)
灰度GS(Grey Scale)
参考文献
[1]丰茂森.遥感图像数字处理[M].北京:地质出版社,1992
[2]Kendall Analysis[M].England:Charles Griffin and Company Limited,1975
[3]Crosta A P, McmMoore of landsat thematic mapper imagery for residual soil mapping in SW Minais Gerrain[A].In:Proceedings of the 7 th(ERIM)Thematic conference: Remote Sensing for Exploration Geology[C].Calgary,1989,1173-1187
[4]Loughlin W component analysis for alteration mapping[A].In:Proceedings of the 8th Thematic conference on Geologic Remote Sensing[C].Denver,USA,1991,293-306
[5]李昌国,张玉君.试用主分量分析方法提取澜沧江兰坪地区铜矿化蚀变遥感信息[J].国土资源遥感,1997,(1):20~30
[6]张玉君,杨建民.基岩裸露区蚀变遥感信息的提取方法[J].国土资源遥感,1998,(2):46~53.
[7]张玉君,杨建民,陈薇.ETM+(TM)蚀变遥感异常提取方法研究与应用——地质依据和波谱前提[J].国土资源遥感,2002,(4):30~36.
[8]冯师颜.误差理论与实验数据处理[M].北京:科学出版社,1964.
[9]居余马,胡金德等.线性代数及其应用[M].北京:中国电视大学出版社,1986.
[10]甘甫平,王润生,杨苏明.西藏Hyperion数据蚀变矿物识别初步研究[J].国土资源遥感,2002,(4):44~50.
THE METHODS FOR EXTRACTION OF ALTERATION ANOMALIES FROM THE ETM+(TM)DATA AND THEIR APPLICATION:METHOD SELECTION AND TECHNOLOGICAL FLOW CHART
Zhang Yu jun,Zeng Zhao ming,Chen Wei
(China Aero Geophysical Survey andRemote Sensing Center for Landand Resources,Beijing100083,China)
Abstract: This is the second paper among the serial articles on the new parameter for the prediction of the mineral resources: alteration RS RM(Ratio Method), SAM(Spectral Angle Mapper)and PCA(Principle Component Analysis)were compared with each other by the analysis of their principles and by the model was selected as the main method and SAM as the supplementary essential concepts of the error theory were standard deviation σ was used as the measurefor anomaly flow chart“De-interfered Anomalous Principal Component Thresholding Technique”was practical application example for Qulong-Jiama anomaly group in Tibet was given in comparison with the hyperspectral results.
Key words: PCA;SAM;RM;Quantitative measure of anomaly slicing;Thresolding
原载《国土资源遥感》,2003,。