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.

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.

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.

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.

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.

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.

Wavefield reconstruction inversion with weaktotal-variation constraint based on first arrival traveltime tomography

The objective function of full waveform inversion is a strong nonlinear function,the inversion process is not unique,and it is easy to fall into local minimum.Firstly,in the process of wavefield reconstruction,the wave equation is introduced into the construction of objective function as a penalty term to broaden the search space of solution and reduce the risk of falling into local minimum.In addition,there is no need to calculate the adjoint wavefield in the inversion process,which can significantly improve the calculation efficiency;Secondly,considering that the total variation constraint can effectively reconstruct the discontinuous interface in the velocity model,this paper introduces the weak total variation constraint to avoid the excessive smooth estimation of the model under the strong total variation constraint.The disadvantage of this strategy is that it is highly dependent on the initial model.In view of this,this paper takes the long wavelength initial model obtained by first arrival traveltime tomography as a prior model constraint,and proposes a weak total variation constrained wavefield reconstruction inversion method based on first arrival traveltime tomography.Numerical experimental results show that the new method reduces the dependence on the initial model,the interface description is more accurate,the error is reduced,and the iterative convergence efficiency is significantly improved.

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.

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.

海水速度结构的OBS透射走时层析反演

近年来,海洋高分辨率地震勘探领域取得了显著进展,对地震资料反演的精度要求日益严苛.为了提升后续处理的效率和准确性,获取更为精确的海水速度分布显得至关重要.然而,传统的走时层析方法在利用有限偏移距的海面拖缆数据时,往往难以达到理想的垂向分辨率,从而限制了海水速度建模的精确性.鉴于此,本文提出了一种基于OBS(海底地震仪)观测系统的透射波走时层析方法.该方法采用伴随状态法走时层析作为反演工具,利用OBS观测系统声压分量中的海水透射波走时信息,对海水速度结构进行精确反演.理论模型和实际资料实验结果表明,该方法能够准确建立海水速度模型,提高海水以下地层的成像效果.

全球非均匀模型中的远震走时层析成像

远震走时层析成像是研究地球内部结构的重要工具.通常,远震走时层析成像使用一维地球参考模型计算研究区域外地震射线.由于地球存在广泛非均匀结构,传统远震走时层析成像在计算射线时忽略了研究区域外非均匀结构对射线的影响,最终导致对研究区域内的成像结果产生影响.为了减小模型外地球非均匀结构对远震走时层析成像的影响,本文利用迎风差分格式数值求解程函方程,得到全球非均匀模型中的远震射线,以此为基础发展了基于全球非均匀模型的远震走时成像方法.为了检验本文成像方法提升成像精度的有效性,选取青藏高原东北缘为研究对象.合成测试和实际数据成像结果显示,在考虑地球非均匀结构以后,成像结果更精确,结果展现更多成像细节.

基于声波测井资料的波动方程走时层析反演速度建模

地震勘探目标逐渐由浅层(<2000 m)转向中深层(2000~3500 m),高精度速度建模是实现精确地震成像的关键手段之一.目前,初至波或早至波速度层析建模方法已较为成熟,然而对深层速度建模更为重要的反射波层析建模方法依然面临走时拾取困难及计算量大等难题.针对当前中深层走时层析速度建模方法面临的难题,本文提出了一种基于声波测井资料的波动方程全波走时层析速度建模方法及策略.该方法首先通过早至波层析反演获得浅层及大尺度速度信息;其次,通过提取初次成像剖面的构造倾角信息,利用声波测井资料进行稀疏反演,进一步获得更高精度的速度模型;最后,以此为初始模型,进行全波层析反演获得最终偏移速度模型.数值结果表明,本文所提出的速度建模方法及策略不仅可以避免传统反射波层析建模方法因采用Born近似反偏移所带来的大计算量问题,且大幅提高了中深层速度建模的精度,进而能有效地提高中深层复杂构造区域的地震成像质量.

海底地震仪与拖揽观测系统的联合走时层析海水速度建模

在海洋地震勘探过程中,如果不考虑海水的非均匀性将会在后续层析成像过程中引入伪影,因此准确的海水速度建模是非常有必要的.海底地震仪(OBS)记录数据中的透射波具有较大的射线角度,能够为海水建模提供垂向分辨率;而拖缆观测可以补充OBS分布稀疏造成的横向欠采样,同时可以定位海底深度并更新OBS位置信息.综合利用OBS和拖缆观测系统可以克服各自的局限性,提高海水速度建模的精度.为此,本文提出了一种基于OBS与拖缆观测系统的海水透射波、海底反射波联合走时层析方法,同时反演海水声速分布和海底深度.二维理论模型实验和实际资料处理结果表明,该方法能够准确反演海水声速分布和海底深度,提高震相走时的匹配程度,为后续成像提供准确的非均匀海水模型.

3D near-surface P-wave velocity structure imaging with Distributed Acoustic Sensing and electric hammer source

Distributed Acoustic Sensing(DAS) is an emerging technique for ultra-dense seismic observation, which provides a new method for high-resolution sub-surface seismic imaging. Recently a large number of linear DAS arrays have been used for two-dimensional S-wave near-surface imaging in urban areas. In order to explore the feasibility of three-dimensional(3D) structure imaging using a DAS array, we carried out an active source experiment at the Beijing National Earth Observatory. We deployed a 1 km optical cable in a rectangular shape, and the optical cable was recast into 250 sensors with a channel spacing of 4 m. The DAS array clearly recorded the P, S and surface waves generated by a hammer source. The first-arrival P wave travel times were first picked with a ShortTerm Average/Long-Term Average(STA/LTA) method and further manually checked. The P-wave signals recorded by the DAS are consistent with those recorded by the horizontal components of short-period seismometers. At shorter source-receiver distances, the picked P-wave arrivals from the DAS recording are consistent with vertical component recordings of seismometers, but they clearly lag behind the latter at greater distances.This is likely due to a combination of the signal-to-noise ratio and the polarization of the incoming wave. Then,we used the Tomo DD software to invert the 3D P-wave velocity structure for the uppermost 50 m with a resolution of 10 m. The inverted P-wave velocity structures agree well with the S-wave velocity structure previously obtained through ambient noise tomography. Our study indicates the feasibility of 3D near-surface imaging with the active source and DAS array. However, the inverted absolute velocity values at large depths may be biased due to potential time shifts between the DAS recording and seismometer at large source-receiver distances.

Solving Traveltime Tomography with Deep Learning

This paper introduces a neural network approach for solving two-dimensional traveltime tomography(TT)problems based on the eikonal equation.The mathematical problem of TT is to recover the slowness field of a medium based on the boundary measurement of the traveltimes of waves going through the medium.This inverse map is high-dimensional and nonlinear.For the circular tomography geometry,a perturbative analysis shows that the forward map can be approximated by a vectorized convolution operator in the angular direction.Motivated by this and filtered backprojection,we propose an effective neural network architecture for the inverse map using the recently proposed BCR-Net,with weights learned from training datasets.Numerical results demonstrate the efficiency of the proposed neural networks.

Application of teleseismic tomography to the study of shallow structure beneath Shizigou in the western Qaidam basin

Teleseismic body wave traveltime tomography is used to inverse the three-dimensional seismic velocity structure beneath Shizigou in the western Qaidam basin. The travel time are picked from the continuous observation data on a small seismic array of stations deployed during 2004-2007. The tomographic results obtained indicate that a NW-trending low velocity anomaly just beneath the target region insert northeastwards with a high dip angle. In the north, northeast and east of the low velocity anomaly, some high-velocity anomalies distribute with the same strike and coverage as those of Shizigou anticline.

Teleseismic P-Wave Tomography of the New Guinea-Solomon Arc System

A P-wave tomographic traveltime inversion was applied to obtain a new model of seismic velocity anomalies beneath the New Guinea-Solomon arc system(PN-SL).The P-wave traveltime data,obtained from the revised International Seismological Center catalog,were recorded by 82 seismic stations in the PN-SL.Under the constraints of the epicenter distance,magnitude,and the number of stations recorded,15009 effective P-wave traveltime data were selected from 2011 teleseismic events.The obtained model showed that the Solomon Sea Plate subducted beneath the New Britain Island along the New Britain Trench at an angle of>70°and that the slab can be traced down to a depth of approximately 800 km.Conversely,we cannot observe a high-velocity anomaly exhibited by the subducted Solomon Sea Plate in the deep mantle at the Trobriand Trench,and the slab stopped at a depth of<200 km.The double subduction of the Solomon Sea Plate strongly modified the subduction patterns of the early subducted Pacific and Australian plates in the mantle along the West Melanesian Trench and the Pocklington Trough,respectively.In addition,the subducted Solomon Sea Plate induced the melting of the upper mantle to form a low-velocity anomaly,which provided the deep dynamic source for the expansion of the Bismarck Sea.Based on the joint consideration of the tomography results and a petrological analysis,the low-velocity anomalies beneath the Solomon Sea and Woodlark Basin are closely related to the early subduction of the Pacific and Australian Plates,respectively.

基于单频走时敏感度核函数的三维波动方程初至层析

地震波走时层析是主要的近地表速度建模技术.对于变速模型,尤其是复杂近地表模型,波动方程能够更加精确地描述地震波的一阶绕射效应,理论上反演精度高于射线(束)类层析.然而,波动方程的数值求解的计算量较射线理论增加了几个数量级,严重制约了波动方程走时层析技术的三维实用化.为解决波动方程走时层析中遇到的计算和存储问题,本文通过引入随机边界条件以及离散Fourier积分,提出了基于单频梯度反演策略的实用化三维波动方程初至层析方法.本文方法的优势在于:在计算和存储方面,相比于基于波场重构算法(需要三次波传播)的梯度计算策略,本文提出的单频梯度计算方法对内存需求极低(相对于波动方程正演,仅需要增加两个单频波场以及一个成像体),且计算量至少减少1/3(仅需要两次波传播且无需处理数值吸收边界).此外,随机边界的引入使得正演算法大幅简化因而更适用于GPU(Graphics Processing Unit,图形处理器)加速.在理论层面,由于影响地震波一阶散射效应的主要区域集中在连接炮检的中心射线邻域的第一菲涅尔带内,随机边界条件的引入降低了随机边界内产生的散射波的空间相关性,使得第一菲涅尔带内的单频梯度相干加强;同时,地震观测系统的多次覆盖特性有助于第一菲涅尔带外的高波数振荡相互干涉.数值实验表明,光滑处理后的单频梯度仍然保留了非常可观的低波数成分,有助于实现背景速度反演.此外,三维复杂模型反演结果验证了本文方法能够较为准确地反演近地表速度结构.相对于传统的波动方程走时层析实现方案,本文方法在显著降低内存需求的同时进一步提高了计算效率,有望将波动方程初至层析技术推向三维实用化.