姚 時，侯 爵3,，黃躍鵬，徐 濤，白志明，高正輝

1 中國科學(xué)院地質(zhì)與地球物理研究所 巖石圈演化國家重點實驗室，北京 100029

2 中國科學(xué)院大學(xué)，北京 100049

3 中國地震局地球物理研究所，北京 100081

4 中國科學(xué)院地球科學(xué)研究院，北京100029

5 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引 言

地震波走時計算在層析成像（徐濤等, 2014,2015; Liang et al., 2016; 張明輝等, 2019; 林吉焱等,2020）、偏移成像（Gray and May, 1994; Sun, 1998;劉國峰等, 2009）和微震定位（Ben-Zion et al.,2015; Inbal et al., 2015）中都發(fā)揮著重要的作用. 目前走時計算的方法可分為射線追蹤類（Julian and Gubbins, 1977; 徐濤等, 2004; Xu et al., 2006; Rawlinson et al., 2008; 李飛等, 2013）和程函方程數(shù)值求解類（Vidale, 1988, 1990; Sethian, 1996, 1999; Rawlinson and Sambridge, 2004; Zhao, 2005; Qian, 2007a,2007b; Treister and Haber, 2016）. 射線追蹤類方法需要先計算射線路徑，然后沿著射線路徑計算走時，在檢波點較多的情況下計算效率較低，且在模型復(fù)雜時會出現(xiàn)陰影區(qū)等問題（Xu et al., 2010, 2014）.程函方程數(shù)值求解類方法有效避免了上述問題，此類方法中最常用的是快速推進法（fast marching method, FMM）和快速掃描法（fast sweeping method, FSM）. FMM （Sethian, 1996, 1999; Rawlinson and Sambridge, 2004; Treister and Haber, 2016）是一種基于網(wǎng)格的差分數(shù)值計算格式，主要包括利用迎風(fēng)差分格式求解局部程函方程獲取走時和利用窄帶技術(shù)模擬波前的傳播過程兩個步驟. FSM（Zhao, 2005; Qian et al., 2007a, 2007b）是一種利用Gauss-Seidel 或者Gauss-Jacobi 迭代格式進行迎風(fēng)差分求解程函方程的迭代方法，首先基于因果關(guān)系將走時場傳播的方向分成有限數(shù)量的組，對于每一組分別利用選定的迭代方法求解非線性逆風(fēng)差分格式離散化后的方程組. FMM 和FSM 對速度變化很大的非均勻介質(zhì)依然有很好的穩(wěn)定性和適用性，均可以準確地計算地震波走時（蘭海強等，2012）. FMM 和FSM 都屬于有限差分類方法，在計算走時的過程中需要求解由每一個震源激發(fā)的走時場，而走時的計算是一個消耗時間和存儲空間的過程，因此，計算走時的時間和內(nèi)存消耗都會隨著震源和網(wǎng)格數(shù)量的增加而增加.

深度學(xué)習(xí)是一種以神經(jīng)網(wǎng)絡(luò)為架構(gòu)、對數(shù)據(jù)進行表征學(xué)習(xí)的算法（Hinton et al., 2006; Goodfellow et al., 2016），目前，深度學(xué)習(xí)方法在許多地震學(xué)問題中已經(jīng)開始發(fā)揮較大的作用（楊旭等，2021），包括以速度模型為輸入預(yù)測介質(zhì)的聲波響應(yīng)（Moseley et al., 2018）、提高黏彈性模擬的速度（DeVries et al., 2017）、使用無監(jiān)督神經(jīng)網(wǎng)絡(luò)解決偏微分方程的正反演問題（Bar and Sochen, 2019）、利用卷積神經(jīng)網(wǎng)絡(luò)進行地震波走時層析成像（Fan and Ying, 2019）、地震波形反演（Sun et al.,2020）、利用卷積神經(jīng)網(wǎng)絡(luò)進行地震波形自動分類與識別（趙明等，2019）和地震數(shù)據(jù)隨機噪聲去除（韓衛(wèi)雪等，2018）等. 此外，還可以向網(wǎng)絡(luò)中加入數(shù)學(xué)物理定律等先驗知識構(gòu)建基于物理信息高效的神經(jīng)網(wǎng)絡(luò)，以解決更復(fù)雜的物理問題（Raissi et al., 2019）. Smith 等（2020）提出了一種名為EikoNet的使用深度神經(jīng)網(wǎng)絡(luò)求解程函方程的方法. EikoNet方法在給定速度模型和其中任意兩點的情況下可以無網(wǎng)格地快速確定兩點之間的走時. 該方法通過在三維空間中采樣生成訓(xùn)練樣本，以給定的速度模型為標簽實現(xiàn)訓(xùn)練過程中對網(wǎng)絡(luò)的優(yōu)化. 并且，在訓(xùn)練過程中，空間梯度都是利用神經(jīng)網(wǎng)絡(luò)的可微性解析計算的. 此外，EikoNet 可以大規(guī)模并行、訓(xùn)練和預(yù)測都可以在GPU 上進行，這樣就使得其計算效率顯著提高，并同時大大降低其在內(nèi)存上的消耗.

本文使用EikoNet 方法和FMM 在均勻速度模型、塊狀速度模型、層狀速度模型和棋盤格速度模型上進行了實驗，對比了二者的效率和精度.

1 方法原理

1.1 EikoNet

圖 1 處理工作流概述.（a）由全連接層和殘差塊組成的神經(jīng)網(wǎng)絡(luò)體系結(jié)構(gòu)，每個殘差塊由3 個全連接層組成，有512 個神經(jīng)元，ELU 激活應(yīng)用于所有隱藏層.（b）T s→r 和 Vr的程函方程總結(jié).（c）在整個三維空間采樣源-接收對構(gòu)建訓(xùn)練數(shù)據(jù)集.（d）通過最小化與預(yù)測和已知速度值相關(guān)的損失函數(shù)進行網(wǎng)絡(luò)訓(xùn)練.（e）通過傳遞用戶定義的源-接收對檢查神經(jīng)網(wǎng)絡(luò)輸出（修改自Smith et al., 2020）Fig. 1 Overview of the processing workflow. (a) Neural network architecture composed of fully connected layers and residual blocks. Each residual block is composed of three fully connected layers with 512 neurons. ELU activations are applied on all hidden layers. (b) Summary of Eikonal equation for Ts→r and Vr. (c) Sampling of source-receiver pairs across the 3-D volume to build the training data set. (d) Network training through the minimization of loss function relating predicted and observed velocity values. (e) Inspection of neural network outputs by passing user-defined source-receiver pairs (modified from Smith et al., 2020)

如圖1 所示，EikoNet 是一個多層的前饋神經(jīng)網(wǎng)絡(luò)，由一系列全連接層組成的殘差塊和其后的非線性單元構(gòu)成，使用pytorch 深度學(xué)習(xí)框架搭建，不同的殘差塊可以使用跳層連接的方式實現(xiàn)，緩解了在深度神經(jīng)網(wǎng)絡(luò)中增加深度帶來的梯度消失問題，并解決了深度神經(jīng)網(wǎng)絡(luò)的退化問題，并且可以通過增加一定的深度來提高神經(jīng)網(wǎng)絡(luò)計算的準確率.

EikoNet 的輸入是6 維的數(shù)據(jù)，輸出是一維的數(shù)據(jù)，其輸入和輸出結(jié)構(gòu)各由兩個線性層組成，維數(shù)變化分別為 6 →32 →512和 512 →32 →1. 其實現(xiàn)大致可分為以下四個步驟：

（1）模型采樣

使用加權(quán)隨機采樣的方法從連續(xù)的三維介質(zhì)中隨機采樣源點和接收點，權(quán)值通過網(wǎng)絡(luò)預(yù)測的速度和真實速度的相對誤差來定義，并將每個點的速度標注在接收點的位置來建立訓(xùn)練數(shù)據(jù)集.

采取的樣本形式為 (

x

,

y

)

，

x

是源點和接收點坐標的集合，

y

為接收點處的速度.

式中，

x

,

y

,

z

表示源點和接收點的坐標，下標s 和r分別代表源點和接收點.

（2）正向傳播

源點和接收點之間經(jīng)分解的走時可以表示為:

式中，

T

=‖→

x

-→

x

‖，表示接收點到源點之間的距離函數(shù)； τ是實際走時和

v

=1的均勻速度模型走時的偏差.

程函方程是一個非線性的一階偏微分方程，其一般形式為：

式中，

s

和

v

分別代表慢度和速度.

將式（3）代入（4）并展開可以得到因式程函方程：

訓(xùn)練一個深度神經(jīng)網(wǎng)絡(luò)

f

來預(yù)測輸入源點和接收點之間的走時偏差. 走時偏差可以表示為：

開始訓(xùn)練之前 τ是未知的，我們可以賦予 τ一個初始值，利用式（5）可以計算得到速度，將此速度定義為預(yù)測速度

v

?，此過程中涉及到的空間梯度都可以利用神經(jīng)網(wǎng)絡(luò)的可微性進行計算.

（3）反向傳播

利用

v

? 和 已知的速度

v

可以定義一個均方誤差損失函數(shù)：

將計算得到的失配函數(shù)反向傳播梯度到神經(jīng)網(wǎng)絡(luò)的參數(shù)，通過迭代更新每個神經(jīng)元的參數(shù)來最小化目標函數(shù)，在此過程中就達到了優(yōu)化網(wǎng)絡(luò)的作用.

（4）結(jié)果輸出

針對每一個速度模型都需要進行一次神經(jīng)網(wǎng)絡(luò)的訓(xùn)練，在本文中設(shè)置網(wǎng)絡(luò)訓(xùn)練的迭代次數(shù)為200，每次訓(xùn)練大約需要2～3 個小時. 神經(jīng)網(wǎng)絡(luò)訓(xùn)練完成后，只需要保存90 M 左右的神經(jīng)網(wǎng)絡(luò)模型，在三維速度模型中給定一定數(shù)量的源點和接收點就可以快速輸出預(yù)測的走時場.

1.2 快速推進法（FMM）

本文使用快速推進法（FMM）計算走時，主要用于與EikoNet 方法的計算結(jié)果進行對比，F(xiàn)MM 方法求解的是一般形式的程函方程式（4）.首先，我們可以使用迎風(fēng)差分格式近似程函方程的時間梯度項，基于網(wǎng)格的數(shù)值格式求解局部程函方程獲取走時.

隨后，使用窄帶技術(shù)模擬波前的傳播過程，如圖2a 所示，窄帶技術(shù)將計算區(qū)域內(nèi)所有網(wǎng)格節(jié)點分為走時計算已完成的完成點（黑色）、已計算走時但還未確定最小值的窄帶點（灰色）和未進行走時計算的遠離點（白色）. 圖2b 展示了窄帶技術(shù)從源點向外拓展的幾個步驟，重復(fù)此過程可以計算得到所有網(wǎng)格點的走時.

1.3 快速掃描法（FSM）

本文使用精細網(wǎng)格因式分解的FSM（周小樂等，2020）計算走時，主要用作參考解對比EikoNet方法和FMM 方法的精度.

圖 2 （a）窄帶技術(shù)原理示意圖. （b）從源點開始的窄帶拓展示例（修改自Rawlinson and Sambridge, 2004）Fig. 2 (a) The principle of the narrow-band method. (b) Example of how the narrow band evolves from a source point(modified from Rawlinson and Sambridge, 2004)

快速掃描法主要有兩個步驟，首先，給所有的網(wǎng)格點賦初始走時值，后續(xù)計算時會被更新；在交替的4個方向上，使用Gauss-Seidel 型迭代法求解離散的非線性方程組，分別按4 個方向掃描.

2 數(shù)值模型實驗

2.1 精度效率對比方法

2.1.1 精度對比方法

本文在不同模型下分別設(shè)置了走時解析解和精細網(wǎng)格下計算的走時數(shù)值解作為走時參考解，通過計算EikoNet 方法和FMM 與走時參考解的相對誤差來衡量精度.

（1）均勻速度模型的走時參考解

在均勻速度模型中，我們使用程函方程的解析解作為走時參考解.

（2）其他速度模型的走時參考解

對于塊狀速度模型、層狀速度模型和棋盤格速度模型，則使用1.3 節(jié)中提到的因式分解的FSM來計算程函方程走時作為走時參考解.

我們采用以下步驟確定其走時參考解：

①給定一個初始的網(wǎng)格間距，使用上述方法計算走時場并保存走時計算結(jié)果；

②不斷縮小網(wǎng)格間距，在每一個網(wǎng)格間距下分別計算一次走時場，將每一次的走時計算結(jié)果與上一次的走時計算結(jié)果進行對比計算誤差；

③如果誤差不再發(fā)生較大變化，我們就將此網(wǎng)格間距下的走時計算結(jié)果當(dāng)作走時參考解，否則，重復(fù)步驟②.

最終，我們選擇了0.05 km 的網(wǎng)格間距下得到的走時計算結(jié)果作為走時參考解.

（3）走時場相對誤差

設(shè)走時參考解為

t

，使用EikoNet 計算得到的走時為

t

，使用FMM 計算得到的走時為

t

，EikoNet的計算誤差為

P

，F(xiàn)MM 的計算誤差為

P

，

P

和

P

的表達式分別定義為：

2.1.2 效率對比方法

在進行走時場計算的過程中，分別記錄EikoNet方法和FMM 的計算耗時來衡量兩種方法的效率.EikoNet 方法的走時計算是在Nvidia GeForce RTX 2080 Ti 上進行的，F(xiàn)MM 的走時計算是在Intel Core i9-9900K 上進行的.

2.2 計算實例

本文設(shè)置了四種速度模型，均為

X

×

Y

×

Z

=20 km×20 km×20 km 的三維模型，設(shè)置了兩種震源位置以驗證EikoNet 方法在不同的震源位置均具有較高計算精度，均勻速度模型、塊狀速度模型和層狀速度模型分別對應(yīng)于均勻介質(zhì)、速度異常體和層狀介質(zhì)三種實際地質(zhì)情況，棋盤格速度模型則是用來檢驗EikoNet 方法在計算過程中適應(yīng)正負速度異常的能力，此外，使用FMM 進行計算時將網(wǎng)格間距設(shè)置為0.1 km.

圖 4 均勻速度模型走時計算相對誤差. （a）、（b）、（c）和（d）、（e）、（f）依次EikoNet 走時計算相對誤差和FMM 走時計算相對誤差的X-Y、X-Z 和Y-Z 切片F(xiàn)ig. 4 The relative error of the traveltime calculation of the homogeneous model. (a), (b), (c) and (d), (e), (f) represent the X-Y, X-Z and Y-Z slices of the relative error of the EikoNet and the relative error of the FMM

EikoNet 神經(jīng)網(wǎng)絡(luò)使用學(xué)習(xí)率為5×10的Adam優(yōu)化算法，殘差元數(shù)量為10，訓(xùn)練集樣本數(shù)量為10，批尺寸為752，每個批次的數(shù)據(jù)尺寸為[752, 6]，每次學(xué)習(xí)迭代200 次，此外，針對每個速度模型都需要進行一次訓(xùn)練.

（1）均勻速度模型

設(shè)置模型的速度為5 km/s，震源放置在（10 km, 10 km, 0.1 km），EikoNet 方法的走時計算結(jié)果如圖3 所示. 由圖4 可以看出EikoNet 方法的走時計算相對誤差基本為0，而FMM 則因為求解的是未經(jīng)因式分解的程函方程，故而存在點源奇異性造成了誤差值達到50%左右的近源誤差，其相對誤差平均值為0.34%. 此外，EikoNet 方法和FMM 的計算耗時分別為5.97 s 和27.53 s，效率比為4.61.

（2）塊狀速度模型

圖 5 EikoNet 方法的走時計算結(jié)果. （a）、（b）、（c）和（d）、（e）、（f）依次代表速度模型和走時的X-Y、X-Z 和Y-Z 切片F(xiàn)ig. 5 The traveltime result of the EikoNet. (a), (b), (c) and (d), (e), (f) represent the X-Y, X-Z and Y-Z slices of the velocity modeland the the travel time

圖 6 塊狀速度模型走時計算相對誤差. （a）、（b）、（c）和（d）、（e）、（f）依次代表EikoNet 走時計算相對誤差和FMM 走時計算相對誤差的X-Y、X-Z 和Y-Z 切片F(xiàn)ig. 6 The relative error of the traveltime calculation of the block model. (a), (b), (c) and (d), (e), (f) represent the X-Y, X-Z and Y-Z slices of the relative error of the EikoNet and the relative error of the FMM

設(shè)置模型的背景速度為5 km/s，在模型的正中心增加了一個大小為8 km×8 km×8 km、速度為7 km/s的塊體，震源放置在（10 km, 10 km, 10 km），EikoNet 方法的走時計算結(jié)果如圖5 所示. 如圖6 所示，EikoNet 方法出現(xiàn)了一些較小的誤差，相對誤差平均值為0.097%，而FMM 除近源誤差外，在三維空間還存在著一些數(shù)值誤差，相對誤差平均值為0.74%. EikoNet 方法和FMM 的計算耗時分別為5.99 s 和28.51 s，效率比為4.76.

（3）層狀速度模型

層狀速度模型共有5 層，每層的深度均為4 km，1—5 層的速度依次為3 km/s、4 km/s、5 km/s、6 km/s和7 km/s，震源放置在（10 km, 10 km, 10 km），EikoNet 方法的走時計算結(jié)果見圖7. 如圖8 所示，EikoNet 方法在

Z

方向速度間斷面上存在一定的誤差，相對誤差平均值為0.078%，F(xiàn)MM 除源點附近的誤差外，在其他區(qū)域內(nèi)也分布著一些比EikoNet方法更大的誤差，相對誤差平均值為0.84%. EikoNet方法和FMM 的計算耗時分別為6.09 s 和30.64 s，效率比為5.03.

圖 7 EikoNet 方法的走時計算結(jié)果. （a）、（b）、（c）和（d）、（e）、（f）依次代表速度模型和走時的X-Y、X-Z 和Y-Z 切片F(xiàn)ig. 7 The traveltime result of the EikoNet. (a), (b), (c) and (d), (e), (f) represent the X-Y, X-Z and Y-Z slices of the velocity model and the travel time

圖 8 層狀速度模型走時計算相對誤差. （a）、（b）、（c）和（d）、（e）、（f）依次代表EikoNet 走時計算相對誤差和FMM 走時計算相對誤差的X-Y、X-Z 和Y-Z 切片F(xiàn)ig. 8 The relative error of the traveltime calculation of the layered model. (a), (b), (c) and (d), (e), (f) represent the X-Y, X-Z and Y-Z slices ofthe relative error of the EikoNet and the relative error of the FMM

（4）棋盤格速度模型

設(shè)置模型的背景速度為5 km/s，在其中增加了速度擾動，擾動的振幅為1 km/s，模型的速度在4～6 km/s 之間變化，震源放置在（10 km, 10 km,0.1 km），圖9 以二維切片的形式展示了EikoNet的走時計算結(jié)果. 從圖10 可以看出，即便在速度正負變化較多棋盤格速度模型中，EikoNet 方法也較好地控制了誤差，相對誤差平均值為0.31%，F(xiàn)MM 誤差分布則與前幾個模型類似，相對誤差平均值為0.58%. EikoNet 方法和FMM 的計算耗時分別為6.06 s 和31.70 s，效率比為5.23.

3 結(jié) 論

本文引入了基于深度神經(jīng)網(wǎng)絡(luò)進行程函方程走時計算的EikoNet 方法，通過速度模型實驗，將其與快速推進法（FMM）進行了計算精度和效率的對比，結(jié)果表明：

圖 9 EikoNet 方法的走時計算結(jié)果.（a）、（b）、（c）和（d）、（e）、（f）依次代表速度模型和走時的X-Y、X-Z 和YZ 切片F(xiàn)ig. 9 The traveltime result of the EikoNet. (a), (b), (c) and (d), (e), (f) represent the X-Y, X-Z and Y-Z slices of the velocity model and the travel time

圖 10 棋盤格速度模型走時計算相對誤差. （a）、（b）、（c）和（d）、（e）、（f）依次EikoNet 走時計算相對誤差和FMM 走時計算相對誤差的X-Y、X-Z 和Y-Z 切片F(xiàn)ig. 10 The relative error of the traveltime calculation of the checkerboard model. (a), (b), (c) and (d), (e), (f) represent the X-Y, X-Z and Y-Z slices ofthe relative error of the EikoNet and the relative error of the FMM from top to bottom

（1）EikoNet 方法可以無網(wǎng)格地快速確定三維空間內(nèi)任意兩點之間的走時，計算速度是FMM 的4～5 倍，計算精度也相對更高.

（2）EikoNet 方法的訓(xùn)練和預(yù)測都是大規(guī)模并行的，計算走時只需要加載訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)參數(shù)即可，無需存儲模型的速度文件，大大降低了內(nèi)存空間.

但是，文中僅使用EikoNet 方法在幾種簡單的速度模型上進行了嘗試，在起伏界面等較為復(fù)雜的速度模型上的適用性還需要進一步研究.

致謝

感謝中國科學(xué)院地質(zhì)與地球物理研究所蘭海強副研究員和周小樂博士提供的因式分解程函方程走時計算結(jié)果. 感謝中國地震局地球物理研究所李永華研究員、張風(fēng)雪研究員以及劉有山副研究員的建設(shè)性意見. 感謝三位匿名審稿人的寶貴建議，對稿件質(zhì)量提升幫助很大.

附中文參考文獻

韓衛(wèi)雪，周亞同，池越. 2018. 基于深度學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)的地震數(shù)據(jù)隨機噪聲去除[J]. 石油物探，57（6）：862-869.

蘭海強，張智，徐濤，等. 2012. 地震波走時場模擬的快速推進法和快速掃描法比較研究[J]. 地球物理學(xué)進展，27（5）：1863-1870.

李飛，徐濤，武振波，等. 2013. 三維非均勻地質(zhì)模型中的逐段迭代射線追蹤[J]. 地球物理學(xué)報，56（10）：3514-3522.

林吉焱，唐國彬，徐濤，等. 2020. 欽杭—武夷山成礦帶上地殼速度結(jié)構(gòu)與基底特征: 萬載—惠安寬角反射/折射地震剖面約束[J].地球物理學(xué)報，63（12）：4396-4409.

劉國峰，劉洪，李博，等. 2009. 山地地震資料疊前時間偏移方法及其GPU 實現(xiàn)[J]. 地球物理學(xué)報，52（12）：3101-3108.

徐濤，徐果明，高爾根，等. 2004. 三維復(fù)雜介質(zhì)的塊狀建模和試射射線追蹤[J]. 地球物理學(xué)報，47（6）：1118-1126.

徐濤，張明輝，田小波，等. 2014. 麗江—清鎮(zhèn)剖面上地殼速度結(jié)構(gòu)及其與魯?shù)?p>M

6.5 級地震孕震環(huán)境的關(guān)系[J]. 地球物理學(xué)報，57（9）：3069-3079.

徐濤，張忠杰，劉寶峰，等. 2015. 峨眉山大火成巖省地殼速度結(jié)構(gòu)與古地幔柱活動遺跡: 來自麗江—清鎮(zhèn)寬角地震資料的約束[J].中國科學(xué)：地球科學(xué)，45（5）：561-576.

楊旭，李永華，蓋增喜. 2021. 機器學(xué)習(xí)在地震學(xué)中的應(yīng)用進展[J].地球與行星物理論評，52（1）：76-88.

張明輝, 劉有山, 侯爵, 等. 2019. 近地表地震層析成像方法綜述[J].地球物理學(xué)進展, 34(1): 54-69.

趙明，陳石，Yuen D. 2019. 基于深度學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)的地震波形自動分類與識別[J]. 地球物理學(xué)報，62（1）：374-382.

周小樂，蘭海強，陳凌，等. 2020. 曲線坐標系因式分解程函方程及其走時計算[J]. 地球物理學(xué)報，63（2）：638-651.