# Combining Higher-Order Statistics and Array Techniques to Pick Low-Energy P-Seismic Arrivals
## 3 Methodology
### Single-Trace Stage  
STS aims to determine P-wave arrival times in the waveforms. It is based on the HOS approach introduced by [5]. The basic assumption is that the sought signal and the background noise are ruled by different physical processes, resulting in distinct statistics. The method utilizes HOS gradient-based characteristic functions. These functions possess the desirable property of yielding null values when computed on time windows containing Gaussian-distributed random processes, such as those consisting solely of noise, while exhibiting non-zero values when a signal is encountered. The onset (hereinafter denoted as \(o\)) of the signal is detected among the samples where the characteristic function significantly deviates from zero. For each seismic trace \(u(t),t\in\mathbb{N}\), we can define the time series \(\gamma_{n}(t)\), representing the n-th order relevant statistic. Indicating with \(E[\cdot]\) the expectation operator, for a generic discrete random process \(\{Y(t),t\in\mathbb{N}\}\) with expected values \(\mu_{t}=E[Y(t)]\), \(\gamma_{n}(t)\) is given by the following:  
\[\gamma_{n}(t)=\hat{\gamma}_{n}(t)-\hat{\gamma}_{n}^{G}, \tag{1}\]  
where \(\hat{\gamma}_{n}(t)\) is the estimation of \(\frac{E[(Y(t)-\mu_{t})^{n}]}{(E[(Y(t)-\mu_{t})^{2}])^{n/2}}\) in a time window of length \(w\) ending at time \(t\), and \(\hat{\gamma}_{n}^{G}\) is the value of \(\hat{\gamma}_{n}\) for a Gaussian variable [31]. Thus, the estimation of \(\gamma_{n}(t)\) for Gaussian variables is zero. Explicitly, for a seismic trace \(u(t)\) we have the following:  
\[\gamma_{n}(t)=\frac{1}{w-1}\frac{\sum_{i=t-w}^{t}(u(i)-\hat{m}_{w})^{n}}{(\sum _{i=t-w}^{t}(u(i)-\hat{m}_{w})^{2})^{n/2}}-\hat{\gamma}_{n}^{G}, \tag{2}\]  
where \(\hat{m}_{w}=\frac{1}{w}\sum_{i=t-w}^{t}u(i)\). We seek the P-arrival in a part of the trace \(u\) enclosed between the event origin time (OT) and OT + \(\bar{t}\). To increase the signal-to-noise ratio (SNR), we use a band-pass filter beforehand on the signals in the frequency range defined by \(f_{1},f_{2}\). Each time series is then scanned by a sliding time window of length equal to \(w\). The values of \(w,f_{1},f_{2},\bar{t}\) depend on the dataset and the geometry of the receivers and will be fixed later (Section 4). An HOS value is computed in each time window \(w\) sliding of one sample at each step, as described in Equation (2), resulting in a new time series: \(h(t)=\gamma_{n}(t)\). A moving average window of three points is used to smooth the \(h(t)\) and mitigate noise fluctuations. In the following, based on previous works [31; 32], we test two _HOS_ values: \(\gamma_{4}\) and \(\gamma_{6}\), which  
Figure 2: Brief outline of the proposed methodology, schematically illustrating the two-stage structure of the HOSA method. For a detailed description, refer to Algorithm 1 or Section 3.  
are, respectively, the 4-th and 6-th order statistics. As noted in [7], the time derivative of \(h(t)\) exhibits a considerably more impulsive character, making it more appropriate as a characteristic function. Consequently, derivatives are computed, generating the time series \(h^{\prime}(t)=\frac{dh(t)}{dt}\). Some typical waveforms, along with their \(h(t)\) and \(h^{\prime}(t)\), are depicted in Figure 3.  
How do we identify P-onsets from the \(h^{\prime}(t)\) functions? Several studies have simply used the maximum of \(h^{\prime}(t)\). However, as explained in [32], this criterion can lead to inaccurate detection in the presence of spurious spikes before the actual P-onset or in the case of very unstable \(h^{\prime}(t)\). This can be seen in Figure 3, where \(h(t)\) and \(h^{\prime}(t)\) have been plotted for three representative cases of seismograms in which the P-onset is clearly impulsive, emergent, and not clearly identifiable, respectively.