Shallow water bathymetry inversion using hybrid multi-band composites resulted from integration of drone-based RGB and multispectral imagery

Shallow bathymetry inversion algorithms have long been applied in various types of remote sensing imagery with relative success. However, this approach requires that imagery with increased radiometric resolution in the visible spectrum is available. The recent developments in drones and camera sensors allow for testing current inversion techniques on new types of datasets. This study explores the bathymetric mapping capabilities of fused RGB and multispectral imagery, as an alternative to costly hyperspectral sensors. Combining drone-based RGB and multispectral imagery into a single cube dataset, provides the necessary radiometric detail for shallow bathymetry inversion applications. This technique is based on commercial and open-source software and does not require input of reference depth measurements in contrast to other approaches. The robustness of this method was tested on three different coastal sites with contrasting seafloor types. The use of suitable end-member spectra which are representative of the seafloor types of the study area and the sun zenith angle are important parameters in model tuning. The results of this study show good correlation (R 2 >0.7) and less than half a meter error when they are compared with sonar depth data. Consequently, integration of various drone-based imagery may be applied for producing centimetre resolution bathymetry maps at low cost for small-scale shallow areas.


Introduction
Shallow water bathymetry retrieval using optical imagery, is a field of ongoing research which has been greatly expanded in recent years. Consequently, the technique of satellite-derived bathymetry has seen significant growth with plentiful applications (Kutser et al., 2020).
There are two major groups of algorithms for deriving bathymetry from optical imagery, namely the empirical and the analytical methods. The empirical algorithms are mostly based on the models suggested by (Lyzenga, 1978) and Stumpf et al. (2003) which are implemented in various contexts (Geyman and Maloof, 2019;Gholamalifard et al., 2013;Ma et al., 2014;Traganos et al., 2018) and rely on the availability of ground-truth depth measurements for model calibration compared to the analytical methods. The empirical methods do not necessarily require absolute radiometric and atmospheric corrections (Gholamalifard et al., 2013;Kibele and Shears, 2016) and depending on model performance, they can be applied on datasets with similar seafloor types (Caballero and Stumpf, 2019), whereas analytical methods account for any seafloor type. The empirical algorithms are suitable for cases when imagery comprises of up to three bands in the visible range of the spectrum.
The analytical algorithms have been developed using in-situ calibrated spectral data which are fitted with radiative-transfer models (Dekker et al., 2011;Lee et al., 1999;Mobley et al., 2005). These algorithms are considered more suitable for imagery consisting of multiple bands (>4) in the visible spectrum, they do not require input of a priori depth information and they account for the inherent optical properties (IOPs) of water and bathymetric uncertainty as well, in contrast with the empirical methods. Ideally, the analytical algorithms perform well on hyperspectral imagery data with optimal quality and with appropriate radiometric and atmospheric corrections applied (Dekker et al., 2011;Castillo-López et al., 2017). A review from Dekker et al. (2011) evaluated the performance of several analytical algorithms for shallow water bathymetry retrieval using airborne hyperspectral data. These algorithms produced good results with low residuals (<~1 meter) for depths up to 10 meters. For greater depths the residual error was increasing with water depth.
The recent developments in drone technology offered tremendous opportunities for the development of novel geospatial applications. Drones are becoming increasingly popular in remote sensing studies since they are low-cost platforms; they provide centimeter-scale spatial resolution that is suitable for observing objects and/or processes in unique detail; they require negligible logistic effort, allowing for frequent deployment on demand, thus increasing the temporal resolution of imagery; they operate in close range without being influenced by clouds or other atmospheric effects (Alevizos, 2019;Román et al., 2021;Rossi et al., 2020).
There have been a few recent studies applying empirical SDB algorithms on drone-based multispectral imagery (Kabiri et al., 2020;Parsons et al., 2018;Rossi et al., 2020) showing relatively good results. Additionally, Chirayath and Earle (2016) developed a pioneering sensor for seafloor mapping providing active compensation for refraction and other optical distortions due to waves on sea surface.
Considering proximal sensing, there are several commercial sensors for drones, such as the MicaSense-RedEdge© dual camera offering 10 bands in the visible and near infrared areas, and other lightweight RGB or hyperspectral cameras which are routinely used on various projects such as environmental mapping (Manfreda et al., 2018), water quality monitoring (Isgró et al., 2021) and intertidal/subtidal habitat mapping (Fallati et al., 2020;Manfreda et al., 2018;Murfitt et al., 2017;Parsons et al., 2018;Román et al., 2021;Rossiter et al., 2020). Ideally, hyperspectral sensors would provide significant input data for shallow bathymetry mapping with analytical methods, due to their enhanced radiometric resolution in the visible range of the spectrum (Parsons et al., 2018). In practice however, their considerable cost and other hardware-related issues, probably do not allow them to be utilized extensively in drone-based projects (Rossiter et al., 2020).
A similar approach with the one presented here was applied by Tait et al. (2019) for classifying microalgae habitats. In their study they combined drone-based RGB and multispectral imagery together in order to enhance the supervised classification accuracy of microalgae taxa but without including bathymetry calculation. Apart from the presented study, to our knowledge there have not been any studies published currently, where shallow-water analytical models are applied on drone-based imagery. With an increasing availability of various camera sensors for drones, there are further opportunities for developing novel approaches regarding shallow bathymetry mapping at landscape scale.
The goal of this study is to provide an experimental alternative solution to costly hyperspectral sensors for shallow bathymetry applications without the need of a priori depth information. The study is based on the novel concept of integrating low-cost, drone-based RGB and multispectral imagery resulting in a mutli-band cube which can be utilized in shallow bathymetry inversion algorithms. This approach requires that both RGB and multispectral datasets have complementary spectral responses, they are radiometrically calibrated and that suitable endmember seafloor cover spectra are available. Therefore, we test the effectiveness of this new approach in three coastal sites with contrasting seafloor types and water qualities. Bathymetry outputs are validated with sonar measurements obtained from an unmanned surface vehicle (USV).

Study areas and fieldwork
In this study, we deployed a DJI Phantom 4 Pro drone equipped with a 1-inch, 20-megapixel CMOS sensor and also with a MicaSense RedEgde-MX © multispectral camera. A synopsis of the drone survey acquisition characteristics for each study area is shown in Table 1.
The following study areas have been selected for assessing the performance of hybrid composite imagery regarding shallow bathymetry retrieval (Fig.1). The first area is the small bay of Lambayanna located in Argolida region (Peloponesse, Greece) which comprises of smooth seafloor covered with medium sand and with some beach-rock outcrops in places. This area has been studied thoroughly in recent years due to its significance as a submerged prehistoric site (Papadopoulos et al., 2021). Particularly, in Papadopoulos et al. (2021 performed an extensive shallow bathymetric survey that took place using an inflatable boat yielding more than 10,000 bathymetry points (Fig.2).
The second area is the Kalamaki beach located in the western part of Chania city (Crete, Greece). It is a semi-enclosed bay with smooth, sandy seafloor covered with scattered pebbles and beachrock in places which are further covered with algae. The third study area, is the Plakias bay located in the southern part of Rethymno region (Crete, Greece). Plakias bay forms a shallow embayment with numerous rocky outcrops and gravelly sand on the seafloor. All depth measurements were acquired with an Ohmex BTX single-beam sonar with an operating frequency of 235 kHz. The sonar is integrated with a Real-Time Kinematics (RTK) GPS sensor for collecting attitude-corrected bathymetry points at 2 Hertz rate. The RTK-GPS measurements provide high spatial accuracy (< 10 cm) which is essential in processing drone-based imagery with a pixel resolution of a few centimetres. Regarding the Kalamaki and Plakias areas, the bathymetric surveys took place with a remotely controlled USV. An overview of depth data from each area is shown on the boxplot graphs in Fig.2.

Pre-processing of drone-based imagery
The diagram in Fig.3 shows the steps of the acquisition, processing and analysis pipeline followed in this study. Both sensors were set to collect nadiral images at one second interval for maximizing overlap between adjacent images and assisting with the mosaicking process. Although the MS sensor records five spectral bands simultaneously (Blue, Green, Red, Red edge and Near infrared) in this study we only considered the Blue, Green and Red bands from the visible spectrum since these are more favourable in optical bathymetry studies (Albert, 2004) and their spectral characteristics are complementary with the P4P bands ( Fig.4; Table 2, Appendix). The MS sensor was integrated with an external Downwelling Light Sensor (DLS-2) module which records sun illumination parameters (i.e.: angle, radiance) that are stored in the imagery metadata. These recordings are required during radiometric correction processing of multispectral imagery in Pix4D © software. In addition, the DLS-2 module provides GPS and attitude information for each acquired image, assisting the georeferencing and mosaicking of processed imagery using the Pix4D © software. In order to convert the pixel values to reflectance values, we acquired images of a spectral calibration panel which is specifically provided for the MS sensor and has a known reflectance coefficient for each band. In this way, the final processed data are suitable for quantitative analysis. The spectral responses of the P4P and the MS sensors are shown in Fig.4. Both RGB and multispectral datasets were processed with Pix4D © software for subsequent radiometric and geometric corrections. Alevizos & Alexakis (2021, in press) suggest that radiometric corrections are also required for drone RGB imagery for optimal results in shallow water bathymetry mapping. The RGB images were adjusted for radial lens distortion using the specific camera model included in the Pix4D © software and individual band mosaics (Red, Green, Blue) were exported. A reference reflectance value of 0.51 was set for all bands (both RGB and MS) considering that this value accounts for all wavelengths in the visible spectrum (MicaSense©, personal communication by email, 03/11/2020). Regarding the images captured with the MS camera, these were corrected for radial lens distortion using the respective camera model provided by the Pix4D © software. The multispectral imagery was radiometrically corrected using the reference spectral panel and the information collected by the on-board sun illumination sensor of the camera. The main points of radiometric calibration procedure regarding the MS sensor are described on the Pix4D © software website (Pix4D © , n.d.).
After the pre-processing stage, both RGB and MS reflectance mosaics were resampled at 20 cm pixel size and stacked together resulting in a six-band composite cube. The cube was converted to ENVI standard format for processing with the open-source WASI-2D software. In this study atmospheric correction of drone based imagery was not necessary. This is due to the fact that the drone surveys took place at significantly low altitude (<150 m) and with optimal weather conditions, so that atmospheric effects on recorded reflectance are totally negligible.

Shallow bathymetry inversion in WASI-2D
The WASI software is one of the few open-source tools for analyzing the spectral properties of aquatic environments. The WASI tool was initially designed for studying the water properties of fresh water environments and it has been applied on a limited number of bathymetry studies, mainly in lake environments so far (Gege, 2014a;Dörnhöfer et al., 2016;Niroumand-Jadidi et al., 2020). The software is based on earlier bio-optical models developed by Albert (2004), Albert & Mobley (2003), and Gege & Albert (2006). WASI has a 2D module which allows for image analysis on per-pixel basis (Gege, 2014b). This is particularly useful for analyzing imagery from multi-or hyper-spectral sensors. Regarding bathymetry retrieval, the WASI tool considers the influence of water-column constituents (IOPs) along with combinations of end-member seafloor reflectance spectra on water-leaving reflectances. For the Lambayanna area we used one of the default endmember spectrum embedded in WASI (Fig.4B). This is a sand type spectrum that it was measured by Pinnel, (2007) at shallow depth (0.5 m) using a submersible RAMSES spectroradiometer at the coast of Bolivar (south Australia).
At the study areas of Kalamaki and Plakias, we considered to replace some default end-member spectra of WASI with more appropriate ones that are similar to the seafloor types that naturally occur there. Thus, we imported the end-member spectra shown in Fig.4B which include rock, turf algae and brown algae, collected underwater with a spectroradiometer and a reference reflectance panel at various coastal locations in the southwest Indian Ocean (Mouquet & Quod (2010).
Suitable initial values of geometric (i.e.: sun zenith angle) and irradiance model parameters ( Table  2, Appendix) are required for accurate fitting of the spectral signatures. Once the model is tuned, the depth and seafloor type are fitted using the least squares method iteratively. The modelled spectral signature showing the lowest residual with the observed signature is used to determine the depth and seafloor type for each pixel. A detailed description of the WASI 2D tool can be found in Gege (2014).
The remote sensing reflectance in WASI is modelled according to the equations of Albert and Mobley (2003) and Albert (2004): The superscript sh indicates shallow water, deep deep water, b bottom, and the symbol λ indicates the wavelength. The first term on the right-hand side is the contribution of water column with depth zB, the second term represents the contribution of the bottom albedo. Light attenuation is described by the attenuation coefficients Kd for down-welling irradiance, KuW for upwelling radiance originating from the water layer, and KuB for upwelling radiance from the bottom surface. These three coefficients are calculated as a function of the sun zenith angle, viewing direction and the concentrations of water constituents using equations also derived by Albert and Mobley (2003) and Albert (2004). Ars,1 and Ars,2 are empirical constants.
The WASI2D algorithm iterates the spectral signatures on per pixel basis trying to fit an optimal spectrum given the constant values of model parameters (Table 3, Appendix). Inverse modeling takes place by approximating the Rrs spectra (of each pixel) with suitable WASI spectra for different depths. The best fit with the observed image spectrum is obtained by minimizing a cost function that calculates the correlation between the Rrs and the WASI spectra. The inversion algorithm employs the absolute difference function in order to identify an optimal set of fit parameters (depth and seafloor type), which minimize the residual of the cost function (Niroumand-Jadidi et al., 2020).

Results
Bathymetry outputs resulted from shallow water inversion in WASI-2D were validated using insitu sonar measurements (Fig.1, Fig.2B). At the Lambayanna area, the WASI inverted depths show very good correlation with sonar data (R 2 =0.82), and have a mean absolute error (MAE) of 0.41 m (Fig.5B). Considering the overall homogeneity of the study area only one end-member spectrum (WASI sand) was used for inversion within the range of 400-700 nm. The initial model parameters (e.g.: chlorophyll-a [CHL-a]; suspended particulate matter [SPM]) were estimated by fitting a few (2-3) single pixels with known depth. The maximum depth estimated by inversion is -5.5 meters. The resulted bathymetry at the Kalamaki area (Fig. 6B) fits very well (R 2 =0.87) with USV sonar measurements and has a MAE of 0.42 m. The end-member spectra of WASI-sand, turf algae and brown algae were selected as more suitable for driving the bathymetry inversion in this area since they correspond with the seafloor cover composition at the Kalamaki area. Inversion took place at the range of 420-600 nm since the pixel values of the MS-red band (660 nm) were too high for the very shallow part (>-1m) and they were not providing useful reflectance values for deeper (<-2m) water. The regression plot in Fig.6B shows that there is not an apparent increase in error with increasing depth. The maximum depth estimated by inversion is -5.5 meters. The bathymetry inversion results for the Plakias area (Fig.7B) show good correlation (R 2 =0.75) with USV sonar measurements and have a MAE of 0.44 m. The bathymetry error seems to be distributed evenly across all depths. Here three end member spectra (rock, brown and turf algae) were applied for inversion since these two are more prominent in the study area. Inversion took place at the range of 420-600 nm since the pixel values of the MS-red band (660 nm) were too high for the very shallow part (>-1m) and they were not providing useful reflectance values for deeper (<-2m) water. The maximum depth estimated by inversion is -7.5 meters. Analysis of the bathymetry residuals (absolute difference between reference depth and inverted depth) shows that residuals with large difference values (>2m) are clearly associated with extremely low values of the respective multispectral bands (Fig.8).

Discussion
This experimental approach was based on leveraging open-software and literature data (spectra) for producing the first-of-its-kind shallow water bathymetry inversion. The results of this study highlight the effectiveness of hybrid multiband composites for proximal sensing of shallow water bathymetry. Optimal bathymetry results were obtained for all three study areas regardless of the optical complexity of the seafloor types.
The combination of RGB and MS imagery provides a cube with improved spectral resolution in the visible spectrum which is suitable for shallow water inversion. The issue of spectral resolution in shallow bathymetry inversion studies has been highlighted by Mobley et al. (2005) and Leiper et al. (2014). The number of spectral bands determines the detail (and hence the accuracy) at which the seafloor reflectance is approximated and assists in estimating the IOPs values. In addition, only bands sensitive to visible wavelengths are suitable for bathymetry inversion since those are transmitted (some more than others) through the water column without being absorbed or scattered. Particularly blue and green wavelengths penetrate deeper into the water column compared to yellow and red wavelengths that are easily absorbed by the first 1-2 meters of water depth (Eugenio et al., 2017;Kobryn et al., 2013). Our findings suggest that any combination of dual sensors (either RGB or MS) with complementary bands in the visible spectrum, are suitable for shallow bathymetry inversion, provided that they are radiometrically calibrated.
It has to be mentioned that it was not necessary to scale the cube reflectance values by pi for the areas of Kalamaki and Plakias. This is a standard practice for obtaining remote sensing reflectance (Rrs) values from atmospherically corrected satellite imagery (Cillero Castro et al., 2020;Dierssen et al., 2003). Considering that the case of satellite acquisition geometry is large (compared to imagery footprint) the behaviour of Earth's surface radiance may be characterized as Lambertian (i.e.: transmitted towards all directions) at this scale (Gao et al., 2009). Consequently, scaling the satellite data by pi is necessary for adjusting their reflectance values so they match the Rrs values of the spectral libraries which are used by shallow bathymetry inversion models.
When a pi factor was applied in Kalamaki and Plakias areas data cubes, the bathymetry inversion was failing. This model behaviour suggests that drone imagery captures Rrs directly and further scaling using pi is unnecessary. This can be explained by the fact that when the sun zenith angle is relatively small (0-50 degrees), seafloor albedo reflected from the water/air interface may be characterized by non-Lambertian reflectance (Garcia, 2015). Therefore it is suggested that drone acquisition during midday time provides optimal conditions for capturing original Rrs values. The same is not valid for satellite imagery due to the scale of the acquisition geometry (i.e.: altitude).
Additionally, the close range of the drone to the seafloor assists in capturing Rrs without the need for applying atmospheric corrections on drone imagery.
Only at the Lambayanna area the cube reflectance values were scaled by pi. This is probably explained by the fact that the drone images were acquired early in the morning (9 am) and the seafloor albedo was reflected in a more Lambertian pattern (above the sea-surface) due to the high sun zenith angle (70 degrees).
End-member spectra are a very important input in shallow bathymetry inversion (Dörnhöfer et al., 2016;Garcia et al., 2018). In the case of Kalamaki beach it seems that although deeper (>-2m) parts of the scene are covered with brown algae, bathymetry error is not increasing over those areas. This highlights the robustness of the input end-member spectra for inferring bathymetry in this area. Thus, it is hypothesized that there is a considerable shape similarity between endmember spectra used from literature and those naturally occurring at the Kalamaki area. At the Lambayanna area, only one end-member spectrum (WASI sand) was sufficient for yielding bathymetry results with low error. This is attributed to the overall homogeneity of the seafloor type in this area. In contrast, at the Plakias area we applied three end-member spectra in order to approximate better the seafloor types occurring in this complex area.
It is suggested that smooth seafloor geometry at the Lambayanna and Kalamaki areas assisted in calculating more accurate bathymetry inversion than the Plakias area. Smooth seafloor provides a more unobstructed and uniform transmission of seafloor albedo resulting in capturing a clear reflectance signal by drone-based imagery. In contrast, the rugged seafloor surface at Plakias is probably responsible for the observed deviations in bathymetry results. In order to examine this hypothesis we classified the bathymetry residuals at the Plakias area according to their magnitude and we extracted their corresponding reflectance values for each band. The reflectance boxplots for each residual class and for all bands are shown in Fig.8. These plots reveal that residuals with maximum deviation from actual depth, are associated with extremely low reflectance values across all bands of the cube. Consequently, it is inferred that seafloor rugosity along with cracks and crevices are blocking sunlight and cast shadows in places, causing significantly lower reflectance values than their surroundings. These obscured areas, result in considerable overestimation of depth by the inversion model which falsely fits them with spectral signatures corresponding to deeper waters.
In addition, the results at the Plakias area showed increased errors over seagrass patches (isolated blue spots in Fig.7A) and therefore were not include in the scatterplot in Fig.7B. This issue was encountered even when a seagrass end-member was used in the model. The reason for this problem may be the fact that the characteristic seagrass spectrum beyond the 650 nm range becomes rapidly absorbed even from one meter depth, leading to confusion with other spectra (e.g.: rocks) corresponding to deeper water. Considering that in Plakias bay, seagrass patches occur in areas with at least 2 meters depth, it is understood that absorption in the water column, leads to depth overestimation.
This model behaviour can be explained by the effect of 'non-uniqueness' that is a common issue hindering bathymetry retrieval from hyperspectral data. This problem has been described by Defoin-Platel and Chami (2007), Mobley et al. (2005) and Garcia (2015) and is related to the natural phenomenon where different combinations of a) optically-active constituents in the water-column and/or b) seafloor albedos can produce the same spectral signature, challenging the accuracy of the final SDB results. In the case of Plakias area, the seagrass patches are fitted with very deep water signatures because they cannot be differentiated by the rest of the seafloor types as their spectral features beyond the 650 nm are absorbed.
Shallow bathymetry analytical models require careful tuning of initial parameters in order to produce accurate results. By selecting a smaller number of suitable parameters, the model is better constrained and requires less computational resources (Gege, 2014b;Niroumand-Jadidi et al., 2020). Thus, when additional information is known (e.g.: depth) for a particular location on the image, then model inversion should be applied on individual pixels first, in order to determine the initial values of model parameters (Table 3, Appendix). Ideally, spectral measurements from a deep water area (where seafloor albedo is negligible) should be collected in order to provide more realistic initial values for CHL-a and SPM (Niroumand-Jadidi et al., 2020). Alternatively, these data could be extracted from coarse resolution satellite imagery or from available literature.

Conclusions
This study shows that the integration of drone-based RGB and multispectral imagery provides promising results in shallow bathymetry inversion. Both imagery types are required to be radiometrically corrected and to have complementary band-widths in the visible range of the spectrum. The multi-band composite cube provides a low-cost alternative to costly hyperspectral commercial sensors, allowing for a wider field of shallow water applications. In addition, the cube is suitable for analysis with physics-based optical models, reducing further the costs for obtaining in-situ bathymetry data for model calibration. Suitable end-member spectra of local seafloor types should be carefully selected as inputs to shallow analytical algorithms for maximizing the quality of bathymetry output. The presented approach works well across various water types and for relatively smooth seafloor, while seafloor roughness found to contribute larger bathymetric errors. This novel method is appropriate for high resolution bathymetry mapping at landscape scale when only drone imagery is available from areas with sufficient water transparency.