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

1 同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院，上海 200092

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

1 引 言

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

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

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

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

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

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

各向異性介質(zhì)中的運(yùn)動(dòng)學(xué)射線(xiàn)方程最早由Cerveny[14］給出.下文簡(jiǎn)單回顧其推導(dǎo)過(guò)程.一般各向異性介質(zhì)無(wú)源的彈性波動(dòng)方程可以表示為

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

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

其中，Γjk=aijklpipl為Christoffel矩陣，aijkl=cijkl／ρ為密度歸一化剛度系數(shù)，pi=?τ／?xi為慢度矢量的各個(gè)分量.方程（3）對(duì)應(yīng)一個(gè)標(biāo)準(zhǔn)的特征值問(wèn)題，且特征值滿(mǎn)足G（pi，xi）=1.（3）式可改寫(xiě)為

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

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

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

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

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

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

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

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

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

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

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

聲學(xué)近似就是假設(shè)沿對(duì)稱(chēng)軸方向qSV波的傳播速度為0，即VS0=0，這樣就可將原始的VTI介質(zhì)彈性波動(dòng)方程及其頻散關(guān)系簡(jiǎn)化.假設(shè)地下介質(zhì)為聲學(xué)介質(zhì)，由VTI介質(zhì)彈性波動(dòng)方程及其頻散關(guān)系可推導(dǎo)出近似的qP波標(biāo)量波動(dòng)方程，進(jìn)而得到相應(yīng)的程函方程和射線(xiàn)方程[17］.研究表明，聲學(xué)近似對(duì)qP波運(yùn)動(dòng)學(xué)特征的負(fù)面影響基本可以忽略[17］.

根據(jù)VTI介質(zhì)qP波的頻散關(guān)系，聲學(xué)近似qP波波動(dòng)方程滿(mǎn)足[17］：

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

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

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

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

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

2.3 兩種各向異性射線(xiàn)追蹤算法對(duì)比

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

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

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

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

表1 兩種射線(xiàn)方程計(jì)算效率對(duì)比Table 1 Computational cost comparison between two ray tracing systems

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

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

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

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

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

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

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

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

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

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

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

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

4 數(shù)值算例

4.1 SEG／HESS VTI模型

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

圖4 SEG／HESS VTI模型數(shù)值算例（a）vP0模型；（b）ε模型；（c）δ模型；（d）考慮各向異性偏移結(jié)果；（e）忽略各向異性偏移結(jié)果.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逆沖模型

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

5 結(jié) 論

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

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

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

圖6 逆沖模型偏移結(jié)果（a）VTI偏移結(jié)果；（b）TTI偏移結(jié)果；（c）平均入射角域共成像點(diǎn)道集；（d）照明傾角域共成像點(diǎn)道集.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ì)射線(xiàn)追蹤方程

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

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

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

當(dāng)TI介質(zhì)對(duì)稱(chēng)軸與垂向存在夾角為ν時(shí)，相速度滿(mǎn)足：

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

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

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

其中，

（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]吳國(guó)忱.各向異性介質(zhì)地震波傳播與成像.東營(yíng)：中國(guó)石油大學(xué)出版社，2005.Wu G C.Seismic Wave Propagation and Imaging in Anisotropic Medium （in Chinese）.Dongying：China Petroleum University Press，2005.

[4]康瑋，程玖兵.橫向各向同性介質(zhì)擬聲波方程及其在逆時(shí)偏移中的應(yīng)用.地球物理學(xué)報(bào)，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.