基于自動(dòng)微分的圖空間最優(yōu)輸運(yùn)全波形反演

唐杰， 孟濤， 劉英昌， 孫成禹

中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 山東青島 266580

0 引言

全波形反演(Full Waveform Inversion，F(xiàn)WI)是一種通過最優(yōu)化算法減少合成地震數(shù)據(jù)與觀測地震數(shù)據(jù)之間的波形信息差異來重建地下介質(zhì)參數(shù)的高分辨率地震勘探方法(Tarantola，1984；劉玉柱等，2020).FWI首先通過正演模擬獲得合成地震記錄，然后利用目標(biāo)函數(shù)來衡量合成地震記錄與觀測地震記錄之間的差異，再通過優(yōu)化算法更新模型參數(shù).FWI具有精確刻畫模型細(xì)節(jié)的潛力，能夠?yàn)閮?chǔ)層解釋等提供可靠依據(jù).FWI除了在時(shí)間域?qū)崿F(xiàn)以外，還可以在頻率域以及Laplace域進(jìn)行，Pratt(1999)研究了頻率域全波形反演方法，提高了反演計(jì)算效率和穩(wěn)定性.Shin和Cha(2008)提出Laplace域全波形反演方法，隨后又研究了Laplace-Fourier域全波形反演(Shin and Cha，2009).

傳統(tǒng)的FWI以L2范數(shù)作為目標(biāo)函數(shù)，利用伴隨狀態(tài)法(Fichtner et al.，2006)計(jì)算梯度，沿著梯度下降的方向迭代更新地下介質(zhì)的模型參數(shù).以L2范數(shù)為目標(biāo)函數(shù)的FWI方法(L2 Full Waveform Inversion，L2-FWI)是一個(gè)高度非線性優(yōu)化問題，且反演過程中為了減小計(jì)算量通常使用局部優(yōu)化算法，當(dāng)初始模型較差時(shí)(或缺少長偏移距透射數(shù)據(jù)時(shí))，反演容易遇到周期跳變(cycle skipping)問題，無法獲得準(zhǔn)確的反演結(jié)果(Hu et al.，2018).此外，實(shí)際地震資料中往往缺乏低頻信息，并含有一定強(qiáng)度的噪聲，也會(huì)使得目標(biāo)函數(shù)在反演過程中陷入局部極小值，影響反演結(jié)果的準(zhǔn)確度.為了解決這一問題，眾多學(xué)者不斷尋找具有更好的凸性和噪聲魯棒性的目標(biāo)函數(shù)，近年來提出的一些目標(biāo)函數(shù)有：互相關(guān)目標(biāo)函數(shù)(Van Leeuwen and Mulder，2010)、包絡(luò)目標(biāo)函數(shù)(Wu et al.，2014；Luo and Wu，2015；胡勇等，2017)、卷積目標(biāo)函數(shù)(Zhang et al.，2016)和匹配濾波目標(biāo)函數(shù)(Zhu and Fomel，2016；Sun and Alkhalifah，2019)等.基于最優(yōu)輸運(yùn)理論的目標(biāo)函數(shù)(Engquist and Froese，2014；Métivier et al.，2016)因其對(duì)于數(shù)據(jù)的時(shí)移變化具有非常好的凸性，得到了廣泛的關(guān)注.最優(yōu)輸運(yùn)問題最早在1781年由Monge提出，1942年俄國數(shù)學(xué)家Kantorovich(1958)給出了一種更實(shí)際的松弛形式，定義了最優(yōu)輸運(yùn)距離，近年來最優(yōu)輸運(yùn)距離已成為數(shù)學(xué)界的一個(gè)研究熱點(diǎn)，被廣泛應(yīng)用于最優(yōu)化、計(jì)算機(jī)視覺和深度學(xué)習(xí)等眾多領(lǐng)域.Engquist和Froese(2014)首次將Wasserstain最優(yōu)輸運(yùn)距離作為目標(biāo)函數(shù)引入波形反演之中，并進(jìn)行了數(shù)學(xué)推導(dǎo)證明(Engquist et al.，2016).Yang等(2018)將二次Wasserstein度量應(yīng)用于FWI反演，在一定程度上緩解了反演過程中的周期跳變問題.Górszczyk等(2021)針對(duì)寬角海底地震儀數(shù)據(jù)研究了一種圖空間最優(yōu)輸運(yùn)(Graph Space Optimal Transport，GSOT)技術(shù)，降低了對(duì)反演初始模型準(zhǔn)確度的要求，取得了較好的實(shí)際應(yīng)用效果.Pladys等(2021)對(duì)最優(yōu)輸運(yùn)距離作為反演目標(biāo)函數(shù)的特性進(jìn)行了對(duì)比評(píng)估，展示了圖空間最優(yōu)輸運(yùn)目標(biāo)函數(shù)在反演中的應(yīng)用潛力.Cuturi(2013)為了克服最優(yōu)輸運(yùn)距離計(jì)算復(fù)雜度較高的缺點(diǎn)，使用熵正則化對(duì)最優(yōu)輸運(yùn)距離進(jìn)行熵約束，將最優(yōu)輸運(yùn)問題轉(zhuǎn)化為凸優(yōu)化問題，得到了Sinkhorn距離(Sinkhorn Distance，SD)，在保證最優(yōu)輸運(yùn)距離一定精度的前提下，節(jié)約了計(jì)算成本，加快了計(jì)算速度.

隨著地震勘探儀器的發(fā)展以及勘探技術(shù)的進(jìn)步，地震勘探已經(jīng)進(jìn)入大數(shù)據(jù)時(shí)代，數(shù)據(jù)科學(xué)和人工智能得到了日益廣泛的應(yīng)用.近年來GPU技術(shù)和深度學(xué)習(xí)技術(shù)與地震勘探領(lǐng)域的相關(guān)工作的結(jié)合日益緊密(桂生等，2017)，Biswas等(2019)使用物理信息引導(dǎo)的卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行疊前和疊后地震反演，較為準(zhǔn)確地預(yù)測了地下介質(zhì)的彈性特征.Zhang和Alkhalifah(2019)研究了一種相約束的彈性全波形反演方法，利用深度學(xué)習(xí)網(wǎng)絡(luò)來估計(jì)相分布，緩解了因初始模型不準(zhǔn)確引起的周期跳變問題.此外，在地震數(shù)據(jù)斷層檢測(Wu et al.，2019)、速度建模(Araya-Polo et al.，2019)、初至拾取(Yuan et al.，2020)、散射體成像(奚先和黃江清，2020)和最小二乘逆時(shí)偏移(Vamaraju et al.，2021)等工作中深度學(xué)習(xí)算法也得到了廣泛的應(yīng)用.自動(dòng)微分(Baydin et al.，2017；Li et al.，2020)是一種基于鏈?zhǔn)揭?guī)則的自動(dòng)梯度求取方法，在深度學(xué)習(xí)中主要用于訓(xùn)練由一系列線性變換和非線性激活函數(shù)組成的神經(jīng)網(wǎng)絡(luò)模型.基于偏微分方程的數(shù)值模擬(Cao and Liao，2015)類似于神經(jīng)網(wǎng)絡(luò)，在正演模擬中，利用模型參數(shù)、震源位置以及震源時(shí)間函數(shù)等得到模擬地震數(shù)據(jù)，然后在反演過程中，通過自動(dòng)微分計(jì)算模型參數(shù)的梯度，并利用深度學(xué)習(xí)中的優(yōu)化算法來反演更新模型參數(shù).基于深度學(xué)習(xí)軟件的自動(dòng)微分法能夠自動(dòng)求取反演過程中梯度，與常規(guī)的伴隨狀態(tài)法求得的梯度具有等價(jià)性(Sambridge et al.，2007；Tan et al.，2010)，而且便于在GPU上加速運(yùn)算.

本文提出了一種基于自動(dòng)微分和圖空間最優(yōu)輸運(yùn)的標(biāo)量波全波形反演方法，針對(duì)L2范數(shù)目標(biāo)函數(shù)存在的不足，將基于最優(yōu)輸運(yùn)理論的Sinkhorn距離作為反演的目標(biāo)函數(shù)，在反演過程中通過自動(dòng)微分自動(dòng)求取梯度，并利用深度學(xué)習(xí)中的Adam優(yōu)化算法更新模型參數(shù).基于Pytorch深度學(xué)習(xí)框架實(shí)現(xiàn)了以圖空間Sinkhorn距離作為目標(biāo)函數(shù)的全波形反演方法(Sinkhorn Distance Full Waveform Inversion，SD-FWI)，采用GPU加速，結(jié)合小批次(Minibatch)提高了反演速度.模型數(shù)據(jù)測試結(jié)果表明，SD-FWI能夠有效地緩解反演過程中的周期跳變問題，并且對(duì)初始子波、初始模型的依賴性較小.

1 原理

1.1 自動(dòng)微分全波形反演

FWI首先通過正演模擬獲得模型的理論計(jì)算波場，在檢波器位置記錄數(shù)據(jù)，得到合成地震記錄，然后通過迭代更新模型減少合成地震記錄和觀測地震記錄的差異來對(duì)地下介質(zhì)模型參數(shù)進(jìn)行反演.考慮二維常密度標(biāo)量波方程控制的全波形反演，具有恒定密度的標(biāo)量波方程滿足(Yang et al., 2018)：

(1)

其中，v(x,z)為波速，u(x,z,t)和S(x,z,t)分別為t時(shí)刻的波場值和震源值.使用有限差分法計(jì)算波場值，簡記為：

ut=v2Δt2(D2ut-1+ft-1)+2ut-1-ut-2,(2)

其中，ut為t時(shí)刻的波場值，D2為Laplace算子，Δt和ft-1分別為時(shí)間步長和t-1時(shí)刻的震源函數(shù)值.

在獲得標(biāo)量波合成記錄后，給定目標(biāo)函數(shù)，通過優(yōu)化算法降低目標(biāo)函數(shù)值，迭代更新模型參數(shù).對(duì)于觀測數(shù)據(jù)dobs，常用的L2范數(shù)目標(biāo)函數(shù)滿足：

(3)

其中，Nt為時(shí)間步長總數(shù)，R為波場限定算子.

FWI需要計(jì)算目標(biāo)函數(shù)相對(duì)于模型參數(shù)的梯度，使得模型參數(shù)能夠以降低目標(biāo)函數(shù)值的方向進(jìn)行更新.常規(guī)FWI方法一般利用伴隨狀態(tài)法(Plessix，2006)來計(jì)算反演過程中的梯度.自動(dòng)微分是一種基于鏈?zhǔn)揭?guī)則的梯度求取方法，首先通過正向模式計(jì)算每一個(gè)節(jié)點(diǎn)的梯度，再通過反向模式使用鏈?zhǔn)椒▌t計(jì)算目標(biāo)函數(shù)對(duì)于特定變量的梯度.自動(dòng)微分能夠利用計(jì)算機(jī)程序自動(dòng)計(jì)算反演過程中的梯度，求得的梯度與常規(guī)的伴隨狀態(tài)法求取的梯度具有等價(jià)性(Tan et al.，2010)，并且便于利用GPU并行運(yùn)算提升反演效率.

如圖1所示是基于自動(dòng)微分的全波形反演流程圖.首先給定反演初始模型，正演模擬得到合成地震記錄，然后通過目標(biāo)函數(shù)衡量真實(shí)地震記錄與合成地震記錄之間的差異.由于數(shù)值模擬波傳播的所有操作都是可微的，在反演過程中可以將自動(dòng)微分法與FWI結(jié)合，利用自動(dòng)微分計(jì)算機(jī)程序自動(dòng)計(jì)算模型參數(shù)的梯度.再不斷迭代直至達(dá)到預(yù)設(shè)的終止條件，最終保存反演結(jié)果.

圖1 反演流程圖Fig.1 Flow chart of inversion

1.2 最優(yōu)輸運(yùn)理論

最優(yōu)輸運(yùn)是一種能夠?qū)⒁粋€(gè)質(zhì)量分布轉(zhuǎn)化為另一個(gè)質(zhì)量分布的有效方法.設(shè)f和g分別為定義在X?d(d為X的維度)上的兩個(gè)概率測度，將兩者看作兩種不同的質(zhì)量分布，最優(yōu)輸運(yùn)問題就是尋找出一種能夠?qū)的質(zhì)量完全輸運(yùn)到g上耗費(fèi)“能量”最少的輸運(yùn)計(jì)劃.為了便于表示，輸運(yùn)計(jì)劃可以用P表示，P∈d×d.f映射到g所有可能的輸運(yùn)計(jì)劃構(gòu)成的集合為：

U(f,g)={P∈d×d|P1d=f,PT1d=g,P≥0}

(4)

f輸運(yùn)到g耗費(fèi)的成本矩陣可以用C表示，Cij=‖fi-gj‖2，‖.‖2為定義在d上的歐幾里得范數(shù).對(duì)于最優(yōu)的輸運(yùn)計(jì)劃，f和g之間的輸運(yùn)距離被稱為Wasserstein距離：

(5)

其中，Pij表示從fi向gj輸運(yùn)的總質(zhì)量.

經(jīng)典最優(yōu)輸運(yùn)距離的缺點(diǎn)是精確求解所需的計(jì)算復(fù)雜度較高(Cuturi，2013)，可以對(duì)Wasserstein距離添加熵約束，將最優(yōu)輸運(yùn)問題轉(zhuǎn)變?yōu)橐粋€(gè)凸優(yōu)化問題，降低計(jì)算成本，具有熵正則化約束的最優(yōu)輸運(yùn)距離也稱為Sinkhorn距離：

(6)

其中，λ是正則化系數(shù).

Sinkhorn距離可以通過Sinkhorn-Knopp定點(diǎn)迭代算法(Knight，2008)進(jìn)行快速求解，相比常規(guī)的Wasserstein距離具有更低的計(jì)算成本(Cuturi，2013).

1.3 圖空間Sinkhorn目標(biāo)函數(shù)

采用最優(yōu)輸運(yùn)距離作為度量函數(shù)時(shí)，被度量的數(shù)據(jù)需要滿足元素非負(fù)和質(zhì)量守恒兩個(gè)要求.由于地震信號(hào)通常是零值附近振蕩的，需要對(duì)地震數(shù)據(jù)進(jìn)行特定的預(yù)處理.目前常用的預(yù)處理方法可分為兩大類：第一類方法是通過規(guī)范化函數(shù)(線性、指數(shù)或Softplus函數(shù))將地震數(shù)據(jù)轉(zhuǎn)化為正值(Engquist and Froese，2014；Engquist et al.，2016；Yang et al.，2018)，但是這類方法缺點(diǎn)是會(huì)破壞最優(yōu)輸運(yùn)目標(biāo)函數(shù)的凸性或產(chǎn)生不可微點(diǎn)；第二類是一種基于1-Wasserstein距離的方法，雖然可以將最優(yōu)輸運(yùn)距離擴(kuò)展到度量非正數(shù)據(jù)上(Métivier et al.，2016)，但是會(huì)增加計(jì)算成本.本文采用的預(yù)處理方法是將地震數(shù)據(jù)轉(zhuǎn)化為圖空間中的點(diǎn)云(Métivier et al.，2019)，該方法可以有效保留地震信號(hào)的幾何形狀和最優(yōu)輸運(yùn)距離目標(biāo)函數(shù)對(duì)于信號(hào)時(shí)移的凸性，并且不會(huì)過多的增加計(jì)算成本.

首先，引入一個(gè)向量t∈K,t=(t1,…,tK)，地震道d(t)的離散圖是一個(gè)由2內(nèi)的K個(gè)點(diǎn)組成的集合，即(t1,d1),…,(tK,dK)，簡記為(t,d).地震數(shù)據(jù)(時(shí)間信號(hào))在時(shí)間離散后可變?yōu)閳D空間中的點(diǎn)云，每個(gè)點(diǎn)包含時(shí)間維度和振幅維度.對(duì)于離散的觀測數(shù)據(jù)(t,dobs)和合成數(shù)據(jù)(t,dcal)，在計(jì)算目標(biāo)函數(shù)時(shí)比較的是兩者的離散圖，而不是數(shù)據(jù)本身，這不僅保留了地震信號(hào)的幾何形狀而且可以確保最優(yōu)輸運(yùn)距離目標(biāo)函數(shù)對(duì)于信號(hào)時(shí)移變化的凸性.

使用圖空間Sinkhorn距離作為全波形反演的目標(biāo)函數(shù)：

(7)

將圖空間最優(yōu)輸運(yùn)Sinkhorn目標(biāo)函數(shù)與深度學(xué)習(xí)中的自動(dòng)微分相結(jié)合，可以方便地利用GPU并行加速運(yùn)算(Peyré and Cuturi，2019)，提高全波形反演的計(jì)算效率.

2 目標(biāo)函數(shù)凸性分析

2.1 Ricker子波測試

利用具有時(shí)移、振幅不同、波形相反和主頻不同的Ricker子波信號(hào)，進(jìn)行四種情況的一維數(shù)值模擬實(shí)驗(yàn)，將Ricker子波f固定為參考信號(hào)，相應(yīng)的時(shí)移變化信號(hào)g的中心沿時(shí)間軸0.3 s移動(dòng)到0.7 s，并以L2范數(shù)作為對(duì)比，對(duì)Sinkhorn距離作為目標(biāo)函數(shù)的凸性進(jìn)行研究分析.

如圖2所示，圖2a是Ricker子波信號(hào)f與其時(shí)移信號(hào)g的波形圖.圖2b是Ricker子波信號(hào)f和其縮放信號(hào)g=0.5f，模擬了波在地下傳播時(shí)發(fā)生的振幅衰減現(xiàn)象.圖2c表示是波形相反的兩個(gè)Ricker子波信號(hào)，用來模擬波在地下傳播時(shí)可能會(huì)遇到不同阻抗差異的地層使得波的極性發(fā)生反轉(zhuǎn)的現(xiàn)象.此外，波在地層的傳播過程中，還會(huì)發(fā)生頻率衰減改變信號(hào)頻譜成分的情況，采用了如圖2d所示的兩個(gè)具有不同主頻的Ricker子波對(duì)這種情況進(jìn)行模擬.以上四種不同情況的歸一化后的L2范數(shù)距離和Sinkhorn距離變化曲線如圖3所示.

圖2 四種不同情況的Ricker子波(a) 時(shí)移； (b) 振幅不同； (c) 波形相反； (d) 主頻不同.Fig.2 Ricker wavelets in four different cases(a) Time shift； (b) Different amplitude； (c) Opposite waveform； (d) Different dominant frequency.

圖3 四種不同情況下的歸一化L2范數(shù)距離和Sinkhorn距離(a) L2范數(shù)距離； (b) Sinkhorn距離.Fig.3 Normalized L2 norm distance and normalized Sinkhorn distance under four different conditions(a) L2 norm distance； (b) Sinkhorn distance.

如圖3a所示，四種情況下的L2范數(shù)距離變化曲線中均具有較多的局部極值，這些局部極值對(duì)于反演將產(chǎn)生非常不利的影響.相比于L2范數(shù)，Sinkhorn距離作為目標(biāo)函數(shù)在四種情況下均可以得到一個(gè)全局極小值(圖3b)，即使信號(hào)之間發(fā)生了時(shí)移、振幅縮放、極性翻轉(zhuǎn)和主頻變化等處理仍然具有更加良好的凸性.

2.2 抗噪性分析

當(dāng)采用實(shí)際數(shù)據(jù)進(jìn)行FWI時(shí)，觀測數(shù)據(jù)中往往會(huì)存在噪聲，因此目標(biāo)函數(shù)對(duì)于噪聲的魯棒性也是衡量目標(biāo)函數(shù)好壞的一個(gè)重要標(biāo)準(zhǔn).如圖4所示，利用圖4a所示的含有噪聲的信號(hào)f_noise和不含噪聲的信號(hào)g進(jìn)行了抗噪性測試分析.圖4b和圖4c分別是L2范數(shù)距離測試結(jié)果和Sinkhorn距離測試結(jié)果，結(jié)合兩者在時(shí)移0.1～0.4 s的局部圖可以看到，在含有噪聲的情況下，L2范數(shù)距離曲線(圖4d)會(huì)受到噪聲影響產(chǎn)生抖動(dòng)，而Sinkhorn距離曲線(圖4e)仍然較為平滑，這初步說明Sinkhorn距離目標(biāo)函數(shù)對(duì)于噪聲具有一定的魯棒性.

圖4 含噪信號(hào)測試(a) 含噪信號(hào)與無噪信號(hào)圖; (b) L2范數(shù)距離; (c) Sinkhorn距離; (d) L2范數(shù)距離0.1～0.4 s的局部圖; (e) Sinkhorn距離0.1～0.4 s的局部圖.Fig.4 Noisy signal test(a) Signal with noise and signal without noise； (b) L2 norm distance； (c) Sinkhorn distance； (d) Local figure of L2 norm distance from 0.1 s to 0.4 s； (e) Local figure of Sinkhorn distance from 0.1 s to 0.4 s

圖5 Marmousi速度模型(a) 真實(shí)速度模型； (b) 初始速度模型.Fig.5 Marmousi velocity model(a) True velocity model； (b) Initial velocity model.

圖6 L2-FWI與SD-FWI的反演結(jié)果對(duì)比(a) L2-FWI結(jié)果； (b) SD-FWI結(jié)果.Fig.6 Comparison of inversion results between L2-FWI and SD-FWI(a) L2-FWI results； (b) SD-FWI results.

3 波形反演測試

3.1 Marmousi模型

3.1.1 目標(biāo)函數(shù)比較

Marmousi模型是地震勘探中用來測試各種正演反演算法的常用模型，其結(jié)構(gòu)較為復(fù)雜，具有大型的斷層、背斜構(gòu)造以及高速地層等地質(zhì)單元.如圖5a所示，使用截取并重采樣的部分Marmousi模型測試Sinkhorn距離作為目標(biāo)函數(shù)在復(fù)雜介質(zhì)模型反演中的適應(yīng)性和穩(wěn)定性.選用的Marmousi模型橫向距離為3 km，縱向距離為1 km，模型網(wǎng)格大小為5 m，對(duì)真實(shí)速度模型進(jìn)行高斯濾波平滑后得到與真實(shí)模型差異較大的初始速度模型(圖5b).

在模型頂部以100 m等間距均勻地放置了30個(gè)震源，并以10 m等間距均勻地放置300個(gè)檢波器.震源是主頻為10 Hz的雷克子波，采樣時(shí)間為2 s，采樣間隔為2 ms.反演測試共迭代400輪，并利用型號(hào)為Tesla T4(15G顯存)的GPU并行加速運(yùn)算，反演總耗時(shí)約為2800 s左右，反演結(jié)果如圖6所示.

圖6a和圖6b分別是L2-FWI的反演結(jié)果和SD-FWI的反演結(jié)果，從上到下分別為兩種方法迭代0、20、100和400次的反演結(jié)果.對(duì)比可知，在迭代200次時(shí)，SD-FWI相比于L2-FWI已經(jīng)得到了較好的反演結(jié)果，對(duì)于模型中的高速地層和斷層反演更加清晰準(zhǔn)確，尤其是對(duì)于模型的深部特征的反演更加準(zhǔn)確.

為了能夠更加清晰直觀地對(duì)比，分別從真實(shí)模型、初始模型、L2-FWI和SD-FWI的反演結(jié)果中提取了水平距離為1 km和2 km處的縱向速度曲線.如圖7所示，對(duì)于Marmousi模型500 m以上的淺部速度，L2-FWI和SD-FWI的反演結(jié)果相差不大，但是對(duì)于500 m以下的深部地層速度，SD-FWI的反演結(jié)果更加準(zhǔn)確.

圖7 L2-FWI反演結(jié)果與SD-FWI反演結(jié)果單道對(duì)比(a) 水平距離為1 km處的縱向速度曲線； (b) 水平距離為2 km處的縱向速度曲線.Fig.7 Single traces comparison between L2-FWI result and SD-FWI result(a) Vertical velocity curves at 1 km horizontal distance； (b) Vertical velocity curves at 2 km horizontal distance.

圖8 不同迭代次數(shù)的反演結(jié)果殘差的波數(shù)譜(a) L2-FWI結(jié)果； (b) SD-FWI結(jié)果.Fig.8 The wave-number spectrums of inversion results residual with different iterations(a) L2-FWI results； (b) SD-FWI results.

圖8展示了L2-FWI和SD-FWI不同迭代次數(shù)的反演結(jié)果與真實(shí)模型的殘差在波數(shù)域的變化趨勢.可以看出SD-FWI對(duì)模型進(jìn)行更新時(shí)具有全局特征，整個(gè)波數(shù)范圍內(nèi)的模型殘差相比于L2-FWI下降速度更快，并且對(duì)于低波數(shù)信息具有更好的敏感性，在重構(gòu)低波數(shù)分量時(shí)比L2-FWI更穩(wěn)健，能夠有效恢復(fù)更多的模型低波數(shù)結(jié)構(gòu)信息.

為了定量地衡量SD-FWI和L2-FWI在反演過程中的反演結(jié)果與真實(shí)模型之間的差異，本文定義了一個(gè)偏差函數(shù)：

(8)

其中，minv為反演得到的速度模型，mtrue為真實(shí)速度模型，Misfit表示反演結(jié)果與真實(shí)模型之間偏差值，偏差值越小說明反演結(jié)果與真實(shí)模型相似性越高.

圖9給出了反演過程中的偏差變化曲線，可以看到SD-FWI相比于L2-FWI，不僅反演過程中前期的偏差下降速率更快，而且后期能夠收斂到更低的值，反演結(jié)果與真實(shí)模型相似性更高.

圖9 反演結(jié)果的偏差變化曲線Fig.9 Misfit curve of inversion results

3.1.2 深度學(xué)習(xí)優(yōu)化算法的選擇

全波形反演本質(zhì)是一個(gè)非線性優(yōu)化問題，優(yōu)化算法會(huì)對(duì)反演過程中的計(jì)算成本、收斂速度和最終反演結(jié)果產(chǎn)生較大的影響.基于深度學(xué)習(xí)框架實(shí)現(xiàn)FWI，除了可以利用自動(dòng)微分對(duì)梯度進(jìn)行求取之外，還可以利用深度學(xué)習(xí)中先進(jìn)的優(yōu)化算法.

深度學(xué)習(xí)中常用的優(yōu)化算法包括Adagrad、RMSprop和Adam等.Adagrad算法會(huì)將所有參數(shù)的每一次迭代的梯度取平方累加后再開方，用全局學(xué)習(xí)率除以這個(gè)數(shù)，作為學(xué)習(xí)率的動(dòng)態(tài)更新量，在迭代前期可以激勵(lì)放大梯度，后期懲罰縮小梯度，能夠?yàn)椴煌淖兞刻峁┎煌膶W(xué)習(xí)率.RMSprop優(yōu)化算法會(huì)設(shè)定權(quán)重變化加速因子與減速因子，在迭代過程中當(dāng)連續(xù)誤差梯度符號(hào)不變時(shí)，采用加速策略，加快訓(xùn)練速度；當(dāng)連續(xù)誤差梯度符號(hào)變化時(shí)，采用減速策略，以期穩(wěn)定收斂.RMSprop優(yōu)化算法步長的更新緩和，在反演的前期可以加速模型參數(shù)更新，后期有利于防止更新幅度過大而跨過最優(yōu)值.Adam優(yōu)化算法(Kingma and Ba，2017)由隨機(jī)梯度下降法發(fā)展而來，是一種基于一階梯度的隨機(jī)目標(biāo)函數(shù)優(yōu)化算法，通過計(jì)算參數(shù)梯度的一階矩估計(jì)和二階矩估計(jì)對(duì)每個(gè)參數(shù)的學(xué)習(xí)速率進(jìn)行動(dòng)態(tài)調(diào)整，可以快速收斂到一個(gè)較好的模型，適合于數(shù)據(jù)或參數(shù)量較大的情況以及非平穩(wěn)目標(biāo)的問題，具有收斂速度快、所需內(nèi)存小等優(yōu)點(diǎn).

為了測試深度學(xué)習(xí)中的優(yōu)化算法在SD-FWI中的表現(xiàn)，對(duì)Adagrad、RMSprop和Adam的不同學(xué)習(xí)率(learning rate，Lr)進(jìn)行了測試和比較，三種優(yōu)化算法不同學(xué)習(xí)率的偏差函數(shù)變化曲線如圖10所示.

比較圖10a、圖10b和圖10c，可知Adagrad、RMSprop和Adam優(yōu)化算法在SD-FWI測試實(shí)驗(yàn)中較為合適的學(xué)習(xí)率分別為100、10和20左右，三種優(yōu)化算法在學(xué)習(xí)率選取合適的情況下，應(yīng)用于SD-FWI均可以取得一定的效果.Adagrad (Lr=100)、RMSprop (Lr=10)和Adam (Lr=20)的反演結(jié)果如圖10所示.當(dāng)然在實(shí)際反演中需要結(jié)合數(shù)據(jù)特征通過測試分析選取合理的學(xué)習(xí)率.

圖10 不同深度學(xué)習(xí)優(yōu)化算法的反演結(jié)果偏差變化曲線(a) Adagrad； (b) RMSprop； (c) Adam； (d) 三種優(yōu)化算法對(duì)比.Fig.10 Misfit curves of inversion results with different deep learning optimization algorithms(a) Adagrad； (b) RMSprop； (c) Adam； (d) Comparison of three optimization algorithms.

如圖11所示，三種優(yōu)化算法在分別選擇其較為合適的學(xué)習(xí)率的情況下，整體上均取得了不錯(cuò)的反演結(jié)果，模型中的地層和主要斷面信息均在一定程度上得到了有效地恢復(fù).這進(jìn)一步表明了深度學(xué)習(xí)中的這三種深度學(xué)習(xí)優(yōu)化算法應(yīng)用于SD-FWI是可行且有效的.

圖11 不同優(yōu)化算法的反演結(jié)果(a) Adagrad (Lr=100)； (b) RMSprop (Lr=10)； (c) Adam (Lr=20).Fig.11 Inversion results of different optimization algorithms

如圖12所示，分別為從真實(shí)模型、初始模型和三種優(yōu)化算法的反演結(jié)果中提取的水平距離為1 km和2 km處的縱向速度變化曲線.可以看到，三種優(yōu)化算法的反演結(jié)果均能較好的恢復(fù)模型的整體速度變化趨勢，相比于Adagrad優(yōu)化算法，Adam與RMSprop優(yōu)化算法對(duì)模型的淺層速度變化趨勢反演更加精細(xì)、抖動(dòng)更小.綜合對(duì)比三種優(yōu)化算法在反演測試過程中最優(yōu)學(xué)習(xí)率的偏差函數(shù)變化曲線(圖10d)可以發(fā)現(xiàn)，采用學(xué)習(xí)率為20的Adam優(yōu)化算法在實(shí)驗(yàn)中的偏差曲線整體下降較快，可以收斂到更低.因此，本文其他測試實(shí)驗(yàn)中的優(yōu)化算法均將選用Adam優(yōu)化算法.

圖12 不同優(yōu)化算法反演結(jié)果單道對(duì)比(a) 水平距離1 km處的縱向速度曲線； (b) 水平距離2 km處的縱向速度曲線.Fig.12 Single traces comparison of inversion results with different optimization algorithms(a) Vertical velocity curves at 1 km horizontal distance； (b) Vertical velocity curves at 2 km horizontal distance.

3.1.3 初始震源子波對(duì)反演結(jié)果的影響

初始震源子波對(duì)FWI的影響也是不可忽視的，F(xiàn)WI通常假設(shè)初始震源子波是已知的，但在實(shí)際地震數(shù)據(jù)反演過程中初始震源子波往往是不準(zhǔn)確的，錯(cuò)誤的初始震源子波可能會(huì)使得反演收斂到局部最小值，甚至得到一個(gè)錯(cuò)誤的反演結(jié)果(Luo et al.，2014).采用如圖13所示兩種不準(zhǔn)確的初始震源子波對(duì)SD-FWI進(jìn)行反演測試：第一種情況如圖13a所示，初始子波與真實(shí)子波能量不同；第二種情況是初始子波與真實(shí)子波的主頻和相位均不同，如圖13b所示，初始子波由一個(gè)與真實(shí)子波主頻不同的雷克子波經(jīng)過翻轉(zhuǎn)和時(shí)移得到，兩個(gè)子波在波形上沒有任何重疊的部分.反演測試模型與參數(shù)設(shè)置均與3.1.1節(jié)相同.反演過程中初始震源子波與速度模型同時(shí)更新，優(yōu)化過程中對(duì)模型參數(shù)和子波的學(xué)習(xí)率設(shè)置不同(由于速度模型與子波的尺度大小不同，本文為速度模型和子波設(shè)置的學(xué)習(xí)率分別為20和0.001).

子波反演結(jié)果分別如圖13a和圖13b中的黑色虛線所示，可以看到初始震源子波能量不準(zhǔn)確時(shí)，子波能量雖然沒有得到恢復(fù)，但是模型仍然可以取得較好的反演結(jié)果(圖14a)；在初始震源子波與真實(shí)子波差距很大的情況下，雖然子波的反演結(jié)果右側(cè)有一些抖動(dòng)，但是與真實(shí)子波對(duì)應(yīng)位置的波形卻得到了較好的恢復(fù).

圖13 不準(zhǔn)確的初始震源子波(a) 不同能量的震源子波； (b) 不同主頻和相位的震源子波.Fig.13 Inaccurate initial seismic wavelets(a) Wavelets with different energy； (b) Wavelets with different dominant frequency and phase.

圖14為使用上述的兩種不準(zhǔn)確震源子波的模型反演結(jié)果，從圖14a可以看出，震源子波能量不同對(duì)SD-FWI的影響較小，反演結(jié)果中斷層的斷面和高速層等顯示清晰.圖14b是使用不同主頻和相位的震源子波的反演結(jié)果，由于使用的震源子波與真實(shí)子波差異巨大，反演結(jié)果中出現(xiàn)了部分偽影干擾.

圖14 不準(zhǔn)確震源子波的反演結(jié)果(a) 不同能量的震源子波的反演結(jié)果； (b) 不同主頻和相位的震源子波的反演結(jié)果. Fig.14 Inversion results of inaccurate seismic wavelets(a) Inversion result of seismic wavelets with different energy； (b) Inversion result of seismic wavelets with different dominant frequency and phase.

3.1.4 噪聲對(duì)反演結(jié)果的影響

如圖15所示，為兩個(gè)含有不同強(qiáng)度高斯隨機(jī)噪聲的單道數(shù)據(jù)，圖15a中記錄的信噪比(SNR)為1 dB，表示強(qiáng)噪聲情況；圖15b中的噪聲較小，信噪比為10 dB.分別使用信噪比SNR=10 dB和SNR=1 dB的含噪聲數(shù)據(jù)重復(fù)前面的FWI實(shí)驗(yàn)，模型選取和參數(shù)設(shè)置與先前實(shí)驗(yàn)相同，Adam優(yōu)化算法的學(xué)習(xí)率為5.含噪數(shù)據(jù)反演測試結(jié)果如圖16所示.

圖15 單道含噪數(shù)據(jù)(a) SNR=1 dB； (b) SNR=10 dB.Fig.15 Single trace of noisy data

圖16 含噪數(shù)據(jù)反演結(jié)果(a) L2-FWI, SNR=10 dB； (b) SD-FWI, SNR=10 dB； (c) L2-FWI, SNR=1 dB； (d) SD-FWI, SNR=1 dB.Fig.16 Inversion results of noisy data

圖16a和圖16b分別是地震數(shù)據(jù)中含有信噪比為10 dB的噪聲情況下L2-FWI和SD-FWI得到的反演結(jié)果.可以看出，SD-FWI相比于L2-FWI方法取得了更好的反演結(jié)果，層位分層較好，斷層的斷面顯示清晰.當(dāng)?shù)卣饠?shù)據(jù)中噪聲強(qiáng)度較高時(shí)(SNR=1)，SD-FWI的反演結(jié)果(圖16d)相比于L2-FWI的反演結(jié)果(圖16c)仍能恢復(fù)Marmousi2模型中的大部分特征.這表明SD-FWI在反演時(shí)具有一定的抗噪性.

3.1.5 一維初始模型反演測試

為了測試本文方法對(duì)于初始模型的依賴性，使用了如圖17a所示的一維初始速度模型，進(jìn)行反演測試(其他反演相關(guān)參數(shù)設(shè)置與上文3.1.1節(jié)相同).

由圖17c可以看出，盡管初始模型與真實(shí)模型偏離較大，SD-FWI相比于L2-FWI(圖17b)仍然得到了較好的反演結(jié)果，地層速度反演較為準(zhǔn)確，斷層的斷面顯示清晰.這表明本文方法在反演過程中對(duì)于精確初始模型的依賴性較小.

圖17 一維初始速度模型測試結(jié)果(a) 一維初始速度模型； (b) L2-FWI結(jié)果； (c) SD-FWI結(jié)果.Fig.17 Test results of 1D initial velocity model(a) 1D initial velocity model； (b) L2-FWI result； (c) SD-FWI result.

3.2 BP2004模型

BP2004模型根據(jù)墨西哥灣西部的一個(gè)復(fù)雜地質(zhì)橫截面建立，具有強(qiáng)烈的非均勻性，由于高阻鹽體的存在會(huì)產(chǎn)生強(qiáng)反射，反演具有一定的難度.為了進(jìn)一步測試SD-FWI方法對(duì)于具有強(qiáng)反射的鹽體和高傾角高速體的反演效果，使用截取并重采樣后的BP2004模型(Billette and Brandsberg-Dahl，2005)進(jìn)行了反演測試.

如圖18a所示，模型的深度為1 km，橫向延伸3 km.反演初始模型選用如圖18b所示的一維初始速度模型.反演測試時(shí)，在模型頂部以100 m等間距均勻地放置30個(gè)震源，并在模型頂部以10 m間距均勻地放置300個(gè)檢波器.震源采用主頻為10 Hz的雷克子波，時(shí)間采樣間隔為2 ms，采樣時(shí)間為2 s.

圖18 BP2004模型反演結(jié)果(a) 真實(shí)速度模型； (b) 一維初始速度模型； (c) L2-FWI結(jié)果； (d) SD-FWI結(jié)果.Fig.18 Inversion results of BP2004 model(a) True velocity model； (b) 1D initial velocity model； (c) L2-FWI result； (d) SD-FWI result.

反演共迭代400輪，如圖18c為L2-FWI的反演結(jié)果，鹽體和高傾角高速體的形狀沒有得到有效恢復(fù)，這是L2范數(shù)作為目標(biāo)函數(shù)在反演過程中陷入局部極小值的典型情況.而在SD-FWI反演結(jié)果(圖18d)中，左側(cè)鹽體得到了較好反演，整體輪廓顯示清晰準(zhǔn)確，且右側(cè)的高傾角高速體的特征也得到了準(zhǔn)確恢復(fù).

如圖19所示，為真實(shí)模型、1D初始速度模型、L2-FWI反演結(jié)果和SD-FWI反演結(jié)果中水平距離0.5 km和2.2 km處的單道速度變化曲線.由圖可知，SD-FWI方法能夠很好地恢復(fù)出模型中鹽體和高傾角高速體的速度變化趨勢.

圖19 BP2004模型反演結(jié)果單道對(duì)比(a) 水平距離為0.5 km處的縱向速度曲線； (b) 水平距離為2.2 km處的縱向速度曲線.Fig.19 Single trace comparison of BP2004 model inversion results(a) Vertical velocity curves at 0.5 km horizontal distance； (b) Vertical velocity curves at 2.2 km horizontal distance.

由圖20可知，L2-FWI的反演結(jié)果與真實(shí)模型之間的偏差逐步增大，反演朝著錯(cuò)誤的方向進(jìn)行，而SD-FWI的反演偏差變化曲線能夠有效下降，得到了有效的迭代更新.

圖20 反演結(jié)果與真實(shí)BP2004速度模型之間的偏差曲線Fig.20 Misfit curves between inversion result and true BP2004 velocity model

4 結(jié)論

本文提出了一種基于自動(dòng)微分和圖空間最優(yōu)輸運(yùn)Sinkhorn距離的標(biāo)量波全波形反演方法，通過數(shù)值模擬分析和模型測試，驗(yàn)證了本文方法的有效性，得到了以下結(jié)論：

(1)圖空間Sinkhorn距離作為目標(biāo)函數(shù)對(duì)信號(hào)時(shí)移和振幅變化具有良好的凸性，能夠緩解反演過程中的周期跳變問題.利用走時(shí)不同、振幅不同、極性相反和主頻不同的Ricker子波信號(hào)進(jìn)行測試，結(jié)果表明圖空間Sinkhorn距離相比于L2范數(shù)作為目標(biāo)函數(shù)在四種情況下均可以得到一個(gè)全局極小值，具有更加良好的凸性.此外Sinkhorn距離相比于L2范數(shù)具有更好的噪聲魯棒性.

(2)本文采用自動(dòng)微分計(jì)算反演過程中的模型參數(shù)梯度，利用深度學(xué)習(xí)中的Adam優(yōu)化算法更新速度模型，并采用Minibatch方法減少內(nèi)存占用和利用GPU加速運(yùn)算.模型測試結(jié)果表明，本文研究方法具有一定的抗噪性，對(duì)于震源子波和初始模型的依賴性較低，可以有效改善FWI結(jié)果并且反演效率較高.

基于自動(dòng)微分計(jì)算梯度有助于研究人員在全波形反演方面嘗試新的想法和進(jìn)行快速實(shí)驗(yàn)測試.本文采用Minibatch技術(shù)能夠?qū)ε跀?shù)和檢波器數(shù)進(jìn)行分批次反演，在一定程度上減小了反演的內(nèi)存需求，當(dāng)然在實(shí)際三維地震數(shù)據(jù)反演應(yīng)用中仍然存在內(nèi)存消耗較大的問題，想要對(duì)三維實(shí)際數(shù)據(jù)進(jìn)行反演，仍需要在大型工作站或計(jì)算機(jī)集群上進(jìn)行，如何優(yōu)化反演過程中的內(nèi)存分配是一個(gè)非常值得探索的重要方向.此外，本文研究的是標(biāo)量波情況下的反演，將本文方法擴(kuò)展到彈性波情況下將會(huì)具有更多實(shí)際應(yīng)用的意義.