Original Wavefield Research Articles

A reciprocity-aware, low-rank regularization for multidimensional deconvolution

A novel factorization-based, low-rank regularization method is introduced to solve multidimensional deconvolution (MDD) problems in the frequency domain. In this approach, each frequency component of the unknown wavefield is represented as a complex-valued square matrix and approximated as the product of a rectangular matrix and its transpose. The benefit of such a parameterization is twofold. First, each frequency matrix is guaranteed to be symmetric. From a physical standpoint, the recovered wavefield inherently adheres to the principle of reciprocity; in other words, the response to a receiver from a virtual source does not change when the two are swapped. Therefore, we refer to this low-rank parameterization as reciprocity-aware in that it respects the underlying physics of wave propagation associated with the MDD solution. Second, the size of the unknown matrix is greatly reduced compared with that of the original wavefield of interest (and halved compared with the existing factorization-based low-rank approximations). We further show that the associated objective function can be successfully optimized using the accelerated proximal gradient algorithm and present a robust strategy to define the initial guess of the solution. Numerical examples on synthetic and field data demonstrate the effectiveness of the proposed method in compressing the retrieved Green’s function while preserving its accuracy.

A fast amplitude preserving three-parameter 3D parabolic Radon transform and its application on multiple attenuation

Seismic wavefields propagate through three-dimensional (3D) space, and their precise characterization is crucial for understanding subsurface structures. Traditional 2D algorithms, due to their limitations, are insufficient to fully represent three-dimensional wavefields. The classic 3D Radon transform algorithm assumes that the wavefield's propagation characteristics are consistent in all directions, which often does not hold true in complex underground media. To address this issue, we present an improved 3D three-parameter Radon algorithm that considers the wavefield variation with azimuth and provides a more accurate wavefield description. However, introducing new parameters to describe the azimuthal variation also poses computational challenges. The new Radon transform operator involves five variables and cannot be simply decomposed into small matrices for efficient computation; instead, it requires large matrix multiplication and inversion operations, significantly increasing the computational load. To overcome this challenge, we have integrated the curvature and frequency parameters, simplifying all frequency operators to the same, thereby significantly improving computation efficiency. Furthermore, existing transform algorithms neglect the lateral variation of seismic amplitudes, leading to discrepancies between the estimated multiples and those in the data. To enhance the amplitude preservation of the algorithm, we employ orthogonal polynomial fitting to capture the amplitude spatial variation in 3D seismic data. Combining these improvements, we propose a fast, amplitude-preserving, 3D three-parameter Radon transform algorithm. This algorithm not only enhances computational efficiency while maintaining the original wavefield characteristics, but also improves the representation of seismic data by increasing amplitude fidelity. We validated the algorithm in multiple attenuation using both synthetic and real seismic data. The results demonstrate that the new algorithm significantly improves both accuracy and computational efficiency, providing an effective tool for analyzing seismic wavefields in complex subsurface structures.

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.

Amplitude-preserving P-/S-wavefield separation with the discontinuous Galerkin method on unstructured meshes

The decoupling of P and S waves is an essential prerequisite for elastic reverse time migration, which effectively reduces crosstalk artifacts, but most wavefield separation algorithms are implemented on uniform rectangular grids. We develop an amplitude- and phase-preserving P- and S-wavefield separation approach on unstructured meshes, which can effectively decompose the original elastic wavefield into P and S wavefields. The isotropic case is considered. With the aid of viscoelastic theory, we choose to attenuate P or S waves and preserve the other wave mode, so as to achieve wavefield decomposition. Viscoelastic wave equations are first reformulated as decoupling wave equations with a selective strong attenuation. We then use the discontinuous Galerkin (DG) method to simulate decoupling P or S wavefield propagation on triangular and tetrahedral meshes. We adopt a quadrature-free DG approach and the arbitrary mesh is mapped into the reference mesh for numerical calculation, where no additional volume and surface integrations are involved. The amplitude and phase information of this vector decomposition agrees with that of the original elastic data. Four numerical examples are used to demonstrate the superior performance of this vector decomposition algorithm. The isotropic example indicates the applicability and correctness of our scheme and the second example displays superiority in handling strong velocity contrasts. The third example exhibits mesh flexibility in dealing with complex structures, such as caves, faults, and undulating surfaces. The last example indicates the effectiveness of our developed algorithm extended to a 3D case.

Free-surface elastic wave injection and wavefield reconstruction with applications to elastic reverse time migration

Accurately injecting and extrapolating multicomponent seismic data is a major step in elastic reverse time migration (RTM) and full-waveform inversion. The theoretical basis of this process is the seismic representation theorem. We propose a numerical approach based on the collocated-grid MacCormack finite-difference (FD) scheme to implement the injection and extrapolation. The MacCormack FD scheme defines wavefield parameters on the same grid. When dealing with the boundary condition in the representation theorem, it is more advantageous than the commonly used staggered-grid FD scheme which defines the stress and velocity on offset grids; thus, interpolations are required to link the boundary values to nearby nodes. In contrast, with the MacCormack scheme, the traction and particle velocity data can be directly injected into the boundary nodes as equivalent concentrated forces, force dipoles, and double couples. The accuracy of this approach is validated by comparing the regenerated wavefield with the original wavefield. Our primary purpose is to develop an injection method for elastic RTM, but this approach also can be used to numerically investigate other interesting issues in the elastic wave injection, including errors and artifacts caused by injection from unclosed boundaries, incomplete boundary values, and uncertain velocity models. In particular, seismic data are usually acquired on the ground surface where the traction-free condition is naturally maintained. We simulate a land acquisition scenario by generating a synthetic multicomponent data set on a traction-free surface. The reconstructed wavefields from velocity data are compared with the original forward-propagated waves to evaluate the accuracy of the boundary treatments. Two multicomponent elastic RTM examples are presented to demonstrate the applications of the proposed injection method. By combining the regenerated receiver-side wavefield and the vectorized imaging condition, the PP, PS, SP, and SS images are calculated.

A Baseline-Free Lamb Wave Damage Localization Method Based on Compressed Sensing

In conventional Lamb wave-based damage detection, estimation is generally performed using the difference between the measured and baseline signals. However, it is difficult to maintain consistent measurement conditions between the measured and baseline signals, likely leading to large measurement errors. This paper proposed a baseline-free Lamb wave-based detection method for locating damage to structures. The damage scattering waves are extracted according to their features in the frequency-wavenumber spectrum and the damage locations are visualized using probability imaging. To reduce the complications of full wavefield acquisition, wavefield data are acquired using piezoelectric transducer (PZT) arrays at sampling rates much lower than the Nyquist frequency, and the original wavefield is reconstructed by applying compressed sensing. Experimental results on an aluminum plate indicate that the proposed method can provide the damage probability at each plate position without requiring a reference signal, suggesting its applicability to damage diagnosis of structures. Moreover, the proposed method achieves a compression ratio of 86.7% compared with the use of Nyquist sampling.

Interpretation of time reversal focal spots based on point-source homogenization

SUMMARYThe time reversal method is based on the backpropagation of seismic waveforms recorded at a set of receivers. When this set forms a closed surface and the elastic properties of the medium are correct, the seismic energy focuses at the source location, creating a focal spot. Such a spot is smooth in space, whereas the original wavefield usually shows a displacement discontinuity at the source. The goal of this paper is to discuss the link between the focal spot and the original source using the concept of homogenized point source. We show that the backpropagated wavefield is equivalent to the sum of two low-wavenumber fields resulting from the homogenization of the original point source. In other words, the homogenized point source is the equivalent force for producing the focal spot. In addition to the demonstration in the general 3-D heterogeneous case, we present some numerical examples in 2-D.

Robust boundary-preserving multi-source least-squares reverse-time migration with combined sparse constraints

Least-squares reverse time migration (LSRTM) can generate a higher-resolution image than other algorithms but with a higher computation cost. Although the conventional source-encoding algorithm can improve the efficiency, it faces the problems of multi-source crosstalk noise and observed-synthetic data mismatching caused by the unpractical assumption of the fixed acquisition setup. In this paper, we propose a hybrid sparsity-constrained multi-source encoding algorithm to overcome the above bottleneck problems. Harmonic wavelets are used as an encoding operator so that the blended multi-source wavefields can be decomposed as each original single-source wavefield based on the orthogonality of trigonometric functions. Due to that the harmonic frequency selection might lead to an insufficient stack of single-frequency components, there would generate background noise on the image, reducing the image quality and convergence rate. Therefore, we apply a combined-sparsity constraint to improve the robustness and convergence rate of this algorithm. Besides, we use the convolution-norm misfit function to eliminate the influence of incorrect source wavelet. Synthetic data examples show the rationality and feasibility of our algorithm.

Elastic reverse time migration angle gathers using a stabilized Poynting vector without zero points within the wave propagation ranges

Generation of PP and PS angle gathers for elastic reverse time migration (ERTM) has important applications in velocity model updating and subsurface parameter inversion. The Poynting-vector-based method provides a cost-effective way to extract angle gathers from ERTM. However, the estimated propagation directions or angles based on Poynting vectors at or near zero points are unstable, and the reflection energy at these points is likely to be placed in an incorrect angle bin. This instability problem severely limits its use in the generation of high-quality angle gathers. To overcome this problem, we add the conventional Poynting vector with an ancillary one of the Hilbert transformed wavefield, a 90° phase-shifted version of the original wavefield. The Poynting vector of the Hilbert transformed wavefield points in the same propagation direction as the conventional Poynting vector, but they have different zero-point distributions. Adding them together allows us to obtain a stabilized Poynting vector which thoroughly eliminates the zero points within the wave propagation range, and thus we obtain a stable algorithm to extract high-quality angle gathers from ERTM. The Hilbert transformed wavefield already exists in the P/S and up/down wavefield decomposition procedure, which has become a necessary step for ERTM to separate overlapping wavefields at reflectors before the computation of Poynting vectors. Results from two synthetic numerical examples indicate the effectiveness of our method in generating high-quality angle gathers for ERTM.

Phase-amplitude-based polarized direct envelope inversion in the time-frequency domain

Large-scale geologic structures with strong contrasts present difficulties for seismic imaging and inversion. Previously, direct envelope inversion (DEI) with ultralow frequency has been proposed to retrieve strong-contrast velocity models. However, the envelope signals lose the instantaneous phase of the original wavefield, resulting in erroneous seismic inversion results. To solve this problem, the instantaneous phase of the wavefield is encoded into envelopes for DEI. In this way, a polarized envelope is created to contain both abundant low frequencies and the instantaneous phase information of the original seismic waveform. Furthermore, to improve the resolution of deep regions underneath the strong reflections, we propose a phase-amplitude-based polarized direct envelope inversion (PA-PDEI) misfit in the time-frequency domain that uses an amplitude factor to emphasize the phase information. This makes it possible to adjust the amplitude influence and boost the low-amplitude seismic signals for deep targets. Tests on a model with a salt body and on field data show that combining the PA-PDEI with the phase-amplitude-based full-waveform inversion (PA-FWI) make it possible to reliably reconstruct velocity models, especially for deep regions, despite large and strong perturbations.

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.

Features of correlation measurements of the parameters of pulsed hyperspectral optical fields using an asymmetric interferometer

Differences in correlation measurements of the parameters of pulsed hyperspectral optical fields using symmetric and asymmetric interferometers are considered. It is shown analytically that the resulting cross-correlation function is sensitive to phase perturbations in the original wave field. The considered setup, which contains a telescopic reflective 4f system of parabolic mirrors in one arm, demonstrates that in the case of an asymmetric interferometer, the presence of aberrations leads to degradation of the reconstructed image, whereas in the case of symmetric interferometers these aberrations do not affect the result.

Modeling the viscoelastic effects in P-waves with modified viscoacoustic wave propagation

Accurate full-waveform seismic modeling is a powerful tool for understanding wave propagation and building subsurface images. However, it can be computationally expensive for viscoelastic media. Viscoacoustic seismic modeling is much cheaper, but at the trade-off of using incomplete physics. We have developed a modified viscoacoustic wave simulation algorithm for modeling the viscoelastic effects of P-waves. The algorithm contains two viscoacoustic forward-modeling steps; the first is the same as the traditional viscoacoustic modeling, whereas the second propagation is generated using a residual error source, which is derived by comparing the viscoacoustic and viscoelastic wave equations in the form of stress-particle velocity formulations. The corrected P-wave particle velocities can be obtained by adding the wavefield from the second simulation step to the original (the first simulation step) viscoacoustic wavefield. Only P-waves are modeled. The overall cost is about twice that of viscoacoustic modeling, but it is significantly less than a viscoelastic propagation because there are fewer calculations, and we can use a coarser grid and larger time steps for the same accuracy. Numerical examples indicate that the P-wave waveforms, after correction, match those from viscoelastic wave modeling better than those from the original viscoacoustic simulation. Our method provides a cost-efficient alternative for approximating the viscoelastic effects in P-wave modeling.

An improved fast converted-wave imaging method

The conventional fast converted-wave imaging method directly uses backward P- and converted S-wavefield to produce joint images. However, this image is accompanied by strong background noises, because the wavefields in all propagation directions contribute to it. Given this issue, we improve the conventional imaging method in the two aspects. First, the amplitude-preserved P- and S-wavefield are obtained by using an improved space-domain wavefield separation scheme to decouple the original elastic wavefield. Second, a converted-wave imaging condition is constructed based on the directional-wavefield separation and only the wavefields propagating in the same directions used for cross-correlation imaging, resulting in effectively eliminating the imaging artifacts of the wavefields with different directions; Complex-wavefield extrapolation is adopted to decompose the decoupled P- and S-wavefield into directional-wavefields during backward propagation, this improves the efficiency of the directional-wavefield separation. Experiments on synthetic data show that the improved method generates more accurate converted-wave images than the conventional one. Moreover, the improved method has application potential in micro-seismic and passive-source exploration due to its source-independent characteristic.

Poynting and polarization vectors based wavefield decomposition and their application on elastic reverse time migration in 2D transversely isotropic media

ABSTRACTWith the progress in computational power and seismic acquisition, elastic reverse time migration is becoming increasingly feasible and helpful in characterizing the physical properties of subsurface structures. To achieve high‐resolution seismic imaging using elastic reverse time migration, it is necessary to separate the compressional (P‐wave) and shear (S‐wave) waves for both isotropic and anisotropic media. In elastic isotropic media, the conventional method for wave‐mode separation is to use the divergence and curl operators. However, in anisotropic media, the polarization direction of P waves is not exactly parallel to the direction of wave propagation. Also, the polarization direction of S‐waves is not totally perpendicular to the direction of wave propagation. For this reason, the conventional divergence and curl operators show poor performance in anisotropic media. Moreover, conventional methods only perform well in the space domain of regular grids, and they are not suitable for elastic numerical simulation algorithms based on non‐regular grids. Besides, these methods distort the original wavefield by taking spatial derivatives. In this case, a new anisotropic wave‐mode separation scheme is developed using Poynting vectors. This scheme can be performed in the angle domain by constructing the relationship between group and polarization angles of different wave modes. Also, it is performed pointwise, independent of adjacent space points, suitable for parallel computing. Moreover, there is no need to correct the changes in phase and amplitude caused by the derivative operators. By using this scheme, the anisotropic elastic reverse time migration is more efficiently performed on the unstructured mesh. The effectiveness of our scheme is verified by several numerical examples.

Attenuation of spatial aliasing by means of stationary-wavelet and slant-stack transforms

Spatial aliasing happens in some seismic data acquisition and affects severely the performance of multichannel data processing and interpretation. Conventional methods (e.g. antialias filtering) remove high-frequency information from data in hand and distort the signals. The procedure we present here produces a unaliased estimation of all frequency components within the original seismic wavefield. The method is benefitted from the advantages of the Wavelet and Radon transforms, for instance the overlap of information between wavelet scales at the same frequency, in decomposing the seismic time series. We demonstrated the performance of the developed dealiasing algorithm in wavelet-Radon domain by applying on synthetic seismic data.

Investigations on flexural wave propagation and attenuation in a modified one-dimensional acoustic black hole using a laser excitation technique

Abstract Acoustic Black Holes (ABHs), as a new type of passive structure for vibration damping enhancement and noise attenuation, have been drawing increasing attentions of many researchers. Due to the difficulty in manufacturing the sharp edges required by the ABH structures, it is important to understand the wave propagation and attenuation process in the presence of damping layers in non-ideal ABHs with a truncated edge. In this paper, an analytical expression of the wave reflection coefficient in a modified one-dimensional ABH is derived and a time-domain experimental method based on a laser excitation technique is used to visualize the wave propagation. In the experimental studies, the flexural waves in the ABH were excited by a scanning pulse laser and measured by a Laser Doppler Vibrometer (LDV). The incident wave and reflected wave were separated from the measured original wave field and the decrease of the wave velocity in the ABH was exhibited. The reflection coefficient was calculated from the ratio of the amplitude of the reflected wave to that of the incident wave for different ABH parameters and different thicknesses of the damping layer. The measured reflection coefficients were used to identify the unknown coefficients in the theoretical formula. The results confirm that there exists an optimal thickness for the damping layer, which leads to the minimum wave reflection. Based on the laser-induced visualization technique and various signal processing and feature extraction methods, the entire process of the wave propagation in a non-ideal one-dimensional ABH structure can be visualized and scrutinized.

Diffraction separation by plane-wave prediction filtering

Seismic data processing typically deals with seismic wave reflections and neglects wave diffraction that affect the resolution. As a general rule, wave diffractions are treated as noise in seismic data processing. However, wave diffractions generally originate from geological structures, such as fractures, karst caves, and faults. The wave diffraction energy is much weaker than that of the reflections. Therefore, even if wave diffractions can be traced back to their origin, their energy is masked by that of the reflections. Separating and imaging diffractions and reflections can improve the imaging accuracy of diffractive targets. Based on the geometrical differences between reflections and diffractions on the plane-wave record; that is, reflections are quasi-linear and diffractions are quasi-hyperbolic, we use plane-wave prediction filtering to separate the wave diffractions. First, we estimate the local slope of the seismic event using planewave destruction filtering and, then, we predict and extract the wave reflections based on the local slope. Thus, we obtain the diffracted wavefield by directly subtracting the reflected wavefield from the entire wavefield. Finally, we image the diffracted wavefield and obtain high-resolution diffractive target results. 2D SEG salt model data suggest that the plane-wave prediction filtering eliminates the phase reversal in the plane-wave destruction filtering and maintains the original wavefield phase, improving the accuracy of imaging heterogeneous objects.

Research Note: Near‐surface layer replacement for sparse data: is interpolation needed?

ABSTRACTNear‐surface problem is a common challenge faced by land seismic data processing, where often, due to near‐surface anomalies, events of interest are obscured. One method to handle this challenge is near‐surface layer replacement, which is a wavefield reconstruction process based on downward wavefield extrapolation with the near‐surface velocity model and upward wavefield extrapolation with a replacement velocity model. This requires, in theory, that the original wavefield should be densely sampled. In reality, data acquisition is always sparse due to economic reasons, and as a result in the near‐surface layer replacement data interpolation should be resorted to. For datasets with near‐surface challenges, because of the complex event behaviour, a suitable interpolation scheme by itself is a challenging problem, and this, in turn, makes it difficult to carry out the near‐surface layer replacement. In this research note, we first point out that the final objective of the near‐surface layer replacement is not to obtain a newly reconstructed wavefield but to obtain a better final image. Next, based upon this finding, we propose a new thinking, interpolation‐free near‐surface layer replacement, which can handle complex datasets without any interpolation. Data volume expansion is the key idea, and with its help, the interpolation‐free near‐surface layer replacement is capable of preserving the valuable information of areas of interest in the original dataset. Two datasets, i.e., a two‐dimensional synthetic dataset and a three‐dimensional field dataset, are used to demonstrate this idea. One conclusion that can be drawn is that an attempt to interpolate data before layer replacement may deteriorate the final image after layer replacement, whereas interpolation‐free near‐surface layer replacement preserves all image details in the subsurface.

Discrete wave equation upscaling

Abstract We present homogenization technique for the uniformly discretized wave equation, based on the derivation of an effective equation for the low-wavenumber component of the solution. The method produces a down-sampled, effective medium, thus making the solution of the effective equation less computationally expensive. Advantages of the method include its conceptual simplicity and ease of implementation, the applicability to any uniformly discretized wave equation in 1-D, 2-D or 3-D, and the absence of any constraints on the medium properties. We illustrate our method with a numerical example of wave propagation through a 1-D multiscale medium and demonstrate the accurate reproduction of the original wavefield for sufficiently low frequencies.