Inversion of Satellite Ocean-Color Data
Robert Frouin
Scripps Institution of Oceanography
La Jolla, California, USA

ADEOS-2 AMSR/GLI Workshop,Tsukuba, Japan, 30 January 2007

Collaborators
Pierre-Yves Deschamps, University of Lille
Hajime Fukushima, Tokai University
Lydwine Gross-Colzy, Capgemini, Toulouse
Hiroshi Murakami, EORC, JAXA
Takashi Nakajima, Tokai University
Bruno Pelletier, University of MontpellierLinear combination of TOA reflectancePrinciple of the Method-The top-of-atmosphere reflectance in selected spectral bands iscombined linearly, after correction for molecular scattering and sun
glint contributions.
-The coefficients of the linear combination minimize the perturbing
effects, due to scattering and absorption by aerosols, and reflection
by the surface. These effects are decomposed into a polynomial or
principal components.
-The spectral bands are selected so that the linear combination issensitive to chlorophyll-a concentration.

Minimization of Perturbing Effects-TOA reflectance
RTOA(λ) = Rm(λ) + Ra(λ) + Rma(λ) + Tm(λ)Ta(λ)Rw(λ)Rc(λ) = RTOA(λ) – Rm(λ)= Ra(λ) + Rma(λ) + Tm(λ)Ta(λ)Rw(λ)= R’(λ) + Tm(λ)Ta(λ)Rw(λ)-Linearly combining the corrected reflectance Rc in spectral
bands centered at λi yields the index
I = ∑i [aiRc(λi)] = ∑i [ai R’(λi)] + ∑i [ai Tm(λi)Ta(λi)Rw(λi)]-To eliminate most of the atmospheric influence on I, one has
to find coefficients ai that fulfill∑i [ai R’(λi)] = 0

-For this, R’(λi) is decomposed in a polynomial or principal
components, i.e.,R’(λi) ≈ ∑j [bjeji]-In general, a satisfactory representation can be obtained with
only a few eigenvectors, ej, since R’ is a smooth function ofwavelength.
-Substituting R’ by its linear expression, we obtain∑i {[ai ∑j[bjeji]} = ∑j {[bj ∑i aieji]} = 0-To satisfy this equation, it is sufficient to have, for each ej∑i [aieji] = 0-This system of linear equations is solved using p wavelengths,
n = p-1 eigenvectors, and a1 = 1.
-Note that the coefficients bj, which vary with geometry andgeophysical conditions do not need to be known.

Figure 1. Application of the linear combination method to simulated
GLI imagery. Perturbing effects are expressed as a polynomial.

Figure 2. Application of the linear combination method to global GLI
data. (Courtesy of H. Murakami, JAXA.)

Figure 3. GLI Version 2 chlorophyll-a concentration product.
(Courtesy of H.Murakami, JAXA.)

11/19/2002

Figure 4. Application of the linear combination method to SeaWiFS imagery.
The perturbing effects are decomposed into principal components.

Retrieval of the Ocean Signal-First, a set of spectral bands is selected in the red and near infrared,
for which the water body can be considered black, except in one of thespectral bands.
-Second, other sets of spectral bands are selected, that progressively
include shorter wavelengths. At each step, only marine reflectance inone spectral band is unknown and therefore estimated.

[λi , i = 1, 2, …,8] = [412, 443, 490, 510, 555, 670, 765, 865]

1: Tm(λ6)Ta(λ6)Rw(λ6) ≈ ∑i [aiRc(λi)] , i = 6, 7, 8
2: Tm(λ5)Ta(λ5)Rw(λ5) ≈ - a6Tm(λ6)Ta(λ6)Rw(λ6)+ ∑i [aiRc(λi)] , i = 5, 6, 7.
.
.
6: Tm(λ1)Ta(λ1)Rw(λ) ≈ - a2Tm(λ2 )Ta(λ2)Rw(λ2) - a3Tm(λ3 )Ta(λ3)Rw(λ3)+ ∑i [aiRc(λi)] , i = 1, 2, 3

670 nm

555 nm

510 nm

490 nm443 nm412 nmFigure 5. Residual perturbing effects of the atmosphere and surface,i.e., ∑i [ai R’(λi)], as a function of aerosol optical thickness at 555 nm.

670 nm

555 nm

510 nm

490 nm443 nm412 nmFigure 6. Comparison of estimated and actual ocean signal, i.e., Tm(λ)Ta(λ)Rw(λ),
when the optical thickness at 550 nm is between 0.1 and 0.2.

670 nm

555 nm

510 nm

490 nm443 nm412 nmFigure 7. Performance statistics for the retrieval of Tm(λ)Ta(λ)Rw(λ).
Bias (solid circles) and standard deviation (open circles) are given as
a function of aerosol optical thickness at 550 nm.

Conclusions-The perturbing influence of the atmosphere and surface is minimized
adequately for each set of wavelengths, except when aerosol loading islarge. The residual effects exhibit a bias increasing with aerosol
optical thickness. The bias can be reduced by taking into account thelast eigenvector of the decomposition into principal components, butjust globally.
-Errors in the estimated ocean signal, i.e., Tm(λ)Ta(λ)Rw(λ), increase withdecreasing wavelength (residual effects at longer wavelengthspropagate) and with increasing aerosol optical thickness. They become
unacceptable when the optical thickness at 550 nm is above 0.3.
-Performance can be improved by optimizing the sets of selected
wavelengths, or by using a knowledge of the aerosol optical thickness,
which can be estimated from the satellite data.

Fields of Nonlinear Regression Models

ProblemTo estimate marine reflectance ρw from top-of-atmospherereflectance ρTOA and angular variables t without knowing theother variables x that influence the radiative transfer in theocean-atmosphere system

Methodology-Explanatory variables (ρTOA) are considered separately from theconditioning variables (t).
-An inverse model is attached to each t, and the attachment iscontinuous, i.e., the solution is represented by a continuum ofparameterized statistical models (a field of non-linear regressionmodels) indexed by t:

ρw = ζt(ρTOA) + ε
where ε is the residual of the modeling.

Methodology (cont.)
Ridge functions, selected for their approximation properties,
especially density, are used to define the statistical models
explaining ρw from ρTOA and t:ζtj(ρTOA) = Σi = 1, …, n cijh(ai. ρTOA + bi)
ρwj = ζtj(ρTOA) + εjwhere ai(t), bi(t), and cij(t) are the model parameters.

Simulated Data Sets
62,000 joint samples of ρTOA and ρw split in two data sets, D0e andD0v, for construction and validation. Noisy versions D1e, D1v, D2e, andD2v generated, by adding 1 and 2% of noise to ρTOA. The noise isdefined by:ρTOAj’ = ρTOAj + νcρTOAj + νucjρTOAj
where νc and νucj are random variables uniformly distributed on theinterval [-ν/200, ν/200], where ν is the total amount of noise inpercent.

Function Field Construction-The free parameters of the field, i.e., the maps ai(t), bi(t), andcij(t), are estimated by multi-linear interpolation on a regular gridcovering the range of t.
-The adjustment is considered in the least-square sense, andminimization of the mean squared error is carried out using astochastic gradient descent algorithm.

Function Field Construction (cont.)-A sufficient number of n = 15 basis functions was selected via
simulations, and three fields of this kind, ζ0, ζ1, and ζ2 wereconstructed for 0, 1, and 2% of noise.
-Since the components ζtj take their values in the same vector space(the vector space spanned by the linear combinations of ridgefunctions), the approach is not equivalent to separate retrievals ona component-by-component basis.

Theoretical Results for GLI

Table 1. Root Mean Squared error (RMS) and Root Mean Squared Relative error(RMSR) for the models ζ0 and ζ1 evaluated on the construction and validation data
sets (D0e and D0v) and on 1% noisy versions of them (D1e and D1v).

Figure 8. Estimated versus expected marine reflectancefor model ζ1 adjusted on 1% noisy data.

Figure 9. Conditional quantiles (of order 0.1, 0.25, 0.5, 0.75, and 0.9) of theresidual ρw error distributions as a function of aerosol optical thickness at
550nm for model ζ1 applied to 1% noisy data.

Figure 10. Conditional quantiles (of order 0.1, 0.25, 0.5, 0.75, and 0.9) of theresidual ρw error distributions as a function of scatering angle for model ζ1applied to 1% noisy data.

Figure 11. ρw(443)/ρw(545) as a function of [Chl-a] for theoretical ρw
and for ρw estimated by ζ1 from 1% noisy data.

Application to SeaWiFS Imagery-Function field methodology tested on SeaWiFS imagery acquired on
day 323 of year 2002 over Southern California.

-ζt2 gives large differences in ρw compared with SeaDAS values,

resulting in 78% difference in chlorophyll-a concentration on average.
-Differences may be explained by large noise level on ρTOA (e.g., 14%at 412 nm), due to RT modeling uncertainties.
-Noise distribution estimated on 2,000 randomly selected pixels of
the imagery, and introduced during the execution of the stochastic
fitting algorithm, yielding function field ζt*.

Figure 12. Marine reflectance ρw estimated by ζ* for SeaWiFS
imagery acquired on day 323 of year 2002 over Southern California.

Figure 13. Marine reflectance ρw estimated by SeaDAS for SeaWiFSimagery acquired on day 323 of year 2002 over Southern California.

Figure 14. Histograms of marine reflectance ρw retrieved by SeaDAS and ζ*.

Figure 15. Marine reflectance spectra retrieved by SeaDAS and ζ*.

Figure 16. [Chl-a] retrieved by SeaDAS and ζ*, fractional difference, andhistograms for SeaWiFS imagery acquired on day 323 of year 2002 overSouthern California. Average difference is 19.6%.

Conclusions
Fields of non-linear regression models emerge as solutions to acontinuum of similar statistical inverse problems. They match well
the characteristics of the remote sensing problem, allowing
separation of the explanatory variables (ρTOA) from the conditioning
variables (t).
The inversion is robust, with good generalization, andcomputationally efficient. The retrievals of ρw are accurate, with anerror uniform over the entire range of ρw values. Situations ofabsorbing aerosols are handled well.
For noise levels up to a few percent, a general noise scheme may beappropriate, but for large noise levels, the noise distribution needs
to be estimated. A plug-in approach may be reasonable.