Advances in controlled‐source seismic imaging

We have witnessed a transit of seismic tomography from ray-based traveltime inversion to wave equation-based waveform inversion during the past two decades. As is widely known, full waveform inversion outperforms traveltime tomography in resolving velocity variations with dimensions comparable and smaller than the dominant wavelength. However, this may not be always true in real practice, mainly due to the presence of unknown data noise, the lack of accurate initial models for the iterative inversion process (including material properties and source mechanisms), and the high demand for computational resources. Further efforts are required to develop techniques for using full waveform contents in seismic imaging studies, especially for the use of high-frequency waveform data (>1 Hz on regional scales). We have attempted to use common-source double-difference traveltime data in wave equation-based adjoint tomography studies of subsurface structures beneath Northeast Japan and Alaska. Because of the high level of waveform similarity, reliable common-source double-difference traveltime data at neighboring seismic stations are measured via cross-correlation approach. Insightful results are obtained. However, it is still computationally prohibitive to model high-frequency data by solving 3-D wave equations, limiting the resolution of seismic images. We admit that there is still a long way to go for the possible wide application of high-frequency full waveform inversion.Traveltime is the most reliable information that can be extracted from raw seismological recordings.  We have developed a new modality of seismic imaging, called adjoint-state traveltime tomography, to unleash the full potential of traveltime in imaging subsurface structures. To avoid potential failure of ray tracing in 3D complex media, isotropic eikonal equation and anisotropic eikonal equations are used to model seismic traveltime field in heterogenous and anisotropic media, and the associated inverse problems are solved by the efficient adjoint state method. No ray tracing is required for the novel adjoint-state traveltime tomography method. Importantly, the sensitivities of traveltime-related objective functions to material parameters can be accurately measured even in complex media. We view adjoint-state traveltime tomography as a new framework for seismic tomography, as it naturally implements various body wave and surface wave tomographic inversions in a very similar way. Good performances of the adjoint-state traveltime tomography method will be reported via various case studies in regions with typically different tectonic settings. It is worth noting that an accompanying software package, TomoATT, is under development.

Velocity anisotropy caused by fine layering has a strong effect on the travel paths at large vertical offsets in crosshole seismic data. This effect must be considered in seismic waveform tomography, so that the inversion procedure can match the waveforms in real data. In this paper, we propose to properly take into account this anisotropic effect in two stages. First, we invert for the anisotropy parameter simultaneously with the (horizontal) velocity in travel time tomography. Then we incorporate this estimated anisotropy parameter model in the waveform simulation during waveform tomography, to refine the velocity image. When considering this anisotropic effect, far-offset data are properly used in the inversion, then both travel time and waveform tomography may have a better ray coverage, producing a velocity image with a higher resolution.

Seismic travel time tomography using surface waves is an effective tool for three-dimensional crustal imaging. Historically, these surface waves are the result of active seismic sources or earthquakes. More recently, however, surface waves retrieved through the application of seismic interferometry have also been exploited. Conventionally, two-step inversion algorithms are employed to solve the tomographic inverse problem. That is, a first inversion results in frequency-dependent, two-dimensional maps of phase velocity, which then serve as input for a series of independent, one-dimensional frequency-to-depth inversions. As such, a set of localized depth-dependent velocity profiles are obtained at the surface points. Stitching these separate profiles together subsequently yields a three-dimensional velocity model. Relatively recently, a one-step three-dimensional non-linear tomographic algorithm has been proposed. The algorithm is rooted in a Bayesian framework using Markov chains with reversible jumps, and is referred to as transdimensional tomography. Specifically, the three-dimensional velocity field is parameterized by means of a polyhedral Voronoi tessellation. In this study, we investigate the potential of this algorithm for the purpose of recovering the three-dimensional surface-wave-velocity structure from ambient noise recorded on and around the Reykjanes Peninsula, southwest Iceland. To that end, we design a number of synthetic tests that take into account the station configuration of the Reykjanes seismic network. We find that the algorithm is able to recover the 3D velocity structure at various scales in areas where station density is high. In addition, we find that the standard deviation of the recovered velocities is low in those regions. At the same time, the velocity structure is less well recovered in parts of the peninsula sampled by fewer stations. This implies that the algorithm successfully adapts model resolution to the density of rays. It also adapts model resolution to the amount of noise in the travel times. Because the algorithm is computationally demanding, we modify the algorithm such that computational costs are reduced while sufficiently preserving non-linearity. We conclude that the algorithm can now be applied adequately to travel times extracted from station–station cross correlations by the Reykjanes seismic network.

Seismic attenuation tomography and its application to rock mass evaluation

Several types of seismological data, including surface wave group and phase velocities, travel times from large explosions, and teleseismic travel time anomalies, have indicated that there are significant regional variations in the upper few hundred kilometers of the mantle beneath continental areas. Body wave travel times and amplitudes from large chemical and nuclear explosions are used in this study to delineate the details of these variations beneath North America. As a preliminary step in this study, theoretical P wave travel times, apparent velocities, and amplitudes have been calculated for a number of proposed upper mantle models, those of Gutenberg, Jeffreys, Lehman, and Lukk and Nersesov. These quantities have been calculated for both P and S waves for model CIT11GB, which is derived from surface wave dispersion data. First arrival times for all the models except that of Lukk and Nersesov are in close agreement, but the travel time curves for later arrivals are both qualitatively and quantitatively very different. For model CIT11GB, there are two large, overlapping regions of triplication of the travel time curve, produced by regions of rapid velocity increase near depths of 400 and 600 km. Throughout the distance range from 10 to 40 degrees, the later arrivals produced by these discontinuities have larger amplitudes than the first arrivals. The amplitudes of body waves, in fact, are extremely sensitive to small variations in the velocity structure, and provide a powerful tool for studying structural details. Most of eastern North America, including the Canadian Shield has a Pn velocity of about 8.1 km/sec, with a nearly abrupt increase in compressional velocity by ~ 0.3 km/sec near at a depth varying regionally between 60 and 90 km. Variations in the structure of this part of the mantle are significant even within the Canadian Shield. The low-velocity zone is a minor feature in eastern North America and is subject to pronounced regional variations. It is 30 to 50 km thick, and occurs somewhere in the depth range from 80 to 160 km. The velocity decrease is less than 0.2 km/sec. Consideration of the absolute amplitudes indicates that the attenuation due to anelasticity is negligible for 2 hz waves in the upper 200 km along the southeastern and southwestern margins of the Canadian Shield. For compressional waves the average Q for this region is > 3000. The amplitudes also indicate that the velocity gradient is at least 2 x 10-3 both above and below the low-velocity zone, implying that the temperature gradient is < 4.8°C/km if the regions are chemically homogeneous. In western North America, the low-velocity zone is a pronounced feature, extending to the base of the crust and having minimum velocities of 7.7 to 7.8 km/sec. Beneath the Colorado Plateau and Southern Rocky Mountains provinces, there is a rapid velocity increase of about 0.3 km/sec, similar to that observed in eastern North America, but near a depth of 100 km. Complicated travel time curves observed on profiles with stations in both eastern and western North America can be explained in detail by a model taking into account the lateral variations in the structure of the low-velocity zone. These variations involve primarily the velocity within the zone and the depth to the top of the zone; the depth to the bottom is, for both regions, between 140 and 160 km. The depth to the transition zone near 400 km also varies regionally, by about 30-40 km. These differences imply variations of 250 °C in the temperature or 6 % in the iron content of the mantle, if the phase transformation of olivine to the spinel structure is assumed responsible. The structural variations at this depth are not correlated with those at shallower depths, and follow no obvious simple pattern. The computer programs used in this study are described in the Appendices. The program TTINV (Appendix IV) fits spherically symmetric earth models to observed travel time data. The method, described in Appendix III, resembles conventional least-square fitting, using partial derivatives of the travel time with respect to the model parameters to perturb an initial model. The usual ill-conditioned nature of least-squares techniques is avoided by a technique which minimizes both the travel time residuals and the model perturbations. Spherically symmetric earth models, however, have been found inadequate to explain most of the observed travel times in this study. TVT4, a computer program that performs ray theory calculations for a laterally inhomogeneous earth model, is described in Appendix II. Appendix I gives a derivation of seismic ray theory for an arbitrarily inhomogeneous earth model.

This thesis addresses the fundamental problem of determining the radial and lateral structure of the earth's interior using body wave observations. My approach is a cautious one, involving detailed analysis of a substantial data set in which I attempt to isolate the contributions from radial and lateral structure in both the upper and lower mantles. This is an elusive task, as I am concerned with the 2 or 3% fluctuations about standard earth models, which produce rather subtle effects on teleseismic signals. However, this is the level of precision to which the three dimensional structure of the earth must be determined if we are to map the dynamical and compositional configuration of the earth. Studies similar to those described here cannot be conducted on a detailed global basis, given the intrinsic limitations due to station and earthquake distribution, so I have repeatedly emphasized the qualitative implications of my results, as they probably provide a representative sampling of the subtle, but significant heterogeneity of the earth. The topics addressed in the following work appear unrelated at first glance, ranging from determination of the detailed shear velocity structure at the base the mantle to variations in attenuation and velocity structure of the upper mantle. However, the data set used throughout is largely the same, and this in itself indicates the need to consider all of the complexity discussed herein for future progress to be made in mapping the three dimensional structure of the earth with a high level of confidence. The first chapter of the thesis presents results of a waveform analysis of transversely polarized SH signals that propagate through the lowermost mantle. The waveforms of these phases show clear interference patterns due to interaction with a discontinuous shear velocity increase about 280 km above the core-mantle boundary. This discontinuity, which has not been detected previously, is manifested in signals sampling three widely separated portions of the lower mantle, and hence is a good candidate for a global radial earth structure. Very detailed inspection of the signals reveals evidence for lateral variations in the depth of the discontinuity, which provides a procedure by which to map the structure of the D region in detail. Analysis of the relative amplitudes of SH and ScSH signals reveals that the velocity gradient above the core boundary is consistent with the smoothly varying gradients in most gross earth models, but evidence for local high velocity gradients at the base of the mantle is detected in ScSV signals. Chapter II presents a travel time analysis of the same data set used in Chapter I, which was motivated by the observation of large amplitude and travel time anomalies in the S and ScS data. An emphasis is placed on isolating the portions of the S and ScS paths which are anomalous. A strong empirical case is made for the existance of localized regions with scale lengths of 1000 to 2000 km and 2% velocity anomalies within the lower mantle at depths from 1000 to 2700 km. The long period signals traversing these regions show as much as a factor of 2 amplitude enhancement or diminution. This result demonstrates that both amplitude and travel time anomalies are induced by lateral structure in the portion of the mantle assumed to be homogeneous in most seismological studies. The third chapter is an analysis of the influence of upper mantle variations in attenuation, velocity structure and receiver structure on the S and ScS signals analyzed in the lower mantle studies. These variations contaminate and complicate the interpretation of any data set used to study deeper earth structure. Along with evidence for very strong and abrupt variations in upper mantle properties, results are presented which indicate the inadequacy of assumptions that are frequently made about the nature of long period body-wave receiver functions.

plex shallow subsurface. This workflow incorporates dispersion curve inversion, traveltime tomography and full waveform inversion into a single framework to optimize the advantages of each. An accurate near-surface velocity is derived using the dispersion curve inversion technique that overcomes limitations of traditional traveltime tomography. This velocity is then incorporated into layer-stripping full waveform or RMO inversion schemes by “freezing” the ultra-shallow part of the model in order to retrieve the deeper velocity structure. We illustrate the workflow on a test data set from Inner Mongolia, China, whose long-offsets and low-frequency content provide the necessary input for both the dispersion curve and full waveform inversion techniques. Using this data set, we demonstrate the efficacy of the integrated methodology in determining both the macro-velocity model and near-surface velocities in structurally complex geological setting.

One of the most important studies of the earth sciences is that of the Earth's interior structure. There are many sources of data for Earth tomography models: flrst-arrival passive seismic data (from the actual earthquakes), flrst-arrival active seismic data (from the seismic experiments), gravity data, and surface waves. Currently, each of these datasets is processed separately, resulting in several difierent Earth models that have speciflc coverage areas, difierent spatial resolutions and varying degrees of accuracy. These models often provide complimentary geophysical information on earth structure (P and S wave velocity structure). Combining the information derived from each requires a joint inversion approach. Designing such joint inversion techniques presents an important theoretical and practical challenge. While such joint inversion methods are being developed, as a flrst step, we propose a practical solution: to fuse the Earth models coming from difierent datasets. Since these Earth models have difierent areas of coverage, model fusion is especially important because some of the resulting models provide better accuracy and/or spatial resolution in some spatial areas and in some depths while other models provide a better accuracy and/or spatial resolution in other areas or depths. The models used in this paper contain measurements that have not only difierent accuracy and coverage, but also difierent spatial resolution. We describe how to fuse such models under interval and probabilistic uncertainty. The resulting techniques can be used in other situations when we need to merge models of difierent accuracy and spatial resolution.

One of the most important studies of the earth sciences is that of the Earth’s interior structure. There are many sources of data for the construction of tomographic models of earth structure. These include the time of the first arriving waves from earthquakes and from man-made sources, measurements of the earth’s gravity field and measurements of the dispersion of surface waves generated from earthquakes. Formally integrating the information derived from multiple types of data sources is an important theoretical and practical challenge. While such combination methods are being developed, as a first step, we propose a practical solution: to fuse the Earth models coming from different data sets. The models used in this paper contain measurements that have not only different accuracy and coverage but also different spatial resolution. We describe how to fuse such models under interval and probabilistic uncertainty. The resulting techniques can be used in other situations when we need to merge models of different accuracy and spatial resolution.

Improvement of seismic tomographic results can be expected if we use not only the first P-arrival times but also incorporate later arrivals with more complicated ray paths. The later arrivals, which are generated by the complexity of the earth structure, provide more information on the media. However, traveltime tomography is a non-linear inverse problem and generally involves tens of thousands of arrival times. Sophisticated techniques for use in complex velocity structures are generally too time-consumming for our applications. An efficient technique for rapidly computing the traveltimes for ray paths through a realistic 2-D/3-D earth model is required.

Body waves traversing the Earth's interior from a seismic source to receivers on the surface carry rich information about its internal structures. Their travel time measurements have been widely used in seismology to constrain Earth's interior at the global scale by mapping the time anomaly along their ray paths. However, picking the travel time of global seismic waves, suitable for studying Earth's fine‐scale structures, requires highly skilled personnel and is often fairly subjective. Here, we report the development of an automatic picker for PKIKP waves traveling through the inner core (IC), especially nearly along Earth's diameters, based on the latest advances in supervised deep learning. A convolutional neural network (CNN) we develop automatically determines the PKIKP onset on vertical seismograms near its theoretical prediction of cataloged earthquakes. As high‐quality manual onset picks of global seismic phases are limited, we employ a scheme to generate a synthetic supervised training data set containing 300,000 waveforms. The PKIKP onsets picked by our trained CNN automatic picker exhibit a mean absolute error of ∼0.5 s compared to 1,503 manual picks, comparable to the estimated human‐picking error. In an integration test, the automatic picks obtained from an extended waveform data set yield a cylindrically anisotropic IC model that agrees well with the models inferred from manual picks, which illustrates the success of this pilot model. This is a significant step closer to harvesting an unprecedented volume of travel time measurements for studying the IC or other regions of the Earth's deep interior.

Detecting and imaging dense non-aqueous phase liquids (DNAPLs) in the subsurface is a challenging problem that is of interest to the Department of Energy, Department of Defense, and several local and state agencies. Current DNAPL detection techniques, such as Direct Push Probe Technologies (DPT) and In-Situ Tracers (IST) have risks and limitations. Cross-well radar (CWR) is a radar-based geophysical technique with low invasiveness for real time monitoring of DNAPLs. This technique uses electromagnetic waves transmitted and received through antennas in the subsurface. The computerized tomography is applied to a pilot-scale experimental facility constructed by the authors (referred as SoilBED) in this research. Cross-tomography data at multiple depths and locations were collected to simulate a 1/100 scaled contaminated soil problem. The experimental results are converted to the time domain signals. The resulting signal travel times are compared with the simulated results by FDTD. The transformed signals are used to prepare signal intensity and travel time tomograms of the SoilBED cross-sectional slices, for background and scattered field to study the signature of the scatterers. The results are able to visualize the dielectric objects in the saturated soil. Introduction Cross-well radar is a potential technique for DNAPL detection. Most of the research and field works done in Cross-well radar and other radar based geophysical techniques deal with a single frequency. Selecting one single frequency may be effective for one non-dispersive medium, but perhaps not for other non-dispersive or most dispersive media. This research covers a broadband of frequency from 400 MHz to 2.2 GHz. Having broadband frequency data helps to either use the entire frequency range and work in frequency domain or select a single frequency or narrower frequency range that seems to have enough detection potential out of the entire observed frequency response. It is important to consider cost and practicality of the system. Establishing practical and cost effective systems require strong theoretical and laboratory background conditions. Therefore, both theoretical and experimental aspects of the research are studied in parallel, to evaluate the feasibility of DNAPL detection using CWR and to study wave propagation behavior and forward problem. Collecting repeatable and reproducible frequency response data is a very challenging process. This difficultly achieved data can be manipulated first or directly used towards image processing. As any research subject, conducting both experimental and theoretical approaches is the most desirable way to advance fundamental knowledge in the field of interest. Frequency response and arrival time can be converted to each other. Thus, frequency domain results can be transformed to time domain and used towards techniques such as travel time tomography and vice versa. The experimental frequency domain results can be converted to CTF (Channel Transfer Functions) and compared to the simulated frequency domain results by the first order Born approximation (Farid 2004, Farid et al. 2005, and Zhan 2005). The FDTD simulation code (Farid 2002 and 2003) is available as the theoretical tool to perform the comparison in time domain and study travel time. If the experimental results are converted to time domain and transformed to travel time, they can be compared to the FDTD simulated results. Travel Time Computed by FDTD To establish a better knowledge of wave propagation in soil, both dry and saturated backgrounds are studied. The sandy soil in the SoilBED facility with 3.88% moisture content is referred to as dry soil.

Seismic bulletins, with trustworthy phase picks, origin times, and source locations are key for regional seismic studies, such as travel-time (TT) tomography, attenuation tomography, and anisotropy studies. To lay the groundwork for such studies in Israel, we revised the seismic bulletin of Israel and the surrounding area and obtained a trustworthy TT data set. From the earthquake and explosion bulletins of the Geophysical Institute of Israel, we compiled a starting data set of about 123,000 earthquakes and explosions that occurred during the past 40 yr. After screening out the poorly recorded events, we were left with a data set of ∼38,000 well-recorded events. We then revised the remaining data set in two consecutive steps. In the first, we reviewed and updated station metadata, including changes in station metadata parameters over time. In the second step, we jointly relocated a list of selected seismic events, using the Bayesian hierarchical location software package (BayesLoc) of Myers et al. (2007) that performs joint relocation of multiple events. We observed striking dissimilarities between the spatial distributions of the newly relocated catalog and the initial locations. Although the depth distribution of the starting catalog is trimodal with peaks at 0, 5, and 10 km, the distribution in this study is unimodal, with a broad peak between 7.5 and 12.5 km. By differencing the observed arrival times and the origin times obtained through relocation with BayesLoc, we obtained a revised TT database that consists of 261,336 Pg, 132,876 Pn, 114,816 Sg, and 60,394 Sn arrivals, from a set of 30,458 jointly relocated seismic sources. We compared prerevision and postrevision TTs as a function of epicentral distance and concluded that the revised data set contains far fewer outliers and inconsistencies than the original data set. The revised TT data set may be used for seismic studies, such as TT tomography, attenuation tomography, and anisotropy studies.

This paper describes the use of computer algebra in the derivation of analytical expressions for travel times of rays at finite offsets for a 3D earth model comprising a few layers where layer thickness, velocity and dip angles are represented symbolically. Expressions for offset and travel times are derived symbolically in terms of the initial ray directions and parameters of the earth model. These expressions and the derivatives of the offsets with respect to ray direction are converted to C program subroutines. Numerical techniques are then used to find travel times for given values of offset, limited only by computer precision. Since derivatives of offsets with respect to initial ray directions can also be calculated directly, a modified Newton’s method was used to give rapid convergence to initial ray directions required for a given value of offset. The travel time is then calculated from the ray directions so derived. This is an alternative to previously used techniques based either on approximations or on purely numerical calculations. Numerical values of model parameters may be readily changed and the calculations repeated for the changed model. Analytical solutions have been derived for a number of multiple and inter-bed ray paths as well as for primary events. Travel time calculations are fast and accurate since they are based on analytical expressions. Shot records, CMP gathers or other geometry can be readily defined. This paper discusses the potential use of these travel time calculations for improving numerical modelling, identifying multiple reflection events, identifying anisotropic effects and determining the variation of reflection point locations in CMP gathers. The potential for, and possible limitation of, extensions to include multi-layer earth models, converted waves and calculation of amplitudes is also discussed.

Despite proven analogues elsewhere in Malaysia, pre-carbonate hydrocarbon targets offshore Sarawak have met with limited drilling success. This has been mainly attributed to the poor seismic imaging beneath complex velocity carbonate reefs and platforms with low signal-to-noise ratio and inadequate illumination. Seismic imaging is challenging due to inconsistent illumination, multiple energy contamination and complex wave propagation. A collaborative study was undertaken to model these challenges and recommend optimum acquisition and signal processing solutions to produce a step change in seismic imaging of the pre-carbonate targets. The study involved building a detailed anisotropic 3D earth model using borehole and seismic data covering more than 600 km2 area to create a robust synthetic seismic dataset. Geological and petrophysical review preceded seismic interpretation and simultaneous AVO inversion results were used to populate a fine scale model. Testing and validation of the model was performed on a data subset to ensured the kinematics and amplitudes were consistent with the actual seismic response. The model was then used to create over 15 billion prestack traces using an anisotropic finite difference wave propagation modelling algorithm. From this data, different survey designs were extracted including; narrow-azimuth (NAZ), multi-azimuth (MAZ), wide-azimuth (WAZ) and full-azimuth (FAZ) circular shooting streamer geometries. These datasets were processed and analyzed to test optimal demultiple and imaging approaches and to understand the relative benefits of different acquisition geometries. The results indicated that acquiring azimuthally rich data improves illumination of the pre-carbonate section. The demultiple testing successfully demonstrates the workflow needed to attenuate both short and long period multiples. Reverse Time Migration (RTM) produced an improved image compared to Kirchhoff pre stack depth migration. This case study demonstrates that large scale 3D finite difference seismic modelling, aimed at addressing specific geological challenges, is an effective tool in understanding the optimum acquisition design and processing workflow prior to beginning a new seismic acquisition. Based on the findings of this collaborative study, PETRONAS has now acquired a new seismic data set in the area targeting pre-carbonate objectives using a 3D full-azimuth circular acquisition.