李松齡

(東北石油大學(xué) 地球科學(xué)學(xué)院，黑龍江 大慶 163318)

0 前 言

勘探地球物理主要的目的是獲得對(duì)地球內(nèi)部的精確成像。近年來(lái)，隨著計(jì)算硬件性能的提升以及GPU(Graphic Processing Unit)等高性能計(jì)算設(shè)備在地震數(shù)據(jù)成像中的成功應(yīng)用，具有高分辨能力的全波形反演(Full Waveform Inversion，F(xiàn)WI)方法受到了地球物理學(xué)家們青睞[1-2]。對(duì)于照明均勻和低頻可用的觀測(cè)數(shù)據(jù)，理論上，F(xiàn)WI能夠同時(shí)恢復(fù)地下構(gòu)造的長(zhǎng)波長(zhǎng)和短波長(zhǎng)信息，相當(dāng)于同時(shí)對(duì)地震數(shù)據(jù)做層析和偏移處理[3-4]。因此，具有高分辨率的全波形反演方法在油氣勘探開(kāi)發(fā)和地球動(dòng)力學(xué)研究等領(lǐng)域中彰顯了更大的研究潛力和發(fā)展前景。

時(shí)間域的FWI方法最先由Tarantola實(shí)現(xiàn)，隨后，Pratt將FWI發(fā)展到頻率域[5]，并發(fā)展了由低頻到高頻的反演策略，Bunk在時(shí)域發(fā)展了由低頻到高頻多尺度反演策略[6]。時(shí)域和頻率域的FWI均有各自的優(yōu)勢(shì)，全頻帶的頻率FWI等同于時(shí)域FWI。頻率FWI在二維情況下優(yōu)勢(shì)巨大，在頻率域僅執(zhí)行一次LU分解，便可實(shí)現(xiàn)多炮的正演，但對(duì)于三維情況下，頻率域正演求解變得非常困難。所以，目前多數(shù)三維FWI方法的實(shí)現(xiàn)常常在時(shí)間域完成。

除了初始模型的問(wèn)題會(huì)干擾FWI的應(yīng)用效果，巨大的計(jì)算量和存儲(chǔ)量同樣也成為影響FWI實(shí)用化推廣的棘手問(wèn)題。在FWI迭代計(jì)算過(guò)程中，反演梯度的計(jì)算與逆時(shí)偏移(Reverse Time Migration，RTM)方法流程類似，因此，采用互相關(guān)成像方法時(shí)，震源波場(chǎng)的存儲(chǔ)量決定著RTM和FWI方法對(duì)計(jì)算硬件存儲(chǔ)性能的要求。同理，針對(duì)RTM方法的存儲(chǔ)策略同樣可借鑒應(yīng)用至FWI方法中。Clapp提出在速度外圍構(gòu)建隨機(jī)速度層作為邊界條件[7]，將邊界波場(chǎng)能量以隨機(jī)噪聲的形式保存于波場(chǎng)內(nèi)部，進(jìn)而通過(guò)波動(dòng)方程的反向計(jì)算即可獲得歷史波場(chǎng)，該方法由于波場(chǎng)存儲(chǔ)量較低，被廣泛應(yīng)用于RTM方法中；Shen等對(duì)上述方法進(jìn)行了改進(jìn)，加強(qiáng)了隨機(jī)邊界對(duì)低頻數(shù)據(jù)的散射強(qiáng)度，并應(yīng)用于波形反演算法中[8]。但隨機(jī)邊界條件產(chǎn)生的隨機(jī)噪聲會(huì)對(duì)RTM及FWI結(jié)果產(chǎn)生干擾，影響最終的成像精度。Clapp提出了一種保存邊界波場(chǎng)的存儲(chǔ)策略[9]，能夠一定程度上的降低波場(chǎng)存儲(chǔ)量，而且可應(yīng)用于吸收邊界條件下的波場(chǎng)重建；隨后，該方法被眾多學(xué)者加以改進(jìn)并應(yīng)用至RTM計(jì)算和GPU加速策略中[10-11]。FENG和WANG等提出了差分階數(shù)的漸變策略[12]，用于進(jìn)行波場(chǎng)重建，但該方法重建的波場(chǎng)存有一定誤差，雖適用于RTM成像，但不適用于基于迭代計(jì)算的FWI方法。Kalita和Alkhalifah提出利用激發(fā)振幅成像方法實(shí)現(xiàn)FWI算法中成像梯度的計(jì)算[13]，能夠大幅降低波場(chǎng)數(shù)據(jù)的存儲(chǔ)量和計(jì)算量，但對(duì)于復(fù)雜構(gòu)造，激發(fā)振幅成像方法往往在多波至路徑的區(qū)域會(huì)產(chǎn)生走時(shí)計(jì)算誤差，影響成像精度。地球物理學(xué)家們目前仍在積極探索針對(duì)FWI方法中波場(chǎng)存儲(chǔ)量的改進(jìn)方法，該方面的研究探討對(duì)提升FWI的理論水平和應(yīng)用價(jià)值十分重要。

在前人研究工作的基礎(chǔ)上，本文提出了一種基于波場(chǎng)重建的時(shí)域FWI方法。借助拉格朗日乘子法，首先推導(dǎo)了FWI的梯度表達(dá)式，并引入解析步長(zhǎng)方法和混合共軛梯度優(yōu)化方法求解該反問(wèn)題。針對(duì)FWI成像梯度計(jì)算所需的巨大存儲(chǔ)量問(wèn)題，我們引入了一種能大幅降低存儲(chǔ)成本的波場(chǎng)重建策略，為時(shí)域FWI的實(shí)用化推廣提供一套切實(shí)可行的理論方案。最后通過(guò)勻速模型和Marmousi2模型的波場(chǎng)重建和FWI測(cè)試驗(yàn)證了本文方法的有效性。

1 方法原理

1.1 全波形反演的梯度

時(shí)空域二階聲波常密度波動(dòng)方程為：

(1)

式中，x是空間坐標(biāo)；xs為炮點(diǎn)空間坐標(biāo)；v是地震波傳播速度；p是地震壓力波場(chǎng)；fs是震源項(xiàng)。定義在檢波點(diǎn)處的模擬的數(shù)據(jù)滿足：

pcal(xr,t;xs)=f(v)

(2)

式中，xr是檢波點(diǎn)的空間坐標(biāo)；pcal為模擬的地震壓力波場(chǎng)；f為波動(dòng)方程正演算子，即方程(1)。

全波形反演(FWI)通過(guò)最小化模擬數(shù)據(jù)和野外觀測(cè)數(shù)據(jù)之間的誤差，推導(dǎo)地下的速度參數(shù)，其目標(biāo)函數(shù)通常定義為模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù)殘差的L2范數(shù)

(3)

式中，ns和ng分別為炮點(diǎn)和檢波點(diǎn)數(shù)目；pobs為野外觀測(cè)的地震數(shù)據(jù)；tmax表示記錄地震波的最大時(shí)間。

該目標(biāo)函數(shù)受到模擬數(shù)據(jù)滿足的波動(dòng)方程的限制，我們采用拉格朗日乘子法法求解該約束性優(yōu)化問(wèn)題，建立拉格朗日乘子符為：

(4)

式中，λ(x,t；xs)為未知的拉格朗日乘子函數(shù)。目標(biāo)函數(shù)的梯度表達(dá)式可以根據(jù)拉格朗日乘子符隨速度的偏導(dǎo)數(shù)得到

(5)

(6)

該伴隨方程與聲波方程(1)具有相同的表達(dá)形式，因此只需要將震源項(xiàng)替換為模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù)的殘差，即可利用相同的數(shù)值求解方法模擬震源波場(chǎng)和伴隨波場(chǎng)。

因此，傳統(tǒng)的聲波全波形反演的流程可以簡(jiǎn)述為：

1)沿時(shí)間正向求解方程(1)，得到正向傳播的震源波場(chǎng)pcal(x,t；xs)，并保存所有時(shí)刻的正傳波場(chǎng)在計(jì)算機(jī)內(nèi)存中。

2)沿時(shí)間反方向求解方程(6)，得到反向傳播的檢波點(diǎn)波場(chǎng)λ(x,t；xs)。

3)從計(jì)算機(jī)內(nèi)存中提取保存的波場(chǎng)，應(yīng)用公式(5)，構(gòu)建FWI的梯度。

4)采用合適的反演算法，計(jì)算優(yōu)化的步長(zhǎng)，更新速度參數(shù)v。

5)重復(fù)1～4步驟，直到目標(biāo)函數(shù)值小于設(shè)定的值或滿足最大迭代次數(shù)。

1.2 全波形反演的解析步長(zhǎng)

上一節(jié)，我們已經(jīng)推導(dǎo)了FWI的梯度，給出了傳統(tǒng)FWI的實(shí)現(xiàn)流程。然而，為了得到合理的反演結(jié)果，最優(yōu)化的搜索步長(zhǎng)也尤其重要，常規(guī)的步長(zhǎng)求取方法有很多，如：固定步長(zhǎng)[14]、拋物線擬合[15]等方法。在每1次的迭代過(guò)程中，拋物線擬合的方法對(duì)每炮數(shù)據(jù)至少需要進(jìn)行兩次正演計(jì)算，才能夠給出合理的步長(zhǎng)，計(jì)算量過(guò)大。因此，本文采用解析步長(zhǎng)的方法[16-17]，該方法僅需對(duì)每炮數(shù)據(jù)額外進(jìn)行一次正演正演計(jì)算即可得到優(yōu)化步長(zhǎng)。

將第k+1次迭代的目標(biāo)函數(shù)進(jìn)行二階泰勒展開(kāi)，可得：

(7)

(8)

(9)

式中，f為式(1)所示的正演算子，ε為試探步長(zhǎng)，它滿足條件

(10)

1.3 混合共軛梯度方法

常規(guī)最速下降方法在迭代尋優(yōu)的過(guò)程，收斂速度非常緩慢，會(huì)帶來(lái)鋸齒折疊效應(yīng)。為了在迭代過(guò)程中，求解一個(gè)合理的下降方向dk，本文我們將一種混合類的共軛梯度方法應(yīng)用到FWI方法中，以加速反演的收斂速度。第k次迭代的下降方向dk可以根據(jù)下列公式得到：

(11)

其中，βk滿足[18-19]：

(12)

在得到共軛梯度方向和優(yōu)化步長(zhǎng)后，速度參數(shù)則能夠迭代地更新，更新公式可以表示為：

vk+1=vk+αkdk

(14)

2 基于波場(chǎng)重建的時(shí)域FWI

在傳統(tǒng)FWI的實(shí)現(xiàn)過(guò)程中，需要對(duì)背景震源波場(chǎng)和伴隨波場(chǎng)進(jìn)行零延遲互相關(guān)來(lái)構(gòu)建梯度。但二者在時(shí)間方向上互逆，因此可預(yù)先保存某一方向所有時(shí)刻的波場(chǎng)值，當(dāng)另一個(gè)波場(chǎng)傳播時(shí)，提取相應(yīng)時(shí)刻保存的波場(chǎng)進(jìn)行互相關(guān)運(yùn)算，但由此產(chǎn)生巨大的存儲(chǔ)需求和I/O(Input/Output)耗時(shí)通常是難以接受的。因此，本文提出一種適用于FWI的逆時(shí)重建震源波場(chǎng)的方法，旨在降低了FWI方法的存儲(chǔ)需求和I/O耗時(shí)。

本文采用高階有限差分法求解波動(dòng)方程(1)：

波動(dòng)方程差分公式

p(x,z,t+Δt)=2p(x,z,t)-p(x,z,t-Δt)+

(15)

如圖1所示，假設(shè)真實(shí)波場(chǎng)范圍為圖1中的A區(qū)域，如果采用空間12階差分計(jì)算格式求解波動(dòng)方程，僅需保存A區(qū)域外側(cè)6層波場(chǎng)數(shù)據(jù)，如圖1中B區(qū)域，即可滿足差分條件，C區(qū)域?yàn)槲者吔鐓^(qū)域。

圖1 邊界存儲(chǔ)示意圖

為了實(shí)現(xiàn)震源波場(chǎng)的逆時(shí)重建，震源波場(chǎng)沿正向傳播時(shí)，我們需要記錄每時(shí)刻的B區(qū)域的波場(chǎng)信息以及最后時(shí)刻的A區(qū)域的波場(chǎng)信息。利用這些保存的信息，我們可以實(shí)現(xiàn)逆時(shí)波場(chǎng)重建：首先通過(guò)波動(dòng)方程的逆向運(yùn)算實(shí)現(xiàn)震源波場(chǎng)的逆時(shí)傳播，同時(shí)，將對(duì)應(yīng)時(shí)刻的B區(qū)域內(nèi)的波場(chǎng)使用正傳過(guò)程中保存的波場(chǎng)更新，以滿足差分運(yùn)算條件。在震源波場(chǎng)逆時(shí)重建過(guò)程中，結(jié)合公式5即可進(jìn)行梯度的構(gòu)建，直到所有時(shí)刻和所有炮完成疊加。其余反演步驟與傳統(tǒng)方法一致，綜上所述，本文提出的基于波場(chǎng)重建的時(shí)域FWI流程如圖2所示。

圖2 基于波場(chǎng)重建的時(shí)域FWI流程

文中首先采用簡(jiǎn)單的均勻各向同性的介質(zhì)進(jìn)行震源波場(chǎng)逆時(shí)重建方法測(cè)試，該模型的水平和垂直方向網(wǎng)格點(diǎn)均為301，沿水平和垂直方向的網(wǎng)格尺寸均為10 m。該模型的縱波傳播速度為2 800 m/s，為了吸收來(lái)自由于邊界截?cái)喈a(chǎn)生的反射波，我們應(yīng)用了完美匹配層吸收邊界(PML)[20]，吸收邊界厚度為50個(gè)網(wǎng)格點(diǎn)。正演時(shí)，震源函數(shù)為主頻20 Hz的雷克子波，該震源位于模型(1.5 km, 1.5 km)，時(shí)間步長(zhǎng)為1 ms。

圖3 均勻各向同性介質(zhì)的波場(chǎng)快照(t=0.15，0.30，0.45，0.60 s)

圖3為均勻各向同性介質(zhì)不同時(shí)刻的波場(chǎng)快照，其中圖3a和圖3b分別顯示了正向傳播和本文方法逆時(shí)重建的震源波場(chǎng)，從左自右的時(shí)刻分別為0.15，0.30，0.45，0.60 s。我們不難發(fā)現(xiàn)，模型區(qū)域內(nèi)的震源波場(chǎng)已經(jīng)被很好地重建出來(lái)，由于PML邊界內(nèi)的波場(chǎng)不參與梯度的計(jì)算，因此并沒(méi)有重建。我們進(jìn)一步對(duì)波場(chǎng)快照抽取單道進(jìn)行對(duì)比，正傳波場(chǎng)和重建波場(chǎng)的單道數(shù)據(jù)對(duì)比如圖4所示，抽取位置在t=0.45 s，x=2 km處。觀察可知，正傳的單道圖4a和逆時(shí)重建的單道圖4b幾乎沒(méi)有差別。兩者的誤差單道圖4c所示，可知單道誤差10-4的數(shù)量級(jí)。

圖4 正傳及逆時(shí)重建的波場(chǎng)(t=0.45s)單道數(shù)據(jù)對(duì)比

為了證明文章中所述方法能夠節(jié)約存儲(chǔ)，我們對(duì)比了傳統(tǒng)全波場(chǎng)存儲(chǔ)的方法和文章所提出方法。假設(shè)，該模型的傳播總時(shí)間間隔為5 000，波場(chǎng)數(shù)據(jù)在GPU端的數(shù)據(jù)類型為4字節(jié)浮點(diǎn)型，全波場(chǎng)存儲(chǔ)的計(jì)算方法為：nx·hz·nt·4byte，本文的方法的存儲(chǔ)為z·(nx+nz)·nt·6·4byte，其中6表示為有限差分方法的空間差分階數(shù)的1/2，nx和nz為模型沿水平和垂直方向的空間網(wǎng)格點(diǎn)數(shù)。存儲(chǔ)消耗對(duì)比如表1所示，這驗(yàn)證了本文所述方法能夠大幅降低存儲(chǔ)。特別地，當(dāng)模型為3D時(shí)，文中所述的方法將大幅降低存儲(chǔ)消耗，加速反演的計(jì)算速度。

表1 勻速模型波場(chǎng)重建存儲(chǔ)消耗對(duì)比

3 模型測(cè)試

為了驗(yàn)證文中所述波場(chǎng)重建方法和全波形反演方法的有效性，分別對(duì)Marmousi2模型進(jìn)行的波場(chǎng)重建測(cè)試和FWI方法測(cè)試。

3.1 波場(chǎng)重建方案測(cè)試

首先，進(jìn)行了基于Marmousi2模型的波場(chǎng)重建測(cè)試，速度模型如圖5所示。模型沿水平和垂直方向的網(wǎng)格點(diǎn)數(shù)分別為500和200，沿水平和垂直方向的網(wǎng)格間距均為10 m，為了壓制了來(lái)自邊界的虛假反射，我們?cè)谀Ｐ退闹芗虞d了厚度為50層PML吸收邊界。震源位置位于水平方向2.5 km，垂直方向0.65 km處，震源函數(shù)為主頻20 Hz的雷克子波，時(shí)間采樣間隔為1 ms。

圖5 Marmousi2 速度模型

分別截取正傳至t=0.20，0.40，0.60，0.80 s時(shí)刻的震源波場(chǎng)，如圖6a所示，當(dāng)?shù)卣鸩▊鞑ブ磷畲髸r(shí)刻5 s時(shí)，進(jìn)行波場(chǎng)逆時(shí)重建，重建所得的t=0.20，0.40，0.60，0.80 s時(shí)刻的波場(chǎng)如圖6b所示。根據(jù)圖6對(duì)比可知，本文所述的逆時(shí)重建方法能夠很好地恢復(fù)出在模型區(qū)域內(nèi)的震源波場(chǎng)。進(jìn)一步對(duì)波場(chǎng)快照抽取單道數(shù)據(jù)進(jìn)行對(duì)比，所選數(shù)據(jù)為t=0.40 s，x=3 km處的波場(chǎng)數(shù)據(jù)，對(duì)比結(jié)果如圖7所示。正傳波場(chǎng)和逆時(shí)重建波場(chǎng)的的單道數(shù)據(jù)分別如圖7a和7b所示，二者并無(wú)肉眼可見(jiàn)差別，對(duì)兩者做相減運(yùn)算，分析所得誤差如圖7c所示，誤差為10-3數(shù)量級(jí)，與實(shí)際有效波振幅相差4個(gè)數(shù)量級(jí)，對(duì)真實(shí)結(jié)果無(wú)明顯影響。本次測(cè)試中，基于Marmousi2模型的波場(chǎng)重建方案和傳統(tǒng)全波場(chǎng)保存方案的存儲(chǔ)消耗對(duì)比如表2所示，相比于全波場(chǎng)保存方案，本文方法能降低約90%的存儲(chǔ)量，效果顯著。

圖6 Marmousi模型的波場(chǎng)快照(t=0.20，0.40，060，0.80 s)

圖7 正傳及逆時(shí)逆時(shí)反傳的波場(chǎng)(t=0.45 s)單道對(duì)比

模型名稱模型大小/m總時(shí)間間隔/s全波場(chǎng)存儲(chǔ)/M邊界存儲(chǔ)/MMarmousi2500×20050001907160

3.2 Marmousi2全波形反演測(cè)試

為了驗(yàn)證本文提出的FWI方法對(duì)復(fù)雜模型的有效性，對(duì)Marmousi2模型進(jìn)行了FWI測(cè)試，精確速度模型如圖5所示。模型的基本信息與波場(chǎng)重建測(cè)試相同。震源和檢波器均置于地面，震源的范圍在0～5 km，炮間距為250 m，共20炮，采用主頻為10 Hz的雷克子波，共500個(gè)檢波器，檢波點(diǎn)間距為10 m。時(shí)間采樣間隔為1 ms，總時(shí)間長(zhǎng)度為5 s。

由于實(shí)際地震勘探無(wú)法獲得精確的速度模型，因此在進(jìn)行FWI測(cè)試時(shí)，首先將圖5所示的精確速度進(jìn)行平滑，如圖8所示，以此作為初始速度模型。不同迭代次數(shù)的反演結(jié)果如圖9所示。隨著迭代次數(shù)的增加，反演結(jié)果逐漸接近真實(shí)速度，對(duì)比不同迭代次數(shù)的反演結(jié)果可進(jìn)一步發(fā)現(xiàn)，初期的反演結(jié)果對(duì)中淺部的速度構(gòu)造更新效果較好，當(dāng)?shù)螖?shù)超過(guò)50次后，深層的速度構(gòu)造也得以更新，趨近真實(shí)速度。

分別提取了水平方向1，2，3，4 km處的反演結(jié)果的單道數(shù)據(jù)，如圖10所示，黑色實(shí)線表示初始速度，藍(lán)色實(shí)線表示真實(shí)速度，綠色和紅色實(shí)線分別表示迭代20次和100次的反演結(jié)果。曲線對(duì)比的結(jié)果與前文所述結(jié)論一致，執(zhí)行20次迭代計(jì)算后，文中所述的方法可以很好地反演出中淺層的速度，隨著迭代次數(shù)的增加，深層構(gòu)造也得到了有效更新。圖11顯示了目標(biāo)函數(shù)隨迭代次數(shù)的下降情況，隨著迭代計(jì)算，理論模擬數(shù)據(jù)與觀測(cè)數(shù)據(jù)誤差逐漸減少，進(jìn)一步驗(yàn)證了本文方法的有效性。

圖8 Marmousi 平滑速度模型

圖9 不同迭代次數(shù)Marmousi的反演結(jié)果

圖10 不同迭代次數(shù)的單道對(duì)比圖

圖11 目標(biāo)函數(shù)隨迭代次數(shù)的下降曲線

4 結(jié) 論

波場(chǎng)數(shù)據(jù)的存儲(chǔ)策略對(duì)全波形反演方法的應(yīng)用要求至關(guān)重要，本文從有限差分的原理出發(fā)，提出了保存局部波場(chǎng)數(shù)據(jù)的波場(chǎng)重建策略，能夠大幅降低全波形反演的波場(chǎng)存儲(chǔ)需求。多種模型試算也表明本文方法的正確性和有效性。

本文方法對(duì)基于復(fù)雜方程(如一階速度應(yīng)力方程和彈性波方程等)的FWI方法或三維FWI均具有一定的借鑒意義，這也是本文下一步的研究重點(diǎn)；本文提出的波場(chǎng)重建策略，僅適用于非損耗介質(zhì)下的波場(chǎng)重建，對(duì)于粘性介質(zhì)的波場(chǎng)重建方法仍有待進(jìn)一步研究。

參考文獻(xiàn)：

[1] Tarantola A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 1984, 49(8): 1259-1266.

[2] Pratt R G. Seismic waveform inversion in the frequency domain Part 1 theory and verification in a physical scale model[J]. Geophysics, 1999, 64(3): 888-901.

[3] Mora P. Inversion = migration + tomography[J]. Geophysics, 1989, 54(12): 1575-1586.

[4] Wang H, Singh S C, Audebert F, et al. Inversion of seismic refraction and reflection data for building long-wavelength velocity models[J]. Geophysics, 2015, 80(2): R81-R93.

[5] Pratt R G, Shin C, Hick G J, et al. Gauss Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophysical Journal International, 1998, 133(2): 341-362.

[6] Bunks C, Saleck F M, Zaleski S, et al. Multiscale seismic waveform inversion[J]. Geophysics, 1995; 60(5): 1457-1473.

[7] Clapp R G. Reverse time migration with random boundarie[J]. Seg Technical Program Expanded Abstracts, 2009, 2809-2813.

[8] Shen X, Clapp R G. Random boundary condition for memory-efficient waveform inversion gradient computation[J]. Geophysics, 2015, 80(6): R351-R359.

[9] Clapp R G. Reverse time migration: Saving the boundaries[J]. Stanford Exploration Project, 2008,136-144.

[10] 王保利, 高靜懷, 陳文超,等. 地震疊前逆時(shí)偏移的有效邊界存儲(chǔ)策略[J]. 地球物理學(xué)報(bào), 2012, 55(7): 412-2421.

[11] Yang P, Gao J, Wang B, et al. RTM using effective boundary saving: A staggered grid GPU implementation[J]. Computers & Geosciences, 2014; 64-72.

[12] FENG Bo, WANG Huazhong. Reverse time migration with source wavefield reconstruction strategy[J]. Journal of Geophysics and Engineering, 2012, 9(1): 69-74.

[13] Kalita M, Alkhalifah T. Efficient full waveform inversion using the excitation representation of the source wavefield[J]. Geophysical Journal International, 2017, 210(3): 581-1594.

[14] Kohn D, De Nil D, Kurzmann A, et al. On the influence of model parametrization in elastic full waveform tomography[J]. Geophysical Journal International, 2012, 191(1): 325-345.

[15] Vigh D, Starr E W. 3D prestack plane-wave, full-waveform inversion[J]. Geophysics, 2008, 73(5): VE135-VE144.

[16] Gauthier O, Virieux J, Tarantola A, et al. Two-dimensional nonlinear inversion of seismic waveforms: Numerical results[J]. Geophysics, 1986, 51(7): 1387-1403.

[17] Pica A, Diet J P, Tarantola A, et al. Nonlinear inversion of seismic reflection data in a laterally invariant medium[J]. Geophysics, 1990; 55(3): 284-292.

[18] Hager W W, Zhang H. A survey of nonlinear conjugate gradient methods[J]. Pacific Journal of Optimization, 2006, 2(1): 35-58.

[19] Yang P, Gao J, Wang B, et al. A graphics processing unit implementation of time-domain full-waveform inversion[J]. Geophysics, 2015, 80(3): F31-F39.

[20] 王守東. 聲波方程完全匹配層吸收邊界[J]. 石油地球物理勘探, 2003, 38(1):31-34.