分形滤波技术在新疆黄山-镜儿泉镍铜成矿带中的应用

娄德波,肖克炎,左仁广 ,丁建华

1)中国地质科学院矿产资源研究所,北京 100037;2)中国地质大学(武汉)资源学院,湖北武汉 430074;3)地质过程与矿产资源国家重点实验室,中国地质大学(武汉),湖北武汉 430074

“横看成岭侧成峰,远近高低各不同,不识庐山真面目,只缘身在此山中”,这是北宋大词人苏轼所做诗篇“题西林壁”中的诗句,它描写庐山变化多姿的面貌,并借景说理,指出观察问题应客观全面,如果主观片面,就得不出正确的结论。地球科学亦是如此,如果想正确地理解地球物理、地球化学数据,也必须在客观地质事实的前提下,全面把握地质过程发生发展的客观规律,选择反映或刻划其本质的技术方法。由于多数地球物理和地球化学测量数据往往是反映地表及地下多种地质体和地质过程的最终叠加结果(Xu et al.,2001),作为研究地质体、地质过程以及矿产勘查的基本手段,地球物理、地球化学场分解和异常的提取已经成为众多学者研究的重要内容(Cheng et al.,2000)。基本的方法有两类,即频率分析和空间分析。频率分析主要是研究值的频率分布,如直方图、Q-Q图、盒式图以及概率图等,主要目的是将样品从不同的总体或者源区分离出来,其基本假设是不同的地质过程产生不同的总体,且这些总体是可以区分的。空间分析主要是从空间方面处理地球化学数据,像使用多种元素进行组合异常分析(马生明等,2011)以及用于产生各种地球物理、地球化学空间分布图的各类插值方法,如克里格、移动平均、径向基函数和张力样条函数等。虽然这些方法可以用来处理局部和区域地球化学数据,但是也有很多不足之处。如通过频率分析方法得到的一个或者两个总体阈值在一些地球化学背景变化强烈的地区不能起到很好的效果;使用频率分析在多个总体严重重叠的区域区分不同总体是困难的;而且这些方法不能表征地球物理、地球化学异常的空间分布和几何特征,而这些对于矿产勘查是尤为重要的。对于空间分析,诸如克里格或者移动平均一般不能分离异常;而且一些方法,诸如移动平均等又通常存在人为确定操作窗口大小和形状等问题,这为地球化学图的编制带来了一些主观因素(Xu et al.,2001)。人们已经做了许多工作来提高这项技术,诸如因子克里格(Goovaerts,1992)、C-A方法(Cheng et al.,1994;徐明钻等,2010)、U统计(Cheng,1999)以及S-A(Cheng et al.,2000;成秋明,2001;张焱等,2011)等。其中新近发展的使用频率和空间信息的分形滤波技术(S-A)在一定程度上克服了上述缺点,具有人为介入较少,适应数据本身由于地质作用和成矿作用所形成的客观存在的各向异性及内部结构复杂性等特点,因此被广泛用于场的分解和信息提取。本文将以新疆东天山黄山-镜儿泉镍铜成矿带Ni元素为例来介绍分形滤波技术。之前,首先简单介绍分形滤波技术的一般原理。

1 分形滤波技术

在矿产预测工作中,同一类型的矿床当处于不同的地质环境、不同的埋藏深度时,其所处地区的地球化学测量数据将表现有很大的不同,这是地质环境和矿床本身综合作用的结果,因此将矿致异常从复杂的背景场中分离出来,显得尤为重要。同一类型矿床的地球化学特征通常具有自相似性,因此基于多重分形理论的分形滤波技术(S-A)便为此问题的解决提供了技术支持。分形滤波技术的基本地质假设是特定地质过程产生的地球化学图具有分形特征,地球化学值的分布与尺度之间遵循一定的幂率关系。以前的分形及多重分形工作表明大部分地质过程产生的图像具有尺度不变性,它通常表现为自相似性或自仿射性,这种性质能够在空间域和频率域中来衡量(Turcotte,1997)。在空间域中,尺度不变性通常与图像的空间几何形状、值得分布有关,几何形态的变化与值的变化相一致。在频率域,尺度不变性主要通过功率谱分布来表现(Lewis et al.,1999)。分形滤波器通过基于频率域中功率谱与“面积”的幂率关系来定义,目的是把功率谱分解成几个有着相同或相似的尺度性质的组分,然后应用该滤波器从原始图中分解背景和异常以及提取其它有用的信息。由于分形滤波技术是在功率谱分析技术和分形浓度-面积技术(C-A)的基础上发展起来的,因此在介绍分形滤波技术之前,首先简单介绍功率谱分析和C-A技术。

1.1 功率谱分析

从功率谱分析的观点来看,一幅图像(如地球化学图)能看作存在于两个域内:空间域和频率域。空间域上的图像或信号可以看作为不同波长的叠加信号,这种图像或信号能够被分解成几个端元组分,每个组分对应着相应的频率范围,且在空间域中复杂的卷积关系在频率域中会变为简单的乘积关系,如公式(1)和(2)。由于傅里叶变换可以把空间域中的图像或信号转化成频率域内的图像或信号的功能,因此可以在空间域内,把二维图像通过傅里叶变换(Dobrin et al.,1988),公式(3):转化进入频率域,以及通过公式(4)生成功率谱。

f(x,y),f1(x,y),f2(x,y)为空间域信号或图像,F(Kx,Ky),F1(Kx,Ky),F2(Kx,Ky)为傅里叶变换后频率域中的函数,⊗代表卷积计算,Kx和Ky是波数,与x和y相对应。

其中,S(Kx,Ky)为功率谱。

函数R(Kx,Ky)和I(Kx,Ky)可以通过乘以一个滤波器函数G(Kx,Ky)来修正,公式如下:

同时,通过反傅里叶变换(Dobrin et al.,1988),如公式(6)可以把分解的经过修正的各种频率重新组合并在空间域中形成新的图像或信号。

1.2 C-A 技术

成秋明等发展了在空间域内分离背景和异常的浓度-面积分形方法(C-A)(Cheng et al.,1994)。它已被用于分析各种各样的地球化学数据,如湖泊沉积物(Cheng et al.,1996),水系沉积物(Cheng et al.,1994)以及岩石样品等等。它主要是在一个浓度和面积取对数的坐标系中,表明在面积和元素浓度值之间存在着一定的幂率关系,即:

C为常数,α为奇异性指数,对于不同的地球化学浓度值区间,有不同的α值,它们可以通过最小二乘法(LS)使用不同的直线段拟合。直线的交叉点及所对应的浓度值作为分界值把浓度值分成几个组分,达到背景和异常分离的目的。

1.3 分形滤波器(S-A)

在Kx−Ky二维平面内,高功率谱值主要分布在平面中心( Kx≈0,Ky≈0)对应着低频范围,随着离中心距离增大,功率谱值随之降低,但不呈严格对应关系,这是因为功率谱的空间分布不仅取决于波数,还取决于功率谱函数本身。在频率域中功率谱值的空间分布可以使用 C-A方法来描述,即将空间自相似性引入对场的滤波和分解数据处理中,不仅考虑了频率的相对大小,而且还顾及了频率的空间分布和功率谱的空间自相似性。人们已经开始从多重分形角度研究尺度不变性,并且已经研究确定,在S-A二维平面内,面积A(≥ S)与功率谱(S)之间存在着幂函数关系。

其中,不同的β值可以在在功率谱和面积取对数的坐标系中,通过最小二乘法(LS)拟合的方式来计算。这样不同的β值根据自相似性将功率谱值分成几类,分别对应着高功率谱值和低功率谱值,这些被不同范围覆盖的面积一般和高频与低频范围相关,但不完全一致,这就是 S-A滤波方法。这样各种滤波器可以被构建,如若获取背景信息,则可使高值S(Kx,Ky)=0,低值S(Kx,Ky)=1;反之,使高值S(Kx,Ky)=1,使低值S(Kx,Ky)=0。

1.4 分形滤波技术的构建过程:

(1)通过使用各种插值方法,如移动平均,克里格等形成地球化学图像;

(2)通过傅里叶变换生成功率谱实部R(Kx,Ky),功率谱虚部I(Kx,Ky)和功率谱S(Kx,Ky)本身的图像;

(3)根据自相似性原理,使用C-A技术把功率谱图像的功率谱值分解成不同的类,分别对应着高、中、低等不同范围的功率谱值;

(4)根据不同的需求,构建不同的滤波器G(Kx,Ky),如为了提取异常可令低值S(Kx,Ky)=1,高值S(Kx,Ky)=0;

(5)通过乘以一个滤波器函数G(Kx,Ky)来修正实部R(Kx,Ky)和虚部I(Kx,Ky);

(6)修正后的实部和虚部通过反傅里叶变换转化回空间域,获得目标图像(如背景图,异常图等)。

2 黄山-镜儿泉镁铁质超镁铁质岩带镍异常识别

2.1 地质背景

图1 黄山-镜儿泉地区地质略图(顾连兴等,2007)Fig.1 Geological sketch map of Huangshan-jing’erquan area,eastern Tianshan (after GU Lian-xing et al.,2007)

黄山-镜儿泉镍铜成矿带位于新疆东天山东部,吐哈盆地南侧,以甘新公路为界,与西段的康古尔地区共同组成了康古尔-镜儿泉构造岩浆岩带,如图1所示,据顾连兴的地质略图修改(顾连兴等,2007)。泥盆-石炭纪期间,觉罗塔格洋沿阿齐克库都克-沙泉子断裂向北俯冲形成弧后盆地,至石炭纪末期,这一带的大洋盆地已经完全碰撞闭合;从早二叠开始,该区进入后碰撞伸展阶段,发生了镁铁质超镁铁质岩和花岗岩类的就位,区内区域变质作用、混合岩化和韧性剪切活动强烈。区内构造主要包括三条近北东东向大断裂带,即北部的康古尔韧性剪切带和南部的雅满苏断裂带和阿奇克库都克大断裂带;区内地层主要为泥盆-石炭系火山岩,其中雅满苏断裂以北梧桐窝子组是一套海相喷发的基性熔岩,绿色-暗绿色细碧岩,干墩组是一套巨大的火山碎屑沉积岩和含炭硅质岩,镁铁质超镁铁质岩体主要侵位于该两套地层当中,断裂以南为雅满苏组,主要为一套双峰式的火山岩;此外,在阿奇克库都克断裂以南,主要分布着前寒武纪古老变质岩系,岩性为大理岩和变质砂岩;区内岩体主要为镁铁质超镁铁质岩类和花岗岩岩类,其中镁铁质超镁铁质类岩体主要分布在康古尔韧性剪切带附近,岩石类型有橄榄岩、辉橄岩、橄辉岩、二辉岩、辉长苏长岩以及辉长岩和闪长岩,岩体规模较小,最大者不超过10 km2,该类岩石中往往产出镍铜硫化物矿床;花岗岩类分布遍布全区,主要包括钙碱性-钾长钙碱性花岗岩、片麻状花岗岩和过铝花岗岩,其中钙碱性-钾长钙碱性花岗岩与斑岩型矿床有关,而过铝花岗岩往往与稀有元素伟晶岩型矿床有关(顾连兴等,2007)。

2.2 镍铜硫化物矿床

在黄山-镜儿泉成矿带上产出有岩浆型镍铜硫化物矿床、斑岩型铜钼矿床和伟晶岩型锂铍矿床,其中岩浆型铜镍硫化物矿床在该区乃至整个东天山地区占据重要的作用,目前已经勘探查明黄山和图拉尔根两个大型矿床、香山、黄山南、二红洼、葫芦、土墩和马蹄等中小型矿床。该类矿床的成矿时代与岩体近似,不同方法获得的年龄在288～269 Ma之间,表明该类矿床是早二叠纪造山后碰撞伸展阶段的产物(毛景文等,2002;秦克章等,2002)。矿床在区域上主要受康古尔-黄山弛张性深大断裂及其次级断裂控制,沿断裂上侵就位的镁铁质超镁铁质岩体是其主要含矿围岩,岩浆熔离作用是主要成矿方式,且各岩相分异程度越高,含矿性越好。主要矿石矿物为磁黄铁矿、镍黄铁矿、黄铜矿和黄铁矿等。

2.3 镍元素S-A技术分析

地球化学异常的空间结构较全面地反映了异常的特征。化探异常与其它类型地质异常相比,除具有化学分析含量数值精度高以外,还能较连续地反应异常的多尺度空间变化模式,因此局部结构模式的定量化有助于理解背景或异常的结构并提供研究地质过程、成矿过程以及异常识别的新线索。基于多重分形理论所定义的分形滤波技术便可以从自相似性的角度来进行地球化学场的分解和异常提取,本文采用中国地质大学(武汉)成秋明教授开发的GeoDAS 4.0软件来进行镍元素的S-A技术分析。

2.3.1 统计特征分析

本次使用的地球化学原始数据来自于新疆烟墩幅(样品性质为水系沉积物或岩屑,样品采集单位为新疆第六地质大队,分析单位为新疆地矿实验研究所,完成时间为 1988—1990年)和梧桐窝子泉幅(样品性质为水系沉积物或岩屑,样品采集单位为新疆物化探大队,分析单位为新疆地矿实验研究所,完成时间为 1987—1990年)1:20万地球化学测量,平均每4 km2含有一个化学分析样品,共计3064个化学分析样品,包括主量、微量等39个元素,其含量主要通过X荧光光谱分析测得。由于镍元素在本地区岩浆型镍铜硫化物矿床中是主要成矿元素,对于镁铁质超镁铁质岩体又具有成矿专属性,因此以镍元素为研究对象,研究镍铜硫化物矿床的成矿过程以及异常分解。首先使用传统的方法分析并了解镍元素数据的统计特征,然后通过分析3064个镍元素样品统计直方图和Q-Q图(图 2和图 3)可知,该区的大部分样品服从正态分布,小部分数据(右端截尾极高值)可能服从分形分布,但是由于高值区样品数据较少,因此在分离异常和背景方面,传统方法(如均值+2倍标准差)效果不是太好。

2.3.2 生成地球化学图

整个研究区 Ni元素浓度分布图通过反距离加权移动平均方法生成网格化数据来实现。窗口尺寸为0.5 km×0.5 km,逆距离滞后指数为2,最小样本数为12。生成的Ni元素空间分布图如图4所示。

2.3.3 生成功率谱图像及构建分形滤波器

图2 Ni元素浓度统计直方图Fig.2 Statistic histogram of Ni element concentrations

图3 Ni元素浓度Q-Q图Fig.3 Q-Q chart of Ni element concentrations

图4 Ni元素值/10-6地球化学图Fig.4 Geochemical map of Ni element values/10-6

通过使用傅里叶变换以及公式(4)可以计算出功率谱值S,功率谱值在Kx−Ky二维平面内的分布如图5所示:

在双对数纸上研究lgA(≥S)与lgS之间的关系,根据图6的线性特征,使用最小二乘法(LS),这些斜率不同的直线段反映了数据中的多重分形。

使用第三个交叉点(从左到右对应的第三条竖线)所对应的功率谱值 44433,根据需要生成背景和异常在频率域中所对应的不规则的滤波器,如图 7所示,其中功率谱值高于 44432的区域对应背景,功率谱值低于44433的区域对应异常。

2.3.4 构建地球化学背景和异常图

图5 功率谱/106分布图Fig.5 Distribution map of power spectra/106

图6 “面积”-功率谱双对数图Fig.6 lg-lg plot showing relationship between “areas” and power-spectra

图7 分形滤波器Fig.7 Filter constructed on the basis of distinct patterns

把使用S-A技术产生的滤波器应用到功率谱的实部R(Kx,Ky)和虚部I(Kx,Ky),然后使用反傅里叶变换将其转换到空间域。根据傅里叶变换的性质,图 8主要是原始区域的高频端元,这可能包括一些局部异常和噪音,图9主要是原始区域的低频部分,主要包括区域的背景。

2.4 地质解释及找矿方向分析

从 Ni元素原始地球化学分布图(图 4)上可以看出:镍铜硫化物矿床在空间上大部分对应着Ni的中等规模的化探异常或弱异常及其边部,甚至某些矿床位于背景区,如二红洼、土墩等,并且总体上异常呈带状分布,与康古尔深大断裂带相对应;而强度最大的地区位于研究区的西北方向,目前为止没有发现任何矿床,在地质图上大致与吐哈盆地相对应,主要被第三系和第四系盖层所覆盖;研究区的东南部(研究区的右下方),整体上处在一个元素含量较低的地区,仅有零星的弱异常出现,虽然到目前为止还没有发现矿床,但是在研究区外,阿齐克库都克大断裂附近已经发现多个镍铜矿床,如白石泉、天宇、天香等,该断裂在研究区内亦有分布。那么造成以上地球化学图分布格局的原因是什么？今后找矿应该到哪些地方才有可能有所突破呢？

大量东天山镍铜硫化物矿床的矿床学、矿床地球化学的研究结果普遍认为,石炭纪晚期至二叠世早期,东天山地区发生大规模的可能与地幔柱有关的强烈的镁铁质超镁铁质岩浆-火山活动(王登红等,2000),它由下地幔底侵上涌,喷溢到地表,可能是整个中亚地区地幔柱活动的一部分。由于受新特提斯洋关闭和青藏高原隆升影响,新生代以来天山-阿尔泰再次大幅隆升,原来可能大面积分布的早二叠世镁铁质火山岩已经被大量剥蚀,但其根部,抑或喷出岩的根部抑或潜火山岩的根部仍然保留着。当前这些沿多条区域性大断裂与镍铜硫化物矿床有关的东西向镁铁质超镁铁质岩体很可能是大面积喷溢熔岩的根部或者补给通道,东天山地区镍铜硫化物矿床的形成和保存模式如图 10所示(Mao et al.,2008)。而研究区西北部出现的强度很大的镍异常可能是二叠纪形成的溢流玄武岩在新生代被剥蚀,由于地形的原因,被运移到地势相对较低的吐哈盆地内,根据物质不灭原理,玄武岩中含有的大量的 Ni元素也被运移至此,因此造成了目前吐哈盆地内第四系沉积物中含有较高的Ni元素。而其它地区由于玄武岩已经被剥蚀或者镁铁质超镁铁质岩体面积相对较小,因此在区域上Ni元素含量较低或者仅仅为弱异常。

图8 Ni元素/10-6异常图Fig.8 Anomaly map of Ni element values/10-6

图9 Ni元素/10-6背景图Fig.9 Background map of Ni element values/10-6

图10 东天山地区后碰撞镍铜硫化物矿床形成和保存模式图(Mao et al.,2008)Fig.10 Model map showing metallogenic process and change of the post-collisional Ni-Cu sulfide deposits in East Tianshan (after Mao et al.,2008)

从根据多重分形滤波技术所获取的 Ni元素异常图(图 8)和背景图(图 9)上可以看出,康古尔大断裂带附近以及研究区东南部地区具有低背景强异常,且异常具有带状分布的特点,这主要是因为这些地区具有多条深大断裂分布,如康古尔塔格大断裂,雅满苏大断裂以及阿奇克库都克大断裂等,虽然玄武岩已经被剥蚀殆尽,但是其根部对应的镁铁质超镁铁质岩体还有可能大量存在(包括一部分隐伏的岩体),因此也是寻找镍铜硫化物矿床的有利地区(如图 8);研究区西北部,即吐哈盆地内,Ni元素呈现出高背景,弱异常的特点,且大部分异常比较分散,无明显的带状特征,这主要是由于该区处在吐哈盆地内,第四系内虽然含有大量的Ni元素,但大多数为风化剥蚀的产物,大都处于分散状态,难以形成规模矿床,因此不利于镍铜硫化物矿床的找寻。

3 结论

1)基于自相似原理,使用频率和空间信息的分形滤波技术(S-A)具有人为介入较少,能适应数据本身由于地质作用和成矿作用所形成的客观存在的各向异性及内部结构复杂性等特点,因此在场的分解和信息提取方面具有广泛的应用前景;

2)在吐哈盆地,Ni元素化学值最高,但经S-A技术分解之后,表现为背景高,异常弱且分散,这主要是因为第四系沉积物是二叠纪玄武岩被风化剥蚀运移至此的结果,因此成矿效果并不理想,不利于镍铜硫化物矿床的找寻;

3)在康古尔深大断裂带以及研究区东南部(尤其是阿奇克库都克断裂带附近),Ni元素化学值总体上中等或较低,但经 S-A技术分解之后,表现为背景低,异常强且呈带状分布,这主要是因为这些地区具有多条深大断裂分布,虽然玄武岩已经被剥蚀殆尽,但是其根部对应的镁铁质超镁铁质岩体还有可能大量存在(包括一部分隐伏的岩体),因此也是寻找镍铜硫化物矿床的有利地区。

致谢:本文在撰写过程中,阅读了中国地质大学(武汉)成秋明教授的大量文章,并学习使用了他的软件GeoDAS 4.0,受益匪浅,在此表示诚挚的感谢！