Exploring the Most Effective Information for Satellite-Derived Bathymetry Models in Different Water Qualities

Abstract

Satellite-derived bathymetry (SDB) is an effective means of obtaining global shallow water depths. However, the effect of inherent optical properties (IOPs) on the accuracy of SDB under different water quality conditions has not been clearly clarified. To enhance the accuracy of machine learning SDB models, this study aims to assess the performance improvement of integrating the quasi-analytical algorithm (QAA)-derived IOPs using the Sentinel-2 and ICESat-2 datasets. In different water quality experiments, the results indicate that four SDB models (the Gaussian process regression, neural networks, random forests, and support vector regression) incorporating QAA-IOP parameters equal to or outperform those solely based on the remote sensing reflectance (R_rs) datasets, especially in turbid waters. By analyzing information gains in SDB, the most effective inputs are identified and prioritized under different water qualities. The SDB method incorporating QAA-IOP can achieve an accuracy of 0.85 m, 0.48 m, and 0.74 m in three areas (Wenchang, Laizhou Bay, and the Qilian Islands) with different water quality. Also, we find that incorporating an excessive number of redundant bands into machine learning models not only increases the demand of computing resources but also leads to worse accuracy in SDB. In conclusion, the integration of QAA-IOPs offers promising improvements in obtaining bathymetry and the optimal feature selection should be carefully considered in diverse aquatic environments.

1. Introduction

Coastal regions are among the most economically significant ecosystems [1]. Accurate shallow water depth data are essential to marine engineering, maritime economics, and integrated coastal zone management [2] and are vital for research on nearshore marine dynamics [3]. Basically, all marine activities and maritime spatial planning require depth data [4]. Since survey vessels may not operate in very shallow waters, unmanned aerial vehicles (UAVs) equipped with acoustic bathymetric instruments [5] and airborne bathymetric lidars [6,7] offer a precise and effective method for depth measurement. However, while shipborne and airborne bathymetric instruments can provide very precise depth estimations, their spatial coverages are limited to regional scales rather than a global scale [8,9].

Satellite-derived bathymetry (SDB) offers a rapid and efficient provision of high-resolution underwater topographic maps and acts as a valuable complement to traditional bathymetric techniques [10]. The principle of estimating water depths from multispectral satellite images is based on the premise that the light penetration through the water column at different wavelengths is a function of seawater characteristics [11,12], and models are constructed under various observational conditions to interpret the relationship between multi-dimensional feature values and water depth. Many reliable methods for deriving bathymetry have been developed over the past fifty years [13]. For empirical algorithms, earlier approaches to SDB relied on simple but effective mathematical models constructed from differential reflectance across various bands, such as Lyzenga’s linear band combination model [14] and Stumpf’s logarithmic band ratio model [15]. In develo** versatile SDB models, machine learning is perceived as an efficacious methodology, leveraging its capacity for fitting to autonomously explore numerical models and proffer optimal solutions [16]. In recent years, various machine learning algorithms, including neural networks [17,18], random forests [19,20,21], and support vector machines [21,22], have been employed in building SDB models.

Conventionally, most of the SDB studies used in situ measured depth data as control (or training) points to construct empirical SDB models, which limited their applications in some remote areas without available in situ data. Many physical-based SDB models that do not rely on in situ data were also proposed [23,24,25,26], but these have the cost that more spectral information is needed. ICESat-2, as the first photon-counting lidar satellite, provides a new means of directly acquiring bathymetric points worldwide in place of in situ bathymetric data. Within the water depths of 30 m, the accuracy of ICESat-2-derived bathymetric points is generally better than 0.6 m [27]. A number of studies have confirmed that ICESat-2 data are compatible with a wide range of hyperspectral and multispectral data sources including Landsat-8, Sentinel-2, Worldview-2, Pleiades-1, and Zhuhai-1 to achieve SDB using only satellite-based datasets [9,28,29,30,31,32,33,34,35,36,37].

Current studies based on empirical algorithms for SDB typically input all available features into machine learning models. In addition to the remote sensing reflectance R_rs that are directly observed by spaceborne optical sensors, some studies also considered the absorption coefficients at specific wavelengths or other inherent optical parameters (IOPs) when training the SDB models [38,39,40]. The absorption and scattering effect of the water column precludes deep penetration of the light and greatly influence the remote sensing reflectance R_rs [41]. In other words, differences in IOPs can lead to significant differences in R_rs even at the same water depth, and thus inputting IOP information can improve the bathymetric accuracy in complex coastal scenarios. However, on one hand, present machine learning studies on SDB are mostly concentrated along coasts or on islands far from urban areas, where the water qualities are very good. The empirical SDB band selection methods established for clear water bodies probably perform poorly in relatively turbid waters [42]. On the other hand, incorporating an excessive number of redundant bands into machine learning models can lead to worse accuracy in SDB and increased computational demand [43]. As a result, identifying the most effective bands under different water quality conditions is well worth investigating. The quasi-analytical algorithm (QAA) is a classic model for IOP retrievals [11], whose accuracy in coastal areas has been verified as reliable in both clear and turbid waters [44,45,46]. Utilizing machine learning models and the QAA algorithm, the effectiveness of input information for SDB models can be evaluated under varying water quality conditions.

This study aims to derive a machine learning SDB model and increase its accuracy by integrating the QAA algorithm, and then quantitatively analyze the most effective input features under different water quality conditions. Specifically, we extract ICESat-2 bathymetric points as ground truth data from three study areas with different water qualities. After preprocessing Sentinel-2 data, we expand the dimensionality of the R_rs dataset from Sentiel-2 to QAA-derived IOPs in different bands. This expansion is then processed by a dimensionality reduction through the wrapper method. Reducing the dimensions of the expanded dataset can be viewed as an evaluation of the most effective input information. We use four machine learning algorithms, i.e., the Gaussian process regression, neural networks, random forests, and support vector regression, to develop their corresponding SDB models. The accuracy and robustness of this new method are subsequently evaluated in three study areas. By assessing information gains, this study presents the information quantities of different features in SDB and gives priority rankings under different water quality conditions.

2. Study Areas and Dataset

Three areas in China are selected with different water qualities (shown in Figure 1). The first area is located in the South China Sea and near the coast of Wenchang [19.7810°N, 111.0021°E]. This coastal area is generally rural with relatively good water quality, which belongs to water type III based on the Jerlov classification [47].

The second area is located in Laizhou Bay [37.2724°N, 119.8064°E] and near the estuary of the Yellow River. This coastal area contains a high sediment load, numerous aquaculture zones and docks, and belongs to water type 7C based on the Jerlov classification. The third area is located in the Qilian Islands [16.9803°N, 112.2386°E] of the South China Sea. This coral reef has very good water quality, which belongs to water type IB based on the Jerlov classification. Figure 1 displays the study area using the Sentinel-2 satellite imagery as the basemap, where the yellow tracks represent the ICESat-2 data used in the experiment.

Table 1. Study area information and Sentinel-2 and ICESat-2 data acquisition time.

Location Center Coordinates	(A) Wenchang				(B) Laizhou Bay			(C) Qilian Islands
	19.7810°N, 111.0021°E				37.2724°N, 119.8064°E			16.9803°N, 112.2386°E
Inherent optical parameters, IOPs (m−1)	a488	bb488		Kd490	a488	bb488	Kd490	a488	bb488	Kd490
	0.09	0.011		0.21	0.23	0.032	0.58	0.03	0.006	0.04
Sentinel-2 image	21 June 2020				9 September 2019			4 April 2023
ICESat-2 points (training and validation)	28 January 2019		20 June 2020		26 October 2019			13 April 2021
	29 April 2019		26 July 2020		21 March 2020			16 July 2021
	22 June 2019		18 June 2021		19 September 2020			14 January 2022
	29 July 2019		24 July 2021		24 December 2020			15 July 2022
	21 September 2019		17 September 2021		17 September 2021			9 January 2023

lists the details of all the data utilized in this study.

2.1. Sentinel-2 Data

The European Space Agency’s (ESA) Copernicus Sentinel-2 mission includes two polar-orbiting satellites carrying the MultiSpectral Instrument (MSI), which provides thirteen available spectral bands. The Sentinel-2 data utilized in this study can be obtained from the European Space Agency (ESA) Sentinel Scientific Data Hub. In this study, we utilize the 10 m resolution bands at 490 nm, 560 nm, and 665 nm, 20 m resolution bands at 865 nm, and 60 m resolution bands at 443 nm. The Sentinel-2 Level-1C top-of-atmosphere reflectance data released by the ESA have undergone geometric correction and radiometric calibration and are geometrically projected in the UTM/WGS84 (Universal Transverse Mercator/World Geodetic System 1984) framework. During the processing of Level-1C datasets to Level-2A, the atmospheric correction is performed using the C2RCC plugin in the SNAP software v9.0.0. Compared to the atmospheric correction by Sen2Cor, C2RCC exhibits better performance in nearshore seawater [35]. Subsequently, the Deglint Operator plugin from the Sen2Cor toolbox is employed for deglint, using NIR (near-infrared) or SWIR (short-wave infrared) to estimate and correct the sun glint in the visible bands. After the above processing, the remote sensing reflectance (R_rs) at 443 nm, 490 nm, 560 nm, 665 nm, and 865 nm are used to construct the SDB dataset. Specifically, the 10 m resolution R_rs at 443 nm and 865 nm are obtained from their 60 m and 20 m resolution R_rs by using the resampling tool in the SNAP software.

2.2. ICESat-2 Data

The payload on ICESat-2 has six laser beams aligned along its trajectory and the cross-track distance is 3.3 km between adjacent beam pairs. A beam pair consists of two laser beams, i.e., the strong beam has a pulse energy that is approximately four times stronger than that of the weak beam. The ICESat-2 Level-2 ATL03 Geolocated Photons product, sourced from https://search.earthdata.nasa.gov/search, accessed on 1 July 2023, is used in this study. The Level-2 ATL03 dataset includes all geolocated photons (including both signal and noise photons) along six tracks (three strong beams and three weak beams). Each photon possesses a unique latitude, longitude, and elevation on the WGS84 ellipsoid reference. Corrections have been conducted to the dataset for the atmospheric delay, solid earth tides, and system pointing biases. However, bathymetry errors caused by surface water fluctuations, instantaneous sea surface heave, etc., are not removed, and this part of the error correction will be introduced in Section 3.

3. Methods

In this study, we employ the density-based spatial clustering of applications with noise (DBSCAN) method to extract surface and bottom photons from ATL03 data and correct the bathymetric errors for the bottom photons. For Sentinel-2 data, by applying the quasi-analytical algorithm (QAA) for four visible light bands of Sentinel-2, the dimensionality of the dataset is expanded to twelve channels. To minimize information redundancy and enhance the efficiency of machine learning models, the wrapper method is then applied to reduce the data dimensionality. Combining the bathymetric points from ICESat-2 and multiple spectral information from Sentinel-2, four machine learning algorithms are utilized to obtain the SDB. To compare and assess the improvement with and without incorporating IOPs, four machine learning models are separately trained using the IOPs+ R_rs dataset or only using R_rs. The effectiveness of these models was evaluated by comparing them with the bathymetric results from ICESat-2. The machine learning techniques used included Gaussian process regression (GPR), neural network (NN), random forest (RF), and support vector regression (SVR). The overall workflow is shown in Figure 2.

3.1. ICESat-2 Bottom Photon Detection

The ATL03 product includes signal and noise photons, and this study employs the DBSCAN method to detect the bottom signal photons. The fundamental assumption of DBSCAN denoising is that the density of signal photons is spatially higher than that of noise photons. Specifically, every 10,000 continuous geolocated photons are segmented in the along-track distance. Within each segment, the cluster density of each photon is evaluated as signal if the photon number exceeds the minimum threshold MinPts with a given radius R_a. To prevent the misidentification of solar background noise as signal in the daytime, R_a for daytime data is set as 1.5 m, whereas for nighttime data, it is set as 2.5 m. The minimum threshold MinPts is adaptively calculated as

$$\begin{array}{c}MinPts=\frac{2S{N}_1-S{N}_2}{\ln \left(\frac{2S{N}_1}{S{N}_2}\right)} \end{array}$$

(1)

where SN₁ represents the expected number of signal and noise photons within a cluster, while SN₂ indicates the expected number of noise photons in the water column [48]. Normally, the vertical layer with the lowest elevation does not include seabed reflection signals, but only contains water backscatter and background noise.

Using the DBSCAN method, sea surface and bottom signal photons are detected from the ATL03 geolocated photons and manually processed to further discard possible noise photons. Then, the bathymetric errors of bottom photons need further correction. In this study, we consider and correct the effects of the water column refraction and forward scattering, sea surface refraction, and sea surface undulations [27,30,35,49]. For each bottom photon, the average elevation of nearby sea surface photons is taken as the instantaneous mean sea surface level, which is then subtracted from the elevation of current bottom photon to determine the instantaneous water depth at that moment. Next, the FES2014 model is used to compensate for water depth variations due to the tidal effect [50,51]. Following the aforementioned procedures, the water depths from the ICESat-2 bathymetric points are obtained.

3.2. QAA for Sentinel-2

In coastal areas, especially near urban areas, water color varies widely by region. Using only remote sensing reflectance for bathymetry may be unreliable due to the significant differences in IOPs. The QAA is a semi-analytical model that adopts the bio-optical model proposed by Lee et al. [11], which calculates the IOPs such as the total absorption (a) and particle backscattering coefficient (b_bp) for various wavelengths using multiple spectral reflectance. Utilizing the QAA provides extra IOP parameters for bathymetry. In the latest version of QAA_v6, the central wavelengths of the two reference bands are 555 nm and 665 nm. In this study, we approximate the bands of 555 nm and 670 nm to that of 560 nm and 665 nm, whose R_rs are provided by the MSI sensor. The reliability of the above substitution has been extensively validated [28,29,30]. The wavelengths used in the QAA are 443 nm, 490 nm, 560 nm, and 665 nm. The 10 m resolution images at 443 nm and 865 nm are resampled by the SNAP software.

QAA_v6 estimates a and b_bp in all other bands by selecting a reference band λ₀, which can be expressed as

$$\begin{array}{c}{b}_{bp}\left(\lambda \right)=b_{bp}\left({\lambda }_0\right){\left(\frac{{\lambda }_0}{\lambda }\right)}^{\eta }\end{array}$$

(2)

$$\begin{array}{c}a\left(\lambda \right)=\left[1-u\left(\lambda \right)\right]\left[b_w\left(\lambda \right)+b_{bp}\left(\lambda \right)\right]/u\left(\lambda \right)\end{array}$$

(3)

a_w and b_bw are the IOPs of pure seawater and are fixed constants. The coefficients of u and η can be calculated in the QAA_v6 equations, which provides the detailed and complete version of the formula and empirical coefficients. For the reference band λ₀, in relatively turbid and clear waters, two different empirical formulas are selected by the magnitude of R_rs(665) in this study in place of R_rs(670). When R_rs(665) ≥ 0.0015 sr⁻¹, the reference band λ₀ is at 665 nm, and a and b_bp at λ₀ at are defined as

$$\begin{array}{c}a\left({\lambda }_0\right)=a\left(665\right)=a_w\left(665\right)+0.39{\left(\frac{R_{rs}\left(665\right)}{R_{rs}\left(443\right)+R_{rs}\left(490\right)}\right)}^{1.14}\end{array}$$

(4)

$$\begin{array}{c}{b}_{bp}\left({\lambda }_0\right)=b_{bp}\left(665\right)=\frac{u\left({\lambda }_0\right)\times a\left({\lambda }_0\right)}{1-u\left({\lambda }_0\right)}-b_{bw}\left(665\right)\end{array}$$

(5)

When R_rs(665) < 0.0015 sr⁻¹, the reference band λ₀ is at 560 nm, and a, b_bp at λ₀ are calculated using the R_rs(665), R_rs(560), R_rs(490), and R_rs(443). All computational steps are detailed in the QAA_v6 Equation Table [52]. Following the application of the QAA algorithm, in addition to the original four band R_rs with 10 m resolution, eight IOP features are obtained, which are presented in

Table 2. Features after dataset expansion by QAA.

Rrs	490 nm	560 nm	665 nm	865 nm
a	443 nm	490 nm	560 nm	665 nm
bbp	443 nm	490 nm	560 nm	665 nm

.

3.3. Dimensionality Reduction

After processing by the QAA algorithm, the Sentinel-2 dataset encompasses 12 features. Some significant autocorrelations may exist among the QAA-derived IOP parameters. Moreover, an excessive number of feature dimensions can lead to an overcomplicated model and reduction in the accuracy of predicted results [43]. As a result, this study utilizes the wrapper method [53] for feature selection to preserve as much of the nonlinear enhancement relationship among IOPs and R_rs at different wavelengths, instead of directly reducing available dimensions using methods like the principal component analysis (PCA). The wrapper method is an iterative feature selection process, which uses the machine learning algorithm to select the optimal subset of features. By choosing the best-performing subset as the starting point for the next iteration, the iterative process continues until a predefined stop** criterion, such as a maximum number of iterations or a certain threshold, is satisfied. The computational process of the wrapper method is depicted in Figure 3.

In this study, a greedy algorithm is adopted as the search strategy, and the iteration process will terminate when no further decrease is obtained in the residual [54,55]. Clearly, the computational steps involved in the wrapper method can be time-consuming. To reduce the computation time, we choose the random forest algorithm as the iterative model among the four machine learning algorithms introduced in Section 3.4.

3.4. SDB Machine Learning Models

Compared with the classic linear band model and band ratio model [14,15], machine learning algorithms, which are capable of handling nonlinearities, high dimensionality, and potential complex interactive effects, can better model the complex relationships between the water depth and IOP/R_rs features specifically in relatively turbid waters. For bathymetric inversion in this study, we plan to use four types, i.e., the GPR, NN, RF, and SVR.

GPR is a non-parametric model that infers continuous distributions over functions for predictions, utilizing a kernel to define covariances among input data [56]. NN is a computational model, which simulates biological neural networks and learns complex data patterns through interconnected nodes with adjustable weights [57]. RF is an ensemble method employing multiple decision trees to enhance prediction accuracy and control over-fitting by aggregating their outcomes [58]. SVR extends SVM for regression, aiming to minimize deviations from actual targets within a specified margin and ensuring robust predictive performance.

Some hyperparameters need to be set when using these machine learning algorithms. Specifically, for GPR, the shape parameter α is set to 1, the length scale of the kernel l is set to 0.075, and the signal standard deviation is set to 2. For NN, we set the learning rate as 0.01, the momentum parameter as 0.9, the maximum number of iterations as 1000, and the iteration error limit as 1 × 10⁻⁶. For RF, we set the number of learners as 300 and defined the minimum leaf size as 8. For SVR, we utilize the radial basis function (RBF) as the kernel function, setting the kernel scale to 1. The box constraint is defined at 2.324; the epsilon is set at 0.2324; and the data undergoes a standardization.

3.5. Validation Metrics

In this study, we use three metrics, i.e., the root-mean-square error (RMSE), the coefficient of determination (R²), and the mean absolute percentage error (MAPE) to assess the prediction accuracy of the developed model, which are shown in Equations (6)–(8).

$$\begin{array}{c}RMSE=\sqrt{\frac{1}{n}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\left({\widehat{h}}_i-h_i\right)}^2}\end{array}$$

(6)

$$\begin{array}{c}{R}^2=\frac{{\displaystyle \sum _{i=1}^n{\left({\widehat{h}}_i-\overline{h}\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(h_i-\overline{h}\right)}^2}}\end{array}$$

(7)

$$MAPE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|\frac{\widehat{h_i}-h_i}{\widehat{h_i}}\right|}$$

(8)

where n represents the number of validation points, $h_i$ is the i-th truth depth, ${\widehat{h}}_i$ is the i-th estimated depth, and $\overline{h}$ is the average truth depth.

4. Results

4.1. Dataset Dimension Reduction and QAA Algorithm Results

As shown in Figure 4, the wrapper method identifies the optimal feature combination for Wenchang as R_rs(560), R_rs(865), a(443), and b_bp(490), and these four features compose the QAA-IOP dataset in the study area of Wenchang. It should be noted that multiple combinations of optimal features may exist, which offer an approximately equivalent accuracy, e.g., the combination of R_rs(490), R_rs (865), a(665), and b_bp(555) for Wenchang. Near Wenchang, the values of a(443) obtained through the QAA algorithm range from 0.04 to 0.18 m⁻¹, with an average of 0.12 m⁻¹. The values of b_bp(490) derived from the QAA algorithm range from 0.005 to 0.025 m⁻¹, with an average value of 0.008 m⁻¹.

A similar situation exists near Laizhou Bay, where multiple feature combinations could result in comparably precise results. As illustrated in Figure 5, near Laizhou Bay, the optimal feature combination obtained from the wrapper method includes R_rs(665), R_rs(865), a(560), and a(665). Near Laizhou Bay, a(560) ranges from 0.15 to 0.74 m⁻¹ with an average of 0.32 m⁻¹, and a(665) ranges from 0.42 to 0.80 m⁻¹ with an average of 0.65. As shown in Figure 6, the optimal feature combination for the Qilian Islands includes R_rs(490), R_rs(560), R_rs(665), and R_rs(865), which means under this very clear water quality, the performance of the QAA-IOP dataset is equivalent to original R_rs dataset.

Due to the differences of water types or qualities, the optimal bathymetric inversion combinations can be quite different. As a result, to construct more accurate SDB models, finding the optimal combination of features is recommended which considers the water optical properties in different study areas.

4.2. Results of Different SDB Models and Training Dataset

By matching the ICESat-2 bathymetry points with Sentinel-2 imagery, 5312 data pairs are obtained near Wenchang, 2451 pairs near Laizhou Bay, and 3612 pairs near the Qilian Islands. Figure 7 shows how the RMSE of each bathymetry model changes as the number of data training pairs increases. As indicated by the trends in Figure 7, with the exception of the Qilian Islands area, the dataset incorporating QAA-IOP considerably outperforms the dataset with only R_rs, which suggests that the application of QAA-IOP can enhance predictive accuracy and performance, especially in worse water quality or less transparency. Near Wenchang, as the number of training pairs increases to 1000, the GPR and NN algorithms training on the QAA-IOP dataset gradually perform better than that on the R_rs dataset. Notably, the GPR algorithm performs well even with only 1000 training pairs, indicating it may perform better under less data conditions. Near Laizhou Bay, when training pairs reach 500, the accuracy of all four QAA-IOP datasets begin to consistently surpass the four R_rs models. As illustrated in Figure 7C-III,D-III, near the Qilian Islands, the convergence rates of RF and SVR are significantly lower than the other two areas.

Overall, near the Wenchang area, all four models tend to converge around 3500 training pairs, whereas the SVR algorithm may continue to improve with further increments in training pairs. Models near Laizhou Bay tend to converge around 2000 training pairs. Models near the Qilian Islands tend to converge around 3000 training pairs, whereas the RF algorithm may continue to improve with further increments in training pairs. As a result, in subsequent experiments, 3500 pairs are randomly selected and used to train the models for Wenchang and the remaining 2600 pairs are allocated for validation. Near Laizhou Bay, 2000 pairs are randomly selected and used for training and 400 for validation. Meanwhile, near the Qilian Islands, 3000 pairs are selected for training and 600 for validation.

Table 3. Comparison of the model training and validation with the R_rs and QAA-IOP datasets.

Study Area	Training Algorithm	RMSE (m) (Training)	R2 (Training)	RMSE (m) (Validation)	R2 (Validation)
Wenchang	GPR with Rrs	0.85	0.90	0.92	0.88
	GPR with QAA-IOP	0.79	0.92	0.80	0.91
	NN with Rrs	0.88	0.90	0.91	0.88
	NN with QAA-IOP	0.80	0.91	0.85	0.91
	RF with Rrs	0.87	0.90	0.92	0.88
	RF with QAA-IOP	0.83	0.91	0.87	0.89
	SVR with Rrs	0.89	0.90	0.91	0.89
	SVR with QAA-IOP	0.89	0.89	0.79	0.91
Laizhou Bay	GPR with Rrs	0.49	0.71	0.54	0.71
	GPR with QAA-IOP	0.44	0.77	0.49	0.74
	NN with Rrs	0.51	0.68	0.55	0.70
	NN with QAA-IOP	0.44	0.78	0.48	0.77
	RF with Rrs	0.50	0.70	0.56	0.68
	RF with QAA-IOP	0.45	0.75	0.52	0.73
	SVR with Rrs	0.52	0.68	0.58	0.67
	SVR with QAA-IOP	0.48	0.74	0.51	0.72
Qilian Islands	GPR with Rrs	0.72	0.96	0.75	0.97
	NN with Rrs	0.73	0.96	0.74	0.96
	RF with Rrs	0.91	0.94	0.93	0.94
	SVR with Rrs	0.98	0.93	0.97	0.94

lists the performance metrics of the model training and bathymetry validation. Near Wenchang, the GPR, NN, and RF models using the QAA-IOP dataset evidently outperform the models training on the R_rs dataset. Specifically, the GPR and NN models exhibit better performances than the RF and SVR models. Near Laizhou Bay, all four machine learning methods utilizing the QAA-IOP dataset obtain better results than those utilizing the R_rs dataset. Notably, the NN model corresponds to the best performance with the least RMSE (0.44 m and 0.48 m) and largest R² (0.78 and 0.77) in training and validation, respectively. Near the Qilian Islands, GPR and NN models also exhibit better performances than the RF and SVR models.

In both Figure 7 and Table 3, considering QAA-IOP in Wenchang and Laizhou Bay can improve the model accuracy. Compared to the results in Wenchang, the QAA-IOP dataset has more improvements in all four models near Laizhou Bay, which indicates that the QAA-IOP is more effective in worse water quality or less transparency. In very clear waters, i.e., the Qilian Islands, the performance of the QAA-IOP dataset is identical to the R_rs dataset. Near Wenchang and Laizhou Bay, apart from differences in training and validation accuracy, models that use the R_rs dataset exhibit a slight tendency towards overfitting. For instance, in the GPR model in Wenchang, the training RMSE is 0.85, while the validation RMSE drops to 0.92. Models using the QAA-IOP dataset demonstrate slightly smaller score variations between training and validation sets, suggesting that QAA-IOP might be helpful for reducing overfitting risks. However, to evaluate overfitting risks more precisely, further experiments on larger datasets and detailed model tuning are required.

Table 3 demonstrates that the neural network (NN) model exhibits the best accuracy among the four algorithms, and thus the NN model is used to generate bathymetry results in the three study areas. Figure 8 presents a scatterplot of the bathymetry estimations from the NN model against the ICESat-2 validation bathymetric points. In Wenchang, in Figure 8A,B, the maximum depth measured by ICESat-2 is 16 m. Incorporating QAA-IOP parameters significantly reduces the RMSE. For depths less than 10 m, the scatter plots are closely clustered around the reference line of y = x. Due to the limited reference data within the 10–13 m range, a slight accuracy reduction is observed with the QAA-IOP method for these depths.

In Laizhou Bay, in Figure 8C,D, the maximum depth measured by ICESat-2 is only 5 m due to the relatively turbid waters. The worse water quality further introduces lower R² values than that in Wenchang on both the R_rs and QAA-IOP datasets. Compared to models trained on the R_rs datasets, the accuracy of models incorporating the QAA-IOP parameters is significantly enhanced, with the regression line almost coinciding with the reference line of y = x and error scatter uniformly distributed on both sides of y = x. In the Qilian Islands, in Figure 8E, the maximum depth measured by ICESat-2 is 20 m. Comparing the results in Figure 8A–E, the accuracy of training results on the R_rs dataset in very clear waters is better than that of other water qualities.

Figure 9 displays the bathymetry maps of Wenchang, Laizhou Bay, and the Qilian Islands generated using the NN models and the QAA-IOP dataset. The bathymetry result in Wenchang indicates that the underwater terrain has many sandy ridges. Near the Yellow River estuary, the bathymetry of Laizhou Bay exhibits a more gradual topographic change. The bathymetry map clearly shows the contours of ridges in the northern part of Laizhou Bay. Furthermore, the navigation channels between Harbor I and Harbor II are distinctly identified within the red dashed box on the Laizhou Bay bathymetry map. The bathymetry result of the Qilian Islands exhibits that the central part of the area is a relatively flat coral reef, with very steep edges around the coral reef. The MAPEs increase as the water qualities worsen in three study areas, i.e., 13.6% in Qilian Islands, 17.5% in Wenchang, and 21.6% in Laizhou Bay.

5. Discussion

5.1. Error Analysis

In this study, the accuracy of the SDB model is mainly influenced by Sentinel-2 imagery and “truth data” from ICESat-2. The aerosols and sun glint from the Sentinel-2 imagery contribute significantly to the errors in bathymetry [59]. During image preprocessing, the C2RCC algorithm effectively mitigates most of the interference from physical factors on spectral data [60,61]. However, in deeper water areas, uncertainties from atmospheric correction and residual errors from sun glint correction may still introduce discrepancies [62]. In the deeper water regions in Wenchang and the Qilian Islands, regularly shaped error fluctuations are observed, which may be caused by residual errors from sun glint correction. The QAA algorithm in optical shallow waters can be influenced by the high reflectivity of sea bottoms, which further leads to an erroneous estimation of the bottom reflectivity and ultimately manifests as an overestimation of the absorption coefficient [63,64].

Although most scattering and refraction errors in ICESat-2 depth measurements have been corrected, related studies suggest that an approximately 0.4 m bathymetric error still exists in depths less than 30 m [6]. This remaining error is due to the randomness or uncertainty of a photon occurrence, as the bottom signal is normally extremely weak, even less than 0.1 photon on average per shot. Specifically, this uncertainty increases when the width of the detection probability curve broadens, e.g., in more turbid waters or reflected from a bottom with a large slope. Additionally, the uneven distribution of training data pairs in different depth ranges introduces errors in the SDB model [65]. For instance, a notably scarce number of samples exist in Wenchang within the 10–13 m depth range, resulting in significant errors in Figure 8B within this range.

5.2. Variations in the Quantities of Information during Feature Selection Process

When employing the wrapper method to identify the optimal feature combination, multiple equivalent optimal combinations can be found. Taking Laizhou Bay as an example, the optimal combination utilized in this study comprises R_rs(665), R_rs(865), a(560), and a(665). However, alternative combinations, such as R_rs(665), R_rs(865), a(490), a(560), and b_bp(665), also achieve nearly identical accuracy levels. To interpret the relationship between the SDB model accuracy and input features, as depicted in Figure 10, we plot the information gains generated by all bands during the SDB process, which is the difference between a certain band and a column of random numbers in conditional entropy [66]. Additionally, we calculate the mutual information among all bands for each area, and the results show that the mutual information between all bands exceeded 0.8. This indicates a very high degree of information overlap among the features.

Based on the information gain and mutual information, it can be inferred that the quantities of information carried by different R_rs and IOPs bands vary. Consequently, a combination of bands, each carrying a lesser quantity of information, can substitute for a single band with a greater quantity of information. Regarding the previously mentioned equivalent combinations for the Laizhou area, the use of a combination of a(490) and b_bp(665) achieves the same information content as a(665), which leads to comparable SDB accuracy. Furthermore, for the example of a feature selection iteration process in Laizhou Bay, R_rs(665), R_rs(865), a(490), a(560), and b_bp(665) are selected in sequence. As new features are incorporated, the RMSE changed in the following sequence: 0.59, 0.55, 0.51, 0.49, and 0.47 m. It is evident that due to the excessive mutual information, the accuracy improvement experiences diminishing returns. When one more feature (such as a(443)) is further added in the model training, the RMSE even increases to 0.48 m, i.e., the performance gets worse. For machine learning models, reducing data dimensionality may be an effective way to enhance accuracy [43]. Therefore, identifying the most effective bands under different water quality conditions is crucial, rather than inputting all available information into the machine learning model.

As shown in Figure 10, for the Wenchang region that belongs to the Jerlov III water type, the diffuse attenuation coefficient (K_d) is primarily contributed by scattering from water and dissolved particulate matter [47]. Hence, incorporating the particulate backscattering coefficient, such as b_bp(490), can effectively improve the accuracy of SDB. In the case of Laizhou Bay, which belongs to the Jerlov 7C water type, the information gains from R_rs at visible wavelengths become significantly lower compared to IOPs. As a result, during the SDB process, the priority of features should be given to the bands that include the effects of Colored Dissolved Organic Matter (CDOM) and chlorophyll, such as a(490) and a(665). For the Qilian Islands, considering IOPs fail to achieve higher SDB accuracy than the R_rs dataset, which may be due to the very clear waters that contain little particles. For regions across all three water quality types, R_rs(856) provides a very high information gain, which might be that the near-infrared (NIR) band supplies essential information for compensating for residuals in the sun glint correction and atmospheric correction.

6. Conclusions

This study presents a new approach to SDB by integrating the QAA with data from the Sentinel-2 and ICESat-2 satellites. With the remote sensing reflectance R_rs from Sentinel-2 and IOPs from QAA as features, we employ ICESat-2 bathymetric data as training data to derive SDB machine learning models. After dimensionality reduction of the features using the wrapper method, we estimate water depths by employing four machine learning models and assess the bathymetric accuracy in three different study areas which correspond to different water qualities.

The results indicate that with the exception of very clear waters (Jerlov IB), incorporating QAA-IOP parameters can significantly enhance bathymetric accuracy, especially in turbid waters. By analyzing information gains in SDB, the most effective inputs are identified and prioritized under different water qualities. Also, incorporating an excessive number of redundant bands into machine learning models not only increases the demand of computing resources but also leads to worse accuracy in SDB. This research evaluates the importance of selecting optimal spectral bands for different water qualities and highlights the potential of machine learning in improving SDB models. Future work should focus on the robustness of these findings across larger datasets and varied coastal environments.