Search

We present a method that skeletonizes the first arriving seismic refractions by machine learning and inverts them for the subsurface velocity model. In this study, first arrivals can be compressed in a low-rank sense with their skeletal features extracted by a well-trained autoencoder. Empirical experiments suggest that the autoencoder’s $1\times 1$ or $2\times 1$ latent vectors vary continuously with respect to the input seismic data. It is, therefore, reasonable to introduce a misfit functional measuring the discrepancies between the predicted and the observed data in a low-dimensional latent space. The benefit of this approach is that an elaborated autoencoding neural network not only refines intrinsic information hidden in the refractions but also improves the quality of inversion for a reliable background velocity model. Numerical tests on both synthetic and field data demonstrate the effectiveness of this method, especially in recovering the low-to-intermediate wavenumber parts of the subsurface velocity distribution. Comparisons are made with the other three relevant methods, the wave-equation travel-time (WT) inversion, the envelope inversion, and the full waveform inversion (FWI). As expected, the cycle skipping problem is alleviated due to the reduction of dimensions of data space. This method outperforms the envelope inversion in resolution, and it is no worse than WT. Moreover, there is no need for careful manual travel-time picking with this methodology. In general, this inversion framework provides an extendable strategy to compress any input data for reconstructing high-dimensional physical parameters.

Information extracted from aerial photographs has been used for many practical applications such as urban planning, forest management, disaster relief, and climate modeling. In many cases, labeling of information in the photo is still performed by human experts, making the process slow, costly, and error-prone. This article shows how a convolutional neural network can be used to determine the location of ground control points (GCPs) in aerial photos, which significantly reduces the amount of human labor in identifying GCP locations. Two convolutional neural network (CNN) methods, sliding-window CNN with superpixel-level majority voting, and superpixel-based CNN are evaluated and analyzed. The results of the classification and segmentation show that both of these methods can quickly extract the locations of objects from aerial photographs, but only superpixel-based CNN can unambiguously locate the GCPs.

We have recast the forward pass of a multilayered convolutional neural network (CNN) as the solution to the problem of sparse least-squares migration (LSM). The CNN filters and feature maps are shown to be analogous, but not equivalent, to the migration Green’s functions and the quasi-reflectivity distribution, respectively. This provides a physical interpretation of the filters and feature maps in deep CNN in terms of the operators for seismic imaging. Motivated by the connection between sparse LSM and CNN, we adopt the neural network version of sparse LSM. Unlike the standard LSM method that finds the optimal reflectivity image, neural network LSM (NNLSM) finds the optimal quasi-reflectivity image and the quasi-migration Green’s functions. These quasi-migration Green’s functions are also denoted as the convolutional filters in a CNN and are similar to migration Green’s functions. The advantage of NNLSM over standard LSM is that its computational cost is significantly less and it can be used for denoising coherent and incoherent noise in migration images. Its disadvantage is that the NNLSM quasi-reflectivity image is only an approximation to the actual reflectivity distribution. However, the quasi-reflectivity image can be used as an attribute image for high-resolution delineation of geologic bodies.

Irregular topography can cause strong scattering and defocusing of propagating surface waves, so it is important to account for such effects when inverting surface waves for shallow S-wave velocity structures. We have developed a 3D surface-wave dispersion inversion method that takes into account the topographic effects modeled by a 3D spectral element solver. The objective function is the frequency summation of the squared wavenumber differences [Formula: see text] along each azimuthal angle of the fundamental mode or higher-order modes of Rayleigh waves in each shot gather. The wavenumbers [Formula: see text] associated with the dispersion curves are calculated using the data recorded along the irregular free surface. Numerical tests on synthetic and field data demonstrate that 3D topographic wave equation dispersion inversion (TWD) can accurately invert for the S-wave velocity model from surface-wave data recorded on irregular topography. Field data tests for data recorded across an Arizona fault demonstrate that, for this example, the 2D TWD model can be as accurate as the 3D tomographic model. This suggests that in some cases, the 2D TWD inversion is preferred over 3D TWD because of its significant reduction in computational costs. Compared to the 3D P-wave velocity tomogram, the 3D S-wave tomogram agrees much more closely with the geologic model taken from the trench log. The agreement with the trench log is even better when the [Formula: see text] tomogram is computed, which reveals a sharp change in velocity across the fault. The localized velocity anomaly in the [Formula: see text] tomogram is in very good agreement with the well log. Our results suggest that integrating the [Formula: see text] and [Formula: see text] tomograms can sometimes give the most accurate estimates of the subsurface geology across normal faults.

To image deeper portions of the earth, geophysicists must record reflection data with much greater source-receiver offsets. The problem with these data is that the signal-to-noise ratio (S/N) significantly diminishes with greater offset. In many cases, the poor S/N makes the far-offset reflections imperceptible on the shot records. To mitigate this problem, we have developed supervirtual reflection interferometry (SVI), which can be applied to far-offset reflections to significantly increase their S/N. The key idea is to select the common pair gathers where the phases of the correlated reflection arrivals differ from one another by no more than a quarter of a period so that the traces can be coherently stacked. The traces are correlated and summed together to create traces with virtual reflections, which in turn are convolved with one another and stacked to give the reflection traces with much stronger S/Ns. This is similar to refraction SVI except far-offset reflections are used instead of refractions. The theory is validated with synthetic tests where SVI is applied to far-offset reflection arrivals to significantly improve their S/N. Reflection SVI is also applied to a field data set where the reflections are too noisy to be clearly visible in the traces. After the implementation of reflection SVI, the normal moveout velocity can be accurately picked from the SVI-improved data, leading to a successful poststack migration for this data set.

We are grateful to the sponsors of the Center for SubsurfaceImaging and Modeling (CSIM) Consortium for their financialsupport. This work was supported by the Natural ScienceFoundation of China (41874134) and Jilin Excellent YouthFund of China (20190103142JH).

Reflection full-waveform inversion (RFWI) can recover the low-wavenumber components of the velocity model along with the reflection wavepaths. However, this requires an expensive least-squares reverse time migration (LSRTM) to construct the perturbation image that can still suffer from cycle-skipping problems. As an inexpensive alternative to LSRTM, we use migration deconvolution (MD) with RFWI. To mitigate cycle-skipping problems, we develop a multiscale reflection phase inversion (MRPI) strategy that boosts the low-frequency data and should only explain the phase information in the recorded data, not its magnitude spectrum. We also use the rolling-offset strategy that gradually extends the offset range of data with an increasing number of iterations. Numerical results indicate that the MRPI + MD method can efficiently recover the low-wavenumber components of the velocity model and is less prone to getting stuck in local minima compared to conventional RFWI.

We present a quasi-elastic wave equation as a function of the pressure variable, which can accurately model PP reflections with elastic amplitude variation with offset effects under the first-order Born approximation. The kinematic part of the quasi-elastic wave equation accurately models the propagation of P waves, whereas the virtual-source part, which models the amplitudes of reflections, is a function of the perturbations of density and Lamé parameters [Formula: see text] and [Formula: see text]. The quasi-elastic wave equation generates a scattering radiation pattern that is exactly the same as that for the elastic wave equation, and only requires the solution of two acoustic wave equations for each shot gather. This means that the quasi-elastic wave equation can be used for true-amplitude linearized waveform inversion (also known as least-squares reverse time migration) of elastic PP reflections, in which the corresponding misfit gradients are with respect to the perturbations of density and the P- and S-wave impedances. The perturbations of elastic parameters are iteratively updated by minimizing the [Formula: see text]-norm of the difference between the recorded PP reflections and the predicted pressure data modeled from the quasi-elastic wave equation. Numerical tests on synthetic and field data indicate that true-amplitude linearized waveform inversion using the quasi-elastic wave equation can account for the elastic PP amplitudes and provide a robust estimate of the perturbations of P- and S-wave impedances and, in some cases, the density. In addition, true-amplitude linearized waveform inversion provides images with a wider bandwidth and fewer artifacts because the PP amplitudes are accurately explained. We also determine the 2D scalar quasi-elastic wave equation for P-SV reflections and the 3D vector equation for PS reflections.

We present a theory for wave-equation inversion of Love-wave dispersion curves, in which the misfit function is the sum of the squared differences between the wavenumbers along the predicted and observed dispersion curves. Similar to inversion of Rayleigh-wave dispersion curves, the complicated Love-wave arrivals in traces are skeletonized as simpler data, namely, the picked dispersion curves in the [Formula: see text] domain. Numerical solutions to the SH-wave equation and an iterative optimization method are then used to invert these dispersion curves for the S-wave velocity model. This procedure, denoted as wave-equation dispersion inversion of Love waves (LWD), does not require the assumption of a layered model or smooth velocity variations, and it is less prone to the cycle-skipping problems of full-waveform inversion. We demonstrate with synthetic and field data examples that LWD can accurately reconstruct the S-wave velocity distribution in a laterally heterogeneous medium. Compared with Rayleigh waves, inversion of the Love-wave dispersion curves empirically exhibits better convergence properties because they are completely insensitive to the P-velocity variations. In addition, Love-wave dispersion curves for our examples are simpler than those for Rayleigh waves, and they are easier to pick in our field data with a low signal-to-noise ratio.

The 2D wave-equation dispersion (WD) inversion method is extended to 3D wave-equation dispersion inversion of surface waves for the shear-velocity distribution. The objective function of 3D WD is the frequency summation of the squared wavenumber [Formula: see text] differences along each azimuth angle of the fundamental or higher modes of Rayleigh waves in each shot gather. The S-wave velocity model is updated by the weighted zero-lag crosscorrelation between the weighted source-side wavefield and the back-projected receiver-side wavefield for each azimuth angle. A multiscale 3D WD strategy is provided, which starts from the pseudo-1D S-velocity model, which is then used to get the 2D WD tomogram, which in turn is used as the starting model for 3D WD. The synthetic and field data examples demonstrate that 3D WD can accurately reconstruct the 3D S-wave velocity model of a laterally heterogeneous medium and has much less of a tendency to getting stuck in a local minimum compared with full-waveform inversion.

Viscoacoustic migration can significantly compensate for the amplitude loss and phase distortion in migration images computed from highly attenuated data. However, solving the viscoacoustic wave equation requires a significant amount of storage space and computation time, especially for least-squares migration methods. To mitigate this problem, we used acoustic reverse time migration (RTM) instead of viscoacoustic migration to migrate the viscoacoustic data and then we correct the amplitude and phase distortion by hybrid deblurring filters in the image domain. Numerical tests on synthetic and field data demonstrate that acoustic RTM combined with hybrid deblurring filters can compensate for the attenuation effects and produce images with high resolution and balanced amplitudes. This procedure requires less than one-third of the storage space and is [Formula: see text] times faster compared with the viscoacoustic migration, but at the cost of mildly reduced accuracy. Here, [Formula: see text] represents the number of iterations used for least-squares migration method. This method can be extended to 3D migration at even a greater cost saving.

Velocity errors in the shallow part of the velocity model can lead to erroneous estimates of the full-waveform inversion (FWI) tomogram. If the location and topography of a reflector are known, then such a reflector can be used as a reference reflector to update the underlying velocity model. Reflections corresponding to this reference reflector are windowed in the data space. Windowed reference reflections are then crosscorrelated with reflections from deeper interfaces, which leads to partial cancellation of static errors caused by the overburden above the reference interface. Interferometric FWI (IFWI) is then used to invert the tomogram in the target region, by minimizing the normalized waveform misfit between the observed and predicted crosscorrelograms. Results with synthetic and field data with static errors above the reference interface indicate that an accurate tomogram can be inverted in areas lying within several wavelengths of the reference interface. IFWI can also be applied to synthetic time-lapse data to mitigate the nonrepeatability errors caused by time-varying overburden variations. The synthetic- and field-data examples demonstrate that IFWI can provide accurate tomograms when the near surface is ridden with velocity errors.

We thanks to the Deep Earth Imaging Future Science Platformof CSIRO for funding and computing resources of CSIRO.

We thank the financial support from the sponsors of the Consortium of Subsurface Imaging and Fluid Modeling (CSIM). We thank KAUST for funding this research. This work is supported by the Natural Science Foundation of China (41874134) and the China Postdoctoral Science Foundation 2106T902503 and 2015M571366.

We have developed a wave-equation traveltime inversion method with multifrequency bands to invert for the shallow or intermediate subsurface velocity distribution. Similar to the classical wave-equation traveltime inversion, this method searches for the velocity model that minimizes the squared sum of the traveltime residuals using source wavelets with progressively higher peak frequencies. Wave-equation traveltime inversion can partially avoid the cycle-skipping problem by recovering the low-wavenumber parts of the velocity model. However, we also use the frequency information hidden in the traveltimes to obtain a more highly resolved tomogram. Therefore, we use different frequency bands when calculating the Fréchet derivatives so that tomograms with better resolution can be reconstructed. Results are validated by the zero-offset gathers from the raw data associated with moderate geometric irregularities. The improved wave-equation traveltime method is robust and merely needs a rough estimate of the starting model. Numerical tests on the synthetic and field data sets validate the above claims.

Two robust imaging technologies are reviewed that provide subsurface geologic information in challenging environments. The first one is wave-equation dispersion (WD) inversion of surface waves and guided waves (GW) for the shear-velocity (S-wave) and compressional-velocity (P-wave) models, respectively. The other method is traveltime inversion for the velocity model, in which supervirtual refraction interferometry (SVI) is used to enhance the signal-to-noise ratio of far-offset refractions. We have determined the benefits and liabilities of both methods with synthetic seismograms and field data. The benefits of WD are that (1) there is no layered-medium assumption, as there is in conventional inversion of dispersion curves. This means that 2D or 3D velocity models can be accurately estimated from data recorded by seismic surveys over rugged topography, and (2) WD mostly avoids getting stuck in local minima. The liability is that WD for surface waves is almost as expensive as full-waveform inversion (FWI) and, for Rayleigh waves, only recovers the S-velocity distribution to a depth no deeper than approximately 1/2 to 1/3 wavelength of the lowest-frequency surface wave. The limitation for GW is that, for now, it can estimate the P-velocity model by inverting the dispersion curves from GW propagating in near-surface low-velocity zones. Also, WD often requires user intervention to pick reliable dispersion curves. For SVI, the offset of usable refractions can be more than doubled, so that traveltime tomography can be used to estimate a much deeper model of the P-velocity distribution. This can provide a more effective starting velocity model for FWI. The liability is that SVI assumes head-wave first arrivals, not those from strong diving waves.

Seismic time-lapse surveys are susceptible to repeatability errors due to varying environmental conditions. To mitigate this problem, we propose the use of interferometric least-squares migration to estimate the migration images for the baseline and monitor surveys. Here, a known reflector is used as the reference reflector for interferometric least-squares migration, and the data are approximately redatumed to this reference reflector before imaging. This virtual redatuming mitigates the repeatability errors in the time-lapse migration image. Results with synthetic and field data show that interferometric least-squares migration can sometimes reduce or eliminate artifacts caused by non-repeatability in time-lapse surveys and provide a high-resolution estimate of the time-lapse change in the reservoir.

Different types of model parameters, such as the P- and S-wave velocities, can be coupled to one another in multiparameter seismic inversion. This coupling effect is not taken into account by the conventional approximation to the Hessian inverse for a single type of parameter, and so it can lead to a slow rate of convergence by an iterative inversion method. To solve this problem, we have developed a multiparameter deblurring filter that approximates the Hessian inverse. This filter takes into account the coupling between different parameters by using stationary local filters to approximate the submatrices of the Hessian inverse for the different types of parameters. The filters are calculated by matching the reference multiparameter migration images to their reference reflectivity models. Numerical tests with elastic migration and inversion indicate that the multiparameter deblurring filter not only reduces the footprint noise, balances the amplitude, and increases the resolution of the elastic migration images, but it also mitigates the crosstalk artifacts caused by the coupling effect. When used as a preconditioner, it also accelerates the convergence rate for elastic inversion.

We present the theory for wave-equation inversion of dispersion curves obtained from guided P waves. The misfit function is the sum of the squared differences between the wavenumbers along the predicted and observed dispersion curves, and the inverted result is a high-resolution estimate of the near-surface P-velocity model. This procedure, denoted as the wave equation dispersion inversion of guided Pwaves (WDG), is valid for near-surface waveguides with irregular layers. It is less prone to the cycle skipping problems of full waveform inversion (FWI) and can sometimes provide velocity models with higher resolution than wave-equation traveltime tomography (WT). The synthetic and field data examples demonstrate that WDG for guided P-waves can accurately reconstruct the P-wave velocity distribution in laterally heterogeneous media.

The theory and practice of multisource full waveform inversion of marine supergathers are described with a frequency-selection strategy. The key enabling property of frequency selection is that it eliminates the crosstalk among sources, thus overcoming the aperture mismatch of marine multisource inversion. Tests on multisource full waveform inversion of synthetic marine data and Gulf of Mexico data show speedups of 4× and 8×, respectively, compared to conventional full waveform inversion. This article is protected by copyright. All rights reserved

A resonant multiple is defined as a multiple reflection that revisits the same subsurface location along coincident reflection raypaths. We have found that resonant first-order multiples can be migrated with either Kirchhoff or wave-equation migration methods to give images with approximately twice the spatial resolution compared with poststack primary-reflection images. We have developed a moveout-correction stacking method to enhance the signal-to-noise ratios of the resonant multiples for superresolution migration. The effectiveness of this procedure is validated by synthetic and field data tests.

We present a scheme for multiscale phase inversion (MPI) of seismic data that is less sensitive than full-waveform inversion (FWI) to the unmodeled physics of wave propagation and to a poor starting model. To avoid cycle skipping, the multiscale strategy temporally integrates the traces several times, i.e., high-order integration, to produce low-boost seismograms that are used as input data for the initial iterations of MPI. As the iterations proceed, lower frequencies in the data are boosted by using integrated traces of lower order as the input data. The input data are also filtered into different narrow frequency bands for the MPI implementation. Numerical results with synthetic acoustic data indicate that, for the Marmousi model, MPI is more robust than conventional multiscale FWI when the initial model is moderately far from the true model. Results from synthetic viscoacoustic and elastic data indicate that MPI is less sensitive than FWI to some of the unmodeled physics. Inversion of marine data indicates that MPI is more robust and produces modestly more accurate results than FWI for this data set.

The theory for the seismic scanning tunneling macroscope is extended from acoustic body waves to elastic body-wave propagation. We show that, similar to the acoustic case, near-field superresolution imaging from elastic body waves results from the O(1/R) term, where R is the distance between the source and near-field scatterer. The higher-order contributions $$R^{-n}$$ for $$n>1$$ are cancelled in the near-field region for a point source with normal stress.

We thank the sponsors for supporting the Consortium of Subsurface Imaging and Fluid Modeling (CSIM). We also thank KAUST for providing funding by the CRG grant OCRF-2014-CRG3-2300. For computer time, this research used the resources of the IT Research Computing Group and the Supercomputing Laboratory at KAUST.

Wave-equation migration velocity analysis (WEMVA) based on subsurface-offset, angle domain, or time-lag common-image gathers (CIGs) requires significant computational and memory resources because it computes higher dimensional migration images in the extended image domain. To mitigate this problem, we have developed a WEMVA method using plane-wave CIGs. Plane-wave CIGs reduce computational cost and memory storage because they are directly calculated from prestack plane-wave migration and the number of plane waves is often much smaller than the number of shots. In the case of an inaccurate migration velocity, the moveout of plane-wave CIGs is automatically picked by a semblance analysis method, which is then linked to the migration velocity update by a connective function. Numerical tests on two synthetic data sets and a field data set validate the efficiency and effectiveness of this method.

We have developed a tutorial for skeletonized inversion of pseudo-acoustic anisotropic vertical symmetry axis (VTI) data. We first invert for the anisotropic models using wave-equation traveltime inversion. Here, the skeletonized data are the traveltimes of transmitted and/or reflected arrivals that lead to simpler misfit functions and more robust convergence compared with full-waveform inversion. This provides a good starting model for waveform inversion. The effectiveness of this procedure is illustrated with synthetic data examples and a marine data set recorded in the Gulf of Mexico.

Full-waveform inversion of seismic data is often plagued by cycle-skipping problems such that an iterative optimization method often gets stuck in a local minimum. To avoid this problem, we simplify the objective function so that the iterative solution can quickly converge to a solution in the vicinity of the global minimum. The objective function is simplified by only using parsimonious and important portions of the data, which are defined as skeletonized data. We have developed a mostly nonmathematical tutorial that explains the theory of wave-equation inversion of skeletonized data. We also demonstrate its effectiveness with examples.

We have developed a methodology for detecting the presence of near-surface heterogeneities by naturally migrating backscattered surface waves in controlled-source data. The near-surface heterogeneities must be located within a depth of approximately one-third the dominant wavelength [Formula: see text] of the strong surface-wave arrivals. This natural migration method does not require knowledge of the near-surface phase-velocity distribution because it uses the recorded data to approximate the Green’s functions for migration. Prior to migration, the backscattered data are separated from the original records, and the band-passed filtered data are migrated to give an estimate of the migration image at a depth of approximately one-third [Formula: see text]. Each band-passed data set gives a migration image at a different depth. Results with synthetic data and field data recorded over known faults validate the effectiveness of this method. Migrating the surface waves in recorded 2D and 3D data sets accurately reveals the locations of known faults. The limitation of this method is that it requires a dense array of receivers with a geophone interval less than approximately one-half [Formula: see text].

Least-squares migration (LSM) can produce images with better balanced amplitudes and fewer artifacts than standard migration. The conventional objective function used for LSM minimizes the L2-norm of the data residual between the predicted and the observed data. However, for field-data applications in which the recorded data are noisy and undersampled, the conventional formulation of LSM fails to provide the desired uplift in the quality of the inverted image. We have developed a least-squares reverse time migration (LSRTM) method using local Radon-based preconditioning to overcome the low signal-to-noise ratio (S/N) problem of noisy or severely undersampled data. A high-resolution local Radon transform of the reflectivity is used, and sparseness constraints are imposed on the inverted reflectivity in the local Radon domain. The sparseness constraint is that the inverted reflectivity is sparse in the Radon domain and each location of the subsurface is represented by a limited number of geologic dips. The forward and the inverse mapping of the reflectivity to the local Radon domain and vice versa is done through 3D Fourier-based discrete Radon transform operators. The weights for the preconditioning are chosen to be varying locally based on the relative amplitudes of the local dips or assigned using quantile measures. Numerical tests on synthetic and field data validate the effectiveness of our approach in producing images with good S/N and fewer aliasing artifacts when compared with standard RTM or standard LSRTM.

We use elastic least-squares reverse time migration (LSRTM) to invert for the reflectivity images of P- and S-wave impedances. Elastic LSRTM solves the linearized elastic-wave equations for forward modeling and the adjoint equations for backpropagating the residual wavefield at each iteration. Numerical tests on synthetic data and field data reveal the advantages of elastic LSRTM over elastic reverse time migration (RTM) and acoustic LSRTM. For our examples, the elastic LSRTM images have better resolution and amplitude balancing, fewer artifacts, and less crosstalk compared with the elastic RTM images. The images are also better focused and have better reflector continuity for steeply dipping events compared to the acoustic LSRTM images. Similar to conventional least-squares migration, elastic LSRTM also requires an accurate estimation of the P- and S-wave migration velocity models. However, the problem remains that, when there are moderate errors in the velocity model and strong multiples, LSRTM will produce migration noise stronger than that seen in the RTM images.

We present a parsimonious wave-equation travel-time inversion technique for refraction waves. A dense virtual refraction dataset can be generated from just two reciprocal shot gathers for the sources at the endpoints of the survey line, with N geophones evenly deployed along the line. These two reciprocal shots contain approximately 2N refraction travel times, which can be spawned into O(N2) refraction travel times by an interferometric transformation. Then, these virtual refraction travel times are used with a source wavelet to create N virtual refraction shot gathers, which are the input data for wave-equation travel-time inversion. Numerical results show that the parsimonious wave-equation travel-time tomogram has about the same accuracy as the tomogram computed by standard wave-equation travel-time inversion. The most significant benefit is that a reciprocal survey is far less time consuming than the standard refraction survey where a source is excited at each geophone location.

Many explorationists think of surface waves as the most damaging noise in land seismic data. Thus, much effort is spent in designing geophone arrays and filtering methods that attenuate these noisy events. It is now becoming apparent that surface waves can be a valuable ally in characterizing the near-surface geology. This review aims to find out how the interpreter can exploit some of the many opportunities available in surface waves recorded in land seismic data. For example, the dispersion curves associated with surface waves can be inverted to give the S-wave velocity tomogram, the common-offset gathers can reveal the presence of near-surface faults or velocity anomalies, and back-scattered surface waves can be migrated to detect the location of near-surface faults. However, the main limitation of surface waves is that they are typically sensitive to S-wave velocity variations no deeper than approximately half to one-third the dominant wavelength. For many exploration surveys, this limits the depth of investigation to be no deeper than approximately 0.5–1.0 km.

We present a wave-equation inversion method that inverts skeletonized seismic data for the subsurface velocity model. The skeletonized representation of the seismic traces consists of the low-rank latent-space variables predicted by a well-trained autoencoder neural network. The input to the autoencoder consists of seismic traces, and the implicit function theorem is used to determine the Fréchet derivative, i.e., the perturbation of the skeletonized data with respect to the velocity perturbation. The gradient is computed by migrating the shifted observed traces weighted by the skeletonized data residual, and the final velocity model is the one that best predicts the observed latent-space parameters. We denote this as inversion by Newtonian machine learning (NML) because it inverts for the model parameters by combining the forward and backward modeling of Newtonian wave propagation with the dimensional reduction capability of machine learning. Empirical results suggest that inversion by NML can sometimes mitigate the cycle-skipping problem of conventional full-waveform inversion (FWI). Numerical tests with synthetic and field data demonstrate the success of NML inversion in recovering a low-wavenumber approximation to the subsurface velocity model. The advantage of this method over other skeletonized data methods is that no manual picking of important features is required because the skeletal data are automatically selected by the autoencoder. The disadvantage is that the inverted velocity model has less resolution compared with the FWI result, but it can serve as a good initial model for FWI. Our most significant contribution is that we provide a general framework for using wave-equation inversion to invert skeletal data generated by any type of neural network. In other words, we have combined the deterministic modeling of Newtonian physics and the pattern matching capabilities of machine learning to invert seismic data by NML.