Ronge－Kutta法在计算复杂构造区域地震波旅行时中的应用

邹长桥 段永红 唐秋惠 徐 佼 陈立波

1 中国地震局地球物理勘探中心，郑州市文化路75号，450002

2 桂林理工大学地球科学学院，桂林市建干路12号，541004

地震波旅行时信息在地震勘探和地震学领域占有重要的位置，在叠前偏移、叠前速度分析、地震走时层析成像及地震定位方面应用广泛［1－2］。计算初至波走时的方法主要有4类：基于高频近似射线理论的最短路径方法、基于程函方程的数值解方法、基于惠更斯原理的波前构建法和频率域波动方程法［3－10］。程函方程数值解法不需要计算射线路径，且计算速度较高，具有较好的稳定性，但是计算精度比较低，不少学者通过引入高阶差分来提高计算精度。波前构建法弥补了程函方程数值解法精度低的问题，稳定性好，但需要在射线网格和规则网格之间互相转换，从而影响了计算效率。频率域波动方程法虽然能适应任意复杂介质模型的计算，但计算精度和效率仍很低。传统的射线方法计算精度较高，稳定性也较好，但是需要应用较多的网格节点，影响了计算效率［11－12］。在构造复杂区域常常会遇到成像效果差或基本不成像的难题，主要原因是大多数地震成像技术和地震波理论都是基于水平地表和层状介质的简单模型得出，使得常规的数据采集、处理和解释技术在复杂构造区域中不再适用［13－14］。

Ronge－Kutta射线追踪法利用4 阶Ronge－Kutta法解射线追踪方程，采用低阶导数的线性组合和逼近高阶导数的方法，避免了求解高阶导数的巨大计算量，能取得很高的计算精度［15－16］。本文拟采用Ronge－Kutta法计算构造复杂区域的地震波旅行时。

1 Ronge－Kutta法

1.1 一阶常微分方程初值问题的求解

设有一阶常微分方程组的初值问题：

根据Ronge－Kutta法的一般式：

其中bi、ci、aij为待定常数，c1＝0，a1j＝0，i＝1，2，…，s－1。由泰勒级数和待定系数法可以得到4阶Ronge－Kutta法形式为：

其中，

1.2 Ronge－Kuttta法解射线方程

给定一个速度模型的初始条件，例如炮点坐标（x0，z0）和射线出射角θ0，根据Ronge－Kuttta法求解射线方程组，其中θ为入射角与z轴间的夹角。假设速度模型选择的参数化方式为梯形网格单元化，则在射线追踪过程中，当射线近似水平方向时，采用下式：

其中，vx为x方向速度，vz为z方向速度。当射线近似垂直方向时，采用下式：

其中，

其中，4阶经典Ronge－Kuttta公式为：

通过式（8）即可求解式（5）、（6），从而得出地震波的射线路径和走时。

2 实 例

速度模型1（图1（a））为层状速度模型，模型速度从上至下随深度递增（连续介质）大致分为4层：第1层介质速度为1 000m／s，第2层为3 000 m／s，第3层为5 000m／s，第4层为7 500m／s，模型在低速层内包含一个高速圆形异常体。如图1（b）所示，炮点位于（0，0），每条射线的出射角已知，观测系统位于地表。图1中的射线路径是采用Ronge－Kutta射线追踪法得到的。由图1（b）可以看出，有些初至波的射线轨迹发生了弯曲、回折现象，对应为地震回折波。这些回折波也表明模型1在第3层的低速层内包含一个圆形高速异常体，这就验证了模型1构建时的猜想是正确的，同时也验证了Ronge－Kutta射线追踪法在复杂构造模型中是可行的。

图2（a）为光滑Marmousi模型（模型2）。可以看出，该模型的速度与深度成非线性，即该模型的介质为不连续变化介质，但底部为高速介质。该速度模型主要有3个不连续的起伏界面。炮点位于（0，0）处，观测系统位于地面，每条射线的出射角已知。由图2（b）得到的Ronge－Kutta射线追踪可知，地震波射线轨迹在均匀介质中为圆弧形，而在异常体左侧发生了回折现象，对应为回折波，右侧发生了弯曲、交叉现象，则表明该异常体为圆形高速异常体。

图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）的速度变化复杂，速度不随深度连续变化，且没有明显的分层界面。该模型有3个断层，在低速层内包含有高速介质，底部的高速介质中也包含有低速介质。图3（b）为炮点在（0，0）处、观测系统位于地表的射线追踪图。从初至波的射线轨迹可以看出，模型3具有3个主反射界面，第1个界面在1 300m 处，第2个界面在1 600m 处，第3个界面在2 400m 处。地震波在第3个界面发生了弯曲、透射现象，表明第3个界面速度不均匀。

图3 Marmousi模型及Ronge－Kutta射线追踪图Fig.3 The Marmousi model and ray tracing of the Ronge－Kutta

3 结果对比

针对前面建立的速度模型，分别作出时距曲线图、等时线图和旅行时图，对比Ronge－Kutta射线追踪法和程函方程有限差分旅行时算法。

一般来说，在进行数值计算时，步长越小、网格越密，计算效果越好。为此，本文按照空间采样定理分别进行z方向和x方向加密。在得到各自方向的网格加密步长后，以不同的权系数进行z轴和x轴的同步加密，将理论走时与加密后得到的走时结果进行比较，最接近于理论走时的网格尺度为最优网格步长。图4为用不同的网格参数对模型1进行网格化后的时距曲线图，实线是程函方程计算的理论值，虚线是不同网格参数计算的值，红色虚线网格步长为0.05，蓝色虚线为0.01，黑色虚线为0.005。图4（a）表示在z轴加密，红色网格数为11×101，蓝色网格数为51×101，黑色网格数为101×101。从图中可以看出，黑色虚线比较接近理论值，因此，z轴加密时，应选择网格步长为0.005和网格数为101×101的网格参数。图4（b）表示在x轴加密，红色网格数为51×21，蓝色网格数为51×101，白色网格数为51×201。从图中可以看出，蓝色虚线比较接近理论值，因此，对x轴加密时，应选择网格步长为0.01和网格数为51×101的网格参数。

图5为速度模型1同时在z轴和x轴网格加密的时距曲线图。实线是程函方程计算得到的理论值，虚线是不同网格参数的计算值。红色虚线的网格步长为0.05，网格数为21×51；蓝色虚线的网格步长为0.01，网格数为51×101；黑色虚线的网格步长为0.005，网格数为101×201。可以看出，黑色虚线比较接近实线理论值，所以对速度模型1进行网格x轴和z轴加密时，应选择网格步长为0.005和网格数101×201的网格参数。

图4 单坐标轴加密的时距曲线Fig.4 The hodograph of encryption in the single Axis

图5 双坐标轴加密的时距曲线Fig.5 The hodograph of encryption in the double axis

图6为模型2 的初至波时距曲线图，3 条曲线表示炮点位于不同深度。蓝色点线是具有因果关系的双平方根求解程函方程得出的结果，黑色实线为有限差分求解程函方程得出的结果。

图6 初至波时距曲线Fig.6 The hodograph of preliminary wave

图7为模型1的旅行时等时线，因为等时线与射线路径是相互垂直的，所以图7的等时线与图1～3的射线路径是相互对应的。从图7也可以看出，模型1的第3层低速层内包含着一个高速异常体。

图8为程函方程计算出的地震波旅行时，炮点位于（0，0）处。可以看出波前在介质中的传播情况，而且从波前的形状可以看出，该模型的中部有异常体。

图7 旅行时等时线Fig.7 The isochron diagram of travel time

图8 程函方程计算旅行时Fig.8 The eikonal equation calculated travel times

4 结 语

通过对几个模型地震初至波和反射波的Ronge－Kutta法射线追踪数值模拟的研究和实例分析，以及Ronge－Kutta射线追踪模型和程函方程计算的时距曲线效果对比，验证了Ronge－Kutta射线追踪法的易实现性和程函方程有限差分旅行时算法的强稳定性。Rong－Kutta射线追踪法在已知出射方向和倾角大、埋藏深、断层和裂隙孔洞分布无规律、速度变化剧烈等构造复杂区域进行地震初至波射线追踪时，能得到较好的模拟效果，精度较高，但是稳定性和计算效率不及程函方程有限差分法。因此，建议在复杂构造地震波旅行时计算中应尽量同时使用两种方法。

［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］孙章庆，孙建国，韩复兴.针对复杂地形的三种地震波走时算法及对比［J］.地球物理学报，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］.西南石油大学学报：自然科学版，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］.西安：长安大学，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］李强，白超英.复杂介质中地震波前及射线追踪综述［J］.地球物 理 学 进 展，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］.地球物理学报，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］.地球物理学进展，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］.长春：吉林大学，2007（Yang Hao.Study on the Fast Computation Technique of Seismic Travel Times with Finite－Difference Method［D］.Changchun：Jilin University，2007）

［14］赵爱华，张中杰，王光杰，等.非均匀介质中地震波走时与射线路径快速计算技术［J］.地震学报，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