近年来，在工程地质勘探领域，人们对于海底浅部地层结构特征成像的精度要求越来越高，这对海洋地震勘探的分辨率提出了更高的要求。基于电火花震源的超高分辨率地震勘探，可获取到宽频带（0～4 kHz）、高分辨率的地震剖面^[1]，而且操作环保、安全、使用成本低，同时也能很好地适应一些狭窄通道、海域环境复杂区域海洋地震数据的采集工作^[2]，因此，目前基于电火花震源的超高分辨率地震勘探得到越来越广泛的关注。骆迪等^[3]利用电火花震源进行浅表层天然气水合物小道距高分辨率探测，垂向精度达到了1 m；等离子体电火花震源目前也已开始应用于海洋工程地震调查，能达到0.3～1 m的分辨率^[4]；荷兰GEO公司研制的Geo-Spark1600深海多电极电火花发射阵列系统，垂直分辨率可达到30 cm^[1]。基于电火花震源的超高分辨率勘探已在我国的海洋区域调查以及工程地球物理勘探中发挥着重要的作用。

基于电火花震源的超高分辨率地震勘探在实际作业时，往往需将多个电火花震源以一定的排列方式组合成枪阵，以克服单一电火花震源能量不足、激发的地震波穿透能力差等问题。枪阵的排列方式对于枪阵子波的波形和分辨率具有重要影响。油气地震勘探中，为获得理想的空气枪枪阵子波，通常在施工前基于波动方程有限差分正演模拟技术对于不同排列方式的枪阵波场进行数值模拟，然后通过对模拟波场的深入分析，从中优选出最佳的枪阵排列方式。常规的有限差分正演模拟技术在电火花枪阵波场的数值模拟中遇到困难，其主要原因是，为保证电火花枪阵子波的超高分辨率，各电火花震源间距一般为厘米级，这意味着基于常规有限差分正演模拟技术的电火花枪阵波场数值模拟需要付出巨大的计算和内存消耗。

不同于常规有限差分数值模拟算法，变网格有限差分数值模拟技术算法通过在模型的不同区域采用不同尺度的网格实现波场模拟，可显著提高计算效率，降低内存消耗，是一种理想的多尺度波场数值模拟技术。变网格有限差分数值模拟技术最早由JASTRAM和BEHLE提出，其通过改变网格步长来实现波场的稳定传播^[5]；随后，WANG和SCHUSTER^[6]提出了插值变网格算法，其在粗细网格过渡区域进行插值计算；李胜军^[7]发展了纵向变网格步长差分算法，进一步节省了计算机内存；朱生旺等^[8]将变网格算法推广到了弹性波数值模拟中；李振春等^[9]提出一种新的交错网格数值模拟方法，提高了常规变网格算法的精度；刘春园等^[10]将变网格算法在碳酸盐储层预测中进行了应用，发现了孔洞密度对地震属性的影响规律；张慧等^[11]提出了时间和空间同时变化的双变网格算法，进一步提高了计算效率，孙成禹等^[12]对其精度进行了分析；郭念民等^[13]和李振春等^[14]将变网格应用到了逆时偏移中，较为准确得对微小构造进行成像；曲英铭等^[15]、SUN等^[16]将变网格算法应用于全波形反演领域中，拓展了变网格算法的应用范围；解闯等^[17]对变网格有限差分正演模拟的虚假反射进行了分析，得到了虚假反射与正演模拟中各种参数之间的关系；姜占东等^[18]利用变网格高阶有限差分算法，高精度地实现了鬼波模拟。

常规变网格数值模拟算法在数值模拟过程中尚存在如下问题：如果粗细网格比的比例过大，就会产生严重的虚假反射，严重影响模拟精度^[17]。因此，常规变网格数值模拟时，为保证模拟精度其粗细网格比不能太大。为进一步提升电火花枪阵波场的数值模拟效率，降低内存消耗，本文在电火花枪阵数值模拟中提出了多重网格策略，即设置多层不同比例的网格，各层网格之间采用小网格比进行过渡，最终实现了小间隔电火花枪阵波场的高效高精度模拟，可为实际基于电火花枪阵的超高分辨率地震勘探提供重要的技术支撑。

三维各向同性介质中的一介速度-应力声波方程组为：

式中：$ K = \rho {v^2} $，$\rho $为地质体密度，kg/m³；

$v$为纵波速度，$ {v}_{x}{\text{、}}{v}_{y}{\text{、}} $${v_{\textit{z}}}$分别为$ x{\text{、}}y{\text{、}}{\textit{z}} $方向上质点的振动速度，m/s；

$p$为应力，Pa；

$t$为时间，s。

利用泰勒公式，对公式(1)进行离散化并整理后，可得：

式中：${{\rm{D}}}_{x}{\text{、}}{{\rm{D}}}_{y}{\text{、}}{{\rm{D}}}_{z}$分别表示对$ x{\text{、}}y{\text{、}}{\textit{z}} $的一阶微分算子；

i、j、k分别为空间x、y、$ {\textit{z}} $方向上的采样点；

$n$为时间采样点号；

$\Delta t$为时间采样间隔。

用变量$g$来表示应力$p$或者速度$ {v}_{x}{\text{、}}{v}_{y}{\text{、}} $${v_{\textit{z}}}$，则$g$关于$x$的偏导数公式为：

式中：${a_n}$为可变交错网格差分系数；

${\Delta _{2n - 1}}$为差分算子。

由公式（3）知，差分系数共$2N$个。这是由于在交错网格中，存在整网格点$k$和半网格点$k + 1/2$这2个对称点。下面对$k$点的差分系数的求取展开介绍。

如图1所示，计算节点与对称点$k$之间的距离用$ {\Delta _i} $表示，${L_k}$表示$k$点处的网格步长，由图1可知，整网格点为对称点时${{\Delta }_i}$为：

将${\Delta _i}$代入公式（3）并整理，可得矩阵^[9]：

求解矩阵（5）即可得到可变交错网格差分系数${a_i}$。

实际数值模拟时，中心波场模拟采用式（2）计算，在边界上还需做特殊处理以消除边界反射。本文采用完全匹配层方法（PML）压制人工边界反射，详见Berenger文章^[19]，这里不再赘述。

图  1  可变交错网格计算节点示意图

Figure 1.  Computer nodes of variable and staggered grid

常规变网格有限差分数值模拟时，当数值模拟参数设置不当时会产生强虚假反射，严重影响波场模拟精度。解闯等^[17]在一维均匀介质假设前提下，推出了常规变网格数值模拟的虚假反射率函数公式：

式中：$\gamma $为虚假反射率，无量纲常数；

$c$为介质速度，m/s；

$L$为网格步长，m；

$\omega $为子波主频，Hz；

$\beta $表示网格比，无量纲常数。

由式（6）可知，虚假反射率与介质的速度、网格步长、子波主频和网格比有关。当模型、震源与网格步长一定时，虚假反射率函数只与网格比有关，其随着网格比的增大而增大。

为了降低虚假反射对波场模拟的影响，本文提出多重网格策略，即在网格剖分时，设置多层不同步长的网格，各层之间采用小网格比，波场逐层过渡（常规变网格与多重网格剖分方式如图2所示），通过多重小网格比网格剖分，可有效降低虚假反射，提高波场模拟的精度，同时其可同样达到降低网格点数、提高计算效率和降低内存消耗的目的。

图  2  网格剖分方式

Figure 2.  Mesh mode of grids

首先以二维模型的数值模拟为例测试均匀网格、常规变网格和多重网格有限差分数值模拟效果。所采用的各向同性、介质均匀的二维模型如图3所示，其横纵向长度均为100 m，密度为3 000 kg/m³，速度为4 000 m/s。利用时间二阶、空间十二阶交错网格进行正演模拟，震源采用雷克子波，如图4所示，其主频为1 000 Hz，时间间隔为0.002 ms，记录时间为8 ms，炮点位置置于模型中间。

图  3  二维模型

Figure 3.  2D model

图  4  雷克子波

Figure 4.  The Ricker wavelet

利用均匀网格对图3二维模型进行模拟时，网格步长设为0.02 m，横纵向分别需要设置5 000个网格（5000×0.02 m=100 m），均匀网格模型剖分示意图如图5a所示。利用常规变网格进行模拟时，网格比设为25，网格步长分别为0.02、0.5 m，常规变网格模型剖分示意图如图5b所示，模型中间红色线条代表细网格，网格步长为0.02 m，红色网格线横纵向均设置100条；模型外围的黑色线条代表粗网格，网格步长为0.5 m，此时横纵向黑色网格线均设置196条（100×0.02 m+196×0.5 m=100 m）。利用多重网格进行模拟时，网格比设为5，网格步长分别为0.02、0.01和0.5 m，多重网格模型剖分示意图如图5c所示，模型中间红色线条代表细网格，网格步长为0.02 m，红色网格线横纵向均设置100条；蓝色线条代表中细网格，实现第一层过渡，网格步长为0.1 m，蓝色网格线横纵向均设置100条；模型外围黑色线条代表粗网格，网格步长为0.5 m，黑色网格线横纵向均设置176条（100×0.02 m+100×0.1 m+176×0.5 m=100 m）。

图  5  模型剖分示意图

Figure 5.  Diagram of model division

8 ms时刻的均匀网格、常规变网格和多重网格数值模拟波场快照如图6所示（由于采用的网格点数不同，波的初至位置存在差异），可以看出，采用均匀网格进行正演模拟时，波场效果最好；采用常规变网格算法，会产生严重的虚假反射；而采用多重网格算法，虚假反射得到有效压制，几乎可以达到和均匀网格同样的精度。

图  6  波场快照

Figure 6.  The wavefront snapshots

上述实验的运行程序基于型号为NVIDIA GT730的GPU卡进行运算。对均匀网格、常规变网格和多重网格进行正演模拟所用的内存消耗以及运行时间如表1所示。

从表1可以看出，均匀网格数值模拟算法所占内存最大（内存占用分别为常规变网格和多重网格的289倍和176倍），运行时间最长（计算时间分别为常规变网格和多重网格的117倍和84倍），因此虽然其模拟精度最高但是其对计算资源性能要求过高，在实际生产中的应用受到限制；常规变网格算法占内存最少，运行时间最快，但其模拟精度难以满足实际生产要求；相对来讲，基于多重网格的数值模拟，其可在较小的内存和计算消耗的前提下，实现多尺度模型的高精度数值模拟，其是较为理想的电火花枪阵数值模拟算法。

本文采用如图7所示的排列方式进行电火花枪阵三维数值模拟。该排列方式中电火花阵列共设置4行，行间距为20 cm，每行40个点震源，点震源之间距离为2 cm。

图  7  电火花阵列

Figure 7.  The array of spark

为了验证本文提出的多重网格算法对于电火花枪阵模拟的可行性，建立一个长、宽、高分别为100、20和30 m的三维水平层状模型（如图8所示，红色五角星代表电火花枪阵）。模型共分为3层，蓝色区域为第1层，深度为0～10 m，速度为1 500 m/s，密度为1 500 kg/m³，灰色区域为第2层，深度为10～20 m，速度为2 000 m/s，密度为2 000 kg/m³，黑色区域为第3层，深度为20～30 m，速度为3 000 m/s，密度为3000 kg/m³。利用时间二阶，空间十二阶差分正演算子进行电火花正演模拟，采用主频为1 000 Hz、时间采样步长为0.002 ms的雷克子波，时间延续长度为40 ms。为模拟海平面反射，模型上边界不加吸收边界，模型其他边界设置为50层PML吸收边界。

图  8  三维模型示意图

Figure 8.  3D model sketch map

观测系统如图9所示，采用单船拖缆的采集方式，图片顶部黑色箭头代表航向，绿线代表拖缆，拖缆共6条，每条拖缆长度为18 m，深度为2 m。每条缆接收道数为100道，道间距为0.18 m，拖缆间距为3.6 m，航向最小偏移距为5.8 m。红色圆点代表震源位置，震源深度为2 m，炮间距为1.98 m。

图  9  观测系统示意图

Figure 9.  Geometry sketch map

在图8所示的三维模型中采用多重网格算法，由震源阵列开始从内向外，设置三重网格，每重网格间网格比例设为3。第1重网格步长为0.02 m，网格点数为50×50×50；第2重网格包围在第1重网格的外面，由第1重网格向外延拓了50个网格点，网格步长变为0.06 m；剩余外层网格设为第3重网格，网格步长为0.18 m。基于这样的多重网格剖分算法，内存消耗约为544.16 M，而若采用常规网格算法（网格步长均为0.02 m）计算，内存消耗约为111.6 G，由此可知，基于实际大模型的电火花震源波场其常规均匀网格数值模拟一般计算资源难以承受。

第1炮的xoz、xoy和yoz方向局部切片分别如图10a、10b和10c所示。由图10可知，基于多重网格数值模拟的波场快照，其波前面清晰，波形连续，未见明显的虚假反射。第1炮的地震记录如图11所示，在图11的炮集记录中，直达波、反射波和多次反射波清晰可见，也未见到明显的虚假反射。由此证明，基于多重网格的有限差分数值模拟技术能够适应实际大模型的电火花枪阵数值模拟。

图  10  三维模型局部波场快照

Figure 10.  Local wavefield snapshots of a 3D model

图  11  单炮地震记录

Figure 11.  The seismic record of one shot

本文针对电火花枪阵波场的小网格数值模拟难题，在常规变网格数值模拟思想的基础上，提出了多重网格数值模拟策略。并通过数值实验对比了均匀网格、常规变网格和多重网格算法的模拟精度、内存消耗和计算效率，最后实现了基于多重网格的电火花枪阵波场的高精度数值模拟。数值模拟实验结果表明，多重网格算法可在较小的内存和计算消耗的前提下，实现电火花阵列的高精度数值模拟，其可为实际电火花枪阵勘探的观测系统设置提供重要的技术支撑。