0 引言

在水下应用中，电磁波起到至关重要的作用，比如包括潜水员与潜水员间通信，水下目标探测，通信以及石油和天然气勘探等.研究表明当水下目标电磁特性不同时，水下电磁响应具有较大电磁差异，这一点在水下勘探中已作为重要的理论依据^[1-3].在水下电磁场的计算中，FDTD方法^[4]由于其独特的时域全波算法特点，越来越受到广大研究人员的关注.在最近的研究中表明FDTD方法适合用于水下电磁场的仿真^[5].

水是有耗的色散媒质，电磁波频率越高在水下传播衰减越快.因此，水下电磁场的远距离传播主要考虑的是低频信号.通常情况下，FDTD计算中剖分空间的网格尺寸δ与入射波长有关，即与计算区域内介质中最小波长λ有关，需满足色散条件δ ≤ λ/10.由于低频信号具有波长长的特点，当利用FDTD仿真水下电磁场时，如果利用较大网格剖分足够大的计算空间时，很难精确地描述激励源；如果采用较小网格剖分计算空间时，会造成计算目标的电尺寸不能取很大.目前解决这一问题的常用方法是非均匀网格技术^[4].通常非均匀网格FDTD算法包括亚网格技术和变网格技术两种.亚网格技术在局部采用均匀细网格，在全局仍然采用均匀粗网格.变网格步长方法在空间上采用渐变步长，配合被计算目标的结构或材料的空间变化，按需调整空间步长.但是，以上两种方法实现的粗细网格尺寸的比值较小，处理过渡区域较为复杂，易在粗细网格过渡区域产生非物理反射，特别是当粗细网格尺寸变大会造成算法不稳定.针对上述问题，文献^[6]提出了双离散空间的思路，但该文章未给出详细的过渡边界处理过程.文献^[7]提出了一种等效源FDTD方法，但给出的过渡边界处理方法只适合于水下单一媒质中电磁场的计算，并且没有考虑甚低频电磁波在高损耗媒质中的传播特性.

针对上述问题，本文提出基于总场-散射场技术的大比例网格变换的总场-散射场源FDTD方法(TSS-FDTD：Total-field Scattered-field Source FDTD).通过与已有解析解的对比，验证了该方法在分析甚低频电磁波在高损耗媒质中传播的有效性和正确性.最后，计算了二维空间水下构造模型分析算例，为海底油气水合物等高阻和金属等低阻物体勘探提供重要的参考.

1 TSS-FDTD方法原理

FDTD方法是一种求解微分Maxwell方程组的时域全波电磁场数值方法.该方法的电场E、磁场H分量在空间和时间上采用交替抽样的离散方式，即每一个E(或H)场周围都有四个H(或E)场分量环绕.应用这种离散方式将含时间变量的Maxwell旋度方程转化为一组差分方程，并在时间轴上逐步推进的求解空间电磁场。

TSS-FDTD方法包括两次独立的FDTD计算过程.如图 1所示，图 1(a)为第一次FDTD计算的近场计算空间，图 1(b)为第二次FDTD计算的全计算空间.计算中，近场空间用细网格剖分，全空间部分用粗网格剖分.近场空间和全空间部分分别由两次独立的FDTD计算处理.即，第一次FDTD计算得到近场空间的场分布，在第二次计算中利用第一次的计算结果得到全空间问题的场分布.两次FDTD计算以总场-散射场变换边界(TSS boundary)作为桥梁.如图 1所示.由于FDTD中电场和磁场是交替抽样离散的，因此变换边界由两个彼此相邻的电场边界(E surface)和磁场边界(H surface)组成.本文采用卷积完全匹配层作为吸收边界^[8].

在第一次计算中建立包围激励源的变换边界，记录变换边界上的电磁场分布，其对应的电磁流可以表示为

${J_s} = n \times {H_s} {J_{ms}} = {E_s} \times n$

(1)

其中,n代表变换面的法线方向.在第二次FDTD计算中，以电场边界为基准，设电场边界以内的区域为散射场区，电场边界以外的区域为总场区.基于等效原理并采用总场散射场技术，利用第一次FDTD计算得到的面电磁流，完成第二次FDTD计算在变换边界处加入入射平面波的过程.

下面以TM波为例对该方法详细地说明，如图 2所示.定义E、H和E′、H′分别为第一次和第二次计算的电磁场节点.第一次计算的网格尺寸为δ₁，第二次计算网格尺寸为δ₂，且设δ₂/δ₁=N.根据Courant条件，两次计算的时间步长为Δt₂/Δt₁=N.第二次计算中，定义总场-散射场边界范围为{i′₀δ₂≤x≤i′_aδ₂，j′₀δ₂≤y≤j′_bδ₂}，对应在第一次计算的范围可表示为{i₀δ₁≤x≤i_aδ₁，j₀δ₁≤y≤j_bδ₁}，即第二次计算中x,y方向棱边上分别有i′_a-i′₀+1，j′_b-j′₀+1个E′_z节点，对应的，在第一次计算棱边上分别有N(i′_a-i′₀) +1，N(j′_b-j′₀) +1个E_z节点.

考虑边界x=i′₀δ₂上的等效电磁流为

${J_s} = n \times H = - x \times ({H_x},{H_y},0) = - {H_y}z,$

(2)

${J_{ms}} = - n \times E = x \times (0,0,{E_z}) = - {E_z}y.$

(3)

由麦克斯韦旋度方程可知

$\nabla \times H\prime = \frac{{\partial D\prime }}{{\partial t}} + J\prime + {J_s},$

(4)

$\nabla \times E\prime = - \frac{{\partial B\prime }}{{\partial t}} - {J_{ms}}.$

(5)

将式(2) 、(3) 代入式(4) 、(5) 可得分量方程

$\frac{{\partial H{\prime _y}}}{{\partial x}} - \frac{{\partial H{\prime _x}}}{{\partial y}} = \varepsilon \frac{{\partial E{\prime _z}}}{{\partial t}} - {H_y},$

(6)

$ - \frac{{\partial E{\prime _z}}}{{\partial x}} = - \mu \frac{{\partial H{\prime _y}}}{{\partial t}} + {E_z}.$

(7)

将式(6) 和(7) 在空间和时间上Yee离散就可以得到总场散射场边界处加入入射源的FDTD迭代公式.下面详细地给出相关的迭代公式.

研究发现，当粗细网格比例N为偶数时的计算公式要比N为奇数时的计算公式复杂，因此在这里限于篇幅，仅考虑N为偶数的情况.考虑y=j′₀δ₂的边界，总场边界设置在电场节点处，边界上的电场E′_zⁿ⁺¹(i′δ₂,j′₀δ₂)属于总场，计算涉及两个H_x节点，其中H′_x^n+1/2(i′δ₂,j′₀δ₂+δ₂/2) 处于散射场区，其余一个磁场节点位于总场区.根据总场散射场技术，需要在磁场节点H′_x^n+1/2 (i′δ₂,j′₀δ₂+δ₂/2) 处扣除入射场值，该入射场值利用第一次计算中的等效面电磁流得到.由于第一次计算的细网格节点和第二次计算的粗网格节点在空间位置上并不完全重合，因此需要在空间和时间上同时采用插值的方式得到该点处的入射场值.通过计算整理可得，磁场节点H′_x^n+1/2(i′δ₂,j′₀δ₂+δ₂/2) 处的入射场值为

$\begin{array}{l} H_{ - in}^{'n + 1/2}(i\prime {\delta _2},j{\prime _0}{\delta _2} + {\delta _2}/2) = \\ \frac{1}{4}[H_x^{N\left( {n + 1/2} \right) - 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N + 1} \right){\delta _1}/2) + \\ H_x^{N\left( {n + 1/2} \right) + 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N + 1} \right){\delta _1}/2) + \\ H_x^{N\left( {n + 1/2} \right) - 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N - 1} \right){\delta _1}/2) + \\ H_x^{N\left( {n + 1/2} \right) + 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N - 1} \right){\delta _1}/2)], \end{array}$

(8)

将上式代入式(6) ，得到电场边界上E′_zⁿ⁺¹(i′δ₂,j′₀δ₂)的迭代方程为

$\begin{array}{l} E_z^{'n + 1}(i\prime {\delta _2},{\rm{ }}j{\prime _0}{\delta _2}) = \\ CA\left( m \right)E_z^{'n}(i\prime {\delta _2},{\rm{ }}j{\prime _0}{\delta _2}) + CB\left( m \right)[\nabla \times H\prime ]_z^{n + 1/2} - \\ \frac{{CB\left( m \right)}}{{4{\delta _2}}}[H_x^{N\left( {n + 1/2} \right) - 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N + 1} \right){\delta _1}/2) + \\ H_x^{N\left( {n + 1/2} \right) + 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N + 1} \right){\delta _1}/2) + \\ H_x^{N\left( {n + 1/2} \right) - 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N - 1} \right){\delta _1}/2) + \\ H_x^{N\left( {n + 1/2} \right) + 1/2}(i{\delta _1},{\rm{ }}{j_0}{\delta _1} + \left( {N - 1} \right){\delta _1}/2)], \end{array}$

(9)

式中[▽×H′]_z^n+1/2为常规FDTD的离散公式，即

$\begin{array}{l} \left[ {\nabla \times H\prime } \right]_z^{n + 1/2} = \frac{1}{{{\delta _2}}}\left[ {H_y^{'n + 1/2}(i{\prime _0}{\delta _2} + {\delta _2}/2,{\rm{ }}j\prime {\delta _2}} \right) - \\ H_y^{'n + 1/2}\left( {i{\prime _0}{\delta _2} - {\delta _2}/2,j\prime {\delta _2}} \right) - \\ H_x^{'n + 1/2}(i{\prime _0}{\delta _2},j\prime {\delta _2} + {\delta _2}/2) + \\ H_x^{'n + 1/2}(i{\prime _0}{\delta _2},j\prime {\delta _2} - {\delta _2}/2)], \end{array}$

(10)

(10) 式中系数CA(m)和CB(m)为迭代公式的系数，如CB(m)为

$CB\left( m \right) = \left[ {\frac{{\Delta t}}{{\varepsilon \left( m \right)}}} \right]/{\left[ {1 + \frac{{\sigma \left( m \right)\Delta t}}{{2\varepsilon \left( m \right)}}} \right]_{m = i\prime dx,j\prime dy}}$

(11)

同理，对于磁场边界y=j′_bδ₂－1/2，计算位于散射场区的H′_x^n+1/2(i′δ₂,j′_bδ₂－δ₂/2) 节点涉及两个电场节点，其中电场节点E′_z(i′δ₂,j′_bδ₂)位于总场区，另外一个电场节点位于散射场区.根据总场散射场技术，为了保证等号两端场量属性的一致性，应在电场节点E′_z(i′δ₂,j′_bδ₂)处减去入射波的场值.如果需要n时刻电场节点E′_z(i′δ₂,j′_bδ₂)处的入射波场值，则需要记录第一次FDTD计算Nn时刻电场节点E_z(iδ₁,j_bδ₁)处的场值.因此，磁场边界y=j′_bδ₂－1/2处H′_x的迭代方程为

$\begin{array}{l} H_x^{'n + 1/2}x(i\prime {\delta _2},{\rm{ }}j{\prime _b}{\delta _2} - {\delta _2}/2) = \\ CP\left( m \right)H_x^{'n - 1/2}(i\prime {\delta _2},{\rm{ }}j{\prime _b}{\delta _2} - {\delta _2}/2) + CQ\left( m \right)\left[ {\Delta \times E\prime } \right]_y^N + \\ CQ\left( m \right)E_z^{Nn}(i{\delta _1},{\rm{ }}{j_b}{\delta _1})/\delta 2, \end{array}$

(12)

上式中CP(m)和CQ(m)为迭代系数.(12) 式中[▽×E′]_yⁿ的计算公式为

$\begin{array}{l} \left[ {\nabla \times E\prime } \right]_y^n = \\ \left[ {E_z^{'n}(i{\prime _0}{\delta _2},{\rm{ }}j\prime {\delta _2}) - E_z^{'n}\left( {i{\prime _0}{\delta _2} + {\delta _2},{\rm{ }}j\prime {\delta _2}} \right)/{\delta _2}} \right]. \end{array}$

(13)

这样，在第二次FDTD计算的交换边界上引入了入射波的场值.需要说明的是在式(8) ～(13) 中两次FDTD计算的空间坐标具有以下关系

$i\prime {\delta _2} = i{\delta _1},j\prime {\delta _2} = j{\delta _1}$

(14)

对于变换边界的其余的三条边，按照相似的步骤进行处理，在这里不一一给出迭代式.

2 模型验证

本算例通过对比解析解和TSS-FDTD的仿真结果验证方法的正确性.数值算例的计算空间如图 3所示.假设求解空间充满均匀分布的水媒质，其电磁参数为ε_r=80，σ=1 S·m^-1.如图 3(a)所示为近场计算空间，图 3(b)所示为全计算空间.在图中虚线表示的为边长10 m的总场-散射场边界.源为线电流脉冲源，位于坐标系原点处.计算参数设定如表 1所示.脉冲源的表达式为

$I\left( t \right) = \exp \left[ { - {{\left( {\frac{{t - T}}{T}} \right)}^2}} \right],T = 600\mu {\rm{s}}.$

(15)

图 4给出了图 3所示观测线上电场的分布.其中图 4(a)给出在不同工作频率时，电场振幅随空间的分布情况.为了验证，在图中也给出了文献^[9]中解析解的计算结果.通过对比图 4中各组曲线可以看出TSS-FDTD(N=10) 计算结果与解析解吻合很好.以解析解的结果为基准，图 4(b)给出了本算例的绝对误差分析.图上显示TSS-FDTD(N=10) 方法的误差很小，并且在所关心的距离范围内TSS-FDTD(N=10) 的绝对误差最大仅为13μV·m^-1，符合工程电磁应用的要求.

本算例在PC HP Pro3380 CPU2.4GHz MEM2.0GB平台运行.表 1也列出了常规FDTD和TSS-FDTD(N=10) 中运行时间和所需内存的对比.通过对比表中的数据可以得到，TSS-FDTD方法比常规FDTD方法所需要的运行时间和内存都降低了，即消耗的运行时间和所需内存分别降低为原来的17.8%和3.9%.TSS-FDTD方法的优势很明显，计算效率得到了提高.

为了进一步地说明TSS-FDTD方法的正确性，提取了频率f=50 Hz的电场幅值的空间分布图，如图 5所示.TSS-FDTD计算中第二次计算的网格尺寸与第一次之比为N=10.图 5(a)给出了第一次计算的幅值分布图，计算空间为10 m×10 m，图中的虚线框为TSS边界.图 5(b)给出了第二次计算的相位分布图，计算空间为60 m×60 m，图中的阴影区域为TSS边界区域.通过对比第一次计算的幅值图和第二次计算的幅值图，可以发现第一次计算的电磁波通过总场-散射场边界的变换，很好地传播到第二次计算的区域.

3 探测水下目标

如图 6所示，模型空间包括三个部分，分别是空气、水和岩层.水的深度为400 m.z方向极化的激励源脉宽为0.2 ms，放置于在y方向，距离水层的上下两个界面200 m处.激励源的位置设为坐标原点.在岩层上表面上方20 m的地方设置12个接收点.接收点的位置沿x方向等间隔均匀分布，间隔为100 m，第一个接收点位于激励源的正下方，接收点以三角符号“Δ”所示.异常体位于岩层中，其尺寸为200 m×1 000 m，距离岩层上表面的距离为20 m.变换边界是以激励源为中心，边长为100 m的正方形闭合曲面，如图 6中虚线所示.应用本文方法计算异常体的电磁响应.本算例中水的本构参数为ε_r=80和σ=0.018 S·m^-1；岩层的本构参数为ε_r=6和σ=0.33 S·m^-1.当岩层中的异常体为高阻体时，其本构参数为ε_r=1和σ=0 S·m^-1；当异常体为低阻体时，其本构参数为ε_r=4和σ=100 S·m^-1.

本算例计算了岩层中异常体为高阻体、低阻体和岩层等三种情况下，12个观测点处接收到的电场强度.表 3给出了常规FDTD和TSS-FDTD(N=10) 两种情况下，计算内存和运行时间的对比情况.从表 3中可以看出，大比例网格变换的TSS-FDTD方法可以节约大量的计算时间和运行内存.计算时间和运行内存的消耗分别减少到原来的0.4%和3.54%，计算空间越大该方法的优势越明显.需要说明的是，水中激励源的脉宽很宽，并且水是一种有耗介质，这就造成如果想得到整个脉冲信号在观测点的场值，就需要较多的时间步.如表 3所示，当使用常规FDTD方法计算整个空间的运行时间将达到734 551 s之多；但如果采用本文方法取N=10，计算整个空间的总的时间仅需要2 927 s，本文方法的优势很明显.

当岩层中异常体为高阻体、低阻体和岩层三种情况时，12个接收点处归一化电场的时序波形图如图 7所示.图 7中虚线、点线和实线分别表示的异常体为高阻体、低阻体和岩层时的电场时序图.通过该图可以看出，当水下异常体的电磁参数发生变化时，观测点处电场的时序图有明显的差异.观测点距离激励源越远，接收到的时域波形脉宽越宽，这是由于水和岩层的色散特性引起的.当异常体为高阻体时，接收的电场峰值最大；其次是岩层情况；最后的是低阻体情况.对比这三种情况可以得出，根据观测点处电场强度和时序图变化情况可以判断异常体是低阻体还是高阻体.为了进一步分析不同异常体的电磁响应情况，下面以异常体为岩层情况的电场时域值作为基准值，分别对异常体为高阻体或低阻体情况的电场时域值进行差值归一，处理后的结果见图 8.图中虚线表示的是异常体为低阻体时的差值归一化结果；实线为异常体是高阻体时的差值归一化结果.

图 8可以看出，当异常体为高阻体时其反射场的值要大于为岩层时的反射场值，差值大于零；当异常体为低阻体时其反射场值小于岩层时的反射场值，差值小于零.这是由于，在激励源的频率范围内，高阻体的本质阻抗要大于岩层的本质阻抗，所以得到图 7中“高阻体—岩层”情况下的反射系数大于零，相位为正；同样，低阻体的本质阻抗要小于岩层的本质阻抗，得到“低阻体—岩层”情况下的反射系数小于零，相位为负.从反射波的波前到达观测点的顺序来看，“高阻体—岩层”情况下反射波的波前要早于“低阻体—岩层”情况.这是由于，电磁波在电导率低的媒质(高阻体)中传播速度要比在电导率高的媒质(低阻体)中传播速度要快，这是由电磁波的传播规律决定的.

4 结束语

提出一种基于等效原理的总场-散射场源大比例网格变换FDTD方法.给出该方法的理论推导及数值结果.利用总场-散射场源实现细-粗网格变换，在近场可以精确地描述源，同时也可模拟远场空间的电磁波传播.通过数值算例可以看出该方法在提高效率情况下仍然可以保持计算的精度.最后，给出水下构造模型电磁场问题的数值仿真，结果看出当水下大地以下存在异常的低阻体或高阻体时，在接收点处时序波形的幅度和到达时间都不一样.通过时序和幅度的差别可以分辨出水下大地中异常体的大致电磁特性.该方法可以为水下低频电磁计算提供有效的方法支持.