胡江濤, 王華忠

同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與反演成像研究組， 上海 200092

?

基于解析時(shí)間波場(chǎng)外推與波場(chǎng)分解的逆時(shí)偏移方法研究

胡江濤, 王華忠

同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與反演成像研究組， 上海 200092

雙程波方程逆時(shí)深度偏移是復(fù)雜介質(zhì)高精度成像的有效技術(shù)，但其結(jié)果中通常包含成像方法引起的噪音和假象，一般的濾波方法會(huì)破壞成像剖面上的振幅，其中的假象也會(huì)給后續(xù)地質(zhì)解釋帶來(lái)困擾.將波場(chǎng)進(jìn)行方向分解然后實(shí)現(xiàn)入射波與反射波的相關(guān)成像能夠有效地消除這類成像噪音，并提高逆時(shí)偏移成像質(zhì)量.波傳播方向的分解通常在頻率波數(shù)域?qū)崿F(xiàn)，它會(huì)占用大量的存儲(chǔ)和計(jì)算資源，不便于在沿時(shí)間外推的逆時(shí)深度偏移中應(yīng)用.本文提出解析時(shí)間波場(chǎng)外推方法，可以在時(shí)間外推的每個(gè)時(shí)間片上實(shí)現(xiàn)波傳播方向的顯式分解，逆時(shí)深度偏移中利用分解后的炮檢波場(chǎng)進(jìn)行對(duì)應(yīng)的相關(guān)運(yùn)算，實(shí)現(xiàn)成像噪音和成像信號(hào)的分離.在模型和實(shí)際數(shù)據(jù)上的測(cè)試表明，相比于常規(guī)互相關(guān)逆時(shí)偏移成像結(jié)果，本文方法能夠有效地消除低頻成像噪音和特殊地質(zhì)構(gòu)造導(dǎo)致的成像假象.

解析時(shí)間波場(chǎng)外推； 波場(chǎng)分解； 噪音消除； 逆時(shí)深度偏移

1 引言

基于雙向波波動(dòng)方程的逆時(shí)偏移方法(Baysal et al.,1983; McMechan, 1983; Whitmore, 1983)克服了Kirchhoff偏移的焦散和單走時(shí)問(wèn)題，單向波波動(dòng)方程的陡傾角成像限制問(wèn)題，成為復(fù)雜介質(zhì)高精度成像的有力工具.借助GPU和多核CPU等高性能計(jì)算設(shè)備的發(fā)展(劉紅偉等，2010； 李博等， 2010； 劉守偉等， 2013)，逆時(shí)偏移方法已經(jīng)被廣泛應(yīng)用于工業(yè)生產(chǎn)中.但其成像結(jié)果通常存在一定的、由成像方法引起的成像噪音和假象，這是由于不對(duì)應(yīng)的波場(chǎng)進(jìn)行相關(guān)成像而引起的.

低頻成像噪音是逆時(shí)偏移成像中最常見的問(wèn)題.在速度模型中存在高波數(shù)的速度擾動(dòng)或反射界面時(shí)，正、反傳波場(chǎng)中存在較強(qiáng)的背向散射(或反射)，低頻成像噪音非常嚴(yán)重.Loewenthal 等(1987)采用平滑速度場(chǎng)的方式來(lái)消除背向散射波場(chǎng)進(jìn)而消除低頻噪音.Fletcher 等(2006)采用壓制反射界面背向散射的方式來(lái)消除低頻噪音.Yoon et和Marfurt (2006)、康瑋和程玖兵(2012)采用波印廷矢量方法進(jìn)行波場(chǎng)傳播方向預(yù)測(cè)，實(shí)現(xiàn)成像結(jié)果的角度切除，進(jìn)而消除低頻噪音.Guitton 等(2007)采用最小平方濾波的方式實(shí)現(xiàn)低頻噪音消除.Zhang和Sun(2009)采用Laplace濾波方法來(lái)進(jìn)行低頻噪音消除，并通過(guò)頻率波數(shù)域關(guān)系實(shí)現(xiàn)了Laplace濾波的振幅校正.Liu 等 (2007， 2011)提出采用波場(chǎng)分解方式來(lái)構(gòu)造逆時(shí)偏移成像條件，實(shí)現(xiàn)入射波和反射波的相關(guān)成像，進(jìn)而達(dá)到消除低頻噪音的目的.杜啟振等(2013)對(duì)上述低頻噪音壓制方法做了詳細(xì)分類并對(duì)比了各種方法的優(yōu)缺點(diǎn)，指出基于波場(chǎng)分解的成像方式在復(fù)雜介質(zhì)中能更好地消除低頻噪音和實(shí)現(xiàn)振幅保真.

除了低頻成像噪音，F(xiàn)ei 等(2010，2014)分析了逆時(shí)偏移對(duì)偏移速度的敏感性，并指出在速度模型精度較差時(shí)，當(dāng)向斜構(gòu)造存在強(qiáng)速度梯度時(shí)，逆時(shí)偏移會(huì)產(chǎn)生較強(qiáng)的假象.產(chǎn)生這類假象的構(gòu)造通常會(huì)出現(xiàn)在巖丘邊界，潛山頂面及火成巖與圍巖的接觸面等地方.他指出在這類構(gòu)造區(qū)域，采用波場(chǎng)分解成像方式能夠有效地消除成像假象.雖然波場(chǎng)分解成像能夠有效地消除成像噪音，但它通常需要在頻率波數(shù)域進(jìn)行實(shí)現(xiàn).在逆時(shí)深度偏移過(guò)程中，波場(chǎng)分解通常需要將炮點(diǎn)和檢波點(diǎn)波場(chǎng)存盤并進(jìn)行關(guān)于時(shí)間的傅里葉變換，而在沿時(shí)間的外推過(guò)程中，時(shí)間維處于數(shù)據(jù)的最慢維，因此一般的波場(chǎng)分離成像條件的I/O量和計(jì)算量非常大.

針對(duì)上述問(wèn)題，本文提出了解析時(shí)間波場(chǎng)外推及波場(chǎng)分解方法.在逆時(shí)深度偏移過(guò)程中，僅利用外推過(guò)程中的、每個(gè)時(shí)間步的解析時(shí)間波場(chǎng)進(jìn)行不同傳播方向的波場(chǎng)的顯式分離，避免了時(shí)間域向頻率域轉(zhuǎn)換的I/O量和計(jì)算量.方向分解后的波場(chǎng)進(jìn)行對(duì)應(yīng)的波場(chǎng)相關(guān)可以有效地壓制成像噪音和假象.模型和實(shí)際數(shù)據(jù)的測(cè)試表明了所提出方法的有效性.

2 逆時(shí)深度偏移及波場(chǎng)分解成像條件

(2)其中，dR為檢波點(diǎn)端接收的炮波場(chǎng)，gx=(gx,gy)為檢波點(diǎn)空間位置，uR為延拓后的檢波點(diǎn)波場(chǎng).

雙向波沿時(shí)間的正向和反向外推可以模擬聲波方程(1)和(2)所能模擬的所有波現(xiàn)象，在偏移速度場(chǎng)比較精確的情況下，理論上可以對(duì)一次和多次散射(和反射)波進(jìn)行偏移成像.這是逆時(shí)深度偏移成為復(fù)雜介質(zhì)成像首選方法技術(shù)的原因.存在的問(wèn)題是計(jì)算代價(jià)大和不對(duì)應(yīng)的波場(chǎng)相關(guān)產(chǎn)生的噪音和假象.一般地，后者需要更改成像條件來(lái)解決.

本文采用有限差分方法求解公式(1)和公式(2)來(lái)實(shí)現(xiàn)波場(chǎng)的延拓.將延拓的波場(chǎng)應(yīng)用如下的零延遲相關(guān)成像條件即可提取地下界面的成像結(jié)果，

(3)

其中，tmax為延拓的最大時(shí)間值.

在速度模型比較準(zhǔn)時(shí)，炮點(diǎn)正向外推和檢波點(diǎn)反向外推波場(chǎng)都會(huì)產(chǎn)生前向和背向散射(和反射)，炮點(diǎn)和檢波點(diǎn)波場(chǎng)中同方向傳播的波場(chǎng)相關(guān)產(chǎn)生低頻噪音，某些反方向傳播的波場(chǎng)會(huì)形成成像假象(Liuetal. 2011;Feietal.，2010).因此可以利用波場(chǎng)的方向分解來(lái)實(shí)現(xiàn)相對(duì)應(yīng)波場(chǎng)的相關(guān)，進(jìn)而壓制低頻噪音和/或成像假象.

Liu等(2007)提出將炮點(diǎn)波場(chǎng)和檢波點(diǎn)波場(chǎng)進(jìn)行上下行波分解：

(4a)

(4b)

其中，上標(biāo)up代表上行波，down代表下行波.需要說(shuō)明的是，Liu等(2007)所指炮檢端的上下行波是從炮點(diǎn)端來(lái)看的，其上下行波的定義與單向波方程的上下行波定義一致.

將公式(4)帶入公式(3)，可得出相關(guān)成像的結(jié)果為

=I1(x)+I2(x)+I3(x)+I4(x),

(5)

其中，I1(x)為炮點(diǎn)下行波場(chǎng)與檢波點(diǎn)上行波場(chǎng)相關(guān)成像結(jié)果，I2(x)為炮點(diǎn)上行波場(chǎng)與檢波點(diǎn)下行波場(chǎng)相關(guān)成像結(jié)果，I3(x)為炮點(diǎn)下行波場(chǎng)與檢波點(diǎn)下行波場(chǎng)相關(guān)成像結(jié)果，I4(x)為炮點(diǎn)上行波場(chǎng)與檢波點(diǎn)上行波場(chǎng)相關(guān)成像結(jié)果.在速度模型精度較高時(shí)，I1(x)和I2(x)分別是炮檢波場(chǎng)中反向傳播的波場(chǎng)的相關(guān)，對(duì)應(yīng)為反射界面(或散射點(diǎn))的像.Fei等(2014)指出在速度模型精度不是太高時(shí)，I2(x)在特殊地質(zhì)構(gòu)造區(qū)域可能導(dǎo)致成像假象.I3(x)和I4(x)分別是炮檢波場(chǎng)中同向傳播波場(chǎng)的相關(guān)，表現(xiàn)為低頻噪音.由公式(5)也可看出，采用波場(chǎng)分解成像條件可以有效地實(shí)現(xiàn)逆時(shí)偏移成像噪音和成像信號(hào)的分離.

3 解析時(shí)間波場(chǎng)外推及波場(chǎng)分解

雖然波場(chǎng)分離成像條件能夠有效地分離逆時(shí)偏移中的信號(hào)和噪音，但是波場(chǎng)分解方法是其實(shí)現(xiàn)的核心.下面以炮點(diǎn)波場(chǎng)為例來(lái)說(shuō)明波場(chǎng)分解過(guò)程，檢波點(diǎn)波場(chǎng)的方向分解與其類似，不再贅述.波場(chǎng)方向分解通常借助傅里葉變換來(lái)實(shí)現(xiàn)，波場(chǎng)的時(shí)間和空間傅里葉變換可以表示為

×e-i(ω t+xkx+yky+zkz)dtdxdydz.

(6a)

波場(chǎng)的上下行波的定義為

(7a)

(7b)

同理可以定義左右行波和前后行波.由公式(7)可以看出波場(chǎng)的方向由頻率和波數(shù)的符號(hào)共同決定，也就是說(shuō)常規(guī)波場(chǎng)分離需要在頻率域進(jìn)行實(shí)現(xiàn).時(shí)間域逆時(shí)偏移實(shí)現(xiàn)常規(guī)波場(chǎng)分離，需要將地下波場(chǎng)存盤，并進(jìn)行關(guān)于最慢維的時(shí)間傅里葉變換，所以其I/O量和計(jì)算量很大.但是時(shí)間域外推中比較便于在波場(chǎng)時(shí)間片上獲取空間波數(shù)信息，所以為了提高計(jì)算效率，Liu等(2011)根據(jù)頻率和波數(shù)的符號(hào)關(guān)系推導(dǎo)出

(8)

其中，

(9a)

(9b)

(9c)

(9d)

由公式(7)和公式(8)可以看出，波場(chǎng)的顯式分離需要由頻率和波數(shù)的符號(hào)共同決定，而為了實(shí)現(xiàn)的高效，需要盡量利用空間波數(shù)的符號(hào)來(lái)進(jìn)行波場(chǎng)分離.那么如果能夠獲取一個(gè)時(shí)間波場(chǎng)，其頻譜僅包含正或負(fù)頻率，那么波場(chǎng)的方向分解就可以依靠空間波數(shù)的符號(hào)來(lái)實(shí)現(xiàn).這樣的時(shí)間信號(hào)在信號(hào)處理中廣泛應(yīng)用，這就是復(fù)信號(hào)，也即解析時(shí)間信號(hào).解析時(shí)間信號(hào)的實(shí)部為信號(hào)本身，而虛部為信號(hào)的Hilbert變換.將解析時(shí)間信號(hào)概念擴(kuò)展到波場(chǎng)上，我們將其稱為解析時(shí)間波場(chǎng).通過(guò)常規(guī)波場(chǎng)傳播方程可以獲得波場(chǎng)本身(即解析時(shí)間波場(chǎng)的實(shí)部)，但是無(wú)法獲得波場(chǎng)的Hilbert變換(即解析時(shí)間波場(chǎng)的虛部).可以對(duì)地下波場(chǎng)進(jìn)行Hilbert變換來(lái)實(shí)現(xiàn)解析時(shí)間波場(chǎng)虛部的構(gòu)造，但Hilbert變換為一個(gè)時(shí)間方向的卷積，這就需要將波場(chǎng)存盤并進(jìn)行卷積運(yùn)算，I/O量同樣會(huì)增加.Hilbert變換為一個(gè)90°相移運(yùn)算，在頻率域是一個(gè)相位濾波器.而波場(chǎng)傳播算子為線性算子，也即傳播方程(1)和(2)的源項(xiàng)和波場(chǎng)成線性關(guān)系(Liuetal., 2006).那么參考波動(dòng)方程相位編碼的實(shí)現(xiàn)(Liuetal., 2006)，可以對(duì)波動(dòng)方程的源項(xiàng)進(jìn)行Hilbert變換.由于波場(chǎng)傳播算子的線性特征，波場(chǎng)也進(jìn)行了90°相位移動(dòng)，這樣就獲得了解析時(shí)間波場(chǎng)的虛部.炮點(diǎn)的解析時(shí)間波場(chǎng)外推可以表示為

(10)

(11)

(12a)

(12c)

(12d)其中，Re代表取實(shí)運(yùn)算.同樣，可以寫出檢波點(diǎn)的上下行波：

(13a)

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

為了測(cè)試方法的有效性，下面采用向斜模型進(jìn)行測(cè)試，如圖1所示.該模型第一層，第二層和第三層速度值分別為2000m·s-1, 2200m·s-1和4500m·s-1.第二層和第三層之間在實(shí)際情況可能為巖丘邊界，火成巖邊界及潛山頂面等地質(zhì)構(gòu)造.模擬數(shù)據(jù)一共100炮，炮間隔為50m，最大偏移距為4000m.采用常規(guī)相關(guān)逆時(shí)偏移成像方法得到的結(jié)果如圖2所示，由圖可以看出，成像結(jié)果受到嚴(yán)重的低頻噪音影響.

圖1 準(zhǔn)確速度模型Fig.1 Exact velocity model

圖2 傳統(tǒng)逆時(shí)偏移成像結(jié)果Fig.2 Result of conventional RTM

下面利用例子來(lái)驗(yàn)證本文分解方法的有效性.圖3為采用本文解析時(shí)間波場(chǎng)外推方法得到的CDP號(hào)為380位置處的解析VSP記錄，其中炮點(diǎn)布設(shè)在CDP號(hào)為300位置處.采用常規(guī)波場(chǎng)外推方法僅能得到該解析VSP記錄的實(shí)部.圖4為常規(guī)實(shí)波場(chǎng)和解析波場(chǎng)在0.5s時(shí)間片上的垂直波數(shù)譜對(duì)比.由圖3可以看出該記錄在0.5s時(shí)間片上僅存在下行波記錄.而常規(guī)實(shí)波場(chǎng)的波數(shù)譜并不能指示這一信息，所以常規(guī)波場(chǎng)分離方法通常在頻率域?qū)崿F(xiàn).但在本文解析時(shí)間波場(chǎng)的波數(shù)譜上則可以清楚地識(shí)別該下行波記錄.利用本文波場(chǎng)分解方法在每個(gè)時(shí)間片上進(jìn)行上下行波分解可以得到分離的上下行波記錄，如圖5所示.

在成像空間應(yīng)用上述波場(chǎng)分解方法將波場(chǎng)分解為各個(gè)方向矢量，然后進(jìn)行成像，得到的炮點(diǎn)下行和檢波點(diǎn)上行成像結(jié)果如圖6a所示，炮點(diǎn)上行和檢波

圖3 解析VSP記錄(a)解析VSP記錄實(shí)部； (b)解析VSP記錄虛部.Fig.3 Analytic VSP record(a) Real part of analytic VSP record; (b) Imaginary part of analytic VSP record.

圖4 波數(shù)譜對(duì)比圖(a)常規(guī)實(shí)波場(chǎng)波數(shù)譜；(b)解析波場(chǎng)波數(shù)譜.Fig.4 Comparison of wavenumber spectrums(a) Wavenumber spectrum of conventional real wavefield; (b) Wavenumber spectrum of analytic wavefield.

圖5 分離的上下行波(a)下行波記錄；(b)上行波記錄.Fig.5 Separated up and down going wavefield(a) Down-going wavefield; (b) Up-going wavefield.

圖6 波場(chǎng)分離成像結(jié)果(a)炮點(diǎn)下行與檢波點(diǎn)上行成像結(jié)果；(b)炮點(diǎn)上行與檢波點(diǎn)下行成像結(jié)果；(c)炮點(diǎn)下行與檢波點(diǎn)下行成像結(jié)果；(d)炮點(diǎn)上行與檢波點(diǎn)上行成像結(jié)果.Fig.6 Results of wavefield separation imaging(a) Result from down going source and up going receiver wavefield; (b) Result from up going source and down going receiver wavefield; (c) Result from down going source and down going receiver wavefield; (d) Result from up going source and up going receiverwavefield.

點(diǎn)下行成像結(jié)果如圖6b所示，炮點(diǎn)下行和檢波點(diǎn)下行成像結(jié)果如圖6c所示，炮點(diǎn)上行和檢波點(diǎn)上行成像結(jié)果如圖6d所示.由圖可以看出，圖6a和圖6b為有效信號(hào)，而圖6c和圖6d為低頻噪音.由圖也可以看出，本文方法可以實(shí)現(xiàn)各個(gè)成像分量的分離.在這種速度模型比較精確的情況，可以疊加圖6a和圖6b結(jié)果作為最終成像結(jié)果，也即等效波場(chǎng)分離成像公式結(jié)果.在這種情況下，等效波場(chǎng)分離成像條件的結(jié)果和本文方法結(jié)果等價(jià).

在實(shí)際情況下，由于速度建模方法對(duì)速度結(jié)構(gòu)的分辨能力限制，通常只能得到一個(gè)相對(duì)平滑的速度模型.將準(zhǔn)確速度模型進(jìn)行平滑來(lái)模擬實(shí)際情況，平滑速度模型如圖7所示.采用常規(guī)相關(guān)逆時(shí)偏移成像，結(jié)果如圖8a所示.由圖可以看出，其結(jié)果存在一定的低頻噪音，還存在一定的成像假象.采用本文方法提取炮點(diǎn)下行和檢波點(diǎn)上行成像，炮點(diǎn)上行和檢波點(diǎn)下行成像的疊加結(jié)果如圖8b所示，該結(jié)果等價(jià)于等效波場(chǎng)分離成像公式得到的結(jié)果.由圖8b可以看出，相比于圖8a，成像假象明顯減少，但是在向斜的中心存在一個(gè)較強(qiáng)的成像假象，且該假象呈現(xiàn)與向斜構(gòu)造相似的形態(tài).這個(gè)假象產(chǎn)生的主要原因?yàn)椋河捎诘诙雍偷谌龑又g存在強(qiáng)速度分界面，因此平滑速度模型中存在強(qiáng)梯度速度過(guò)渡帶，當(dāng)炮點(diǎn)和檢波點(diǎn)波場(chǎng)經(jīng)過(guò)這一區(qū)域時(shí)，都會(huì)產(chǎn)生向上回傳的波場(chǎng)，而由于向斜獨(dú)特的地質(zhì)形狀，導(dǎo)致炮點(diǎn)上行波和檢波點(diǎn)下行波場(chǎng)都匯聚到向斜中心，進(jìn)而在應(yīng)用成像條件時(shí)產(chǎn)生假象.由此，可以看出這類假象是由炮點(diǎn)上行和檢波點(diǎn)下行波場(chǎng)進(jìn)行相關(guān)所導(dǎo)致的.采用本文方法分別提取炮點(diǎn)下行和檢波點(diǎn)上行成像如圖8c所示，炮點(diǎn)上行和檢波點(diǎn)下行成像結(jié)果如圖8d所示.由圖可以看出，圖8c為有效信號(hào)，而圖8d則是圖8b所包含的假象.在實(shí)際情況下，這種假象可能導(dǎo)致巖丘邊界、潛山面及火成巖邊界的錯(cuò)誤解釋.由圖8c和圖8d也可以看出在這種情況下，本文方法相比于等效波場(chǎng)分離成像條件能夠更好地分離各個(gè)成像分量，進(jìn)而消除假象.

圖7 平滑速度模型Fig.7 Smooth velocity model

圖8 成像結(jié)果對(duì)比(a)傳統(tǒng)成像結(jié)果；(b)等效波場(chǎng)分離成像結(jié)果；(c)炮點(diǎn)下行與檢波點(diǎn)上行成像結(jié)果；(d)炮點(diǎn)上行與檢波點(diǎn)下行成像結(jié)果.Fig.8 Imaging results comparison(a) Result of conventional RTM; (b) Result from equivalent wavefield imaging condition; (c) Result from down going source and up going receiver wavefield; (d) Result from up going source and down going receiver wavefield.

圖9 不同偏移方法成像結(jié)果的對(duì)比(a)相關(guān)成像結(jié)果；(b)Laplacian濾波方法；(c)波場(chǎng)方向分解方法.Fig.9 Comparison of imaging results from different methods(a) Result of correlation imaging condition; (b) Result from Laplacian filter; (c) Result from wavefield separation method.

下面對(duì)比本文方法與常規(guī)濾波方法在噪音消除方面的表現(xiàn).Laplacian濾波方法由于其高效性而被廣泛應(yīng)用于低頻噪音消除中，因此，這里將本文方法與該方法進(jìn)行對(duì)比.測(cè)試的模型為一個(gè)兩層速度模型，速度分界面在1500m處，第一層速度為2000m·s-1，第二層速度為2500m·s-1.由該模型速度分界面的速度差異可以得出：在該模型上，低頻成像噪音強(qiáng)度較弱，同時(shí)成像噪音僅分布在成像子波之上.這樣，利用相關(guān)成像方法得到的成像子波就可近似地作為其他成像方法在成像振幅保真度方面的參考.圖9展示了相關(guān)成像方法、Laplacian濾波方法以及本文波場(chǎng)分解方法的結(jié)果.由圖可以看出，Laplacian濾波方法和本文波場(chǎng)分解方法都較好地對(duì)低頻背景噪音進(jìn)行了去除.但是對(duì)比二者能量關(guān)系可以發(fā)現(xiàn)，Laplacian濾波方法顯著地改變了成像振幅.為了更清楚地進(jìn)行對(duì)比，下面抽取各種方法在CDP號(hào)為300位置處的單道成像結(jié)果進(jìn)行對(duì)比，如圖10所示.由圖10可以看出，常規(guī)相關(guān)方法由于在反射層之上受到低頻背景噪音的影響，其成像振幅小于零.Laplacian濾波方法和本文方法的結(jié)果在反射層之上的成像振幅趨近于零，由此可以看出它們都較好地實(shí)現(xiàn)了噪音的消除.但是進(jìn)一步對(duì)比成像子波可以發(fā)現(xiàn)，本文波場(chǎng)分解方法的成像子波與相關(guān)方法更好地貼合，而Laplacian濾波方法的成像子波能量則存在顯著差異.由此可以看出本文方法在成像子波保真性方面具有較好的表現(xiàn).

圖10 不同成像結(jié)果的單道對(duì)比(紅色代表相關(guān)成像結(jié)果；綠色代表Laplacian濾波結(jié)果；藍(lán)色代表波場(chǎng)方向分解結(jié)果)Fig.10 Trace comparison of imaging result from different methods

圖11 實(shí)際數(shù)據(jù)偏移速度模型Fig.11 Migration velocity model of real data

圖12 實(shí)際數(shù)據(jù)成像結(jié)果對(duì)比(a) 傳統(tǒng)逆時(shí)偏移結(jié)果； (b) 波場(chǎng)分離成像結(jié)果.Fig.12 Imaging results comparison of real data(a) Result of conventional RTM; (b) Result of wavefield separation.

下面采用某地區(qū)實(shí)際數(shù)據(jù)進(jìn)行方法測(cè)試，偏移采用的速度模型為層析速度反演結(jié)果，如圖11所示.在速度模型中可以看出較強(qiáng)的速度分界面，因此在常規(guī)逆時(shí)偏移結(jié)果(如圖12a)上存在較強(qiáng)的低頻噪音.采用本文方法提取炮點(diǎn)下行和檢波點(diǎn)上行成像，炮點(diǎn)上行和檢波點(diǎn)下行成像的結(jié)果，如圖12b所示.由圖可以看出，低頻噪音得到了很好的消除，同時(shí)產(chǎn)生低頻噪音區(qū)域的構(gòu)造更加清楚.

5 結(jié)論

有效的炮點(diǎn)及檢波點(diǎn)波場(chǎng)分解及相應(yīng)波場(chǎng)的相關(guān)成像是高質(zhì)量逆時(shí)深度偏移的核心.本文提出了解析時(shí)間波場(chǎng)外推及波場(chǎng)分解方法，能夠在逆時(shí)深度偏移外推過(guò)程中的每個(gè)時(shí)間片上僅應(yīng)用空間傅里葉變換實(shí)現(xiàn)波傳播方向的有效分解，進(jìn)而實(shí)現(xiàn)相應(yīng)波場(chǎng)的相關(guān)成像，得到高質(zhì)量的逆時(shí)深度偏移成像結(jié)果，低頻噪音和成像假象被有效壓制.相比于頻率波數(shù)域波場(chǎng)分離方法，該方法具有較高的實(shí)現(xiàn)效率，與常規(guī)的相關(guān)+低頻濾波的逆時(shí)深度偏移方法相比，增加的計(jì)算量較少.當(dāng)然，該方法也可以應(yīng)用于方位角度成像道集的生成，為后續(xù)的基于逆時(shí)深度偏移的速度層析反演和AVA分析與反演提供高質(zhì)量的工具.

Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration.Geophysics, 48(11): 1514-1524.

Du Q Z, Zhu Y T, Zhang M Q, et al. 2013. A study on the strategy of low wavenumber noise suppression for prestack reverse-time depth migration.ChineseJ.Geophys. (in Chinese), 56(7): 2391-2401, doi: 10.6038/cjg20130725.Fei T W, Luo Y, Aramco S, et al. 2010. De-blending reverse-time migration. ∥SEG Expanded Abstracts, 3130-3134.

Fei T W, Luo Y, Qin F H. 2014. An endemic problem in reverse-time migration. ∥SEG Expanded Abstracts, 3811-3815.

Fletcher R P, Fowler P J, Kitchenside P, et al. 2006. Suppressing unwanted internal reflections in prestack reverse-time migration.Geophysics, 71(6): E79-E82.Guitton A, Kaelin B, Biondi B. 2007. Least-squares attenuation of reverse-time migration artifacts.Geophysics, 72(1): S19-S23.

Kang W, Cheng J B. 2012. Methods of suppressing artifacts in prestack reverse time migration.ProgressinGeophys. (in Chinese), 27(3): 1163-1172.Li B, Liu H W, Liu G F, et al. 2010. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU.ChineseJ.Geophys. (in Chinese), 53(12): 2938-2943, doi: 10.3969/j.issn.0001-5733.2010.12.017.

Liu F Q, Hanson D W, Whitmore N D, et al. 2006. Toward a unified analysis for source plane-wave migration.Geophysics, 71(4): S129-S139.

Liu F Q, Zhang G Q, Morton S A, et al. 2007. Reverse-time migration using one-way wavefield imaging condition. ∥SEG Expanded Abstracts, 26: 2170-2174.

Liu F Q, Zhang G Q, Morton S A, et al. 2011. An effective imaging condition for reverse-time migration using wavefield decomposition.Geophysics, 76(1): S29-S39.

Liu H W, Li B, Liu H, et al. 2010. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation.ChineseJ.Geophys. (in Chinese), 53(7): 1725-1733, doi: 10.3969/j.issn.0001-5733.2010.07.024.

Liu S W, Wang H Z, Chen S C, et al. 2013. Implementation strategy of 3D reverse time migration on GPU/CPU clusters.ChineseJ.Geophys. (in Chinese), 56(10): 3487-3496, doi: 10.6038/cjg20131023.

Loewenthal D, Stoffa P L, Faria E L. 1987. Suppressing the unwanted reflections of the full wave equation.Geophysics, 52

(7): 1007-1012.

McMechan G A. 1983. Migration by extrapolation of time-dependent boundary values.GeophysicalProspecting, 31(3): 413-420.

Whitmore N D. 1983. Iterative depth migration by backward time propagation. ∥SEG Expanded Abstracts, 382-385.

Yoon K, Marfurt K J. 2006. Reverse-time migration using the Poynting vector.ExplorationGeophysics, 37(1): 102-107.

Zhang Y, Sun J. 2009. Practical issues in reverse time migration: true amplitude gathers, noise removal and harmonic source encoding.FirstBreak, 26(1): 29-35.

附中文參考文獻(xiàn)

杜啟振, 朱釔同, 張明強(qiáng)等. 2013. 疊前逆時(shí)深度偏移低頻噪聲壓制策略研究. 地球物理學(xué)報(bào), 56(7): 2391-2401, doi: 10.6038/cjg20130725.

康瑋, 程玖兵. 2012. 疊前逆時(shí)偏移假象去除方法. 地球物理學(xué)進(jìn)展, 27(3): 1163-1172, 10.6038/j.issn.1004-2903.2012.03.041.

李博, 劉紅偉, 劉國(guó)峰等. 2010. 地震疊前逆時(shí)偏移算法的CPU/GPU實(shí)施對(duì)策. 地球物理學(xué)報(bào), 53(12): 2938-2943, doi: 10.3969/j.issn.0001-5733.2010.12.017.

劉紅偉, 李博, 劉洪等. 2010. 地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn). 地球物理學(xué)報(bào), 53(7): 1725-1733, doi: 10.3969/j.issn.0001-5733.2010.07.024.

劉守偉, 王華忠, 陳生昌等. 2013. 三維逆時(shí)偏移GPU/CPU機(jī)群實(shí)現(xiàn)方案研究. 地球物理學(xué)報(bào), 56(10): 3487-3496, doi: 10.6038/cjg20131023.

(本文編輯 胡素芳)

Reverse time migration using analytical time wavefield extrapolation and separation

HU Jiang-Tao, WANG Hua-Zhong

WPI,SchoolofOceanandEarthScience,TongjiUniversity,Shanghai200092，China

Two-way wave equation based reverse time migration is powerful in imaging complex area. However, its result usually contains imaging noise which poses challenge to seismic interpretation. The imaging amplitude would be distorted by simple noise filter. Wavefield separation imaging condition is an effective way to suppress the imaging noise and enhance the imaging quality. Conventionally, the wavefield separation is achieved in frequency wavenumber domain which cannot be conveniently implemented in time domain reverse time migration. We propose the analytical wavefield extrapolation method which can efficiently separate wavefield into its directional components during the time domain extrapolation. We use this method in reverse time migration and separate both the source and receiver wavefield. And the imaging noise and signal can be separated after applying the imaging condition to the separated wavefield.The wavefield propagation direction is usually defined by temporal and spatial Fourier transform. In the Fourier domain, the wavefield propagation direction is defined by the sign of frequency and spatial wavenumber. In order to efficiently apply the wavefield separation imaging condition in reverse time migration, we extend the analytical time signal into the wavefield and call it analytical time wavefield. Since the analytical time wavefield contains only the positive frequency component, the wavefield propagation direction can be defined by the sign of spatial wavenumber. To avoid the I/O cost in generating the analytical wavefield at every imaging point, we propose an analytical wavefield propagation equation based on the linear relation between the source term and wavefield in wave equation. We solve the proposed equation by finite difference method. Then we separate the source and receiver wavefield into their up and down going components and apply the imaging condition to the separated wavefield. Four imaging components (i.e., imaging component from up going source and down going receiver wavefield, imaging component from down going source and down going receiver wavefield, imaging component from down going source and up going receiver wavefield and imaging component from up going source and up going receiver wavefield) are effectively separated. Then the imaging noise and signal are separated.We use the syncline model and real data to test the proposed method. Numerical example on syncline model shows that the proposed method can effectively separate imaging component from down going source and up going receiver wavefield, imaging component from up going source and down going receiver wavefield, imaging component from down going source and down going receiver wavefield and imaging component from up going source and up going receiver wavefield. When the migration velocity is not accurate, the correlation imaging condition and equivalent wavefield separation imaging condition would generate imaging artifact in central area of the syncline. And the proposed method can eliminate the artifact. Numerical example on real data shows that this method can generate imaging result free from low frequency noise.Wavefield separation imaging condition can effectively separate the imaging noise and signal in reverse time migration and enhance the imaging quality. The proposed analytical wavefield extrapolation method can separate the source and receiver wavefield into their up and down going components. The imaging noise can be suppressed by applying the imaging condition to the separated wavefield. We want to further use this method to generate common angle imaging gather in reverse time migration and provide input data for AVA and migration velocity inversion.

Analytical time wavefield extrapolation; Wavefield separation; Noise elimination; Reverse time migration

973項(xiàng)目(2011CB201002)，國(guó)家自然科學(xué)基金(41374117)以及國(guó)家重大專項(xiàng)(2011ZX05005-005-008HZ,2011ZX05006-002,2011ZX05023)資助.

胡江濤，男，1987年生，博士研究生，主要研究方向?yàn)榈卣鸩▊鞑ヅc成像方法及應(yīng)用.E-mail:huidaojiaxiang@126.com

10.6038/cjg20150822.

10.6038/cjg20150822

P631

2015-01-08，2015-07-22收修定稿

胡江濤, 王華忠.2015.基于解析時(shí)間波場(chǎng)外推與波場(chǎng)分解的逆時(shí)偏移方法研究.地球物理學(xué)報(bào)，58(8):2886-2895,

Hu J T, Wang H Z. 2015. Reverse time migration using analytical time wavefield extrapolation and separation.ChineseJ.Geophys. (in Chinese),58(8):2886-2895,doi:10.6038/cjg20150822.