基于波場重建的時域全波形反演

李松齡

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

0 前 言

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

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

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

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

1 方法原理

1.1 全波形反演的梯度

時空域二階聲波常密度波動方程為：

(1)

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

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

(2)

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

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

(3)

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

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

(4)

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

(5)

(6)

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

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

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

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

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

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

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

1.2 全波形反演的解析步長

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

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

(7)

(8)

(9)

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

(10)

1.3 混合共軛梯度方法

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

(11)

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

(12)

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

vk+1=vk+αkdk

(14)

2 基于波場重建的時域FWI

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

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

波動方程差分公式

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

(15)

如圖1所示，假設(shè)真實波場范圍為圖1中的A區(qū)域，如果采用空間12階差分計算格式求解波動方程，僅需保存A區(qū)域外側(cè)6層波場數(shù)據(jù)，如圖1中B區(qū)域，即可滿足差分條件，C區(qū)域為吸收邊界區(qū)域。

圖1 邊界存儲示意圖

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

圖2 基于波場重建的時域FWI流程

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

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

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

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

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

表1 勻速模型波場重建存儲消耗對比

3 模型測試

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

3.1 波場重建方案測試

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

圖5 Marmousi2 速度模型

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

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

圖7 正傳及逆時逆時反傳的波場(t=0.45 s)單道對比

模型名稱模型大小/m總時間間隔/s全波場存儲/M邊界存儲/MMarmousi2500×20050001907160

3.2 Marmousi2全波形反演測試

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

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

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

圖8 Marmousi 平滑速度模型

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

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

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

4 結(jié) 論

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

本文方法對基于復(fù)雜方程(如一階速度應(yīng)力方程和彈性波方程等)的FWI方法或三維FWI均具有一定的借鑒意義，這也是本文下一步的研究重點；本文提出的波場重建策略，僅適用于非損耗介質(zhì)下的波場重建，對于粘性介質(zhì)的波場重建方法仍有待進(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] 王保利, 高靜懷, 陳文超,等. 地震疊前逆時偏移的有效邊界存儲策略[J]. 地球物理學(xué)報, 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.