## 2 Methodology  
### Influence of Attenuation Media on Depth-Domain Source Wavelets  
In viscoelastic theory, the degree of media attenuation is usually measured by the quality factor \(Q\), which is defined as follows:  
\[Q=\frac{2\pi E}{\Delta E}, \tag{1}\]  
where \(E\) stands for the energy of a seismic wave at a reference point and \(\Delta E\) is the energy loss after propagating one wavelength (or a period). The smaller the \(Q\), the more substantial the degree of attenuation, and conversely, when \(Q\) tends to infinity, the media can be regarded as elastic. Let us start with the analytical solution of the homogeneous viscoacoustic media to investigate the influence of the attenuation media on the waveform of the source wavelet. Assume that a homogeneous viscoacoustic medium has a velocity value of \(v_{0}\) and a quality factor value of \(Q_{0}\). When a time-domain source wavelet is excited at point \(x_{0}\), its analytical solution at point \(x\) is the convolution of a source wavelet and a three-dimensional (3D) Green function in this medium. The analytical solution can be expressed as follows:  
\[s(t)=f(t)*\widetilde{g}_{0,3D}(x|x_{0},t), \tag{2}\]  
where \(s(t)\) is the analytical solution, \(f(t)\) represents the source wavelet, and \(\widetilde{g}_{0,3D}(x|x_{0},t)\) stands for the 3D Green function. The 3D Green function is written as follows [38]:  
\[\widetilde{g}_{0,3D}(x|x_{0},t)=\frac{1}{4\pi|x-x_{0}|}\delta\bigg{(}t-\frac{ |x-x_{0}|}{v_{0}}\bigg{)}, \tag{3}\]  
where \(\delta\) represents the Dirac delta function and \(|x-x_{0}|\) denotes the distance between the source and receiver points. The result of the Fourier transform of Equation (2) can be rewritten as follows:  
\[s(\omega)=F(\omega)\cdot g_{0,3D}(x|x_{0},\omega), \tag{4}\]  
where \(s(\omega)\), \(F(\omega)\), \(g_{0,3D}(x|x_{0},\omega)\) denote the analytical solution, source wavelet, and 3D Green function in the frequency domain, respectively. The definition of the frequency-domain 3D Green function is as follows:  
\[g_{0,3D}(x|x_{0},\omega)=\frac{e^{-i\widetilde{k}|x-x_{0}|}}{4\pi|x-x_{0}|}, \tag{5}\]  
where \(\widetilde{k}\) is the complex wavenumber. The expression of \(\widetilde{k}\) is as follows [39]:  
\[\widetilde{k}=\frac{1}{v}\Big{[}\omega-i\tan\Big{(}\frac{\pi\gamma}{2}\Big{)} |\omega|\Big{]}, \tag{6}\]  
where \(v\) represents complex velocity and \(\gamma\) denotes the attenuation factor. Their definitions are as follows:  
\[\left\{ \begin{array}{c}\gamma=\frac{1}{\pi}\tan^{-1}\Big{(}\frac{1}{Q_{0}} \Big{)}\\ v=v_{0}\Big{|}\frac{\omega}{\omega_{0}}\Big{|}^{\gamma}\end{array} \right., \tag{7}\]  
where \(\omega_{0}\) is the reference angular frequency (in radians per second) and its corresponding velocity is \(v_{0}\). Therefore, substituting \(v_{0}\) and \(Q_{0}\) into Equation (7) leads to the numerical complex velocity and attenuation factor. Then the numerical complex wavenumber and the 3D frequency-domain Green function can be obtained. Substituting Equation (5) and the spectrum of the source wavelet into Equation (4) yields the frequency-domain analytical solution. Moreover, the time-domain analytical solution of the wavelet can be obtained by the inverse Fourier transform of \(s(\omega)\). In the depth domain, the complex velocity should be a function of the wavenumber:  
\[\left\{ \begin{array}{c}\gamma=\frac{1}{\pi}\tan^{-1}\left(\frac{1}{Q_{0}}\right) \\ v=v_{0}\left|\frac{k}{k_{0}}\right|^{\gamma}\end{array} \right., \tag{8}\]  
where \(k\) is the wavenumber, \(k_{0}\) denotes the reference angular wavenumber (in radians per meter), and its corresponding velocity is \(v_{0}\).