A three-dimensional immersed boundary method for accurate simulation of acoustic wavefields with complex surface topography

Abstract

1. Introduction

The impact of irregular topography is crucial in seismic processing, particularly in the areas of migration imaging and full-waveform inversion (Tarantola 1984, Beasley and Lynn 1989, Reshef 1991, Li et al. 2016, Borisov et al. 2018, Li et al. 2022). In 3D numerical simulations of seismic waves incorporating real topography, the Earth model approximates the irregular surface as a staircase structure when discretized using cuboidal grids. This approximation leads to artefacts and traveltime errors in events related to the surface (Lombard and Piraux 2004). These errors can affect seismic imaging and full-waveform inversion. Therefore, it is necessary to develop a stable and accurate wavefield simulation method that can handle true 3D surface topography in land seismic exploration.

Wavefield simulation employs various methods, including finite-difference methods (FDM) as discussed by Alterman and Karal (1968), Virieux (1986), Levander (1988), Yao et al. (2016, 2018a, 2018b), Mulder and Huiskes (2017), Cao and Chen (2018), and Li et al. (2020); pseudo-spectral methods (PSM) as described by Kosloff and Baysal (1982); finite-element methods (FEM) detailed by Marfurt (1984), Liu et al. (2014), and Meng and Fu (2017); and spectral-element methods (SEM), as explored by Kelly et al. (1976), Komatitsch and Vilotte (1998), and Oral et al. (2019). Each has its advantages; for instance, FEM and SEM adeptly handle the irregular surface by using triangular/tetrahedral and quadrilateral/hexahedral mesh designs and are beneficial for simulating low-velocity layers. However, they are more computation-intensive than FDM and PSM (Alterman and Karal 1968, Kosloff and Baysal 1982, Marfurt 1984, Virieux 1986).

FDM is widely used in computational seismology due to its simplicity, and high computational efficiency. Therefore, it is particularly favoured for implementing reverse-time migration (RTM) and full-waveform inversion. However, the use of conventional rectangular/cuboid grids in FDM can result in artefacts such as seismic scattering and diffractions when dealing with an undulated topography. These artefacts can degrade the quality of migration results, posing challenges for accurate subsurface imaging. The artefacts can be mitigated by using finer grids (Hayashi et al. 2001, Bohlen and Saenger 2006). Numerical experiments have demonstrated that at least 60 sampling points per minimum wavelength are required to accurately simulate the seismic wavefield with an irregular surface, resulting in a significant increase in computation (Hayashi et al. 2001, Bohlen and Saenger 2006, Bleibinhaus and Rondenay 2009). By contrast, Tarrass et al. (2011) applied optimized filters to a conventional grid to address the irregular topography issue.

To address the limitations of traditional FDM, various schemes with curvilinear grids and coordinate systems have been explored. Curvilinear grid-based FDM aligns the discretized grid more closely with the irregular surface, reducing ‘staircase’ artefacts (Puente et al. 2014). These methods introduce additional partial derivative terms in the wave equations, increasing computational cost and complexity. In addition, a skewed grid is prone to becoming unstable (Fornberg 1988, Jastram and Tessmer 1994, Sun et al. 2016, 2017, Li et al. 2019, Wang et al. 2019, Li and Qu 2022, Liu 2023). Rao and Wang (2013, 2018) employed the summation-by-parts finite-difference method to tackle instabilities in forward modelling in the case of a significantly skewed grid. Qu et al. (2022a, 2022b) proposed a new method for solving continuous wavefield operators in a 3D curvilinear domain, improving accuracy and efficiency in wavefield modelling with topography and achieving good results in RTM.

To achieve free surface boundary conditions with curved grids, Zhang et al. (2012) proposed the traction image method with velocity–stress wave equations. Appelö and Petersson (2007) presented a boundary-modified difference operator for second-order displacement–stress wave equations.

In Cartesian FDM, numerical scattering occurs due to discretizing the topographic free surface using a staircase-like approach. Approaches such as the embedded-boundary method have been used to address this issue in seismic wave modelling (Tseng and Ferziger 2003, Zhao 2010, Gao et al. 2015), but stability can become a concern after a long simulation time (Kreiss and Petersson 2006, Lombard et al. 2008, Mittal et al. 2008). Methods such as the arbitrary accuracy derivatives approach have been used to manage curved surfaces with staggered-grid finite difference (FD), as discussed by Almuhaidib and Toksöz (2015). Hu (2016) used interpolation techniques to calculate the wavefield at ghost points, which are directly aligned and orthogonal to the surface, followed by iterative extrapolation from these inside fictitious points to ghost points. However, both methodologies entail extrapolation within the forward modelling framework, which may induce instability and reduce the accuracy of the FD wavefield modelling.

To overcome these obstacles, Li et al. (2020) introduced a refined immersed boundary method, which uses iterative symmetric interpolation for the acoustic wavefield simulation in the presence of irregular topography. Initially developed for two-dimensional applications, this methodology is expanded to incorporate a three-dimensional immersed boundary method (3DIBM) in this paper. It computes spatial coordinates for ghost and mirror points in 3D space, employing 3D spatial symmetric interpolation iteratively. This notably includes the Kaiser-windowed sinc interpolation technique for wavefield interpolation. This approach not only reduces diffraction waves but also enhances the accuracy of traveltimes for free surface reflections, offering an improvement over traditional methods such as the vacuum method.

2. Method

2.1. Forward modelling with 3DIBM

The acoustic wave equation can be expressed as

where |$p({{\bf x}},t)$| signifies the pressure component, V represents velocity, |$\rho $| denotes density, |$s(t)$| designates the source signature, and |$\delta ({{\bf x}} - {{{{\bf x}}}_s})$| symbolizes the Dirichlet delta function, which is employed for injecting the source at the specified location |${{{{\bf x}}}_s}$|⁠.

The Earth's free surface |$\Omega $| is a rigid reflective boundary in geophysics, vital for accurate seismic wavefield simulations. The free surface boundary condition can be expressed as

where |${{\bf n}}$| represents the surface's normal direction at the free surface |$\Omega $|⁠. The free surface boundary condition in Equation (2) implies that the traction is zero in the surface's normal direction. In acoustic wave simulations, the pressure component p is equal to the sum of normal stresses and equal to zero at the surface, which can be expressed as

The second-order wave equation, denoted as Equation (1), can be resolved by employing the FD using regular mesh grids. The free surface, demarcating the boundary between the air and subsurface media, exhibits a reflection coefficient of −1. For a horizontal surface, the free surface condition for a flat surface is satisfied by directly antisymmetric mirroring of the wavefields above and below the surface, maintaining zero pressure at the boundary.

However, challenges emerge with undulating terrains in regular mesh-based FD simulation, as the irregular surfaces may not align with grid nodes, resulting in compromised accuracy. To address this issue, the immersed boundary method (IBM) is introduced as a pivotal solution for implementing the free boundary condition of an undulating surface in seismic wavefield simulations. Figure 1 illustrates how the 3DIBM enables the seismic simulation with a sloped surface in 3D space. The 3D surface is immersed into the entire 3D meshes. Points located on integer nodes above the free surface are integral to FDM computation and are designated as ghost points. These ghost points, situated within the air area above the surface, and their corresponding mirror points inside the model exhibit symmetry relative to the 3D surface.

Figure 1.

Schematic illustration of the IBM applied in a 3D spatial domain. The coloured surface represents a 3D sloped topography. An unfilled circle symbolizes a ghost point situated on an integer node above the 3D surface within the air region, whereas the filled black circle marks the corresponding mirror point situated on a fractional node inside the model. A black line orthogonal to the 3D free surface intersects both the ghost and mirror points pair. Additionally, the grey dots denote integer grid points engaged in 3D interpolation.

Open in new tabDownload slide

By employing pressure antisymmetric mirroring, IBM circumvents the limitations posed by regular cuboid meshes. In each time step of the FD process, determining the wavefield at the ghost point is essential to facilitate the FD computation near the surface boundary. The wavefield value to ghost points will be set to the opposite value to its corresponding mirror points, and enables the implicit implementation of the free surface boundary condition, ensuring the surface pressure component equals zero.

However, it is noteworthy that mirror points are commonly situated on the fractional grids. Thereby, precise calculation of the wavefield at these mirror points in 3D space is a pivotal aspect of the 3DIBM.

2.2. Computing the wavefield at mirror points in 3D space

Figure 1 illustrates that mirror points are situated on non-integer spatial grids. This situation presents difficulties for traditional FDM that rely on cuboid grids. This section explains how to compute wavefields at mirror points in 3D subsurface areas, which are often on fractional spatial grids. To address this, we examined several interpolation techniques, including Bilinear, Lagrange, Kaiser-windowed sinc, etc. Based on our tests, Kaiser-windowed sinc interpolation method is the most accurate one. It applies a Kaiser window to truncate the sinc function into a finite length (Hicks 2002, Oppenheim and Schafer 2010, Li et al. 2020). We adapted this method for use in 3D spaces. The implementation includes four steps as follows:

• Mirror point centring: the mirror point is positioned at the centre of the 3D sinc function. As the sinc function is discretized into an odd number serial, we select four integer sample points in each unilateral direction, resulting in a total of nine points for interpolation in one space direction. Therefore, in total 729 (9 × 9 × 9) integer points are used in the 3D interpolation for each mirror point.

  The sinc function is an even function defined as

  $$\begin{eqnarray} {\mathop{\rm{sinc}}\nolimits} (x) = \frac{{\sin (\pi x)}}{{\pi x}}, \end{eqnarray}$$

  (4)

  where x is the distance between the mirror point and the integer grid points.
• Calculation of spatial interpolation coefficients: in step (i), nine sample points are used to truncate the continuous sinc function. To mitigate the truncation effect, we use the Kaiser function for nine sample points to smooth the truncation of the sinc function. The Kaiser function is defined as

  $$\begin{eqnarray} {\mathop{\rm{Kaiser}}\nolimits} (x) = \left\{ {\begin{array}{@{}*{2}{l}@{}} {\frac {{{I}_0}\left( {6.31\sqrt {1 - {{{\left(\displaystyle {\frac{x}{4}} \right)}}^2}} } \right)} {{{I}_0}(6.31)}}, &\quad { - r \le x \le r}\\ 0\quad \quad \quad \quad \quad \ \ , &\quad {{\rm{otherwise}}} \end{array}} \right., \end{eqnarray}$$

  (5)

  where |${{I}_0}$| denotes the first kind of zero-order Bessel function, r is the half-width of the Kaiser window. In our calculation, it is set to four to achieve a good trade-off between efficiency and accuracy (Hicks 2002). Figure 2 illustrates them for a 2D case of the xy plane. For 3D interpolation, coefficients are cumulatively determined across an additional dimension.
• Wavefield interpolation at mirror points: the interpolation formula for wavefield at every mirror point in 3D space is as

  $$\begin{eqnarray} &&{{P}_{3d \hbox{-} \textit{interp}}}\left( {x,y,z} \right) = \sum\limits_{ix = - 4}^4 {\sum\limits_{iy = - 4}^4 {\sum\limits_{iz = - 4}^4}}\\ &&\quad \times {{{\left( {\frac{{\sin (\pi d)}}{{\pi d}} \cdot {\mathop{\rm{Kaiser}}\nolimits} (d) \cdot P\left( {ix,iy,iz} \right)} \right)} } } , \end{eqnarray}$$

  (6)

  where |$ix$|⁠, |$iy$|⁠, |$iz$| are the integer grid numbers of the 3D area surrounding the mirror point, d is the distance between the integer grid points |$( {ix,iy,iz} )$| and the mirror point |$( {x,y,z} )$|⁠, |$P( {ix,iy,iz} )$| denotes the pressure component at a specified integer node around the mirror point during the current modelling timestep. |${{P}_{3d \hbox{-} \textit{interp}}}( {x,y,z} )$| represents the pressure component at the mirror point after the application of 3D Kaiser-windowed sinc interpolation within the spatial domain.

  The primary computational cost is associated with calculating the Kaiser-windowed sinc function. To enhance efficiency for multi-shot 3D simulations, we propose pre-computing these interpolation coefficients with preset subgrid intervals. For instance, with a 20 subgrid point precision, we precompute 9261 (=21³) sets of coefficients and approximate a mirror point with the closest subgrid point. As an illustrative example, consider the spatial grids of 10 m in three directions (x, y, and z) in spatial, the subgrid interval is 0.5 m, which is sufficient for exploration seismic wavefield modelling. The pre-computation strategy helps reduce the cost of 3D sinc interpolation significantly.

• Iterative wavefield update at mirror points: the interpolation process may involve points in the air layer above the 3D surface boundary, indicated as ghost points, where the wavefield value is unknown. We initially set them to zero. After each interpolation iteration, the wavefield at the ghost point is adjusted to the opposite value of their corresponding mirror point, which could be expressed as follows:

  $$\begin{eqnarray} p_{{ghost}}^{n + 1}\,\, = - p_{{mirror}}^n, \end{eqnarray}$$

  (7)

  where |$p_{{ghost}}^{n + 1}\,$| signifies the pressure component at ghost points during the |$n + 1$|th iteration, |$p_{{mirror}}^n$| indicates the pressure component at the mirror points during the nth interpolation iteration, which is equal to |${{P}_{3d \hbox{-} \textit{interp}}}( x )$| in Equation (6). This iteration continues until the wavefield at ghost points converges, usually within 10 iterations in 3D cases. So that the ghost point can participate in calculating the next FD time step near the 3D surface.

Figure 2.

The 2D sinc interpolation coefficients for xy plane. The colour surface is the continuous 2D sinc function, the yellow dots are the discrete sinc coefficients, and the grey surface is the Kaiser window.

Open in new tabDownload slide

2.3. Computing the particle velocity component at the surface

Geophones record the particle vibration velocity component for land surveys at or close to the surface. Additionally, the surface topography is typically not flat, resulting in the placement of receivers on the fractional grids. To obtain the land seismic records, the following first-order Equation (8) is employed to convert the near-surface pressure component into particle vibration velocity:

where |${{\bf v}}$| denotes the vector of vibration velocity, comprising three components |${{v}_x}$|⁠, |${{v}_y}$| and |${{v}_z}$|⁠, corresponding to the x, y, and z directions of a three-dimensional model, respectively. Subsequently, the 3D interpolation method outlined in Section 2.2 can be applied to calculate the particle velocity at the locations of the receivers.

The flowchart in Fig. 3 outlines the sequential steps involved in 3DIBM for forward modelling in seismic wavefield simulations. It begins with the initialization of the 3D model and grids. The next step is identifying ghost points, which are necessary for implementing FD stencil near the free surface boundary. Following this, mirror points are determined (Fig. 1), essential for accurate wavefield interpolation. Then, the normal FD calculations are performed below the surface to simulate wave propagation. Then, the wavefield at mirror points is interpolated using the Kaiser-windowed sinc interpolation, ensuring precision even on fractional spatial grids. This is followed by an iterative process where the wavefield at the ghost points is updated until it converges. Once converged, the particle velocity at the surface is computed, crucial for converting near-surface pressure components into observable data. Finally, the simulated seismic data velocity components are output and saved.

Figure 3.

The flowchart of 3DIBM for forward modelling in seismic wavefield simulations.

Open in new tabDownload slide

3. Numerical examples

To demonstrate the effectiveness of our 3DIBM in this paper, we conducted modelling on two 3D models with irregular topography.

3.1. The 3D sloped model

To evaluate the efficacy of the proposed 3DIBM, we employed a 3D homogenous model characterized by an inclined surface and a uniform propagation P-wave velocity of 3000 m/s, as depicted in Fig. 4. The 3D surface plane of this model passes through the three points with coordinates (0, 0, 300), (0, 1500, 500), and (1500, 1500, 700), hence the model is a true 3D topography, rather than a pseudo-3D model that remains uniform in the y-axis direction. The model is discretized into a grid of 151 × 151 × 151 cells, with a spatial resolution of 10 m across the x, y, and z axes. Seismic wavefield propagation is simulated using a 30 Hz dominant frequency Ricker wavelet. The spatial sampling interval is 10 m, and the modelling time sampling interval is 1 ms. The seismic source, denoted by a yellow dot in Fig. 4, and the receiver arrays are uniformly distributed within the model at a depth of 700 m along the plane defined by x = 750 m and y = 750 m.

Figure 4.

A 3D homogeneous model featuring sloped surface topography. The yellow dot signifies the source location, which is located at the coordinate (x = 750 m, y = 750 m, depth = 750 m), while the black dashed lines denote the positions of two receiver arrays.

Open in new tabDownload slide

Figure 5 illustrates the spatial configuration of two types of point near the slope surface in a 3DIBM. The coloured surface in the figure represents a non-flat, sloped topography. It distinguishes between two kinds of point: ‘ghost points’, which are shown as coloured, unfilled circles arranged in four layers, and ‘mirror points’, which are represented by solid dots and correspond to the ghost points. To compute the wavefield at a mirror point, interpolation is required near the surface boundary in three-dimensional space, as described in Equation (6). This process necessitates a minimum of four layers of ghost points to ensure sufficient data for half the interpolation window radius. In simulations that use a FD scheme with a spatial order greater than eight, the number of required ghost point layers increases to half the order of the FD scheme. The illustration also shows that ghost points are aligned with the integer grid coordinates in 3D space, while mirror points are generally found at fractional grid coordinates. An effective absorbing condition is applied around the model except for the surface (Yao et al. 2018b, Zhang and Yao 2024).

Figure 5.

(a) Enlarged view of the configuration of mirror and ghost points within the model depicted in Fig. 4. The coloured surface represents the free surface. Unfilled circles of various colours denote the four layers of ghost points, while filled circles represent their corresponding mirror points. Lines perpendicular to the free surface connect each ghost point to its mirror point. (b) A zoomed-in section of (a) with a spatial discrete grid overlay. Ghost points are located at the integer grid endpoints while mirror points are on the fractional grid within the meshes.

Open in new tabDownload slide

For comparative analysis, we conducted simulations of wavefields utilizing the traditional vacuum method with the irregular free surface, a prevalent approach in wavefield modelling as noted by Oprsal and Zahradnik (1999). The two wavefield snapshots show that when seismic waves propagate to the free surface and are reflected, the vacuum method generates many erroneous diffractions related to reflected waves in the free surface (Fig. 6a and c). Conversely, the wavefield snapshots depicted in Fig. 6b and d exhibit clarity, unequivocally demonstrating the capability of the 3DIBM in mitigating the erroneous diffractions induced by the ‘staircase-like’ configuration of the free surface.

Figure 6.

Snapshots of the 3D pressure wavefield utilizing the vacuum method (a, c) versus the 3DIBM (b, d), captured at 180 and 250 ms. Black arrows highlight the distinctions in wave propagation between the two methods.

Open in new tabDownload slide

Figure 7 presents seismic records obtained using both the vacuum method and 3DIBM. The vacuum method results in pronounced erroneous diffractions, whereas the 3DIBM produces clear, artefact-free direct and reflected waves. It is worth noting that the seismic records received by the two perpendicular receiver lines are identical because the receiver lines are symmetrical relative to the dipping line of the surface. This means that the 3DIBM is applicable to true 3D cases with irregular surface topography. Furthermore, the cyan dashed lines in Fig. 7 draw the precise traveltime curves of reflections from the free surface.

Figure 7.

Seismic profiles of the pressure wavefield comparing the vacuum method (a, c) with the 3DIBM (b, d). Cyan dashed lines illustrate the accurate traveltime of reflections from the free surface. The receiver array in (a, b) is along the line in y = 750 m, while the receiver array in (c, d) is along the line in x = 750 m.

Open in new tabDownload slide

A more detailed seismogram comparison is shown in Fig. 8. The peak of surface reflection calculated by the 3DIBM in red waveforms has a very high consistency with the theoretical traveltime indicated by the blue vertical bars. At the same time, there are almost no false diffracted waves after free surface reflected waves. By contrast, the vacuum method does not accurately represent the sloped surface, resulting in erroneous traveltime and artificial diffractions.

Figure 8.

Seismograms corresponding to records from Fig. 7a and b, located at the line of y = 750 m. Blue vertical bars mark the theoretical traveltimes for reflections from the free surface.

Open in new tabDownload slide

3.2. The 3D ‘dual-complexity’ geological model of Tarim Basin

Figure 9 showcases a 3D ‘dual-complexity’ model, designed for the exploration of the Tarim Basin in Western China, as reported by Li et al. (2023). This model was established to research seismic wavefield simulation and imaging caused by complex surface topography and complex subsurface structures.

Figure 9.

A ‘double-complex’ 3D geological model of Tarim Basin showcasing (a) the velocity model and (b) an enlarged view of the 3D surface topography. The source location is marked by a yellow dot, which is located at the coordinate (x = 1275 m, y = 1950 m, depth = 338 m), and black dashed lines delineate the boundaries of the surface topography.

Open in new tabDownload slide

The 3D model, as depicted in Fig. 9a, is discretized into grids measuring 701 by 301 by 351 in the x, y, and depth dimensions, respectively, with a consistent spacing interval of 15 metres between each grid. We used a Ricker wavelet as the seismic source, which has a main frequency of 15 Hz and is initiated at the surface. Figure 10 displays two perpendicular slices of the wavefield through the source point. In these simulations, the 3DIBM clearly depicts the surface reflected waves, whereas the vacuum method results show these reflections are contaminated by artificial diffractions.

Figure 10.

3D seismic wavefield snapshots using the vacuum method (a, c, e) and their 3DIBM counterparts (b, d, f), with different recording times indicated in each subfigure. The wavefield slices are positioned at x = 1275 m and y = 1950 m, traversing the source location marked by the yellow dot. Black arrows and dashed-line boxes identify differences in wave propagation between the two methods.

Open in new tabDownload slide

Figures 11 and 12 display the shot gathers, showcasing three-component modelling results. The continuity of seismic events in the profiles simulated from 3DIBM is significantly better than that of the vacuum method, showcasing the mitigation of artificial diffractions and, hence clearer reflection signals. This affirms the 3DIBM's improvement in simulation quality, particularly for irregular topography, making it a better 3D forward modelling engine for land data inversion and imaging.

Figure 11.

Seismic profiles were recorded at the receiver line of x = 1050 m, displaying particle vibration velocity components from the vacuum method (a, c, e) and 3DIBM (b, d, f). Black arrows and dashed-line boxes highlight discrepancies in the 3D modelling outcomes.

Open in new tabDownload slide

Figure 12.

Seismic profiles were recorded at the receiver line of x = 1800 m, presenting particle vibration velocity components from the vacuum method (a, c, e) and 3DIBM (b, d, f). Black arrows and dashed-line boxes identify differences in the 3D modelling results.

Open in new tabDownload slide

To further verify the accuracy of the simulation results, we carried out the 3D Kirchhoff pre-stack depth migration (PSDM) based on ray theory for the middle line of x = 2250 m. Figure 13 shows the P-wave velocity section and the corresponding reflection coefficient of the target line.

Figure 13.

The P-wave velocity profile (a) and its corresponding reflection coefficient profile (b) along the target line. Black dashed lines denote the surface topography.

Open in new tabDownload slide

In the stage of simulating synthetic data, we fired seismic sources at the first integer grid beneath the surface at a horizontal interval of 75 m, a total of 141 sources along the target line. For recording seismic data, we evenly placed all the receivers on the free surface with a horizontal interval of 15 m and a total of 701 receivers for one receiver line. We took the target line at x = 2250 m as the centre line and set up a total of 11 receiver lines to record 3D seismic data. The source wavelets are Rickers with 15 and 20 Hz dominant frequencies. The vacuum method and the 3DIBM are used to simulate the records. Thus, there are four sets of seismic data. Figure 14 shows the CMP gathers used for Kirchhoff PSDM imaging. We can see that the noise of the 3DIBM simulated data is significantly less than that of the vacuum method simulated data. These differences will be reflected in further imaging results. The simulation process of the total 141 shots reflects the stability of the proposed 3DIBM.

Figure 14.

The CMP gathers of the four sets of seismic data. The first row is two data sets with a 15 Hz Ricker wavelet, while the second row is for 20 Hz. (a) and (c) The CMP gathers generated by the vacuum method. (b) and (d) The CMP gathers from the 3DIBM. The red arrows indicate the data difference.

Open in new tabDownload slide

We used the small-scale smoothed true P-wave velocity for the 3D Kirchhoff PSDM migration, and set the maximum migration aperture to 6000 m. The imaging datum is defined as the true surface topography. After PSDM imaging, the migrated common reflection points (CRP) gathers are shown in Fig. 15. As can be seen, the CRP gathers of the 3DIBM are superior to the vacuum method in terms of noise level and event continuity.

Figure 15.

CRP gathers at CDP = 4650 m for the vacuum method (a, c) and 3DIBM (b, d). Differences induced by the two forward modelling approaches are indicated with black arrows. The modelling method and source wavelet frequency are specified in the top right corner of each panel.

Open in new tabDownload slide

The final migrated stack sections after PSDM are shown in Fig. 16. Same migration parameters and pre-stack mute are used in all migrations. We see that the imaging results of the 3DIBM simulated data are clearer and more accurate than those of the vacuum simulated data at both dominant frequencies. The imaging results of data generated by the vacuum method have obvious structural artefacts. By contrast, the imaging results of 3DIBM modelled data are improved from shallow to deep layers. It is foreseeable that this imaging improvement will be more obvious for more complex undulating surfaces.

Figure 16.

The stacked imaging sections of the target line after PSDM. (a) and (e) The vacuum method imaging results with dominant frequencies of 15 and 20 Hz, respectively. (c) and (g) The 3DIBM imaging results at these frequencies. The second column showing (b), (d), (f), and (h) are enlarged views of the regions within the black dashed boxes in (a), (c), (e), and (g), respectively, highlighting the disparities produced by the two modelling techniques.

Open in new tabDownload slide

We believe that choosing Kirchhoff PSDM better demonstrates the accuracy of the 3DIBM simulation operator compared to RTM. This is because RTM typically employs the same forward engine for both forward and backward simulation and imaging. At the same time, Kirchhoff PSDM provides cross-validation from a ray perspective, thereby confirming the accuracy of the 3DIBM forward modelling algorithm.

4. Conclusion

This study introduces a stable 3DIBM tailored for modelling acoustic wavefields with irregular topography using cuboidal meshes and second-order acoustic equation. Incorporating 3D Kaiser-windowed sinc interpolation and iterative symmetric interpolation in 3D space, the 3DIBM precisely computes wavefield values at mirror points. Two numerical examples demonstrate that the proposed 3DIBM is effective in mitigating artificial diffractions resulting from incorrect model discretization to the irregular surface and in improving reflection traveltime accuracy compared to traditional vacuum methods. Notably, the surface reflects wavefield traveltime aligns closely with theoretical expectations, and the 3DIBM effectively reduces the distortion of reflection events by strong diffractions in complex 3D models. PSDM tests underscore its value as a forward modelling tool for 3D land data inversion and imaging.

Funding

This research was supported by the basic research project of CNPC (grant no. 2023ZZ05-02).

Conflict of interest statement: None declared.

Data Availability

The data underlying this paper is available by contacting the corresponding author.

References