## 2 Method
### Internal Multiple Prediction Using High-Order Born Modeling  
Here, the internal multiple prediction method we use is the high-order Born modeling method. This method obtains the predicted internal multiples in a forwarding process. The high-resolution velocity model and reflectivity model that are needed in the forwarding process can be provided and updated from LSRTM. This method provides separated forward-propagated internal multiples for LSRTM imaging. In the theory of scattering series, we first get the wave equation of the actual velocity model and the reference velocity model:  
\[A\mathbf{U}=-\boldsymbol{\updelta}\]  
\[A_{0}\mathbf{U}_{0}=-\boldsymbol{\updelta} \tag{5}\]  
where \(\mathbf{U}\) and \(\mathbf{U}_{0}\) is wavefield of the actual model and the reference model, \(A\) and \(A_{0}\) are the differential operators of the actual model and the reference model and satisfy:  
\[A=\frac{1}{v(\mathbf{x})^{2}}\frac{\partial^{2}}{\partial t^{2}}+\
abla^{2}\]  
\[A_{0}=\frac{1}{v_{0}(\mathbf{x})^{2}}\frac{\partial^{2}}{\partial t^{2}}+ \
abla^{2} \tag{6}\]  
The \(\mathbf{U}\) can be expressed as the wavefield of the reference reflectivity adding the scatter wavefield:  
\[\mathbf{U}=\mathbf{U}_{0}+\mathbf{U}_{s} \tag{7}\]  
where \(\mathbf{U}_{s}\) is the scattering wavefield. Let the two equations of Equation (5) be subtracted, and let \(B=A-A_{0}=m(x,z)\frac{\partial^{2}}{\partial t^{2}}\), where \(m(\mathbf{x})=\frac{1}{v^{2}(\mathbf{x})}-\frac{1}{v_{0}^{2}(\mathbf{x})}\) is the reflectivity model following, we can have:  
\[A_{0}\mathbf{U}_{s}=B\mathbf{U} \tag{8}\]  
Equation (8) describes how the scatter wavefield can be computed, However, We cannot get the value of \(\mathbf{U}\) in real cases, so we use \(\mathbf{U}_{0}\) as an approximate replacement of \(\mathbf{U}\). We therefore have:  
\[A_{0}\mathbf{U}_{s}=B\mathbf{U}_{0} \tag{9}\]  
We use Equation (9) to compute \(\mathbf{U}_{s}\). By continuously substituting \(\mathbf{U}_{s}\) into Equation (7) and substituting Equation (7) into Equation (8) to compute \(\mathbf{U}_{s}\), we can get the internal multiple predicted by the high-order Born modeling method. The derivation is shown in Appendix A and the method is expressed as follows:  
\[\left\{ \begin{array}{l}A_{0}\mathbf{U}_{s,1}=\boldsymbol{\upgamma}_{0}B \mathbf{U}_{0,\gamma}\\ A_{0}\mathbf{U}_{s,2}=\boldsymbol{\upgamma}_{1}B\mathbf{U}_{s,1,\lambda}\\ A_{0}\boldsymbol{\upphi}=\boldsymbol{\upgamma}_{2}B\mathbf{U}_{s,2,\gamma} \end{array} \right. \tag{10}\]  
where \(\mathbf{U}_{s}\) is the scatter wavefield and can be written as the series \(\mathbf{U}_{s}=\mathbf{U}_{s,1}+\mathbf{U}_{s,2}+\mathbf{U}_{s,3}+\ldots\), \(\mathbf{U}_{s,l}\) is the portion of \(\mathbf{U}_{s}\) that is ith order. The subscripts \(\gamma\) and \(\lambda\) represent the down- and upgoing components, and \(\boldsymbol{\upPhi}\) is the first-order internal multiple we need to calculate.\(\boldsymbol{\upgamma}\) is the zeroing matrix and is expressed as:  
primary or multiple by given different weighting according to the observed data and the predicted internal multiple. Here, we set:  
\[W_{pri}(x,t)= \begin{cases}\alpha|d_{obsn}(x,t)/d_{cmn}(x,t)|\left|d_{cmn}(x,t) \right|>\alpha|d_{obsn}(x,t)|\\ 1\left|d_{cmn}(x,t)\right|\leq\alpha|d_{obsn}(x,t)|\end{cases}  \tag{19}\]  
\[W_{mul}(x,t)= \begin{cases}\beta|d_{cmn}(x,t)/d_{obsn}(x,t)|\left|d_{obsn}(x,t )\right|>\beta|d_{cmn}(x,t)|\\ 1\left|d_{obsn}(x,t)\right|\leq\beta|d_{cmn}(x,t)|\end{cases}  \tag{20}\]  
where \(W_{pri}(x,t)\) and \(W_{mul}(x,t)\) represent the value of \(\mathbf{W}_{pri}\) and \(\mathbf{W}_{mul}\) when the coordinate is \((x,t)\), \(d_{cmn}(x,t)\) represents the normalized predicted internal multiples, \(d_{obsn}(x,t)\) represents the normalized observed data, and \(0\leq\alpha\leq 1\), \(0\leq\beta\leq 1\) is a threshold parameter. Through introducing the weighting matrix, we can select the wave we need.