TY - JOUR
AU - An, Shengpei
AB - Abstract Numerical methods have been widely applied to simulate seismic wave propagation. However, few studies have focused on internal multiples modeling. The formation mechanism and response of internal multiples are still unclear. Therefore, we develop a weighted-optimised-based internal multiples simulation method under 3D conditions. Using a one-way wave equation and full-wavefield method, the different-order internal multiples are computed numerically in a recursive manner. The traditional Fourier finite-difference (FFD) method has low numerical accuracy in a horizontal direction. A globally optimised FFD (OFFD) method is used to improve the lateral propagation accuracy of the seismic waves. Meanwhile, we adopt an adaptive variable-step technique to improve computational efficiency. The 3D internal multiples modeling technique is capable of calculating the different-order multiple reflections in complex structures. We use the present method to simulate internal multiples in several models. Theoretical analyses are consistent with the numerical results. Numerical examples demonstrate that the 3D internal multiples modeling technique has superior performance when adapting to lateral velocity changes and steep dip. This also implies that our method is fit for the simulation of internal multiples propagation in a 3D complex medium and can assist in identifying the internal multiples from full-wavefield data. internal multiples, 3D numerical modeling, Fourier finite-difference, wave propagation, optimisation 1. Introduction Multiple reflections are very troublesome in the seismic data process, and harm imaging and interpretation. Due to the existence of a series of high- and low-velocity layers, strong internal multiples will be generated, especially deep areas where the energy of the multiple reflections may conceal the energy of the primaries. Similar structures are common in western China, such as the Jurassic coal seams in Tarim Basin and Yili Basin, and the Ordovician mudstone-shale layers in Sichuan Basin. Many suppression methods have been developed, including the inverse scattering series (Weglein et al. 1997; Innanen 2017), feedback-iterative-based surface-related multiples elimination (Berkhout & Verschuur 1997; Dragoset et al. 2010), the estimation of primaries by sparse inversion (Van Groenestijn & Verschuur 2009; Ypma & Verschuur 2013), the inversion-based method (Luo et al. 2003; Wang 2004, 2007), virtual seismic events method (Ikelle 2006; Liu et al. 2018), Marchenko method (Wapenaar et al. 2014; da Costa Filho et al. 2017) and the seismic interferometry method (Wapenaar & Fokkema 2006; An et al. 2017). Also, some scholars find that the multiples have smaller reflection angles and their traveltime is longer than that of primaries (Berkhout & Verschuur 2011). Multiples are also widely applied for imaging (Berkhout & Verschuur 2003; Liu et al. 2015; Nath & Verschuur 2020). However, whether in multiples suppression or imaging, it is impossible to explain the formation mechanism and characteristics of multiple reflections. As a consequence, understanding the response mechanism of seismic multiple reflections through numerical modeling becomes very important, especially in land seismic exploration. The internal multiples modeling method can separate different-order internal multiples directly, and the phase and amplitude of wavefields are consistent with the geological properties. Kennett (1979) first proposed an internal multiples modeling technique based on the reflectivity method, but that is unable to adapt to lateral changes in the medium. Covey et al. (1989) presented a ray-path tracing scheme with high computational efficiency to simulate internal multiples. However, it cannot deal with complex structures. Berkhout (2014) introduced a full-wavefield forward method that can model the wavefield in stages and handle relatively complex geological models. Verschuur and Davydenko applied the full-wavefield method to synthetic and real data imaging (Verschuur et al. 2016; Davydenko & Verschuur 2018, 2019). Kuang et al. (2020) developed an internal multiples modeling technique based on adaptive step-length-varying wavefield extrapolation, which is applied to generate internal multiples in a specific layer. The one-way wave equations are usually used to simulate the primaries, and corresponding methods include the phase-shift method (PSM) (Gazdag 1978; Gazdag & Sguazzero 1984), split-step Fourier method (SSFM) (Stoffa et al. 1990) and Fourier finite-difference (FFD) method (Ristow & Ruhl 1994; Zhu et al. 2018). The conventional FFD method adds a finite-difference correction term in SSFM, which improves the accuracy in the medium with steep dips and lateral velocity changes. Wang (2001) and Zhang et al. (2008) extended the Fourier finite-difference method to 3D conditions. The conventional FFD method maintains enough accuracy in the depth direction, but the accuracy in the lateral direction is insufficient. The primaries recorded by receivers are accurate at the near offset, but there exist errors at the far-offset seismic traces. For internal multiples, the traveling path is longer in the lateral direction than one of the primaries, so the signals at far-offset traces deviate from the real position more seriously, especially in the case of steep dip angles. For this reason, the conventional FFD method is not applicable to extrapolate the internal multiples. In this paper, inspired by the high-order equation reduction (Ma 1983) and weighted optimisation coefficient (Lee & Suh 1985), we present a high-order globally optimised FFD (OFFD) method to simulate the 3D internal multiple reflections. The presented method is qualified to simulate internal multiples in a 3D medium with high steep dip angles and strong lateral change. Also, the adaptive variable-step technique is introduced to further enhance computational efficiency. This paper is organised into three sections. First, we introduce the basic theories of the related methods. Then, we present a homogeneous model to verify the numerical accuracy and simulate the internal multiples in other complex models. Finally, brief conclusions for implementing the proposed method are summarised. 2. Methodology 2.1. 3D FFD scheme The one-way wave equation for a 3D heterogeneous medium in the frequency domain is formulated in a compact form: $$\begin{equation} \frac{{\partial \hat{u}\left( {\omega ;x,y,z} \right)}}{{\partial z}} = i{k_z}\hat{u}\left( {\omega ;x,y,z} \right), \end{equation}$$(1) where |$\hat{u}( {\omega ;x,y,z} )$| denotes the wavefield in the frequency domain, the notations |$x$|⁠, |$y$| and |$z$| denote the three directions of the coordinate axis in the space domain, respectively, |$\omega $| is the angular frequency, |$i$| is the imaginary number and |${k_z} = \sqrt {{{\omega ^{2}}/{v^{2}}} + {\partial ^{2}/{\partial {x^2}}} + {{\partial ^{2}}/\partial {y^2}}}$| (⁠|$v$| denotes medium velocity) represents the square-roots operator. Since the square-roots operator is not directly used to compute the wavefield, we must solve it at reference velocity |${v_0}$|⁠. Let |${k_z} = {k_{{z_0}}} + ( {{k_z} - {k_{{z_0}}}})$|⁠, where |${k_{{z_0}}}$| denotes the square-roots operator in reference velocity, and carry out Taylor series and continued-fraction expressions; so equation (1) is further split into the following three cascaded equations (Wang 2001; Zhang et al. 2008): $$\begin{equation} \frac{{\partial {{\hat{u}}_1}\left( {\omega ;x,y,z} \right)}}{{\partial z}} = i{k_{{z_0}}}{\hat{u}_1}\left( {\omega ;x,y,z} \right), \end{equation}$$(2) $$\begin{equation} \frac{{\partial {{\hat{u}}_2}\left( {\omega ;x,y,z} \right)}}{{\partial z}} = i\omega \Delta l{\hat{u}_2}\left( {\omega ;x,y,z} \right), \end{equation}$$(3) and $$\begin{equation}\frac{{\partial {{\hat{u}}_3}\left( {\omega ;x,y,z} \right)}}{{\partial z}} = \frac{{b\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)}}{{1{\rm{\ +\ }}a\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)}}{\hat{u}_3}\left( {\omega ;x,y,z} \right),\end{equation}$$(4) where equation (2) represents the phase-shift operator (Gazdag 1978). Equation (3) represents the slowness correction term (Staff et al. 1990), |$\Delta l = 1/v - 1/v_{0}$| denotes the slowness, equation (4) represents the finite-difference correction term (Ristow & Ruhl 1994; Huang & Michael 2000) and |$a$| and |$b$| are constant coefficients. Actually, equations (2)–(4) are equivalent to $$\begin{equation}{k_z} \approx {k_{{z_0}}} + \omega \Delta l + \frac{{b\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)}}{{1{\rm{\ +\ }}a\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)}}.\end{equation}$$(5) Equation (5) is regarded as a conventional FFD method. Compared with PSM and SSFM, the conventional FFD operator enhances the maximum propagation angle in wavefield modeling. But it guarantees that the angles within 75o are correct within the error range of only 1% (Huang & Michael 2000). 2.2. High-order weighted-optimised extrapolation operator The conventional FFD method uses only a second-order finite-difference correction term. When implementing Taylor series expansion for |${k_z}$|⁠, the higher-order terms are truncated. Moreover, if the high-order term is added to approximate the square-roots operator |${k_z}$|⁠, the discretisation for the high-order term is quite difficult and causes numerical instability. In view of all these arguments, a globally optimised FFD operator is presented to solve this problem. We add the high-order correction term into the conventional FFD operator and decompose it into a series of cascade second-order equations (Ma 1983). The least-squares scheme proposed by Lee & Suh (1985) is adopted to calculate the optimisation coefficients. The specific optimisation scheme is shown next. First, define |${u^2}$|⁠: as $$\begin{equation}{u^2} = - \frac{{{v^2}\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)}}{{{\omega ^2}}} = \frac{{{v^2}\left( {k_x^2 + k_y^2} \right)}}{{{\omega ^2}}} = {\sin ^2}\alpha, \end{equation}$$(6) where |${k_x}$| and |${k_y}$| denote the wavenumbers in horizontal and vertical directions, respectively, and |$\alpha $| represents the propagation angle with the vertical direction. Define |$r = \frac{{{v_0}}}{v}$|(⁠|$r \le 1$|⁠) where |$r$| denotes the inhomogeneity of medium, and |$\varepsilon $| denotes the error between the theoretical operator and approximate operator that reads: $$\begin{equation}\varepsilon = \frac{{\sqrt {1 - {u^2}} }}{{\left( {1 - r} \right)}} - \frac{{\sqrt {1 - {r^2}{u^2}} }}{{r\left( {1 - r} \right)}} + \frac{1}{r} + \left[ {\frac{{\sum\limits_{i = 1}^m {{a_i}{{\left( {{u^2}} \right)}^i}} }}{{1 - \sum\limits_{i = 1}^m {{b_i}{{\left( {{u^2}} \right)}^i}} }}} \right],\end{equation}$$(7) where |${a_i}$| and |${b_i}$| represent the optimisation coefficients in the ith second-order equation. Introduce |$A = \frac{{\sqrt {1\ - \ {u^2}} }}{{( {1\ -\ r } )}} - \frac{{\sqrt {1\ -\ {r^2}{u^2}} }}{{r( {1\ -\ r } )}} + \frac{1}{r}$|⁠, and minimise the integral of equation (7), so we can get $$\begin{equation}J = \int_{0}^{\varphi }{{\left\{ {A + \left[ {\frac{{\sum\limits_{i = 1}^m {{a_i}{{\left( {{u^2}} \right)}^i}} }}{{1 - \sum\limits_{i = 1}^m {{b_i}{{\left( {{u^2}} \right)}^i}} }}} \right]} \right\}\left( \theta \right)d\theta }},\end{equation}$$(8) where |$\varphi $| denotes the maximum propagation angle. Equation (8) indicates that the cumulative error is smallest within the largest propagation angle |$\varphi $|⁠. However, it is very difficult to solve equation (8) directly. To solve it, we adopt the weighted optimisation method to redefine a new weighted-optimised error as: $$\begin{equation}\varepsilon ^{\prime}{\rm{ = }}\sum\limits_{i = 1}^m {{a_i}{{\left( {{u^2}} \right)}^i}} {\rm{\ +\ }}A\left[ {1 - \sum\limits_{i = 1}^m {{b_i}{{\left( {{u^2}} \right)}^i}} } \right].\end{equation}$$(9) Using the least-squares method, we can construct a series of standard equations. The solutions of the equations represent the optimisation coefficients in equation (7). The specific formulations are constructed as follows: $$\begin{eqnarray}\sum\limits_{i = 1}^m {{a_i}} \int{{{{\left( {{u^2}} \right)}^{i + j}}d\theta }} &-& \sum\limits_{i = 1}^m {{b_i}} \int{{{{\left( {{u^2}} \right)}^{i + j}}\cdot Ad\theta }}\nonumber\\ &&= {-} \int{{A\centerdot {{\left( {{u^2}} \right)}^j}d\theta}},\end{eqnarray}$$(10-a) $$\begin{eqnarray}\sum\limits_{i = 1}^m {{a_i}} \int{{{{\left( {{u^2}} \right)}^{i + j}}\cdot Ad\theta }} &-& \sum\limits_{i = 1}^m {{b_i}} \int{{{{\left( {{u^2}} \right)}^{i + j}}\cdot {A^2}d\theta }}\nonumber\\ &&= {-} \int{{{A^2}\cdot {{\left( {{u^2}} \right)}^j}d\theta }}.\quad\quad\end{eqnarray}$$(10-b) Finally, we obtain the high-order optimised FFD (OFFD) operator, which reads: $$\begin{equation}k_z^{new} \approx {k_{{z_0}}} + \omega \Delta l + \sum\limits_i^m {\frac{{{b_i}\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)}}{{1{\rm{\ +\ }}{a_i}\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)}}} ,\end{equation}$$(11) where |$k_z^{new}$| represents the new FFD operator. To compare the property of the presented method and the traditional FFD method, the accuracy is quantitatively investigated by defining the relative error, which reads: $$\begin{equation}\varepsilon _{FFD}^m\left( \alpha \right) = \frac{{\left| {{\varepsilon ^m}} \right|}}{{r\sqrt {1 - {{\sin }^2}\left( \alpha \right)} }},\end{equation}$$(12) where |$\varepsilon _{FFD}^m$| represents the relative error of the 2m-order optimised FFD operator. Figure 1 is the comparison of the relative errors for the optimised FFD method. The dashed line in the horizontal direction represents the relative error accepted within 1%. We can observe that the maximum propagation angle for the fourth-order optimisation reaches 87o, and the maximum propagation angle above the sixth-order optimisation method reaches 90o. It shows that the high-order optimised FFD method possesses higher accuracy than the conventional FFD method, because the high-order FFD operator adds more finite-difference terms to make the truncation error smaller. In general, the fourth- or sixth-order optimised FFD operator can satisfy the accuracy requirements for the internal multiples modeling. Following this, we use the sixth-order optimised FFD operator to simulate internal multiples in different media. Figure 2 denotes the relative error curves of the FFD operators in various inhomogeneous mediums. We find that as the inhomogeneity |$r$| becomes smaller, the relative error gradually decreases and the dip angle becomes larger within the error range of 1%. Figure 1. Open in new tabDownload slide The relative error for different-order optimised FFD operators with the propagation angle (⁠|$r = 0.5$|⁠). Figure 1. Open in new tabDownload slide The relative error for different-order optimised FFD operators with the propagation angle (⁠|$r = 0.5$|⁠). Figure 2. Open in new tabDownload slide Relative error with the propagation angle in inhomogeneous media (the smaller the |$r$|⁠, the stronger the inhomogeneity). Figure 2. Open in new tabDownload slide Relative error with the propagation angle in inhomogeneous media (the smaller the |$r$|⁠, the stronger the inhomogeneity). By comparing with figures 1 and 2, we have a deeper conclusion that even if the inhomogeneity is weak in the low-order case, the maximum propagation angle cannot reach 90o, which is caused by the defect of the traditional FFD method itself. The high-order optimisation proposed by this paper can make up for this defect. Meanwhile, the maximum propagation angle in high-order cases is hardly affected by the inhomogeneity. 2.3. Internal multiples modeling The full-wavefield migration was first proposed by Berkhout and Verschuur (2011), which regards seismic data as full-wavefield data. Based on full-wavefield migration or inversion, we can obtain the underground structure. The full-wavefield method provides a good guide for the 3D internal multiples modeling. The generation mechanism and characteristics of internal multiples have been studied in 2D and 3D conditions (Berkhout 2014; Davydenko & Verschuur 2018). We propose a high-precision 3D internal modeling technique using a globally optimised Fourier finite-difference propagation operator. For 3D internal multiples modeling, a 3D model is split into a series of 2D planar thin plates along the vertical direction (z axis) (as shown in figure 3), and the multiples modeling is done recursively during each round trip. The wavefield propagation of each thin plate just needs to know the previous thin plate, instead of the entire model. The basic principle in each round trip is described as follows. Figure 3. Open in new tabDownload slide A sketch of the 3D domain decomposition. Figure 3. Open in new tabDownload slide A sketch of the 3D domain decomposition. As discussed by Berkhout (2014), wavefields are divided into downward and upward waves according to the vertical direction. The wavefield relationship in the propagation process is shown in figure 4. The superscript |$n$| indicates the nth interface, the subscript |$u$| denotes the upward wave and the subscript |$d$| denotes the downward wave. |$A_d^n$| and |$A_u^n$| indicate the downward and upward wavefields above the |${z_n}$| interface, respectively. |$B_d^n$| and |$B_u^n$| represent the downward and upward wavefields below the |${z_n}$| interface, respectively. Also, we introduce a propagation operator |$W( {{z_n},{z_{n - 1}}} )$| from the depth level |${z_{n - 1}}$| to the neighboring depth level |${z_n}$|⁠, which is the high-order optimised FFD operator described previously. Figure 4a shows the wavefield relationship of the downward extrapolation and figure 4b shows the wavefield relationship of the upward extrapolation. In every depth level |${z_n}$|⁠, the wavefield must satisfy: $$\begin{equation}B_d^n = T_d^nW\left( {{z_n},{z_{n - 1}}} \right)B_d^{n - 1} + R_u^nB_u^n,\end{equation}$$(13-a) $$\begin{equation}A_u^n = T_u^nW\left( {{z_n},{z_{n + 1}}} \right)A_u^{n + 1} + R_d^nA_d^n,\end{equation}$$(13-b) where |$T$| denotes transmission matrix and |$R$| denotes reflection matrix (assume that only pure P-wave reflection and transmission occur, satisfy the formulations |${R_p} = ( {{\rho _2}{v_2} - {\rho _1}{v_1}} )/( {{\rho _2}{v_2} + {\rho _1}{v_1}} )$| and |${T_p} = ( {2{\rho _1}{v_1}} )/( {{\rho _2}{v_2} + {\rho _1}{v_1}} )$|⁠, where |$\rho $| denotes mass density). Figure 4. Open in new tabDownload slide Wavefield relationship in the extrapolation process. (a) The downward extrapolation. (b) The upward extrapolation. Figure 4. Open in new tabDownload slide Wavefield relationship in the extrapolation process. (a) The downward extrapolation. (b) The upward extrapolation. During each round trip, based on the recursive iteration of equation (13a and b), the expressions of the downward and upward waves in the depth level |${z_n}$| are obtained: $$\begin{equation}A_d^n = W\left( {{z_n},{z_m}} \right)\left( {S_d^m + R_d^mA_d^m + R_u^mB_u^m} \right),\end{equation}$$(14-a) $$\begin{equation}B_u^n = W\left( {{z_n},{z_m}} \right)\left( {S_u^m + R_d^mA_d^m + R_u^mB_u^m} \right),\end{equation}$$(14-b) where |$S$| denotes the primary source and |$R_d^mA_d^m + R_u^mB_u^m$|can be understood as the secondary source excited at the |${z_m}$| interface, and it satisfies |$T_d^m = R_d^m + I$| and |$T_u^m = R_u^m + I$| (⁠|$I$| is the identity matrix). According to equation (14a and b), the wavefield continuously propagates and scatters. Finally, signals are recorded by the observation system on the surface. Equation (14a and b) only represents the computational process of the upward and downward wavefields in one round trip. The stepwise modeling of the seismic wavefield is a continuous iterative procedure, and the corresponding-order wavefield can be obtained by controlling the number of round trips. In this paper, the primary is considered as zero-order internal multiples. The iterative procedure for internal multiples modeling is shown in figure 5. Through the first-round trip, the zero-order internal multiples are computed (as shown in figure 5a). The dots represent the secondary sources excited during the extrapolation process. Each additional round trip can obtain the higher-order internal multiples. Green paths in figure 5b lead to the first-order internal multiples, and blue paths in figure 5c lead to the second-order internal multiples. To obtain higher-order internal multiples, more round trips are required. Based on this theoretical research, the high-precision modeling for 3D internal multiple reflections can be realised numerically. Figure 5. Open in new tabDownload slide Propagation paths of the internal multiples in different round trips (asterisk* represents the seismic source, and the inverted triangle represents receiver). (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 5. Open in new tabDownload slide Propagation paths of the internal multiples in different round trips (asterisk* represents the seismic source, and the inverted triangle represents receiver). (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. 2.4. Adaptive wavefield extrapolation scheme The FFD extrapolation operator is obtained by an approximate expansion, and the truncated errors are bound to exist. A large-step extrapolation will reduce the recursive error and save the computational cost than multiple small-step extrapolations (Liu et al. 2006). Due to the advantages of large-step extrapolation, the adaptive variable-step wavefield extrapolation technique is applied to improve the efficiency before without influencing accuracy. According to the complexity of the underground medium, a large step is applied in a simple region and a small step is applied in complex structures with lateral varying. The formulation of adaptive extrapolation is simplified as: $$\begin{equation}{S_n} = \left\| {{V_n} - {V_{n - 1}}} \right\|,\end{equation}$$(15) where |${S_n}$| denotes the nth layer factor in the model and |${V_n}$| represents the velocity vector at the nth layer. When |${S_n} = 0$|⁠, it means that the velocity vector in the nth layer is the same as the one in the |$( {n - 1} )$|th layer, and the wavefields in the n- and |$( {n - 1} )$|th layers are computed by one step. When |${S_n} \ne 0$|⁠, the velocity vector in the nth layer is regarded as a single extrapolation step. 3. Numerical experiments To investigate its accuracy, we construct homogeneous and graben models to test the proposed method in 2D situations. Then, the 3D internal multiples modeling method is applied to layered and inhomogeneous medium models to further investigate its validity. The 20-Hz Ricker wavelet is used as the excitation source. 3.1. Impulse responses In this section, a homogeneous medium model with a horizontal length of 6 km and a depth of 3 km is considered. The model velocity is 4 km s−1, and the background velocity is 2 km s−1, which denotes the inhomogeneity |$r = 0.5$|⁠. The source is located at the coordinate (0 km, 3 km). Figure 6 shows snapshots of wavefields for different-order optimised FFD operators at 0.5 s. The yellow dashed line indicates the location of the wavefront at 0.5 s under theoretical conditions. We can see that the accuracy of the optimised FFD method is significantly improved with an increase of the orders. The maximum propagation angle for the second-order FFD operator reaches 75o, and the maximum propagation angle for the fourth-order FFD operator reaches 87o. The sixth-order FFD method is close to 90o. The numerical results are consistent with the previous theoretical analysis, and this verifies that the proposed optimisation has higher accuracy and overcomes the shortcomings of the traditional FFD method. A dip filter scheme is used to eliminate the heart-shaped evanescent wave noise generated in the extrapolation process (Claerbout 1985). Based on the high-order propagation operator, the internal multiples are numerically simulated. Figure 6. Open in new tabDownload slide Snapshots of wavefields at 0.5 s. Medium velocity is |$v$| = 4 km s−1 and a reference medium with |${v_0}$| = 2 km s−1 was used to obtain these impulse responses. The yellow line represents a semicircle. (a) Second-order OFFD. (b) Fourth-order OFFD. (c) Sixth-order OFFD. (d) Eighth-order OFFD. Figure 6. Open in new tabDownload slide Snapshots of wavefields at 0.5 s. Medium velocity is |$v$| = 4 km s−1 and a reference medium with |${v_0}$| = 2 km s−1 was used to obtain these impulse responses. The yellow line represents a semicircle. (a) Second-order OFFD. (b) Fourth-order OFFD. (c) Sixth-order OFFD. (d) Eighth-order OFFD. 3.2. 2D graben model We use the optimised FFD method to simulate the internal multiples. A reflectivity model with a horizontal length of 9.6 km and a depth of 4 km is shown in figure 7. The density is uniformly 2 g cm−3. A source is excited at the location (0 km, 4.8 km). The recorded length is 3.5 s, the time interval is 5 ms and the trace interval is 0.005 km. The model size is 801 × 1921 with a 0.005-km grid size. According to the given parameters, we implement numerically internal multiples modeling. Figure 7. Open in new tabDownload slide Reflectivity structure of the graben model. Figure 7. Open in new tabDownload slide Reflectivity structure of the graben model. Figure 8 shows the different-order multiples seismic records obtained by the sixth-order optimised FFD. Figure 9 is the comparisons of multiples seismic records obtained by different numerical methods, where figure 9a is the full-wavefield record obtained by the finite-difference method, figure 9b shows the first-order internal multiples calculated by the second-order traditional FFD and figure 9c shows the first-order internal multiples calculated by the sixth-order optimised FFD. We can see that the internal multiples events obtained by the FFD are consistent with the ones obtained by the finite difference (indicated by the black arrow in figure 9). Moreover, we further compare the seismic traces at near and far offsets and find that at near offset (see figure 10), the seismic signals obtained by the traditional and sixth-order optimised FFD operators are consistent with the finite-difference signals. Whereas, at far offset (see figure 11), the seismic signals obtained by the second-order traditional FFD method have a certain delay, the sixth-order optimised FFD does not. This is caused by the insufficient propagation accuracy of low-order FFD in the lateral direction. Our method can make up for this shortcoming. Figure 8. Open in new tabDownload slide Seismic records for internal multiples. (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 8. Open in new tabDownload slide Seismic records for internal multiples. (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 9. Open in new tabDownload slide Comparison of seismic records with different numerical methods. (a) Finite difference. (b) First-order internal multiples with traditional FFD. (c) First-order internal multiples with 6th-order optimised FFD. Figure 9. Open in new tabDownload slide Comparison of seismic records with different numerical methods. (a) Finite difference. (b) First-order internal multiples with traditional FFD. (c) First-order internal multiples with 6th-order optimised FFD. Figure 10. Open in new tabDownload slide Comparisons of Single channel for three numerical methods at near offset (trace 1). (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 10. Open in new tabDownload slide Comparisons of Single channel for three numerical methods at near offset (trace 1). (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 11. Open in new tabDownload slide Comparisons of Single channel for three numerical methods at far offset (trace 961). (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 11. Open in new tabDownload slide Comparisons of Single channel for three numerical methods at far offset (trace 961). (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. To further investigate the characteristics of internal multiples, we superimposed different data volumes to get stacked profiles (as shown in figure 12). The shot interval is 20 m, and a total of 421 shots are implemented. Figure 12a is the stacked data obtained by finite difference, and figure 12 parts b–d show the stacked profiles of the zeroth- to second-order internal multiples. The black arrows refer to the stacked data of zeroth-order internal multiples corresponding to the finite-difference stacked data, the orange arrows refer to the stacked data of first-order internal multiples and the white arrow refers to the stacked data of second-order internal multiples. The internal multiples superposition profile corresponding to different orders can be matched with the profile of the finite difference, which proves the correctness of the numerical simulation. It also validates that the real seismic data is a superposition of different-order multiple reflections. Figure 12. Open in new tabDownload slide Stacked profiles based on (a) finite-difference method, (b) zeroth-order internal multiples, (c) first-order internal multiples, and (d) second-order internal multiples. (different colored arrows indicate the positions of the corresponding-order multiples). Figure 12. Open in new tabDownload slide Stacked profiles based on (a) finite-difference method, (b) zeroth-order internal multiples, (c) first-order internal multiples, and (d) second-order internal multiples. (different colored arrows indicate the positions of the corresponding-order multiples). 3.3. 3D five-layer medium model The 3D real seismic data are very complicated. Modeling the seismic waves in 3D structures can help real seismic data to accurately image them (Pu et al. 2015; Wang et al. 2017; Qi et al. 2020). Here, we use a 3D model to study the formation mechanism and response of the internal multiple reflections and help the identification of the internal multiples in real seismic data. A 3D five-layer media model is selected, the velocities from top to bottom are 3, 1.5, 3.5, 1.5 and 2.5 km s−1. The model contains two low-velocity layers, which can generate strong internal multiples. And the reflectivity is shown in figure 13. The source is located at the point (0 km, 2 km, 2 km). The model size is 301 × 401 × 401 grid points and the mesh size is 0.01 km. The receiver interval in the in-line direction is 0.02 km, the receiver interval in the cross-line direction is 0.2 km and medium density is 2 g cm−3. Figure 13. Open in new tabDownload slide 3D reflectivity model. Figure 13. Open in new tabDownload slide 3D reflectivity model. We use the 3D modeling technique to simulate the internal multiples propagation and get the related wavefields and synthetic seismograms. Figures 14 and 15 show the snapshots of the downward and upward wavefields of the zeroth-order multiples at a time level |$T\,\,{\rm{ = }}\,\,0.35$|⁠, 0.45, 0.55 and 0.65 s, respectively. The wavefields in the first-round trip propagate steadily in the computational domain. Figure 16 is the simulated seismograms of the first- and second-order internal multiples. We can observe that the events of the different phases are clear and the different-order multiples are accurately separated. Figure 14. Open in new tabDownload slide Snapshots of downward wavefields for zeroth-order internal multiples at a time level (a) |$T\,\,{\rm{ = }}\,\,0.35$| s, (b) |$T\,\,{\rm{ = }}\,\,0.45$| s, (c) |$T\,\,{\rm{ = }}\,\,0.55$| s, and (d) |$T\,\,{\rm{ = }}\,\,0.65$| s. Figure 14. Open in new tabDownload slide Snapshots of downward wavefields for zeroth-order internal multiples at a time level (a) |$T\,\,{\rm{ = }}\,\,0.35$| s, (b) |$T\,\,{\rm{ = }}\,\,0.45$| s, (c) |$T\,\,{\rm{ = }}\,\,0.55$| s, and (d) |$T\,\,{\rm{ = }}\,\,0.65$| s. Figure 15. Open in new tabDownload slide Snapshots of upward wavefields for zeroth-order internal multiples at a time level (a) |$T\,\,{\rm{ = }}\,\,0.35$| s, (b) |$T\,\,{\rm{ = }}\,\,0.45$| s, (c) |$T\,\,{\rm{ = }}\,\,0.55$| s, and (d) |$T\,\,{\rm{ = }}\,\,0.65$| s. Figure 15. Open in new tabDownload slide Snapshots of upward wavefields for zeroth-order internal multiples at a time level (a) |$T\,\,{\rm{ = }}\,\,0.35$| s, (b) |$T\,\,{\rm{ = }}\,\,0.45$| s, (c) |$T\,\,{\rm{ = }}\,\,0.55$| s, and (d) |$T\,\,{\rm{ = }}\,\,0.65$| s. Figure 16. Open in new tabDownload slide Simulated seismograms. (a) First-order one-shot record. (b) Second-order one-shot record. Figure 16. Open in new tabDownload slide Simulated seismograms. (a) First-order one-shot record. (b) Second-order one-shot record. To further investigate the properties of internal multiples, we compare the 3D multiples technique with the finite-difference method in this model. Figure 17 shows the comparison of the seismograms, which are calculated by the FFD method and finite-difference method, respectively, whereas figure 17 parts a and c are the sum of the multiples from 0 to 2 orders at near-offset and far-offset traces, respectively. The orange arrow denotes direct wave. Figure 18 is the zoomed-in view of figure 17. We can clearly observe that the synthetic seismograms are in good agreement with the ones computed by the finite difference. The phase of multiples is consistent with the finite-difference solution (pointed by the green and red arrows). These results validate the correctness and applicability of our method. Meanwhile, we quantitatively compare the relationship between the results of the multiples modeling method and the finite-difference results. When there is no normalisation, the amplitude of the different-order multiples is about 1.1 times that of the finite-difference results at near offset, and the amplitude is close to twice that of the finite-difference results at far offset. The travel times of the two methods are the same. Many trials are conducted to quantitatively investigate the reflections, similar results can be obtained. Figure 17. Open in new tabDownload slide Comparison of seismograms in the 3D model, which are calculated by the FFD method (left) and finite-difference method (right), respectively. Here, (a) and (b) denote the near-offset seismic traces, and (c) and (d) denote the far-offset seismic traces. Figure 17. Open in new tabDownload slide Comparison of seismograms in the 3D model, which are calculated by the FFD method (left) and finite-difference method (right), respectively. Here, (a) and (b) denote the near-offset seismic traces, and (c) and (d) denote the far-offset seismic traces. Figure 18. Open in new tabDownload slide Zoomed-in view. Here, (a) and (b) denote the near-offset seismic traces, and (c) and (d) denote the far-offset seismic traces. Figure 18. Open in new tabDownload slide Zoomed-in view. Here, (a) and (b) denote the near-offset seismic traces, and (c) and (d) denote the far-offset seismic traces. Also, we compare the computational efficiency of the different methods in the current model. The result is shown in figure 19. The adaptive modeling technique is the most efficient, and the finite-difference method is much more expensive. Because the finite-difference method must be applied on each grid point. The adaptive technique can automatically adjust the grid step size to compute the wavefields according to the simplicity of the model. It shows that the 3D internal multiples simulation method can effectively compute internal multiples in 3D media. Figure 19. Open in new tabDownload slide Comparison of the computational efficiency of different methods. Figure 19. Open in new tabDownload slide Comparison of the computational efficiency of different methods. 3.4. 3D inhomogeneous medium model In the final model, we verify the applicability of this method in 3D transverse heterogeneous media. The reflectivity structure is shown in figure 20, which is heterogeneous in the cross-line and in-line directions. The model size is 101 × 401 × 401 grid points and the grid size is 0.01 km. On the surface (0 km, 2 km, 2 km), a source is excited. The time length is 3 s with a 2 ms sampling interval. Figure 21 is the synthetic seismic records of internal multiples. We can clearly observe the different events of the internal multiples. This also indicates that the proposed method can effectively model the internal multiples in an inhomogeneous medium. A comparison between the finite difference and the proposed method is presented in figure 22. Seismic records calculated by finite difference are considered to be the standard reference result. The results of the two methods at near offset are relatively consistent. However, the amplitudes and waveforms are significantly different at far offset. Figure 20. Open in new tabDownload slide Reflectivity model in a 3D inhomogeneous medium. Figure 20. Open in new tabDownload slide Reflectivity model in a 3D inhomogeneous medium. Figure 21. Open in new tabDownload slide Synthetic seismograms computed by the optimised FFD operator. (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 21. Open in new tabDownload slide Synthetic seismograms computed by the optimised FFD operator. (a) Zeroth-order internal multiples. (b) First-order internal multiples. (c) Second-order internal multiples. Figure 22. Open in new tabDownload slide A comparison between the finite-difference results and the results of the multiples modeling method at (a) near (trace 1) and (b) far (trace 151) offsets. Figure 22. Open in new tabDownload slide A comparison between the finite-difference results and the results of the multiples modeling method at (a) near (trace 1) and (b) far (trace 151) offsets. 4. Conclusions In this paper, we develop a high-order globally optimised FFD operator to simulate the internal multiple reflections in a 3D medium because the conventional FFD method is not enough accurate in a lateral direction. To solve this problem, we add several finite-difference correction items to the conventional FFD operator for enhancing the numerical accuracy. Moreover, we apply the new FFD operator and full-wavefield method to simulate the 3D internal multiples. Through the wavefield extrapolation, a 3D model is divided into a series of thin plates along the vertical direction so that the simulation of 3D internal multiples can be realised numerically. The comparisons with the finite-difference method verify the correctness and practicability of our method; and the numerical results show that our method has superior performance to adapt to lateral velocity changes and steep dip, and the seismic record is a superposition of different-order multiple reflections. In addition, the 3D simulation technique may be applied to identify internal multiples noises in real seismic data and pseudo-stratigraphic structures in post-stack data. Acknowledgements All authors thank the National Key Research and Development Program of China (2018YFA0702503, 2018YFA0702505) and the National Natural Science Foundation of China (41674122) for supporting this research. Furthermore, all authors appreciate the constructive comments and suggestions on the manuscript by the editors and three anonymous reviewers. Conflict of interest statement. None declared. References An S.P. , Hu T.Y., Liu Y.M, Peng G.X., Liang X.H., 2017 . Automatic first-arrival picking based on extended super-virtual interferometry with quality control procedure , Exploration Geophysics , 48 , 124 – 130 . Google Scholar OpenURL Placeholder Text WorldCat Berkhout A.J. , 2014 . Review paper: an outlook on the future of seismic imaging, Part I: forward and reverse modelling , Geophysical Prospecting , 62 , 911 – 930 . Google Scholar OpenURL Placeholder Text WorldCat Berkhout A.J. , Verschuur D.J., 1997 . Estimation of multiple scattering by iterative inversion; Part 1, Theoretical considerations , Geophysics , 62 , 1586 – 1595 . Google Scholar OpenURL Placeholder Text WorldCat Berkhout A.J. , Verschuur D.J., 2003 . Transformation of multiples into primary reflections , 73rd SEG Annual International Meeting, Expanded Abstract , 1925 – 1928 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Berkhout A.J. , Verschuur D.J., 2011 . Full wavefield migration, utilizing surface and internal multiple scattering , 81th SEG Annual International Meeting Expanded Abstract , 3212 – 3216 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Claerbout J.F. , 1985 . Imaging the Earth's Interior . Blackwell Scientific Publications, Inc. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Covey J.D. , Hron F., Peacock K.L., 1989 . On the role of partial ray expansion in the computation of ray synthetic seismograms for layered structures , Geophysical Prospecting , 37 , 907 – 923 . Google Scholar OpenURL Placeholder Text WorldCat da Costa Filho C.A. , Meles G.A., Curtis A., 2017 . Elastic internal multiple analysis and attenuation using Marchenko and interferometric methods , Geophysics , 82 , Q1 – Q12 . Google Scholar OpenURL Placeholder Text WorldCat Davydenko M. , Verschuur D.J., 2018 . Including and using internal multiples in closed-loop imaging field data examples , Geophysics , 83 , R297 – R305 . Google Scholar OpenURL Placeholder Text WorldCat Davydenko M. , Verschuur D.J., 2019 . Using the full wavefield both in FWI & wavefield tomography , 89th SEG Annual International Meeting, Expanded Abstract , 1239 – 1243 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Dragoset B. , Verschuur E., Moore I., Bisley R., 2010 . A perspective on 3D surface-related multiple elimination , Geophysics , 75 , A245 – A261 . Google Scholar OpenURL Placeholder Text WorldCat Gazdag J. , 1978 . Wave equation migration with the phase-shift method , Geophysics , 43 , 1342 – 1351 . Google Scholar OpenURL Placeholder Text WorldCat Gazdag J. , Sguazzero P., 1984 . Migration of seismic data by phase shift plus interpolation , Geophysics , 49 , 124 – 131 . Google Scholar OpenURL Placeholder Text WorldCat Huang L.J. , Michael C., 2000 . Globally optimised Fourier finite-difference migration method , 70th SEG Annual International Meeting, Expanded Abstract , 802 – 805 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Ikelle L.T. , 2006 . A construct of internal multiples from surface data only: the concept of virtual seismic events , Geophysical Journal International , 164 , 383 – 393 . Google Scholar OpenURL Placeholder Text WorldCat Innanen K.A. , 2017 . Time- and offset-domain internal multiple prediction with nonstationary parameters , Geophysics , 82 , V105 – V116 . Google Scholar OpenURL Placeholder Text WorldCat Kennett B.L.N. , 1979 . Theoretical reflection seismograms for elastic media , Geophysical Prospecting, 27 , 301 – 321 . Google Scholar OpenURL Placeholder Text WorldCat Kuang W.K. , Hu T.Y., Duan W.S., Li J.S., Li Y.D., 2020 . Modeling inter-layer multiples based on adaptive step-length-varying wavefield extrapolation , Chinese Journal of Geophysics (in Chinese) , 63 , 2043 – 2055 . Google Scholar OpenURL Placeholder Text WorldCat Lee M.W. , Suh S.Y., 1985 . Optimization of one-way wave equations , Geophysics , 50 , 1634 – 1637 . Google Scholar OpenURL Placeholder Text WorldCat Liu H. , Yuan J.H., Chen J.B., Shou H., Li Y.M., 2006 . Theory of large-step wavefield depth extrapolation , Chinese Journal of Geophysics (in Chinese) , 49 , 1779 – 1793 . Google Scholar OpenURL Placeholder Text WorldCat Liu J.H. , Hu T.Y., Peng G.X., Cui Y.F., 2018 . Removal of internal multiples by iterative construction of virtual primaries , Geophysical Journal International , 215 , 81 – 101 . Google Scholar OpenURL Placeholder Text WorldCat Liu Y. , Hu H, Xie X.B., Zheng Y., Li P., 2015 . Reverse time migration of internal multiples for subsalt imaging , Geophysics , 80 , S175 – S185 . Google Scholar OpenURL Placeholder Text WorldCat Luo Y. , Kelamis P.G., Wang Y.H., 2003 . Simultaneous inversion of multiples and primaries: Inversion versus subtraction , The Leading Edge , 22 ( 9 ), 814 – 891 . Google Scholar OpenURL Placeholder Text WorldCat Ma Z.T. , 1983 . Splitting algorithm for high-order equation migration , Chinese Journal of Geophysics (in Chinese) , 26 , 377 – 389 . Google Scholar OpenURL Placeholder Text WorldCat Nath A. , Verschuur D.J., 2020 . Imaging with surface-related multiples to overcome large acquisition gaps , Journal of Geophysics and Engineering , 17 , 742 – 758 . Google Scholar OpenURL Placeholder Text WorldCat Pu Y.T. , Yin X.Y., Cao D.P., Zhang C., 2015 . The small-scale forward modeling method for large models based on the wave equation , Journal of Geophysics and Engineering , 12 , 988 – 995 . Google Scholar OpenURL Placeholder Text WorldCat Qi Z.P. , Li X., Fan K.R., Zhi Q.Q., 2020 . The three-dimensional finite element forward modelling of complex excitation source NMR , Journal of Geophysics and Engineering , 17 , 127 – 137 . Google Scholar OpenURL Placeholder Text WorldCat Ristow D. , Ruhl T., 1994 . Fourier finite-difference migration , Geophysics , 59 , 1882 – 1893 . Google Scholar OpenURL Placeholder Text WorldCat Stoffa P.L. , Fokkema J.T., Freire R., Kessinger W.P., 1990 . Split-step Fourier migration , Geophysics , 55 , 410 – 421 . Google Scholar OpenURL Placeholder Text WorldCat Van Groenestijn G.J. , Verschuur D.J., 2009 . Estimating primaries by sparse inversion and application to near-offset data reconstruction , Geophysics , 74 , A23 – A28 . Google Scholar OpenURL Placeholder Text WorldCat Verschuur D.J. , Staal X.R., Berkhout A.J., 2016 . Joint migration inversion: simultaneous determination of velocity fields and depth images using all orders of scattering , The Leading Edge , 35 , 1037 – 1046 . Google Scholar OpenURL Placeholder Text WorldCat Wang W.Z. , Hu T.Y., Song J.Y., Li Y.D., Li J.S., Zhang Y., 2017 . An optimized scheme of dispersion suppression for elastic-wave variable-order rotated staggered-grid forward modeling , Journal of Geophysics and Engineering , 14 , 1624 – 1638 . Google Scholar OpenURL Placeholder Text WorldCat Wang Y. , 2001 . ADI plus interpolation: accurate finite-difference solution to 3D paraxial wave equation , Geophysical Prospecting , 49 , 547 – 556 . Google Scholar OpenURL Placeholder Text WorldCat Wang Y. , 2004 . Multiple prediction through inversion: a fully data-driven concept for surface-related multiple attenuation , Geophysics , 69 , 547 – 553 . Google Scholar OpenURL Placeholder Text WorldCat Wang Y. , 2007 . Multiple prediction through inversion: theoretical advancements and real data application , Geophysics , 72 , V33 – V39 . Google Scholar OpenURL Placeholder Text WorldCat Wapenaar K. , Fokkema J., 2006 . Green's function representations for seismic interferometry , Geophysics , 71 , SI33 – SI46 . Google Scholar OpenURL Placeholder Text WorldCat Wapenaar K. , Thorbecke J., Van Der Neut J., Broggini F., Slob E., Snieder R., 2014 . Marchenko imaging , Geophysics , 79 , WA39 – WA57 . Google Scholar OpenURL Placeholder Text WorldCat Weglein A.B. , Gasparotto F A, Carvalho P.M., Stolt R.H., 1997 . An inverse-scattering series method for attenuating multiples in seismic reflection data , Geophysics , 62 , 1975 – 1989 . Google Scholar OpenURL Placeholder Text WorldCat Ypma F.H.C. , Verschuur D.J., 2013 . Estimating primaries by sparse inversion, a generalized approach , Geophysical Prospecting , 61 , 94 – 108 . Google Scholar OpenURL Placeholder Text WorldCat Zhang J.H. , Wang W.M., Fu L.Y., Yao Z.X., 2008 . 3D Fourier finite-difference migration by alternating-direction-implicit plus interpolation , Geophysical Prospecting , 56 , 95 – 103 . Google Scholar OpenURL Placeholder Text WorldCat Zhu F. , Huang J.P., Yu H., 2018 . Least-squares Fourier finite-difference pre-stack depth migration for VTI media , Journal of Geophysics and Engineering , 15 , 421 – 437 . Google Scholar OpenURL Placeholder Text WorldCat © The Author(s) 2021. Published by Oxford University Press on behalf of the Sinopec Geophysical Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
TI - 3D adaptive internal multiples modeling based on globally optimised Fourier finite-difference method
JF - Journal of Geophysics and Engineering
DO - 10.1093/jge/gxab022
DA - 2021-07-27
UR - https://www.deepdyve.com/lp/oxford-university-press/3d-adaptive-internal-multiples-modeling-based-on-globally-optimised-2DsKHi7fUV
SP - 429
EP - 445
VL - 18
IS - 4
DP - DeepDyve
ER -