物理驱动深度学习波动数值模拟方法及应用1)

2023-02-25王进廷李小军

力学学报 2023年1期

关键词：波场元法波动

陈苏丁毅孙浩赵密王进廷李小军,,3)

* (北京工业大学城市与工程安全减灾教育部重点实验室,北京 100124)

† (中国人民大学高瓴人工智能学院,北京 100872)

** (清华大学水沙科学与水利水电工程国家重点实验室,北京 100084)

†† (中国地震局地球物理研究所,北京 100081)

引言

近场波动既是地震学和地震工程等领域所关心的关键科学问题,也是涉及大型工程地震输入和抗震分析的重大工程问题.在过去的几十年里,波动模拟主要采用有限差分法[1-2]、有限元法[3-4]、谱元法[5-7]等,将连续变量的解析问题转化为离散变量的数值问题.特别是谱元法作为较新的数值算法,不仅具备处理复杂几何问题的灵活性,还具有伪谱法的高精度性和快速收敛性,在波动模拟中取得了广泛应用.为了模拟地震波在无限域条件下的传播特性,各类数值方法对计算模型的空间离散方法和人工边界的处理方式有很高的要求.求解结果的优劣依赖于网格的剖分,对于复杂波动问题,可能会导致巨大的计算和存储;其次,所引入的人工边界的精度和稳定性问题一直以来都是研究近场波动数值模拟的一个重要问题[8-10].近场波动在数学上可归结为偏微分方程的初、边值问题,近年来,基于数学物理方法前沿成果的近场波动数值模拟成为热点.

随着数据和计算资源的爆炸式增长,数据驱动的机器学习在图像、自然语言处理等方面取得了很多革命性的成果.高度依赖数据驱动的方法,缺乏物理可解释性、易陷入过拟合以及可获取数据的稀疏性等问题,限制了这类方法在诸多学科和工程领域中的应用.人工智能科学范式(AI for science)已成为新的交叉领域,将物理原理与数据相融合,以建立更通用、解释性更强的机器学习模型,得到了广泛的关注,比如Deep Ritz 方法[11]、深度Galerkin 方法[12],PINN 方法[13]等.其中,Raissi 等[13]提出了物理信息神经网络(physics-informed neural networks,PINN),将偏微分方程以及初边值条件融入到神经网络的损失函数设计中,在求解偏微分方程的相关正问题和反问题方面获得了很好的效果.PINN 方法的基本框架理念是由Lagaris 等[14]在20 世纪90 年代提出的,而Raissi 等[13]结合深度学习技术(如自动微分)实现并扩展了应用场景.

PINN 方法在流体动力学[15-19]、热传导[20-21]、地球物理学[22-25]、固体力学[26-29]、生物医学[30]等多个领域中获得了应用(表1).综述文献[31-34] 对PINN 方法及其应用作了较为详细的介绍.

表1 PINN 方法在不同领域中的应用Table 1 Application of PINN method in different fields

已有少数研究者将物理信息神经网络方法应用于波动方程的求解以及全波形反演.比如Moseley等[23]设计了一个物理信息神经网络求解复杂二维声学介质中的压力响应,该网络只使用解的前几个时间步作为训练数据,可模拟分层介质、Marmousi p 波速度模型中由波动方程引起的各种物理现象.Rao 等[28]提出一种不依赖于任何标记数据的物理信息神经网络来建模弹性动力学问题,成功对弹性波在有限域中的传播进行了模拟.Rasht-Behesht等[24]提出了适用于二维声波波动方程的PINN算法,通过将传感器数据加入网络进行训练,可在前向求解的同时实现全波形反演(FWI).在前人工作的基础上,本文提出了一种基于物理信息神经网络的近场波动模拟方法,探索了迁移学习在波动问题中的应用以及复杂起伏地表边界在PINN 中的实现方法.

本文从物理信息深度学习的数学理论出发,系统探讨该方法应用于近场波动数值模拟的相关问题,包括模型及损失函数构建以及优化器模式等.通过多组数值试验对所提方法应用于均质场地、空间不均匀及复杂地形场地时的有效性进行验证.

1 动力系统物理驱动深度学习原理

通过对时空偏微分方程的初、边值问题来解析物理驱动深度学习方法的设计原理.设u(x,t)满足如下形式的偏微分方程

式中,u(x,t) 为偏微分方法的解,D (·) 为带参数 λ的微分算子,x是域Ω ∈Rd(d=1,2,···,n)中的空间变量,域边界 ∂ Ω由Dirichlet 边界 ∂ ΩD和Neumann 边界 ∂ ΩN组成.t∈[0,T]为时间变量.BD(t),BN(t),I0(x),I1(x)分别代表偏微分方程初边值问题中的Dirichlet和Neumann 边界条件以及初始条件函数.

众多动力系统均可采用上述形式来建模.通过给定物理模型的初始条件 I0(x),I1(x)和边界条件BD(t),BN(t),采用解析或数值方法求解偏微分方程得到任意时、空点的解值量u(x,t).从函数逼近论的角度,具有单一隐藏层的神经网络可精确逼近任何线性/非线性连续函数或运算符[35],因此可以将神经网络看作非线性函数逼近器.Raissi 等[13]考虑通过训练全连接神经网络 N N(x,t;θ)来逼近偏微分方程的解,其中 θ表示网络中的可训练参数.深度学习现已实现强形式的自动微分技术,并集成在包括Pytorch 及Tensorflow 等开发平台中,并获取神经网络所逼近的方程解的各阶偏导数,实现将PDE 以及初、边值残差作为正则项加入损失函数中.

对于上述的时空偏微分初、边值问题,通过设计如下损失函数 L (θ)来架构并训练神经网络

式中,Lp(θ),Lbc(θ),Lic(θ)分别为PDE 残差损失项、边界条件损失项和初始条件损失项.在全域 Ω,Dirichlet边界 ∂ ΩD和Neumann 边界 ∂ ΩN选取的时空采样点个数分别为Np,Nd,Nn,其中时间t∈[0,T].参与计算初始条件损失项的时空采样点个数为Nic.使神经网络逼近的解收敛到偏微分方程真解,需选择足够数量的时空采样点,本文采用Sobol 序列算法[36]进行时空采样.

以无源项一维波动(一维波动方程)为例,介绍物理驱动深度学习的实现方法.

初始条件

边界条件

式中,波速c=1.上述波动方程解析解为

使用Sobol 序列算法分别生成PDE 残差损失项、初始条件损失项、边界条件损失项的时空采样点,Np=5000,Nd=400,Nn=400,采样点的空间分布如图1 所示.采用3 个隐藏层、每层100 个神经元的全连接神经网络.定义PINN 的预测结果与解析解之间的相对 L2范数误差为

图1 无源项一维波动问题: 不同损失项采样点时空分布Fig.1 One-dimensional fluctuation problem of passive term:spatiotemporal distribution of sampling points for different loss terms

使用Adam 优化器训练40000 步,各个损失项以及相对误差收敛过程如图2 所示.PINN 的预测解与解析解、有限差分方法计算的数值解对比如图3所示.其中,有限差分方法中空间、时间离散均采用标准的二阶中心差分格式,网格大小分别为0.01 m和0.005 s.PINN 的预测解、有限差分解与解析解之间的绝对误差最大值分别为0.008 m 和0.007 m.图3表明: 经过训练的神经网络已具备较强逼近波动方程真解的能力.

图2 无源项一维波动问题: 神经网络训练过程中不同损失项的演化Fig.2 One-dimensional fluctuation problem of passive terms: evolution of different loss terms during neural network training

图3 无源项一维波动问题: 解析解与PINN 预测的结果对比Fig.3 One-dimensional fluctuation problem with passive term: comparison between analytical solution and PINN prediction

2 物理信息深度学习波动数值模拟

2.1 方法介绍

以二维波动问题阐明物理驱动深度学习方法求解波动方程的实现过程.二维标量波方程为

式中,u(x,z,t) 为出平面波动位移,c为介质物理波速.f(x,z,t) 为外荷载,令f≡0,并通过给定初始波场等效施加外力.

在求解二维波动方程的PINN 框架中,指定网络的输入为时空坐标X=(x,z,t),输出为神经网络逼近的PDE 解,即位移场Y=u(x,z,t).采用深度学习自动微分方法,计算得到,并构建损失函数.本文采用谱元法计算得到的t1,t2时刻的两个早期位移场作为初始条件,第一个波场快照 U1约束了震源的位置和形状,第二个波场快照 U2约束波动传播方向.垂直地表的应力为零的边界条件可通过损失函数中的Neumann 边界条件 ∇u(x,z=zmax,t)=0表示.包含PDE 残差损失项 Lp(θ)、初始条件损失项 Lic(θ)、自由边界条件损失项 Lbc(θ)的损失函数L(θ)可表示为

图4 求解二维波动方程的物理信息神经网络架构图Fig.4 Physical information neural network architecture diagram for solving the 2 D wave equation

2.2 方法验证

以无限均匀介质中的内源波动问题为例,验证方法的可行性.计算模型如图5 所示,模型所在区域为宽600 m、高600 m 的矩形,介质物理波速c=400 m/s.

图5 无限均匀介质内源波动计算模型Fig.5 Calculation model of endogenous fluctuation in infinite uniform medium

使用谱元方法计算得到t1=0.18 s,t2=0.2 s 的位移场作为初始条件参与训练,其余时刻的结果作为真解对PINN 的计算结果进行精度对比.谱元法中使用二阶显式Newmark 时间步长格式,时间步长为0.1 ms,计算总时长为0.9 s.PINN 中的波场时刻从第一个初始条件的时间t1=0.18 s 算起,波场训练的最终时刻为T=0.72s.采用四阶勒让德谱单元对模型在 1 00×100的网格上进行空间离散.PINN 方法具备稀疏数据下的高精度学习能力,仅使用两个时刻的50 × 50 的波场作为初始条件进行训练,即可达到与谱元方法相关系数超过0.99 的精度.

在模型四边各5 个单元厚度的区域边界使用完美匹配层[37](PML)以模拟点源在无限介质中的传播.输入波位移时程采用主频为20 Hz 的高斯源时间函数,源位置处于模型中心x=300m,z=300m.采用Sobol 序列算法在整个域中生成用于计算PDE残差损失项的采样点和初始条件采样点,空间分布如图6 所示.采样点数Np=10 000,Nic1=2500,Nic2=2500.

图6 无限均匀介质波动问题不同损失项时空采样点分布Fig.6 Spatio-temporal sampling points distribution of different loss terms for infinite uniform medium fluctuation problem

使用5 个隐藏层,每层30 个神经元的全连接网络逼近公式(7)的解.使用学习率为6×10−3的Adam优化器和L-BFGS 优化器分别训练10000 步.优化组合增强了PINN 的全局搜索和局部调优能力,从而得到高精度的结果.模型训练后,可对任意分辨率下任意时空点的波动方程解及解的各阶偏导进行预测.定义某一时刻PINN 方法与谱元方法模拟波场之间的决定系数R2为

图7 PINN 方法与谱元方法模拟波场之间的相对 L2 范数误差与决定系数 R2Fig.7 Relative L2 norm error and determination coefficient R2 between the results of PINN method and spectral element method

图8 表明: PINN 方法对波场的预测结果与谱元的结果基本一致,可以准确捕捉到波动在无限介质中的传播特性.训练中并未对边界施加任何约束,內域的波动也可实现边界透射.

图8 无限均匀介质波动问题: 谱元法与PINN 模拟的波场快照对比Fig.8 Wave field snapshot of spectral element method and PINN simulation

构建不同初始条件进行测试,发现经过训练的神经网络具有在不同初始条件下进行泛化的能力.仅改变上例中的源位置为x=200m,z=200m,使用谱元法生成t=0.18 s,t=0.2s 的位移场作为新的初始条件输入网络参与构成损失函数.在源位置为x=300m,z=300m 训练的神经网络可直接对新的初始条件进行泛化,得到高精度预测结果.PINN 方法与谱元方法模拟波场之间的相对 L2范数误差以及决定系数R2如图9 所示.图10 给出了其中四个时刻模拟波场对比,图中的时间均为与第一个初始条件相隔的时间.

图9 PINN 方法使用迁移学习与谱元法模拟波场之间的相对 L2 范数误差与决定系数 R2Fig.9 Relative L2 norm error and determination coefficient R2 between the results of PINN method using transfer learning and spectral element method

图10 谱元法与PINN 方法使用迁移学习得到的波场快照对比Fig.10 Wave field snapshots obtained using spectral element method and PINN method using transfer learning

由图9 可知,超出训练时间范围后,相对 L2范数误差出现较快的上升段、决定系数R2出现较快的下降段.在t=1.02 s 波即将完全透出的时刻,仍与谱元方法保持0.99 以上的高相关性,由此可见波场训练最终时刻为0.72 s 的物理驱动神经网络具有在新的初始条件下进行一定时间外推的能力.模拟结果验证了迁移学习在近场波动问题中应用的可行性.

3 物理驱动深度学习波动模拟应用

3.1 空间不均匀介质波动模拟

为进一步验证PINN 方法处理复杂介质波场问题的有效性,考虑空间不均匀介质中二维波动问题.模型所在区域为宽1000 m、高1000 m 的矩形,在波速c=1000m/s 的基础上,添加四个二维高斯混合的复杂波速分布.计算模型及介质物理波速分布如图11 所示,波速分布范围为500 m/s～ 1000 m/s.

图11 空间不均匀介质物理波速分布Fig.11 Physical wave velocity distribution in spatially inhomogeneous media

使用谱元方法计算t=0.15 s,t=0.175s 的位移场作为初始条件参与训练,其余时刻的结果与PINN 方法结果进行对比.每个初始条件的采样点为50×50共2500 个采样点.谱元法中二阶显式Newmark时间步长格式的时间步长为0 μs,计算总时长为0.7 s.PINN 中的计算时间从第一个初始条件的时间t1=0.15s 算起,训练波场的最终时刻为T=0.55s.

训练由两部分构成,首先使用学习率为6×10−4的Adam 优化器训练10000 个epoch,实现全域收敛,再调用L-BFGS 优化器实现优化收敛,训练步数为40000 步.

PINN 方法与谱元方法模拟波场之间的相对 L2范数误差以及决定系数R2如图12 所示.图13 对比了不同时刻PINN 预测和谱元方法模拟的结果,图中的时间为与第一个初始条件相隔的时间.结果表明: PINN 方法可实现非均匀介质中的高精度波动模拟.

图12 PINN 方法与谱元方法模拟波场之间的相对 L2 范数误差与决定系数 R2Fig.12 Relative L2 norm error and determination coefficient R2 between the results of PINN method and spectral element method

图13 空间不均匀介质波动问题: 谱元法与PINN 模拟的波场快照Fig.13 Spatially inhomogeneous media fluctuation problem: wave field snapshots simulated by spectral element method and PINN

3.2 起伏地形波动模拟

物理驱动神经网络实现波动模拟的一个关键问题是自由边界条件的加载方法.自由地表条件理论上概念明确,有限元法、谱元法可自动满足,地表起伏增加了求解与建模的难度[38].本文采用PINN 方法开展相关工况验证.

构建左边界长度为450 m、底边界长度为1200 m计算模型,介质物理波速c=1000 m/s,如图14 所示.

图14 起伏地形内源波动问题计算模型Fig.14 Computational model for endogenous fluctuation problems in undulating terrain

使用谱元方法分别计算了无限介质中t=0.09s,t=0.1 s 的位移场作为初始条件、起伏地表下t=0.3s内的位移场作为精度验证.两种计算模型均采用四阶勒让德谱单元在 1 00×100的网格上进行空间离散,其中每个初始条件使用 5 0×50共2500 个采样点.吸收边界仍使用5 个单元厚度的完美匹配层(PML).输入波位移时程采用主频为20 Hz 的高斯源时间函数,源位置为x=250 m,z=750m.采用Sobol 序列算法在整个域中生成用于计算PDE 残差损失项和自由边界条件损失项的采样点,采样点数分别为Np=10 000,Nbc=2500.采样点的空间分布如图15 所示.

图15 起伏地形内源波动问题不同损失项时空采样点分布Fig.15 Distribution of spatial and temporal sampling points for different loss terms for endogenous fluctuation problems

使用5 个隐藏层,每层30 个神经元的全连接网络.调用学习率为6×10−3的Adam 优化器训练10000个epoch,再使用L-BFGS 优化器训练40000 步.PINN 方法与谱元方法模拟波场之间的相对 L2范数误差以及决定系数R2如图16 所示.图17 展示了三个时刻的PINN 预测和谱元方法模拟结果对比,图中的时间为与第一个初始条件相隔的时间.

图16 PINN 方法与谱元方法模拟波场之间的相对 L2 范数误差与决定系数 R2Fig.16 Relative L2 norm error and determination coefficient R2 between the results of PINN method and spectral element method

图17 起伏地表内源波动问题: 谱元法与PINN 模拟的波场快照Fig.17 Internal fluctuation problem of undulating surface: wave field snapshots of spectral element method and PINN simulation