Distributed acoustic sensing (DAS) has emerged as a transformative technology for high-resolution seismic exploration. However, compared to conventional geophysics, the vertical seismic profile (VSP) data obtained by DAS has a relatively low signal-to-noise ratio (SNR). This limitation stems primarily from the transient fiber deployment in casing operations, which leads to suboptimal fiber-ground coupling and creates characteristic zigzag noise patterns that degrade signal fidelity. This study presents a feature picking principal component analysis (FPPCA) framework for adaptive zigzag noise suppression. The method includes four stages: power spectral density estimation with multilevel spectral smoothing to preserve critical mid-low frequency components, design of frequency-adaptive hanning windows targeting harmonic sidelobe suppression, bandwidth parameter optimization guided by dominant frequency localization to prevent spectral aliasing during feature extraction, and PCA based noise separation using cumulative variance thresholds derived from localized zigzag features. Validation dataset comprising one synthetic and one field DAS-VSP datasets demonstrates the framework’s ability to maintain broadband signal integrity while achieving spectrally consistent noise attenuation.

The presence of internal multiples significantly impacts seismic data processing, particularly in conventional imaging, where it introduces notable image artifacts. Existing methods for suppressing internal multiples primarily focus on the prestack data domain, often involving trade-offs between accuracy and computational cost. Some strategies have been proposed to mitigate image artifacts in the image domain, offering a more cost-effective solution. To address the high cost of suppressing internal multiples in the prestack data domain and eliminate the impact of artifacts in the image domain, this paper investigates imaging of data containing both primaries and internal multiples using reverse time migration based on adaptive matching filter. This study begins with a theoretical framework for reverse time migration of data containing both primaries and internal multiples, followed by an analysis of the generation of image artifacts in the imaging process. A strategy for artifact suppression based on adaptive matching filter is then proposed. Numerical and field data examples are applied to validate the feasibility and effectiveness of the proposed method. This method can eliminate image artifacts in the image domain instead of high-cost suppression of internal multiples in the prestack data domain, offering significant practical advantages.

Seismic high-resolution processing plays a pivotal role in enhancing the accuracy of velocity analysis and subsurface imaging in seismic exploration. However, both low and high frequencies are often attenuated due to Earth filtering and diffraction effects. Conventional compressed sensing-based frequency expansion methods, while effective in broadening the seismic frequency spectrum, are typically applied to individual traces. This approach, however, overlooks lateral spatial relationships, resulting in increased noise susceptibility and poor lateral continuity. In this study, we propose a multi-trace compressed sensing frequency expansion technique that is constrained by geological features. This approach improves lateral spatial relationships during inversion, thereby enhancing noise robustness and lateral continuity. To accurately estimate time shifts between adjacent noisy traces, we introduce an improved noise-robust dynamic time warping (Noise-DTW) algorithm. These time shifts, combined with a differential operator, are used to construct a geological feature operator for multi-channel inversion. The estimated geological feature operator is incorporated as a regularization term into the L1-norm-based sparse inversion objective function, which is solved using the Alternating Direction Method of Multipliers (ADMM). The effectiveness of the proposed method is validated through both modeling and real seismic data. The results demonstrate significant improvements in seismic data resolution, while also enhancing the lateral continuity of the imaging results.

The theory of compressed sensing (CS), offering a novel perspective beyond the Nyquist-Shannon sampling theorem, has garnered significant attention from geophysicists. It enables the reconstruction of ideal high-density seismic signals using cost-effective random undersampling, provided that the seismic signals can be sparsely represented. However, the near-surface noise challenges the sparse representation of seismic reflection signals. To suppress near-surface noise and achieve a better sparse representation of seismic reflection signals within the CS framework, this article develops a method for CS-based near-surface noise suppression, centered on 3-D conical Radon transform (CRT) combined with the CS theory. First, this study introduces the 3-D CRT into seismic signal processing, exploiting the “conical surface” characteristic of near-surface noise in 3-D CMP records. A method for extracting seismic signals and suppressing noise in the 3-D conical Radon domain is proposed. Subsequently, an adaptive filtering strategy is integrated with the 3-D CRT to develop a method for adaptive near-surface noise suppression in the 3-D conical Radon domain. Finally, combining CS theory, adaptive filtering strategies, and the 3-D Radon transform, this article proposes a method for adaptive near-surface noise suppression in the 3-D conical Radon domain based on CS. Numerical examples demonstrate that the integration of the CS random undersampling framework, adaptive filtering strategies, and the 3-D CRT significantly mitigates the impact of near-surface noise on seismic reflection signal reconstruction within the CS framework. This effectively addresses the issue of the 3-D CRT’s unsuitability for near-surface noise suppression within the CS framework.

The presence of internal multiples significantly impacts seismic data imaging. Due to the subtle differences in velocity and energy between primaries and internal multiples, the task of eliminating internal multiples poses a significant challenge, often leading to inaccuracies or high costs using current methods. Considering that the internal multiples are also real reflections from interfaces containing valuable structural information, reverse time migration (RTM) can be used to handle their imaging but with artifacts introduced. Therefore, to avoid the high costs associated with internal multiple suppression and to address the impact of internal multiples on imaging, our study focuses on the joint imaging of primaries and multiples using RTM, aiming to mitigate the influence of internal multiples artifacts during the imaging process. First, the generation of the artifacts in joint imaging of primaries and internal multiples is analyzed. Then, an imaging constraint strategy is proposed for the images of separated upgoing and downgoing wavefields based on similarity analysis to address the image artifacts, according to the different behavior of the separated images. Finally, the feasibility and effectiveness of the proposed method are verified by numerical examples. The proposed method can realize the elimination of image artifacts in the image domain instead of suppression of internal multiples in the prestack data domain, thereby reducing calculation costs. In addition, the migration of internal multiples may offer potential advantages in imaging complex cases.

Ground-roll attenuation is a challenging seismic processing task in land seismic surveys. The ground-roll coherent noise with low frequency and high amplitude seriously contaminates the valuable reflection events, corrupting the quality of seismic data. The transform-based filtering methods leverage the distinct characteristics of the ground roll and seismic reflections within the transform domain to attenuate the ground-roll noise. However, the ground roll and seismic reflections often share overlaps in the transform domain, making it challenging to remove ground-roll noise without attenuating useful reflections. We propose to apply a conditional diffusion denoising probabilistic model (c-DDPM) to attenuate the ground-roll noise and recover the reflections efficiently. We prepare the training dataset using the finite-difference modeling method and the convolution modeling method. After the training process, the c-DDPM can generate the clean data given the seismic data as the condition. The ground roll obtained by subtracting the clean data from the seismic data might contain some residual reflection energy. Thus, we further improve the c-DDPM to allow for generating the clean data and ground roll simultaneously. We then demonstrate the feasibility and effectiveness of our proposed method by using the synthetic data and the field data. The methods based on the local time-frequency (LTF) transform and U-Net are also applied to these two examples for comparison with our proposed method. The test results show that the proposed method performs better in attenuating the ground-roll noise from the seismic data than the LTF and U-Net methods.

Land seismic data are often obscured by noise, severely affecting the accuracy of subsequent seismic imaging and interpretation. Dictionary learning (DL) is an effective method for noise suppression. However, finding a fast DL method that is suitable for weak signals and can suppress multi-scale strong noise is still a hot topic. In this paper, we introduce a noise suppression method that combines DL with fast adaptive empirical mode decomposition (FABEMD). We leverage the advantages of FABEMD in multi-scale signal decomposition, along with the efficient sparse representation capabilities of DL, to achieve noise suppression for low signal-to-noise ratio seismic signals. We group bi-dimensional intrinsic mode functions based on their cross-correlation coefficients and train dictionaries for components using the sequential generalization K-means method, enhancing computational efficiency and adaptability. Numerical examples using both synthetic and field data validate the practicality and versatility of the proposed method, indicating its improved performance in denoising compared tof-xEMD, BEMD, and traditional DL methods.

As seismic exploration focuses on deep and ultra-deep hydrocarbon targets, seismic data are characterized by weak reflection signals and extremely low signal-to-noise ratio (SNR). Although weak reflections help delineate deep geological structures, the low SNR presents challenges for traditional denoising methods. We propose a high-dimensional adaptive progressive seismic denoising (APSD) method to enhance the SNR of deep weak reflection signals. Instead of using a global noise variance, we estimate local noise variances using a 3-D Laplacian mask based on the local characteristics of seismic data at different locations for nonstationary seismic signals. It employs local noise variance to calculate a Gaussian bilateral kernel function to estimate high-amplitude noise in the time domain and low-amplitude noise in the frequency domain. We further extend this algorithm to three dimensions by adjusting the parameters of the 3-D kernel function during the iterative process. Numerical examples of synthetic and field data demonstrate the feasibility and adaptability of the proposed method. Its comparison with the dictionary learning method and optimal damped rank reduction methods shows that the proposed method can significantly improve the SNR of deep reflection signals and is a good tool for processing deep and ultra-deep seismic data.

Suppressing random noise and improving the signal-to-noise ratio of seismic data are of great significance for subsequent high-precision processing. As one of the most popular denoising methods, the algorithms based on deep learning usually utilize the network’s poor generalization ability for denoising processing. As a result, such methods usually face the problem of high training set construction costs and high computational costs. In addition, the widely used network UNet’s processing of multiscale features is limited to its U-shaped processing path. Therefore, we proposed a denoising method based on unsupervised learning, including a new denoising strategy based on spatial correlation, the corresponding training method, and a new network M-ResUNet. The training method allows networks to be trained directly on the test area and conducts independent training and processing for subtest areas with different geological characteristics. This not only effectively solves the problems caused by the method using generalization ability for denoising but also improves the denoising accuracy to a certain extent. In addition, M-ResUNet breaks through the limitation of the U-shaped processing path and effectively improves the training effect by controlling the bias degree of the network output to features with different scales. Compared with the traditional denoising methods and networks, the results’ test on synthetic and field data indicates that the proposed denoising method has superior performance in random noise attenuation, and M-ResUNet is effective.

Building subsurface models is a very important but challenging task in hydrocarbon exploration and development. The subsurface elastic properties are usually sourced from seismic data and well logs. Thus, we design a deep learning (DL) framework using Vision Transformer (ViT) as the backbone architecture to build the subsurface model using well log information as we apply full waveform inversion (FWI) on the seismic data. However, training a ViT network from scratch with limited well log data can be difficult to achieve good generalization. To overcome this, we implement an efficient self-supervised pre-training process using a masked autoencoder (MAE) architecture to learn important feature representations in seismic volumes. The seismic volumes required by the pre-training are randomly extracted from a seismic inversion, such as an FWI result. We can also incorporate reverse time migration (RTM) image into the seismic volumes to provide additional structure information. The pre-training task of MAE is to reconstruct the original image from the masked image with a masking ratio of 75%. This pre-training task enables the network to learn the high-level latent representations. After the pre-training process, we then fine-tune the ViT network to build the optimal mapping relationship between 2D seismic volumes and 1D well segments. Once the fine-tuning process is finished, we apply the trained ViT network to the whole seismic inversion domain to predict the subsurface model. At last, we use one synthetic data set and two field data sets to test the performance of the proposed method. The test results demonstrate that the proposed method effectively integrates seismic and well information to improve the resolution and accuracy of the velocity model.

In marine seismic exploration, the ocean bottom cable can receive both reflected P- and S-wave information from sub-seabed structures. Recently, a water-land dual-detector observation system has been developed to record acoustic pressure fields using water detectors and elastic displacement wave fields using land detectors. These water-land dual-detector data sets are commonly used to suppress multiples. However, the conventional elastic least-squares reverse time migration (LSRTM) based on single elastic wave equations cannot simultaneously utilize the ocean bottom dual-sensor (OBD) data. To solve the acoustic-elastic joint inverse problem by simultaneously using observed pressure and displacement records in the OBD data, we propose an OBD-data-based joint acoustic and decoupled-elastic LSRTM (JADE-LSRTM). This method constructs a new joint acoustic and decoupled-elastic misfit function, acoustic-elastic wave backward-adjoint operators and demigration operators in the curvilinear system, and gradient directions with respect to P- and S-wave velocity. In these acoustic-elastic coupled wavefield propagation operators, acoustic equations are used to calculate forward-propagated and backward-propagated pressure wavefields in the seawater based on water detector data sets, while decoupled elastic equations are applied to produce the displacement wavefields in the underlying elastic medium, with a higher computational efficiency than the traditionally combining individual acoustic and decoupled-elastic wavefield operators. Comparing to the conventional elastic LSRTM, the proposed method enables the use of acoustic data, which is crucial for OBD data. Two numerical examples demonstrate that the proposed curvilinear coordinated JADE-LSRTM based on water-land dual-detector data can produce accurate images in P- and S-component with higher efficiency and accuracy.

Compared with traditional adjoint migration, the least-squares migration (LSM) can effectively mitigate the unbalanced illumination and limited resolution associated with finite acquisition apertures, complex overburden structures and band-limited records. Data-domain LSM needs many times of Born modeling and adjoint migration to converge to a good solution, which is still challenging for large-scale 3D model under current computational capacity. To reduce computational cost and produce high-quality images, we directly approximate the Hessian inverse in the image-domain LSM using a new framework of generative adversarial networks (GANs). The migrated images, source illumination and migration velocity model are used as input data for the GANs, and the ground-truth reflectivity is utilized as the label data to train the network. Directly applying conventional GAN framework to implement the image-domain LSM leads to dislocated reflection events and incorrect images. To overcome this issue, we develop a new GAN framework that is more suitable for the Hessian approximation of image-domain LSM, which is named as LsmGANs. In the new framework, we use a max-pooling instead of convolution to downsample the feature maps to capture horizontal and vertical variations of reflectors. This enables us to map reflection events to correct location in downsampling. To address the lateral discontinuity of events in the predicted image from conventional GANs, we further apply multiple transform layers to strengthen feature transformation to guide Hessian approximation. Finally, we add the skip connection in the transform layer to enhance the information exchange of the feature channels and avoid the gradient vanishing problem to improve image resolution. Assembling predicted patches to construct a whole reflectivity image is a key step in the neural-network-based LSM. We investigate four strategies using different overlapping ratio and window functions to assemble the LSM patches and observe that less overlapping produces more patch-edge artifacts and the partition of unit with a Gaussian window has the best performance. Numerical experiments for synthetic and field data show that the proposed LsmGAN method can produce high-quality images with balanced amplitudes, reduced artifacts and improved resolution.

In viscoelastic media, reverse time migration (RTM) requires accurate Q value for energy compensation and phase correction. Full waveform inversion (FWI) has the potential to provide accurate Q value. In the wide frequency band seismic data, the quality factor Q almost does not change with frequency, which is called constant Q model, but its computational complexity is large. In the narrow frequency band data, when the memory variable order L is assumed to be 1 in the standard linear solid (SLS) model, the attenuation simulation effect of the subsurface Q value is acceptable, and the calculation amount is small. In this paper, we derive the back propagation wave field equation and gradient formula of full waveform inversion of viscoelastic medium based on the SLS model with memory variable order L = 1, and analyzes the effectiveness and feasibility of this method. In viscoelastic medium, low velocity body will produce low frequency noise in seismic record, which will lead to heterogeneous gradient energy. This low frequency noise is similar to the low frequency noise produced near the source of shallow strata. This will reduce the accuracy of the inversion results. The gradient preconditioning method using the pseudo-Hessian operator can make the gradient energy more balanced. Therefore, we introduce the pseudo-Hessian operator in gradient preconditioning and derives the gradient preconditioning formula for the full waveform inversion in viscoelastic media, which solves the problem of the non-uniform gradient energy. Examples show that the method can suppress the non-uniform gradient energy in the full waveform inversion of viscoelastic media.

Conventional full-waveform seismic inversion (FWI) tries to estimate a subsurface model that can accurately predict surface records by minimizing an${L} _{\mathbf {2}}$-norm misfit between observed and synthetic data. If the initial model is far away from the true model, the cycle-skipping issue might occur and the${L} _{\mathbf {2}}$-norm-based FWI produces a spurious model. To mitigate this problem, we present a novel FWI scheme using a high-dimensional local coherence misfit function. A 2-D/3-D window is first used to extract local seismic waveform from the common-shot gathers. Then, we apply a normalized cross-correlation to measure the coherence of local synthetic and observed records, which is used as the misfit to iteratively update the subsurface velocity model. The new misfit function enhances the contribution of phase fitting while reducing the amplitude contribution, which helps to increase the tolerance of FWI to an inaccurate initial velocity model. In addition, the computation of local waveform coherence along the temporal and spatial axes can adaptively balance the adjoint source amplitudes for strong near-offset reflections and weak far-offset refractions, which improves the low- wavenumber updates. Numerical experiments for synthetic and field data demonstrate that the proposed FWI scheme has a better tolerance to inaccurate starting models and is less sensitive to cycle-skipping issues compared with the conventional FWI method.

Anisotropic elastic wave vector decomposition for common-shot multicomponent data is a critical step in multicomponent seismic exploration. It is more difficult than the separation for elastic wavefield snapshots because the spatial derivatives along vertical direction cannot be calculated directly. In this letter, we propose an effective workflow to implement elastic wave vector decomposition for common-shot multicomponent data based on the pure-wave equation in transversely isotropic (TI) media. First, we extrapolate the multicomponent data backward in time to a reference surface using an anisotropic elastic wave equation, implement wave-mode separation by the projection of the extrapolated anisotropic elastic wavefields onto the polarization vectors of P- and S-wave-modes, and record the separated pure-wave data at the reference depth. Second, we extrapolate the separated anisotropic pure-wave data forward in time using the anisotropic pure-wave equations, perform the vector decomposition by the projection of the extrapolated pure-wave wavefields onto the polarization vectors of P- and S-modes, and record the vector anisotropic pure-wave data at the recording surface. This method can produce vector pure-mode data from the anisotropic medium with relatively high accuracy. Synthetic examples show the feasibility of this method.

In viscoelastic full waveform inversion, there is interference between the gradients of the P- and S- wave velocities, leading to the degradation of the inversion accuracy. Therefore, the interference problem in viscoelastic media needs to be solved beforehand. Even though the P-S separation method has been demonstrated to be effective in reducing the crosstalk in elastic media, little attention has been paid to its efficacy in viscoelastic media. We incorporate the P-S separation method in the full waveform inversion of viscoelastic media (QEFWI) and derive a gradient formula for the QEFWI. Numerical tests on the QEFWI demonstrate that the proposed method can successfully reduce the interference problem between the gradients of P- and S- wave velocities, resulting in significant improvement of the viscoelastic full waveform inversion.

Internal multiples bring great challenges for conventional seismic data imaging and interpretation, in which primaries only are regarded as signals. It is of great importance to suppress the internal multiples before conventional imaging, which is computationally expensive. Considering that internal multiples are also real reflections from the underground interfaces and reverse time migration (RTM) can realize the imaging of multiples. Instead of suppression before imaging, we propose to do the joint imaging of primaries and internal multiples based on RTM. However, the artifacts caused by the internal multiples should be suppressed to guarantee the precision. First, the imaging condition based on up-going and down-going extrapolated wavefield separation is introduced to improve the imaging accuracy. Then, the generation mechanism of image artifacts is analyzed, it concludes that nonphysical wave paths caused by the missing boundary wavefields under the surface result in the infidelity reversed-time backward-extrapolation, which is the main reason for the artifacts. Next, we propose to apply the modeling boundary wavefields with a match based on pseudomultichannel matching filter to compensate the missing boundary wavefields, which can avoid the image artifacts and realize high-precision joint imaging. Finally, the feasibility and effectiveness of the proposed method are verified by numerical examples and field data.

Prismatic waves carry steeply dipping structural information that primaries cannot contain. Therefore, prismatic waves are separately used in some migration methods to improve the illumination and imaging effect on steeply dipping structures. Least-squares reverse time migration of prismatic waves (LSRTM-P) can produce high-resolution images with improved steeply dipping structures. However, viscoelasticity exists widely on the Earth, which poses great difficulty for imaging. The effect of attenuation on prismatic waves is difficult to be compensated when conducting LSRTM-P because prismatic waves have three propagation paths. To overcome this problem, a$Q$-compensated LSRTM ($Q$-LSRTM)-P method is proposed by deriving$Q$-compensated forward-propagated operators and backward-propagated adjoint operators of prismatic waves, which compensates for$Q$attenuation along all the three propagation paths of prismatic waves. The proposed$Q$-LSRTM-P is conducted to update the image after applying the conventional$Q$-LSRTM. Besides, the proposed method can be adapted to the irregular surface media. Numerical examples on two synthetic and a field datasets verify that our method can produce better imaging results with clearer steeply dipping structures, higher signal-to-noise ratio (SNR), higher resolution, and more balanced amplitude than noncompensated LSRTM-P and conventional$Q$-LSRTM.

Curvilinear-grid-based least-squares reverse time migration (LSRTM) can produce an accurate image of complex subsurface structures. However, a huge amount of computational cost of LSRTM makes it difficult in real data applications, especially in 3-D cases. We propose a wavefield continuation operator in a new curvilinear-$\tau $domain to make the sampling space in the vertical direction to be uniform by stretching and compressing the low- and high-velocity zones, respectively. An objective function based on a conical wave encoding and student’s$T$distribution is constructed to improve the computational efficiency and robustness of LSRTM. The gradient formula is derived based on the objective function, and to correct the unbiased random estimation error, the idea of random optimization is introduced to obtain a weighted gradient. Demigration and adjoint wave equations in the curvilinear-$\tau $domain are derived to calculate the synthetic records and the backward-propagated wavefields. Numerical examples on synthetic and field datasets suggest that the proposed 3-D LSRTM method produces better images with a higher signal-to-noise ratio (SNR), more improved resolution, and more balanced amplitude than the conventional 3-D LSRTM and greatly improves the computational efficiency of 3-D LSRTM.

The acquisition footprint, finite-frequency source, and unbalanced subsurface illumination make it difficult for traditional adjoint-based migration to produce a high-quality image for complex structures. By fitting reflection events with linearized simulation data, least-squares migration (LSM) can iteratively incorporate the effects of the Gauss–Newton Hessian (GNH) to produce high-quality depth profiles. However, high computational costs of forward and adjoint simulations limit the LSM applications in production. In this study, we present an efficient approximation approach for the GNH and utilize it as a preconditioner for the misfit gradient of the LSM to accelerate its convergence. The analytic solution of the GNH in homogeneous media reveals that the columns of the GNH are local spatial functions. Based on this observation, we design a space-wavenumber filter to approximate the GNH for heterogeneous media, which can be efficiently computed with the S-transform and spectral division. The mixed-domain property of this filter allows it to automatically take the wavenumber dependence of the GNH into account and, therefore, helps to improve spatial resolution. Numerical examples demonstrate that the approximated GNH can considerably improve the image quality and speed up the convergence of the LSM.

The process of seismic multiple removal in oil seismic exploration is crucial for the imaging of underground structures with primary reflections. The inclusion of prediction and subtraction in the multiple removal method requires adaptive subtraction to remove the complex differences between the true and modeled multiples. The traditional adaptive subtraction method is generally expressed as a linear regression (LR) problem. In this article, we introduce U-net, a popular deep learning tool, to represent the complex differences between the true and modeled multiples in a nonlinear relationship. Thus, we present adaptive subtraction as a non-LR problem. The modeled multiples and full recorded seismic response with multiples and primaries are used as the input and labels to train U-net. The proposed U-net method is able to avoid over-fitting of the primaries due to the sufficient number of 2-D data windows for the training of U-net, as well as the network parameter regularization and the L1 norm minimization constraint on the primaries. The proposed U-net method attains 20.5 dB and 2.5 dB improvement in the signal-to-noise ratio (SNR) using the first synthetic data and second synthetic (Sigsbee2B) data set, respectively, compared with the traditional LR method, and a qualitative improvement for a real data set test.

In the deep-marine environment, seismic data contain multiples and are seriously affected by$Q$attenuation. Multiples have been used in migration to image the shadow zones and improve the resolution. However, the effect of attenuation on multiples is more serious than primaries because multiples have longer propagation paths. Therefore, we compensate the forward-propagated source-side and backward-propagated receiver-side wavefields along all the propagation paths of multiples. To fully use the primaries and multiples, we construct an objective function of least-squares reverse time migration of joint primaries and multiples (LSRTM-J) to update the imaging results by jointly using primaries and different-order multiples. In practice in a seawater medium, seismic waves can hardly be affected by attenuation. To decrease the computational cost, we divide the medium in the deep-marine environment into an acoustic medium part and a viscoacoustic medium part and derive the acoustic–viscoacoustic coupled compensated forward continuation operator, compensated adjoint operator, attenuated demigration operator, and gradient formula of joint primaries and multiples. To eliminate the severe scattering and diffracted noise caused by the strong-reflected irregular seabed interface, we mesh the velocity and$Q$models into curvilinear grids to perfectly match the seabed structure and realize the proposed viscoacoustic LSRTM-J in the curvilinear domain. Numerical examples on two typical models and a real data test suggest that the proposed method produces images with high SNR, high resolution, balanced amplitude, clear imaging structures, and strong deep region energy, and the total computational cost is the least of the other four conventional methods.