界面二次源波前扩展法全局最小走时射线追踪技术

以Moser方法为代表的最短路径射线追踪算法可以快速稳定地获得整个追踪区域的全局最小走时和路径，但它存在两个缺陷：一是射线大多由折线呈锯齿状相连，长度和位置偏离真实射线路径；二是在低变速区容易出现射线路径多值现象．本文提出的界面二次源波前扩展法全局最小走时射线追踪技术（以下简称界面源法）旨在解决上述两个问题．不同于Moser方法，界面源法只在物性分界面上设置子波源点，子波出射射线可以到达任何不穿越物性界面而直接到达的空间点和界面离散点，在均匀块体内或层内地震波以精确的射线路径传播．显然，界面源法的子波出射方向数远远大于传统方法，算法的追踪误差主要由界面离散引起的，因此，界面源法很好地解决了Moser法存在的问题，大大提高了追踪的精度．同时，由于界面源法的子波源点数远远小于Moser法，因而效率也很高．模型实算证实了该算法的高效性．     

基于分区多步快速行进法下3D起伏层状介质中多震相射线追踪

快速行进法(FMM)是一种求解程函方程数值解计算网格节点走时,然后向后处理进行射线追踪的方法.为了求取任意起伏界面下高精度多震相的地震走时与相应的射线路径,本文采用任意起伏地表条件下的的三维不等距上行差分公式结合分区多步计算技术实现了三维复杂层状起伏介质中多震相(透射、反射、转换波)地震走时的计算,利用上行有限差分公式逐次进行射线路径的追踪,并且通过与较为成熟的不规则最短路径法(ISPM)对比,验证了本算法的计算精度和有效性.数值模拟实例和对比结果表明该算法具有较高的计算精度,数值计算稳健,能灵活处理含任意三维起伏界面模型中多震相地震走时及相应射线路径的追踪问题.     

计算最小走时和射线路径的界面网全局方法

用慢度分块均匀正方形模型将介质参数化，仅在正方形单元的边界上设置计算结点，这些结点构成界面网．根据Ｈｕｖｓｅｎｓ和Ｆｅｒｍａｔ原理，由不断扩张、收缩的波前点扫描代替波前面搜索，在波前点附近点的局部最小走时计算中对波前点之间的走时使用双曲线近似，通过比较确定最小走时和相应的次级源位置，记录在以界面网点位置为指针的３个一维数组中．借助这些数组通过向源搜索可计算任意点（包括界面网以外的点）上的全局最小走时和射线路径．这一方法不受介质慢度差异大小限制，占内存少，计算速度较快，适于走时反演和以Ｍａｓｌｏｖ射线理论为基础的波场计算．     

三维复杂介质中转换波走时快速计算

复杂介质中转换波走时计算是多波勘探地震学中重要内容之一.本项研究利用惠更斯原理和费玛原理,获得了三维复杂介质中转换波快速计算的改进型最小走时树方法.其中,在保证精度的条件下,为了提高三维转换波走时计算效率,首先对初至波最小走时树基本算法进行了改进.本方法通过将转换波分为上、下行波分别进行射线追踪以实现三维转换波走时的快速计算.模型计算表明,方法的计算速度快,而且稳定性强,对多波地震勘探具有较大的应用价值.     

分区多步最短路径极值法多值多次反射波追踪

基于网格单元扩展的射线追踪算法,如:较为流行的有限差分解程函方程法和最短路径法均是建立在费马(最小走时)原理基础上的射线追踪算法,只能进行单值(最小走时)的多次反射波的追踪.然而当介质速度反差较大或存在复杂反射界面(如:常见的向斜、背斜、透镜体、塌陷构造等)时,将出现波前的自相缠绕(即蝴蝶结现象),相应的地震射线则为多...     

波前构建法研究现状

波前构建法是一种能够快速计算地震波走时、射线路径及振幅的算法,该算法考虑从整个波场发出的射线族的追踪计算,为几何射线理论在勘探地震学中的应用带来了决定性的变化,使得能够同时计算成千上万务射线以及附在上面的走时和振幅成为可能.本文在收集、整理国内外研究现状的基础上,对波前构建法目前的应用介质模型、求解射线追踪方程组的算法...     

一种改进的线性走时插值射线追踪算法

线性走时插值法（LTI）在走时的计算中,由于射线方向考虑不全,计算得到的节点走时不一定最小,导致追踪的射线路径无法满足最小走时.针对这一问题,本文提出了一种改进的射线追踪算法,通过采用多方向的循环计算,得到所有计算节点的最小走时,使追踪到的射线路径能真正满足最小走时,以确保射线追踪的精度.模拟实验结果表明,在介质速度变化剧烈的结构中,该算法与传统的LTI算法相比,有效地提高了射线追踪的精度.     

有序波前重建法的射线追踪

建立了一种新的计算最小走时和射线路径的方法——有序波前重建法. 文中算法按照波前面的实际扩展顺序外推计算走时，采用以计算点为中心的走时计算策略，直接记录计算点获取最小走时的前一节点坐标，同步计算最小走时和射线路径，得到一种全局算法. 该方法具有原理简单、易于实现、不受介质速度差异大小限制、计算速度快等优点. 数值实验表明有序波前重建法具有较高的计算精度和运行效率.     

利用三维弯曲界面的反演方法重建唐山震区莫霍界面三维构造形态

研究利用反射波的走时反演三维弯曲界面和介质层速度的计算方法．各层界面利用分段非完全三次多项式描述．正问题采用一种快速的三维射线追踪方法．反演过程采用变阻尼最小二乘法．数值模拟结果表明，解很快收敛到真模型．处理了通过唐山震区的实测资料，重建了该震区莫霍界面三维构造形态，并揭示了该区域构造与地震活动的关系．      

基于图形结构的三维射线追踪方法

在地震层析成像研究中,为了克服最小走时射线路径追踪方法存在的问题,对该方法计算过程中的关键步骤进行了改进．在节点走时的计算中引入Ｂｒｅｓｅｎｈａｍ画线算法;在最小走时节点查寻中,结合使用快速排序算法与插入排序算法,替代以往方法中多采用的堆排序算法;所采用的节点设置方式,可以引入速度界面,还可以实现反射波射线追踪．模型计算证明,改进的最小走时射线路径方法具有精度高,速度快的特点,所提出的三维空间反射波射线追踪算法简便易行。     

A new 3-D ray tracing method based on LTI using successive partitioning of cell interfaces and traveltime gradients

We present a new method of three-dimensional (3-D) seismic ray tracing, based on an improvement to the linear traveltime interpolation (LTI) ray tracing algorithm. This new technique involves two separate steps. The first involves a forward calculation based on the LTI method and the dynamic successive partitioning scheme, which is applied to calculate traveltimes on cell boundaries and assumes a wavefront that expands from the source to all grid nodes in the computational domain. We locate several dynamic successive partition points on a cell's surface, the traveltimes of which can be calculated by linear interpolation between the vertices of the cell's boundary. The second is a backward step that uses Fermat's principle and the fact that the ray path is always perpendicular to the wavefront and follows the negative traveltime gradient. In this process, the first-arriving ray path can be traced from the receiver to the source along the negative traveltime gradient, which can be calculated by reconstructing the continuous traveltime field with cubic B-spline interpolation. This new 3-D ray tracing method is compared with the LTI method and the shortest path method (SPM) through a number of numerical experiments. These comparisons show obvious improvements to computed traveltimes and ray paths, both in precision and computational efficiency.     

PARABOLIC AND HYPERBOLIC PARAXIAL TWO-POINT TRAVELTIMES IN 3D MEDIA1

The 4 × 4 T -propagator matrix of a 3D central ray determines, among other important seismic quantities, second-order (parabolic or hyperbolic) two-point traveltime approximations of certain paraxial rays in the vicinity of the known central ray through a 3D medium consisting of inhomogeneous isotropic velocity layers. These rays result from perturbing the start and endpoints of the central ray on smoothly curved anterior and posterior surfaces. The perturbation of each ray endpoint is described only by a two-component vector. Here, we provide parabolic and hyperbolic paraxial two-point traveltime approximations using the T -propagator to feature a number of useful 3D seismic models, putting particular emphasis on expressing the traveltimes for paraxial primary reflected rays in terms of hyperbolic approximations. These are of use in solving several forward and inverse seismic problems. Our results simplify those in which the perturbation of the ray endpoints upon a curved interface is described by a three-component vector. In order to emphasize the importance of the hyperbolic expression, we show that the hyperbolic paraxial-ray traveltime (in terms of four independent variables) is exact for the case of a primary ray reflected from a planar dipping interface below a homogeneous velocity medium.     

Eigenrays in 3D heterogeneous anisotropic media,Part I: Kinematics

We present a new ray bending approach, referred to as the Eigenray method, for solving two‐point boundary‐value kinematic and dynamic ray tracing problems in 3D smooth heterogeneous general anisotropic elastic media. The proposed Eigenray method is aimed to provide reliable stationary ray path solutions and their dynamic characteristics, in cases where conventional initial‐value ray shooting methods, followed by numerical convergence techniques, become challenging. The kinematic ray bending solution corresponds to the vanishing first traveltime variation, leading to a stationary path between two fixed endpoints (Fermat's principle), and is governed by the nonlinear second‐order Euler–Lagrange equation. The solution is based on a finite‐element approach, applying the weak formulation that reduces the Euler–Lagrange second‐order ordinary differential equation to the first‐order weighted‐residual nonlinear algebraic equation set. For the kinematic finite‐element problem, the degrees of freedom are discretized nodal locations and directions along the ray trajectory, where the values between the nodes are accurately and naturally defined with the Hermite polynomial interpolation. The target function to be minimized includes two essential penalty (constraint) terms, related to the distribution of the nodes along the path and to the normalization of the ray direction. We distinguish between two target functions triggered by the two possible types of stationary rays: a minimum traveltime and a saddle‐point solution (due to caustics). The minimization process involves the computation of the global (all‐node) traveltime gradient vector and the traveltime Hessian matrix. The traveltime Hessian is used for the minimization process, analysing the type of the stationary ray, and for computing the geometric spreading of the entire resolved stationary ray path. The latter, however, is not a replacement for the dynamic ray tracing solution, since it does not deliver the geometric spreading for intermediate points along the ray, nor the analysis of caustics. Finally, we demonstrate the efficiency and accuracy of the proposed method along three canonical examples.     

Interval velocity and thickness estimate from wide-angle reflection data

A method to estimate interval velocities and thickness in a horizontal isotropic layered medium from wide-angle reflection traveltime curves is presented. The method is based on a relationship between the squared reflection traveltime differences and the squared offset differences relative to two adjacent reflectors. The envelope of the squared-time versus offset-difference curves, for rays with the same ray parameter, is a straight line, whose slope is the inverse of the square of the interval velocity and whose intercept is the square of the interval time. The method yields velocity and thickness estimates without any knowledge of the overlying stratification. It can be applied to wide-angle reflection data when either information on the upper crust and/or refraction control on the velocity is not available. Application to synthetic and real data shows that the method, used together with other methods, allows us to define a reliable 1D starting model for estimating a depth profile using either ray tracing or another technique.     

速度差对井间走时层析成像质量的影响及其改进方法

在井间地震初至走时层析成像中，随着相邻地质体速度差的增大，使得射线分布不均匀，以及网格剖分不合适，导致层析成像结果不理想.物理和数值模型的井间走时层析成像表明：当速度差超过300%时，层析结果畸变较大；在300%～150%之间时，层析结果较好；低于150%时，层析结果好.在此基础上，提出了井间多尺度初至走时层析成像方法，即对同一模型采用多种网格剖分来同时进行层析成像，以获得研究区域的速度图像.数学和物理模型的井间多尺度走时层析结果表明：该方法很好地兼顾了层析成像的分辨率和精度，极大地改善了井间地震层析成像的质量.即使速度差超过30%，其多尺度的层析结果仍然较好.因此，这种方法具有实际应用的潜力.     

Simultaneous elastic parameter inversion in 2-D/3-D TTI medium combined later arrival times

Traditional traveltime inversion for anisotropic medium is, in general, based on a “weak” assumption in the anisotropic property, which simplifies both the forward part (ray tracing is performed once only) and the inversion part (a linear inversion solver is possible). But for some real applications, a general (both “weak” and “strong”) anisotropic medium should be considered. In such cases, one has to develop a ray tracing algorithm to handle with the general (including “strong”) anisotropic medium and also to design a non-linear inversion solver for later tomography. Meanwhile, it is constructive to investigate how much the tomographic resolution can be improved by introducing the later arrivals. For this motivation, we incorporated our newly developed ray tracing algorithm (multistage irregular shortest-path method) for general anisotropic media with a non-linear inversion solver (a damped minimum norm, constrained least squares problem with a conjugate gradient approach) to formulate a non-linear inversion solver for anisotropic medium. This anisotropic traveltime inversion procedure is able to combine the later (reflected) arrival times. Both 2-D/3-D synthetic inversion experiments and comparison tests show that (1) the proposed anisotropic traveltime inversion scheme is able to recover the high contrast anomalies and (2) it is possible to improve the tomographic resolution by introducing the later (reflected) arrivals, but not as expected in the isotropic medium, because the different velocity (qP, qSV and qSH) sensitivities (or derivatives) respective to the different elastic parameters are not the same but are also dependent on the inclination angle.     

Moveout approximation for vertical seismic profile geometry in a 2D model with anisotropic layers

I introduce a new explicit form of vertical seismic profile (VSP) traveltime approximation for a 2D model with non‐horizontal boundaries and anisotropic layers. The goal of the new approximation is to dramatically decrease the cost of time calculations by reducing the number of calculated rays in a complex multi‐layered anisotropic model for VSP walkaway data with many sources. This traveltime approximation extends the generalized moveout approximation proposed by Fomel and Stovas. The new equation is designed for borehole seismic geometry where the receivers are placed in a well while the sources are on the surface. For this, the time‐offset function is presented as a sum of odd and even functions. Coefficients in this approximation are determined by calculating the traveltime and its first‐ and second‐order derivatives at five specific rays. Once these coefficients are determined, the traveltimes at other rays are calculated by this approximation. Testing this new approximation on a 2D anisotropic model with dipping boundaries shows its very high accuracy for offsets three times the reflector depths. The new approximation can be used for 2D anisotropic models with tilted symmetry axes for practical VSP geometry calculations. The new explicit approximation eliminates the need of massive ray tracing in a complicated velocity model for multi‐source VSP surveys. This method is designed not for NMO correction but for replacing conventional ray tracing for time calculations.     

THE INTERPRETATION OF TRAVELTIME FIELDS IN REFRACTION SEISMOLOGY*

We consider multiply covered traveltimes of first or later arrivals which are gathered along a refraction seismic profile. The two-dimensional distribution of these traveltimes above a coordinate frame generated by the shotpoint axis and the geophone axis or by the common midpoint axis and the offset axis is named a traveltime field. The application of the principle of reciprocity to the traveltime field implies that for each traveltime value with a negative offset there is a corresponding equal value with positive offset. In appendix A procedures are demonstrated which minimize the observational errors of traveltimes inherent in particular traveltime branches or complete common shotpoint sections. The application of the principle of parallelism to an area of the traveltime field associated with a particular refractor can be formulated as a partial differential equation corresponding to the type of the vibrating string. The solution of this equation signifies that the two-dimensional distribution of these traveltimes may be generated by the sum of two one-dimensional functions which depend on the shotpoint coordinate and the geophone coordinate. Physically, these two functions may be interpreted as the mean traveltime branches of the reverse and the normal shot. In appendix B procedures are described which compute these two functions from real traveltime observations by a least-squares fit. The application of these regressed traveltime field data to known time-to-depth conversion methods is straightforward and more accurate and flexible than the use of individual traveltime branches. The wavefront method, the plus-minus method, the generalized reciprocal method and a ray tracing method are considered in detail. A field example demonstrates the adjustment of regressed traveltime fields to observed traveltime data. A time-to-depth conversion is also demonstrated applying a ray tracing method.     

一种新的预条件伴随状态法初至波走时层析

伴随状态法初至波走时层析是基于最优化理论的一种层析成像方法,该方法不必进行射线追踪,用两次正演的计算量便可以获得梯度,具有计算效率高、内存占用小等优点,但是其一阶方向在初始模型或观测孔径不理想的情况下往往无法获得正确的反演结果,而二阶方向的实现又比较困难且费时.在伴随状态法的基础上,将走时差替换为定值,再次进行反演,便可以得到类似于射线密度的矩阵,用该矩阵的逆可以方便地进行预条件.基于该方法,本文提出了一种简单易行的预条件伴随状态法初至波走时层析的实现方法.理论模型和实际资料处理结果都表明,该方法既保留了伴随状态法初至波走时层析的优点,又可以克服一阶方向的局限,获得良好的反演效果.