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.     

Interpretation of shallow refraction seismic data by reflection/refraction tomography

A new algorithm for tomographic inversion of traveltimes of reflected and refracted seismic waves is developed. The inversion gives interface configurations and velocity distributions in layers. The important features of the algorithm are: (a) the inclusion of shot time delays in the list of unknown parameters; (b) the regularization is applied in such a way that the most probable model is characterized by the similarity of neighbouring interfaces. As the problem under consideration is non-linear, several iterations are necessary in order to obtain the final model. In the case of a very inexact initial model, a 'layer-by-layer' inversion strategy is recommended as a first inversion step. The inversion program is supplied with a user interface, thanks to which one can: (a) pick interactively and identify seismic traveltimes; (b) build and edit depth/velocity models; and (c) display calculated traveltime curves and compare them with picked traveltimes as well as with the original seismic sections. The efficiency of the inversion software developed is illustrated by a numerical example and a field example in which shallow seismic data are considered. Application to wide-aperture reflection/refraction profiling (WARRP) data is also possible.     

NON-NORMAL INCIDENCE INVERSION: EXISTENCE OF SOLUTION*

An inverse problem is one in which the parameters of a model are determined from measured seismic data. Important to the solution of inverse problems is the issue of whether or not a solution exists. In this paper we show, in a constructive manner, that a solution does exist to the specific inverse problem of determining the parameters of a horizontally stratified, lossless, isotropic and homogeneous layered system that is excited by a non-normal incidence (NNI) plane wave. Mode conversion between P- and S-waves is included. We develop a seven-step layer-recursive procedure for determining all of the parameters for layer j. These parameters are P-wave and S-wave velocities and angles of incidence, density, thickness, traveltimes, and reflection- and transmission-coefficient matrices. Downward continuation of data from the top of one layer to the top of the next lower layer is an important step in our procedure, just as it is in normal incidence (NI) inversion. We show that, in order to compute all parameters of layer j, we need to (and can) compute some parameters for layer j+ 1. This is a non-causal phenomenon that seems to be necessary in NNI inversion but is not present in NI inversion.     

METHOD FOR DETERMINATION OF VELOCITY AND DEPTH FROM SEISMIC REFLECTION DATA1

The estimation of velocity and depth is an important stage in seismic data processing and interpretation. We present a method for velocity-depth model estimation from unstacked data. This method is formulated as an iterative algorithm producing a model which maximizes some measure of coherency computed along traveltimes generated by tracing rays through the model. In the model the interfaces are represented as cubic splines and it is assumed that the velocity in each layer is constant. The inversion includes the determination of the velocities in all the layers and the location of the spline knots. The process input consists of unstacked seismic data and an initial velocity-depth model. This model is often based on nearby well information and an interpretation of the stacked section. Inversion is performed iteratively layer after layer; during each iteration synthetic travel-time curves are calculated for the interface under consideration. A functional characterizing the main correlation properties of the wavefield is then formed along the synthetic arrival times. It is assumed that the functional reaches a maximum value when the synthetic arrival time curves match the arrival times of the events on the field gathers. The maximum value of the functional is obtained by an effective algorithm of non-linear programming. The present inversion algorithm has the advantages that event picking on the unstacked data is not required and is not based on curve fitting of hyperbolic approximations of the arrival times. The method has been successfully applied to both synthetic and field data.     

Estimating the elastic parameters of anisotropic media using a joint inversion of P-wave and SV-wave traveltime error

A method is presented to estimate the elastic parameters and thickness of media that are locally laterally homogeneous using P‐wave and vertically polarized shear‐wave (SV‐wave) data. This method is a ‘layer‐stripping’ technique, and it uses many aspects of common focal point (CFP) technology. For each layer, a focusing operator is computed using a model of the elastic parameters with which a CFP gather can be constructed using the seismic data. Assuming local homogeneity, the resulting differential time shifts (DTSs) represent error in the model due to anisotropy and error in thickness. In the (τ?p) domain, DTSs are traveltimes Δτ that connect error in layer thickness z, vertical slowness q, and ray parameter p. Series expansion is used to linearize Δτ with respect to error in the elastic parameters and thickness, and least‐squares inversion is used to update the model. For stability, joint inversion of P and SV data is employed and, as pure SV data are relatively rare, the use of mode‐converted (PSV) data to represent SV in the joint inversion is proposed. Analytic and synthetic examples are used to demonstrate the utility and practicality of this inversion.     

Determination of a shallow velocity–depth model from seismic refraction data by coherence inversion1

Seismic refractions have different applications in seismic prospecting. The travel- times of refracted waves can be observed as first breaks on shot records and used for field static calculation. A new method for constructing a near-surface model from refraction events is described. It does not require event picking on prestack records and is not based on any approximation of arrival times. It consists of the maximization of the semblance coherence measure computed using shot gathers in a time window along refraction traveltimes. Time curves are generated by ray tracing through the model. The initial model for the inversion was constructed by the intercept-time method. Apparent velocities and intercept times were taken from a refraction stacked section. Such a section can be obtained by appling linea moveout corrections to common-shot records. The technique is tested successfully on synthetic and real data. An important application of the proposed method for solving the statics problem is demonstrated.     

Exact traveltime computation in multi-layered transversely isotropic media with vertical symmetry axis

An approach to calculate the accurate ray paths and traveltimes in multi-layered VTI media (transversely isotropic media with a vertical symmetry axis) is proposed. The expressions of phase velocity, group velocity and Snell’s law used for computation are all explicit and exact. The calculation of ray paths and traveltimes for a given ele-mentary wave is equivalent to that of a transmission problem which is much easier to be treated with the formulae proposed. In the section of numerical examples, the proce...     

Near-surface layer traveltime inversion: a synthetic example1

The near-surface layer is modelled as a constant-velocity layer with varying thickness. The base of the layer is described by a B-spline curve. The optimum model is calculated by minimizing, with respect to the model parameters, the difference between traveltimes predicted by the model and those observed in the data. Once a model has been produced, corrections that are dependent on the raypath geometry through the near-surface layer can be calculated. The effect of the near-surface layer is normally considered to be consistent at each shot or geophone station for all traveltimes arriving at that location (the surface-consistent approximation). This assumption linearizes the problem, allowing timeshifts to be calculated and the traveltimes corrected to a chosen datum, representing static corrections. The single correction at each point is an averaged correction, based on an assumption that is particularly inaccurate in the presence of lateral variations of velocity or thickness of the surface layer, in the presence of large surface layer velocities or in the presence of a thick surface layer. The method presented considers the non-linear relationship between data and model explicitly, hence the correction that is dependent on the raypath. Linearization removes this dependence and reduces the problem to a surface-consistent approximation. The method is applied to synthetic data calculated from a model with surface layer variations. Comparisons are made between the corrected data resulting from the method described here and the conventional surface-consistent approach. From these results it becomes apparent that the near-surface layer inversion method presented here can reproduce accurate models and correct for near-surface layer effects in cases where conventional methods encounter difficulties. Additionally the method can be readily extended to 3D.     

Common-reflection-point time migration

Since the early days of seismic processing, time migration has proven to be a valuable tool for a number of imaging purposes. Main motivations for its widespread use include robustness with respect to velocity errors, as well as fast turnaround and low computation costs. In areas of complex geology, in which it has well-known limitations, time migration can still be of value by providing first images and also attributes, which can be of much help in further, more comprehensive depth migration. Time migration is a very close process to common-midpoint (CMP) stacking and, more recently, to zero-offset commonreflection- surface (CRS) stacking. In fact, Kirchhoff time migration operators can be readily formulated in terms of CRS parameters. In the nineties, several studies have shown advantages in the use of common-reflection-point (CRP) traveltimes to replace conventional CMP traveltimes for a number of stacking and migration purposes. In this paper, we follow that trend and introduce a Kirchhoff-type prestack time migration and velocity analysis algorithm, referred to as CRP time migration. The algorithm is based on a CRP operator together with optimal apertures, both computed with the help of CRS parameters. A field-data example indicates the potential of the proposed technique.     

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.     

各向异性TI介质qP反射波走时层析成像

地震走时层析成像是反演地层各向异性参数分布的有效方法,但是关于地震各向异性介质走时层析成像的研究并不多,其技术远远没有达到成熟的阶段.在野外数据采集时,地表反射波观测方式相对井间和垂直地震剖面观测方式的成本更低,利用qP反射波走时反演各向异性参数具有更加广泛的实用价值.本文实现的TI介质地震走时层析成像方法结合了TI介质反射波射线追踪算法、走时扰动方程和非线性共轭梯度算法,它可以对任意强度的TI介质模型进行反演,文中尝试利用qP反射波走时重建TI介质模型的参数图像.利用qP反射波对层状介质模型和块状异常体模型进行走时反演,由于qP波相速度对弹性模量参数和Thomsen参数的偏微分不同,所以可以分别反演弹性模量参数和Thomsen参数.数值模拟结果表明:利用qP反射波可以反演出TI介质模型的弹性模量参数与Thomsen参数,不同模型的走时迭代反演达到了较好的收敛效果,与各向同性介质走时反演结果相比较,各向异性介质走时反演结果具有较好的识别能力.     

工程地震钻孔柱状图的微机绘图系统

根据工程地震工作的需要,建立了适用于工程地震的岩土柱状图形库,编制了一套工程地震钻孔柱状图的微机绘图系统.该系统可直接生成集各岩土层参数、标贯值、岩土柱状图及剪切波速直方图等内容于一体的工程地震钻孔柱状图.所形成的图件规范、美观,内容完整、精确,大大提高了工作效率.     

ESTIMATION OF ELASTIC PARAMETERS FROM AVO EFFECTS IN THE TAU-P DOMAIN1

A quantitative AVO algorithm suitable for media with slow lateral parameter variations is developed. The method is based on a target-oriented inversion scheme for estimation of elastic parameters in a locally horizontally stratified medium. The algorithm uses band-limited PP reflection coefficients in the τ-p domain to estimate P- and S-wave velocities, densities and layer thicknesses. To obtain these reflection coefficients, a pre-processing involving the Radon transform and multiple attenuation is necessary. Furthermore, a macromodel for the velocities above the target zone must be found prior to the inversion. Various inversion tests involving synthetic data with white Gaussian noise and modelling errors that are likely to occur in conjunction with real data have been performed. In general, the inversion algorithm is fairly robust, since it is able to reproduce the main features of the reference model: main interface locations and relative contrasts in the three unknown layer parameters are recovered. From a test combining the effect of source directivity, one thin layer and 20% white Gaussian noise, it was found that neglect of the source directivity in the inversion caused the largest errors in the estimates. This indicates that it is very important either to eliminate the source directivity in a preprocessing step, or to take the directivity into account in the present algorithm. Despite these problems it was concluded that the inversion algorithm was able to reproduce the main features of the reference model.     

利用联合反演技术进行反射地震的波速成象

本文介绍了根据反射地震数据进行波速成象的一种方法,其基础为多种反演技术的综合。由于要求的波速图象C（x,z）具有间断性,除利用走时数据T（x,t）外,在地层比较水平的情况下,还利用了均方根速度V（x,t）和统计子波W（t）的数据来成象。计算机层析成象过程分为三步:首先重做速度分析,取得与初次反射走时一致的均方根速度数据;然后用反射走时与均方根速度联合反演对应分析道的层速度和界面深度;最后由联合反演结果和反射面走时求波速图象函数的数字化版。文中还给出了波速成象方法在我国西北某沉积盆地上的应用及验证结果。     

P‐wave stacking‐velocity tomography for VTI media

A major complication caused by anisotropy in velocity analysis and imaging is the uncertainty in estimating the vertical velocity and depth scale of the model from surface data. For laterally homogeneous VTI (transversely isotropic with a vertical symmetry axis) media above the target reflector, P‐wave moveout has to be combined with other information (e.g. borehole data or converted waves) to build velocity models for depth imaging. The presence of lateral heterogeneity in the overburden creates the dependence of P‐wave reflection data on all three relevant parameters (the vertical velocity V_P0 and the Thomsen coefficients ε and δ) and, therefore, may help to determine the depth scale of the velocity field. Here, we propose a tomographic algorithm designed to invert NMO ellipses (obtained from azimuthally varying stacking velocities) and zero‐offset traveltimes of P‐waves for the parameters of homogeneous VTI layers separated by either plane dipping or curved interfaces. For plane non‐intersecting layer boundaries, the interval parameters cannot be recovered from P‐wave moveout in a unique way. Nonetheless, if the reflectors have sufficiently different azimuths, a priori knowledge of any single interval parameter makes it possible to reconstruct the whole model in depth. For example, the parameter estimation becomes unique if the subsurface layer is known to be isotropic. In the case of 2D inversion on the dip line of co‐orientated reflectors, it is necessary to specify one parameter (e.g. the vertical velocity) per layer. Despite the higher complexity of models with curved interfaces, the increased angle coverage of reflected rays helps to resolve the trade‐offs between the medium parameters. Singular value decomposition (SVD) shows that in the presence of sufficient interface curvature all parameters needed for anisotropic depth processing can be obtained solely from conventional‐spread P‐wave moveout. By performing tests on noise‐contaminated data we demonstrate that the tomographic inversion procedure reconstructs both the interfaces and the VTI parameters with high accuracy. Both SVD analysis and moveout inversion are implemented using an efficient modelling technique based on the theory of NMO‐velocity surfaces generalized for wave propagation through curved interfaces.     

基于HAFMM的无射线追踪跨孔雷达走时层析成像

本文使用最小二乘线性迭代反演方法对跨孔雷达直达波初至时数据进行反演,每次迭代过程中,用有限差分法求解走时程函方程,并用高精度快速推进方法(HAFMM)进行波前扩展,通过追踪波前避免了进行射线追踪.为了验证该方案,我们对三组合成数据进行了测试,分析了单位矩阵算子、一阶差分算子和拉普拉斯算子等三种不同模型参数加权算子对模型的约束和平滑效果;讨论了FMM和HAFMM对反演精度的影响;测试了LSQR,GMRES和BICGSTAB等三种矩阵反演算法的反演效果.此外,我们还对一组野外实测数据进行了反演,对比了基于本方案以及基于平直射线追踪和弯曲射线追踪的走时层析成像反演效果.对比分析结果表明,使用拉普拉斯算子和HAFMM进行反演能较好地进行目标体重建,而三种矩阵反演方法对反演效果的影响差别不大;并且通过对波前等时线图的分析可以定性地判断异常体的性质和位置;而在对实测数据目标体的重建上,本方案能达到甚至优于弯曲射线算法的重建效果.     

2D多尺度非线性地震速度成像

将遗传算法和单纯形算法相结合,得到了一种高效、健全的2D混合地震走时反演方法.把速度场划分为不同的空间尺度,定义网格节点上的速度作为待反演参数,采用双三次样条函数模型参数化,正问题采用有限差分走时计算方法,反问题采用多尺度混合反演方法.首先在较大的空间尺度内反演,然后减小空间尺度,将大尺度的反演结果作为次一级尺度反问题的初始模型,再进行混合反演,如此类推逐次逼近全局最优解.一个低速度异常体的数值模拟试验和抗走时扰动试验表明该方法是有效和健全的.我们将该方法应用到青藏高原东北缘阿尼玛卿缝合带东段上部地壳速度结构研究中,并与前人的成果进行了对比.     

Statistical Travel-time Tomography in Terms of Stacking Velocity

— Velocity evaluation is a key step in seismic analysis. The covariance of the true velocity field must be known when interpolating or simulating velocities from well measurements using geostatistical methods. In addition, inversion procedures often require information pertaining to this covariance. Traditionally it has been taken to be the covariance of stacking velocities. We present a simple example to show that this approximation can lead to significant errors. Better methods, such as those of Touati (1996) and Iooss (1998), use the variance of prestack picked travel times as a function of offset to infer that of the velocities. In this paper we extend their results on the estimation of the covariance of the reflected traveltimes, and obtain an explicit expression for the covariance of the square of the stacking slowness as a function of the covariance of the velocities. Although we are not able to invert the formula analytically to yield an explicit estimator for these parameters, the results obtained using it furnish a good and quick estimation of the velocity's covariance. This is illustrated with synthetic examples.     

ELASTIC PARAMETER ESTIMATION BY COHERENCY OPTIMIZATION1

We present a method for estimating P- and S-velocities within defined layers (macromodel), using only kinematic properties (i.e. traveltimes) of the wavefield. The method does not require identification of mixed-mode events on prestack or post-stack data. After obtaining a V_p-depth model by coherency inversion, S-velocities are determined by coherency optimization along computed traveltime curves for mixed-mode events on prestack data. Since the method does not involve any dynamic wavefield computations, a simple ray-tracing algorithm is used to solve the forward problem. The simplicity of the scheme, together with the ability to apply it locally, makes it highly suitable for interactive use. Results of this method may be used to detect Poisson's ratio anomalies within or between layers and may serve as an initial model for more complicated elastic inversion algorithms.