考虑发射波形的瞬变电磁三维模型降阶快速正演算法

周建美，刘文韬，齐彦福，李貅*，戚志鹏，鲁凯亮,刘正银,邢涛

1 长安大学地质工程与测绘学院，西安 710054 2 长安大学地球物理场多参数综合模拟实验室(中国地球物理学会重点实验室)，西安 710054 3 山东省交通规划设计院，济南 250031 4 北京探创资源科技有限公司，北京 100071

0 引言

瞬变电磁法作为一种重要的电磁勘探方法(底青云等，2019;Xue et al.,2020)，广泛应用于金属矿产勘探(Di et al.,2018;底青云等，2020)、煤田灾害探测(王若等,2018;Jiang et al.，2019;Chang et al.，2019)、水文地球物理研究(Danielsen et al.，2003)和隧道超前地质预报(Xue et al.,2018)等领域.理论上瞬变电磁发射波形设计为一个阶跃波，接收不含一次场的地下感应涡流产生的纯二次场信号(Nabighian and Macnae，1991).由于仪器制造上的限制(Davis and Macnae，2008)或探测目标的需要(殷长春等，2015)，实际中采用多种发射波形，包括半正弦波、三角波、梯形波(Liu，1998)和其他复杂波形(Davis and Macnae，2008;Chen et al.，2015).

发射波形对瞬变电磁响应的影响得到了广泛的研究.Liu(1998)、陈曙东等(2012)利用自由空间中回线作为目标体，推导了不同波形的解析解，并研究了不同发射波形对机载瞬变电磁响应的影响.Chen等(2015)开发了一个多波形(MULTIPULSE)系统，该系统在半周期内间断发射一个高功率半正弦波和一个低功率梯形波，保证探测深度的同时提供较高的近表面分辨率.殷长春等(2015)以半正弦波和梯形波为例，系统研究了几种典型地下目标情况下机载TEM方法的测量能力，并提供了参数组合以实现最佳勘探效果.不同发射波形的瞬变电磁响应不同，因此在开发瞬变电磁三维正演算法时需要考虑发射波形的影响(Zeng et al.，2019).

计算全波形瞬变电磁正演响应的方法主要有以下几种：(1)褶积算法(殷长春等，2013)，对阶跃响应与发射波形的时间导数进行褶积，可以得到任意发射波形的全时瞬变电磁响应.该方法可以直接利用高效的阶跃响应正演算法，但该算法需要获得精确的电流二阶导数，实际中发射电流数据采样率有限，难以获得电流二阶导数的精确值，严重影响了褶积结果的计算精度(齐彦福等，2017)；(2)时域有限差分法(Commer and Newman,2004；孙怀凤等，2013)，直接对连续发射电流波形进行采样并将其包含到控制方程式中，从而得到全时响应.该方法不需要求解大规模线性方程组，内存需求相对较低，但为了保持其数值稳定性，只能采用非常小的时间步长，从而需要较长的计算时间(孙怀凤等，2013)；(3)后退Euler法(Um et al.，2010；齐彦福等，2017；Zeng et al.，2019)，对发射波形的处理方法与时域有限差分法类似，虽然该方法是无条件稳定的，但仍需要谨慎选取时间步长，以在计算效率和准确性之间取得有效的平衡(Li et al.,2018).

对于三维阶跃响应计算，基于Krylov子空间投影的模型降阶算法(Zaslavsky et al.，2011;Börner et al.，2015;Zhou et al.，2018；Qiu et al.，2019)比时域有限差分法和后退Euler法具有更高的计算效率(Zaslavsky et al.，2011;Börner et al.，2015).结合直接法求解器，一次时间域正演的模型降阶算法的计算量可以减少为一次系数矩阵分解和上百次矩阵回代(Zhou et al.,2018).最近Zhou等(2020)将模型降阶算法扩展到了全波形瞬变电磁正演响应的计算，通过指数求积法则(Hochbruck and Ostermann,2010)，将on-time响应和off-time响应分别表示为关于空间离散系数矩阵的phi函数和指数函数，采用基于位移逆Krylov子空间投影的模型降阶算法(Zhou et al.,2018)实现了高精度正演计算.但该算法在计算on-time时间内每个时刻的正演响应时都需要构建一个新的位移逆Krylov子空间；而为了保证on-time响应的正演精度，一般需要在on-time范围内计算上百个时刻的响应，导致该算法在on-time时间内的正演响应计算非常耗时.

由于仪器设备的限制，目前实际工作中一般只对off-time数据进行处理解释，因此对于实际需求而言，只需要精确模拟任意发射波形的off-time响应即可.Smith和Neil(2013)指出，on-time发射波形非精确采样，也可以得到精确的off-time响应.Zhou等(2020)同样发现，在一定范围内改变on-time发射波形采样时刻数量，对off-time响应的精度没有显著影响.由于模型降阶算法(Zhou et al.,2020)on-time响应的计算时间与on-time范围内的采样时刻数量近似成正比，因此减少on-time范围内的采样时刻数量，就能够减少正演计算时间.基于频谱分析理论，时间域信号是频率域信号的叠加，而频率域信号的幅值随时间快速衰减，因此off-time某一时刻的响应只和特定频率范围内的信号有关.基于这一原则，本文提出了一种定量化评价on-time采样波形精度的方法，能够在保证off-time响应精度的同时，显著的减少on-time采样时刻，最大程度的减少正演计算时间.本文首先简述全波形瞬变电磁模型降阶算法；然后给出定量化评价on-time采样波形精度的方法；最后通过典型发射波形的正演计算验证本文算法的精度和计算速度.

1 正演算法

1.1 控制方程离散

全波形瞬变电磁法对应的时间域Maxwell方程，忽略位移电流：

(1)

式中，E是电场强度矢量，B是磁感应强度矢量，t是时间，σ是电导率，μ是磁导率，s是源项，考虑如图1所示的任意发射波形的电流源，在on-time时间段，0

图1 任意发射波形示意图Fig.1 Schematic diagram of the arbitrary transmitting waveform

采用简单的自然边界条件(Haber,2014)，利用拟态有限体积方法(周建美等，2018)对方程(1)进行空间离散，消去电场，得到关于磁场的空间离散控制方程：

(2)

(3)

式中，tk=tk-1+hk，k=1,2,…,n,hk是时间步长,t0=0,b0=0.

对于on-time时间段，公式(3)中g不为零.利用指数求积法则(Hochbruck and Ostermann,2010)对公式(5)中的积分项进行离散化，并引入函数φj(z),j=0,1,2,…，满足φj(0)=1/j!且具有递归关系φ0(z)=exp(z)，φj+1(z)=[φj(z)-φj(0)]/z，则on-time磁场响应可以表示为(Zhou et al.,2020)：

bk=bk-1+hk(-Abk-1+gk-1)+hkφ2(-hkA)

×[-hkA(-Abk-1+gk-1)+gk-gk-1]，

0

(4)

对于off-time时间段，公式(3)中g为零.为了避免在求解每个时刻的响应时都重新构建子空间，统一以b(toff)为起点(Börner et al.,2015;Zhou et al.,2018)，则off-time时间段任意时刻磁场的求解不再依赖于上一个时刻的磁场，简化后的公式(3)可以转化为：

bk=e-(tk-toff)Ab(toff),tk≥toff，

(5)

综合公式(4)和(5)，全时域磁场响应可表示为：

(6)

其中r=φ2(-hkA)[-hkA(-Abk-1+gk-1)+gk-gk-1]，q=e-(tk-toff)Ab(toff).可以看出，公式(6)中大部分的计算量集中于r和q的求解.它们可以写成统一的矩阵与向量乘积的函数形式

y=f(A)v，

(7)

对于on-time 时段，

f(A)=φ2(-hkA),

v=-hkA(-Abk-1+gk-1)+gk-gk-1,

0

(8)

对于off-time时段，f(A)=e(toff-tk)A,v=b(toff),tk≥toff.

(9)

1.2 位移逆Krylov子空间算法

采用基于位移逆Krylov子空间投影的模型降阶算法 (van den Eshof and Hochbruck,2006;Botchev,2016;Zhou et al.,2018)求解公式(7)的近似解，主要计算量是构建合适的Krylov子空间.对于给定的位移量γ和子空间阶数m，采用Arnoldi 方法(Zhou et al.,2018)构建Krylov子空间Km((I+γA)-1,v).根据子空间的正交基Vm和上Hessenberg矩阵Hm+1,m得到公式(11)的模型降阶解为(Hochbruck and Ostermann,2010)：

(10)

其中，Hm,m是上Hessenberg矩阵Hm+1,m中除去第m+1列，e1=[1,0,…,0]T.公式(10)中矩阵的维度m远远小于原系统中系数矩阵A的维度，可以快速求解.注意到公式(7)中向量v在off-time时段是常数，而在on-time时段随着发射时间变化，因此on-time时段的计算比off-time时段的计算复杂.为了得到整体高效的算法，对这两个时段分别采用不同的子空间近似策略进行求解.

在off-time时段，把公式(10)代入公式(9)中，得到off-time时段的磁场响应：

(11)

在on-time时段，把公式(10)代入公式(8)中，而后代入公式(6)，可得on-time时段的磁场响应：

bk≈bk-1+hk(-Abk-1+gk-1)+hkVmkφ2

(12)

1.3 频谱分析

在on-time时段进行密集时刻采样，可以得到精确的全波形正演响应(Zhou et al.,2020).但在on-time时段每个计算时刻tk需要重新构建子空间Kmk((I+γkA)-1,vk)，导致on-time响应的计算效率低下.如果只关注off-time响应，则on-time采样时刻可以显著减少，从而减小总的正演计算时间.

由傅里叶变换理论，任意一个时间域信号，都可以分解为无穷多个不同频率的正弦信号的叠加(胡广书,2003)，因此可以从频率域信号出发，分析接收的时间域信号的特点.如图1所示，瞬变电磁发射波形是一个连续非周期信号.实际工作中对发射波形进行离散化采样，得到时间域中一个离散非周期信号.对该信号进行快速傅里叶变换,即可得到该发射波形对应的频谱分布(胡广书,2003).下面从频谱分析角度，根据一维介质中电磁场的分布特征，定量化的给出on-time时段的波形采样规则.

(1)根据电磁场的传播理论，忽略位移电流，一维介质中频率域磁场的通解为(汤井田和何继善，2005)：

B=B0e-z/δeiφ，

(13)

Vp=ωδ.

(14)

对于瞬变电磁法，在发射源关断后t时刻接收到的时间域信号中，频率为f的场的相位传播距离为

z1=Vpt=2πfδt.

(15)

(2)常见发射波形的频谱中，高频信号的幅值相比低频信号的幅值小(牛之琏，2007)，对于测量的时间域信号的贡献有限，可以进一步将比较波形的频谱范围限定为大于最大幅值Amax的1%的频率范围.本文令Ac=0.01Amax为截止振幅.

基于上述两条原则，可以在保证off-time响应计算精度的同时，极大的减少on-time采样时间的数量.

2 数值实例

通过对比基于位移逆Krylov子空间模型降阶的全波形算法(Zhou et al.,2020)(简称为SAI)和本文结合频谱分析的位移逆Krylov子空间快速算法(简称为SAI-SA)计算典型的半正弦波、梯形波和VTEM波的正演响应，测试本文算法在计算off-time响应的精度和速度.SAI算法对on-time波形进行密集时刻采样(Zhou et al.,2020)，SAI-SA算法对on-time波形进行基于频谱分析的稀疏时刻采样.本文计算设备为32 G内存、四核主频3.6 GHz的Intel i7 CPU的台式电脑.

2.1 半正弦波

首先采用均匀半空间模型，通过与一维解析解结果(殷长春等，2013)对比，验证本文算法的计算精度.半空间模型的电阻率为10 Ωm，空气层电阻率设置为10-6Ωm，采用中心回线装置的航空装置，发射线圈是边长6.66 m的正八边形，飞行高度为30 m.发射波形设置为单位电流强度、脉冲宽度为4 ms的半正弦波(图2a中黑色曲线所示)，接收off-time时间段[0.01,10]ms的dbz/dt响应.

图2 半正弦波形和采样数nton=6的采样波形对比(a)半正弦波形(黑色)和采样波形(灰色)；(b)半正弦波形(黑粗线)和采样波形(灰粗线)频谱分布，虚线对应截止频率，黑色细实线对应截止振幅；(c)采样波形的频谱分布与半正弦波形的频谱分布的比值.Fig.2 Comparison of half-sine waveform and sampling waveform with nton=6(a)Half-sine waveform (black)and sampling waveform (gray);(b)Spectrum of half-sine waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and half-sine waveform.

图3 半正弦波形和采样数nton=7的采样波形对比(a)半正弦波形(黑色)和采样波形(灰色)；(b)半正弦波形(黑粗线)和采样波形(灰粗线)频谱分布，虚线对应截止频率，黑色细实线对应截止振幅；(c)采样波形的频谱分布与半正弦波形的频谱分布的比值.Fig.3 Comparison of half-sine waveform and sampling waveform with nton=7(a)Half-sine waveform (black)and sampling waveform (gray);(b)Spectrum of half-sine waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and half-sine waveform.

图4 半正弦波形和采样数nton=8的采样波形对比(a)半正弦波形(黑色)和采样波形(灰色)；(b)半正弦波形(黑粗线)和采样波形(灰粗线)频谱分布，虚线对应截止频率，黑色细实线对应截止振幅；(c)采样波形的频谱分布与半正弦波形的频谱分布的比值.Fig.4 Comparison of half-sine waveform and sampling waveform with nton=8(a)Half-sine waveform (black)and sampling waveform (gray);(b)Spectrum of half-sine waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and half-sine waveform.

接下来采用采样数nton=7的SAI-SA算法计算半正弦波的正演响应，同时采用密集采样数nton=100的SAI算法计算的结果(Zhou et al.,2020)作为三维高精度结果进行比较.结果如图5所示，图5a为三维正演结果，图5b为三维正演结果相对一维解析解的误差，其中实线为一维解析解，圆圈为采样数nton=100的SAI算法的计算结果，叉点为采样数nton=7的SAI-SA算法的计算结果.由图5可知，密集采样的SAI算法能够得到on-time和off-time时间段内的高精度正演结果；而稀疏采样的SAI-SA算法的正演结果虽然在on-time时间段的误差很大，但在off-time时间段内与密集采样的SAI正演结果基本一致，相对一维解析解的相对误差均小于5%.

图5 半正弦波模型的dbz/dt响应(a)及误差对比(b)实线为一维解析解，圆圈为采样nton=100的SAI算法的正演结果，叉点为采样nton=7的SAI-SA算法的正演结果.Fig.5 dbz/dt response (a)and error comparison (b)of the half-sine waveform modelThe line is 1D analytical solution,circle is result of SAI algorithm with nton=100 and cross is result of SAI-SA algorithm with nton=7.

表1给出了密集采样的SAI算法和稀疏采样的SAI-SA算法的计算时间对比.由表1可知，采用密集采样nton=100的SAI算法，由于每次计算on-time响应时需要构建一次新的Krylov子空间，导致总的系数矩阵回代次数高达1622，总的正演时间为217.1 s；而采用稀疏采样nton=7的SAI-SA算法，计算on-time响应只需要构建7次新的Krylov子空间，从而总的系数矩阵回代次数减少为171，总的正演时间减少为39.5 s，计算速度提高为密集采样算法的5.5倍.

表1 SAI-SA算法和SAI算法计算半正弦波模型的时间对比Table 1 Solution time of SAI-SA algorithm and SAI algorithm for a half-sine waveform model

2.2 梯形波

接下来考察梯形波的正演.采用与半正弦发射波形相同的半空间模型和发射接收装置系统.发射波形设置为单位电流强度，脉冲宽度为4 ms、上升沿和下降沿均为0.2 ms的梯形波(图6a中黑色曲线)，接收off-time时间段[0.01,10]ms的dbz/dt响应.利用SAI-SA算法计算off-time响应.采用与半正弦发射波形相同的空间离散网格.先确定on-time时间段的采样波形.对发射波形采用线性等间隔采样，显然当采样数nton=20(不包括初始时刻t=0)时，采样波形与原梯形波完全重合(图6a中灰色圆圈)，频谱分布也完全一致(图6b).注意到梯形波的规则形状，如果采用线性不等间隔采样，采样数nton=3(不包括初始时刻t=0)，既可以保证采样波形与原梯形波完全重合(图7a)，频谱分布也完全一致(图7b).

图6 梯形波形和采样数nton=20的采样波形对比(a)梯形波形(黑色)和采样波形(灰色)；(b)梯形波形(黑粗线)和采样波形(灰粗线)频谱分布，虚线对应截止频率，黑色细实线对应截止振幅；(c)采样波形的频谱分布与梯形波形的频谱分布的比值.Fig.6 Comparison of trapezoid waveform and sampling waveform with nton=20(a)Trapezoid waveform (black)and sampling waveform (gray);(b)Spectrum of trapezoid waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and trapezoid waveform.

图7 梯形波形和采样数nton=3的采样波形对比(a)梯形波形(黑色)和采样波形(灰色)；(b)梯形波形(黑粗线)和采样波形(灰粗线)频谱分布，虚线对应截止频率，黑色细实线对应截止振幅；(c)采样波形的频谱分布与梯形波形的频谱分布的比值.Fig.7 Comparison of trapezoid waveform and sampling waveform with nton=3 (a)Trapezoid waveform (black)and sampling waveform (gray);(b)Spectrum of trapezoid waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and trapezoid waveform.

但在on-time时间段采用线性不等间隔采样时，需要重新考虑最优化位移量γon的选取.当on-time时间段采用线性等间隔采样，SAI算法的最优化位移量γon与采样间隔成正比(Zhou et al.,2020).当在on-time时间段采用线性不等间隔采样时，每个采样间隔对应一个最优化位移量，但采用多个不同位移量会增加系数矩阵的LU分解次数，从而增大计算时间.由于位移量对采样间隔不灵敏(van den Eshof and Hochbruck,2006；Zhou et al.，2018)，位移量γon选取可以采用两种简单的方法：(1)选取某个采样间隔对应的最优化位移量作为总的最优化位移量，比如以第一个时间间隔对应的最优化位移量作为总的最优化位移量，对于采样数nton=3，采用SAI-SA算法，求解on-time响应需要构建的3次Krylov子空间，子空间阶数分别为m1=19，m2=55，m3=19；(2)选取off-time时间段对应的最优化位移量γoff作为总的最优化位移量，对于采样数nton=3，采用SAI-SA算法，求解on-time响应构建的3次Krylov子空间的阶数分别为m1=16，m2=30，m3=16.两种方法需要构建的Krylov子空间阶数相差不大，对应计算量也相差不大，而且都表现为大的采样间隔对应大的子空间阶数.但选取γoff作为总的最优化位移量时，求解on-time时间段正演响应需要LU分解的系数矩阵与求解off-time时间段正演响应需要LU分解的系数矩阵相同，因此可以在求解过程中减少1次系数矩阵的LU分解，表现出显著的计算优势.因此本文选取γon=γoff来计算线性不等间距稀疏采样波形的正演响应.

接收off-time时间段[0.01,10]ms仍然离散为31个对数等间隔时刻.分别采用密集采样数nton=100的SAI算法(Zhou et al.,2020)、采样数nton=20的线性等间隔采样的SAI-SA算法、采样数nton=3的线性不等间隔采样的SAI-SA算法和一维解析算法(殷长春等，2013)计算梯形波模型的正演响应.结果如图8所示，图8a为正演计算结果，图8b为三维正演结果相对一维解析解的误差，其中实线为一维解析解，圆圈为nton=100的SAI算法的计算结果，三角形为nton=20的SAI-SA算法的计算结果，叉点为nton=3的SAI-SA算法的计算结果.由图8可知，nton=100的密集采样的SAI算法能够得到on-time和off-time时间段内的高精度正演结果；而nton=20和nton=3的稀疏采样的SAI-SA算法虽然在on-time时间段的误差很大，但由于nton=20和nton=3的稀疏采样的频谱分布与nton=100的密集采样的频谱分布是相同的，都与原梯形波完全重合，因此三种采样得到的off-time时间段内的正演响应是完全相同的，相对一维解析解的相对误差均小于5%.

图8 梯形波模型的dbz/dt响应(a)及误差对比(b)实线为一维解析解，圆圈为采样nton=100的SAI算法的正演结果，叉点为采样nton=3的SAI-SA算法的正演结果，三角形为采样nton=20的SAI-SA算法的正演结果.Fig.8 dbz/dt response (a)and error comparison (b)of the trapezoid waveform modelThe line is 1D analytical solution,circle is result of SAI algorithm with nton=100,cross is result of SAI-SA algorithm with nton=3,and triangle is result of SAI-SA algorithm with nton=20.

表2给出了密集采样的SAI算法和稀疏采样的SAI-SA算法的计算时间对比.由表2可知，采用密集采样nton=100的SAI算法，由于总的系数矩阵回代次数高达1546，总的正演时间为213 s；采用线性等间隔稀疏采样nton=20的SAI-SA算法，总的系数矩阵回代次数减少为303，总的正演时间减少为56.5 s，计算速度提高为密集采样算法的3.8倍；采用线性不等间隔稀疏采样nton=3的SAI-SA算法，不仅总的系数矩阵回代次数减少为171，而且减少了一次LU分解，因此总的正演时间减少为36.5 s，计算速度提高为密集采样算法的5.8倍.

表2 SAI-SA算法和SAI算法计算梯形波模型的时间对比Table 2 Solution time of SAI-SA algorithm and SAI algorithm for a trapezoid waveform model

2.3 VTEM波

进一步采用复杂VTEM波的三维模型(齐彦福等，2017)的正演响应验证本文算法的有效性.三维模型如图9所示，电阻率为100 Ωm的均匀半空间中存在一个电阻率为1 Ωm的低阻异常体，异常体上顶面距离地表50 m，异常体尺寸为100 m×100 m×200 m，空气层电阻率设置为106Ωm，在三维目标体的正上方采用与半正弦发射波形相同的发射接收装置系统.

图9 三维模型Fig.9 3-D model

图10 VTEM波形和采样数nton=26的采样波形对比(a)VTEM波形(黑色)和采样波形(灰色)；(b)VTEM波形(黑粗线)和采样波形(灰粗线)频谱分布，虚线对应截止频率，黑色细实线对应截止振幅；(c)采样波形的频谱分布与VTEM波形的频谱分布的比值.Fig.10 Comparison of VTEM waveform and sampling waveform with nton=26 (a)VTEM waveform (black)and sampling waveform (gray);(b)Spectrum of VTEM waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and VTEM waveform.

接下来采用采样数nton=26的SAI-SA算法计算VTEM发射波形的正演响应，同时采用密集采样数nton=200的SAI算法计算的结果(Zhou et al.,2020)作为三维高精度结果进行比较.结果如图11所示，图11a为三维正演结果，其中实线为三维隐式时间步迭代算法的计算结果(齐彦福等，2017)，圆圈为采样数nton=200的SAI算法的计算结果，叉点为采样数nton=26的SAI-SA算法的计算结果.图11b为nton=26的SAI-SA算法的计算结果相对nton=200的SAI算法的计算结果的误差.由图11可知，与半正弦波模型正演和梯形波模型正演结果一致，密集采样的SAI算法能够得到on-time和off-time时间段内的高精度正演结果；而稀疏采样的SAI-SA算法的正演结果虽然在on-time时间段的误差很大，但在off-time时间段内与密集采样的SAI正演结果基本一致，相对密集采样的正演结果的相对误差小于5%.

图11 VTEM波模型的dbz/dt响应(a)实线为三维BDF2算法的正演结果，圆圈为采样nton=200的SAI算法的正演结果，叉点为采样nton=26的SAI-SA算法的正演结果;(b)为nton=26的SAI-SA算法的正演响应相对nton=200的SAI算法的正演响应的误差.Fig.11 dbz/dt response of the VTEM waveform model(a)The line is result of 3D BDF2 algorithm,circle is result of SAI algorithm with nton=200,cross is result of SAI-SA algorithm with nton=26;(b)Relative error of response of SAI-SA algorithm with nton=26 to response of SAI algorithm with nton=200.

表3 SAI-SA算法和SAI算法计算VTEM波模型的时间对比
Table 3 Solution time of SAI-SA algorithm and SAI algorithm for VTEM waveform model

表3给出了密集采样的SAI算法和稀疏采样的SAI-SA算法的计算时间对比.由表3可知，采用密集采样nton=200的SAI算法，由于总的系数矩阵回代次数高达3114，总的正演时间为903.9 s；采用线性不等间隔稀疏采样nton=26的SAI-SA算法，不仅总的系数矩阵回代次数减少为436，而且减少了一次LU分解，因此总的正演时间减少为154.4 s，计算速度提高为密集采样算法的5.9倍，在保证off-time响应计算精度的同时，表现出非常好的加速效果.

2.4 不同发射波形的探测能力分析

为了比较不同发射波形的勘探能力，参考Commer和Newman(2004)，设计如图12的三维垂直接触带模型.空气电阻率为106Ωm，地表下方是厚度为50 m、电阻率为100 Ωm的覆盖层，覆盖层下方由两部分不同电阻率的垂直接触带构成，电阻率分别为500 Ωm和800 Ωm.垂直接触带中间存在一个电阻率为2 Ωm的三维复杂形状的低阻体，沿x方向的长度为400 m，沿y方向的宽度为100 m，沿z方向的厚度为500 m，其具体位置信息如图12所示.采用中心回线的航空装置，发射线圈是边长6.66 m的正八边形，飞行高度为30 m.采用直立六面体网格对计算区域进行剖分，最小网格为10 m，网格扩大因子为1.36，三个方向的网格数为66×76×60.

SAI算法(nton=200)SAI-SA算法(nton=26)求解次数求解时间(s)求解次数求解时间(s)矩阵分解228.11 13.9矩阵回代3114833.6436116.0搜索位移量γ7.57.5求解方程(7)和(8)23110.9574.8其他23.812.2总计903.9154.4

图12 三维垂直接触带模型Fig.12 3D model of a conductor at a vertical contact

为了便于比较，采用具有相同的发射能量(发射电流与时间轴围成的面积相同)的半正弦波和梯形波形，其中半正弦波的持续时间为4 ms，发射电流为1 A；梯形波的上升沿和下降沿均为0.2 ms，持续时间为4 ms，发射电流为0.58 A.图13为半正弦波和梯形波的off-time时段[0.01,50]ms时间范围内∂bz/∂t响应的多测道图.从图13可以看出：半正弦波和梯形波的off-time时段∂bz/∂t响应均能对地下低阻异常体有明显的反映，在x=-25 m区域，对应模型图12中低阻异常体最浅部位置，多测道曲线呈下凹特征，显示明显的低阻异常.

图13 不同发射波形off-time阶段dbz/dt响应的多测道图(a)半正弦波；(b)梯形波.Fig.13 Profiles of off-time dbz/dt responses for different transmitting waveforms(a)Half-sine wave;(b)Trapezoid wave.

图14为x=-25 m处正上方半正弦波与梯形波在off-time阶段的∂bz/∂t响应曲线，其中虚线为电阻率100 Ωm的半空间的正演响应曲线，实线为图12所示三维模型的正演响应曲线，灰色为半正弦波的响应曲线，黑色为梯形波的响应曲线.由图14可知，半正弦波和梯形波的响应曲线对异常体均有明显的反映：早期响应相比半空间响应的数值大，显示为低阻异常；晚期响应相比半空间响应的数值小，则表现为高阻地层的影响.对比图14中的曲线幅值可知，在激发能量相同时，早期的梯形波响应相比半正弦波响应的幅值大，因此在相同噪声条件下，梯形波可以得到信噪比更好的早期响应；但晚期的梯形波响应曲线和半正弦波响应曲线基本重合，说明发射波形对晚期响应影响小.

图14 x=-25 m处不同发射波形off-time阶段的dbz/dt响应对比灰色为半正弦波响应，黑色为梯形波响应.Fig.14 Comparison of off-time dbz/dt responses at x=-25 m for different transmitting waveformsGray for half-sine wave and black for trapezoid wave.

3 结论

本文采用位移逆Krylov子空间模型降阶算法结合频谱分析理论实现了考虑发射波形的瞬变电磁off-time响应的三维快速正演.利用指数求积法则，将on-time响应和off-time响应表示为统一的矩阵与向量乘积的函数形式，采用位移逆Krylov子空间模型降阶算法实现有效求解.为了减少计算时间，从频谱分析角度，定量化的给出了on-time时段的波形采样规则，减少了on-time波形采样数，从而在保证off-time响应计算精度的同时，极大的减少了总的正演计算时间.

半正弦波模型、梯形波模型以及复杂的VTEM波模型的正演计算结果表明，对于规则变化的发射波形，比如半正弦波，采用线性等间隔采样，能够实现好的加速效果；对于具有典型拐点的发射波形，比如梯形波和VTEM波，采用线性不等间隔采样，能够实现更好的加速效果.三种发射波形的模型正演结果，相比密集采样的常规正演结果(Zhou et al.,2020)，在保证off-time响应计算精度的同时，都实现了5倍以上的加速，表现出非常稳定的加速效果.

线性不等间隔采样能够最大化的实现加速效果，但对于复杂的发射波形，比如VTEM，需要进行手动选取采样点，很难快速直接的得到最优化的线性不等间隔采样波形.未来进一步工作考虑引入有效的自适应选取线性不等间隔采样波形的方法.