CN104408021A - 一种电偶源三维时域有限差分正演成像方法 - Google Patents

Abstract

本发明涉及一种电偶源三维时域有限差分正演成像方法，其包括：将高斯脉冲加载到电偶极子源上；对海洋空气、海水和海底大地这三分空间均建立麦克斯韦方程组和本构方程；对三分空间的棱柱体模型均进行均匀网格剖分；假定每个剖分得到的网格上电导率和磁导率固定不变，根据剖分得到的网格，采用时域有限差分法得到海水和海底大地的差分方程；采用解析解对海洋空气的麦克斯韦方程组进行处理，计算得到海面以上空气中的电磁场；对剖分空间的边界条件进行处理并设定稳定性条件；结合对边界条件的处理结果和设定的稳定性条件，对建立的差分方程进行求解，得到海水和海底大地中任意时刻电磁场的分布。

Description

一种电偶源三维时域有限差分正演成像方法

技术领域

本发明涉及一种海洋物理电磁勘探方法，特别是关于一种电偶源三维时域有限差分正演成像方法。

背景技术

海洋可控源数值模拟技术不论是在频域还是在时域，主流的方法主要有有限差分法、有限元法和积分方程法三大类。其中，积分方程法是一个非常有效的电磁模拟计算方法，基本思想是将空间电导率分成两部分，背景电导率和异常电导率，背景电导率主要用于求解格林函数，异常电导率即积分区域电导率。积分方程法与其它数值模拟方法相比，最大的优势在于只需要计算异常区域，而不需要对全区域进行计算；其局限性在于目前只能得到简单层状空间的格林函数，所以难于对复杂异常区域模型积分，并且难于对复杂多异常体模型进行计算。积分方程法在地球物理中最早的成功应用见于1975年，由Gerald W.Hohman提出。后来经过S.Zhdanov、Sam C.Ting、Wannamaker P E、Gregory A.Newman、殷长春、Masashi Endo和陈桂波等人的不断丰富和发展，现在已经非常成熟，可以模拟非均匀各向异性介质，即只要能得到介质的电导率、磁导率及介电常数对空间坐标的函数就可以模拟介质内电磁场的传播。也可以进行海洋电磁大尺度模拟。在海洋可控源电磁数值模拟中，2D电磁模拟不适合描述复杂模型的异常特征，所以很少使用，通常都采用3D电磁模拟。近年来，3D电磁模拟算法发展很快，最引人注目的和使用最多的要数犹他大学Zhdanov等人的研究成果(Weiss and Constable,2006；Hursan,and Zhdanov,2002；Zhdanov et al.,2000)。Ueda和Zhdanov(2005)采用一种准线性近似算法进行了3D高精度CSEM(ControlledSource Electromagnetic Methods，可控源电磁法)积分方程法模拟，该算法结合实际数据可以得到很好的结果，但是依据该算法对多个源模拟时却非常耗时。Zhdanov和Wan(2005)采用3D积分方程算法计算了一个非均匀背景下的复杂模型，在岩丘和火成岩等基底已知信息的约束下，可以准确的模拟出油气储藏的清晰图像。Zhdanov和Wan(2005)又用该方法进行了快速CSEM数据偏移成像。结果显示该算法对海洋电阻率异常体成像是很有效的。Zhdanov和Yoshioka(2005)做了MCSEM(MarincControlled Source Electromagnetic Methods，海洋可控源电磁法)数据的3D交互反演，获得很好的效果。但是计算必须在大型机上进行。积分方程法对于多目标、非均匀介质模拟存在明显缺陷，背景场难求，对计算机CPU和内存的要求非常高。而FDTD方法(Finite Difference Time Domain，时域有限差分方法)在达到同样计算精度的同时对计算机CPU和内存的要求低得多。积分方程法模拟目前大多在频率域内进行，转到时间域也是一个大的挑战。

有限元法也是一种非常有效的电磁模拟计算技术，其对复杂模型和复杂边界的处理非常有效。有限元法在地球物理中的应用是非常广泛的，1999年Zhdanov公布了2.5D有限元法在频域时域的正演模拟方法；国内罗延钟、底青云、沈金松、王若，王妙月、熊彬和王华军等对有限元法的研究做了大量工作，取得了许多实用成果。目前，海洋电磁中有限元法使用最好的是由Scripps海洋研究所李予国等人提出的海洋电磁2.5维频域自适应有限元正反演算法，应用效果很好，而且也实现了程序的并行化。但3D有限元模拟难度较大，目前还没有实用化报道。一定意义上有限元法模拟的精度是最高的，可以根据介质的几何变化自适应的剖分空间，但是对CPU和内存的要求很高，因为剖分越密，精度越高，内存需求呈几何倍增。特别是三维有限元模拟是一个非常大的挑战。目前大都在频率域内研究，要想转到时间域就更困难了。

与积分方程法和有限元法相比，本发明所采用的时域有限差分法(FDTD)数学理论没有前两者复杂，虽然需要在整个模拟空间进行剖分，但其具有天然的并行性，非常适合大尺度模拟计算，而且很容易适应介质特性，理论上可以模拟任意复杂介质，最主要的是它可以给出全空间全时段电磁场的分布。FDTD在低频地球物理勘探中最早见于Hohman(1984、1985、1989、1993)等人的研究，他们针对地面可控源电磁勘探用FDTD对2D、3D的模型分别进行了研究，得到了比较好的结果Oristagl io&Hohmann1984给出了2D(二维)FDTD解，Adhidjaja&Hohmann 1989给出了3D(三维)FDTD解，Tislliwang 1993年进一步做了研究，将FDTD成功用于地面瞬变电磁法勘探及井地电磁探测中，对计算机要求不高的情况下，给出了均匀半空间中低阻异常的电磁响应，与解析解相比结果令人满意。2002年以来国内也有学者研究此方法，闫述2002，2004，2005，2007，2009做了2D、3D的研究，并且在2007年论述了低频FDTD要解决的几个重要问题，如激励源的加入和模拟，边界条件处理，数值稳定性等。在有耗介质中使用3D FDTD模拟，史红锫的硕士论文中研究了采用线性电流源代替瞬变电流的电磁响应，得到了跟解析解相比较满意的结果。史红锫将Tislliwang等的方法做了一些改进，不需再求解初始值，而是直接将线性源电流加入方程，但其并没有详细论证和叙述电流的波形及频谱适用的范围。李貅、薛国强、石显新等发表论文针对陆地煤层地质瞬变电磁响应，分别做了2D、3D的研究。2002年闫玉波也做了2D的研究，之后刘云2012年做了2D长线源FDTD的研究，以上这些研究全部是在地面上做的，共同特征是将源置于地面，并且采用长线源和磁偶极子。前人的研究充分说明FDTD在低频地球物理勘探数值模拟中应用的有效性。但是存在不足，(1)不能适应地面起伏地形；(2)激励源加入没有得到很好解决，发射信号频率也不可控；(3)未涉及海洋电磁的模拟。

发明内容

针对上述问题，本发明的目的是提供一种采用高斯脉冲作为发射电流信号的海洋电磁三维时域有限差分正演成像方法。

为实现上述目的，本发明采取以下技术方案：一种电偶源三维时域有限差分正演成像方法，其包括以下步骤：1)将高斯脉冲加载到电偶极子源上，电偶极子源的源电流I(t)为：

$I\left(t\right)=\exp (-\frac{{(t-t_0)}^2}{{\mathrm{\τ}}^2}),$

式中，t为时间，t₀为时移常数，τ为脉冲宽度因子，取τ＝0.1s；2)对海洋空气、海水和海底大地这三分空间均建立如下麦克斯韦方程组：

$\left\{\begin{array}{c}\▿\×E(r,t)=-\frac{\∂B(r,t)}{\∂t}\\ \▿\×H(r,t)=J(r,t)\\ \▿\·B(r,t)=0\\ \▿\·J(r,t)=0\end{array}\right.,$

对海洋空气、海水和海底大地均建立如下本构方程：

$\left\{\begin{array}{c}B(r,t)=\mathrm{\μ}\left(r\right)H(r,t)\\ J(r,t)=\mathrm{\σ}\left(r\right)E(r,t)\end{array}\right.,$

式中，B(r,t)、H(r,t)、E(r,t)分别表示磁感应强度、磁场强度、电场强度。μ(r)表示介质中磁导率，σ(r)表示介质中电导率，J(r,t)表示介质中的传导电流密度，r表示测点到原点的距离，3)对海洋空气、海水和海底大地这三分空间的棱柱体模型均进行均匀网格剖分；4)假定每个剖分得到的网格上电导率和磁导率固定不变，根据步骤3)剖分得到的网格，采用时域有限差分法对海水和海底大地这两部分空间的麦克斯韦方程组均进行差分处理，得到差分方程；5)根据步骤4)得到的海水和海底大地中的差分方程，采用解析解对海洋空气的麦克斯韦方程组进行处理，计算得到海面以上空气中的电磁场；6)为得到步骤4)中建立的差分方程的稳定解，对剖分空间的边界条件进行处理并设定稳定性条件；7)结合步骤6)对边界条件的处理结果和设定的稳定性条件，对步骤4)建立的差分方程进行求解，得到海水和海底大地中任意时刻电磁场的分布。

所述步骤1)中，高斯脉冲在时域中的脉冲宽度取6τ，其频谱宽度取4τ/5。

所述步骤2)中，对海洋空气、海水和海底大地棱柱体模型进行均匀网格剖分时，最大空间步长为45m。

所述步骤4)中，采用时域有限差分法对海水和海底大地这两部分空间的麦克斯韦方程组均进行差分处理，得到海水和海底大地中的电场强度E(r,t)的差分方程为：

$\begin{array}{c}{E}_x^{n+1}(i+1/2,j,k)=\frac{2\mathrm{\γ}-{\mathrm{\Δt}}_n\mathrm{\σ}(i+1/2,j,k)}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i+1/2,j,k)}{E}_x^n(i+1/2,j,k)\\ +\frac{{2\mathrm{\Δt}}_n}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i+1/2,j,k)}[\frac{B_z^{n+1/2}(i+1/2,j+1/2,k)-B_z^{n+1/2}(i+1/2,j-1/2,k)}{({\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j-1})/2}\\ -\frac{B_y^{n+1/2}(i+1,j,k+1/2)-B_y^{n+1/2}(i+1/2,j,k-1/2)}{({\mathrm{\Δz}}_j+{\mathrm{\Δz}}_{j-1})/2}]\end{array},$

$\begin{array}{c}{E}_y^{n+1}(i,j+1/2,k)=\frac{2\mathrm{\γ}-{\mathrm{\Δt}}_n\mathrm{\σ}(i,j+1/2,k)}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j+1/2,k)}{E}_y^n(i,j+1/2,k)\\ +\frac{{2\mathrm{\Δt}}_n}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j+1/2,k)}[\frac{E_x^{n+1/2}(i,j+1/2,k+1/2)-E_x^{n+1/2}(i,j+1/2,k-1/2)}{({\mathrm{\Δz}}_k+{\mathrm{\Δz}}_{k-1})/2}\\ -\frac{B_z^{n+1/2}(i+1/2,j+1/2,k)-B_z^{n+1/2}(i-1/2,j+1/2,k)}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i-1})/2}]\end{array},$

$\begin{array}{c}{E}_z^{n+1}(i,j,k+1/2)=\frac{2\mathrm{\γ}-{\mathrm{\Δt}}_n\mathrm{\σ}(i,j,k+1/2)}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j,k+1/2)}{E}_z^n(i,j,k+1/2)\\ +\frac{{2\mathrm{\Δt}}_n}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j,k+1/2)}[\frac{B_y^{n+1/2}(i+1/2,j,k+1/2)-B_y^{n+1/2}(i-1/2,j,k+1/2)}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i-1})/2}\\ -\frac{B_x^{n+1/2}(i+1/2,j,k+1/2)-B_x^{n+1/2}(i,j-1/2,k+1/2)}{({\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j-1})/2}]\end{array},$

式中，γ表示虚拟介电常数，Δt_n为时间步长，Δ为空间步长，μ为磁导率；E_x、E_y和E_z分别为电场强度E(r,t)在X轴、Y轴和Z轴的分量；海水和海底大地中的磁感应强度B(r,t)的差分方程为：

$\begin{array}{c}{B}_x^{n+1/2}(i,j+1/2,k+1/2)=B_x^{n-1/2}(i,j+1/2,k+1/2)\\ -\frac{{\mathrm{\Δt}}_{n-1}+{\mathrm{\Δt}}_n}{2}\×[\frac{E_z^n(i,j+1,k+1/2)-E_z^n(i,j,k+1/2)}{{\mathrm{\Δy}}_j}\\ -\frac{E_y^n(i,j+1/2,k+1)-E_y^n(i,j+1/2,k)}{{\mathrm{\Δz}}_k}]\end{array},$

$\begin{array}{c}{B}_y^{n+1/2}(i+1/2,j,k+1/2)=B_y^{n-1/2}(i+1/2,j,k+1/2)\\ -\frac{{\mathrm{\Δt}}_{n-1}+{\mathrm{\Δt}}_n}{2}\×[\frac{E_x^n(i+1/2,j,k+1)-E_x^n(i+1/2,j,k)}{{\mathrm{\Δz}}_k}\\ -\frac{E_z^n(i+1,j,k+1/2)-E_z^n(i,j,k+1/2)}{{\mathrm{\Δx}}_i}]\end{array},$

$\begin{array}{c}{B}_z^{n+1/2}(i+1/2,j+1/2,k)=B_z^{n+1/2}(i+1/2,j+1/2,k+1)\\ +{\mathrm{\Δz}}_k[\frac{B_x^{n+1/2}(i+1,j+1/2,k+1/2)-B_x^{n+1/2}(i,j+1/2,k+1/2)}{{\mathrm{\Δx}}_i}\\ +\frac{B_y^{n+1/2}(i+1/2,j+1,k+1/2)-B_y^{n+1/2}(i+1/2,j,k+1/2)}{{\mathrm{\Δy}}_j}]\end{array},$

式中，B_x、B_y和B_z分别为磁感应强度B(r,t)在X轴、Y轴和Z轴的分量；海水和海底大地中的磁场强度H(r,t)的差分方程为：

$\left\{\begin{array}{c}{H}_x(i,j+1/2,k+1/2)=B_x(i,j+1/2,k+1/2)/\mathrm{\μ}(i,j+1/2,k+1/2)\\ H_y(i+1/2,j,k+1/2)=B_y(i+1/2,j,k+1/2)/\mathrm{\μ}(i+1/2,j,k+1/2)\\ H_z(i+1/2,j+1/2,j+1/2,k)=B_z(i+1/2,j+1/2,k)/\mathrm{\μ}(i+1/2,j+1/2,k)\end{array}\right.,$

式中， $\mathrm{\μ}(i,j+1/2,k+1/2)=\frac{{\mathrm{\Δx}}_{i-1}\mathrm{\μ}(i-1,j,k)+{\mathrm{\Δx}}_i\mathrm{\μ}(i,j,k)}{{\mathrm{\Δx}}_{i-1}+{\mathrm{\Δx}}_i},$ i表示各网格点在X轴方向的位置，j表示各网格点在Y轴方向的位置，k表示各网格点在Z轴方向的位置，i＝i_min,i_min+1,…,i_max；j＝j_min,j_min+1,…,j_max；k＝k_min,k_min+1,…,k_max；H_x、H_y和H_z分别为磁场强度H(r,t)在X轴、Y轴和Z轴的分量。

所述步骤5)中，采用解析解计算得到海面以上空气中的电磁场，其具体包括以下步骤：I)对海面以上的电磁场进行二维傅里叶变换，从空间域到波数域；对磁感应强度B(r,t)在Z轴的分量B_z做二维傅里叶变换，得到波数域将波数域乘以操作因子计算得到波数域将波数域乘以操作因子计算得到波数域II)计算波数域中海面以上半空间步长处水平方向的磁感应强度；由磁感应强度B(r,t)在自由空间中遵循矢量Laplace方程：以及B_z的2D空间傅里叶变换得到波数域方程组：

式(1)和式(2)中，表示B_x的2D空间傅里叶变换，表示B_y的2D空间傅里叶变换；k_x表示相对于x坐标的波数域变量，k_y表示相对于y坐标的波数域变量；利用式(1)和式(2)在海平面向上连续计算得到在自由空间中和的值：

将式(3)和式(4)分别代入式(1)和式(2)中，得到：

式(3)～式(6)中，exp表示以自然对数e为底的指数函数；III)对空气中的磁感应强度进行二维傅里叶逆变换，从波数域到空间域；对进行二维空间傅里叶反变换，得到空间域的B_x，对进行二维空间傅里叶反变换，得到空间域的B_y；IV)在海面以下采用时域有限差分法计算边界处的切向电场；通过海面上B_z及海面以上半空间步长处的B_x、B_y，用差分方程计算海面上切向电场E_x、E_y。

所述步骤6)中，对剖分空间的边界条件进行处理，其具体包括：(I)对海洋空气与海水之间的边界进行处理；利用电磁场边界条件，介质面上磁场切向连续，磁感应强度法向连续的基本规律，麦克斯韦方程组隐含着两个连续性边界条件，电磁场在穿过不连续的介质面时遵循以下边界条件：电场和磁场的切向分量连续，即E_1ta＝E_2ta，H_1ta＝H_2ta，ta表示切向；电通量密度和磁感应强度法向分量连续，即D_1n＝D_2n，B_1n＝B_2n，n表示法向，或者ε₁E_1n＝ε₂E_2n，μ₁H_1n＝μ₂H_2n；(II)对剖分空间的前边界、剖分空间的后边界、剖分空间的左边界、剖分空间的右边界以及海底大地下无穷远处的边界进行处理；剖分空间的前、后、左、右边界以及海底大地下无穷远处的边界都采用狄利克雷吸收边界，看成在无穷远处，电磁场的值衰减为0；(III)在介质的边界面上，采用电场的切向连续和磁感应强度的法向连续对海水与海底大地之间的边界进行处理。

所述步骤6)中，设定的稳定性条件为：

${\mathrm{\Δt}}_{\max }=\mathrm{\α}{\left(\frac{{\mathrm{\μ}}_{\min }\mathrm{\σt}}{6}\right)}^{1/2}{\mathrm{\Δ}}_{\min },$

式中，系数α的取值范围是0.1～0.2，Δt_max为最大时间步长，Δ_min为最小空间步长，Δ_min取决于计算的精度要求，μ_min为最小介质磁导率，σ为介质电导率。

本发明由于采取以上技术方案，其具有以下优点：1、本发明采用高斯电流脉冲作为发射信号，其低频信号丰富，信号频谱易于控制；针对不同探测深度，所用高斯电流脉冲的有效频率也不同，可以灵活调整高斯电流脉冲的宽度；高斯电流脉冲用于正演模拟，其n阶导数连续，能够很好地保持迭代方程的稳定性；高斯脉冲在硬件技术上更容易实现。2、本发明可以模拟任意复杂导电介质中低频电磁波的传播特征，展示任意时刻空间电磁场的分布，既可以观察其总场分布，也便于观察其二次场响应；同时本发明也适用于地面电磁勘探数值模拟，适用范围广。基于以上优点，本发明可以广泛应用于海洋电磁勘探、陆地电磁勘探和航空电磁勘探的三维正演成像中。

附图说明

图1是本发明采用的高斯脉冲的波形图及其频谱图；其中，图(a)是高斯脉冲的波形图，图(b)是高斯脉冲的频谱图；

图2是海洋空气、海水和海底大地的空间网格剖分示意图；其中，图(a)是海洋空气的空间网格剖分示意图，图(b)是海水的空间网格剖分示意图，图(c)是海底大地的空间网格剖分示意图；

图3是FDTD中的Yee元胞图；

图4是海洋空气与海水之间的边界处理流程图；

图5是常规海洋电磁三维构造勘探模型；

图6是海底地下有不同电导率三维异常体时E_x二次场20道时序波形堆积图；其中，图(a)是棱柱体为高阻体时的时序图，图(b)是棱柱体为低阻体时的时序图；

图7是海底地下有不同电导率三维体时MVO曲线比较；其中，图(a)是填充高阻异常与背景场的MVO曲线对比，图(b)是填充低阻异常与背景场的MVO曲线对比；

图8是海底地下有高阻三维体时E_x在不同时刻的波场快照；其中，图(a)～图(d)分别是1.03s、2.17s、3.72s和5.69s四个时刻全空间的波场快照；

图9是海底地下有低阻三维体时E_x在不同时刻的波场快照；其中，图(a)～图(d)分别是1.03s、2.17s、3.72s和5.69s四个时刻全空间的波场快照。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明电偶源三维时域有限差分正演成像方法包括以下步骤：

1)将高斯电流脉冲加载到电偶极子源上

将如图1所示的最高频率为10Hz的高斯脉冲加载到电偶极子源上，其中电偶极子源的源电流I(t)为：

$I\left(t\right)=\exp (-\frac{{(t-t_0)}^2}{{\mathrm{\τ}}^2})---\left(1\right)$

式中，t为时间，t₀为时移常数，τ为脉冲宽度因子，一般取τ＝0.1s。

为使高斯脉冲在初始时刻近似为零，高斯脉冲在时域中的脉冲宽度取6τ，则时移常数t₀取t₀＝3τ。高斯脉冲的频谱宽度取

2)对海洋空气、海水和海底大地这三分空间均建立如下麦克斯韦方程组：

$\left\{\begin{array}{c}\▿\×E(r,t)=-\frac{\∂B(r,t)}{\∂t}\\ \▿\×H(r,t)=J(r,t)\\ \▿\·B(r,t)=0\\ \▿\·J(r,t)=0\end{array}\right.---\left(2\right)$

对海洋空气、海水和海底大地均建立如下本构方程：

$\left\{\begin{array}{c}B(r,t)=\mathrm{\μ}\left(r\right)H(r,t)\\ J(r,t)=\mathrm{\σ}\left(r\right)E(r,t)\end{array}\right.---\left(3\right)$

式中，B(r,t)、H(r,t)、E(r,t)分别表示磁感应强度、磁场强度、电场强度。μ(r)表示介质中磁导率，σ(r)表示介质中电导率，J(r,t)表示介质中的传导电流密度，r表示测点到原点的距离，

3)对三分空间进行网格剖分；

对海洋空气、海水和海底大地这三分空间的棱柱体模型均进行均匀网格剖分。如图2所示，建立空间直角坐标系O-XYZ，其中X轴水平向右为正，Y轴垂直于纸面向外为正，Z轴竖直向下为正。将一定空间内的海洋空气、海水和海底大地均抽象为一棱柱体模型，在空间直角坐标系O-XYZ的第一象限分别对海洋空气、海水和海底大地棱柱体模型进行均匀网格剖分，最大空间步长为45m。其中，i表示各网格点在X轴方向的位置，j表示各网格点在Y轴方向的位置，k表示各网格点在Z轴方向的位置，i＝i_min,i_min+1,…,i_max；j＝j_min,j_min+1,…,j_max；k＝k_min,k_min+1,…,k_max。

根据海洋电磁环境中空气和导电介质中的色散确定最小网格步长。

对于三维空间模拟，如果最小空间步长Δ_min和波长λ满足以下关系：

$\left\{\begin{array}{c}{\mathrm{\Δ}}_{\min }\≤\frac{\mathrm{\λ}}{N}\\ N\≥12\end{array}\right.---\left(4\right)$

则差分近似带来的色散将非常小。

I)空气中的色散

三分空间的空气中波速为光速，根据波长λ、频率f和光速c三者之间的关系λf＝c，按照发射的高斯脉冲最高频率为10Hz，计算得到波长的最小值为λ_min＝3×10⁷m。将λ_min＝3×10⁷m代入式(4)，则最小空间步长Δ_min可选为2.5×10⁶m。空气中不存在色散，因此可以不再考虑空气中的色散问题。

II)导电介质中的色散

导电介质中波长λ的表达式为：

$\mathrm{\λ}=2\mathrm{\π}\sqrt{\frac{2}{\mathrm{\μ\ω\σ}}}---\left(5\right)$

式中，μ为介质磁导率，ω为角频率，σ为介质电导率。

将海水的电导率3.3s/m、磁导率4π×10^-7H/m以及高斯脉冲的最高频率10Hz代入式(5)中，计算得到波长的最小值为λ_min＝550.5m。将λ_min＝550.5m代入式(4)，则最小空间步长Δ_min可选为45m。对于高斯脉冲电流来说，高频出现在中期，数值色散更可能发生在电磁响应的中期，引起波形畸变，或者使脉冲展宽。

4)假定每个剖分得到的网格上电导率和磁导率固定不变，根据步骤3)剖分得到的网格，采用时域有限差分法对海水和海底大地这两部分空间的麦克斯韦方程组均进行差分处理，得到差分方程；

如图3所示，将剖分后的每一个网格均看作一个Yee元胞，其中电场处在元胞的棱边上，磁场处在元胞各个面的中心，这样每一个电场分量周围有四个磁场分量，每一个磁场分量周围也有四个电场分量，电场环绕磁场，磁场环绕电场，这样整个剖分空间都是由此元胞构成，因而整个剖分空间的电磁场分量都可以由元胞的棱边和面表示出来。其中电场和磁场相差半个时间步长。整个剖分空间的电磁场每一个时间步长刷新一遍，遵循电磁波在空间的传播规律，涡旋电场产生涡旋磁场，涡旋磁场产生涡旋电场，随着时间的推移，电磁波从偶极源向周围传播开来。

根据Yee元胞中反映的电场与磁场之间的关系，对步骤2)建立的麦克斯韦方程组进行时空离散，得到海水和海底大地中的电场强度E(r,t)的差分方程均为：

$\begin{array}{c}{E}_x^{n+1}(i+1/2,j,k)=\frac{2\mathrm{\γ}-{\mathrm{\Δt}}_n\mathrm{\σ}(i+1/2,j,k)}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i+1/2,j,k)}{E}_x^n(i+1/2,j,k)\\ +\frac{{2\mathrm{\Δt}}_n}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i+1/2,j,k)}[\frac{B_z^{n+1/2}(i+1/2,j+1/2,k)-B_z^{n+1/2}(i+1/2,j-1/2,k)}{({\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j-1})/2}\\ -\frac{B_y^{n+1/2}(i+1,j,k+1/2)-B_y^{n+1/2}(i+1/2,j,k-1/2)}{({\mathrm{\Δz}}_j+{\mathrm{\Δz}}_{j-1})/2}]\end{array}---\left(6\right)$

$\begin{array}{c}{E}_y^{n+1}(i,j+1/2,k)=\frac{2\mathrm{\γ}-{\mathrm{\Δt}}_n\mathrm{\σ}(i,j+1/2,k)}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j+1/2,k)}{E}_y^n(i,j+1/2,k)\\ +\frac{{2\mathrm{\Δt}}_n}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j+1/2,k)}[\frac{E_x^{n+1/2}(i,j+1/2,k+1/2)-E_x^{n+1/2}(i,j+1/2,k-1/2)}{({\mathrm{\Δz}}_k+{\mathrm{\Δz}}_{k-1})/2}\\ -\frac{B_z^{n+1/2}(i+1/2,j+1/2,k)-B_z^{n+1/2}(i-1/2,j+1/2,k)}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i-1})/2}]\end{array}---\left(7\right)$

$\begin{array}{c}{E}_z^{n+1}(i,j,k+1/2)=\frac{2\mathrm{\γ}-{\mathrm{\Δt}}_n\mathrm{\σ}(i,j,k+1/2)}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j,k+1/2)}{E}_z^n(i,j,k+1/2)\\ +\frac{{2\mathrm{\Δt}}_n}{2\mathrm{\γ}+{\mathrm{\Δt}}_n\mathrm{\σ}(i,j,k+1/2)}[\frac{B_y^{n+1/2}(i+1/2,j,k+1/2)-B_y^{n+1/2}(i-1/2,j,k+1/2)}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i-1})/2}\\ -\frac{B_x^{n+1/2}(i+1/2,j,k+1/2)-B_x^{n+1/2}(i,j-1/2,k+1/2)}{({\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j-1})/2}]\end{array}---\left(8\right)$

式(6)～式(8)中，γ表示虚拟介电常数，其中Δt_n为时间步长，Δ为空间步长，E_x、E_y和E_z分别为电场强度E(r,t)在X轴、Y轴和Z轴的分量。

海水和海底大地中的磁感应强度B(r,t)的差分方程均为：

$\begin{array}{c}{B}_x^{n+1/2}(i,j+1/2,k+1/2)=B_x^{n-1/2}(i,j+1/2,k+1/2)\\ -\frac{{\mathrm{\Δt}}_{n-1}+{\mathrm{\Δt}}_n}{2}\×[\frac{E_z^n(i,j+1,k+1/2)-E_z^n(i,j,k+1/2)}{{\mathrm{\Δy}}_j}\\ -\frac{E_y^n(i,j+1/2,k+1)-E_y^n(i,j+1/2,k)}{{\mathrm{\Δz}}_k}]\end{array}---\left(9\right)$

$\begin{array}{c}{B}_y^{n+1/2}(i+1/2,j,k+1/2)=B_y^{n-1/2}(i+1/2,j,k+1/2)\\ -\frac{{\mathrm{\Δt}}_{n-1}+{\mathrm{\Δt}}_n}{2}\×[\frac{E_x^n(i+1/2,j,k+1)-E_x^n(i+1/2,j,k)}{{\mathrm{\Δz}}_k}\\ -\frac{E_z^n(i+1,j,k+1/2)-E_z^n(i,j,k+1/2)}{{\mathrm{\Δx}}_i}]\end{array}---\left(10\right)$

$\begin{array}{c}{B}_z^{n+1/2}(i+1/2,j+1/2,k)=B_z^{n+1/2}(i+1/2,j+1/2,k+1)\\ +{\mathrm{\Δz}}_k[\frac{B_x^{n+1/2}(i+1,j+1/2,k+1/2)-B_x^{n+1/2}(i,j+1/2,k+1/2)}{{\mathrm{\Δx}}_i}\\ +\frac{B_y^{n+1/2}(i+1/2,j+1,k+1/2)-B_y^{n+1/2}(i+1/2,j,k+1/2)}{{\mathrm{\Δy}}_j}]\end{array}---\left(11\right)$

式(9)～式(11)中，B_x、B_y和B_z分别为磁感应强度B(r,t)在X轴、Y轴和Z轴的分量。

根据磁感应强度B(r,t)和磁导率μ计算得到海水和海底大地中的磁场强度H(r,t)的差分方程为：

$\left\{\begin{array}{c}{H}_x(i,j+1/2,k+1/2)=B_x(i,j+1/2,k+1/2)/\mathrm{\μ}(i,j+1/2,k+1/2)\\ H_y(i+1/2,j,k+1/2)=B_y(i+1/2,j,k+1/2)/\mathrm{\μ}(i+1/2,j,k+1/2)\\ H_z(i+1/2,j+1/2,j+1/2,k)=B_z(i+1/2,j+1/2,k)/\mathrm{\μ}(i+1/2,j+1/2,k)\end{array}\right.---\left(12\right)$

式中， $\mathrm{\μ}(i,j+1/2,k+1/2)=\frac{{\mathrm{\Δx}}_{i-1}\mathrm{\μ}(i-1,j,k)+{\mathrm{\Δx}}_i\mathrm{\μ}(i,j,k)}{{\mathrm{\Δx}}_{i-1}+{\mathrm{\Δx}}_i},$ H_x、H_y和H_z分别为磁场强度H(r,t)在X轴、Y轴和Z轴的分量。

5)根据步骤4)得到的海水和海底大地中的差分方程，采用解析解对海洋空气的麦克斯韦方程组进行处理，计算得到海面以上空气中的电磁场，其具体包括以下步骤：

I)对海面以上的电磁场进行二维傅里叶变换，从空间域到波数域；

对磁感应强度B(r,t)在Z轴的分量B_z做二维傅里叶变换，得到波数域将波数域乘以操作因子计算得到波数域将波数域乘以操作因子计算得到波数域

II)计算波数域中海面以上半空间步长处水平方向的磁感应强度；

在准静态场条件下，磁感应强度B(r,t)在自由空间中遵循矢量Laplace(拉普拉斯)方程：

${\▿}^{2}B=0---\left(13\right)$

满足方程(13)的磁场，其水平分量可以由垂直分量导出。由式(13)可以推出如下波数域方程组：

式(14)和式(15)中，表示B_x的2D空间傅里叶变换，表示B_y的2D空间傅里叶变换，表示B_z的2D空间傅里叶变换；k_x表示相对于x坐标的波数域变量，k_y表示相对于y坐标的波数域变量。

利用式(14)和式(15)在海平面向上连续计算得到在自由空间中和的值。

将式(16)和式(17)分别代入式(14)和式(15)中，得到：

式(16)～式(19)中，exp表示以自然对数e为底的指数函数。

III)对空气中的磁感应强度进行二维傅里叶逆变换，从波数域到空间域；

对进行二维空间傅里叶反变换，得到空间域的B_x，对进行二维空间傅里叶反变换，得到空间域的B_y。

IV)在海面以下采用时域有限差分法计算边界处的切向电场；

通过海面上B_z及海面以上半空间步长处的B_x、B_y，用差分方程计算海面上切向电场E_x、E_y。

6)为得到步骤4)中建立的差分方程的稳定解，对剖分空间的边界条件进行处理并设定稳定性条件，其具体包括：

I)对剖分空间的边界条件进行处理

整个剖分空间包括六个外边界和一个内边界，其中六个外边界包括海洋空气与海水之间的边界、剖分空间的前边界、剖分空间的后边界、剖分空间的左边界、剖分空间的右边界和海底大地下无穷远处的边界，内边界为海水与海底大地之间的边界。

(I)如图4所示，对海洋空气与海水之间的边界进行处理。

利用电磁场边界条件，介质面上磁场切向连续，磁感应强度法向连续的基本规律，麦克斯韦方程组隐含着两个连续性边界条件，电磁场在穿过不连续的介质面时遵循以下边界条件：

电场和磁场的切向分量(水平分量)连续，即E_1ta＝E_2ta，H_1ta＝H_2ta，ta表示切向。

电通量密度和磁感应强度法向分量(垂直分量)连续，即D_1n＝D_2n，B_1n＝B_2n，n表示法向，或者ε₁E_1n＝ε₂E_2n，μ₁H_1n＝μ₂H_2n。

在差分离散中，利用电磁场在穿过不连续的介质面时遵循的边界条件来确定交错网格上电磁场的连续性。海面以下电磁场根据差分方程(6)～(12)计算得到，海面以上的电磁场根据方程(14)～(19)的解析解计算得到，海面以下电磁场和海面以上的电磁场通过边界条件进行连接。

(II)对剖分空间的前边界、剖分空间的后边界、剖分空间的左边界、剖分空间的右边界以及海底大地下无穷远处的边界进行处理；

剖分空间的前、后、左、右边界以及海底大地下无穷远处的边界都采用狄利克雷吸收边界，看成在无穷远处，电磁场的值衰减为0。

(III)对海水与海底大地之间的边界进行处理；

FDTD原则上可以模拟任意介质。海水和海底大地都属于低阻介质，在介质的边界面上采用电场的切向连续和磁感应强度的法向连续处理，迭代方程不变，通过界面点上的电阻率和磁导率加权平均便可自然过渡，也就是差分方程不变，只是调整界面处介质电阻率和电导率值即可。

如求取电导率：

$\left\{\begin{array}{c}\mathrm{\σ}(i+1/2,j,k)=\mathrm{\σ}(i,j,k)w(i,j,k)+\mathrm{\σ}(i,j-1,k)w(i,j-1,k)\\ +\mathrm{\σ}(i,j,k-1)w(i,j,k-1)+\mathrm{\σ}(i-1,j,k-1)w(i-1,j,k-1)\\ \mathrm{\σ}(i,j+1/2,k)=\mathrm{\σ}(i,j,k)w(i,j,k)+\mathrm{\σ}(i-1,j,k)w(i-1,j,k)\\ +\mathrm{\σ}(i,j,k-1)w(i,j,k-1)+\mathrm{\σ}(i-1,j,k-1)w(i-1,j,k-1)\\ \mathrm{\σ}(i,j,k+1/2)=\mathrm{\σ}(i,j,k)w(i,j,k)+\mathrm{\σ}(i,j-1,k)w(i,j-1,k)\\ +\mathrm{\σ}(i-1,j,k)w(i-1,j,k)+\mathrm{\σ}(i-1,j-1,k)w(i-1,j-1,k)\end{array}\right.---\left(20\right)$

式中，w是由磁环穿起来的四个网格对环内电导率的贡献率。这里假定在整个磁环区域内电场及其对时间的导数是不变的，每一磁场分量在其所在的磁环边上磁场是恒定的。

以σ(i+1/2,j,k)为例，电场处在一个磁环中，电导率为磁环所穿过的四个网格电导率的加权平均值，其中

$w(i,j,k)=\frac{{\mathrm{\Δy}}_j{\mathrm{\Δz}}_k}{({\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j-1})({\mathrm{\Δz}}_k+{\mathrm{\Δy}}_{k-1})},$

$w(i,j-1,k)=\frac{{\mathrm{\Δy}}_{j-1}{\mathrm{\Δz}}_k}{({\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j-1})({\mathrm{\Δz}}_k+{\mathrm{\Δy}}_{k-1})},$

而磁导率则是相邻两个网格的磁导率加权平均值，例如X方向的磁导率可以表示为：

$\mathrm{\μ}(i,j+1/2,k+1/2)=\frac{{\mathrm{\Δx}}_{i-1}\mathrm{\μ}(i-1,j,k)+{\mathrm{\Δx}}_i\mathrm{\μ}(i,j,k)}{{\mathrm{\Δx}}_{i-1}+{\mathrm{\Δx}}_i}---\left(21\right)$

实际计算中，空间的介质参数都用式(20)和式(21)计算得到，如果是均匀介质，加权平均后还是其实际值。

II)为使步骤4)建立的差分方程的解在空间任意一点和任意时刻都一致趋于原麦克斯韦方程组的解，设定稳定性条件：

${\mathrm{\Δt}}_{\max }=\mathrm{\α}{\left(\frac{{\mathrm{\μ}}_{\min }\mathrm{\σt}}{6}\right)}^{1/2}{\mathrm{\Δ}}_{\min }---\left(22\right)$

式中，系数α的取值范围是0.1～0.2，Δt_max为最大时间步长，Δ_min为最小空间步长，Δ_min取决于计算的精度要求，μ_min为最小介质磁导率，σ为介质电导率。

7)结合步骤6)对边界条件的处理结果和设定的稳定性条件，对步骤4)建立的差分方程进行求解，得到海水和海底大地中任意时刻电磁场的分布。

实施例：如图5所示的常规海洋电磁三维构造勘探模型中存在海洋空气、海水、海底大地这三分空间，其中海洋空气、海水、海底大地的电导率分别为0s/m、3.7s/m和0.333s/m。海水深380m，在海底地下160m深处存在一个棱柱体，其长为1000m，宽为300m，高为50m。距离海底50m处的海水中设置一发射机,发射机发送高斯脉冲电流；在海底设置一接收机，用于接收电场在X、Y、Z三个方向的分量以及磁场在X、Y、Z三个方向的分量，发射机与接收机之间的距离为0m～2000m，采用本发明电偶源三维时域有限差分正演方法计算脉冲响应。

如图6所示，绘出20道二次场时序，分两组比较，分别为高阻异常与背景(异常棱柱体重填充围岩，其电阻率为0.333s/m)以及低阻异常与背景，其中图(a)表示有高阻异常时的二次场时序，时间同相轴中间向左弯，这是因为相对于海底大地，电磁波在高阻体中的传播速度大，相同距离，从发射点到接收点的传播时间短；图(b)表示有低阻异常体时的二次场时序，时间同相轴中间向右弯，这是因为相对于海底大地，电磁波在低阻体中的传播速度小，相同距离，从发射点到接收点的传播时间长。

如图7所示，3.72s时刻的归一化幅度偏移距曲线图中分别对高阻异常与背景以及低阻异常与背景进行了比较，主要用于观察当海底地下有不同3D电阻率异常体时，幅度随偏移距的变化。异常体为高阻体和低阻体时分别与背景场比较，归一化MVO(Magnitude Versus Offset幅度偏移距曲线)差别明显。

如图8所示，测得当海底地下有3D高阻异常时的不同时刻电场E_x的分布情况，图(a)～图(d)分别是1.03s、2.17s、3.72s和5.69s四个时刻全空间的波场快照，坐标单位为网格数，网格大小为45m,X方向300网格，Y方向50网格，Z方向200网格。

如图9所示，坐标单位为网格数，网格大小为45m,X方向300网格，Y方向50网格，Z方向200网格。测得当地下有3D低阻异常时不同时刻电场E_x的分布情况，图(a)～图(d)分别是1.03s、2.17s、3.72s和5.69s时全空间的波场快照，可以清楚地看到电偶极子激发的电场在海洋三分空间传播过程中遇到低阻3D异常体时的响应变化。异常反应主要集中在异常体四个顶点。

3D异常体跟2D异常体的响应有很大的不同，到后期，能量主要集中在异常体的角顶点处。这个也符合静电平衡原理，处于静电平衡的导体，内部没有自由移动的电荷，电荷主要集中在曲度大的棱角和尖端突起点上。充分说明二维模拟在构造方向介质性质相同近似有时候是不符合实际规律，而三维模拟更能反应介质真实的电磁响应规律，开展三维模拟有着非常重要意义。

本发明对硬件要求低、计算速度快且稳定性好。与有限元法和积分方程法相比，本发明可以快速的对电磁波在海洋介质中的传播规律进行正演成像，且成像结果与实际吻合。采用用高斯脉冲作为电流信号，正演模拟结果波形清晰，这对研究电磁波在海洋介质中的传播规律有很大帮助。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和方法步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (7)

1.一种电偶源三维时域有限差分正演成像方法，其包括以下步骤：

1)将高斯脉冲加载到电偶极子源上，电偶极子源的源电流I(t)为：

$I\left(t\right)=\exp (-\frac{{(t-t_0)}^2}{{\mathrm{\τ}}^2}),$

式中，t为时间，t₀为时移常数，τ为脉冲宽度因子，取τ＝0.1s；

2)对海洋空气、海水和海底大地这三分空间均建立如下麦克斯韦方程组：

$\left\{\begin{array}{c}\▿\×E(r,t)=-\frac{\∂B(r,t)}{\∂t}\\ \▿\×H(r,t)=J(r,t)\\ \▿\·B(r,t)=0\\ \▿\·J(r,t)=0\end{array}\right.,$

对海洋空气、海水和海底大地均建立如下本构方程：

$\left\{\begin{array}{c}B(r,t)=\mathrm{\μ}\left(r\right)H(r,t)\\ J(r,t)=\mathrm{\σ}\left(r\right)E(r,t)\end{array}\right.,$

式中，B(r,t)、H(r,t)、E(r,t)分别表示磁感应强度、磁场强度、电场强度，μ(r)表示介质中磁导率，σ(r)表示介质中电导率，J(r,t)表示介质中的传导电流密度，r表示测点到原点的距离， $r=\sqrt{x^2+y^2+z^2};$

3)对海洋空气、海水和海底大地这三分空间的棱柱体模型均进行均匀网格剖分；

4)假定每个剖分得到的网格上电导率和磁导率固定不变，根据步骤3)剖分得到的网格，采用时域有限差分法对海水和海底大地这两部分空间的麦克斯韦方程组均进行差分处理，得到差分方程；

5)根据步骤4)得到的海水和海底大地中的差分方程，采用解析解对海洋空气的麦克斯韦方程组进行处理，计算得到海面以上空气中的电磁场；

6)为得到步骤4)中建立的差分方程的稳定解，对剖分空间的边界条件进行处理并设定稳定性条件；

7)结合步骤6)对边界条件的处理结果和设定的稳定性条件，对步骤4)建立的差分方程进行求解，得到海水和海底大地中任意时刻电磁场的分布。

2.如权利要求1所述的一种电偶源三维时域有限差分正演成像方法，其特征在于：所述步骤1)中，高斯脉冲在时域中的脉冲宽度取6τ，其频谱宽度取4τ/5。

3.如权利要求1所述的一种电偶源三维时域有限差分正演成像方法，其特征在于：所述步骤2)中，对海洋空气、海水和海底大地棱柱体模型进行均匀网格剖分时，最大空间步长为45m。

4.如权利要求1或2或3所述的一种电偶源三维时域有限差分正演成像方法，其特征在于：所述步骤4)中，采用时域有限差分法对海水和海底大地这两部分空间的麦克斯韦方程组均进行差分处理，得到海水和海底大地中的电场强度E(r,t)的差分方程为：

$\begin{array}{c}{E}_x^{n+1}(i+1/2,j,k)=\frac{2\mathrm{\γ}-\mathrm{\Δ}{t}_n\mathrm{\σ}(i+1/2,j,k)}{2\mathrm{\γ}+\mathrm{\Δ}{t}_n\mathrm{\σ}(i+1/2,j,k)}{E}_x^n(i+1/2,j,k)\\ +\frac{2\mathrm{\Δ}{t}_n}{2\mathrm{\γ}+\mathrm{\Δ}{t}_n\mathrm{\σ}(i+1/2,j,k)}[\frac{B_z^{n+1/2}(i+1/2,j+1/2,k)-B_z^{n+1/2}(i+1/2,j-1/2,k)}{(\mathrm{\Δ}{y}_j+{\mathrm{\Δy}}_{j-1})/2}\\ -\frac{B_y^{n+1/2}(i+1/2,j,k+1/2)-B_y^{n+1/2}(i+1/2,j,k-1/2)}{(\mathrm{\Δ}{z}_j+\mathrm{\Δ}{z}_{j-1})/2}]\end{array},$

$\begin{array}{c}{E}_y^{n+1}(i,j+1/2,k)=\frac{2\mathrm{\γ}-\mathrm{\Δ}{t}_n\mathrm{\σ}(i,j+1/2,k)}{2\mathrm{\γ}+\mathrm{\Δ}{t}_n\mathrm{\σ}(i,j+1/2,k)}{E}_y^n(i,j+1/2,k)\\ +\frac{2\mathrm{\Δ}{t}_n}{2\mathrm{\γ}+\mathrm{\Δ}{t}_n\mathrm{\σ}(i,j+1/2,k)}[\frac{E_x^{n+1/2}(i,j+1/2,k+1/2)-E_x^{n+1/2}(i,j+1/2,k-1/2)}{(\mathrm{\Δ}{z}_k+\mathrm{\Δ}{z}_{k-1})/2}\\ -\frac{B_z^{n+1/2}(i+1/2,j+1/2,k)-B_z^{n+1/2}(i-1/2,j+1/2,k)}{(\mathrm{\Δ}{x}_i+\mathrm{\Δ}{x}_{i-1})/2}]\end{array},$

$\begin{array}{c}{E}_z^{n+1}(i,j,k+1/2)=\frac{2\mathrm{\γ}-\mathrm{\Δ}{t}_n\mathrm{\σ}(i,j,k+1/2)}{2\mathrm{\γ}+\mathrm{\Δ}{t}_n\mathrm{\σ}(i,j,k+1/2)}{E}_z^n(i,j,k+1/2)\\ +\frac{2\mathrm{\Δ}{t}_n}{2\mathrm{\γ}+\mathrm{\Δ}{t}_n\mathrm{\σ}(i,j,k+1/2)}[\frac{B_y^{n+1/2}(i+1/2,j,k+1/2)-B_y^{n+1/2}(i-1/2,j,k+1/2)}{(\mathrm{\Δ}{x}_i+\mathrm{\Δ}{x}_{i-1})/2}\\ -\frac{B_x^{n+1/2}(i+1/2,j,k+1/2)-B_x^{n+1/2}(i,j-1/2,k+1/2)}{(\mathrm{\Δ}{y}_j+\mathrm{\Δ}{y}_{j-1})/2}\end{array},$

式中，γ表示虚拟介电常数，Δt_n为时间步长，Δ为空间步长，μ为磁导率；E_x、E_y和E_z分别为电场强度E(r,t)在X轴、Y轴和Z轴的分量；

海水和海底大地中的磁感应强度B(r,t)的差分方程为：

$\begin{array}{c}{B}_x^{n+1/2}(i,j+1/2,k+1/2)=B_x^{n-1/2}(i,j+1/2,k+1/2)\\ -\frac{\mathrm{\Δ}{t}_{n-1}+\mathrm{\Δ}{t}_n}{2}\×[\frac{E_z^n(i,j+1,k+1/2)-E_z^n(i,j,k+1/2)}{\mathrm{\Δ}{y}_j}\\ -\frac{E_y^n(i,j+1/2,k+1)-E_y^n(i,j+1/2,k)}{\mathrm{\Δ}{z}_k}]\end{array},$

$\begin{array}{c}{B}_y^{n+1/2}(i+1/2,j,k+1/2)=B_y^{n-1/2}(i+1/2,j,k+1/2)\\ -\frac{\mathrm{\Δ}{t}_{n-1}+\mathrm{\Δ}{t}_n}{2}\×[\frac{E_x^n(i+1/2,j,k+1)-E_x^n(i+1/2,j,k)}{\mathrm{\Δ}{z}_k}\\ -\frac{E_z^n(i+1,j,k+1/2)-E_z^n(i,j,k+1/2)}{\mathrm{\Δ}{x}_i}\end{array},$

$\begin{array}{c}{B}_z^{n+1/2}(i+1/2,j+1/2,k)=B_z^{n+1/2}(i+1/2,j+1/2,k+1)\\ +\mathrm{\Δ}{z}_k[\frac{B_x^{n+1/2}(i+1,j+1/2,k+1/2)-B_x^{n+1/2}(i,j+1/2,k+1/2)}{\mathrm{\Δ}{x}_i}\\ +\frac{B_y^{n+1/2}(i+1/2,j+1,k+1/2)-B_y^{n+1/2}(i+1/2,j,k+1/2)}{\mathrm{\Δ}{y}_j}\end{array},$

式中，B_x、B_y和B_z分别为磁感应强度B(r,t)在X轴、Y轴和Z轴的分量；

海水和海底大地中的磁场强度H(r,t)的差分方程为：

$\left\{\begin{array}{c}{H}_x(i,j+1/2,k+1/2)=B_x(i,j+1/2,k+1/2)/\mathrm{\μ}(i,j+1/2,k+1/2)\\ H_y(i+1/2,j,k+1/2)=B_y(i+1/2,j,k+1/2)/\mathrm{\μ}(i+1/2,j,k+1/2)\\ H_z(i+1/2,j+1/2,k)=B_z(i+1/2,j+1/2,k)/\mathrm{\μ}(i+1/2,j+1/2,k)\end{array}\right.,$

式中， $\mathrm{\μ}(i,j+1/2,k+1/2)=\frac{\mathrm{\Δ}{x}_{i-1}\mathrm{\μ}(i-1,j,k)+\mathrm{\Δ}{x}_i\mathrm{\μ}(i,j,k)}{\mathrm{\Δ}{x}_{i-1}+\mathrm{\Δ}{x}_i},$ i表示各网格点在X轴方向的位置，j表示各网格点在Y轴方向的位置，k表示各网格点在Z轴方向的位置，i＝i_min,i_min+1,…,i_max；j＝j_min,j_min+1,…,j_max；k＝k_min,k_min+1,…,k_max；H_x、H_y和H_z分别为磁场强度H(r,t)在X轴、Y轴和Z轴的分量。

5.如权利要求4所述的一种电偶源三维时域有限差分正演成像方法，其特征在于：所述步骤5)中，采用解析解计算得到海面以上空气中的电磁场，其具体包括以下步骤：

I)对海面以上的电磁场进行二维傅里叶变换，从空间域到波数域；

对磁感应强度B(r,t)在Z轴的分量B_z做二维傅里叶变换，得到波数域

将波数域乘以操作因子计算得到波数域将波数域乘以操作因子计算得到波数域

II)计算波数域中海面以上半空间步长处水平方向的磁感应强度；

由磁感应强度B(r,t)在自由空间中遵循矢量Laplace方程：▽²B＝0以及B_z的2D空间傅里叶变换得到波数域方程组：

式(1)和式(2)中，表示B_x的2D空间傅里叶变换，表示B_y的2D空间傅里叶变换；k_x表示相对于x坐标的波数域变量，k_y表示相对于y坐标的波数域变量；

利用式(1)和式(2)在海平面向上连续计算得到在自由空间中和的值：

将式(3)和式(4)分别代入式(1)和式(2)中，得到：

式(3)～式(6)中，exp表示以自然对数e为底的指数函数；

III)对空气中的磁感应强度进行二维傅里叶逆变换，从波数域到空间域；

对进行二维空间傅里叶反变换，得到空间域的B_x，对进行二维空间傅里叶反变换，得到空间域的B_y；

IV)在海面以下采用时域有限差分法计算边界处的切向电场；

通过海面上B_z及海面以上半空间步长处的B_x、B_y，用差分方程计算海面上切向电场E_x、E_y。

6.如权利要求4所述的一种电偶源三维时域有限差分正演成像方法，其特征在于：所述步骤6)中，对剖分空间的边界条件进行处理，其具体包括：

(I)对海洋空气与海水之间的边界进行处理；

利用电磁场边界条件，介质面上磁场切向连续，磁感应强度法向连续的基本规律，麦克斯韦方程组隐含着两个连续性边界条件，电磁场在穿过不连续的介质面时遵循以下边界条件：

电场和磁场的切向分量连续，即E_1ta＝E_2ta，H_1ta＝H_2ta，ta表示切向；

电通量密度和磁感应强度法向分量连续，即D_1n＝D_2n，B_1n＝B_2n，n表示法向，或者ε₁E_1n＝ε₂E_2n，μ₁H_1n＝μ₂H_2n；

(II)对剖分空间的前边界、剖分空间的后边界、剖分空间的左边界、剖分空间的右边界以及海底大地下无穷远处的边界进行处理；

剖分空间的前、后、左、右边界以及海底大地下无穷远处的边界都采用狄利克雷吸收边界，看成在无穷远处，电磁场的值衰减为0；

(III)在介质的边界面上，采用电场的切向连续和磁感应强度的法向连续对海水与海底大地之间的边界进行处理。

7.如权利要求4所述的一种电偶源三维时域有限差分正演成像方法，其特征在于：所述步骤6)中，设定的稳定性条件为：

$\mathrm{\Δ}{t}_{\max }=\mathrm{\α}{\left(\frac{{\mathrm{\μ}}_{\min }\mathrm{\σt}}{6}\right)}^{1/2}{\mathrm{\Δ}}_{\min },$

式中，系数α的取值范围是0.1～0.2，Δt_max为最大时间步长，Δ_min为最小空间步长，Δ_min取决于计算的精度要求，μ_min为最小介质磁导率，σ为介质电导率。