張生強，劉春成，韓立國，楊小椿

1.吉林大學地球探測科學與技術學院，長春 130026

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

0 引言

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

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

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

1 方法原理

頻率域波動方程全波形反演是利用波形，即波場的振幅、頻率以及相位信息，將數據殘差反傳波場與初始模型的入射波場互相關，更新地下介質參數，進而對地下介質物性參數成像的方法。其具有揭示復雜地質背景下構造與巖性細節(jié)信息的潛力?；诓▌臃匠痰娜ㄐ畏囱荼旧硎莻€強非線性反問題，解這類問題的經典方法是首先將問題線性化，再通過逐步迭代最終確定一個使理論計算值與觀測記錄較好擬合的反演結果。從數學的角度來講，全波形反演可以對應地分為2個問題：大型稀疏矩陣線性方程組的求解以及非線性局部最優(yōu)化問題。

1.1 頻率域正演

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

1.1.1 Helmholtz方程

對時域二階標量聲波方程作關于時間的Fourier變換可得到頻率域標量聲波方程，即Helmholtz方程［13］：

其中：▽2為 Laplace算子；A＝－▽2－k2為Helmholtz算子；u為地震波場；k＝ω/v為波數，ω為角頻率，v為介質速度；S為震源項。

采用LU分解法求解方程（1）即可求得地震波場u。地震波場u和模型參數m之間的關系是非線性的，可以通過算子D（D為m和u之間的非線性函數）來定義：

1.1.2 同時激發(fā)震源

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

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

超炮道集的構建可由下式表達：

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

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

1.2 全波形反演

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

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

其中，＊表示共軛轉置算子。

1.2.1 波恩近似及反演問題的線性化

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

在初始速度模型v0點附近對目標函數進行二階Taylor－Lagrange展開：

其中：整數M表示速度模型v中的元素個數；O（v3）表示誤差項。式（7）兩端同時對速度模型參數vl求導并寫成矩陣形式：

當目標函數的一階導數為0時，則目標函數在v0附近達到極小值，同時得到速度模型擾動的表達式：

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

1.2.2 L－BFGS算法

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

其中：函數C∈R，是實值函數，稱C（v）為問題（10）的目標函數；集合Ω為問題的定義域。

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

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

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

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

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

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

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

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

第5步，若k＝k＋1，轉第2步。

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

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

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

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

2 數值模擬試驗

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

引入擬合誤差來定量刻畫速度反演的正確程度，擬合誤差W定義為

其中：N為模型總道數；M為每道采樣點數；v′ij為速度反演結果；vij為理論速度值；i、j則分別為采樣點號和道號。

2.1 Marmousi模型

通過對 Marmousi縱波速度模型（圖2a）的反演，來驗證并分析筆者提出的全波形反演方法。模型網格為128×384，網格間距24m。原始Marmousi速度模型大尺度平滑后的速度設定為初始速度模型（圖2b）。對初始速度模型從低頻到高頻逐次反演10個頻帶，得到最終速度模型（圖2c）。

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

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

2.2 高速楔形體模型

利用信噪比（SNR）對壓制隨機噪聲的質量進行評價，定義式［17］如下：

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

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

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

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

2.3 逆沖斷層模型

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

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

3 結論與展望

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

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

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

3）介質中的隨機擾動對高分辨率地震勘探有較大影響，筆者嘗試了用全波形反演方法反演隨機介質，將隨機場理論引入到全波形反演中，充分考慮了地下介質的隨機特性，使反演結果更符合實際地質情況。反演算法成功地將介質的速度擾動特征反演出來，且分辨率較高。

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

（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］韓利，韓立國，李翔，等.二階聲波方程頻域PML邊界條件及頻域變網格步長并行計算［J］.吉林大學學報：地球科學版，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變換及其在地震數據去噪和插值中的應用研究［D］.長春：吉林大學，2010.Zhang Yahong.Focal Transformation Research and Application in Seismic Data Denoising and Interpolation［D］.Changchun：Jilin University，2010.

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