# Combining Higher-Order Statistics and Array Techniques to Pick Low-Energy P-Seismic Arrivals
## 3 Methodology
### Multi-Channel Stage  
The second part of the algorithm involves a multi-channel stage (MCS) that includes the analysis of multiple channels of an array of receivers. This stage exploits the array disposition information to enhance the reliability of the previous onsets that passed the multi-threshold criterion. In order to formalize this process, we introduce the following notation. Suppose that a seismic event is registered by the receivers \(r_{ij}\), resulting in the traces \(u_{ij}\), where the indices \(i\) and \(j\) run, respectively, over the number of arrays and the number of receivers of array \(i\). Let us assume that for all or a portion of the traces, \(u_{ij}\) STS is assigned an onset of \(o_{ij}\). For each array \(R_{i}=\{r_{i1},r_{i2},r_{i3}...\}\), we denote the set of onsets as \(O_{i}=\{o_{ij_{1}},o_{ij_{2}},o_{ij_{3}}...\}\). To evaluate the consistency of the arrivals, we first perform a clustering analysis of the elements of \(O_{i}\). In general, a clustering procedure is an unsupervised learning technique [33] used to group data points together, based on some similarity criterion, into subsets or clusters that are homogeneous or well separated [34]. In this study, we employed hierarchical agglomerative clustering [35; 36; 37; 38; 39] to build groups of concordant onsets. Hierarchical agglomerative clustering begins by partitioning the dataset into single dad (\(\{o_{ij}\}\) in our case), each representing a cluster. A value of dissimilarity is then computed for each pair of initial clusters. Subsequently, the clustering procedure iteratively merges the current pair of mutually closest clusters into a new object, computing new dissimilarities until there is one final cluster left. The user can stop the fusion of the objects at a certain level of dissimilarity or when a desired number of clusters is reached. Here, we adopt the complete linkage distance [37] as a measure of dissimilarity between two fused objects (clusters) of onsets \(C_{k}=\{o_{ik_{1}},o_{ik_{2}}...\}\) and \(C_{l}=\{o_{il_{1}},o_{il_{2}}...\}\), defined as follows:  
\[d(C_{k},C_{l})=\max_{o_{ik}\in C_{k},o_{il}\in C_{l}}|o_{ik}-o_{il}|. \tag{5}\]  
When clusters \(C_{k}\) and \(C_{l}\) are merged into a new cluster \(C^{\prime}\), the value \(d(C_{k},C_{l})\) will correspond to its diameter, that is, the maximum dissimilarity (distance) between two entities within \(C^{\prime}\)[40; 41]. A critical point of hierarchical agglomerative clustering is to establish the optimal number of final clusters or, equivalently, the maximum allowed dissimilarity between fused objects [34]. Given the nature of our objects (seismic arrival times), we fixed the maximum intra-cluster distance as the largest physically possible time difference between P-wave arrival times at receivers of array \(i\), which we indicate as \(\Delta t^{max}_{i}\). This parameter is related to the maximum spatial distance \(D^{\text{max}}_{i}\) between two generic receivers \(r_{il},r_{il^{\prime}}\) of the same array and the P-wave propagation velocity \(v_{p}\) in the following way:  
\[\Delta t^{max}_{i}=\frac{D^{\text{max}}_{i}}{v_{p}}+\delta, \tag{6}\]  
where \(\delta\) is a tolerance term, which takes into account any possible uncertainty or fluctuation in the considered physical quantities. Thus, we interrupted the clustering procedure at the fusion step, after which a new cluster with a diameter larger than \(\Delta t^{max}_{i}\) was created. In this way, the clustering step returns a partition of \(O_{i}\) in a set of different clusters \(\{C_{k}\}\) with the largest intra-cluster distance less than or equal to \(\Delta t^{max}_{i}\). Figure 4 illustrates some possible scenarios arising from the clustering procedure in a schematic way. Once the clusters of onsets at array \(i\) are formed, a selection of suitable onsets is performed as follows. Let \(N_{i}\) be the total number of onsets in \(O_{i}\), which is, in general, less than or equal to 10; that is, the number of receivers for our arrays.