基于位移逆Krylov子空...磁法三维数值模拟和响应特征_董岩.pdf

第 37 卷第 1 期2023 年 2 月现代地质GEOSCIENCEVol.37No.1Feb.，2023DOI:10.19657/j.geoscience.1000 8527.2022.075基于位移逆 Krylov 子空间的全波形航空瞬变电磁法三维数值模拟和响应特征董岩，谭捍东，付兴(中国地质大学(北京)地球物理与信息技术学院，北京100083)收稿日期:2022-06-30;改回日期:2022-10-10。基金项目:国家自然科学基金项目(41830429);山西省重点研发计划项目(202102080301001)。作者简介:董岩，男，硕士研究生，1998 年出生，地球探测与信息技术专业，主要从事电法勘探算法研究。Email:dongyan_。通信作者:谭捍东，男，教授，博士生导师，1966 年出生，地球探测与信息技术专业，主要从事电法勘探理论及应用研究。Email:thd 。摘要:为了研究复杂地电模型的航空瞬变电磁法全波形响应特征，需要开发考虑发射波形的三维数值模拟算法。本研究基于非结构四面体网格和位移逆 Krylov 子空间(Shift-and-Invert Krylov，简称 SAI Krylov)方法，采用基于电偶极子离散的场源处理方法模拟场源，在时间域进行计算实现了全波形航空瞬变电磁法矢量有限元三维数值模拟。使用均匀半空间模型在阶跃波、半正弦波、三角波和梯形波激发下的全波形解析解、VTEM 实际激发波形的后推欧拉算法计算结果，检验了本研究开发的数值模拟算法的正确性。设计地表起伏异常体模型，计算和分析了航空瞬变电磁响应特征。开发的基于位移逆 Krylov 子空间的全波形航空瞬变电磁法三维数值模拟算法适合模拟复杂地电模型的响应，具有较高的计算精度。关键词:航空瞬变电磁法;三维正演;全波形;位移逆 Krylov 子空间;矢量有限元中图分类号:P631文献标志码:A文章编号:1000 8527(2023)01 0074 10Three-dimensional Full-time Airborne TEM Modeling with Shift-and-InvertKrylov Subspace Method and esponse Characteristic AnalysisDONG Yan，TAN Handong，FU Xing(School of Geophysics and Information Technology，China University of Geosciences，Beijing100083，China)Abstract:To study the fulltime response characteristics of airborne TEM with complex geoelectric models，it isnecessary to develop a 3D forward algorithm with transmitted waveform consideration Based on the unstructuredtetrahedral mesh and the Shift-and-Invert Krylov subspace(SAI Krylov)method，we achieve the fulltime 3Dedge finite-element simulation of airborne transient electromagnetic，by using the source processing technology-based dipole discretization technique We simulate the shape of the transmitting loop by topography，and thenthe instantaneous current pulse technology is used to simulate the transmitting waveform The numerical simula-tion algorithm developed here is verified by the analytical solutions of the uniform half-space model，under theexcitation of step wave，half-sine wave，triangle wave，and trapezoidal wave，and by the calculation results ofthe 3D FETD algorithm of the VTEM actual transmitting waveform The anomaly spherical model with upliftedterrain is designed，and the characteristics of airborne transient electromagnetic response are calculated and ana-lyzed The fulltime 3D numerical simulation algorithm of airborne TEM(based on SAI Krylov subspace)is suit-able to calculate the response of complex geoelectric models，and is highly accurateKey words:airborne TEM;3-D modeling;fulltime waveform;SAI Krylov subspace;edge finite-element0引言航空瞬变电磁法(Airborne Transient Electro-magnetic Method，简称 ATEM)是一种基于航空平台的瞬变电磁探测方法，具有速度快、成本低、通行性好的优点，可有效克服地面条件限制，在海洋、地形起伏剧烈、大面积森林或植被覆盖等地面通行困难的区域开展勘测工作，广泛应用于矿产、地下水勘探等领域，应用前景广阔1 3。ATEM 理想的发射波形是阶跃波，由于硬件的限制实际的 ATEM 仪器系统很难实现理想的发射波形4，通常会存在关断延时。三角波、半正弦波等发射波形在航空瞬变电磁数据采集中被广泛使用，有必要开发全波形的 ATEM 三维正演算法，并研究全波形 ATEM 电磁场的响应特征。殷长春等5 研究了层状介质的全波形航空瞬变电磁响应计算问题，通过频时转换首先得到层状介质的航空电磁系统的阶跃波响应，再将阶跃波响应与发射电流强度对时间的二阶导数进行褶积，得到了层状介质的全波形航空瞬变电磁响应。这一做法对于实际发射电流的全波形计算存在困难，因为实际发射电流采样率有限，难以可靠计算电流强度对时间的二阶导数值。孙怀凤等6 研究了三维介质的全波形航空瞬变电磁响应计算问题，将源的电流密度加入麦克斯韦方程中，采用有限差分法实现了全波形瞬变电磁三维数值模拟算法。时间离散采用显式差分格式，为了保证计算结果的准确性步长要足够小，导致计算时间相对较长。齐彦福等7 应用矢量有限元和后推 Euler方法研究了全波形航空瞬变电磁的三维正演。对比显式方法，隐式差分格式无条件稳定8，但是步长的选取仍要谨慎。为了提高三维正演效率，众多学者将 Krylov子空间方法引入到电磁场的计算中。Druskin、Knizhnerman 等9 10 最先将有理多项式方法引入时间域电场的直接求解。Brner 等11 应用位移逆Krylov 子空间(Shift-and-Invert Krylov，简称 SAIKrylov)。子空间方法得到频域解，然后通过频时转换得到时间域的解。Druskin 等12 将有理 Krylov子空间投影应用于求解麦克斯韦方程的三维扩散问题，并优化了有理子空间的偏移参数，使计算结果不再受限于模型电导率的差异。Brner 等13 通过有理 Krylov 子空间方法求得频域解，然后用快速 Hankel 变换得到时域解;并提出了一种新的替代优化方法来选择有理 Krylov 子空间的极点参数，该方法可以在一个确定的迭代次数内收敛，而不受网格尺寸和电导率结构的影响，但频时转换的精度受到频点个数影响。周建美等14 提出了频率域单极点模型降阶方法，使用拟态有限体积法实现了多频三维电磁问题的快速计算。邱长凯15 等提出了一种有理 Krylov 子空间的指数加权偏移参数优化方法，使得有理 Arnoldi 近似在瞬变衰减晚期具备更高的精度，结合矢量有限元离散得到了航空瞬变电磁的阶跃波响应。ZHOU 等16 开发了一种新的时间积分方法，引入二阶 函数，通过构造新的指数函数形式得到了全波形瞬变电磁响应的表达式，基于结构化网格有限体积离散和 SAI Krylov 子空间方法实现了全波形瞬变电磁三维数值模拟。为了更好地模拟起伏地表和复杂地电模型，本文基于非结构网格和 SAI Krylov 子空间方法，系统研究和实现了全波形航空瞬变电磁法三维数值模拟算法，验证了本文算法的精度，设计地表起伏异常体模型，计算和分析了其航空瞬变电磁响应特征。1基于 SAI Krylov 子空间的 ATEM 全波形三维数值模拟1.1ATEM 电磁场满足的微分方程和边界条件在研究区域 内，忽略位移电流和磁导率差异的影响，ATEM 电磁场满足以下微分方程17:e(r，t)=0h(r，t)t(1a)h(r，t)=e(r，t)t+js(r，t)(1b)其中，e(r，t)和 h(r，t)分别为位置 r 处 t 时刻的电场强度和磁场强度，为电导率，js(r，t)为源电流密度，0为真空中的磁导率。根据(1a)和(1b)式可以得到电场满足的扩散方程为:10 e(r，t)+e(r，t)t=js(r，t)t(2)在研究区域 的边界上，采用狄利克雷边界条件，取边界上的切向电场为零。e(r，t)|=0(3)ATEM 常用的激发波形有阶跃波、半正弦波、梯形波等。本文研究 ATEM 全波形数值模拟方法，全波形的定义见图 1，包含供电时间段(on-time)57第 1 期董岩等:基于位移逆 Krylov 子空间的全波形航空瞬变电磁法三维数值模拟和响应特征和断电时间段(off-time)。发射线框内的电流 I 在on-time 时段非零，而在 off-time 时段为零，toff为on-time 和 off-time 的交点。图 1任意发射波形示意图Fig.1Schematic diagram of the arbitrary transmitting waveform1.2ATEM 全波形三维矢量有限元数值模拟采用矢量有限元法18 求解微分方程(2)。采用非结构四面体对整个研究区域进行剖分，如图 2所示。电场定义在四面体单元的棱边上(图 3)。则第 k 个单元内的电场 ek(r，t)使用矢量基函数计算得到，具体表达式为:ek(r，t)=6j=1ekj(t)nkj(r)(4)其中，ekj(t)为第 k 个单元第 j 条棱边上的切向电场强度，nkj(r)为单元矢量基函数，j=1，2，6。在整个研究区域 采用伽辽金法对微分方程(2)进行离散化、计算单元矩阵、整体合成等处理(具体细节参考文献 19)，得到有限元控制方程:Ae(t)+Bde(t)dt=J(t)(5)图 2研究区域非结构网格剖分示意图Fig.2Schematic diagram of unstructured gridschematics in the study area图 3四面体单元棱边电场分布图Fig.3Edge electric field distribution(of tetrahedral ele-ment)diagram其中，A 为刚度矩阵，B 为质量矩阵，J(t)为外加电流源项，e(t)为所有剖分单元棱