The atmospheric correction algorithm employed by the NASA Ocean Biology Processing Group requires an assumption of negligible water-leaving reflectance in the near-infrared region of the spectrum. For waters where this assumption is not valid, an optical model is used to estimate near-infrared water-leaving reflectance. We describe this optical model as implemented for the sixth reprocessing of the SeaWiFS mission-long time-series (September 2009). Application of the optical model resulted in significant reductions in the number of negative water-leaving reflectance retrievals in turbid and optically complex waters, and improved agreement with in situ chlorophyll-a observations. The incidence of negative water-leaving reflectance retrievals at 412 nm was reduced by 40%, while negative reflectance at 490 nm was nearly eliminated.
©2010 Optical Society of America
The retrieval of accurate geophysical data products (e.g., spectral remote sensing reflectance Rrs(λ), and concentrations of the phytoplankton pigment chlorophyll-a, Ca) from satellite-borne radiometers requires the use of an ‘atmospheric correction’ algorithm . A critical component of this atmospheric correction algorithm is estimation of spectral aerosol reflectance, ρa(λ), which are subtracted from the total measured signal as part of the process of determining Rrs(λ). The algorithm implemented by the NASA Ocean Biology Processing Group (OBPG)  requires an assumption of negligible water-leaving reflectance in the near infrared (NIR) region of the spectrum (i.e., Rrs(NIR) = 0 sr−1). The NIR bands of the NASA Sea-viewing Wide Field-of-view Sensor (SeaWiFS) are centered on 765 and 865 nm. With this ‘black-pixel’ assumption, the measured top-of-atmosphere (TOA) reflectance in two NIR bands can be used to estimate both the magnitude and spectral dependence of ρa(NIR). Unfortunately, this assumption is rarely valid for waters with significant particulate (e.g., algal and mineral) backscattering . If Rrs(NIR) can be modeled, however, it can be removed from the TOA signal, allowing the atmospheric correction to proceed in estimating ρa(NIR) as if Rrs(NIR) were negligible.
The OBPG implemented a bio-optical model to estimate Rrs(NIR) as part of the third reprocessing of the SeaWiFS mission in May 2000 (R2000) [2,3]. This approach required an initial run through the atmospheric correction process (ignoring Rrs(NIR) > 0 sr−1) to estimate Ca. This preliminary Ca was used to estimate spectral particulate backscattering, bbp(NIR), which was used in turn to reconstruct Rrs(NIR). With this modeled Rrs(NIR) removed from the TOA signal, the atmospheric correction process was repeated, Ca was recalculated, the process was iterated upon until convergence in Ca was achieved [2,4]. With the fourth SeaWiFS reprocessing in July 2002 (R2002), this Ca-driven algorithm was replaced with a reflectance-based algorithm . The OBPG switched from  to  after observing that  depressed ratios of ρa(NIR) in optically complex waters. This often resulted in the selection of spectrally flat aerosol models within the atmospheric correction process, which artificially depressed the final Ca retrievals. The initial estimates of Ca were unrealistically high in these complex water masses, which led to over-estimates of bbp(NIR) and, subsequently, Rrs(NIR). The reflectance-based approach  was less prone to such over-correction. In this revised model, Rrs(λ) from an initial run through the atmospheric correction process and empirical estimates of the absorption coefficients for particles and dissolved materials at 670-nm were used to estimate bbp(670) via an ocean surface reflectance model . A spectral scattering function  was then used to derive bbp(NIR) from bbp(670). Finally, Rrs(NIR) was reconstructed from bbp(NIR), the atmospheric correction process was repeated, and the process was iterated upon until convergence in Rrs(NIR) was achieved.
Despite such efforts to account for non-negligible Rrs(NIR), the atmospheric correction algorithm continued to demonstrate sub-optimal performance in highly scattering waters, as manifest by negative Rrs(λ), particularly in the blue region of the spectrum, and erroneous spectral dependence in visible (VIS) part of the spectrum, with the latter leading to inaccurate Ca retrievals. To address these issues, a re-evaluation of the implementation of  was initiated prior to the sixth SeaWiFS reprocessing in September 2009 (R2009) and several modifications were made to the Rrs(NIR) model, which are presented below. First, an inconsistency in the relationship between Rrs(λ) and the inherent optical properties (IOP) was eliminated. Second, alternative estimates of the spectral dependence of bbp(λ) and of the absorption coefficients for particles and dissolved materials at 670 nm were adopted. Finally, a revised iteration scheme was implemented.
2. NIR reflectance model development
2.1 Inversion of remote-sensing reflectance
Remote sensing reflectance is proportional to the IOPs, specifically the marine absorption and backscattering coefficients, through:
G(λ) is a function of the illumination conditions, water constituents, and sea-state, bb(λ) is the total backscattering coefficient [ bbw(λ) + bbp(λ)], and a(λ) is the total absorption coefficient [ aw(λ) + ap(λ) + ag(λ)]. The additional subscripts w, p, and g indicate specific contributions by water, particles, and dissolved material, respectively.
The NIR reflectance model  implemented for R2002 and used through SeaWiFS reprocessing 5.2 in 2007 (R2007) suffered from an inconsistent application of Eq. (1). First, an approximation of this relationship, where G(λ) was fixed at 0.51 and X(λ) was simplified to bb(λ) / a(λ), was used to obtain near-infrared IOPs from Rrs(VIS). The quadratic form of , where G(λ) = [g1 + g2 X(λ)], g1 = 0.0949, and g2 = 0.0794, was later used to obtain Rrs(NIR) from the modeled near-infrared IOPs. While a minor inconsistency, using dissimilar definitions for G(λ) and X(λ) was found to have significant effects in highly scattering waters.
Furthermore, the approach was inconsistent with the methods of  that are used elsewhere (and independently) in the atmospheric correction process to normalize Rrs(λ) for bi-directional effects. In , G(λ) = f(λ) / Q(λ), where f(λ) is a function of the illumination conditions, water constituents, and sea-state, Q(λ) is the ratio of spectral upwelling irradiance to spectral upwelling radiance, and ratios of f(λ) to Q(λ) are pre-tabulated from extensive radiative transfer simulations. For the Rrs(NIR) model implemented in R2009, G(λ) = f(λ) / Q(λ), using the look up tables (LUTs) of , for both the conversion of retrieved Rrs(VIS) to IOPs and the conversion of modeled IOPs to Rrs(NIR). While these LUTs were originally generated for waters whose light field is predominantly influenced by phytoplankton, they proved suitable for use in this model, as will be demonstrated in Section 4. A full discussion of the applicability of  to all water types is beyond the scope of this paper.
2.2 Spectral backscattering
Given the relationship in Eq. (1), a measured Rrs(λ0), and knowledge of a(λ0), an estimate of bb(λ0) can be obtained, where λ0 is any reference wavelength (e.g., 670 nm). This retrieved bb(λ0) can then be extrapolated to the NIR through a spectral model. The spectral dependence of bb(λ) can be expressed as:
In Eq. (2)b, the η term is typically in the range of [0, 2], with values around 2 observed in open ocean waters and turbid waters approaching zero . In contrast, the model employed in R2002-R2007  used a spectral scattering relationship for Eq. (2)a developed by :Eq. (2)b was adopted. Rather than relying on a fixed η value, however, the empirical relationship developed in  was adopted:
The NASA bio-Optical Marine Algorithm Data set (NOMAD)  includes measurements of bb(VIS) at a number of wavelengths spanning a range of values from 0.0007 to 0.0133 m−1. Unfortunately, this data set does not include measured bb(NIR). Nevertheless, this data set was used to evaluate the various forms of Eq. (2) and 3 using a reference wavelength of 670 nm and a retrieval wavelength of 555 nm. The formulation of  had the minimal bias and median percent difference relative to in situ measurements (Table 1 , Fig. 1 ). and performed better than the linear, spectral scattering function .
2.3 Estimating absorption
In the NIR, a(λ) is dominated by the absorption due to water, aw(λ) [10,11], such that absorption coefficients for particles, ap(NIR), and dissolved material, ag(NIR), can be ignored. However, a reasonable estimation of the absorption coefficient is crucial for the inversion of Rrs(λ) to retrieve bb(670), where ap(670) and ag(670) can be a significant portion of the total. In , the estimate of a(670) involved two functions, one based on Ca  and the other on a green-to-red reflectance ratio . In the modified Rrs(NIR) model, this was updated to a single, Ca-based relationship derived from NOMAD  (Type II linear regression using spectroscopy only, limited to Ca > 0.2 mg m−3; Fig. 2 ):
3. Iteration scheme
Since the retrieval of Rrs(λ) in highly scattering waters requires knowledge of Rrs(NIR), and estimation of Rrs(NIR) requires knowledge of Rrs(VIS), the Rrs(NIR) algorithm must be implemented iteratively. The iteration begins by assuming Rrs(NIR) = 0 sr−1, so that initial Rrs(λ) can be derived and a first estimate of Rrs(NIR) can be achieved. The iteration converges when Rrs(λ) changes by less than 2%. In the majority of cases, the method converges after 3-4 iterations, but the algorithm allows up to 10 iterations.
If, however, the initial atmospheric correction results in negative or an otherwise non-physical spectral distribution for Rrs(λ) (criteria based on in situ observations), the iteration is re-initialized assuming ρa(NIR) = 0 (i.e., all reflectance except that from Rayleigh scattering and Sun glint originates from the water mass). This is effectively the opposite extreme from the initial condition of Rrs(NIR) = 0 sr−1. The iteration process is then allowed to proceed to convergence as previously described. For the rare cases where the iteration fails to converge after this reset, an additional iteration is forced, again with an assumption of ρa(NIR) = 0. This pixel is then flagged with an ‘atmospheric correction warning’ to alert users to the questionable nature of the pixel, although a qualitatively useful retrieval may still result.
The controls on this iteration scheme differ slightly from previous implementations , especially with regards to the how the iteration is reset, which is currently more robust than that of past methods. It should also be noted that the Rrs(NIR) is forced to zero if the initial iteration results in a Ca < 0.3 mg m−3. This is done to avoid the addition of noise into clear-water pixels through the estimation of otherwise negligible Rrs(NIR) contributions. Furthermore, to avoid the introduction of transitional artifacts, the Rrs(NIR) estimation is linearly weighted from 0 to 1 for 0.3 ≤ Ca ≤ 0.7 mg m−3 (Fig. 3 ). This Ca-based weighting of Rrs(NIR) was first introduced in R2002 and is effectively unchanged for R2009.
A long-term satellite and in situ data set was recently compiled to study Ca variability in Chesapeake Bay . These data were used to evaluate the modifications to the NIR correction described above, as implemented for the most recent SeaWiFS reprocessing. Specifically, the NIR correction implemented with R2002 (and used through R2007) was compared to the method described in this work and employed for R2009. The satellite and in situ data acquisition and processing methods are described in .
The revised correction reduced the percentage of negative Rrs(412) for the time series of SeaWiFS retrievals over the Bay by nearly half (18.5% vs. 31.9%). The incidence of negative Rrs(490) was reduced to negligible levels (Table 2 ). The satellite estimates of Ca showed a substantial improvement in the agreement with ground truth measurements in the turbid and highly productive regions of the Middle and Upper Bay (Fig. 4 ).The revised correction also improved reflectance retrievals relative to coincident in situ measurements. Following the method of , and using a global validation data set compiled by the OBPG, we compared Rrs(λ) retrievals for SeaWiFS after limiting the data set to where the NIR correction was used in the processing of the satellite data (e.g., Ca ≥ 0.3 mg m−3). Improvement in the satellite agreement with in situ Rrs was evident for all bands, with the greatest improvement seen for the 412 and 443 nm bands (Table 3 ).
The revised Rrs(NIR) model that is presented in this paper has been implemented in the operational processing of satellite ocean color sensor data by the NASA Ocean Biology Processing Group. This revised correction eliminated a number of inconsistencies present in the previous version. A variable estimate for the spectral slope of the backscattering coefficient was implemented in place of a fixed linear relationship. A simplified estimate of absorption in the red region of the visible spectrum was also implemented. This resulted in a significant reduction in the incidence of negative Rrs for satellite retrievals in the 412 to 490 nm range, while improving the agreement of satellite-derived Rrs(λ) and Ca with in situ measurements.
The results and discussion presented here have been focused on ocean color retrievals from the SeaWiFS sensor. It should be noted, however, that the same algorithm has been applied in the processing of MODIS-Aqua and other ocean color sensor data sets archived by NASA’s Ocean Biology Processing Group.
The approach presented, while a significant improvement over the previous implementation, cannot correct all cases where the black pixel assumption is invalid. The use of empirical optical models (e.g. spectral backscattering, particulate and dissolved absorption) limits the applicability of this algorithm to waters that are similar to those over which these empirical models were developed. In addition, the practical implementation of the iteration scheme can result in conditions of non-convergence. Fortunately, the optical domains covered by the data sets used in the generation of these algorithms are quite diverse and the cases of convergence failure are small. Therefore, these limitations do not detract from the benefits realized by the use of this algorithm.
References and Links
1. H. R. Gordon and M. Wang, “Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS: a preliminary algorithm,” Appl. Opt. 33(3), 443–452 (1994). [CrossRef] [PubMed]
2. D. A. Siegel, M. Wang, S. Maritorena, and W. Robinson, “Atmospheric correction of satellite ocean color imagery: the black pixel assumption,” Appl. Opt. 39(21), 3582–3591 (2000). [CrossRef]
3. C. R. McClain, E. Ainsworth, R. Barnes, R. Eplee, Jr., F. Patt, W. Robinson, M. Wang, and S. W. Bailey, “SeaWiFS Postlaunch Calibration and Validation Analyses, Part 1,” NASA Tech. Memo. 206892, National Aeronautics and Space Administration, Goddard Space Flight Center, Greenbelt, MD (2000).
4. A. Morel, “Optical modeling of the upper ocean in relation to its biogenous matter content (case 1 waters),” J. Geophys. Res. 93(C9), 10749–10768 (1988). [CrossRef]
5. R. P. Stumpf, R. A. Arnone, J. R. W. Gould, P. M. Martinolich, and V. Ransibrahmanakul, “A partially coupled ocean-atmosphere model for retrieval of water-leaving radiance from SeaWiFS in coastal waters,” in Patt, F.S., et al., 2003: Algorithm Updates for the Fourth SeaWiFS Data Reprocessing. NASA Tech. Memo. 206892, National Aeronautics and Space Administration, Goddard Space Flight Center, Greenbelt, MD (2003).
6. H. R. Gordon, O. B. Brown, R. H. Evans, J. W. Brown, R. C. Smith, K. S. Baker, and D. K. Clark, “A semianalytic radiance model of ocean color,” J. Geophys. Res. 93(D9), 10909–10924 (1988). [CrossRef]
7. R. W. Gould Jr, R. A. Arnone, and P. M. Martinolich, “Spectral dependence of the scattering coefficient in case 1 and case 2 waters,” Appl. Opt. 38(12), 2377–2383 (1999). [CrossRef]
8. NIR, html, NIR Correction, http://oceancolor.gsfc.nasa.gov/REPROCESSING/SeaWiFS/R4/NIR.html
9. A. Morel, D. Antoine, and B. Gentili, “Bidirectional reflectance of oceanic waters: accounting for Raman emission and varying particle scattering phase function,” Appl. Opt. 41(30), 6289–6306 (2002). [CrossRef] [PubMed]
11. R. M. Pope and E. S. Fry, “Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements,” Appl. Opt. 36(33), 8710–8723 (1997). [CrossRef]
12. Z. Lee, R. A. Arnone, C. Hu, P. J. Werdell, and B. Lubac, “Uncertainties of optical parameters and their propagations in an analytical ocean color inversion algorithm,” Appl. Opt. 49(3), 369–381 (2010). [CrossRef] [PubMed]
13. P. J. Werdell and S. W. Bailey, “An improved in-situ bio-optical data set for ocean color algorithm development and satellite data product validation,” Remote Sens. Environ. 98(1), 122–140 (2005). [CrossRef]
14. A. Bricaud, A. Morel, M. Babin, K. Allali, and H. Claustre, “Variations of light absorption by suspended particles with the chlorophyll a concentration in oceanic (case 1) waters: analysis and implications for bio-optical models,” J. Geophys. Res. 103(C13), 31033–31044 (1998). [CrossRef]
15. P. J. Werdell, S. W. Bailey, B. A. Franz, L. W. Harding Jr, G. C. Feldman, and C. R. McClain, “Regional and seasonal variability of chlorophyll-a in Chesapeake Bay as observed by SeaWiFS and MODIS- Aqua,” Remote Sens. Environ. 113(6), 1319–1330 (2009). [CrossRef]
16. S. W. Bailey and P. J. Werdell, “A multi-sensor approach for the on- orbit validation of ocean color satellite data products,” Remote Sens. Environ. 102(1-2), 12–23 (2006). [CrossRef]