## 2 Methodology
### Source-Time Function Estimation  
The source-time function \(S(t)\) absorbing origin time and event magnitude is usually estimated by the Gaussian and Ricker wavelet [8; 37], which is calculated as follows:  
\[S(t)=-\sqrt{\pi}\omega_{0}\exp(-\frac{\omega_{0}^{2}}{4}(t-\tau_{0})^{2}) \tag{6}\]  
where, \(\tau_{0}\) is the time value corresponding to wavelet central time; \(\omega_{0}\) denotes the controlled dominant frequency, which is related with specified sensor frequency response. While [PERSON] [48] proposed a kind of Gaussian family wavelet based on fractional-order derivatives, whose frequency domain form is expressed as follows:  
\[\Phi^{(u)}(\omega)=\left(\frac{u}{2}\right)^{-u/2}\frac{\omega^{u}}{\omega_{0} ^{u}}\exp(-\frac{\omega^{2}}{\omega_{0}^{2}}+\frac{u}{2})\times\exp(-i\omega \tau_{0}+i\pi\left(1+\frac{u}{2}\right)) \tag{7}\]  
where, \(u\) stands for the fractional-order of time differential. The fractional-order Gaussian function can change fractional order \(u\) and \(w_{0}\) at the same time, thus, it has more freedom than the pure Gaussian or Ricker wavelet. However, for quite dissimilar waveform data, there are still some limitations in wavefield modeling using the two parameters (\(u\), \(w_{0}\)) fractional-order Gaussian function as source-time function. In order to further improve fitting performance for complex waveforms, we put forward a new approach to obtain the source-time function based on observed waveforms themselves, namely the DD-SVD-FD wavelet estimation. The detailed steps are shown as follows: firstly, each recording waveform containing the first half period of P-phase arrival is cut by the semiautomatic windowing method proposed by [PERSON] et al. [6]: (1) A two passes, four-poles, 200 Hz lowpass Butterworth filter was conducted to the recorded signal. (2) We manually picked a rough time \(T_{1}\) before the P-phase arrival and a rough time \(T_{2}\) after the first peak amplitude; (3) The \(T_{2}\) is automatically extended to the next zero amplitude time \(T_{3}\); (4) The windowed waveform corresponds to the filtered time series in the time interval \([T_{1},T_{3}]\). The waveforms are aligned based on the cross correlation time difference between different windowed waveforms. Then, the SVD average of windowed and aligned signals is taken as a basis function, where SVD is a classical signal processing approach based on coherence between different data traces. It can enhance the signal-to-noise ratio and suppress incoherent background noise [49]. The aligned and windowed data traces can be represented by a data matrix \(S(\mathbf{x}_{r},t_{i})_{M\times N}\) (dimension is \(M\times N\)), with \(M\) traces and \(N\) time samples per trace (usually \(M<N\)). Assuming that the data matrix \(S(\mathbf{x}_{r},t_{i})_{M\times N}\) has a rank \(R<M\), the SVD of \(S(\mathbf{x}_{r},t_{i})_{M\times N}\) is a orthogonal expansion of the data space:  
\[S(\mathbf{x}_{r},t_{i})_{M\times N}=\sum_{i=1}^{R}\sigma_{i}u_{i}\tau_{i}{}^{T} \tag{8}\]  
where, \(\sigma_{i}\) is the \(i\)th nonzero singular value, and \(u_{i}\) and \(v_{i}\) are the corresponding left and right singular eigenvectors with dimensions \(M\) and \(N\), respectively. The product \(\sigma_{i}u_{i}\sigma_{i}{}^{T}\) is a \((M\times N)\) matrix with rank one. If the windowed dataset has a high degree of trace-to-trace coherence, then the most coherent or average trace can be constructed from the largest singular eigenvector with the largest eigenvalue. The SVD acts as a data-driven, lowpass filter by keeping the main character of the dataset and rejecting highly incoherent noises. Finally, time scaling factor and fractional-order are introduced to estimate the frequency-time domain characteristics and obtain the source-time function in the frequency-time domain as (Equation (9)). \[\left\{ \begin{array}{l}\overline{S}(w)=\int_{-\infty}^{+\infty}S(t)e^{-iwt} \mathrm{d}t\\ \overline{S}(at)=S(w/a)/a\end{array} \right.