自适应物理驱动的深度学习在地震波建模中的应用复杂地形Self-adaptive physics-driven deep learning for seismic wave modeling in complex topography

摘要

解决散射波场问题是地震学和地震工程领域的关键科学问题。近年来，物理信息神经网络（PINNs）的发展在提升地震建模和反演的灵活性与效能方面展现了巨大潜力。受自适应物理信息神经网络（SA-PINNs）的启发，我们引入了一个用于复杂地形下地震波建模的框架。相关的理论模型构建以一维（1D）波动方程为例进行了演示。通过将SA-PINNs与由谱元法（SEM）生成的稀疏初始波场数据相结合，我们对二维（2D）SH波传播进行了数值模拟，涵盖了无限/半无限域和弧形峡谷/山地地形等典型情况。针对复杂的散射波场，在SA-PINNs中引入了基于时域分解的序列学习方法，以提高网络的可扩展性和解的精度。通过将SA-PINNs方法计算的位移地震图与SEM计算的结果进行对比，验证了该方法在复杂地形下模拟波传播的准确性和可靠性。结果表明，SA-PINNs具有无网格和精细化仿真的优势，并且能够实现如自由表面和侧边界波场传播等数值模拟条件。

0. 引言

当地震波遇到复杂地形时，会产生反射、折射和散射等复杂的传播特性，从而引发地形效应，导致地面运动的放大或衰减。为了揭示地形效应，关于地表地形对地震地面运动影响的理论研究已经通过解析和数值模拟方法（Trifunac, 1972；Zhang et al., 2017）进行。过去几十年里，地震波传播的数值模拟主要采用了有限差分法（Boore et al., 1971；Virieux, 1986, 1984）、有限元法（Ichimura et al., 2009）、谱元法（Komatitsch 和 Tromp, 2002a, b）和边界元法（Liu et al., 2022），这些方法将具有连续变量的解析问题转化为使用离散变量求解的数值问题。为了在无限域条件下模拟地震波的传播特性，各种数值方法对计算模型的空间离散化以及吸收边界条件（ABC）的处理提出了较高要求。解的结果质量取决于网格离散化、构成模型的选择以及输入模式，这可能会导致复杂波形模拟的收敛性和稳定性不足（Bielak et al., 2005）。引入的ABC的准确性和稳定性在波传播数值模拟研究中一直是重要问题（Xing et al., 2021a, b；Clayton 和 Engquist, 1977）。地震波传播模拟可数学化简为偏微分方程（PDE）的初始/边值问题。近年来，使用深度神经网络（DNN）结合物理原理求解PDE的方法得到了快速发展。随着数据量和计算资源的爆炸性增长，数据驱动的机器学习在图像和自然语言处理方面取得了革命性的成果。然而，数据驱动方法对数据的高度依赖、物理解释性缺乏、容易过拟合以及可访问数据的稀缺性限制了这些方法在许多学科和工程领域的应用。科学领域的人工智能（AI）范式作为一个新兴的跨领域研究方向逐渐兴起，它将物理原理与数据相结合，构建出更为通用和解释性强的机器学习模型。实例包括Deep Ritz（E 和 Yu, 2018）、Deep Galerkin（Sirignano 和 Spiliopoulos, 2018）以及PINNs（Raissi et al., 2019）。Raissi等人（2019）提出了PINNs，将PDE及其初始/边值条件纳入到DNN的损失函数设计中，并在PDE的正问题和逆问题求解中取得了良好效果。PINNs方法的基本框架概念由Lagaris等人（1998）在1990年代提出，Raissi等人（2019）通过使用自动微分的深度学习技术实现并扩展了应用场景。

Raissi等人（2019）提出的连续PINNs算法无缝地结合了数据项和控制方程项，在揭示由已知PDE控制的物理现象时展现了极大的灵活性，因此被广泛应用于PDE的求解。PINNs方法已被应用于多个领域，如流体动力学（Sun et al., 2020；Guo et al., 2022；Rao et al., 2020；Jin et al., 2020；Raissi et al., 2020；Jagtap et al., 2022b, a；Mao et al., 2020）、地球物理学（Smith et al., 2020；Moseley et al., 2006；Rasht-Behesht et al., 2022；Song 和 Alkhalifah, 2022；bin Waheed et al., 2021）、固体力学（Shukla et al., 2020；Rao et al., 2021；Haghighat et al., 2021；Shukla et al., 2022）（见表1）。相关文献综述（Cai et al., 2021a；Cuomo et al., 2022；Karniadakis et al., 2021）提供了关于PINN方法及其应用的更详细描述。

与其他广泛研究的PDE相比，求解波动方程更具挑战性，因为它的解具有多尺度特性、波形几何扩散等。已有少量研究将PINNs应用于波动方程及其全波形反演。例如，Moseley等人（2006）设计了一个神经网络来求解复杂二维（2D）声学介质中的压力响应，仅使用了解的前几个时间步作为训练数据。Rao等人（2021）提出了一个无需标签数据的PINNs框架来模拟弹性动力学问题，并成功模拟了有限域中弹性波的传播。Rasht-Behesht等人（2022）提出了一种适用于二维声波方程的PINNs算法，通过将传感器数据添加到网络中进行训练，同时前向求解，实现了全波形反演（FWI）。然而，一些学者已将PINNs框架引入到频域中的波动方程求解中。通过解Helmholtz方程并考虑简单的背景波场，神经网络得到了频域中的散射波场（Song 和 Alkhalifah, 2022；Alkhalifah et al., 2021；Song et al., 2021；Wu et al., 2022）。

然而，PINNs方法通过重新定义解为损失函数的最小值来获取DNN的权重，这可能导致近似不准确或收敛不完全。许多研究试图改善PINNs的收敛性或求解效率，例如自适应采样方法（Wu et al., 2023；Peng et al., 2022；Das 和 Tesfamariam, 2022）、自适应激活函数（Jagtap et al., 2020a, b, 2022c；Jagtap 和 Karniadakis, 2022），以及缓解高频难以学习的问题（Wang et al., 2021, 2022；Huang 和 Alkhalifah, 2022a, b）。PINN方法的主要局限之一是当计算域较大时，DNN训练的高计算成本。为此，Jagtap等人基于基本PINNs方法提出了空间域分解（cPINNs）（Jagtap et al., 2020c）以及更通用的时空域分解（XPINNs, APINNs）（J. 和 Karniadakis, 2020；Hu et al., 2022b, a）。除了可能的并行实现（Shukla et al., 2021），这些域分解方法可以根据我们对解行为的先验知识，灵活选择每个子域的DNN和超参数。一些研究（Shin et al., 2020；Mishra 和 Molinaro, 2023；De Ryck et al., 2023）已证明并估计了由不同PDE近似引起的（扩展）PINNs的泛化误差上限。

除了这些对PINN的改进，McClenny和Braga-Neto（2022）提出的自适应物理信息神经网络（SA-PINNs）利用了一种软注意机制，使DNN能够自主学习并专注于解决难以收敛的解域，且可以针对每个配点自适应训练权重。在本文中，我们使用SA-PINNs解决一维（1D）和二维（2D）波动方程在考虑复杂地形下的无限域和半无限域问题。研究目标是探索该方法在地震工程领域广泛研究的地形效应中的应用可能性。时间域内的分解是逐步进行的，从早期有限的时间范围逐步学习到整个时间域的解，遵循时间依赖系统中的时间因果关系。讨论了SA-PINNs在二维波传播模拟中的自由表面边界条件“软约束”方法，并通过数值实验验证了该方法在弧形峡谷和山地地形中的有效性。

1. 自适应物理信息神经网络（SA-PINNs）

通过空间-时间偏微分方程（PDE）的初始/边值问题介绍SA-PINNs的设计原理。令满足以下形式的PDE：

边界条件为：

初始条件为：

其中，是PDE的解，是带参数的微分算子，是定义在域中的空间变量，是时间变量。表示一组边界算子，包括Dirichlet、Neumann、Robin或周期边界条件，为对应的边界条件，和表示初始条件函数。

许多动态系统都可以通过这种形式建模。给定初始条件和边界条件，可以通过解析或数值方法得到解。从函数逼近理论的角度看，神经网络具有单隐层时，可以准确逼近任意线性或非线性连续函数或算子（Chen 和 Chen, 1995）；因此，神经网络可以视为一种非线性函数逼近器。Raissi等人（2019）通过训练深度神经网络（DNN）来逼近PDE的解，其中表示DNN中的可训练参数。

对于上述时间相关的PDE初始/边值问题，SA-PINNs通过设计损失函数来进行训练，具体形式为：

其中，、、分别为PDE残差、边界条件和初始条件损失项，、、是可调超参数，需先验定义，、、是对应每个损失项的自适应权重。

其中，和分别表示在区域和边界处选择的时空配点的数量，时间。初始条件损失项中的配点数量为。为了确保DNN逼近的解收敛于PDE的解，必须选择足够数量的配点，本研究中采用Sobol序列算法（Sobol, 1967）进行时空采样。

SA-PINNs训练神经网络的一个关键特征是同时最小化损失函数关于DNN权重的梯度，同时最大化自适应权重。设，优化目标为：

梯度下降/上升算法更新如下：

其中，为训练DNN参数的学习率，为自适应权重训练的学习率。在本研究中，我们通过PyTorch（Paszke等人, 2019）深度学习框架实现了梯度上升/下降算法。网络参数的更新可以直接通过Adam优化器（Kingma 和 Ba, 2015）实现。自适应参数的梯度上升可以通过在每个训练步骤更新为来实现。在Adam优化器之后，自适应权重固定，网络参数进一步使用LBFGS优化器（Liu 和 Nocedal, 1989）优化，以获得高精度解。

我们以一维波动方程为例，介绍SA-PINNs的实现：

初始条件为：

边界条件为：

波速。上述波动方程的精确解可以通过分离变量法得到：

在本文中，通过对解变量和相应方程的量纲分析，确定了各损失项的可调超参数。我们还对网络超参数进行了调整，发现用于训练网络的超参数处于精度的可接受范围内。PDE残差、初始条件和边界条件损失项的配点数量分别为，，，配点的空间分布如图1所示。我们设置了，。采用了一个具有四个隐藏层且每层包含30个神经元的全连接神经网络。本文的源代码和相关数据将在发表后公开，链接为：https://github.com/norery/SA-PINN-SH2d。

DNN预测结果与解析解的相对误差定义如下：

其中，为PINNs/SA-PINNs预测的波场，为由精确解给出的参考波场。用于对比的配点数为。

使用两个Adam优化器以5e-3的学习率分别实现网络参数的梯度下降和自适应权重的梯度上升，训练步数为10,000步。PINNs/SA-PINNs预测结果与精确解的相对误差收敛过程如图1所示。训练2,000步时，PINNs/SA-PINNs的相对误差分别为0.56和0.13；训练8,000步时，分别为0.46和0.02。

图2显示了每经过2,000步Adam优化器训练后SA-PINNs在PDE损失函数项中的自适应权重的时空分布。在训练初期，自适应权重在区间[0, 1]内随机初始化。随着训练迭代的进行，权重值逐渐增加，并在时空中实现了自适应分布，使SA-PINNs能够更加关注具有较大权重的区域的网络参数优化。

图3展示了经过10,000步Adam优化器训练后，PINNs和SA-PINNs的预测解与精确解以及通过有限差分法计算的数值解的对比。在使用LBFGS优化器训练5,000步后，PINNs和SA-PINNs预测解与精确解的最大绝对误差分别为0.1和0.01。与PINNs相比，SA-PINNs能够更快逼近波动方程的解，并且解的精度提高了一个数量级，这反映了SA-PINNs更快的收敛速度和更高的逼近精度。

2. SA-PINNs中二维波传播的数值模拟方法

2.1 方法

在使用PINNs模拟二维或三维（3D）波传播时，研究了几种源项引入的方法。Song等人（2021）通过基于背景波场解散射压力波场，以消除源附近波场的奇异性和复杂性。Rao等人（2021）使用圆形波源来在有限域内施加等效的地震源。Moseley等人（2006）通过使用早期几个时刻的波场训练神经网络，并在中途加入控制方程损失项来模拟声波的传播。Rasht-Behesht等人（2022）使用SPECFEM2D计算的早期两个时刻的波场作为初始条件来模拟声波问题。

本研究基于使用谱元法（SEM）计算的两个早期位移波场作为初始条件来模拟震源。第一个波场快照限制源的位置和形状，第二个波场快照限制波传播的方向。

二维SH波传播问题的控制方程可以表示为：

其中，是位移场函数，是物理波速。令，因此通过给定初始波场来等效施加力。

用于求解二维波动方程的SA-PINNs架构如图4所示，其中表示通过DNN逼近的PDE解，即。和是分别使用SEM计算的两个早期时刻和的波场。神经网络的输入是时空坐标，输出是由DNN逼近的PDE解，即位移场。通过自动微分法计算、、，并构建损失函数。

在本研究中，选择了tanh激活函数作为非线性激活函数。DNN中的可训练参数使用Xavier初始化，波动方程归一化方法参照Rasht-Behesht等人（2022）进行。Sobol序列算法用于生成用于PDE残差损失项的时空配点。附录中总结了本研究使用的DNN超参数。

在SEM中，采用时间步长为秒的二阶显式Newmark时间步格式。使用SEM计算和时刻的位移场作为参与训练的初始条件，其余时间的结果作为真实值来比较SA-PINNs的精度。SA-PINNs中的波场时间从主初始条件开始计时。模型在空间上使用四阶Legendre谱元法离散化，在SA-PINNs中每个初始条件仅使用50×50个配点（共2500个）。在边界（不包括自由表面）外使用了完美匹配层（PML）（Komatitsch 和 Tromp, 2003）以实现边界处外向波的快速衰减。

2.2 无限域SH波传播模拟

以无限均匀介质中的波传播为例，验证该方法的可行性。计算模型如图5所示，位于一个宽度和高度均为1500米的矩形内。剪切波速设定为 m/s。输入波为主频为10 Hz的Ricker源时函数，源位置位于模型中心点米和米处。在模型的纵向方向上等间距布置了41个位移传感器，位置为米。SEM计算的总持续时间为1.0秒，初始条件时间为秒，秒。SA-PINNs训练波场的最终时间为秒。用于求解二维无限均匀介质的损失函数仅包含PDE残差和初始条件损失项，即

其中，，。使用学习率为的Adam优化器训练10,000步，然后使用LBFGS优化器继续训练20,000步以获得更高精度的结果。

如图5所示，SA-PINNs和SEM计算了在0.775秒内41个纵向等距分布传感器的位移地震图。图6展示了SA-PINNs模拟的四个波场快照与SEM的对比，所示时间均从主要初始条件起始。图6表明，SA-PINNs对波场的预测结果与谱元法一致，能够准确捕捉无限介质中地震波的传播特性。

3. SA-PINNs在复杂地形下的地震波建模

3.1 半无限域SH波传播模拟

使用PINNs进行波传播模拟时，一个关键问题是如何实现自由表面边界条件。在有限元法和谱元法中，地震波方程弱形式中包含的边界积分项被直接设为零，自由边界条件可以方便地实现。然而，高精度模拟大规模复杂地形需要大量的存储和计算。使用PINNs求解均匀介质波传播问题时，传输边界条件会自动满足，自由表面的反射通过“软约束”边界条件来约束。

计算模型如图7所示，位于一个宽1500米、高1000米的矩形内。剪切波速设定为 m/s。输入波为主频为15 Hz的Ricker源时函数，源位置位于模型中心点米和米处。在模型的纵向方向上等间距布置了41个传感器，位置为米。

SEM计算的总持续时间为0.45秒，初始条件时间为秒，秒。SA-PINNs训练波场的最终时间为0.3秒。当顶部边界为地表时，应力自由条件可以通过Neumann边界条件在损失函数中表示，此处米。因此，损失函数包括PDE残差、初始条件和边界条件损失项，即

其中，，，。使用学习率为的Adam和LBFGS优化器分别训练10,000个周期，之后使用LBFGS优化器进行微调。SA-PINNs预测的41个纵向等距分布传感器的位移地震图与SEM结果进行了比较，如图7所示。

图8展示了SA-PINNs模拟的四个波场快照与SEM的对比。将无限/半无限域缩减为截断计算域，并减少或消除人工边界处反射波的干扰，是确保波动模拟效率和准确性的关键。通常使用人工吸收层（例如PML）、多重传输公式和其他人工边界来避免波反射问题。这些方法引入了复杂的数学情境和可能的不稳定问题。图8显示，损失函数中的“软约束”自由边界损失项能够表示地表反射，并满足侧边和底部边界处外向波的传输。

3.2 弧形峡谷/山地SH波传播模拟

Raissi等人（2019）提出的基础PINNs方法训练DNN同时预测整个时空域（即在所有空间和时间点预测）。在长时间序列信号中，特征容易丢失。Wang等人（2022）通过指出连续PINNs框架可能违反时间因果性并容易收敛到错误解，揭示了隐式偏差。在传统的波动方程数值算法中，需要使用适当的时间积分方法逐步积分求解离散节点的运动方程，确保在时刻的解在时刻之前获得。而PINNs框架同时预测整个时空域的波场，不遵循这种时间优先原则。

为了解决这个问题，最近的研究（Ren等人，2022；Wight和Zhao，2021；Krishnapriyan等人，2021）通过引入连续训练策略，改进了PINNs在复杂动态系统中的学习能力。复杂峡谷或山地地形上的地震波传播会产生散射波形问题。本节介绍了一种用于时域分解的连续训练策略，使DNN能够学习散射地震波的复杂行为。例如，如果在时间域内求解PDE，PDE的解在时间段（，且）内依次学习。持续训练的时间段的可训练参数被用作DNN在时间段的初始参数，并使用LBFGS优化器进行微调。实际上，网络的顺序训练旨在实现更好的训练收敛性。在此过程中，我们保持网络规模及其他参数设置一致。在第一阶段，我们在考虑较短波传播持续时间的情况下预训练网络。这可以视为网络可训练参数的“热启动”搜索。随后，我们逐渐考虑更长的持续时间，并继承已训练的参数重新训练网络。我们发现，这种顺序训练策略明显提高了训练收敛性。唯一可能出现的问题是网络对于初始小时间间隔可能过大，随着域的扩展而变得更好。否则，较小的网络对于初始域更为合适，但对于较大域可能过小，无法捕捉细节。一种可以帮助提高对域大小变化适应性的解决方案是神经元分裂（Huang和Alkhalifah，2022a）。

弧形峡谷和山地模型如图9所示，弧形模型参数为原点（1500米，1500米），半径米；山地模型参数为原点（750米，900米），半径米，剪切波速 m/s。输入波形使用主频为15 Hz的Ricker源时函数，源位置为（1500米，500米）和（750米，500米）。传感器沿模型纵向等距分布，水平坐标分别为1000米和400米，每个位置布置41个传感器。

在一个长1500米、宽1000米的矩形域内，使用SEM计算了秒和秒的位移场，这些位移场可以直接用作山地案例中的初始条件（）。经过简单的坐标变换后，可以用作峡谷案例的初始条件。此外，SEM用于计算弧形峡谷/山地地形在秒内的位移场，以验证SA-PINNs方法的精度。模型在空间上使用四阶Legendre谱元法离散化，网格为150×150，除顶部边界外，在三侧使用厚度为五节点的完美匹配层（PML）。

在弧形峡谷案例中，时域分解使用的时间跨度为0.15秒，即在秒、秒、秒和秒生成时间配点用于损失函数优化。弧形山地案例使用更细粒度的时域分解，因为山顶的地震波能量增强，形成更复杂的散射波场。时域分解使用固定时间跨度，时间采样间隔为秒、秒、秒、秒、秒和秒。

使用SA-PINNs和SEM计算的传感器位移地震图如图9所示，表明两种方法模拟的波场之间具有高度相关性。为了研究地形效应导致的地震波放大和衰减效应，计算了弧形峡谷、山地地形和自由场不同位置的地面峰值位移，如图10a和10b所示。地表位移最大值在空间上对称分布，峡谷和山地中心及边缘外侧显示出放大效应，边缘内部显示出衰减效应。

图11和图12展示了弧形峡谷/山地地形案例中六个时刻SA-PINNs预测结果与SEM模拟波场的对比。可以看出，在采用细粒度时域分解的连续训练策略后，DNN能够捕捉到微小的散射波细节，并获得高精度的时序波场模拟结果。差异图表明误差随时间不会增加，这是良好的现象，可能是由于波场振幅随时间减小，但也表明域分解在这些例子中效果良好。

4. 结论

在本研究中，我们将波传播数值模拟的基本理论与SA-PINNs相结合，建立了包含一维和二维波动方程的模拟系统。通过与多种解析解或数值模拟结果的比较，验证了该方法在无限/半无限域以及复杂地形波传播问题中的有效性。解决正问题的能力对于PINNs至关重要，因为它为逆问题提供了基础。在此，我们展示了PINNs/SA-PINNs在复杂地形中正向建模地震波的能力，为解决地震反演问题奠定了基础。提出的方法具有以下特点：

1. 避免点源奇异性：在模拟二维地震波传播时，SA-PINNs通过结合稀疏初始波场数据，能够避免点源奇异性。SA-PINNs中引入的基于时间域分解的顺序学习方法，可应用于复杂地形波传播问题的数值模拟。

2. 自适应权重与收敛速度：与PINNs方法相比，SA-PINNs可以为每个配点获取自适应权重，并且在求解波动方程时具有快速收敛和高精度的优势。

3. 边界传输与吸收边界条件：SA-PINNs在均匀介质中模拟波传播问题时，无需添加吸收边界的“软约束”来实现边界传输。由于SA-PINNs直接约束强形式的波动方程，因此可以轻松地将连续吸收边界条件嵌入到损失函数中，从而提高训练效率和精度。