Spectral features of dual-frequency multibeam echosounder data for benthic habitat mapping

Automatic methods of seafloor mapping are still in their early stage of development, despite the technical progress made in recent years. A serious imperfection is the limited types of predictor features available for seabed classification. It is therefore desirable to introduce new class of spectral features to benthic habitat mapping. In this study, we introduced eight spectral features of a rough seafloor surface that were indicative of better seabed classification. We compared them with traditional secondary features, like terrain variables and textural features. The suitability of 48 variables was tested, and the most important features were identified. The selected variables were used to perform a supervised object-based image analysis using four machine learning algorithms. We found that backscatter was the strongest predictor, followed by several spectral features from bathymetry that appeared more predictive than bathymetry itself. The highest overall accuracy of predictive model reached approximately 86% using the support vector machine classifier. The innovative results of this study suggest further application of the spectral features for predictive benthic habitat mapping, including research based on multi-frequency multibeam echosounder datasets. The utilisation of spectral features derived from bathymetry provide an important step towards more accurate maps of benthic habitats and seabed sedi-


Introduction
The ocean floor is the least explored surface of Earth. At present, it is estimated that less than 15% of the seafloor has been mapped in detail. On the other hand, the surfaces of the Moon and Mars have been mapped in significantly greater detail (Jones, 1999). Global initiatives, such as Seabed 2030 of the General Bathymetric Chart of the Oceans group, aim to change this state of knowledge by mapping the entire seabed by the year 2030 (Mayer et al., 2018).
Apart from side-scan based seafloor analyses, remote sensing measurements of the seabed surface have often employed multibeam echosounder (MBES) systems. Originally, MBES equipment was designed to collect measurements of the seafloor bathymetry, which allowed for the generation of digital elevation models (DEMs) of the seabed. In the early 1990s, an MBES was developed that could measure the backscattering strength from the seafloor based on its corresponding properties (Lamarche and Lurton, 2017). Recent recommendations have suggested concurrent acquisition of bathymetry and seafloor backscatter strength (or a related variable), as well as further generation of georeferenced grids of bathymetry and co-registered backscatter mosaics (Schimel et al., 2018).
term 'spectral' has been used as a descriptor of the rough seafloor surface.
Recent benthic habitat mapping studies have underlined the utilisation of different features of MBES bathymetry and backscatter data (Diesing et al., 2016;Lecours et al., 2016;Held and Schneider von Deimling, 2019). Special attention has been paid to the geomorphometric analysis of bathymetry (Goff and Jordan, 1988;Wilson et al., 2007;Micallef et al., 2012;Li et al., 2016;Diesing and Thorsnes, 2018;Gafeira et al., 2018;Lucieer et al., 2018), textural analysis of backscatter (Montereale- Gavazzi et al., 2017;Prampolini et al., 2018;Samsudin and Hasan, 2017), and multi-scale analysis, including the geographic context of MBES datasets (Lecours et al., 2015;Misiuk et al., 2018;Porskamp et al., 2018). For example, Conti et al. (2019) used different feature categories (from bathymetry, texture, optical properties, and object-based shape) for the classification of a cold-water coral mound based on MBES and a high-resolution video mosaic. Additionally, Janowski et al. (2018a) extracted sixteen features of bathymetry and nine object-based features of backscatter to execute the mapping of seabed sediments in the Polish coastal area of the southern Baltic Sea. Furthermore, Rattray et al. (2015) introduced wave exposure as an oceanographic predictor variable for the mapping of high energy temperate reefs in Victoria, Australia.
Recent reviews of literature emphasise the need for new features for benthic habitat mapping (Diesing et al., 2016). It is known that environmental factors, such as light penetration in the water column, primary productivity, hydrodynamics, temperature, salinity, and oxygen concentration determine the distribution of habitats on the seabed (Brown et al., 2011). However, modelled features representing these factors are typically generated in spatial and temporal scales that are very different (e.g. often significantly larger) from those of MBES datasets, thereby limiting their applicability. On the other hand, secondary features extracted from MBES bathymetry and backscatter are directly related to the spatial and temporal extent of their counterparts. The development of new features may allow for a better understanding of the environmental processes occurring on and/or influencing the seabed, and proper application may increase predictive power and improve classification accuracy.
Spectral parameters have been successfully used to classify sediment types using single beam echosounder registrations (Tegowski et al., 2003). Previous research on spectral features of MBES bathymetry has indicated their utility for a detailed description of roughness and seafloor geomorphology, as well as the classification of seafloor sediments. Additionally, they have been used for benthic habitat mapping using Principal Components Analysis to reduce correlated data and the Fuzzy C-means clustering algorithm with a declared number of three classes (Tegowski et al., 2018). In this study, we presented and evaluated eight spectral features derived from bathymetry and applied them for the classification of the seabed using object-based image analysis (OBIA). These parameters originated from two-dimensional fast Fourier transformation (2D FFT). Lyons et al. (2002) described one of the first applications of this method for seabed characterisation with high-resolution, in which the photogrammetric method (stereoscopic photograph) was utilised; a three-dimensional model of the bottom surface was generated using this approach. Application of the 2D FFT allowed for the spatial distribution of the power spectral density of the surface heights to be obtained. The same technique has been applied in several other studies (e.g. Briggs et al., 2005). Moreover, this method was improved upon and applied to the analysis of high-resolution bathymetry from modern hydroacoustic measurements, including the application of a MBES (e.g. Cazenave et al., 2008;Lefebvre et al., 2009). Schönke et al. (2017) used 2D Fourier transform to describe the microroughness of the seabed based on underwater laser line scanning in the southeastern North Sea.
Classification of the seabed substrata and benthic habitats is one of the main tasks necessary for the spatial planning of the marine environment. In addition to providing crucial information for the establishment of Marine Protected Areas, such actions are within the main aims of Descriptor 6 of the Marine Strategy Framework Directive 2008/56/EC, which is related to seafloor integrity. In general, these actions assume the development of standardised methods for seabed mapping and monitoring. Although recent studies include proposals for the working procedures of benthic habitat mapping or the development of habitat classification schemes, diversity of specific environmental conditions (e.g. depth or sediment types), causes that they are still only valid for strict spatial areas (Strong et al., 2019).

Multi-frequency multibeam echosounder studies
A recent trend in benthic habitat mapping is the use of multi-frequency MBES data. Bottom backscattering strength registered by echosounder strongly depends on the frequency of the emitted signal, its true incidence angle, seabed roughness, and geo-acoustic properties of the seafloor. The dependency of backscatter strength on frequency has been observed in laboratory and field studies with various responses from different sediment types exhibited (e.g. Jackson et al., 1986;Urick, 1983). Recent habitat mapping studies have emphasised the use of multi-frequency MBES datasets for better discrimination between seabed types (Feldens et al., 2018;Gaida et al., 2018;Janowski et al., 2018b;Fakiris et al., 2019).
High-frequency pulses allow for the detecting smaller objects and seabed structures; however, such pulses are strongly attenuated, thus limiting the sonar range. Low-frequency signals are not as attenuated; they can penetrate deeper into the sediments below the seafloor, but they are less sensitive to small features and weak boundaries with a slight change in acoustic impedance, such as the boundary between water and mud. Overall, acoustic images (especially of seafloor sediments) that are recorded at several frequencies often provide more information with respect to the physical and biological characteristics of seabed habitats compared with that of a single-frequency (Feldens et al., 2018;Gaida et al., 2018;Janowski et al., 2018b;Fakiris et al., 2019Finer sediments, such as sands and silts, are more sensitive to acoustic frequencies than coarser sediments, such as gravel, shells, or boulders (Jackson et al., 1986;Williams et al., 2009;Hefner et al., 2010;Gaida et al., 2018).
In this study, we focused on the spectral features derived from bathymetry, which was considered independent of the frequency. Most applicable for this approach were MBES, which delivered bathymetry, as well as co-registered and geolocated backscatter of the seabed. In this study, we did not focus on multi-spectral MBES, although we used several frequencies to enlarge our feature space. Our objectives were as follows: (1) to introduce eight spectral features of a rough seafloor surface, (2) to evaluate the importance of spectral features for benthic habitat mapping, and (3) to classify benthic habitats and estimate the accuracy of classification, including the input from spectral features.

Study site
The study site was situated in a shallow area off the Polish coast of the southern Baltic Sea. The water depth ranged between 3.8 and 20.1 m below sea level (Fig. 1). The site was in direct proximity to Slowinski National Park and partially located within a Natura 2000 area, which protects marine areas up to 10 m below sea level. In the research area, six classes of benthic habitats were present, including flat areas of very fine sand with traces of worm burrows (VFS), sands with ripple marks (S), sandy gravels or gravelly sands (SG_GS), boulders covered with a large concentration of Mytilus trossulus bivalves (B), boulders covered with Mytilus trossulus and large patches of red algae (R), and artificial structures (A), such as a shipwreck located in the centre of the study area (Kendzierska, 2009;Tegowski et al., 2009;Janowski et al., 2018b).
The seabed consisted of valleys and crests with depths of approximately 2 m and lengths of dozens of metres to approximately 180 m. The main geomorphologic structure was a large 1100 m × 500 m moraine outcrop located in the centre of the area. Filled with glacial tills, it was covered with areas comprising numerous boulders with red algae vegetation. Such a hard substratum has not been typical along the Polish coast of the Baltic Sea; however, it allowed for conditions necessary for the development of unique benthic communities, such as Mytilus trossulus bivalves and Furcellaria lumbricalis or Polysiphonia fucoides red algae. Previous studies of this site have confirmed its high ecological relevance (Kendzierska, 2009;Tegowski et al., 2009). The area surrounding the moraine shoal appeared inhomogeneous and consisted of sands of different grain sizes and occasional gravel admixtures. The bio-and geodiversity of this area has made it highly suitable for the evaluation of non-invasive research methods of the seafloor, including the determination of different acoustic characteristics.

Data acquisition and processing
MBES datasets were acquired using a NORBIT iWBMS STX system (manufactured by Norbit ASA: Po box 1858, Lade 7440, Trondheim, Norway) mounted on a portable pole on the 'Zelint' research vessel. The MBES device was manufactured especially for use in shallow marine areas and typically reserved for hydroacoustic measurements from 0.2 to 160 m below sea level. At a maximum frequency of 400 kHz, the receiving beam width was 0.9°x 0.9°and allowed for the collection of 512 beams. The MBES had an integrated WaveMaster (manufactured by Applanix: 85 Leek Crescent, Richmond Hill, ON Canada, L4B 3B3) Global Navigation Satellite System/Inertial Navigation System that was supported by Real Time Kinematic/Global Positioning System corrections for precise positioning and altitude measurements. Using the Polish Active Geodetic Network -European Position Determination System NAWGEO service (www.asgeupos.pl), we received real-time positioning with an accuracy of 3 cm horizontally and 5 cm vertically. In this study, the influence of acoustic absorption on the recorded signals was initially ignored; however, it was considered during post-processing. To fulfil our research purposes, the frequency was set to either 150 or 400 kHz, and the swath range covered 150-160°. The maximum ping rate for both frequencies was 20 Hz, and we applied a 200 μs (for 150 kHz) and 500 μs (for 400 kHz) modulated chirp with a bandwidth of 6 and 80 kHz, respectively. Surveys were designed in respect to the systematic collection of five sound velocity profiles. A constant vessel speed 2.83-3.09 m/s was maintained.
The MBES datasets were processed using QPS Qimera 1.6.3 and Fledermaus Geocoder Toolbox 7.8.4 software, which allowed for bathymetry and backscatter data processing, cleaning, and mosaicking. After registration, a patch test was applied. A bathymetric grid with a pixel size of 0.25 m × 0.25 m was calculated for both frequencies.
Because we did not find any significant differences in the MBES measurements recorded at 150 and 400 kHz, we combined the frequencies to obtain the bathymetry from a dense point cloud. The Qimera software allowed for the manual cleaning of any outliers and/or acoustic spikes. Backscatter grids were generated based on beam time series (snippets) with resolutions of 0.75 m and 0.5 m for 150 kHz and 400 kHz, respectively, using a mosaicking method with an Angle Varying Gain (AVG) correction included in the Geocoder engine (Fonseca et al., 2009). The AVG method has been commonly utilised for the correction of MBES angular dependency and to obtain a normalised seafloor backscatter dataset. We applied the default settings of the AVG, which were 'flat' (mode), 'blend' (mosaicking style), and '300' (size of processing window). The flat mode was responsible for reducing backscatter signal noise and smoothing fine variations. The blend mosaicking style was responsible for the management of overlapping MBES swaths. This allowed for the blending of the pixels along the nadir track line of the vessel with other overlapping pixels (Schimel et al., 2018). The window size corresponded to a specific number of consecutive MBES pings considered for AVG correction (e.g. see Parnum and Gavrilov, 2011). All the MBES datasets were extracted as surface floating point files in a Universal Transverse Mercator (zone 33 N) projected coordinate system.
We applied the general workflow for benthic habitat mapping developed by Janowski et al. (2018b) to the MBES data. Hence, we extracted statistical and geomorphometric features of bathymetry, for instance, slope and vertical ruggedness measure (Sappington et al., 2007), bathymetric position index (Wilson et al., 2007), as well as textural features of backscatter, such as different types of grey level cooccurrence matrices (Haralick et al., 1973). Additionally, we derived firstand second-order 2D spectral parameters of bathymetry. A list of all extracted features is presented in Table 1. Features 1-23 were created using the Benthic Terrain Modeler toolbox for ArcGIS (Walbridge et al., 2018), features 24-39 were calculated using algorithms coded in MATLAB, and features 40-48 were created using OBIA workflows in Trimble eCognition software (Janowski et al., 2018b). When possible, we tested various sizes of rectangular moving windows or scales of image-based objects (Table 1), which enabled us to perform multi-scale analysis of geospatial datasets to a certain extent (Misiuk et al., 2018).

Bathymetric 2D spectral parameters
Utilising 2D FFT of the bathymetric grid allowed us to generate eight spectral parameters, which will be discussed in respect to their predictive power for the description of seafloor geomorphology and classification. The eight spectral parameters were zero-order spectral moment (m 0 ), second-order spectral moment (m 2 ), mean frequency (ω 0 ), spectral width (ν 2 ), spectral skewness (γ s ,), spectral skewness defined for central moments (γ s_centr ), quality factor (Q-factor), and fractal dimension (Dfft).

Considerations of 2D seabed spectral parameters
Assuming that the height values of the bottom surface are normally distributed, and the surface is isotropic, the power spectral density can  be expressed as the following power function (Jackson et al., 1986;Jackson and Richardson, 2007): x y is the wave vector of surface inequalities, γ 2 is the exponent of the spectrum, and w 2 is the spectral power of the rough seabed surface expressed in cm 4 . Both the γ 2 and w 2 spectral parameters comprehensively characterise the scale and degree of surface roughness in that they are basic parameters of the physical models of sound scattering at the bottom (APL, 1994). Measurements of the γ 2 using different techniques, such as stereophotography, laser scanning, acoustic scanning, and mechanical stylus scanning, indicated that the value of this parameter ranges from 2.4 to 3.9 (in most instances), whereas the mean value is 3.25 (APL, 1994). The values of the parameters w 2 and γ 2 can be determined from the Fourier spectrum of a rough surface, and such a method was adopted in this study.
The 2D normalised bathymetric cross spectrum s(x,y) can be represented by Fourier transformation as follows: , x y i π K x K y 2 ( , ) x y where K x and K y are the spatial wave numbers expressed in the cycle•m −1 . The result of the transformation is a spatial spectrum of the surface height. The 2D FFT requires subtraction of the average height and trend removal to avoid spectral leakage. Further reduction of the spectral leakage effect requires additional multiplication of the transformed bathymetric surface by a function of the discrete prolate spheroidal sequences spectral window, or in our case, a first-order parameter in which NW (Slepian bandwidth parameter) = two window widths. Our algorithm was performed on the basis of a moving window with the dimensions of 20 m × 20 m or 35 m × 35 m with a 90% overlap. After executing the 2D FFT, we extracted one-dimensional (1D) cross-sections spectra from 0°to 180°(every 5°) for which we calculated spectral parameters.
To find the specific features of the tested surface of the bottom, 1D spectra were parameterized, and the results were averaged. For each of the 37 spectra, spectral moments m r (Clough and Penzien, 1975) were calculated as follows: where r is the order of the moment, ω is the circular frequency, and S (ω) is the density of the power spectrum.
The mean frequency is defined as the following equation: The spectral width is calculated according to the following equation: It is defined as the concentration of the spectral energy density around the mean frequency. The higher the value of average frequency, the lower the parameter value. Furthermore, the parameter value is higher for a wider spectrum and lower when the opposite is true. A very sensitive parameter for changes in the shape of the surface is the spectral skewness defined for central moments (Davidson and Louglin, 2000), which is calculated as follows: Another parameter based on 1D spectral density of power is the Dfft. If the rough surface of the bottom has fractal properties, the ratio between the spectrum S(f) and the frequency f takes the form of an exponential relation for the frequency interval f (Mandelbrot, 1982) as follows: where K is a constant, and β is the exponent of the power function. The spectrum slope is calculated by linear regression. The Dfft is defined as follows: Additionally, the Q-factor, which is a combination of spectral moments, can be also calculated with the following equation: The Q-factor is a measure of spectrum peak 'sharpness'. For the 1D spectra obtained in this way, the eight spectral parameters defined above were calculated. Whereas spectral parameters created with a moving window size of 20 m × 20 m had a resulting pixel resolution of 2 m × 2 m, parameters generated with the larger window size (35 m × 35 m) had a pixel resolution of 3.5 m × 3.5 m.

Ground-truth data acquisition and processing
Ground-truth samples were acquired with a remotely operated vehicle (ROV) and a Van Veen grab sampler. Samples were retrieved during three surveys on 7 September 2018, 20-23 November 2018, and 21-25 January 2019. Based on previous research in this area, as well as backscatter acoustic characteristics, the locations were carefully chosen in a representative way to capture all the properties of the seabed (Tegowski et al., 2009). ROV video recordings were collected in more  than half the sites, allowing for the investigation of the locations and surrounding areas. Sediment samples were collected from 46 sites, and they were classified using the methods of Folk and Ward (1957), as well as the Wentworth method (Wentworth, 1922). From the 57 ground-truth samples, 29 were chosen for the training of supervised classifiers, and 28 were used to test the classification performance. The split was performed based on random single splitting. Our split for the training/test samples was targeted to obtain better prediction performance than that of model fitting. Fig. 2 presents the locations of the ground-truth samples used for training and validation with reference to the MBES backscatter. Note that the A benthic habitat class is omitted in Fig. 2 because of its presence in only one specific site. It was, therefore, classified manually at the conclusion of the supervised classification process, based on the exact location of the shipwreck visible on the MBES bathymetry and ROV video datasets.

Image analysis for predictive habitat mapping
To evaluate the importance of individual features we applied the Boruta feature selection algorithm (Kursa and Rudnicki, 2010) based on the random forest (RF) machine learning algorithm by Breiman (2001). The algorithm belongs to the wrapper feature selection methods that evaluate the performance of a certain model after searching for all possible feature selections. Typically, wrapper methods aim to minimise prediction error, and because of this, they belong to common minimal-optimal feature selection methods as well (Kursa, 2016). The wrapper is implemented in R software using the 'Boruta' library. The wrapper method iteratively evaluates sets of different input features and calculates a Z-score, which is indicative of feature importance. Each evaluation is done by the introduction of other irrelevant features that are treated as a reference for the assessment of the original features.
The Z-score is calculated based on the RF method during the training of the classifier (Breiman, 2001). Based on the feature importance measure, feature selection is performed iteratively, successively removing irrelevant features. To exclude tentative (unallocated) features, the maximum number of Boruta iterations was set to 5000. Although our aim was to identify all relevant features, including weakly relevant ones (Nilsson et al., 2007), we allowed the possibility of refinement if some were correlated. To remove highly correlated features, a correlation analysis was performed in the R software using a 'caret' package. Features with an absolute Pearson's correlation of 0.75 or higher were removed.
In this study, we used Trimble eCognition software to conduct an OBIA. This image processing technique was developed in the 2000s to manage an increasing number of high-resolution remote sensing images containing larger amounts of heterogenous information (Blaschke, 2010). Through a multiresolution segmentation (MS) algorithm, the OBIA merged similar pixels of an image into groups of uniform shapes and sizes (Benz et al., 2004). MS had various parameters that we defined and tested to generate meaningful image objects. Whereas the colour parameter corresponded to the relative values of the MBES backscatter intensity, the associated parameter (shape) was related to the ratio between compactness and smoothness. Compactness referred to the ratio between the segment border length and the square root of the pixel count within. Smoothness was related to the ratio between the border length of the segment and its bounding box (Benz et al., 2004). Both weighted pairs of parameters were determined with values of 0.1 to 0.9, and the total value of each pair was 1. The MS parameters of shape and compactness were defined as 0.1 and 0.5, respectively. We also tested 1-20 scale parameters of MS that were responsible for the termination of the merging process of the image objects. Segmentation was performed to delineate image objects based equally on two Fig. 6. Results of correlation matrix for features selected in this study; backscatter 150 kHz (bs150khz); backscatter 400 kHz (bs400khz); spectral skewness defined for central moments (20; SpSk_c_2); spectral moment m 0 (35; m 0 _3_5); spectral moment m 2 (35; m 2 _3_5); bathymetry (bath); Q-factor (35; Q_3_5); fractal dimension (35; Dfft_3_5); spectral width (35; sp_w_3_5); fractal dimension (20; Dfft_2); spectral skewness (35; SpSk_3_5). backscatter derivatives, namely, 150 kHz and 400 kHz.
Similar to other benthic habitat mapping studies, a few supervised classification approaches were tested to generate predictive outcomes based on ground-truth samples (Diesing et al., 2014;Hasan et al., 2012;Montereale Gavazzi et al., 2016). In this study, they included k-nearest neighbours (KNN), classification and regression trees (CART), RF, and support vector machines (SVM). We implemented these algorithms, which were available in the eCognition software. The KNN classified a specified object (query point) by a certain number (K) of known training samples that were located at the nearest neighbour around the query point. Euclidean distances (between the object and each instance) were calculated in the feature space to estimate the influence area of the neighbours. The KNN classification algorithm has been described in detail by Bremner et al. (2005). CART makes classification rules by recursively partitioning the data into increasingly homogenous groups. The algorithm created a decision tree that was associated with a system of questions and answers, thereby allowing the determination of the final classification (Breiman et al., 1984). RF was the machine learning method used for classification, regression, and other tasks, which consisted of constructing multiple decision trees that generated the class dominant or predicted average of individual trees (Breiman, 2001). SVM based on the machine learning technique used an algorithm that transformed datasets into a multidimensional feature space to find the appropriate boundary between them. Data points were called vectors, and the vectors that supported border selection were called support vectors. Machine learning models that use support vectors have been called SVM (Cortes and Vapnik, 1995).
We used validation ground-truth samples to create error matrices and calculate accuracy assessment statistics (Foody, 2002). These included the accuracy of the user and producer (Congalton, 1991;Story and Congalton, 1986), overall accuracy, and kappa index of agreement (Cohen, 1960).

Spectral parameters of bathymetry
The processed MBES bathymetry was presented in Fig. 1, and the co-registered backscatter for both frequencies was shown in Fig. 2.  Figs. 3 and 4 show the bathymetry and its eight derived spectral parameters for moving windows of 20 m × 20 m and 35 m × 35 m, respectively.
Visual insight for the created parameters demonstrated that several could match the geomorphologic features of bathymetry, especially a few spectral parameters (e.g. m 0 , m 2 , and spectral skewness defined for central moments) redrawn with specific features such as valleys and crests (Figs. 3 and 4). Moreover, the comparison of MBES backscatter (Fig. 2) showed a rough similarity to seabed areas of strong absorption or backscattering (visible as Dfft).

Ground-truth data processing
Six benthic habitat classes were determined, and they included VFS, S, SG_GS, B, R, and, A. The extensive identification of the benthic habitats with reference to MBES backscatter in this area was described by Janowski et al. (2018b).

Feature selection
The results of the Boruta feature selection are presented in Fig. 5. The algorithm performed 612 iterations and confirmed 11 features as important. The most important feature in this study was 400 kHz backscatter; this was followed by 150 kHz backscatter, Dfft (35), Dfft (20), spectral skewness defined for central moments (20), Q-factor (35), spectral skewness (35), bathymetry, spectral moment m 0 (35), spectral width (35), and spectral moment m 2 (35). The Boruta results indicated that certain spectral parameters were of greater significance than bathymetry, from which they were derived. This was especially visible in the Dfft parameter, which had a importance score approximately twice as great as that of bathymetry. It was also noteworthy that all the other extracted features (including the geomorphometric, statistical, and textural features) of the MBES bathymetry and backscatter were not considered important. Our initial results suggested the relevance of multi-frequency MBES; in other words, the first and second most important features were backscatter collected at different frequencies. The correlation analysis removed six highly correlated features, as shown in Fig. 6. The retained features were 400 kHz backscatter, bathymetry, spectral skewness defined for central moments (20), Q-factor (35), and spectral moment m 2 (35).
From the twenty scales of MS and four methods of supervised classification, the best classification performance was found with MS 8 and the SVM classifier. We adopted the following properties of the SVM classifier: radial-basis function kernel with C factor 100 and gamma 0.1. We created three predictive models using the following sets of features: (1) all relevant; (2) uncorrelated; and (3) only primary features (backscatter 400 kHz and bathymetry). The predictive benthic habitat maps generated using this approach are shown in Fig. 7.
The error matrix and accuracy assessment of the predictive habitat mapping method are presented in Table 2. Based on the validation subset of the ground-truth samples, the model with all relevant features confirmed a high performance that achieved a prediction accuracy of 86% and a Kappa index of agreement of 0.82. The second model, which considered only uncorrelated features, also achieved high accuracy; however, in comparison with the previous map, it misclassified one validation sample. The reference model without spectral features had an overall accuracy of 64% and a Kappa index of agreement of 0.55. Taking the accuracy of the user and producer into consideration, the two best-performing models were in reasonable agreement for specific classes, such as VFS, S, and R. We performed McNemar's chi-squared test for the statistical significance of differences in overall accuracy between the three models (Foody, 2004). The test result for differences between all relevant and uncorrelated features was 0.0. It means that the difference between these two models is statistically insignificant at the 5% level of significance. The McNemar's chi-squared test for differences between all relevant and only primary features models was 4.17 with p-value = .04, while the same test between uncorrelated and only primary features was 3.2 with p-value = .07. Mentioned results mean that there was a significant difference in the accuracy between all relevant and only primary features models and lack of significant difference between uncorrelated and only primary features models at the 5% level of significance.

Discussion
In this study, we introduced eight spectral features of a rough seafloor surface. The significance of the spectral features was evaluated and expressed using an importance score. We built and estimated the accuracy of three models of benthic habitat mapping, thereby demonstrating that the majority of introduced spectral features (i.e. seven out of eight) could improve the predictive power of supervised classifiers.
This study emphasised the importance of spectral parameters derived from bathymetry for predictive benthic habitat mapping based on multi-frequency MBES measurements. We did not observe significant differences in the bathymetry between both datasets (150 and 400 kHz). However, moderate differences existed in the backscatter of both frequencies, thus supporting the usefulness for a multi-frequency approach. We assumed that consistency in the bathymetry gathered with different frequencies was valid for sandy and gravelly sediments, as well as hard substrates. However, substantial depth differences in softer sediments could occur when significant acoustic penetration occurs (Schneider von Deimling et al., 2013). Furthermore, the beam resolution was linked with the frequency using one acoustic array with lower beam resolution when using lower frequencies. This could have affected the bathymetric results presented in this study. A visual comparison of the spectral parameters (particularly Dfft) indicated a high similarity between certain features of the multibeam backscatter datasets. Remarkably, the Boruta results showed that the spectral parameters of bathymetry, in general, had a greater significance than bathymetry itself. If the spectral parameters could match certain types of seabeds, they could be very useful for benthic habitat mapping at times when only MBES bathymetry is available. This study highlighted that the implementation of these spectral parameters could significantly improve supervised classification and benthic habitat mapping.
Other research of the Rowy Site demonstrated high applicability of the KNN and RF methods of classification (Janowski et al., 2018b). On the other hand, in this study, we obtained the best prediction performance with the SVM technique. Additional ground-truth samples were available for this study, thereby doubling the amount in relation to previous research (Janowski et al., 2018b). Increasing the number of samples provided a more realistic reference. Therefore, the predictions presented in this study could be considered as more robust.
The spectral analyses of surfaces, including Dfft, have already been applied in the analyses of seafloor data (Goff et al., 1999;Wilson et al., 2007) for morphological description. However, geomorphometry has been applied to the terrestrial environment very intensely (e.g. Sofia et al., 2016). Spectral analysis of the land surface has been used by geomorphologists; for example, Hutchinson and Gallant (2000) explored the usefulness of using numerical geomorphometric methods in terrain shape analysis. In the marine realm, studies using the autocovariance function of multibeam bathymetry successfully characterised the widths of morphological structures such as abyssal hills and continental slope canyons (Goff and Jordan, 1988;Goff, 2001). However, it was not determined as to what the size of the processing window of the spectral parameters should be for useful classification. This issue should be investigated further. The analysis of the spectral features, such as lidar data, from other sources of DEMs would be another topic worth exploring.
An MBES currently allows for the registration of absolute (calibrated) backscattering strength values; however, uncalibrated backscattering is still the most commonly used of these two measurement types. There is a need for calibrated systems with accurate hydroacoustic measurements so that data from different registrations can be compared. Because backscatter is dependent upon frequency, it can also be a disadvantage when compiling various datasets. In turn, the bathymetric spectral features represent absolute values, and therefore, they are reliable and easy to compare with measurements from other datasets. Spectral parameters are generally not dependent on the operating frequency of the MBES when the effects of sediment penetration can be excluded. However, taking such spectral features into account requires a high quality MBES bathymetry dataset and precise motion compensation. Any vessel motion artefacts can interfere, or 'leak', into the spectral features when they are not compensated. However, modern motion compensation systems work well correcting these errors. The dataset presented in this study was recorded on a vessel with a length of 8 m and width of 3 m at a sea state between 2 and 4. Despite unfavourable weather conditions for such a small vessel, the surveys yielded valuable results.
Because our study site was characterised by complex geomorphological features, we can ascertain that the presented method of predictive benthic habitat mapping could be especially valuable in other areas with diverse morphology (e.g. reefs). Considering a broader perspective, the spectral analysis of seafloor bathymetry could provide new insight into the analyses of DEMs of other sources, such as gravity models (Smith and Sandwell, 1997), which would allow for the exploration and interpretation of large scale complex geomorphological features, including volcanic structures, seamounts, or mid-ocean ridges.

Declaration of Competing Interest
None.