An Effective Acoustic Impedance Imaging Based on a Broadband Gaussian Beam Migration

: The estimation of the subsurface acoustic impedance (AI) model is an important step of seismic data processing for oil and gas exploration. The full waveform inversion (FWI) is a powerful way to invert the subsurface parameters with surface acquired seismic data. Nevertheless, the strong nonlinear relationship between the seismic data and the subsurface model will cause nonconvergence and unstable problems in practice. To divide the nonlinear inversion into some more linear steps, a 2D AI inversion imaging method is proposed to estimate the broadband AI model based on a broadband reﬂectivity. Firstly, a novel scheme based on Gaussian beam migration (GBM) is proposed to produce the point spread function (PSF) and conventional image of the subsurface. Then, the broadband reﬂectivity can be obtained by implementing deconvolution on the image with respect to the calculated PSF. Assuming that the low-wavenumber part of the AI model can be deduced by the background velocity, we implemented the AI inversion imaging scheme by merging the obtained broadband reﬂectivity as the high-wavenumber part of the AI model and produced a broadband AI result. The developed broadband migration based on GBM as the computational hotspot of the proposed 2D AI inversion imaging includes only two GBM and one Gaussian beam demigraton (Born modeling) processes. Hence, the developed broadband GBM is more efﬁcient than the broadband imaging using the least-squares migrations (LSMs) that require multiple iterations (every iteration includes one Born modeling and one migration process) to minimize the objective function of data residuals. Numerical examples of both synthetic data and ﬁeld data have demonstrated the validity and application potential of the proposed method. The bottom four subﬁgures ( e – h ) are the corresponding spectrum of top four subﬁgures in local wavenumber domain. It is clearly seen that the conventional result has a bandlimited spectrum ( e ), while the spectrum of the corrected image ( g ) has been recovered greatly, and it is very close to the true one ( h ).


Introduction
As the consumption of oil and gas resources increases, the target of seismic exploration is shifting from the large-scale structural reservoirs to small-scale stratigraphic reservoirs [1]. Acoustic impedance (AI), as a basic physical parameter, is defined as the product of density and velocity, which plays a key role to connect the surface seismic data to the subsurface model [2,3]. Therefore, the broadband AI imaging is crucial for seismic interpretation and beneficial to reservoir predication [4,5].
The typical method of AI estimation is to obtain the subsurface reflectivity first and then to calculate the AI. A recursive formula expressed by the reflectivity of normal incident is commonly used for AI estimation [6,7], followed by a prepared reflectivity from the conventional seismic migration. However, the recursive algorithm is sensitive to noise and often causes accumulated errors and inversion instabilities [8]. Furthermore, the recursive algorithm is usually applied to 1D seismic traces, but the spatial coherence of the subsurface model cannot be incorporated into the AI inversion [9,10]. To overcome the limitations of the 1D AI inversion, the full waveform inversion (FWI) [11][12][13] has been developed to estimate the broadband velocity of AI model in 2D and 3D cases. However, the strong nonliner relationship between the seismic data and the subsurface model will cause nonconvergence and unstable problems in practice.
The travel time-based velocity tomography and the conventional migrations are widely used linear inversion methods [14,15], and the full subsurface AI model can also be divided into two parts, the low-wavenumber, and the high-wavenumber parts, for two different linear inversion methods, respectively [16]. The low-wavenumber part is the background model, which can be obtained by tomography [17]. The subsurface reflectivity is related to the high-wavenumber part of AI, which can be inverted by seismic migrations [18]. The merging of the estimated low-wavenumber part and the high-wavenumber part of the AI model often leads to an intermediate wavenumber gap, and the intermediate wavenumber components are usually poorly reconstructed [14]. In recent years, the broadband, wide-azimuth, and high-density seismic data acquisition technologies have been proposed for broadband and true-amplitude seismic imaging [1,19], offering possibilities to bridge the intermediate wavenumber gap.
The least-squares migration (LSM) is a powerful method for seismic imaging of the broadband reflectivity [20][21][22]. Different LSM imaging schemes can be developed based on different seismic wave propagators. The ray-based Kirchhoff migrations method [23,24] performs efficiently but suffers accuracy problems in areas with complex subsurface. The two-way wave-based least-squares reverse time migration (LSRTM) is known as the most accurate imaging method. However, it is also known as the most expensive LSM for computation cost [25,26]. Gaussian beam is a one-way wave propagator that inherits the flexibility of ray and the accuracy of wave equation, and which has an extensive application in seismic modeling and migration [27][28][29][30][31][32]. Gaussian beam migration (GBM) can be extended for LSM by minimizing the objective function of residuals between the simulated data and the observed data [33,34]. However, the iterative fitting in data domain is a computationally intensive job, usually requiring more than 10 iterations for an acceptable result [35]. There are several noniterative approaches toward the similar goal of LSM by approximating the Hessian matrix and calculating its inverse matrix. Illumination effects linked with Hessian Matrix have been proved to be very important in retrieving the broadband reflectivity [36,37]. From the perspective of illumination analysis, the noniterative broadband reflectivity imaging can employ Kirchhoff migrations [38,39], finite difference (FD) migrations [40,41] and beam-based migrations [19]. The point spread function (PSF), as the localized elements of the full Hessian matrix, can be linked to the illumination [42,43] and produce a broadband reflectivity by a PSF-based deconvolution [44].
In this paper, we propose an effective AI imaging approach based on a broadband reflectivity imaging profile under the assumption that the background velocity or AI has been obtained by tomography. Firstly, we adopted GBM to implement migration process and calculated the PSF with the prepared background velocity. Then, we used the PSFs to deconvolve the image of GBM to produce a broadband subsurface reflectivity as the high-wavenumber part of AI model. Finally, we implemented the full AI inversion imaging by merging the obtained broadband reflectivity and the low-wavenumber part of AI model that can be derived from background velocity under assumption of constant density.

The Principle of Gaussian Beam Migration
Gaussian beam is an asymptotic solution of wave equation in a coordinate system centered on a ray under the high-frequency approximation, which combines the merits of both wave equation and ray theory. The frequency domain solution [28] of the Gaussian beam of a scalar wave equation can be expressed as where s and q T = (q 1 , q 2 ) refer to the coordinates along and perpendicular to the centeredray in the ray-centered coordinates, w is the circular frequency and s 0 is the initial point of ray. The real functions v(s) and τ(s) refer to velocity and travel time along the path of centered-ray, respectively, which can be computed by kinematic ray tracing. The complex functions P(s) and Q(s) refer to the dynamic ray parameters, which can be computed by dynamic ray tracing. For GBM, the Green's function can be decomposed into a summation of Gaussian beams [29,45] G(x, x s , ω) ≈ iω 2π dp s x dp s where u GB is an individual Gaussian beam in Cartesian coordinates, x refers to the coordinate of subsurface space, x s refers to the location of the source on the surface and p s = (p s x , p s y , p s z refers to the slowness vector in the initial emergence direction of the local plane-wave. For implementing the pre-stack GBM, we use the classic imaging condition of cross-correlation as below: where G(x, x s , ω) is the background Green's function with the source, and G * (x, x r , ω) is the conjugate of background Green's functions with the receiver. Both can be calculated using Equation (2), describing wave propagation from source and receiver locations to the image point, respectively. D(x r , x s , ω) is the recorded pre-stack data in frequency domain, and x r represents the coordinate of the receiver.

The Broadband Reflectivity Estimation Based on Point Spread Function
The pre-stack GBM image calculated using Equation (3) can be related to the subsurface reflectivity through a PSF in the local wavenumber domain [42] I(x, k) = P(x, k)R(x, k) (4) where k refers to the local wavenumber, I(x, k) refers to the imaging result of conventional migration, P(x, k) presents the PSF and R(x, k) refers to the reflectivity of subsurface model. All quantities in Equation (4) are dual domain functions, which are in local wavenumber domain defined in a small area adjacent to x. We see that the PSF serves as a filter that blurs the reflectivity and converts it into a distorted image. If the localized PSF is available, the conventional imaging result can be corrected to recover the reflectivity, as below: where R c refers to the obtained broadband image after correction. The calculation of PSF, a localized expression the full Hessian matrix P(x, k), is the key step of broadband reflectivity estimation. In this study, a novel scheme for PSF calculating the PSF by using the GBM is developed based on the method of Chen and Xie [46] and Kang et al. [47]. Note that Equation (4) is a multiplication operation in the local wavenumber domain, which is equivalent to convolution operation in the spatial domain.
where, the local coordinate x k refers to the inverse FFT of local wavenumber k. The PSF P(x, x k ) equals to the conventional imaging result I(x, x k ) if the reflectivity model is replaced by a model with scatters (every single scatter presents a delta function in mathematical aspect). From the perspective of seismic migration, the PSF P(x, x k ) can be obtained straightforward by imaging the seismic scattering wave of the scatter-point model. The Born modeling [48] is generated on the scatter-point model to simulate the original pre-stack data: where D s is the simulated data, and m s (x) is scatter-point model, which includes many scatter points (perturbation of delta functions) embedded into the background model. If the Green's functions of Gaussian beam Equation (2) are prepared, we can obtain the pre-stack data D s by the integration of Equation (7). Then, we can calculate the PSF P(x, x k ) in space domain by applying D s to Equation (3). In order to correct the image in wavenumber domain, the local FFT and the inverse FFT are generated before and after correction, and a stable factor is applied to the denominator of Equation (5) for a stabilized division.
Since the PSF includes all illumination and resolution effects from the acquisition layout, the geologic structure of subsurface and seismic wavelet, etc., and the division operation in Equation (5) tend to reverse those effects. From the prospective of LSM, the inverse of PSF can not only recover the imaging amplitude for subsurface structures but also brings a deblurring effect to the image to enhance its resolution, achieving a broadband imaging result.

The Impedance Inversion Based on the Estimated Reflectivity
The reflectivity of normal incidence is the high-wavenumber part of AI, which is important for seismic interpretation [49,50]. The 1D reflectivity of normal incidence can be expressed by AI as where ρ i+1 , v i+1 and Z i+1 = ρ i+1 v i+1 are density, velocity and AI in the (i + 1) th layer, respectively, and r i , ρ i , v i and Z i = ρ i v i are reflectivity, density, velocity and AI in the i th layer, respectively. If the absolute value of the reflection coefficients is small, this nonlinear relationship between the reflectivity and the AI can be reformulated according to [51]: We adopt a = (ln Z)/2, and the liner relationship between 1D reflectivity and a can be derived by where, r refers to the 1D reflection coefficients series of all reflectors and a refers to the corresponding series of half natural log of AI. ∇ z refers to vertical one-order difference operator. By considering the spatial coherence between traces of the AI model, we incorporate the 2D seismic profiles and employ the forward model that solves for multitraces simultaneously. Equation (10) can be extended to the 2D case according to [9]: where R = (r 1 , r 2 , . . . , r M ) refers to a matrix that denotes the reflectivity imaging profile of normal incidence, the matrix A = (a 1 , a 2 , . . . , a M ) denotes the corresponding profile of half natural log of AI and ∇ n refers to vertical one-order difference operator in the normal direction of the reflector. The TV regularization is introduced to preserve the blocky characteristics of structure, and the AI profile can be estimated as follows [9]: Energies 2021, 14, 4105

of 12
Here, A low is the given background of A, and L refers to the low-pass filter. A low = LA is a boundary constraint condition in this AI inversion scheme. In this scheme, the broadband reflectivity takes the role of seaming the wavenumber gap between the lowwavenumber part of AI and the high-wavenumber part. In the case of 2D model, the TV norm in Equation (12) can be written as where ∇ x is the horizontal one-order difference operator. The constrained inversion problem Equation (12) cannot be solved directly, as the TV regularization term is nondifferentiable. The Bregman iteration inversion method [52] and the Shrinkage-Thresholding method [53] can be introduced to solve this problem.

Numerical Examples
The Sigsbee2a model was designed according to the geological background near the Sigsbee escarpment, containing several normal and thrust faults that separate the sedimentary deposits [54]. The left part of the Sigsbee2a model was chosen to test the proposed method, which is parameterized by 361 × 601 mesh points in the 2D model and 8 m sampling intervals in both horizontal and vertical directions. The true impedance and reflectivity model are shown in Figure 1a,b, respectively. The FD-based Born modeling was used for the original pre-stack data, and a Ricker wavelet with the main spectrum in 25 Hz was selected for the modeling. The pre-stack data set contains 150 shots that are evenly distributed on the surface with a 32 m sampling interval. Each common shot point gathered has 201 receivers arrayed with offset from −1600 m to 1600 m at a 16 m spacing. A sorted common shot gather that is located at 2880 m, and the zero offset gather are shown in Figure 1c,d, respectively. we also zoomed out one PSF to show its spectrum in detail, as the red rectangle shows in Figure 2b. We found that the PSFs in shallower positions have superior amplitude than the deeper ones, and the PSFs in the middle have broader wavenumber coverage than the sider ones, which reveals that the shallower and middle positions have better seismic illumination and coverage.   To produce the PSFs for correction, we inserted 493 point scatterers in the background (migration) velocity with a 168 m sampling interval in both horizontal and vertical directions. The pre-stack diffracted data can be generated by using the Born modeling method on the designed point scatterers model, as calculated by Equation (7). Based on the obtained diffracted data in Equation (7), we implemented the conventional GBM to produce all subsurface PSFs at once. The PSFs in space domain are shown in Figure 2a, which contain the main blurring effects of the conventional migration as a localized Hessian matrix of the inversion imaging system. One PSF in Figure 2a is zoomed in to show its shape in details, but it still has not directly revealed seismic illumination. Therefore, we implemented the FFT to obtain the spectrum of the PSFs in local wavenumber domain, which are illustrated in Figure 2b. In addition to illuminating all PSFs in the whole model, we also zoomed out one PSF to show its spectrum in detail, as the red rectangle shows in Figure 2b. We found that the PSFs in shallower positions have superior amplitude than the deeper ones, and the PSFs in the middle have broader wavenumber coverage than the sider ones, which reveals that the shallower and middle positions have better seismic illumination and coverage.  Next, we implemented conventional GBM to produce the subsurface image and then performed the correction to produce a broadband reflectivity using Equation (5). The image from the conventional GBM and the corrected image using the proposed method are shown in Figure 3a,b, respectively. By comparing conventional GBM image with the corrected reflectivity image, we found that the image after correction has more balanced amplitudes in both the shallower and the deeper parts of the model, as well as more abundant wavenumber contents. We chose a local part of the model to show the quantitative im- Next, we implemented conventional GBM to produce the subsurface image and then performed the correction to produce a broadband reflectivity using Equation (5). The image from the conventional GBM and the corrected image using the proposed method are shown in Figure 3a,b, respectively. By comparing conventional GBM image with the corrected reflectivity image, we found that the image after correction has more balanced amplitudes in both the shallower and the deeper parts of the model, as well as more abundant wavenumber contents. We chose a local part of the model to show the quantitative improvement of the PSF based correction. Figure 4a  The results indicate that the spectrum of the corrected image has been recovered greatly, and it is very close to the true spectrum.
Owing to the broadband reflectivity that was obtained beforehand, the absolute AI imaging by using the boundary constraint inversion scheme can be implemented as Equation (12). The constant density of 2500 kg/m 3 is used in this inversion scheme as Equation (12), and the low wavenumber part of AI is deduced by the background migration velocity. As shown in Figure 5a,b for comparison, we implemented this scheme by using conventional GBM image and the corrected image with PSF, respectively. By comparing the absolute AI results from the corrected image and the conventional GBM image, the AI inversed by the corrected image shows better consistency with the true AI model of Figure 1a. For further comparison between the inversed AI results from conventional GBM and the proposed method, we implemented the FFT to local areas of both results. Figure 5c,d show the spectrum of the areas marked by red rectangles; Figure 5e shows the spectrum of the corresponding true AI model in Figure 1a. We can see that the inverted impedance of the proposed method has a wider spectral band than the conventional one, and it has similar spectral pattern to the true spectrum. For quantitative comparison of the results in space domain, we measured the relative errors of the inverted AI profiles from both conventional GBM and the proposed broadband GBM. The comparison results in Figure 6a,b show that the relative errors of the proposed method are almost all below 5%, and the relative errors of traditional method are greater than 10% in many places. For further comparison of absolute value of the inverted AI, we extract two traces at 1840 m and 3280 m from Figure 5a,b, and the extracted traces are indicated in Figure 7a,b, respectively. The green, blue and red dashed lines represent the inverted impedance results from the conventional GBM, the proposed method and true reflectivity, respectively. It is obvious that the AI inversed by the proposed method is very close to the true AI, and the AI inversed by conventional GBM is far away from the true AI.   It is clearly seen that the conventional result has a bandlimited spectrum (e), while the spectrum of the corrected image (g) has been recovered greatly, and it is very close to the true one (h).
Owing to the broadband reflectivity that was obtained beforehand, the absolute AI imaging by using the boundary constraint inversion scheme can be implemented as Equation (12). The constant density of 2500 kg/m 3 is used in this inversion scheme as Equation (12), and the low wavenumber part of AI is deduced by the background migration velocity. As shown in Figure 5a,b for comparison, we implemented this scheme by using conventional GBM image and the corrected image with PSF, respectively. By comparing the absolute AI results from the corrected image and the conventional GBM image, the AI inversed by the corrected image shows better consistency with the true AI model of Figure  1a. For further comparison between the inversed AI results from conventional GBM and the proposed method, we implemented the FFT to local areas of both results. Figure 5c,d   It is clearly seen that the conventional result has a bandlimited spectrum (e), while the spectrum of the corrected image (g) has been recovered greatly, and it is very close to the true one (h).
Owing to the broadband reflectivity that was obtained beforehand, the absolute AI imaging by using the boundary constraint inversion scheme can be implemented as Equation (12). The constant density of 2500 kg/m 3 is used in this inversion scheme as Equation (12), and the low wavenumber part of AI is deduced by the background migration velocity. As shown in Figure 5a,b for comparison, we implemented this scheme by using conventional GBM image and the corrected image with PSF, respectively. By comparing the absolute AI results from the corrected image and the conventional GBM image, the AI in-  relative errors of traditional method are greater than 10% in many places. For further comparison of absolute value of the inverted AI, we extract two traces at 1840 m and 3280 m from Figure 5a,b, and the extracted traces are indicated in Figure 7a,b, respectively. The green, blue and red dashed lines represent the inverted impedance results from the conventional GBM, the proposed method and true reflectivity, respectively. It is obvious that the AI inversed by the proposed method is very close to the true AI, and the AI inversed by conventional GBM is far away from the true AI.  (a,b), respectively, and (e) shows the spectrum of the same position of the true AI model in Figure 1a. We can see that the inverted impedance of the proposed method has a wider spectral band than the conventional one, and it has similar spectral pattern to the true one.  (a,b), respectively, and (e) shows the spectrum of the same position of the true AI model in Figure 1a. We can see that the inverted impedance of the proposed method has a wider spectral band than the conventional one, and it has similar spectral pattern to the true one.
show that the relative errors of the proposed method are almost all below 5%, and the relative errors of traditional method are greater than 10% in many places. For further comparison of absolute value of the inverted AI, we extract two traces at 1840 m and 3280 m from Figure 5a,b, and the extracted traces are indicated in Figure 7a,b, respectively. The green, blue and red dashed lines represent the inverted impedance results from the conventional GBM, the proposed method and true reflectivity, respectively. It is obvious that the AI inversed by the proposed method is very close to the true AI, and the AI inversed by conventional GBM is far away from the true AI.  Figure 1a. We can see that the inverted impedance of the proposed method has a wider spectral band than the conventional one, and it has similar spectral pattern to the true one. To further certify the application potential of this method, we also chose a field seismic data set for testing. The whole field data came from China, and we sorted 91 shots for testing. All sorted shots are located on the surface from 0 to 5500 m with a 60 m interval, and each shot has 75 receivers bilaterally located from offset 40 m to 3000 m at a spacing of 40 m. Figure 8a shows the background AI model from the industry, which is obtained by velocity analysis and tomography module of commercial software. Figure 8b shows the broadband reflectivity of the field data deduced by the proposed GBM. Figure 8c is the inverted AI through merging the obtained broadband reflectivity and background AI model. It can be found in Figure 8c that the detailed structural features of the inverted AI are consistent with the reflectivity profile of Figure 8b, and the macro features are consistent with the background AI of Figure 8a. This test demonstrates that the proposed AI inversion imaging scheme has the capability to merge the obtained broadband reflectivity and the given background AI and produce the absolute broadband AI model effectively. To further certify the application potential of this method, we also chose a field seismic data set for testing. The whole field data came from China, and we sorted 91 shots for testing. All sorted shots are located on the surface from 0 to 5500 m with a 60 m interval, and each shot has 75 receivers bilaterally located from offset 40 m to 3000 m at a spacing of 40 m. Figure 8a shows the background AI model from the industry, which is obtained by velocity analysis and tomography module of commercial software. Figure 8b shows the broadband reflectivity of the field data deduced by the proposed GBM. Figure 8c is the inverted AI through merging the obtained broadband reflectivity and background AI model. It can be found in Figure 8c that the detailed structural features of the inverted AI are consistent with the reflectivity profile of Figure 8b, and the macro features are consistent with the background AI of Figure 8a. This test demonstrates that the proposed AI inversion imaging scheme has the capability to merge the obtained broadband reflectivity and the given background AI and produce the absolute broadband AI model effectively.

Discussions
Today, the wide-azimuth and high-fidelity seismic datasets enable the high dimensional inversion imaging to estimate the subsurface parameters with higher resolution. The classic FWI is a popular high dimensional nonlinear inversion method from seismic data to subsurface model but usually suffers the "cycle-skipping" problem [55]. In this study, we focused on the broadband reflectivity under the linear LSM inversion frame and then estimated the broadband AI model by incorporating the background AI model. In other words, we divided the full model into different parts and tried to invert different parts by the linear inversion separately.
We developed the GBM to implement the LSM, because the GBM has better imple- Figure 8. The background AI model of a field data (a), the broadband reflectivity generated from the proposed imaging method with PSF correction (b) and the inverted absolute AI result of the proposed AI inversion (c).

Discussions
Today, the wide-azimuth and high-fidelity seismic datasets enable the high dimensional inversion imaging to estimate the subsurface parameters with higher resolution. The classic FWI is a popular high dimensional nonlinear inversion method from seismic data to subsurface model but usually suffers the "cycle-skipping" problem [55]. In this study, we focused on the broadband reflectivity under the linear LSM inversion frame and then estimated the broadband AI model by incorporating the background AI model. In other words, we divided the full model into different parts and tried to invert different parts by the linear inversion separately.
We developed the GBM to implement the LSM, because the GBM has better implementation efficiency and superiority for broadband imaging than the FD-based migrations. Specifically, the computational cost increased nearly 1600% when FD-based methods are used to propagate seismic wave with frequency from 30 Hz to 70 Hz, while beam-based method only requires less than a 90% increment in the same case [19]. Moreover, the GBM is a one-way propagator that does not have a back-scattering wave in propagation, and the GBM-based method has no problem of crosstalk noises in PSF calculation due to interferometry of scatter-points [46,55].
In this study, the low-wavenumber part of the model is assumed to be known beforehand, because this part is commonly inverted by travel time inversion [56], which is linear inversion and demonstrates higher stability than FWI. We used the constant density model in all the tests in this paper, because density is not sensitive to travel time, and it always couples with velocity in AI inversion. Although the acoustic wave equation with constant density is commonly used in the industry for impedance inversion as the geometric position of subsurface reflectors can be inverted without density [10,16,57], we are still in the transition from AI inversion to the elastic impedance inversion. Since well-logging is an effective way to get the density directly, the elastic impedance inversion can be developed to adapt spatial-variable density model in our future work.

Conclusions
As an alternative of the nonlinear inversion imaging, we propose an effective AI inversion imaging method based on a broadband GBM imaging. The extended GBM method is a linear imaging under LSM framework, which produces the broadband reflectivity as the high-wavenumber part of the AI model. The numerical tests demonstrate that the bandwidth of the reflectivity is important for inverting the detailed texture of the AI model, and the background AI should be included to guarantee the correct value of the absolute AI model.
The proposed broadband imaging method includes only one Gaussian beam Born modeling and two GBM processes, which is more efficient than the conventional LSM that requires multiple iterations of Born modeling and gradient calculation (migration). Benefitting from the implementation flexibility of beam-based method with complex geometry, this new method can be developed to process the AI inversion in the case of rugged topography straightforward. In addition, the developed method can be extended to 3D cases with less computational cost by using GBM than FD-based migrations.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.