时间域航空电磁法2.5维有限元模拟

引用本文

强建科, 周俊杰, 满开峰. 时间域航空电磁法2.5维有限元模拟[J]. 物探与化探, 2015,39(5): 1059-1062.
QIANG Jian-Ke, ZHOU Jun-Jie, MAN Kai-Feng. Synthetic study of 2.5-D ATEM based on finite element method[J]. Geophysical and Geochemical Exploration, 2015,39(5): 1059-1062.
Doi:10.11720/wtyht.2015.5.29 复制到剪切板

Permissions

《物探与化探》编辑部

时间域航空电磁法2.5维有限元模拟

强建科¹, 周俊杰², 满开峰¹

1. 中南大学地球科学与信息物理学院,湖南长沙 410083

2.核工业北京地质研究院,北京 100029

作者简介: 强建科(1967-),男,副教授,主要从事电磁法正演模拟与反演成像研究工作。

收稿日期: 2014-11-08

基金: 国家自然科学基金项目(41174104,41472301)

摘要

采用三角网格剖分的有限元法,研究了2.5维航空瞬变电磁法正演模拟问题。利用时频变换数值方法将时间域电磁场转换到拉氏域,再利用傅里叶变换将三维问题降维变为2.5维问题,然后由有限元法求解得到拉氏域二维电磁场,逆拉氏变换后得到时间域航空瞬变响应。为了回避正演模拟中总感应磁场在场源处的奇异性问题,采用异常场算法,场源响应通过在微分方程中施加背景电磁场实现。由于瞬变电磁信号具有较大的动态范围,而且需要经过两次正、逆拉氏变换和傅里叶变换,每个环节的计算精度和速度要严格控制在较高的水平上,否则积累误差会非常大。模型计算表明均匀大地和层状大地模型解析解与数值解吻合很好。这证明该算法是正确可行的,可作为研究二维复杂地质体的方法手段。

关键词: 航空瞬变电磁; 2.5维正演; 有限单元法; 异常场法; 响应特征

中图分类号:P631 文献标志码:A 文章编号:1000-8918(2015)05-1059-04 doi: 10.11720/wtyht.2015.5.29

Synthetic study of 2.5-D ATEM based on finite element method

QIANG Jian-Ke¹, ZHOU Jun-Jie², MAN Kai-Feng¹

1. Shool of Central South University, Central South University, Changsha 410083, China

2. Beijing Research Institute of Uranium Geology,Beijing 100029, China

Abstract

Based on the triangular mesh dissection of the finite element method, this paper presents the 2.5-D airborne transient electromagnetic method forward modeling. Firstly the time domain electromagnetic partial equations are converted into laplace domain with the numerical time-frequency transform algorithm, then the 2.5-D expressions are obtained from 3-D case by adopting the fourier transform. The airborne data are acquired from the results of the digital inverse laplace transform of the 2-D electric and magnetic components in laplacian domain, which is computed using finite element method.Instead of the conventional algorithm, the anomalous field method is employed throughout the finite element process, so that the source singularity of total field is avoided. the source impact is embodied in the differencial equation by appling the background electromagnetic field item.The computing accuracy and efficiency would be managed rigorously due to the wide dynamic range of the transient singnal, as well as the two-time inverse laplace and fourier transforms. Otherwise the accumulative error will soar acceptably. The synthetic model experiments shows that the numeric modeling solutions match the analytical results of homogeneous and layered erath responses well, which also proves the effectiveness of the algorithm.

Keyword: ATEM; 2.5-D forward modeling; finite element method; anomalous field algorithm; response characteristic

Show Figures

近年来, 航空瞬变电磁法ATEM(Airborne Transient Electromagnetic)作为一种快速勘查方法得到了较快发展, 其优势是探测范围更广, 测深效率更高, 如茂密森林、崎岖山地、沼泽、湖区、浅海等地区, 都是航空瞬变电磁的用武之地, 而常规物探方法却无能为力。目前广泛应用于生产实践中的是基于电导率— 深度成像技术和一维层状大地反演, 如“ 镜像场源法” 或“ 烟圈理论” ^[1]、最大电流层模型^[2]、薄板逼近法 $\begin{matrix} ^{[34]} \end{matrix}$ 和伪层半空间模型^[5]等近似解释法, 后来发展了一维层状大地反演 $\begin{matrix} ^{[612]} \end{matrix}$ 。由于一维模型的局限性, 往往无法解释或错误解释二、三维地质体上的ATEM数据, 给野外勘查造成较大的损失。研究高维时间域航空电磁响应的难度较大, 进展缓慢。

Stoyer等^[13]于1976年首先研究了2.5D时间域磁偶源电磁响应正演模拟问题, 尽管模拟速度很慢, 时间很长, 但2.5D问题更符合实际野外情况^[14]。另一种方法是先求解多个频点的频率域电磁响应, 再通过Hankel数字滤波算法把频率域电磁响应转换到阶跃函数的时间域电磁响应, 由此间接实现2.5D时间域航空电磁正演模拟^[15]。近年来, 时间域航空电磁在高维正反演方面得到了很大发展 $\begin{matrix} ^{[1620]} \end{matrix}$ , 如Wilson等发展了2.5维和3维正反演技术, 提高了对实测数据的解释水平。在国内, 自2003年罗延钟等人开始研究时间域航空电磁方法以来, 已有不少学者进行了一维、二维、三维正演模拟 $\begin{matrix} ^{[2131]} \end{matrix}$ 和一维反演方面的工作 $\begin{matrix} ^{[912]} \end{matrix}$ , 但还远不能进行生产应用, 该领域的研究工作还需长期进行。

目前研究2.5D航空瞬变电磁法存在的最大问题是精度不够高, 速度太慢。此次研究借鉴了地面瞬变电磁的一些技术, 提高了场源附近的计算精度, 通过优化算法提高了计算速度。

1 理论

ATEM与常规TEM的区别在于发射源的位置不同, 此次研究采用异常场法, 以避开场源直接影响。设瞬变电磁场可分为背景场与异常场两部分:

$\begin{matrix} E = E_{b} + E_{a}, H = H_{b} + H_{a} (1) \end{matrix}$

式中E、H分别代表电场和磁场, 下标a代表异常(abnormal), 下标b代表背景(background)。这里异常是指目标模型下场的响应与背景模型下场的响应之差。根据Maxwell方程组, 在直角坐标系下, 假定y为走向方向, z为垂直向下方向, 拉氏傅氏域异常电磁场所满足的偏微分方程为

$\begin{matrix} \begin{matrix} \{\begin{matrix} - im E_{zx} - \frac{\partial E_{ys}}{\partial z} = - μs H_{xs} \\ \frac{\partial E_{xs}}{\partial z} + \frac{\partial E_{zs}}{\partial x} = - μs H_{ys} \\ \frac{\partial E_{ys}}{\partial x} + im E_{xs} = - μs H_{zs} \end{matrix}, \\ \{\begin{matrix} - im H_{zs} - \frac{\partial H_{ys}}{\partial z} = (εs + σ) E_{xs} + J_{xs} \\ \frac{\partial H_{xs}}{\partial z} - \frac{\partial H_{zs}}{\partial x} = (εs + σ) E_{ys} + J_{ys} \\ \frac{\partial H_{ys}}{\partial x} + im H_{xs} = (εs + σ) E_{zs} + J_{zs} \end{matrix}, (2) \end{matrix} \end{matrix}$

式中i为虚数单位, ε 代表介电常数, σ 代表介质电导率, μ 代表磁导率, s代表拉氏域变量, m代表傅氏域变量, J代表场源电流密度矢量, 其表达式为

$\begin{matrix} J_{s} = σ_{s} E_{b}, (3) \end{matrix}$

由此可见, 异常场方程组中的场源项是由背景场量及异常电导率共同体现的, 与实际场源无关。因此, 可假定实际场源在某一简单模型(如均匀大地或者层状大地)下产生的场为背景场, 求解该背景场所满足的偏微分方程, 即可得到各电磁分量的解, 再用方程(2)求解得到异常场, 最终用式(1)得到总感应场的数值解。

从拉氏傅氏域电磁场微分方程组(2)出发, 消去电磁x、z分量, 可得到y分量的偏微分方程组, 与之对应的Ritz泛函极值问题为

$\begin{matrix} \begin{matrix} J (E, H) = \underset{Ω}{\int \int} {a_{1} [{(\frac{\partial E}{\partial x})}^{2} + {(\frac{\partial E}{\partial z})}^{2}] + \\ a_{3} [{(\frac{\partial H}{\partial x})}^{2} + {(\frac{\partial H}{\partial z})}^{2}] + a_{4} E^{2} + a_{5} H^{2} + \\ 2 a_{2} (\frac{\partial E}{\partial x} \frac{\partial H}{\partial z} - \frac{\partial E}{\partial z} \frac{\partial H}{\partial x}) + 2 b_{4} E_{b} E + \\ 2 b_{1} (\frac{\partial E_{b}}{\partial x} \frac{\partial E}{\partial x} - \frac{\partial E_{b}}{\partial z} \frac{\partial E}{\partial z}) + 2 b_{3} (\frac{\partial H_{b}}{\partial x} \frac{\partial H}{\partial x} - \frac{\partial H_{b}}{\partial z} \frac{\partial H}{\partial z}) + \\ 2 b_{2} (\frac{\partial H_{b}}{\partial z} \frac{\partial H}{\partial x} - \frac{\partial H_{b}}{\partial x} \frac{\partial E}{\partial z})} + 2 b_{2} (\frac{\partial E_{b}}{\partial x} \frac{\partial H}{\partial z} - \frac{\partial E_{b}}{\partial z} \frac{\partial H}{\partial x}) dΩ + \\ \frac{1}{2} \int_{l} η^{H} H^{2} + η^{E} E^{2}) dl \Rightarrow \min 。 (4) \end{matrix} \end{matrix}$

式中:Ω 代表求解区域, l代表求解区域边界, E代表拉氏傅氏域异常电场y分量, H代表拉氏傅氏域异常磁场y分量与虚数单位i的乘积, a_n、b_n(n=1, 5)、η ^E、η ^H为一些系数。该泛函取得极值的条件为

$\begin{matrix} δJ (E, H) = 0 。 (5) \end{matrix}$

利用有限元法可求解该式。求解思路为:先将求解区域进行单元剖分, 在每个单元做插值后对泛函积分; 再对各单元求和, 将连续函数的泛函离散成节点上的函数值的泛函; 再根据泛函极值条件推出各节点处应满足的线性方程组, 进而解得各节点函数值; 之后, 再利用

$\begin{matrix} H_{zs} = - a_{1} \frac{\partial E}{\partial x} - a_{2} \frac{\partial H}{\partial z} (6) \end{matrix}$

可解得垂直异常磁场分量, 式中参数意义与前文相同。

2 模型计算

根据上述算法, 用Fortran语言编制了采用规则三角网格剖分的有限元计算程序, 并分析典型2.5维模型上中心回线装置的ATEM响应。

图1a是均匀半空间模型和层状介质模型的正演模拟结果。均匀半空间电阻率ρ =100 Ω · m, 层状模型电阻率为:ρ ₁=100 Ω · m, ρ ₂=10 Ω · m, ρ ₃=100 Ω · m; 层厚度h₁=50 m, h₂=100 m; 矩形发射线框边长a=b=20 m, 发射电流I=1 A, 接收线框等效面积为50 m², 高度与发射线框相同; 吊舱高度Δ h分别为0 m、10 m、20 m、30 m和80 m。

	Figure Option View Download New Window
	图1 感应电动势随收发线圈高度的衰减曲线

从图1b中看出, 均匀模型的感应电动势随飞行高度呈对数值线性衰减, 在1~1 000 μ s之间幅值衰减明显, 但在晚期, 信号衰减很小。对于层状介质而言, 瞬变响应总体特征与均匀半空间相似, 但在0.1~10 ms之间, 感应电动势明显增大, 出现向上拱起的现象, 证明该算法正确反映了第二层低电阻率特征。

复杂模型测试中, 图2a展示了倾斜板状低阻异常体模型, 目标体电阻率为0.1 Ω · m, 围岩电阻率100 Ω · m, 上顶板距地面6 m, 垂直方向延伸27 m, 水平方向偏移18 m, 飞行高度10 m, 矩形发射线框边长a=b=20 m, 发射电流I=1 A, 接收线框等效面积为50 m²。

	Figure Option View Download New Window
	图2 倾斜板状低阻异常体多测道剖面模拟结果

由多测道剖面图2b看出, 感应电动势曲线发生严重畸变, 而且极大值由早期测道(1.29 μ s)逐渐从右向左下方偏移, 在时间窗口129.42~514.53 μ s区段出现几个极大值, 然后峰值逐渐下降。感应电动势的这种变化正好反映出倾斜板状体倾斜的方向, 极大值出现在倾斜板的重心上方, 说明2.5D正演模拟程序能够在复杂模型上得到正确的瞬变响应。

3 结论

基于拉氏傅氏域的异常电磁场变分问题, 利用异常场求解思路实现了航空瞬变电磁法的2.5维有限元正演模拟, 回避了场源问题, 从而降低了求解难度。通过求解均匀模型的电磁场, 为有限元加载了场源信息。模型计算表明, 2.5D航空瞬变电磁响应不仅与一维模型解析解吻合, 而且能明显反映局部异常体的瞬变响应, 从多测道剖面图上还能够获得异常体的形状、产状等性质特征, 进一步证明本算法正确、稳定可行。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Nabighian M N. Quasi-static transient response of a conducting half space: An approximate representation[J]. Geophysics, 1979, 44(9): 17001705. [本文引用:1]
[2]	Fullagar P K. Generation of conductivity-depth parasections from coincident loop and in-loop TEM data[J]. Exploration Geophysics, 1989, 20: 4345. [本文引用:1]
[3]	Keating P, Crossley D. The inversion of time-domain airborne electromagnetic data using the plate model[J]. Geophysics, 1990, 55(6): 705711. [本文引用:1]
[4]	Zhdanov M S, Pavlov D, Ellis R. Localized S-inversion of time domain electromagnetic data[J]. Geophysics, 2002, 67(4): 11151125. [本文引用:1]
[5]	Huang H, Rudd J. Conductivity-depth imaging of helicopter-borne TEM data based on a pseudolayer half-space model[J]. Geophysics, 2008, 73(3): F115F120. [本文引用:1]
[6]	Huang H, Palacky G J. Damped least-squares inversion of time-domain airborne EM data based on singular value decomposition[J]. Geophysical Prospecting, 1991, 39: 827844. [本文引用:1]
[7]	Paterson N R, Grant F S, Misener D J. Development of computer software for geophysical interpretation, in Pye. E. G. , and Barlow R. B. ,Eds. , Exploration technology development program for the Board of Industrial Leadership and Development: Summary of research1981 —1983[J]. Ontario Geol. Surv. miscellaneous paper , 1983, 115: 5361. [本文引用:1]
[8]	Sattel D. Inverting airborne electromagnetic(AEM) data with Zohdy's Method[J]. Geophysics, 2005, 80: 7785. [本文引用:1]
[9]	李永兴, 强建科, 汤井田. 航空时间域电磁法一维正反演研究[J]. 地球物理学报, 2010, 53 (3): 751759. [本文引用:1]
[10]	强建科, 李永兴, 龙剑波. 航空瞬变电磁数据一维Occam反演[J]. 物探化探计算技术, 2013, 35(5): 501505. [本文引用:1]
[11]	罗勇. 时间域航空电磁一维正反演研究[D]. 成都: 成都理工大学, 2013. [本文引用:1]
[12]	毛立峰, 陈斌, 吕东伟. 测高数据不准时的直升机航空瞬变电磁一维反演方法理论研究[J]. 物探与化探, 2011, 35(3): 402405. [本文引用:1]
[13]	Stoyer C H, Greenfield R L. Numerical solutions of the response of a two dimensional earth to an oscillating magnetic dipole source[J]. Geophysics, 1976, 41: 519530. [本文引用:1]
[14]	Sugeng F, Raiche A, Rijo L. Comparing the time-domain EM response of 2-D and elongated 3-D conductors excited by a rectangular loop source[J]. Geomag and Geoelec, 1993, 45: 873885. [本文引用:1]
[15]	Raiche, A P. Modelling the time-domain response of AEM systems[J]. Expl Geophys, 1998, 29: 103106. [本文引用:1]
[16]	Annetts D, Raiche A, Sugeng F. The time-domain electromagnetic response of wedge-like structures[C]//ASEG 17th Geophysical conference and exhibition, Sydney 2004. [本文引用:1]
[17]	Wilson G, Raiche A, Sugeng F. Practical 3D AEM inversion using 2. 5D modeling[C]//AESC 2006, Melbourne, Australia, 2006. [本文引用:1]
[18]	Wilson G, Raiche A, Sugeng F. 2. 5D inversion of airborne electromagnetic data[J]. Exploration Geophysics, 2006, 37: 363-371. [本文引用:1]
[19]	Viezzoli A, Auken E, Munday T. Spatially constrained inversion for quasi 3D modeling of airborne electromagnetic data: an application for environmental assessment in the Lower Murray Region of South Australia[J]. Exploration Geophysics, 2009, 40: 173-183. [本文引用:1]
[20]	Meng Y, Li W, Zhdanov M, et al. 2. 5-D electromagnetic forward modeling in the time and frequency domain using the finite-element method[C]//SEG 69th Annual Meeting Extended abstracts, USA, 1999. [本文引用:1]
[21]	罗延钟, 张胜业, 王卫平. 时间域航空电磁法一维正演研究[J]. 地球物理学报, 2003, 46(5): 719-724. [本文引用:1]
[22]	林君, 稽艳鞠, 于生宝. ATEM-Ⅱ型瞬变电磁系统的实际应用[J]. 物探与化探, 2004, 28(6): 496-499. [本文引用:1]
[23]	王卫平, 王守坦. 吊舱式直升机频率域电磁系统在北京密云红光铁矿的勘查效果[J]. 物探与化探, 2006, 30(5): 420-426. [本文引用:1]
[24]	刘桂芬. 回线源层状大地航空瞬变电磁场的理论计算[D]. 长春: 吉林大学, 2008. [本文引用:1]
[25]	周俊杰. 航空瞬变电磁2. 5维有限元正演模拟[D]. 长沙: 中南大学, 2011. [本文引用:1]
[26]	强建科, 罗延钟, 汤井田, 等. 航空瞬变电磁法关断电流斜坡响应的计算[J]. 地球物理学进展, 2012. 2, 27(1): 345351. [本文引用:1]
[27]	许洋铖, 林君, 李肃义, 等. 全波形时间域航空电磁响应三维有限差分数值计算[J]. 地球物理学报, 2012, 55(6): 21052114. [本文引用:1]
[28]	龙剑波, 强建科. 航空瞬变电磁一维正演的连分式算法研究, [J]. 湖南科技大学学报: 自然科学版, 2013, 28(4): 8691. [本文引用:1]
[29]	殷长春, 黄威, 贲放. 时间域航空电磁系统瞬变全时响应正演模拟[J]. 地球物理学报, 2013, 56(9): 31533162. [本文引用:1]
[30]	王宇航. 吊舱式时间域直升机航空电磁2. 5维正演及响应曲线分析[D]. 成都: 成都理工大学, 2013. [本文引用:1]
[31]	田培培. 时域航空回线源的三维异常体电磁响应探测分辨率研究[D]. 长春: 吉林大学, 2013. [本文引用:1]