Listening for rain: Principal component analysis and linear discriminant analysis for broadband acoustic rainfall detection

Rainfall is an essential input to weather prediction models. More than 75% of the Earth's precipitation occurs over the ocean (e.g., Schmitt, 1995). Weather satellites can assist with rain detection and estimation at sea, but satellite measurements have large uncertainties (e.g., Tian and Peters-Lidard, 2010; Vinogradova and Ponte, 2017). These uncertainties reduce the accuracy of models that use satellite data. Furthermore, traditional rain gauges are impractical at sea, where equipment may be subject to turbulent surface conditions and human interference. For these reasons, meteorologists and climatologists seek alternate rain detection methods.

Rainfall at sea can be detected acoustically, using hydrophones. Rain falling on the ocean creates distinct acoustic spectra based on the physics of individual raindrops striking the water (e.g., Medwin , 1992; Nystuen, 1996; see also Fig. 1). The features of these spectra can be exploited by acoustic rain detection algorithms (Nystuen, 1986). Acoustic rain detection algorithms have been a topic of interest for decades. In the past, available computing power and storage space limited algorithm complexity, especially for algorithms intended to run on remote moorings.

Ma and Nystuen (2005) (MN05) proposed a rain and drizzle detection algorithm that is a decision tree comparing a few narrowband “discriminant frequencies.” The MN05 algorithm is computationally simple, and other recent papers (Pensieri , 2015; Yang , 2015) proposed similar decision tree rain detectors for different datasets still based on a few narrowband frequencies. However, these decision tree algorithms suffer from high rates of false alarms for a given detection rate. The MN05 algorithm is especially vulnerable in environments with high levels of noise within the discriminant frequency bands.

Recently, Trucco (2022b) (TR22) compared the performance of four common machine learning algorithms at detecting rain from broadband acoustic data. Their data set comprised of more than 18 000 h of acoustic data sampled at 100 kHz. The recording site was in the Mediterranean Sea and the hydrophone was at a depth of 36 m in water 1200 m deep and 80 km offshore. The data also includes hourly rainfall and average wind speed measured by a co-located meteorological station 10 m above sea level. The input to the machine learning algorithms were vectors of the 64 bin acoustic spectra with the same variable bandwidths above and below 3 kHz used in Ma and Nystuen (2005) without any dimension reduction preprocessing. However, TR22 analyzed 1 h average spectra, not the shorter records in prior rainfall detection studies. Additionally, the use of 1-h intervals affect the interpretation of rainfall rates, since “low” hourly average rainfall rates are more likely to contain bursts of high-intensity rain as compared to shorter intervals. Trucco (2023) notes that rain/wind predictions made over short intervals are significantly less accurate than predictions made using 1-h averages.

TR22 compared their results to six prior decision tree rain detectors for narrowband frequencies (Ma and Nystuen, 2005; Nystuen , 2015; Vagle , 1990). They found that Random Forest detectors performed best of the broadband methods (Detection Probability $P D = 70.8 %$ for False Alarm Probability $P F A = 1 %$ with an 0.1 mm/h minimum rainfall rate, and $P D = 92.1 %$ for $P F A = 1 %$ with a 1.0 mm/h minimum rainfall rate). They also found that all of the modern broadband algorithms outperformed the narrowband decision trees (or rule-based algorithms¹). The MN05 algorithm performed the best of the narrowband decision trees on the TR22 data set, with $P D = 52.1 %$ for $P F A = 1 %$ after compensating for hydrophone sensitivity bias. In their study, the MN05 algorithm performed only slightly worse than the linear discriminant analysis (LDA).

Like TR22, this study hypothesizes that the use of broadband spectral information makes rain detection more reliable. Additionally, this study investigates the effect of applying dimension reduction preprocessing to broadband acoustic spectra. The rainfall detector in this study operates on 3 min time windows in shallow water, in contrast to the one-hour averaging window of Trucco (2022b) in 1200 m water. The algorithm presented here has two main parts:

1. Principal Component Analysis (PCA), an algorithm that reduces the dimensionality of the broadband spectral information (Bishop, 2007, Chap. 12.1). It is in contrast to the frequency-bin reduction achieved by smoothing and decimation in TR22. Section III B contains details of the PCA.

2. LDA, a simple machine learning algorithm (Bishop, 2007, Chap. 4.1). LDA makes a binary decision of whether or not rain is falling in a given 3-minute window of time. LDA is advantageous in this experiment because it does not require much data to train. Section III C contains details of the LDA.

This PCA/LDA algorithm detects a larger fraction of rain events than the MN05 algorithm when both algorithms are tested against the same shallow-water acoustic data.² In this context, PCA does not appear to be significantly beneficial as compared to smoothing-and-decimation, but it is a valid method of reducing data dimensionality (see Sec. IV C). Section III describes the PCA/LDA algorithm in detail. Section IV gives the results of using PCA/LDA for rain detection, and compares them to the results of applying MN05's algorithm to the same input data.

The PCA/LDA algorithm is a binary hypothesis test for whether or not it is raining at a given time. LDA is a simple supervised machine learning algorithm. Like all supervised machine learning algorithms, LDA requires a labeled training data set that contains both input data and ground truth outputs. The trained LDA predicts the rain-state of a set of test data. The training and test data sets consist of simultaneously recorded acoustic data and weather data. The weather data provides ground truth. The experimental site is a pier in a saltwater inlet in New Bedford, MA (Fig. 2). The time period of the experiment is from August 6 to December 16, 2021. The PCA/LDA algorithm and the MN05 algorithm use the same data, so the experimental setup described here applies to both algorithms.

This section has three subsections: Collection and processing of acoustic data (Sec. II A), collection and classification of weather data (Sec. II B), and partitioning of the data into training and test sets (Sec. II C).

The acoustic recordings are continuous, except for brief maintenance periods. The data is stereo, but the rain detection experiment in this paper uses only the shallow 1 m hydrophone. The deeper hydrophone is only used to identify and reject transients and outliers.

A stereo Zoom H6 recorder (Zoom North America, Hauppauge, NY) recorded the underwater acoustic field measured by two HTI-96-MIN hydrophones (High Tech, Inc., Long Beach, MS), sampling at 44.1 kHz with an anti-aliasing lowpass filter at 22.05 kHz.³ The hydrophones are deployed at approximately 1 and 2 m below the low tide mark on a concrete piling supporting the pier at the UMass Dartmouth School of Marine Science and Technology. The mean low tide water depth at this location is approximately 3 m, with an average high tide depth of 4.3 m (resulting in a maximum water depth of about 2.3 and 3.3 m for the hydrophones).

The recorder produces.wav files up to 2.1 GB long (3 h 22 min for stereo data, with 16 bits per sample). The hydrophone runs on battery power to isolate the power supply from the seawater pumps located on the pier, reducing susceptibility to electrical noise.

The algorithm begins by considering data from the .wav files in 3-min windows. It splits each window into blocks of 2¹⁶ samples (⁠ $≈$ 1.5 s). The blocks have a 2-sided overlap of 10% and taper with a 10% Tukey taper (Blackman and Tukey, 1958). There are K = 127 blocks per window, with about half a block dropped at the end of each window. The algorithm computes the median acoustic power for the blocks in the window, then discards blocks with power greater than this median power. This discarding excludes loud transients. The time-domain data in the undiscarded blocks transform into frequency-domain power spectral densities (PSDs). There are $( K − 1 ) / 2 = 63$ PSDs per 3-min window.

The native frequency resolution of these PSDs (44.1 kHz/2¹⁶ FFT bins = 0.67 Hz/bin) is excessive for the requirements of rainfall detection.⁴ Applying the Blackman and Tukey (1958; Sec. B.17) frequency smoothing and decimation approach with a decimation factor of D = 32 reduces the dimensionality of the data while changing the frequency resolution to a more tractable 21.5 Hz/bin. This bin-width of 21.5 Hz/bin and decimation factor of D = 32 is used by the PCA. (To compare PCA to a pure smoothing-decimation method of dimensionality reduction, other values of D are used; see Sec. IV C.)

With D = 32, there are 1024 one-sided frequency bins, spanning 0–22.05 kHz. Highpassing at 1 kHz reduces the bin-count to F = 978 bins. These reduced-frequency-resolution PSDs also have excessive time resolution, about 1.5 s per PSD. The Bartlett (1948) periodogram averaging reduces the time resolution. After applying this procedure, there is one estimated PSD per window. Optionally, this PSD can be converted into physical units of dB re $μ Pa 2$/Hz, based on the sensitivity of the hydrophones.

Each estimated PSD is a “datapoint” when considered alongside the weather which occurs during the associated audio recording. A set of about 900 datapoints appears in Fig. 1.

Ground-truth weather data comes from a weather station on the pier (see Fig. 2). The weather station consists of a Rainlog-Rainew 2.0 tipping bucket rain gauge and a PortLog weather station unit.⁵ The Rainew is the primary source of rain data, and it provides rainfall rate estimates at 1-min intervals. Backup rain data is obtained from the PortLog at 10-min intervals. The PortLog also provides 10-min average windspeeds and air temperature measurements.  Appendix A describes how the raw weather data is smoothed to estimate rainfall rates and associated wind speeds.  Appendix A describes the weather estimation algorithm. The time-stamps of the estimated weather data are chosen to align with the time-stamps of the PSD estimates obtained in Sec. II A. The weather estimates serve as the ground-truth labels for the acoustic data. Together these PSD and weather estimates form a datapoint.

The weather estimates determine the classification of each datapoint for training and validation of the PCA/LDA algorithm. These classifications are used to partition the data into training and test sets. There are four rain classes: Dry (0 mm/h), Damp (0–1 mm/h), Light Rain (1–8 mm/h), and Heavy Rain (>8 mm/h). There are also three wind classes: Calm (<4 m/s), Breezy (4–8 m/s), and Heavy Wind (>8 m/s). A datapoint has both a rain and wind class. For example, a datapoint with a rainfall rate of 3 mm/h and a windspeed of 6 m/s would be classified as Light Rain + Breezy. Table I shows the counts of datapoints of each class in the experimental dataset.

The rain classes derive from the physical properties of raindrops striking the water. Rain generally contains a mix of drop sizes, which relates to the rainfall rate. Drops of different sizes create different sounds when they strike the water. Large, loud drops dominate Heavy Rain. These large drops have a relatively flat audio spectrum. Light Rain has a higher percentage of small drops. Small drops produce sound by a mechanism that depends on the angle at which the drop strikes the water. Consequently, wind has a large influence on the sound of Light Rain (Medwin , 1992; Nystuen, 1996). The “Damp” classification applies to precipitation that is below the chosen threshold of rain detection for this experiment (1 mm/h). Such precipitation mildly affects the audio spectrum, according to experimental observations.

Breaking waves generated by wind create their own characteristic PSD by stirring up waves on the water's surface, as studied by Knudsen (1948). However, the spectral interaction between rain and wind is complicated (Schwock and Abadi, 2021). The audio spectrum of windy rain depends on the rain rate, the amount of time the wind has been blowing, and (in a cove) the local topography. Whitecaps form on the water's surface at about 4 m/s (Wenz, 1962). This is the windspeed threshold between “Calm” and “Breezy.” The threshold between “Breezy” and “Heavy Wind” is more arbitrary.⁶ Sufficiently heavy wind causes a layer of bubbles to form beneath the surface of the water (Farmer and Lemon, 1984), attenuating sound from rain. Furthermore, wind-driven waves may affect the experimental site in unpredictable ways, given the shallow depth and scarcity of Heavy Wind datapoints. The windspeed threshold for “Heavy Wind” lies at 8 m/s, since wind over 8 m/s may make rain detection unreliable in this dataset.

A datapoint's rain class and wind class, taken together, are its weather class.

PCA/LDA requires the dataset to be partitioned into training and test sets. MN05 does not required training data in the same way PCA/LDA does but uses the training data to calibrate the levels. A given choice of data partition is called a “trial.” To estimate both means and standard deviations for the performance metrics, there are 100 trials in this experiment. PCA/LDA and MN05 both run on each of these trials independently.

The partition is constrained by representation and continuity requirements: The weather conditions in the test set must be well-represented by weather conditions in the training set, and the datapoints in each set must be reasonably contiguous in time. For these reasons, datapoints are grouped into “clips” of contiguous data prior to defining any trials. The partitioning process of each trial assigns clips of points to the training and test sets, rather than assigning individual points. Ultimately, each trial places about half of the usable datapoints in the test set and half in the training set, but the exact count of points in each set varies from trial to trial due to the variable length of each clip. For details of the clipping algorithm; see  Appendix B.

Excluding noisy datapoints is an accepted first step to analyzing passive acoustic data (Yang , 2015). A datapoint is “usable” if it meets certain conditions:

1. The air temperature must be above 0 °C, since the experiment is not intended to detect freezing precipitation.⁷

2. The point is not contaminated by intermittent biologic noise that can create false positives for rain detection. See  Appendix D, Item 2.

3. The point is not contaminated by intermittent pumphouse noise, which can also create false positives for rain detection. See  Appendix D, Item 4.

4. The point is not contaminated by an unknown tidally-correlated noise that appeared in late fall. This noise is identified by an algorithm that seeks a particular moiré-like pattern of interference at higher frequencies. See  Appendix D, Item 3.

Figure 1(c) marks usable and unusable datapoints. There are 37 344 usable datapoints in the total dataset, so a training or test set contains about 18 700 ± 900 datapoints. Table I shows the number of usable datapoints in each weather class.

The rain detection algorithm consists of two modular stages: PCA for dimension reduction, and LDA for binary hypothesis testing. Section III A gives a concise list of steps to implement PCA/LDA. For details on PCA, see Sec. III B. For details on LDA, see Sec. III C.

Steps 2–4 are part of the PCA. Steps 5–7 are part of the LDA. If running multiple trials of PCA/LDA, these steps should be implemented for each trial. The experiment in this paper uses 100 trials.

For each trial of PCA/LDA:

1. Partition the data into training and test sets, as described in Sec. II C.

2. Perform a singular value decomposition (SVD) on the spectral data in the training set. The result is a set of principal components (eigenvectors), which span the power spectra present in the training set.

3. Project the spectral data in the training set against the most important principal components found in Step 2. The result is a projection vector for each datapoint in the training set [Eq. (1)].

4. Project the spectral data in the test set against the most important principal components found in Step 2. The result is a vector of scalar values for each datapoint in the test set. (This is the last step in the PCA part of the algorithm.)

5. Use the projection vectors from Step 3 and the associated weather station data to train an LDA to detect rain.

6. Use the projection vectors from Step 4 as input to the LDA to attempt to detect rain at each datapoint in the test set. The result is a binary decision vector of about whether or not rain is present at each datapoint in the test set.

7. Compare the decision vector from Step 6 to the ground truth from the weather station. The result is a set of detection statistics for this trial.

After all the trials run, record the mean and standard deviation of the detection statistics over all trials. These are the performance metrics for the PCA/LDA rainfall detection algorithm (Table II).

This study investigates the use of PCA in reducing the dimensionality of spectral data. This section describes the PCA. For a comparison of PCA to smoothing-decimation, see Sec. IV C.

PCA expresses each spectrum as a weighted sum of a few eigenvectors (i.e., the “principal components”). Principal components are the left singular vectors of the training dataset X, where the columns of X are the PSDs obtained in Sec. II A, with each PSD representing 3 min of audio data. The training data matrix X is (⁠ $F × D$⁠), $D$ is the number of training datapoints, varying from trial to trial depending on the random partitioning of the data. F is the number of frequency bins remaining after highpass filtering the data at 1 kHz. In this study, F = 978.

Factoring $X = E Σ V †$ with an SVD yields the principal components as the columns of E (Golub and Reinsch, 1970). Sorting the eigenvalues on the diagonal of Σ in descending order also orders the columns of E in order of decreasing importance in explaining the data in X. The first N columns of E form an orthonormal basis for the N-dimensional subspace (or linear manifold) best representing the training data, and thus hopefully the test data as well. Each eigenvector can also be interpreted as a PSD for an uncorrelated acoustic source in the data (i.e., the principal components). This SVD factorization to find the principal components is Step 2 in Sec. III A.

Figure 3 shows a scatterplot of the first two eigenspace projections $c → 1 , c → 2$ for the entire dataset. Each point in the scatterplot represents a datapoint, with the colors used to represent weather data. Figure 3 also contains an LDA illustration, which will be discussed in Sec. III C.

Figure 4 shows the four most important eigenvectors of X, as determined by taking the mean and standard deviation of the $E [ : , 1 : 4 ]$ over 100 trial realizations of X. Most eigenvectors vary little from trial to trial; the translucent shading indicating the standard deviation is barely visible at this scale.

The first 2 eigenvectors in Fig. 4 have a noticeable relation to the physics of rain:

1. Eig 1 (⁠ $e → 1$⁠) is a mostly flat, positive-definite spectrum, representing the overall power of a PSD. Heavy rain creates loud PSDs, and therefore heavy rain produces high values of c₁. This can be seen in Fig. 3.

2. Eig 2 (⁠ $e → 2$⁠) is related to the acoustics of both rain and wind. According to Medwin (1992) and others, small raindrops are associated with greater power above 10 kHz, whereas wind is associated with greater power at lower frequencies (see Wenz, 1962, and others). Most rainy datapoints have positive projections against $e → 2$⁠. Meanwhile, most dry, windy datapoints have negative projections against $e → 2$⁠. Negative projections of $e → 2$ imply a boost in power at lower frequencies, such as that produced by wind. Figure 3 shows every datapoint's projection $c → 2$ against Eig 2.

Note that near-constant noise is not detrimental to the PCA/LDA's detection rates. For example, it is unnecessary to bandstop filter the spikes at 6.5 and 18 kHz. Bandstop filtering these spikes actually reduces rainfall detection rates, since it reduces the broadband information available to the algorithm. Likewise, the low-frequency pumphouse noise in Eigs 3 and 4 is part of a near-constant background noise. Highpassing the PSDs to cut this noise reduces PCA/LDA's detection rates.

Let Y be the test dataset, which contains the usable datapoints which are not in the training dataset X. Y projects into the eigenspace the same way as X, i.e., $C test = E [ : , 1 : N ] T Y$⁠. Projecting the test data is Step 4 of Sec. III A.

LDA is a simple, easily trained machine learning algorithm for binary hypothesis tests. LDA learns to detect rain from the projection coefficient matrix of the training set C [Eq. (1)], then is tested on the projection coefficient matrix of the test set C_test. The LDA follows Bishop (2007, Chap. 4). Training the LDA requires estimating both a weight vector $w →$ and a best decision threshold $Th$⁠. For a given datapoint with index d, the decision is implemented as $w → T C [ : , d ] ≷ Th$⁠. Interpreted geometrically, $w →$ and $Th$ divide the space of all projection vectors C into two classes (Rain and Dry) with an N – 1 dimensional hyperplane as a decision boundary between the classes. Figure 3 illustrates a notional decision boundary projected onto a line on the c₁ vs c₂ scatterplot.

The LDA sets a minimum rainfall rate as either 1.0 mm/h or 0.1 mm/h; below this rate, a datapoint is not considered “rainy.” When the minimum rate is 1.0 mm/h, rainfall below 0.4 mm/h is considered a dry datapoint,⁸ and when the minimum rate is 0.1 mm/h, rainfall below 0.05 mm/h is considered a dry datapoint. In both cases, there is a boundary zone between rainy and dry datapoints (0.4–1.0 and 0.05–0.10 mm/h, respectively). Datapoints in the boundary zone are excluded from the training and testing statistics. This exclusion zone clearly segregating dry and rainy data allows LDA to train better. Excluded points do not contribute to the calculation of detection or false alarm probabilities during either training or testing, but they pass through the LDA because this implementation of LDA uses median filtering (see Secs. III C 1 and III C 2, and  Appendix C). Note that the Dry/Rain/Excluded classification scheme used by the LDA is distinct from the rain classes used during the partitioning process (Secs. II B and II C). The list of classifications of all points in the training set is a vector $G →$⁠, and $G →$ is considered to be ground truth by the LDA. The equivalent ground truth vector for the test set is called $G → test$⁠.

This experiment trains the LDA both with and without supplementary data from the wind gauge. The inclusion of wind data improves algorithm performance (compare Tables II and III). However, PCA/LDA can function whether or not wind data is available.

When using wind gauge data, the LDA trains and tests on different wind classes separately. The LDA begins by partitioning C into three submatrices C^w, where each C^w contains the columns of C whose wind class is $w ∈${Calm, Breezy, Heavy Wind}. The LDA runs once for each wind class, using only the relevant C^w for training and testing. For example, the LDA uses only C^Breezy when training or testing rain detection in Breezy conditions. C^w is (⁠ $N × D w$⁠), where $D w$ is the count of datapoints in the training or test set with wind class w. Later, during testing, C_test will be partitioned the same way as C. The same logic applies to the ground truth classification vectors $G →$ and $G → test$⁠: When superscripted as $G → w$⁠, the w indicates a partition of $G →$ by wind class.

To avoid assuming a prior distribution on rain events, the LDA detector optimizes the Neyman-Pearson criterion for hypothesis tests, maximizing the probability of detection P_D subject to a constraint on false alarm probability P_FA (Van Trees, 1968, Sec. 2.2.1). Given that Table I indicates that only 3.2% of the measurements contain rainfall, this experiment chooses a false alarm tolerance of $P F A ⪅ 1$%, the same tolerance used in Trucco (2022b). Training the LDA for each wind class separately makes it difficult to set the overall false alarm rate precisely.

This subsection assumes that the LDA trains separately for each wind class w. As noted previously, it is also possible to ignore the wind classes and train with all of the data combined into a single wind class. The LDA's training algorithm produces 2 outputs for each w: The trained weight vector $w → w$ and the best threshold value $T h w$⁠.

Let $D rain w$ be the count of Rain datapoints in C^w, i.e., the count of datapoints d in $G →$ such that (⁠ $G d ∈ Rain$⁠). $D dry w$ is the count of Dry points. These datapoints have projections $C rain w$ and $C dry w$⁠, of size (⁠ $N × D rain w$⁠) and (⁠ $N × D dry w$⁠), respectively. $C rain w$ and $C dry w$ are mutually exclusive subsets of C^w. However, they are not partitions of C^w, because C^w will also contain Excluded points if $G →$ contains Excluded points.

The LDA training proceeds as follows for each wind class:

1. Take the means of $C rain w$ and $C dry w$ over all points in each wind class: That is, $⟨ C rain w ⟩ = ( 1 / D rain w ) ∑ d = 1 D rain w C rain , d w$⁠, and the same construction is used for $⟨ C dry w ⟩$⁠. Then create zero-mean versions of $C rain w$ and $C dry w$⁠, i.e., $C ¯ rain w = C rain w − ⟨ C rain w ⟩$⁠, with the same construction used for $C ¯ dry w$⁠.

2. Create a “whitening matrix” S^w for the combined zero-mean Rain and Dry data.

   $S w = C ¯ rain w ( C ¯ rain w ) T + C ¯ dry ( C ¯ dry w ) T .$
3. Create an $( N × 1 )$ “weight vector” $w → w$ from the whitened difference of the means of the rain classes:

   $w → w = ( S w ) − 1 ( ⟨ C rain w ⟩ − ⟨ C dry w ⟩ ) .$
4. Use $w →$ to create a ‘detection statistic’ $Z → w$ for all $D w$ points in the combined data set C^w (both Rain and Dry):

   $Z → w = ( w → w ) T C w .$

   (2)

   Optionally, $Z → w$ may then be smoothed by a continuous median filtering process to reduce false alarms (see  Appendix C). Figure 5 plots an example of a detection statistic (which has been smoothed). Like C^w and $G → w , Z → w$ may contain Excluded points, which are useful to the median filter. However, Excluded points are ignored by the LDA's performance metrics.

5. Find a “best threshold value” $T h w$ for $Z → w$ such that the binary hypothesis test,

   $H → w = Z → w Dry Rain ≷ T h w ,$

   (3)

   maximizes rain detection probability without exceeding the chosen false alarm tolerance FA_tol. Detection probability (P_D) and false alarm probability (P_FA) are calculated by comparing $H → w$ to $G → w$ at all indices d such that (⁠ $G d w ∈$ Rain XOR Dry). The best threshold $T h w$ should produce a binary hypothesis vector $H → w$ that yields $P D w = max ( P D w ) | P F A w ≤ F A tol$⁠.

There are multiple valid ways to seek $T h w$⁠. This experiment sought $T h w$ by looping through K = 100 possible threshold values $t h k w$⁠, in order to collect data on the P_FA vs P_D produced by different thresholds (i.e., the Receiver Operating Characteristic, ROC, of the rain detector; see Fig. 6). The values of $t h k w$ were distributed to match the nonlinear distribution of values in $Z → w$⁠, as determined by a nearest-neighbor estimation of the probability density of $Z → w$ (Bishop, 2007, Chap. 2.5).

Like $Z → w , H → w$ is $( D w × 1 )$⁠, and it may contain Excluded points. These Excluded points are used by an optional binary median filter on $H → w$⁠, which is applied prior to comparison with $G →$ to reduce false alarms (see the  Appendix C).

The LDA's performance as a rain detector is estimated by the testing phase. Testing the LDA should be done on each wind class separately if the LDA is trained on separate wind classes w. This subsection assumes separate wind classes. Testing uses the best threshold $T h w$ and weight vector $w → w$ found during training.

Some of the steps of the testing phase are analogous to steps in the training phase. However, it is not necessary to create a whitening matrix (Training Step 2) or weight vector (Training Step 3) during testing, since the information from these training steps is already incorporated in the weight vector $w → w$⁠.

1. Using the weight vectors $w → w$ from Training Step 3, create a detection statistic $Z test w$ for the test set $C test w$⁠, analogously to Training Step 4. If continuous median filtering was used on $Z → w$ during training, continuous median filtering should also be applied to $Z → test w$ immediately after this step.

2. Apply the hypothesis test from Training Step 5, using the best threshold value $T h w$⁠. The result is the binary vector $H → test w$⁠, describing a hypothesis for Rain/Dry for every point in $C test w$⁠. If binary median filtering was used on H^w during training, binary median filtering should also be applied to $H test w$ immediately after this step.

3. Merge the $H → test w$ for each wind class into one vector $H → test$⁠. There is no analogy to this step in training. An optional binary median filter is applied to this merged $H → test$⁠.

4. Compare $H → test$ to the ground truth logic vector for the test set, $G → test$⁠, at all datapoints d such that (⁠ $G test , d ∈$ Rain XOR Dry). This comparison yields the overall probability of detection (P_D) and false alarm (P_FA). A similar comparison can be used to calculate Precision, i.e., the probability that if the detector goes off, it is a real detection.

The LDA trains and tests independently for each trial in the experiment, giving an independent value for P_D, P_FA, and Precision for each trial. Table II of Sec. IV shows the mean and standard deviations of these statistics over 100 trials.

This section compares the performance of the PCA/LDA detector with the narrowband binary tree detector in MN05 (Ma and Nystuen, 2005), Sec. IV B. It also examines how the LDA performs when PCA is used, versus a pure smoothing-and-decimation approach to dimensionality reduction (Sec. IV C).

As described in Sec. III, the PCA/LDA algorithm bases decisions on 3-min PSDs, while the MN05 algorithm decides based on 4 × 10.24 ms PSDs spread across roughly 15 s. The first two metrics focus on rainfall events for each 3 min window. The first is a ROC quantifying the probability of detection (P_D) vs the probability of false alarm (P_FA) (Van Trees, 1968, Sec. 2.2). The probability of detection P_D is the fraction of the true rainfall events detected by the algorithm, while the probability of false alarms P_FA is the fraction of dry events that are erroneously reported as containing rainfall. Both P_D and P_FA vary with the detection threshold $T h w$⁠, so varying $T h w$ sweeps out a curve of the available P_D vs P_FA for the algorithm. The second performance metric was a comparison of the Precision and Recall for the rainfall detector (Bishop, 2007). Precision is the fraction of estimated detections that are true rainfall events, while Recall is the same as P_D. The third metric considers the total rainfall across the experiment window, i.e., what fraction of the total rainfall over four months was detected by the algorithm.

The PCA/LDA detector attempts to set the detection threshold $T h w$ to target a false alarm rate of approximately 1%. This choice for P_FA results from the MN05 implementation on the dataset, yielding a fair comparison of the algorithms.

Table II reports the performance for the PCA/LDA rainfall detector averaged across 100 Monte Carlo trials, with two different minimum rainfall rates (1.0 and 0.1 mm/h). As described previously, Calm and Breezy conditions use separate LDAs. Heavy Wind conditions are omitted from the experiment due to insufficient sample size of these conditions.

For comparison, Table III reports the detection results using 1.0 mm/h minimum rain rate and a single LDA for all wind conditions (including Heavy Wind). This simulates the performance of PCA/LDA if the detector does not have access to wind data. Lack of wind data degrades detection (and recall) performance by about 13%, but precision degrades by only 2%. The LDA without wind data also detects about 10% less of the total rain volume.

Processing the same dataset with MN05 establishes a benchmark for the performance of the PCA/LDA detector. Because the original MN05 algorithm assumed a deep water environment, the shallow water environment required some small modifications of the MN05 algorithm.⁹ The binary decision tree rules for narrowband discriminant frequencies of the original MN05 detector are unaltered. The performance reported for the MN05 implementation in this study used the same rainfall thresholds as the PCA/LDA results to allow a fair comparison: Dry $≤$ 0.4 mm/h, Rainy $≥$ 1.0 mm/h, and excluding 0.4–1.0 mm/h. These categories differ from the hard binary decision at 0.4 mm/h in the original paper (Ma and Nystuen, 2005).

MN05 refines prospective detections through three successive stages: initial detection, noise removal, and continuity check. The noise removal and continuity check stages reduce false alarms, but at the expense of decreasing the detection performance. The final results after the continuity check with $P F A ≈ 1 %$ form the primary basis for comparison with the PCA/LDA detector, as the new algorithm's threshold was chosen to match this false alarm rate.

Table V summarizes the detection performance of all three stages of the MN05 detector on the same datasets as processed by the PCA/LDA detector in Tables II and III. The narrowband MN05 algorithm only detects 31.5 ± 4.7% of the rain events with a $P F A = 1.0 ± 0.1 %$⁠. The fraction of the total rainfall detected is 51.6 ± 6.2%. Both by detection probability and fraction of total rainfall detected, the LDA/PCA detector performs better than the MN05 detector at a comparable false alarm rate.

To test the value of PCA, PCA/LDA's performance should be compared to the performance of an LDA detector that does not use PCA, but which has a similar dimensionality of data. Call this detector the High-Decimation LDA.

To create a High-Decimation LDA, the acoustic data of Sec. II A is reprocessed with a decimation factor of D = 512, instead of the D = 32 used in Sec. II A. This yields a set of acoustic PSDs X₆₄ with 64 one-sided frequency bins, or 61 bins after highpassing. These bins have a width of 344.5 Hz. Because X₆₄ already has low dimensionality, there is no need to use PCA to further reduce dimensionality. Therefore, X₆₄ is used to train the LDA directly. [By contrast, the PCA/LDA trains the LDA on the principal component projections C of a set of PSDs X with 1024 frequency bins and a bin width of 21.5 Hz. See Eq. (1) and Sec. III C]. Testing of the High-Decimation LDA proceeds similarly to training.

For the best comparison of High Decimation LDA to PCA/LDA, training and testing of the High-Decimation LDA uses the same data points as the PCA/LDA. Additionally, PCA/LDA is rerun with N = 61 principal components, so that the data-dimensionality of High Decimation LDA and PCA/LDA is identical. The results are shown in Table IV.

The performance metrics of the High Decimation LDA are similar to the performance metrics of PCA/LDA when data dimensionality is identical. Therefore, although PCA is a valid method of reducing data dimensionality, it is not significantly beneficial to a rainfall detection algorithm based on LDA.

The rain detection algorithm consists of two modular components, PCA and LDA. LDA is advantageous in this experiment because it does not require much data to train. LDA is the simplest linear binary hypothesis algorithm, yet the PCA/LDA algorithm has significantly better rain detection rates than the narrowband MN05 algorithm. This result indicates that the broadband acoustic data is valuable in detecting rainfall and that PCA is a valid method of capturing this data (although not the only valid method). Replacing LDA with a more-advanced algorithm may lead to even better detection rates.

Both broadband LDA detectors perform better than the MN05 narrowband decision tree algorithm on the same datasets, reinforcing the value of broadband acoustic data in detecting rain. The combined PCA/LDA detector yields a $P D = 72.6 ± 5.6 %$ for a minimum rain rate of 0.1 mm/h, significantly outperforming the standalone LDA results reported by TR22 (⁠ $P D = 58.3 ± 5.3 %$⁠), and slightly outperforming TR22's Random Forest method (⁠ $P D = 70.8 ± 5.4 %$⁠) for this minimum rain rate.

However, it would be naive to read too much into comparisons with TR22, since the datasets differ in almost every detail. Physical differences include water depth (3–4 m vs 1200 m), hydrophone depth (1 vs 36 m), and recording environment (urban coastal North Atlantic Ocean vs offshore Mediterranean Sea). The temporal differences are even more critical, namely, the use of 3 min windows instead of 1 h windows. As an example, a 1-h window with an average rainfall rate of 1.0 mm/h may correspond to a downpour of 20 mm/h in a 3-min windowing scheme, if the rest of the hour is dry. Likewise, a single 3-min window with a 1.0 mm/h rainfall rate will correspond to an hourly average rainfall rate of only 0.05 mm/h, if the rest of the hour is dry. TR22's results are therefore not directly comparable to the results of this study.

This study chose 3-min windows for datapoints due to the high levels of intermittent noise at the experimental site. In particular, motorboats passing the site were extremely loud for about 1.5 min at a time. By using 3-min windows, datapoints with particularly bad intermittent noise-contamination could be removed from the dataset, while quieter datapoints between the intermittent noises were allowed to remain. Acoustic rain detection in noisy areas may benefit from using short time windows such as those used in this study. The use of short windows does not entirely rule out the possibility of compounding quieter datapoints into longer intervals. A discussion of compounding strategies is beyond the scope of this paper but is featured in Trucco (2022a) in regard to acoustic wind prediction.

The PCA/LDA algorithm presented here uses N = 56 eigenvectors, representing 88% of the power of the PSDs in the dataset.¹⁰ Note, however, that this is 88% of the power, not 88% of the information which is useful for discriminating rainy points from dry points. The useful information content of each eigenvector is not known. Therefore, the choice of these 56 eigenvectors is somewhat arbitrary. Using more eigenvectors is not necessarily better, since the LDA must “learn” how rain affects all N eigenvectors. The optimal N for a given dataset depends on the amount of training data available.

False alarms tended to occur repeatedly at the same datapoints for different trials. Often, these false alarms occur adjacent to actual rain events (Fig. 7). This suggests that many “false alarms” may be due to two causes. First, there may be errors in the ground truth data: The rain gauge has 1-min time resolution and is subject to quantization error, so rainfall rates must be estimated rather than measured directly (see the  Appendix A). Second, the hydrophones may detect approaching or departing rain signals a few minutes before or after rain falls directly on the measurement site. This is a limitation of the equipment and experimental site: Namely, the experimental site is a shallow cove which can channel acoustic echoes from distant locations. By contrast, the rain gauge has a very small, precise collection area which can only detect local rain. Both the PCA/LDA algorithm and the MN05 implementation frequently suffered false alarms adjacent to actual rain events.

Wind plays an important role in the performance of acoustic rainfall detection. As noted in the results, training separate LDAs for calm and windy conditions yields significant improvements in performance. Unfortunately, the dataset contains few rain events during periods of Heavy Wind. The 78% rain detection rate in Table II only applies for average windspeeds up to 8 m/s.

The background noise levels are around 40 dB re $μ Pa 2 / Hz$ even during quiet periods (see Fig. 1). The experimental site is a shallow, rocky-bottomed urban cove; therefore, a high background noise level is reasonable.

Adaptability is an advantage of using a machine learning algorithm like PCA/LDA algorithm: The algorithm is capable of ignoring noise which is unrelated to rain, provided that the noise is well-represented in the training data.

Intermittent noise poses a more difficult problem since intermittent noise is not necessarily well-represented in the training set. Some intermittent noises can be removed by excluding loud outlying blocks of audio data (Sec. II A). For rainfall detectors placed in the open ocean, intermittent anthropogenic noise is less dominant than in urban littoral environments like the experimental site. Therefore, excluding loud anthropogenic outliers results in a better estimate of the algorithm's performance in the open ocean. Some intermittent noise is biologic, and intermittent biologic noise can be expected in the open ocean. Spatial filtering may be useful for removing some of these intermittent biologic noises ( Appendix D, Item 2), since acoustic rainfall detection relies upon noise originating from the surface of the water.

The PCA/LDA algorithm is modular, and either PCA or LDA may be replaced by other algorithms. For example, the LDA may be replaced by a decision algorithm based on support vector machines, as presented by the authors of this study (Berg , 2022).

This study found PCA to be of negligible benefit, as compared to a method of dimensionality reduction based purely on smoothing-and-decimation. This may be due to the adaptability of machine learning (ML) methods: ML learns important spectral features on its own, so intermediary processing steps like PCA may not have much to offer ML-based rainfall detection. Although this is a null result for PCA, awareness of null results may be useful in guiding future research. Additionally, although PCA was not found to be especially beneficial here, it was not detrimental. PCA may have utility when it is not possible to train an ML algorithm on site-specific data: The principal components of relatively clean training data may be more robust against the introduction of new site-specific noise than the smoothed-and-decimated-spectra. However, investigating the robustness of PCA across multiple sites was outside the scope of this study.

This paper presents a rain detection algorithm which exploits broadband acoustic data. PCA is used to reduce the dimensionality of this data. The primary principal components of the acoustic data relate to the physics of rain and wind interacting with the ocean. LDA decides whether rain is falling in a given 3-min window.

The PCA/LDA algorithm is a simple broadband detection algorithm, yet it performs better than the narrowband MN05 algorithm on the same data. The improved performance indicates that the broadband acoustic spectrum contains valuable information for detecting rainfall. While PCA is a valid means of extracting this broadband information, it does not offer significant performance gains for LDA in this experiment, as compared to broadband spectral data which has been smoothed and decimated to reduce its dimensionality.

This project was sponsored by the Department of the Navy, Office of Naval Research Grant No. N00014-20-1-2849. Alan Andonian designed and created the mounting bracket for the hydrophones. The authors wish to thank Forrest Kennedy and Jennifer L. Benson for their diving services installing the hydrophones.

The weather measurement equipment consists of a Rainlog-Rainew 2.0 tipping-bucket rain gauge and a separate PortLog weather station (both made by Rainwise, Boothwyn, PA, rainwise.com). The Rainew reports rainfall at 1-min intervals, and the PortLog reports both rainfall and windspeed at 10-min intervals. The minimum rainfall measurable by these devices is 0.254 mm (0.01 in) per reporting interval, corresponding to a single tip of the collection bucket. The rainfall measurements must be scaled to get an equivalent rate in mm/h (the standard units used in meteorology and prior acoustic rainfall work).

Both devices are subject to quantization error. However, the Rainew has a more pronounced quantization error due to its shorter reporting interval. Let $r raw [ t ]$ (mm/min) be the Rainew's raw estimate of the rainfall rate. The minimum nonzero value for $r raw [ t ]$ is 0.254 mm/min (15.24 mm/h), which is Heavy Rain. Periods of Light Rain will appear in $r raw [ t ]$ as a mix of zeros and 0.254 mm/min. By contrast, the Portlog's minimum nonzero rainfall rate is 1.524 mm/h, close to the minimum threshold of rain detection for this experiment. The Portlog's quantization error is therefore less concerning, especially since the Portlog's rain data were supplemental to the Rainew's data (as discussed later in this section).

To obtain reasonable estimates of Light Rain from the Rainew, it is necessary to use a smoothing algorithm to mitigate the quantization error. The smoothing algorithm used here is similar to the smoothing algorithm used in the original MN05 paper. Both the PCA/LDA and the implementation of MN05 presented here use the same smoothed rain data. The steps of smoothing $r raw [ t ]$ are as follows:

1. Convolve $r raw [ t ]$ with a Hanning window, $h [ t ]$⁠. The impulse response $h [ t ]$ is seven timesteps (minutes) long, inclusive of endpoint zeros, and is normalized so that $∑ t = 1 7 h [ t ] = 1$⁠. The convolution yields a smoothed rain time-series $r smooth [ t ] = r raw [ t ] ⋆ h [ t ]$⁠. Convolution with a symmetric 7-point impulse response introduces a constant group delay of 3 minutes between $r raw [ t ]$ and $r smooth [ t ]$ (Oppenheim , 1999, Sec. 5.1). In addition, the rain accumulation in the tipping bucket is effectively convolution of a theoretical instant rain rate with a 1-minute rectangular window, introducing an additional 0.5-minute group delay. The cascade of the bucket integration and smoothing creates a total delay of 3.5 minutes between the rain falling into the tipping bucket and the smoothed time series $r smooth [ t ]$ that must be accounted for when comparing the acoustic data with the meteorological ground truth.

2. Let $r ̂ ( τ )$ be the estimated rainfall rate in mm/h at a point τ in continuous time, where τ has units of minutes, and it is associated with the discrete time index $t$⁠. Then $r ̂ [ τ − ( 1 / 2 ) ] = 60 min/h × r smooth [ t + 3 ] mm/min$⁠, accounting for both the 3-minute delay caused by convolution and the 0.5 minute delay associated with the accumulation of rainfall within the tipping bucket.

The times $[ τ − ( 1 / 2 ) ]$ may not correspond precisely with the time-stamps of the acoustic PSDs used by the experiment. If necessary, interpolate values for $r ̂ ( τ )$ at the desired τ.

Smoothing with the Hanning window is effectively a lowpass filter. Extracting the rainfall estimates every 3 min out of the 1 min time series effectively downsamples the rainfall signal. Lowpass filtering followed by downsampling in time improves the signal-to-noise ratio of the rainfall data corrupted by quantization noise (Oppenheim , 1999, Sec. 4.9).

The accuracy of $r ̂ ( τ )$ is limited by both the amplitude quantization error and the time resolution of the Rainew. The Rainew recorder's clock does not store seconds, so the time-stamps of the Rainew data should only be regarded as accurate to within one minute. The Rainew recorder's clock was reset every few week when data were collected, and no drift in the time-stamps was observed over the collection period.

The PortLog system, which supplements rainfall data, operates on a separate clock not under the control of this experiment. Cross-correlating the rain data from the PortLog with the smoothed Rainew data found that the PortLog clock was advanced from the Rainew clock by 12 min. All weather data from the PortLog is corrected by shifting it back in time by 12 min. This 12 min shift is most relevant for the wind data since the PortLog is the only source of wind data. Rain data from the PortLog was only used directly when the Rainew data were unavailable, notably for 2 days in August 2021. The PortLog rain data is interpolated to the time-stamps of the acoustic PSDs without additional smoothing since the 10 min reporting interval of the PortLog is equivalent to smoothing minute-by-minute data by a 10-min rectangular window.

Windspeed data from the PortLog is the average windspeed of the past 10 min. This data is the best estimate of the instantaneous windspeed 5 min before the (corrected) time-stamp of the PortLog. This estimated windspeed is then interpolated to estimate the windspeeds at the time-stamps τ of the datapoints. Note that detailed wind data is unnecessary for this experiment. The windspeed estimate is only used to classify the wind state of a datapoint as Calm (< 4 m/s), Breezy (4 to 8 m/s), or Heavy Wind (> 8 m/s).

Datapoints must be combined into “clips” of temporally neighboring points prior to being partitioned into training and test sets at each trial (see Sec. II C). This clipping occurs only once, and its results are used for all trials. Clipping assists with satisfying two constraints of the training and test sets:

1. “Representation requirement”: The rain and wind conditions present in the test set should be well-represented by rain and wind conditions in the training set.

2. “Continuity requirement”: To the extent possible, datapoints in the training and test sets should be contiguous in time. This is because the presence or absence of rain is correlated from minute to minute, and both PCA/LDA and MN05 use this correlation to improve accuracy via median filtering (see  Appendix C, which is equivalent to MN05's continuity check). Additionally, contiguous datapoints give the most representative performance of both PCA/LDA and MN05, since both algorithms are intended for long deployments in the open ocean.

Requirement 1 implies that clips should be defined based on the weather present at one or more of the datapoints within the clip. Since the rarest weather does not occur in many datapoints, the clipping algorithm must prioritize defining clips based on datapoints with rare weather. The clipping algorithm begins by selecting datapoints with the rarest weather classes to serve as anchors (‘primary points’) for clips. A clip inherits its weather class from the weather class of its primary point. The weather class of a clip does not impact any process other than the partitioning process.

In this experiment, (Heavy Rain + Heavy Wind) was the rarest weather class, so the eight datapoints with this weather class were the first datapoints to be selected as primary points. Figure 8 shows a set of illustrative dummy data containing three primary points (A, B, C). Point A is a (Heavy Rain + Heavy Wind) datapoint, and so the algorithm selects Point A first. The algorithm will define a clip around A (that is, Clip A) before it selects any other primary points like B (Light Rain + Calm) or C (Dry + Breezy).

Primary points like Point A should be grouped with neighboring points, so that the clip contains a reasonable amount of chronologically contiguous data. This is done to meet Requirement 2. For this experiment, a clip can contain up to ten neighboring points on each side of the primary point. A neighboring point will be assigned to the clip if:

1. The neighbor point is usable by the experiment (see Sec. II C). In Fig. 8, unusable datapoints are marked in black.

2. The clip has not been “filled” with its allotted ten temporal neighbors. In Fig. 8, Clip A has been filled to the left of Point A: It contains its maximum of ten usable neighbors to the left.¹¹

3. The neighbor point is not more than about 10 min away from the nearest other point in the clip. A clip is allowed to have small time gaps in it (up to 3 datapoints) because some datapoints are unusable. In Fig. 8, Clip A contains a 1-point time-gap seven points to the left of Point A, but this gap is too small to truncate Clip A. However, Clip A is unable to expand more than five points to the right, because it is truncated by a time-gap of four unusable datapoints.

4. The neighbor point is not currently assigned to another clip. In Fig. 8, Clip A is completely defined before the algorithm begins to define a clip around Point C. Consequently, Clip C cannot expand more than 3 points to the right, because the points further to the right are already assigned to Clip A.

If two adjacent clips inherit the same weather class, the algorithm merges these clips into a single clip. This merging ensures that points from a brief extreme weather event are not assigned to separate clips. As a side effect of the merging, clips with common weather classes may be several hours long. The median timespan of a clip is 62 min. Therefore, this clipping procedure ensures that Requirement 2 is met.

After the algorithm defines clips around all the primary points in one weather class, it moves on to the next rarest weather class. In this way, rare weather classes are assured adequate representation in clips. Table VI shows the counts of clips by their inherited weather classes. These clips are defined only once; they do not vary from trial to trial.

After all clips are defined, they are available to be assigned to the training and test sets for each trial. The partitioning algorithm randomly assigns different clips to the training and test sets for each trial.¹² The partitioning algorithm assigns approximately 50% of the clips in each weather class to the test set, and 50% to the training set. In this way, Requirement 1 is met.

The LDA has modifications which are not described in the main text.

Median filtering is the PCA/LDA detector's version of MN05's continuity check. It reduces false alarm rates significantly, so the LDA always uses it (although technically, it is optional). The LDA implements median filtering two ways. The first median filter is continuous: It smooths the continuous detections statistic $Z → w$ (Step 4 in Sec. III C) with a central median filter with window length 3. The second median filter is binary: It processes the binary decision vector $H → w$ (Step 5 in Sec. III C). Each binary decision is set equal to the majority of the three decisions at that datapoint, the previous data point, and the next datapoint. This is effectively a sliding window version of the common M-of-N binary detector (Abraham, 2019, Sec. 10.3.1) for M = 2 and N = 3.

The median filters reduce false alarms significantly. They motivate the “continuity requirement” of  Appendix B. Data in the test set is median-filtered the same way as data in the training set, with one exception: If the LDA uses separate wind classes w, the binary median filter can be applied once more during the testing phase, after all $H → w$ have been merged into a single $H →$⁠.

To optimize performance of PCA/LDA, the LDA must be able to use median filtering on neighboring datapoints of different wind classes. In most regards, the LDA considers datapoints of different wind classes separately (Sec. III C), and uses a different detection statistic $Z → w$ and best-threshold $T h w$ for each wind class w. However, datapoints with different wind classes are interleaved in time, causing the data in each wind class to have many temporal “edges.” In order to avoid transient errors caused by these edges, a length-3 median filter must incorporate one point beyond each edge of each wind class w into $Z → w$⁠. These extra points outside the primary wind class w are called “shoulder points.” In the continuous median filter, shoulder points ensure the selection of the best threshold $T h w$ for a given $Z → w$⁠.¹³ In the binary median filter, shoulder points are used on the binary decision vector $H → w$ for each wind class in training and testing. Shoulder points are not necessary when median-filtering the decision vector $H →$ after merging the wind classes.

The PCA/LDA algorithm can tolerate continuous, self-similar noise sources, but does not work well with loud intermittent noises that are common at the experimental site (a shallow urban cove). To the extent feasible, a cleanup algorithm removes datapoints that are contaminated by these intermittent noises, prior to PCA/LDA analysis or MN05 analysis.

There are four significant sources of intermittent noise in the raw data:

1. Brief, unpredictable transient noise is mitigated by the use of order statistics. See Sec. II A.

   Additionally, for noise that is closer to one hydrophone than the other, there may be a poor correlation between the power spectra of the two hydrophones. Therefore, times with poor power correlation are marked as noisy by the cleanup algorithm, and excluded from the experiment.

2. There is nightly biologic noise in summer which appears to be from cusk eels, a common species in northeastern coastal waters. This noise can be automatically identified and excluded by using the two hydrophones as a simple direction-finding array (Fig. 9).

   The two hydrophones are located 1 m apart, corresponding to a half-wavelength of approximately 765 Hz in summer seawater. The eels are very loud at this frequency. When the eels are not active, most of the sound at 765 Hz appears to originate directly above or below the hydrophones, indicating either reflections or sound conduction through the leg of the pier. However, when the eels are active, much of this sound originates in the middle of the water column where the eels swim.

   The cleanup algorithm marks a datapoint as cusk-eel-contaminated if the direction-of-arrival of 765 Hz sound has a large horizontal component on a summer evening (when the cusk eels are most active). This is a form of spatial filtering.

3. Fall data sometimes included a “moiré” spectrogram pattern (Fig. 10) of unknown cause. The moiré pattern has distinct ripples between 16 and 22 kHz. These ripples result in the PSD having high variance between these frequencies. The moiré can be identified by looking at the variance of each PSD across these frequencies. The cleanup algorithm assumes that PSD datapoints with high variance in this frequency range are contaminated by moiré.

4. The pumphouse machinery is more active at some times than others (Fig. 10), and was particularly active during the autumn months. This activity is not dictated by weather conditions.

   The high-activity cycle of the pumphouse is common, loud, and has consistent spectral characteristics at lower frequencies. The cleanup algorithm identifies times of high pumphouse activity by seeking datapoints whose lower PSD frequencies are highly correlated with other datapoints, or datapoints which contain particular spectral patterns at these frequencies. The cleanup algorithm excludes these datapoints.