黃光南,鄧居智,李紅星,李澤林,張　華,王安東

(1.東華理工大學(xué)核技術(shù)應(yīng)用教育部工程研究中心,江西南昌 330013;2.中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430074;3.中國石油大學(xué)(北京)油氣資源與探測(cè)國家重點(diǎn)實(shí)驗(yàn)室,北京 102249)

?

非均勻節(jié)點(diǎn)網(wǎng)格TI介質(zhì)反射波射線追蹤研究

黃光南1,2,3,鄧居智1,李紅星1,李澤林1,張華1,王安東1

(1.東華理工大學(xué)核技術(shù)應(yīng)用教育部工程研究中心,江西南昌 330013;2.中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430074;3.中國石油大學(xué)(北京)油氣資源與探測(cè)國家重點(diǎn)實(shí)驗(yàn)室,北京 102249)

摘要:地下介質(zhì)不僅具有各向異性性質(zhì),還存在一定程度的復(fù)雜地質(zhì)構(gòu)造。運(yùn)用均勻節(jié)點(diǎn)網(wǎng)格射線追蹤方法計(jì)算旅行時(shí)要求網(wǎng)格單元的劃分非常小以達(dá)到較高的旅行時(shí)精度,這種做法會(huì)產(chǎn)生大量的網(wǎng)格單元從而降低計(jì)算效率?；诰鶆蚬?jié)點(diǎn)網(wǎng)格算法,研究了非均勻節(jié)點(diǎn)網(wǎng)格TI介質(zhì)反射波射線追蹤方法,采用較大的網(wǎng)格單元,通過在網(wǎng)格單元的每條邊增加次一級(jí)節(jié)點(diǎn)來提高旅行時(shí)計(jì)算的精度與效率。首先討論了qP,qSV和qSH波的群速度計(jì)算方法,然后給出了非均勻節(jié)點(diǎn)網(wǎng)格TI介質(zhì)反射波射線追蹤算法,最后分別求取了3種不同對(duì)稱軸傾角的TI介質(zhì)起伏地層模型qP,qSV和qSH反射波的射線路徑和旅行時(shí)。結(jié)果表明:①非均勻節(jié)點(diǎn)網(wǎng)格射線追蹤算法能夠適用于復(fù)雜TI介質(zhì)模型;②相同模型中不同波模式具有不同的反射路徑和旅行時(shí);③相同波模式在不同對(duì)稱軸傾角模型中也具有不同的反射路徑和旅行時(shí)。

關(guān)鍵詞：起伏地層模型;非均勻節(jié)點(diǎn)網(wǎng)格;各向異性TI介質(zhì);射線追蹤算法

各向異性介質(zhì)射線追蹤可以模擬地震波在復(fù)雜地下介質(zhì)中的傳播，為實(shí)際地震資料解釋提供有效的技術(shù)手段，因而被廣泛應(yīng)用于各向異性地震記錄合成、各向異性參數(shù)反演、各向異性地震偏移成像等技術(shù)領(lǐng)域。國外學(xué)者對(duì)各向異性介質(zhì)射線追蹤技術(shù)的研究較早[1-2]。如CERVENY[3]推導(dǎo)了不均勻各向異性介質(zhì)射線走時(shí)和振幅的計(jì)算方法,為實(shí)施射線追蹤算法奠定了理論基礎(chǔ);CERVENY等[4]針對(duì)非均勻、弱各向異性介質(zhì)情形,提出了無需射線追蹤的線性化走時(shí)計(jì)算方法;GAJEWSKI等[5]將打靶法射線追蹤技術(shù)成功擴(kuò)展至地層起伏變化的三維各向異性介質(zhì)模型,通過數(shù)值模擬計(jì)算得到了射線路徑、旅行時(shí)和射線振幅信息;SHEARER等[6]利用線性梯度算法刻畫各向異性介質(zhì)模型,通過求解多項(xiàng)式方程得到射線路徑與旅行時(shí),利用數(shù)值模擬結(jié)果展示了相關(guān)理論和算法的有效性。

國內(nèi)最早計(jì)算各向異性旅行時(shí)的方法基于時(shí)距曲線方程。例如:張文生等[7]推導(dǎo)了水平層狀各向異性介質(zhì)走時(shí)近似公式,通過數(shù)值模擬P波和SV波時(shí)距曲線,驗(yàn)證了近似表達(dá)式的有效性;苑書金等[8]利用泰勒級(jí)數(shù)法推導(dǎo)了P波和SV波在短、中、長(zhǎng)排列的反射波旅行時(shí)公式,通過數(shù)值模擬驗(yàn)證了這種計(jì)算方法的有效性和可行性。后來,國內(nèi)學(xué)者逐漸提出了算法相對(duì)簡(jiǎn)單的各向異性射線追蹤方法。例如:鄧懷群等[9]根據(jù)群速度和相速度的關(guān)系,結(jié)合各向異性射線參數(shù)表達(dá)式,推導(dǎo)出了射線追蹤的射線參數(shù)表達(dá)式;孔選林等[10]利用逐段迭代射線追蹤方法計(jì)算簡(jiǎn)單層狀各向異性模型的地震波走時(shí),有效地模擬了P波和SV波的時(shí)距曲線;熊金良等[11]利用水平層狀介質(zhì)時(shí)距曲線方程和二分法射線追蹤算法計(jì)算層狀各向異性模型的方位旅行時(shí),當(dāng)?shù)叵麓嬖诹严稌r(shí),不同方位的時(shí)距曲線之間存在較大差別;李建國等[12]將試射法射線追蹤算法應(yīng)用于VTI介質(zhì)理論模型研究,利用VSP觀測(cè)系統(tǒng)計(jì)算了上、下行P波的射線路徑和旅行時(shí);郝奇等[13]基于二維層狀VTI介質(zhì)模型推導(dǎo)了qP波和qSV波的精確Snell定律公式,歸納了VTI介質(zhì)qP波和qSV波的射線追蹤算法流程,通過水平和起伏界面模型的反射波射線追蹤數(shù)值模擬,證明了該算法的準(zhǔn)確性和高效性。

上述地球物理學(xué)者研究各向異性射線追蹤時(shí),大多采用地層對(duì)稱軸傾角垂直的VTI介質(zhì),因?yàn)樗亲顬楹?jiǎn)單的一種各向異性介質(zhì)情形,利用它推導(dǎo)的射線追蹤算法可以避免相對(duì)復(fù)雜的計(jì)算。如QIAN等[14]通過建立群速度和相速度之間的顯式關(guān)系式尋找有限差分走時(shí)場(chǎng)的外推方向,方便地求出各向異性介質(zhì)qP波的初至走時(shí)信息,通過二維和三維VTI介質(zhì)數(shù)值模擬實(shí)驗(yàn),證明了該方法具有較高的數(shù)值精度和運(yùn)算效率。劉玉柱等[15]實(shí)現(xiàn)了VTI介質(zhì)射線追蹤算法,采用兩步法反演兩個(gè)Thomsen參數(shù),針對(duì)異常體模型的反演取得了較好的數(shù)值效果。SCHNEIDER[16]利用擾動(dòng)理論來求解P波的程函方程,當(dāng)各向異性強(qiáng)度在一定范圍內(nèi)時(shí),計(jì)算得到的各向異性旅行時(shí)具有較高的精度。然而,自然界地層形成之后往往會(huì)伴隨一系列的地殼運(yùn)動(dòng),導(dǎo)致地層對(duì)稱軸傾角發(fā)生變化;因此,野外實(shí)際地層大多數(shù)屬于TTI介質(zhì)情形。ZHOU等[17]將最短路徑射線追蹤算法拓展至各向異性TTI介質(zhì)情形,計(jì)算了傾斜界面模型3種體波的射線路徑和旅行時(shí)。WANG[18]針對(duì)各向異性介質(zhì)彎曲射線法產(chǎn)生的高度非線性問題修改了牛頓迭代法,用于求解穩(wěn)定非線性問題,并利用井間觀測(cè)方式反演了各向異性介質(zhì)模型參數(shù)。趙后越等[19]結(jié)合均勻節(jié)點(diǎn)網(wǎng)格和非均勻節(jié)點(diǎn)網(wǎng)格算法實(shí)現(xiàn)了起伏地表各向異性TTI介質(zhì)的射線追蹤,理論模型試算結(jié)果表明該方法具有較高的數(shù)值精度。QIN等[20]根據(jù)惠更斯原理推導(dǎo)了各向異性介質(zhì)初至波走時(shí)的計(jì)算方法,并將該方法應(yīng)用于幾種各向異性TTI介質(zhì)模型,通過與波動(dòng)方程所得旅行時(shí)相比,說明它是一種穩(wěn)定和精確的各向異性射線追蹤算法。WANG等[21]針對(duì)復(fù)雜各向異性TTI介質(zhì),利用基于網(wǎng)格的射線追蹤算法,反演得到了模型的各向異性參數(shù)剖面。KUMAR等[22]針對(duì)各向異性TTI介質(zhì),提出了一種直接計(jì)算P波初至走時(shí)的方法,并將它與克?；舴蚍e分偏移算法相結(jié)合,極大地改善了偏移成像的精度。PRATT等[23]運(yùn)用彎曲射線追蹤算法反演井間各向異性TTI介質(zhì)的參數(shù)分布和地層的對(duì)稱軸傾角。LOU[24]將求解程函方程的快速匹配方法用于TTI介質(zhì)的初至旅行時(shí)求取,數(shù)值實(shí)例表明它是一種穩(wěn)定和精確的初至旅行時(shí)計(jì)算方法。本文研究非均勻節(jié)點(diǎn)網(wǎng)格各向異性TI介質(zhì)反射波射線追蹤算法,它不僅適用于任意對(duì)稱軸傾角的地層模型,也適用于起伏地層模型,其運(yùn)用條件更加貼近于實(shí)際地質(zhì)情形。

1各向異性介質(zhì)群速度計(jì)算方法

各向同性介質(zhì)射線追蹤算法基于速度網(wǎng)格模型,而各向異性TI介質(zhì)射線追蹤算法基于群速度網(wǎng)格模型。ZHOU等[25]提出了兩種計(jì)算各向異性介質(zhì)群速度的方法:特征值法和特征向量法。對(duì)于二維各向異性模型,根據(jù)彈性模量參數(shù){c11,c13,c33,c44,c66},利用特征值法可以得到群速度在水平與垂直方向的表達(dá)式為:

(1)

(2)

其中,?vm/??是相速度表達(dá)式對(duì)?的偏導(dǎo)數(shù),計(jì)算公式如下:

(3)

其中,

(4)

式中:v1,v2和v3分別代表qP波、qSV波和qSH波速度;c11,c13,c33,c44,c66為彈性模量參數(shù)。DALEY等[27]給出了各向異性介質(zhì)的相速度表達(dá)式,即克里斯托弗矩陣的特征值:

(5)

P和Q的表達(dá)式為:

Q=Q1Q2-Q3

(6)

其中Q1,Q2和Q3的表達(dá)式為:

Q1=c44cos2?+c11sin2?

Q2=c33cos2?+c44sin2?

(7)

公式(2)至公式(7)中,角度?為相慢度向量n=(sinθ,cosθ)和各向異性介質(zhì)對(duì)稱軸傾角方向ez=(sinθ0,cosθ0)之間的夾角。

2非均勻節(jié)點(diǎn)網(wǎng)格TI介質(zhì)反射波射線追蹤算法

與各向同性介質(zhì)旅行時(shí)計(jì)算方法類似,各向異性介質(zhì)也可以利用線性積分沿著射線路徑計(jì)算旅行時(shí),計(jì)算表達(dá)式為:

(8)

式中:x為空間坐標(biāo)向量,r0是射線方向從xA到xB的單位向量,ds為沿射線路徑的線段。由(8)式可知旅行時(shí)計(jì)算需要確定群速度U(x,θ0,r0)和射線路徑R(x)。第1節(jié)給出了群速度計(jì)算方法,這里討論如何計(jì)算旅行時(shí)。根據(jù)費(fèi)馬原理,在空間位置xB處的旅行時(shí)計(jì)算表達(dá)式為:

(9)

其中,ΩB是模型內(nèi)xB點(diǎn)的鄰域。將(9)式與最短路徑射線追蹤算法結(jié)合,可以得到各向異性射線追蹤算法。下面介紹各向異性TI介質(zhì)非均勻節(jié)點(diǎn)網(wǎng)格反射波射線追蹤算法的實(shí)現(xiàn)步驟。

1) 劃分非均勻節(jié)點(diǎn)網(wǎng)格參數(shù)模型。固定網(wǎng)格單元大小,將二維各向異性模型m(x)劃分成由主節(jié)點(diǎn)組成的網(wǎng)格模型:

(10)

其中,Nx,Nz分別代表x和z軸方向的網(wǎng)格個(gè)數(shù)。這樣模型就可以用很多矩形網(wǎng)格組成,每個(gè)矩形網(wǎng)格的4個(gè)角被定義為主節(jié)點(diǎn)。根據(jù)這種網(wǎng)格劃分方法,通常需要將模型劃分為非常小的網(wǎng)格以獲得較高的旅行時(shí)計(jì)算精度。采用BAI等[28]提出的非均勻節(jié)點(diǎn)網(wǎng)格算法將模型劃分為較大的網(wǎng)格單元,通過在網(wǎng)格單元的每條邊增加次一級(jí)節(jié)點(diǎn)來提高旅行時(shí)計(jì)算的精度與效率;將網(wǎng)格模型轉(zhuǎn)化成群速度模型Um(xk,θ0,r0),m=1,2,3分別對(duì)應(yīng)qP波、qSV波和qSH波。如果次一級(jí)的網(wǎng)格節(jié)點(diǎn)數(shù)為N2,那么網(wǎng)格模型的節(jié)點(diǎn)總數(shù)為N1+N2。另外,如果炮點(diǎn)、檢波點(diǎn)位置和網(wǎng)格模型的節(jié)點(diǎn)位置不重合,那么必須增加炮點(diǎn)和檢波點(diǎn)位置節(jié)點(diǎn)。因此整個(gè)模型包括3類網(wǎng)格節(jié)點(diǎn):主節(jié)點(diǎn)、次節(jié)點(diǎn)、炮點(diǎn)和檢波點(diǎn)位置對(duì)應(yīng)的節(jié)點(diǎn)。

2) 初始化走時(shí)場(chǎng)。對(duì)炮點(diǎn)所在節(jié)點(diǎn)的走時(shí)賦初值0,其它網(wǎng)格節(jié)點(diǎn)的走時(shí)賦初值無窮大(例如107),從炮點(diǎn)所在節(jié)點(diǎn)開始,計(jì)算炮點(diǎn)周圍節(jié)點(diǎn)的走時(shí),再從相鄰節(jié)點(diǎn)計(jì)算其它節(jié)點(diǎn)的走時(shí),并記錄相鄰節(jié)點(diǎn)的走時(shí)最小值。

3) 非均勻節(jié)點(diǎn)網(wǎng)格初至波走時(shí)計(jì)算。單個(gè)網(wǎng)格內(nèi)部?jī)蓚€(gè)節(jié)點(diǎn)之間的旅行時(shí)可以利用如下公式計(jì)算:

(11)

4) 非均勻節(jié)點(diǎn)網(wǎng)格反射波旅行時(shí)計(jì)算。如果各向異性模型存在復(fù)雜的地層界面,那么地震波會(huì)在界面處產(chǎn)生反射波。可以用一系列密集的網(wǎng)格節(jié)點(diǎn)來代表地層界面,地層界面節(jié)點(diǎn)與前述3類節(jié)點(diǎn)組合起來,可以得到一個(gè)新的節(jié)點(diǎn)集合,它包括主節(jié)點(diǎn)、次節(jié)點(diǎn)、反射界面節(jié)點(diǎn)和炮、檢點(diǎn)位置節(jié)點(diǎn)。每個(gè)地層界面對(duì)應(yīng)的反射波旅行時(shí)τ的計(jì)算表達(dá)式為:

(12)

其中,τS(xint)和τR(xint)分別是炮點(diǎn)S與檢波點(diǎn)R到界面節(jié)點(diǎn)位置xint的初至走時(shí),Ωint代表反射界面節(jié)點(diǎn)的集合。在反射波旅行時(shí)計(jì)算過程中,分別假設(shè)炮點(diǎn)S與檢波點(diǎn)R為震源點(diǎn)位置,進(jìn)行兩次正演走時(shí)計(jì)算,然后以(12)式為準(zhǔn)則尋找相應(yīng)的界面反射節(jié)點(diǎn),從而求取反射波的旅行時(shí)與射線路徑。值得注意的問題是:在計(jì)算反射波旅行時(shí)時(shí),群速度模型只與反射界面以上的非均勻節(jié)點(diǎn)網(wǎng)格模型有關(guān),也就是說,在計(jì)算τS(xint)和τR(xint)時(shí),只需要截取反射界面以上非均勻節(jié)點(diǎn)網(wǎng)格模型參與運(yùn)算,因此可以提高算法的整體運(yùn)算效率。

5) 非均勻節(jié)點(diǎn)網(wǎng)格反射波射線路徑求取。當(dāng)炮點(diǎn)和檢波點(diǎn)分別完成走時(shí)計(jì)算后,從非均勻節(jié)點(diǎn)網(wǎng)格模型中得到入射節(jié)點(diǎn)系列,根據(jù)入射節(jié)點(diǎn)(它包括主節(jié)點(diǎn)、次節(jié)點(diǎn)和界面節(jié)點(diǎn))的序號(hào)以反向追蹤的形式求取連接檢波點(diǎn)與炮點(diǎn)之間的射線路徑,從而得到共炮點(diǎn)道集的射線路徑。根據(jù)這種方法可以獲得所有炮點(diǎn)對(duì)應(yīng)的檢波點(diǎn)的射線路徑。

3數(shù)值模擬

利用非均勻節(jié)點(diǎn)網(wǎng)格射線追蹤算法對(duì)長(zhǎng)度800m,深度500m的起伏地層各向異性模型進(jìn)行了射線追蹤數(shù)值模擬。該模型含有兩個(gè)起伏界面和一個(gè)水平界面,第一層介質(zhì)的彈性模量參數(shù)為c11=9.08,c13=2.98,c33=7.53,c44=2.27,c66=3.84,第二層介質(zhì)的彈性模量參數(shù)為c11=20.3,c13=9.58,c33=22.3,c44=8.35,c66=11.35,第三層介質(zhì)的彈性模量參數(shù)為c11=13.86,c13=4.31,c33=10.93,c44=3.31,c66=4.34。炮點(diǎn)位于地表400m處,檢波點(diǎn)共33個(gè),第1個(gè)檢波點(diǎn)位于地表0處,第33個(gè)檢波點(diǎn)位于800m處,檢波點(diǎn)間距為25m(參見圖1)。

圖1給出了網(wǎng)格模型的主節(jié)點(diǎn)和次節(jié)點(diǎn)分布。為了方便圖形顯示和易于理解,這里假設(shè)模型的網(wǎng)格間距為50m,那么主節(jié)點(diǎn)數(shù)為187個(gè);網(wǎng)格的每條邊增加1個(gè)次一級(jí)節(jié)點(diǎn),那么次一級(jí)節(jié)點(diǎn)數(shù)為346個(gè),主次節(jié)點(diǎn)的總數(shù)為533個(gè)。由于網(wǎng)格間距較大,計(jì)算效率相對(duì)較高;同時(shí),網(wǎng)格模型的主次節(jié)點(diǎn)數(shù)較大,算法的旅行時(shí)計(jì)算精度也相對(duì)較高。圖2 為qP波在模型內(nèi)部的波前等值線圖,與均勻節(jié)點(diǎn)網(wǎng)格算法產(chǎn)生的波前等值線圖相比,由于它采用了大量的次一級(jí)網(wǎng)格節(jié)點(diǎn),其波前等值線不夠光滑,但是它具有較高的數(shù)值精度。為了提高旅行時(shí)的數(shù)值精度,模型的網(wǎng)格間距在水平與垂直方向上均為5m,在網(wǎng)格單元的每一條邊增加7個(gè)次一級(jí)網(wǎng)格節(jié)點(diǎn),這樣一個(gè)網(wǎng)格單元的主次節(jié)點(diǎn)數(shù)為32個(gè)。這里研究第一層和第二層介質(zhì)對(duì)稱軸傾角分別為0,45°和90°(VTI,TTI和HTI)時(shí),qP,qSV與qSH反射波的射線路徑和旅行時(shí)分布。圖3至圖5 是地層對(duì)稱軸傾角為0時(shí),qP,qSV與qSH反射波的射線路徑和旅行時(shí)分布。圖6至圖8是地層對(duì)稱軸傾角為45°時(shí),qP,qSV與qSH反射波的射線路徑和旅行時(shí)分布。圖9至圖11是地層對(duì)稱軸傾角為90°時(shí),qP,qSV與qSH反射波的射線路徑和旅行時(shí)分布。比較這3組各向異性模型射線追蹤結(jié)果可以看出,對(duì)于同一各向異性模型,不同波模式(qP,qSV與qSH波)具有不同的反射波射

圖1　模型內(nèi)部的非均勻節(jié)點(diǎn)網(wǎng)格剖分(白色實(shí)心圓代表主節(jié)點(diǎn);白色空心圈代表次節(jié)點(diǎn);紅色菱形點(diǎn)代表界面節(jié)點(diǎn))

線路徑和旅行時(shí)分布(圖3,圖4和圖5);當(dāng)?shù)貙訉?duì)稱軸傾角不同時(shí),同一波模式的反射波射線路徑和旅行時(shí)分布也不相同(圖3,圖6和圖9)。

圖2　模型內(nèi)部的波前等值線(qP波)

圖3　地層對(duì)稱軸傾角為0時(shí)qP波的射線路徑(a)和旅行時(shí)分布(b)

圖4　地層對(duì)稱軸傾角為0時(shí)qSV波的射線路徑(a)和旅行時(shí)分布(b)

圖5　地層對(duì)稱軸傾角為0時(shí)qSH波的射線路徑(a)和旅行時(shí)分布(b)

圖6　地層對(duì)稱軸傾角為45°時(shí)qP波的射線路徑(a)和旅行時(shí)分布(b)

圖7　地層對(duì)稱軸傾角為45°時(shí)qSV波的射線路徑(a)和旅行時(shí)分布(b)

圖8　地層對(duì)稱軸傾角為45°時(shí)qSH波的射線路徑(a)和旅行時(shí)分布(b)

圖9　地層對(duì)稱軸傾角為90°時(shí)qP波的射線路徑(a)和旅行時(shí)分布(b)

圖10　地層對(duì)稱軸傾角為90°時(shí)qSV波的射線路徑(a)和旅行時(shí)分布(b)

圖11　地層對(duì)稱軸傾角為90°時(shí)qSH波的射線路徑(a)和旅行時(shí)分布(b)

4結(jié)束語

本文研究的非均勻節(jié)點(diǎn)網(wǎng)格TI介質(zhì)射線追蹤技術(shù)可以在較少的網(wǎng)格單元數(shù)情況下,同時(shí)兼顧數(shù)值精度和運(yùn)算效率。首先求取各向異性群速度模型,然后結(jié)合群速度模型和最短路徑射線追蹤算法求取射線路徑和旅行時(shí),最后根據(jù)入射節(jié)點(diǎn)在非均勻節(jié)點(diǎn)網(wǎng)格模型內(nèi)的序號(hào)反向追蹤,求取連接炮點(diǎn)和檢波點(diǎn)之間的射線路徑。由于qP,qSV和qSH波的群速度表達(dá)式不同,因此需要分別求取這3種波模式的射線路徑和旅行時(shí)。對(duì)于同一各向異性模型,不同波模式具有不同的反射波射線路徑和旅行時(shí)。對(duì)于不同對(duì)稱軸傾角模型,同一種波模式的反射波射線路徑和旅行時(shí)也不相同。作為一種正演技術(shù),這種非均勻節(jié)點(diǎn)網(wǎng)格TI介質(zhì)反射波射線追蹤算法在復(fù)雜地層各向異性參數(shù)反演和各向異性地震偏移等領(lǐng)域具有廣泛的應(yīng)用前景。

致謝:感謝阿拉伯聯(lián)合酋長(zhǎng)國石油大學(xué)周兵教授對(duì)本文給予的指導(dǎo)和幫助。

參考文獻(xiàn)

[1]BABICH V M.Ray method for the computation of the intensity of wave fronts in elastic inhomogeneous anisotropic medium:problems of the dynamic theory of propagation of seismic waves[M].Saint Petersburg:Leningrad University Press,1961:36-46

[2]CERVENY V.Seismic ray theory[M].Cambridge:Cambridge University Press,2001:234-400

[3]CERVENY V.Seismic rays and ray intensities in inhomogenous anisotropic media[J].Geophysical Journal of the Royal Astronomical Society,1972,29:1-13

[4]CERVENY V,FIRBAS P.Numerical modeling and inversion of traveltimes of seismic body waves in inhomogeneous anisotropic media[J].Geophysical Journal of the Royal Astronomical Society,1984,76:41-51

[5]GAJEWSKI D,PSENCIK I.Computation of high-frequency seismic wavefields in 3-D laterally inhomogeneous anisotropic media[J].Geophysical Journal of the Royal Astronomical Society,1987,91:383-411

[6]SHEARER P M,CHAPMAN C H.Ray tracing in anisotropic media with a linear gradient[J].Geophysical Journal International,1988,94(3):575-580

[7]張文生,何樵登.橫向各向同性介質(zhì)中的反射波時(shí)距曲線[J].石油物探,1997,36(2):15-24

ZHANG W S,HE Q D.Reflection time-distance curves in transversely isotropic media[J].Geophysical Prospecting for Petroleum,1997,36(2):15-24

[8]苑書金,董敏煜,魏修成,等.橫向各向同性介質(zhì)中的反射波旅行時(shí)分析[J].石油地球物理勘探,2001,36(4):488-494

YUAN S J,DONG M Y,WEI X C,et al.Analysis of reflection travel-time in transversely isotropic media[J].Oil Geophysical Prospecting,2001,36(4):488-494

[9]鄧懷群,劉雯林,趙正茂.橫向各向同性介質(zhì)中縱波和轉(zhuǎn)換橫波的快速射線追蹤方法[J].石油物探,2000,39(4):1-11

DENG H Q,LIU W L,ZHAO Z M.Fast ray tracing method for compressional and conveted waves in transversely isotropic media[J].Geophysical Prospecting for Petroleum,2000,39(4):1-11

[10]孔選林,李錄明,羅省賢,等.各向異性介質(zhì)中地震波射線正演[J].物探化探計(jì)算技術(shù),2008,30(3):178-184

KONG X L,LI L M,LUO S X,et al.Seismic wave ray forward in anisotropy medium[J].Computing Techniques for Geophysical and Geochemical Exploration,2008,30(3):178-184

[11]熊金良,劉洋,侯伯剛.任意各向異性介質(zhì)方位旅行時(shí)正演[J].石油地球物理勘探,2005,40(3):300-304

XIONG J L,LIU Y,HOU B G.Azimuth travel time forward modeling in arbitrary anisotropic media[J].Oil Geophysics Prospecting,2005,40(3):300-304

[12]李建國,李彥鵬,郭曉玲.VTI介質(zhì)試射射線追蹤[J].石油地球物理勘探,2010,45(4):491-496

LI J G,LI Y P,GUO X L.VTI medium test firing ray tracing[J].Oil Geophysics Prospecting,2010,45(4):491-496

[13]郝奇,何樵登,王德利,等.二維層狀VTI介質(zhì)的梯形塊狀建模及qP/qSV反射波運(yùn)動(dòng)學(xué)射線追蹤[J].石油地球物理勘探,2009,44(6):662-670

HAO Q,HE Q D,WANG D L,et al.Trapezoid block modeling and qP/qSV reflection wave kinematic ray tracing in 2D layered VTI media[J].Oil Geophysical Prospecting,2009,44(6):662-670

[14]QIAN J,SYMES W W.Finite-difference quasi-P traveltimes for anisotropic media[J].Geophysics,2002,67(1):147-155

[15]劉玉柱,王光銀,董良國,等.VTI介質(zhì)多參數(shù)聯(lián)合走時(shí)層析成像方法[J].地球物理學(xué)報(bào),2014,57(10):3402-3410

LIU Y Z,WANG G Y,DONG L G,et al.Joint inversion of VTI parameters using nonlinear traveltime tomography[J].Chinese Journal of Geophysics,

2014,57(10):3402-3410

[16]SCHNEIDER W A.Linearization of the P-wave eikonal equation for weak vertical transverse isotropy[J].Geophysics,2003,68(3):1075-1082

[17]ZHOU B,GREENHALGH S A.Ray path and travel time computations for 2D transversely isotropic media with dipping symmetry axes[J].Exploration Geophysics,2006,37(2):150-159

[18]WANG Y H.Seismic ray tracing in anisotropic media:a modified Newton algorithm for solving highly nonlinear systems[J].Geophysics,2014,79(1):T1-T7

[19]趙后越,張美根.起伏地表?xiàng)l件下各向異性地震波最短路徑射線追蹤[J].地球物理學(xué)報(bào),2014,57(9):2910-2917

ZHAO H Y,ZHANG M G.Tracing seismic shortest path rays in anisotropic medium with rolling surfacing[J].Chinese Journal of Geophysics,2014,57(9):2910-2917

[20]QIN F,SCHUSTER G.First-arrival travel time calculation for anisotropic media[J].Geophysics,1993,58(9):1349-1358

[21]WANG X X,TSVANKIN I.Ray-based gridded tomography for tilted transversely isotropic media[J].Geophysics,2013,78(1):C11-C23

[22]KUMAR D,SEN M K,FERGUSON R J.Travel time calculation and prestack depth migration in tilted transversely isotropic media[J].Geophysics,2004,69(1):37-44

[23]PRATT R G,CHAPMAN C H.Travel time tomography in anisotropic media (Ⅱ):application[J].Geophysical Journal International,1992,109(1):20-37

[24]LOU M.Travel time calculation in 3D TTI media by the fast marching method[J].68thEAGE Conference and Exhibition,2006:H034

[25]ZHOU B,GREENHALGH S A.On the computation of elastic wave group velocity for a general anisotropic medium[J].Journal of Geophysics and Engineering,2004,1(3):205-215

[26]CRAMPIN S.A review of wave motion in anisotropic and cracked elastic-media[J].Wave Motion,1981,3(4):343-391

[27]DALEY P,HRON F.Reflection and transmission coefficients for transversely isotropic media[J].Bulletin of the Seismological Society of America,1977,67(3):661-675

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

(編輯：戴春秋)

Reflected wave ray tracing in TI medium based on the nonuniform node meshes

HUANG Guangnan1,2,3,DENG Juzhi1,LI Hongxing1,LI Zelin1,ZHANG Hua1,WANG Andong1

(1.EngineeringResearchCenterofNuclearTechnologyApplication,MinistryofEducation,Nanchang330013,China;2.HubeiSubsurfaceMulti-scaleImagingKeyLaboratory(SMIL),ChinaUniversityofGeosciences,Wuhan430074,China;3.StateKeyLaboratoryofPetroleumResourcesandProspecting,ChinaUniversityofPetroleum,Beijing102249,China)

Abstract:The actual stratum not only has anisotropic property,but also has complicated geological structure to some extent.The efficiency of uniform grid ray tracing method is poor because there are too many cells caused by the refining grid algorithm.The ray tracing algorithm based on the nonuniform node meshes can obtain much higher numerical resolution and computational efficiency even with less primary node number and some secondary node number.Firstly,the calculation of the group velocities for qP,qSV and qSH waves are discussed in this article.Then,the reflected wave ray tracing in TI media is realized by combining these group velocities with the ray tracing algorithm based on the nonuniform node meshes.Finally,three undulating-layered TI models with different symmetry axes are built,the raypaths and traveltime of the qP,qSV and qSH reflected waves are computed for these models respectively.The numerical simulation results suggest that:①the ray tracing algorithm based on the nonuniform node meshes has applicability for such a complex TI media;②three wave-modes have different raypaths and traveltimes for the same model;③the same wave-mode also has different raypaths and traveltime for the model with different symmetry axes.

Keywords:undulating-layered model,nonuniform node meshes,anisotropic TI media,ray tracing algorithm

文章編號(hào)：1000-1441(2016)01-0025-08

DOI:10.3969/j.issn.1000-1441.2016.01.004

中圖分類號(hào)：P631

文獻(xiàn)標(biāo)識(shí)碼：A

基金項(xiàng)目：國家科技重大專項(xiàng)(2011ZX05024-001-02)、國家自然科學(xué)基金(41504095,41004048,41364004,41104074,41304097)、國家科技支撐計(jì)劃(2011BAB04B03)、東華理工大學(xué)博士科研啟動(dòng)基金(DHBK2013212)、核技術(shù)應(yīng)用教育部工程研究中心基金(HJSJYB2015-9)、中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室基金(SMIL-2015-10)和江西省教育廳基金(GJJ14476)項(xiàng)目聯(lián)合資助。

作者簡(jiǎn)介：黃光南(1983—),男,博士,講師,主要從事地震速度層析成像和地震數(shù)字處理方法研究。

收稿日期：2015-05-04;改回日期：2015-09-27。

黃光南,鄧居智,李紅星,等.非均勻節(jié)點(diǎn)網(wǎng)格TI介質(zhì)反射波射線追蹤研究[J].石油物探,2016,55(1):-32

HUANG Guangnan,DENG Juzhi,LI Hongxing,et al.Reflected wave ray tracing in TI medium based on the nonuniform node meshes[J].Geophysical Prospecting for Petroleum,2016,55(1):-32

This research is financially supported by the National Science and Technology Major Project of China (Grant No.2011ZX05024-001-02),National Natural Science Foundation of China (Grant Nos.41504095,41004048,41364004,41104074,41304097),National Science and Technology Supported Program (Grant No.2011BAB04B03),Doctoral Research Foundation of East China University of Technology (Grant No.DHBK2013212),Foundation of the Engineering Research Center of Nuclear Technology Application,Ministry of Education (Grant No.HJSJYB2015-9),Foundation of Hubei Subsurface Multi-scale Imaging Key Laboratory,China University of Geosciences (Grant No.SMIL-2015-10),Foundation of Education Department of Jiangxi Province (Grant No.GJJ14476).