Research Article  Open Access
Wendian Lai, Zhongping Lee, Junwei Wang, Yongchao Wang, Rodrigo Garcia, Huaguo Zhang, "A Portable Algorithm to Retrieve Bottom Depth of Optically Shallow Waters from TopOfAtmosphere Measurements", Journal of Remote Sensing, vol. 2022, Article ID 9831947, 16 pages, 2022. https://doi.org/10.34133/2022/9831947
A Portable Algorithm to Retrieve Bottom Depth of Optically Shallow Waters from TopOfAtmosphere Measurements
Abstract
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 Landsat8 and ICESat2 as examples, we present an approach to retrieve directly from the topofatmosphere (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 Rayleighcorrected 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 Landsat8 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 twoband 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 [3–6].
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 highcost, timeconsuming, and limited to the areas of vessels or airplanes can reach, with the resulted data from space hardly to form a highresolution 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, 10–12]. 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 [12–14]. Although its accuracy and precision cannot match that from Lidar, the imaging capability offers an efficient and costeffective 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, airsea 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 abovewater radiometry [4, 12, 15–17]. 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 [17–19] do not require known depths for algorithm development, they require adequate spectral bands in order to separate watercolumn 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 Satellite2 (ICESat2), highprecision 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 spatialresolution imagers, such as Landsat8 (30 m resolution), Sentinel2 (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 [21–23]. One typical example is the TwoBand 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 [26–29]. 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 ICESat2 and imager data (using OLI of Landsat8 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 topofatmosphere 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 topofatmosphere reflectance () as the input, along with a MLP for the estimation of . More importantly, we further applied the trained NN algorithm to Landsat8 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, highresolution, bathymetry map of global OSW waters.
2. Data
2.1. Landsat8 Satellite Image
While many multiband imagers can be used for the empirical retrieval of , we used Landsat8 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 (NIRSWIR)(3)Topofatmosphere radiance from Landsat8 OLI is well calibrated [33, 34], thus adequate for the calculation of Rayleighcorrected reflectance ()
In addition, Landsat8 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 [36–38] where is the total reflectance at topofatmosphere, 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 Landsat8 OLI measurements were also generated using SeaDAS. For these products, pixels designated as lowquality were removed based on the standard Level2 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. ICESat2 Data
ICESat2 was launched on September 15, 2018, and carries ATLAS, a photoncounting 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 acrosstrack 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 ICESat2 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 (https://nsidc.org/data/atl03#) [40].
The procedure to obtain the bottom depth from ICESat2 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 WGS84 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 WGS84 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 densitybased 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 lowdensity 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 alongtrack 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 alongtrack direction due to tide or other hydrodynamic processes (Figure 1(b)). According to the photon geolocation algorithm for ATL03 of ICESat2, photons that are characterized as likely signals are binned in approximately 20 m alongtrack segments [40]. Within each segment, the water depth is calculated as the difference between the height of bottomreflected photons and the mean height of surfacereflected photons.(5)Correcting the impact of light speed: photons from ICESat2 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 lightcorrection 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 Landsat8 and ICESat2
Assuming the variation of bottom depth (after tidal correction) is negligible within a short period, a time constraint of ±2 weeks between Landsat8 OLI and ICESat2 measurements was used to matchup radiometric and Lidar data [25]. Meanwhile, the corrected bottom depth of ICESat2 was adjusted to match the tidal cycle of the Landsat8 overpass, where the classical tidal harmonic analysis model [44] was applied to calculate tide information of targeted locations.
To match the measurements between ICESat2 and OLI, we first located the ICESat2 track within the OLI image. For a point of ICESat2 corrected bottom depth, an OLI pixel was first selected based on the closest Euclidean distance. Since the footprint of ICESat2 along the ground track is 0.7 m, while the spatial resolution of OLI is 30 m, there are many ICESat2 points within a Landsat8 pixel. The weighted average of ICESat2 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 ICESat2 points within the OLI pixel, is the ICESat2 bottom depth from Section 2.2, and and are the tide heights at the time of the ICESat2 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 ICESat2 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 Landsat8 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.0010.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 (Landsat8 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 ICESat2, 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 ICESat2 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 Landsat8 and ICESat2 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.

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 neuralnetworkbasedTOAalgorithmforbathymetry, NNTOAB. Figures 6(a) and 6(b) present schematic charts of the two MLP systems. Fundamentally, MLP is a forwardstructured 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 easytouse interface of Tensorflow [53], was used to train the two MLPs. Tensorflow is an opensource 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 userfriendliness [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 Landsat8 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.

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, NNTOAB achieved an and MARD of 0.96 and 8.8%, respectively, for ranging between ~0 and 25 m, suggesting high predictability of NNTOAB. 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 dropoffs 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, NNTOAB 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 Landsat8 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 highresolution (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 ICESat2 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 NNTOAB 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 NNTOAB 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 Landsat8 and ICESat2 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 BandRatio Algorithm
The empirical twoband 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 Landsat8 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).

As the evaluation of NNTOAB 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 NNTOAB can obtain highly accurate estimation of from OLIderived (430–2300 nm) and that NNTOAB has high portability, implying that reliable could be obtained from new OLI acquisitions without modifying NNTOAB. One of the likely reasons for the promising results lies in the simultaneous use of the visibleSWIR domain where traditionally twoseparate 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 NIRSWIR 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. NNTOAB does not explicitly separate the two processes for the retrieval of from , but since all spectral information in the visibleSWIR 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 topofatmosphere radiance in the bluegreen 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 NNTOAB from OLI measurements are somewhat surprising, as Landsat8 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 HOPEbased 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 NNTOAB 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 Landsat8 OLI images of various latitudes, the evaluation of NNTOAB will be expanded; and, as all algorithm development, the present NNTOAB will likely be revised and updated to cover midtohigh latitude regions. In addition, as the focus of this effort is on the development and evaluation of , the confidence at pixel level [25] of NNTOAB 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 Landsat8 OLI and ICESat2, we show that bottom depth of OSW can be retrieved with high accuracy from the Rayleighcorrected topofatmosphere reflectance spectrum after a NNbased workflow (NNTOAB) was successfully trained with a large number and a wide range of data. This NNTOAB shows promising portability to other shallow regions not included in the training, thus overcoming a longstanding limitation of traditional empirical algorithms for that works best with data/image used in algorithm training. With data from ICESat2 and highresolution imagery growing by day, we see the generation of accurate, highresolution, 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 Landsat8 OLI data are publicly available and can be downloaded from https://glovis.usgs.gov/app. The ICESat2 data are also publicly available through the National Snow and Ice Data Center (NSIDC). The geolocated photon data (ATL03) are found online (https://nsidc.org/data/atl03), 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 Landsat8 and ICESat2 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.
Acknowledgments
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 Landsat8 OLI L1 imagery, NASA’s Ocean Biology Processing Group (OBPG) for the SeaDAS processing software, and the NASA ICESat2 team for providing the data used in this study.
References
 T. Kutser, I. Miller, and D. L. Jupp, “Mapping coral reef benthic substrates using hyperspectral spaceborne images and spectral libraries,” Estuarine, Coastal and Shelf Science, vol. 70, no. 3, pp. 449–460, 2006. View at: Google Scholar
 R. P. Stumpf, K. Holderied, and M. Sinclair, “Determination of water depth with highresolution satellite imagery over variable bottom types,” Limnology and Oceanography, vol. 48, no. 1, pp. 547–556, 2003. View at: Google Scholar
 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
 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
 M. F. Karim and N. Mimura, “Impacts of climate change and sealevel 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
 B. B. Barnes, R. Garcia, C. Hu, and Z. Lee, “Multiband 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
 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
 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
 J. L. Irish and T. E. White, “Coastal engineering applications of highresolution Lidar bathymetry,” Coastal Engineering, vol. 35, no. 12, pp. 47–71, 1998. View at: Publisher Site  Google Scholar
 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
 B. J. Babbel, C. E. Parrish, and L. A. Magruder, “ICESat2 elevation retrievals in support of satellitederived bathymetry for global science applications,” Geophysical Research Letters, vol. 48, no. 5, 2021. View at: Publisher Site  Google Scholar
 D. Traganos, D. Poursanidis, B. Aggarwal, N. Chrysoulakis, and P. Reinartz, “Estimating SatelliteDerived Bathymetry (SDB) with the Google Earth Engine and Sentinel2,” Remote Sensing, vol. 10, no. 6, p. 859, 2018. View at: Publisher Site  Google Scholar
 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
 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
 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
 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.
 A. G. Dekker, S. R. Phinn, J. Anstee et al., “Intercomparison of shallow water bathymetry, hydrooptics, 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
 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
 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
 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
 Y. Ma, N. Xu, Z. Liu et al., “Satellitederived bathymetry using the ICESat2 Lidar and Sentinel2 imagery datasets,” Remote Sensing of Environment, vol. 250, article 112047, 2020. View at: Publisher Site  Google Scholar
 N. Thomas, A. P. Pertiwi, D. Traganos et al., “Spaceborne cloudnative satellitederived bathymetry (SDB) models using ICESat2 and Sentinel2,” Geophysical Research Letters, vol. 48, no. 6, 2021. View at: Publisher Site  Google Scholar
 I. Caballero and R. P. Stumpf, “Retrieval of nearshore bathymetry from Sentinel2A and 2B satellites in South Florida coastal waters,” Estuarine, Coastal and Shelf Science, vol. 226, article 106277, 2019. View at: Publisher Site  Google Scholar
 Y. Liu, J. Zhao, R. Deng et al., “A downscaled bathymetric mapping approach combining multitemporal Landsat8 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
 Z. Lee, M. Shangguan, R. A. Garcia et al., “Confidence measure of the shallowwater 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
 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
 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
 Z. Lee, J. Sandidge, and M. Zhang, “OceanColor 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
 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
 I. Caballero and R. P. Stumpf, “Atmospheric correction for satellitederived bathymetry in the Caribbean waters: from a single image to multitemporal approaches using Sentinel2A/B,” Optics Express, vol. 28, no. 8, pp. 11742–11766, 2020. View at: Google Scholar
 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
 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
 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
 J. Wei, Z. Lee, R. Garcia et al., “An assessment of Landsat8 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
 B. A. Franz, S. W. Bailey, N. Kuring, and P. J. Werdell, “Ocean color measurements with the Operational Land Imager on Landsat8: implementation and evaluation in SeaDAS,” Journal of Applied Remote Sensing, vol. 9, no. 1, article 096070, 2015. View at: Publisher Site  Google Scholar
 J. Wang, Z. Lee, D. Wang, S. Shang, J. Wei, and A. Gilerson, “Atmospheric correction over coastal waters with aerosol properties constrained by multipixel observations,” Remote Sensing of Environment, vol. 265, article 112633, 2021. View at: Publisher Site  Google Scholar
 H. R. Gordon and M. Wang, “Retrieval of waterleaving 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
 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
 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
 T. A. Neumann, A. Brenner, D. Hancock et al., ATLAS/ICESat2 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
 M. Ester, H. P. Kriegel, J. Sander, and X. Xu, “A densitybased 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
 C. E. Parrish, L. A. Magruder, A. L. Neuenschwander, N. ForfinskiSarkozi, M. Alonzo, and M. Jasinski, “Validation of ICESat2 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
 Y. Chen, Y. Le, D. Zhang, Y. Wang, Z. Qiu, and L. Wang, “A photoncounting LiDAR bathymetric method based on adaptive variable ellipse filtering,” Remote Sensing of Environment, vol. 256, article 112326, 2021. View at: Publisher Site  Google Scholar
 R. Pawlowicz, B. Beardsley, and S. Lentz, “Classical tidal harmonic analysis including error estimates in MATLAB using TTIDE,” Computers & Geosciences, vol. 28, no. 8, pp. 929–937, 2002. View at: Google Scholar
 Z. Lee, C. Hu, R. Arnone, and Z. Liu, “Impact of subpixel variations on ocean color remote sensing products,” Optics Express, vol. 20, no. 19, pp. 20844–20854, 2012. View at: Google Scholar
 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
 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.
 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
 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
 C. Nwankpa, W. Ijomah, A. Gachagan, and S. Marshall, “Activation functions: comparison of trends in practice and research for deep learning,” 2018, https://arxiv.org/abs/1811.03378. View at: Google Scholar
 D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by backpropagating errors,” Nature, vol. 323, no. 6088, pp. 533–536, 1986. View at: Google Scholar
 D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014, https://arxiv.org/abs/1412.6980. View at: Google Scholar
 M. Abadi, P. Barham, J. Chen et al., “Tensorflow: a system for largescale 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
 A. Géron, Handson machine learning with ScikitLearn, Keras, and TensorFlow: concepts, tools, and techniques to build intelligent systems, O'Reilly Media, 2019.
 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
 E. C. Franklin, J. S. Ault, S. G. Smith et al., “Benthic habitat mapping in the Tortugas region, Florida,” Marine Geodesy, vol. 26, no. 12, pp. 19–34, 2003. View at: Google Scholar
 I. Caballero, R. P. Stumpf, and A. Meredith, “Preliminary assessment of turbidity and chlorophyll impact on bathymetry derived from Sentinel2A and Sentinel3A satellites in South Florida,” Remote Sensing, vol. 11, no. 6, p. 645, 2019. View at: Google Scholar
 H. R. Gordon, J. W. Brown, and R. H. Evans, “Exact Rayleigh scattering calculations for use with the Nimbus7 coastal zone color scanner,” Applied Optics, vol. 27, no. 5, pp. 862–871, 1988. View at: Google Scholar
Copyright
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).