A functional regression analysis of vessel source level measurements from the Enhancing Cetacean Habitat and Observation (ECHO) database

The international shipping lanes in the coastal waters of southern British Columbia, Canada directly overlap with the critical habitat of the endangered Southern Resident killer whales, who use sound to navigate, communicate, and, most importantly, echolocate to find prey. Fisheries and Oceans Canada identified underwater noise from shipping as a key threat to the recovery of this endangered species (DFO, 2017). Indeed, vessel noise is often the main driver of acoustic habitat quality within coastal areas of importance to sensitive marine mammal species (Erbe et al., 2012; Williams et al., 2013; Erbe et al., 2019). In 2014, the Enhancing Cetacean Habitat and Observation (ECHO) program was established by the Vancouver Fraser Port Authority to better understand and reduce the effects of shipping on endangered whales with a focus on underwater noise (VFPA, 2016). In conjunction with partners at Transport Canada (Ottawa, ON), Ocean Networks Canada (Victoria, BC), and JASCO Applied Sciences (Victoria, BC), a project was initiated in 2015 to develop a real-time method for measuring and reporting vessel source levels for deep sea vessels transiting to the Port of Vancouver (Hannay et al., 2016).

Between September 2015 and January 2020, more than 20 000 vessel source level measurements for vessels of opportunity were collected on behalf of the ECHO program, referred to here as the ECHO database. These data were collected using three different hydrophone stations (Fig. 1), approximately conforming to the ANSI (2014) and ISO (2016) standards for vessel radiated noise level (RNL) and source level measurements. A dataset of this size provides an unprecedented opportunity for statistical investigation of correlations between underwater noise emissions measurements and readily available vessel design and operational characteristics.

A previous study (MacGillivray et al., 2019) identified that slowing vessel speed is an effective method for reducing underwater radiated noise from commercial vessels. Vessel slowdowns can provide an immediate benefit to underwater noise reduction but may not always be safe and feasible. The key to reducing overall underwater noise for the Southern Resident killer whale and other aquatic species worldwide will be reducing vessel noise at the source, likely through a combination of operational and design measures (Southall et al., 2017). This study was undertaken to investigate correlations between design and operational characteristics to develop a statistical model to predict underwater radiated noise and identify the key drivers most affecting noise predictions. Understanding these relationships can inform or enhance underwater noise modelling efforts and assist in designing future noise reduction efforts.

Although engineering controls for vessel underwater radiated noise are well developed (Spence and Fischer, 2017), the reasons for the large observed differences in underwater radiated noise between similar types of vessels remain unclear (Leaper and Renilson, 2012). While it has been proposed that quieting efforts should be focused on the noisiest vessels, the first step must be to identify which design and operating characteristics are associated with the noisiest vessels (Leaper et al., 2014). Recent European research programs have highlighted the need for more noise measurements of commercial vessels to fill current data gaps (Brooker and Humphrey, 2015; Audoly et al., 2016). Similarly, a recent horizon scanning exercise found agreement between researchers, policy makers, and other stakeholders that research into the characteristics of marine vessels which are associated with underwater noise should be prioritized (McWhinnie et al., 2017). The current study was undertaken in an effort to help fill in these data gaps using the large ECHO source level dataset.

This article describes the development of a statistical model that was used to associate vessel design and operating characteristics with measurements from the ECHO database. The statistical model was used to identify how frequency-dependent source levels varied with nine predictor variables in six different categories: cruise ships, container ships, bulk carriers and general cargo ships, tankers, tugs, and vehicle carriers. Section II describes the data collection and analysis methods and the development of the statistical model. Section III describes the results of applying the statistical model to the ECHO dataset. Section IV includes a discussion of which model parameters were most influential, sources of uncertainty, comparisons with prior work, and a discussion of remaining data gaps. Finally, Sec. V provides a summary of the main conclusions of this work.

Understanding underwater radiated noise from marine vessels is a priority study area for the ECHO program. Since September 2015, the ECHO program has accumulated a comprehensive database of source level measurements for vessels of opportunity transiting along the international shipping lanes near the Port of Vancouver (Canada). Source level data were collected at three sites within the Salish Sea (Fig. 1) in water depths ranging from 173 to 250 m using a combination of cabled hydrophone nodes and autonomous seabed recorders. The ECHO database was automatically collected and analyzed using JASCO's ShipSound™ system (Victoria, BC) and includes measurements of many commercial vessel types, which were identified through correlation of acoustic measurements to the automated identification system [AIS; by Maritime Mobile Service Identity (MMSI) or International Maritime Organization (IMO) number]. Each vessel measurement in the ECHO database includes the monopole source level (MSL), RNL, and associated metadata, including closest point of approach (CPA) to the station, speed through water (STW), vessel draft, navigational information, surface current velocity, and weather parameters. The measurements analyzed for the present work were collected from September 2015 through January 2020.

The MSL and RNL were calculated in standard decidecade bands from 10 to 63 000 Hz. The RNL was calculated assuming spherical spreading propagation loss (PL; i.e., $PL=20×log10r$⁠) as per the ANSI standard, whereas the MSL was calculated using a frequency-dependent PL model based on the numerical solution of the acoustic wave equation, which accounts for the effect of the environment on sound transmission. At 5 kHz and below, PL was calculated using a split-step Padé parabolic equation model (Collins, 1993), modified to treat shear wave reflection losses in the seabed (Zhang and Tindle, 1995). Above 5 kHz, PL was calculated using an image-method model (Brekhovskikh and Lysanov, 2003), which accounts for surface and seabed reflection coefficients and frequency-dependent absorption (François and Garrison, 1982). The RNL and MSL metrics incorporated seawater absorption into the PL calculation above 5 kHz, which is not specified in the ANSI standard for calculating the RNL. The incorporation of absorption into these calculations was necessary because vessel CPA distances often differed from the ANSI S12.64 (ANSI, 2014) specified distance. The average PL in each decidecade band was based on the mean propagation factor calculated at 50 frequencies, which were spaced logarithmically between the minimum and maximum band limits.

Calculation of the MSL depends on the choice of the source depth, which was taken to be half the actual vessel draft at the time of measurement. The PL was smoothed by assuming the source depth had a Gaussian distribution in a manner similar to that used by Wales and Heitmeyer (2002), where the standard deviation was taken to be 30% of the mean source depth. It should be noted that the value of the MSL is sensitive to the assumed monopole source depth at low frequencies, and the monopole source depth was not a directly measured dimension or otherwise determined quantity (although recent work by Tollefsen and Dosso, 2020, suggests that it may be estimated using arrays of hydrophones and matched-field inversion). It should be noted that the ISO 17208–2 standard for vessel source level measurement in deep water (ISO, 2019a,b) assumes that the source depth is equal to 70% of the vessel draft rather than the 50% used here, but this standard was published after data collection for the ECHO program began. Differences in assumed monopole source depth may lead to systematic differences in reported source levels between studies.

While the RNL and MSL measurements are collectively referred to here as source levels, the RNL is actually an affected source level (i.e., a source level measurement affected by surface and seabed reflections). Only the MSL strictly corresponds to the ISO standard definition of a source level (ISO, 2017). More details regarding the methods used by ShipSound for calculating source levels of marine vessels are given in Hannay et al. (2016) and MacGillivray et al. (2019).

This research project was limited to commercial vessels in the following six categories:

• bulk carriers and general cargo ships (hereafter referred to as “bulkers”),

• container ships,

• cruise ships (i.e., passenger vessels greater than 100 m in length, excluding ferries),

• tankers,

• tugs, and

• vehicle carriers.

These six categories were based on the vessel types recorded by the local pilotage authority and are used by the ECHO program for their targeted underwater noise reduction initiatives (Burnham et al., 2021; Trounce et al., 2021). Previous studies of vessel source level measurements (MacGillivray et al., 2019) indicated clear differences in the operational and noise profiles of these different vessel categories. As such, these categories were maintained as separate in the statistical analysis with the intent to develop a robust source level prediction model for each category. These categories were also consistent with those used by a previous vessel underwater radiated noise study in the same region (Veirs et al., 2016).

Prior to the inception of the this study, the ECHO program convened a meeting of their Acoustic Technical Committee (ATC), an advisory group of acousticians, bio-acousticians, naval architects, and engineers, to identify how their database of source level measurements could be used to better understand factors driving underwater radiated noise from vessels (VFPA, 2018). The goal of the ATC meeting was to identify those design characteristics that were expected to be correlated with vessel underwater radiated noise and recommend methods for investigating their influence. Following the meeting, the ECHO program sought additional information on those vessel design characteristics recommended by the ATC, although it should be noted that many design characteristics (particularly those related to propeller design and resilient mounting of machinery) were unobtainable from publicly available sources.

Records from the ECHO database were joined to two additional databases to provide information on vessel design characteristics and operating conditions for the noise correlations study:

• a database of vessel design characteristics provided by Lloyd's List Intelligence (London, UK), including variables describing dimensions, displacement, propulsion, and nominal operating conditions, and hereafter referred to as the LLI database; and

• a database of transit logs providing records of actual vessel draft at the time of transit from the Pacific Pilotage Authority (Vancouver, BC) as recorded by on-duty pilots, and hereafter referred to as the PPA database.

Each measurement in the ECHO database was matched to records from the LLI and PPA databases based on the IMO number whenever possible. The IMO number is a seven-digit code that uniquely identifies large cargo vessels (>300 gross tons) and large passenger vessels (>100 gross tons). In cases where an IMO number was unavailable or was recorded incorrectly, records were instead matched on the basis of MMSI or by vessel name. IMO numbers, MMSI numbers, and vessel names in the ECHO database were obtained from the AIS as broadcast at the time of measurement. The data from all three databases were merged into a single vessel noise database for subsequent analysis. Vessel subtypes from the LLI database were cross-checked against the six ECHO categories to verify that the ship classifications were consistent between the different databases.¹

Only measurements that passed a manual quality review were retained for subsequent statistical analysis. Vessels that could not be matched to an entry in the LLI database were excluded from the analysis. A total of 9880 accepted measurements of 3188 unique vessels met these criteria and were retained in the merged vessel noise database.

In developing the statistical model, predictor variables were taken to be those variables that may have influenced underwater noise emissions. Two types of predictor variables were considered in this study:

• Design: variables related to the design characteristics of a marine vessel. These variables are assumed to remain constant between repeat measurements of the same vessel; and

• Operational: variables related to the operational characteristics of a marine vessel. These variables may change between repeat measurements of the same vessel.

A set of suitable candidate predictors were identified from all of the data sources captured in the merged vessel noise database (Table I). These candidate predictors were evaluated during the model development to determine which variables should be retained for the final statistical model (see Sec. II E).

A manual review was performed to clean the merged database and remove outlier values. Invalid data were corrected when possible (e.g., using online databases) or flagged as missing if they could not be corrected. A small subset of the dataset contained sonar-like signals (e.g., from echosounders) in the 16–31.5 kHz frequency range. Source levels in frequency bands affected by this issue were flagged as missing data [i.e., value set to not available (NA)] because none of the available design data in the LLI database pertained to sonar equipment and, furthermore, the measurement geometry was not designed to sample narrow-beam sound emissions from hull-mounted sonar transducers. The cleaned version of the merged database was used for all of the subsequent statistical analyses.

The merged database contained missing (NA) values where information for a specified predictor was unavailable in the LLI and PPA databases. The percentage of missing data was calculated for each of the candidate predictors, broken down by the vessel category (Table II). Source levels were treated as missing (NA) when background noise levels were within 3 dB of received signal levels during a vessel measurement, following the ANSI S12.64 standard (ANSI, 2014). Additionally, the source levels for the 40, 50, and 63.1 kHz decidecade bands were missing for a subset of the measurements collected in the Georgia Strait location as the sampling rate was lower at this location (64 kHz). The missingness (or absence) was generally greatest at the lowest and highest frequencies for all of the vessel categories. Imputation procedures, based on non-missing data for similar measurements, were used to fill in missing predictor and source level values for the functional regression analysis.¹

Physical models were used to capture the effect of water currents, wind, and source-receiver geometry on measured source levels. The magnitudes and directions of wind and water currents were captured by ShipSound at the time of measurement and stored in the ECHO database. The windspeed and direction data were obtained from nearby Environment Canada weather stations in Haro Strait, Georgia Strait, and Boundary Pass (Environment Canada, 2020). Where possible, ocean current data were obtained from a combination of acoustic Doppler current profiler (ADCP) measurements collocated with the underwater listening station (ULS) nodes in Georgia Strait and Boundary Pass. Where direct ocean current measurements were unavailable (e.g., for the Haro Strait hydrophones), ocean current data at the Haro Strait ULS were obtained from the WebTide tidal prediction model (version 0.7.1), provided by Fisheries and Oceans Canada (Bedford Institute of Oceanography, 2015).

For each measurement, the STW vector was computed as the difference between the speed over the ground vector (from AIS) and ocean current vector. STW was used in the statistical analysis to capture in a single predictor the effect of the water current and vessel speed on underwater radiated noise. Thus, the magnitude of the resulting STW vector was used as a predictor in the subsequent statistical analysis. The speeds over the ground and ocean currents were not used as predictors because they are implicitly included in the STW calculation. The effect of wind on vessel source levels was captured using a dimensionless wind resistance factor, K_w, which depended on the speed and course over the ground of the vessel and magnitude and direction of the wind at the time of measurement (see the  Appendix).

While the measurement procedure was established to generate repeatable source levels and minimize the influence of source-receiver geometry on measured source levels, differences in the source-receiver geometry between measurements may nonetheless have a residual influence on the measured source levels. Such influence comes about naturally in the RNL metric, which does not explicitly account for the effect of the environment on each measurement. Such influence may also come about in the MSL metric due to imperfect characterization of the environment in the PL model or mismatch in the actual versus assumed source depth. To control for these residual effects, the surface grazing angle was included as a predictor in the statistical analysis to control for differences in sampling geometry between measurements and sites. The surface grazing angle was calculated from the horizontal distance of the CPA of the vessel (⁠$x0$⁠) and depth of the hydrophone (h) such that

Besides capturing $x0$ and h in a single statistical predictor, another advantage of using the surface grazing angle is that it also captures residual cancellation of radiated sound by the sea-surface (i.e., the Lloyd-mirror effect), which is not fully accounted for in the RNL calculation.

Three dimensionless derived quantities were also calculated from groups of predictors in the merged vessel noise database and evaluated for predictive significance. The first derived quantity was the block coefficient, which was the ratio of the volume of displacement (⁠$V, in cubic metres)$ to the product of the breadth (⁠$B$⁠), length overall (⁠$L$⁠), and draft (⁠$d$⁠),

Note that the block coefficient could only be calculated for the summer draft (i.e., as a static value) because the true displacement depends on the actual draft in a manner that depends on the hull design. The second derived quantity was a dimensionless speed, which was the ratio of the actual STW (⁠$v$⁠) to the design speed of the vessel (⁠$v0)$⁠, such that

The third derived quantity was the dimensionless draft, which was the ratio of the actual draft (d) to the summer draft (⁠$d0$⁠),

These three derived quantities were evaluated using a univariate analysis to determine whether they were more strongly correlated with vessel source levels than the constituent quantities on their own.

Many vessels in the ECHO vessel noise database had more than 1 measurement, with 161 vessels (of the 3188 total) having 10 or more measurements. Repeat measurements are valuable when they capture the same vessel under different operating conditions, but they can also bias the analysis by weighting the result toward the most frequently sampled vessels. To balance these competing effects, repeat vessel measurements were randomly subsampled (without replacement) so that they were included only when the following three operating conditions were substantially different: STW, actual draft, and wind resistance. To perform the subsampling, repeat measurements of a vessel were binned according to these three variables and a single measurement was randomly selected from each bin, up to a maximum of eight randomly selected measurements per vessel. The following bin widths were used for the subsampling procedure:

• STW: bin width equal to 20% of vessel design speed;

• actual draft: bin width equal to 20% of vessel summer draft; and

• wind resistance: bin width equal to 100 = (10 m/s windspeed)²/(1 m²/s²).

The subsampling ensured that different operating conditions were captured in the statistical analysis without biasing the result too heavily to any single vessel.

Exploratory data analysis methods were used to examine the interrelationships among the ship characteristics and operational parameters and noise across frequency bands. Bivariate scatterplots (i.e., X-Y plots) and correlation coefficients were created for pairs of numerical variables to identify trends and outliers. Density plots (i.e., smoothed histograms) were used to assess the distributions of numerical variables. These exploratory plots were also used for identifying which variables should be transformed for subsequent regression analysis. A logarithmic transformation was identified as the most suitable transformation for most predictor variables because source levels also measure radiated noise on a logarithmic (i.e., decibel) scale.

Not all of the predictors in the merged vessel noise database were retained for the regression analyses. Variables with a large amount of missing data or redundant variables were removed. Including redundant predictor variables causes collinearity that can lead to numerical instability and inaccurate regression coefficient estimates. Additionally, fitting too many predictors can lead to overfitting the data and reducing the generalizability of the statistical model. Some predictors were not immediately relevant to underwater radiated noise, and others did not exhibit a strong correlation with measured source levels, hence, they were dropped from the statistical model.

Variable selection for developing the final model was performed using a combination of statistical analysis and expert knowledge. Functional regression analysis (see Sec. II F) was used to investigate linear trends between individual (transformed) variables and source levels across decidecade bands. Those variables that had weak correlations were omitted, particularly if another variable captured a closely related design characteristic of a vessel. A short list of design predictors was reviewed in consultation with a team of subject-matter experts to identify which ones were to be retained, developing a functional regression model with multiple predictors.

The final functional regression models simultaneously analyzed the correlation of the subset of the predictors with vessel source levels. Functional regression models were developed in a forward-stepwise fashion by incrementally adding predictors to the model and evaluating the coefficient of determination (R²) versus frequency. Operational parameters were added first, followed by design parameters. Predictors were retained in the final model if they increased the R² value over a range of frequencies.

Functional data analysis is a modern statistical technique (Ramsay and Silverman, 2005) used to analyze the characteristics of curves or profiles. Functional regression analysis provides a means for carrying out regression analysis when the outcome or predictor variable values are curves rather than individual data points. For example, Ainsworth et al. (2011) used functional data analysis techniques to relate daily river flow patterns to annual salmon returns. This approach identified the river flow patterns and seasons with the strongest relationships to salmon return rates. In the present context, functional regression analysis is used to model how source level patterns across sound frequencies are affected by ship characteristics and operational parameter variables.

The advantage of functional data analysis is that it captures the information in the entire curve (i.e., the entire source level versus frequency profile). Rather than looking at individual correlations of ship characteristics and operational parameter variables with source levels in distinct frequency ranges, functional regression analysis simultaneously models the relationship of ship characteristics with the entire source level curve. Thus, functional regression models identify which ship characteristics and operational parameters best predict noise and how these relationships change across frequency bands.

Functional regression analysis is an extension of standard regression analysis. For each observation, the outcome variable value (or predictor variable values) can be a curve rather than a single number. Standard regression analysis leads to a single regression coefficient for each predictor. On the other hand, functional regression analysis avoids the need to run multiple regression analyses on noise levels aggregated across frequency bands. It provides simultaneous estimates in the form of a regression coefficient function for each predictor variable. This function indicates which frequencies are most correlated with a predictor variable and the direction of the relationship across all frequencies. For example, a predictor variable may have a positive relationship with source levels at low frequencies, no relationship with source levels at mid frequencies, and a negative relationship with source levels at high frequencies.

Similar to standard regression analysis, functional regression analysis models can include a single predictor or multiple-predictor variables (single-variable and multiple-variable cases). Multiple-predictor functional regression models show the relationship between predictors and vessel source levels (MSL and RNL) in decidecade bands.

The functional regression model with multiple covariates is

where $yi(f)$ is the source level for observation i at frequency f, $α(f)$ is the intercept, x_i is the vector of covariate values for observation i, β(f) is the vector of regression slopes, one for each covariate, and ε_i(f) is the error for observation i. The noise and regression slopes are modelled as smooth functions across the frequency range. Therefore, the slope estimates can be viewed as a curve where the x axis is frequency and the y axis is the regression slope. Therefore, the regression estimates can be viewed as a curve for each covariate with frequency on the x axis and slope on the y axis. These curves show the direction and strength of each covariate across all frequencies.

Functional regression modelling was performed using the R programming language (version 4.0.5; R Core Team, 2020) using the fRegress function in the fda package (version 5.1.9; Ramsay et al., 2020). Smoothing was applied to the curves using a b-spline basis evaluated on the range of observed frequencies. Lambda, which specifies the amount of smoothing, was chosen to minimize the average generalized cross-validation (GCV) criterion.

Plots of the decidecade band measurements showed that the average source level versus frequency curves exhibited a similar pattern for all six vessel categories (Fig. 2): noise emissions below 100 Hz were dominated by a broad hump, centered around 50 Hz, and noise emissions above 100 Hz steadily decayed with frequency at a rate of approximately −8.2 dB per decade. Average noise emissions were generally the highest for container vessels, which was consistent with previous observations by McKenna et al. (2012). Average noise emissions were generally lowest for cruise vessels. Although the average source level trend was smooth for each category, there was a significant amount of variability between measurements, particularly at the lower and higher ends of the measured frequency range. Furthermore, many individual measurements featured strong spikes, corresponding to narrowband or tonal features in their source level spectrum. It is interesting to note that similar narrowband tones also were observed in measurements of container vessels in Santa Barbara Channel by McKenna et al. (2013).

Several notable observations resulted from the exploratory analysis of the merged vessel noise database:

• A range of operational and measurement conditions were sampled for each vessel category (Fig. 3). The spread in STW values and the bi-modal distributions of some categories were due to the inclusion of measurements from the Haro Strait and Boundary Pass slowdown periods in the data set (see MacGillivray et al., 2019; Burnham et al., 2021).

• Scatterplots and density plots indicated that logarithmic transformations were appropriate for linearizing observed trends in most predictor variables. Exceptions were for variables containing negative or zero values (e.g., wind resistance) and categorical variables (e.g., engine type). The surface grazing angle was not subjected to a logarithmic transformation as it spanned a small range of values.

• Several design characteristics were dominated by a single value, which made them unsuitable for statistical analysis. For example, the overwhelming majority of vessels used conventional diesel propulsion (only seven vessels used diesel-electric or gas turbine propulsion), fixed-pitch propellers (97.3% of vessels), and had a single two-stroke engine (96.7% of vessels). The homogeneity of these variables meant that meaningful correlations could not be inferred from the data.

• In general, design characteristics of tugs and cruise vessels were considerably different from those of the other categories.

Univariate correlation coefficients were used to identify which predictors had the strongest correlation with vessel source levels over a wide range of frequencies. In general, the strongest univariate correlations were observed for the two main operational parameters: STW and actual draft. Other parameters were observed to have weaker correlations with vessel source levels. Notably, the derived variables described by Eqs. (2)–(4) (block coefficient, speed fraction, and draft fraction) had lower univariate correlation coefficients than their component variables. This suggested that including these derived variables in a multiple-predictor model would not provide any more explanatory power than including their component variables. Furthermore, logarithmic transformation of these derived variables would simply transform them to sums in the functional regression model [Eq. (5)]. Thus, the derived variables were excluded from consideration in the final functional regression model.

Four design parameters, all related to vessel dimensions, were found to be strongly correlated with each other: length overall (LOA.LLI), gross tonnage (Gross.LLI), moulded breadth (Breadth.Moulded.LLI), and displacement (Displacement.LLI; Table III). Due to their strong collinearity, these variables could not be added together during stepwise development of the final functional regression model. To address this issue, length overall (LOA.LLI) was selected as the primary size-related design parameter during initial model development. The remaining size-related variables were held back and added only after all of the other design variables were added. This was done as a final check to assess whether they provided additional improvement to the final model.

While the regression results were generally different for each vessel category, incremental analysis of the stepwise R² values showed that the operational variables (STW, actual draft) explained the largest fraction of source level variability. The design parameters, by contrast, generally explained much less of the variability. In particular, the additional size-related variables did not generally improve the fit of the model to the data, therefore, they were excluded. Thus, the final model included nine predictors and had the following form:

where $yif$ was the estimated MSL (⁠$LS)$ for measurement i, and the predictors included in the model were as follows (with the database variable name indicated in parentheses):

1. Actual draft, $d$ (Actual.Vessel.Draft);

2. STW, $v$ (STW);

3. wind resistance, $Kw$ (Wind.Resistance);

4. surface grazing angle, $φ$ (Surface.Angle);

5. length overall, $L$ (LOA.LLI);

6. nominal main engine revolutions per minute (RPM), $n$ (Main.Engine.RPM.LLI);

7. total main engine power, $P$ (Main.Engine.kW.LLI);

8. design speed, $v0$ (Speed.LLI); and

9. vessel age, $a$ (Vessel.Age).

The model of Eq. (6) was fit to the ECHO dataset separately for each of the six vessel categories by minimizing the model errors, $εi(f)$⁠, as a function of frequency.

The ranges of the various design parameters used to develop the final model were summarized in box-and-whisker plots for each category (Fig. 4). Functional regression analysis was applied to the full set of MSL measurements using the nine predictors listed above. This resulted in sets of $βf$ functions that described the frequency-dependent slope of the trend between (transformed) predictors and decidecade band source levels in each vessel category (Fig. 5). The regression functions indicated how strongly the source levels at each frequency correlated with each predictor (either positive or negative). Confidence intervals (95%) were calculated for the regression functions (Fig. 6). The correlation was not considered statistically significant where the confidence intervals crossed zero. Note that functional regression analysis was also applied to the RNL data.¹ This showed that within the same category regression functions for a given predictor were generally very similar for the RNL and MSL with the notable exception of actual draft at frequencies below approximately 100 Hz. The reason for this latter difference is that the draft is directly involved in the MSL calculation, whereas this is not the case for the RNL (see Sec. IV A).

Note that the unusual trends observed for cruise vessels for some predictors were likely caused by the small number of measurements included in this category (only 136 were included after subsampling for repeat measurements). Unlike other categories, cruise vessels were extremely heterogeneous, hence, large deviations from the mean could be due to a single vessel. More data would likely be needed to better understand the relationships between predictors and source levels for this category. Tugs also exhibited clearly distinct trends from the other categories. This may reflect the fact that tugs had very different design characteristics from the other types of vessels included in the analysis (see Fig. 4) and there was a much greater fraction of missing design data for this category (see Table II).

The unexplained variability in the functional regression model was quantified by calculating the residuals of the statistical model in different frequency ranges (Fig. 7). Recall that the residuals are the differences between the estimated and measured source levels [cf. Eq. (5)] such that

The root mean square error (RMSE) of the residuals provides an approximate estimate of the prediction error in decibels of the functional regression model versus frequency (see Sec. IV C) such that

where N is the total number of measurements. This analysis showed that the model residuals were approximately normally distributed with RMSE ranging from 3.3 to 8.0 dB, depending on frequency range, and with a mean RMSE of 5.1 dB when averaged across the vessel category and frequency range (Table IV). Approximately normal distributions of the residuals for a population of ships with RMSE in the range of 3–7 dB was similarly observed in earlier measurements by Scrimger and Heitmeyer (1991) and Wales and Heitmeyer (2002). Depending on the frequency band and vessel category, the functional regression model was generally able to explain 25%–50% of the variance in the observed source level measurements in the ECHO data set (Fig. 8). This was reflected by the fact that the model was able to accurately reproduce the broadscale features of source level profiles (i.e., the smooth MSL versus frequency trend) but not the fine-scale features of the profiles.¹ The fine-scale features—particularly the narrowband spikes—were not reproduced by the model because they did not follow a predictable trend between different vessels. Furthermore, the functional regression model intentionally smooths the regression coefficients across frequencies to constrain the beta coefficients to be similar at neighboring frequencies and minimize random error in the predictions. This was also consistent with the observation of Simard et al. (2016) that “[source level] variability appears to be an intrinsic property of a diverse fleet.”

The regression coefficient functions for the six vessel categories were used to create a spreadsheet implementation of the functional regression model that can be used to predict the source level versus frequency curve for a vessel, given its set of operational and design parameters.¹ It should be noted that the surface grazing angle does not represent the directivity of the source level with respect to the sea-surface. As discussed in Sec. II C, this parameter was intended to control (and correct) for any systematic errors in the estimated source level due to differences in the measurement geometry. Thus, when using the functional regression model in a predictive sense, the intent is to use a nominal fixed value of the surface grazing angle to represent the mean measurement condition (i.e., 30 deg, corresponding to the preferred ISO 17208–1 measurement geometry).

To aid in the interpretation of the regression functions, predicted source level plots (Fig. 9) were created to show the influence of individual predictors on source levels for an average vessel in each category. That is, for each vessel type, the median of each covariate (vessel characteristic) value was taken to represent a “typical vessel.” These plots showed the predicted source levels obtained by changing the value of a single predictor (color-coded) while holding all other predictors at fixed average values. For predictors having more than 200 possible values in the data, 200 values were randomly selected, along with the minimum and maximum value. These predicted source level plots were used to gauge the influence of the different predictors on vessel source levels.

The statistical analysis showed that STW and actual vessel draft (i.e., the two primary operational parameters) were generally the most influential predictors of source level for all of the vessel categories. Wind resistance had only a marginal influence on source level, which appeared to be statistically significant only at low frequencies (below approximately 100 Hz). Higher STW was associated with higher underwater radiated noise for all of the vessel groups. It had the greatest influence at high frequencies (>1000 Hz), where cavitation dominates, and the smallest influence at intermediate frequencies (100–1000 Hz), where machinery noise generally dominates at speeds below cavitation inception. This observation was consistent with results obtained during the 2017 ECHO slowdown study in Haro Strait (MacGillivray et al., 2019). This trend is to be expected for conventional fixed-pitch propellers, which is the predominant type of propulsion employed by vessels in this this data set, but it is not necessarily the case for controllable-pitch propellers (Baudin and Mumm, 2015; Traverso et al., 2015; McIntyre et al., 2021).

Actual vessel draft also had a strong influence on vessel source levels, primarily at frequencies above 100 Hz. As with STW, the strongest influence of actual vessel draft was at high frequencies (>1000 Hz), where cavitation dominates the noise spectrum. Although actual draft was correlated with reduced MSL below 100 Hz, this correlation was attributed primarily to the linear dependence between monopole intensity and source depth (which is proportional to draft) at these frequencies. This is because the radiated sound field is dipole rather than monopole in nature when the sound source is within a quarter wavelength of the sea-surface, as is typically the case for surface vessels at low frequency (Ross, 1976). Inspection of the corresponding RNL regression coefficients, which are more closely related to the dipole intensity at the source, suggested that actual vessel draft had a negligible or slightly negative influence on radiated noise.¹ Increasing draft likely increases underwater radiated noise because it increases hydrodynamic drag and puts more of the hull area in contact with water.

Each of the six vessel categories shared a distinct set of design characteristics (see Fig. 4), and the functional regression analysis likewise showed that each category exhibited distinct source level versus frequency patterns and trends. Whereas the design characteristics were generally less influential on vessel source levels than the operational parameters, it was nonetheless possible to use the functional regression model to rank the influence of the different design characteristics for each vessel category. Rankings were not provided for the cruise category because the statistical significance of the results was questionable given the relatively small sample size. Similarly, rankings were not provided for vehicle carriers because their design characteristics did not exhibit a sufficient range of variation to demonstrate clear trends.¹ For the remaining four vessel categories, the functional regression model suggested the following trends:

• For bulkers and tankers, length overall and main engine RPM were the two most influential design characteristics, followed by main engine power and vessel age. Increased engine RPM was associated with higher source levels at all frequencies. The influence of length overall was frequency dependent: greater length was associated with higher source levels at low frequency (<100 Hz) and lower source levels at high frequency (>100 Hz). Older vessels in these categories tended to have higher source levels than newer vessels. Main engine power and design speed appeared to have a secondary influence in these categories.

• For container ships, length overall, main engine power, and design speed were the most influential design characteristics. Greater length overall was associated with higher source levels at all frequencies, whereas greater engine power was associated with lower source levels at all frequencies. Main engine RPM and vessel age did not appear to be very influential for this category.

• For tugs, main engine RPM was the only design parameter that appeared to have a clearly significant trend with source levels. Other design characteristics (length overall included) had weak correlations with tug source levels and their regression coefficient functions were not clearly significant over a wide range of frequencies. Thus, design characteristics for tugs were difficult to associate with underwater noise emissions despite the large number of measurements in the data set.

• In all of the categories, higher design speed was associated with lower source levels below 100 Hz and above 1000 Hz (i.e., at frequencies where cavitation dominates). A possible explanation for this trend is that higher design speed was associated with increased propeller cavitation inception speed (see Sec. IV E).

When interpreting these results, it is important to note that the statistical methods employed by this study only had the ability to examine correlation not causation. The analysis was also limited by the sampling methods inherent to the dataset, which were collected for vessels of opportunity transiting in and out of the Port of Vancouver (i.e., not in a manner that controlled for design parameters). For example, the functional regression with multiple predictors showed that larger main engine power was associated with lower underwater radiated noise for the container and vehicle carrier categories. This result is difficult to interpret unless one also considers that main engine power was strongly correlated with length overall in this category, and greater length is more strongly associated with higher radiated noise. The underlying relations of length, engine power, and radiated noise are very difficult for a statistical analysis to pull apart because the underlying predictors are often strongly correlated.

The relationship of design predictors with source levels was generally weakest for the tug category. This could be explained, in part, by the fact that tugs had more missing design information in the database than other categories (see Table II). It could also be related to the fact that radiated tug noise depends on the operating mode of the vessel when transiting. For example, a tug performing escort duties (not towing or pushing) would likely have different noise emissions than a tug engaged in towing or part of an articulated-tug-barge unit. Information was not available in the ECHO database on whether tugs were involved in towing or pushing, which could affect their radiated noise emissions.

Like any model, predictions of the functional regression model developed in this study have an associated error or uncertainty. If this model is intended to be used in a predictive fashion, it is important to quantify the uncertainty of its predictions. The functional regression model itself captures the smoothed trend of the measurements, but there remains a component of the data variability that is unexplained by the model. In statistical terms, this residual variance can be used to estimate the uncertainty of the model predictions by treating it as a random variable and calculating a prediction interval. Functional regression is a relatively new statistical technique and, therefore, methods for calculating prediction intervals for functional regression models are not widely available in popular statistics packages. Nonetheless, methods for calculating prediction intervals for simpler linear regression models are well understood and can be used as a guide for estimating approximate uncertainties. In linear regression, the standard error of the prediction is proportional to the Pythagorean sum of the RMSE and standard error of the predicted value (Devore, 1995). In the present case, the RMSE term generally dominates the sum when the desired prediction is within the range of the covariate data that was used to develop the model.¹ Thus, we expect that the RMSE (see Table IV) gives a good approximation to the model uncertainty when the desired prediction is not extrapolated too far outside the range of the vessel designs comprising the ECHO dataset (see Figs. 3 and 4). Uncertainties in vessel noise predictions remain an important topic for future research.

Early studies based on measurements of World War II (WWII)-era vessels recognized that vessel speed and size were strong predictors of underwater radiated noise (Ross, 1976). A straightforward empirical method for estimating source levels, which emerged from this early research, was to adjust a reference spectrum according to an assumed power law dependence on speed and length (Breeding et al., 1996). This approach has been refined by several investigators using more recent source level data sets to incorporate class-specific spectrum levels and different speed and length trends (Chion et al., 2017; Jiang et al., 2020; MacGillivray and de Jong, 2021). Nonetheless, such simple models do not address the question as to what drives differences between similar vessels.

Three recent studies applied statistical methods to investigate a wider set of design and operational parameters using radiated noise datasets for vessels of opportunity. Simard et al. (2016) applied various statistical techniques to analyze MSL measurements of 255 commercial vessels measured in the St. Lawrence seaway, collected according to the ANSI S12.64 (ANSI, 2014) methodology. Their analysis found statistically significant trends with ship length, breadth, draft, and speed, but did not find a significant correlation with wind speed. They presented three different possible source level models with varying numbers of predictors and polynomial trends that furthermore varied with frequency in the range 20–12 000 Hz. However, the forms of these statistical models were rather different than those employed for the current study, and the influence of the various covariates was difficult to directly compare to the trends obtained using the functional regression model.

McKenna et al. (2013) used a general additive model (GAM) to investigate trends of radiated noise with operational conditions, design characteristics, and oceanographic data for 593 radiated noise level measurements from the Santa Barbara Channel. They reported significant, frequency-dependent trends with vessel speed and length in the range of 20–1000 Hz, as well as significant trends with wind and wave conditions. They did not observe a significant influence of draft, although the authors suggested that this could have been due to the lack of variation in the loading conditions of the vessels, which were measured purely in the outbound direction from the Port of Los Angeles (the ECHO database measured inbound and outbound vessels from the Port of Vancouver under a wide range of loading conditions; see Fig. 3). A prior study by McKenna et al. (2012) also identified vessel type as an important predictor of radiated noise, which was a key underwater radiated noise predictor in the functional regression model. Once again, however, the GAM model employed by McKenna et al. (2013) was quite different than the functional regression model from the current study; therefore, it was difficult to directly compare the observed trends.

Finally, a recent study by Chion et al. (2019) applied a generalized linear mixed model (GLMM) to analyze publicly available vessel source level datasets, including a set of over 2100 radiated noise measurements from Veirs et al. (2016). The GLMM analysis identified significant trends with speed, breadth, and ship class. However, the reported relationships were obtained only for broadband radiated noise and were, thus, not easily compared to the frequency-dependent source level trends obtained using the functional regression model. These comparisons suggest that one possibility for future research would be to conduct a parametric analysis using the statistical models of McKenna et al. (2013), Simard et al. (2016), and Chion et al. (2019) to evaluate the consistency of the observed trends from different datasets reported by different investigators. While absolute source levels are likely to be different between studies, due to a lack of standardization in measurement methods (Chion et al., 2019; MacGillivray and de Jong, 2021), it may nonetheless be possible to compare trends in underwater radiated noise with changes in predictor values between studies.

Reducing underwater noise from shipping requires an understanding of which factors are driving underwater radiated noise from marine vessels. When considering ship design from a noise control perspective, the most important factors to consider include the following (Bahtiarian, 2017; Spence and Fischer, 2017; Wittekind and Schuster, 2017):

• Factors influencing propeller noise and cavitation inception, including propeller design characteristics, propeller placement, and uniformity of wake outflow; and

• factors influencing transmission of machinery noise through the hull, including use of quiet machinery, vibration isolation of machinery, and acoustic insulation of machinery compartments.

No parameters describing these factors were part of the LLI database nor were they available from other publicly available sources. It is, furthermore, unlikely that such information could easily be obtained for thousands of vessels of opportunity. However, machinery sound and vibration are influenced by main engine power and ship speed. The exposed hull area is related to length and draft. Propellers are designed for maximum efficiency at the vessel service speed, but this generally comes at the cost of reduced cavitation inception speed (Gray and Greeley, 1980). Thus, the trends identified by the functional regression analysis from this study may point to which design factors are most important for controlling radiated noise for specific frequency bands or vessel classes. Future studies could investigate the relationships between these underlying design factors and functional regression model predictions from this study. The results of this study could also be used to refine alternative vessel noise modelling approaches, which are based on the detailed engineering parameters considered by naval architects, such as the methods of Wittekind (2014) or Audoly and Rizzuto (2015).

One limitation of the approach taken in this study was that some categorical information related to engine design that was present in the LLI database (e.g., engine model and make) did not easily fit into the framework of functional regression analysis. This is because categorical variables cannot be ranked or represented as a continuously varying quantity. Such information might be amenable to future analysis if additional data related to the engine design (and noise and vibration) were to be obtained, e.g., from the equipment vendors. Alternatively, categorical data on engine design might be investigated using descriptive statistics or by focusing on a smaller subset of the data.

The largest data gap identified during this study was the lack of publicly available information related to propeller design and cavitation inception speed for commercial vessels. Propeller cavitation is the dominant source of underwater radiated noise for most surface ships. Therefore, efforts to control radiated noise must address propeller cavitation. However, data on propeller design characteristics are not available in public databases and must often be obtained directly from shipbuilders (who may consider such information proprietary). Alternatively, it may be possible to estimate cavitation inception speed based on other design parameters as in Jalkanen et al. (2018). Future work could focus on obtaining data on propeller design and cavitation inception data for a subset of vessels in the ECHO dataset to determine their influence on the radiated noise.

Functional regression analysis of a large database of 9880 unique vessel noise measurements found that two operational parameters (STW and actual draft) were the most influential predictors of source levels for six categories of vessels. Within each category, five design characteristics were generally less influential on vessel source levels than the operational parameters. Of the design characteristics, vessel size (represented via length overall) was the design parameter with the strongest correlation to underwater radiated noise for three categories of vessels. Other design parameters that were investigated (main engine RPM, main engine power, design speed, and vessel age) had weaker but nonetheless statistically significant correlations with underwater radiated noise. In all of the categories, higher design speed was correlated with reduced source levels at frequencies dominated by cavitation. Wind resistance had a minor influence on source levels, which was nonetheless significant below 100 Hz. The explanatory power of the model (i.e., the coefficient of determination, R²) was largest for those vessel categories in which vessels were measured while operating under a wide range of speed and draft conditions. Depending on the frequency band and category, the functional regression model was generally able to explain 25%–50% of the variance in the observed source level measurements in the ECHO data set. The standard deviation of the residual model errors was 5.1 dB when averaged over vessel category and frequency band, which reflected the approximate uncertainty of the model predictions.

Some of the influential factors identified in this study, such as vessel size, speed, and class, were also identified by earlier statistical studies (McKenna et al., 2013; Simard et al., 2016; Chion et al., 2019), however differences in data collection and analysis methods made it difficult to directly compare the reported trends to earlier studies. Thus, one possible avenue for future investigation would be to validate the functional regression model from this work against other source level datasets. Future work could also focus on associating the influential parameters identified during this study with more detailed ship design characteristics used by naval architects. The largest data gap identified by the present study was the lack of publicly available information regarding propeller design characteristics, which are undoubtedly of importance in determining cavitation noise, a dominant source of vessel underwater radiated noise. The ECHO program is continuing its data collection, research, and outreach efforts to address these and other important topics aimed at reducing underwater radiated noise from vessels in the marine environment.

The authors would like to acknowledge the Vancouver Fraser Port Authority (Vancouver, BC) and Transport Canada (Ottawa, ON) for funding the deployment of underwater listening stations for the collection and analysis of vessel source level measurements. The ECHO program's ATC provided valuable advice and guidance to the design of the study, and Michael Bahtiarian (Acentech, Cambridge, MA) provided additional expertise on the selection of potentially influential design parameters to be included in the statistical analysis. Many technical and support staff at JASCO Applied Sciences Ltd. (Victoria, BC) contributed to the data collection and analyses projects that made this study possible.

The cross section of a vessel's hull that sits above the water line experiences an aerodynamic drag force that depends on the speed of airflow around the hull (i.e., apparent wind speed) and the cross-sectional hull area that is exposed to the airflow. In naval architecture literature, the propulsive power required to overcome this drag force is taken to be proportional to a parameter known as the heading coefficient (C_γ; Lewis, 1988), which depends on the relative wind heading with respect to the direction of travel of a vessel. It is positive (representing a resistive force) when the apparent wind direction is toward the bow of the ship and negative (representing a driving force) when the apparent wind direction is toward the stern of the ship (Fig. 10). For calculating wind resistance, we use the following piecewise polynomial approximation to the heading coefficient from Fig. 33 of Lewis (1988):

where $θ$ is the direction of the apparent wind vector (⁠$0°≤θ≤180°$⁠) with respect to the heading of the vessel such that $θ=0°$ is a pure headwind and $θ=180°$ is a pure tailwind. The apparent windspeed vector is calculated from the vector difference

where $vw$ is the windspeed vector with respect to ground, $vsog$ is the speed over ground vector of the vessel, and $vw′$ is the apparent windspeed vector in the reference frame of the vessel. The direction of the apparent windspeed vector is equal to

where $ĥ$ is a unit vector representing the vessel heading.

A dimensionless wind resistance factor (K_w), assumed to be proportional to the power required to overcome the aerodynamic drag force, was then calculated from the product of the square of the apparent wind speed and heading coefficient such that

This wind resistance factor was used as a predictor in the statistical analysis rather than separate wind speed and direction as it better reflected the increase in propulsion power (and, thus, associated noise and vibration) required to overcome wind-induced drag forces.