郑 纬,朱谷昌,,张远飞,石菲菲,王冬寅,李 红

(1.中南大学地学与环境工程学院,长沙410083;2.有色金属矿产地质调查中心,北京100012)

0 引言

矿化蚀变信息是地质演化过程中的地质异常事件,是矿床在生成、演化全过程中的成矿条件表现[1]。由于矿床与地质体的各种异常属性相关联,同时这些异常属性又与相对均一的地质体背景属性所有偏离和差别,矿化蚀变信息就是包含在地质体综合信息中的各种属性异常微弱信息,人们必须采用各种可能的方法去捕捉和分辨地质体各种属性的异常,才能更好地预测矿床。遥感地质方法正是通过识别遥感信息中的矿化蚀变信息进行成矿预测的一种技术手段。

多层次分离提取技术就是以岩石和矿物的电磁波特征反射谱带作为提取岩石蚀变带信息的理论基础,通过数学建模和最佳变量的选择,运用各种数学方法和图像增强手段,多层次地从遥感信息中逐步剔除背景信息(即干扰信息),一次次地分离提取出矿化蚀变信息(即目标特征信息),最后得到包括目标特征信息的图像[2]。

本文在“遥感蚀变信息多层次分离提取”的主技术框架下,通过遥感图像特征诊断,把复杂的遥感蚀变信息提取转化为遥感图像的背景、干扰与蚀变信息3个主要对象的分析;然后,应用光谱数据空间(即特征空间)的几何结构分析对背景、干扰与蚀变信息3个主要对象的聚类结构进行研究与空间定位,特别是对干扰因素进行识别与划分[3-5],在分析的基础上再进行菲律宾宿务岛地区矿化蚀变信息的检测与提取。

1 宿务岛矿区地质背景与数据库

工作区位于位于菲律宾宿务岛北部宿务省宿务市 Kansi村,距宿务市北约 80 km,面积约16.92 km2。矿区地形为低山区,炎热潮湿,植被为杂草和树木混合。

菲律宾群岛属于西太平洋新生代岛弧火山岩带,区内地层以新生界为主,广泛分布第三系中性火山岩,其基底为前侏罗系,分布于菲律宾群岛中部和西部局部地区。

菲律宾岛弧区火山活动早期以古新世及渐新世为主,主要分布在东部大断裂带的东缘;晚期以中新世及上新世为主,主要分布在大断裂的两侧。侵入岩以闪长岩为主,还有花岗闪长岩、石英二长岩和正 长岩等。可分为2期:早期为35～60 Ma(渐新世),与块状硫化物和斑岩铜(钼)矿有关;晚期为1.5～2.5 Ma,与富金的斑岩铜矿有关。火山热液型金银矿与晚期中(基)性火山岩有关。超基性岩主要分布于大断裂西侧,东侧较少。

本文选用了1景 TM数据,即113/53(时相为1992-06-29)。以 TM为代表的Landsat类多光谱遥感数据,其地面覆盖范围宽,空间分辨率、光谱分辨率能满足区域地质调查及地质信息提取的要求,数据形式易于增强处理,多波段优化组合的图像信息丰富,是中小比例尺区域遥感地质调查的理想数据源。

2 矿化蚀变图像特征分析与提取

2.1 矿化蚀变的多波段数据背景-异常子空间分析

根据多层次分离提取技术思路,首先对遥感图像的背景-异常子空间进行划分[6]。通过遥感多波段数据的背景端元数M值的估计方法,能够基本上确定出背景-异常子空间。背景包括岩石、土壤等地物,同时也包含大的干扰(如植被区、云区、水系等)。小的干扰是通过光谱数据的空间结构分析来识别。

背景端元数M值的估计是通过赤池信息准则(AIC准则)或最小描述长度准则(MDL准则)来实现的。

本文根据 Crosta法先选择 F(铁化)因子的TM1,TM3,TM4,TM5等4个波段组合计算,得到标准差与方差分数(即贡献率)(表1)、M值的估计函数 Q(M)值(表 2)。

分析表1的方差分数,前2个主成分的贡献率已达到95.711%,表明主要背景(含大干扰)地物为前2个主成分;第3和第4主成分的合计贡献率仅为4.29%,主要为一些小干扰与可能的铁化蚀变异常信息。主成分4的贡献率比较低,结合Q(M)值可以预计这里包含的除了噪声之外,还有小干扰或弱信息。

M是主成分的本征数,背景(含大干扰)的端元数目应该等于M+1。分析表2中Q(M)值的变化会发现,这些值具有2个台阶,第一个台阶值范围是0.845～0.853,表明主要存在3类背景地物(含大干扰);第二个台阶值范围是0.853～1.000,表明小干扰与可能的蚀变异常信息属于台阶变化。这同前面的方差贡献率分析是一致的。

表1 Crosta法TM1,TM3,TM4,TM5的标准差与贡献率Table 1 Standard division and contribution rate of TM1,TM3,TM4 and TM5 calculated with Crosta method

表2 M值的估计函数 Q(M)值Table 2 The estimated value of function Q(M)of M value

表3 Crosta法TM1,TM4,TM5,TM7的标准差与贡献率Table 3 Standard division and contribution rate of TM1,TM4,TM5 and TM7 calculated with Crosta method

选择 H(泥化)因子 TM1,TM4,TM5,TM7等4个波段组合计算,得到标准差与方差分数(即贡献率)(表3)、M值的估计函数 Q(M)值(表4)。

分析表3的方差分数,前2个主成分的贡献率已达到95.188%,表明主要背景(含大干扰)地物为前2个主成分;第3和第4主成分的合计贡献率仅为4.812%,这2个主成分含有一些小干扰与可能的泥化蚀变异常信息。主成分4的贡献率比较低,结合Q(M)值可以预计这里包含的除了噪声之外,还有小干扰或弱信息。

表4 M值的估计函数 Q(M)值Table 4 The estimated value of function Q(M)of M value

M是主成分的本征数,背景(含大干扰)的端元数目应该等于M+1。分析表4中Q(M)值的变化会发现,这些值的变化具有2个台阶,第一个台阶值范围是0.711～0.752,表明主要存在3类背景地物(含大干扰);第二个台阶值范围是0.752～1.000,表明小干扰与可能的蚀变异常信息属于这个台阶。这同前面的方差贡献率分析是一致的。

结合遥感图像分析,发现该地区植被覆盖面很广、干扰强烈,而且已构成主背景地物,水体面积也较大,而其他岩石裸露区则是小背景地物。很明显,一些小干扰与可能的蚀变异常信息会“湮没”在这2类背景地物之中。

2.2 基于背景异常分析的矿化蚀变异常的特征波段选择

在对遥感图像数据进行背景-异常子空间分析后,通过回归偏度曲线分析和地物光谱反射特征来选择信息提取的特征波段。

回归偏度曲线是用来描述线性回归的二维散点图的图形不对称性,回归偏度较大,则预示着有不同于背景的异常或干扰信息存在[7]。

图1是由剔除了水体影响后的7个波段的回归偏度曲线。图中显示出 TM4,TM5,TM3,TM6和TM1等5个波段的偏度较大,说明它们载荷了较多的背景、干扰或蚀变异常信息。由地物波谱的物理机制可知,偏度最大的 TM4反映的是植被信息;第2位的 TM5既有植被信息、岩石信息,还可能载荷有泥化蚀变信息;第3位的TM3包含了铁化蚀变信息;与 TM5形成反向偏度的 TM7应该是由泥化蚀变信息引起的;同时注意到 TM3与 TM1的偏度反差不大,TM3与 TM4的偏度反差主要是由植被造成的,即 TM4高、TM3低。

由上述的分析得知,TM4,TM5,TM6,TM3和TM1是该地区重要的特征波段。所以波段选择方案为:①剔除植被干扰的波段:TM4和 TM5;②岩石与裸露区识别的波段:TM5;③铁化蚀变异常信息识别与提取波段:(偏度曲线反映不明显,待光谱空间结构分析后再定);④泥化蚀变异常信息识别与提取:TM5和 TM7。

2.3 矿化蚀变的波谱空间几何结构分析

在选定特征波段后,对遥感数据在波段组合中的光谱空间几何结构进行分析,从而确定不同蚀变信息提取的波段选择。

遥感的多(高)光谱数据是一种多元数据集合,每1个像元代表1个波谱矢量。所以,图像多元数据集合在高维空间中形成1个点阵,这个空间点阵具有一定的几何结构。

多(高)光谱数据集合属高维空间中的点阵,高维空间的点阵几何体形态是无法直接观测的。但是,从低维空间入手,通过分析二维(2个波段)、三维(3个波段)数据的散点图,可以探讨高维情况下的数据点阵的结构及在空间中的聚类分布特征。

以图2为例来说明光谱数据的二维散点图所包含的丰富信息:

(1)二维散点图又称二维直方图,它是二维变量(波段)联合概率密度分布的几何表达,反映了二维变量数据空间的聚类结构。

图1 TM数据7个波段的回归偏度曲线Fig.1 Partial regression curve of TM data of the 7 bands

图2 光谱数据的二维散点图Fig.2 Planar scatter plot of spectral data

(2)二维散点图也是2个变量主成分分析的几何解释,其中 PC1为第1主成分轴,PC2为第2主成分轴。对于多波段(3个波段以上)主成分分析,比如Crosta定向主成分分析法,可以通过多个二维散点图的观测去重构脑海中的高维光谱数据点阵几何结构。

(3)若约定图形坐标系统的Bx(横坐标)表示具有吸收峰的波段,比如 TM1或 TM7,而By(纵坐标)表示具有反射峰的波段,如 TM3或 TM5。根据几何意义,坐标平面内的点P(bx,by),位于对角线(图2中绿色线)上方的数据点必定有by＞bx,即(by/bx)＞1.0;反之,对角线下方的点必定有by

图3 TM3,TM4二维散点图Fig.3 Planar scatter plot of band TM3 and TM4

(4)对于比值运算,例如,TM3/TM1,则等同于TM3为by,TM1为bx,且有 TM3/TM1=by/bx=tgθ,这里tgθ为坐标平面点P(bx,by)与原点P(0,0)联线的斜率。

基于二维散点图的上述重要性,所以,它将作为对图像光谱数据空间几何结构特征分析的重要技术手段。

本区遥感数据的主要波段基于二维散点图的空间几何结构分析如下。

图3的背景主轴是沿着 TM4波段方向的,所以,该工作区的主要背景地物是植被。其他地物(含蚀变信息)仅构成一个次椭圆。

图4表明,散点信息全在对角线下方,从物理意义而言,表明 TM3与 TM1的比值未能反映出铁化异常蚀变信息。

结合图3与图4,说明该地区的3价铁(Fe3+)的吸收谱带可能主要在0.90μm(TM4(0.76～0.90)),而不在 0.45μm(TM1(0.45～0.52))。TM3与 TM4能反映出铁化异常蚀变信息。

图5为 TM5与 TM7的二维散点图,由于植被与泥化蚀变信息均为 TM5值高、TM7值低,所以主要包含植被信息的背景中心在对角线上方,说明植被的干扰还是比较严重的。

图6是 TM5与 TM6的二维散点图,利用 TM6与TM5的比值来提取硅化信息。

图4 TM3,TM1二维散点图Fig.4 Planar scatter plot of band TM3,TM1

以上二维散点图所反映的主要地物在光谱空间的聚类结构,表明该地区提取蚀变异常信息主要受植被影响最大,其次是云区干扰。蚀变异常信息未能形成自身的“聚类”形态,而是与干扰一起构成“点群”,说明蚀变异常信息在空间上分布较分散,大片异常不多,是“小”而“碎”的异常较多。

经分析认为,本区的植被、水体及云区是3种主要的干扰地物。

2.4 矿化蚀变图像的识别与提取

通过对遥感图像光谱数据的综合分析,选择基于Crosta法的F因子与 H因子分析方法,根据不同蚀变信息提取的波段选择,以自行开发的软件对工作区遥感图像进行矿化蚀变信息提取。

(1)铁化信息提取。根据 Crosta法先选择 F(铁化)因子的 TM1,TM3,TM4,TM5等4个波段组合计算,得到特征向量矩阵(表5)。

经分析比较,第3主成分(PCA3)的载荷 TM3和 TM5为正值,TM1和 TM4为负值,即 TM3与TM1,TM4载荷符号相反。所以,PCA3应该属于铁化蚀变信息主成分。

表5 研究区TM 1,3,4,5波段的特征向量Table 5 Characteristic vector of band TM1,3,4,5 of the study area

在提取软件中对滤波处理后的PCA3主成分影像进行密度分割,通过分析最优分割段内离差平方总和随分割段数变化曲线发现,当分割段数达到14后,曲线趋于平衡,因此14为合理分割段数,分割段数14的段内离差平方总和为28 108 033.71,最优分割区间为:{255,238}{237,219}{218,200}{199,183}{182,162}{161,145}{144,127}{126,109}{108,93}{92,78}{77,62}{61,47}{46,1}{0,0}。

分析对比认为,第1、第2分割段(灰度级为219～255)为铁化蚀变信息异常区,对该灰度段作进一步的3段最优密度分割,并分别赋以红、黄、绿色,该灰度段以外的灰度级作为背景处理,赋为白色,得到铁化蚀变信息异常图(图7)。

(2)泥化信息提取。选择 H(泥化)因子的TM1,TM4,TM5,TM7等4个波段组合计算,得到特征向量矩阵(表6)。

由分析比较,第4主成分(PCA4)的载荷 TM5为负值,TM7为正值,即 TM5与 TM7的载荷符号相反。所以,PCA4应该属于泥化蚀变信息主成分。根据PCA4取反向的数值,TM5与 TM7提取的泥化蚀变信息图像见图8。具体软件操作与铁化提取类似,采用最优彩色分割对泥化信息进行提取。

表6 研究区TM 1,4,5,7波段的特征向量Table 6 Characteristic vector of band TM1,4,5,6 of the study area

图7 TM1,3,4,5的 PCA3主成分提取的铁化信息图像Fig.7 The extracted iron-based information Images of principal component PCA3 of TM1,3,4,5

(3)硅化信息提取。根据前文的分析,这里利用ETM6/ETM5来提取硅化信息。

2.5 矿化蚀变信息提取结果分析

从图7、图8、图9中可以看出主要蚀变区集中于图幅的南部,但在图8与图9中图幅的东北部泥化、硅化均表现较弱,但却有一致性。

(1)结合构造解译可见,图幅南部的铁化、泥化、硅化均有所反映,它正好与宿务岛火山弧中的一个环形构造带相吻合,在环形构造的中心即是托莱多斑岩铜矿(大型)的产出部位。托莱多斑岩铜矿具有良好的蚀变分带:中央为黑云母的钾化带,其外侧为硅化较强的石英绢云母化带,再外则是泥化带,最外侧是青磐岩化带。钾化带与青磐岩化带的铁质都比较高。因为钾化带的黑云母、青磐岩化中的绿泥石和黄铁矿均富含铁质,所以会有较强的铁化蚀变显示。在石英绢云母化带内,SiO2是主要的成分,必然有硅化的显示,而泥化带的主要成分是高岭石、伊例石、叶腊石等,故而出现泥化反映。

(2)在图幅的东北角也有一个较小的环形构造和一组放射状断裂,其南侧还有几个次火山岩岩筒,这里的铁化较弱,但泥化、硅化却较明显,虽然强度 不高,但铁化、泥化、硅化也有吻合性,可能也是具有铜矿化的一种蚀变标志,可以列为找矿标志。

(3)东南角大片的泥化与硅化可能为滨海沙滩,那里可能有大量的石英和泥质的海滩,他们并不是矿化的反映。

3 结语

将蚀变信息提取结果与地质矿产资料相比较,发现区内己知的矿床(点)和矿化蚀变带均在所提取的蚀变信息异常内;通过野外实地调查表明,采用的矿化蚀变图像识别与提取技术是有效和可靠的。

[1] 赵鹏大,陈永清,刘吉平,等.地质异常成矿预测理论与实践[M].武汉:中国地质大学出版社,1999.

[2] 张远飞,朱谷昌,吴德文.地质矿产调查的遥感蚀变信息多层次分离提取技术与应用[C]∥中国遥感应用协会2007年学术年会论文集,2007.

[3] 张远飞,吴德文,朱谷昌,等.遥感蚀变信息检测中背景与干扰问题研究[J].国土资源遥感,2008(2):22-26.

[4] 张远飞,杨自安,朱谷昌,等.遥感图像蚀变信息检测中的光谱数据空间结构分析[J].遥感信息,2009(1):3-9.

[5] 张远飞,吴健生.基于遥感图像提取矿化蚀变信息[J].有色金属矿产与勘查,1999,8(6):604-606.

[6] 张远飞,杨自安,张普斌,等.高(多)光谱数据的背景-异常子空间模型研究[J].地球信息科学学报,2009,11(3):283-285.

[7] 张远飞,吴德文,张艮中,等.高光谱数据的波段序结构分析与应用研究[J].国土资源遥感,2010(1):30-38.

[8] 张玉君,曾朝铭,陈薇.ETM+(TM)蚀变遥感异常提取方法研究与应用——方法选择和技术流程[J].国土资源遥感,2003(2):44-49.

[9] 王建平.基于遥感的河南卢氏西部地区蚀变信息提取与分析[J].地球信息科学学报,2007(6):111-115.

[10] 李智勇,郁文贤,匡纲要,等.基于高维几何特性的高光谱异常检测算法研究[J].遥感技术与应用,2003(6):379-383.

[11] 李智勇,匡纲要,郁文贤,等.基于高光谱图像主成分分量的小目标检测算法研究[J].红外与毫米波学报,2004(4):286-290.

[12] 方洪宾,李志中.遥感化探信息综合分析在地质找矿中的应用研究[J].国土资源遥感,1998(4):33-36.

[13] 陈彩芬,舒宁.SAR影像与 TM影像的几种融合处理方法[J].国土资源遥感,1999(4):53-57.

[14] 吴德文,朱谷昌,吴健生,等.青海芒崖地区岩石光谱特征分析及应用[J].国土资源遥感,2001(4):28-34.

[15] 赵元洪,张福祥.波段比值的主成分复合在热液蚀变信息提取中的应用[J].国土资源遥感,1991(3):12-26.