The current study focuses on a seafloor area (37 km²) adjacent to the Doce River mouth, situated on the eastern Brazilian continental shelf, in the state of Espírito Santo. The Doce River is one of the main river deltas recognized along the Brazilian eastern coast (Dominguez, 2006), and is a major sediment source to the adjacent continental shelf (Bastos et al., 2015). In terms of sedimentary processes, Quaresma et al. (2015) characterized a deltaic lobe extending to 30 m depth, with an accumulation of fine sediments (>75% mud) in a main depocenter south from the river mouth. Towards offshore (depths >30 m), sandy facies are predominant (Figure 1)—a trend also observed and corroborated by Vieira et al. (2019). Though useful, previous research in this area has been conducted at a coarse spatial resolution, and fine scale geomorphological and substrate characterization are still lacking. The MBES and sediment sample data sets presented here were analyzed to fill this knowledge gap.

A multispectral MBES dataset was acquired to provide high resolution information on the geomorphology and substrate distribution of the study area. Bathymetry and backscatter were acquired using an R2Sonic 2024 echosounder, with the sonar head deployed through a moon pool in the steel-hulled survey vessel, “Santa Edwiges” (19.6 m in length). The MBES was configured to collect data by sequentially pinging at 170, 280, 400, and 700 kHz (true frequency of 697.674 kHz) operating frequencies in equiangular mode. The 700 kHz frequency did not operate during the entire survey due to the effects of depth and water turbidity on data quality. Specifications such as power, pulse length, gain, and spreading were held constant. A 90° angular sector and 256 beams were used for all frequencies except 700 kHz, which is limited to a 70° sector. Because the receiver array is flat, the size of the receive beam grows with distance from nadir (i.e., the lower the grazing angle, the larger the beam footprint; Figure 1). The transmit array is round and does not suffer the same issue.

The MBES system was paired with a POS MV Wave Master Inertial Navigation System (INS), with differential positioning. During the survey, SVP casts were deployed every 3 hours using a Valeport Mini configured to collect sound velocity, salinity, temperature, and pressure. These data are essential to ensure the correction of sound velocity effects and to calculate and apply the absorption coefficient. Application of absorption coefficients for each frequency was the only radiometric correction applied during the data acquisition—all other corrections (e.g., beam width) were applied during post processing. All systems were integrated using QPS QINSy 8.18.3 for data acquisition, and the survey lines were planned with ∼30% overlap. The dataset was assessed during the survey to ensure data quality using QPS Qimera 2.0, and post-processing was carried out using QPS Qimera and QPS FMGT 7.9.5 (Fledermaus Geocoder Toolbox, hereafter FMGT).

In compliance with the International Hydrographic Organization, bathymetric processing consisted of manually cleaning erroneous soundings and tidal correction, which was performed using data extracted from an ADCP moored adjacent to the study area. The dataset was filtered for each operating frequency in QPS Qimera. Digital bathymetric models (DBM) for each frequency were exported at a 2 m horizontal resolution.

Although there is still no unified standard for processing backscatter, the multispectral data were processed following recommendations in Lamarche and Lurton (2018) and Schimel et al. (2018). The main steps involve the frequency filtering, correction of acquisition parameters for each frequency (gain, transmit power, pulse length, beam width) and AVG correction. Backscatter mosaics were exported at 1 m resolution. All backscatter processing was completed using FMGT.

Elements of both image-based and angular range approaches were used to analyze backscatter data and explore the multispectral response of the seabed. The methodological framework of the analysis is summarized in Figure 2, outlining each approach and its expected results.

The large amount of information from the different analysis approaches (Multispectral Backscatter Analysis) was summarized spatially using density-based clustering. Unlike many clustering methods that rely on measures of multivariate distance to identify clusters, these approaches identify clusters using estimates of multivariate data density (Hahsler et al., 2019; Hashler and Piekenbrock, 2021). Density-based clustering offers several advantages compared to many other approaches, such as the automatic identification of cluster numbers, rejection of outliers that do not fall within any dense regions, and clustering of data of arbitrary multivariate shape (Ester et al., 1996; Kriegel et al., 2011). Ordering Points to Identify Clustering Structure (OPTICS) is a method that can additionally identify clusters of varying density (Ankerst et al., 1999). OPTICS is an extension of the DBSCAN algorithm (Ester et al., 1996), which facilitates data exploration through “reachability plots”—one-dimensional representations of the “reachability distance” between points (Hahsler et al., 2019). These can be used to inform the clustering outcome, and enable additional flexible solutions such as hierarchical density-based clustering. Despite the widespread application of unsupervised approaches for seabed mapping and classification, and the potential benefits offered by density-based clustering, it has received relatively little attention in this field (but see Le Quilleuc et al., 2021).

Results from both image-based and angular range approaches were used for clustering, with the exception of the frequency of 700 kHz, which has a smaller coverage area and was not considered for this analysis. These layers included AVG backscatter mosaics for each frequency (170, 280, and 400 kHz), GLCM variance for each frequency (170, 280, and 400 kHz), and backscatter values from four different angular sectors (1–10°, 11–20°, 21–30°, 31–40°) for each frequency (170, 280 and 400 kHz), totaling 18 rasters (3 AVG mosaics, 3 rasters from GLCM results, and 12 rasters from ARA). We applied the methods presented in Che Hasan et al. (2012) in order to aggregate the acoustic data to object-based image segments. The three-band RGB backscatter mosaic was segmented using the “mean shift” algorithm (Yizong, 1995; Comaniciu and Meer, 2002) in ArcGIS Pro using the maximum spatial and spectral detail, but a minimum segment size of 2000 pixels. The mean shift is a “mode-seeking” algorithm that assigns data points to clusters based on their position relative to the local maxima of a multivariate kernel density estimate. The values of the 18 rasters were stacked and averaged over the resulting segments, which allowed for continuous estimates of the angular backscatter values across the full extent of the study area, and a reduction of data volume by several orders of magnitude, producing 6902 segments (step 2 of Figure 3). The values for each segment (18 acoustic variables) were converted to matrices in R, and principal component analysis (PCA) was applied to eliminate collinearity and reduce dimensionally. The first two components explained 95.4% of the total variance and were retained for clustering.

Principal components of the image-based segments were clustered using OPTICS to identify acoustically distinct seabed types across the study area. OPTICS was implemented using the DBSCAN package in R (Figure 3; Ester et al., 1996, Campello et al., 2013; Hahsler et al., 2019). The reachability distance was first plotted using the parameters (the density search radius) and (a threshold number of points used for reachability calculations) to visualize the data density and estimate the number of clusters. The reachability plot suggested a hierarchical data structure, and two sets of clusters were extracted at coarse and fine levels of detail. These indicated broad data groupings, which contained smaller sub-clusters at a higher level of detail. The clustering results were exported to GIS format to generate benthic map predictions.

After preliminary analysis of the backscatter mosaics, samples of surficial sediments were collected at 12 sites within the survey area (Figure 1). Sediments were collected using a Van Veen grab sampler in July 2021 (nearly 2 years after the acoustic data survey). To account for the temporal offset between the collection of sediment samples and acoustic data, we verified the differences and the general trend from the sediment distribution (Quaresma et al., 2015; Vieira et al., 2019). The grain size analysis provides results as percentage gravel, coarse sand, medium sand, fine sand, silt and clay, which were used to ground-truth results from the multiple backscatter analysis approaches and density-based clustering.

The DBM outputs, which contain all processed depths from all frequencies combined, ranged from 10 to 32 m, and revealed a predominantly flat bottom with locally steep features. These were evident in maps of local slope and roughness, and ridges oblique to the coast were visible at the deep end of the profile (Figure 4).

The spatial distribution of depth differences was also investigated. DBM surfaces generated for each operating frequency were subtracted from each other using raster algebra to visualize the distribution of depth differences across the study area (Figure 5). In general, regions close to higher slope values (steep feature and the edges of oblique ridges) show higher bathymetric difference between frequencies. Differences between 170 and 280 kHz were, on average, 0.1 m greater than the differences between 280 and 400 kHz, and differences between 170 and 400 kHz were greater still. The differences between 700 and 170 kHz were the greatest observed, whereas differences between 400 and 700 kHz were three times lower. The bathymetric differences (which can reach up to 0.2 m) observed between the frequencies largely conformed to expectation given the current understanding of substrate-frequency response and the bottom type in the study area.

The angular range analysis provides results that can be explored using geophysical and empirical approaches, including physical model comparison and empirical parameters to distinguish seafloor types, and also angular dependence curves of different substrates (Lurton and Lamarche, 2015). The angular backscatter values were interpolated spatially using Inverse Distance Weighting (IDW) in ArcGIS Pro to produce grids of 1 × 1 m cells (Figure 11). The results demonstrated that, for most of the area, decreasing backscatter was associated with increased incidence angle and frequency. The exception occurs over the high-backscatter feature, which marks a region with high scattering across the entire swath.

The 1–10° sector shows high backscatter values throughout the area and similar results for all frequencies. Compared to results presented for the lower frequencies at 11–20° sector, the 400 kHz frequency seems to achieve a better power of seafloor discrimination for the shallower region. Three distinct regions are clearly visible using the 21–30° sector. Additionally, the 400 kHz data within this sector may suggest that the shallow region is composed of different material than the elongated feature, even though they appear similar at other frequencies and angular sectors. Lower backscatter was observed using the outermost angular sector over much of the study area—primarily at 400 kHz—but high values were still observed at the elongated feature.

Angular response curves were extracted for all frequencies to characterize the acoustic response of the seabed at each sample site (Figure 12). The angular response curves can be broadly described in three main groups. The curves at sample sites Doce 01 and Doce 02 show similar shapes for all frequencies, with decreasing backscatter values as the frequency increases (with exception of 700 kHz). The shape shows low backscatter level loss for the first 20° from the nadir, followed by a greater decrease in backscatter with increasing incidence angle. The curves extracted at Doce 07, 09, 11, and 12 are comparatively flat, with high backscatter values and very slight differences between frequencies. The third group is characterized by a high backscatter level loss for all frequencies from nadir to the outer beams, and lower values at 400 kHz.

The combination of PCA and density-based clustering enabled synthesis of a large number of backscatter data layers to inform seabed classification. PCA reduced the 18 backscatter-derived layers to two principal components that explained the majority of variance (>95%) in the dataset, suggesting substantial collinearity among the variables. The OPTICS reachability plot from these principal components (Figure 13) enabled identification of clusters at two hierarchical levels. Low reachability distances, represented as valleys in the plot, indicate dense data regions, which are separated by “peaks” representing regions of low data density. First, two broad, well-defined clusters were evident at a reachability distance of 0.88 (Figure 13A). These results agree with general trends observed for many of the backscatter layers (Figures 6, 10, 11), in which the study area is separable into distinct areas of low and high backscatter intensity and variance. Cluster 1, which occurred across the full extent of the study area (Figure 14A), could be further divided at a reachability distance of 0.48 to produce four sub-clusters. Unlike the broader Cluster 1, the four sub-clusters show a distinct gradient from the coast (Figure 14B), which was not immediately apparent from visual analysis of the various backscatter layers. Outliers appear in regions of extreme high and low backscatter values.

The results from grain size analysis (Figures 14C,D) show that samples Doce 01 and Doce 02 are composed primarily of fine sand. Doce 03 and Doce 04 share a similar mud content and both stations are close to bathymetric anomalies (high slope and high ruggedness sites). Samples Doce 05, Doce 08 and Doce 10 (Figure 14F) had an increased mud content, while sample sites at the high backscatter features (Doce 07, 09, 11, 12) - e.g., Figure 14E - had higher proportions of medium and coarse sand fractions, with low amounts of mud.

Although there not enough samples to estimate the classification accuracy, the classes seem to distinguish benthic clusters through the multivariate combination of the results. Table 1 summarizes the main properties for each cluster defined.

Cluster 1 contained the greatest number of sediment samples and can be characterized as muddy fine sand. Cluster 2 is medium sand, with a smaller coarse sand fraction. Cluster 1a contains no ground truthing validation, but seems to indicate a transition between clusters 1 and 2. Cluster 1b seems to have characteristics close to cluster 1a, likely with higher mud content, or morphological characteristics that provide a soft backscatter response. Cluster 1c also does not have sediment samples, but it can be interpreted as similar to cluster 1b, according to the angular sector maps (Figure 11). Cluster 1d has two samples and can be defined as fine sand. Some areas were left unclassified but may comprise two additional substrate types—one with higher coarse sand/gravel contents (sample Doce 07) and the other with extremely low backscatter values, indicating a muddy bottom. It is possible that the latter has undergone modifications during the time between the collection of acoustic data and ground truthing.

Highly dimensional input data resulting from multiple analysis approaches can be difficult to synthesize spatially, especially given sparse ground truth information, yet density-based clustering provides an objective means to identify well-defined seabed units under such circumstances. Here, backscatter mosaics (Figure 6) and the composite RGB (Figure 9) potentially suggest two or three seabed classes, while the angular response curves at the 12 sample sites (Figure 12) suggest at least three classes. GLCM products and spatial analysis of different angular sectors further confound a unified manual seabed segmentation solution. The combination of image-based segmentation, dimensionality reduction, and density-based clustering provide an objective solution to synthesizing these data spatially to identify distinct seabed types. Similar solutions have applied PCA for dimensionality reduction of backscatter data followed by k-means clustering for unsupervised classification (e.g., Alevizos et al., 2018), and other approaches towards similar ends have included combinations of self-organizing maps and fuzzy clustering (Chakraborty et al., 2015), and Bayesian probability estimation (Amiri-Simkooei et al., 2009; Simons and Snellen, 2009). Several properties of density-based clustering are desirable in an exploratory context though, including automatic identification of number of clusters, enhanced data visualization and hierarchical clustering solutions, and rejection of outliers that are dissimilar from other well-defined clusters. These qualities are lacking in many other common approaches, which may be well-suited to classifying gradational or poorly defined seabed types (e.g., k-means, ISO, agglomerative clustering).

Quantifying the advantages of using multispectral backscatter data as a predictor for different seafloor types is beyond the scope of the current study, but the findings suggested that the multifrequency acoustic data provided greater discrimination of muddy and fine sand sediments than coarser sediments, aligned with findings from other studies (Gaida et al., 2018; Brown et al., 2019; Costa, 2019). Part of the data identified as noise by the OPTICS algorithm is located in the region of highest backscatter, where the sediment shows the higher gravel and coarse sand content. Conversely, multiple distinct clusters were identifiable in muddy sediments. We believe this to be a property of the acoustic data, which were highly variable and inhomogeneous at the high backscatter regions. Further investigations using multispectral backscatter for a greater range of seabed types is necessary to determine the potential of this technology for substrate discrimination. This is also the first application of density-based clustering methods in this context as far as we know, and further comparison with other clustering methods in various environments would be of great interest.

Despite limitations associated with the use of uncalibrated backscatter data (Lurton and Lamarche, 2015; Schimel et al., 2018; Malik et al., 2019), the results provided basic information for interpretation of different seafloor types, and provided useful information on the acoustic response of the seabed. For most areas, the median uncalibrated backscatter values from the mosaics for all frequencies were low (slightly higher for lower frequency), which is characteristic of soft muddy deposits (Brown et al., 2019; Diesing et al., 2020). It is also worth noting that the characteristics of sub-bottom deposits can influence the acoustic response, especially for lower frequencies (Jackson and Richardson, 2007; Williams et al., 2009; Feldens et al., 2018). The acoustic response from sandy features (both for image-based and ARA results) is largely influenced by interface scattering, but could be further investigated by observing the influence of survey azimuth (Lurton et al., 2018), presence of bedforms, the variety of microtopography, and even more ground truthing to check lateral variability.

The use of multispectral technology offered increased seafloor discriminatory power based on the different frequency response of the seafloor across the frequency range, achieving a benthic prediction in agreement with earlier studies in this region (Quaresma et al., 2015; Vieira et al., 2019), and making it possible to improve the seafloor classification. Quaresma et al. (2020) found lower seabed density values in the shallow region close to the mud depocenter (Figure 1), and high density values (>1400 kg/m³) in regions with depths greater than 30 m. Here, the lower uncalibrated backscatter values appear close to the steep feature (i.e., high slope–see Figure 4) and on the edge of the sandy elongated feature, where the bathymetric differences (Figure 5) between the frequencies reached up to 10 cm, indicating a mud accumulation partially driven by the seafloor morphology. In these regions, the lower frequency was presumably detecting the sub-surface sediment characteristics, while the higher frequency reflected primarily off the surface sediments of the seafloor, potentially indicating a thick, muddy deposit in this area. Although results from the bathymetric data (e.g., slope, ruggedness, bathymetric differences) played a secondary role compared to the backscatter here, they can serve as important predictors for seabed classification (Walbridge et al., 2018; Lucatelli et al., 2019), and some derivatives may be more important than bathymetry itself (Trzcinska et al., 2020).

The analysis of the composite band image proved to be an interesting tool to visualize the results in an integrated way, making it possible to observe regions where all frequencies indicated similar characteristics for distinct seabed types (for example, darker tones indicating more muddy bed, and whiter tones indicating more sandy sediment). Nevertheless, the observation of the RGB mosaic also allowed to recognize differences in acoustic responses in some regions, such as in the redder region near to the sample sites Doce 01 and Doce 02–indicating higher backscatter values for the lower frequency, and a lighter bluish region near to sample site Doce 11)–coloring that indicates higher backscatter values for the combination of the higher frequency bands (280 and 400 kHz). Quantifying the gain due the use of such methodology is beyond the scope of this work, but the result brings new information about the spectral signature of certain types of bottoms, and helps to consolidate the potential to improve the distinction and classification of the seafloor using multispectral backscatter (Costa, 2019).

Results from the GLCM analysis provide information that may aid in interpreting geoacoustic properties of the seabed. Jackson and Briggs (1992) found that the scattering interface may be a dominant factor controlling backscatter in sandy bottoms, while volume homogeneity is increasingly important in bottoms of higher mud (and lower medium-coarse sand) content. Here, the highest GLCM variance values were observed at the oblique and elongated sandy features, indicating higher uncalibrated backscatter values due to mosaic grey-level values—particularly at the 280 kHz operating frequency (Figure 10). Hughes Clarke (2015) also noted an increase in backscatter intensity over sand sheets with increasing frequency, yet retained uncertainty regarding textural classification, as the bottom ruggedness measurement depends on the acoustic wavelength. Closer to the river mouth, differences in GLCM variance may be associated with other variables that influence the backscatter response, such as the thickness and homogeneity of the muddy subsurface deposit. In addition to informing classification, the GLCM variance contributed to visual interpretations, providing results that assist with backscatter analysis and classification of the seabed (Samsudin and Hasan, 2017; Runya et al., 2021).

While the 30–60° angular range is ideal for many seabed mapping applications (Lucieer et al., 2018), the utility of angular sector may change with frequency, and therefore, the extraction of backscatter values for different angular sectors may provide increased information with which to inform seabed classification. The 1–10° sector contains the most specular reflection, and here, appeared to show low contrast between different seabed types. The 11–20° data provided useful information at the fine sandy region, primarily at the higher frequencies. The angular response curves for this region show low backscatter loss for the first 20°, followed by a higher backscatter loss with increasing incidence angle. The 21–30° angles allowed for distinguishing bottom types using the lower frequencies. The 400 kHz data at these angles allows for visually distinguishing differences in bottom type close to the river mouth, which were not apparent at lower angles. The 31–40° sector showed high contrast between the sandy feature and the rest of the region. While the distinction between fine sand and muddy regions was not always clear using image-based methods, the use of backscatter at different angular sectors and frequencies enabled their discrimination, further indicating the benefits of combined analysis approaches for enhanced seafloor discrimination (Che Hasan et al., 2012; Ierodiaconou et al., 2018).

The ARA results appeared to correspond to the composition of sediment samples. Those with a large proportion of mud showed a greater rate of backscatter decrease with incidence angle than samples describing coarser sediments (Figure 12), and the greatest differences between frequencies. This trend was most apparent from the 400 kHz data, which is consistent with lower uncalibrated backscatter values mapped by the image-based results and with the acoustic influence of the muddy deposit homogeneity. Gaida et al. (2018) observed an increasingly flat angular response curve at 400 kHz in soft sediments compared to lower frequencies, which they attributed to an increased sensitivity to roughness at higher frequencies. Here, the angular response curves observed for sandy sample sites (Doce 07, 09, 11, and 12) exhibited flatter shape, similar to acoustic model curves (Applied Physics Laboratory, 1994). Combining with the other results, it can be concluded here that the flatter the angular response curve, the recognition of medium/coarse sandy bottoms is better observed in the image textural analysis, corroborating the greater influence of the scattering interface. Deriving consistent trends between backscatter and frequency over a range of seabed types remains difficult, yet Trzcinska et al. (2021) summarized backscatter values and angular response curves extracted from different studies over a range of seabed types. Multispectral backscatter facilitated this relationship with the results presented in the angular response curves: higher backscatter values were observed for higher frequencies in fine sand bottoms here, which appeared in a similar shape to the very fine sand in Trzcinska et al. (2021). Higher backscatter values for 280 kHz on a medium sand bottom here also appeared similar to their findings for sandy gravel (Trzcinska et al., 2021). For muddy bottoms, this relationship is made difficult by the suspected presence of a heterogeneous subsurface. For a better understanding of the results, it would certainly be interesting to include some subsurface analysis (e.g., sub bottom profiler or even subsurface sediment samples) and observe the ARA results at smaller angular intervals with more ground truthing.

The seafloor variability within the swath is one of the difficulties related to the angular range analysis, but the use of RGB image segmentation before clustering the variables of both approaches was an effective way to combine the methods and explore the dataset more adequately. The combination of such approaches applied to multispectral backscatter is still scarce in the scientific literature, but showed potential to better explore the advantages of each method using multifrequency dataset, even not involving bathymetric predictors.