Seabed Mapping in Coastal Shallow Waters Using High Resolution Multispectral and Hyperspectral Imagery

Coastal ecosystems experience multiple anthropogenic and climate change pressures. To monitor the variability of the benthic habitats in shallow waters, the implementation of effective strategies is required to support coastal planning. In this context, high-resolution remote sensing data can be of fundamental importance to generate precise seabed maps in coastal shallow water areas. In this work, satellite and airborne multispectral and hyperspectral imagery were used to map benthic habitats in a complex ecosystem. In it, submerged green aquatic vegetation meadows have low density, are located at depths up to 20 m, and the sea surface is regularly affected by persistent local winds. A robust mapping methodology has been identified after a comprehensive analysis of different corrections, feature extraction, and classification approaches. In particular, atmospheric, sunglint, and water column corrections were tested. In addition, to increase the mapping accuracy, we assessed the use of derived information from rotation transforms, texture parameters, and abundance maps produced by linear unmixing algorithms. Finally, maximum likelihood (ML), spectral angle mapper (SAM), and support vector machine (SVM) classification algorithms were considered at the pixel and object levels. In summary, a complete processing methodology was implemented, and results demonstrate the better performance of SVM but the higher robustness of ML to the nature of information and the number of bands considered. Hyperspectral data increases the overall accuracy with respect to the multispectral bands (4.7% for ML and 9.5% for SVM) but the inclusion of additional features, in general, did not significantly improve the seabed map quality.


Introduction
Coastal ecosystems are essential because they support high levels of biodiversity and primary production, but their complexity and high spatial and temporal variability make their study particularly challenging.Seagrasses are extremely important marine angiosperms (flowering plants) with a worldwide distribution.Seagrass meadows are among the most productive ecosystems in the world, which help protect the shoreline from soil erosion, serve as a refuge area for other species, and absorb carbon from the atmosphere [1,2].Thus, seagrasses are essential, and their preservation in a sustainable manner needs the appropriate management tools.In this sense, satellite remote sensing is a cost-effective solution that has many advantages, compared to traditional techniques, like airborne photography with photo-interpretation or in-situ measurements (binomic maps from oceanographic ships).This way, satellite remote sensing is becoming a fundamental technology for the monitoring of benthic habitats (e.g., seagrass meadows) in shallow waters, as it provides periodic and synoptic data at different spatial scales and spectral resolutions [3].
Seafloor mapping using satellite remote sensing is a complex and challenging task, as optical bands have limited water penetration capability and the best channels to reach the seafloor (shorter wavelengths) suffer from higher atmospheric distortion.Hence, the signal recorded at the sensor level coming from the seabed is very low, even in clear waters [4,5].Towards the goal of mapping benthic habitats at high spatial resolution and achieving a reasonable accuracy, the use of hyperspectral (HS) imagery can be considered as an alternative to multispectral (MS) data.Unfortunately, high spatial hyperspectral sensors onboard satellites are not yet available and, in consequence, high-resolution data from airborne or drone HS sensors are the only options to collect HS data to map complex benthic habitats environments.
To map the seafloor, the use of high-resolution remote sensing is promising but requires the application of different geometric and radiometric corrections.Specifically, the removal of the atmospheric absorption and scattering and the sunglint effect over the sea surface are essential preprocessing steps.In addition, the water column disturbance can be corrected; however, it is a very complex issue in coastal areas due to the variability of the scattering and absorption in the water column, the bottom type, and the water depth [6].
Regarding the removal of the atmospheric effects, correction approaches can be basically grouped into physical radiative transfer models and empirical methods exclusively considering information obtained from the image scene itself [7].Many scene-based empirical approaches have been developed to remove atmospheric effects from multispectral and hyperspectral imaging data [8][9][10][11].Concerning the physical models, they are more advanced, complex and based on simulations of the conditions of the atmosphere from its physical-chemical characteristics and the day and time of acquisition of the image.At the present time, there are a number of model-based correction algorithms, for example MODerate resolution atmospheric TRANsmission (MODTRAN), Atmosphere CORrection Now (ACRON), Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH), High-accuracy Atmospheric Correction for Hyperspectral Data (HATCH), Atmospheric and Topographic CORrection (ATCOR), or Second Simulation of a Satellite Signal in the Solar Spectrum (6S) [12][13][14].Some of these algorithms include more advanced features, such as spectral smoothing, topographic correction, and adjacency effect correction.
On the other hand, the removal of sunglint is necessary for the reliable retrieval of bathymetry and seafloor mapping in shallow-water environments.Deglinting techniques have been developed for low-resolution open waters and also for high-resolution coastal applications [15].In general, algorithms use the near-infrared (NIR) channel to eliminate sunglint assuming that water reflectivity in the NIR band is negligible [16].This assumption is usually correct, except when turbidity is high, or the seabed reflectance is important, which can occur in very shallow areas [6].
Concerning the water column correction, Lyzenga [17] proposed the depth invariant index (DII), an image-based method to decrease the water column attenuation effect.This correction technique has been applied in previous works, due to its simplicity, with different degrees of success [18][19][20][21].On the other hand, in the last decades, some radiative transfer models have been proposed, but they are more complex, and the difficulty of accurately measuring some in-situ water parameters can limit their applicability [22][23][24][25].
Once preprocessing algorithms have been applied, classification techniques can be used to generate the seabed maps.Classification is one of the most active areas of research in the field of remotely sensed image processing.For example, the classification of hyperspectral imagery is a challenging task because of the imbalance among the high dimensionality of the data and the limited amount of available training samples, as well as the implicit spectral redundancy.For this reason, specific approaches have been developed, like random forests, support vector machines (SVMs), deep learning or logistic regressions [26].Unmixing techniques have also attracted the attention of the hyperspectral community.Unmixing algorithms separate the pixel spectra into a collection of constituent pure spectral signatures, named endmembers, and the corresponding set of fractional abundances, representing the percentage of each endmember that is present in the pixel [27].
Recent research to create seabed maps using remote sensing imagery has been mainly devoted to map coral reefs [28][29][30][31][32][33][34] or seagrass meadows [3,18,[35][36][37][38][39][40].Commonly, these studies address very shallow, clear and calm waters, and very dense vegetal species (i.e., Posidonia oceanica).As a continuation of our preliminary study [14], in this work, hyperspectral and multispectral imagery have been used to compare the benefits of each type of data to map the seafloor in a complex coastal area where submerged green aquatic vegetation meadows have low density, are relatively located at considerable depths (5 to 20 meters), where the sea surface is usually not completely calm due to persistent local surface winds, and, consequently, where very few bands reach an acceptable signal-to-noise ratio.Hence, a thorough analysis has been performed to obtain a robust methodology to produce accurate benthic habitat maps.To achieve this goal, different corrections, object-oriented, and pixel-based classification approaches have been considered, and diverse feature extraction strategies have also been tested.In summary, contributions are presented regarding the best correction techniques, feature extraction methods and classification approaches in such a challenging scenario.Moreover, a comparative assessment of the benefits of satellite multispectral and airborne hyperspectral imagery is included to map the seafloor in complex coastal zones.

Study Area
The Maspalomas Natural Reserve (Gran Canaria, Spain) is an important coastal-dune ecosystem covering approximately 4 km 2 and with a high touristic pressure, visited by more than 2 million people each year.The marine vegetation in the coastal fringe is basically composed of seagrass beds.The most abundant seagrass species is Cymodocea nodosa but, more recently, the Caulerpa prolifera green algae has also become dominant in this area.Figure 1 shows the geographic location of Maspalomas.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 21 specific approaches have been developed, like random forests, support vector machines (SVMs), deep learning or logistic regressions [26].Unmixing techniques have also attracted the attention of the hyperspectral community.Unmixing algorithms separate the pixel spectra into a collection of constituent pure spectral signatures, named endmembers, and the corresponding set of fractional abundances, representing the percentage of each endmember that is present in the pixel [27].
Recent research to create seabed maps using remote sensing imagery has been mainly devoted to map coral reefs [28][29][30][31][32][33][34] or seagrass meadows [3,18,[35][36][37][38][39][40].Commonly, these studies address very shallow, clear and calm waters, and very dense vegetal species (i.e., Posidonia oceanica).As a continuation of our preliminary study [14], in this work, hyperspectral and multispectral imagery have been used to compare the benefits of each type of data to map the seafloor in a complex coastal area where submerged green aquatic vegetation meadows have low density, are relatively located at considerable depths (5 to 20 meters), where the sea surface is usually not completely calm due to persistent local surface winds, and, consequently, where very few bands reach an acceptable signalto-noise ratio.Hence, a thorough analysis has been performed to obtain a robust methodology to produce accurate benthic habitat maps.To achieve this goal, different corrections, object-oriented, and pixel-based classification approaches have been considered, and diverse feature extraction strategies have also been tested.In summary, contributions are presented regarding the best correction techniques, feature extraction methods and classification approaches in such a challenging scenario.Moreover, a comparative assessment of the benefits of satellite multispectral and airborne hyperspectral imagery is included to map the seafloor in complex coastal zones.

Study Area
The Maspalomas Natural Reserve (Gran Canaria, Spain) is an important coastal-dune ecosystem covering approximately 4 km 2 and with a high touristic pressure, visited by more than 2 million people each year.The marine vegetation in the coastal fringe is basically composed of seagrass beds.The most abundant seagrass species is Cymodocea nodosa but, more recently, the Caulerpa prolifera green algae has also become dominant in this area.Figure 1

Multisensor Remotely Sensed Data
A flight campaign was performed on June 2, 2017 and data were collected at 2.5 m resolution with the Airborne Hyperspectral Scanner (AHS).The AHS sensor (developed by ArgonST, USA) is operated by the Spanish Aerospace Institute (INTA) onboard the CASA 212-200 Paternina.AHS incorporates an 80-band imaging radiometer covering the range from 0.43 to 12.8 µm.In our study, only the first 20 channels were selected, covering the visible and near-infrared (NIR) spectrum from 0.434 to 1.015 µm with 12-bits of radiometric resolution [41].Additionally, a WorldView-2 (WV-2) image collected on January 17, 2013 was considered for this study.Its sensor has a radiometric resolution of 11 bits and a spatial resolution, at the nadir, of 1.8 m for the 8 multispectral bands (0.40-1.04 µm).The panchromatic wideband was not used because this channel only provides information about the seabed in the very first meters of depth and, in consequence, pansharpening algorithms are not effective.
Table 1 includes the spectral characteristics of both sensors and Figure 2a,b show the Worldview and AHS imagery for the area of interest processed in this work, respectively.

In-Situ Measurements
In-situ data were acquired simultaneously to the AHS campaign.A total of 6 transects were performed, measuring the bathymetry using an ecosounder Reson Navisound 110 and recording images of the seafloor using two different video cameras (Neptune and Go Pro Hero 3+).Precise geographic and temporal information was provided by a differential GPS receiver model Trimble DSM132.Ten additional sites visited during 2015 were also used in the analysis, providing bathymetry and video records from the Go Pro camera.The variation of the sea level due to tides and waves was obtained from a nearby calibrated tide gauge.
Figure 2c presents the real ship transects during the 2017 campaign, as well as the sites monitored in 2015 (marked by yellow dots).Isobaths are also included, and the sites and equidistant transects are perpendicular to the shore to get the maximum information of the area at different depths up to 20 m.
Finally, to assess the accuracy of the benthic maps derived from HS and MS imagery, apart from the accurate information from the in-situ transect and sampling sites, a reference map providing global information for the whole area was desirable.Unfortunately, very limited cartography is available and the map considered (Figure 2d) is the reference map available of the Maspalomas coast, which is part of a Spanish coastline eco-cartographical study from a series of maritime engineering and marine ecology studies structured in a GIS [42].
This information was used just as a coarse reference, but the quantitative validation only considered the precise information measured during the campaigns of 2015 and 2017.
The data used in this study correspond to different years (2013 to 2017).However, the seabed of Maspalomas is in a fairly stable area with the exception of Punta de la Bajeta (marked in Figure 2a), which is the zone with the greatest topographic and sedimentary variability due to storms from the southwest.
The seafloor images recorded during the field campaign (see Figure 3 for examples) show the complexity to discriminate habitat classes because they are usually mixed.Especially, submerged green aquatic vegetation meadows grow on sandy bottoms, and they have low density and mainly live between 5 and 25 m depth.In addition, rocks are partially covered by algae making automated classifications more challenging.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 21 This information was used just as a coarse reference, but the quantitative validation only considered the precise information measured during the campaigns of 2015 and 2017.
The data used in this study correspond to different years (2013 to 2017).However, the seabed of Maspalomas is in a fairly stable area with the exception of Punta de la Bajeta (marked in Figure 2a), which is the zone with the greatest topographic and sedimentary variability due to storms from the southwest.
The seafloor images recorded during the field campaign (see Figure 3 for examples) show the complexity to discriminate habitat classes because they are usually mixed.Especially, submerged green aquatic vegetation meadows grow on sandy bottoms, and they have low density and mainly live between 5 and 25 m depth.In addition, rocks are partially covered by algae making automated classifications more challenging.

Mapping Methodology
The overall processing protocol to generate the seabed maps is presented in Figure 4. Different inputs to the classification algorithm were analyzed to select the best methodology.Following, there is a detailed description of the different steps involved.

•
Multispectral Satellite Imagery DigitalGlobe owns and operates a world-class constellation of high resolution, high accuracy Earth imaging satellites.The acquired WV-2 image is the level 2 ortho-ready product that has the geometric correction implemented and a horizontal accuracy specification of 5 meters or better [43].
In coastal areas, radiometric and atmospheric corrections have proven to be a crucial step in the processing of high-resolution satellite images due to the low signal-to-noise ratio at sensor level.A comparative evaluation of advanced atmospheric models (FLAASH, ATCOR, and 6S) in the coastal area of Maspalomas using high-resolution WorldView-2 data showed that the 6S algorithm achieved the highest accuracy when the corrected reflectance was compared to field spectroradiometer data (RMSE of 0.0271 and bias of −0.0217) [14].Hence, we have applied the 6S correction model but adapted to the WorldView-2 spectral response and to the particular scene geometry and date of the image.6S is a radiative transfer model that generates the constants ( , ) to estimate the surface (BOA: Bottom Of Atmosphere) reflectance by [14,44,45]: corrected to consider the adjacency effect by: where is the radiance measured by the sensor (TOA: Top Of Atmosphere), is the initially corrected surface reflectivity by the 6S model, is the surface reflectivity taking into account the

Mapping Methodology
The overall processing protocol to generate the seabed maps is presented in Figure 4. Different inputs to the classification algorithm were analyzed to select the best methodology.Following, there is a detailed description of the different steps involved.Moreover, specular reflection of solar radiation is a serious disturbance for water quality, bathymetry, and benthic mapping in shallow-water environments.We applied the method suggested by References [15,16] to remove sunglint.Therefore, regions of the image having sunglint were selected, preferably areas of deep water.For each visible channel, all the pixels from these regions were included in a linear regression between the NIR against each visible band.Therefore, the reflectance of each pixel in the visible band i ( , ) can be deglinted ( , ) using the following

Multisensor Imagery Corrections
• Multispectral Satellite Imagery DigitalGlobe owns and operates a world-class constellation of high resolution, high accuracy Earth imaging satellites.The acquired WV-2 image is the level 2 ortho-ready product that has the geometric correction implemented and a horizontal accuracy specification of 5 meters or better [43].
In coastal areas, radiometric and atmospheric corrections have proven to be a crucial step in the processing of high-resolution satellite images due to the low signal-to-noise ratio at sensor level.A comparative evaluation of advanced atmospheric models (FLAASH, ATCOR, and 6S) in the coastal area of Maspalomas using high-resolution WorldView-2 data showed that the 6S algorithm achieved the highest accuracy when the corrected reflectance was compared to field spectroradiometer data (RMSE of 0.0271 and bias of −0.0217) [14].Hence, we have applied the 6S correction model but adapted to the WorldView-2 spectral response and to the particular scene geometry and date of the image.6S is a radiative transfer model that generates the constants (x a , x b and x c ) to estimate the surface (BOA: Bottom Of Atmosphere) reflectance by [14,44,45]: corrected to consider the adjacency effect by: where L TOA is the radiance measured by the sensor (TOA: Top Of Atmosphere), ρ BOA is the initially corrected surface reflectivity by the 6S model, ρ BOA is the surface reflectivity taking into account the adjacency effect, τ o di f and τ o dir are the diffuse and direct transmittances and ρ is the average reflectivity contribution from the pixel background [44].Moreover, specular reflection of solar radiation is a serious disturbance for water quality, bathymetry, and benthic mapping in shallow-water environments.We applied the method suggested by References [15,16] to remove sunglint.Therefore, regions of the image having sunglint were selected, preferably areas of deep water.For each visible channel, all the pixels from these regions were included in a linear regression between the NIR against each visible band.Therefore, the reflectance of each pixel in the visible band i (ρ BOA,i ) can be deglinted ρ DG BOA,i using the following equation: where b i is the slope of the regression line for band i, ρ BOA,N IR is the reflectance of the NIR channel and MIN NIR corresponds to the minimum reflectance value of the NIR image.Improvements in the glint removal algorithm were performed [46] because, for WorldView-2, all the sensor bands are not recording the energy at the same precise time.In addition, the deglinting process, in the presence of considerable waves, alters the spectral content of the image.To overcome this inconvenience, a histogram matching was applied to statistically equalize each channel after and before deglinting.Finally, to remove the foam caused by the waves (whitecaps), pixels achieving reflectance values above a threshold were replaced by interpolation [46].
Water column correction is a very complex matter due to the variability of the bottom type, water depth, and water attenuation (scattering and absorption in the water column).Lyzenga proposed a simpler depth invariant index [17] but, in the last decades, different water models have been developed.In this work, a radiative transfer model for coastal waters was implemented [14].The assumption is that the water leaving radiance is not only caused by the water column (IOPs: water Inherent Optical Properties), but it is also affected by the seafloor albedo and its corresponding bathymetry.In particular, the Lee et al. [47] exponential expression was used, that is an improved version of the formulation proposed by Reference [22].The model can be expressed by: where r m rs (λ) is the modeled reflectivity, r rs,∞ (λ) corresponds to the reflectivity below the sea surface for deep waters generated by the IOPs, ρ alb (λ) represents the seafloor reflectivity (albedo), z stands for the bathymetry, k d is the water diffuse attenuation coefficient, µ sw s and µ sw v allow the corrections in the sun and sensor trajectories for the downwards and upwards directions, and, finally, D c u and D b u are light diffusion ascending factors due to the water column and the bottom reflectivity, respectively.
As proposed in Reference [48], the Fully Constrained Linear Unmixing method (FCLU) was used to model the seabed albedo.The albedo is obtained as the sum of the products of the abundances of the p pure benthic elements (ab i ) by their albedo for each wavelength (em i (λ)): Due to the limitations of the sensitivity of the multispectral sensors, we decided to use the three most significant and separable spectra of the seabed coverage in the area (sand, rocks and green vegetation).The radiative transfer modeling was adapted to multispectral sensors by integrating the result of the monochromatic model for the bandwidths of the 6 first channels of WV-2.

• Airborne Hyperspectral Imagery
Regarding the AHS geometric accuracy, the airborne inertial system achieves a final angular precision of 0.008 • for the roll and the pitch, and 0.015 • for the true heading; and a 12 channels GPS receiver provides trajectory location with accuracies of 5 to 10 cm [41].
For the hyperspectral data, atmospheric and illumination corrections were performed by INTA using the ATCOR4 model [41,49].ATCOR4 is the ATCOR specific version for hyperspectral airborne data, and it is based on the MODTRAN-5 radiative transfer code.
The multispectral channels of the AHS sensor are somewhat narrower than the WV2 channels, providing, in the visible range, 12 multispectral channels instead of the 6 of the WV2.The greater number of channels implies greater spectral information and, therefore, an assumed greater sensitivity in the classification of the benthic classes.In any case, the methodology carried out for the elimination of the specular solar brightness, and the radiative transfer modeling is identical to that applied to WV-2.
Finally, land and deep water masks were applied to the corrected bands after both sensors have been homogenized in terms of resolution and have been spatially aligned using a large database of singular and well-distributed ground control points.Seafloor is difficult to properly monitor over 20 m depth in this area using satellite or airborne remote sensing imagery; therefore, a mask was applied at the 20 m isobath.

Feature Extraction
In addition to the spectral channels provided by each sensor, we also studied the inclusion of additional information to check the possible increase in classification performance.Therefore, we considered adding, to the multispectral or hyperspectral bands, components after Principal Component Analysis (PCA), Independent Component Analysis (ICA) and Minimum Noise Fraction (MNF) transforms; textural features to enrich the spatial information, and abundance maps extracted from linear unmixing techniques.Next, the feature extraction techniques are explained in more detail.

• Image Transforms
In hyperspectral imagery, the high number of narrow bands requires an increase in the number of training pixels to maintain a minimum statistical confidence and functionality for classification.This problem, known as the Hughes' phenomenon [50] or the curse of dimensionality, can be addressed by overcoming data redundancy.Several transforms have been proposed in the last decades to extract reliable information, reducing redundancy and noise.Traditional feature-extraction techniques, mainly applied to reduce the dimensionality of hyperspectral data, are PCA, ICA, and MNF [51].
PCA uses an orthogonal transformation to convert a set of correlated bands into a new set of uncorrelated components [52].PCA is frequently applied for dimensionality reduction because it retains most of the information of the original data in the first principal components.In consequence, computation times decrease, and the Hughes' phenomenon is avoided, and, accordingly, the reduction of the considered components allows obtaining more precise thematic maps [53].
ICA decomposes the set of image bands into a linear combination of independent source signals [54].ICA not only decorrelates second-order statistics, like PCA does but also decreases higher-order dependencies, generating a new set of components as independent as possible.It is an alternative approach to PCA for dimensionality reduction.
MNF Rotation is a linear transformation to segregate noise in the data and to reduce the computational requirements for subsequent processing [55].This MNF transformation orders the new components so that they maximize the signal-to-noise ratio, rather than the information content [56].MNF, apart from reducing the dimensionality, can be used to remove noise from data by performing a forward transform, by determining which bands contain the coherent information, by examining the images and eigenvalues, and then applying the inverse transform using only the appropriate bands, or adding as well the filtered or smoothed noisy bands.

• Texture Maps
We analyzed the inclusion of spatial information into the classification to improve the separability between the classes.There is a wide variety of texture measurements and, in this experiment, the parameters used were derived from the Gray-Level-Co-occurrence Matrix (GLCM) [57].The main idea of the GLCM measurements is that the texture information contained in an image is based on the adjacency relationship between gray levels in the image.The relationship of the occurrence frequencies between pixel pairs can be calculated reliably for specific directions and distances between them.
After a preliminary analysis, we observed that most of the texture maps, derived from the co-occurrence matrix, were very similar (variance, entropy, dissimilarity, etc.).Therefore, to avoid the inclusion of redundant information, only the mean and variance parameters were finally used.Instead of applying the GLCM to each original band and, consequently, to have a large dataset of redundant texture information, the GLCM was only applied to the first component of the PCA and MNF transforms.

• Abundance Maps
In general, many pixels in the image represent a mixture of spectral signatures from different classes (e.g., seagrass and sand).In this context, mixing models estimate the contribution (abundance) of each endmember (pure class) to the total reflectance of each pixel.Linear mixing models provide adequate results in a large number of applications [27], and while non-linear models can be more precise, they require detailed information about the geometry and physical properties of the objects, which is not usually available and, thus, hinders its usefulness.
There are different strategies for the selection of pure pixels.In this work, as supervised classification techniques were applied, the training regions selected for each class were used to get the spectral signatures of each endmember.
After the application of unmixing techniques, an abundance map for each class to be discriminated was generated.Each map represents the percentage of contribution of a class to the total reflectance value of the pixel.In consequence, maps have values between 0 and 1 and, if a pixel is a mixture of different classes, the abundance of each class is obtained.
ML classification is one of the most common supervised techniques used with remotely sensed data [51].ML considers that each class can be modeled as a normal distribution, allowing to describe that class by a Gaussian probability function, from its vector of means and covariance matrix.ML assigns the pixel to the class that maximizes the probability function.
The SAM classifier [52] compares the similarity between two spectra from their angular deviation, assuming that they form two vectors in an n-dimensional space (n being the number of bands).This algorithm measures the similarity as a function of the angle both vectors form in such a space.
SVM [60] is a machine learning algorithm to discriminate two classes by fitting an optimal hyperplane to separate the training samples of each class.The samples closest to the decision boundary are the so-called support vectors.SVM has been efficiently applied to classify both linear and nonlinearly separable classes applying a kernel function into a higher dimensional space, whose new data distribution allows better fitting of a linear hyperplane.Although deep learning approaches are becoming popular [61], they require large training datasets, and that is a great inconvenience in many operational applications.In addition, a recent assessment comparing advanced classification methods (SVMs, random forests, neural networks, deep convolutional neural networks, logistic regression-based techniques, and sparse representation-based classifiers) demonstrated that SVM is widely used because of its accuracy, stability, automation, and simplicity [26].After a detailed review of SVM literature [18,[62][63][64] and many tests conducted using the SVM algorithm with high-resolution imagery [65,66], the Gaussian radial basis function kernel K x i , x j was selected for the SVM classification and their parameters were properly adjusted: For the segmentation and merging steps in the object based classification (OBIA), after testing different combinations, using diverse features and number of bands, AHS channels 1 to 3 and WV-2 channel 1 were used.During the segmentation process, the image is divided into homogeneous regions according to several parameters (band weights, scale, color, shape, texture, etc.) defined by the operator, with the objective of creating the suitable object borders.We tested a multiresolution segmentation approach [67] and an algorithm based on watershed segmentation and merging stages [68].An over-segmented image was preferred, and as soon as a suitable segmentation was attained, the classifier was applied.
Finally, the classification accuracy for each possible combination of input data was estimated using independent test regions of interest (ROIs) located in the image and computing the kappa coefficient, the confusion matrix, and its derived measures.

Results and Discussion
After the correction of each dataset (see Figure 4), three supervised classifiers were applied to different combinations of input data.All the analysis was performed for HS and MS imagery (AHS and WV-2, respectively) at the complex area of Maspalomas.
As seafloor reflectivity is very weak, and following the steps of Figure 4, precise preprocessing algorithms were applied to correct limitations in the sensor calibration, solar illumination geometry, viewing effects, as well as the atmospheric, sunglint, and water column disturbances.In this sense, geometric, radiometric and atmospheric corrections were performed.As specified, 6S was selected to model the atmosphere and to remove the absorption and scattering effects in the multispectral image [14], and ATCOR4 for the hyperspectral data [49].This selection took into account the results of a previous validation campaign comparing real sea surface reflectance recorded by a field spectroradiometer (ADS Fieldspec 3) and the reflectance estimated for WV-2 data using different models.Next, deglinting algorithms were applied to eliminate the solar glint and whitecaps over both datasets.Finally, the seafloor albedo was generated applying the radiative transfer model described in Section 2.4.1.As shown in Figure 5, for a small area, the improvement is considerable, especially for the AHS data, as some areas were severely affected by de-sunglint.
After the correction of each dataset (see Figure 4), three supervised classifiers were applied to different combinations of input data.All the analysis was performed for HS and MS imagery (AHS and WV-2, respectively) at the complex area of Maspalomas.
As seafloor reflectivity is very weak, and following the steps of Figure 4, precise preprocessing algorithms were applied to correct limitations in the sensor calibration, solar illumination geometry, viewing effects, as well as the atmospheric, sunglint, and water column disturbances.In this sense, geometric, radiometric and atmospheric corrections were performed.As specified, 6S was selected to model the atmosphere and to remove the absorption and scattering effects in the multispectral image [14], and ATCOR4 for the hyperspectral data [49].This selection took into account the results of a previous validation campaign comparing real sea surface reflectance recorded by a field spectroradiometer (ADS Fieldspec 3) and the reflectance estimated for WV-2 data using different models.Next, deglinting algorithms were applied to eliminate the solar glint and whitecaps over both datasets.Finally, the seafloor albedo was generated applying the radiative transfer model described in Section 2.4.1.As shown in Figure 5, for a small area, the improvement is considerable, especially for the AHS data, as some areas were severely affected by de-sunglint.A preliminary analysis was performed to find the most suitable corrected imagery to address the classification problem.Specifically, images obtained after the different pre-processing steps were assessed to identify the more reliable data source for the mapping production.The following thematic classes were considered: sand (yellow), rocks (brown) and Cymodocea or Caulerpa (green).Using the information from the ship transects and sampling sites (Figure 2c), sets of training and validation regions were generated including regions of each class at five-meter step depths from 5 to 20 m.Approximately, 3000 and 6000 pixels per class were selected for the training and test ROIs, respectively.The class pair separability (Jeffries-Matusita distance [52]) in the bands ranges from 1.218 to 1.693 for WV-2 and between 1.802 and 1.985 for AHS.These values corroborate a better discrimination capability of HS data as more spectral richness is available.
Table 2 presents the results of applying the three supervised classifiers to the data after the atmospheric, sunglint, and water column corrections.The same independent training and validation regions of interest were used in all the experiments.As expected, the airborne hyperspectral imagery allows a better classification than the satellite multispectral data (92.01%with respect to 88.66%).It can be appreciated that the best overall accuracy was achieved after the deglinting step.The water column removal did not improve the seafloor mapping, even after applying a complex radiative A preliminary analysis was performed to find the most suitable corrected imagery to address the classification problem.Specifically, images obtained after the different pre-processing steps were assessed to identify the more reliable data source for the mapping production.The following thematic classes were considered: sand (yellow), rocks (brown) and Cymodocea or Caulerpa (green).Using the information from the ship transects and sampling sites (Figure 2c), sets of training and validation regions were generated including regions of each class at five-meter step depths from 5 to 20 m.Approximately, 3000 and 6000 pixels per class were selected for the training and test ROIs, respectively.The class pair separability (Jeffries-Matusita distance [52]) in the bands ranges from 1.218 to 1.693 for WV-2 and between 1.802 and 1.985 for AHS.These values corroborate a better discrimination capability of HS data as more spectral richness is available.
Table 2 presents the results of applying the three supervised classifiers to the data after the atmospheric, sunglint, and water column corrections.The same independent training and validation regions of interest were used in all the experiments.As expected, the airborne hyperspectral imagery allows a better classification than the satellite multispectral data (92.01%with respect to 88.66%).It can be appreciated that the best overall accuracy was achieved after the deglinting step.The water column removal did not improve the seafloor mapping, even after applying a complex radiative model.Even providing adjusted water IOPs and bathymetry values, the modeling of the background albedo by linear mixing of benthic classes in this complex area does not seem adequate for the subsequent classification.The very low reflectivity of the coastal bottom, which usually contributes less than 1% of the radiation observed by the sensor, produces errors in the adjustment of the abundances of the modeled pure benthic elements.Clearly, the model considered has to be further improved.As indicated, the water column modelling in coastal areas is complex and depends on the water quality parameters, as well as the bathymetry and the type of seabed.For this reason, this preprocessing is not always considered and some studies demonstrate that better results are not always achieved [18,20,40].Finally, regarding the classification algorithm, SVM is the most appropriate approach for AHS but Maximum Likelihood works better with WV-2.
Figure 6 shows examples of seafloor maps generated for the AHS sensor, with SVM, and for the WV-2 data using the ML classifier.Comparing the results with the reference benthic map (Figure 2d) and the available video records from the ship transects, higher accuracy can be noted for AHS and using the imagery after the atmospheric and sunglint correction (middle row).Excessive amount of submerged vegetation is identified for WV-2 and some rocks incorrectly appear on the right side when these pixels should be labeled as vegetation.Figure 6 shows examples of seafloor maps generated for the AHS sensor, with SVM, and for the WV-2 data using the ML classifier.Comparing the results with the reference benthic map (Figure 2d) and the available video records from the ship transects, higher accuracy can be noted for AHS and using the imagery after the atmospheric and sunglint correction (middle row).Excessive amount of submerged vegetation is identified for WV-2 and some rocks incorrectly appear on the right side when these pixels should be labeled as vegetation.To improve the previous seabed cartography, a detailed feature extraction and classification assessment was only performed using the preprocessed data after the atmospheric and sunglint correction stages.
As stated in Section 2.4.2, to improve the benthic maps, additional information was obtained using feature extraction techniques.In particular, PCA, ICA and MNF were applied to the corrected spectral bands.In the analysis, the classifier performance was assessed including the complete new set of components after these transforms and, in addition, the best components were also tested discarding noisy bands.Figure 7 shows the first bands of each transform, as well as the original spectral bands as a reference (the remaining bands were not displayed as they are too noisy).Regarding the spectral channels, we can appreciate that only shorter wavelengths (first bands) can reach the seafloor and, in consequence, even dealing with hyperspectral data only a few channels are really valuable to map benthic habitats up to a depth of 20 m.On the other hand, PCA, ICA and MNF provide useful information in the first four components.The true color image and false color To improve the previous seabed cartography, a detailed feature extraction and classification assessment was only performed using the preprocessed data after the atmospheric and sunglint correction stages.
As stated in Section 2.4.2, to improve the benthic maps, additional information was obtained using feature extraction techniques.In particular, PCA, ICA and MNF were applied to the corrected spectral bands.In the analysis, the classifier performance was assessed including the complete new set of components after these transforms and, in addition, the best components were also tested discarding noisy bands.Figure 7 shows the first bands of each transform, as well as the original spectral bands as a reference (the remaining bands were not displayed as they are too noisy).Regarding the spectral channels, we can appreciate that only shorter wavelengths (first bands) can reach the seafloor and, in consequence, even dealing with hyperspectral data only a few channels are really valuable to map benthic habitats up to a depth of 20 m.On the other hand, PCA, ICA and MNF provide useful information in the first four components.The true color image and false color composites using the first three components are also included, and it is possible to check the worse behavior of ICA and the noise removal effect of MNF.
assessment was only performed using the preprocessed data after the atmospheric and sunglint correction stages.
As stated in Section 2.4.2, to improve the benthic maps, additional information was obtained using feature extraction techniques.In particular, PCA, ICA and MNF were applied to the corrected spectral bands.In the analysis, the classifier performance was assessed including the complete new set of components after these transforms and, in addition, the best components were also tested discarding noisy bands.Figure 7 shows the first bands of each transform, as well as the original spectral bands as a reference (the remaining bands were not displayed as they are too noisy).Regarding the spectral channels, we can appreciate that only shorter wavelengths (first bands) can reach the seafloor and, in consequence, even dealing with hyperspectral data only a few channels are really valuable to map benthic habitats up to a depth of 20 m.On the other hand, PCA, ICA and MNF provide useful information in the first four components.The true color image and false color composites using the first three components are also included, and it is possible to check the worse behavior of ICA and the noise removal effect of MNF.Additional textural parameters and abundance maps after unmixing were also inputted to the classifiers as auxiliary information.
Table 3 summarizes the AHS and WV-2 accuracy results of each classifier for the following input combinations: • Spectral bands after atmospheric and sunglint corrections.Additional textural parameters and abundance maps after unmixing were also inputted to the classifiers as auxiliary information.
Table 3 summarizes the AHS and WV-2 accuracy results of each classifier for the following input combinations:

•
Spectral bands after atmospheric and sunglint corrections.

•
Components after the application of three-dimensionality reduction techniques (PCA, ICA, and MNF).The complete dataset and a reduced number of bands or components were both tested.

•
Abundance maps of each class after the application of linear unmixing techniques.

•
Texture information (mean and variance) extracted from the first PCA/MNF component.Pixel-based classification was applied to the previous options and, finally, object-based classification was applied to the spectral bands.
With respect to the sensors, we can appreciate that AHS provides better accuracy than WV-2, as expected, mainly due to the availability of additional bands and a better radiometric resolution.Specifically, a major improvement is attained for SVM (mean accuracy increase of 9.5%) than for ML (4.7% average increase).
Concerning the classification algorithms, SAM did not work properly because, even being more insensitive to variations of the bathymetry, classes are spectrally overlapped, and only very few bands are useful due to the water column attenuation.SVM is the algorithm achieving the best accuracy, but the simpler and faster ML demonstrates good performance and, in many cases, better than SVM (average results in Table 3 confirm it).Actually, Figure 8 presents the comparative performance of both classifiers, and it can be appreciated that ML is more robust, providing more stable results regardless of the input information used or the number of bands considered.Specifically, the standard deviation (averaged for AHS and WV-2) of the overall accuracy for the different combinations is 2.9% for ML and 6.4% for SVM.Pixel-based classification was applied to the previous options and, finally, object-based classification was applied to the spectral bands.
With respect to the sensors, we can appreciate that AHS provides better accuracy than WV-2, as expected, mainly due to the availability of additional bands and a better radiometric resolution.Specifically, a major improvement is attained for SVM (mean accuracy increase of 9.5%) than for ML (4.7% average increase).
Concerning the classification algorithms, SAM did not work properly because, even being more insensitive to variations of the bathymetry, classes are spectrally overlapped, and only very few bands are useful due to the water column attenuation.SVM is the algorithm achieving the best accuracy, but the simpler and faster ML demonstrates good performance and, in many cases, better than SVM (average results in Table 3 confirm it).Actually, Figure 8 presents the comparative performance of both classifiers, and it can be appreciated that ML is more robust, providing more stable results regardless of the input information used or the number of bands considered.Specifically, the standard deviation (averaged for AHS and WV-2) of the overall accuracy for the different combinations is 2.9% for ML and 6.4% for SVM.PCA and MNF perform much better than ICA, but the improvement is, in general, negligible with respect to the original bands.Also, the reduction of the number of bands/components to avoid the Hughes' phenomenon is basically not increasing the classification accuracy except for ML and the hyperspectral data.The number of training pixels for each class is high enough (3000), and that could be a possible explanation.PCA and MNF perform much better than ICA, but the improvement is, in general, negligible with respect to the original bands.Also, the reduction of the number of bands/components to avoid the Hughes' phenomenon is basically not increasing the classification accuracy except for ML and the hyperspectral data.The number of training pixels for each class is high enough (3000), and that could be a possible explanation.
The application of unmixing techniques before the classification did not improve the accuracy due to the small number of bands actually available.It can be appreciated the degraded performance of SVM when only the three abundances are considered in the classification scheme.
Finally, texture information is a feature that can be included in the final methodology as precision values in some circumstances increase the performance.Specifically, the improvement is more evident for SVM and using the texture information provided by the first component of PCA.For ML, texture generally does not provide a better accuracy of the benthic map.
It is important to highlight that results obtained by the object-based classification techniques (OBIA) are not always the best.Basically, OBIA only provides superior performance than pixel-based techniques for the SVM algorithm.However, results are quite dependent on the type of segmentation considered.
In general, the overall accuracies for ML and SVM are high as few classes are considered and the validation pixels chosen to numerically assess accuracy were selected in clear and central locations of each seabed type.In any case, the relative results between the different classifiers and input combinations displayed in Table 3 are reliable.Figure 9 includes an example of the AHS and WV-2 segmentation for a specific area.AHS provides more detailed information and, in consequence, the number of objects increases.
For ML, texture generally does not provide a better accuracy of the benthic map.
It is important to highlight that results obtained by the object-based classification techniques (OBIA) are not always the best.Basically, OBIA only provides superior performance than pixel-based techniques for the SVM algorithm.However, results are quite dependent on the type of segmentation considered.
In general, the overall accuracies for ML and SVM are high as few classes are considered and the validation pixels chosen to numerically assess accuracy were selected in clear and central locations of each seabed type.In any case, the relative results between the different classifiers and input combinations displayed in Table 3 are reliable.Figure 9 includes an example of the AHS and WV-2 segmentation for a specific area.AHS provides more detailed information and, in consequence, the number of objects increases.Figure 10 compares the best pixel-based seafloor maps generated by ML and SVM for the AHS image.A majority filter of 5 × 5 window size was applied to remove the salt and pepper effect.Both maps are very accurate, but ML overestimates vegetation (green) in some specific areas, while SVM the rocks (brown) in others.Finally, Figure 11 shows the best maps for each sensor obtained using the object-based classification with the SVM algorithm.Results are similar and, in general, match the available eco-cartographic map included in Figure 2d, except for the western side of the area.In any case, as indicated, this map was just considered a coarse reference and really ship transects T1 and T2 in Figure 2c demonstrate the existence of vegetation meadows in that area, in agreement with AHS and WV-2 maps.It is also important to highlight that a vulnerable and complex ecosystem was studied where the density of submerged green aquatic vegetation beds is quite low and, therefore, there is a considerable mixture of sand and plant contributions in each pixel of the image.These methodologies will be shortly applied to generate precise benthic maps of natural protected ecosystems in other vulnerable coastal ecosystems.In addition, these will be applied to hyperspectral imagery recorded from drone platforms with the goal of discriminating between the different vegetation species.Figure 10 compares the best pixel-based seafloor maps generated by ML and SVM for the AHS image.A majority filter of 5 × 5 window size was applied to remove the salt and pepper effect.Both maps are very accurate, but ML overestimates vegetation (green) in some specific areas, while SVM the rocks (brown) in others.Finally, Figure 11 shows the best maps for each sensor obtained using the object-based classification with the SVM algorithm.Results are similar and, in general, match the available eco-cartographic map included in Figure 2d, except for the western side of the area.In any case, as indicated, this map was just considered a coarse reference and really ship transects T1 and T2 in Figure 2c demonstrate the existence of vegetation meadows in that area, in agreement with AHS and WV-2 maps.It is also important to highlight that a vulnerable and complex ecosystem was studied where the density of submerged green aquatic vegetation beds is quite low and, therefore, there is a considerable mixture of sand and plant contributions in each pixel of the image.

Conclusions
A comprehensive analysis was performed to identify the best input dataset and to obtain a robust classification methodology to generate accurate benthic habitat maps.The assessment considered pixel-based and object-oriented classification methods in shallow waters using hyperspectral and multispectral data.
A vulnerable and complex coastal ecosystem was selected where the submerged green aquatic vegetation meadows to be classified are located at depths between 5 and 20 meters and have low These methodologies will be shortly applied to generate precise benthic maps of natural protected ecosystems in other vulnerable coastal ecosystems.In addition, these will be applied to hyperspectral imagery recorded from drone platforms with the goal of discriminating between the different vegetation species.

Conclusions
A comprehensive analysis was performed to identify the best input dataset and to obtain a robust classification methodology to generate accurate benthic habitat maps.The assessment considered pixel-based and object-oriented classification methods in shallow waters using hyperspectral and multispectral data.
A vulnerable and complex coastal ecosystem was selected where the submerged green aquatic vegetation meadows to be classified are located at depths between 5 and 20 meters and have low density, implying the availability of very few spectral channels with information and a considerable mixing of spectral contributions in each image pixel.
Appropriate and improved atmospheric and sunglint correction techniques were applied to the HS and MS data.Next, a water radiative transfer model was also considered to remove the water column disturbances and to generate the seafloor albedo maps.A preliminary analysis was performed to identify the most suitable preprocessed imagery to be used for seabed classification.Three different supervised classifiers (maximum likelihood, support vector machines, and spectral angle mapper) were tested.
A detailed analysis of different feature extraction methods was performed with the goal to increase the discrimination capability of the classifiers.To our knowledge, the effect of three rotation transforms to generate benthic maps was assessed for the first time.Texture parameters were, as well, added to check whether spatial and context information improve classifications.Finally, the inclusion of abundance maps for each cover, obtained by the application of linear unmixing algorithms, was also considered but, given the small number of spectral bands actually reaching the seafloor, results were not fully satisfactory.The best results were produced by SVM and the OBIA approach.However, to generate benthic habitat maps, the simple ML has shown an excellent performance and superior

Figure 1 .
Figure 1.Maspalomas area: (a) Geographic location; (b) Panoramic view of the study area (scale is approximate).
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 21 adjacency effect, and are the diffuse and direct transmittances and ̅ is the average reflectivity contribution from the pixel background [44].

Figure 4 .
Figure 4. Flowchart of the processing methodology to generate benthic maps.

Figure 7 .
Figure 7. AHS color composite (RGB for the original bands and the 3 first components for the transforms) and the first bands for each transform: (a) original bands; (b) Principal Component Analysis (PCA); (c) Independent Component Analysis (ICA); (d) Minimum Noise Fraction (MNF).

Figure 7 .
Figure 7. AHS color composite (RGB for the original bands and the 3 first components for the transforms) and the first bands for each transform: (a) original bands; (b) Principal Component Analysis (PCA); (c) Independent Component Analysis (ICA); (d) Minimum Noise Fraction (MNF).

Figure 8 .
Figure 8. ML and SVM overall accuracies for both sensors and the different input combinations (LU: Linear Unmixing, B: Bands, Text: Texture.The number of input bands appears in parentheses).

Figure 8 .
Figure 8. ML and SVM overall accuracies for both sensors and the different input combinations (LU: Linear Unmixing, B: Bands, Text: Texture.The number of input bands appears in parentheses).

Figure 9 .
Figure 9. Objects after the segmentation for rocky and sandy areas of the seafloor: (a) AHS; (b) WV-2.

Figure 9 .
Figure 9. Objects after the segmentation for rocky and sandy areas of the seafloor: (a) AHS; (b) WV-2.

Figure 10 .
Figure 10.AHS pixel-based classification (majority 5 × 5): (a) ML using bands 1 to 8; (b) SVM using the bands plus the texture of the first principal component.

Table 2 .
Overall accuracy (%) of Maximum Likelihood (ML), Support Vector Machine (SVM) and Spectral Angle Mapper (SAM) for Airborne Hyperspectral Scanner (AHS) and Worldview-2 (WV-2), and the different input combinations after each correction stage (AC: Atmospheric correction, SC: Sunglint Correction and WCC: Water Column Correction.Best accuracies marked in bold).

Table 2 .
Overall accuracy (%) of Maximum Likelihood (ML), Support Vector Machine (SVM) and Spectral Angle Mapper (SAM) for Airborne Hyperspectral Scanner (AHS) and Worldview-2 (WV-2), and the different input combinations after each correction stage (AC: Atmospheric correction, SC: Sunglint Correction and WCC: Water Column Correction.Best accuracies marked in bold).

Table 3 .
Overall accuracy (%) of ML, SVM, and SAM for AHS and WV, and the different input combinations (Best accuracies marked in bold).