基于物理信息驅動神經網絡的三維初至波旅行時計算方法

都國寧，譚軍*,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深度神經網絡在三維速度模型中計算初至波旅行時的方法。經不同速度模型實驗，結果表明本文所提方法求解的旅行時比目前常用的有限差分法計算結果更準確或結果相似，并且適用于復雜的三維速度模型。相對于傳統深度學習算法，本文構建的神經網絡是基于物理模型驅動的，即在損失函數中加入了三維程函方程，使計算結果更符合物理規律。此外，在給定三維速度模型中選擇一部分網格點作為震源位置進行訓練學習后，即可輸入任意震源位置迅速求解地震波旅行時，在多震源應用中表現出有限差分法不可比擬的效率優勢。