三分量定源瞬变电磁解释技术及其在金属矿区的实验
智庆全1, 武军杰1,2, 王兴春1, 陈晓东1, 杨毅1, 张杰1, 邓晓红1
1.中国地质科学院 地球物理地球化学勘查研究所,河北 廊坊 065000
2.长安大学 地质工程与测绘学院,陕西 西安 710054

作者简介: 智庆全(1987-),男,助理工程师,硕士,主要从事电磁法研究工作。E-mail:zhiqingquan@igge.cn

摘要

目前的瞬变电磁仪器已经为三分量瞬变电磁解释提供了硬件基础,但瞬变电磁的实用解释技术仍局限于对垂直分量的解释或一些简单近似技术。为充分利用三分量数据,进一步提高定源瞬变电磁方法的分辨率,针对定源瞬变电磁装置给出三分量全域视电阻率定义和三分量联合一维反演的方法。在进行全域视电阻率定义时,利用迭代求取反函数值的方法替代直接求取反函数的过程;引入三分量瞬变电磁数据的权重以构建联合三分量数据的目标函数,利用可行域法实现三分量联合一维瞬变电磁反演;最后将三分量解释算法分别应用在模型计算和金属矿的实际勘查资料上,以检验其有效性。模型试验和实际资料处理结果表明,三分量解释方法能够利用更多的瞬变电磁数据,获得更高分辨率的处理结果,且具有较快的计算速度。

关键词: 定源瞬变电磁; 全域视电阻率; 三分量联合反演; 金属矿
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2016)04-0798-06 doi: 10.11720/wtyht.2016.4.25
Three-component interpretation technique of fixed source TEM and its experimental application in metallic ore district
ZHI Qing-Quan1, WU Jun-Jie1,2, WANG Xing-Chun1, CHEN Xiao-Dong1, YANG Yi1, ZHANG Jie1, DENG Xiao-Hong1
1. Institute of Geophysical and Geochemical Exploration, Langfang 065000, China
2. School of Geology Engineering and Geomatics, Chang'an university, Xi'an 710054, China
Abstract

The TEM instruments have provided a hardware foundation for three-component TEM interpretation; nevertheless, the practical interpretation technology is still confined to the vertical component or simple approximation techniques. To make full use of the three-component data and improve the resolution of the fixed source TEM, this paper put forwards the three component full domain apparent resistivity definition and one-dimensional three-component joint inversion method. For the full domain apparent resistivity definition, the authors use the iterative method to calculate the inverse function value instead of directly calculating the inverse function. The authors introduce the three-component weights of TEM data to build the joint objective function, and use the feasible region method to implement three-component joint inversion. Finally, the three-component interpretation algorithm is applied to the model calculation and the actual exploration data of the metallic mine to verify its effectiveness. Model test and actual data processing results show that the three-component interpretation method can use more TEM data and obtain processing results with higher resolution, thus achieving a fast calculation speed.

Keyword: fixed TEM; full domain apparent resistivity; three-component joint inversion; metallic mine

定源瞬变电磁(或称大回线瞬变电磁)具有探测深度大、可以在框内外进行多分量采集、工作效率高的特点, 尤其适合深部的金属矿产资源勘查, 是一种应用广泛的瞬变电磁探测装置。目前已有多种仪器可实现三分量的瞬变电磁数据采集, 很多工程项目也采集了多分量数据, 但是在解释技术方面, 仍主要采用基于中心点垂直分量的电阻率解释、一维反演等手段。由于定源回线和中心回线装置的响应规律不同, 这种解释方法不能利用水平分量数据, 同时也是比较近似、粗糙的。

考虑到这种近似处理技术所带来的问题, 有必要针对定源瞬变电磁装置, 研究非中心点、多分量的正演和解释方法。目前, 在正演方面已发表的成果主要有:国外学者Poddar[1]、Raiche[2]等采用电偶极子叠加和线积分的方法, 给出了层状大地上多边形回线源的瞬变电磁场的算法; 国内学者薛国强[3]、翁爱华[4]等对利用电偶极子和磁偶极子叠加方法计算大回线源瞬变电磁响应进行了研究, 并进行了精度和特征分析。在解释方法的研究上, 席振铢等[5]通过对场的响应规律进行分析, 得出了水平分量对异常体位置较敏感的结论, 并通过水平分量与垂直分量的比值函数进行了三分量数据解释; 武军杰 67、戚志鹏 89、翁爱华 1011等研究了定源瞬变电磁的三分量视电阻率定义方法, 通过理论计算证明了瞬变电磁三分量解释方法能综合利用更多的信息, 获得更高的解释分辨率; 戚志鹏、李貅等[9]研究了综合利用三分量瞬变电磁数据进行一维反演的方法, 结果表明三分量联合反演能快速、准确地获得反演模型。这些研究证明了三分量联合解释方法能充分利用各分量数据, 进而提高定源瞬变电磁解释精度的优势, 为定源装置的三分量解释方法的实用化打下了基础。

前人的研究结果表明, 定源回线的水平分量在确定异常体位置等方面具有一些重要的优点。综合利用定源回线的垂直分量和水平分量进行三分量瞬变电磁解释, 能充分利用更多的瞬变电磁场信息, 提高解释效果。鉴于此, 笔者结合理论模型进一步研究定源瞬变电磁的三分量视电阻率定义和三分量联合反演方法, 并通过实际资料处理分析三分量解释方法在金属矿区的应用效果。

1 定源瞬变电磁全域三分量视电阻率定义方法

一般定义电磁法视电阻率的方式为:根据均匀半空间模型上电磁响应随电阻率的变化函数来确定其反函数, 利用反函数来计算视电阻率[12]。然而在定源回线装置中, 瞬变电磁响应随着电阻率的变化规律是十分复杂的, 其反函数求取几乎不可能实现, 仅仅能表示为一种隐函数形式。在以前的视电阻率定义研究中, 往往通过早、晚期或者近、远区近似将响应函数近似为简单函数, 分别定义早、晚期或近、远区视电阻率, 以在较小计算量的前提下获得近似的视电阻率值。由于这种近似条件受到地电断面特征的控制, 对于实测数据而言, 何种情况下采用早期或晚期(近区或远区)近似是很难确定的, 在不满足近似条件的情况下使用分区或分期视电阻率定义公式, 很可能会导致视电阻率剖面畸形, 造成解释工作上的误判。因此, 文中采用一种无需分期或分区的全域(指全期和全区)视电阻率定义方法。

根据电磁场理论, 由发射源激发的瞬变电磁响应可以统一地表达为

B=f(x, y, z, t, ρ, a, I0), (1)

其中, B为瞬变电磁响应(衰减电压或磁场值), xyz为测点坐标, t为采样时间点, ρ 为均匀半空间的电阻率, I0为阶跃电流的大小, a为回线几何参数, f代表观测场值关于装置参数(发射回线半径、发射电流、观测位置等)和模型参数的函数。将装置参数和观测参数看作函数中的参量, 当装置固定不变时, 略去这些参量不写, 即有

B=f(ρ)(2)

如果能求取上式的反函数, 即可获得瞬变电磁任意分量的视电阻率定义式。考虑到在大回线装置中, 上式的形式一般较为复杂, 反函数难以显式表达, 采取一种等价方法求取其反函数值。取

g(ρ)=B-f(ρ)B2, (3)

则求取式(2)的反函数值等价于求解方程

g(ρ)=0(4)

假定ρ 0g(ρ )零点ρ s的一个近似点, 将g(ρ )在ρ 0处作Taylor展开, 并忽略二阶以上高阶项, 进行移项整理, 即有

ρ=ρ0+g(ρ)-g(ρ0)g'(ρ0), (5)

式(5)即为迭代求取视电阻率值的基本格式。数学上已经证明, 此迭代格式是线性收敛的。

这种迭代的视电阻率定义方法考虑了时间、测点空间坐标等参数, 无需进行分期或分区处理, 属于全域的视电阻率定义方法。通过数值迭代, 可以求得瞬变电磁响应任意分量或其组合的视电阻率值。对于迭代过程中可能出现的“ 双解” 、“ 无解” 和“ 偏解” 现象[13], 根据迭代残差和迭代电阻率值予以剔除, 并通过插值方法进行“ 补全” , 以保证视电阻率曲线或剖面的完整性。

2 定源瞬变电磁三分量一维反演

通过反演技术可以根据观测数据来获得测量区域的电性分布情况, 从而提供较为准确的电阻率— 深度剖面辅助地质解释。对于地层分布具有近似层状规律的测区, 一维反演方法能快速、准确地给出反演剖面, 是目前较为实用的反演方法。文中将通过引入三分量权重的方法, 综合利用定源回线装置的三个分量进行一维反演, 以期利用更多的测量信息获得更为准确的反演结果。

瞬变电磁响应和地电模型之间的关系是非线性的, 为获得较实用的反演速度, 一般使用线性化技术进行瞬变电磁数据的反演。线性化反演的一个重要环节是目标函数设置, 合理的目标函数可以保证反演过程的收敛速度, 实现对真实模型的最佳拟合。

从最佳拟合观测数据的角度出发, 将目标函数设置为实测数据与拟合数据之差, 并利用方差控制各数据对于目标函数的权重, 即取目标函数为

Φ=rTCmr, (6)

其中r= dobs-dpredobs, Cm为方差矩阵, 数据向量d= BxtBytBzt

在进行定源装置测量时, 由于瞬变电磁场的特性随着位置坐标变化而不同, 水平分量数据和垂直分量数据往往在量级上存在差异, 甚至某些分量数据为零。这种差异对于目标函数值的影响是很大的, 进行反演时有必要加以考虑。这里采用引入三分量权重的办法, 即重新取数据向量

d=BxtBytBztC3, (7)

其中C3为三分量的权重。若某分量的幅值不足数据向量模值的0.1%, 或者场分量受到干扰影响较大时, 可将其相应权重设置为0, 该分量不参与反演计算。

不同分量和时间道上的瞬变电磁响应存在量级差异, 考虑到这种差异可能引起的奇异性, 在数据拟合时采用对数比例。此时目标函数应进行一定修改, 取

r=(lndobs-lndpre)C3, (8)

相应的目标函数修改为

Δ=eΦ/m-1(9)

在反演过程中, 使用修正后的目标函数。

瞬变电磁场随电阻率的变化近似满足双对数关系, 在线性坐标上, 较低电阻率的变化引起的电磁响应变化比较剧烈, 而较高电阻率的变化引起的电磁响应变化不大。为保证拟合过程的稳定性, 在进行拟合时同样将反演电阻率向量取为其对数。实际岩石的电阻率值一般分布在[10-1, 104] Ω · m范围内, 文中的反演均将电阻率值限制在此范围内, 以减小反演过程的多解性。对于目标函数的最小化, 采用可行域方法。

3 H型地电模型视电阻率定义和一维反演算例

首先设计H型地层进行视电阻率定义和反演试算, 以研究前文算法的正确性和可靠性。设计模型参数为:第一层电阻率值为100 Ω · m, 厚度为200 m; 第二层电阻率值为10 Ω · m, 厚度为100 m; 第三层电阻率值为100 Ω · m。发射参数设置为:方形发射回线边长为600 m, 发射电流为1 A, 关断时间为500 μ s。观测参数为:接收面积1 m2, 观测时间为 0.086~7.93 ms, 观测时间道按对数间隔分布, 共27道; 以发射回线框中心为原点, 在y=0的主剖面上布置观测点, 测点在x轴上分布于-800~800 m之间, 点距为100 m。利用数字滤波算法计算其瞬变电磁响应, 然后利用前文所述的视电阻率定义方法进行计算, 获得x分量、z分量的视电阻率剖面(图1、图2)。

图1 H型模型x分量视电阻率剖面

图2 H型模型z分量视电阻率剖面

从视电阻率断面图中可以看出, x分量和z分量视电阻率剖面大致上反映了模型电性在深度上的分布情况。x分量视电阻率剖面反映的低阻层在纵向上有较大的延伸, 这种特征在框内较框外更为明显, 使得低阻层在视电阻率剖面上呈现出了“ 外凸” 的现象。z分量视电阻率剖面对于低阻层在深度上的定位更为准确, 下延的特征较不明显, 但在框内时反映出的异常体厚度较大, 而在框外较小。在视电阻率剖面上, 部分点的电阻率值较大, 超出了100 Ω · m, 这是由全域视电阻率定义中的“ 偏解” 问题引起的。虽然对不合理电阻率值进行了剔除, 这种影响仍然造成了视电阻率剖面的畸变。因此, 相比于中心回线装置, 定源回线的视电阻率剖面对于真实电性分布的反映更为粗略, 仅仅只能作为解释工作的一种参考。

图3为联合使用x分量和z分量瞬变电磁数据进行一维反演的结果。初始模型取为电阻率值为100 Ω · m的均匀半空间。从图中可以看出, 通过联合反演所得到的深度和电阻率参数几乎完全和设计模型一致, 仅仅在电阻率界面部分稍有差异。数据的整体相对拟合误差为3.2%, 单点反演时间约为 1 min。这说明, 文中所给出的一维三分量联合反演方法是快速、有效的。

图3 H型模型x-z分量联合反演剖面

图4为仅使用z分量数据进行单分量反演的结果。初始模型仍取为电阻率为100 Ω · m的均匀半空间, 反演的整体拟合误差为2.0%。与联合使用x-z分量的反演结果相比, 仅使用z分量数据的反演剖面仍能比较准确地给出低阻层的位置参数, 但在电阻率交界处出现了电阻率稍高的低阻区域, 高、低阻界面相对模糊。二者的对比表明, 联合使用垂直分量和水平分量的反演方法能更好地实现异常体定位和确定电阻率参数, 具有更高的分辨率。

图4 H型模型z分量反演剖面

4 金属矿勘查的有效性实验

青海夏日哈木铜镍矿区位于青海省东昆仑西段祁漫塔格地区[14], 地处柴达木盆地西南缘, 是东昆仑成矿带新发现的一处超大型岩浆熔离型铜镍硫化物矿床。矿体多呈厚大似层状、透镜状, 少数漏斗状和不规则状矿体呈上悬式位于岩体中上部, 沿走向方向, 矿体中间厚、品位高, 向两侧趋于尖灭。在HS26异常区L9~L17线开展了三分量定源瞬变电磁法的野外试验, 测线与地质勘探线重合, 方向为近南北向, 点号由南向北递增。此处以L9线为例, 分析三分量解释技术在金属矿勘查上的应用效果。发射框边长为600 m× 600 m, 发射框中心点在x方向上对应16号测点, 在y方向上向东偏离测线160 m, 发射电流13 A, 关断时间为500 μ s。测点点距为50 m, 接收使用PEM探头(10 kHz), 接收面积3 850 m2。实测中采用右手坐标系, 定义z分量向下为正, x沿测线大号点方向为正, y方向垂直测线向西为正。L9线地质剖面及工作布置如图5所示。

图5 夏日哈木矿区L9线测点布置及剖面示意

图6~图8分别为利用本文算法计算获得的xyz分量视电阻率剖面。从图中可以看出, 在100~400 m深度范围内存在一低阻体, 在xz分量视电阻率剖面上, 该低阻的位置大致对应地质剖面上偏右侧的铜镍矿体; 在y分量视电阻率剖面上, 低阻体的形态类似于低阻地层, 向两侧延展较大, 但随点号减小, 电阻率值稍有增加, 表现出尖灭的趋势。总体而言, 三分量视电阻率剖面粗略地反映了矿体的形态和位置。

图6 夏日哈木矿区L9线x分量视电阻率剖面

图7 夏日哈木矿区L9线y分量视电阻率剖面

图8 夏日哈木矿区L9线z分量视电阻率剖面

图9为联合使用三分量数据进行一维反演所得到的反演剖面, 其中红线部分为矿体轮廓。反演共用时37 min, 总体的数据相对拟合误差为 16.2%。对比地质剖面图容易看出, 反演剖面较为准确地反映了两个矿体的位置和形态。反演剖面上存在两个向上“ 拱起” 的低阻异常, 分别对应着地质剖面上的两个矿体(图中红色框范围), 其中右侧的矿体较大, 距离地表较近, 在反演剖面上反映出的低阻异常区域范围较宽泛, 电阻率值较低; 而左侧矿体较小, 埋深相对右侧矿体较大, 在反演剖面上反映出的低阻异常区域较小, 与背景的电阻率差异相对右侧矿体小一些。在两个矿体之间的连接部分, 反演剖面上的低阻异常形态呈现了“ 下凹” 的趋势, 这和已知地质剖面的矿体形态是一致的。由于使用一维近似, 连接部分的厚度相对已知地质剖面偏大。同样由于一维近似的原因, 在矿体两侧的尖灭部分, 反演剖面上仍出现了具有一定厚度的低阻异常。

图9 夏日哈木矿区L9线反演剖面

5 结论

1) 利用迭代方法求取瞬变电磁响应的反函数, 可以进行无需分期、分区的全域瞬变电磁视电阻率定义。但由于瞬变电磁场的复杂特性, 全域视电阻率定义方法存在“ 无解” 和“ 偏解” 的问题。“ 无解” 和“ 偏解” 的视电阻率值可根据迭代残差和迭代电阻率值予以剔除, 并通过插值方法进行“ 补全” , 以保证视电阻率曲线或剖面的完整性。模型计算结果表明, 由于场的特性随几何位置而不同以及插值补全的影响, 定源回线的全域视电阻率剖面对于真实电性分布的反映更为粗略(如相比于中心回线的全期视电阻率剖面), 仅仅只能作为解释工作的一种参考。

2) 通过引入不同分量的权重因子可以有效地联合瞬变电磁三分量数据进行反演。模型计算结果表明, 在合理设计目标函数、约束电阻率搜索范围的基础上, 联合使用定源瞬变电磁的多分量数据进行反演, 能获得更准确、更高分辨率的反演结果。

3) 利用本文算法进行实际资料处理的结果表明, 定源装置的三分量全域视电阻率剖面能粗略地反映地下低阻体的形态和位置; 联合使用三分量数据的一维反演剖面能较为准确地反映出低阻体的位置和形态, 且具有较快的计算速度。

The authors have declared that no competing interests exist.

参考文献
[1] Poddar M. A rectangular loop source of current on multilayered earth[J]. Geophysics, 1983, 48(1): 107-109. [本文引用:1]
[2] Raiche A P. Transient electromagnetic field computations for polygonal loops on layered earths[J]. Digital Library Home, 2012, 52(6): 785-793. [本文引用:1]
[3] 周楠楠, 薛国强, 李梅芳, . 基于电偶极子近似的多边形回线源瞬变电磁响应[J]. 煤田地质与勘探, 2011, 39(4): 49-54. [本文引用:1]
[4] 翁爱华, 刘云鹤, 陈玉玲, . 矩形大定源层状模型瞬变电磁响应计算[J]. 地球物理学报, 2010, 53(3): 646-650. [本文引用:1]
[5] 席振铢, 刘剑, 龙霞, . 瞬变电磁法三分量测量方法研究[J]. 中南大学学报: 自然科学版, 2010, 41(1): 272-276. [本文引用:1]
[6] 武军杰, 张杰, 王兴春, . 基于等效磁偶极子的定源回线瞬变响应计算方法及视电阻率定义[J]. 煤田地质与勘探, 2013, 41(3): 68-71. [本文引用:1]
[7] 武军杰, 王兴春, 邓晓红, . 定源回线瞬变电磁 x 分量视电阻率计算方法[J]. 物探与化探, 2012, 36(4): 684-687. [本文引用:1]
[8] 戚志鹏, 李貅, 郭文波, . 瞬变电磁水平分量视电阻率定义[J]. 煤炭学报, 2011 (s1): 88-93. [本文引用:1]
[9] 戚志鹏, 智庆全, 李貅, . 大定源瞬变电磁三分量全域视电阻率定义与三分量联合反演[J]. 物探与化探, 2014, 38(4): 742-749. [本文引用:1]
[10] 翁爱华, 陆冬华, 刘国兴. 利用连分式定义瞬变电磁法全区视电阻率研究[J]. 煤田地质与勘探, 2003, 31(3): 56-59. [本文引用:1]
[11] 张成范, 翁爱华, 孙世栋, . 计算矩形大定源回线瞬变电磁测深全区视电阻率[J]. 吉林大学学报: 地球科学版, 2009, 39(4): 755-758. [本文引用:1]
[12] 殷长春, 朴化荣. 电磁测深法视电阻率定义问题的研究[J]. 物探与化探, 1991, 15(4): 290-299. [本文引用:1]
[13] 蒋邦远. TEM中心回线电动势求解全区视电阻率的实用特性[J]. 物探与化探, 2015, 39(3): 530-536. [本文引用:1]
[14] 杜玮, 凌锦兰, 周伟, . 东昆仑夏日哈木镍矿床地质特征与成因[J]. 矿床地质, 2014, 33(4): 713-726. [本文引用:1]