基于有限元数值模拟的直流电阻率激电找矿探测技术

马德锡, 谭捍东, 张志勇, 陈亮, 皮进军

(1.中国地质调查局地球物理调查中心, 廊坊 065000; 2.中国地质大学地球物理与信息技术学院, 北京 100083; 3.东华理工大学地球物理与测控技术学院, 南昌 330013)

在金属非金属矿产勘查方面,直流电法是应用最广泛的物探方法,尤其是激发极化法的作用更加显著,据《中国矿床发现史·物化探卷》统计,中国早期发现的贵金属和多金属矿床中,激发极化法是最有效的方法,适用于热液型、块状硫化物型、细脉浸染型、石英脉型等多种类型矿床的探测,至近年,国内外诸多学者还在持续开展激发极化法找矿研究及应用,刘明等[1]基于Dias 模型进行一维层状介质正演,研究天然场源激电异常参数。涂永明[2]研究了激电异常与金多金属矿化的关系。郝海强等[3]、赵伟峰等[4]、赵林林[5]、蒙凯等[6]研究了激发极化法在铜、锌、锰、银、钼、钨等多金属矿找矿中高极化异常与矿化的对应关系及找矿效果,然而金属矿地质条件的复杂性和勘查任务的多样性,使得地质目标体的异常解释变得困难[7-10]。马德锡等[11]、侯宇健[12]、潘北斗[13]利用多种方法综合解释、电阻率极化率联合反演、精细反演等方法解释异常,从而降低多解性,为了更准确识别不同异常所反映的地质目标体,并有效指导工程验证,20世纪70—90年代初,中国部分地球物理学者开展了水槽、土槽、薄水层、电阻网络及导电纸等物理模拟试验研究[14-16],取得了一定的认识成果,极大推动了中国直流电阻率激电方法进步及在找矿方面的应用,但不规则的地质体往往难以构建物理模型并进行物理模拟,而且物理模拟实验也无法满足野外无限半空间的地质条件。90年代末计算机技术的发展极大地推动了数值模拟技术的进步,近年来数值模拟由于成本低、效率高,在许多行业得到了广泛应用[17-20],尤其是电法勘探研究领域,有限元数值模拟以其独特的优势成为研究地质体异常的主要手段,诸多地质地球物理学者从不同方面进行了研究,孔重阳等[21]、戴前伟等[22]、彭艳华等[23]、李晓斌等[24]分别从稳定电流场满足的微分方程、非网格有限元方法、长电极线和电源不同模式等出发研究数值模拟对于成果解释的应用,但是针对电法找矿探测地质解释的异常精确判别、定量解释方面研究的系统性不够、代表性不强,尤其是基于有限元数值模拟,按照反演异常解释目标体,构建目标体模型通过正演求取异常响应,从而与反演异常对比来分析并提高反演解释准确性和定量化方面相比较少[25-28],基于此,现系统分析直流电阻率激电找矿探测中地质目标体的几种基本情况,重点用长方体低阻高极化模型、直立板状体低阻高极化模型、首层低阻弱极化模型、首层高阻弱极化模型、首层不连续高阻弱极化模型、正地形(凸)模型、负地形(凹)模型来代表独立规则矿(化)体、断层构造、低阻屏蔽、高阻盖层、地表不均匀体、山峰、山谷等进行有限元数值模拟及异常响应对比研究,从而指导直流电阻率激电找矿勘查实践,并将有限元数值模拟引向电法找矿探测生产实践。

1 直流电阻率激电数值模拟研究基础

代表独立矿化体的低阻高极化体、代表断层构造(矿脉)的低阻高极化板状体、低阻屏蔽、高阻盖层及复杂地形条件影响等在找矿实践中较为常见,也是直流激电找矿探测需要重点分析和解释的问题,图1列举了几种常见情况的电阻率激电找矿探测断面实例,图1(a)中,低阻高极化实例是青海一个矿区探测结果反演断面图,断面中部低阻和高极化率组合异常推测可能是矿化体;图1(b)中,低阻高极化实例是新疆一个矿区探测结果反演断面图,倾伏低阻带对应有高极化率异常;图1(c)中,断层构造实例是新疆一个矿区探测成果反演断面图,近直立的低阻带,推测为低阻断裂,同时在低阻带边部有高值极化率异常,推测为矿化断裂;图1(d)中,低阻屏蔽是内蒙古一个矿区探测成果反演断面图,显示低阻的沙层覆盖造成严重的低阻屏蔽,电阻率和充电均无明显异常显示;图1(e)中,正(凸)地形实例是内蒙古另一矿区探测成果反演断面图,显示剖面中部地形高差超过70 m,地形给深度异常产生明显影响;图1(f)中,负(凹)地形实例是青海一个矿区探测成果反演断面图,地形变化剧烈,主要的山谷负地形落差超过60 m,电阻率低阻异常形态因受地形严重影响,地质解释有一定的难度。

图1 直流电阻率激电探测异常典型实例

数值模拟在工程、物理、工业等领域解决不适定问题方面有广泛的应用,电法探测及反演解释本质也是基于先验信息不断寻优解的过程,基于探测信息的反演也具有极大的不适定性,作为有限元数值模拟其根本性数学基础没有变,直流电阻率激电有限元数值模拟其物理基础满足有限元数值模拟条件[29-30]。直流电法电场分布在地下无限半空间中,地下电场电位以供电点为中心向外逐渐变小,根据点源电场直流电阻率问题,电位及电流在地下满足如下边值方程[31-32]。

(1)

式(1)中:σ为电导率;φ为电位;I为供电电流;δ(·)为狄拉克函数;A为供电点;Ω为研究区;Γs和Γ∞分别为Ω的地表边界和地下无穷远边界;n为边界Γs外法线方向;r为源点到计算点的距离。研究点电源供电的二维地电断面,正演算法主要是求解傅里叶变换域中的二维偏微分方程,点源二维边值方程如式(2)。

(2)

式(2)中:k为波数;K0、K1分别为零阶、一阶第二类贝塞尔函数;U为傅氏域电位。因此,解析方程也表明电场在一定范围之外的电流密度很小对研究的问题影响可忽略[33],正演计算附加狄利克雷边界条件(Dirichlet)对正演模型边界进行约束[34-35],因此用有限半空间加边界条件来代替无限半空间,即用有限网格单元代替无限半空间并通过数值计算的方法将连续分布的电场计算为离散网格节点的电场值问题[36]。因为有限元法在“同样”网格离散化的情况下可以产生更精确的正演模拟解,而且模拟的地质体异常形态更好[37],因此选择用有限元数值模拟方法研究直流电法勘探的金属矿异常特征,正演计算采用Cholesky分解(CD),对于多电极的情况(>20),该方法的计算是稳健和稳定的,而且占用内存小[38],因此基于有限元数值模拟方法对直流电阻率激电找矿探测技术进行研究。

电法勘探的主要依据是目标体与周围介质的电性差异,电性主要是导电性和激发极化特性,这也是直流电阻率法和激发极化法用于找矿探测的基础,通常用电阻率和极化率(充电率)来表征这两种性质,这两种属性参数也是直流电法找矿勘查的基础参数指标,具体计算公式如下。

(3)

(4)

(5)

式中:ρs为视电阻率;ηs为极化率;Ms为充电率;k为装置系数;ΔU为电位差;ΔU2(t)为二次电位差;ΔU1+2为总电位差;t为供电时间。一般金属矿体因其含有金属矿物而表现为良导特性,相对于围岩表现为低电阻率,同时金属矿物由于激发极化效应而表现为高极化率(充电率),所以本文有限元数值模拟中,在构建先验模型时,拟用低阻高极化属性模型代表电阻率激电找矿探测中矿(化)目标体。

2 规则模型有限元数值模拟及异常响应

2.1 独立矿(化)体模型有限元数值模拟

有限元数值模拟在离散网格构建方面,即要保证建模精度同时还不能使CPU占用时间太长,因此横向网格设置为两个电极之间为2个分区,网格剖分单元为矩形网格,网格单元宽度相等,厚度增量系数、深度系数均为1.1,网格总数110×13。独立矿化体模型为长方形,背景电阻率设为700 Ω·m、充电率设为10 mV/m、矿化体模型电阻率设为100 Ω·m、充电率设为200 mV/m,矿化体模型相对背景的电性特征为低阻高极化特征,具体独立矿体的二维有限元剖分网格断面和异常体模型如图2(a)和图2(b)所示。

图2 相对围岩为低阻高极化的独立矿(化)体模型及数值模拟异常响应断面图

图2(c)和图2(d)中独立矿化体模型数值模拟异常响应断面图中,相对周围介质为低阻高极化的独立矿化体模型对应的是低阻高极化异常,电阻率和充电率异常均表现为椭圆形,椭圆长轴方向与长方形模型的长轴方向一致,异常宽度与模型一致,而厚度与模型差异较大,充电率异常相比电阻率异常向深部偏移较多。数值模拟结果说明直流电法对独立的规则矿化体的宽度及上表面(或顶板)反映准确,而对于矿化体下界面(底板)反映不准确。故在反演解释时,根据异常推断矿化体宽度、埋深、延深等空间属性时,实际矿化体下部延深(下边界)推断不能单独依据异常下边界,需要综合分析电阻率激电异常的叠加关系,并用其他资料提供约束,否则会得出错误的推断解释结果。

2.2 直立板状断层(矿脉)模型数值模拟

断层使原始地层及岩石结构破坏,常产生较多的次级断裂和裂隙,这些断裂不仅成为矿物运移和储存的通道,且常因含水或充填物而具有低电阻率高极化率的地球物理特征,在找矿勘查中是重要的地质和地球物理找矿标志。另外,水相比砂砾、沙、粉砂岩等电阻率低,且随着矿化度的提高,水和岩石的电阻率均会下降而趋于零[39],所以含矿(矿化)断层其电阻率更低,同时能够引起高的极化率异常,故用低阻高极化的直立板体模型代表断层(矿脉)做数值模拟,构建初始模型如图3(a)和图3(b)所示,背景电阻率为700 Ω·m,充电率为20 mV/m,断层模型电阻率为20 Ω·m,充电率为150 mV/m[40-41]。

图3 相对周围介质为低阻高极化的陡立断层(矿脉)模型及数值模拟异常响应断面图

直立板状断层(矿脉)数值模拟异常响应如图3(c)和图3(d)所示,电阻率和充电率异常形态均和模型相似,浅部异常宽度和模型一致,深部逐渐变宽呈等腰三角形;14 m以下,电阻率异常值比模型电阻率值高,充电率异常值和模型充电率值一致,由数值模拟异常响应形态和值的变化特征得出,应用直流电阻率激电对断层(矿脉)进行勘探时,成果解释需要综合分析异常整体结构特点、梯度变化趋势和异常中心位置,依据电阻率异常确定断层(矿脉)中心位置、产状及宽度等形态属性,依据充电率确定断层与矿化的产出关系及和矿脉的含矿性,必要时可利用数值模拟辅助进行解释推断。

3 特殊地质条件有限元数值模拟及异常响应

3.1 “低阻屏蔽”数值模拟

沼泽湿地、盐碱地、松散沙层及刚下过雨的地表会形成良导层,位移电流较大,易形成地表电流密度层,有时还会增加直流电法的电磁耦合效应,使电阻率数据发生畸变,这种效应称为“低阻屏蔽”,为此采用有限元数值模拟对其进行研究。假设背景电阻率为700 Ω·m、充电率为10 mV/m;设良导层电阻率为20 Ω·m、充电率为50 mV/m,深部矿化体模型的电阻率为100 Ω·m、充电率为200 mV/m,即矿化体模型相比围岩为低阻高极化体,构建正演模型如图4(a)和图4(b)所示。“低阻屏蔽”数值模拟异常响应如图4(c)和图4(d)所示,从电阻率异常断面[图4(c)]可看出低阻盖层下的低阻体没有显示异常,图4(d)显示矿化体模型有对应的异常,但异常发生畸变,形态范围与实际模型相差较大,数值模拟研究表明地表良导层对深部矿化体电阻率影响严重,而其弱极化效应对矿化体的充电率影响较小。所以在湿地沼泽区、盐碱地、沙层覆盖区、下过雨等情况使用电阻率法探测没有效果,但单独开展激发极化法开展找矿探测是可行的[42-43]。

图4 “低阻屏蔽”地质条件矿化体模型及数值模拟异常响应断面图

3.2 高阻盖层数值模拟

地表灰岩风化层、基岩出露及土壤腐殖质等会形成高阻层或局部高阻层,这又使得直流电法电流损耗加剧,影响深部数据采集。构建如图5(a)和图5(b)的高阻盖层及矿化体模型,设盖层电阻率为2 000 Ω·m、充电率为10 mV/m,围岩背景电阻率为700 Ω·m、充电率为50 mV/m,矿化体模型电阻率为100 Ω·m、充电率为200 mV/m,即矿化体模型相比围岩为低阻高极化体,以此模型采用有限元数值模拟研究直流电法中地表高阻盖层对深部矿化体的影响[44-45]。数值模拟异常响应如图5(c)和图5(d)所示,异常断面图显示均匀的高阻盖层对深部矿化体影响弱,电阻率和充电率异常的形态范围均和实际模型接近。

图5 高阻盖层和矿化体模型及数值模拟异常响应断面图

图6为地表不均匀局部高阻盖层模型及数值模拟异常响应断面图,图6(d)中矿化体模型异常响应发生畸变,异常形态和产状与模型差异较大,而且充电率异常向深部被放大,所以地表不均匀分布的局部高阻层会给直流电阻率激电带来明显影响。

图6 不均匀局部高阻盖层和矿化体模型及数值模拟异常响应断面图

3.3 地形起伏条件数值模拟

为研究起伏地形对电阻率激电野外探测的影响,构建如图7(a)和图7(b)所示的山峰地形条件和低阻矿化体模型断面及图8(a)和图8(b)山谷地形条件和低阻矿化体模型断面[46-47]。山峰高差约20 m,山谷高差约12 m,数值模拟异常响应断面如图7(c)和图7(d)所示。通过对地形起伏条件有限元数值模拟分析得出,地形对电阻率和充电率异常均有一定程度影响,使得矿化目标体异常发生倾斜,尤其是充电率异常较为明显,异常倾伏具有背离山峰、朝向山谷的特征,所以在复杂地形区依据直流电阻率激电异常进行矿化体地质解释时,单纯依据异常形态容易误判矿化体的产状,需要根据山峰和山谷对矿化体的影响及异常体上的规律,对分布位置及产状等做必要的矫正,从而得出更加准确的地质推断解释结果。

图7 山峰地形和矿化体模型及有限元数值模拟异常响应断面图

图8 山谷地形和矿化体模型及有限元数值模拟异常响应断面图

4 实测数据反演异常及数值模拟异常响应对比分析

在直流电法找矿探测中,实测数据反演异常地质解释常难以定量,尤其是地表没有指示标志的深部异常定性都很困难。图9(a)和图9(b)为内蒙古一个金矿区直流电阻率激电找矿探测反演断面图,图9(c)和图9(d)为地质解释断面图。电阻率断面[图9(a)]主要呈现高阻(>700 Ω·m),浅部存在相对高阻岩体(>2 000 Ω·m),在断面横向300～310 m有电阻率低阻带(<120 Ω·m),西倾近直立,推断为断层;充电率断面[图9(b)]三处异常对应在电阻率断面的低阻与高阻过渡带附近,强度大(>260 mV/m)套合好的充电率异常分布在断层与上盘围岩接触带位置,据充电率异常特点及与低阻异常的位置关系推断异常由矿化体引起。

图9 实测数据反演电阻率极化率异常断面图及反演解释断面图

由此,以实测数据反演异常并结合地质解释推断开展有限元数值模拟及异常响应分析,首先基于电性结构断面图构建有限元数值模拟模型如图10(a)和图10(b),电阻率模型断面[图10(a)]中的红色部分模拟高阻岩体,电阻率值在5 000～10 000 Ω·m,浅蓝色模拟该低阻断裂,电阻率为120 Ω·m。充电率模型断面[图10(b)]中有4个独立的矿化体模型,模型中心充电率值预设为350 mV/m,边部设置为150 mV/m。岩体、断层、矿化体模型数值模拟异常响应如图10(c)和图10(d)所示,电阻率断面[图10(c)]中,对比反演解释异常断面[图9(c)],有限元数值模拟后浅部高阻不均匀岩体(模型)没有得到清晰的异常反映,低阻断层有较明显异常,产状及中心位置等与反演解释(模型)较一致,深部宽度有放大。充电率异常断面[图10(d)]三处异常与反演解释异常(模型)位置对应,但与反演解释矿化体形态不完全相同,没有出现虚假异常[48-49]。通过有限元数值模拟异常响应与反演解释异常地质体对比分析,结果表明直流电阻率激电测量结果所推断的不均匀岩体、断层和矿化地质体都是存在的,但属性方面和数值模拟先验模型有部分差异,所以借助有限元数值模拟方法,可以对直流电法成果推断解释的地质异常体进行验证,对难于定性的尤其是隐伏的地质体推断解释具有辅助作用[50-52]。

图10 实测数据反演异常的数值模拟模型及数值模拟异常响应断面图

5 结论

(1)有限元数值模拟具有比物理模拟更强大的实现能力,且高效方便;基于有限元数值模拟技术,可对直流电阻率激电找矿探测异常解释中构建任意形状地质目标体模型,通过异常响应对比分析不断先验模型,从而得出使地质推断解释结果更接近地质目标体实际的位置、形态及属性,有助于定量解释和三维地质建模。

(2)地质体模型有限元数值模拟异常响应与实际地质体属性特征具有一致性;因此可以根据实测数据反演异常灵活构建各种地质体模型,用有限元数值模拟进行验证,从而实现对探测目标体定性定量的解释。

(3)规则矿化体模型有限元数值模拟分析表明,直流电阻率激电反演结果中规则异常是规则的地质体引起的;依据低阻高极化的组合异常可以推断解释矿体空间位置和顶板埋深,但底板深度或底部延深不易判断。

(4)直流电阻率激电反演结果中陡倾的连续的条带状低阻高极化组合异常,通常是由控矿断层或者矿脉所引起的,如果异常在深部变宽,可能与反演初始模型有关,是否代表断层或者矿脉在深部规模变化,需要仔细研究异常等值线的变化趋势来确定。

(5)地表良导层、高阻盖层、高阻不均匀体、地形起伏变化等条件会使深部目标体电阻率、极化率异常发生畸变;尤其是地表存在的良导层会掩盖下部探测目标的电阻率异常,但激发极化法受低阻层影响较小。高阻盖层不会对深部目标异常造成影响;模型异常响应与探测目标体位置一致,但形态并不完全相同。深部矿化体受起伏地形影响造成异常假性产状,通常倾向山峰、背离山谷。