Research on the Improved Least Squares Reverse Time Migration Imaging Method

. An improved Least Squares Reverse Time Migration (LSRTM) method is proposed in the paper, which can effectively improve convergence speed and imaging accuracy. Firstly, the key techniques in the implementation of LSRTM are discussed. Secondly, a condition factor is introduced in the iteration process of conjugate gradient method. Finally, the imaging effect and performance of the algorithm are analyzed. The experiment results indicate that it can speed up the convergence speed and improve the convergence accuracy, so as to improve the imaging effect. Compared with the conventional LSRTM, the data residual of improved LSRTM can be reduced by about 5%.


Introduction
With the deepening of oil and gas exploration, the least squares migration method based on the inversion idea is optimized to find the reflection coefficient model that best matches the observed seismic data, correct the amplitude errors in the migration results, and further improve the solution of seismic imaging, so as to better perform lithologic reservoir imaging and reservoir parameter inversion [1][2][3].
The idea of least squares migration inversion was first proposed by Lebras et al. [4], which is based on the method of data fitting to obtain the optimal imaging value through successive iterations. Nevertheless, the wavefield propagation operator cannot accurately predict the actual data, which will greatly increase the mount of iterative inversion calculation, thereby resulting in the convergence becomes quite difficult [5][6]. In order to improve the convergence performance and obtain higher imaging quality, a large number of scholars have carried out lots of academic research work in related fields. Dai et al. [7][8] proposed LSRTM for synchronized source data, which employs a defuzzification filter as a preconditioning operator to speed up the convergence. Li Zhenchun et al. [9], Liu Jinyu et al. [10], and FAN et al. [11] added constraints to least squares reverse time migration. Li Chuang [12] et al. adopted bidirectional lighting as a precondition and employed a fast high-pass filtering method to preprocess the gradient. Li Qingyang et al. [13] saved the calculation cost of least squares reverse time migration and improved the calculation efficiency by optimizing the phase encoding. Wang Liankun et al. [14] suppressed migration artifacts by constructing cross-correlation error functionals, thereby improving the fidelity and resolution of imaging results.
Fang Xiuzheng et al. [15] adopted inverse scattering imaging conditions to effectively suppress the lowfrequency noise in the gradient, so that the convergence speed during the inversion process was accelerated and the imaging quality was significantly improved. Tian Wenbin et al. [16] used a guided filtering operator to suppress the plane wave least squares migration artifacts, meanwhile, improved the convergence efficiency and imaging accuracy. Ke Xuan et al. [17] introduced GPU computing to least squares reverse time migration, which has a significant acceleration effect. Yang Jidong [18] incorporated anisotropy into the two-dimensional elastic wave equation and established an elastic least squares reverse time migration method for multi-component data in isotropic media. Fang Jinwei [19] applies the objective function based on convolution to LSRTM, and relaxes the requirement for accurate wavelet, thereby improving the robustness of the inversion in the presence of noise.
The current paper is devoted to a novel approach called the improved least squares reverse time migration imaging. The condition factor is added to the conjugate gradient, which to effectively settle the issues of illposed inversion and huge iterative calculation as well as improving imaging accuracy and convergence speed.

Basic principles of LSRTM
In order to obtain the best offset result matching the recorded data, the least squares migration problem is equivalent to optimizing a specific objective function, that is, to minimize the data errors between forward modeling and observation. Therefore, the objective function can be expressed as: Where obs d represents the observation data, m manifests the model parameter, and L is the inverse migration operator. Here, we employ reverse time demigration for data reconstruction [20], and the formula can be described as follows: Equation (2) is the wave field forward modeling, where 0 u is the background wave field. Equation (3)

Condition factor for conjugate gradient
Generally speaking, it is difficult to directly settle the objective function (9). Therefore, gradient-based solutions such as steepest descent method and conjugate gradient method are usually employed. Then the gradient expression can be written as follows: Where * L is the adjoint operator of L . In the actual solution process, the residuals of demigration data and observation data are obtained firstly, and then perform reverse time migration according to the residual of the data.
The calculation cost of least squares reverse time migration is extremely high, except for the larger calculation amount of the wavefield propagation operator of the two-way wave equation, another more important reason is the slower iterative convergence speed. Compared with the steepest descent method, the conjugate gradient method is an inversion iterative algorithm with faster convergence speed, and its search direction is a combination of the negative gradient direction and the previous search direction. The formula can be established as follows: Where k g is the gradient, k g is the search direction, and  is the correction factor of the search direction k g . Different conjugate gradient methods has different  . Famous conjugate gradient methods include HS (Hestenes-Stiefel), FR(Fletcher-Reeves), PRP(Polak-Ribiere -Polyak), DY(Dai-Yuan), etc. Among them, the convergence of PRP method is relatively weak and the numerical performance is better, while that of FR method is conversely. Therefore, this paper focuses on the PRP conjugate gradient method to obtain the solution of the objective function. By introducing parameters to adjusts the descent direction of the gradient by to meet the sufficient descent, thereby accelerating the convergence speed, saving calculation costs, and improving imaging accuracy [21]. The improvement formula can be obtained as follows: where PRP   is the modified conjugate gradient correction factor, k d is the conjugate gradient, k is the number of iterations, and k  is a self-adjusting parameter, which can automatically adjust the parameter  .

Experiments and Results
The size of the model is 501 sampling points horizontally and 301 sampling points vertically, as shown in Figure 2 (a). Seismic records are obtained by using time 2 order space 8 order finite difference forward modeling. In the forward simulation, the vertical and horizontal grid sizes were both 10m, the time step was 1ms, the recording duration was 2.5s, and the main frequency of wavelet was 20Hz. A total of 24 guns were recorded, starting at 100m and spacing was 200m. There were 501 receiving shots per gun and 2500 sampling points per shot.
To discuss the effectiveness of the improvement scheme, conventional LSRTM and Improved LSRTM were tested. As shown in Figure 1 (c), the offset velocity is smoothed by the true velocity. The imaging results are shown in Fig. 1 (d), Figure 1 (e) and Figure 1 (f), respectively. As shown in Figure 1 (d), the energy of conventional LSRTM imaging is not balanced, the source effect (red arrow) and noise are obvious in the upper layer, and the energy of deep area (blue arrow) is not balanced. Figure 1 (e) Comparing with Figure 1 (d), it shows that the source effect is significantly reduced and the noise suppression is obvious after Improved LSRTM. Furthermore, the energy in deep regions is compensated, and the effect is better than LSRTM, which proves the validity of the improved CG method. Figure 2 compares the convergence curves of conventional LSRTM and Improved LSRTM. It can be seen that the convergence speed of the residual is accelerated and can converge to a lower value, which verifies the effectiveness of the algorithm.

Conclusions
In this paper, an improved LSRTM is proposed. On this basis, through the imaging test of the model, the following conclusions are obtained： (l) Improved LSRTM method can improve the image quality and has better amplitude preservation and resolution in depth and on both sides, which is beneficial to the identification of deep structures; (2) Improved LSRTM method has faster convergence speed and better convergence accuracy than the conventional LSRTM in both simple and complex constructions.
Although the LSRTM with improved conjugate gradient method proposed in this paper can improve the imaging effect to some extent, there is still much room for further development to improve its computing efficiency and imaging accuracy. The improved method should be further optimized and combined with the multi-source method and coding technology to make the LSRTM image better and more efficient.