A Semi-Empirical Chlorophyll-a Retrieval Algorithm Considering the Effects of Sun Glint, Bottom Reflectance, and Non-Algal Particles in the Optically Shallow Water Zones of Sanya Bay Using SPOT6 Data

Abstract

Chlorophyll-a (Chl-a) concentration retrieval is essential for water quality monitoring, aquaculture, and guiding coastline infrastructure construction. Compared with common ocean color satellites, land observation satellites have the advantage of a higher resolution and more data sources for retrieving the concentration of Chl-a from optically shallow waters. However, the sun glint (R_sg), bottom reflectance (R_b), and non-algal particle (NAP) derived from terrigenous matter affect the accuracy of Chl-a concentration retrieval using land observation satellite image data. In this paper, we propose a semi-empirical algorithm based on the remote sensing reflectance (R_rs) of SPOT6 to retrieve the Chl-a concentration in Sanya Bay (SYB), considering the effect of R_sg, R_b, and NAP. In this semi-empirical algorithm, the Cox–Munk anisotropic model and radiative transfer model (RTM) were used to reduce the effects of R_sg and R_b on R_rs, and the Chl-a concentration was retrieved by the Chl-a absorption coefficient at 490 nm (a_phy(490)) to remove the effect of NAP. The semi-empirical algorithm was in the form of Chl-a = 43.3[a_phy(490)]^1.454, where a_phy (490) was calculated by the total absorption coefficient and the absorption coefficients of each component by empirical algorithms. The results of the Chl-a concentration retrieval show the following: (1) SPOT6 data are available for Chl-a retrieval using this semi-empirical algorithm in oligotrophic or mesotrophic coastal waters, and the accuracy of the algorithm can be improved by removing the effects of R_sg, R_b, and NAP (R² from 0.71 to 0.93 and root mean square error (RMSE) from 0.23 to 0.11 ug/L); (2) empirical algorithms based on the blue-green band are suitable for oligotrophic or mesotrophic coastal waters, and the algorithm based on the blue-green band difference Chl-a index (DCI) has stronger anti-interference in terms of the effects of sun glint and bottom reflectance than the algorithm based on the blue-green ratio (BGr); (3) in the case of ignoring R_sg unrelated to inherent optical properties (IOPs), NAP is the biggest interference factor when >9.5 mg/L and the effect of bottom reflectance should be considered when the water depth (H) <5 m in SYB; and (4) the inherent optical properties of the waters in SYB are dominated by NAP (Chl-a = 0.2–2.6 ug/L and NAP = 2.2–30.1 mg/L), and the nutrients are concentrated by enclosed terrain and southeast current. This semi-empirical algorithm for Chl-a concentration retrieval has the potential to monitor Chl-a in oligotrophic and mesotrophic coastal waters using other land observation satellites (e.g., Landsat8 OLI, ASTER, and GaoFen2).

1. Introduction

Though coastal zones only cover 2% of the world’s land area, almost two-thirds of urban settlements with a population higher than 5 million live in these areas [1]. In particular, the optically shallow regions of coastal waters have the most diverse marine life and are significant for aquaculture and tourism. Monitoring the concentration of chlorophyll-a (Chl-a) in coastal waters by using ocean color satellites (e.g., SeaWiFs, MODIS, SeaWiFS, and MERIS) is one of the critical methods for water quality monitoring [2,3]. However, the retrieval of Chl-a concentration from optically shallow waters requires satellites with a higher spatial resolution. Therefore, the multi-spectral land observation satellites (e.g., SPOT6, Landsat8 OLI, ASTER, and GaoFen2) may have the a potential to retrieve Chl-a concentration from optically shallow waters.

According to the Chl-a absorption characteristics and band sets of satellites, selection of the best bands and different algorithms is a continuously hot research topic. For Case-1 waters without terrigenous matter inputs (i.e., sediments and colored dissolved organic matter (CDOM)) or resuspended sediments, the optical properties of the ocean waters are tightly subordinated to the abundance of pigmented algal cells [4,5]. According to the optically characteristics of oceanic waters, the bio-optical algorithms based on the blue-green band-ratio (BGr) are suitable for oligotrophic waters [6,7,8]. According to the absorption peak of Chl-a near 440 nm, the blue band at 443 nm was adopted in the early BGr algorithm [9]. To minimize the interference of the strong absorption of CDOM in the blue band, using in the algorithms based on BGr (e.g., SeaWiFS OC3, OC4, and MODIS OC3M), the blue band was red-shifted from 440 nm to 490 nm or 510 nm [10,11,12]. Though the Chl-a retrieval algorithm for Case-1 waters is relatively mature and there are several algorithms working showing how to estimate from Case-2 waters, the band selection for Chl-a retrieval is still a difficult task in waters with complex constituents [13,14].

With the inflow of terrigenous matter from the runoff and the resuspension of bottom sediment, non-algal particle (NAP) is the main interfering factor for Chl-a retrieval in optically complex waters (i.e., Case 2 waters). The quasi-analytical algorithm (QAA) [15] based on the visible bands and the empirical or semi-empirical algorithm based on the near-infrared (NIR) bands (e.g., the three-band algorithm or the four-band algorithm) are effective in reducing the effects of CDOM and NAP interference in eutrophic waters [16,17,18]. However, for most land observation satellites (e.g., SPOT6, Landsat8 OLI, ASTER, and GaoFen2), the spectral resolution and the lack of NIR band near 705 nm do not meet the requirements of the QAA algorithm and algorithms based on the NIR band. However, for mesotrophic or oligotrophic waters, the low reflectance may cause a reduction in the signal-to-noise ratio (SNR) and, based on the BGr, both MERIS OC3E and OC4E perform better than algorithms based on NIR bands [14]. However, the spectral resolution and the lack of NIR band near 705 nm of land observation satellites (e.g., SPOT6, Landsat8 OLI, ASTER, and GaoFen2) do not meet the requirements of the QAA algorithm. The empirical or semi-empirical algorithms based on the near-infrared (NIR) bands are effective in reducing the effects of CDOM and NAP interference in eutrophic waters [19,20,21]. For land observation satellites, the lack of NIR band near 705 nm limits the ability of Chl-a retrieval in Case-2 waters. However, for mesotrophic or oligotrophic waters, the low reflectance may cause a reduction in the SNR, and algorithms based on the BGr, both MERIS OC3E and OC4E, perform better than algorithms based on NIR bands [14]. The two-rand-ratio (2Br) algorithm using Landsat8 OLI performs well in Case-2 waters [22,23,24]. Therefore, algorithms based on visible bands have the potential to retrieve the concentration of Chl-a in coastal waters.

For optically shallow waters, the acquisition of remote sensing reflectance or water-leaving reflectance (R_rs) is limited by the effects of boundary conditions, that is, sun glint and bottom reflectance. The contribution of sun glint to R_rs (R_sg) mainly varies with the sun–satellite relative geometric relationship [25,26]. In the images of common ocean color sensors (e.g., MODIS, SeaWiFS, and MERIS), the trace of R_sg can be easily found because medium spatial resolution satellites have a wide field of view [27,28]. The view angle of sensors for coastal zone observation with high spatial resolution changes little, and sun glint may be ignored when the R_sg is lower than the water-leaving reflectance [29]. On the same remote sensing image, the R_sg changes with the wavelength owing to the variations of the water refractive index, and the difference of R_sg in the range of 350–2500 nm is approximately 10% [30,31]. The R_sg may cause a lot of errors, especially in clear waters with low reflectance. The method of removing sun glint by the near-infrared (NIR) band does well in oceanic waters with a low concentration of NAP, but it is not appropriate for coastal waters with a high concentration of NAP. The Cox–Munk analytical anisotropic model has previously been used to calculate the value of sun glint, but its applicability in areas with low R_sg still needs to be validated.

The water-leaving reflectance contributed by the bottom (R_bw) depends on the bottom reflectance and the water column attenuation coefficient. The type of bottom is the main factor for bottom reflectance, because the slope of the seafloor changes little in coastal zones. Bottom types mainly include sand, rock, mud, coral, and seagrass, which all result in an overestimation of the concentration of Chl-a because the reflectance of the green band is higher than that of the blue band. The radiative transfer model (RTM) has previously been used to calculate the R_bw [32,33,34,35], but the attenuation of the water column and the bottom reflectance need to be used as the input parameters.

Compared with ocean color satellites with moderate spatial resolution (e.g., MODIS, SeaWiFS, and MERIS), multi-spectral land observation satellites (e.g., SPOT6, Landsat8 OLI, ASTER, and GaoFen2) have the advantage of a higher resolution in reducing the effect of bottom. We chose SPOT6 data as an example. Although SPOT6 was not designed to be an ocean color sensor, SPOT6 image data have been widely used in land vegetation monitoring [36,37,38] and have the potential to be applied to Chl-a concentration retrieval in coastal waters. According to the band sets of SPOT6 (blue: 455–525 nm; red: 530–590 nm) and a_phy (shown in Figure A1), the blue and red bands reflect the absorption and reflection characteristics of Chl-a, respectively. Therefore, the blue-green band of SPOT6 data may have the potential to retrieve the concentration of Chl-a in oligotrophic or mesotrophic waters. However, it has a shortcoming in that broad-band spectral data are more sensitive to suspension than hyper-spectral data [39,40,41].

For Chl-a retrieval using broad-band spectral data in optically shallow waters of coastal waters, the main problems are the effects of sun-glint, bottom reflectance, and NAP. To solve these problems, we proposed a semi-empirical Chl-a retrieval algorithm combining the Cox–Munk anisotropic model [25,26], the RTM [32,33,34,35], and a bio-optical model. This semi-empirical algorithm has the potential to retrieve Chl-a using other broad-band spectral data without NIR band in mesotrophic and oligotrophic waters.

2. Materials

2.1. Study Area

Sanya Bay (SYB) is in the southeast of Hainan Island, the second-largest island in China (Figure 1), and lies in a tropical monsoon climate [42]. The coastline of Hainan Island mainly includes sandy, rocky, biological, and artificial coastlines. The coastline of the study area is a sandy coastline, and nearly no pollution has been found in this area. Most of the coastal waters of SYB are open water areas with large waves, and the mean monthly wave height is 0.87–1.87 m [43]. The water depth range is 0–35 m and there are many kinds of benthos and fishes in this study area.

The waters in the study area mainly include three types: ultraoligotrophic, oligotrophic, and mesotrophic waters. The water types vary with the ocean currents and the distance to the shore. The sea waves pick up large amounts of terrigenous particles from shallow waters and send them to the deep-water regions. The turbidity is very high within 2 m of the water depth, and the waters are as clear as oceanic waters when the depth is greater than 30 m. The West Island mitigates the impact of the northeast current, which allows the southwest side of the island to retain a large amount of nutrients and sediment. On the southwest side of the West Island, the clear water and the wide shallow water region provide the best conditions for coral growth. Overall, the coastal waters of SYB are relatively clear and the optical properties are dominated by terrigenous suspended mineral particles.

2.2. SPOT6 Data

SPOT6 data have three visible and near-infrared bands (wavelengths: 450–520 nm; 530–590 nm; 625–695 nm; and 760–890 nm), and the spatial resolution is 6 m. In the application of ocean color remote sensing, the visible bands, and the NIR bands are used for monitoring phytoplankton and non-algal particulars, respectively. Compared with other common ocean color satellites, the higher spatial resolution can help to analyze the changes of Chl-a in coastal zones. According to the band sets of SPOT6 (blue: 455–525 nm; red: 530–590 nm) and a_phy (shown in Appendix A.1, Figure A1), the blue and red bands reflect the absorption and reflection characteristics of Chl-a, respectively. The image of SPOT6 was acquired on 21 January 2013.

2.3. In Situ Measured Dataset

The in situ measured parameters include wind speed, water depth, and water column parameters (i.e., Chl-a, NAP, and CDOM) (as shown in

Table 1. The measured water column parameters for the model input include the absorption coefficients, the water constituent concentrations, and the water depths. NAP, non-algal particle.

	Wavelength (nm)	aphy (m−1)	aNAP (m−1)	aCDOM (m−1)	aw (m−1)	bb*NAP (m2 g−1)	Chl-a (ug/L)	NAP (mg/L)	H (m)
mean	490	0.034	0.1	0.037	0.019	0.0156	0.587	5.565	7.803
	560	0.011	0.076	0.013	0.071	0.0150
	660	0.017	0.033	0.003	0.410	0.0141
max	490	0.171	0.562	0.059	-		1.952	23.570	15.000
	560	0.039	0.218	0.021	-
	660	0.063	0.151	0.005	-
min	490	0.023	0.02	0	-		0.205	2.670	1.0
	560	0.006	0.012	0	-
	660	0.011	0.006	0	-

). The in-situ samples were collected on 23 January 2013. The wind speed data were obtained from the meteorological station of the China Meteorological Administration, and the acquisition time was at 11:00 on 1 January 2013. The water depth-measured data used the single-point sonar data measured by the Guangzhou Marine Geological Survey in November 2014. The Chl-a concentration was obtained through the extraction of a filtered sample (0.45 μm What-man GF/F) with 90% acetone solution and the determination of the a_phy with a Shimadzu UV-2600PC spectrophotometer. The CDOM absorption spectra were determined using a Shimadzu UV-2600PC spectrophotometer, fitted with a 1 cm quartz cuvette, in the spectral region between 200 and 800 nm at 1 nm intervals. The water samples for the CDOM measurement were filtered through a 0.7 μm pre-combusted Whatman GF/F filter, and then through a 0.22 μm Millipore membrane cellulose filter. The absorption coefficients of the NAP were measured by a UV2600 PC spectrophotometer based on the filter pad methodology, and the membrane was 0.7 μm Whatman GF/F [23,44]. The measured absorption coefficients of Chl-a, CDOM, and NAP are shown in Figure 2.

The specific backscattering coefficient (b_b*) as input parameter for RTM can be calculated by the concentration of Chl-a and NAP, as well as their specific scattering coefficient (b*_phy and b*_NAP) and backscattering ratio (B_phy and B_NAP). The mean value of specific backscattering coefficient (b*) at 532 nm measured in mineral-dominated waters was 0.519 m² g⁻¹ [45] and the mean value of b*_NAP after removing b*_phy was 0.4 m² g⁻¹ [46]. The B_NAP at 488 nm ranged from 0.006 to 0.03 in Hudson Bay (Canada) with a relatively closed terrain and that measured in the English Channel with open terrain was 0.04 [46]. The b_b*_phy ranges from 0.04 for macrophytoplankton to 0.005 for microphytoplankton [47]. For the total attenuation coefficient in coastal waters with Chl-a < 2 ug/L, b_b*_phy was a negligible quantity. The specific backscattering coefficient for NAP (b_b* _NAP) was calculated by B_NAP and b* _NAP from Neil’s datasets measured by HOBI Labs HS–2 backscattering meter (B_NAP = 0.04 and b* _NAP = 0.4 m² g⁻¹ at 488 nm) [46].

3. Remote Sensing Reflectance Algorithm

In optically deep waters, the remote sensing reflectance or leaving water reflectance (R_rs) is defined by the ratio of water-leaving radiance (L_w) to downward irradiance (E_d) [48]. However, R_rs in optically shallow waters is affected by bottom. The available Chl-a concentration retrieval algorithm is based on the acquirement of the accurate R_rs (i.e., the remote sensing reflectance without the effect of the bottom). Although ocean color remote sensing has been studied for over half a century, the acquisition of accurate R_rs in optically shallow waters is still a challenge for Chl-a concentration retrieval [2,9,33]. After accurate atmospheric correction, the reflectance gathered by a sensor just above the water surface is referred to as apparent remote sensing reflectance (R_ap). In oceanic waters, R_ap is the sum of the sun glint (i.e., the specular reflectance (R_sg)) and the reflectance contributed by the water column (R_cw), and the reflectance in the sun glint areas may be far larger than R_cw. Although there are a number of novel models for sun glint calculation, the most common approach of the Cox–Munk model is still widely used [25,26,27,28]. In optically shallow waters, bottom reflectance is another boundary condition for the retrieval of R_rs, and the reflectance contributed by the bottom (R_bw) can be calculated using the RTM [32,33,34,35]. Assuming atmospheric correction is available, in order to retrieve R_rs or R_cw, knowledge of the propagation of light within a water body is necessary, and the boundary conditions and inherent optical properties (IOPs) are the main research objects (as shown in Figure 3).

The satellite-measured top-of-atmosphere (TOA) reflectance can be expressed as follows:

$${\rho }_{\mathrm{TOA}}\left(\mathsf{\lambda }\right)=t_{\mathrm{gas}}\left(\mathsf{\lambda }\right)\left(t\right(\mathsf{\lambda })·{\rho }_{\mathrm{t}}(\mathsf{\lambda })+{\rho }_{\mathrm{r}}(\mathsf{\lambda })+{\rho }_{\mathrm{A}}(\mathsf{\lambda }\left)\right)$$

(1)

where t_gas(λ), t(λ), ρ_t(λ), ρ_r(λ), and ρ_A(λ) are the gaseous transmittance, the total diffuse atmospheric transmission, the irradiance reflectance just above the water surface, the molecular Rayleigh scattering reflectance, and the aerosol reflectance, respectively. The term ρ_t(λ) can be calculated by the fast line of sight atmospheric analysis of spectral hypercube (FLAASH) module to eliminate the effect of the atmosphere. According to R_rs(λ) = ρ_t(λ)/π, the general model of R_rs above the water surface is as follows:

R_rs(IOPs, H) = R_ap(SGR, wind, IOPs, H, R_b) − R_sg(SGR, wind) − R_bw(IOPs, H, R_b)

(2)

where H, SGR, winds, and R_b represent the water depth, the sun–satellite geometric relationship, the wind (i.e., speed and aspect), and the bottom reflectance, respectively. The symbols used in this paper are shown in

Table 2. The symbols and abbreviations for the parameters used in this paper.

Symbols and Abbreviations	Description	Units
aw	Absorption coefficient of pure water	m−1
aphy	Absorption coefficient of algal pigments	m−1
aNAP	Absorption coefficient of non-algal pigments	m−1
aCDOM	Absorption coefficient of CDOM	m−1
a	Absorption coefficient of the total (=aw + aphy + aNAP + aCDOM)	m−1
b	Scattering coefficients	m−1
B	Backscattering ratio
bb	Backscattering coefficients	m−1
Chl-a	Chlorophyll-a concentration	ug/L
β(ψ,λ)	Volume scattering function (VSF)
f/Q	Water column bidirectional factor	sr−1
D	Bottom slope	deg
D′	Bottom aspect	deg
H	Water depth	m
tgas(λ)	Gaseous transmittance	m−1
t(λ)	Total diffuse atmospheric transmission	m−1
ρt(λ)	Irradiance reflectance just above the water surface
ρA(λ)	Aerosol reflectance
ρr(λ)	Molecular Rayleigh scattering reflectance
Rrs	Remote sensing reflectance	sr−1
Rap	Apparent remote sensing reflectance	sr−1
Rsg	Surface specular reflectance	sr−1
Rcw	Remote sensing reflectance from water-column scattering	sr−1
Rbw	Remote sensing reflectance from bottom reflectance	sr−1
θ0	Solar zenith angle above water surface	deg
θ	View angle above the water surface	deg
θ0′	Subsurface solar zenith angle	deg
θ′	Subsurface view zenith angle	deg
θn	Water surface slope gradient	deg
μ0	Cosine of θ0
μ	Cosine of θ
μn	Cosine of θn
tu	Water–air transmittance	sr−1
td	Air–water transmittance	sr−1
L	Optical path-elongation factor	m−1
k	Beam attenuation coefficient (=a + bb)	m−1
φ	Relative azimuth angle	de
Θ	Sun–satellite relative phase angle	deg
σu	Mean square slope in an upwind direction	deg
σc	Mean square slop in the crosswind direction	deg
λ	Wavelength	nm
W	Wind speed	m/s
△φ	Wind direction	deg
n	Water refractive index (≈1.34 in seawater)
SYB	Sanya Bay
CDOM	Coloured dissolved organic matter	m−1
Chl-a	Chlorophyll-a
CZCS	Coastal zone color scanner
IOPs	Inherent optical properties
MODIS	Moderate resolution imaging spectroradiometer
2Br	Two-band-ratio
BGr	Blue-green band-ratio
DCI	Difference Chl-a index
NAP	Suspended mineral particle	mg/L
NIR	Near-infrared
NRr	NIR-red band-ratio
RGr	Red-green band–ratio
RSRF	Relative spectral response function
RTM	Radiative transfer model
SeaWiFS	Sea-viewing wide field of view sensor
SAA	Solar azimuth angle	deg
SGR	Sun–satellite geometric relationship	deg
SZA	Solar zenith angle	deg
VZA	View zenith angle	deg
VAA	View azimuth angle	deg

.

3.1. Sun Glint Effect

The view zenith angle and the view azimuth angle of traditional ocean color satellites (e.g., CZCS, MODIS, and SeaWiFS) are designed to avoid the sun glint. For land observation sensors (e.g., SPOT6, Landsat8, Aster, and GaoFen2), the acquired images may be located in sun glint areas, and so R_sg needs to be removed.

The water column has strong absorption in the NIR band of SPOT6 (i.e., 825 nm); therefore, the approach of wiping off the sun glint with the NIR band has been used in remote sensing reflectance retrieval of oceanic waters [30]. However, this simple method is not available for coastal waters with high concentrations of suspended matter. Thus, we chose the physical Cox–Munk anisotropy model [25,26] to calculate the sun glint of SPOT6.

The Cox–Munk model is the most classical model for simulating sun glint, and the anisotropic form of the Cox–Munk model is still the most widely used method based on the research of Zhang [27]. Sun glint may exist in the image of satellites in large view zenith angle (VZA) [28], but the sun glint pollution on the image of land observation satellites is easily ignored because there are no obvious sun glint distribution features when the view angle of each pixel in an image changes little [49].

3.2. Water-Leaving Reflectance Contributed by Bottom Reflectance (R_bw)

For optically shallow waters, the effect of R_bw is a non-negligible factor. R_bw is a function that depends on the bottom reflectance (R_b), the diffuse attenuation coefficient, and the water depth. According to the RTM [32], R_bw is written as follows:

$$R_{\mathrm{bw}}=\frac{t_{\mathrm{u}}·t_{\mathrm{d}}}{n^2}·R_{\mathrm{b}}·{\mathrm{e}}^{-\left(k_{\mathrm{u}}+k_{\mathrm{d}}\right)H}$$

(3)

where t_u, t_d, n, k_u, k_d, and H are the water–air transmittance, the air–water transmittance, the refractive index, the upwelling diffuse attenuation coefficient for the bottom, the downwelling diffuse attenuation coefficient for the bottom, and the water depth, respectively. Diffuse attenuation coefficients (k_u and k_d) depend on the SGR, a, and b_b, and can be simulated by HydroLight [35] (described in Appendix A). The R_b depends on the SGR underwater surface and bottom albedo, and can be calculated by an empirical bottom bidirectional reflectance distribution fucntion (BRDF) model derived from datasets of measured bottom remote sensing reflectance [49]:

$$R_{\mathrm{b}}({\theta }_0^{\mathrm{s}},{\phi }^{\mathrm{s}},r_{\mathrm{b}})=x_1+x_2·{\theta }_0^{\mathrm{s}}+(x_3+x_4·{\theta }_0^{\mathrm{s}})·\exp (-(x_5+x_6·{\theta }_0^{\mathrm{s}})·{\phi }^{\mathrm{s}})·r_{\mathrm{b}}$$

(4)

where the parameters ${\theta }_0^{\mathrm{s}}$, φ^s, and r_b represent the incident angle at the slope, the relative azimuth angle at the slope, and the bottom albedo, respectively. The empirical parameters, x_i, are based on six types of bottoms, and the parameters of grapestone sand (similar to the bottom type in Sanya) were adopted. The x_1–6 were −0.341, 0.00594, −0.853, 0.0436, 0.00884, and 0.000696, respectively. The parameter of ${\theta }_0^{\mathrm{s}}$ is the function of the bottom terrain and the SGR underwater surface:

$${\theta }_0^{\mathrm{s}}=\frac{1}{2}\arccos (\cos {\theta }_0^{\mathrm{w}}\cos {\theta }^{\mathrm{w}}+\sin {\theta }_0^{\mathrm{w}}\sin {\theta }^{\mathrm{w}}\cos ({\phi }^{\mathrm{s}}))-D$$

(5)

where ${\theta }_0^{\mathrm{w}}$, θ^w, and D represent the incident angle underwater subsurface, the view angle underwater surface, and the bottom slope, respectively; the parameter φ^s = φ − D′, where D′ is the bottom aspect.

4. Chl-a Concentration Retrieval Algorithm

The effects of sun glint, bottom reflectance, and NAP have a great impact on Chl-a concentration retrieval in optically shallow zones of coastal waters using broad-band satellites. In this paper, a semi-empirical algorithm is proposed to reduce the effects of these limitations, and SPOT6 data serve as the reflectance data. The general flow chart of the Chl-a concentration retrieval process is shown in Figure 4.

4.1. Review of the Chl-a Retrieval Algorithm

According to the IOPs in different types of waters, there are mainly three kinds of empirical or semi-empirical algorithms for Chl-a retrieval: two-band ratio algorithms (e.g., BGr and red-green band ratio, RGr), spectral band difference algorithms (e.g., Fluorescence Line Height (FLH), color index (CI), and maximum chlorophyll index (MCI)). The selection of bands and methods for Chl-a retrieval needs to consider the positions of the Chl-a absorption peak and the interferences from CDOM and NAP. The absorption peaks of Chl-a are mainly located near 440 nm, 660 nm, and 705 nm, and the absorption peak of Chl-b (near 490 nm) has been used to replace 440 nm to reduce the effect of CDOM. The excellent correlation between Chl-a and chlorophyll-b (Chl-b) ensures the feasibility of this replacement. The algorithms based on an absorption peak at 705 nm are proposed to reduce the effect of NAP.

The principle of the QAA algorithm is to obtain the Chl-a absorption coefficient by the IOPs calculated by the reflectance, and the principle of the empirical or semi-algorithms are based on the spectral features of waters with different Chl-a concentrations. On the basis of the band settings of SPOT6 data (red band: 450–520 nm; green band: 530–590 nm), the algorithms of BGr, CI, and RGr were adopted to estimate the concentration of Chl-a in this paper. For reducing the effect of sun glint and bottom reflectance, we also tried another form of the blue-green band algorithm to retrieve the Chl-a concentration, and this algorithm is called the difference Chl-a index (DCI). The DCI algorithm is an empirical algorithm based on the difference of reflectance between blue band and green band. The algorithms of BGr, CI, and DCI are based on an absorption peak at blue band and RGr is based on an absorption peak at red band. The common Chl-a retrieval algorithms are shown in

Table 3. The common algorithms and spectral bands (in nanometers) for Chl-a concentration retrieval.

Algorithm	Sensor	Band	Water Type	Reference
BGr	OLI	483, 562	Eutrophic inland	[23,24]
RGr	Sentinel-2	665, 560	Eutrophic inland	[13]
OC3M	MODIS	443, 488, 555	Case-2 coastal	[12]
QAA	SeaWiFS	443, 490, 550, 667	Case-2 coastal	[10,11]
FLH	MODIS	667, 678, 746	Eutrophic coastal	[50]
CI	MERIS	443, 555, 670	Oligotrophic	[51]
MCI	MERIS	681, 709, 753	Eutrophic coastal	[52,53]

. For reducing the effect of NAP on empirical algorithms, we propose an improved semi-empirical algorithm for Chl-a concentration algorithm and describe in the next section.

4.2. Improved Chl-a Concentration Algorithm

Chl-a concentration retrieval using remote sensing methods is based on the relationship between IOPs and the remote sensing reflectance. For optically deep waters, R_rs is the function of IOPs and the SGR [54,55]. However, for optically shallow waters, R_rs is also a function of water depth. On the basis of the RTM, the R_rs can be written as follows [32,35]:

$$R_{\mathrm{rs}}(\mathit{IOPs},f,Q)=\frac{t_{\mathrm{u}}·t_{\mathrm{d}}}{n^2}·\frac{f}{Q}·r_{\mathrm{rs}}\left(a,b_{\mathrm{b}}\right)·\left[1-\exp \left(-\left(k_{\mathrm{u}}+k_{\mathrm{d}}\right)·H\right)\right]$$

(6)

where r_rs and the parameters r_rs, f/Q, a, and b_b are the subsurface remote sensing reflectance, the bidirectional factor for reflectance just above the water surface, the total absorption coefficient, and the total backscattering coefficient, respectively. The polynomial, $\left[1-\exp \left(-\left(k_{\mathrm{u}}+k_{\mathrm{d}}\right)·H\right)\right]$ in Equation (5), is the effect of the water depth and needs to be removed. The parameter f/Q is the function of the SGR, the winds, and the water constituents, and it can be calculated by the experimental function simulated by the HydroLight dataset [35]. In the absence of the value of f/Q, the subsurface remote sensing reflectance (r_rs) can be calculated by the QAA when the solar zenith angle (SZA) is larger than 20° [15]:

$$r_{\mathrm{rs}}=\frac{R_{\mathrm{rs}}\left(\mathsf{\lambda }\right)}{0.52+1.7R_{\mathrm{rs}}\left(\mathsf{\lambda }\right)}$$

(7)

The relationship between r_rs and the ratio of the total backscattering coefficient to the total absorption coefficient u(λ) can be written as follows:

$$u(\mathsf{\lambda })=\frac{-{\mathrm{g}}_0+\sqrt{{\left({\mathrm{g}}_0\right)}^2+4{\mathrm{g}}_1·r_{\mathrm{rs}}\left(\mathsf{\lambda }\right)}}{2{\mathrm{g}}_1}$$

(8)

The ratio of u(λ₁) to u(λ₂), (u_ratio) can be used to retrieve the concentration of Chl-a instead of the ratio of r_rs(λ₁) to r_rs(λ₂). According to the relationship between IOPs and u(λ), u(λ) = b_b(λ)/(a(λ) + b_b(λ)), the u_ratio can be written as follows:

$$u_{\mathrm{ratio}}=\frac{u\left({\mathsf{\lambda }}_1\right)}{u\left({\mathsf{\lambda }}_2\right)}=\frac{b_{{b}}\left({\mathsf{\lambda }}_1\right)}{b_{{b}}\left({\mathsf{\lambda }}_2\right)}·\frac{a\left({\mathsf{\lambda }}_2\right)}{a\left({\mathsf{\lambda }}_1\right)}$$

(9)

The main constituents affecting IOPs are the water molecules, Chl-a, NAP, and CDOM, and the parameters of a and b_b can be written as follows [56,57]:

$$a=a_{\mathrm{w}}+a_{\mathrm{phy}}^{*}·\mathrm{Chl}-\mathrm{a}+a_{\mathrm{NAP}}^{*}·\mathrm{NAP}+a_{\mathrm{CDOM}}$$

(10a)

$$b=b_{\mathrm{b}\mathrm{w}}+b_{\mathrm{bphy}}^{*}·\mathrm{Chl}-\mathrm{a}+b_{\mathrm{bNAP}}^{*}·\mathrm{NAP}$$

(10b)

The pure seawater backscattering (b_bw) is constant and much less than the particle backscattering (b_bp = b_bNAP + b_bphy) [58,59], and the particle backscattering coefficient (b_p) changes little in the visible and NIR bands [35,57,60,61,62]. To reduce the model complexity, u_ratio can be converted to the following [8]:

$$u_{\mathrm{ratio}}\approx \frac{a\left({\mathsf{\lambda }}_2\right)}{a\left({\mathsf{\lambda }}_1\right)}$$

(11)

To reduce the effect of a_NAP, a_w, and a_cdom, Equations (8) and (10b) can be merged and written as follows:

$$a_{\mathrm{phy}}({\mathsf{\lambda }}_2)=\frac{u\left({\mathsf{\lambda }}_1\right)}{u\left({\mathsf{\lambda }}_2\right)}·a\left({\mathsf{\lambda }}_1\right)-[a_{\mathrm{w}}({\mathsf{\lambda }}_2)+a_{\mathrm{NAP}}({\mathsf{\lambda }}_2)+a_{\mathrm{CDOM}}({\mathsf{\lambda }}_2)]$$

(12)

For oceanic waters or phytoplankton-dominated waters, a_phy is approximately equal to the particle absorption coefficient (a_p). For coastal waters of SYB in winter, a_NAP is an important parameter, whose proportion is over 80% in the a_p (as shown in Table 1). We thus propose the semi-empirical algorithm to retrieve the concentration of Chl-a by absorption of a_phy. The parameters of a_w, a_NAP, a_CDOM, and a are unknown. Therefore, a_w was adopted by the measured data from Pope and Fry [63]; a_NAP was calculated by the semi-analytical model proposed by Tassan, as shown in Appendix A.3; and a_CDOM is calculated by an empirical algorithm as shown in Appendix A.4. The Chl-a concentration was estimated by the empirical algorithm (e.g., DCI, BGr, and RGr). The parameter a(λ) is the sum of a_w, a_NAP, a_CDOM, and a_phy (as shown in Equation (10a)), and the semi-empirical algorithm uses the parameters of a, a_NAP, a_CDOM, and a_w to calculate a_phy (as shown in Equations (11) and (12)).

5. Results

5.1. SPOT6 Data Atmospheric Correction

The atmospheric correction of remote sensing data is the first step in Chl-a retrieval and may remove over 70%–90% of the signal measured by the sensor in the visible bands [27,28,64]. The fast line of sight atmospheric analysis of spectral hypercube (FLAASH) module is a widely used atmospheric correction module for multi-band remote sensing data and is based on MODTRAN4 (moderate spectral resolution atmospheric transmittance algorithm and computer model) [64,65,66]. FLAASH has been successfully used for the atmospheric correction of SPOT6 image data [67,68], and the input parameters are shown in

Table 4. Input parameters in fast line of sight atmospheric analysis of spectral hypercube (FLAASH).

Sensor Type	SPOT 6	Atmospheric Model	Tropical
Sensor altitude	695 km	Aerosol model	Maritime
Ground elevation	0 m	Initial visibility	30 km
Flight data	21 January 2013	View zenith angle	22°
Flight time	10:58	View azimuth angle	135°

. The meteorological data were retrieved from the National Meteorological Information Center: http://data.cma.cn/data/cdcindex/cid/6b2344874db4bbaf.html. As we did not have in situ reflectance measurements above the water surface in the SYB, we collected the reflectance of the wet and dry sand for calibration of the atmospheric correction (as shown in Figure 5). We also applied the FLAASH module to the Landsat8 OLI (imaging time: 15 June 2013) and used the in situ reflectance data from coastal waters in Qinghuangdao (sampling in 119.80° E, 39.95° N on 15 June 2013) to verify the calibration of the atmospheric correction (as shown in Figure 6). The reflectance data were measured by the FieldSpec 3 Spectroradiometer (America, ASD Company, Alpharetta, GA, USA).

5.2. Parameter u(λ) Correction

The a_phy was calculated by u(λ), as shown in Equation (12). Therefore, u(λ) was calculated by R_ap, R_rs, and r_rs to show the effects of sun-glint and bottom. In Figure 7, the green circles, the red diamonds, and the blue rectangles represent the u(λ) calculated by R_ap, R_rs, and r_rs, respectively. The effects of sun glint (as shown in green circles) and bottom reflectance (as shown in blue rectangles) increase the value of u(λ). The effect of sun glint is larger than the bottom reflectance, and the error caused by sun glint may be over 20%. After model correction, the u(λ) calculated by the r_rs of the SPOT6 data and those calculated by the in situ data have good agreement.

5.3. Comparison of the Different Algorithms

For SPOT6 data, there is no available analytical or semi-analytical algorithm, so we tried different empirical algorithms to retrieve the concentration of Chl-a. According to the absorption characteristics of Chl-a and the relative response function (RSRF) of SPOT6, the BGr, RGr, CI, and DCI algorithms [8,10,12,13,51] were applied to the remote sensing reflectance processed by atmospheric correction using SPOT6 data. The algorithm correlation coefficients are shown in Figure 8.

Although the correlation between RGr and Chl-a concentration is better than the other three algorithms, RGr and Chl-a concentration should be positively correlated rather than negatively correlated based on the absorption characteristics of Chl-a. NAP and RGr are positively correlated, and the NAP signal masks the Chl-a concentration signal in waters with a high concentration of NAP using the RGr algorithm. Thus, we applied the BGr and DCI on R_rs (i.e., R_ap without the effect of sun glint) to retrieve the concentration of Chl-a (Figure 9).

In the same image, the sun glint decreases as the wavelength gets longer owing to the variation of the water refractive index, but the sun glint changes less than 5% in the visible bands [30]. If the sun glint is small, the sun glint effects on BGr or DCI can be ignored, and the DCI has stronger resistance to sun glint (as shown in Figure 8a,d and Figure 9). The value of the SPOT6 image sun glint calculated by the Cox–Munk model is 0.01 in the blue band; thus, the effect of the sun glint is very small for the DCI algorithm.

The effect of the bottom on Chl-a concentration retrieval in optically shallow waters depends on the bottom reflectance and the attenuation coefficient. The attenuation coefficients in the visible bands are less than those in the NIR bands, and the effect of R_bw on the blue and green bands needs to be considered. Better correlation coefficients were obtained by the BGr and DCI algorithms without bottom reflectance, and a relatively good correlation coefficient (R² = 0.75) was obtained by the DCI algorithm (as shown in Figure 10). The accuracy of global Chl-a concentration retrieval based on the MODIS algorithm reached R² = 0.85, and thus the accuracy of the DCI algorithm still needs to be improved.

For the BGr and DCI algorithms, the effect of NAP may be the most important factor for Chl-a concentration retrieval, because the absorption of NAP in the blue and green bands is much larger than the absorption of Chl-a. We applied the semi-analytical algorithm based on a_phy (490) to reduce the effect of NAP, and the accuracy of algorithm was improved (R² = 0.91 as shown in Figure 11). The retrieval results using a_phy (660) inverted by the semi-analytical algorithm were not as good as those using a_phy (490), which may be because the wavelength range of the SPOT6 red band (625–695 nm) contains not only a strong Chl-a spectral absorption, but also a weak Chl-a spectral reflection (as shown in Appendix A.1).

According to the comparison of different algorithms (as shown in Figure 8, Figure 9, Figure 10 and Figure 11), the algorithms based on BGr, DCI, and a_phy (490) are available. We compared the Chl-a retrieval results by different reflectance using different algorithms to analyze the performance of the semi-empirical algorithm proposed in this paper (as shown in Figure 12). This semi-empirical algorithm can be used to reduce the effects of boundary conditions (i.e, sun glint and bottom) and NAP.

5.4. Statistical Analysis

A number of statistical methods were used in this paper to compare the different algorithms, the coefficient of determination (R²), the root mean square error (RMSE), the mean absolute percentage error (MAPE), and the relative error (RE).

$${\mathrm{R}}^2=1-\frac{\sum {\left(y_{\mathrm{measured}}-y_{\mathrm{modeled}}\right)}^2}{\sum {\left(y_{\mathrm{measured}}-\overline{y_{\mathrm{modeled}}}\right)}^2}$$

(13)

$$\mathrm{MAPE}=\frac{1}{m}\sum 100\left|\frac{y_{\mathrm{measured}}-y_{\mathrm{modeled}}}{y_{\mathrm{measured}}}\right|$$

(14)

$$\mathrm{RMSE}=\sqrt{\frac{1}{m}{\left(y_{\mathrm{measured}}-y_{\mathrm{modeled}}\right)}^2}$$

(15)

$$\mathrm{RE}=\left|\frac{y_{\mathrm{measured}}-y_{\mathrm{modeled}}}{y_{\mathrm{measured}}}\right|\times 100\%$$

(16)

where m is the sample number, y_measured is the observed value, y_modeled is the predicted value, and $\overline{y_{\mathrm{modeled}}}$ is the expectation of the observed value, where y is the Chl-a concentration calculated by different algorithms using different reflectance of R_ap, R_rs, and r_rs. The empirical coefficients of the algorithms in this paper are shown in

Table 5. The statistics of algorithms in this paper. RMSE, root mean square error; MAPE, mean absolute percentage error.

Algorithm	Data	Empirical Coefficients	R2	RMSE	MAPE
BGr	Rap (λ)	y = 67.85x2 − 163.9x + 99.31	0.72	0.2255	41.72%
		Y = −10x + 12.35	0.68	0.2565	66.23%
	Rrs (λ)	y = 42.98x2 − 108.9 + 68.18	0.74	0.2243	36.28%
		y = −7.66x + 9.859	0.7	0.2444	50.46%
	rrs (λ)	y = 47.85x2 − 116.9 + 71.78	0.75	0.2223	40.93%
		y = −7.57x + 9.5	0.72	0.2422	41.56%
DCI	Rap (λ)	y = −8901x2 − 1042x + 2.061	0.76	0.2151	36.09%
		y=-245x + 2.507	0.71	0.2277	45.92%
	Rrs (λ)	y = −8901x2 − 1025x + 2.04	0.76	0.2151	36.04%
		Y = −245x + 2.475	0.71	0.2277	45.76%
	rrs (λ)	y = −4478x2 − 558.2x + 1.843	0.77	0.2141	36.15%
		y = −131.3x + 2.057	0.75	0.2253	38.97%
aphy (490)	rrs (λ)	y = 43.3x1.454	0.91	0.11	21.51%

.

5.5. Sun Glint Effect

We used polarization coordinates to show the sun glint distribution and the radius and polar angle present the view zenith angle and relative azimuth angle (solar direction = 180°), respectively. Although the SPOT6 image is not located in the center of the sun glint (red circle represents the view direction of SPOT6, as shown in Figure 13), it is a non-negligible signal (≈0.01 sr⁻¹) because R_rs is very small. The scope of the sun glint zone expands as the wind speed increases, and the SPOT6 image will locate in the sun glint zone if the wind speed is strong enough. Figure 13 shows the sun glint distribution, and Figure 14 shows that the reflectance at 825 nm is close to the R_sg calculated by the Cox–Munk model when the concentration of NAP is low, which proves that the Cox–Munk model is still reliable when the image is on the edge of the sun glint area with a low value.

In our study, R_rs ranged from 0.03 to 0.05 sr⁻¹ in the visible bands and the proportion of R_sg in R_rs ranged from approximately 20% to 30%. Hence, we simulated the sensitivity of the two-band algorithm (including BGr, RGr, and DCI) to R_sg to analyze the effect of sun glint on Chl-a concentration retrieval (Figure 15). The remote sensing reflectance (R_ap) at 490 nm and 560 nm was simulated by Equations (1) and (5), and the parameters NAP = 1 mg/L and Chl-a = 0.3–10 were the model input parameters. In order to control the variables, the effects of bottom reflectance and NAP were not considered in this analysis.

As the value of the sun glint goes up, the ratio of the two-band reflectance approaches 1 and the variations of the two-band ratio (2Br) in different Chl-a concentrations are reduced. This phenomenon reduces the sensitivity of R_sg to the variation of Chl-a concentration and causes the 2Br algorithm to lose its ability to retrieve Chl-a. Relative to the 2Br, the DCI has better resistance ability to sun glint. However, we also need in situ data to correct the algorithm because the R_sg increases the value of DCI.

According to the equation of BGr ($BGr=\frac{R_{\mathrm{rs}}\left(490\right)+R_{\mathrm{sg}}\left(490\right)}{R_{\mathrm{rs}}\left(560\right)+R_{\mathrm{sg}}\left(560\right)}$), the relationship of BGr will not be affected by sun glint when BGr = 1. The Chl-a retrieval results using the value of BGr = 1 against the measured Chl-a data may help to verify whether the retrieval result is affected by sun glint.

5.6. Bottom Effect

To study the effect of the bottom on Chl-a concentration retrieval, we subtracted the results retrived by r_rs (i.e., the remote sensing reflectance after removing the effect of bottom) from those by R_rs (i.e., the remote sensing reflectance affected by the bottom). Figure 16 shows that the error caused by the effect of bottom is approximately 0–0.3 ug/L. The bottom in shallow waters causes larger Chl-a concentratio retrieval results because the water-leaving reflectance contributed by a sandy bottom in the green band is larger than that of the blue band. The error spatial distritbution shows that error in the waters with H < 2 m is low, which may be beacause the effect of the bottom not only changes with water different depths, but also in terms of turbidity. The waters near the coast contain a high concentration of NAP picked up by the sea waves. To analyze the relationship between errors caused by bottom reflectance and H, we select three hatch lines (blue lines in Figure 16). Figure 17 shows that error between H = 0 and 3.8 m is exponential rise and that H > 3.8 m is exponential decline, and the error maximum occurs at H = 3.8 m.

According to the RTM, the remote sensing reflectance contributed by R_bw depends on the bottom albedo, the water constituent, and the wavelength. A bottom of sandy coastline is mainly fine-grained sand (Figure 18a), while that of the rocky coastline is mainly coarse quartz arkose (Figure 18b), and the bottom of corals include health corals and bleaching corals donated by Chen’s research [69] (Figure 18c,d). We collected the spectra of different types of bottoms to analyze the contribution of R_bw to R_rs.

The spectra of fine-grained sand and coarse quartz arkose are similar, and the spectra of stones are slightly higher (Figure 18a,b). Compared with sand and coarse quartz arkose, the coral spectrum has a strong absorption at 675 nm, and the coral reflectance is relatively flat until 675 nm and then rises sharply after 675 nm (Figure 18c,d). The bottom reflectance generally increases as the wavelength increases in the visible and near-infrared bands, and we analyzed R_bw in different types of bottoms and at different water depths and wavelengths. The input parameters included the attenuation coefficient, H, λ, and SGR. The IOP data were averaged from the in situ samples, and the VZA, SZA, and view azimuth angle (VAA) were obtained by SPOT6 data. The effect of the bottom on the remote sensing reflectance is shown in Figure 19.

Although the blue light has the best penetration ability, the maximum of R_bw occurs in the red band. The R_bw in the NIR band is the lowest, the water depths that affect the remote sensing reflectance are less than 2 m, and the near-infrared band has the potential to retrieve the parameters of the water column in optically shallow waters. The green lines show the water depths when the ratio of R_bw to remote sensing reflectance is equal to 5%, and the variations of the water depths on the green lines are independent of bottom types. In coastal waters, the R_bw depends on the bottom albedo when the water depth is less than 2 m. In fact, the error caused by bottom reflectance with H < 2 m is low owing to the high concentration of NAP. Therefore, the R_bw mainly depends on the attenuation of the water column, and we chose four types of waters to analyze the effect of R_bw on 2Br. In coastal waters with a low degree of pollution, the water constituent is mainly affected by terrigenous matter, and we classified coastal waters into four types based on the following [23,70,71]: (1) ultraoligotrophic waters with Chl-a < 0.3 ug/L and NAP < 3 mg/L; (2) oligotrophic waters with Chl-a = 0.3–1 ug/L and NAP = 3–5 mg/L; (3) mesotrophic waters with Chl-a = 1–5 ug/L and NAP = 5–10 mg/L; and (4) eutrophic waters with Chl-a = 5–10 ug/L and NAP = 10–20 mg/L. We simulated the two-band reflectance ratio and difference by water depths and water constituents to analyze the sensitivity of the parameters (Figure 20).

The R_bw varies with the water column attenuation coefficient and bottom albedo, and it is difficult to retrieve the Chl-a concentration without information of the R_bw. Figure 20 shows that the DCI is less affected by the bottom than BGr and RGr, and the optical depth for DCI is less than 2 m in clear waters (type 1 waters) and 1 m in turbid waters (type 4 waters). According to the empirical algorithms, the R_bw increases the Chl-a concentration by approximately 0–4.2 ug/L and 0–1.2 ug/L using BGr and DCI, respectively.

5.7. The Effect of NAP on Spatial Distribution

To study the effect of NAP on the distribution of the Chl-a concentration, we separated the results between those retrieved by the r_rs and those by a_phy (490) calculated using the semi-empirical algorithm (as shown in Figure 21). Similar to the effect of the bottom, NAP also causes larger retrieval results because the absorption coefficient of NAP in the blue band is larger than that in the green band. Figure 21 shows that the errors caused by NAP range from approximately 0.05 to 0.45 ug/L. To analyze the effect of NAP on the Chl-a retrieval result, we calculated the concentration of NAP by the empirical algorithm using NIR band (as show in Appendix A.3). Figure 22 shows the distribution of NAP concentration in SYB, and the range is approximately 2.2–30.1 mg/L. Shallow waters contain more NAP, and the waters with H < 5 m are most turbid (NAP > 25 mg/L). Open waters contain low NAP, and the waters with H > 35 m are clear (NAP < 5 mg/L). We also selected the results from three hatch lines in Figure 19 to analyze the relationship between the error caused by NAP and the concentration of NAP. Figure 23 shows that the effect of NAP is positively correlated with the concentration of NAP.

The two-band algorithm based on the visible bands is widely used in Chl-a concentration retrieval, but it is largely affected by NAP. Figure 24 shows the effect of NAP on the recognition ability of two-band algorithms of Chl-a. Figure 21a shows that the BGr algorithms have a stronger anti-interference ability for NAP than that of DCI and RGr in waters with NAP < 5 mg/L, and the recognition ability of these three algorithms of Chl-a is not good in waters with NAP > 5 mg/L. RGr and DCI lose their recognition ability of Chl-a when the concentration of NAP is larger than 9.7 mg/L and 12 mg/L, respectively. Thus, we used the semi-empirical algorithm in this paper to reduce the effect of NAP.

5.8. Chl-a Retrieval of SYB

The Chl-a concentration of SYB was calculated using a_phy (490) with the semi-empirical algorithm (Figure 25). A high concentration of Chl-a (>1.7 ug/L) is located in the areas of water depth (H) <5 m, and a medium concentration (0.8–1.7 ug/L) distributes in the areas of 5 m < H < 10 m, while the concentration is low (<0.8 ug/L) in the areas of H > 15 m. In the southwest of SYB, the Chl-a concentration is lower than 0.5 ug/L in the areas of H > 20 m, and the water type is similar to oceanic water. The concentration of Chl-a is concentrated between 0.5 and 1.1 ug/L (as shown in Figure 26), and the areas with the highest concentration of Chl-a (>2.5 ug/L) are concentrated in harbor fishing ports, docks, and shores, while the Chl-a concentration is lowest at water depths deeper than 30 m. The highest concentration of Chl-a occurs in the east of the Luhuitou Peninsula owing to the winter current flows from the northwest to the southeast and the barrier of the peninsula. The Chl-a concentration in the east of the study area is larger than that in the west, which is caused by the terrain and the current. The West Island and the Luhuitou Peninsula block the impact of the current and protect the coral area in the west of the island and peninsula. On the whole, the water in SYB is clear, and tourism activities such as diving are more suitable for economic activities in SYB than aquaculture in winter.

6. Discussion

We compared the results using the algorithm in this paper with the MODIS Ocean Color3 (OC3M) algorithm, which has been widely used for oceanic and coastal waters [72]. Figure 27 shows that the algorithm used in this paper is better than that of OC3M, and the RE in oligotrophic waters (Chl-a < 1 ug/L) is greater than that in mesotrophic waters. Therefore, we needed to analyze the effect of interference factors (i.e., sun glint, bottom, and NAP) on the Chl-a retrieval algorithm. According to the water samples collected in the study area, the water types are oligotrophic and mesotrophic. Hence, we compared four common empirical algorithms to choose the best algorithm as the start of the algorithm in this paper.

Figure 7 shows that the effect of R_sg on u(λ) is larger than that of bottom reflectance, and thus it is necessary to restrain the sun glint. We used the Cox–Munk model to calculate the R_sg and Figure 13 shows that the angle of the satellite observation may be affected by sun glint when the wind is strong enough. Figure 14 also shows that the R_rs of clear water in the NIR bands can replace the R_sg calculated by the Cox–Munk model, and Figure 8 shows that DCI has the ability to restrain the sun glint effect. Figure 9 shows that, after removing the effect of R_sg, both DCI and BGr have a relatively good coefficient with the Chl-a concentration.

For optically shallow waters, R_bw is non-negligible. We collected six types of bottoms and analyzed the effect of R_bw on R_rs. Figure 16 shows that R_bw in the red band is higher than that in the blue band, and so R_bw causes the abnormal rise in Chl-a concentration using BGr or DCI. Figure 9 also shows that DCI is better than BGr when R_rs is affected by the bottom. R_bw decreases as H gets deeper, and is lower than 5% of R_rs when H > 4.5 m. In fact, the maximum R_bw does not occur where the water depth is the lowest. Large amounts of re-suspended sediment near the coastline increase the turbidity. Figure 16 and Figure 17 show that retrieval is most affected by bottoms at H = 3.8 m and the biggest error is approximate 0.3 ug/L.

There are more data sources from satellites for land observation than there are for coastal zone observation, and land observation satellites have the potential to be applied in coastal zones. For Chl-a retrieval in oligotrophic or mesotrophic waters, algorithms based on the blue-green band may be more suitable than algorithms based on the NIR bands [14]. Therefore, multi-spectral land observation satellites with the blue-green band can also be used for Chl-a concentration retrieval. NAP is a big barrier for Chl-a concentration retrieval using observation satellites. Thus, we used a semi-empirical algorithm to reduce the effect of NAP and to calculate a_phy to retrieve the concentration of Chl-a. Figure 11 and Figure 12 show that the Chl-a concentration retrieval results calculated by the semi-empirical algorithm are better than those calculated by the empirical algorithms. Figure 21 also shows that the errors caused by NAP range from 0.05 to 0.45 ug/L and Figure 23 shows that the errors caused by NAP have a good correlation with the concentration of NAP.

In order to study which parameter has a greater influence on Chl-a inversion, we simulated the errors caused by sun glint, bottom reflectance, and NAP. We set the sand spectrum, Chl-a = 1 ug/L, H = 3–7 m, and R_sg = 0.01 sr⁻¹ as the input parameters. Figure 28 shows that the error caused by bottom reflectance is lower than 0.1 ug/L when H > 5 m, and the error caused by sun glint in our study is approximately 0.11 ug/L. NAP is the biggest interference factor when the concentration of NAP > 9.5 mg/L. Figure 25 also shows that the concentration of NAP in the bay mouth is higher than 15 mg/L, and is lower than 5 mg/L in H > 30. Hence, the R_bw of H < 5 m should be considered and the NAP is the biggest interference factor in the semi-enclosed bay.

In winter, the Luhuitou Peninsula on the east side of SYB blocks out suspended matter brought by the southeast current, and the semi-enclosed terrain keeps in more suspended matter and provides more nutrients for phytoplankton near the shore. However, the water of SYB is not highly nutritive in winter (Figure 25). According to the Chl-a retrieval algorithm (i.e., Chl-a = A·a_phy^B), the empirical coefficients of A and B change with the $a_{\mathrm{phy}}^{*}$, which changes with the types of phytoplankton. The Lee’s dataset taken over the shallow waters along the west Florida shelf, Key West, Fla., and the Bahamas showed that $a_{\mathrm{phy}}^{*}$(490) ranged from 0.03 to 0.17 m² mg⁻¹ [32], and the Neil’s dataset measured from Irish Sea and the Bristol Channel shows that $a_{\mathrm{phy}}^{*}$(490) ranges from 0.02 to 0.22 m² mg⁻¹ [46]. Therefore, the change in the types of phytoplankton may be another source of error in this algorithm.

7. Conclusions

We proposed a semi-empirical Chl-a retrieval algorithm for SPOT6 data as a demonstration for multi-spectral data, and we applied it to the waters in SYB. A comparison of the performance of this semi-empirical algorithm with that of the OC3M algorithm, and a concurrent measurement of the concentration of Chl-a from 16 surface water points in SYB, demonstrated the appropriateness of the semi-empirical algorithm to estimate the Chl-a concentration in coastal waters using SPOT6 data (R² = 0.93 with an RMSE = 0.11 ug/L and MAPE = 21%).

Through a comparison of the Chl-a concentration retrieval results using different algorithms and the sensitivity of R_sg, R_bw, and NAP to two-band algorithms (e.g., BGr, RGr, CI, and DCI), algorithms that use the blue-green band are available for oligotrophic and mesotrophic coastal waters. DCI is a better empirical algorithm than BGr because it has a stronger anti-interference ability in terms of sun glint and bottom reflectance. To improve the accuracy of the empirical algorithm, we calculated the total absorption coefficient and the absorption coefficient of each component by empirical algorithms, and then created the Chl-a concentration retrieval semi-empirical algorithm in the form of Chl-a = A[a_phy(490)]^B, where A = 43.3 and B = 1.454. The empirical coefficients of A and B can be calculated by the in situ data of the Chl-a concentration or $a_{\mathrm{phy}}^{*}$(490), which changes with the types of phytoplankton.

In the case of ignoring R_sg unrelated to IOPs, NAP is the biggest interference factor when the concentration of NAP > 9.5 mg/L and the effect of bottom reflectance should be considered when H < 5 m in SYB. The errors caused by NAP and the concentration of NAP are positively and linearly correlated, and the biggest error caused by the bottom occurs at H = 2.5 m.

In SYB, the Chl-a concentration retrieved by SPOT6 data ranges from 0.2 to 2.6 ug/L, and the concentration of NAP ranges from 2.2 to 30.1 mg/L. In general, the IOPs of coastal waters in Sanya are dominated by NAP. The distribution of the Chl-a concentration retrieval results shows that the terrain and the current affect the distribution of the Chl-a concentration. The West Island and the Luhuitou Peninsula block the impact of the southeast winter current and protect the phytoplankton in the coral areas. In general, the Chl-a concentration in shallow waters is higher than that in deep waters, and terrigenous matter may be the main source of nutrients for phytoplankton.

Compared with the remote sensing data from ocean color satellites, the data from land observation satellites have the advantage of higher spatial resolution and more data sources. Though most land observation sensors have the disadvantage of low spectral resolution and lack of NIR bands around 705 nm, the data of land observation satellites have potential for application to oligotrophic or mesotrophic coastal waters. The semi-empirical algorithm based on SPOT6 reflectance data can be used for Chl-a concentration retrieval in optically shallow areas of coastal waters after removing the effects of sun glint, bottom reflectance, and NAP, and it also has the potential to monitor Chl-a in optically shallow waters with a high NAP concentration using other broad-band remote sensing data without the red-edge band (e.g., Landsat8 OLI, ASTER, and GaoFen2).