地球物理学报  2020, Vol. 63 Issue (2): 715-725   PDF    
三维时间域航空电磁有理Krylov正演研究
邱长凯1,2, 殷长春2, 刘云鹤2, 张博2, 任秀艳2, 齐彦福3, 蔡晶2     
1. 中国地质调查局发展研究中心, 北京 100037;
2. 吉林大学地球探测科学与技术学院, 长春 130026;
3. 长安大学地质工程与测绘学院, 西安 710054
摘要:常规的三维时间域航空电磁模拟通常采用隐式步长方法进行时间离散,需要几次矩阵分解和上百次右端源项回带,计算效率较低.为了提高正演计算效率,本文提出使用有理Krylov方法求解时间域电场扩散方程.首先使用非结构四面体网格进行空间离散,采用Nédélec矢量基函数近似四面体单元内的电场;然后基于有限元离散给出矩阵指数和矢量乘积表示的电场显式解;最后采用有理Arnoldi算法构造Krylov子空间内的正交基函数并进一步求解矩阵指数与矢量的乘积,直接得到任意时刻的电场解向量,避免步长离散过程.此外,本文还提出一种指数加权偏移参数优化方法,使得有理Arnoldi近似在瞬变衰减晚期具备更高的精度,从而降低Krylov子空间阶数并提高计算效率.通过和层状模型解析解的对比验证了有理Krylov方法的精度.针对三维异常体模型使用全局网格和局部网格剖分并和其他数值方法比较,进一步说明了有理Krylov方法的有效性.
关键词: 航空电磁      时间域      三维正演      有限元法      有理Krylov方法     
3D time-domain airborne electromagnetic forward modeling using the rational Krylov method
QIU ChangKai1,2, YIN ChangChun2, LIU YunHe2, ZHANG Bo2, REN XiuYan2, QI YanFu3, CAI Jing2     
1. Development and Research Center, China Geological Survey, Beijing 100037, China;
2. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China;
3. School of Geological Engineering and Geomatics, Chang'an University, Xi'an 710054, China
Abstract: Traditional Three-Dimensional (3D) time-domain airborne Electromagnetic (EM) modeling utilizes an implicit time discretization technique to calculate the EM responses in the time domain,which requires several matrix decompositions and hundreds of forward and backward substitutions. This reduces largely the modeling efficiency. To speed up the 3D forward modeling for time-domain Airborne Electromagnetic (AEM),we propose to solve the diffusion equation for the electric field using the rational Krylov method. The spatial discretization is completed by the unstructured tetrahedral grids,where lowest order Nédélec first kind vector basis functions are adopted to approximate the electric field inside each element. Then we derive the electric field solutions directly in terms of the product of a matrix exponential function with a vector. The rational Krylov method is used to solve the matrix function,where the orthogonal basis vectors in the rational Krylov subspace are constructed using the rational Arnoldi algorithm. Finally,the electric field is calculated for any desired time by the rational Arnoldi approximation,which avoids explicit or implicit time stepping. In addition,an exponential weight is incorporated into the rational Arnoldi approximation errors when optimizing the shift parameters,resulting in small approximation errors at late times. This further reduces the size of the rational Krylov subspace and improves the computational efficiency. We verify the correctness and accuracy of the proposed rational Krylov method by comparing the forward modeling results for a homogeneous half-space model with the semi-analytic solution. Furthermore,we prove the effectiveness of the presented rational Krylov method by comparing our solutions with time-domain finite-element and finite-volume solutions for a 3D abnormal model using global and local grids. The simulations of AEM responses for two typical 3D abnormal models demonstrate that the rational Krylov method can speed up the time-domain AEM forward modeling for as up to 7 times while keeping high modeling accuracy.
Keywords: Airborne electromagnetic    Time-domain    3D forward modeling    Finite-element method    Rational Krylov method    
0 引言

航空电磁法由于采用飞行装置作为搭载平台,特别适合地形复杂地区能源和资源勘查(殷长春等,2015孙栋华等,2017宁媛丽等,2019).时间域航空电磁法由于发射功率大、勘探深度大(赵越等,2017武欣等,2019),特别适合勘探深部隐伏矿体.经过70多年的发展,时间域航空电磁法已经在矿产、油气、环境工程和水资源勘查等领域发挥积极应用(殷长春等,2015).时间域航空电磁数据处理解释技术也由近似成像和反演(Macnae et al., 1991, 1998)、一维横向约束反演(殷长春等,2016)和拟二、三维空间约束反演(Viezzoli et al., 2008王琦等,2016殷长春等,2018)等发展到三维反演解释(Cox et al., 2010Haber and Schwarzbach, 2014Liu and Yin, 2016Liu et al., 2016).考虑到正演是反演的基础,一个高效(快速、高精度)的正演算法可为有效的反演算法提供前提.因此,快速高精度的三维正演算法十分重要.

目前主流的三维时间域航空电磁正演模拟算法主要分为两种,即频时转换和时域直接模拟.频时转换算法首先计算宽频带电磁响应,进而采用正余弦变换(殷长春等,2013)或GS变换(Stehfest, 1970朴化荣和殷长春,1987)将频率域场变换到时间域.频时转换法的优势在于直接使用比较成熟的频率域正演模拟算法,计算相对简单.因此,目前的三维反演大多使用积分方程法(Cox et al., 2010)或有限差分法(Liu and Yin, 2016)计算频域响应,然后变换到时间域得到时间域正演响应.其不足之处在于为了保证转换算法的计算精度,一般需要计算几十个频率域响应,计算效率较低.

时域直接模拟方法采用时间步长离散并通过迭代方式求解不同时刻的电磁响应,分为显式和隐式两种方法.显式方法通过前几个(取决于差分阶数)时刻的电磁响应计算当前时刻的场,而隐式方法需要求解当前时刻和向后时刻的方程.时域有限差分法是一种典型的显式步长方法.Wang和Hohmann(1993)首先提出使用修正的Du Fort-Frankel方法(Du Fort and Frankel, 1953)进行三维时间域电磁模拟.Commer和Newman(2004)使用交错网格和时域有限差分法并行求解电性源电磁响应.孙怀凤等(2013)将瞬变电磁回线源直接加入到安培环路定理中,并考虑发射电流下降沿,实现三维瞬变电磁时域有限差分正演.在时域有限差分法中,由于虚拟位移电流项(Chew, 1995)的引入,时间步长必须满足Courant稳定性条件, 因此受网格尺寸和最低电阻率的约束,初始步长非常小(Wang and Hohmann, 1993),造成效率相对较低.

相比于显式方法,隐式方法有时是无条件稳定的(Ascher and Petzold, 1998),进而可以使用较大的步长,因此节省针对每一时间步长求解方程的需求.后推欧拉(Backward Euler, BE)方法是一种典型的隐式方法,由于其无条件稳定(Ascher and Petzold, 1998; Um et al., 2010),能够压制早期的高频误差(Haber et al., 2002),近十年逐渐在时间域电磁模拟中受到重视.Haber等(2002)首先提出使用后推欧拉和有限体积法求解电磁势问题.Um等(2010)将直接求解技术引入到后推欧拉方法中求解时间域海洋电磁模拟问题.Yin等(2016a, b)提出使用对称正定电导率张量和任意源项离散技术模拟复杂介质中航空电磁响应,并研究了各向异性和地形对时间域航空电磁的影响.Zhang等(2018)基于法向电流和切向磁场连续的混合误差估计算子开展时间域航空电磁自适应模拟.虽然直接求解技术能够提高时域电磁模拟速度,然而大量右端项回带仍十分耗时.特别地,对于后推欧拉方法,初始步长和步长增长因子选择不当会导致计算误差(Li et al., 2018).使用较少的步长能够提高时间域电磁模拟的速度,却无法保证计算精度.

本文研究使用有理Krylov方法求解矩阵函数,避免时间步长离散.有理Krylov方法是一种Krylov子空间投影算法,因其子空间是由有理函数作为基函数而得名.其基本思想为,构建m维子空间ϕ,在ϕ中寻找矩阵函数f(M)b的一个近似fmfm可看成是f(M)b在子空间ϕ中的一个投影.由于有限元离散的矩阵通常为大型稀疏矩阵,而Krylov子空间算法避免直接计算函数f(M),并在构建子空间时充分利用矩阵的稀疏性进行矩阵向量计算,极大提高计算效率.Krylov子空间算法在电磁模拟问题中已获得应用.Druskin和Knizhnerman(1988)首先使用谱Lanczos方法求解三维电磁非稳态问题.而后,Druskin和Knizhnerman(1994)又提出使用多项式Krylov方法进行频域和时域电磁模拟.Börner等(2015)使用有理Krylov方法开展地面瞬变电磁模拟研究.周建美等(2018)结合拟态有限体积离散和有理Krylov模型降阶算法求解多频可控源电磁的快速计算问题.

本文引入非结构有限元和有理Krylov方法实现时间域航空电磁快速正演模拟.非结构有限元的优势在于能够剖分任意航空电磁发射源,实现在发射源、接收点和电性分界面处的局部加密.考虑到时间域航空电磁感应电动势随时间近似指数衰减,在选择Krylov子空间偏移参数时,本文使用指数权重函数对有理Arnoldi近似误差进行约束,提高有理Arnoldi近似在晚期时间道的精度.相比于均一误差优化方法(Börner et al., 2015),本文算法能够显著降低Krylov子空间的阶数,进而提高计算效率.本文首先介绍离散电场扩散方程的有限元算法,推导矩阵函数和向量乘积表示的电场解向量,着重阐述有理Krylov方法的基本原理以及偏移参数优化方法;进而通过与一维半解析解进行对比,验证有理Krylov方法的精度;最后针对典型三维地电模型,通过和直接时域有限元法比较,验证有理Krylov方法的有效性.

1 正演理论 1.1 控制方程和有限元离散

在航空电磁法中位移电流对电磁响应的影响可以忽略(Yin and Hodges, 2005),因此在有源区域内电磁场满足的麦克斯韦方程为(Jin, 2014):

(1a)

(1b)

其中,e(r, t)和h(r, t)分别是t时刻的电场和磁场,r是位置矢量,μ是磁导率,σ是电导率,Js(r, t)是外加电流源项.对(1a)式取旋度并代入(1b)式,可得时间域电场扩散方程为

(2)

本文使用矢量有限元方法对(2)式进行离散.为此,定义残差矢量R(r, t)为

(3)

对模拟区域进行非结构四面体网格剖分,并使用矢量基函数(Nédélec, 1980)近似单元内电场,即:

(4)

其中Ni(r)是矢量基函数,ei(t)是t时刻单元内第i个棱边的切向电场值.Nédélec矢量基函数自动满足电场切向分量连续且散度为零的条件.按照伽辽金方法,选择权重函数为矢量基函数,并使得残差R(r, t)最小可得:

(5)

使用格林定理对(5)式第一项进行分部积分,将(4)式代入(5)式,得到:

(6)

式中,ABS分别为整体刚度矩阵、质量矩阵及源项,在单元内的局部形式为

(7a)

(7b)

(7c)

对于一阶单元,AeBe可直接通过单元积分得到.源项Se的离散参考Yin等(2016a).

1.2 边界条件和初值条件

使用非结构网格可以在电磁场梯度大的区域(比如电性分界面)很好地实现局部加密,同时网格尺度可快速增加到扩边区域.为保证电磁场解的唯一性,在外边界上施加狄利克雷边界条件:

(8)

其中,n为外边界面的法向向量.

本文假设航空电磁系统发射波形为阶跃波.根据电磁感应定律,t=0时刻,发射线圈中的恒定磁场在空中和地下均不产生电场,因此电场的初值条件为

(9)

1.3 有理Krylov方法

由(6)式可得矩阵指数函数和向量乘积表示的电场解向量为(Börner et al., 2015):

(10)

其中,矩阵函数为ft(z)=exp(-tz).式(10)中的b向量可以通过Bb=S求解.

由于有限元离散的矩阵AB均为大型稀疏矩阵,直接计算矩阵函数ft(M)b非常耗时,而有理Krylov子空间算法可以用来构建类似矩阵函数的有效近似.为此,首先定义关于质量矩阵B的正交内积和范数为

(11)

式中H代表共轭转置.

给定矩阵M和向量b,以及m阶有理函数rm=pm/qm,且有理多项式qm是非奇异的,则m+1阶有理Krylov子空间定义为

(12)

时,有理Krylov子空间

(13)

其中,I是单位矩阵,分母qm的零点ξj为定义Krylov子空间的偏移参数.有理Krylov方法的核心在于有理Arnoldi近似算法.其本质是使用Gram-Schmidt方法(Ruhe, 1994)构建一系列有理Krylov子空间内的正交基v(累计m+1个).记起始基函数在第j步(jm)计算第j+1个关于B矩阵正交的向量vj+1,即:

(14)

其中,且hi, j>0使得

将基函数v写成矩阵形式为

(15)

式中,N为有限元离散的未知数个数,m为偏移参数个数(等于子空间阶数减1).则由基函数Vm+1可得矩阵函数ft(M)b的有理Arnoldi近似为

(16)

其中,.显而易见,N×N阶的矩阵指数函数ft(M)降阶为(m+1)(m+1)阶的矩阵函数ft(Mm+1).由于m远小于N,计算复杂度极大降低,降阶后的问题非常容易计算.

在构建正交基的过程中,求解m次线性代数方程组(14)式最耗时,主要计算量在于偏移参数ξj改变时的矩阵分解和m次右端项回带.当ξj保持不变时,线性代数方程组系数矩阵不变,此时若采用直接求解器求解,则矩阵分解结果可以重复利用,大大提高计算效率.

偏移参数ξj是构建正交基Vm+1的重要参数,如何选择ξj是有理Krylov方法的关键.有理Arnoldi近似误差存在最大值为(Börner et al., 2015):

(17)

其中,[αβ]是矩阵M的奇异值范围,rm(z)是f(z)的有理近似.因此,必定存在一组偏移参数使得有理Arnoldi近似误差小于设定值.根据Börner等(2015),当偏移参数ξj的总个数m保持不变时,l=2个不同的偏移参数重复m/2次具有最高的计算精度.另外,考虑到航空电磁场衰减很快,晚期和早期电磁响应的幅值差异可达8至9个数量级,有理Arnoldi近似在晚期应具有更小的绝对误差,以保证晚期道计算精度.由Nabighian (1991),感应电动势在晚期的衰减近似满足,因此本文将有理Arnoldi误差乘以权重函数w(t)=t2.5作为最优化问题的目标函数,以保证有理Krylov近似在晚期误差小、精度高.偏移参数优化问题可以表述为:寻找一系列ξj以及有理函数,使得加权误差最小,即:

(18)

本文针对进行偏移参数优化,取l=2,将不同偏移参数重复20次,累计产生m=40个偏移参数,生成41阶有理Krylov子空间.使用有理Arnoldi近似求解(18)式,给定矩阵为对角阵,为单位矢量,使用全局搜索得到加权误差最小的偏移参数分别为:=-104.20,其近似误差如图 1所示.

图 1 41阶有理Arnoldi近似误差 Fig. 1 Errors of rational Arnoldi approximation of order 41
2 算例分析

下面本文分别针对均匀半空间和三维异常体模型验证有理Krylov方法的计算精度和效率.计算环境为配备Intel (R) Core (TM) i5-4590 CPU @ 3.30 GHz中央处理器和16 GB内存的个人计算机.航空电磁系统发射线圈为外接圆半径是15.35m的正十二边形,匝数为1匝,发射波形为阶跃波形,电流为1安培;接收线圈位于发射线圈右上方,水平距离为10 m,垂直距离为20 m.

2.1 精度验证

首先,本文设计一个均匀半空间模型,通过与一维半解析解比较来验证有理Krylov方法的精度和有效性.均匀半空间和空气的电阻率分别为50 Ωm和106 Ωm.发射源中心位于(0, 0, -30)m处,接收点位于(10, 0, -50)m.考虑到航空电磁系统的影响范围(footprint)有限(Yin et al., 2014),将四面体网格剖分的核心区域设为1800 m×1800 m×1500 m.为了满足狄利克雷边界条件,在xyz方向各扩边20 km.考虑到发射源附近电磁场变化剧烈,空气与地表电性对比度大,因此在源和地表附近均对网格进行了局部加密,网格尺寸随着远离发射源而逐渐增加,最终生成13701个节点和80660个四面体单元,产生94622个未知数,四面体网格截面图如图 2所示.

图 2 均匀半空间模型四面体网格截面图 Fig. 2 Section of tetrahedral grid at y=0 m for a homogeneous half-space model

图 3给出有理Krylov法和使用Airbeo(Raiche et al., 2007)计算得到的一维半解析解对比结果.其中,一维半解析解是由频率域磁场响应经过余弦变换得到的.由图可以看出,有理Krylov算法结果和一维结果吻合很好,从早期到晚期相对误差均小于5%,证明了有理Krylov方法具有很高的计算精度.

图 3 均匀半空间模型有理Krylov结果和一维半解析解的比较 (a) dBz/dt响应;(b)相对误差. Fig. 3 Comparison of rational Krylov solution and 1D semi-analytic solution (a) dBz/dt response; (b) Relative error.

为了比较有理Krylov方法和后推欧拉方法的计算效率,本文针对后推欧拉方法设计三组时间步长,如图 4所示.其中起始步长均为10-4 ms,时间步长增长因子分别为2.67、4.25和7.60,步长增大次数分别为7、5和4次,总步长数分别为800(Scheme 1)、510(Scheme 2)和200(Scheme 3)步.

图 4 后推欧拉法时间步长随时间变化 Fig. 4 Temporal step varying with time of 3 backward Euler schemes

图 5表示针对均匀半空间模型三组时间步长计算得到的后推欧拉法计算结果和一维半解析解的相对误差.其中,步长1和步长2的计算误差满足精度要求,而时间步长3由于增长太快导致误差很大.对比图 3图 5的结果表明,通过加权偏移参数优化后,41阶有理Arnoldi近似可以得到和510至800个步长的后推欧拉方法相当的计算精度.表 1给出有理Krylov方法和三组时间步长的后推欧拉法的计算效率对比,可以看出有理Krylov法相对使用800个步长的后推欧拉方法计算效率提高了13.15倍,相对于使用510个步长的后推欧拉方法计算效率提高了8.28倍.

图 5 后推欧拉方法结果和一维半解析解的相对误差 Fig. 5 Relative errors of 3 backward Euler solutions against 1D semi-analytic solution
表 1 有理Krylov方法和后推欧拉法计算时间对比 Table 1 Computational time for rational Krylov method and 3 backward Euler schemes
2.2 全局网格和局部网格的比较

图 6给出一个三维水平板状体模型,本文分别使用全局网格和局部网格进行离散,研究有理Krylov法对于全局网格和局部网格的计算效率.水平板大小为200 m×200 m×50 m,顶部埋深50 m.均匀半空间电阻率为100 Ωm,水平板电阻率为1 Ωm.分别模拟y=0 m剖面上-400~400 m范围内共41个测点的时间域航空电磁响应,点距为20 m.如图 7a所示,剖分全局网格时,在所有发射源和测点附近均进行了网格加密,共生成50162个节点和314153个四面体单元,产生364461个未知数.相比之下,在进行局部网格剖分时,参见图 7b,只在某个源或测点附近进行网格加密,大大降低了网格数量,平均每套网格约生成105026个单元,产生121994个未知数.将本文有理Krylov法分别和基于有限元离散(使用和本文相同的网格, FETD)(Yin et al., 2016a)和有限体积离散(使用规则六面体网格, FVTD)(Ren et al., 2017)的后推欧拉方法的计算结果进行比较以进一步验证本文算法的精度.

2.2.1 计算精度对比

图 8分别给出利用有理Krylov法、时间域有限元法(Yin et al., 2016a)和时间域有限体积法(Ren et al., 2017)计算得到的航空电磁响应.由图可以看出,三组计算结果吻合很好,地下良导板状体反映明显.在扩散早期,感应电流聚集在地表,dBz/dt无异常;随着扩散时间的增加,感应电流被低阻水平板吸引,异常体上方dBz/dt响应相对变小;当感应电流穿过异常体时,低阻板内电流密度增大且占据主导地位,因此接收机处磁场变大,异常体上方出现明显的单峰异常;在扩散的晚期,感应电流越过异常体继续向下和向外扩散,电磁响应异常逐渐消失.

图 6 三维水平板状体模型 Fig. 6 Plate model in a homogeneous half-space
图 7 (a) 全局网格; (b)第21个源的局部网格.模型参数见图 6 Fig. 7 Comparison of global mesh (a) and 21st local mesh (b) for the model in Fig. 6
图 8 有理Krylov和时域有限元及时域有限体积计算结果对比 (a)全局网格;(b)局部网格. Fig. 8 Comparison of dBz/dt responses by rational Krylov method (solid line), finite-element (cycle) and finite-volume time-domain methods (cross) (a) Global grid; (b) Local grid.
2.2.2 计算效率对比

表 2给出有理Krylov法和后推欧拉法的计算时间对比.从表可以看出,对于后推欧拉方法,计算效率随着时间步长数增加而降低;对于全局和局部网格,有理Krylov方法计算效率相对于使用800个时间步长的后推欧拉方法分别提高3.96和9.15倍.结合图 8表 2可以看出,有理Krylov法和局部网格技术不仅能满足计算精度要求,而且能显著提高计算效率.因此,后面的典型异常体正演模拟中本文均采用局部网格技术对模型进行离散.

表 2 水平板模型有理Krylov法和后推欧拉法计算时间对比 Table 2 Computational time of rational Krylov method and 3 backward Euler schemes for the plate model
2.3 典型异常体模型电磁响应

本节对典型的球体和双倾斜板模型的时间域航空电磁响应进行模拟,进一步检验基于四面体网格离散的有理Krylov正演方法的有效性.

2.3.1 球体模型

球体模型如图 9所示,半径为40 m,中心位于地下60 m处.均匀半空间和球体的电阻率分别为100 Ωm和1 Ωm.测线位于球体正上方,发射线圈高度为30 m,接收线圈高度为50 m,测点距为10 m.使用局部网格技术对模型进行离散,网格的平均棱边数约为158000个.

图 9 三维球体模型 Fig. 9 Spherical model in a homogeneous half-space

图 10给出球体模型的dBz/dt响应.整体上异常形态与图 8给出的板状体响应特征相似,响应峰值位于球体上方.在0.1 ms之前,由于扩散深度较浅,dBz/dt响应只有微弱的异常.随着感应电流向下扩散,异常电流在低阻球体内聚集,在0.2 ms至0.9 ms时,异常体内的电流占主导地位,dBz/dt呈现正异常,且在一定时刻达到最大值;随着时间的增加,感应电流逐渐远离球体,在1.8 ms之后dBz/dt异常逐渐消失.

图 10 低阻球体模型dBz/dt响应 Fig. 10 dBz/dt responses of the spherical model in Fig. 9
2.3.2 双倾斜板模型

双倾斜板状体模型如图 11所示.板状体水平距离为60 m,顶板埋深均为20 m,顶面中心分别位于(-30, 0, 20) m和(30, 0, 20) m处,每个板的x方向宽度为20 m,y方向宽度为100 m,z方向垂直高度为100 m,倾角为60°.均匀半空间的电阻率为100 Ωm,两个低阻板的电阻率均为1 Ωm.

图 11 双倾斜板体模型 Fig. 11 Dual dipping plate model in a homogeneous half-space

图 12给出双倾斜板模型的dBz/dt响应.整体响应特征同前,但由于板的倾斜特征异常形态相对复杂.早期(0.088 ms之前)感应电流聚集在地表,无法分辨两个低阻板,dBz/dt呈单峰异常;随时间推移,电磁场向深部扩散,感应电流集中于良导板状体内且占据主导地位,异常呈现双峰特征,表征两个异常体的存在.随着电磁场向深部进一步扩散,异常响应再次呈现单峰,幅值减弱直至最后消失.

图 12 双倾斜板模型dBz/dt响应 Fig. 12 dBz/dt responses for the duel dipping plate model
2.3.3 计算效率比较

表 3给出有理Krylov法和后推欧拉法的计算时间对比.从表 3可以看出,当采用局部网格时,对于球体模型和双倾斜板模型,有理Krylov法分别比使用800个时间步长的后推欧拉方法提速8.03和7.37倍,同时比使用200个时间步长的后推欧拉方法也提速近2.80和2.70倍.由此,本文提出的有理Krylov法能够显著地提高时间域航空电磁的模拟速度.

表 3 典型异常体模型有理Krylov方法和后推欧拉方法计算时间对比 Table 3 Computational time of rational Krylov method and 3 backward Euler schemes for typical models
3 结论

本文使用有理Krylov法和非结构有限元成功实现时间域航空电磁三维正演模拟.在优化偏移参数时,使用指数权重函数对有理Arnoldi近似误差进行加权,使得在晚期道有理Arnoldi近似的绝对误差更小,从而降低Krylov子空间阶数,提高计算效率.和均匀半空间模型的半解析解对比验证了本文有理Krylov法的有效性.针对不同的异常体模型对本文算法和其他算法的计算结果进行对比表明,本文提出的有理Krylov方法在保证计算精度的同时可极大提高计算效率,说明有理Krylov法是进行时间域航空电磁正演模拟快速有效的方法.

致谢  感谢曼彻斯特大学数学学院Stefan Güttel博士在文章准备过程中提供的帮助.
References
Ascher U M, Petzold L R. 1998. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Philadelphia, PA: SIAM: Society for Industrial and Applied Mathematics.
Börner R U, Ernst O G, Güttel S. 2015. Three-dimensional transient electromagnetic modelling using rational Krylov methods. Geophysical Journal International, 202(3): 2025-2043. DOI:10.1093/gji/ggv224
Chew W C. 1995. Waves and Fields in Inhomogenous Media. New York: Wiley-IEEE Press.
Commer M, Newman G. 2004. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources. Geophysics, 69(5): 1192-1202. DOI:10.1190/1.1801936
Cox L H, Wilson G A, Zhdanov M S. 2010. 3D inversion of airborne electromagnetic data using a moving footprint. Exploration Geophysics, 41(4): 250-259.
Druskin V, Knizhnerman L. 1994. Spectral approach to solving three-dimensional Maxwell's diffusion equations in the time and frequency domains. Radio Science, 29(4): 937-953.
Druskin V L, Knizhnerman L A. 1988. Spectral differential-difference method for numeric solution of three-dimensional nonstationary problems of electric prospecting. Izvestiya, Earth Physics, 24(8): 641-648.
Du Fort E C, Frankel S P. 1953. Stability conditions in the numerical treatment of parabolic differential equations. Mathematical Tables and Other Aids to Computation, 7(43): 135-152. DOI:10.2307/2002754
Haber E, Ascher U, Oldenburg D W. 2002.3D forward modelling of time domain electromagnetic data.//SEG Technical Program Expanded Abstracts 2002. Society of Exploration Geophysicists: 641-644.
Haber E, Schwarzbach C. 2014. Parallel inversion of large-scale airborne time-domain electromagnetic data with multiple OcTree meshes. Inverse Problems, 30(5): 055011. DOI:10.1088/0266-5611/30/5/055011
Jin J M. 2014. The Finite Element Method in Electromagnetics. 3rd ed. Hoboken, NJ: Wiley-IEEE Press.
Li J H, Lu X S, Farquharson C G, et al. 2018. A finite-element time-domain forward solver for electromagnetic methods with complex-shaped loop sources. Geophysics, 83(3): E117-E132. DOI:10.1190/GEO2017-0216.1
Liu Y H, Yin C C. 2016. 3D inversion for multipulse airborne transient electromagnetic data. Geophysics, 81(6): E401-E408. DOI:10.1190/GEO2015-0481.1
Liu Y H, Yin C C, Ren X Y, et al. 2016. 3D parallel inversion of time-domain airborne EM data. Applied Geophysics, 13(4): 701-711. DOI:10.1007/s11770-016-0581-x
Macnae J, King A, Stolz N, et al. 1998. Fast AEM data processing and inversion. Exploration Geophysics, 29(2): 163-169.
Macnae J C, Smith R, Polzer B D, et al. 1991. Conductivity-depth imaging of airborne electromagnetic step-response data. Geophysics, 56(1): 102-114. DOI:10.1190/1.1442945
Nabighian M N. 1991. Electromagnetic Methods in Applied Geophysics: Volume 2, Application, Parts A and B. Tulsa: Society of Exploration Geophysicists.
Nédélec J C. 1980. Mixed finite elements in R3. Numerische Mathematik, 35(3): 315-341. DOI:10.1007/BF01396415
Ning Y L, Zhou Z Y, Jiang M Z, et al. 2019. Airborne anomaly characteristics and prospecting intention of the upper reaches of Xiagalai'aoyi river Pb-Zn deposit in Heilongjiang. 34(3): 1074-1080, doi: 10.6038/pg2019CC0153.
Piao H R, Yin C C. 1987. Calculation of transient E.M. sounding using the Gaver-Stehfest inverse Laplace transform method . Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 9(4): 295-302.
Raiche A, Sugeng F, Wilson G. 2007. Practical 3D EM inversion? the P223F software suite. ASEG Extended Abstracts, 2007(1): 1-5.
Ren X Y, Yin C C, Liu Y H, et al. 2017. Efficient modeling of time-domain AEM using finite-volume method. Journal of Environmental and Engineering Geophysics, 22(3): 267-278. DOI:10.2113/JEEG22.3.267
Ruhe A. 1994. Rational Krylov algorithms for nonsymmetric eigenvalue problems.//Golub G, Luskin M, Greenbaum A eds. Recent Advances in Iterative Methods. New York, NY: Springer, 149-164.
Stehfest H. 1970. Algorithm 368:Numerical inversion of Laplace transforms. Communications of the ACM, 13(1): 47-49.
Sun D H, Li H Y, Jiang M Z, et al. 2017. Applications of time domain airborne electromagnetic and aeromagnetic in Wulonggou gold deposit exploration area. Progress in Geophysics (in Chinese), 32(6): 2533-2544. DOI:10.6038/pg20170634
Sun H F, Li X, Li S C, et al. 2013. Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time. Chinese Journal of Geophysics (in Chinese), 56(3): 1049-1064. DOI:10.6038/cjg20130333
Um E S, Harris J M, Alumbaugh D L. 2010. 3D time-domain simulation of electromagnetic diffusion phenomena:A finite-element electric-field approach. Geophysics, 75(4): F115-F126. DOI:10.1190/1.3473694
Viezzoli A, Christiansen A V, Auken E, et al. 2008. Quasi-3D modeling of airborne TEM data by spatially constrained inversion. Geophysics, 73(3): F105-F113. DOI:10.1190/1.2895521
Wang Q, Zhang Q, Yu S B, et al. 2016. Quasi 2D holistic inversion of fixed-wing time-domain electromagnetic data. Progress in Geophysics (in Chinese), 31(3): 1173-1180. DOI:10.6038/pg20160334
Wang T, Hohmann G W. 1993. A finite-difference, time-domain solution for three-dimensional electromagnetic modeling. Geophysics, 58(6): 797-809. DOI:10.1190/1.1443465
Wu X, Xue G Q, Fang G Y. 2019. Development of helicopter-borne transient electromagnetic in China. Progress in Geophysics (in Chinese), 34(4): 1679-1686. DOI:10.6038/pg2019CC0312
Yin C C, Hodge G. 2005. Influence of displacement currents on the response of helicopter electromagnetic systems. Geophysics, 70(4): G95-G100. DOI:10.1190/1.1993710
Yin C C, Huang W, Ben F. 2013. The full-time electromagnetic modeling for time-domain airborne electromagnetic systems. Chinese Journal of Geophysics (in Chinese), 56(9): 3153-3162. DOI:10.6038/cjg20130928
Yin C C, Huang X, Liu Y H, et al. 2014. Footprint for frequency-domain airborne electromagnetic systems. Geophysics, 79(6): E243-E254. DOI:10.1190/GEO2014-0007.1
Yin C C, Qi Y F, Liu Y H, et al. 2016a. 3D time-domain airborne EM forward modeling with topography. Journal of Applied Geophysics, 134: 11-22. DOI:10.1016/j.jappgeo.2016.08.002
Yin C C, Qi Y F, Liu Y H. 2016b. 3D time-domain airborne EM modeling for an arbitrarily anisotropic earth. Journal of Applied Geophysics, 131: 163-178. DOI:10.1016/j.jappgeo.2016.05.013
Yin C C, Qiu C K, Liu Y H, et al. 2016. Weighted laterally-constrained inversion of time-domain airborne electromagnetic data. Journal of Jilin University (Earth Science Edition) (in Chinese), 46(1): 254-261. DOI:10.13278/j.cnki.jjuese.201601302
Yin C C, Zhang B, Liu Y H, et al. 2015. Review on airborne EM technology and developments. Chinese Journal of Geophysics (in Chinese), 58(8): 2637-2653. DOI:10.6038/cjg20150804
Yin C C, Zhu J, Qiu C K, et al. 2018. Spatially constrained inversion for airborne EM data using quasi-3D models. Chinese Journal of Geophysics (in Chinese), 61(6): 2537-2547. DOI:10.6038/cjg2018K0559
Zhang B, Yin C C, Ren X Y, et al. 2018. Adaptive finite element for 3D time-domain airborne electromagnetic modeling based on hybrid posterior error estimation. Geophysics, 83(2): WB71-WB79. DOI:10.1190/GEO2017-0544.1
Zhao Y, Xu F, Li X. 2017. Review on time-domain AEM system and applied potential. Progress in Geophysics (in Chinese), 32(6): 2709-2716. DOI:10.6038/pg20170656
Zhou J M, Liu W T, Liu H, et al. 2018. Research on rational Krylov subspace model order reduction algorithm for three-dimensional multi-frequency CSEM modelling. Chinese Journal of Geophysics (in Chinese), 61(6): 2525-2536. DOI:10.6038/cjg2018L0311
朴化荣, 殷长春. 1987. 利用G-S逆拉氏变换法计算瞬变测深正演问题. 物化探计算技术, 9(4): 295-302.
宁媛丽, 周子阳, 江民忠, 等. 2019. 黑龙江下噶来奥伊河上游铅锌多金属矿航空电磁异常特征及找矿意义. 地球物理学进展, 34(3): 1074-1080. DOI:10.6038/pg2019CC0153
孙栋华, 李怀渊, 江民忠, 等. 2017. 时间域航空电磁、航磁在五龙沟金矿整装勘查区中的应用研究. 地球物理学进展, 32(6): 2533-2544. DOI:10.6038/pg20170634
孙怀凤, 李貅, 李术才, 等. 2013. 考虑关断时间的回线源激发TEM三维时域有限差分正演. 地球物理学报, 56(3): 1049-1064. DOI:10.6038/cjg20130333
王琦, 张琼, 于生宝, 等. 2016. 固定翼时间域航空电磁数据拟二维整体反演. 地球物理学进展, 31(3): 1173-1180. DOI:10.6038/pg20160334
武欣, 薛国强, 方广有. 2019. 中国直升机航空瞬变电磁探测技术进展. 地球物理学进展, 34(4): 1679-1686. DOI:10.6038/pg2019CC0312
殷长春, 黄威, 贲放. 2013. 时间域航空电磁系统瞬变全时响应正演模拟. 地球物理学报, 56(9): 3153-3162. DOI:10.6038/cjg20130928
殷长春, 邱长凯, 刘云鹤, 等. 2016. 时间域航空电磁数据加权横向约束反演. 吉林大学学报(地球科学版), 46(1): 254-261. DOI:10.13278/j.cnki.jjuese.201601302
殷长春, 张博, 刘云鹤, 等. 2015. 航空电磁勘查技术发展现状及展望. 地球物理学报, 58(8): 2637-2653. DOI:10.6038/cjg20150804
殷长春, 朱姣, 邱长凯, 等. 2018. 航空电磁拟三维模型空间约束反演. 地球物理学报, 61(6): 2537-2547. DOI:10.6038/cjg2018K0559
赵越, 许枫, 李貅. 2017. 时间域航空电磁系统回顾及其应用前景. 地球物理学进展, 32(6): 2709-2716. DOI:10.6038/pg20170656
周建美, 刘文韬, 刘航, 等. 2018. 多频可控源电磁法三维有理函数Krylov子空间模型降阶正演算法研究. 地球物理学报, 61(6): 2525-2536. DOI:10.6038/cjg2018L0311