An Expanded Three Band Model to Monitor Inland Optically Complex Water Using Geostationary Ocean Color Imager (GOCI)

Due to strict spectral band requirements, the three-band (TB) chlorophyll-a concentration (C chla) estimation algorithm cannot be applied to GOCI image, which has great potential in frequently monitoring inland complex waters. In this study, the TB algorithm was expanded and applied to GOCI data. The GOCI TB algorithm was subsequently calibrated using an in-situ dataset which contains 281 samples collected from 17 inland lakes in China between 2013 and 2020. MERIS TB and GOCI band ratio (BR) models were selected as comparisons to assess the proposed model. The results showed that the proposed GOCI TB model has similar accuracy with MERIS TB model and overperformed GOCI BR model. The root mean square error (RMSE) of the GOCI TB, MERIS TB, and GOCI BR algorithms are 14.212 μg/L, 12.096 μg/L, and 20.504 μg/L, respectively. The mean absolute percentage error (MAPE) (when C chla is larger than 10 μg/L) of the three models were 0.377, 0.250, and 0.453, respectively. Similar conclusion could be drawn from a match-up dataset containing 40 samples. Finally, a simulation experiment was carried out to analyze the robustness of the models under various total suspended matter concentration (C TSM) conditions. Both the in-situ validation and simulation experiment indicated that the GOCI TB factor could effectively eliminate the optical influence of C TSM. Furthermore, the broader spectral range requirement of GOCI TB model made it proper for many other multispectral sensors such as Sentinel two Multispectral Instrument (S2 MSI), Moderate Resolution Imaging Spectroradiometer (MODIS) (onboard the Terra/Aqua satellite), and Visible Infrared Imaging Radiometer Suite (VIIRS) (onboard the National Polar-orbiting Partnership satellite). Compared with the GOCI BR algorithm, the GOCI TB algorithm has stronger stability, better accuracy, and greater potential in practice.

GOCI is the world's first geostationary ocean color multispectral system, with medium spatial resolution (500 m) and very high temporal resolution (1 h). Its high frequency refresh rate provides great potential for inland complex water monitoring (Huang et al., 2015a;Guo et al., 2020). However, in previous studies of remotely estimate of C chla using GOCI images, researchers tend to use empirical models such as BR model, with indistinct mechanisms (Huang et al., 2014;Bao et al., 2015;Huang et al., 2015b). For GOCI data, an accurate and explainable C chla estimation algorithm is urgently needed.
In this study, two questions are addressed. First, could the GOCI BR algorithm correctly describe C chla patterns in optically complex inland waters? Then, is there a clear theory foundation to apply the TB algorithm to GOCI? Focusing on these two questions, the objectives of this paper are 1) extend the TB model to make it proper for GOCI images; 2) calibrate and validate the proposed model in highly turbid inland waters in China using a dataset covering a long time series and a large spatial scale; and 3) assess the robustness and potential of the proposed model.

In-situ Data
Field measurements from 12 cruises, including 281 points, were conducted to calibrate and validate the TB and BR algorithms. The sampling areas covered several turbid productive inland waters in China (Figure 1), including Hongze Lake, Taihu Lake, Dongting Lake, Datong Lake, Changhu Lake, Honghu Lake, Huanggai Lake, Wushan Lake, Poyang Lake, Liangzi Lake, Huangda Lake, and Cihu Lake. The sampling stations are illustrated in Figure 1 and their associate details are listed in Table 1. All measurements and water samples were obtained within 6 h of solar noon (13:00 GMT +8). At each station, the following parameters were measured: remote sensing reflectance (R rs (λ), sr −1 ), C chla , C TSM , and absorption of phytoplankton (a ph (λ)). In 281 samples, 188 of them were randomly selected as calibration dataset and the left 93 samples were used to validate the performance of the C chla estimation models. The in-situ measured R rs and C chla were used for model calibration and validation. In-situ measured C TSM was used to discuss the sensitivity of the models.

Measurement of R rs
R rs (λ) measurements were conducted by means of an ASD FieldSpec spectroradiometer, which has a spectral range of 350-1,050 nm at increment of 1.5 nm. The spectral resolution was interpolated into 1 nm after measurement. According to the above-water measurement method described in the Ocean Optical Protocols (Mueller et al., 2003), the above-water measurement method was used to measure the radiance spectra of the reference panel, water, and sky, respectively. At each site, ten spectra were collected, from which abnormal ones were eliminated and valid ones retained and averaged. Specific observation geometry was applied to effectively avoid the interference of a ship with the water surface and the influence of direct solar radiation during the measurement (Tang et al., 2004). Finally, R rs (λ) was derived via the following equation (Tang et al., 2004;Le et al., 2009) Where L t is the total radiance received from the water surface; L sky is the radiance from sky; L p is the simultaneously observed radiance of the reference gray board. In this process, skylight reflectance at the air-water surface (r) was taken as 2.2% for calm weather, 2.5% for wind speed of up to 5 ms −1 , and 2.6-2.8% for wind speed of about 10 ms −1 (Tang et al., 2004). The reflectance of the gray diffuse board (ρ p ) had been accurately corrected to be 30% in the factory, before we used it to carry out these field experiments.

Measurement of Water Constituents' Concentration
Water samples were collected from the water surface (<20 cm) and kept in a cooler with ice. The fraction of each sample was used to measure concentrations of water constituents. Water samples were filtered through 0.7 μm Whatman GF/F glass fibre filters that had been combusted at 550°C for 4 h. The glass fibre filters were dried at 105°C for 4 h. C TSM was obtained by measuring the difference in weights between combusted and dried glass fibre filters. Water samples for measuring C chla were filtered with GF/C filters (Whatman). Chlorohphyll-a was extracted with ethanol (90%) at 80°C for 6 h in darkness and then analyzed spectrophotometrically at 750 and 665 nm with a correction for phaeopigments using a spectrophotometer (Shimadzu UV-3600) (Chen et al., 2006).

Measurement of a ph (λ)
The water samples were filtered and analyzed by a spectrophotometer (Shimadzu UV-3600) to obtain a ph (λ) using the quantitative filter technique (Cleveland and Weidemann, 1993). First, the water samples were filtered through GF/C (Whatman) glass fibre filters to obtain TSM, and the absorbance of TSM was measured using a spectrophotometer. Pigments were removed from the water samples with NaClO and the water samples were refiltered to obtain tripton. The absorbance of tripton was obtained from these glfibre filters using a spectrophotometer. Data processing to calculate a ph (λ) from the absorbance of tripton and TSM, respectively, was performed as described by Huang et al. (2011).

Satellite Data
Level-1b GOCI images covering Taihu Lake and Hongze Lake were downloaded from the Korea Ocean Satellite Centre (http:// kosc.kordi.re.kr/). For the algorithm validation, the sampling times for the in-situ match-up data were ±0.5 h for the GOCI transit time. There were 40 match-up points (11 from Taihu Lake, Aug. 2013, five from Taihu Lake, Aug. 2019, 10 from Hongze Lake, Apr. 2019, and 14 from Hongze Lake, Nov. 2020) for GOCI images. These images were vector masked to remove land and islands after geometric correction. The GOCI atmospheric correction was carried out using an improved land target-based atmospheric correction method (Guanter et al., 2010;Liu et al., 2015).

Expanding of the TB Algorithm
The TB algorithm that developed by Dall'Olmo and Gitelson (Dall'olmo and Gitelson, 2005;Dall'olmo and Gitelson, 2006) is based on the following relationship between C chla and R rs : where R rs is a function of the inherent absorption (a(λ)) and scattering (b b (λ)) properties of the medium, according to the basic radiative transfer equation (Gordon et al., 1988).
Absorption a(λ) can be separated into absorption related to a ph , a d , a CDOM , and pure water (a w ), while b b (λ) is the measurement of total backscattering. f/Q is depended on Sun zenith angle (Morel and Gentili, 1993), and can be approximated to be 0.0945 (Gordon et al., 1988), and t/n 2 = 0.54 (Austin, 1974;Clark, 1981).
The difference between the reciprocal reflectance R −1 (λ 1 ) and R −1 (λ 2 ) is approximated by based on the following assumptions: (a) b b is spectrally invariant between λ 1 and λ 2 ; (b) a ph (λ 1 ) >> a ph (λ 2 ); (c) a d (λ 1 ) + a CDOM (λ 1 ) ≈ a d (λ 2 ) + a CDOM (λ 2 ). Then, the third band, λ 3 is included to remove the effect of b b . λ 3 is chosen to be in the near infrared (NIR) wavelengths, where reflectance by a ph , a d , and a CDOM is minimal: The a ph could be extracted by combining Eqs 4, 5: Eq. 6 is the basic formula of the traditional TB algorithm, where λ 1 is in the spectral region in which chlorophyll-a shows maximum absorbance (λ 1 = 660-690 nm). λ 2 is in the spectral region in which chlorophyll-a shows minimum absorbance. The absorption of tripton and CDOM at λ 2 are very close to those at λ 1 (λ 2 = 700-730 nm). λ 3 is located in the spectral region such that R rs (λ 3 ) is minimally affected by absorption of water constituents (λ 3 = 740-760 nm). Therefore, in multispectral sensors, only thin-band multispectral sensors like MERIS and S3 OLCI were widely used. According to the previous description, MERIS centre wavelengths for TB algorithm are λ 1 = 681nm, λ 2 = 708nm, and λ 3 = 753 nm. To simplify the discussion, Eq. 6 could be written as: According to some previous studies (Smyth et al., 2006;Chen et al., 2014), the ratio of absorption in different wavelength is fixed: What's more, the GOCI band five is located in 660 nm (Table 1), which is nearby band 6 (680 nm). Then, the following equation can be established: Therefore, the position of λ 2 in the standard TB algorithm could be moved from 710 to 660 nm to fit the GOCI band setting. Eq. 4 could be written as: Then, the expanded TB algorithm could be expressed as the following equation by combing Eq. 5 and Eq. 10 in the assumption that ε aph (680,660) is fixed:

Accuracy Assessment
MAPE and RMSE were used to indicate errors in the estimated values, which could be calculated through the following equations: Where n is the number of samples; y i is the measured value; y' i is the estimated value. In practical, MAPE indexes under different C chla levels fluctuated dramatically. To better reflect the model performance, MAPE low (only the samples that measured C chla less than 10 μg/L were included) and MAPE high (only the samples that measured C chla greater or equal to 10 μg/L were included) were calculated separately.

Simulation Experiment
The main theoretical advantage of TB model is it will eliminate the optical influence of C TSM in NIR spectral range. Therefore, we analyzed this influence based on a bio-optical model calibrated by Huang et al. (2011). In the bio-optical model, R rs was expressed as a function of C TSM , C chla , and a CDOM (λ 0 ) (Appendix A). To reveal the influence of C TSM to C chla estimation models, we simulated R rs spectra between 400 and 800 nm under different C chla and C TSM combinations. C chla and C TSM were ranging from 1 μg/L to 200 μg/L, and 1 mg/L to 200 mg/L, respectively. a CDOM (λ 0 ) was set to 0.5. Then, we defined a variable named model factor difference (Δ factor) to evaluate the stability of model factors to C TSM . The expression of Δ factor is as follows:  (753)), GOCI TB (GTB) (i.e., [R -1 rs(680)-R -1 rs (660)] R rs (745)) factor, and GOCI BR (GBR) factor (i.e., R rs (745)/R rs (680)). All the factors were calculated from R rs at a specific C chla and C TSM level. Δ factor represents the stability of the current factor to C TSM . A Δ factor close to 0 means the current factor is not sensitive to C TSM at the current C chla level.

R rs Spectra Characteristics
The average R rs (Figure 2A) of 15 cruises exhibited wide variability in terms of both magnitude and shape. Spectra of Wushan Lake and Cihu Lake were not plotted because we only collected two samples in these two lakes. This is inadequate to generate the correlation line. For eutrophicationic water, such as Huanggai Lake (Jan. 2018) and Taihu Lake (Aug. 2016), strong phytoplankton absorption formed the R rs peaks around 570, 700, and 815nm, and the trough near 675 nm. For some other cruises such as Changhu Lake (Jan. 2018), its high turbidity and light eutrophication determined the stronger reflectance in red spectral range and weaker trough near 675 nm (Figure 2A). For relatively clean waters that contain less C TSM and C chla , like Liangzi Lake (Jan. 2018), the average R rs curve have small magnitude and weak phytoplankton characters in red and NIR spectral ranges. The correlation coefficients between R rs and C chla at each wavelength ( Figure 2B) also indicated the complexity of the samples. Datasets with significant trough at 680 nm and peak at 700 nm exhibited similar correlation curves. However, for lakes with high turbidity and low eutrophication level, like Changhu Lake, the correlation line looks flat: peak in 700 nm disappeared because high turbidity weakened the R rs sensitivity to C chla .
In conclusion, the complex optical properties of inland water require a C chla estimation model that could eliminate the optical influence of other water color constituents, especially C TSM.

a ph Spectra Characteristics
The average a ph spectrum (Figure 3) contains three peaks. The first one is at 440 nm, which is caused by chlorophyll-a, chlorophyll-b, chlorophyll-c, carotenoid, and other accessory pigments. The second peak at 625 nm is caused by phycocyanin. The last one around 675 nm is mainly caused by the absorption of chlorophyll-a (Haardt and Maske, 1987;Dekker, 1993).
In order to investigate the rationality of the core assumption of the GOCI TB model, i.e., the stability of ε aph , we calculated linear correlation coefficients between a ph at 680 nm and other spectral bands (i.e., line r (680, λ) in Figure 3). Overall, the r (680, λ) line had a similar shape with the average a ph line. This means the main factor that affects inter-band correlation of a ph is its magnitude. In detail, the correlationship became stronger in the spectral range from 400 to 500 nm. Then, a trough appears in 550 nm. In  the wavelength range between 600 and 690 nm, r (680, λ) is entirety greater than 0.9. The strongest correlation appears in 660-690 nm, where the r (680, λ) line is above 0.98. For the bands beyond 690 nm, with the increasing of wavelength, the correlation coefficients rapidly decrease together with the a ph spectrum.
The relationship between a ph (660) and a ph (680) (Figure 4) of all the available in-situ samples denote that when the intercept is set to 0, the linear fitting R 2 is 0.965. This means a ph (660) can explain more than 96% changes of a ph (680) in our dataset. Therefore, the assumption in Eq. 9, i.e., ε aph (680,660) is fixed, is considered to be reasonable in this study (ε aph (680,660) = 1.23).

Calibration and Validation of the Models
Three algorithms were calibrated and validated in this study, they are: MERIS TB algorithm (based on Eq. 7), GOCI TB algorithm (based on Eq. 11), and GOCI BR algorithm R rs (745)/R rs (680) (Huang et al., 2014;Bao et al., 2015) was selected as the model factor. The whole in-situ dataset was randomly separated into 188 calibration samples and 93 validation samples. The parameters of the three models were determined in the calibration dataset by a simple linear fitting method ( Figures 5A,C In the validation dataset, RMSE of the MERIS TB, GOCI TB, and GOCI BR models are 12.943 μg/L, 13.313 μg/L, and 22.613 μg/L, respectively. MAPE low samples are 8. 089,7.725,and 29.898,respectively. MAPE high are 0.433,0.463,and 0.915,respectively (Figures 5B,D,F). The results indicated that GOCI TB algorithm performs better than GOCI BR algorithm and similar with MERIS TB algorithm. Wavelengths in the brackets represent specific band centers of MERIS and GOCI. Table 2 shows the atmospheric correction results of the 40 match-up points. The average MAPE and RMSE over all spectral bands are 0.216 and 4.578 × 10 −3 , respectively. Generally, the atmospheric correction yields satisfactory accuracy, especially for red and NIR bands.

Performance in Match-Up Samples
The C chla estimation results from GOCI TB and BR models represent that, generally, the performance of the two models are similar with those of in-situ dataset ( Figure 5). RMSE of BR and TB models are 12.42 μg/L and 9.90 μg/L, respectively. This means that GOCI TB model improved the accuracy for about 25%. MAPE of BR and TB models are 0.94 and 0.60, respectively. The improvement was about 56%. This suggesting that the GOCI TB model could improve more estimation accuracy in low C chla conditions.

C chla Maps
We applied GOCI TB and BR models to the preprocessed GOCI images and yielded the C chla maps ( Figure 6). Generally, the GOCI TB and GOCI BR C chla maps have similar spatial distribution: high C chla values appeared in Zhushan Bay, Meiliang Bay, north of Gonghu Bay, and south west lake. Low C chla region distributed in the center lake. Notably, an opposite trend appeared between GOCI TB C chla series and GOCI BR C chla series: from 8:30 to 15:30, C chla is getting larger in GOCI TB C chla maps and smaller in the GOCI BR yielded ones.

Error Distribution Analysis
TB model can eliminate the optical influence of total suspended matter in NIR spectral range. To visualize this elimination, MAPE of each sample was calculated for all the three models and plotted in the C chla -C TSM space (Figure 7). The results indicated that, generally, MERIS TB and GOCI TB have similar MAPE distribution. GOCI TB algorithm performs slightly better in low C chla (C chla < 10 μg/L), high C TSM (C TSM > 30 mg/L) region, where GOCI BR result exhibits significant error ( Figure 5E and Figure 7C).
To understand the failure of the GOCI BR algorithm in low C chla , high C TSM condition, we carried out a further analysis. Based on GOCI band settings, the linear correlation coefficients of BR factor with both C chla and C TSM were calculated (Figure 8). What interesting is, the correlation map of C TSM has a similar distribution with that of C chla . Even the overall correlation coefficients of C TSM are lower than C chla , the significant correlation between BR factors and C TSM at 680, 745, and 865 nm indicated the unneglectable influence of suspended matter in BR model. Meanwhile, the correlation coefficients between GOCI TB factors and C TSM was 0.286, which is obviously lower than that of BR factors (higher than 0.5 in Figure 8B). This comparison suggesting that compared with GOCI BR factor, GOCI TB factor could effectively eliminate the impact of TSM and yield a more reliable C chla estimation result.

Theoretical Analysis of the Model Sensitivity to C TSM
The simulation experiment results will be discussed as follows: For each C chla level, we calculated Δ factor under all C TSM  conditions to yield a sensitivity line in Figure 9. For an ideal factor that could completely erase the impact of C TSM , the sensitivity line will be y = 0 (the dashed lines in Figure 9), which means no matter how much TSM is in the water, the spectral factor kept the same. Therefore, the distance of Δ factor line to the dashed line reflects the stability of the current factor.
The results (Figure 9) indicated that the three factors have different stable C chla ranges. Generally, the influence of TSM to C chla spectra factors are increasing together with C TSM . However, For MERIS TB factor, the most stable C chla range is about 70 μg/L ( Figure 9B). High C chla sensitivity and low C chla sensitivity lines equally distributed around the dashed line. GOCI TB factor is insensitive to C TSM when C chla is smaller than 40 μg/L ( Figure 9C). With the increasing of C TSM , offset of the high C chla sensitivity line became larger. Δ GBR has a linear relationship with C TSM . For lower C chla , the factor difference is higher, which indicates that the GOCI BR factor is unstable in low C chla conditions.  In conclusion, the MERIS TB and GOCI TB factors perform similar: they both have stable C chla ranges, which are around 76 μg/L and less than 40 μg/L, respectively. GOCI BR algorithm did not show a stable region: with the increasing of C TSM , the factor became more and more sensitive to TSM. The simulation results can partly explain the model performances in section 4.1.

Comparison of C chla Maps Yielded From GOCI TB and BR Models
Water constituents time series of specific water body has been widely reported (Le et al., 2013;Feng et al., 2014;Palmer et al., 2015). Meanwhile, we proved that different models have different performance under varied conditions (section 3.4). Then, how will this difference affect the spatial-temporal distribution of C chla maps? To answer this question, we calculated the correlation coefficients between GOCI TB and BR yielded C chla maps pixel by pixel ( Figure 10A). A positive correlation coefficient denotes that C chla yielded from TB and BR algorithms have similar temporal changes, and vice versa. The results indicated that in high C chla regions, like Zhushan Bay, Meiliang Bay, Gonghu Bay, and south west lake, two map series showed strong consistency. The correlation coefficients are generally higher than 0.8. However, a large negative correlation region in the center lake means that the two series have totally opposite C chla trends in this region. C TSM sensitivity discussed in section 4.2 can explain this phenomenon: resuspention in the center lake increased C TSM and the instability of GOCI BR model. The histogram of correlation coefficient map ( Figure 11B) demonstrated that the two models yielded opposite trends in a considerable region of Taihu Lake. This will definitely influence the statistical analysis in long-term water C chla monitoring missions.

Influence of Position and Width of Band λ 2 to GOCI TB Model
Our expansion to the TB model is focusing on λ 2 . Following the original ideal of the TB model, we analyzed the effect of band λ 2 to the GOCI TB model by calculating the band-by-band correlation coefficients between C chla and GOCI TB factor.
We fixed the positions of λ 1 and λ 3 and changed band λ 2 from 600 to 740 nm to yield the correlation coefficient line (Figure 11). The widest and highest correlation peak (the first peak in Figure 11) appears at around 710nm, where is close to the theoretical optimal FIGURE 8 | Correlation maps of GOCI BR factors to (A) C chla and (B) C TSM . Marks "*", "**", and "***" represent for p value less than 0.05, 0.01, and 0.001, respectively. band position in previous researches (Dall'olmo and Gitelson, 2005;Le et al., 2011). Therefore, the first peak is the best choice when the sensor could provide proper spectral bands. It is also worth noting that the second peak exists at around 660 nm ( Figure 11). In fact, this peak had been reported by Dall'olmo and Gitelson, 2005 andLe et al. (2009) researches. This suggesting that 660 nm is always a potential selection for TB models.
To further discuss the impact of bandwidth to the GOCI TB model, we resampled R rs (660) to 5, 15, 25, and 35 nm width, respectively. Then, we calculated the GOCI TB factor based on the resampled R rs (660) and their correlation coefficients with C chla (Figure 11). The result denoted that, even the bandwidth was set to 35 nm, the GOCI TB factor still have high correlation coefficient with C chla (larger than 0.87). This means that the GOCI TB factor is not sensitive to the width of band λ 2 .

Potential of the Expanded TB Algorithm
An obstacle of applying traditional TB algorithm is band λ 2 , which located in range of 690-710 nm (Dall'olmo and Gitelson, 2005;Dall'olmo and Gitelson, 2006). By a strict theoretical derivation, the expanded TB model suggests that when the sensor set no bands in 690-710 nm, λ 2 can be moved to 660-690 nm. The expanded TB algorithm has both clear mechanism and satisfactory performance. In fact, apart from GOCI, the proposed model can be also applied to multispectral sensors like S2 MSI, MODIS (onboard the Terra/Aqua satellite), and VIIRS (onboard the National Polar-orbiting Partnership satellite) ( Figure 12).
Therefore, for MODIS (250/500 m band settings), VIIRS, and S2 MSI, their TB and best-fitted BR models were calibrated using the in-situ dataset ( Table 3). In general, all the sensors perform better using TB model. Especially for MODIS and S2 MSI. For    For VIIRS, the band width of λ 2 is so broad ( Figure 12) that it is partly out of the second high correlation peak ( Figure 11). Even though, TB model still improved the performance slightly. Conclusively, the expanded TB algorithm provides a wider choice for accurately monitoring inland complex water.

CONCLUSION
In this research, the TB algorithm is expanded for GOCI image to remotely estimate C chla for optically complex water. The GOCI TB model was calibrated and validated by a comprehensive insitu measured dataset and GOCI images. By comparing with MERIS TB and GOCI BR algorithms, the following conclusions can be drawn: 1) In the expanded TB algorithm, λ 2 is located in the spectral range between 680 and 700 nm or 660-680 nm. The core assumption of the expanded TB algorithm is reasonable. The validation result indicated that GOCI TB algorithm have similar accuracy with MERIS TB algorithm. 2) In highly turbid waters, GOCI BR factor cannot explain the change of C chla . Meanwhile, the GOCI TB algorithm successfully eliminates b b at NIR spectral range and performs well in high C TSM conditions.
3) Compared with the BR algorithm, GOCI TB algorithm can yield more reasonable C chla map for water environment monitoring and analyzing. 4) The expanded TB algorithm is proper for other multispectral sensors like S2 MSI, MODIS, and VIIRS. It provides wider choice for the future inland water remote monitoring missions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
YG performed the data and analysis. YG and CH wrote the paper. CH designed the research. YL, CD, and YL helped to design the experiments. WC contributed to the model development. LS and GJ helped to refine the algorithm.