## 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(x_{R},x_{S},\omega)\) 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 \(P_{i}(x_{R},x_{S},\omega)\) of \(i\)th interface and the primary \(\overset{\rightarrow}{P_{i}}(x_{R},x_{S},\omega)\) generated by the interface below this interface. Then,  
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. Figure 2: Velocity spectrum of the different original CMP trace sets; (**a**) velocity spectrum of CMP1; (**b**) velocity spectrum of CMP2. correlating the primary \(\overline{P_{l}}(x_{R},x_{S},\omega)\) generated by other interfaces below this interface with the time reversal \(P_{l}^{H}(x_{R},x_{S},\omega)\) of the primary event \(P_{l}(x_{R},x_{S},\omega)\) to build the virtual event \(P_{iv}(x_{R},x_{S},\omega)\) of this interface (Figure 4a). Finally, the convolution operation is carried out by using the constructed virtual event \(P_{iv}(x_{R},x_{S},\omega)\) and the primary \(\overline{P_{l}}(x_{R},x_{S},\omega)\) generated by the interface below this interface to construct the internal multiples about this interface (Figure 4b). Namely,  
\[P_{iv}(x_{R},x_{S},\omega)=\sum_{x}\overline{P_{l}}(x,x_{S},\omega)P_{l}^{H}( x_{R},x,\omega) \tag{1}\]  
\[M_{i}(x_{R},x_{S},\omega)=\sum_{x}P_{iv}(x_{R},x,\omega)\overline{P_{l}}(x,x_{ 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 \(\omega\) 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. \(P_{iv}(x_{R},x_{S},\omega)\) represents the constructed post-stack virtual events. \(M_{i}(x_{R},x_{S},\omega)\) 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 \(P_{v}(x_{R},x_{S},\omega)\) for the correlation of delay wave \(\overline{P}(x,x_{S},\omega)\) and lead wave \(P^{H}(x_{R},x,\omega)\). Then, the constructed virtual events and primary convolution is used to obtain the internal multiples of the post-stack-related layers \(M(x_{R},x_{S},\omega)\), 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}(x_{R},x_{S},\omega)=P(x_{R},x_{S},\omega)-M(x_{R},x_{S},\omega), \tag{3}\]  
where \(P_{0}(x_{R},x_{S},\omega)\) is the data of internal multiple suppression and \(P(x_{R},x_{S},\omega)\) 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(x_{R},x_{S},\omega)=P_{0}(x_{R},x_{S},\omega)\}_{0}AP(x_{S},x_{R},\omega) \tag{4}\]  
where \(A\) is the surface factor.