The least-squares error objective function can be expressed by the acoustic wave equations as shown in Equations 1, 2: E m = 1 2 ∥ L m − d ∥ 2 , (1) 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 being m = L T L − 1 L T d , (2) where T denotes the transpose of a matrix. The Hessian matrix operator, L T L , 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 L T L − 1 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: g k = L T L m k − d , (3) μ k = g k T g k L g k T L g k , (4) m k + 1 = m k − μ k g k , (5) where g k is the gradient at the k th iteration and μ k 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.

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: ∇ 2 P x , z , ω + w 2 v 2 P x , z , ω = − S ω , (6) 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 ∂ 2 P ∂ x 2 and ∂ 2 P ∂ z 2 : ∂ 2 P ∂ x 2 ≈ α 0 Δ x 2   a P m + 1 , n − 2 P m , n + P m − 1 , n + 1 − a 2 P m + 1 , n + 1 − 2 P m , n + 1 + P m − 1 , n + 1 + 1 − a 2 P m + 1 , n − 1 − 2 P m , n − 1 + P m − 1 , n − 1 + α 1 Δ x 2   a P m + 2 , n − 2 P m + 1 , n + P m , n + 1 − a 2 P m + 2 , n + 1 − 2 P m + 1 , n + 1 + P m , n + 1 + 1 − a 2 P m + 2 , n − 1 − 2 P m + 1 , n − 1 + P m , n − 1 + α 1 Δ x 2   a P m , n − 2 P m − 1 , n + P m − 2 , n + 1 − a 2 P m , n + 1 − 2 P m − 1 , n + 1 + P m − 2 , n + 1 + 1 − a 2 P m , n − 1 − 2 P m − 1 , n − 1 + P m − 2 , n − 1 + α 2 Δ x 2   a P m + 1 , n + 1 − 2 P m , n + 1 + P m − 1 , n + 1 + 1 − a 2 P m + 1 , n + 2 − 2 P m , n + 2 + P m − 1 , n + 2 + 1 − a 2 P m + 1 , n − 2 P m , n + P m − 1 , n + α 2 Δ x 2   a P m + 1 , n − 1 − 2 P m , n − 1 + P m − 1 , n − 1 + 1 − a 2 P m + 1 , n − 2 P m , n + P m − 1 , n + 1 − a 2 P m + 1 , n − 2 − 2 P m , n − 2 + P m − 1 , n − 2 + α 3 Δ x 2   a P m + 2 , n + 1 − 2 P m + 1 , n + 1 + P m , n + 1 + 1 − a 2 P m + 2 , n + 2 − 2 P m + 1 , n + 2 + P m , n + 2 + 1 − a 2 P m + 2 , n − 2 P m + 1 , n + P m , n + α 3 Δ x 2   P m + 2 , n − 1 − 2 P m + 1 , n − 1 + P m , n − 1 + 1 − a 2 P m + 2 , n − 2 P m + 1 , n + P m , n a + 1 − a 2 P m + 2 , n − 2 − 2 P m + 1 , n − 2 + P m , n − 2 + α 3 Δ x 2   a P m , n + 1 − 2 P m − 1 , n + 1 + P m − 2 , n + 1 + 1 − a 2 P m , n + 2 − 2 P m − 1 , n + 2 + P m − 2 , n + 2 + 1 − a 2 P m , n − 2 P m − 1 , n + P m − 2 , n + α 3 Δ x 2   a P m , n − 1 − 2 P m − 1 , n − 1 + P m − 2 , n − 1 + 1 − a 2 P m , n − 2 P m − 1 , n + P m − 2 , n + 1 − a 2 P m , n − 2 − 2 P m − 1 , n − 2 + P m − 2 , n − 2 . (7)

Therefore, ω 2 v 2 P ≈ ω 2 v 2 d 0 P m , n + d 1 P m + 1 , n + P m − 1 , n + d 2 P m , n + 1 + P m , n − 1 + d 3 P m + 1 , n + 1 + P m − 1 , n + 1 + P m + 1 , n − 1 + P m − 1 , n − 1 + d 4 P m + 2 , n + P m − 2 , n + d 5 P m , n + 2 + P m , n − 2 + d 6 P m + 2 , n + 1 + P m − 2 , n + 1 + P m + 2 , n − 1 + P m − 2 , n − 1 + d 7 P m + 1 , n + 2 + P m − 1 , n + 2 + P m + 1 , n − 2 + P m − 1 , n − 2 + d 8 P m + 2 , n + 2 + P m − 2 , n + 2 + P m + 2 , n − 2 + P m − 2 , n − 2 , (8) where α i , a , b , d i represent the weighting coefficients, and they are d 0 + 2 d 1 + 2 d 2 + 4 d 3 + 2 d 4 + 2 d 5 + 4 d 6 + 4 d 7 + 4 d 8 = 1 , α 0 + 2 α 1 + 2 α 2 + 4 α 3 = 1 . By substituting Equations 7–9 into Equation 6, we obtain the improved 25-point difference scheme for the acoustic wave: A u m − 2 , n − 2 + B u m − 1 , n − 2 + C u m , n − 2 + B u m + 1 , n − 2 + A u m + 2 , n − 2 + D u m − 2 , n − 1 + E u m − 1 , n − 1 + F u m , n − 1 + E u m + 1 , n − 1 + D u m + 2 , n − 1 + G u m − 2 , n + H u m − 1 , n + I u m , n + H u m + 1 , n + G u m + 2 , n + D u m − 2 , n + 1 + E u m − 1 , n + 1 + F u m , n + 1 + E u m + 1 , n + 1 + D u m + 2 , n + 1 + A u m − 2 , n + 2 + B u m − 1 , n + 2 + C u m , n + 2 + B u m + 1 , n + 2 + A u m + 2 , n + 2 = 0 . (9)

The values of the coefficients in the equations are as shown in Equation 10 A = 1 − a 2 α 3 Δ x 2 + 1 − b 2 α 3 Δ z 2 + ω 2 v 2 d 8 , B = 1 − a 2 α 2 Δ x 2 + a − 1 α 3 Δ x 2 + 1 − b 2 α 2 Δ z 2 + b α 3 Δ z 2 + ω 2 v 2 d 7 , C = a − 1 α 2 Δ x 2 + 1 − a α 3 Δ x 2 + b α 2 Δ z 2 + 1 − b α 3 Δ z 2 + ω 2 v 2 d 5 , D = 1 − a 2 α 1 Δ x 2 + a α 3 Δ x 2 + 1 − b 2 α 1 Δ x 2 + b − 1 α 3 Δ x 2 + ω 2 v 2 d 6 , E = 1 − a 2 α 0 Δ x 2 + a − 1 α 1 Δ x 2 + a α 2 Δ x 2 − 2 a α 3 Δ x 2 + 1 − b 2 α 0 Δ z 2 + b α 1 Δ z 2 + b − 1 α 2 Δ z 2 − 2 b α 3 Δ z 2 + ω 2 v 2 d 3 , F = a − 1 α 0 Δ x 2 + 1 − a α 1 Δ x 2 − 2 a α 2 Δ x 2 + 2 a α 3 Δ x 2 + b α 0 Δ z 2 + 1 − b α 1 Δ z 2 − 2 b α 2 Δ z 2 + 2 b − 1 α 3 Δ z 2 + ω 2 v 2 d 2 , G = a α 1 Δ x 2 + 1 − a α 3 Δ x 2 + b − 1 α 1 Δ z 2 + 1 − b α 3 Δ z 2 + ω 2 v 2 d 4 , H = a α 0 Δ x 2 − 2 a α 1 Δ x 2 + 1 − a α 2 Δ x 2 + 2 a − 1 α 3 Δ x 2 + b − 1 α 0 Δ z 2 − 2 b α 1 Δ z 2 + 1 − b α 2 Δ z 2 + 2 b α 3 Δ z 2 + ω 2 v 2 d 1 , I = − 2 a α 0 Δ x 2 + 2 a α 1 Δ x 2 + 2 a − 1 α 2 Δ x 2 + 2 1 − a α 3 Δ x 2 − 2 b a 0 Δ z 2 + 2 b − 1 a 1 Δ z 2 + 2 b a 2 Δ z 2 + 2 1 − b a 3 Δ z 2 + d 0 w 2 v 2 . (10)

Next, we solve for the weighting coefficients of the improved 25- point difference scheme. The dispersion relation is given by Equation 11 V ph v = { − { [ a 0 2 α + 2 1 − α C A − 1 + a 1 2 α + 2 1 − α C B − 2 A + 1 + a 2 4 α A − 1 C + 2 1 − α A D − D + A − 1 + a 3 [ 4 α B − 2 A + 1 C + 2 1 − α B D − 2 A D + D + B − 2 A + 1 ] ] + r 2 [ a 0 2 b + 2 1 − b A C − 1 + a 1 4 b C − 1 A + 2 1 − b B C − B + C − 1 + a 2 2 b + 2 1 − b A D − 2 C + 1 + a 3 4 b D − 2 C + 1 A + 2 1 − b B D − 2 B C + B + D − 2 C + 1 ] / [ 2 π G ( d 0 + 2 d 1 A + 2 d 2 C + 4 d 3 A C + 2 d 4 B + 2 d 5 D + 4 d 6 B C + 4 d 7 A D + 4 d 8 B D ) 1 / 2 . (11)

Similarly, the weighting coefficients can be obtained by solving the equations using the conjugate gradient method.

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

LSRTM based on safe Anderson acceleration type I.

Input: Choose an initial reflectivity model m 0 , memory parameter M ≥ 1 , regularization constant θ ¯ , τ , β ∈ 0 , 1 , safety protection constants D > 100 and 0 < ε < 0.01 , and fixed-point iteration operator f . In the first iteration, calculate the data residual L m 0 − d o b s , obtain the gradient ∂ E ∂ m = L T L m 0 − d o b s , update the reflectivity f H ∼ m 0 = m 0 − μ ∂ E ∂ m , and output f H ∼ m 0 as the reflectivity for the first iteration. Calculate g 0 = m 0 − f H ∼ m 0 , initialize the regularization function H 0 = I , set accumulation constants M A = N A = 0 , r ^ 1 = , and U ¯ = g 0 2 2 .

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.

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.

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.

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.

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.