單程和雙程波動方程疊前深度偏移方法

田東升，王云專，李義鵬，石 穎，柯 璇，李婷婷，劉淑芬

（1.東北石油大學 地球科學學院，黑龍江 大慶 163318； 2.北華航天工業(yè)學院 電子工程系，河北 廊坊 065000）

單程和雙程波動方程疊前深度偏移方法

田東升1，王云專1，李義鵬2，石 穎1，柯 璇1，李婷婷1，劉淑芬1

（1.東北石油大學 地球科學學院，黑龍江 大慶 163318； 2.北華航天工業(yè)學院 電子工程系，河北 廊坊 065000）

疊前深度偏移是獲得地下構造映像的有效手段，而基于波動方程的疊前深度偏移方法對速度橫向變化劇烈的地層有更好的適應性.分析基于單程波方程的相移法、相移加插值法、頻率—空間域有限差分法、傅里葉有限差分法和基于雙程波方程的逆時偏移方法，借助于地塹模型與鹽丘模型，測試5種逆時偏移方法成像復雜構造的精度和適應性.結果表明，基于波動方程的疊前深度偏移方法可實現(xiàn)橫向變速地下構造成像，相比于基于雙程波方程的逆時偏移方法，單程波方程方法對垂直斷層等高陡傾角構造成像有局限性；逆時偏移方法對垂直斷層、鹽丘下邊界等復雜構造可以清晰成像，輔以精確的地層速度，逆時偏移方法在地震資料成像領域中有廣闊的發(fā)展和應用前景.

疊前深度偏移；相移法；相移加插值；頻率—空間域有限差分；傅里葉有限差分；逆時偏移

0 引言

隨著油氣勘探目標日趨復雜，具體表現(xiàn)在斷塊小、傾角陡、縱橫向速度變化劇烈等方面，常規(guī)偏移方法很難達到地震資料處理要求，加之巨大的數(shù)據(jù)處理量，亟需研究高精度、高效率的偏移算法，為精細地質(zhì)構造解釋及儲層識別提供重要依據(jù).疊前深度偏移是獲得精確地下構造的有效途徑.在數(shù)學解法上，疊前深度偏移分為兩類：基于波動方程積分解的射線偏移和基于波動方程微分解的波動方程疊前深度偏移.波動方程疊前深度偏移較射線偏移在處理橫向變速問題上具有更強的適應性，不存在射線偏移法中成像點的多值走時問題.其中波動方程疊前深度偏移分為單程波法和雙程波法.

20世紀70年代，Claerbout J F[1]提出應用波動方程進行偏移，采用有限差分法求解得到單程波動方程15°近似公式，該方法在主傳播方向小范圍內(nèi)具有較好的效果，對寬方位地震波傳播模擬并不理想，尤其對于地震數(shù)據(jù)中包含的水平和陡傾角反射信息的偏移效果不明顯.Gazdag J最初提出的頻率—波數(shù)域相移法[2]僅適用于對地下橫向速度不發(fā)生變化的介質(zhì)進行成像；進而提出的相移加插值方法[3]可以適應存在橫向速度變化的介質(zhì)，但需要頻繁計算參考波場，效率較低[4].隨后，人們將有限差分法與傅里葉法結合研究，提出混合偏移方法[5-6]，對橫向速度變化強烈的介質(zhì)成像效果良好.王玉學等[7]推導上行波方程的兩種高階近似表達式.馮鳳萍等[8]研究加吸收層的三維45°上行波方程的隱式差分格式.然而，所有單程波偏移方法存在局限性，即單程波算子在成像大角度傳播的波時將發(fā)生相位改變和振幅減弱的現(xiàn)象，無法對陡傾角進行成像；另外，單程波法也無法成像回轉(zhuǎn)波.

采用雙程波動方程進行偏移能很好地適應劇烈的橫向速度變化，可以有效解決復雜地質(zhì)體成像問題，最典型的方法就是逆時偏移方法.該方法最早由Whitmore D N等[9]提出，最初應用于處理疊后資料.逆時偏移方法簡單、易于實現(xiàn)，不對波動方程做任何近似，從而不存在傾角的限制，可以對透射波、多次波、繞射波等進行成像.近些年，隨著計算機硬件技術的迅猛發(fā)展及勘探要求的日益提高，逆時偏移方法的研究也從疊后走向疊前[10-11]，從二維走向三維.GPU加速計算技術的引入[12]和噪音壓制策略的研究[13]推動逆時偏移技術的發(fā)展.同時，如全波形反演的精準速度建模方法的研究[14]也加速波動方程疊前深度偏移的研究進程.

筆者首先闡述單程和雙程波動方程疊前深度偏移方法的基本原理，分析不同方法優(yōu)缺點，通過地塹模型和二維鹽丘模型進行成像測試，分析不同方法成像復雜構造的精度和適應性，為高精度疊前深度偏移方法的工業(yè)化應用提供依據(jù).

1 單程波方程疊前深度偏移方法

1.1 相移和相移加插值法

在縱波勘探中，假設地下介質(zhì)為均勻各向同性介質(zhì)，并且介質(zhì)密度恒定，可以用聲波方程描述地震波的傳播，二維形式為

式中：v為地下介質(zhì)速度；W 為t時刻（x，z）位置處的波場值；t為時間；x，z為空間方向.式（1）在頻率—波數(shù)域的解析解為

式中：kx為水平波數(shù)；kz為垂直波數(shù)；d z為深度延拓步長；w為頻率.

通過標量波動方程的頻散關系得到上、下行波分解的形式表達式，將它代入式（2）并求偏導數(shù)，得到頻率—波數(shù)域單程波方程為

式（3）為相移法的基本原理.相移法的優(yōu)點在于穩(wěn)定性，對網(wǎng)格間距沒有要求，但是僅能對橫向速度無變化的介質(zhì)進行偏移；因此也決定它無法對復雜地質(zhì)體準確成像.

為了能更好地適應橫向速度變化，引入相移加插值法，將式（2）分解成2個獨立的公式：

式中：v′為v 的近似，v′≠v（x，z）.

利用式（4）對每道進行時移；再根據(jù)式（5）使用v（x，z）的最小值和最大值進行相移；最后對兩個相移波場的模數(shù)和相位角進行線性插值，得到最終的波場信息.相移加插值法較相移法有很大的改進，可對橫向速度變化緩慢的介質(zhì)進行成像[4].

1.2 頻率—空間域有限差分法

對式（3）二階近似后分離成兩部分，可得頻率—空間域有限差分法表達式：

式（6）為繞射項方程，可以通過有限差分求解，實現(xiàn)對繞射波的收斂.式（7）可看作折射項方程，在頻率—空間域表示為

式中：Dxx為x方向二階導數(shù)；Dz為z方向一階導數(shù).

通過有限差分法或相移法計算式（8），可以完成橫向速度變化的時差校正[15].對偏移傾角不大的地質(zhì)構造，該方法成像效果較理想.

1.3 混合偏移法

混合偏移法是將有限差分法與傅里葉法相結合，主要包括裂步傅里葉法和傅里葉有限差分法.裂步傅里葉法基于擾動理論，將速度場分離成背景速度項和擾動項.該方法首先在頻率—波數(shù)域以背景速度延拓，而后在頻率—空間域利用表示速度橫向變化的擾動項進行相移；在成像橫向速度變化強烈的介質(zhì)時，成像效果不理想.

傅里葉有限差分法以裂步傅里葉法為基礎，在波場延拓過程中，引入自適應有限差分算子.由于在選取偏移速度時，采用介質(zhì)速度v和參考速度vr偏移得到的結果不一致，根據(jù)二者偏移誤差d，可得傅里葉有限差分法表達式：

第Ⅰ部分是在頻率—波數(shù)域進行的相位移算子；第Ⅱ部分是裂步傅里葉法在頻率—空間域進行的相位移算子；第Ⅲ部分為有限差分算子.當?shù)叵陆橘|(zhì)為層狀介質(zhì)時，式（9）只保留第Ⅰ部分，該算法變?yōu)橄嘁品?；當?shù)叵陆橘|(zhì)速度橫向變化劇烈時，該算法變?yōu)橛邢薏罘址╗16].因此，傅里葉有限差分法是裂步傅里葉法和頻率—空間域有限差分法的混合偏移方法，兼顧二者各自的優(yōu)點.

2 雙程波方程疊前深度偏移方法

相對于單程波方程偏移算法，通過求解雙程波方程得到波場傳播信息的典型方法是疊前深度逆時偏移方法，實現(xiàn)步驟主要包括：（1）震源處波場沿時間軸正向延拓，保存各個時刻波場信息；（2）檢波點處波場沿時間軸反向延拓；（3）對同一時刻的兩個波場進行成像，完成單炮逆時深度偏移.對各炮的成像結果進行疊加，得到疊前逆時偏移結果.

在逆時偏移計算中，對式（1）進行高階有限差分，空間方向2N階差分格式[17]為

時間方向差分格式為

逆時偏移的成像條件主要有激發(fā)時刻成像、互相關成像和除法成像.目前最常用的成像條件是互相關成像條件，即

式中：s（x，z，t）、r（x，z，t）分別為震源波場與檢波點波場.

采用Robert G Clapp提出的隨機邊界條件[18]，無需存儲波場正傳過程中所有時刻的波場，節(jié)省大量的存儲空間；利用GPU/CPU協(xié)同并行技術加速逆時偏移計算[19]，提高算法的計算效率；由于逆時偏移成像結果中含有大量低頻噪音，通過拉普拉斯算子法對噪音進行壓制[20-21].不同波動方程疊前深度偏移方法優(yōu)缺點見表1.

表1 波動方程疊前深度偏移方法優(yōu)缺點Table 1 Comparison of different wave equation prestack depth migration approaches

3 計算精度分析

為測試單程和雙程波動方程疊前深度偏移方法對復雜構造的偏移效果，分別成像地塹模型和二維鹽丘模型，分析不同算法的成像精度和適應性.

3.1 地塹模型

地塹模型網(wǎng)格點數(shù)為200×100，網(wǎng)格大小為5 m×5 m，上層介質(zhì)速度為2 000 m/s，下層速度為3 000 m/s，激發(fā)震源采用Ricker子波，主頻為40 Hz，震源與檢波器始于地面最左端，炮間距20 m，共50炮，每炮200道接收，道間距為5 m.

地塹速度模型見圖1（a）；采用頻率—空間域有限差分算法得到的偏移結果見圖1（b）、（c），其中圖1（c）將正演中的棱柱波切除.由圖1（b）、（c）可知，棱柱波對成像結果存在影響（圖中方框位置）.另外，成像剖面中上層存在較大頻散，是由算法本身造成的（圖中箭頭指示位置）.采用傅里葉有限差分算法得到的偏移結果見圖1（d）、（e），其中圖1（e）將正演中的棱柱波切除.由圖1（d）、（e）可知，傅里葉有限差分法對棱柱波不能準確歸位，因此對垂直斷層成像效果不理想.采用裂步傅里葉法得到的偏移結果見圖1（f）.由圖1（f）可知，下部同相軸不清晰，邊界無法識別.基于單程波方程的3種算法對水平分層可有效成像，但無法成像垂直斷層.

采用逆時偏移方法經(jīng)拉普拉斯算子法壓制低頻噪音后得到的偏移結果見圖1（g）.由圖1（g）可知，成像剖面清晰，水平分層明顯，垂直斷層得到很好歸位（圖中橢圓位置）.相對于單程波算法，逆時偏移方法可以對高陡傾角準確成像，即對橫向速度變化劇烈的介質(zhì)成像效果較好.

圖1 不同疊前深度偏移算法成像地塹模型Fig.1 Graben model imaging result by different prestack depth migration approaches

3.2 二維鹽丘模型

二維鹽丘模型網(wǎng)格點數(shù)為150×649，網(wǎng)格大小為24.384 m×24.384 m；震源采用Ricker子波，主頻為18 Hz；共325炮，每炮176道接收，炮間距為48.768 m.

二維鹽丘速度模型見圖2（a），采用有限差分算法得到的偏移結果見圖2（b）.由圖2（b）可知，鹽丘模型整體成像不清晰，上部斷層不明顯，上邊界較模糊，下邊界幾乎無法識別.采用相移加插值法得到的偏移結果見圖2（c），采用傅里葉有限差分算法得到的偏移結果見圖2（d）.由圖2（c）、（d）可知，兩者對橫向變速的地質(zhì)體成像有明顯改善，上部斷層成像不清晰，下部邊界不準確.總體上，在單程波法偏移中，相移加插值算法成像精度較高，但耗時最長.采用疊前逆時偏移方法，經(jīng)拉普拉斯算子法壓制低頻噪音后得到的成像結果見圖2（e）.由圖2（e）可知，與單程波偏移結果相比，逆時偏移得到的二維鹽丘剖面成像清晰，上部斷層得到很好歸位（圖中方框位置），鹽丘下邊界同相軸明顯（圖中橢圓位置），下部陡傾角得到清晰成像，成像剖面質(zhì)量良好.

圖2 不同疊前深度偏移算法成像結果Fig.2 Imaging result by different prestack depth migration approaches

4 結論

（1）基于波動方程的疊前深度偏移方法不受高頻近似和多值走時的影響，對速度橫向變化劇烈的地層有很好的適應性.

（2）單程波方程方法在成像大角度傳播的地震波場時，存在相位改變和振幅削弱的問題；另外，因無法成像回轉(zhuǎn)波，導致難以對陡傾角構造成像.

（3）基于雙程波的逆時偏移法可有效成像各種地震波，利用精準的地層速度，對橫向變速地層和高陡構造地層均能高精度成像，但是存儲尤其是3D數(shù)據(jù)存儲問題還有待于解決.

[1] Claerbout J F.Toward a unified theory of reflector mapping[J].Geophysics，1971，36（3）：467-481.

[2] Gazdag J.Wave equation migration with the phase-shift method[J].Geophysics，1978，43（7）：1342-1351.

[3] Gazdag J，Sguazzero P.Migration of seismic data by phase shift plus interpolation[J].Geophysics，1984，49（2）：124-131.

[4] 張釙，李幼銘，劉洪.幾類疊前深度偏移方法的研究現(xiàn)狀[J].地球物理學進展，2000，15（2）：30-39.

Zhang Po，Li Youming，Liu Hong.The situation of several prestack depth migration methods[J].Progress in Geophysics，2000，15（2）：30-39.

[5] Stoffa P L，F(xiàn)okkema J T.Split-step Fourier migration[J].Geophysics，1990，55：410-421.

[6] Ristow D，Ruhl T.Fourier finite-difference migration[J].Geophysics，1994，59（12）：1882-1893.

[7] 王玉學，周瑞芬，張廷全.三維上行波的高階近似[J].大慶石油學院學報，2003，27（2）：120-122.

Wang Yuxue，Zhou Ruifen，Zhang Tingquan.High order approximation of 3D one way wave equation[J].Journal of Daqing Petroleum Institute，2003，27（2）：120-122.

[8] 馮鳳萍，周瑞芬，劉暢.加吸收層的三維45°上行波方程的隱式差分格式[J].大慶石油學院學報，2005，29（4）：98-100.

Feng Fengping，Zhou Ruifen，Liu Chang.Implicit difference scheme of 45 degree up-going wave equation with an absorbing layer[J].Journal of Daqing Petroleum Institute，2005，29（4）：98-100.

[9] Whitmore D N.Iterative depth imaging by back time propagation[C].53rd Annual International Meeting，SEG Expanded Abstracts，1983：382-385.

[10] Yoon K，Marfurt K J，Starr W.Challenges in reverse-time migration[C].74th Annual International Meeting，SEG Expanded Abstracts，2004，23（1）：1057-1060.

[11] 劉紅偉，李博，劉洪，等.地震疊前逆時偏移高階有限差分算法及GPU實現(xiàn)[J].地球物理學報，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[J].Chinese Journal of Geophysics，2010，53（7）：1725-1733.

[12] 李博，劉紅偉，劉國峰，等.地震疊前逆時偏移算法的CPU/GPU實施對策[J].地球物理學報，2010，53（12）：2938-2943.

Li Bo，Liu Hongwei，Liu Guofeng，et al.Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J].Chinese Journal of Geophysics，2010，53（12）：2938-2943.

[13] Zhang Yu，Sun J.Practical issues of reverse time migration：True amplitude gathers，noise removal and harmonic-source encoding[J].First Break，2009，26，29-35.

[14] 楊午陽，王西文，雍學善，等.地震全波形反演方法綜述[J].地球物理學進展，2013，28（2）：0766-0776.

Yang Wuyang，Wang Xiwen，Yong Xueshan，et al.The review of seismic full waveform inversion method[J].Progress in Geophysics，2013，28（2）：0766-0776.

[15] 馬淑芳，李振春.波動方程疊前深度偏移方法綜述[J].勘探地球物理進展，2007，30（7）：153-161.

Ma Shufang，Li Zhenchun.Review of wave equation prestack depth migration methods[J].Progress in Exploration Geophysics，2007，30（7）：153-161.

[16] 鄭曉.起伏地表單程波疊前深度偏移方法研究[D].北京：中國地質(zhì)大學，2011，34-43.

Zheng Xiao.Method study of rugged topography prestack one-way equation depth migration[D].Beijing：China Unversity of Geosciences，2011，34-43.

[17] 石穎，柯璇，田東升，等.復雜構造地震數(shù)據(jù)疊前逆時偏移方法[J].數(shù)學的實踐與認識，2013，43（10）：206-213.

Shi Ying，Ke Xuan，Tian Dongsheng，et al.Seismic data prestack reverse time migraion approach on complicated construction[J].Mathematics in Practice and Theory，2013，43（10）：206-213.

[18] Robert G Clapp.Reverse time migration with random boundaries[C].79th Annual International Meeting，SEG Expanded Abstracts，2009：2809-2813.

[19] 石穎，陸加敏，柯璇，等.基于GPU并行加速的疊前逆時偏移方法[J].東北石油大學學報，2012，36（4）：111-115.

Shi Ying，Lu Jiamin，Ke Xuan，et al.Prestack reverse time migration based on GPU parallel accelerating algorithm[J].Journal of Northeast Petroleum University，2012，36（4）：111-115.

[20] 劉紅偉，劉洪，鄒振.地震疊前逆時偏移中的去噪與存儲[J].地球物理學報，2010，53（9）：2171-2180.

Liu Hongwei，Liu Hong，Zou Zhen.The problems of denoise and storage in seismic reverse time migration[J].Chinese Journal of Geophysics，2010，53（9）：2171-2180.

[21] 石穎，王建民，田東升，等.應用于疊前逆時偏移的拉普拉斯算子去噪法[J].大慶石油地質(zhì)與開發(fā)，2012，31（5）：154-157.

Shi Ying，Wang Jianmin，Tian Dongsheng，et al.Denoise approach by Laplace operator applied in prestack reverse time migration[J].Petroleum Geology and Oilfield Development in Daqing，2012，31（5）：154-157.

TE132.1

A

2095-4107（2014）04-0039-06

2013-02-12；

任志平

黑龍江省教育廳科學技術研究項目（12511025）

田東升（1989-），男，碩士研究生，主要從事地震資料處理方面的研究.通訊作者：石 穎，E-mail：shiyingdqpi@163.com

DOI 10.3969/j.issn.2095-4107.2014.04.006