## Appendix A  
The scattering wavefield can be written as the series:  
\[\mathbf{U}_{s}=\mathbf{U}_{s,1}+\mathbf{U}_{s,2}+\mathbf{U}_{s,3}+\ldots \tag{11}\]  
\(\mathbf{U}_{s,i}\) is the portion of \(\mathbf{U}_{s}\) that is \(i\)th order. And we can write Equation (9) as:  
\[A_{0}\mathbf{U}_{s,1}=B\mathbf{U}_{0} \tag{12}\]  
Use \(\mathbf{U}_{s,1}\) as an approximate replacement of \(\mathbf{U}_{s}\). Substitute \(\mathbf{U}_{s}\) into Equation (7) and substitute Equation (7) into Equation (8), and we can have:  
\[A_{0}\mathbf{U}_{s}=B(\mathbf{U}_{0}+\mathbf{U}_{s,1})=A_{0}\mathbf{U}_{s,1}+B \mathbf{U}_{s,1} \tag{13}\]  
So we can get:  
\[A_{0}\mathbf{U}_{s,2}=B\mathbf{U}_{s,1} \tag{14}\]  
Use \(\mathbf{U}_{s,1}\) and \(\mathbf{U}_{s,2}\) approximate replacement of \(\mathbf{U}_{s}\). Substitute \(\mathbf{U}_{s}\) into Equation (7) and substitute Equation (7) into Equation (8) again, we can have:  
\[A_{0}\mathbf{U}_{s,3}=B\mathbf{U}_{s,2} \tag{15}\]  
The \(\mathbf{U}_{s,3}\) represents the wavefield experiencing scattering for three times. The first-order internal multiple \(\mathbf{\Phi}\) is also the wavefield that experienced scattering for three times and contained in the \(\mathbf{U}_{s,3}\). The propagation of the first-order internal multiple is down-up-down, so we can get the first order internal multiple \(\mathbf{\Phi}\) when \(\mathbf{U}_{s,3}\) meets the depth condition of \(z_{1}<z_{2},z_{3}<z_{2}\), where \(z_{1},z_{2},z_{3}\) is the depth of the first, second, and third scatter point. Here we do wavefield decomposition to meet the depth condition. Then we can get:  
\[\left\{ \begin{aligned} \mathbf{A}_{0}\mathbf{U}_{s,1}= \mathbf{BU}_{0,\gamma}\\ \mathbf{A}_{0}\mathbf{U}_{s,2}=\mathbf{BU}_{s,1,\lambda}\\ \mathbf{A}_{0}\mathbf{\Phi}=\mathbf{BU}_{s,2,\gamma}\end{aligned} \right. \tag{16}\]  
The subscripts \(\gamma\) and \(\lambda\) indicate that the propagation of the wavefield is down and up. Subscripts 1 and 2 indicate that the scatter wavefields is first and second order. In Figure 1\(u\) is the wavefield, and the subscripts \(r\) indicate that the propagation of the wavefield is rightward, and \(u^{n+1}\) is generated by \(u^{n}\) according to Equation (16). Through Figure 2 we see that \(u^{n+1}_{r}\) makes the virtual source of \(u^{n}\) still exist in \(u^{n+1}\). In other words, there is overlap of \(u^{n}\) and \(u^{n+1}\). In this case, the depth condition can only reach \(z_{1}\leq z_{2},z_{3}\leq z_{2}\). To solve this problem, we need not only to split the up- and downgoing wavefields but also to introduce the zeroing matrix to remove the overlap of \(u^{n}\) and \(u^{n+1}\). Then we can get:  
\[\left\{ \begin{aligned} \mathbf{A}_{0}\mathbf{U}_{s,1}= \mathbf{\Psi}_{0}\mathbf{BU}_{0,\gamma}\\ \mathbf{A}_{0}\mathbf{U}_{s,2}=\mathbf{\Psi}_{1}\mathbf{BU}_{s,1, \lambda}\\ \mathbf{A}_{0}\mathbf{\Phi}=\mathbf{\Psi}_{2}\mathbf{BU}_{s,2, \gamma}\end{aligned} \right. \tag{17}\]  
* [PERSON] and [PERSON] (1992) [PERSON]; [PERSON] Suppression of multiple reflections using the Radon transform. _Geophysics_**1992**, _57_, 386-395. [CrossRef]
* [PERSON] (1998) [PERSON] Wave equation prediction and removal of interbed multiples. In Proceedings of the 68 th Annual International Meeting, SEG Expanded Abstracts, New Orleans, LA, USA, 13-18 September 1998; pp. 1527-1530. * [PERSON] and [PERSON] (2005) [PERSON]; [PERSON] Removal of internal multiples with the common- focus-point (CFP) approach: Part 2--Application strategies and data examples. _Geophysics_**2005**, _70_, V61-V72. [CrossRef]
* [PERSON] (2006) [PERSON] A construct of internal multiples form surface data only: The concept of virtual seismic events. _Geophysics_**2006**, _164_, 383-393. * [PERSON] (2009) [PERSON] Scattering diagrams in seismic imaging: More insights into the construction of virtual events and internal multiples. _J. Appl. Geophys._**2009**, _67_, 150-170. [CrossRef]
* [PERSON] et al.