Bathymetry of the Coral Reefs of Weizhou Island Based on Multispectral Satellite Images

Abstract

Shallow water depth measurements using multispectral images are crucial for marine surveying and mapping. At present, relevant studies either depend on the use of other auxiliary data (such as field water depths or water column data) or contain too many unknown variables, thus making these studies suitable only for images that contain enough visible wavebands. To solve this problem, a Quasi-Analytical Algorithm (QAA) approach is proposed in this paper for estimating the water depths around Weizhou Island by developing a QAA to estimate the diffuse attenuation coefficients and simplifying the parameterization of the bathymetric model. The approach contains an initialization sub-approach and a novel global adjustment sub-approach. It is not only independent of other auxiliary data but also greatly reduces the number of unknowns. Experimental results finally demonstrated that the Root Mean Square Errors (RMSEs) were 1.01 m and 0.77 m for the ZY-3 image and the WorldView-3 (WV-3) image, respectively, so the approach is competitive to other QAA bathymetric methods. Besides, the global adjustment sub-approach was also seen to be superior to common smoothing filters: if the Signal-to-Noise Ratio (SNR) is as low as 42, i.e., ZY-3, it can smooth the water depths and improve the accuracies, otherwise can avoid the over-smoothing of water depths.

1. Introduction

Shallow seawater depth measurements are considered to be very important, as they are directly relevant to the applications of environmental management, exploration, and research [1,2,3,4]. For example, the natural development of coral reefs can strongly influence the physical structure of their environment. Information associated with water depth is fundamental for characterizing and discriminating between different kinds of coral reef habitats, such as patch reefs, spur-and-groove systems around reef fronts, and sea grass beds [1]. The knowledge of the water depth can also facilitate estimations of the bottom albedo, which can improve the habitat mapping [5]. Therefore, shallow water depth measurements have consistently represented an important component of marine surveying and mapping.

The bathymetry of the seafloors has usually been determined using conventional echosounders [6,7,8,9]. However, as the echosounders are usually fixed in ships [7,10], such echosounder methods are believed to be expensive and time-consuming, especially for shallow waters [3,11]. One reason is that the scan width of the echosounders will become narrow as the water becomes shallow. Another reason is that some remote and dangerous areas are not suitable for echosounder methods, as massive hidden reefs can make such areas unreachable to ships.

Hence, we focus on developing an approach to improve the performances of the multispectral satellite image-based bathymetric methods, since it is currently easy to collect multispectral satellite images. Besides, if multispectral satellite images are applied, then we no longer need to ship to the locations where the water depths will be measured. To test the approach, we chose the coral fringing reefs of the Weizhou Island in the Beibu Gulf, northern South China Sea, as our study area.

Specifically, Weizhou Island is located at 21°00′–21°05′N, 109°00′–109°10′E, approximately 36 miles from Beihai, Guangxi, as shown in Figure 1. Its coverage area of Weizhou Island is approximately 25 km². According to Wang et al. [12], the annual average Sea Surface Temperature (SST) and salinity values around Weizhou Island are about 24.55 °C and 31.9‰, respectively, thus it is ideal for coral growth even though the seawater transparency varies only from 3.0 m to 10.0 m [12]. Coral reefs have developed around Weizhou Island since the mid-Holocene, and these reefs cover an area of approximately 6–8 km² [13,14]. These coral reefs have long been the focus of studies analyzing the responses of coral to global warming and human disturbances, because they are located in a relatively high-latitude area and are heavily influenced by anthropogenic activities [12,14,15,16].

Within this study area, we expect to improve remote sensing data-based solutions for multispectral images to attain relatively accurate and reliable water depth estimations that may be applied to remote, large, or dangerous areas in the future. In addition, another practical motivation of this paper is to provide bathymetric results that will enhance future research of the effects of ambient environmental conditions on the coral reef of Weizhou Island itself.

2. Related Works

In this paper, we mainly focus on band combination methods and semi-analytic methods to make the proposed approach as effective, practical, and convenient as possible. Some related works are summarized as follows.

Based on the law that light is attenuated exponentially with depth in the water column, as follows:

$$L_i=L_{si}+k_i{r}_{Bi}\exp (-{\kappa }_{i}fz)$$

(1)

Lyzenga et al. [17,18] argued that the ratio of bottom reflectance between two bands is constant over a given scene for all bottom types; they demonstrated that errors caused by the presence of different bottom types could be corrected with two bands. Accordingly, the two band and n-band combination algorithms were then proposed for the water property-independent and bottom-reflectance-independent water depth extractions. Similarly, Clark et al. [19] introduced a linear multiband method for shallow water depth extraction, which was tested using bands 1 and 2 of Landsat TM imagery acquired in the vicinity of Isla de Vieques. To further develop an empirical solution for the presence of variable bottom types, Stumpf et al. [1] improved the linear multiband combination methods with a ratio of reflectances with only two tunable parameters.

Although band combination methods have been in use since the early 1970s, until now, they are still found to be useful and applicable because of their simple forms and efficient characteristics. Doxani et al. [20] evaluated the effectiveness of several band combination bathymetry algorithms via using a high spatial and spectral resolutions of WorldView-2 (WV-2) image. Yuzugullu and Aksoy [21] compared the linear regression model with a non-linear regression model to predict the depths in Lake Eymir based on WV-2 multispectral satellite images; their results showed that the non-linear regression model predicted these depths slightly better than the linear model. Papadopoulou et al. [3] used Lyzenga’s band combination methods [17,18] and high-resolution IKONOS-2 image to create digital bathymetric maps of the coastal area of Nea Michaniona, Thessaloniki, in northern Greece. Mohamed et al. [22] improved the linear multiband method by using a fitting algorithm called Ensemble Learning (EL) and obtained better performance and accuracy than those obtained by conventional PCA and generalized linear models. Many other related studies [2,23,24,25] were also found to focus on the improvements, comparisons, and applications of these algorithms.

Even though band combination methods are simple enough to have been widely applied by many previous studies, they are obviously highly dependent on field data, which are used for the derivation of the model coefficients. As a result, semi-analytic models were further developed and applied to overcome the disadvantages of the band combination methods by considering such factors as the Inherent Optical Properties (IOPs) of water, bottom types, and in-water light transfer.

Assuming that the water quality and bottom substratum are homogeneous, Bierwirth [26] outlined a method that unmixed the exponential influence of depth in each pixel by employing a mathematical constraint, which was applied in the analysis of Landsat TM data from Hamelin Pool in Shark Bay, Western Australia. Kutser et al. [27] have proved that the bottom type and water depth of the Australian Great Barrier Reef (GBR) can be simultaneously mapped by matching the Hyperion at-sensor radiance data to the simulated at-sensor spectral library. Lafon et al. [28] proposed a semi-empirical bathymetric methodology based on the Hydrolight Radiative Transfer code, which was calibrated using in situ bottom reflectance and effective attenuation coefficient measurements and then applied to determine the depth from SPOT images and to develop topographic maps of the tidal inlet of Arcachon. Lee et al. [29,30,31] also developed a semi-analytical method for shallow water depth estimation using an inversion-optimization approach, which can simultaneously derive the water depth and water column properties from the hyperspectral data in coastal waters. Brando et al. [4] further enhanced the physics-based inversion-optimization approach that was developed by Lee et al. [29,30,31]. They further considered the concentrations of optically active constituents in the water column, different types of substratum covers, and the contribution of every substratum to the remote sensing signal in their enhanced approach. With Lee’s semi-analytic forward model for shallow water remote sensing reflectance, Hedley et al. [32] was able to map water depths efficiently via the introduction of an Adaptive Look-Up Tree (ALUT). Knudby et al. [11] further applied this ALUT method to Landsat 8 data for producing the bathymetric maps for shallow waters in Canada. Recently, Hedley et al. [33] further presented a method for mapping the water depths and the leaf area index (LAI) of seagrasses according to the idea that variations in reflectance is caused by leaf length, leaf position, sediment coverage on leaves, water depth, and solar zenith angle, etc. Sylvain et al. [34] also presented a novel statistical semi-analytic method for mapping water depth and water quality using hyperspectral remote-sensing data, which is based on local maximum likelihood (ML) estimation performed on subimages. The authors assumed that deep water bottom is locally flat so that the spatial correlation between neighboring pixels can be considered sufficiently during the water depth estimation. This method is able to provide multi-resolution maps, where the resolution depends on local water depth.

Semi-analytic methods are usually regarded to be independent of field water depth measurements and having a wider scope of applications than band combination methods. However, many semi-analytic methods were originally developed for high resolution hyperspectral images, e.g., when Dekker et al. [35] reviewed and compared the performances of an empirical and five different analytic/semi-analytic methods, airborne hyperspectral data (Ocean PHILLS and CASI-2) was applied in the experiments. Thus, it is still currently difficult to declare that semi-analytic methods are able to completely solve all bathymetric problems, including those associated with the remote sensing images themselves, especially for the multispectral satellite images that are currently being collected.

For example, according to a radiative transfer model that can connects the water-leaving reflectance of the optically deep waters $r_{\infty }$, the bottom reflectance $r_b$, the diffuse attenuation coefficient of the water column $K_d$, and the water depth $H$ to the water-leaving reflectance $r_w$ as follows [36,37]:

$$r_w(\lambda )=r_{\infty }(\lambda )\left\{1-\exp \left[-2K_d(\lambda )H\right]\right\}+r_b(\lambda )\exp (-2K_d(\lambda )H),$$

(2)

Hu et al. [37] developed a semi-analytic method to derive the coastal water depths of the South China Sea islands using Landsat TM and ETM+ images. In this method, the diffuse attenuation coefficient $k_d(\lambda )$ of the water column is estimated from the adjacent optically deep waters by using the empirical model of MODIS or Landsat data for Case I water, as shown in Figure 2.

On the one hand, when this approach is applied to other types of satellite images, another empirical model need to be built for the estimation of the diffuse attenuation coefficient by using the optically deep water pixels selected from the surrounding areas. Otherwise, when other auxiliary datasets, e.g., MODIS, are needed for the estimation of the diffuse attenuation coefficient, we must ensure that the gap between the acquisition time of the multispectral satellite image and that of the auxiliary data is small enough. On the other hand, this model ignores the influence of the subsurface solar zenith angle ${\theta }_s$, the subsurface view angle ${\theta }_v$, the variations of the bottom reflectances, as well as the differences between the diffuse attenuation coefficients for downwelling and upwelling light.

A more accurate and typical radiative transfer model was developed by Lee et al. [29]. The expression constructs the relationships between the subsurface remote-sensing reflectance for optically deep waters ${\rho }_{\infty }$, the bottom irradiance reflectance ${\rho }_b$, the water depth $H$, and the subsurface remote-sensing reflectance $\rho$ as follows:

$$\rho ={\rho }_{\infty }\left\{1-\exp \left[-(\frac{1}{\cos {\theta }_s}+\frac{D_u^c}{\cos {\theta }_v})K_{d}H\right]\right\}+\frac{{\rho }_b}{\pi }\exp \left[-(\frac{1}{\cos {\theta }_s}+\frac{D_u^b}{\cos {\theta }_v})K_{d}H\right],$$

(3)

where $D_u^c$ and $D_u^b$ are the optical path-elongation factors for scattered photons from the water column and the bottom, respectively, $K_d$ is the diffuse attenuation coefficient. Some important semi-analytic bathymetric methods have been developed mainly based on the optimization of this model [30,38].

For example, the semi-analytic bathymetric model developed by Eugenio et al. [38] can be summarized as shown in Figure 3. Specifically, the optimum values of $a_{\varphi }(440)$, $a_g(440)$, $b_{bp}(400)$, $Y$, and $H$ are determined by minimizing the differences between the observed subsurface remote-sensing reflectance and the estimated reflectance, where $a_{\varphi }(440)$ is the attenuation coefficient for phytoplankton at 440 nm, $a_g(440)$ is the attenuation coefficient for gelbstoff at 440 nm, $b_{bp}(400)$ is the backscattering coefficient for particles at 400 nm, $Y$ is the spectral shape that depends on the particulate shape and size, $B$ is the bottom albedo at 550 nm, and $H$ is the water depth. This approach requires researchers to first initialize the $H$ parameter by applying the ratio algorithm, and then, to solve the nonlinear Model (3) using the Levenberg–Marquardt method.

Because the ratio algorithm is used for the initialization, this approach can barely function without the incorporation of field water depth data. In addition, many multispectral satellite images can provide only 3–4 visible wavebands, except for a few types of multispectral satellites (such as WV-2/3 can provide 6 visible wavebands). Hence, the model can become ill-posed because it contains six unknowns for each image pixel, which is usually more than the number of the visible wavebands.

To meet the requirements of most current multispectral satellite images, this paper develops a bathymetric approach on the basis of Model (3). In contrast to the band combination methods and the previously proposed semi-analytic methods, this study mainly highlights the following aspects: (a) making the approach independent of other auxiliary data, such as field water depth measurements and the production of MODIS diffuse attenuation coefficients; (b) reducing the number of unknowns so that the approach can satisfy most current multispectral satellite images; and (c) refining the estimated global water depths to smooth the final water depths and improve the accuracies. The specific details of this approach are discussed in the following sections.

3. Materials and Methods

3.1. Overview of the Approach

In this study, a Quasi-Analytical Algorithm (QAA) is firstly applied to some adjacent optically deep waters to estimate the diffuse attenuation coefficient $K_d(\lambda )$ and the optical path-elongation factors $D_u^c$ and $D_u^b$. As a result, the number of unknowns in the semi-analytic Model (3) is reduced from six ($a_{\varphi }(440)$, $a_g(440)$, $b_{bp}(400)$, $Y$, $B$, and $H$) to two ($B$ and $H$). Then, the shallow water depths can be further estimated by using an optimization without the diffuse attenuation coefficients being empirically estimated at first. This choice can also make the method completely independent of field water depth measurement data. Thereafter, we further attempt to propose a global adjustment method to refine the base of the estimated water depths depending on the following two conditions: (a) most of the water depths observed between a pixel and its adjacent pixels should be similar to each other; and (b) if the resolution of the image is high (being better than 5.8 m for the experiments), most of the water depths of the pixels that are close to the transient coastline should be approximately zero. In previous studies, these global adjustments have rarely been taken into consideration. Figure 4 shows the framework of the approach in detail.

3.2. Materials and Preprocessing

3.2.1. Experimental Data

Our experimental data consists of ZY-3 and WorldView-3 (WV-3) multispectral satellite images with 4 and 8 bands, respectively, as shown in Figure 5.

The ZY-3 satellite was launched on 9 January 2012 and represents the first civil high-resolution optical transmission surveying and mapping satellite of China. It is equipped with four optical cameras, including a nadir panchromatic Time Delay and Integration (TDI) CCD camera with a ground resolution of 2.1 m, two back/forward sight panchromatic TDI CCD cameras with a ground resolution of 3.6 m, and a nadir multispectral camera with a ground resolution of 5.8 m. The ZY-3 images are geo-referenced in the World Geodetic System (WGS-84). The geo-positioning accuracy is higher than 14.0 m without ground control after undergoing the on-orbit geometric calibration [39]. In addition, the images have a signal-to-noise ratio (SNR) of approximately 42 [40,41]. The ZY-3 multispectral image chosen as one of the remote sensing data sources for these experiments was acquired on 24 August 2015, at a corresponding tide height of 3.660 m. The satellite azimuth angle and zenith angle are 17.0° and 6.6°, respectively. The sensor band information of the multispectral sensor is presented in

Table 1. Band information of the ZY-3 multispectral satellite image.

Bands	Wavelength (nm)	Radiometric Calibration Coefficients ($\mathbf{W}·{\mathbf{m}}^{-2}·\mathbf{s}{\mathbf{r}}^{-1}·\mathbf{\mu }{\mathbf{m}}^{-1}$)
		Gain	Offset
Blue	450–520	0.233	0
Green	520–590	0.2162	0
Red	630–690	0.1789	0
Near-Infrared	770–890	0.1949	0

.

The WV-3 satellite was launched by DigitalGlobe Inc., Longmont, CO, USA on 13 August 2014. The WV-3 satellite sensor is a high-resolution commercial satellite sensor operating at an altitude of 617 km. It provides panchromatic, multispectral, short-wave infrared, and CAVIS images with resolutions of 0.31 m, 1.24 m, 3.7 m, and 30 m, respectively. The WV-3 multispectral image used in our experiments was acquired on 1 October 2015, at a corresponding tide height of 2.305 m. The satellite azimuth and zenith values are 218.3° and 21.1°, respectively. The geo-positioning accuracy is predicted to be higher than 3.5 m CE90 without ground control [42]. The WV-3 multispectral images consist of 8 different wave bands (the Coastal, Blue, Green, Yellow, Red, RedEdge, Near-IR1, and Near-IR2 wave bands). The wavelength and the radiometric calibration coefficient information can be all found the IMD file of the WV-3 data.

Moreover, when we performed ecological investigations of the coral reefs around Weizhou Island in October 2015, we measured the water depths at 45 points around the survey stations, as shown in Figure 6. The water depths were measured using an SM-5A hand-held portable depth sounder with an accuracy of 0.1 m, and the horizontal positions were mapped using a Magellan eXplorist 610 hand-held GPS navigation instrument with an accuracy of 3 m. Therefore, the measured water depths can be matched to the satellite data by using the measured horizontal positions and the geo-referenced information of the images.

In this study, these water depths were used to determine the accuracy of the proposed bathymetric approach after the influence of the tidal levels was eliminated. More specifically, all the measured water depths are subtracted by the corresponding tide heights before the accuracy assessment, as well as the water depths estimated from the experimental images. However, when the water depths estimated from the two experimental images are compared with each other, we selected the tide level of the ZY-3 image as a reference, and then corrected water depths estimated from the WV-3 image to this reference. As the tide height corresponding to the ZY-3 image is greater than the one corresponding to the WV-3 image, negative water depths can be then avoided in the comparisons by choosing such a reference.

3.2.2. Data Preprocessing

To implement the bathymetric approach, four main preprocessing steps of the multispectral imagery must be performed prior to obtaining the water depth estimation; these steps include the radiometric calibration, the atmospheric correction, the sunglint correction, and the subsurface remote sensing reflectance estimation.

Radiometric Calibration

The radiometric calibration is performed to transform the Digital Numbers (DN) into the Top-Of-Atmosphere (TOA) radiance $L$.

For the WV-3 image, this calibration is performed with the following equation:

$$L=Gain·DN·\frac{ACF}{EBW}+Offset,$$

(4)

where $ACF$ and $EBW$ are the abscalfactor and the effective bandwidth from the IMD file of the WV-3 data, respectively, and $Gain$ and $Offset$ are both provided by data from the DigitalGlobe website [43].

For the ZY-3 image, the calibration is easy to implement by using the radiometric calibration coefficients shown in Table 1 as follows:

$$L=Gain·DN+Offset.$$

(5)

Atmospheric Correction

The atmospheric correction aims to transform the TOA radiance into the top-of-canopy reflectance ${\rho }_{TOC}$. Similar to the study of Eugenio et al. [38], the atmospheric correction can also be performed by using the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer code [44,45,46]. Specifically, the ${\rho }_{TOC}$ can be estimated by using the following equations:

$$\{\begin{array}{c}y=x_a·L-x_b\\ {\rho }_{TOC}=\frac{y}{1+x_c·y}\end{array},$$

(6)

where $x_a$, $x_b$, and $x_c$ are the atmospheric correction coefficients calculated by the 6S code.

However, the aerosol optical thickness is also needed to run the 6S code, in addition to the spectral response functions of the sensor bands, the view zenith and azimuth angles, the sun zenith and azimuth angles, the date and time of the image acquisition, and the latitude and longitude of the center of the image scene.

In accordance with Tian et al. [47] and Zheng et al. [48], we first select the minimum TOA radiance of the adjacent optically deep waters to be the dark-object in this study. Then, we change the aerosol optical thickness from 0.0 to 2.0 in intervals of 0.001 and run the 6S code for each of the aerosol optical thicknesses. Consequently, the radiance of the dark-object is transformed into different values of ${\rho }_{TOC}$ corresponding to different aerosol optical thicknesses by using Equation (6), which allows us to select the aerosol optical thicknesses that cause the ${\rho }_{TOC}$ of the near-infrared bands (Near-IR1 and Near-IR2 for WV-3 and Near-Infrared for ZY-3) to become as close as possible to zero (which is the optimal aerosol optical thickness).

Thereafter, the optimal aerosol optical thickness is then input into the 6S code to obtain the atmospheric correction coefficients. Therefore, the ${\rho }_{TOC}$ image can finally be calculated by using Equation (6). The atmospheric correction approach is briefly presented in Figure 7.

Sunglint Correction and Subsurface Remote Sensing Reflectance Estimation

Because we actually use the remote sensing reflectance $R_{rs}$ rather than the TOC reflectance ${\rho }_{TOC}$ in further processing, ${\rho }_{TOC}$ must be further transformed into the remote sensing reflectance. As the 6S code can also directly provide the values of solar irradiance $E_{dir}$ and atmospheric diffuse irradiance $E_{dif}$, we can approximate the remote sensing reflectance $R_{rs}$ using the following equation [49,50,51,52,53]:

$$R_{rs}=\frac{{\rho }_{TOC}}{\pi }-\frac{r_{sky}·E_{dif}}{\pi (E_{dir}+E_{dif})},$$

(7)

where $r_{sky}$ is the reflectivity of the gas-water interface that can be calculated using the Fresnel reflection formula as follows [49,50,51]:

$$r_{sky}=\frac{1}{2}\left\{{\left[\frac{\sin ({\theta }_i-{\theta }_v)}{\sin ({\theta }_i+{\theta }_v)}\right]}^2+{\left[\frac{\tan ({\theta }_i-{\theta }_v)}{\tan ({\theta }_i+{\theta }_v)}\right]}^2\right\},$$

(8)

where ${\theta }_i$ and ${\theta }_v$ are the view zenith angle and the corresponding refraction angle, respectively.

On the other hand, sunglint is usually caused by the specular reflection of solar radiation on the slopes of the waves of surface winds. Thus, it is difficult to avoid and may seriously influence the bathymetry results, especially for images with high spatial resolution. Fortunately, a sunglint correction algorithm has already been developed by Hochberg et al. [54], and it was improved and used to correct the sunglint of high-spatial-resolution images for the estimation of water depth and diffuse attenuation coefficients by Eugenio et al. [38,55]. Thus, to perform the sunglint correction in this study, we simply applied Eugenio’s method to the remote sensing reflectance image. For convenience, we defined the sunglint-corrected remote sensing reflectance corresponding to the remote sensing reflectance $R_{rs}$ as $R_{rs}^{\prime }$.

The bathymetric Model (3) provides only the subsurface remote-sensing reflectance $\rho$, rather than the above-surface remote sensing reflectance $R_{rs}^{\prime }$. Therefore, we further defined the subsurface remote sensing reflectance obtained from $R_{rs}^{\prime }$, according to the relationship between the above-surface and subsurface remote sensing reflectances provided by Lee et al. [29,30], as follows:

$$\rho =\frac{R_{rs}^{\prime }}{0.518+1.562R_{rs}^{\prime }}.$$

(9)

3.3. Bathymetric Model and Its Solution

3.3.1. Reduce the Number of Unknowns

Many current multispectral satellite images can provide only 3–4 visible wavebands; for example, the ZY-3 imagery in this study provides only three visible wavebands (the blue, green, and red bands). As a result, we must reduce as many unknown variables of the bathymetric model as possible. Otherwise, the bathymetric model may become ill-conditioned, and the approach will end in failure.

According to the work of Prieur and Sathyendranath [56] and Lee et al. [29], the phytoplankton pigment absorption coefficient at 440 nm, $a_{\varphi }(440)$, can be further parameterized using the chlorophyll-a concentration $C$ as follows:

$$a_{\varphi }(440)=0.06C^{0.65}.$$

(10)

As illustrated in Figure 3, when Model (3) is used, the backscattering coefficients $b_{bp}(\lambda )$ of the particles can be parameterized, using $b_{bp}(400)$ and $Y$, as follows:

$$b_{bp}(\lambda )=b_{bp}(400){\left(\frac{400}{\lambda }\right)}^Y.$$

(11)

To reduce the number of unknown variables, $b_{bp}(\lambda )$ is further updated using the chlorophyll-a concentration $C$ (as shown in Appendix A and Appendix B) as follows:

$$b_{bp}(\lambda )=b_{bp}(550){\left(\frac{550}{\lambda }\right)}^Y,$$

(12)

where $Y=0.67875$, and

$$b_{bp}(550)=0.0111C^{0.62}.$$

(13)

To summarize, the model described in Figure 3 can now be updated using Equations (10), (12), and (13), and a constant of $Y=0.67875$. As a result, the number of unknown variables is reduced to only four: $H$, $B$, $C$, and $a_g(440)$.

To further reduce the number of unknowns, we note that

$${\rho }_{\infty }(\lambda )=g_{0}u(\lambda )+g_1{u}^2(\lambda )$$

(14)

can be rewritten as

$$u(\lambda )=\frac{-g_0+\sqrt{g_0^2+4g_1{\rho }_{\infty }(\lambda )}}{2g_1}.$$

(15)

Therefore, we can further separate the updated semi-analytic bathymetric model into two parts, as shown in Figure 8. In other words, we can solve the bathymetric problem with the following two steps: (a) according to Part 2 of the updated semi-analytic bathymetric model, determine the optimal values of $a_g(440)$ and $C$, so that the differences between the values of the estimated ${\rho }_{\infty }(\lambda )$ and the subsurface remote sensing reflectance of the optical deep waters are as small as possible, thus allowing the diffuse attenuation coefficient $K_d(\lambda )$ to be estimated; and (b) according to Part 1 of the updated semi-analytic bathymetric model, determine the optimal values of $H$ and $B$ for every pixel, so that the differences between the values of the estimated subsurface remote sensing reflectance $\rho$ and the observed values calculated by Equation (9) are as small as possible. We will illustrate these two steps in the following two sections, respectively.

3.3.2. Solve the Diffuse Attenuation Coefficient

Eugenio et al. [55] adapted the QAA for WV-2 multispectral images to estimate the values of phytoplankton attenuation $a_{\varphi }$, gelbstoff and detritus attenuation $a_g$, and particle backscattering $b_{bp}$. In this study, we present a similar approach to the adjacent optically deep waters according to Part 2 of the updated semi-analytic bathymetric model (Figure 8). In contrast to Eugenio et al. [55], our model is parameterized using $a_g(440)$ and $C$, which can be applied to multispectral images that have no fewer than 2 visible bands.

Specifically, to solve the parameters of $a_g(440)$ and $C$, we first change $a_g(440)$ from 0 m⁻¹ to 0.5 m⁻¹ using intervals of 0.001 m⁻¹, then change $C$ from 0.1 mg/m³ to 5.0 mg/m³ using intervals of 0.1 mg/m³, and calculate the corresponding subsurface remote sensing reflectance of each pair of $a_g(440)$ and $C$ using Part 2 of the updated semi-analytic bathymetric model (Figure 8). Then, the initialization of the values of $a_g(440)$ and $C$ are selected by making the sum of squared differences between the calculated subsurface remote sensing reflectance and the observed reflectance to be as close to zero as possible. Thereafter, the Levenberg–Marquardt method is used to solving for the optimal values of $a_g(440)$ and $C$. Finally, the average diffuse attenuation coefficient $K_d(\lambda )$ is then estimated from these optimal $a_g(440)$ and $C$ values, according to Part 2 of the updated semi-analytic bathymetric model (Figure 8). The estimated $K_d(\lambda )$ here is thus assumed to be known and is input into Part 1 of the model to complete the following steps of the proposed approach.

The effectiveness of the frame of the QAA approach has already been proved by Eugenio et al. [55], so we only need to further argue that the estimated diffuse attenuation coefficient here is approximately constant in the study area for each of the two experimental images.

To prove this, several deep-water regions were first designated manually from the two experimental images. Then, 100 pixels were randomly selected from the designated regions. The diffuse attenuation coefficients corresponding to these pixels were finally estimated by using our approach. The statistical analyses of these results are shown in

Table 2. Statistics of the Diffuse Attenuation Coefficients (WV-3).

Bands	Coastal	Blue	Green	Yellow	Red	RedEdge
Mean	0.0482	0.0529	0.0714	0.2776	0.4289	1.3665
Std. Deviation	0.00248	0.00201	0.00112	0.00091	0.00150	0.00073
Skewness	0.209	0.296	0.373	0.390	0.426	0.358
Kurtosis 1	0.006	0.669	0.069	−0.025	0.432	−0.265

¹ The Kurtosis has already been reduced by 3.

and

Table 3. Statistics of the Diffuse Attenuation Coefficients (ZY-3).

Bands	Blue	Green	Red
Mean	0.0371	0.0677	0.4143
Std. Deviation	0.00947	0.00329	0.00265
Skewness	−0.111	0.481	0.784
Kurtosis 1	0.446	−0.050	0.088

¹ The kurtosis has already been reduced by 3.

, as follows.

In Table 2 and Table 3, the standard deviation of the diffuse attenuation coefficient is comparatively small, as are the values of skewness and kurtosis. This means that the diffuse attenuation coefficients estimated from most deep-water pixels are close to their means. In other words, the diffuse attenuation coefficients obtained from both the WV-3 image and the ZY-3 image prove to be approximately constants.

3.3.3. Initialize the Water Depths

The bottom irradiance reflectance ${\rho }_b(\lambda )$ is usually expressed as follows [30,37,38]:

$${\rho }_b(\lambda )=B{\rho }_N(\lambda )$$

(16)

where ${\rho }_N(\lambda )$, which is usually assumed to be known in most bathymetric approaches, $B$ is the spectral shape of the reflectance normalized at the reference wavelength. For example, Lee et al. [30] and Eugenio et al. [38] used a 550-nm-normalized sand-albedo shape, and Hu et al. [37] tried to derive ${\rho }_N(\lambda )$ from the spectral reflectance values of sand, algae, dark sediment, and their mixtures.

To estimate water depths using Part 1 of the updated semi-analytic bathymetric model (Figure 8), we must also specify the shape of ${\rho }_N(\lambda )$ prior to running our bathymetric approach. Equation (3) illustrates that if the water depth is close to zero, the subsurface remote sensing $\rho$ will be approximately equal to

$${\rho }_0(\lambda )=\frac{{\rho }_b(\lambda )}{\pi }.$$

(17)

Therefore, we manually select several typical subsurface remote sensing reflectance shapes located close to the transient coastline as references. We assign these subsurface remote sensing reflectance shapes values of ${\rho }_0^{(i)}(\lambda )$ ($i=1,2,\cdots ,n_{\rho }$), where $n_{\rho }$ is the number of selected reference shapes. Considering that the water depths of the selected shapes are not actually equal to zero, we parameterized the bottom irradiance reflectance shapes according to Equation (17), as follows:

$${\rho }_b^{(i)}(\lambda )=B\pi {\rho }_0^{(i)}(\lambda )$$

(18)

or

$${\rho }_N^{(i)}(\lambda )=\pi {\rho }_0^{(i)}(\lambda ).$$

(19)

Thereafter, in principle, we can estimate pixel-by-pixel water depths according to Part 1 of our updated semi-analytic bathymetric model (Figure 8) using the Levenberg–Marquardt algorithm. However, the initialization of the water depth must be performed prior to applying the Levenberg–Marquardt algorithm. In this study, the initial values of $H$ and $B$ are estimated using the following steps:

(a) For each value of ${\rho }_N^{(i)}(\lambda )$ ($i=1,2,\cdots ,n_{\rho }$), change $H$ from 0 m to 40 m using intervals of 0.5 m, change $B$ from 0.5 to 1.5 using intervals of 0.01, and calculate the subsurface remote sensing reflectance ${\rho }_{HB}^{(i)}(\lambda )$ corresponding to every pair of $H$ and $B$ according to Part 1 of the model, as shown in Figure 8.

(b) For each value of $i$, build a K-D tree [57] based on the calculated subsurface remote sensing reflectance ${\rho }_{HB}^{(i)}(\lambda )$, and define this K-D tree as $KdTre{e}^{(i)}$, where K-D tree is known as a space-partitioning data structure for organizing points in a k-dimensional space, and it was useful for us to efficiently search the closest reflectance of the observed one from these calculated ones ${\rho }_{HB}^{(i)}(\lambda )$ in the following step.

(c) For each observed value of subsurface remote sensing reflectance $\rho (\lambda )$ (calculated by Equation (9)), determine its closest corresponding reflectance value from $KdTre{e}^{(i)}$, define the result as ${\rho }_{HB\_\min }^{(i)}(\lambda )$, and define the corresponding square distance as $D\left[\rho (\lambda ),{\rho }_{HB\_\min }^{(i)}(\lambda )\right]$.

(d) For each observed value of subsurface remote sensing reflectance $\rho (\lambda )$, determine the minimum distance from the set of $\left\{D_i|D_i=D\left[\rho (\lambda ),{\rho }_{HB\_\min }^{(i)}(\lambda )\right],i=1,2,\cdots ,n_{\rho }\right\}$, define the result as $D_{i\_\min }$, and define the corresponding normalized reflectance spectral shape as ${\rho }_N{}^{(i\_\min )}(\lambda )$.

As a result, the value of ${\rho }_N{}^{(i\_\min )}(\lambda )$ is selected as the normalized reflectance spectral shape corresponds to the observed $\rho (\lambda )$. Furthermore, based on such normalized reflectance spectral shape, the values of $H$ and $B$ corresponding to $D_{i\_\min }$ are selected to initialize the Levenberg–Marquardt algorithm. Finally, the result of the Levenberg–Marquardt algorithm is then used as the initialization to perform the following global adjustment.

3.4. Global Adjustment of Water Depths

In contrast to most bathymetric methods, we further propose a global adjustment method to refine the estimated water depths, for trying to smooth the water depths and improve the accuracies.

As mentioned in Section 3.1, the principle is based on the following two conditions: (a) most water depths should be close to each other for any two adjacent pixels; and (b) if the pixel size is small, most water depths of the pixels that are close to the transient coastline should be approximately zero. The specific computational model is described as follows.

First, we define the water depth $H$ and the parameter $B$ of a given pixel $i$ as $H_i$ and $B_i$, respectively, and define the corresponding normalized reflectance spectral shape (see Section 3.3.3) as ${\rho }_{N,i}(\lambda )$. As a result, Condition (a) can be described as follows:

$$H_i-H_j\to 0,$$

(20)

where $j$ represents one pixel that is adjacent to pixel $i$. In addition, if pixel $i$ is located close to the transient coastline, then Condition (b) can be described as follows:

$$H_i\to 0.$$

(21)

Therefore, by considering Part 1 of the updated semi-analytic model (Figure 8) together with Equations (20) and (21), we can construct the following objective function to refine the estimated water depths:

$$L(H_i,B_i,i\in I)={\displaystyle \sum _{i\in I}{({\rho }_{o,i}-{\rho }_i)}^T}({\rho }_{o,i}-{\rho }_i)+W_{\Delta }{\displaystyle \sum _{i\in I}{\displaystyle \sum _{j\in N_i}(H_i-H_j}}{)}^2+W_0{\displaystyle \sum _{i\in H_0}{H}_i^2}$$

(22)

where $I$ represents the region of the image that is used to perform the bathymetric approach, $N_i$ is the adjacent pixel set of $i$ in region $I$, and $H_0$ is the pixel set where the water depths are close to zero. In this study, $H_0$ consists of pixels that are located close to the transient coastline (at a distance of less than two pixels); ${\rho }_{o,i}$ is the observed subsurface remote sensing reflectance vector calculated by Equation (9); and ${\rho }_i$ is the subsurface remote sensing reflectance vector expressed by Equation (3), which is dependent on the unknown parameters $H_i$ and $B_i$. $W_{\Delta }$ and $W_0$ are the weights of the corresponding items of the objective function. $W_{\Delta }$ can usually be determined according to the observed local variations in the water depths; if the actual variation in the water depth is large, then $W_{\Delta }$ should be set to a relatively small positive number. To avoid the over-smoothing of the water depths, we set $W_{\Delta }$ to a small number in these experiments, i.e., $W_{\Delta }=0.001$. $W_0$ is used to constrain the water depths of the pixels around the transient coastline to be as close to zero as possible; thus, in this study, $W_0=2.0$.

The first item of the objective function causes the results to fit the radiative transfer Model (3) as closely as possible. The second item is a smoothing item, which can reduce the influence of noise on the final water depth results. When it is difficult to determine the pixels with water depths of approximately zero, the third item of the objective function can be ignored in the global adjustment approach.

In principle, all values of the water depth $H_i$ and the parameter $B_i$ can be globally adjusted by minimizing the proposed objective function. Specifically, we started at the values of water depth $H_i$ and parameter $B_i$ that were estimated in Section 3.3.3 and employed the trust–region–reflective method proposed by Coleman and Li [58,59] to optimize the objective function in this study.

Note that in case of some pixels that are close to the transient coastline but are not really zeroes in depth, we should discard such pixels from $H_0$ before performing the global adjustment.

4. Results

We have already applied the proposed bathymetric approach to the two experimental images that were shown in Section 3.2.1. Accordingly, the results will be illustrated and discussed in the following sections.

4.1. Visualization of the Water Depths

The estimated water depths before and after the global adjustment are shown in Figure 9 (ZY-3) and Figure 10 (WV-3). To facilitate these comparisons, the water depths shown in Figure 10 have already been corrected to refer to the tide level of the ZY-3 image.

4.2. Evaluations of the Accuracies before and after the Global Adjustment

Since water depths at 45 points around the Weizhou Island have been measured during October 2015, the accuracy of the proposed bathymetric approach can be then evaluated by comparing the estimated water depths with the measured ones. These comparisons are shown in Figure 11, and the corresponding statistics are illustrated in

Table 4. Statistics of the differences between the estimated and measured water depths ¹ (Unit: m)

Satellite	Before Global Adjustment			After Global Adjustment
	MAX	MAE	RMSE	MAX	MAE	RMSE
ZY-3	2.587	0.97	1.22	2.265	0.81	1.01
WV-3	2.11	0.58	0.75	2.27	0.56	0.77

¹ Note: MAX represents the maximum value of the absolute errors between the estimated and measured water depths, MAE represents the mean value of the absolute errors, and RMSE represents the root mean square error, which is calculated using the following equation:

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{(H_i^m-H_i^e)}^2}}{n}}$$

(23)

where $H_i^m$ is the measured water depth and $H_i^e$ is the water depth estimated using the proposed approach.

.

It should be note that although the gap between the time that the ZY-3 image was obtained and the water depth measurements were acquired was approximately 40 days, the water bottom around Weizhou Island changed little during this time due to the following three reasons: (a) no very strong tropical storms passed through Weizhou Island during this time; (b) the island is piled up under the water when Quaternary basalt magma gushed out, and now has a relatively stable structure [60]; and (c) the distributions of the water depths estimated from the two images were seen to be similar to each other, as is seen in the visualization of Section 4.1.

4.3. Comparisons between the Water Depths Estimated from the Two Different Images

In contrast to other works, we compared the differences between the water depths estimated from the ZY-3 and the WV-3 images for further testing the performances of the proposed approach. Specifically, we randomly selected 200 points from the two images and quantitatively compared the results obtained from these two images. The results are shown in Figure 12. When the statistics are further estimated, we find that the maximum value of the absolute errors (MAX) is 4.69 m, the mean of the absolute errors (MAE) is 1.02 m, and the root mean square error (RMSE) is 1.24 m, where the errors are calculated by using the differences between the water depths estimated from the ZY-3 and the WV-3 images.

5. Discussion

5.1. Discussion on the Visualization of the Water Depths

We can see from Figure 9 and Figure 10 that the water depths are slightly smoother after the global adjustment, especially for the ZY-3 image. To more carefully compare the water depths obtained before and after the global adjustment, we further examined an arbitrary transect of water depths shown in Figure 13.

When comparing only the water depths obtained before and after the global adjustment for the ZY-3 image and the WV-3 image, the water depths along the transect are smoother after performing the global adjustment, as shown in Figure 13.

Note that the ZY-3 image has a lower SNR, and this may be the reason why the smoothing effect of the ZY-3 image was more significant than the WV-3 image.

On the other hand, the following aspects are also illustrated in Figure 13: (a) when the water depths are less than 5 m, the water depths estimated from the WV-3 image are slightly greater than those estimated from the ZY-3 image, although most of these differences are less than approximately 1 m; (b) when the water depths range from 5 m to 10 m, the water depths obtained from the two different images agree well with each other; and (c) when the water depths are more than 10 m, the differences between the two images may reach up to 2 m (see the dotted circle in Figure 13).

When the water depths are lower than 5 m (i.e., about 2 m without tide corrections), the proposed approach may be more susceptible to the errors of the spectral shape of the selected bottom reflectance. This is regarded as the reason why the differences between the water depths estimated from the two images become slighter and slighter as the water depths increase until about 10 m, as shown in Figure 13.

As can also be seen in Figure 13, before the global adjustment, the local variations in water depth increase with increasing water depth. This is due to the decreases in the ratio between the radiations reflected by the water bottom and those coming from the water column. Thus, the influences of the noise become increasingly significant, so that the local variations of water depths become large, as do the differences between the water depths estimated from the two images as seen in the dotted circle of Figure 13.

Besides, we note that the ZY-3 and the WV-3 images are acquired at 11:23 a.m., 24 August 2015, with the sun zenith angle of 22.1°, and 11:37 a.m., 1 October 2015, with the sun zenith angle of 27.6°, respectively. Hence, the WV-3 image may correspond to a weaker direct solar radiation intensity than the ZY-3 image. This can also cause the local variations of the water depths estimated from the WV-3 image being greater than the ZY-3 image at the places where the water depths are deeper than 10 m, as shown in the dotted circle of Figure 13.

5.2. Discussion on the Accuracy Assessment of the Estimated Water Depths

Except for the visualization results, we focus more on the quantitative statistics of the differences between the estimated and measured water depths.

As shown in Figure 11, the water depths estimated for the ZY-3 image appear to be closer to the measured depths after undergoing the global adjustment sub-approach, but this trend is not remarkable for the WV-3 image. This phenomenon is further confirmed by the statistics shown in Table 4: the MAE and RMSE changed from 0.97 m and 1.22 m before the global adjustment to 0.81 m and 1.01 m after undergoing the global adjustment for the ZY-3 image, while from 0.58 m and 0.75 m to 0.56 m and 0.77 m for the WV-3 image.

When Figure 13 was further reviewed, we can clearly find that the proposed global adjustment sub-approach is able to adjust some possible large outliers to the relatively reasonable values, as shown in the dotted boxes of Figure 13. At the same time, we also find that if the variation of the adjacent water depths is not remarkable, the global adjustment sub-approach can hardly influence the final water depth estimation. Accordingly, we speculate that the reasons for the mentioned phenomenon may include the following aspects:

(a) The SNR of the ZY-3 image is lower than that of the WV-3 image; therefore, the influences of noise were more serious for ZY-3 than for WV-3.

(b) The band number of the ZY-3 image is considerably less than that of the WV-3 image; thus, the bathymetry of the ZY-3 image is much more sensitive to noise than that of the WV-3.

(c) The smoothing item of the objective function (see Equation (22)) works well when the noises are sensitive or the influences of the noises are serious.

These factors indicate that the global adjustment sub-approach can not only improve the accuracy when the influences of the noises are remarkable, but can also avoid the over-smoothing of the water depths when the results are not sensitive to the noises. This is what common smoothing filters are difficult to realize.

On the other hand, if we focus only on the statistics after the global adjustment as shown in Table 4, we note that the proposed approach have exceeded those of Hu et al. [37] and Eugenio et al. [38] on RMSEs, as shown in

Table 5. Comparisons of the RMSEs obtained using different methods.

Author/Method	Location	Data	RMSE (Unit: m)
Hu et al. (2014)	Pratas Atoll	Landsat ETM+	1.6
	North Danger Reefs	Landsat ETM+	1.1
Eugenio et al. (2015)	Granadilla area of Tenerife Island	WV-2	1.20
	Corralejo area of Fuerteventura Island	WV-2	1.94
Our Approach	Weizhou Island	ZY-3	1.01
	Weizhou Island	WV-3	0.77

. Although these RMSEs are calculated by using different images with different resolutions, we are more inclined to believe that they are comparable to each other. The reason is that we can see from Table 5 that Hu et al. [37] obtained better RMSEs (1.1 m and 1.6 m) with 30 m resolution of Landsat ETM+ images than those of Eugenio et al. [38] (1.20 m and 1.94 m) with 2.0 m resolution of WV-2 images, while the bathymetric model used by Hu et al. [37] (see Equation (2)) can be seen as a simplification or an approximation of the one used by Eugenio et al. [38] (see Equation (3)). Actually, we used a ZY-3 multispectral image with a resolution of 5.8 m and a WV-3 multispectral image with a resolution of 1.24 m in the experiments, and obtained the RMSEs of 1.01 m and 0.77 m respectively, which were less than those of Hu et al. [37] and Eugenio et al. [38].

When we further compare the proposed approach with those of Hu et al. [37] and Eugenio et al. [38], we can see that the proposed approach is not only better than other methods on accuracies as mentioned previously, but also has the following advantages:

(a) This approach can extract water depths directly from the multispectral image itself and does not require the use of extra MODIS data, water depth data, or water column data, except for the evaluation of its results;

(b) This approach has fewer variables than many other current methods, which means that it can be used to analyze most current multispectral satellite images, even some images that contain only 2–3 visible wavebands.

Note that the proposed approach reduced the unknown variables through parameterizing the radiative transfer model using only two unknown parameters (see Section 3.3.1) instead of six of Eugenio et al. [38] (see Equation (3) in Section 2). This seems to increase the model errors for the proposed approach, and may lead to the statistics of the differences between the estimated and measured water depths becoming worse, but this is not true as shown in Table 4 and Table 5. Hence, we believe that the parameterization is reasonable without causing significant model errors in this study.

In addition, Figure 12 further demonstrates that these comparisons between the water depths estimated from the ZY-3 and the WV-3 images are similar to the results derived from the transect shown in Figure 13; the results of statistical analysis indicate that the value of the root mean square error is 1.24 m (see Section 4.3).

In fact, we can also estimate the RMSE of the differences between the water depths estimated from the ZY-3 and WV-3 images according to the results shown in Table 4 as follows:

$$RMS{E}_{diff}=\sqrt{RMS{E}_{ZY-3}^2+RMS{E}_{WV-3}^2}=\sqrt{{1.01}^2+{0.77}^2}=1.27(m).$$

(24)

Such value is relatively small and very close to 1.24 m. This indicates that the water depths estimated from the ZY-3 and WV-3 images are consistent to each other, and then is able to further enhance our confidence in the effectiveness and the accuracies of the proposed approach.

5.3. Discussion on Supplementary Notes and Future Work

We note that the global adjustment sub-approach depends slightly on the assumption that most water depths should be similar to each other between one pixel and its adjacent pixels. This may indicate that this method involves an implicit assumption that the detected sea floor elevations are continuous and vary gradually. Thus, when the sea floor is not sufficiently continuous, this refinement may become inaccurate or may even fail. Fortunately, when the weight of the smooth item (the second item of the objective function, Equation (22)) is set as a small value ($W_{\Delta }=0.001$), we find that the smooth item mainly works when the noises are sensitive or their influences are serious, as discussed in Section 5.2. In other words, when the noises of the images are significantly sensitive or the influences of the noises are serious, the global adjustment approach works well on improving the water depth results, otherwise the global adjustment approach is able to avoid the over-smoothing of the water depths.

As a result, we believe that the proposed global adjustment sub-approach can be applied to estimate the water depth of many types of sea floors. Specifically, the proposed approach has great potential for the marine surveying and mapping of large, remote, or dangerous areas in the future, particularly for further research on the bathymetries of typical coral reefs.

Thus, more different images and regions are needed to further test and improve this approach in our next works. Moreover, we also plan to improve this approach by further considering the variations in the properties of the water column between shallow and deep waters in future studies that may introduce the global adjustment sub-approach to smooth the water-column properties.

6. Conclusions

In this paper, a QAA method is firstly explored to estimate the diffuse attenuation coefficients, and then the parameterization of bathymetric model is further simplified. Based on these two aspects, we have developed a QAA approach that includes sub-approaches of initialization and a novel global adjustment to estimate the water depths around Weizhou Island from multispectral images. One advantage is that the approach is able to reduce the number of unknowns and is independent of auxiliary datasets. Another advantage is that the global adjustment sub-approach can not only smooth the final water depths and improve the accuracies, but can also avoid the over-smoothing of water depths. Experimental results have finally illustrated that the water depths estimated from the ZY-3 and WV-3 images are approximately consistent to each other, and yield EMSEs of 1.01 m and 0.77 m, respectively. Therefore, the proposed approach is thought to be competitive to other QAA bathymetric methods.