## 2 Method
### Internal Multiple Prediction Using High-Order Born Modeling
The weighting matrix is set to a value of 1 at the position of the wave we need and a minimum value at the position of the wave we do not need. Where the value of \(\mathbf{W}_{pri}\) is calculated as 1, the useful content is the primary and the internal multiples could be suppressed; where the value of \(\mathbf{W}_{mul}\) is calculated as 1, the useful content are internal multiples and the primary could be suppressed. By building a division formula, we can set the values of \(\mathbf{W}_{pri}\) and \(\mathbf{W}_{mul}\) as we want. Both LSRTM and the multiple prediction method are also based on the Born modeling theory, so they can use each other's parameters. Thus, the calculation of the scatter wavefield requires LSRTM to provide an accurate reflectivity \(\mathbf{m}\). However, There may be a difference in phase between actual reflectivity and the reflectivity \(\mathbf{m}\) provided by LSRTM. Thus, the internal multiples we compute may differ from the observed data in the phase and position. Using the computed weighting matrix to suppress them may not completely suppress the internal multiples. Thus, we use the Hilbert transforms in order to use the low frequency component of the predicted internal multiples and observed data to compute \(\mathbf{W}\). The mismatch of the predicted internal multiples occurs mainly in the high frequency part, so introducing the Hilbert transform can reduce the influence of the mismatch. Then Equations (19) and (20) can be changed as:  
\[W_{pri}(x,t)= \begin{cases}\alpha|\mathrm{H}[d_{obsn}(x,t)]|/|\mathrm{H}[d_{ cmn}(x,t)]|\left|\mathrm{H}[d_{cmn}(x,t)]\right|>\alpha|\mathrm{H}[d_{obsn}(x,t)]| |\\ 1\left|\mathrm{H}[d_{cmn}(x,t)]\right|\leq\alpha|\mathrm{H}[d_{obsn}(x,t)]| |\end{cases}  \tag{21}\]  
\[W_{mul}(x,t)= \begin{cases}\beta|\mathrm{H}[d_{cmn}(x,t)]|/|\mathrm{H}[d_{obsn} (x,t)]|\left|\mathrm{H}[d_{obsn}(x,t)]\right|>\beta|\mathrm{H}[d_{cmn}(x,t)]| |\\ 1\left|\mathrm{H}[d_{obsn}(x,t)]\right|\leq\beta|\mathrm{H}[d_{cmn}(x,t)]| \end{cases}  \tag{22}\]  
Because the cross-correlation of the nth-order internal multiples in the forward- propagated wavefield and the backward-propagated wavefield of the nth-order internal multiples can create two correct images. Thus, the equation of the gradient should contain two terms related to the internal multiple. According to the second part and third part of equation \(\mathbf{A}_{0}\mathbf{U}_{s,2}=\mathbf{\Psi}_{1}\mathbf{B}\mathbf{U}_{s,1, \lambda}\), \(\mathbf{A}_{0}\mathbf{\Phi}=\mathbf{\Psi}_{2}\mathbf{B}\mathbf{U}_{s,2,\gamma}\), we can get the gradient:  
\[ \begin{split}&\mathrm{g}(\mathbf{x})\\ &=\int\limits_{0}^{\mathrm{max}}u_{0}^{*}(\mathbf{x},t,\mathbf{x}_{ \gamma})\frac{\partial^{2}}{\partial t^{2}}u_{0}(\mathbf{x},t,\mathbf{x}_{s}) dt\\ &+\int\limits_{0}^{\mathrm{max}}u_{1}^{*}(\mathbf{x},t,\mathbf{x}_{ \gamma})\frac{\partial^{2}}{\partial t^{2}}u_{1,\lambda}(\mathbf{x},t,\mathbf{ x}_{s})dt\\ &+\int\limits_{0}^{\mathrm{max}}u_{1}^{*}(\mathbf{x},t,\mathbf{x}_{ \gamma})\frac{\partial^{2}}{\partial t^{2}}u_{2,\gamma}(\mathbf{x},t,\mathbf{ x}_{s})dt\end{split}  \tag{23}\]  
where, \(u_{1,\lambda}(\mathbf{x},t,\mathbf{x}_{s})\) is the upward of the primary and \(u_{2,\gamma}(\mathbf{x},t,\mathbf{x}_{s})\) is the downward of the first-order internal multiple. To realize it, we have:\[(\frac{1}{v_{0}(\mathbf{x})^{2}}\frac{\partial^{2}}{\partial t^{2}}-\
abla^{2})u_{0 }^{*}(\mathbf{x},t,\mathbf{x}_{s})=W_{pri}W_{pri}[d_{abs}(\mathbf{x},T-t)-d_{cal} (\mathbf{x},T-t)] \tag{24}\]  
\[(\frac{1}{v_{0}(\mathbf{x})^{2}}\frac{\partial^{2}}{\partial t^{2}}-\
abla^{2})u _{1}^{*}(\mathbf{x},t,\mathbf{x}_{s})=W_{mut}W_{mut}[d_{abs}(\mathbf{x},T-t)-d_ {cal}(\mathbf{x},T-t)] \tag{25}\]  
Thus, we can conduct the LSRTM process using the weighting matrix to realize the internal multiples, as Figure 2 shows.