基于物理信息驱动神经网络的三维初至波旅行时计算方法

都国宁，谭军*,2,3，宋鹏,2,3，解闯，王绍文

(1.中国海洋大学海洋地球科学学院，山东青岛 266100； 2.青岛国家海洋科学与技术实验室，山东青岛 266100； 3.中国海洋大学海底科学与探测技术教育部重点实验室，山东青岛 266100)

0 引言

精确求取初至波旅行时对于实现静校正[1-2]、基尔霍夫积分偏移[3]、旅行时反演[4-9]、地震定位[10]、偏移成像[11]等具有重要意义。传统求取初至波旅行时方法采用射线追踪类算法，包括试射法[12-14]和弯曲法[15-16]，但计算效率较低。由于复杂速度模型中两点之间的射线路径存在多种可能，计算极易陷入局部收敛，为了克服射线追踪类算法的缺陷，程函方程类算法得以发展。Vidale[17]基于盒式扩张的思想提出采用有限差分法求解二维程函方程，但这种算法并不完全符合波前传播规律，在因果性方面存在缺陷，由于忽略了初至能量可能迂回传播的情况，致使该算法存在不稳定性，不能取得真正的全局最小旅行时，当网格间距较小时还会带来巨大的计算量； Qin等[18]使用波前扩张理论改进了有限差分算法的因果性；VanTrier等[19]使用迎风有限差分算子提高了算法的稳定性；Sethian[20]提出快速推进法(Fast Marching Method,FMM)，利用逆风差分格式求解局部程函方程,采用窄带延拓重建旅行时波前,利用堆选排技术保存旅行时,将最小旅行时放在堆的顶部，显著缩短了寻找极小值的时间; Qian等[21]将快速清理法(Fast Sweeping Method,FSM)用于二维程函方程求解，进一步提高了有限差分算法的计算效率，其主要思想是基于因果关系将旅行时波场传播的方向分成有限个组,对每一组分别利用Gauss-Seidel迭代方法求解非线性逆风差分格式离散化后的方程组。一些学者提出因式分解形式的程函方程[22-23]，解决了震源奇异性问题，进一步提高了有限差分算法的准确度。Vidale[24]还将有限差分算法扩展到三维，此后便成为求取三维旅行时最常用的方法; Soukina等[25]将三维有限差分算法应用于各向异性介质旅行时的求取；Hole等[26]使用三维有限差分算法求解反射波旅行时。相对于二维模型，求取三维地震旅行时对计算精度和效率要求更高，而高精度的高阶有限差分会导致巨大的计算成本[27-28]。如何在求解三维程函方程时平衡精确度与计算成本是求解旅行时亟待解决的问题。

深度学习技术的蓬勃发展为解决以上问题拓展了新途径。20世纪90年代，诸多学者提出了用神经网络求解偏微分方程的想法[29-30]。Sirignano等[31]提出使用全连接神经网络(Fully Connected Neural Network,FCNN)求解偏微分方程的无网格化方法；Tompson等[32]使用卷积神经网络提高了稀疏线性偏微分方程的求解效率，得到Navier-Stokes方程的数值解。针对传统神经网络没有考虑偏微分方程本身携带的物理信息并且缺少物理意义可解释性的缺陷，Raissi等[33]提出物理信息驱动的神经网络(Physics-informed Neural Network，PINN)。与传统深度学习网络不同，PINN在利用FCNN的函数拟合功能实现偏微分方程近似的基础上，在神经网络训练过程中加入偏微分方程和实际物理条件的约束，使求解偏微分方程的结果更具实际物理意义。近年来，PINN在地球物理领域得到了广泛的应用。Xu等[34]将PINN应用于速度模型的反演；Karimpouli等[35]使用PINN求解波动方程进行正演模拟；Song等[36]将PINN用于VTI介质频率域声波方程的求解；Waheed等[37]使用PINN求解二维程函方程并将其应用于层析成像，结果显示PINN在二维初至波旅行时的求取中表现出高于传统方法的效率和准确度。

本文基于PINN实现了三维程函方程的高效、高精度求解。模型实验结果显示，基于PINN的三维初至波旅行时计算方法相对于传统的有限差分法有更高的精度和更高的计算效率。

1 程函方程的有限差分数值解

三维形式的程函方程为

(1)

式中：T(x,y,z)为(x,y,z)处的地震波旅行时；v(x,y,z)为(x,y,z)处的速度。

为了便于表示，式(1)可写为

(2)

(3)

为求解三维程函方程，Vidale[24]提出了有限差分算法，将三维速度模型剖分为如图1所示的若干正方体网格。假设已知A点旅行时和各网格点的速度，且网格间距为h，则有

图1 三维差分网格示意图

TB=TA+hsB

(4)

(5)

(6)

式(4)～式(6)分别为一维、二维、三维有限差分算子，使用以上算子对所有网格点进行差分计算，即可求得整个计算区域的旅行时。

有限差分算法凭借较高的精度和计算效率得到了广泛的应用，但该方法存在不可避免的震源奇异性问题[28]。为此，需要将程函方程转化为因式分解形式[29-30]，即将待求解旅行时T(x,y,z)分解为两个因式，有

T(x,y,z) =T0(x,y,z)τ(x,y,z)

(7)

式中：T0(x,y,z)是被指定的已知因式；τ(x,y,z)是需要求解的旅行时未知因式。

将式(7)代入式(2)即可得到因式分解形式的程函方程

T02(x,y,z)|τ(x,y,z)|2+

τ2(x,y,z)|T0(x,y,z)|2+

2T0(x,y,z)τ(x,y,z)×

=s2(x,y,z)τ(xs,ys,zs)=1

(8)

式中τ(xs,ys,zs)为震源位置的未知因式。

T0(x,y,z)用以下解析式指定

T0(x,y,z)=s(xs,ys,zs)×

(9)

式中s(xs,ys,zs)为震源处的慢度。

因式分解形式的程函方程有效解决了震源奇异性问题，但三维有限差分算法计算精度有限，而具有更高计算精度的高阶有限差分的计算成本也随之提升[27-28]。因此，本文尝试将深度学习与旅行时计算相结合，提出兼具精度与计算效率的基于PINN的地震旅行时求取方法。

2 基于PINN求取地震旅行时

2.1 PINN简介

由通用近似算法[38]可知，一个输入层有n个神经元、输出层有m个神经元的神经网络可以用来表示任意一个多维非线性函数u：Rn→Rm，而偏微分方程的建模过程也是寻找满足约束条件的非线性函数，两者具有相通之处。得益于在深度神经网络(DNN)中广泛使用的自动微分技术，在设计神经网络的损失函数时，可以融入偏微分方程中微分形式的约束条件，以获得包含物理模型约束的神经网络，即PINN。PINN的基本结构是通过FCNN近似表示一个函数，再利用自动微分技术求出偏微分方程残差和初边值残差约束，并作为正则项添加至损失函数中，最后利用梯度下降法等优化算法获得神经网络权重参数和偏微分方程物理参数。

(10)

式中σ(·)表示激活函数。常用的激活函数是Sigmoid、Tanh和ReLu等。

自动微分的反向模式就是反向传播算法的一般化，其思路是根据计算路径从后向前计算，依次得到对每个中间变量节点的偏导数，直至自变量节点处。在每个节点处根据其后续节点计算导数值，整个过程对应于多元复合函数求导时从最外层逐步向内侧求导。自动微分不涉及差分近似误差，因此能够较准确地计算导数，从而求出程函方程残差和初边值残差约束作为损失函数，再通过合适的优化算法进行神经网络参数的迭代优化，当损失函数达到最小时即实现了程函方程的精确求解。

2.2 基于PINN的三维程函方程求解

利用神经网络近似表示函数映射关系求解三维程函方程，定义一个损失函数，使训练集中的因式分解形式程函方程的残差最小。求解过程主要包括：

(1)利用FCNN近似旅行时未知因式τ(x,y,z)；

(2)定义包含了三维程函方程的损失函数，并在随机分配的网格上进行采样，构建用于神经网络训练的数据集；

(3) 通过DNN的自动微分算法计算τ(x,y,z)相对于空间坐标的偏导数；

(4) 选择合适的优化器，通过更新网络参数来最小化损失函数。

为了便于描述，本文假设震源(xs,ys,zs)处有τ(xs,ys,zs)=1。用多层深度神经网络Rτ近似(x,y,z)位置的旅行时未知因式τ(x,y,z)，即

(11)

激活函数在神经网络参数的优化中也起重要作用。许多研究证明，局部自适应激活函数比传统的固定激活函数具有更好的学习能力[40]，因此除了在最后一层使用ReLu线性激活函数外，所有的隐藏层都使用了局部自适应反正切(locally adaptive arctangent)函数，如此可在每个神经元的激活函数中使用可扩展参数，改变激活函数的斜率以提高网络的学习效率。

使用均方差(MSE)构建损失函数

(12)

式中：I表示所有采样点的集合；NI为采样点的数量；X*表示采样点位置坐标；L为因式分解形式程函方程的残差，其表达式为

s2(x,y,z)

(13)

(14)

在PINN的优化过程中，损失函数中不同组成部分的收敛效率存在差异[41]。针对损失函数各部分的权值分配方法有两种，即自适应权重和固定权重。Waheed等[37]证明了在PINN的学习过程中，自适应权重分配方法能够有更佳的收敛效率。因此，本文采取自适应权重分配策略，即通过反向传播的梯度值调整分配给损失函数中的不同项的权重，从而提高收敛效率。

基于PINN求解三维程函方程的流程如图2所示，(x*,y*,z*)表示训练集中随机选择的各点坐标。将这些点的坐标输入到随机初始化的神经网络中进行训练，并将各点对应的速度v(x*,y*,z*)、已知旅行时参数T0(x*,y*,z*)及空间导数T0(x*,y*,z*)作为网络的已知信息，通过最小化损失函数的过程实现对神经网络参数的迭代优化。神经网络训练完成后，输入待求位置的坐标(x,y,z)，使用神经网络预测并输出未知旅行时因式再与T0(x,y,z)相乘，即可求得最终的旅行时

图2 基于PINN求解三维程函方程的流程图

(15)

3 数值模拟实验与讨论

为了研究本文所提基于PINN的三维旅行时计算方法的适用性，设计均匀速度模型、水平层状等速度梯度模型、倾斜层状等速度梯度模型及局部三维Marmousi模型,开展模拟实验，并将本文方法与三维有限差分法计算结果进行了比较。所有模型实验均在Intel(R) Core(TM) i7-4720HQ的CPU上进行，模型参数见表1。

表1 模型参数统计表

3.1 均匀模型

均匀模型中的初至波解析旅行时可以通过距离与速度之比直接求出。截取x=3 m、y=14 m、z=36 m三个剖面分别展示PINN计算的旅行时结果与有限差分法计算结果及解析值的对比。如图3所示，在三个剖面中PINN计算结果都与理论旅行时相差无几，而有限差分求解结果则在接近模型边界的个别位置相对于解析值出现了偏差。

图3 均匀模型3种方法旅行时计算结果对比(★为震源位置，下同)

图4展示了上述三个剖面应用有限差分法(图4上)与本文方法(图4下)的计算结果相对于解析旅行时的绝对误差。可以看到，基于PINN求解的旅行时计算结果绝对误差较小。

图4 均匀模型有限差分法(上)、本文方法(下)计算结果相对于解析旅行时的绝对误差

本次模型实验中，采用PINN方法训练过程共耗费约19 min，训练完成后对待求旅行时进行预测仅需3 s，而使用有限差分法进行相同操作则需要约5 min。虽然在神经网络的训练过程耗时较多，但三维旅行时的应用通常需要计算多个震源的旅行时，随着震源数量的增加，本文提出的方法会展现更明显的效率优势。

3.2 水平层状等速度梯度模型

对于具有恒定速度梯度的模型，解析旅行时的方程为[42]

(16)

式中：T′(x,y,z)是从震源(xs,ys,zs)到某个网格点(x,y,z)的解析旅行时；gx、gy、gz分别表示沿x、y、z三个方向的速度梯度分量。

创建的水平层状等速度梯度模型在x方向的剖面如图5所示，设gx=gy=0、gz=5 s-1，速度在z方向以5 m/s的梯度递增。

图5 水平层状等速度梯度模型(x=1 m剖面)

取x=10 m、y=20 m、z=35 m三个剖面展示PINN求取旅行时的结果与有限差分法计算结果及解析值的对比(图6)。可以看到，在三个剖面中PINN计算结果都与理论旅行时几乎相等，而有限差分求解结果则在震源附近的个别位置相对于解析值出现了细微偏差，并且该偏差沿同一方向影响了整个速度模型的旅行时计算结果。

图6 水平层状等速度梯度模型3种方法旅行时计算结果对比

图7为上述三个剖面的有限差分算法计算结果(图7上)和本文所用方法计算结果(图7下)相对于解析旅行时的绝对误差。由图可见，基于PINN求解的旅行时绝对误差仍然小于有限差分法计算结果。

图7 水平层状模型有限差分法(上)、本文方法(下)计算结果相对于解析旅行时的绝对误差

本次模型实验中，PINN方法训练过程共耗费约25 min，训练完成后对整个速度模型的旅行时预测耗时4 s，而有限差分法计算旅行时约耗时5 mim。本文方法在多炮计算效率方面展现出巨大潜力。

3.3 倾斜层状等速度梯度模型

为了证明本文方法的泛化性以及在多震源计算中的效率优势，建立一个倾斜层状等速度梯度模型，x、y、z三个方向的速度变化梯度分别为gx=1 s-1、gy=8 s-1、gz=5 s-1。选择如图8三个剖面中的三角形所示的48个位置作为训练集震源，生成旅行时训练集数据；以不属于训练集的位置(25 m，25 m，25 m)作为测试集的震源。

图8 倾斜层状等速度梯度模型

选取x=15 m、y=30 m、z=20 m三条剖面展示PINN法求取旅行时的结果与有限差分法计算结果以及解析值的对比。如图9所示，在三个剖面中PINN法计算结果都与理论旅行时几乎完全拟合，而有限差分求解结果在对角线方向相对于解析值出现了较明显的偏差。

图9 倾斜层状等速度梯度模型3种方法旅行时计算结果对比

图10为上述三个剖面中有限差分法计算结果(图10上)与本文方法计算结果(图10下)相对于解析旅行时的绝对误差。可见在较复杂的速度模型中，基于PINN求解的旅行时的精度仍然优于有限差分法计算结果。

图10 倾斜层状模型有限差分法(上)与本文方法(下)计算结果相对于解析旅行时的绝对误差

为了验证多震源计算的效率优势，在y=11 m、y=12 m、y=13 m三个剖面上选取123个网格点作为震源，且这三个剖面上的点均与训练集中的震源点无重合。使用本文方法和有限差分法分别计算所有震源的旅行时，图11为截取(11 m,17 m,19 m)、(12 m,18 m,21 m)、(12 m,37 m,5 m)、(13 m,28 m,16 m)四个点作为震源时，x=20 m剖面上的旅行时计算结果。在该模型试验123炮的计算中，PINN方法训练、预测过程共耗时约68 min，而有限差分法计算旅行时共耗时约458 min，PINN方法体现出了相当明显的效率优势。

图11 多震源位置3种方法旅行时计算部分计算结果对比(x=20 m剖面)

3.4 三维Marmousi模型

选取x=40 m、y=60 m、z=50 m三条剖面，将二维Marmousi模型中速度变化十分复杂的部分拼接为如图12所示的局部三维Marmousi模型。

图12 倾斜层状等速度梯度模型

图13展示了训练完成的神经网络的旅行时预测结果与有限差分法计算结果的对比，可以看到，在三个剖面中PINN计算结果都与有限差分法计算结果基本吻合。

图13 三维Marmousi模型2种方法旅行时计算结果对比

图14为上述三个剖面中有限差分法计算结果与本文所用方法计算结果之间的绝对误差。由于无法计算Marmousi模型的解析旅行时，故而无法进行精度的对比，但是通过本文方法与有限差分法所求结果的绝对误差，可以验证PINN方法即使对于复杂构造模型也具有较高的稳定性。另外，该模型试验中有限差分法计算旅行时用时约536 min，而PINN方法训练、预测过程共耗时约145 min，在复杂模型中同样表现出明显的效率优势。

图14 三维Marmousi模型中有限差分法与本文方法旅行时求解结果的绝对误差

通过数值模拟实验，证明了PINN算法具有效率上的绝对优势，且在简单速度模型中应用精度也有所提高。但对于复杂模型则无法验证所求解旅行时的准确度，只能说明PINN算法能够得到与有限差分算法相似的结果。

4 结束语

本文提出一种基于PINN深度神经网络在三维速度模型中计算初至波旅行时的方法。经不同速度模型实验，结果表明本文所提方法求解的旅行时比目前常用的有限差分法计算结果更准确或结果相似，并且适用于复杂的三维速度模型。相对于传统深度学习算法，本文构建的神经网络是基于物理模型驱动的，即在损失函数中加入了三维程函方程，使计算结果更符合物理规律。此外，在给定三维速度模型中选择一部分网格点作为震源位置进行训练学习后，即可输入任意震源位置迅速求解地震波旅行时，在多震源应用中表现出有限差分法不可比拟的效率优势。