李興旺 白超英*② 李曉玲

（①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系，陜西西安710054；②長(zhǎng)安大學(xué)計(jì)算地球物理研究所，陜西西安710054；③國(guó)家知識(shí)產(chǎn)權(quán)專利局專利審查協(xié)作河南中心，河南鄭州450000）

·地震模擬·

四面體單元剖分下三維各向異性TI介質(zhì)中多次波射線追蹤

李興旺①白超英*①②李曉玲③

（①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系，陜西西安710054；②長(zhǎng)安大學(xué)計(jì)算地球物理研究所，陜西西安710054；③國(guó)家知識(shí)產(chǎn)權(quán)專利局專利審查協(xié)作河南中心，河南鄭州450000）

采用不規(guī)則四面體單元做模型剖分，實(shí)現(xiàn)了三維復(fù)雜層狀各向異性TI介質(zhì)中多次（透射、反射及轉(zhuǎn)換）波的追蹤計(jì)算，彌補(bǔ)了現(xiàn)行算法主要采用規(guī)則網(wǎng)格計(jì)算多次波的不足。根據(jù)介質(zhì)的5個(gè)彈性參數(shù)及其對(duì)稱軸傾角和走向，計(jì)算主節(jié)點(diǎn)上qP、qSV和qSH波群速度，由此線性插值得到次級(jí)節(jié)點(diǎn)的群速度值；然后采用不規(guī)則最短路徑算法，并結(jié)合分區(qū)多步計(jì)算技術(shù)實(shí)現(xiàn)多次波追蹤計(jì)算。均勻各向異性介質(zhì)中計(jì)算出的數(shù)值解與解析解吻合，驗(yàn)證了算法的有效性。復(fù)雜層狀模型數(shù)值模擬結(jié)果表明：采用不規(guī)則四面體單元可有效擬合起伏地表及地下不規(guī)則速度間斷面，同時(shí)該算法適用于任意傾角的TI介質(zhì)，且不受各向異性強(qiáng)度限制。

三維各向異性TI介質(zhì) 改進(jìn)型最短路徑算法 多次波射線追蹤 四面體單元 分區(qū)多步計(jì)算技術(shù)

1 引言

地球內(nèi)部地層普遍存在各向異性現(xiàn)象，地震波在其中傳播時(shí)，運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)屬性均會(huì)發(fā)生明顯改變，最為突出的是地震波速度的各向異性，即地震波速度不但與空間位置有關(guān)，而且與傳播方向有關(guān)。當(dāng)?shù)貙拥母飨虍愋圆豢珊雎詴r(shí)，采用基于各向同性的方法來(lái)處理地震數(shù)據(jù)，將帶來(lái)較大的計(jì)算誤差，降低資料的解釋精度，嚴(yán)重時(shí)可能得到錯(cuò)誤的解釋。因此，研究各向異性介質(zhì)中地震波的傳播規(guī)律顯得尤為重要。早期研究大多針對(duì)具有垂直對(duì)稱軸的VTI介質(zhì)以及水平對(duì)稱軸的HTI介質(zhì)［1-5］。然而，實(shí)際地層更為復(fù)雜，例如褶皺和上沖等作用使地層發(fā)生傾斜、變形，介質(zhì)對(duì)稱軸也將發(fā)生變化，即不再是垂直或水平。此時(shí)，簡(jiǎn)單的VTI或HTI介質(zhì)的假設(shè)難以滿足地震勘探的實(shí)際要求。因此，有必要研究具有傾斜對(duì)稱軸的TTI介質(zhì)模型。

以往各向異性介質(zhì)中的射線追蹤，主要是在弱各向異性的假設(shè)條件下，采用傳統(tǒng)算法（如打靶法或彎曲法）進(jìn)行初至波追蹤計(jì)算［6-11］。若介質(zhì)較為復(fù)雜或討論強(qiáng)各向異性介質(zhì)時(shí)，上述方法將不再適用（出現(xiàn)“偏靶”或陷入局部解）。近年來(lái)出現(xiàn)的基于網(wǎng)格（或單元）波前擴(kuò)展的算法（有限差分解程函方程法、最短路徑算法、波前構(gòu)造法等），由于其計(jì)算的穩(wěn)健性、全局解、精度高、耗時(shí)少等優(yōu)點(diǎn)，得到了廣泛應(yīng)用［12-17］。但此類研究依然局限于對(duì)初至波的討論。另外，多數(shù)算法采用規(guī)則網(wǎng)格（或單元）進(jìn)行模型參數(shù)化，不易處理起伏地表和地下復(fù)雜波阻抗界面，所以當(dāng)模型復(fù)雜時(shí)存在較大計(jì)算誤差。為了彌補(bǔ)以上缺陷，同時(shí)討論一般各向異性介質(zhì)中多次波的追蹤計(jì)算，已實(shí)現(xiàn)分區(qū)多步不規(guī)則最短路徑算法下的2D／3D各向異性TTI介質(zhì)中多次波的追蹤計(jì)算（簡(jiǎn)稱 Multistage ISPM_TTI）［18］，以及三角網(wǎng)格剖分下二維各向異性TTI介質(zhì)中多次波的追蹤計(jì)算（簡(jiǎn)稱 Multistage TSPM_TTI）［19］。本文則是將四面體單元剖分下各向同性介質(zhì)中的多次波計(jì)算方法（簡(jiǎn)稱 Multistage TSPM_3D）［20，21］推廣至一般各向異性TTI介質(zhì)中，從而實(shí)現(xiàn)含起伏地表（包括地下不規(guī)則界面）三維各向異性TTI介質(zhì)中多次波地震射線的追蹤計(jì)算。該算法采用不規(guī)則四面體單元做模型剖分，能更精確地?cái)M合實(shí)際地學(xué)模型，適用于一般各向異性介質(zhì)，主要難點(diǎn)是網(wǎng)格剖分和TTI介質(zhì)各方向上相、群速度的計(jì)算，包括qSV波在奇異點(diǎn)出現(xiàn)多值時(shí)的處置。

2 TTI介質(zhì)中相速度和群速度的計(jì)算

各向異性介質(zhì)中相、群速度發(fā)生分離（圖1），射線沿群速度方向傳播，故將Multistage TSPM_3D算法推廣到各向異性介質(zhì)時(shí)，關(guān)鍵步驟是計(jì)算群速度。

圖1 三維各向異性TTI介質(zhì)相、群速度示意圖

為了計(jì)算TTI介質(zhì)的相、群速度，通常的做法是首先得到VTI介質(zhì)中相速度的精確表達(dá)式，然后再經(jīng)過(guò)坐標(biāo)旋轉(zhuǎn)得到TTI介質(zhì)的相、群速度計(jì)算式。Daley等［22］給出了VTI介質(zhì)中相速度的精確表達(dá)式，它僅依賴于以下5個(gè)彈性參數(shù)

式中：c1、c2、c3分別表示qP、qSV 和qSH 波的相速度；υ為入射角；（a11，a13，a33，a55，a66）是介質(zhì)的五個(gè)彈性參數(shù)，且隨位置坐標(biāo)變化。

根據(jù)相、群速度的關(guān)系［23］（下標(biāo)m=1，2，3分別對(duì)應(yīng)三種類型的波，即qP、qSV和qSH），可得

上述討論得到了VTI介質(zhì)中群速度的計(jì)算方法，為了計(jì)算有任意傾斜對(duì)稱軸的TTI介質(zhì)，需將介質(zhì)對(duì)稱軸旋轉(zhuǎn)一定角度（θ0，φ0）。對(duì)于某特定的相速度角（θ，φ），存在

將其代入式（1）～式（3），即可得到TTI介質(zhì)中群速度值，其中相位的傳播方向?yàn)?/p>

射線的傳播方向是

另外，特殊（奇異值）點(diǎn)上的qSV波會(huì)出現(xiàn)三重值，此時(shí)宜選擇群速度值中的最小值進(jìn)行計(jì)算［14］。

3 四面體單元剖分下的最短路徑算法原理

3.1 散點(diǎn)集四面體單元網(wǎng)的生成

近十幾年來(lái)，散點(diǎn)四面體剖分一直是眾多學(xué)者研究和關(guān)注的焦點(diǎn)。到目前為止出現(xiàn)了不少算法：Dwyer［24］提 出 Voronoi圖 法；Radovitzky 等［25］則采用逐點(diǎn)插入法；郭際元等［26］給出了三種建立四面體網(wǎng)格的算法；孔娟華等［27］則是在Delaunay四面體的基礎(chǔ)上提出了遞推生成算法。本文將二維三角網(wǎng)格生長(zhǎng)算法［20］推廣到三維情形，形成的四面體滿足空球準(zhǔn)則，很好地解決了四面體不相容問(wèn)題。圖2給出了四面體網(wǎng)格剖分的示意圖，為了成圖清晰，僅給出模型3個(gè)表面（未標(biāo)出次級(jí)節(jié)點(diǎn)）的分布情況。

3.2 模型參數(shù)化及界面離散

對(duì)于三維模型，按照各區(qū)界面的起伏情況采用不規(guī)則四面體單元對(duì)模型進(jìn)行網(wǎng)格剖分（圖2），四面體單元的4個(gè)角點(diǎn)為主節(jié)點(diǎn)（圖3a），四個(gè)面上插入等間距的次級(jí)節(jié)點(diǎn)（圖3b），四面體內(nèi)部無(wú)節(jié)點(diǎn)（炮點(diǎn)和檢波器可位于其內(nèi)）。根據(jù)速度界面的具體情況，將其離散化，界面上離散點(diǎn)的間距應(yīng)與四面體單元面上次級(jí)節(jié)點(diǎn)的間距相當(dāng)。速度采樣在主節(jié)點(diǎn)上，次級(jí)節(jié)點(diǎn)和界面離散點(diǎn)的速度值可通過(guò)線性插值獲得［20，21］

圖2 含起伏界面的三維層狀模型中反射、透射、轉(zhuǎn)換波計(jì)算示意圖

圖3 四面體內(nèi)速度插值示意圖

式中：v j表示i節(jié)點(diǎn)所在單元各主節(jié)點(diǎn)的速度值；uj為體積坐標(biāo)。圖3a中P點(diǎn)在四面體T1T2T3T4內(nèi)的體積坐標(biāo)u1可定義為：u1=V PT2T3T4／V T1T2T3T4，其中V PT2T3T4、V T1T2T3T4分別表示四面體PT2T3T4和T1T2T3T4的體積，對(duì)于u2、u3、u4，依此類推。

3.3 波前計(jì)算和射線追蹤原理

分區(qū)多步計(jì)算方法是追蹤多次波的核心。結(jié)合模型的起伏地表及波阻抗界面的幾何形態(tài)、位置，將模型劃分成多個(gè)相對(duì)獨(dú)立的層狀起伏區(qū)域。在射線追蹤過(guò)程中，可根據(jù)多次波的類型進(jìn)行分區(qū)多步計(jì)算，即確定每一步計(jì)算時(shí)的計(jì)算區(qū)域和調(diào)用相應(yīng)的速度模型。由于每一步計(jì)算僅針對(duì)確定區(qū)域，規(guī)避了“無(wú)用區(qū)域”（射線不需要通過(guò)區(qū)域），在很大程度上減少計(jì)算時(shí)間和內(nèi)存。對(duì)于大型三維模型，該優(yōu)勢(shì)尤為突出。

具體來(lái)講，自震源所在的單元開(kāi)始計(jì)算各節(jié)點(diǎn)的最小旅行時(shí)，該單元所有節(jié)點(diǎn)的旅行時(shí)計(jì)算結(jié)束后，進(jìn)行波前擴(kuò)展，直到該區(qū)域內(nèi)所有節(jié)點(diǎn)均獲得最小旅行時(shí)。此時(shí)波前停留在該區(qū)的速度離散界面上，同時(shí)保存速度界面上各離散點(diǎn)的旅行時(shí)和相應(yīng)的射線路徑。然后按照所要追蹤的射線類型選擇下一個(gè)區(qū)域進(jìn)行計(jì)算，直至追蹤到檢波器為止。圖2所示的多次波射線路徑（紅線）為P1SV2P2SV2P2P1（這里規(guī)定上標(biāo)代表上行波，下標(biāo)代表下行波，數(shù)字1、2表示分區(qū)號(hào)，另外，為了方便起見(jiàn)，這里P表示qP；SV表示qSV；SH表示qSH），采用上述原理可獲得更為復(fù)雜多次波的射線路徑和相應(yīng)旅行時(shí)，技術(shù)細(xì)節(jié)可參見(jiàn)文獻(xiàn)［18-20］。

3.4 節(jié)點(diǎn)旅行時(shí)計(jì)算方法

由于單元網(wǎng)格間距較小，相鄰點(diǎn)間速度差異不大，對(duì)于i點(diǎn)附近任意節(jié)點(diǎn)j處的旅行時(shí)（可僅限于相同單元，也可僅限于相鄰單元，參見(jiàn)文獻(xiàn)［20，21］，即采用一階或二階Forward Star計(jì)算技術(shù)），可由以下公式（兩點(diǎn)間距離除以平均群速度）近似計(jì)算

式中：i是波前面上當(dāng)前的次級(jí)震源；j為將要計(jì)算的節(jié)點(diǎn)；Cj是j鄰近節(jié)點(diǎn)的集合；是自炮點(diǎn)到達(dá)i節(jié)點(diǎn)射線的最小旅行時(shí)，其中m=1、2、3分別對(duì)應(yīng)qP、qSV和qSH 三類波；D（X i，X j）是節(jié)點(diǎn)i與節(jié)點(diǎn)j之間的距離；和分別是i和j節(jié)點(diǎn)的群速度值。

當(dāng)模型參數(shù)（a11，a13，a33，a55，a66，θ0，φ0）不隨空間坐標(biāo)變化時(shí)，即為均勻各向異性介質(zhì)，U m只是（θ′，φ′）的函數(shù)。由公式（5）可得到等時(shí)面解析解，即Δt時(shí)刻沿（θ′，φ′）方向的一點(diǎn)到坐標(biāo)原點(diǎn)的距離R（θ′，φ′），可確定該時(shí)刻等時(shí)面上的一個(gè)波前點(diǎn)（θ′，φ′，R（θ′，φ′））。 當(dāng)θ′、φ′分 別 取 遍 ［0°，180°］、［0°，360°］各值時(shí)，可以得到完整的 Δt時(shí)刻的等時(shí)面（解析解）

4 實(shí)例模擬和結(jié)果分析

4.1 群速度與相慢度角和射線角的關(guān)系

為直觀顯示TI介質(zhì)中群速度隨相慢度角、射線角的變化情況，選擇一組彈性參數(shù)（a11=5.2，a13=0.93，a33=4.0，a55=1.0，a66=1.0），分別針對(duì)VTI（θ0=0°），HTI（θ0=90°），TTI（θ0=45°）介質(zhì)進(jìn)行了試算。圖4左圖顯示群速度隨相慢度角的變化情況，右圖為群速度隨射線角的變化情況（從上到下θ0依次為0°、45°、90°）。從圖4可見(jiàn)，qSV波群速度隨射線角變化時(shí)出現(xiàn)了多重值現(xiàn)象，即同一射線角對(duì)應(yīng)多個(gè)群速度值，這是Multistage TSPM_3D算法不能處理的，故應(yīng)選擇最小速度值進(jìn)行計(jì)算以保證速度的連續(xù)性。而群速度與相慢度角的關(guān)系則是單一的，并沒(méi)有出現(xiàn)多重值現(xiàn)象。同時(shí)，注意到qSH波群速度為常數(shù)，這是因?yàn)閍55=a66=1.0。

圖4 對(duì)稱軸角度不同時(shí)三種波（qP、qSV和qSH）的群速度隨相慢度角（左）及射線角（右）變化曲線

4.2 TTI介質(zhì)群速度分布

采用Thomsen［28］給出的一組TTI模型參數(shù)（a11=25.7，a13=15.2，a33=15.4，a55=4.2，a66=9.0，θ0=45°和φ0=315°），計(jì)算得到t=1.0s時(shí)qP、qSV、qSH三種波等時(shí)波前的解析解（或稱群速度分布圖，圖5a～圖5c）和數(shù)值解（圖5d～圖5f）。從圖5中可以看出，三類波群速度分布與該點(diǎn)和對(duì)稱軸之間的夾角有關(guān)，在垂直于對(duì)稱軸的平面內(nèi)均為常量，但由于介質(zhì)的各向異性，三類波都有著各自獨(dú)特的波前面或群速度空間分布形態(tài)，尤其是qSV波的解析波前發(fā)生了分裂，沿對(duì)稱軸方向及其垂向均出現(xiàn)了三重值（或重疊）現(xiàn)象。由于數(shù)值解僅使用qSV波最小速度進(jìn)行計(jì)算，其結(jié)果并無(wú)波前分裂及三重值現(xiàn)象發(fā)生。

圖5 TTI介質(zhì)中qP、qSV、qSH波（從左到右）在t=1.0s時(shí)刻的等時(shí)波前

4.3 計(jì)算精度分析

為驗(yàn)證算法的有效性，選取Thomsen的Clay Shale模型（a11=15.1，a13=1.6，a33=10.8，a44=3.1，a66=4.3）構(gòu)建均勻各向異性模型。模型尺寸為10km×10km×5km，炮點(diǎn)位于x軸中間點(diǎn)（5km，0，0）處。圖6左右兩部分分別給出了（θ0=0°，φ0=0°）和（θ0=45°，φ0=0°）時(shí)，x-z和x-y剖面上qP、qSV和qSH的等時(shí)波前，圖中黑色為解析解，紅色為數(shù)值解。為清晰起見(jiàn)，圖中數(shù)值解時(shí)間間隔是解析解的一半?？梢?jiàn)因群速度差異，三類波都有各自獨(dú)特的波前面，且波前隨介質(zhì)對(duì)稱軸傾斜方向不同而變化，亦證明群速度大小與介質(zhì)對(duì)稱軸的傾角有關(guān)。同時(shí)還發(fā)現(xiàn)qSV數(shù)值解中未出現(xiàn)解析解中的蝴蝶結(jié)（多重值）現(xiàn)象，這是由于數(shù)值算法無(wú)法使用多個(gè)群速度值來(lái)進(jìn)行射線追蹤造成的（只取了多重值中的最小群速度值進(jìn)行計(jì)算）。除了蝴蝶結(jié)處外，其余部分?jǐn)?shù)值解與解析解相吻合，甚至在蝴蝶結(jié)附近也非常一致。

4.4 復(fù)雜TTI介質(zhì)中多次波追蹤

本文所提算法可用于復(fù)雜TTI介質(zhì)中多次波追蹤，為使展示結(jié)果清晰明了，數(shù)值模擬時(shí)采用層狀模型。選擇如圖7所示的三維層狀起伏模型：模型尺寸為80km×80km×50km；地表起伏，最大起伏幅度為2km；在25km深度處有一個(gè)凹陷反射界面，最大凹陷幅度為4km；底部?jī)A斜反射界面位于47～50km深度處。炮點(diǎn)位于（40km，40km，10km）的位置（圖7中黑色圓點(diǎn)），31個(gè)檢波器以等間距（8.0km）分布在模型地表的三條邊上（圖7中灰色圓點(diǎn)）。第一層模型參數(shù)為：a11=15.1，a13=1.6，a33=10.8，a55=3.1，a66=4.3，θ=φ=0°第二層模型參數(shù)為：a11=6.3，a13=2.25，a33=4.51，a55=1.0，a66=1.5，θ0=45°，φ0=315°。圖7a展示了一次反射的情形，震相分別為P1P1、P1P2P2P1，即自炮點(diǎn)激發(fā)的P波經(jīng)第一、二層起伏界面的反射波。圖7b給出了多次透射、反射及轉(zhuǎn)換波的射線路徑：第一個(gè)震相為P1SV1P1，炮點(diǎn)激發(fā)的P波上行至地表，經(jīng)反射并轉(zhuǎn)換成SV波后下行至第一層反射界面，再次反射并轉(zhuǎn)換成P波上行至檢波器；第二個(gè)震相為P1SV1P2SV2P1，炮點(diǎn)激發(fā)的P波上行至地表，經(jīng)反射轉(zhuǎn)換為SV波后下行至第一層反射界面，透射且轉(zhuǎn)換成P波后繼續(xù)下行至第二層反射界面，然后反射轉(zhuǎn)換成SV波上行至第一層反射界面，最后透射轉(zhuǎn)換成P波上行至檢波器。

圖6 基于三維Clay Shale模型的qP、qSV、qSH波在x-z和x-y剖面上的等時(shí)波前

圖7 一次反射波P1 P1（藍(lán)線）和P1 P2 P2 P1（紅線）射線路徑（a）、多次反射（轉(zhuǎn)換）震相P1 SV1 P1（藍(lán)線）和P1 SV1 P2 SV2 P1（紅線）射線路徑（b）

5 結(jié)論

本文將用于三維各向同性介質(zhì)中的Multistage TSPM_3D算法推廣至三維各向異性介質(zhì)中，實(shí)現(xiàn)了三維TTI介質(zhì)中多次（透射、反射、轉(zhuǎn)換）波的追蹤計(jì)算（簡(jiǎn)稱：Multistage TSPM_3D_TTI算法）。該算法適用于任意傾角的TTI介質(zhì)，且不受各向異性強(qiáng)度（即無(wú)弱各向異性的假設(shè)）的限制，可應(yīng)對(duì)各類復(fù)雜介質(zhì)；由于模型剖分時(shí)采用了不規(guī)則四面體單元，對(duì)起伏地表和速度間斷面的擬合更加精確，能夠處理各種復(fù)雜層狀模型；同時(shí)，分區(qū)多步計(jì)算技術(shù)的采用，可實(shí)現(xiàn)復(fù)雜多次波的追蹤計(jì)算。數(shù)值模擬驗(yàn)證了該算法的有效性和準(zhǔn)確性，表明該算法具有較高的計(jì)算精度。下一步研究的重點(diǎn)則是利用這些多震相旅行時(shí)數(shù)據(jù)進(jìn)行各向異性參數(shù)（其中包括傾斜對(duì)稱軸）的反演。

［1］ ?erveny V.Seismic rays and rays intensities in inhomogeneous anisotropic media.Geophysical Journal of Royal Astronomic Society，1972，29（1）：1-13.

［2］ Babich V M.Ray method of calculating the intensity of wavefronts in the case of a heterogeneous，anisotropic，elastic medium.Geophysical Journal International，1994，118（2）：379-383.

［3］ Cardarelli E，Cerreto A.Ray tracing in elliptical anisotropic media using the linear traveltime interpolation method applied to traveltime seismic tomography.Geophysical Prospecting，2002，50（1）：55-72.

［4］ Qian J，Symes W W.Finite-difference quasi-P traveltimes for anisotropic media.Geophysics，2002，67（1）：147-155.

［5］ 馬德堂，朱光明.橫向各向同性介質(zhì)中初至波的旅行時(shí)計(jì)算.石油地球物理勘探，2006，41（1）：26-31.Ma Detang，Zhu Guangming.Calculating traveltime of the first arrival in transversely isotropic medium.OGP，2006，41（1）：26-31.

［6］ Grechka V Y，Mc Mechan G A.3-D tow-point ray tracing for heterogeneous，weakly transversely isotropic media.Geophysics，1996，61（6）：1883-1894.

［7］ Psencik I，F(xiàn)arra V.First-order ray tracing for qP waves in inhomogeneous weakly anisotropic Media.Geophysics，2005，70（6）：D65-D75.

［8］ ?erveny V，F(xiàn)irbas P.Numerical modelling and inversion of traveltimes of seismic body waves in inhomogeneous anisotropic media.Geophysical Journal Royal Astronomical Society，1984，76（1）：41-51.

［9］ Gajewski D，Psencik I.Computation of high-frequency seismic wavefields in 3-D laterally inhomogeneous anisotropic media.Geophysical Journal Royal Astronomical Society，1987，91（2）：383-411.

［10］ Shearer P M，Chapman C H.Ray tracing in anisotropic media with a linear gradient.Geophysical Journal International，1988，94（3）：575-580.

［11］ Farm V.First-order ray tracing for qS waves inhomogeneous weakly anisotropic media.Geophysical Journal International，2005，161（2）：309-324.

［12］ 張中杰，何樵登.各向異性介質(zhì)中折射波及其走時(shí)曲線研究.西安石油學(xué)院學(xué)報(bào)，1992，7（3）：1-6.Zhang Zhongjie and He Qiaoden.The study of refracted wave and its traveltime curve in anisotropic medium.Journal of Xi’an Petroleum Institute，1992，7（3）：1-6.

［13］ 趙愛(ài)華，張美根，丁志峰.橫向各向同性介質(zhì)中地震走時(shí)模擬.地球物理學(xué)報(bào)，2006，49（6）：1762-1769.Zhao Aihua，Zhang Meigen，Ding Zhifeng.Seismic traveltime computation for transversely isotropic media.Chinese Journal of Geophysics，2006，49（6）：1762-1769.

［14］ Zhou B，Greenhalgh S.On the computation of elastic wave group velocities for a general anisotropic medium.Journal of Geophysics and Engineering，2004，1（3）：205-215.

［15］ Zhou B，Greenhalgh S.Shortest path＇ray tracing for most general 2D／3D anisotropic media.Journal of Geophysics and Engineering，2005，2（1）：54-63.

［16］ Zhou B，Greenhalgh S.Raypath and traveltime computations for 2D transversely isotropic media with dipping symmetry axes.Exploration Geophysics，2006，37（2）：150-159.

［17］ 馬德堂，朱光明，范廷恩.二維TTI介質(zhì)中初至波旅行時(shí)的搜索算法.石油地球物理勘探，2011，46（5）：710-714.Ma Detang，Zhu Guangming，F(xiàn)an Ting’en.Search algorithm of first arrival traveltime in 2D TTI medium.OGP，2011，46（5）：710-714.

［18］ Bai C Y，Huang G J，Li X L et al.Ray tracing of multiply transmitted／reflected／converted waves in 2D／3D layered anisotropic TTI media and application to crosswell traveltime tomography.Geophysical Journal International，2013，195（2）：1068-1087.

［19］ 李曉玲，白超英，胡廣義.起伏層狀TI介質(zhì)中多次波射線追蹤.石油地球物理勘探，2013，48（6）：924-931.Li Xiaoling，Bai Chaoying，Hu Guangyi.Multiple ray tracing in an undulated layered TI media.OGP，2013，48（6）：924-931.

［20］ Bai C Y，Li X L，Tang X P.Seismic wavefront evolution of multiply reflected，transmitted，and converted phases in 2D／3D triangular cell model.Journal of Seismology，2011，15（4）：637-652.

［21］ Bai C Y，Li X，Wang Q L，Peng J B.Multiple arrival tracking within irregular triangular or tetrahedral cell model.Journal of Geophysics and Engineering，2012，9（1）：29-38.

［22］ Daley P，Hron F.Reflection and transmission coefficients for transversely isotropic media.Bulletin of the Seismological Society of America，1977，67（3）：661-675.

［23］ Berryman J G.Long-wave elastic anisotropy in transversely isotropic media.Geophysics，1979，44（5）：896-917.

［24］ Dwyer R A.Higher-dimensional voronoi diagram in linear expected time.Discrete and Computational Geometry，1991，6（4）：64-67.

［25］ Radovitzky R，Ortiz M.Tetrahedral mesh generation based on node insertion in crystal lattice arrangements and advancing-front-Delaunay triangulation.Computer Methods in Applied Mechanics and Engineering，2000，187（3-4）：543-569.

［26］ 郭際元，龔俊芳.由三維離散數(shù)據(jù)生成四面體網(wǎng)格算法研究.中國(guó)地質(zhì)大學(xué)學(xué)報(bào)（地球科學(xué)），2002，27（3）：271-273.Guo Jiyuan，Gong Junfang.Algorithms of producing tetrahedral network from three dimensional dispersed data.Journal of China University of Geosciences（Earth Science），2002，27（3）：271-273.

［27］ 孔娟華，鄭江濱.一種三維離散點(diǎn)數(shù)據(jù)生成結(jié)構(gòu)四面體算法.計(jì)算機(jī)工程與科學(xué)，2009，31（1）：35-37.Kong Juanhua，Zheng Jiangbin.An algorithm of generating unstructured tetrahedrons from 3D discrete points.Computer Engineering and Science，2009，31（1）：35-37.

［28］ Thomsen L.Weak elastic anisotropy.Geophysics，1986，51（10）：1954-1966.

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.008

李興旺，白超英，李曉玲.四面體單元剖分下三維各向異性TI介質(zhì)中多次波射線追蹤.石油地球物理勘探，2017，52（1）：48-55.

1000-7210（2017）01-0048-08

*陜西省西安市長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系，710054。Email：baicy＠chd.edu.cn

本文于2015年11月10日收到，最終修改稿于2016年12月12日收到。

本項(xiàng)研究受國(guó)家重大科技專項(xiàng)子課題“海上斜井井間地震資料成像處理技術(shù)及應(yīng)用研究”（2011ZX05024-001-03）、長(zhǎng)安大學(xué)中央高?；究蒲袠I(yè)務(wù)費(fèi)課題（310826150006）資助。

（本文編輯：朱漢東）

李興旺 博士研究生，1989年生；2012年畢業(yè)于長(zhǎng)安大學(xué)地球物理學(xué)專業(yè)，獲理學(xué)學(xué)士學(xué)位；2014年畢業(yè)于該校固體地球物理學(xué)專業(yè)，獲理學(xué)碩士學(xué)位；現(xiàn)在該校攻讀地球探測(cè)與信息技術(shù)專業(yè)博士學(xué)位，近期主要致力于各向異性介質(zhì)中地震射線追蹤及菲涅耳體旅行時(shí)層析成像研究。