1. 引言
重磁地球物理数据具有深穿透性,且通常含有地壳深部关键地质结构和矿化信息 [1] [2]。而这些关键信息由于传统信息提取数学模型(如傅氏变换等)的局限性而没有得到有效地提取和应用,非线性技术为提取深部地质矿化信息,进一步精确构建深部地质结构和地质体模型提供了新的机遇 [3]。
地球科学的复杂性主要涉及地质过程的多期多阶段性及其结果的叠加性,记录这一复杂过程和叠加现象的数据集往往具有非线性结构和非平稳特征,使得诸如地质统计学和傅立叶变换等常用数据处理方法,严格意义上,并不适用于处理非线性和非平稳数据 [4] - [9]。为此,Huang等 [10] [11] [12] [13] 建立一种普适的信号分析技术,称为Hilbert-Huang变换(Hilbert-Huang Transform, HHT)。该方法包括两部分算法,即,经验模分解(Empirical mode decomposition, EMD)和Hilbert谱分析,二者有机结合定量刻画非线性和非平稳过程。该方法不依赖于核函数的先验选择,而是分解信号到起源于连续性极值的本征振荡模。该方法的创立与发展为基于非平稳数据研究致矿地质异常事件与过程开辟一个新领域。
经验模分解(EMD)可将复杂数据集分解为一系列有限个数的本征模态函数(intrinsic mode function, IMF),这一系列的IMF可以通过筛分过程(sifting process)得到,通过筛分过程还能够从原始数据中得到局部的相对高通滤波振动。该方法适用于地震数据处理 [14],噪声滤波 [15],航空重力数据一维处理 [16]。Nunes等 [17] [18] 发展二维EMD (Bi-EMD,简称BEMD)方法,使之适用于处理二维图像;沈滨等 [19] 将BEMD方法应用于纹理分割及图像瞬时频率估计。目前,该方法已被广泛应用于提取区域和局部重磁致矿异常和地球化学元素组合致矿异常信息提取,取得显著效果 [3] [20] [21] [22]。近年来,许多学者应用多维分型模型(multi-fractal mode)提取重磁异常亦取得明显进展 [23] [24] [25] [26]。
本文尝试应用BEMD方法,对鲁西地块1:20万重力数据进行分解,结合其地质矿化特征,试图揭示鲁西地块深部地质结构及其与金矿化分布的关系,为金矿深部成矿预测与评价提供科学依据。
2. 二维EMD基本原理及算法实现
EMD方法 [10] 能够自适应地将一组数据分解为不同频率的组成——本征模态函数(intrinsic mode function, IMF)。通过筛分,可以得到从高通滤波到低通滤波的一系列的IMF分量,而每一个IMF分量包含以下性质:1) IMF分量有相同的过零点数和极值点个数或者差1;2) 在任一点,由极大值点构成的上包络和由极小值点够成的下包络的平均值趋于0。筛分过程要保证数据的极值点个数至少为2。EMD分解可以将数据分解为有限个的IMF分量和剩余分量,不同的IMF分量代表了不同的振动模式(oscillatory mode),而剩余分量描述了数据的整体趋势。
2.1. 二维EMD的基本原理
在二维空间中,对一维EMD方法进行推广,实现了二维的EMD的分解(BEMD)。二维EMD方法过程与一维EMD类似,主要的区别在于二维EMD方法矩阵中极值点的确定和包络面的拟合比一维EMD方法更复杂(后面将详细说明)。令
为待分解的二维数据,通过二维筛分过程,可将二维数据分解为有限个的二维IMF分量(BIMF),分别代表二维数据的不同频率(尺度)的结构特征,按照频率的由高到低,依次为
,则有:
(1)
其中,
为第i个二维IMF分量,
为剩余分量。在滤波过程中,与一维类似,可以设计不同的滤波器
,
,
,分别用于高通、带通和低通滤波。也可以选择性的选取某些(个)反映特定频率(尺度)结构特征的二维IMF分量作为滤波结果。
(2)
(3)
(4)
换句话说,我们能通过长周期分量相加获得低通滤波结果,通过短周期分量相加获得高通滤波结果,通过选择性周期分量相加获取带通滤波结果。这就是该方法能够用于信号分解和致矿异常信息提取的方法学基础。
此外需设定筛分过程的停止条件,这里主要通过限定SD的大小来得到,SD由两个相邻的筛分过程的结果得到 [13] [17] [18],在二维中,设定:
(5)
SD的设定对BIMF有影响,如果SD的值较小,BIMF的个数将增加,筛选过程也更为细致,但运算时间会有所增加。但在二维中,
的取值大小还没有一个标准,往往通过经验来判断和设定,亦可借助局部均值矩阵
的期望值、方差获取。
2.2. 算法实现
二维EMD的流程步骤与一维EMD类似,主要是将一维的极值求取和包络线的拟合变为二维的极值求取与二维包络面的拟合。二维EMD的算法流程(图1)如下:
为待分解的二维数据,
1) 初始化:令
,并
(j为二维IMF的标号);
2) 提取第j个IMF分量
① 令
,
② 查找二维数据
的极大(小)值,
③ 计算上(下)包络
(
),
④ 计算均值
,
⑤ 令
,
⑥ 验证结束条件,
⑦ 循环执行步骤②~⑥,直到满足停止条件为止跳出,而此时的
即为第j个二维IMF分量
,
3) 剩余量
,
4) 令
,循环步骤2-4,直到剩余量
的极值数小于2为止。
3. 二维EMD方法在鲁西地块重力致矿异常信息提取中的应用
前已述及,重磁信息的一个显著特点是具有“透视性”,它不仅能够反映浅部的地质现象,而且通过对重磁场的分解,还能够获取深部地质结构信息。重磁信息的另一个特点是具有多解性,这是因为通常我们获取的重磁数据测量的是不同规模、不同深度、不同形态和不同磁性(密度)地质体组合的叠加场。这就要求我们根据地质体场的性质和特点,借助于信息处理技术实现叠加场的分解,尽可能使场与地质体一一对应,并结合地质矿化信息等约束条件,获取目标信息,最终达到解决地质与找矿疑难问题之目的 [27] [28]。
Figure 1. Flow chart of the empirical mode decomposition
图1. 经验模分解法流程图(根据Chen等,2017) [3]
3.1. 鲁西地块地质矿化特征
鲁西地区系指沂沭断裂以西的基岩隆起区。区内众多的金矿点和一些金矿床初步划分为三种成因类型:1) 变质热液金矿床;2) 与侵入岩有关的热液金矿床;3) 火山热液金矿床(图2)。与金矿成矿有关的主要因素如下:
a) 断裂:以NE向沂沭断裂及NW向断裂为主的断裂系统。在沂沭断裂带内或其附近的中生代火山岩盆地,一些金矿点受其控制呈NE向分布。在沂沭断裂带以西,NW向断裂系统控制了前寒武纪变质岩系、古生代碳酸盐岩沉积盆地、中生代火山岩盆地和中生代侵入体以及金矿床、矿点的分布。因此,NE向沂沭断裂带及NW向断裂系统构成了研究区的一级控矿地质因素 [28]。
b) 地层:本区前寒武纪结晶基底由表壳岩(5%)和深成侵入岩(95%)组成。表壳岩主要包括中太古代(>3 000 Ma)沂水岩群和新太古代(2800 Ma)泰山岩群。沂水岩群分布于沂水断裂带,其主要岩性为斜长二辉麻粒岩、二辉角闪斜长片麻岩和斜长角闪岩以及含紫苏磁铁石英岩等一套深变质岩系;其原岩为超镁铁质–镁铁质熔岩、凝灰岩及泥质砂岩夹硅铁质岩石。泰山岩群分布于沂沭断裂带以西的鲁西隆起区,其主要岩性为斜长角闪岩,透闪片岩、角闪黑云变粒岩夹铁闪磁铁石英岩;其原岩为超镁铁质–镁铁质熔岩(科马提岩)、凝灰岩、泥质粉砂岩和硅铁质岩;泰山岩群是我国保存最好、发育最完整的典型新太古代绿岩带 [29]。
c) 中生代岩体:与岩体有关的金矿床是本区金的主要工业矿床类型。与金矿有关的侵入体皆为中生代浅成斑状杂岩体,成岩时代为190~125 Ma [30] [31] [32],发育于隐伏基底(这里指被古生代沉积碳酸盐岩覆盖的前寒武结晶基底区)区。金矿床(矿点)通常环绕侵入体分布(图2)。这类金的工业矿床的形成归结于矿源(前寒武纪结晶基底),热源(中生代侵入岩),有利围岩(古生代沉积碳酸盐岩)和成矿构造等成矿要素的最佳匹配 [28]。
Ary:沂水岩群;Art:泰山岩群;TTG:前寒武纪变质深成杂岩;Pz:古生代沉积碳酸盐岩;Mz:中生代碎屑沉积岩。1:中生代侵入岩;2:中生代火山岩盆地;3:韧性剪切带;4:断裂;5:变质热液金矿床;6:侵入岩型热液金矿床;7:火山热液金矿床;8:城镇
Figure 2. Regional outline map of gold mineral resources in the western Shandong uplift block
图2. 鲁西地块地质矿产图(据Chen et al, 2001修编) [28]
3.2. 鲁西地块重力致矿异常信息提取
论文研究所用重力数据来源于山东地质矿产局第物化探大队1980年代实施的1:20万地面综合物探调查中的重力测量结果。其测量网距为2 km × 2 km,工作总精度为±2.32 g.μ,控制面积约48,000 km2。鲁西地块出露地层密度参数由高至低依次为:泰山岩群(2.70~2.90 g/cm3,平均值2.81 g/cm3)→寒武–奥陶系(2.64~2.73 g/cm3,平均值2.70 g/cm3)→石炭–二叠系(2.58~2.67 g/cm3,平均值2.61 g/cm3)→侏罗–白垩系(2.39~2.54 g/cm3,平均值2.51 g/cm3)。侵入岩密度参数由高至低依次为:前寒武纪花岗质侵入岩(TTG岩系:2.59~2.84 g/cm3,平均值2.75 g/cm3)→中生代花岗岩类(2.59~2.71 g/cm3,平均值2.65 g/cm3) [20]。
原始重力数据图像(图3)是不同深度、不同密度地质单元及其地质结构的综合反映。为揭示鲁西地块深部地质结构及其对金矿化系列形成与分布的控制,对其1:20万重力数据(图3)进行二维EMD分解,获得4个二维IMF分量,以及剩余量
,
。各个二维IMF分量代表二维数据的不同频率的结构特征。BIMF1、BIMF2、BIMF3和BIMF4按照频率从高到低产生的,对于相同的局部区域BIMF1带通滤波频率通常高于BIMF2的频率。在鲁西1:20万重力数据的处理过程中,设定SD = 0.02。
Figure 3. The original gravity data image compiled by the gravity data surveyed at scale 1:200000 (the original gravity data from the second geological survey team of Shandong bureau of geology & mineral resource, 1991). The legends for geology and the gold deposits see in Fig.3, AB and CD lines mark respectively two one dimension gravity sections at perpendicular to each other for orthogonality assessment of IMFs
图3. 鲁西地块原始重力数据图像(数据来源山东地矿局物化探大队,1991),地质图例同图3,AB和CD线段是对分解的不同的IMF重力分量进行正交性检验的一维重力剖面位置
3.3. IMF分量正交性检验
一维IMF满足正交性,虽然还不能在理论上进行证明,但实验结果证明一维IMF分量是正交的 [10]。假设各个一维IMF分量是正交的,则:
,(
) (6)
中等式右边的第二项等于0 (或趋近于0),其中
,
,
(剩余分量)。在实际计算中,用IO来表示正交性指标,
(7)
或者可以计算两个不同IMF分量之间的
值:
(8)
Figure 4. The IMF components of one dimensional gravity data in NE trend of section (AB)
图4. 鲁西地块北东向(AB)一维重力数据分解的各个IMF分量
Figure 5. The IMF components of one dimensional gravity data in NW trend of section (CD)
图5. 鲁西地块北西向(CD)一维重力数据分解的各个IMF分量
这里,T是总的时间区间,为获得IMF的最优分解,IO应该非常小,典型地,其值介于0.01和0.001之间是可以接受的。如果信号长度很短,IO的值会增大至0.05左右。线性分离系统需满足正交性,对于一些特殊的信号,Huang [10] 认为虽然分量之间不再满足正交性,但由EMD分解到的IMF分量依然是有意义的。分别详细研究图4中北东向(AB)和北西向(CD) 2个方向的二维IMF分量的一维数据结果(图4,图5)。
并对这两个方向的IMF分量进行正交性检验(表1),结果表明两个方向的IMF分量近似满足正交性。
Table 1. Orthogonality assessment of IMFs
表1. IMF重力数据分量的正交性检验
4. 结果与讨论
在上述分解的4个分量中,其中IMF4 (图6)和IMF3 (图7)具有较明确地质意义。
Figure 6. BIMF4 image for the western Shandong uplift block, the legends for the gold deposits and their geology can be seen in Fig.3
图6. 鲁西地块IMF4重力分量图像,地质矿产图例同图3
Figure 7. IMF3 image decomposed from t original gravity data for the western Shandong uplift block, the legends for the gold deposits and their geology see in Fig. 3
图7. 鲁西地块BIMF3重力分量图像,地质矿产图例见图3
4.1. IMF4和IMF3重力分量图像的地质意义
图6是低通滤波重力分量(IMF4)图像。该图表明,该重力分量值的高低分布与地表岩石的密度没有对应关系,反映了鲁西地块深部区域性地幔隆起和凹陷。在该图像上,沿沂沭断裂呈北东向分布的高重力分量异常带可能反映沿临沂—沂水分布一地幔隆(或称基底隆起)起,简称临沂幔隆(IUY),重力IMF4分量值分布−40~80 μm。沂沭断裂带基本上沿该幔隆分布,其中火山岩盆地和幔隆两侧的金矿床亦呈北东向分布(图6)。
分布于沂沭断裂以西的IMF4低重力分量异常区反映了地幔凹陷特征。该区存在两处低重力分量异常区,分别对应两处地幔凹陷。其一是分布在平邑以南呈北西向分布的凹陷,简称平邑幔陷(IDP),重力IMF4分量值分布−310~−190 μm,对应于该凹陷地表出露的地质单元是由太古代TTG和少量变质岩系构成的泰山岩群。在平邑凹陷北西侧分布有呈北西向展布的平邑–费县中生代火山岩盆地、燕山期鲁西杂岩体以及相关的大型归来庄金矿床等 [33] [34],在其南东侧分布有苍山–龙宝山花岗岩带以及相关的金矿床。其二是分布于蒙阴–新泰断裂以北以莱芜为中心呈弧(扇)形分布的地幔凹陷,简称莱芜幔陷(IDL);其重力IMF4分量值分布−310~−190 μm。对应地表出露的地质单元内弧为太古代TTG岩系 [35],外弧为古生代碳酸盐岩。幔陷内自北向南发育三条呈弧形分布的中生代岩浆岩(侵入岩和火山岩)带,且不同程度的发育金矿化(图7)。此外,垂直于弧形幔陷的轴线发育放射性断裂。
总之,幔隆和幔陷从整体上控制研究区基本地质格局以及中生代与火成岩有关的金矿床和相关火成岩(侵入岩和火山岩)体的区域分布。
图7是带通滤波重力分量(IMF3)图像。该图表明,该重力分量值的高低分布与地表岩石的密度没有对应关系,可能反映在一级区域性地幔隆陷基础上发育二级区域性地幔隆陷 [3] [36]。在整个研究区这种二级区域性幔隆和幔陷呈带状相间分布。在研究区东部,临沂–沂南–沂水幔隆西侧分布有相同展布方向和近似规模的孟良崮–棱背岭–沂山–唐吾幔陷,在其南东侧发育临沭–莒南幔陷。在研究区中西部,存在4对呈北西向分布的地幔隆陷。自南至北依次为张庄–滕州幔陷;向城–石井–水泉–南辛幔隆;罗庄–郑城–四海山幔陷;临沂–费县–平邑–泗水幔隆;蒙山–宫里–泰安幔陷;蒙阴–新泰–徂徕–柳埠–济南幔隆;雁翎关–跺庄幔陷;辛庄–腰关–阎家峪幔隆。构造线的分布与幔隆和幔陷的展布方向具有高度一致性。绝大多数中生代火山岩及其相关的金矿床和铜金矿床分布于幔隆区,而绝大多数中生代侵入岩及其相关的金矿床和铜金矿床则分布于幔隆的边缘。
4.2. 金矿化及其相关火成岩形成的构造环境
一些学者 [37] [38] [39] [40] 的研究资料表明:中生代太平洋板块向欧亚板块俯冲,导致鲁西地区的整体隆起;伴随郯庐断裂大规模的左行平移,在其西侧形成了一系列NW-NNW向断裂和陆内断陷盆地,在以张应力为主的构造环境下,出现了大量火山喷发和浅成侵入活动,从而形成了一系列火山岩盆地,和浅成次火山斑状杂岩体。这种发育中生代火成岩的以张应力为主的构造环境恰是一级幔陷中的局部隆起所致。王先美等 [41] [42] 由出露于沂沭断裂带、鲁西地体、鲁东地体的中生代地层、岩浆岩,结合断裂活动年代学、区域地质等资料分析,将沂沭断裂带晚中生代构造演化划分出距今约160 Ma、130~110 Ma、90~80 Ma等3个关键时期,并分别与左行压剪、左行张剪、右行压剪构造活动相对应。晚中生代北北东走向的沂沭断裂带与北西向断裂系的共轭属性及其构造活动方式,是在华北地块周边不同板(地)块、不同时期的汇聚与离散等相对运动,以及地球内部不同圈层的物质与能量交换的背景下发生的。沂沭断裂带早期左行压剪与北西向断裂系的右行压剪构造活动的动力来自华北、华南地块的陆–陆碰撞;沂沭断裂带早白垩世左行张剪运动与北西向断裂系同期右行张剪活动是Izanagi板块俯冲与华北岩石圈拆沉引起上地壳快速伸展的联合效应;Kula板块晚白垩世向西的快速俯冲可能构成华北地块挤压隆升的主要动力,沂沭断裂带与北西向断裂系则相应表现出右行压剪与左行压剪构造活动。
某些学者 [43] [44] [45] [46] 根据同位素定年,微量元素,稳定同位素和稀土元素等分别对鲁西铜石、铜井、莱芜、龙宝山、莲子旺、粟园等杂岩体进行了较系统研究。其结果表明,其成岩时代分布112.5~189.8 Ma,众数值是110~126 Ma,这一时期与前述的沂沭断裂左行张剪相对应,亦与早白垩世中国东部岩浆、成矿等作用最为强烈的时期(120~130 Ma)相吻合 [47] [48] [49] [50] [51],亦很可能对应岩石圈减薄的高峰时期 [52]。
林景仟等 [44] 通过对鲁西地区中生代火成活动的研究认为:鲁西地区的火成活动远与太平洋板块对欧亚板块的俯冲导致沂沭断裂的活动有密切联系。此时,在鲁西地块上北西向的断裂频频活动,该组断裂的强烈拉张最早发生在早侏罗世,发育于鲁西的南部,生成了高钾钙碱性–高钾碱性的岩石组合(如铜石杂岩体),至早白垩世地壳拉张活动加剧,南部区域的北西向断裂切割深度大,从上地幔的深层生成了碱性岩浆,北部区域北西向断裂亦切割到上地幔,生成了钙碱性的辉长质岩浆,它与上部壳层产生的熔体发生了混染,形成了二辉闪长质岩浆(较高水条件下生成了角闪闪长质岩浆),二长质岩浆是在混入了更多的酸性熔融体的情况下产生的。在此时期,近沂沭断裂带的地段,早白垩世较早阶段生成的辉长质岩浆发生混染生成闪长质岩浆,上部的壳层则产生了花岗质岩浆,后一组岩石多在沂沭断裂带近旁分布,它们亦沿北西向断裂活动,分布于较远离沂沭断裂带的部位,如滕州桑村岩体、蒙阴的虎头崖岩体等。
形成杂岩体的岩浆最初可能起源于上地幔,在其上升过程中不同程度地同化了太古代基底岩层,使部分Au、Cu、F、B等物质进入岩浆。随着岩浆的进一步侵位,冷凝结晶(局部伴有交代作用),成矿组分(Au、Cu等)和矿化剂(F、B等)在富碱质的热液中进一步富集形成液态矿源,并在岩体内外有利的构造部位沉淀成矿,最终形成火成岩体(侵入岩和火山岩)和金与金–铜矿化 [53] [54] [55]。
鲁西在中生代除发育花岗质岩石及其相关的金(铜)矿床外,还发育基性火成岩 [56]。典型玄武岩发现于费县方城一带,K-Ar年龄为119~125 Ma [56] [57]。鲁西中生代火山岩的主体为钾玄质火山岩,蒙阴盆地钾玄质火山岩的40Ar-39Ar年龄变化于114.8~124.3 Ma之间 [58]。谭东娟和林景仟 [59] 测得济南辉长质杂岩体的40Ar-39Ar年龄为115 Ma,山东省地质矿产局(1991)测得邹平盆地侵入青山组火山岩的辉长质岩体 40Ar-39Ar年龄为128~130 Ma [60]。它们皆属于燕山晚期岩浆活动的产物。
上述研究表明:中生代太平洋板块向欧亚板块俯冲导致的郯庐断裂大规模的左行平移引起鲁西地幔隆陷,地幔隆陷导致的拉张和挤压作用引起的构造岩浆活动是中生代金、铜成矿的有利背景,亦是该类矿床找矿的远景地段。
5. 结论
二维经验模分解(BEMD)方法被成功地应用于鲁西地块重力数据分解,获取反映鲁西地块不同尺度深部地幔隆陷及其与金矿化空间分布关系的IMF4和IMF3重力分量图像。
1) IMF4低通滤波重力分量图像显示鲁西地块存在1处一级区域性地幔隆起,临沂幔隆(I)和两处一级区域性地幔凹陷,平邑幔陷(IDP)和莱芜幔陷(IDL)。几乎所有金矿床及其相关的中生代火成岩(侵入岩和火山岩)皆沿幔隆或幔陷的边缘分布。此外,研究区北部近于平行分布的三个弧形花岗岩带的形成与分布受莱芜弧形幔陷的控制,且不同程度的发育金矿化。
2) IMF3带通滤波重力分量图像,揭示在一级区域性地幔隆陷的基础上发育的二级区域性地幔隆陷特征。在整个研究区这种二级区域性幔隆和幔陷呈带状相间分布。在研究区东部,临沂–沂南–沂水幔隆西侧分布有相同展布方向和近似规模的孟良崮–棱背岭–沂山–唐吾幔陷,在其南东侧发育临沭–莒南幔陷。在研究区中西部,分布4对呈北西向分布的地幔隆陷。自南至北,依次为张庄–滕州幔陷;向城–石井–水泉–南辛幔隆;罗庄–郑城–四海山幔陷;临沂–费县–平邑–泗水幔隆;蒙山–宫里–泰安幔陷;蒙阴–新泰–徂徕–柳埠–济南幔隆;雁翎关–跺庄幔陷;辛庄–腰关–阎家峪幔隆。研究区构造线的分布与二级幔隆和幔陷的展布方向具有高度一致性。金矿床及其相关的中生代岩浆岩(侵入岩和火山岩)绝大多数分布于二级幔隆区或其边缘。
3) 中生代太平洋板块向欧亚板块俯冲导致的沂沭断裂大规模的左行平移导致鲁西地块在一级地幔凹陷的基础上发育二级区域性地幔隆陷。张应力为主的构造环境中形成的中生代火成岩恰是二级幔隆所致;深部地质结构对金矿化的形成与分布的控制是通过二级(局部)地幔隆起对与金–铜矿化有关的中生代火成岩的控制而实现的。因此,发育于二级地幔隆起及其边缘的中酸性岩浆岩区是找寻金–铜矿床的远景地段。
基金项目
国家自然科学基金(批准号:41972312,41672329),国家重点研发计划课题(批准号:2016YFC0600509)资助。
NOTES
*第一作者。
#通讯作者。