0 引言

多通道瞬变电磁法是由英国爱丁堡大学提出的(Wright et al., 2001)，它一次布设多道采集电极，通过沿测线逐次移动发射源的方式达到测深的目的.该方法在陆上、海上油田的油气水识别中取得了成功，申请了专利，受到了国内外关注(Ziolkowski et al, , 2007; Anderson et al, , 2008;汤井田和罗维斌，2008；武欣等，2015)，国家重大科研装备项目的相关科研人员对多通道瞬变电磁法从仪器到方法进行了全面研究.

多通道瞬变电磁法的工作方式与传统瞬变电磁法不同.传统瞬变电磁法多采用线框做为发射源(牛之琏，2007)，发射有一定占空比的方波，接收三个方向的二次磁场，因此信号比较弱，易受天然、人文电磁信号的干扰.多通道瞬变电磁法中常采用伪随机系列编码源连续激发，同时在近源处监测发射信号(认为发射信号已知)，在资料处理中，用发射信号与多通道接收设备记录的时间域信号进行相关辨识，获得接收点处的大地脉冲响应.

综上所述，多通道瞬变电磁法有两个特色：一是装置形式与多维反射地震勘探和高密度电法类似，一次性布设接收电极，通过沿测线移动发射源实现测深，记录点为发射源与接收点的中点；二是发射源为伪随机不间断信号源，接收一次电场与二次电场的叠加信号，因此接收信号的强度明显大于传统瞬变电磁法的二次场信号.

国内外学者对多通道瞬变电磁法研究较多，国外的文献多集中在方法机理与应用研究方面(Wright et al., 2002；Ziolkowski et al., 2007, 2010).国内学者多把研究重点放在了相关辨识及对噪声的压制上，对辨识的方法理论做了充分的探讨，汤井田和罗维斌(2008)分析了逆重复m序列的波形特征、相关函数及频谱特性，讨论了对各种干扰的压制方法及测量精度，薛国强等(2015)分析了MTEM伪随机编码发射的响应特性和方法的关键技术，武欣等(2015)研究了m序列伪随机编码源电磁响应的精细辨识问题.在电磁场偏移成像技术上，国内学者也做了尝试，智庆全(2015)对MTEM波场变换与偏移成像方法进行了研究，钟华森等(2016)对多道瞬变电磁法虚拟波场提取技术做了专门研究.在正演模拟方面，齐彦福等(2015)利用方波响应移位叠加和电流导数与阶跃响应褶积两种方法实现理论m序列和实际发射波形的全时正演模拟，我们也进行了伪随机编码源的信号合成研究(王若等, 2016)，但这些都只局限于一维模拟，对于较复杂的模型，还需要做多维正演模拟研究.本文试图开展多通道瞬变电磁法的2D有限元大地脉冲响应的正演模拟研究.

本文的目的是采用2D频率域有限元方法和程序获得2D时间域大地脉冲响应，也就是已知电性结构模型，数值正演计算理论大地脉冲响应，并研究如何用大地脉冲响应有效地表征目标体的存在，以此来探讨多通道瞬变电磁法的应用效果.

1 2D有限元正演模拟方法的实现 1.1 二维电场的边值问题

若电磁场的时间因子为e^－iωt，则含源二维电场的边值问题可写为

(1)

其中，，q=σ－iωε，ε和μ分别为介质中的介电常数和磁导率，在本文中，ε和μ分别取为空气中介电常数和磁导率ε₀和μ₀，ω为角频率，σ为介质的电导率.u为x方向的电场E_x，g=－J_cx，为x方向的发射电流密度，k₀和k分别为空气中的波数和地层中的波数，Γ₁为内部边界，DABC为空中边界，DEFC为地下外边界，空中边界和地中边界共同构成计算区域Ω.

1.2 有限元方法

应用变分理论，与边值问题(1)式等价的变分为

(2)

将计算区域剖分成系列矩形单元，用双线性插值函数作为形函数(i=1, 2, 3, 4)，每个单元的场值u可用四个角点的场值u_i与形函数N_i表示为

(3)

将(3)式代入(2)式，经推导可得2D频率域有限元方程：

(4)

其中，K为有限元的系数矩阵，U为电场矢量，P为源项，具体表达式可参考相关文献(王若等，2006).用LU分解法解方程组，即可得到各节点的电场值.通过公式可得到各个节点上的磁场.

1.3 宽频带2D有限元结果的修正

我们已用2D有限元方法实现了可控源音频大地电磁法(CSAMT)的正演模拟，并取得了精度较高的结果(王若等，2006).相较于CSAMT方法的频带范围(0.25~8192 Hz)，多通道瞬变电磁法中的2D有限元方法需要的频带范围要宽得多，带宽从10^-6 Hz到10⁶Hz，因此，把原程序代码移植于多通道瞬变电磁法前，需要对其带宽范围内的频率域计算精度进行二次检验.这里针对100 Ωm的均匀半空间模型，分别用2D有限元方法及利用源叠加的解析公式(朴化荣，1990；方文藻等，1993)计算了地面电磁场，二者对比结果如图 1所示.

图 1a中的电场频谱图有三条曲线，实心圆点为2D有限元计算的宽频带谱图，空心圆圈为用相同程序计算的CSAMT频谱范围的谱图，黑色实线为用解析公式计算的宽频带谱图，图 1b为磁场谱图.从图中可以看出，在比CSAMT稍宽一些的频带范围内，2D有限元的计算结果与解析结果基本吻合，但对电场的频谱来说，在宽频带的高频(>8192 Hz)与低频(< 0.003 Hz)部分，2D有限元的计算结果显著偏离了解析解.磁场在低频部分，与解析解偏离时的频点稍高一些.

通过分析，发现造成2D有限元正演程序在宽频带的高频和低频与解析解偏离的原因互不相同.高频时2D有限元方法的误差应是数值频散造成的，在网格剖分固定时，频率越高，频散越严重.低频时的误差应是边界影响引起的，在低频时，电磁波波长较长，传播较远，衰减较慢，即便将边界设置得离源很远，随着传播时间的加长，边界上的电磁场值将不趋近于零，边界条件不再适用，因此导致低频时2D电磁场结果偏离解析解，频率越低，偏离程度越高.

根据上面的分析，可以通过加密网格降低数值频散影响，并可将边界设置得更远以保证低频场在边界上的值趋近于零.但是这种解决方案无疑会大大增加网格的剖分数量，对于有限元方法来说，即便是在2D条件下实施，也会引起内存及计算时间的大量消耗，因此，需寻找更好的解决方法.根据电磁波的传播理论，在本文所说的高频时，电磁波的能量主要集中于浅层，甚至是表层中，满足远区条件，因此可以用远区的近似公式来代替有限元法获得高频部分的解答；同样，对于低频，电磁波的传播满足近区传播规律，用近区公式来代替有限元法计算低频的解，便可得到整个宽频带内的电磁场解.采用这种方法对2D有限元法进行修正，可以不显著增加有限元的计算时间及内存，便可得到较正确的解.

本文采用第二种方法对2D有限元结果进行修正.在近区，电磁场修正公式为

(5)

(6)

在远区，电磁场修正公式为

(7)

(8)

其中，P_E为电偶极距，r为收发距，θ为源与发收距矢径的夹角，ρ₁为地下介质的电阻率.

修正后的结果如图 2所示.由图 2可见，通过校正后的2D有限元电磁场计算结果与解析解吻合程度较高，说明所用的校正方法可行.

1.4 远区与近区阈值的讨论

本文涉及到远区与近区的校正，但不同模型的远、近区频率范围会不同，而如何界定远区与近区的频率范围，是一个难题.在频率域有源电磁法中，远区与近区的划分通常用|k₁r|的值来判断(朴化荣, 1990).理论上说，在均匀半空间中，当为均匀半空间的波数，r为收发距)时，定义为近区，当时，定义为远区.

根据远区与近区的定义，当地下介质一定时，其划分不只与频率有关，还与收发距有关.本文的修正包括高频(远区)资料修正与低频(近区)资料的修正，对于高频电磁场资料，如前所述，计算误差由数值频散引起，因此，用频率作为判断标准比用|k₁r|更为合理.本文将频率大于等于10000 Hz作为修正范围，采用这个频率阈值，针对文中的模型，当r>200 m时，都满足的条件，r越大，远区条件越会得到满足，因此，从远区意义上说，频率阈值选为10000 Hz也是合理的；对于低频(近区)电磁场资料的修正，先用成熟软件(针对不同背景，本文用自编源叠加的均匀半空间电磁场响应程序或Utah积分方程软件)计算出某一收发距处的背景场，与用有限元计算的背景场进行比较，找到二者在低频部分产生分歧的频率点，用此频率点计算出相应收发距处的|k₁r|值，把该值做为低频阈值，然后用近场公式对低频场资料进行修正.经过计算，并通过和已有程序计算结果比对，表明用这种方法得到的阈值可以应用到该模型所有收发距的低频校正中.

注意计算层状背景的|k₁r|及应用(5)式时，电导率为最后一层的电导率，这是考虑到低频时，电磁波穿透较深，且具有体积效应，针对本文模型(见3.1节中图 4)，低频时对电磁波有主要贡献的是最后一层介质.

2 瞬变电磁脉冲响应 2.1 频时变换

频时变换方法有多种，如正弦变换、余弦变换、G-S变换、Euler算法、Talbot算法等，前面两种属于逆傅氏变换的变种，后面几种属于逆拉氏变换(Li et al., 2016).正弦变换和余弦变换可以通过折线法和数字滤波法(李锋平等, 2016; Anderson et al., 2008; Kauffman, 2014; 朴化荣, 1990; 李大俊等, 2014)来实现.通过研究发现：折线法是近似方法，其计算精度受区间分割数的影响较大(李大俊等, 2014)，对于复杂介质可能会产生误差；G-S变换虽然计算速度快，但需要4倍精度来实现，对计算机的要求太高；而余弦变换不需要对变量定义太高的精度，而且计算精度也较高(昌彦君等，2000).因此，本文选择余弦变换作为频时变换的手段.

余弦变换的表达式为

(9)

其中，t为时间，ω为角频率，F(ω)为与频率有关的核函数.选用合适的函数对求取余弦变换的滤波权系数，利用滤波系数可将(9)写为

(10)

其中，C_i为汉克尔变换的滤波权系数，x为时间的自然对数，A_i=A₀+ilne，A₀=－30.3025124，e= 1.221402758160169834，i为1~281之间的正整数，N₁和N₂是与被积函数的性质和滤波截止误差有关的正整数.

本文利用了方文藻等(1993)的滤波权系数，并参考公开的余弦变换子程序，根据公式(10)将程序适当修改，使之成为适合于2D正演计算调用的形式.

2.2 大地脉冲响应

频率域的计算结果经余弦变换得到的是阶跃响应，采用余弦变换的另一个原因是因为阶跃响应比较容易验证，获得阶跃响应后再通过对阶跃响应求取时间的导数即可得到时间域的脉冲响应.

图 2中的频率域结果所对应的阶跃响应与大地脉冲响应如图 3所示.为了检验频时变换的准确性，将源叠加频率域解析解经余弦变换得到的阶跃响应与脉冲响应数值解(下面称为1D数值解)也示于图中，通过二者对比发现，2D有限元结果经余弦变换后，得到的阶跃响应与1D数值解结果基本吻合，电场的吻合度要高于磁场的吻合度，经过对时间求导后，磁场与1D数值解的差异进一步缩小.为了突出电场的脉冲响应，将图 3c的纵轴用算术坐标表示，从图 3c看出，由2D有限元得到的计算结果早期有两个缺失点，出现在t < 5×10^-4 s，这两个点跳变到显示范围之外.其跳变的原因，可能是空气波与数值滤波算法精度的合并影响，也不排除高频修正后的误差带来的影响，但对时间域方法来说，更多关注的是1 ms后的结果，因此本文对此不做更深入的讨论.总体来说，将频域2D有限元结果作高低端频率特性校正后，再用余弦变换法获得2D时间域的结果是可行的.

3 多通道瞬变电磁法数值模拟实例

电磁法对于低阻矿脉都有一定的响应.CSAMT法对顶面埋深较浅、近垂直方向展布的矿脉有较好的分辨能力(Wang et al., 2010)，但对有一定埋深、水平展布的矿脉的分辨力较差，这是由于CSAMT法用频率域资料做解释，而CSAMT法使用频点数常常不够密，水平矿脉引起的异常容易被体积效应圆滑，难于分辨.为此，设计一个顺层展布的2D矿脉模型，利用本文的模拟方法，查验多道瞬变电磁法通过大地脉冲响应场特征分辨这类矿脉的能力.

3.1 矿脉模型与观测装置

首先，设计一个三层模型，第一层模拟第四系松散沉积层，层厚30 m，第二层模拟经历变质作用的沉积层，层厚400 m，第三层模拟花岗岩基底，将此三层模型称为背景模型.然后，在上述背景模型的第二层中，设计一个水平展布的矿脉，水平方向延伸500 m，厚度为40 m，其顶界面埋深为300 m，称为存在异常体的模型.根据经验及文献中不同岩石的电阻率值，将各层及矿体的电阻率进行了赋值，见图 4所示.

若将第一次发射源的位置定为坐标原点，则矿脉两端在地表的投影分别为2000 m和2500 m.在地面布设一条长度为4400 m的观测剖面，位置从100 m到4500 m，测点间距设定为100 m.与2D反射地震勘探的野外装置和高密度电法中的装置类似，采用电性源多次发射的方式，第一次发射源位于0 m点，沿测线移动发射源，每次移动200 m，移动至4600 m结束发射，共计发射24次.每移动一次源，整条剖面上的所有测点记录一次信号，并将源检位置的中点作为记录位置，若以偏移距(源检距)的负数为纵轴，以地面的记录位置为横轴，则得到如图 5所示的观测装置图.

3.2 频率域结果(CSAMT响应)

该模型单个发射源(源位于原点)的频率域响应如图 6所示.在总场断面图 6a和图 6b上，含矿异常的分辨率极低，很难定位矿体异常的位置.但用总场减去背景(层状半空间)场，会看到明显的纯异常场，其横向位置基本准确.

3.3 多通道瞬变电磁时间域结果

选取4个点来展示大地脉冲响应随时间的变化趋势(单个发射源的结果，源位于原点)，这四个点分别距源1000、2000、2300 m和3500 m，偏移距为2000 m的点为矿体左边界在地面的投影，而偏移距为2300 m的点为接近矿体中心的点在地面的投影.大地脉冲响应结果如图 7所示，红线为不存在矿脉时的三层地层的脉冲响应，离散点为含矿脉时大地的脉冲响应，可以看出，虽然二者在第一个脉冲极值处有幅值差异，但若无背景场的对比，从脉冲响应曲线中无法判断矿脉的存在.

因此，和频率域的结果类似，计算二者的偏差，如图 7中黑实线所示，可以看出，二者的偏差突出了矿脉的影响，矿脉引起的异常在其上方最大(见图 7c).

3.4 讨论

从图 6和图 7可以看出，对于有一定埋深的异常体来说，无论从频率域的总场还是时间域的大地脉冲响应，单个发射源(源位于原点)都无法突出矿体的响应，要想突出异常，必须去掉相应背景介质的响应，而在实际工作中，背景介质的电性分布是未知的，所以，应该寻求更好的方式来直观地展示异常.多通道瞬变电磁法的观测装置为多样化地、有效地展示观测结果提供了便利与可能.

考虑到背景为层状介质时，等偏移距的电磁波在各个点的传播速度一致，也就是说，偏移距相同时，每个记录位置相同时间道的穿透深度一致.若有异常体存在，会对层状介质的大地脉冲响应产生干扰，引起异常，因此，可以用记录位置作为横轴，绘制大地脉冲响应的等时剖面曲线，探究矿体的影响.图 8是偏移距为300 m、记录时间为100 ms时三层介质包含和不包含矿体时得到的等时剖面，可以看出，当矿体不存在时，等时剖面曲线是一条直线，而当矿体存在时，可以看到很明显的异常场存在，曲线低值异常的中心对应着矿体的中心在地面的投影(x=2250 m).由此可见，用等时剖面曲线做为显示图件，可以突出矿体异常响应.

3.5 不同偏移距的电场等时曲线特征

为了研究图 4模型中的矿体在图 5所示的观测装置下的异常响应特征，将不同偏移距的等时曲线显示于同一张图中，如图 9所示.图 9给出了4张图片，分别展示了记录时间是0.01、1、100 ms和1 s时5个偏移距的等时剖面曲线.

在图 9a中，t=0.01 ms时，电磁波的扩散仍处于早期，而且第一层为低阻层，所以电磁波穿透较浅，给出的所有偏移距的等时曲线均为直线，说明电磁波还未抵达矿体所在的位置.大地脉冲响应记录的是电场随时间的变化率，偏移距越大，大地脉冲响应幅值越小，说明电场随时间的变化越小.

随着时间的推移，电磁波继续向地下扩散，矿体对大地脉冲响应的扰动逐渐显现.当观测时间为1 ms时，除了偏移距2000 m外，其他4个偏移距的等时曲线均显示了较好的异常.结合图 7b分析，图 7b给出的是偏移距2000 m时的大地脉冲响应，从7b可以看出，当记录时间是1 ms时，矿体引起的异常恰巧几乎为零.因此在图 9b中，偏移距是2000 m的等时曲线是一条直线是合理的.

图 9c是记录时间为100 ms时的等时曲线，其中，偏移距小于等于500 m时，矿体引起的异常均很明显，而偏移距是1000 m时的黄色曲线却趋于成直线，这是因为随着时间的推移，电磁波的扩散范围增大，体积效应使得矿体引起的异常场减小，各记录点电场随时间的变化率越来越趋于相同所致.图 9d中的异常均进一步减小.

在图 9b、9c中，红色曲线代表的100 m偏移距大地脉冲响应等时曲线幅值偏低，由于偏移距小，电磁波的体积效应小，其值反映的主要是第一层的低阻层耦合了矿体的响应，故导致整体幅值较低.当电磁波传播到晚期时，因为体积效应，矿体及第一层低阻引起的异常逐渐被平滑，如图 9d所示.

由图 9可以看出，从早期到晚期，在等时剖面图上，矿体的影响在合适的偏移距处都能分辨，这说明除了对特别早的早期结果矿体不产生影响外，在接近矿体埋深的偏移距上，矿体对各个时间道都有影响.由此说明用等时剖面图展示多通道瞬变电磁法异常体的影响是非常有效的.

3.6 电场等时断面特征

考虑到对有一定埋深且顺层产出的矿体，用等时曲线表示能突出异常的特点，为了充分利用多通道瞬变电磁法所获得的大量资料，抽取出图 9中4个时间点记录的大地脉冲响应资料，形成记录点-偏移距等时断面图，以便直观地了解矿体异常的空间特征，如图 10所示.

从图 10可以看出，在t=1 ms和t=100 ms的等时断面图中，当偏移距小于1500 m时，在2000~ 2500 m之间的记录点处，可以看到矿体引起的扰动.而在t=0.01 ms的早期和t=1 s的晚期，更多体现的是地层的信息，图 9d中矿体引起的异常被幅值更高的地层信息所掩盖.

3.7 磁场的等时曲线特征

在已有文献中(Wright，2002)，只应用了电场资料，而没有用到磁场资料.本文对磁场大地脉冲响应进行了探讨.从图 3可以看出，磁场的大地脉冲响应和电场的大地脉冲响应曲线形态完全不同，对磁场大地脉冲响应的等时资料进行抽取并成图，如图 11所示.为了和电场的等时曲线对比，图 11也给出了4个记录时间的等时曲线，总体来说，矿体引起的磁场异常比电场异常小，且不对称，随着偏移距的增大，不对称性更加明显.在早期(t=0.01 ms)，和电场相同，电磁波还未扩散至矿体，所以无异常；在晚期(t=1 s)，矿体的影响是存在的，如偏移距为500 m时，大地脉冲响应比偏移距为1000 m时幅值略小，应是矿体的影响.

图 11给出的是源在接收点小号点方向时的等时曲线，对于多通道瞬变电磁法来说，还存在源在接收点大号点方向的情况，若将前者的偏移距用正值表示，则后者的偏移距为负值.由于图 11中磁场正偏移距的等时曲线不对称，那么负偏移距的等时曲线也将不对称，针对图 11b中偏移距为1000 m和2000 m的等时曲线，把相应负偏移距的等时曲线绘制在一幅图中探查矿体引起的异常，如图 12所示.

从图 12可以看出，虽然正负偏移距的磁场大地脉冲响应等时曲线是不对称的，但二者呈镜像关系，而且二者的交叉点对应着矿体异常的中心位置，这就为多通道瞬变电磁法应用磁场信息提供了可能.

总体来说，磁场大地脉冲响应等时曲线显示的异常相对电场来说小，且不对称，所以前人不用磁场资料是有一定理论依据的，但是，如果同时应用磁场大地脉响应等时资料的正负偏移距信息，矿体中心的异常位置仍能很好地定位出来.

4 结论及建议

通过本文的研究，得到以下认识：

(1) 针对原有程序在2D宽频带正演模拟高低频的错误结果进行了校正，得到了基本正确的宽频带结果.

(2) 本文提出的电场分量的等偏移距等时曲线与等时断面可以清晰地分辨出有一定埋深且顺层产出的矿脉引起的异常，而且能准确地分辨出矿体异常中心的位置，说明多通道瞬变电磁法对有一定埋深且顺层产出的矿体有较好的分辨能力.

(3) 如果同时应用正负偏移距磁场分量的信息，也能较好地定位矿体中心的异常位置，因此，磁场信息也有一定的利用价值.

(4) 大地脉冲响应的早期仍存在部分跳点，虽然不影响1 ms后的结果，但仍需进一步研究与改善.

致谢  审稿专家给本文提出了建设性意见及建议，对本文的提升起到了至关重要的作用，在此表示感谢.