频率域航空电磁法人机交互式二维正演研究
1.
2.
The human-computer interaction two dimensional forward research of the frequency domain airborne electromagnetic method
1.
2.
责任编辑: 沈效群
收稿日期: 2017-09-21 修回日期: 2017-11-12 网络出版日期: 2018-08-05
基金资助: |
|
Received: 2017-09-21 Revised: 2017-11-12 Online: 2018-08-05
Fund supported: |
|
作者简介 About authors
谭林(1985-),男,硕士研究生,主要从事航空物探和电磁法研究工作。 。
研究了频率域航空电磁法的电磁场理论、边界条件、有限元单元法求解等内容,重点介绍了网格剖分原理,实现了人机交互网格剖分软件以及二维有限单元数值模拟算法。通过二维数值模拟计算,了解了圆柱体模型的响应曲线特征,提高了频率域航空电磁数据的解释水平。
关键词:
This paper researched the electromagnetic field theory of the frequency domain airborne electromagnetic method, the boundary conditions, the finite element method and so on. Mainly introduced the mesh generation principle, completed the human-computer interaction mesh generation software and two dimensional finite element numerical simulation algorithm. Understand the characteristics of response of the cylinder model by the simulation calculation. Improved the interpretation level of the frequency domain airborne electromagnetic data.
Keywords:
本文引用格式
谭林, 周道卿, 宁墨奂, 秦溱.
TAN Lin, ZHOU Dao-Qing, NING Mo-Huan, QIN Zhen.
0 引言
航空电磁法是基于岩石电性和磁性的差异,利用电磁感应原理,以飞行器为观测载具,实现地球物理探测的勘探方法。2002年,中国国土资源航空物探遥感中心引进了加拿大IMPULSE系统,先后在我国江西武夷山、广东河源、内蒙乌达和吉林东部地区开展了试生产工作[1],并且取得了初步的地质效果。该系统拥有水平共面、垂直同轴2组线圈,每组线圈同时发射和接收3个频率电磁信号,大大提高了信息采集量和工作效率。由于引进的时间不长,针对其采集的数据资料进行处理解释的方法有限,目前对于资料的解释仅限于一维反演,尚没有针对该系统的完整实用化的二维正反演软件系统,资料的解释水平相对滞后[1]。因此,有必要针对该系统特点进行二维正演研究,了解典型二维地电模型的航空电磁法的响应特征,提高航空电磁法的解释水平。
1 频率域航空电磁法二维正演理论
1.1 电磁场方程
对于航空电磁法,假设时间变化关系为eiωt,忽视位移电流的影响,总电场和磁场由麦克斯韦方程描述[4]:
式中:Mp是外加的磁性源,σ是二维电导率,μ0是自由空间的磁导率。对于点电磁源,Mp很难由一个离散公式准确地表达,求解二次场的Es和Hs更有利些。二次场的Es和Hs满足[4]:
式中:σa=σ-σ0,是异常电导率;σ0是正常或背景电导率;与σ0对应的一次电场Ep可以数值积分形式进行求解。
图1
沿走向方向对y做傅里叶变换,将式(2)从频率—空间域变换到频率—波数域中。式中“^”表示频率—波数域的结果。变换后经整理得到的频率—波数域的电磁场方程为[4]
其中,
式中:
因此,通过方程组(4)求解出
1.2 边界条件
为了通过式(4)求解电磁场分布,还需要在研究区域边界上给出特定的边界条件,该边界条件是否有效将直接影响数值模拟算法计算结果的精度以及计算效率。在实际情况中,由于电磁场是以指数形式衰减的,电磁场在距离源点无限远处衰减为零,如果在有限的范围内假设它衰减为零或某个特殊值,通常会得到错误或误差大的计算结果;如果将有限单元网格延伸到几个趋肤深度以外,必将增加网格单元数,从而导致矩阵方程变的太大而不好求解,影响计算效率。
考虑到我们用有限单元法要模拟的是二次场的分布,二次场是局部异常体产生的,其强度弱于一次场,衰减快;另外,航空电磁法是采用高频工作,场源是磁偶极子,场的影响和分布范围有限,因此,对于足够大的剖分区域,可以给出在剖分区域边界处二次场值为零的边界条件。
综上所述,频率域航空电磁法二维数值模拟对应的边值问题可以表示为:
1.3 有限单元法求解
1.3.1 有限单元法计算公式
微分方程(6)可以利用基于加权残差概念的有限元法求解。假设场V(x,z)满足公式(6),在边界∂Ω为狄利克雷边界条件的Ω区域中,它以LV=f来表示,其中L为微分算子,f为指定的源函数。区域Ω被剖分成NL个二维单元,第k个单元以坐标为(xi,zi)(i=1,2,…,n(k))的n(k)节点所定义。
V在第k个单元的解是节点值vi与在第k个单元的一组形函数ψi(i=1,2,…,n(k))的近似。Vk通过下式获得[5]:
利用这种V的近似,加权残差方程可以组合成一个矩阵方程[5]:
其中,
式中:U、V分别是网格节点上待求的
求解矩阵方程(8),可得到频率—波数域中的结果,再通过傅里叶反变换,即可得到频率域的结果。在进行傅里叶反变换过程中,为保证计算结果的精度,需合理选取波数。通过对均匀半空间和层状介质进行试算,在一个数量级内对数等间隔选取4个波数,波数分布范围为(0.000 1,1)[2]。在这些波数处由有限单元法计算对应的频谱值,再经过样条插值得到正弦或余弦变换系数处的频谱值,可以保证傅里叶反变换的计算精度。
1.3.2 方程组求解
在有限元分析中,用单元定位向量来集成整体刚度矩阵的方法虽然简单, 但在计算过程中需要存储整体刚度矩阵,即使用变带宽方法,仍然需要一个很大的矩阵来存储这些数据,占用很大的内存,这样解题的规模就受到了限制。在实际消元的过程中,不必等所有的单元全部组装进整体刚度矩阵后再进行,只要将与此结点相关的所有单元(即连接此结点的单元)组装进整体刚度矩阵就可以进行消元了。已完成消元的把消元结果写入文件,其占用的内存再存放其他数据, 解方程时再从文件中读出消元结果求结点值。这样组装单元刚度矩阵与消元就可以交替进行,这种方法叫做波前法。
采用有限元法求解时,往往离散模型划分的单元和相应的结点很多,得到的求解方程的阶数一般都很高,系数矩阵往往不能全部进入计算机内存。波前法是利用外存解决计算机内存容量不够时较好的解法。
波前法的特点是:刚度矩阵K和场源向量S不按自然编号进入内存,而按计算时参加运算的顺序排列;在内存中保留尽可能少的一部分K和S中的元素。计算过程可简单介绍如下:
1) 按单元顺序扫描计算单元刚度矩阵及等效结点载荷列阵,并送入内存进行集成。
2) 检查哪些自由度已集成完毕,将集成完毕的自由度作为主元,对其他行、列的元素进行消去修正。
3) 对其他行列元素完成消元修正后,将主元有关K和S中的元素移到计算机外存。
4) 重复1~3步骤,将全部单元扫描完毕。
5) 按消元顺序,由后向前依次回代求解。
2 网格剖分原理及人机交互式实现
有限单元法数值模拟需要将待求解区域进行合理地网格剖分,网格剖分质量直接影响电磁法二维正演的效率及精度,网格剖分越细,越接近真实,精度越高,但带来的问题是计算效率降低。剖分的合理性在电磁法二维正演问题中至关重要,本次在自动化离散二维地电模型的基础上,引入了人机交互式剖分的方法,使得网格剖分更加高效,更加合理。
2.1 剖分原理
首先通过自动化的方式对模型进行离散化,生成离散化网格;然后根据需要进行网格的精细化剖分;最后根据边界条件进行边界扩展处理,形成剖分最终的有限单元网格。
2.1.1 模型的离散化
离散的基本思路是从模型各线段的起始点开始,根据给定的网格单元(dx,dy),当内部域面积大于整个单位域面积的二分之一时,将整个单位域划为内部域,否则划分为外部单位域。图2为异常体的剖分示意,离散化时存在4种情况。
图2
图3
第二种情况:S1位于左边,如图4所示。根据第一种情况推理,有
图4
第三种情况:S1位于上边,如图5所示。根据第一种情况推理,有
图5
第四种情况:S1位于右边,如图6所示。根据第一种情况推理,有
图6
从模型的首个线段开始,判断该线段处于的单元网格属于哪一种情况,再利用式(5)~式(8)即可判定此单位域划分为内域还是外域,然后依次判断下一线段直到最后一个线段结束为止。离散后的结果如图7所示。
图7
2.1.2 网格的精细化剖分
上述自动化的离散化方法仅完成模型的离散化,为降低计算误差,还需进行精细化剖分。精细化剖分是在上述自动离散化形成的网格基础上,利用人机交互式方法进行,以满足各种复杂地电模型的需要。
2.2 人机交互式实现
为实现人机交互功能,采用了QtGui图形用户界面库的QGraphics View Framework框架[6],以 QGraphicsItem类为基类,设计组成网格的子类:
1) GeoBlockItem类:网格单元类,以颜色代表此单元的电阻率值。
2) GeoLineItem类:网格线类,该类数据成员中包括了该网格线一系列的GeoBlockItem对象。
以QGraphicsView类为基类,设计GeoDivisionView类:该类显示网格,且实现了人机交互功能,包括网格单元的电阻率值修改,网格加密、抽稀等。
图8
图9
3 计算实例
以圆柱体模型为例(图10)。该圆柱体截面直径为80 m,顶面埋深20 m,电阻率50 Ω·m;围岩介质电阻率500 Ω·m。计算参数设置为:垂直同轴装置,飞行高度30 m,收发距6.5 m;O点为坐标原点。
图10
图11
图12
1) 虚部大于实部;随着频率变低,异常数值和异常幅度逐渐减小。
2) 异常剖面曲线以原点O对称。据此可判断圆柱体中心水平位置。
3) 高频(f=21 750 Hz)时,实部和虚部均为单峰曲线,单峰曲线两侧出现幅度不大的极小值;中频(f=4 350 Hz)时,实部和虚部均为单峰曲线,单峰曲线两侧不出现极小值;低频(f=870 Hz)时,虚部为单峰曲线,实部为双峰曲线,但异常幅度不大。
4 结论和建议
通过引入人机交互式方法,提升了网格剖分的效率,并实现了频率域航空电磁法人机交互式二维有限单元数值模拟算法,提高了频率域航空电磁数据的解释水平。今后将继续对人机交互式软件进行改进,解决数值模拟算法的效率问题,提高频率域航空电磁法二维正演的实用性,为二维反演奠定基础。
参考文献
一种矩形网格自动剖分及加密算法
[J].
带地形的可控源音频大地电磁法二维正演
[J].
2.5维电磁响应的有限元模拟与波数取值研究
[J].
DOI:10.3969/j.issn.1001-1749.2008.02.013
URL
利用有限元法实现了任意方向偶 极子源在二维介质中频率域电磁响应的数值模拟,研究了波数取值对模拟结果的影响。通过对构造走向的Fourier变换,将全三维电磁问题,转化为一系列二 维问题,并在波数域求解,极大地减小了计算工作量,导出了波数域耦合适用于二维电性介质中任何方向电或磁偶极子响应计算的电磁场方程。针对每个给定的波 数,上述耦合电磁场方程用等参有限元方法在x-z平面内求解。采用Fourier逆变换,将波数域解积分,得到空间域电磁场。针对电磁模拟计算中,源点的 奇异性,采用具有一定面积的伪δ函数表达源电流分布,使数值解精度得以提高。另外,采用等参有限元,使地下复杂地质体得到准确表达。利用不同波数值对均匀 介质与层状介质的模拟结果与解析解的对比,验证了算法的正确性与精度。利用层状介质模型的解析解与数值计算结果的对比,分析了波数的优化取值范围及取值点 数对数值模拟结果的影响,考察了算法对非均匀介质的适应性。
复杂介质有限元法2.5维可控源音频大地电磁法数值模拟
[J].
DOI:10.3321/j.issn:0001-5733.2004.04.026
URL
利用有限元 2 .5维可控源音频大地电磁法 (简称CSAMT)数值模拟方法 ,对 10 0Ωm均匀半空间介质中有限长度的电偶极源产生的电场、磁场及视电阻率、相位特征进行了数值模拟 ,研究了场的空间变化规律 .在一个象限中 ,场的特征存在双叶现象 ,当收发距大于 4个趋肤深度时 ,电阻率较接近介质真实的电阻率 ,这些结果为观测系统和收发距的选择提供了依据 .波数域场的特征表明 ,低波数对源的贡献占较大的比例 .有限元法 2 .5维CSAMT数值模拟的优势在于能较准确地获得复杂介质结构的波场特征 .本文结合直立异常体、倾斜异常体及断陷模型对CSAMT电阻率、相位剖面特征及频率曲线特征的可靠性进行了研究 .数值模拟结果直观地给出了异常体的剖面异常形态 .通过对比研究异常体的剖面异常形态和半空间场的特征进一步说明本文方法和软件在模拟复杂介质结构场特征时是可靠的 .这为认识观测数据 ,指导反演解释提供了较好的依据
航空电磁法的发展现状
[J].
DOI:10.6053/j.issn.1001-1412.2006.1. 009
URL
航空物探方法具有低成本快速度,从宏观上探测大区域地球物理场的优势。在我国航空电磁法相对于航磁而言发展较为缓慢,但航空电磁法的物理原理又使它具有其他航空物探方法所不能替代的功能和优势(如寻找地下水、环境调查等)。文章从航空电磁观测系统、航空电磁数据处理和解释以及航空电磁应用等方面阐述了2000年以来国内外航空电磁的最新发展动态,在此基础上分析总结了航空电磁未来发展的趋势和重点。
有限元法在大地电磁测深正演计算中的应用与改进
[J].
航空电磁法的应用现状和前景
[J].本文介绍了航空电磁法近年来在应用方面的进展情况,指出了航空电磁/电阻率填图的方法技术问题及其解决途径,并概述了航电未来的发展前景:
频率域线源大地电磁法有限元正演模拟
[J].
DOI:10.3321/j.issn:0001-5733.2006.06.035
URL
本文介绍了频率域线源大地电磁法有限元正演模拟的研究结果,在外边界上统一应用适合于人工源 的一阶吸收边界条件来形成边值问题,可减小基于平面波假设造成的人为截断边界的影响,程序编辑中设计了两个二维数组分别存储总体系数矩阵的非零元素和在总 体结点编号中的位置,使内存占用量减少,且物理意义明确,方便用高斯-赛德尔等迭代法解有限元方程时调用.采用视δ函数模拟线源,提高了解方程组的稳定 性.最后通过对1个简单模型和1个复杂模型的模拟,证明所用的方法对异常体能够有明显的反映,说明了该方法的可靠性和有效性.
Electromagnetic methods in applied geophysics
[J].
/
〈 |
|
〉 |
