0 引 言

随着国民经济的快速发展，我国对油气资源的需求飞速增长，油气勘探已从过去的东部，浅层，陆相，构造油气藏向西部，深层，海相，地层油气藏转变.地震勘探中声波速度这一单一岩石物性的局限性逐渐凸显，因此业界越来越强调综合地球物理探勘技术在油气勘探中的作用(刘光鼎，2005；唐建伟等，2008；张春贺等，2011；王学军等，2014).电磁法因其能提供地层电阻率信息，在客观上为区分油气储层(高阻异常)和围岩(低阻背景)提供了可能，近年来受到石油工业界的青睐(曹中林等，2006；罗维斌和汤井田，2008；汤井田等，2008；何展翔和余刚，2008；李予国和Constable，2010；沈金松等，2012).

在众多电磁新方法中，多通道瞬变电磁方法(MTEM)因其对高阻薄层的较强分辨能力，高效的施工方式，以及类地震的数据采集和解释方式而日益受到关注.自Ziolkowski等人(2007a)开始提出MTEM的概念以来，MTEM方法已由陆地走向浅海，接收站由传统的TEM接收方式发展为专用多通道接收机接收和海洋电缆接收，发射波形从最初的方波改进为伪随机码(PRBS)(D’Arienzo et al.，2010；Ziolkowski et al.，2010；Ziolkowski et al.，2011).随着MTEM方法在油气勘探中的成功应用(Ziolkowski et al.，2007b；Ziolkowski et al.，2010)发展适应于MTEM发射源特点的数值工具已成为一项紧迫的任务.

本文在前人的工作基础上(Wang and Hohmann，1993；Commer and Newman，2004；许洋铖等，2012；孙怀凤等，2013)，研究了一种既能适应复杂表层地电模型，又能适应MTEM发射波形的有限差分正演模拟算法.和Wang等人的算法相比，本文有四点不同：第一，初始值是通过解一个三维直流电法正演模拟问题获得的，可以适应复杂表层的地电模型；第二，差分格式利用阶跃响应和脉冲响应的数值关系对二者进行交错递推，可同时获得脉冲响应和阶跃响应而不增加计算量；第三，用Kirchhoff积分延拓代替了基于静磁场假设的地空边界延拓公式，提高了边界条件精度；第四，通过褶积获得最终响应，能很好地适应PRBS源的特点.文章先通过均匀半空间模型验证了所构造差分格式的正确性.在此基础上计算了两种较为复杂模型下的高阻油气薄层的响应.

1 多通道瞬变电磁正演模拟 1.1 MTEM的基本特征

MTEM方法采用接地电流源向地下发射电磁波，在发射源轴向(inline)用阵列式布置的接收机进行标量或者张量观测.标量观测一般只测量inline方向电场分量(Ex)，张量观测可以既测量Ex，Ey两个水平电场，还可以根据需要测量磁场分量.野外工作时一般采取多道观测的方式，即不同接收点 X _R同时记录时域波形e_x(k，X _R，t; X _S，l).与此同时记录发射源发出的电流波形i(k，X _S，t).Ward和Hohmann(1988)等人认为接收电场是发射电流波形i(k，X _S，t)和大地脉冲响应g(X _R; X _S，l，t)的线性褶积.考虑到接收机的系统响应和噪声，观测到的电场表达式为(Ziolkowski et al.，2007b)

其中ir_R(t)是接收机的系统响应，n(k，X _R，t)是观测点的噪声水平.k代表观测次数，一般为了提高信噪比，需在时域进行多次观测并叠加. X _S是发射电极中心点的位置，X _R是接收电极中心点的位置.l是发射电极长度.大地脉冲响应既反映了地下介质电学特性，又包含了观测系统参数.在确定的观测系统下可以通过它来推断地下介质的电学参数.它是MTEM方法中最为关键的概念，也是MTEM数据处理解释的基础.本文正是先计算大地脉冲响应，然后通过(1)式获得接收电场波形的.

1.2 时间离散格式

由于单位脉冲函数是阶跃函数的导数，根据褶积性质可知单位脉冲响应是阶跃响应的时间导数(Ward and Hohmann，1988).若 e=e(X，t)，h=h(X，t)和b =μ h 分别表示电场强度，磁场强度和磁感应强度的阶跃响应，则$\dot b = \frac{{\partial b}}{{\partial t}},\dot h = \frac{{\partial h}}{{\partial t}},\dot e = \frac{{\partial e}}{{\partial t}}$就是他们的脉冲响应.时间离散时，将 e和$\dot b$ 定义在半时间节点上，将 $\dot e$ 定义在时间的整数节点上，从麦克斯韦方程组出发可得到以下时间递推格式为

其中$\tilde \varepsilon $ 是虚拟介电常数.根据Wang和Hohmann(1993)，离散格式(2a)-(2c)必须满足Courant-Friedrichs-Lewy条件才能保持稳定.此时每个时间步的取值为 这里要求α≥3.μ_min是介质的最小磁导率，Δ_min是空间离散后的最小空间步长.可以看出该格式的稳定性受时间步长，空间步长，介质的电学特性的共同约束.

如(2)式所示，计算过程中阶跃响应和脉冲响应相互耦合，交错递推，最终和同时求出二者而不增加额外的计算.递推过程如图 1所示.初始化电场及其导数后，首先利用向前差分获取t=$\frac{1}{2}$Δt₀时刻的电场强度(①步). 然后可由(2b)求出该时间节点上的磁感应强度的脉冲响应(②步)，之后可以利用(2c)把电场的脉冲响应递推到下一个时间步(③步)，最后利用(2a)把电场强度递推到下一个时间半节点上(④步).循环②③④步就可以递推出整个计算区域磁感应强度，磁场强度和电场强度的脉冲响应和电场强度阶跃响应随时间演化过程.

空间离散时采用交错网格，把电场强度定义在网格棱的中点上，把磁感应强度时间导数定义在网格面的中心上，电场强度的时间导数和电场强度的位置重合.其格式类似于Wang和Hohmann(1993)，这里不再赘述.

为增加数值模拟过程中的稳定性，模拟过程分别对磁感应强度的时间导数和电流密度的散度进行了约束(Commer and Newman，2004).由于这里的磁感应强度的时间导数是脉冲响应，磁感应强度时间导数的散度为零和磁感应强度本身的散度为零的意义是相同的，即磁场的力线是无源可循的.一般在电磁法勘探中大地可视为良导体，根据连续性方程可获得对电流密度的散度约束条件.(4a)式约束时可和(2b)式混合使用，即磁感应强度时间导数的水平分量可利用(2b)式求解，其垂直分量利用(4a)式根据水平分量求解出来(Commer and Newman，2004).(4b)式的约束可通过欧姆定律转成对电场的约束，同样可以通过和(2a)式的混合使用来实现，公式为

1.3 源的加载

为提高计算效率，文章直接计算大地脉冲响应，然后利用(1)式合成发射任意波形时的响应.在计算脉冲响应时，本文将非齐次麦克斯韦方程中源的处理问题转换成齐次方程的初值问题.先通过解一个直流电法正演问题获得阶跃响应的初值，然后根据麦克斯韦方程组利用导数关系获得脉冲响应的初值.为获取阶跃响应的初值，需要解(5)式这个泊松方程为

其中 Δ $\phi $= e _dc.进行三维直流电法正演模拟时同样采用了有限差分进行了数值模拟，主要采用了文献(Lowry et al.，1989；Spitzer，1995；Zhao and Yedlin，1996；Wu et al.，2003)中介绍的手段技巧，这里不再赘述.得到直流电法的电场值 e _dc之后，可通过下式获得MTEM数值模拟中的初始值 $\dot b$，$\dot h$，$\dot e$，也就是脉冲响应的初值，公式为 其中 e _s是电流源 j _s本身的电场，可通过欧姆定律求得. $\dot b$ ⁰，$\dot h$ ⁰，$\dot e$ ⁰分别是磁感应强度，磁场强度和电场强度的时间导数(即脉冲响应)的初始值.可以看出(6b)和(6c)是麦克斯韦方程的拟静态近似.

1.4 边界条件

在处理地-地边界时，我们采用了高斯吸收边界.其思想是电磁波传播到边界后按照单边高斯函数缓慢过渡到第一类边界条件.这相当于在边界区加载了一个低通滤波器，高频电磁波很快被吸收，低频电磁波缓慢衰减.虽然这种吸收边界条件并不能完全消除边界处理想导体的吸引效应，但是其较单一的第一类边界条件，计算精度有所提高，且编程简单，特别适合并行实现.

由于电磁场在空气层中满足波动方程，地-空边界处理上，我们可利用地震勘探中常用的波场延拓公式来提高向上延拓的精度.本文采用的是Kirchhoff积分延拓公式，表达式为

其中cosθ=z/R，R=$\sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {y - {y_0}} \right)}^2} + {z^2}} $，是从空中一点(x，y，z)到地面某点(x₀，y₀，0)的距离.是偏移孔径.c=1/$\sqrt {\mu \varepsilon } $为光速.若采用静磁场假设，则可认为光速c为无穷大，此时(7)式就退化成Oristaglio和Hohmann(1984)给出的静磁场假设下波数域磁场延拓公式的空间域直接表达(Adhidjaja and Hohmann，1989).

2 模型验证

为检验所构造的的差分格式及边界条件的正确性，文章 用均匀半空间模型对算法进行了验证.所采用的解析解来源于Ward和Hohmann(1988).使用的模型为电阻率为30 Ω·m 的均匀半空间.发射电流为20安培，发射极距为1000 m，发射方向为X轴正向. Z轴垂直地面向下.模型剖分采用不均匀网格剖分模式——靠近源的区域采用均匀细网格，远离源的区域网格按照等比数列逐渐增大.

图 2是6个不同偏移距处inline方向电场的脉冲响应的计算结果，横纵坐标轴均用的是对数坐标.曲线形态均表现出先快速上升后缓慢下降的趋势(其中偏移距为475 m的曲线由于时间轴的截取只显示出了下降部分).在上升阶段，偏移距越小的接收点曲线上升的速度越快，曲线越陡峭，而偏移距大的接收点曲线则更为平滑，这是由于偏移距越小，电磁波的高频成分越丰富，该接收点的脉冲响应越接近于Δ(t)函数.峰值时刻则明显和偏移距的大小成正相关.下降阶段特别是在晚期，不同偏移距的曲线的衰减形态基本保持一致，这主要是因为晚期的信号以低频为主，低频信号的长波长特点使得偏移距相差不是特别大的接收点电磁场值相差不大.图 3是用公式(1)通过褶积合成的MTEM响应.所使用的电流波形为周期为1秒的5阶PRBS，大地脉冲响应是图 2中偏移距为2975 m处接收点处的数值模拟结果.图 4是使用不同的虚拟介电常数时模拟的2975 m偏移距处的inline方向电场脉冲响应的对比结果.

从图 2中可以看出电场脉冲响应在早期和解析解之间有一定误差，这主要是因为我们的差分格式(2c)考虑了介电常数，而解析解没有考虑这一项.由于在早期的时候频率较高，介电常数所代表的波动项的能量比较强，而由公式(3)可以看出我们使用的虚拟介电常数，在早期的取值相对较大(早期Δt_n较小)，所以在早期的计算误差相对于晚期较大.这种误差的变化在图 3中表现为尖峰处出现误差较大，这是因为尖峰处对应于PRBS的符号翻转时刻，以高频成分为主，而图 2中高频成分主要出现在早期.从图 2中都可以看出越靠近发射源误差越大，越远离发射源误差越小.这主要有两方面的原因：首先近源高频成分丰富，介电常数的影响较大；其次由于模拟的是总场，源的奇异性造成的误差也表现出在源附近大，在远离源的地方小的特征.虚拟介电常数对计算结果的影响在图 4中表现的更为直接.可以看出虚拟介电常数越大，计算结果和解析解之间的误差越大.当虚拟介电常数较小时，这种误差主要体现在早期.当虚拟介电常数足够大时，波动性已经完全不可忽略，计算结果和解析解在各个时刻都会出现很大的差距.计算结果表明当把虚拟介电常数控制的较小时，所构造的差分格式是稳定可靠的.

3 高阻油气薄层的响应

为探究MTEM方法在油气薄层探测中的应用效果，分别构造了层状背景下的油气藏模型和复杂高阻盐丘附近的油气藏模型，进行了数值模拟，分析模型中高阻油气薄层的响应.

3.1 层状背景中的薄层油气藏

使用的模型如图 5所示.该4层层状模型源于Hördt等(2000)的参考模型，它是一个陆地沉积盆地油气模型.电阻率为200 Ω·m的油气藏的顶层埋深为1000 m，厚度为50 m，沿X方向宽1000 m，沿Y方向(走向)长1500 m.发射电流为20安培，发射电极距为1000 m，偶极方向为X轴正向.偶极中心位于图 5中坐标位置为(0，0，0)处.首选计算了没有油气藏时层状模型的大地脉冲响应，然后以此作为背景场计算了由薄层油气藏引起的二次场.图 6和图 7分别展示了不同偏移距处的inline方向电场分量和crossline方向磁感应强度分量的异常场大小.可以看出在油气藏的正上方异常场最大，远离油气藏，异常场逐渐减小.离油气藏越近，异常出现时间越早，也越早出现极大值.这可以理解为背景场的电磁波到达油气藏时发生散射.根据惠更斯原理，油气藏相当于一个二次场源，向外发射散射电磁波场.离油气藏越近，越先接收到散射场，且路径相对较短，异常场的数值越大.inline方向电场的最大异常值有10 mv/km，处于接收机动态范围之内，且异常出现的时间为一毫秒到几百毫秒，理论上可以通过MTEM方法辨识出油气藏的存在.

3.2 复杂高阻盐丘背景中的薄层油气藏

由于空气波的影响，由于在在浅水环境中(海水深度<300 m)，电磁法，特别是频率域电磁法在勘探油气薄层时面临着困难(李予国和Constable，2010；殷长春等，2012；张建国等，2013).为研究浅水区域，利用MTEM进行油气勘探的行性，文章构筑了一个浅水环境的盐丘模型.该盐丘模型是从SEG地震盐丘模型转化而来.根据(Meju et al.，2003；Wirianto et al.，2010)的经验，可以通过一个近似的函数关系把纵波速度转换成电阻率，从而构造出一个电性模型.Meju等采用的转换公式为

得到盐丘电阻率模型(盐丘电阻率约为100 Ω·m，背景电阻率相对较小并缓慢变化，大部分在10 Ω·m以下)后，在模型的正西方向(X轴负向)设置了一高阻油气藏，油气藏的大小为1000 m×2000 m×200 m(X，Y，Z)，油气藏的电阻率为150 Ω·m.海水深度约为150 m到200 m.为了尽可能模拟高阻盐丘对油气藏勘探地干扰，我们将油气藏设置在盐丘的邻接的位置，将发射源置于盐丘的另一侧. 以发射源的中心位置为原点建立直角坐标系，则油气藏的中心位置为(-2750 m，0 m，810 m). 盐丘分布的区域范围约为 -3000 m<x<0 m，-2000 m<y<2000 m，z<3000 m.这里收发距和油气藏的埋深满足2~4倍关系，这里主要是考虑到目前PGS使用的托缆长度约为1500 m，加上发射源和最小偏移距直接一般相距1000 m左右，发射源长度为500 m刚好是发射源中心到拖缆最后一个接受点的距离.图 8展示了该模型对数电阻率的切片.

图 9是SEG盐丘模型正演模拟获得的inline方向电场的响应.发射源为周期为1 s的7阶PRBS码，发射电流为1000 A.发射源的长度为500 m，共激发了5 s.采样速率为24 k.从上至下的偏移距依次为525 m，1025 m，1525 m，2025 m，2525 m，3025 m.由于模型的背景电阻率非常低，对电磁波高频成分衰减的非常严重，所以图中场值的变化已经和发射源的电流波形差距非常大.

图 10是从带波形的场值中提取的inline方向电场的单位脉冲响应的剖面.可以看出图中有两条同相轴.第一条同相轴同横轴平行，走时基本不变，认为它是空气波.空气波的传播路径路径为从发射点垂直向上传播，达到海面后，沿海面以光速传播至接收点正上方，之后垂直向下传播，到达接收点(殷长春等，2012).考虑到实际勘探中的偏移距相对光速很小，所以空气中传播时间可以忽略不计，空气波到达接收点的时间基本上与偏移距无关，主要取决于发射机和接收点到海面的距离(Hunziker et al.，2011).第二条同相轴包括了从低阻背景中扩散到接收点的直达波和被高阻体(包括盐丘和油气藏)反射回来的反射波，和层状介质的情况相同，它们是叠加在一起的.空气波同相轴的幅值明显强于有效同相轴，且在近道，二者是重合在一起的.但偏移距超过1000 m时，二者之间有明显的走时差异，这时可以利用走时上的差异将二者分开，这是MTEM法和频率域海洋电磁法相比最大的优势之一.

图 11和12分别是电场inline方向分量总场脉冲响应和垂直方向分量总场脉冲响应的三维透视图.前者的观测时刻为2 s.后者的观测时刻为0.5 s.从电场inline方向分量的透视图上可以清楚的观测到在盐丘的中心位置场值增大，而在外边界，场值降低，这是因为电阻率分界面处，电流法向分量连续，电场的法向分量不连续.根据欧姆定律则有：电阻率大的一侧场值大，电阻率小的一侧场值小.水平切面相当于对电场做了一个水平投影，所以水平切面上更有利于观测电场的水平分量的变化；而相反垂直切面更有利于电场垂直分量的观测.此外，在电场inline方向分量的透视图上，在油气藏的由边界上也可以观测到同样的现象.从电场垂直方向分量的透视图上可以见到介质背景电阻率近似的垂直的分层情况，在油气藏和盐丘的中心都有非常明显的变化.特别是在油气藏的位置，出现了非常强烈的场值增强，且变化的边界与油气藏的边界吻合度很高.

图 13是偏移距不同的2个观测道处电场inline方向分量脉冲响应的总场和背景场对比曲线.曲线形态上表现为两个波峰，前一个波峰在不同偏移距处对应的走时基本保持不变，且振幅随着偏移距的增加基本呈三次方衰减，这是空气波基本的特征，判断为空气波，后一个波峰则振幅上小的多，且随着偏移距则增加衰减的速度要快于空气波，是对勘探有利的有效波.但可以看出空气波的振幅比有效波强的多.

图 14为电场inline方向分量二次场大地脉冲响应的道时剖面和归一化的道时剖面.如果认为0.3 mv/km的绝对异常和10%的相对异常是比较有利的观测条件的话，则在偏移距为-3000 m到-4000 m(接收点位于发射点的负X方向，定义偏移距为负)之间是本油气模型比较有利的观测位置，这个位置也刚好处于油气藏在海平面上投影的位置.

4 结 论

文章构造了一种适合于MTEM正演计算的新的差分格式.通过解直流电法问题获取初始值，提高了初始值对复杂模型的适应性；用Kirchhoff积分公式代替静磁场假设下的地空边界条件，提高了地空边界的延拓精度.并通过对模型的计算，获得了以下认识：

(1)MTEM方法对高阻薄层油气藏有非常好的分辨能力.在陆上，埋深为1000 m，厚度为50 m的油气薄层可以通过MTEM这种新的电磁技术探测出来.在海洋环境中，MTEM方法特别适合使用拖缆系统进行走航观测.这种观测方式有两个优点，首先数据采集的工作效率特别高；其次，它能克服海洋浅水环境中空气波支配着电磁响应的问题.当把海洋浅表拖拽观测方式应用于MTEM方法时，浅水区电磁勘探的问题可以很好的得到解决.

(2)异常场最大值出现在异常体的正上方，为此勘探过程中接收机应布置于异常体的正上方.偏移距的选择要均衡空气波的影响和接收信号的信噪比，一般而言二者之间是矛盾对立的.最佳偏移距的选择不仅和区域的地质背景有关，也与异常体的电阻率，埋深，大小等因素有关，需要根据测区的具体情况确定.

(3)电场的垂直分量对高阻油气藏有着更好的灵敏度和分辨率，在条件允许的条件下(如井中勘探)，可以尝试观测该分量以增加对油气藏的辨识度.

致 谢 文章作者在进行MTEM有限差分正演模拟的研究过程中得到了Tsili Wang 老师，孙怀凤老师，王妙月老师的热情帮助和指导，并与吴学银同学，严君同学，李展辉同学进行了很多有益的探讨，在此一并表示感谢.