居興國,郭 愷,劉定進

(中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

?

基于相速度的TTI介質(zhì)射線追蹤方法研究

居興國,郭 愷,劉定進

(中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

TTI介質(zhì)射線追蹤決定了TTI介質(zhì)建模、成像的效率和精度。在程函方程法射線追蹤的基礎(chǔ)上,提出了基于相速度的TTI介質(zhì)射線追蹤新策略,該策略以Hamilton公式為基礎(chǔ)推導(dǎo)出相速度法射線追蹤系統(tǒng),再代入群速度和出射角將其轉(zhuǎn)換為TTI介質(zhì)射線追蹤微分方程,利用龍哥庫塔解微分方程得到每一步延拓的坐標和出射角。相速度法射線追蹤在保證精度和穩(wěn)定性的前提下,提高了計算效率,所得出的參數(shù)可用于后期層析迭代反演等。模型試算結(jié)果表明基于相速度的TTI介質(zhì)射線追蹤方法有效且實用。

TTI介質(zhì);TTI建模;射線追蹤;程函方程;相速度

地下介質(zhì)普遍存在各向異性特征。在地層接近水平或各向異性程度較弱的地區(qū),各向同性或VTI介質(zhì)假設(shè)條件能夠得到較好的地下介質(zhì)成像效果。但是,如果各向異性程度較強,并且地層傾角較大,基于各向同性和VTI介質(zhì)的假設(shè)則無法得到地下介質(zhì)的高精度成像結(jié)果。此種情況下通常會出現(xiàn)成像道集無法拉平、構(gòu)造模糊或失真、成像位置和深度存在誤差等現(xiàn)象,給后續(xù)油藏描述和儲層預(yù)測帶來困難,從而降低鉆井成功率[1-4]。因此,隨著地震勘探的不斷發(fā)展以及勘探目標越來越復(fù)雜,各向同性和VTI介質(zhì)的假設(shè)條件已經(jīng)不能夠滿足勘探精度的需求,需要研究更加符合地下地質(zhì)特征的TTI介質(zhì)假設(shè)條件下的建模成像技術(shù)以提高復(fù)雜構(gòu)造的勘探精度,而TTI射線追蹤是TTI建模成像的關(guān)鍵技術(shù),它決定了建模成像的精度和效率,是復(fù)雜構(gòu)造精確成像研究的重要內(nèi)容之一。

目前各向異性介質(zhì)射線追蹤方法主要有[5]:打靶法、波前構(gòu)建法、插值法、最短路徑法和程函方程法等。打靶法基于漸進射線理論[6],適用于一般各向異性介質(zhì),但對于TTI這種復(fù)雜介質(zhì)則會出現(xiàn)陰影區(qū)問題。波前構(gòu)建法是根據(jù)當前的波前利用射線理論估計新的波前,白海軍等[7]雖然利用群速度和相速度表達的射線追蹤方程避免了復(fù)雜的特征值求解問題,但是實現(xiàn)方案較為復(fù)雜。插值法基于費馬原理[8-10],分為向前確定走時和向后確定射線路徑兩個過程,在各向異性介質(zhì)中往往得不到精確的出射點位置。最短路徑法基于惠更斯原理和費馬原理[11-13],能同時求取震源到模型所有節(jié)點的走時和相應(yīng)的射線路徑,但只適用于各向異性介質(zhì)初至波射線追蹤。VIDALE[14]基于擴張波前的思想建立了經(jīng)典的程函方程射線追蹤系統(tǒng),之后,ALKHALIFAH[15]提出了各向異性程函方程,并在此基礎(chǔ)上建立了各向異性射線追蹤系統(tǒng),該方法穩(wěn)定性強、精度高,是VTI介質(zhì)射線追蹤最主要的方法,但是在TTI介質(zhì)中,程函方程則非常繁瑣,無論是對于推導(dǎo)射線追蹤方程還是對于射線追蹤效率,影響都很大,因此該方法不適用于TTI介質(zhì)。

本文在前人研究的基礎(chǔ)上,利用Hamilton相速度公式推導(dǎo)出相速度射線追蹤方程,再將TTI介質(zhì)群速度和相速度公式代入相速度射線追蹤方程得到TTI介質(zhì)射線追蹤方程,然后利用射線參數(shù)與相角的關(guān)系將射線追蹤方程中射線參數(shù)與時間的偏微分公式轉(zhuǎn)換為射線出射角與時間的偏微分公式,最后采用龍哥庫塔法求解射線追蹤系統(tǒng)的微分方程組,從而得到射線追蹤每一步延拓的坐標、出射角和時間。這種方法相比較于程函方程法,公式簡單、物理意義明確,運算效率高,并且解微分方程得到的是出射角而不是程函方程方法中的射線參數(shù),可以提供后期層析反演所需要的數(shù)據(jù),更有利于運算。

1 方法原理

1.1 TTI介質(zhì)相速度和群速度

圖1為TTI介質(zhì)相角與群角示意圖,其中,TTI介質(zhì)中出射角(射線方向與z軸的夾角)為θ,對稱軸傾角(對稱軸與z軸的夾角)為θc,在本文公式中記為θ′,射線方位角(射線方向水平投影與0方位的夾角)為φ,對稱軸方位角(對稱軸水平投影與0方位的夾角)為φc,在本文公式中記為φ′,則TTI介質(zhì)中相角γ可表示為:

(1)

圖1 TTI介質(zhì)相角與群角示意

通過相角可以得到TTI介質(zhì)相速度公式:

(2)

其中:

式中:v為地震波相速度;vP0是沿對稱軸方向的縱波速度;ε和δ是Thomsen各向異性參數(shù)。

BERRYMAN提出射線速度可以表示為波動矢量和相速度的函數(shù),并給出了射線速度公式[16]:

(3)

式中:vG表示地震波射線速度,即群速度;v是地震波相速度;波動矢量k=kxI+kyJ+kzK;I,J,K代表直角坐標系三個方向的單位向量。

吳國忱等[17]以公式(3)為基礎(chǔ),推導(dǎo)出的TTI介質(zhì)三維群速度公式為:

(4)

1.2 基于相速度的射線追蹤

基于相速度的射線追蹤方法是從Hamilton相速度公式出發(fā)推導(dǎo)運動學(xué)射線追蹤系統(tǒng),然后將TTI介質(zhì)相速度和群速度代入該系統(tǒng),得到TTI介質(zhì)射線追蹤方程,此時的方程是坐標、射線參數(shù)和旅行時的微分方程,解此方程就會得到當前點的坐標、射線參數(shù)及傳播時間,可以用作計算下一個傳播點的初始值。但是射線參數(shù)只是中間變量,與速度和時間并無聯(lián)系,無法直接作為初始值使用,因此利用射線參數(shù)與角度的關(guān)系進一步將TTI介質(zhì)射線追蹤方程變?yōu)樽鴺?、出射角和旅行時的微分方程,解方程得到的坐標、出射角及傳播時間可以直接作為初始值進行計算。通過該微分方程就可以從出射點一直追蹤到檢波點,計算出一條射線曲線及旅行時信息?；谙嗨俣鹊纳渚€追蹤方法具體推導(dǎo)過程如下。

首先從Hamilton相速度公式出發(fā)[18]:

(5)

得到相應(yīng)的運動學(xué)射線追蹤方程:

(6)

式中:xi是坐標;T是旅行時;τ是時間步長;K是Hamilton相速度;p是射線參數(shù);v是TTI介質(zhì)相速度;vG是TTI介質(zhì)群速度。

將TTI介質(zhì)相速度((2)式)和群速度((4)式)公式代入射線追蹤系統(tǒng)((6)式),可以得到基于相速度的三維TTI介質(zhì)射線追蹤方程:

(7)

將射線參數(shù)與出射角的關(guān)系式px=sinθcosφ/v,py=sinθsinφ/v,pz=cosθ/v代入公式(7),即可得到三維TTI介質(zhì)射線追蹤方程出射角形式:

(8)

其中:

[cosθcosθ′cos(φ-φ′)-sinθsinθ′]+

2sinθcosθsin2(φ-φ′)

sinθcosθ′sin(φ-φ′)+

2sin2θsin(φ-φ′)cos(φ-φ′)

公式(8)即為基于相速度的TTI介質(zhì)三維射線追蹤方程,它是以等時間為步長進行延拓的,即每一點旅行時一樣。先計算每一步延拓點的坐標,然后采用相速度對坐標的偏導(dǎo)數(shù)計算出射角?；谙嗨俣鹊腡TI介質(zhì)三維射線追蹤方程形式簡單,方程中的變量物理意義明確,并且只計算5個參數(shù),運算穩(wěn)定、計算效率高。

最后用龍哥庫塔法求解微分方程組(公式(8)),可以得到射線追蹤每一步延拓的坐標、出射角和時間。龍哥庫塔公式為:

(9)

式中:Yn和Yn+1分別是n次和n+1次延拓得到的值;h是延拓步長;K為龍哥庫塔四階計算值;tn是n次延拓的時間。

2 模型試算

2.1 正演模擬驗證

在驗證方法的精度和效率之前,首先驗證方法的正確性。將TTI聲波正演模擬產(chǎn)生的單炮記錄與TTI射線追蹤得到的旅行時曲線進行疊合,如果完全重合,說明本文方法所得到的旅行時信息正確;如果不完全重合,說明本文方法存在問題,需要修改。

設(shè)計一個雙層模型,模型大小為橫向4000m,縱向4000m,網(wǎng)格大小為橫向10m,縱向10m,網(wǎng)格數(shù)量為橫向401個,縱向401個。模型第一層速度vP0=3000m/s,各向異性參數(shù)ε=0.1,δ=0.1,對稱軸傾角θ=45°,第二層速度vP0=3250m/s,各向異性參數(shù)ε=0.12,δ=0.12,對稱軸傾角θ=25°,反射界面深度為2000m,炮點位于地表(2000m,10m)處,兩組檢波器分別放置在深度為1990m和3990m兩個層位上，檢波器間距均為10m。將TTI聲波正演所得單炮記錄(黑白同相軸)與相速度法得到的旅行時曲線(綠色曲線)進行疊合,如圖2所示,其中圖2a是1990m處的疊合圖,圖2b是3990m處的疊合圖。

圖2 旅行時與正演記錄疊合顯示a 1990m; b 3990m

由圖2a可以看出,在第一層均勻模型中,利用本文提出的基于相速度的射線追蹤法得到的旅行時曲線與聲波正演模擬產(chǎn)生的單炮記錄吻合程度非常高,說明本文方法在均勻模型中的應(yīng)用是正確的。由圖2b 可以看出,射線從地表炮點出發(fā),在兩種介質(zhì)分界面發(fā)生了一次透射,最終到達3990m處的檢波點,利用基于相速度的射線追蹤法得到的旅行時曲線與聲波正演模擬產(chǎn)生的單炮記錄吻合程度同樣非常高,說明該方法在復(fù)雜介質(zhì)中依然準確。

2.2 二維模型試算

通過二維模型驗證方法的精度和效率。程函方程法的精度很高,但是運算速度較慢,這里分別采用恒速度模型和變速度模型繪制射線路徑和等時線(波前面)圖,并通過等時線與程函方程法進行精度對比,如果兩者等時線重合率高,說明基于相速度的射線追蹤法的精度同樣很高;如果兩者等時線出現(xiàn)了較大的偏差,則需要進一步驗證哪種方法精度更高。再采用一個較為復(fù)雜的二維模型利用密集射線追蹤驗證基于相速度的射線追蹤法的計算效率,如果相速度法的效率明顯高于程函方程法,說明相速度法更具有實用性;如果兩種方法的計算效率相當,則說明相速度法與程函方程法具有相同的實用性。

設(shè)計一個二維模型,模型大小為橫向4000m,縱向4000m,網(wǎng)格大小為橫向10m,縱向10m,網(wǎng)格數(shù)量為橫向401個,縱向401個。圖3中的模型速度vP0=2000m/s,各向異性參數(shù)ε=0.3,δ=-0.1,對稱軸傾角θ=45°,炮點在模型正中間,坐標為(2000m,2000m)。圖4中的模型速度vP0=(2000+0.5z)m/s,各向異性參數(shù)ε=0.3,δ=-0.1,對稱軸傾角θ=45°,炮點在地表中間位置,坐標為(2000m,0)。

圖3和圖4分別為恒速度場、變速度場射線路徑與等時線以及基于相速度的射線追蹤法與程函方程法所得等時線疊合顯示圖。圖3a和圖4a為基于相速度的射線追蹤法得到的射線路徑和等時線,背景為速度場。從等時線可以看出,各向異性情況下波前面不再是圓形,而是橢圓形,并且隨著對稱軸的傾斜發(fā)生旋轉(zhuǎn),變速情況下射線發(fā)生了彎曲,出現(xiàn)了回轉(zhuǎn)現(xiàn)象。圖3b和圖4b為相速度法和程函方程法得到的等時線疊合圖,其中黑色為程函方程的等時線,紅色為相速度的等時線。因為兩種方法的等時線完全重合,所得到的旅行時信息一致,所以圖中只能看到紅色的相速度等時線,說明相速度法射線追蹤的效果和精度跟程函方程法一致,驗證了本文方法的有效性。

圖3 恒速度場(速度2000m/s)射線路徑與等時線(a)以及相速度法與程函方程法所得等時線疊合顯示(b)

圖4 變速度場射線路徑與等時線(a)以及相速度法與程函方程法所得等時線疊合顯示(b)

采用復(fù)雜的二維TTI模型進一步驗證本文方法的有效性和運算效率。圖5a為速度場,圖5b為ε場,圖5c為δ場,圖5d為傾角場,圖5e為射線路徑圖,圖5f為等時線圖。該模型包含高速鹽丘和高陡構(gòu)造等復(fù)雜介質(zhì),射線在淺層出現(xiàn)了回轉(zhuǎn)現(xiàn)象,在鹽丘頂部出現(xiàn)了透射和全反射,等時線也體現(xiàn)了射線的傳播方向和傳播速度。

表1為該TTI模型采用兩種方法的計算時間,射線密度X°/180°表示射線出射角范圍是0～180°,方向向下,每隔X°出射一條射線,旅行時單位為s。我們采用3組射線密度(10°/180°,2°/180°和0.5°/180°)分別測試程函方程法和相速度法射線追蹤的用時。從表1可以看出,隨著射線密度的增大,兩種方法的用時都是成倍數(shù)增長,射線密度越大,兩者用時差距也越大,當射線密度為0.5°/180°時,相速度射線追蹤法用時比程函方程法少了12.62s,可見,相速度法的計算效率明顯優(yōu)于程函方程法,這種優(yōu)勢在三維模型等大規(guī)模數(shù)據(jù)運算時更為明顯。

通過模型試算發(fā)現(xiàn),雖然相速度法射線追蹤與程函方程法射線追蹤具有相同的精度和效果,都能夠為建模和成像提供準確的射線路徑和旅行時信息,但是相速度法的計算效率遠超程函方程法,能夠在很大程度上縮短TTI建模和成像的時間,比程函方程法更適用于大規(guī)模實際數(shù)據(jù)的處理。

圖5 復(fù)雜二維TTI介質(zhì)的參數(shù)場、射線路徑和等時線a 速度場; b ε場; c δ場; d 傾角場; e 射線路徑; f 等時線

射線密度10．0°/180°2．0°/180°0．5°/180°程函方程法用時/s1．459．2736．01相速度法用時/s0．925．7623．39

2.3 三維模型試算

通過三維模型試算驗證方法在三維空間中的射線追蹤路徑和旅行時信息,以及等時線在空間中的旋轉(zhuǎn)情況。設(shè)計一個三維模型,模型大小為橫1000m,寬1000m,深度1000m,x,y,z方向網(wǎng)格大小均為10m,x,y,z方向網(wǎng)格數(shù)量均為101個,炮點在模型正中間,坐標為(500m,500m,500m)。

圖6中,模型速度vP0=2000m/s,各向異性參數(shù)ε=0.3,δ=-0.1,對稱軸傾角θ=45°,方位角φ=45°;圖7中,模型速度vP0=(2000+0.5z)m/s,各向異性參數(shù)ε=0.3,δ=-0.1,對稱軸傾角θ=45°,方位角φ=45°。圖6a和圖7a分別為模型的速度場;圖6b 和圖7b為射線路徑圖;圖6c和圖7c為等時線圖。由圖6b和圖7b可以看出,恒速度場中的射線呈直線傳播,變速度場中的射線呈彎曲狀,出現(xiàn)了回轉(zhuǎn)現(xiàn)象;由圖6c和圖7c可以看出,波前面不但垂向上旋轉(zhuǎn)了45°(傾角),水平方向也旋轉(zhuǎn)了45°(方位角)。三維模型測試結(jié)果說明相速度射線追蹤方法在三維空間中也能取得正確有效的射線路徑和旅行時信息。

圖6 三維模型恒速度試算結(jié)果a 速度場; b 射線路徑; c 等時線

圖7 三維模型變速度試算結(jié)果a 速度場; b 射線路徑; c 等時線

3 認識與結(jié)論

本文在程函方程法TTI介質(zhì)射線追蹤的基礎(chǔ)上,提出了一種簡潔、高效的基于相速度的TTI介質(zhì)射線追蹤方法。該方法以Hamilton相速度為基礎(chǔ)推導(dǎo)出了射線追蹤系統(tǒng),并利用群速度公式和出射角公式將其轉(zhuǎn)換為TTI介質(zhì)出射角格式,得到基于相速度的TTI介質(zhì)射線追蹤方程,其公式簡單、物理意義明確,并且所得出參數(shù)能夠用于層析反演等后續(xù)處理。二維、三維模型試算結(jié)果表明,本文提出的相速度法射線追蹤的精度較高、效果較好,能提供準確的射線路徑和旅行時信息,計算效率較高,非常適用于大規(guī)模實際數(shù)據(jù)處理。

在實際生產(chǎn)中,射線追蹤可用于射線類偏移和建模,本文所述射線追蹤方法只反映透射波、回轉(zhuǎn)波和折射波,不能夠反映反射波,因此,在實際使用過程中需要確定反射點位置和反射界面的傾角和方位角,再遵循反射定律向炮點和檢波點進行射線追蹤,從而完成一個炮檢對的射線路徑模擬。由于無法模擬反射波,因此也無需考慮邊界效應(yīng),但是,必須考慮全反射現(xiàn)象,這種現(xiàn)象是我們不想看到的,為了避免發(fā)生這種現(xiàn)象,在射線追蹤之前一般需要將模型進行平滑,消除界面的突變帶來的不穩(wěn)定性,再進行射線追蹤計算。

[1] 張千祥,王德利,周進舉.二維TTI介質(zhì)的純P波波動方程數(shù)值模擬[J].石油物探,2015,54(5):485-492 ZHANG Q X,WANG D L,ZHOU J J.Acoustic wave equation numerical simulation for pure P-wave in 2D TTI medium[J].Geophysical Prospecting for Petroleum,2015,54(5):485-492

[2] 郭立鵬,楊勤勇,李振春,等.復(fù)雜各向異性介質(zhì)初至波射線追蹤[J].石油物探,2016,55(1):18-24 GUO L P,YANG Q Y,LI Z C,et al.First arrival ray tracing in complex anisotropic medium[J].Geophysical Prospecting for Petroleum,2016,55(1):18-24

[3] 黃光南,鄧居智,李紅星,等.非均勻節(jié)點網(wǎng)格TI介質(zhì)反射波射線追蹤研究[J].石油物探,2016,55(1):25-32 HUANG G N,DENG J Z,LI H X,et al.Reflected wave ray tracing in TI medium based on the nonuniform node meshes[J].Geophysical Prospecting for Petroleum,2016,55(1):25-32

[4] 郭愷,王鵬燕,婁婷婷.基于P-SV波分離的VTI介質(zhì)射線追蹤方法[J].地球科學(xué)進展,2015,30(9):1028-1033GUO K,WANG P Y,LOU T T.A new ray tracing method for VTI medium based on separated P-SV waves[J].Advances in Earth Science,2015,30(9):1028-1033

[5] 趙愛華,張美根,丁志峰.橫向各向同性介質(zhì)中地震波走時模擬[J].地球物理學(xué)報,2006,49(6):1762-1769 ZHAO A H,ZHANG M G,DING Z F.Seismic traveltime computation for transversely isotropic media[J].Chinese Journal of Geophysics,2006,49(6):1762-1769

[6] CERVENY V.Seismic ray theory[M].Cambridge:Cambridge University Press,2001:148-153

[7] 白海軍,孫贊東,王學(xué)軍.基于波前構(gòu)建法的TTI介質(zhì)射線追蹤[J].石油地球物理勘探,2011,46(增刊1):1-6 BAI H J,SUN Z D,WANG X J.Ray tracing in TTI media based on wavefront construction method[J].Oil Geophysical Prospecting,2011,46(S1):1-6

[8] 鄧懷群,劉雯林.橫向各向同性介質(zhì)中地震波旅行時的計算[J].石油地球物理勘探,2000,35(4):508-516DENG H Q,LIU W L.Computation of travel timein transversely isotropic media[J].Oil Geophysical Prospecting,2000,35(4):508-516

[9] 馬德堂,朱光明.橫向各向同性介質(zhì)中的初至波旅行時計算[J].石油地球物理勘探,2006,41(1):26-31 MA D T,ZHU G M.Computation of seismic firstbreak travel time in transversely isotropic media[J].Oil Geophysical Prospecting,2006,41(1):26-31

[10] 馬德堂,朱光明,范廷恩.二維TTI介質(zhì)中初至波旅行時的搜索算法[J].石油地球物理勘探,2011,46(5):710-714 MA D T,ZHU G M,FAN Y E.2D TTI medium first arrival travel time search algorithm[J].Oil Geophysical Prospecting,2011,46(5):710-714

[11] MOSER T J.Shortest path calculation of seismic rays[J].Geophysics,1991,56(1):59-67

[12] ZHOU B,GREENHALGH S A.Shortest path ray tracing for most general 2D/3D anisotropic media[J].Journal of Geophysics and Engineering,2005,2(1):54-63

[13] BAI C,GREENHALGH S A,ZHOU B.3D ray tracing using a modified shortest path method[J].Geophysics,2007,72(4):T27-T36

[14] VIDALE J.Finite-difference calculation of travel times[J].Bulletin of the Seismological Society of America,1988,78(6):2062-2076

[15] ALKHALIFAH T.An acoustic wave equation for anisotropic media[J].Geophysics,2000,65(4):1239-1250

[16] BERRYMAN J G.Long-wave elastic anisotropy in transversely isotropic media[J].Geophysics,1979,44(5):896-917

[17] 吳國忱,梁鍇,戚艷平.三維TTI介質(zhì)相速度和群速度[J].地球物理學(xué)進展,2009,24(6):2097-2105 WU G C,LIANG K,QI Y P.Phase velocity and group velocity in 3D TTI media[J].Progress in Geophysics,2009,24(6):2097-2105

[18] 段鵬飛,程玖兵,陳三平,等.TI介質(zhì)局部角度域射線追蹤與疊前深度偏移成像[J].地球物理學(xué)報,2013,56(1):269-279 DUAN P F,CHENG J B,CHEN S P,et al.Local angle-domain ray tracing and prestack depth migration in TI medium[J].Chinese Journal of Geophysics,2013,56(1):269-279

(編輯:朱文杰)

Research on a ray tracing method for TTI medium based on phase velocity

JU Xingguo,GUO Kai,LIU Dingjin

(SinopecGeophysicalResearchInstitute,Nanjing211103,China)

The ray tracing for TTI medium determines the efficiency and precision of modeling and imaging in TTI medium.In this paper we proposed a new strategy named phase velocity ray tracing on the basis of eikonal equation ray tracing.This strategy established phase velocity ray tracing system based on Hamilton formula.Through substituting group velocity and emergence angle we convert the system to ray tracing differential equation for TTI medium.Then we use Runge-Kutta method to solve the differential equation and get coordinates and emergence angles for every continuation.Phase velocity ray tracing method improved the operation efficiency with high accuracy and stability,which is benefit to the later tomography iterative inversion.The results of numerical tests prove the correctness and practicality of this method.

TTI medium,TTI modeling,ray tracing,Eikonal equation,phase velocity

2016-06-16;改回日期:2016-08-24。

居興國(1965—),男,高級工程師,現(xiàn)在主要從事地震資料處理技術(shù)研究工作。

國家高技術(shù)研究發(fā)展計劃(863計劃)課題(2011AA060303)資助。

P631

A

1000-1441(2017)02-0171-08

10.3969/j.issn.1000-1441.2017.02.002

This research is financially supported by National High-tech R&D Program (863 Program) of China (Grant No.2011AA060303).

居興國,郭愷,劉定進.基于相速度的TTI介質(zhì)射線追蹤方法研究[J].石油物探,2017,56(2):-178

JU Xingguo,GUO Kai,LIU Dingjin.Research on a ray tracing method for TTI medium based on phase velocity[J].Geophysical Prospecting for Petroleum,2017,56(2):-178