Comparison of Aerosol Reflectance Correction Schemes Using Two Near-Infrared Wavelengths for Ocean Color Data Processing

This paper reanalyzes the aerosol reflectance correction schemes employed by major ocean color missions. The utilization of two near-infrared (NIR) bands to estimate aerosol reflectance in visible wavelengths has been widely adopted, for example by SeaWiFS/MODIS/VIIRS (GW1994), OCTS/GLI/SGLI (F1998), MERIS/OLCI (AM1999), and GOCI/GOCI-II (A2016). The F1998, AM1999, and A2016 schemes were developed based on GW1994; however, they are implemented differently in terms of aerosol model selection and weighting factor computation. The F1998 scheme determines the contribution of the most appropriate aerosol models in the aerosol optical thickness domain, whereas the GW1994 scheme focuses on single-scattering reflectance. The AM1999 and A2016 schemes both directly resolve the multiple scattering domain contribution. However, A2016 also considers the spectrally dependent weighting factor, whereas AM1999 calculates the spectrally invariant weighting factor. Additionally, ocean color measurements on a geostationary platform, such as GOCI, require more accurate aerosol correction schemes because the measurements are made over a large range of solar zenith angles which causes diurnal instabilities in the atmospheric correction. Herein, the four correction schemes were tested with simulated top-of-atmosphere radiances generated by radiative transfer simulations for three aerosol models. For comparison, look-up tables and test data were generated using the same radiative transfer simulation code. All schemes showed acceptable accuracy, with less than 10% median error in water reflectance retrieval at 443 nm. Notably, the accuracy of the A2016 scheme was similar among different aerosol models, whereas the other schemes tended to provide better accuracy with coarse aerosol models than the fine aerosol models.


Introduction
In the last few decades, ocean observations using visible (VIS) to near-infrared (NIR) wavelength satellite imagery have successfully retrieved oceanic environmental information over a large spatial and temporal range [1][2][3][4][5][6][7][8].Atmospheric correction is a primary process in satellite ocean color missions, in which surface water reflectance, ρ wn , is extracted from the top of atmosphere (TOA).The quantity of ρ wn is generally less than 10% of the TOA reflectance (ρ TOA ); thus, atmospheric reflectance and transmittance should be computed as accurately as possible.Ignoring the surface-reflectance from sun-glint and whitecaps, ρ TOA at each wavelength (λ) can be described as follows [9,10]: where ρ r is the Rayleigh multiple-scattering (molecular multiple-scattering) reflectance in the absence of aerosols, and ρ a + ρ ra is the aerosol multiple-scattering reflectance in the presence of air molecules (denoted ρ am ).Terms t s d and t v d are the diffuse transmittances of the atmosphere from the sun to the sea surface and from the sea surface to the sensor, respectively.
Traditional atmospheric correction processes first compute ρ r (λ) within 1% error using a radiative transfer simulation [11][12][13][14][15].Following subtraction of the Rayleigh reflectance contribution, ρ am (NIR) can be estimated by assuming that the NIR water-leaving reflectance is negligible, i.e., ρ wn (NIR) = 0, due to the relatively strong light absorption of the waterbody or so-called "black pixel" assumption.ρ am (VIS) is then extrapolated from satellite-observed ρ am (NIR) by inversion of the radiative transfer simulation.
The inverse processes for aerosol correction are generally separated into two steps: aerosol model selection, and subsequent spectral extrapolation of ρ am (VIS) from ρ am (NIR).In the model selection step, two similar aerosol models, and the contributions of their weighting factors, are extrapolated from ρ am (NIR).ρ am (VIS) is then computed from ρ am (NIR) using the aerosol information estimated in the extrapolation step.
Schemes GW1994, F1998, and AM1999, but not A2016, were compared and evaluated for a coarse maritime aerosol model with a relative humidity (RH) of 80%, and presented in a report by the International Ocean Colour Coordinating Group (IOCCG) [21].However, the IOCCG comparison focused more on differences in atmospheric correction systems than on differences of the aerosol correction schemes.Since each scheme used different aerosol models and different radiative transfer code and different Mie scattering code, it is difficult to determine which factors account for the different performances between the schemes.
In this paper, we first revisit and compare the above-mentioned schemes, including the A2016 scheme, and then present our results, in which the performance of each algorithm was examined using a simulated TOA dataset generated by radiative transfer simulations [22][23][24].Our analysis was performed using aerosol models with a size distribution including fine to coarse aerosols, whereas the IOCCG study only looked at coarse aerosols in their models.For clarification of the comparison, we used the same code and model for the radiative transfer, Mie scattering calculations, and aerosol models for both the implementation and the analysis result.

The Two-NIR Aerosol Correction Methods Used by the Various Schemes
This section describes various two-NIR aerosol correction schemes.The correction schemes and look-up tables were developed simultaneously to allow comparison (as summarized in Table 1).Fundamentally, F1998, AM1999, and A2016 were developed based on GW1994; the principle differences between the listed schemes lie in the details of the aerosol model and weighting factor determination.Specifically, the GW1994 scheme uses the spectral ratio of the aerosol single-scattering reflectance of two NIR bands, and the F1998 scheme uses the spectral ratio of the aerosol optical thickness (τ a ) of two NIR bands}.Both the AM1999 and A2016 schemes apply the spectral ratio of aerosol multiple-scattering reflectance of two NIR bands; however, A2016 is slightly modified to consider the spectrally varying weighting factor [19,20].In the AM1999 scheme, the path reflectance increment ratio term (ρ am + ρ r )/ρ r by aerosol concentration is utilized, whereas other schemes use ρ am .

GW1994 Scheme
The GW1994 scheme determines aerosol model contributions in the single-scattering domain.Therefore, GW1994 first converts ρ am into single-scattering aerosol reflectance (ρ as ) for two NIR bands (hereafter, the shorter wavelength NIR band is denoted by NIR S , and the longer wavelength NIR band is denoted by NIR L ) for all i-th candidate aerosol models (M i ), as follows: where F m→s is the empirical function for converting ρ as to ρ am for specific λ, M i , and Θ.The Θ is the scanning geometry variable that is a combination of solar zenith angle, satellite zenith angle, and relative azimuth angle values [9,10,21,25].
ρ as values are then used to calculate the single-scattering reflectance ratio (ε as ), used for selecting the appropriate aerosol models from among the candidate models (M i ): Each M i has a theoretical eigenvalue of ε as (denoted as ε M as ) derived from the analytical single-scattering reflectance model, determined by the single-scattering albedo of M i and the phase function.The GW1994 scheme selects two similar aerosol models, M L and M H , by comparing the representative ε as (denoted as ε r as ) to ε M as , as follows: The representative ε as (NIR S , NIR L ) can be approximated by averaging over the ε M as of all candidate models, as the conversion function F m→s is less sensitive to the aerosol model for NIR wavelengths: Then, the weighting factors w M L and w M H for the two selected aerosols (M L and M H ) can be estimated by the fraction of the linear distance between ε r as and ε M as : The ρ as (VIS) for the two selected models can be extrapolated spectrally from ρ as (NIR) based on ε M as (NIR, VIS).Having derived the values w M L and w M H for M L and M H , respectively, the aerosol multiple-scattering reflectance (ρ am ) in the VIS bands is given by: where F s→m is the inverse function of F m→s [9,10,21,25].

F1998 Scheme
The F1998 scheme determines the contributions of the aerosol models to the τ a domain [16], and has the advantage that the model-wise single-scattering reflectance ratios are often unevenly distributed over the single-scattering space [26].For model determination, F1998 uses the inter-band τ a ratio of two NIR bands, as follows: where ε τa is the ratio of τ a of two wavelengths.Therefore, the F1998 scheme first converts ρ am into τ a for all candidate models, M i , using an empirical conversion function (F m→τ ) between ρ am and τ a , in contrast to Equation ( 2), as follows: The F m→τ function can be expressed in several ways [18,19,27,28], and F1998 uses a third-order polynomial function to describe the empirical relationship [16].
Theoretically, each M i has a scan geometry-independent eigenvalue of ε M τa that can be expressed by the ratio of the aerosol extinction coefficient K ext (M i , λ), as follows: In a similar way to Equation ( 4), the contributions of M L and M H can be determined by comparing the representative ε τa (denoted as ε r τa ) to ε M τa , as follows: To enhance the accuracy in determining the two aerosol models (M H and M L ) and their corresponding weighting factors (w M L and w M H ), the F1998 scheme uses a weighted-average value of ε τa for all M i : In the last step, ρ am for all VIS bands can be computed similarly to Equation ( 8), as follows: where F τ →m is the inverse function of F m→τ [16].

AM1999 Scheme
In the AM1999 scheme, the atmospheric path reflectance ratio to Rayleigh reflectance (ρ χ ) is used for the aerosol estimation, whereas other general aerosol correction schemes compute ρ r and ρ am separately, as follows: According to Antoine and Morel [18,28], ρ χ (λ 1 ) for M i can be modeled from ρ χ (λ 2 ) in a similar way to F1998, using a ρ χ versus τ a relationship in which Equation ( 10) and ( 18) takes the following form: where F χ→τ is the conversion function from ρ χ to τ a , expressed as a quadratic equation, and F τ →χ is the inverse function of F χ→τ [18,19,28].The value ρ M χ is the modeled ρ χ for candidate model M i .To select the two-closest aerosol models and determine the corresponding weighting factor, ρ χ (NIR S ) for all candidate aerosol models ρ M χ (M i , NIR S ) are first computed using a quadratic expression describing the relationship of ρ χ with τ a and K ext (M i , λ) [28].Then, M L and M H are selected by comparing ρ χ (NIR S ) observed by satellite to the ρ M χ (M i , NIR S ) of all candidate models, as follows: The mixing ratio w M H is then derived directly in the multiple-scattering domain: Using the two models (M L and M H ) and their derived weighting factors (w M L and w M H , respectively), ρ χ for all VIS bands can be computed as follows: where ρ M χ (M L , λ) and ρ M χ (M H , λ) are spectrally extrapolated from ρ χ (NIR L ) [29].

A2016 Scheme
To select the appropriate aerosol models and perform spectral extrapolation of their reflectance from the NIR to VIS bands, the A2016 method directly estimates the reflectance fraction of the two-closest models in the multiple-scattering domain, without considering the single-scattering domain.It uses the spectral relationships between multiple-scattering aerosol reflectance and the different wavelengths expressed by polynomial functions, whereas τ a changes to M i , θ v , θ s , and φ sv , according to the following equation: where D is the degree of the polynomial (Table 2) and ρ M am (M i , λ) is the theoretically computed ρ am (λ) for model M i , considering the geometries and band pairs.The term c n represents the constants of the polynomial equation for models M i , θ s , θ v , and φ sv .To estimate ρ am (VIS) using A2016, as explained above, ρ M am (M i , NIR S ) for all candidate aerosol models can be calculated using Equation (26).Then, the most similar models, M L and M H , can be determined by comparison of the observed ρ am (NIR S ) and ρ M am (M i , NIR S ) of all candidate models, as follows: Therefore, two aerosol models, M L and M H , contribute to ρ am (NIR), as described below: The reflectance fraction at NIR (w M L and w M H ) can be directly calculated by solving the quadratic formula without any residual errors [20].
Extending Equation ( 26) by considering the wavelength-dependent reflectance fraction (w M L and w M H ), ρ am (VIS) can be derived as follows:

Simulation Dataset for the Evaluation
The four schemes for aerosol multiple-scattering reflectance estimation were evaluated using a simulation dataset generated by a vector radiative transfer code [22][23][24].Simulations were carried out for three aerosol models, excluding candidate aerosol models, and 24 scan geometries: θ s = 0 • , 25  , and 60 • ; and φ sv = 60 • and 120 • .The ρ wn for VIS wavelengths was modeled [29] for chlorophyll-a concentrations of 0.1, 0.3, and 1.0 mg/m 3 ; however, the NIR ρ wn was assumed to be negligible, to satisfy the black pixel assumption [30].Three aerosol models based on the work of Shettle and Fenn [31] were used for the evaluation: (1) the maritime model, with an RH of 80% (M80) that represents a coarse-size distribution of aerosol particles; (2) the coastal model, with an RH of 80% (C80); and (3) the tropospheric model, with an RH of 90% (T90) to represent a fine-scale aerosol particle distribution [31].The input parameters for the simulations are summarized in Table 3.We excluded validation data when ρ am (865 nm) was greater than 0.027, based on the cloud masking threshold of the SeaWiFS Data Analysis System (SeaDAS) [32].To avoid uncertainties arising from the sun-glint effect, cases in which the sun-glint reflectance at the surface exceeded 0.001 were removed.

Results and Discussion
This section describes the atmospheric correction results and intermediate aerosol parameters over Case-1 waters for the three aerosol models specified, after integrating the four atmospheric correction schemes.For quantitative analysis of the four schemes, we used the following statistical parameters: error (∆), relative percentage error (RPE), mean absolute percentage error (MAPE), median absolute percentage error (Med.), and root mean square error (RMSE), as follows: where K is the total number of matched pairs, and v n t and v n e are the true and derived values of the nth matched entry, respectively.The analysis was performed after replacing the four aerosol reflectance correction schemes within the atmospheric correction algorithm.The following nine aerosol models were used as candidates for the evaluation: oceanic model with RH 99% (O99), maritime model with RH 50, 70, 90, and 95% (M50, M70, M90, and M95, respectively), coastal model with RH 50 and 70% (C50 and C70, respectively), and tropospheric model with RH 50 and 80% (T50 and T80, respectively) [31].We used the same radiative transfer code [22][23][24] for both the evaluation data and candidate aerosol look-up tables.
The results are summarized as box-and-whisker plots with median, minimum, maximum, and quartile values (Figures 1 and 2). Figure 1a-c shows the errors (∆) in aerosol reflectance retrieval at 443, 555, and 660 nm, respectively, for the four schemes.Figure 1d,e show the errors (∆) in τ a at 865 nm and the Ångström exponent for 443 nm relative to 865 nm, respectively.Figure 2a-c show the RPE in ρ wn retrieval at 443, 555, and 660 nm, respectively.Table 4 summarizes the statistical results for ρ wn retrieval, for the four schemes and three aerosol types.Note that A2016 uses Equation (10) for retrieval of the τ a and Ångström exponent, because it does not calculate τ a itself.
parameters: error (Δ), relative percentage error (RPE), mean absolute percentage error (MAPE), median absolute percentage error (Med.), and root mean square error (RMSE), as follows: RPE (%) 100 ( ) where K is the total number of matched pairs, and vn t and vn e are the true and derived values of the nth matched entry, respectively.The analysis was performed after replacing the four aerosol reflectance correction schemes within the atmospheric correction algorithm.The following nine aerosol models were used as candidates for the evaluation: oceanic model with RH 99% (O99), maritime model with RH 50, 70, 90, and 95% (M50, M70, M90, and M95, respectively), coastal model with RH 50 and 70% (C50 and C70, respectively), and tropospheric model with RH 50 and 80% (T50 and T80, respectively) [31].We used the same radiative transfer code [22][23][24] for both the evaluation data and candidate aerosol look-up tables.
The results are summarized as box-and-whisker plots with median, minimum, maximum, and quartile values (Figures 1 and 2). Figure 1a-c shows the errors (Δ) in aerosol reflectance retrieval at 443, 555, and 660 nm, respectively, for the four schemes.Figure 1d,e show the errors (Δ) in τa at 865 nm and the Å ngström exponent for 443 nm relative to 865 nm, respectively.Figure 2a-c show the RPE in ρwn retrieval at 443, 555, and 660 nm, respectively.Table 4 summarizes the statistical results for ρwn retrieval, for the four schemes and three aerosol types.Note that A2016 uses Equation ( 10) for retrieval of the τa and Å ngström exponent, because it does not calculate τa itself.The median value of ∆ρ am (443 nm) for the four schemes falls between −0.0011~+0.0007(Figure 1).The GW1994, F1998, and AM1999 schemes tend to produce more errors, with underestimation of ρ am retrieval for T90 aerosol models compared to the other models due to the underestimation of the aerosol optical thickness and Ångström exponent.In the A2016 scheme, on the contrary, ρ am retrieval errors are more significant, with overestimation for M80 due to the overestimation of aerosol optical thickness and Ångström exponent.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 14 The median value of Δρam(443 nm) for the four schemes falls between −0.0011~+0.0007(Figure 1).The GW1994, F1998, and AM1999 schemes tend to produce more errors, with underestimation of ρam retrieval for T90 aerosol models compared to the other models due to the underestimation of the aerosol optical thickness and Å ngström exponent.In the A2016 scheme, on the contrary, ρam retrieval errors are more significant, with overestimation for M80 due to the overestimation of aerosol optical thickness and Å ngström exponent.The ρwn estimation results showed that the atmospheric correction accuracy for the four schemes was acceptable, with median APE values of ρwn of <5.1% and <3.5% at 443 nm and 555 nm, respectively.Although the RMSE of ρwn at 660 nm was smaller than that of ρwn at 443 nm, the MAPE of ρwn(660 nm) had a higher value than that of ρwn(443 nm), as the ρwn to aerosol reflectance ratio at 660 nm is significantly lower than that at 443 nm, due to relatively smaller water reflectance to aerosol reflectance ratio by relatively stronger water absorption at 660 nm.Similar to the aerosol estimation results, the accuracy of the ρwn estimation for the three schemes GW1994, F1998, and AM1999 improved for the coarse-sized aerosol models, in which aerosol reflectance had less of a multiplescattering effect.The A2016 scheme showed relatively similar accuracy to the three aerosol models, in which the MAPE ranged from 2.7% to 4.1%.The ρwn estimation results showed that the atmospheric correction accuracy for the four schemes was acceptable, with median APE values of ρwn of <5.1% and <3.5% at 443 nm and 555 nm, respectively.Although the RMSE of ρwn at 660 nm was smaller than that of ρwn at 443 nm, the MAPE of ρwn(660 nm) had a higher value than that of ρwn(443 nm), as the ρwn to aerosol reflectance ratio at 660 nm is significantly lower than that at 443 nm, due to relatively smaller water reflectance to aerosol reflectance ratio by relatively stronger water absorption at 660 nm.Similar to the aerosol estimation results, the accuracy of the ρwn estimation for the three schemes GW1994, F1998, and AM1999 improved for the coarse-sized aerosol models, in which aerosol reflectance had less of a multiplescattering effect.The A2016 scheme showed relatively similar accuracy to the three aerosol models, in which the MAPE ranged from 2.7% to 4.1%.

40.0%, Color scale for mean absolute percentage error (MAPE).
The ρ wn estimation results showed that the atmospheric correction accuracy for the four schemes was acceptable, with median APE values of ρ wn of <5.1% and <3.5% at 443 nm and 555 nm, respectively.Although the RMSE of ρ wn at 660 nm was smaller than that of ρ wn at 443 nm, the MAPE of ρ wn (660 nm) had a higher value than that of ρ wn (443 nm), as the ρ wn to aerosol reflectance ratio at 660 nm is significantly lower than that at 443 nm, due to relatively smaller water reflectance to aerosol reflectance ratio by relatively stronger water absorption at 660 nm.Similar to the aerosol estimation results, the accuracy of the ρ wn estimation for the three schemes GW1994, F1998, and AM1999 improved for the coarse-sized aerosol models, in which aerosol reflectance had less of a multiple-scattering effect.The A2016 scheme showed relatively similar accuracy to the three aerosol models, in which the MAPE ranged from 2.7% to 4.1%.

Note and Summary
Two-NIR-band-based aerosol reflectance correction schemes have been widely employed.We implemented and tested four such aerosol correction schemes: SeaWiFS/MODIS/VIIRS (GW1994), OCTS/GLI/SGLI (F1998), MERIS/OLCI (AM1999), and GOCI/GOCI-II (A2016).The GW1994 scheme determines the contribution of the aerosol models in the single-scattering domain after conversion from aerosol multiple-to single-scattering reflectance.F1998 determines the contribution to the aerosol optical thickness domain and then uses a weighted average approach to enhance the calculation accuracy of the weighting factors, whereas GW1994 uses an average value.AM1999 determines the contribution in the multiple-scattering domain and then uses linear distance to calculate the weighting factor.There are methods [28,34] that apply an approach similar to that of AM1999 to determine the contribution based on the linear distance of the reflectance in the multiple-scattering domain; these methods are expected to offer comparable performance to that of AM1999.A2016 is similar to AM1999 in terms of selecting aerosol models in the multiple-scattering domain, except that the weighting scheme for A2016 is enhanced to consider the wavelength-dependent weighting factor.
In this study, we intercompared four aerosol correction schemes in an assimilation.However, discrepancies between aerosol models and the actual aerosols can introduce more significant errors than the inaccuracies associated with the inverse scheme.Two-NIR aerosol correction schemes that we evaluated rely on the assumption of non-or less-absorbent aerosol optical properties, although actual aerosols originating from land can absorb more light.Alternative aerosol models are being developed that provide a more realistic representation of the optical properties of aerosols [35]; however, these models are also based on non-or less-absorbent models.All these schemes will significantly underestimate the ρ wn for strongly absorbing aerosols that contain a soot component, as discussed by the IOCCG study [21].To determine absorbance properties future aerosol correction schemes will require additional aerosol information from other wavelengths, i.e., from the near-ultraviolet regime [17] or from polarization [36].Errors in variables such as wind speed [37], trace gases [38][39][40][41], and air pressure [38] are higher at higher solar and sensor zenith angles.The atmospheric correction accuracy is also impacted by the system vicarious calibration [42][43][44].The radiance calibration requirement for atmospheric correction is <1%.However, the visible calibration gain can be varied by more than 1% due to different NIR intercalibration.Further work on the NIR calibration is needed to improve this.
The most recent atmospheric correction algorithms yield errors of less than 5% [2,45], which satisfies the accuracy requirement for ocean colour observations in the open ocean.However, more accurate atmospheric correction algorithms are required to observe diurnal changes from geostationary platforms such as GOCI, since geostationary measurements are made over a wide range of solar zenith angles which significantly increases the atmospheric correction uncertainty.Therefore, further investigation on the error sources from all relevant parameters is needed to improve atmospheric correction algorithms.

Figure 1 .
Figure 1.Errors in aerosol retrieval for the four schemes, presented as box-and-whisker plots.Aerosol reflectance differences by scheme for (a) the blue band (443 nm); (b) the green band (555 nm); and (c) the red band (660 nm); (d) aerosol optical thickness difference by scheme for the green band; and (e) aerosol Ångström exponent for 443 nm relative to 865 nm.

Figure 2 .
Figure 2. Water reflectance (ρ wn ) retrieval accuracy of the (a) blue; (b) green; and (c) red bands in the atmospheric correction algorithm after integrating the four aerosol correction schemes.Chlorophyll-a estimation error based on OC3 [33] is plotted in (d).

Figure 2 .
Figure 2. Water reflectance (ρwn) retrieval accuracy of the (a) blue; (b) green; and (c) red bands in the atmospheric correction algorithm after integrating the four aerosol correction schemes.Chlorophylla estimation error based on OC3 [33] is plotted in (d).

Table 1 .
Aerosol reflectance correction schemes based on the black pixel assumption.

Table 3 .
Summary of the input parameters for the simulations.

Table 4 .
Statistical summary of atmospheric correction accuracy.