高斯射線束疊前深度偏移成像研究

謝 飛，李 佩，黃中玉，3，魏修成，3，朱成宏

（1.中國(guó)石油化工股份有限公司石油勘探開(kāi)發(fā)研究院，北京 100083；2.中國(guó)地質(zhì)大學(xué)（北京），北京 100083；3.中國(guó)石油化工集團(tuán)公司多波地震重點(diǎn)實(shí)驗(yàn)室，北京 100083）

在復(fù)雜介質(zhì)條件下，疊前深度偏移技術(shù)能提供很好的構(gòu)造成像結(jié)果。Kirchhoff積分疊前深度偏移方法以其靈活的目標(biāo)處理功能，高效的數(shù)值計(jì)算和對(duì)空間非規(guī)則采樣數(shù)據(jù)的適應(yīng)能力等特點(diǎn)，在石油工業(yè)界得到了廣泛應(yīng)用。然而，對(duì)于由復(fù)雜構(gòu)造引起的多值走時(shí)、陰影、焦散等問(wèn)題，采用Kirchhoff偏移方法成像效果往往不如單程波方法。高斯束偏移方法是Kirchhoff偏移方法的改進(jìn)，既具備Kirchhoff偏移方法的靈活性和適應(yīng)性，又能解決焦散、多值走時(shí)等問(wèn)題，其成像效果堪比單程波算子，同時(shí)還能克服單程波算子的傾角限制，成像回轉(zhuǎn)波［1－6］。

Hill［7］提出了疊后高斯束偏移方法，通過(guò)局部的傾斜疊加將波場(chǎng)分解為高斯束。對(duì)每個(gè)高斯波束偏移成像、疊加得到最終的偏移結(jié)果。Hale［8－9］詳細(xì)分析了Kirchhoff偏移、傾斜疊加偏移和高斯束偏移的關(guān)系，指出高斯束偏移與Kirchhoff偏移、傾斜疊加偏移相比的優(yōu)點(diǎn)，并給出了具體的計(jì)算公式，分析了其計(jì)算效率。

Hill［10］給出了疊前高斯束偏移的方法。在共偏移距偏移中，通過(guò)對(duì)所有偏移距射線參數(shù)進(jìn)行掃描找到使總走時(shí)虛部最小的炮點(diǎn)射線參數(shù)和檢波點(diǎn)射線參數(shù)的組合。這種掃描可以在稀疏的偏移成像點(diǎn)網(wǎng)格上執(zhí)行，然后將值再內(nèi)插到較密集的成像網(wǎng)格點(diǎn)上，因此具有較高的計(jì)算效率和精度。

Gray［11］指出由Hill提出的共偏移距偏移方法對(duì)于較規(guī)則的觀測(cè)系統(tǒng)是合適的，但是對(duì)于其它觀測(cè)系統(tǒng)的記錄，或者是地表高程變化大、近地表速度變化劇烈時(shí)不利于成像。Gray提出炮集高斯束偏移更適合陸地起伏地表的地震數(shù)據(jù)。

1 高斯束偏移方法原理

假設(shè)地震波場(chǎng)滿足標(biāo)量波動(dòng)方程：

則可將由r′點(diǎn)源激發(fā)傳播到r點(diǎn)的地震波場(chǎng)表示為瑞雷積分形式：

Hill［10］的研究表明，點(diǎn)源波場(chǎng)的格林函數(shù)可用高斯射線束的積分形式表示：

其中，uGB（r，r′，p′；ω）為在r′點(diǎn)激發(fā)、方向?yàn)閜′的射線束。因此，當(dāng)將（2）式中的格林函數(shù)用（3）式的高斯束形式表示時(shí)，就可以得到利用高斯束表示的地震波場(chǎng)。

當(dāng)射線束的中心不在r′點(diǎn)而在鄰近的r0點(diǎn)時(shí)，由于在高斯束發(fā)射的初始位置波束是特定傳播方向的平面波波場(chǎng)，Hill［10］引入相移因子exp［－iωp′（r′－r0）］，則r′點(diǎn)到r點(diǎn)的格林函數(shù)表示為

在克?；舴蚍e分公式（2）中使用格林函數(shù)（4）時(shí)，應(yīng)該將公式（2）中的積分范圍劃分為小的空間積分窗，并且保證在計(jì)算格林函數(shù)（4）時(shí)，空間積分窗的中心r0點(diǎn)與積分窗內(nèi)的任意點(diǎn)r′相距不遠(yuǎn)，能保持初始高斯束的平面波性質(zhì)，為此，Hill［10］引入了如下高斯函數(shù)：

其中，ωl是參考頻率是射線束寬度；a是網(wǎng)格點(diǎn)最小間距。并得到了地震數(shù)據(jù)局部?jī)A斜疊加的頻率域表達(dá)形式：

利用Claerbout［12］的相關(guān)成像原理，將炮點(diǎn)波場(chǎng)和檢波點(diǎn)波場(chǎng)均按照克希霍夫積分公式向下延拓，并將兩個(gè)波場(chǎng)進(jìn)行互相關(guān)可得到偏移成像公式：

將偏移成像公式（7）中的格林函數(shù)用（4）式的高斯束表示，并結(jié)合高斯函數(shù)（5）式，就可以得到以高斯束形式表達(dá)的偏移成像公式：

其中，C0是與高斯函數(shù)相關(guān)的系數(shù)。

Hale［8］對(duì)克?；舴?、傾斜疊加和高斯束3種偏移方法進(jìn)行了詳細(xì)比較，揭示了三者間的聯(lián)系與區(qū)別，認(rèn)為在時(shí)間域做高斯束偏移要比在頻率域效率高，并給出了單個(gè)高斯束對(duì)地下成像貢獻(xiàn)的時(shí)間域表達(dá)式：

其中，AR和AI分別對(duì)應(yīng)復(fù)值振幅的實(shí)部和虛部；τR和τI分別為復(fù)值走時(shí)的實(shí)部和虛部。復(fù)值振幅和復(fù)值走時(shí)均可通過(guò)高斯束射線追蹤得到。px為射線參數(shù)，代表傾斜疊加形成的平面波傳播方向；（τR，τI）為高斯窗內(nèi)地震數(shù)據(jù)的傾斜疊加結(jié)果；（τR，τI）為（τR，τI）的Hilbert變換。（τR，τI）可由下式得到：

其中，Bj（ω，px）為第j個(gè)高斯窗內(nèi)通過(guò)傾斜疊加得到的頻率域結(jié)果，其具體表達(dá)式為

式中：C為與高斯窗函數(shù)相關(guān)的振幅系數(shù)；Fj（ω，kx）為地震記錄fj（x，t）的二維Fourier變換。

關(guān)于高斯窗函數(shù)的空間間隔與波束射線參數(shù)采樣間隔的關(guān)系，Hill［7，10］和Hale［8］均進(jìn)行了詳細(xì)的討論，我們采用了Hill［7，10］給出的參數(shù)選取準(zhǔn)則。

2 基于高斯射線束的角道集

角度域共成像點(diǎn)道集在AVA 分析和偏移速度分析方面有著非常廣闊的應(yīng)用前景。郭朝斌等［13］和李振春等［14］給出了利用粗網(wǎng)格走時(shí)信息計(jì)算角度域共成像點(diǎn)道集的方法。我們利用高斯束中心射線所攜帶的豐富的速度、方向信息，便捷直觀地給出了角度域共成像點(diǎn)道集的計(jì)算方法。

在中心射線坐標(biāo)系下，某一網(wǎng)格點(diǎn)處的走時(shí)可由對(duì)應(yīng)的中心射線通過(guò)二階泰勒展開(kāi)得到［15］（圖1）：

其中，τ0為中心射線走時(shí)；M可由動(dòng)力學(xué)射線追蹤得到。

利用坐標(biāo)變換，將中心射線坐標(biāo)系轉(zhuǎn)換為直角坐標(biāo)系，（13）式變?yōu)?/p>

其中，

式中：θ表示射線方向與Z軸的夾角。由（14）式對(duì)走時(shí)分別計(jì)算X方向和Z方向的偏導(dǎo)數(shù)，得到慢度px和pz：

其中，

得到成像點(diǎn)處的水平和垂直慢度后，我們可以分別計(jì)算出來(lái)自炮點(diǎn)的入射波場(chǎng)的射線角度φs 和來(lái)自射線束中心點(diǎn)的反射波場(chǎng)的射線角度φr：

在分別計(jì)算出φs 和φr 后，可通過(guò)計(jì)算二者之間的夾角來(lái)得到孔徑角。對(duì)于常規(guī)縱波勘探而言，反射角γ為孔徑角的一半，即

圖1 由中心射線計(jì)算成像點(diǎn)走時(shí)和射線方向原理示意圖解

3 模型測(cè)試

我們首先利用Marmousi國(guó)際標(biāo)準(zhǔn)模型對(duì)基于高斯束偏移的角度域共成像點(diǎn)道集算法進(jìn)行了測(cè)試。

圖2給出了Mamousi速度模型。圖3給出了采用本文方法對(duì)Mamousi模型數(shù)據(jù)進(jìn)行偏移得到的結(jié)果，可見(jiàn)成像效果良好。在輸出時(shí)我們選擇的角度范圍為［－60°，60°］，但由于觀測(cè)系統(tǒng)為右邊放炮左邊接收，所以在［0，60°］基本沒(méi)有能量，因此，我們只顯示［－60°，0］的角度域道集（圖4）。圖4a對(duì)應(yīng)圖3中最左邊的三角形位置，該處構(gòu)造相對(duì)簡(jiǎn)單，呈現(xiàn)層狀結(jié)構(gòu)，對(duì)應(yīng)的角道集也被拉平；圖4b對(duì)應(yīng)圖3中間的三角形位置，圖4c對(duì)應(yīng)圖3中最右邊的三角形位置，這兩處對(duì)應(yīng)的構(gòu)造復(fù)雜，橫向變速劇烈，角道集也同樣得到了拉平。由此可見(jiàn)，即使在如此復(fù)雜的構(gòu)造條件下，算子的成像精度依然很高，保證了最終高質(zhì)量的成像結(jié)果。

圖2 Marmousi模型

我們?cè)O(shè)計(jì)了一個(gè)包含直立斷層的復(fù)雜模型，該模型尺寸水平方向?yàn)?7 000m，垂直方向?yàn)? 000m，其速度場(chǎng)如圖5所示。利用有限差分技術(shù)模擬炮集記錄，觀測(cè)系統(tǒng)為中間放炮兩邊接收，檢波點(diǎn)間隔為50m，炮點(diǎn)間隔為100m，偏移距范圍為［－6 000m，6 000m］，共110炮，第1炮炮點(diǎn)位置在8 000m處。

圖5 包含直立斷層的復(fù)雜構(gòu)造速度模型

圖6a和圖6b分別給出了采用二階廣義屏偏移方法和高斯束偏移方法對(duì)該復(fù)雜模型的炮集記錄進(jìn)行偏移得到的結(jié)果。對(duì)比圖6a和圖6b可以看出，兩種方法整體上均能對(duì)構(gòu)造有比較好的成像效果，但是在陡傾角區(qū)域，高斯束偏移方法成像效果更好（圖6中紅色箭頭、紅色方框、紅色橢圓所示位置），二階廣義屏偏移方法無(wú)法對(duì)直立斷層進(jìn)行成像，而高斯束偏移方法能夠清晰、準(zhǔn)確地對(duì)直立斷層進(jìn)行成像。

圖7給出了圖6中紅色方框內(nèi)偏移結(jié)果的局部放大顯示。對(duì)比圖7a和圖7b可見(jiàn)，廣義屏算法在陡傾角構(gòu)造部位引入了較強(qiáng)的偏移噪聲，而且沒(méi)有將能量歸位到正確的位置上；而高斯束偏移結(jié)果信噪比高，同時(shí)將能量歸位到了正確的傾角構(gòu)造位置上。但由于觀測(cè)系統(tǒng)的原因，高斯束方法也未成像出中間的一段陡傾角構(gòu)造。圖8給出了圖6中紅色橢圓內(nèi)偏移結(jié)果的局部放大顯示。對(duì)比圖8a和圖8b可見(jiàn)，對(duì)于陡傾角構(gòu)造，高斯束的偏移結(jié)果要比廣義屏算法能量聚焦更好。

同時(shí)，我們對(duì)上述復(fù)雜構(gòu)造模型采用二階廣義屏偏移和高斯束偏移的計(jì)算效率進(jìn)行了對(duì)比，其結(jié)果如表1所示?？梢?jiàn)，高斯束偏移的計(jì)算效率要明顯高于二階廣義屏偏移的計(jì)算效率。

表1 高斯束偏移與二階廣義屏偏移有關(guān)參數(shù)

4 應(yīng)用效果分析

我們將高斯束疊前深度偏移算法應(yīng)用于西亞某地區(qū)的實(shí)際資料成像處理。研究區(qū)地質(zhì)條件復(fù)雜，速度橫向變化劇烈，存在高陡構(gòu)造，而且原始資料品質(zhì)差，經(jīng)過(guò)前期處理的偏移輸入道集如圖9a所示。由圖9a可見(jiàn)，該區(qū)地震資料信噪比較低，尤其是在中深層，有效信號(hào)被噪聲淹沒(méi)。圖9b為建立的深度域速度模型。圖10a和圖10b分別給出了對(duì)研究區(qū)地震資料進(jìn)行克?；舴蚱坪透咚故频慕Y(jié)果；圖11a和圖11b分別給出了CMP50和CMP250處高斯束偏移生成的角度域共成像點(diǎn)道集。對(duì)比圖10a和圖10b可見(jiàn)，高斯束偏移結(jié)果信噪比較高，同相軸連續(xù)性更好，特別是在陡傾角區(qū)域（圖中橢圓區(qū)域），高斯束偏移方法實(shí)現(xiàn)了對(duì)高陡構(gòu)造的優(yōu)質(zhì)成像；在深層由于輸入道集信噪比低，克希霍夫方法抗噪能力相對(duì)較弱，成像效果較差，而高斯束偏移剖面中深層信噪比較高，同相軸連續(xù)性較好。

圖11 角度域共成像點(diǎn)道集

5 結(jié)束語(yǔ)

高斯束偏移方法保留了射線方法靈活高效、無(wú)傾角限制的優(yōu)勢(shì)，同時(shí)克服了焦散、陰影區(qū)、多值走時(shí)等問(wèn)題；基于射線中心坐標(biāo)系的動(dòng)力學(xué)射線追蹤得到的慢度信息可直接用于計(jì)算共成像點(diǎn)角度道集。復(fù)雜模型和實(shí)際數(shù)據(jù)的測(cè)試結(jié)果驗(yàn)證了高斯束偏移方法在陡傾角構(gòu)造、低信噪比資料成像方面的優(yōu)勢(shì)。

本研究?jī)H限于各向同性介質(zhì)的構(gòu)造成像，在振幅保持以及各向異性高斯束偏移方面有待更深入的研究。

［1］Virgilio M，Simone R，Daniele C.Gaussian beam migration：prestack，common－shot，TTI，true amplitude in angular domain［J］.Expanded Abstracts of 78thAnnual Internat SEG Mtg，2008，2236－2239

［2］Cockshott I.Specular beam migration—a low cost 3D prestack depth migration［J］.Expanded Abstracts of 76thAnnual Internat SEG Mtg，2006，2634－2637

［3］Gao F，Zhang P，Dirks V.Fast beam migration—a step toward interactive imageing［J］.Expanded Abstracts of 75thAnnual Internat SEG Mtg，2005，2470－2474

［4］ˇCervenyV，Popov M M，Pˇsenˇcík I.Computation of wave fields in inhomogeneous media—gaussian beam approach［J］.Geophysical Journal Royal Astronomical Society，1982，70：109－128

［5］ˇCervenyV，Pˇsenˇcík I.Gaussian beams in elastic 2－D laterally varying layered structures［J］.Geophysical Journal Royal Astronomical Society，1984，78：65－91

［6］Costa C A，Raz S，Kosloff D.Gaussian beam migration［J］.Expanded Abstracts of 59thAnnual Internat SEG Mtg，1989，1169－1171

［7］Hill N R.Gaussian beam migration［J］.Geophysics，1990，55（11）：1416－1428

［8］Hale D.Migration by the Kirchhoff，slant stack，and Gaussian beam methods［R］.Colorado：CWP Annual Project Review Meeting，1992

［9］Hale D.Computational aspects of Gaussian beam migration［R］.Colorado：CWP Annual Project Review Meeting，1992

［10］Hill N R.Prestack Gaussian beam depth migration［J］.Geophysics，2001，66（4）：1240－1250

［11］Gray S H.Gaussian beam migration of common－shot records［J］.Geophysics，2005，70（4）：S71－S77

［12］Claerbout J F.Imaging the earth's interior［M］.Oxford：Blackwell Scientific Publication，1985：376－378

［13］郭朝斌，李振春，岳玉波.高斯束成像技術(shù)及其應(yīng)用［J］.石油物探，2011，50（1）：38－44，58 Guo C B，Li Z C，Yue Y B.Gaussian beam migration and its application［J］.Geophysical Prospecting for Petroleum，2011，50（1）：38－44，58

［14］李振春，岳玉波，郭朝斌，等.高斯波束共角度保幅深度偏移［J］.石油地球物理勘探，2010，45（3）：360－365 Li Z C，Yue Y B，Guo C B，et al.Gaussian beam common angle preserved－amplitude migration［J］.Oil Geophysical Prospecting，2010，45（3）：360－365

［15］ˇCerveny V.Seismic ray theory［M］.New York：Cambridge University Press，2001：234－400