李興旺 白超英*② 李曉玲

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

·地震模擬·

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

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

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

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

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

1 引言

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

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

2 TTI介質中相速度和群速度的計算

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

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

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

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

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

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

將其代入式（1）～式（3），即可得到TTI介質中群速度值，其中相位的傳播方向為

射線的傳播方向是

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

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

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

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

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

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

圖2 含起伏界面的三維層狀模型中反射、透射、轉換波計算示意圖

圖3 四面體內速度插值示意圖

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

3.3 波前計算和射線追蹤原理

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

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

3.4 節(jié)點旅行時計算方法

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

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

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

4 實例模擬和結果分析

4.1 群速度與相慢度角和射線角的關系

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

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

4.2 TTI介質群速度分布

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

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

4.3 計算精度分析

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

4.4 復雜TTI介質中多次波追蹤

本文所提算法可用于復雜TTI介質中多次波追蹤，為使展示結果清晰明了，數(shù)值模擬時采用層狀模型。選擇如圖7所示的三維層狀起伏模型：模型尺寸為80km×80km×50km；地表起伏，最大起伏幅度為2km；在25km深度處有一個凹陷反射界面，最大凹陷幅度為4km；底部傾斜反射界面位于47～50km深度處。炮點位于（40km，40km，10km）的位置（圖7中黑色圓點），31個檢波器以等間距（8.0km）分布在模型地表的三條邊上（圖7中灰色圓點）。第一層模型參數(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，即自炮點激發(fā)的P波經(jīng)第一、二層起伏界面的反射波。圖7b給出了多次透射、反射及轉換波的射線路徑：第一個震相為P1SV1P1，炮點激發(fā)的P波上行至地表，經(jīng)反射并轉換成SV波后下行至第一層反射界面，再次反射并轉換成P波上行至檢波器；第二個震相為P1SV1P2SV2P1，炮點激發(fā)的P波上行至地表，經(jīng)反射轉換為SV波后下行至第一層反射界面，透射且轉換成P波后繼續(xù)下行至第二層反射界面，然后反射轉換成SV波上行至第一層反射界面，最后透射轉換成P波上行至檢波器。

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

圖7 一次反射波P1 P1（藍線）和P1 P2 P2 P1（紅線）射線路徑（a）、多次反射（轉換）震相P1 SV1 P1（藍線）和P1 SV1 P2 SV2 P1（紅線）射線路徑（b）

5 結論

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

［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］ 馬德堂，朱光明.橫向各向同性介質中初至波的旅行時計算.石油地球物理勘探，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］ 張中杰，何樵登.各向異性介質中折射波及其走時曲線研究.西安石油學院學報，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］ 趙愛華，張美根，丁志峰.橫向各向同性介質中地震走時模擬.地球物理學報，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介質中初至波旅行時的搜索算法.石油地球物理勘探，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介質中多次波射線追蹤.石油地球物理勘探，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)格算法研究.中國地質大學學報（地球科學），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］ 孔娟華，鄭江濱.一種三維離散點數(shù)據(jù)生成結構四面體算法.計算機工程與科學，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介質中多次波射線追蹤.石油地球物理勘探，2017，52（1）：48-55.

1000-7210（2017）01-0048-08

*陜西省西安市長安大學地質工程與測繪學院地球物理系，710054。Email：baicy＠chd.edu.cn

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

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

（本文編輯：朱漢東）

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