New Feature Classes for Acoustic Habitat Mapping—A Multibeam Echosounder Point Cloud Analysis for Mapping Submerged Aquatic Vegetation (SAV)

Abstract

A new method for multibeam echosounder (MBES) data analysis is presented with the aim of improving habitat mapping, especially when considering submerged aquatic vegetation (SAV). MBES data were acquired with 400 kHz in 1–8 m water depth with a spatial resolution in the decimeter scale. The survey area was known to be populated with the seagrass Zostera marina and the bathymetric soundings were highly influenced by this habitat. The depth values often coincide with the canopy of the seagrass. Instead of classifying the data with a digital terrain model and the given derivatives, we derive predictive features from the native point cloud of the MBES soundings in a similar way to terrestrial LiDAR data analysis. We calculated the eigenvalues to derive nine characteristic features, which include linearity, planarity, and sphericity. The features were calculated for each sounding within a cylindrical neighborhood of 0.5 m radius and holding 88 neighboring soundings, on average, during our survey. The occurrence of seagrass was ground-truthed by divers and aerial photography. A data model was constructed and we applied a random forest machine learning supervised classification to predict between the two cases of “seafloor” and “vegetation”. Prediction by linearity, planarity, and sphericity resulted in 88.5% prediction accuracy. After constructing the higher-order eigenvalue derivatives and having the nine features available, the model resulted in 96% prediction accuracy. This study outlines for the first time that valuable feature classes can be derived from MBES point clouds—an approach that could substantially improve bathymetric measurements and habitat mapping.

1. Introduction

1.1. Multibeam Echosounder Evolution

Multibeam echosounders (MBES) were invented in the 1960s with the first commercial devices for deep water measurements entering the market in the 1970s [1,2,3]. In addition to military purposes, multibeam’s objective was to perform accurate and efficient depth soundings, and to contribute to a better interpretation of the seabed that came along with the broader coverage, three-dimensional (3D) surface generation, and visualization that the multibeam can offer. Within the last two decades, great progress has been made in antenna array design, signal processing, and bottom detection, as well as in positioning and motion sensing, which enables highly accurate and high resolution depth measurements to be obtained today, especially at shallow water depths [4]. Early multibeam systems were able to form 16 beams only across a 45-degree swath, resulting in 16 soundings or points per ping. With the progress in digital signal processing and storage capabilities, modern systems are able to acquire significantly more data; today, we deal with swath angles up to 160° and the invention of the Mills-cross array led to an increased sounding density with typical values between 256 and 512 beams. Recently, multi-swath systems with sector-coded frequency modulation have been introduced, which has doubled the number of soundings possible. So-called “near-field focusing”, leading to narrower beams in the sonar near field, is another improvement [5]. With sophisticated dualhead systems, the swath coverage can be increased and one ping cycle can produce up to 1600 soundings. With ping rates exceeding 50 Hz in very shallow water, this translates into 80.000 points per second and it provides a high sounding density per square meter. Such high point density holds superior information regarding grid data; however, this information is often erased when dealing with calculated DTMs.

In the last few years, it has become possible to record water column data from all MBES systems as beam backscatter time series. In general, such recordings are valuable but they remain computationally very expensive and difficult to handle [6,7,8], and they are merely recorded by the authorities and hydrographic agencies. Bottom detection algorithms (BDA) with multi-detect capabilities [9] present a compromise as they do not neglect water column imaging (WCI) but still operate economically. In most cases, hydrographic databases are still missing full WCI recordings, and they are beyond the scope of this study for practical reasons.

1.2. MBES Seafloor Classification and Habitat Mapping

Today, MBES present powerful devices for seafloor classification and habitat mapping using both bathymetry and co-located backscatter, with the latter sometimes being referred to as MBES snippet backscatter when many samples are calculated within one beam. There was little ambition from the manufacturers until a decade ago to make extensive use of such data. One reason was that the surveying market did not require such data. Backscatter was considered as ‘nice to have’ but not adding value to nautical charting. This has changed over recent years, which is partly due to increasing user interest and the application of MBES for seafloor classification and habitat mapping.

Various techniques that also integrate MBES backscatter, e.g., for angular range analyses, have been proposed, which makes MBES very valuable seafloor classification tools [10,11,12]. In addition to the bathymetry data itself, many first and second order geomorphometric derivatives, such as slope, aspect ratio, bathymetric positioning index, and higher order derivatives [13] are calculated for better habitat assessments [14]. Micallef et al. [15] combined such bathymetric features with textuaral backscatter analysis of MBES data to achieve a semi-automated method for habitat mapping. Recent work clearly outlines that increasing the number of features (gridded bathymetry, backscatter, and derivatives) can highly improve habitat modeling; however, Stephens and Diesing [16] showed that the results of supervised classification approaches can be further improved when only using a subset of the most predictive features. Even more discriminative power for habitat mapping can be expected from the arising multispectral MBES capabilities [17,18,19,20].

There is no doubt that many types of submerged aquatic vegetation (SAV), e.g., eelgrass or seagrass can be acoustically detected. Detailed lab studies showed the dedicated backscattering response from the seagrass’ tissue in the range of 0.5–2.5 kHz for three different seagrass species [21]. A recent review study by Gumusay et al. [22] discusses seagrass detection with various acoustic systems between ~1 and 1000 kHz.

1.3. Point Cloud Analysis

Point cloud (PCL) analysis is quite common for terrestrial research in the LIDAR community, where one important task is to separate between the ground and non-ground returns (i.e., returns from vegetation, buildings, etc.) [23]. A commonly used and very promising approach is “eigen-features”—values that are computed from eigenvalues of the covariance matrix of the local neighborhood [24]. There are several other approaches, such as height-based approaches [25,26]. Yan et al. provides an overview of terrestrial LiDAR classifications methods [27].

LiDAR and MBES are based on the same working principle: the measurement of angle and two way travel time of a transmitted electromagnetic/acoustic pulse. In both cases, the returns are stored as raw and ungridded PCLs. However, in marine sciences, MBES PCLs are normally gridded before further analysis takes place, which may come along with information loss.

1.4. Submerged Aquatic Vegetation

SAV in the Ocean can occur down to a depth of 200 m in the photic zone and it occurs in coastal areas all over the globe. The coastal areas of the Western Baltic Sea (at water depths between 0.6 and 7.6 m) are frequently populated with the eelgrass Zostera marina [28]. Seagrass meadows are one of the most important habitats on Earth in terms of their various ecological functions: they act as nursery grounds for fish and a food source for marine animals [29,30], function as a sediment trap for coastal erosion mitigation and as a CO₂ sink [31,32], and they can even improve the water quality [33]. However, while seagrass improves water quality, the growth of the seagrass itself is dependent on the water quality; in particular, the maximum distribution depth of seagrass strongly depends on the light permeability of the water column, which in turn depends on the water eutrophication. Consequently, the maximum distribution depth of seagrass is an indicator for the water quality [34].

1.5. MBES Zostera Marina Survey and Study Site Introduction

We employed a hydrographic MBES as our device of choice in the coastal area of the Baltic Sea. In this region, water depths between 1 and 7 m are frequently populated with the eelgrass Zostera marina. The seafloor at the study site mainly consists of fine grained sediments ranging from clay to coarse sand and it has gentle slopes with along-shore sand bars.

1.6. Study Objectives

In this study, we claim the meaningfulness of MBES PCLs analysis, which has, until now, not been applied for MBES data, except for some special cases, such as least depth wreck detection and seabed gas seepage [6,35,36]. We propose that new features that are directly derived from ungridded, native PCLs provide additional discriminatory power for habitat mapping, as these PCLs usually hold more information than downsampled gridded data.

LiDAR and MBES using similar operation principles and the separation of returns from the crown of trees and returns from the ground in LiDAR data is furthermore a similar problem to the separation of seagrass and seafloor echoes in MBES data. Consequently, we adapted and tested the segmentation and classification methods that were developed by the LiDAR community for the analysis of MBES PCLs.

Therefore, the aim of this study is twofold. Firstly, we will introduce PCL analysis for MBES data and we adapted classification methods from the terrestrial community for the application in marine science, in particular for SAV mapping. Secondly, we integrated aerial imagery for the training of a data driven machine learning model for SAV detection.

2. Methods

2.1. MBES Survey and Ground-Truthing

The data were collected during a survey with the small vessel Klaashahn from the Institute for Baltic Sea Research (IOW) from 4th to 8th September 2017 in a coastal area of the Western Baltic Sea, which is located off Heidkate, Germany (Figure 1). The Klaashahn is a low draft vessel able to operate in waters as shallow as about 1 m water depth. The MBES data were recorded using Hypack 2017 and a Norbit iWBMSe that was pole-mounted on the portside of the vessel at 0.45 m depth. This system has an integrated Applanix SurfMaster inertial navigation and attitude system. The MBES device was operated with a 80 kHz width chirp at a mean frequency of 400 kHz, forming 256 beams with individual beam resolution of 0.9°, and in equiangle mode, with an opening angle of 150°. The upper gate of the bottom detector was placed between 1 and 1.5 m above the seafloor to ensure that the seagrass canopy was within its search range. The sound velocity profiles for raytracing were collected with an AML Oceanographic Base X sound velocity probe. The average survey speed was 5 kn.

Research divers from the scientific diving center of the University of Kiel performed ground-truthing at the same time. Additionally, punctual optical ground-truthing was done with an underwater video camera, equipped with a CMOS image sensor with 1000 TV lines, during the survey. The captured frozen images have a resolution of 720 × 288 pixels. The aerial images of this area were recorded during a LiDAR campaign on the 1st September 2017. These images have a resolution in the decimeter range and they were used for an aided interpretation of the shallow MBES data and for the training of the supervised data driven machine learning SAV model.

2.2. MBES Data Processing and Feature Generation

We used the open-source software MB-System by MBARI (Monterey Bay Aquarium Research Institute) for the common MBES post-processing steps, like patch test calibration, raytracing, tidal reduction, and outlier removal of the MBES data [37]. The cleaned ungridded three-dimensional (3D) PCL was then exported as xyz-ASCII in cartesian UTM 32 North coordinates with reference to the mean sea level, for further processing.

Next, we defined neighborhoods for all points in the PCL. We chose a cylindrical neighborhood with a radius of 0.5 m, as it enables feature computation that is based on the height distribution [38], resulting in an average neighborhood of 88 points. The neighborhoods of each point in the PCL were computed with the very fast open source machine learning library mlpack [39].

To calculate the following features, we developed a C++ program for computation. The matrix operations and eigenvalue calculations within this program were performed using the open-source library Eigen version 3.3, which is a popular library that is used, amongst others, by Google’s TensorFlow and at the ATLAS experiment at CERN.

The first feature we calculated was the vertical displacement dz between each point of the PCL and the lowest point in its neighborhood. We first computed the covariance matrix Cov of each point neighborhood to derive the eigen-features:

$$Cov=\frac{{\displaystyle \sum }\left(A_i-\overline{A_i}\right)\cdot \left(A_j-\overline{A_j}\right)}{N-1}$$

(1)

where A_i contains the N observations in the neighborhood and A_i holds the mean of the respective dimension. Subsequently, we computed the 3 eigenvalues λ_i of the covariance matrices Cov:

$$Cov\cdot v_i={\lambda }_i\cdot v_i$$

(2)

where v_i are the corresponding eigenvectors. Assuming that the eigenvalues are sorted as λ₁ ≥ λ₂ ≥ λ₃ ≥ 0, they can be used to calculate the following eigen-features [40]:

$$\begin{array}{l}Linearity: L_{\lambda }=\frac{{\lambda }_1-{\lambda }_2}{{\lambda }_1}\\ Planarity: P_{\lambda }=\frac{{\lambda }_2-{\lambda }_3}{{\lambda }_1}\\ Sphericity: S_{\lambda }=\frac{{\lambda }_3}{{\lambda }_1}\\ Omnivariance: O_{\lambda }= \sqrt[3]{{\lambda }_1{\lambda }_2{\lambda }_3}\\ Anisotropy: A_{\lambda }=\frac{{\lambda }_1-{\lambda }_3}{{\lambda }_1}\\ Change of curvature: C_{\lambda }=\frac{{\lambda }_3}{{\lambda }_1+{\lambda }_2+{\lambda }_3}\end{array}$$

(3)

The first three of these features describe the local geometry [41]. The linearity is a good feature for the recognition of lined up points, but it also increases at the borders of the PCL. The flat regions in the PCL have high planarity, whereas sphericity detects high-curvature points [24]. In addition to these eigen-features and dz, we computed features based on a local plane fitted to the neighborhoods. We used the RANSAC (RANdom SAmple Consensus) algorithm, which can better handle the outliners than a normal least square fit, to fit these planes. The RANSAC algorithm randomly chooses s points (s = 3 for a plane) from the dataset and then estimates the model (plane). Subsequently, it computes the distances of all the datapoints to the model. If a distance is greater than a defined threshold δ, then the corresponding datapoint is an outliner and erased. All other datapoints are saved. This process is repeated until the biggest set of datapoints, the so-called “consensus set”, is found. The consensus set is then used to fit the model via a least square fit. We used δ = 0.1 m for our computations, which is about the ten-fold SONAR range accuracy of the MBES that we used. These planes were used to derive the distance dp between the queued point and the plane, the summed distances dsum of all points in the neighborhood and the plane, and the angle φ between the normal of the local plane and the vertical.

In total, we computed nine features for each point in the MBES PCL. We used a random forests (RF) supervised machine learning approach to classify, and differentiate between, returns from SAV and from the seafloor. The RFs are ensemble-based classifiers that build multiple decision trees on the subsets of the input variables [42]. The user, depending on the number of variables and classes, can select the number of decision trees. The class predicted by most decision trees is output as the result. The RFs were set up using the mlpack library [39].

These computations on ungridded MBES PCLs require high computational power. Consequently, we restricted the test of our new approach to the subsamples of different vegetation densities. The sub-areas were identified in the aerial images (Figure 1): SA1 contains a dense seagrass meadow, SA2 comprises sandy seafloor with distributed seagrass patches (Figure 1, inlets), and SA3 consists of mainly unvegetated seafloor.

We used sub-area SA2 for the training and testing of the RFs, since it is comprised of dense seagrass patches and unvegetated seafloor. For this purpose, we randomly distributed sample points with a minimum distance of 1 m over the area of SA2 with QGIS 3.6 Noosa. These points were then assigned to the existing points in the PCL and manually labeled as either seafloor or SAV return using the aerial imagery and underwater photographs. Subsequently, we randomly split these sample points into a training and test dataset using mlpack_preprocess_split, resulting in 2932 (80%) training points and 732 (20%) test points (Figure 2a,b).

Afterwards, we trained several RFs with different feature ensembles and number of trees and the training dataset. The RF with the best performance (ratio between correctly assigned points in the test dataset and the total number of test points) was then used to classify the PCLs of SA1, SA2, and SA3 into SAV and seafloor returns.

Subsequently, we used the complete PCLs to compute a DSM (digital surface model), and only the points that were classified as seafloor to compute a DTM (digital terrain model). The computations of DSM and DTM were performed using a splines-with-tension algorithm in gmt surface (gmt version 5.4.3; [43]). Afterwards, we computed a canopy height model (CHM) by subtracting the DTM from the DSM [26] in QGIS 3.6 Noosa.

3. Results

3.1. Seafloor, Habitat Description, and Acoustic Bottom Mis-Detection

Figure 1 presents the edited and gridded bathymetry of our survey. The seabed is covered by sand with sand bars aligned parallel to the coastline. The seabed is frequently populated with the seagrass Zostera marina as visible in the aerial photography (Figure 1) and ground-truthed by scuba divers (Figure 1, inset). The maximum height of the seagrass was determined by the divers to be 0.6 m. Occasionally, blue mussels occur within the seagrass meadows and, in some places, the seagrass is accompanied by red algae up to 0.1–0.2 m high. No cobbles or boulders were reported in this area. Comparison between the aerial images and the shallow MBES data highlights that acoustic depth anomalies quantitatively match with areas that are covered by seagrass, even in places where the seagrass occurs in bunches with loose leaf density. Thorough inspection of PCL shows that the acoustic energy was not completely backscattered by the canopy of the seagrass, but sufficient energy was obviously received from the interior of the meadows, and even from the seafloor (Figure 3), which prompts bottom misdetection of various heights between real seafloor and seagrass canopy by the BDA. This resulted in acoustic shadowing and rougher bathymetry in regions that are covered by seagrass when compared to uncovered regions, as visible in the point cloud (Figure 3), and indicated by the bathymetric DTM (Figure 1) and computed terrain roughness index (measure for elevation differences of neighboring cells [44]; not shown).

3.2. Calculated Features from Point Clouds

Nine features were derived from the native PCL of SA2 after Equation (3) and they are shown in Figure 4. The highest values of Linearity L_λ mainly correspond with the edges of the PCL and with gaps in the dataset, while the values of Planarity P_λ are a maximum, where the seafloor is flat and uncovered. The canopies and edges of the SAV patches are associated with lower P_λ values. In contrast, sphericity S_λ is clearly enhanced at the canopy of the SAV patches and it is zero for the flat unvegetated seafloor. dz is clearly influenced by the seagrass fields; contrastingly, the values of dp bear much less resemblance to the seagrass. φ mainly delineates the border of the SAV patches and has similar values (around 0°) for both the tops of the patches and the unvegetated seafloor. The feature O_λ again clearly responds on the SAV patches, whereas A_λ also shows their extent to a lesser degree. Last but not least, C_λ again clearly highlights SAV. Figure 4 highlights that P_λ, S_λ, dz, φ, O_λ, and C_λ are the features that show the best response to SAV.

3.3. SAV Identification in PCLs with RF

Several RF runs with different tree numbers and feature ensembles were performed. We found that increasing the amount of trees above 20 for ensemble sizes between six to nine features, and above 10 for ensemble sizes between three and five features did not affect the training accuracy.

Table 1. Performance of RFs for SAV classification in PCL against manually flagged SAV returns, for SA2 and different feature ensembles, ordered by performance from top to bottom.

Feature Ensemble:	RF Performance [%]:
Lλ, Pλ, Sλ, dz, dp, φ, Oλ, Aλ, Cλ	96.1
Pλ, Sλ, dz, φ, Oλ, Cλ	95.9
Lλ, Pλ, Sλ, dz	95.6
Lλ, Pλ, Sλ	88.5
Lλ, Pλ, Sλ, Oλ, Aλ, Cλ	88.2
Sλ, dp, φ	86.9
Lλ, Pλ, Sλ, φ	86.5

shows the test outcomes for the different feature ensembles on SA2. Together, the most common eigenfeatures L_λ, P_λ, and S_λ are able to correctly predict 88.5% of the manually-classified SA2. The feature ensemble (P_λ, dz, φ, O_λ, C_λ), which by eye has the highest response to SAV, has a prediction accuracy of 95.9% with regard to our manual classification. However, the RF build-up that includes all nine features achieved the best performance with 96.0% prediction accuracy. Table 1 also shows that it is not only the number of features that is decisive for the classification result, but the choice of meaningful features.

However, we used the RF, trained with all nine features, to classify the PCLs of SA1, SA2, and SA3, since it has the best test performance. The classified SAV returns of SA2 are superimposed on the corresponding aerial image in Figure 5c,d, which illustrates the high agreement between the two. Figure 6 shows the computed CHMs. The maximum SAV canopy heights vary between 0.49 m (SA3) and 0.77 m (SA1), which is in good agreement with ground-truthing by divers who reported a maximum height of 0.6 m, but could only do sparse, pointwise sampling.

4. Discussion

In this study, we introduced a PCL analysis of MBES soundings for SAV mapping similar to that which has been successfully conducted for terrestrial LiDAR data [26]. Our adaption of LiDAR processing methods is based on the fact that LiDAR and MBES are similar in their measurement principle, i.e., both based on the time of flight or two way travel time and angle measurement of a signal pulse to provide coordinate triplets (x,y,z) of the target, sometimes combined with intensity (LiDAR) or backscatter strength (MBES).

Diesing et al. have reviewed the similarities between terrestrial and marine remote sensing and habitat mapping methodology itself [14], but only for image based and gridded data. To the best of our knowledge, there is no marine study yet that uses native PCLs for the investigation SAV or for habitat detection and classification, as presented here. The very recent review study regarding SAV mapping by Gumusay et al. [22] lists several devices for seagrass detection and mapping, amongst which are single beam echosounders, acoustic doppler current profilers (ADCP), side scan sonars (SSS), MBES, and even low-cost recreational fishfinder-SSS. However, the authors only describe image-based approaches (SSS- and MBES-BS-mosaics, echograms, and WCI data) and features that are derived from DTMs, but do not mention the use of native MBES generated PCLs.

4.1. Comparison to other Remote Sensing SAV Methods

In a recent study, Kruss et al. [45] used backscatter mosaics from a MBES to classify different seafloor types and identify SAV. Additionally, they used the near nadir beams (+/−10°) to estimate the macro-algae height. However, they only tested their approach for estimating the SAV height, and not on the outer beams. If this approach does not work on the outer beams, it would imply that denser track lines are required, as is already the case when using MBES in shallow water. Contrastingly, we did not observe an angle-dependence of our approach and it seems to work on the outer beams, so MBES trackline planning does not have to be adapted in order to apply the approach that is outlined in this study.

Fakiris et al. suggests the possibility of increasing the trackline spacing for economic habitat mapping in shallow water, by using a combination of SSS and MBES [19]. They used two SSSs with different frequencies and MBES backscatter to increase spectral information for an improved habitat classification. The tracklines were planned with gaps between the MBES swaths. The improved classification is thereby performed by MBES-derived features, which are extrapolated using segmented SSSs between the tracklines. Fakiris et al. [19] used their approach for habitat classification and mapping, however their approach is not able to assess SAV canopy heights.

The use of WCI also yields possibilities for the detection of seagrass, as shown by Hamana and Komatsu [46], who developed a method for real-time classification of seagrass on flat bottoms using a MBES. However, the disadvantages of this technique are that not all MBES devices support WCI, it requires huge amounts of memory, and it is difficult to analyse the outer beams due to elongated echo responses near the seabed.

Backscatter analysis is performed on mosaics (equi-distant gridded data), which enables statistical analysis that improves habitat classification, such as grey level co-occurrence matrix [47]. These backscatter analyses are normally carried out on SSS data, hence do not provide depth and SAV height information, and are therefore not suitable for biomass estimations.

As SAV growth only occurs within the photic zone, optical approaches may also be applicable for its mapping. These optical approaches can make use of satellite, aircraft, or unmanned aerial vehicle-borne aerial images. In a recent study, Poursanidis et al. [48] used WorldView-2 data for the seafloor classification and mapping of Posidonia oceanica in Greece. The authors were also able to determine the maximum growing depth of this seagrass species. One disadvantage of such optical approaches is that their success strongly depends on the water turbidity. Its application for the mapping of deeper seagrass in turbid waters, like the Baltic Sea, is therefore not feasible. Additionally, satellite data have lower resolutions than near-range hydro-acoustic data and are hence not suitable to determine small-scale changes of seagrass distribution that may be caused by coastal protection measures, such as sand nourishment, wave breaker, or harbor construction activities.

4.2. Robustness of Point Cloud Analyses from an Environmental Perspective

MBES point clouds ideally represent the seabed, although it should be noted that the respective BDA producing them may get confused by acoustic features, such as SAV, which vertically rise from the seabed (Figure 3 and Figure 7).

In the Baltic Sea, the investigated seagrass Zostera marina grows on sandy or muddy sediments, with no apparent preference for either of the two substrates (pers. comm. T. Reusch, F. Weinberger, P. Schubert). The morphology of sedimentary settings in the Baltic Sea is generally gentle and vertically-aligned features, such as SAV, physically and acoustically stand out from the flat seabed by up to 0.68 m (Figure 3, ground-truthed by divers). We suggest four possible candidates as vertically-aligned features on an at-least-decimeter scale that might be prone to being picked by the BDA as “the bottom”: (a) benthic fish schools, (b) cobbles/boulders, (c) true morphology (such as sand bars), and (d) other SAV, where (b) and (c) would be considered to be part of a DTM. Additionally, BDA bottom mis-detections may also arise from acoustic interference, bottom reverberation, uncompensated water column anomalies, or special seabed geoacoustic properties, like fluid mud [49], very soft sediment, or gas [50,51].

4.3. Possible Improvement of MBES Point Cloud Analyses from a Technical Perspective

The evaluation of PCLs and calculation of the presented features are both very susceptible to outliers [24]. As the features are calculated from overlapping swaths, it is important to “iron out the creases” that come along with improper calibration especially in shallow water [52]. First, we suggest that special care should be taken during data acquisition regarding motion compensation, in order to obtain very high quality measurements that are a prerequisite for a successful PCL analysis. Secondly, the data need to be thoroughly patch-tested, because outliers are not smoothed by any gridding operation, and overlapping soundings with depth offsets would immediately impact the calculation of features such as dz.

We did not observe systematic dependence of misdetections on the beam angle, and thus the prevailing type of bottom detection, i.e., inner beam amplitude or outer beam zero phase crossing detection [53]. This means that the BDA reliably identifies SAV irrespective of amplitude or phase detection being applied for bottom detection. For future studies, we suggest BDAs may need environmental training to account for SAV and other environmental settings that can cause misdetections (Figure 7).

The data suffered in some bad weather and we sailed our survey without the centimeter-accurate real time kinematic GPS, which is usually available in coastal seas today. The motion sensor that was used is not designed for ocean swell, and we attribute some of the data anomalies to the motion sensor performance. We assume that the most limiting factor for the data is that it was acquired without a multi-detect mode [9,54] in the BDA, as it was not yet available at the time of collection. Since our model’s prediction provides very good results with this data, we expect even better performance with a survey optimized, as discussed previously. PCL density could also be improved with multi-swath or larger antennas with narrower beams, synthetic aperture sonar technology, and lower survey speeds. Of course, the inclusion of the backscatter might be a very promising task to be considered in future studies, but it was left out in this study in order to demonstrate what is feasible only using PCL analysis.

We used a desktop computer with a 64 bit Intel^® Core^TM 4 core CPU and 32 GB RAM for the computations that are presented in this study. With this computer, it was not possible to compute the neighborhoods for the complete MBES dataset visible in Figure 1, using the range_search tool of mlpack [39] due to a RAM overflow. However, we are not sure if this was really a computer issue or a software issue because of uneconomic memory usage within the code. This problem may be overcome by using an state-of-the-art computer with more RAM or by using another software for neighborhood analysis.

4.4. Model Adaptations and Algorithmic Issues, and Performance

The search radius r for feature calculation was set to 0.5 m in the xy plane, resulting in approximately 88 soundings in each neighborhood cylinder. While three points in the neighborhood are sufficient for the computation of a local plane and derived features, the choice of the neighborhood radius is a crucial point when computing features that are based on the covariance matrix [40]. The main problem is that the neighborhood radius is manually chosen and cannot account for different point density in the PCL [41]. This second reason is particularly important when analysing MBES PCLs that cover a wider depth range. To deal with this issue, Blomley et al. [40] suggest the use of a weighted covariance instead of the normal one, and Demantke et al. [41] developed an algorithm for finding the best neighborhood size. For further application of the MBES PCL features, it is certainly reasonable to test these approaches, especially if the MBES dataset comprises several depth ranges. While the calculated features are not depth-dependent, the size of the neighborhood may need to be adjusted to the mean depth.

The performance of a RF inherently depends on the accuracy of the user-defined training and test cases. In our study, it was difficult to find the exact transition between the seagrass and seafloor returns in the MBES PCL when manually flagging the seagrass points. This was especially the case at the border of the seagrass patches and meadows, where the transition from the seafloor to the seagrass is smooth. It became even more difficult when the local bathymetry was uneven. Consequently, we recommend the use of aerial images to support the manual classification of the training and test cases. Different operators will more likely-than-not find different borderlines around the seagrass patches and meadows [55]. With this in mind, we did not expect a 100% detection rate when conducting the RF test runs. There are also other supervised and unsupervised classification methods, such as k-means, support vector machine, neural networks, and naive bayes [22], which we have not tested here, but may achieve better results.

5. Conclusions and Outlook

SAV is commonly found in the photic zone; however, optic remote sensors struggle to detect it in deeper waters. Therefore, acoustic methods are still required in order to map the spatial distribution of SAV in the ocean.

To date, most acoustic seafloor classification work has concentrated on backscatter and bathymetry data grids, and derived features thereof (mainly first and second order DTM derivatives). Common analysis is performed with grey level co-occurrence matrices, segmented images, or modelled or heuristic angular range behavior. Here, we introduce another acoustic feature class, and make use of the strength of MBES point clouds, in order to identify SAV extent and canopy height. This is reasonable, because all of the previously described features are likely to be frequency-dependent and as the backscattering strength of SAV is not yet known in detail.

Here, we introduce PCL analyses for MBES data as an alternative to working with gridded bathymetry. Our approach is based on the spatial properties of the neighborhoods of each point within the PCL, by computing the so-called “eigen-features” and features based on a plane, fitted to the neighborhoods. Subsequently, we classified the PCL points into SAV and seafloor returns with a trained random forest on the features. With this approach, we achieved accuracies of up to 96.1%. As this was the first approach of this kind, and, given a non-ideal dataset, we suggest that even better results are possible with this method. The classified PCL can be split into SAV and true seafloor returns to compute a canopy height and a terrain model. Thus, our new approach is not only suitable for SAV detection, but also for the assessment and monitoring of SAV biomass using MBES. How inclusion of the backscatter can reliably improve SAV detection and biomass estimates remains a matter of ongoing research.

We would also like to emphasize that we have only used a regular state-of-the-art MBES and freely available open source software, and therefore our approach does not require any special devices or commercial software.

Our approach shows, in principle, that any normal MBES dataset could be reprocessed in terms of SAV mapping. Hence. we foresee a mutual benefit between the seabed habitat mapping community and the hydrographic surveying and nautical chart production community, with the goal of better understanding how MBES seabed mapping interacts with SAV.

The application of the PCL approach to multi-detect MBES data should, in future, result in better DTMs, even in denser seagrass meadows. Incorporating MBES backscatter as an additional feature would also be advantageous, as backscatter is the most prominent and significant feature in side-scan-based seafloor classification. Additionally, the response of the introduced MBES PCL features to other SAV species should be evaluated.

Our adaption of the point cloud analyses for MBES soundings in this study also highlights the mutual benefit of exchanging post-processing strategies between the LiDAR and MBES data.