Efficient least-squares imaging of a complex mountainous area in China using LsmGANs

1 Introduction

High-precision seismic imaging is essential for revealing complex subsurface structures (Miller et al., 1998). Seismic migration, the core technique for achieving this goal, repositions reflected wavefields recorded at the surface to their true subsurface impedance interfaces and thereby establishes a spatial correspondence between seismic responses and geological structures. Early migration methods, grounded in diffraction summation, provided foundational structural information for interpretation. Driven by the increasing demand for high-resolution imaging in complex geological environments, the technology has progressed from Kirchhoff migration (Bleistein and Gray, 2001) and beam migration (Yang and Zhu, 2018) through one-way wave equation migration (Naghadeh and Riahi, 2013) to reverse-time migration (RTM) (Whitmore, 1983).

Kirchhoff migration is an efficient ray-based method that handles steeply dipping structures (Gray et al., 2001; Hokstad, 2000). However, it suffers from amplitude singularities near caustics and struggles to accurately model multipathing wavefields in complex geology. While efficient traveltime compression can alleviate the computational burden of large traveltime tables in 3D prestack Kirchhoff migration (Alkhalifah, 2011), it does not directly mitigate the migration artifacts and uneven illumination caused by multipathing. Beam migration can image multipathing arrivals without dip limitations (Yang et al., 2018), but its primary limitation is the inaccurate simulation of finite-frequency effects in regions with strong lateral velocity contrasts (Gray, 2005; Yang et al., 2022b). One-way wave equation migration is constrained by the directional propagation approximation and fails to reliably image steeply dipping reflectors (Zhang et al., 2003). In contrast, RTM uses a two-way propagator that places no fundamental restrictions on velocity structure, reflector dip, or wavefield type, and it is widely deployed in complex settings such as overthrust faults, subsalt sediments, and buried hills (Liu et al., 2011; Yang et al., 2021).

Mathematically, the migration methods described above are essentially adjoint operators of the forward modeling operator. Consequently, when acquisition geometry is irregular, frequency bandwidth is limited, or subsurface illumination is uneven, these methods rarely yield satisfactory images. Building on the fundamental concept of least-squares seismic inversion introduced by Bamberger et al. (1982), Lailly (1983) and Tarantola (1984) formulated seismic migration as a linear inverse problem, establishing the theoretical framework of least-squares migration (LSM). Compared with conventional adjoint migration, LSM effectively suppresses acquisition footprints caused by irregular geometry, attenuates migration artifacts, enhances amplitude fidelity, and produces higher-quality images.

Since its inception, LSM has been extended considerably, progressing from Kirchhoff migration (Wang et al., 2013) to one-way wave-equation migration (Kuehl and Sacchi, 2002) and RTM (Dai and Schuster, 2013; Yang et al., 2019), from acoustic media (Sun et al., 2016) to anelastic anisotropic media (Ren et al., 2017), and from conventional geophone data (Kaplan et al., 2010) to distributed acoustic sensing data (Chen Ning et al., 2019). More recently, LSM has also been extended to pseudo-inverse extended Born formulations for true-amplitude imaging and migration velocity analysis, which better account for lateral heterogeneity (Alali et al., 2020). LSM can be formulated in either the data domain or the image domain. Data-domain methods iteratively update the reflectivity to minimize the misfit between synthetic and observed data, avoiding explicit construction of the Hessian matrix (Kühl and Sacchi, 2003). However, without proper preconditioning, this iterative inversion remains computationally expensive and converges slowly. Image-domain methods directly compute the Hessian matrix and apply its inverse to the migrated image, but this incurs substantial storage costs (Yang et al., 2022a). Consequently, although LSM improves imaging resolution and amplitude fidelity, its high computational and memory demands remain major obstacles for large-scale applications.

As a data-driven approach, deep learning (DL) (LeCun et al., 2015) with multi-layer neural architectures can construct nonlinear mappings by mining large volumes of training data (Araya-Polo et al., 2019; Yang and Ma, 2019). In recent years, it has been applied to LSM to enhance imaging efficiency and accuracy. Cheng et al. (2020) used a convolutional neural network (CNN) to automatically identify faults and imposed this structural information on the LSM gradient to constrain the image of complex-structures; Liu et al. (2020) investigated the physical meaning of feature maps and designed a CNN architecture for sparsity-promoting LSM; Vamaraju et al. (2021) adopted mini-batch and adaptive moment estimation to implement CNN-based least-squares reverse-time migration with reduced computational time; Kaur et al. (2020) and Huang et al. (2022) employed cycle-consistent generative adversarial networks (GANs) to enhance image quality; Torres and Sacchi. (2022) combined traditional least-squares reverse-time migration with DL to avoid explicit regularization and obtain accurate reflectivity in only a few iterations. Sun et al. (2023) approximated the inverse Hessian by training a GAN to learn a mapping from RTM images to true reflectivity models and the trained network is then applied to the RTM images to improve imaging quality. More recently, Alfarhan et al. (2025) applied a learned Hessian to full waveform inversion (FWI), training a deep neural network to learn a deblurring mapping between Hessian-blurred and unblurred FWI gradients, so that the network acts as an approximate inverse-Hessian operator and improves the FWI update direction.

Supervised DL methods offer a highly promising route to enhancing image quality while reducing computational cost by learning a mapping that approximates the inverse of the Hessian matrix. Most existing DL-based LSM approaches rely on conventional GAN architectures, which can, however, produce inaccurate reflector positioning and migration artifacts. To overcome these limitations, we develop a refined GAN framework tailored for image-domain LSM, termed LsmGANs. The proposed framework has been validated on synthetic data, where it effectively suppresses migration artifacts and produces high-quality reflectivity images.

Although such DL-based LSM methods have achieved remarkable success on synthetic models, their application to field data remains challenging due to the domain gap between training and target datasets, noise contamination under complex near-surface conditions, and the difficulty of efficiently handling three-dimensional data. To address these issues, this study uses a complex mountainous area of the northeastern Sichuan Basin as a test site and systematically evaluates the performance of LsmGANs on real seismic data. Several key practical issues are identified and clarified, providing guidance and reference for the large-scale deployment of DL-based LSM in similar exploration settings.

2 Methods

2.1 Principle of conventional LSM

Imaging seeks to reconstruct the subsurface reflectivity model m from recorded data d. As shown in Equation 1, this relationship can be linearized and expressed in matrix formwhere L is the forward modeling operator. Conventional migration uses adjoint forward modeling L T as the imaging operator to yield an approximation m mig, as shown in Equation 2:

d = Lm , (1)

m mig = L T d . (2)

While adjoint imaging correctly positions most reflections, it does not preserve true reflection amplitudes and suffers from limited resolution. LSM addresses this by treating imaging as an inverse problem, aiming to minimize an objective function. LSM implementations are broadly classified into data-domain LSM and image-domain LSM.

As shown in Equation 3, data-domain LSM directly inverts Equation 1 by minimizing the data misfit:where d obs denotes the observed data, d cal denotes the simulated record. The problem is typically solved using iterative optimization algorithms. Each iteration requires one forward simulation and one migration. As shown in Equation 4, the update can be expressed aswhere m i + 1 is the updated model, m i is the model at iteration i, α is a step-length parameter, and ∇ J m i is the gradient of the data-misfit objective function. Although data-domain LSM rigorously honors the observed data and can, in principle, yield the most accurate reflectivity estimate, its computational cost is prohibitive for large-scale problems, as each iteration involves expensive wave-equation solves.

J m = 1 2 d cal − d obs 2 = 1 2 Lm − d obs 2 , (3)

m i + 1 = m i − α ∇ J m i , (4)

Image-domain LSM offers a computationally cheaper alternative by working entirely in the image domain. It is based on the relation between the adjoint image in Equation 2 and the desired reflectivity, which involves the Hessian matrix H = L T L. As shown in Equation 5, the ideal least-squares solution can be expressed as

m = L T L − 1 L T d = H − 1 m mig . (5)

Compared with data-domain LSM, image-domain LSM circumvents the repeated wave-equation simulations and step-length searches required in iterative data-domain inversion. This makes the approach more efficient and adaptable for seismic imaging applications.

2.2 Improved GAN architecture

A GAN consists of a generator G and a discriminator D that compete in an adversarial training process: the generator G learns to produce samples that resemble real data, while the discriminator D learns to distinguish these generated samples from actual ones. This competition drives both components to improve until they reach a stable equilibrium.

The generator G synthesizes samples that emulate the distribution of the input data, while the discriminator D evaluates their authenticity by outputting a probability between 0 and 1, where 1 denotes high similarity to genuine samples and 0 indicates a clear distinction. Trained in an adversarial manner, the generator G learns to capture the underlying data distribution and produces new samples from random noise, whereas the discriminator functions as a binary classifier. The generator optimizes its parameters to maximize the probability of the discriminator being deceived. This adversarial process reaches a unique equilibrium when the generator G reproduces the true data distribution and the discriminator D converges to a constant output of 0.5. As shown in Equation 6, the alternating optimization of both components is expressed as:where E · represents the expectation operator; p d a t a x is the distribution of real samples; p z is the distribution of the input noise; x denotes a real sample; z is the random noise vector; and G z is the synthetic sample produced by the generator. The architecture of the GAN is illustrated in Figure 1.

min G max D V D , G = E x ∼ P d a t a x log ⁡ D x + E x ∼ P z log 1 − D G z , (6)

During training, the network weights Θ G and Θ D are updated alternately using the stochastic gradient descent method. As shown in Equation 7, at each iteration, the discriminator D is updated using the gradient calculated by maximizing the misfit function, as follows:

g GAN D = ∇ Θ D 1 n ∑ i = 1 n log ⁡ D x i + log 1 − D G z i . (7)

In Equation 8, the generator G is updated using the gradient calculated by minimizing the misfit function, as follows:where n is the mini-batch size.

g GAN G = ∇ Θ G 1 n ∑ i = 1 n log 1 − D G z i , (8)

It is worth noting that the standard GAN framework described above is trained solely with the adversarial loss, which encourages the generator to match the distribution of the target data but does not explicitly enforce pixel-wise consistency. For the task of seismic reflectivity inversion, the “real data” in the standard GAN framework shown in Figure 1 corresponds to the true reflectivity model. Relying only on the adversarial loss may lead to structural distortions and inaccurate recovery of weak illumination zones. To further constrain the generator and force the output to match the true reflectivity x i, an L 1 reconstruction loss is incorporated. As shown in Equation 9, the gradient of the generator with respect to the L1 reconstruction loss is calculated by minimizing the misfit function, as follows:where x i denotes the true reflectivity corresponding to the i-th input, and · 1 represents the L 1 norm. Consequently, the overall optimization objective for the generator G combines the adversarial loss and the reconstruction loss. As shown in Equation 10, the total misfit function for updating the generator is formulated as:where λ is a weighting hyperparameter that balances the contributions of the adversarial loss and the reconstruction loss. In this study, λ is empirically set to 100 to ensure structural fidelity during early training while allowing the adversarial loss to refine the textural details. By minimizing this composite objective, the network learns to produce high-fidelity reflectivity models with suppressed migration artifacts and enhanced weak-illumination recovery.

g rec G = ∇ Θ G 1 n ∑ i = 1 n G z i − x i 1 , (9)

g total G = g GAN G + λ g rec G , (10)

Both the generator G and the discriminator D in a GAN are multi-layer neural networks built from convolutional and fully connected layers. Their architectures are shown in Figures 2a,b. The generator G follows a U-shaped architecture consisting of an input layer (layer 1), an encoder block (layers 2 and 3), a transform block (layers 4 and 5), a decoder block (layers 6 and 7), and an output layer (layer 8). The input layer contains a 7 × 7 convolutional layer, normalization, and rectified linear unit (ReLU) activation. The larger kernel size in the input layer can enrich the basic features extracted from the input image. Both the encoder and decoder blocks contain two 2-D 3 × 3 convolutional layers followed by normalization and ReLU activation. In the encoder block, we use the convolutional layers with a stride of 2 to downsample the size of the input. The number of channels starts at 32 and doubles after each downsampling operator. In the decoder block, we use the up-convolution with a stride of 2 to upsample feature maps back to the same size as the input image. The number of channels in the decoder block starts at 128 and halves after each upsampling layer. The transform block contains one residual structure. The residual structure consists of two 3 × 3 convolutional layers: the first is followed by the normalization and ReLU activation, whereas the second is followed by the normalization operator. The number of channels and the feature map dimensions remain within the residual structure. The output layer uses the Tanh activation to overcome the centrosymmetry problem associated with Sigmoid activation. In the U-shaped architecture, the encoder block first gradually reduces the spatial dimension of the input image. Then, the extracted high-resolution features from the encoder block are mapped from one domain to another through the transform block. Finally, the decoder block and output layer recover the spatial dimensions of transformed feature maps and output the reflectivity image. The discriminator D is composed of one side of the U-Net and different colors denote different operators.

In the conventional GAN framework, downsampling is typically performed by a convolution with a stride of 2. This operation inevitably discards a portion of the feature maps and can lead to mislocated reflection events and lateral discontinuities when applied to migrated seismic images. To overcome these limitations, we introduce several architecture modifications tailored to the characteristics of seismic data, resulting in an improved framework termed LsmGANs (Figure 3). The key innovations are threefold. First, we replace the strided convolution with a combination of a stride-1 convolution followed by max-pooling. The stride-1 convolution is designed to preserve input information, while max-pooling retains the strongest local responses during downsampling. This design is motivated by the geophysical observation that the peak amplitude of a seismic wavelet, rather than its full waveform, carries the critical travel-time information that controls reflector positioning. By preserving local maxima during downsampling, max-pooling is expected to help anchor reflection events to their correct spatial locations and reduce the sensitivity of the network to small spatial variations of reflectors. Second, we employ multiple stacked transform layers in the generator, rather than a single one. Each transform layer incrementally expands the receptive field, which may enable the network to better reconcile information across adjacent traces and enforce lateral coherence. This progressive refinement is motivated by the iterative processing workflows commonly used in conventional seismic imaging. The stacked transform layers are designed to enhance the lateral continuity of the predicted reflectors, an improvement that is difficult to achieve through max-pooling modifications alone. Third, we introduce skip connections between corresponding encoder and decoder levels. In image-domain LSM, the migration velocity model provides an essential low-wavenumber constraint that should be preserved across the full network depth. Skip connections provide direct pathways that are intended to prevent this physical prior from being attenuated in deeper layers, while simultaneously alleviating the gradient-vanishing problem. Consequently, the network is expected to learn a more robust nonlinear mapping and yield improved image quality.

These three modifications form a complementary design. Max-pooling is intended to provide correctly located events, multiple transform layers are designed to promote lateral continuity, and skip connections are incorporated to help retain the velocity prior throughout the network. Although the individual contributions of each of the components are not separately quantified, a comparison between the full LsmGANs framework and the original GAN baseline indicates that the combined integration of these modifications leads to improved event placement, enhanced lateral coherence, and greater training stability.

2.3 Image-domain LSM with LsmGANs

The Image-Domain LSM solution in (5) requires computing the Hessian inverse H − 1, a key operator for seismic resolution. However, the computation is expensive, especially for large-scale data. We therefore propose the LsmGANs network, a novel GAN-based architecture that learns a complex mapping from a standard migration image to the true reflectivity. Equation 5 states that the true reflectivity m is obtained by applying the inverse Hessian H − 1 to the migrated image m mig. The Hessian encodes the spatially variant illumination and blurring imposed by the acquisition geometry and velocity model. In a supervised learning framework, training a network to map m mig to m using a large number of paired examples drives the learned mapping to approximate H − 1, because the network aims to approximate H − 1, which serves as the ideal operator for this transformation. However, under limited illumination or complex acquisition geometries, H is often singular or ill-conditioned, where a strictly unique H − 1 may not exist. The supervised learning framework therefore learns a pseudo-inverse or regularized inverse to mitigate these issues. The velocity model is supplied as an additional input to inform the network of the kinematic factors that dominate the structure of H, thereby constraining the approximation to physically admissible solutions. The network learns a data-driven approximation to the action of a regularized inverse Hessian within the distribution of the training data. Compared with a traditional GAN, LsmGANs incorporates three key adaptations tailored for seismic imaging: max-pooling for structure-preserving downsampling, multiple stacked transform layers for enhanced lateral continuity, and skip connections for improved gradient flow. Together, these modifications address the dislocation and discontinuity artifacts that commonly arise in conventional GAN-based image-domain LSM. As shown in Equation 11, under the DL framework, the LSM solution can be written aswhere m mig = L T d obs is the input migrated image, m is the reflectivity model, and G , ; Θ G denotes the parameterized network with trainable weights Θ G. The inputs to network G are the migrated image m mig and migration velocity V 0, and the output is the reflectivity image m. Here, V 0 provides essential structural and kinematic constraints that help produce more reliable reflectivity estimates. Figure 4 illustrates the workflow of learning a data-driven approximation to the action of a regularized inverse Hessian with the proposed LsmGANs framework.

m = G m mig , V 0 ; Θ G , (11)

Our workflow takes two primary geophysical inputs: the migrated seismic image m mig, which encodes the raw reflectivity structure, and the initial velocity prior V 0, which provides low-wavenumber constraints for regularization. These inputs are fed into the generator G, which learns the nonlinear mapping G m mig | V 0 to produce a synthetic reflectivity estimate m *. This synthetic estimate is then paired with the velocity prior to form the conditional input m * , V 0, while the discriminator D evaluates this pair against the ground-truth counterpart m , V 0, where m denotes the reference reflectivity model. Conditioned on V 0, D learns to distinguish between real and synthetic data, ensuring that the generated outputs remain geophysically consistent with the prior velocity model. Through adversarial optimization, the generator minimizes the discriminator’s ability to distinguish synthetic from true LSM results, while the discriminator maximizes its classification accuracy.

3 Synthetic data tests

3.1 Construction of the 2D training dataset

DL methods require a large and diverse training set to learn effective features. To train the LsmGANs network, we construct a training dataset from the Marmousi model, which contains complex faults with steep dips and strong velocity variations that challenge conventional imaging approaches. The migration velocity models are obtained by smoothing the true models with a Gaussian filter. The network labels are the true reflectivity models calculated as R = v t r u e − v m i g / v m i g. We compute the seismograms by directly solving the acoustic wave equation with a free-surface boundary condition at the top and absorbing conditions at the lateral and bottom edges. The migrated images are computed by RTM. Rather than adding synthetic random noise, we introduce realistic noise extracted from the field data to the migrated images. Specifically, ambient noise is isolated from the pre-first-arrival window of the raw field shot gathers, where no coherent reflection signal is present. This field-derived noise preserves the natural spectral characteristics, spatial correlation, and amplitude statistics of the acquisition environment, providing a more representative challenge for the network than artificially generated Gaussian or uniform noise. The extracted noise is scaled to achieve a target signal-to-noise ratio and added to the noise-free synthetic migrated images. For reproducibility without access to the original field data, comparable results can be obtained using Gaussian white noise with a signal-to-noise ratio (SNR) of approximately 5 dB. Each resulting noisy migrated image, together with its corresponding velocity model and the noise-free reference, constitutes one training pair. During training, the network takes the noisy migrated image and migration velocity as inputs, and outputs the true reflectivity. The network weights are updated to learn a data-driven approximation to the action of a regularized inverse Hessian within the training-data distribution. Although GANs are conventionally trained on large datasets, the LsmGANs framework converges stably with a moderate number of training samples. This data efficiency arises from the constrained nature of the image-domain LSM task. Unlike unconditional image generation where a GAN must synthesize diverse content from random noise, image domain LSM requires the network only to learn a residual correction to the input migrated image, which already captures the dominant structural framework of the subsurface. The generator does not invent geological structures from scratch but refines an existing geophysical image by removing spatially variant blur and compensating for illumination loss. This substantially reduces the learning burden and lowers the data volume required for stable convergence. For the Marmousi model, we generate 405 pairs of 128 × 128 input and 128 × 128 output patches, divided into 305 pairs for training, 50 for validation, and 50 for testing.

3.2 Testing on the marmousi model

Both the GAN baseline and the proposed LsmGANs share the same training and inference configuration. Specifically, the framework is implemented using the PyTorch DL library on a single NVIDIA Quadro RTX 5000 GPU with 16GB of memory. The model is trained using the Adam optimizer with a learning rate of 1 × 10 − 4 and a batch size of 1. In this section, the Marmousi model (Figure 5a) is used to evaluate the performance of the proposed LsmGANs. The model is discretized on a grid 344 × 1536 with an interval of 10 m. Each 128 × 128 patch therefore covers approximately one-third of the depth extent and 1–12th of the lateral extent, providing sufficient local structural context such as multiple reflectors and faults for the network to learn dip-invariant features. The migration velocity shown in Figure 5b is obtained by applying a Gaussian filter with a window size of 25 × 25m to the true model. Figures 5c,d show the original RTM image and the RTM image with noise added, respectively. The true reflectivity model is shown in Figure 6a and is used as the reference for comparison. After 100 training epochs, the GAN result in Figure 6b exhibits visibly curved deep reflectors and limited spatial resolution, indicating that the network has not yet learned a stable mapping. The LsmGANs prediction in Figure 6c shows noticeable improvement over the GAN counterpart, yet several reflection events remain laterally misaligned by approximately 2–3 sample points, residual migration noise is pervasive, and the gas reservoir in the upper-left region remains poorly resolved. These artifacts indicate that 100 epochs are insufficient for the network to learn a stable and physically accurate mapping. After 300 training epochs, the GAN result (Figure 6d) still suffers from limited resolution, noticeable background noise, and visibly curved deep reflectors that deviate from the true model. In contrast, the LsmGANs result (Figure 6e) shows marked improvement: most events are correctly positioned, lateral continuity is restored, noise is substantially suppressed, and the deep reflections appear sharper and structurally consistent with the reference. The LsmGANs prediction at 300 epochs provides the closest match to the true reflectivity in terms of both structural fidelity and amplitude preservation.

Figure 7 presents a progressive comparison across the three highlighted regions. The noisy RTM image (first column) serves as the input baseline. After 100 epochs, the GAN prediction (second column) shows only marginal improvement over the RTM input, while the LsmGANs prediction (third column) demonstrates noticeable gains yet still exhibits residual misalignment and noise. After 300 epochs, the GAN result (fourth column) remains limited in resolution and structural accuracy. In contrast, the LsmGANs result (fifth column) consistently resolves both flat reflectors and steeply dipping faults with markedly improved structural continuity and amplitude fidelity, approaching the true reflectivity (sixth column) more closely than any other prediction. This progressive comparison demonstrates that LsmGANs not only converges faster but also achieves a substantially higher performance ceiling than the standard GAN under identical training conditions.

Table 1 provides a comprehensive quantitative evaluation on the Marmousi model, comparing the reconstruction performance of the proposed LsmGANs against both the conventional baseline (RTM) and the standard GAN framework. The conventional RTM method serves as a physical imaging baseline, yielding an SSIM of 0.3738 and a PSNR of 18.43 dB. In contrast, the proposed LsmGANs framework demonstrates superior performance across all tested epochs. Notably, as training progresses from 100 to 300 epochs, the LsmGANs model exhibits a steady optimization trajectory, improving the SSIM from 0.7655 to 0.9143 and the PSNR from 30.27 dB to 34.11 dB. When compared to the standard GAN after 300 epochs of training, which achieved an SSIM of 0.6578 and a PSNR of 28.00 dB, LsmGANs surpasses it by a substantial margin, with a relative enhancement of approximately 39.0% in SSIM and a 6.11 dB gain in PSNR. Most critically, compared to the physical RTM baseline, the proposed LsmGANs at 300 epochs achieves a staggering relative SSIM improvement of nearly 144.6% and an absolute PSNR gain of 15.68 dB. Crucially, all reported SSIM and PSNR values are computed on the final fully reconstructed continuous sections, rather than on individual overlapping patches, thereby preventing any metric inflation caused by highly correlated subsets. These quantitative metrics unequivocally validate the efficacy of the LsmGANs architecture in refining spatial and structural details and suppressing noise, substantially outperforming both traditional physics-based approaches and standard deep generative networks.

TABLE 1. Quantitative comparison of RTM, GAN, and LsmGANs on the Marmousi model.

4 Field data tests

4.1 Overview of the complex mountainous survey area

We apply the trained LsmGANs to a field dataset acquired in the M1 block of the Tongnanba gas field, located in a structurally complex mountainous region of the northeastern Sichuan Basin, China. The Tongnanba field is situated at the intersection of the Micangshan thrust belt and the Dabashan arcuate nappe belt. Subjected to Yanshanian and Himalayan orogenies, the strata were intensely deformed into a dominant ENE-trending elongated anticline belt with pervasive fracturing. The primary reservoir intervals are the Xu-2 and Xu-4 Members of the Upper Triassic Xujiahe Formation, which consist of thin interbedded sandstones and shales deposited in braided-river delta systems, with individual sand bodies frequently close to or below the seismic tuning thickness.

These geological characteristics impose several significant challenges for seismic migration imaging. First, the target reservoir is buried at considerable depth, causing substantial attenuation of high-frequency seismic energy and broadening of the effective source wavelet prior to migration. Second, the tightly folded anticline architecture produces steeply dipping flanks where bedding dips frequently exceed 40°. Combined with rugged surface topography and strong lateral velocity contrasts across the fault zones, these steep dips generate highly uneven subsurface illumination. In the deeper portions of the anticline and adjacent to the Mujialiang fault zone, the conventional RTM image exhibits weak amplitude fidelity, reduced spatial resolution, and pronounced migration artifacts. Third, the thin interbedded reservoir layers fall below or near the tuning thickness, making individual sand bodies difficult to resolve in the migrated image and posing a direct challenge to the resolution enhancement capability of any imaging method. These intertwined difficulties establish the survey area as a rigorous and representative test for evaluating the amplitude compensation and resolution improvement offered by LsmGANs.

4.2 Construction and reconstruction of the 3D testing dataset

The original field dataset is a 3D volume of dimensions 1701 × 681 × 281, corresponding to depth, inline, and crossline directions, respectively. To adapt the data to our 2D LsmGANs network and construct the test set, the volume is sliced along the crossline direction into a series of inline sections, each of dimensions 1701 × 681. Each extracted 2D section is then subdivided into 128 × 128 patches to match the network input size. Along each section, this results in 26 patches in the depth direction and 10 patches in the inline direction, yielding 260 patches per section. With a total of 281 inline sections, the entire 3D volume is partitioned into 281 × 260 patches of 128 × 128 samples, expanding the original volume into a comprehensive set of 2D patches for validation. After the network predicts each 128 × 128 patch, the local outputs are first reassembled into complete 2D inline sections using an overlapping tiling strategy. The sliding step between adjacent patches is set to be smaller than the patch size, so that neighboring patches share a region of overlap in both the depth and inline directions. Each 2D inline section is reconstructed from 26 patches in the depth direction and 10 patches in the inline direction, yielding 260 patches per section. During reconstruction, a two-dimensional Hanning window is applied to each patch before placement, and a pixel-wise accumulator records the weighted sum and the total weight at each grid point. The standard stride between adjacent patches is fixed at 64 pixels, which corresponds to a 50% overlap in both the depth and inline directions. To accommodate the volume boundaries without resorting to zero-padding, an adaptive boundary strategy is adopted. The starting indices of the terminal patches are adaptively shifted so that the final patch aligns exactly with the data boundaries, yielding reduced strides of 37 pixels in the depth direction and 41 pixels in the inline direction. The final value at each pixel is obtained by dividing the accumulated weighted sum by the total weight. This edge-weighted averaging scheme effectively blends the overlapping regions, eliminates boundary discontinuities within each section. Once all 281 inline sections are reconstructed, they are stacked sequentially along the crossline direction to form the complete 3D volume.

4.3 Analysis of field application results

We apply the LsmGANs framework to a land field dataset. Figure 8a shows the 3D migration velocity model, and Figure 8b displays the 3D prestack depth migration (PSDM) volume. Following the slicing and tiling strategy established earlier, both volumes are partitioned into numerous 128 × 128 patches. These velocity and PSDM patches are then fed into the trained network, and the resulting local reflectivity images are reassembled to produce the final output.

We selected inline 100 and crossline 300 as representative examples. Figures 9a,b show the corresponding slices extracted from the 3D migration velocity model. These two orthogonal profiles together capture the spatial velocity heterogeneity and structural complexity of the deep reservoirs.

As shown in Figure 9a, the inline 100 profile exhibits strong lateral velocity variations and steeply dipping strata. A prominent low-velocity anomaly of approximately 4.5 km/s appears at the right flank of the shallow section, likely associated with near-surface gas chimneys or unconsolidated fluid-bearing formations. Beneath this zone, tilted sedimentary beds with irregular interfaces suggest active tectonic deformation. At depth, an undulating high-velocity basement exceeding 6.0 km/s poses a challenging imaging target due to strong wavefield distortions and energy attenuation caused by the overlying complex overburden. In contrast, the crossline 300 profile shown in Figure 9b displays relatively horizontal and gently undulating layers with smooth lateral velocity gradients. A localized low-velocity lens is observed near the shallow left edge, while the mid-to-deep section transitions gradually from sedimentary strata to a deep high-velocity regime, lacking the steep structures seen in the inline profile. This laterally continuous setting provides a favorable benchmark for evaluating the algorithm’s ability to preserve lateral coherence in flat-layered geological environments. The combination of these two profiles, one structurally complex and the other gently stratified, enables a rigorous assessment of the proposed framework across contrasting geological scenarios, demonstrating its robustness for practical 3D imaging applications.

Figure 10a shows the inline 100 slice extracted from the PSDM volumes. To match the coverage of the original PSDM, a mask was applied to the network prediction. The corresponding LsmGANs network prediction is shown in Figure 10b. Although the predicted result exhibits sharper reflection events and moderately improved resolution relative to the input, a clear asymmetry emerges: most reflections on the right are barely recovered, whereas those on the left are well reproduced. The dominant frequency content is comparable on both sides, ruling out a frequency-dependent cause. A closer inspection reveals a systematic difference in structural orientation: the left side is dominated by upward-dipping reflectors, while the right side consists predominantly of downward-dipping structures. Because the Marmousi training dataset contains predominantly upward-dipping configurations, the network fails to generalize to the downward-dipping patterns encountered in this field section.

This directional bias reveals a fundamental limitation of the original training strategy. The Marmousi model, despite its structural complexity, does not uniformly represent all possible reflector orientations. A network trained on this dataset therefore inherits an implicit preference for the dominant dip directions and systematically underperforms on structures dipping in underrepresented directions. To address this, we augmented the original Marmousi training patches by applying a combination of classical geometric transformations, including horizontal flipping, vertical flipping, and rotations of 90 °, 180 °, and 270 °, thereby expanding the dataset by a factor of six. These augmentations were originally adopted to enrich the training set and improve the network’s generalization to field data with diverse structural dips. However, it should be noted that vertical flipping may violate the fundamental geological principle of increasing velocity with depth. While the augmented dataset proved effective in our experiments, we recommend that future applications prioritize physically admissible transformations and retain only horizontal flipping to ensure consistency with subsurface velocity-depth relationships. The network was first trained from scratch on this augmented dataset for 300 epochs, allowing it to learn a dip-invariant feature representation. It was then fine-tuned on an additional set of patches representing 10% of the original training data, which were selected to contain densely spaced reflection events. This fine-tuning step functions as targeted transfer learning: the augmented pretraining provides a strong, dip-robust initialization, and the subsequent fine-tuning adapts this general representation to the specific amplitude and frequency characteristics of the field data, which often differ from those of the synthetic training data. The final predictions after data augmentation and fine-tuning are presented in Figures 10c,d, where Figure 10c corresponds to the GAN baseline and Figure 10d to the proposed LsmGANs framework. Compared with the initial result in Figure 10b, the fine-tuned LsmGANs result in Figure 10d successfully recovers the right-side reflectors to a considerable extent, demonstrating the effectiveness of the fine-tuning strategy. Figure 10c serves as the baseline for comparison. The GAN result still suffers from limited resolution and noticeable residual noise. In contrast, the LsmGANs result in Figure 10d achieves markedly improved structural continuity, more balanced amplitude characteristics, and more effective noise suppression, suggesting the superior imaging capability of the proposed framework under the same training conditions.

We selected these three representative regions spanning the shallow, intermediate, and deep sections to systematically evaluate imaging performance across distinct geological settings and SNR conditions. The red box, located at the shallow-to-intermediate depth, features steeply dipping structures with relatively strong reflection amplitudes, providing a benchmark for assessing spatial resolution and steep-dip imaging accuracy. The yellow box, at intermediate depth, contains laterally continuous layered strata and tests the algorithm’s ability to maintain coherent reflection events across adjacent traces. The blue box, situated in the deep section, represents a challenging environment where wavefield attenuation and absorption of high-frequency energy result in weak signals and severe amplitude loss, serving as a stringent test of the algorithm’s robustness in recovering faint reflections under low signal-to-noise ratio conditions.

Figure 11 presents enlarged views of the three regions marked in Figures 10a–d. From left to right, the columns correspond to PSDM, initial LsmGANs, GAN after fine-tuning, and LsmGANs after fine-tuning. The vertical arrangement of these regions enables a direct comparison of the spatial alignment and amplitude characteristics across depths. The initial LsmGANs prediction (second column) offers a modest improvement over the PSDM image by enhancing subsalt illumination and attenuating migration noise, although the right-side reflectors remain poorly recovered. The fine-tuned GAN result (third column) still suffers from limited resolution and noticeable residual noise. In contrast, the fine-tuned LsmGANs result (fourth column) exhibits substantial improvements across all three regions. Structural conformability is restored in the steeply dipping shallow section, lateral continuity is reinforced in the intermediate layered strata, and weak deep reflections are recovered with visually more balanced amplitude characteristics. This consistent enhancement from the shallow to the ultra-deep section confirms the effectiveness of the proposed augmentation and fine-tuning strategies and demonstrates the superior imaging capability of the LsmGANs framework under the same training conditions.

Figure 12a shows the crossline 300 slice extracted from the 3D PSDM volume. A mask was applied to the network prediction to match the coverage of the original PSDM. The initial LsmGANs prediction in Figure 12b achieves improved resolution relative to the PSDM input. The fine-tuned GAN result in Figure 12c serves as a baseline for comparison and still exhibits limited resolution with residual noise. In contrast, the fine-tuned LsmGANs result in Figure 12d delivers the best overall performance. Because the crossline section is dominated by nearly flat reflectors, the improvement appears modest in structural terms. Nevertheless, subtle enhancements are visible, with several reflection events appearing sharper and more clearly defined, confirming the superior imaging capability of the proposed framework under the same training conditions.

The crossline profile is dominated by subhorizontal, laterally continuous layers at intermediate depths. Three regions, marked by red, yellow, and blue boxes, are selected at different lateral positions along the crossline direction to evaluate horizontal continuity and stability within this depth window. The red and yellow boxes target robust reflection events at depths of approximately 7–11 km, where the SNR is favorable and the primary challenge lies in preserving lateral coherence. The blue box extends to roughly 9–13 km to further assess the algorithm’s sensitivity to weak signals where attenuation becomes significant. These laterally distributed regions complement the structural complexity analysis performed on the inline profile, providing a focused examination of the algorithm’s ability to maintain reflector continuity and suppress migration artifacts in relatively flat, depth-varying stratigraphic settings.

The zoomed-in views for three highlighted regions in Figures 12a–d are shown in Figure 13. The first column shows the PSDM images, and the second column presents the initial LsmGANs predictions, which improve upon PSDM by suppressing migration noise. The third column displays the GAN results after fine-tuning, serving as a baseline for comparison. The third column displays the GAN results after fine-tuning, serving as a baseline for comparison. The fourth column shows the LsmGANs results after fine-tuning, which achieve improved resolution and more effective noise suppression over both the initial prediction and the fine-tuned GAN baseline. In the red and yellow regions at intermediate depths, the enhanced LsmGANs result preserves lateral continuity of the strong reflection events with greater structural conformity. In the deeper blue region, where the original signal is intrinsically weak, faint reflections are recovered with improved amplitude balance and spatial resolution beneath the salt. This consistent enhancement across laterally distributed regions, each posing distinct imaging challenges, confirms the robustness and depth-adaptive capability of the proposed method.

To evaluate the computational efficiency of the proposed method, we compared the training and inference times of the GAN baseline and LsmGANs on both synthetic and field data. Both models were trained for 300 epochs. For the Marmousi dataset, the training set was augmented by a factor of six through geometric transformations, as detailed in Section 4.3. As summarized in Table 2, LsmGANs consistently exhibits slightly lower training and inference times than the GAN baseline, while both methods achieve rapid inference on the order of a few seconds per image or slice. These results indicate that the improved imaging quality of LsmGANs is not obtained at the expense of computational efficiency.

TABLE 2. Comparison of training and inference times between the GAN baseline and the proposed LsmGANs.

5 Conclusion and discussion

Conventional image-domain LSM is constrained by the prohibitive computational cost of explicitly constructing and inverting the Hessian matrix. To address this bottleneck, this study introduced LsmGANs, a GAN-based framework that learns a data-driven approximation to the action of a regularized inverse Hessian within the training-data distribution. The proposed architecture incorporates three key modifications to the standard GAN: max-pooling for structure-preserving downsampling, multiple stacked transform layers for enhanced lateral continuity, and skip connections for improved gradient flow and velocity-prior preservation. These adaptations effectively mitigate the dislocated reflection events and lateral discontinuities commonly observed in conventional GAN-based LSM methods.

The LsmGANs framework is explicitly conditioned on the migration velocity model, which provides an essential low-wavenumber constraint that helps guide the solution toward geologically plausible reflectivity estimates. On the Marmousi model, quantitative evaluation confirms the effectiveness of the proposed design: at 300 training epochs, LsmGANs achieves an SSIM of 0.9143 and a PSNR of 34.11 dB, substantially outperforming both the standard GAN and the conventional RTM baseline. Numerical experiments further demonstrate that LsmGANs consistently improves structural continuity, amplitude fidelity, noise suppression, and spatial resolution on the widely used Marmousi synthetic model. Field-data tests in a structurally complex mountainous area of the northeastern Sichuan Basin validate the practical viability of the method. Compared with the GAN baseline under identical data augmentation and fine-tuning conditions, LsmGANs achieves superior structural continuity and noise suppression on both inline and crossline field profiles. Transfer learning experiments confirm that the Marmousi-trained network retains a useful degree of generalization when applied to unseen datasets, and data augmentation with geometric transformations combined with targeted fine-tuning. The geometric data augmentation effectively enhances the network’s robustness to diverse structural dips, while the targeted fine-tuning further improves its adaptability to the amplitude and frequency characteristics of field data. Together, these strategies enhance the network’s overall generalizability and robustness during transfer to field data. Notably, the data efficiency of the framework arises from the constrained nature of the image-domain LSM task: the network learns a residual correction to the input migrated image rather than synthesizing geological structures from scratch, substantially reducing the training data required for stable convergence. Once trained, the network generates a high-quality reflectivity image in seconds, offering a computationally efficient alternative to iterative inversion schemes. This case study provides a useful reference for the large-scale deployment of DL-based LSM in complex exploration settings.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Tariq Alkhalifah, King Abdullah University of Science and Technology, Saudi Arabia

Reviewed by: Swarajya Madhuri Rayavarapu, Gayatri Vidya Parishad College of Engineering, India

Zhongming Liu, Jiangxi Normal University, China