基于第二類切比雪夫多項式直接展開的疊前深度偏移方法

關鍵詞： 頻率—波數域，波場延拓，第二類切比雪夫展開，疊前深度偏移中圖分類號：P61 文獻標識碼：A DOI：10. 1810/j. cnki. issn. 1000-7210. 2024028

Prestack depth migration method based on direct expansion of the second kind of Chebyshev polynomials

ZHANG Qi1，2，LI Zhenchun1，2，QU Yingming1，2，YANG Jing3，ZHANG Lei3tate Key Laboratory of Deep Oil and Gas，China University of Petroleum（East China），Qingdao，Shandong 266580，China；2. School of Geosciences，China University of Petroleum（East China），Qingdao，Shandong 266580，China；3. Data Processing Center，Geological Research Institute，BGP Inc. ，CNPC，Zhuozhou，Hebei ，China）

Abstract：The second kind of Chebyshev expansion has been widely used in numerical approximations since it has stricter constraints on the objective function compared to the first kind of Chebyshev expansion，resulting in smaller error accumulation. Starting from the frequency?domain acoustic wave equation，this paper introduces the second kind of Chebyshev expansion into the Fourier one ?way wave continuation operator to achieve wave? field continuation. On this basis，with the use of the Born approximation theory and cross ? correlation imaging conditions，prestack depth migration imaging is realized based on the direct expansion of the second kind of Che? byshev polynomials. This method utilizes the accurate velocity at each location during depth propagation，re? ducing errors compared to the conventional one?way wave continuation method which uses the mode of depth ? layer background velocity combined with perturbation velocity. This achieves higher accuracy of the wavefield continuation operator under strong lateral velocity variations. The effectiveness of the method is verified through wavefield snapshot analysis and model imaging tests. It has significantly improved imaging accuracy compared to conventional prestack depth migration methods for one?way wave.

Keywords：frequency?wavenumber domain，wavefield continuation，second kind of Chebyshev expansion， prestack depth migration

張琦，李振春，曲英銘，等 . 基于第二類切比雪夫多項式直接展開的疊前深度偏移方法［J］. 石油地球物理勘探，2025，60（3）：710?717.ZHANG Qi，LI Zhenchun，QU Yingming，et al. Prestack depth migration method based on direct expansion ofthe second kind of Chebyshev polynomials［J］. Oil Geophysical Prospecting，2025，60（3）：710?717.

0 引言

地震偏移成像是在一定的數學物理模型（聲介質，彈性介質等）基礎上，應用相應的地球物理理論，將地面觀測到的多次覆蓋數據反傳，消除地震波的傳播效應，得到地下介質的模型圖像。偏移成像方法按照其地球物理理論的不同大致可以分為兩類：以Kirchhoff 積分法為代表的射線類偏移成像方法，以及以逆時偏移（RTM）、傅里葉有限差分偏移（FFD）等為代表的波動方程類偏移成像方法。Kirchhoff 積分法依賴于地震射線方法計算地震波的旅行時，計算效率很高且對觀測系統具有良好的適應性，但它受限于理論特性，當地下速度橫向變化強烈、構造復雜時，Kirchhoff 積分法偏移成像的準確性會大幅降低，難以對地下構造準確成像。相比之下，以波動方程為基礎的、未做高頻近似假設的偏移成像方法能夠更逼近地震波在地下的真實傳播情況，使地表接收的繞射、散射和反射信息更準確地歸位到真實的構造位置[1]。全波波動方程方法雖然能夠模擬出全波場信息，但在偏移剖面上會出現低頻干擾噪聲，處理大型數據時對內存和算力的需求較高。隨著圖形處理器（GPU）技術的不斷發展，引入 GPU 加速的 RTM 計算效率大幅提高，但是在大型生產資料的處理中仍然面臨諸多計算方面的挑戰[2-3]。單程波偏移成像方法通過沿波傳播的主要方向分解完整的波動方程，與 Kirchhoff積分法等射線類偏移方法相比，避免了對波動方程繁瑣的逼近假設，具有更高的成像精度；與逆時偏移方法相比，單程波動方程計算效率更高，可節省計算機內存空間。發展并完善單程波偏移理論對偏移成像技術的綜合應用具有重要意義。

近年來，學者們基于單程波理論實現了多種單程波疊前深度偏移算法。1971 年 Clearbout[4]首先通過有限差分實現了基于波動理論的偏移方法；Gazdag[5]、Stolt[6]隨后提出了頻率—波數域單程波偏移方法；為應對橫向速度變化劇烈的情況，在頻率—空間域（F-X）或者混合域（F-X、F-K）實現波場外推的成像方法也相繼出現[7-10]。You 等[11]提出一種基于雙檢模式的波場延拓模式，進一步提高了單程波算子的成像角度；Terekhov[12]提出一種三維 Laguerre單程波偏移方法，具有較好應用效果；Song 等[13]推導了切比雪夫直接展開形式的各階算子。單程波偏移相關理論得到了長足的發展和完善[14-19]。

單程波波場外推主要包括兩種方式：一是基于有 限 差 分 （FD） 法[4，20-22]，二 是 采 用 傅 里 葉 方法[5，7-10，23]。基于 FD 的方法雖然可以適應速度的橫向變化，但是存在雙程分裂誤差、數值頻散等問題[24]；基于偽譜法思想的傅里葉方法可以有效避免數值頻散和雙程分裂誤差。傳統的單程波傅里葉方法存在的主要問題是在大傾角和速度變化劇烈的情況下成像效果較差。對此，一種處理方式是將介質近似為均勻速度的背景場加上速度變化產生的擾動場，因此相繼提出了裂步傅里葉（SSF）、傅里葉有限差分（FFD）、廣義屏（GSP）等單程波方法[5-10]。利用這類方法可以得到較為準確的波場，但是精度受速度擾動的影響較大，并且這種影響也無法通過增加高階項來消除。另一處理方式是無需引入速度擾動場，而是通過各類近似手段將波動方程直接展開成沿深度延拓的單程波算子，延拓中直接使用各位置的速度進行波場延拓，這類方法在介質橫向速度劇烈時具有更強的適應能力。

為了對復雜構造區域進行準確成像，利用切比雪夫展開實現的波場延拓算法不斷發展。傳統的切比雪夫單程波延拓算法類似于裂步傅里葉方法，通過時移校正和切比雪夫高階校正適應復雜速度變化區域，但這類方法存在與 SSF、FFD 相同的問題。Song等[13]推導了一種傅里葉單程波傳播算子的切比雪夫直接展開形式，在實現單程波場傳播時不需要引入分層背景速度。本文從頻率域二維聲波波動方程出發，詳細推導了波動方程引入第二類切比雪夫多項式的展開形式，再將其與 De Hoop 等[10]提出的歸一化算子結合，得到第二類 4 階切比雪夫直接展開的單程波延拓算子。算子由兩部分組成：第一部分為第二類切比雪夫展開常數項對應的頻率空間域的初步延拓；第二部分為包含第二類切比雪夫展開高階項與歸一化算子共同構成的波數頻率域的校正項。最后，基于Born 近似理論與互相關成像條件實現了一種新型單程波疊前深度偏移算法，通過波場快照分析和模型測試驗證方法的準確性和優勢。

1 第 二 類 切 比 雪 夫 波 場 延 拓 算 子推導

頻率—波數域中二維聲波方程的表達式為

式中： k_x 是水平方向的波數； P（k_x，z，ω） 是頻率—波數域波場， ω 表示角頻率； 是波場傳播速度。垂直方向波數 k_z 與 k_x 滿足

式中 R=vk_x/ω ，其取值范圍為 [-1，1] 。以每層的參考常速度代入式（1），可以得到相移法單程波延拓算子，但是該算子只適用于均勻層狀介質或者橫向速度變化不大的介質。但是，對于實際情況，波傳播速度 v 往往是隨空間位置變化的。對于存在垂直或者橫向速度變化的情況，一種解決方法是通過添加各個空間位置的速度擾動項，如裂步傅里葉方法和傅里葉有限差分方法；另一種解決方法是通過對式（2）進行近似展開，如泰勒級數展開、Padé 近似展開、切比雪夫展開等。

第二類切比雪夫展開一直廣泛應用于數值近似，相比于第一類切比雪夫展開，對目標函數的約束更為嚴格，具有更小的誤差累積。若函數 f（λ）∈ [-1，1]，連續且有限，函數的第二類切比雪夫多項式形式為

式中 U_n（λ） 是由遞歸關系定義的第二類切比雪夫多項式，具有如下關系和初值

展開式的系數 a_n 為

當 時，系數方程級數展開的奇階系數為零，下文以第二類 2 階切比雪夫展開（ （n=4） ）為例，推導并得到相關的系數，表達式為

求得展開系數后，對應的多項式值為

得到對應的切比雪夫展開系數后，第二類 2 階和 4 階切比雪夫展開單程波波動方程可分別表示為

為了在數值實現中避免溢出，延拓過程中需要對展開式中的高階項進行正則化處理[10]

式中 P 和 q 分別表示復數的實部和虛部。結合式（9）的波動方程展開形式和式（10）的正則化處理，基于

第二類 4 階切比雪夫多項式直接展開的單程波延拓算子數值實現包括兩個步驟：首先，在頻率—空間域開展初步延拓

然后，在頻率—波數域開展基于切比雪夫展開的高階項校正

{F⁺[P₁（x，z+Δz;ω）]+ξ}

式中： P（x，z;ω） 為深度 z 的波場； P₁（x，z+Δz;ω） 為經初步延拓的波場； Δz 為波場延拓的單位深度；Re（?） 和 Im（?） 分別表示求實部和虛部； ζ 是一個很小的數（例如 0. 0001），作用是避免分母為零時出現不穩定性； F⁺ 和 F^- 分別為 x 方向的正、反傅里葉變換； c 為基于第二類 4 階切比雪夫展開的高階校正項； d 為正則化處理算子； P₂（x，z+Δz;ω） 為延拓到下一層深度的波場。

圖1 第一類和第二類切比雪夫2 階和4 階延拓算子的波場快照

為了對第一類和第二類切比雪夫直接展開得到的單程波延拓算子的特征進行比較，分別求取第一類（CF-I-2）和第二類切比雪夫 2 階（CF-II-2）和 4 階延拓算子（CF-I-4，CF-II-4）的波場快照，如圖 1 所示。所用模型為二維均勻速度場，網格數為 501× 251，網格間距為 10m ，震源采用 30Hz 主頻的雷克子波，模型速度為 4000m/s ，左、右象限分別為對應方法波場快照的左、右部分。由圖 1a 可見：第一類切比雪夫延拓算子在 2 階情況下，在 CF-I-2 箭頭處，波場位置離理論波場位置存在較大偏差；4 階情況下，在 CF-I-4 箭頭處，波場存在較大位置偏差和強振幅。由圖 1b 可見：第二類切比雪夫延拓算子的成像角度范圍雖然略小于第一類，但第二類切比雪夫延拓算子的波場與理論波場較第一類更為接近，不同角度的振幅都十分均勻，整體上誤差明顯更低。

2 疊前深度偏移成像

基于切比雪夫直接展開的單程波疊前偏移算法是根據 Born 近似理論建立的。根據 Born 近似理論，頻率域一次反射波可以表示為[24]

式中：下標 r 和 s 分別代表檢波點和震源； V（x） 表示含有散射體的模型空間； U_r（x_r|x_s，ω） 表示在檢波點 x_r 記錄的一次反射波； U_s（x|x_s，ω） 表示由震源位置 σ_xs 傳播到散射點 x 的入射波場； m（x） 是地下位置 x 點的反射系數，即偏移成像的目標； G（x|x_s，ω） 是波從震源位置 ?_xs 傳播到散射點 x 的 Green 函數；G（x_r|x，ω） 是由散射點 x 傳播到接收點 x_r 的 Green函數； f（ω） 是頻率域震源子波。將式（15）表示為矢量矩陣形式

U_r=Lm

式中： m 為反射系數矩陣； L 為正演算子。正演過程即震源子波由地下模型空間 ?_m 向數據空間的投影過程。 L 與偏移算子 M 具體的表達形式如下

3 模型試算

為了驗證本文基于第二類切比雪夫直接展開的傅里葉疊前偏移成像的效果，分別選取雙層模型和橫向速度變化劇烈的 Marmousi 模型進行偏移成像測試。偏移過程使用的合成數據由前文推導的單程波方法正演得到，成像條件為互相關成像。

雙層模型參數見圖 2a 。模型空間剖分網格數為 201×101 ，網格間距為 20m ，炮點位于第 101 個網格點，共 1 炮，震源采用 20Hz 主頻的雷克子波。兩類切比雪夫展開的單炮成像結果見圖 2b 和圖 2c 。由圖可見，兩類切比雪夫展開方法都能準確歸位平層，但第二類切比雪夫展開疊前偏移相比第一類的優勢主要體現在層位整體振幅的均勻性 （ 圖 2b、圖 2c 中箭頭處），這與前文波場快照（圖 1a 中箭頭處）得出的結論是一致的。

Marmousi 模型參數見圖 3。模型空間剖分網格數為 681×381 ，網格間距為 20m ，炮點位于第10～670 個網格點，共 31 炮，炮間距為 440m ，震源采用 30Hz 主頻的雷克子波。每炮有 681 道接收，道間距為 20m ，每道記錄時間為 6 s，采樣間隔為1ms 。

圖 4 為 Marmousi 模型四種單程波延拓算子（SSF、FFD、CF-II-2 和 CF-II-4）得到的疊前偏移圖像。由圖可見，這四種方法都能夠較準確地歸位平層區域，但是在模型中央速度變化劇烈的區域成像效果差別較大。基于第二類切比雪夫展開的兩種疊前偏移方法（圖 4c、圖 4d）得到的反射軸能量更加均勻、層位更加清晰準確，尤其是復雜區域的成像效果明顯優于 SSF 和 FFD 的偏移方法（圖 4a、圖 4b）。CF-II-4 雖然在部分區域反射軸的能量低于 CF-II-2，但是整體能量收斂更好，反射軸的分辨率更高，且在復雜層位區域不存在 FFD 和 CF-II-2偏移圖像中的噪聲干擾（圖 4b、圖 4c），因此 CF-II-4的綜合成像效果是這四類方法中最優的。當然，由于傅里葉類偏移方法的計算量主要集中于各類域變換（時間—頻率域，空間—波數域）過程中的傅里葉變換，從式（12）可知，切比雪夫展開階數越高則需要進行域變換的項也越多，計算量也會相應增大。

圖2 雙層模型及兩類切比雪夫展開的單炮成像結果

圖 3 Marmousi 速度模型

為了更清晰地對比前述四種方法的成像效果，對 Marmousi 模型中速度變化劇烈的兩個區域進行局部放大（圖 3 中的區域 A 和區域 B），結果見圖 5和圖 6。在區域 A，中央位置的高陡地層、右上區域的復雜平層及底部的層位交錯位置是三個主要的成像難點。由圖 5a 可見，SSF 延拓算子對中央位置的高陡地層和右上角位置的復雜平層成像十分模糊；

圖4 Marmousi 模型不同單程波方法成像結果

圖5 Marmousi 模型區域A 不同方法偏移結果

由圖 5b 可見，FFD 延拓算子對右上角位置的復雜平層成像效果較差；在底部層位交錯位置，FFD 和CF-II-2 延拓算子均存在因反射波未準確歸位引起的噪聲干擾（圖 5b、圖 5c），出現層位信息失準現象；由圖 5d 可見，CF-II-4 延拓算子對應的偏移剖面上，區域 A 中三個成像難點部分成像清晰，反射波準確歸位，更好地反映了地下構造信息。此外，在一些弱反射層上，CF-II-2 的成像軸能量高于 CF-II-4，這是由于在部分相位角范圍內 2 階算子的精度高于 4 階算子[11]。

圖 6 展示的是細小傾斜層位密集排列的區域 B的成像剖面。可以進一步看出 CF 延拓算子對復雜速度層位的適應力比 SSF 和 FFD 更強，對區域中傾斜層位的刻畫也更為準確。在區域 B，CF-II-2 成像中同相軸比 CF-II-4 更為清晰、能量更強，CF-II-2整體的成像效果更好，但是 CF-II-4 對部分層位具有更高的分辨率。可以發現，區域中部分較薄層位置的上、下邊界距離較近，在 CF-II-4 成像剖面上得到較好的歸位，但在 CF-II-2 算子成像剖面上，上、下界面混疊在一起，難以區別（圖 6c、圖 6d 中箭頭位置）。

圖 6 Marmousi 模型區域 B 偏移結果

4 結論

本文推導了第二類切比雪夫直接展開的單程波延拓算子，基于 Born 近似理論實現了一種單程波疊前偏移成像新方法。通過波場分析和模型試算，得到了如下結論。

（1）CF-I、CF-II 波場延拓算子的波場傳播情況和單層模型成像結果對比結果表明，第二類切比雪夫傅里葉方法的整體精度超過對應階數的第一類切比雪夫方法。

（2）Marmousi 模型疊前偏移結果表明 CF-II 算子對復雜區域具有良好的適應能力，相較于傳統的SSF、FFD 偏移方法，該方法能夠得到更好的偏移成像結果。在實際應用中，可以綜合考慮成像需求和計算成本，選擇使用 2 階或更高階的 CF-II 偏移方法。

因此，基于第二類切比雪夫多項式直接展開的單程波疊前偏移對橫向速度變化劇烈的復雜介質具有更高的成像精度，能克服數值頻散和雙向分裂誤差，具有較高的計算效率。

參 考 文 獻

［1］ 李振春. 地震偏移成像技術研究現狀與發展趨勢[J]. 石油地球物理勘探， 2014， 49（1）： 1?21. LI Zhenchun. Research status and development trends for seismic migration technology[J]. Oil Geophysical Prospecting， 2014， 49（1）： 1?21.

［2］ 段心標 . 應用 GPU 的傅里葉有限差分逆時偏移[J]. 石油地球物理勘探， 2020， 55（5）： 1039?1046. DUAN Xinbiao. Fourier finite ? difference reverse time migration using GPU[J]. Oil Geophysical Prospecting， 2020， 55（5）： 1039?1046.

［3］ 楊禮勝， 黃金強， 高國超， 等. 應用球型差分模板的 低秩有限差分法純qP 波逆時偏移[J]. 石油地球物理 勘探， 2023， 58（5）： 1101?1114. YANG Lisheng， HUANG Jinqiang， GAO Guochao， et al. Pure qP ? wave reverse time migration using low rank finite difference method with spherical difference stencil[J]. Oil Geophysical Prospecting， 2023， 58（5）： 1101?1114.

［4］ CLEARBOUT J F. Toward a unified theory of reflec? tor mapping[J]. Geophysics， 1971， 36（3）： 467?481.

［5］ GAZDAG J. Wave equation migration with the phase ? shift method[J]. Geophysics， 1978， 43（7）： 1342?1351.

［6］ STOLT R H. Migration by Fourier transform[J]. Geophysics， 1978， 43（1）： 23?48.

［7］ STOFFA P L， FOKKEMA T J T， DE LUNA FREIRE R M， et al. Split ? step Fourier migration[J]. Geophysics， 1990， 55（4）： 410?421.

［ 8 ］ RISTOW D， RUHL T. Fourier finite ? difference mi? gration[J]. Geophysics， 1994， 59（12）： 1882?1893.

［9］ WU R， DE HOOP M V. Accuracy analysis and nu? merical tests of screen propagators for wave extrapola? tion[C]. Proceedings of SPIE， 1996： 196?209.

［10］ DE HOOP M V， LE ROUSSEAU J H， WU R. Generalization of the phase ? screen approximation for the scattering of acoustic wave[J]. Wave Motion， 2000， 32（1）： 43?70.

［11］ YOU J， WU R， LIU X， et al. Two?way wave equa? tion?based depth migration using one?way propagators on a bilayer sensor seismic acquisition system[J]. Geo? physics， 2018， 83（3）： S271?S278.

［12］ TEREKHOV A V. A three?dimensional Laguerre one? way wave equation solver[J]. Applied Numerical Mathematics， 2022， 173（3）： 380?394.

［13］ SONG H， ZHANG J， ZOU Y. Direct expansion of Fourier extrapolator for one ? way wave equation using Chebyshev polynomials of the second kind[J]. Geo? physics， 2022， 87（2）： S63?S73.

［14］ LAKOBA T I. Study of instability of the Fourier split? step method for the massive Gross?Neveu model[J]. Journal of Computational Physics， 2020， 402（1）： 109100.

［15］ TANG C， FU L， PAN W， et al. Optimized pseudo? Padé Fourier migrator in terms of propagation angles [J]. IEEE Access， 2020， 8： 32054?32065.

［16］ YU Y， LI Y， WU X， et al. Enhancing one?way wave equation?based migration with deep learning[J]. Geo? physics， 2022， 88（1）： WA105?WA114.

［17］ ESLAMINA M， ELMELIEGY A M， GUDDATI M N. Improved least?squares migration through double? sweeping solver[J]. Geophysics， 2023， 88（3）： S131 ? S141.

［18］ YOU J， PAN N， LIU W， et al. Efficient wavefield separation by reformulation of two ? way wave ? equation depth?extrapolation scheme[J]. Geophysics， 2022， 87（4）： S209?S222.

［19］ YOU J， LIU W， HUANG X， et al. Q ?compensated wavefield depth extrapolation ? based migration using a vicoacoustic wave equation[J]. Geophysics， 2023， 89（1）： S1?S14.

［20］ WANG Y. ADI plus interpolation： Accurate finite?dif? ference solution to 3D paraxial wave equation[J]. Geo? physical Prospecting， 2001， 49（5）： 547?556.

［21］ TEREKHOV A V. The Laguerre finite difference one? way equation solver[J]. Computer Physics Communica? tions， 2017， 214： 71?82.

［22］ TEREKHOV A V. The stabilization of high?order multistep schemes for the Laguerre one?way wave equa? tion solver[J]. Journal of Computational Physics， 2018， 368（1）： 115?130.

［23］ WU R. Wide?angle elastic wave one?way propagation in heterogeneous media and an elastic wave complex ? screen method[J]. Journal of Geophysical Research， 1994， 99： 751?766.

［24］ 朱峰 . VTI 介質單程波最小二乘偏移方法研究[D]. 青島： 中國石油大學（華東）， 2017. ZHU Feng. The Research of Least ?squares One ?way Wave Migration in VTI Media[D]. Qingdao： China University of Petroleum （East China）， 2017.

（本文編輯：劉海櫻）

作 者 簡 介

張琦 博士研究生，1998 年生；2021年獲中國石油大學（華東）勘查技術與工程專業學士學位；現在中國石油大學（華東）攻讀地質資源與地質工程專業博士學位，主要研究方向是單程波疊前偏移和高性能地震偏移成像技術。