Abstract
Least-squares reverse time migration (LSRTM) is widely used in seismic imaging for high-resolution subsurface imaging, particularly in complex geological structures. This technique helps reveal detailed subsurface features that are crucial for fields such as oil and gas exploration and geotechnical studies. However, the iterative nature of LSRTM and its reliance on the least-squares approach result in high computational costs, making it challenging for large-scale applications. To address this challenge, this article proposes a safe Anderson-type-I LSRTM, built upon an enhanced 25-point finite difference scheme. This method incorporates a coefficient-optimized 25-point frequency-domain finite difference scheme, alongside Powell regularization, restart checking, and safety protection steps, which are applied to Anderson acceleration type I in order to improve stability and accelerate convergence. Model tests demonstrate that the proposed safe Anderson type-I LSRTM, based on the improved 25-point finite difference scheme, results in faster data residual convergence, higher imaging signal-to-noise ratio, superior resolution, clearer imaging of the in-phase axis, and a closer match between the imaging and the true reflection coefficient model, compared to the steepest descent method, conjugate gradient method, and limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) method. This method significantly enhances the practical feasibility of LSRTM for large-scale, high-resolution seismic imaging.
Introduction
In recent years, seismic imaging technology has played a crucial role in oil exploration, geological surveys, and earthquake monitoring. As a key method in seismic imaging, least-squares reverse time migration (LSRTM) has been widely used for subsurface structure interpretation and oil and gas exploration. LSRTM employs the least-squares approach, based on reverse time migration, to iteratively refine the calculated data so that it better matches the observed data. This method can generate high-precision images in complex media. Due to its strong amplitude-preserving capability, LSRTM can produce high-resolution images in complex subsurface imaging, attracting significant attention in both academia and industry, and has seen successful application in practical production (Gray, 1997). The LSRTM method can be categorized into three main types (Li and Qu, 2022; Qin et al., 2025; Wang et al., 2023; Li et al., 2024): least-squares ray-based migration, least-squares one-way wave migration, and least-squares two-way wave migration (Cole and Karrenbach, 1992). However, applying LSRTM to large-scale data, especially in three-dimensional cases, requires substantial computational effort due to the complexity of the structures.
Improving the computational efficiency of LSRTM can be approached in four ways: The first method is encoded in LSRTM, which aims to suppress the interference of crosstalk noise in imaging. Several LSRTM encoding methods have been proposed to achieve this, such as polarity encoding (Sun et al., 2022), frequency division encoding (Huang and Schuster, 2012), plane wave encoding (Dai and Schuster, 2013), amplitude encoding (Hu et al., 2016), and random encoding. The second method involves gradient preconditioning LSRTM: due to the complexity of inverting the Hessian matrix, direct solutions are not feasible. Typically, the Hessian matrix is approximated, and the inverse of this approximated Hessian matrix is used as a preconditioner to increase the convergence speed of the iterations. The Hessian matrix, under assumptions of infinite aperture and high-frequency approximations, can be treated as a diagonal matrix. Plessix and Mulder (2004) introduced and analyzed four diagonal Hessian matrices in the computation of frequency-domain migration. Valenciano et al. (2006) simplified the Hessian matrix to a diagonal matrix for gradient preconditioning. Aoki and Schuster (2009) used filters, estimated from two migration results, as preconditioners, that is, deblurring preconditioners, to accelerate the convergence of the inversion. Ren H et al. (2013) provided mathematical and physical insights into the Hessian matrix in seismic inversion imaging and discussed two simplified Hessian operators suitable for LSRTM and full waveform inversion. Gao et al. (2020) proposed using the Kronecker product summation approximation of the Hessian matrix to reduce the computational load of LSRTM. The third method utilizes regularization-constrained LSRTM: appropriate regularization constraints can improve the convergence speed of LSRTM. Two main types of regularization methods are commonly used. The first type applies constraints to common imaging point gathers, such as smoothness constraints (Kuehl and Sacchi, 2003) and dip constraints (Liu et al., 2013). The second type constrains the imaging results by denoising the signal, including Radon transform constraints (Dutta et al., 2015), Seislet domain sparse constraints (Dutta and Schuster, 2015), curvelet domain sparse constraints (Herrmann and Brown, 2009), prior information constraints, L1 sparse constraints (Wu et al., 2016), and structural constraints (Liu et al., 2012). The fourth method is improved gradient update methods for LSRTM: traditional LSRTM typically employs the steepest descent method or the conjugate gradient method for gradient updates. However, these conventional methods struggle to meet the increasing demands of large data volumes and high computational efficiency. To increase the computational efficiency of LSRTM, various new gradient update methods have been proposed, such as LSRTM based on LBFGS (Liu et al., 2016), LSRTM based on Anderson acceleration (Zhang et al., 2025; Huang et al., 2022), LSRTM based on the Adam optimizer (Vamaraju et al., 2021; Wang et al., 2022), and LSRTM based on stochastic quasi-Newton methods (Farias et al., 2023).
LSRTM is applicable not only in the time domain but also in the frequency domain (Shin et al., 2001; Kim et al., 2011; Zhang et al., 2025). Ren R et al. (2013) utilized Green’s functions of sources and receivers to effectively compute shot record data and gradients under the Born approximation, exploring the feasibility of implementing LSRTM in the frequency domain. They demonstrated that frequency-domain LSRTM is more efficient than time-domain LSRTM. Compared to time-domain LSRTM, frequency-domain LSRTM offers several notable advantages. First, for the implementation of frequency-domain LSRTM, Green’s functions in the frequency domain can be solved in parallel using LU factorization (Jo, 1996a). Second, in the frequency domain, the second-order time derivatives can be conveniently represented using angular frequency, which facilitates computation. Suitable Helmholtz equation solvers can also improve the efficiency of computing the LSRTM gradient in the frequency domain. Additionally, the velocity field values in frequency-domain LSRTM can be expressed as complex numbers, which accounts for the attenuation effects of seismic waves in the frequency domain (Operto et al., 2007). However, because the source observation data of actual records are typically collected as time series, implementing frequency-domain LSRTM requires transforming the time-domain observation data into frequency-domain data using the Fourier transform. This transformation is not conducive to the application of large-scale real data in frequency-domain LSRTM. Furthermore, the implementation of frequency-domain LSRTM typically involves matrix operations, which are memory-intensive and place significant demands on computer memory, thereby limiting the practical application of frequency-domain LSRTM.
Lysmer and Drake applied the finite element method to seismic wave forward modeling, marking the formal introduction of frequency-domain forward modeling (Lysmer et al., 1972). Pratt employed the traditional second-order difference scheme for seismic wave numerical simulation and demonstrated through model tests that at least 13 grid points per wavelength are required to ensure accuracy (Pratt and Worthington, 1990). In the same year, Luo and Schuster (1990) proposed the staggered grid method; however, its large memory consumption and severe dispersion phenomena hindered the widespread application of frequency-domain forward modeling. It was not until Jo introduced a new coordinate system by rotating the coordinate axes and weighted the wavefield values obtained through finite difference discretization under two different coordinate axes, thus optimizing the nine-point difference scheme, that the grid point requirement was reduced to at least four points per wavelength, mitigating dispersion and reducing memory usage (Jo et al., 1996b). Subsequently, Hustedt solved the one-dimensional velocity-stress equation using the rotated coordinate system (Hustedt et al., 2004). Following this, various researchers rotated the coordinate axes at different angles, extending the method to 17-point, 25-point, and other different schemes (Jo, 1996a; Cao and Chen, 2012; Liu et al., 2013; Shin and Sohn, 1998).
In this article, we introduce Powell regularization, restart checks, and safety protection steps into the Anderson acceleration type-I method to address the instability and oscillation of data residuals in the Anderson acceleration type-I method, thereby improving its computational efficiency. The enhanced Anderson acceleration type-I method is then applied to frequency-domain LSRTM with the improved 25-point difference scheme, effectively accelerating the convergence rate of data residuals in LSRTM.
Theory
Gradient update methods for LSRTM
The least-squares error objective function can be expressed by the acoustic wave equations as shown in Equations 1, 2:where E is the misfit function, L is a linear modeling operator, m is a parameter field, and d is the observed data. The imaging result of LSRTM can be obtained using the equation above, with the least-squares solution beingwhere T denotes the transpose of a matrix. The Hessian matrix operator, , can also be referred to as the blurring integral kernel operator. The Hessian matrix is generally defined as the second-order partial derivative of the objective function with respect to the parameter vector. Because the inverse of the Hessian operator is sensitive to noise and data errors, and computing the inverse of the Hessian operator is computationally intensive, it is usually not directly computed. Instead, indirect iterative methods, such as the steepest descent method and the conjugate gradient method, are employed to approximate the inverse of the Hessian operator (Broyden et al., 1973). The Anderson acceleration algorithm discussed in this article is also an indirect iterative method for approximating the inverse of the Hessian operator. The Anderson acceleration algorithm is used to speed up fixed-point iterations. In LSRTM, it is generally used to accelerate the fixed-step steepest descent method. It is not used to accelerate the conjugate gradient method because the conjugate gradient method updates the gradient and step size based on the previous iteration, starting from the second iteration, which does not meet the fixed-step requirement of fixed-point iterations. Therefore, Anderson acceleration is mainly used to accelerate the fixed-step steepest descent method. The LSRTM using the steepest descent method is composed of the following equations:where is the gradient at the th iteration and is the step size at the kth iteration. The LSRTM using the steepest descent method first obtains the gradient by backpropagating the residuals through reverse time migration (RTM) using Equation 3. Then, the step size is calculated using Equation 4, and finally, the reflectivity is updated using Equation 5. In the fixed-step steepest descent method used in this article, the step size is calculated during the first iteration and fixed for subsequent iterations. Instead of recalculating the step size for each iteration, Anderson acceleration is employed to speed up the update of the reflectivity by utilizing past iteration results. This approach significantly reduces the computation time by eliminating the need to calculate the step size.
Coefficient-optimized 25-point difference scheme
To address the substantial computational cost associated with solving large-scale linear algebraic equation systems, Sun et al. (2022) proposed an improved 25-point difference scheme for elastic waves. Building upon this method, we introduce it into the acoustic wave equation. In the current optimization scheme, we employ different finite difference (FD) formulas to approximate spatial derivatives while managing the mass acceleration term by distributing the mass acceleration at the central point across all grid points in the FD scheme.
The frequency-domain acoustic wave equation is as follows:where P denotes an acoustic pressure field, v denotes the velocity, w denotes the frequency, S denotes the source term, and x and z are the spatial coordinates. We apply the weighted average method to discretize the partial derivatives and :
Therefore,where represent the weighting coefficients, and they are , . By substituting Equations 7–9 into Equation 6, we obtain the improved 25-point difference scheme for the acoustic wave:
The values of the coefficients in the equations are as shown in Equation 10
Next, we solve for the weighting coefficients of the improved 25- point difference scheme. The dispersion relation is given by Equation 11
Similarly, the weighting coefficients can be obtained by solving the equations using the conjugate gradient method.
LSRTM based on safe Anderson acceleration type I
Applying Anderson acceleration to LSRTM can speed up inversion convergence (Zhang et al., 2025); however, the classical scheme may become unstable. To eliminate the instabilities of the type-I Anderson acceleration and simultaneously reduce its computational cost, we incorporate Powell regularization (Supplementary Appendix A), a restart check (Supplementary Appendix B), and a safeguarding step (Supplementary Appendix C) into the type-I Anderson framework, yielding a safe first-type Anderson acceleration. This safe variant is then embedded in LSRTM to accelerate the convergence of data residuals and to improve the overall efficiency of LSRTM. The detailed implementation is summarized in Algorithm 1.
Algorithm 1
Input: Choose an initial reflectivity model , memory parameter , regularization constant , safety protection constants and , and fixed-point iteration operator . In the first iteration, calculate the data residual , obtain the gradient , update the reflectivity , and output as the reflectivity for the first iteration. Calculate , initialize the regularization function , set accumulation constants , , and .
fork = 2, … , Maximum number of iterations or number of iterations to reach suitable
convergence do
Step one: ifdo
update
else do
update
end if
Calculate the data residual , obtain the gradient ,
update the reflectivity ; calculate
, ; increment by one.
Step two: ifdo
Save and , output as the reflectivity for
this iteration; increment by one.
else do
Set ; save ; recalculate the data residual
, recalculate the gradient , and update the reflectivity
again; save and output as the
reflectivity for this iteration.
end if
Save ; Calculate
Step three: if k = 2 do
Set
else do
Set
end if
Step four: ifordo
reset , , ,
end if
Step five: update , calculate ,
, update
end for
Numerical examples
To test the applicability of safe Anderson acceleration type-I LSRTM in the frequency domain, experiments were conducted using the Canadian foothills model (Figure 1). The parameters of the Canadian foothills model, including source frequency, shot point and receiver distribution, memory variables for safe Anderson acceleration type I, and limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) (Bautista et al., 2020), were kept consistent. The sampling frequency interval was set to 0.5 Hz, with 1,000 sampling points, and the maximum sampling frequency was 500 Hz. Using the Canadian foothills model, Born forward modeling was performed. Figures 2a–c show wavefield snapshots at 45 Hz, 60 Hz, and 75 Hz with a shot point at (2010 m, 10 m). These figures illustrate the propagation process of the acoustic wavefield at different frequency slices. Figure 3 presents the shot records of the Born forward modeled data for 20 shots. Figure 4 compares the data residual convergence curves of frequency-domain LSRTM using the steepest descent method, conjugate gradient method, LBFGS (Bautista et al., 2020), and safe Anderson acceleration.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
The figure shows that as CPU time increases, the data residuals for all methods gradually decrease, indicating that the algorithms are continuously optimizing and approaching the optimal solution. However, there are significant differences in the convergence rates among the different methods. First, the safe AA-type-I method, represented by the purple dotted line, shows the fastest convergence speed within 8,500 s of CPU time. This means that the LSRTM method with safe Anderson acceleration can achieve lower data residuals more quickly, demonstrating its significant advantage in accelerating data residual convergence. In contrast, the steepest descent (SD) method, conjugate gradient (CG) method, and LBFGS method have relatively slower convergence rates. The SD method, represented by the solid cyan line, has the slowest convergence rate throughout the process, while the CG method (dashed blue line) and the LBFGS method (dashed red line) are faster than the SD method but still not as fast as the safe AA-type-I method. This difference indicates that the safe Anderson acceleration-based LSRTM can more effectively accelerate data residual convergence in both time and frequency domains than traditional methods. This not only means that the method has an advantage in computational efficiency but also suggests that it may have higher accuracy and stability when dealing with complex inversion problems.
Figure 5 displays the frequency-domain LSRTM images generated using four different optimization methods: SD, CG, LBFGS, and safe Anderson acceleration after 8,500 s of CPU time, along with their similarity to the Canadian foothills model reflectivity model shown in Figure 1c. The similarity between the generated images and the model reflectivity was calculated using Madagascar software, a widely recognized tool for geophysical data processing.
FIGURE 5
The detailed comparison shown in Figure 6 provides a closer look at the imaging quality across different methods. Figures 5, 6 show that under the same computational time constraints, the imaging resolution of frequency-domain LSRTM based on safe Anderson acceleration is higher than that of the other methods. This is particularly evident within the area enclosed by the dashed red circle, where the imaging of the stratum layers is significantly better when using safe Anderson acceleration. In the images produced by the steepest descent method, conjugate gradient method, and LBFGS, the stratum within the dashed red circle is poorly defined and lacks clarity. These methods struggle to accurately capture the fine details of the subsurface structure, leading to a less accurate representation of the geological features. In contrast, the frequency-domain LSRTM based on safe Anderson acceleration effectively restores the stratum layers, providing a clearer and more detailed image of the subsurface. This enhanced imaging capability of safe Anderson acceleration-based LSRTM is crucial for accurate geological interpretation and can significantly improve the reliability of subsurface exploration and resource assessment. The ability to resolve fine details within a given time frame is a testament to the computational efficiency and effectiveness of the safe Anderson acceleration method in the context of frequency-domain LSRTM. Furthermore, the similarity comparison shows that the frequency-domain LSRTM image based on safe Anderson acceleration has a higher similarity to the Canadian foothills model reflectivity model than the other methods. This confirms that the imaging quality of LSRTM based on safe Anderson acceleration is superior to the other methods in both the time and frequency domains under the same time constraints.
FIGURE 6
The Marmousi model was adopted to evaluate the noise sensitivity of the proposed method. The model is 2.0 km deep and 3.68 km wide, discretized with a uniform grid spacing of 10 m in both directions. Figures 7a,b show the P-wave velocity and the corresponding smoothed P-wave velocity models, respectively. A Ricker wavelet with a dominant frequency of 20 Hz was used as the source-time function. The acquisition geometry consisted of 19 shots evenly distributed from (10 m, 10 m) to (3610 m, 10 m) with a 200 m interval, and 368 receivers deployed on the surface from (10 m, 10 m) to (3680 m, 10 m) with a 10 m spacing. The recording parameters were a 0.5 m sampling interval and a maximum recording time of 2.5 s. Figures 8a,b show the observed data of the 10th shot using our coefficient-optimized 25-point scheme and the Sun et al. (2022) 25-point scheme, respectively. As can be seen from the figures, our coefficient-optimized 25-point scheme has no dispersion, with clear and continuous reflection axes and high accuracy. In contrast, the Sun et al. (2022) 25-point scheme still has obvious dispersion. Therefore, our coefficient-optimized 25-point scheme is suitable for complex simulations and has high accuracy. To examine the performance of LSRTM under noisy conditions, random noise corresponding to a signal-to-noise ratio (SNR) of 0.14 was added to the synthetic data. Figure 8c displays the noisy shot gather of the 10th shot, where the effective reflections are almost completely buried in noise. Figure 9 presents the images at iteration 8,000 obtained with SD, CG, LBFGS (Bautista et al., 2020), and the proposed AA-type-I LSRTM schemes, together with their structural similarity (SSIM) indices against the true Marmousi reflectivity model. The proposed method produces images with higher SNR and better resolution even under strong noise contamination, and its SSIM is consistently higher than those of the competing algorithms. The proposed method produces images with higher SNR and better resolution even under strong noise contamination, and its SSIM is consistently higher than those of the competing algorithms. Figure 9 shows that the images obtained using the AA-type-I LSRTM scheme are of superior quality than those from the SD, CG, and LBFGS methods. Figures 9a–d reveal a significant difference in image quality among the different methods. The AA-type-I LSRTM scheme (Figure 9d) clearly outperforms the SD (Figure 9a), CG (Figure 9b), and LBFGS (Figure 9c) schemes in terms of imaging quality. Particularly in the presence of high noise levels, the AA-type-I LSRTM scheme is effective in restoring subsurface structures, whereas the other methods lose a significant amount of detail during the imaging process. Figures 9e–h further show the differences in the imaging results among the various methods. The AA-type-I LSRTM scheme (Figure 9h) shows clearer stratigraphic structures and higher contrast in the imaging results, while the other methods exhibit varying degrees of blurriness and noise interference. This indicates that the AA-type-I LSRTM scheme has stronger robustness and higher imaging quality when dealing with noisy data. Additionally, the SSIM indices in Figure 9 confirm the imaging quality advantage of the AA-type-I LSRTM scheme. The SSIM index is a metric that measures the structural similarity between two images, with values ranging from 0 to 1, where higher values indicate greater structural similarity between images. The figure shows that the SSIM index of the AA-type-I LSRTM scheme is significantly higher than those of the other methods, further demonstrating its advantage in imaging quality.
FIGURE 7
FIGURE 8
FIGURE 9
Real data
To evaluate the proposed scheme on field records, we applied LSRTM to a land survey acquired over an 18 km × 4 km area. The migration velocity shown in Figure 10 was obtained through tomographic inversion. Figure 11 displays the first five shot gathers collected by 240 receivers spaced 40 m apart along the surface. The acquisition geometry includes 264 shots, and each trace is recorded for 4.0 s with a 2 m sampling interval. Prior to migration, the gathers are conditioned by attenuating random noise, eliminating direct arrivals, and discarding dead traces. Figure 12a presents the results of reverse time migration (RTM) imaging, while Figures 12b,c show the imaging results of our method and the LBFGS method at the same calculation time. Figure 12a presents the results of reverse time migration (RTM) imaging, while Figures 12b,c show the imaging results of our method and the LBFGS method at the same calculation time, respectively. Figure 12 illustrates that our method delivers a crisper subsurface image, with better resolved reflectors and more continuous structural dips, whereas the LBFGS yields a slightly smeared result. Our method still faithfully recovers the main geological features. Specifically, the RTM imaging results shown in Figure 12a depict the basic shape of the subsurface structure but may lack detail in certain areas. In contrast, the image generated by our method in Figure 12b displays clearer stratigraphic interfaces and more intricate geological structures along both the depth and distance axes. This indicates that our method is more effective at preserving geological information when dealing with complex subsurface structures, thus providing more accurate imaging results. The imaging results from the LBFGS method shown in Figure 12c also demonstrate the subsurface structure, but compared to our method, the reflectors and structural dips in the image are less distinct and exhibit some degree of smearing. This blurriness may be due to the limitations of the LBFGS method in handling noise or complex geological structures.
FIGURE 10
FIGURE 11
FIGURE 12
Conclusion
The proposed safe Anderson acceleration type-I least-squares reverse time migration (LSRTM), based on an improved 25-point finite difference scheme, effectively increases both computational efficiency and image quality. By incorporating an enhanced 25-point frequency-domain finite difference scheme, along with Powell regularization, restart checks, and safety protection steps in Anderson acceleration type I, the method achieves significant improvements in stability and convergence speed. Model tests demonstrate that this approach results in faster data residual convergence, a higher imaging signal-to-noise ratio and resolution, clearer in-phase axis imaging, and greater similarity between the imaging results and the true reflection coefficient model, outperforming traditional methods such as the steepest descent method, conjugate gradient method, and the LBFGS method. The proposed method mitigates, to a certain extent, the explosive weights, oscillations, or even overflow caused by the exploding condition number of the least-squares coefficient matrix and achieves a markedly faster convergence. Nevertheless, the memory footprint remains heavy: every step requires solving an m × m linear least-squares problem, so both memory usage and floating-point operations grow linearly with m, turning large-scale problems into a bottleneck.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
HY: Formal Analysis, Methodology, Resources, Writing – original draft. GL: Formal Analysis, Investigation, Writing – original draft. SC: Formal Analysis, Validation, Writing – review and editing. SJ: Formal Analysis, Investigation, Writing – review and editing. GZ: Formal Analysis, Investigation, Validation, Writing – review and editing. QZ: Methodology, Supervision, Writing – review and editing. JS: Investigation, Writing – review and editing. YQ: Methodology, Supervision, Writing – review and editing.
Funding
The author(s) declared that financial support was received for this work and/or its publication. The program of China National Petroleum Corporation (2023ZZ14YJ05), Youth Science and Technology Special Programme of CNPC (no. 2024DQ03013), and Opening Foundation of State Key Laboratory of Deep Oil and Gas (no. SKLDOG2024-KFYB-08).
Conflict of interest
Authors HY, GL, SC, SJ, GZ, QZ, and JS were employed by PetroChina.
The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2025.1665343/full#supplementary-material
References
1
AokiN.SchusterG. T. (2009). Fast least squares migration with a deblurring filer. Geophysics74 (6), WCA83–WCA93. 10.1190/1.3155162
2
BautistaK. D.PetersF. C.de SouzaR. V.MansurW. J. (2020). L-BFGS based least-squares reverse time migration for reflectivity imaging. XL Ibero-Latin Am. Congr. Comput. Methods Eng.2 (02).
3
BroydenC. G.DennisJr J. E.MoréJ. J. (1973). On the local and superlinear convergence of quasi-newton methods. IMA J. Appl. Math.12 (3), 223–245. 10.1093/imamat/12.3.223
4
CaoS. H.ChenJ. B. (2012). A 17-point scheme and its numerical implementation for high-accuracy modeling of frequency-domain acoustic equation. Chin. J. Geophys. (In Chinese)55 (10), 3440–3449. 10.6038/j.issn.0001-5733.2012.10.027
5
ColeS.KarrenbachM. (1992). Least-squares kirchhoff migration, 101–110.
6
DaiW.SchusterG. T. (2013). Plane-wave least-squares reverse-time migration. Geophysics78 (4), S165–S177. 10.1190/geo2012-0377.1
7
DuttaG.SchusterG. T. (2015). “Sparse least-squares reverse time migration using seislets,” in SEG Technical Program Expanded Abstracts 2015. (Tulsa, OK: Society of Exploration Geophysicists), 4232–4237.
8
DuttaG.GiboliM.WilliamsonP. (2015). “Least-squares reverse time migration with factorization-free preconditioning,” in SEG Technical Program Expanded Abstracts 2015. (Tulsa, OK: Society of Exploration Geophysicists), 4270–4275.
9
FariasF.SouzaA. A.PestanaR. C. (2023). Applying a stochastic quasi-newton optimizer to least-squares reverse time migration. Computers and Geosciences171, 105292. 10.1016/j.cageo.2022.105292
10
GaoW.MatharuG.SacchiM. D. (2020). Fast least-squares reverse time migration via a superposition of kronecker products. Geophysics85 (2), S115–S134. 10.1190/geo2019-0254.1
11
GrayS. H. (1997). True-amplitude seismic migration: a comparison of three approaches. Geophysics62 (3), 929–936. 10.1190/1.1444200
12
HerrmannF. J.BrownC. R.ErlanggaY. A.MoghaddamP. P. (2009). Curvelet-based migration preconditioning and scaling. Geophysics74 (4), A41–A46. 10.1190/1.3124753
13
HuJ.WangH.FangZ.ZhangJ. (2016). Efficient amplitude encoding least-squares reverse time migration using cosine basis. Geophysical Prospecting64 (6), 1483–1497. 10.1111/1365-2478.12356
14
HuangY.SchusterG. T. (2012). Multisource least-squares migration of marine streamer and land data with frequency-division encoding. Geophysical Prospecting60 (4), 663–680. 10.1111/j.1365-2478.2012.01086.x
15
HuangC.QuY.LiZ. (2022). Least-squares reverse time migration using the safe Type-Ⅰ Anderson Accaleration. 83rd EAGE Annual Conference and Exhibition2022 (1), 1–5. 10.3997/2214-4609.202210435
16
HustedeB.OpertoS.VirieuxJ. (2004). Mixed grid and ataggered grid finite-difference methods for frequency-domain acoustic wave modelling. Geophys. J. Int.157 (3), 1269–1296. 10.1111/j.1365-246X.2004.02282.x
17
JoC. H. (1996a). An optimal 9-point finite-differnce frequency-space 2-D scalar wave extrapolator. Geophysics61 (2), 529–537. 10.1190/1.1443979
18
JoC. H.ShinC.SuhJ. H. (1996b). An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator. Geophysics61 (2), 529–537. 10.1190/1.1443979
19
KimY.MinD. J.ShinC. (2011). Frequency-domain reverse-time migration with source estimation. Geophysics76 (2), S41–S49. 10.1190/1.3534831
20
KuehlH.SacchiM. D. (2003). Least-squares wave-equation migration for AVP/AVA inversion. Geophysics68 (1), 262–273. 10.1190/1.1543212
21
LiZ.QuY. (2022). Research progress on seismic imaging technology. Petroleum Science19 (1), 128–146. 10.1016/j.petsci.2022.01.015
22
LiJ.QuY.HuangC.LiY.LiZ. (2024). Ocean bottom dual-sensor Q-compensated elastic least-squares reverse time migration based on viscoacoustic and separated-viscoelastic coupled equations. Geophysics89 (3), S155–S173. 10.1190/geo2023-0433.1
23
LiuG. C.ChenX. H.SongJ. Y.RuiZ. H. (2012). A stabilized least-squares imaging condition with a structure constraint. Applied Geophysics9 (4), 459–467. 10.1007/s11770-012-0358-9
24
LiuL.LiuH.LiuH. W. (2013). Optimal 15-point finite difference forward modeling in frequency-space domain. Chinese Journal of Geophysics (In Chinese)56 (2), 644–652. 10.6038/cjg20130228
25
LiuY. J.LiZ. C.WuD. (2013). The research on local slope constrained least-squares migration. Chinese Journal of Geophysics (In Chinese)56 (3), 1003–1011. 10.6038/cjg20130328
26
LiuY.TengJ.XuT.BaiZ.LanH.BadalJ. (2016). An efficient step-length formula for correlative least-squares reverse time migration. Geophysics81 (4), S221–S238. 10.1190/geo2015-0529.1
27
LuoY.schusterG. T. (1990). “Wave-equation traveltime waveform inversion,” in SEG Technical Program Expanded Abstracts 1990. (Tulsa, OK: Society of Exploration Geophysicists), 1223–1228.
28
LysmerJ.DrakeL. A.BoltB. (1972). A finite element method for seismology. Methods in Computational Physics11, 181–216. 10.1016/b978-0-12-460811-5.50009-x
29
OpertoS.VirieuxJ.AmestoyP.L’ExcellentJ. Y.GiraudL.AliH. B. H. (2007). 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: a feasibility study. Geophysics72 (5), SM195–SM211. 10.1190/1.2759835
30
PlessixR. E.MulderW. A. (2004). Frequency-domain finite-difference amplitude-preserving migration. Geophysical Journal International157, 975–987. 10.1111/j.1365-246x.2004.02282.x
31
PrattR. G.WorthingtonM. (1990). Inverse theory applied to multi-source cross-hole tomography. Part 1:Acoustic wave-equation method. Geophysical Prospecting38 (3), 287–310. 10.1111/j.1365-2478.1990.tb01846.x
32
QinC.GaoF.WuK. (2025). An overview of research progress on wave-equation depth migration for seismic imaging. CT Theory and Applications, 1–15. 10.15953/j.ctta.2025.251 (in Chinese).
33
Ren HH.WangH.ChenS. (2013). Least-squares reverse time migration in frequency domain using the adjoint-state method. Journal of Geophysics and Engineering10 (3), 035002. 10.1088/1742-2132/10/3/035002
34
Ren RR.HuangG. H.WangH. Z. (2013). A research on the Hessian operator in seismic inversion imaging. Chinese Journal of Geophysics56 (7), 2429–2436. 10.6038/cjg20130725
35
ShinC. S.SohnH. J. (1998). A frequency-space 2-D scalar wave extrapolator using extend-ed 25-point finite-difference operater. Geophysics63 (1), 289–296. 10.1190/1.1444323
36
ShinC.JangS.MinD. J. (2001). Improved amplitude preservation for prestack depth migration by inverse scattering theory. Geophysical Prospecting49 (5), 592–606. 10.1046/j.1365-2478.2001.00279.x
37
SunW.LiZ.QuY. (2022). The 3D conical Radon transform for seismic signal processing. Geophysics87 (5), V481–V504. 10.1190/geo2021-0278.1
38
ValencianoA. A.BiondiB.GuittonA. (2006). Target-oriented wave-equation inversion. Geophysics71 (4), A35–A38. 10.1190/1.2213359
39
VamarajuJ.VilaJ.Araya-PoloM.DattaD.SidahmedM.SenM. K. (2021). Minibatch least-squares reverse time migration in a deep-learning framework. Geophysics86 (2), S125–S142. 10.1190/geo2019-0707.1
40
WangY.HuangC.QuY.LiM.LiJ. (2023). Velocity-adaptive irregular point spread function deconvolution imaging using X-shaped denoising diffusion filtering. IEEE Trans. Geosci. Remote. Sens.61, 5916808. 10.1109/TGRS.2023.3303184
41
WangS. W.SongP.TanJ. (2022). The least-squares reverse time migration with gradient optimization based on QHAdam. Chinese Journal of Geophysics (In Chinese)65 (7), 2673–2680. 10.6038/cjg2022P0364
42
WuD.YaoG.CaoJ.WangY. (2016). Least-squares RTM with L1 norm regularisation. Journal of Geophysics and Engineering13 (5), 666–673. 10.1088/1742-2132/13/5/666
43
ZhangQ.LiZ.QuY.YangJ.ZhangL. (2025). Prestack depth migration method based on direct expansion of the second kind of Chebyshev polynomials. Oil Geophysical Prospecting60 (3), 710–717. 10.13810/j.cnki.issn.1000-7210.20240328
Summary
Keywords
improved 25-point frequency differential format, least-squares reverse time migration, Powell regularization, restart checking, safety protection, safe Anderson type I
Citation
Yang H, Lei G, Cheng S, Jian S, Zhao G, Zeng Q, Sui J and Qu Y (2026) Safe Anderson type-I least-squares reverse time migration based on a coefficient-optimized 25-point difference scheme. Front. Earth Sci. 13:1665343. doi: 10.3389/feart.2025.1665343
Received
14 July 2025
Revised
19 December 2025
Accepted
29 December 2025
Published
24 February 2026
Volume
13 - 2025
Edited by
Sanjit Kumar Pal, IIT(ISM) Dhanbad, India
Reviewed by
Jincheng Xu, Southern University of Science and Technology, China
Jianguang Han, Chinese Academy of Geological Sciences (CAGS), China
Updates
Copyright
© 2026 Yang, Lei, Cheng, Jian, Zhao, Zeng, Sui and Qu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qingcai Zeng, zqc69@petrochina.com.cn; Yingming Qu, quyingming@upc.edu.cn
Disclaimer
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.