How to analyse overlapping sounds in the marine environment using supervised multi-label classification

Passive acoustic monitoring (PAM) is a non-invasive tool that can be used for long-term monitoring of the underwater environment^1,2,3. In addition to naturally occurring sounds in the ocean, human activity contributes to the complex underwater soundscape, e.g. through vessel noise, active sonar and pile driving⁴. PAM is frequently employed to study cetaceans, whose vocalisations provide the opportunity to identify the presence of species, or species groups, in an area, to analyse specific behaviours^5,6,7, to estimate abundance^8,9,10 and assess the impact of human activities^11,12,13.

The majority of cetacean species are highly vocal^14,15, utilising sound for social purposes and all odontocetes (toothed whales) also use it to sense their environment through echolocation^16,17. Cetacean vocalisations are usually broadly categorised as: burst pulses and whistles, which are used for social interactions, and echolocation clicks used for sensing their environment^17,18.

The increasingly common use of PAM has led to a growing need for automated detection systems capable of distinguishing sounds and advances in machine learning have greatly improved this capability^19,20. One common strategy is to divide a time series into a set of segments of fixed duration and compute a single label for each segment. Some methods are simple detectors, implementing binary classification^21,22,23 whereas others seek to perform single-label multi-class classification^24,25,26. In an acoustic environment where multiple audio events co-occur within a single segment, multi-label classification (MLC) systems, in which multiple labels can be assigned to one segment, are a natural tool.

In terrestrial acoustics, MLC systems have been frequently used for analysing overlapping sounds from different sources^27,28,29. However, there is a lack of corresponding research in marine bioacoustics, where the interactions of anthropogenic sounds and the vocalisations of aquatic fauna are widely studied^4,30,31. This is partly due to challenges such as the limited availability of labelled datasets and spectral overlap between anthropogenic and biophonic sound sources, as well as the greater complexity of underwater background geophonic noise. When designing architectures for MLC, there is a spectrum of approaches that aim to balance the use of shared network elements with overall network complexity. At one end lies binary relevance (BR), where independent binary classifiers are constructed for each class^32,33,34,35. Whilst at the other end are the Conventional Multi-Label Classifiers (CMLC), that employ a single shared network^27,28,36,37.

CMLC systems are optimised to learn general features across all classes²⁷. However, this generalisation may compromise performance compared to the tailored features learned in class-specific BR networks³⁷. BR systems, while capable of creating bespoke features for individual tasks, being trained independently cannot exploit inter-class relationships such as feature co-occurrences³². Consequently, in the field of MLC, intermediate architectures that lie between the extremes of BR and CMLC are often implemented using multi-task learning (MTL) strategies, which aim to balance feature sharing with task-specific specialisation^33,38,39,40. Two primary forms of MTL are commonly used⁴¹. The first is hard parameter sharing (HPS), where some, or all, of the feature extraction process is performed by a single network shared across all tasks⁴². The second is soft parameter sharing (SPS), which maintains distinct layers for the feature extraction process but facilitates information exchange between them⁴¹. There is a wide variety of methods for realising such information exchange^43,44. In this study, we use neural discriminative dimensionality reduction (NDDR) layers⁴⁵ as one such mechanism to integrate them within an SPS MTL framework, which we refer to as NDDR MTL system.

Class-specific binary detectors (BR) often achieve strong performance^32,34,35,37, although their computational requirements do not scale well as the number of classes increases, since a new network is required for each new class. Moreover, because the classifiers operate independently, BR cannot exploit correlations or naturally co-occurring structure between classes, which one might expect to degrade performance. Conversely, BR models are more flexible, since classes can be added or removed from a problem without necessarily retraining.

Accordingly, we compare the performance of four systems, BR, CMLC, HPS MTL, and NDDR MTL, which span the range of MLC approaches, for underwater sound classification by using data from the west of Scotland.

A comparison of the performance of MLC systems

The performance of any classifier should be made with reference to some baseline non-informative system⁴⁶. For MLC with an imbalance in the data, then one approach is to consider a system that predicts zero prevalence for all classes across audio segments⁴⁶. The resulting baseline accuracy using our validation dataset is 58.2% and baseline exact match is 31.4%.

Table 1 shows the accuracy and exact match scores for 4 systems considered. The performance metrics are computed for each of the folds with the mean and standard deviation across the 4 folds presented in the table along with the baseline scores. It is evident that all systems substantially outperform the baseline. Overall, the 4 systems offer similar levels of performance, with small differences in mean values (within 0.03). For the single network approaches (i.e. excluding BR), introducing partial layer sharing (as in HPS MTL and NDDR MTL) provides small improvements in average performance, consistent with previous MTL studies^40,41,47,48. The BR system achieves the highest mean performance across both metrics; however, the differences relative to the other systems are small, which shows all the systems have comparable performances under the present experimental setup, despite the larger overall model size of BR.

System performances under varying event concurrency

In multi-label sound classification tasks, different combinations of concurrent events can cause varying levels of difficulty^40,49. There are at least two factors affecting this performance, one is the number of different event types occurring concurrently within a sample (we refer to this as the concurrency level) and the other is the number of instances that a particular combination occurs in the dataset. One might hypothesise that performance on samples containing many different event types is poorer than on samples with fewer event types and, further, that performance on rarely occurring combinations is poorer than on combinations that occur frequently. The following results aim to explore these hypotheses.

To study the effect of event concurrency, we computed the macro averaged exact match values for each concurrency level, the results are shown in Fig. 1a for each of the networks. The data were grouped based on their concurrency level and the exact match scores were averaged across these groups. Exact match values of the systems for each unique combination are given in Fig. S1.

There is no consistent trend in the performance with increasing concurrency level (Fig. 1a). Whilst the best performance is obtained for the concurrency level of zero (no event present, i.e. ambient noise), there not a monotonic decrease in performance from that point: the next best performing case (for 3 of the 4 classifiers) is for concurrency level 4. This does not support the hypothesis that higher concurrency levels lead to poorer performance.

To shed light on the reasons for this, Fig. 1b shows the number of events in the dataset for each combination and the concurrency levels are indicated. Of the possible 32 combinations, only 18 are represented in this dataset. The correlation between presence of the event types leads to an uneven distribution of data across concurrency levels. For instance, because three of the classes are associated with the presences of dolphins, they frequently co-occur in some combinations, this is one reason why the higher concurrency levels (3 and 4) are more common in the dataset. This greater representation coincides with better system performances, especially evident at a concurrency level 4. Suggesting that the prevalence of combinations in the dataset is a more important factor affecting performance than the concurrency level.

To further evidence of the impact of the amount of data on system performance is provided by Fig. 2, which shows the exact match scores for each combination, averaged across the networks, against how many times that combination is represented in the validation dataset. There is a clear positive correlation between classification performance and the number of audio segments for specific combinations present in the data. Highlighting the critical role of the availability of sufficient audio segments for each combination in enabling accurate MLC.

One notable outlier in Fig. 2 is the label vector [C,0,0,V,S], which corresponds to the co-occurrence of delphinid clicks, vessel noise and sonar signals. Within our dataset, delphinid clicks naturally co-occur with delphinid vocalisations more frequently than with sonar signals. Furthermore, sonar signals, particularly echosounders, which typically occupy a similar frequency band to delphinid tonal, introduce similar spectral and temporal patterns as delphinid whistles. As a result, systems often (falsely) associate these sonar signals with delphinid tonals, leading to false insertions.

Confusion patterns across audio events

Identifying confusion patterns between audio events provides insights into the underlying reasons for system difficulties^26,50. In MLC, confusion must be examined at the combination level, as multiple events can co-occur within the same segment. To analyse these patterns, we computed a normalised confusion matrix showing the conditional probability of predicting a given combination for each ground truth combination (Fig. 3).

The confusion matrix was computed by first identifying all segments belonging to each ground truth combination across all systems and folds, pooling their predictions and counting how frequently each possible combination was predicted. These counts were then divided by the number of occurrences of the ground truth combination. To reduce visual clutter, only results for the commonly occurring ground truth combinations, specifically those which occurred more than 100 times in the training dataset, are shown. Further, we removed rarely occurring predictions, i.e. columns in which the prediction probabilities never exceeded 0.1 for any ground truth combination, so that the results shown contain only the interpretable and non-negligible prediction patterns.

The resulting confusion matrix has a strong leading diagonal, those values corresponding to the probability of an exact match for the given combination. This matrix provides greater insight regarding the outlier [C,0,0,V,S], highlighted in Fig. 2, 51% of this combination were predicted with combinations that include delphinid tonal (T), with the majority (36%) of those being predicted as [C,T,0,V,S], confirming the high rate of false positives of delphinid tonals. Further, the systems often get confused when sonar does not co-occur with tonals, for the combinations [0,T,0,V,S] and [C,T,0,V,S], the majority of the segments were correctly identified.

Delphinid vocalisations, delphinid clicks (C), tonal and burst pulses (B) frequently trigger insertions of one another. For example, 28% of [0,T,B,V,0] segments were predicted as [C,T,B,V,0], 14% of [C,0,B,V,0] were shifted to [C,T,B,V,0], and for [0,T,0,V,S] the systems inserted delphinid clicks in 40% of the cases and burst pulses in 11% of the cases. These interactions are also evident in false negatives, where one delphinid component suppresses another. For example, 33% of [C,T,0,0,0] samples were missing clicks and 12% were missing tonals. Similarly, 25% of [C,0,B,V,0] samples failed to detect burst pulses. This is mainly because delphinid vocalisations often naturally co-occur in varying combinations, which can confuse the systems when identifying event patterns.

While the insertion of tonals in the presence of sonar is a noticeable issue, the opposite pattern does not occur. When tonals are present, sonar is rarely inserted (there is only 4% insertion of sonar in total when delphinid tonal is present in the segments). Similarly, when sonar and tonals co-occur, the systems miss both the tonal and sonar together in 5% of the cases for both [0,T,0,V,S] and [C,T,0,V,S]. Therefore, the systems’ ability to identify sonar is not affected by the presence of tonals, in contrast to the presence of sonar causing the systems to incorrectly insert tonals. A prominent insertion of sonar occurs for ambient noise, where 17% of ambient noise was predicted as sonar. The main reason is the presence of sonar segments that propagate from long distances (especially sonar signals occurring without vessel noise), which exhibit very faint spectral patterns, causing systems identifying some ambient noise segments as sonar.

Vessel noise does not show a strong confusion pattern and is mostly identified correctly. The only prominent miss rate is 7% for [C,0,B,V,0]. In other cases, vessel noise is not missed on its own but only as part of broader misclassifications.

Segment-based PR curve analysis (per-class performance analysis)

To analyse the per-class performances of the systems and variabilities in the performances across the folds, we computed precision-recall (PR) curves with area under the curve (AUC) values (Fig. 4). All the systems achieve AUC values above 84.0% across the audio events, exceeding 90.0% when delphinid tonals are excluded and reaching above 98.0% for vessel noise. The greater variability in performance is observed for delphinid tonal class. This was caused by the uneven distribution of echosounder signals, particularly validation datasets for fold 1 and fold 2 have high representation of this signal type, causing a high false positive (FP) rate. A sharp drop in precision at low recall values, driven by these FP, yields the low AUC values.

Across all systems, vessel noise is the most accurately classified event with AUCs exceeding 98.0%, with the lowest variability. This low variability highlights this class’s predictability, indicating its less complex acoustic characteristics regardless of the system type. The factors contributing to this is that vessel noise is particularly, as a persistent tonal sound, more consistent across the audio segments, its being well-represented along with delphinid clicks (at least one of these two audio events are present in ~45.0% of the entire dataset), and the MS-PCEN input representation used in this study providing higher resolution at lower frequencies. This shows all systems abilities in identifying discriminative and well-represented events.

System model size effect

Experiments were conducted to understand the effect of the models’ size by systematically adjusting the number of trainable parameters across a range, from approximately 5 thousand to 12 million. This was realised by varying the number of convolutional kernels and recurrent units, while keeping all other architectural components and training procedures consistent with those described in Systems and Model Training Settings.

It is important to note that there are two common strategies to compare the performance of BR with other MLC systems in such experiments: (1) normalising BR’s total trainable parameters³⁷ or (2) keeping each BR network the same as that in the other systems^35,51. We followed the second approach because each BR network operates independently, meaning it does not share representations across labels like the other systems. Reducing BR’s total parameter count would not align with its independent-task structure, as each classifier learns separately without benefiting from shared features. Hence, we compare the systems as they naturally are, considering both total parameter counts and network structures.

The average accuracy values with min-max error bars across the four folds are shown in Fig. 5, where error bars illustrate the spread of performance across folds for each model size. As model size increases, all systems show improved performance and the performance differences between systems diminish. A key observation is that after 410 K trainable parameters, all systems reach a performance plateau. The differences between the 5 K and 410 K configurations are as follows: BR shows a 5.0% change in accuracy; NDDR MTL changes by 8.4%; HPS MTL by 10.7%; and CMLC by 13.8%. These results indicate that BR is particularly robust to variations in model size.

When comparing systems, BR and NDDR MTL generally achieve higher accuracy than CMLC and HPS MTL across nearly all model sizes. Notably, performance differences between systems become more apparent under reduced model capacity, where BR consistently maintains higher accuracy compared to the other approaches. The performance variability, as indicated by min-max error bars, shows the differences in the systems robustness. In general, variability increases as model size decreases. Overall, CMLC system has the highest variability among the systems. For example, at the smallest model size (5 K), CMLC experiences a sharp drop in accuracy (6.8%) compared to the next configuration (10 K), with the widest min-max spread of 20.0%. This shows the CMLC system’s sensitivity to the data distribution across folds and may struggle in domain adaptation or cross-site classification tasks.

This study provides an effective methodology for analysing marine polyphonic soundscapes that contain both biophonic and anthropogenic signals using MLC systems. By systematically comparing four supervised deep learning systems (CMLC, BR, HPS MTL, and NDDR MTL), our work highlights the strengths and limitations of each approach under realistic classification scenarios. We evaluated performance with respect to the number of concurrent event types present, the diversity of co-occurring event combinations, and model size.

The BR system has previously shown success in classifying a number of distinct sonar signals³⁷. Our results extend this finding to a more complex underwater task involving concurrent biophonic and anthropogenic audio events. Overall, the BR system achieved the best MLC performance, closely followed by NDDR MTL. The success of these two systems shows the value of class-specific network structures for accurately classifying co-occurring underwater sound events.

Although BR produced the best classification performance, its computational properties introduce practical limitations for large-scale PAM applications. Each BR classifier has training and inference times comparable to the full CMLC network because it uses the same underlying architecture (defined in ‘Methods’). BR’s total computational cost increases almost linearly with the number of classes; in our five-label task, this resulted in an approximately five-fold increase in inference time relative to the integrated systems (about four times longer than NDDR MTL). This limitation is expected to become increasingly critical as the number of target classes grows (e.g. applications involving larger acoustic taxonomies). For example, extending the current framework to 15–20 audio events would require a proportional increase in the number of independent BR classifiers, which substantially increases computational cost, whereas integrated MLC systems would maintain a single network with shared layer(s). The model size experiments provide additional insight into this trade-off. When the total number of parameters is reduced, BR exhibits the smallest decrease in performance, even at relatively small model sizes (e.g. the 5k configuration, corresponding to 25k total parameters for the system). This shows that the class-specific nature of BR classifiers allows them to retain strong performance under limited model capacity. However, this robustness at the classifier level does not completely offset the fundamental scalability limitation of BR at the system level, since achieving comparable per-class capacity requires maintaining one independent network per audio event, and normalising the total BR system size implies proportionally smaller binary classifiers. At the same time, we note that an important practical advantage of BR that requires training only an additional binary classifier when introducing a new audio event, whereas integrated systems require reconfiguring and retraining the entire network. On the other hand, integrated systems generate predictions for all classes with a single forward pass, which makes them computationally efficient when operating on long-duration PAM datasets or when scaling to a larger number of event types, particularly given that the performance gap between NDDR-MTL (10k model size) and BR (single binary classifier with 205k parameters) remains below 5% accuracy.

A key insight from our study is that the diversity of co-occurrence types and their representation in the training dataset have a greater impact on classification performance than the number of co-occurring classes alone. This contrasts with prior studies showing monotonic decreases in performance with increased concurrency levels^40,49, which used artificially mixed datasets. Our dataset only includes naturally co-occurring signals, and our findings emphasise the importance of co-occurrence frequency in training data. This is mainly because artificially mixed signals often lack coherent acoustic patterns, as events are combined without preserving their original environmental or behavioural contexts. Such combinations can result in unnatural overlaps where spectral components of distinct signals mask or distort each other, which makes them harder to disentangle. As a result, increasing the concurrency level in artificially mixed datasets consistently decreases the networks’ ability to distinguish between them.

We further observed consistent misclassification of delphinid tonal signals in the presence of high-energy echosounders occurring within a similar frequency band. This issue was not observed for sonar signals in lower frequency bands or for lower energy echosounders operating within the same frequency range. These findings support the hypothesis that automated systems may misclassify sonar signals as delphinid whistles, especially in regions where both signal types co-occur. Given that whistles play critical roles in identifying whistling animal and studying species behaviour^52,53, such misclassifications could hinder ecological analyses. Therefore, care must be taken in data collection to avoid regions dominated by such anthropogenic signals or to take additional measures to improve classifier robustness.

These observations can be further interpreted in terms of how the systems learn from the data. In the evaluated systems, each class is associated with its own output and loss, and therefore, the models are primarily trained to learn class-specific spectro-temporal features from the input representations. However, the results also show that model behaviour is influenced by the co-occurrence structure present in the dataset, as demonstrated in both the confusion patterns and the representation of different combinations in the data. These effects appear to vary across architectures. The BR system relies on independent classifiers and does not model shared representations across classes, whereas integrated systems (CMLC, HPS MTL, and NDDR MTL) incorporate shared or partially shared features. Lower performance of CMLC for the delphinid tonal class, as observed in the PR curve analysis, suggests that fully shared representations may not sufficiently capture the spectro-temporal characteristics required to distinguish acoustically similar events, potentially due to feature interference when events share similar spectro-temporal patterns. This effect appears to be more evident in CMLC, which also shows greater sensitivity to reduced model capacity in the model size experiments. In contrast, architectures that allow either class-specific feature learning (BR) or controlled feature sharing (e.g. NDDR MTL) appear better able to distinguish between spectro-temporally similar events.

Beyond the system-level findings, this study also highlights the broader value of multi-label approaches for marine bioacoustics. Many ecological questions require understanding not only the presence of individual sound types, but also how biophonic and anthropogenic signals co-occur within the same acoustic scene. Unlike single-label or binary detection frameworks, MLC systems can characterise these polyphonic conditions by providing a more realistic and diverse representation of marine soundscapes. Such capabilities are important for future ecological analyses, particularly in scenarios where the types or combinations of anthropogenic signals are unknown a priori. The framework established in this study can therefore support subsequent applications, including population- or behaviour-level studies, by providing a foundation for analysing complex acoustic mixtures. At Tolsta, the frequent co-occurrence of vessel noise with delphinid vocalisations shows that both sound sources are persistently present in the same acoustic scenes. This may suggest a regular overlap between anthropogenic activity and cetacean presence in the region, rather than short-term avoidance behaviour. Future work may extend this approach to new regions and/or incorporate additional event categories.

Our choice of Tolsta, a site reflective of Western Scotland²⁶, enabled this study to capture a range of co-occurring biophonic and anthropogenic sounds. However, system generalisation to other locations remains untested. Domain adaptation and cross-site validation have gained prominence as problems in underwater sound classification^50,54 and extending our study to include other datasets will be essential to fully assess generalisability. Nevertheless, our findings, particularly the observed variability in the CMLC system under different folds, already hint at potential differences in how well each system may generalise across acoustic environments. Additionally, although we evaluated five sound types (three delphinid vocalisations, vessel noise, and sonar), the framework is adaptable to additional biophonic, anthropogenic, and geophonic sources (e.g. baleen whale vocalisations, seismic airgun pulses, rain), contingent on annotated data availability.

We also note that we curated the full 101k-segment dataset to allow a controlled and fair comparison of the four systems, as class imbalance is known to adversely affect both training stability and the interpretation of results^55,56. The original dataset (Table 2) is heavily dominated by vessel noise, which would have masked system-level differences and made the comparison largely reflective of dataset bias rather than model behaviour. The curated dataset therefore had a specific methodological purpose, as outlined in Data Preparation and Experimental Setup: reducing extreme imbalance, improving the representation of rare but ecologically meaningful combinations, and ensuring that all systems were trained under the same experimental conditions (same input representation, training settings, and balanced label distribution). Although loss-function adjustments (e.g. weighted loss) can also address imbalance, such approaches are not straightforward for multi-label problems because segments may contain multiple events and class frequencies are not independent. Developing imbalance-aware loss strategies is an important direction for future work, but lies beyond the scope of the present study, which focuses on comparing system behaviour under controlled and reproducible conditions. We note that system performance and relative comparisons may vary when evaluated on the full imbalanced dataset.

The overall performance and stability of both BR and NDDR multi-task learning (NDDR MTL) systems for different numbers of concurrent event types, variability representation of label combinations in data and model sizes suggest they are better suited for real-world polyphonic analysis. BR’s success supports the utility of class-specific modelling in class-imbalanced, complex environments, while NDDR MTL’s competitive performance highlights the promise of SPS MTL systems. Together, these findings provide a foundation for MLC in marine bioacoustics, particularly in complex and ecologically significant situations.

Multi-labelled data

This study uses a subset of the PAM data collected during the COMPASS project (EU INTERREG). The data are collected over January, April, June and November in 2019 from a single site, Tolsta (58.392° N, −6.008° W) (Fig. 6). The COMPASS mooring comprises a single omni-directional broadband recorder (Ocean Instruments SoundTrap 300 HF) positioned in 100 m water depth⁵⁵. The recording protocol followed a 20/40 min on/off duty cycle at a sampling rate, ({f}_{s}), of 96 kHz.

Multi-labelling process

The MLC is based on the analysis of 3 s segments extracted from each 20-min recording. Five labels are considered: delphinid tonal, delphinid clicks, delphinid burst pulses, vessel noise and sonar, where sonar is used to represent anthropogenic, modulated, narrow-band signals, including sources such as echosounders and pingers. A subset of the available recordings was labelled based on visual inspection of the spectrograms with linear and mel frequency scales (2048 sample FFT size with 0.75 overlap) and audio playback, while also considering data outside of the 3 s segment to provide context. More than one label can be applied to a segment. There were several challenges faced when labelling which were resolved by using some basic criteria. Delphinid clicks were distinguished from other transient sounds by their short duration and based on only having energy above 5 kHz. In rare instances delphinid whistles and sonar could be challenging to separate, usually whistles could be distinguished by their higher degree of frequency and amplitude modulation.

All primary labelling was conducted by the lead author. To ensure reliability, labelled data was independently cross-checked by a second analyst. Any uncertain or ambiguous segments identified during the cross-check were excluded from the labelled dataset if there was no consensus between analysts; in total, fewer than 0.1% of all segments were removed through this process.

Labels were encoded as a 5-digit binary vectors for each 3 s segment, using the order [Delphinid Clicks, Delphinid Tonal, Delphinid Burst Pulses, Vessel Noise, Sonar], where ‘1’ indicates presence and ‘0’ absence of that class. For example, a vector of [1, 0, 1, 1, 0] would indicate the presence of clicks, burst pulses, and vessel noise, with no tonal and sonar. The binary vector [0, 0, 0, 0, 0] represents the ambient noise class.

A total of 101,535 segments (approximately 85 h of data) were annotated, yielding 192,209 individual labels (Table 2). The breakdown clearly shows that the dataset is dominated by vessel noise, which occurs roughly 2.5 times more frequently than the next most common event type, delphinid clicks. Several seasonal patterns are also apparent: delphinid vocalisations are most prevalent during the winter months (January and November); sonar detections increase in April, coinciding with an annual naval exercise in the region; and in November, no segments were labelled as ambient noise due to the near-constant presence of other classes, primarily those associated with delphinids. Four multi-labelled exemplars are given in Fig. 7.

Input representation

Systems benefit from using an input representation that is robust to changes in the background noise^23,26,57. Herein, we adopt per-channel energy normalisation (PCEN)^58,59 to provide dynamic range compression and stabilise signal energy levels to the Mel-spectrogram (MS), resulting in what we refer to as the MS-PCEN⁵⁵.

MS-PCEN is calculated for each 3 s audio segment based on a short-time Fourier transform (STFT) using a 2048-point (21 ms) Hanning window with 0.75 overlap. The MS is computed using a Mel-filter bank with 64 bands. PCEN is then applied as described in refs. ^58,59: The PCEN hyperparameters we used were: offset, (delta =0.05), gain parameter, (alpha =0.98), scaling factor, (r=0.5), numerical stability parameter (epsilon =1.4), and smoothing coefficient, ({rm{gamma }}=0.967) following⁵⁵. The resulting representation contains 64 Mel bands and 559 time frames for each 3 s segment. MS-PCEN is provided as a three-channel input by repeating the single-channel spectrogram. Although MS-PCEN can also be represented as a single-channel input, the three-channel format was retained to remain consistent with the formulation applied in previous work⁵⁵, yielding a final segment dimension of 64 × 559 × 3.

Systems background

Four different MLC systems are considered for analysing concurrent (polyphonic) acoustic data. Rather than focusing on the design or optimisation of new deep learning models, our primary objective is to examine how different MLC frameworks operate and perform. All the systems in this study share a common foundation based on a Convolutional Recurrent Neural Network (CRNN) architecture, which is widely used for similar tasks^27,28,37,60 and is detailed in the following.

In our study, CRNN architecture consists of three convolutional blocks followed by three recurrent blocks. Each convolutional block in the architecture comprises 2D convolutional layer followed by Batch Normalisation (BN)⁶¹. A Rectified Linear Unit (ReLU) activation function is used to introduce non-linearity⁶², 2D non-overlapping max pooling (MP) is used, followed by a dropout layer with a rate of 0.25 to help mitigate overfitting⁶³. We use 5 × 5 kernels across all convolutional layers, zero-padding the inputs to the convolutional layers and adopt pooling sizes of [(5, 1), (4, 1), (2, 1)] respectively. All the CNN outputs are reshaped by using time distributed flatten layers to ensure each timestep’s feature map is individually flattened for processing by the subsequent recurrent layers.

The recurrent blocks incorporate 128 bidirectional gated recurrent units (GRUs). Each GRU uses a tanh activation function, with a dropout rate of 0.2 applied to the input connections. The outputs of bidirectional GRUs are combined by using averaging⁶⁴. A temporal MP (TMP) layer is applied before the final dense output layer, which uses sigmoid activation to implement binary logistic regression function(s) (BLRF) for MLC task. This architecture is used uniformly across all systems.

System architectures

Here we provide details the four MLC architectures considered, these are illustrated in Fig. 8 and summarised in Table 3 with training and inference times. All the systems considered are based on networks of 4.1 M trainable parameters, the BR system employs K (=5) such networks, except in the System Model Size Effect experiment, where the number of trainable parameters is systematically varied.

The CMLC network shares a single CRNN structure across all labels. A dense layer is then used to produce the output vector, which contains one element for each class. The convolutional layers in CMLC use 160, 192 and 512 kernels. The network is trained to output predictions in the form of the K element binary label vector representing absences (0) or presences (1) for each of the K labels.

The BR system consists of K independent binary classifiers, each trained to detect a particular audio event. The output predictions of these K classifiers are concatenated to form the output vector. Each binary classifier in the BR system uses the same architecture used in the CMLC network but with a single output unit. Note that this means that the BR system has roughly K times as many parameters as in the corresponding CMLC network.

Instead of building separate models for each class (as in BR) or using completely shared network layers (as in CMLC), MTL architectures balance shared and task-specific components within a unified network. We implement two common types of MTL structures, which are examples of hard and SPS. The distinction between these approaches lies in the degree to which features are shared across classes during learning.

The HPS MTL system employs convolutional blocks that are shared across all classes, with distinct recurrent blocks used within each class. In our implementation, the shared convolutional layers have 96, 128 and 192 kernels, after which each class has its own bidirectional GRU layers.

Within the SPS MTL system, each class employs its own convolutional and recurrent blocks, but the blocks can exchange information within the network. We use 64 kernels across the convolutional layers to keep the size of the network similar to the other networks considered. For feature exchange, we employ NDDR layers with shortcuts in the convolutional blocks⁴⁵. At the end of each class-specific convolutional block, feature maps from all blocks are concatenated along the channel dimension and passed through a BN layer for stabilisation before being processed by the NDDR layers. The concatenation combines discriminative features from all classes, which results in feature maps with K times the number of channels compared to a single block. NDDR layers use 1 × 1 convolutional layers after BN layer to fuse features across blocks by also reducing the channel dimension to match a single output. ({l}_{2}) weight decay is applied to the 1 × 1 convolutional weights, with a moderate regularisation factor of 0.01 and each NDDR layer employs class-specific diagonal matrix initialisers to prioritise class-relevant features (values used are: 0.6 for self, 0.1 for others).

After each NDDR block, bilinear interpolation is used to ensure the outputs of all the layers are resized to match the spatial dimensions of the output of the final NDDR layer⁴⁵. For each class, resized outputs from the NDDR layers are concatenated along the channel dimension. Skip connections are constructed individually for each class, and the resultant feature maps, after concatenation, are subsequently reduced to the original channel dimension by an additional 1 × 1 convolutional layer with the same configuration as that NDDR layers, with the exception of diagonal initialiser. In this case, the initialiser assigns a weight of 0.6 to the final NDDR layer and 0.2 to the two earlier layers, which ensures that skip connections emphasise recent information while retaining contributions from prior layers. These convolutional layers are followed by class-specific recurrent layers to process each classification path accordingly. This SPS MTL structure with NDDR layers using skip connections is referred to as NDDR MTL in this study. Both HPS MTL and NDDR MTL systems produce K binary outputs, with each class-specific path terminating in a dense layer with one output unit. This design yields binary label vectors across the segments, same as the CMLC system.

Model training settings

A training batch size of 16 is used throughout. An early stopping criterion was not used in this study as it can lead to some system networks stopping earlier than others, which potentially causes incomplete training and inconsistencies due to the criterion’s dependence on validation (or test) error curves that frequently exhibit multiple local minima. In preliminary work, the error curves were monitored to determine a suitable number of epochs, the value used was 100, which was selected as it allowed all the networks to converge.

We use adaptive moment estimation (Adam) optimiser with exponential learning rate decay⁶⁵. The exponential decay schedule used an initial rate of 0.001 and a decay rate of 0.75 at every 90 steps with staircase function for stable training^26,50,55. The binary cross-entropy loss function is applied for all the systems. All the initialisation is Glorot uniform⁶⁶ for the weights and zero initialisation is used for bias for all the layers in the systems, except for that NDDR layers that use diagonal initialiser for the weights.

All experiments are implemented in Python using Keras library⁶⁷ and executed on the cloud environment Google Colab Pro+ framework with 52 GB system RAM and 15360 MiB of Tesla T4 GPU memory.

Data preparation and experimental setup

The imbalance in the labelled dataset (see ‘Multi-labelling process’) potentially causes significant challenges in system evaluation, as severe overrepresentation of a certain class can lead to biased training and testing. This imbalance is not readily addressed in MLC problems as it is with single class classification.

To mitigate this issue, we applied a systematic data arrangement strategy, for more details see Supplementary Material D. The aim of this process is not to eliminate imbalance entirely, but to substantially reduce extreme dominance effects (particularly vessel noise) and to provide the representation of rare event combinations while preserving seasonal trends. The process was carried out stepwise by prioritising underrepresented events and combinations. The initial steps targeted rare combinations, such as delphinid burst pulses co-occurring with sonar signals. Then, segments with sonar and burst pulses without vessel noise were selected to reduce the risk of over-representing vessel noise. Next, delphinid tonal and click signals were included, with priority given to segments that did not contain vessel noise.

Following this, the dataset was expanded to include segments with high co-occurrence rates among the events. This ensured that the added segments contributed to a balanced increase in the counts of multiple events. Ambient noise segments were also included to provide a realistic baseline and ensure that the dataset adequately represents the full soundscape. Although ambient noise is presented as one of six audio event classes in our datasets, it is not treated as a separate class during the training or inference, it is simply a record of when no class was detected. The curated dataset incorporates 15,520 audio segments, the distribution of the audio events within is shown in Fig. 9.

The curated data was randomly split into 4-folds. The number of training segments in each fold is 2904 ± 1, with 976 ± 1 segments in the validation data. This is arranged so that the percentage of labels in each fold and between the training and validation datasets is kept approximately consistent. Specifically the percentage of label in the ambient, delphinid tonal, clicks, burst pulse, vessel and sonar classes are, respectively, roughly 13.2%, 18.0%, 22.7%, 11.8%, 22.7% and 11.6% (the differences between the folds are <0.02% and the differences between the training and validation datasets are all <0.2%). For the number of audio segments audio events are present across the folds, see Table S3.

All training and validation experiments reported in this study are conducted in the curated dataset described above. The full data mentioned in ‘Multi-labelling process’ is not used directly for model training but is presented to characterise the original data distribution (Table 2) and to motivate the data arrangement strategy. Accordingly, the reported results correspond to validation performance, as no independent hold-out test set was used.

Evaluation metrics

In MLC systems, the system outputs are vectors of predicted probabilities across the labels. Evaluation is carried out by comparing predicted label vectors and target labels in the validation datasets. In MLC, a prediction can be partially correct, meaning it may correctly predict some, but not all, of the true labels. To measure the performance of MLC systems, we use accuracy and exact match⁴⁶.

Accuracy measures the proportion of correctly predicted individual class labels to the total number of labels in the label vectors across the audio segments.

where ({bf{y}}=left{{{bf{y}}}_{s}right}) is the set of target labels in the validation data comprising N segments, ({{bf{y}}}_{s}) is the binary vector for the ({s}^{{rm{th}}}) segment indicating the presence or absence of each of the K classes. ({y}_{s}^{left(lright)}) represents the ({l}^{{rm{th}}}) element of the label vector ({{bf{y}}}_{s}) corresponding to the label for the ({l}^{{rm{th}}}) class in segment s. The ‘hat’, as in ({hat{y}}_{s}^{left(lright)}), denotes predictions generated from the systems. The Iverson bracket (left[cdot right]) is the indicator function returning 1 if its argument is true and 0 otherwise.

An exact match defines the proportion of segments where the predictions for all the classes match the target labels exactly. It can be defined as:

The product over all classes ensures that the value is 1 only when all labels are predicted correctly for the segment; otherwise, it is 0. The exact match and accuracy can be related to each other if the accuracy is class independent, call that accuracy α and decisions for each class in a segment are assumed to be statistically independent, then the exact match score is given by ({alpha }^{K}). For consistency, a prediction probability of 0.5 is applied to the system outputs when making classification decisions for computing both accuracy and exact match scores.

Performance is also evaluated by segment-based precision and recall (PR) curves⁶⁸, which are computed by averaging PR curves across the folds, also the maxima and minima across the folds can be computed to assess variability. The area under the curve (AUC, aka average precision) for each class is also computed by averaging across the folds.