盧寶坤，張劍鋒，蒲泊伶

1 中國科學(xué)院地質(zhì)與地球物理研究所，北京 100029

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

1 引 言

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

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

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

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

2 基于單程波的疊前時(shí)間偏移方法

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

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

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

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

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

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

對(duì)（3）式用穩(wěn)相點(diǎn)原理[17］求解，可求得漸進(jìn)解為

（4）式中

是炮點(diǎn) （xs，ys，0）下傳到成像點(diǎn) （x，y，T）處的走時(shí).

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

由（9）式求解pr進(jìn)而求得，代入（5）式可以求走時(shí)，代入可以求取幾何擴(kuò)散振幅補(bǔ)償項(xiàng)，代入可以求取透射衰減振幅補(bǔ)償項(xiàng).

3 速度梯度對(duì)走時(shí)、角度、幅值影響的定量分析

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

3.1 速度梯度對(duì)走時(shí)影響的定量分析

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

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

3.2 速度梯度對(duì)角度影響的定量分析

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

由（11）、（12）式可以看出，可以直接由炮點(diǎn)、接收點(diǎn)以及像點(diǎn)的水平面坐標(biāo)獲得射線參數(shù)其計(jì)算效率優(yōu)于迭代求解方法[14］.

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

其中，px，，py分別是x，y方向的射線參數(shù)，x0，y0是像點(diǎn)坐標(biāo)，v是成像點(diǎn)的層速度，θz是入射波與z軸的夾角.

可以看出求取角度時(shí)采用的是成像點(diǎn)處的速度v，而不是像文獻(xiàn)[11-13］那樣采用均方根速度.

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

圖3 入射波與z軸夾角對(duì)比實(shí)線為解析解，短劃線為本文方法解，灰線為文獻(xiàn)[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 速度梯度對(duì)振幅影響的定量分析

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

4 角度域疊前時(shí)間偏移

4.1 入射角求取

對(duì)（13）式應(yīng)用余弦定理進(jìn)一步得到：

圖4 振幅補(bǔ)償項(xiàng)對(duì)比實(shí)線代表幾何擴(kuò)散補(bǔ)償項(xiàng)加射線參數(shù)精確值計(jì)算的透射損失補(bǔ)償項(xiàng)，短劃線代表幾何擴(kuò)散補(bǔ)償項(xiàng)加射線參數(shù)近似值計(jì)算的透射損失補(bǔ)償項(xiàng)，灰線代表只有幾何擴(kuò)散補(bǔ)償項(xiàng).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.利用炮點(diǎn)和成像點(diǎn)的幾何關(guān)系，可進(jìn)一步得到

這樣可解得：

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

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

其中h是半偏移距，進(jìn)一步求得入射角為

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

4.2 角度域疊前時(shí)間偏移角道集獲得

獲得較準(zhǔn)確的走時(shí)、振幅值和入射角后，就可以利用波動(dòng)方程深度偏移的反褶積成像條件[21］進(jìn)行角度域疊前時(shí)間偏移計(jì)算.若設(shè)震源是一時(shí)間脈沖，則對(duì)一個(gè)地震道有成像結(jié)果：

式中F′（t）是f（ω）對(duì)應(yīng)的時(shí)域函數(shù)的一階導(dǎo)數(shù)，As和ts是炮點(diǎn) （xs，ys，0）至成像點(diǎn) （x，y，T）的幅值和走時(shí)，Ar和tr是求得的檢波點(diǎn)至成像點(diǎn)的幅值和走時(shí).

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

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

4.3 角度域疊前時(shí)間偏移流程

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

（1）對(duì)疊前地震數(shù)據(jù)求導(dǎo)數(shù)，計(jì)算走時(shí)以及幅值表和偏移孔徑表；

（2）在橫向均勻的速度場下對(duì)每道數(shù)據(jù)計(jì)算脈沖響應(yīng)，對(duì)每個(gè)CDP點(diǎn)將各道數(shù)據(jù)的脈沖響應(yīng)按數(shù)據(jù)的偏移據(jù)分類、部分疊加，形成CRP偏移距道集；

（3）對(duì)各CRP道集按橫向均勻的疊加速度做反動(dòng)校正、然后作動(dòng)校正并將這一速度作為該點(diǎn)新的疊加速度；

（4）在新的疊加速度下計(jì)算脈沖響應(yīng)，對(duì)同一個(gè)CRP點(diǎn)計(jì)算入射角，把同一入射角數(shù)據(jù)疊加，形成共成像點(diǎn)角道集.

（5）對(duì)同一CRP角道集數(shù)據(jù)作剩余動(dòng)校正，拉平同相軸.

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

5 數(shù)值算例

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

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

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

表1 反射界面參數(shù)Table 1 Parameters of three reflections

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

6 實(shí)際海上數(shù)據(jù)算例

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

7 結(jié) 語

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

圖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，F(xiàn)oster D.Common angle imaging conditions for pre-stack depth migration.70thAnnual International Meeting，SEG，Expanded Abstracts.2000：830-833.

[3]Sava P C，F(xiàn)omel S.Angle-domain common-image gathers by wavefield continuation methods.Geophysics，2003，68（3）：1065-1074.

[4]陳生昌，曹景忠，馬在田.疊前地震數(shù)據(jù)的平面波深度偏移法.地球物理學(xué)報(bào)，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小波束疊前深度偏移的角度域共成像道集.地球物理學(xué)報(bào)，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]王西文，劉文卿，王宇超等.共反射角疊前偏移成像研究及應(yīng)用.地球物理學(xué)報(bào)，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.波動(dòng)方程偏移成像陰影的照明補(bǔ)償.地球物理學(xué)報(bào)，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]張輝，王成祥，張劍鋒.基于單程波方程的角度域照明分析.地球物理學(xué)報(bào)，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]張江杰，張劍鋒.波動(dòng)方程偏移速度建模：一種直接反演方法.地球物理學(xué)報(bào)，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疊前深度偏移速度分析方法.地球物理學(xué)報(bào)，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]程玖兵，王楠，馬在田.表驅(qū)三維角度域Kirchhoff疊前時(shí)間偏移成像方法.地球物理學(xué)報(bào)，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疊前時(shí)間偏移角度道集.地球物理學(xué)報(bào)，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.