**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 km2, 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.

**Figure 1.** Well A synthetic record calibration and well seismic profile.

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.

**Figure 2.** Velocity spectrum of the different original CMP trace sets; (**a**) velocity spectrum of CMP1; (**b**) velocity spectrum of CMP2.

**Figure 3.** Velocity spectrum of the different CMP gathers after multiple suppression by the Radon transform (corresponds to Figure 2); (**a**) velocity spectrum of CMP1; (**b**) velocity spectrum of CMP2.

### **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 *P*(*xR*, *xS*, *ω*) 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 *Pi*(*xR*, *xS*, *ω*) of *i*th

interface and the primary <sup>−</sup> *Pi*(*xR*, *xS*, *ω*) generated by the interface below this interface. Then, correlating the primary <sup>−</sup> *Pi*(*xR*, *xS*, *ω*) generated by other interfaces below this interface with the time reversal *P<sup>H</sup> <sup>i</sup>* (*xR*, *xS*, *ω*) of the primary event *Pi*(*xR*, *xS*, *ω*) to build the virtual event *Piv*(*xR*, *xS*, *ω*) of this interface (Figure 4a). Finally, the convolution operation is carried

out by using the constructed virtual event *Piv*(*xR*, *xS*, *<sup>ω</sup>*) and the primary <sup>−</sup> *Pi*(*xR*, *xS*, *ω*) generated by the interface below this interface to construct the internal multiples about this interface (Figure 4b). Namely,

$$P\_{iv}(\mathbf{x}\_{\mathbb{R}}, \mathbf{x}\_{\mathbb{S}^\star} \boldsymbol{\omega}) = \sum\_{\mathbf{x}} \bar{P}\_i(\mathbf{x}, \mathbf{x}\_{\mathbb{S}^\star} \boldsymbol{\omega}) P\_i^H(\mathbf{x}\_{\mathbb{R}}, \mathbf{x}\_{\mathbb{R}} \boldsymbol{\omega}) \tag{1}$$

$$M\_i(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) = \sum\_{\mathbf{x}} P\_{\text{iv}}(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}, \omega) \bar{P}\_i(\mathbf{x}, \mathbf{x}\_{\mathcal{S}}, \omega) \tag{2}$$

where *x* represents any point on the surface, *R* is the receiver point, *S* is the source, *H* 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 *i*th layer. *Piv*(*xR*, *xS*, *ω*) represents the constructed post-stack virtual events. *Mi*(*xR*, *xS*, *ω*) is the internal multiple of the constructed post-stack data.

**Figure 4.** Construction of the virtual events and the internal multiples; (**a**) construction of the virtual events; (**b**) construction of the internal multiples. There pentagram represents convolution, the black line and black arrow represent the actual detectable seismic wave paths, red line and red arrow is used to assist in constructing virtual events and internal multiples.

Figure 4 shows the construction process of virtual events and related internal multiples; constructing post-stack virtual events *Pv*(*xR*, *xS*, *<sup>ω</sup>*) for the correlation of delay wave <sup>−</sup> *P*(*x*, *xS*, *ω*) and lead wave *PH*(*xR*, *x*, *ω*). Then, the constructed virtual events and primary convolution is used to obtain the internal multiples of the post-stack-related layers *M*(*xR*, *xS*, *ω*), 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:

$$P\_0(\mathbf{x}\_{\mathbb{R}\prime}\mathbf{x}\_{\mathbb{S}\prime}\omega) = P(\mathbf{x}\_{\mathbb{R}\prime}\mathbf{x}\_{\mathbb{S}\prime}\omega) - M(\mathbf{x}\_{\mathbb{R}\prime}\mathbf{x}\_{\mathbb{S}\prime}\omega),\tag{3}$$

where *P*0(*xR*, *xS*, *ω*) is the data of internal multiple suppression and *P*(*xR*, *xS*, *ω*) is the original seismic data. According to the estimation of multiple scattering by iterative inversion [25,26], the multiple can be written as follows:

$$M(\mathbf{x}\_{\mathbb{R}}, \mathbf{x}\_{\mathbb{S}\prime}\omega) = P\_{\mathbb{0}}(\mathbf{x}\_{\mathbb{R}}, \mathbf{x}\_{\mathbb{S}\prime}\omega) \{ \_{0}AP(\mathbf{x}\_{\mathbb{S}\prime}\mathbf{x}\_{\mathbb{R}\prime}\omega) \tag{4}$$

where *A* is the surface factor. Substituting Formula (4) into Formula (3):

$$P\_0(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) = P(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) - P\_0(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) \, \_0AP\_0(\mathbf{x}\_{\mathcal{S}}, \mathbf{x}\_{\mathcal{R}}, \omega) \tag{5}$$

In terms of the Neumann series, Formula (5) can be written as follows:

$$\{P\_0(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega)\}^{(n)} = P(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) - \{P\_0(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega)\}\_{0}^{(n-1)} A^{(n)} P(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega),\tag{6}$$

where *n* represents the number of iterations and the expression of surface factor *A* is as follows:

$$\begin{array}{l} A^{(n)} = \left( \left\{ P\_0(\mathbf{x}\_{R'}, \mathbf{x}\_{S'} \boldsymbol{\omega}) \right\} \right)^{(n-1)}\_{0} \left( P(\mathbf{x}\_{R'}, \mathbf{x}\_{S'} \boldsymbol{\omega}) \right) \\ \quad - \left\{ P\_0(\mathbf{x}\_{R'}, \mathbf{x}\_{S'} \boldsymbol{\omega}) \right\}^{\mathrm{tr}} \left( P(\mathbf{x}\_{R'}, \mathbf{x}\_{S'} \boldsymbol{\omega}) \right)^{-1} \end{array} \tag{7}$$

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:

$$\{M(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega)\}^{(n+1)} = \{P\_0(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega)\}\_0^{(n)} A^{(n)} P(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) \tag{8}$$

The multiples after *n* iterations are expressed as follows:

$$\begin{aligned} \left\{ M(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) \right\}^{(n)} &= \left\{ P\_{0}(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) \right\}^{(n-1)}\_{0} \left[ \left\{ P\_{0}(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) \right\} \right]^{(n-2)}\_{0} \left[ \right. \end{aligned} \tag{9}$$
 
$$\times \left( P(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) - \left\{ P\_{0}(\mathbf{x}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{S}}, \omega) \right\} \right)^{(n-1)} \tag{9}$$

Bringing Formula (7) into Formula (9), the internal multiple model after *n* iterations is as follows:

$$\begin{array}{rcl}\{M(\mathbf{x}\_{\mathcal{R}},\mathbf{x}\_{\mathcal{S}},\omega)\}^{(n)} &=& \{P\_{0}(\mathbf{x}\_{\mathcal{R}},\mathbf{x}\_{\mathcal{S}},\omega)\}\_{0}^{(n-1)} (\{P\_{0}(\mathbf{x}\_{\mathcal{R}},\mathbf{x}\_{\mathcal{S}},\omega)\}\_{0}^{(n-1)})^{T} \\ & \times [\{P\_{0}(\mathbf{x}\_{\mathcal{R}},\mathbf{x}\_{\mathcal{S}},\omega)\}\_{0}^{(n-2)} (\{P\_{0}(\mathbf{x}\_{\mathcal{R}},\mathbf{x}\_{\mathcal{S}},\omega)\}\_{0}^{(n-2)})^{T}]^{-1} \\ & \times (P(\mathbf{x}\_{\mathcal{R}},\mathbf{x}\_{\mathcal{S}},\omega) - \{P\_{0}(\mathbf{x}\_{\mathcal{R}},\mathbf{x}\_{\mathcal{S}},\omega)\}\_{0}^{(n-1)})^{T}) \end{array} \tag{10}$$

where *T* represents matrix transpose.

The primary focus of this article is on the suppression of internal multiples in poststack 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.
