吳 彥,馬 玥,劉玉金,駱 毅

(1.沙特阿美北京研發(fā)中心,北京100102;2.沙特阿美EXPEC高級研究中心,沙特阿拉伯達(dá)蘭)

全走時反演及其應(yīng)用

吳 彥1,馬 玥1,劉玉金1,駱 毅2

(1.沙特阿美北京研發(fā)中心,北京100102;2.沙特阿美EXPEC高級研究中心,沙特阿拉伯達(dá)蘭)

波動方程速度建模方法是利用地震數(shù)據(jù)或成像道集中的旅行時信息或運動學(xué)特征自動反演出背景速度模型,但是已有方法普遍面臨一個嚴(yán)重的問題,即反演過程無法避免受振幅信息的影響。為此,提出了一種基于波動方程的旅行時反演方法,稱之為全走時(即完全依賴于走時信息)反演,從地震數(shù)據(jù)中自動估計出運動學(xué)意義上正確的速度模型,其核心思想是使反演完全依賴于走時信息,從而防止反演受到振幅信息的干擾。基于速度擾動只產(chǎn)生旅行時變化這一假設(shè),介紹了數(shù)據(jù)域和成像域兩種全走時反演方法,分別對應(yīng)透射波和反射波走時反演。全走時反演不需要準(zhǔn)確的初始速度模型和地震數(shù)據(jù)中的低頻信息。理論模型和實際資料測試結(jié)果表明,全走時反演在常速度初始模型的情況下也能得到令人滿意的反演結(jié)果。

波動方程;速度建模;全自動;旅行時反演

全波形反演(FWI)根據(jù)觀測地震數(shù)據(jù)構(gòu)建地下速度模型,利用了地震波波形中的所有信息,包括振幅和相位信息。常規(guī)的全波形反演[1-3]和反射波波形反演[4-5]都是根據(jù)最小化觀測數(shù)據(jù)和正演數(shù)據(jù)之間的波形差異,沿著波路徑反傳數(shù)據(jù)誤差,得到模型更新量。FWI最大的問題是相減型目標(biāo)函數(shù)與速度模型呈強(qiáng)非線性關(guān)系。只要預(yù)測數(shù)據(jù)和觀測數(shù)據(jù)中存在相位差異大于半個周期的情況,FWI就會遇到嚴(yán)重的周波跳躍問題。在這種情況下,反演容易陷入局部極小值。從小偏移距到大偏移距進(jìn)行FWI是解決周波跳躍問題的一種實用化方法;另外,從低頻到高頻進(jìn)行FWI也能夠緩解這一問題。為了避免觀測波形和正演波形之間存在太大的旅行時差異,仍然需要一個盡可能準(zhǔn)確的初始速度模型。

相對于波形,旅行時和長波長速度變化之間更符合線性關(guān)系,因此工業(yè)界常采用旅行時層析反演地下結(jié)構(gòu)信息。地震層析只利用了旅行時信息,但是能夠提供滿足構(gòu)造成像要求的速度模型,也即運動學(xué)意義上正確的速度模型。射線層析是工業(yè)界最常用的旅行時反演方法,該方法通過旅行時信息反演得到平滑的偏移速度模型。在射線層析中,需要在偏移后的共成像點道集拾取剩余時差,然后將拾取的剩余時差沿著射線路徑反傳得到速度更新量。射線層析方法計算效率高，當(dāng)?shù)卣鸩▊鞑ミ^程可用射線理論近似時,利用射線層析反演能夠得到合理的速度模型[6-8]。但是,在復(fù)雜地質(zhì)構(gòu)造條件下,由于射線多路徑的影響及射線理論高頻近似的局限,射線層析方法很難反演出正確的速度模型。

通過引入有限頻地震波傳播算子,波動方程旅行時層析可以克服射線層析的多路徑問題。在處理運動學(xué)頻散和采集孔徑有限的數(shù)據(jù)時,波動方程層析優(yōu)于射線層析[7]。波動方程層析中一個關(guān)鍵的步驟是提取正演數(shù)據(jù)和觀測數(shù)據(jù)之間的旅行時差異。在數(shù)據(jù)域中,最直接的旅行時時差估算方法是在波形數(shù)據(jù)中拾取兩個相關(guān)的波形,然后選擇觀測波形和計算波形互相關(guān)最大時所對應(yīng)的時移量[9]。在成像域中,旅行時時差可以在共成像點道集(CIGs)中拾取同相軸的剩余時差或剩余深度差[10-11]。上述旅行時反演方法的目標(biāo)函數(shù)是最小化拾取的相對旅行時時差。盡管工業(yè)界存在多種自動拾取技術(shù),但是在多數(shù)情況下,為了降低自動拾取帶來的誤差,質(zhì)量控制和人工干預(yù)仍然必不可少。

為了避免拾取帶來的誤差,VAN LEEUWEN等[12]采用互相關(guān)加權(quán)范數(shù)作為目標(biāo)函數(shù),該目標(biāo)函數(shù)隱式地度量正演數(shù)據(jù)和觀測數(shù)據(jù)之間的旅行時差,不需要做任何拾取,就可以進(jìn)行全自動反演。但是,由于觀測數(shù)據(jù)中不同反射波到達(dá)時之間容易發(fā)生串?dāng)_,因此這個目標(biāo)函數(shù)對反射波數(shù)據(jù)容易出現(xiàn)多走時問題。一些成像域的速度建模方法,如波動方程偏移速度分析[13],可以通過偏移成像克服這個問題。這類方法通過最小化成像擾動,使得成像空間中的共成像點道集拉平或者地下偏移距道集聚焦在零偏移距。在波動方程偏移速度分析類方法中,微分相似度最優(yōu)化(DSO)方法將震源波場及接收波場互相關(guān)之后的結(jié)果乘以地下偏移距作為誤差函數(shù)[14]。DSO通過隱式度量地下偏移距道集的聚焦性或者角道集的拉平程度自動更新速度模型。但是,常規(guī)的DSO梯度往往出現(xiàn)明顯的近垂直的條帶狀噪聲,這類噪聲會降低反演的收斂速度[15-18]。通過改進(jìn)成像剩余量的提取方法能夠避免這種噪聲,比如水平收縮算子提取無限小的成像擾動[18],但是這類方法的速度更新量已經(jīng)不是原目標(biāo)函數(shù)的梯度。2009年,ZHANG等[19]提出了一種基于旅行時的波動方程反演方法,該方法反傳延遲炮記錄,并通過最大化反傳波場在震源位置處對應(yīng)激發(fā)時間的能量,更新速度模型,它只依賴于地震數(shù)據(jù)中的走時信息,而對地震數(shù)據(jù)中存在的振幅誤差不敏感。2011年,LUO等[20]提出反褶積型波動方程反演,該方法最小化正演數(shù)據(jù)和觀測數(shù)據(jù)之間加權(quán)反褶積的范數(shù),全自動,對振幅誤差不敏感,而且對帶限(或非脈沖)震源函數(shù)也不敏感。2014年,WARNER等[21]提出自適應(yīng)波形反演(AWI),該方法首先利用維納濾波器匹配預(yù)測數(shù)據(jù)和觀測數(shù)據(jù),然后通過最小化維納濾波器非零延遲系數(shù)的能量,更新速度模型。如果維納濾波器和輸入數(shù)據(jù)長度相同,AWI和反褶積型波動方程反演類似。實際上,設(shè)計一個正演數(shù)據(jù)和觀測數(shù)據(jù)每個地震道都完全匹配的維納濾波器十分困難。

本文介紹一種全自動的波動方程旅行時反演方法,該方法能夠利用旅行時信息反演得到運動學(xué)意義上準(zhǔn)確的速度模型。首先介紹透射波和反射波情形下全走時反演(FTI)的基本原理,并分析FTI的梯度和Born近似下的梯度存在的差異,然后通過理論模型和實際資料測試,驗證FTI方法的正確性和有效性。

1 基本原理

1.1 透射波FTI

首先介紹數(shù)據(jù)域透射波FTI的基本原理。為方便起見,本文假定:①地下模型為常密度聲介質(zhì)模型；②假定模型擾動Δs只產(chǎn)生地震波場的時間變化Δτ?；诘冖趥€假設(shè),地震波同相軸可以完全由旅行時來表征。令Pobs(xr,t;xs)和Pcal(xr,t;xs)分別表示炮點位于xs,接收點位于xr,時刻為t的觀測地震道和計算地震道。定義兩道之間的互相關(guān)函數(shù):

(1)

式中:τ為觀測地震道和計算地震道之間的時移量。目標(biāo)函數(shù)定義為[12]:

(2)

式中:τ∈[-T,T]是線性加權(quán)函數(shù),用于壓制互相關(guān)c(xr,τ;xs)在非零時移處的能量,其中,T是估算觀測地震道和計算地震道之間旅行時時差的最大值。最小化目標(biāo)函數(shù)((2)式)意味著選擇速度模型使得互相關(guān)能量聚焦在零時移處。與LUO等[9]提出的互相關(guān)波動方程旅行時反演方法相比,該方法不需要顯式拾取互相關(guān)時移量,而是自動計算出來,且對振幅譜的差異也更加不敏感[12]。

經(jīng)過一系列的公式推導(dǎo)[22],目標(biāo)函數(shù)((2)式)關(guān)于模型的Born梯度為:

(3)

這里,g(x′,t;xr)為格林函數(shù)，表示單位源在地表xr處激發(fā)后在地下x′處所觀測的波場。Born近似表明速度擾動和波場擾動之間呈線性關(guān)系。波場擾動包含了振幅和相位信息。但是,互相關(guān)目標(biāo)函數(shù)((2)式)用于度量相對時移量,和基于Born近似的梯度公式((3)式)的意義不一致,這會導(dǎo)致常規(guī)的Born近似最優(yōu)化算法無法獲得目標(biāo)函數(shù)所對應(yīng)的正確梯度,使得反演過程收斂很慢或者不收斂。

目標(biāo)函數(shù)隱式地最小化拾取觀測和計算的透射波到達(dá)時差。當(dāng)梯度和目標(biāo)函數(shù)意義一致時,目標(biāo)函數(shù)的這種誤差度量方式才能產(chǎn)生有意義的模型更新量。

假設(shè)模型擾動Δs只產(chǎn)生波傳播的旅行時變化Δτ。,對于單個炮點-檢波點對,當(dāng)慢度模型從s0變化到s0+Δs,計算地震道Pcal(xr,t;xs)s0僅存在時間變化Δτ,即:

(4)

利用公式(4)并經(jīng)過一系列公式推導(dǎo)[22],得到目標(biāo)函數(shù)((2)式)關(guān)于模型的FTI梯度:

(5)

與Born近似下的梯度公式((3)式)相比,FTI梯度僅包含旅行時信息,這和目標(biāo)函數(shù)的誤差度量方式一致。注意,對單個炮點-檢波點對,伴隨源(或FTI反傳剩余量)為:

(6)

基于速度擾動只產(chǎn)生旅行時變化的假設(shè),可以將公式中的Pcal(xr,t;xs)替換為Pobs(xr,t+Δτ;xs),其中，Δτ為合成地震道匹配觀測地震道的時移量。伴隨源可以改寫為:

(7)

實際上,滿足方程(7)的條件是地震子波是Delta函數(shù),這樣對任意的τ≠Δτ都滿足c(xr,τ;xs)=0。在此情況下,積分項中只有τ=Δτ時不為0。和直接拾取時移量Δτ不同,伴隨源是通過最大化互相關(guān)函數(shù)c(xr,τ;xs)計算出“最佳”的時移量。

1.2 反射波FTI

反射波FTI在成像空間可以由透射波的波路徑來表征,其目標(biāo)是將擴(kuò)展成像結(jié)果聚焦在反射層深度上。地震波場由平滑的背景速度模型和反射系數(shù)正演得到。這種方式預(yù)測的地震數(shù)據(jù)包含所有的一次反射波,而不含任何多次波。

震源波場Ps(x,t;xs)和反向傳播得到的接收波場Pr(x,t;xs)在成像點x位置處的時移互相關(guān)[23]為:

I(x,τ)=?Ps(x,t+τ;xs)Pr(x,t;xs)dtdxs

(8)

其中,τ(x,xs)∈[-T,T]是時移量。沿著時移坐標(biāo)軸方向?qū)ζ瞥上襁M(jìn)行擴(kuò)展,可以線性模擬速度擾動產(chǎn)生的時移量。時移量彌補(bǔ)了速度誤差并使得成像結(jié)果在非零時移處形成反射能量。正確模型下的成像結(jié)果聚焦在零時移處。為了自動度量反射能量在不同深度位置的時移量,定義目標(biāo)函數(shù)為:

(9)

時移量τ用于壓制非零時移處的能量。最小化上述目標(biāo)函數(shù)意味著選擇速度模型使得擴(kuò)展成像結(jié)果聚焦在零延遲時刻。

對于DSO方法,速度更新量(梯度)表示為:

(10)

公式(10)的詳細(xì)推導(dǎo)過程參考文獻(xiàn)[14]。在原來的DSO方法中,擴(kuò)展成像道集定義在地下偏移距域。為簡單起見,這里僅分析時移成像道集。地下空移量和時移量可以相互轉(zhuǎn)換[24]。與公式(5)類似,DSO梯度包含了振幅和相位信息,在Born近似下無法解耦。但是目標(biāo)函數(shù)度量的是相對時移量,因此,DSO梯度和目標(biāo)函數(shù)的意義不一致。FEI等[15]指出反射層截斷或畸變會對DSO梯度引入噪聲,降低收斂速度。SHEN等[18]提出了一種地下偏移距域水平收縮方法,可以產(chǎn)生無噪聲的速度更新量。該方法通過水平校正算子無限小地改進(jìn)當(dāng)前成像結(jié)果。成像擾動則是校正后的和原來的成像結(jié)果之差。盡管由此生成的速度更新量不會產(chǎn)生噪聲,但是它已經(jīng)不是原來DSO目標(biāo)函數(shù)的梯度。

和透射波情形類似,FTI假定模型擾動只產(chǎn)生時移成像道集聚焦能量的時間變化,即:

(11)

和數(shù)據(jù)域推導(dǎo)類似,目標(biāo)函數(shù)(9)式關(guān)于慢度的梯度公式[22]為:

(12)

其中,分母E為:

(13)

旅行時剩余量Δτ度量的是聚焦性,可以在時移成像道集中拾取得到。由于成像道集I(x,τ)只聚焦在Δτ位置,因此公式(12)在任何τ≠Δτ處生成的波路徑可以忽略,即:

(14)

(14)式不需要拾取Δτ,而是采用積分自動提取Δτ對應(yīng)的主要能量。

2 數(shù)值試驗

2.1 Wedge模型

首先采用二維理論模型對本文方法進(jìn)行測試。圖1a是真實速度模型。采用線性Born正演方法合成反射地震數(shù)據(jù)。震源是主頻為10Hz的雷克子波,共正演126炮,炮間距為80m,每炮以固定接收排列方式接收數(shù)據(jù),共251道,道間距為40m。初始速度模型為常速(2000m/s)模型。圖1b為FTI反演迭代10次后的速度模型。對比圖1a和圖1b可以看出,反演結(jié)果較好地恢復(fù)了真實速度模型的變化趨勢。圖2各個子圖分別顯示初始速度模型、反演速度模型和真實速度模型下的疊前深度偏移成像結(jié)果和地下偏移距道集。對比不同速度模型下的成像結(jié)果和道集可以看出,采用反演后的速度模型,反射層基本偏移到真實的空間位置,地下偏移距道集聚焦在零偏移距處。理論模型測試結(jié)果表明:FTI從常速度模型開始反演,得到的背景速度模型,可以滿足構(gòu)造成像的精度要求。

圖1 真實速度模型(a)和FTI反演速度模型(b)

2.2 實際資料

采用某二維實際地震數(shù)據(jù)對本文方法進(jìn)行測試。該實際資料共801炮,炮間距25m,每炮以固定接收排列方式接收數(shù)據(jù),共1026道,道間距25m。采用常速模型(4500m/s)作為初始速度模型進(jìn)行反演。圖3是FTI反演迭代10次后的速度模型。圖4a和圖4b分別為初始速度模型下的疊前深度偏移成像剖面和地下偏移距道集。圖5a和圖5b分別為FTI反演速度模型下的疊前深度偏移成像剖面和地下偏移距道集。對比圖4和圖5可以看出,采用FTI反演得到的速度模型,可以改善成像質(zhì)量,提高偏移成像道集的聚焦程度。

圖2 對Wedge模型數(shù)據(jù)用不同速度模型偏移的成像剖面(左)和地下偏移距成像道集(右)a 初始速度模型;b 反演速度模型;c 正確速度模型

圖3 FTI反演迭代10次后的速度模型

圖4 對實際數(shù)據(jù)用初始速度模型偏移的成像剖面(a)和地下偏移距成像道集(b)

圖5 對實際數(shù)據(jù)用反演速度模型偏移的成像剖面(a)和地下偏移距成像道集(b)

3 結(jié)論

FTI是一種利用旅行時信息構(gòu)建運動學(xué)意義上正確速度模型的反演方法，其核心思想是減少振幅對旅行時反演的影響。FTI可以在數(shù)據(jù)域也可以在成像域?qū)崿F(xiàn),分別用于接收端的透射波和反射層上的反射波旅行時反演。數(shù)據(jù)域全走時反演的目標(biāo)函數(shù)是觀測和計算的透射波數(shù)據(jù)之間互相關(guān)的線性加權(quán)范數(shù)。成像域FTI的目標(biāo)函數(shù)是時間或者空間延遲擴(kuò)展成像道集的線性加權(quán)范數(shù)。兩個目標(biāo)函數(shù)都可以隱式地度量旅行時相對誤差?；谒俣葦_動只產(chǎn)生波場的旅行時變化這一假設(shè),本文分別介紹了透射波和反射波情形下的FTI梯度(波路徑)。與Born近似意義下的梯度相比,FTI梯度只包含旅行時信息,和目標(biāo)函數(shù)極小化旅行時時差的意義一致。FTI不需要準(zhǔn)確的初始速度模型和地震數(shù)據(jù)中的低頻信息。理論模型和實際資料測試結(jié)果表明,FTI具有很好的收斂性,采用十分簡單的初始速度模型,在很少的迭代次數(shù)內(nèi)也能收斂到滿足構(gòu)造成像精度的要求。

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

[2] PRATT R G,SHIN C,HICKS G J.Gauss-Newton and full Newton methods in frequency domain seismic waveform inversion[J].Geophysical Journal International,1998,133(2):341-362

[3] VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[4] XU S,WANG D,CHEN F,et al.Inversion on reflected seismic wave[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-7

[5] BROSSIER R,OPERTO S,VIRIEUX J.Velocity model building from seismic reflection data by full-waveform inversion[J].Geophysical Prospecting,2015,63(2):354-367

[6] MENG Z,VALASEK P A,WHITNEY S A,et al.3D global tomographic velocity model building[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:2379-2382

[7] WOODWARD M,NICHOLS D,ZDRAVEVA O,et al.A decade of tomography[J].Geophysics,2008,73(5):VE5-VE11

[8] JONES I F.Tutorial:velocity estimation via ray-based tomography[J].First Break,2010,28(2):45-52

[9] LUO Y,SCHUSTER G T.Wave-equation traveltime inversion[J].Geophysics,1991,56(5):645-653

[10] YANG T,SAVA P.3D image-domain wavefield tomography using time-lag extended images[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:4816-4821

[11] ZHANG S,SCHUSTER G,LUO Y.Wave-equation reflection traveltime inversion[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2705-2710

[12] VAN LEEUWEN T,MULDER W A.A correlation-based misfit criterion for wave-equation traveltime tomography[J].Geophysical Journal International,2010,182(3):1383-1394

[13] SAVA P,BIONDI B.Wave-equation migration velocity analysis I:theory[J].Geophysical Prospecting,2004,52(6):593-607

[14] SHEN P,SYMES W W.Automatic velocity analysis via shot profile migration[J].Geophysics,2008,73(5):VE49-VE59

[15] FEI W,WILLIAMSON P.On the gradient artifacts in migration velocity analysis based on differential semblance optimization[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:4071-4076

[16] VYAS M,TANG Y.Gradients for wave-equation migration velocity analysis[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:4077-4081

[17] ALBERTIN U.An improved gradient computation for adjoint wave equation reflection tomography[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:3969-3973

[18] SHEN P,SYMES W W.Horizontal contraction in image domain for velocity inversion[J].Geophysics,2015,80(3):R95-R110

[19] ZHANG Y,WANG D.Traveltime information-based wave-equation inversion[J].Geophysics,2009,74(6):WCC27-WCC36

[20] LUO S,SAVA P.A deconvolution-based objective function for wave-equation inversion[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2788-2792

[21] WARNER M,GUASCH L.Adaptive waveform inversion:theory[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:1089-1093

[22] LUO Y,MA Y,WU Y,et al.Full-traveltime inversion[J].Geophysics,2016,81(5):R261-R274

[23] SAVA P,FOMEL S.Time-shift imaging condition in seismic migration[J].Geophysics,2006,71(6):S209-S217

[24] SAVA P,FOMEL S.Angle-domain common-image gathers by wavefield continuation methods[J].Geophysics,2003,68(3):1065-1074

(編輯：顧石慶)

Full-traveltime inversion and its application

WU Yan1,MA Yue1,LIU Yujin1,LUO Yi2

(1.AramcoBeijingResearchCenter,AramcoAsia,Beijing100102,China;2.EXPECAdvancedResearchCenter,SaudiAramco,Dhahran,SaudiArabia)

Many previously published wave-equation based methods,which attempt to automatically invert traveltime or kinematic information in seismic data or migrated gathers for smooth velocities,suffer a common and severe problem,that is the inversion are involuntarily and unconsciously hijacked by amplitude information.To overcome this problem,we present a new wave-equation based traveltime inversion methodology,referred to as full traveltime inversion (FTI),to automatically estimate a kinematically accurate velocity model from seismic data.The key idea of FTI is to make the inversion fully dependent on traveltime information,thus prevent amplitude interferenceduring the inversion.Under the assumption that velocity perturbations cause only traveltime changes,we derive the FTI method in both data domain and image domain which are applicable to transmitted arrivals and reflected waves,respectively.FTI does not require an accurate initial velocity model or the low frequencyseismicdata.Synthetic and field data tests demonstrate that FTI produces satisfactory inversion results,even using constant velocity models as initials.

wave equation,velocity modeling,automatic,traveltime inversion

2016-10-17;改回日期：2016-11-28。

吳彥(1985—),男,博士,主要從事地震速度建模及反演方法等研究。

馬玥(1984—),女,博士,主要從事超分辨率地震成像、壓縮感知、速度建模及反演方法等研究。

P631

A

1000-1441(2017)01-0050-07

10.3969/j.issn.1000-1441.2017.01.006