瑞雷波頻散曲線的深度學(xué)習(xí)反演方法

張志厚， 石澤玉， 馬寧*， 王虎， 喬中坤， 趙思為，姚禹， 趙明浩， 葉志虎

1 西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院， 成都 611756 2 西南交通大學(xué)，高速鐵路線路工程教育部重點實驗室， 成都 610031 3 浙江工業(yè)大學(xué)理學(xué)院，浙江省量子精密測量重點實驗室， 杭州 310023 4 中鐵二院成都地勘巖土工程有限責(zé)任公司， 成都 610000

0 引言

瑞雷波是縱橫波在地面相互干涉形成并沿地表傳播的一種地震面波，其具有能量強、衰減慢、信噪比高等優(yōu)勢，被廣泛用于淺地表及中深部速度結(jié)構(gòu)的調(diào)查(Miller et al., 1999; Du et al., 2019; Tang et al., 2019; Berg et al., 2020; Wang et al., 2021)等.諸多學(xué)者分別從瑞雷波場模擬(Xu et al., 2007; Zeng et al., 2011; 胡書凡等, 2021)、頻散曲線提取(Hu et al., 2020; Xi et al., 2021; Dai et al., 2021; 易佳等, 2021)、頻散曲線反演(Xia et al., 1999; 陳祥和孫進忠, 2006; Lu et al., 2007; Picozzi and Albarello, 2007; 于東凱等, 2018; Mi and Xia，2021; Poormirzaee and Fister, 2021)以及波形反演(Pan et al., 2019; Zhang and Alkhalifah, 2019)等方面進行研究，使得瑞雷波研究與應(yīng)用得到了極大的發(fā)展.頻散曲線反演是瑞雷波勘探的重要環(huán)節(jié)之一(Lu et al., 2007)，也是獲取近地表橫波速度的關(guān)鍵步驟(Xia et al., 1999; 夏江海等, 2015).然而，同其他大多數(shù)地球物理反演優(yōu)化問題類似，頻散曲線反演具有高度的非線性、多參數(shù)及多極值等特點(于東凱等, 2018)，一般線性反演方法效果取決于初始模型的選擇和偏導(dǎo)數(shù)矩陣計算的準確性；非線性全局最優(yōu)化算法如粒子群(Song et al., 2012; Poormirzaee，2016)、遺傳(Picozzi and Albarello, 2007)、模擬退火(Pei et al., 2007; Lu et al., 2016)、布谷鳥搜索(Poormirzaee and Fister, 2021)以及多種全局優(yōu)化方法改進或結(jié)合算法(Sun et al., 2017; Lei et al., 2019)等存在收斂速度慢、精度低等問題(于東凱等, 2018).

深度學(xué)習(xí)(Deep Learning, DL)是利用各類型的深層次神經(jīng)網(wǎng)絡(luò)模型來模擬人腦結(jié)構(gòu)，該過程是信息的分布式存儲和并行協(xié)同處理一個非線性系統(tǒng)的過程，盡管每個神經(jīng)元或網(wǎng)絡(luò)層的結(jié)構(gòu)和功能非常簡單，但其組合時卻能完成復(fù)雜的非線性映射.深度學(xué)習(xí)在過去10年間為數(shù)據(jù)科學(xué)領(lǐng)域帶來了一場革命，如在圖像分割(Krizhevsky et al., 2012)、目標檢測(Ren et al., 2017)、語音識別(Hinton et al., 2012)、機器翻譯(Sutskever et al., 2014)、情感分類(Zhao et al., 2019)以及圍棋游戲(Silver et al., 2016)等取得了突破性進展，使許多之前被認為無法實現(xiàn)的任務(wù)成為了可能(Bronstein et al., 2021).由于DL輸入到輸出的復(fù)雜映射關(guān)系難以明確解構(gòu)，因此其屬于典型“黑盒”方法.相比統(tǒng)計學(xué)上大量使用的“數(shù)據(jù)模型”，DL網(wǎng)絡(luò)架構(gòu)聚焦于尋找一種預(yù)測的算法，因此被稱為“算法模型”(Diaconis and Mosteller, 1989)，或稱為一種“網(wǎng)絡(luò)結(jié)構(gòu)映射”，其本質(zhì)是一個函數(shù)估計問題：在訓(xùn)練集上給定某些未知函數(shù)的輸出，試圖從某些假設(shè)函數(shù)類別中找到一個函數(shù)F，該函數(shù)可以很好地擬合訓(xùn)練數(shù)據(jù)，并可以預(yù)測出先前未見過的輸入對應(yīng)的輸出(Bronstein et al., 2021).并且根據(jù)萬能近似定理(Hornik et al., 1989)，DL可以以理想的精度近似任意復(fù)雜非線性映射函數(shù)(LeCun et al., 2015)，有助于解決地球物理正反演問題(Jin et al., 2017)，以及DL與地球物理反演方法都具有類似的優(yōu)化算法(Kim and Nakata, 2018).DL的顯著優(yōu)勢也促使了諸多學(xué)者將其引入到地球物理領(lǐng)域中，無論是對原始數(shù)據(jù)進行預(yù)處理(Wang et al., 2019; Xu et al., 2020; 俞若水等, 2020; Dai et al., 2021; Guo, 2021；諸峰等，2022),還是對勘探成果的自動化解譯(Wu et al., 2019; Xiong et al., 2018; Shi et al., 2019; 何易龍等，2022；張志厚等，2022)，以及在地球物理“端到端”的反演中(Puzyrev, 2019; Liu et al., 2020; Li et al., 2020a,b; 張志厚等, 2021a，b)等都取得了較好效果.

DL除去利用已訓(xùn)練好的網(wǎng)絡(luò)模型進行反演預(yù)測本身之外，其主要包括樣本數(shù)據(jù)集的構(gòu)建和網(wǎng)絡(luò)結(jié)構(gòu)的設(shè)計兩大重要內(nèi)容.

在數(shù)據(jù)集構(gòu)建方面：DL樣本數(shù)據(jù)主要來源于真實數(shù)據(jù)和人工合成數(shù)據(jù).相關(guān)學(xué)者探索了基于精細化標注真實數(shù)據(jù)訓(xùn)練的語義分割模型(Bai et al., 2018)，然而真實樣本數(shù)據(jù)標注耗時、昂貴，且對人為錯誤敏感，并不能滿足地球物理反演的需求；人工合成方法可以獲取大量的樣本數(shù)據(jù)，足量代表性數(shù)據(jù)可極大提高訓(xùn)練模型的泛化性能，如在DL的重磁反演(張志厚等, 2021a，b)、電磁反演(Puzyrev, 2019; Li et al., 2020b)、直流電阻率反演(Liu et al., 2020)、地震反演(Li et al., 2020a)中取得了較理想的成果，但以上內(nèi)容缺乏樣本數(shù)據(jù)集構(gòu)建方法的系統(tǒng)性分析.

在網(wǎng)絡(luò)結(jié)構(gòu)精心設(shè)計方面：近10年來，以ImageNet為代表的大型、高質(zhì)量數(shù)據(jù)集和不增長的計算資源(GPU)使我們能夠設(shè)計各種被用于此類大型數(shù)據(jù)集的網(wǎng)絡(luò)結(jié)構(gòu)(Bronstein et al., 2021)，如卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Networks，CNN)(Krizhevsky et al., 2012)、深度信念網(wǎng)絡(luò)(Deep Belief Networks, DBN)(Hinton et al., 2006)、對抗神經(jīng)網(wǎng)絡(luò)(Generative Adversarial Nets, GAN)(Goodfellow et al., 2014)、長短時記憶網(wǎng)絡(luò)(Long Short-term Memory Nets, LSTM)(Hochreiter and Schmidhuber, 1997)等.雖然很難理解以上每種網(wǎng)絡(luò)結(jié)構(gòu)之間存在的關(guān)系，但有學(xué)者提出了幾何深度學(xué)習(xí)的概念(Bronstein et al., 2017)，從而嘗試將以上網(wǎng)絡(luò)進行統(tǒng)一，并利用對稱先驗和尺度分離原理推導(dǎo)出當(dāng)下較為典型的網(wǎng)絡(luò)架構(gòu)，如由平移對稱推導(dǎo)出了CNN，由置換不變性導(dǎo)出的圖神經(jīng)網(wǎng)絡(luò)(Graph Neural Networks, GNN)、DeepSets和Transformer網(wǎng)絡(luò)，由時間扭曲不變形導(dǎo)出的LSTM，以及由規(guī)范對稱性導(dǎo)出的計算機圖形和視覺中使用的Intrinsic Mesh CNN等.瑞雷波不同頻率相速度反映不同深度的橫波速度，表明其具有一定的空間對應(yīng)性，頻散數(shù)據(jù)本質(zhì)也是一種時間序列的變換，說明了瑞雷波頻散數(shù)據(jù)同時具備空間對應(yīng)性和時間序列性.而CNN可以通過卷積算子提取頻散曲線的空間特征，LSTM用于頻散序列數(shù)據(jù)的預(yù)測，且LSTM在處理梯度爆炸與消失問題方面相比其他網(wǎng)絡(luò)結(jié)構(gòu)具有更好的優(yōu)勢(Kong et al., 2019).

基于此，本文首先針對DL所依賴的樣本數(shù)據(jù)集的足量多樣屬性，展開了淺地表速度結(jié)構(gòu)所具備遍歷性和有序性的討論與分析，并由此提出了約束馬爾科夫決策的樣本數(shù)據(jù)構(gòu)建方法；其次，基于不同網(wǎng)絡(luò)結(jié)構(gòu)的各自優(yōu)勢，針對瑞雷波頻散曲線具備的空間性和序列性，提出了一種CNN銜接LSTM的DL網(wǎng)絡(luò)結(jié)構(gòu)(CNN-LSTM)用于瑞雷波頻散曲線反演；最后，對本文所提方法進行理論數(shù)據(jù)和實測數(shù)據(jù)的檢驗，以此來證實文中方法的可行性、有效性和適用性.

1 反演樣本數(shù)據(jù)集構(gòu)建

近地表速度結(jié)構(gòu)數(shù)據(jù)和速度結(jié)構(gòu)一一對應(yīng)的頻散曲線數(shù)據(jù)組成了樣本數(shù)據(jù)對.訓(xùn)練樣本數(shù)據(jù)的多樣性和足量性與DL方法所具備的泛化性能密切相關(guān)，其決定方法的實用性，然而采用單純的遍歷方法所構(gòu)建的近地表速度結(jié)構(gòu)模型一是計算成本高，二是不具有真實地下速度結(jié)構(gòu)的典型性，三是造成部分大量樣本數(shù)據(jù)的價值密度低，四是無法通過樸素的方式進行學(xué)習(xí).因此，近地表速度結(jié)構(gòu)模型的科學(xué)構(gòu)建顯現(xiàn)的尤為重要.

1.1 近地表速度結(jié)構(gòu)遍歷假說

“世上沒有兩片完全相同的樹葉”，看似簡單的自然現(xiàn)象，竟然可以引起人們對宇宙萬物的深思.46億年地球演化過程從沒有簡單的重復(fù)，物種的多樣性，自然的多元性，都在隨時隨地證明自然系統(tǒng)的各態(tài)遍歷屬性，使得遍歷假說也許是打開非線性科學(xué)的一把鑰匙(楊文采, 2008).毫無疑問，近地表地層速度結(jié)構(gòu)作為地球演化的表層系統(tǒng)也具有其內(nèi)在遍歷性.伴隨著地球的演化過程，近地表速度結(jié)構(gòu)的演化軌道也將力圖遍歷其相空間各種可能的體制和狀態(tài).

各態(tài)遍歷假說闡述了一個“時間平均”與“空間平均”的等價關(guān)系，即在一個充分大空間內(nèi)，同一區(qū)域不同時間上的各個可能狀態(tài)遍歷了同一時間不同區(qū)域的所有可能狀態(tài)，反之亦然(Zhan, 1999).那么通過地球演化視角來看，同一區(qū)域近地表地層速度結(jié)構(gòu)最終覆蓋其全部的相空間，可以通過同一時間不同區(qū)域的速度結(jié)構(gòu)狀態(tài)來替代.因此，可以利用不同區(qū)域的速度結(jié)構(gòu)來組成深度學(xué)習(xí)所需的樣本數(shù)據(jù)集，遍歷速度結(jié)構(gòu)可用式(1)表達：

vk=random(vmin,vmax)，k=1,2,3,…,n，

(1)

式中，vk表示第k層的速度，n表示一維層數(shù)，random為隨機采樣函數(shù)，vmin與vmax為隨機采樣空間的界限.

1.2 速度結(jié)構(gòu)演化特征的有序性

自然系統(tǒng)在演化中還有力圖保持自身特征形態(tài)的慣性(楊文采, 2008)，使其隨時空變化有序.如向日葵、松果籽、鳳梨呈現(xiàn)出的斐波那契數(shù)列；斑馬的皮毛、蝴蝶的紋理以及花朵的形狀展現(xiàn)出規(guī)則的圖案，都在隨時隨地驗證自然的有序性.自然過程的慣性也是大自然的最佳選擇，如河流演變是在流域內(nèi)既定的地形、植被與地質(zhì)等條件下，歷經(jīng)千萬年的時空隨機降水，通過地上地下的匯聚而形成的，其流動能耗最低.綿延的群山，地層的形變，巖石的自然裂縫，同一區(qū)域不同時間的近地表速度結(jié)構(gòu)等也可以認為是大自然的慣性使然.對于近地表的層狀速度結(jié)構(gòu)，其慣性特點為速度隨深度的增加而增大(圖1e)，這是由于地層隨著深度的增加，其自重應(yīng)力也隨著增加，在力的作用下，其速度也將發(fā)生變化，可以用式(2)表述：

vk+1=vk+Δvk+1∝γkρkghk，

(2)

式中ρk、hk分別為第k層的密度與層厚，g為重力加速度，γk為大于零的系數(shù).而不會出現(xiàn)隨深度增加而減少(圖1b).即便是針對不同巖性的情況也很少出現(xiàn)圖1a所示的震蕩速度結(jié)構(gòu)，但極有可能出現(xiàn)圖1c、d所示的速度結(jié)構(gòu)，并且類似于圖1e、f所示的速度結(jié)構(gòu)較為普遍.

1.3 約束馬爾科夫決策的速度模型

基于近地表速度結(jié)構(gòu)的遍歷性與有序性，文中將馬爾科夫鏈引入到速度建模中.為了便于表述，將各層狀模型速度隨深度由淺到深的變化看成時間譜序列，每一層對應(yīng)一個時間，每一層的速度對應(yīng)該時間的狀態(tài).速度結(jié)構(gòu)由淺到深建模過程視為一個隨機過程{v(t),t∈T}，其對隨深度層增加的有限時序t1

P{v(tn)=vn|v(tn-1),…,v(t1)}

=P{v(tn)=vn|v(tn-1)}.

(3)

這種隨深度增加而變化的速度結(jié)構(gòu)特性稱為速度層的馬爾科夫過程無后效應(yīng)，具有該性質(zhì)的v(t)稱為馬爾科夫過程.{v(t),t∈T}為馬爾科夫鏈(Markov Chain, MC)，它表明當(dāng)前地層的速度狀態(tài)僅依賴于上一層速度，且下一層速度僅依賴于當(dāng)前層的速度.

在條件概率P{v(tn)=j|v(tn-1)=i}中，v(tn)=j表示在tn時刻過程處于狀態(tài)j，故條件概率P{v(tn)=j|v(tn-1)=i}表示過程在時刻tn-1處于狀態(tài)i條件下，時刻tn過程變化到狀態(tài)j的概率(Ching and Ng, 2006).它相當(dāng)于隨機移動質(zhì)點在時刻tn-1的狀態(tài)為i，經(jīng)過一步隨機移動變化到狀態(tài)為j的概率，并記該條件下的概率為Pi,j:

Pi,j=P{v(tn)=j|v(tn-1)=i}，

(4)

圖1 速度結(jié)構(gòu)模型示意圖Fig.1 Schematic diagram of the velocity structure model

(5)

稱為此過程的一步轉(zhuǎn)移概率矩陣.

由于淺地表速度結(jié)構(gòu)普遍具有底層速度大于上層速度(圖1e、f所示)的特征，即下一層速度大于當(dāng)前層速度為大概率事件.劉紅帥等(2010)統(tǒng)計了不同地區(qū)8769個點的土層剪切波速度，結(jié)果表明隨著深度增加，剪切波速呈近線性增大的關(guān)系.因此，本文速度結(jié)構(gòu)建模的MC如圖2所示，下層波速大于當(dāng)前層波速并呈一定近線性增大為大概率事件，相比當(dāng)前層波速下層波速顯著增大以及減小或緩增為小概率事件.

圖2 速度結(jié)構(gòu)建模的馬爾科夫鏈示意圖Fig.2 Schematic diagram of Markov Chain for velocity

1.4 速度模型構(gòu)建的基本步驟

下面給出約束馬爾科夫決策速度模型構(gòu)建的基本步驟：

步驟1：本文將待正反演地層劃分為10層，并令第一層的厚度h1=2 m，隨后各層的厚度隨深度的增加而遞增，表達式為：

hk+1=1.10×hk,k=1,2,3,…,9，

(6)

步驟3：在約束速度范圍內(nèi)隨機選擇首層橫波波速：

(7)

方案一：

約束范圍內(nèi)隨機選擇第k+1層橫波波速與第k層橫波波速的近線性增加比：

(8)

(9)

方案二：

約束范圍內(nèi)隨機選第k+1層的橫波波速：

(10)

其中，min(a,b)為取a和b的最小值函數(shù).

方案三：

約束范圍內(nèi)隨機選第k+1層的橫波波速：

(11)

其中，max(a,b)為取a和b的最大值函數(shù).

步驟5：根據(jù)文獻(Tang et al., 2020)計算各層近似縱波波速：

(12)

并依據(jù)文獻(Ikeda et al., 2017)計算各層的密度：

(13)

其中，k=1,2,3,…,10,half.

1.5 樣本數(shù)據(jù)隨機抽樣實例

速度模型頻散曲線正演計算方法的選擇是人工合成高效、準確性樣本數(shù)據(jù)集的主要內(nèi)容，其最終決定了DL反演方法的實用性.瑞雷波頻散曲線的正演方法有Thomson-Haskell法(Haskell，1953)、δ矩陣法(Pestel and Leckie, 1963)、Schwab-Knopoff 法(Schwab and Knopoff, 1972)以及快速矢量傳遞法(凡友華等，2002)等.Thomson-Haskell法在高頻段會出現(xiàn)數(shù)值溢出等問題；Schwab和Knopoff(1972)采用行列變換使頻散函數(shù)公式更為簡潔，獲得了瑞雷波頻散曲線快速正演算法；快速矢量傳遞算法(凡友華等，2002)解決了高頻數(shù)值精度丟失等問題.Dunkin(1965)將系數(shù)傳遞矩陣方法應(yīng)用于瑞雷波頻散方程的推算過程，并克服了數(shù)值精度缺失等問題，且該算法在高頻段計算穩(wěn)定.考慮本文研究的正演頻散曲線用于淺表層速度結(jié)構(gòu)的反演，對應(yīng)頻散曲線的高頻相速度，因此文中采用Dunkin(1965)方法計算瑞雷波頻散曲線.

本文將速度模型和其對應(yīng)的頻散曲線組成樣本數(shù)據(jù)對，通過約束MC構(gòu)建大量的速度模型，從而獲得本文監(jiān)督學(xué)習(xí)所需的樣本數(shù)據(jù)集，其中，τ等于0.8和0.9時分別創(chuàng)建了4萬和10萬個樣本數(shù)據(jù)，總共包含14萬個樣本數(shù)據(jù)對.隨機抽樣了12組樣本數(shù)據(jù)如圖3所示，可以看出速度模型中包含速度層的正常遞增、突然變大、突然變小、某一層變大或變小的情況.表明通過本文提出的約束MC決策過程可以構(gòu)建具有一定代表性的樣本數(shù)據(jù)集.

圖3 抽樣樣本數(shù)據(jù)對(a) 橫波速度模型； (b) 頻散曲線.Fig.3 Randomly selected sampling data(a) Transversal wave velocity model; (b) Dispersion curve.

圖4 速度概率密度圖Fig.4 Velocity probability density map

圖4為統(tǒng)計了14萬個樣本數(shù)據(jù)中各層速度區(qū)段(速度區(qū)段范圍為10 m·s-1)的分布概率，可以發(fā)現(xiàn)：(1)概率密度較大的速度段整體呈現(xiàn)一帶狀分布區(qū)，整體上橫波速度隨深度的增加而近線性增大，說明這種情況下的速度結(jié)構(gòu)占主體優(yōu)勢，符合速度結(jié)構(gòu)自然“慣性”特點；(2)存在速度突然變大及突然變小的情況，其速度分布分別位于兩個“淺藍色”的帶狀區(qū)，符合近地表速度層的突變特點；(3)首層速度為一約束范圍的隨機取值，不存在突變情況.并且在約束范圍內(nèi)速度分布較為均一，說明利用“擲骰子”方法可以有效選取約束范圍具有遍歷屬性的速度.綜上所述，表明約束MC的樣本數(shù)據(jù)集構(gòu)建方法具有一定的有效性.

2 瑞雷波頻散曲線反演方法

2.1 頻散曲線反演網(wǎng)絡(luò)總體結(jié)構(gòu)

本文提出了一種名為CNN-LSTM的網(wǎng)絡(luò)結(jié)構(gòu)用于瑞雷波頻散曲線的反演，其流程圖如圖5所示，該網(wǎng)絡(luò)包含了輸入(Input)、編碼(Encoder)、解碼(Decoder)及輸出(Output)4個部分，其中編碼部分包含了3個局部特征學(xué)習(xí)模塊(Local Feature Learning Block, LFLB)和1個LSTM模塊，解碼部分是1個全連接層.該網(wǎng)絡(luò)旨在從頻散曲線中學(xué)習(xí)深層次特征，因此，每個LFLB模塊的卷積層和池化層都是一維的.隨著網(wǎng)絡(luò)深度的增加，卷積核通道數(shù)分別為64、128、256；該網(wǎng)絡(luò)結(jié)構(gòu)的頂層為softmax分類器，用于近地表速度結(jié)構(gòu)的分類.

圖5 一維CNN-LSTM網(wǎng)絡(luò)總體架構(gòu)圖LFLB為編碼部分的局部特征學(xué)習(xí)模塊，LSTM為編碼階段的全局學(xué)習(xí)模塊，64、128和256表示卷積通道數(shù)；Dispersion Curves為頻散曲線；Fully Connection為全連接；Shear Wave Velocity為剪切波速.Fig.5 One-dimensional CNN-LSTM network structureLFLB is the local feature learning module of the coding part, LSTM is the global feature learning module of the coding part, The number of convolution channels is 64, 128 and 256.

當(dāng)一維向量表示的頻散曲線輸入到網(wǎng)絡(luò)時，LFLB可以學(xué)習(xí)到局部特征.這些特征重構(gòu)后，從最后一個LFLB模塊輸出再輸入到LSTM層中，然后從這些輸入的局部層次特征中學(xué)習(xí)全局信息的相關(guān)性.因此LSTM的輸出包含了局部特征信息和全局信息的相關(guān)性.進而再將學(xué)習(xí)到的特征傳遞到全連接層，最后再輸入到分類器中.

2.2 頻散曲線局部特征學(xué)習(xí)模塊(LFLB)

用LFLB來學(xué)習(xí)頻散曲線的局部特征，每個LFLB模塊包含了兩個CNN層，兩個批處理歸一化層(Batch Normalization, BN)，兩個線性激活函數(shù)層(ReLU)，一個最大池化層(Max Pooling)，如圖6所示.LFLB的核心是卷積和池化，其中卷積算子具有權(quán)值共享性和局部空間連通性(Shi et al., 2019)，重點學(xué)習(xí)頻散曲線與速度結(jié)構(gòu)之間的非線性關(guān)系.BN層對每批卷積層激活結(jié)果進行歸一化處理，可提高深層次網(wǎng)絡(luò)的優(yōu)良性和魯棒性.池化層不僅可以抑制噪聲提高模型的魯棒性，同時還可以對數(shù)據(jù)進行降維，進一步提高計算速度.最大池化是常用的非線性函數(shù)，可將輸入特征劃分為不重疊區(qū)域，提取每一子區(qū)域的最大特征(Ciregan et al., 2012).ReLU進行非線性關(guān)系的激活，其主要優(yōu)勢是稀疏性，能夠減少梯度消失的可能性，還具有非常好的收斂性能和較低的計算成本(張志厚等，2021a).

圖6 局部特征學(xué)習(xí)模塊框圖Convolution為卷積算子，其大小為3，步長為1，跨度為1；BN為批處理歸一化算子；ReLU為線性整流函數(shù)；Max Pooling為一維最大池化算子，其大小為3，步長為1.Fig.6 Block diagram of local feature learning moduleConvolution operator size is 3, step is 1, span is 1; BN is the batch normalization operator;ReLU is a linear rectification function; The size of one-dimensional maximal pooling operator is 3, step is 1.

局部特征學(xué)習(xí)可以根據(jù)不同的任務(wù)進行不同的組合，主要體現(xiàn)在卷積和池化的各種參數(shù)上.卷積層主要對局部特征進行提取，它會通過輸入數(shù)據(jù)的長度與卷積內(nèi)核進行卷積，然后通過計算核元素與輸入元素之間的點積來生成特征映射.其中，一維卷積計算表達式為：

(14)

式中，x(n)為輸入信號，z(n)為輸出結(jié)果，w(n)長度為l的一維卷積核，文中卷積核初始化采用隨機方式.

然后將卷積特征輸入到BN層，再輸入到激活函數(shù)中.BN是一種數(shù)據(jù)的變換計算，其主要作用將特征的均值歸一到0，將方差歸一到1，這樣有助于加快模型的收斂速度，減小“梯度彌散”.當(dāng)歸一化特征之后再輸入到激活函數(shù)層后，輸出特征可以表示為：

(15)

(16)

其中，E[x]與D[x]分別為x的均值和方差，ε是很小的正數(shù)，α是尺度因子，β是平移因子.ReLU(x)為線性激活函數(shù)，表達式為：

ReLU(x)=max(x,0).

(17)

上述卷積、批處理歸一化和線性激活執(zhí)行兩次之后，再將特征輸入到最大池化層.池化執(zhí)行信息的非線性降采樣功能，同時也可以增加網(wǎng)絡(luò)結(jié)構(gòu)的感受野.

2.3 頻散曲線全局特征學(xué)習(xí)模塊

LSTM是一種循環(huán)神經(jīng)網(wǎng)絡(luò)體系結(jié)構(gòu)，可以進行全局序列信息相關(guān)性學(xué)習(xí)(Hochreiter and Schmidhuber, 1997).因此它被疊加在LFLB上，從提取的局部特征序列中學(xué)習(xí)全局信息的依賴關(guān)系.LSTM有4個類型的網(wǎng)絡(luò)組件，包含輸入門、輸出門、遺忘門和一個具有自循環(huán)連接單元(圖7).LSTM單元在每個時間步長t處的更新用等式(18)—(24)來描述(Hochreiter and Schmidhuber, 1997)：

圖7 LSTM模塊結(jié)構(gòu)圖黃色虛線框內(nèi)為細胞狀態(tài)更新(Ct)；紅色虛線框內(nèi)為遺忘門(ft)；藍色虛線框內(nèi)為輸入門(it)；綠色虛線框內(nèi)為輸出門(ot)；σg為sigmoid函數(shù)，tanh是雙曲正切函數(shù)；⊙是阿達瑪積；⊕表示相加.Fig.7 LSTM module structureThe cell status update (Ct) is in the yellow dashed box; The forgotten door (ft) is in the red dashed box; The input door (it) is in the blue dashed box; The output door (ot) is in the green dashed box; σg is the sigmoid function; tanh is the hyperbolic tangent function; ⊙ is the Hadamard product; ⊕means addition.

(18)

it=σg(Wi[ht-1,xt]+bi)，

(19)

ot=σg(Wo[ht-1,xt]+bo)，

(20)

(21)

(22)

xt=ot⊙σgCt，

(23)

ht=ot⊙tanh (Ct)，

(24)

式中，Ct為LSTM當(dāng)前細胞狀態(tài)更新值，xt為輸入值，ht為輸出值，W為參數(shù)矩陣，b為偏置向量，ft、it和ot分別為遺忘門、輸入門和輸出門變量，σg為sigmoid函數(shù)，tanh是雙曲正切函數(shù)，⊙是阿達瑪積，?表示相加，[·,·]表示兩個變量拼接.

2.4 全連接模塊及速度分類器

將LSTM模塊輸出的特征再輸入到全連接層，全連接層可以表示為：

zl=bl+zl-1·wl，

(25)

式中變量上標l和l-1分別對應(yīng)輸入和輸出的特征指標.

softmax是整個網(wǎng)絡(luò)結(jié)構(gòu)的速度分類器，它利用輸入的特征進行預(yù)測，該速度分類器是logistic回歸多分類問題的一種推廣，表達式為：

zi=∑jhjWji，

(26)

(27)

其中，zi為softmax分類器的輸入，hj為倒數(shù)第二層的激活函數(shù)，Wji為倒數(shù)第二層與分類器連接的權(quán)重.

3 模型實驗

3.1 模型訓(xùn)練

通過前文所述的方法共構(gòu)建的樣本數(shù)據(jù)為140000對.將樣本數(shù)據(jù)包含了訓(xùn)練、測試以及驗證數(shù)據(jù)，比例為18∶1∶1，即126000個數(shù)據(jù)作為訓(xùn)練，驗證數(shù)據(jù)和測試數(shù)據(jù)各有7000個.訓(xùn)練輪數(shù)可以自由設(shè)定，圖8所示為頻散曲線深度學(xué)習(xí)網(wǎng)絡(luò)的誤差(輪數(shù)為100)，可以看出，訓(xùn)練誤差與驗證誤差的變化保持同步，模型訓(xùn)練良好，未出現(xiàn)過欠擬合.

圖8 訓(xùn)練誤差及驗證誤差Fig.8 Training loss and validation loss

3.2 抽取測試數(shù)據(jù)反演結(jié)果

為了測試本文所提方法的反演效果，將測試集中的7000組數(shù)據(jù)全部進行反演，共耗時5.267 s，平均每組數(shù)據(jù)反演耗時0.75 ms.從測試集中隨機抽取了一些反演結(jié)果，如圖9所示.從圖中可以看出，CNN-LSTM反演方法無論是對正常地層遞增的速度模型(圖9a)，還是速度異常夾層(圖9d、e)，以及速度突變層(圖9b、c、f、g)都具有較好的反演效果，反演的橫波速度值和理論模型基本一致.同時，理論速度模型的頻散曲線與反演結(jié)果的頻散曲線擬合度極高(圖9h—n).為了進一步檢驗本文方法的穩(wěn)定性，將圖9所示的7組模型的頻散曲線加入5%的隨機噪聲，然后獲得其反演結(jié)果(圖10)，可以看出，反演結(jié)果能夠完全反映出地下速度的變化趨勢及其分層特性.

圖9 測試集中抽取的理論值與反演結(jié)果對比(a)—(g) 分別為7組理論模型與反演結(jié)果對比圖； (h)—(n) 分別為(a)—(g)理論模型與反演結(jié)果的頻散曲線圖. 圖中紅色與藍色曲線分別代表理論值與反演結(jié)果.Fig.9 Comparison of theoretical values and inversion results of randomly selected data pairs(a)—(g) Correspond to 7 groups of theoretical models and inversion results respectively; (h)—(n) are the dispersion curves of (a)—(g) theoretical model and inversion results respectively. The red and blue curves in the figure represent the theoretical value and inversion result respectively.

圖10 加噪聲數(shù)據(jù)的反演結(jié)果與理論模型對比(a)—(g) 分別為7組加噪聲數(shù)據(jù)反演結(jié)果與理論模型對比； (h)—(n) 分別為(a)—(g)理論模型與加噪聲反演結(jié)果的頻散曲線圖. 圖中紅色與藍色曲線分別代表理論值與反演結(jié)果.Fig.10 Comparison of inversion results and theoretical models for noise-contained data(a)—(g) Correspond to 7 groups of noise-contained data inversion results and theoretical models; (h)—(n) are the dispersion curves of (a)—(g) theoretical model and noise-contained data inversion results respectively. The red and blue curves in the figure represent the theoretical value and inversion result respectively.

為了客觀評價反演效果，本文分別計算了理論數(shù)據(jù)反演結(jié)果與含噪數(shù)據(jù)反演結(jié)果的誤差，如表1所示，無噪聲情況下，反演速度模型的相對誤差(Relative Error, RE)最大為0.09，頻散曲線擬合的平均相對誤差(Mean Relative Error, MRE)最大為0.056；含噪聲情況下，反演速度模型的相對誤差最大為0.15，頻散曲線擬合的平均相對誤差最大為0.089.可以看出，本文所提方法的計算精度高、魯棒性強.

表1 反演結(jié)果的誤差Table 1 Error of inversion results

總體而言，CNN-LSTM反演方法具有較高的精度與反演速度.另外，基于DL的反演方法，在抗噪能力上有一定的特點，并與訓(xùn)練輪數(shù)相關(guān).訓(xùn)練輪數(shù)少，相應(yīng)的反演精度相對低，但魯棒性就越強，即網(wǎng)絡(luò)的抗噪能力越強；訓(xùn)練輪數(shù)多，相應(yīng)的反演精度就高，但魯棒性會變?nèi)?，即DL泛化性能會降低.因此，如何平衡訓(xùn)練輪數(shù)與魯棒性之間的關(guān)系是一個值得深思的問題.

3.3 二維模型試驗

為了進一步測試本文所提方法的效果，采用二維邊坡模型(Mi et al., 2020)來檢驗方法的有效性.該模型簡化為一個傾斜界面(圖11a)，在25 m處界面深度為9.2 m，界面與水平方向的夾角為10°，該模型可以看作是由從左至右的兩層模型拼接而成，深度逐漸減小.首先采用多道面波分析(Multichannel Analysis of Surface Waves, MASW)數(shù)值計算方法(Mi et al., 2017)獲得該模型不同采集剖面中點的瑞雷波頻散速度(采用有限差分方法模擬13條合成多道記錄，在有限差分法中，空間網(wǎng)格大小為0.2 m，時間步長為0.02 ms，震源是一個位于自由表面20 Hz的雷克子波，延遲60 ms，模擬剖面道間距1 m，剖面長度24 m，約為最大勘探深度的2倍，檢波器道數(shù)為25，炮點位于剖面左側(cè)，最小偏移距為30 m, 炮間距2 m，采用滾動采集模式).然后利用本文所提方法對所有頻散曲線進行一維反演，隨后通過空間插值方案(Miller et al., 1999)對倒置一維模型進行網(wǎng)格化，再合成相應(yīng)的偽二維結(jié)果(圖11b),可以看出反演結(jié)果與真實界面位置具有較好的對應(yīng)性.

圖11 二維模型及反演結(jié)果(a) 模型示意圖(Mi et al., 2020)； (b) 反演橫波速度偽二維圖.圖中黑色虛線為界面位置.Fig.11 Two-dimensional model and inversion result(a) Illustration of a syntheticmodel (Mi et al., 2020); (b) Pseudo-two-dimensional VS Section of the inversion results, the dashed line in (b) represents the dipping interface.

4 實際資料應(yīng)用

將本文方法應(yīng)用于同震偏移(2008年汶川MW7.9地震)勘探的瑞雷波數(shù)據(jù)中.同震偏移的確定對于理解斷裂帶力學(xué)機制具有重要的意義(Scholz, 2019; Wang et al., 2021).2008年汶川MW7.9地震發(fā)生在青藏高原東緣龍門山逆沖帶；該斷裂帶主要由三條北東向逆沖斷層組成，分別為茂縣—汶川斷裂、江油—灌縣斷裂和映秀—北川斷裂.

作為發(fā)震斷裂之一的江油—灌縣斷裂同震垂向位移量為1～2 m(Xu et al., 2009; Zhang et al., 2010).本文選擇位于該斷裂的白鹿鎮(zhèn)同震地表破裂帶開展了瑞雷波勘探，旨在探測淺表層第四紀堆積物的速度結(jié)構(gòu)，為同震偏移的“局部吸收”效應(yīng)提供科學(xué)佐證(Wang et al., 2021).瑞雷波頻散曲線DL反演結(jié)果如圖12所示(數(shù)據(jù)采集采用滾動模式，最小偏移距1 m，道間距2 m，共24道，炮間距2 m，震源為錘擊震源，圖12a中水平距離10 m處的單炮地震記錄及對應(yīng)頻散曲線見附錄中圖A1).白鹿鎮(zhèn)同震地表破裂帶較好保存在原白鹿中學(xué)一帶，該部位主要發(fā)育有兩級階地，階地高差約4 m(圖13，圖14a).該處的同震破裂帶走向約北東向，表現(xiàn)為逆沖型破裂，其中在T2階地白鹿中學(xué)操場產(chǎn)生約1.8 m的垂向位錯量(BLL2，圖12a，圖14a、b).而在操場北東方向的T1階地白鹿河左岸一帶，同震垂直位錯量達到約2.5 m (BLL1，圖12b，圖14a、d)，Ran等(2010)在該觀測點以北測量一小路被同震垂向位錯約2.4 m(圖14a、c).因此，T1階地的同震垂向位錯量比T2階地大39%.考慮到操場和白鹿河岸旁觀測點之間的距離為185m(圖13b)，獲得兩個場地之間的同震垂向位錯量變化率為3.8 m·km-1.此外，橫跨T2階地操場地表破裂的瑞雷波勘探剖面顯示，第四系沉積物在斷層上盤的視厚度為9.6 m (BLS2，圖12a)，這與白鹿中學(xué)操場觀測點處開挖的大型探槽(圖14a、e)揭示的第四系地層厚度具有一致性，表明瑞雷波的深度學(xué)習(xí)反演方法具有較好的可靠性.同時本文反演T1階地觀測點處在斷層上盤的第四系沉積物視厚度為6.4m(BLS1，圖12b)，上述成果表明第四系松散地層可能在同震地表位錯變形的局部化方面具有重要的作用，能夠一定程度上吸收同震位錯變形量(Wang et al., 2021).

圖12 白鹿鎮(zhèn)地表破裂的瑞雷波反演結(jié)果(a) T2垂直偏移約1.8 m，斷層上盤顯示第四系沉積厚度約9.6 m； (b) T1垂直偏移約2.5 m，斷層上盤顯示第四系沉積，厚度約6.4 m.Fig.12 Rayleigh wave inversion results of surface ruptures at Bailu(a) The second terrace was vertically offset by about 1.8 m and shows Quaternary deposits with an apparent thickness of about 9.6 m at the hanging wall of the fault; (b) The first terrace was vertically offset by about 2.5 m and shows Quaternary deposits with an apparent thickness of about 6.4 m at the hanging wall of the fault.

圖13 白鹿鎮(zhèn)地形(a) 高分辨率圖像(谷歌地球)； (b) 地形解釋，紅線表示2008年MW7.9汶川地震有關(guān)的地表破裂，藍線和綠線分別表示瑞雷波剖面和地形調(diào)查剖面，“T”為階地.Fig.13 Landform at Bailu(a) A high-resolution image (Google Earth); (b)Landform interpretation. Red lines denote surface ruptures associated with the 2008 MW7.9 Wenchuan earthquake. Blue and green short lines show Rayleigh wave profiles and landform survey profiles, respectively. “T” marks terrace.

圖14 白鹿鎮(zhèn)地表破裂帶(a) 無人機拍攝的白鹿鎮(zhèn)照片； (b) T2上兩棟校舍之間的操場垂直偏移約1.8 m； (c) T1上的小路垂直偏移約2.4 m； (d) 白鹿河左岸T1面垂直偏移約2.5 m； (e) 在操場的地表破裂處挖掘了一條溝渠(Ran et al., 2019)研究表明次級階地的地層單元主要為砂質(zhì)黏土，含粒徑為0.2～5 cm的小礫石.黃色虛線表示階地和坡狀沉積的邊界，藍線和綠線分別表示地震剖面的地形調(diào)查剖面，紅線表示地表破裂與2008年有關(guān)MW7.9汶川地震，白色虛線表示與地震相關(guān)的突變斷層；圖片(b)、(c)、(e)分別攝于2010年、2008年和2009年.Fig.14 Surface fracture zone at Bailu(a) Photo ofBailu taken by drone; (b) The playground between two school buildings on the second terrace (T2) was vertically offset by about 1.8 m; (c) The small road on the first terrace (T1) was vertically offset by about 2.4 m; (d) Surface of T1 at the left bank of the Bailu River was vertically offset by about 2.5 m; (e) A trench was excavated across the surface rupture at the playground (Ran et al., 2019). The trench shows that the stratigraphic units of the second terrace mainly consist of sandy clay containing small gravels with diameters of 0.2～5 cm. Yellow dashed lines represent boundaries of terrace and sloperelated deposit. Blue and green short lines demote seismic profiles and landform survey profiles, respectively. “T” marks terrace. Red lines denote surface ruptures associated with the 2008 MW7.9 Wenchuan earthquake, and white dashed lines demote the coseismic fault scarps related to the earthquake. Photos (b)、(c)、(e) were taken in 2010, 2008, and 2009, respectively.

圖A1 (a)單炮地震記錄; (b)圖(a)對應(yīng)頻散圖 Fig.A1 (a) Single shot seismic record; (b) Dispersion curve corresponding to (a)

5 結(jié)論與討論

本文提出一種基于瑞雷波頻散曲線的DL反演新方法.DL主要包含了兩大重要的內(nèi)容，一是樣本數(shù)據(jù)，二是網(wǎng)絡(luò)結(jié)構(gòu).在樣本數(shù)據(jù)集構(gòu)建方面：本文針對地球物理反演的近地表速度結(jié)構(gòu)展開討論，闡述了近地表速度結(jié)構(gòu)的演化軌道也將力圖遍歷其相空間各種可能的體制和狀態(tài)，具有其各態(tài)遍歷性，以及其在演化中還有力圖保持自身特征形態(tài)的有序性特征.基于速度結(jié)構(gòu)的遍歷性和有序性，本文提出了約束馬爾科夫決策的近地表速度結(jié)構(gòu)模型，隨后對速度模型進行正演獲得了具有代表性的樣本數(shù)據(jù)對.隨機抽樣試驗表明了約束MC的樣本數(shù)據(jù)集構(gòu)建方法具有一定的合理性.在網(wǎng)絡(luò)結(jié)構(gòu)搭建方面：本文提出了一種CNN-LSTM的混合網(wǎng)絡(luò)結(jié)構(gòu)，該網(wǎng)絡(luò)結(jié)構(gòu)包含了3個局部特征學(xué)習(xí)模塊和1個長短時記憶層，局部特征學(xué)習(xí)模塊主要包含了雙卷積層和最大池化層，用于局部特征提取和層次關(guān)聯(lián)學(xué)習(xí)，LSTM從學(xué)習(xí)到的局部特征中繼續(xù)學(xué)習(xí)特征的長期依賴關(guān)系.這種混合網(wǎng)絡(luò)結(jié)構(gòu)能夠結(jié)合CNN和LSTM各自優(yōu)勢，取長補短，最終對頻散曲線的序列數(shù)據(jù)進行反演.

模型試驗表明本文方法具有較高的計算精度和良好的魯棒性，且計算效率極高，通常不到1 ms.以及將本文方法應(yīng)用于同震偏移勘探的瑞雷波數(shù)據(jù)反演中，其結(jié)果與實地地震地質(zhì)條件具有很好的一致性，表明本文反演方法具有較好的可信度.

本文DL的反演方法，在抗噪能力上有一定的特點，并與訓(xùn)練輪數(shù)相關(guān).訓(xùn)練輪數(shù)少，相應(yīng)的反演精度相對低，但魯棒性就越強，同時也對應(yīng)網(wǎng)絡(luò)結(jié)構(gòu)的泛化性能強；訓(xùn)練輪數(shù)多，相應(yīng)的反演精度就高，但魯棒性會變?nèi)?，即DL泛化性能會降低.因此，一方面如何平衡訓(xùn)練輪數(shù)與魯棒性之間的關(guān)系是一個值得深思的問題；另一方面，如何進一步同時提高計算精度和網(wǎng)絡(luò)的泛化性能，需要更深更合理的網(wǎng)絡(luò)結(jié)構(gòu)以及更加完備的樣本數(shù)據(jù)集.然而，這勢必會引起網(wǎng)絡(luò)復(fù)雜度與計算量增加.因此，在隨后的工作中，課題組將進一步對樣本數(shù)據(jù)構(gòu)建以及網(wǎng)絡(luò)結(jié)構(gòu)精修進行系統(tǒng)性的研究，以此更加符合真實地質(zhì)情況的應(yīng)用需求，從而使得頻散曲線的深度學(xué)習(xí)反演方法更加普適化.

致謝感謝匿名審稿專家對論文提出的寶貴修改建議，感謝浙江大學(xué)宓彬彬副研究員提供的二維模型及對應(yīng)的頻散速度數(shù)據(jù)！

附錄