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

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

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

1 引言

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

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

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

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

2 RTM中回轉波假象成因分析

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

(1)

式中：x為地下成像點位置；t為炮檢點傳播到成像點的時間；s(t,x)和r(t,x)分別為時、空域炮檢點波場；Tmax為最大記錄時間；I(x)為成像結果。RTM的實現是基于雙程波動算子，故在波場外推時會產生反射波和透射波。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)表示震源波場下行波與檢波點波場上行波互相關成像結果；I2(x)表示震源波場上行波與檢波點波場下行波互相關成像結果；I3(x)表示震源波場下行波與檢波點波場下行波互相關成像結果；I4(x)表示震源波場上行波與檢波點波場上行波互相關成像結果。在當前反射波偏移成像框架下,I1(x)和I2(x)對應震源波場和檢波點波場反向傳播的波場進行互相關，都可對地下反射界面進行成像。具體而言：I1(x)類似于單程波偏移的成像條件，結果即對應一次反射波成像；I2(x)即對應能量較一次波弱的多次波成像；I3(x)和I4(x)對應震源波場和檢波點波場同向傳播的波場進行互相關，即表現為低頻假象，一般通過拉普拉斯濾波可較徹底地消除。當震源波場和檢波點波場通過強速度梯度或復雜反射界面時，將產生向上回傳的波場(圖1)，此時I2(x)則表現為一次回轉波產生的假象，其能量強于多次波成像，且難以被去除。

圖1 回轉波假象產生示意圖

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

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

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

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

3 用波長平滑算子去除回轉波假象

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

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

(5)

(6)

式中：t表示一個波長內計算網格點到加權點的直線旅行時；T表示對應的周期；p為速度場的維度。

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

為了進一步分析高斯平滑與波長相關平滑對波

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

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

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

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

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

4 數值算例

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

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

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

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

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

5 結論與討論

在地下強速度梯度區域，基于互相關成像條件的RTM會產生回轉波假象，該假象區別于典型的低頻噪聲，有較強的能量且不容易去除。本文通過波場分解成像條件分析回轉波假象產生原因，認為該假象是緣于震源波場上行波與檢波點波場下行波的互相關，即強梯度帶的回轉波在傳統相關成像條件下產生。據此提出一種基于波長平滑算子的回轉波假象去除方法：通過波長相關平滑算子進行速度平滑，在波傳播過程中壓制回轉波的產生，進而消除回轉波假象。該方法區別于通過波場分解法消除回轉波成像假象，可以很好地對高陡構造進行成像且不增加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] 楊午陽，張厚柱，撒利明等. 逆時偏移關鍵問題探討. 石油地球物理勘探，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] 王鵬飛，何兵壽.基于行波分離的三維彈性波矢量場點積互相關成像條件.石油地球物理勘探，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] 葉月明，郭慶新，孫小東等.多次波對逆時偏移構造成像的影響.石油地球物理勘探，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實現.地球物理學報，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 協同三維疊前逆時偏移方法研究與應用.石油地球物理勘探，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，Faria 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，Fowler 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.