基于波長平滑算子的逆時偏移回轉(zhuǎn)波假象去除方法

唐永杰 宋桂橋 劉少勇 顧漢明* 嚴 哲

(①中國地質(zhì)大學(武漢)地球物理與空間信息學院，湖北武漢 430074； ②地球內(nèi)部多尺度成像湖北省重點實驗室，湖北武漢 430074； ③中國石化油田勘探開發(fā)事業(yè)部，北京 100728)

1 引言

隨著地震勘探目標越來越復雜、越來越精細，基于射線理論和單程波方程理論的地震偏移技術(shù)已難以滿足業(yè)界需求[1-3]； 而采用雙程波動方程的逆時偏移(RTM)具有成像精度高、相位準確、適應(yīng)于變速介質(zhì)和陡傾構(gòu)造成像的特點[4，5]，并且可對繞射波、多次波、棱柱波等進行成像[6-8]，因此，RTM已成為目前復雜構(gòu)造成像的常用方法之一。

RTM方法自20世紀80年代被提出以來[9-11]，主要面臨計算效率和偏移假象兩大問題的挑戰(zhàn)。隨著計算機技術(shù)的飛速發(fā)展以及各種RTM優(yōu)化算法和并行策略的提出[12-20]，計算效率問題基本得到了解決。針對偏移假象問題，傳統(tǒng)相關(guān)成像條件中來自炮檢點同向傳播波場產(chǎn)生的低頻噪聲受到廣泛關(guān)注，并涌現(xiàn)出多種消除此類低頻噪聲的方法。Loewenthal等[21]采用平滑速度方式消除背向散射波場進而壓制低頻噪聲； Mulder等[22]通過帶通濾波器對偏移后的成像結(jié)果進行濾波處理以消除低頻假象； Valenciano等[23]利用反褶積成像條件壓制低頻噪聲； Fletcher等[24]提出采用方向性衰減方法壓制層間反射從而去除RTM中的低頻假象； Yoon等[25]將Poynting矢量引入成像條件以消除低頻噪聲進而提高成像質(zhì)量； Guitton等[26]采用最小平方濾波法去除低頻噪聲； Liu等[27]先將全波場分解為各自的單程波分量，再選取合適分量進行互相關(guān)成像； Zhang等[28]通過切除角度域成像道集中的大角度部分，再做角度疊加，以此壓制低頻假象； Liu等[29]提出一種基于Hilbert變換的波場分解法，該方法不需將波場顯式分解，也不需存儲整個波場，只需選擇性地對有效成像路徑進行互相關(guān)從而達到去除低頻噪聲的目的； Rocha等[30]提出基于能量范數(shù)的成像條件來消除低頻噪聲。

除了常見的低頻假象外，當偏移速度梯度較大或反射界面較復雜時，RTM結(jié)果中會產(chǎn)生回轉(zhuǎn)波假象，該假象不能像低頻噪聲那樣通過簡單的濾波手段去除。Fei等[31]提出可將源、檢波場進行分解，僅利用震源波場的下行波與檢波點波場的上行波做互相關(guān)，從而消除該假象。此類方法采用類似于單程波偏移的成像條件時會對高陡傾構(gòu)造成像產(chǎn)生一定影響，且震源(檢波點)波場的顯式分解需在頻率—波數(shù)域進行，其計算較復雜。

本文引入上、下行波分離成像條件分析RTM中產(chǎn)生假象的原因，認為回轉(zhuǎn)波假象是在強速度梯度帶由震源波場上行波與檢波點波場下行波互相關(guān)產(chǎn)生，進而提出一種基于波長平滑算子的回轉(zhuǎn)波假象去除方法。該方法避免波場分解，通過引入波長平滑算子，使波場在強速度梯度處的傳播特征更接近于真實速度中波場傳播特性，進而壓制回轉(zhuǎn)波假象，因此可較好地兼顧效率和高陡構(gòu)造成像。

2 RTM中回轉(zhuǎn)波假象成因分析

逆時偏移過程包含炮檢波場的正反傳播和成像條件施加兩部分，其中常用的互相關(guān)成像條件可表示為

(1)

式中：x為地下成像點位置；t為炮檢點傳播到成像點的時間；s(t,x)和r(t,x)分別為時、空域炮檢點波場；Tmax為最大記錄時間；I(x)為成像結(jié)果。RTM的實現(xiàn)是基于雙程波動算子，故在波場外推時會產(chǎn)生反射波和透射波。Liu等[27]提出將波場進行上、下行波分解，則震源波場s(t,x)可分解為上行波su(t,x)和下行波sd(t,x)，檢波點波場r(t,x)可分解為上行波ru(t,x)和下行波rd(t,x)，即炮檢點波場可表示為

s(t,x)=su(t,x)+sd(t,x)

(2)

r(t,x)=ru(t,x)+rd(t,x)

(3)

式中 u(上行)和d(下行)代表波的傳播方向。將式(2)和式(3)代入式(1)，可得

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

(4)

式中：I1(x)表示震源波場下行波與檢波點波場上行波互相關(guān)成像結(jié)果；I2(x)表示震源波場上行波與檢波點波場下行波互相關(guān)成像結(jié)果；I3(x)表示震源波場下行波與檢波點波場下行波互相關(guān)成像結(jié)果；I4(x)表示震源波場上行波與檢波點波場上行波互相關(guān)成像結(jié)果。在當前反射波偏移成像框架下,I1(x)和I2(x)對應(yīng)震源波場和檢波點波場反向傳播的波場進行互相關(guān)，都可對地下反射界面進行成像。具體而言：I1(x)類似于單程波偏移的成像條件，結(jié)果即對應(yīng)一次反射波成像；I2(x)即對應(yīng)能量較一次波弱的多次波成像；I3(x)和I4(x)對應(yīng)震源波場和檢波點波場同向傳播的波場進行互相關(guān)，即表現(xiàn)為低頻假象，一般通過拉普拉斯濾波可較徹底地消除。當震源波場和檢波點波場通過強速度梯度或復雜反射界面時，將產(chǎn)生向上回傳的波場(圖1)，此時I2(x)則表現(xiàn)為一次回轉(zhuǎn)波產(chǎn)生的假象，其能量強于多次波成像，且難以被去除。

圖1 回轉(zhuǎn)波假象產(chǎn)生示意圖

為進一步驗證該回轉(zhuǎn)波假象是由I2(x)(震源波場上行波與檢波點波場下行波互相關(guān)結(jié)果)分量產(chǎn)生，選取一個兩層模型進行測試。模型上、下層速度分別為2000m/s和4000m/s(圖2a)，其偏移速度為該兩層模型做高斯平滑后速度模型(圖2b)。

基于該模型，將震源波場和檢波點波場分解為各自的上、下行波，分別選取對應(yīng)的波場分量進行成像(圖3)。圖3a表示式(4)中I1(x)和I2(x)相加的成像結(jié)果，可見當偏移速度具有較強的速度梯度時，成像結(jié)果中真實反射界面上方出現(xiàn)了假的同相軸(紅色箭頭)； 圖3b表示I1(x)震源波場下行波與檢波點波場上行波互相關(guān)成像結(jié)果，對應(yīng)反射界面的成像； 圖3c表示I2(x)震源波場的上行波與檢波點波場的下行波互相關(guān)成像結(jié)果，對比可知I2(x)成像結(jié)果即是圖3a中回轉(zhuǎn)波假象。

圖2 兩層介質(zhì)速度模型 (a)準確速度； (b)高斯光滑速度

圖3 采用高斯光滑偏移速度得到的不同分量偏移結(jié)果 (a)I1(x)+I2(x)； (b)I1(x)； (c)I2(x)

3 用波長平滑算子去除回轉(zhuǎn)波假象

通過以上分析可知，在采用互相關(guān)成像條件進行RTM成像時，在強速度梯度處會產(chǎn)生回轉(zhuǎn)波，相關(guān)成像結(jié)果表現(xiàn)為類似于一次反射波成像結(jié)果的假同相軸，該假象沒有典型的低頻噪聲特征，不能通過簡單疊后濾波方法去除。若不考慮多次波成像，且引入波場分解成像條件，僅利用式(4)中I1(x)(震源波場下行波與檢波點波場上行波互相關(guān)結(jié)果)分量成像即可較好地回避此假象[31]。但該方法需對波場進行顯式分解，而一般的波場分解法宜在頻率—波數(shù)域進行，需將整個波場進行存盤，并進行關(guān)于時間最慢維的高維傅里葉變換，其計算量和存儲量都較大； 另一方面該方法也損失了RTM相對于單程波偏移的優(yōu)勢，因其僅是對波場進行上、下行波分離，則對陡傾角構(gòu)造的成像有一定的局限性。

在強變速介質(zhì)中，對速度場做簡單平滑(如高斯平滑)會人為引入強速度梯度帶而產(chǎn)生回轉(zhuǎn)波，此時常規(guī)相關(guān)成像條件不能對回轉(zhuǎn)波進行正確成像。如果控制波場在傳播過程中只產(chǎn)生符合波傳播規(guī)律的透射和反射波場，消除回轉(zhuǎn)波，則在互相關(guān)成像時該回轉(zhuǎn)波假象就能被壓制或消除?；诓ㄩL平滑算子的速度平滑法在平滑處理時考慮地震波長因素，使波場更接近于在真實速度中的傳播特性，以壓制在強速度梯度處產(chǎn)生的向上回傳波場，進而壓制或消除由此產(chǎn)生的回轉(zhuǎn)波假象。在波長相關(guān)平滑中，每一網(wǎng)格點的平均速度為距離該點一個波長內(nèi)所有點速度的加權(quán)平均。對于某一網(wǎng)格點，其波長相關(guān)平滑速度為

(5)

(6)

式中：t表示一個波長內(nèi)計算網(wǎng)格點到加權(quán)點的直線旅行時；T表示對應(yīng)的周期；p為速度場的維度。

與采用固定窗口的高斯平滑算子相比，本文平滑算子通過引入波長相關(guān)的加權(quán)函數(shù)對速度場進行波長相關(guān)平滑，使波場在界面處更符合其物理傳播特征。為凸顯二者區(qū)別，分別采用兩種平滑算子對上述兩層速度模型進行平滑。圖4a～圖4c依次對應(yīng)準確速度、高斯平滑后速度和波長相關(guān)平滑后速度縱向變化三種情形。對比發(fā)現(xiàn)圖4b中界面上下速度梯度一致，而圖4c中界面上下速度變化率與波長相關(guān)。

為了進一步分析高斯平滑與波長相關(guān)平滑對波

圖4 不同速度的縱向變化 (a)準確速度； (b)高斯平滑速度； (c)波長平滑速度

傳播的影響，利用有限差分算子針對該模型做波場模擬。圖5為采用主頻為30Hz震源、以不同平滑速度在不同時刻的波場快照。其中紅色虛線代表真實速度模型中波傳播的準確旅行時。

通過對比波場快照可見：圖5a和圖5d中紅色虛線與真實速度中的波前面相吻合； 圖5b和圖5e中高斯平滑速度場的波前較準確時間更快(紅色箭頭)，這是由于高斯平滑在界面處速度變化率是對稱的，波場通過平滑區(qū)域所用時間更短；圖5c和圖5f中波長相關(guān)平滑速度場的波前與紅色虛線表示的準確旅行時一致，這是由于波長相關(guān)平滑考慮了波長因素，平滑后保持了強速度界面變化特征，使波場更接近于準確速度的傳播特性，該特性確保不會因平滑而過度降低逆時偏移精度，同時兼顧達到去除回轉(zhuǎn)波假象的目的。

圖5 基于不同平滑速度模型的波場快照 (a)、(b)和(c)分別為真實速度、高斯平滑速度和30Hz波長相關(guān)平滑速度的0.7s波場快照； (d)、(e)和(f)為對應(yīng)的0.85s波場快照

從上述試算結(jié)果及其分析可知，波長平滑算子能改善地震波在復雜介質(zhì)中的傳播，這是在逆時偏移互相關(guān)成像中壓制回轉(zhuǎn)波假象的基礎(chǔ)。

4 數(shù)值算例

采用洼陷模型測試本文提出的波長平滑算子對RTM中回轉(zhuǎn)波假象的去除效果。圖6a所示準確速度模型的尺寸為5000m(橫)×2000m(縱)。對該模型進行有限差分模擬的參數(shù)為：子波是主頻為20Hz的雷克子波，道距為10m， 炮檢距范圍是-5000～5000m，采樣間隔為1ms，采樣長度為2.5s，炮間距為50m，共101炮。圖6b和圖6c分別為高斯平滑速度場和20Hz波長平滑速度場。圖7a為采用高斯平滑速度獲得的RTM成像結(jié)果，可見洼陷兩側(cè)出現(xiàn)假同相軸(紅色箭頭)； 圖7b為采用20Hz波長平滑速度獲得的RTM成像結(jié)果，易見洼陷兩側(cè)假象已消除。

采用復雜鹽丘模型[31]進一步測試本文方法對復雜地質(zhì)體適用性。圖8a所示準確速度模型的尺寸為20km(橫)×7.8km(縱)。模型中部為火山巖侵入構(gòu)造，該鹽丘具有不規(guī)則邊界，其速度高于圍巖。對該模型進行有限差分模擬的參數(shù)為：子波是主頻為20Hz的雷克子波，道距為20m，炮檢距范圍是-10～10km，采樣間隔為4ms，采樣長度為10s，炮間距為80m，共251炮。圖8b和圖8c分別為高斯平滑速度場和20Hz波長平滑速度場。分別采用準確速度場、高斯平滑速度場和20Hz波長平滑速度場進行逆時偏移成像(圖9)。在采用準確速度場得到的偏移結(jié)果(圖9a)中可見鹽丘頂部不存在強速度梯度帶，不產(chǎn)生回轉(zhuǎn)波，即鹽丘頂界面并未產(chǎn)生假象；而在采用高斯平滑偏移速度得到的偏移結(jié)果(圖9b)中看到鹽丘頂部出現(xiàn)了假同相軸(紅色箭頭)，該假象頻帶與有效信號頻帶相近，且不能通過濾波去除； 在采用波長相關(guān)平滑速度場得到的偏移結(jié)果(圖9c)中，可見鹽丘頂部回轉(zhuǎn)波假象(相對圖9b的高斯平滑速度偏移結(jié)果)被消除，這有助于對鹽丘邊界、潛山面和火成巖邊界的地震解釋，而且鹽丘兩側(cè)陡構(gòu)造成像效果也更好(紅色橢圓)；另外，采用波長平滑速度的偏移結(jié)果(圖9c)在成像精度上并不明顯低于準確速度模型成像結(jié)果(圖9a)。

圖6 不同速度場的洼陷模型 (a)準確速度場； (b)高斯平滑速度場； (c)20Hz波長平滑速度場

圖7 不同偏移速度獲得的洼陷模型成像結(jié)果 (a)高斯平滑速度； (b)波長平滑速度

圖8 不同速度場的鹽丘模型 (a)準確速度場； (b)高斯平滑速度場； (c)20Hz波長平滑速度場

5 結(jié)論與討論

在地下強速度梯度區(qū)域，基于互相關(guān)成像條件的RTM會產(chǎn)生回轉(zhuǎn)波假象，該假象區(qū)別于典型的低頻噪聲，有較強的能量且不容易去除。本文通過波場分解成像條件分析回轉(zhuǎn)波假象產(chǎn)生原因，認為該假象是緣于震源波場上行波與檢波點波場下行波的互相關(guān)，即強梯度帶的回轉(zhuǎn)波在傳統(tǒng)相關(guān)成像條件下產(chǎn)生。據(jù)此提出一種基于波長平滑算子的回轉(zhuǎn)波假象去除方法：通過波長相關(guān)平滑算子進行速度平滑，在波傳播過程中壓制回轉(zhuǎn)波的產(chǎn)生，進而消除回轉(zhuǎn)波假象。該方法區(qū)別于通過波場分解法消除回轉(zhuǎn)波成像假象，可以很好地對高陡構(gòu)造進行成像且不增加RTM的計算復雜度。

[1] Yoon K，Marfurt K J，Starr W.Challenges in reverse-time migration.SEG Technical Program Expanded Abstracts,2004，23：1057-1060.

[2] Zhu J，Lines L R.Comparison of Kirchhoff and reverse time migration methods with application to prestack depth imaging of complex structures.Geophysics，1998，63(4)：1166-1176.

[3] Sava P，Hill S J.Overview and classification of wavefield seismic imaging methods.The Leading Edge，2009，28(2)：170-183.

[4] 楊午陽，張厚柱，撒利明等. 逆時偏移關(guān)鍵問題探討. 石油地球物理勘探，2016，51(6)：1251-1262. Yang Wuyang，Zhang Houzhu，Sa Liming et al.Some issues about reverse time migration.OGP，2016，51(6)：1251-1262.

[5] 王鵬飛，何兵壽.基于行波分離的三維彈性波矢量場點積互相關(guān)成像條件.石油地球物理勘探，2017，52(3)：477-483. Wang Pengfei，He Bingshou.Vector field dot product cross-correlation imaging based on 3D elastic wave separation.OGP，2017，52(3)：477-483.

[6] 葉月明，郭慶新，孫小東等.多次波對逆時偏移構(gòu)造成像的影響.石油地球物理勘探，2014，49(3)：523-529. Ye Yueming，Guo Qingxin，Sun Xiaodong et al.Multiples influence on structure imaging in reverse time migration.OGP，2014，49(3)：523-529.

[7] Deeks J，Lumley D.Prism waves in seafloor canyons and their effects on seismic imaging.Geophysics，2015，80(6)：S213-S222.

[8] Weglein A B.Multiples：Signal or noise? Geophysics，2016，81(4)：V283-V302.

[9] Whitmore D N.Iterative depth migration by backward time propagation.SEG Technical Program Expanded Abstracts，1983，2：382-385.

[10] Baysal E，Kosloff D D，Sherwood J W C.Reverse time migration.Geophysics，1983，48(11)：1514-1524.

[11] McMechan G A.Migration by extrapolation of time-dependent boundary values.Geophysical Prospecting，1983，31(3)：413-420.

[12] Fricke J R.Reverse-time migration in parallel：A tutorial.Geophysics，1988，53(9)：1143-1150.

[13] Harris C E，McMechan G A.Using downward continu-ation to reduce memory requirements in reverse-time migration.Geophysics，1992，57(6)：848-853.

[14] Symes W.Reverse time migration with optimal check-pointing.Geophysics，2007，72(5)：SM213-SM221.

[15] 劉紅偉，李博，劉洪等.地震疊前逆時偏移高階有限差分算法及GPU實現(xiàn).地球物理學報，2010，53(7)：1725-1733. Liu Hongwei，Li Bo，Liu Hong et al.The algorithm of high order finite difference prestack reverse time migration and GPU implementation.Chinese Journal of Geophysics，2010，53(7)：1725-1733.

[16] Anderson J E，Tan L，Wang D.Time-reversal check pointing methods for RTM and FWI.Geophysics，2012，77(4)：S93-S103.

[17] 張慧，蔡其新，秦廣勝等.基于GPU并行加速的逆時偏移成像方法.石油地球物理勘探，2013，48(5)：707-710. Zhang Hui，Cai Qixin，Qin Gunagsheng et al.GPU-accelerated reverse time migration and its application.

OGP，2013，48(5)：707-710.

[18] 唐祥功，匡斌，杜繼修等.多GPU 協(xié)同三維疊前逆時偏移方法研究與應(yīng)用.石油地球物理勘探，2013，48(6)：910-914. Tang Xianggong，Kuang Bin，Du Jixiu et al.3D pre-stack reverse time migration based on multiple-GPU collaboration.OGP，2013，48(6)：910-914.

[19] Clapp R G.Seismic processing and the computer revolution.SEG Technical Program Expanded Abstracts，2015，34：4832-4837.

[20] Yang P，Brossier R，Virieux J.Wavefield reconstruction by interpolating significantly decimated boundaries.Geophysics，2016，81(5)：T197-T209.

[21] Loewenthal D，Stoffa P L，F(xiàn)aria E L.Suppressing the unwanted reflections of the full wave equation.Geophysics，1987，52(2)：1007-1012.

[22] Mulder W，Plessix R.One-way and two-way wave-equation migration.SEG Technical Program Expanded Abstracts,2003，22：881-884.

[23] Valenciano A，Biondi B.2D deconvolution imaging condition for shot-profile migration.SEG Technical Program Expanded Abstracts，2003，22：1059-1062.

[24] Fletcher R，F(xiàn)owler P，Kitchenside P.Supressing unwanted internal reflections in prestack reverse-time migration.Geophysics，2006，71(6)：E79-E82.

[25] Yoon K，Marfurt K J.Reverse-time migration using the poynting vector.Exploration Geophysics，2006，37(1)：102-107.

[26] Guitton A，Kaelin B，Biondi B.Least-square attenuation of reverse-time migration artifacts.SEG Technical Program Expanded Abstracts，2006，25：2348-2352.

[27] Liu F，Zhang G Q，Morton S A.Reverse-time migration using one-way wavefield imaging condition.SEG Technical Program Expanded Abstracts，2007，26：2170-2174.

[28] Zhang Y，Sun J.Practical issues of reverse time mi-gration：true amplitude gathers，noise removal and harmonic-source encoding.ASEG Extended Abstracts，2009，1-5.

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

[30] Rocha D，Tanushev N，Sava P.Acoustic wavefield imaging using the energy norm.Geophysics，2016，81(4)：S151-S163.

[31] Fei T W，Luo Y，Yang J et al.Removing false images in reverse time migration：The concept of de-primary.Geophysics，2015，80(6)：S237-S244.