角度域疊前時間偏移：非均勻介質中的單程波方法

盧寶坤，張劍鋒，蒲泊伶

1 中國科學院地質與地球物理研究所，北京 100029

2 中國石油勘探開發研究院，北京 100083

1 引 言

振幅隨入射角變化（AVA）技術廣泛應用于油氣儲層反演中，預測巖石物性，提高了儲層預測的準確性.工業生產中通常用振幅隨偏移距變化（AVO）技術來替代AVA技術或者把偏移距通過簡單的幾何計算影射成角度進行AVA分析.偏移距與入射角度只有在地下是水平層狀介質才具有簡單的幾何影射關系，因此當地下介質比較復雜時，直接采用上述方法計算出的入射角存在較大的誤差.

近年來，用于速度分析以及AVA分析的角度域道集（ADCIG）引起了學者廣泛的關注，ADCIG相對于偏移距共成像點道集（ODCIG）可以很大程度上減少假象，同時可以提供地下的入射角信息.Berkhout[1］提出利用地面平面波道集的偏移成像得到平面波參數域共成像點道集的方法；Moster和Foster[2］提出利用雙平方根波場得到地下成像點的偏移距kh域波場，然后由投影定理把kh轉化為偏移距射線參數ph，由此獲得成像點處的ph域共成像點道集.Sava和Formel[3］利用單程波偏移成像的雙平方根方程，并結合傾斜疊加，生成了波動方程角度域成像道集.國內學者也在角度域偏移成像方法上進行了研究[4-6］，并且把角度域偏移方法應用于照明分析[7-8］以及速度分析方面[9-10］.這些研究多是基于波動方程深度偏移的，要么數據體龐大且內存需要大，要么計算量大而且存在不穩定性問題.

與疊前深度偏移方法相比，疊前時間偏移方法（PSTM）可以對一類傾角、斷層較為復雜但速度橫向變化不是很劇烈的構造較好成像，具有較高的計算效率.只是需要疊加速度，即可簡單地通過速度掃描等方式得到恰當的速度模型，回避了疊前深度偏移方法難以解決的困難，因此采用疊前時間偏移技術輸出共成像點的入射角道集，用于AVA分析是較好的選擇.Fomel和Prucha[11］討論了角度域疊前時間偏移的基本方法，其入射角計算是通過偏移距、成像深度（時間）以及偏移速度等參數經過簡單的幾何關系來表達的.Zheng[12］采用基于時間域直射線近似Kirchhoff疊前時間偏移角度域成像方法，在高分辨率裂隙中得到應用.其走時計算采用雙平方根算子，入射角計算直接采用偏移距、炮點和接收點分別到成像點的距離通過余弦定理得到；Perez等[13］在角度域疊前偏移中加入波形校正對偏移中的波形拉伸給予校正.其角度計算是分別計算炮點和接收點到像點的走時向量，然后求取二者的夾角，取其一半作為入射角度.上述方法計算角度時使用的是偏移速度，而不是成像點的層速度，計算的角度比實際角度小；程玖兵等[14］提出走時基于射線追蹤的表驅動Kirchhoff疊前時間偏移角度域算法，走時考慮射線彎曲效應，取得了好的成像效果；由于射線參數是逐層遞推獲得的，因此射線參數和入射角度計算比較準確.鄒振等[15］在假設水平射線參數不變的基礎上，采用炮點到像點走時對炮點到像點的水平距離求導或者采用接收點到像點走時對像點到接收點的距離求導來獲得射線參數.走時計算考慮射線彎曲效應，相對直射線的方法[11-13］入射角角度范圍增大.

本文從“真幅值”單程波相移法[16］出發，直接求取成像點處的射線參數，進而獲得入射角；計算入射角采用成像點處的層速度，保證了入射角的精度.由于入射角是以度為計量單位，因此需要的層速度不需要達到深度偏移成像的速度精度，可以通過如Dix公式等簡單方法獲得，速度建模簡單易行；由于該方法是從“真幅值”單程波方程出發，較好的補償了地震波傳播的幾何擴散效應和透射衰減，因而很好的保留了地震波的振幅特征，其計算結果用于疊前反演可以得到較準確的效果；求解入射角時充分利用走時計算產生的數據，同時對反余弦計算提前建表儲存，計算時間比單獨疊前時間偏移成像所需計算時間增加有限.因此本文方法既可以較精確計算入射角，又可以兼顧計算效率，可以方便的應用于工業處理流程.

2 基于單程波的疊前時間偏移方法

由于生成角道集的目的是服務于疊前反演，因此幅值問題尤為重要.常規的單程波算子僅能考慮地震波的幾何擴散，未能考慮介質速度變換導致的透射衰減，因而不能很好的解決幅值問題.本文從“真幅值”的單程波方程出發研究疊前時間偏移方法，推導出走時、振幅以及入射波與反射波的夾角計算公式，進而獲得具有更準確幅值的角道集.

針對介質的非均勻性，沿用波動方程疊前深度偏移的處理方法，首先假定介質是橫向不變的，求得地震波傳播的走時、幅值后，再允許速度橫向變換，用這一解代表非均勻介質的解；擬微分算子理論的研究結果表明，這一解確實反應了地震波傳播的主要特征，并且現行的波動方程深度偏移的實踐應用證明了該處理方法的有效性.

在三維層狀介質情況下，波數-頻率域的真幅值的單程波方程可以表示為

式中，kx和ky是水平波數，ω是頻率，kz=是橫向均勻的介質速度，是需要求解的波數-頻率域的波場，與常規單程波方程相比，增加了一個因子

式中，T是用走時表達的時間深度，f（ω）是炮點或檢波點的時域信號的傅里葉變換，px和py分別是水平射線參數，有px=kx／ω，py=ky／ω，v0表示地層第一層的層速度.

對（2）式做空間傅里葉反變換，得到空間-頻率域波場為

對（3）式用穩相點原理[17］求解，可求得漸進解為

（4）式中

是炮點 （xs，ys，0）下傳到成像點 （x，y，T）處的走時.

引入新的變量pr，定義則（7）和（8）式可轉化為單變量的方程：

由（9）式求解pr進而求得，代入（5）式可以求走時，代入可以求取幾何擴散振幅補償項，代入可以求取透射衰減振幅補償項.

3 速度梯度對走時、角度、幅值影響的定量分析

在走時計算時，通常采用均方根速度代替層速度，使用級數展開到偏移距的二次項，這種近似不能很好的描述含有縱向速度梯度的地層，尤其對于大偏移距數據.下面用常速度梯度模型來定量分析速度梯度對地震波的走時、入射角以及振幅的影響.設計了如下速度模型A：速度沿縱向變化，速度梯度保持常量為1200m／s2，選擇時間深度為2s（此時深度為6400m）的成像點作為分析點.在上述速度模型下考慮到入射波與反射波關于成像點對稱，我們可以選擇炮點的下行波來進行分析.

3.1 速度梯度對走時影響的定量分析

在模型A下，從炮點正傳到成像點再反傳到接收點的走時是炮點正傳到成像點的2倍.圖1a中實線是按照公式（5）計算的走時（精確解）隨偏移距的變化曲線，虛線是按照走時近似公式[18-20］

計算的走時（記為Tvrms），其中為成像點的坐標向量，s、g分別表示震源與接收點的坐標向量，τξ為成像點對應的雙程走時深度，Vrms為均方根速度.圖1b是二者的差.可以看出隨著偏移距的增大，二者的差距變大，說明速度梯度的作用變大，因此對于大偏移距數據進行走時計算時，應該考慮速度梯度的影響.同時我們選擇偏移距為10000m的數據做動校正（NMO）計算，此時偏移距是成像點深度大約1.5倍.圖2中實線是未做動校正前的波形，虛線是動校正后的波形，可以看出經過動校正后，波形的周期明顯變大，產生了拉伸，在實際資料處理中往往做切除處理.這時再對圖1中偏移距小于10000m的數據進一步分析，可以看出采用近似公式（10）計算的走時誤差小于15ms，此時我們可以采用走時近似公式（10）替代解析公式（5）進行走時計算.

3.2 速度梯度對角度影響的定量分析

采用均方根速度Vrms代替介質速度，（5）式近似變化為

由（11）、（12）式可以看出，可以直接由炮點、接收點以及像點的水平面坐標獲得射線參數其計算效率優于迭代求解方法[14］.

獲得了成像點處的射線參數（（11）、（12）式）后，就可以進一步求取炮點正傳波場與軸的夾角.由Snell定理通過（11）、（12）式可以得到如下的關系：

其中，px，，py分別是x，y方向的射線參數，x0，y0是像點坐標，v是成像點的層速度，θz是入射波與z軸的夾角.

可以看出求取角度時采用的是成像點處的速度v，而不是像文獻[11-13］那樣采用均方根速度.

采用解析公式θz=sin-1（v＊p）獲得入射波與z的夾角（記為Exact），結果如圖3實線所示，p是射線參數精確值；同時分別采用（14）式以及文獻[5-7］的方法計算夾角近似值.其中圖3中短劃線是式（14）的計算值（記為Angle-Vint），點線是采用文獻[11-13］方法的計算值（記為 Angle-Vrms）.可以看出采用式（14）計算的角度，其與精確解（實線）相比誤差很小；而文獻[11-13］不考慮成像點處的速度，其角度計算結果（Angle-Vrms）誤差較大.因此進行角度計算時要考慮介質速度的影響.由于射線參數的精確解在疊前時間偏移中往往不容易得到，所以實際計算中可以選擇（14）式計算角度，本文選擇該公式用于角度計算.

圖3 入射波與z軸夾角對比實線為解析解，短劃線為本文方法解，灰線為文獻[5-7］的解.Fig.3 Comparison of angle between incident wave and zaxis The solid line denotes exact，dash line denotes that presented by this paper（Angle-Vint）and gray line denotes that presented by reference[5-7］（Angle-Vrms）.

3.3 速度梯度對振幅影響的定量分析

綜上所述，本文在角度域疊前偏移成像中采用如下的計算方法：在走時計算時采用走時近似公式（10）計算；夾角計算時采用成像點處的層速度；采用均方根速度獲得射線參數近似值進而得到振幅幾何擴散補償項和透射衰減補償項.

4 角度域疊前時間偏移

4.1 入射角求取

對（13）式應用余弦定理進一步得到：

圖4 振幅補償項對比實線代表幾何擴散補償項加射線參數精確值計算的透射損失補償項，短劃線代表幾何擴散補償項加射線參數近似值計算的透射損失補償項，灰線代表只有幾何擴散補償項.Fig.4 Comparison of amplitude compensation factors The solid line denotes amplitude compensation for geometric divergence plus for loss of refraction obtained by exact rayparameters（Spread+P-exact），dash line denotes amplitude compensation for geometric divergence plus that for loss of refraction obtained by approximate ray-parameters（Spread+P-vrms）and gray line denotes compensation for geometric divergence only（Spread only）.

式中，定義ρb=v／vrms.利用炮點和成像點的幾何關系，可進一步得到

這樣可解得：

同理可解得反傳波場的三個方向余弦，利用它們之間的點積，可求得兩者間的夾角.令

可求得入射波與界面法線的夾角的余弦值為

其中h是半偏移距，進一步求得入射角為

從（17）式可以看出求解入射角采用的是層速度，而不是均方根速度，據此計算的入射角較準確.由于是以度為計量單位，因此入射角計算需要的層速度模型精度比深度偏移成像低，可以通過如Dix公式等簡單方法獲得，大大節省了速度建模的時間.

4.2 角度域疊前時間偏移角道集獲得

獲得較準確的走時、振幅值和入射角后，就可以利用波動方程深度偏移的反褶積成像條件[21］進行角度域疊前時間偏移計算.若設震源是一時間脈沖，則對一個地震道有成像結果：

式中F′（t）是f（ω）對應的時域函數的一階導數，As和ts是炮點 （xs，ys，0）至成像點 （x，y，T）的幅值和走時，Ar和tr是求得的檢波點至成像點的幅值和走時.

公式（18）成像條件表明：對任一道地震道和任一成像點，對地震信號求一階導數，依據成像點處的偏移速度，計算總走時tr+ts和權系數Ar／As，在地震信號的一階導數上拾取ts+tr時刻的值并乘以權系數，即得到該地震道在該成像點的偏移幅值.同時由公式（17）可以求得成像點處的入射角γ，這樣就可將該成像點的偏移幅值記錄在入射角γ的位置上.將所有地震道的脈沖響應結果按入射角排列、部分疊加，就可以獲得角道集；如果把同一個像點不同角度的偏移數據疊加，也可以獲得最終的偏移結果，由此可以看出：角道集計算可以在疊前時間偏移成像過程中獲得，不需要另外計算過程.

由公式（18）可以看出，求取入射角時，可以充分利用走時計算中ps、pr、T＊Vrms等已知量，同時對cos-1提前建表，ρb提前計算并儲存，這樣就可以大大節省入射角計算所需要的計算時間，使其更好的應用于工業生產.

4.3 角度域疊前時間偏移流程

角度域疊前時間偏移可以同時獲得角度域道集與偏移距道集，其偏移流程如下：

（1）對疊前地震數據求導數，計算走時以及幅值表和偏移孔徑表；

（2）在橫向均勻的速度場下對每道數據計算脈沖響應，對每個CDP點將各道數據的脈沖響應按數據的偏移據分類、部分疊加，形成CRP偏移距道集；

（3）對各CRP道集按橫向均勻的疊加速度做反動校正、然后作動校正并將這一速度作為該點新的疊加速度；

（4）在新的疊加速度下計算脈沖響應，對同一個CRP點計算入射角，把同一入射角數據疊加，形成共成像點角道集.

（5）對同一CRP角道集數據作剩余動校正，拉平同相軸.

（6）如果需要偏移成像結果，可以通過對角道集疊加獲得.

5 數值算例

為了驗證本文提出的角度域疊前時間偏移的效果，模擬了海上勘探數據模型，其波阻抗模型如圖5所示.該模型水平方向距離3300m，深度為1650m，包含三個低速帶低波阻抗砂體，模型的主要構造為大角度的正斷層，最高角度可以達到60°.由于時間偏移的適用條件所限，盡管設計的構造非常復雜，但是模型的速度橫向變化平緩.從圖6偏移成像剖面上可以看出，其界面成像較好，斷點及砂體端點收斂比較清晰；其能量對比反映了界面上下的波阻抗關系.說明該方法可以很好的對模型數據成像，且很好的保留了反射波的振幅信息.

為了獲得共成像點角道集的振幅隨入射角的變化關系，特選取圖5箭頭所示位置的角道集進行分析，如圖7所示.可以看出角道集上反射波同相軸已經拉平，從上到下有效角度范圍變小.其中標示1的地層角度可以達到40°以上.讀取圖7箭頭所示三個位置的振幅，同時采用聲波反射系數隨入射角的變化公式（19）[22］進行反射系數理論值計算：

其中，Rpp（θ1）代表反射界面的反射系數，VP1，VP2分別是界面上下兩層的P波速度，ρ1，ρ2分別是界面上下兩層的密度.VP是波阻抗.θ1，θ2分別是界面的入射角和反射角，二者滿足Snell定律.其與圖7三個反射界面對應的地層參數如表1所述.計算結果如圖8所示，其中實線代表從圖7中標志序號的同相軸讀出的振幅值，虛線是按照公式（19）計算的理論值，二者分別進行歸一化處理.可以看出兩者很好的吻合，說明該方法計算的角道集有很好的保幅效果，基于這一角道集進行疊前反演是可以得到好的反演結果.

表1 反射界面參數Table 1 Parameters of three reflections

圖7 圖5箭頭所示位置處的角道集Fig.7 The angle-gather of the station in Fig.5marked by an arrow

6 實際海上數據算例

針對我國南海某實際三維地震資料，采用本文提出的疊前時間偏移角道集計算方法進行了實際計算.該數據體共3700炮，總數據量21G，其中道距13.33m，采樣2ms；成像CDP間隔6.665m，成像深度5s.選取一條Inline方向的成像剖面用于分析，如圖9所示.可以看出資料的信噪比較高，構造形態清晰.分別抽取了入射角度為8°和16°剖面，如圖10和圖11所示.可以看到二者的構造形態與圖9一致，特別是三者圖中橢圓所圈位置的斷層形態一致，說明入射角角度計算正確.為了進一步分析單個角道集的特征，特抽取6km處的角道集，如圖12所示.可以看出反射波同相軸已經拉平，且連續性較好.入射角有效范圍大，角度數據豐富.其淺層（1.5s）以上部分有效角度可以達到40°，可以很好的用于后續的AVA反演.

7 結 語

本文從“真幅值”單程波方程出發，構建了可直接生成角度域成像道集的疊前時間偏移方法與角度域疊前偏移流程.在角度以及振幅計算上充分考慮了速度梯度的影響，因此可以獲得較準確的入射角度；較好的補償了地震波傳播的幾何擴散效應和透射損失，因此可以得到較準確的振幅.可以在疊前時間偏移成像過程完成角道集的計算，入射角計算充分利用走時計算的過程結果，計算效率高.理論模型、實際資料偏移疊加結果以及抽取的角道集驗證了本文方法的有效性.

圖12 6km處的角道集Fig.12 The angle-gather of 6km

（References）

[1]Berkhout A J.Pushing the limits of seismic imaging，Part II：Integration of pre-stack migration，velocity estimation，and AVO analysis.Geophysics，1997，62（3）：954-969.

[2]Mosher C，Foster D.Common angle imaging conditions for pre-stack depth migration.70thAnnual International Meeting，SEG，Expanded Abstracts.2000：830-833.

[3]Sava P C，Fomel S.Angle-domain common-image gathers by wavefield continuation methods.Geophysics，2003，68（3）：1065-1074.

[4]陳生昌，曹景忠，馬在田.疊前地震數據的平面波深度偏移法.地球物理學報，2003，46（6）：821-826.Chen S C，Cao J Z，Ma Z T.A method of plane wave depth migration for pre-stack seismic data.Chinese J.Geophys.（in Chinese），2003，46（6）：821-826.

[5]陳凌，吳如山，王偉君.基于Gabor-Daubechies小波束疊前深度偏移的角度域共成像道集.地球物理學報，2004，47（5）：876-884.Chen L，Wu R S，Wang W J.Common angle image gathers obtained from Gabor-Daubechies beamlet pre-stack depth migration.Chinese J.Geophys.（in Chinese），2004，47（5）：876-884.

[6]王西文，劉文卿，王宇超等.共反射角疊前偏移成像研究及應用.地球物理學報，2010，53（11）：2723-2738.Wang X W，Liu W Q，Wang Y C，et al.Research and application of common reflection angle domain pre-stack migration.Chinese J.Geophys.（in Chinese），2010，53（11）：2732-2738.

[7]陳生昌，馬在田，Wu R S.波動方程偏移成像陰影的照明補償.地球物理學報，2007，50（3）：844-850.Chen S C，Ma Z T，Wu R S.Illumination compensation for wave equation migration shadow.Chinese J.Geophys.（in Chinese），2007，50（3）：844-850.

[8]張輝，王成祥，張劍鋒.基于單程波方程的角度域照明分析.地球物理學報，2009，52（6）：1606-1614.Zhang H，Wang C X，Zhang J F.One-way wave equation based illumination analysis in angle domain.Chinese J.Geophys.（in Chinese），2009，52（6）：1606-1614.

[9]張江杰，張劍鋒.波動方程偏移速度建模：一種直接反演方法.地球物理學報，2011，54（3）：835-844.Zhang J J，Zhang J F.Wave equation based migration velocity estimation：a direct inversion approach.Chinese J.Geophys.（in Chinese），2011，54（3）：835-844.

[10]杜啟振，李芳，侯波等.角度域彈性波Kirchhoff疊前深度偏移速度分析方法.地球物理學報，2011，54（5）：1327-1339.Du Q Z，Li F，Hou B，et al.Angle-domain migration velocity analysis based on elastic-wave Kirchhoff pre-stack depth migration.Chinese J.Geophys.（in Chinese），2011，54（5）：1327-1339.

[11]Fomel S，Prucha M.Angle-gather time migration：Stanford.Exploration Project Report SEP 100，1999，141-151.

[12]Zheng Y.Seismic Azimuthal Anisotropy and Fracture Analysis from PP Reflection Data.Calgary：University of Calgary，2006.

[13]Perez G， Marfurk K J.Improving lateral and vertical resolution of seismic images by correcting for wavelet stretch in common-angle migration.Geophysics，2007，72（6）：C95-C104.

[14]程玖兵，王楠，馬在田.表驅三維角度域Kirchhoff疊前時間偏移成像方法.地球物理學報，2009，52（3）：792-800.Cheng J B，Wang N，Ma Z T.Table-driven 3-D angledomain imaging approach for Kirchhoff pre-stack time migration.Chinese J.Geophys.（in Chinese），2009，52（3）：792-800.

[15]鄒振，劉洪，劉紅偉.Kirchhoff疊前時間偏移角度道集.地球物理學報，2010，53（5）：1207-1214.Zou Z，Liu H，Liu H W.Common angle gathers based on Kirchhoff pre-stack time migration.Chinese J.Geophys.（in Chinese），2010，53（5）：1207-1214.

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

[17]Bleistein N，Cohen J K，Stochwell J W Jr.Mathematics of Multidimensional Seismic Imaging，Migration and Inversion.New York：Springer，2001.

[18]Yilmaz O.Pre-stack Partial Migration.Stanford：Stanford University，1979.

[19]Bancroft J C，Geiger H D，Margrave G F.The equivalent offset method of pre-stack time migration.Geophysics，1998，63（6）：2042-2053.

[20]Biondi B.3DSeismic Imaging.Tulsa：Society of Exploration Geophysicists，2006.

[21]Zhang J F，Wapenaar K.High-resolution depth imaging with sparseness-constrained inversion.Geophysical Prospecting，2006，54（1）：49-62.

[22]Smith G C，Gidlow P M.Weighted stacking for rock property estimation and detection of gas.Geophysical Prospecting，1987，35（9）：993-1014.