A Novel Atmospheric Correction for Turbid Water Remote Sensing

Abstract

For the remote sensing of turbid waters, the atmospheric correction (AC) is a key issue. The “black pixel” assumption helps to solve the AC for turbid waters. It has proved to be inaccurate to regard all water pixels in the SWIR (Short Wave Infrared) band as black pixels. It is necessary to perform atmospheric correction in the visible bands after removing the radiation contributions of water in the SWIR band. Here, the modified ACZI (m-ACZI) algorithm was developed. The m-ACZI assumes the spatial homogeneity of aerosol types and employs the BPI (Black Pixel Index) and PIFs (Pseudo-Invariant Features) to identify the “black pixel”. Then, the radiation contributions of waters in the SWIR band are removed to complete the atmospheric correction for turbid waters. The results showed that the m-ACZI had better performance than the SeaDAS (SeaWiFS Data Analysis System) -SWIR and the EXP (exponential extrapolation) algorithm in the visible band (sMAPE < 30.71%, RMSE < 0.0111 sr⁻¹) and is similar to the DSF (Dark Spectrum Fitting) algorithm in floating algae waters. The m-ACZI algorithm is suitable for turbid inland waters.

1. Introduction

Atmospheric correction (AC) is a key step of ocean color data processing. The “black pixel” assumption is commonly applied to AC for waters (including open ocean, coastal, and inland waters) [1,2,3,4,5,6,7]. The radiance of the “black pixel” is assumed to be dominated by atmospheric radiance, and the contribution from water can be neglected; hence, the AC of satellite images can be completed using the “black pixel” assumption [6,8]. This assumption depends on the fact that the water radiation is dominated by the strong absorption of pure water. Early on, waters in near-infrared (NIR) bands were determined to fit this assumption, and the water-leaving reflectance (ρ_w) (or the remote sensing reflectance (R_rs)) in the visible bands was derived by removing the aerosol scattering through extrapolation from the NIR bands [6,8,9].

However, follow-up studies showed that the contribution of turbid water in NIR bands could not be neglected [10,11]. The short-wave infrared (SWIR) bands with a longer wavelength and stronger absorption were considered to fit the “black pixel” assumption for turbid waters [9,12,13,14]. Wang et al. (2019) found that the ρ_w of algae waters was strong in SWIR bands, meaning these bloom waters no longer met the “black pixel” assumption [15]. In Lake Taihu, significant water signals in SWIR bands were detected when algal blooms occurred (Figure 1). Therefore, the contribution of waters in SWIR bands should be removed before the estimation of aerosol scattering in visible bands. Wang et al. (2019) developed the BPI (black pixel index) to obtain the “black pixel” from turbid waters and proposed the ACZI algorithm (atmospheric correction algorithm based on a zero assumption for the short-wave infrared band) for inland turbid waters. The water pixels in the 2201 nm band were considered black pixels by the ACZI algorithm for Landsat 8 images, and aerosol scattering in the visible bands was derived from the assumption of the zero contribution of waters in the 2201 nm band [12]. Figure 1 shows that the Rayleigh-corrected reflectance (ρ_rc) in SWIR bands is very different in water areas without algal blooms; it is incorrect that all water pixels in the SWIR band are treated as black pixels. It is necessary to obtain the “black pixel” with a more accurate identification method.

Here, we modified the identification method of the “black pixel” in the ACZI algorithm [12], and the m-ACZI algorithm was developed. In the m-ACZI algorithm, the water contribution in SWIR bands was removed using the combination of BPI [12] and PIFs (pseudo-invariant features) [16,17], and accurate aerosol scattering in SWIR bands was derived. This scheme was evaluated in turbid inland waters, with its performance also compared with other AC algorithms.

Figure 1. (a) False color image (R: SWIR band, G: NIR band, B: Blue band), (b) Algal bloom map derived by FAI (floating algae index) [18,19], and the Rayleigh-corrected reflectance (ρ_rc) in SWIR bands (c,d) and its histograms from Lake Taihu with algal bloom on 30 January 2021. (e) Average ρ_rc from three boxes (Z1, Z2, and Z3) without algal blooms.

Figure 1. (a) False color image (R: SWIR band, G: NIR band, B: Blue band), (b) Algal bloom map derived by FAI (floating algae index) [18,19], and the Rayleigh-corrected reflectance (ρ_rc) in SWIR bands (c,d) and its histograms from Lake Taihu with algal bloom on 30 January 2021. (e) Average ρ_rc from three boxes (Z1, Z2, and Z3) without algal blooms.

2. Dataset and Methods

2.1. Satellite Data

The cloudless Landsat8-OLI (Operational Land Imager) Level-1 data were from the USGS (United State Geological Survey) and were processed to the R_rs data via the AC algorithms (details in Section 2.3.2). In this work, three OLI images of Lake Taihu (dates: 11 May 2017, 27 May 2017, and 21 December 2017) were matched with the in situ data to evaluate the performance of the AC algorithms. The time window of data matching was set as ±3 h of the Landsat8 overpass. The AC-driven R_rs data were obtained by averaging a 3 × 3 pixel (the coefficient of variation for these valid pixels was < 10%) area surrounding the sample point location [20,21,22,23]. A total of 71 sampling points were matched (11 match-up points from 11 May 2017; 22 match-up points from 27 May 2017; and 38 match-up points from 21 December 2017) (Figure 2).

2.2. Field Spectral Data

The R_rs was measured using the spectrometer (ASD FieldSpec 4) following NASA (National Aeronautics and Space Administration) Ocean optics protocols [24] (Figure 3):

$$R_{rs}=\frac{L_w}{E_d}=\frac{\left(L_t-\sigma {·L}_{sky}\right)·{\rho }_p}{\pi ·L_p}$$

(1)

where L_t is the total water-leaving radiance; ρ_p is the reflectance; L_g is the radiance of the reference panel; L_sky is the sky radiance; L_w is the water-leaving radiance; and E_d is the downwelling plane irradiance. Additionally, the water surface reflectance factor σ was assumed to be 0.028 [25].

2.3. Methods

2.3.1. The m-ACZI Algorithm

The black pixel of the ACZI algorithm was identified using the combination of BPI (Equation (2)) [12] and FAI (the floating algae index) [18,19], and this algorithm assumed that all water pixels in the 2201 nm band accord with the “black pixel” hypothesis. However, the ρ_w in SWIR bands was significant and non-negligible when algal blooms occurred (Figure 1). In this study, the m-ACZI algorithm changed ACZI’s black pixel identification method of the black pixels of ACZI to the BPI-PIF combination, and these black pixels are applied to remove the water contribution in SWIR bands. After this stage, the following process is the same as the ACZI algorithm [12]. The BPI was developed by Wang et al. (2019a) [12] to identify the black pixel from turbid waters (Equation (2)). Schott et al. (1988) [17] considered that urban construction land has pseudo-invariant features (PIFs) and proposed the PIF mask method. The BPI and PIF mask followed the methods described by Wang et al. (2019a) [12] and Schott et al. (1988) [17], respectively.

$$BPI=\frac{\left|{\rho }_{rc}\left(655\right)-{\rho }_{rc}\left(561\right)\right|}{{\rho }_{rc}\left(655\right)-{\rho }_{rc}\left(865\right)}$$

(2)

The PIFs are with features that do not change their reflectivity properties drastically over time in a scene [17]. The surface reflectance of PIFs is similar for the two dates. Hence, the difference between the Rayleigh-corrected reflectance (ρ_rc) of PIFs from two dates can be considered as the difference between two values of aerosol reflectance (ρ_a) of PIF. The black pixels are pixels consisting only of atmospheric radiance (Rayleigh scattering and aerosol scattering). If a pixel is a true “black pixel”, the difference between its ρ_rc for two dates is the same as the difference between its ρ_a for the same two dates (Equations (5)–(7)). Furthermore, theoretically, if the difference in the ρ_rc of a water pixel between two dates in the SWIR band is equal or approximate to the difference in the ρ_rc of PIFs between two dates, this water pixel can be taken as a “black pixel”. Therefore, the characteristics of PIFs can be used to more accurately obtain the BPIs for atmospheric correction. The flowchart of the m-ACZI algorithm is shown in Figure 4.

(1) The ρ_rc and transmittance (t) are obtained from a LUT (lookup table) generated for all OLI bands using 6SV (the vector version of the Second Simulation of the Satellite Signal in the Solar Spectrum) according to sun and sensor geometry [26,27]. The ρ_rc can be derived from the reflectance of the atmosphere (ρ_t). The sea surface reflectance from whitecaps and sunglint were ignored [1,2,28,29].

$${\rho }_{rc}={\rho }_t-{\rho }_r$$

(3)

For land pixels:

$${\rho }_{rc}={\rho }_a+t{\rho }_s$$

(4)

where ρ_s is the surface reflectance.

For a water pixel:

$${\rho }_{rc}={\rho }_a+t{\rho }_w$$

(5)

(2) Two masks are made by the BPI for two date images (T1 and T2) [4,12], and the overlapped pixels in two masks (BPI mask) are retained. The difference in ρ_rc (SWIR) for two date images (D^BPI) is calculated and masked by the BPI mask.

$$\begin{gathered} D^{BPI}={\rho }_{{rc}_{T1}}^{BPI}\left(SWIR\right)-{\rho }_{{rc}_{T2}}^{BPI}\left(SWIR\right) \\ ={(\rho }_{a_{T1}}^{BPI}\left(SWIR\right)+{t_{T1}\rho }_{w_{T1}}^{BPI}\left(SWIR\right))-{(\rho }_{a_{T2}}^{BPI}\left(SWIR\right)+{t_{T2}\rho }_{w_{T2}}^{BPI}\left(SWIR\right)) \end{gathered}$$

(6)

(3) The PIF mask is extracted according to Schott et al. (2016) [16]. The construction land features are considered as PIFs, and this extraction is applied over the construction area within the image (Figure 5). For PIFs, the ρ_s is invariable, and the difference in ρ_rc (SWIR) for two date images (D^PIF) is masked by the PIF mask.

$$\begin{gathered} D^{PIF}={\rho }_{{rc}_{T1}}^{PIF}\left(SWIR\right)-{\rho }_{{rc}_{T2}}^{PIF}\left(SWIR\right) \\ ={(\rho }_{a_{T1}}^{PIF}\left(SWIR\right)+t_{T1}{\rho }_{s_{T1}}^{PIF}\left(SWIR\right))-{(\rho }_{a_{T2}}^{PIF}(SWIR)+{t_{T2}\rho }_{s_{T2}}^{PIF}(SWIR\left)\right) \\ ={\rho }_{a_{T1}}^{PIF}\left(SWIR\right)-{\rho }_{a_{T2}}^{PIF}\left(SWIR\right)+\Delta t{\rho }_s^{PIF}\left(SWIR\right) \end{gathered}$$

(7)

Δt will make ρ_S^PIF residual in D^PIF and affect the matching of D^PIF and D^BPI; for this, we analyze and discuss the impact of Δt on D^PIF (1609) (details in Section 4.4). The results showed that when Δt was less than 0.0107, the D^PIF change caused by Δt was about 2.70%, which is within the allowable error range.

(4) The ρ_a^BPI and ρ_a^PIF are similar in one OLI image based on the assumption of the spatial uniformity distribution of aerosol types. The D^BPI matches the mode number of D^PIF (with ±10% variation) to obtain the final BPI-PIF mask. The black pixels were identified using the BPI-PIF mask (Example: Figure 6). The aerosol scattering from black pixels is used to remove the water contribution for nonblack pixels.

$$D^{BPI}\cong D^{PIF}$$

(8)

(5) The atmospheric correction for all water pixels is completed using the same steps as the ACZI algorithm [12], and R_rs (or ρ_w) data are obtained.

Described below is the process of developing the m-ACZI algorithm.

(1)

We extracted the PIF mask according to Schott et al., 1988 [17] and Concha and Schott 2016 [16]. For Lake Taihu, the 2201 nm band with a threshold of 0.05 was used to obtain the water mask, and the threshold of the NIR/RED band ratio was set to 1.3 for obtaining the vegetation mask. The previous two masks were combined using logic. AND. gate, creating a “PIF mask” that rejects the vegetation and water and only accepts the construction land features.

(2)

In this step, we calculated the ρ_rc and t using the 6SV-LUT from ACOLITE (atmospheric correction for OLI lite).

(3)

Based on step 2, we adopted the BPI [12] to obtain the BPI mask, and the threshold of the BPI was set to 0~0.1.

(4)

Following steps 1, 2, and 3, concerning t values, two images with similar t values were selected to obtain their PIF mask and BPI mask, respectively. Then, these masks were applied to the ρ_rc (SWIR) images, respectively, to obtain the PIF and BPI images. We subtracted two of each image type to obtain the D^PIF and D^BPI images.

(5)

We counted the average value of D^PIF and used the average D^PIF with 10% variation as the judgment value for D^BPI, i.e., D^BPI pixels greater than the average D^PIF with 10% variation were discarded.

(6)

Following the ACZI process [12], we counted the median value from the final retained D^BPI pixels, and this median value was taken as ρ_a on two SWIR bands. The aerosol scattering ratio (ε), consisting of these two ρ_a(SWIR) values, was applied to all water pixels, and aerosol scattering in the visible and near-infrared bands was obtained using an exponential extrapolation method, thus completing the final atmospheric correction process to obtain R_rs.

2.3.2. Other Atmospheric Correction Algorithms

The SeaDAS-SWIR Algorithm

For turbid waters, the “black pixel” assumption in NIR bands is invalid because of the contribution of phytoplankton, detritus, and suspended sediment to water backscatter [2,10]. Research from Hale et al. (1973) indicated that the water in SWIR bands has stronger absorption than that in NIR bands [30], and the turbid water in the former meets the requirement of the “black pixel” assumption [13]. Therefore, Wang and Shi (2005) [10] developed an AC algorithm for turbid water based on extrapolation from SWIR bands and integrated this AC algorithm into SeaDAS (SeaWiFS Data Analysis System, referred to as SeaDAS-SWIR in this paper; more details can be found in Wang and Shi (2005) [10]).

The ACOLITE (atmospheric correction for OLI lite; this study worked with version 20211112) tool was developed for coastal and inland waters by the Remote Sensing and Ecosystem Modelling (REMSEM) team. It performs this process using dark spectrum fitting (DSF) and exponential extrapolation (EXP) [14,26,31,32,33].

The EXP Algorithm

The idea of the EXP algorithm is consistent with the SeaDAS-SWIR algorithm, both of which are based on exponential extrapolations of the ratio of multiple aerosol scattering in SWIR bands, but the EXP algorithm does not deal with aerosol LUTs (lookup tables) [14]. The EXP algorithm only uses an exponential extrapolation method to estimate aerosol scattering in visible bands and NIR bands [28].

The DSF Algorithm

The DSF algorithm is presented for meter-scale-resolution optical satellite imagery [26,31]. Similar to other AC algorithms, this algorithm was also developed based on the “black pixel” concept [26]. However, the difference is that it defines black pixels as pixels with zero surface reflection in one scene image, including land pixels and water pixels. These black pixels are employed to estimate the atmospheric radiance for other pixels. It can be seen that the DSF depends on two assumptions: (1) the spatial distribution of the atmosphere is uniform in an image; (2) at least one pixel whose surface reflection is about zero exists in one scene image [26,31].

2.3.3. Algorithm Accuracy Analysis

The root-mean-square error (RMSE), symmetric mean absolute percentage error (sMAPE), and Bias were used to determine the differences between the in situ data and the AC-driven data:

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left(R^m-R^i\right)}^2}{n}}$$

(9)

$$sMAPE=\frac{{\sum }_{i=1}^n\frac{\left|R^m-R^i\right|}{\left|R^m+R^i\right|}}{n}\times 2\times 100\%$$

(10)

$$Bias=\frac{{\sum }_{i=1}^n\left[\frac{\left(R^m-R^i\right)}{R^i}\right]}{n}\times 100\%$$

(11)

where Rⁱ is the in situ R_rs data, R^m is the OLI-R_rs data via atmospheric correction, and n is the number of match-up pairs.

3. Results

3.1. Assessment of the m-ACZI Algorithm

As Figure 7 and Figure 8 show, the m-ACZI raised the R_rs estimation compared to the ACZI algorithm, so the histograms of R_rs generally moved and concentrated to high value areas (Figure 8). Notably, the m-ACZI algorithm did not increase the R_rs values for all water pixels, but mainly for algal bloom pixels (green pixels in the RGB true-color images) (Figure 2 and Figure 7). It can also be observed that the m-ACZI algorithm reduced the R_rs estimation of the non-algal bloom pixels in the 865 nm band (Figure 7), which were consistent with the characteristics of the in situ R_rs (Figure 3).

Most of the scatterplots of the ACZI were below the 1:1 line, and part of the scatterplots derived in the 865 nm band were above the 1:1 line (Figure 9a). The scatterplots of the m-ACZI were distributed around the 1:1 line, with several significant off-group points in the 865 nm band (Figure 9b). In general, compared with the scatterplots of the ACZI, the scatterplots of the m-ACZI were more concentrated near the 1:1 line. The ACZI-driven R_rs was underestimated in all visible bands. The R_rs derived by the m-ACZI was also underestimated in the 561 nm and 655 nm bands, but overestimated in two BLUE bands (Figure 9 and

Table 1. Statistics of the accuracy measures for the ACZI and the m-ACZI algorithms using in situ R_rs measurements. The BLUE number represents the minimal statistical value.

Method		443	483	561	655	865
m-ACZI	RMSE(sr−1)	0.0063	0.0062	0.0111	0.0082	0.0228
	sMAPE(%)	30.71	23.01	21.41	21.27	48.31
	Bias(%)	52.07	7.59	−16.48	−11.65	88.83
ACZI	RMSE(sr−1)	0.0070	0.0076	0.0146	0.0102	0.0203
	sMAPE(%)	35.88	34.07	29.16	29.10	57.48
	Bias(%)	−8.91	−12.32	−22.70	−16.18	94.16

). The minimum sMAPE was obtained by the m-ACZI in the 655 nm band, followed by the m-ACZI in the 561 nm band. Similarly, the minimum RMSE was provided by the m-ACZI in the 483 nm band, followed by the m-ACZI in the 443 nm band. However, the minimum bias was obtained from the data derived by the ACZI in the 443 nm band (Table 1). Overall, from RMSE and sMAPE, the m-ACZI was superior to the ACZI in all visible bands, especially in the 561 nm and 655 nm bands; the performance of the m-ACZI algorithm was greatly improved. However, these 2 algorithms still failed in the NIR band (RMSE > 0.0203 sr⁻¹; sMAPE > 48.31%; bias > 88.83%). From the scatterplots of the 865 nm band in Figure 9, these 2 AC algorithms have overestimations. The assessment results indicate that the m-ACZI algorithm is generally better than the ACZI algorithm.

3.2. Comparison with other AC Algorithms

Only 54 match-up sample points were used to evaluate the SeaDAS-SWIR algorithm because some water pixels were masked by this algorithm (the “proc_ocean” option of SeaDAS was set to “2-force all pixels to be processed as ocean”). In Figure 10, the scatterplots of the m-ACZI algorithm were more concentrated, while those of other algorithms were discrete, especially in the 443 nm and 483 nm bands. The scatterplot distributions of the m-ACZI and the DSF were very similar in the 561 nm and 655 nm bands and were closer to the 1:1 line than others (Figure 10). The statistical results of the accuracy also supported these findings (

Table 2. Statistics of the accuracy measures for the m-ACZI, SeaDAS, DSF, and EXP algorithms using in situ R_rs measurements. The BLUE number represents the minimal statistical value.

Method		443	483	561	655	865
m-ACZI	RMSE(sr−1)	0.0063	0.0062	0.0111	0.0082	0.0228
	sMAPE(%)	30.71	23.01	21.41	21.27	48.31
	Bias(%)	52.07	7.59	−16.48	−11.65	88.83
SeaDAS-SWIR	RMSE(sr−1)	0.0113	0.0100	0.0108	0.0104	0.0134
(n = 54)	sMAPE(%)	45.11	32.70	20.65	27.33	76.35
	Bias(%)	33.32	13.36	−1.38	4.37	126.68
EXP	RMSE(sr−1)	0.0060	0.0070	0.0150	0.0105	0.0217
	sMAPE(%)	34.12	27.76	32.79	28.29	55.56
	Bias(%)	2.50	−10.31	−27.00	−19.72	90.38
DSF	RMSE(sr−1)	0.0058	0.0061	0.0098	0.0080	0.0244
	sMAPE(%)	32.27	24.34	17.68	21.80	63.00
	Bias(%)	34.22	4.94	−4.46	0.01	113.55

). The assessment results from the m-ACZI and the DSF were very similar: the RMSE difference between them was less than 0.0008 in visible bands, and both were more accurate than the SeaDAS and the EXP algorithms. It should be noted that all the AC algorithms failed in the NIR band (RMSE > 0.0134 sr⁻¹; sMAPE > 49.31%; Bias > 88.83%), and the m-ACZI algorithm showed the best comprehensive performance.

3.3. Assessment of AC Algorithms for Different Water Types

According to the water classification method of Sun et al. (2012), waters were classified into two types: floating bloom (FB) and un-floating bloom (uFB) [15,34]. Figure 11 shows the scatterplots of R_rs between in situ data and the AC algorithm in different water types, with detailed statistics presented in

Table 3. Statistics of the accuracy measures for atmospheric correction algorithms for un-floating blooms using in situ R_rs measurements. The BLUE number represents the minimal statistical value.

Method		443	483	561	655	865
m-ACZI	RMSE(sr−1)	0.0058	0.0058	0.0087	0.0072	0.0127
(n = 53)	sMAPE(%)	25.63	20.49	18.97	17.86	48.62
	Bias(%)	21.32	3.93	−14.76	−10.45	91.30
ACZI	RMSE(sr−1)	0.0071	0.0073	0.0109	0.0085	0.0112
(n = 53)	sMAPE(%)	31.77	25.94	23.13	22.59	59.54
	Bias(%)	−6.01	−6.81	−18.20	−10.54	102.18
SeaDAS-SWIR	RMSE(sr−1)	0.0113	0.0101	0.0107	0.0103	0.0123
(n = 48)	sMAPE(%)	44.17	31.74	20.02	26.57	74.39
	Bias(%)	30.91	11.20	−2.26	2.61	119.49
EXP	RMSE(sr−1)	0.0058	0.0068	0.0126	0.0093	0.0123
(n = 53)	sMAPE(%)	26.63	24.18	30.28	24.10	55.98
	Bias(%)	−3.30	−10.77	−25.36	−17.42	93.74
DSF	RMSE(sr−1)	0.0055	0.0058	0.0075	0.0068	0.0144
(n = 53)	sMAPE(%)	26.92	21.81	15.16	19.31	65.51
	Bias(%)	8.72	3.16	−1.85	1.95	117.79

and

Table 4. Statistics of the accuracy measures for atmospheric correction algorithms for floating bloom using in situ R_rs measurements. The BLUE number represents the minimal statistical value.

Method		443	483	561	655	865
m-ACZI	RMSE(sr−1)	0.0078	0.0072	0.0162	0.0107	0.0396
(n = 18)	sMAPE(%)	45.68	30.42	28.61	31.29	47.40
	Bias(%)	142.60	18.39	−21.53	−15.18	81.53
ACZI	RMSE(sr−1)	0.0067	0.0086	0.0222	0.0140	0.0355
(n = 18)	sMAPE(%)	51.47	58.00	46.92	48.25	51.41
	Bias(%)	−19.90	−28.52	−35.93	−32.81	70.54
SeaDAS-SWIR	RMSE(sr−1)	0.0111	0.0097	0.0115	0.0110	0.0204
(n = 6)	sMAPE(%)	52.71	40.43	25.69	33.33	92.05
	Bias(%)	52.55	30.56	5.67	18.46	184.26
EXP	RMSE(sr−1)	0.0067	0.0075	0.0205	0.0134	0.0376
(n = 18)	sMAPE(%)	56.19	38.30	40.20	40.63	54.32
	Bias(%)	19.58	−8.93	−31.81	−26.47	80.47
DSF	RMSE(sr−1)	0.0067	0.0070	0.0147	0.0107	0.0416
(n = 18)	sMAPE(%)	48.00	31.78	25.10	29.13	55.60
	Bias(%)	109.31	10.17	−12.16	−5.71	101.07

. All AC algorithms had more discrete scatterplots in FB than in uFB, especially in the two blue bands, where the bias of the algorithm increased significantly (Table 3 and Table 4). Table 3 shows that the assessment results of the m-ACZI and the DSF algorithms were superior to other algorithms in the uFB. In particular, the m-ACZI and the DSF algorithms showed satisfactory performance in the 561 nm and 655 nm bands (RMSE < 0.0087 sr⁻¹; sMAPE < 19.31%), which was sufficient to support the quantitative remote sensing research for turbid waters [5]. Compared with those in the uFB, the accuracy results of the AC algorithms in FB were reduced (Table 3 and Table 4). The DSF algorithm had the best performance in the FB waters, followed by the m-ACZI algorithm. The SeaDAS-SWIR algorithm masked most of the FB pixels, with only 6 in situ points (18 in situ points in total) participating in the evaluation. The performance of these AC algorithms was not satisfactory, especially in the two blue bands, as was that of the m-ACZI algorithm. These AC algorithms (except the ACZI algorithm) overestimated the water-leaving reflection in the 443 nm and 483 nm bands (Table 4). A possible reason was that the m-ACZI algorithm for the FB underestimated the aerosol scattering in the SWIR band, which was further exacerbated by exponential extrapolation. This finally resulted in a serious underestimation of the aerosol scattering and an overestimation of the water-leaving reflection in the 443 nm and 483 nm bands.

4. Discussion

The m-ACZI algorithm depends on two prerequisite assumptions: (1) the spatial distribution of aerosol types is uniform; (2) the black pixels from two date images overlap. These are the determining conditions for the success of the m-ACZI algorithm.

4.1. The Assumption of the Spatially Uniformity Distribution of Aerosol Types

Solving the aerosol scattering of nonblack pixels relies on the assumption of the spatially uniform distribution of aerosol types for the m-ACZI algorithm. This assumption has been recognized by a great amount of research [1,4,14,26,28]. For example, the DSF and EXP adopted this assumption for the atmospheric correction process. Taking the Lake Taihu image without algal blooms (obtained on 2 February 2016) as an example (remove cloud pixels) (Figure 12), the spatial distributions of ρ_rc (1609) and ρ_rc (2201) were uniform, and the standard deviations of ρ_rc (1609) and ρ_rc (2201) were 0.0026 (average: 0.0218) and 0.0016 (average: 0.0137), respectively. The standard deviation of ρ_rc (1609)/ρ_rc (2201) was 0.0546 (average: 1.5944) (Figure 12). This suggests that Lake Taihu is consistent with the assumption of the spatially homogeneous distribution of aerosol types.

4.2. The Black Pixels from Two Date Images Overlap

The m-ACZI algorithm searches and obtains the “black pixel” based on the combination of BPI and PIF for two date images, which requires the BPI and PIF pixels of the two periods to overlap. The PIF pixels are from urban construction land, which is invariable in a year. Therefore, the PIF pixels of the two periods should overlap. The occurrence probability of BPI pixels is calculated in Lake Taihu based on cloudless images from 2013 to 2018 (a total of 36 images) (Figure 13), showing more than 30% in a large area of central Lake Taihu. This means that there is a very high probability of obtaining overlapping BPI pixels from three images, which is sufficient to support the m-ACZI algorithm. Landsat 8-OLI has accumulated a large number of images to complete the m-ACZI algorithm.

4.3. The m-ACZI Algorithm Depends on Pure Pixels

The m-ACZI algorithm needs PIF pixels to identify the “black pixel” because they do not drastically change their reflectivity properties over time. The reflectivity properties of PIFs are derived from construction land, and mixing from other variable land types (such as water, forest land, and agricultural land) weakens the invariant reflectivity properties of PIFs. The theoretical basis of the m-ACZI algorithm requires minimizing the uncertainty caused by a mixture of pixels, and the high spatial resolution can reduce the uncertainty of a mixture of pixels [35]. From this point of view, the m-ACZI algorithm depends on pure pixels.

4.4. The Impact of Δt on D^PIF

To analyze the influence of the Δt on D^PIF, the t was calculated using 6SV-LUT (from ACOLITE software (version 20211112)) based on cloudless images from 2013 to 2018 (a total of 36 images). In the 1609 nm band, the variation range of t was 0.9981 to 0.9543, the average t was 0.9632, the standard deviation (SD) of t was 0.0111 (Figure 14), the range of Δt was 0.04375 to 0.00002, the average Δt was 0.00971, and the SD of Δt was 0.01255. In the 2201 nm band, the variation range of t was 0.9663 to 0.8982, the average t was 0.9149, the SD of t was 0.0165 (Figure 14), the range of Δt was 0.06810 to 0.00003, the average Δt was 0.01656, and the SD Δt was 0.01699. The variation in t and Δt in the 1609 nm band was less than in 2201 nm; this band was selected as the judgment band for D^PIF ≈ D^BPI. Suzhou City (the red box in Figure 15a) near Lake Taihu was taken as an example to illustrate the influence of the Δt on D^PIF (1609). Furthermore, to ensure the independence of verification, ρ_s was driven by the ACOLITE. To reduce the error caused by the ACOLITE-AC algorithm, the average ρ_s^PIF was considered as the real ρ_s^PIF (Figure 15c). This ρ_s^PIF deviated from the true value, but it can still be considered as a reference for analyzing the influence of Δt on D^PIF. Since there were no images that completely corresponded to the average Δt value, we selected the corresponding two images of the approximate average Δt for analysis. The maximum Δt (Δt = 0.04375) in the 1609 nm band was from 2 images from 12/29/2014 and 05/14/2018, and the Δt = 0.01065 (approximate average Δt) in the 1609 nm band was from 2 images from 12/29/2014 and 06/15/2018. As shown in Figure 15, when Δt = 0.04375, the average Δtρ_s^PIF was 0.0058, the median Δtρ_s^PIF was 0.0055, the average Δtρ_s^PIF/D^PIF was 14.85%, and the median Δtρ_s^PIF/D^PIF was 10.84%. When Δt = 0.01065, the average Δtρ_s^PIF was 0.0014, the median Δtρ_s^PIF was 0.0013, the average Δtρ_s^PIF/D^PIF was 4.62%, and the median Δtρ_s^PIF/D^PIF was 2.68%. It can be seen that Δt can introduce an error into the algorithm, but when Δt is less than 0.01065 (this was a common occurrence), the error it introduced was small and within the allowable error range. We propose selecting two images with similar t values for the m-ACZI algorithm, and this t value can be pre-estimated by the LUT.

5. Conclusions

The “black pixel” assumption plays a key role in providing atmospheric (especially aerosol) scattering information to solve the atmospheric correction process for water. We modified the identification method for the “black pixel” in the ACZI algorithm as a combination of BPI and PIF and provided the m-ACZI algorithm. The m-ACZI algorithm assumes the spatial homogeneity of aerosol types and employs BPI and PIF to identify the “black pixel”, which together make this algorithm suitable for inland water (such as lakes, rivers, and reservoirs). Compared with the ACZI algorithm, m-ACZI had better performance in the visible bands, especially in the 561 nm and 655 nm ones. The evaluation results showed that the performance of the m-ACZI algorithm is superior to the SeaDAS-SWIR, ACZI, and EXP algorithms in waters with algal blooms. This algorithm is a useful attempt to improve the accuracy of AC for turbid water and provides a new scheme for obtaining black pixels from turbid inland waters.