Elastic reverse time migration based on nested triangular mesh in 2D isotropic media (2024)

Search for other works by this author on:

Article history

Abstract

1. Introduction

Reverse time migration (RTM) (Baysal 1983, McMechan 1983, Whitmore 1983) extrapolates wavefield using the full wave equation. This technique is advantageous due to its absence of dip limitations and its capability to image complicated geological structures distinguishing itself from other migration algorithms. This method is a pivotal technique in the field of seismic imaging. Recently, elastic RTM (ERTM) has become increasingly feasible due to advancements in seismic acquisition and computational capabilities. Moreover, as the Earth is an inhom*ogeneous elastic medium, ERTM, in comparison to single-component RTM, can provide more accurate images with a clear physical interpretation.

To date, two approaches have been proposed for ERTM. One method utilizes Helmholtz decomposition (Sun etal. 2004, Yan and Sava 2008, Du etal. 2012, Duan and Sava 2015, Liu 2019; Lu etal. 2020) to decompose the multi-component wavefields into P- and S-waves during wavefield extrapolation. However, this technique changes the amplitude and phase of the wavefield. Moreover, when dealing with converted PS and SP images, correction for polarity reversals is required. To address this issue, an alternative ERTM scheme by extrapolating the decoupled elastic wave equation (Xiao and Leaney 2010, Wang and McMechan 2015, Wang etal. 2016, Yong etal. 2016, Du etal. 2017, Yang etal. 2018, Shi etal. 2019, Zhang and Shi 2019, Zhou etal. 2019, Zhang etal. 2020, Bao etal. 2022, Pan etal. 2022) has been proposed. In this approach, the P- and S-wave vector decomposition preserves amplitude and phase properties, eliminating the need for correcting polarity reversals in PS and SP images. In this study, we adopt the latter strategy.

Surface topography brings challenge to the propagation of elastic waves, for example, the free surface conditions, the conversion of different waves, and the scattered waves in the near surface. Therefore, an appropriate forward modeling algorithm is essential for ERTM. In recent decades, various wave equation methods have been widely utilized to tackle surface topography challenges. These methods can be broadly categorized into several groups, including finite difference method, finite-element method, spectral element method, and boundary element method, etc. Each of these methods have its advantages and disadvantages. In the numerical simulation of surface, the appropriate method is often selected according to the characteristics of the simulated underground structure. In the forward modeling of an undulating surface, mesh discretization is employed, and meshes are generally categorized into two types: structured and unstructured. The so-called structural mesh is simply one in which all the internal points in the grid area have the same number of adjacent units. The generation techniques of structural mesh include conventional rectangular grid, curvilinear grid (Zhang and Chen 2006; Zang etal. 2021), body-fitted coordinate method, etc. By contrast, an unstructured grid (Zhang and Liu 1999) means that internal points within the grid region do not share the same number of adjacent units.

In finite-element calculations, due to the non-diagonal nature of the mass matrix, the process of matrix inversion requires significant computational and storage resources. However, the grid method transforms the mass matrix into a diagonal matrix, significantly reducing storage requirements. Moreover, it simplifies the process of matrix inversion, maintaining the advantages of finite-element methods in accurately depicting the complex surface and adaptive mesh refinement. It also facilitates the handling of free surface conditions. Moreover, the grid method combines the advantages of finite differences by employing Gauss's divergence theorem, converting the solution of surface integrals into line integral solutions. In the computation of stress for triangular elements, this method efficiently utilizes linear interpolation algorithms for triangles, resulting in high computational efficiency.

In this study, we used the grid method to simulate elastic wave equation in the region with rugged topography. For the workflow of ERTM, there are two main steps: forward-propagating the source wavefield in a migration velocity model, and back-propagating the recorded wavefield through the same velocity mode. For the grid method, if we all use small meshes to simulate wavefield propagation for the two steps, the computational cost is high. To address this problem, we propose the following strategies: for the forward and backward simulations of the source wavefields, the high-order large meshes are used; for the backward-propagated receiver wavefield, the nested low-order small and high-order large meshes are used.

The organization of this paper is as follows. First, the methodology of grid method is presented. Second, the implementation of ERTM based on the nested unstructured meshes is demonstrated. Finally, numerical tests on three synthetic models demonstrate the validity of the proposed method.

2. Methods

The grid method (Zhang and Liu 1999) is a numerical technique founded on an irregular and unstructured mesh, providing the capability to accurately model acoustic or elastic wave propagation in complicated structures, featuring the interface and surface topographies. This approach is capable of simulating wave propagation in regions with strong variations in velocities. In this section, we will provide a concise overview of the grid method.

2.1. The decouple elastic wave equation based on the grid method

The 2D first-order velocity-stress elastic equations can be expressed as

where |${v_x}$| and |${v_z}$| are the horizontal and vertical components of particle-velocity, |${\tau _{xx}}$| and |${\tau _{zz}}$| are the horizontal and vertical components of stress tensor, |${\tau _{xz}}$| is the shear stress, ρ is the density, t is the time, and λ and μ are the Lamé parameters. To solve Equation(1), we integrate both sides of this equation around each spatial point (taking the node k shown in Fig.1 as an example), and obtain the following expressions by discretizing (Zhang and Liu 1999):

where m denotes the number of triangular elements around the target point k, M_k denotes the parameter about the total mass sum of all triangles around point k, and b_k and c_k denote coefficients. Considering the triangle ijk shown in Fig.1 as an example, the first-order spatial derivatives |$\frac{{\partial {v_x}}}{{\partial x}}$| and |$\frac{{\partial {v_x}}}{{\partial z}}$| can be derived as follows:

where Δ denotes the area of triangle ijk, |${b_i} = \frac{{{z_k} - {z_j}}}{2}$|⁠, |${c_i} = \frac{{{x_k} - {x_j}}}{2}$| and the remaining coefficients can be obtained similarly using the values of coordinates for points i, j, and k in the same way. The calculation process of components |$\frac{{\partial {v_z}}}{{\partial x}}$| and |$\frac{{\partial {v_z}}}{{\partial z}}$| is similar to Equation(3). Equations(2) and (3) provide a means to solve Equation(1) using the grid method.

As the implementation of ERTM in this study is based on the decoupled equations, we need to make some changes based on the original Equation(1) (Lu and Liu 2018):

where |${v_p}$| and |${v_s}$| denote the P- and S-wave velocities, respectively, |${\bf{v}} = ({v_x},{v_z})$|⁠, |${{\bf{v}}^{\bf{P}}} = ({v_{px}},{v_{pz}})$|⁠, and |${{\bf{v}}^{\bf{S}}} = ({v_{sx}},{v_{sz}})$|⁠. We summarize the whole wavefield extrapolation process as follows: first, compute the velocity components |${\bf{v}} = ({v_x},{v_z})$| and the stress tensor components |$\tau = ({\tau _{xx}},{\tau _{zz}},{\tau _{xz}})$| based the grid method using Equation(1). Second, compute the scalar wavefield θ and vector wavefield ω using Equation (4). Finally, compute the vector P- and S-wavefields |${{\bf{v}}^{\bf{P}}}$| and |${{\bf{v}}^{\bf{S}}}$|⁠.

2.2. Mesh generation and wavefield extrapolation

Here, we will explain the mesh generation algorithm and the nested mesh framework proposed in this study. The geometric information of undulating terrain or velocity interfaces is taken as constraints. Subsequently, using the Delaunay algorithm with locally adaptive discretization intervals, a mesh is generated. For a given Voronoi polygon, its centroid can be defined as

The reason for introducing Voronoi polygons first is that they are duals of Delaunay triangles. In the mesh generation process, the density function |$\rho (x)$| is defined as

where |${v_{\max }}$| is the maximal migration velocity, and |$\tilde v(x,z)$| is the local migration velocity after smoothing determines the spatial sampling step size. Unfortunately, the Delaunay algorithm does not guarantee the quality of the generated mesh. We will employ the centroid Voronoi tesselations (CVTs) algorithm, taking the mesh generated by Delaunay as the initial solution. If the centroids of a set of Voronoi polygons coincide with the vertices of the corresponding Delaunay triangles, then we can refer to these as CVTs. Local velocity is used as the density function to calculate the global moment of inertia. The objective function is iteratively optimized to generate a higher-quality mesh, which is written as

To optimize this objective function, the Lloyd algorithm is employed, which is an explicit iterative algorithm for constructing CVT. This algorithm continuously moves the centroids of Voronoi polygons, making the objective function F monotonically decrease until certain convergence criteria are satisfied. In this way, velocity-adaptive, higher-quality grids with closer-to-equilateral triangles can be generated.

Low-order equal-divided meshes effectively capture the fluctuations of complex interfaces due to their small scale. Conversely, high-order equalized meshes (as shown in Fig.2) offer a higher acceleration ratio, yet their larger scale presents challenges in accurately describing detailed surface features. For ERTM, to simulate the forward and backward propagation of seismic waves in complex topography media, we adopt the following strategies:

First, based on the surface elevation changes, an intersection line is automatically fitted within the shallow layer, as shown in Fig.3, representing the overall trend of elevation variation. For the receiver wavefield, first-order equal division is applied to the upper section of the boundary line to describe complex surface undulations, as shown in Fig.4.

Second, in the lower part of the intersection line and all part of the shot wavefield, the 16 high-order equal divisions are adopted to take advantage of the high acceleration ratio. The nodes of low- and high-order meshes overlap on the intersection line so the nodes on the intersection line can exchange information explicitly to update the wavefield.

CUDA has many successes applications in the field of exploration, and the computing efficiency is greatly improved. However, compared to the conventional finite difference methods based on structured grids, using grid method based on unstructured meshes directly accelerated by graphics processing units (GPU) is relatively difficult, largely due to the scattered data read–write pattern inherent. As a result, the effectiveness of GPU acceleration for numerical algorithms based on a double mesh is proposed (Yang etal. 2017), which effectively addresses the discontinuity issues in GPU memory access. A double mesh consists of a large mesh and small meshes within it, as illustrated in Fig.2. The large one is partitioned based on dispersion requirements and constraints imposed by the velocity interface. On the other hand, the small ones are automatically generated through structured 16-order equal division. Exploiting the similarity and regular arrangement of these small grids allows for reduced memory storage, as well as the consolidation of necessary storage reads, enhancing the GPU's acceleration performance.

2.3. Imaging condition

The implementation process of ERTM mainly consists of three steps: forward extrapolation of the source wavefield, reverse extrapolation of the receiver wavefield, and applying imaging conditions. The cross-correlation imaging condition for ERTM is as follows (Liu 2019):

where |${I_{PP}}(x,z)$| and |${I_{PS}}(x,z)$| are the PP and PS migration stack profiles, respectively, |${S_P}(x,z,t)$| is the seismic source wavefield, and |${R_P}(x,z,t)$| and |${R_S}(x,z,t)$| are the P- and S-wave receiver wavefields, respectively. The benefit of RTM technology is its capability to image the entire wavefield without limitations on imaging angles. Therefore, large-angle low-frequency noise can also be included in the migration results. In the implementation of imaging conditions, the angles |${D_{PP}}({{\hat{\bf n}}_s},{{\hat{\bf n}}_g},\alpha )$| and |${D_{PS}}({{\hat{\bf n}}_s},{{\hat{\bf n}}_d},\alpha )$| calculated by Poynting vector are restricted to eliminate low-frequency noise.

2.4. Numerical implementation

The proposed scheme that handles ERTM based on the grid method can be performed as follows:

• Input the model parameters and generate the shot and receiver meshes, then store the mesh partitioning parameters.

• Forward the source wavefield using the grid method based on the large mesh.

• Backward propagate the receiver wavefield with the observed data on a nested mesh, and cross-correlate the source and receiver wavefields on the large mesh.

• Repeat 2–3 times for all shots, and stack the migration results from multiple shots to generate the final migration results.

3. Numerical experiments

Here, we use three synthetic datasets to test the effectiveness of our ERTM procedure, and the forward simulation data are generated using the grid method. We use the perfectly matched layer boundary condition for forward modeling and RTM, and the receivers are located within the computational domain. Through applying the proposed procedure, the efficiency of the original small mesh partitioning has improved from several minutes to several tens of seconds. The forward simulation and imaging processes have also been accelerated from the hourly level to the minute level.

3.1. Three-layer model

First, we use a three-layer model with topography to test the proposed method, as shown in Fig.5. This model consists of 1447 points in the horizontal direction, and 801 points in the vertical direction spaced at intervals of 10m. There are 140 explosive sources placed at the model surface with an interval of 100 m, and the synthetic seismic recorded data for each shot were generated with a total of 723 receivers positioned on the surface, uniformly spaced at 20-m intervals. The source function employed in this model was a Ricker wavelet with a peak frequency of 20 Hz. The time sampling interval is 5 ms, and the total recording time is 4 s. Figure6 illustrates a synthetic seismic record generated by the grid method with a shot positioned at the distance of 7.5km. Figure7 parts a and b display the PP and PS migration profiles, respectively. Note that our method provides clear image results of this simple topography model.

Table1 illustrates the generation time costs of the nested mesh proposed in this study compared to the traditional low-order mesh generation. It is evident from the table that the efficiency of our approach is significantly higher. Consequently, the efficiency of wavefield simulation corresponding to this approach is expected to improve accordingly. Wavefield simulations on low-order meshes cannot be accelerated using GPU, whereas our nested mesh approach can be accelerated using GPU.

3.2. Hess model

Second, we employ the modified Hess model with irregular surface to further test the validity of our method. The velocity model parameters are shown in Fig.8, which include high-velocity salt structure, complex fault, and undulating surfaces, and we set the density to a fixed constant number. The model's computational area is 36.16 km long and 11 km wide, with spatial step size of 10 m in both the x- and z-directions. The source is loaded as a Riker wavelet with a peak frequency of 20 Hz, and 360 shots are placed on the model surface to generate the seismic data with intervals of 100 m. At the same time, there are 2000 receivers with intervals of 15 m. The time sampling interval is 5 ms, and the total recording time is 8 s. Figure9 presents a synthetic seismic record generated by the grid method with the nested triangular mesh, and the shot is positioned at a distance of 17 km. Figure10a and b present the PP and PS migration results of the modified Hess model, which clearly show the complex high-velocity salt structure, fault, and rugged surface features. Correspondingly, Fig.11 shows the imaging results without considering the influence of surface topography. The comparison of the two figures shows that the imaging structures of our method is consistent with that of the flat terrain imaging structures, indicating the effectiveness of the proposed method. Similarly, Table2 demonstrates the superiority of the proposed approach.

3.3. Foothill model

Finally, we evaluate the effectiveness of our approach using the Foothill model, which features intricate surfaces and strong velocity variations. The P-wave velocity of the Foothill is depicted in Fig.12, the S-wave velocity is equal to V_p/2, and we set the density to a constant number. The lateral extent of the Foothill model is 25km, and the depth range of the model is 10km. The maximum elevation difference is 1.4km, and the grid spacings are set to 15 and 10m in the x- and z-directions, respectively. The seismic recording time is 8s with an interval of 0.5ms. A Riker wavelet with a peak frequency of 20Hz is loaded to the source. The synthetic records consist of 278 common-shot gathers with a shot spacing of 90m, received by 833 receivers, and the offset ranges from 15 to 3600m. Figure13 illustrates the migrated stacked profiles for the Foothill model, we can see that the images of the main structures are very clear and continuous, the rugged surface is accurately represented and the imaging depth is accurate. This numerical demonstration highlights the capability of ERTM based on the grid method to produce high-resolution imaging results.

4. Conclusions

ERTM based on the unstructured meshes is a crucial tool in imaging the underground structures with complex rugged topography. However, the computational complexity and efficiency become the important factors that restrict the applications of this method. To depict rugged topography, very fine meshes are needed for generation, which leads to low computational efficiency. In this study, we proposed to adopt two kinds of different mesh combined with GPU acceleration to implement ERTM. Several surface layers are discretized by small fine meshes to describe topography, and remaining layers discretized by large coarse meshes such as the flat surface model. During wavefield propagation, for the forward and backward propagation of a source wavefield, large meshes are used, and for the backward propagation of a receiver wavefield, nested small and large meshes are used. In this way, the complexity of mesh generation and the computational cost of wavefield propagation are greatly reduced. Through numerical examples, we demonstrated the validity and effectiveness of the proposed method.

Conflict of interest statement: None declared.

Funding

This research was supported by the University Stability Support Program of Shenzhen-General Project (grant no. 20231128113233002), the National Science Foundation of China (grant nos. 42104119 and 42104106) and Mathematical Model and Machine Learning Algorithm of Seismic Surface Waves (grant no. 2021KTSCX145).

Data availability

Data associated with this research are available and can be obtained by contacting the corresponding author.

References

© The Author(s) 2024. Published by Oxford University Press on behalf of the SINOPEC Geophysical Research Institute Co., Ltd.

Issue Section:

Research article

Email alerts

Citing articles via

More from Oxford Academic