Current Advances in Uncertainty Estimation of Earth Observation Products of Water Quality

Remote sensing data over a water body are related to the physical and biological properties of water constituents through inherent optical properties (IOPs). These IOPs characterize the absorption and scattering of the water column and are used as proxies to water quality variables. The scientific procedure to derive IOPs from ship/space borne remote sensing data can be divided into three steps: iforward modeling, relates the radiometric data to the IOPs of the water column; iiparametrization, defines the minimal set of IOPs whose values completely characterize the observed radiance; iiiinversion, derives the values of IOPs, and hence water quality variables, from radiometric data.


Introduction
Remote sensing data over a water body are related to the physical and biological properties of water constituents through inherent optical properties (IOPs).These IOPs characterize the absorption and scattering of the water column and are used as proxies to water quality variables.The scientific procedure to derive IOPs from ship/space borne remote sensing data can be divided into three steps: i-forward modeling, relates the radiometric data to the IOPs of the water column; ii-parametrization, defines the minimal set of IOPs whose values completely characterize the observed radiance; iii-inversion, derives the values of IOPs, and hence water quality variables, from radiometric data.
Reliable methods for uncertainty quantification of earth observation (EO) products of IOPs are important for sensor and algorithm validation, assessment, and operational monitoring.High accuracy in both observations and algorithms may reduce considerable ranges of errors.EO derived IOPs, however, have an inherent stochastic component.This is due to the dynamic nature of aquatic biogeophysical quantities, intrinsic fluctuations, model approximations, correction schemes, and inversion methods.Due to stochasticity of the measurements, as well as model approximations and inversion ambiguity, the retrieved IOPs are not the only possible set that caused the observed spectrum (Sydor et al., 2004).Instead, many other IOPs sets may be derived.Each of these sets has an unknown probability of being the derived product.The probability distribution of the estimated IOPs provides, therefore, all the necessary information about the variability and uncertainties of derived IOPs.
Generally, uncertainty assessment of EO-data falls under one of two methods, namely analytical deterministic or stochastic methods.
Deterministic methods are based on gradient techniques and have been used to asses the uncertainty of IOPs as derived from EO-data.Duarte et al. (2003) analyzed the sensitivity of the observed remote sensing reflectance due to variable concentrations of water constituents.Maritorena & Siegel (2005) employed a deterministic technique for consistent merging of different products using their uncertainties.Wang et al. (2005) performed a detailed study on the uncertainties of model inversion related to fluctuations in each of the IOPs and their spectral shapes.Salama et al. (2009) studied the uncertainty of model-inversion using the gradient-based method.They found that the derived IOPs are linearly related to their errors.Lee et al. (2010) used analytical derivative of the quasi-analytical algorithm (Lee et al., 2002, QAA) to estimate the uncertainty of IOPs as derived from QQA.On the other hand, Salama et al. (2011) developed a gradient based method to estimate the accuracy of a specific model-parameterizations setup.The

Ocean color model inversion
Remote sensing reflectance, the ratio of radiance to irradiance, above the water surface Rs w can be related to the inherent optical properties (IOPs) using the ocean color model of Gordon et al. (1988): (1) Where Rs w (λ) is the remote sensing reflectance leaving the water surface at wavelength λ; g i are constants taken from Gordon et al. (1988); t and n w are the sea−air transmission factor and water index of refraction, respectively.Their values are taken from literatures (Gordon et al., 1988;Lee, 2006;Maritorena et al., 2002).The parameters b b (λ) and a(λ) are the bulk backscattering and absorption coefficients of the water column, respectively.The light field in the water column is assumed to be governed by four optically significant constituents, namely: water molecules, phytoplankton green pigment chlorophyll-a (Chl-a), colored dissolved organic matter (CDOM) and detritus/suspended particulate matter (SPM).The absorption and backscattering coefficients are modeled as the sum of absorption and backscattering from water constituents: b b (λ)=0.5bw (λ)+ηb spm (λ). ( Where the subscripts on the right hand side of equations ( 2) and (3) denote water constituents: water w; phytoplankton green pigment ph; lumped absorption effects of CDOM and detritus dg and suspended particulate matter spm.η is the backscattering fraction, its value is estimated from Petzold's "San Diego harbor" scattering phase function as η ∼ 0.018 (Petzold, 1977).
The absorption and scattering coefficients of water molecules, a w (λ) and b w (λ), are assumed to be constant.Their values are obtained from Pope & Fry (1997) and Mobley (1994), respectively.The total absorption of phytoplankton pigments a ph (λ) is approximated as in Lee et al. (1998), a ph (λ) ≃ a 0 (λ)a ph (440)+a 1 (λ)a ph (440) ln a ph (440), where a 0 (λ) and a 1 (λ) are statistically derived coefficients of Chl-a, their values are taken from Lee et al. (1998).
The absorption effects of detritus and colored dissolved organic matter (CDOM) are combined due to the similar spectral signature (Maritorena et al., 2002) and approximated using the model of Bricaud et al. (1981), where s is the spectral exponent of combined effects of detritus and CDOM.The scattering coefficient of SPM b spm (λ) is parameterized as a single type of particles with a spectral dependency exponent y (Kopelevich, 1983): Equation ( 1) is inverted to derive five parameters from the IOCCG data set and three parameters from the NOMAD data set.The derived parameters are called the set of IOPs and expressed in a vector notation as iop.The exponents s and y are assumed to be unknown (Salama et al., 2009) and are derived from the IOCCG data set as: iop = a ph (440), a dg (440), b spm (550), s, y .
The numerical inversion is carried out using the constrained Levenberg-Marquardt Algorithm (LMA) (Press et al., 2002), where the constraints are set such that they guarantee positive and physically meaningful values: between 0 and 100 m −1 for a ph (440), a dg (440) and b spm (550), between 0 and 2.5 for y and between 0 and 0.03 for s.Optimization is started using the initial values of Lee et al. (1999) and s = 0.021 nm −1 and y = 1.7.Maximum number of iteration is set equal to 100.

Description
The uncertainty in the derived IOPs is attributed to the infinitesimal change of radiance in equation ( 1) as, where ∆Rs w (λ) represents the radiometric uncertainty at the wavelength λ; w ph , w dg ,w spm are the partial derivatives of Rs w with respect to the derived IOPs.Equation ( 8) represents an over determined linear set of equations that can only be solved if the radiometric uncertainty is known in at least n wavelengths, with n being the number of derived IOPs.

231
Current Advances in Uncertainty Estimation of Earth Observation Products of Water Quality

www.intechopen.com
Analytical expressions of partial derivatives in (8) are listed hereafter.To simplify the notations let us define the ratio ω as, where c b = b b (λ)+a(λ).The partial derivative w ph is, where ζ ph is the spectral dependency of Chla, The parameters j i are j 1 = −g 1 and j 2 = −2g 2 c b .The term w dg is expressed as, The partial derivative w spm is expressed as, where Based on the above theoretical formulation in equation ( 8), Lee et al. (2010) obtained the uncertainty of IOPs using the quasi analytical algorithm (Lee et al., 2002) and a prior information on the radiometric errors.Salama et al. (2011), on the other hand, proposed a method that produces a single (or ensemble) uncertainty measure for the collective errors in the derived IOPs relative to the radiometric uncertainty without the need for model inversion or prior information on the radiometric errors.In addition, the method provides the optimum accuracy which can be achieved by a model-parametrization setup.The method of Salama et al. (2011) is self-contained and is directly applicable to existing satellite based IOP products, we therefore, brief this method hereafter.

Ensemble uncertainty of IOPs
Applying Taylor series approximation of the second moment on equation (8) gives: Where σ 2 r (λ) is the radiometric variance and σ 2 ph (440), σ 2 dg (440), and σ 2 spm (550) are the variances of the derived IOPs.The covariance terms in equation( 14) is assumed to be zero, i.e. the IOPs are mutually independent.Knowledge on the radiometric uncertainty is now 232 Earth Observation www.intechopen.comavoided by dividing both sides of equation ( 14) by the radiometric variance, with The ensemble uncertainty of IOPs per radiometric error, Ψ(λ), is derived from equation ( 15) by normalizing both sides by the squared sum of partial derivatives and taking its square-root: . ( 16)

Detailed uncertainty of IOPs
Based on equation(8), Bates & Watts (1988) devised an elegant method to quantify the uncertainties for each derived IOPs as, Where IOP i± is the upper "+" and lower "-" bounds of the derived IOP; W is the matrix of partial derivatives; σ is the standard deviation of residuals between measured and model best-fit radiances; t(N − m, α/2) is the upper quantile for a Student's t distribution with N − m degrees of freedom.N is the number of bands and m is the number of unknowns.R is the upper triangle matrix of QR decomposition of the jacobian matrix.equation ( 17) has widely been used to estimate the error of derived IOP (Salama et al., 2009;Van Der Woerd & Pasterkamp, 2008).The derivative term in equation ( 17), can be approximated as being the gradient of equation ( 1) with respect to the derived IOPs and is computed for model-best-fit to the observation.This approximation is derived as follows.
Observed remote sensing reflectance can be approximated as being the sum of the model best-fit Rs m (λ) and its deviations from the observed one ǫ(λ): The term Rs m (λ) is obtained from fitting the model in equation (1) to the radiometric observation of ocean color or/and field sensors.model goodness-of-fit, measurements and atmospheric noises.For simplicity this term is assumed to be nearly independent the derived IOPs.The derivative of (18), with respect to the derived values, can then be written as: By definition of the least square minimization that was used to derive model-best-fit Rs m (λ), we have: Equation ( 19) can then be reduced to: The simplification in equation ( 21) implies that the gradient of measured remote sensing reflectance can be approximated by the gradient of the model in (1) which can easily be computed as in equation ( 21).

Description
In this section we summarize the method of Salama & Stein (2009) as it is the only stochastic method published so far in the field of ocean color.
Salama and Stein used prior information to obtain plausible ranges of the IOPs.These ranges are used in a log-normal distribution to generate a first-estimate of the probability distribution (PD) of the IOPs.This first-estimate PD is called the prior PD of the IOPs.The method, explained hereafter, uses the prior PD to converge to a "posterior" probability distribution that better describes the IOPs.
Prior information is obtained from known radiometric errors in Rs w and model-inversion intrinsic errors.Radiometric errors are: (i) noise equivalent radiance of the sensor and (ii) error in aerosol optical thickness.Sensor equivalent radiance is known from sensor specifications and post-launch calibrations.Model approximation and inversion-accuracy can be quantified by evaluating the performance of the employed ocean color model against measurements and radiative transfer simulations.Atmospheric error, due to variation in aerosol optical thickness, can be evaluated from available measurements or by using standard atmospheric correction models.The error estimate algorithm will follow sequential steps as detailed hereafter.
An initial estimate of the confidence interval around water remote sensing reflectance can be computed using the method of (Bates & Watts, 1988, pp.59, cf. 1.36 ) or available knowledge on plausible fluctuations for model, noise and atmospheric residual respectively.The upper and lower bounds of this interval are then inverted to derive the corresponding two sets of IOPs iop u , iop l .These sets with the derived iop obs from the water remote sensing reflectance, hereafter will be called the IOP-triplet: (iop l , iop obs , iop u ) and denoted as ω.The value log iop obs is assumed to approximate the mean of a first-estimate, i.e. prior, probability distribution (PD) of IOPs in the logarithmic space.The prior PD is first elicited using the IOP-triplet and prior knowledge on the log-normal shape of the IOPs as explained in

235
Current Advances in Uncertainty Estimation of Earth Observation Products of Water Quality www.intechopen.comsection (4.2).The posterior probability distribution, or our gain in information, is then inferred by maximizing the expected utility (Bernardo, 1979;Carlin & Polson, 1991) as explained in section(4.3).

Prior probability distribution
Estimating the IOP-triplet, iop l , iop obs and iop u , is the first step towards deriving the prior probability distribution of the IOPs.The use of flat or improper priors, e.g.uniform distribution, may invalidate the derivation of the posterior probability (Goutis & Robert, 1998).According to the maximum entropy principle (Jaynes, 1957a;b) a proper prior probability distribution should have the maximum entropy provided by the IOP-triplet.However applying the maximum entropy principle on the information provided by the IOP-triplet will give the probability values of iop l , iop obs and iop u but not the whole probability distribution P(iop); for more detail one may consult Jaynes (1968).To overcome this limitation, in data values, we introduce the following method to elicit the prior distribution of IOPs assuming that they are log-normally distributed.The log-normal assumption is based on Campbell's work (Campbell, 1995) who pointed out that, in general, marine bio-geophysical quantities follow a log-normal distribution i.e. their log transform has a Gaussian distribution.
The IOP-triplet is first transformed to the log space, allowing us to use a Gaussian distribution to simulate the PD of IOPs.Second we assume that log iop obs approximates the mean of the prior PD of the IOPs.The Gaussian distribution of the IOPs can be standardized to a N(0,1) distribution, i.e. normal distribution with zero mean and unity standard deviation.
The standard Gaussian variate for log iop u is, where α u is a sample drawn from the N(0,1) that corresponds to iop u .The parameters log iop obs and σ are the expectation and the standard deviation of the population.From equation ( 22) and the second set in the IOP-triplet iop l we can establish the ratio, and for convenience we set log iop u > log iop l .The standardization of the IOPs distribution allows us to use the N(0,1) random number generator to simulate values of α as in equation ( 22).The ratios of these random values are also computed and compared to the ratio of the IOP-triplet in equation ( 23).The best fit allocates the two values α u and α l , hence the standard deviation of the prior distribution can be computed from equation ( 22).The prior probability distribution of the IOPs, is now known: N(log iop obs ,σ), i.e. a Gaussian distribution with log iop obs mean and σ standard deviation.
Random values (1000) are generated from the N(0,1) distribution such that they satisfy an imposed acceptance-rejection condition.This condition requires that the ratio in equation ( 23) defines a unique ordered pair of α.This is to enable the use of a simple searching method with a fast convergence to the best-fit ratio.The uniqueness in this sense implies that the squared difference between the computed ratio, from the IOP-triplet, and the best-fit is a global minimum resolvable by the searching method and the used computer processor.Three look-up tables (LUTs) are then created from the generated values.These LUTs correspond to 236 Earth Observation www.intechopen.comthe following three scenarios: The generated N(0,1) values are, first, subdivided into two sets containing positive and negative values.The ratios of the first and second LUTs are, then, computed from the ordered descending sets as; x i /x i+1 .The third LUT is generated from all possible combinations of the unordered positive and negative sets.This will results in ratio values between 0 and 1, > 1 and < 0 for the first, second and third LUT respectively.The ratio in equation ( 23) is first estimated from IOP-triplet.Based on the values of this triplet (equation 24) a lookup table is selected and searched to find the best-fit value to the computed ratio (equation 23).This best-fit is found either by direct search or interpolated.One of the corresponding pair is then used in equation ( 22) to compute the standard deviation of the prior PD P(iop).

Posterior probability distribution
In section (4.2) we derived a proper prior distribution of the IOPs.This first-estimate, i.e. prior distribution, is converged to a posterior distribution that better describes the IOPs using the concept of Entropy.Entropy is a numerical measure of error associated with probability distribution of derived IOPs or any hydrological parameter (Singh, 1998).For a population with N sets of IOPs it is expressed as the Shannon entropy (Shannon, 1948): where P(iop) is the prior probability distribution (PD) of the derived set of IOPs iop.
If we design a function D that measures the information, e.g.equation ( 25), between the prior and the posterior PD, then we can derive the posterior PD such that it maximizes the expected information to be gained in D (Bernardo, 2005;Christakos, 1990).In other words, maximizing the function D will maximize the gained information from the posterior PD (Bernardo, 1979).The Kullback-Leibler divergence (Kullback & Leibler, 1951), or cross-entropy, belongs to this type of utility functions (Johnson & Geisser, 1985).It measures the divergence between the posterior P(iop|ω) and the prior P(iop) probability distribution as: where P(iop|ω) is the posterior probability of iop given the IOP-triplet ω.Equation ( 26) can be rewritten in view of equation ( 25) as: where H{P(iop|ω),P(iop)} is expressed as: Maximizing the cross-entropy in equation ( 26) or the corresponding expression in ( 27) is equivalent to minimizing the entropy (uncertainty) of the posterior probabilities distribution, i.e. maximizing gained information.The errors can then be estimated from the reconstructed posterior probability distribution of IOPs P(iop|ω).
The posterior probability distribution is inferred by maximizing the utility function, i.e.Kullback-Leibler divergence (equation 26).The maximum is found by iteration through a sequential updating of the posterior using the prior parameters mean μ and variance σ 2 (Rubinstein & Kroese, 2004).The corresponding log-normal mean m and variance v are computed as Kendall & Stuart (1987): The following steps describe the algorithm, as implemented, to derive the posterior PD P(iop|ω): 1. From the water remote sensing spectrum estimate the initial radiometric confidence interval using the method of (Bates & Watts, 1988, pp.59, cf. 1.36) or prior information on atmospheric and noise-induced radiometric fluctuations.2. Invert the ocean color model in equation ( 1) to derive the IOPs from the water remote sensing spectrum and the upper and lower bounds.This will results in three sets of IOPs: iop l , iop obs , iop u ; IOP-triplet.3. Based on the order of this IOP-triplet allocate the suitable LUT using equation ( 24). 4. Search for the best-fit ratio calculated from equation ( 23). 5. Use equation ( 22) to estimate the standard deviation of the prior PD. 6. Use the standard deviation and log iop obs to generate the prior PD. 7. Use initial values of the mean and standard deviation to generate n Monte Carlo samples of PD. 8. Select the population that have the maximum Kullback-Leibler divergence (equation 26), and update the initial values.9. Repeat step 7 to 8 till convergence.10.Update the prior PD with the resulting posterior PD (from the pervious step: 9), and iterate steps 7 to 10 till convergence.
The convergence is defined by a threshold as follow.Keep track of the best ten candidates which maximize equation (26).The system converges if the variance of these ten values is less than 10 −4 .

Inter-comparison between deterministic and stochastic methods
The inter-comparison between the deterministic method, described in Section (3), and the stochastic method, described in Section (4), is carried out using two data sets.The first, is radiative transfer simulations of synthetic IOPs obtained from the International Ocean Color Coordination Group (IOCCG), report-5 (Lee, 2006, IOCCG data set).The second consists of concurrent observations from the Sea viewing Wide Field-of-view Sensor (SeaWiFS) and measured inherent and apparent optical properties, retrieved from the NASA bio-Optical Marine Algorithm Data set (NOMAD) Version 1.3 (Werdell & Bailey, 2005, SeaWiFS matchup data set).

IOCCG
IOCCG data set (Lee, 2006) of synthesized IOPs and their radiative transfer simulations at 30 • sun zenith angle are used to inter compare the results of the deterministic and stochastic methods.IOCCG simulated spectra, between 400 nm and 750 nm at 10 nm interval, are inverted using the ocean color model in equation ( 1) to derive five variables.These variables are: Chlorophyll-a absorption at 440 nm a ph (440), detritus and CDOM absorption at 440 nm a dg (440) and their spectral dependency s, SPM scattering at 550 nm b spm (550) and SPM spectral dependency y, as shown in equation ( 7).
The standard deviation of the posterior PD represents the error/confidence of the derived value iop obs .The deviation of the posterior PD from known IOPs is measured using root-mean-square of errors (RMSE).These two values, RMSE and standard deviation, are related through the bias, i.e the actual difference between derived and measured IOPs. Figure (2) shows estimated errors, expressed as standard deviation using equation ( 30), against the known root-mean-square of errors (RMSE).The actual RMSE is estimated from the posterior PD and the known IOPs.The reproduced errors for the IOPs other than a ph (440) have a high accuracy with r 2 values between 0.77 and 0.96.Estimated errors of a ph (440) have the lowest r 2 and n values.It is worth noting that the determinacy method of Bates & Watts (1988) generally underestimates model-errors of the IOPs with lower r 2 values than the presented stochastic method.This is apparent at an almost threefold difference for the error values of a ph (440).On the other hand, the stochastic method has a tendency to overestimate the errors of the IOPs with a better fit and improved capability, in the sense that it can be applied to populations of any bio-geophysical variable.

NOMAD
Due to the limited number of available visible bands in this data set we reduced the number of unknowns to three only.The first three IOPs in equation ( 7) are derived from SeaWiFS spectra using the ocean color model (equation 1) and the constrained LMA technique.The values of s and y are set to 0.021 nm −1 and 1.7 respectively.The actual RMSE values are computed from the posterior PD and measured IOPs.The total error on derived IOPs is estimated by applying the stochastic method using (Bates & Watts, 1988, pp.59, cf. 1.36) radiometric confidence interval.The estimated errors are expressed as standard deviation using equation ( 30) and plotted against RMSE values in figure (3).The reproduced total error values are strongly correlated to the known RMSE values with r 2 between 0.67 and 0.9 and >90% of valid retrievals.Estimated errors from the deterministic technique (Bates & Watts, 1988), however, did not correspond to the actual values of RMSE.
Errors are computed for the ocean color model and SeaWiFS visible bands centered at [412,443,490,510,555,670] nm.The average values of the derived standard deviation are 1.7802, 1.1431 and 1.6177 m −1 for a ph (440), a dg (440) and b spm (550), respectively.

Uncertainty sources 6.1 Description
The total remote sensing reflectance received at the sensor altitude can be written as the sum of several components (Gordon, 1997): The subscript of the remote sensing reflectance Rs represents the contribution from: (i) the atmosphere (path), i.e. air molecules and aerosol multiple scattering; (ii) sea-surface (sfc); and (iii) water (w).T(λ) is the diffuse transmittance.
The contribution of air molecules, i.e. the Rayleigh scattering, to the atmospheric path is well described in terms of geometry and atmospheric pressure (Gordon et al., 1988).The contribution of sea-surface reflectance Rs sfc can be estimated using the probabilistic formulations of Cox & Munk (1954) and ancillary data on wind field.Gaseous transmittance can be calculated from ancillary data on ozone and water vapor concentrations using the transmittance models of Goody (1964) and Malkmus (1967).For viewing angles < 60 • the diffuse transmittance T is weakly dependent on aerosol and can be approximated following Gordon et al. (1983).Following the aforementioned approximations will basically leave two unknowns; the aerosol and the water remote sensing reflectance.In other words, the errors in Rs w can be attributed to errors in aerosol estimation and any noise in the sensor, i.e. noise equivalent radiance (NER).
Radiometric errors in Rs w , beside to model-inversion intrinsic errors, will accumulate and propagate to the IOPs during the retrieval.The total error of the derived IOPs can therefore be decomposed into three major components, namely model-inversion error, sensor noise and error in aerosol estimation.These errors are originated by various mechanisms during the processing chain of ocean color data as explained hereafter.
Each error component, x, will be expressed as the variance σ 2 x of IOPs caused by this error x.The subscript x will be replaced by inv, ner and a to represent the contribution of model-inversion, noise equivalent radiance and aerosol, respectively.

Model-inversion error, σ 2 inv :
The employed approximations in the forward-model (equation 1) may not precisely describe the optical processes that have caused the observed signal (Zaneveld, 1994).Moreover, the numerical technique used for inversion provides an ambiguous solution, i.e. the derived IOPs are not unique (Sydor et al., 2004).These assumptions and ambiguity will generate error that is, at the one hand, inherent to the employed ocean color forward model and, on the other hand, dependent on the accuracy of the inversion scheme which could be related to the optical complexity of the water.Model-inversion error is quantified as a lumped sum of errors due to the approximation in (1), the parametrization of IOPs and inversion and abbreviated as model error.

Noise equivalent radiance, σ 2 ner :
Noise equivalent radiance (NER) depends on sensor specifications and performance over time, i.e. sensor degradation.This fluctuation could either increase or decrease the observed remote sensing reflectance and could also be wavelength dependent or random.The effects of NER is inversely proportional to the value of signal-to-noise ratio.Sensor degradation, i.e. sensitivity losses over time, will cause decrease in the signal-to-noise ratio of the sensor leading to low signal reading.Low signal can also be observed over clear water at the near infrared part of the spectrum or over turbid water, with high CDOM, detritus and Chl-a contents, at the blue part of the spectrum.The propagated error from NER to IOPs will therefore be dependent on sensor specification, sensor degradation over time, water turbidity and observing wavelength.Atmospheric correction errors are, generally, caused by unknown aerosol type and optical thickness (AOT).The residual signals from atmospheric correction will have spectral and spatial dependency.The spectral dependency is due to the error about the aerosol type e.g.absorbing aerosol, while the spatial dependency is, on the one hand, related to the error about AOT spatial variations and, on the other hand, to water turbidity (Hu et al., 2004).It is assumed that aerosol optical thickness has a higher spatial variability than aerosol type, so that aerosol type can be assumed to be known and homogenous.Within the validity of this assumption, the residual signals from atmospheric correction will be caused by errors in estimating the aerosol optical thickness.

Decomposition
The total error of the derived IOPs, expressed as the variance σ 2 t σ 2 inv , σ 2 ner , σ 2 a , is thus described as a function of the three error components, σ 2 inv , σ 2 ner and σ 2 a .Assuming that this function is continuous in its variables, we can approximate it by a first order Taylor series as: where, σ 2 t0 is the value of the function σ 2 t (0, 0, 0).According to the assumption that the total error is caused by three components, the value of σ 2 t0 is negligible, i.e. σ 2 t0 ≃ 0. In other words, if we have perfect measurements, accurate atmospheric correction and exact model parameterizations and inversion then the total error on the derived IOPs will be negligible.The total error of the derived IOPs can thus be approximated as a weighted sum of the individual error components as: where the weights w inv , w ner , and w a are the partial derivatives in equation ( 32).The functionality in σ 2 t , however, is commonly unknown and it is therefore difficult to find proper estimates of the weights w inv , w ner and w a .An intuitive approach would be setting all the weights in equation ( 33) to unity and check its validity: Figure ( 4) depicts the relationship between the sum of the righthand side of equation ( 34) and the total error on the derived IOPs.On the X axis is the total error of the IOPs as calculated from all possible error sources σ 2 t .We then calculated each error component apart and summed their variances in the Y axis as: σ 2 a + σ 2 ner + σ 2 a .As anticipated from equation ( 33) there is a linear relationship between the actual variance and the linear sum of individual variances with R 2 values above 0.75 for the absorption coefficients of Chl-a and detritus-CDOM.The value of R 2 decreases to 0.69 for SPM scattering and 0.64 for the total absorption.The dispersion value as measured with RMSE is large for all IOPs.The results in figure (4) indicate that the linear sum in equation ( 34) is an acceptable approximation to the total variance.Due to the large values of RMSE in figure (4), the computed relative contribution should be treated with caution.
While model-induced error can directly be estimated from the techniques described in Brad (1974) and Bates & Watts (1988), noise and atmospheric-induced errors should be inferred from the available information.This information forms the prior knowledge that we will use in the following section to derive the error of the IOPs.Prior information is obtained from known sensor's noise, variation in aerosol optical thickness and ocean-color model's approximations and inversion accuracy.

IOCCG
The noise is estimated based on NER values of the Medium Resolution Imaging Spectrometer (MERIS) (Doerffer, 2008;Hoogenboom & Dekker, 1998).The variation in aerosol optical thickness (AOT) is set to be ±0.02.This value is estimated from the variation of recorded aerosol optical thickness by a newly calibrated sunphotometer (CIMEL) and cloud free condition (Holben et al., 2000).The values of aerosol optical thicknesses are obtained from sunphotometer measurements situated at (51.225 N, 2.925 E) at the 8th of June 2006.The atmospheric paths are estimated with radiative transfer computation (Vermote et al., 1997) using maritime aerosol model with a nadir looking sensor at 30 • sun-zenith and 203 • sun-azimuth angles.
The relative contribution of model, noise and atmospheric errors are shown in table (1) and quantified for each of the derived IOP as follow.First we computed the total error, i.e. the total error in Rs w (λ) is due to aerosol estimation and sensor noise, inversion error will add 244 Earth Observation www.intechopen.comup during the inversion.The same step is repeated for each error source in three steps: (i) model error is estimated from the error-free Rs w (λ); (ii) atmospheric-induced error σ 2 a is computed from Rs w (λ) that contains errors due to aerosol estimation only; (iii) noise error is calculated from Rs w (λ) that contains sensor noise only.Note that model error will add up during the inversion in the last two steps.Now we can use equation ( 34) to estimate the relative contribution of each error component to the total error of the IOPs.
Errors due to atmospheric correction are the major source of errors in the derived IOPs.Imperfect atmospheric correction, due to the variability of aerosol optical thickness, is responsible for more than 50% of the total error and up to 82%.One fifth of the total errors on derived IOPs (except for the SPM scattering: one tenth) is attributed to noise-error.Model-error has the lowest contribution (≈7%) to the total error on derived b spm (550) values, but it has a significant contribution (≈ 16%) to y.This can be attributed to the assumed parametrization.On the one hand, the absorption of other constituents than water molecules is negligible at the near infrared (NIR) which will cause stability (one-to-one relation) in the derived SPM scattering coefficient, leading to a significant contribution from the atmosphere at the NIR region.On the other hand, the error in b spm (λ) will decrease towards the NIR region due to the assumed exponential spectral dependency.In general, model-induced errors are large for the spectral shape coefficients y and s.Note that the spectral shape of chlorophyll-a absorption is imbedded in coefficient a ph (440  1.The average relative contribution (%) of error components on IOCCG data set.

NOMAD
The total error on estimated IOPs from the NOMAD data is derived from the values presented in figure (3).Model induced errors are subtracted from the total error using equation ( 33) to deduce atmospheric and noise-induced errors.The results are shown as percentages in table (2).Main uncertainty is due to atmospheric and noise-induced errors for a ph (440) and b spm (550), while model inversion is the main source of error to a dg (440) in this data set.These results are within the validity of the linear assumption expressed in equation ( 33

Model-sensor error table
The linear sum of individual variances in equation ( 34) can describe about 70%, the value of R 2 , of the total variance of the IOPs.This linearization of the total variance is a simple yet effective approach.It allows us to estimate the relative contribution of the different error components to the total error budget on the IOPs.The relative contribution of model, noise and atmospheric errors to the total error budget using IOCCG data set are 20-40%, 10-25% and 40-80%, respectively.Model-induced errors, due to approximation and inversion, are inherent to the derived IOPs and inversely proportional to model-inversion degree-of-freedom, while atmospheric-induced errors are the major contributor to the total error budget on IOPs.These results are for assumed levels of noise and atmospheric fluctuations.This suggests that error table can be generated for specific model, sensor and range of IOPs.This model-sensor error table can serve as a benchmark to estimate the atmospheric-induced errors in the derived IOPs.The merit of this argument is based on the fact that the computations of model and noise-induced errors can be quantified using water radiative transfer simulations, for a specific range of IOPs, and known sensor's NER.The magnitude of these errors are in principle known for the ocean color model and the used sensor.An example of such a table for the MERIS sensor and the ocean color model is shown in table (3).This table is computed from table (1) for the MERIS visible bands centered at [412,443,490,510,560,620,665,708,778] nm, i.e. we simply reduced the spectral bands of IOCCG data set to fit those of MERIS.Table (3) shows that the reduced number of spectral bands for MERIS setup has increased model contribution to the total error approximately two fold.This will reduce noise and atmospheric contribution to the total error, since the relative contributions of all error components should sum to 1.Note that for weak radiometric signals, the lower bound might end up to negative values which will lead to further reduction in the number of bands (negative values are set to zero).This approach is demonstrated for ocean color observations obtained from NOMAD data set.Model and noise-induced errors are simulated from the IOCCG data set and subtracted from the total error of IOPs estimated from the NOMAD data set.The simulation is carried out simply by selecting IOCCG wavelengths that correspond to NOMAD spectral set-up.The simplicity of this approach can pose a limitation on the accuracy of equation ( 34).On the one hand, the method shows that model approximation and inversion are main contributors, ≈57%, to the total error of a dg (440).On the other hand, the presented stochastic method quantified these errors most efficiently.Atmospheric and noise-induced errors are significant for a ph (440) and b spm (550).This may suggests that model-induced errors are better quantified with the current method.

Spectral propagation of errors and error correlation
The presented errors of IOPs were for two wavelengths: 440 nm for the absorption coefficients and 550 nm for the scattering coefficient, as defined by equation ( 7).We can use the parameterizations in equations ( 4), ( 5) and ( 6) to derive analytical description of error propagation to other wavelengths.Here below we provide an analytical derivation of error propagation and numerical examples for two wavelengths one at the blue, 400 nm, and the other at the red, 680 nm.
The errors of the IOPs will propagate to shorter and longer wavelengths following the parameterizations in equations ( 4), ( 5) and (6).For example, the error in b spm (550) has two components; one in b spm (550) itself and the other in the spectral shape y.Using the parametrization in equation ( 6) we will have: Carrying the derivation of the palatial derivatives, equation ( 35) can be written as: In this exercise we will neglect the error in the wavelength, i.e. ∆λ ≈ 0 and we will show that the derivative ∂/∂λ is negligible.
Equations (37,38) show that the error in SPM scattering coefficient at the blue wavelength is larger than that at the red wavelength if the relative error in the scattering coefficient satisfies the condition: In a similar approach we can quantify the propagated errors of a dg (440) to other wavelengths: If we assume the value s = 0.021 nm −1 and take our reference bands to be the blue (400 nm) and red (680 nm) wavelengths we will have, with λ in meter: ∆a dg (400)=2.316∆adg (440)+92.654× 10 −9 a dg (440)∆S (41) ∆a dg (680)=6.47 × 10 −3 ∆a dg (440) − 1.554 × 10 −9 a dg (440)∆S (42) The wavelength variation term ∂/∂λ is also negligible.It takes the values −4.86 × 10 −11 a dg (440)∆λ and −1.36 × 10 −13 a dg (440)∆λ for the blue and the red wavelengths, respectively.The error at the blue will be larger than that at the red if the relative error of a dg (440) satisfies the condition (from equations 41 and 42): The parametrization of Chl-a absorption is based on the tabulated values a 0 and a 1 , see equation( 4).These tabulated values are taken to be constant per wavelength, i.e. a ph (λ) is function of a ph (440) only.The error in a ph (440) will propagate to other wavelengths following the derivative of equation ( 4): For the two reference bands, 400 nm and 680 nm, we will have: From equations (45,46) it can be shown that the error at the blue band is larger than that at the red if the following condition is satisfied: log a ph (440)∆a ph (440) < −1.562 The analytical expressions in equations ( 36), ( 40) and (44) show that the errors are related to absolute values of the IOPs.Therefore, the three error components are expected to be correlated to water turbidity, and hence to each others.The results of the numerical examples also demonstrate that the errors of b spm (λ) and a dg (λ) will be larger at the blue than that at the red if the relative errors of b spm (550) and a dg (440) satisfy equations ( 36) and (40), respectively.Whereas the error in a ph (440) will propagate to other wavelengths following equation (44) and will be larger at the blue if the condition in (47) is satisfied.

Advantages and limitations of error estimation methods
Estimated errors from the deterministic method (Bates & Watts, 1988) did not correspond to the actual values of RMSE.This is due to the atmospheric and noise radiometric fluctuations.These fluctuations are imbedded in the observed signal and do not vary with IOPs values, i.e. different response function.Their large fluctuations may cause an ill-conditioned Jacobian matrix that produces erroneous estimates, see (Bates & Watts, 1988, pp.59, cf. 1.36 ) application produces a single (or ensemble) uncertainty measure for the collective errors in the derived IOPs relative to the radiometric uncertainty without the need for model inversion or prior information on the radiometric errors.
Error decomposition exercise shows that atmospheric and NER induced errors can be better quantified when prior knowledge is available.This is important for ocean color band ratio or single band algorithms, e.g.(Austin & Petzold, 1981;Salama et al., 2004).These algorithms are empirical in nature, i.e.Jacobian matrix is not available.In this case, deterministic methods to derive the error are not applicable.In contrary, the presented stochastic method is generic and can be applied to quantify the error of any derived bio-geophysical parameter regardless of the used derivation method.This is true if, beside to the derived quantity, two other values are known a priori so the IOP-triplet can be constructed.
The prior values were inferred from the quantiles of the populations.In practice this information is not available but it could be estimated from historical measurements or high temporal observations.The later, high temporal sampling, can be realized using sensors on board of geostationary satellites to quantify marine bio-geophysical parameters.For instance, the visible band of the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board of the METEOSAT second generation satellite (MSG) can be used to quantify the concentrations of SPM (Neukermans et al., 2008).With MSG 15 minutes of repeated sampling cycle, the stochastic method can be applied on three consecutive acquisitions, i.e. each 45 minutes, to produce SPM concentration and related-error maps.This error map provides vital input to the recently developed SPM assimilation model (Eleveld et al., 2008).Moreover it can be used, as weights, for ocean color products merging (Pottier et al., 2006).This generality aspect of the presented stochastic method expands its applicability to different fields other than ocean color.For example, Velde van der et al. ( 2008) developed a basis for Synthetic Aperture Radar (SAR)-based soil moisture downscaling methodologies.
One limitation of the presented stochastic method is the choice of the acceptance-rejection method.Although it facilitates the search for a unique pair of N(0,1) values, the derived σ become sensitive to the ratio in (23), i.e. sensitive to the lower and upper pair (iop u ,iop l )in the IOP-triplet.This may caused the 7∼10% failure to reproduce the values of the standard deviation.This can be attributed to the small values of α ≪ 1 which produce large values of σ.These large values will further be magnified by equation (30).
Using equation (Bates & Watts, 1988, pp.59, cf. 1.36) to estimate the total error as a linear sum of all other error components is another limitation.Atmospheric or noise radiometric fluctuations can be interpreted, by model inversion, as high/low IOPs values with high goodness-of-fit.Using the same reasoning, bad fit to very complex signal (turbid water with high SPM, CDOM and Chl-a contents) can be attributed to atmospheric and sensor noise errors, although the observed signal might be error-free.
Model-sensor error tables were simulated from IOCCG data set without accounting for sensor's band width and response function.A more detailed simulations that includes band width, response function of the sensor and a specific range of the IOPs should be carried out to establish a more accurate model-sensor error tables.
Although we showed that equation ( 34) is an acceptable approximation to the total variance, the computed relative contribution of errors should be treated with caution.

Summary and future developments
In this chapter we reviewed the recent advances in uncertainty estimation of the earth observation products of water quality.Both deterministic and stochastic methods are presented and their results are inter-compared.The stochastic method is more appropriate to estimate actual errors of ocean color derived products than the deterministic methods, however, it is still limited to few studies and as the deterministic approach requires prior information.The uncertainties could be decomposed only if additional information is provided a priori.Using a simple exercise it was shown that atmospheric-induced errors are major contributors to the total error of IOPs whereas model-induced errors are inherent to the derived IOPs depending on the used derivation method and number of spectral bands.
The error in this chapter was estimated as the difference between ground truth measurement and satellite derived products.Direct matching between earth observation data and just above the water field measurements imbed, however, an inherent scale difference.This scale difference between in-situ observation and a pixel of ocean color satellite is at least three to four orders of magnitude for nadir match-up sites and much larger for off-nadir ones.This huge scale difference, means that point measurement is sampling a tiny fraction of the water body which is observed by a satellite pixel.Few studies were carried out to address the scale difference between point and aerospace measurements directly.Most of these studies have used re-sampling to smooth out the scale differences in the match-up sites, see (Bailey & Werdell, 2006;Bissett et al., 2004;Harding Jr. et al., 2005;Hu et al., 2000).For example, Hyde et al. (2007) applied a correction algorithm to SeaWiFS products of chlorophyll-a to overcome the mismatch which was partially due to sampling size differences.Although this assumption of spatial homogeneity have resulted in good matches for most open ocean matchup data (Carder et al., 2004;Garcia et al., 2005;Karl & Lukas, 1996;McClain et al., n.d.), it lowers the percentage of usable match-up points considerably (Hooker & McClain, 2000;Mélin et al., 2005) and should be avoided for productive waters (Chang & Gould, 2006;Darecki & Stramski, 2004;Harding Jr. et al., 2005).Salama & Su (2010;2011), used the differences between the earth observation products and in situ data to quantify the sub satellite pixel spatial viabilities using both the deterministic and stochastic methods, respectively and neglecting the error.In principle the mismatch between earth observation derived products and in situ measured quantities is attributed to the scale difference and errors due to noise, correction and retrieval accuracy.Current uncertainty estimation methods do not consider the spatial dependency of errors and their relationships to the actual distribution of IOPs.Understanding the spatial characteristics of errors is necessary to resolve the smallest sub-scale variability of the IOPs.This aspect should be investigated in the future to define spatial-thresholds of measurable physical processes based on their errors.Moreover, the dependency of both deterministic and stochastic methods on the radiometric uncertainties limit their accuracy and application to cases where such data are available with an acceptable degree of confidence.A self-consistent and operational method is still required to estimate the uncertainties of IOPs without additional inputs or assumptions on the radiometric fluctuations.

Acknowledgment
The authors would like to thank NASA Ocean Biology Processing Group and individual data contributors for maintaining and updating the SeaBASS database.

Fig. 1 .
Fig. 1.Time series of ensemble-uncertainty of IOPs at 440 nm relative to the sum of derived IOPs.

) 237 Current
Advances in Uncertainty Estimation of Earth Observation Products of Water Quality www.intechopen.com

) 239 CurrentFig. 2 .
Fig. 2. Derived versus known errors of the IOPs estimated from the IOCCG data set for: (a) Chl-a absorption at 440 nm; (b) CDOM and detritus absorption at 440 nm; (c) SPM scattering at 550 nm; and (d) the total absorption at 440 nm.The data on the plots are log transformed.The coefficients of determination r 2 s and r 2 d are for stochastic and deterministic method respectively.

Fig. 3 .
Fig. 3. Derived versus known errors of the IOPs estimated from the NOMAD data set for: (a) Chl-a absorption at 440 nm; (b) CDOM and detritus absorption at 440 nm; (c) SPM scattering at 550 nm; and (d) the total absorption at 440 nm.The data on the plots are log transformed.The coefficients of determination r 2 s and r 2 d are for stochastic and deterministic method respectively.

243CurrentFig. 4 .
Fig. 4. Sum of variances versus the total variance of the IOCCG data set for: (a) Chl-a absorption at 440 nm; (b) absorption of detritus and CDOM at 440 nm; (c) SPM scattering at 550 nm; and (d) total absorption coefficient at 440 nm.
) and the imposed values of s and y.
. It should, nevertheless, be emphasized that the deterministic method is a well established technique to estimate retrieval errors.It can be used for the quantification of the combined accuracy of ocean color models and the parameterizations of IOPs, or model-parametrization setup.Its