TI介質(zhì)局部角度域射線追蹤與疊前深度偏移成像

段鵬飛，程玖兵＊，陳三平，何光明

1 同濟大學海洋與地球科學學院，上海 200092

2 中國石油川慶鉆探工程有限公司地球物理勘探公司，成都 610213

1 引 言

為了改善地震波照明與成像、噪聲壓制以及油氣儲層的預測和描述，長偏移距、寬方位的地震數(shù)據(jù)采集越來越多，以往針對有限偏移距、窄方位地震數(shù)據(jù)的各向同性介質(zhì)假設與在此基礎上的地震成像方法明顯受到了挑戰(zhàn).忽略各向異性帶來的誤差會引起反射波歸位不準，繞射波收斂不徹底，能量不聚焦，對長偏移距、寬方位數(shù)據(jù)尤為突出.這就要求地震偏移成像與速度模型建立方法考慮速度各向異性，否則花費巨大代價觀測到的地震資料得不到合理的利用.大量的觀測與研究發(fā)現(xiàn)，許多頁巖與薄互層沉積地層等在地震波長尺度下均可等效成TI介質(zhì).而這些地層在全球范圍內(nèi)分布非常廣泛.因此，在考慮各向異性地震波成像問題時，VTI介質(zhì)和TTI介質(zhì)是最常用的等效模型.

疊前深度偏移是強橫向非均勻介質(zhì)復雜構造成像與速度模型建立依賴的關鍵技術.其算法實現(xiàn)要么基于射線理論，如Kirchhoff偏移和高斯束偏移，要么基于波動理論，如單程波方程深度延拓偏移和雙程波方程逆時延拓偏移.近10多年來，各向異性介質(zhì)深度偏移方法也得到了極大的發(fā)展，先后出現(xiàn)了 TI介質(zhì) Kirchhoff偏移[1］、高斯束偏移[2］、單程波方程偏移[3］與逆時偏移[4］等深度域成像方法.盡管波動方程偏移存在精度上的優(yōu)勢，射線理論基礎上的偏移方法因其在靈活性、面向局部目標的成像能力以及計算成本等優(yōu)勢，在復雜構造成像尤其是速度模型建立過程中得到廣泛應用.目前主要地震數(shù)據(jù)處理軟件中的深度域偏移速度模型構建都仍以Kirchhoff偏移作為引擎.

在復雜介質(zhì)條件下，即使偏移速度是合理的，傳統(tǒng)的偏移距域和炮域共成像點道集都可能存在假象干擾.為此，近十幾年來人們一直在致力于研究射線理論或波動理論基礎上的角度域成像方法.基于射線理論和廣義拉冬變換（GRT），de Hoop[5］提出了共散射角偏移／反演理論.隨后Xu等[6］和Brandsberg-Dahl等[7］提出了Kirchhoff疊前深度偏移共散射角成像方法，并用于地震成像和速度模型建立過程中.Bleistein[8］也系統(tǒng)地闡述了共散射角成像／反演理論.最近，Koren等[9］在局部角度域成像理論框架下提出了方向型和反射型共成像點道集的產(chǎn)生方法及其用途.Cheng等[10］提出了適應TI與方位各向異性介質(zhì)的方位保真局部角度域疊前時間偏移成像方法并展示了在儲層成像與描述中的初步應用.

局部角度域Kirchhoff疊前深度偏移算法的核心在于穩(wěn)健、快速地計算地震射線的走時與方向信息.在傳統(tǒng)Kirchhoff疊前深度偏移過程中，程函方程有限差分解法[11］與波前重建算法[12］被廣泛用于走時表的計算[1，13］.然而，對局部角度域成像與反射走時層析基礎上的偏移速度分析而言，射線追蹤算法顯得更有吸引力，因為它除了計算走時，還可以顯式地得到射線路徑及其方向信息.不過，傳統(tǒng)的各向異性介質(zhì)射線追蹤方程是以剛度系數(shù)而不是Thomsen參數(shù)表示的[14］，不方便數(shù)值計算，效率也較低[14-16］.為此，一些學者重新推導了Thomsen參數(shù)表征的射線方程.例如，Alkhalifah[17］基于聲學近似推導出了VTI介質(zhì)的射線方程.Zhu[15-16］等推導了一種基于相速度、適用于一般各向異性介質(zhì)的射線方程.

本文根據(jù)強橫向非均勻各向異性介質(zhì)地震波成像與偏移速度模型建立的需要，研究基于射線理論的TI介質(zhì)局部角度域疊前深度偏移成像方法.首先討論與對比兩種各向異性射線追蹤算法，然后論述局部角度域疊前深度偏移成像原理及其算法實現(xiàn)，最后借助國際上通用的理論模型檢驗算法的可靠性.

2 各向異性介質(zhì)射線追蹤方法

各向異性介質(zhì)中的運動學射線方程最早由Cerveny[14］給出.下文簡單回顧其推導過程.一般各向異性介質(zhì)無源的彈性波動方程可以表示為

其中，ui為位移，cijkl為介質(zhì)的剛度系數(shù)，i，j，k，l=1，2，3，ρ為密度，t為走時.其頻率域表達式為

其中，ω為圓頻率.在零階射線理論中，（2）式的近似解可寫成uk（xi，ω）=Uk（xi）eiωτ（xi），其中Uk（xi）和τ（xi）分別是射線上的振幅和走時.把這個解帶入（2）式，當ω→ ∞ 時，可得到Christoffel方程：

其中，Γjk=aijklpipl為Christoffel矩陣，aijkl=cijkl／ρ為密度歸一化剛度系數(shù)，pi=?τ／?xi為慢度矢量的各個分量.方程（3）對應一個標準的特征值問題，且特征值滿足G（pi，xi）=1.（3）式可改寫為

其中，gk是單位特征向量（即極化矢量）.（4）式兩邊同乘以gj，結合gkgk=1，可得到

對于程函τ（xi），方程（4）是一個非線性一階偏微分方程.這個程函方程可通過漢密爾頓方程求解，進而表示成一般各向異性介質(zhì)的運動學射線追蹤方程組[14］：

方程組等式最右側的函數(shù)非常復雜，計算不但費時，且需要在射線追蹤的每一步求解特征值問題.此外，方程組（6a）與（6b）用剛度系數(shù)來描述介質(zhì)的彈性性質(zhì)，這與實際地震資料處理中通常用Thomsen參數(shù)的情況不一致.為此，文中將討論兩種不需用剛度系數(shù)表征的各向異性射線方程.

2.1 相速度表示的各向異性介質(zhì)射線方程

為了克服剛度系數(shù)表示的射線方程的復雜性及其計算上的麻煩，Zhu[15-16］等重新推導了各向異性介質(zhì)中的運動學射線追蹤方程.根據(jù)文獻[14］，沿xi方向的群速度可表示為Vi=aijklplgjgk.于是，方程（6a）改寫為

其中，VGi為群速度對空間坐標xi的導數(shù).考慮到（4）式中特征值G及其偏導數(shù)?G／?xi都是pi的齊次方程，容易得到v2=G（xi，ni），故而有：

與各向同性介質(zhì)程函方程推導思路一樣，將平面波解帶入方程（10）可推導出VTI介質(zhì)的程函方程：

其中，ni為單位慢度矢量，v=v（xi，ni）為相速度.將（8）式帶入方程（6b）并聯(lián)立（7）式得到：

由于群速度可通過相速度計算得出，因此（9a）與（9b）式就組成了相速度表示、適應一般各向異性介質(zhì)的射線方程組.這樣就回避了傳統(tǒng)各向異性射線追蹤過程中每一步都要計算的特征值問題.注意，由于空間矢量x與單位慢度矢量n都是相速度方程v=v（xi，ni）的獨立變量，因此（9b）式右邊對相速度求偏導數(shù)時，其隱函數(shù)的鏈式求導中不依賴ni，只需對相速度表達式中出現(xiàn)的與空間坐標xi有關的參數(shù)對xi求導即可.

在許多地質(zhì)條件下，受構造運動或其它因素影響，一些橫向各向同性地層大多數(shù)情況下都非水平層狀，其對稱軸通常與垂向存在一定的夾角.這時，采用TTI模型來描述速度各向異性就更合理.附錄A給出了VTI與TTI介質(zhì)中相速度表示的射線追蹤方程的具體形式.

2.2 TI介質(zhì)聲學近似意義下的射線方程

聲學近似就是假設沿對稱軸方向qSV波的傳播速度為0，即VS0=0，這樣就可將原始的VTI介質(zhì)彈性波動方程及其頻散關系簡化.假設地下介質(zhì)為聲學介質(zhì)，由VTI介質(zhì)彈性波動方程及其頻散關系可推導出近似的qP波標量波動方程，進而得到相應的程函方程和射線方程[17］.研究表明，聲學近似對qP波運動學特征的負面影響基本可以忽略[17］.

根據(jù)VTI介質(zhì)qP波的頻散關系，聲學近似qP波波動方程滿足[17］：

其中VP0為qP波垂直速度，VNMO為NMO速度，η為反橢圓系數(shù)，且與Thomsen參數(shù)ε與δ存在如下關系：設η=0和VNMO=VP0，（12）式就退化成各向同性介質(zhì)的程函方程.通過特征值方法可進一步推導出描述射線路徑的常微分方程組.為此，將（12）式改寫為如下形式：

其中，τ代表沿著射線的走時，i對應x，y和z分量.該方程組描述了VTI介質(zhì)聲學近似意義下的射線路徑、走時及傳播方向信息.可以看出，式（14）與式（9）形式上非常相似，且在算法框架上也基本一致.只不過式（9）在理論上是準確的，而式（14）代表的射線追蹤算法存在聲學近似引入的誤差.

TTI介質(zhì)與VTI介質(zhì)并沒有物理上的本質(zhì)區(qū)別，若在沿傾斜對稱軸的坐標系下考查TTI介質(zhì)波的傳播問題則與VTI介質(zhì)完全等價.為了推導TTI介質(zhì)的程函方程，將VTI介質(zhì)程函方程進行坐標旋轉(zhuǎn).假設x為標準坐標系，按對稱軸的傾角ν與方位角α旋轉(zhuǎn)后的傾斜坐標系記為x′，則坐標變換對應的雅可比矩陣B可以表示成：

于是由（12）式得到三維TTI介質(zhì)聲學近似意義下的程函方程（見附錄B）.設方位角滿足α=0，則得到其二維形式：

將（16）式改寫為漢密爾頓方程 F（x，y，z，px，py，pz，ν）=0，通過特征值方法仍可將聲學近似意義下TTI介質(zhì)的射線方程表示成（14a）與（14b）的形式.只不過這時對稱軸傾角ν也為空間坐標xi的函數(shù)，因此（16）式對空間坐標xi求導時，還需對傾角ν求導數(shù).

2.3 兩種各向異性射線追蹤算法對比

由于局部角度域疊前深度偏移成像算法既需要射線的走時信息，還需要其方向信息，因此在偏移之前按上述原理編寫射線追蹤算法創(chuàng)建走時與起飛（相）角的數(shù)值表.基于國際上通用的二維TTI推覆體模型，對比兩種算法的計算精度與效率.如圖1（a—d）分別顯示了該模型的垂直qP波速度vP0、Thomsen參數(shù)ε與δ以及對稱軸傾角ν.該模型橫向與縱向采樣點分別為900和200，采樣間隔均為10m.由于聲學近似意義下的射線追蹤方程不是完全精確的，因此本文將基于它計算的射線路徑、走時和文中另一種精確的各向異性介質(zhì)射線追蹤算法進行對比.圖2a與圖2b中紅色射線路徑與走時曲線為相速度表示的射線方程計算結果，而藍色曲線為聲學近似射線方程計算結果.兩種結果非常接近，表明聲學近似射線追蹤算法精度仍然比較高，適用于角度域疊前深度偏移算法中走時與起飛角等相關參數(shù)的計算.

圖1 逆沖模型（a）vP0模型；（b）ε模型；（c）δ模型；（d）ν模型.Fig.1 Numerical example on overthrust model（a）vP0model；（b）εmodel；（c）δmodel；（d）νmodel.

圖2 局部角度域射線追蹤圖（a）射線路徑；（b）走時曲線圖.Fig.2 Local angle domain ray tracing（a）Raypath；（b）Traveltime.

基于同樣的模型，表1列出了兩種算法在VTI與TTI兩種介質(zhì)情況下計算走時表與起飛角度表所需的CPU時間.對于VTI介質(zhì)情況，將原TTI介質(zhì)模型對稱軸傾角取為0.從表中可以看出，VTI介質(zhì)聲學近似射線追蹤算法效率非常高，但當考慮對稱軸傾角時，其計算成本提高了數(shù)倍，這主要是因為從VTI介質(zhì)到TTI介質(zhì)的坐標變換引入了太多的額外計算量.而基于相速度的射線追蹤算法考慮對稱軸傾角增加的計算成本不到10%.

表1 兩種射線方程計算效率對比Table 1 Computational cost comparison between two ray tracing systems

3 地震波局部角度域成像原理

如果從地下成像點的視角來觀察地震波場，成像點處有兩種波場，即入射波與散射波（包括反射與繞射波）.在高頻漸近意義下，入射和散射波前面分別具有各自的格林函數(shù)屬性，如走時、幾何擴散因子、走時梯度或慢度矢量等.如圖3所示，三維情況下，入射慢度矢量ps和散射慢度矢量pr共同描述了散射點m處波的傳播方向特征.入射與散射慢度矢量之和pm稱為照明矢量.根據(jù)地震勘探的需要，可用兩類、四個角度共同定義局部傳播方向[18］.第一類是描述入射與散射（包括繞射和反射）方向特征的兩個角度，即入射角γ（散射張角θ的一半）和散射方位角（即局部入射與散射慢度所在平面的方位角）φ.偏移速度分析、AVA分析／反演就是考察和利用時差或振幅隨這兩個角度的變化.第二類是描述局部照明方向的兩個角度，即照明矢量的傾角?與方位角φ.基于射線理論，這四個角度參數(shù)可由走時的空間梯度計算得到[10］.

圖3 成像點處地震波局部角度特征示意圖Fig.3 Local angle characteristics of a selected ray pair at a subsurface image point

設入射射線的起飛角為βs，方位角為αs，散射射線的起飛角為βr，方位角為αr，很容易得入射射線與散射射線的單位慢度矢量：

根據(jù)矢量運算法則，四個局部角度參數(shù)分別滿足[10］：

式中x、y與z分別代表沿坐標軸的單位矢量，其中y指向正北方向并作為定義方位角的參考方向，pmz為照明矢量的垂向分量.可見，只要根據(jù)起飛角及其方位角計算得到入射與散射慢度矢量，就可根據(jù)上述方程求取四個局部角度參數(shù).

在局部角度域進行射線追蹤時，從地下成像點以起飛角βs（或βr）與方位角αs（或αr）等間隔向上發(fā)射一簇射線到達地表各觀測點，將這些不同方向起飛射線的走時與角度信息保存在數(shù)值表中.在局部角度域成像時，根據(jù)炮點-成像點-接收點關系，在數(shù)值表中讀取計算好的樣本通過插值獲得實際射線路徑的走時、起飛角及其方位角，進而按公式（17a）至（18d）從射線路徑的局部角度參數(shù) （βs，αs；βr，αr）轉(zhuǎn)換成局部角度域成像需要的角度參數(shù)（φ，γ；φ，?）.

根據(jù)Kirchhoff積分偏移原理，地面觀測地震記錄各個時刻振幅所對應的偏移脈沖響應按空間位置疊加起來就得到地下構造圖像.事實上，脈沖響應曲面任意一點都與可能的特定射線路徑相對應，且在各點都具有其局部角度屬性.常規(guī)的成像結果相當于不同傳播方向波場成像值的某種平均，而局部角度域成像就是要在疊前偏移過程中保留這些像的局部方向信息.為了降低計算成本，一般在完全疊加成像數(shù)據(jù)與炮檢距域共成像點道集之外，僅額外輸出入射角度（φ，γ）域或者照明角度（φ，?）域的共成像點道集.入射角度域共成像點道集適用于偏移速度分析、成像振幅隨入射角或方位角變化（AVA／AVAZ）分析.照明角度域共成像點道集可用于提取地層走向與傾角屬性，也可用于繞射波的分離與成像等[10］.

根據(jù)Bleistein等[8］的建議，本文把三維入射角度域Kirchhoff疊前深度偏移成像公式寫成：

其中I（x，γ，φ）為入射角度域成像結果，A（x，xs）與A（x，xr）分別對應入射波與反射波格林函數(shù)中的振幅項，W（x，xr，xr）為與該成像點處的覆蓋次數(shù)成反比的加權系數(shù)，D（x，xr，xr）為濾波后的地震數(shù)據(jù)，其表達式為

式中u（xr，xr，ω）表示地表觀測數(shù)據(jù)，τ（x，xr，xr）為入射射線與散射射線走時之和，K（x，pm，γ，φ）為KMAH參數(shù)，它代表射線經(jīng)過焦散處引起相位反轉(zhuǎn)的次數(shù).振幅系數(shù)A（x，xs）與A（x，xr）以及 KMAH參數(shù)通常需由動力學射線追蹤得到.為降低計算成本，本文數(shù)值算法中的振幅系數(shù)采用了均勻介質(zhì)或橫向均勻介質(zhì)格林函數(shù)振幅來替代.如果不考慮隨入射方位角的變化，可將（19）式改寫為：

其中，I（x，?，φ）為照明角度域成像結果.

4 數(shù)值算例

4.1 SEG／HESS VTI模型

圖4（a、b、c）展示了SEG／HESS二維VTI模型垂直qP波速度、Thomsen參數(shù)ε和δ.模型左邊有一個被各向異性巖層包圍的鹽丘，右側為一個斷層面，鹽丘兩側與斷層較為陡峭.射線走時與局部角度計算采用VTI介質(zhì)聲學近似射線追蹤算法.圖4d為本文VTI介質(zhì)Kirchhoff疊前深度偏移結果，鹽丘邊緣與斷層都得到了很好的聚焦成像；克服了各向同性介質(zhì)Kirchhoff疊前深度偏移（圖4e）中反射歸位不準、繞射未完全收斂、能量不聚焦等問題.圖5a展示了該模型偏移后的平均入射角域共成像點道集，其中入射角范圍為0°～60°.與傳統(tǒng)的偏移距域共成像點道集（圖5b）相比，角度域共成像點道集能量更聚焦，波形拉伸效應的影響也弱得多.類似地，圖5c與圖5d分別展示了忽略各向異性影響的情況下得到的入射角域與偏移距域共成像點道集.可以看出，這兩種成像道集大多數(shù)同相軸均未拉平，都有上翹的現(xiàn)象.借助這些剩余曲率信息便能夠進行各向異性偏移速度分析.

圖4 SEG／HESS VTI模型數(shù)值算例（a）vP0模型；（b）ε模型；（c）δ模型；（d）考慮各向異性偏移結果；（e）忽略各向異性偏移結果.Fig.4 Numerical example on SEG／HESS VTI model（a）vP0model；（b）εmodel；（c）δmodel；（d）Considered the effects of anisotropy；（e）Ignored the effects of anisotropy.

4.2 TTI逆沖模型

下面以國際上通用的逆沖模型（圖1）為例，測試TTI介質(zhì)局部角度域疊前深度偏移成像算法.該模型數(shù)據(jù)正演模擬時底部放有一水平反射界面用于測試上覆各向異性介質(zhì)對該界面反射地震波傳播與成像的影響.射線走時與局部角度計算采用前文介紹的相速度表征的TTI介質(zhì)射線追蹤算法.首先對合成數(shù)據(jù)進行VTI介質(zhì)Kirchhoff疊前深度偏移（圖6a），發(fā)現(xiàn)由于忽略對稱軸傾角給走時計算帶來明顯誤差，進而導致該反射界面未能準確成像.圖6b顯示了TTI介質(zhì)Kirchhoff疊前深度偏移結果，此時底部的水平界面得到了正確成像，觀測孔徑限制引起的未完全收斂與疊加掉的繞射能量也減弱了許多.圖6c為TTI介質(zhì)局部角度域Kirchhoff疊前深度偏移成像獲得的平均入射角域共成像點道集.由于合成數(shù)據(jù)偏移距范圍較小，考慮的入射角范圍為0°～25°.圖6d展示了TTI介質(zhì)局部角度域Kirchhoff疊前深度偏移成像獲得的照明傾角域共成像點道集.可見，偏移后的反射波為似雙曲線形狀，其中能量最強的頂點所對應的照明傾角與反射界面的傾角一致，而偏移后的繞射波則具有明顯不同的形態(tài).關于照明角度域共成像點道集物理含義和實際用途的討論可參見文獻[18］，[9］和[10］.

5 結 論

本文基于改進的射線追蹤算法實現(xiàn)了一種面向TI介質(zhì)的局部角度域疊前深度偏移成像方法.對比研究的兩種改進的射線追蹤方法為局部角度域成像提供射線的走時與局部方向信息，它們均不再采用剛度系數(shù)而是相速度或Thomsen參數(shù)來表示各向異性速度模型.其中聲學近似射線追蹤方法是基于VTI介質(zhì)聲學近似qP波波動方程與程函方程推導的，本文通過坐標旋轉(zhuǎn)將其擴展到了TTI介質(zhì).今后可以照此思路推導諸如正交各向異性等更復雜各向異性介質(zhì)的聲學近似射線方程.從推覆體模型數(shù)值試驗看，聲學近似射線追蹤算法精度與文中另一種準確的射線追蹤算法非常接近，但從計算效率考慮更適合VTI介質(zhì).

根據(jù)地震波在成像點處的局部方向特征，文中采用局部角度域成像空間的脈沖響應疊加，實現(xiàn)了入射角度域和照明角度域Kirchhoff疊前深度偏移成像算法.文中二維VTI與TTI模型數(shù)值試驗表明，該算法能夠給各向異性偏移速度分析提供強有力的支撐.這種局部角度域成像方法的優(yōu)點還在于用到了成像孔徑內(nèi)所有波場數(shù)據(jù)而不是某些按地面炮檢距與方位角分選的數(shù)據(jù)子集，成像振幅更合理地反映了目的層“原位”的隨入射角及其方位變化的帶限反射系數(shù)信息.本文方法與常規(guī)Kirchhoff疊前深度偏移一樣，在處理多波至、焦散等問題是存在缺陷.基于與本文TI介質(zhì)運動學射線追蹤對應的動力學射線追蹤算法，可在本文理論框架下發(fā)展局部角度域高斯束偏移成像方法，進而有效地解決多波至與焦散問題，提高在復雜非均勻各向異性介質(zhì)中的成像精度.這部分工作擬另文介紹.

致 謝 感謝SEG提供文中使用的理論模型數(shù)據(jù).

圖6 逆沖模型偏移結果（a）VTI偏移結果；（b）TTI偏移結果；（c）平均入射角域共成像點道集；（d）照明傾角域共成像點道集.Fig.6 TTI data migration results（a）VTI migration；（b）TTI migration；（c）Average incident-angle CIGs；（d）Illumination-dip domain CIGs.

附錄A 相速度表示的TI介質(zhì)射線追蹤方程

前文給出了Zhu等[15-16］提出的相速度表示的一般各向異性介質(zhì)的射線方程.在VTI與TTI介質(zhì)中的具體形式，根據(jù)前文，有如下各向異性射線追蹤方程組：

其中群速度可由相速度表示[19］，

而VTI介質(zhì)相速度滿足[19］：

當TI介質(zhì)對稱軸與垂向存在夾角為ν時，相速度滿足：

于是在（A1）中采用這種適用于TTI介質(zhì)的相速度公式即可.如果令相速度公式中的f=1，上述射線方程也就退化成聲學近似形式.

附錄B：TTI介質(zhì)聲學近似程函方程

TTI介質(zhì)聲學近似意義下的程函方程有如下形式：

其中，

（References）

[1]Kumar D，Sen M K，F(xiàn)erguson R J.Traveltime calculation and prestack depth migration in tilted transversely isotropic media.Geophysics，2004，69（1）：37-44.

[2]Zhu T F，Gray S H，Wang D L.Prestack Gaussian-beam depth migration in anisotropic media.Geophysics，2007，72（3）：S133-S13.

[3]吳國忱.各向異性介質(zhì)地震波傳播與成像.東營：中國石油大學出版社，2005.Wu G C.Seismic Wave Propagation and Imaging in Anisotropic Medium （in Chinese）.Dongying：China Petroleum University Press，2005.

[4]康瑋，程玖兵.橫向各向同性介質(zhì)擬聲波方程及其在逆時偏移中的應用.地球物理學報，2012，55（3）：1033-1045.Kang W，Cheng J B.Pseudo-acoustic wave equations for reverse-time migration in TI media.Chinese J.Geophys.（in Chinese），2012，55（3）：1033-1045.

[5]de Hoop M V，Bleistein N.Generalized radon transform inversions for reflectivity in anisotropic elastic media.Inverse Problems，1997，13（3）：669-690.

[6]Xu S，Chauris H，Lambaré G，et al.Common-angle migration：A strategy for imaging complex media.Geophysics，2001，66（6）：1877-1894.

[7]Brandsberg-Dahl S，Ursin B，de Hoop M V.Seismic velocity analysis in the scattering angle／azimuth domain.Geophysical Prospecting，2003，51（4）：295-314.

[8]Bleistein N，Gray S H.A proposal for common-openingangle migration／inversion.Center for Wave Phenomena，Colorado School of Mines.Research Report CWP-420，2002.

[9]Koren Z，Ravve I，Ragoza E，et al.Full-azimuth angle domain imaging.SEG Technical Program Expanded Abstracts，2008：2221-2225.

[10]Cheng J B，Wang T F，Wang C L，et al.Azimuth-preserved local angle-domain prestack time migration in isotropic，vertical transversely isotropic and azimuthally anisotropic media.Geophysics，2012，77（2）：S51-S64.

[11]Vidale J E.Finite-difference calculation of traveltimes in three dimensions.Geophysics，1990，55（5）：521-526.

[12]Schneider W A，Ranzinger K A，Balch A H，et al.A dynamic programming approach to first-arrival traveitlme computation in media with arbitrarily distributed velocities.Geophysics，1992，57（1）：39-50.

[13]Farra V.Ray tracing in complex media.Journal of Applied Geophysics，1993，30（1-2）：55-73.

[14]Cerveny V.Seismic Ray Theory.Cambridge：Cambridge University Press，2001.

[15]Zhu T F，Gray S H，Wang D.Kinematic and dynamic raytracing in anisotropic media-aphase-velocity formulation.Theory and application：75th Annual International Meeting，SEG，ExpandedAbstracts，2005，96-99.

[16]Zhu T F，Gray S H，Wang D L.Prestack gaussian-beam depth migration in anisotropic media.Geophysics，2007，72（3）：S133-S138.

[17]Alkhalifah T.An acoustic wave equation for anisotropic media.Geophysics，2000，65（4）：1239-1250.

[18]Audebert F，F(xiàn)roidevaux P，Racotoarisoa H，et al.Insights into migration in the angle domain.SEG Technical Program Expanded Abstracts，2002：1188-1191.

[19]Tsvankin I.Seismic Signatures and Analysis of Reflection Data in Anisotropic Media.New York：Elsevier，2001.