## 2 Methods
### Different-Order Multiple Crosstalk Artifacts and Illumination Analysis
Equation (5) formulates the following modeling scheme for the downward wavefield:  
\[D_{F}(z_{m})=\sum_{n=0}^{m-1}\mathbf{W}_{D}(z_{m},z_{n})\delta S(z_{n}), \tag{5}\]  
where \(D_{F}(z_{m})\) is the downward forward propagation wavefield at the depth level \(z_{m}\), \(\delta S(z_{n})\) represents the scattered wavefield at every grid point positioned at the depth level \(z_{n}\), matrix \(W_{D}(z_{m},z_{n})\) represents a cumulative downward propagation operator that brings the scattered wavefield related to the depth level \(z_{n}\) to the depth level \(z_{m}\), which is formulated as follows:  
\[\mathbf{W}_{D}(z_{m},z_{n})=\mathbf{W}_{D}(z_{m},z_{m-1})\mathbf{W}_{D}(z_{m-1},z_{m-2}) \times\cdots\times\mathbf{W}_{D}(z_{n+1},z_{n}). \tag{6}\]  
Equation (7) formulates the following modeling scheme for the upward wavefield:  
\[U_{F}(z_{m})=\sum_{n=m+1}^{N}\mathbf{W}_{U}(z_{m},z_{n})\delta S(z_{n}), \tag{7}\]  
where \(U_{F}(z_{m})\) is the upward forward propagation wavefield at the depth level \(z_{m}\), matrix \(\mathbf{W}_{U}(z_{m},z_{n})\) represents a cumulative upward propagation operator that brings the scattered wavefield related to the depth level \(z_{n}\) to the depth level \(z_{m}\), which is formulated as follows:  
\[\mathbf{W}_{U}(z_{m},z_{n})=\mathbf{W}_{U}(z_{m},z_{m+1})\mathbf{W}_{U}(z_{m+1},z_{m+2}) \times\cdots\times\mathbf{W}_{U}(z_{n-1},z_{n}) \tag{8}\]  
The scattering operator \(\delta S(z_{m})\) depends on the reflectivity operator \(\mathbf{R}(z_{m})\), the downward wavefield \(D_{F}(z_{m})\), and the upward wavefield \(U_{F}(z_{m})\),  
\[\delta S(z_{m})=\mathbf{R}^{\cup}(z_{m})D_{F}(z_{m})+\mathbf{R}^{\cap}(z_{m})U_{F}(z_{ m}), \tag{9}\]  
where \(\mathbf{R}^{\cup}\) and \(\mathbf{R}^{\cap}\) are the reflectivity matrices from above and below, \(\mathbf{R}^{\cup}=-\mathbf{R}^{\cap}\). The formula for calculating the \((i+1)\)th-order scattering wavefield \(M^{i}_{illu}\) two-way wavefield illumination is as follows:  
\[M^{i}_{illu}(z)=\sum_{\omega}\sum_{nshot}\Bigl{(}D^{i}_{F}(z,\omega)\Bigl{[}D ^{i}_{F}(z,\omega)\Bigr{]}^{*}+U^{i}_{F}(z,\omega)\Bigl{[}U^{i}_{F}(z,\omega) \Bigr{]}^{*}\Bigr{)}, \tag{10}\]  
where \(nshot\) is the number of shots. We use a simple four-layer model, and the velocity model and density model are shown in Figure 3. We then employed the FWMOD algorithm to generate primary waves and multiples of different orders. We performed migration imaging on each wavefield separately and obtained the crosstalk artifact and correct imaging results for each order's scattering (Figure 4). Figure 4 displays an image of the primary migration of full-wavefield data, as described in Equation (3). Figure 4 illustrates an image of the primary migration of primary data, as described in Equation (1). Figure 4 demonstrates the crosstalk artifact of the first-order surface-related multiple and first-order internal multiple, respectively, as shown in Equation (4) \(U^{1}_{B}D^{0^{*}}_{F}\). Figure 4 shows the correct imaging of the first-order surface-related multiple, corresponding to \(U^{1}_{B}D^{1^{*}}_{F}\) in Equation (4). Figure 4 shows the correct imaging of the second-order surface-related multiple, corresponding to \(U^{2}_{B}D^{2^{*}}_{F}\) in Equation (4). Figure 3: Four-layer model: (**a**) P-wave velocity model; (**b**) density model. Within the context of multiple imaging, both the surface-related multiples and internal multiples are incorporated during the source forward propagation process. The surface-related multiples and internal multiples can be regarded as new virtual sources, allowing for illumination analysis to be conducted on the primaries and different orders of multiples (Figure 5). The use of normalized amplitude is employed to represent the illumination intensity, achieved by dividing by the average energy of the first depth grid point illuminated by the primary wave. Figure 5a illustrates only primary wavefield illumination. Figure 5b shows the illumination of the primary wavefield and the internal multiple wavefield.