Deep embedded clustering of coral reef bioacoustics^a) Free

Deep learning, a powerful, recent subfield of machine learning,¹ excels at learning representations of large amounts of data and often outperforms shallower machine learning methods.^2–8 Deep convolutional networks have been particularly effective for image classification due to their scalable capacity. In ocean bioacoustics, machine learning has become an effective signal detection and classification tool^2,9–25 with research encompassing both supervised methods,^9–13,18 in which reliable labeled data are present, and unsupervised methods^16,17,19,26 when labeled data are limited or unavailable. Deep convolutional learning of marine mammal signals in the time-frequency domain has shown promising results for supervised detection and classification.^15,20–25 For other ocean bioacoustic signals, particularly marine fishes, just a few studies have considered deep learning classification.^26–28 This study proposes a variant of deep convolutional learning dubbed deep clustering as a method for time-frequency representation learning and classification of unlabeled, automatically detected signals in a cacophonous coral reef environment, characterized by significant presence of spatially dense calling fish.

Coral reefs are amongst the most biodiverse ecosystems in the ocean but are under threat from global climate change, overfishing, and pollution.^29–32 Passive acoustic studies allow non-invasive study of reef ecology over continuous time scales of days to years, and offer a complement to measurements that traditionally have been collected through direct observation by SCUBA divers or snorkelers, taking point measurements in time.³³ Coral reef ambient biological sound, or soundscapes, are an emerging topic of interest in the coral reef scientific community. Reef sound has been linked to relative ecosystem health, abundance of both reef building coral and fleshy macroalgae, and fish density.^34,35 Healthy reef sounds have also been shown to enhance larval recruitment,³⁶ suggesting that reef soundscapes are not just a by-product of biological activity but an integral part of complex ecosystem function. The biological components of reef soundscapes comprise of inadvertent sounds from organism activities such as feeding, locomotion,³⁷ and photosynthesis,³⁸ and acoustic communication such as the fish calls discussed here. Acoustic classification of individual fish calls informs our understanding of spatial and temporal movement, species assemblage, and response to human activities.

Acoustic classification of tropical marine fishes, such as damselfish (family Pomacentridae), has been improved through passive acoustic field experiments^39–41 but lacks established terminology across studies.⁴² To reduce the labeling burden, a few studies have considered automatic classification of fish calls by utilizing machine learning tools.^18,26 Unsupervised methods such as Gaussian mixture models (GMM) showed improved detection of fish chorusing events compared to the conventional spectral energy detector.²⁶ Deep neural networks, including convolutional networks (CNN), recurrent networks (RNN),²⁷ and sparse autoencoders (SAE)⁴³ have achieved higher classification accuracy of grouper calls than a sparse feature classification method.²⁸ 

In this study, deep clustering is applied to automatically detected, unlabeled signals, combining the benefits of deep learning with the flexibility of clustering. Due to the complexity of the dataset, clustering is restricted to two distinct signals: low-frequency fish pulse trains (fish) and segments of whale song units (whale). Standard GMM-based clustering and other conventional clustering of handpicked features are conducted for comparison.

1. Spectral and time domain features were manually chosen, or handpicked, based on their observed relation to coral reef fish calling and on studies of fish calling spectral and temporal properties.^39–41 Clustering of events using handpicked features provides physical intuition about the event signals but is difficult due to the varying feature properties. Thus, the feature extraction and classification steps are combined into one algorithm using deep learning.

2. Fixed-length spectrograms were used in deep embedded clustering (DEC),⁴⁴ a deep-learning image compression algorithm that produces an accurate image reconstruction. The DEC latent features were directly clustered with GMM, then the DEC was trained further using a joint clustering loss to encourage deep cluster formation among the latent features.^45–47 

Simulated signals were used to compare the limitations of the two approaches before applying them to data recorded on a Hawaiian coral reef between February and March 2020.⁴⁸ 

In Sec. II, the clustering methods are overviewed and metrics are presented for measuring classification success. Section III outlines the simulated signals used for method comparison. Experimental data collection and automatic signal detection on a Hawaiian coral reef in February–March 2020 are discussed in Sec. IV along with discussion of handpicked feature extraction. A comparison of clustering results on the simulated datasets are presented in Sec. V A, with an analysis of the effect of the DEC latent dimension. Finally, Sec. V B presents DEC, GMM, and conventional clustering of two signal types for the experimentally detected signals. Section VI summarizes the approach and discusses challenges associated with the methods and dataset.

Clustering methods are unsupervised frameworks for categorizing data according to their similarities.⁴⁹ This section describes DEC for learning the latent feature vectors and forming clusters directly from high-dimensional data. For comparison, the GMM, K-means, and hierarchical agglomerative clustering algorithms, common methods that rely on approximations of the input feature vector properties, are also presented.

DEC is a modified convolutional autoencoder, a deep feature learning method,⁵⁰ that uses a joint loss function to transform a set of N input data, $X={X1,…,XN},Xn∈ℝd1×d2$⁠, into a set of latent feature vectors, $Z={z1,…,zN},zn∈ℝP$⁠. The output is a reconstruction of the input data, $X̂n∈ℝd1×d2$⁠, estimated using the latent feature vectors. The DEC structure consists of two stacked networks (Fig. 1): the encoder network, which maps input data into a lower-dimensional latent space, and the decoder network, which reconstructs an approximation of the input from the latent space. For a given latent feature dimension, P, DEC learns representative latent features using the mean squared reconstruction error of the inputs,

Inspired by recent applications to spectrograms of unlabeled seismic events,^46,47 the compressive architecture downsamples the input using convolutions of stride length 2. The rectified linear unit (ReLU) was used to transform the outputs at each layer. The neural network structure⁴⁶ was adapted by reducing the convolutional filter size and adding convolutional layers to downsample the input.

The DEC weights were pretrained for 1000 epochs using the MSE of fixed-length spectrogram images containing N_f frequencies and N_t time samples, $Xn∈ℝNf×Nt$⁠, and their reconstructions, $X̂n∈ℝNf×Nt$⁠. The Adam optimizer⁵¹ was used with learning rate of $10−3$⁠. The spectrogram length N_t corresponds to a signal duration of 0.5 s = N_tdT and was motivated by the signals of interest in the coral reef soundscape. Events longer than 0.5 s were clipped to length to ensure a fixed size.

Each latent feature vector $zn$ is assumed to belong to one of K clusters, C_k, with means $μk∈ℝP, k=1,…,K$⁠. The latent features $zn$ were clustered with GMM to initialize the deep clustering, with means $μk, k=1,…,K$⁠. These clustering results are reported as “GMM, latent features” for the simulated and experimental data.

Last, the DEC latent feature vectors were updated by incorporating the Kullback-Leibler (KL) divergence into the training loss. DEC was trained for an additional 20 epochs using the joint clustering/reconstruction loss function^44,45

where Eq. (4) is the empirically estimated Student's t-distribution of the latent feature vectors around each cluster mean and Eq. (5) further penalizes points that are distant from cluster centers.⁴⁵ The objective of additional training is to learn latent feature vectors that are both disjoint in feature space and representative of the inputs. The weights for each loss term in Eq. (2) are from Ref. 45. Details of the model are given in Fig. 1 and Table I. The DEC was written in keras⁵² using tensorflow.⁵³ 

GMM aim to partition a set of N feature vectors, $X={x1,…,xN}$ into K clusters, C_k, $k=1,…,K$⁠, assuming the kth cluster is Gaussian with prior probability π_k, mean $μk$⁠, and covariance $Σk$⁠.⁵⁴ C_k is a nonempty set. If the number of clusters, K, is assumed a priori, then the Gaussian parameters $θk=(πk,μk,Σk)$ can be iteratively estimated.

GMM employs the expectation-maximization (EM) algorithm to estimate $θk$ from $X$ via the complete data likelihood.⁵⁴ EM is an alternating algorithm that updates the weighted posterior probability (responsibility, r_nk) and the prior, mean, and covariance by⁵⁴ 

for $xn∈X, n=1,…,N$⁠. $θk={πk,μk,Σk}$ represents the estimated parameters of the kth cluster at the tth iteration. In GMM, the data likelihood p is Gaussian,

where $∝$ represents proportionality. EM is iterated until π_k, $μk$⁠, and $Σk$ converge for all k. The estimated covariance matrix can become ill-conditioned if there are fewer than K true clusters.

The cluster priors π_k, means $μk$⁠, and covariance $Σk,k=1,…,K$ were initialized for 100 starts using the K-means++ algorithm (see  Appendix A).⁵⁵ 

K-means⁵⁶ is an approximation to the EM algorithm that partitions a set of data-derived feature vectors $X={x1,…,xN}$ into K clusters, C_k, $k=1,…,K$⁠, where C_k is a nonempty set. The number of clusters K must be assumed a priori. The K-means algorithm is derived assuming the feature vectors belonging to the kth cluster are drawn from a Gaussian distribution with covariance $Σ=σ2I$ and prior probability $πk=1/K$⁠, where $I$ is the identity matrix,⁵⁷ 

for $n=1,…,N$⁠. $|Ck|$ denotes the cardinality or size of the kth cluster.

The cluster means $μj, j=1,…,K$ were initialized for 100 starts using the K-means++ algorithm ( Appendix A).⁵⁵ 

K-means is commonly used due to its computational speed and interpretability, but its validity is subject to the isotropic covariance assumption and the equality of the cluster priors. If the true number of classes differs from K, the clusters will be incorrectly estimated.

Agglomerative hierarchical clustering,^49,54 also called bottom-up clustering, partitions a set of N data-derived feature vectors, $X={x1,…,xN}$⁠, into K clusters C_k, $k=1,…,K$⁠, by grouping the most similar data at each step. Hierarchical clustering successively merges nearby clusters until the stop criterion is achieved. In this case, the stop criterion is met when K or fewer clusters remain, where K must be set by the practitioner.

Agglomerative clustering does not require the number of clusters to be assumed in advance. As demonstrated with classification of dolphin echolocation clicks,^16,17 bottom-up methods allow for detailed clusters to be sequentially merged as additional information or labels become available.

In agglomerative clustering, the kth cluster at the mth iteration is represented as $Ckm, km=1,…,N−m+1$⁠. The number of clusters is decreased to N – m by merging the two closest clusters that satisfy

$d(Ci,Ci′)$ is a metric measuring the distance between all points in sets i and $i′$⁠. This agglomerative process is repeated $N−K+1$ times, until there are at most K clusters C_k remaining, with $k=1,…,K$⁠.

The distance metric d chosen here is from Ward's method.⁵⁸ Ward's method measures the incremental sum-of-squares resulting from merging two sub-clusters.⁵⁹ For clusters C_j and C_k with cluster means $μj$ and $μk$ that were merged for form $Ck′$ with mean $μk′$⁠, the Ward's method distance is

$|Cj|$ and $|Ck|$ are the cardinality of cluster C_j and C_k. Ward's method has been empirically successful,⁶⁰ but has been shown to perform worse in cases with unequal-sized clusters.⁶¹ Single-linkage, centroid, and complete-linkage methods were also considered, but had low accuracy.

The performance of the handpicked feature clustering and DEC was measured by

where t_i is the true label and $t̂i$ is the categorical prediction. The indicator function $I(x)=1$ if x is true and $I(x)=0$ otherwise. Here, whale (label 0) was the true category. Fish were indicated by label 1. The label 2 was used in a scenario to indicate overlapping fish/whale signals. Precision (positive predictive value or PPV) measured the ratio of correctly predicted whale to total predicted whale. Recall (hit rate or true positive rate) measured the ratio of correctly classified whale to the true total. Higher metrics correspond to improved performance, with perfect performance given by accuracy, precision, and recall all equal to 1.

A set of simulated coral reef bioacoustic signals, motivated by observed signals from a coral reef soundscape near Hawaii, was simulated to compare conventional clustering of handpicked features to DEC.

One hundred simulation trials were conducted, with 10 000 signals simulated for each trial. The signals were simulated as timeseries sampled at 1 kHz and mimicked whale and fish signals (Fig. 2). A spectrogram image of each signal was generated using a 256–pt FFT with 90% overlap. The fixed-length spectrograms were created by clipping all signals to 0.5 s.

The DEC latent feature dimension was optimized over a range of values (Fig. 7). The optimal DEC performance was compared to clustering with automatically extracted handpicked features. Then, a third category containing an overlapping fish and whale signals was included, with 33% of the total samples belonging to each category.

Quadratic FM sweeps mimic parts of humpack whale song units (whale). The equation for the instantaneous frequency of a quadratic sweep is⁶² 

where t is the continuous time variable, f₀ is the initial sweep frequency, $Δf$ is the total bandwidth, and T is the duration of the signal. The signal is an FM upsweep when $β>0$⁠, an FM downsweep when $β<0$ and a tonal when β = 0. The phase $Φ(t)$ of the time domain signal is found by integrating the instantaneous frequency⁶² 

where t₀ is the signal start time. For simulation, t is discretized at 0.001 s (sampling frequency of 1 kHz).

Timeseries of impulses, or pulse trains, simulated fish pulse calling (fish). The pulses were a set of N_p superimposed Gaussian-modulated sinusoids⁶³ spaced $Δt$ apart

where τ is the half-power pulse width and f_c is the center frequency.

The signal parameters were varied randomly for each sample (Table II). Pulse width and center frequency were fixed to achieve a representative pulse structure. The number of pulses and spacing were drawn from experimentally estimated distributions. The duration, initial frequency, and total bandwidth of the FM sweep were drawn at uniform random from a range of realistically observed values. All signals were centered within the 0.5 s spectrogram and assigned a random delay of within ±0.1 s.

White noise was added to the simulated signals using a fixed signal-to-noise ratio (SNR),

where $σs2$ is the signal power and $σn2$ is the noise power. The SNR range of each signal was determined from the experimental spectrograms during manual labeling. The SNR in dB was estimated as the difference of the peak signal power to the median power of the background.

The signal power was estimated as the bandwidth-normalized mean power over the signal duration,⁶⁴ 

Three handpicked features were extracted: peak frequency (25), kurtosis (26), and number of timeseries peaks (27). Duration, median power, and cross-sensor coherence were excluded due to limitations of the fixed simulation parameters. Then, GMM, K-means, and hierarchical clustering were applied to the feature matrix to discover K = 2 clusters (fish or whale) or K = 3 (fish, whale, or both).

Three directional autonomous seafloor acoustic recorders (DASARs) were deployed adjacent to a coral reef westward of the island of Hawaii. The DASARs, labeled N, W, and S for the north-most, middle, and south-most sensors, measured pressure and horizontal particle velocity with x– and y–components oriented at orthogonal compass directions. The array was roughly oriented N-S with inter-sensor spacing about 15 m (Fig. 3).

The DASARs recorded continuously for 7 days with a sampling rate of 1 kHz. On February 25, 2020, the dominant soundscape contributors below 500 Hz were reef fish, humpback whales, and motor noise from transiting surface boats, with boat noise occurring predominantly during daylight and fish calling most pronounced during dusk. The data were processed in 5 min chunks to account for DASAR clock drift. For each chunk, the pressure spectrograms on the three sensors were cross-correlated and aligned to within a time bin.⁶⁵ The spectral density power spectrogram, denoted as matrix $X∈ℝNf×Nt$⁠, was generated using a 256–point FFT with 90% overlap and Hanning window and normalized to have units of $log$ μPa Hz⁻¹, with $dt=0.026$ s and $df=3.91$ Hz.

A directional event detector⁶⁵ was developed to utilize the DASARs' directional capability by combining two DASARs ( Appendix B). For this study, the detector parameters were $Δθ=90°$⁠, T = 120 Hz, $Msep=1$(⁠$Msepdt=0.0256$ s), and $Tmax=2$ s. Detected events were localized⁶⁶ to ensure that their signal had sufficient bandwidth for feature extraction. Detections for which the localization algorithm failed to converge were discarded. The remaining events were spatially filtered within a 100 m by 100 m box from DASAR S (Fig. 3). A total of 92 736 potentially localizable detections within 100 m were kept for further analysis, on average about 1 detection per second.

The handpicked features (Table III, Fig. 4) are time-frequency properties known to relate to fish call type including power, duration, and peak frequency^39–41 as well as timeseries estimates of impulsive noise.^18,67

If $tn,1$ and $tn,2$ represent the absolute start and end of the nth detected signal in seconds, the event duration is

The median power and peak frequency of the nth detection were extracted from the spectrogram, [Fig. 4(a)–4(c)], with

The spectral features extracted at DASAR S, closest to the majority of detected signals (Fig. 3), were used for clustering.

The pressure timeseries of the nth detected signal $yn∈ℝM$ was extracted between $tn,1$ and $tn,2$ (Fig. 4) and contained M time samples. For the experimental data, the vector sensor x– and y–velocity during the nth detected signal, $vn,x∈ℝM$ and $vn,y∈ℝM$⁠, were used to create a beamformed pressure timeseries, $yb∈ℝM$⁠, for improved detection SNR,

where $Z0=ρc$ is the impedance in water with density ρ (kg m³) and sound speed c (m s⁻¹).⁴⁸ $θ̂$ is the signal azimuth estimated during the directional detector processing.⁶⁵ For the simulated data, $vn,x=vn,y=0$⁠.

Three metrics were chosen for timeseries extraction: kurtosis, number of peaks, and cross-sensor coherence. The kurtosis is a ratio of moments and has recently been applied to the task of differentiating impulsive from non-impulsive sounds,^18,67

$σ2$ is the variance and the overbar represents the arithmetic mean.

The number of peaks and cross-sensor coherence were extracted from the Hilbert transform $H()$ of the beamformed pressure timeseries (Fig. 4)

Here, the number of peaks is the number of local maxima with at least 5 dB prominence relative to the standard deviation

where (a, b) is an interval in $Yb, maxminjYb[j]$ indicates the highest trough in (a, b), and σ is the standard deviation of Y_b.

Last, for the experimental detections, the normalized correlation coefficient of the timeseries envelope across DASARs was computed to measure the spatial coherence of the signal propagation,

where $YbN$ and $YbS$ are the Hilbert transform of the beamformed timeseries at the North and South DASARs. The measure of coherence assumes that a good quality biological signal will be received as similar timeseries on two spatially separated sensors.

Figure 5 shows the number of events detected along with extracted feature median, 10%, and 90% levels for every 15 min. The number of peaks in the timeseries has been shown to be an indicator of fish species and call context for Hawaiian reef fish.^40,41 In this study, fish were most common during nighttime [Fig. 5(g)], with pulse trains peaking during the evening chorus after nautical twilight (19:15 HST). The evening chorus corresponded to a visible increase in median power, event duration, and number of time peaks and a visible decrease in the 90th percentiles of peak frequency and kurtosis.

The handpicked features were examined for separability for K = 2 equal-sized clusters (fish, whale) and K = 3 clusters (fish, whale, both). Figure 6 shows the feature values colored by true signal type for one simulation trial. The whale and fish signals overlapped in their peak frequency and in the number of automatically extracted temporal peaks. The signals were most strongly separated in kurtosis, with whale signals having very low kurtosis.

For two clusters, all clustering methods assumed K = 2. When the cluster sizes were equal, all clustering methods achieved high accuracy (Table IV). For the unequal cluster scenario with 25% whale and 75% fish (2500/7500), GMM achieved much higher accuracy than either K-means or hierarchical agglomerative clustering on the handpicked features (Table IV). These simulation results are supported by a supplementary example of the clustering algorithm assumptions ( Appendix A), which demonstrates that GMM has fewer prior assumptions and is thus more flexible. The classification accuracy of all methods for three equal-sized clusters (Table IV), with whale, fish, and a cluster containing a temporally overlapping whale and fish, was lower than for two equal-sized clusters due to the signal overlap. GMM is capable of estimating overlapping clusters and achieved the best accuracy and recall on the handpicked features of the overlapping signal scenario.

The effect of the DEC latent feature dimension was examined across 100 random simulation trials (Fig. 7). The DEC model, motivated by a successful architecture for classifying seismic spectrograms,⁴⁶ was retrained and tested on all trials at varying latent feature vector dimensionality (P) for the K = 2 equal-sized clusters and K = 3 cluster scenarios (Fig. 7). When the latent dimension was low, e.g., P < 8, all features were likely to be clustered in a single cluster. Above P > 8, DEC accuracy was consistent. Values of P = 10 and P = 15 were chosen for the K = 2 and K = 3 cluster scenarios.

Figure 8 shows t-SNE visualizations from a representative trial of the latent features clustered with GMM [Figs. 8(b), 8(e), and 8(h)], and the DEC latent features and cluster predictions after 20 additional training epochs using (2) [Figs. 8(c), 8(f), and 8(i)]. Additional training with (2) increased the accuracy of DEC clustering prediction from 96% to 99.2% in the case of 2 equal sized clusters and from 70.8% to 78.2% in the case of 3 clusters. In the scenario containing 2 unequal-sized clusters, the additional DEC training reduced the accuracy from 97.6% to 78.3%. In that case, high recall and low precision (Table IV) indicate that DEC overpredicted the smaller cluster. GMM clustering of the latent features learned from spectrograms achieved comparable accuracy to GMM clustering of handpicked features in the non-overlapping scenarios with K = 2 (Table IV). For the K = 3 scenario, the highest classification confusion for the spectrogram-based deep learning was between fish and the overlapping signal cluster. Figure 9 shows that fish signals with large bandwidth and duration may dominate the spectral signature, which may lead to increased misclassification.

The clustering methods were applied to a subset of 4000 unlabeled detections randomly selected from the Hawaii 2020 experiment. Each detection contained one or more directional signals. All clustering methods assumed K = 2. Then, for post-clustering analysis, labels of whale/no whale (only fish) were manually assigned to all 4000 samples, based on the signal within the detection window. About two-thirds were labeled as no whale and contained only fish.

Figure 10 shows the clustering results using t-SNE representations of the learned latent feature vectors [Figs. 10(a)–10(c)] and the DEC latent features after additional training [Fig. 10(d)]. Clustering of the P = 6–dimensional handpicked feature vectors did not correspond well to the manual labels, with K-means achieving the highest accuracy of 67%. GMM clustering of the latent features learned from the signal spectrograms achieved an accuracy of 78% compared to the manual labels, which was the highest accuracy of all methods examined. GMM of latent features had high precision and low recall (Table V).

In the experimental signals, additional training of the DEC latent features reduced the classification accuracy [Fig. 10(d)]. Similar to the simulated unequal-sized cluster scenario, DEC obtained higher recall with lower precision (Table V).

Spectrograms correctly classified by DEC are shown in Fig. 11, with the reconstructed spectrograms shown before and after additional DEC training using (2). These demonstrate that whale signals were primarily identified by their narrow bandwidth and temporal extent, whereas fish signals were identified as broadband. The timeseries in Fig. 11 demonstrate the magnitude variation between events that was not attributed to signal type, motivating the normalization of the spectrograms.

A deep clustering approach was presented for interpreting unlabeled coral reef bioacoustic detections. This approach leverages deep learning and clustering, motivated by recent improvements in classification accuracy of fish calling using unsupervised detection and deep neural network classifiers.^18,26,27,43

Clustering of simulated fish (Gaussian pulses) and whale (FM sweeps) demonstrated that unique signal categories could be clustered with GMM of deep latent features learned from spectrograms, DEC clustering,^45,46 or with clustering of handpicked features motivated by studies of coral reef fish.^39–41 GMM applied to either deep latent features or handpicked features achieved the highest accuracy when there were two clusters of unequal size, and GMM of handpicked features had the highest accuracy when the clusters contained overlapping signals. Deep clustering learned features directly from spectrograms and was successful in separating fish and whale clusters. DEC trained with a joint clustering loss function increased the classification accuracy in simulations compared to clustering the deep latent features with GMM, except when the cluster sizes were significantly different. For unequal-sized simulated clusters, applying GMM to the learned latent features demonstrated higher classification accuracy.

Broadband bioacoustic events detected on a Hawaiian coral reef in February–March 2020 using a directional detector were analyzed with deep clustering and clustering of handpicked features. A labeled subset of these detections with whale/no whale indicated that about two-thirds of the detections contained primarily fish and one-third contained a whale signal. GMM with K = 2 applied to latent features learned from the signal spectrograms achieved the highest accuracy and tended to overpredict the larger cluster containing fish. After further training DEC with a joint clustering loss, the DEC recall increased but DEC classification accuracy was reduced, indicating that whale was overpredicted. DEC reconstructions of the input spectrograms demonstrate that the learned features are representative of spectral features identified by manual labelers. Clustering of handpicked features with K = 2 on the experimental, manually labeled data achieved lower accuracy, precision, and recall.

These results demonstrate that deep clustering is a promising method for classifying unlabeled bioacoustic signals with distinct spectral signatures. In simulated and experimental studies, both GMM applied to deep latent features and DEC trained with a joint clustering loss achieved competitive classification accuracy for a range of latent feature dimensions, without requiring the selection of handpicked features. Our results indicate that the accuracy of handpicked feature clustering depends strongly on the feature properties and the choice of clustering method. Finally, the choice of clustering algorithm is an important consideration for applications that are subject to signal clusters of unequal size.

Thank you to Greeneridge Sciences for providing the DASAR sensors and to Alex Conrad for assisting with DASAR post-processing. Thank you to Richard Walsh for assistance with hardware deployment and recovery. This work was supported by the Office of Naval research under Award No. N00014-18-1-2065.

For a continuous data variable, x, that is Gaussian with mean $μk$ and covariance $Σk$⁠, the posterior probability of the cluster is given by⁵⁷ 

where C_k and C_j represent the kth and jth classes and $|Σk|$ is the determinant of $Σk$⁠. The boundary between the two classes occurs when there is an equal probability of x belonging to either class,

The general solution for x in Eq. (A3) is a multi-dimensional parabola, which simplifies to a linear boundary if the distributions have a shared covariance such that $Σk=Σj∀j,k$⁠. The solution in Eq. (A3) is extensible to K > 2 classes by considering the joint distributions of all classes (see Bishop⁵⁷ Chap. 4.2.1 for details).

Three 2D Gaussian distributions, each with N = 10 000 points, were used to simulate overlapping clusters (Fig. 12). The true cluster means were $μ1=(0,2), μ2=(10,$−$8)$⁠, and $μ3=(21,3)$⁠. The covariance of the first dataset was

with $σx2=σy2=3$⁠. For the first dataset [Fig. 12(a)], there were no off diagonal covariance terms.

The second dataset [Fig. 12(b)] was generated by rotating the data counterclockwise at θ, with

with $σx2=6, σy2$ = 3. The three clusters were rotated by $θ=120°,25°,$ and $0°$⁠. The non-isotropic covariance is a result of correlation between variables and indicates feature correlations. A non-isotropic covariance matrix of a feature vector indicates correlation between features.

GMM assumes that the clusters are Gaussian distributed around a cluster mean. The K-means algorithm further assumes the Gaussians have the same diagonal covariance and uniform prior probability. Thus, GMM and K-means perform best for the first dataset [Figs. 12(a), 12(c), 12(e)] with identical, diagonal covariance for all clusters. By incorporating intercluster distance through the Ward metric, hierarchical agglomerative identifies 3 similar clusters.

When the cluster covariances are not of the form $σ2I$⁠, K-means and agglomerative clustering provide solutions that are suboptimal to the true class boundaries. GMM achieves better accuracy [Figs. 12(b), 12(d), and 12(f)] but does not converge to the maximum likelihood solution due to estimation errors in the covariances. All three methods provide insight into the cluster memberships despite invalid class assumptions, but GMM is the most versatile clustering method.

K-means++ is an algorithm proposed to efficiently initialize parameters for Gaussian clustering methods, including GMM and K-means. K-means++ improves the runtime of the clustering algorithms and the quality of the final solution.⁵⁵ For a set of data samples $X={x1,…,xN}$⁠, the K-means++ algorithm follows:⁵⁵ 

1. Initialize the priors $πk=1/K, k=1,…,K$⁠.

2. Initialize the covariance matrices as $Σ=σ2I$⁠, with $σ2=$ var $[X]$⁠.

3. Select the first cluster mean as a random data sample, $μ1=xl, l∈U(1,N)$⁠.

4. Assign all data samples to the cluster with the nearest mean. We denote that the mth sample belongs to cluster C_i as $xm∈Ci$⁠, satisfying $d(xm,μi)≤d(xm,μp), i,p=1,…,j, j≤K$⁠. Note that when j = 1 (first step), $xm∈C1∀m$⁠, i.e., all data samples belong to one cluster.

5. Select the j + 1-th cluster mean at random from the remaining data samples with probability

   $p(xn=μj+1)=d2(xn,μj+1)∑xm∈Cid2(xm,μi), n=1,…,N−j,n≠l.$

   (A6)
   d is the Mahalanobis distance,

   $d(xi,xj)=(xi−xj)TΣ−1(xi−xj).$

   (A7)

   Equation (A6) ensures that each mean is selected with a probability proportional to its distance from all existing means.

6. Repeat steps 4 and 5 until K cluster means are chosen.

The initialized parameters define $θ0=(πk,μk,Σk), k=1,…,K$ for the first E-step of the GMM algorithm in (6). For K-means, the prior and covariance are assumed fixed, and the means are updated according to Eq. (10).

For data with more than two dimensions, $xn∈ℝP$ for P > 2, clusters may be visualized by applying dimensionality reduction. In this study, 2D t-stochastic neighbor embedding (t-SNE)⁶⁸ was used to visualize P–dimensional features.

The similarity of one point, $xi∈ℝP$⁠, to another point, $xj∈ℝP$⁠, is found from the conditional probability that the points are neighbors within a Gaussian density with mean $xi$⁠,⁶⁸ 

for N points i, j = $1,…,N$⁠. The neighborhood of $xi$⁠, as determined by σ_i, is defined implicitly in terms of the perplexity (Fig. 13),^68,69

where $Pi=Σjpj|i$⁠, and H is the Shannon entropy. The optimal value of σ_i in Eq. (A8) for each point is solved with a binary search for a given value of perplexity.⁶⁹ 

Then, a set of projected data is randomly initialized with zero-mean Gaussians of low variance,⁶⁸ $yi∈N(0,10−4I), yi∈ℝ2$⁠. The $yi$ are iteratively updated until the Kullback-Leibler (KL) divergence between the conditional probability in Eq. (A8) and the Student's t-distribution of the projected data is minimized,

In contrast to classic SNE which uses only Gaussians, t-SNE's use of the Student t-distribution further penalizes outliers in the projected data,⁶⁸ resulting in a more compact representation. As shown in Fig. 13, the value of perplexity should be varied according to user preference to obtain the desired visualization.

A directional event detector⁶⁵ was developed to utilize the DASARs' directional capability by combining two DASARs. The complex spectrograms of the x– and y–particle velocity, matrices $Vx∈ℂNf×Nt$ and $Vy∈ℂNf×Nt$⁠, were generated identically to the spectrogram X with units m s⁻¹.

The active intensity, a measure of in-plane energy, was used to determine the noise directionality

where the asterisk is the complex conjugate and $ℜ$ the real component. a tan² is a piecewise function that computes the elementwise angle between the elements of two matrices, with domain $(0°,359°)$ defined counterclockwise from the y axis (⁠$0°$ = North). A is therefore the time-frequency representation of compass directionality. In the following, matrix $AN$ is called the azigram⁷⁰ for DASAR N, likewise for $AW$ and $AS$⁠.

The event detector makes the following signal assumptions:

1. An event arrives from a constant azimuthal sector for each DASAR.

2. Target events are broadband below 500 Hz. The minimum required bandwidth was set with an empirical threshold (see Appendix in Thode et al.⁶⁵).

The detection algorithm is demonstrated in Fig. 14. First, the azigrams for the north- and southmost DASARs, $AN$ and $AS$⁠, were used to create binary maps $BN$ and $BS$ of time-frequency points within a fixed azimuthal sector [Fig. 14(a)],

and likewise for $BS$ [Fig. 14(b)]. $I$ is the elementwise indicator function, with $I(true)=1$⁠. Binary maps were generated for all combinations of azimuthal sectors $θN,$$θS∈([0,Δθ]T,[Δθ/2,3Δθ/2]T,…,[(360°−Δθ),360°])$⁠. Here, $Δθ=90°$⁠. Next, overlapping signals on both DASARs were discovered by creating a combined map [Fig. 14(c)], which was summed across frequency to determine the detection timeseries,

d measures the bandwidth of an event. Event start and end times were determined for $dj>T, j=1,…,Nt$ for threshold T. Events separated by less than $Msepdt$ were merged and events longer than $Tmax$ were removed.