## 1 Introduction  
In seismic data, there are the primary and multiples. In most cases, the multiples cause artifacts in imaging. Thus, multiple removal methods are an important part of seismic imaging. However, multiples contain information about the subsurface structure, and the imaging of multiples can be a good complement to imaging. As a result, much research has focused on using multiples in imaging. To reduce the artifacts caused by multiples in imaging, many multiple removal methods have been proposed. Previously, multiple removal methods used the difference in periodicity and travel time between the primary and multiple waves; these methods included the normal moveout (NMO) method, common midpoint stack, Wiener filtering, predictive deconvolution, f-k filter, and Radon transformation [1; 2; 3; 4; 5; 6; 7; 8]. Currently, the main internal multiple suppressing methods are based on the wave equation. There are two effectively multiple removal methods that have been widely used. Internal multiple elimination (IME) methods predict internal multiples based on convolution and correlation [9; 10; 11; 12; 13]. This method isn't rely on the model and predicts internal multiples layer-by-layer. The inverse scattering series (ISS) internal multiple attenuation methods is derived from the Lippmann-Schwinger equation [14; 15; 16; 17; 18; 19; 20; 21; 22]. By computing the component which satisfied the depth condition of its third term in the wavenumber domain, the internal multiples can be predicted. The ISS can predict the multiple precisely but need several Fourier transforms which make its computational cost high. In addition, many other methods have been proposed in recent years. [PERSON] published a model-driven method in 2008 which predicts the internal multiples by using several extrapolations of the data [23]. Marchenko multiple elimination is derived from the coupled Marchenko equation and it is a model independent [24]. It can estimate the surface multiple and internal multiple, and use them in imaging at the same time, but it needs to be improved in complex structure. [PERSON] proposed a method called full wave migration (FWM) [25; 26; 27] which uses one-way wave equation to estimate the primary, surface multiples and internal multiples and uses them in imaging. Reverse time migration (RTM) was proposed by [PERSON] and developed thereafter [28; 29; 30; 31; 32]. Least squares migration (LSM) was first proposed by [PERSON] in 1984 [33]. [PERSON] et al. [34] proposed LSRTM by combining LSM and RTM in 2011. LSRTM is a method of obtaining the underground reflectivity distribution by solving the optimization problem. LSRTM therefore has more balanced amplitude and less crosstalk than RTM. Many studies have been proposed to remove multiples or use the multiples as useful information in LSRTM. [PERSON] and [PERSON] used both primaries and multiples to image in LSM [35]. [PERSON] et al. proposed the RTM of multiples (RTMM) method in 2011 [36; 37]. By cross-correlating the backward-propagated and forward-propagated recorded data on the surface, RTMM can realize the image of the primary and free-surface multiples. [PERSON] and [PERSON] proposed the LSRTM of multiples (LSRTMM) method based on RTMM in 2014 [38] to make the amplitude more balanced and to have crosstalk artifacts than RTMM. There are also applications in vertical seismic profile (VSP) ghosts and the viscoacoustic media of LSRTMM [39; 40]. Joint LSRTM (JLSRTM) was proposed by [PERSON] et al. in 2015 [41]. Through constructing a modified modeling generator, JLSRTM can jointly image the primary and free-surface multiples and increase the illumination of a complex structure and reduce the artifacts. [PERSON] et al. have applied JLSRTM to the plane wave domain and achieved a good effect [42]. [PERSON] et al. have applied JLSRTM to viscoacoustic media [43]. [PERSON] proposed LSRTM using controlled-order multiples in 2016 [44].