## 2 Methods
### Multi-Order Wavefield Global Normalization Imaging Condition  
The deconvolution imaging condition \(I_{deco}\) was proposed by [PERSON] [44].  
\[I_{deco}(x)=\sum_{\omega}\sum_{nshot}\frac{U_{B}(x,\omega)}{D_{F}(x,\omega)}, \tag{11}\]  
where \(U_{B}(x,\omega)\) denotes the backward-propagating upward wavefield at the receiver, and \(D_{F}(x,\omega)\) denotes the forward-propagating downward wavefield at the source.  
The cross-correlation imaging condition \(I_{corr}\) is as follows:  
\[I_{corr}(x)=\sum_{\omega}\sum_{nshot}U_{B}(x,\omega)D_{F}{}^{*}(x,\omega). \tag{12}\]  
Compared to the cross-correlation imaging condition, the deconvolution imaging condition has the instability of a zero denominator. Many methods have been proposed to improve this, and Equation (13) \(I_{sm-deco}\) is the well-known smoothing deconvolution imaging condition [16] of adding a smoothing window \(\langle\rangle\) with a small number \(\varepsilon\). This can ensure the stability of the algorithm, but the validity of its value is limited, to some extent, by the selection of the number \(\varepsilon\). In primary imaging, the source's forward wavefield only has a primary wavefield, while in multiple imaging, additional virtual source wavefields are involved, making the selection of an appropriate \(\varepsilon\) value more difficult.  
Figure 5: Normalized amplitude of two-way illumination of different order wavefield: (**a**) primary illumination; (**b**) primary and internal multiple illumination; (**c**) primary and surface-related multiple illumination; (**d**) full-wavefield (primary and surface-related multiple and internal multiple) illumination.  
\[I_{sm-deco}(x)=\sum_{\omega}\sum_{nshot}\frac{U_{B}(x,\omega)D_{F}{}^{*}(x,\omega) }{\langle U_{B}(x,\omega)D_{F}{}^{*}(x,\omega)\rangle+\varepsilon}. \tag{13}\]  
Equation (13) can be further improved to a global normalization imaging condition \(I_{gl}\) (Equation (14)). For each single-frequency component, the denominator term considers the forward wavefields of all the shot points, and its physical meaning is the energy distribution of all the shots in the subsurface space. The advantage of this imaging condition is that, since it calculates each single-frequency component, it retains the resolution-improving characteristic of the deconvolution imaging condition while providing additional illumination compensation, promoting the subsequent least-squares migration.  
\[I_{gl}(x)=\sum_{\omega}\frac{\sum\limits_{nshot}U_{B}(x,\omega)D_{F}{}^{*}(x, \omega)}{\sum\limits_{nshot}U_{B}(x,\omega)D_{F}{}^{*}(x,\omega)}. \tag{14}\]  
Considering that multiples can provide broader illumination, Equation (14) can be further improved by adding illumination compensation for multiple orders of surface-related multiples and internal multiples to the denominator term, resulting in a global-normalized imaging condition for multiple-order wavefields \(I_{mgl}\) (Equation (15)), which can be used for both multiple and primary imaging frameworks.  
\[\begin{array}{ll}I_{mgl}(x)&=\sum\limits_{\omega}\frac{\sum\limits_{nshot}U_{ B}(x,\omega)D_{F}^{*}(x,\omega)}{\sum\limits_{nshot}^{D}(x,\omega)D_{F}^{*}(x, \omega)}\\ &=\sum\limits_{\omega}\frac{\sum\limits_{nshot}U_{B}(x,\omega)\left(D_{F}^{* }(x,\omega)+D_{F}^{*}(x,\omega)+\cdots\right)}{\sum\limits_{nshot}\left(D_{F}^ {*}(x,\omega)+D_{F}^{*}(x,\omega)+\cdots\right)}\end{array}. \tag{15}\]