The constant density acoustic wave equation is ( ∇ 2 − 1 v 2 ( X ) ∂ 2 ∂ 2 t ) p = δ ( X- X s ) δ ( t ) (1) where p is the wave field function, v(X) is the velocity at the X position, and the X _s and X are the position of the source and the geophone, respectively.

The actual velocity field can be composed of normal field and disturbance: 1 v 2 ( X ) = 1 v 0 2 ( X ) ( 1 − α ( X ) ) (2)

By expanding Taylor’s Eq. 2 at v ₀ and removing the higher order term, we obtain the following result: 1 v 0 2 − 2 Δ v v 0 3 = 1 v 0 2 ( 1 − α ) (3) where α = 2 Δ v v 0 , α(X) is the disturbance. After applying Eq. 2 into (1), and transforming the equation into Fourier frequency domain, we have [ ∇ 2 + ω 2 v 0 2 ( 1 − α ( X ) ) ] P ( X , X S , ω ) = δ ( X − X S ) (4) where P ( X , X S , ω ) = ∫ − ∞ ∞ p ( X , X S , t ) e − jωt dt .

By expanding Eq. 4, ( ∇ 2 + ω 2 v 0 2 ) P ( X , X S , ω ) = δ ( X − X S ) + ω 2 v 0 2 α ( X ) P ( X , X S , ω ) (5)

The total observed wavefield is the sum of the incident field and the scattered field: P ( X , X S , ω ) = P 0 ( X , X S , ω ) + P S ( X , X S , ω ) (6)

Applying Eq. 6 into (5) and decomposing it into two formulas: ( ∇ 2 + ω 2 v 0 2 ) P 0 ( X , X S , ω ) = δ ( X − X S ) (7) ( ∇ 2 + ω 2 v 0 2 ) P S ( X , X S , ω ) = ω 2 v 0 2 α ( X ) P ( X , X S , ω ) (8) where P ₀ is the background field and P _s is the scattering disturbance field. The total wavefield can be written as P ( X , X S , ω ) = P 0 ( X , X S , ω ) + ∭ V P 0 ( X , X ′ , ω ) ω 2 v 0 2 α ( X ) P ( X , X ′ , ω ) dX ′ (9) where V is the target region of velocity variation. Eq. 9 is the Lippmann–Schwinger integral formula. Assuming α is small, the scattering wave field under Born approximation (Liu, 2008) is expressed by: P S ( X , X S , ω ) = ∭ V P 0 ( X , X ′ , ω ) ω 2 v 0 2 α ( X ) P ( X , X ′ , ω ) dX ′ (10)

Defining m ( X ) = 2 Δ v v 0 3 to replace α(X) in Eq. 10, the wavefield in Eq. 10 can be obtained from Eqs. 11 and 12 as ( ∇ 2 + ω 2 v 0 2 ) P 0 ( X , X S , ω ) = δ ( X , X S ) (11) ( ∇ 2 + ω 2 v 0 2 ) P S ( X , X S , ω ) = ω 2 m ( X ) P 0 ( X , X S , ω ) (12)

The Born forward operator is represented by vector matrix: d = Lm (13) where m is the matrix form of migration profile or reflection coefficient model, d is the matrix form of simulation data, and L is the Born approximate forward operator matrix. The calculation of scattering wavefield can be obtained by forward simulation of Eqs. 11 and 12.

Conventional RTM can be expressed as m 0 ( x , z ) = L T D (14) where m ₀ is the RTM profile, and L ^T is the approximate migration operator. There are some errors when replacing the migration operator with the transpose of the forward modeling operator. To minimize the difference between the simulated data and the field data, the error function is defined as f ( m ) = 1 2 ‖ Lm − D ‖ 2 (15)

After taking partial derivative with respect to m , g = ∂ f ( m ) ∂ m = L T ( Lm − D ) (16)

When the gradient g is zero, the optimal solution of the least-squares problem is obtained: m = ( L T L ) − 1 L T D (17) where L ^T L is the Hessian matrix. Because it is so large and difficult to obtain, the gradient is gradually close to zero by iteration to avoid getting the inverse of Hessian matrix.

The cross-correlation conjugate gradient method for solving Eq. 15 can be expressed as (Huang et al., 2016): g ( k + 1 ) = L T [ L m ( k ) − D ] β = g ( k + 1 ) g ( k + 1 ) g ( k ) g ( k ) z ( k + 1 ) = g ( k + 1 ) + β z ( k ) α = [ z ( k + 1 ) ] T g ( k + 1 ) [ L z ( k + 1 ) ] T L z ( k + 1 ) m ( k + 1 ) = m ( k ) − α z ( k + 1 ) (18)

The gradient g expansion based on the steepest descent method can be expressed as g = ∑ t ∂ S s ( t , x , z ) ∂ t R res ( t , x , z ) (19) where β is the correction factor of the conjugate gradient method, α is the update step, z is the conjugate gradient, g is the steepest descent gradient, D is the field data, S s is the source forward wave field, and R res is the backward wavefield of the residual error between Born approximate forward data and field data.

It can be seen from Eq. 18 that the difficulty of conjugate gradient LSRTM method lies in gradient calculation and Born approximate forward modeling. d res = L m ( k ) − D (20) d _res is the residual of Born approximate forward data and field data. The calculation of the steepest descent gradient is similar to the cross-correlation of conventional RTM, except that conventional RTM is the cross-correlation of field data and forward wavefield, while the steepest descent gradient is the cross-correlation of the forward modeling wave field. In the conventional cross-correlation RTM, because of the use of an imprecise migration operator, the unbalanced wave field energy affects the imaging results. When it is close to the source and geophone with strong energy, the signal may become blurred (Yang et al., 2018). To weaken the influence of energy imbalance, the imaging conditions of source-normalized cross-correlation RTM for compensating underground illumination were proposed. Considering that the preconditioner in LSRTM of conjugate gradient method is difficult to obtain and cannot approximate the Hessian matrix well, the source effect has a great impact on the gradient. Here, the normalization is used to improve the calculation process of the steepest descent gradient and weaken the influence of the source effect.

The source-normalized cross-correlation imaging condition of the RTM source is m 0 ( x , z ) = ∑ t S s ( t , x , z ) R s ( t , x , z ) ∑ t ( S s ( t , x , z ) ) 2 (21) where R _S is the reverse propagation field of seismic record.

Likewise, the source-normalized steepest descent gradient formula of LSRTM is g ˜ ( k + 1 ) = ∑ t ∂ S s ( t , x , z ) ∂ t R res ( t , x , z ) ∑ t ( ∂ S s ( t , x , z ) ∂ t ) 2 (22) ∂ S s ( t , x , z ) ∂ t is the first-order partial derivative of source forward propagation field, which is very important for obtaining the zero-phase imaging profile (Yao and Wu, 2015). Eq. 18 can be written as g ˜ ( k + 1 ) = L T [ L m ( k ) − D ] ‖ L T [ L m ( k ) − D ] ‖ s β ˜ = g ˜ ( k + 1 ) g ˜ ( k + 1 ) g ˜ ( k ) g ˜ ( k ) z ˜ ( k + 1 ) = g ˜ ( k + 1 ) + β ˜ z ˜ ( k ) α ˜ = [ z ˜ ( k + 1 ) ] T g ˜ ( k + 1 ) [ L z ˜ ( k + 1 ) ] T L z ˜ ( k + 1 ) m ( k + 1 ) = m ( k ) − α ˜ z ˜ ( k + 1 ) (23) where β ˜ is the normalized conjugate gradient correction factor, α ˜ is the update step, z ˜ is the normalized conjugate gradient, and g ˜ = L T [ Lm − D ] ‖ L T [ Lm − D ] ‖ s is the source-normalized steepest descent gradient. The LSRTM imaging is realized through the iterative calculation of the above equation.

The real velocity model and the Gaussian smoothed velocity model are shown in Figures 1A,B, which is used for RTM and LSRTM background velocity field. The lateral length of the model is 1,250 m and the depth is 1,000 m. It consisted of 250 traces with a trace interval of 5 m and 11 shots with a shot interval of 125 m. The receivers are fixed, and the shot point moves at equal intervals.

It can be seen that conventional cross-correlation RTM results contain low-frequency noise (Figure 2A), coupled with strong near-surface energy caused by the source effect, resulting in unbalanced wave field energy, weak deep illumination, and unclear imaging. The LSRTM can solve the problems in RTM imaging. With an increase in the number of iterations, the data of Born forward simulation are gradually approaching the field data, and the migration profile is also gradually close to the reflection coefficient profile. As can be seen from the red arrows in Figures 2B,C, after 20 iterations, the noise in the shallow part of the migration profile gradually disappears, and the energy of the profile becomes more balanced, yielding clearer deep structure imaging. Under the same iteration times, the results of LSRTM processing by conjugate gradient normalization method are better than that by the conventional conjugate gradient method.

To eliminate low-frequency noise, Laplace filtering is performed on the profiles processed by RTM, cross-correlation conjugate gradient LSRTM, and conjugate gradient normalized LSRTM. From the comparison of red arrows in Figures 2D–F, it can also be seen that the amplitude of migration profile obtained by conjugate gradient LSRTM is more balanced than that obtained by RTM, and the imaging results of deep structure are better than that processed by RTM, while the migration profile obtained by conjugate gradient normalized LSRTM is better than that processed by the cross-correlation conjugate gradient LSRTM and RTM in both amplitude equalization and deep illumination (yellow arrow).

To verify the applicability of the conjugate gradient normalized LSRTM for complex model, imaging experiments were carried out on the Marmousi model (Figure 3). The velocity model and the Gaussian smoothed velocity model are shown in Figures 3A,B. The lateral length of the model is 4,000 m and the depth is 2,495 m. It consisted of 650 traces with a trace interval of 5 m and 14 shots with a shot interval of 250 m. The receivers were fixed, and the shot point moved at equal intervals. The seismic records were obtained using finite-difference forward modeling (Figure 3C shows the single-shot record).

The experimental work of the complex model was carried out on a workstation using the Intel(R) Xeon(R) Silver 4210R CPU @ 2.40 ghz, 128 GB memory, 64-bit operating system, and an X64-based processor. The graphics card was NVIDIA GeForce RTX3090 with 24 GB of video memory. In this computing environment, both the conventional LSRTM and conjugate gradient normalized LSRTM take approximately 3 min per iteration.

The real and Gaussian smoothed velocity field of Marmousi are shown in Figures 4A,B, which is used for the background velocity field of the RTM and LSRTM migration. From the migration results shown in Figures 4A,B, both cross-correlation RTM and normalized cross-correlation RTM will be contaminated by low-frequency noise (indicated by the yellow arrow). The low-frequency noise in the shallow part of the cross-correlation RTM is more serious. Normalized cross-correlation RTM is better in suppressing the noise in the shallow part, but there will be a small amount of low-frequency noise in the deep part. The main reason for this is that although the source-normalized imaging conditions suppress the shallow strong energy and enhance the deep illumination, they also enhance the wave field energy in the continuation process. The LSRTM can eliminate the low-frequency noise very well. It can be seen from Figures 4C,D (yellow arrow) that the low-frequency noise in the shallow part is obviously eliminated by least-squares processing, and as the number of iterations increases, the noise will continue to weaken. By comparing the low-frequency noise suppression results of the two imaging conditions, under the same number of iterations, the conjugate gradient normalized LSRTM is more significant for low-frequency suppression. In contrast, strong shallow low-frequency noise affects the imaging of underground structures by RTM. Figures 4E,F (red arrow) clearly show that reverse-time migration under different imaging conditions does not achieve accurate imaging of underground structures, some events cannot reflect accurate structure information well, and there are residual low-frequency noises in shallow parts. After the LSRTM processing, it can be seen from Figures 4G,H (red arrows) that the LSRTM can well eliminate low-frequency noise and realize accurate imaging of underground structures, and the overall amplitude of the profile is more balanced. Compare the underground illumination of the LSRTM under two different imaging conditions, it can be seen from the two enlarged images Figures 4I,J (corresponding to the two red dashed rectangular boxes in Figures 4G,H), the conjugate gradient normalized LSRTM is significantly stronger for the deep illumination than the result of the cross-correlation LSRTM. From the comparison of the whole section, we can also see that the conjugate gradient normalized LSRTM is better than the conventional LSRTM for noise suppression in the shallow part, illumination of the deep part, and imaging of the whole structure.

The spectrum of Figure 4G is shown in Figure 5A, and the spectrum of Figure 4H is shown in Figure 5B. Comparing Figures 5A,B, it can be seen that the conjugate gradient normalized LSRTM has a wider frequency band and more information. This also shows the superiority of conjugate gradient normalized LSRTM in the spectrum.

To compare the convergence speed of conjugate gradient normalized LSRTM and conventional conjugate gradient cross-correlation LSRTM, Figure 6 shows the comparison between the seismic records of Born approximate simulation and the observation records under the same number of iterations for the same shot in the simple model.

Figure 6A shows the sixth observation seismic record. The sixth shot simulation records of Born approximate forward modeling after 20 iterations of conventional conjugate gradient cross-correlation LSRTM, and conjugate gradient normalized LSRTM are shown in Figures 6B,D. We calculated the difference between the simulation record and observation record of the two LSRTM methods (Figures 6C,E). From the comparison of Figures 6C,E, it can be seen that under the same number of iterations, the Born approximate forward modeling record using the conjugate gradient normalized LSRTM is closer to the observation.

The normalized error reduction curve of the observation record and the simulation is shown in Figure 6F. The abscissa is the number of iterations of the two LSRTM methods, and the ordinate is the relative error. The error formula can be expressed as E r = | R o − R s | R o , (24) where E _r is the relative error, R _o denotes the observation record, and R _s denotes the simulation record. Figure 6F shows that compared to the conventional LSRTM, the conjugate gradient normalized LSRTM converges faster, and the residual error will eventually converge to a lower level. Therefore, it can be seen from the above results that the conjugate gradient normalized LSRTM converges faster than the conventional LSRTM, and its residual error is smaller after 20 iterations.

To test the adaptability of the method to the field data, two-dimensional land-based real data were imaged using the conjugate gradient normalized LSRTM. The velocity model of the field data and single-shot record are presented in Figures 7A,B. The velocity field of these data has a horizontal length of 32,050 m and a depth of 7,560 m. The data consisted of 328 shots with the shot interval of 100 m. The time sampling was 6 s with the sampling rate at 4 ms. Figure 8A shows the result of conventional normalized cross-correlation RTM filtering of field data, and Figure 8B shows the result of 20 iterations of conjugate gradient normalized LSRTM filtering. The conjugate gradient normalized LSRTM is superior to the conventional normalized cross-correlation RTM method, in both the suppression of shallow low-frequency noise and the illumination of deep structures. The energy of the shallow and deep layers of the method in this study is more balanced, and the profile imaging is better, especially for the middle and shallow imaging (as shown in the red box). Thus, it is better than the conventional normalized cross-correlation RTM method.