Bathymetry Retrieval From Spaceborne Multispectral Subsurface Reflectance

A few scholars have developed the models for retrieval of water depth from subsurface reflectance of multispectral images to avoid the influences of sun glitter. However, the models are only suitable for case I water. For this reason, this study proposes a bathymetry retrieval model using subsurface reflectance for both case I and case II water. The model first corrects the water surface reflectance image and then converts it into a subsurface reflectance image, and the subsurface reflectance image is used as the water depth retrieval image. Landsat 8 images were taken for experiments in case 1 water and case 2 water, and two water areas, Weizhou Island, Guangxi, China, and Molokai Island, Hawaii, USA, were used to verify the proposed model. The experimental results showed that the proposed model reduced the root-mean-squared error of the retrieved water depth in the Weizhou and Molokai areas from 3.113 to 2.903 m and 4.239 to 3.653 m, respectively, i.e., improve accuracy of water depth at 6.75% and 13.82% for Weizhou and Molokai areas, respectively. Therefore, the results demonstrate that the proposed model using subsurface reflectance can significantly improve the accuracy of bathymetry retrieval via spaceborne multispectral images.

Semiempirical bathymetric model: The bathymetry of satellite images is obtained by combining theoretical models and empirical parameters, which are divided into single-, dual-, and multiband models. Seawater irradiance or remote sensing reflectance is divided into two parts: deep water and bottom reflection. The single-band model only needs to regress two empirical parameters of the model to invert the water depth. Because only one band is utilized, it is required that the band has a relatively high transmissivity to water; therefore, the green band is generally used. The disadvantage of these models is that only one band of bathymetric information is utilized and, thus, information from the other bands is neglected, which increases the error of the method. Lyzenga [35], Paredes and Spero [36], and Spitzer and Dirks [37] proposed a two-band log-linear bathymetric model in which the optical properties of the water in the study area were assumed to be homogeneous and the reflectance of the two bands was considered constant over different substrate types, and a log-linear model was established in which the retrieval effect was much better than that of the single-band model. van Hengel and Spitzer [38] proposed that according to the above assumptions, the retrieval methods will have certain errors. If the water body under study has inhomogeneous water optics, then the suspended matter in the water body will be unstable and change with time, while the seabed topography will not change with time. Bathymetry retrieval can be performed using multitemporal satellite data with short intervals, and the stability and accuracy of bathymetric retrieval results can be improved by averaging the multitemporal retrieval results. Su et al. [39] found that the previous bathymetric retrieval studies used a single regression coefficient in the entire image, although in actual conditions, the type of substrate and water quality in each scene will vary with space; thus, these authors proposed a method of chunking the study area and assigning weights to the adjacent bathymetric reference points according to the This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ geographic location of the pixels in the geographic area or local area, with the points close to the centroid point having higher weights than those far from the centroid point, thus changing the bathymetric retrieval model parameters and solving the problem of inhomogeneous seafloor types and water quality. This method improves the reliability and accuracy of retrieval. Figueiredo et al. [40] suggested that the water depth retrieval model of the entire image should not use the same set of coefficients; rather, different retrieval coefficients should be used for different seabed types and water qualities so that the retrieval depth can be closer to the real water depth. However, they and Su et al. [39] adopted a different approach. Figueiredo et al. [40] performed Tikhonov regularization on the original model and tested the modified model in Lisbon, Portugal, and the results show that this method can improve the accuracy of the water depth retrieval.
Semianalytical bathymetric model: This model was proposed by Lee et al. [41] in 1999 and does not require actual water depth data, and it is currently the most widely used hyperspectral remote sensing bathymetric model. The model takes advantage of the high information capacity of the hyperspectral image to invert the inherent optical properties and water depth of the water body. Based on this model, Lee et al. [42] developed a quasi-analytical algorithm (QAA). These authors suggested that the water quality of different sea areas is very different and the accuracy of the simulated backscattering and absorption coefficients of water bodies is very low, which leads to large errors. Thus, the retrieval accuracy will be improved if the two are calculated directly. Compared with the results of field measurements, the retrieval results of the QAA algorithm were very similar to the real data. Li et al. [43] introduced semianalytic methods into the stump model to retrieve the water depth under the condition of case I water (Case I water is defined as a concentration of phytoplankton high compared with that of other particles in water [44]) in 2019. In this method, the parameters of the stump model are established using absorption characteristics of chlorophyll-a. Planet Dove satellite images were used as experiments image for verification of the method. The major disadvantages of this method are its uncertainty of the model parameters due to the different components under the various water environment, such as chlorophyll-a and seawater yellow substance, which leads to bathymetric retrieval error. Chen et al. [45] also presented the application of subsurface reflectance (the subsurface layer is a theorical concept, which is usually defined as just below the surface of the water) for retrieval of the water depth for case I water in 2019. The GeoEye-1 and Gaofen-2 images are used to verify the accuracy of the bathymetry by a semianalytic model. However, the semianalytic model is commonly applied in retrieval of bathymetry on the basis of the hyperspectral images. The multispectral images have a longer band width than the hyperspectral images have, and the central wavelength between the hyperspectral and multispectral images is not completely consistent, which also leads to bathymetric retrieval error.
As analyzed above, most of the models proposed by authors, such as Li et al. [43] and Chen et al. [45], only can be used in case I water for bathymetric retrieval. In fact, not all of water belong to case I water, i.e., many water is in case II water (Case II water not only has phytoplankton, has also inorganic particles in the water, and inorganic particles dominate in the water [44]). Therefore, this article proposes the development of a new multispectral-based bathymetric retrieval model using subsurface multispectral reflectance for both case I and case II water. The rest of this article is organized as follows. Section II describes the proposed bathymetric retrieval model that uses subsurface reflectance. Section III presents the experimental results and analysis. Finally, Section IV concludes this article. Fig. 1 explains how water surface reflection causes errors during the establishment of a multispectral-based bathymetric retrieval model. When sunlight reaches the water surface, it produces sun glitter on the water surface. The satellite receives reflectance radiation, which consists of both sun glitter and seafloor reflection information, resulting in an error in the multispectral water depth retrieval model. However, if the subsurface layer is used for the multispectral water depth retrieval model, this type of error should significantly be reduced.

II. BATHYMETRIC RETRIEVAL THROUGH SUBSURFACE REFLECTANCE
A Landsat 8 satellite image was taken as an example to explain the method of establishing the subsurface-based water depth retrieval model below.
First, because raw surface reflectance images contain noise that causes errors during conversion into the subsurface reflectance process, image noise must be spatially smoothed. For this reason, many scholars have suggested that the near-infrared Landsat 8 Band 5 should be used to smooth out other band images since water bodies have strong absorption at Band 5. This finding implies that Band 5 theoretically has a very low reflectance associated with low and weak noise. The same noises in all bands of Landsat 8 were considered to have approximately the same value. Therefore, the noise in the other bands can be removed by subtracting the image values from Band 5. Thus, smoothing of all raw images is performed by subtracting from where λ i is the ith band (i = 2, 3, 4), R raw rs (λ i ) is the raw surface reflectance image at Band i(i = 2, 3, 4), and R s rs (λ i ) is the surface reflectance image after smoothing the Band i (i = 2, 3,4). This means that images at Bands 2, 3, and 4 are smoothed using (1).
In addition, the raw image at Band 5 also contains noise. To obtain a "true" value of Band 5, the "noiseless" image at Band 5 can be calculated as follows [46]: where R n rs (λ 5 ) is the "noiseless" surface reflectance image at Band 5 and R s rs (λ 4 ) is the surface reflectance image after smoothing the image at Band 4.
After smoothing the noise of the images at Bands 2 and 3 using (1), the correction of surface reflectance for images at Bands 2 and 3 can be carried out as follows [47]: The surface reflectance R c rs (λ i ) carries water depth information, including the optical properties of the water body and bottom reflectance. According to Lee et al. [48], the relationship between the surface reflectance and subsurface reflectance can be expressed by the following equation: where r rs (λ i ) is the subsurface reflectance image of Band i (i = 2, 3). Rewriting (4), the subsurface reflectance can be obtained by where the symbols are the same as those above. With (5) and the semiempirical model proposed by Zhou [34], water depth retrieval based on subsurface reflectance can be expressed as follows: where r sh rs (λ i )(i = 2, 3) is the subsurface reflectance in the shallow water of Band i (i = 2, 3), r dp rs (λ i )(i = 2, 3) are the minimum values of subsurface reflectance in the deep water area in Band i (i = 2, 3) image, the deep water area is selected by tailoring from a big study area away from the coast under the same water quality environment and n is the magnification factor, which is set to a value of 10 000 for Landsat 8 images.

III. EXPERIMENTS AND VALIDATION
A. Two Study Areas 1) Study Area 1: Study Area 1 is located on the west side of Weizhou Island, Guangxi, China, at geographic coordinates from 21°01 47.18 N to 21°03 25.45 N and 109°03 22.46 E to 109°05 9.7 E. The study area is a square water area of 3.06 km×3.06 km (see Fig. 2). There are four main types of substrate, with black volcanic bedrock, dark magnesian sand, light-colored felsic sand, and coral in the south. The maximum water depth in this area was 22.96 m, and the average water depth was 14.45 m. The study area is case II water.
The RESON Seabat 7125 multibeam bathymetric system was used to measure the bathymetry over an area of approximately 24 km 2 on western Weizhou Island (see Fig. 3) in November 2019. The data collected by the multibeam were processed into point data using PDS2000 and HIPS&SIPS software and then  processed and cropped to obtain bathymetry data of the study area using ArcGIS (see Fig. 4). The accuracies of the measured water depth and positioning are listed in Table I. The error of the sounding points is listed in Table II [49]. The bathymetric data have been made the tide correction.
2) Study Area 2: Study Area 2 was selected from a sea area on the west side of Molokai Island, Hawaii, USA, with geographic coordinates from 21°08 37.28 N to 21°11 20.12 N and 157°14 52.91 W to 157°18 11.34 W, as shown in Fig. 5. The maximum water depth of this area is 37.39 m, and the average depth is 18.82 m. Molokai Island is the least polluted island in the Hawaiian Islands, and its water quality is clear; thus, it is an excellent water area for bathymetry retrieval by satellite image. The study area is case I water.
The CZMIL airborne LiDAR was used to conduct the bathymetry in a shallow water area around Molokai Island in October 2013. The specifications of the CZMIL are listed in   Table III [50]. In this study, a part of the west coast of Molokai Island was selected as the experimental area (see Fig. 5), and the bathymetric data obtained from the CZMIL LiDAR system were transformed into DEM data using ArcGIS interpolation method with a resolution of 30 m (see Fig. 6). The bathymetric data have been made the tide correction.
With Tables I and III, although the accuracy of the bathymetric data obtained by CZMIL LiDAR System is higher than that obtained by the multibeam measurements, the accuracy for the two datasets can meet the requirement for the verification of the bathymetric data retrieved by multispectral satellite images in this article.
3) Ten Known Reference Points (KRPs): The retrieval of water depth using the model in this study requires the selection of ten points at different water depths that are uniformly distributed throughout the region, which are referred to as KRPs, from the "true" water depth data. Thus, ten KRPs were selected in both  study areas in this study, and the distribution of the ten KRPs is depicted in Figs. 5 and 7.

4) Check Points and Check Lines:
To verify the effect of the retrieved bathymetry, ten check points distributed at various water depths were selected in this study (see Figs. 5 and 7). Their values are shown in Table IV (Study Area 1) and V (Study Area 2). In addition, four bathymetric check lines marked Lines 1 to 4 were selected in Study Area 1 (see Fig. 4), and three bathymetric check lines, marked Lines 1-3, were selected in Study Area 2 (see Fig. 6), which were evenly distributed over the study areas.   Image preprocessing included atmospheric correction using the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) Model. Then, the region of interest within the test field was used to clip the image corresponding to the same size of the test field, as shown in the red boxes of Figs. 3 and 6.

1) Experiments in Study Area 1: Weizhou Island:
The bathymetric retrieval using the proposed model (6) was performed using the following steps. The first step is to calculate the required surface reflectance. The correction of the surface reflectance at Bands 2 and 3 was obtained using (1)- (3). The results are shown in Fig. 7(a) for Band 3 and Fig. 7(b) for Band 2.
The second step is to calculate the subsurface reflectance r rs , which can be obtained using (5). The results are shown in Fig. 8(a) for Band 3 and Fig. 8(b) for Band 2.   The difference value between the subsurface reflectance and the raw surface reflectance of Study Area 1 is shown in Fig. 9(a) for Band 3 and Fig. 9(b) for Band 2.
Finally, the water depth is calculated. By substituting the calculated subsurface reflectance into (6), the parameters (a 0 and a i (i = 2, 3)) required in (6) were regressed according to the reference points selected in Fig. 4. The water depth retrieved using the proposed model for the study area is shown in Fig. 13(a).
2) Experiments in Study Area 2: Molokai Island: The bathymetric retrieval using the proposed model (6) was performed using the following steps. The first step is to calculate the required surface reflectance. The correction of the surface reflectance at Bands 2 and 3 was obtained using (1)- (3). The results are shown in Fig. 10(a) for Band 3 and Fig. 10(b) for Band 2.  The second step is to calculate the subsurface reflectance, which can be obtained using (5). The results are shown in Fig. 11(a) for Band 3 and Fig. 11(b) for Band 2.
The difference value between the subsurface reflectance and the raw surface reflectance of Study Area 2 is shown in Fig. 12(a) for Band 3 and Fig. 12(b) for Band 2.
Finally, water depth is calculated. By substituting the calculated subsurface reflectance into (6), the parameters [a 0 and a i (i = 2, 3)] required in (6) were regressed according to the reference points selected in Fig. 6. The water depth retrieved using the proposed model for the study area is shown in Fig.  16(a).

1) Study Area 1. Weizhou Island Water Areas:
a) Accuracy evaluation of the bathymetry results using the deepest water depth data: To compare the accuracy of the bathymetry results using the proposed model, the conventional model presented by Lyzenga (1985) was used to retrieve the water depth. The results are shown in Fig. 13(b).
As observed in Fig. 13, the deepest water depths retrieved by the conventional model and the proposed model [i.e., (6)   uses subsurface reflectance that can effectively reduce the sun glitter interference and ensure high accuracy of retrieval water depth. With the resulted data, the contour lines in Fig. 13(a) are relatively smooth. In contrast, while the surface reflectance applied by the conventional model cannot reach this effect. Thus, the proposed model can significantly improve the accuracy of water depth measurements. b) Accuracy evaluation of the bathymetry results using checkpoints: The ten check points on Weizhou Island are shown in Table IV. Table V presents a comparison of the error at the checkpoints between the two models. The maximum absolute error from the conventional model was 5.14 m, whereas that from the proposed model was 4.42 m relative to the "true" water depth. The maximum absolute error occurs at check point IV, where the "true" water depth is 12.6 m, while the retrieved water depth is 17.02 m. Therefore, the maximum error in the water depth from the proposed model decreased by 14.01% relative to that of the conventional model. The mean error of the water depths retrieved by the proposed model was 2.20 m, while that for the conventional model was 4.84 m. Therefore, the mean error of the water depth from the proposed model decreased by 33.99% relative to that of the conventional model. The root-mean-squared error (RMSE) relative to "true" water depth by the proposed model was 2.16 m, while that from the conventional model was 2.84 m. Therefore, the RMSE of the water depth from the proposed model decreased by 23.94% relative to that of the conventional model. Therefore, the accuracy of the bathymetry results retrieved by the proposed model was improved by 23.94% relative to that of the conventional model through the verification of ten check points.  c) Accuracy evaluation of the bathymetry results using check lines: Fig. 14(a)-(d) shows the profiles of Lines 1 to 4 in Fig. 4. The topography obtained by the two models in the same profile clearly shows that the trends of the two are roughly the same, with a greater water depth on the right side and a smaller water depth on the left side, both conforming to the shape of the true seafloor. In most cases, the subsurface reflectance model is much closer to the "true" water depth. The proposed model reduces the retrieved water depth error at most by 3.2 m at Line 1; 3.8 m at Line 2; 3.1 m at Line 3; and 5.6 m at Line 4 compared with the conventional model. In the topography of the retrieval, the seafloor terrain obtained by the conventional model forms a jagged shape, and the depth difference between the retrieval of two adjacent pixels is too large, as shown in Fig. 14, where the maximum depth difference between two adjacent pixels can reach 5 m. Compared with the seafloor terrain in the retrieval of the subsurface reflectance model and conventional model, the fluctuation is smaller and the depth of the retrieval of adjacent pixels is closer to the real seafloor terrain. The depth of adjacent pixels is closer to the real seafloor topography, which effectively improves the retrieval stability.
The above analysis shows that the use of the model proposed in this article can effectively reduce the errors caused by the water surface. This means that the model proposed in this study can achieve higher accuracy than the conventional model does in practical applications. d) Accuracy evaluation of the bathymetry results using scatter plots: To further compare the accuracy distributions of water depth measurement, the scatter plots of all the points from the proposed model and the conventional model with respect to the bathymetry measured ("true" value) using multibeam device at Study Area 1, Weizhou Island, are depicted in Fig. 15. As observed from Fig. 15, it can be found that the proposed model has stronger compactness than the conventional model does. Moreover, when the water depth was greater than 15 m, the correlation between the water depth obtained by the proposed model and the "true" value was much stronger than that of the conventional model does. This implies that the water surface had a big influence on the surface reflectance observed by the satellite, which shifted the bottom reflection information when passing through the water surface. By correcting the water surface, the surface reflectance was converted into subsurface reflectance and the influence of the water surface on the outgoing sunlight was weakened or even eliminated, thus achieving the purpose of improving the accuracy of the bathymetry retrieval.
In addition, at a water depth range of 10-15 m, both the proposed model and the conventional model had certain deviations. The conventional model had a large retrieval range, strong fluctuation, and low robustness, whereas the proposed model appeared to overestimate the water depth but had strong robustness. The proposed model replaced the surface reflectance with subsurface reflectance and eliminated the fluctuations generated at the water surface in the water depth range of 10-15 m, which increased the stability of the retrieval results. The overestimation of the water depth at 10-15 m in this model was caused by the linear function used in the fitting. As observed from Figs. 15 and 18, when the water depth is less than 5 m, the proposed model and the conventional model are dispersion (see Part A in Fig. 15 and Part A in Fig. 18). When the water depth reaches around 10 m, the proposed model and the conventional model become not dispersion (see Part B in Fig. 15 and Part B in Fig. 18). When the water depth reaches around 15 m and more, the proposed model and conventional model are become aggregation (see Part C in    Fig. 18). This phenomenon demonstrated that it is common for bathymetric retrieval to produce dispersion in shallow waters (<5 m).
However, in Weizhou Island, it is slightly dispersion relative to that in Molokai Island. It may have been due to the fact that the influence of seafloor reflectivity on pixels becomes larger in shallow water areas and the reflective from seafloor substrate is strong and the turbidity of the seawater.
In addition, the R 2 and RMSE from the conventional model achieved 0.63 and 3.113 m, and 0.68 and 2.903 m from the proposed model. The RMSE from the proposed model relative to the "true" value was 0.21 m, which is lower than that from the conventional model relative to the "true" water depth. This means that the accuracy increases 6.75%.
2) Study Area 2. Molokai Island Water Areas: a) Accuracy evaluation of the bathymetry results using the deepest water depth data: To compare the accuracy of bathymetry using the proposed model, the conventional model presented by Lyzenga (1985) was used to retrieve the water depth, and the results are shown in Fig. 16(b).
As observed from Fig. 16, the deepest water depth retrieved by the conventional model and the proposed model (i.e., (6)) in Molokai Island was approximate −33.71 m and −35.90 m, respectively (see the symbol, "star" in Fig. 16). The deepest water depth was −37.39 m. Therefore, the errors from the conventional model and the proposed model were 3.68 m and 1.49 m, respectively. Thus, the proposed model can significantly improve the accuracy of water depth measurements using spaceborne multispectral retrieval.
b) Accuracy evaluation of the bathymetry results using checkpoints: The ten check points on Molokai Island are shown in Table VI. Table VII presents a comparison of the error at the checkpoints between the two models. The maximum absolute error of the conventional model was 7.50 m, whereas the maximum absolute error of the proposed model was 5.16 m. Therefore, the maximum error in the water depth from the proposed model decreased by 31.20% relative to that of the conventional model. The mean error of the water depths retrieved by the proposed model was 2.81 m, while that of the conventional model was 3.88 m. Therefore, the mean error of the water depth from the proposed model decreased by 27.58% relative to that of the conventional model. The RMSE relative to the "true" water depth by the proposed model was 3.21 m, while it was 4.24 m  Fig. 6 are shown in Fig. 17, and they were evenly distributed in the experimental area. The topography obtained by both models in the same water depth profile clearly showed that the trends of both were roughly the same and the water depth gradually increased from left to right, which was the same trend as the results of LiDAR scanning. The proposed model reduced the retrieved water depth error at most by 3.1 m at Line 1; 2.1 m at Line 2; and 3.4 m at Line 3 compared with the conventional model. In most cases, the proposed model was closer to the "true" water depth. d) Accuracy evaluation of the bathymetry results using scatter plots: To further compare the accuracy distribution of bathymetry in Study Area 2, Molokai area, Fig. 18 depicts the scatter plot of all points of the proposed model and the conventional model relative to the bathymetry ("true" values) measured using the Lidar device in Study Area 2. As observed from Fig. 18, the proposed model was more compact. The water depth obtained by the proposed model was overall closer to the "true" water depth in the range of 10-30 m, and the proposed model was more accurate than the conventional model and showed a certain degree of improvement for the underestimated water depth of conventional when the water depth was greater than 30 m, but the error of retrieving water depths above 30 m becomes larger. From the experimental results, it can be discovered that the proposed model in case I can measure water depths greater than 30 m. The similar results can also be demonstrated by Ceyhun and Yalçın [51] and Poursanidis et al. [52]. Overall, the proposed model performed better than the conventional model.
In addition, the R 2 and RMSE of the conventional model were 0.84 and 4.239 m, and the R 2 and RMSE of the proposed model were 0.87 and 3.653 m (see Fig. 18), respectively. The RMSE from the proposed model was 0.586 m lower than that from the conventional model relative to the "true" water depth. This means that the accuracy from the proposed model increases 13.82% when compared with the conventional model. This is because the application of subsurface reflectance can reduce the information that interferes with bathymetric retrieval, such as surface sun glitter, so the accuracy of bathymetric retrieval can be improved.
The above analysis indicates that the proposed model produces smaller errors and improves the accuracy of the bathymetric retrieval of satellite images.
And compared with the scatter plot of Weizhou Island, the scatter plot of Molokai Island is more compact overall, which is because the water quality of Molokai Island is better than that of Weizhou Island, and the better water quality can ensure the stability of the retrieval.

IV. CONCLUSION
Previous models for retrieval of water depth from subsurface reflectance of multispectral images are used to avoid the influences of sun glitter. However, the models are only suitable for case I water. For this reason, this study proposes a bathymetry retrieval model using subsurface reflectance for both case I and case II water. The innovative idea is based on the fact that sunlight passing through the water surface produces a certain bias; however, available correction methods do not correct this bias.
Two experimental areas located on Weizhou Island (case II water) in Guangxi and Molokai Island (case I water) in Hawaii were used to verify the accuracy achievable using the proposed model. The RMSE of the water depth retrieved by the model proposed in this article was reduced by 0.21 m and 0.586 m relative to the "true" water depth measured by a multibeam device and an airborne LiDAR in the two study areas. The accuracy of the water depth retrieved by the proposed model in this study was 6.75% and 13.82% higher relative to the conventional model (Lyzenga model) in the two study areas, respectively. Therefore, bathymetric retrievals using the proposed model based on Landsat 8 data can significantly increase the accuracy of water depth measurements.