Multi-Frequency, Multi-Sonar Mapping of Shallow Habitats—Efficacy and Management Implications in the National Marine Park of Zakynthos, Greece

Abstract

In this work, multibeam echosounder (MBES) and dual frequency sidescan sonar (SSS) data are combined to map the shallow (5–100 m) benthic habitats of the National Marine Park of Zakynthos (NMPZ), Greece, a Marine Protected Area (MPA). NMPZ hosts extensive prairies of the protected Mediterranean phanerogams Posidonia oceanica and Cymodocea nodosa, as well as reefs and sandbanks. Seafloor characterization is achieved using the multi-frequency acoustic backscatter of: (a) the two simultaneous frequencies of the SSS (100 and 400 kHz) and (b) the MBES (180 kHz), as well as the MBES bathymetry. Overall, these high-resolution datasets cover an area of 84 km² with ground coverage varying from 50% to 100%. Image texture, terrain and backscatter angular response analyses are applied to the above, to extract a range of statistical features. Those have different spatial densities and so they are combined through an object-based approach based on the full-coverage 100-kHz SSS mosaic. Supervised classification is applied to data models composed of operationally meaningful combinations between the above features, reflecting single-sonar or multi-sonar mapping scenarios. Classification results are validated against a detailed expert interpretation habitat map making use of extensive ground-truth data. The relative gain of one system or one feature extraction method or another are thoroughly examined. The frequency-dependent separation of benthic habitats showcases the potentials of multi-frequency backscatter and bathymetry from different sonars, improving evidence-based interpretations of shallow benthic habitats.

1. Introduction

The cumulative human-induced threats on marine biodiversity and ecosystem goods and services are increasing worldwide while habitat loss and degradation, pollution, overfishing and climate change effects are classified among the most important threats that affect the marine realm [1,2]. Yet, the successful and effective protection and conservation of the marine ecosystem largely depends on our knowledge regarding the spatial extent, the geographical range and the ecological characteristics of the resource or the habitat of interest [3]. Thus, accurate and high-resolution habitat mapping can be considered as a prerequisite in spatial marine management and conservation planning. In recent years the urgent need for accurate habitat maps of the seafloor has been largely acknowledged within the framework of the relative legislation such as the European Habitats Directive (92/43/EEC) and Marine Strategy Framework Directive (MSFD; 2008/56/EC). Likewise, marine spatial planning, systematic conservation planning, biodiversity and marine resources management also demand good and comprehensive spatio-temporal knowledge of the seascape and the habitat types [4,5].

Swath mapping echosounder is the preferred and most efficient way to rapidly map benthic habitats at all ranges and depths [6,7,8]. A significant part of the survey efforts rests on acquiring acoustic backscatter intensity of the seafloor with sidescan sonars (SSS) and, more recently, multibeam echosounders (MBES). Variations in surveying approaches limit repeatability and comparability of the mapping results, but MBES offers the opportunity of building standardized procedures for automated seafloor mapping [9,10]. There is a growing literature that makes use of MBES for benthic habitat mapping [11,12], taking advantage of the wide spectrum of information that it provides, including co-registered backscatter and bathymetry and angular backscatter intensity. Bathymetric and backscatter texture features are long being used for quantitative analysis of the seabed components [13,14], but also Backscatter Angular Response Analysis is gaining popularity over the years and tends to be a popular procedure for seafloor characterization [15,16,17,18,19]. Among the disadvantages of traditional MBES is the relatively low resolution that they provide, although newer systems (e.g., interferometric sonars) tend to compete with SSS systems in spatial sample density in shallow waters. The above limits severely the backscatter texture content of the recorded seafloor, while angular responses are traditionally analyzed in regard to a full swath strip [9,17], limiting the spatial density of the extracted statistical parameters. In addition, MBES systems lack in seafloor coverage in shallow waters, as their swath width is proportional (3–4 times) to the depth. Significant technology advances of MBES rest on standardization of the collection and processing of calibrated backscatter data and the introduction of multispectral backscatter [9,20,21].

The acoustic signature of seagrass, and in particular Posidonia oceanica, is well known, e.g., [22,23], modulated by the efficiency of photosynthesis [24], i.e., related to the time of day and season at which the seagrass was mapped and how individual responses from algae combine over meadows of different densities [25]. It is also affected by gas content [26] and general health [27], with clear differences between frequencies and imaging angles [28]. Information necessary to reliably identify seagrass, reefs and thin or thick sedimentary layers is limited by the angle of imaging, the spatial resolution and other intrinsic parameters such as beam aperture and pulse length but most importantly by the frequency e.g., [29]. However, the frequency dependent visibility of benthic habitats, including Posidonia, is not yet known and experiments on multi-frequency dependencies on benthic habitats is a recent field of evolving research, with recent examples being [30,31,32,33]. As a general rule, the choice of imaging frequency not only directly influences the resolution of the system and the time it takes to map an area, but it is also related to the acoustic response (backscatter intensity and texture) of a seafloor type.

Development of multispectral MBES, although highly promising, is still in its infancy and therefore multi-frequency datasets are still mostly acquired by successive seafloor scans of the same area with different transducers, which is generally operationally inefficient or even impossible [21,33], and/or on using multiple MBES systems on a single survey [20]. On the contrary, the majority of modern SSS have a dual-frequency capacity and so they can offer a multispectral perspective to seafloor mapping [34,35]. Recently, SSS systems with the capability of simultaneously producing three frequencies have been developed and demonstrated [36]. Additionally, SSS systems offer much higher resolution backscatter images of the seafloor, suitable for detailed habitat mapping through texture analysis [11,25,37,38,39]. SSS can acquire data from very high grazing angles, not limited to specific beam angles as MBES are. This makes it more suitable for shallow waters. What is considered an important drawback of SSS is its lack of secondary information regarding bathymetry and angular dependency of returning echoes. This way, SSS measurements, as opposed to MBES, cannot be accurately corrected for radiometric and geometric artifacts [40], and it is also hard to investigate angular backscatter dependencies for habitat discrimination.

As MBES advances lead to vastly improved maps of benthic habitats, MBES tends to replace SSS for seafloor mapping. Although there is work comparing these two systems [11], there seems to be a research gap concerning if and how SSS and MBES could be used jointly to complement each other in terms of spatial and spectral resolution and mapping coverage, especially in shallow habitats. The present work uses a single frequency MBES and a dual frequency SSS spanning three frequencies (100, 180 and 400 kHz) during a single survey. Texture analysis is applied to each frequency’s respective backscatter mosaic, thus extracting a range of textural features from each. In addition, the MBES bathymetry is processed to extract bathymetric features, while well accepted approaches are applied for the extraction of statistics on angular backscatter response curves (ARCs). The above features constitute the basis for building supervised seabed classifications. All the pieces of the above dataset needed to be analyzed under a common scale, achieving stationary information in all extents of the area. This common scale was achieved through a segmentation process on the 100 kHz mosaic (being the only full seabed coverage dataset), which provided the means for an Object Based Image Analysis (OBIA). This way, regardless of the different data densities and non-stationarity of some data sources, all their features can be joined into common minute regions covering the full seafloor scene.

This work is an effort to exploit the synergistic use of MBES and SSS as complementary sources of information for improved benthic habitat mapping with emphasis on Mediterranean seagrass habitats. By weighing their relative importance for habitat discrimination, we identify optimal mapping strategies, informing future surveys and extending interpretation of the huge amount of existing datasets. At the same time, a high-detail habitat map of a very important Mediterranean MPA, the National Marine Park of Zakynthos (NMPZ), is presented for the first time.

The remainder of this paper is organized as follows: Section 2 describes the ecological and geological settings of the study area (NMPZ); Section 3 describes the data acquisition materials and methods and the detailed benthic habitat cartographic work derived by the expert interpretation of the available data. Section 4, Section 5, Section 6 and Section 7 present the data processing methods, regarding the feature extraction and supervised seafloor classification approaches. Finally, Section 8 provides an in-depth analysis of the results along with a discussion for general guidelines and future perspectives.

2. Study Area

The National Marine Park of Zakynthos (NMPZ) (Figure 1) represents the first marine protected area (MPA) in the Mediterranean for the protection of sea turtles, including no-take zones. It was established in 1999 with the primary goal of protecting and managing one of the most important nesting rookeries of the loggerhead sea turtles Caretta caretta in the region [41]. The MPA mainly encompasses four habitat types that are all included in the EU Habitats Directive 92/43/EEC; namely, Posidonia oceanica (EU habitat code—EUhc: 1120—“Posidonia beds”), rocky reefs (EUhc: 1170—“Reefs”), Sandbanks (EUhc 1110), soft substrates (EUhc 119A—“Unvegetated soft bottoms” and 119B—“Vegetated soft bottoms”) and marine caves (EUhc: 8330) westernmost of the NMPZ [42]. Previous research efforts on marine habitats in the MPA of the NMPZ were rather limited [42,43]; they predominantly mapped seagrass, especially Posidonia oceanica meadows [44,45]. The offshore geology of the NMPZ area is characterized by diapiric salt intrusions (salt domes) forming extensive reefs, up to 20 m higher than the soft bottom, as shown in the chirp SBP stratigraphy of Figure 2. As in many cases worldwide, Zakynthos’ salt diapirs are associated with a wider hydrocarbon province, apparent in various free hydrocarbon seeps in the area off and on shore [46]. The same SBP profiles, on the northern shallow parts of NMPZ, clearly show sandbanks, as well as soft sediments on top of two sedimentary basins of significant thickness (>25 m), formed between the salt dome outcrops, on its deeper parts. Posidonia oceanica beds are typically found between 10 m and 35 m deep, and mostly develop on shallow sandbanks (Figure 2).

3. Data Acquisition and Manual Habitat Mapping

The NMPZ has been mapped in a multi-platform manner (Figure 1), encompassing acoustic systems ranging from SSS and MBES to a chirp Sub-Bottom Profiler (SBP). An extensive supplementary visual seafloor inspection survey, using a towed underwater camera, offered the means for ground-truthing and validation of the mapped seafloor components. During the coupled SSS-MBES survey, 36 survey lines, with a total length of 262 km, were carried out using an Edgetech 4200 SP dual-frequency SSS, transmitting simultaneously at 100 and 400 kHz and a dual-head MBES Elac Nautic Seabeam 1185, transmitting at 180 kHz, covering a total area of 84 km² (see Figure 3). To meet MBES survey requirements, vessel motion was acquired using an SMC IMU-108 MRU system while sound velocity (SV) in the water column was recorded using Valeport MiniSVS as a keel SV probe and Valeport MIDAS SVP as an SV profiler, realizing down-casts on a frequent basis.

The MBES data was acquired through Hypack 2016, while Elac’s HydroStar software was used as the inertial beam-forming, navigation and attitude system. A Real Time Kinematics (RTK) GPS was used to obtain a 10 cm lateral positioning accuracy. To avoid changes in source levels, pulse lengths, receiver gains and directivity patterns, all above parameters were set constant throughout the survey. Automatic Gain Control (AGC) and the Time Varied Gain (TVG) were turned off. Pulse length and frequency were constant through the swath, with no sector dependence, and automatic adjustment of the pulse length as a function of depth was set off to 1 ms. The swath coverage sector was changing between 130° and 150° depending on depth (130° at areas deeper than 50 m) with an along- and across-ship beam of 1.5°, eventually producing 106 or 126 equiangular soundings per swath respectively, corresponding to data point densities ranging from 0.4 to 3.2 points/m.

Adapting to the needs of a time-limited field-work schedule and to the inherent low coverage of MBES in shallow waters, as it is proportional to depth, full coverage data was only feasible through the 100 kHz SSS dataset (see Figure 3), whose slant range was set to 200 m. The lower the frequency the lower the sound attenuation in the water column and thus low frequency SSS can achieve full coverage seabed surveys in a fraction of time than with higher frequency SSS. Given the above, the effective slant range of the 400 kHz SSS was found to be 100 m, thus leaving systematic data gaps between the successive survey-lines, corresponding to about 50% of the surveyed area (see Figure 3). The equivalent mean data coverage of the MBES was about 70%, ranging from <30% shallower than −20 m and >90% deeper than −30 m. A 20-km long ground-truth survey (see Figure 1) was concurrently run using a SeaViewer 950 analog towed camera with an attached GoPro to validate acoustic bottom types. The ground-truth survey was planned so that all distinct seafloor types, having been recognized by expert interpreters, were visually inspected (see Figure 4). Due to the lack of an underwater positioning system, adequate positioning of the underwater imagery was performed at post-processing using landmarks, collecting salient bottom features, evident both in the video and in the backscatter mosaic, calculating the layback of the tow-camera to the vessel and accounting for the mid-distances using standard spline interpolation.

Careful interpretation of all the available marine geoacoustical and ground-truth datasets led to a detailed manual habitat map of the NMPZ, on the basis of the Natura habitat classification scheme (Figure 4). The Natura classes considered are the Posidonia oceanica beds (Natura code category 1120), “Sandbanks” (1110), “Reefs” (1170), “Unvegetated soft bottoms” (119A) and “Vegetated soft bottoms” (119B). Vegetated soft bottoms, in the case of NMPZ, corresponds to Cymodocea nodosa beds. Posidonia beds have further been separated into low relief (“meadows”) and high relief (“matte”) (Figure 4). The division of Natura code category 1120 was done because, although Posidonia meadows and matte are morphological deviations of the same species, they constitute totally different habitat types, both in backscatter texture and in ecological importance.

4. Multifrequency Backscatter Mosaics and Texture Analysis Features

4.1. Mosaicking MBES (180 kHz) and Dual Frequency (100 and 400 kHz) SSS Backscatter data

In this work, we used both a textural and mean angular response approach. The textural approach is less sensitive to absolute backscatter calibration [20] but it heavily relies on radiometric and geometric corrections to produce consistent backscatter responses across the full swaths. Thus, these corrections were considered both for the MBES and SSS backscatter mosaics. MBES backscatter data were derived from the “sidescan” option of the Elac system, i.e., digitization of the received beams from the receive beamformer into 2048 samples. The backscatter values have been normalized to an arbitrarily chosen reference level to remove the angular dependencies, here set to the middle-range incidence angle of 45°.

Acoustic returns from MBES produce intense geometric and radiometric artefacts in the backscatter mosaics, which in turn affect the accuracy of facies delineation. To identify and correct radiometric and geometric MBES backscatter artifacts, we used the Geocoder tool, part of the Hypack 2016 suite. Geocoder, e.g., [18,40,47], uses advanced slant range correction, beam pattern compensation, speckle removal, overlap and feathering options. This leads to a drastically improved MBES mosaic of the NMPZ. Individual MBES time series were slant-range corrected based on the bathymetric model of the survey area, replacing the flat bottom assumption with an estimation of the average slope in the across-track direction of each ping. The effective incident angle (based on the bathymetric model) was corrected using Lambert’s law, taking in account the changes in volume scattering and slope. A residual beam pattern correction was removed by calculating a moving average of angular responses (in our case 100 pings). Overlapping used a data quality indicator, attributing low values to samples closer to and further from the nadir and high values to mid-ranges, allowing selection of the highest quality samples. A speckle noise filter and a feathering algorithm were used to reduce nadir artifacts and the seam artifact between overlapping acquisition lines. All the above processes led to a drastically improved mosaic (see Figure 3), with a pixel size of 50 × 50 cm, suitable for producing spatially consistent results through any image segmentation or classification approach. Further details about the processes implemented in Geocoder can be found in [48]. In [40], the importance of Geocoder’s acoustic backscatter corrections describe for texture-based seafloor characterization. No absolute dB values are given in the backscatter mosaic, as recommended by the GeoHab community, e.g., [9].

The raw SSS records underwent radiometric and geometric corrections using the Isis software (Triton Imaging Inc., Capitola, CA, USA). Beam-pattern correction and ping energy level normalization were followed with slant range correction and ping by ping spatial registration. Mosaicking of both frequencies of the SSS data used the Triton Map software (Triton Imaging Inc., Capitola, CA, USA), producing mosaic with 50 cm × 50 cm pixel size (see Figure 3).

In Figure 3 the final seafloor backscatter mosaics at each individual frequency are shown, where clear differences in area coverage are evident between them. A closer comparison between the three individual frequencies is later presented in Figure 11, revealing clear differences in backscatter texture and intensities between the various bottom types, especially regarding the seagrass cover.

4.2. Image Texture Features

Various texture-based sonar seafloor classification software have been popular among researchers during the recent years, including academic ones, such as TexAn [49,50,51] and SonarClass [52], and commercial ones, such as QTC Side and Multi-view [53,54], Geocoder [18] and others [55,56]. Image texture relies on the relative positions of grey levels within an image neighborhood, and describes the roughness or the smoothness, the variability or the homogeneity and the repeatability or the randomness of different areas on the swath sonar backscatter [51]. In the present work, a set of 11 textural features (Figure 5) was extracted through SonarClass from each one of the three backscatter mosaics.

SonarClass is a Matlab tool used with success in various benthic habitat mapping and target detection works [25,37,40,47,52,57,58,59,60]. SonarClass employs three feature extraction algorithms, based namely on first order grey-level statistics, metrics on grey level co-occurrence matrices (GLCMs) originally defined in [61] and 2D power spectrum specifications, all fully described in [62] (Figure 5). The textural features were extracted from each backscatter mosaic using 30 × 30-m windows, with a 20 m sliding step, leading to 11 feature maps with 10 m spatial resolution each. The 30 × 30 m patch size was decided after visual interpretation of the main image textures in the mosaics, so that primitive texture elements from different bottom types (e.g., seagrass patches, minute reefs, sand ripples) could completely fit into a single image patch. Smaller image patches would lead to incomplete analysis of bottom textures while bigger patches would increase the likelihood that different bottom types are included in a single patch.

The nadir patches were excluded from analysis as they are affected by strong acoustic artifacts [7,8], making them inappropriate for use in automated classification processes. Likewise, transmitting two frequencies with the same system but with one at half the ping rate of the other, induced a persistent interference, expressed as a parallel to the nadir and at a distance equal to the high frequency (400 kHz) slant range linear strip noise (Figure 6). In this work, automated localization of both above linear artifacts was performed using Independent Component Analysis (ICA) decomposition of the feature vector. ICA holds the advantage of separating irregularities in the image scene, usually captured in just one independent component [58]. The selection of the IC corresponding to the artifact relied on visual inspection of the ICs. The theory and principles of ICA have first been applied in SSS images in [58], to which readers are directed for further details. Figure 6 shows an example of artifact exclusion through ICA decomposition.

5. MBES Bathymetric Features

Post processing and registration of the data was achieved using Hypack 2016 (Hysweep module) software package. The data processing included the following stages:

• Data Review: The sonar position and attitude available with the navigation survey lines were closely examined and the sonar data were compensated for heave-pitch-roll, tide-and-draft as well as for sound velocity profile information. Tidal variations were estimated using the Real Time Kinematics (RTK) GPS receiver, installed on the research vessel. True Heave, pitch and roll corrections were acquired through the MRU.
• Swath-by-swath editing: Automated geometric filters removed outliers present as spikes in the data.
• CUBE: Developed by UNH-CCOM [63], CUBE (Combined Uncertainty and Bathymetric Estimator) provides a near-automated editing of multibeam data, allowing for a rapid turn-around of data. CUBE is an error-model based, direct DTM generator, that estimates the depth plus a confidence interval directly on each node point of a bathymetric grid. The output DTM had a 2 m pixel size.

Bathymetric features were extracted with the Benthic Terrain Modeler (BTM) toolbox for ArcMap 10.1 [64]. The BTM toolbox contains a set of tools that allow users to create grids of bathymetric position index (BPI), standardized BPIs, slope, and terrain ruggedness from an input bathymetric data set. In this study, 10 bathymetric features were extracted (see Figure 5), namely depth, slope, aspect, Eastness, Northness, ruggedness (at two scales: 9 m and 21 m) and BPI (at three scales: 25 m, 50 m and 150 m). The BPI algorithm compares each cell’s elevation to the mean elevation of the surrounding cells within a rectangle, annulus (donut shaped area) or circle. The ruggedness is the ratio of the surface area to the planar area and it effectively captures variability in slope and aspect into a single measure. The above-mentioned scales, corresponding to simple circular neighborhoods, were decided upon trial and error, under the restrictions posed in minimum and maximum effective scales by the data resolution and the wide gaps in the shallower areas (higher radius would result in extensive no-data areas) and trying to hardly follow a factorial sequence multiplied by the minimum scale to be analyzed.

6. MBES Backscatter Angular Response Features

Backscatter angular response approaches use the angular response curves (ARCs) to extract statistical parameters that describe their overall structure and sometimes even directly link them to modelled responses of certain sediment granulometries. This work followed the approach of [65]. ARC statistics are extracted from seafloor patches, defined as average per 2.5° angular bin (range of grazing angles) of several consecutive sonar pings, set here to 30. The stacked angular responses are then divided in three angular ranges: Near, Far and Outer. Finally, slopes and intercepts are extracted for each range per channel (port and starboard) (Figure 7a). In our case ARCs are calculated directly from the recorded raw data files and they are corrected for transmission losses (spreading and absorption) and for the volume scattering due to the ensonified footprint following [66].

Factor Analysis application to the ARCs, using the consecutive stacks of pings as cases and the respective angular bins as variables, clearly revealed 3 angular ranges with highly correlated responses. Those ranges, defined as Near: 0–25°, Far: 25–45° and Outer: 45–65° (Figure 7b), nearly validates the work in [65], which defined them as Near: 0–25°, Far: 25–55° and Outer: 55–85°.This way, a calibration is realized, especially in the Far and Outer ranges, so that those ranges are better adapted to our dataset. Following the above findings, slopes and intercepts where extracted from the ranges identified by Factor Analysis application, finally constructing an ARC derived feature vector of 6 features (Figure 5).

7. Supervised Classification

7.1. Facing Coverage and Resolution Inconsistencies: 100 kHz SSS Mosaic Segmentation for Object Based Classification

Object Based Image Analysis (OBIA), has become a standard practice for gaining detailed delineations of habitats through swath sonar images. It involves similar pixels being grouped into connected segments (objects), classified using some intrinsic information. It has mostly been applied to MBES data [12,15,16,18,56,67,68], where the objects are extracted through segmentation of MBES backscatter mosaics and the available information (bathymetric, backscatter and angular backscatter responses or acoustic modeling parameters) is averaged within each segment. Segmentation practices are also used in SSS analyses using backscatter information for habitat mapping [38,47], and more recently for SBP data [47]. In this study, object-oriented analysis offers the means to successfully combine MBES and dual-frequency SSS data into an integrated classification procedure, although the dataset has systematic gaps throughout the surveyed area and there are different resolutions between one system or another. The 100-kHz SSS system was the only one with full coverage of the survey–area; its mosaic thus formed the basis for object extraction through image segmentation techniques. The average value of each derivative within each segment constitute a new feature vector, composed of all 49 features considered in this study, fully covering the survey-area. It is assumed that at least some data derivative points from each available system are located within each segment, to assure that the size of the segments should be bigger than the data gaps that may exist between any pair of consecutive survey-lines. Thus, any segmentation technique used should allow for scale adjustment.

To meet the goal of image segmentation with selectable scales, the Segmentation by Weighted Aggregation, as implemented in the SWA tool [69] was used. This method constructs a pyramid of graphs, which represents aggregates of pixels of similar properties. At fine levels these aggregates represent regions of homogenous intensities while at coarser levels coherent regions of more complex third order image texture measures. The 100-kHz SSS greyscale mosaic was down-sampled to 1000 × 700 pixels in order to make the segmentation algorithm run efficiently. The intensity “fine” and “average” coefficients were set to 8 and 5 respectively, the intensity variance minimum scale was set to 3 with a coefficient equal to 0.5 and the 6th scale was chosen as the most suitable for the needs of this work. Figure 8 shows the result of the segmentation process overlaid on and compared side-by-side to the 100-kHz SSS mosaic. Examples showing the average response of selected texture, terrain and angular response analysis features within each segment are shown in Figure 3.

7.2. Data Models: Systems’ Features Fusion and Feature Selection

Each data model is a group of features associated with a single system, a single data analysis method or a meaningful combination of the two. The choice of features in each data model was driven by a meaningful and operationally realistic data acquisition and processing concept. For instance, data models m1–m3 regard using a single dual frequency SSS, where each frequency or both are used for textural analyses. Data models m4–m8 respect a data acquisition concept in which a single MBES is used, but with the application of different levels of data processing i.e., terrain, backscatter texture or backscatter angular response analysis, a combination of the above or all together. Data model m9 concerns both SSS and MBES used together, but with just their backscatter mosaics being analyzed, which is a concept adapted to exploring any possible rate of improvement when increasing the spectral bands of the backscatter texture analysis. Finally, data model m10 assumes that both MBES and SSS are used together and all available feature extraction methods are applied to them. For each group, a feature selection algorithm is used to detect the features (no matter which system they are derived from) that offer the best classification accuracy possible (see Section 7.3). A wrapper feature selection method [70] has been chosen, evaluating the feature sets in each data model, by using a learning scheme (Naïve Bayes here—see Section 7.3) that offer the best information gain. In this way, the most important features in each group are determined to reduce the overall complexity of the feature space that in turn will lead to producing simpler decision rules during the supervised classification process. Feature selection is an attractive option to validate which feature extraction methods offer better habitat discriminations in each concept and weight their relative importance per data model as well as per habitat type. The final data models considered after the feature selection stage are presented in Figure 5, organized by feature extraction method and swath mapping system. It is important to mention that at least one derivative from each feature extraction method has been selected through this process, implying that each one is offering information gain in the habitat discrimination process and they are all worth being considered.

7.3. Supervised Classification: Training, Validation and Accuracies

Training of the classifiers was performed using the averaged features (49 in total) from a neighborhood of 100 m radius around each ground-truth data point. More precisely, the underwater videos have been converted into georeferenced snapshots every 100 m intervals, which constituted the ground-truth points after expert interpretation of the visible bottom types.

The validation dataset was composed of the averaged features per segment with each segment characterized after the most frequent bottom class included in its extents, as derived by the detailed manual habitat map presented in Figure 4. This object-averaged version of the manual habitat map actually is the best possible classification that could be achieved through the given mapping scale, governed by the size of the segments, and so it is used as the basis for evaluating the classification accuracy achieved through supervised classification of each data model.

A simple Naïve Bayes classifier was selected to avoid over-segmentation (overfitting) of the feature space and to extract simple decision rules. Although using Random Forest (RF), Support Vector Machines (SVM), or Neural Network Classifiers (NNC) or other sophisticated classifiers could potentially produce higher classification precisions, Naïve Bayes produces much more generalized decision rules, which constitute a fair baseline for comparison between different data models, avoiding any feature space over-segmentation. Simple decision rules need well separated classes in the Euclidean space, which is a desired hypothesis sought to be elaborated in this work. Classification was performed on the raw derivative points without performing any object-wise averaging. Objects were considered on a later stage, during which the majority class of the derivative points falling inside each segment was considered to construct the final classification maps per data model.

Classification accuracy assessment per data model was achieved through the Cohen’s Kappa coefficient, estimated by the confusion matrix of the predicted and the manual classes for all segments. Classification accuracy assessment per Habitat type was performed using the “recall” measure. This is to unveil which system or which system’s derived components worked better for each habitat type. This is important to exploit the inherent capacity of each system or feature extraction method towards the discrimination of certain bottom types.

8. Results

8.1. Angular Backscatter Responses

Figure 9 shows a characteristic ARC per Habitat type, which actually is the average of all ARC points (middle swath points per side) falling inside each habitat type region of the manual map in Figure 4. Most ARCs are well separated, with the most acoustically reflective substrates (P. oceanica beds and reefs) following a Lambertian distribution (isotropic scattering), and the soft bottoms (sandbanks and unvegetated soft bottoms) showing strong specular components with a slight slope across the ARC. Vegetated soft bottoms (C. nodosa beds) presents an intermediate ARC, owed to a soft bottom with sparse and low vegetation of relatively higher reflectivity. The volume scattering in this seabed type plays an important role in the ARC’s shape, as the extensive rhizomes just below the seafloor, supporting the seagrass leaves, are strong scatterers at 180 kHz and similar frequencies, (e.g., [24]), which penetrates a few centimeters in the sediment (e.g., [71,72]). Yet, the two different sub-categories of Posidonia oceanica beds, low-relief meadows and high-relief matte, have almost identical ARCs. This is interesting as, although their textural appearance and relief are totally different in the backscatter and bathymetry data, it is their physical properties that governs the shape of the ARCs, i.e., both are variations of the same species. Almost the same happens with “Sandbanks” and “Unvegetated soft bottoms”, which also have very similar ARCs. However, they can still be separated at Far ranges.

8.2. Frequency Dependent Separation of Seagrasses

Frequency separation of habitats is a key field of research, aiming to introduce multi-spectral satellite approaches to seafloor mapping. In this work, multi-spectral backscatter data came from the multi-source swath-sonar systems used. The SSS system had a dual-frequency capacity producing instantaneously co-registered maps at 100 kHz and 400 kHz. This offered the option to explore the level of information revealed with each frequency. Figure 10 shows such an example in Area A in Figure 1, where P. oceanica, C. nodosa and sandy beds laying in proximity are ensonified using 100 kHz and 400 kHz frequencies. Those two frequencies are further combined to produce a false-color composite, with green and blue channels corresponding to 100 kHz and 400 kHz respectively. Comparably higher backscatter at 100 kHz is shown in green shades in the pseudo-color image, higher 400 kHz intensities in blue and similar intensities in grey. The maximum benefit of using both frequencies is that the C. nodosa beds are remarkably better separated. This is because the slightly submerged rhizomes of C. nodosa exhibits a much higher contrast at 100 kHz than at 400 kHz, as lower-frequency signals penetrate some centimeters in the sediments [71,72]. This intensity difference between the color channels highlights the respective areas in the composite image. Conversely, soft sediments exhibit higher backscatter strength at 400 kHz than at 100 kHz, because of increased volume scatter (due to penetration) in the latter, thus producing blue shades in the respective regions. Finally, Posidonia oceanica, being a strong acoustic reflector with an isotropic scattering behavior, produces similar backscatter intensities in both frequencies, leading to grey shades in the pseudo-color image.

SSS backscatter mosaics cannot be directly compared to MBES because of the lack of accurate co-registration. Moreover, the different pulse geometries and lengths between the two systems affect the volume scattering components of the seabed and thus a false-color image would not truly reflect the seabed response just on different frequencies but also to other transmit parameters. Nonetheless, each one of the three available backscatter mosaics exhibits different levels of separation between the various bottom types, and so they all worth to be included in a robust seafloor classification scheme, as it was realized in this work. Figure 11 shows how different habitat types P. oceanica meadows, C. nodosa beds and Sandbanks seem when mapped using the different transmit frequencies and swath mapping sonars.

8.3. Relative Importance of MBES and SSS Data Models for Acoustic Habitat Mapping. Allies or Not?

Figure 12 presents the classification accuracy assessment results of the ten data models. An 85% Kappa coefficient was achieved in data model m10, i.e., when all MBES and SSS features have been used, for which the classification map is shown in Figure 13. When only MBES features were used (data model m8), the classification accuracy reached about 77% while when only SSS features (regarding both 100 and 400 kHz frequencies) were used (data model m3), it was 10 percent lower, about 67%. The above supports the rationale of this work that was to fill the gap concerning if and how SSS and MBES could complement each other towards improved habitat mapping in shallow habitats, as in the NMPZ case, a single-system survey would offer results from 10% to 20% less accurate. MBES backscatter information alone (either mosaic texture or angular range analyses) have been shown to provide the worse results (data models m5 (texture): 39% and m7 (angular responses): 50%). However, if combined with the auxiliary bathymetric analysis features, the above offered the best single-system automated habitat mapping results (data models m8). Regarding the contribution of multi-spectra backscatter mosaics to the classification performance through the texture information, it was shown that individual frequencies provided very poor results, with Kappa coefficients around or below 50% (data models m1, m2 and m5). Combination of the 100 kHz and 400 kHz SSS mosaic texture features increased classification accuracy more than 15% in comparison to the single frequency results (data model m3: 67%). Using all three mosaics allowed for a slight improvement of less than 5% over the latter data model (data model m9: 70%). Finally, the contribution of the angular response analysis features to the overall performance of the classifier was very low and more precisely about 3%, as indicated by comparing data models m6 and m8 that concerns MBES single-system scenarios, the former without and the latter with backscatter angular responses consideration.

8.4. Classification Accuracy per Habitat Type. Which System is Better for Each Shallow Habitat Type?

Figure 14 presents the classification accuracy of each individual Natura habitat type when SSS and MBES systems are used alone (data models m3 (SSS) and m8 (MBES)) or together (data model m8). MBES is further exploited concerning the added gain when just using bathymetry (data model m4), bathymetry and backscatter texture (data model m6) or a full feature set of bathymetry, backscatter texture and angular responses (data model m8). As a rule, MBES was found to be the best single-system choice for separating most habitat types, with high superiority in rough, high relief ones, i.e., “Reefs”. However, dual-frequency SSS (data model m3) was better for textured seabed types e.g., “Posidonia beds—matte” or featureless “Sandbanks” and “Unvegetated soft bottoms”. “Vegetated soft bottoms” (C. nodosa beds) were better separated by MBES (data model m8). This latter observation seems odd, as discussed in Section 8.2, as C. nodosa beds have a very characteristic acoustic signature in dual-frequency SSS records. However, it is the only high-intensity/low-ruggedness seabed type, and as such the MBES bathymetry proved equally important to the backscatter texture. “Posidonia beds—low relief” (meadows) were hard to separate using either swath-system, as they were confused locally with other high intensity bottom types (e.g., “Vegetated soft bottoms”, “Posidonia beds—matte” and “Reefs”). Nonetheless, the multi-source swath sonars scenario (data model m10), with conjunction of statistics on MBES and dual frequency SSS data, proved to provide the best possible classification results in almost all bottom types, with significantly better performance in “Unvegetated soft bottoms” and “Posidonia beds—low relief”, worse in “Reefs”, and equally good in “Posidonia beds—matte”.

9. Discussion

A thorough evaluation process clearly showed that automated acoustic habitat mapping can be highly accurate when using multi-platform swath-sonars for shallow habitat mapping. MBES alone proved capable of separating most seafloor types, but it suffered from low coverage in shallow waters. This made lower-frequency SSS a very important complementary system, as its transmit geometries allow for full coverage even in very shallow habitats. We take advantage of the very low grazing angles that low-frequency SSS can achieve to maximize the scanned area. At the same time, we gain dense bathymetric information along with additional backscatter data from the MBES. This combines into four layers of information; three backscatter mosaics at different frequencies and bathymetry. Although SSS and MBES data cannot be precisely co-registered, how important this is for habitat mapping studies needs to be compared with the desired resolution. For each dataset, the accuracy in the positions of sonar measurements on the ground must ideally be smaller than their spatial resolution. In our case, where quite large seafloor segments define the baseline spatial scale for the analysis, positional accuracies are always comparatively higher.

Regarding the system selection per mapping subject and area, SSS offers the fastest possible mappings in terms of survey time, retaining wide swaths in shallow depths at lower frequencies and it preserves high textural content in the acquired images. The above makes it suitable for shallow habitats with highly textured patchy seagrasses on a plain seafloor. At the same time, common dual-frequency SSS offer the opportunity for a multi-spectral perception of the seafloor. While low frequencies show the highest overall contrast between different seafloor types, a consideration of all frequencies permits for an improved interpretation of subtle sea-floor features. Multi-spectra backscatter proved to be advantageous for the separation of seabed types where the volume scattering versus the isotopic scattering is of importance, i.e., when there is some shallow layering just below the seafloor as in the case of sparse seagrasses with extended rhizomes (e.g., Cymodocea nodosa). However, due to the lack of auxiliary depth and incidence angle information, traditional SSS offers less capabilities for the application of a full range of geometric and radiometric corrections and the classification of habitats in rough terrains (e.g., reefs and steep slopes) is less accurate. In the contrary, despite the lower resolution and coverage that common MBES provide, they have been proven to offer the best trade-off between habitat discrimination and survey cost, as it is a many in one system. Given the recent advances in MBES Mills-Cross technology, offering snippet imagery with centimeter-accuracy and wideband operation e.g., [32], one can get much more predictive results than the ones derived by the older Elac system used in this study. One should consider though that it would have been impossible for any MBES to provide a full coverage mapping without the object-oriented analysis that was made feasible through the full coverage low frequency SSS mosaic. The above arguments imply that while MBES takes the lead in habitat mapping works at present, SSS still remains an invaluable tool for higher resolutions and coverages in shallow waters that is always worth employing.

Regarding the habitat composition in NMPZ, measurements show that “Posidonia beds” cover 27% (22.5 km²) of the NMPZ seafloor, C. nodosa beds (“Vegetated soft bottoms”) 6% (4.8 km²), “Reefs” 9% (7.8 km²) and soft sediments 58% (48 km²), highlighting the ecological importance of the area. For the first time, a very detailed benthic habitat map was produced of this area, with classifications suitable for any future spatial planning or monitoring works in the NMPZ.

Effective protection and conservation of the marine realm and the goods and services that it provides to humanity, should follow an ecosystem-based management (EBM) approach [73,74]. However, in order to successfully reach the EBM objectives detailed, accurate and easy to use habitat maps are required by both the scientists for the needs of integrated ecological assessments and the decision makers for marine spatial planning and management purposes, respectively [75,76,77]. Therefore, the approach that we propose here for high accuracy habitat mapping through automated procedures may facilitate the enhancement of management effectiveness of MPAs in protecting species and habitats in the EBM context. In this respect, fine scale habitat mapping, in conjunction with the spatial information on human activities, can further enhance our knowledge in managing the conflicts between human uses and conservation objectives while mitigating the cumulative impacts of stressors [77,78] both inside and outside the MPAs. Yet such information is valuable in identifying, designing and prioritizing networks of MPAs at various spatial scales [79,80]. Still, the thorough knowledge of the seafloor topography through accurate habitat maps can also be of use in a fisheries management context [81].

It is it well known that seagrasses and especially Posidonia oceanica meadows are declining in the Mediterranean at alarming rates due to human activities and climate change and thus they have become one of the main targets of protection and management [82]. Yet in the Mediterranean Sea the Natura 2000 marine sites have been designated under this protection status mainly due to the presence of the protected seagrass Posidonia oceanica (priority habitat—Habitat Type 1120: “P. oceanica beds”) [5]. Thus, the accurate mapping of this habitat-forming species will allow precise detection of any change in the distribution range size, demographic status and mortality events of P. oceanica meadows through space and time and consequently to assist conservation efforts.

10. Conclusions

A fully automated shallow habitat classification approach was proposed, making use of multi-source, high-resolution swath sonars, for improved seafloor classifications in a total area of 84 km². We exploit the efficacy of the combined use of MBES and SSS for gaining more acoustic backscatter spectral bands in a single survey, as well as dense bathymetric and angular backscatter responses information. Integration of MBES and SSS is an operationally efficient practice, quite common in marine geophysical surveys, to gain maximum information from a single survey. Such a multi-platform approach can offer at least three acoustic backscatter spectral bands on a single survey with varying data resolutions and mapping scales, supplemented by bathymetry and angular backscatter response information. Combination of the above systems proved to offer a good trade-off between survey effort and information gain, mainly in shallow waters. The poor textural content and coverage in shallow waters of the MBES was supplemented by the higher equivalents of SSS, offering an improved classification of the bottom classes. As a rule, MBES although providing poorer coverage in shallow waters, it provided the best single-system seafloor classification accuracies, but it needed to be combined with an object-oriented analysis made available through the full coverage SSS mosaic. On the other end, individual frequency backscatter mosaics, no matter which system they were derived from, were unable to provide adequate habitat discrimination. This was much improved when all backscatter frequencies were used in synergy, but it was still below the accuracy levels achieved with the fully integrated data model, combining both SSS and MBES data derivatives. NMPZ, as an MPA of high local importance, being the first no-take zone in Greece, stood as a perfect case study to showcase the above findings. It is a shallow, high diversity MPA, hosting extensive fields of the protected phanerogams P. oceanica and C. nodosa that have been fully and in detail mapped, offering the background for local management stakeholders to set the standards for an optimized marine spatial planning of the area.