Get Our e-AlertsSubmit Manuscript
Journal of Remote Sensing / 2022 / Article

Research Article | Open Access

Volume 2022 |Article ID 9831947 |

Wendian Lai, Zhongping Lee, Junwei Wang, Yongchao Wang, Rodrigo Garcia, Huaguo Zhang, "A Portable Algorithm to Retrieve Bottom Depth of Optically Shallow Waters from Top-Of-Atmosphere Measurements", Journal of Remote Sensing, vol. 2022, Article ID 9831947, 16 pages, 2022.

A Portable Algorithm to Retrieve Bottom Depth of Optically Shallow Waters from Top-Of-Atmosphere Measurements

Received04 Sep 2021
Accepted15 Dec 2021
Published03 Feb 2022


Bottom depth () of optically shallow waters can be retrieved from multiband imagery, where remote sensing reflectance () are commonly used as the input. Because of the difficulties of removing the atmospheric effects in coastal areas, quite often, there are no valid from satellites for the retrieval of . More importantly, the empirical algorithms for are hardly portable to new measurements. In this study, using data from Landsat-8 and ICESat-2 as examples, we present an approach to retrieve directly from the top-of-atmosphere (TOA) data. It not only bypasses the requirement to correct the effects of aerosols but also shows promising portability to areas not included in algorithm development. Specifically, we use Rayleigh-corrected TOA reflectance () in the 443–2300 nm range as input, along with a multilayer perceptron (), for the retrieval of . More than 78,000 matchup points of and (0–25 m) were used to train , which resulted in a Mean Absolute Percentage Difference (MARD) of 8.8% and a coefficient of determination () of 0.96. This was further applied to Landsat-8 data of six regions not included in the training phase, generating MARD and values of 8.3% and 0.98, respectively. In contrast, a conventional two-band ratio algorithm with as the input generated MARD and values of 31.6% and 0.68 and significantly fewer retrievals due to failures in atmospheric correction. These results indicate a breakthrough of algorithm portability of in sensing of optically shallow waters.

1. Introduction

Nearshore shallow water environments such as coral reefs, seagrass, and kelp beds are among the most socioeconomically important and productive ecosystems in the world; its monitoring is an important task of many government agencies [1, 2]. In addition to monitoring changes across bottom substrates in such ecosystems, one desired parameter is the bottom depth, as it is important not only for navigation but also for studies of coastal processes and management of coastal events ranging from monitoring of storm surge to site selection of wind farms [36].

Conventionally, the measurement of bottom depth (or bathymetry) is carried out using multibeam echosounder. In more recent years, with the advancement of laser technology, more measurements of bottom depth (, units: m) are taken by airborne or spaceborne Lidar. While these methods and systems provide a high precision product [7, 8], they are high-cost, time-consuming, and limited to the areas of vessels or airplanes can reach, with the resulted data from space hardly to form a high-resolution bathymetric map [9].

With the operation of remote sensing imagers (or passive remote sensing) and that an optically shallow bottom can alter the spectral signature of the radiance emerging from the water surface, many radiometric data have been applied to estimate [2, 4, 1012]. Based on the imaging characteristics of such measurements, bathymetry can be derived from a variety of satellite imaging sensors that have a range of spatial resolutions [1214]. Although its accuracy and precision cannot match that from Lidar, the imaging capability offers an efficient and cost-effective way to obtain an map of large spatial regions, with its spatial resolution dictated by the sensors.

Passive sensors, especially those onboard satellite platforms, capture, after calibration, the total radiance () defined as the sum of radiance contributions from the atmosphere, air-sea interface, water column, and bottom substrate. Sophisticated algorithms are thus required to derive from . To this end, numerous algorithms have been developed since the 1970s for the estimation of from above-water radiometry [4, 12, 1517]. Generally, these algorithms can be divided into empirical and semianalytical approaches, with empirical algorithms purely data driven, while semianalytical algorithms are based on radiative transfer theory. Although semianalytical algorithms [1719] do not require known depths for algorithm development, they require adequate spectral bands in order to separate water-column and bottom information [20], a prerequisite met by few current satellite image sensors.

With the recent launch and operation of the Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2), high-precision bottom depth ( hereafter) of optically shallow waters (OSW) can be well produced from data obtained by the Advanced Topographic Laser Altimeter System (ATLAS). In the meantime, many satellites with good spatial-resolution imagers, such as Landsat-8 (30 m resolution), Sentinel-2 (20 m resolution), and Coastal Zone Imager (50 m resolution), have been taking routine measurements of coastal regions, thus it is logical to fuse both and imager data to generate maps of bottom depth through empirical regressions [2123]. One typical example is the Two-Band Ratio Algorithm (TBRA) developed about two decades ago [2] that uses the ratio of remote sensing reflectance () at two bands for the estimation of , with the empirical coefficients of TBRA derived for each image [3, 24]. However, the use of two spectral bands runs into nonuniqueness issues that results in large uncertainties particularly when bottom substrates and/or water properties change significantly within a scene [25].

With the advancement in computer technology, neural network (NN) is a powerful approach for accurate bathymetry retrievals. Compared with traditional band ratio schemes, NN is flexible in taking a wide range of inputs. At present, many NN algorithms with as the input have been developed for the estimation of H [2629]. However, for nearshore shallow regions, such an approach is sensitive to atmospheric correction algorithms [30]. When there are no valid due to failures in atmospheric correction [31, 32], there will be no product for those locations.

In this study, we take advantage of the availability of both from ICESat-2 and imager data (using OLI of Landsat-8 as a demonstration), combined with NN (more specifically, multilayer perceptron (MLP)), to show a promising scheme to obtain map. In particular, to avoid issues related to failures in atmosphere correction, we used top-of-atmosphere data, similar to Kutser et al. [1], for the retrieval of . Rather than creating a specifically tailored hyperspectral database [1], where its applicability is likely limited to a unique environment, we used Rayleigh corrected top-of-atmosphere reflectance () as the input, along with a MLP for the estimation of . More importantly, we further applied the trained NN algorithm to Landsat-8 OLI images of different regions not included in the training phase and obtained highly consistent compared with , which highlights a breakthrough in portability compared to the traditional empirical algorithms. With the increased data and imager measurements, we foresee the generation of accurate, high-resolution, bathymetry map of global OSW waters.

2. Data

2.1. Landsat-8 Satellite Image

While many multiband imagers can be used for the empirical retrieval of , we used Landsat-8 OLI data due to its following features: (1)A spatial resolution of 30 m provides detailed features of coastal regions where bottom depth may be highly heterogeneous(2)Seven spectral bands cover the 430–2300 nm range, which contains spectral information of water optical properties, depth, and bottom substrate (visible domain) as well as atmospheric aerosols (NIR-SWIR)(3)Top-of-atmosphere radiance from Landsat-8 OLI is well calibrated [33, 34], thus adequate for the calculation of Rayleigh-corrected reflectance ()

In addition, Landsat-8 OLI data can be processed using the SeaDAS package (v7.5.3) [35], where [31, 32] of OLI measurements can be adequately generated. Here, is defined as [3638] where is the total reflectance at top-of-atmosphere, defined as . Here is the extraterrestrial solar irradiance [39], the solar zenith angle, the gas transmittance, and the reflectance from Rayleigh scattering. Both and can be calculated with SeaDAS using satellite ancillary information.

The SeaDAS cloud albedo value (default: 0.018) used to mask pixels with cloud contamination was increased to 0.03 to avoid masking bright shallow water pixels at the consequence of including pixel contaminated with thin clouds. with a value less than 0.0 or greater than 1.0 at any band was also masked. To compare the performance of the traditional TBRA scheme, images of of these Landsat-8 OLI measurements were also generated using SeaDAS. For these products, pixels designated as low-quality were removed based on the standard Level-2 quality flags included in SeaDAS, such as atmospheric correction failure (ATMFAIL), land pixel (LAND), probable cloud or ice contamination (CLDICE), very high or saturated observed radiance (HILT), and strong sun glint contamination (HIGLINT).

2.2. ICESat-2 Data

ICESat-2 was launched on September 15, 2018, and carries ATLAS, a photon-counting laser altimeter. The green (532 nm) laser pulse from ATLAS illuminates three beam pairs on the surface, with the distance between beam pairs about 3 km in the across-track direction [40]. In each beam pair, the ratio of the transmit energy of the beams is approximately 1 : 4, which are called weak and strong beams, respectively [40]. The laser pulse of each beam along the ground track is 0.7 m. The photon heights of each beam pair are very similar, and in order to reduce the degree of data redundancy, the ICESat-2 global geolocated data ATL03 from the strong beam (or week beam in the case of a lot of noise in the strong beam data) was used in this study, which contains photon height, time, geodetic latitude and longitude, and Geoid height for individual photons above the WGS84 ellipsoid ( [40].

The procedure to obtain the bottom depth from ICESat-2 photons is as follows [11]: (1)Invalid photon removal: within the ATL03 dataset, a quality parameter of “confidence” (with a range of -2 to 4) is provided to classify each photon as likely signal or noise. Here, as other studies [21, 40], we consider the pixels having a “confidence” value from 0 to 4 to be valid.(2)Geoid system conversion: the photon height reported in the ATL03 dataset is WGS-84 ellipsoidal height (); thus, each track of raw photons on the water surface is not parallel to Geoid. As such, the photon height was transformed from the WGS-84 system to the Geoid system, , where is the Geoid height (distance between Geoid model and ellipsoid), and is available within ATL03.(3)Outlier removal: the raw ATL03 photon data are contaminated with a significant amount of noise (see Figure 1(a)). In this work, a density-based spatial clustering of applications with noise (DBSCAN) method [41] was employed, which groups points that are closely packed together (points with many nearby neighbors) and marks outliers that lie alone in low-density regions (whose nearest neighbors are too far away) [21, 41]. Similarly, as in Babbel et al. [11], we used the occurrence of photons in a rectangular window (length is the along-track distance and width is the elevation) to remove outliers. Basically, for any photon at the center of this rectangle, if the number of total photons within the rectangle is below a threshold, it is treated as noise. In this study, we found that a window size of 10 m (distance along the track) by 1 m (elevation) with a threshold of 7 photons is suitable, where photons from sea surface and those from sea bottom are clearly separated (see Figure 1(b)). Sliding the rectangular window along the track then generates a line of apparent bottom depth.(4)Calculating the geometric depth of bottom photons: although the Geoid system conversion is applied, it can be clearly seen that photons reflected from the water’s surface have a slope in the along-track direction due to tide or other hydrodynamic processes (Figure 1(b)). According to the photon geolocation algorithm for ATL03 of ICESat-2, photons that are characterized as likely signals are binned in approximately 20 m along-track segments [40]. Within each segment, the water depth is calculated as the difference between the height of bottom-reflected photons and the mean height of surface-reflected photons.(5)Correcting the impact of light speed: photons from ICESat-2 propagate from air to water and then water to air. Thus, there is a change of light speed between the two media that results in a systematic bias in the calculated bottom depth. The light-correction algorithm of Parrish et al. [42] was applied to remove this effect from the bottom depth obtained in #4 above, where the uncertainty of the corrected depth is extremely low (a few centimeters) [21, 42, 43].

2.3. Data Matchup between Landsat-8 and ICESat-2

Assuming the variation of bottom depth (after tidal correction) is negligible within a short period, a time constraint of ±2 weeks between Landsat-8 OLI and ICESat-2 measurements was used to matchup radiometric and Lidar data [25]. Meanwhile, the corrected bottom depth of ICESat-2 was adjusted to match the tidal cycle of the Landsat-8 overpass, where the classical tidal harmonic analysis model [44] was applied to calculate tide information of targeted locations.

To match the measurements between ICESat-2 and OLI, we first located the ICESat-2 track within the OLI image. For a point of ICESat-2 corrected bottom depth, an OLI pixel was first selected based on the closest Euclidean distance. Since the footprint of ICESat-2 along the ground track is 0.7 m, while the spatial resolution of OLI is 30 m, there are many ICESat-2 points within a Landsat-8 pixel. The weighted average of ICESat-2 pixels located within an OLI pixel was calculated following Lee et al. [45], and this average (represented as hereafter) is considered the “true” bottom depth corresponding to an OLI pixel. This is calculated as where is the number of ICESat-2 points within the OLI pixel, is the ICESat-2 bottom depth from Section 2.2, and and are the tide heights at the time of the ICESat-2 and OLI acquisitions, obtained from the classical tidal harmonic analysis model [44].

2.4. Training Dataset
2.4.1. Depth Dataset

A reliable NN system requires a large training dataset with a wide range of water optical and bottom properties. In this study, the ICESat-2 and OLI data over the Great Bahama Bank (GBB) and the Cay Sal Bank (CSB) (yellow boxes in Figure 2) were selected to train the NN algorithm. The Landsat-8 acquisition dates are listed in Table 1. To ensure enough representation, a cloud coverage of less than 20% was applied for the CSB data, while this cloud threshold is 10% for GBB. A total of 72,178 pairs of and for GBB and 10,728 pairs for CSB were obtained following the steps described above. In contrast, and due to failure in atmosphere correction, the matchup between and was 30,549 pairs for GBB and 5,499 pairs for CSB, which are about 50% of that between and .



The bottom depth in the GBB and CSB area ranges ~0–15 m and ~5–25 m, respectively (Figure 2). The bottom substrates of GBB is spatially and temporally variable, which includes seagrass, microalgae, brown macroalgae and sand [46], while the major benthic substrates of CSB are seagrass and sand [47]. Figure 2 also includes the histogram of from OLI, which is in a range of ~0.001-0.2, with most of the data having (equivalent aerosol optical density at this wavelength is ~0.4), indicating a wide range of aerosol loadings.

2.4.2. Dataset for Optically Deep and Optically Shallow Waters

We classified water types into optically deep (ODW) and optically shallow (OSW) waters, where a separate NN was developed first for this segmentation, as empirical algorithms (explicit empirical functions or implicit neural networks) for can only be applied to OSW. The data used for the NN training came from known ODW (Landsat-8 OLI measurements in Massachusetts Bay (2019/04/17), Chesapeake Bay (2019/11/25), Yangtze River (2020/12/22), and Tongue of The Ocean (TOTO) (2018/10/14)) and OSW environments (Figure 3). OSW is determined by having valid from ICESat-2, while ODW is based on the literature and visual inspection. From these OLI acquisitions, 80,000 spectra were randomly selected from the ODW regions, labeled as “ODW,” and 80,000 groups of were selected from the OSW regions and labeled as “OSW.” Figure 4 shows examples of spectra of optically deep and optically shallow waters.

2.5. Evaluation Dataset

We extend the evaluation of the NN performance with data from six different regions across the Caribbean and South China Sea, covering a wide range of aerosol loadings, bottom depths, benthic types, and water properties that were not included in the training dataset. Figure 5 presents true color images and ICESat-2 tracks of these six regions that include the Little Bahama Bank (LBB), Dry Tortugas (DT), the Bight of Acklins (BOA), the Caicos Bank (CB), Xisha Island (XI), and Dongsha Island (DI) (Figures 5(a)–5(f)). The acquisition time of Landsat-8 and ICESat-2 data used in this effort is included in Table 2. These regions represent a variety of geomorphic zones (reef crest, patch reef, lagoon, and open ocean) and benthic types (coral reef, seagrass, and sand), which provide a completely independent evaluation of the algorithms developed.

AreaLat/LonLandsat-8 dateICESat-2 dateMatchup points

LBB (a)26.44N/-77.48E2019/07/132019/07/181064
DT (b)24.66N/-82.85E2019/09/042019/09/2446
BOA (c)22.53N/-74.09E2020/09/282020/09/28289
CB (d)21.59N/-71.96E2020/12/102020/12/161317
XI (e)16.23N/111.69E2019/04/202019/04/2111
DI (f)20.69N/116.81E2019/05/172019/10/2170

3. Methods

3.1. Neural Network Systems

An image processing workflow consisting of two MLP systems was developed for deriving from . First, a pixel is classified as OSW or ODW using , where the depth of pixels classed as OSW is derived with . We termed this workflow as neural-network-based-TOA-algorithm-for-bathymetry, NNTOA-B. Figures 6(a) and 6(b) present schematic charts of the two MLP systems. Fundamentally, MLP is a forward-structured artificial neural network, which is capable of handling nonlinear separable problems, easy to be implemented, scalable, fast, high stability, and efficient in the learning process [48].

Three units constitute an MLP architecture: an input layer, several hidden layers, and an output layer. The number of neurons in the input layer is determined by the input data (the number of spectral bands). The MLP networks for both and have the same input layer containing seven neurons (the total number of OLI bands). The number of hidden layers and the number of neurons in each layer need to be determined for optimal results. After trial and error, two hidden layers with 32 and 16 neurons (Figure 6(a)) are found to work the best for , while three hidden layers with 128, 32, and 16 neurons (Figure 6(b)) were optimal for . The output of each neuron is converted to the input for the neurons of the next layer through a nonlinear activation function. The Rectified Linear Unit (ReLU) was employed as the activation function of the hidden layers [49], which largely avoids gradient disappearance. The number of neurons in the output layer is determined by the desired output (the depth), which is one neuron for both and . In this layer, the sigmoid activation function was used for binary classification of water types [50], while the linear activation function was used for water depth prediction.

Backpropagation, short for “error backpropagation,” is a common procedure used in conjunction with optimization methods to train NN [51]. This procedure computes the gradient of the loss function (mean squared error for and accuracy for ) for all the reweightings in the network, and this gradient is fed back to the optimization method to update the weights and biases to minimize the loss function using the learning rate () as a step. The adaptive moment estimation (Adam) of minibatch gradient descent, with regular hyperparameters as , was used for optimization calculations [52]. The learning rate of 0.001 and batch size of 128 performed well in this process. The training phase is finally achieved when the iterations stop (2,000 epochs in our case) and the loss function converges.

The Keras library, which is a powerful and easy-to-use interface of Tensorflow [53], was used to train the two MLPs. Tensorflow is an open-source library for numerical computation [53], and its interface Keras provides a Python interface to NN. Creating NN models on Keras is quite straightforward because of its emphasis on implementing user-friendliness [54].

3.2. Accuracy Assessment

To verify whether the trained model is over- or underfitted, we subdivided the (GBB and CSB) training dataset (Section 2.4.1) into subsets of 95/5 training/testing. Changing from 95/5 to 80/20 did not impact the results.

To measure the deviation or difference of the estimated from (represented as hereafter), we adopt the coefficient of determination () between any two sets of data, the Root Mean Square Difference (RMSD) and the Mean Absolute Relative Difference (MARD), between and as metrics. They are defined as follows: where is the total number of pairs used in the analyses.

4. Results

4.1. Optically Deep vs. Optically Shallow Waters

First, we evaluated the classification accuracy using the test subset. An accuracy of 100% was reached, where no ODW was labeled as OSW, and vice versa. To verify its applicability to a wider range of environments, we applied to the six independent images from different regions, with classified ODW and OSW presented in Figure 7. Visually, the OSW pixels in these regions, which have distinctively different colors (see Figure 5) and spectra (see Figure 4), are separated from the ODW pixels. Further, using pixels having as an independent reference of OSW, the accuracy of correctly identifying OSW by is in general greater than 94% for these images (see Table 3), as a few OSW pixels were classified as ODW due to low signal from the bottom or complex topography in the Landsat-8 OLI images. Nevertheless, these evaluations indicate that is effective in separating OSW from ODW, although future upgrades and revisions are always necessary in algorithm development.

AreaAccuracy (%)

LBB (a)100.0
DT (b)97.8
BOA (c)100.0
CB (d)94.3
XI (e)100.0
DI (f)98.6

4.2. Bottom Depth
4.2.1. Training Data

Figure 8(a) shows plotted against of the 78,698 subset for training . For these multitemporal images with varying water properties, aerosols, and bottom characteristics, NNTOA-B achieved an and MARD of 0.96 and 8.8%, respectively, for ranging between ~0 and 25 m, suggesting high predictability of NNTOA-B. There are a few (<0.1% of the total data) outliers in , which are likely due to imprecise matchup between the footprints of the two satellite sensors and the sharp drop-offs of near cliffs (see Section 4.2.3), as well as photon outliers that were not properly filtered. These are not avoidable when measurements from two satellites of different characteristics are compared with each other.

Further, NNTOA-B was evaluated using the test subset (4,141 points), i.e., 5% of the training data, with compared with showing in Figure 8(b). Like the training subset, most of the data points are around the 1 : 1 line, with as 0.92 and MARD as 8.6% for in the range of 0–25 m. If the few (<0.1%) outliers are removed, the and MARD become 0.98 and 8.1%. These results further support the capability of to estimate using (430–2300 nm).

4.2.2. Bottom Depth of GBB and CSB

Using the trained and , we processed eight Landsat-8 OLI images in the GBB and CSB and compiled a shallow water bathymetry map from these OLI measurements (see Figure 9(a)). For each image, the water depth () is the result of minus the tidal height at the observation time; i.e., they all are normalized to the mean sea level. If one pixel has more than one product from OLI, an averaged was calculated from the available products. This bathymetry map, a composite from eight OLI images, provides a complete and high-resolution (30 m) detail of the GBB area never achieved before. Compared with the depth (300 m resolution) generated from 187 MERIS images over the GBB [46], the bathymetry patterns of both maps are very consistent ( is 0.89 between the two images), along with a mean absolute relative difference of 14.7% (see Figure 9(b)). In view that the MERIS product did not use any ICESat-2 information and that different data sources with different lifespan and fundamentally different algorithms used for the generation of the maps, this comparison indicates highly consistent products from the two independent determinations.

4.2.3. Depth of Independent Dataset

The key to any remote sensing algorithm is its applicability to new measurements or new environments that are not included in the training phase of algorithm development. We evaluated the performance of NNTOA-B with the six independent OLI images obtained from regions not covered by the training data. The resulted bottom depths are shown in Figure 10, while Figure 11 shows scatterplots between and . Across all six locations, agrees with excellently (Figure 11). Overall, the value is 0.98, along with RMSD of 0.5 m and MARD of 8.3% for these data. Since these image data were not included in the training, these results show a portability achieved by NNTOA-B that was missing in previous empirical algorithms for bathymetry, at least for the data included in this effort.

Among the six locations, better statistics were found for LBB ( ranging between 0.5 and 9 m, Figure 11(a)) and CB ( ranging between 0.5 and 20 m, Figure 11(d)). In part, this could be a result of much more data points, so a few outliers due to matchup issues have less impact on the statistical evaluation; or the fact that these regions are closer to the data used to train . For BOA ( ranging from 0.5 to 3 m, Figure 11(c)), the estimated was systematically deeper by ~0.2 m. This is likely due to the nearly closed topography of this region, where an accurate tide correction is difficult [55].

For the other three regions (DT, XI, and DI), the statistics between and are not as good as that of the three regions discussed above. In part, it is because of the fewer number of matchups for these regions: 46 for DT (Figure 11(b)), 11 for XI (Figure 11(e)), and 70 for DI (Figure 11(f)), a result of small terrain and high cloud cover in the OLI images. Therefore, a few outliers in the matchup data can impact the statistics significantly. In addition, because of the small terrain and subpixel variations (e.g., XI; see Figure 10(e)), it is difficult to achieve precise matchup between Landsat-8 and ICESat-2 measurements, where the footprint of OLI is 30 m. Further, the spatial variation in DT area is quite complex [56] (also see Figure 10(b)); thus, spatial heterogeneity will also contribute to the slightly degraded statistics for these regions. Nevertheless, although the training of is completely independent of these regions and environments, the overall and MARD for these three regions are still 0.85 and 19.0%, indicating very good estimation of on average by .

For XI and DI, located in the Eastern Hemisphere, there is a similar trend of underestimation for depths greater than 7 m. The reason being that the water depth changes rapidly within an OLI pixel as it is at the edge of the reefs (see Figure 12 for an example), making it extremely difficult to obtain a representative bottom depth corresponding to the OLI spectrum of such a pixel, even though the weighted average of water depth was used [45]. Because shallower bottom contributes more to the measured total signal within a pixel, remotely sensed bottom depth from water color will in general tilt to an underestimation if there are strong depth variations within a pixel [45].

4.3. Performance of a Band-Ratio Algorithm

The empirical two-band ratio algorithm (TBRA) using as the input for the estimation of is often applied to obtain maps of bottom depth from multiband imagery [2, 21, 57]. As first proposed by Stumpf et al. [2], TBRA takes the form of

Here, and for Landsat-8 OLI are 482 and 561 nm; can be fixed at 1000, while values of and are empirical coefficients to be derived from a training dataset [2, 10, 21, 57].

For the matchup points used to train , we also compiled a dataset containing and , but the total numbers are much less (36,048 matchups) due to issues with atmosphere correction in shallow coastal waters. Also, 95% of this dataset was used to obtain the values of ; the other 5% used to verify the empirical algorithm. Figure 13(a) compares with of the model training subset, while Figure 13(b) compares with of the model test subset. For these datasets, TBRA obtained , RMSD and MARD of 0.75, 1.5 m and 28% for the training subset and 0.73, 1.5 m, and 30% for the test subset, respectively, which are significantly worse compared to the performance of for these measurements (see Figure 11 and Table 4).

MethodDatasetRMSD (m)MARD (%)Data number

Training subset0.960.78.878,698
Test subset0.920.98.64,141
TBRATraining subset0.751.528.034,246
Test subset0.731.530.01,802

As the evaluation of NNTOA-B to independent data, we also applied this newly updated TBRA to the six regions independent of algorithm development (see Figure 14). Despite the accurate retrievals for CB (see Figure 14(d)), the other five regions show quite large uncertainties, which include LBB where the environment is quite similar with that of GBB and CSB used to update the coefficients of TBRA. These tests and evaluations highlight, as indicated in previous studies [2, 25], that the algorithm coefficients () of TBRA are highly region and data dependent. Thus, for each scene, a new set of algorithm coefficients are required for more reliable derivation. Furthermore, because it uses as the input, significantly less coverage in the products is obtained compared to using (430–2300 nm) as the input. These results emphasize the significantly improved applicability of for the remote sensing of bottom depth. In addition, as pointed out in Lee et al. [25], the application of TBRA requires a step to separate OSW from ODW; otherwise, erroneous depths can be obtained.

5. Discussion

The above tests and evaluations indicate that NNTOA-B can obtain highly accurate estimation of from OLI-derived (430–2300 nm) and that NNTOA-B has high portability, implying that reliable could be obtained from new OLI acquisitions without modifying NNTOA-B. One of the likely reasons for the promising results lies in the simultaneous use of the visible-SWIR domain where traditionally two-separate routines have been used in the data processing of satellite ocean color measurements; the first being atmospheric correction and the second the retrieval of water and/or bottom properties. Atmospheric correction focuses on the spectral information in the NIR-SWIR range, but the second component (for properties below the sea surface) focuses on the spectral information in the visible range, although all come from the same or spectrum. NNTOA-B does not explicitly separate the two processes for the retrieval of from , but since all spectral information in the visible-SWIR is employed in this algorithm, the information related to aerosol is included (where is as high as 0.2, see Figure 2), although no aerosol product is generated with . It thus avoided the challenges of obtaining highly accurate for nearshore shallow waters [30]. It is important to emphasize that the contributions from Rayleigh scattering can be accurately calculated [58] and removed. Since Rayleigh contributions make a significant part of the top-of-atmosphere radiance in the blue-green domain, which is also the spectral window that solar radiation best reaches the bottom for most natural waters, a spectrum has significantly enhanced water signal (see Figure 4) than a spectrum where contributions of Rayleigh scattering are still included.

The excellent results of retrieval by NNTOA-B from OLI measurements are somewhat surprising, as Landsat-8 OLI has four bands (with wide bandwidth) in the visible (400–700 nm) domain, which contradicts the conclusion of an earlier study [20] that it is better to have 15 or so bands in the visible domain for shallow water remote sensing. This difference lies in the approaches used, where Lee and Carder [20] used a semianalytical algorithm (HOPE) to simultaneously solve for unknowns of water and bottom properties, a scheme does not assume any ranges or values of the targeted parameters, while relies on “learning” of provided pairs of and ; thus, the solution of from is in a predetermined parameter space, rather than the widely open space that HOPE-based semianalytical algorithms take. Fundamentally, an spectrum of OSW is a function of four variables [18]: absorption and backscattering coefficients, bottom reflectance, and . With limited to a predetermined space from the learning processing, four bands in the visible domain of OLI imagery appear to have sufficient information for accurate retrieval of by a NN, although the accuracy may not necessarily be the same for the retrieval of the optical properties such as backscattering coefficient or bottom reflectance.

It is necessary to point out that, although the results are greatly promising, this study focused on OSW in the tropical and subtropical regions, where water properties are generally clear, and the range of (865) is 0–0.2, along with marine aerosols. It remains to be seen how the present NNTOA-B performs for OSW in mid to high latitude waters, where both water and atmospheric properties could be very different from those encountered here. With the availability of more and Landsat-8 OLI images of various latitudes, the evaluation of NNTOA-B will be expanded; and, as all algorithm development, the present NNTOA-B will likely be revised and updated to cover mid-to-high latitude regions. In addition, as the focus of this effort is on the development and evaluation of , the confidence at pixel level [25] of NNTOA-B estimated was not calculated. This is also because that due to the significantly expanded data pool for training (>80,000), it is computationally demanding to obtain a confidence score following the scheme developed in Lee et al. [25], where improvements in its computation efficiency are beyond the scope of this study.

6. Conclusions

In this effort, with collocated measurements from both Landsat-8 OLI and ICESat-2, we show that bottom depth of OSW can be retrieved with high accuracy from the Rayleigh-corrected top-of-atmosphere reflectance spectrum after a NN-based workflow (NNTOA-B) was successfully trained with a large number and a wide range of data. This NNTOA-B shows promising portability to other shallow regions not included in the training, thus overcoming a long-standing limitation of traditional empirical algorithms for that works best with data/image used in algorithm training. With data from ICESat-2 and high-resolution imagery growing by day, we see the generation of accurate, high-resolution, bathymetric product of global OSW in the horizon. Such data product will not only help economic activities and monitoring of coastal shallow environments but also provide valuable input for the estimation of both water column properties and substrate types from multiband imagery. Remote sensing of OSW is now entering a new era.

Data Availability

The Landsat-8 OLI data are publicly available and can be downloaded from The ICESat-2 data are also publicly available through the National Snow and Ice Data Center (NSIDC). The geolocated photon data (ATL03) are found online (, with descriptions can be found in the cited reference of Neumann et al. [40].

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Authors’ Contributions

W Lai conducted the study and drafted the manuscript; Z Lee conceptualized the study and finalized the manuscript; J Wang helped processing Landsat-8 and ICESat-2 data; Y Wang helped neural networks and data analyzing; R Garcia helped in analyzing results and revising the manuscript; H Zhang contributed to the concept and helped in revising the manuscript.


Financial supports from the Chinese Ministry of Science and Technology through the National Key Research and Development Program of China (#2016YFC1400905) and the National Natural Science Foundation of China (#41941008, #41890803, and #41830102), the Joint Polar Satellite System (JPSS) funding for the NOAA ocean color calibration and validation (Cal/Val) project, and the University of Massachusetts Boston are greatly appreciated. The authors would like to thank USGS for the distribution of Landsat-8 OLI L1 imagery, NASA’s Ocean Biology Processing Group (OBPG) for the SeaDAS processing software, and the NASA ICESat-2 team for providing the data used in this study.


  1. T. Kutser, I. Miller, and D. L. Jupp, “Mapping coral reef benthic substrates using hyperspectral space-borne images and spectral libraries,” Estuarine, Coastal and Shelf Science, vol. 70, no. 3, pp. 449–460, 2006. View at: Google Scholar
  2. R. P. Stumpf, K. Holderied, and M. Sinclair, “Determination of water depth with high-resolution satellite imagery over variable bottom types,” Limnology and Oceanography, vol. 48, no. 1, pp. 547–556, 2003. View at: Google Scholar
  3. T. Kutser, J. Hedley, C. Giardino, C. Roelfsema, and V. E. Brando, “Remote sensing of shallow waters–a 50 year retrospective and future directions,” Remote Sensing of Environment, vol. 240, article 111619, 2020. View at: Publisher Site | Google Scholar
  4. D. R. Lyzenga, N. P. Malinas, and F. J. Tanis, “Multispectral bathymetry using a simple physically based algorithm,” IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 8, pp. 2251–2259, 2006. View at: Publisher Site | Google Scholar
  5. M. F. Karim and N. Mimura, “Impacts of climate change and sea-level rise on cyclonic storm surge floods in Bangladesh,” Global Environmental Change, vol. 18, no. 3, pp. 490–500, 2008. View at: Publisher Site | Google Scholar
  6. B. B. Barnes, R. Garcia, C. Hu, and Z. Lee, “Multi-band spectral matching inversion algorithm to derive water column properties in optically shallow waters: an optimization of parameterization,” Remote Sensing of Environment, vol. 204, pp. 424–438, 2018. View at: Publisher Site | Google Scholar
  7. J. L. Irish and W. J. Lillycrop, “Scanning laser mapping of the coastal zone: the SHOALS system,” ISPRS Journal of Photogrammetry and Remote Sensing, vol. 54, no. 2, pp. 123–129, 1999. View at: Publisher Site | Google Scholar
  8. M. Doneus, I. Miholjek, G. Mandlburger et al., “Airborne laser bathymetry for documentation of submerged archaeological sites in shallow water,” in Underwater 3D Recording and Modeling (ISPRS TC V, CIPA), Piano di Sorrento, Italy, 2015. View at: Google Scholar
  9. J. L. Irish and T. E. White, “Coastal engineering applications of high-resolution Lidar bathymetry,” Coastal Engineering, vol. 35, no. 1-2, pp. 47–71, 1998. View at: Publisher Site | Google Scholar
  10. J. Li, D. E. Knapp, S. R. Schill et al., “Adaptive bathymetry estimation for shallow coastal waters using Planet Dove satellites,” Remote Sensing of Environment, vol. 232, article 111302, 2019. View at: Publisher Site | Google Scholar
  11. B. J. Babbel, C. E. Parrish, and L. A. Magruder, “ICESat-2 elevation retrievals in support of satellite-derived bathymetry for global science applications,” Geophysical Research Letters, vol. 48, no. 5, 2021. View at: Publisher Site | Google Scholar
  12. D. Traganos, D. Poursanidis, B. Aggarwal, N. Chrysoulakis, and P. Reinartz, “Estimating Satellite-Derived Bathymetry (SDB) with the Google Earth Engine and Sentinel-2,” Remote Sensing, vol. 10, no. 6, p. 859, 2018. View at: Publisher Site | Google Scholar
  13. P. Jagalingam, B. Akshaya, and A. V. Hegde, “Bathymetry mapping using Landsat 8 satellite imagery,” Procedia Engineering, vol. 116, pp. 560–566, 2015. View at: Publisher Site | Google Scholar
  14. Z. Lee, K. L. Carder, R. F. Chen, and T. G. Peacock, “Properties of the water column and bottom derived from Airborne Visible Infrared Imaging Spectrometer (AVIRIS) data,” Journal of Geophysical Research: Oceans, vol. 106, no. C6, pp. 11639–11651, 2001. View at: Publisher Site | Google Scholar
  15. D. R. Lyzenga, “Remote sensing of bottom reflectance and water attenuation parameters in shallow water using aircraft and Landsat data,” International Journal of Remote Sensing, vol. 2, no. 1, pp. 71–82, 1981. View at: Publisher Site | Google Scholar
  16. F. Polcyn, W. L. Brown, and I. Sattinger, The Measurement of Water Depth by Remote Sensing Techniques, Michigan Univ Ann Arbor Inst of Science and Technology, 1970.
  17. A. G. Dekker, S. R. Phinn, J. Anstee et al., “Intercomparison of shallow water bathymetry, hydro-optics, and benthos mapping techniques in Australian and Caribbean coastal environments,” Limnology and Oceanography: Methods, vol. 9, no. 9, pp. 396–425, 2011. View at: Publisher Site | Google Scholar
  18. Z. Lee, K. L. Carder, C. D. Mobley, R. G. Steward, and J. S. Patch, “Hyperspectral remote sensing for shallow waters: 2. Deriving bottom depths and water properties by optimization,” Applied Optics, vol. 38, no. 18, pp. 3831–3843, 1999. View at: Google Scholar
  19. V. E. Brando, J. M. Anstee, M. Wettle, A. G. Dekker, S. R. Phinn, and C. Roelfsema, “A physics based retrieval and quality assessment of bathymetry from suboptimal hyperspectral data,” Remote Sensing of Environment, vol. 113, no. 4, pp. 755–770, 2009. View at: Publisher Site | Google Scholar
  20. Z. Lee and K. L. Carder, “Effect of spectral band numbers on the retrieval of water column and bottom properties from ocean color data,” Applied Optics, vol. 41, no. 12, pp. 2191–2201, 2002. View at: Google Scholar
  21. Y. Ma, N. Xu, Z. Liu et al., “Satellite-derived bathymetry using the ICESat-2 Lidar and Sentinel-2 imagery datasets,” Remote Sensing of Environment, vol. 250, article 112047, 2020. View at: Publisher Site | Google Scholar
  22. N. Thomas, A. P. Pertiwi, D. Traganos et al., “Space-borne cloud-native satellite-derived bathymetry (SDB) models using ICESat-2 and Sentinel-2,” Geophysical Research Letters, vol. 48, no. 6, 2021. View at: Publisher Site | Google Scholar
  23. I. Caballero and R. P. Stumpf, “Retrieval of nearshore bathymetry from Sentinel-2A and 2B satellites in South Florida coastal waters,” Estuarine, Coastal and Shelf Science, vol. 226, article 106277, 2019. View at: Publisher Site | Google Scholar
  24. Y. Liu, J. Zhao, R. Deng et al., “A downscaled bathymetric mapping approach combining multitemporal Landsat-8 and high spatial resolution imagery: demonstrations from clear to turbid waters,” ISPRS Journal of Photogrammetry and Remote Sensing, vol. 180, pp. 65–81, 2021. View at: Publisher Site | Google Scholar
  25. Z. Lee, M. Shangguan, R. A. Garcia et al., “Confidence measure of the shallow-water bathymetry map obtained through the fusion of Lidar and multiband image data,” Journal of Remote Sensing, vol. 2021, article 9841804, pp. 1–16, 2021. View at: Publisher Site | Google Scholar
  26. B. Ai, Z. Wen, Z. Wang et al., “Convolutional neural network to retrieve water depth in marine shallow water area from remote sensing images,” Ieee Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 13, pp. 2888–2898, 2020. View at: Publisher Site | Google Scholar
  27. S. Liu, L. Wang, H. Liu, H. Su, X. Li, and W. Zheng, “Deriving bathymetry from optical images with a localized neural network algorithm,” IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 9, pp. 5334–5342, 2018. View at: Publisher Site | Google Scholar
  28. Z. Lee, J. Sandidge, and M. Zhang, “Ocean-Color Inversion: A Combined Approach by Analytical Solution and Neural Networks,” in Ocean Remote Sensing and Imaging II, San Diego, California, United States, 2003. View at: Google Scholar
  29. J. C. Sandidge and R. J. Holyer, “Coastal bathymetry from hyperspectral observations of water radiance,” Remote Sensing of Environment, vol. 65, no. 3, pp. 341–352, 1998. View at: Publisher Site | Google Scholar
  30. I. Caballero and R. P. Stumpf, “Atmospheric correction for satellite-derived bathymetry in the Caribbean waters: from a single image to multi-temporal approaches using Sentinel-2A/B,” Optics Express, vol. 28, no. 8, pp. 11742–11766, 2020. View at: Google Scholar
  31. M. Wang, “A refinement for the Rayleigh radiance computation with variation of the atmospheric pressure,” International Journal of Remote Sensing, vol. 26, no. 24, pp. 5651–5663, 2005. View at: Publisher Site | Google Scholar
  32. M. Wang, “The Rayleigh lookup tables for the SeaWiFS data processing: accounting for the effects of ocean surface roughness,” International Journal of Remote Sensing, vol. 23, no. 13, pp. 2693–2702, 2002. View at: Google Scholar
  33. E. Micijevic, M. O. Haque, and N. Mishra, “Radiometric Calibration Updates to the Landsat Collection,” in Earth Observing Systems XXI, San Diego, California, United States, 2016. View at: Google Scholar
  34. J. Wei, Z. Lee, R. Garcia et al., “An assessment of Landsat-8 atmospheric correction schemes and remote sensing reflectance products in coral reefs and coastal turbid waters,” Remote Sensing of Environment, vol. 215, pp. 18–32, 2018. View at: Publisher Site | Google Scholar
  35. B. A. Franz, S. W. Bailey, N. Kuring, and P. J. Werdell, “Ocean color measurements with the Operational Land Imager on Landsat-8: implementation and evaluation in SeaDAS,” Journal of Applied Remote Sensing, vol. 9, no. 1, article 096070, 2015. View at: Publisher Site | Google Scholar
  36. J. Wang, Z. Lee, D. Wang, S. Shang, J. Wei, and A. Gilerson, “Atmospheric correction over coastal waters with aerosol properties constrained by multi-pixel observations,” Remote Sensing of Environment, vol. 265, article 112633, 2021. View at: Publisher Site | Google Scholar
  37. H. R. Gordon and M. Wang, “Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS: a preliminary algorithm,” Applied Optics, vol. 33, no. 3, pp. 443–452, 1994. View at: Google Scholar
  38. F. Steinmetz, P.-Y. Deschamps, and D. Ramon, “Atmospheric correction in presence of sun glint: application to MERIS,” Optics Express, vol. 19, no. 10, pp. 9783–9800, 2011. View at: Google Scholar
  39. G. Thuillier, M. Hersé, P. C. Simon, H. Mandel, and D. Gillotay, “Observation of the solar spectral irradiance from 200 nm to 870 nm during the ATLAS 1 and ATLAS 2 missions by the SOLSPEC spectrometer,” Metrologia, vol. 35, no. 4, p. 689, 1998. View at: Google Scholar
  40. T. A. Neumann, A. Brenner, D. Hancock et al., ATLAS/ICESat-2 L2A Global Geolocated Photon Data, Version 3, NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado USA, 2020. View at: Publisher Site
  41. M. Ester, H. P. Kriegel, J. Sander, and X. Xu, “A density-based algorithm for discovering clusters in large spatial databases with noise,” in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, Portland, Oregon, 1996. View at: Google Scholar
  42. C. E. Parrish, L. A. Magruder, A. L. Neuenschwander, N. Forfinski-Sarkozi, M. Alonzo, and M. Jasinski, “Validation of ICESat-2 ATLAS bathymetry and analysis of ATLAS's bathymetric mapping performance,” Remote Sensing, vol. 11, no. 14, p. 1634, 2019. View at: Publisher Site | Google Scholar
  43. Y. Chen, Y. Le, D. Zhang, Y. Wang, Z. Qiu, and L. Wang, “A photon-counting LiDAR bathymetric method based on adaptive variable ellipse filtering,” Remote Sensing of Environment, vol. 256, article 112326, 2021. View at: Publisher Site | Google Scholar
  44. R. Pawlowicz, B. Beardsley, and S. Lentz, “Classical tidal harmonic analysis including error estimates in MATLAB using T-TIDE,” Computers & Geosciences, vol. 28, no. 8, pp. 929–937, 2002. View at: Google Scholar
  45. Z. Lee, C. Hu, R. Arnone, and Z. Liu, “Impact of sub-pixel variations on ocean color remote sensing products,” Optics Express, vol. 20, no. 19, pp. 20844–20854, 2012. View at: Google Scholar
  46. R. A. Garcia, Z. Lee, B. B. Barnes, C. Hu, H. M. Dierssen, and E. J. Hochberg, “Benthic classification and IOP retrievals in shallow water environments using MERIS imagery,” Remote Sensing of Environment, vol. 249, article 112015, 2020. View at: Publisher Site | Google Scholar
  47. A. W. Bruckner, J. Kerr, G. Rowlands, A. Dempsey, S. J. Purkis, and P. Renaud, Khaled bin Sultan Living Oceans Foundation Atlas of Shallow Marine Habitats of Cay Sal Bank, Panoramic Press, Great Inagua, Little Inagua and Hogsty Reef, Bahamas, 2014.
  48. M. I. Jordan and T. M. Mitchell, “Machine learning: trends, perspectives, and prospects,” Science, vol. 349, no. 6245, pp. 255–260, 2015. View at: Google Scholar
  49. G. E. Dahl, T. N. Sainath, and G. E. Hinton, “Improving Deep Neural Networks for LVCSR Using Rectified Linear Units and Dropout,” in 2013 IEEE international conference on acoustics, speech and signal processing, Vancouver, BC, Canada, 2013. View at: Google Scholar
  50. C. Nwankpa, W. Ijomah, A. Gachagan, and S. Marshall, “Activation functions: comparison of trends in practice and research for deep learning,” 2018, View at: Google Scholar
  51. D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” Nature, vol. 323, no. 6088, pp. 533–536, 1986. View at: Google Scholar
  52. D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014, View at: Google Scholar
  53. M. Abadi, P. Barham, J. Chen et al., “Tensorflow: a system for large-scale machine learning,” in 12th {USENIX} symposium on operating systems design and implementation ({OSDI} 16), pp. 265–283, Savannah, GA, USA, 2016. View at: Google Scholar
  54. A. Géron, Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow: concepts, tools, and techniques to build intelligent systems, O'Reilly Media, 2019.
  55. G. D. Egbert, A. F. Bennett, and M. G. Foreman, “TOPEX/POSEIDON tides estimated using a global inverse model,” Journal of Geophysical Research: Oceans, vol. 99, no. C12, pp. 24821–24852, 1994. View at: Google Scholar
  56. E. C. Franklin, J. S. Ault, S. G. Smith et al., “Benthic habitat mapping in the Tortugas region, Florida,” Marine Geodesy, vol. 26, no. 1-2, pp. 19–34, 2003. View at: Google Scholar
  57. I. Caballero, R. P. Stumpf, and A. Meredith, “Preliminary assessment of turbidity and chlorophyll impact on bathymetry derived from Sentinel-2A and Sentinel-3A satellites in South Florida,” Remote Sensing, vol. 11, no. 6, p. 645, 2019. View at: Google Scholar
  58. H. R. Gordon, J. W. Brown, and R. H. Evans, “Exact Rayleigh scattering calculations for use with the Nimbus-7 coastal zone color scanner,” Applied Optics, vol. 27, no. 5, pp. 862–871, 1988. View at: Google Scholar

Copyright © 2022 Wendian Lai et al. Exclusive Licensee Aerospace Information Research Institute, Chinese Academy of Sciences. Distributed under a Creative Commons Attribution License (CC BY 4.0).

 PDF Download Citation Citation
Altmetric Score