Ronge－Kutta法在計(jì)算復(fù)雜構(gòu)造區(qū)域地震波旅行時(shí)中的應(yīng)用

鄒長(zhǎng)橋 段永紅 唐秋惠 徐 佼 陳立波

1 中國(guó)地震局地球物理勘探中心，鄭州市文化路75號(hào)，450002

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

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

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

1 Ronge－Kutta法

1.1 一階常微分方程初值問(wèn)題的求解

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

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

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

其中，

1.2 Ronge－Kuttta法解射線方程

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

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

其中，

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

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

2 實(shí) 例

速度模型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)包含一個(gè)高速圓形異常體。如圖1（b）所示，炮點(diǎn)位于（0，0），每條射線的出射角已知，觀測(cè)系統(tǒng)位于地表。圖1中的射線路徑是采用Ronge－Kutta射線追蹤法得到的。由圖1（b）可以看出，有些初至波的射線軌跡發(fā)生了彎曲、回折現(xiàn)象，對(duì)應(yīng)為地震回折波。這些回折波也表明模型1在第3層的低速層內(nèi)包含一個(gè)圓形高速異常體，這就驗(yàn)證了模型1構(gòu)建時(shí)的猜想是正確的，同時(shí)也驗(yàn)證了Ronge－Kutta射線追蹤法在復(fù)雜構(gòu)造模型中是可行的。

圖2（a）為光滑Marmousi模型（模型2）?？梢钥闯?，該模型的速度與深度成非線性，即該模型的介質(zhì)為不連續(xù)變化介質(zhì)，但底部為高速介質(zhì)。該速度模型主要有3個(gè)不連續(xù)的起伏界面。炮點(diǎn)位于（0，0）處，觀測(cè)系統(tǒng)位于地面，每條射線的出射角已知。由圖2（b）得到的Ronge－Kutta射線追蹤可知，地震波射線軌跡在均勻介質(zhì)中為圓弧形，而在異常體左側(cè)發(fā)生了回折現(xiàn)象，對(duì)應(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個(gè)斷層，在低速層內(nèi)包含有高速介質(zhì)，底部的高速介質(zhì)中也包含有低速介質(zhì)。圖3（b）為炮點(diǎn)在（0，0）處、觀測(cè)系統(tǒng)位于地表的射線追蹤圖。從初至波的射線軌跡可以看出，模型3具有3個(gè)主反射界面，第1個(gè)界面在1 300m 處，第2個(gè)界面在1 600m 處，第3個(gè)界面在2 400m 處。地震波在第3個(gè)界面發(fā)生了彎曲、透射現(xiàn)象，表明第3個(gè)界面速度不均勻。

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

3 結(jié)果對(duì)比

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

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

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

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

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

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

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

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

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

圖7 旅行時(shí)等時(shí)線Fig.7 The isochron diagram of travel time

圖8 程函方程計(jì)算旅行時(shí)Fig.8 The eikonal equation calculated travel times

4 結(jié) 語(yǔ)

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

［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］孫章慶，孫建國(guó)，韓復(fù)興.針對(duì)復(fù)雜地形的三種地震波走時(shí)算法及對(duì)比［J］.地球物理學(xué)報(bào)，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ìn)的利用程函方程計(jì)算旅行時(shí)的方法［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］陳可洋.地震波旅行時(shí)計(jì)算方法及其模型試驗(yàn)分析［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］孟憲海，金穎，李吉?jiǎng)?，?基于三角域快速行進(jìn)法的地震波走時(shí)計(jì)算［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］彭直興，沈忠民.基于雙線性插值的三維地震波旅行時(shí)計(jì)算［J］.西南石油大學(xué)學(xué)報(bào)：自然科學(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］.西安：長(zhǎng)安大學(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］李強(qiáng)，白超英.復(fù)雜介質(zhì)中地震波前及射線追蹤綜述［J］.地球物 理 學(xué) 進(jìn) 展，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］朱金明，王麗燕.地震波走時(shí)的有限差分法計(jì)算［J］.地球物理學(xué)報(bào)，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］蘭海強(qiáng)，張智，徐濤，等.地震波走時(shí)場(chǎng)模擬的快速推進(jìn)法和快速掃描法比較研究［J］.地球物理學(xué)進(jìn)展，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］楊昊.有限差分法地震波走時(shí)計(jì)算的快速算法研究［D］.長(zhǎng)春：吉林大學(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ì)算技術(shù)［J］.地震學(xué)報(bào)，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