## 3 Results  
### Four-Layer Model: Finite Difference Data  
First, we used a simple four-layer model to verify the algorithm's capability (Figure 3). The velocities of the model from top to bottom were 1500, 2050, 2300, and 2600 m/s, and the densities were 1, 2, 3.5, and 2 g/cm\({}^{3}\). In this model, the third layer exhibited density inversion, forming two strong impedance interfaces above and below, leading to the development of internal multiples between the second and third layers. The parameters of this model do not directly correspond to the real lithology. The main objective was to use a simple four-layer model to visually demonstrate the interference artifacts of different-order multiples in imaging and the suppression effects of LSFWM. The model size was 701 \(\times\) 401, with horizontal and vertical sampling of 10 m and 5 m, respectively. Figure 6 displays the reflection data of this experiment, which have been modeled by an acoustic finite difference modeling scheme. In Figure 6a, seismic data for a common-shot gather are represented, where the seismic records were obtained by fixing the source (keeping a specific shot point location) and arranging receivers to record the seismic signals. In Figure 6b, a common-offset section is depicted, where multiple shots are considered, the shot points are not fixed, and the distance from the shot point to the receiver is constant. The difference between the two lies in the different selections of the entire seismic matrix. The black arrows indicate the primaries. The solid white arrows and the dashed white arrows represent some surface-related multiples and internal multiples, respectively. A total of 101 shots were evenly distributed between 0 and 7 km, with 701 traces per shot and a trace interval of 10 m. The source wavelet was a Ricker wavelet with a dominant frequency of 15 Hz, and both the source and receiver depths were at 5 m. The simulated data, obtained using free surface conditions, contain primaries, surface-related multiples, and internal multiples. The receiver sampling interval was 4 ms, and the total recording time was 3 s. Due to the practical impossibility of obtaining a true velocity model, the migration velocity model (Figure 7) was derived by applying Gaussian smoothing to the true velocity model (Figure 3a). Figure 6: Full-wavefield data of layer model computed with the finite difference method: (**a**) common\(-\)shot gather; (**b**) common\(-\)offset gather. The black arrow is primary wave, the solid white arrow is surface-related multiple, the dashed white arrow is internal multiple. The full-wavefield data were migrated by employing three imaging conditions (Equations (12), (13), and (15)) (Figure 8). By comparing the yellow rectangles in Figure 8a-c, where the black arrows represent the true reflectors, it is evident that the coherency axis resolution is higher in Figure 8b,c compared to Figure 8a. Additionally, the energy in Figure 8c is stronger than in Figure 8b. This observation suggests the deconvolution imaging condition has higher resolution compared to the cross-correlation imaging condition. In addition, the multi-order wavefield normalization condition preserves the high-resolution benefit of the deconvolution imaging condition while offering supplementary deep energy compensation. This is attributed to the fact that, for each frequency, the summation of all the shots for the downward wavefield can function as an auxiliary source illumination condition for the imaging, thus facilitating amplitude compensation in deeper regions. Figure 8: Migration results under three imaging conditions: (**a**) cross-correlation; (**b**) deconvolution; (**c**) global normalization. The black arrow is primary correct imaging, the yellow box indicates the area of interest. Figure 7: Migration velocity model (Gaussian smoothing of the true velocity model). Figure 9 shows the suppression results for multiple interferences in least-squares full-wavefield migration. Figure 9a displays the least-squares primary migration.