Abstract In mudstones, undercompaction inhibits grain rearrangement, preserving pores and generating anomalous petrophysical responses. These anomalies typically manifest in geophysical exploration as abnormal acoustic and resistivity logs. However, the microstructural differences between undercompacted and normally compacted conditions remain unclear, and the mechanisms linking pore preservation to petrophysical property anomalies are not well understood. Neogene mudstones from the Bohai Bay Basin serve as a case study, integrating well log analysis, experimental research, and digital rock physics methods to investigate pore structure and petrophysical property responses in mudstones under different compaction conditions. Although well log analysis exhibits non-uniqueness, rock physics experiments conclusively demonstrate that these anomalies are related to mudstone undercompaction. Results reveal that undercompacted mudstones exhibit significantly higher porosity (average 20.54%) than normally compacted mudstones (average 14.26%), along with a higher clay content and a well-connected network of interaggregate pores. The primary cause of their low resistivity is the combined effect of this interconnected macro-pore system and high clay content. Density-acoustic transit time crossplots confirm that overpressure in the study area originates from undercompaction rather than fluid expansion or diagenesis. Digital rock simulations quantitatively reveal that porosity has approximately twice the effect of clay content on resistivity reduction, and further demonstrate that undercompacted mudstones generate low acoustic impedance anomalies similar to sandstone reservoirs but are distinguished by significantly higher Vp/Vs ratios, providing a robust seismic discriminant. These findings clarify the pore and mineral controlled origins of petrophysical anomalies in undercompacted mudstones and offer practical criteria for their accurate identification in geophysical exploration, thereby supporting improved reservoir interpretation and drilling hazard assessment.
Abstract The Majiaoba area, on the northwestern margin of the Longmen Mountains within the Sichuan Basin, exposes a complete Silurian–Triassic succession. Silurian–Devonian strata are widespread, Carboniferous–Triassic formations are well developed, and metamorphic or igneous rocks are rarely exposed at the surface. With the exception of the weakly magnetic Lower Triassic Feixianguan Formation, most units are essentially non-magnetic. Despite this, the region exhibits high-amplitude magnetic anomalies (up to ~103 nT) whose genesis remains uncertain, and basin-scale magnetic highs in central-western Sichuan are commonly ascribed to rifting even though their driving mechanisms and lateral continuity are debated. Two high-precision profiles across Majiaoba were acquired, and gravity and magnetic datasets were inverted independently using a minimum-support, adaptive iterative focusing algorithm selected under multi-criteria evaluation, yielding high-resolution images of subsurface structure beneath both lines. Rock-magnetic measurements on Feixianguan samples, combined with the inversion results, reveal a pronounced vertical offset between the principal magnetic source and the surface-exposed reddish mudstone-shale of the Feixianguan Formation. Along profile L2, a closed high-susceptibility body is resolved at ~400–600 m depth. The anomaly pattern is inconsistent with a purely shallow sedimentary origin and is best explained by mafic (diabase) intrusions. At the regional scale, the deep magnetic source is interpreted as a product of sustained activity related to the large-scale magmatic episode along the western margin of the Yangtze Block during the Late Permian, with magmatism in the northern Longmen Mountain segment plausibly persisting into the Early Triassic (Feixianguan time). Considering gravity–magnetic anomaly patterns across northwestern Sichuan, the magmatic system was likely widespread rather than local in extent.
Abstract In Antarctica, more than 98% of the continent is covered by ice sheets. Seismic exploration is recognized as one of the most effective methods for high-resolution imaging of ice sheets and the underlying bedrock. However, conducting active-source seismic surveys remains highly challenging due to the extreme environmental conditions and limited logistical support. Passive seismic exploration, which does not require active sources, has therefore emerged as a promising approach for efficient seismic investigations on Antarctic ice sheets. However, its imaging capability under field conditions still requires direct validation against active-source results. This study provides a field-scale validation of passive seismic imaging by directly comparing passive results with co-located active-source results acquired along a survey line on the ice sheet in the Larsemann Hills, East Antarctica. Two seismic processing workflows are developed for surface-wave and body-wave imaging. For surface waves, least-squares noise suppression and beamforming-based stationary-phase selection are used to improve the quality of reconstructed virtual shot gathers. Multichannel analysis of surface waves (MASW) is conducted on both active and passive datasets, and the resulting dispersion characteristics and inverted shear-wave velocity models show strong consistency. For body-wave imaging, conventional reflection processing techniques used for active-source data are adapted for virtual shot gathers retrieved from passive recordings. The stacked passive seismic profile reproduces the main reflection events observed in the active-source profile, particularly the ice-bed interface and shallow reflections within the bedrock. These comparative results demonstrate that, in the Larsemann Hills region of Antarctica, passive seismic methods provide an environmentally sustainable, cost-effective, and reliable alternative to active seismic surveys for large-scale imaging of ice sheets and shallow bedrock structures.
Abstract Wave-induced local fluid flow at different scales, which contributes most to the wave velocity and attenuation, is determined not only by the compressibility of rock but also by its permeability properties. To investigate the potential relationship between permeability and seismic attributes, a wave propagation theory for fluid-saturated media with infinituple-porosity is extended to include rock textures (inclusions) with scale-dependent permeability. The effects of permeabilities, radii and elastic moduli of inclusions on wave dispersion and attenuation were first investigated. The results show that the stiffness moduli of the inclusions dominate, while for the same moduli, the permeability of the inclusions also exerts an influence at seismic frequencies. To validate the model, the results are compared with laboratory measurements on tight sandstones (1 Hz-106 Hz) and field measurements on marine sediments (50 Hz-400 kHz). The agreement confirms that the proposed model can be used to analyse and interpret the observed wave propagation and attenuation in a wide frequency range.
Abstract Scholte waves are surface waves that propagate along a fluid-solid interface and arise from the interference of P- and SV-waves. The conventional dispersion equations for surface waves are typically derived under the assumption of horizontal interfaces, and cannot be applied to inclined strata. The presence of inclined multi-layered seabed introduces substantial mathematical challenges in deriving dispersion equations for Scholte waves. To address this, a natural stratum coordinate system and a seismic wave propagation coordinate system are established. Coordinate rotation transformation matrices are derived to convert the seismic wave equations and boundary conditions between these two coordinates. In this natural stratum coordinate system, the dispersion equations of Scholte waves are derived with the presence of inclined multi-layered seabed. In addition, the dispersion curves and wavefields of Scholte waves are calculated for theoretical models with several inclined interfaces. The numerical results demonstrate that both the dipping angles and azimuths substantially influence the dispersion characteristics of Scholte waves. The effects of dipping angles on the dispersion curves are larger than those of azimuths. These derived equations can be further applied to high precision inversion of dispersion curves and to predict the S-wave velocities, dipping angles, and azimuths of shallow seabed with complex fluid-solid interfaces.
Abstract Poststack seismic acoustic impedance inversion is of vital significance in reservoir prediction. Nonetheless, conventional poststack seismic inversion approaches linearize the relationship between seismic response and acoustic impedance, which restricts their adaptation to weak reflection interfaces only. The ensemble smoother with multiple data assimilation (ESMDA) has emerged as a powerful tool for generating a credible set of prior ensemble members by matching simulated seismic responses to available observations. Yet, in standard implementations of ESMDA, single-trace inversion methods often encounter spatial discontinuities and instability. This is mainly because of the independence among prior ensemble members at each trace. Generally, ESMDA consists of two parts: generation of prior ensemble members and update of ensemble members. To enhance the spatial continuity of the ESMDA inversion results, stratigraphy information is incorporated into the prior ensemble members through stratigraphy-guided sequence Gaussian simulation. In this way, spatial continuity is introduced into the prior ensemble members, and inversion is performed to update these models. Tests on both synthetic and field data demonstrate the effectiveness of the proposed ESMDA inversion scheme. It can provide high-resolution inversion results with robust spatial continuity.
Abstract 3D time-domain electromagnetic (TEM) simulations typically employ sequential time stepping, requiring early-time calculations before late times, reducing parallelizability compared to frequency-domain methods. To address this issue, a decomposition method is developed to parallelize forward modeling of TEM time-stepping for arbitrarily complex waveforms. Based on the theory of survey decomposition and scale matching principle, our approach computes the time-domain responses of each time channel separately in parallel. Computational efficiency comes from two ideas. First, TEM fields diffuse at later times so that each time channel can be stepped at a step length adapted to its temporal scale. Second, exact simulation of the entire on-time waveform is not always necessary – early off-time channels only need a narrow portion of the pulse width before the turn-off to be precisely described; late off-time channels need the time integration of the on-time waveform to account for the total energy transmitted. Our approach first empirically determines a characteristic step length (δt) for a particular time channel by exploiting those two properties. Then, a step-off response is obtained by stepping at a constant δt; next, the discrete impulse response is calculated by differentiating the step-off response with an impulse width δt. Finally, the discrete impulse response is convolved with the effective portion of the transmitter waveform discretized by δt. By experimenting with a variety of TEM waveforms and many random models, we obtain a set of empirical parameters for δt and the effective pulse width for practical use. 3D TEM examples using the VTEM and HELITEM waveform have demonstrated that multi-scale and parallel time-stepping consumes a small fraction of what would be required by the conventional sequential method. The decoupled computation enables accelerated wide-band TEM time-stepping modeling in massively parallel environments by eliminating inter-channel dependencies similar to the frequency-domain methods.
Abstract Full waveform inversion (FWI) provides high-resolution velocity models by exploiting the complete information embedded in seismic waveforms. However, the advantages of FWI come with cycle skipping when sufficient information is not available to build an accurate initial model. To integrate non-seismic data into seismic FWI, this study establishes a direct conversion method, leveraging the correlation between multi-physical models. The diffusive nature of electromagnetic (EM) data makes it an independent, low-cost information source that complements seismic low-frequency data. Unlike joint inversion approaches, our flexible FWI framework initializes FWI using EM resistivity inversion results. This novel approach extracts low-frequency skeleton structure from electromagnetic data to construct an initial FWI model. By design, it avoids strong assumptions about petrophysical correlation and structural alignment and does not require complex multi-physics objective functions, resulting in a straightforward implementation. Our three-stage implementation begins by inverting EM data to generate a blurred resistivity image. Pixels in this image are then clustered into models, each assumed to have a constant velocity. To prevent cycle skipping, a preliminary FWI is performed to roughly fit seismic data by determining velocities for these models. Finally, a full-bandwidth FWI refines the preliminary velocity model, which already incorporates EM and low-frequency seismic information, aiming for the highest possible resolution. Applied to the Marmousi model, our approach outperformed conventional seismic-only FWI methods in terms of data fitting and computational performance. Using EM to warm-up FWI is theoretically similar to employing an initial velocity model by smoothing the true model, because EM surveys can be regarded as a low-pass filter of the subsurface structure, as EM surveys physically average subsurface structures. Our findings underscore the importance of incorporating non-seismic data in seismic imaging and offer a robust workflow for joint multi-physical data acquisition and analysis.
Abstract 3D seismic data acquisition, processing, and analysis have become increasingly prevalent in modern geophysical exploration. However, various types of coherent noise, especially linear noise, can cause serious interference and effective noise reduction techniques are critical. Most deep learning-based denoising methods rely on converting 3D seismic volumes into 2D slices and then reassembling them, which disrupts the spatial continuity and structural integrity of the data. To address this issue, we propose a novel linear noise suppression network capable of directly processing 3D seismic volumes, effectively preserving spatial features and continuity. Our proposed approach, the 3D Deformable Convolution-Transformer U-Net (3D-DCTU) network, integrates 3D deformable convolutional layers with a Transformer architecture. The deformable convolutional layers adaptively adjust sampling locations through learned offsets to capture spatial variations in seismic data, while the Transformer's self-attention mechanism captures long-range dependencies. This combination enables both the extraction of local features and the modeling of global information. To enhance perceptual quality of the denoised results, the network incorporates the gradient difference loss and mean squared error as optimization targets. We evaluated the proposed network using both synthetic data and field data from a region in western China. The results suggest promising denoising performance, particularly for linear noise attenuation, indicating a viable technical pathway for 3D seismic data denoising.
ABSTRACT Geologic feature analysis from seismic images plays a vital role in geologic assessments, resource exploration, natural disaster prediction and assessment, and carbon capture and storage. Traditional deep learning-based geologic feature recognition methods often suffer from a high demand for labeled training samples and limited model generalization. In recent years, foundational models pretrained on large-scale unlabeled data through self-supervised learning have gained significant attention in the computer vision community owing to their strong generalization and robustness, and initial attempts have been made to extend such models to seismic image analysis. Nevertheless, most existing efforts rely on 2D slices and focus primarily on general seismic tasks, with little emphasis on the specialized optimization needed for identifying diverse geologic features. This study proposed a dual pretraining semi-supervised framework to develop SS-UNETR, a foundational model that leverages a 3D Swin Transformer backbone. SS-UNETR was designed to concurrently segment multiple geologic features in seismic volumes through multitask learning. The training procedure consisted of two sequential stages. In the first stage, the encoder of SS-UNETR was pretrained on 75,651 field seismic images via three proxy tasks (image inpainting, rotation prediction, and contrastive learning) to learn underlying representations and basic patterns without explicit labeling. In the second stage, SS-UNETR was fine-tuned using 2860 synthetic seismic images with corresponding geologic annotations to adapt the model to geologic feature segmentation. Experimental results demonstrated that SS-UNETR excelled in analyzing geologic features of seismic images compared to methods based on seismic attributes and neural networks trained for specific tasks. SS-UNETR achieved effective simultaneous identification of multiple geologic features, even when these features exhibited highly similar seismic responses, for such features as faults, channels, and caves. Furthermore, SS-UNETR outperformed conventional deep learning models without pretraining and general-purpose seismic foundational models in few-shot geologic feature identification tasks. These results indicated that SS-UNETR can serve as a valuable tool in the industry to enhance the efficiency and accuracy of geologic interpretation, reduce the burden of manual labeling, and reduce the carbon footprint of model training.
ABSTRACT SS-wave (SH–SH and SV–SV waves) data can more accurately calculate fracture parameters, image gas-cloud areas, and invert S-wave velocity and density. However, due to high acquisition costs and a low signal-to-noise ratio, high-quality SS-wave seismic data are limited in field seismic exploration. In the laboratory, previous experimental studies on the fracture parameters and geophysical properties of reservoirs were mostly based on 1D ultrasonic S-wave transmission, which failed to study amplitude variation with offset in SS-wave reflection gathers. A new physical modeling method on a solid surface was systematically introduced, and 2D reflection wavefields of P–P, SV–SV, SH–SH, and converted P–SV waves acquired using the method were analyzed. Unlike traditional physical modeling on water surface that only captures PP-wave reflections, this new method acquired nine-component seismic reflection gathers featuring stronger amplitudes, although it resulted in higher noise levels. F–K filtering effectively improved the signal-to-noise ratio (SNR) of multiwave data. Notable observations included the following: the amplitude of P-wave reflections varied considerably across different wavefields, whereas frequency changes were relatively small; the signal types in the SV–SV wavefield were diverse, including P-wave, SV-wave, and converted waves; in contrast, the SH–SH wavefield primarily consisted of SH-wave reflections with higher amplitudes. The amplitude variation with angle characteristics of SH-wave reflections influenced by porosity parameters was also discussed and analyzed, revealing that these characteristics varied across different porosity ranges. The experimental results could contribute to advancing the application of SH-wave data in reservoir characterization and enhancing the understanding of the characteristics of various signals within SV- and SH-wave-excited wavefields.
ABSTRACT Conventional acoustic impedance inversion methods have long faced technical bottlenecks such as inaccurate wavelet estimation and strong dependence on initial models. Although existing deep learning approaches can partially alleviate these problems, they often compromise model simplicity and training efficiency while introducing new challenges such as limited generalizability and heavy reliance on labeled data. To overcome these limitations, a lightweight inversion framework that tightly integrates physics-driven and data-driven paradigms was proposed. The physics-driven component adopted a neural network architecture largely consistent with traditional modeling processes, enabling direct optimization of physically meaningful parameters through backpropagation, thereby avoiding the construction of excessively complex inverse operators. Meanwhile, regularization methods were introduced to enforce geologic prior knowledge (that is, the “layered geologic model” assumption) on the network parameters, improving the spatial continuity of the reconstructed impedance models. The data-driven component used an enhanced 2D U-Net integrated with class activation mapping to generate accurate reference models from sparse well-log data. Tests on synthetic and field data sets demonstrated the advantages of the proposed method: (1) a physically interpretable network design; (2) strong robustness to noise and reduced dependence on training data; (3) higher accuracy and better spatial continuity compared with conventional and purely data-driven methods. The method provided a new perspective for addressing long-standing challenges in seismic impedance inversion.
ABSTRACT Full-waveform inversion (FWI) is a high-resolution seismic inversion technique widely used in oil and gas exploration. Traditional FWI uses the l2 norm measurement to minimize the misfit between observed and predicted seismic data. However, when the background velocity is inaccurate or the seismic data lack low-frequency components, the conventional FWI suffers from cycle skipping, leading to inaccurate inversion results. To address this issue, a multiscale structural similarity index measure (SSIM) objective function was introduced for FWI. Anisotropic total variation regularization with the lp quasi-norm (ATpV) was also incorporated to further improve the accuracy of FWI. Multiscale SSIM extracts multiscale structural features of seismic data in time and space dimensions. These features can reduce the risk of cycle skipping and improve the stability of FWI. Additionally, ATpV applies structural constraints to the velocity gradients, which helps suppress artifacts and preserve the sharp boundaries of geologic formations. Automatic differentiation (AD) was used to efficiently and stably optimize this novelly introduced FWI objective function. Experiments on synthetic and field seismic data demonstrated that the proposed method accurately characterizes complex subsurface velocity structures, even when the background velocity is crude, the data lack low-frequency components, or contain noise.
Abstract Accurate prediction of the boundaries and filling degree of fault-controlled fracture–caves is critical for the effective exploration of fault-controlled reservoirs. An optimized method for predicting the filling degree of such reservoir bodies was developed, targeting the Maokou Formation in the Luzhou area of the southern Sichuan Basin. Based on one-dimensional wave equation forward modeling, gradient structure tensor attributes were computed with standard deviations σ1 = 0.6 and σ2 = 1.6, using λ2 as the eigenvalue to delineate boundaries. A threshold-based and normalized fault-controlled fracture–cave probability volume was generated as a constraint for prestack seismic prediction. AVA (Amplitude Versus Angle) forward modeling verified the feasibility of filling degree prediction using prestack data, and original gathers were optimized. Prestack simultaneous inversion yielded elastic parameter volumes. A rock physics interpretation panel was constructed to classify filling degrees, enabling volume estimation for different categories. Results demonstrate that the gradient structure tensor effectively identified boundaries, with a threshold >0.19 indicating fracture–cave development. AVA responses varied by filling degree: fully filled bodies showed no AVA features, partially filled exhibited Class I or II, and unfilled showed Class III characteristics. The Vp/Vs ratio and P-impedance served as effective discriminators. Fracture-caves primarily occurred on both sides of the faults, with unfilled and partially filled ones concentrated near large faults and intersecting smaller ones, comprising 48.81% of the total volume. The prediction results were validated by dynamic verification wells, achieving an 83.3% concordance rate. The integration of multi-source seismic data and a progressive, constraint-driven prediction strategy enables a significant improvement in the accuracy and efficiency of predicting the filling degree in fault-controlled fracture–cave systems.
Abstract Full waveform inversion (FWI) is a powerful technique for building high-resolution subsurface models, but it is fundamentally ill-posed. Traditional FWI cannot quantify uncertainties arising from noise, modeling errors, the intrinsic high nonlinearity, and other sources. In Bayesian inference, Markov chain Monte Carlo (MCMC) methods facilitate direct sampling from the posterior distribution to resolve these concerns. However, MCMC approaches are slow to converge and inefficient for exploration in high-dimensional spaces due to the complex posterior distribution in FWI. In this study, we develop an uncertainty quantification framework for Bayesian FWI based on underdamped Langevin dynamics, which adds momentum and inertia to sample trajectories, allowing for more efficient posterior exploration. We implement two splitting schemes, Strang splitting and OBABO splitting (named after the sequence of operator applications in underdamped Langevin dynamics), to discretize the stochastic differential equations (SDE). These splitting methods divide the Langevin dynamics into simpler components, each of which can be solved more accurately, and can then be recombined. This reduces discretization errors and improves stability, allowing more efficient sampling. We validate the proposed approach through numerical experiments on two examples. The results show that our methods can effectively explore high-probability regions of the posterior distribution, enabling reliable estimation of posterior means, variances, and marginal probability densities for uncertainty quantification. The Strang splitting scheme exhibits a slightly faster convergence rate and lower computational overhead. Overall, at an acceptable computational cost, our methods achieve rapid convergence, provide robust uncertainty quantification, and yield posterior statistics that are largely insensitive to the choice of initial model.
ABSTRACT One of the primary tasks of electrical imaging logging in oil-based mud (OBM-EIL) is fracture evaluation. At present, there was no research on the fracture aperture calculation of OBM-EIL similar to that of electrical imaging logging in water-based mud. First, numerical modeling was used to analyze the fracture responses of OBM-EIL under the influence of multiple parameters such as the background formation resistivity, fracture resistivity, fracture aperture, and fracture dip angle. Then, a strategy for determining the fracture aperture was proposed, which included establishing multiparameter response databases, adopting a stepping strategy, and developing calculation models of fracture apertures based on backpropagation neural networks. Finally, several tests were conducted to assess the effectiveness of fracture aperture calculation models. The results indicated that fracture apertures in the range of 0.01–100 mm could be accurately calculated, and the determination coefficient could reach 97.1%. The fracture resistivity and the fracture dip angle were essential parameters for the fracture aperture calculation. The lack of background formation resistivity and gap distance has a slight effect on the fracture aperture calculation. These findings will greatly enhance the researchers’ confidence in the fracture aperture calculation using the proposed method, even when some of these parameters might be acquired through inversion or machine learning methods.
Abstract Density interfaces are discontinuous boundaries characterized by significant density contrasts within Earth's interior, such as the seafloor, sedimentary basement, and the Moho. Gravitational observations can be used to infer the geometry of these undulating interfaces. Current and future satellite missions provide global high-quality gravity, gravity-gradient, and gravity-curvature data, enabling local- to global-scale studies and mitigating the limitations of terrestrial measurements. Density interface inversion is intrinsically non-linear, and conventional linearization inevitably introduces bias, particularly in regions with pronounced interface relief. To address this limitation, non-linear forward and inverse models under a spherical approximation are considered, and a higher?order regularization framework is developed to approximate the nonlinear solution. Synthetic tests based on theoretical models, together with applications to Moho recovery beneath the Tibetan Plateau, demonstrate superior performance relative to conventional Tikhonov regularization, indicating that higher-order contributions are non-negligible when the interface undulations are moderate. The application also reveals limitations of gravity-curvature Moho inversion, as the refined gravity-curvature disturbances are strongly influenced by short-wavelength signals associated with imperfect stripping reductions, resulting in weak sensitivity to Moho geometry. Overall, the proposed higher-order regularization provides a robust and effective framework for nonlinear density interface inversion and offers improved algorithmic support for further applications.
ABSTRACT Seismic wavefield simulation serves as an essential step for high-resolution seismic imaging. Traditional methods like the finite-difference method (FDM) require dense grid meshing for complex velocity models. In 3D frequency-domain wavefield simulation, FDM needs to invert a large impedance matrix. This process is extremely slow and consumes substantial computational resources. As an emerging mesh-free paradigm for forward modeling, physics-informed neural networks (PINNs) adopt physical principles as constraints to optimize the network and have demonstrated significant potential in simulating seismic wavefields for complex models. However, hindered by the inherent spectral bias of neural networks, PINNs struggle to capture the high-frequency spatial oscillations of complex-valued seismic wavefields for the Helmholtz equation, leading to decreased accuracy and slow convergence. To overcome this challenge, a new form of the Helmholtz equation parameterized by amplitude and unwrapped phase for PINNs was proposed. Consequently, the frequency-domain wavefield was represented by amplitude and unwrapped phase, which was significantly smoother, resulting in a significantly improved convergence in training PINNs. Numerical experiments on 2D and 3D complex models demonstrated that, compared with the conventional Helmholtz equations based on complex-valued wavefield representation, the proposed method significantly enhanced convergence speed and wavefield simulation accuracy, even when simple fully connected networks were used.
ABSTRACT For seismic characterization of reservoirs, it is necessary to extract the time-frequency characteristics hidden in time-domain data for effective use in processing and interpretation. In the time-frequency analysis algorithm of the chirp-modulated W-transform (CW-transform), the spread of the analysis kernel was the primary factor controlling the resolution of the power spectral density. A multiresolution CW-transform was proposed to improve the localization of the power spectral density and extend the capability of time-frequency algorithms in delineating reservoirs and thin layers. Multiplication of CW-transforms with different resolutions was also proposed to produce a highly localized and well-resolved power spectral density. The resultant spectrum can be interpreted as the geometric mean of the spectra estimated with various scaling factors. Compared with the arithmetic mean, this approach exhibits superior robustness and noise resistance, thereby achieving higher resolution and greater energy concentration. The application of statistical properties of the time-varying spectrum was investigated to indicate the presence of gas-bearing reservoirs. Experiments were conducted to demonstrate that the multiresolution CW-transform can successfully estimate a time-frequency spectrum with high energy concentration, showing great potential for application in the seismic characterization of complex reservoirs.
ABSTRACT Reverse time migration (RTM) requires numerically solving the partial differential wave equation because analytical solutions are infeasible. A significant challenge in numerical methods arises from inaccuracies in derivative approximation, making the numerical wave velocity frequency-dependent and causing numerical dispersion. It is an unphysical artifact that degrades modeling and imaging, particularly at higher frequencies and over time. Avoiding numerical dispersion requires finer spatial grids, which substantially increase computational costs to achieve high-resolution imaging results. The modified nearly analytical discretization (MNAD) method reduces numerical dispersion by incorporating additional analytical relations through simultaneous numerical solutions of the wavefield and its spatial gradient fields, using them in higher-order derivative approximations, and improving spatial derivative estimation via energy conservation optimization. MNAD was introduced for RTM in large-scale studies, where leveraging compact stencils and coarser spatial and temporal grids enables high-resolution imaging with substantially lower computational and memory costs compared to conventional finite-difference (FD) methods. Furthermore, adjoint-state imaging was enhanced with a novel data boundary condition interpolation using MNAD gradient fields, mitigating aliasing effects in recovered images from data recorded at half the Nyquist rate, enabling imaging with fewer sources or receivers, and alleviating acquisition costs. Synthetic experiments validated the method’s performance in modeling and imaging on coarser grids than FD methods and in maintaining stability over longer propagation times. Furthermore, the MNAD-based RTM application to ocean-bottom seismometer (OBS) data in a large-scale study confirmed its capability to achieve high-resolution images with reduced computational costs. Finally, imaging with data sampled at half the Nyquist rate highlighted the potential of the proposed approach for minimizing acquisition costs without sacrificing resolution and suffering from aliasing. These findings affirmed MNAD as a robust and efficient alternative to FD methods for large-scale, high-resolution imaging and offer significant advantages in computation, storage, and acquisition efficiency.
