复合电偶源瞬变电磁张量测量
朱伟国1, 席振铢2, 韩雪1
1.广东省地质物探工程勘察院,广东 广州 510800
2.中南大学 地球科学与信息物理学院,湖南 长沙 410083

作者简介: 朱伟国(1987-),男,硕士学位,助理工程师,从事地球物理勘查工作。

摘要

复合电偶源瞬变电磁张量测量能够有效解决单一电偶源在复杂三维电性结构体的瞬变场响应失真的问题。采用复合电偶源瞬变电磁张量的坐标不变量 P2研究均匀大地介质和一维典型层状介质,计算结果表明,坐标不变量 P2参数不论在早期或晚期都能准确反映一维介质的电性规律。对于复杂的三维地质体,采用坐标不变量 P2以及张量电阻率的极大值 ρmax和极小值 ρmin组成的旋转椭圆表征三维电性结构的特征,用低阻长方体模型进行正演模拟,计算结果表明:不同时间窗口坐标不变量 P2的响应以及极大值 ρmax和极小值 ρmin的旋转椭圆,不但可以准确反映出地下三维结构的变化情况,而且能够直观地勾画出异常区域边界的变化范围,这为后期的资料解译提供了理论指导。

关键词: 复合电偶源; 张量电阻率; 坐标不变量; 旋转椭圆
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2016)03-0534-07 doi: 10.11720/wtyht.2016.3.15
A study of electric dipole transient electromagnetic tensor survey
ZHU Wei-Guo1, XI Zhen-Zhu2, HAN Xue1
1. Geologic & Geophysical Engineering Exploration Institute of Guangdong, Guangzhou 510800, China
2. Central South University, Changsha 410083, China
Abstract

It is difficult to detect the complicated 3D structures by only using a single grounded source, as it could only generate single direction current vector and will have distortion sometimes. Therefore, the authors use the electric dipole sources transient electromagnetic tensor method, which is suitable for detecting 3D structure. In combination with numerical models, the authors present the data with the tensor apparent resistivity It is found that in uniform half-space and layered earth, the coordinate invariant P2 has a good representation for the change of underground media. In 3D model, the coordinate invariant P2 and ellipse are formed by the coordinate invariant ρmin and ρmax, which are calculated for 3D resistivity models. The authors hold that profiles of LOTEM measurements across a 3D cube structure can be used to accurately reflect the changes of the 3D structure and sketch the anomaly boundary. All of these results provide theoretical guidance for further TEM data interpretation.

Keyword: multiple electric dipole sources; resistivity tensor; coordinate invariant; ellipse

电偶源瞬变电磁法又称为长偏移距瞬变电磁法(LOTEM), 它是利用接地电偶极子向地下供电, 在供电间歇期观测二次场变化情况的一种时间域电磁法 13。LOTEM一直都是采用单一的偶极发射装置作为供电场源, 但是这种发射场源仅能在地下建立单一的电流矢量, 适用于简单的层状介质, 难以适应探测复杂的三维结构体, 而实际地质情况往往是复杂的, 因此, 考虑采用多场源电偶源瞬变电磁法来满足实际需求。

瞬变电磁法最常用的解释方法就是把观测场值转换成视电阻率曲线, 然后进行反演解释, 而视电阻率定义的好坏与否取决于能否真实地反映出地下介质的变化情况。瞬变电磁法视电阻率定义是参照直流电法的视电阻率定义, 虽然含义相同, 但是半空间场值与半空间电阻率存在着复杂的隐性函数关系, 难以用解析法推导出视电阻率与场之间的显式反函数, 只能使用各种近似定义方法或精确定义近似计算的方法 46。早期、晚期视电阻率虽然计算方法简单, 但是在不满足条件时存在一段不确定的区域, 无法用具体函数表示, 严重影响数据解释。为了更连续、更全面地解释多个场源产生的张量数据, 文中引入张量视电阻率 67作为解释参数。

外国学者Bibby首次提出张量视电阻率的定义, 并且引入到多场源偶极子直流电法中, 数值模拟结果表明张量视电阻率在解释地下结构时有效地避免了标量视电阻率引起的歧义 78; Bibby、Hohmann研究应用视电阻率张量解释三维多源偶极子直流电阻率数据, 数值模拟得出多源测量结合直流张量视电阻率的解释在许多方面优于单极— 偶极法[9]; Li、Pedersen将多场源和张量视电阻率应用到直流电法中, 克服了单场源带来的场源依赖性[10]; Caldwell引入张量视电阻率到时间域电磁法, 通过数值模拟来分析三维地下介质的响应特征 1112。在国内, 关于电偶源瞬变电磁法, 国内学者研究了LOTEM从数值模拟到实际应用, 但均是采用单一场源, 且观测的信号也仅限于磁场的变化信号, 几乎没有涉及电场以及多源装置的情况 1314

在前人研究的基础上, 进一步开展多场源电偶源瞬变电磁法的研究。电偶源瞬变电磁张量测量的场源至少由两个不同极化方向的电偶源组成, 通过获取不同极化方向的电场信号增加信息量, 减少单一信号在探测复杂结构下的不确定性; 同时, 引入张量视电阻率到时间域电磁法中, 结合各个坐标不变量不同特性, 全面分析三维地质体张量测量的数据, 以期为实现瞬变电磁法张量测量提供一定的理论基础。

1 基本理论
1.1 张量视电阻率的定义

标量视电阻率定义常用于共线装置, 广泛应用在物探数据解释中, 然而在多个场源解释三维结构时, 标量视电阻率的定义会产生假异常而造成歧义[8]。为了更好地解决上述问题, 引出张量视电阻率的概念。

如图1所示, 当两个电偶极子源ABCD向地下供电, 根据欧姆定律的微分形式, 可得到大地表面视电阻率张量ρ a(t)的关系式 78

E(t)=ρa(t)J, (1)

式中:E(t)是水平总电场向量; J是电流密度向量; ρ a(t)是一个与时间相关的二阶张量, 它包含4个相互独立的非零元素, 在笛卡尔坐标系下可表示为

ρa(t)=ρxx(t)ρxy(t)ρyx(t)ρyy(t); (2)

而电流密度J是一个与时间无关、与电偶极子的几何形状以及场源电流I大小相关的向量。如图1所示, 在P点所测得的电流密度水平分量可表示为

JAB=I(ra/ra3-rb/rb3)/2π(3)

接地电偶源在零秒时刻向地下供入阶跃电流, E(t)表示在t时刻接收的总电场向量。该总电场是一次场和二次场的叠加, 随着时间推移, 二次场慢慢衰减, 总电场在晚期会趋向直流张量电阻率法对应的电场值。

根据式(1)仅能获取单个场源的电场信号, 不足以求取张量视电阻率的各个元素, 必须使用两个不同方向的电偶极子(如图1所示相互垂直的场源ABCD), 获取相应的电流密度以及总电场向量, 就能计算出视电阻率张量

ρa(t)=1(Jx1Jy2-Jy1Jx2)×Ex1(t)Jy2-Ey1(t)Jx2Ey1(t)Jx1-Ex1(t)Jy1Ex2(t)Jy2-Ey2(t)Jx2Ey2(t)Jx1-Ex2(t)Jy1, (4)

上式中的角标1和角标2, 分别指与场源ABCD相关的变量。

图1 张量电偶源场源布置

1.2 张量不变量

Bibby已经对张量视电阻率的坐标不变量的基本特性有过详细介绍, 该特性也可以应用到时间域电磁法, 下面对其特点进行简单介绍。

张量视电阻率矩阵各个值会随着坐标的变化而改变, 这会对数据解释带来不便, 因此, 参照大地电磁法阻抗张量的坐标不变量[15], 从张量视电阻率矩阵推出3个最基本的旋转不变量 78, 根据式(2)可推出这三个坐标不变量的表达式

P1=trace(ρa)/2=(ρxx+ρyy)/2, (5)P2=det(ρa)2=(ρxxρyy-ρxyρyx)2, (6)P3=(ρxy-ρyx)/2(7)

从上面公式可知, 3个坐标不变量P1P2P3分别由视电阻率矩阵的迹、行列式、差变化得出, 均具有反映地下电阻率信息的特征, 基于不变量的组合变化也可用于数据解释中, 其特征将在实例中进行讨论。

1.3 电阻率旋转椭圆

在任意时刻t, 瞬态视电阻率张量可以用相应的旋转椭圆 67来表述。旋转椭圆表示的是总场视电阻率沿电场向量方向随场源变化的轨迹。假设场源偶极子旋转, 总场视电阻率沿所测电场方向的轨迹就是一个椭圆, 其长轴和短轴的长度分别对应最大和最小总场的电阻率(ρ man, ρ min), 表达式如下:

ρmax=P12+P32+P12+P32-P22, (8)ρmin=P12+P32-P12+P32-P22; (9)

由于椭圆面积是由坐标不变量变化而成的, 因此它也具有坐标不变量的性质:

πρmaxρmin=πP22(10)

如图2所示, 当电场向量在方位角ϕ =(α -β )时, 总场视电阻率达到最大值; 而电场向量在方位角ϕ =(α -β +π /2)时, 总场视电阻率达到最小值。其中α 是对称主轴的方位角; β 表示电场向量与电流密度方向的夹角, 表达式分别为

tan(2α)=(ρxy+ρyx)(ρxx-ρyy), (11)tan(2β)=P3/P1(12)

2 均匀半空间的P2特征

在均匀半空间中, 张量视电阻率可用关于场源对称的极轴坐标系表示。在极轴坐标系中, 张量非主轴分量为0。在均匀半空间电阻率为ρ 的视电阻率张量可简化为二阶矩阵

ρa(t)=ρrrρρθrρθθ=1-f(r/δ)2001+f(r/δ)13

式中:rθ 分别表示极性坐标轴的径向和切向分量, δ 是与距离相关的变量, 函数f(r/δ )是与时间相关的电场准静态误差函数[16]

通过张量视电阻率的定义, 计算出电阻率为200 Ω · m均匀半空间以下坐标不变量值, 与T.G Caldwell的计算结果[12]进行对比, 结果见表1

表1 坐标不变量P2在半空间的对比

模型计算表明, 文中通过定义计算出的不变量P2与前人计算结果基本一致, 相对误差最大仅为 0.7%, 并且在整个时间间隔基本保持常数不变, 接近半空间电阻率。可见, 在均匀半空间下, 相比于传统视电阻率区分早、晚期, 张量不变量P2在有效时间很好地反映出了半空间电阻率。

3 层状介质的P2特征

当地下介质的电阻率随着深度的变化而变化, 且接地线源可近似为电偶极子时, 视电阻率张量矩阵的求取可近似为关于发射装置的圆对称问题。瞬态视电阻率张量在极轴坐标系下是关于发射装置对称的, 于是通过变换可以得出瞬态视电阻率张量是

关于对角线对称的。而根据张量不变量的特性, 当地下介质是一维的情况下(P3=0), 推导出的一维介质视电阻率矩阵是对称的。

对A型、K型、H型、Q型4种三层地电模型进计算, 场源为两个相互垂直的电偶极子, 长度均为1 km, 发射电流为10 A, 中心位于坐标原点; 在x=y=7 000 m处获取其水平总电场, 并用电流密度进行归一化得到不变量参数。各模型中, h1=500 m, h2是变化的, 分别为100、500、1 000 m; A型电阻率模型:ρ 1=10 Ω · m, ρ 2=100 Ω · m, ρ 3=1 000 Ω · m; H型模型:ρ 1=100 Ω · m, ρ 2=10 Ω · m, ρ 3=1 000 Ω · m; K型模型:ρ 1=10 Ω · m, ρ 2=100 Ω · m, ρ 3=10 Ω · m; Q型模型:ρ 1=1 000 Ω · m, ρ 2=100 Ω · m, ρ 3=10 Ω · m。图3给出了4种地电模型张量不变量P2以及早、晚期视电阻率随中间层厚度变化所得到的视电阻率曲线。

从图3可知, 在满足极限条件下, 早期视电阻率响应能反映出浅部的真实电阻率, 晚期视电阻率趋近于深部的电阻率, 但在不满足极限条件下, 早晚期视电阻率都无法反映出中间层的电性结构。而坐标不变量P2不分早、晚期在全区段有具体的定义, 基本反映出层状介质电阻率的变化情况。相比于早期视电阻率, 坐标不变量P2在早期准确地反映出浅部的真电阻率; 在晚期接近深部的电阻率, 受收发距的影响, 最后趋近于某一稳定值。根据4种三层模型的响应效果可见, 坐标不变量的响应曲线在中间低阻层的反映效果要好于中间高阻层。

图3 三层介质视电阻率曲线

4 三维介质的张量特征

实际中的地下结构比较复杂, 一般大地结构是三维地质体, 因此有必要讨论三维介质的张量特征。当地下介质是三维时, 任意接收点需要获取至少两个不同方向的场源才能全面表示电磁场响应, 因此, 视电阻率张量的引入使得数据分析更加直观。文中采用积分方程法 1722来模拟三维结构体。

4.1 数值模拟

如图4所示的三维地电模型, 在电导率为σ b的层状大地中存在一个电导率为σ =σ b+Δ σ b的异常体, 异常区域为D。假定大地中的磁导率是自由空间中的磁导率μ , 并忽略大地中的位移电流, 取时谐因子e-iω t。根据电磁场理论, 三维地电模型中电磁场满足麦克斯韦方程

图4 层状介质的三维模型

其中 jΔσaD区域的异常电流, 表达式为

    jΔσa=ΔσaE,  rD0, rD(16)

从上述公式可知, 模型的电场可以表示为由背景电导率σ b产生的背景电场Eb及由异常电导率Δ σ a产生的异常电场 EΔσa的叠加。其总场表达式为

E=Eb+EΔσa, (17)

其中, Eb是场源在层状大地产生的一次场, 可通过数值方法求取。根据张量格林函数, 可求出异常在rj处的异常场 EΔσa:

EΔσa=DG˙E(rj|r)·ΔσaE(r)dv;  r, rjD(18)

因此, 只要能求出区域D任意位置的异常场, 就可以得到总场响应

E(rj)=Eb+DG˙E(rj|r)·ΔσaE(r)dv,  r, rjD(19)

4.2 三维长方体模型

如图5所示, 假设电阻率为100 Ω · m的均匀半空间大地中存在一个大小为2 000 m× 2 000 m× 600 m、顶部埋深为600 m的三维长方体, 其电阻率为10 Ω · m, 以对数形式在1 ms到10 s之间取25个时间点。场源ABCD的中心均位于坐标原点, 距离异常中心水平距离约为4 km, 其长度均为 1 km, 发射电流为1 A; 其中测量点距为500 m, 线距为500 m, 均布置在异常区域水平地面, x方向的接收范围为-11 000~-3 000 m, y方向的接收范围为-3 500~3 500 m。

图6所示为坐标不变量P2在不同时刻的剖面以及相应视电阻率旋转椭圆的叠加。在1 ms的时刻, 由于扩散时间短, 电磁波未完全到达异常区域所在深度, 此时旋转椭圆处于变化初期, 在异常区域上方的椭圆与区域外基本一致, 主轴均沿径向指向场源; 随着扩散时间的增加, 在31.6 ms时, 不变量P2响应基本反映到异常区域的大小, 也较清晰地反映出异常的范围, 同时旋转椭圆开始总体沿圆的趋势变化, 并且在异常区域边界变化最大, 基本垂直异常边界, 总体指向异常体中心点。在100 ms后, 剖面曲线基本不再变化, 最低电阻率为30 Ω · m左右, 椭圆极化图也基本呈圆形状, 在异常区域范围指向异常边界。此后, 由于受场源收发距及电流大小的影响, 在瞬变电磁张量测量后期近似于直流电法的张量测量。通过不变量P2在不同时刻的剖面切片结合旋转椭圆的变化, 不仅能够准确地勾画出异常的大小变化以及边界的范围, 还可以直观地反映出模型不同时刻电阻率的变化情况。

图5 三维低阻长方体正演模型示意

图6 坐标不变量P2的时间剖面与旋转椭圆叠加

5 结论

相比于LOTEM法, 电偶源瞬变电磁法张量测量采用多场源来激发信号, 观测总电场信号, 并且辅以张量视电阻率作为解释参数, 联合多场源数据反映地下结构。文中首先分析了张量视电阻率在半空间下的特性, 结果表明张量视电阻率不分早期、晚期, 在全时间段内都有很好的定义, 所推导的不变量P2与真值误差小, 克服了常规标量视电阻率需要区分早晚期而带来了不确定性; 在层状介质的模型下, 不变量P2能真实地反映出层状介质电阻率的变化情况, 并且响应曲线在中间低阻层的反映效果要好于中间高阻层; 通过模拟三维模型得出视电阻率张量在不同时刻的剖面图以及相应的旋转椭圆, 结果表明, 不变量P2可重构地下介质的主要特征, 旋转椭圆对于地下介质的变化边界比较灵敏, 在可探测范围内, 呈垂直于物性边界趋势, 可见, 旋转椭圆为三维数据解释提供了一种直观、有效的方式。

The authors have declared that no competing interests exist.

参考文献
[1] 牛之琏. 时间域电磁法原理[M]. 长沙: 中南大学出版社, 2007: 8-10. [本文引用:1]
[2] 朴化荣. 电磁测深法原理[M]. 北京: 地质出版社, 1990: 10-27. [本文引用:1]
[3] Strack K M. Exploration with deep transient electromagnetics[M]. Amsterdam: Elsevier, 1992. [本文引用:1]
[4] Spies B R, Eggers D E. The use and misuse of a pparent resistivity in electromagnetic methods[J]. Geophysics, 1986, 51(7): 1462-1471. [本文引用:1]
[5] 严良俊, 胡文宝, 陈清礼. 长偏移距瞬变电磁测深的全区视电阻率求取及快速反演方法[J]. 石油地球物理勘探, 1999, 34(5): 532-538. [本文引用:1]
[6] 白登海, Maxwell A Meju, 卢健, 等. 时间域瞬变电磁法中心方式全程视电阻率的数值计算[J]. 地球物理学报, 2003, 46(5): 697-704. [本文引用:1]
[7] Bibby H M. The apparent resistivity tensor[J]. Geophysics, 1977, 42(6): 1258-1261. [本文引用:1]
[8] Bibby H M. Analysis of multiple source bipole dipole resistivity surveys using the apparent resistivity tensor[J]. Geophysics, 1986, 51(4): 972-983. [本文引用:1]
[9] Bibby H M, Hohmann G W. Three dimensional interprettation of multiple-source bipole-dipole data using the apparent resistivity tensor[J]. Geophysical Prospecting, 1993, 41: 697-723. [本文引用:1]
[10] Li X, Pederson L B. Controlled source tensor magnetitellurics[J]. Geophysics, 1991, 56(9): 1456-1461. [本文引用:1]
[11] Caldwell T G, Bibby H M. The instantneous apparent resistivity tensor: a visualization scheme for LOTEM electric field measurements[J]. Geophysical Journal International, 1998, 135: 817-834. [本文引用:1]
[12] Caldwell T G, Bibby H M, Brown C. Controlled source apparent resistivity tensors and their relationship to the magnetotelluric impedance tensor[J]. Geophysical Journal International, 2003, 151: 755-770. [本文引用:1]
[13] 陈明生. 电偶源瞬变电磁测深研究(一)——基本原理[J]. 煤田地质与勘探, 1999, 27(2): 55-59. [本文引用:1]
[14] 陈明生. 电偶源瞬变电磁测深研究(五)——实测感应电压转换成垂直磁场[J]. 煤田地质与勘探, 1999, 27(5): 54-57. [本文引用:1]
[15] Szarka L, Menvielle M. Analysis of rotational invariants of the magneto telluric impedance tensor[J]. Geophysical Journal International, 1997, 129: 133-142. [本文引用:1]
[16] 米萨克N 纳比吉安. 勘察地球物理电磁法(第一卷)[M]. 赵经祥, 王艳君, 译. 北京: 地质出版社, 1992: 7-40. [本文引用:1]
[17] Xiong Z H. Electromagnetic fields of electrical dipoles embedded in a stratified anisotropic earth[J]. Geophysics, 1989, 54: 1643-1646. [本文引用:1]
[18] Xiong Z H. Electromagnetic modeling three-dimensional structures by the methods of system iterations using equation[J]. Geophysics, 1992, 57(12): 1556-1561. [本文引用:1]
[19] Xiong Z H, Tripp A. A block iterative algorithm for 3-D electromagnetic modeling using integral equations with symmetrized substructures[J]. Geophysics, 1995, 60: 291-295. [本文引用:1]
[20] Zhsanov M S, Fang S. Quiasi linear approximation in 3D electromagnetic modeling[J]. Geophysics. 1996, 61(3): 646-665. [本文引用:1]
[21] Hursan G, Zhsanov M S. Contraction integral equation method in three-dimensional electromagnetic modeling[J]. Radio Science, 2002, 37(6): 1089-1101. [本文引用:1]
[22] 张辉, 李桐林, 董瑞霞, . 利用高斯求积和连分式展开计算电磁张量格林函数积分[J]. 地球物理学进展, 2004, 20(3): 667-670. [本文引用:1]