An Approach to Minimize Atmospheric Correction Error and Improve Physics-Based Satellite-Derived Bathymetry in a Coastal Environment

Abstract

Physics-based radiative transfer model (RTM) inversion methods have been developed and implemented for satellite-derived bathymetry (SDB); however, precise atmospheric correction (AC) is required for robust bathymetry retrieval. In a previous study, we revealed that biases from AC may be related to imaging and environmental factors that are not considered sufficiently in all AC algorithms. Thus, the main aim of this study is to demonstrate how AC biases related to environmental factors can be minimized to improve SDB results. To achieve this, we first tested a physics-based inversion method to estimate bathymetry for a nearshore area in the Florida Keys, USA. Using a freely available water-based AC algorithm (ACOLITE), we used Landsat 8 (L8) images to derive per-pixel remote sensing reflectances, from which bathymetry was subsequently estimated. Then, we quantified known biases in the AC using a linear regression that estimated bias as a function of imaging and environmental factors and applied a correction to produce a new set of remote sensing reflectances. This correction improved bathymetry estimates for eight of the nine scenes we tested, with the resulting changes in bathymetry RMSE ranging from +0.09 m (worse) to −0.48 m (better) for a 1 to 25 m depth range, and from +0.07 m (worse) to −0.46 m (better) for an approximately 1 to 16 m depth range. In addition, we showed that an ensemble approach based on multiple images, with acquisitions ranging from optimal to sub-optimal conditions, can be used to estimate bathymetry with a result that is similar to what can be obtained from the best individual scene. This approach can reduce time spent on the pre-screening and filtering of scenes. The correction method implemented in this study is not a complete solution to the challenge of AC for satellite-derived bathymetry, but it can eliminate the effects of biases inherent to individual AC algorithms and thus improve bathymetry retrieval. It may also be beneficial for use with other AC algorithms and for the estimation of seafloor habitat and water quality products, although further validation in different nearshore waters is required.

1. Introduction

Bathymetric information from satellite data is of fundamental importance in optically shallow waters, where the seafloor is visible from space and the water-leaving radiance (L_w) is influenced by reflection off the seafloor. Such information, in the form of maps of water depth, is essential for a wide variety of purposes including offshore activities (e.g., pipeline laying), resource management (e.g., fishery), and defense operations (e.g., navigation). Traditional bathymetric charts are based on soundings obtained during hydrographic surveys. However, as ship-borne surveys are costly and time-consuming, and many shallow-water environments are highly dynamic, it is impossible to survey all areas of interest, and the difficulty in accessing shallow and remote areas means that in practice, up-to-date data are typically only available for limited areas (harbors and main navigation corridors). Airborne light detection and ranging (LiDAR) Bathymetry (ALB) systems, such as CZMIL (Coastal Zone Mapping and Imaging LiDAR) [1], LADS MK 3 (Laser Airborne Depth Sounder MK 3) [2], and EAARL-B (Experimental Advanced Airborne Research LiDAR B) [3] can also be used to map water depth. With these techniques, a vertical accuracy of about ± 15 cm in shallow water is possible [4], although accuracy is affected by turbidity and the LiDAR system. While precise bathymetric mapping of water depth to about 20–70 m depth can be achieved with airborne LiDAR [5,6], costs associated with these systems are relatively high, thus limiting their application over large, or remote, areas.

Passive optical satellite remote sensing can also be used to map bathymetry, typically known as satellite-derived bathymetry (SDB), based on the relationship between the color of a shallow-water area and the depth of the water. SDB can be implemented using empirical or physics-based methods. The empirical methods are based on the simple premise that a statistical relationship can be established between water depth and the remotely sensed radiance of a water body, using regression or similar analysis [2,7,8,9,10]. Thus, all empirical approaches require coincident in-situ data on water depth for calibration; ideally, these data should be up to date and have good geographic and depth distribution. Empirical approaches assume that the inherent optical properties (IOPs) of the water, as well as seafloor spectral reflectance, do not vary across the image, and therefore, the results may contain large errors and require manual editing when this is not the case. A key advantage of empirical approaches is the ability to retrieve water depth relatively easily, but their reliance on calibration from coincident field observations means that they cannot be used for systematic regional and global mapping and monitoring. Physics-based methods instead estimate bathymetry on per-pixel basis through the inversion of a radiative transfer model (RTM). As such, they do not assume uniform IOPs and seafloor reflectance, nor do they rely on coincident depth data for calibration. In addition to bathymetry, seafloor reflectance and water IOPs, which can be used to infer substrate and water quality respectively, can be simultaneously retrieved, and per-pixel uncertainties of all these parameters, including water depth, can also be determined. While originally developed for and tested on airborne hyperspectral imagery, physics-based methods for SDB have also been demonstrated for multispectral satellite sensors [11,12,13,14]. Physics-based methods can be implemented using either look-up tables (LUTs) [15,16] or semi-analytical optimization methods [17,18]. In the first case, a database of remote sensing reflectance (R_rs) spectra is built from an RTM provided with a range of values for water depth, spectral seafloor reflectance, water column optical properties (absorption and backscattering coefficients), and known environmental conditions such as sun angle and wind speed. For the retrieval of parameters (water depth, water IOPs, and seafloor reflectance) in each image pixel, a search is then performed to find the R_rs in the LUT that best matches the one observed in the pixel. With semi-analytical optimization methods, the radiative transfer equation is used to estimate water depth by iterative optimization of the same parameters. In both methods, the best match between modeled and observed reflectance is determined using a least squares or similar matching technique.

Despite the advantages of physics-based methods, a substantial challenge is that they rely on precise estimates of absolute radiometry, typically in the form of R_rs or L_w. Unlike other optical remote sensing applications, including the empirical approaches to satellite-derived bathymetry, physics-based retrieval algorithms may perform very poorly if R_rs is incorrectly estimated, and high-quality R_rs data from a robust atmospheric correction (AC) is essential for accurate physics-based water depth estimation. Accordingly, a variety of AC algorithms have been developed for ocean color (OC) products retrieval such as bathymetry, and several studies have validated their performance against in situ data. For example, Pahlevan et al. [19] validated R_rs produced from different AC schemes in the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) Data Analysis System (SeaDAS) with in situ data from the AERONET-OC network. Likewise, Doxani et al. [20] assessed the performance of different AC methods and validated their R_rs with match-up datasets over both land and water surfaces in an AC inter-comparison exercise. Warren et al. [21] evaluated the accuracy of a wide range of freely available AC processors by comparing them to reference R_rs data from different coastal and inland waters. Similarly, in a more recent AC exercise, Zhang and Hu [22] also analyzed an AC algorithm, comparing its R_rs images with those measured over a few sites from the AERONET-OC stations. Collectively, these studies demonstrated that accurate AC remains a challenge for OC remote sensing where precise R_rs data are needed. Therefore, it is important to explore ways by which errors in AC outputs, and their effect on the products derived from them, can be minimized. One way to address some of the problems posed by imprecise AC is to assess and quantify the impacts of environmental variables on AC accuracy and then account for this in the atmospherically corrected image. In an earlier study [23], four publicly available AC processors (2 land-based and 2 water-based) for deriving the R_rs in coastal waters were compared and validated with 54 R_rs match-up datasets from AERONET-OC stations. The study revealed that biases from ACOLITE and SeaDAS, two of the state-of-the-art AC algorithms, are influenced by environmental variables. In this study, we demonstrated the potential of Landsat 8 (L8) data for SDB in US coastal waters and assessed the performance of a commonly used and publicly available water-based AC algorithm (ACOLITE [24]) for physics-based SDB. To minimize the effect of imperfect AC on the bathymetry retrieval, we further used a correction factor to improve the original atmospherically corrected image from ACOLITE. Using a set of 9 images, SDB estimates from these two AC procedures were then compared with LiDAR-derived bathymetry of the area. Lastly, we used an ensemble approach to produce SDB of the study area using all the corrected images.

2. Study Sites and Imagery

2.1. Study Sites

The Florida Keys is a series of islands that extend from the southern end of Florida, USA, to the south–southwest. Their nearshore shallow waters include coral reef tracts, patch reefs, bank reefs, seagrass meadows, and unvegetated hard and soft bottom. This site was chosen because of its relatively clear waters, the good knowledge of seabed features, and availability of LiDAR-derived depth data for validating SDB estimates of water depth. The benthic environment of the section of the Florida Keys used in this study is dominated by extensive seagrass beds, with some patches of reef and unconsolidated sediments. Figure 1 shows this area with the distribution of bathymetric LiDAR data used for validating the SDB estimates in this study.

2.2. Satellite Data

Nine L8 images (Figure 2) from the Florida Keys, acquired during both optimal and near-optimal conditions for SDB, were downloaded from the archive of the United States Geological Survey after visually inspecting all available images from May 2013 to May 2019. L8 OLI (Operational Land Imager) collects visible, Near Infrared (NIR) and Short-wave infrared (SWIR) spectral band imagery at 30 m spatial resolution. In addition to the improved positional accuracy of 14 m, compared to 50 m for its predecessors in the Landsat series, L8 includes coastal and aerosol (433–453 nm) and blue (450–515 nm) bands for coastal and bathymetric mapping [25,26].

2.3. LiDAR Data

To validate the SDB estimates, a bathymetry topographic digital elevation model (DEM) was acquired from the National Oceanic and Atmospheric Administration (NOAA) National Centres for Environment Information (NCEI) coastal LiDAR archive. The LiDAR data collection was conducted in December 2014 over South Florida and the Florida Keys as part of efforts by NOAA to study sea level rise and coastal flooding impacts on US coasts. Several LiDAR sources including topographic and bathymetric LiDAR sensors were used to develop and create a suite of tiled bathymetric–topographic DEMs for South Florida and the Florida Keys [27]. A portion of the DEM tiles covering the study site (Figure 1) was retrieved from the Office of Coastal Managements Data Access Viewer [28] where all DEM data are archived. The DEMs, with a vertical accuracy of approximately 0.5 m, are referenced vertically to the North American Vertical Datum of 1988. Horizontal positions were provided in geographic coordinates and referenced to the North American Datum of 1983 [29,30]. A portion of the collection covering the Florida Keys coastal area was referenced to mean sea level and resampled from 0.3 m to 30 m to match the spatial resolution of L8.

3. Methodology

3.1. Data Preprocessing

3.1.1. Atmospheric Correction

We implemented two types of AC methods for water depth retrieval: (1) we used ACOLITE to process L8 images into R_rs values (henceforth Rrs_raw) and (2) then applied a correction factor to reduce errors in the original ACOLITE output and create new corrected R_rs values (henceforth Rrs_corrected). ACOLITE [24], specifically designed for AC over water surfaces, is an AC method that estimates L_w by simulating contributions from molecular (Rayleigh) and particulate (aerosol) scattering using a 6SV-based LUT [31]. Based on Ruddick et al. [32], aerosol reflectance is estimated by determining a per-tile aerosol type (or epsilon) from the ratio of reflectances in two bands over water pixels where L_w can be assumed to be zero. Then, the epsilon is used to extrapolate the observed aerosol reflectance to the visible bands to remove atmospheric contributions. ACOLITE was originally designed for processing L8 images, but it has been modified and updated to also process Sentinel-2 data [33]. Furthermore, the most recent version, which can be adapted to commercial sensors such as Pleiades, contains an additional AC scheme (now the default setting) called the dark spectrum fitting (DSF) algorithm, as well as a sun glint correction scheme [34]. In this study, ACOLITE (version 20170113.0) was used to produce all R_rs images, which are the direct input into the bathymetry algorithm. The default SWIR option (1609 and 2201 nm band combination) was implemented for all images. This band combination takes advantage of the longest-wavelength SWIR band, where water absorption is the highest. In a previous study [23], in which a range of AC algorithms were compared and validated against in situ L_w from 14 AERONET-OC stations, statistically significant relationships were demonstrated between errors in ACOLITE’s R_rs estimates for L8′s 443 nm and 482 nm bands and three environmental variables: Solar Zenith Angle (SZA), Aerosol Optical Thickness (AOT) at 865 nm (AOT₈₆₅), and wind speed (u₁₀); probable but statistically non-significant relationships were also demonstrated for the 561 nm and 655 nm bands. Using multiple linear regression, we therefore derived a set of coefficients that were used to estimate the error of ACOLITE’s R_rs estimates for each of those four bands in each image, as a function of SZA, AOT₈₆₅, and wind speed. Then, each of the four bands used for depth retrieval in this study was corrected using Equation (1):

Rrs_corrected = Rrs_raw − (a + b*SZA + c*AOT₈₆₅ + d*u10)

(1)

where Rrs_corrected and Rrs_raw are the R_rs images with and without correction, respectively; and a, b, c, and d are coefficients obtained through fitting a linear model to the data from Ilori et al. [23]. SZA was obtained from the metadata of each L8 scene. AOT₈₆₅ was processed and obtained using the l2gen processor in the SeaDAS software, and an average value used for each image was calculated by randomly sampling multiple pixels over the area of the study site. Wind speed data were obtained from the National Centers for Environmental Prediction Reanalysis project [35], where 6 h global wind speed estimates are archived.

Table 1. Environmental parameter variables for each image. AOT₈₆₅: Aerosol Optical Thickness (AOT) at 865 nm, SZA: Solar Zenith Angle.

Scene Date (dd/mm/yyyy)	SZA (Degrees)	AOT865	u10 (m/s)
01/12/2013	50.36	0.081	5.29
05/01/2015	52.79	0.088	1.07
26/01/2017	50.13	0.083	2.49
28/12/2017	52.98	0.076	6.45
13/01/2018	52.14	0.142	3.11
14/02/2018	45.63	0.12	4.84
02/03/2018	40.66	0.11	3.21
01/02/2019	49.02	0.122	3.74
05/03/2019	39.74	0.143	4.67

presents the value of each environmental parameter for each image used in this study.

3.1.2. Sun Glint Correction

As sun glint correction is not inherently part of the ACOLITE version used in this study, we implemented the NIR method [36] to remove specular reflection off the sea surface for images where glint was visually obvious. This method assumes that for optically deep areas (where radiation reflected from the seafloor has a negligible influence on L_w), any remaining NIR signal after AC must be due to sea surface reflection. Thus, glint intensity and removal is performed by establishing a linear relationship between the NIR and visible bands over an optically deep area in the image, and that relationship is then used across all water pixels to reduce R_rs for the visible bands to its assumed glint-free value.

3.1.3. Estimation of Noise Equivalent Reflectance

Bathymetry model inversion based on least squares optimization techniques is generally sensitive to environmental noise [37,38]; thus, high environmental noise may make images unsuitable for bathymetry extraction. The noise-equivalent difference in reflectance, NE∆R_rs (sr⁻¹), is a measure of image noise, with contributions from the sensor (e.g., instrument degradation) and the environment (e.g., variability in atmosphere and water surface state) [37,38]. The NE∆R_rs can be used to assess the suitability of a satellite imagery for aquatic remote sensing applications. For example, it has been used to determine the suitability of the Compact Airborne Spectrographic Imager (CASI) for benthic mapping [11]. Therefore, following AC, we estimated the NE∆R_rs (sr⁻¹) [39] by calculating the band-wise standard deviation of R_rs from a 33 × 33-pixel window over a homogeneous optically deep area using Equation (2) [40]. This approach assumes that any observed spectral variations in the selected area is due to noise; thus, selected pixels must be as homogenous as possible for a robust standard deviation estimate. Ideally, the NE∆R_rs should be lower than 0.00025 sr⁻¹ in each of the visible bands [41], which was the case for all nine images used in this study.

Table 2. The noise equivalent difference in reflectance (NE∆R_rs), computed from a kernel of 33 × 33 pixels from an optically deep and homogeneous area, for each image used in this study.

Scene Dates (dd/mm/yyyy)	Band 1	Band 2	Band 3	Band 4
01/12/2013	0.000200	0.000154	0.000096	0.000061
05/01/2015	0.000136	0.000108	0.000084	0.000063
26/01/2017	0.000092	0.000072	0.000057	0.000042
28/12/2017	0.000151	0.000105	0.000081	0.000053
13/01/2018	0.000111	0.000103	0.000069	0.000047
14/02/2018	0.000157	0.000129	0.000108	0.000063
02/03/2018	0.000126	0.000110	0.000069	0.000043
01/02/2019	0.000148	0.000127	0.000100	0.000063
05/03/2019	0.000086	0.000081	0.000059	0.000042

shows the per-band value obtained for each of the 9 images used in this study.

$$\mathrm{NE}∆\mathrm{R}{\mathrm{rs}}_{ }=\mathsf{\sigma }{\mathrm{R}}_{\mathrm{rs}}$$

(2)

where $\mathsf{\sigma }$R_rs is the standard deviation in each band over an as homogeneous as possible area of optically deep water within the image.

3.1.4. Parameterization of Environmental Properties

To implement the physics-based approach to SDB, values of optical properties and substratum spectral reflectance that are representative of the environment in question are needed. Water inherent optical properties (IOP) (P₄₄₀, G₄₄₀, and X₅₅₀) parameterization for forward modeling for each site was based on assessment from Level 3 OC products from the Visible Infrared Imaging Radiometer Suite Visible Infrared Imaging Radiometer Suite (VIIRS) Generalized Inherent Optical Property (GIOP) algorithms [42]. P₄₄₀ is the phytoplankton absorption coefficient at 440 nm, G₄₄₀ is the absorption of gelbstoff and detrital materials coefficient at 440 nm, and X₅₅₀ is the particulate backscattering of suspended particles coefficient at 550 nm. Using parameter values obtained from these OC products, ranges of values for each parameter were determined by observing the lowest and highest parameter values for all dates from GIOVANNI, which is an online visualization tool for OC products [43]. Then, values slightly lower and higher than the observed lowest and highest values, respectively, were chosen (

Table 3. Parameter ranges used for forward modeling.

P	G	X	Z
0.006–0.04	0.004–0.04	0.0005–0.006	0.1–25

). As part of the inversion model, seafloor reflectance spectra are also needed. We used two seafloor spectra (Figure 3), based on the area’s benthic description [44]. Depth (Z), which was also needed for forward modeling, was set to 0.1 and 25 m with the understanding that depth penetration greater than 25 m would be difficult.

3.1.5. Forward Modeling of Remote Sensing Reflectance

To derive water depth, we applied a modified version of the semi-analytical inversion model of Lee et al. [17,18] to the atmospherically corrected images. In this inversion scheme, the sub-surface remote sensing reflectance, r_rs, (the ratio of upwelling radiance to downwelling irradiance just below the surface) is related to absorption (a) and backscattering properties ($b_b$) of the water column, the seafloor reflectance (ρ), and water depth (H). For nadir-viewing satellites, the model can be expressed as:

$$R_{rs }\approx \frac{0.5r_{rs}}{1-1.5r_{rs}}$$

(3)

where r_rs, the subsurface remote-sensing reflectance, is expressed as:

$$\begin{gathered} r_{rs}\left(a,b_{b, }H, \rho \right)\approx \left(0.084+0.170u\right)u\left(\begin{array}{c}1-exp\left\{-\left[\frac{1}{\cos \left({\theta }_w\right)}+\frac{1.03\sqrt{1+2.4u}}{\cos \left({\theta }_v\right)}\right]kH\right\}\end{array}\right)+ \\ \frac{\rho }{\pi }exp\left\{-\left[\frac{1}{\cos \left({\theta }_w\right)}+\frac{1.04\sqrt{1+5.4u}}{\cos \left({\theta }_v\right)}\right]kH\right\} \end{gathered}$$

(4)

$$u = \frac{b_b}{\mathrm{a}+ b_b }$$

(5)

$$k = \mathrm{a} + {\mathrm{b}}_{\mathrm{b}}$$

(6)

where ${\theta }_w$ and ${\theta }_v$are the sub-surface solar zenith and sub-surface sensor viewing angles, respectively. Absorption (a) and backscattering coefficients ($b_b$) are functions of (1) the absorption coefficient of phytoplankton at 440 nm, P; (2) the absorption coefficient of colored dissolved materials at 440 nm, G; and (3) the backscattering coefficient of suspended particles at 550 nm, X. These are expressed as:

$$a= a_w+ P{a}^{*}{ }_{phy }+G{e}^{-0.015\left(\lambda -440\right)}$$

(7)

$$b_b= b_{bw}+X{\left[\frac{550}{\lambda }\right]}^Y$$

(8)

where $a_w$and $b_{bw}$ are the absorption and backscattering coefficients of pure water, respectively [45], $a^{*}{ }_{phy }$is the specific absorption of coefficient of phytoplankton (normalized to a value of 1.0 at 440 nm), λ is the center wavelength, and Y is the spectral shape that depends on the particulate shape and size.

While Lee’s inversion model uses the albedo of only one key benthic substrate (sand), our model includes a parameterization to set the seafloor reflectance as a linear mix of the two bottom types (i.e., sand and algae; [46]). To forward model the R_rs as a function of water depth, water quality parameters, and the seafloor reflectance, the adaptive look-up table (ALUT) method [11,16] was implemented, which ensures efficient construction and search through the table. In this approach, an LUT consisting of the modeled R_rs values of L8 bands 1–4, seafloor reflectance (Figure 3), water optical properties (absorption and scattering characteristics of water), and water depths of the optically shallow zone of the area in question (Table 3) is constructed. With realistic minimum and maximum values of all environmental parameters in the table, the LUT construction is optimized by using a hierarchical structure to efficiently cover the range of expected R_rs values while minimizing the under- or over-sampling of spectrally similar regions of environmental space, which is common with discretization by regular intervals in conventional forward models. For example, a small change in depth in shallow water areas will cause a significant/large change in R_rs, and will typically result in under-sampling if the depth parameter is discretized by regular intervals. Likewise, oversampling may occur in deep water areas where a small change in R_rs is expected [11,16]. To identify the best parameter values to be included in the hierarchy in the ALUT technique, discretization is based on an evenly sampled spectral space and not on an evenly sampled parameter space (e.g., water depth). This approach requires bounded ranges for all the modeled parameters, for which we used the value ranges in Table 3.

3.1.6. Inversion of Remote Sensing Reflectance

Model inversion was subsequently performed using the binary space partitioning (BSP) approach [11,16], as described in Knudby et al. [12]. Briefly, this technique subdivides the LUT created during forward modeling into different nodes. First, the BSP splits the whole LUT into two (left and right child nodes) and subsequently subdivides the nodes into a partitioning tree, which facilitates the optimization of the per-pixel LUT search. After model inversion and depth retrieval, water depths were corrected for tidal height at the time of each image acquisition using tidal height estimates obtained from Oregon State University’s tide prediction service [47].

Depth cannot be accurately estimated for optically deep pixels (i.e., where reflection off the seafloor is a negligible component of the radiance measured by the sensor); thus, we estimated a depth threshold value for each scene to distinguish between optically deep and optically shallow areas. These threshold values were obtained by calculating the mean minus 2 standard deviations of depth predictions within a 33 × 33 window of homogeneous optically deep water. Results are reported for each scene both a) for the full depth range, and b) for depths from the surface and down to these scene-specific thresholds.

3.2. Validation of Depth Estimates

Validation of depth estimates from the two AC procedures was performed by comparing the estimated depths to the LiDAR data. The number of depth estimates used for validation (

Table 4. Summary validation statistics for satellite-derived bathymetry (SDB) estimates (SDB_raw and SDB_corrected) for the two depth limits in this study. Bold letters in the root-mean-squared-error (RMSE) column indicate where an observed difference between SDB_raw and SDB_corrected estimates is more than 0.1 m.

Scene Date dd/mm/yyyy	Max Depth (m)	RMSE (m) (SBDraw/SBDcorrected)	Bias (SBDraw/SBDcorrected)	R2 (SBDraw/SBDcorrected)	Number of Validation Points
01/12/2013	25	1.96/1.66	0.83/0.84	0.66/0.79	3102
	13.80	1.80/1.50	0.85/0.73	0.58/0.71	2988
05/01/2015	25	2.02/2.03	−0.17/−0.46	0.71/0.71	3317
	11.23	1.84/1.88	−0.46/−0.71	0.53/0.53	3050
26/01/2017	25	1.34/1.25	0.49/0.09	0.85/0.85	3311
	12.03	1.21/0.95	0.57/0.23	0.81/0.85	3142
28/12/2017	25	1.87/1.39	1.16/0.41	0.74/0.81	3148
	17.54	1.85/1.39	1.18/0.42	0.79/0.81	3086
13/01/2018	25	2.21/2.04	1.45/1.40	0.66/0.76	2804
	15.21	2.19/1.89	1.48/1.35	0.56/0.74	2721
14/02/2018	25	1.86/1.72	−1.15/−1.04	0.77/0.79	3316
	14.81	1.80/1.69	−1.12/−1.01	0.71/0.74	3247
02/03/2018	25	2.03/1.65	0.71/0.68	0.61/0.74	3268
	16.84	1.98/1.62	0.69/0.67	0.53/0.71	3226
01/02/2019	25	1.35/1.05	0.55/0.62	0.83/0.92	3255
	14.61	1.33/1.02	0.54/0.64	0.76/0.90	3187
05/03/2019	25	1.36/1.27	0.76/0.53	0.89/0.89	3312
	16.39	1.32/1.25	0.79/0.55	0.85/0.83	3263

) varied between the nine images due to differences in the number of pixels for which depth was successfully estimated, as pixels that did not pass the AC’s internal quality checks (e.g., due to clouds), pixels with negative depths, and pixels that were visually impacted by boats, wake, or cloud shadows were eliminated prior to validation. Based on the remaining pixels, we used the coefficient of determination (R²), RMSE (root-mean-squared-error) (Equation (9)) and bias (Equation (10)) to compare the accuracy of the uncorrected and corrected SDB estimates with the LiDAR datasets. The RMSE is used to measure the accuracy of the estimated depth values; and bias is used to indicate overestimation (positive value) or underestimation (negative value):

$$\mathrm{RMSE} = \sqrt{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(x^{est}-x^{obs}\right)}^2}$$

(9)

$$\mathrm{bias} = \frac{{{\displaystyle \sum }}_{i=1}^n\left(x^{est}-x^{obs}\right)}{n}$$

(10)

where n is the number of observations, and x^est and x^obs are the estimated and measured depths, respectively. Values closer to zero for both error metrics indicate a better result. SDB obtained with Rrs_raw and Rrs_corrected are hereafter referred to as SDB_raw and SDB_corrected, respectively. These summary statistics (R², RMSE, and bias) were calculated both for the full depth range and for depths ranging from the surface down to the per-scene depth threshold (Table 4).

4. Result and Discussion

Scatterplots showing water depth estimates produced from both Rrs_raw and Rrs_corrected images and the LiDAR depth measurements are shown in Figure 4, and summary statistics (R², RMSE, and bias) are listed in Table 4.

The RMSE values for SDB_raw/SDB_corrected estimates range from 1.35/1.05 m to 2.21/2.04 m for the full depth range, and from 1.21/0.95 m to 2.19/1.89 m when applying the per-scene depth threshold. These values are broadly comparable with other SDB studies (e.g., [10,13,14,41,48]). For example, Dekker et al. [48], who compared one empirical and five physics-based approaches to bathymetry mapping using hyperspectral imagery, reported RMSE values between 0.86 (best) and 4.71 m (worst) for depths less than 13 m for two clear tropical waters in the Bahamas and eastern Australia, suggesting that our results are typical of what should be expected from a physics-based bathymetry method. It is worth mentioning here that impacts from recent hurricanes over parts of the Florida Keys, notably in 2016 and 2017, have resulted in an average increase of approximately +0.3 m in seafloor elevation over different habitat types [49,50]. Such changes were not accounted for in this study and may have had a small effect on the results, although it is worth noting that the best SDB estimate is from 2019 [Figure 4h], after the hurricanes.

Figure 4 shows that accuracy decreases with depth for both SDB_raw and SDB_corrected, particularly beyond approximately 15 m where the proportion of the measured signal originating from reflection at the seafloor becomes very small. In general, for depths shallower than 15 m, SDB_corrected points cluster more tightly around the 1:1 line that do the SDB_raw points.

4.1. Effects of Image Conditions on Depth Accuracy

4.1.1. Turbidity

Out of the nine Rrs images we applied the correction factor to, eight corrections resulted in negative RMSE changes when considering the two depth limits used in this study, with reductions ranging from 0.09 to 0.48 m (for the full depth range) and 0.07 to 0.46 m (for the per-scene depth threshold). For both depth limits, only one correction resulted in increased RMSE (i.e., the image from 01/05/2015, see Table 4) (RMSE values for SDB_raw and SDB_corrected will hereafter be referred to as RMSE_raw and RMSE_corrected, respectively). For this image, accurate depth estimates were not possible beyond approximately 14 m (Figure 4b), regardless of correction, and RMSE_corrected increased marginally by 0.01 m for the full depth range and by 0.04 m for the per-scene depth threshold (i.e., 0–11.23 m). A visual inspection of this image shows sediment plumes in the study area (Figure 5b), which suggests that turbidity contributed to an underestimation of water depth for both SDB_raw and SDB_corrected [14], and the image is of marginal use for SDB, regardless of correction. Similarly, one of the two corrected images with the lowest RMSE reduction for the full depth range also has what looks to be a silt plume emerging from nearshore channels in the southwestern portion of the area for which depth was calculated (Figure 5a). For this image, RMSE value marginally decreased by 0.09 m for the full depth range and by 0.26 m for the per-scene depth threshold (i.e., 1 to 12.03 m).

4.1.2. Glint

With RMSE_raw values of 2.21 m and 2.19 m for the full depth range and the per-scene depth threshold, respectively, the image from 13/01/2018 produces the poorest SDB_raw and SDB_corrected estimates out of the nine images. As shown in Table 4, this image has the highest RMSE and positive bias values. A visual inspection of this image (Figure 5c) indicates the presence of a moderate glint. Glint correction was not performed, as the image did not show any noticeable improvement after the initial testing. Likewise, for the image from 14/02/2018 (Figure 5d), the high negative bias values for both depth ranges, regardless of correction (Table 4), may be attributed to residual sun glint in addition to light turbidity. While an attempt was made to de-glint this image as described in Section 3.1.2, given that the NIR-based de-glinting method [36] implemented (1) relies on manual selection of deep-water pixels to estimate glint contribution, (2) assumes that there are glint-free pixels among those selected [51], and (3) assumes a homogenously low R_rs(NIR) across all water pixels, failure to meet these conditions may have resulted in the observed residual glint. For example, R_rs(NIR) may be non-negligible in glint-free but very shallow or turbid waters, or where reflective vegetation such as seagrass is close to the upper water column [52]. For these two images, the correction produces slightly reduced RMSE values (Figure 4e,f), substantially so for their respective depth thresholds (Table 4).

4.2. Effect of Wind Speed and SZA on SDB Performance

Out of the eight scenes whose SDB performance improved with the correction for the R_rs images, greater corrections were done for five scenes (i.e., Figure 4a,c–f) acquired with SZA > 49° (Table 1), and one scene (i.e., Figure 4g) acquired during high wind speed (3.21 m/s) (Table 1), when considering the per-scene depth thresholds. Likewise, when considering the full depth range, greater corrections were also done for the same number of scenes under similar environmental conditions, with the exception of the image from 26 January 2017, whose change in RMSE is −0.09 m (Table 4). It should be noted that the most noticeable RMSE reduction for SDB_corrected images for both depth limits considered in this study (i.e., Figure 4d) was observed for the scene with the highest SZA and highest wind speed. This is supporting evidence for the existence of a relationship between ACOLITE’s overestimation of R_rs in the first two bands of L8 and these two environmental variables [23], and it gives an idea of the magnitude of its impact on SDB performance. A recent study [53] also found a dependency between AC retrieval accuracy and wind speed in coastal waters. While there is strong evidence to conclude that the correction factor used in this study lowers RMSE values for images with high wind speed, it should be noted that wind speed data used in this study come from 6-h estimates in reanalysis model, and thus, they have their own uncertainty. For this reason, more testing may be needed for a firmer conclusion about the relationship between ACOLITE’s error and wind speed.

4.3. Bathymetry Estimates at Different Depth Ranges

Figure 6 shows the performance of SDB_raw and SDB_corrected for each image, binned to 5 m depth increments for the 1–25 m depth range. The accuracy of SDB estimates decreases with increasing depth for both SDB_raw and SDB_corrected. While higher RMSE values should be expected at deeper depths due to the diminishing signal from seafloor reflectance, it should be noted that the number of LiDAR points for validation is also smaller at deeper depths, leading to increased uncertainty around the RMSE values reported at these depths.

4.4. SDB Estimates Using an Ensemble Approach

Most SDB studies are based on a single image for each study area, with researchers typically selecting the best available image using visual inspection [14]. Our results indicate that this may not be a robust approach. To illustrate the problem, we invite readers to visually inspect the nine images used in this study (Figure 2) and identify the one that looks most suitable for SDB. Then, proceed to Table 4 to see if it was indeed the one that produced the best results, as measured by RMSE, bias, or R². An informal test among our colleagues, all of whom work on OC remote sensing, suggests that it is not easy to identify the best scene. However, a unique advantage of optical remote sensing is the repetitive acquisition of images over the same area. This is especially important for SDB, where the suitability of a given image is determined by transient environmental factors, such as cloud and aerosols, sea surface state, and turbidity [54]. We explored one way of taking advantage of the multiple images available for the study area by testing an ensemble approach in which we calculated the per-pixel median depth value of all nine corrected images (i.e., SDB_corrected) used in this study. Then, we compared the resulting depth estimates with those obtained using the best individual image from the analysis in Section 4.1 (i.e., the image from 01/02/2019, see Table 4). Figure 7 shows that the results produced by the ensemble are very similar to those obtained with the best individual image. SDB estimates up to approximately 15 m are similar to the SDB_corrected estimates from the 01/02/2019 image, as are the RMSE values for the 1–5, 5–10, and 10–15 m depth ranges (Figure 8). Outliers are noticeably reduced in the ensemble result when compared to any of the sub-optimal images that were also included in its calculation, suggesting that the use of median depth is effective in eliminating noise in the ensemble. Thus, the ensemble approach eliminates the need for the selection of a single best image, while producing SDB results of similar accuracy. In this context, it is noteworthy that one of the two best images (i.e., image from 05/03/2019) is also the one with the lowest NE∆R_rs in bands 1, 2, and 4, as well as the second-lowest NE∆R_rs in band 3 (Table 2). This suggests that one effective way to pre-screen images, either for a single best image approach or to determine which images should be included in an ensemble, could be to estimate NE∆R_rs and select those images with the lowest values across the visible bands.

5. Conclusions

In this study, we demonstrated the use of Landsat 8 data for physics-based SDB in US coastal waters. A state-of-the-art AC method (ACOLITE) was used to convert per-pixel radiometric units to R_rs, and an RTM was inverted to estimate water depth, which was compared to airborne LiDAR validation data. The results showed that ACOLITE can be used to produce SDB from imagery that is free of conditions such as clouds, glint, sediments plumes, boats and wakes, with an accuracy (RMSE 1.34 m for 1–25 m depth range and 1.21 m for approximately 1–15 m depth range) comparable to that reported from empirical and physics-based SDB elsewhere. To account for ACOLITE’s known overestimation of R_rs for Landsat 8′s coastal and blue bands, we applied a correction factor, which was calculated as a function of solar zenith angle, aerosol optical thickness, and wind speed, to obtain a corrected set of R_rs images. This correction further improved bathymetry estimates for eight of the nine scenes we tested, with the resulting changes in bathymetry RMSE ranging from +0.01 m (worse) to −0.48 m (better) for a 1 to 25 m depth range, and from + 0.04 m (worse) to −0.46 m (better) for an approximately 1–16 m depth range. Using a total of nine Landsat 8 images, we showed that the correction factor improved SDB results, both on average (ΔRMSE = −0.22 m) and for the best single image (ΔRMSE = −0.30 m) for the 1–25 m depth range. SDB improvements from application of the correction factor were the greatest for images acquired at a high solar zenith angle and at high wind speeds, where ACOLITE is known to have the greatest bias. The correction method demonstrated in this study can be implemented with any appropriate AC algorithm. Finally, we demonstrated that an ensemble approach based on multiple images, with acquisitions ranging from optimal to sub-optimal conditions, can be used to derive bathymetry with a result that is similar to what can be obtained from the best individual image. This is important because it is rarely visually obvious which of several images is best for SDB, and the ensemble approach can be automated to reduce time spent on pre-screening and filtering scenes, and it can potentially also reduce the amount of missing pixels caused by clouds and cloud shadows encountered in any single image. Automating SDB will ultimately facilitate the efficient and operational use of the globally available L8 (and other multispectral) datasets.