張生強(qiáng)，劉春成，韓立國(guó)，楊小椿

1.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

2.中國(guó)海洋石油研究總院，北京 100027

0 引言

在地震勘探的每一個(gè)環(huán)節(jié)中，速度均起著至關(guān)重要的作用。目前沒(méi)有任何直接手段可以獲得地下介質(zhì)的速度分布，通常只能利用地表觀測(cè)到的地震資料來(lái)確定地層速度。用于確定地層速度的基本方法有速度分析、射線層析成像和疊前（逆時(shí)）偏移等，但得到的結(jié)果均比較粗糙。全波形反演可以利用疊前地震記錄，迭代反演出高精度的地下介質(zhì)速度。它不僅能夠顯示地下介質(zhì)的構(gòu)造信息，還能顯示地下介質(zhì)的速度變化。因此，研究全波形反演方法具有極其廣泛的意義和價(jià)值。

全波形反演技術(shù)經(jīng)歷了漫長(zhǎng)的發(fā)展過(guò)程。20世紀(jì)80年代，Tarantola［1－2］提出了基于廣義最小二乘反演理論的時(shí)間域全波形反演方法，對(duì)近20多年多維地震反演理論的發(fā)展產(chǎn)生了深遠(yuǎn)的影響。Marfurt［3］明確指出，對(duì)于含有大量炮點(diǎn)的波動(dòng)方程模擬問(wèn)題，使用頻率域有限元或有限差分法是最有效的手段。為提高計(jì)算效率，20世紀(jì)80年代末90年代初，Pratt等［4－7］將這一思想與 Tarantola的全波形反演方法相結(jié)合，發(fā)展了頻率域聲波和彈性波全波形反演方法。頻率域全波形反演中各個(gè)頻率數(shù)據(jù)是相互獨(dú)立的，可以并行計(jì)算，相對(duì)于時(shí)域全波形反演具有計(jì)算高效、數(shù)據(jù)選擇靈活等優(yōu)勢(shì)。頻率采樣間隔選擇具有一定的靈活性，早期的頻率域反演通常使用等間隔頻率采樣。Song等［8］首先將低頻數(shù)據(jù)波形反演結(jié)果作為相鄰高頻數(shù)據(jù)波形反演的初始模型以提高穩(wěn)定性。該方法隨后被廣泛應(yīng)用于全波形反演研究中，并被稱為串行反演策略，如Sourbier等［9］應(yīng)用該方法，有效地防止了求解過(guò)程陷入目標(biāo)函數(shù)的局部極小值。Pratt［10］隨后提出將有效頻段內(nèi)頻率域數(shù)據(jù)由低頻到高頻分組，每組內(nèi)各頻率成分同時(shí)參與反演（即組內(nèi)并行反演策略），將較低頻分組反演結(jié)果作為相鄰較高頻率分組數(shù)據(jù)反演的初始模型（即組間串行反演策略），有利于壓制噪音和反演假象。

全波形反演計(jì)算成本高、效率低、內(nèi)存需求大，這些都是限制其發(fā)展的主要原因。筆者提出在頻率多尺度全波形反演中將L－BFGS（limited memory BFGS）算法與同時(shí)激發(fā)震源技術(shù)相結(jié)合的方法來(lái)改善這一現(xiàn)狀。同時(shí)激發(fā)震源技術(shù)可以大大地提高波形反演過(guò)程中正演模擬的計(jì)算效率。求解全波形反演問(wèn)題的優(yōu)化算法有很多，如最速下降法、共軛梯度法和擬牛頓法等，Nocedal［11］提出的 L－BFGS數(shù)值優(yōu)化方法是求解優(yōu)化問(wèn)題最有效的擬牛頓方法之一。2007年Qin Zhiwei（Tony）［12］對(duì)共軛梯度法、BFGS算法以及2種有限內(nèi)存優(yōu)化算法進(jìn)行了對(duì)比，證明了L－BFGS算法是一種占用內(nèi)存小、收斂速度快、計(jì)算精度高的算法，適用于求解大規(guī)模地球物理優(yōu)化問(wèn)題。但是由于全波形反演是一個(gè)強(qiáng)非線性反問(wèn)題，針對(duì)這個(gè)問(wèn)題L－BFGS算法可能會(huì)收斂到局部極值，得不到全局最優(yōu)解。為了克服這個(gè)缺點(diǎn)，筆者使用頻率多尺度反演方法，即將有效頻段內(nèi)頻率域數(shù)據(jù)由低頻到高頻分組，采取組內(nèi)并行反演、組間串行的反演策略。這樣可以有效地避免收斂于局部極值，通常都能夠?qū)ふ业饺肿顑?yōu)值。最后通過(guò)數(shù)值模擬試驗(yàn)對(duì)這種頻率多尺度全波形反演方法及其抗噪能力和反演速度擾動(dòng)能力進(jìn)行了測(cè)試，收到了很好的成效。

1 方法原理

頻率域波動(dòng)方程全波形反演是利用波形，即波場(chǎng)的振幅、頻率以及相位信息，將數(shù)據(jù)殘差反傳波場(chǎng)與初始模型的入射波場(chǎng)互相關(guān)，更新地下介質(zhì)參數(shù)，進(jìn)而對(duì)地下介質(zhì)物性參數(shù)成像的方法。其具有揭示復(fù)雜地質(zhì)背景下構(gòu)造與巖性細(xì)節(jié)信息的潛力?；诓▌?dòng)方程的全波形反演本身是個(gè)強(qiáng)非線性反問(wèn)題，解這類問(wèn)題的經(jīng)典方法是首先將問(wèn)題線性化，再通過(guò)逐步迭代最終確定一個(gè)使理論計(jì)算值與觀測(cè)記錄較好擬合的反演結(jié)果。從數(shù)學(xué)的角度來(lái)講，全波形反演可以對(duì)應(yīng)地分為2個(gè)問(wèn)題：大型稀疏矩陣線性方程組的求解以及非線性局部最優(yōu)化問(wèn)題。

1.1 頻率域正演

頻率域正演是頻率域全波形反演的基礎(chǔ)，正演的精度和效率在很大程度上決定了反演的精度和效率。本文研究了二階常密度聲波方程頻域正演。目前，頻率域全波形反演中使用的正演方法包括有限差分法、有限元法和有限體積法三大類，筆者采用有限差分法（9點(diǎn)有限差分格式）求解頻域聲波方程，即將時(shí)空域的微分方程轉(zhuǎn)化為頻率－空間域的大型稀疏矩陣線性方程組問(wèn)題，并采用LU分解法求解。該方法的優(yōu)點(diǎn)在于可以重復(fù)使用已分解的矩陣因子，提高了求解相同頻率不同震源的正演計(jì)算效率。

1.1.1 Helmholtz方程

對(duì)時(shí)域二階標(biāo)量聲波方程作關(guān)于時(shí)間的Fourier變換可得到頻率域標(biāo)量聲波方程，即Helmholtz方程［13］：

其中：▽2為 Laplace算子；A＝－▽2－k2為Helmholtz算子；u為地震波場(chǎng)；k＝ω/v為波數(shù)，ω為角頻率，v為介質(zhì)速度；S為震源項(xiàng)。

采用LU分解法求解方程（1）即可求得地震波場(chǎng)u。地震波場(chǎng)u和模型參數(shù)m之間的關(guān)系是非線性的，可以通過(guò)算子D（D為m和u之間的非線性函數(shù)）來(lái)定義：

1.1.2 同時(shí)激發(fā)震源

全波形反演中，正演模擬最佳震源選擇的一個(gè)重要標(biāo)準(zhǔn)是多炮模擬的有效性。2008年Berkhout［14］提出了混合激發(fā)采集，即不同空間位置的多個(gè)震源按照隨機(jī)線性編碼方式激發(fā)組成超級(jí)炮。筆者在全波形反演過(guò)程中選擇了混合激發(fā)采集策略（同時(shí)激發(fā)震源技術(shù)，如圖1），主要原因有二：1）提高聲波有限差分算法多炮正演模擬的計(jì)算效率，使用同時(shí)激發(fā)震源技術(shù)來(lái)代替每一炮獨(dú)立的正演，這樣可以在一個(gè)炮記錄里同時(shí)模擬多炮正演；2）可以用更多的炮來(lái)對(duì)地下進(jìn)行照明。

圖1 混合激發(fā)采集（同時(shí)激發(fā)震源）方法示意圖Fig.1 Schematic diagram of blended acquisition（simultaneous sources）method

超炮道集的構(gòu)建可由下式表達(dá)：

式中，P指包含多個(gè)單炮記錄共頻道集Pj的道集集合?；旌喜杉捎秒S機(jī)線性編碼，不同震源按照一定的隨機(jī)時(shí)間延遲激發(fā)。而對(duì)于同時(shí)激發(fā)震源，也采用隨機(jī)線性編碼，只是不同震源按照相同的時(shí)間激發(fā)。

當(dāng)使用單炮正演時(shí)，每炮正演耗時(shí)為tseq，共Ns炮，則正演總用時(shí)為Nstseq。若對(duì)這些震源采用隨機(jī)線性編碼，按照相同的時(shí)間激發(fā)，組成同時(shí)激發(fā)震源，此時(shí)正演耗時(shí)為tsim。在這2種情況下，由于tsim?Nstseq，所以同時(shí)激發(fā)震源技術(shù)大大地提高了正演模擬的計(jì)算效率［15］。

1.2 全波形反演

對(duì)于地震勘探中的每一個(gè)震源－檢波器，定義波場(chǎng)殘差Δd為檢波器位置處的觀測(cè)波場(chǎng)dobs與正演模擬計(jì)算的波場(chǎng)dcal（m）之間的差值，即

式中：模型m表示地下離散計(jì)算域中的一些物理參數(shù)，本文中這些模型參數(shù)是指定義在每個(gè)數(shù)值網(wǎng)格節(jié)點(diǎn)上的P波速度，下文將用v來(lái)代替m表示速度模型。這時(shí)需要定義一個(gè)目標(biāo)函數(shù)C（v）用來(lái)衡量波場(chǎng)觀測(cè)值與計(jì)算值之差。全波形反演的目的就是通過(guò)迭代更新使C（v）達(dá)到最小，最終求出v。根據(jù)最小二乘原理，目標(biāo)函數(shù)為

其中，＊表示共軛轉(zhuǎn)置算子。

1.2.1 波恩近似及反演問(wèn)題的線性化

根據(jù)波恩近似理論，更新后的速度模型v可以表示為速度模型初值v0與速度模型擾動(dòng)Δv之和，即

在初始速度模型v0點(diǎn)附近對(duì)目標(biāo)函數(shù)進(jìn)行二階Taylor－Lagrange展開(kāi)：

其中：整數(shù)M表示速度模型v中的元素個(gè)數(shù)；O（v3）表示誤差項(xiàng)。式（7）兩端同時(shí)對(duì)速度模型參數(shù)vl求導(dǎo)并寫(xiě)成矩陣形式：

當(dāng)目標(biāo)函數(shù)的一階導(dǎo)數(shù)為0時(shí)，則目標(biāo)函數(shù)在v0附近達(dá)到極小值，同時(shí)得到速度模型擾動(dòng)的表達(dá)式：

其中：目標(biāo)函數(shù)的一階導(dǎo)矩陣為梯度，表示在v0點(diǎn)處目標(biāo)函數(shù)最速上升的方向；目標(biāo)函數(shù)的二階導(dǎo)矩陣為Hessian矩陣，表示目標(biāo)函數(shù)在v0點(diǎn)附近的彎曲度。求解方程（9）的方法通常稱為牛頓方法。地震波場(chǎng)和速度模型之間的關(guān)系是非線性的（u＝D（v）），在全波形反演中，反演需要迭代多次才能收斂到目標(biāo)函數(shù)的極小值。

1.2.2 L－BFGS算法

考慮如下優(yōu)化問(wèn)題：

其中：函數(shù)C∈R，是實(shí)值函數(shù)，稱C（v）為問(wèn)題（10）的目標(biāo)函數(shù)；集合Ω為問(wèn)題的定義域。

求解最優(yōu)化問(wèn)題（10），通常采用迭代算法，其迭代格式如下：

其中：步長(zhǎng)α（k）由線搜索策略確定；d（k）為搜索方向。在眾多的優(yōu)化方法中，Nocedal［11］提出的BFGS方法因其穩(wěn)定的數(shù)值效果和快速收斂性而被公認(rèn)為是求解優(yōu)化問(wèn)題（10）最有效的擬牛頓算法。該方法在構(gòu)造搜索方向時(shí)只需要利用目標(biāo)函數(shù)及其一階導(dǎo)數(shù)的信息，避免了Hessian矩陣的計(jì)算，減少了計(jì)算量，并且保持超線性收斂的優(yōu)點(diǎn)。它的基本思想是用某個(gè)給定的對(duì)稱正定矩陣H（k）近似于Hessian矩陣的逆，即通過(guò)校正公式逐步修正矩陣H（k），使其越來(lái)越近似于［▽2C（v（k））］－1。該算法中Hessian矩陣的逆近似H（k）的修正公式［11］為

但BFGS方法的一個(gè)顯著缺點(diǎn)是需要在內(nèi)存中存儲(chǔ)Hessian矩陣的逆近似，當(dāng)需要求解的問(wèn)題規(guī)模很大時(shí)，需要消耗很多的計(jì)算機(jī)內(nèi)存，因此不適合求解大規(guī)模優(yōu)化問(wèn)題（例如地球物理問(wèn)題）。針對(duì)該問(wèn)題，Nocedal［11］又提出了一種計(jì)算高效的 L－BFGS方法，其主要思想是：在當(dāng)前迭代點(diǎn)v（k）處，對(duì)給定的初始對(duì)稱正定矩陣和非負(fù)整數(shù)m，利用前面m步的信息，對(duì)進(jìn)行修正m次產(chǎn)生H（k）。與標(biāo)準(zhǔn)的BFGS方法不同，L－BFGS算法無(wú)需存儲(chǔ)矩陣H（k），僅存儲(chǔ)m＋1個(gè)向量組就能計(jì)算H（k＋1）。

利用式（12）導(dǎo)出L－BFGS算法第k＋1次迭代的修正矩陣H（k＋1）的表達(dá)式：

綜上所述，求解優(yōu)化問(wèn)題（10）的L－BFGS算法步驟如下。

第1步，給定初始速度模型v0∈Ω、初始對(duì)稱正定矩陣H0、非負(fù)整數(shù)m、誤差限ε＞0，令k＝0。

第3步，使用反向傳播線搜索方法，確定步長(zhǎng)α（k）＞0滿足下面的 Wolfe－Powell條件：

第4步，令m＝min｛k＋1，m｝，取由公式（13）確定

第5步，若k＝k＋1，轉(zhuǎn)第2步。

在 L－BFGS算法步驟中，計(jì)算H（k）g（k）的一個(gè)高效公式見(jiàn)文獻(xiàn)［11］。Liu和 Nocedal［16］證明，在適當(dāng)?shù)臈l件下，L－BFGS方法在求解一致凸函數(shù)極小化問(wèn)題時(shí)具有全局收斂性，大量的數(shù)值試驗(yàn)結(jié)果表明，此算法非常適合求解大規(guī)模優(yōu)化問(wèn)題。在實(shí)際計(jì)算中，可以根據(jù)問(wèn)題規(guī)模的大小以及計(jì)算機(jī)容許的內(nèi)存，通過(guò)選擇適當(dāng)?shù)膍來(lái)控制所需的存儲(chǔ)量。通常情況下，m的取值為3～20［16］。

1.2.3 頻率多尺度串行反演方法

頻率域全波形反演開(kāi)始前需要確定選用數(shù)據(jù)的最低和最高頻率。低頻數(shù)據(jù)對(duì)于恢復(fù)大尺度宏觀速度模型極為重要，因此反演最低頻率通常都設(shè)置在數(shù)據(jù)有效頻段內(nèi)能使用的最低頻率附近。高頻數(shù)據(jù)則對(duì)恢復(fù)速度模型細(xì)節(jié)信息極為重要，最大頻率通常由反演目標(biāo)、有效數(shù)據(jù)頻段、頻散誤差優(yōu)化和硬件條件等因素綜合確定。筆者使用的頻率多尺度串行反演方法是將有效頻段內(nèi)頻率域數(shù)據(jù)由低頻到高頻分組，每組內(nèi)各頻率成分同時(shí)參與反演，將較低頻分組反演結(jié)果作為相鄰較高頻率分組數(shù)據(jù)反演的初始模型，有利于壓制噪音和反演假象。

筆者在詳細(xì)介紹了同時(shí)激發(fā)震源方法、L－BFGS算法、頻率多尺度串行反演方法以及它們優(yōu)點(diǎn)的基礎(chǔ)上，提出了在頻率多尺度全波形反演中將LBFGS數(shù)值優(yōu)化算法與同時(shí)激發(fā)震源技術(shù)相結(jié)合的方法，以期能在計(jì)算效率和存儲(chǔ)量方面有所改進(jìn)。其實(shí)現(xiàn)流程是：首先需要在選定反演頻率的初始頻帶下給定一個(gè)初始速度模型，使用同時(shí)激發(fā)震源作正演模擬，利用觀測(cè)數(shù)據(jù)和正演模擬數(shù)據(jù)的差值構(gòu)造目標(biāo)函數(shù)；然后使用L－BFGS優(yōu)化算法對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，并求出目標(biāo)函數(shù)取得極小值時(shí)對(duì)應(yīng)的速度模型；接下來(lái)將該速度模型作為相鄰較高頻帶數(shù)據(jù)反演的初始速度模型，使用上述相同的方法求出該頻帶對(duì)應(yīng)的速度模型；這樣反復(fù)進(jìn)行頻帶循環(huán)計(jì)算，直到達(dá)到最高頻帶循環(huán)次數(shù)為止，最后得出的速度模型即為最終反演得到的速度模型。

2 數(shù)值模擬試驗(yàn)

為了測(cè)試本文提出的全波形反演方法的性能，分別對(duì)復(fù)雜Marmousi模型、高速楔形體模型及逆沖斷層模型進(jìn)行了速度反演、抗噪能力研究和反演速度擾動(dòng)能力研究。在各模型的試算中選用了以下相同的參數(shù)：震源子波為雷克子波，主頻12Hz；接收時(shí)間長(zhǎng)度為2.4s；頻率采樣間隔為0.42Hz，反演了2.92～25.41Hz間的10個(gè)頻帶，共55個(gè)頻率采樣點(diǎn)，其中每個(gè)頻帶包含10個(gè)頻率，每個(gè)頻帶最多迭代計(jì)算30次；單炮震源數(shù)等于相應(yīng)模型橫向網(wǎng)格點(diǎn)數(shù)，并利用這些單炮組合成20個(gè)同時(shí)激發(fā)震源；檢波器個(gè)數(shù)為相應(yīng)模型橫向網(wǎng)格點(diǎn)數(shù)；檢波器和震源均位于模型的最頂部。

引入擬合誤差來(lái)定量刻畫(huà)速度反演的正確程度，擬合誤差W定義為

其中：N為模型總道數(shù)；M為每道采樣點(diǎn)數(shù)；v′ij為速度反演結(jié)果；vij為理論速度值；i、j則分別為采樣點(diǎn)號(hào)和道號(hào)。

2.1 Marmousi模型

通過(guò)對(duì) Marmousi縱波速度模型（圖2a）的反演，來(lái)驗(yàn)證并分析筆者提出的全波形反演方法。模型網(wǎng)格為128×384，網(wǎng)格間距24m。原始Marmousi速度模型大尺度平滑后的速度設(shè)定為初始速度模型（圖2b）。對(duì)初始速度模型從低頻到高頻逐次反演10個(gè)頻帶，得到最終速度模型（圖2c）。

圖2 Marmousi模型全波形反演Fig.2 Full waveform inversion of Marmousi model

對(duì)于該模型，采用10個(gè)頻帶的單炮正演一次所需要的時(shí)間大約為85s，同時(shí)激發(fā)震源（384炮）正演一次所需時(shí)間大約為700s；因此，單炮震源正演完384炮所需的時(shí)間32 640s，要遠(yuǎn)遠(yuǎn)大于同時(shí)激發(fā)震源正演所需要的時(shí)間。所以，相對(duì)于常規(guī)單炮全波形反演，筆者的全波形反演方法在計(jì)算效率上得到了很大提高，并且在計(jì)算過(guò)程中可以明顯地發(fā)現(xiàn)在對(duì)計(jì)算機(jī)內(nèi)存的占用方面也得到了很大的改善。對(duì)比圖2a和圖2c可以看出，用筆者的全波形反演方法通過(guò)有限頻帶反演得到的速度模型與實(shí)際模型非常接近，反演效果非常好。利用式（14）計(jì)算最終反演結(jié)果與實(shí)際Marmousi模型的擬合誤差較?。╓＝0.095 9），定量地說(shuō)明了反演得到的速度精度較高。另外，在研究過(guò)程中發(fā)現(xiàn)，初始頻段的選取和初始速度模型的選取是一種矛盾的關(guān)系：即若選取的初始頻段頻率較高，反演方法對(duì)初始速度模型的要求就比較高，初始速度模型比較準(zhǔn)確時(shí)才能得到比較好的反演結(jié)果；相反，若選取的初始頻段頻率較低，反演方法對(duì)初始速度模型的要求降低，比較差一點(diǎn)的初始速度模型也能得到比較好的反演結(jié)果。

2.2 高速楔形體模型

利用信噪比（SNR）對(duì)壓制隨機(jī)噪聲的質(zhì)量進(jìn)行評(píng)價(jià)，定義式［17］如下：

式中：S0為原始模型正演地震記錄；S1為含噪聲地震記錄或壓制噪聲后的地震記錄。

為了檢驗(yàn)全波形反演方法對(duì)噪聲干擾的適應(yīng)程度，同時(shí)也為了問(wèn)題的說(shuō)明方便，筆者利用簡(jiǎn)單高速楔形體對(duì)含噪聲模型進(jìn)行試算。實(shí)際模型是在均勻介質(zhì)中包含一個(gè)高速的楔形體，均勻介質(zhì)和高速楔形體的縱波速度分別為3 500m/s和4 000m/s（圖3a），模型網(wǎng)格為160×320，網(wǎng)格間距20m。以3 500m/s的均勻背景速度場(chǎng)作為初始速度模型，將加噪聲的地震記錄數(shù)據(jù)用于筆者的全波形反演方法中，圖3b為最終反演結(jié)果。圖4a是對(duì)截取的第10個(gè)同時(shí)激發(fā)震源作用于原始模型的時(shí)間域正演地震記錄加入高斯隨機(jī)白噪生成的含噪地震記錄，SNR為11.147 3dB。圖4b為加噪數(shù)據(jù)最終反演結(jié)果的同時(shí)激發(fā)震源時(shí)間域正演地震記錄。從圖3可以看到，反演的楔形體很清楚，與實(shí)際速度模型非常接近。利用式（14）計(jì)算最終反演結(jié)果與實(shí)際速度模型的擬合誤差較?。╓＝0.024 1）；從而定量地說(shuō)明了該反演方法反演得到的速度精度很高。由圖4可以明顯看出，隨機(jī)噪聲得到了壓制。根據(jù)式（15）計(jì)算可知，信噪比提高為22.251 8dB。所以筆者提出的全波形反演方法具有一定的抗噪能力，對(duì)于提高地震資料信噪比也是比較有效的。

圖3 高速楔形體模型全波形反演Fig.3 Full waveform inversion of high－speed wedge model

圖4 高速楔形體模型同時(shí)激發(fā)震源時(shí)間域正演地震記錄Fig.4 Simultaneous source time domain forward seismic records of high－speed wedge model

2.3 逆沖斷層模型

在經(jīng)典的波場(chǎng)傳播研究中，通常將實(shí)際介質(zhì)簡(jiǎn)化為均勻介質(zhì)或?qū)訝罹鶆蚪橘|(zhì)，常規(guī)的地震勘探工作即建立在該假設(shè)的基礎(chǔ)上。然而實(shí)際上，小尺度范圍內(nèi)的速度和密度的擾動(dòng)廣泛存在于地殼之中，并通過(guò)向各個(gè)方向上散射入射波的能量，共同影響人工地震數(shù)據(jù)的質(zhì)量。所以，通常建立的均勻介質(zhì)模型并不符合實(shí)際地下情況，需要對(duì)這些小尺度的異常加以考慮。在地震勘探中，介質(zhì)的隨機(jī)擾動(dòng)可以理解為產(chǎn)生隨機(jī)噪聲的一種來(lái)源，所以它像隨機(jī)噪聲一樣，對(duì)勘探精度有很大影響。為了考慮這些細(xì)小但是無(wú)法忽略的異常，筆者使用全波形反演方法對(duì)只具有速度擾動(dòng)的隨機(jī)介質(zhì)逆沖斷層模型進(jìn)行試算。自相關(guān)函數(shù)、尺度、方差和粗糙系數(shù)是隨機(jī)介質(zhì)的4個(gè)要素［18］。介質(zhì)隨機(jī)擾動(dòng)可以分為3種情況：1）以隨機(jī)介質(zhì)為目標(biāo)體；2）目標(biāo)體被隨機(jī)介質(zhì)遮擋；3）上述兩種情況的復(fù)合。筆者選用的是第三種情況。

逆沖斷層縱波速度模型如圖5a所示，模型網(wǎng)格為187×801，網(wǎng)格間距25m。在圖5a模型上加入相關(guān)長(zhǎng)度為200m×20m、方差為3%、粗糙系數(shù)為2的Von Karman型隨機(jī)介質(zhì)，生成只具有速度擾動(dòng)的隨機(jī)介質(zhì)模型（如圖5b所示，為方便對(duì)比，淺部450m以上沒(méi)有加入速度擾動(dòng)）。對(duì)圖5b模型做大尺度平滑后得到初始速度模型（圖5c），并利用筆者提出的全波形反演方法進(jìn)行反演，得到最終速度模型（圖5d）。從圖5d可以看到，反演得到的最終速度模型很清楚，與具有速度擾動(dòng)特性的實(shí)際速度模型（圖5b）非常接近。利用式（14）計(jì)算最終反演結(jié)果與實(shí)際速度模型的擬合誤差W＝0.036 0，也定量說(shuō)明了反演得到的速度精度較高，同時(shí)驗(yàn)證了筆者所提出的全波形反演方法具有反演速度擾動(dòng)（隨機(jī)介質(zhì)）能力。

3 結(jié)論與展望

1）模型試算的結(jié)果表明：基于L－BFGS算法的頻率多尺度全波形反演可以通過(guò)有限的頻帶反演出精度很高的速度；同時(shí)也證明了筆者提出的將LBFGS數(shù)值優(yōu)化算法與同時(shí)激發(fā)震源技術(shù)相結(jié)合的方法有效地減少了內(nèi)存開(kāi)銷，提高了計(jì)算效率，適用于求解全波形反演問(wèn)題。

2）分析了隨機(jī)噪聲對(duì)筆者提出的全波形反演方法的影響，通過(guò)模型試算可知，該方法具有一定的抗噪能力，表現(xiàn)出較好的可靠性、有效性和穩(wěn)定性。

圖5 逆沖斷層模型全波形反演Fig.5 Full waveform inversion of thrust fault model

3）介質(zhì)中的隨機(jī)擾動(dòng)對(duì)高分辨率地震勘探有較大影響，筆者嘗試了用全波形反演方法反演隨機(jī)介質(zhì)，將隨機(jī)場(chǎng)理論引入到全波形反演中，充分考慮了地下介質(zhì)的隨機(jī)特性，使反演結(jié)果更符合實(shí)際地質(zhì)情況。反演算法成功地將介質(zhì)的速度擾動(dòng)特征反演出來(lái)，且分辨率較高。

4）全波形反演是構(gòu)造高分辨率地層速度模型的有效手段，隨著計(jì)算機(jī)軟硬件技術(shù)的推進(jìn)、地震模擬技術(shù)的進(jìn)步和大規(guī)模優(yōu)化計(jì)算技術(shù)的發(fā)展，全波形反演必將在不遠(yuǎn)的將來(lái)被廣泛地應(yīng)用于生產(chǎn)實(shí)際。

（References）：

［1］Tarantola A.Inversion of Seismic Reflection Data in the Acoustic Approximation［J］.Geophysics，1984，49（8）：1259－1266.

［2］Tarantola A.A Strategy for Nonlinear Elastic Inversion of Seismic Reflection Data［J］.Geophysics，1986，51（10）：1893－1903.

［3］Marfurt K J.Accuracy of Finite－Difference and Finite－Element Modeling of the Scalar and Elastic Wave E－quation［J］.Geophysics，1984，49（5）：533－549.

［4］Pratt R G.Frequency－Domain Elastic Wave Modeling by Finite Differences：A Tool for Crosshole Seismic Imaging［J］.Geophysics，1990，55（5）：626－632.

［5］Pratt R G，Worthington M H.Inverse Theory Applied to Multi－Source Cross－Hole Tomography：Part 1：Acoustic Wave－Equation Method［J］.Geophysical Prospecting，1990，38（3）：287－310.

［6］Pratt R G.Inverse Theory Applied to Multi－Source Cross－Hole Tomography：Part 2：Elastic Wave－Equation Method［J］.Geophysical Prospecting，1990，38（3）：311－329.

［7］Pratt R G，Shin C，Hicks G J.Gauss－Newton and Full Newton Methods in Frequency－Space Seismic Waveform Inversion［J］.Geophysical Journal International，1998，133（2）：341－362.

［8］Song Z M，Williamson P R，Pratt R G.Frequency－Domain Acoustic－Wave Modeling and Inversion of Crosshole Data：Part 2：Inversion Method，Synthetic Experiments and Real－Data Results［J］.Geophysics，1995，60（3）：796－809.

［9］Sourbier F，Operto S，Virieux J，et al.FWT2D：A Massively Parallel Program for Frequency－Domain Full－Waveform Tomography of Wide－Aperture Seismic Data：Part 1：Algorithm［J］.Computers and Geosciences，2009，35（3）：487－495.

［10］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.

［11］Nocedal J.Updating Quasi－Newton Matrices with Limited Storage［J］.Mathematics of Computation，1980，35：773－782.

［12］Qin Z W（Tony）.The Relationships Between CG，BFGS，and Two Limited－Memory Algorithms［J］.Electronic Journal of Undergraduate Mathematics，2007，12：5－20.

［13］韓利，韓立國(guó)，李翔，等.二階聲波方程頻域PML邊界條件及頻域變網(wǎng)格步長(zhǎng)并行計(jì)算［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2011，41（4）：1226－1232.Han Li，Han Liguo，Li Xiang，et al.PML Boundary Conditions for Second－Order Acoustic Wave Equations and Variable Grid Parallel Computation in Frequency－Domain Modeling［J］.Journal of Jilin University：Earth Science Edition，2011，41（4）：1226－1232.

［14］Berkhout A J.Changing the Mindset in Seismic Data Acquisition［J］.The Leading Edge，2008，27（7）：924－938.

［15］Neelamani R，Krohn C E，Krebs J R，et al.Efficient Seismic Forward Modeling Using Simultaneous Random Sources and Sparsity［J］.Geophysics，2010，75（6）：WB15－WB27.

［16］Liu D C，Nocedal J.On the Limited Memory BFGS Method for Large Scale Optimization［J］.Mathematical Programming，1989，45（1）：503－528.

［17］張亞紅.Focal變換及其在地震數(shù)據(jù)去噪和插值中的應(yīng)用研究［D］.長(zhǎng)春：吉林大學(xué)，2010.Zhang Yahong.Focal Transformation Research and Application in Seismic Data Denoising and Interpolation［D］.Changchun：Jilin University，2010.

［18］王恩利.基于各向異性的復(fù)合介質(zhì)彈性波場(chǎng)模擬與特征分析［D］.長(zhǎng)春：吉林大學(xué)，2008.Wang Enli.The Forward Modeling and Character A－nalysis for Complex Media Based on Anisotropic Media［D］.Changchun：Jilin University，2008.