鄒長橋 段永紅 唐秋惠 徐 佼 陳立波

1 中國地震局地球物理勘探中心，鄭州市文化路75號，450002

2 桂林理工大學(xué)地球科學(xué)學(xué)院，桂林市建干路12號，541004

地震波旅行時信息在地震勘探和地震學(xué)領(lǐng)域占有重要的位置，在疊前偏移、疊前速度分析、地震走時層析成像及地震定位方面應(yīng)用廣泛［1－2］。計算初至波走時的方法主要有4類：基于高頻近似射線理論的最短路徑方法、基于程函方程的數(shù)值解方法、基于惠更斯原理的波前構(gòu)建法和頻率域波動方程法［3－10］。程函方程數(shù)值解法不需要計算射線路徑，且計算速度較高，具有較好的穩(wěn)定性，但是計算精度比較低，不少學(xué)者通過引入高階差分來提高計算精度。波前構(gòu)建法彌補了程函方程數(shù)值解法精度低的問題，穩(wěn)定性好，但需要在射線網(wǎng)格和規(guī)則網(wǎng)格之間互相轉(zhuǎn)換，從而影響了計算效率。頻率域波動方程法雖然能適應(yīng)任意復(fù)雜介質(zhì)模型的計算，但計算精度和效率仍很低。傳統(tǒng)的射線方法計算精度較高，穩(wěn)定性也較好，但是需要應(yīng)用較多的網(wǎng)格節(jié)點，影響了計算效率［11－12］。在構(gòu)造復(fù)雜區(qū)域常常會遇到成像效果差或基本不成像的難題，主要原因是大多數(shù)地震成像技術(shù)和地震波理論都是基于水平地表和層狀介質(zhì)的簡單模型得出，使得常規(guī)的數(shù)據(jù)采集、處理和解釋技術(shù)在復(fù)雜構(gòu)造區(qū)域中不再適用［13－14］。

Ronge－Kutta射線追蹤法利用4 階Ronge－Kutta法解射線追蹤方程，采用低階導(dǎo)數(shù)的線性組合和逼近高階導(dǎo)數(shù)的方法，避免了求解高階導(dǎo)數(shù)的巨大計算量，能取得很高的計算精度［15－16］。本文擬采用Ronge－Kutta法計算構(gòu)造復(fù)雜區(qū)域的地震波旅行時。

1 Ronge－Kutta法

1.1 一階常微分方程初值問題的求解

設(shè)有一階常微分方程組的初值問題：

根據(jù)Ronge－Kutta法的一般式：

其中bi、ci、aij為待定常數(shù)，c1＝0，a1j＝0，i＝1，2，…，s－1。由泰勒級數(shù)和待定系數(shù)法可以得到4階Ronge－Kutta法形式為：

其中，

1.2 Ronge－Kuttta法解射線方程

給定一個速度模型的初始條件，例如炮點坐標(biāo)（x0，z0）和射線出射角θ0，根據(jù)Ronge－Kuttta法求解射線方程組，其中θ為入射角與z軸間的夾角。假設(shè)速度模型選擇的參數(shù)化方式為梯形網(wǎng)格單元化，則在射線追蹤過程中，當(dāng)射線近似水平方向時，采用下式：

其中，vx為x方向速度，vz為z方向速度。當(dāng)射線近似垂直方向時，采用下式：

其中，

其中，4階經(jīng)典Ronge－Kuttta公式為：

通過式（8）即可求解式（5）、（6），從而得出地震波的射線路徑和走時。

2 實 例

速度模型1（圖1（a））為層狀速度模型，模型速度從上至下隨深度遞增（連續(xù)介質(zhì)）大致分為4層：第1層介質(zhì)速度為1 000m／s，第2層為3 000 m／s，第3層為5 000m／s，第4層為7 500m／s，模型在低速層內(nèi)包含一個高速圓形異常體。如圖1（b）所示，炮點位于（0，0），每條射線的出射角已知，觀測系統(tǒng)位于地表。圖1中的射線路徑是采用Ronge－Kutta射線追蹤法得到的。由圖1（b）可以看出，有些初至波的射線軌跡發(fā)生了彎曲、回折現(xiàn)象，對應(yīng)為地震回折波。這些回折波也表明模型1在第3層的低速層內(nèi)包含一個圓形高速異常體，這就驗證了模型1構(gòu)建時的猜想是正確的，同時也驗證了Ronge－Kutta射線追蹤法在復(fù)雜構(gòu)造模型中是可行的。

圖2（a）為光滑Marmousi模型（模型2）。可以看出，該模型的速度與深度成非線性，即該模型的介質(zhì)為不連續(xù)變化介質(zhì)，但底部為高速介質(zhì)。該速度模型主要有3個不連續(xù)的起伏界面。炮點位于（0，0）處，觀測系統(tǒng)位于地面，每條射線的出射角已知。由圖2（b）得到的Ronge－Kutta射線追蹤可知，地震波射線軌跡在均勻介質(zhì)中為圓弧形，而在異常體左側(cè)發(fā)生了回折現(xiàn)象，對應(yīng)為回折波，右側(cè)發(fā)生了彎曲、交叉現(xiàn)象，則表明該異常體為圓形高速異常體。

圖1 速度模型1及Ronge－Kutta射線追蹤圖Fig.1 The first velocity model and ray tracing of the Ronge－Kutta

圖2 速度模型2及Ronge－Kutta射線追蹤圖Fig.2 The second velocity model and ray tracing of the Ronge－Kutta

圖3（a）所示Marmousi模型（模型3）的速度變化復(fù)雜，速度不隨深度連續(xù)變化，且沒有明顯的分層界面。該模型有3個斷層，在低速層內(nèi)包含有高速介質(zhì)，底部的高速介質(zhì)中也包含有低速介質(zhì)。圖3（b）為炮點在（0，0）處、觀測系統(tǒng)位于地表的射線追蹤圖。從初至波的射線軌跡可以看出，模型3具有3個主反射界面，第1個界面在1 300m 處，第2個界面在1 600m 處，第3個界面在2 400m 處。地震波在第3個界面發(fā)生了彎曲、透射現(xiàn)象，表明第3個界面速度不均勻。

圖3 Marmousi模型及Ronge－Kutta射線追蹤圖Fig.3 The Marmousi model and ray tracing of the Ronge－Kutta

3 結(jié)果對比

針對前面建立的速度模型，分別作出時距曲線圖、等時線圖和旅行時圖，對比Ronge－Kutta射線追蹤法和程函方程有限差分旅行時算法。

一般來說，在進行數(shù)值計算時，步長越小、網(wǎng)格越密，計算效果越好。為此，本文按照空間采樣定理分別進行z方向和x方向加密。在得到各自方向的網(wǎng)格加密步長后，以不同的權(quán)系數(shù)進行z軸和x軸的同步加密，將理論走時與加密后得到的走時結(jié)果進行比較，最接近于理論走時的網(wǎng)格尺度為最優(yōu)網(wǎng)格步長。圖4為用不同的網(wǎng)格參數(shù)對模型1進行網(wǎng)格化后的時距曲線圖，實線是程函方程計算的理論值，虛線是不同網(wǎng)格參數(shù)計算的值，紅色虛線網(wǎng)格步長為0.05，藍色虛線為0.01，黑色虛線為0.005。圖4（a）表示在z軸加密，紅色網(wǎng)格數(shù)為11×101，藍色網(wǎng)格數(shù)為51×101，黑色網(wǎng)格數(shù)為101×101。從圖中可以看出，黑色虛線比較接近理論值，因此，z軸加密時，應(yīng)選擇網(wǎng)格步長為0.005和網(wǎng)格數(shù)為101×101的網(wǎng)格參數(shù)。圖4（b）表示在x軸加密，紅色網(wǎng)格數(shù)為51×21，藍色網(wǎng)格數(shù)為51×101，白色網(wǎng)格數(shù)為51×201。從圖中可以看出，藍色虛線比較接近理論值，因此，對x軸加密時，應(yīng)選擇網(wǎng)格步長為0.01和網(wǎng)格數(shù)為51×101的網(wǎng)格參數(shù)。

圖5為速度模型1同時在z軸和x軸網(wǎng)格加密的時距曲線圖。實線是程函方程計算得到的理論值，虛線是不同網(wǎng)格參數(shù)的計算值。紅色虛線的網(wǎng)格步長為0.05，網(wǎng)格數(shù)為21×51；藍色虛線的網(wǎng)格步長為0.01，網(wǎng)格數(shù)為51×101；黑色虛線的網(wǎng)格步長為0.005，網(wǎng)格數(shù)為101×201。可以看出，黑色虛線比較接近實線理論值，所以對速度模型1進行網(wǎng)格x軸和z軸加密時，應(yīng)選擇網(wǎng)格步長為0.005和網(wǎng)格數(shù)101×201的網(wǎng)格參數(shù)。

圖4 單坐標(biāo)軸加密的時距曲線Fig.4 The hodograph of encryption in the single Axis

圖5 雙坐標(biāo)軸加密的時距曲線Fig.5 The hodograph of encryption in the double axis

圖6為模型2 的初至波時距曲線圖，3 條曲線表示炮點位于不同深度。藍色點線是具有因果關(guān)系的雙平方根求解程函方程得出的結(jié)果，黑色實線為有限差分求解程函方程得出的結(jié)果。

圖6 初至波時距曲線Fig.6 The hodograph of preliminary wave

圖7為模型1的旅行時等時線，因為等時線與射線路徑是相互垂直的，所以圖7的等時線與圖1～3的射線路徑是相互對應(yīng)的。從圖7也可以看出，模型1的第3層低速層內(nèi)包含著一個高速異常體。

圖8為程函方程計算出的地震波旅行時，炮點位于（0，0）處?？梢钥闯霾ㄇ霸诮橘|(zhì)中的傳播情況，而且從波前的形狀可以看出，該模型的中部有異常體。

圖7 旅行時等時線Fig.7 The isochron diagram of travel time

圖8 程函方程計算旅行時Fig.8 The eikonal equation calculated travel times

4 結(jié) 語

通過對幾個模型地震初至波和反射波的Ronge－Kutta法射線追蹤數(shù)值模擬的研究和實例分析，以及Ronge－Kutta射線追蹤模型和程函方程計算的時距曲線效果對比，驗證了Ronge－Kutta射線追蹤法的易實現(xiàn)性和程函方程有限差分旅行時算法的強穩(wěn)定性。Rong－Kutta射線追蹤法在已知出射方向和傾角大、埋藏深、斷層和裂隙孔洞分布無規(guī)律、速度變化劇烈等構(gòu)造復(fù)雜區(qū)域進行地震初至波射線追蹤時，能得到較好的模擬效果，精度較高，但是穩(wěn)定性和計算效率不及程函方程有限差分法。因此，建議在復(fù)雜構(gòu)造地震波旅行時計算中應(yīng)盡量同時使用兩種方法。

［1］Cerveny V.Seismic Ray Theory［M］.Cambridge：Cambridge University Press，2001

［2］Zhang Z J，Wang G J，Teng J W，et al.CDP Mapping to Obtain the Fine Structure of the Crust and Upper Mantle from Seismic Sounding Data：An Example for the Southeastern China［J］.Physics of the Earth and Planetary Interiors，2000，122（1－2）：133－146

［3］孫章慶，孫建國，韓復(fù)興.針對復(fù)雜地形的三種地震波走時算法及對比［J］.地球物理學(xué)報，2012，55（2）：560－568（Sun Zhangqing，Sun Jianguo，Han Fuxing.The Comparison of Three Schemes for Computing Seismic Wave Travel Times in Complex Topographical Conditions［J］.Chinese Journal of Geophysics，2012，55（2）：560－568）

［4］李文杰，魏修成，寧俊瑞，等.一種改進的利用程函方程計算旅行時的方法［J］.石油地球物理勘探，2008，43（5）：589－594（Li Wenjie，Wei Xiucheng，Ning Junrui，et al.A Method Using Improved Eikonal Equation to Compute Travel Time［J］.Oil Geophysical Prospecting，2008，43（5）：589－594）

［5］陳可洋.地震波旅行時計算方法及其模型試驗分析［J］.石油物 探，2010，49（2）：154－157（Chen Keyang.Compution Method for Seismic Wave Travel－Time and Its Model Experiment Analysis［J］.GPP，2010，49（2）：154－157）

［6］孟憲海，金穎，李吉剛，等.基于三角域快速行進法的地震波走時計算［J］.軟件，2011，32（11）：36－42（Meng Xianhai，Jin Yin，Li Jigang，et al.Seismic Travel－Time Computation Using Triangulated Fast Marching Method［J］.Software，2011，32（11）：36－42）

［7］彭直興，沈忠民.基于雙線性插值的三維地震波旅行時計算［J］.西南石油大學(xué)學(xué)報：自然科學(xué)版，2008，30（5）：85－92（Peng Zhixing，Shen Zhongmin.Computation of Three－Dimensional Seismic Travel Time Based on Bilinear Interpolation［J］.Journal of Southwest Petroleum University：Science ＆Technology Edition，2008，30（5）：85－92）

［8］陳超群.逐段迭代法射線追蹤三維地震道集記錄正演模擬［D］.西安：長安大學(xué)，2007（Chen Chaoqun.The Forward Simulation of Iteration Segment by Segment Ray－Tracing in 3D Seismic Gather［D］.Xi’an：Chang’an University，2007）

［9］李強，白超英.復(fù)雜介質(zhì)中地震波前及射線追蹤綜述［J］.地球物 理 學(xué) 進 展，2012，27（1）：92－104（Li Qiang，Bai Chaoying.Review on Seismic Wave Front and Ray Tracing in Complex Media［J］.Progress in Geophysics，2012，27（1）：92－104）

［10］朱金明，王麗燕.地震波走時的有限差分法計算［J］.地球物理學(xué)報，1992，35（1）：86－92（Zhu Jinming，Wang Liyan.Finite Difference Calculation of Seismic Travel Times［J］.Chinese Journal of Geophysics，1992，35（1）：86－92）

［11］蘭海強，張智，徐濤，等.地震波走時場模擬的快速推進法和快速掃描法比較研究［J］.地球物理學(xué)進展，2012，27（5）：1 863－1 870（Lan Haiqiang，Zhang Zhi，Xu Tao，et al.A Comparative Study on the Fast Marching and Fast Sweeping Methods in the Calculation of First－Arrival Travel Time Field［J］.Progress in Geophysics，2012，27（5）：1 863－1 870）

［12］Fomel S.Travel Time Computation with the Linearized Eikonal Equation［J］.SEP Report，1997，94：123－131

［13］楊昊.有限差分法地震波走時計算的快速算法研究［D］.長春：吉林大學(xué)，2007（Yang Hao.Study on the Fast Computation Technique of Seismic Travel Times with Finite－Difference Method［D］.Changchun：Jilin University，2007）

［14］趙愛華，張中杰，王光杰，等.非均勻介質(zhì)中地震波走時與射線路徑快速計算技術(shù)［J］.地震學(xué)報，2000，22（2）：151－157（Zhao Aihua，Zhang Zhongjie，Wang Guangjie，et al.The Fast Computation Technique of Seismic Wave Travel Times and Ray Paths in Inhomogeneous Media［J］.Acta Seismologica Sinica，2000，22（2）：151－157）

［15］Ettrich N，Gajewski D.Wave Front Construction in Smooth Media for Prestack Depth Migration［J］.Pure and Applied Geophysics，1996，148（3／4）：481－502

［16］Bulant P，Klime?L.Interpolation of Ray Theory Traveltimes within Ray Cells［J］.Geophysical Journal International，1999，139（2）：273－282