物理驅(qū)動(dòng)深度學(xué)習(xí)波動(dòng)數(shù)值模擬方法及應(yīng)用1)

陳 蘇 丁 毅 孫 浩 趙 密 王進(jìn)廷 李小軍,,3)

* (北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室,北京 100124)

? (中國(guó)人民大學(xué)高瓴人工智能學(xué)院,北京 100872)

** (清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100084)

?? (中國(guó)地震局地球物理研究所,北京 100081)

引言

近場(chǎng)波動(dòng)既是地震學(xué)和地震工程等領(lǐng)域所關(guān)心的關(guān)鍵科學(xué)問(wèn)題,也是涉及大型工程地震輸入和抗震分析的重大工程問(wèn)題.在過(guò)去的幾十年里,波動(dòng)模擬主要采用有限差分法[1-2]、有限元法[3-4]、譜元法[5-7]等,將連續(xù)變量的解析問(wèn)題轉(zhuǎn)化為離散變量的數(shù)值問(wèn)題.特別是譜元法作為較新的數(shù)值算法,不僅具備處理復(fù)雜幾何問(wèn)題的靈活性,還具有偽譜法的高精度性和快速收斂性,在波動(dòng)模擬中取得了廣泛應(yīng)用.為了模擬地震波在無(wú)限域條件下的傳播特性,各類數(shù)值方法對(duì)計(jì)算模型的空間離散方法和人工邊界的處理方式有很高的要求.求解結(jié)果的優(yōu)劣依賴于網(wǎng)格的剖分,對(duì)于復(fù)雜波動(dòng)問(wèn)題,可能會(huì)導(dǎo)致巨大的計(jì)算和存儲(chǔ);其次,所引入的人工邊界的精度和穩(wěn)定性問(wèn)題一直以來(lái)都是研究近場(chǎng)波動(dòng)數(shù)值模擬的一個(gè)重要問(wèn)題[8-10].近場(chǎng)波動(dòng)在數(shù)學(xué)上可歸結(jié)為偏微分方程的初、邊值問(wèn)題,近年來(lái),基于數(shù)學(xué)物理方法前沿成果的近場(chǎng)波動(dòng)數(shù)值模擬成為熱點(diǎn).

隨著數(shù)據(jù)和計(jì)算資源的爆炸式增長(zhǎng),數(shù)據(jù)驅(qū)動(dòng)的機(jī)器學(xué)習(xí)在圖像、自然語(yǔ)言處理等方面取得了很多革命性的成果.高度依賴數(shù)據(jù)驅(qū)動(dòng)的方法,缺乏物理可解釋性、易陷入過(guò)擬合以及可獲取數(shù)據(jù)的稀疏性等問(wèn)題,限制了這類方法在諸多學(xué)科和工程領(lǐng)域中的應(yīng)用.人工智能科學(xué)范式(AI for science)已成為新的交叉領(lǐng)域,將物理原理與數(shù)據(jù)相融合,以建立更通用、解釋性更強(qiáng)的機(jī)器學(xué)習(xí)模型,得到了廣泛的關(guān)注,比如Deep Ritz 方法[11]、深度Galerkin 方法[12],PINN 方法[13]等.其中,Raissi 等[13]提出了物理信息神經(jīng)網(wǎng)絡(luò)(physics-informed neural networks,PINN),將偏微分方程以及初邊值條件融入到神經(jīng)網(wǎng)絡(luò)的損失函數(shù)設(shè)計(jì)中,在求解偏微分方程的相關(guān)正問(wèn)題和反問(wèn)題方面獲得了很好的效果.PINN 方法的基本框架理念是由Lagaris 等[14]在20 世紀(jì)90 年代提出的,而Raissi 等[13]結(jié)合深度學(xué)習(xí)技術(shù)(如自動(dòng)微分)實(shí)現(xiàn)并擴(kuò)展了應(yīng)用場(chǎng)景.

PINN 方法在流體動(dòng)力學(xué)[15-19]、熱傳導(dǎo)[20-21]、地球物理學(xué)[22-25]、固體力學(xué)[26-29]、生物醫(yī)學(xué)[30]等多個(gè)領(lǐng)域中獲得了應(yīng)用(表1).綜述文獻(xiàn)[31-34] 對(duì)PINN 方法及其應(yīng)用作了較為詳細(xì)的介紹.

表1 PINN 方法在不同領(lǐng)域中的應(yīng)用Table 1 Application of PINN method in different fields

已有少數(shù)研究者將物理信息神經(jīng)網(wǎng)絡(luò)方法應(yīng)用于波動(dòng)方程的求解以及全波形反演.比如Moseley等[23]設(shè)計(jì)了一個(gè)物理信息神經(jīng)網(wǎng)絡(luò)求解復(fù)雜二維聲學(xué)介質(zhì)中的壓力響應(yīng),該網(wǎng)絡(luò)只使用解的前幾個(gè)時(shí)間步作為訓(xùn)練數(shù)據(jù),可模擬分層介質(zhì)、Marmousi p 波速度模型中由波動(dòng)方程引起的各種物理現(xiàn)象.Rao 等[28]提出一種不依賴于任何標(biāo)記數(shù)據(jù)的物理信息神經(jīng)網(wǎng)絡(luò)來(lái)建模彈性動(dòng)力學(xué)問(wèn)題,成功對(duì)彈性波在有限域中的傳播進(jìn)行了模擬.Rasht-Behesht等[24]提出了適用于二維聲波波動(dòng)方程的PINN算法,通過(guò)將傳感器數(shù)據(jù)加入網(wǎng)絡(luò)進(jìn)行訓(xùn)練,可在前向求解的同時(shí)實(shí)現(xiàn)全波形反演(FWI).在前人工作的基礎(chǔ)上,本文提出了一種基于物理信息神經(jīng)網(wǎng)絡(luò)的近場(chǎng)波動(dòng)模擬方法,探索了遷移學(xué)習(xí)在波動(dòng)問(wèn)題中的應(yīng)用以及復(fù)雜起伏地表邊界在PINN 中的實(shí)現(xiàn)方法.

本文從物理信息深度學(xué)習(xí)的數(shù)學(xué)理論出發(fā),系統(tǒng)探討該方法應(yīng)用于近場(chǎng)波動(dòng)數(shù)值模擬的相關(guān)問(wèn)題,包括模型及損失函數(shù)構(gòu)建以及優(yōu)化器模式等.通過(guò)多組數(shù)值試驗(yàn)對(duì)所提方法應(yīng)用于均質(zhì)場(chǎng)地、空間不均勻及復(fù)雜地形場(chǎng)地時(shí)的有效性進(jìn)行驗(yàn)證.

1 動(dòng)力系統(tǒng)物理驅(qū)動(dòng)深度學(xué)習(xí)原理

通過(guò)對(duì)時(shí)空偏微分方程的初、邊值問(wèn)題來(lái)解析物理驅(qū)動(dòng)深度學(xué)習(xí)方法的設(shè)計(jì)原理.設(shè)u(x,t)滿足如下形式的偏微分方程

式 中,u(x,t) 為偏 微分方法 的解,D (·) 為帶參數(shù) λ的微分算子,x是域? ∈Rd(d=1,2,···,n)中的空間變量,域邊界 ? ?由Dirichlet 邊界 ? ?D和Neumann 邊界 ? ?N組成.t∈[0,T]為 時(shí)間變量.BD(t),BN(t),I0(x),I1(x)分別代表偏微分方程初邊值問(wèn)題中的Dirichlet和Neumann 邊界條件以及初始條件函數(shù).

眾多動(dòng)力系統(tǒng)均可采用上述形式來(lái)建模.通過(guò)給定物理模型的初始條件 I0(x),I1(x)和邊界條件BD(t),BN(t),采用解析或數(shù)值方法求解偏微分方程得到任意時(shí)、空點(diǎn)的解值量u(x,t).從函數(shù)逼近論的角度,具有單一隱藏層的神經(jīng)網(wǎng)絡(luò)可精確逼近任何線性/非線性連續(xù)函數(shù)或運(yùn)算符[35],因此可以將神經(jīng)網(wǎng)絡(luò)看作非線性函數(shù)逼近器.Raissi 等[13]考慮通過(guò)訓(xùn)練全連接神經(jīng)網(wǎng)絡(luò) N N(x,t;θ)來(lái)逼近偏微分方程的解,其中 θ表示網(wǎng)絡(luò)中的可訓(xùn)練參數(shù).深度學(xué)習(xí)現(xiàn)已實(shí)現(xiàn)強(qiáng)形式的自動(dòng)微分技術(shù),并集成在包括Pytorch 及Tensorflow 等開(kāi)發(fā)平臺(tái)中,并獲取神經(jīng)網(wǎng)絡(luò)所逼近的方程解的各階偏導(dǎo)數(shù),實(shí)現(xiàn)將PDE 以及初、邊值殘差作為正則項(xiàng)加入損失函數(shù)中.

對(duì)于上述的時(shí)空偏微分初、邊值問(wèn)題,通過(guò)設(shè)計(jì)如下?lián)p失函數(shù) L (θ)來(lái)架構(gòu)并訓(xùn)練神經(jīng)網(wǎng)絡(luò)

式中,Lp(θ),Lbc(θ),Lic(θ)分別為PDE 殘差損失項(xiàng)、邊界條件損失項(xiàng)和初始條件損失項(xiàng).在全域 ?,Dirichlet邊界 ? ?D和Neumann 邊界 ? ?N選取的時(shí)空采樣點(diǎn)個(gè)數(shù)分別為Np,Nd,Nn,其中時(shí)間t∈[0,T].參與計(jì)算初始條件損失項(xiàng)的時(shí)空采樣點(diǎn)個(gè)數(shù)為Nic.使神經(jīng)網(wǎng)絡(luò)逼近的解收斂到偏微分方程真解,需選擇足夠數(shù)量的時(shí)空采樣點(diǎn),本文采用Sobol 序列算法[36]進(jìn)行時(shí)空采樣.

以無(wú)源項(xiàng)一維波動(dòng)(一維波動(dòng)方程)為例,介紹物理驅(qū)動(dòng)深度學(xué)習(xí)的實(shí)現(xiàn)方法.

初始條件

邊界條件

式中,波速c=1.上述波動(dòng)方程解析解為

使用Sobol 序列算法分別生成PDE 殘差損失項(xiàng)、初始條件損失項(xiàng)、邊界條件損失項(xiàng)的時(shí)空采樣點(diǎn),Np=5000,Nd=400,Nn=400,采樣點(diǎn)的空間分布如圖1 所示.采用3 個(gè)隱藏層、每層100 個(gè)神經(jīng)元的全連接神經(jīng)網(wǎng)絡(luò).定義PINN 的預(yù)測(cè)結(jié)果與解析解之間的相對(duì) L2范數(shù)誤差為

圖1 無(wú)源項(xiàng)一維波動(dòng)問(wèn)題: 不同損失項(xiàng)采樣點(diǎn)時(shí)空分布Fig.1 One-dimensional fluctuation problem of passive term:spatiotemporal distribution of sampling points for different loss terms

使用Adam 優(yōu)化器訓(xùn)練40000 步,各個(gè)損失項(xiàng)以及相對(duì)誤差收斂過(guò)程如圖2 所示.PINN 的預(yù)測(cè)解與解析解、有限差分方法計(jì)算的數(shù)值解對(duì)比如圖3所示.其中,有限差分方法中空間、時(shí)間離散均采用標(biāo)準(zhǔn)的二階中心差分格式,網(wǎng)格大小分別為0.01 m和0.005 s.PINN 的預(yù)測(cè)解、有限差分解與解析解之間的絕對(duì)誤差最大值分別為0.008 m 和0.007 m.圖3表明: 經(jīng)過(guò)訓(xùn)練的神經(jīng)網(wǎng)絡(luò)已具備較強(qiáng)逼近波動(dòng)方程真解的能力.

圖2 無(wú)源項(xiàng)一維波動(dòng)問(wèn)題: 神經(jīng)網(wǎng)絡(luò)訓(xùn)練過(guò)程中不同損失項(xiàng)的演化Fig.2 One-dimensional fluctuation problem of passive terms: evolution of different loss terms during neural network training

圖3 無(wú)源項(xiàng)一維波動(dòng)問(wèn)題: 解析解與PINN 預(yù)測(cè)的結(jié)果對(duì)比Fig.3 One-dimensional fluctuation problem with passive term: comparison between analytical solution and PINN prediction

2 物理信息深度學(xué)習(xí)波動(dòng)數(shù)值模擬

2.1 方法介紹

以二維波動(dòng)問(wèn)題闡明物理驅(qū)動(dòng)深度學(xué)習(xí)方法求解波動(dòng)方程的實(shí)現(xiàn)過(guò)程.二維標(biāo)量波方程為

式中,u(x,z,t) 為出平面波動(dòng)位移,c為介質(zhì)物理波速.f(x,z,t) 為外荷載,令f≡0,并通過(guò)給定初始波場(chǎng)等效施加外力.

在求解二維波動(dòng)方程的PINN 框架中,指定網(wǎng)絡(luò)的輸入為時(shí)空坐標(biāo)X=(x,z,t),輸出為神經(jīng)網(wǎng)絡(luò)逼近的PDE 解,即位移場(chǎng)Y=u(x,z,t).采用深度學(xué)習(xí)自動(dòng)微分方法,計(jì)算得到,并構(gòu)建損失函數(shù).本文采用譜元法計(jì)算得到的t1,t2時(shí)刻的兩個(gè)早期位移場(chǎng)作為初始條件,第一個(gè)波場(chǎng)快照 U1約束了震源的位置和形狀,第二個(gè)波場(chǎng)快照 U2約束波動(dòng)傳播方向.垂直地表的應(yīng)力為零的邊界條件可通過(guò)損失函數(shù)中的Neumann 邊界條件 ?u(x,z=zmax,t)=0表示.包含PDE 殘差損失項(xiàng) Lp(θ)、初始條件損失項(xiàng) Lic(θ)、自由邊界條件損失項(xiàng) Lbc(θ)的損失函數(shù)L(θ)可表示為

圖4 求解二維波動(dòng)方程的物理信息神經(jīng)網(wǎng)絡(luò)架構(gòu)圖Fig.4 Physical information neural network architecture diagram for solving the 2 D wave equation

2.2 方法驗(yàn)證

以無(wú)限均勻介質(zhì)中的內(nèi)源波動(dòng)問(wèn)題為例,驗(yàn)證方法的可行性.計(jì)算模型如圖5 所示,模型所在區(qū)域?yàn)閷?00 m、高600 m 的矩形,介質(zhì)物理波速c=400 m/s.

圖5 無(wú)限均勻介質(zhì)內(nèi)源波動(dòng)計(jì)算模型Fig.5 Calculation model of endogenous fluctuation in infinite uniform medium

使用譜元方法計(jì)算得到t1=0.18 s,t2=0.2 s 的位移場(chǎng)作為初始條件參與訓(xùn)練,其余時(shí)刻的結(jié)果作為真解對(duì)PINN 的計(jì)算結(jié)果進(jìn)行精度對(duì)比.譜元法中使用二階顯式Newmark 時(shí)間步長(zhǎng)格式,時(shí)間步長(zhǎng)為0.1 ms,計(jì)算總時(shí)長(zhǎng)為0.9 s.PINN 中的波場(chǎng)時(shí)刻從第一個(gè)初始條件的時(shí)間t1=0.18 s 算起,波場(chǎng)訓(xùn)練的最終時(shí)刻為T(mén)=0.72s.采用四階勒讓德譜單元對(duì)模型在 1 00×100的網(wǎng)格上進(jìn)行空間離散.PINN 方法具備稀疏數(shù)據(jù)下的高精度學(xué)習(xí)能力,僅使用兩個(gè)時(shí)刻的50 × 50 的波場(chǎng)作為初始條件進(jìn)行訓(xùn)練,即可達(dá)到與譜元方法相關(guān)系數(shù)超過(guò)0.99 的精度.

在模型四邊各5 個(gè)單元厚度的區(qū)域邊界使用完美匹配層[37](PML)以模擬點(diǎn)源在無(wú)限介質(zhì)中的傳播.輸入波位移時(shí)程采用主頻為20 Hz 的高斯源時(shí)間函數(shù),源位置處于模型中心x=300m,z=300m.采用Sobol 序列算法在整個(gè)域中生成用于計(jì)算PDE殘差損失項(xiàng)的采樣點(diǎn)和初始條件采樣點(diǎn),空間分布如圖6 所示.采樣點(diǎn)數(shù)Np=10 000,Nic1=2500,Nic2=2500.

圖6 無(wú)限均勻介質(zhì)波動(dòng)問(wèn)題不同損失項(xiàng)時(shí)空采樣點(diǎn)分布Fig.6 Spatio-temporal sampling points distribution of different loss terms for infinite uniform medium fluctuation problem

使用5 個(gè)隱藏層,每層30 個(gè)神經(jīng)元的全連接網(wǎng)絡(luò)逼近公式(7)的解.使用學(xué)習(xí)率為6×10?3的Adam優(yōu)化器和L-BFGS 優(yōu)化器分別訓(xùn)練10000 步.優(yōu)化組合增強(qiáng)了PINN 的全局搜索和局部調(diào)優(yōu)能力,從而得到高精度的結(jié)果.模型訓(xùn)練后,可對(duì)任意分辨率下任意時(shí)空點(diǎn)的波動(dòng)方程解及解的各階偏導(dǎo)進(jìn)行預(yù)測(cè).定義某一時(shí)刻PINN 方法與譜元方法模擬波場(chǎng)之間的決定系數(shù)R2為

圖7 PINN 方法與譜元方法模擬波場(chǎng)之間的相對(duì) L2 范數(shù)誤差與決定系數(shù) R2Fig.7 Relative L2 norm error and determination coefficient R2 between the results of PINN method and spectral element method

圖8 表明: PINN 方法對(duì)波場(chǎng)的預(yù)測(cè)結(jié)果與譜元的結(jié)果基本一致,可以準(zhǔn)確捕捉到波動(dòng)在無(wú)限介質(zhì)中的傳播特性.訓(xùn)練中并未對(duì)邊界施加任何約束,內(nèi)域的波動(dòng)也可實(shí)現(xiàn)邊界透射.

圖8 無(wú)限均勻介質(zhì)波動(dòng)問(wèn)題: 譜元法與PINN 模擬的波場(chǎng)快照對(duì)比Fig.8 Wave field snapshot of spectral element method and PINN simulation

構(gòu)建不同初始條件進(jìn)行測(cè)試,發(fā)現(xiàn)經(jīng)過(guò)訓(xùn)練的神經(jīng)網(wǎng)絡(luò)具有在不同初始條件下進(jìn)行泛化的能力.僅改變上例中的源位置為x=200m,z=200m,使用譜元法生成t=0.18 s,t=0.2s 的位移場(chǎng)作為新的初始條件輸入網(wǎng)絡(luò)參與構(gòu)成損失函數(shù).在源位置為x=300m,z=300m 訓(xùn)練的神經(jīng)網(wǎng)絡(luò)可直接對(duì)新的初始條件進(jìn)行泛化,得到高精度預(yù)測(cè)結(jié)果.PINN 方法與譜元方法模擬波場(chǎng)之間的相對(duì) L2范數(shù)誤差以及決定系數(shù)R2如圖9 所示.圖10 給出了其中四個(gè)時(shí)刻模擬波場(chǎng)對(duì)比,圖中的時(shí)間均為與第一個(gè)初始條件相隔的時(shí)間.

圖9 PINN 方法使用遷移學(xué)習(xí)與譜元法模擬波場(chǎng)之間的相對(duì) L2 范數(shù)誤差與決定系數(shù) 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 方法使用遷移學(xué)習(xí)得到的波場(chǎng)快照對(duì)比Fig.10 Wave field snapshots obtained using spectral element method and PINN method using transfer learning

由圖9 可知,超出訓(xùn)練時(shí)間范圍后,相對(duì) L2范數(shù)誤差出現(xiàn)較快的上升段、決定系數(shù)R2出現(xiàn)較快的下降段.在t=1.02 s 波即將完全透出的時(shí)刻,仍與譜元方法保持0.99 以上的高相關(guān)性,由此可見(jiàn)波場(chǎng)訓(xùn)練最終時(shí)刻為0.72 s 的物理驅(qū)動(dòng)神經(jīng)網(wǎng)絡(luò)具有在新的初始條件下進(jìn)行一定時(shí)間外推的能力.模擬結(jié)果驗(yàn)證了遷移學(xué)習(xí)在近場(chǎng)波動(dòng)問(wèn)題中應(yīng)用的可行性.

3 物理驅(qū)動(dòng)深度學(xué)習(xí)波動(dòng)模擬應(yīng)用

3.1 空間不均勻介質(zhì)波動(dòng)模擬

為進(jìn)一步驗(yàn)證PINN 方法處理復(fù)雜介質(zhì)波場(chǎng)問(wèn)題的有效性,考慮空間不均勻介質(zhì)中二維波動(dòng)問(wèn)題.模型所在區(qū)域?yàn)閷?000 m、高1000 m 的矩形,在波速c=1000m/s 的基礎(chǔ)上,添加四個(gè)二維高斯混合的復(fù)雜波速分布.計(jì)算模型及介質(zhì)物理波速分布如圖11 所示,波速分布范圍為500 m/s～ 1000 m/s.

圖11 空間不均勻介質(zhì)物理波速分布Fig.11 Physical wave velocity distribution in spatially inhomogeneous media

使用譜元方法計(jì)算t=0.15 s,t=0.175s 的位移場(chǎng)作為初始條件參與訓(xùn)練,其余時(shí)刻的結(jié)果與PINN 方法結(jié)果進(jìn)行對(duì)比.每個(gè)初始條件的采樣點(diǎn)為50×50共2500 個(gè)采樣點(diǎn).譜元法中二階顯式Newmark時(shí)間步長(zhǎng)格式的時(shí)間步長(zhǎng)為0 μs,計(jì)算總時(shí)長(zhǎng)為0.7 s.PINN 中的計(jì)算時(shí)間從第一個(gè)初始條件的時(shí)間t1=0.15s 算 起,訓(xùn)練波 場(chǎng)的最 終時(shí)刻 為T(mén)=0.55s.

訓(xùn)練由兩部分構(gòu)成,首先使用學(xué)習(xí)率為6×10?4的Adam 優(yōu)化器訓(xùn)練10000 個(gè)epoch,實(shí)現(xiàn)全域收斂,再調(diào)用L-BFGS 優(yōu)化器實(shí)現(xiàn)優(yōu)化收斂,訓(xùn)練步數(shù)為40000 步.

PINN 方法與譜元方法模擬波場(chǎng)之間的相對(duì) L2范數(shù)誤差以及決定系數(shù)R2如圖12 所示.圖13 對(duì)比了不同時(shí)刻PINN 預(yù)測(cè)和譜元方法模擬的結(jié)果,圖中的時(shí)間為與第一個(gè)初始條件相隔的時(shí)間.結(jié)果表明: PINN 方法可實(shí)現(xiàn)非均勻介質(zhì)中的高精度波動(dòng)模擬.

圖12 PINN 方法與譜元方法模擬波場(chǎng)之間的相對(duì) L2 范數(shù)誤差與決定系數(shù) R2Fig.12 Relative L2 norm error and determination coefficient R2 between the results of PINN method and spectral element method

圖13 空間不均勻介質(zhì)波動(dòng)問(wèn)題: 譜元法與PINN 模擬的波場(chǎng)快照Fig.13 Spatially inhomogeneous media fluctuation problem: wave field snapshots simulated by spectral element method and PINN

3.2 起伏地形波動(dòng)模擬

物理驅(qū)動(dòng)神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)波動(dòng)模擬的一個(gè)關(guān)鍵問(wèn)題是自由邊界條件的加載方法.自由地表?xiàng)l件理論上概念明確,有限元法、譜元法可自動(dòng)滿足,地表起伏增加了求解與建模的難度[38].本文采用PINN 方法開(kāi)展相關(guān)工況驗(yàn)證.

構(gòu)建左邊界長(zhǎng)度為450 m、底邊界長(zhǎng)度為1200 m計(jì)算模型,介質(zhì)物理波速c=1000 m/s,如圖14 所示.

圖14 起伏地形內(nèi)源波動(dòng)問(wèn)題計(jì)算模型Fig.14 Computational model for endogenous fluctuation problems in undulating terrain

使用譜元方法分別計(jì)算了無(wú)限介質(zhì)中t=0.09s,t=0.1 s 的位移場(chǎng)作為初始條件、起伏地表下t=0.3s內(nèi)的位移場(chǎng)作為精度驗(yàn)證.兩種計(jì)算模型均采用四階勒讓德譜單元在 1 00×100的網(wǎng)格上進(jìn)行空間離散,其中每個(gè)初始條件使用 5 0×50共2500 個(gè)采樣點(diǎn).吸收邊界仍使用5 個(gè)單元厚度的完美匹配層(PML).輸入波位移時(shí)程采用主頻為20 Hz 的高斯源時(shí)間函數(shù),源位置為x=250 m,z=750m.采用Sobol 序列算法在整個(gè)域中生成用于計(jì)算PDE 殘差損失項(xiàng)和自由邊界條件損失項(xiàng)的采樣點(diǎn),采樣點(diǎn)數(shù)分別為Np=10 000,Nbc=2500.采樣點(diǎn)的空間分布如圖15 所示.

圖15 起伏地形內(nèi)源波動(dòng)問(wèn)題不同損失項(xiàng)時(shí)空采樣點(diǎn)分布Fig.15 Distribution of spatial and temporal sampling points for different loss terms for endogenous fluctuation problems

使用5 個(gè)隱藏層,每層30 個(gè)神經(jīng)元的全連接網(wǎng)絡(luò).調(diào)用學(xué)習(xí)率為6×10?3的Adam 優(yōu)化器訓(xùn)練10000個(gè)epoch,再使用L-BFGS 優(yōu)化器訓(xùn)練40000 步.PINN 方法與譜元方法模擬波場(chǎng)之間的相對(duì) L2范數(shù)誤差以及決定系數(shù)R2如圖16 所示.圖17 展示了三個(gè)時(shí)刻的PINN 預(yù)測(cè)和譜元方法模擬結(jié)果對(duì)比,圖中的時(shí)間為與第一個(gè)初始條件相隔的時(shí)間.

圖16 PINN 方法與譜元方法模擬波場(chǎng)之間的相對(duì) L2 范數(shù)誤差與決定系數(shù) R2Fig.16 Relative L2 norm error and determination coefficient R2 between the results of PINN method and spectral element method

圖17 起伏地表內(nèi)源波動(dòng)問(wèn)題: 譜元法與PINN 模擬的波場(chǎng)快照Fig.17 Internal fluctuation problem of undulating surface: wave field snapshots of spectral element method and PINN simulation

4 結(jié)論與展望

本文結(jié)合波動(dòng)數(shù)值模擬基本理論與物理驅(qū)動(dòng)深度學(xué)習(xí),建立了包括一維無(wú)源、二維波動(dòng)模擬方法.通過(guò)與理論解、譜元計(jì)算結(jié)果等對(duì)比,驗(yàn)證了方法應(yīng)用于均質(zhì)場(chǎng)地、空間不均勻及復(fù)雜地形場(chǎng)地的有效性.得到的主要結(jié)論與展望如下.

(1) 物理驅(qū)動(dòng)深度學(xué)習(xí)PINN 方法可實(shí)現(xiàn)波動(dòng)模擬,并可結(jié)合稀疏初始波場(chǎng)數(shù)據(jù),實(shí)現(xiàn)高精度泛化.方法具備無(wú)網(wǎng)格、精細(xì)化,波場(chǎng)透射等優(yōu)勢(shì).物理驅(qū)動(dòng)神經(jīng)網(wǎng)絡(luò)僅需存儲(chǔ)網(wǎng)絡(luò)參數(shù),對(duì)硬件存儲(chǔ)要求較低.

(2) 針對(duì)典型工況,訓(xùn)練形成的神經(jīng)網(wǎng)絡(luò)具有在不同初始條件的泛化能力,結(jié)合遷移學(xué)習(xí),可顯著提高網(wǎng)絡(luò)的訓(xùn)練效率.

(3) 通過(guò)“軟約束”可在PINN 方法中嵌入應(yīng)力型及位移型等不同邊界條件.由于物理驅(qū)動(dòng)深度學(xué)習(xí)計(jì)算損失函數(shù)的靈活性,可展望PINN 方法與有限差分、偽譜法等多類型數(shù)值算法開(kāi)展耦合,以實(shí)現(xiàn)對(duì)空間大尺度、高精度復(fù)雜波動(dòng)問(wèn)題的模擬.