Subsurface Velocity Structure Research Articles

Diffraction-angle filtering of gradient for acoustic full-waveform inversion

Seismic full-waveform inversion (FWI) estimates subsurface velocity structures by reducing data misfit between observed and modeled data. Simultaneous matching of transmitted and reflected waves in seismic FWI causes different updates of different wavenumber components of a given model depending on the diffraction angle between incident and diffracted rays. Motivated by the inverse-scattering imaging condition and elastic FWI, we have applied a diffraction-angle filtering technique in acoustic FWI, which enables us to separate transmission and reflection energy in the partial derivative wavefields. Diffraction-angle filtering is applied to the virtual source, which is the model parameter perturbation acting as a source for the partial derivative wavefields. Diffraction-angle filtering consists of two diffraction-angle filters (DAFs), DAF-I and DAF-II. The former is derived from the particle acceleration of the incident wavefields and suppresses energies at either small or large diffraction angles by changing the sign of the weighting factor. DAF-I is exactly identical to the conventional inverse scattering approach. The DAF-II is derived from the artificial shear strain of the incident P-wave and additionally suppresses energies at intermediate diffraction angles. With this mechanism, we can design various types of diffraction-angle filtering to control the updates of wavenumber components of the misfit gradient with respect to the P-wave velocity. Using the synthetic Marmousi-II data and real ocean-bottom seismic data from the North Sea, we determined that diffraction-angle filtering enables us to control low-, intermediate-, and high-wavenumber components of the gradient direction.

Imaging Strong Lateral Heterogeneities Across the Contiguous US Using Body‐To‐Surface Wave Scattering

AbstractBody‐to‐surface wave scattering, originated from strong lateral heterogeneity, has been observed and modeled for decades. Compared to body waves, scattered surface waves propagate along the Earth's surface with less energy loss and, thus, can be observed over a wider distance range. In this study, we utilize surface waves converted from teleseismic SH or Sdiff wave incidence to map strong lateral heterogeneities across the entire contiguous United States. We apply array‐based phase coherence analysis to broadband waveforms recorded by the USArray Transportable Array and other permanent/temporary networks to detect coherent signals that are associated with body‐to‐surface wave scattering. We then locate the source of the scattering by back‐propagating the beamformed energy using both straight‐ray and curved‐ray approximations. Our results show that the distribution of scatterers correlates well with known geological features across the contiguous United States. Topographic/bathymetric relief along the continental slope off the Pacific Border is the major source of scattering in the western United States. On the other hand, sedimentary basins, especially their margins, are the dominant scatterers in the central United States. Moho offsets, such as the one around the periphery of the Colorado Plateau, are also a strong contributor to scattering, but isolating their effect from that of other near‐surface structures without any additional constraints can be complicated. Finally, we demonstrate the possibility of using scattered surface waves to constrain subsurface velocity structures, as complementary to conventional earthquake‐ or ambient noise‐based surface wave tomography.

Estimation of shallow subsurface S-wave velocity structure in urban area based on inversion of apparent dispersion curve

AbstractDue to its efficiency, convenience, non-destructive nature and strong anti-interference capability, the microtremor survey method (MSM) has found wide applications in urban geological surveys. The spatial autocorrelation method is diffusely applied to extract the dispersion curves from microtremor signals, which needs to satisfy the assumption that the energy of the fundamental Rayleigh wave is dominant. However, for layered media containing a layer with a significant low- or high-velocity contrast, this assumption is distinctly incorrect for certain frequency ranges. We present a processing methodology comprising the extraction and inversion of the apparent dispersion curves based on extended spatial autocorrelation method and fast simulated-annealing algorithm. We analyse synthetic microtremor signals for three selected geological models, and then compare the S-wave velocity structures estimated from their inversions with the actual models. Subsequently, a filed data example is given to detect the shallow stratigraphic structures in Guangzhou city, China, in which the new MSM was used. The estimated two-dimensional S-wave velocity model provided an accurate description of the thickness and depth of the strata in the study area, based on a priori information. Moreover, the S-wave velocity structures estimated from the MSM and the results from the drilling match very well, indicating that MSM is a reliable geophysical technique in urban geological surveys. Combined with available borehole information, MSM can be a very robust and effective method for detecting the shallow three-dimensional velocity structures in an urban area.

Ambient noise wavefield decomposition and shear-wave velocity retrieval in the south of Tehran, Iran and in the Colfiorito basin, Italy

Knowledge about ambient noise wavefield, through decomposition into different participant waves and determination of precise singularities in ellipticity curves, could be very beneficial to improve the obtained subsurface velocity structures. In this paper, two array processing methods; WaveDec and high-resolution Rayleigh three-component beamforming (RTBF) and four single-station methods; Horizontal-to-Vertical Spectral Ratio (HVSR) or H/V, HVTFA (H/V using time frequency analysis), RayDec (Random Decrement Technique) and DOP-E (Degree Of Polarization) are applied to decompose ambient noise wavefield and to extract dispersion and ellipticity curves of surface waves. Two sites with different geological and geophysical conditions are chosen to examine the results of the measurements. The first selected site is in the south of Tehran, Iran, as a site with a low impedance contrast structure. The second one is in the Colfiorito basin, Italy, as a major impedance contrast structure site. Wavefield decomposition resulted in Rayleigh and Love wave separation and also mode distinction. The decomposed modes of surface waves and ellipticity curves of Rayleigh waves are used in a joint inversion process. The retrieved shear-wave velocity profiles are in close agreement with previous studies in the mentioned sites.

Sinkhole detection with 3D full seismic waveform tomography

Sinkhole collapse may result in significant property damage and even loss of life. Early detection of sinkhole attributes (buried voids, raveling zones) is critical to limit the cost of remediation. One of the most promising ways to obtain subsurface imaging is 3D seismic full-waveform inversion. For demonstration, a recently developed 3D Gauss-Newton full-waveform inversion (3D GN-FWI) method is used to detect buried voids, raveling soils, and characterize variable subsurface soil/rock layering. It is based on a finite-difference solution of 3D elastic wave equations and Gauss-Newton optimization. The method is tested first on a data set constructed from the numerical simulation of a challenging synthetic model and subsequently on field data collected from two separate test sites in Florida. For the field tests, receivers and sources are placed in uniform 2D surface grids to acquire the seismic wavefields, which then are inverted to extract the 3D subsurface velocity structures. The inverted synthetic results suggest that the approach is viable for detecting voids and characterizing layering. The field seismic results reveal that the 3D waveform analysis identified a known manmade void (plastic culvert), unknown natural voids, raveling, as well as laterally variable soil/rock layering including rock pinnacles. The results are confirmed later by standard penetration tests, including depth to bedrock, two buried voids, and a raveling soil zone. Our study provides insight into the application of the 3D seismic FWI technique as a powerful tool in detecting shallow voids and other localized subsurface features.

On the differences between horizontal-to-vertical spectral ratios caused by earthquakes and ambient noise—A case study of vertical-array observations in Northern China

Horizontal-to-vertical spectral ratio (HVSR) of both earthquake and ambient noise are applied to perform inversions of the subsurface velocity structure of a certain site, however, their mechanisms differ significantly. A comparative study on earthquake HVSR (EHV) and noise HVSR (NHV) for the same site is rare, hence, herein, a case study is conducted on the Xiangtang vertical-array observation station, located in the city of Tangshan, northern China. The earthquake and noise records observed at this array station are used to investigate the spectral characteristics of EHV and NHV. Differences between the EHV and NHV are compared, and the possible causes for the differences are explored. Further, the mechanisms of the EHV and NHV are discussed, and the suggestions about HVSR inversion are offered. The results of this study show that the EHV and NHV significantly differ in the spectral shape and amplitude in the medium-to-high frequency band. Typically, the amplitude of the observed EHV in the frequency band larger than two times of the predominant frequency (f0) is higher than that of the observed NHV; the mechanisms of the EHV and NHV show clear differences. For frequencies larger than 2f0, the observed EHV is better explained by the theoretical body-wave HVSR, and the observed NHV is better explained by the theoretical full-wave HVSR. Therefore, when performing an inversion to obtain subsurface velocity structure using HVSR, one must carefully distinguish the data sources and choose the correct forward method to obtain more reliable results.

Identification of reservoir characteristics based on seismic velocity tomography using microearthquake method in “X” geothermal field

This research studies microearthquake (MEQ) data analysis in “X” geothermal field. The data was recorded in two phases; before and during drilling, with a total number of 2863 events recorded. The quality assessment method was applied using wadati diagram analysis in order to select the best quality data for advance processing. The result shows that the MEQ events are having trend of Vp/Vs ratio around 1.76. The MEQ data then processed using tomography inversion technique to obtain the image of subsurface seismic velocity structure. The results of seismic tomography inversion shows anomaly pattern which consistent with the lithology and geological structure. The results also shows a good agreement with the interpretation of resistivity section from magnetotelluric (MT) data, as well as from geochemical data. The region infered as reservoir zone are observed with dominantly high Vp/Vs anomaly which indicates incompressible characteristics that may relate to existence of water content, while the low Vp/Vs anomaly which indicates steam cap existence are not observed in this region. This pattern of anomaly leads to inference that the reservoir characteristics are dominated with liquid phase of fluid.

SCoda and Rayleigh Waves From a Decade of Repeating Earthquakes Reveal Discordant Temporal Velocity Changes Since the 2004 Sumatra Earthquake

AbstractTemporal changes in the subsurface seismic velocity structure reflect the physical processes that modulate the properties of the media through which seismic waves propagate. These processes, such as healing of the surface damage zone and deep crustal deformation, are described by similar functional forms and operate on similar timescales, making it difficult to determine which process drives the observed changes. We examine earthquake‐induced velocity changes using the measured lag‐time time seriesτ(t) of the repeating earthquake sequences since the 2004Mw9.2 Sumatra and 2005Mw8.6 Nias earthquakes. TheScoda velocity changes (δVS, equivalent to−τS) recover steadily during the 2005–2015 period. The Rayleigh wave velocity changes (δVLR, or−τLR) undergo transient recovery, followed by a strongδVLRreduction in late 2007.δVSrecovery is most likely driven by deep processes, whereas the temporal breaks inδVLRrecovery in 2007 mostly reflect surface damage and healing induced by the strong ground motions of the 2004Mw9.2, 2005Mw8.6, 2007Mw8.4, andMw7.9 Bengkulu and 2008Mw7.3 Simeulue earthquakes. The observed differences between the temporal variations inδVSandδVLRcan distinguish deep processes from healing of the surface damage zone.

Velocity model estimation by means of Full Waveform Inversion of transmitted waves: An example from a seismic profile in the geothermal areas of Southern Tuscany, Italy

We propose an FWI strategy that makes use of transmitted waves as input data and utilizes both global and local optimization methods to estimate the P-wave velocity model of the subsurface. We envisage that our approach may be applicable to difficult seismic land data, like those from geothermal areas characterised by complex geological structures. As a test case, we considered the CROP/18A seismic reflection profile that crosses the geothermal field of Larderello (southern Tuscany, Italy). The aim is to estimate the P-wave velocity model down to a few kilometres depth below the surface that could be used as complementary information to the standard seismic reflection image which, in this case, does not show interpretable reflections in a range of depths accessible to industrial drillings.One innovative aspect of the inversion we propose with respect to conventional FWI approaches is its independence of a starting model that, ideally, should reproduce the true long wavelength velocity structure of the subsurface and that may be rather difficult to obtain in case of low quality data and complex geology. We lessen the dependence on knowledge of a suitable starting model by performing a sequence of two inversions. First, we employ a genetic-algorithm (GA) based inversion, a global optimisation method that does not require any specific starting model, resulting in a long wavelength, low-resolution velocity model. This model then becomes the starting model for a second FWI, driven by a local optimisation algorithm, aimed at bringing in the fine details of the subsurface velocity structure. The reliability of the final model is checked by comparing observed and predicted waves for many common shot gathers along the seismic line and through the matching between the velocities measured by check shots in two nearby wells and the FWI velocities in the same locations. Many details of the velocity field, likely related to metamorphic and igneous formations, become apparent and may complement the interpretation of the standard reflection image. From these results, it appears that the use of transmitted waves and of the FWI approach discussed here may effectively improve the information for geophysical interpretation of challenging seismic land data, such as those that characterise many areas of geothermal exploration.

Ambient Vibration Analysis on Large Scale Arrays When Lateral Variations Occur in the Subsurface: A Study Case in Switzerland

The ambient vibration analysis is a non-invasive and low-cost technique used in site characterization studies to reconstruct the subsurface velocity structure. Depending on the goal of the research, the investigated depth ranges from tens to hundreds of meters. In this work, we aimed at investigating the deeper contrasts within the crust and in particular down to the sedimentary-rock basement transition located at thousands of meters of depth. To achieve this goal, three seismic arrays with minimum and maximum interstation distances of 7.9 m and 26.8 km were deployed around the village of Schafisheim. Schafisheim is located in the Swiss Molasse Basin, a sedimentary basin stretching from Lake Constance to Lake Geneva with a thickness ranging from 800 to 900 m in the north to 5 km in the south. To compute the multimodal dispersion curves for Rayleigh and Love waves and the Rayleigh wave ellipticity angles, the data were processed using two single-station and three array processing techniques. A preliminary analysis of the inversion results pointed out a good agreement with the fundamental modes of Rayleigh and Love waves used in the inversion and a quite strong disagreement with the higher modes. The impossibility to explain at the same time most of the dispersion curves was interpreted as the co-existence, within the investigated area, of portions of the subsurface with different geophysical properties. The hypothesis was confirmed by the Horizontal-to-Vertical spectral analysis (H/V) which indicated the presence of two distinguished areas. The observation allowed a new interpretation and the identification of the Rayleigh and Love wave fundamental modes and of the S-wave velocity profiles to be reconstructed for each investigated zone. It results in two S-wave velocity profiles with similar velocities down to 15 km deferring only in their shallow portions due to the occurrence of a low velocity zone at a depth of 50–150 m at the centre of the investigated area.

The deep subsurface structure beneath Lusi and the adjacent volcanic chain inferred from local travel-time tomography

Since 2006, the Lusi mud eruption is ongoing in the northeast of Java Island, Indonesia. We perform a 3D seismic tomography of the area surrounding Lusi and the adjacent Arjuno-Welirang volcanic complex using P-wave arrivals of 797 local earthquakes recorded on 31 seismic stations for two years between 2015 and 2016. Resolution tests show that we can reasonably well constrain the subsurface velocity structure beneath Lusi down to a depth of 20 km. The results of our tomography reveal that Lusi has a shallow plumbing system that is interconnected with the neighboring Penanggungan volcano. Their joint system reaches down to a depth of 20 km. The Arjuno-Welirang volcanic complex has its own deep-rooted plumbing system that might be connected to Penanggungan volcano and Lusi along the Watukosek Fault zone. The connection of Lusi with this large reservoir volume elucidates the longevity of the Lusi eruption.

Crustal Structure Beneath KLNI Station in Lombok Island Based on Teleseismic Receiver Function

Lombok Island is part of the Sunda-Banda arc system which is one of the most active tectonic regions in the world. In Lombok area there are two subduction zones namely Java trench in the south and Flores back arc thrust in the north. Information on crustal structure is needed to explain the tectonic development and geodynamic interaction of the area. The receiver function is one of the techniques used in the study to identify subsurface and wave velocity structures based on earthquake data regardless of the source mechanism. We analyzed receiver function of teleseismic events (epicenter distance 30-90°) and magnitude more than 6 (M>6) to estimate crustal thickness, Vp/Vs ratio and wave velocity structure beneath KLNI station using H-κ stacking and inversion method. We obtained that the crustal thickness and Vp/Vs ratio beneath KLNI station is about 21.38 km and 2.05. The high Vp/Vs ratio might be an indication of partial melting related to the upwelling of hot asthenosphere material through the subduction slab of Java Trench. The depth of the low velocity zone beneath KLNI station is at 10-18 km which is estimated due to the presence of geothermal

The effect of body waves on phase-velocity determined by the spatial autocorrelation (SPAC) method, evaluated using full-wave modelling

Body waves may affect phase velocity obtained from microtremor array surveys in some rare cases. Fitting theoretical phase velocities based on a surface-wave theory to observed phase velocities affected by body waves would therefore result in distorted images of subsurface S-wave velocity structure. In this study, we present a method for the theoretical calculation of phase velocities in which the full-wave field (i.e. a wavefield including not only surface waves but also body waves) is taken into account. In numerical experiments conducted in this study, in which we considered the full-wave field, we generated synthetic microtremors by randomly distributing point vibration sources on the surface of a horizontally stratified velocity model. We then determined the phase velocities by applying the spatial autocorrelation (SPAC) method to the synthetic vertical-component wave data. The phase-velocity dispersion curve thus obtained exhibited a shape with a clear peak, with a peak value (peak phase velocity) exceeding the S-wave velocity of a bedrock in the model, which was not explainable with a surface-wave (Rayleigh-wave) theory. We conducted systematic numerical experiments and clarified the following two features of the peak phase velocity: (1) the peak phase velocity becomes large as the contrast of the S-wave velocities between the surface layer and the bedrock, or the P-to-S-wave velocity ratio (related to the Poisson’s ratio) in the surface layer gets large, and (2) the frequency at which peak phase velocity occurs (peak frequency) lies in the vicinity of the S-wave resonance frequency of the ground. Both the peak phase velocity and the peak frequency were theoretically reproduced by the calculation method that we propose in this study, based on a SPAC method modified to consider the full-wave field. These results imply the possible improvement in the accuracy of microtremor array survey analysis for velocity-structure inference, by applying a full-wave theory to the peak phase velocity.

Sub-surface structures and site effects extracted from ambient noise in metropolitan Guangzhou, China

In this paper, ambient noise is used to investigate near-surface structures and site effects in metropolitan Guangzhou. We deployed 94 short period stations across Guangzhou area, and later a dense linear array across the Shougouling Fault (SF). Using more than one month's continuous data, we invert three-dimensional shear wave velocity structures of Guangzhou area via ambient noise tomography. Results show low velocity near the Guangzhou-Conghua Fault (GCF), Shougouling Fault (SF), and beneath the alluvial plain in the south region. Meanwhile, high-velocity anomalies are beneath the mountain area in northeastern region. Moreover, we obtain the sediment thickness and sub-surface shear wave velocity structures around the Shougouling Fault (SF) by horizontal to vertical spectral ratio (HVSR) method. HVSR results show a significant shift in the thickness of the sedimentary layer across the SF. Shear wave velocity derived from HVSR curves has a consistent trend of variation with the sediment thickness. Our results provide a better understanding of sub-surface structures of metropolitan Guangzhou, and can be served as a reference model for geological disaster migration prediction in the city.

InversionNet: An Efficient and Accurate Data-Driven Full Waveform Inversion

Full-waveform inversion problems are usually formulated as optimization problems, where the forward-wave propagation operator $f$ maps the subsurface velocity structures to seismic signals. The existing computational methods for solving full-waveform inversion are not only computationally expensive, but also yields low-resolution results because of the ill-posedness and cycle skipping issues of full-waveform inversion. To resolve those issues, we employ machine-learning techniques to solve the full-waveform inversion. Specifically, we focus on applying the convolutional neural network~(CNN) to directly derive the inversion operator $f^{-1}$ so that the velocity structure can be obtained without knowing the forward operator $f$. We build a convolutional neural network with an encoder-decoder structure to model the correspondence from seismic data to subsurface velocity structures. Furthermore, we employ the conditional random field~(CRF) on top of the CNN to generate structural predictions by modeling the interactions between different locations on the velocity model. Our numerical examples using synthetic seismic reflection data show that the propose CNN-CRF model significantly improve the accuracy of the velocity inversion while the computational time is reduced.

What can P-wave polarization direction data tell us about subsurface velocity structures?

SUMMARY As seismic waves propagate in the Earth, the directions of particle motions are affected by the media that they encounter, and thus seismic wave polarization direction carries the information on the media. So far there remains unclear about what can be inferred from the P-wave polarization direction data. For clarifying it, we discuss the mapping relation between polarization direction and velocity distribution. It is found that the velocity model cannot be derived uniquely from the polarization direction data. By analysing the relation between slowness vectors of the seismic ray at the source and the receiver, we find that relative velocity gradient is the physical quantity that describes the capability to deflect seismic rays in a continuous medium. The equation describing the relation between polarization direction and relative velocity gradient is given. For imaging relative velocity gradients, we derive the calculation formula for the partial derivative of polarization direction with respect to velocity gradient parameters. Synthetic experiments are conducted. The test results demonstrate that the absolute velocity model cannot be recovered from P-wave polarization direction data, but the relative velocity gradient model can. Polarization direction tomography gives a way to build gradient maps for the geometric characteristic of the subsurface velocity structures.

Modeling of Subsurface Velocity Structures from Seismic Bedrock to Ground Surface in the Tokai Region, Japan, for Broadband Strong Ground Motion Prediction

For sophistication of strong ground motion prediction in terms of disaster mitigation, one of the principal issues is to model subsurface velocity structures so that characteristics of earthquake ground motions can be reproduced in the broadband range 0.1 Hz to 10 Hz. In recent years, subsurface structures have been modeled in sedimentary layers on seismic bedrock for a few regions of Japan, in a national project. In this study, subsurface velocity structures were modeled from seismic bedrock to the ground surface for the Tokai region. These models were constructed in accordance with the subsurface velocity structure modeling scheme published by the Headquarters for Earthquake Research Promotion. To begin with, initial models were constructed based on existing bore-hole data, geological information, etc. Next, they were improved based on results of microtremor explorations which had been conducted in recent years. It was found that the new model had different characteristics to the conventional model. This paper will present the modeling process and characteristics of distribution maps for velocity structures and amplification index.

Passive Adjoint Tomography of the Crustal and Upper Mantle Beneath Eastern Tibet With a W2‐Norm Misfit Function

AbstractFull waveform tomography is an effective method to obtain high‐resolution subsurface velocity structures. It is, however, difficult for the conventional L2‐norm‐based inversion to reach the global minimum because of cycle‐skipping problem. In this study, we investigated the feasibility of using the quadratic Wasserstein metric distance (W2 norm) as the misfit function for passive adjoint tomography. We first derived equations of the Fréchet gradient and adjoint source under W2‐norm misfit function. We then conducted numerical experiments to illustrate the effectiveness of the proposed method in avoiding cycle skipping. We finally applied the adjoint tomography to eastern Tibet and obtained a 3‐D velocity model of the lithosphere. The adjoint tomography revealed two low‐velocity channels beneath the NE and SE margins of the Tibetan plateau, which were also observed by previous studies. Our results are consistent with the lower crust flow model that was proposed to explain the deformation occurring at the two margins.

Comparing performance of multi-frequency bands Occam’s receiver function inversion to standard linearized receiver function inversion

Receiver function (RF) is a characterized waveform that sensitive to the sub-surface velocity structure near the seismic station. In the geophysical problem, an inversion scheme has been used for determining the structure model that corresponded to the observed data. Because RF inversion is a non-linear problem, the standard method such as linearized inversion suffered in the non-uniqueness of the results. Thus, using a stochastic optimization algorithm (e.g. Monte Carlo or Generic Algorithm) is the most interesting trend for solving the problem. However, their computation cost is high and the efficiency is strongly depended on user setting. In this work, Occam’s inversion algorithm, which is popular in the electromagnetic method, and multi-frequency bands receiver function have been integrated to improve the efficiency of receiver function inversion. The propagation method has been applied for the calculation of seismogram in this work. The synthetic cases, which are thin low velocity zone model (TLVZ), broad low velocity zone model (BLVZ), thin high velocity zone model (THVZ) and broad high velocity zone model (BHVZ), were used to compare the performance of this newly implemented algorithm and the standard linearized inversion. Single frequency band (SFO) and multi-frequency bands Occam’s inversion (MFO) have been tested for each case. For, TLVZ, THVZ and BHVZ, the MFO inversion provided better fitting results than SFO in most cases. Comparing to linearized inversion of CPS330, the MFO results still provide better performance. For the smooth low velocity structure in BLVZ model, both linearized and MFO can recover a significant structure of the true model. More concisely, the linearized inversion can recover the absolute velocity of the true model better than MFO. By the way, the MFO provides the best data fitting for the BLVZ. In summary, the MFO receiver function inversion that implemented in this study provides a new improved tool for a seismologist to investigating the subsurface structure.

Ray‐based reflection traveltime tomography using approximate stationary points

ABSTRACTReflection traveltime tomography has been used to describe subsurface velocity structures which, in practice, can be used as a background or initial model for pre‐stack depth migration or full waveform inversion. Conventional reflection traveltime tomography is performed by solving an optimization problem based on a ray‐tracing method. As a result, reflection traveltime tomography requires heavy computational efforts to carry out ray tracing and solve a large matrix equation. In addition, like most data‐domain tomography methods, reflection traveltime tomography depends on initial guesses and suffers from non‐uniqueness and uncertainty of solutions. In this research, we propose a deterministic ray‐based reflection traveltime tomography method by applying seismic interferometry. This method does not suffer from the non‐uniqueness problem and does not require a priori information on subsurface media. By adding a virtual layer (whose properties are known) on top of the real surface and applying convolution‐type interferometry, we approximately determine the stationary points (i.e., incident raypaths in the virtual layer). Then, we generate reflection points for a range of assumed velocities and estimate the velocity by considering the number of reflection points and the traveltime difference between the observed and calculated data. The reflection surface can then be recovered by using the estimated velocity. Once the first target layer is resolved, we can recover the whole media by recursively applying the same method to the lower layers. Numerical examples using surface seismic profile data for homogeneous‐layer (with a low‐velocity layer) and inhomogeneous‐layer models and real field data experiments on the Congo data set demonstrate that our method can successfully recover the velocities and depths of subsurface media without initial guesses. However, our method has some limitations for multi‐layer models because the method does not have sufficient reflection points for the deeper layers.