Applying principles of metrology to historical Earth observations from satellites

Level 0 (L0) is the raw telemetry downlinked by a ground receiving station, comprising a mix of raw scientific observations from sensors together with engineering data from sensors and the platform. Transforming L0 data to scientifically useful products is a complex engineering task, and is a responsibility of the satellite operator.

The raw sensor data in L0 are in the form of binary numbers, referred to as 'counts' or 'digital numbers'. Each count is recorded using typically 10 to 16 bits. Such digitisation is coarse compared to laboratory metrology, and is driven by band-width limitations in downlinking data. Binary representation of the raw sensor values places a fundamental lower limit on the uncertainty present in calibrated radiances and derived geophysical estimates. 10-bit digitisation corresponds to ~0.1% resolution of the range.

Part of L0 processing involves estimating the satellite orbit and calculating, from the attitude of the platform in space, the origin of the measured radiances projected onto Earth's surface. This is known as geo-locating the satellite imagery. Orbit estimation can generally be improved retrospectively, so near-real time ('operational') satellite data may be superseded for non-time-critical applications by delayed mode and/or reprocessed data, using improved orbit information.

Level 1 (L1) datasets contain calibrated radiances or counts together with calibration parameters to map the counts into radiance. Auxiliary data locate the radiances in time, latitude and longitude on the Earth's surface. The satellite and solar azimuthal and zenith angles of the measured radiances are often included for convenience. Varying amounts of information from onboard calibration processes and engineering data (e.g. instrument temperatures, scanning rate, power consumption, etc) are provided. L1 data files are generally organised in image layers, one per channel, that have a common set of pixels. L1 data for microwave imagers can be an exception, since the footprint of an observation on the ground generally varies dramatically with microwave frequency for a fixed antenna size.

In scanning sensors, the individual measured radiances comprise roughly non-overlapping pixels that combine to create an image. For sensors on low Earth orbit, across-track scanning may be used, in which the image is obtained by successively scanning the Earth perpendicularly to the platform motion. In other cases, the scanning arrangement may be conical, in which cases pixels may, for convenience, be re-arranged (e.g. by nearest neighbour sampling) onto a rectilinear image grid. Pixels in an image may be referred to by their line (l) and element (e) indices which define the pixel position.

L1 radiances are typically calculated using an equation similar to:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{L}_{c,l,e}}={{a}_{0}}+{{a}_{1}}C_{c,l,e}^{{\rm E}}+{{a}_{2}}{{\left( C_{c,l,e}^{{\rm E}} \right)}^{2}}+0\nonumber \end{align} \tag{ 2.1 }$$

where: ${{L}_{c,l,e}}$ is the calculated Earth radiance for channel c of the pixel at image co-ordinate (l, e); $C_{c,l,e}^{{\rm E}}$ is the Earth count for the pixel recorded by the sensor; ${{a}_{0}},{{a}_{1}}\,{\rm and}\,{{a}_{2}}$ are calibration parameters, with ${{\boldsymbol{a}}^{{\rm E}}}=\left[ {{a}_{0}},{{a}_{1}},{{a}_{2}} \right]$. Uncertainty in the calculated radiance arises in part from uncertainties in the values of the quantities on the right hand side. In general, there are also other effects expected to have zero mean that contribute uncertainty in the calculated radiance; the '+  0' is included as a reminder of this, and will be discussed further in section 3.2.1. (Note that equation (2.1) is not representative of the calculation of radiance spectra from interferometers.)

In general, the calibration parameters are determined by in-flight calibration data plus other information; section 4 gives an example of the calibration parameters for a specific sensor. Sensors are generally reasonably linear, so ${{a}_{2}}$ is often small. ${{a}_{1}}$ is the 'gain' of the sensor. In-flight calibration systems typically have two reference points, such as a dark and bright target, to estimate changes in gain over time in flight. With two references, characterisation of non-linearity needs external information, such as pre-flight characterisation or in-flight characterisation against another sensor.

Changes in the gain and offset are to be expected as the sensor's space environment changes and the sensor degrades. An important source of change in the space environment is any precession of the orbit plane relative to the Sun (changing local solar time) over the mission lifetime. Material properties of components of the in-flight calibration system may degrade in time. The uncertainties associated with measured radiances will therefore evolve, and the stability of measurement may not be calculable from the satellite data alone.

Typically, 3 to 10 years after launch, the sensor will fail or the platform will be decommissioned. Multi-decadal datasets are built from sensors on a series of missions. Ideally, the sensors in a series would have identical spectral response and missions would overlap by a year or more. In practice, nominally equivalent channels have significant differences in SRFs and the overlaps between missions are not fully controllable.

Metrologists may be surprised that uncertainty estimates are not routinely included in L1 products under current practices in EO. A metrological approach to L1 products would include adopting the principle that every measured value in an L1 product should have associated context-specific uncertainty information (provided per datum if necessary). Although this principle (for Level 1 and higher products) and a set of associated guidance was endorsed by CEOS (Committee on Earth Observation Satellites) in 2010 as part of the Quality Assurance Framework for Earth Obervation (QA4EO) http://qa4eo.org/docs/QA4EO_guide.pdf this is still in the process of adoption at all space agencies.

Sensor channels are chosen to provide differential sensitivity to aspects of the state of the Earth's surface and atmosphere that are of interest. This differential sensitivity provides the information which is exploited in 'retrieval', which is inverse estimation from radiances of one or more geophysical variables. Retrieval algorithms used in EO are highly varied: analytic and implicit; empirical and physics-based; simple and complex; using only the observed radiances and employing significant external information.

The defining characteristic of L2 products is that retrieved geophysical variables are presented on the L1 image co-ordinates, i.e. at the highest spatial resolution possible.

The retrieval of L2 variables, $\boldsymbol{z}$, can often be written in a form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{z}=g\left( \boldsymbol{y},{{\boldsymbol{S}}_{\boldsymbol{\varepsilon }}},{{\boldsymbol{z}}_{\boldsymbol{a}}},{{\boldsymbol{S}}_{\boldsymbol{a}}},\boldsymbol{b} \right)+0.\nonumber \end{align} \tag{ 2.2 }$$

Here, z is the 'retrieved state' for a given pixel at line-element (l, e); this 'state vector' may contain a single variable (at the surface or a particular atmospheric level), a vertical profile of a single variable, or a set of various variables. Along with other lower-case bold variables, it is a column vector. $g$ is the measurement function for the L2 product. This may be a simple, analytic equation, or considerably more complex. The ${{n}_{c}}$ values of radiance used for the retrieval are in the 'observation vector' $\boldsymbol{y}={{\left[ {{L}_{1,l,e}},~\ldots ,{{L}_{{{n}_{c}},l,e}}~ \right]}^{{\rm T}}}$. ${{\boldsymbol{S}}_{\boldsymbol{\varepsilon }}}$ is an evaluation of the error covariance matrix of the measured radiances. ${{\boldsymbol{S}}_{\boldsymbol{\varepsilon }}}$ would ideally be based on uncertainty information available in the L1 product, but presently this is rarely, if ever, the case—an example of the need for improved metrology of Earth observations. ${{\boldsymbol{z}}_{\boldsymbol{a}}}$ is a prior estimate of the state with estimated error covariance ${{\boldsymbol{S}}_{\boldsymbol{a}}}$. ${{\boldsymbol{S}}_{\boldsymbol{\varepsilon }}}$, ${{\boldsymbol{z}}_{\boldsymbol{a}}}$ and ${{\boldsymbol{S}}_{\boldsymbol{a}}}$ are explicitly present in some retrieval algorithms such as optimal estimation (Rogers 2000), in which case they can be explicitly used to evaluate the error covariances of the retrieved state. In other cases, they are implicit and unevaluated. All retrieval algorithms have additional parameters used in retrieval, b. This can be as simple as a set of weights for combining radiances, or could be tens of thousands of spectroscopic parameters embedded within a radiative transfer model (or 'forward model') used in the inversion. More generically, one could regard all the terms ${{\boldsymbol{S}}_{\boldsymbol{\varepsilon }}},{{\boldsymbol{z}}_{\boldsymbol{a}}},{{\boldsymbol{S}}_{\boldsymbol{a}}},\boldsymbol{b}$ distinguished above as retrieval parameters, and write the retrieval equation as $\boldsymbol{z}=\boldsymbol{g}\left( \boldsymbol{y},\boldsymbol{b} \right)+0$. The '+  0' indicates that not all aspects of the inverse problem are necessarily captured by terms in the measurement function; some uncertainty is associated with those aspects of the inverse problem, even though their net effect is assumed to be zero mean. The symbols used here are somewhat adapted from widespread EO usage (e.g. Rogers (2000)) in that we use $z$ for the state vector rather than $x$.

The expression in equation (2.2) does not encompass all possibilities. The retrieval for pixel (l, e) may use radiances in other pixels in the image, to exploit spatial coherence in state. Classification or clustering of pixels may be performed prior to retrieval, for example to distinguish whether a water body or a particular class of land cover is present; different retrievals may then be applied to pixels according to their class or cluster.

Users of L2 products are presently accustomed to there being no or limited uncertainty information available in many L2 products, and perhaps relying on validation results reported in the literature. Users of L2 products often interrogate pixel-level quality indicators for indications about which values are more trustworthy than others. Quality indicators have different meanings in different L2 products. They may indicate the likely magnitude of uncertainty, or reflect the degree of confidence in the classification of the pixel, or a contextual judgement about the validity (in some sense) of the retrieval, or a combination of such factors. A metrological approach to L2 products would include adopting the principle that every geophysical estimate should have associated with it context-specific uncertainty information. EO experts working on several ECVs in the European Space Agency Climate Change Initiative (Hollmann et al 2013) have reached consensus with regards to best practice for geophysical products (Merchant et al 2017), and that consensus is consistent with this view.

The locations of geophysical observations in L2 do not exactly repeat between orbits, which is not always convenient. L2 product data volumes may also be large (many terabytes). Therefore L2 data are used where maximum spatial resolution is needed, but for applications tolerant of lower spatiotemporal sampling, L3 products are generated that are easier to use and have smaller volumes.

L3 products are gridded products, made by aggregating L2 values in space and/or time on a regular space-time grid. Aggregation involves the averaging of the L2 data obtained within the space-time volume of a given cell. The averaging may be weighted: for example, a weight of 0 may be applied to data whose L2 quality indicator is lower than that of the best quality data available within the cell. The formula used to obtain the gridded product value is usually a relatively simple equation, similar to:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle z \right\rangle =\sum\limits_{m}{{{w}_{m}}{{z}_{m}}+0}\nonumber \end{align} \tag{ 2.3 }$$

where the m index runs over all the L2 measured values within the latitude-longitude-time interval of the cell, w_m is the weight of the pixel value z_m, and $\left\langle z \right\rangle$ is the cell value in the L3 product, which represents the estimated mean value over the measurand across a spatio-temporal cell.

The uncertainty in the '+  0' term here includes sampling uncertainty (e.g. Bulgin et al (2016)), also known as representativity. Geophysical variability is continuous in space and time, and the L2 values sample this variability at certain, usually irregular, locations and/or times within the cell. Assuming the sampling is independent of the variability, the resulting sampling error has an expectation of zero, but can be drawn from an error distribution whose width is not negligible compared to the uncertainty in $\left\langle z \right\rangle$ arising directly from the uncertainties in the z_m.

A metrological approach to L3 products would involve the following principles: when creating aggregated products, propagate uncertainty appropriately, accounting for error correlation; clearly define the intended measurand represented by an aggregated observation and estimate the uncertainty from representativity/sampling accordingly. General solutions to these challenges are yet to be established and applied in practice.

Some applications require gridded data that are complete (gap free in space and time). An example would be where a climate or other environmental model is run with one or more variables being prescribed (i.e. imposed rather than calculated within the model). This may be done to explore the response of other variables to the 'real' history of the prescribed variable. A common example is running the atmospheric circulation of a global climate model while prescribing a spatially complete sea surface temperature field as a boundary condition.

L2 and L3 products are not gap free, and therefore L4 production requires additional processing involving spatio-temporal interpolation. L2 and L3 data from many satellite sensors may be used as inputs, and L4 analyses (as they are called) may also ingest values measured in situ. There are many possibilities, and the metrology of this level of processing will not be explored further in this paper.

The Guide to the Expression of Uncertainty in Measurement (GUM 2008; hereafter 'the GUM'), and its supplements provide guidance on how to express, determine, combine and propagate uncertainty. The GUM and its supplements are maintained by the JCGM (Joint Committee for Guides in Metrology), a joint committee of all the relevant standards organisations and the International Bureau of Weights and Measures, the BIPM.

The GUM defines the uncertainty of measurement as a:

parameter, associated with the result of a measurement, that characterizes the dispersion of the values that could reasonably be attributed to the measurand.

Uncertainty is non-negative. The 'standard uncertainty' is the measurement uncertainty expressed as a standard deviation. Evaluations of uncertainty in this paper refer throughout to standard uncertainty.

The error of measurement is defined in the GUM (section 2.2.19) as the:

result of a measurement minus a true value of the measurand

Since the true value of the measurand is unknown, the measurement error is unknown. However, by analysing the sources of errors in an EO measurement, the uncertainty in the measured value can be evaluated. This is the topic of section 3.2 below.

Estimating the uncertainty involves complex considerations, because any measurement is subject to several different sources of error, or 'effects'. The approach to uncertainty analysis presented below is to evaluate the uncertainty arising from different effects and combine these to determine the uncertainty associated with a measured value.

The measured value in each pixel of an EO image is the result of a sequence of steps and transformations. In transforming from raw data (L0) to calibrated radiances (L1), many measured values relevant to calibration measurements may be combined. In transforming L1 radiances to L2 geophysical products (figure 2), radiances from different spectral bands may be used. In L2 to L3  +  processing, data across different pixels are then combined. Where correlation exists between errors in different measured values (different wavelengths and/or pixels), this error correlation needs to be considered in the uncertainty analysis.

The GUM defines systematic and random errors, concepts that are widely used in science.

Random errors are errors manifesting independence: the error in one instance is in no way predictable from knowledge of the error in another instance. A complication arises in EO imagery when one instance of a parameter in the radiance measurement function is used in the calculation of the Earth radiance across many pixels. That component of the error in the radiance image is then correlated across pixels, even though the originating effect is random. Put another way, the originating random error contributes errors with a particular structure to the image.

Systematic errors are those that could in principle be corrected for if we had sufficient information to do so: that is, they arise from unknowns that could in principle be estimated rather than from chance processes. All systematic errors in EO are structured in that there is a pattern of influence on multiple data. They include, but are not limited to, effects that are constant for a significant proportion of a satellite mission—i.e. biases, for which the structure is a simple error in common.

When considering EO imagery, it is can be useful to categorise effects primarily according to their cross-pixel error correlation properties, as independent, structured or common effects.

3.1.2.1. Independent errors.

3.1.2.2. Structured errors.

3.1.2.3. Common errors.

3.1.2.4. Notes on usage of 'error' and 'uncertainty'.

Uncertainty analysis in the GUM begins with modelling the measurement, i.e. linking the measurand to the input quantities from which it is derived. A generic measurement model for a L1 radiance would be

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Y=f\left( {{X}_{1}},~{{X}_{2}},\ldots ,\boldsymbol{A} \right)+\Delta \nonumber \end{align} \tag{ 3.1 }$$

where: ${{X}_{1}},~{{X}_{2}},\ldots$ are input quantities; A is the vector of calibration parameters (which are also input quantities but are usefully distinguished); and $\Delta$ is an input quantity introduced to represent any inadequacy of the function f  to represent all phenomena that affect the measurand.

The equations, such as equation (3.1), used to populate L1 products, evaluate the measurand (radiance) using estimates of the input quantities. In the GUM, the convention is for estimates to be represented with the lower-case characters corresponding to the quantities written in upper case. Equation (2.1) is then seen as a particular case of the expression by which the measurand is estimated:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y=f\left( {{x}_{1}},{{x}_{2}},\ldots ,\boldsymbol{a} \right)+\delta \nonumber \end{align} \tag{ 3.2 }$$

where the input estimates include the recorded sensor counts, etc. This clarifies the meaning of the '  +  0' term previously introduced: 0 is our best estimate of $\delta$, which is the expectation of $\Delta$ (assuming we are using the best measurement model we can formulate).

The uncertainty in the measured value is derived from the evaluation of the uncertainty in each input estimate (or, strictly, from the distribution of possible values of ${{X}_{i}}$ given ${{x}_{i}}$). The uncertainty in all input estimates, including calibration parameters and $\delta$ is relevant. Evaluation of the uncertainty in $y$ means propagation through the measurement model of these uncertainties (or, strictly, distributions).

The GUM and its supplements describe both the 'Law of Propagation of Uncertainty' (hereafter, LPU) and Monte Carlo methods of uncertainty propagation.

The LPU propagates standard uncertainties for the input quantities through a locally-linear first-order Taylor series expansion of the measurement function to obtain the standard uncertainty associated with the estimate, $y$, of the measurand. Higher order approximations can be applied if necessary.

Monte Carlo methods approximate the input probability distributions by finite sets of random draws from those distributions and propagate the sets of input values through the measurement function to obtain a set of output values regarded as random draws from the probability distribution of the measurand. The output values are then analysed statistically, for example to obtain expectation values, standard deviations and error covariances. The measurement function in this case need not be linear nor written algebraically. Steps such as inverse retrievals and iterative processes can be addressed in this way. The input probability distributions can be as complex as needed, and can include distributions for digitised quantities, which are very common in EO, where signals are digitised for on-board recording and transmission to ground.

Monte Carlo methods can provide information about the shape of the output probability distribution for the measurand, deal better with highly non-linear measurement functions and with more complex probability distributions, and can be the only option for models that cannot be written algebraically. However, they are computationally more expensive, which is an important consideration with the very high data volumes of EO.

Often uncertainty analyses will use a combination of Monte Carlo methods and the LPU, for example by using Monte Carlo methods to determine the uncertainty for a particular quantity, which is then used as an input to LPU in a subsequent uncertainty analysis. In this section, we discuss and illustrate the LPU, while acknowledging the role of Monte Carlo methods.

The LPU for a single-variable measurand evaluated from ${{n}_{j}}$ input quantities with values, ${{\boldsymbol{x}}^{{\rm T}}}=~\left[ {{x}_{1}},\ldots {{x}_{{{n}_{j}}}} \right]$, can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{u}^{2}}\left( y \right)={{\boldsymbol{c}}^{{\rm T}}}\boldsymbol{S}(\boldsymbol{x})\boldsymbol{c}\nonumber \end{align} \tag{ 3.3. }$$

for a generic measurement function of the form $y=f\left( \boldsymbol{x} \right)$. Here, ${{\boldsymbol{c}}^{{\rm T}}}=[{{c}_{1}},\ldots ,{{c}_{{{n}_{j}}}}]$ contains the sensitivity coefficients of the measurand with respect to the input quantities. The input error covariance matrix, $\newcommand{\e}{{\rm e}} \boldsymbol{S}(\boldsymbol{x})=\left[ \begin{array}{@{}cccccccccccccccccccc@{}} {{u}^{2}}\left( {{x}_{1}} \right) & u({{x}_{1}},{{x}_{2}}) & \cdots \\ u({{x}_{2}},{{x}_{1}}) & {{u}^{2}}\left( {{x}_{2}} \right) & \cdots \\ \vdots & \vdots & \ddots \end{array} \right]$, is square. Each row and each column corresponds to one input quantity. The elements represent the error covariance associated with pairs of input quantities, and, down the diagonal, the error variance (uncertainty squared) of each input quantity.

The sensitivity coefficients are the first-order derivatives of the measurement function evaluated at the estimates of the input quantities and denoted by ${{c}_{j}}=\frac{\partial f}{\partial {{x}_{j}}}$. They can be analytically calculated when the measurement function can formally be differentiated, or can be numerically approximated using finite-difference methods, or can be experimentally determined.

Using the symbol conventions per satellite processing level in figure 2, consider a single-pixel L2 product calculated using a function $g({{Y}_{1}},\ldots {{Y}_{{{n}_{c}}}},\boldsymbol{B})$ which combines measured Earth radiances, ${{y}_{c}},$ in ${{n}_{c}}$ spectral channels and uses estimates, b, of the retrieval parameters. Then, the LPU takes the form ${{u}^{2}}\left( z \right)={{c}^{{\rm T}}}\boldsymbol{S}(\boldsymbol{y})\boldsymbol{c}$, where ${{\boldsymbol{y}}^{{\rm T}}}=\left[ {{y}_{1}},\ldots {{y}_{{{n}_{c}}}},{{\boldsymbol{b}}^{{\rm T}}} \right]$. The vector of sensitivity coefficients relates the Earth radiance values in different spectral channels and the retrieval parameters to the retrieved L2 variable (see section 2.2.3):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\boldsymbol{c}}^{{\rm T}}}=\left[ \frac{\partial g}{\partial {{y}_{1}}},~\ldots ,~\frac{\partial g}{\partial {{b}_{1}}},\ldots \right]\nonumber \end{align} \tag{ 3.4 }$$

and ${{\left[ \boldsymbol{S}(\boldsymbol{y}) \right]}_{1,1}}={{u}^{2}}\left( {{y}_{1}} \right)$, ${{\left[ \boldsymbol{S}(\boldsymbol{y}) \right]}_{1,2}}=u\left( {{y}_{1}},{{y}_{2}} \right),~$ etc. The derivatives in equation (3.2) are evaluated at the estimates (measured values) of the corresponding input quantities.

The LPU can be extended to multiple measurands (output quantities) (GUM-102 2011). Consider the case of retrieval of L2 variables from radiances, $\boldsymbol{Z}=\boldsymbol{g}\left( \boldsymbol{Y} \right)$, where Z consists of several geophysical quantities of interest obtained jointly from a common set of multi-channel radiances, Y. The propagated uncertainties associated with values of z are the contents of the following error covariance matrix:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{S}=\boldsymbol{CUR}{{\boldsymbol{U}}^{{\rm T}}}{{\boldsymbol{C}}^{{\rm T}}}.\nonumber \end{align} \tag{ 3.5 }$$

Here, $\boldsymbol{UR}{{\boldsymbol{U}}^{{\rm T}}}$ is the error covariance matrix for the input quantities, expressed in terms of a diagonal matrix of input standard uncertainties, U, and a matrix of error correlation coefficients between input quantities, $\newcommand{\e}{{\rm e}} \boldsymbol{R}=\left[ \begin{array}{@{}cccccccccccccccccccc@{}} 1 & {{r}_{1,2}} & \cdots \\ {{r}_{2,1}} & 1 & \cdots \\ \vdots & \vdots & \ddots \end{array} \right]$. The sensitivity coefficient matrix has elements ${{\left[ \boldsymbol{C} \right]}_{n,m}}=\frac{\partial {{g}_{n}}}{\partial {{y}_{m}}}$, i.e. $\boldsymbol{C}=\frac{\partial \boldsymbol{g}}{\partial \boldsymbol{y}}$, evaluated at the estimates of the corresponding input quantities. The output error covariance matrix, S, contains the error covariance of the estimates of the measurands. Even if none of the input terms has correlated errors (and thus R is diagonal), S will generally still be non-diagonal, if (as is likely) the sensitivity matrix includes terms for those output quantities that depend on the same input quantities. A summation of the impact of all error effects is implicit in the summation over all the input quantities.

A subtly different restatement of the LPU will prove useful in section 3.3, where the contribution to the overall error covariance matrix will be calculated one-effect-at-a-time and then summed over the effects explicitly. This approach will be adopted for considering the error covariances between different pixels in an image and between different channels at a given pixel location. This situation differs from the multiple-measurand case above in that now each pixel/channel radiance (output value) is the same function of its own set of input values—i.e. ${{\boldsymbol{y}}^{{\rm T}}}=\left[\,f\left( {{\boldsymbol{x}}_{1}} \right),f\left( {{\boldsymbol{x}}_{2}} \right),\ldots \right]$. In this case, ${{\left[ \boldsymbol{C} \right]}_{n,m}}=0$ for $n\ne m$, and C is square and diagonal. For a single effect, $k$, applying to measurement function term, $j$, the contributed error covariance is ${{\boldsymbol{C}}_{j}}{{\boldsymbol{U}}_{k}}{{\boldsymbol{R}}_{k}}\boldsymbol{U}_{k}^{{\rm T}}\boldsymbol{C}_{j}^{{\rm T}}$. This looks the same as equation (3.5), but the meanings and structures of the matrices are a particular case. These points will be explained more fully below.

The uncertainty analysis we describe here is centred on the measurement function for calibrated radiances, the data of L1 products. Typically, some input quantities are directly determined (by measurement, parameterisation or data analysis), while others are determined through their own measurement function of other input quantities. Thus, there can be a hierarchical aspect to uncertainty analysis. To help organise a measurement-function centred analysis of uncertainty, it can be useful to prepare a graphical representation in the form of an uncertainty analysis diagram (figure 3).

Zoom In Zoom Out Reset image size

We find that most radiance measurement function terms are sensitive to one or more effects, and the uncertainty due to each effect needs to be estimated. Quantities not subject to effects include mathematical and physical constants, and agreed reference quantities. In some cases, contributing effects may be estimated separately and then combined, and in other cases, all the effects operative on a single quantity may only be jointly estimated. Many effects are possible in EO: examples at L1 include detector noise, temperature gradients across reference targets, stray light ingress and temperature sensitivity of electronics. The '+  0' term in a radiance measurement function represents effects relating to the assumptions underlying the form of the measurement function, e.g. that the sensor response is linear or quadratic in underlying radiance. Example effects for a particular sensor are discussed in section 4. Uncertainty analysis starts at the terminations of the uncertainty analysis diagram and propagates uncertainty through each measurement function, in sequence through the hierarchy as necessary.

Uncertainty analysis presumes that the result of a measurement has been corrected for all recognized systematic effects, and that every effort has been made to identify such effects (GUM 2008 section 3.2.4). This means that the measurand will be as accurate as possible given the current state of knowledge. When we perform analyses at the effects level we need to decide whether the effect could be fully, or partially corrected for, and if so we should apply the correction. The residual uncertainty after correction is smaller. Performing uncertainty analysis therefore often has the co-benefit of improving EO data.

Each effect causes uncertainty, meaning it gives rise to an unknown error in each measured value of radiance. Errors may be independent, in common, or have intermediate degrees of correlation. Three different forms of correlation need to be understood to characterise L1 uncertainty for users:

• The error correlation between different quantities in the measurement function. Such error correlations arise from effects in common between different quantities in the measurement function, and where quantities are determined from a common fit to data.
• The error correlation between image pixels for a given spectral band, which may arise when different image pixels have used some of the same input data in the calculation.
• The error correlation between spectral bands for a given pixel, which may arise when different image pixels have used some of the same input data in the calculation.

It is advantageous for uncertainty analysis to formulate the measurement function to avoid the first form of correlation, so that effects are defined by construction to be independent. The LPU can then be expressed as a sum of covariance matrices over effects. Effect independence can be achieved by using hierarchical analysis to pre-combine terms that have a common error component. For example, consider calculation of a reflectance value that depends on the satellite view angle. The satellite view angle may be calculated from longitude and latitude values, both of which are sensitive to a common time error. We can construct the measurement function in terms of the satellite view angle and consider the effect of time error on that angle. A second example arises when there is error covariance between empirical calibration parameters that were obtained in a joint fitting procedure. We can treat these parameters as a single, vector input quantity, and handle the error correlation using their estimated error covariance.

Error correlation between pixels and spectral bands also arises from errors in common. It is important to distinguish 'common effects' from 'common errors'. The same effect might be present in multiple measurements (for example, all measured radiances will be sensitive to noise), but a common error only arises if the same instance (value) of a quantity is used. Earth-view radiances in an image come from different measurements for each channel and pixel, and therefore a different error instance in each radiance. In contrast, a single view of an on-board calibration target may be used to calibrate several lines of pixels, and therefore there is a common error instance for all those pixels. Similarly, any error in characterising the temperature of a black-body calibration target affects all the calibrated thermal channels.

Since our uncertainty analysis is based on metrologically-independent effects by construction of the measurement function, we can find the error covariance matrix for a vector of Earth radiance values (e.g. radiance values for different spectral channels and/or for different image pixels) as a sum over effect-specific error covariance matrices:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{S}=\sum\limits_{j}{\sum\limits_{k|j}{{{\boldsymbol{C}}_{j}}{{\boldsymbol{U}}_{k}}{{\boldsymbol{R}}_{k}}\boldsymbol{U}_{k}^{{\rm T}}\boldsymbol{C}_{j}^{{\rm T}}}.}\nonumber \end{align} \tag{ 3.6 }$$

This is written as a sum over input quantities, $j$, in the measurement model. The sum is over the associated effects, $k|j$, i.e. $k$ given $j$. We write it this way to make explicit that there may be multiple effects for a given input quantity, although, since effects influence only one quantity by construction, we could simply sum over effects. The sensitivity coefficient matrix, ${{\boldsymbol{C}}_{j}}$, is diagonal, because each radiance is separately determined from its own evaluation of its measurement model and any quantities in common are addressed in the correlation matrix (see below). Each row/column refers to one corresponding radiance in the vector measurand (i.e. per spectral band or per image pixel). The units of the sensitivity coefficient are the units of radiance divided by the units of the input quantity, $j$, and the derivatives are evaluated at the estimate of the quantity.

The uncertainty matrix ${{\boldsymbol{U}}_{k}}$ is a diagonal matrix of the same dimension, with the standard uncertainty associated with the effect down the diagonal. These uncertainties are in the units of the input quantity, j .

The correlation coefficient matrix ${{\boldsymbol{R}}_{k}}$ shows the error correlation for the effect $k$ between the evaluations of the measurement function for the different radiances. Typically, all the radiances have the same form of measurement function, with the same terms. Some of these terms have an independent value for each evaluation of radiance, as when each radiance is estimated from the corresponding Earth counts. In this case, there is no correlation, and ${{\boldsymbol{R}}_{k}}$ is the identity matrix. Other terms may have the same value for each evaluation of radiance, in which case ${{\boldsymbol{R}}_{k}}$ is the matrix of ones. In more complex cases, terms may have different values whose errors are neither fully correlated nor uncorrelated, in which case ${{\boldsymbol{R}}_{k}}$ will have a complex, non-diagonal form.

We choose to decompose the error covariance matrix for the term into uncertainty and correlation matrices, ${{\boldsymbol{U}}_{k}}{{\boldsymbol{R}}_{k}}\boldsymbol{U}_{k}^{{\rm T}}$, because correlation structures are often fixed even if uncertainty varies. Not all effects lend themselves to this treatment, and in such cases the error covariance matrix, ${{\boldsymbol{S}}_{k}}$, may be better evaluated directly. In this case, the contribution to the error covariance matrix for that effect would be ${{\boldsymbol{C}}_{j}}{{\boldsymbol{S}}_{k}}\boldsymbol{C}_{j}^{{\rm T}}$.

The description above has considered a general set of radiances constituting a vector measurand. Typically, we are interested in one of three specific cases: radiance values across the elements of a line for a particular channel; radiance values across the lines for a particular element in a particular channel; and radiance values for a given pixel (line and element) across the channels measured by the sensor. In the first case, the result is the cross-element error covariance matrix for the given channel, $\boldsymbol{S}_{e}^{l}$. In the second case, it is the cross-line error covariance matrix, $\boldsymbol{S}_{l}^{e}$. The third case gives the cross-channel error covariance matrix, $\boldsymbol{S}_{c}^{\,p}$. The indices in this notation operate as follows. The subscript indicates the dimension across which the error covariances are calculated: thus the cross-channel matrix has subscript, $c$, for example. The superscript indicates a dimension that is not averaged across. Thus, the cross-element error covariance is evaluated for a particular line, $l$. Since it would never make sense to average error covariance matrices across channels, there is no need to have $c$ in the superscript: that is implicit. To calculate a cross-element error covariance matrix representative of many lines in an image, ${{\boldsymbol{S}}_{e}}$, we would need to evaluate $\boldsymbol{S}_{e}^{l}$ for many different lines, and then average these. Denoting by ${{\left\langle x \right\rangle }_{l}}$ the average of $x$ over an adequate set of lines, ${{\boldsymbol{S}}_{e}}={{\left\langle \boldsymbol{S}_{e}^{l} \right\rangle }_{l}}$. Further, $\boldsymbol{S}_{c}^{\,p}$ indicates that a cross-channel error covariance matrix is calculated at a particular pixel, $p=\left( l,e \right)$. Such a notation is required because there are many possible covariance matrices that can be calculated and used, and they need to be uniquely indexed.

To make this description more concrete we can consider the cross-channel error covariance matrix calculation for a three-channel sensor subject to two effects. Let the first effect, $k=1$, be fully correlated, because a common value, ${{x}_{1}}$, is used in the calculation of all three channel radiances. ${{x}_{1}}$ might quantify the state of a common calibration target, for example. Then, evaluating at a pixel, $p=P$, we have:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}lllccccccccccccccccc@{}} & \boldsymbol{S}_{c}^{\,p=P,k=1}=\boldsymbol{C}_{c}^{\,p=P,k=1}\boldsymbol{U}_{c}^{\,p=P,k=1} \nonumber \\ &\boldsymbol{R}_{c}^{\,p=P,k=1}\boldsymbol{U}{{_{c}^{\,p=P,k=1}}^{T}}\boldsymbol{C}{{_{c}^{\,p=P,k=1}}^{T}} \nonumber \\ &= \left( \begin{array}{cccccccccccccccccccc} \frac{\partial {{y}_{1}}}{\partial {{x}_{1}}} & 0 & 0 \nonumber \\ 0 & \frac{\partial {{y}_{2}}}{\partial {{x}_{1}}} & 0 \nonumber \\ 0 & 0 & \frac{\partial {{y}_{3}}}{\partial {{x}_{1}}} \nonumber \end{array} \right)\left( \begin{array}{cccccccccccccccccccc} {{u}_{{{x}_{1}}}} & 0 & 0 \nonumber \\ 0 & {{u}_{{{x}_{1}}}} & 0 \nonumber \\ 0 & 0 & {{u}_{{{x}_{1}}}} \nonumber \end{array} \right)\left( \begin{array}{cccccccccccccccccccc} 1 & 1 & 1 \nonumber \\ 1 & 1 & 1 \nonumber \\ 1 & 1 & 1 \nonumber \end{array} \right) \nonumber \\ &{{\left( \begin{array}{cccccccccccccccccccc} {{u}_{{{x}_{1}}}} & 0 & 0 \nonumber \\ 0 & {{u}_{{{x}_{1}}}} & 0 \nonumber \\ 0 & 0 & {{u}_{{{x}_{1}}}} \nonumber \end{array} \right)}^{{\rm T}}}{{\left( \begin{array}{cccccccccccccccccccc} \frac{\partial {{y}_{1}}}{\partial {{x}_{1}}} & 0 & 0 \nonumber \\ 0 & \frac{\partial {{y}_{2}}}{\partial {{x}_{1}}} & 0 \nonumber \\ 0 & 0 & \frac{\partial {{y}_{3}}}{\partial {{x}_{1}}} \nonumber \end{array} \right)}^{{\rm T}}}. \nonumber \\ \end{array}\nonumber \end{align} \tag{ 3.7 }$$

This fully correlated effect ($\boldsymbol{R}_{c}^{\,p=P,k=1}=\boldsymbol{J}$, the all-ones matrix) has the same uncertainty for each channel, but differing sensitivity coefficients. An effect, $k=2$, that has no correlation between channels, would have a correlation matrix $\boldsymbol{R}_{c}^{\,p=P,k=2}=\boldsymbol{I}$. In that case the uncertainty may be different for each channel: the values of the equivalent term in the measurement function for each channel may be separately quantified, with a different uncertainty in each case. A cross-channel error covariance that is representative for many pixels in an image could be calculated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\boldsymbol{S}}_{c}}={{\left\langle \sum\limits_{k}{\boldsymbol{S}_{c}^{\,p,k}} \right\rangle }_{p}}.\nonumber \end{align} \tag{ 3.8 }$$

Analogous steps would be necessary to find the error covariance matrix for Earth radiance values obtained in different spatial pixels, such as across the elements of a line. The matrices in that case would have a dimension equal to the number of pixels of interest. The correlation coefficient matrix for a single effect, $\boldsymbol{R}_{p}^{k}$, would represent the correlation coefficients between the errors in a particular term across that set of pixels.

We previously introduced (section 3.3.1) the uncertainty analysis diagram, in order graphically to illustrate (i) that uncertainty analysis may have a hierarchical aspect, and (ii) the different paths by which effects may influence the measurand. The measurement function is centrally placed in the uncertainty analysis diagram, and the analysis of EO we present here is likewise 'measurement-function centred'. The uncertainty analysis diagram is an aid to the organisation and documentation of metrologically traceable measurement-function centred uncertainty analysis of complex EO sensors. For each term of the measurement function, the origin of uncertainties is identified. On the branches of the diagram, the sensitivity coefficients are placed that relate the uncertainty in terms of the uncertainties in the measured values. The end points of the diagram identify the complete set of effects analysed.

There are choices about how the uncertainty propagation is in practice to be calculated, and these choices are reflected in the structuring of the uncertainty analysis diagram. Careful thought is needed to balance the need for an understandable measurement function while avoiding too many hierarchical levels of measurement function and constructing each effect as metrologically independent (section 3.3.1). In documentation and in the processing code, it is necessary to record the origin of each term in the measurement function and the data used in processing.

Each effect requires several aspects to be quantified and documented. A useful aid to organising this effort is to synthesise the knowledge of each effect in its own 'effects table'. The effects table records the information needed about this effect to perform a full uncertainty analysis and propagation calculations. An effects table suitable for EO metrology is shown as table 2.

Table 2. Table for codifying the correlation structure in L1 radiances.

Table descriptor		How this is codified	Notes
Name of effect		A name for each source of error	Names should be unique
Measurement function term affected		Name and standard symbol of affected term	Usually an effect will only affect a single term, by construction
Instruments in the series affected		Identifier of the specific sensors/satellite platform for which this effect matters	Some effects may arise only for specific sensors in a series
Correlation type and form	Cross-element	One of a number of error correlation forms arising in EO	Codifies the form of spatio-temporal error correlation, addressing all the relevant dimensions that are relevant.
	Cross-line
	Between images/orbits
	Over time
Correlation scale	Cross-element	In units of pixels, lines, time or images/orbits—what is the scale of the correlation form?	The scale parameter has a meaning specific to each correlation form
	Cross-line
	Between images/orbits
	Over time
Channels/bands	List of channels/bands affected	Channel names in standard form	Some effects may only affect certain channels
	Error correlation coefficient matrix	A matrix (or specification to calculate the matrix)	This is $\boldsymbol{R}_{c}^{k}$
Uncertainty	PDF shape	Functional form of estimated error distribution for the term	Often Gaussian, or rectangular
	Units	Units in which PDF shape is expressed	Units of term, or can be as percentage etc
	Magnitude	Dispersion of PDF as standard uncertainty	May point to data quantifying uncertainty per pixel/channel
Sensitivity coefficient		Value(s) or equation for of sensitivity of measurand to term	Yields measurand units when multiplied with term standard uncertainty

First, the table documents uniquely the name of the source of error and which term is affected. Even within a series of nominally similar sensors, a particular effect may not arise in every case, so the affected sensors must be identified.

Next, the functional dependence ('form') of spatio-temporal error correlation is codified. The form of error correlation must be characterised across several 'dimensions': between elements, between lines, between images/orbits and over time. If correlation does not arise over one of these dimensions, the corresponding correlation form is recorded as independent for the given effect. Across a dimension where correlation does arise, the error correlation structure can take various forms, depending on the origin of the common component of error. For example, there will be an uncertainty associated with the measurement of an on-board calibration target. On-board calibration may be performed at intervals, meaning that the error in calibration result affects all lines observed between calibration cycles. The error correlation then has a rectangular form, being 1 for the lines that use that calibration result, and 0 for lines prior and after that refer to the result of an independent calibration cycle. If the calibration results are smoothed using a boxcar running average over several calibration cycles, which introduces a triangular form of error correlation.

Several forms of error correlation arise often in EO, and in the effects table we use defined names (e.g. rectangular, triangular, bell-shaped) to refer to these. A correlation scale parameter is defined for each of these forms that specifies the range of elements, lines, images or time over which the error correlation is non-zero.

The effects table then characterises the error across channels, specifying the cross-channel correlation matrix, which is often either I or J, as discussed in section 3.3.3.

The estimated shape and dispersion of the error distribution are tabulated next. The uncertainty may be as simple as a constant per channel for a given effect, or may be specified by an equation or parameterisation, or may be a layer of data per pixel per channel. The same may be true of the sensitivity coefficients.

Various methods are available for estimating uncertainties and sensitivity coefficients.

Uncertainty can be evaluated from the standard deviation or the Allan deviation (Allan 1966) of repeated measured values, where available. Uncertainty may also be evaluated using simulation, the results of pre-flight tests, or characterised behaviour of similar instruments. Particularly for historical satellite sensors, we may not have all the information ideally needed, so pragmatic estimates must be made.

Where there is an analytic measurement function, a sensitivity coefficient is most easily determined as its partial derivative. Where it is not possible to evaluate a partial derivative, sensitivity coefficients can be determined through perturbation of the values of the inputs of a software function, or from experimental (e.g. pre-flight) data, for example that showing the sensitivity of the instrument gain to temperature.

EO missions are typically described by a set of documentation. To undertake a metrological analysis of a historical EO sensor, the available documentation should be used to understand the sensor operation, including hardware and procedures for in-flight calibration, as illustrated for AVHRR below.

In the IR, the AVHRR uses two different types of solid state detector. An InSb detector is used for its 3.7 µm channel, which has a close-to-linear response apart from minor circuit non-linearity. The HgCdTe detectors used for the 11 µm and 12 µm channels have a distinct non-linear response. The AVHRR on-board calibration system uses an internal calibration target which has a quoted emissivity of 0.985 140. Note that this emissivity value is purely theoretical and is considered to be significantly more uncertain than implied by the number of significant figures quoted. Its value is, however, known to be significantly less than 1, since reflected components (proportional to 1 minus emissivity) are seen in the calibration signal, including solar contamination (e.g. Cao et al (2004)).

The AVHRR is a scanning radiometer that collects images as a sequence of scan lines at right angles to the direction of travel of the satellite over the ground. The in-flight calibration procedure makes ten measurements (as counts) per scan line when viewing the internal calibration target (ICT) and ten measurements per scan line of a space view. Four platinum resistance thermometers (PRTs) measure the temperature of the ICT and allow an estimate of the spectral radiance of the ICT to be made, using Planck's Law and the quoted emissivity. From the spectral radiance, the channel-integrated spectral radiance (usually referred to as simply 'radiance') of the ICT is calculated by integrating spectral radiance across the (assumed known) spectral response function (SRF) of a given channel. This ICT radiance coupled with the counts for the ICT and space-view enable the instrument gain to be estimated. The gain is combined with other calibration parameters to convert the counts observed when looking at the Earth into Earth-view radiance. Radiance is often then expressed as brightness temperature, by inverting the channel-integrated Planck function numerically.

The calibration algorithms of the AVHRR have evolved during the series of missions from 1978 to present. Originally, no account was made for the non-linearity of the 11 µm and 12 µm HgCdTe detectors. Subsequent calibration schemes used a range of look-up tables or correction terms for this non-linearity. The current operational calibration is based on Walton et al (1998) and involves calculation of a linear radiance which is then corrected using a quadratic correction. This formulation is intrinsically problematic (Mittaz et al 2009). Further, using parameters derived from pre-launch data within this formulation gives rise to significant biases (Wang and Cao 2008, Mittaz and Harris 2011).

An improved radiance measurement function, modified from Mittaz and Harris (2011), is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & {{L}_{{\rm E}}}={{a}_{0}}+\frac{{{a}_{1}}{{L}_{{\rm T}}}-{{a}_{2}}\hat{C}_{{\rm T}}^{2}}{{{{\hat{C}}}_{T}}}{{{\tilde{C}}}_{{\rm E}}}+{{a}_{2}}\tilde{C}_{{\rm E}}^{2}+0 \nonumber \\ & {{{\hat{C}}}_{{\rm T}}}=\frac{1}{n}\sum{\left( {{C}_{{\rm S},i}}-{{C}_{{\rm T},i}} \right)} \nonumber \\ & {{{\tilde{C}}}_{{\rm E}}}={{C}_{{\rm S}}}-{{C}_{{\rm E}}} \nonumber \\ \end{array}\nonumber \end{align} \tag{ 4.1 }$$

where ${{L}_{{\rm E}}}$ is the Earth radiance, ${{L}_{{\rm T}}}$ is the estimated radiance of the ICT, ${{\hat{C}}_{{\rm T}}}$ is the counts observed when looking at the ICT averaged over a number of scanlines (in practice a space-view count C_s minus a target count) and ${{C}_{{\rm E}}}$ is the Earth-view count (in practice a space-view count minus the Earth-view count) and ${{a}_{1}},{{a}_{2}},{{a}_{3}}$ are calibration parameters that have physical interpretations. We assume that this quadratic form represents the non-linearity adequately though this will be discussed in section 4.2.5. For the AVHRR the dark, space-view count, is set to be higher than the signal count for operational reasons.

An uncertainty analysis diagram for the AVHRR measurement function is shown in figure 4. It shows the sources of uncertainty from their origin (shown on the outside portions of the diagram) through to the uncertainty in the final measurand. The sources of uncertainty range from independent effects with a physical origin, such as detector noise sources, to structured effects associated with the estimate of the internal calibration target (ICT) radiance to common effects related to the harmonisation process. A full discussion of all effects is beyond the scope of this paper, and this complete diagram is included mainly to illustrate the scale of a full uncertainty analysis for a relatively simple instrument. At the stage of constructing an uncertainty analysis diagram, it is also best practice to include all potential sources of uncertainty, even those suspected to be negligible. For a metrologically-traceable uncertainty estimate, the neglect of any contributing term should be justified and documented.

Zoom In Zoom Out Reset image size

Note also that the diagram includes some intermediate processes that themselves contribute to the total uncertainty budget. For example, the ICT radiance is derived from measurements of the ICT temperature from four PRTs. The uncertainty in the ICT radiance is a combination of the uncertainty in the ICT temperature (affected by PRT noise, PRT bias, digitisation and the representativity of this measurement to the ICT temperature seen by the detector, and the issue of thermal gradients across the ICT), uncertainty in the spectral response function, and approximations from use of look-up tables in conversions ('discretisation'). The ICT radiance is also a source of cross-channel error correlations, because the same ICT temperature is used to calibrate all three IR channels and therefore any error will be common to all channels.

In the upper half of the diagram are the components related to harmonisation. Harmonisation involves the re-estimation of the calibration parameters for consistency across the satellite series by using sensor-to-sensor collocated observations ('matches'). Many of the sources of uncertainty listed here are repeated from the lower portion of the diagram, because the harmonisation fit depends on a set of measured counts, each subject to the same effects. The matches have their own uncertainties associated with the quality of the collocation. A detailed analysis of uncertainties in sensor-to-sensor matches has been undertaken by Hewison (2013).

The diagram also includes the sensitivity coefficients on each branch that enable the source uncertainties to be propagated to the measurand.

The uncertainty analysis diagram documents in a visual form the full set of error effects to be analysed and their connections to the measurement function. For each effect represented on the diagram, we use an effects table (section 3.4.2) to aid analysis and documentation.

In the following sections we use specific examples to illustrate the use of the effect tables and consider effects including noise in the space-view counts (C_s), uncertainties in the estimation of the calibration target temperature and the uncertainty due to imperfect characterisation of the HgCdTe detectors (which is an uncertainty associated with the term '+  0' of the radiance measurement function). The 'counts' example illustrates structured errors, the ICT temperature example illustrates channel-to-channel errors and the 'non-quadratic' example illustrates the metrological principle that effects need to be considered even when there is no direct way of quantifying them.

Uncertainty in the space-view counts can be derived from the telemetry stored in the AVHRR Level 1b file. Error correlation arises because, as is commonly done in EO, the space-view counts are averaged over a number of scanlines to reduce noise in the calibration values. Table 3 is the corresponding effects table.

Table 3. Effects table for noise in space view counts.

Table descriptor		How this is codified
Name of effect		Space view count uncertainty
Affected term in measurement function		Space view count (${{{C}}_{{\rm S}}}$)
Channels/bands		Channel 3B, 4 and 5 (3.7 µm, 11 µm and 12 µm)
Correlation type and form	Within scanline (pixels)	Rectangular absolute
	From scanline to scanline (scanlines)	Triangular relative
	Between images/orbits (orbits)	None
	Between channels/bands	None
Correlation scale	Within scanline (pixels)	(−∞, +∞)
	From scanline to scanline (scanlines)	(−N, +N) where N is AVHRR dependent
	Between images/orbits (orbits)	None
	Between channels/bands	None
Uncertainty PDF shape		Gaussian
Uncertainty units		Counts
Uncertainty magnitude		Estimated from Allan deviation from space views accumulated over a complete orbit
Sensitivity coefficient		$\frac{\partial L_E}{\partial C_S} = \frac{a_1 L_T}{\hat{C}_T} \left(1-\frac{\tilde{C}_E}{\hat{C}_T}\right)+a_2(\tilde{C}_E-\hat{C}_T)$

The effects table shows that there are two dimensions of error correlation that have to be dealt with appropriately. First, since the calibration parameters are recalculated each scanline, the space-view counts error is the same for all pixels across a given scanline, and therefore has a correlation of 1 across the whole scanline, which is a 'rectangular absolute' systematic correlation. Note that while the space-view counts error is in common, the sensitivity coefficient differs for the different pixels across the scanline, since the Earth-view count is present in the analytic expression for the derivative. The corresponding radiance errors across the scanline therefore have a predictable structure (determined by the sensitivity coefficient variation) scaled by an unknown random error (since the space counts error is random). Second, because multiple (N) scanlines are averaged to create the calibration coefficients, there is a correlated error in the running-mean space-view counts between scanlines used within the average. The running mean uses a constant weight, and so the form of the correlation is triangular (the correlation decreases for scanlines increasingly far apart) with a base width of  ±N scanlines (beyond which the correlation is zero). N varies between different AVHRR sensors and can vary from approximately 20 scanlines up to ~100 scanlines.

The measurement function for the AVHRR contains the radiance of the ICT used to calibrate the instantaneous instrument gain. This quantity is itself calculated from the ICT temperature, and that temperature is itself subject to different sources of uncertainty. The ICT radiance is therefore one of the intermediate effects. One of the sources of uncertainty for the ICT radiance is having to estimate the ICT radiant temperature based on temperatures measured at four points only using platinum resistance thermometers attached to the ICT baseplate. The uncertainty then arises because the ICT is subject to large (>1 K), time-variable thermal gradients, so the relationship between the true ICT radiant temperature and these four point temperature measurements is not simple. During the satellite lifetime, the radiant temperature was estimated operationally by the arithmetic mean of the four temperatures. This approach is simple and would give reasonable estimates for an ICT in a well-regulated environment. However, due to the thermal gradients across its surface as well as stray light issues during times of rapid variations in the instrument's thermal state the use of a simple average will lead to significant errors. Figure 5 shows two different sensors, one where this effect is thought to be large (the AVHRR on NOAA-12) and one where it is likely to be small (NOAA-17). Previous studies of this effect on the earlier AVHRRs (e.g. Trishchenko et al (2002)) have commented on the range of variability seen in the PRT measurements (see figure 5) and that these could be used as some kind of broad estimate of the uncertainty on the ICT temperature. This proposed approach cannot be considered as metrologically robust, since it involves no attempt to remove the major source of error.

Zoom In Zoom Out Reset image size

In a metrological approach, we remove known sources of systematic error (such as errors in the ICT temperature estimate) to the best of our knowledge, and quantify the remaining uncertainty. To do this, we have to consider the process giving rise to the effect (variations in the thermal gradients across the instrument) and how the error will propagate into measured radiances.

Four PRT temperatures are inadequate to determine much detail about the complex thermal gradients across the surface of the ICT. We can, however, exploit understanding of the behaviour of the 3.7 µm channel to constrain the problem. This channel uses an InSb detector which is known to be linear for the radiance levels observed from Earth observation. The 3.7 µm channel gain is therefore expected to be constant around the orbit. We can then reasonably infer that significant deviations from a constant gain reflect errors in the ICT temperature estimate. A direct mapping from gain error to ICT temperature error constrained by the PRT measurements themselves can then be made (full details of the procedure will be published elsewhere). Figure 6 shows data for NOAA-12 from the same orbit as shown in figure 5 (left panel). The left hand plot of figure 6 shows the original and corrected 3.7 µm channel gain. The right-hand plot shows the ICT temperature estimates made using the simple mean of the PRTs (original) and corrected by analysing the 3.7 µm channel gain. The corrected gain has lower variance, suggesting a substantial reduction in the ICT temperature error and the remaining variability in the 3.7 µm channel gain can be used to estimate the ICT temperature uncertainty. The improved ICT temperature can then also be used to calibrate the 11 µm and 12 µm channels, where constant gain is not expected. This metrological approach estimates ICT temperature with reduced systematic errors, and provides a method to evaluate the remaining uncertainty, for propagation to uncertainty in measured radiances.

Zoom In Zoom Out Reset image size

This effect is also one that will introduce channel-to-channel error correlation because any remaining error in the ICT temperature will be present in the calibration of all the infrared channels. This effect has never previously been characterised for the AVHRR. For any given pixel we can use equation (3.7) to estimate a cross-channel error covariance matrix. An example is shown below (table 4) together with the related correlation matrix for a single pixel of NOAA-12. The cross-channel error correlation between the 11 and 12 µm channels is large, and this is significant for retrieval of geophysical quantities using these channels. The effect table for the ICT thermal gradient effect is shown as table 5. We note that currently for the AVHRR one of the table entries (correlation type and form between images/orbits) is labelled as 'Unknown' because we do not currently know what this is apart from a general statement that we might expect close in time orbits to have similar correlation structures and those far apart in time having less correlation. This will be looked at in the future.

Table 4. Correlation and Covariance matrices for a single pixel of NOAA-12 data showing the correlation and covariances introduced by the ICT temperature error effect on the 3.7 µm, 11 µm and 12 µm channels. The unit of the covariance matrix is K².

IR channels correlation matrix	IR channels covariance matrix
$\newcommand{\e}{{\rm e}} \left( \begin{array}{@{}cccccccccccccccccccc@{}} 1.0 & 0.0209 & 0.0209 \\ 0.0209 & 1.0 & 0.9958 \\ 0.0209 & 0.9958 & 1.0 \end{array} \right)$	$\newcommand{\e}{{\rm e}} \left( \begin{array}{@{}cccccccccccccccccccc@{}} 0.1074 & 0.0017 & 0.0016 \\ 0.0017 & 0.0726 & 0.0689 \\ 0.0016 & 0.0689 & 0.0714 \end{array} \right)$

Table 5. Effect table for the error in the ICT temperature estimate.

Table descriptor		How this is codified
Name of effect		ICT Temperature Error
Affected term in measurement function		${{L}_{{\rm T}}}$
Channels/bands		Channel 3, 4 and 5 (3.7 µm, 11 µm and 12 µm)
Correlation type and form	Within scanline (pixels)	Rectangular absolute systematic
	From scanline to scanline (scanlines)	Triangular relative systematic
	Between images/orbits (orbits)	Unknown at the moment
	Between channels/bands	Across the IR channels
Correlation scale	Within scanline (pixels)	(−∞, +  ∞)
	From scanline to scanline (scanlines)	(−N, +N) where N is AVHRR sensor specific
	Between images/orbits (orbits)	None
	Between channels/bands	None
Uncertainty PDF shape		Gaussian
Uncertainty units		Kelvin
Uncertainty magnitude		Estimated from corrected 3.7 µm channel gain variance
Sensitivity coefficient		$\frac{\partial {{L}_{{\rm E}}}}{\partial {{T}_{{\rm ICT}}}}=\frac{{{a}_{1}}}{{{{\hat{C}}}_{{\rm T}}}}\frac{\partial {{L}_{{\rm T}}}}{\partial {{T}_{{\rm ICT}}}}$

The effects tables shown as table 6 represent effects related to the assumption that the AVHRR HgCdTe detector can be modelled as a quadratic function of counts with a constant non-linear coefficient. This assumption is built into the AVHRR measurement equation, and therefore any error in this assumption is an effect on the  +0 term of the measurement equation. The origin of the non-linearity is dominated by Auger recombination, which is itself related to the lifetime of semiconductor carriers. The exact details of the non-linear behaviour is related to doping levels, lattice defects, etc, in the detector and is not expected to be strictly quadratic. We cannot determine this effect in detail from the pre-launch measurements and do not have access to the original manufacturer's data, and so we cannot directly characterise the errors introduced by the quadratic assumption. So, we must make a Type-B uncertainty estimate, based on expert judgement and independent knowledge. In this case, the latter consists of a numerical model of HgCdTe detectors of the geostationary operational environmental satellite (GOES) Imager (Bicknell 2000). The model determines the Auger recombination lifetimes of the carriers and derives variations in the predicted voltage seen for a given input photon flux. We tuned the model parameters to match the photon fluxes and magnitude of non-linearity seen in the AVHRR sensors. The top two plots of figure 7 shows that the modelled errors in the brightness temperature from approximating the detector behaviour with a quadratic form are of order millikelvin, which is small.

Table 6. Effect table for the assumptions regarding detector non-linearity.

Table descriptor		How this is codified
Name of effect		Non-quadratic non-linearity
Affected term in measurement function		+0
Channels/bands		Channel 4 and 5 (11 µm and 12 µm)
Correlation type and form	within scanline (pixels)	Rectangular absolute systematic
	from scanline to scanline (scanlines)	Rectangular absolute systematic
	between images/orbits (orbits)	Rectangular absolute systematic
	between channels/bands	Rectangular absolute systematic
Correlation scale	within scanline (pixels)	(−∞, +∞)
	from scanline to scanline (scanlines)	(−∞, +∞)
	between images/orbits (orbits)	(−∞, +∞)
	between channels/bands	None
Uncertainty PDF shape		Gaussian
Uncertainty units		Radiance
Uncertainty magnitude		Type-B estimate from detector modelling
Sensitivity coefficient		$\frac{\partial {{R}_{E}}}{\partial +0}=1$

Zoom In Zoom Out Reset image size

While the size of the uncertainties added from the constant quadratic non-linear model used in the measurement equation have been shown to be small compared to other effects, it is important to understand that without the modelling studies undertaken here we would not have a good estimate as to the impact of these particular effect and therefore would not be able to claim metrological traceability in the sense used in this paper. If we had either ignored these effects or not estimated the size of uncertainty in a justifiable manner we would have an incomplete uncertainty model.

We have simulated AVHRR BTs for one day's worth of full-swath orbital data covering a scene temperature range from 200 K to 300 K and generated at each simulation point 1000 realisations of the simulated 'measured' BT to compare with the 'true' value inputted to the simulation.

Figure 8 shows the simulated error distributions for a cold (200 K) and a hot (300 K) scene for the 11 µm and 12 µm channels. The error distributions are not Gaussian: digitisation causes each distribution to comprise a few separated peaks. The degree of separation depends on the scene temperature because a single count corresponds to a smaller BT increment at higher radiances. This effect is random between pixels, since the detector noise that is digitised is random between pixels. The spread of each peak reflects the effect of the calibration uncertainty, including noise in the calibration measurements.

Zoom In Zoom Out Reset image size