Non‐uniqueness with refraction inversion – the Mt Bulga shear zone

The tau‐p inversion algorithm is widely employed to generate starting models with many computer programs that implement refraction tomography. However, this algorithm can frequently fail to detect even major lateral variations in seismic velocities, such as a 50 m wide shear zone, which is the subject of this study. By contrast, the shear zone is successfully defined with the inversion algorithms of the generalized reciprocal method. The shear zone is confirmed with a 2D analysis of the head wave amplitudes, a spectral analysis of the refraction convolution section and with numerous closely spaced orthogonal seismic profiles recorded for a later 3D refraction investigation. Further improvements in resolution, which facilitate the recognition of additional zones with moderate reductions in seismic velocity, are achieved with a novel application of the Hilbert transform to the refractor velocity analysis algorithm. However, the improved resolution also requires the use of a lower average vertical seismic velocity, which accommodates a velocity reversal in the weathering. The lower seismic velocity is derived with the generalized reciprocal method, whereas most refraction tomography programs assume vertical velocity gradients as the default. Although all of the tomograms are consistent with the traveltime data, the resolution of each tomogram is comparable only with that of the starting model. Therefore, it is essential to employ inversion algorithms that can generate detailed starting models, where detailed lateral resolution is the objective. Non‐uniqueness can often be readily resolved with head wave amplitudes, attribute processing of the refraction convolution section and additional seismic traverses, prior to the acquisition of any borehole data. It is concluded that, unless specific measures are taken to address non‐uniqueness, the production of a single refraction tomogram that fits the traveltime data to sufficient accuracy does not necessarily demonstrate that the result is either correct, or even the most probable.     

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.     

P‐wave velocity distribution in basalt flows of the Enni Formation in the Faroe Islands from refraction seismic analysis

The main objective of this work is to establish the applicability of shallow surface‐seismic traveltime tomography in basalt‐covered areas. A densely sampled ～1300‐m long surface seismic profile, acquired as part of the SeiFaBa project in 2003 ( Japsen et al. 2006 ) at Glyvursnes in the Faroe Islands, served as the basis to evaluate the performance of the tomographic method in basalt‐covered areas. The profile is centred at a ～700‐m deep well. V_P, V_S and density logs, a zero‐offset VSP, downhole‐geophone recordings and geological mapping in the area provided good means of control. The inversion was performed with facilities of the Wide Angle Reflection/Refraction Profiling program package ( Ditmar et al. 1999 ). We tested many inversion sequences while varying the inversion parameters. Modelled traveltimes were verified by full‐waveform modelling. Typically an inversion sequence consists in several iterations that proceed until a satisfactory solution is reached. However, in the present case with high velocity contrasts in the subsurface we obtained the best result with two iterations: first obtaining a smooth starting model with small traveltime residuals by inverting with a high smoothing constraint and then inverting with the lowest possible smoothing constraint to allow the inversion to have the full benefit of the traveltime residuals. The tomogram gives usable velocity information for the near‐surface geology in the area but fails to reproduce the expected velocity distribution of the layered basalt flows. Based on the analysis of the tomogram and geological mapping in the area, a model was defined that correctly models first arrivals from both surface seismic data and downhole‐geophone data.     

Kinematical characteristics of the factorized velocity model

The factorized velocity model that incorporates both vertical heterogeneity and constant anisotropy is one of the complicated analytical models used in seismic data processing and interpretation. In this paper, I derive the analytic equations for offset, traveltime and relative geometrical spreading for the quasi‐compressional (qP‐) waves that can be used for modelling and inversion of the traveltime parameters. I show that the presence of anelliptic anisotropy usually dominates over the vertical heterogeneity with respect to the non‐hyperbolicity of the factorized velocity model.     

Assessing uncertainty in refraction seismic traveltime inversion using a global inversion strategy

To analyse and invert refraction seismic travel time data, different approaches and techniques have been proposed. One common approach is to invert first‐break travel times employing local optimization approaches. However, these approaches result in a single velocity model, and it is difficult to assess the quality and to quantify uncertainties and non‐uniqueness of the found solution. To address these problems, we propose an inversion strategy relying on a global optimization approach known as particle swarm optimization. With this approach we generate an ensemble of acceptable velocity models, i.e., models explaining our data equally well. We test and evaluate our approach using synthetic seismic travel times and field data collected across a creeping hillslope in the Austrian Alps. Our synthetic study mimics a layered near‐surface environment, including a sharp velocity increase with depth and complex refractor topography. Analysing the generated ensemble of acceptable solutions using different statistical measures demonstrates that our inversion strategy is able to reconstruct the input velocity model, including reasonable, quantitative estimates of uncertainty. Our field data set is inverted, employing the same strategy, and we further compare our results with the velocity model obtained by a standard local optimization approach and the information from a nearby borehole. This comparison shows that both inversion strategies result in geologically reasonable models (in agreement with the borehole information). However, analysing the model variability of the ensemble generated using our global approach indicates that the result of the local optimization approach is part of this model ensemble. Our results show the benefit of employing a global inversion strategy to generate near‐surface velocity models from refraction seismic data sets, especially in cases where no detailed a priori information regarding subsurface structures and velocity variations is available.     

Sensitivity analysis and application of time‐lapse full‐waveform inversion: synthetic testing and field data example from the North Sea,Norway

Time‐lapse refraction can provide complementary seismic solutions for monitoring subtle subsurface changes that are challenging for conventional P‐wave reflection methods. The utilization of refraction time lapse has lagged behind in the past partly due to the lack of robust techniques that allow extracting easy‐to‐interpret reservoir information. However, with the recent emergence of the full‐waveform inversion technique as a more standard tool, we find it to be a promising platform for incorporating head waves and diving waves into the time‐lapse framework. Here we investigate the sensitivity of 2D acoustic, time‐domain, full‐waveform inversion for monitoring a shallow, weak velocity change (?30 m/s, or ?1.6%). The sensitivity tests are designed to address questions related to the feasibility and accuracy of full‐waveform inversion results for monitoring the field case of an underground gas blowout that occurred in the North Sea. The blowout caused the gas to migrate both vertically and horizontally into several shallow sand layers. Some of the shallow gas anomalies were not clearly detected by conventional 4D reflection methods (i.e., time shifts and amplitude difference) due to low 4D signal‐to‐noise ratio and weak velocity change. On the other hand, full‐waveform inversion sensitivity analysis showed that it is possible to detect the weak velocity change with the non‐optimal seismic input. Detectability was qualitative with variable degrees of accuracy depending on different inversion parameters. We inverted, the real 2D seismic data from the North Sea with a greater emphasis on refracted and diving waves’ energy (i.e., most of the reflected energy was removed for the shallow zone of interest after removing traces with offset less than 300 m). The full‐waveform inversion results provided more superior detectability compared with the conventional 4D stacked reflection difference method for a weak shallow gas anomaly (320 m deep).     

模拟退火方法在三维速度模型地震波走时反演中的应用

采用块状建模以及三角形拼接的界面描述方式,并通过立方体速度网格线性插值获得块体内部的速度分布。正演过程中采用逐段迭代射线追踪方法计算三维复杂地质模型中的射线走时,并采用模拟退火方法进行了三维模型中的地震波走时反演研究。模型测试结果表明,使用的射线追踪和走时反演算法有效。     

3D P‐wave traveltime computation in transversely isotropic media using layer‐by‐layer wavefront marching

Subsurface rocks (e.g. shale) may induce seismic anisotropy, such as transverse isotropy. Traveltime computation is an essential component of depth imaging and tomography in transversely isotropic media. It is natural to compute the traveltime using the wavefront marching method. However, tracking the 3D wavefront is expensive, especially in anisotropic media. Besides, the wavefront marching method usually computes the traveltime using the eikonal equation. However, the anisotropic eikonal equation is highly non‐linear and it is challenging to solve. To address these issues, we present a layer‐by‐layer wavefront marching method to compute the P‐wave traveltime in 3D transversely isotropic media. To simplify the wavefront tracking, it uses the traveltime of the previous depth as the boundary condition to compute that of the next depth based on the wavefront marching. A strategy of traveltime computation is designed to guarantee the causality of wave propagation. To avoid solving the non‐linear eikonal equation, it updates traveltime along the expanding wavefront by Fermat's principle. To compute the traveltime using Fermat's principle, an approximate group velocity with high accuracy in transversely isotropic media is adopted to describe the ray propagation. Numerical examples on 3D vertical transverse isotropy and tilted transverse isotropy models show that the proposed method computes the traveltime with high accuracy. It can find applications in modelling and depth migration.     

Piecewise‐continuous velocity inversion

We suggest a new method to determine the piecewise‐continuous vertical distribution of instantaneous velocities within sediment layers, using different order time‐domain effective velocities on their top and bottom points. We demonstrate our method using a synthetic model that consists of different compacted sediment layers characterized by monotonously increasing velocity, combined with hard rock layers, such as salt or basalt, characterized by constant fast velocities, and low velocity layers, such as gas pockets. We first show that, by using only the root‐mean‐square velocities and the corresponding vertical travel times (computed from the original instantaneous velocity in depth) as input for a Dix‐type inversion, many different vertical distributions of the instantaneous velocities can be obtained (inverted). Some geological constraints, such as limiting the values of the inverted vertical velocity gradients, should be applied in order to obtain more geologically plausible velocity profiles. In order to limit the non‐uniqueness of the inverted velocities, additional information should be added. We have derived three different inversion solutions that yield the correct instantaneous velocity, avoiding any a priori geological constraints. The additional data at the interface points contain either the average velocities (or depths) or the fourth‐order average velocities, or both. Practically, average velocities can be obtained from nearby wells, whereas the fourth‐order average velocity can be estimated from the quartic moveout term during velocity analysis. Along with the three different types of input, we consider two types of vertical velocity models within each interval: distribution with a constant velocity gradient and an exponential asymptotically bounded velocity model, which is in particular important for modelling thick layers. It has been shown that, in the case of thin intervals, both models lead to similar results. The method allows us to establish the instantaneous velocities at the top and bottom interfaces, where the velocity profile inside the intervals is given by either the linear or the exponential asymptotically bounded velocity models. Since the velocity parameters of each interval are independently inverted, discontinuities of the instantaneous velocity at the interfaces occur naturally. The improved accuracy of the inverted instantaneous velocities is particularly important for accurate time‐to‐depth conversion.     

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

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

Refraction traveltime and amplitude corrections for very near- surface inhomogeneities

The performance of refraction inversion methods that employ the principle of refraction migration, whereby traveltimes are laterally migrated by the offset distance (which is the horizontal separation between the point of refraction and the point of detection on the surface), can be adversely affected by very near‐surface inhomogeneities. Even inhomogeneities at single receivers can limit the lateral resolution of detailed seismic velocities in the refractor. The generalized reciprocal method ‘statics’ smoothing method (GRM SSM) is a smoothing rather than a deterministic method for correcting very near‐surface inhomogeneities of limited lateral extent. It is based on the observation that there are only relatively minor differences in the time‐depths to the target refractor computed for a range of XY distances, which is the separation between the reverse and forward traveltimes used to compute the time‐depth. However, any traveltime anomalies, which originate in the near‐surface, migrate laterally with increasing XY distance. Therefore, an average of the time‐depths over a range of XY values preserves the architecture of the refractor, but significantly minimizes the traveltime anomalies originating in the near‐surface. The GRM statics smoothing corrections are obtained by subtracting the average time‐depth values from those computed with a zero XY value. In turn, the corrections are subtracted from the traveltimes, and the GRM algorithms are then re‐applied to the corrected data. Although a single application is generally adequate for most sets of field data, model studies have indicated that several applications of the GRM SSM can be required with severe topographic features, such as escarpments. In addition, very near‐surface inhomogeneities produce anomalous head‐wave amplitudes. An analogous process, using geometric means, can largely correct amplitude anomalies. Furthermore, the coincidence of traveltime and amplitude anomalies indicates that variations in the near‐surface geology, rather than variations in the coupling of the receivers, are a more likely source of the anomalies. The application of the GRM SSM, together with the averaging of the refractor velocity analysis function over a range of XY values, significantly minimizes the generation of artefacts, and facilitates the computation of detailed seismic velocities in the refractor at each receiver. These detailed seismic velocities, together with the GRM SSM‐corrected amplitude products, can facilitate the computation of the ratio of the density in the bedrock to that in the weathered layer. The accuracy of the computed density ratio improves where lateral variations in the seismic velocities in the weathered layer are known.     

Quasi-2D hybrid joint inversion of seismic and geoelectric data

Many joint inversion schemes use 1D forward modelling in the integrated interpretation of various geophysical data. In extending the joint inversion approach to the investigation of 2D structures, the discretization of the model parameters and the appropriate choice of the forward‐modelling procedure play a very important role. In this paper, a hybrid seismic–geoelectric joint inversion method is proposed for the investigation of 2D near‐surface geological structures. The electric and seismic models are coupled together through the use of common boundaries between the adjacent layers. Assuming a 2D model composed of homogeneous layers with curved boundaries, a fast ray‐tracing algorithm is used for the calculation of refraction seismic traveltime data. In the geoelectric forward modelling, a locally 1D approximation is used. The boundary surfaces are written in the form of series expansion; the inversion algorithms are formulated for the expansion coefficients and the petrophysical parameters as unknowns. Two versions of the inversion method are proposed: in versions A and B, interval‐wise constant functions and Chebyshev polynomials are, respectively, used as basis functions of the series expansion. The versions are tested by means of synthetic and in situ measured data. The tests show that both methods are stable and accurate.     

THE RESOLUTION OF NARROW LOW-VELOCITY ZONES WITH THE GENERALIZED RECIPROCAL METHOD1

The depth to the surface of a refractor and the seismic velocity within the refractor are very often intimately related. In the shallow environment, increased thicknesses of weathering occur in areas of jointing, shearing or lithological variations, and these zones of deeper weathering can have lower subweathering refractor velocities. This association is important in geotechnical investigations and in the measurement of weathering thicknesses and sub-weathering velocities for statics corrections for reflection seismic surveys. Algorithms, which employ forward and reverse traveltime data and which explicitly accommodate the offset distance through the process known as refraction migration, are necessary if detailed structure on a refractor and rapid lateral variations of the seismic velocity within it are to be resolved. These requirements are satisfied with wavefront construction techniques, Hales’ method and the generalized reciprocal method (GRM). However, these methods employ refraction migration in fundamentally different manners. Most methods compute an offset distance with an often imprecise knowledge of the seismic velocities of the overlying layers. In contrast, the GRM uses a range of offset distances from less than to greater than the optimum value, with the optimum value being selected with a minimum-variance criterion. The approach of the GRM is essential where there are undetected layers and where there are rapid variations in the depth to a refractor and the seismic velocity within it. In the latter situations the offset distance necessary to define the seismic velocities can differ considerably from the value required to define depths. The efficacy of the GRM in resolving structure and seismic velocity is demonstrated with three model studies and two field examples.     

Misfit functionals in Laplace‐Fourier domain waveform inversion,with application to wide‐angle ocean bottom seismograph data

In seismic waveform inversion, non‐linearity and non‐uniqueness require appropriate strategies. We formulate four types of L2 normed misfit functionals for Laplace‐Fourier domain waveform inversion: i) subtraction of complex‐valued observed data from complex‐valued predicted data (the ‘conventional phase‐amplitude’ residual), ii) a ‘conventional phase‐only’ residual in which amplitude variations are normalized, iii) a ‘logarithmic phase‐amplitude’ residual and finally iv) a ‘logarithmic phase‐only’ residual in which the only imaginary part of the logarithmic residual is used. We evaluate these misfit functionals by using a wide‐angle field Ocean Bottom Seismograph (OBS) data set with a maximum offset of 55 km. The conventional phase‐amplitude approach is restricted in illumination and delineates only shallow velocity structures. In contrast, the other three misfit functionals retrieve detailed velocity structures with clear lithological boundaries down to the deeper part of the model. We also test the performance of additional phase‐amplitude inversions starting from the logarithmic phase‐only inversion result. The resulting velocity updates are prominent only in the high‐wavenumber components, sharpening the lithological boundaries. We argue that the discrepancies in the behaviours of the misfit functionals are primarily caused by the sensitivities of the model gradient to strong amplitude variations in the data. As the observed data amplitudes are dominated by the near‐offset traces, the conventional phase‐amplitude inversion primarily updates the shallow structures as a result. In contrast, the other three misfit functionals eliminate the strong dependence on amplitude variation naturally and enhance the depth of illumination. We further suggest that the phase‐only inversions are sufficient to obtain robust and reliable velocity structures and the amplitude information is of secondary importance in constraining subsurface velocity models.     

三维复杂地形近地表速度估算及地震层析静校正

在地表一致性模型的基础上提出一种可适用于宽线剖面、弯曲测线、传统的二维和目前广泛使用的三维地震观测.在地形及近地表低降速带地质结构复杂的探区，低降速带厚度及速度估算的精度是静校正处理的关键.本研究根据三维地震观测的初至走时数据，利用最小平方与QR分解相结合的算法，在三维空间重建近地表低降速带速度模型，根据重建速度模型实现了静校正长波长分量与短波长分量的同步计算.分析了复杂的近地表低降速带模型初至波的性质，在观测值的自动拾取以及理论值的计算中充分考虑了可能成为初至波的直达波、折射波和反射波的利用，提高了低降速带速度模型反演的精度.在初至走时观测数据的拾取中，本研究采用分形算法克服了初至波波形差异以及折射波相位反转导致的拾取误差，实现了三维初至拾取的大规模全自动化运算.在射线路径与初至波理论走时的计算中，本研究采用一种计算量与模型复杂程度无关的三维射线追踪方法，该方法以最小走时射线路径保证了与观测数据有同等意义的初至波的射线追踪及理论走时的计算.野外实际资料的处理结果表明了方法的有效性.     

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.     

Solving the complex near‐surface problem using 3D data‐driven near‐surface layer replacement

In many land seismic situations, the complex seismic wave propagation effects in the near‐surface area, due to its unconsolidated character, deteriorate the image quality. Although several methods have been proposed to address this problem, the negative impact of 3D complex near‐surface structures is still unsolved to a large extent. This paper presents a complete 3D data‐driven solution for the near‐surface problem based on 3D one‐way traveltime operators, which extends our previous attempts that were limited to a 2D situation. Our solution is composed of four steps: 1) seismic wave propagation from the surface to a suitable datum reflector is described by parametrized one‐way propagation operators, with all the parameters estimated by a new genetic algorithm, the self‐adjustable input genetic algorithm, in an automatic and purely data‐driven way; 2) surface‐consistent residual static corrections are estimated to accommodate the fast variations in the near‐surface area; 3) a replacement velocity model based on the traveltime operators in the good data area (without the near‐surface problem) is estimated; 4) data interpolation and surface layer replacement based on the estimated traveltime operators and the replacement velocity model are carried out in an interweaved manner in order to both remove the near‐surface imprints in the original data and keep the valuable geological information above the datum. Our method is demonstrated on a subset of a 3D field data set from the Middle East yielding encouraging results.     

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

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

The use of low frequencies in a full‐waveform inversion and impedance inversion land seismic case study

Velocity model building and impedance inversion generally suffer from a lack of intermediate wavenumber content in seismic data. Intermediate wavenumbers may be retrieved directly from seismic data sets if enough low frequencies are recorded. Over the past years, improvements in acquisition have allowed us to obtain seismic data with a broader frequency spectrum. To illustrate the benefits of broadband acquisition, notably the recording of low frequencies, we discuss the inversion of land seismic data acquired in Inner Mongolia, China. This data set contains frequencies from 1.5–80 Hz. We show that the velocity estimate based on an acoustic full‐waveform inversion approach is superior to one obtained from reflection traveltime inversion because after full‐waveform inversion the background velocity conforms to geology. We also illustrate the added value of low frequencies in an impedance estimate.     

Seismic traveltime inversion of 3D velocity model with triangulated interfaces

Seismic traveltime tomographic inversion has played an important role in detecting the internal structure of the solid earth. We use a set of blocks to approximate geologically complex media that cannot be well described by layered models or cells. The geological body is described as an aggregate of arbitrarily shaped blocks, which are separated by triangulated interfaces. We can describe the media as homogenous or heterogeneous in each block. We define the velocities at the given rectangle grid points for each block, and the heterogeneous velocities in each block can be calculated by a linear interpolation algorithm. The parameters of the velocity grid positions are independent of the model parameterization, which is advantageous in the joint inversion of the velocities and the node depths of an interface. We implement a segmentally iterative ray tracer to calculate traveltimes in the 3D heterogeneous block models. The damped least squares method is employed in seismic traveltime inversion, which includes the partial derivatives of traveltime with respect to the depths of nodes in the triangulated interfaces and velocities defined in rectangular grids. The numerical tests indicate that the node depths of a triangulated interface and homogeneous velocity distributions can be well inverted in a stratified model.