Improving modelling and inversion in refraction seismics with a first-order Eikonal solver

A first-order Eikonal solver is applied to modelling and inversion in refraction seismics. The method calculates the traveltime of the fastest wave at any point of a regular grid, including head waves as used in refraction. The efficiency, robustness and flexibility of the method give a very powerful modelling tool to find both traveltimes and raypaths. Comparisons with finite-difference data show the validity of the results. Any arbitrarily complex model can be studied, including the exact topography of the surface, thus avoiding static corrections. Later arrivals are also obtained by applying high-slowness masks over the high-velocity zones. Such an efficient modelling tool may be used interactively to invert for the model, but a better method is to apply the refractor-imaging principle of Hagedoorn to obtain the refractors from the picked traveltime curves. The application of this principle has already been tried successfully by previous authors, but they used a less well-adapted Eikonal solver. Some of their traveltimes were not correct in the presence of strong velocity variations, and the refractor-imaging principle was restricted to receiver lines along a plane surface. With the first-order Eikonal solver chosen, any topography of the receiving surface can be considered and there is no restriction on the velocity contrast. Based on synthetic examples, the Hagedoorn principle appears to be robust even in the case of first arrivals associated with waves diving under the refractor. The velocities below the refractor can also be easily estimated, parallel to the imaging process. In this way, the model can be built up successively layer by layer, the refractor-imaging and velocity-mapping processes being performed for each identified refractor at a time. The inverted model could then be used in tomographic inversions because the calculated traveltimes are very close to the observed traveltimes and the raypaths are available.     

Depth-Recursive Tomography of the Bohemian Massif at the CEL09 Transect—Part A: Resolution Estimates and Deblurring Aspects

The refraction CEL09 profile from the CELEBRATION 2000 project intersects the main terranes of the Bohemian Massif in the NW–SE direction: the Saxothuringian, the Teplá-Barrandian, the Moldanubian and the Moravo-Silesian. In its easternmost part, it crosses the Western Outer Carpathians overthrust westward onto the Bohemian Massif. Only the first 450 km were surveyed with the densest deployment of shot points providing data suitable for a reliable geological interpretation. The first-arrival depth-recursive tomography was applied here to derive a P-wave velocity image of the upper and middle crust (Part A). The proper interpretation of the obtained velocity features is the subject of the accompanying paper (Part B). The attained resolution in the velocity image is shown to be superior as compared with the previous CEL09 models based also on the more uncertain later arrivals of reflection waves. The applied DRTG (depth-recursive tomography on grid) method is based on a regular network of refraction grid rays generated for iteratively updated starting models. Only the distinct first arrivals with minimum uncertainty are used for the DRTG inversions to yield the maximum resolution. Thanks to the full control of the data fit by the grid rays used, the statistical lateral resolution could be determined at single grid depths for the chosen confidence levels. Thus, the lateral sizes of the anomalies that can be yet resolved are determined in dependence on their depths and their velocity excesses. The defocusing of the imaged features is studied on the basis of the spatial responses to spike excitation. The calculated spatial responses also allowed the edge smearing of the velocity anomalies to be assessed. Special attention is paid to the imaging of low-velocity zones that are usually suppressed by the smoothing measures used in standard tomographic methods. An improvement was achieved if the smoothing was suggested with regard to the occurrence of the low-velocity zones repeatedly appearing in higher iterations. The gained deblurring effect concerns both the negative and positive anomalies as documented on the velocity features interpreted in the accompanying paper.     

SEISMIC REFRACTION AND SCREENING BY THIN HIGH-VELOCITY LAYERS *

The screening effect of thin, relatively shallow high-velocity layers often presents considerable problems in seismic exploration. Such layers prevent the greater part of the seismic energy from travelling to greater depths and introduce additional refraction arrivals, confusing the seismogram still further. In order to investigate both the screening and refractive properties of high-velocity layers, scale-model experiments have been made over a wide range of layer-thickness/ wavelength ratios (0.05 < d/λ < 2) for suitably chosen material contrasts. The results may be summarised as follows. Refraction arrivals from thin layers in the field may be recognised by their relatively rapid amplitude decay. Furthermore, the “echeloning”-effect observed for refraction first arrivals may be due to the presence of a (thin) layered structure. Since the apparent refraction velocity varies with d/λ when d/λ < 1, differences between vertical well-log velocities and velocities observed along the surface may be expected, making time/depth conversion using surface velocity data inaccurate. Transmission of elastic energy may be expected, if anywhere, only near the shotpoint, at small geophone offset, and for relatively thin screens (d/λ < 0.1). The transmitted signal shape is then independent of the layer thickness. This transmitted energy may be registered either in a reflection set-up with geophones near the shotpoint, or in long-distance refraction work. Three possibilities are offered for overcoming the screening effect of thin high-velocity layers: Use longer-wavelength signals Apply short-spread reflection shooting Apply long-distance refraction shooting The experimental results obtained in scale-model arrangements of such set-ups confirm the potentialities of these methods.     

STATIC CORRECTIONS FOR SHEAR WAVE SECTIONS*

The computation of static corrections requires information about subsurface velocities. This information can be obtained by different methods: surface wave analysis, short refraction lines, downhole times, uphole times and first arrivals from seismograms. For pure shear waves generated by SH sources the analysis of first arrivals from seismograms combined, if necessary, with short refraction lines has proved to be most accurate and economic. A comparison of first-arrival plots from P- and S-wave surveys of the same line measured in areas of unconsolidated sediments in northern Germany illustrates the characteristic differences between the two velocity models. P-waves show a marked velocity increase at the water table from about 600 to 1800 m/s. S-wave velocities of the same strata increase gradually from about 100 to 400 m/s. As a consequence, S-wave models are vertically and laterally more complex and, in general, show no significant velocity increase at a defined boundary as P-wave models do. Therefore, other suitable correction levels with specific velocities must be chosen. A comparison of “tgd-corrections” (correction time between geophone position and datum level) for P- and S-waves in areas of unconsolidated sediments shows that their ratio is different from the P-/S-velocity ratio for the respective correction level because of the greater depth of the S-wave refractor. Therefore, P- and S-waves are influenced by different near-surface anomalies, and time corrections calculated for both wave types are largely independent.     

The problem of velocity inversion in refraction seismics: some observations from modelling results

The applicability of seismic refraction profiling for the detection of velocity inversion, which is also known as a low-velocity layer (LVL), is investigated with the aid of synthetic seismogram computations for a range of models. Our computational models focus on the inherent ambiguities in the interpretation of first-arrival time delays or 'skips' in terms of LVL model parameters. The present modelling results reveal that neither the measure nor even the existence of a shadow zone and/or a time shift (skip) in first arrivals is necessarily indicative of an LVL. Besides attenuation effects, the cap-layer velocity gradient is a critical parameter, determining the termination point of the cap-layer diving wave and thus the time skip.
We suggest that shallow LVLs can be delineated more reliably by traveltime and amplitude modelling of coherent phases reflected from their top and bottom boundaries, often clearly observed in the pre- and near-critical ranges in seismogram sections of refraction profiling experiments with a close receiver spacing. We demonstrate the applicability of this approach for a field data set of a refraction profile in the West Bengal Basin, India. The inferred LVL corresponds to the Gondwana sediments underlying the higher-velocity layer of the Rajmahal Traps. This interpretation is consistent with the data from a nearby well in the region.     

2D inversion of refraction traveltime curves using homogeneous functions

A method using simple inversion of refraction traveltimes for the determination of 2D velocity and interface structure is presented. The method is applicable to data obtained from engineering seismics and from deep seismic investigations. The advantage of simple inversion, as opposed to ray‐tracing methods, is that it enables direct calculation of a 2D velocity distribution, including information about interfaces, thus eliminating the calculation of seismic rays at every step of the iteration process. The inversion method is based on a local approximation of the real velocity cross‐section by homogeneous functions of two coordinates. Homogeneous functions are very useful for the approximation of real geological media. Homogeneous velocity functions can include straight‐line seismic boundaries. The contour lines of homogeneous functions are arbitrary curves that are similar to one another. The traveltime curves recorded at the surface of media with homogeneous velocity functions are also similar to one another. This is true for both refraction and reflection traveltime curves. For two reverse traveltime curves, non‐linear transformations exist which continuously convert the direct traveltime curve to the reverse one and vice versa. This fact has enabled us to develop an automatic procedure for the identification of waves refracted at different seismic boundaries using reverse traveltime curves. Homogeneous functions of two coordinates can describe media where the velocity depends significantly on two coordinates. However, the rays and the traveltime fields corresponding to these velocity functions can be transformed to those for media where the velocity depends on one coordinate. The 2D inverse kinematic problem, i.e. the computation of an approximate homogeneous velocity function using the data from two reverse traveltime curves of the refracted first arrival, is thus resolved. Since the solution algorithm is stable, in the case of complex shooting geometry, the common‐velocity cross‐section can be constructed by applying a local approximation. This method enables the reconstruction of practically any arbitrary velocity function of two coordinates. The computer program, known as godograf , which is based on this theory, is a universal program for the interpretation of any system of refraction traveltime curves for any refraction method for both shallow and deep seismic studies of crust and mantle. Examples using synthetic data demonstrate the accuracy of the algorithm and its sensitivity to realistic noise levels. Inversions of the refraction traveltimes from the Salair ore deposit, the Moscow region and the Kamchatka volcano seismic profiles illustrate the methodology, practical considerations and capability of seismic imaging with the inversion method.     

COMPUTING FIELD STATICS WITH THE HELP OF SEISMIC TOMOGRAPHY*

Field static corrections in general need be applied to all onshore seismic reflection data to eliminate the disturbing effects a weathering layer or near-surface low velocity zone has on the continuity of deep seismic reflections. The traveltimes of waves refracted at the bottom of the low velocity zone (or intermediate refracting interfaces) can often be observed as first breaks on shot records and used to develop a laterally inhomogeneous velocity model for this layer, from which the field static corrections can then be obtained. A simple method is described for computing accurate field statics from first breaks. It is based on a linearization principal for traveltimes and leads to the algorithms that are widely and successfully applied within the framework of seismic tomography. We refine an initial model for the low velocity layer (estimated by a standard traveltime inversion technique) by minimizing the errors between the observed first arrivals on field records and those computed by ray theory through an initial model of the low velocity layer. Thus, one can include more lateral velocity variations within the low velocity layers, which are important to obtain good field static corrections. Traditional first break traveltime inversion methods cannot, in general, provide such refined velocity values. The technique is successfully applied to seismic data from the Amazon Basin. It is based on a simple model for the low velocity layer that consists of an undulating earth surface and one planar horizontal refractor overlain by a laterally changing velocity field.     

最短路径算法下二维层状介质中多次波追踪

在改进后的最短路径算法（MSPM）中引入分区多步计算技术实现了二维层状起伏介质中的多次透射、反射及转换波波前传播的数值模拟,以及相应的走时和射线路径的跟踪计算.其原理是将二维复杂层状模型按速度界面分成若干个独立的计算区域,采用分步计算技术进行多次波的跟踪计算.基于多次波是通过速度界面的简单入射、透射、反射及转换波按一定规律的不同组合,因此可实施分区多步计算技术.通过某一上、下层界面的透射（或透射转换）波实际上是由上层得到的下行波加上由该界面透射的下行波组成,若为转换波则使用不同的速度模型;而经过某一界面的反射（或反射转换）波实际上是由某层内计算得到的下行（或上行）波再加上由该界面反射的上行（或下行）波组成.这样即可得到分区独立计算,并通过速度界面分步连接达到跟踪多次波的目的.计算结果表明MSPM算法下的分区多步计算技术具有单步SPM算法中的诸多优点,即:算法简单、数值计算稳键、计算精度高、速度快及全球解等,因而是解决多次波跟踪计算行之有效的方法.     

Delineation of Subtrappean Mesozoic Sediments in Deccan Syneclise, India, Using Travel-Time Inversion of Seismic Refraction and Wide-Angle Reflection Data

2-D shallow velocity structure is derived by travel-time inversion of the first arrival seismic refraction and wide-angle reflection data along the E–W trending Narayanpur–Nandurbar and N–S Kothar–Sakri profiles, located in the Narmada–Tapti region of the Deccan syneclise. Deccan volcanic (Trap) rocks are exposed along the two profiles. Inversion of seismic data reveals two layered velocity structures above the basement along the two profiles. The first layer with a P-wave velocity of 5.15–5.25 km s^?1 and thickness varying from 0.7–1.5 km represents the Deccan Trap formation along the Narayanpur–Nandurbar profile. The Trap layer velocity ranges from 4.5 to 5.20 km s^?1 and the thickness varies from 0.95 to 1.5 km along the Kothar–Sakri profile. The second layer represents the low velocity Mesozoic sediments with a P-wave velocity of 3.5 km s^?1 and thickness ranging from about 0.70 to 1.6 km and 0.55 to 1.1 km along the E–W and N–S profiles, respectively. Presence of a low-velocity zone (LVZ) below the volcanic rocks in the study area is inferred from the travel-time ‘skip’ and amplitude decay of the first arrival refraction data together with the prominent wide-angle reflection phase immediately after the first arrivals from the Deccan Trap formation. The basement with a P-wave velocity of 5.8–6.05 km s^?1 lies at a depth ranging from 1.5 to 2.45 km along the profiles. The velocity models of the profiles are similar to each other at the intersection point. The results indicate the existence of a Mesozoic basin in the Narmada–Tapti region of the Deccan syneclise.     

Various Approaches to Forward and Inverse Wide-Angle Seismic Modelling Tested on Data from DOBRE-4 Experiment

The interpretation of seismic refraction and wide angle reflection data usually involves the creation of a velocity model based on an inverse or forward modelling of the travel times of crustal and mantle phases using the ray theory approach. The modelling codes differ in terms of model parameterization, data used for modelling, regularization of the result, etc. It is helpful to know the capabilities, advantages and limitations of the code used compared to others.This work compares some popular 2D seismic modelling codes using the dataset collected along the seismic wide-angle profile DOBRE-4, where quite peculiar/uncommon reflected phases were observed in the wavefield.The ~505 km long profile was realized in southern Ukraine in 2009, using 13 shot points and 230 recording stations. Double P_MP phases with a different reduced time (7.5–11 s) and a different apparent velocity, intersecting each other, are observed in the seismic wavefield. This is the most striking feature of the data. They are interpreted as reflections from strongly dipping Moho segments with an opposite dip. Two steps were used for the modelling. In the previous work by Starostenko et al. (2013), the trial-and-error forward model based on refracted and reflected phases (SEIS83 code) was published. The interesting feature is the high-amplitude (8–17 km) variability of the Moho depth in the form of downward and upward bends. This model is compared with results from other seismic inversion methods: the first arrivals tomography package FAST based on first arrivals; the JIVE3D code, which can also use later refracted arrivals and reflections; and the forward and inversion code RAYINVR using both refracted and reflected phases. Modelling with all the codes tested showed substantial variability of the Moho depth along the DOBRE-4 profile. However, SEIS83 and RAYINVR packages seem to give the most coincident results.     

Hagedoorn's plus-minus method: the beauty of simplicity

Hagedoorn developed the plus-minus method with the objective of providing the exploration world with a simple and rapid approximation of Thornburgh's wavefront reconstruction method. Straightforward adding and subtracting produces depths and velocities at all geophones where the refracted waves from a set of reciprocal shots are first arrivals. Even though seismic refraction has not kept up with the revolutionary advance of reflection technology during the past decades, the method still has a lot to offer, especially in shallow engineering, and environmental, groundwater and sea-bed surveying. Its strong points are the ease of application, the low costs, the effectiveness in the shallow zone, the unique ability to provide detailed velocity information on the deepest refractor and the capability for producing parameters for layer and rock characterization.Because of its simplicity, the plus-minus method is ideally suited for real-time processing of refraction data in the field, thus monitoring the data quality and the optimal shot configuration. Today, such field processing is easily performed with a laptop and dedicated software. matlab ^® is a software package that enables tailor-made processing, offering a combination of a programming language, visualization tools and a large library of ready-made functions. This paper presents a plus-minus program developed in matlab and illustrates its application with a case study in Yemen, where seismic refraction was used in a regional groundwater study. Here, refraction provided not only a detailed section across the recharge area of a coastal plain but also the additional information needed to reduce ambiguity in the interpretation of the vertical electrical soundings made in the area.     

THE CORRELATION REFRACTION METHOD AS APPLIED TO WEATHERED ZONE STUDIES IN A GRANITE TERRAIN*

The combined observation of first and later arrivals in shallow seismic refraction surveys, particularly on hard rock terrains, is discussed. Details of experimental weathered-zone investigations by the correlation refraction method in a granite terrain (i.e. field procedure, seismograms obtained, plotting of the data, and identification of the waves are presented). Complete travel time data and interpreted subsurface sections of a few test refraction surveys are included. In one instance the interpreted results of normal and converted refracted wave data have been tested by drilling at three points along a 220 m long profile.     

3D seismic residual statics solutions derived from refraction interferometry

We apply interferometric theory to solve a three‐dimensional seismic residual statics problem to improve reflection imaging. The approach calculates the static solutions without picking the first arrivals from the shot or receiver gathers. The static correction accuracy can be significantly improved by utilising stacked virtual refraction gathers in the calculations. Shots and receivers may be placed at any position in a three‐dimensional seismic land survey. Therefore, it is difficult to determine stationary shots and receivers to form the virtual refraction traces that have identical arrival times, as in a two‐dimensional scenario. To overcome this problem, we use a three‐dimensional super‐virtual interferometry method for residual static calculations. The virtual refraction for a stationary shot/receiver pair is obtained via an integral along the receiver/shot lines, which does not require knowledge of the stationary locations. We pick the maximum energy times on the interferometric stacks and solve a set of linear equations to derive reliable residual static solutions. We further apply the approach to both synthetic and real data.     

APPLICATION OF REFRACTION SEISMICS TO SUBSALT TECTONIC PROBLEMS IN A DEEP SALTDOME BASIN*

Despite the use of CDP and digital methods the Zechstein base is still the deepest horizon in the vast salt-dome basin of Central Europe for which continuous information can be obtained by reflection seismics. Thus in North-western Germany, in addition to reflection seismics, the refraction seismic method has been increasingly used for a reliable survey of deeper horizons. The first part of the paper deals with the investigation of the various possibilities and limitations of refraction seismics with regard to the investigation of Pre-Zechstein layers in a basin with a tectonically very complicated overburden. The recording techniques specially developed for continuous profiling of the desired refraction seismic arrivals and the data processing methods are described. The main problems of interpretation are then discussed, in particular with regard to depth representation. The advantages and disadvantages of the various methods, e.g. Gardner's, Hales' and Wyrobek's, and of the wave-front method, are compared. On account of the tectonically complicated overburden Thornburgh's wave-front method proved to be the most useful. In a further section the various possibilities for velocity determinations are mentioned, e.g. Wyrobek's determination of the overburden velocity, for which the wave-front method automatically furnishes the necessary corrections to a deep datum. Finally, some examples are given for the results obtained, including some incidental information on the deeper crust.     

Determination of Dips and Depths of Near Surface Layers by Radon Transform

The refracted arrivals on seismic shot records have long been recognized as an efficient tool for the computation of detailed near-surface information. In this paper, a new concept of refraction static, which is based on the Radon transform and avoids the tedious process of picking first arrival times, is proposed. This method is particularly suitable when a rough near-surface problem necessitates the utilization of numerous shallow refraction data for the one reflector case. Quasi-linearity of refractors and a constant velocity medium are assumed within the shooting range. Synthetic and real cases have been tested to evaluate the performance of the method. The result is revealed to be satisfactory. Comparison of the synthetic model with the results obtained through the Radon transform reveals a very good accuracy for the proposed method.     

基于交叉梯度约束的地震初至纵波与瑞雷面波联合反演

本文提出了一种初至纵波(P波)与瑞雷面波的交叉梯度联合反演策略.通过对初至P波进行全波形反演可以获得近地表P波速度结构;通过对仅含瑞雷面波信息的地震数据转换到频率-波数域进行加窗振幅波形反演(Windowed-Amplitude Waveform Inversion,w-AWI)可获得近地表横波(S波)速度结构.在二者反演的目标函数中均加入P波速度和S波速度的交叉梯度作为正则化约束项,使得在反演过程中P波速度和S波速度相互制约,相互约束,从而实现对地震初至P波与瑞雷面波的联合反演.数值模拟结果表明交叉梯度联合反演可以提高S波速度反演分辨率,而P波速度反演结果并没有得到提高.实际资料的反演结果表明,交叉梯度联合反演能够获得更加可信的近地表速度结构.     

岩土层速度结构回折波探测原理与方法

岩上层波速有随深度增加的总趋势。用浅层折射法记录到的折射波均属回折波。本文采用垂向速度梯度为不同正的常数多层水平层状介质模型，推导出回析波正反演公式，有简单解析表达式，据此解释的实测剖面，与验证孔波速测井资料相当吻合。表明在岩土层保持天然结构和横向比较均匀地区，回折波探测方法有良好应用前景。     

基于多震相联合反演江苏地区一维速度模型研究

利用地震走时数据,采用联合反演方法获取了江苏地区的一维P波速度模型。与仅采用初至波走时的传统天然地震走时获取方法相比,该方法充分利用了大量存在的续至波参与反演,能有效改进中下地壳的反演能力。针对地震震相目录中常存在震相标识错误的问题,采用的自动判别筛选震相方法能最大限度提高数据走时的精度,可以对不同震相进行有效区分。与其他常用一维速度模型相比,本文反演的模型对Pg、Pn震相走时拟合效果最佳,残差最小。当所用走时数据拥有较高定位精度时,该反演方法能为研究区三维速度结构成像和地震定位提供较可靠的一维速度模型。     

REFRACTION LE LONG DES BANCS MINCES RAPIDES ET EFFET D'ECRAN POUR LES MARQUEURS PROFONDS *

Refraction along thin high velocity layers and along basement is investigated in two cases. a) high velocity layer just on the basement. b) high velocity layer higher above. Period and attenuation of refracted waves are givers as a function of the layer thickness H. Refracted arrivals along thin high velocity layers are visible at significant distances if the layer thickness is not smaller than A/6, where A is the longitudinal wavelength in high velocity medium. The pseudoperiod is proportional to the layer thickness H. The attenuation at large distance follows an x^-ne^-k1x law, where n is close to I and k₁ is inversely proportional to H. Refracted arrivals along the basement are observable even in the case of thin high velocity layers situated in the overburden; their intensity is smaller and their pseudo-period larger than when no layer exists in the overburden. The intensity of the basement arrival decreases and the pseudoperiod increases with increasing layer thickness. The pseudoperiod and the attenuation of refracted arrivals along high velocity layers and along the basement are also highly dependent on acoustic contrasts. Both arrivals from a high velocity layer and from the basement can be recorded simultaneously, provided the frequency spectrum of the seismic chain is sufficiently broad. In all cases layer arrivals show a character very different from basement arrivals.     

Combined structure and velocity stacks via the tau–p transform

Reflection and refraction data are normally processed with tools designed to deal specifically with either near- or far-offset data. Furthermore, the refraction data normally require the picking of traveltimes prior to analysis. Here, an automatic processing algorithm has been developed to analyse wide-angle multichannel streamer data without resorting to manual picking or traveltime tomography. Time–offset gathers are transformed to the tau–p domain and the resulting wavefield is downward continued to the depth–p domain from which a velocity model and stacked section are obtained. The algorithm inputs common-depth-point (CDP) gathers and produces a depth-converted stacked section that includes velocity information. The inclusion of long-offset multichannel streamer data within the tau–p transformation enhances the signal from high-velocity refracted basalt arrivals. Downward continuation of the tau–p transformed wavefield to the depth–p domain allows the reflection and refraction components of the wavefield to be treated simultaneously. The high-slowness depth–p wavefield provides the velocity model and the low-slowness depth–p wavefield may be stacked to give structural information. The method is applied to data from the Faeroe Basin from which sub-basalt velocity images are obtained that correlate with an independently derived P-wave model from the line.