1. Introduction
The increasing complexity of geological structures necessitates advancements in oil and gas exploration and development technologies, making multiple suppression an essential and challenging problem to address. The presence of multiples interferes with the identification of primary reflections, reducing the signal-to-noise ratio of seismic data, compromising the accuracy and reliability of seismic imaging, and potentially leading to erroneous interpretations and source rock investigations. Consequently, effective suppression of multiples is crucial for enhancing seismic data quality and facilitating oil and gas exploration and development. In northwest China, there are abundant resources of oil and natural gas. However, the complexity of its geological features and strong reflective interfaces in underground media, such as rock mounds, coal seams, basalt, etc., can generate strong internal multiples. The understanding of internal geological structures is constrained. Therefore, the demand for effective suppression of seismic multiples is very urgent.
Seismic multiples are generally classified into two categories based on the location of reflective layers: surface-related multiples and internal multiples. Surface-related multiples exhibit relatively stronger energy, particularly in offshore seismic data. Consequently, numerous researchers have focused on predicting and suppressing surface-related multiples, resulting in the maturation of suppression techniques for this category. In 1992, Verschuur et al. employed the concept of series expansion to decompose the seismic wavefield into primaries and multiples of all orders, effectively suppressing surface-related multiples [
1]. This method obviates the need for subsurface geological model information, is a data-driven method, delivers superior results in field data processing, and exhibits high computational efficiency. Many geophysicists have expressed interest in and conducted research on this approach, leading to its widespread adoption in the exploration and development field. In 2009, Van Groenestijn and Verschuur proposed using sparse inversion to estimate primaries, eliminating the requirement for adaptive matched subtraction and near offset extrapolation [
2,
3]. This advancement significantly optimized the traditional surface-related multiple elimination (SRME) methods. In 2015, Ma et al. developed a suppression method for 3D surface-related multiples [
4]. This technique assumes that the multiple contribution trace set corresponding to each seismic record in 3D seismic data is hyperbolic, enabling multiple predictions and suppression in 3D seismic data through sparse inversion. In 2022, Zhu et al. introduced a least squares datum continuation method for suppressing surface-related multiples within the least squares inversion theoretical framework [
5]. This approach targets seafloor-separated up- and down-wave records, and, without multiple separation or fractional extraction, iterative inversion yields pre-stack seismic data from virtual observations on the seafloor, enhancing the resolution of effective signals and improving seismic imaging quality.
In comparison to surface-related multiples, internal multiples exhibit smaller differences in frequency, stacking velocity, and normal moveout (NMO) correction, making their suppression more challenging. In 1998, Jakubowicz introduced a data-driven method for internal multiple suppression [
6]. This approach directly employs surface observed seismic data to construct internal multiples, enhancing computational efficiency. However, the method necessitates the accurate selection of primary reflection events from seismic data, which is difficult to implement in complex field data. In 2005, Berkhout and Verschuur proposed an internal multiple prediction method based on common focus technology, which is grounded in wavefield extrapolation [
7,
8]. This method is better suited to complex geological conditions but relies on a velocity model, limiting its application in field data. In 2006, Weglein et al. presented the inverse scattering series method based on the point scattering model [
9]. As a data-driven, wave equation-based method, it does not require prior information such as velocity models. However, the method involves substantial computation and is only effective for near offsets, making it insufficient for current data processing demands. In 2013, Ypma and Verschuur introduced an internal multiple suppression method based on sparse inversion [
10]. This approach minimizes the damage caused by adaptive matched subtraction to primary signals and effectively protects effective waves. Nevertheless, the method still entails considerable computation and yields unstable wavelets, complicating large-scale field data processing. In 2015, inspired by Marchenko’s imaging method, Meles et al. proposed a self-focusing-based Marchenko internal multiple suppression method [
11]. This technique does not require accurate velocity models, only an inaccurate macro velocity model to forward direct waves. Subsequently, direct waves and original data are used to construct the up- and down-going Green functions for relevant virtual source points, generating the associated internal multiples. However, this method demands significant computation due to the need for virtual source points at various depths, making it difficult to apply widely to large scale, complex structured seismic data. In 2018, Zhang and Staring proposed a one-step internal multiple suppression method based on Marchenko’s self-focusing approach, also known as the data domain internal multiple suppression method [
12]. This technique does not require macro velocity models for direct wave estimation and directly outputs primaries without internal multiples. However, this method imposes high requirements on the observation system, and its large-scale application in field data remains limited. In 2020, Zhang et al. applied the Marchenko internal multiple suppression method to field seismic data processing, achieving favorable results [
13]. Nonetheless, in low signal-to-noise ratio scenarios, energy may not be sufficiently focused, rendering the Marchenko internal multiple suppression method ineffective. In 2022, Zhang et al. compared the self-focusing-based Marchenko method for internal multiple suppression with the data domain method [
14]. In the same year, Zhang et al. introduced a multiple suppression method based on self-attention convolutional auto encoders [
15]. This method can reduce artificial cost, reduce dependence on unknown prior information, and improve data processing efficiency.
In 2006 and 2009, Ikelle proposed a virtual seismic event construction method to predict internal multiples [
16,
17]. This approach does not require prior information, such as velocity models or subsurface structures, enabling accurate internal multiple predictions. However, this method requires sequential extraction of primaries to construct internal multiples generated by relevant layers, which increase the computational complexity. In 2013, Wu et al. utilized a multi-channel L
1 norm adaptive matching algorithm to suppress multiples, considering the differences in amplitude and phase between virtual event method-predicted internal multiples and actual internal multiples, achieving favorable results in model data [
18]. In 2018, Liu et al. introduced an adaptive virtual event method to address the excessive dependency on the matching algorithm in the traditional virtual event method, enhancing internal multiple suppression capabilities [
19]. Although researchers have improved and modified the virtual event method to address various issues, the traditional virtual event method still requires sequential internal multiple constructions on subsurface interfaces, incurring high computational costs for seismic data derived from complex structures and making primary extraction difficult. To overcome this challenge, we propose an iterative virtual event method. Compared to the traditional virtual event method, the iterative virtual event technique incorporates an iterative calculation process, mitigates the influence of false frequencies generated by data space convolution on predicted multiple models, and significantly improves computational efficiency. When applied to actual onshore post-stack seismic data from the southwest depression of the Tarim Basin, our method effectively suppresses internal multiples that are challenging to eliminate in pre-stack traces, emphasizes the energy of deep effective waves, and further enhances the quality of seismic profiles.
The paper is organized as follows. After the introduction, we provide an overview of the study region. Subsequently, we present the theory of virtual events for internal multiple eliminations, followed by synthetic and field data experiments that validate the effectiveness of our approach. Finally, we draw conclusions based on detailed analyses and discussions.
2. The Overview of the Study Region
The Southwest Depression of the Tarim Basin, situated in the southwestern part of the basin, is bordered by the South Tianshan orogenic belt to the north and the West Kunlun orogenic belt to the south. Encompassing an area of approximately 100 × 100 km
2, it is the largest sub-basin within the Tarim Basin and a favorable area for oil and gas exploration in Cambrian Ordovician carbonate rocks. The traps formed during the Himalayan movement predominantly consist of effective hydrocarbon source rocks distributed in the Bachu fault uplift and Maigaiti slope area (Qu et al., 2000) [
20]. After extensive research, a consensus has emerged concerning the Cambrian subsalt of the Maigaiti slope: the Maigaiti slope contains Ordovician Carboniferous oil and gas sources, while drilling through the Cambrian south of Bachu revealed no Cambrian source rock, indicating the development of Cambrian subsalt source rock in the slope area. The Middle-Lower Cambrian serves as an important reservoir and regional cap rock for the area’s deep carbonate rocks and is a significant horizon for hydrocarbon source rock development (Cui et al., 2017; Zhang et al., 2015) [
21,
22]. Its seismic sequence features three strong reflective interfaces within the Middle Cambrian. The energy at the bottom boundary of the Lower Cambrian is not strong, whereas the bottom of the upper salt section displays a continuous and stable strong wave crest. A noticeable phase change occurs at the bottom of the lower salt section, and the sedimentary pattern from the Early Cambrian to Middle Cambrian exhibits significant changes.
The three strong reflective interfaces in the Middle Cambrian provide favorable conditions for internal multiple formations. However, internal multiples can distort seismic waveforms, decrease the signal-to-noise ratio, and affect effective wave identification, thereby complicating seismic processing, diminishing seismic imaging authenticity and reliability, influencing hydrocarbon source rock studies, and increasing the difficulty of oil and gas exploration and development in the northwest region. To better understand and identify the internal multiples in the study area, well-seismic calibration was performed using existing seismic and logging information, as illustrated in
Figure 1. The well-seismic calibration results reveal that the bottom interface of the upper salt section at the Middle Cambrian is stable. The lithology of the lower salt section of the Middle Cambrian comprises pure gypsum salt rock with a strong wave crest at the bottom, while the lower salt section contains gypsum-bearing dolomite and exhibits a trough at the bottom boundary. The bottom boundary of the Lower Cambrian primarily appears as a trough, and the time thickness of the Lower Cambrian varies significantly across different drilling wells, as does the mismatch between well synthetic records and seismic records. Subsequently, velocity spectra were generated for multiple identification and confirmation, as depicted in
Figure 2.
Figure 2 displays the velocity spectra of two distinct CMP gathers, revealing numerous low-velocity energy clusters. Following NMO correction of CMP gathers using primary velocity, downward-bending events indicate the presence of a large number of multiples in the original seismic data, necessitating multiple suppression processing.
Considering the characteristics of multiples in the study area, we initially employed the Radon transform to suppress internal multiples in the pre-stack data. However, the Radon transform cannot entirely eliminate all internal multiples present in the data, and this method can only suppress some long-period multiples. Multiples interference still exists under the Cambrian salt, as shown by the red circle in the velocity spectrum in
Figure 3. Additionally, there are significant differences in energy and frequency of multiples in the horizontal direction. It is challenging to suppress short-path internal multiples in the pre-stack data, necessitating the selection of an appropriate multiple suppression method from the post-stack data for this purpose. The iterative method for suppressing internal multiples based on virtual events is capable of accurately predicting internal multiples without requiring prior information, such as velocity models and underground structures, while maintaining high computational efficiency. Consequently, we opt for the iterative method to suppress internal multiples in the seismic data for the study area.
3. Internal Multiple Elimination Using the Theory of Virtual Events
The iterative virtual event internal multiple suppression approach extends the traditional virtual event internal multiple suppression method and avoids constructing internal multiple from every layer, which greatly improves computational efficiency [
23,
24]. Seismic virtual events are not directly recorded in standard seismic data acquisition, but their existence can help us construct internal multiples with scattering points at the sea surface [
16,
17]. Therefore, the total seismic wave field data
received from the surface can be used to successfully construct internal multiples. The process of predicting internal multiples mainly includes two parts. For the internal multiples related to the first underground interface, it is necessary to first extract the primary event
of
ith interface and the primary
generated by the interface below this interface. Then, correlating the primary
generated by other interfaces below this interface with the time reversal
of the primary event
to build the virtual event
of this interface (
Figure 4a). Finally, the convolution operation is carried out by using the constructed virtual event
and the primary
generated by the interface below this interface to construct the internal multiples about this interface (
Figure 4b). Namely,
where
represents any point on the surface,
is the receiver point,
is the source,
represents complex conjugation, and
represents angular frequency. The subscripts represent the number of layers, namely, the lower subscript
i in Formulas (1) and (2) represents the
ith layer.
represents the constructed post-stack virtual events.
is the internal multiple of the constructed post-stack data.
Figure 4 shows the construction process of virtual events and related internal multiples; constructing post-stack virtual events
for the correlation of delay wave
and lead wave
. Then, the constructed virtual events and primary convolution is used to obtain the internal multiples of the post-stack-related layers
, where the red star represents the convolution operation.
Assuming that the median of seismic data does not include surface-related multiples but only internal multiples, the effective wave after internal multiples suppression is as follows:
where
is the data of internal multiple suppression and
is the original seismic data. According to the estimation of multiple scattering by iterative inversion [
25,
26], the multiple can be written as follows:
where
is the surface factor. Substituting Formula (4) into Formula (3):
In terms of the Neumann series, Formula (5) can be written as follows:
where
n represents the number of iterations and the expression of surface factor
is as follows:
In order to obtain more abundant internal multiple information, the multiple iteration theory can be used for multiple models. The multiples after iteration
n + 1th can be expressed as follows:
The multiples after
n iterations are expressed as follows:
Bringing Formula (7) into Formula (9), the internal multiple model after
n iterations is as follows:
where
T represents matrix transpose.
The primary focus of this article is on the suppression of internal multiples in post-stack profiles. In the post-stack profile, firstly, we can employ the tracing of horizons method to obtain the primary events. Secondly, by shifting the window up and down, we can acquire the primary data volume and other primary data volume below this primary used to construct the virtual events. Subsequently, we input the seismic data obtained, the iterative virtual event internal multiple suppression technique initially utilizes the virtual events internal multiple methods, such as Formulas (1) and (2), to construct the internal multiple models for the corresponding layer. Finally, through iterative updates, such as Formula (10), a more comprehensive multiple model is obtained, ultimately leading to the suppression of the multiple results. This approach not only enhances the computational efficiency, accuracy, and applicability but also reduces the requirements for the original input data, thereby improving imaging accuracy, which can make it widely applicable in the processing of field data.
5. Field Example
Per the interpreter’s request, the primary focus is on suppressing internal multiples below 6 s to observe the original seismic profile more clearly. The seismic profile below 4 s is shown in
Figure 9. From
Figure 9a, it can be seen that the events below 6 s are similar in wave shape to the strong energy primaries, and the energy and frequency of these events are significantly different in the horizontal direction. Based on the relevant characteristics of multiples and multiples identification methods, these events are determined as internal multiples generated by the three sets of strong reflection interfaces of Middle Cambrian. Proof by facts, these internal multiples are challenging to suppress in pre-stack data using the filtering methods, such as the Radon transform method. Therefore, on top of the pre-stack suppression of multiples, we suppressed those internal multiples that are difficult to suppress in post-stack data using the iterative virtual events internal multiples suppression method.
Figure 9 shows the profile comparison before and after the internal multiple suppression.
Figure 9a displays the partial stack profile before internal multiple suppression, while
Figure 9b shows the corresponding stack profile after the internal multiple suppression. As seen in
Figure 9b, the internal multiples below 6 s are mostly suppressed, and the horizontal difference in energy and frequency of the events below 6 s is smaller after the internal multiples are suppressed, particularly at the position indicated by the red arrow. To more clearly observe the effect of internal multiple suppression before and after, the part in the yellow box in
Figure 9 is enlarged, as shown in
Figure 10. In
Figure 10, the event, which is similar to the upper strong reflection waveform with the opposite phase, periodic appearance, and significant lateral variation, is essentially suppressed. This is consistent with the analysis and identification results of multiples discussed earlier (
Figure 2 and
Figure 3). The results demonstrate that the iterative virtual event multiple suppression method can effectively suppress the internal multiples in field data.
To further analyze the suppression effect and accuracy of internal multiples, some sections before and after internal multiple suppression are selected for spectrum analysis, as shown in
Figure 11.
Figure 11a is the partial superimposed section before internal multiple suppression,
Figure 11b is the partial superimposed section after internal multiple suppression, and
Figure 11c shows the frequency distribution curve before and after internal multiple suppression. It can be seen from the figure that after suppressing internal multiples, the frequency spectrum is significantly broadened, and the spectral energy of the seismic profile is increased. This indicates that while using the iterative virtual event technology to suppress the internal multiples, the effective wave is protected, and the resolution of seismic data is improved. The method can provide reliable seismic data for subsequent studies of hydrocarbon source rocks.