mFAST:A MATLAB toolbox for ocean bottom seismometer refraction first-arrival traveltime tomography

First-arrival seismic traveltime tomography(FAST)is a well-established technique to estimate subsurface velocity structures.Although several existing open-source packages are available for first-arrival traveltime tomography,most were written in compiled languages and lack sufficient extendibility for new algorithms and functionalities.In this work,we develop an open-source,selfcontained FAST package based on MATLAB,one of the most popular interpreted scientific programming languages,with a focus on ocean bottom seismometer refraction traveltime tomography.Our package contains a complete traveltime tomography workflow,including ray-tracing-based first-arrival traveltime computation,linearized inversion,quality control,and high-quality visualization.We design the package as a modular toolbox,making it convenient to integrate new algorithms and functionalities as needed.At the current stage,our package is most efficient for performing FAST for two-dimensional ocean bottom seismometer surveys.We demonstrate the efficacy and accuracy of our package by using a synthetic data example based on a modified Marmousi model.

First-arrival traveltime tomography for monitoring the excavation damaged zone in the Horonobe Underground Research Laboratory

Subsurface excavation results in the formation of a zone called excavation damaged zone(EDZ)around the tunnel wall.An EDZ is a major concern in the field of high-level radioactive waste disposal because it may act as a flow path after the closure of a repository.In this study,first-arrival traveltime tomography was repeatedly conducted on the EDZ at a depth of 350 m in the Horonobe Underground Research Laboratory.However,the acquired data was highly affected by the support structure on the drift wall.For proper visualization of the EDZ,information about the structure was incorporated into the inversion by modifying the model constraint.The synthetic study showed that the approach reproduced the EDZ in the model without the artifacts.The method was applied to field data,and the EDZ around the drift was detected.The inversion was extended to a time-lapse inversion to trace the changes in P-wave velocity in the EDZ.The synthetic study demonstrated that temporal changes in the P-wave velocity distribution could be detected.Data obtained from 12 surveys under open-drift conditions were analyzed by time-lapse inversion.The results indicated that the EDZ did not undergo sealing or evolution at the site for approximately seven years.

An adaptive finite-difference method for seismic traveltime modeling based on 3D eikonal equation

3D eikonal equation is a partial differential equation for the calculation of first-arrival traveltimes and has been widely applied in many scopes such as ray tracing,source localization,reflection migration,seismic monitoring and tomographic imaging.In recent years,many advanced methods have been developed to solve the 3D eikonal equation in heterogeneous media.However,there are still challenges for the stable and accurate calculation of first-arrival traveltimes in 3D strongly inhomogeneous media.In this paper,we propose an adaptive finite-difference(AFD)method to numerically solve the 3D eikonal equation.The novel method makes full use of the advantages of different local operators characterizing different seismic wave types to calculate factors and traveltimes,and then the most accurate factor and traveltime are adaptively selected for the convergent updating based on the Fermat principle.Combined with global fast sweeping describing seismic waves propagating along eight directions in 3D media,our novel method can achieve the robust calculation of first-arrival traveltimes with high precision at grid points either near source point or far away from source point even in a velocity model with large and sharp contrasts.Several numerical examples show the good performance of the AFD method,which will be beneficial to many scientific applications.

Joint tomographic inversion of first-arrival and reflection traveltimes for recovering 2-D seismic velocity structure with an irregular free surface

Irregular surface flattening,which is based on a boundary conforming grid and the transformation between curvilinear and Cartesian coordinate systems,is a mathematical method that can elegantly handle irregular surfaces,but has been limited to obtaining first arrivals only.By combining a multistage scheme with the fast-sweeping method(FSM,the method to obtain first-arrival traveltime in curvilinear coordinates),the reflected waves from a crustal interface can be traced in a topographic model,in which the reflected wavefront is obtained by reinitializing traveltimes in the interface for upwind branches.A local triangulation is applied to make a connection between velocity and interface nodes.Then a joint inversion of first-arrival and reflection traveltimes for imaging seismic velocity structures in complex terrains is presented.Numerical examples all perform well with different seismic velocity models.The increasing topographic complexity and even use of a high curvature reflector in these models demonstrate the reliability,accuracy and robustness of the new working scheme;checkerboard testing illustrates the method's high resolution.Noise tolerance testing indicates the method's ability to yield practical traveltime tomography.Further development of the multistage scheme will allow other later arrivals to be traced and used in the traveltime inversion.

A review of the influencing factors on teleseismic traveltime tomography

Teleseismic traveltime tomography is an important tool for investigating the crust and mantle structure of the Earth.The imaging quality of teleseismic traveltime tomography is affected by many factors,such as mantle heterogeneities,source uncertainties and random noise.Many previous studies have investigated these factors separately.An integral study of these factors is absent.To provide some guidelines for teleseismic traveltime tomography,we discussed four main influencing factors:the method for measuring relative traveltime differences,the presence of mantle heterogeneities outside the imaging domain,station spacing and uncertainties in teleseismic event hypocenters.Four conclusions can be drawn based on our analysis.(1)Comparing two methods,i.e.,measuring the traveltime difference between two adjacent stations(M1)and subtracting the average traveltime of all stations from the traveltime of one station(M2),reveals that both M1 and M2 can well image the main structures;while M1 is able to achieve a slightly higher resolution than M2;M2 has the advantage of imaging long wavelength structures.In practical teleseismic traveltime tomography,better tomography results can be achieved by a two-step inversion method.(2)Global mantle heterogeneities can cause large traveltime residuals(up to about 0.55 s),which leads to evident imaging artifacts.(3)The tomographic accuracy and resolution of M1 decrease with increasing station spacing when measuring the relative traveltime difference between two adjacent stations.(4)The traveltime anomalies caused by the source uncertainties are generally less than 0.2 s,and the impact of source uncertainties is negligible.

Finite-difference calculation of traveltimes based on rectangular grid

To the most of velocity fields, the traveltimes of the first break that seismic waves propagate along rays can be computed on a 2-D or 3-D numerical grid by finite-difference extrapolation. Under ensuring accuracy, to improve calculating efficiency and adaptability, the calculation method of first-arrival traveltime of finite-difference is de- rived based on any rectangular grid and a local plane wavefront approximation. In addition, head waves and scat- tering waves are properly treated and shadow and caustic zones cannot be encountered, which appear in traditional ray-tracing. The testes of two simple models and the complex Marmousi model show that the method has higher accuracy and adaptability to complex structure with strong vertical and lateral velocity variation, and Kirchhoff prestack depth migration based on this method can basically achieve the position imaging effects of wave equation prestack depth migration in major structures and targets. Because of not taking account of the later arrivals energy, the effect of its amplitude preservation is worse than that by wave equation method, but its computing efficiency is higher than that by total Green′s function method and wave equation method.

Finite-difference calculation of traveltimes based on rectangular grid

To the most of velocity fields, the traveltimes of the first break that seismic waves propagate along rays can be computed on a 2-D or 3-D numerical grid by finite-difference extrapolation. Under ensuring accuracy, to improve calculating efficiency and adaptability, the calculation method of first-arrival traveltime of finite-difference is de- rived based on any rectangular grid and a local plane wavefront approximation. In addition, head waves and scat- tering waves are properly treated and shadow and caustic zones cannot be encountered, which appear in traditional ray-tracing. The testes of two simple models and the complex Marmousi model show that the method has higher accuracy and adaptability to complex structure with strong vertical and lateral velocity variation, and Kirchhoff prestack depth migration based on this method can basically achieve the position imaging effects of wave equation prestack depth migration in major structures and targets. Because of not taking account of the later arrivals energy, the effect of its amplitude preservation is worse than that by wave equation method, but its computing efficiency is higher than that by total Green′s function method and wave equation method.

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.

2-D crustal structure of Changbaishan-Tianchi volcanic region determined by seismic traveltime inversion

2-D velocity structure and tectonics of the crust and upper mantle is revealed by inversion of seismic refraction and wide-angle reflection traveltimes acquired along the profile L1 in the Changbaishan-Tianchi volcanic region. It is used in this study that seismic traveltime inversion for simultaneous determination of 2-D velocity and interface structure of the crust and upper mantle. The result shows that, under Changbaishan-Tianchi crater, there exists a low-velocity body in the shape of an inverted triangle, and the crustal reflecting boundaries and Moho all become lower by a varying margin of 2-6 km, forming a crustal root which is assumed to be the Changbaishan-Tianchi volcanic system. Finally, we make a comparison between our 2-D velocity model and the result from the studies by using trial-and-error forward modeling with SEIS83.

A new scheme for fast calculation of seismic traveltimes and ray paths in heterogeneous media

A new method based on Huygens' principle and Fermat's principle is presented to calculate seismic traveltime and ray-paths in this paper. The algorithm can be easily understood and programmed, and can be applied to heteroge neous media in order to eliminate its disadvantage of slowness, we have improved the basic algorithm to speed its calculation to practical stage without the loss of its accuracy.

Upper crustal structure under Jingtai–Hezuo profile in Northeastern Tibet from topography-dependent eikonal traveltime tomography

The Northeastern Tibetan plateau records Caledonian Qilian orogeny and Cenozoic reactivation by continental collision between the Indian and Asian plates. In order to provide the constraint on the Qilian orogenic mechanism and the expansion of the plateau,wide-angle seismic data was acquired along a 430 km-long profile between Jingtai and Hezuo. There is strong height variation along the profile,which is dealt by topography flattening scheme in our crustal velocity structure reconstruction. We herein present the upper crustal P-wave velocity structure model resulting from the interpretation of first arrival dataset from topography-dependent eikonal traveltime tomography. With topography flattening scheme to process real topography along the profile,the evenness of ray coverage times of the image area（upper crust）is improved,which provides upper crustal velocity model comparable to the classic traveltime tomography（with model expansion scheme to process irregular surface）. The upper crustal velocity model shows zoning character which matcheswith the tectonic division of the Qaidam-Kunlun-West Qinling belt,the Central and Northern Qilian,and the Alax blocks along the profile. The resultant upper crustal P-wave velocity model is expected to provide important base for linkage between the mapped surface geology and deep structure or geodynamics in Northeastern Tibet.

Reflection-based traveltime and waveform inversion with second-order optimization

Reflection-based inversion that aims to reconstruct the low-to-intermediate wavenumbers of the subsurface model, can be a complementary to refraction-data-driven full-waveform inversion(FWI), especially for the deep target area where diving waves cannot be acquired at the surface. Nevertheless, as a typical nonlinear inverse problem, reflection waveform inversion may easily suffer from the cycleskipping issue and have a slow convergence rate, if gradient-based first-order optimization methods are used. To improve the accuracy and convergence rate, we introduce the Hessian operator into reflection traveltime inversion(RTI) and reflection waveform inversion(RWI) in the framework of second-order optimization. A practical two-stage workflow is proposed to build the velocity model, in which Gauss-Newton RTI is first applied to mitigate the cycle-skipping problem and then Gauss-Newton RWI is employed to enhance the model resolution. To make the Gauss-Newton iterations more efficiently and robustly for large-scale applications, we introduce proper preconditioning for the Hessian matrix and design appropriate strategies to reduce the computational costs. The example of a real dataset from East China Sea demonstrates that the cascaded Hessian-based RTI and RWI have good potential to improve velocity model building and seismic imaging, especially for the deep targets.

Finite difference calculation of traveltime on non-orthogonal grid

Finite difference methods have been widely employed in solving the eikonal equation so as to calculate traveltime of seismic phase. Most previous studies used regular orthogonal grid. However, much denser grid is required to sample the interfaces that are undulating in depth direction, such as the Moho and the 660 km discontinuity.Here we propose a new finite difference algorithm to solve the eikonal equation on non-orthogonal grid（irregular grid）.To demonstrate its efficiency and accuracy, a test was conducted with a two-layer model. The test result suggests that the similar accuracy of a regular grid with ten times grids could achieve with our new algorithm, but the time cost is only about 0.1 times. A spherical earth model with an undulant660 km discontinuity was constructed to demonstrate the potential application of our new method. In that case, the traveltime curve fluctuation corresponds to topography. Our new algorithm is efficient in solving the first arrival times of waves associated with undulant interfaces.

A Neural Network for Weighted Least－Squares Criteria of Traveltime Tomography

TH network has been applied successfully to linear optimum and some quadratic optimum problems.This paper discusses how to determine the gain function of the net amplifiers and choose appropriate parameters to solve weighted least-squares problems.To test the performance of the net,it is used for seismic traveltime tomography in computer simulation. Two media of different contrast are taken in simulation.The experimental results show that if the parameters are determined appropriately, the performance of the network is good and the results are close to the ideal ones.

Using Traveltime Tomography to Reconstruct Electrical Conductivity

A new method for reconstructing electrical conductivity distribution from electromagnetic (EM) data by using traveltime tomography is presented in this paper. Diffusive EM fields can be mathematically transformed to wavefields defined in a time like variable. The transform uniquely relates a field satisfying a diffusion equation in time, or in frequency, to an integral of the corresponding wavefield. This paper first transforms numerically calculated transient magnetic fields to wavefields. Traveltime data from a source to the receivers are estimated from the transformed wavefields. Then an iterative reconstruction algorithm is used to obtain the slowness distribution of a medium. This algorithm is an improved ART algorithm taking account of bending ray paths. The slowness distribution is transformed to electrical conductivity distribution according to their relation. The simulation result is presented in this paper.

Source-independent wave-equation based microseismic source location using traveltime inversion

A new source location method using wave-equation based traveltime inversion is proposed to locate microseismic events accurately. With a sourceindependent strategy, microseismic events can be located independently regardless of the accuracy of the source signature and the origin time. The traveltime-residuals-based misfit function has robust performance when the velocity model is inaccurate. The new Fréchet derivatives of the misfit function with respect to source location are derived directly based on the acoustic wave equation, accounting for the influence of geometrical perturbation and spatial velocity variation. Unlike the mostly used traveltime inversion methods, no traveltime picking or ray tracing is needed.Additionally, the improved scattering-integral method is applied to reduce the computational cost. Numerical tests show the validity of the proposed method.

An Improved Iterative Reconstruction Algorithm for Traveltime Tomography

Traveltime tomography is a technique to reconstruct acoustic, seismic, or electromagnetic wave-speed distributions from first arrival traveltime data. The ray paths that should be used for tomographic techniques strongly depend on the wave-speed distribution. In this paper, a new method is proposed for finding out the ray paths from Fermat's principle, that means the traveltime of the ray path should be a minimum value. The problem of finding out the ray path is actually an optimum problem. Our new method uses the idea to find out the shortest path in a weighted directed graph to solve the problem. The ray paths found out by this method are used in the iterative reconstruction algorithm. Computer simulation result produced by this reconstruction algorithm is better than that by the conventional ones. It also shows that the new algorithm is effective with good convergency and stability.

Imagingcomplex near-surface structures in Yumen oil field by joint seismic traveltime and waveform inversion

The first-arrival traveltime tomography is a standard approach for near-surface velocity estimation.However,it cannot resolve complex near-surface structures and will produce a smooth velocity model with low resolution.Early arrival waveform inversion is a robust tool for imaging the near surface structures,but it requires a good initial model to avoid cycle skipping between the predicted and observed data.Furthermore,waveform inversion requires substantial computation efforts.Therefore,we present joint seismic traveltime and waveform inversion method,and we expect the joint inversion method retains the advantages of both traveltime inversion and full waveform inversion and overcomes their respective drawbacks at the same time.The objective function includes both the traveltime and waveform misfit.At each iteration,the traveltimes are calculated by wavefront raytracing,and the waveforms are computed using a finite-difference method.The nonlinear optimization problem is solved by the conjugate gradient method.We apply the joint inversion method to study complex near-surface area where shallow overthrust and rugged topography present a significant challenge for applying traveltime inversion and waveform inversion alone.We test synthetic data to verify the advantages of the joint method,and then apply the method to a 2Ddataset acquired in Yumen Oil field,China.The inversion results suggest that the joint traveltime and waveform inversion helps constrain the very shallow velocity structures and also resolve complex overthrust with large velocity contrasts.

Wave equation tomographyin baseband with phase correction from the first arrival traveltimes

We use the first arrival traveltime to correct the phase distortion in a nonlinear wave equation inversion scheme.This improves the precision of tomographic reconstruction of a velocity structure with large variations and helps solve the ill-posed problem of wave equation inversion.When the variation of the velocity distribution is large,general non-linear wave equation inversions are very ill-posed and for such strong nonlinear we can not obtain a correct inversion.One of main reasons is that the calculated and observed phase of the wavefield differs greatly if the initial model is far from the true model.This leads to highly mismatched phase between the calculated and the observed wave field.This is so-called"Cycle Skipping"problem in the full waveform inversion.The phase mismatch is even more pronounced if a high operating frequency is employed in order to increase resolution.To address this problem,we use the first arrival to"demodulate"the wave field in the frequency domain with a goal of restoring the phase of wave field.Then we minimize an objective function consisting of so called"demodulated"wave field to solve wave equation inversion problem.In this way,we find that the inversion is much improved,and when the velocity perturbation in a complicated model reaches 35%,we can still obtain a good inversion.A computer simulation shows that our method is very robust for acoustical wave inversion with good reconstruction precision.

反射波法隧道超前地质预报地震波走时模拟研究

隧道施工安全常常由于复杂的地质情况面临着巨大的挑战,反射波法地震勘探具有探测深度大、抗干扰能力强等优点,是确保隧道施工安全的重要超前地质预报手段。反射地震波走时计算与地震数据处理中的纵横波分离、速度建模、偏移成像等多个步骤息息相关,直接影响着隧道地质预报的结果。本文将仿照反射波法隧道超前地质预报的实际施工区域进行速度建模,以求解程函方程为基础,结合快速推进法(FMM)中的迎风有限差分方法以及窄带扩展技术实现针对隧道模型的初至地震波走时模拟。并通过对隧道模型中反射界面位置信息的提取,构建反射波初始窄带以及结合窄带扩展技术实现针对隧道模型反射波走时的计算。通过数值分析以及复杂隧道模型计算表明,文中给出的反射波法隧道超前地质预报地震波走时模拟具有良好的计算精度以及较强的适应能力。