Inland Water Atmospheric Correction Based on Turbidity Classification Using OLCI and SLSTR Synergistic Observations

Atmospheric correction is an essential prerequisite for obtaining accurate inland water color information. An inland water atmospheric correction algorithm, ACbTC (Atmospheric Correction based on Turbidity Classification), was proposed in this study by using OLCI (Ocean and Land Color Instrument) and SLSTR (Sea and Land Surface Temperature Radiometer) synergistic observations for the first time. This method includes two main steps: (1) water turbidity classification by the GRA index (GRAdient of the spectrum index); and (2) atmospheric correction by synergistic use of OLCI and SLSTR images. The algorithm was validated with 72 in situ sampling sites in Lake Erhai, Lake Hongze, and Lake Taihu, and compared with other atmospheric correction methods, i.e., C2RCC (Case 2 Regional Coast Colour processor), MUMM (Management Unit of the North Seas Mathematical Models), FLAASH (Fast Line-of-sight Atmospheric Analysis of Hypercubes), POLYMER (POLYnomial based algorithm applied to MERIS), and BPAC (Bright Pixel Atmospheric Correction). The results show that (1) the GRA index performed better than the proposed turbidity classification indices, i.e., the Diff (spectral difference index) and the Tind (turbid index), in inland lakes by using the reflectance peak at 1020 nm in clean water; (2) the synergistic use of OLCI and SLSTR performed feasibly for atmospheric correction, and the ACbTC algorithm achieved full-band average values of the mean absolute percentage error (MAPE) = 29.55%, mean relative percentage error (MRPE) = 13.98%, and the root mean square of error (RMSE) = 0.0039 sr−1, which were more reliable than C2RCC, MUMM, FLAASH, POLYMER, and BPAC; and (3) the synergistic use of the 17th band (865 nm) on OLCI and the 5th band (1613 nm) on SLSTR are suitable for clean inland lakes, while both the 5th band (1613 nm) and 6th band (2250 nm) on SLSTR are advisable for the turbidity.


Introduction
The monitoring of inland lakes is becoming gradually more critical due to the increase of anthropogenic pressures and the effects of climate change [1,2]. Compared to in situ measurements, Figure 1. The 12 lakes in the plateau and the plains are annotated using red crosses in (a) Lake Qinghai, Lake Erhai, Lake Chenghai, Lake Dianchi, Lake Yangzonghai, Lake Fuxian, and Lake Xingyun in the plateaus; Lake Poyang, Lake Zhuhu, Lake Hongze, Lake Gehu, and Lake Taihu in the plains. Moreover, three sets of satellite-synchronous in situ stations are presented using red dots, including (b) Lake Erhai (N = 14); (c) Lake Hongze (N = 23), where CZH (Chengzihu), LHB (Lihe Bay) and EHZ (Eastern part of Lake Hongze) are the three regions of Lake Hongze; (d) Lake Taihu (N = 35), where ZSB (Zhushan Bay), MLB (Meiliang Bay), GHB (Gonghu Bay), XKB (Xukou Bay), and ETH (the eastern part of Lake Taihu) are the five regions of Lake Taihu. Dots in black are stations out of the 6 h time window referring to the image time. The 12 lakes in the plateau and the plains are annotated using red crosses in (a) Lake Qinghai, Lake Erhai, Lake Chenghai, Lake Dianchi, Lake Yangzonghai, Lake Fuxian, and Lake Xingyun in the plateaus; Lake Poyang, Lake Zhuhu, Lake Hongze, Lake Gehu, and Lake Taihu in the plains. Moreover, three sets of satellite-synchronous in situ stations are presented using red dots, including (b) Lake Erhai (n = 14); (c) Lake Hongze (n = 23), where CZH (Chengzihu), LHB (Lihe Bay) and EHZ (Eastern part of Lake Hongze) are the three regions of Lake Hongze; (d) Lake Taihu (n = 35), where ZSB (Zhushan Bay), MLB (Meiliang Bay), GHB (Gonghu Bay), XKB (Xukou Bay), and ETH (the eastern part of Lake Taihu) are the five regions of Lake Taihu. Dots in black are stations out of the 6 h time window referring to the image time.

Image and Preprocessing
Nine full-resolution (300 m) OLCI images (OLCI Level 1B, the Instrument Processing Facility (IPF) version is 6.07) and nine corresponding full-resolution (500 m) SLSTR images (SLSTR Level 1B, the IPF version is 6.14) were acquired from the ESA (European Space Agency) data hub (https: //scihub.copernicus.eu/) (see Table 2 for details). The image pre-processing, including image reading, sub-setting, collocation, re-projecting, and radiance to reflectance transferring were performed using SNAP 6.0 (Sentinel Application Platform, http://step.esa.int/main/toolboxes/snap/), which is the toolbox designed to process and analyze the Sentinel series images. Specifically, the Sentinel-3 products were read with per-pixel geo-codings instead of using tie-points, and their vicarious calibration coefficients were set as the default (1.00 for each band) in this study. To match with available OLCI and SLSTR observations, all data were resampled to a 300 m pixel size and were converted to Universal Transverse Mercator (UTM) projection. Then, the 18 top-of-atmosphere (TOA) reflectance bands were extracted from the OLCI images of which 761. 25,764.375, 767.5, 900, and 940 nm were excluded to avoid atmospheric absorption [45]. The geolocation information including SZA (solar zenith angle), SAA (solar azimuth angle), VZA (view zenith angle), VAA (view azimuth angle), and total column water vapor and ozone were obtained from the SLSTR images. Table 2. List of the image data used in this study. The satellite-synchronous in situ sampling measurements were observed in Lake Taihu, Lake Hongze, and Lake Erhai are in bold font. However, it is necessary to discuss the difference of the two sensors for the radiation obtaining ability and the image registration because they are the critical prerequisite for the synergistic use of OLCI and SLSTR [33]. The two sensors have three similar band pairs in the VISR to NIR: OL06-560 nm and SL01-555 nm, OL08-665 nm and SL02-659 nm, and OL17-865 nm and SL03-865 nm (Figure 2a). Afterwards, the central regions of the lakes were selected from Lake Erhai, Lake Honge, and Lake Taihu. Their corresponding TOA reflectance was extracted from these regions (a total of 70,901 pixels) for scatter plot analysis. The results are shown in Figure 2b-d. Despite the small difference in the spectral respond function between OLCI and SLSTR (Figure 2a), three groups of points are concentrated near the 1:1 line, indicating the synergistic use of OLCI and SLSTR is technically feasible. nm and SL01-555 nm, OL08-665 nm and SL02-659 nm, and OL17-865 nm and SL03-865 nm ( Figure  2a). Afterwards, the central regions of the lakes were selected from Lake Erhai, Lake Honge, and Lake Taihu. Their corresponding TOA reflectance was extracted from these regions (a total of 70,901 pixels) for scatter plot analysis. The results are shown in Figure 2b-d. Despite the small difference in the spectral respond function between OLCI and SLSTR (Figure 2a), three groups of points are concentrated near the 1:1 line, indicating the synergistic use of OLCI and SLSTR is technically feasible.
The misregistration of the OLCI and SLSTR images is reduced by the collocation tool from 5~7 pixels to 1~2 pixels based on the visual interpretation [29]. Moreover, the scatter plot results prove that the collocation has inhibited the misregistration between the images. The dark pixel assumption atmospheric correction is based on the average of regional aerosol signals so that the atmospheric correction of OLCI is less affected by the mismatch [18,46]. A few vertical stripes in the scatter plots were caused by resampling of SLSTR which is similar to that of previous studies [47], but the effect is not significant due to its small quantity. Finally, the two registered SWIR-band images of SLSTR (i.e., SL05-1613 nm and SL06-2250 nm) are combined with OLCI into a single image for the synergistic operation. Figure 2. (a) Spectral response function of OLCI and SLSTR at OL06-560 nm, OL08-665 nm, OL17-865 nm, SL01-555 nm, and SL02-659 nm and SL03-865 nm. The solid and dotted lines represent bands of OLCI and SLSTR, respectively. (b-d) Pixel-by-pixel scatter plots between the TOA reflectance of OL06-560 nm and SL01-555 nm, OL08-665 nm and SL02-659 nm, and OL17-865 nm and SL03-865 nm, respectively. The solid and dotted lines represent the 1:1 lines and the regression lines, respectively, and the color bars represent the density of each dot. (a) Spectral response function of OLCI and SLSTR at OL06-560 nm, OL08-665 nm, OL17-865 nm, SL01-555 nm, and SL02-659 nm and SL03-865 nm. The solid and dotted lines represent bands of OLCI and SLSTR, respectively. (b-d) Pixel-by-pixel scatter plots between the TOA reflectance of OL06-560 nm and SL01-555 nm, OL08-665 nm and SL02-659 nm, and OL17-865 nm and SL03-865 nm, respectively. The solid and dotted lines represent the 1:1 lines and the regression lines, respectively, and the color bars represent the density of each dot.
The misregistration of the OLCI and SLSTR images is reduced by the collocation tool from 5~7 pixels to 1~2 pixels based on the visual interpretation [29]. Moreover, the scatter plot results prove that the collocation has inhibited the misregistration between the images. The dark pixel assumption atmospheric correction is based on the average of regional aerosol signals so that the atmospheric correction of OLCI is less affected by the mismatch [18,46]. A few vertical stripes in the scatter plots were caused by resampling of SLSTR which is similar to that of previous studies [47], but the effect is not significant due to its small quantity. Finally, the two registered SWIR-band images of SLSTR (i.e., SL05-1613 nm and SL06-2250 nm) are combined with OLCI into a single image for the synergistic operation.

In Situ Measurement Data
Four field cruises (a total of 72 measurement sites shown in Figure 1) were conducted on the dates shown in Table 2 over Lake Taihu (n = 35), Lake Hongze (n = 23), and Lake Erhai (n = 14), Remote Sens. 2018, 10, 1002 6 of 29 respectively, to examine the atmospheric correction algorithms in this study. Each cruise measured the remote sensing reflectance, R rs (λ), and collected water samples in the same date of satellite images obtaining (Table 2). To obtain high-quality in-situ data at each sample site, sample collection and optical radiometric measurements were performed under clear skies, with no wind or low winds, between 9:00 and 16:00 local time. The longitude and latitude for each site were obtained using a global positioning system receiver, and surface water samples were collected at a depth of 50 cm using Niskin bottles that were frozen at −20 • C for laboratory analysis of the Chla and TSM concentrations [27,48].
The R rs was measured by the FieldSpec spectroradiometer (Analytical Spectral Devices, Inc., Boulder, CO, USA) with a wavelength range from 350-1075 nm. The water spectra were measured according to the above-water method [49], which has been proven to be reliable in previous studies [6,27,36,50,51]. The radiance of the water surface, sky, and a 30% gray board were measured over the deck of an anchored ship. The remote sensing reflectance was extracted by the following equation: where L sw , L sky , L p , and ρ p represent the total radiance, the sky diffuse radiance and the radiance and reflectance of the gray board, respectively; and r demotes skylight reflectance at the air-water surface, and its value is affected by water-surface roughness caused by wind speed (2.2% for calm weather, 2.5% for <5 m/s wind speed, and 2.6-2.8% for 10 m/s wind speed).

The Basic Theory of the Algorithm
The total signal measured from the satellite sensor, ρ t (λ), can be decomposed into several terms as described in Equation (2): where ρ r (λ), ρ a (λ), ρ ar (λ), ρ g (λ), ρ wc (λ), and ρ b (λ) represent the contributions of Rayleigh scattering, aerosol scattering, the coupling term between air molecules and aerosols, sun glint off the water's surface, water whitecaps, and bottom material, respectively. ρ w (λ) is the water leaving reflectance and that is the desired quantity in water remote sensing. T and t represent the direct and diffuse transmittances of the atmospheric column, respectively [11]. The coupling term ρ ar (λ) is zero in the single-scattering case, in which photons are only scattered once, and can be ignored as the amount of multiple scattering is small [46,52]. The sun-glitter-contaminated observations are generally avoided by the observational geometry [15,53]. Moreover, since the light can hardly transmit to the bottom of the lake in turbid inland water, the contribution of the bottom can also be ignored. The whitecap can be negligible due to the inland lake wind speed (the dominant factor of ρ wc (λ) is much smaller). Further, the OLCI could mitigate the sun-glint contamination by tilting cameras in a westerly direction. After ignoring the above contributions, ρ t (λ) in Equation (2) can be rewritten as below: The remaining signal is only contributed by the aerosol scattering and the water-leaving reflectance after the Rayleigh correction which can be performed by Seadas-LUTs (Look-Up Tables of  SeaWiFS Data Analysis System) the aerosols at two wavelengths have a proportional relationship by a large number of aerosol model simulation results [12], and the atmospheric correction factor, ε, can be defined in the following method: where ρ a (λ i ) and ρ a (λ j ) are, respectively, the aerosol reflectances at short wavelength, λ i , and long wavelength, λ j , which satisfies the dark pixel assumption, usually in the NIR or SWIR band. C is the constant under a specific aerosol model, assuming that the aerosol model is relatively stable at the inland lake scale [53].
Since the water-leaving reflectance is almost zero at the dark pixel, the Rayleigh corrected (RC) reflectance ρ t (λ) − ρ r (λ) is similar to ρ a (λ) and Equation (4) can be converted to: After calculating the median value of the C in all the dark pixels, as it is less sensitive to outliers than a regression analysis or average values [19], the aerosol reflectance, ρ a (λ), in the visible channels can be extrapolated as: To calculate an accurate parameter C, it is essential to select the suitable correction bands and to choose the dark pixels. Usually, the dark pixels are considered as the water pixels with low reflectance due to the characteristics of pure water [19]. In general, more turbid waters need the longer correction bands [54]. The parameter C will become greater when shorter bands are selected in turbid waters where the water signal could not be negligible. Furthermore, the ρ a calculated by Equation (6) will become greater at short wavelengths resulting in overcorrection for aerosols [55]. However, if the selected bands are too long, the estimated aerosol reflectance might be insufficient. Therefore, it is inappropriate to fix the specific correction bands in Equation (5) for lakes under different turbid conditions. The calculation of parameter C in this study will be completed by respective band pairs based on the proposed water turbidity detection method.
To assess the performance of the proposed ACbTC algorithm, the other atmospheric correction algorithms are compared in this study, which are FLAASH [27,56], MUMM [10], C2RCC [57], POLYMER [7], and BPAC [14]. Due to the similar band settings with MERIS, the MUMM_alpha and MUMM_gamma coefficients were set as 1.817 and 1.0, which is recommended by [17,58]. Meanwhile, the MUMM_epsilon was calculated from the Rayleigh-corrected reflectance over the region of clean pixels for each image [10,18]. Please see more details in the references [10,17,18,58]. As the standard OLCI algorithm, the BPAC results could be downloaded from the EUMETSAT (European Organisation for the Exploitation of Meteorological Satellites, https://codarep.eumetsat.int). The system vicarious calibration (SVC) coefficients used in POLYMER (https://www.hygeos.com/polymer) were set as 1.0 in keeping with other algorithms [59,60].

The Processing Procedure of the ACbTC
Processing procedure of the ACbTC algorithm in this study is mainly divided into three steps ( Figure 3): First, the image was pre-processed according to Section 2.2.1; second, the Rayleigh correction was conducted by Seadas LUTs and 6SV for OLCI and SLSTR, respectively; and, lastly, the aerosol reflectance was calculated based on turbidity detection by the GRA index using the dark pixel method [11] according to Section 3.1.

Rayleigh Correction for the OLCI and SLSTR
The Rayleigh correction for OLCI is performed by Seadas LUTs to facilitate the atmospheric correction in Seadas (i.e., MUMM). However, the latest version of the Seadas does not support the Rayleigh correction for SLSTR. Thus, the bands of SL05-1613 nm and SL06-2250 nm were Rayleigh corrected using the 6SV model. To reduce the number of calculations, several points were selected to calculate the Rayleigh scattering, which represents the whole area: (1) The smallest rectangle containing the lake area was divided into several squares of the same size (the number of squares depends on the computing ability and the number in this study was selected as nine), and the center pixel of each square was set as the representative point (shown in the Figure 4a). (2) The corresponding auxiliary information of each point were extracted from the SLSTR dataset.
(3) The 6SV model was compiled for each point to simulate the scattering process of the Rayleigh molecules in advance. Then, the Rayleigh reflectance spectrum was constructed and it was assigned to represent this square area. (4) In each square, the reflectance of the top of the atmosphere was subtracted from the Rayleigh reflectance to complete the Rayleigh scattering correction.
To analyze the difference between the two Rayleigh correction methods, the Rayleigh correction was performed by using both methods for OLCI bands. The corresponding Rayleigh correction results of 72 in situ points were selected for linear analysis (Figure 4b,c). The results show that the linear correlation coefficient becomes closer to 1 with the increasing wavelength and the results from Seadas and 6SV are similar at longer wavelengths. Due to the weak Rayleigh scatter effect at long wavelengths [45], the 6SV can be used to supplement Seadas at SL05-1613 nm and SL06-2250 nm.

Aerosol Reflectance Calculation by Turbidity Detection
To calculate the aerosol reflectance, the dark pixel should be selected, and the parameter C should be estimated depending on the zero water contribution.

Rayleigh Correction for the OLCI and SLSTR
The Rayleigh correction for OLCI is performed by Seadas LUTs to facilitate the atmospheric correction in Seadas (i.e., MUMM). However, the latest version of the Seadas does not support the Rayleigh correction for SLSTR. Thus, the bands of SL05-1613 nm and SL06-2250 nm were Rayleigh corrected using the 6SV model. To reduce the number of calculations, several points were selected to calculate the Rayleigh scattering, which represents the whole area: (1) The smallest rectangle containing the lake area was divided into several squares of the same size (the number of squares depends on the computing ability and the number in this study was selected as nine), and the center pixel of each square was set as the representative point (shown in the Figure 4a). (2) The corresponding auxiliary information of each point were extracted from the SLSTR dataset.
(3) The 6SV model was compiled for each point to simulate the scattering process of the Rayleigh molecules in advance. Then, the Rayleigh reflectance spectrum was constructed and it was assigned to represent this square area. (4) In each square, the reflectance of the top of the atmosphere was subtracted from the Rayleigh reflectance to complete the Rayleigh scattering correction.

Rayleigh Correction for the OLCI and SLSTR
The Rayleigh correction for OLCI is performed by Seadas LUTs to facilitate the atmospheric correction in Seadas (i.e., MUMM). However, the latest version of the Seadas does not support the Rayleigh correction for SLSTR. Thus, the bands of SL05-1613 nm and SL06-2250 nm were Rayleigh corrected using the 6SV model. To reduce the number of calculations, several points were selected to calculate the Rayleigh scattering, which represents the whole area: (1) The smallest rectangle containing the lake area was divided into several squares of the same size (the number of squares depends on the computing ability and the number in this study was selected as nine), and the center pixel of each square was set as the representative point (shown in the Figure 4a). (2) The corresponding auxiliary information of each point were extracted from the SLSTR dataset.
(3) The 6SV model was compiled for each point to simulate the scattering process of the Rayleigh molecules in advance. Then, the Rayleigh reflectance spectrum was constructed and it was assigned to represent this square area. (4) In each square, the reflectance of the top of the atmosphere was subtracted from the Rayleigh reflectance to complete the Rayleigh scattering correction.
To analyze the difference between the two Rayleigh correction methods, the Rayleigh correction was performed by using both methods for OLCI bands. The corresponding Rayleigh correction results of 72 in situ points were selected for linear analysis (Figure 4b,c). The results show that the linear correlation coefficient becomes closer to 1 with the increasing wavelength and the results from Seadas and 6SV are similar at longer wavelengths. Due to the weak Rayleigh scatter effect at long wavelengths [45], the 6SV can be used to supplement Seadas at SL05-1613 nm and SL06-2250 nm.

Aerosol Reflectance Calculation by Turbidity Detection
To calculate the aerosol reflectance, the dark pixel should be selected, and the parameter C should be estimated depending on the zero water contribution. To analyze the difference between the two Rayleigh correction methods, the Rayleigh correction was performed by using both methods for OLCI bands. The corresponding Rayleigh correction results of 72 in situ points were selected for linear analysis (Figure 4b,c). The results show that the linear correlation coefficient becomes closer to 1 with the increasing wavelength and the results from Seadas and 6SV are similar at longer wavelengths. Due to the weak Rayleigh scatter effect at long wavelengths [45], the 6SV can be used to supplement Seadas at SL05-1613 nm and SL06-2250 nm.

Aerosol Reflectance Calculation by Turbidity Detection
To calculate the aerosol reflectance, the dark pixel should be selected, and the parameter C should be estimated depending on the zero water contribution.
(1) Dark pixel selection Some studies determine the dark pixels directly by defining thresholds [19]. However, this method is limited by spatial objects and is not suitable for groups of lakes. In this study, histogram statistics of water pixels are used to ensure that the darkest (cleanest) pixels can be found in each study area. The specific three steps are as follows: (a) Exclude the boundary between water and land using the Canny filter [61] to avoid the mixed pixel contamination; (b) Undeterred by the complex and diverse geographical environments of the images, more stringent taxonomies are established based on NDWI [62] to separate the water pixel from the ice, snow and clouds. The thresholds are determined by the strict visual interpretation in this study [29]: (c) The intersection of the two 10th percentiles of OL17-865 nm and SL05-1613 nm [8,18,19] is determined as the threshold for dark pixels selection.
(2) Turbidity detection by the GRA index The band pairs need to be determined when Equation (6) is used to calculate the aerosol reflectance. For this purpose, the GRA index was proposed to detect the water turbidity and, then, the different band pairs were fixed for the clean or turbid waters.
A large amount of Rayleigh corrected images show that there is a clear morphological spectrum difference between clean and turbid water in the NIR-SWIR range. Thus, the GRA index was promoted to identify the turbid water by calculating the difference between the gradient of OL18-885 nm to OL21-1020 nm and that of OL18-885 nm to SL05-1613 nm. The GRA can be calculated by Equation (8) as follows: where R rs (885), R rs (1020), and R rs (1613) are the Rayleigh corrected remote sensing reflectance of the OL18-885 nm and OL21-1020 nm and SL05-1613 nm, respectively. The numbers 885, 1020, and 1613 are the center wavelengths (nm) of the above corresponding bands.
To determine the threshold of GRA for turbidity classification, statistical analysis of the histograms were performed by visually interpreting 5000 pixels of both turbid and clean lakes. It is reasonable to set −0.07 as the threshold of the GRA to efficiently distinguish clean or turbid water in this study ( Figure 5). The overlap of columns for Turbid and Clean with the limited quantity, which accounts for~0.98%. Additionally, the histogram statistics of all pixels of the study lakes also have a clear distinction at −0.07. Finally, we temporarily define −0.07 as the threshold of the GRA. Thus, a pixel with a GRA < −0.07 is considered as turbid water. Otherwise, it belongs to the clean water category.
(3) Turbidity detection by the GRA index For turbid waters, the assumption of water reflectance approaching zero is invalid when the wavelength is longer than 1000 nm, so the calculation of aerosol reflectance is required for SL05-1613 nm and SL06-2250 nm, while, for clean waters, the combination is OL17-865 nm and SL05-1613 nm since the assumption is valid in these wavelengths (see the related procedure in Equations (4)- (6)).
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 29 Figure 5. Histogram distribution of the GRA index of water pixels where the label "All pixels" represents all the study lakes in this study, and the label "Clean" and "Turbid" are, respectively, the clean and turbid water pixels that we manually interpreted referenced from water quality parameters and spectral curves.

Accuracy Assessment
The statistical analysis considered the comparison between the in-situ measured Rrs vs. OLCIderived Rrs. Six popular indices, namely, the MAPE, MRPE, RMSE, slope, intercept, and Pearson correlation coefficient, were used for assessing the accuracy of atmospheric correction. The MAPE, MRPE, and RMSE are defined as Equations (9)-(11), respectively:

MAPE=
where xi and yi are the in situ measured and OLCI-derived reflectance for the ith sample site and N is the number of samples.

Visual Inspection of Atmospherically Corrected Images
The original images of false color and the Rrs of the selected bands corrected by ACbTC of Lake Erhai, Lake Hongze and Lake Taihu are shown in Figure 6. The presented six bands are common water color bands, of which 442.5 nm is the characteristic band of Colored Dissolved Organic Matter (CDOM) where the absorption of CDOM and Non-Algal Particulate (NAP) predominates [63]; 560 nm and 708.75 nm are, respectively, the first and second reflection peak positions of Chla; 620 nm is the position of the absorption peak of phycocyanin; and 778.75 nm is commonly used for the estimation of TSM in turbid inland lakes [27]. These six bands can describe the Rrs of the entire OLCI bands and provide a qualitative analysis for atmospheric correction. The area of cyanobacteria . Histogram distribution of the GRA index of water pixels where the label "All pixels" represents all the study lakes in this study, and the label "Clean" and "Turbid" are, respectively, the clean and turbid water pixels that we manually interpreted referenced from water quality parameters and spectral curves.

Accuracy Assessment
The statistical analysis considered the comparison between the in-situ measured R rs vs. OLCI-derived R rs . Six popular indices, namely, the MAPE, MRPE, RMSE, slope, intercept, and Pearson correlation coefficient, were used for assessing the accuracy of atmospheric correction. The MAPE, MRPE, and RMSE are defined as Equations (9)-(11), respectively: where x i and y i are the in situ measured and OLCI-derived reflectance for the ith sample site and N is the number of samples.

Visual Inspection of Atmospherically Corrected Images
The original images of false color and the R rs of the selected bands corrected by ACbTC of Lake Erhai, Lake Hongze and Lake Taihu are shown in Figure 6. The presented six bands are common water color bands, of which 442.5 nm is the characteristic band of Colored Dissolved Organic Matter (CDOM) where the absorption of CDOM and Non-Algal Particulate (NAP) predominates [63]; 560 nm and 708.75 nm are, respectively, the first and second reflection peak positions of Chla; 620 nm is the position of the absorption peak of phycocyanin; and 778.75 nm is commonly used for the estimation of TSM in turbid inland lakes [27]. These six bands can describe the R rs of the entire OLCI bands and provide a qualitative analysis for atmospheric correction. The area of cyanobacteria blooms in Lake Taihu (Figure 6d0), highlighted in the yellow rectangle, was masked to improve the interpretability of the R rs images from the western part to ZSB. blooms in Lake Taihu (Figure 6d0), highlighted in the yellow rectangle, was masked to improve the interpretability of the Rrs images from the western part to ZSB.
From the studies [26,29,64] and our in situ measurement data (Table 1), Lake Erhai is the cleanest among the three lakes but Lake Hongze and Lake Taihu are all dominated by suspended matter. Additionally, Lake Taihu is dominated both by suspended matter and algae which is more eutrophic in summer than Lake Hongze [63]. The maximum Rrs values of Lake Hongze and Lake Taihu are, respectively, 681.25 nm and 560 nm, which correspond to the fact that they are dominated by optically-active substances, while their minimum values are 778.75 nm and 681 nm due to the strong absorption of pure water and pigment, respectively. These results mean that the atmosphericallycorrected Rrs at different bands corresponds to the water quality status according to previous studies [26,29,64]. As shown in Figure 6a1-a6, the ACbTC-corrected Rrs values of the six bands in Lake Erhai range from 0.002 to 0.01 sr −1 . The image of Rrs(442.5) presents high values (~0.004 sr −1 ) in the southern area. However, the images of the other bands present high values in the northern region, of which the highest is Rrs(560) with the broadest coverage. The high Rrs occurs because the northern part of Lake Erhai is the area where the algal bloom is most likely to occur [26,64] and Rrs(560) is sensitive to Chla concentration. With the increasing wavelengths, the Rrs of Lake Erhai tends to be spatially homogenous (~0.002 sr −1 ) over most of the area in 708.75 and 778.75 nm. The magnitude and spatial distribution of Rrs in Lake Erhai after ACbTC are the same as those of the previous studies [64].
The images of Lake Hongze obtained in 7 December 2016 and 18 May 2018 are shown in Figure  6b0-b6,c0-c6, respectively. Lake Hongze (Figure 6b0,c0) is the most turbid of the three lakes, but it is visually cleaner in west LHB, which is mainly covered by wetlands [42]. After the atmospheric correction, the images of the two seasons both show strong spatial heterogeneity at the six bands. The west LHB has the lowest full-band value, close to 0.012 sr −1 , while the center is the most turbid at 681.25 nm, approaching 0.05 sr −1 . Among them, the spatial distribution heterogeneity of Rrs(778.75) is the strongest, mainly due to the dominant absorption of pure water in these wavelengths, the weak are presented at (a0-d0) using the same stretching method in the study areas. The atmospheric correction results by the ACbTC algorithm, i.e., Lake Erhai-20170419, Lake Hongze in winter (Lake HongzeW)-20161207, Lake Hongze in spring (Lake HongzeS)-20170518, and Lake Taihu-20170724, for which (a1-d1) R rs (442.5), (a2-d2) R rs (560), (a3-d3) R rs (620), (a4-d4) R rs (681.25), (a5-d5) R rs (708.75), and (a6-d6) R rs (778.75) in units of sr −1 in the corresponding lakes. The highlighted yellow rectangle is masked to improve the interpretability of the R rs images from the western part to ZSB.
From the studies [26,29,64] and our in situ measurement data (Table 1), Lake Erhai is the cleanest among the three lakes but Lake Hongze and Lake Taihu are all dominated by suspended matter. Additionally, Lake Taihu is dominated both by suspended matter and algae which is more eutrophic in summer than Lake Hongze [63]. The maximum R rs values of Lake Hongze and Lake Taihu are, respectively, 681.25 nm and 560 nm, which correspond to the fact that they are dominated by optically-active substances, while their minimum values are 778.75 nm and 681 nm due to the strong absorption of pure water and pigment, respectively. These results mean that the atmospherically-corrected R rs at different bands corresponds to the water quality status according to previous studies [26,29,64].
As shown in Figure 6a1-a6, the ACbTC-corrected R rs values of the six bands in Lake Erhai range from 0.002 to 0.01 sr −1 . The image of R rs (442.5) presents high values (~0.004 sr −1 ) in the southern area. However, the images of the other bands present high values in the northern region, of which the highest is R rs (560) with the broadest coverage. The high R rs occurs because the northern part of Lake Erhai is the area where the algal bloom is most likely to occur [26,64] and R rs (560) is sensitive to Chla concentration. With the increasing wavelengths, the R rs of Lake Erhai tends to be spatially homogenous (~0.002 sr −1 ) over most of the area in 708.75 and 778.75 nm. The magnitude and spatial distribution of R rs in Lake Erhai after ACbTC are the same as those of the previous studies [64].
The images of Lake Hongze obtained in 7 December 2016 and 18 May 2018 are shown in Figure 6b0-b6,c0-c6, respectively. Lake Hongze (Figure 6b0,c0) is the most turbid of the three lakes, but it is visually cleaner in west LHB, which is mainly covered by wetlands [42]. After the atmospheric correction, the images of the two seasons both show strong spatial heterogeneity at the six bands. The west LHB has the lowest full-band value, close to 0.012 sr −1 , while the center is the most turbid at 681.25 nm, approaching 0.05 sr −1 . Among them, the spatial distribution heterogeneity of R rs (778.75) is the strongest, mainly due to the dominant absorption of pure water in these wavelengths, the weak absorption of particulates, and the strong particulate scatter effect, which can efficiently represent the distribution of the TSM concentration [42].
Among the six bands of the ACbTC-corrected results of Lake Taihu, the R rs (560) is locally the highest, approximately 0.032 sr −1 in MLB, GHB, and XKB, showing that the dominant optically-active substance of these parts of Lake Taihu are algal particles, which is consistent with previous investigations [3,6,21,53]. In the images of R rs (620) and R rs (681.25), the value of XKB is high, which is consistent with the water quality characteristics dominated by suspended matter [65]. The lower R rs (442.5),~0.01 sr −1 , is observed near the area of the algal bloom, which is mainly affected by the strong absorption of pigment particles. Additionally, the high values appearing in 708.75 and 778.75 nm show the characteristic spectrum of algal blooms. In general, Lake Taihu (Figure 6d0) is cleaner in the west to the central area, while more turbid in the eastern part, including XKB and ETH. Additionally, the lower R rs of the six bands in the center of Lake Taihu is consistent with the observation of Wang et al. [3] where the dark pixels were selected in this study.

Comparison with In Situ Measurements
To quantitatively assess the accuracy of the proposed algorithm, the in situ measured and atmospherically-corrected R rs values were compared. Due to the relatively stable condition of inland lakes [66], the time window between in situ measurements and the satellite overpass was relaxed to ±6 h to obtain a sufficient number of coincident matchup pairs [50,67]. To avoid sub-pixel light contamination, the OLCI-derived remote sensing reflectance (R rs ) values were extracted using a 3 × 3 average pixel-box. The coefficient of variation (CV) of each pixel-box was less than 3% to ensure their homogeneity [51]. Additionally, the sampling sites within two pixels from the shore were excluded to avoid land-adjacent effects based on prior knowledge [29]. To match with the OLCI data, the in situ spectra were convolved to represent the spectra of OLCI using a method referred from the literature [27]. Figure 7 presents the scatter plots of the in situ measured and OLCI-derived R rs values. Overall, the scattered dots are concentrated near the 1:1 line. Notably, the scattering dots of 400 and 412.5 nm show a flat distribution that means the ACbTC algorithm fails to distinguish the different water conditions. For these two bands, the OLCI-derived and ASD-measured CV are, respectively, approximately 4.5 and 26% (Lake Taihu), 8.0 and 22% (Lake Hongze), and 28 and 40% (Lake Erhai). We attribute this to the following two aspects: (1) the spatial heterogeneity of the Rayleigh and aerosol scattering mainly occurs in the shortwave bands, while the assumption of our algorithm is incapable of that [45]; and (2) the aerosol correction parameter, C, in the dark pixel method is calculated by fitting the reference bands through the exponential function, so the error is relatively large when extrapolated to the short wavelengths, as shown in Equation (6), leading to inaccurate aerosol estimation [11]. However, scattering dots begin to show a more uniform distribution from 490 nm onward, indicating that the OLCI-derived results begin to reflect the difference between sampling measurements. Additionally, the overcorrection occurs between 500 and 700 nm when the measured R rs is greater than 0.035 sr −1 , especially in the more turbid areas (e.g., Lake Hongze). The overcorrection may be due to the non-aerosol signal remaining in the reference bands, resulting in the overestimation of the aerosol reflectance. However, these overcorrections are within an acceptable range, and there are similar correction results in these bands in previous studies based on the dark pixel method [18,55,68].
The relationship of the ACbTC results between the total data and each lake is shown in Figure 8. Note that only MAPE, MRPE, and RMSE values were analyzed in this part due to the apparent optical deviation between the study lakes ( held steady from 510 to 708.75 nm, increase gradually from 708.75 nm, and reach the maximum at 885 nm. There are four MAPE local maximums at 400, 560, 753.75, and 885 nm which are mainly caused by Lake Erhai and Lake Taihu, as shown in Figure 8a. However, Lake Taihu, Lake HongzeW, and Lake HongzeS maintain the low MAPE values from 400 to 708.75 nm and Lake Taihu performs worse than Lake Hongze. The MRPE values of both Lake Erhai and Lake Taihu show very high, from 400 to 442.5 nm, and the local minimum at 560 nm, which indicates the undercorrection and overcorrection by ACbTC, respectively. The MRPE values of Lake Hongze are better than the other lakes. The total RMSE value at 400-442.5 nm is contributed from Lake Erhai and Lake Taihu. The two causes of this high error in the NIR of Lake Erhai and Lake Taihu are as follows: (1) the NIR water signal in Lake Erhai is weak and susceptible to land-adjacent effects, causing the algorithm to ineffectively eliminate the signal of aerosols [31]; and (2) dominated by suspended matter and algae, Lake Taihu presents considerable spatial heterogeneity in the optical properties in the NIR. Furthermore, Lake Taihu may not satisfy the assumption of the dark pixel method of spatial average for its relatively large area [69]. These two factors may contribute to the poor atmospheric correction performance of Lake Taihu in the NIR. high, from 400 to 442.5 nm, and the local minimum at 560 nm, which indicates the undercorrection and overcorrection by ACbTC, respectively. The MRPE values of Lake Hongze are better than the other lakes. The total RMSE value at 400-442.5 nm is contributed from Lake Erhai and Lake Taihu.
The two causes of this high error in the NIR of Lake Erhai and Lake Taihu are as follows: (1) the NIR water signal in Lake Erhai is weak and susceptible to land-adjacent effects, causing the algorithm to ineffectively eliminate the signal of aerosols [31]; and (2) dominated by suspended matter and algae, Lake Taihu presents considerable spatial heterogeneity in the optical properties in the NIR. Furthermore, Lake Taihu may not satisfy the assumption of the dark pixel method of spatial average for its relatively large area [69]. These two factors may contribute to the poor atmospheric correction performance of Lake Taihu in the NIR.

Comparison with Other Atmospheric Algorithms
To assess the performance of the proposed ACbTC algorithm, the atmospheric correction results of the FLAASH, MUMM, C2RCC, POLYMER, and BPAC were compared in this section from the spectral shape and the magnitude of R rs . Figure 9 shows a stricter comparison for analysis of the spectral shape between the in situ R rs values (n = 9) and different atmospheric correction results for a match-up time window of less than 1 h [18,70,71]. Figure 10 presents the comparison of ACbTC with other atmospheric algorithms in MAPE, MRPE, RMSE, slope, r, and intercept (n = 72).

Comparison with Other Atmospheric Algorithms
To assess the performance of the proposed ACbTC algorithm, the atmospheric correction results of the FLAASH, MUMM, C2RCC, POLYMER, and BPAC were compared in this section from the spectral shape and the magnitude of Rrs. Figure 9 shows a stricter comparison for analysis of the spectral shape between the in situ Rrs values (N = 9) and different atmospheric correction results for a match-up time window of less than 1 h [18,70,71]. Figure 10 presents the comparison of ACbTC with other atmospheric algorithms in MAPE, MRPE, RMSE, slope, r, and intercept (N = 72). As shown in the Figure 9a-c, the Rrs values derived by FLAASH are far lower than those of the in situ Rrs values from 400 to 708.75 nm, indicating that FLAASH overcorrects the Rrs in Lake Erhai. The same situation occurred in Lake Hongze (Figure 9d-f) and Lake Taihu (Figure 9g-i). However, FLAASH works well in NIR bands in the sediment-dominated Lake Hongze validating that its application for TSM/SSD estimation bands (nearing NIR) in turbid inland lakes is suitable [27,72]. This improvement leads to a gradual decrease of MAPE and MRPE of FLAASH after 708.75 nm four strict synchronization samples match very well with the in situ Rrs values. Although there is a relative bias in Figure 9c,e, the results of ACbTC are still optimal compared to those of the other algorithms. For three synchronization samples of Lake Taihu, ACbTC has an insufficient correction at short wavelengths as well, but it is better than that of MUMM because its better performance of the spectral shape in NIR. Meanwhile, the MAPE, MRPE, and RMSE of ACbTC are similar with MUMM, but the regression results of the former is better than the latter, as shown in Figure 10d-f. As shown in the Figure 9a-c, the R rs values derived by FLAASH are far lower than those of the in situ R rs values from 400 to 708.75 nm, indicating that FLAASH overcorrects the R rs in Lake Erhai. The same situation occurred in Lake Hongze (Figure 9d-f) and Lake Taihu (Figure 9g-i). However, FLAASH works well in NIR bands in the sediment-dominated Lake Hongze validating that its application for TSM/SSD estimation bands (nearing NIR) in turbid inland lakes is suitable [27,72]. This improvement leads to a gradual decrease of MAPE and MRPE of FLAASH after 708.75 nm (Figure 10a,b). Although the R rs value trend of C2RCC is similar to that of the in situ R rs , these results present the most severe deviation (the undercorrection in Lake Erhai and the overcorrection in Lake Hongze and Lake Taihu) and the erroneous correction of the characteristic bands, such as 708.75 nm (Figure 9d-f). The MAPE and MRPE results show that the C2RCC is less effective than FLAASH (Figure 10a,b), the overall bias is the highest of all (Figure 10c), and the regression results also present the worst among all of the algorithms (Figure 10d-f). This is caused by the fact that C2RCC is predetermined by the neural network (NN) method, with the training data mostly simulated by the Hydrolight model [73] and collected from European water [57] whose water constituents and bio-optical features were considerably different from those of the lakes evaluated in this study [57]. MUMM performs better than C2RCC and FLAASH algorithms and has a smaller deviation with the in situ R rs , even in Lake Erhai (Figure 9a-c). However, MUMM has an overcorrection in Lake Hongze as well and fails to reduce the atmospheric effect in the short wavelength bands in Lake Taihu. The deviation (MAPE, MRPE, and RMSE) of MUMM maintains stability from 500 to 700 nm and the correction result still performs better than C2RCC and FLAASH (Figure 10a-c) showing the stability of MUMM [18]. The error mainly comes from the short bands (<500 nm) which also results in the poor regression outcomes (Figure 10d,e). The overcorrection of POLYMER is observed in all lakes (Figure 10b) but it is reduced in Lake Erhai from 510 to 620 nm (Figure 9a-c). In addition, the overcorrection in Lake Hongze and Lake Taihu becomes more serious compared to C2RCC (Figure 9d-i). This overcorrection indicates that POLYMER may be suitable for clean inland lakes rather than extremely turbid inland lakes. The MAPE of POLYMER is similar to C2RCC after 620 nm, but improves at short wavelengths (Figure 10a) showing more stability than C2RCC. The regression results of POLYMER performs better than C2RCC and FLAASH. In contrast to POLYMER, the BPAC shows serious overcorrection in Lake Erhai (Figure 9a-c). However, the results in Lake Hongze and Lake Taihu are, respectively, similar to those of MUMM (Figure 9d-f) and POLYMER (Figure 9g-i) in spectral shape. The MAPE of BPAC is greater than 80% in the whole wavelength (Figure 10a) and the MRPE is less than −80% (Figure 10b), showing the overcorrection of BPAC in the whole wavelength. This indicates that the IOP models of BPAC might be not suitable for inland lakes in this study. The ACbTC algorithm performs the best, as shown in Figure 9a,b,d,f, and the spectra of the four strict synchronization samples match very well with the in situ R rs values. Although there is a relative bias in Figure 9c,e, the results of ACbTC are still optimal compared to those of the other algorithms. For three synchronization samples of Lake Taihu, ACbTC has an insufficient correction at short wavelengths as well, but it is better than that of MUMM because its better performance of the spectral shape in NIR. Meanwhile, the MAPE, MRPE, and RMSE of ACbTC are similar with MUMM, but the regression results of the former is better than the latter, as shown in Figure 10d-f.
In summary, FLAASH, C2RCC and BPAC obtained the most inferior results, showing significant differences in the magnitude and spectral shape of the corrected R rs results, especially before the NIR. Among them, FLAASH slightly improved after the NIR, mainly because of its excellent performance for Lake Hongze. POLYMER performed more stably than FLAASH and C2RCC, but there was a ubiquitous overcorrection phenomenon, especially in turbid lakes. Furthermore, both MUMM and ACbTC showed excellent and stable correction effects, but the latter was more accurate for the overall state of the spectrum, showing better statistical results.

The Usability of OL21-1020 nm for Distinguishing Turbid Water from Clean Water
The R rs of OL21-1020 nm was used to detect water turbidity by Equation (8) depending on the weak peak of clean inland water. However, is this band useful for distinguishing between clean and turbid water? By extracting R rs values from nine Sentinel-3 images (Table 2), a weak peak of Rayleigh corrected R rs at the 1020 nm band (denoted as the 1020 nm peak) was observed in clean inland water, such as Lake Erhai, Lake Qinghai, and Lake Zhuhu, as shown in Figure 11. Thus, this peak might be an indicator for the identification of these two types of water. To check this hypothesis, a simple normalization index, T, was constructed and applied to the study area for identification. The formula of the T index is: where R rsRC (885) and R rsRC (1020) denote the Rayleigh corrected R rs in OL18-885 nm and OL21-1020 nm, respectively. When T ≥ 0, the study area belongs to the clean inland waters; otherwise, it is turbid water.
The T values of the surveyed lakes are shown in Figure 11 and are sliced where T ≥ 0. In the turbid lakes, the T values are less than zero in the central region of lakes except in the nearshore region, while, in the clean water, the T values are greater than zero in the central region with a decreasing trend from nearshore region to the lake center. This difference indicates that R rs (1020) is useful in distinguishing the clean water from the central region of turbid water. However, the nearshore region of the turbid water may be confused with the clean water since the T value is greater than zero there as well (as shown the spectrum of Lake Hongze_1 in Figure 11).
Why is the 1020 nm peak present in the clean water and the nearshore region of turbid water? Is this peak caused by the water signal or a non-water signal? To answer these questions, the absorption and backscattering coefficient of pure water was first analyzed. The absorption spectrum of pure water presents the first absorption peak at 972 nm and an apparent absorption valley at 1071 nm [74,75], while the backscattering coefficient is consistent after 800 nm that means the peak may be caused by other components in the water or by the effect of the shore reflectance spectrum. Since there are very low concentrations of Chla and TSM in the clean water (Table 1), the peak is mainly caused by the impact of land stray light [76][77][78][79], which is also called the land adjacent effect [80][81][82]. Therefore, images of the T index can also characterize the impact of land stray light on the water signal. The land stray light has a different effect on the reflectance of 1020 nm in all the images surveyed in this study, mainly due to the geo-location of view, the atmospheric condition, the lake area, land cover type, and ground elevation [80]. The sag ponds, such as Lake Fuxian, Lake Erhai, and Lake Dianchi, are always surrounded by high altitude mountains whose dense vegetation cover and snow cover can contaminate the 1020 nm signal of the water surface [76,77,83]. As shown in Figure 12, the T values of Lake Erhai, relatively small in area, shows a clear spatial distribution which gradually decreases from the nearshore to the central region. This result indicates that the contamination of land stray light at 1020 nm decreases with longer distances from the shore. Although there are no high mountains around the Lake Qinghai, the contamination from the vegetation and rocks still exist, as shown in Figure 12.
The plain lakes, such as Lake Taihu, Lake Hongze, and Lake Poyang, were mainly surrounded by vegetation, wetland, and urban land. However, are suffered from the stray light as well, as shown in Figure 12. For turbid water, previous studies have proved that the backscattering of sediments in highly turbid water is not close to zero even in the vicinity of 1071 nm [30,31,84]. Therefore, the higher reflectance at 885 and 1020 nm will weaken the influence of land stray light when the distance is far from the land. This phenomenon is present in Lake Hongze_2, as shown in Figure 11. However, when near the lakeshore, the peak appears to be due to the stronger stray light caused by the land cover (Lake Hongze_1). in Figure 12. For turbid water, previous studies have proved that the backscattering of sediments in highly turbid water is not close to zero even in the vicinity of 1071 nm [30,31,84]. Therefore, the higher reflectance at 885 and 1020 nm will weaken the influence of land stray light when the distance is far from the land. This phenomenon is present in Lake Hongze_2, as shown in Figure 11. However, when near the lakeshore, the peak appears to be due to the stronger stray light caused by the land cover (Lake Hongze_1). Figure 11. Six typical spectral curves from the different water bodies. (a) Lake Hongze_1 and Lake Hongze_2 represent the pixels near the shore and away from the shore, respectively. The spectral curves of Lake Erhai, Lake Taihu, Lake Qinghai, and Lake Zhuhu are located in the corresponding lakes. To eliminate differences in magnitude between different turbid conditions, all spectral curves were normalized based on Rrs(2250) with their distributions displayed in the histogram (b).  Figure 11. Six typical spectral curves from the different water bodies. (a) Lake Hongze_1 and Lake Hongze_2 represent the pixels near the shore and away from the shore, respectively. The spectral curves of Lake Erhai, Lake Taihu, Lake Qinghai, and Lake Zhuhu are located in the corresponding lakes. To eliminate differences in magnitude between different turbid conditions, all spectral curves were normalized based on R rs (2250) with their distributions displayed in the histogram (b). Figure 11. Six typical spectral curves from the different water bodies. (a) Lake Hongze_1 and Lake Hongze_2 represent the pixels near the shore and away from the shore, respectively. The spectral curves of Lake Erhai, Lake Taihu, Lake Qinghai, and Lake Zhuhu are located in the corresponding lakes. To eliminate differences in magnitude between different turbid conditions, all spectral curves were normalized based on Rrs(2250) with their distributions displayed in the histogram (b).

Wavelength (nm)
The discussion of the T index shows: first, the signal of 1020 nm is contaminated by land stray light no matter in clean or turbid water; second, turbid water with strong backscatter at 1020 nm will inhibit stray light effect, especially when the water pixels are far from the land; and, third, the clear water is vulnerable to form the 1020 nm peak due to the weaker backscattering. Therefore, the 1020 nm peak could be an indicator to discriminate between clean and turbid water excluding pixels near the shore.

The Effect of the GRA Index to Distinguish the Water Turbidity
As described in Section 5.1, the T index could not discriminate between clean inland water and turbid nearshore water. Thus, the improvement is required to distinguish the water turbidity more accurately. In fact, the R rs (885) varies significantly in inland lakes [29,63]. The high R rs (885) values of turbid water are mainly affected by the backscattering of the particles while the relatively low R rs (885) values of clean water are affected primarily by the absorption of the pure water. Meanwhile, the R rs value at 1613 nm is relatively low in the different waters. Four typical water spectra are shown in Figure 13 of which the slope of the turbid nearshore water (spectrum d) is flatter than that of the clean water and the central turbid water from 885 to 1613 nm. Therefore, Equation (8) can efficiently distinguish the turbid water with the combination of the 1020 nm peak detection (GRA 1 ) and the decreasing gradient (GRA 2 ).
To evaluate this method, the results of the GRA were compared with other methods, such as the Diff and Tind indices. The Diff index was proposed to classify the turbidity of near-shore water by calculating the difference between the logarithmic gradients of R rsRC (470) and R rsRC (660) and that of R rsRC (470) and R rsRC (1240) [85,86]. It is assumed that Diff is less than 0 for clean water pixels. The calculation bands of Diff were modified slightly according to the OLCI bands in this study, which are 490, 665, and 1020 nm. The Tind index uses the combinations of MODIS measured radiance or reflectance in the short visible, NIR and SWIR bands to detect the water turbidity area [8,23]. The corresponding bands of OLCI + SLSTR for Tind are 753.75, 1020, and 2250 nm. The pixels of Tind > 1.3 are considered as the turbid (see Equations (13) and (14), respectively): As shown in Figures 14 and 15, the turbid or clean pixels of the study lakes are distinguished using Diff, GRA, and Tind. The results of plateau lakes are shown in Figure 14, of which only Lake Dianchi is considered to be a turbid lake (Table 1) noting that the pixels in the southern part of Lake Dianchi were masked due to the cyanobacteria bloom. Diff mainly misclassified clean lakes as turbid, such as Lake Yangzonghai, Lake Xingyun, and Lake Chenghai, and locally misclassified the central part of Lake Fuxian and Western Lake Qinghai. Tind performed better than Diff. However, the noise of Tind, which comes from the SL05-1613 nm or SL06-2250 nm in clean water, or the stray light, mainly from the land surface, may affect the judgment in some parts of a lake, such as the western part of Lake Dianchi, Lake Fuxian, and Lake Qinghai. Finally, the GRA proposed in this study can not only effectively classify the whole lake, but can also resist the stray light effect under certain conditions.
In Figure 15, lakes in the middle and lower reaches of the Yangtze River are all sediment-dominated except for Lake Zhuhu ( Table 1). The three indices have less misjudgment in general, but more misjudgment at the local level. Among the classification results of Lake Hongze, Tind is the worst and Diff is second, but the GRA is the best. Although central part of Lake Poyang was distinguished as turbid water by Diff and GRA, pixels of the main channel, in the northern region of Lake Poyang, was mistakenly identified as clean by Tind. For Lake Taihu and Lake Gehu, better discrimination can be found in Diff and GRA, in general, although some misjudgment appeared in the eastern region of Lake Taihu. However, there were more misjudged pixels by Tind in Lake Taihu, especially in the nearshore regions.
In general, it is accurate to discriminate water turbidity using Tind, which was also confirmed in the previous research on its application in coastal waters [8]. However, Tind occasionally misjudged the lake turbidity on the local scale, especially for the nearshore pixels. It can be explained that Tind may not be directly applicable to inland lakes though it performed well in the coastal regions [3,53]. Diff performed well in turbid lakes, but it was vulnerable to misjudge the whole lake in clean water. Moreover, the Diff index has poor resistance to hillside reflections and submerged plants in the inland lakes. However, based on the spectral morphology (Figure 13), the GRA performed better than Diff and Tind because it is supported by the characteristics of the 1020 nm peak of clean water. These findings indicate that the method of turbidity classification based on the spectral morphology is simple, but effective.  Figure 13. Four typical water spectra at 885, 1020, and 1613 nm: a and b denote the clean water in the lake center and nearshore regions, while c and d denote the turbid water in the lake center and the nearshore regions, respectively. Figure 13. Four typical water spectra at 885, 1020, and 1613 nm: a and b denote the clean water in the lake center and nearshore regions, while c and d denote the turbid water in the lake center and the nearshore regions, respectively.

Figure 14.
The results of water turbidity classification of Lake Dianchi, Lake Yangzonghai, Lake Fuxian, Lake Xingyun, Lake Erhai, and Lake Chenghai, which are all plateau lakes in Southwestern China, and Lake Qinghai is located on the Qinghai-Tibet Plateau. Blue and yellow represent clean and turbid waters, respectively.
In Figure 15, lakes in the middle and lower reaches of the Yangtze River are all sedimentdominated except for Lake Zhuhu ( Table 1). The three indices have less misjudgment in general, but 800 900 1000 1100 1200 1300 1400 1500 1600 1700 Wavelength (nm) Figure 14. The results of water turbidity classification of Lake Dianchi, Lake Yangzonghai, Lake Fuxian, Lake Xingyun, Lake Erhai, and Lake Chenghai, which are all plateau lakes in Southwestern China, and Lake Qinghai is located on the Qinghai-Tibet Plateau. Blue and yellow represent clean and turbid waters, respectively.

The Necessity of Using Different Bands for Atmospheric Correction in Turbid and Clean Inland Lakes
To analyze the necessity of band selections for the ACbTC algorithm, six bands were selected, which are the commonly used atmospheric correction bands, from the OLCI and SLSTR: OL16-778.75 nm, OL17-865 nm, OL21-1020 nm, SL05-1613 nm, and SL06-2250 nm [5][6][7]18,53]. Based on the discussion of Sections 5.1 and 5.2, the OL21-1020 nm band was excluded when paired with the NIR band, as the Rayleigh-corrected R rs (1020) value in clean water is sometimes even higher than the R rs in the NIR. Finally, six pairs were generated: OL16 + OL17, OL16 + SL05, OL17 + SL05, OL21 + SL05, OL21 + SL06, and SL05 + SL06, and the atmospheric correction was performed based on these six pairs. Correction results were compared with the in situ measurements and analyzed using negative (the percent of negative values), MAPE, MRPE, and RMSE, as shown in Figure 16. Meanwhile, the full-band average values of these indices are shown in Table 3.
For the OL16 + OL17, there are negative values in both the short and near-infrared bands, especially for Lake Erhai where the match-up points from 400 to 665 nm are all negative, resulting in poor MAPE, MRPE and RMSE values, as shown in Figure 16b-d. These figures indicate that the band combination of NIR bands, 778.75 and 865 nm, will overcorrect even in clean inland water due to the non-negligible water signal. Moreover, the overcorrection caused by OL16 + OL17 is the most serious in this study with the MRPE far less than zero, which results in higher MAPE and RMSE values. However, the correction results are improved when the band of SLSTR is added to the combination, such as OL16 + SL05 and OL17 + SL05. The most considerable improvement of correction is Lake Erhai, but the turbid water, especially Lake Hongze-W and Lake Taihu, are still overcorrected, which means longer bands are required for atmospheric correction in turbid lakes. Although OL21 + SL05 and OL21 + SL06 perform better than OL16 + SL05 and OL17 + SL05, the combinations including OL21-1020 nm result in the critical overcorrection (even the negative results after NIR) in Lake Erhai. The two combinations perform better for turbid water, such as Lake Hongze and Lake Taihu, but are still inferior to the combination of SL05 + SL06, probably due to a more significant amount of stray light from the water at OL21-1020 nm. This result shows again that the impact of the stray light on the accuracy of atmospheric correction for inland lakes should be considered, and weakening its influence is an essential prerequisite for accurate corrections [82].
Remote Sens. 2018, 10, x FOR PEER REVIEW 20 of 29 more misjudgment at the local level. Among the classification results of Lake Hongze, Tind is the worst and Diff is second, but the GRA is the best. Although central part of Lake Poyang was distinguished as turbid water by Diff and GRA, pixels of the main channel, in the northern region of Lake Poyang, was mistakenly identified as clean by Tind. For Lake Taihu and Lake Gehu, better discrimination can be found in Diff and GRA, in general, although some misjudgment appeared in the eastern region of Lake Taihu. However, there were more misjudged pixels by Tind in Lake Taihu, especially in the nearshore regions. Figure 15. The results of water turbidity classification of Lake Hongze in the winter and the spring, as well as Lake Poyang, Lake Gehu and Lake Taihu. Blue and yellow represent clean and turbid waters, respectively.
In general, it is accurate to discriminate water turbidity using Tind, which was also confirmed in the previous research on its application in coastal waters [8]. However, Tind occasionally misjudged the lake turbidity on the local scale, especially for the nearshore pixels. It can be explained that Tind may not be directly applicable to inland lakes though it performed well in the coastal Figure 15. The results of water turbidity classification of Lake Hongze in the winter and the spring, as well as Lake Poyang, Lake Gehu and Lake Taihu. Blue and yellow represent clean and turbid waters, respectively. study, the analysis here is only a simple numerical simulation of these phenomena. More in-depth research requires more synchronously in-situ measurement of atmospheric and water surface, especially in the SWIR bands.
As shown in Table 3, OL17 + SL05 is a reasonable band combination for clean inland lakes such as Lake Erhai, whereas SL05 + SL06 is recommended for the more turbid inland lakes such as Lake Hongze and Lake Taihu.   As shown in Figure 16a,n, Lake Erhai and Lake Taihu present negative results in the NIR while the proportion of Lake Taihu is smaller than that of Lake Erhai, nearing 5%. To analyze this phenomenon, the aerosol reflectance is numerically simulated based on the 6SV model. The atmospheric model was defined by total column water vapor and ozone, the aerosol model was assumed as "Continental", and the aerosol optical depth (AOD) ranged from 0.1 to 1.9 to simulate different aerosol conditions [87]. According to the information provided in Sections 3.1 and 5.1, it is considered that there are three main components in ρ RC , which are the water signal (ρ w ), the aerosol signal (ρ a ), and the stray light signal (ρ S ). The simulated reflectance of water surface was referred to the average of in situ R rs in the nearshore of Lake Taihu. The proportion of the stray light, ρ S /(ρ r + ρ a ), was assumed as a constant 0.1 reference from previous studies [81,88]. When the dark pixel assumption is satisfied, ρ RC = ρ a ; otherwise, ρ RC = ρ w + ρ a + ρ S . Therefore, Equation (6) changes to Equations (15) and (16). Finally, their difference calculated by Equation (17) presents the overcorrection (∆ρ a > 0) or undercorrection (∆ρ a < 0) for the aerosol, respectively: where subscripts T and M denote the theoretical and measured values, respectively [23]. The band combinations including OL16-778.75 nm, OL17-865 nm and OL21-1020 nm (shown in Figure 17a-e) present the overcorrection in short bands which corresponds to Figure 16p. As the reference wavelength increases, the ∆ρ a generally decreases showing that the overcorrection is relieved. However, SL05 + SL06 shows an undercorrection which is relieved with the increasing wavelength ( Figure 17f). The simulation results show once again that OL16 + OL17, OL16 + SL05, OL17 + SL05, OL21 + SL05 and OL21 + SL06 will overestimate the aerosol signal due to the non-negligible water signal and the land stray light, which universally result in the overcorrection. In addition, the points with negative values mainly distributed in the center of Lake Taihu where the water is much cleaner (Figure 6d5) which could explain the overcorrection and the negative corrected values in Figure 16n. The changes of AOD are less sensitive to these band combinations. On the contrary, SL05 + SL06 will cause the undercorrecton which becomes more serious when AOD in 0.5~1.5. However, there are still more interferences in the actual earth observation by satellites. In this study, the analysis here is only a simple numerical simulation of these phenomena. More in-depth research requires more synchronously in-situ measurement of atmospheric and water surface, especially in the SWIR bands.
As shown in Table 3, OL17 + SL05 is a reasonable band combination for clean inland lakes such as Lake Erhai, whereas SL05 + SL06 is recommended for the more turbid inland lakes such as Lake Hongze and Lake Taihu.

Limitations of ACbTC in Applications
Although the ACbTC algorithm has certain advantages compared with other atmospheric correction algorithms, the algorithm itself has certain limitations. The algorithm must be evaluated for applicability of the atmospheric correction.
First, the performance of ACbTC at 400 and 412.5 nm in the dataset of this study was not ideal, with total MAPE values of 43.73% and 41.16%, respectively (Figure 7). Lake Erhai is the primary source of error (shown in Figure 8a) whose MAPE values at 400, 412.5, and 442.5 nm reached 100.21%, 92.41%, and 59.29%, respectively. This means the ACbTC algorithm failed to obtain relatively accurate atmospheric correction results at short wavelengths in clean inland lakes where R rs (778.75) is lower than 0.001 sr −1 .

Limitations of ACbTC in Applications
Although the ACbTC algorithm has certain advantages compared with other atmospheric correction algorithms, the algorithm itself has certain limitations. The algorithm must be evaluated for applicability of the atmospheric correction.
First, the performance of ACbTC at 400 and 412.5 nm in the dataset of this study was not ideal, with total MAPE values of 43.73% and 41.16%, respectively (Figure 7). Lake Erhai is the primary source of error (shown in Figure 8a) whose MAPE values at 400, 412.5, and 442.5 nm reached 100.21%, Figure 17. The differences between measured and theoretical aerosol reflectance, ∆ρ a , for the different band pairs via 6SV simulation. The x-axis represents the wavelength from 400 to 2250 nm. The color bar indicates ∆ρ a in sr −1 . The y-axis represents the AOD, from 0.1 to 1.9, with an interval value of 0.1.
Second, the performance of ACbTC in 753.75-885 nm is not good enough for Lake Taihu, dominated by both suspended sediment and algal particulate matter (Figure 8), and is worse than MUMM (Figure 9g-i), but better in 400-490 nm. This indicates that the atmospheric correction for inland lakes, such as Lake Taihu, is more suitable at short wavelengths using ACbTC, and at NIR wavelengths using MUMM.
Third, for the two seasons of data for Lake Hongze in this study, ACbTC performed best. The MRPE value over the whole band remained at~17.10% and was stable between 510-778.75 nm ( Figure 8). Therefore, the use of ACbTC for atmospheric correction is recommended for inland lakes, such as Lake Hongze (lakes that are dominated by suspended sediments).
Fourth, the upcoming products of Sentinel-3 (such as L1C) are required for synergistic use of OLCI and SLSTR. On the one hand, although the collocation tool of SNAP was used in this study to eliminate the image misregistration of the two sensors, the remaining mismatched pixels (ca. two-pixel distance) will impact the turbidity classification, water pixels extraction, and atmospheric correction [33,89,90]. The future work should discuss the sensitivity of the algorithm to the image misregistration. On the other hand, the first version of the OLCI data was reported to show bias in the blue and NIR so that the SVC is necessary [59,60] for Sentinel-3 products.

Conclusions
An improved atmospheric correction method based on the turbidity classification (ACbTC) algorithm using the synergy of OLCI and SLSTR was presented in this study for inland lakes in China.
(1) The GRA index is proposed to distinguish turbid inland water from clean water, which depends on the feature of the 1020 nm peak of R rs in clean inland lakes and the decreasing slope from 885 to 1613 nm, especially for clean inland lakes that cover relatively small areas. Through the qualitative analysis of 12 inland lakes and the existing observation data, the threshold of the GRA index was determined as −0.07 in this study. The classification results of the GRA index performed better than the Diff and Tind indices for inland lakes by constraining the land adjacency effect at 1020 nm on the similar spectral shapes of turbid lakes.
(2) The synergistic use of OLCI and SLSTR was verified to provide additional atmospheric correction bands, SL05-1613 nm and SL06-2250 nm, for turbid inland lakes to satisfy the dark pixel assumption. The incorrect use of atmospheric correction bands in turbid inland lakes may cause an overcorrection in the short bands due to the invalid assumption and the stray light effect. Thus, OL17 + SL05 and SL05 + SL06 are recommended for the clean and turbid inland waters, respectively. Additionally, OL21-1020 nm is unsuitable for atmospheric correction of inland lakes because it is a vulnerable band easily affected by land stray light. (3) The ACbTC method calculated the aerosol reflectance based on turbidity identification.
Compared to the in situ measurements, the ACbTC algorithm achieved full-band average values of MAPE = 29.55%, MRPE = 13.98%, and RMSE = 0.0039 sr −1 , which were more accurate than the C2RCC, MUMM, FLAASH, POLYMER, and BPAC. The results indicate that both the synergistic use of OLCI and SLSTR and the water turbidity classification are essential for accurate atmospheric correction in inland lakes when using Sentinel-3 data.