大定源回线TEM地空系统全域视电阻率定义
赵越1, 王祎鹏1,2, 李貅1
1.长安大学 地质工程与测绘学院,陕西 西安 710054
2.中航勘察设计研究院有限公司,北京 100098
通讯作者:李貅(1958-),男,吉林长春人,博士,教授,研究方向为电磁法理论与应用。

作者简介:赵越(1989-),女 ,江西九江人,博士研究生,主要从事瞬变电磁地空系统、航空时间域电磁技术研究。

摘要

介绍了一种快速有效的瞬变电磁地空系统全域视电阻率定义方法。利用磁场垂直分量与水平分量的时间域响应定义地空系统全域视电阻率,实现全空间、全时域视电阻率定义。通过深入分析均匀半空间中地空系统的电磁场表达式,发现地空系统的时间域磁场响应可以用与地面类似的多项式表示,利用视电阻率平移算法,实现全域视电阻率定义。理论模型计算表明,该算法无需迭代,速度快且精度高。最后,对比分析了磁场强度法与感应电动势法定义视电阻率对于异常体的分辨能力。

关键词: 瞬变电磁地空系统; 视电阻率; 多分量; 全域
中图分类号:P631.3 文献标志码:A 文章编号:1000-8918(2015)02-0352-06
The Definition of Full-domain Apparent Resistivity Based on Air-ground Transient Electromagnetic Data
ZHAO Yue1, WANG Yi-Peng1,2, LI Xiu1
1.College of Geology Engineering and Geomatics, Chang 'an University, Xi'an 710054,China
2.AVIC Geotechnical Engineering Institute Co.,Ltd., Beijing 100098,China
Abstract

A fast and efficient numerical calculation method of full-domain apparent resistivity for AGTEM has been developed in this paper. We defined the full-domain apparent resistivity using the vertical and horizontal component of time domain magnetic field, in this way we can achieve full space, all time domain apparent resistivity definition. We analyzed thoroughly the Air-ground transient electromagnetic field analytic expression in homogeneous half space, and discovered that the Air-ground time domain response can be expressed in similar polynomial with the central vertical response, use the apparent resistivity translation algorithm, we can achieve full-domain apparent resistivity definition. Model calculation results show that the apparent resistivity definition method does not require iterative, with quicker speed and higher precision. Finally, compared the magnetic field and the voltage based apparent resistivity for resolution of abnormal body.

Keyword: Air-ground transient electromagnetic method system; Apparent resistivity; Multi-component; Full-domain

瞬变电磁地空系统是将电性或磁性发射源放置在地表, 利用无人机或无人飞艇携带探头在空中进行信号采集, 采用全域、高密度、扫面式的三维测量的装置。瞬变电磁地空系统结合了航空电磁法工作效率高和地面电磁法采集信号信噪比高的优点, 工作成本低, 操作安全且工作方式灵活多变, 是一种极具潜力的电法勘探新方法, 特别适合在山区、无人区、沼泽等地区进行详细的地质调查工作。

瞬变电磁地空系统这种工作方式的原型最早由Misac N. Nabighia[1]提出, 接着在日本与加拿大不断发展。20世纪90年代, 具有代表性的有FLAIRTEM系统和TerraAir系统[2, 3, 4], 之后R.S.Smith[5]在对地空系统、地面TEM及航空TEM进行对比之后, 发现对于深部目标体的探测, 地空系统优于航空TEM系统。对于电法勘探而言, 视电阻率参数仍然是目前广泛应用的参数, 通过引入适当的视电阻率定义可以使视电阻率响应曲线能够很好地反映出地下电性界面参数的变化情况。目前, 对于全区视电阻率定义问题, 已有许多学者做了大量的相关研究工作[6, 7, 8, 9], 但大多是基于中心回线单分量进行研究, 对于非中心点多分量视电阻率定义问题鲜有研究。事实上, 研究非中心点处多分量视电阻率定义问题是非常有意义的, 尤其是对于地空系统, 其工作方式通常是进行扫面测量, 即测量的是全空间内的电磁响应数据, 为了既保证工作效率, 又能全面利用采集数据, 有必要对全时域、全空间内的视电阻率进行定义。由此, 笔者针对一维均匀层状介质进行正演, 分析大定源地空系统非中心点处垂直分量与水平分量磁场强度的曲线变化规律, 并以此为基础, 提出了一种基于平移算法的地空系统全域(即全空间, 全时域)视电阻率定义方法。

1 时间域正演响应计算

以大定源装置为例(图1), Rx、Tx分别为接收回线与发送回线, r为偏移距, 接收高度为Z。通过推导得到均匀半空间大定源地空系统频率域电磁场表达式:

Eφ=-iωμFr=iωμI0a0λλ+u1eλzJ1(λa)J0(λr)d, Hr=2Frz=I0a0u1λλ+u1eλzJ1(λa)J1(λr),

Hz=-1rrrFr=I0a0λ2λ+u1eλz·J1(λa)J0(λr), (1)

式中:Eφ 为电场强度; HrHz分别为磁场强度的水平分量与垂直分量; I0为供电电流强度; a为圆形发射回线半径, 若发射回线是边长为L的方形回线时, 可通过两者面积相等来计算其等效半径(a=L/ π); r为偏移距; u1为瞬变电磁场参数, u1= λ2+k12, k12=-iω σ μ 0; σ 为均匀半空间电导率; μ 0 为真空磁导率。

图1 地空系统模型装置示意

从以上公式中可以看出, 在非中心点处地空系统的电磁场响应均为双重贝塞尔函数的积分, 由于该函数的强振荡与慢衰减特性难以应用通常的数值积分算法进行计算, 笔者采用贝塞尔函数展开法对其进行计算。这种方法不仅精度高[10] , 并且速度快:将含有双重贝塞尔函数的积分式子中的一个贝塞尔函数写进核函数, 则剩下的式子用线性数字滤波法进行计算。令

Fr=u1λλ+u1eλzJ1(λr), Fz=λ2λ+u1eλzJ0(λr),

则有 Hr=-I0a0FrJ1(λa), 2Hz=I0a0FzJ1(λa)3

频率域磁场和时间域磁场的对应关系为[11]:

H(t)=2π0ImH(ω)ωcos(ωt)H(t)=H0-2π0ReH(ω)ωsin(ωt)Hz(t)t=-2π0ImHz(ω)·sin(ωt)Hz(t)t=-2π0ReHz(ω)·cos(ωt)(4)

对于时间域磁场的求取, 根据频谱分析理论, 先在频率域中求解给定场源的电磁场, 然后通过傅氏反变换得到相应的时间域电磁场。

2 时间域电磁响应的平移性质

对于地面装置, 当在回线源中心点处接收时, 磁场垂直分量的时间域响应解析式为[12]

Bz(t)=I0μ02af(u)=I0μ02a1-3u2Φ(u)+2π3ue-u22, (5)Φ(u)=2π0u(t)e-t2/2dt; u=aμ02ρt, f(u)=1-3u2Φ(u)+2π3ue-u2/2

Φ (u)称为概率积分, f(u)为Bz的核函数, 又称为归一化响应函数。

可以将地表非中心点处电磁场的垂直分量和水平分量用与中心点处垂直分量相似的多项式来表示[13]。从式(1)中可以看出, 与地面装置响应表达式不同的是地空系统电磁响应的表达式中多了一个e的指数项eλ z, 相当于多了一个系数项(当接收高度与偏移距确定时, 指数项相当于一个常数)。同理, 可以将地空系统非中心点垂直分量与水平分量用于中心点垂直分量相似的多项式

Bp(t)=I0μ02acp1-cp21u2Φ(u)+2πcp31ue-u2/26

来表示, 式中的cp1、cp2、cp3为未知系数(系数的选取与接收高度和偏移距有关); p=z和p=r分别对应于磁场的垂直分量和水平分量。为了进一步讨论方便, 将式(2)改写成以下形式

Bp(t)(σ, a, t)=I0μ02afp(u)

王华军[9]指出均匀半空间的瞬变响应曲线具有平移伸缩的特性, 提出了一种半空间全区视电阻率的直接解法, 该方法速度快, 精度高。以此为出发点, 讨论该方法在地空系统全域视电阻率计算中的适用性。

σ* =, t* =Kt, K为常数, 则有

u* =aμ0σ* 2t* =aμ02Kt* =aμ0σt=u, Bp(t)(σ* , a, t* )=I0μ02afp(u* )=I0μ02afp(u)=Bp(t)(σ, a, t),

上式可改写成

Bp(t)(, a, Kt)=Bp(t)(σ, a, t)

由此可知, 在观测装置、接收高度、偏移距不变的条件下, 若均匀半空间的电导率变为原来的K倍, 则将观测时间延长至原来的K倍时, 观测到的瞬变响应的磁场值等于原瞬变响应的磁场值。这说明地空系统均匀半空间磁场具有平移性质。

3 大定源地空系统全区视电阻率定义
3.1 均匀半空间地空系统磁场响应特征

为了更加直观地的展现地空系统均匀半空间响应的平移伸缩性质, 将在同一高度接收, 位于同一偏移距处电导率分别为σ 的理论瞬变响应曲线绘制于单对数坐标下(图2), 模型采用方形回线发射, 线框边长200 m , 发射电流I0=1 A , 接收高度Z=-100 m 。从图中不难看出, 其中任意一条曲线必可由另外一条曲线平移得到, 这个曲线的平移轨迹曲线称为平移曲线, 计算两者斜率可得

k=Bp(t)(, a, Kt)-Bp(t)(σ, a, t)log(Kt)-logt=0logK=0

由此可知, 平移曲线是一条斜率为0平行于x轴的直线。

图2 地空系统均匀半空间中不同偏移距磁场垂直分量曲线的平移特性

图2显示, 在框内r=50 m处观测到的Bz 随时间的变化曲线是一个单调减函数, 这与中心点处的变化规律是一致的[14]; 在r=300 m(框外200 m)处观测到的Bz随时间的变化曲线并不是单调的, 其在早期为负值, 到晚期转为正值, 这种现象可以通过“ 烟圈效应” 来解释[15], 符号改变的时间可以通过 {t}μs=0.121{r}2m{ρ}Ω·m计算得到, 即符号改变的时间与偏移距和地下电阻率有关。

3.2 地空系统平移算法的具体实现

引入文献[9]中的平移截距概念, 对于地空系统磁场而言, 由于平移曲线平行于x轴, 可将平移截距定义为

b(t)=Bp(t)(σ, a, t),

即可认为此平移截距为磁场值。所有平移截距构成的曲线称为平移截距曲线。从图2中可知, 地空系统均匀半空间的平移截距曲线可分为两类。

(1)在中心点及框内偏移距较小处, 平移截距曲线为单值函数, 可以直接利用插值办法算出视电阻率的唯一解。

(2)框外偏移距较大且发生“ 极性反转” 现象时, 平移截距曲线为上凸的双值曲线, 其解并不是唯一的, 在转折点处存在“ 无解” 以及“ 双解” 现象。参照文献[9]中的处理方法, 以截距最大值对应的时间点tm为界, 将截距曲线分为左右两支, 然后让它们分别与理论平移截距曲线的左右两支相对应, 则每段曲线都是单值的。这样既可消除多解性的问题, 也可采用插值的办法来直接计算, 速度极快。

3.3 视电阻率定义的中解的存在性与唯一性问题

中心回线方式的磁场强度响应随时间变化的曲线是一个单值函数, 然而, 在大偏移距情况下, 磁场强度随时间的变化曲线却并不是单调的。野外工作时并不是均匀半空间的情况, 而是复杂的地电情况, 在观测参数固定的情况下, 实测截距参数曲线的极大值与理论截距参数曲线的极大值并非完全相等, 可能出现以下三种情况:① bobsmax(t)< bmax(t), 即出现视电阻率值的“ 双解” 现象, 导致全域视电阻率曲线并不连续, 给进一步解释带来困难; ② bobsmax(t)=bmax(t), 这是最理想的情况, 电阻率有唯一解; ③ bobsmax(t)> bmax(t), 超出了均匀半空间条件下所能达到的极大值, 出现“ 无解” 的现象, 导致测深视电阻率曲线不连续。

关于“ 双解” 和“ 无解” 问题, 主要存在于电阻率的过渡段, 是由于这一段的响应变化不单调造成的。如何解决这类问题, 已有许多学者做了大量的工作。苏朱刘[7]等在中心回线的全期视电阻率定义中引入了虚拟全区视电阻率和虚拟半径的概念, 利用感应电动势的最大值对观测感应电动势进行校正, 保证视电阻率在全时段有解, 但是该方法要求实测的响应曲线是完整的, 这就使得这种方法存在一定的局限性。白登海等[8] 借助瞬变场参数把整个瞬变过程分为早期阶段、早期到晚期的转折点、晚期阶段, 然后分别迭代求取早期、晚期视电阻率的精确值, 进而通过转折点构成一条完整的全程视电阻率曲线, 然而迭代方法速度慢, 且精度难以保证。王华军[9]指出可以通过在实测平移截距极大点附近设一个过渡区, 而过渡区内的视电阻率值, 通过插值方法得到, 以保证获得平滑、连续的全区视电阻率曲线。综上所述, 对于本文中过渡区的处理, 为保证全时段视电阻率曲线的连续性以便于成图和解释, 笔者沿用插值方法计算过渡区段的视电阻率值。

4 理论模型算例分析

依据上述算法, 针对水平三层地层几种典型模型进行计算, 并与感应电动势平移算法进行对比, 分析这两种方式对于异常体的反映情况。

模型参数设置如下:采用方形回线发射, 发射回线边长100 m, 发射电流I0=1 A, 接收高度Z=-100 m, 各层电阻率为ρ i(i=1, 2, …, n), 相应层厚为 hi(i=1, 2, …, n)。

4.1 磁场强度平移算法多分量模型验证

目前, 瞬变电磁法的研究方法, 不管是理论上还是在实际应用中, 大都只是局限在磁场的垂直分量上。但是, 随着对物探精度要求的提高, 单一对磁场垂直分量的反演解释已经不能满足需要, 因此研究瞬变电磁法多分量解释技术意义在于用以提高其综合解释精度, 增强其探测复杂地质构造的能力。

图3、图4分别为不同模型、不同偏移距下利用磁场垂直分量与水平分量求得的地空系统全域视电阻率曲线。从图4中可以看出, 视电阻率曲线在形态上很好地反应了各电性层参数, 对于中间层的反应, 良导层的效果比高阻层的效果更好; 其次, 垂直分量与水平分量的视电阻率曲线形态基本一致; 最后, 从图中可以看出, 曲线形态与偏移距无关, 即无论是垂直分量还是水平分量, 全域视电阻率都不受偏移距的影响。

图3 H型模型磁场强度平移法全域视电阻率曲线随偏移距变化情况

图4 K型模型磁场强度平移法全域视电阻率曲线随偏移距变化情况

4.2 磁场强度法与电动势法视电阻率曲线特征

以三层模型为例, 分别给出磁场强度与感应电动势视电阻率曲线, 分析两类方法视电阻率曲线表现特征, 并进一步分析其对于异常体的分辨能力。就此引入在数学领域经常使用的相对变化率的概念, 定义相对异常幅度:

相对异常幅度=计算异常-围岩背景值模型异常-围岩背景值

对于低阻异常体的分辨, 选用三层H型地电模型, 接收高度为100 m, 偏移距为30 m , 以垂直分量为例(表1、图5)。对于高阻异常体的分辨, 选用K型地电模型, 见表1、图6。

表1 H型、K型全区视电阻率相对异常幅度对照

图5 H型模型平移法全域视电阻率曲线对比

图6 K型模型平移法全域视电阻率曲线对比

通过以上两类模型曲线可以看出:

(1)两种视电阻率的定义方式都对异常体有所反应, 对于低阻异常体的反应明显优于高阻异常体, 相对感应电动势法, 磁场强度定义的视电阻率曲线形态简单, 易于判断曲线类型, 并且不存在Overshoot(当地下电阻率由高阻突然变为低阻时, 在界面附近出现的一个向高阻震荡的极大值)或Undershoot(当地下电阻率由低阻突然变为高阻时, 在界面附近出现的一个向低阻震荡的极小值)[16]

(2)中间层越厚, 视电阻率曲线的异常表现越明显, 相对异常幅度也越大, 并且感应电动势法的Overshoot现象也越明显。

(3)从表1、表2中可以看出, 两种方法计算得到的相对异常幅度值差异不大, 感应电动势的分辨能力稍强。

5 结论

受文献[9]的启发, 提出了一种适合瞬变电磁地空系统的全空间全域视电阻率定义方法, 实现在时间上不分早晚、在距离上不分远近的全域视电阻率定义。根据文中的模型计算结果和讨论, 总结出以下几点结论。

(1)利用磁场强度的平移性质对全域视电阻率进行定义, 使用该方法无需迭代, 计算速度极快, 大大提高了效率, 为地空系统大数据量的解释工作提供了基础。

(2)大定源地空系统全域视电阻率曲线能准确反映地下介质的电性分布, 垂直分量与水平分量的全区视电阻率曲线形态一致, 并且曲线形态均不受偏移距的影响。

(3)通过对比磁场强度法与电动势法视电阻率曲线, 发现利用感应电动势定义视电阻率存在不同程度的假极值现象, 特别是H型地层假极值现象尤为明显, 而磁场强度法则不存在此种现象, 能够更加真实反应地下的电性特征。但是, 用感应电动势定义视电阻率的分辨率更高, 对于异常体的反映能力更好, 实际工作中应根据工作需要选择合适方法。

The authors have declared that no competing interests exist.

参考文献
[1] Nabighian, Misac N, ed. Electromagnetic methods in applied geophysics[M]. SEG Books, 1988, 2. [本文引用:1]
[2] Verma S K, Mogi T, Allah S A. Response characteristics of GREATEM system considering a half-space model[C]//20th IAGA WG1. 2 Workshop on electromagnetic Induction in the earth, 2010. [本文引用:1]
[3] Smith R S, Annan A P, McGowan P D. A comparison of data from airborne, semi-airborne, and ground electromagnetic systems[J]. Geophysics, 2001, 66(5): 1379-1385. [本文引用:1]
[4] Mogi T, Tanaka Y, Kusunoki K, et al. Development of grounded electrical source airborne transient EM (GREATEM)[J]. Exploration Geophysics, 1998, 29(2): 61-64. [本文引用:1]
[5] 嵇艳鞠, 王远, 徐江, . 无人飞艇长导线源时域地空电磁勘探系统及其应用[J]. 地球物理学报, 2013, 56(11): 36-40. [本文引用:1]
[6] 殷长春, 朴化荣. 电磁测深法视电阻率定义问题的研究[J]. 物探与化探, 1991, 15(4): 290-299. [本文引用:1]
[7] 苏朱刘, 胡文宝. 中心回线方式瞬变电磁测深虚拟全区视电阻率和一维反演方法[J]. 石油物探, 2002, 41(2): 216-221. [本文引用:2]
[8] 白登海, Meju M A, 卢健, 等. 时间域瞬变电磁法中心方式全程视电阻率的数值计算[J]. 地球物理学报, 2003, 46(5): 697-704. [本文引用:2]
[9] 王华军. 时间域瞬变电磁法全区视电阻率的平移算法[J]. 地球物理学报, 2008, 51(6): 1936-1942. [本文引用:3]
[10] 许洋铖, 林君, 嵇艳鞠, . 航空时间域电磁法回线源有限差分初始场计算[J]. 电波科学学报, 2010(2): 259-264. [本文引用:1]
[11] 李貅. 瞬变电磁测深的理论与应用[M]. 西安: 陕西科学技术出版社, 2002. [本文引用:1]
[12] 朴化荣. 电磁测深法原理[M]. 北京: 地质出版社, 1990. [本文引用:1]
[13] Sun H, Li X, Li S, et al. Multi-component and multi-array TEM detection in karst tunnels[J]. Journal of Geophysics and Engineering, 2012, 9(4): 359. [本文引用:1]
[14] Raiche A P. Comparison of apparent resistivity functions for transient electromagnetic methods[J]. Geophysics, 1983, 48(6): 787-789. [本文引用:1]
[15] 李金铭. 地电场与电法勘探[M]. 北京: 地质出版社, 2005. [本文引用:1]
[16] Raiche A P, Spies B R. Coincident loop transient electromagnetic master curves for interpretation of two-layer earths[J]. Geophysics, 1981, 46(1): 53-64. [本文引用:1]