Least squares reverse time migration imaging with illumination preconditioned based on improved PRP conjugate gradients

Least squares reverse time migration (LSRTM) imaging is the one of the most accurate methods for migration imaging at present, and Polak–Ribiere–Polyak conjugate gradient (PRPCG) for LSRTM has the good numerical performance but weak convergence, so we construct an optimization factor to improve the iteration direction of the gradient, which can automatically generate a sufficient descent direction. The improved PRPCG (IPRPCG) can reduce the data residual values and speed up the iteration. And the illumination preconditioned (IP) operator is employed to IPRPCG-LSRTM which solves the problem of low resolution due to the insufficient iterative gradient information. In this paper, the experiments show that the imaging results of the proposed method (IPRPCG-IP-LSRTM) is improved greatly in detail characterization and events continuity, the iterative curve converged faster significantly, and the normalized data residual was reduced by 6.55% on average, which improved the accuracy of migration imaging effectively.


Least squares reverse time migration imaging with illumination preconditioned based on improved PRP conjugate gradients
Xiaodan Zhang 1* , Rui Li 1 , Lin Cui 1 , Dongxiao Liu 1 , Guizhong Liu 2 & Zhiyu Zhang 3 Least squares reverse time migration (LSRTM) imaging is the one of the most accurate methods for migration imaging at present, and Polak-Ribiere-Polyak conjugate gradient (PRPCG) for LSRTM has the good numerical performance but weak convergence, so we construct an optimization factor to improve the iteration direction of the gradient, which can automatically generate a sufficient descent direction.The improved PRPCG (IPRPCG) can reduce the data residual values and speed up the iteration.And the illumination preconditioned (IP) operator is employed to IPRPCG-LSRTM which solves the problem of low resolution due to the insufficient iterative gradient information.In this paper, the experiments show that the imaging results of the proposed method (IPRPCG-IP-LSRTM) is improved greatly in detail characterization and events continuity, the iterative curve converged faster significantly, and the normalized data residual was reduced by 6.55% on average, which improved the accuracy of migration imaging effectively.
With the development of oil and gas exploration, it poses greater challenges to traditional migration imaging methods for underground structure imaging 1,2 .In order to improve imaging accuracy and amplitude fidelity, researchers proposed least squares reverse time migration (LSRTM) 3,4 .LSRTM often uses Conjugate Gradient (CG) algorithm to find the optimal solution.The classical CG methods include the Fletcher-Reeves (FR) , Polak-Ribiere-Polyak (PRP) , Hestenes Stiefel (HS), and Conjugate Descent (CD) method.Their different iteration directions and step sizes will show different convergence rates 5,6 .Therefor, it is an important topic in the application research of LSRTM to research the new algorithm of iteration directions to achieve global convergence 7,8 .For the problem of insufficient illumination in deep subsurface imaging [9][10][11] , which leads to insufficient iterative gradient information resulting in low imaging resolution and slow convergence of iterative curves, illumination preprocessing has been developed and has attracted extensive research in academia [12][13][14] .
As for the development of conjugate gradient algorithm and illumination preprocessing, Liu et al. 15 performed time-domain full waveform inversion for eight versions of the CG method, and numerical experiments showed that the PRP conjugate gradient (PRPCG) method is more efficient.Kim et al. 16 proposed direction-oriented wavefield imaging, which compensates for possible illumination effects during acquisition and can produce highfidelity depth profiles and correct imaging of complex wavepaths.Yu et al. 17 developed a 3D nonlinear conjugate gradient traveltime inversion method using the PRP conjugate gradient method for solving constrained damping least squares problems, and model tests showed that the method can substantially improve the inversion resolution.Chen et al. 18 proposed the use of a wave illumination compensation method for the possible existence of offset shadows in seismic wave illumination, which can eliminate the offset imaging shadows and improve the computational efficiency.Wu et al. 19 introduced the reweighted regularized conjugate gradient method of PRP formula to constrain the inversion for the "skin effect" problem in gravity data inversion, which improved the computational efficiency and accuracy of the inversion.Sun et al. 20 developed a stabilization compensation operator that attenuates poor illumination of subsurface structures and does not amplify high frequency noise in the data.Hu et al. 21derived the gradient of the target generalized function on the model parameters from the Love wave fluctuation equation combined with the PRP conjugate gradient algorithm, and model tests verified that the method can improve the computational efficiency.Du et al. 22 used log-absolute error functions for AVO inversion and used a new spectral PRP conjugate gradient method in iterations to solve large-scale optimization problems, and then they combined a smooth nonconvex regularization method with adaptive individual weight gain and used a PRP conjugate gradient method to minimize the objective function 23 , recently they used a smooth L1 parametrization as the loss function and used a new spectral PRP conjugate gradient algorithm to optimize the inversion, and proposed a robust AVO inversion algorithm based on generalized nonconvex dictionary learning 24 .
The conventional PRPCG method is considered to have good stability and numerical performance, but it may fall into infinite loops near local minima during iteration, making the objective function unable to converge.Therefore, to ensure the inversion stability and improve the inversion results, this paper first proposed an improved PRPCG-LSRTM (IPRPCG-LSRTM) imaging method, which constructed an optimization factor to improve the iteration direction to make it fall sufficiently and can improve the iteration speed.Subsequently, a bidirectional illumination preconditioning (IP) operator was used for IPRPCG-LSRTM, and the IPRPCG-IP-LSRTM was proposed, which can solve the low-resolution problem caused by insufficient iterative gradient information.

Methods
Principles of LSRTM.LSRTM is based on traditional reverse time migration imaging and combines inversion ideas to perform migration imaging on seismic profiles.It is the migration result under the framework of least squares inversion, which corrects the migration imaging process and results.The migration principle is shown in Fig. 1: Figure 1 is the migration principle diagram under the constant velocity model.The profile that has not been migration is 0C1D1 , and the migration profile is 0CD .On the unmigrated profile 0C1D1 , the reflected wave gen- erate from point C directly below point E is observed at point A and plotted at point C1 .The true reflection point is in the upward tilt direction of its apparent position.After using the IPRPCG-IP-LSRTM method proposed in this article, the reflection layer segment 0C1D1 can be regressed to its real position CD.
Born approximation shows that: where s is the source point, x is the underground reflection point, and r is the detection point, D(r, t|s, 0) is the recorded seismic data model, m(x) is the reflectivity model, W(t) is the source, G(r, t|x, 0) is Green function from source to underground reflection point, G(r, t|s, 0) is Green function from underground reflection point to detection point.The matrix of Born approximation is expressed as: where L is the Born forward operator, also known as the reverse migration operator.Since L −1 is difficult to achieve, the conventional reverse time offset is simplified using L T instead of L −1 .The expression is as follows: where L T denotes the adjoint operator, also known as the migration operator.
This substitution simplifies the process of solving the reverse time migration, but it has an impact on the resolution and amplitude preservation of the migration imaging.To mitigate the negative effects of this effect, the idea of least square inversion is introduced to construct the objective function for fitting and minimizing the objective function: (1) where d obs is the observed seismic data.The minimization of Eq. ( 4) is to minimize the derivative of the objective function with respect to the model m: The equation can be solved as follows: where L T L is the Hessian matrix, forming the Hessian matrix is challenging, and computing the inverse of this matrix is also very difficult, especially for large 3D imaging applications.Therefore, an iterative algorithm is used to solve the above objective function.
In this paper, data reconstruction of simulated seismic data can be achieved by these equations as follows: Equation ( 7) is the wave field forward modeling, where v 0 (x) is the background velocity, u 0 (x s , x, t) is the background wave field, δ(x − x s ) is the Dirac function of the source.Equation ( 8) uses the background wave field and model parameter m as the secondary source, the disturbed wave field δu(x s , x r , t) is obtained.It can be seen that the solution of the disturbance wave field requires two finite difference simulations.The numerical simulation process of the Finite difference method is shown in Appendix 1. Equation ( 9) records the disturbance wave field at the detection point x r to obtain simulated data d(x s , x r , t).

Principles of IPRPCG-IP-LSRTM.
In solving the objective function Eq. ( 4), to avoid directly solving the inverse of the Hessian matrix, the gradient method is often used for multiple indirect iterations.The conjugate gradient method has the characteristics of a small amount of calculation, high accuracy, good stability, and low storage requirements.Due to the different conjugate gradient parameters β , it can be divided into FR, PRP, HS, CD, etc.Among them, PRPCG is recognized as one of the best numerical performance methods 25 .When the algorithm generates a small step, the gradient of this objective function is close to the gradient of the previous objective function, and the search direction will automatically approach the negative gradient direction, which can effectively avoid the generation of successive small steps.Therefore, this article uses the PRPCG to find the least squares optimal solution, and the calculation formula for the parameter β is: where g k = ∇f (x k ) represents the gradient operator of the objective function f (x) .However, the convergence of the PRPCG is relatively weak, and it may loop infinitely far from the global optimal solution, and produce a search direction that makes the objective function rise resulting in diverging scenario and stopping further computations 26 .Therefore, a parameter β containing an optimization factor σ is constructed in this paper to solve this problem.This factor can modify the iteration direction of β to reduce the residual values and improve the iteration speed, as shown in Eq. ( 11): To decrease the value of the β k , the optimization factor σ k must satisfy σ k > 1 .so the construction of σ k is shown as follows: In Eq. ( 12), µ k is defined as the limiting factor of σ k .The limiting factor µ k is related to the current gradient and the last gradient, and we construct a function µ k , which is given as µ k = f (g k−1 , g k ) .According to the theory of global convergence and the direction of the fastest descending gradient, the expression of µ k is constructed as Eq. ( 13): where τ is a constant greater than or equal to one.The limiting factor µ k is the fine-tuning quantity of β k , which modifies the optimization factor σ k through the last gradient value and the current gradient value based on ensuring the decline of β k .Then, the optimized conjugate gradient parameter β k is given as follows: , the following equation holds: It can be seen that the optimization factor σ k can automatically adjust parameter , that is, the constructed σ k can adjust the descent direction of the gradient to automatically generate a sufficient descent direction, search along the feasible descent direction, find the feasible point to lower the objective function, determine the appropriate moving step and accelerate the convergence speed of the objective function.
In LSRTM, the gradient is the reverse time migration imaging result of the residual between simulated seismic data and observed seismic data, which also requires two finite difference calculations to calculate the background and residual wave fields: where q(x s , x, t) is the residual wave field, the background wave field is a forward continuation in the time direc- tion, and the residual wave field is a backward continuation in the time direction.Therefore, the gradient can be obtained using the following formula: The calculation of gradient is the reverse time migration of data residuals, which can be expressed in the following matrix form: The search direction generated by the IPRPCG satisfies the sufficient descent condition and the global convergence theorem.The convergence analysis of the algorithm is shown in Appendix 2. The specific process of using IPRPCG to find the optimal solution for least-squares reverse time migration (IPRPCG-LSRTM) is shown in Table 1.
IPRPCG not only retains the excellent performance of traditional PRPCG, but also makes the search direction satisfy the condition that it is always the descent direction of the objective function, which has full descent, good global convergence performance and numerical performance.The accurate descent direction can make the data residuals smaller, so it has a great improvement in the computational accuracy and imaging resolution.
As the conventional migration operator is only the transpose of the forward operator, the illumination of some underground structures, especially deep and high-speed salt bodies, will be uneven, or even unable to image 27 .The Hessian operator can characterize the distribution of underground lighting energy, and its Inverse matrix ( 14) Table 1.The solution process of IPRPCG-LSRTM.
The input: m = 0 : Initial value of reflection model ϕ > 0 : Threshold for iterative data residuals, or k max : The maximum number of iterations, its initial value k = 0 The output: m : Image results of LSRTM Step 1: Obtain the gradient g of the objective function: Step 2:Take conjugate gradient optimization parameters β , conjugate gradient d and the iteration step α: Step 3: Update the model: .However, it is difficult to calculate the inverse of the full Hessian matrix under the current computing power.Considering that the Hessian matrix is a matrix with dominant principal diagonal elements, its diagonal elements can be used to approximate the Hessian operator, and then its Inverse matrix can be used as a precondition operator in the implementation of LSRTM 29 , which can also improve the imaging resolution.Therefore, this paper introduces the Illumination precondition operator P = H −1 in the optimization parameters β k and conjugate gradient d k , whose approximation matrix is : where H(a) is the approximately linear part of the Hessian operator, Re{} means taking the real part for the com- plex numbers in parentheses, w is the angular frequency, x s denotes the location of the source point, r denotes the location of the detector point, x denotes the location of any point in the subsurface, f s (w) is the frequency domain expression of the source , G(x s , x, w) is the Green's function of the source, and G(x, r, w) is the Green's function of the checkpoint.
Then the optimized parameters ε k and conjugate gradient h k with IP are shown follows: Compared with the CG method, the proposed method IPRPCG-IP only adds a precondition operator and some dot product calculation at the beginning of the calculation, which can be ignored compared with LSRTM, but it can greatly improve the migration imaging effect.The overall implementation flow chart of IPRPCG-IP-LSRTM proposed in this paper is shown in Fig. 2.

Experimental results and analysis
Noise resistance test.Sigsbee2B truncation model was used for testing the noise resistance of PRPCG-LSRTM, IPRPCG-LSRTM, and IPRPCG-IP-LSRTM imaging.The velocity models were shown in Fig. 3. Figure 3a was the real velocity model, and Fig. 3b was the background velocity model, which was obtained by real velocity model smoothing.The real reflection coefficient model was shown in Fig. 3c, and the detailed test parameters were shown in Table 2.
PRPCG-LSRTM, IPRPCG-LSRTM and IPRPCG-IP-LSRTM imaging methods were employed on Sigsbee2B truncation model, and the three methods had been used for 30 iterations.The imaging results were shown in Fig. 4, where Fig. 4a was the PRPCG-LSRTM imaging results, Fig. 4b was the IPRPCG-LSRTM imaging results, and Fig. 4c was the IPRPCG-IP-LSRTM imaging results.From the above three migration imaging results, it can be seen that the PRPCG-LSRTM imaging energy was unbalanced and more serious low-frequency noise still remained.Comparing with PRPCG-LSRTM, IPRPCG-LSRTM method imaging of Fig. 4b, the noise in the imaging profile was suppressed, the regional energy was compensated, and the tomographic part was clearer, and the effect was better than that of PRPCG-LSRTM, which proved the effectiveness of IPRPCG-LSRTM.IPRPCG-IP-LSRTM imaging of Fig. 4c compared with IPRPCG-LSRTM of Fig. 4b, the low-frequency noise was further suppressed, the profile amplitude was more balanced, the signal-to-noise ratio was higher, the wavefield energy was more balanced, and the imaging quality was improved significantly, which can be seen in the red and blue circle of Fig. 4a-c.
To evaluate the image quality of the experimental results more rigorously, Table 3 quantitatively compared the imaging results of the sigsbee2B truncated model using Peak Signal to Noise Ratio (PSNR) and Structural Similarity (SSIM).From the data in the table, it can be seen that the imaging results of IPRPCG-LSRTM iteration 30 times increased PSNR and SSIM values compared to PRPCG-LSRTM iteration 30 times, indicating a significant improvement in imaging quality.However, IPRPCG-IP-LSRTM iteration 30 times has the highest PSNR and SSIM values, indicating the highest image quality generated.
It can be seen from the Fig. 5 that the three methods were stable from the 15th iterations.IPRPCG-IP-LSRTM (red) converged the fastest in the three methods, which improved the computational efficiency.The residual value of IPRPCG-IP-LSRTM was the smallest, so the imaging result was the closest to the real structure, the accuracy of IPRPCG-IP-LSRTM imaging was the best.
Table 4 showed the comparison of the program running time of PRPCG-LSRTM, IPRPCG-LSRTM, and IPRPCG-IP-LSRTM with 30 iterations.It can be seen that the program running time of IPRPCG-LSRTM with 30 iterations was the shortest, which improved the computational efficiency by about 24.73% compared to PRPCG-LSRTM.After adding IP, the program running time of IPRPCG-IP-LSRTM was extended compared to IPRPCG-LSRTM due to the addition of a preconditioned operator and some point multiplication calculations at the beginning of the calculation.However, compared to PRPCG-LSRTM, the program running time of IPRPCG-IP-LSRTM with 30 iterations was significantly reduced, and the computational efficiency was improved by about 23.59%.
To compare the amplitude fidelity of PRPCG-LSRTM, IPRPCG-LSRTM, and IPRPCG-IP-LSRTM, singlechannel imaging data at a distance of 2400m (indicated by the yellow dashed line) were extracted from the profiles shown in Figs.3c, and 4a-c.for comparison, as shown in Fig. 6.From the figure, it can be seen that the effect of IPRPCG-LSRTM (green) was better than that of PRPCG-LSRTM (blue), but they only had good amplitude fidelity in the shallow layer, and there was a significant difference in the true reflection coefficient (black dashed line) in the middle and deep layers.And the amplitude of IPRPCG-IP-LSRTM (red) was closer to the real reflection coefficient, especially in the middle and deep layers, indicating that IPRPCG-IP-LSRTM could better compensate for the amplitude and illumination of energy attenuation caused by geometric diffusion or absorption attenuation in the middle and deep layers, thus it had higher amplitude preservation.
Artificially adding random noise to Sigsbee2B truncation model to form a data volume with low signal-tonoise ratio, which can be seen in Fig. 7.
Figure 7a was the original shot data and Fig. 7b was the one of the shots with low signal-to-noise data.PRPCG-LSRTM and IPRPCG-IP-LSRTM were employed on the low signal-to-noise data.
Figure 8 was the imaging results of the Low signal-to-noise ratio shot data of Sigsbee2B.Figure 8a was the imaging results of PRPCG-LSRTM with 30 iterations, and Fig. 8b was the imaging results of IPRPCG-IP-LSRTM with 30 iterations.We can be seen that both contained migration artifacts caused by the noise, but the low-frequency noise was suppressed better by IPRPCG-IP-LSRTM in Fig. 8b, and the imaging energy was more balanced, and the imaging resolution was higher.
Accuracy and convergence test.Marmousi model was employed to test the convergence and accuracy of RTM, PRPCG-LSRTM, and IPRPCG-IP-LSRTM proposed in this paper, which can be shown in Fig. 9. Figure 9a was Marmousi's real velocity model.The velocity range was shown as the scale.Figure 9b was the background velocity model obtained by Marmousi's real velocity model smoothing.The detailed test parameters were shown in Table 5.
The migration imaging results of the above methods were shown in Fig. 10. Figure 10a showed the RTM imaging results of Marmousi.It can be seen that there was the serious low-frequency noise in the imaging results, the shallow imaging was week and the deep area imaging was not clear.Figure 10b showed the imaging result   www.nature.com/scientificreports/after Laplacian filtering based on RTM, and it can be seen from the figure that there was the serious noise in shallow surface area.Figure 10c-f showed the migration imaging results of 10, 20, 30 and 40 iterations based on PRPCG-LSRTM.It can be seen that with the increase in the number of iterations, the imaging result become better and better, the source effect decreased, the energy become more balanced, the imaging on both sides and deep layers become clearer, the event energy was restored, the noise was suppressed, and the imaging resolution was higher.10f, the imaging energy is more uniform and the imaging of details was clearer.And it can be seen from Fig. 10f, g that the IPRPCG-IP-LSRTM migration imaging result with 30 iterations was better than PRPCG-LSRTM with 40 iterations, especially the areas in the red boxes.
Figure 10g showed the imaging results of IPRPCG-IP-LSRTM with 30 iterations when the parameter τ = 4 , and Fig. 10h showed the imaging results of IPRPCG-IP-LSRTM with 30 iterations when the parameter τ = 2 .After experiments with parameter τ in the optimization factor σ , it was found that the best results are obtained when parameter τ = 4 of Marmousi model.Comparing Fig. 10g, h, it can be seen that the imaging of depth boundary with less data and deeper layers of Fig. 10g was clearer, with more balanced energy and higher resolution.
We compared the magnified details in the red box (1) of Fig. 10f, g, they were shown in Fig. 11. Figure 11a showed the enlarged detail image of the imaging results in the red box (1) of Figs.10f, and 11b showed the enlarged image of the imaging results in the red box (1) of Fig. 10g.Magnifying the imaging details in the red circle, it can be seen that the geological structure (horizon) of PRPCG-LSRTM had a fracture and no imaging, while the layer structure (event) was continuous and clear in Fig. 11b.It can be seen that the IPRPCG-IP-LSRTM method proposed in this paper can improve the resolution of the migration imaging result effectively.
Table 6 showed the image quality evaluation of the Marmousi model's imaging results using PSNR and SSIM indicators.It can be seen from the table that the image quality of RTM had improved after Laplace filtering, and as the number of iterations increases from 10 to 40, the imaging effect of PRPCG-LSRTM also improved.The imaging quality of the IPRPCG-IP-LSRTM proposed in this article was already better at 30 iterations than   Figure 12 showed a comparison of the convergence curves of normalized data residuals using the Barzilai Borwein algorithm for least squares reverse time migration (Barzilai-Borwein-LSRTM) based on the Marmousi model, compared to PRPCG-LSRTM and IPRPCG-IP-LSRTM with 40 iterations.From the figure, it can be seen that the normalized data residuals of Barzilai-Borwein-LSRTM decreased rapidly during 1-5 iterations, and the effect of Barzilai-Borwein-LSRTM was similar to that of PRPCG-LSRTM during 10-25 iterations.At 26 to 40 iterations, the normalized data residuals of Barzilai-Borwein-LSRTM were smaller than those of PRPCG-LSRTM.However, from beginning to end, the normalized data residuals of Barzilai-Borwein-LSRTM and PRPCG-LSRTM were not lower than those of IPRPCG-IP-LSRTM.This meant that the convergence accuracy of the IPRPCG-IP-LSRTM proposed in this paper was the highest, and the imaging data was closer to real data.
Table 7 showed the comparison of the computational efficiency of IPRPCG-IP-LSRTM and PRPCG-LSRTM.From the data in the table, it can be seen that the normalized residual values of IPRPCG-IP-LSRTM decreased faster than PRPCG-LSRTM, with an average of about 6.55%.Additionally, the program running time of IPRPCG-IP-LSRTM was shorter than PRPCG-LSRTM, with an average improvement in computational efficiency of about 22.45%.This demonstrated the superiority and effectiveness of IPRPCG-IP-LSRTM in terms of imaging accuracy and descent speed.
For observing the imaging details, single-trace data of migration imaging was extracted.Figure 13 showed the single-channel amplitude comparison between the imaging results of the 30th iteration of PRPCG-LSRTM, IPRPCG-IP-LSRTM, and the real reflection coefficient at a distance of 1200 m.From the figure, it can be seen that the imaging effect of IPRPCG-IP-LSRTM was better than PRPCG-LSRTM, especially when there was a high-speed layer at a depth of 2400-2600 m, the amplitude curve of IPRPCG-IP-LSRTM could better match the true reflection coefficient.The amplitude energy of IPRPCG-IP-LSRTM was stronger than that of PRPCG-LSRTM.At the depth of 2700-3000 m, the amplitude energy of PRPCG-LSRTM was greatly reduced, while the amplitude energy of IPRPCG-IP-LSRTM was still very strong, closer to the amplitude energy of the real reflection coefficient, as shown in the green dashed circle in Fig. 13.www.nature.com/scientificreports/

Discussions
PRPCG had not been widely used because of its good numerical performance and relatively weak convergence, which may produce a search direction that makes the objective function rise and thus fail to converge.With the work done, we proposed a method that improved the PRPCG and incorporated the IP operator for the combination.Compared with earlier work, Zhao et al. 30 used the original PRP conjugate gradient method to solve the sparse joint inverse objective generalization, and we solved the problem that the PRP method may not converge globally.Regarding the earlier finding of Jin et al. 31 that effective tectonic information could not be obtained due to uneven illumination when performing offset imaging in complex tectonic zones in the subsurface, we used IP operators to compensate for illumination and improve image resolution.In addition, recent work by Rong et al. 32 showed that it was possible to use deep learning to obtain low-illumination regions of geological models and recovered complex geological formations in the region.

Conclusions
The advantages of IPRPCG-IP-LSRTM proposed in this paper were that it can effectively improve the migration imaging resolution and computational efficiency.The method started from the gradient descent direction, constructed the optimization factor, improved the conjugate gradient parameters, added the illumination preconditioning operator, accelerated the convergence speed, reduced the normalized data residual values, and effectively improved the computational efficiency and imaging accuracy.
According to the experiments in this paper, the method IPRPCG-IP-LSRTM proposed can reduce the average residuals of normalized data by about 6.55%, which verified the effectiveness and feasibility of the algorithm.The disadvantage was that due to the huge computational power of LSRTM, a suitable number of iterations should be selected when using the IPRPCG-IP-LSRTM method in practical applications.Finally, combining our method with deep learning and intelligent optimization was our further research.

Appendix 1: Numerical simulation of finite difference method
The two-dimensional constant density acoustic Wave equation can describe the propagation of acoustic waves in Seismic waves in underground media: where u is the wave field function u(x, z, t) , v is the medium velocity v(x, z) , x and z are the spatial coordinate components, x represents the horizontal direction, z represents the vertical direction, and t is the temporal coordinate component.
The finite difference method is used to Discretization the numerical solution of Eq. ( 28).The wave field space should first be Discretization.Therefore, the spatial model needs to be meshed, as shown in Fig. 14: If the grid spacing after Discretization in x and z directions is x and z , and t is the time sampling step, x = i x , z = j z , t = n t can be obtained.Where i , j , and n are integers.Assuming u k i,j represents the wave field value at time k at point (i, j) , expand (x, z) at time k using Taylor's formula: Expand u k−1 i,j at time k at point (i, j) using Taylor's expansion: Add Eqs. ( 29) and (30) above, omit high-order small quantities, and organize the second-order time partial differential of point (i, j) and time k as follows: Similarly, the second-order space partial differential at point (i, j) and time k can be obtained as follows: The above equation replaces the derivative of the second-order partial differential with the difference of the wave field values of discrete grid points.By substituting Eqs.(31), (32), and (33) into Eqs.( 28), (34) can be obtained: Similarly, the difference scheme of the second-order Time derivative and the space 2N order acoustic equa- tion is as follows: Assuming that the objective function f (x) is lower bounded on the level set Ŵ = x ∈ R n f (x) ≤ f (x 1 ) , it is known that the sequence f k is monotonically bounded and thus converges by combining with Eq. (45), i.e., lim k→∞ f k+1 is a constant.Assume that f (x) is differentiable in some neighborhood of Ŵ and that ∇f (x) satisfies the Lipschitz condition, i.e., there exists ∀x, y ∈ and a constant Z > 0 such that: g(x) − g y ≤ Z x − y .Combining Eq. ( 46
with a main frequency of 20 Hz excited a total of 60 shots, with the first shot at a position of 38.1 m and a shot interval of 76.2 m.Each shot received a total of 600 smoothed background velocity model was shown in Fig. 3b for details Objective function solving IPRPCG-IP method Vol:.(1234567890)Scientific Reports | (2023) 13:13623 | https://doi.org/10.1038/s41598-023-40578-8

Figure 5 .
Figure 5.Comparison of residual convergence curves for 30 iterations of data.

Figure 6 .
Figure 6.Comparison of single-channel imaging curves at a distance of 2400 m.

Figure 7 .
Figure 7.The 30th shot data without and with noise.

Figure 13 .
Figure 13.Comparison of single-channel imaging curves at a distance of 1200 m.

Table 2 .
Test parameters for sigsbee2B truncation model.

Table 3 .
Image quality evaluation of sigsbee2B truncation model imaging results.Figure 10g showed the imaging results of IPRPCG-IP-LSRTM with 30 iterations.Compared with PRPCG-LSRTM imaging results with 40 iterations of Fig.

Table 4 .
Comparison of program running times of different algorithms.

Table 5 .
Test parameters for Marmousi model.

Table 6 .
Image quality evaluation of Marmousi model imaging results.

generated by different methods Evaluation indicators PSNR/dB SSIM
Comparison of residual convergence curves for normalized data.

Table 7 .
Comparison of computational efficiency.