時(shí)間二階積分波場的全波形反演

陳生昌， 陳國新

浙江大學(xué)地球科學(xué)學(xué)院， 杭州　310027

?

時(shí)間二階積分波場的全波形反演

陳生昌， 陳國新

浙江大學(xué)地球科學(xué)學(xué)院， 杭州310027

通過對波場的時(shí)間二階積分運(yùn)算以增強(qiáng)地震數(shù)據(jù)中的低頻成分，提出了一種可有效減小對初始速度模型依賴性的地震數(shù)據(jù)全波形反演方法—時(shí)間二階積分波場的全波形反演方法.根據(jù)散射理論中的散射波場傳播方程，推導(dǎo)出時(shí)間二階積分散射波場的傳播方程，再利用一階Born近似對時(shí)間二階積分散射波場傳播方程進(jìn)行線性化.在時(shí)間二階積分散射波場傳播方程的基礎(chǔ)上，利用散射波場反演地下散射源分布，再利用波場模擬的方法構(gòu)建地下入射波場，然后根據(jù)時(shí)間二階積分散射波場線性傳播方程中散射波場與入射波場、速度擾動(dòng)間的線性關(guān)系，應(yīng)用類似偏移成像的公式得到速度擾動(dòng)的估計(jì)，以此建立時(shí)間二階積分波場的全波形迭代反演方法.最后把時(shí)間二階積分波場的全波形反演結(jié)果作為常規(guī)全波形反演的初始模型可有效地減小地震波場全波形反演對初始模型的依賴性.應(yīng)用于Marmousi模型的全頻帶合成數(shù)據(jù)和缺失4 Hz以下頻譜成分的缺低頻合成數(shù)據(jù)驗(yàn)證所提出的全波形反演方法的正確性和有效性，數(shù)值試驗(yàn)顯示缺失4 Hz以下頻譜成分?jǐn)?shù)據(jù)的反演結(jié)果與全頻帶數(shù)據(jù)的反演結(jié)果沒有明顯差異.

散射波方程; 時(shí)間二階積分; 增強(qiáng)低頻; Born近似; 全波形反演

1　引言

利用地震數(shù)據(jù)獲取地下速度分布是地震勘探中的一個(gè)經(jīng)典問題，同時(shí)又是地震數(shù)據(jù)反演方法理論研究中的難點(diǎn).地震波場與地下速度模型之間是一種非線性的函數(shù)關(guān)系，利用觀測的地震波場反演地下速度模型的分布屬于典型的非線性反演問題.對于非線性反演問題的求解通常需要給定一個(gè)初始速度模型，由給定的速度模型和震源函數(shù)以及相應(yīng)的初邊值條件可計(jì)算出與速度模型對應(yīng)的計(jì)算波場，進(jìn)而獲得觀測波場與計(jì)算波場間的殘差波場.如果殘差波場小于設(shè)定的閾值，則認(rèn)為給定的速度模型就是利用地震波場反演出的地下速度模型.如果殘差波場大于設(shè)定的閾值，則需要對給定的速度模型進(jìn)行修正.根據(jù)殘差波場信息對初始速度模型進(jìn)行修正主要兩類方法，1)是非線性最優(yōu)化方法；2)是線性化迭代反演方法.第1類的非線性最優(yōu)化方法又包含全局最優(yōu)化方法和局部最優(yōu)化方法，由于得到計(jì)算波場的波場模擬計(jì)算量巨大，通常采用局部最優(yōu)化方法對速度模型進(jìn)行修正，這也就是目前廣泛研究和應(yīng)用的全波形反演方法.基于局部最優(yōu)化方法的全波形反演自上世紀(jì)80年代初提出以來(Lailly，1983；Tarantola，1984，1986,1987)，就一直是地震波場反演研究的熱點(diǎn)，并逐漸成為了地震勘探中的地下速度建模方法(Virieux and Operto，2009).本文把基于局部最優(yōu)化方法的全波形反演稱為常規(guī)全波形反演.第二類的全波形反演的線性化迭代反演方法主要基于地震波傳播的散射理論(Chew，1990)，利用漸近展開的方法(如Taylor展開、一階Born近似、Raytov近似等方法)把非線性的全波形反演問題線性化為迭代的線性反演問題，即在迭代中利用線性反演方法求解初始速度模型的修正.在數(shù)據(jù)擬合的最小平方準(zhǔn)則下，基于局部最優(yōu)化方法的全波形反演與基于線性化迭代反演方法的全波形反演是等價(jià)的(Vogel，2002).

常規(guī)全波形反演方法僅考慮地震波場與地下速度模型間的非線性關(guān)系，而沒有考慮地震波場傳播的物理過程，因此常規(guī)的全波形反演方法可視為一種純粹基于計(jì)算數(shù)學(xué)而發(fā)展起來的方法.在地震波場的傳播中，散射波場是最基本的波場，其它如反射、折射波場等是散射波場在高頻近似下的退化.地震散射波場傳播的物理過程包括，1)震源激發(fā)的入射波場在背景速度模型中傳播；2)入射波場與相對于背景速度的速度擾動(dòng)的相互作用產(chǎn)生散射波場，散射波場再與速度擾動(dòng)的相互作用產(chǎn)生高次散射波場；3)散射波場在背景速度模型中傳播.如果把常規(guī)全波形反演中的初始速度模型視為散射理論中散射波場傳播方程的背景速度模型，并對散射波場的傳播方程進(jìn)行線性化，則上述的散射波場的產(chǎn)生與傳播過程可為地震波場的全波形反演方法研究提供思路.

基于上述的認(rèn)識，本文根據(jù)散射理論中散射波場的傳播方程，首先推導(dǎo)出時(shí)間二階積分散射波場的傳播方程，然后把一階Born近似條件下的線性化迭代反演方法與散射波場的傳播過程相結(jié)合，構(gòu)建基于線性化迭代方式的時(shí)間二階積分散射波場的全波形反演方法.對波場的時(shí)間二階積分運(yùn)算相當(dāng)于一種低通濾波可以增強(qiáng)數(shù)據(jù)中已有的低頻成分，不同于Laplace、Laplace-Fourier域全波形反演方法、歸一化積分法和包絡(luò)反演方法中利用各種不同變換來提取地震數(shù)據(jù)中的低頻信息，利用時(shí)間二階積分增強(qiáng)地震數(shù)據(jù)中低頻成分的方法是來自于散射波場的傳播方程，它具有明確的數(shù)學(xué)物理含義.在時(shí)間二階積分散射波場全波形反演中，把常規(guī)全波形反演中的初始速度模型視為散射波場傳播中的背景速度模型，把觀測波場與由初始速度模型合成得到的計(jì)算波場之間的殘差波場視為速度擾動(dòng)產(chǎn)生的散射波場；在推導(dǎo)出的時(shí)間二階積分散射波場傳播方程的基礎(chǔ)上，利用逆?zhèn)鞑セ虬殡S傳播的方法由地表的時(shí)間二階積分殘差波場重建出地下散射源的分布，同時(shí)利用震源波場正向傳播的方法獲取地下入射波場分布，并根據(jù)時(shí)間二階積分散射波場線性傳播方程中散射波場與入射波場、速度擾動(dòng)間的線性關(guān)系，再利用類似偏移中的成像公式獲得初始速度模型的修正.最后再把時(shí)間二階積分散射波場全波形反演得到的結(jié)果作為常規(guī)全波形反演的初始速度模型，對地震波場進(jìn)行常規(guī)的全波形反演.通過Marmousi模型的全頻帶合成數(shù)據(jù)和缺低頻合成數(shù)據(jù)的數(shù)值試驗(yàn)驗(yàn)證本文所提出的全波形反演方法的正確性和有效性.

2　方法原理

眾所周知，正演是反演的基礎(chǔ)，因此我們在這節(jié)首先推導(dǎo)出時(shí)間二階積分散射波場的傳播方程，并分析時(shí)間二階積分波場的含義，然后再構(gòu)建基于時(shí)間二階積分散射波場傳播方程的全波形反演方法.

2.1時(shí)間二階積分散射波場的傳播方程

本文考慮單P波速度的全波形反演.給定初始(背景)速度模型v0(x)和震源函數(shù)s(t)，入射波場ui(x,t)的傳播方程有

(1-1)

(1-2)把速度模型表示為v(x)=v0(x)+δv(x)，δv(x)為速度擾動(dòng)，也稱為初始速度模型v0(x)的修正.根據(jù)散射理論，散射波場δu(x,t)的傳播方程有

(2)

方程(2)右邊的波場u(x,t)為與速度模型v(x)對應(yīng)的總波場；δu(x,t)也稱為與擾動(dòng)速度模型δv(x)有關(guān)的擾動(dòng)波場.下文把方程(2)中的近似等號寫為等號，把方程(2)變換到頻率域，有

(3)

(4)

則方程(3)可寫為

(5-1)

把方程(5-1)變換到時(shí)間域，有

(5-2)

根據(jù)表達(dá)式(4)(對該式的詳細(xì)解釋在下一小節(jié))可知，波場δui(x,t)是散射波場δu(x,t)的時(shí)間二階積分，我們稱之為時(shí)間二階積分散射波場，方程(5-1)、(5-2)分別為時(shí)間二階積分散射波場傳播方程的頻率域形式和時(shí)間域形式.

把一階Born近似條件應(yīng)用于方程(5-2)，得到

(6)

(7)式中，T0為波場的最大記錄時(shí)間；Ω為速度擾動(dòng)δv(x)的分布空間.令，

(8)

m(x,t)可視為產(chǎn)生時(shí)間二階積分散射波場的散射源，它是與δv(x)和ui(x,t)的分布、強(qiáng)弱有關(guān)的.利用式(8)，公式(7)可寫為

(9)

式(9)的積分是時(shí)間-空間域的褶積.如果已知左邊的δui(x,t)求右邊積分式中的m(x′,t′)，則意味求解一個(gè)第一類Fredholm線性積分方程.對式(9)兩邊做時(shí)間Fourier變換，有

(10)

(11)

2.2波場的時(shí)間二階積分

(12)

積分運(yùn)算相當(dāng)于一種低通濾波，因此對波場進(jìn)行式(12)所表示的頻率域運(yùn)算有助于增強(qiáng)波場的低頻成分.把式(12)進(jìn)行時(shí)間Fourier反變換，有

(13)

F-1表示時(shí)間Fourier反變換.因此時(shí)間二階積分波場ui(t)相對于原始波場u(t)其低頻成分得到了有效增強(qiáng).

假定波場u(t)為主頻等于15 Hz的Ricker子波，如圖1中的藍(lán)線所示，對應(yīng)的時(shí)間二階積分波場ui(t)如圖1中的紅線所示，它們的振幅譜如圖2所示，圖2中的藍(lán)線為Ricker子波的振幅譜，紅線為Ricker子波時(shí)間二階積分的振幅譜.由圖1的時(shí)間信號和圖2的振幅譜可看出，Ricker子波的二階積分信號相對于Ricker子波表現(xiàn)出明顯的低頻化，通過二階積分運(yùn)算使Ricker子波中的低頻成分得到了增強(qiáng)，二階積分信號的能量主要集中在低頻段，強(qiáng)化了低頻信息在數(shù)據(jù)中的作用.

圖1　主頻為15 Hz的Ricker子波及其時(shí)間二階積分：藍(lán)線為Ricker子波；紅線為Ricker子波的時(shí)間二階積分Fig.1　The ricker wavelet with peak-frequency 15Hz and its second-order time integral. The blue-line denotes ricker wavelet. The red-line denotes the second-order time integral of ricker wavelet

圖2　主頻為15 Hz的Ricker子波及其時(shí)間二階積分的振幅譜；藍(lán)線為Ricker子波的振幅譜，紅線為Ricker子波時(shí)間二階積分的振幅譜Fig.2　The amplitude spectra of ricker wavelet with peak-frequency 15 Hz and its second-order time integral. The blue-line denotes the amplitude spectra of ricker wavelet. The red-line denotes the amplitude spectra of second-order time integral of ricker wavelet

圖3為Marmousi模型的單炮合成記錄，圖4為對圖3中的地震波場進(jìn)行時(shí)間二階積分運(yùn)算得到的時(shí)間二階積分波場，對比圖3、4可看出，圖4中波場的波形變化相對于圖3中的波形表現(xiàn)緩慢，圖4中的波場和圖3中的波場相比表現(xiàn)出明顯的低頻化.圖5為圖3中波場的包絡(luò)信號圖，由于信號的包絡(luò)處理和所得到的包絡(luò)信號的數(shù)學(xué)物理意義不明確，圖5所示的包絡(luò)信號僅是一種能量集中在淺部、變化非常緩慢的非負(fù)值信號.圖4中的時(shí)間二階積分波場信號是受波場傳播方程控制，是一種具有明確數(shù)學(xué)物理意義的地震波場.

2.3時(shí)間二階積分波場的全波形反演

由于地震波場是速度模型的非線性函數(shù)，利用觀測的地震波場數(shù)據(jù)uobs(x,t)反演地下速度模型v(x)屬于非線性反演問題，本文采用線性化迭代反演方法進(jìn)行求解.

給定一個(gè)初始速度模型v0(x)，并計(jì)算出與之對應(yīng)的計(jì)算波場ucal(x,t)，把觀測波場uobs(x,t)與ucal(x,t)間的殘差波場，即

圖3　Marmousi模型的單炮合成記錄Fig.3　The synthetic single-shot record of Marmousi model

圖4　圖3所示的單炮記錄的時(shí)間二階積分Fig.4　The second-order time integral of the synthetic single-shot record in Fig.3

圖5　圖3所示的單炮記錄的包絡(luò)Fig.5　The envelope of the synthetic single-shot record in Fig.3

(14)

視為由地下真實(shí)速度模型與給定速度模型之間速度差異，δv(x)=v(x)-v0(x)，引起的散射波場.計(jì)算散射波場δu(x,t)的頻率域、時(shí)間域二階積分散射波場，即

(15)

(16)

式(15)中的F表示時(shí)間Fourier變換.

(17)

或

(18)

(19)

(20)

得到了速度擾動(dòng)δv(x)的解估計(jì)后，我們就可以對給定的初始速度模型v0(x)進(jìn)行修正，得到新的速度模型v1(x)，即v1(x)=v0(x)+αδv(x)，其中α為修正步長，然后再把v1(x)作為下一步反演的初始速度模型，進(jìn)行迭代反演.

(21)

利用最小二乘方法得到算子方程(21)的最小二乘解，即

(22)

式中G-1為散射波場傳播算子G的最小二乘逆算子，有

(23)

式(22)所表示的運(yùn)算為在時(shí)間-空間域?qū)Σ▓靓膗i進(jìn)行逆?zhèn)鞑?由于這種逆?zhèn)鞑ゲ粌H計(jì)算量巨大，而且還很不穩(wěn)定，常用G的伴隨算子Gt代替逆算子，即G-1≈Gt，因此式(22)改寫為

(24)

式(24)所表示的運(yùn)算為在時(shí)間-空間域?qū)Σ▓靓膗i進(jìn)行伴隨(逆時(shí))傳播.這種伴隨傳播不僅計(jì)算量較小，而且還穩(wěn)定.

把式(23)代入式(22)，再把式(22)代入式(17)就可得到速度擾動(dòng)的解估計(jì)，并把各個(gè)共炮道集數(shù)據(jù)的反演結(jié)果疊加，有

(25)

也可以把式(23)代入式(22)，再把式(22)代入式(19)得到速度擾動(dòng)的解估計(jì)，并把各個(gè)共炮道集數(shù)據(jù)的反演結(jié)果疊加，有

(25-1)

把式(24)代入式(19)，同樣也把各個(gè)共炮道集數(shù)據(jù)的反演結(jié)果疊加，有

(26)

在上述式中，xs表示共炮道集的炮點(diǎn)位置；ui(xs,x,t)表示初始速度模型下正向傳播的震源波場.式(25)、(25-1)、(26)為在時(shí)間域求解δv(x)的公式.

上述的時(shí)間域反演工作也同樣可在頻率域?qū)崿F(xiàn)，把頻率域的方程(10)寫為算子形式，有

(27)

同樣利用最小二乘方法得到算子方程(27)的最小二乘解，即

(28)

(29)

(30)

把式(29)代入式(28)，再把式(28)代入式(18)就可得到速度擾動(dòng)的解估計(jì)，并把各個(gè)共炮道集數(shù)據(jù)的反演結(jié)果疊加，有

(31)

把式(29)代入式(28)，再把式(28)代入式(20)得到速度擾動(dòng)的解估計(jì)，并把各個(gè)共炮道集數(shù)據(jù)的反演結(jié)果疊加，有

(31-1)把式(30)代入式(20)，同樣也把各個(gè)共炮道集數(shù)據(jù)的反演結(jié)果疊加，有

(32)

把正向傳播的震源波場和二階積分波場逆?zhèn)鞑サ木唧w形式代入式(31)、(31-1)，分別有，

(33)

(33-1)

把正向傳播的震源波場和二階積分波場逆時(shí)傳播的具體形式代入式(32)，有，

(34)

利用上述的時(shí)間二階積分波場的全波形反演雖然可以有效地減小對初始速度模型的依賴性，但時(shí)間二階積分波場相對于原始波場存在明顯的低頻化，因此時(shí)間二階積分波場全波形反演得到的反演結(jié)果分辨率受到一定的限制，為了充分地利用地震數(shù)據(jù)中的全部信息，以得到高分辨率的全波形反演結(jié)果，我們把時(shí)間二階積分波場全波形反演得到的反演結(jié)果作為常規(guī)全波形反演的初始速度模型，再進(jìn)行常規(guī)的全波形反演.

圖6　Marmousi模型的速度分布Fig.6　The velocity distribution of Marmousi model

3　數(shù)值試驗(yàn)

為了驗(yàn)證本文所提出的時(shí)間二階積分波場全波形反演(英文簡寫為IntI)的正確性和有效性，我們用Marmousi模型的合成地震數(shù)據(jù)進(jìn)行試驗(yàn).試驗(yàn)中Marmousi模型的網(wǎng)格化參數(shù)：nx=369，nz=125，dx=25m，dz=24m.圖6為Marmousi模型的速度分布圖.試驗(yàn)用的合成地震數(shù)據(jù)共有92炮，炮間距100m，每炮369道，道間距25m.時(shí)間采樣長度5s，采樣間隔2ms.震源為主頻8Hz的Ricker子波.我們合成了2套地震數(shù)據(jù)用于試驗(yàn)，一套是全頻帶地震數(shù)據(jù)；另一套是完全缺失4Hz以下頻率成分，并以4～5Hz為過渡帶的缺低頻地震數(shù)據(jù).圖7為試驗(yàn)用的常梯度初始速度模型.

3.1全頻帶地震數(shù)據(jù)試驗(yàn)

對于合成得到的全頻帶地震數(shù)據(jù)全波形反演試驗(yàn)，我們首先用圖7所示的常梯度速度模型作為時(shí)間二階積分波場全波形反演的初始速度模型，在反演中應(yīng)用時(shí)間域公式(26)進(jìn)行速度修正量的計(jì)算.圖8為應(yīng)用時(shí)間二階積分波場全波形反演得到的反演結(jié)果，對比圖6、8可看出，由時(shí)間二階積分波場反演得到的速度模型大體結(jié)構(gòu)正確，但分辨率較低.我們再把圖8所示的速度模型作為常規(guī)全波形反演(英文簡寫為FWI)的初始速度模型，圖9為用圖8作初始速度模型的常規(guī)全波形反演結(jié)果.作為比較，我們還用圖7所示的常梯度速度模型作為常規(guī)全波形反演的初始速度模型，圖10為用圖7作初始速度模型的常規(guī)全波形反演結(jié)果.對比圖6、9、10可看出，對于同樣的常梯度初始速度模型，時(shí)間二階積分波場全波形反演加常規(guī)全波形反演(IntI+FWI)的反演結(jié)果要明顯優(yōu)于單純的常規(guī)全波形反演的反演結(jié)果.為了更仔細(xì)對比上述兩種反演結(jié)果的細(xì)節(jié)，圖11為Marmousi模型中第85、235網(wǎng)格點(diǎn)兩種反演結(jié)果的速度-深度曲線對比圖，圖中綠線為真實(shí)速度，黑線為初始速度，紅線為常規(guī)全波形反演的結(jié)果，藍(lán)線為時(shí)間二階積分波場全波形反演加常規(guī)全波形反演的結(jié)果.由圖11中的速度-深度曲線可看出，時(shí)間二階積分波場全波形反演加常規(guī)全波形反演能很好地重建出真實(shí)速度模型.

圖7　用于反演的常梯度初始速度模型Fig.7　The initial velocity model with constant gradient for inversion

圖8　時(shí)間二階積分波場全波形反演(IntI)的反演結(jié)果Fig.8　The inversion result by full waveform inversion of the second-order time integral wavefield

圖9　時(shí)間二階積分波場全波形反演加常規(guī)全波形反演(IntI+FWI)的反演結(jié)果Fig.9　The inversion result by full waveform inversion of the second-order time integral wavefield plus conventional full waveform inversion

圖10　常規(guī)全波形反演(FWI)的反演結(jié)果Fig.10　The inversion result by conventional full waveform inversion

圖11　第85、235網(wǎng)格點(diǎn)兩種反演結(jié)果的速度-深度曲線對比：上圖為第85網(wǎng)格點(diǎn)；下圖為235網(wǎng)格點(diǎn).圖中綠線為真實(shí)速度， 黑線為初始速度， 紅線為常規(guī)全波形反演的結(jié)果， 藍(lán)線為時(shí)間二階積分波場全波形反演加常規(guī)全波形反演的結(jié)果Fig.11　Comparisons of velocity-depth curves of the two inversion results at 85th grid-point and 235th grid-point: the up-panel is for 85th grid-point; the bottom-panel is for 235th grid-point. The green line denotes the true velocity, the black line denotes the initial velocity. The red line denotes the velocity from conventional FWI. The blue line denotes the velocity from full waveform inversion of the second-order time integral wavefield plus conventional FWI

3.2缺低頻地震數(shù)據(jù)試驗(yàn)

為了進(jìn)一步驗(yàn)證本文所提出的全波形反演方法對于缺低頻地震數(shù)據(jù)的反演效果，我們進(jìn)行了對完全缺失4 Hz以下頻率成分，并以4～5 Hz為過渡帶的缺低頻地震數(shù)據(jù)全波形反演試驗(yàn).在試驗(yàn)中，我們首先用圖7所示的常梯度速度模型作為時(shí)間二階積分波場全波形反演的初始速度模型，在反演中仍然應(yīng)用時(shí)間域公式(26)進(jìn)行速度修正量的計(jì)算.圖12為應(yīng)用時(shí)間二階積分波場全波形反演得到的反演結(jié)果.我們再把圖12所示的速度模型作為常規(guī)全波形反演的初始速度模型，圖13為用圖12作初始速度模型的常規(guī)全波形反演結(jié)果.作為比較，我們還用圖7所示的常梯度速度模型作為初始速度模型對缺低頻地震數(shù)據(jù)進(jìn)行了單純的常規(guī)全波形反演和包絡(luò)反演加常規(guī)全波形反演(EI+FWI).圖14為用圖7作初始速度模型的常規(guī)全波形反演結(jié)果，圖15為用圖7作初始速度模型的包絡(luò)反演加常規(guī)全波形反演的反演結(jié)果.

對比圖6、13、14、15可看出，對于同樣的常梯度初始速度模型，時(shí)間二階積分波場全波形反演加常規(guī)全波形反演的反演結(jié)果要明顯優(yōu)于單純的常規(guī)全波形反演的反演結(jié)果和包絡(luò)反演加常規(guī)全波形反演的反演結(jié)果，包絡(luò)反演加常規(guī)全波形反演的反演結(jié)果次之，單純的常規(guī)全波形反演的反演結(jié)果最差.為了更仔細(xì)地對比上述三種反演結(jié)果的細(xì)節(jié)，圖16為Marmousi模型中第85、235網(wǎng)格點(diǎn)三種反演結(jié)果的速度-深度曲線對比圖，圖中綠線為真實(shí)速度，黑線為初始速度，深紅線為常規(guī)全波形反演的結(jié)果，淺紅線為包絡(luò)反演加常規(guī)全波形反演的結(jié)果，藍(lán)線為時(shí)間二階積分波場全波形反演加常規(guī)全波形反演的結(jié)果.由圖16中的速度-深度曲線可看出，即使對于缺低頻的地震數(shù)據(jù)，時(shí)間二階積分波場全波形反演加常規(guī)全波形反演也能很好地重建出真實(shí)速度模型.

圖12　時(shí)間二階積分波場全波形反演(IntI)的反演結(jié)果Fig.12　The inversion result by full waveform inversion of the second-order time integral wavefield

圖13　時(shí)間二階積分波場全波形反演加常規(guī)全波形反演(IntI+FWI)的反演結(jié)果Fig.13　The inversion result by full waveform inversion of the second-order time integral wavefield plus conventional full waveform inversion

由全頻帶地震數(shù)據(jù)的時(shí)間二階積分波場全波形反演得到的結(jié)果圖8與缺低頻地震數(shù)據(jù)的時(shí)間二階積分波場全波形反演得到的結(jié)果圖12的對比，以及全頻帶地震數(shù)據(jù)的時(shí)間二階積分波場全波形反演加常規(guī)全波形反演得到的結(jié)果圖9與缺低頻地震數(shù)據(jù)的時(shí)間二階積分波場全波形反演加常規(guī)全波形反演得到的結(jié)果圖13的對比，我們可看出對于缺低頻地震數(shù)據(jù)應(yīng)用時(shí)間二階積分波場全波形反演方法進(jìn)行反演可得到與全頻帶地震數(shù)據(jù)應(yīng)用時(shí)間二階積分波場全波形反演方法進(jìn)行反演基本一致的反演結(jié)果，這表明時(shí)間二階積分波場全波形反演方法對缺失低頻地震數(shù)據(jù)的有效性.圖17為反演結(jié)果圖9、13中第85、235網(wǎng)格點(diǎn)的速度-深度曲線對比圖，圖中綠線為真實(shí)速度，黑線為初始速度，紅線為缺低頻地震數(shù)據(jù)的反演結(jié)果，藍(lán)線為全頻帶地震數(shù)據(jù)的反演結(jié)果，由圖17的對比可看出缺失4 Hz以下頻譜成分?jǐn)?shù)據(jù)的反演結(jié)果與全頻帶數(shù)據(jù)的反演結(jié)果沒有明顯差異.

圖14　常規(guī)全波形反演(FWI)的反演結(jié)果Fig.14　The inversion result by conventional full waveform inversion

圖15　包絡(luò)反演加常規(guī)全波形反演(EI+FWI)的反演結(jié)果Fig.15　The inversion result by envelope inversion plus conventional full waveform inversion

圖16　第85、235網(wǎng)格點(diǎn)三種反演結(jié)果的速度-深度曲線對比：上圖為第85網(wǎng)格點(diǎn)；下圖為235網(wǎng)格點(diǎn).圖中綠線為真實(shí)速度，黑線為初始速度，深紅線為常規(guī)全波形反演的結(jié)果，淺紅線為包絡(luò)反演加常規(guī)全波形反演的結(jié)果，藍(lán)線為時(shí)間二階積分波場全波形反演加常規(guī)全波形反演的結(jié)果Fig.16　Comparisons of velocity-depth curves of the three inversion results at 85th grid-point and 235th grid-point: the up-panel is for 85th grid-point; the bottom-panel is for 235th grid-point. The green line denotes the true velocity. The black line denotes the initial velocity, the dark-red line denotes the velocity from conventional FWI. The light-red line denotes the velocity from envelope inversion plus conventional FWI. The blue line denotes the velocity from full waveform inversion of the second-order time integral wavefield plus conventional FWI

圖17　第85、235網(wǎng)格點(diǎn)全頻帶地震數(shù)據(jù)與缺低頻地震數(shù)據(jù)反演結(jié)果的速度-深度曲線對比：上圖為第85網(wǎng)格點(diǎn)；下圖為235網(wǎng)格點(diǎn).圖中綠線為真實(shí)速度，黑線為初始速度，紅線為缺低頻地震數(shù)據(jù)的反演結(jié)果，藍(lán)線為全頻帶地震數(shù)據(jù)的反演結(jié)果Fig.17　Comparisons of velocity-depth curves of inversion result by full band seismic data and inversion result by lack low frequencies seismic data at 85th grid-point and 235th grid-point. The up-panel is for 85th grid-point; the bottom-panel is for 235th grid-point. The green line denotes the true velocity, the black line denotes the initial velocity, the red line denotes the velocity frominversion result by lack low frequencies seismic. The blue line denotes the velocity from inversion result by full band seismic data

同時(shí)給出時(shí)間域和頻率域的公式推導(dǎo)，主要是便于對本文方法的理解和與現(xiàn)有一些FWI公式的對比.本文反演計(jì)算是在時(shí)間-空間域進(jìn)行的，由給出的公式可知同樣可以在頻率-空間域進(jìn)行反演計(jì)算.

4　結(jié)論

基于局部極值問題和周期跳現(xiàn)象在低頻地震數(shù)據(jù)全波形反演中相對不敏感的認(rèn)識，由散射波場傳播方程推導(dǎo)出了時(shí)間二階積分散射波場的傳播方程，以及一階Born近似下的時(shí)間二階積分散射波場的線性傳播方程.把推導(dǎo)出的時(shí)間二階積分散射波場線性傳播方程應(yīng)用于非線性反演問題的線性化迭代反演，并結(jié)合時(shí)間二階積分散射波場線性傳播方程中散射波場與入射波場和速度擾動(dòng)間的線性關(guān)系，建立了一種首先反演地下散射源分布，然后再求解速度擾動(dòng)的時(shí)間二階積分波場全波形反演方法.對地震波場的時(shí)間二階積分可有效地增強(qiáng)波場中已有的低頻成分，強(qiáng)化了地震數(shù)據(jù)中的低頻信息在反演中的作用，這是本文方法與由低頻到高頻的分頻反演方法的不同之處，也是其長處，使本文提出的時(shí)間二階積分波場全波形反演方法對初始模型的依賴性減弱.不同于包絡(luò)全波形反演和歸一化積分全波形反演等方法中的數(shù)據(jù)變換方法，本文的波場時(shí)間二階積分運(yùn)算是來自于描述散射波場傳播的數(shù)學(xué)物理方程，因此得到的時(shí)間二階積分波場具有清楚的數(shù)學(xué)物理含義.把時(shí)間二階積分波場全波形反演方法與常規(guī)全波形反演方法相結(jié)合，不僅可以得到高分率的反演結(jié)果，還能在缺低頻地震數(shù)據(jù)的全波形反演中發(fā)揮作用.本文的數(shù)值試驗(yàn)顯示時(shí)間二階積分波場全波形反演方法應(yīng)用于缺失4 Hz以下頻譜成分?jǐn)?shù)據(jù)和全頻帶數(shù)據(jù)得到的兩個(gè)反演結(jié)果沒有明顯的差異.

致謝感謝審稿專家和編輯部的支持.

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

Chauris H, Donno D, Calandra H. 2012. Velocity estimation with the normalized integration method.∥ 74thEAGE Conference and Exhibition incorporating EUROPEC 2012. SPE, EAGE, W020.

Chew W C. 1990. Waves and Fields in Inhomogeneous Media. New York: Dover Publications.

Chi B, Dong L, Liu Y. 2013. Full waveform inversion based on envelope objective function.∥ 75thEAGE Conference and Exhibition incorporating SPE EUROPEC 2013. EAGE, TU-P04-09.

Chi B X, Dong L G, Liu Y Z. 2014. Full waveform inversion method using envelope objective function without low frequency data.JournalofAppliedGeophysics, 109: 36-46.Donno D, Chauris H, Calandra H. 2013. Estimating the background velocity model with the normalized integration method.∥ 75thEAGE Conference and Exhibition incorporating SPE EUROPEC 2013. EAGE, Tu-07-04.

Lailly P. 1983. The seismic inverse problem as a sequence of before stack migrations.∥ Conference on Inverse Scattering, Theory and Applications. Philadelphia, PA: Society of Industrial and Applied Mathematics, 206-220. Pratt R G, Shin C, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362. Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part I: Theory and verification in a physical scale model.Geophysics, 64(3): 888-901.

Shin C, Min D J. 2006. Waveform inversion using a logarithmic wavefield.Geophysics, 71(3): R31-R42.

Shin C, Pyun S, Bednar J B. 2007. Comparison of waveform inversion, Part 1: Conventional wavefield vs. logarithmic wavefield.GeophysicalProspecting, 55(4): 449-464.

Shin C, Cha Y H. 2008. Waveform inversion in the Laplace domain.GeophysicalJournalInternational, 173(3): 922-931.

Shin C, Cha Y H. 2009. Waveform inversion in the Laplace-Fourier domain.GeophysicalJournalInternational, 177(3): 1067-1079.

Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies.Geophysics, 69(1): 231-248.

Sirgue L. 2006. The importance of low frequency and large offset in waveform inversion.∥ 68thEAGE Conference & Technical Exhibition. SPE, EAGE, A037.

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

Tarantola A. 1986. A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics, 51(10): 1893-1903.

Tarantola A. 1987. Inverse problem theory: Methods for data fitting and model parameter estimation. Amsterdam: Elsevier Science Publ. Co., Inc.

Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6): WCC1-WCC26.

Vogel C R. 2002. Computational methods for inverse problems. Philadelphia, PA, US: Society for Industrial and Applied Mathematics.

Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model.Geophysics, 79(3): WA13-WA24.

(本文編輯劉少華)

Full waveform inversion of the second-order time integral wavefield

CHEN Sheng-Chang， CHEN Guo-Xin

SchoolofEarthSciences，ZhejiangUniversity，Hangzhou310027,China

We proposed a new full waveform inversion (FWI) method, namely full waveform inversion of the second-order time integral wavefield, by enhancing low-frequency components of seismic data with a second-order time integration to seismic wavefield, which can efficiently reduce the initial model dependence of FWI. According to the propagation equation of scattering wavefield in scattering theory, we derived a propagation equation for the scattering wavefield with second-order time integral, and used the leading order Born approximation for the linearization to the propagation equation. Based on this equation, using the scattering wavefield to inverse the distribution of scattering sources in subsurface, and using wavefield modeling to construct the incident wavefield, and according to the linear relationship between the scattering wavefield and incident wavefield, velocity perturbation in the linear propagation equation for the scattering wavefield with second-order time integral, we applied a formula similar to the imaging formula of migration to estimate velocity perturbation, and established an iterative inversion method for the full waveform inversion of second-order integral wavefield. Applying the inversion result of full waveform inversion of second-order integral wavefield as the initial velocity model for the conventional FWI can efficiently reduce the initial model dependence of FWI. Numerical tests using synthetic data of the Marmousi model demonstrate the validity and feasibility of the proposed method. The final results of the new method can obtain much improved results than the conventional FWI. Furthermore, to test the independence to the seismic frequency-band, we used a low-cut source wavelet (cutting from 4Hz below) to generate the synthetic data. The inversion results by our new method show no appreciable difference from the full-band source results.

Propagation equation of scattering wave; Second-order time integral; Enhancement low-frequency components; Born approximation; Full waveform inversion

10.6038/cjg20161021.

國家自然科學(xué)基金項(xiàng)目(41074133，41374001)資助．

陳生昌，男,1965年生，教授，博士生導(dǎo)師．主要從事勘探地球物理和計(jì)算地球物理研究. E-mail: chenshengc@zju.edu.cn

10.6038/cjg20161021

P631

2015-11-10, 2016-03-18收修定稿

陳生昌， 陳國新. 2016. 時(shí)間二階積分波場的全波形反演. 地球物理學(xué)報(bào)，59(10):3765-3776,

Chen S C， Chen G X. 2016. Full waveform inversion of the second-order time integral wavefield.ChineseJ.Geophys. (in Chinese),59(10):3765-3776,doi:10.6038/cjg20161021.