基于深度學習的接收函數橫波速度預測

楊庭威，曹丹平*，杜南樵，崔榮昂，南方舟，徐亞，梁策

1 中國石油大學(華東)，地球科學與技術學院，青島 266580 2 中國科學院地質與地球物理研究所，中國科學院油氣資源研究重點實驗室，北京 100029 3 中國科學院大學，北京 100049 4 中國科學院地球科學研究院，北京 100029

0 引言

自Phinney(1964)提出接收函數的基本概念以來，接收函數方法便成為研究地球內部圈層間斷面、速度結構的重要方法之一(Phinney,1964;Ammon et al.,1990;Leahy et al.,2012;Shibutani et al.,1996;Zhu and Kanamori,2000).Langston(1977)通過對徑向分量和垂直分量的譜分解(反褶積)，從遠震事件的長周期體波中提取了接收函數(Langston,1977).Owens等(1984)將接收函數的方法擴展到寬頻，并利用遠震P波接收函數的波形擬合方法來反演地殼速度，發展了基于遠震接收函數的線性化反演方法(Owens et al.,1984).盡管線性化方法在某些情況下得到的地殼模型可以很好地擬合觀測到的接收函數波形，但結果通常對初始模型有很強的依賴性，常常收斂到局部最優解而非全局最優解(Sambridge,1999a；任駿聲和沈旭章，2015).相比之下，基于馬爾科夫鏈的蒙特卡羅方法可以避免局部極小值，逼近全局最優解(Sambridge,1999b;Vinnik et al.,2004;Zheng et al.,2005).但該方法需要通過大量的正演計算來采樣反問題的后驗概率密度分布，使其在實際應用中受到了一定限制.

線性化反演方法和蒙特卡羅方法都是模型驅動的，遠震接收函數反演的數據量較大，每次迭代都需要進行大量正演計算.而利用深度學習從大量的遠震接收函數數據中自動提取其復雜的特征，預測臺站下方橫波速度，避免全局反演解的不穩定性，并克服傳統方法中擬合接收函數迭代運算的大運算量缺陷，大幅提升了反演的速度、精度和效率.其中卷積神經網絡通過共享權重，設置相對較少的參數即可構建層數更深的網絡，其提取特征的能力較以往的人工神經網絡有很大提升(劉帥師等,2016;Krizhevsky et al.,2017;安鵬和曹丹平,2018;Puzyrev,2019;Zhu and Beroza,2019)，非常適用于地震波識別與定位以及接收函數反演等這類特征提取與參數優化的問題(Van der Baan and Jutten,2000;Kong et al.,2016;Ross et al.,2018；段剛，2021).在反演中只需將徑向接收函數作為輸入參數，經過卷積層、下采樣層和全連接層的前向傳播和反向傳播，網絡訓練結束后便能直接得到橫波速度，這對于反演巖石圈速度從模型驅動、方程反演朝著自動化、智能化方向發展具有十分重要的意義.

本文基于卷積神經網絡，設計了利用接收函數預測橫波速度的方法.在研究中利用全球模型參數的合成接收函數與高質量臺站觀測接收函數，建立一個與全球地殼結構復雜程度高度一致的樣本庫集合，提高卷積神經網絡的泛化能力.最后，我們將訓練好的神經網絡應用到被動源海底地震儀觀測數據上，對琉球海溝地區的地殼結構進行深入研究.

1 卷積神經網絡預測方法

基于卷積神經網絡對接收函數波形進行訓練并預測地下橫向波速特征，其核心是構建合理的網絡結構，建立具備泛化能力的樣本集進行測試.本文引入U-net的核心結構，其跳層的連接可以減少神經網絡中梯度爆炸或者消失的問題，并且U-net中多層次的特征信息可以從不同尺度挖掘接收函數和橫波速度之間的關系，因而選用U-net提取接收函數復雜的特征，預測臺站下方橫波速度.

基于U-net卷積神經網絡結構的基本技術框架如圖1所示，將遠震P波徑向接收函數的波形作為輸入參數，對應的橫波速度結構作為樣本標簽.左圖每個綠色盒對應一個多通道特征圖.通道的數量顯示在盒子的頂部.黃色盒代表復制的特征圖.不同顏色箭頭表示不同的操作.最先提出的U-net是用于醫學領域的細胞分割(Ronneberger et al.,2015)，輸入輸出均為二維數據，為了適應地震反演的實際問題，我們將U-net拓展為一維.可以更好適用于轉換輸入數據和輸出數據的維度，例如輸入數據是不同震中距的多道接收函數為二維數據,而輸出數據是臺站下方一維的橫波速度結構為一維數據.

圖1 基于U-Net的卷積神經網絡結構Fig.1 Architecture of U-Net-based convolutional neural network

為提高該網絡結構的特征提取能力，本文在U-net網絡上進行局部改進，將最后卷積層的輸出經過Flatten操作后傳入全連接層.加深了神經網絡的深度，也就使深層的神經元可識別更為抽象的信息，同時利用dropout操作讓某神經元的激活值以一定的概率p停止工作，防止過擬合，使模型泛化性更強.從而將反演問題簡單化，精確化.

1.1 卷積神經網絡的前向傳播

池化層的前向傳播過程是以卷積層提取的特征作為輸入傳到下采樣層，并在下采樣層進行池化操作，降低數據的維度，避免過擬合.網絡采用的池化方法為均值池化，表示取局部接受域中的平均值.

卷積神經網絡前向傳播公式為：

x=f(u)，

(1)

u=W*x+b，

(2)

1.2 卷積神經網絡的反向傳播

本文用來衡量真實標簽值與預測結果之間誤差情況的目標函數為：

(3)

其中f(ul(w,b)為網絡最后一層的輸出，V為標簽，w和b分別是權值和偏置，在Adam算法中統一用θ來表示.

本文利用Adam算法更新網絡權值(Kingma and Ba,2017)，具體過程如下:

(1)設置參數.步長ε設為0.001;矩估計的指數衰減速率ρ1和ρ2在區間[0,1)內(默認設為0.9和0.999);用于數值穩定的常數δ(默認設為10-8)；

(2)初始化參數.初始化一階和二階矩變量s=0，r=0;

2 接收函數預測橫波速度的樣本庫構建

樣本庫是深度學習的核心基礎，構建高質量的樣本庫是提高深度學習泛化能力的基礎和關鍵，本文從合成接收函數樣本及實際數據的接收函數樣本兩方面出發提出構建樣本訓練集的優化方案，建立與全球地殼結構復雜程度高度一致的樣本庫集合.

2.1 合成接收函數樣本

本文采用一維水平層狀各向同性介質的反射率方法合成理論地震波形，該技術可以用較少的耗時合成與實際數據較為接近的地震波形，再使用最大熵反褶積法從三分量波形數據中提取接收函數(Levin and Park,1997a,b;Li et al.,2010).最大熵譜反褶積是以熵極大作為互相關函數和自相關函數的外推準則，來求取接收函數.該方法在外推計算過程中，其反射系數總是小于1，以確保最大熵反褶積的穩定性，并對數據時窗外的數據做更加合理的假設，以獲得高分辨率的接收函數(吳慶舉等,2003).同一模型的合成數據包括20道徑向接收函數，以此模擬不同地震事件.正演計算采用2.5的高斯系數，在P到達之后使用32 s的時間窗，從而保證接收函數包含直達波和主要轉換波的信息.

為提高網絡的泛化能力，建立一個與全球地殼結構復雜程度高度一致的樣本庫.為此我們基于Crust1.0模型進行合成接收函數樣本的構建工作.在Crust1.0全球模型中每隔5經緯度取一個點，共2592個模型數據，不同深度的橫波速度分布如圖2所示.并將介質中的每一層細分為1 km厚的區間，用隨機高斯噪聲擾動每一薄層的橫波速度，擾動均值為0，標準差為0.1 km·s-1，利用所得結果正演射線參數分布在0.042～0.080 s·km-1內的20道接收函數(該范圍內接收函數質量較好)，再添加隨機噪聲，使其與真實接收函數更為接近.

全球地殼各處厚度不一，大陸部分平均厚度為37 km；而海洋部分平均厚度則只有約7 km.高山、高原部分地殼最厚可達70 km；裂谷、海溝部分地殼最薄低至1.5 km.樣本中包含了從1.5 km到70 km厚度的地殼模型.樣本集的地殼厚度分布如圖2g所示，其中地殼厚度大于60 km或小于3 km樣本不足，因此本文只討論地殼厚度為5～50 km范圍內時該方法的有效性.將地殼模型進行編號，用隨機數字生成器隨機選取地殼速度模型，與低速體模型、高速體模型、高低混合模型隨機組合形成408個復雜模型.基于模型參數正演射線參數分布在0.042～0.080 s·km-1內的20道接收函數，并添加隨機噪聲(圖3).

圖2 地殼橫波速度結構的水平切片(a)5 km;(b)10 km;(c)15 km;(d)25 km;(e)35 km;(f)45 km；(g)地殼模型分布(數據集來源于Crust1.0).Fig.2 Horizontal slices of crustal shear wave velocity structures at depths(a)5 km;(b)10 km;(c)15 km;(d)25 km;(e)35 km;(f)45 km；(g)Crustal model distribution (dataset from Crust1.0).

圖3 合成接收函數樣本展示(a)低速體模型;(b)合成接收函數;(c)添加隨機噪聲的接收函數.Fig.3 Sample set of synthetic receiver functions(a)Low-speed model;(b)Synthetic receiver functions;(c)Synthetic receiver functions with random noise.

2.2 觀測接收函數樣本

除合成接收函數樣本外，本文進一步基于高質量的實際地震觀測數據進行接收函數提取，并將其作為樣本集的組成之一進行訓練.在JAMSTEC網站上申請了2013年11月23日到2014年2月18日琉球海溝附近30個OBS海底地震儀的連續波形數據.為此對實際地震數據先進行預處理，包括格式轉換和地震識別.首先是把原始記錄SEED格式數據轉換成較為通用的SAC格式，然后根據哈佛大學的CMT目錄截取地震震級在5.5級以上，震中距在30°～90°的三分量地震波形;陸地地震數據方面，在IRIS網站申請美國地區10個固定臺震級在5.5級以上，震中距在 30°～90°的三分量地震波形數據.

從遠震波形中提取出射線參數在0.042～0.080 s·km-1范圍內的接收函數.每個臺站分別挑選出20道質量較高的接收函數，進行可逆跳躍馬爾可夫鏈蒙特卡羅接收函數反演臺站下方的橫波速度結構,作為網絡的標簽(圖4).樣本集兼顧了陸地和海洋不同的地震數據，以提高網絡對不同類型觀測接收函數數據進行橫波速度預測的泛化能力.

圖4 觀測接收函數樣本展示(a)、(c)可逆跳躍馬爾可夫鏈蒙特卡羅接收函數反演的橫波速度;(b)、(d)觀測接收函數.Fig.4 Sample set of observed receiver functions(a)and (c)Shear wave velocities inverted by a reversible jumping Markov chain Monte Carlo；(b)and (d)Observed receiver functions.

3 方法測試與實驗

3.1 輸入樣本的選取

隨著Crust系列全球地殼模型精度的不斷提高，當前最為詳細，分辨率最高的全球地殼模型是Crust1.0地殼模型，其地形特征的多樣性和完整性非常適用于作為神經網絡的訓練集.同時，在全球模型數據的基礎上，本研究加入了隨機生成的陸殼高原、洋殼沉積層以及高低速異常體等復雜模型以及實際的觀測數據，以提高訓練數據集的多樣性，增強模型的泛化能力.因此，將全球模型中2592個合成數據、308個復雜模型的合成數據以及16個觀測數據(包括陸地臺站8個，OBS數據8個)，共計2916個樣本作為網絡的訓練數據集.剩下100個復雜模型的合成數據以及4個觀測數據，共計104個樣本作為網絡的測試數據集，用于評估最終模型的泛化能力.

3.2 測試集數據檢驗

對網絡進行3萬次訓練學習，耗時方面，利用NVIDIA 2080 GPU訓練耗時只有10 h，這很大程度上由計算機硬件配置、學習系數的數據量所決定.用訓練好的網絡模型進行橫波速度預測，在i5四核CPU上預測橫波速度耗時僅在5 s左右.在測試集中計算四個具有代表性的模型：洋殼模型、高原模型、高速模型和低速模型來驗證方法的有效性及準確性.

圖5為隨機選取測試集中四個代表性模型的預測結果，可見四個結果預測橫波速度均與模型真實值基本一致.圖5a為洋殼模型的預測結果，黑色實線表示模型的真實橫波速度隨深度變化的曲線，紅色虛線表示利用網絡預測的橫波速度隨深度變化的曲線，可見對5 km沉積層厚度以及12 km的Moho面深度符合較好，預測誤差不超過0.2%.圖5b為陸地高原模型的預測結果，對3 km沉積層厚度以及40 km的Moho面深度符合較好，預測誤差不超過0.6%.圖5c為高速體模型的預測結果，在17 km處存在8.5 km高速體厚度，以及41 km的Moho面深度符合較好，預測誤差不超過0.5%.圖5d為低速體模型的預測結果，在8 km處存在9 km高速體厚度，以及42 km的Moho面深度符合較好，預測誤差不超過0.4%.預測結果均具有很好的符合度.

圖5 測試集中合成接收函數預測的橫波速度結果(a)洋殼模型的預測結果;(b)陸地高原模型的預測結果;(c)高速體模型的預測結果;(d)低速體模型的預測結果.Fig.5 Shear wave velocity models predicted by synthetic receiver functions in test dataset(a)Oceanic crust model；(b)Plateau model；(c)High-speed model；(d)Low-speed model.

為更好地檢驗網絡的泛化性，在測試集中加入地殼厚度為1～70 km的70個不同樣本如圖6a所示，并計算預測結果與模型均方誤差隨地殼厚度分布變化，以用于評估網絡的可靠性.圖6c測試結果顯示：在地殼厚度為5～60 km范圍內預測誤差較小，而該范圍內地殼厚度已基本能覆蓋實際觀測臺站下方的地殼厚度.由于訓練集中鮮有地殼厚度小于5 km大于60 km的樣本，因此該范圍內地殼厚度預測誤差較大.除此之外，在測試集中再加入橫波速度為0.5～3.5 km·s-1的300個不同的速度異常體樣本如圖6b所示，其預測結果與模型均方誤差隨異常體速度分布如圖6d所示：當低速異常值大于1 km·s-1時預測誤差較小，表明網絡在預測低速異常幅值不低于1 km·s-1時仍具有較高分辨率.

圖6 測試集中合成接收函數預測結果評估(a)地殼厚度為1～70 km模型；(b)異常體橫波速度為0.5～3.5 km·s-1模型；(c)預測結果與模型均方誤差隨地殼厚度分布；(d)預測結果與模型均方誤差隨異常體速度分布.Fig.6 Evaluation of synthetic receiver functions in test dataset(a)Models with crustal thickness 1～70 km；(b)Abnormal body of shear wave velocity as 0.5～3.5 km·s-1；(c)MSE of prediction versus crustal thickness；(d)MSE of prediction versus velocity of anomalous body.

訓練集中單個模型的20道不同的合成接收函數其射線參數是等間隔分布在0.042～0.080 s·km-1范圍內如圖7a所示，旨在利用神經網絡學習到不同震中距的接收函數的平均特征來預測臺站下方的橫波速度結構.而實際觀測到的單臺接收函數的震中距是隨機分布的，特別是流動觀測數據很難得到震中距規則分布的接收函數.因此在測試集中加入射線參數隨機分布在0.042～0.080 s·km-1內的20道含噪接收函數如圖7b所示.同時，為檢驗模型的穩健性，在測試集中加入1000個射線參數隨機分布的含噪接收函數樣本，預測結果與模型均方誤差如圖7d所示，可以看出其誤差均在允許范圍內.

圖7 測試集中射線參數隨機分布的合成接收函數預測結果評估(a)射線參數規則分布的合成接收函數；(b)射線參數隨機分布的合成接收函數;(c)射線參數規則分布的合成接收函數與射線參數隨機分布的合成接收函數預測結果；(d)隨機分布的接收函數樣本預測誤差分布.Fig.7 Evaluation of synthetic receiver functions with random distribution of ray parameters in test dataset(a)Synthetic receiver functions with regular distribution of ray parameters;(b)Synthetic receiver functions with random distribution of ray parameters;(c)Prediction results of synthetic receiver functions with regular distribution of ray parameters and synthetic receiver functions with random distribution of ray parameters;(d)MSE of prediction versus test dataset number.

進一步對實際數據進行測試，隨機選取測試集中OBS觀測接收函數樣本以及陸地臺站觀測接收函數樣本來檢驗網絡在預測橫波速度的泛化能力.

圖8為隨機選取測試集中海底地震儀和陸地地震臺站兩種觀測類型的接收函數樣本的預測結果;可見兩種觀測接收函數預測的橫波速度結果均與可逆跳躍馬爾可夫鏈蒙特卡羅接收函數反演的橫波速度基本一致，橫波速度間斷面所反映的莫霍面深度與H-κ疊加結果也較為一致.

圖8 測試集中觀測接收函數預測的橫波速度結果(a)OBS觀測接收函數預測的橫波速度結果；(b)OBS觀測接收函數的H-κ疊加結果;(c)陸地臺觀測接收函數預測的橫波速度結果；(d)陸地臺觀測接收函數的H-κ疊加結果.Fig.8 Shear wave velocity predicted by observed receiver functions in test set(a)By the receiver function observed on the OBS；(b)Result of H-κ stacking for receiver function observed on the OBS；(c)Shear wave velocity predicted by the receiver function observed on land station；(d)Result of H-κ for the receiver function observed on the land station.

當僅用實際樣本訓練的網絡時，網絡學習到的特征只能針對此研究區域的橫波速度預測，且由于實際樣本較少容易出現過擬合現象.如圖9所示，預測的橫波速度與實際橫波速度點距離回歸線較遠，觀測數據與合成數據預測結果均不太理想.當僅利用Crust1.0模型的合成數據訓練網絡時，網絡能較好預測合成數據和實際數據的橫波速度，相較于僅用實際訓練的預測精度有所提升.如圖10所示，預測橫波速度與實際橫波速度點相對集中于回歸線中，其中合成數據預測結果有明顯的提升.當樣本集中同時包括合成數據和實際數據訓練網絡時，網絡能學習到Crsut1.0中更豐富的模型特征，同時加入少量的實際數據訓練，能使網絡更適用于實際數據的預測，增強網絡的泛化能力.如圖11所示，預測橫波速度與實際橫波速度點基本落在回歸線上，實際數據與合成數據的預測精度均有所提升.

圖9 訓練集中僅有實際數據網絡預測的橫波速度與實際橫波速度統計(a)OBS觀測數據預測的橫波速度與實際橫波速度相關性；(b)合成數據預測的橫波速度與實際橫波速度相關性.Fig.9 Statistics of predicted shear wave velocities and actual shear wave velocity from trained set only with observed samples(a)Correlation between shear wave velocities predicted by OBS observational data and actual shear wave velocities；(b)Correlation between shear wave velocities predicted by synthetic data and actual shear wave velocities.

圖10 訓練集中僅有合成數據網絡訓練預測的橫波速度與實際橫波速度統計(a)OBS觀測數據預測的橫波速度與實際橫波速度相關性；(b)合成數據預測的橫波速度與實際橫波速度相關性.Fig.10 Predicted shear wave velocity and actual shear wave velocity statistics trained only with synthetic samples(a)Correlation between shear wave velocities predicted by OBS observation data and actual shear wave velocities;(b)Correlation between shear wave velocities predicted by synthetic data and actual shear wave velocities.

圖11 訓練集中有實際數據及合成數據網絡預測的橫波速度與實際橫波速度統計(a)OBS觀測數據預測的橫波速度與實際橫波速度相關性；(b)合成數據預測的橫波速度與實際橫波速度相關性.Fig.11 Statistics of predicted shear wave velocity and actual shear wave velocity from trained data set with observed samples and synthetic samples(a)Correlation between shear wave velocities predicted by OBS observational data and actual shear wave velocities；(b)Correlation between shear wave velocities predicted by synthetic data and actual shear wave velocities.

上述結果表明，基于全球模型數據以及部分質量較高的觀測接收函數建立樣本庫能有效提高網絡的泛化性，同時對比接收函數H-κ疊加結果驗證了利用卷積神經網絡預測的橫波速度間斷面具有較高的可靠性.

4 實際數據應用

海洋地殼結構研究能提供大洋內部結構的認識，基于海底地震儀的觀測應用為獲得海洋地殼結構提供了良好的數據基礎.特別是近年來被動源OBS觀測陣列的布設為接收函數研究提供了良好的數據支撐，其研究成果可為研究精細海洋地殼結構、海底擴張和板塊運動提供重要的地球物理依據(Akuhara et al.,2016;Audet,2016;郭衍龍等,2016;Liu et al.,2020).因此，本文選取琉球海溝地區作為研究區進行實際應用檢驗，設計利用卷積神經網絡識別被動源接收函數的波形特征來預測琉球海溝附近臺站下方橫波速度結構，并通過速度間斷面進一步獲取研究區域的莫霍面深度.

為驗證基于卷積神經網絡的接收函數橫波速度預測方法對地殼結構求取的準確性和可靠性，利用在JAMSTEC網站上申請的30個OBS臺站數據進行海底地殼結構的研究.琉球海溝是菲律賓海板塊向下俯沖到歐亞板塊的一個俯沖帶構造，具有強烈的構造活動，是研究太平洋西岸匯聚型板塊邊界的重要研究區域(Ali and Hall,1995;Hall et al.,1995;Sdrolias et al.,2004;Yumul et al.,2009;Arai et al.,2016).圖12中紅色線表示板塊的邊界，黑色線為斷層，白色箭頭表示板塊俯沖的方向，被動源OBS臺站所處位置由北向南從琉球島弧過渡到琉球海溝南部的斜坡帶處.由于3臺OBS數據無法獲取，本文僅選取研究區域27個數據質量較高的OBS，并提取20道接收函數作為輸入參數，經過訓練好的卷積神經網絡輸出橫波速度，并求取橫波速度在深度方向的梯度，從中獲得研究區域OBS下方莫霍面深度.

圖12 琉球海溝OBS地理位置Fig.12 Geographical locations of OBS in the Ryukyu Trench

為驗證卷積神經網絡結果的可靠性，結合接收函數H-κ疊加結果與Crust1.0模型進行對比，如圖13所示.圖14為利用卷積神經網絡預測研究區海域的橫波速度所計算的莫霍面深度分布圖，呈現出由北向南莫霍面深度從26 km到12 km逐漸減小的趨勢，以及地殼厚度的北厚南薄趨勢，明顯地分為三級階梯，與接收函數H-κ疊加結果以及Crust1.0莫霍面深度分布基本一致.通過對比，驗證了利用卷積神經網絡預測的橫波速度間斷面具有較高的可靠性.由于具有更高的分辨率，利用卷積神經網絡得到莫霍面深度分布能夠更為細致地刻畫研究區的地殼結構變化，同時滿足了準確度，精細度.

圖13 研究區域莫霍面深度分布(a)H-κ疊加結果莫霍面深度分布；(b)Crust1.0莫霍面深度分布圖.Fig.13 Moho depth distribution in the study area(a)Moho depth distribution for H-κ stacking；(b)Map showing Crust1.0 Moho depth distribution.

圖14 基于卷積神經網絡預測的橫波速度所計算的莫霍面深度分布圖Fig.14 Moho depth distribution based on the shear wave velocity predicted by convolutional neural network

5 結論

基于卷積神經網絡方法本文設計了一種接收函數橫波速度的預測方法，利用深度學習從大量的遠震接收函數數據中自動提取其復雜的特征，預測臺站下方橫波速度，避免全局反演解的不穩定性，并克服傳統方法中擬合接收函數迭代運算的大運算量缺陷，大幅提升了反演的速度、精度和效率.本研究基于Crust1.0模型構建模型接收函數樣本，基于陸地及海底地震高質量觀測實際數據構建實際接收函數樣本，建立與全球地殼結構復雜程度高度一致的樣本庫集合，增強了網絡學習的泛化能力.經過測試集對模型及實際數據的檢驗均取得了良好的橫波速度預測效果.進一步將該方法用于琉球海溝地區被動源海底地震數據接收函數研究，獲得了該地區海洋地殼結構特征，結果表明該地殼莫霍面在12～26 km范圍內變化，整體莫霍面特征較Crust1.0模型更為精細，為認識琉球海溝地區深部結構提供了可靠參考.

致謝感謝審稿專家提出的寶貴意見.IRIS、Japan Agency for Marine-Earth Science and Technology為本研究提供地震數據.