地球物理学报  2019, Vol. 62 Issue (10): 3950-3963   PDF    
三维激发极化法反演研究及其在金属矿勘探中的应用
张刚1,2,3, 吕庆田2,3     
1. 中国地质大学(北京)地球物理与信息技术学院, 北京 100083;
2. 中国地质科学院地球深部探测中心, 北京 100037;
3. 中国地质科学院地球物理地球化学勘查研究所, 河北廊坊 065000
摘要:激发极化法在金属矿、硫化矿等资源勘探方面应用较广.随着勘探设备与计算机硬件的发展,所采集到的数据具有观测数据量大、布极方式多样等特点.针对实测数据的特点,我们研究完成了基于并行技术的激发极化法对数反演算法,该算法具有如下特点:(1)通过压缩存储技术和并行技术的集成实现了可以处理大数据特征的反演算法;(2)利用对数反演来约束每次计算得到新模型的电导率恒为正,充电率在0到1之间变化,从而保证反演的稳定性和可靠性.本文首先设计了两组对比模型进行试验,通过对不同区块数据采用不同加权的方式来减弱噪声对反演结果影响的效果;其次,采用并行技术提高了反演的计算速度,并利用理论模型分析了不同电极装置对反演分辨率的影响.最后,在甘肃某金属矿开发前景区利用激发极化法开展了中梯装置的采集工作,利用加权后的观测数据反演推断出了测区金属矿开发靶区的大致位置及分布特征.
关键词: 三维激发极化法      对数反演      并行技术      反演分辨率      金属矿勘探     
Study on three-dimensional inversion of induced polarization data and applications in mineral exploration
ZHANG Gang1,2,3, LÜ QingTian2,3     
1. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China;
2. SinoProbe Center, Chinese Academy of Geological Sciences, Beijing 100037, China;
3. Institute of Geophysical and Geochemical Exploration, Chinese Academy of Geological Sciences, Hebei Langfang 065000, China
Abstract: The induced polarization tomography has been extensively used in mineral exploration,investigation of sulfide ore,and other resources detection. With the development of exploration instruments and computer hardware,the amount of observed data in field surveys is getting larger and electrode array is more diverse. Based on the characteristics of observed data,we develop an inversion strategy for the induced polarization data that integrates the parallel technology and log-barrier method. The inversion strategy has the following characteristics:(1) The inversion strategy can handle large amounts of observed data with the compression storage and parallel technology; (2) Using the log-conductivity and log-barrier methods constrain the inverse model more reasonable that the resistivity is guaranteed to be positive and the chargeability always varies between 0 and 1. We firstly conduct two groups of comparative experiments with synthetic models. The experiments aim at controlling the effects of the noise in inversion result with giving the different weighting coefficients to the datum acquired in different areas that simulate the actual situation in the field. Second,the parallel technology is adopted to improve the inversion calculation speed. Meanwhile,we study the effects of the type of electrode array on the inversion resolution with synthetic models. Last,we employed the induced polarization tomography with the central gradient array in a potential area of mineral exploration in Gansu province,the inverted maps from the inversion with different weighting coefficients shows that the approximate location and subsurface three-dimensional resistivity and chargeability distribution of the mineral developing target in the survey area are delineated.
Keywords: Three-dimensional induced polarization inversion    Log-barrier inversion    Parallel technology    Inversion resolution    Mineral exploration    
0 引言

激发极化法可分为两种:一种是观测随时间变化在稳定电流激发下电场的激电效应,被称为时间域激发极化法.另一种是观测在交变电流作用下的电场随频率变化的激电效应,被称为频率域激发极化法.通过近半个世纪以来的实践证明,激发极化法是勘查各类金属矿产的重要方法之一,特别是在勘查电阻率与围岩相差不大的浸染状金属矿床.此外,在寻找地下水、天然气和探测石油等方面也均有一定的成效,取得了较好结果.焦方谦等(2013)在西伯利亚与哈萨克斯坦—准噶尔这两大板块之接合部位的陆缘裂谷带附近采用激发极化法圈定激电异常区,并成功区分矿与非矿异常.周洪祥等(2013)在某地区利用直流激发极化法获取测区的低电阻率异常带同时伴随高极化率异常反映,结合物性测定结果以及地面地质情况,发现了铅锌矿开发远景区.林方丽等(2016)在长江中下游多金属成矿带与华南成矿带之间的江南造山带上的乌溪矿区利用激电方法同时结合多种物探方法构建了研究区的成矿动力学模型,推测了研究区成矿机制,揭示了矿区的成矿潜力.周多等(2016)采用激电中梯与激电测深两种方法进行金属矿产勘查,圈定多个激电异常,并发现了矽卡岩型铜多金属矿开发的远景区.董和平和王勇(2017)利用激电中梯法和激电对称四极法来研究位于亚洲最大的石墨成矿带上的靶区深部电性特征,成功圈定成矿有利部位,并推断地层层位和构造特征.潘燕兵等(2017)利用激发极化法成功圈定撒哈拉沙漠金矿床的两条与黄铁矿化关系密切的金矿体.孙仁斌等(2017)通过激电中梯扫面和测深工作,发现了高极化率异常与锌钨多金属矿体具有良好的空间对应关系,大致圈定了与成矿有关的金属硫化物富集体范围.张晓东等(2017)利用激电测量所获取的视极化率异常解释了控矿构造,为实施钻探工程提供了强有力证据.赵后越(2017)在河北省石墨矿成矿远景重点勘探区利用激发极化法开展工作,联合地质资料进行解译,成功圈定两处有利异常,并推断了与异常相关的两处构造.

激电反演理论主要依据Seigel(1959)提出的视极化率概念,通过直流电阻率法正演假设没有激发极化效应条件下的电位差或者视电阻率,然后正演存在激发极化效应条件下的等效视电阻率,最后通过简单计算可以得到视极化率.这种方法的优势在于将直流电阻率法反演和激发极化反演联系到了一起,激发极化法反演方法首先是通过直流电阻率法反演计算得到地下电导率模型,将该电导率模型设为“真实”地下电导率模型进行第二步的激电数据的反演计算获取地下充电率模型.早期,线性反演方法以Pelton等(1978)Sasaki(1982)为代表,通过求解线性矩阵方程获取第一步直流法电阻率反演电阻率结果,进而获取地下充电率模型.Constable等(1987)提出最小构造方法减弱了多参数引起的超定问题带来的反演固有的非唯一性问题,最终得到最优解可以较好的拟合观测值.La Brecque(1991)将该方法应用在跨孔观测方式的成像试验中.Beard等(1996)提出了近似激电反演方法解决了目标体与围岩电阻率差异较小的模型反演.Ye等(2014)采用非结构化三角形网格进行剖分,且利用对偶加权后验误差估计指导网格自动细化过程,并依据Seigel(1959)理论实现激发极化法2.5维自适应有限元正演模拟算法.李兆祥等(2015)提出了基于交叉梯度约束的时间域激电数据二维同步反演策略,实现了交叉梯度约束的电阻率和极化率二维同步反演算法.

野外激电法的工作采集的数据具有观测数据量大、布极方式多样等特点.针对实测数据的特点,为了减少计算量以达到加快正演速度的目的,对稀疏矩阵进行压缩对角存储的处理,该方法可以方便地应用到各种迭代法中,从而使得用迭代法求解系数矩阵为带状阵的线性方程组时,内存空间得到节省(邹桂红等, 2013;吴小平等, 1998).随着高性能并行机的发展,有大量学者对并行计算进行了大量的研究(Iwashita and Shimasaki, 2002; Iwashita et al., 2005; Loke et al., 2010; Schwarzbach et al., 2005).

本文主要参考Li和Oldenburg(2000)反演算法原理并针对野外实测数据的特点,具体研究实现了能够处理大数据量的并行方法及能够约束每次反演迭代结果在合理范围的对数方法于一体的反演算法.并且通过模型试验研究针对不同区域数据使用不同的加权系数的方式来减弱噪声对反演结果的影响,最后将该方法应用于甘肃某金属矿开发远景区进行调查研究,取得了较好效果.

1 反演算法

稳定电流场条件下满足的偏微分方程,如(1)式所示:

(1)

其中,ϕσ是不考虑地下介质存在充电率前提下的电位(差);I是电流;δ是狄利克雷函数;rrs分别是测点与供电点的位置;σ是电导率.

Dey和Morrison(1979)对偏微分方程(1)进行区域积分后再进行离散化处理,即可以获取正演的离散的矩阵方程.通过有限差分方法离散后的矩阵方程具有大型、带状、稀疏、对称及正定特征,针对这些特征的方程通常采用不完全乔列斯基分解的共轭梯度方法(ICCG)求解使得迭代计算更加高效(Meijerink and Van Der Vorst, 1977; Poole and Ortega, 1987).

假设研究区区域Ω中存在一个点电源,地面处Γ0(z=0)的电位满足诺依曼边界条件,区域其他边界(Γ)处的电位满足混合边界条件,则电场满足的边值条件表示如下:

(2)

其中,n为边界外法向量,r为源点到边界的距离,θnr的夹角.混合边界条件不仅在边界处保持电位的连续性,也消除了对地下边界虚源所产生的反射作用.采用混合边界条件进行正演模拟得到的结果精度较高(Dey and Morrison, 1979).

那么根据公式(1)和(2)可以将地下无激发极化效应情况下的地电场电位公式简化为如式(3)所示,其中为正演算子.

(3)

当地下介质存在充电率为η时,那么此时的总电位ϕη与不存在充电率条件下的ϕσ计算公式有所不同.根据Seigel公式可以得到充电率条件下的电导率为σ*=σ(1-η),这样我们可以通过两次正演计算出视极化率ηa.

(4)

(5)

激发极化法反演的目标函数为

(6)

其中,μ为正则化因子,控制反演过程中数据拟合和模型拟合的权重;ΨdΨm分别为数据约束和模型约束,u为模型充电率的最大值.

结合Oldenburg和Li(1994)提出的非线性反演方法来实现三维激发极化法数据的反演与解释,其反演的雅克比矩阵如式(7).

(7)

其中,σ=σb, j(1-η).通常采用Rodi方法来计算雅克比矩阵G与任意向量x的乘积Gx以及雅克比矩阵的转置GT与任意向量y的乘积GTy,因此加快了反演计算速度.

三维激发极化法反演算法流程包括首先进行视电阻率(电位)数据反演得到地下电阻率模型,为使地下电阻率模型每次迭代计算结果恒为正值(电导率为非负值),我们拟采用log-conductivity方法(Oldenburg and Li, 1994;吴小平和徐果明, 2000)对地下模型电导率进行约束.

,那么根据(8)式可知,每次反演迭代后计算结果恒为正.

(8)

(8) 式可以表示为

(9)

那么雅克比矩阵变为

(10)

(11)

在完成第一步的电阻率反演后,利用反演所得三维电阻率模型作为初始模型进行视极化率数据的反演,为了约束地下充电率范围在0~1之间,我们采用log-barrier方法(Li and Oldenburg, 2000; Sastry et al., 2011; Sastry et al., 2014; Wright, 2001)进行视极化率数据的约束反演.具体实现步骤如下:

(1) 设定初始充电率模型m和正则化因子μ,利用公式(12)计算参数λ的初始值

(12)

其中,Ψd=‖Wd(dpre-dobs)‖2Ψm=‖Wm(m-m0)‖2分别为反演目标函数中的数据约束和模型约束,dpre为利用反演获取的新模型正演计算得到的预测数据,dobs为实测数据.μ的取值采用广义交叉验证(GCV)方法(Golub et al., 1979).

(2) 反演迭代计算新的模型矩阵Δm

(13)

其中,X=diag{m1, …, mM},Y=uI-Xe为单位矩阵.

求解模型修正量的流程如下:

(ⅰ)第k次正演得到数据拟合差矩阵dk-dobs,其中k=1, 2, …, Nmax

(ⅱ)计算公式(13)的右端项:

(ⅲ)共轭梯度(CG)迭代初始化各参数:

(ⅳ)CG迭代计算得到新的模型修正量Δmk

(ⅴ)CG迭代终止,存储模型修正量Δmk和计算并存储第k+1次反演新模型mk+1mk+1=mkmk.

(3) 确定最大模型修正量

(14)

(4) 更新模型和λ参数

(15)

其中,Gill等(1991)γ取值的范围是(0.99,0.999);Li和Oldenburg(2000)通过大量反演得到较好的经验γ取值0.925,也是本文中采用的取值.

(5) 返回第2步并迭代计算直至反演收敛至可接受精度内.

通过log-conductivity及log-barrier方法将每次反演迭代计算所获取的地下电阻率模型和充电率模型参数值控制在合理范围内.

针对任意四极装置观测条件下的直流电阻率法雅克比矩阵G的推导如下:令三维网格节点数为L,网格单元个数为n.

对于正问题有如下关系式:

(16)

其中,A为耦合系数矩阵;φi, j, k为空间网格中第(i, j, k)个节点上电位;S为由剖分网格体的几何参数、电导率及相邻七点的正常电位值相关的计算值.

供电极AB在观测点MN处激发的电位差为

(17)

其中,aman为M、N点的空间位置编号.

偶极源条件下M点的偏导数矩阵为

(18)

可以推导出任意四极装置条件下的偏导数矩阵为

(19)

于是得到

(20)

对于激电数据反演过程中雅可比矩阵的处理与电阻率反演类似,采用Rodi提出的方法来计算雅克比矩阵G与任意向量x的乘积Gx以及雅克比矩阵的转置GT与任意向量y的乘积GTy,从而达到加快了反演计算速度的目的.

随着观测设备及计算机硬件的发展,我们可以采集到更多的数据和拥有更优异硬件的计算机,同时对反演算法的计算能力和速度的要求也在不断提升.基于此,并行技术越来越广泛地应用于地球物理数值模拟和反演计算中来.并行计算是一种一次可执行多个指令的算法,目的是提高计算速度,及通过扩大问题求解规模,解决大型而复杂的计算问题.所谓并行计算可分为时间上的并行和空间上的并行.时间上的并行就是指流水线技术,而空间上的并行则是指用多个处理器并发的执行计算.对于多源的激发极化法反演而言,对正反演模块进行并行化处理是较好的选择.基于Visual Studio.NET开发平台,利用静态类Parallel中的三种方法For、ForEach、Invoke对正演模块和反演模块中相关算法部分进行优化实现并行集成.激发极化反演层析成像符合并行化计算的要求,激发源在不同位置发射信号时的正演计算是相互独立的事件,可以将不同源正演计算交予计算机不同线程完成;同时在反演计算模块中,对反演计算速度影响较大的是偏导数矩阵(Jacobi)的计算,可以实现不同源的Jacobi矩阵交予不同的计算机线程完成.并行化处理的详细流程如图 1所示.

图 1 正反演模块并行化流程图 Fig. 1 Numerical simulation and inversion calculation with parallel processing
2 理论模型试验 2.1 权重因子对噪声压制效果分析

反演的非唯一性由众多因素引起,噪声造成的观测误差就是其中之一.通过对观测数据加权的方式使得不同区域数据参与反演计算的权重不同,从而达到压制噪声对反演结果的影响.设定观测数据的权重因子来改变信噪比高和信噪比低的数据在反演过程中的参与度,从而达到对噪声的压制的效果.设计低阻高极化阶梯状具有同源性目标体的电阻率模型和充电率模型如图 2a所示,目标体电阻率为40 Ωm,充电率为15%;背景围岩的电阻率为100 Ωm,充电率为1%.供电方式采用中梯装置如图 2b所示.

图 2 理论模型(a)和中梯装置(b) 黑色圆点为供电极(A、B),黑色叉为测量极(M、N). Fig. 2 Synthetic model (a) with low-resistivity (40 Ωm) and the central gradient array (b) The black dots and crosses stand for the current (A, B) and potential electrodes (M, N), respectively.

对不加入高斯白噪声的数据但对不同区域数据加权后进行反演计算得到结果如图 3所示.数据加权方式为对目标异常体在地表投影所覆盖观测点数据加权1.0或0.2,与之对应将非覆盖区数据加权0.2或1.0.反演结果表明,当异常目标体在地表投影区域观测数据权重大于其他区域时,反演结果(图 3e3f)最优;而目标体投影区域数据参与计算的权重较低时,电阻率反演结果(图 3g)可以圈定出两到三层深度低阻块体的位置,而极化率反演结果(图 3h)则只能识别低阻目标体的浅层边界位置.结果表明,离目标异常体近的观测数据对目标体更敏感,在每次反演迭代计算中对模型修正量的贡献最大;反之,远离目标体的数据对每次迭代计算中模型修正量的贡献较小.所以当目标投影区覆盖数据权重大于其他区域数据时,反演效果最好.

图 3 理论模型主剖面(y=5 m)电阻率图(a)和充电率图(b). (c)—(h)为未加入高斯白噪声时对观测数据取不同权重因子时反演结果 其中(c)和(d)对所有观测数据权重均取1.0,(e)和(f)对目标异常体在地表投影覆盖区数据取权重为1.0而其他区域数据取0.2,(g)和(h)则目标异常体在地表投影覆盖区数据取权重为0.2而其他区域数据取1.0. Fig. 3 The main section at y=5 m of synthetic resistivity (a) and chargeability (b) model. (c) to (h) are the inverted models using the pure data without adding the Gaussian white noise (c) and (d) are the inverted models using the observed data with a diagonal weighting matrix of 1.0. (e) and (f) are the inverted models using the observed data with weighting 1.0 in the target projection area and weighting 0.2 in the other areas. (g) and (h) are the inverted models using the observed data with weighting 0.2 in the target projection area and weighting 1.0 in the other areas.

其次,依据野外实际施工情况对模拟观测数据加入不同比例的噪声.即供电极距较远时,观测数据的信噪比相对较低;反之,观测数据信噪比较高.据此原理,对观测数据加入高斯白噪声比例见表 1.然后利用加权最优的方式进行反演计算,分析与讨论数据加权反演对噪声压制效果.

表 1 不同供电极距所观测数据加入不同比例高斯白噪声 Table 1 Adding different proportions of Gaussian white noise in the observed data set based on the different polar distance

图 4是选择最优加权方式和不加权条件下反演结果比较图.不加权时反演电阻率结果(图 4c)在边界处出现多余的低阻假构造;充电率结果(图 4d)显示目标体底界面发散较大难以约束其内边界.加权后反演电阻率结果(图 4e)表明多余构造依然存在,但是相较于图 4c有所减少,说明加权后数据反演压制了噪声对反演结果的影响;充电率结果(图 4f)对目标体底界面约束有所改善.

图 4 理论模型主剖面电阻率图(a)和充电率图(b). (c)—(f)为对观测数据加入高斯白噪声后,分别采用和不采用最优权重因子时的反演结果 其中(c)和(d)为不采用权重因子直接反演结果,(e)和(f)为通过加入权重因子方式减小信噪比低数据的权重实现约束反演的结果. Fig. 4 The main sections of synthetic resistivity (a) and chargeability (b) model. (c) to (f) are the inverted models from the inversions with or without adding the optimal weight coefficients for the observed data with Gaussian white noise (c) and (d) are the inverted models from the directly inversion without adding the weight coefficients. (e) and (f) are the inverted models from the constrained inversion with adding the weight coefficients.
2.2 并行模块对计算效率优化分析

为了研究并行化正反演模块的计算效率,分别利用串行化和并行化的正反演模块进行计算并记录反演结束时所用时间.然后通过控制计算机线程数来执行反演计算可以得到不同线程条件下的运行时.计算机的配置为Z820 WorkStation Interl (R) Xeon (R) CPU E5-2620 v2 @2.1 GHz RAM16GB 64位操作系统.

反演模型设定:数据空间观测数据量为815400,模型空间网格块体数为60000(即包含:60000个块体电导率参数和60000个充电率参数).反演参数的设定为:CG迭代10次,INV迭代20次.分别利用并行算法和串行算法来执行反演计算,统计每次计算运行时间并分析并行算法模块对反演计算效率的提升.

通过对CPU的24线程采用逐渐关闭的方式统计反演运行时间和CPU占用率得到如表 2的结果和图 5.根据表 2图 5可以得到:(1)随着线程的增加并行和串行计算时间都有所减少,但是并行算法对反演计算效率的提升更显著(采用24线程时,并行方法大约是串行方法的3.07倍计算效率),从图 5中可以更直观地看出计算效率的提升情况;(2)并行算法中CPU都能被充分利用,而串行算法则只是在单线程中计算从而使得反演计算效率低.

表 2 基于并行与串行算法的激电反演的计算时间和CPU占用率对比 Table 2 Running time and CPU utilization of the induced polarization data inversion with the parallel and sequential algorithm
图 5 并行与串行算法实现反演计算所需时间及CPU占用率对比 Fig. 5 Running time and CPU utilization for the inversion calculation with the parallel and sequential algorithm
2.3 电极装置对分辨率的影响分析

为了研究电极装置对反演结果的影响,我们建立如图 6a所示的模型.在1000 m×1000 m×500 m的模型空间中,建立z轴垂直向下的笛卡尔坐标系.

图 6 横向与纵向上分布不同物性目标体模型 Fig. 6 The synthetic model with targets with different physical properties layout in horizontal and vertical lateral distribution

理论模型的反演试验主要目的是研究电极装置对反演分辨率的影响.为了研究横向及纵向分辨率,理论模型设计如下:在地表(z=0 m)放置宽度分别为100 m、200 m、300 m,厚度均为100 m的不同电阻率分布的三个异常体,同时在深度为200 m放置两个水平方向等尺度大小且厚度均为100 m的高阻和低阻目标体,在深度为400 m处相同位置放置大小一致的低阻和高阻异常体(图 6).目标体的物性参数如图 6所示,背景围岩为100 Ωm,充电率为1%,异常体分别为低阻高极化(ρ=10 Ωm, η=30%)与高阻高极化体(ρ=200 Ωm, η=30%).设计电极装置包括两极装置(图 7a)、三极装置(图 7b)及四极装置(图 7c).

图 7 电极装置图 (a)两极装置(A-M);(b)三极装置(A-MN);(c)四极装置(A-MN-B). Fig. 7 Electrode arrays (a) Pole-pole array; (b) Pole-dipole array; (c) Dipole-dipole array.

三维反演结果如图 8所示,其中(a)、(b)、(c)为不同电极装置数据反演电阻率结果在深度z=50 m、250 m、450 m时y=500 m的剖面曲线;(d)、(e)、(f)为不同电极装置数据反演充电率结果在深度z=50 m、250 m、450 m时y=500 m的剖面曲线.反演结果表明:(1)纵向对比反演不同深度上电阻率结果和极化率结果可以看出,反演浅层(z=50 m)异常体效果明显优于深部(z=450 m)异常体,说明三种装置采集数据均对浅层异常体较敏感而对深部异常体不如浅部敏感;(2)无论电阻率反演结果还是极化率反演结果,横向分辨率比纵向分辨率高;(3)深度为450 m时,三极装置和四极装置所采集数据反演效果优于两极装置,对位于x=350~450 m之间的高阻体边界能够识别出来,而三种装置数据反演结果对边部区域异常体边界的刻画不如对研究区中心位置异常体;(4)图 8e能够较明显地看出极化率反演结果对低阻高极化体圈定较好而对高阻高极化体的识别能力较差.

图 8 不同电极装置数据反演的电阻率和充电率结果 (a)、(b)、(c)为不同电极装置数据反演电阻率结果在深度z=50 m、250 m、450 m时y=500 m的剖面曲线;(d)、(e)、(f)为不同电极装置数据反演充电率结果在深度z=50 m、250 m、450 m时y=500 m的剖面曲线. Fig. 8 Inverted resistivity and chargeability results with different types of electrode arrays (a), (b), and (c) are inverted resistivity profile of y=500 m at depths 50 m, 250 m and 450 m, respectively; (d), (e), and (f) are inverted chargeability profile of y=500 m at depths 50 m, 250 m and 450 m, respectively.
3 应用实例 3.1 地质背景

研究区甘肃某银铅锌多金属矿田位于北山造山带中西部(图 9),大地构造位置属塔里木古陆东北部敦煌地块北缘双鹰山早古生代裂谷型被动陆缘带内的中元古生代裂谷中.北部为一早古生代被动边缘裂谷裂陷带,其中夹杂有微陆块和中元古带裂谷;南部为一古生代多期裂谷带,二者之间以红柳园微陆块相隔.敦煌陆块主要由太古宇敦煌群深变质岩构成,西以阿尔金断裂为界与塔里木陆块为邻,东南部以黑河断裂与阿拉善陆块相隔.早期学者认为测区银铅锌矿床类型为矽卡岩-高中温热液型;后来推断为火山-沉积热液改造型和岩浆热液型,后又有矽卡岩型矿床的观点.2000年之后,大多研究者认为多金属矿田为喷流-沉积矿床类型(杨建国等, 2010;代文军, 2010).

图 9 (a) 测区位置及成矿带构造单元划分示意略图;(b)激发极化法电极分布图;(c)测量仪器;(d)发射机 Fig. 9 Location of the survey area and the geological setting of the core region (a). (b) is the electrode arrangement. (c) and (d) are the measuring equipment and the transmitter, respectively
3.2 激电数据解释

地质背景资料提供了测区部分铅锌多金属矿区地质环境、矿床成因类型等信息.圈定和约束了金属矿开发前景区的调查范围,且喷流-沉积矿床类型多呈低阻高极化异常体,也通过钻孔岩性分析资料论证了该信息.利用激电方法对金属矿开发前景区进行高分辨率的调查工作可以对测区金属矿分布进行更细致的解释工作是较好的选择.

激电仪器(EDJD-1)如图 9c9d所示,该仪器的电压测量范围为±30 V,电压测量精度为±0.1%,电压最高采样分辨率为0.1 μV,视极化率测量精度为±0.2%,50 Hz工频干扰抑制为80 dB.激电观测方式如图 9b所示,观测电极装置采用四极装置(中梯装置).电极距为50 m,测线之间距离为100 m.根据观测数据空间可以设计反演模型空间尺度为800 m×1200 m,三维模型空间网格块体剖分最小单元尺度为50 m.

供电极AB距离分别为1000 m、1500 m、2000 m、2500 m和3000 m时所观测的原始等效视电阻率和视极化率等值线如图 10所示.对比视电阻率和视极化率结果图,圈定出的相对低阻异常和相对高极化异常带共有2处,其中南部区域异常带规模较大,而北部区域的异常带规模相对较小.低阻区和高极化区对应较好的重叠区在成果图中利用虚线圈定出来.

图 10 视电阻率(a—e)和视极化率(f—j)等值线图 (a—e)和(f—j)分别为供电极AB=1000 m、1500 m、2000 m、2500 m、3000 m时视电阻率和视极化率图.其中虚线为推断出的金属矿开发前景区. Fig. 10 Apparent resistivity (a—e) and chargeability (f—j) maps From (a) to (e) and (f) to (j) stand for the apparent resistivity and polarization maps with AB=1000 m, 1500 m, 2000 m, 2500 m, 3000 m, respectively. The area within the dashed black line box infers to the potential area for mineral exploration.

根据图 10所示的视电阻率和视极化率等值线图,可以推断金属矿位于虚线所圈定的两片低阻高极化区域.我们对低阻高极化区域所覆盖的观测数据和其他区域数据分别进行加权1.0和0.2后进行反演计算得到如图 11所示的三维电阻率和充电率分布结果.反演结果的浅部与视电阻率和视极化率等值线结果图具有较好一致性,虚线圈定区域吻合较好.尤其是在北部的低阻区域与高极化区域对应较好,而南部的低阻区域与高极化空间分布对应不是很好,推断为受到近地表低阻覆盖层影响所致,需要更进一步的调查工作进行确认.

图 11 反演结果图.左:电阻率三维分布;右:充电率三维分布 虚线为利用视电阻率图和视极化率图推断的金属矿开发前景区,红色实线为反演结果对应的低阻高激化区域. Fig. 11 Inverted resistivity (left) and chargeability (right) models The area within the dashed line in both models illustrates the potential area for mineral exploration with apparent resistivity and polarization maps. The red solid line infers to be the potential area for mineral exploration with inverted maps.
4 结论

本文利用有限差分方法实现了能够保证反演稳定性的对数方法与能够提供反演速度的并行技术集于一体的三维激电反演算法.首先通过理论模型测试得到结论如下:

(1) 对不同区域观测数据加入不同权系数,使对目标体敏感区数据参与反演的权重大于远离目标体(或目标研究区)的数据,可以更好的压制噪声对反演结果的影响;

(2) 利用理论模型分别模拟采用并行方法与串行方法时的反演计算所消耗的时间和CPU利用率,对比两种方法计算结果可知,在Z820 WorkStation Interl (R) Xeon (R) CPU E5-2620 v2 @2.1 GHz RAM16GB 64位操作系统的计算机上并行计算可以提高反演计算效率3.07倍,且当计算机为多核时每次迭代计算都能够充分利用计算机的线程数;

(3) 利用模拟四极装置和三极装置观测所获取电阻率和极化率数据反演结果对异常体目标体圈定较好;

(4) 极化率反演结果对高阻高极化体的圈定不如低阻高极化体好.

最后将激发极化法反演方法应用到甘肃实测数据中,通过观测视电阻率和视极化率判断出目标体的大致位置.对目标区覆盖数据进行加权反演计算结果更加精细的刻画和判断了金属矿开发前景区的位置,弥补了直接通过视电阻率和视极化率等值线图能够圈定金属矿在水平地表的位置投影而难以判断在深度上的信息.

致谢  感谢中国地质科学院地球深部探测中心、中国地质科学院地球物理地球化学勘查研究所和中国地质大学(北京)研究团队成员的积极支持和帮助.
References
Beard L P, Hohmann G W, Tripp A C. 1996. Fast resistivity/IP inversion using a low-contrast approximation. Geophysics, 61(1): 169-179. DOI:10.1190/1.1443937
Constable S C, Parker R L, Constable C G. 1987. Occam′s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3): 289-300. DOI:10.1190/1.1442303
Dai W J. 2010. Metallogenic study of Huaniushan Au-Ag-Pb-Zn Deposit, Beishan area, Gansu province. Geology and Mineral Resources of South China (in Chinese), 2010(3): 25-33. DOI:10.3969/j.issn.1007-3701.2010.03.004
Dey A, Morrison H F. 1979. Resistivity modeling for arbitrarily shaped three-dimensional structures. Geophysics, 44(4): 753-780. DOI:10.1190/1.1440975
Dong H P, Wang Y. 2017. The research on graphite ore deposit by using induced polarization method in Xinhure, Wulatezhongqi. Chinese Journal of Engineering Geophysics (in Chinese), 14(4): 421-426. DOI:10.3969/j.issn.1672-7940.2017.04.007
Gill P, Murray W, Ponceleon D, et al. 1991. Solving reduced KKT systems in barrier methods for linear and quadratic programming. Technical Report SOL, 91(7): 1-26.
Golub G H, Heath M, Wahba G. 1979. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21(2): 215-223. DOI:10.2307/1268518
Iwashita T, Shimasaki M. 2002. Algebraic multicolor ordering for parallelized ICCG solver in finite-element analyses. IEEE Transactions on Magnetics, 38(2): 429-432. DOI:10.1109/20.996114
Iwashita T, Shimasaki M, Lu J W. 2005. Parallel ICCG solvers for a finite-element eddy-current analysis on heterogeneous parallel computation environment. // IEEE/ACES International Conference on Wireless Communications and Applied Computational Electromagnetics, 2005. Honolulu, HI, USA: IEEE.
Jiao F Q, Zhao X S, Chen C. 2013. The application of the metallic factor to the interpretation of IP anomalies. Geophysical and Geochemical Exploration (in Chinese), 37(5): 848-852. DOI:10.11720/j.issn.1000-8918.2013.5.17
La Brecque D J. 1991. IP tomography. //61st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 413-416.
Li Y G, Oldenburg D W. 2000. 3D inversion of induced polarization data. Geophysics, 65(6): 1931-1945. DOI:10.1190/1.1444877
Li Z X, Tan H D, Fu S S, et al. 2015. Two-dimensional synchronous inversion of TDIP with cross-gradient constraint. Chinese Journal of Geophysics (in Chinese), 58(12): 4718-4726. DOI:10.6038/cjg20151232
Lin F L, Wang G J, Yang X Y. 2016. Application of comprehensive electromagnetic study in deep mineralization mechanism—A case study of the Wuxi polymetallic ore deposit, south Anhui. Chinese Journal of Geophysics (in Chinese), 59(11): 4323-4337. DOI:10.6038/cjg20161132
Loke M H, Wilkinson P B, Chambers J E. 2010. Parallel computation of optimized arrays for 2-D electrical imaging surveys. Geophysical Journal International, 183(3): 1302-1315. DOI:10.1111/j.1365-246X.2010.04796.x
Meijerink J A, Van Der Vorst H A. 1977. An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. Mathematics of Computation, 31(137): 148.
Oldenburg D W, Li Y G. 1994. Inversion of induced polarization data. Geophysics, 59(9): 1327-1341. DOI:10.1190/1.1443692
Pan Y B, Chen Y Y, Wang J W. 2017. Application of IP intermediate gradient method in vein type gold mine of desert. World Nonferrous Metal (in Chinese), 4(17): 107-110.
Pelton W H, Rijo L, Swift C M. 1978. Inversion of two-dimensional resistivity and induced-polarization data. Geophysics, 43(4): 788-803. DOI:10.1190/1.1440854
Poole E L, Ortega J M. 1987. Multicolor ICCG methods for vector computers. SIAM Journal on Numerical Analysis, 24(6): 1394-1418. DOI:10.1137/0724090
Sasaki Y. 1982. Automatic interpretation of induced polarization data over two-dimensional structures. Memoirs of the Faculty of Engineering, Kyushu University, 42(1): 59-74.
Sastry S P, Shontz S M, Vavasis S A. 2011. A log-barrier method for mesh quality improvement. //Proceedings of the 20th International Meshing Roundtable. Berlin: Springer, 329-346.
Sastry S P, Shontz S M, Vavasis S A. 2014. A log-barrier method for mesh quality improvement and untangling. Engineering with Computers, 30(3): 315-329. DOI:10.1007/s00366-012-0294-6
Schwarzbach C, Brner R U, Spitzer K. 2005. Two-dimensional inversion of direct current resistivity data using a parallel, multi-objective genetic algorithm. Geophysical Journal International, 162(3): 685-695. DOI:10.1111/j.1365-246X.2005.02702.x
Seigel H O. 1959. Mathematical formulation and type curves for induced polarization. Geophysics, 24(3): 547-565. DOI:10.1190/1.1438625
Sun R B, Chu L X, Zhao Y J, et al. 2017. Application of the time domain induced polarization method to a Zn-W polymetallic deposit in Xianghuang Banner, Inner Mongolia. Geology and Exploration (in Chinese), 53(3): 519-527.
Wright S J. 2001. On the convergence of the Newton/log-barrier method. Mathematical Programming, 90(1): 71-100. DOI:10.1007/PL00011421
Wu X P, Xu G M, Li S C. 1998. The calculation of three-dimensional geoelectric field of point source by incomplete cholesky conjugate gradient method. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 41(6): 848-855.
Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method. Chinese Journal of Geophysics (in Chinese), 43(3): 420-427.
Yang J G, Zhai J Y, Yang H W, et al. 2010. Metallotectonics and prospection of the huaniushan exhalogene Gold-Silver-Lead-Zinc deposit in Beishan, Gansu Province. Geotectonica et Metallogenia (in Chinese), 34(2): 246-254. DOI:10.3969/j.issn.1001-1552.2010.02.011
Ye Y X, Li Y G, Deng J Z, et al. 2014. 2.5D induced polarization forward modeling using the adaptive finite-element method. Applied Geophysics, 11(4): 500-507. DOI:10.1007/s11770-014-0455-z
Zhang X D, Fang J, Zhang D Y, et al. 2017. The application of the induced polarization method to the replacement resources exploration of the Dongxi gold deposit. Geophysical and Geochemical Exploration (in Chinese), 41(3): 445-451. DOI:10.11720/wtyht.2017.3.08
Zhao H Y. 2017. Application of high-power induced polarization to exploration of graphite ore in Zhangjiakou. Chinese Journal of Engineering Geophysics (in Chinese), 14(5): 546-551. DOI:10.3969/j.issn.1672-7940.2017.05.007
Zhou D, Chen A X, Dong Z M, et al. 2016. Application of the induced polarization method in exploration of a skarnized copper polymetallic deposit in Benxi, Liaoning province. Geology and Exploration (in Chinese), 52(4): 688-694.
Zhou H X, Peng Q Y, Jiang Y L. 2013. Applying the transient electromagnetic methods and the induced polarization method to searching for lead-zinc deposits. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 35(6): 664-667. DOI:10.3969/j.issn.1001-1749.2013.06.06
Zou G H, Liang H Q, Yin H D. 2013. Application frame using compressed diagonal storage format in iterative methods. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 35(1): 107-111. DOI:10.3969/j.issn.1001-1749.2013.01.17
代文军. 2010. 甘肃北山花牛山金银铅锌矿床成因探讨. 华南地质与矿产, 2010(3): 25-33. DOI:10.3969/j.issn.1007-3701.2010.03.004
董和平, 王勇. 2017. 乌拉特中旗新忽热石墨矿的激发极化法研究. 工程地球物理学报, 14(4): 421-426. DOI:10.3969/j.issn.1672-7940.2017.04.007
焦方谦, 赵新生, 陈川. 2013. 激电异常解释中金属因素的应用. 物探与化探, 37(5): 848-852. DOI:10.11720/j.issn.1000-8918.2013.5.17
李兆祥, 谭捍东, 付少帅, 等. 2015. 基于交叉梯度约束的时间域激发极化法二维同步反演. 地球物理学报, 58(12): 4718-4726. DOI:10.6038/cjg20151232
林方丽, 王光杰, 杨晓勇. 2016. 综合电磁法在矿区深部成矿机制中的应用研究——以皖南乌溪多金属矿区为例. 地球物理学报, 59(11): 4323-4337. DOI:10.6038/cjg20161132
潘燕兵, 陈媛媛, 王君伟. 2017. 激电中梯测量在沙漠地区脉型金矿勘查中的应用. 世界有色金属, (17): 107-110.
孙仁斌, 楚丽霞, 赵绎钧, 等. 2017. 时间域激发极化法在内蒙古镶黄旗某锌钨多金属矿的应用. 地质与勘探, 53(3): 519-527.
吴小平, 徐果明, 李时灿. 1998. 利用不完全Cholesky共轭梯度法求解点源三维地电场. 地球物理学报, 41(6): 848-855. DOI:10.3321/j.issn:0001-5733.1998.06.016
吴小平, 徐果明. 2000. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报, 43(3): 420-427. DOI:10.3321/j.issn:0001-5733.2000.03.016
杨建国, 翟金元, 杨宏武, 等. 2010. 甘肃花牛山喷流沉积型金银铅锌矿床控矿因素与找矿前景分析. 大地构造与成矿学, 34(2): 246-254. DOI:10.3969/j.issn.1001-1552.2010.02.011
张晓东, 方捷, 张定源, 等. 2017. 激发极化法在东溪金矿接替资源勘查中的应用. 物探与化探, 41(3): 445-451. DOI:10.11720/wtyht.2017.3.08
赵后越. 2017. 大功率激电在张家口地区石墨矿勘查中的应用. 工程地球物理学报, 14(5): 546-551. DOI:10.3969/j.issn.1672-7940.2017.05.007
周多, 陈安霞, 董再民, 等. 2016. 激发极化法在辽宁本溪某矽卡岩型铜多金属矿勘查中的应用. 地质与勘探, 52(4): 688-694.
周洪祥, 彭乾云, 江玉乐. 2013. 瞬变电磁法与直流激发极化法在寻找铅锌矿中的应用. 物探化探计算技术, 35(6): 664-667. DOI:10.3969/j.issn.1001-1749.2013.06.06
邹桂红, 梁华庆, 尹洪东. 2013. 迭代法中压缩对角存储的应用框架. 物探化探计算技术, 35(1): 107-111. DOI:10.3969/j.issn.1001-1749.2013.01.17