Wavefield Decomposition Method Research Articles

The qP- and qSV-wave separation in 2D vertically transversely isotropic media and their applications in elastic reverse time migration

Compressional and shear wave separation is an important step in anisotropic elastic reverse time migration (ERTM). With separated wavefields, we can generate images between different wave types and reveal more subsurface physical properties, as well as remove unwanted crosstalk and improve image quality. Traditional Helmholtz decomposition based on divergence and curl operations for isotropic media cannot be directly extended to vertically transversely isotropic (VTI) media. Wave-mode separation methods, similar to the nonstationary spatial filter and Poynting vector, cannot preserve the phase and amplitude information of the original coupled wavefield, thus limiting their application in ERTM. Currently, the anisotropic wavefield decomposition methods include the wavenumber-domain approach, the low-rank approximation, and other separation approaches based on different approximations, e.g., weak or elliptical anisotropy. These methods are either computationally intensive or involve large separation errors, especially in models with strong heterogeneities and anisotropy. The VTI wavefield separation operator is essentially defined in a mixed space-wavenumber domain. We introduce scalar operators to transfer operations from this mixed space-wavenumber domain to the space domain. The local wave propagation direction in the scalar operators is estimated using the Poynting vector, which is highly computationally efficient in the space domain. When we apply the wave separation method to ERTM, we not only obtain the vectorized qP and qSV waves but also retrieve the scalar qP wavefield. The scalar and vector wavefields preserve the amplitude and phase information in the original coupled wavefield. For anisotropic ERTM, we suggest using the scalar imaging condition to generate the PP image and the magnitude- and sign-based vector imaging condition to produce the PS image, both having higher image accuracy than the dot-product imaging condition. Numerical examples are used to validate our anisotropic wavefield separation method and the related ERTM workflow.

Reverse time migration angle gathers in acoustic anisotropic media using direction vectors

The application of direction vectors in the generation of reverse time migration (RTM) angle gathers in complex acoustic anisotropic media often encounters three main challenges: not pointing to the phase-velocity direction (PVD) of the Poynting vector, inaccuracy due to overlapping wavefields, and instability due to zero points of the direction vector. In general anisotropic media, the standard Poynting vector indicates the group-velocity direction (GVD), whereas reflection and transmission phenomena rely on the PVD. Anisotropy introduces discrepancies between the GVD and the PVD. To overcome this issue, we use the so-called PVD vector to directly calculate the PVD from anisotropic wavefields, eliminating the need for an approximated conversion from the GVD to the PVD. To mitigate the inaccuracy problem, we apply the Hilbert transform-based wavefield decomposition method to separate overlapping wavefields into their up/down components, and then we calculate the PVDs using the separated wavefields. To tackle the instability problem, we incorporate the additionally simulated quadrature wavefield during the wavefield decomposition procedure. By combining the direction vector of the quadrature wavefield with that of the original wavefield, we can eliminate the zero points and thus obtain a stabilized PVD vector. With those problems solved or alleviated, we present a scheme for the generation of anisotropic RTM angle gathers in complex areas. Two numerical examples using synthetic data sets demonstrate our method’s effectiveness.

Seismic reverse time migration based on biaxial wavefield decomposition

Seismic images generated by reverse time migration (RTM) are often contaminated by low-wavenumber artifacts. Decomposition of up- and downgoing waves during propagation can effectively reduce these artifacts. However, decomposition along the vertical axis is not sufficient to eliminate the migration artifacts. We introduce biaxial wavefield decomposition to decompose seismic wavefields in the lateral and vertical directions. After decomposition, there are 16 crosscorrelation terms between different components, and there are only four crosscorrelation terms that can be used separately to image flat layers and tilted interfaces. Within these four terms, the remaining artifacts can be easily identified and reduced by other methods, and the resolution can be further improved by adjusting the weights of the four terms. The decomposition method based on the analytic wavefield via the Hilbert transform has been shown to agree with the viscoacoustic wave equation with fractional spatial derivatives. Therefore, this biaxial wavefield decomposition method is suitable for decomposing the wavefields for the attenuation-compensated RTM. The decomposed and recombined migration results can better image a complex subsurface structure with high resolution.

Converted-wave reflection traveltime inversion with free-surface multiples for ocean-bottom-node data

Advancements in ocean-bottom-node (OBN) technologies have enabled the recording of high-quality multicomponent seismic data, which can provide crucial information about subsurface properties through elastic seismic imaging. However, imaging multicomponent data is challenging due to the requirement for P- and S-wave velocity models. Conventional processing relies on estimating the S-wave velocity model using a primary converted wave (C wave) because the pure S-wave primaries are not available if using standard air-gun sources. When a free surface is present, the C waves recorded with geophones can be contaminated by source-side P-wave multiples, especially water-layer reverberations, which are difficult to remove even using state-of-the-art wavefield decomposition methods. These C-wave multiples can cause significant mismatches between synthetic and observed data during reflection traveltime or waveform inversion if only primary P-to-S (PS) conversions are simulated. To address these issues and reliably build an S-wave velocity model for PS imaging, we develop a C-wave reflection traveltime inversion method that takes source-side free-surface multiples into account. The traveltime misfit of C-wave data is extracted using dynamic image warping to ensure a better match between the synthetic and observed data. The functional gradient of the objective function is derived using the adjoint state method. Based on sensitivity kernel analyses, a convenient gradient preconditioning strategy is developed to robustly separate the contribution of primary and multiple C waves during backpropagation. This strategy can effectively mitigate the high-wavenumber crosstalk associated with nonphysical wavepaths in the gradient and thus allow more accurate updating for the S-wave velocity model. The synthetic data and shallow-water OBN data from the East China Sea indicate that this approach can reliably recover the S-wave macrovelocity structures for PS imaging.

Crosswell frequency-domain reverse time migration imaging with wavefield decomposition

Abstract Crosswell seismic technique can provide high-resolution imaging and monitoring of the subsurface. However, compared to the surface seismic, crosswell data contain more complex wave components, which increases the difficulty of seismic processing and migrations. Conventional acoustic reverse time migration (RTM) mainly uses the cross-correlation of the source forward and the receiver backward wavefields. The redundant information generated by cross-correlation may undermine the imaging reliability in the crosswell models. Thus, we develop a novel wavefield decomposition imaging condition and only used cross-correlation information of incident and reflection waves in the same propagation directions, which eliminates the artifacts generated from the cross-correlation of wave information unrelated to reflections in the crosswell image. We perform RTM in the frequency domain to maintain efficiency in multi-shot problems. A mono-frequency wavefield decomposition method is applied to separate and process the seismic data. The forward and backward wavefields are reclassified into the up- and downpropagating components. The L1 norm is introduced to enhance the robustness of the proposed imaging method. We then use this method to synthesize data from layering models and analyze the imaging results generated from each pair of cross-correlations using source and receiver wavefields. Results show that the cross-correlation information belonging to the same propagation contributes most to the crosswell image, and the other information always generates migration noises. Moreover, we apply the proposed method to a real-field dataset. Processing results validate the effectiveness of the proposed means for eliminating false events in the crosswell models and improving image quality.

Hierarchical wave-mode separation in the poroelastic medium using eigenform analysis

In the elastic medium, the scalar and vector P- and S-waves decomposition has been extensively studied and some strategies can be extended to the poroelastic medium to extract P- and S-wavefields. However, there are three propagation modes in the poroelastic medium in Biot’s theory, namely, a fast P wave, a slow P wave, and an S wave. Because the propagation characteristic of a slow P wave is different from that of a fast P wave and S wave, the wavefield separation methods in the elastic medium cannot be directly applied to the poroelastic medium to produce a complete wave-mode separation. Based on the eigenform analysis, we have developed a hierarchical wavefield decomposition method to completely separate S waves and fast and slow P waves in the poroelastic medium. Using the Helmholtz decomposition, we first compute scalar and vector potential wavefields to separate P and S waves. Then, a cross-product operator is proposed to decompose fast and slow P waves based on their different polarization directions. To produce correct amplitudes and phases, we apply another cross-product operator and an amplitude correction term to the separated wavefields. Three numerical examples demonstrate that our method can produce accurate fast P-wave, slow P-wave, and S-wave separation results, and the decomposed fast and slow P waves have the same phases and amplitudes as the P-wave potential wavefields.

Suppressing migration noise in reverse time migration of vertical seismic profiles by multiple stacking estimation

Reverse time migration (RTM) has advantages of imaging steep dip structures but may inevitably produce migration swings when applied to vertical seismic profiles (VSPs). These artifacts often are not suppressed by stacking over shots. The noise is mainly generated from the limited coverage of VSP acquisition geometry, which leads to uneven superposition. The noise is processed by the imaging condition and then overlaps with the real layers, generating artificial images with different dip angles that are difficult to distinguish. The difference is that the real layers can be enhanced by the redundancy of multifold data, whereas the noise is usually generated by a smaller subset of the data. A VSP RTM processing flow is developed to suppress this type of noise. This method constructs the dip angle gather and estimates the redundancy of the coverage in the dip domain. In the stratigraphic dip domain, the false structure and the real structure are located at different dip angles. The real dip is determined by identifying the center of the stacking number. Then, the noise is suppressed by attenuating the imaging information of the other dip angles. The efficient wavefield decomposition method based on Hilbert transform-based instantaneous wavenumber method is adopted to determine the wave propagation direction. Compared with the Poynting vector method and the instantaneous wavenumber method, this approach alleviates the influence of multiwave arrivals and is suitable for complex media, ensuring the accuracy of image decomposition. The whole process is data-driven and does not require prior information. Theoretical and field examples indicate that this method can significantly reduce migration noise and improve the image quality of VSP RTM.

Angle-Weighted Reverse Time Migration With Wavefield Decomposition Based on the Optical Flow Vector

Reverse time migration (RTM) is based on the two-way wave equation, so its imaging results obtained by conventional zero-lag cross-correlation imaging conditions contain a lot of low-wavenumber noises. So far, the wavefield decomposition method based on the Poynting vector has been developed to suppress these noises; however, this method also has some problems, such as unstable calculation of the Poynting vector, low accuracy of wavefield decomposition, and poor effect of large-angle migration artifacts suppression. This article introduces the optical flow vector method to RTM to realize high-precision wavefield decomposition for both the source and receiver wavefields and obtains four directions of wavefields: up-, down-, left-, and right-going. Then, the cross-correlation imaging sections of one-way propagation components of forward- and back-propagated wavefields are optimized and stacked. On this basis, the reflection angle of each imaging point is calculated based on the optical flow vector, and an attenuation factor related to the reflection angle is introduced as the weight to generate the optimal stack images. The tests of theoretical model and field marine seismic data illustrate that compared with the conventional RTM with wavefield decomposition based on the Poynting vector, the angle-weighted RTM with wavefield decomposition based on the optical flow vector proposed in this article can achieve wavefield decomposition for both the source and receiver wavefields and calculate the reflection angle of each imaging point more accurately and stably. Moreover, the proposed method adopts angle weighting processing, which can further eliminate large-angle migration artifacts and effectively improve the imaging accuracy of RTM.

A high-efficiency wavefield decomposition method based on the Hilbert transform

Wavefield decomposition can be used to extract effective information in reverse time migration and full-waveform inversion. The wavefield decomposition methods based on the Hilbert transform (HTWD) and the Poynting vector (PVWD) are the most commonly used. The HTWD needs to save the wavefields at all time steps or introduce additional numerical simulation, which increases the computational cost. PVWD cannot handle multiwave arrivals, and its performance is poor in complex situations. We have developed an efficient wavefield decomposition method based on the Hilbert transform (EHTWD). EHTWD constructs two wavefields to replace the original wavefield and the wavefield after the Hilbert transform. The first wavefield is obtained by using the dispersion relation to modify the frequency components. The other wavefield is obtained by time difference approximation. Therefore, there is a 90° phase change between the two wavefields. In EHTWD, we only need two wavefields at different moments, which avoids the need for additional numerical simulation. EHTWD is also suitable for wavefield decomposition in arbitrary directions. Compared to HTWD, the computational complexity can be greatly reduced with the decrease of the number of imaging time slices. The numerical examples of wavefield decomposition demonstrate that our method can realize wavefield decomposition in any direction. The examples of imaging decomposition and real data also indicate that EHTWD suppresses imaging noise effectively.

An improved approach for normal-incidence sound transmission loss measurement using a four-microphone impedance tube

The main aim of this study is to introduce an improved method for determining the sound properties of acoustic materials which is more precise than the common wavefield decomposition method and simpler than the common transfer matrix method. In the first part of the article, a group of formulae for calculating sound transmission loss is represented by combining the wavefield decomposition and transfer matrix methods. Subsequently, a formula for calculating sound absorption coefficients is derived from these formulae by definition. Furthermore, the present formulae are validated by comparing the experimental results achieved with the present formulae and those results obtained by other methods recorded in published articles. Eventually, it is demonstrated that the method can accurately measure the sound insulation performance of materials and the sound absorption properties of limp and lightweight materials.

Stable and High-Efficiency Attenuation Compensation in Reverse-Time Migration Using Wavefield Decomposition Algorithm

Q -compensated reverse time migration ( Q -RTM) can compensate seismic attenuation caused by the anelastic behavior of subsurface media. Although, the traditional Q -RTM has high computational efficiency, it is instable because the high frequency or wavenumber ambient noise is exponentially boosted during forward and backward seismic wavefield propagation. The existing stable Q -RTM method costs twice as much computing time and memory compared to the traditional Q -RTM. In this letter, we propose a new Q -RTM method to address the above issues simultaneously. First, a new viscoacoustic wave equation is derived based on a wavefield decomposition method to obtain the velocity-dispersion-only and viscoacoustic wavefields efficiently. Then, a theoretical framework of stable and high-efficiency Q -RTM method is proposed based on the velocity-dispersion-only and viscoacoustic wavefields. The synthetic example shows that the new stable Q -RTM results match well with the reference images (without attenuation images). Moreover, the field data images also demonstrate the stability and high-efficiency of our proposed Q -RTM method.

Passive multiple reverse time migration imaging based on wave decomposition and normalized imaging conditions

With the development of seismic exploration, passive-source seismic data has attracted increasing attention. Ambient noise passive seismic sources exists widely in nature and industrial production. Passive seismic data is important in logging while drilling (LWD), large-scale structural exploration, etc. In this paper, we proposed a passive multiple reverse time migration imaging (PMRTMI) method based on wavefield decomposition and normalized imaging conditions method. This method differs from seismic interferometry in that it can use raw passive seismic data directly in RTM imaging without reconstruction of virtual active gather, and we use the wavefield decomposition method to eliminate the low frequency noise in RTM. Further, the energy normalized imaging condition is used in full wavefield decomposition, which can not only enhance the image quality of both edge and deep information but also overcome the wrong energy problem caused by uneven distribution of passive sources; furthermore, this method exhibits high efficiency. Finally, numerical examples with the Marmousi model show the effectiveness of the method.

Prestack elastic RTM for VTI media using vector wavefield decomposition and vector imaging conditions

ABSTRACTTo date, studies of elastic reverse time migration (RTM) have undergone much improvement, but mainly focus on isotropic media. Because anisotropic media is widespread, it is necessary to explore the application of elastic RTM in anisotropic media. We extend the wavefield decomposition method, which is based on decoupled propagation, and three vector imaging conditions to transversely isotropic media with vertical symmetry isotropy (VTI). Two of the imaging conditions are based on the excitation amplitude (EA) and the third is based on source-normalized cross-correlation. First, the wavefield decomposition method is extended to VTI media. This is then tested in a two-layer model. The results show that this extension cannot decompose P- and S-waves perfectly in VTI media; some weak residuals remain. However, the results of a simple model test show that no obvious crosstalk is generated by these weak residuals. Finally, a Hess VTI model is adopted to test the adaptive use of this method in complex media. Many subsurface structures can be clearly recognized in the migrated result, for example, the high-velocity rock body, a fault and two low-velocity interlayers. Compared with PP images, converted PS images have many merits, such as clearer imaging of the anisotropic body, higher resolution and a wider migration aperture. We conclude that the vector decomposition method and three vector imaging conditions can be applied to prestack elastic RTM for VTI media and satisfactory results obtained.

Enhancing low-wavenumber components of full-waveform inversion using an improved wavefield decomposition method in the time-space domain

Full-waveform inversion (FWI) is a successful technique that attempts to build high-resolution velocity models by minimizing the residuals between observed and calculated data after a series of iterations. However, its successful application still relies on the quality of the estimated initial velocity model, as well as on the design of a suitable hierarchical workflow. If the initial model is far away from the true model, the widely used gradient-based FWI will inevitably get trapped into meaningless local minima. Given that the properties of the low-wavenumber components can help construct the initial model and overcome cycle skipping, we propose a low-wavenumber update FWI method based on the Hilbert transform called the hybrid method. This method provides a new form of FWI in the time-space domain that appears to be capable of enhancing the low-wavenumber updates. In this method, the updates of low wavenumbers and high wavenumbers are divided into two parts. In the first step, a gradient formula based on the Hilbert transform plays a key role in updating the low-wavenumber components. The gradient is formed by cross-correlation of upgoing and downgoing waves of source and residual wavefields. Instead of using the wavefield decomposition method in the frequency-wavenumber domain with high computational complexity, we apply a novel method to extract the upgoing and downgoing wavefields in the time-space domain. In addition to the acceptable computational burden, the memory is reduced greatly. The second step is the standard FWI, and its task is to update the high-wavenumber components. Numerical tests on an anomaly model and the Marmousi model demonstrate the effectiveness of the hybrid method for extracting low-wavenumber components even for noisy data or when low-frequency data are missing.

Reverse-time migration using multidirectional wavefield decomposition method

Reverse-time migration has attracted more and more attention owing to the advantages of high imaging accuracy, no dip restriction, and adaptation to complex velocity models. Cross-correlation imaging method is typically used in conventional reverse-time migration that produces images with strong low-frequency noise. Wavefield decomposition imaging can suppress such noise; however, some residual noise persists in the imaging results. We propose a 2D multidirectional wavefield decomposition method based on the traditional wavefield decomposition method. First, source wavefields and receiver wavefields are separated into eight subwavefields, respectively. Second, cross-correlation imaging is applied to selected subwavefields to produce subimages. Finally, the subimages are stacked to generate the final image. Numerical examples suggest that the proposed method can eliminate the low-frequency noise effectively and produce high-quality imaging profiles.

Array analysis of the seismic wavefield of long‐period events and volcanic tremor at Arenal volcano, Costa Rica

AbstractWe use wavefield decomposition methods (time domain cross correlation and frequency domain multiple‐signal classification) to analyze seismic data recorded by a dense, small‐aperture array located 2 km West of Arenal volcano, Costa Rica, and operated during 2.5 days. The recorded wavefield is dominated by harmonic tremor and includes also spasmodic tremor and long‐period (LP) events. We find that the initial stages of LP events are characterized by three different wave arrivals. These arrivals propagate with similar back azimuths pointing to the volcano summit (∼80°N) and increasing apparent slowness of 0.4, 1.1, and 1.7 s/km. Spasmodic tremors cannot be regarded as coherent signals. On the contrary, harmonic tremors are highly coherent, characterized by the stability of the apparent slowness vector estimates. Apparent slowness lays in the range 1–2 s/km. Back azimuths point in the general direction of the volcano but with a large variability (40–120°N). Nevertheless, there are long‐term variations and evidences of multiple simultaneous components in the harmonic tremor wavefield. These observations suggest that LP events and tremor are generated in a shallow source area near the volcano summit, although they do not share exactly the same source region or source processes. The tremor source is located in the shallowest part of the plumbing system, beneath the lava crust. This dynamic region is subject to complex fluctuations of the physical conditions. Degassing events at different locations of this region might generate variable seismic radiation patterns. The effects of topography and heterogeneous shallow structure of the volcano may amplify these variations and produce the wide directional span observed for volcanic tremor. On the other hand, the LP source seems to be more repeatable. LP events are likely triggered by fragmentation of the fluid flow in a slightly deeper portion of the volcanic conduits.

Wavefield decomposition based on acoustic reciprocity: Theory and applications to marine acquisition

In marine seismic acquisition, the free surface generates seismic events in our recorded data that are often categorized as noise because these events do not contain independent information about the subsurface geology. Ghost events are considered as such noise because these events are generated when the energy generated by the seismic source, as well as any upgoing wavefield propagating upward from the subsurface, is reflected downward by the free surface. As a result, complex interference patterns between up- and downgoing wavefields are present in the recorded data, affecting the spectral bandwidth of the recorded data negatively. The interpretability of the data is then compromised, and hence it is desirable to remove the ghost events from the data. Rayleigh’s reciprocity theorem is used to derive the relevant equations for wavefield decomposition for multisensor and single-sensor data, for depth-varying and depth-independent recordings from marine seismic experiments using a single-source or dual-source configuration. A comparison is made between the results obtained for a 2D synthetic example designed to highlight the strengths and weaknesses of the various acquisition configurations. It is demonstrated that, using the proposed wavefield decomposition method, multisensor data (measurements of pressure and particle velocity components, or multidepth pressure measurements) allow for optimal wavefield decomposition as independent measurements are used to eliminate the interference patterns caused by the free surface. Single-sensor data using constant-depth recordings are found to be incapable of producing satisfactory results in the presence of noise. Single-sensor data using a configuration with depth-varying measurements are able to deliver better results than when constant-depth recordings are used, but the results obtained are not of the same quality when multisensor data are used.

Simulation of surface waves in porous media

SUMMARY We present a novel numerical algorithm for the simulation of seismic wave propagation in porous media, which is particularly suitable for the accurate modelling of surface wave-type phenomena. The differential equations of motion are based on Biot's theory of poro-elasticity and solved with a pseudospectral approach using Fourier and Chebyshev methods to compute the spatial derivatives along the horizontal and vertical directions, respectively. The time solver is a splitting algorithm that accounts for the stiffness of the differential equations. Due to the Chebyshev operator the grid spacing in the vertical direction is non-uniform and characterized by a denser spatial sampling in the vicinity of interfaces, which allows for a numerically stable and accurate evaluation of higher order surface wave modes. We stretch the grid in the vertical direction to increase the minimum grid spacing and reduce the computational cost. The free-surface boundary conditions are implemented with a characteristics approach, where the characteristic variables are evaluated at zero viscosity. The same procedure is used to model seismic wave propagation at the interface between a fluid and porous medium. In this case, each medium is represented by a different grid and the two grids are combined through a domain-decomposition method. This wavefield decomposition method accounts for the discontinuity of variables and is crucial for an accurate interface treatment. We simulate seismic wave propagation with open-pore and sealed-pore boundary conditions and verify the validity and accuracy of the algorithm by comparing the numerical simulations to analytical solutions based on zero viscosity obtained with the Cagniard-de Hoop method. Finally, we illustrate the suitability of our algorithm for more complex models of porous media involving viscous pore fluids and strongly heterogeneous distributions of the elastic and hydraulic material properties.

A general wave decomposition formula for the measurement of normal incidence sound transmission loss in impedance tube

Two types of general methods can be found in the literature for the determination of the normal incidence sound transmission loss (nSTL) of acoustical elements. The first one is based on the transfer matrix (TM) approach, and the second one is based on the wavefield decomposition (WD) theory. From all the techniques proposed in the literature, the general TM methods (two-load or two-source location) are the only methods yielding the exact nSTL of an acoustical element without any assumptions on its symmetry and on the termination (i.e., the load). Except for the case of an anechoic termination, there is no method based on the WD theory which yields exact nSTL. This paper presents a general WD method to measure the exact nSTL of an acoustical element without any assumptions on its symmetry and on the termination. Similar to general TM methods for non-symmetrical elements, four microphones and two loads will be required. As a first validation of the method, symmetrical and non-symmetrical porous materials are investigated. Results are discussed and compared with some existing methods and with the classical two-load method. A perfect agreement is found with the classical two-load method.

Subsonic near-surface P-velocity and low S-velocity observations using propagator inversion

Detailed knowledge of near-surface P- and S-wave velocities is important for processing and interpreting multicomponent land seismic data because (1) the entire wavefield passes through and is influenced by the near-surface soil conditions, (2) both source repeatability and receiver coupling also depend on these conditions, and (3) near-surface P- and S-wave velocities are required for wavefield decomposition and demultiple methods. However, it is often difficult to measure these velocities with conventional techniques because sensitivity to shallow-wave velocities is low and because of the presence of sharp velocity contrasts or gradients close to the earth's free surface. We demonstrate that these near-surface P- and S-wave velocities can be obtained using a propagator inversion. This approach requires data recorded by at least one multicomponent geophone at the surface and an additional multicomponent geophone at depth. The propagator between them then contains all information on the medium parameters governing wave propagation between the geophones at the surface and at depth. Hence, inverting the propagator gives local estimates for these parameters. This technique has been applied to data acquired in Zeist, the Netherlands. The near-surface sediments at this site are unconsolidated sands with a thin vegetation soil on top, and the sediments considered are located above the groundwater table. A buried geophone was positioned 1.05 m beneath receivers on the surface. Propagator inversion yielded low near-surface velocities, namely, 270 ± 15 m/s for the compressional-wave velocity, which is well below the sound velocity in air, and 150 ± 9 m/s for the shear velocity. Existing methods designed for imaging deeper structures cannot resolve these shallow material properties. Furthermore, velocities usually increase rapidly with depth close to the earth's surface because of increasing confining pressure. We suspect that for this reason, subsonic near-surface P-wave velocities are not commonly observed.