董岩等：基于位移逆Krylov子空间的全波形航空瞬变电磁法三维数值模拟和响应特征

　　基于位移逆Krylov子空间的全波形航空瞬变电磁法三维数值模拟和响应特征

　　董岩,谭捍东,?付兴

　　中国地质大学(北京)地球物理与信息技术学院,北京 100083

　　摘? ? 要

　　为了研究复杂地电模型的航空瞬变电磁法全波形响应特征,需要开发考虑发射波形的三维数值模拟算法。本研究基于非结构四面体网格和位移逆Krylov子空间(Shift-and-Invert Krylov,简称SAI Krylov)方法,采用基于电偶极子离散的场源处理方法模拟场源,在时间域进行计算实现了全波形航空瞬变电磁法矢量有限元三维数值模拟。使用均匀半空间模型在阶跃波、半正弦波、三角波和梯形波激发下的全波形解析解、VTEM实际激发波形的后推欧拉算法计算结果,检验了本研究开发的数值模拟算法的正确性。设计地表起伏异常体模型,计算和分析了航空瞬变电磁响应特征。开发的基于位移逆Krylov子空间的全波形航空瞬变电磁法三维数值模拟算法适合模拟复杂地电模型的响应,具有较高的计算精度。

　　关键词

　　航空瞬变电磁法;?三维正演;?全波形;?位移逆Krylov子空间;?矢量有限元

　　0 引? ? 言

　　航空瞬变电磁法(Airborne Transient Electromagnetic Method,简称ATEM)是一种基于航空平台的瞬变电磁探测方法,具有速度快、成本低、通行性好的优点,可有效克服地面条件限制,在海洋、地形起伏剧烈、大面积森林或植被覆盖等地面通行困难的区域开展勘测工作,广泛应用于矿产、地下水勘探等领域,应用前景广阔[1-3]。

　　ATEM理想的发射波形是阶跃波,由于硬件的限制实际的ATEM仪器系统很难实现理想的发射波形[4],通常会存在关断延时。三角波、半正弦波等发射波形在航空瞬变电磁数据采集中被广泛使用,有必要开发全波形的ATEM三维正演算法,并研究全波形ATEM电磁场的响应特征。

　　殷长春等[5]研究了层状介质的全波形航空瞬变电磁响应计算问题,通过频时转换首先得到层状介质的航空电磁系统的阶跃波响应,再将阶跃波响应与发射电流强度对时间的二阶导数进行褶积,得到了层状介质的全波形航空瞬变电磁响应。这一做法对于实际发射电流的全波形计算存在困难,因为实际发射电流采样率有限,难以可靠计算电流强度对时间的二阶导数值。孙怀凤等[6]研究了三维介质的全波形航空瞬变电磁响应计算问题,将源的电流密度加入麦克斯韦方程中,采用有限差分法实现了全波形瞬变电磁三维数值模拟算法。时间离散采用显式差分格式,为了保证计算结果的准确性步长要足够小,导致计算时间相对较长。齐彦福等[7]应用矢量有限元和后推Euler方法研究了全波形航空瞬变电磁的三维正演。对比显式方法,隐式差分格式无条件稳定[8],但是步长的选取仍要谨慎。

　　为了提高三维正演效率,众多学者将Krylov子空间方法引入到电磁场的计算中。Druskin、Knizhnerman等[9-10]最先将有理多项式方法引入时间域电场的直接求解。B?rner等[11]应用位移逆Krylov子空间(Shift-and-Invert Krylov,简称SAI Krylov)。子空间方法得到频域解,然后通过频时转换得到时间域的解。Druskin等[12]将有理Krylov子空间投影应用于求解麦克斯韦方程的三维扩散问题,并优化了有理子空间的偏移参数,使计算结果不再受限于模型电导率的差异。B?rner等[13]通过有理Krylov子空间方法求得频域解,然后用快速Hankel变换得到时域解;并提出了一种新的替代优化方法来选择有理Krylov子空间的极点参数,该方法可以在一个确定的迭代次数内收敛,而不受网格尺寸和电导率结构的影响,但频时转换的精度受到频点个数影响。周建美等[14]提出了频率域单极点模型降阶方法,使用拟态有限体积法实现了多频三维电磁问题的快速计算。邱长凯[15]等提出了一种有理Krylov子空间的指数加权偏移参数优化方法,使得有理Arnoldi近似在瞬变衰减晚期具备更高的精度,结合矢量有限元离散得到了航空瞬变电磁的阶跃波响应。ZHOU 等[16]开发了一种新的时间积分方法,引入二阶φ函数,通过构造新的指数函数形式得到了全波形瞬变电磁响应的表达式,基于结构化网格有限体积离散和SAI Krylov子空间方法实现了全波形瞬变电磁三维数值模拟。

　　为了更好地模拟起伏地表和复杂地电模型,本文基于非结构网格和SAI Krylov子空间方法,系统研究和实现了全波形航空瞬变电磁法三维数值模拟算法,验证了本文算法的精度,设计地表起伏异常体模型,计算和分析了其航空瞬变电磁响应特征。

　　1 基于SAI Krylov子空间的ATEM全波形三维数值模拟

　　1.1 ATEM电磁场满足的微分方程和边界条件

　　在研究区域Ω内,忽略位移电流和磁导率差异的影响,ATEM电磁场满足以下微分方程[17]:

　　其中,e(r, t)和h(r, t)分别为位置r处t时刻的电场强度和磁场强度,σ 为电导率,js(r, t)为源电流密度,μ0为真空中的磁导率。

　　根据(1a)和(1b)式可以得到电场满足的扩散方程为:

　　在研究区域Ω的边界上,采用狄利克雷边界条件,取边界上的切向电场为零。

　　ATEM常用的激发波形有阶跃波、半正弦波、梯形波等。本文研究ATEM全波形数值模拟方法,全波形的定义见图1,包含供电时间段(on-time)和断电时间段(off-time)。发射线框内的电流I在on-time时段非零,而在off-time时段为零,toff为on-time和off-time的交点。

　　

　　图1???任意发射波形示意图

　　1.2 ATEM全波形三维矢量有限元数值模拟

　　采用矢量有限元法[18]求解微分方程(2)。采用非结构四面体对整个研究区域进行剖分,如图2所示。电场定义在四面体单元的棱边上(图3)。则第k个单元内的电场ek(r, t)使用矢量基函数计算得到,具体表达式为:

　　其中,ekj(t)为第k个单元第j条棱边上的切向电场强度,

　　(r)为单元矢量基函数,j=1,2,…,6。

　　

　　图2???研究区域非结构网格剖分示意图

　　

　　图3???四面体单元棱边电场分布图

　　求解此线性方程组就可得到异常场电势Ψ,将其由解析解计算的一次场电势相加得到所有节点的总场电势,即所有测点电势,从而求出所有测点处的视电阻率值。

　　在整个研究区域Ω采用伽辽金法对微分方程(2)进行离散化、计算单元矩阵、整体合成等处理(具体细节参考文献[19]),得到有限元控制方程:

　　其中,A为刚度矩阵,B为质量矩阵,J(t)为外加电流源项,e(t)为所有剖分单元棱边上待求解的电场组成的向量。

　　(5)式可以进一步表示为:

　　一阶微分方程(6)式的解可以表示成以下形式:

　　其中,ti=ti-1+Δt,ei为ti时刻的电场,ei-1为ti-1时刻的电场。

　　在off-time阶段,S(t)项为0,电场的迭代公式为:

　　off-time时间段任意时刻ti电场响应的求解可以简化为统一以断电时刻toff为起点进行计算,而不再基于上一个时刻

　　。设断电时刻toff对应的电场为eoff,(9)式可以表示为:

　　ei=exp(-(ti-toff)B-1A)eoff,ti>toff

　　1.3 基于电偶极子离散的场源处理方法

　　由于ATEM实际数据采集过程中,发射线圈多为多边形回线,且很难保证位置的稳定,所以为了更加精确地模拟场源的形状和位置,采用基于偶极子离散的场源处理方法[20],将发射线圈近似为若干段首尾相连的电偶极子。假设发射回线穿过某个四面体单元,把穿过四面体的一段源线看作一个电偶极子,l为电偶极子的长度,rs为电偶极子空间位置的中点,v为发射电流的方向,电偶极子存在于第k个单元,则第k个单元的外加电流源项单元矩阵的具体表达式为:

　　1.4 初值条件

　　其中:N表示观测数据总个数,V表示剖分单元的尺寸。在非结构化网格中,每个四面体单元的尺寸无法具体算出,所以采取V为模型总尺寸除去剖分单元总个数得出的单元平均体积这种形式。需要注意的是,与正则化因子(拉格朗日因子)的选择相比,正则化项对δ的变化不敏感[21,24],所以采用平均体积这种做法对反演影响很小。

　　ATEM使用磁性源。当发射波形为理论阶跃波时,断电前发射线圈平稳供电,研究区域只存在稳定磁场而电场为零;当发射波形为其它波形如半正弦波、三角波、梯形波的时候,电场的计算从未供电时刻开始,研究区域不存在电场。综上两种情况,电场的初值条件为:

　　1.5 SAI Krylov子空间方法

　　在电场迭代公式(8)和(10)中,涉及矩阵指数函数φ2(-ΔtB-1A)和exp(-ΔtB-1A)与向量乘法的计算问题。由于矩阵的规模较大,常规计算技术和硬件条件难以对其直接快速求解。所以引入SAI Krylov子空间这一投影算法对矩阵指数函数和向量的乘法进行降阶求解。统一矩阵指数函数与向量乘法的形式为y=f(M)r, M为大型稀疏矩阵,r为向量,f为任意函数。SAI Krylov子空间方法基本思想是将目标向量y

　　从公式(13)可以看出,问题由原来的求解几十上百万阶的矩阵指数函数和向量的乘法问题f(M)r转变成了求解f(ξ(H-1m，m-I))这一几十至几百阶的矩阵指数函数的问题,后者的求解将非常容易。

　　在(8)和(10)式的求解过程中,B-1和A的矩阵乘无法直接得到,所以要先定义正交矩阵I和质量矩阵B的内积和范数为[13]:

　　其中,yH表示为y的共轭转置。若M=B-1A,引入正交矩阵I和质量矩阵B的内积和范数后,(13)式将变为:

　　使用Ruhe[21]的Gram-Schmidt正交化算法来构造Krylov子空间Km,并且得到(15)式中子空间的基Vm和H-1m，m，生成过程如图4所示。

　　

　　图4???有理Arnoldi近似

　　设(8)式中的-ΔtB-1A(-B-1Aei-1+Si-1)+Si-Si-1=r1, 令(10)式中的eoff=r2, 将(15)式代入到(8)式和(10)式中,可得到off-time和on-time任意时刻电场求解的最终表达式为:

　　对于t≤toff时段ei的求解,每次的迭代因为不同的时间间隔Δt对应不同的偏移参数ξon,为了避免由于Δt变化导致的多次矩阵分解,所以在on-time整个阶段取相同的迭代步长Δt和相同的偏移参数ξon。需要注意的是on-time时段的电场为迭代求解,下一时刻的电场ei的求解依赖于上一时刻的电场ei-1以及两个时刻的源项Si-1,这也导致了(16)式中的r1在随着迭代而变化,所以每次迭代都需要重新生成一次Krylov子空间。若从0时刻到toff时刻按时间间隔Δt共取non段,则有toff=Δt×non,所以整个on-time阶段构建Krylov子空间共需要1次大型稀疏对称矩阵分解和mon次右端项回代。

　　对于t>toff时段ei的求解,ti为toff时刻后的任意时刻,因为r2等于eoff,所以r2为固定值,因此只需要生成一次Krylov子空间便可以完成任意时刻ti对应的电场ei的求解。需要注意的是,当发射波形为阶跃波时,on-time时段的电场为0,eoff为断电后瞬间单元棱边的电场强度,令toff=0,按照(15)式中t>toff部分的算法便可以求解阶跃波断电后任意时刻的电场。

　　1.6 ATEM接收点处感应电动势的计算

　　计算得到研究区域所有棱边上的电场值后,任意接收点处的感应电动势dbz/dt由接收点所在单元的棱边上的电场以及相应的矢量插值函数计算,表达式如(17)所示:

　　其中,robs为接收点位置,ei(t)和︿n(robs)(robs)分别为接收点所在单元第i条棱边的电场和矢量插值函数。

　　2 数值算例

　　2.1 精度验证

　　为了检验基于SAI Krylov子空间方法的全波形响应数值模拟算法的计算精度,分别计算了均匀半空间模型阶跃波、半正弦波、三角波和梯形波(图5)的全波形响应感应电动势dbz/dt并与解析解进行了对比。

　　

　　图5???发射波形示意图

　　均匀半空间模型(图6)的电阻率为10?Ω·m,空气的电阻率为1×106Ω·m,发射源为单匝的十二边形回线,边长为7.59?m,距离地面高度为30?m,接收点位于发射线圈中点。设置研究区域的外边界沿坐标轴到中心点的距离为50?km,模拟空间的尺寸为100?km×100?km×100?km。使用开源代码Tetgen[22]对研究区域进行非结构四面体剖分,由于场源附近的电场变化比较剧烈,对测点以及源附近的网格进行局部加密,考虑到地面和空气界面的电阻率差异比较大,所以发射线圈下面对应的300?m附近的空气-地面界面进行面加密,最终生成127322个单元,20540个节点,147993条棱边。

　　

　　图6???均匀半空间模型

　　2.1.1 阶跃波

　　t=4 ms时发射电流瞬间关断(图5),计算断电后(0.01 ms,10 ms)区间范围内的感应电动势dbz/dt,设置子空间的阶数moff=77,偏移参数ξoff为1.6x104[23]。

　　SAI Krylov数值模拟结果与ATEM全波形响应的均匀半空间解析解[6]的对比结果如图7所示,基于SAI Krylov子空间方法计算的感应电动势dbz/dt与解析解的计算结果误差在3%以内,验证了三维数值模拟算法的计算精度。

　　

　　图7???阶跃波SAI Krylov数值模拟结果和解析解的比较

　　(a)dbz/dt响应;(b)相对误差

　　2.1.2 半正弦波/三角波

　　计算半正弦波和三角波(图5)的ATEM响应。由于两种波形的on-time时长相同并且on-time期间发射电流的一阶导数总不为零,所以两种波形计算时采用相同的参数。将on-time时段共分为100段,即令non等于100,设置子空间的阶数mon=24,偏移参数ξon为1/0.026Δt[24],计算断电后秒区间范围(0.01 ms,10 ms)内的感应电动势dbz/dt, 设置断电后的子空间的阶数moff为108, 偏移参数ξon为3.8x104[16]。

　　由图8计算结果可见:对于半正弦波,除了感应电动势正负号变号时刻的时间采样点,SAI Krylov子空间算法与解析解的误差绝大部分都在3%以内。对于三角波,SAI Krylov子空间算法在感应电动势正负号变号的时间采样点以及toff时刻附近误差较大,off-time的响应误差在4%以内。

　　

　　图8???半正弦波、三角波SAI Krylov计算结果和解析解的比较

　　(a)半正弦波dbz/dt响应;(b)半正弦波相对误差;(c)三角波dbz/dt响应(d)三角波相对误差

　　2.1.3 梯形波

　　计算梯形波(图5)的ATEM响应。?on-time阶段共4?ms,上升沿和下降沿为0.2?ms,发射线圈峰值电流为1?A,平稳供电持续时间为3.6?ms,off-time共10?ms。由于上升沿和下降沿较陡,电流的一阶导数相对较大,而平稳供电期内,发射回线不产生变化电场,电流的一阶导数为0,导致在第一个上升沿结束后总场变化较剧烈,考虑到梯形波的复杂情况,将on-time时段分为200段处理,令non等于200,子空间的阶数mon=24,偏移参数ξon为1/0.026Δt,对于off-time时段,计算断电后秒区间范围(0.01?ms,10?ms)的感应电动势,偏移参数与2.1.2小节中相同。

　　SAI Krylov的计算结果与解析解的计算结果如图9所示,两者之间的相对误差除了个别跳变点外在5%以内,on-time阶段在上升沿、下降沿与平稳供电期的交点处由于感应电动势发生了变号所以误差较大,off-time的感应电动势误差整体低于5%。

　　

　　图9???梯形波SAI Krylov结果和解析解的比较

　　(a)dbz/dt响应;(b)相对误差

　　分析4种波形的off-time响应(图10)可以看出:当线框的峰值发射电流相同时,阶跃波的二次场响应幅值最高,梯形波次之,然后是半正弦波,三角波幅值最小,表明同等条件下阶跃波的勘探能力更强一些。

　　

　　图10???阶跃波、三角波、半正弦波、梯形波的off-time?dbz/dt响应

　　2.2 VTEM实际发射波形

　　为了验证算法对复杂波形的ATEM响应的计算能力,对VTEM实际发射电流的ATEM响应进行计算。建立一个立方异常体模型(图11),异常体的尺寸为100?m×100?m×200?m,电阻率为1?Ω·m,设置地下背景电阻率为100?Ω·m。发射线圈设置为边长=6?m的八边形回线,接收点位于发射线圈中点。网格边界的设置和加密策略与2.1小节的模型相同。共计生成160305个四面体单元,25749个节点,186234条棱边。SAI Krylov子空间算法的相关参数与上文的梯形波相同。

　　

　　图11???立方异常体模型

　　将VTEM实际发射电流图(12(a)的ATEM响应的计算结果与后推欧拉算法(FETD)计算的结果[7]对比(图12),可以看出,由于VTEM实际波形非常复杂,所以on-time阶段的ATEM响应的变化也非常复杂,但最终两种算法的计算结果之间拟合非常好,说明本文算法具备准确计算复杂波形ATEM响应的能力。

　　

　　图12???VTEM实际发射波形及其ATEM响应

　　(a)VTEM实际发射波形;(b)SAI Krylov子空间方法和FETD计算的dbz/dt响应

　　2.3 复杂模型的ATEM响应

　　为了研究复杂模型的ATEM响应特征,建立如图13所示的复杂模型。地面有一个高度为30?m、直径为150?m的隆起地形;空气的电阻率为106Ω·m,地下背景的电阻率为100?Ω·m;在地形正下方存在一个半径为50?m、顶部埋深为50?m、电阻率为1?Ω·m的球体。发射线框为单匝的边长为10?m的十二边形回线,距离地面高度为50?m,发射波形为半正弦波,峰值发射电流为1?A。接收点位于发射线框中点,在x=-400?m至x=400?m范围内取点距为10?m,共有81个测点。在源的附近以及地表进行加密,最终生成490693 个单元,569478条棱边。

　　

　　图13???隆起地形球体模型示意图

　　图14给出了复杂模型半正弦波断电后不同时刻的dbz/dt剖面曲线。在0.03?ms以前,异常电流集中在隆起的地形上,呈现负异常,负异常左右两个峰的坐标与模型中地形隆起的边缘一致,大约在x=75?m和x=-75?m处。在0.03?ms至1.9?ms,异常电流聚集在球体,呈现正异常,而后随着时间的推移,感应电动势异常逐渐消失。

　　

　　图14???半正弦波off-time阶段不同时刻的dbz/dt剖面曲线

　　为了更好地分析地形对dbz/dt响应的影响,建立了水平地表均匀半空间模型和水平地表球形异常体模型,这两种模型的背景电阻率、球体几何参数和电阻率值均与图13所示的复杂模型相同。三种模型断电后0.255?ms的感应电动势计算结果如图15所示,从图中可以看出,不加地形的球体模型的响应对比均匀半空间的响应中间有明显的正异常,两边远离异常体的位置感应电动势的幅值开始逐渐接近均匀半空间的感应电动势,直到观测点的坐标大于330?m,异常体对观测值不再有影响,dbz/dt的幅值与均匀半空间相同。加了隆起地形后,0.255?ms时间道的幅值又略微低于水平地表球体模型的响应,原因可能是地形隆起导致的负异常叠加到了低阻球体的正异常上。计算结果表明,地形的起伏会影响ATEM感应电动势的分布,在ATEM数据处理以及反演时不能忽略地形的影响。

　　

　　图15???均匀半空间、水平地表球体、隆起地形球体断电后0.255?ms的dbz/dt响应

　　3 结 论

　　本文采用非结构四面体网格对空间进行离散,基于矢量有限元和SAI Krylov子空间方法实现了全波形航空瞬变电磁三维数值模拟算法。

　　SAIKrylov子空间方法求解ATEM响应时所需要的矩阵分解次数以及时间步迭代次数非常少。尤其是在off-time时段,只需要1次矩阵分解和多次右端项回代便可以生成Krylov子空间然后得到off-time所有待求时间点的dbz/dt响应。

　　通过计算和对比分析带地形的和不带地形的低阻球体模型ATEM响应dbz/dt,地形的起伏会影响ATEM感应电动势的分布,在ATEM数据处理以及反演时不能忽略地形的影响。

　　所开发的三维数值模拟算法计算精度高,可以更好地模拟复杂模型的全波形航空瞬变电磁响应,为带地形的全波形航空瞬变电磁三维反演奠定了基础。

　　致 谢

　　感谢长安大学地质工程与测绘学院周建美副教授、齐彦福博士在本文算法完成过程中提供的帮助。

　　参考文献：

　　[1]殷长春, 任秀艳, 刘云鹤, 等.航空瞬变电磁法对地下典型目标体的探测能力研究[J]. 地球物理学报, 2015, 58(9):3370-3379.

　　[2]赵越, 李貅, 王祎鹏, 等.三维起伏地形条件下航空瞬变电磁响应特征研究[J]. 地球物理学报, 2017, 60(1):383-402.

　　[3]武欣, 薛国强, 方广有.中国直升机航空瞬变电磁探测技术进展[J]. 地球物理学进展, 2019, 34(4):1679-1686.

　　[4]齐彦福. 复杂介质中时间域航空电磁数据仿真技术研究[D]. 长春:吉林大学, 2017.

　　[5]殷长春, 黄威, 贲放.时间域航空电磁系统瞬变全时响应正演模拟[J]. 地球物理学报, 2013, 56(9):3153-3162.

　　[6]孙怀凤, 李貅, 李术才, 等.考虑关断时间的回线源激发 TEM 三维时域有限差分正演[J]. 地球物理学报, 2013, 56(3):1049-1064.

　　[7]齐彦福, 殷长春, 刘云鹤, 等.基于瞬时电流脉冲的三维时间域航空电磁全波形正演模拟[J]. 地球物理学报, 2017, 60(1):369-382.

　　[8]UM E S, HARRIS J M, ALUMBAUGH D L.3D time-domain simulation of electromagnetic diffusion phenomena: A finite-element electric-field approach[J]. Geophysics, 2010, 75(4):F115-F126.

　　[9]DRUSKIN V, KNIZHNERMAN L.On application of the Lanczos method to solution of some partial differential equations[J]. Journal of Computational and Applied Mathematics, 1994, 50(1-3):255-262.

　　[10]DRUSKIN V, KNIZHNERMAN L.Extended Krylov subspaces: approximation of the matrix square root and related functions[J]. SIAM Journal on Matrix Analysis and Applications, 1998, 19(3):755-771.

　　[11]RALPH-UWE B, ERNST O G, SPITZER K.Fast 3-D simulation of transient electromagnetic fields by model reduction in the frequency domain using Krylov subspace projection[J]. Geophysical Journal International, 2008, 173(3):766-780.

　　[12]DRUSKIN V, KNIZHNERMAN L, ZASLAVSKY M.Solution of Large Scale Evolutionary Problems Using Rational Krylov Subspaces with Optimized Shifts[J]. SIAM Journal on Scientific Computing, 2009, 31(5):3760-3780.

　　[13]B?RNER R U, ERNST O G, GüTTEL S.Three-dimensional transient electromagnetic modelling using rational Krylov methods[J]. Geophysical Journal International, 2015, 202(3):2025-2043.

　　[14]周建美, 刘文韬, 刘航, 等.多频可控源电磁法三维有理函数Krylov子空间模型降阶正演算法研究[J]. 地球物理学报, 2018, 61(6):2525-2536.

　　[15]邱长凯, 殷长春, 刘云鹤, 等.三维时间域航空电磁有理Krylov正演研究[J]. 地球物理学报, 2020, 63(2):715-725.

　　[16]ZHOU J, LIU W, LI X, et al.3-D Full-Time TEM modeling using shift-and-invert Krylov subspace method[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(10):7096-7104.

　　[17]WANG T, HOHMANN G W.A finite-difference, time-domain solution for three-dimensional electromagnetic modeling[J]. Geophysics, 1993, 58(1):1646.

　　[18]NEDELEC J C.Mixed finite elements in 3[J]. Numerische Mathematik, 1980, 35(3):315-341.

　　[19]王建国. 电磁场有限元方法[M]. 西安: 西安电子科技大学出版社, 1998:164-198.

　　[20]YIN C, QI Y, LIU Y, et al.3D time-domain airborne EM forward modeling with topography[J]. Journal of Applied Geophysics, 2016, 134:11-22.

　　[21]RUHE A.The rational Krylov algorithm for nonsymmetric eigenvalue problems III: Complex shifts for real matrices[J]. BIT Numerical Mathematics, 1994, 34(1): 165-176.

　　[22]SI H. TetGen, a quality tetrahedral mesh generator and 3D delaunay triangulator[J/OL]. http://tetgen.berlios.de.2007.

　　[23]刘文韬. 基于拟态有限体积与有理 Krylov 空间方法的瞬变电磁三维正演模拟研究[D]. 西安:长安大学, 2019.

　　[24]JASPER V, HOCHBRUCK M.Preconditioning Lanczos approximations to the matrix exponential[J]. SIAM Journal on Scientific Computing, 2006, 27(4):1438-1450.