Radiance Uncertainty Characterisation to Facilitate Climate Data Record Creation

Abstract

The uncertainty in a climate data records (CDRs) derived from Earth observations in part derives from the propagated uncertainty in the radiance record (the fundamental climate data record, FCDR) from which the geophysical estimates in the CDR are derived. A common barrier to providing uncertainty-quantified CDRs is the inaccessibility to CDR creators of appropriate radiance uncertainty information in the FCDR. Here, we propose radiance uncertainty information designed directly to facilitate estimation of propagated uncertainty in derived CDRs at full resolution and in gridded products. Errors in Earth observations are typically highly structured and complex, and the uncertainty information we propose is of intermediate complexity, sufficient to capture the main variability in propagated uncertainty in a CDR, while avoiding unfeasible complexity or data volume. The uncertainty and error correlation characteristics of uncertainty are quantified for three classes of error with different propagation properties: independent, structured and common radiance errors. The meaning, mathematical derivations, practical evaluation and example applications of this set of uncertainty information are presented.

1. Introduction

Climate data records (CDRs) are generated from Earth observations within several programmes [1,2,3,4] in support of societal and scientific needs to quantify change and variability in Earth’s climate system. The uncertainty in the trends and variations of climate variables in CDRs is arguably as important to understand as the trends and variations themselves, in order to avoid incorrect inferences or optimistic levels of confidence in conclusions drawn from the data [5,6,7,8]. Recently, greater efforts have been directed to CDR generation [9,10], and specifically to quantifying uncertainty in CDRs [11,12,13]. Consideration has been given to good practice, drawing on experience across a range of essential climate variables [14], a major conclusion being that in cases where uncertainty in the CDR is highly variable, uncertainty evaluations should be available per datum. This implies the need to understand and model the dominant processes causing uncertainty in the CDR.

Climate data records from Earth observation are derived by a series of transformations of data. Raw telemetry is processed by a ground segment to ‘level-1’ datasets, comprising radiances or counts with sufficient information to convert to calibrated radiances. Note that radiances may then be expressed as brightness temperature or reflectance in level-1 products, and in this paper, we use “radiance” generically to encompass all such quantities. It also needs to be acknowledged that some level-1 data comprise signals of a different nature, as in the cases of altimetry and gravimetry. In level-1 data, observations are also located in time and space. A level-1 dataset of sufficient length, consistency and stability to support the derivation of a CDR may be referred to as a fundamental climate data record (FCDR). From the level-1 dataset, often using retrieval techniques involving multi-channel inferences, estimates of geophysical quantities are derived, and the full resolution dataset of the retrieved variable is described as a CDR at ‘level 2’. For some applications, regularly gridded data are preferred, often at a lower spatio-temporal resolution, and level-2 data may be averaged on such a regular grid to generate a CDR at ‘level-3’. Further processing to other forms of higher-level dataset may also be done.

At each transformation from level 1 to level 3 and beyond, there is both (i) propagation of uncertainty from the lower-level data to the higher, and (ii) introduction of new sources of uncertainty by virtue of the parameters, assumptions and/or limitations of the method of data transformation [15]. The relative importance of these two sources of uncertainty in a CDR will vary, depending on the nature of the radiance data and of the retrieval process. Generally, rigorously quantifying CDR uncertainty requires propagation of radiance uncertainty through the transformations to higher levels.

In practice, a large proportion of level-1/FCDR datasets contain either no or very limited radiance uncertainty information, meaning that CDR creators have no practical means to evaluate propagated radiance uncertainty. In part, this situation arises because it is not obvious what selection of radiance uncertainty information is needed to adequately support CDR creation, given that providing comprehensive uncertainty information would entail infeasible data volumes and/or complexity (see below).

The focus of this paper is to address the question, “What radiance uncertainty information should be provided to support the characterisation of uncertainty in derived CDRs?” The proposed approach is intended to support propagation of uncertainty to levels 2 and 3, which means that inter-channel and spatio-temporal error correlations must be included. The key contributions of the paper are (i) to propose a set of summary information on errors and their correlations that it is feasible and sufficient to provide, and (ii) to summarise the mathematics and computing techniques required to calculate this uncertainty information. References [15,16,17,18,19] illustrate the working out of these concepts in particular cases.

Although facilitating CDR generation is our focus, we note that radiance uncertainty information at level 1 is also directly relevant to environmental re-analysis using data assimilation e.g., [20]. In data assimilation, the ideal approach is to assimilate level-1 observations. An evaluation of the observation-to-model error covariance associated with FCDRs is required by the assimilation, and the observation-related component of this is directly addressed by the methods described in this paper. Data assimilation may also help validate uncertainty information in FCDRs and CDRs [21], although this topic is beyond the present scope.

The remainder of the paper proceeds as follows. In Section 2, we describe the steps a CDR creator would want to be able to take to propagate and quantify radiance uncertainty in their CDR. From these steps, we infer the uncertainty information that should be provided at level 1 in an FCDR. In Section 3, we outline the understanding necessary to be able to characterise the radiance uncertainty in a single datum, based on understanding sources of error in the instrument. We emphasise that CDR creators should not need such detailed knowledge of the instrument: providing quantification and understanding of level-1 uncertainty should be the responsibility of space agencies, in our view. In Section 4, we present the theory and formulae for turning knowledge of per-datum instrument uncertainty into summary uncertainty information suitable for inclusion in FCDRs, ready for CDR creators to use. Practical evaluation of the formulae requires care, and so these considerations are presented in Section 5. Section 6 briefly illustrates what results can look like in a real case without extensive discussion of the sensor or error sources addressed; this section focuses on the meaning and interpretation of the uncertainty information shown. The concluding Section 7 includes summational remarks regarding the assumptions and implications of the ideas presented, as well as some closing discussion.

2. Requirements for Propagation of Uncertainty to a Climate Data Record

Typically, estimation of a geophysical variable for inclusion in a CDR (i.e., at level-2/3) involves several pre-processing steps followed by a multi-channel retrieval (inference of the variable’s value from a combination of radiance values at different wavelengths). The pre-processing steps may include pixel classification (including target selection or cloud screening) and selection of the valid channel set for the retrieval (which may vary between day and night scenes, for example). In the retrieval step that is applied to the selected pixels, a set of $n_c$ observations $\mathit{y}={\left[\begin{array}{cc}\begin{array}{cc}{y}_1& \cdots \end{array}& y_{n_c}\end{array}\right]}^{\mathrm{T}}$ is exploited to infer a geophysical variable, $z$. In creating a level-2 CDR, this process of inference usually operates on single pixels. Let the retrieval be

$$z=g\left(\mathit{y},\mathbf{\beta }\right)+0$$

(1)

where $g$ is the retrieval function which may or may not be expressible analytically, and $\mathbf{\beta }$ is the set of retrieval parameters used (which may include externally sourced information). The term $+0$ obviously does not change the retrieved value $z$, but is explicitly included here to remind us that there are three components to the uncertainty in the retrieved value:

• The propagation of radiance uncertainty through $g$ to $z$.
• The propagation of parameter uncertainty through $g$ to $z$.
• The uncertainty, “$u\left(0\right)$”, introduced into $z$ by any approximations and inherent ambiguities of the retrieval that are not directly traceable from uncertainty in $\mathit{y}$ or $\mathbf{\beta }$.

This paper is concerned with methods to make the first component listed above straightforward for creators of CDRs to calculate, by defining the uncertainty characterisation that they need.

The propagated radiance standard uncertainty can readily be calculated for a given set of radiances $\mathit{y}$ if the radiance error covariance matrix, ${\mathit{S}}_{ϵ}$, is known:

$$u\left(z\right)=\sqrt{{\mathit{c}}^{\mathrm{T}}{\mathit{S}}_{ϵ}\mathit{c}},$$

(2)

where $\mathit{c}={\left[\begin{array}{cc}\begin{array}{cc}\frac{\partial g}{\partial y_1}& \cdots \end{array}& \frac{\partial g}{\partial y_{n_c}}\end{array}\right]}^{\mathrm{T}}$. Equation (2) is the law of propagation of uncertainty [22], which expresses (i) a first-order approximation in which a small error in, say, $y_1$ propagates linearly to an error in $z$ in proportion to the derivative $\frac{\partial g}{\partial y_1}$, and (ii) that independent errors combine as the sum of squares. Readers unfamiliar with the matrix form of (2) can readily verify that in the case where ${\mathit{S}}_{ϵ}$ is diagonal it corresponds to a textbook expression that may be more familiar: $u\left(z\right)=\sqrt{{\displaystyle \sum }_i{\left(\frac{\partial g}{\partial y_i}u\left(y_i\right)\right)}^2}$. The matrix form additionally allows for the convenient propagation of correlated errors through any off-diagonal terms in ${\mathit{S}}_{ϵ}$. An off-diagonal term takes the value $r_{c{c}^{\prime }}u\left(y_c\right)u\left(y_{c^{\prime }}\right)$, where $r_{c{c}^{\prime }}$ is the correlation coefficient of errors arising between radiance in channels $c$ and $c^{\prime }$. When errors are independent, $r_{c{c}^{\prime }}=0$ for all $c\ne c^{\prime }$, and the error covariance matrix is diagonal.

The answer to the question of this paper, “What radiance uncertainty information should be provided to support the characterisation of uncertainty in derived CDRs?”, thus might seem quite simple: “Provide the radiance error covariance matrix”.

If we assume that it is understood how to estimate ${\mathit{S}}_{ϵ}$ for a given set of radiances, there remains a practical problem with this answer: ${\mathit{S}}_{ϵ}$ is not a constant matrix. The radiance uncertainty typically arises from several error-causing effects, the relative importance of which may change with instrument state and with the scene characteristics. The degree of cross-channel error correlation is also likely to differ between effects. Thus, the uncertainties, $u\left(y_c\right)$, and the correlation structures, $r_{c{c}^{\prime }}$, usually vary. On the other hand, it is not practical and affordable to provide ${\mathit{S}}_{ϵ}$ per pixel from the point of view of data volumes [14]. Thus, there is a need to find a useful and efficient approximate representation, preferably one that reflects most of the variability in ${\mathit{S}}_{ϵ}$. That need is addressed in this paper.

Turning to level-3 CDRs, a gridded CDR can generally be considered to contain weighted combinations, $\langle z\rangle$, of $n_z$ selected retrieved values,$\mathit{z}$, from the full resolution CDR:

$$\langle z\rangle =h\left(\mathit{z},\mathbf{\gamma }\right)+0.$$

(3)

The process of selecting which pixels contribute to $\langle z\rangle$ includes identifying which pixels lie within the required space-time cell for a given level-3 datum, and perhaps a process to select the best pixels to include. The parameters, $\mathbf{\gamma }$, may include weights, length-scales, correlation coefficients, etc. As at level 2, there are three categories of uncertainty in $\langle z\rangle$:

• The propagation of retrieved-value uncertainty through $h$ to $\langle z\rangle$.
• The propagation of gridding parameter uncertainty through $h$ to $\langle z\rangle$.
• The uncertainty, $u\left(0\right)$, introduced by incomplete sampling across the space-time domain which the level-3 datum represents [23]. The uncertainty from incomplete sampling is sometimes referred to as “standard error in the mean”, and reflects observations errors, geophysical variability within the space-time cell, and correlations in these.

Depending on the nature and assumptions of the gridding, the gridding parameter uncertainty, the second bullet above may be justifiably neglected. For example, if the aim is to average evenly across the cell area and the fields of view of the contributing pixels are well quantified, then the weight of each pixel in the cell average may be determined with negligible uncertainty.

Consider the propagation of the retrieved-value uncertainty at pixel level to the cell-mean value. The propagated standard uncertainty is:

$$u(\langle z\rangle )=\sqrt{{\mathit{c}}^{\mathrm{T}}{\mathit{S}}_{\mathit{z}}\mathit{c}},$$

(4)

where the sensitivity coefficients are now $\mathit{c}={\left[\begin{array}{cc}\begin{array}{cc}\frac{\partial h}{\partial z_1}& \cdots \end{array}& \frac{\partial h}{\partial z_{n_z}}\end{array}\right]}^{\mathrm{T}}$, and ${\mathit{S}}_{\mathit{z}}$ is the error covariance matrix of the selected pixels being averaged within the cell. Continuing to assume that the retrieved values are derived one-pixel-at-a-time, what is the structure of ${\mathit{S}}_{\mathit{z}}$? The simplest possibility arises if the errors in the retrieved pixels are independent, in which case ${\mathit{S}}_{\mathit{z}}$ is diagonal and the uncertainty in the cell average diminishes as the familiar $1/\sqrt{n}$. This is rarely, if ever, the case. Retrieval parameters are likely to be strongly related for nearby pixels, and generally so too will be their errors. Noise may be the largest uncertainty in radiances, and although noise may often be independent between pixels, the noise frequency spectrum can be such that there is inter-pixel error correlation [24]. On-board calibration cycles for imagers are once-per-scan-line or less frequent, and calibration errors are correlated across scans and across the duration of the calibration cycle. Error correlations mean that the uncertainty $u(\langle z\rangle )$ will be greater than the uncertainty expected on the basis of “$1/\sqrt{n}$”. Making a robust estimate of the level-3 uncertainty requires knowledge of error covariance between retrieved values across the cell (level 2), which in turn is in part determined by spatial error covariance in level-1 radiances—i.e., by cross-pixel correlations in radiance errors.

In summary, for propagation of radiance uncertainty to the CDR at level 2 and level 3, estimates of cross-channel and cross-pixel error covariance are needed, and should be available as part of the FCDR. The next two sections briefly describe how such summary error covariance information can be developed from an understanding of the sensor.

3. Methods 1—Radiance Uncertainty Analysis

There are various sources of information from which a good understanding of radiance uncertainty can be developed. Satellite sensors are extensively analysed during design, build and pre-flight characterisation [24]. On-board calibration systems provide information that enables quantification of sensor noise and related calibration uncertainty [25,26], as well as characterisation of the noise spectrum and cross-pixel error correlation [18]. Specifications and pre-flight characterisation of sensor features such as calibration of target properties [26], detector non-linearity [27], sensitivity to instrument temperatures [28] or quantisation [29,30] can be used to infer estimates of uncertainty in the calculated radiance from a range of sources of errors. These are general comments, and it is beyond the scope here to give a comprehensive account of the sources of information that underpin an uncertainty analysis of radiance for a given sensor or sensors in general. A lengthier account of principles of uncertainty analysis applied to Earth observation radiances, with examples, is presented in [15].

The focus of this paper is to show how the outcome of such an uncertainty analysis can be condensed into summary information about radiance uncertainty that is then useful in CDR development. For the sake of argument, therefore, we assume that an uncertainty analysis has been undertaken for a given sensor, resulting in quantification of an adequately comprehensive set of $n_k$ sources of error (“effects”), which will be indexed as $k=1\dots n_k$. Let the radiance in a given channel, $c$, for a certain pixel, $p$, be calculated as

$$y=f\left(x_1,\dots ,x_{n_j}\right)+0$$

(5)

where $x_j$ is one of the $n_j$ input quantities that appear in the calculation of the FCDR radiance, such as sensor counts, calibration target counts, calibration parameters and measurements of instrument state. The function $f$ is the measurement function for radiance. Quantification of an effect, $k$, means that we:

• have identified which input quantity (index $j$) in the radiance measurement function the effect influences;
• can evaluate the corresponding standard uncertainty $u_k\left(x_j;c,p\right)$ associated with that input quantity, which may be a constant for all $p$, may change for every $p$, or be in-between these extremes;
• can estimate the correlation in error in $x_j$ between any pair of channels $\left(c,c^{\prime }\right)$; this correlation will usually lie between 1, where the same value of $x_j$ is used in the measurement equation for both channels, and 0, where the values are independent;
• can estimate the error correlation in error in $x_j$ between any pair of pixels $\left(p,p^{\prime }\right)$ of the image in a given channel; and
• can quantify $\partial f/\partial x_j$ (the sensitivity coefficient of the measurand to the input quantity), either analytically or by perturbation.

Once this set of information has been developed across the key error-causing effects, the required radiance uncertainty information can be calculated, as derived in the next section of this paper.

Here, we now define how the information listed above, i.e., the result of radiance uncertainty analysis, is encapsulated quantitatively, in various matrices.

First, we define the identity matrix, $\mathit{I}$, and the matrix of ones, $\mathit{J}$. The former describes the error correlation between quantities that are independent. The latter describes the error correlation between quantities subject to errors that are in common (or in fixed proportions). We have:

$$\mathit{I}=\left[\begin{array}{ccc}1& 0& \cdots \\ 0& 1& \cdots \\ ⋮& ⋮& \ddots \end{array}\right];\mathit{J}=\left[\begin{array}{ccc}1& 1& \cdots \\ 1& 1& \cdots \\ ⋮& ⋮& \ddots \end{array}\right]$$

(6)

Second, we introduce notation for error correlation matrices. There are several such matrices that we will need to distinguish carefully, and thus the notation is necessarily detailed. All error correlation matrices will be represented as ${\mathit{R}}_{i_1,\left\{k\right\}}^{i_2,i_3\dots }$. The subscripts and superscripts have meanings according to their position, as follows.

• $i_1$ is an index across a dimension of the satellite image, and indicates the dimension across which the matrix expresses the error correlations. This subscript can take the following values: $c$, the channel dimension, indicating a cross-channel error correlation matrix; $e$, the image element dimension, indicating a cross-element error correlation matrix; $l$, the image line dimension; and $p$, used to index pixels generically (across both lines and elements), as when indexing a selection of pixels within a particular L3 grid-cell.
• The superscripts $i_2,i_3$, etc., identify indices (line, elements, input quantities or effects) over which the matrix does not average or integrate: $p=P$ in the superscript position indicates the error covariance matrix is evaluated for a specific pixel, $P$; $p$ alone means we are referring to evaluation for single pixels generically, without specifying which pixel particularly; $k=K$ in this position indicates that the matrix refers to one specific effect, $K$; and so on.
• $\left\{k\right\}$, if present, indicates a sub-set of effects whose combined effect is quantified by the error correlation matrix; absence of this optional label indicates all effects have been combined; the possible sub-sets of effects will be defined in Section 4.

The partition between superscripts and subscripts is, in summary, as follows: the subscripts indicate the aspects the covariance matrix does include, while the superscripts express aspects across which the matrix does not average or integrate.

For example, ${\mathit{R}}_c^{p,k=K}$ refers to a cross-channel error correlation matrix (subscript $c$) for one identified effect, $K$ (superscript $k=K$), evaluated at a single, arbitrary image pixel, $p$. We use the same labelling convention for corresponding error covariance matrices, such that ${\mathit{S}}_p^{p,j}$ would designate a cross-pixel error covariance matrix (subscript $p$) from all effects affecting the $j$^th measurement equation input quantity, evaluated for a single, arbitrary image-pixel.

Third, we introduce notation for uncertainty. The uncertainty of the measured value in pixel $p=\left(l,e\right)$ in channel $c$ arising from the combination of all effects $k$ is:

$$u\left(y;c,l,e\right)=\sqrt{{\displaystyle \sum }_j{\displaystyle \sum }_{k|j}{c}_j{\left(x_j\left(c,l,e\right)\right)}^2{u}_k{\left(x_j\left(c,l,e\right)\right)}^2}$$

(7)

where the notation of the summation makes explicit that each input quantity can be affected by more than one effect, hence the double sum over all $k|j$ (effect $k$ given input quantity $j$) and $j$ (input quantities). The sensitivity coefficients, $c_j=\partial f/\partial x_j$, change with the input quantity but not the effect, whereas each effect contributes its own uncertainty, $u_k\left(x_j\right)$, in the relevant input quantity. All these effects are likely to vary between channels, and may vary across lines and/or elements in the image. Using the same labelling convention as error correlation matrices, the matrix ${\mathit{U}}_{i_1}^{i_2\dots }$ is a diagonal matrix with the diagonal elements calculated across channels, lines, elements or pixels (as indicated by $i_1$) for the specific instances indicated by the superscript(s), $i_2\dots$. For example,

$${\mathit{U}}_c^{k,l,e}=\left[\begin{array}{ccc}{u}_k\left(x_j\left(c=1,l,e\right)\right)& 0& \cdots \\ 0& u_k\left(x_j\left(c=2,l,e\right)\right)& \cdots \\ ⋮& ⋮& \ddots \end{array}\right]$$

(8)

shows an uncertainty matrix with diagonal elements being the uncertainty in the $j$^th input quantity arising from one of the effects, $k$, relevant to that input quantity.

Lastly, we define the sensitivity matrix, in a similar way to uncertainty. Thus, for example,

$${\mathit{C}}_c^{j,l,e}=\left[\begin{array}{ccc}{c}_j\left(x_j\left(c=1,l,e\right)\right)& 0& \cdots \\ 0& c_j\left(x_j\left(c=2,l,e\right)\right)& \cdots \\ ⋮& ⋮& \ddots \end{array}\right]$$

(9)

shows a sensitivity matrix of the sensitivities to the $j$^th input quantity of the radiance measurement equation across different channels, for a single choice of image line and element.

4. Methods 2—Equations to Define Summary Uncertainty Information

A conclusion of Section 2 was that propagation of radiance uncertainty to the CDR at level 2 and level 3 requires estimates of the cross-channel and cross-pixel error covariances to be provided at level 1. In this section, the equations necessary to quantify these covariances are developed, starting with calculation of the cross-channel error covariance matrix.

Considering a single effect, $k$, at an arbitrary pixel, $p$, in an image, the cross-channel radiance error covariance matrix is equal to:

$${\mathit{S}}_c^{p,k}={\mathit{C}}_c^{p,j}{\mathit{U}}_c^{p,k}{\mathit{R}}_c^{p,k}{\mathit{U}}_c^{p,k}{\mathit{C}}_c^{p,j}$$

(10)

where $j$ is the index of the input quantity in which the effect causes error. The form of this is analogous to the earlier equations for propagation of uncertainty, with the error covariance on the right-hand side (representing the uncertainty that is to be propagated to the radiance) has been decomposed into a diagonal matrix of uncertainty values, ${\mathit{U}}_c^{p,k}$, and the error correlation matrix, ${\mathit{R}}_c^{p,k}$. The total error covariance matrix from all the analysed effects is the sum of the individual error covariances:

$${\mathit{S}}_c^p={\displaystyle \sum }_j{\displaystyle \sum }_{k|j}{\mathit{C}}_c^{p,j}{\mathit{U}}_c^{p,k}{\mathit{R}}_c^{p,k}{\mathit{U}}_c^{p,k}{\mathit{C}}_c^{p,j}$$

(11)

Given the results of radiance uncertainty analysis, this channel covariance matrix could be provided for every pixel, in principle. However, while the radiance data scales as $n_c$, the matrix scales as $n_c^2$, which is impractical for users, since FCDRs are already often multi-terabyte or petabyte datasets. An average estimate can be made and provided, by averaging over an adequate sample of pixels:

$${\mathit{S}}_c=\langle {\mathit{S}}_c^p\rangle {}_p$$

(12)

where ${\langle \rangle }_p$ is an operator to average over a sample of pixels. (In line with the nomenclature introduced earlier, the absence of a superscript in ${\mathit{S}}_c$ indicates that all effects have been combined and that there is averaging over pixels.)

Provision of ${\mathit{S}}_c$ would be an advance over present practice, but can be further improved to allow CDR producers the ability to reflect some of the variability in uncertainty across an image. (Where data are provided as whole orbits, for example, the instrument state may change and radiance uncertainty may vary significantly.) One approach is to provide, instead of ${\mathit{S}}_c$, the average radiance error correlation matrix per image, ${\mathit{R}}_c$, which can be combined with a data layer quantifying per pixel uncertainty. The approximation assumed in this approach is that the error correlations are relatively constant, whereas uncertainty significantly changes. Any covariance matrix, being positive definite and symmetric, can be represented in the form $\mathit{S}=\mathit{U}\mathit{R}{\mathit{U}}^{\mathrm{T}}$. To see this, let element $\left(m,n\right)$ of the matrix, $\mathit{S}$, be ${[\mathit{S}]}_{m,n}^{}$, etc. Make the following identifications:

$${\left[\mathit{U}\right]}_{m,n}={\delta }_{m,n}\sqrt{{\left[\mathit{S}\right]}_{m,n}};{\left[\mathit{R}\right]}_{m,n}=\frac{{\left[\mathit{S}\right]}_{m,n}}{\sqrt{{\left[\mathit{S}\right]}_{m,m}{\left[\mathit{S}\right]}_{n,n}}}$$

(13)

where ${\delta }_{m,n}$ is zero when $m\ne n$ and one when $m=n$. Calculating $\mathit{U}\mathit{R}{\mathit{U}}^{\mathrm{T}}$ recovers $\mathit{S}$. This is the principle on which an overall error correlation matrix can be calculated for a combination of effects. Thus, if ${\mathit{U}}_c$ is the diagonal matrix whose diagonal elements are the square roots of the diagonal elements of ${\mathit{S}}_c$, the required average error correlation matrix is:

$${\mathit{R}}_c={\mathit{U}}_c^{-1}{\mathit{S}}_c{\mathit{U}}_c^{-1}$$

(14)

The use of ${\mathit{R}}_c$, combined with uncertainty estimates across the image, will be described later.

The same equations apply to calculating error covariance across elements and lines. For example, the representative cross-element error correlation matrix (for a given channel) is defined by:

$${\mathit{R}}_e={\mathit{U}}_e^{-1}{\mathit{S}}_e{\mathit{U}}_e^{-1};{\mathit{S}}_e={\langle {\displaystyle \sum }_j{\displaystyle \sum }_{k|j}{\mathit{C}}_e^{l,j}{\mathit{U}}_e^{l,k}{\mathit{R}}_e^{l,k}{\mathit{U}}_e^{l,k}{\mathit{C}}_e^{l,j}\rangle }_l$$

(15)

and equivalent expressions for cross-line error correlations can be obtained by swapping the $e$ and $l$ indices. The dimensions of ${\mathit{R}}_l$ and ${\mathit{R}}_e$ are, respectively, the squares of the number of along-track lines and the number of elements squared ($n_l^2\mathrm{and}{n}_e^2$); this is often likely to represent an impractical additional data volume. If further data reduction is required, the approximation may be made that error correlations vary spatially as a function of element or line separation only, i.e., as a function of $\left|e-e^{\prime }\right|$ and $\left|l-l^{\prime }\right|$, respectively. If we write, for example, ${\mathit{R}}_e=\left[\begin{array}{cccccc}1& r_{1,2}& r_{1,3}& r_{1,4}& r_{1,5}& \ddots \\ r_{1,2}& 1& r_{2,3}& r_{2,4}& r_{2,5}& \ddots \\ r_{1,3}& r_{2,3}& 1& r_{3,4}& r_{3,5}& \ddots \\ r_{1,4}& r_{2,4}& r_{3,4}& 1& r_{4,5}& \ddots \\ r_{1,5}& r_{2,5}& r_{3,5}& r_{4,5}& 1& \ddots \\ \ddots & \ddots & \ddots & \ddots & \ddots & \ddots \end{array}\right]$, then the average error correlation value for pixels on a given line separated by a difference in element index equal to ${\Delta }_e$ can be calculated using:

$${\left[{\mathit{r}}_e\right]}_{{\Delta }_e}={\displaystyle \sum }_{e=1}^{n_e-{\Delta }_e}{r}_{e,e+{\Delta }_e}/(n_e-{\Delta }_e)$$

(16)

with a similar equation for the cross-line average error correlation. ${\mathit{r}}_e$ is a vector containing $n_e$ elements evaluating this equation for pixel separations from ${\Delta }_e=0$ (for which ${\left[{\mathit{r}}_e\right]}_{{\Delta }_e}=1$) to ${\Delta }_e=n_e-1$. From this, an ‘average’ ${\mathit{R}}_e$ can be constructed by substituting ${\left[{\mathit{r}}_e\right]}_{\left|e-e^{\prime }\right|}$ for $r_{e,e^{\prime }}$.

To summarise the implications of the above mathematical development, methods have been proposed to calculate:

• a representative cross-channel error correlation matrix (affordable volume-wise)
• a representative cross-line error correlation function (from which an approximate error correlation matrix etc can be inferred where the matrix itself is unaffordable volume-wise), and
• a representative cross-element error correlation function

from which all the necessary error covariance matrices for uncertainty propagation to the CDR can be generated using the uncertainty provided for each pixel.

The development above is framed in terms of a single set of uncertainty and error correlations being able to represent all effects. This is a strong approximation and does not enable the correlation structure to respond to relative changes in uncertainty from effects with contrasting error structures. (Such relative changes in uncertainty can arise, for example, when an instrument’s gain changes and the balance of pure noise versus calibration uncertainty may change.)

To lessen the degree of the approximation, we can subdivide effects into families of similar error correlation structure, providing the set of information above for each family. Such a subdivision, of course, increases complexity and data volume, so a balance needs to be struck. A traditional split is to partition “random” and “systematic” error sources. We have found that this is not adequate to describe the errors in a satellite image, because there are important intermediate effects that cause errors that are in common or partially correlated across certain groups of pixels (and thus are neither random nor systematic in the usual sense). The minimum complexity that reflects this situation is a tripartite division of effects. A tripartite approach to uncertainty information has previously been adopted in thermal remote sensing of sea and land surface temperatures at level 2/CDR e.g., [19], and here we discuss how it can also apply at level 1 (and support the tripartite calculations for CDRs).

We categorise each of the effects, $k$, into one subset according to their correlation structure. The set of independent effects, $\left\{i\right\}$, contains effects where there is no (or negligible) correlation of error between pixels. White noise from statistical fluctuations in detectors or amplifier circuits on the sensor would fall into this category. At the other extreme would be common effects, $\left\{m\right\}$, which cause the same unknown error in all pixels (at least approximately on some large scale). The uncertainty from common effects can be characterised across an image by a single uncertainty estimate. Between these extremes are structured effects, $\left\{s\right\}$, which can be either random or systematic in their physical origin, whose defining characteristic is that they cause errors that are partially shared across at least some pixels in the image. This tripartite classification is discussed in detail in [15].

Equations (10) to (15) can all be applied to independent and structured effects separately. For example, the modification of eq. 11 for independent effects is:

$${\mathit{S}}_{c,i}^p={\displaystyle \sum }_j{\displaystyle \sum }_{i|j}{\mathit{C}}_c^{p,j}{\mathit{U}}_c^{p,i}{\mathit{U}}_c^{p,i}{\mathit{C}}_c^{p,j}$$

(17)

which uses the fact that ${\mathit{R}}_c^{p,i}=\mathit{I}$ for independent effects by definition. (An extension is made to the notation, with the second subscript, $i$, of ${\mathit{S}}_{c,i}^p$ indicates the effects class if it is a subset of all.) Covariances may be added, and ${\mathit{S}}_c^p={\mathit{S}}_{c,i}^p+{\mathit{S}}_{c,s}^p+{\mathit{S}}_{c,m}^p$, while for a particular pixel the total uncertainty is the root sum of squares of the three component uncertainties.

5. Methods 3—Practical Computation

Calculating the “CURUC” equations presented in the previous section can be practically challenging because of the size and number of dimensions of the arrays involved. A literal implementation of the equations with each matrix in its own array is inefficient with regard to computer memory and the number of required floating point operations. This section presents the principles of a more efficient formulation of the equations for practical computation, and illustrates the principles with reference to specific Python functions that implement these.

Consider a single orbit of data from a radiometer with $n_c$ channels, $n_l$ scanlines, $n_e$ elements per scanline. Without any subsampling, the array containing the cross-line correlation matrix for each element and channel, ${\mathit{R}}_l^{e,k=K}$, has dimensions $\left[n_c,n_e,n_l,n_l\right]$ for a single effect, $k=K$. To store this array for $n_k$ effects requires $n_k{n}_c{n}_e{n}_l^2$ elements. Independent and common effects need not be calculated via CURUC equations, since their calculation simplifies, so the CURUC calculation can be restricted to the number of structured effects, $n_s$, giving $n_s{n}_c{n}_e{n}_l^2$ correlation coefficients to be stored.

Turning to the sensitivity matrices, the simplest approach is to store the sensitivities per structured effect. (This means some stored values are repeated, since the sensitivity coefficients relate to input quantities of the measurement equation rather than effects, and there can be more than one effect per input quantity.) A similar comment applies to the uncertainty matrices.

The calculation of ${\mathit{S}}_l^{e,s}={\mathit{C}}_l^{e,s}{\mathit{U}}_l^{e,s}{\mathit{R}}_l^{e,s}{\mathit{U}}_l^{e,s}{\mathit{C}}_l^{e,s}$ for all channels, effects and elements requires $4n_c{n}_s{n}_e$ matrix multiplications of $n_l\times n_l$ matrices. Optimised matrix multiplication algorithms exist and are implemented in numerical packages; for example, the Python NumPy library uses the Basic Linear Algebra Subprograms (BLAS), which implements the Coppersmith–Winograd algorithm [31]. Even so, the scale of the matrix computations here necessitates further optimisation: a literal implementation of the CURUC equations with matrix multiplication routines would be too slow for practical use.

In our implementation of the CURUC equations, we take advantage of $\mathit{C}$ and $\mathit{U}$ being diagonal. This property allows us to rewrite the matrix multiplication as an elementwise multiplication as follows. In the matrix multiplication ${\mathit{M}}_3={\mathit{M}}_1{\mathit{M}}_2$, where all matrices are square $n_M\times n_M$, the entries of ${\mathit{M}}_3$ are by definition given by ${\left[{\mathit{M}}_3\right]}_{v,w}={\displaystyle \sum }_{\chi =1}^{n_M}{\left[{\mathit{M}}_1\right]}_{v,\chi }{\left[{\mathit{M}}_2\right]}_{\chi ,w}$. If ${\mathit{M}}_1$ is diagonal, ${\left[{\mathit{M}}_1\right]}_{v,w}=0$ when $v\ne w$, such that ${\left[{\mathit{M}}_3\right]}_{v,w}={\left[{\mathit{M}}_1\right]}_{v,v}{\left[{\mathit{M}}_2\right]}_{v,w}$. Defining ${\mathit{d}}_{M_1}$ as the column vector containing the diagonal values of ${\mathit{M}}_1$, we can construct the $\left[n_M\times n_M\right]$ matrix ${\mathit{D}}_{M_1}=\left[{\mathit{d}}_{M_1}{\mathit{d}}_{M_1}\dots \right]$, such that we have ${\left[{\mathit{M}}_3\right]}_{v,w}={\left[{\mathit{D}}_{M_1}\right]}_{v,w}{\left[{\mathit{M}}_2\right]}_{v,w}$. Thus, we have rewritten the matrix multiplication as an elementwise multiplication, which can be parallelised much more easily. In practice, the final step can be calculated without actually constructing ${\mathit{D}}_{M_1}$ (in NumPy this process is called broadcasting). By storing only the diagonal values, the dimensionality of ${\mathit{C}}_l^{e,s}$ and ${\mathit{U}}_l^{e,s}$ is reduced by one and their size by a factor $n_l$.

In the calculation of ${\mathit{S}}_l^{e,s}$, ${\mathit{S}}_e^{l,s}$, ${\mathit{S}}_c^{p,s}$, we recognise that ${\mathit{U}}_l^{e,s}$, ${\mathit{U}}_e^{l,s}$, and ${\mathit{U}}_c^{p,s}$ are all arrays containing the same set of uncertainties, differing only in order. The arrays of the diagonals of ${\mathit{U}}_e^{l,s}$ and ${\mathit{U}}_c^{p,s}$ can be constructed as views of the same memory as ${\mathit{U}}_l^{e,s}$, avoiding duplication. The same applies to the diagonal sensitivity matrices, ${\mathit{C}}_c^{p,s}$, etc.

Even with those optimisations, a typical orbit file may remain too large for calculating the complete arrays as described above. The final optimisation is to approximate the averages across lines and elements, by sub-sampling these dimensions. For example, an AVHRR GAC orbit has $n_e=409$and $n_l\approx 12000$. For $n_c=5$ channels and $n_s=5$structured effects, the full set of ${\mathit{R}}_l^{e,s}$ alone would have $1.5\times {10}^{12}$ elements, which would require 5.9 TB of RAM to store in single IEEE 754 floating point precision; peak memory use in the calculation of ${\left[{\mathit{r}}_e\right]}_{{\Delta }_e}$ would be several times higher. We reduce the pixels across which the covariance matrices are averaged to every 10th element and every 50th line to reduce the RAM requirement to 236 MB. The validity of the sampling approximation for the estimated uncertainty information has to be verified for each case, and depends on the length scale of the effects and the data volume of the sensor. An alternative to sampling would be to calculate the length scales for shorter segments rather than full orbits, but that limits the calculation of correlation to scales shorter than the selected segment.

The optimisation described here has been implemented as a set of Python routines, built using the numpy [32] and numexpr libraries.

6. Results—Example from HIRS

We illustrate the results of applying the methods described above with the example of the High-Resolution Infrared Sounder (HIRS). The HIRS is a 20-channel radiometer with a heritage dating to 1975 [33], and it has been used in numerical weather prediction and analysis, and for climate data records [34]. A companion paper [18] describes the uncertainty properties of HIRS data in more detail. Figure 1 shows an extract of HIRS data, for channel 7 which measures radiance around 13.34 μm, which means it detects infra-red emission from clouds and carbon dioxide in the atmosphere. Also shown are the calculations of the uncertainty attributed to independent, structured and common effects.

In Figure 1, the uncertainty fields show variation with respect to the image brightness temperature. Uncertainty expressed in brightness temperature tends to be larger for the colder parts of the image where high cloud tops are present, since the non-linearity of the relationship between radiance and brightness temperature (Planck function) means a radiance error of a given size corresponds to a greater error in brightness temperature. This is the reason that the common uncertainty has a spatial structure in Figure 1: the common error is constant in radiance. The mean value of the uncertainty from independent effects differs significantly between calibration cycles. (The stripes of missing radiance along the orbit are the times when HIRS undertakes a calibration cycle.) This component of uncertainty characterises the noise of the instrument. The noise is estimated from the variability of the signal when looking at a uniform calibration target. The noisiness of instruments does change over time, and inferring the noise from each calibration cycle is appropriate.

Figure 2 shows the cross-element error covariance matrix from structured effects, ${\mathit{S}}_e^{l_A,s}$, evaluated for the scan line, $l_A$, that includes pixel ‘A’ in Figure 1. The variation in uncertainty reflected in this matrix is associated with the differential impact of structured errors on pixels with different brightness temperatures. The brightness temperature uncertainty tends to be greater for the colder pixels (which are on the right-hand side of ‘A’ in Figure 1, and which correspond to element numbers from 0 to 28 in Figure 2). The error covariance matrix has no obvious diagonal, reflecting the fact that the dominant origin of the uncertainty from structured effects is uncertainty in calibration, which tends to cause errors of the same sign in all radiances.

Figure 3 shows the cross-line error covariance for elements with the same element number as pixel ‘A’ in Figure 1. The scanline number 0 corresponds to the north (top) of the image. The variations in error covariance are again dominated by brightness temperature variations in the image, although the partial decorrelation of errors between different calibration cycles can also be seen (e.g., in the block-diagonal effect obvious as a boundary between row/column 6 and 7 in Figure 3).

The matrices shown in Figure 2 and Figure 3 have not been averaged across the image to give a representative estimate. Since a significant part of their evident structure relates to variations in the pixel brightness temperature, it follows that the error correlation matrix is more stable, and, as described earlier, could be used in conjunction with the per-pixel uncertainty fields to estimate the error covariance for any given pair of pixels more precisely. However, it is not practical for a full orbit to include a cross-line error correlation matrix, and the further approximation using Equation (16) must be used. Figure 4 presents the cross-line error correlation function thus derived for the image of Figure 1, not just for pixel ‘A’, but sampled and averaged across the whole image.

HIRS data are often used in multi-channel atmospheric sounding and for data assimilation, and so the cross-channel error correlations are important. Figure 5 shows an example obtained for the image in Figure 1. The matrix shows the error correlation between each pair of the 19 HIRS infra-red channels for independent effects. For this sensor, errors that are independent between pixels in an image (spatially) turn out to be partially correlated between channels, and so the off-diagonal terms of the cross-channel error correlation matrix are non-zero. There are two groups of channels that have cross-channel correlated noise between the channels in each group, but not between the groups, as can be seen from the block-diagonal form.

7. Conclusions

Uncertainty in derived climate data records and geophysical products comes in part from uncertainty in the satellite radiances used in the geophysical retrieval. Information about radiance uncertainty is generally not present in radiance products (level 1 data), compelling Earth observation practitioners doing geophysical retrievals to invent radiance uncertainty estimates or use other ad hoc methods to attribute uncertainties to results (if they do so at all). This situation falls short of best scientific and metrological practice, in which uncertainties at all stages are understood, characterised and propagated through to the final results. Although much useful science and many applications can be done without rigorous end-to-end treatment of uncertainty, where significant policy and liability questions rest on the integrity of observations, the trustworthiness and traceability of uncertainty estimates must be demonstrated.

However, providing uncertainty information at level 1 is not trivial. Radiance uncertainties are not simple. Uncertainty varies significantly from pixel-to-pixel, errors are correlated across the spatio-temporal domain in complex patterns that have their roots in instrument design and operation, and traditional consideration of “random vs. systematic errors” does not enable adequate characterisation of errors in satellite imagery. Moreover, solutions have to be practical for data producers and users, while still enabling a proper evaluation of uncertainty in derived products to be obtained. The full error covariance information for a multi-channel satellite image, even if we can estimate it, cannot be communicated to users because of the data volume implied.

This paper sets out one set of approaches that address the non-trivial aspects of providing users of radiance products with uncertainty information that is as simple as possible, but not simpler than necessary. Many details have inevitably been presented, so it is worth reviewing here the high-level concepts of the proposed solution.

First, the new uncertainty information available to users within a product following the methods presented here includes:

• Per-pixel, per-channel magnitude of radiance uncertainty
• Per-product, per-channel cross-element and cross-line radiance error correlation coefficients, and
• Per-product cross-channel radiance error correlation

Second, this uncertainty information is presented for a tri-partite classification of uncertainty components that adequately respects the complexity of errors in satellite imagery. The components describe the uncertainty that arises from errors that are independent (uncorrelated between pixels of a given channel), structured (patterned across an image) and common (systematic across a whole image).

Third, this set of uncertainty information gives users of radiance data a new capability to propagate uncertainty in a standardised, rigorous way through their retrieval process to derive uncertainty in their geophysical products. This includes where (i) the retrieval is multi-channel and therefore sensitive to any cross-channel correlations in errors, and (ii) where the target products aggregate in space and time and simplistic treatment radiance errors as “noise” would generally give under-estimates of uncertainty in the results. Propagation of the uncertainty information can be done using analytic or Monte Carlo methods [22].

Fourth, although some increase in data volume is necessary, it is minimised. (i) Uncertainty information need not be held to the same numerical precision as the radiances themselves, and so the per pixel provision of uncertainty information can be done such that the data volume is no more than doubled. (ii) Providing correlation functions as a function of pixel separation and the cross-channel error covariance matrices per radiance product has a small associated cost in data volume.

The framework presented above allows derivation of summary uncertainty (error covariance) information in radiance data based on systematic combination of instrumental error sources—a “bottom-up” approach. Is it reasonable to expect such an uncertainty analysis to become common practice? Space agencies already determine much of the necessary understanding of instrumental errors during instrument design, build and testing. Taking this instrument understanding forward into radiance uncertainty information for level-1 products, using the framework presented here or one similar, is therefore a value-adding extension of a much greater effort of characterisation that is already undertaken by agencies.

This work has been undertaken within the project Fidelity and Uncertainty in Climate data records from Earth Observation (FIDUCEO, www.fiduceo.eu), and the level 1 products being prepared at time of writing in that project have additional features, where relevant, to provide convenience to users. We will re-calibrate the satellite radiances to improve multi-sensor series. The products gather radiance, radiance uncertainty, viewing geometry and geolocation data in a self-described, compressed netCDF format. Lastly, at the scale of level-1 data collections, duplicate data are suppressed and orbits are consolidated. A standardised format for data across different sensor series allows tools to assist users with handling the uncertainty information to be designed, which will be the subject of a future paper.

Earth observation provides the largest and most globally descriptive set of information about the changing state of Earth’s environments. Proper interpretation of and inference from satellite data should be done in the context of well quantified uncertainties, to avoid wrong conclusions and bad decision-making. Practical delivery of the uncertainty information required presents a significant challenge and the methods in this paper are presented as candidate approaches to address this challenge.