0 引言

全波形反演技术是当今地震勘探领域的研究热点^[1-9]。然而，目前该技术的研究和应用主要针对以声波方程为基础的纵波地震数据。众所周知，实际的地下介质是弹性甚至是黏弹性的，相比于传统的声波勘探，弹性波勘探可获得更丰富的矢量波场信息，可为油气勘探提供丰富的构造和岩性成像信息，是油气地震勘探技术的重要发展方向。在此背景下，研究弹性波全波形反演对于丰富弹性波地震处理技术、推动弹性波勘探技术的发展具有重要的意义。

虽然弹性波全波形反演的理论基础与声波全波形反演类似，但其计算量却大许多。一方面，求解弹性波波动方程的计算成本和存储成本比声波方程高数倍；另一方面，由于横波速度相对较低，为避免数值频散需要采用更小的计算网格，这无形中又大幅增加了计算量。因此，即便在计算机硬件高度发达的今天，庞大的计算成本依然是3D弹性波全波形反演技术面临的严峻挑战，也是该技术一直未能广泛应用的主要原因之一。针对全波形反演技术的计算效率不高的问题，现有的方法主要有两类：一类是借助相位编码等手段^[10-13]；另一类则是借助高性能计算设备及优化算法^[14-16]。前者主要通过减少同时反演炮集的个数提高计算效率，其主要限制是需要假设固定的观测系统以及会引入不同炮之间的串扰噪声，因此目前仅在海底电缆(Ocean Bottom Cable, OBC)及海底节点(Ocean Bottom Node, OBN)等符合要求的观测数据上取得较好的效果，对于地面观测数据则应用较少；后者则是一种更通用的加速方式，需要借助CPU、GPU计算集群等具有多核或重核的计算设备，借助MPI、CUDA等并行编译环境实现。3D弹性波全波形反演是一种时间成本与空间复杂度均很大的运算过程，需要同时解决计算时间过长和内存占用过大的问题。在实际应用中需要根据实际设备和具体目标优化算法。

本文针对具有分布式内存的高性能计算集群(PC-Cluster)，基于MPI与OpenMP并行编译环境下发展了基于区域分解的3D弹性波全波形反演并行算法。即结合MPI的粗粒度并行与OpenMP的细粒度并行，在充分利用PC-Cluster计算能力的前提下降低内存的占用空间，提高了全波形反演算法在大尺度模型下的计算可行性。本文首先简要介绍了3D弹性波全波形反演的理论基础，然后对基于区域分解的并行加速算法进行分析，最后通过模型试算验证方法的正确性。

2 3D弹性波全波形反演 2.1 3D波动方程数值模拟方法

3D速度—应力形式的波动方程可表示为

$ \left\{\begin{array}{l} \frac{\partial \boldsymbol{u}}{\partial t}=\boldsymbol{A} \boldsymbol{u}+\boldsymbol{s} \\ \boldsymbol{u}=\left(\begin{array}{lllllllll} v_x & v_y & v_z & \tau_{x x} & \tau_{y y} & \tau_{z z} & \tau_{x y} & \tau_{x z} & \tau_{y z} \end{array}\right)^{\mathrm{T}} \\ \boldsymbol{s}=\left(\begin{array}{lllllllll} s_{f x} & s_{f y} & s_{f z} & s_{\tau_{x x}} & s_{\tau_{y y}} & s_{\tau_{z z}} & s_{\tau_{x y}} & s_{\tau_{x z}} & s_{\tau_{y z}} \end{array}\right)^{\mathrm{T}} \end{array}\right. $

(1)

式中：u为矢量波场；s为震源矢量；t为时间；v_x、v_y、v_z分别为位移速度在x、y、z坐标轴的分量；τ_xx、τ_yy、τ_zz为正应力分量；τ_xy、τ_xz、τ_yz为剪切应力分量；s_fx、s_fy、s_fz分别为集中力源在x、y、z坐标轴的分量；s_τxx、s_τyy、s_τzz为正应力源分量；s_τxy、s_τxz、s_τyz为剪切应力源分量；矩阵A中包含了弹性参数与空间偏导数项，其具体形式为

$\boldsymbol{A}=\left(\begin{array}{ccccccccc}0 & 0 & 0 & \frac{1}{\rho} \frac{\partial}{\partial x} & 0 & 0 & \frac{1}{\rho} \frac{\partial}{\partial y} & \frac{1}{\rho} \frac{\partial}{\partial z} & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{\rho} \frac{\partial}{\partial y} & 0 & \frac{1}{\rho} \frac{\partial}{\partial x} & 0 & \frac{1}{\rho} \frac{\partial}{\partial z} \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\rho} \frac{\partial}{\partial z} & 0 & \frac{1}{\rho} \frac{\partial}{\partial x} & \frac{1}{\rho} \frac{\partial}{\partial y} \\ (\lambda+2 \mu) \frac{\partial}{\partial x} & \lambda \frac{\partial}{\partial y} & \lambda \frac{\partial}{\partial z} & 0 & 0 & 0 & 0 & 0 & 0 \\ \lambda \frac{\partial}{\partial x} & (\lambda+2 \mu) \frac{\partial}{\partial y} & \lambda \frac{\partial}{\partial z} & 0 & 0 & 0 & 0 & 0 & 0 \\ \lambda \frac{\partial}{\partial x} & \lambda \frac{\partial}{\partial y} & (\lambda+2 \mu) \frac{\partial}{\partial z} & 0 & 0 & 0 & 0 & 0 & 0 \\ \mu \frac{\partial}{\partial y} & \mu \frac{\partial}{\partial x} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \mu \frac{\partial}{\partial z} & 0 & \mu \frac{\partial}{\partial x} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \mu \frac{\partial}{\partial z} & \mu \frac{\partial}{\partial y} & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right)$

(2)

式中：λ和μ为介质的拉梅系数；ρ为介质的密度。

通过对比分析不同方法的计算效率与精度，进而选取基于交错网格的有限差分正演方法求解上述偏微分方程。所选取的时间差分精度为2阶、空间差分精度为8阶，同时使用卷积完全匹配层吸收边界条件，自由边界条件采用的是牵引力镜像法。3D交错网格有限差分的具体形式如图 1所示，其中不同的图形代表不同的参量。为提高数值模拟精度，在半网格点上利用算术平均求取密度，利用调和平均得到拉梅系数^[17]。

2.2 全波形反演算法

弹性波全波形反演的目标函数可定义为

$ O(\boldsymbol{m})=\left\|\boldsymbol{d}_{\mathrm{cal}}-\boldsymbol{d}_{\mathrm{obs}}\right\|_2^2 $

(3)

式中：m为反演模型参数；d_cal为正演地震数据；d_obs为观测地震数据。式(3)的目的是寻求最优的模型参数m，使正演数据d_cal与观测数据d_obs尽可能匹配。若利用梯度类方法对上述反演问题进行求解，则模型更新量为

$ \delta \boldsymbol{m}=-\alpha \boldsymbol{P} \nabla_\boldsymbol{m} O(\boldsymbol{m}) $

(4)

式中：P为预处理算子；α为更新步长；$\nabla_\boldsymbol{m}$为梯度向量。

将式(1)代入式(3)，通过伴随状态法可求得目标函数相对于介质模型参数的梯度向量^[18-19]，对于密度的梯度为

$ \frac{\partial O(\boldsymbol{m})}{\partial \rho}=-\int_{\mathrm{G}} \int_{t_0}^{t_1}\left(v_x^{\mathrm{w}} \frac{\partial v_x}{\partial t}+v_y^{\mathrm{w}} \frac{\partial v_y}{\partial t}+v_z^{\mathrm{w}} \frac{\partial v_z}{\partial t}\right) \mathrm{d} t \mathrm{~d} X $

(5)

式中：G代表模型空间；t₀、t₁分别代表起始和终止计算时间；X代表空间位置；上标“w”代表伴随波场分量。目标函数对于λ的梯度为

$ \frac{\partial O(\boldsymbol{m})}{\partial \lambda}=\int_{\rm{G}} \int_{t_0}^{t_1}\left[\left(\tau_{x x}^{\mathrm{w}}+\tau_{y y}^{\mathrm{w}}+\tau_{z z}^{\mathrm{w}}\right) \frac{\frac{\partial \tau_{x x}}{\partial t}+\frac{\partial \tau_{y y}}{\partial t}+\frac{\partial \tau_{z z}}{\partial t}}{(3 \lambda+2 \mu)^2}\right] \mathrm{d} t \mathrm{~d} X $

(6)

目标函数对于μ的梯度为

$ \frac{\partial O(\boldsymbol{m})}{\partial \mu}=\int_{\mathrm{G}} \int_{t_0}^{t_1}\begin{aligned} \left\{\begin{array}{l} \int \tau_{z z}^{\mathrm{w}} \frac{\partial \tau_{x z}}{\partial t} \frac{1}{\mu^2}+\tau_{y z}^{\mathrm{w}} \frac{\partial \tau_{y z}}{\partial t} \frac{1}{\mu^2}+\tau_{x y}^{\mathrm{w}} \frac{\partial \tau_{x y}}{\partial t} \frac{1}{\mu^2}+ \\ \frac{3}{2}\left[\frac{\left(\tau_{x x}^{\mathrm{w}}+\tau_{y y}^{\mathrm{w}}+\tau_{z z}^{\mathrm{w}}\right)\left(\frac{\partial \tau_{x x}}{\partial t}+\frac{\partial \tau_{y y}}{\partial t}+\frac{\partial \tau_{z z}}{\partial t}\right)}{(3 \lambda+2 \mu)^2}+\frac{\left(\tau_{x x}^{\mathrm{w}}-\tau_{z z}^{\mathrm{w}}\right)\left(\frac{\partial \tau_{x x}}{\partial t}-\frac{\partial \tau_{z z}}{\partial t}\right)}{(2 \mu)^2}+\mathrm{d} t \mathrm{~d} X\right. \\ \left.\frac{\left(\tau_{x x}^{\mathrm{w}}-\tau_{y y}^{\mathrm{w}}\right)\left(\frac{\partial \tau_{x x}}{\partial t}-\frac{\partial \tau_{y y}}{\partial t}\right)}{(4 \mu)^2}+\frac{\left(\tau_{y y}^{\mathrm{w}}-\tau_{z z}^{\mathrm{w}}\right)\left(\frac{\partial \tau_{y y}}{\partial t}-\frac{\partial \tau_{z z}}{\partial t}\right)}{(4 \mu)^2}\right] \end{array}\right] \\ \end{aligned} $

(7)

式(5)~式(7)是在3D各向同性介质中，利用一阶速度—应力方程时，目标函数对于密度和拉梅常数的梯度，利用求导链式法则得到的不同模型参数化情况下的梯度。由于在求取梯度的过程中需要分别计算正向传播波场和波场残差的反向传播波场，因此此过程是全波形反演中计算量最大的部分，也是算法提高效率的关键。

3 基于区域分解的加速策略及并行实现 3.1 基于区域分解的正演算法

区域分解算法是处理大尺度模型较常用的一种方法^[20-21]，适用于波动方程有限差分数值求解。如图 2所示，通过将原有模型分解为多块，并将每一块利用不同的CPU进行计算，可在提高计算效率的同时降低数值模拟对内存的需求。然而，由于不同块的计算并不是完全独立的，需要在计算过程中进行块与块之间的信息传递，这在一定程度上限制了该方法的加速比。在计算能力一定的情况下，该方法能够有效降低算法对单个节点内存的消耗。

高性能计算集群单个节点上往往存在多个计算核心，且它们往往共享同一内存。此时，若仅利用MPI进行并行加速，那么单个节点内部的CPU同样需要进行MPI通信，这在一定程度上会降低计算效率。结合面向分布式内存的MPI和面向共享内存的OpenMP编程环境，对单个节点内的多个计算核心利用OpenMP进行加速，对不同节点之间则利用MPI进行通信。如此，可在最大化利用计算资源的同时减少节点内部的通信，从而提高算法的计算效率。基于有限差分的3D弹性波正演和反演均是计算密集型的算法，当模型较大时，巨大的内存占用是一个重要问题。结合区域分解算法可有效提高计算效率、降低内存占用。

在3D弹性波场的正演模拟中，利用区域分解进行加速的关键是同一时刻不同块波场的数据交换以及相关的前、后期处理工作。其中，块与块间波场的数据交换可利用MPI进行点对点通信实现；处理工作包括观测系统、模型、正演数据的分解与组合、节点的分配等。为提升算法的灵活性和稳定性，本文将所有节点进行重新分组，将同一炮正演中的所有节点组成新的通信子以实现不同节点间数据交换等操作。

为测试区域分解算法对计算效率的影响，利用空间尺寸为160×160×160的模型进行波场正演，并对比不同加速策略的加速程度，具体结果如表 1所示。其中MPI算法对每个分块只用1个计算核心进行计算，而MPI+OpenMP算法用8个计算核心参与计算。由表可知，利用区域分解算法可以大幅提高计算效率，特别是在对比仅利用MPI算法时，由于单核计算的负载更大，因而加速比更为明显。当利用MPI+OpenMP的算法时，每一块对应的计算负载要弱，此时节点之间的通信对加速比的影响更大。虽然利用MPI的(2，2，2)区域分解与利用MPI+OpenMP的(1，1，1)区域分解均利用了8个计算核心进行计算，但MPI+OpenMP算法的加速比更高，这同样是由于不同节点之间的通信延迟造成的。

3.2 基于区域分解的反演算法

由于反演算法相对于正演更加复杂，将模型分解为多块后，波场的正向和反向传播、不同节点间的数据传递、不同炮对应梯度的叠加等均影响并行算法的优化程度。上述算法的优化主要针对不同类型的处理需求，并通过对节点进行重新分组实现。

多炮情况下，针对全波形反演的区域分解策略如图 3所示。除了全局通信子MPI_COMM_WORLD外，还建立了对应当前炮的通信子MPI_COMM_shot和对应每一块的通信子MPI_COMM_block。其中，MPI_COMM_shot包含当前炮所需计算核心的总和，主要用于当前炮内部的同步、数据传递等运算操作；MPI_COMM_block主要包含模型分块后相同位置对应计算核心的总和，主要用于模型分解后数据的分发、多炮梯度的归约等运算。通过将不同类型的节点进行分组，可以有效解决区域分解后不同节点间数据传递的复杂性。反演算法的具体流程如图 4。

4 模型试算及效果分析

利用简单透射模型进行试算以验证本文方法的正确性。该模型在水平方向有100×100个网格点，在深度方向有120个网格点，在不同方向上网格间距均为10 m。在均匀背景速度模型的基础上引入5%的速度扰动，扰动体为250 m×250 m×250 m的规则立方体。图 5展示了模型的纵波波场，纵横波速度比为1.5(由于横波构造形态和纵波完全相同，在此并未展示)。由于在震源点和检波点(反向传播时)处的波场为奇异的，因此在深度方向两侧10个网格点的范围内没有加入速度扰动，在反演过程中不对其进行更新。如图 5所示，在深度为50 m处设置12个震源，检波点设置在深度为1150 m的网格点上。震源子波是主频为15 Hz的Ricker子波，时间采样间隔为1 ms，观测的三分量波场时长为1.5 s。

反演过程利用了25个节点，其中24个为计算节点。每个节点具有8个CPU计算核心，利用OpenMP进行细粒度并行计算。每2个节点利用区域分解计算同一炮数据，利用MPI进行粗粒度并行计算。选取纵波初始速度为3000 m/s，横波初始速度为2000 m/s，40次迭代之后的纵、横波速度反演结果如图 6和图 7所示。由图 6可见，反演结果与真实速度模型基本吻合，不论是扰动体的构造还是速度值均得到了很好反演。图 6上还可以看到纵波与横波的反演结果并不完全相同，横波分辨率更高，这种差异性在图 7中的2D切片中更清楚。图 7中可以明显看到，相对于横波的反演结果，纵波反演结果照明范围较窄。这是由于纵波的透射角相对于横波要大，当纵波透射到模型外部时，横波依然能够到达模型底部的接收点。图 8为目标函数的迭代收敛曲线，经40次迭代后可以将目标函数降为初始值的5%左右，从侧面印证了正演波场和记录波场良好匹配。

5 结论与认识

本文基于PC-Cluster并行计算环境，针对3D弹性波全波形反演庞大的计算量以及内存占用等问题，提出了基于区域分解的弹性波全波形反演方法。通过细粒度的OpenMP与粗粒度的MPI并行环境有机结合，达到最大化利用并行计算设备、提高计算效率的目的。并行算法主要涉及MPI通信子的自定义创建、不同节点之间的数据交换等。数值测试表明，为了使计算效率最大化，需要尽可能地利用节点的计算能力以减少节点间通信在整个计算过程中所占的比例。最后通过简单透射模型的测试，验证了本文方法的正确性及高效性。虽然本文采用较小的模型进行测试，但文中相关方法和策略可以适用于大规模的3D弹性波全波形反演。