## 4 Discussion  
In Equation (11b), we can see that the summation of the time derivative and the spatial gradient images had an extra term of \(\cos^{2}\theta\) compared to the time derivative image. When we designed the small-angle imaging condition in Equation (13), we mentioned that the exponential term approximated to \(e^{-\text{asin}^{2}\theta}\). Since this exponential weighting was applied at each time step, we can consider the proposed imaging condition as the time derivative imaging condition weighted by the term \(\cos^{2}\theta e^{-\text{asin}^{2}\theta}\). Figure 10 shows the analysis of this term by choosing different parameter \(\alpha\), compared to the term \(\cos^{2}\theta\). In the figure, we can see that a larger \(\alpha\) attenuated the large reflections more rapidly than a small \(\alpha\); on the other hand, the attenuation curve was smoother for a small \(\alpha\), which generated less artifacts than a large \(\alpha\). In all our examples in the previous section, we empirically chose \(\alpha=5\) based on a compromise between the preservation of accuracy, and artifacts reduction. Figure 11 shows the sensitivity analysis of the inversion results using different \(\alpha\) (\(\alpha=0,5,10,20\)) for the layered model (Figure 4). In the figure, we can see that for \(\alpha=0\), the velocity and impedance perturbation cannot be distinguished, while for a larger \(\alpha\) (\(\alpha=10\) and \(\alpha=20\)), the shallow part is problematic. In future studies, we may consider the functions of different forms to attenuate large reflections while still preserving the imaging quality. Figure 10: Analysis of different angle-dependent functions. Figure 9: Comparison of filtered true impedance perturbation (red) and estimated results (blue) at three locations: (**a**) X = 1900; (**b**) X = 2200; and (**c**) X = 3000. In all our numerical examples, broadband sources were used for synthetic data generation. Broad bandwidth played an important role in seismic inversion. Increasing the bandwidth at the high-frequency end improved seismic resolution, and adding more low frequencies in the data reduced the side lobes of the seismic wavelet [18]. Figure 12 demonstrates the impact of low frequencies on the inversion result. For the Sigsbee2b model, we applied a band-pass filter (\(f_{1}=3\) Hz, \(f_{2}=5\) Hz, \(f_{3}=56\) Hz, \(f_{4}=60\) Hz) to remove low frequencies in the input data, and then estimated the impedance perturbation using our proposed method. The estimated impedance perturbations with and without low frequencies are shown in Figures 12a and 12b, respectively. Figure 12a is the same as Figure 6c, and is repeated here for comparison. Figure 12c shows the trace comparison of Figure 12a,b at the horizontal location of 3500 m. The results demonstrate that low frequencies in the data affect the long-wavelength structures in the image. Low frequencies are crucial for both velocity and impedance inversions [8], and retaining low frequencies in the recorded data depends on broadband seismic acquisition and data processing. Figure 11: Estimated relative impedance perturbation for the layered model using different \(\alpha\) in the weighting function in Equation (14): (**a**) \(\alpha=0\); (**b**) \(\alpha=5\); (**c**) \(\alpha=10\); and (**d**) \(\alpha=20\). Figure 12: Comparison of estimated impedance perturbation with (**a**) and without (**b**) low frequencies in input data. (**c**) Trace comparison of estimated impedance perturbation with (blue) and without (green) low frequencies at the horizontal location of 3500 m.  
As in conventional RTM, the proposed method assumes primaries only in the data, so de-multiple is a necessary preprocessing step in the case of strong multiples in the data. Even though our method was derived from the framework of amplitude-preserving RTM, the inversion was based on the Born approximation (Equation (4)). Therefore, transmission losses due to the overburden effects were also a challenge for our proposed method.