Influence of background heterogeneity on traveltime shifts for compacting reservoirs

Compaction induced by pore‐pressure decrease inside a reservoir can be monitored by measuring traveltime shifts of reflection events on time‐lapse seismic data. Recently we introduced a perturbation‐based formalism to describe traveltime shifts caused by the 3D stress‐induced velocity field around a compacting reservoir. Application of this method to homogeneous background models showed that the offset variation of traveltime shifts is controlled primarily by the anisotropic velocity perturbations and can provide valuable information about the shear and deviatoric stresses. Here, we model and analyse traveltime shifts for compacting reservoirs whose elastic properties are different from those of the surrounding medium. For such models, the excess stress is influenced primarily by the contrast in the rigidity modulus μ across the reservoir boundaries. Synthetic examples demonstrate that a significant (25% or more) contrast in μ enhances the isotropic velocity perturbations outside the reservoir. Nevertheless, the influence of background heterogeneity is mostly confined to the reservoir and its immediate vicinity and the anisotropic velocity changes are still largely responsible for the offset dependence of traveltime shifts. If the reservoir is stiffer than the host rock, the background heterogeneity reduces anisotropic velocity perturbations inside the reservoir but increases them in the overburden. As a result, in this case, the magnitude of the offset variation of traveltime shifts is generally higher for reflections from interfaces above the reservoir. We also study compaction‐induced stress/strain and traveltime shifts for a stiff reservoir embedded in a softer layered model based on velocity profiles from the Valhall Field in the North Sea. Despite producing discontinuities in strain across medium interfaces, horizontal layering does not substantially alter the overall behaviour of traveltime shifts. The most pronounced offset variation of traveltime shifts is observed for overburden events recorded at common midpoints close to the reservoir edges. On the whole, prestack analysis of traveltime shifts should help better constrain compaction‐induced velocity perturbations in the presence of realistic background heterogeneity.     

Influence of subsurface injection on time‐lapse seismic: laboratory studies at seismic and ultrasonic frequencies

Seismic monitoring of reservoir and overburden performance during subsurface CO₂ storage plays a key role in ensuring efficiency and safety. Proper interpretation of monitoring data requires knowledge about the rock physical phenomena occurring in the subsurface formations. This work focuses on rock stiffness and elastic velocity changes of a shale overburden formation caused by both reservoir inflation induced stress changes and leakage of CO₂ into the overburden. In laboratory experiments, Pierre shale I core plugs were loaded along the stress path representative for the in situ stress changes experienced by caprock during reservoir inflation. Tests were carried out in a triaxial compaction cell combining three measurement techniques and permitting for determination of (i) ultrasonic velocities, (ii) quasistatic rock deformations, and (iii) dynamic elastic stiffness at seismic frequencies within a single test, which allowed to quantify effects of seismic dispersion. In addition, fluid substitution effects connected with possible CO₂ leakage into the caprock formation were modelled by the modified anisotropic Gassmann model. Results of this work indicate that (i) stress sensitivity of Pierre shale I is frequency dependent; (ii) reservoir inflation leads to the increase of the overburden Young's modulus and Poisson's ratio; (iii) in situ stress changes mostly affect the P‐wave velocities; (iv) small leakage of the CO₂ into the overburden may lead to the velocity changes, which are comparable with one associated with geomechanical influence; (v) non‐elastic effects increase stress sensitivity of an acoustic waves; (iv) and both geomechanical and fluid substitution effects would create significant time shifts, which should be detectable by time‐lapse seismic.     

Analysis of time‐lapse travel‐time and amplitude changes to assess reservoir compartmentalization

Fluid depletion within a compacting reservoir can lead to significant stress and strain changes and potentially severe geomechanical issues, both inside and outside the reservoir. We extend previous research of time‐lapse seismic interpretation by incorporating synthetic near‐offset and full‐offset common‐midpoint reflection data using anisotropic ray tracing to investigate uncertainties in time‐lapse seismic observations. The time‐lapse seismic simulations use dynamic elasticity models built from hydro‐geomechanical simulation output and a stress‐dependent rock physics model. The reservoir model is a conceptual two‐fault graben reservoir, where we allow the fault fluid‐flow transmissibility to vary from high to low to simulate non‐compartmentalized and compartmentalized reservoirs, respectively. The results indicate time‐lapse seismic amplitude changes and travel‐time shifts can be used to qualitatively identify reservoir compartmentalization. Due to the high repeatability and good quality of the time‐lapse synthetic dataset, the estimated travel‐time shifts and amplitude changes for near‐offset data match the true model subsurface changes with minimal errors. A 1D velocity–strain relation was used to estimate the vertical velocity change for the reservoir bottom interface by applying zero‐offset time shifts from both the near‐offset and full‐offset measurements. For near‐offset data, the estimated P‐wave velocity changes were within 10% of the true value. However, for full‐offset data, time‐lapse attributes are quantitatively reliable using standard time‐lapse seismic methods when an updated velocity model is used rather than the baseline model.     

Traveltime approximation in vertical transversely isotropic layered media

The well‐known asymptotic fractional four‐parameter traveltime approximation and the five‐parameter generalised traveltime approximation in stratified multi‐layer transversely isotropic elastic media with a vertical axis of symmetry have been widely used for pure‐mode and converted waves. The first three parameters of these traveltime expansions are zero‐offset traveltime, normal moveout velocity, and quartic coefficient, ensuring high accuracy of traveltimes at short offsets. The additional parameter within the four‐parameter approximation is an effective horizontal velocity accounting for large offsets, which is important to avoid traveltime divergence at large offsets. The two additional parameters in the above‐mentioned five‐parameter approximation ensure higher accuracy up to a given large finite offset with an exact match at this offset. In this paper, we propose two alternative five‐parameter traveltime approximations, which can be considered extensions of the four‐parameter approximation and an alternative to the five‐parameter approximation previously mentioned. The first three short‐offset parameters are the same as before, but the two additional long‐offset parameters are different and have specific physical meaning. One of them describes the propagation in the high‐velocity layer of the overburden (nearly horizontal propagation in the case of very large offsets), and the other characterises the intercept time corresponding to the critical slowness that includes contributions of the lower velocity layers only. Unlike the above‐mentioned approximations, both of the proposed traveltime approximations converge to the theoretical (asymptotic) linear traveltime at the limit case of very large (“infinite”) offsets. Their accuracy for moderate to very large offsets, for quasi‐compressional waves, converted waves, and shear waves polarised in the horizontal plane, is extremely high in cases where the overburden model contains at least one layer with a dominant higher velocity compared with the other layers. We consider the implementation of the proposed traveltime approximations in all classes of problems in which the above‐mentioned approximations are used, such as reflection and diffraction analysis and imaging.     

Anisotropic 4D seismic response inferred from ultrasonic laboratory measurements: A direct comparison with the isotropic response

Pore-pressure depletion causes changes in the triaxial stress state. Pore-pressure depletion in a flat reservoir, for example, can be reasonably approximated as uniaxial compaction, in which the horizontal effective stress change is smaller than the vertical effective stress. Furthermore, the stress sensitivity of velocities can be angle-dependent. Therefore, time-lapse changes in reservoir elastic anisotropy are expected as a consequence of production, which can complicate the interpretation of the 4D seismic response. The anisotropic 4D seismic response caused by pore-pressure depletion was investigated using existing core velocity measurements. To make a direct comparison between the anisotropic 4D seismic response and the isotropic response based only on vertical velocities, pseudoisotropic elastic properties were utilized, and the two responses were compared in terms of a dynamic rock physics template. A comparison of the dynamic rock physics templates indicates that time-lapse changes in reservoir elastic anisotropy have a noticeable impact on the interpretation of 4D seismic data. Changes in anisotropy as a result of pore-pressure depletion cause a time-lapse amplitude variation with offset response as if there is a reduction in V_P/V_S (i.e., pseudoisotropic V_P/V_S decreases), although the vertical V_P/V_S increases. The impact of time-lapse changes in anisotropy on the amplitude variation with offset gradient was also investigated, and the time-lapse anisotropy was found to enhance changes in the amplitude variation with offset gradient for a given case.     

Review Paper: A semi‐empirical model of strain sensitivity for 4D seismic interpretation

We reformulate the original model of Hatchell and Bourne and Røste, Stovas and Landrø that couples fractional velocity change to subsurface strain via a fundamental constant R. The new model combines elastic compressibility of a dual‐porosity system for a sand–shale mixture with horizontal planes of inter‐granular weakness. The majority of observed R‐factor magnitudes from post‐stack 4D seismic data in both the reservoir and overburden can thus be explained. R is predicted to depend strongly on lithology and also initial strain state. The model is also extended to predict the observed angle‐dependence of time‐lapse time‐shifts from pre‐stack data. An expression for the gradient of time‐shift with incidence angle is obtained in terms of the background V_P/V_S, and also the ratio of tangential to normal compliances B_T/B_N representing loss or creation of inter‐granular coupling. If accurately estimated from data, this compliance ratio can be used as an additional parameter to assess the post‐production state of the overburden. It is concluded that whilst R remains the over‐arching parameter controlling the magnitude of time‐shifts measured from 4D seismic data, B_T/B_N is a subtler parameter that may also prove of future value.     

Velocity/interface model building in a thrust belt by tomographic inversion of global offset seismic data

Between September and November 1999, two test seismic lines were recorded in the southern Apennine region of southern Italy using the global offset technique, which involves the acquisition of a wide offset range using two simultaneously active seismic spreads. One consisted of a symmetrical spread moving along the line, with geophone arrays every 30 m and a maximum offset of 3.6 km. The other one consisted of fixed geophone arrays every 90 m with a maximum offset of 18 km. This experimental acquisition project was carried out as part of the enhanced seismic in thrust belt (ESIT) research project, funded by the European Union, Enterprise Oil and Eni-Agip. An iterative and interactive tomographic inversion of refraction/reflection arrivals was carried out on the data from line ESIT700 to produce a velocity/interface model in depth, which used all the available offsets. The tomographic models allowed the reconstruction of layer interface geometries and interval velocities for the target carbonate platform (Apula) and the overburden sequence. The value of this technique is highlighted by the fact that the standard approach, based on near-vertical reflection seismic and a conventional processing flow, produced poor seismic images in both stack and migrated sections.     

Review Paper: Post-stack 4D seismic time-shifts: Interpretation and evaluation

A review and analysis of post-stack time-lapse time-shifts has been carried out that covers published literature supplemented by in-house datasets available to the authors. Time-shift data are classified into those originating from geomechanical effects and those due to fluid saturation changes. From these data, conclusions are drawn regarding the effectiveness of post-stack time-shifts for overburden and reservoir monitoring purposes. A variety of field examples are shown that display the range and magnitude of variation for each class of application. The underlying physical mechanisms creating these time-shifts are then described, and linked to a series of generic and field-specific rock physics calculations that predict their magnitudes. These calculations serve as a guide for practitioners wishing to utilize this information on their own datasets. Conclusions are drawn regarding the reliability of this attribute for monitoring purposes, and the extent to which further development is required and how it should be reported by authors.     

Sensitivity of time-lapse seismic to reservoir stress path

The change in reservoir pore pressure due to the production of hydrocarbons leads to anisotropic changes in the stress field acting on the reservoir. Reservoir stress path is defined as the ratio of the change in effective horizontal stress to the change in effective vertical stress from the initial reservoir conditions, and strongly influences the depletion‐induced compaction behaviour of the reservoir. Seismic velocities in sandstones vary with stress due to the presence of stress‐sensitive regions within the rock, such as grain boundaries, microcracks, fractures, etc. Since the response of any microcracks and grain boundaries to a change in stress depends on their orientation relative to the principal stress axes, elastic‐wave velocities are sensitive to reservoir stress path. The vertical P‐ and S‐wave velocities, the small‐offset P‐ and SV‐wave normal‐moveout (NMO) velocities, and the P‐wave amplitude‐versus‐offset (AVO) are sensitive to different combinations of vertical and horizontal stress. The relationships between these quantities and the change in stress can be calibrated using a repeat seismic, sonic log, checkshot or vertical seismic profile (VSP) at the location of a well at which the change in reservoir pressure has been measured. Alternatively, the variation of velocity with azimuth and distance from the borehole, obtained by dipole radial profiling, can be used. Having calibrated these relationships, the theory allows the reservoir stress path to be monitored using time‐lapse seismic by combining changes in the vertical P‐wave impedance, changes in the P‐wave NMO and AVO behaviour, and changes in the S‐wave impedance.     

Optimization of spread configuration for seismic data acquisition through numerical modeling in tectonically complex areas: a case study from Badarpur anticline, Cachar, India

Cachar area, a part of the Assam–Arakan frontal fold belt in the eastern part of India is characterized by a series of narrow anticlines which are dissected by a number of faults including thrusts, and broad synclines. Out of several exposed structures, Badarpur anticline is most important because of its major share of hydrocarbon strikes in Cachar. However, the seismic image of this area so far has been extremely poor for deciding locations for exploratory drilling. The problems in obtaining a good seismic image in this area include proper spread configuration including far and near offsets to capture reflected signals from the complex subsurface, variation of charge-sizes and depths along a line, large variation in statics due to rapid changes in elevation and near surface properties. The presence of coherent noise due to scattering and proper migration are major processing problems. This paper deals with the problem of spread configuration. For areas of complex subsurface, spread configuration obtained through field experiment at one spot may be quite different from other spots, requiring a large number of costly experiments. In this work, therefore, we used numerical modeling to derive an optimum spread configuration. For this, we chose a line passing through four wells on Badarpur anticline, and prepared a depth model using available seismic, surface geological and well data. Continuous synthetic common shot records were generated through oblique incidence ray tracing and were processed to obtain stacks with various spread configurations. Comparison of these stacks of synthetic data shows that end-on shooting configuration is inadequate for acquiring a good seismic image of the complex subsurface. From the present study it is evident that a split–spread configuration with near and far offsets of 50 m and 2400 m, respectively, can bring out optimum seismic images in Badarpur and other anticlines in this area.     

Lithology-dependent seismic anisotropic amplitude variation with offset correction in transversely isotropic media

Analysis of amplitude variation with offset is an essential step for reservoir characterization. For an accurate reservoir characterization, the amplitude obtained with an isotropic assumption of the reservoir must be corrected for the anisotropic effects. The objective is seismic anisotropic amplitude correction in an effective medium, and, to this end, values and signs of anisotropic parameter differences (Δδ and Δε) across the reflection interfaces are needed. These parameters can be identified by seismic and well log data. A new technique for anisotropic amplitude correction was developed to modify amplitude changes in seismic data in transversely isotropic media with a vertical axis of symmetry. The results show that characteristics of pre-stack seismic data, that is, amplitude variation with offset gradient, can be potentially related to the sign of anisotropic parameter differences (Δδ and Δε) between two layers of the reflection boundary. The proposed methodology is designed to attain a proper fit between modelled and observed amplitude variation with offset responses, after anisotropic correction, for all possible lithofacies at the reservoir boundary. We first estimate anisotropic parameters, that is, δ and ε, away from the wells through Backus averaging of elastic properties resulted from the first pass of isotropic pre-stack seismic inversion, on input data with no amplitude correction. Next, we estimate the anisotropic parameter differences at reflection interfaces (values and signs of Δδ and Δε). We then generate seismic angle gather data after anisotropic amplitude correction using Rüger's equation for the P-P reflection coefficient. The second pass of isotropic pre-stack seismic inversion is then performed on the amplitude-corrected data, and elastic properties are estimated. Final outcome demonstrates how introduced methodology helps to reduce the uncertainty of elastic property prediction. Pre-stack seismic inversion on amplitude-corrected seismic data results in more accurate elastic property prediction than what can be obtained from non-corrected data. Moreover, a new anisotropy attribute (ν) is presented for improvement of lithology identification.     

Moveout approximation for horizontal transversely isotropic and vertical transversely isotropic layered medium. Part I: 1D ray propagation‡

Anisotropy in subsurface geological models is primarily caused by two factors: sedimentation in shale/sand layers and fractures. The sedimentation factor is mainly modelled by vertical transverse isotropy (VTI), whereas the fractures are modelled by a horizontal transversely isotropic medium (HTI). In this paper we study hyperbolic and non‐hyperbolic normal reflection moveout for a package of HTI/VTI layers, considering arbitrary azimuthal orientation of the symmetry axis at each HTI layer. We consider a local 1D medium, whose properties change vertically, with flat interfaces between the layers. In this case, the horizontal slowness is preserved; thus, the azimuth of the phase velocity is the same for all layers of the package. In general, however, the azimuth of the ray velocity differs from the azimuth of the phase velocity. The ray azimuth depends on the layer properties and may be different for each layer. In this case, the use of the Dix equation requires projection of the moveout velocity of each layer on the phase plane. We derive an accurate equation for hyperbolic and high‐order terms of the normal moveout, relating the traveltime to the surface offset, or alternatively, to the subsurface reflection angle. We relate the azimuth of the surface offset to its magnitude (or to the reflection angle), considering short and long offsets. We compare the derived approximations with analytical ray tracing.     

非零偏VSP弹性波叠前逆时深度偏移技术探讨

非零偏VSP地震资料是一种多分量资料,处理非零偏VSP资料,弹性波叠前逆时深度偏移技术无疑是最适合的处理技术.本文从二维各向同性介质的弹性波波动方程出发,研究了对非零偏VSP资料进行叠前逆时深度偏移的偏移算法,讨论了逆时传播过程中的边值问题和数值频散问题及其相应的解决方案;采用求解程函方程计算得到地下各点的地震波初至时间作为成像时间,实现了非零偏VSP资料的叠前逆时深度偏移.最后进行了模型试算和非零偏VSP地震资料的试处理,结果表明该方法不受地层倾角限制,较适用于高陡构造地区或介质横向速度变化较大地区的非零偏VSP地震资料处理.     

转换波四参数速度分析方法在k71地区的应用

3-D converted-wave data were acquired using digital MEMS （micro-electromechanical system） three component （3C） sensors in the alternating sand and shale sequence in the overburden of the Shengli Ken-71 area. This gives rise to serious non-hyperbolic moveout effects in the converted-wave data due to both the asymmetrical ray path and anisotropic effects. Conventional velocity analysis and moveout correction based on isotropic methods do not flatten reflections events. Here, we use a four-parameter theory to evaluate these effects and process the data. These four parameters include the PS converted wave stacking velocity （Vc2）, the vertical velocity ratio （Y0）, the effective velocity ratio （Yeff）, and the anisotropy parameter （xoff）, The method utilizes the moveout information at different offsets to estimate the different parameters and ensures that the events are properly aligned for stacking, As a result, this four-parameter theory leads to an improvement in imaging quality and correlation between the P-waves and converted-waves.     

REFLECTION-WINDOW MAPPING OF SHALLOW BEDROCK*

A field study was undertaken to evaluate the effectiveness of the high frequency seismic reflection technique for mapping of shallow and irregular bedrock. Bedrock reflections were obtained using a hammer source with both in-line and common offset field layouts. The recording equipment included 12-channel enhancement seismographs, 28 Hz vertical geo-phones and a microcomputer. The latter increased the overall versatility of the seismic system. Field sites for this study are typical of the geological settings of the tin mining areas of Malaysia. The topographical ‘lows’ of the irregular bedrock control the localization of tin ore. The subsurface geology consists of a thin low velocity layer (± 300 m/s) overlying the compact overburden (± 1700 m/s) which in turn lies on bedrock. This paper discusses various criteria for designing an optimum window for obtaining usable reflections between the first arrival and the leading edge of the ground roll cone. Detailed mapping of the overburden and the bedrock interface by the reflection method can be useful in delineating areas for exploratory drilling and for optimum planning of mining operations.     

Enabling affordable omnidirectional subsurface extended image volumes via probing

Image gathers as a function of subsurface offset are an important tool for the inference of rock properties and velocity analysis in areas of complex geology. Traditionally, these gathers are thought of as multidimensional correlations of the source and receiver wavefields. The bottleneck in computing these gathers lies in the fact that one needs to store, compute, and correlate these wavefields for all shots in order to obtain the desired image gathers. Therefore, the image gathers are typically only computed for a limited number of subsurface points and for a limited range of subsurface offsets, which may cause problems in complex geological areas with large geologic dips. We overcome increasing computational and storage costs of extended image volumes by introducing a formulation that avoids explicit storage and removes the customary and expensive loop over shots found in conventional extended imaging. As a result, we end up with a matrix–vector formulation from which different image gathers can be formed and with which amplitude‐versus‐angle and wave‐equation migration velocity analysis can be performed without requiring prior information on the geologic dips. Aside from demonstrating the formation of two‐way extended image gathers for different purposes and at greatly reduced costs, we also present a new approach to conduct automatic wave‐equation‐based migration‐velocity analysis. Instead of focusing in particular offset directions and preselected subsets of subsurface points, our method focuses every subsurface point for all subsurface offset directions using a randomized probing technique. As a consequence, we obtain good velocity models at low cost for complex models without the need to provide information on the geologic dips.     

An estimation method for effective stress changes in a reservoir from 4D seismics data

Advances in seismics acquisition and processing and the widespread use of 4D seismics have made available reliable production‐induced subsurface deformation data in the form of overburden time‐shifts. Inversion of these data is now beginning to be used as an aid to the monitoring of a reservoir's effective stress. Past solutions to this inversion problem have relied upon analytic calculations for an unrealistically simplified subsurface, which can lead to uncertainties. To enhance the accuracy of this approach, a method based on transfer functions is proposed in which the function itself is calibrated using numerically generated overburden strain deformation calculated for a small select group of reference sources. This technique proves to be a good compromise between the faster but more accurate history match of the overburden strain using a geomechanical simulator and the slower, less accurate analytic method. Synthetic tests using a coupled geomechanical and fluid flow simulator for the South Arne field confirm the efficacy of the method. Application to measured time‐shifts from observed 4D seismics indicates compartmentalization in the Tor reservoir, more heterogeneity than is currently considered in the simulation model and moderate connectivity with the overlying Ekofisk formation.     

Wide-angle linear forward modelling of synthetic seismograms

We address the issue of linearity and scale dependence in forward modelling of seismic data from well logs, for large ray parameters, wide angles or large offsets. We present a forward model, within the context of seismic‐to‐well matching, that is linearized in the elastic properties of the earth. This model preserves linearity at large ray parameters and can handle fine‐layering effects such as induced anisotropy. Starting from a low‐contrast small‐ray‐parameter model, we extend it to a large‐ray‐parameter model by fully linearizing the elastic‐property contrasts. Overall linearity of the forward model is extended by partitioning the compressional‐wave and shear‐wave velocity fields into two fundamental scales: a kinematic scale that governs wavefield propagation effects and a dynamic scale that governs wavefield scattering effects. This analysis reveals that the standard practice in forward modelling of strongly filtering the ratio of compressional‐wave velocity to shear‐wave velocity is well founded in the underlying physics. The partitioning of the velocity fields also leads naturally to forward modelling that accounts fully for stretch effects, to resolution of the angle‐of‐incidence versus ray‐parameter dichotomy in seismic‐amplitude analysis, and to full accounting for induced anisotropy and dispersion effects due to fine‐layering of isotropic media. With the onset of routine long‐offset acquisition and the compelling need to optimize asset management in order to maximize reserve recovery, this forward model recognizes the physics of seismic wave propagation and enables a more complete exploitation of amplitude information in pre‐critical seismic data.     

地球介质自组织性对地震波走时和振幅的影响

本文采用随机介质描述地球内部大尺度背景场上存在的小尺度不均匀性和自组织结构;文中没有采用传统的层状网格结构,提出了块状结构来描述复杂的二维自组织结构的方法,分别以高斯型、指数型和von Karman型等自相关函数描述各向同性和各向异性非均匀分布自组织特征;采用射线追踪分析了不同分布特征自组织结构对地震波运动学和动力学特征参数的影响;结果表明,由于地球内部介质的自组织性存在,地震波射线轨迹可能发生明显的畸变; 不同偏移距处,地震反射振幅减弱或增强; 自组织结构从高斯型到von Karman型,在小尺度上表现更大的非均匀性,因此走时和振幅表征依次更强的平均效应.     

A fast and robust two-point ray tracing method in layered media with constant or linearly varying layer velocity

A fast and robust method for two-point ray tracing in one-dimensional layered media is presented. This method is applicable to layered models with constant or linearly varying isotropic layer velocity. For given model properties and source and receiver positions, a ray path can be uniquely determined once its ray parameter (i.e. horizontal slowness) is known. The ray parameter can be obtained by numerically solving the nonlinear offset (i.e. source–receiver horizontal distance) equation using Newton's method, which generally works well at near and mid offsets. However, Newton's method becomes hard to converge at large offsets due to the oversensitivity of offset to ray parameter. Based on the analysis of the characteristic of the offset equation, a modified ray parameter is proposed and used to replace the generic ray parameter in numerical calculation. Numerical experiments show that the iteration process becomes stable and converges rapidly with the modified ray parameter. Moreover, a rational function that asymptotically approximates the shape of the offset equation is introduced for obtaining good initial estimates of the modified ray parameter. Numerical tests show that this method is robust in any situation, and an accurate ray parameter can be obtained within two or three iterations for a wide range of model velocity structure and source–receiver distance. Furthermore, the proposed two-point ray tracing method is easy to implement.