1. Introduction
The literature on the optics of natural waters before 1970 suggested that the diffuse reflectance of the water contained information regarding the water itself and its constituents, i.e., mineral sediments as well as dissolved and suspended organic material (e.g., see Jerlov [
1,
2]). However, it was the pioneering work of Clarke et al. [
3] that demonstrated that the upward spectral radiance measured above the sea contained information related to the concentration of chlorophyll
a in the water. Chlorophyll
a is the photosynthetic pigment contained in phytoplankton, which constitute the first link in the marine food chain. This direct observation led to the development of several aircraft imaging radiometers for viewing coastal (and ocean) waters and then to the first space-borne imager specifically designed for viewing the coasts (and ocean): the Coastal Zone Color Scanner (CZCS). Clarke et al. [
3] also demonstrated that the information was degraded by the backscattering of the atmosphere between the surface and the sensor which they referred to as “air light”. Regarding this “air light”, at the end of their paper, they commented [
3] “If such interference can be eliminated or identified and allowed for, spectroscopic procedures from aircraft (and perhaps from satellites) will be of great value in the rapid investigation of oceanic conditions, including conditions important for biological productivity.” The effect of air light was further examined by Hovis and Leung [
4] from altitudes as high as 16 km, where they observed that as much as 90% of the measured radiance actually resulted from scattering of sunlight between the surface and the sensor. This scattered light contained no information regarding the water. The removal of the air light radiance and the estimation of the radiance backscattered
from within the water—the “water leaving radiance”—from that measured by the sensor is usually referred to as “atmospheric correction”.
In this historical review I will describe the events leading up to the basic algorithm used for atmospheric correction of the modern U.S. sensors that followed the prototype CZCS, namely SeaWiFS, MODIS, and VIIRS (Gordon and Wang [
5]). I devote much of the paper to this particular algorithm because of its widespread operational use, and refer to it as “GW94”. I compare GW94 to the OCTS and MERIS algorithms, showing that they are essentially the same. I also show how (and why) the modifications to the algorithm were effected as a new understanding of the nature of aerosols over the oceans was obtained. A constraint on the original algorithm was that the aerosol exhibited only weak absorption compared to its scattering, and it often fails when strongly absorbing aerosols are present in the atmosphere, e.g., Saharan dust. Methods of dealing with strongly absorbing aerosols are described as well. Noting that the algorithm was developed for Case 1 waters, i.e., waters for which the optical properties all more-or-less covary with the concentration of chlorophyll
a [
6,
7], and it often fails when applied to Case 2 waters (e.g., coastal waters), I will examine some efforts to extend the algorithm into the coastal regime; however, as this special issue of
Remote Sensing is centered on atmospheric correction over such waters, the examination will be brief.
The review is divided into time periods that are quite naturally specific to individual sensors: 1970–1980, CZCS development; 1980–1990; CZCS operation; 1990–2000, SeaWiFS/MODIS preparation (and operation); and 2000 to roughly 2005 (near the end of my active involvement in atmospheric correction development), algorithm enhancements and improvements. However, some work started before, but published after 2005, is described for completeness.
It is natural that in outlining the evolution of events in which one has participated, the result will reflect the author’s perception of their import at the time. As such, there are likely contributions to the subject that, while important in their own right, had little impact on the development of the correction algorithm. I apologize for their omission.
2. 1970–1980: Algorithm Conception and CZCS Development
It is important to note that up to and through the time period in which the CZCS was designed and fabricated there were no actual algorithms relating the color of the water (i.e., the spectrum of upwelling radiance just beneath the surface) to water constituents. The first such algorithm was that of Morel and Prieur [
7] in which the ratio of the subsurface diffuse reflectance at 440 to that at 550 nm was related experimentally and theoretically to the “pigment concentration” (
, the concentration of chlorophyll
a and its degradation product phaeophytin
a) in the water. This was closely followed by similar algorithms, based on a data set from coastal and ocean waters around the U.S. [
8,
9]. In addition, it is clear that in the design of CZCS little or no attention was paid to the problem of removing, or in some way accounting for, the masking effects of the air light. After its formation, the CZCS Experiment Team (
Table 1) was in effect handed the specifications for the instrument and asked to “do what you can with this to make it work”.
The work of Clarke, et al. [
3] and Hovis and Leung [
4] had clearly shown that the color of the water, i.e., the spectral radiance
exiting the water, was significantly modified, when observed from high altitude, by scattering and absorption by the atmosphere between the surface and the sensor. Thus, one of the principal tasks set forth by the CZCS Experiment Team was to develop a method to account for the atmospheric effects. The other was to establish a quantitative relationship between the color of the water and the concentration of chlorophyll
a. I volunteered to take the lead in the former.
Believing that aerosols presented the most significant problem for atmospheric correction, and applying the adage “one man’s signal is another man’s noise”, I started by looking at efforts to retrieve aerosol concentrations from satellite-measured radiance, principally from the LANDSAT MSS sensor. Griggs [
10] compared LANDSAT (then ERTS-1) radiances to the aerosol optical depth
(at 550 nm) over water surfaces and found a near-linear relationship between the two for the small range of
observed in the study. This suggested that given the aerosol optical depth (and possibly other parameters associated with the aerosol) and the
known scattering by the air itself, it might be a simple matter to estimate the radiance at the sensor. Certainly this should be possible if a full multiple scattering solution to the radiative transfer equation were developed; however, with the enormous amount of data from the CZCS (one CZCS two-minute “scene” contained approximately 1.6
pixels with 4 usable wavelengths), full radiative transfer calculations at each pixel were clearly impossible. In addition, the aerosol properties required for such calculations are not known
a priori. Thus, a simpler method was required to determine the atmospheric effect. The fact that the radiances for LANDSAT seemed to relate linearly to aerosol optical depth suggested that single scattering might be a usable approximation for the radiative transfer. Briefly, the (scalar) radiative transfer equation for the radiance
L in a plane parallel medium illuminated from the top (
) by a beam of irradiance
propagating in a direction
is
where
is a unit vector directed
into the medium at the top,
is a unit vector in the direction of propagation of the radiance,“•” indicates the scalar product,
is the optical depth,
is the wavelength,
is the single-scattering albedo (scattering coefficient ÷ extinction coefficient),
P is the scattering phase function, and
is a small solid angle around the direction
. The usual single-scattering solution to this equation involves simply leaving out the integral term and solving the remaining inhomogeneous differential equation. If the atmosphere is modeled as a plane parallel homogeneous slab of thickness
, the radiance exiting the top of the atmosphere (when no radiance enters either the top or bottom of the atmosphere) is
where
and
(note that
for radiance
exiting the top of the slab), and
and
are the azimuth angles for
and
, respectively. This is often referred to as the single-scattering radiance reflected from the slab; however, as it contains terms of order higher than the first in
, i.e., the exponential term, it is not really single scattering. True single scattering is obtained by expanding the exponential and retaining the lowest order term in
. The result is
Often it is desirable to work with a normalized radiance
(referred to as “reflectance” here) rather than radiance. The reflectance is defined through
. (Note that
is dimensionless; however, the defining equation yields units of inverse solid angle,
. Since a solid angle is just the ratio of two areas, as a plane angle is the ratio of two lengths, and hence dimensionless, I prefer to leave
unit less; however, if the reader is uncomfortable with my convention, then the unit
should be attached to all reflectances.) Then,
If we modify the model to account for the presence of a flat specularly reflecting water surface at
(but ignore the radiance
exiting the water from below), it is easy to show that Equation (
3) becomes
where
, and
is the Fresnel reflectance of the air-water interface for an incident angle (in the air) of
. How accurate is this approximation?
Figure 1 shows the reflectance of a slab with Rayleigh scattering characteristic of Earth’s atmosphere at 412 nm (Rayleigh optical thickness of 0.312) as a function of both the solar zenith angle
and the exiting angle
in the perpendicular plane (
).
The accuracy of Equation (
4) in reproducing the results of exact (scalar) theory is truly remarkable. It suggested that this equation could be used to estimate the Rayleigh-scattering contribution to the radiance at the top of the atmosphere with an error of no more than about
(note that the blue band on CZCS has
nm, where the optical thickness is about 30% lower than that in
Figure 1).
Figure 2 shows a similar computation when aerosols are included in the slab. In this case, the aerosols are characteristic of a maritime aerosol with optical thickness of 0.2 at 412 nm. The computation was effected by computing the contributions from Rayleigh scattering (
and aerosol scattering (
) separately and then adding them (
).
The error is now larger than in
Figure 1, but so is the total optical thickness (>0.5). The larger error for small
and small
owes to the large forward scattering by the aerosol, which is reflected from the lower surface. This suffers too little attenuation upon propagating to the upper boundary (zero attenuation in fact). In reality, this particular error in radiance is not very important as the Sun’s glitter pattern for a roughened (by the wind) surface would likely overwhelm the radiance reflected from the atmosphere. Thus, it appeared that, within an error of only a few percent, one could write
for the normalized radiance exiting the atmosphere, where
and
are computed individually using Equation (
4). Note, this agrees with the observation [
10] that LANDSAT radiances varied linearly with aerosol concentration (i.e., aerosol optical depth). It was taken as the working equation for the analysis of CZCS data [
12], i.e.,
where
is the total reflectance measured at the sensor,
is the radiance exiting the water converted to reflectance and
t is its diffuse transmittance [
13] (see below).
It is easy to underestimate the importance of the observation that single-scattering could be used to provide a reasonably accurate estimate of the radiance at the sensor that was due to the atmosphere. One need only consider the computing power that was available at that time. When CZCS was launched (November 1978) 1 MHz workstations with 1MB of memory were about a decade in the future. In fact some of the track line comparisons of ship measurements and satellite estimates of the concentration of chlorophyll we published (e.g., Gordon et al. [
14]) were carried out by the author on a desktop computer with a single line display and 64 kB of memory (Hewlett Packard 9825). Thus, any complex computation of the atmosphere’s influence on the sensor’s radiance was out of the question. Mainframe computers at that time were too expensive to be dedicated to image processing. Actually, it was not until the early to mid 80s that even the single-scattering computation of
could be effected on a pixel-by-pixel basis in an image processing environment (this application employed supermini computers such as the VAX 11/780).
There are many parameters in Equations (
4) and (
5) that were unknown in the application:
,
P, and
for the aerosol in Equation (
4); and
,
t, and
in Equation (
5). These had to be estimated/determined from the imagery itself. In order to make up for the unknown quantities, several assumptions were required. The CZCS had four bands with sufficient sensitivity to be useful for ocean color. Their band centers were at 443, 520, 550, and 670 nm. The first assumption was that
, the “water-leaving reflectance” in the red band, was negligible, i.e.,
. This is a good assumption for Case 1 waters with chlorophyll concentrations <0.5–1 mg/m
3 (the oceanic chlorophyll concentration has a broad distribution from 0.06 to about 0.3 mg/m
3 [
15] with a global median around 0.13 mg/m
3 and higher values at higher latitudes). Then Equation (
5) provided
, as
could be readily computed from Equation (
4). Note that in Equation (
4) the angles are independent of wavelength, so that only
,
P, and
vary with
. It was then assumed that
and
P for the aerosol are independent of wavelength, so only the aerosol optical thickness
was spectrally dependent. Since
is directly proportional to
, the aerosol optical depth, one has
which shows that the water contribution to the total reflectance can be retrieved from
Note that a term is included on the right-hand side above even though it was assumed . This is for later use in situations in which .
As
can be computed directly and
, the only unknown on the right hand side of Equation (
7) is the spectral variation of
. For many environmental aerosols
where the (usually) positive number
n is called the Ängstrom exponent. For a maritime aerosol,
. Thus, at this point all that was required to find
in the blue and green bands was
n. In more modern notation, the multiplier of
is called
, and Equation (
7) is usually written
and for completeness, in the single-scattering approximation
is given by
where
,
P, and
refer to the associated
aerosol properties.
Once
is determined, it remained to estimate
t. In the absence of scattering in the atmosphere, the radiance that exited the water from a pixel being viewed by the sensor can arrive at the sensor by a single path: precisely in the direction of the sensor. This radiance is attenuated by a factor
. However, in the presence of atmospheric scattering radiance leaving the water in a different direction (i.e., from a
different pixel) can also reach the sensor by being scattered into it. The later radiance depends on the angular distribution of radiance leaving the surface. If this radiance is totally diffuse (independent of the exiting direction) and the water is horizontally homogeneous, then for a pure Rayleigh-scattering atmosphere, and in the single-scattering approximation, the diffuse transmittance is given by (see [
11])
where the “*” is to remind us that this formula is restricted to horizontal homogeneity of the water
and to a totally diffuse
at the surface. (The water-leaving radiance distribution is far from diffuse; however, the upward radiance
just below the surface of the water is only weakly dependent on direction. My preference then was to compute
assuming the upwelling radiance
beneath the surface is totally diffuse, as in [
13,
14].) The dependence of
on
owes to photons exiting the surface, backscattering toward the surface, and then reflecting from the surface in the direction of the sensor. The aerosol contribution to
t is much smaller than the Rayleigh contribution because of the strong forward scattering exhibited by aerosols. It usually results in <1% change in
t and was ignored. Thus, letting
, all that remained for determining
from
was the determination of
.
In the initial validation of the CZCS-derived pigment concentration [
16],
was determined by direct ship-based measurements of
at the surface (i.e., at one pixel in the image). The derived
-values were then used throughout the entire image. The uniformity of
was justified by examination of the image provided in
Figure 3. The image of
was derived by using an estimate for
that was valid for the low-chlorophyll concentrations found in the interior of the Gulf of Mexico at the time of the image, Fall, 1978 (see later discussions of
as a function of the chlorophyll concentration).
Figure 4 provides a histogram of the resulting values of
and
for the pixels within the rectangle drawn on
Figure 3. Note that
varies by
while
varies by about
, and some of the variation in
is due to variation in
, albeit small. Thus, it was taken as a working hypothesis at that time that
was constant over an image (or sub-image), even in the presence of variations in the aerosol concentration. This is nearly equivalent to assuming that the aerosol size
frequency distribution and refractive index are uniform over an image, even when the aerosol concentration varies (
nearly because
will undergo some variation because of its dependence on the Sun-viewing geometry).
In sum, for atmospheric correction of the initial imagery from CZCS [
16,
17], the principle assumptions were: (1) the single-scattering approximation was sufficient to estimate
; (2)
; (3) the parameter
was approximately constant over a CZCS scene; and (4) the upward radiance just beneath the water surface was totally diffuse.
Ship measurement of
contemporaneous with the satellite overpass was carried out to determine
, and these values were used for the entire image. The results suggested that the pigment concentration (
) calculated using the radiance ratios retrieved from the CZCS, i.e.,
was valid to within a factor of two. It should be noted that these relationships for
[
8,
9,
16] were developed from the sparce in situ measurements without regard to the Case 1–Case 2 distinction due to the paucity of surface measurements at that time. It is also important to note that within the literature, relationships such as these were dynamic in the sense that the coefficients, i.e., 0.5 and 1.3 in the first, were subject to variation as new in situ data were obtained, and as stratification of in-situ data between Case 1 and Case 2 waters became possible.
4. 1990–2000: Preparation for Improved Instruments
Even before the launch of CZCS it was understood that, although the sensor was more-or-less optimized for observing the principal variations in ocean color, i.e., those resulting from variations in the concentration of chlorophyll-like phytoplankton pigments, it was not well designed for atmospheric correction [
39]. That is, there were no spectral bands for which the sensor recorded
only radiance resulting from atmospheric and surface interactions. In fact, the radiance exiting in the red band is negligible only in low to moderate chlorophyll waters, and hardly ever near the coasts. However, for bands in the near infrared (NIR, 700–1000 nm), it is usually the case that
, except in regions with high suspended particle load (non-living or living). The CZCS had a wide (100 nm) spectral band at 750 nm that had a radiometric sensitivity similar to a corresponding LANDSAT band, but that radiometric sensitivity was too low to be of value for atmospheric correction. The opportunity of increasing the sensitivity of this band electronically was presented to the author after CZCS had been integrated with the NIMBUS satellite; however, as this involved considerably risk, the decision was made not to effect it. A working group [
39] proposed in the summer of 1978 that ocean color sensors following CZCS have at least two spectral bands in the NIR with which the effect of the atmosphere on the radiance at the sensor could be assessed and extrapolated to the visible bands. Such spectral bands have been placed on all post-CZCS sensors, with the MODIS and VIIRS (and some proposed) sensors having additional bands that are still further into the infrared (the short-wave infrared, SWIR), although their presence was justified for other purposes.
4.1. Required Accuracy for Atmospheric Correction
In the case of CZCS, where there were no bio-optical algorithms to use as a guide, the goal was to perform atmospheric correction as best as one could (given the meager computational resources). The goal of the second-generation SeaWiFS sensor was stated (SeaWiFS Prelaunch Report V1 [
40]) “To achieve ... water-leaving radiances to within 5% absolute, and chlorophyll a concentration to within 35% over the range of 0.05–50.0 mg m
−3”. However, the description of this goal was incomplete. For example, at 670 nm
, and retrieving a very small radiance to within 5% is clearly impossible. The correctly stated goal was to retrieve the water-leaving radiance at 443 nm with an accuracy of 5%
in oligotrophic waters. But even this more complete description of the goal still depends on the definition of “oligotrophic.” The author took that to mean a pigment concentration of about 0.05 mg/m
3.
Figure 5 shows that
for this pigment concentration, so I took the goal of atmospheric correction to be
. Linear regression of the log-transformed data in
Figure 5 provides
and when
mg/m
3 the standard error of estimate in
, using the second form for regression, is about 27%. Using this we find
so for the above-stated goal of atmospheric correction the atmospherically-induced error would correspond to an error in
of approximately 18% at
mg/m
3. Combining this with the error in the bio-optical algorithm, in a root-sum-of-squares (RSS) addition of the two, predicts a total error in
of approximately 32% for
mg/m
3; however, for
mg/m
3 the atmospheric correction error approximately doubles and the combined RSS error increases to about 45%. In order to achieve an uncertainty of less than 35% in an RSS sense over the range
mg/m
3, one would need
. With the final SeaWiFS algorithm, this goal was achievable sometimes, but not always; however the goal of
appeared to be achievable under most circumstances as long as the aerosol is non-absorbing or weakly absorbing. Antoine and Morel [
41] obtained similar criteria using the requirement that a sensor (MERIS in their case) be capable of distinguishing 30 equally spaced classes of
(in their paper they used the chlorophyll
a concentration,
C, rather than
) in the range
mg/m
3.
4.2. An Extended Single-Scattering Algorithm
It was important to implement the equations for atmospheric correction, in the single scattering approximation, including these additional bands in the NIR as an initial substitute for the yet-to-be-developed multiple scattering algorithm, which was expected to be orders of magnitude more computer intensive.
Consider two spectral bands in the NIR positioned at
and
, where “
s” and “
ℓ” stand for “short” and “long,” respectively. For SeaWiFS,
nm and
nm. Thus, since
for
and
, the equivalent of Equation (
9) is
with
where
,
, and
refer to the associated
aerosol optical properties. Writing Equation (
23) for
, we have
which would be a known quantity. It should be noted that in Equations (
23) and (
25),
is provided in the form of look-up-tables (LUTS) prepared from
exact computations (including polarization) of the radiance exiting a
purely Rayleigh scattering atmosphere bounded by a
totally absorbing ocean with a flat Fresnel reflecting surface, i.e., the multiple-scattering analog to Equation (
4). In the case of SeaWiFS, an estimate of the atmospheric pressure was to be available and would modify
in a manner prescribed by Equation (
19).
The basic idea of the algorithm was to estimate
given
. In the case of a non-absorbing aerosol for which the phase function is independent of wavelength, this is particularly simple (in this single-scattering approximation):
or
, and the spectral variation of the aerosol optical depth alone determines the extrapolation into the visible. An example we considered earlier was when the aerosol satisfies Ängstrom’s law (
), in which case
determines
n, and then
. Unfortunately, the phase function of most aerosols depends on
, so this simple example (used here only as an illustration) is not realized in practice, although it was
universally assumed in the earlier algorithms (e.g., [
14,
22,
23]). Wang and Gordon [
42] studied the efficacy of Equations (
23)–(
25) when the realistic aerosol models developed by Shettle and Fenn [
36] were used to characterize the aerosol. They found that for these aerosol models
where
k is a constant that depends on the particular model and Sun-viewing geometry. It can be determined via
where
is found from Equation (
25). Using simulations, they concluded that (1) if the aerosol was non-absorbing (or weakly absorbing, e.g., maritime), (2) if
, (3) if the estimated
(
) and (4) if the resulting value of
is reduced by 4.6%, then
is usually less than 0.002. The required reduction in
by 4.6% is consistent with the approximate correction of Gordon and Castaño (1987) [
34] for the effects of multiple scattering. The extended single-scattering algorithm was to be a “fall-back” algorithm if its multiple-scattering counterpart proved too costly in terms of computational resources.
4.3. The SeaWiFS Multiple Scattering Atmospheric Correction Algorithm
In 1989, the author and M. Wang began working on an extension of the algorithm to include the effects of multiple scattering in a less ad hoc manner than Gordon and Castaño [
34]. They started with the Deschamps et al. [
33] observation that there was significant contribution to
from photons that scatter from both molecules and aerosol particles, i.e.,
in Equation (
20) or the replacement of
by
. To avoid confusion regarding what refers to single scattering and what refers to multiple scattering, we now change the notation slightly. Let the single-scattering reflectance exiting the atmosphere after scattering by the aerosol henceforth be indicated by
, where
with the subscript “
a” indicating “aerosol”, and subscript “
s” indicating “single scattering.” Note that, in a given Sun-viewing geometry,
is directly proportional to the aerosol optical thickness
. The success of the single-scattering algorithm at low values of
suggested that it should be mimicked as closely as possible in a multiple-scattering version. That was the path followed. The only link between the aerosol and its physical properties that can be determined from
are the values of
and
, both of which contain the effects of multiple scattering. Because
is proportional to
, it is clear that for a given wavelength, Sun-viewing geometry, and aerosol model (say the
jth), one can write
as a function of
. An example is
where
,
and
are constants that depend on the model and the Sun-viewing geometry. In the single-scattering approximation,
and
. For now we use the following notation for the relationship designated in Equation (
27),
where the big curly brackets mean
for the
jth model given as a function of
for the
jth model. (Later I will simplify this notation to
or
for the
jth model, where the
g’s are the polynomials in
.) Equation (
28) can be inverted to yield
as a function of
for the
jth model, i.e.,
The relationships in Equations (
28) and (
29) must be established by way of rigorous solutions to the radiative transfer equation. These in turn require models for the aerosol and the structure of the atmosphere. Gordon and Wang [
5] simplified the vertical structure of the atmosphere to just three layers: (1) next to the Fresnel-reflecting, flat, water surface, was a layer containing
only aerosols of a specific type and characterized by the optical thickness
; (2) a layer above the aerosol containing only Rayleigh-scattering molecules; and (3) a layer above the second displaying only absorption to account for Ozone (and which can be neglected in the radiative transfer computations). The bulk water itself was considered to be totally absorbing, so no radiance crosses the water surface from below. Calculations showed that, when the aerosol was non- or weakly absorbing, the radiance exiting the top of the atmosphere was insensitive to its vertical structure [
43,
44], so this atmospheric structure should yield radiances that are sufficently accurate to establish the relationships in Equations (
28) and (
29) as long as realistic aerosol models are employed.
The aerosol was modeled according to Shettle and Fenn [
36]. Their aerosol models consisted of two components: a log-normal size distribution of small particles (modal diameter < 1 μm) and a similar distribution of large particles (modal diameter > 1 μm). Over the oceans the small fraction represented either the background aerosol or an aerosol transported from land. In contrast, the large-sized component represented locally generated aerosol (e.g., by breaking waves). Each component swells with increasing humidity and its refractive index decreases (but always remains greater than that of pure water). The background aerosol has a refractive index (dry) of approximately 1.53 with a modest absorption, while the large-sized component has an index close to water and virtually no absorption. Gordon and Wang [
5] used the Maritime models in SF79 at four values of relative humidity (
, 70, 90, and 99%) and an additional model with
was used to test the final algorithm. These are referred to as M50, M70, M90, M99, and M80. They also used models without the large-sized component and referred to these as terrestrial or tropospheric aerosols with the label “T” (i.e., T50, T70, etc.). Finally, they added a class of aerosols that had half of the large fraction compared to the maritime models and referred to these as coastal models (C50, C70, etc.). Summarizing, Gordon and Wang [
5] used aerosol models with a background aerosol and a varying component representing the locally-generated marine contribution.
Full multiple scattering solutions to the (scalar) radiative transfer equation (using the successive order of scattering method) were carried out for each aerosol model (twelve in total) at eight aerosol optical depths (0.05 to 0.8), eight wavelengths (412 to 865 nm) and 33 values of the solar zenith angle (
to
in 2.5 deg increments), for a total of more than 25,000. Carrying out these simulations was a large undertaking, as each solution required about one hour on a 1 MHz workstation. It was effected by utilizing a large number of workstations at the laboratory (RSMAS) during hours when they would otherwise be idle. The result of each computation was twenty Fourier coefficients (in viewing azimuth relative to the Sun,
) of
at 100 viewing angles and the solar zenith angles above, i.e.,
where
is the viewing angle. The functional
was then established as a power series in
via least-squares fitting. This provided look-up-tables from which
could be estimated given
. Examples of the relationship in Equation (
27) are provided in
Figure 6 for aerosols T50 and M99 which have the smallest and largest aerosol particles of the collection of models. The near linearity and similarity of the curves for
and 865 nm suggests that the first term in Equation (
27) is dominant and the multiple scattering effects are approximately the same at both NIR wavelengths. These multiple scattering results were used in the following manner.
For a given pixel we determine the measured values of
at
and
which we denote as
and
, using the assumption that
at both wavelengths. Then, for the
jth model, we can find
i.e., the values of
at
and at
that would be valid if the actual aerosol were identical to that for the
jth model. For the
jth model then, we can find
the value of
if the
jth aerosol model were correct. As mentioned above, since multiple-scattering effects are similar in both NIR bands, e.g.,
Figure 6, the retrieved value of
will be close to the true value, independent of the aerosol model. That is, if we assume several (
N) different models (“candidates”) for the aerosol, but one model say the
ith is correct, each model, e.g., the
jth, will yield an
that will be close in value to the true
. Since the
’s resulting from each of the trial aerosol models are close to one another, rather than pick a single model to estimate this quantity, as it is likely none of the aerosol models are correct, it seemed more reasonable to estimate
through
where
is the value of
derived by assuming that the
jth aerosol model is correct. To obtain a more robust estimate of
, the two models with
values farthest (above and below) the mean are rejected, and a new average is computed omitting these two models. This is executed several times until only 4 models remain, the average of which gives the final value of
. Such an estimate will be most closely bracketed by two of the models, i.e., one with a slightly smaller value and one with a slightly larger value. Call these
and
, respectively. Then
so
Thus, we have two values for . The procedure usually adopted was to assume that falls between and in the same proportion as falls between and . Some kind of assumption such as this is required, but this one is not universally true, in particular when the aerosols are strongly absorbing. The N aerosol models are called “candidate” models, and to the extent that the candidate models were similar to the actual aerosol, this particular assumption worked quite well.
This procedure yielded
, where
t is the diffuse transmittance from the surface to the sensor. Recall that earlier algorithms estimated
t by assuming
was totally diffuse and the atmosphere was aerosol free, i.e.,
. However, admitting that
t depends on the aerosol concentration as well, and assuming that the upwelling radiance just
beneath the water surface was totally diffuse, Gordon and Wang [
5] used the principle of reciprocity (Yang and Gordon [
45]) to compute
for each aerosol model. Then, since
, i.e.,
,
In the Gordon and Wang algorithm, the diffuse transmittance was derived from and in the same manner as and . Finally, from and one can derive and , respectively. The value can then be estimated from these by assuming it falls between and in the same proportion as falls between and .
Tests of the algorithm with simulated data for
, 40, and 60 deg. and
and 45 deg. suggested that the algorithm met the criterion that
for most cases when the true aerosol was SF97’s M80, C80, or T80 (not among the candidate models). Using relationships similar to those in Equations (
12) and (
21), it was found that for pigment concentrations
or 1 mg/m
3, for 95% of the cases tested
. Thus, the algorithm met the goals specified for the SeaWiFS sensor.
4.4. The OCTS and MERIS Multiple Scattering Atmospheric Correction Algorithms
During this time period one ocean color sensor was launched, OCTS in 1996 by NASDA (Japan), and another was in the planning phase, MERIS by ESA. Atmospheric correction algorithms that were similar in concept to the GW94 algorithm were developed for both sensors. The OCTS algorithm developed by Fukushima et al. [
46] (F98) was virtually identical to GW94 in that it employed nine of the SF79 aerosol models (and one Asian dust model, which was strongly absorbing) with radiative transfer (RT) calculations to account for multiple scattering. They did not specify the vertical structure used in the computations, although it is largely irrelevant for the nine SF79 models, which are weakly absorbing.
The principle difference from GW94 was that F98 developed RT look-up-tables for the expansion of
as a function of
, i.e.,
, where
h is a polynomial function of
, as opposed to the tables of
utilized in GW94. Note that since
for a given Sun-viewing geometry, the
g-tables can be converted to
h-tables and vice versa. F98 chose
nm rather than 765 nm to avoid the problem of having to correct the 750 nm band for absorption by atmospheric Oxygen [
47]. This choice restricts use of the algorithm to the open ocean because of the requirement that
(see
Section 5.1).
To operate the algorithm they first inverted the
-
relationship in the red and NIR for each aerosol model:
Then the ’s were averaged (or combined in a weighted average) to find the “derived” estimate of . Finally, the two “best” models among candidates were found based on comparison of the derived with the individual ’s, and the value of was determined at each by interpolation between the two “best” models based on their deviations from . (Recall GW94 used a similar procedure replacing the ’s with ’s.) The performance of this algorithm with simulated data was similar to that in GW94.
In the case of MERIS there were three bands that could be used for atmospheric correction: 705, 775, and 865 mn; however, the basic algorithm for correction utilized the 775 and 865 nm bands. The algorithm developed by Antoine and Morel [
41] was also similar to that of GW94, in that it utilized aerosol models, but differed in more details than the F98 version.
Using radiative transfer theory they confirmed that the “path” reflectance (
) varied almost linearly with the aerosol optical thickness. In their algorithm they chose to develop
, rather than
alone, as a function of the aerosol optical depth, i.e.,
where
is a polynomial in
. Comparing this to GW94, where
was expanded as a function of
, i.e.,
it is easy to see that
Thus, for any aerosol model and Sun-viewing geometry the two representations contain exactly the same information (as they must), i.e., for a given aerosol model and geometry could be derived given and vice versa in much the same manner as could be derived from .
As in GW94, Antoine and Morel [
41] (AM99) developed look-up-tables of
based on radiative transfer computations for an atmosphere bounded by a Fresnel reflecting surface, below which the water was totally absorbing. They used the Monte Carlo method of solving the transfer equation (radiance uncertainty about 2%). The vertical structure of the atmosphere was more elaborate than that of GW94, consisting of fifty 1 km thick layers with the molecular scattering properties assigned to each layer based on the air density. Ozone absorption was assigned to each layer based on the models of Elterman [
48]. The aerosol was assumed to reside in three homogeneous layers (0–2 km, 2–12 km and 12–50 km), but unlike GW94, they assigned (different) aerosols to all three layers.
The layer nearest the surface (the marine boundary layer) contained locally generated aerosols with properties as modeled by SF79. The middle layer (the free troposphere) contained the “background” aerosol that had weak absorption (similar to the SF79 Tropospheric models) and a fixed optical thickness of 0.025 at 550 nm. The upper layer (the stratosphere) contained the non-absorbing stratospheric aerosol composed mostly of sulfuric acid and water and had an optical thickness of 0.005 at 550 nm. These combined with the variable marine boundary layer constituted what AM99 referred to as an “assemblage.” The marine boundary layer could be composed of the SF97’s M70, M80, M90, or M99 aerosols, and the choice of one of these models and the other two (fixed) layers led to one of four “standard assemblages.” Other “assemblages” contain different components, usually strongly absorbing, in the marine boundary layer and/or the free troposphere, e.g., dust or the SF79 Urban models.
The standard assemblages were first used to perform an atmospheric correction as follows. First, the measured is used to determine the set , one for each model within an assemblage, e.g., the jth. Then, the properties of the jth model provides , one value for each model. Next, evaluated at 775 nm is used to generate , which is then compared with the measured counterpart . Finally, in a manner similar to GW94, the two best models, i.e., those that bracket most closely the observed , are used in the correction. The values of are then estimated by mixing the prediction of the two models in the same proportion as the excursions of the two “best” ’s from . If the correction is successful i.e., the resulting water-leaving reflectances pass certain tests (requirements) for Case 1 waters, then the procedure is finished. If the retrieved water-leaving reflectances fail the tests, then the procedure is repeated with other (absorbing) assemblages until success is obtained. Simulations with pseudo-data for cases with weakly-absorbing aerosols suggest that this algorithm has approximately the same accuracy as that of GW94. An interesting difference between AM99 and the earlier two algorithms is the GW94 and F98 had to average the results of the calculations of ’ and to obtain and , respectively, for comparison with the individual models, while in AM99, the values of , to be compared to the model results, came directly from the measured reflectance.
The principal differences between AM99 and GW94 are that (1) AM99 placed aerosols of different types as well as Rayleigh-scattering molecules in all layers with the variable aerosol component in the lowest layer, while GW94 had a single (variable) aerosol type only in the lowest layer, with all of the Rayleigh-scattering molecules in an (aerosol-free) layer above the aerosols; (2) AM99 modeled Ozone absorption using a continuous profile in the generation of the LUTs (discretized into 50 layers, as was atmospheric density for Rayleigh scattering), while GW94 did not include Ozone in their LUTs, but added it during the atmospheric correction process as a single layer at the top of the atmosphere; and (3) AM99 generated their LUTs using Monte Carlo methods with estimated (statistical) errors of the order of 2% while GW94 generated theirs with the successive order of scattering method with much smaller error (order of 0.1%). (With regard to (1) above, it should be noted that many of the GW94 models included both a pure maritime component (locally generated by wave breaking) and a terrestrial component; however, in contrast to AM99, these two components were in the same (lower) layer.) The procedures used in the three atmospheric correction algorithms are summarized in
Table 2.
One should note that in the cases of weakly- or non-absorbing aerosols, e.g., the standard assemblages, the vertical structure of the aerosol is irrelevant at the accuracy required for atmospheric correction. While, for the strongly-absorbing assemblages, e.g., those with dust or the SF97 Urban models, the correction will fail in their presence unless the assumed vertical structure closely approximates the actual structure of the atmosphere (see
Section 4.5).
There were two visible-NIR sensors developed mostly for land and aerosol applications during this era that could also be used for ocean color observations. These instruments were based on the new concept of measuring the angular distribution of radiance exiting the Earth’s atmosphere and as such had the potential to determine the water-leaving radiance in several directions, albeit at slightly different times (up to 7 min). MISR [
49] provided measurements of the radiance exiting the atmosphere at nine viewing directions distributed fore and aft of the satellite’s orbit track, plus the nadir, and at four wavelengths (nominally 440, 550, 670, and 865 nm). POLDER [
50] was capable of measuring the exiting radiance over a continuous range of viewing directions up to about 50 deg. from nadir, at wavelengths similar to MISR, but with the additional feature that the polarization of the radiance could be determined at 443, 670, and 865 nm as well. The atmospheric correction algorithms described in earlier sections could also be used with these sensors by applying them directly with
and
nm. However, the availability of multi-directional reflectance provides further possibilities. In the case of MISR, Wang and Gordon [
51,
52] showed that the angular variation in
measured by MISR, when compared to that produced by modeled aerosols, e.g., SF97, could be used to obtain the aerosol properties, and Gordon [
44] built on this to show that atmospheric correction could be effected using a single band (865 nm, for which
) by using the
angular distribution of
to select SF79 aerosol models. Atmospheric correction based on a single band in the NIR is attractive because of its potential for application to Case 2 waters. Simulations showed that the efficacy of such an algorithm was slightly better than the SeaWiFS-type algorithm for Northern Hemisphere Summer geometries (small solar zenith angle), but much worse for Northern Hemisphere Winter geometries (large solar zenith angle). Note that such an algorithm requires that the phase function of at least one of the candidate aerosol models must closely match that of the actual aerosol extant in the atmosphere, i.e., a much stronger constraint than that of the GW94-type algorithm.
Frouin et al. [
53] developed an approach with POLDER that was similar to that of GW94, but with some exceptions. They proposed to utilize the band at 670 nm along with
nm and
nm. In essence, for each pixel they computed
where
i is an index over the set of viewing directions extant for the given pixel, and
w is a weighting factor (set equal to 1 for the first instrument). This was then compared to its counterpart for RT simulations using the SF79 models to select the two “best” models, which are mixed by interpolation similar to GW94, with similar results [
53]. To the author’s knowledge, the polarization capabilities of POLDER have not been used to assist in the estimation of
.
Comparison of the performance of the SeaWiFS/MODIS, OCTS, MERIS, and POLDER algorithms operating on a synthetic data set can be found in IOCCG Report No. 10 [
54].
4.5. Absorbing Aerosols
Aerosols that have strong absorption, either neutral with wavelength, e.g., carboniferous aerosols from coal-fired power plants, or selective with wavelength such as desert dust with a significant Iron content (strong absorption in the blue), pose significant challenges for atmospheric correction. The first, and probably the most important, challenge is that all algorithms discussed so far utilized the observation that the radiance exiting a weakly-absorbing atmosphere can be considered to be independent of the aerosol vertical structure, at least at the accuracy required for atmospheric correction. This is no longer true when the aerosols are absorbing. An example of this is provided in
Figure 7 for the Shettle and Fenn [
36] Urban model at relative humidity 80% (U80). For this model
–0.74 for
–865 nm, respectively, and also,
. The behavior seen in the figure, a progressive decrease in
as a fixed amount of aerosol is mixed higher and higher into the atmosphere, is a characteristic behavior of an absorbing aerosol and demonstrates the strong dependence on vertical structure. One should also note from
Figure 7 that measurement of
in the NIR provides little or no information regarding the vertical structure even when the physical properties of the aerosol are known. However, if the aerosol is suspected to be absorbing, and candidate models are restricted to those with similar optical properties (e.g., U50, U70, U90, and U99, when the actual aerosol is U80),
Figure 8 shows that an acceptable atmospheric correction can be obtained with the Gordon and Wang [
5] algorithm,
but only when the actual vertical structure is close to that of the candidate models, for which all of the aerosol is in a thin layer at the surface. Thus, it appears that atmospheric correction in the presence of an absorbing aerosol is possible if (1) the general nature of the aerosol is known, i.e., in the
Figure 8 example the aerosol type is known to be Urban, and (2) the vertical distribution of the aerosol is known to be similar to that in the candidate models.
How does one know that absorbing aerosols are present? The only way seems to be the appearance of unexpected results following the usual atmospheric correction, which assumes the aerosol is non-absorbing. For example, consider the U480 case in
Figure 7. If one assumed that the aerosol was non-absorbing, then
would be approximately the extrapolation from that at 765 and 865 nm in
Figure 7, or approximately 0.026; however, the true value would be about 0.014, so the estimated value is 0.012 too high, or
would be 0.012 too low.
Figure 5 shows that if the true value of
were 0.1 mg/m
3, then an error of this magnitude would lead to a retrieved
mg/m
3, so one manifestation of the presence of the absorbing aerosol might be a much higher-than-expected
. Another manifestation would be a retrieved
, that could occur in waters with even moderate
.
Another method of detecting atmospheric correction problems such as the unknown presence of absorbing aerosols is based on the right panel of
Figure 5, which shows that
varies within a restricted range of values. The Antoine and Morel [
41] algorithm actually used
, which shows even less variation than
, to detect failure of the atmospheric correction using their standard assemblages, i.e., retrieved values of
falling outside the range of acceptable values. When failure was detected, atmospheric correction was carried out using another set of assemblages tailored to the nature of the particular situation, e.g., dust, urban pollution, etc. Thus, their procedure for the detection of absorbing aerosols was based on the fact that one has an
a priori expectation of the range of water-leaving reflectance values in some cases. That is, one has a model, either mathematical or experimental, of what to
expect for
for a given location, e.g., for Case 1 waters.
Carrying this expectation further, Gordon, Du, and Zhang [
55] proposed a method of solution, called the “spectral matching algorithm” (SMA), that
required for all values of
to fit a particular semi-analytic bio-optical model [
21]. In their algorithm, first the bio-optical model provided
as a function of
and a particle scattering parameter
, i.e., it was used to find
. Using the same approximation as in Equation (
30), these provided trial values of
for the
ith aerosol model and waters characterized by particular values of
and
. Next, using the usual assumption that
,
for the
ith aerosol model was determined from
, and this yielded
for the other spectral bands, i.e., the value
would have if the
ith aerosol model were correct. Finally, the residual
where
n is the number of visible wavelengths, was examined. The value of
was computed for each aerosol model and discrete set of ocean parameters. Note that, if for a particular parameter combination
and
and a particular aerosol model,
, then an exact solution of the problem has been found, i.e., the top-of-atmosphere
is consistent with the water parameters and the atmospheric model. This is most unlikely, so normally the set of parameters
i,
, and
that yielded the smallest
would be chosen as the best, i.e., the solution to the problem. However, knowing that it is also unlikely that the “correct” model is one of the set of candidates, averaging over some number of the best retrievals (i.e., retrievals with the lowest values of
) to obtain the retrieved ocean and aerosol parameters seemed more reasonable. Gordon, Du, and Zhang [
55] averaged over the best ten.
To test this algorithm they generated pseudo data using the M80, C80, T80, and U80 aerosol models for the atmosphere and the Gordon et al. [
21] water-leaving reflectance model. Note that the same water model is used to generate the pseudo data and to operate the SMA; so this represented the optimum situation, i.e., a correct bio-optical model. Sixteen aerosol models are used as candidates: M, C, T and U models with RH = 50, 70, 90, and 99%. Note the aerosol models used to generate the pseudo data were not part of the candidate set. For this test, the vertical structure of the aerosol in the generated pseudo data is the same as that in the candidate models—again, an optimum situation. The test demonstrated that with the above limitations—correct water model and correct aerosol vertical distribution—the SMA performed very well. Of particular importance was the observation that the
algorithm had no difficulty determining the presence of strongly absorbing aerosols and performed as well as the Gordon and Wang [
5] algorithm when the aerosol was non-absorbing. In a real-world operation of the SMA, the major stumbling block was that one would have to incorporate information regarding the aerosol vertical structure when it was determined that the aerosol was absorbing.
Chomko and Gordon [
56] proposed a different approach, but similar in spirit, to effecting atmospheric correction both in the presence of absorbing and non-absorbing aerosols. Rather than using the realistic aerosol models of SF79, they assumed a simple power-law size distribution:
with fixed values of
,
, and
. The aerosol was assumed to be homogeneous in nature (all particles have the same composition), and was completely specified by
K,
,
, and
, where the refractive index of the (assumed-spherical) particles is
. Thus,
. (Note, for a given set
,
, and
,
.) They used the same water model as the SMA, so
. For SeaWiFS this provided eight equations of the form
the set of which were “solved” by standard optimization techniques, hence the name spectral optimization algorithm (SOA). Applying this to simulated “data” developed using the SF79 aerosol models T80, C80 and M80, they found good retrievals of
(range,
mg/m
3), showing that accurate aerosol size distributions are not required to effect atmospheric correction; however, on the contrary, as one would expect, the derived aerosol optical thickness showed large errors. Excellent retrievals were also obtained for strongly absorbing aerosols (U80) but, again, only when the correct vertical distribution of the aerosols was known and employed in the algorithm.
The SOA was tested using SeaWiFS imagery from the Middle Atlantic Bight and the Sargasso Sea by Chomko and Gordon [
57]. They examined imagery from two days (DOY 279 and 281, 1997), one of which had a turbid atmosphere and the other a clear atmosphere. They found good continuity between the derived pigment concentrations over the two days, and a complete decoupling between the atmosphere and the ocean patterns for a relatively turbid atmosphere, but for the clear day water patterns were clearly evident in the retrieved aerosol products. This atmosphere-ocean coupling was significantly reduced by using a more sophisticated water model (Chomko et al. [
58]). In that study, the water was modeled using radiative transfer along with the absorption coefficients of phytoplankton (
) and colored detrital material (CDM,
), and the backscattering by phytoplankton (
). These parameters were further related to wavelength through
where
C is the concentration of chlorophyll
a,
is the chlorophyll-specific absorption coefficient of phytoplankton,
S characterized the spectral variation of CDM,
specified the spectral variation of phytoplankton backscatter, and
nm. The parameters
,
S, and
were optimized for Case 1 waters by Maritorena et al. [
59] using in-situ data. Thus, the unknown parameters of the water model were
C,
, and
, i.e.,
. When this model was tested with SeaWiFS imagery, a more complete decoupling between the water and the atmosphere was obtained, and good agreement was found between
and its proxy measured contemporaneously by airborne lidar, and between the derived
C and that from standard SeaWiFS processing. To apply such an algorithm to an absorbing aerosol that was distributed higher in the atmosphere, rather than in a layer at the surface, would require look-up-tables tailored to each likely aerosol vertical distribution. The required resolution in the vertical structure for effecting atmospheric correction is yet to be determined.
It should be pointed out that the models in this section that relate
to the absorption and backscattering coefficients of the water’s constituents are approximate in that the resulting upward radiance distribution is taken to be uniform, an assumption that must be invoked because the scattering phase function of the constituents in the water is not used in the calculation. A more realistic model has been presented by Morel, Antoine, and Gentili [
60], in which the angular distribution of
is directly related to
C. This model required a full solution to the radiative transfer equation (including the Raman scattering) and the results were given in the form of lookup tables. They are for Case 1 waters only, and in application to atmospheric correction, have been used in the SeaWiFS implementation to remedy failures of the “black pixel” approximation (see
Section 5.1).
4.6. Second-Order Effects
There are many effects that were omitted in the development of atmospheric correction algorithms as their influence was thought to be much smaller than the effect of aerosols, e.g., they were second-order effects. Among these were whitecaps from breaking waves, Sun glint cause by surface waves, departure from the plane-parallel atmosphere assumption due to the curvature of the Earth, and polarization of light in the ocean-atmosphere system. There were of course issues that were due to idiosyncrasies of specific sensors, e.g., in the case of SeaWiFS interference due to the Oxygen “A” band absorption near 765 nm, non-ideal spectral response functions, or excessive instrument response to the polarization of
as in SeaWiFS and MODIS. Although important, these sensor-specific issues will not be discussed here (e.g., see [
47,
61,
62]).
4.6.1. Reflectance Augmentation from Whitecaps
Whitecaps from breaking waves introduce additional radiance reflected from the surface over and above the water-leaving radiance. Their effect is small but can exceed
in the red and NIR, and therefore interfere with atmospheric correction [
63]. The fraction of the water surface covered by whitecaps is strongly dependent on wind speed (
, where
to 3.5, and
W is the surface, or ten-meter-height, windspeed). Thus, their reflectance is strongly dependent on
W as well. (A complete discussion of early studies aimed at relating whitecap coverage to wind speed and other parameters, such as atmospheric stability, etc., can be found in Stramska and Petelski [
64]). Most studies attempted to relate the whitecap coverage (in % of surface covered) to the wind speed, and then to the augmentation of reflectance by multiplying the coverage by their effective reflectance (e.g., Koepke [
65]). Whitecap coverage was usually determined from ship-board photographs of the water surface [
65,
66], but some measurements were obtained radiometrically [
67,
68,
69]. Both determinations of coverage are difficult and the “derived” increased reflectance has significant variance (standard deviation of coverage ∼ coverage). Frouin et al. [
67] showed that sea foam displays a clear decrease in reflectance from the visible to the NIR. Moore et al. [
69] gave the radiometrically-determined increase in reflectance due to whitecaps as
, where
S is the spectral variation being unity for all bands in the visible and 0.95 and 0.85, respectively for 750 and 860 nm. In contrast, from the Stramska and Petelski [
64] photographically-determined windspeed-coverage relationship and the average reflectance of whitecaps (22% [
65]), the reflectance was given by
, where
, 0.765, and 0.645 at 670, 765, and 865 nm, respectively, (and unity elsewhere) was taken from Frouin et al. [
67], and
is the wind speed 10 m above the surface. The Stramska and Petelski [
64] formulation was used to estimate and remove the whitecap contribution from
in the SeaWiFS/MODIS processing.
4.6.2. Sun Glint
Sun glint, the specular reflection of direct solar radiation from the water surface to the sensor, can either be a zeroth, first, or second order effect. Looking near the specular image of the Sun, i.e., the position of the Sun in an image for a flat water surface, results in a huge radiance that exceeds all other terms in the radiance budget and can even saturate the sensor. This region, the radiance and spatial size of which depends on the wind speed (Cox and Munk [
70]) was always simply masked from processing (not processed). Further away from the specular image the radiance is still large and, although it could be estimated given the wind speed and direction using the Cox and Munk [
70] surface slope probability density, the relationship between the probability density and wind speed was not well enough known to be able to predict the radiance accurately, so this region was masked as well. In fact, the masks were so conservative that Sun glint was generally ignored. There were many Sun glint studies in the early 2000’s aimed mostly in assessing the glint and reducing the size of the glint-masked regions, e.g., Wang and Bailey [
71] implemented a correction using the Cox-Munk slope distribution that operated when the glitter reflectance was predicted to be less than about 0.005–0.01.
4.6.3. Polarization of Radiance in Ocean-Atmosphere System
In the Gordon and Wang [
5] algorithm, the Rayleigh scattering contribution
was computed using full
vector radiative transfer theory, while the aerosol look-up-tables for the relationships in Equations (
28) and (
29) were generated using
scalar radiative transfer theory. Understanding that ignoring polarization in radiative transfer can result in errors of a few percent, Gordon [
44], in a restricted study, showed that the algorithm could perform somewhat better (applied to pseudo data) when vector theory was employed, but did not implement it for SeaWiFS. In a later study, Wang [
72], recomputed the entire set of tables using vector theory. Use of the new tables showed significant improvement when the algorithm was operated with synthetic data; however, after vicarious recalibration, a procedure that is required whenever an algorithm is changed, there was a negligible difference in the derived products, i.e., the
’s. In part, this result was a demonstration of the power of vicarious calibration and in part a demonstration of the robustness of the algorithm concept. When a new set of candidate aerosol models were introduced [
73] (
Section 5.4 below) the computation of the aerosol look-up-tables was carried out using a vector radiative transfer procedure.
4.6.4. Curvature of Earth’s Surface
The GW94 atmospheric correction algorithm (
Section 4.3) was based on radiative transfer (single and multiple scattering) within an assumed plane-parallel geometry (PPG) for the atmosphere as opposed to the correct spherical-shell geometry (SSG). Although one would expect the PPG to be an excellent approximation, its applicability will clearly break down when the solar zenith angle becomes large, i.e., at high latitudes. Gordon and Ding [
74] tested the applicability of the PPG by creating pseudo data in SSG and then using the GW94 algorithm (based on PPG) to effect atmospheric correction. They found that for
deg the algorithm performed as well as with pseudo data created in PPG; however, for
deg, significant deviations were observed. In addition, they discovered that the additional error could be almost eliminated completely by computing
in SSG rather than in PPG. They also presented a simplified computation scheme for
in SSG based on the observation [
75] that the ratio of the true
to that computed assuming single scattering was
approximately the same in PPG and SSG.
6. Concluding Thoughts
There is little doubt that the success of the CZCS led to the development of ocean color as a significant component of research related to living marine resources and to climate change on a global scale. At present count there have been 29 ocean color sensors following CZCS placed in space and 8 in the planning stage (IOCCG website). Its success led directly to the follow-on sensors SeaWiFS, OCTS, and MERIS. It seems reasonable to assert that the success of the CZCS followed directly from the ability to execute atmospheric correction to an accuracy that allowed application of the bio-optical algorithms that retrieved the pigment concentration in the global oceans to an accuracy well within a factor of two. This ability resulted partly from a good measure of serendipity.
The choice of the CZCS spectral bands, although well chosen for retrieving pigments given
, was inadequate for atmospheric correction, a fact already known prior to launch [
39]. Apparently there was no consideration given to atmospheric correction in the design of the CZCS, only to ways of possibly minimizing the atmospheric contamination, e.g., using band ratios
at various wavelengths [
109]. The red band, that ultimately was central to atmospheric correction in the open ocean, was actually included on CZCS for the detection of chlorophyll at high concentrations and had significant water-leaving radiance only at the very high chlorophyll concentrations that one might find in the coastal zone, estuaries, and lakes. (Recall that chlorophyll
a has an absorption band at 670 nm, so any enhanced radiance there had to be due to the backscattering of plankton containing the pigment, their associated detritus, and whatever other kinds of particulates are present–a very indirect and insensitive method of measuring chlorophyll.) The NIR band on CZCS was intended for surface vegetation and was too insensitive for oceanic application. Interestingly, the excess radiance in the red band, although minimal, prevented atmospheric correction in high pigment areas, e.g., the coastal zones, so we had the conundrum that the
Coastal Zone Color Scanner performed poorly near the coasts, but, in the end, very well in the open ocean. It is interesting to note that the CZCS deficiencies led most of the oceanographic community (except the experiment team,
Table 1) in the prelaunch era to consider the CZCS project a waste of precious funding. Had the red band not been included, atmospheric correction on the open ocean may not have been possible, and the CZCS would likely have been considered a failure. As the time between the launches of CZCS and SeaWiFS was 19 years, one is led to wonder how much that time interval would have expanded if CZCS atmospheric correction had failed. But it did not fail, the CZCS mission proved that chlorophyll (actually pigments) could be measured on a global scale and, concomitantly, an indication of primary productivity on a global scale, and the CZCS truly became the “proof-of-concept” mission.
What are the optical characteristics of the atmosphere that led to the possibility of atmospheric correction of CZCS? First and foremost is the fact that the atmosphere is basically optically thin, meaning that single-scattering of sunlight provides a reasonable approximation to the atmosphere’s contribution to
(
Figure 1 and
Figure 2). [It also means that one can image through the atmosphere, an obvious requirement for space-borne ocean color sensors.] In particular, molecular scattering, the largest contribution to
in the blue (the most important spectral band for estimating
in the open ocean,
Figure 5, left panel), could be computed using single-scattering theory, as it was in the initial validation studies [
14,
16]. It enabled atmospheric correction on a pixel-by-pixel basis to be carried out
with the computational resources available at that time. Second, because radiative transfer in the atmosphere was close to single scattering, Rayleigh and aerosol scattering were independent (
) so
([
10], and Equation (
4)). This meant that
at different wavelengths were proportional to one another (Equation (
5)), and therefore could be determined from one another. Third, although the proportionality “constant”
(Equation (
10)) depends on the aerosol phase function, the fact that the dependence is weak, especially for marine aerosols, which are the main aerosol in the open ocean, meant that
, which can be accurately modeled by Ängstrom’s law. Thus, if the aerosol’s Ängstrom exponent remained constant over a scene (a constant aerosol “type,” i.e., uniform aerosol size
frequency distribution and refractive index) and could be determined in one part of the scene, then
could be determined from
. [Although not a characteristic of the atmosphere, a characteristic of Case 1 waters that facilitated this process was the “clear water” radiance concept [
18].] The presence of the red band (670 nm), where
for most of the world’s oceans, facilitated the determination of
, from which
could be found leading to
for use in the bio-optical algorithms (Equation (
12)). These characteristics of the atmosphere (and “clear” water), enabled the arguably “rapid,” correction of imagery obtained early in the mission and preparation and publication of initial results in a prestigious journal [
16,
17], fending off critics. In addition, the fact that
for a marine aerosol (e.g., M99) enabled processing of the full CZCS archive without the cumbersome process of locating clear water to determine
in several sub-scenes of each of the approximately 68,000 images [
35,
38].
What would happen if one of the links in the above chain were broken, e.g., the weak dependence of the aerosol phase function on wavelength. The manifestation of failure of this particular link would be that
determined from one portion of the image, e.g., “clear water,” could not be “carried over” to another portion of the image (without explicit knowledge of the phase functions, Equation (
10)), and the resulting “corrected” image with the assumed constant
would contain easily-identifiable atmospheric features (from spatial scale differences between oceanic and atmospheric features)—a failure of atmospheric correction in particular and the CZCS mission in general.
I cannot over stress the importance of the fact that atmospheric radiative transfer, at the level of accuracy that was required for CZCS, is well-approximated by single scattering. In early 1979, when CZCS imagery became available, the image processing system (the Atmospheric and Oceanic Image Processing System—AOIPS at NASA/GSFC) available for our use with CZCS imagery was powered by a DEC PDP 11/55 mini computer (circa 1976). Although state-of-the-art, it could only perform simple manipulation of the various spectral bands (image planes) of an image. Initially, I computed the Rayleigh reflectance
in scan coordinates (line#–pixel#) using a Univac 1106 mainframe (circa 1965) computer, which would be a “toy” by today’s standards. In this form
could be simply subtracted from
to form
using AOIPS. Then using
determined by assuming
and
, the desired
was found (Equation (
9), with
). The first CZCS image to be atmospherically corrected in this manner (or any manner) is provided in
Figure 12. The correction was based on measured
’s at a location in the upper portion of the image, enabling
to be determined there, and then used throughout the sub-scene. It was photographed from the AOIPS monitor screen by the author as atmospheric correction of the 443 nm band of CZCS took place, line-by-line. The time to compute this simple correction for a single scan line (512 pixels in the sub-scene) was of the order of ten seconds. As I watched the correction take place, it was as if a veil were slowly being pulled downward off the image—and it was exhilarating. Were the algorithm to require multiple scattering computations, it is hard to envisage (1) how imagery could be processed in a fashion timely enough to forestall designation of the CZCS mission as a failure, and (2) how more than a few CZCS scenes could have been processed in the 1980s. “Fortunately, by the late 1980’s, with the introduction of the DEC VAX line of super mini-computers, and software developed at the University of Miami and installed at NASA/GSFC, it was possible to process the entire CZCS archive of images (∼68,000 images, each with 1968 pixels/line and 970 lines) and demonstrate the potential of ocean color imagery [
35,
38].
Comparison of the radiometric performance of SeaWiFS with CZCS showed that a correction algorithm for SeaWiFS had to be significantly more accurate than the simple CZCS algorithm, and in fact required a full multiple-scattering algorithm. Although the spectral bands in the NIR facilitated atmospheric correction and indeed a single-scattering version using these bands was quite accurate [
42], it was not accurate enough. M. Wang and I begin working on a multiple-scattering version (GW94) in 1989 when SeaWiFS was approved by NASA. Although multiple scattering was required, success with the single-scattering algorithm suggested that we should try and use the same structure. Early tests were carried out with aerosol models for which both the phase function and the extinction coefficients were independent of wavelength to approximate a marine aerosol (the aerosol was also non-absorbing). Results were encouraging, suggesting that a correction that was roughly an order of magnitude better than CZCS could be effected [
110,
111]. However, when spectral variations of the aerosol optical properties expected for realistic aerosol models, e.g., SF79, were incorporated into the algorithms and its pseudo data, the resulting accuracy was reduced by a factor of two to three—fortunately, still accurate enough for SeaWiFS.
Implementation of the algorithm for operational SeaWiFS processing posed problems. First, the generation of the LUTs containing the expansion coefficients
a,
b and
c in Equation (
27) was a significant computational problem. The computations required solving the radiative transfer equation at eight wavelengths, 12 aerosol models, 33 values of the solar zenith angles, and eight values of the aerosol optical depth for a total of 25,344 simulations. Our radiative transfer code (successive order of scattering) required one hour per simulation on a DECstation 3100, so the total time required for a single work station was approximately 2.9 years. The actual computations were carried out in about one month using a large number of work stations at RSMAS during times when they would otherwise be idle (thanks to Robert Evans and James Brown for developing a scheduler that effected this task). Second, the LUTs had to be available for rapid access for pixel-by-pixel processing, i.e., stored in the computer’s main memory. Because our successive-order code used Fourier decomposition in the azimuth, rather than store tables of
a,
b and
c for a set of azimuth angles, we stored the first 15 of their Fourier coefficients. The final LUTs contained over 100 Mb of data: 3 × 15 Fourier coefficients for each case of 8
’s, 100 viewing angles, 33 solar zenith angles, and 12 aerosol models. Systems with such large memory were not readily available at that time. But, by launch (1997) the systems with sufficient memory were available and the GW94 algorithm could be operated on SeaWiFS data in a timely manner. Lastly, the prototype algorithm was only capable of processing three to four pixels per second on a DECstation 3100. This was too slow to even consider processing imagery in real time, a goal for the SeaWiFS mission; however, by launch, processors were fast enough to process the imagery efficiently, and by the end of the mission, the entire SeaWiFS data set could be reprocessed in a matter of weeks. Saved by Moore’s law again.
A final development that became an indispensable development tool for SeaWiFS and for future sensors was the SeaDAS computer program developed and maintained at NASA/GSFC which was freely available to all, and allowed processing of ocean color imagery without the significant complications of navigation, calibration, etc., using a personal computer in the comfort of one’s office. As I can personally attest, the availability of the source code for this program made it possible to try new approaches to atmospheric correction with a minimum of effort [
56,
57,
58,
79,
82,
93,
94].
Here, I have tried to describe the ocean color atmospheric correction algorithm and its development within the context of the time, i.e., the available computational resources. Through a good measure of luck and some insight, the CZCS and follow-on sensors worked sufficiently well to enable ocean color remote sensing to achieve respectability, and become a significant player in global marine-biological and climate-related studies. Additional reflections on the journey from CZCS to SeaWiFS can be found in Gordon [
112] and Acker [
113].