都國寧，譚軍*,2,3，宋鵬,2,3，解闖，王紹文

(1.中國海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島 266100； 2.青島國家海洋科學(xué)與技術(shù)實驗室，山東青島 266100； 3.中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點實驗室，山東青島 266100)

0 引言

精確求取初至波旅行時對于實現(xiàn)靜校正[1-2]、基爾霍夫積分偏移[3]、旅行時反演[4-9]、地震定位[10]、偏移成像[11]等具有重要意義。傳統(tǒng)求取初至波旅行時方法采用射線追蹤類算法，包括試射法[12-14]和彎曲法[15-16]，但計算效率較低。由于復(fù)雜速度模型中兩點之間的射線路徑存在多種可能，計算極易陷入局部收斂，為了克服射線追蹤類算法的缺陷，程函方程類算法得以發(fā)展。Vidale[17]基于盒式擴張的思想提出采用有限差分法求解二維程函方程，但這種算法并不完全符合波前傳播規(guī)律，在因果性方面存在缺陷，由于忽略了初至能量可能迂回傳播的情況，致使該算法存在不穩(wěn)定性，不能取得真正的全局最小旅行時，當(dāng)網(wǎng)格間距較小時還會帶來巨大的計算量； Qin等[18]使用波前擴張理論改進了有限差分算法的因果性；VanTrier等[19]使用迎風(fēng)有限差分算子提高了算法的穩(wěn)定性；Sethian[20]提出快速推進法(Fast Marching Method,FMM)，利用逆風(fēng)差分格式求解局部程函方程,采用窄帶延拓重建旅行時波前,利用堆選排技術(shù)保存旅行時,將最小旅行時放在堆的頂部，顯著縮短了尋找極小值的時間; Qian等[21]將快速清理法(Fast Sweeping Method,FSM)用于二維程函方程求解，進一步提高了有限差分算法的計算效率，其主要思想是基于因果關(guān)系將旅行時波場傳播的方向分成有限個組,對每一組分別利用Gauss-Seidel迭代方法求解非線性逆風(fēng)差分格式離散化后的方程組。一些學(xué)者提出因式分解形式的程函方程[22-23]，解決了震源奇異性問題，進一步提高了有限差分算法的準確度。Vidale[24]還將有限差分算法擴展到三維，此后便成為求取三維旅行時最常用的方法; Soukina等[25]將三維有限差分算法應(yīng)用于各向異性介質(zhì)旅行時的求取；Hole等[26]使用三維有限差分算法求解反射波旅行時。相對于二維模型，求取三維地震旅行時對計算精度和效率要求更高，而高精度的高階有限差分會導(dǎo)致巨大的計算成本[27-28]。如何在求解三維程函方程時平衡精確度與計算成本是求解旅行時亟待解決的問題。

深度學(xué)習(xí)技術(shù)的蓬勃發(fā)展為解決以上問題拓展了新途徑。20世紀90年代，諸多學(xué)者提出了用神經(jīng)網(wǎng)絡(luò)求解偏微分方程的想法[29-30]。Sirignano等[31]提出使用全連接神經(jīng)網(wǎng)絡(luò)(Fully Connected Neural Network,FCNN)求解偏微分方程的無網(wǎng)格化方法；Tompson等[32]使用卷積神經(jīng)網(wǎng)絡(luò)提高了稀疏線性偏微分方程的求解效率，得到Navier-Stokes方程的數(shù)值解。針對傳統(tǒng)神經(jīng)網(wǎng)絡(luò)沒有考慮偏微分方程本身攜帶的物理信息并且缺少物理意義可解釋性的缺陷，Raissi等[33]提出物理信息驅(qū)動的神經(jīng)網(wǎng)絡(luò)(Physics-informed Neural Network，PINN)。與傳統(tǒng)深度學(xué)習(xí)網(wǎng)絡(luò)不同，PINN在利用FCNN的函數(shù)擬合功能實現(xiàn)偏微分方程近似的基礎(chǔ)上，在神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程中加入偏微分方程和實際物理條件的約束，使求解偏微分方程的結(jié)果更具實際物理意義。近年來，PINN在地球物理領(lǐng)域得到了廣泛的應(yīng)用。Xu等[34]將PINN應(yīng)用于速度模型的反演；Karimpouli等[35]使用PINN求解波動方程進行正演模擬；Song等[36]將PINN用于VTI介質(zhì)頻率域聲波方程的求解；Waheed等[37]使用PINN求解二維程函方程并將其應(yīng)用于層析成像，結(jié)果顯示PINN在二維初至波旅行時的求取中表現(xiàn)出高于傳統(tǒng)方法的效率和準確度。

本文基于PINN實現(xiàn)了三維程函方程的高效、高精度求解。模型實驗結(jié)果顯示，基于PINN的三維初至波旅行時計算方法相對于傳統(tǒng)的有限差分法有更高的精度和更高的計算效率。

1 程函方程的有限差分數(shù)值解

三維形式的程函方程為

(1)

式中：T(x,y,z)為(x,y,z)處的地震波旅行時；v(x,y,z)為(x,y,z)處的速度。

為了便于表示，式(1)可寫為

(2)

(3)

為求解三維程函方程，Vidale[24]提出了有限差分算法，將三維速度模型剖分為如圖1所示的若干正方體網(wǎng)格。假設(shè)已知A點旅行時和各網(wǎng)格點的速度，且網(wǎng)格間距為h，則有

圖1 三維差分網(wǎng)格示意圖

TB=TA+hsB

(4)

(5)

(6)

式(4)～式(6)分別為一維、二維、三維有限差分算子，使用以上算子對所有網(wǎng)格點進行差分計算，即可求得整個計算區(qū)域的旅行時。

有限差分算法憑借較高的精度和計算效率得到了廣泛的應(yīng)用，但該方法存在不可避免的震源奇異性問題[28]。為此，需要將程函方程轉(zhuǎn)化為因式分解形式[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]。因此，本文嘗試將深度學(xué)習(xí)與旅行時計算相結(jié)合，提出兼具精度與計算效率的基于PINN的地震旅行時求取方法。

2 基于PINN求取地震旅行時

2.1 PINN簡介

由通用近似算法[38]可知，一個輸入層有n個神經(jīng)元、輸出層有m個神經(jīng)元的神經(jīng)網(wǎng)絡(luò)可以用來表示任意一個多維非線性函數(shù)u：Rn→Rm，而偏微分方程的建模過程也是尋找滿足約束條件的非線性函數(shù)，兩者具有相通之處。得益于在深度神經(jīng)網(wǎng)絡(luò)(DNN)中廣泛使用的自動微分技術(shù)，在設(shè)計神經(jīng)網(wǎng)絡(luò)的損失函數(shù)時，可以融入偏微分方程中微分形式的約束條件，以獲得包含物理模型約束的神經(jīng)網(wǎng)絡(luò)，即PINN。PINN的基本結(jié)構(gòu)是通過FCNN近似表示一個函數(shù)，再利用自動微分技術(shù)求出偏微分方程殘差和初邊值殘差約束，并作為正則項添加至損失函數(shù)中，最后利用梯度下降法等優(yōu)化算法獲得神經(jīng)網(wǎng)絡(luò)權(quán)重參數(shù)和偏微分方程物理參數(shù)。

(10)

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

自動微分的反向模式就是反向傳播算法的一般化，其思路是根據(jù)計算路徑從后向前計算，依次得到對每個中間變量節(jié)點的偏導(dǎo)數(shù)，直至自變量節(jié)點處。在每個節(jié)點處根據(jù)其后續(xù)節(jié)點計算導(dǎo)數(shù)值，整個過程對應(yīng)于多元復(fù)合函數(shù)求導(dǎo)時從最外層逐步向內(nèi)側(cè)求導(dǎo)。自動微分不涉及差分近似誤差，因此能夠較準確地計算導(dǎo)數(shù)，從而求出程函方程殘差和初邊值殘差約束作為損失函數(shù)，再通過合適的優(yōu)化算法進行神經(jīng)網(wǎng)絡(luò)參數(shù)的迭代優(yōu)化，當(dāng)損失函數(shù)達到最小時即實現(xiàn)了程函方程的精確求解。

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

利用神經(jīng)網(wǎng)絡(luò)近似表示函數(shù)映射關(guān)系求解三維程函方程，定義一個損失函數(shù)，使訓(xùn)練集中的因式分解形式程函方程的殘差最小。求解過程主要包括：

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

(2)定義包含了三維程函方程的損失函數(shù)，并在隨機分配的網(wǎng)格上進行采樣，構(gòu)建用于神經(jīng)網(wǎng)絡(luò)訓(xùn)練的數(shù)據(jù)集；

(3) 通過DNN的自動微分算法計算τ(x,y,z)相對于空間坐標的偏導(dǎo)數(shù)；

(4) 選擇合適的優(yōu)化器，通過更新網(wǎng)絡(luò)參數(shù)來最小化損失函數(shù)。

為了便于描述，本文假設(shè)震源(xs,ys,zs)處有τ(xs,ys,zs)=1。用多層深度神經(jīng)網(wǎng)絡(luò)Rτ近似(x,y,z)位置的旅行時未知因式τ(x,y,z)，即

(11)

激活函數(shù)在神經(jīng)網(wǎng)絡(luò)參數(shù)的優(yōu)化中也起重要作用。許多研究證明，局部自適應(yīng)激活函數(shù)比傳統(tǒng)的固定激活函數(shù)具有更好的學(xué)習(xí)能力[40]，因此除了在最后一層使用ReLu線性激活函數(shù)外，所有的隱藏層都使用了局部自適應(yīng)反正切(locally adaptive arctangent)函數(shù)，如此可在每個神經(jīng)元的激活函數(shù)中使用可擴展參數(shù)，改變激活函數(shù)的斜率以提高網(wǎng)絡(luò)的學(xué)習(xí)效率。

使用均方差(MSE)構(gòu)建損失函數(shù)

(12)

式中：I表示所有采樣點的集合；NI為采樣點的數(shù)量；X*表示采樣點位置坐標；L為因式分解形式程函方程的殘差，其表達式為

s2(x,y,z)

(13)

(14)

在PINN的優(yōu)化過程中，損失函數(shù)中不同組成部分的收斂效率存在差異[41]。針對損失函數(shù)各部分的權(quán)值分配方法有兩種，即自適應(yīng)權(quán)重和固定權(quán)重。Waheed等[37]證明了在PINN的學(xué)習(xí)過程中，自適應(yīng)權(quán)重分配方法能夠有更佳的收斂效率。因此，本文采取自適應(yīng)權(quán)重分配策略，即通過反向傳播的梯度值調(diào)整分配給損失函數(shù)中的不同項的權(quán)重，從而提高收斂效率。

基于PINN求解三維程函方程的流程如圖2所示，(x*,y*,z*)表示訓(xùn)練集中隨機選擇的各點坐標。將這些點的坐標輸入到隨機初始化的神經(jīng)網(wǎng)絡(luò)中進行訓(xùn)練，并將各點對應(yīng)的速度v(x*,y*,z*)、已知旅行時參數(shù)T0(x*,y*,z*)及空間導(dǎo)數(shù)T0(x*,y*,z*)作為網(wǎng)絡(luò)的已知信息，通過最小化損失函數(shù)的過程實現(xiàn)對神經(jīng)網(wǎng)絡(luò)參數(shù)的迭代優(yōu)化。神經(jīng)網(wǎng)絡(luò)訓(xùn)練完成后，輸入待求位置的坐標(x,y,z)，使用神經(jīng)網(wǎng)絡(luò)預(yù)測并輸出未知旅行時因式再與T0(x,y,z)相乘，即可求得最終的旅行時

圖2 基于PINN求解三維程函方程的流程圖

(15)

3 數(shù)值模擬實驗與討論

為了研究本文所提基于PINN的三維旅行時計算方法的適用性，設(shè)計均勻速度模型、水平層狀等速度梯度模型、傾斜層狀等速度梯度模型及局部三維Marmousi模型,開展模擬實驗，并將本文方法與三維有限差分法計算結(jié)果進行了比較。所有模型實驗均在Intel(R) Core(TM) i7-4720HQ的CPU上進行，模型參數(shù)見表1。

表1 模型參數(shù)統(tǒng)計表

3.1 均勻模型

均勻模型中的初至波解析旅行時可以通過距離與速度之比直接求出。截取x=3 m、y=14 m、z=36 m三個剖面分別展示PINN計算的旅行時結(jié)果與有限差分法計算結(jié)果及解析值的對比。如圖3所示，在三個剖面中PINN計算結(jié)果都與理論旅行時相差無幾，而有限差分求解結(jié)果則在接近模型邊界的個別位置相對于解析值出現(xiàn)了偏差。

圖3 均勻模型3種方法旅行時計算結(jié)果對比(★為震源位置，下同)

圖4展示了上述三個剖面應(yīng)用有限差分法(圖4上)與本文方法(圖4下)的計算結(jié)果相對于解析旅行時的絕對誤差?？梢钥吹?，基于PINN求解的旅行時計算結(jié)果絕對誤差較小。

圖4 均勻模型有限差分法(上)、本文方法(下)計算結(jié)果相對于解析旅行時的絕對誤差

本次模型實驗中，采用PINN方法訓(xùn)練過程共耗費約19 min，訓(xùn)練完成后對待求旅行時進行預(yù)測僅需3 s，而使用有限差分法進行相同操作則需要約5 min。雖然在神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過程耗時較多，但三維旅行時的應(yīng)用通常需要計算多個震源的旅行時，隨著震源數(shù)量的增加，本文提出的方法會展現(xiàn)更明顯的效率優(yōu)勢。

3.2 水平層狀等速度梯度模型

對于具有恒定速度梯度的模型，解析旅行時的方程為[42]

(16)

式中：T′(x,y,z)是從震源(xs,ys,zs)到某個網(wǎng)格點(x,y,z)的解析旅行時；gx、gy、gz分別表示沿x、y、z三個方向的速度梯度分量。

創(chuàng)建的水平層狀等速度梯度模型在x方向的剖面如圖5所示，設(shè)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求取旅行時的結(jié)果與有限差分法計算結(jié)果及解析值的對比(圖6)?？梢钥吹剑谌齻€剖面中PINN計算結(jié)果都與理論旅行時幾乎相等，而有限差分求解結(jié)果則在震源附近的個別位置相對于解析值出現(xiàn)了細微偏差，并且該偏差沿同一方向影響了整個速度模型的旅行時計算結(jié)果。

圖6 水平層狀等速度梯度模型3種方法旅行時計算結(jié)果對比

圖7為上述三個剖面的有限差分算法計算結(jié)果(圖7上)和本文所用方法計算結(jié)果(圖7下)相對于解析旅行時的絕對誤差。由圖可見，基于PINN求解的旅行時絕對誤差仍然小于有限差分法計算結(jié)果。

圖7 水平層狀模型有限差分法(上)、本文方法(下)計算結(jié)果相對于解析旅行時的絕對誤差

本次模型實驗中，PINN方法訓(xùn)練過程共耗費約25 min，訓(xùn)練完成后對整個速度模型的旅行時預(yù)測耗時4 s，而有限差分法計算旅行時約耗時5 mim。本文方法在多炮計算效率方面展現(xiàn)出巨大潛力。

3.3 傾斜層狀等速度梯度模型

為了證明本文方法的泛化性以及在多震源計算中的效率優(yōu)勢，建立一個傾斜層狀等速度梯度模型，x、y、z三個方向的速度變化梯度分別為gx=1 s-1、gy=8 s-1、gz=5 s-1。選擇如圖8三個剖面中的三角形所示的48個位置作為訓(xùn)練集震源，生成旅行時訓(xùn)練集數(shù)據(jù)；以不屬于訓(xùn)練集的位置(25 m，25 m，25 m)作為測試集的震源。

圖8 傾斜層狀等速度梯度模型

選取x=15 m、y=30 m、z=20 m三條剖面展示PINN法求取旅行時的結(jié)果與有限差分法計算結(jié)果以及解析值的對比。如圖9所示，在三個剖面中PINN法計算結(jié)果都與理論旅行時幾乎完全擬合，而有限差分求解結(jié)果在對角線方向相對于解析值出現(xiàn)了較明顯的偏差。

圖9 傾斜層狀等速度梯度模型3種方法旅行時計算結(jié)果對比

圖10為上述三個剖面中有限差分法計算結(jié)果(圖10上)與本文方法計算結(jié)果(圖10下)相對于解析旅行時的絕對誤差?？梢娫谳^復(fù)雜的速度模型中，基于PINN求解的旅行時的精度仍然優(yōu)于有限差分法計算結(jié)果。

圖10 傾斜層狀模型有限差分法(上)與本文方法(下)計算結(jié)果相對于解析旅行時的絕對誤差

為了驗證多震源計算的效率優(yōu)勢，在y=11 m、y=12 m、y=13 m三個剖面上選取123個網(wǎng)格點作為震源，且這三個剖面上的點均與訓(xùn)練集中的震源點無重合。使用本文方法和有限差分法分別計算所有震源的旅行時，圖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剖面上的旅行時計算結(jié)果。在該模型試驗123炮的計算中，PINN方法訓(xùn)練、預(yù)測過程共耗時約68 min，而有限差分法計算旅行時共耗時約458 min，PINN方法體現(xiàn)出了相當(dāng)明顯的效率優(yōu)勢。

圖11 多震源位置3種方法旅行時計算部分計算結(jié)果對比(x=20 m剖面)

3.4 三維Marmousi模型

選取x=40 m、y=60 m、z=50 m三條剖面，將二維Marmousi模型中速度變化十分復(fù)雜的部分拼接為如圖12所示的局部三維Marmousi模型。

圖12 傾斜層狀等速度梯度模型

圖13展示了訓(xùn)練完成的神經(jīng)網(wǎng)絡(luò)的旅行時預(yù)測結(jié)果與有限差分法計算結(jié)果的對比，可以看到，在三個剖面中PINN計算結(jié)果都與有限差分法計算結(jié)果基本吻合。

圖13 三維Marmousi模型2種方法旅行時計算結(jié)果對比

圖14為上述三個剖面中有限差分法計算結(jié)果與本文所用方法計算結(jié)果之間的絕對誤差。由于無法計算Marmousi模型的解析旅行時，故而無法進行精度的對比，但是通過本文方法與有限差分法所求結(jié)果的絕對誤差，可以驗證PINN方法即使對于復(fù)雜構(gòu)造模型也具有較高的穩(wěn)定性。另外，該模型試驗中有限差分法計算旅行時用時約536 min，而PINN方法訓(xùn)練、預(yù)測過程共耗時約145 min，在復(fù)雜模型中同樣表現(xiàn)出明顯的效率優(yōu)勢。

圖14 三維Marmousi模型中有限差分法與本文方法旅行時求解結(jié)果的絕對誤差

通過數(shù)值模擬實驗，證明了PINN算法具有效率上的絕對優(yōu)勢，且在簡單速度模型中應(yīng)用精度也有所提高。但對于復(fù)雜模型則無法驗證所求解旅行時的準確度，只能說明PINN算法能夠得到與有限差分算法相似的結(jié)果。

4 結(jié)束語

本文提出一種基于PINN深度神經(jīng)網(wǎng)絡(luò)在三維速度模型中計算初至波旅行時的方法。經(jīng)不同速度模型實驗，結(jié)果表明本文所提方法求解的旅行時比目前常用的有限差分法計算結(jié)果更準確或結(jié)果相似，并且適用于復(fù)雜的三維速度模型。相對于傳統(tǒng)深度學(xué)習(xí)算法，本文構(gòu)建的神經(jīng)網(wǎng)絡(luò)是基于物理模型驅(qū)動的，即在損失函數(shù)中加入了三維程函方程，使計算結(jié)果更符合物理規(guī)律。此外，在給定三維速度模型中選擇一部分網(wǎng)格點作為震源位置進行訓(xùn)練學(xué)習(xí)后，即可輸入任意震源位置迅速求解地震波旅行時，在多震源應(yīng)用中表現(xiàn)出有限差分法不可比擬的效率優(yōu)勢。