An Algorithm for In-Flight Spectral Calibration of Imaging Spectrometers

Abstract

Accurate spectral calibration of satellite and airborne spectrometers is essential for remote sensing applications that rely on accurate knowledge of center wavelength (CW) positions and slit function parameters (SFP). We present a new in-flight spectral calibration algorithm that retrieves CWs and SFPs across a wide spectral range by fitting a high-resolution solar spectrum and atmospheric absorbers to in-flight radiance spectra. Using a maximum a posteriori optimal estimation approach, the quality of the fit can be improved with a priori information. The algorithm was tested with synthetic spectra and applied to data from the APEX imaging spectrometer over the spectral range of 385–870 nm. CWs were retrieved with high accuracy (uncertainty <0.05 spectral pixels) from Fraunhofer lines below 550 nm and atmospheric absorbers above 650 nm. This enabled a detailed characterization of APEX’s across-track spectral smile and a previously unknown along-track drift. The FWHMs of the slit function were also retrieved with good accuracy (<10% uncertainty) for synthetic spectra, while some obvious misfits appear for the APEX spectra that are likely related to radiometric calibration issues. In conclusion, our algorithm significantly improves the in-flight spectral calibration of APEX and similar spectrometers, making them better suited for the retrieval of atmospheric and surface variables relying on accurate calibration.

1. Introduction

Imaging spectroscopy from satellite and airborne platforms is playing an increasingly important role in Earth sciences to retrieve information on surface and atmospheric variables [1,2]. The retrieval of atmospheric and surface variables relies critically on accurate spectral calibration to properly account for the influence of spectrally-varying absorption and scattering in the atmosphere and at the surface [3,4,5]. In particular, retrievals that exploit narrow spectral features require highly accurate determinations of center wavelength positions (CW) and instrument slit function parameters (SFP). For example, spectral shifts as small as 0.01 spectral pixels can result in noticeable biases in atmospheric trace gases retrieved by remote sensing [6]. Damm et al. [7] showed that also the retrieval of chlorophyll fluorescence in the spectrally-narrow absorption lines of atmospheric oxygen is very sensitive to spectral shifts. In-flight spectral calibration (IFC) methods have been developed to detect changes in spectral calibration occurring during instrument operation. These methods determine CWs and SFPs using an artificial light source together with a filter or taking advantage of sharp signatures in the observed spectrum caused by Fraunhofer lines or atmospheric absorption (e.g., [8,9,10]).

One airborne observatory with a comprehensive infrastructure for on-ground and in-flight calibration is the Airborne Prism Experiment (APEX) [2,11]. APEX is an airborne dispersive push broom imaging spectrometer developed by a Swiss-Belgium consortium on behalf of the European Space Agency (ESA) [2]. It measures backscattered solar irradiance in the visible and near-infrared (VNIR) and in the short-wave infrared (SWIR) with a high spatial and spectral resolution (see

Table 1. Selected characteristics of the APEX imaging spectrometer in its two channels [2].

Characteristic	VNIR Channel	SWIR Channel
Spectral range	372–1015 nm	940–2540 nm
Spectral bands (mode)	114 (binned), 334 (unbinned)	198 (nominal)
Spectral sampling interval	0.45–7.5 nm	5.0–10.0 nm
Spectral resolution (FWHM)	0.86–15.0 nm	7.4–12.3 nm
Center wavelength accuracy	≤0.2 nm for a single flight
Spatial pixels (across track)	1000
Field of view	28.10${}^{\circ }$
Instantaneous field of view	0.028${}^{\circ }$
Spatial resolution	2.5 m at 5000 m above ground
Signal to noise (SNR)	625 (average over 50% reflecting target, Sun zenith: 24.4${}^{\circ }$)

for spectrometer characteristics). The instrument is calibrated regularly in ground calibration facilities [12,13], and an in-flight calibration facility can be used before and after data acquisition to improve the spectral calibration using a quartz tungsten halogen lamp and four spectral filters [14]. D’Odorico et al. [14] determined APEX’s spectral calibration parameters for different flight altitudes, i.e., different temperature and pressure conditions, using oxygen and water absorption lines, thus focusing on the near- and short-wave infrared. They also determined the sensor distortion in the across-track direction, known as spectral smile. Although the in-flight calibration system provides a unique capability to assess sensor performance in flight, the limited spectral range and sampling density hinder detailed characterizations of spectral non-uniformities. To improve the spectral calibration in the visible, Fraunhofer lines have been used to determine APEX’s CWs and SFPs [3,15]. The suggested approach exploits measured radiances and is based on differential optical absorption spectroscopy (DOAS) as implemented, for example, in the QDOAS software package [16]. Using this approach, Popp et al. [3] concluded that the full width at half maximum (FWHM) of the slit function was larger in the nitrogen dioxide (${\mathrm{NO}}_2$) fitting window (470–510 nm) compared to the laboratory calibration without being able to provide an explanation for this observation. In the current Level 1 product, CWs in the VNIR channel are determined from known spectral signatures in the oxygen A-band at 760 nm, the water vapor band at 840 nm and from laboratory measurements of the sensitivity of wavelength shifts due to pressure changes in the spectrometer. The CW uncertainty was estimated at 0.2 nm [13], but a recent re-analysis suggests a larger uncertainty of at least 0.4 nm. SFPs are taken from the laboratory measurements [13]. Although shifts in CWs during flight have already been identified based on IFC facility measurements and DOAS analyses, there has been no algorithm that consistently determines such shifts, as well as changes in SFPs across the entire wavelength range with high accuracy.

In this paper, we present a comprehensive and flexible in-flight spectral calibration algorithm that uses Fraunhofer lines and atmospheric trace gas absorption lines (see Figure 1) to retrieve CWs and SFPs within a single fit across the visible and near-infrared spectral range. CWs, SFPs, spectral albedo and a radiance offset are parametrized by cubic Hermite splines (C-splines), which in contrast to low-order polynomials used in DOAS approaches, are suitable for wide spectral windows. The algorithm is able to incorporate a priori information, such as IFC facility measurements and laboratory calibration, to regularize the retrieval where Fraunhofer lines and atmospheric absorbers do not provide sufficient information. Our particular objectives are to: (i) describe the algorithm and its application to APEX-like instruments; (ii) evaluate its performance using synthetic and real measurement data; and (iii) discuss performance, limitations and potential future improvements of the algorithm.

2. Data

2.1. Synthetic Data

To analyze the performance of the algorithm, synthetic spectra were generated for realistic observation scenarios with a radiative transfer model and convolved with an APEX-like slit function. High-resolution radiance spectra were simulated with the libRadtran radiative transfer model (Version 2.0-beta) [18,19] for an airborne nadir-viewing spectrometer operating at an altitude of 5000 m. The spectral range was 385–900 nm with a 0.001-nm spectral resolution. Solar irradiance was taken from the Kitt Peak Flux Atlas 2005 [20]. Atmospheric profiles were mid-latitude summer profiles [21]. Absorber optical depth profiles of water vapor (H${}_2$O), oxygen (O${}_2$), carbon dioxide (CO${}_2$) and methane (CH${}_4$) were pre-computed with the atmospheric radiative transfer simulator (ARTS, Version 2.2) [22] using the HITRAN 2012 database [23]. Optical depth profiles of broad-band absorbers were pre-calculated using temperature-dependent absorption cross-sections of nitrogen dioxide (NO${}_2$) (Vandaele et al. [24]), ozone (O${}_3$) (Serdyuchenko et al. [25]) and tetraoxygen (O${}_4$) (Thalman and Volkamer [26]) from laboratory measurements. The pre-computed optical depth profiles were used by libRadtran. Atmospheric scattering and absorption were simulated with libRadtran’s default aerosol parametrization [27] and Rayleigh scattering [28]. Surface reflectance was modeled assuming Lambertian reflection and was taken from the ASTER spectral library for “green grass”, “asphalt paving” and “metal roofing” (Figure 2c) [29]. The model parameters are summarized in

Table 2. Parameters used for the simulation of synthetic spectra.

Parameter
Atmosphere	mid-latitude summer
Atmospheric absorber (near infrared)	${\mathrm{H}}_2\mathrm{O}$, ${\mathrm{O}}_2$, ${\mathrm{CO}}_2$, ${\mathrm{CH}}_4$ from HITRAN2012 [23]
Atmospheric absorber (visible)	${\mathrm{NO}}_2$ [24], ${\mathrm{O}}_3$ [25] and ${\mathrm{O}}_4$ [26]
Atmospheric extinction	libRadtran default Rayleigh and Mie scattering, no Raman scattering
Solar reference spectrum	Kitt Peak Flux Atlas 2005 [20]
Solar zenith angle	23${}^{\circ }$
Spectral range	385–900 nm (0.001-nm resolution)
Surface elevation	0 m
Surface reflectance	“asphalt paving”, “green grass” and “metal roofing” from ASTER
	spectral library [29]
Instrument altitude	5000 m
Instrument viewing zenith angle	0${}^{\circ }$ (nadir)

.

The high-resolution spectra were convolved with an APEX-like slit function. As a reference, we used CWs and FWHMs from the APEX Level 1 processor. The APEX slit function is to a very good approximation a Gaussian curve with FWHM increasing from 0.90–8.25 nm over the selected spectral range [13] (Figure 2b). APEX CW positions can shift due to pressure changes of the nitrogen gas in the spectrometer kept at 200 hPa above ambient pressure. The pressure affects the index of refraction of the nitrogen gas and, thus, the dispersion at the prism, resulting in a wavelength shift. The dependency of the shift on pressure was determined in the laboratory. The in-flight pressure is measured with an accuracy of about 15 hPa. Figure 2a shows the wavelength shift caused by a pressure change of only 15 hPa. For the synthetic spectra, CW positions were shifted from the reference using the wavelength pressure sensitivity measured in the laboratory (Figure 2a). APEX FWHMs were also measured in the laboratory (Figure 2b) and do not show any significant sensitivity to pressure changes. Nonetheless, FWHMs as large as twice the laboratory values were found for in-flight spectral calibration near 500 nm [3,15]. To simulate this effect, reference FWHMs were scaled by a factor between 0.5 and 2.0 for the tests. Finally, we added a radiance offset to represent potential issues with dark current calibration or stray light in the instrument. The offset was simulated as a parabolic curve:

$$O\left(i\right)=-0.0002i(i-316)\mathrm{mW}·{\mathrm{m}}^{-2}·{\mathrm{nm}}^{-1}·{\mathrm{sr}}^{-1}$$

(1)

where i is the band number of the spectral pixel. The radiance offset has a maximum of 5 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ near 508 nm, i.e., about 10% of the radiance simulated over “green grass”.

2.2. Real Data from the APEX Imaging Spectrometer

APEX measurements analyzed here were obtained during an aircraft campaign (altitude: 6646 m above sea level) over Zurich-Opfikon, Switzerland, on 30 August 2013 around noon. Data were acquired in the so-called “unbinned mode”, which allows measuring at the highest possible spectral resolution offered by the instrument (Table 1) to facilitate atmospheric applications. Figure 3a shows the true color composite of this test dataset consisting of 1000 × 1000 pixels in the across- and along-track direction, respectively. The dataset covers diverse surface types, including urban areas, forest and farmland, making it suitable for testing the algorithm. Because of detector saturation in the NIR bands, about 40% of the pixels had to be removed, mainly over farmland. The signal saturation in the NIR was intentionally accepted to increase the signal-to-noise ratio in the visible bands used for NO${}_2$ remote sensing. To increase the signal-to-noise ratio, 40 × 20 pixels (about 100 × 100 m${}^2$) were spatially binned resulting in 25 × 50 allocated pixels (Figure 3b). The 25 lines in along-track direction are marked in Figure 3. The numbering of lines is reversed with respect to the across-track position, i.e., Line 0 corresponds to across-track Positions 960–999. For this dataset, the measurement precision ${\sigma }_{\epsilon }$ of a single pixel was estimated by the APEX Level-1 processor to 1.31 ± 0.24 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$. The uncertainty is reduced by spatial binning depending on the degree of randomness in the instrument noise.

3. Spectral Calibration Algorithm: General Framework and Application to APEX

3.1. General Framework

The spectral calibration algorithm is designed without inferring a specific instrument to make it as generally applicable as possible for different types of spectrometers. It comprises two main components, a forward model and a retrieval module.

3.1.1. Forward Model

The core of the algorithm is a simple forward model that calculates the spectral radiance measured by a downward viewing spectrometer without requiring a radiative transfer solver. The spectral radiance $L\left(\lambda \right)$ is solar irradiance $E_0\left(\lambda \right)$ that is scattered and absorbed in the atmosphere and reflected at the surface. Ignoring scattering of light in the direction of the instrument, the at-sensor radiance can be calculated as (e.g., Liou [30], Chapter 7):

$$L\left(\lambda \right)=\rho \left(\lambda \right){\mu }_0{E}_0\left(\lambda \right)e^{-\tau \left(\lambda \right)}+N\left(\lambda \right)$$

(2)

where $\rho \left(\lambda \right)$ is the surface reflectance, ${\mu }_0$ is the cosine of the solar zenith angle, $\tau \left(\lambda \right)$ is the optical depth of the atmosphere and $N\left(\lambda \right)$ is the measurement error. Equation (2) calculates the continuous radiance spectrum. In practice, the discrete radiance needs to be computed on a high-resolution wavelength grid with indices $i_{\mathrm{hr}}$. The forward model convolves the high-resolution radiance spectrum $L\left(i_{\mathrm{hr}}\right)$ with the instrument slit function and projects the radiances onto the low-resolution pixels $i_{\mathrm{lr}}$ of the spectrometer:

$$L\left(i_{\mathrm{lr}}\right)={\left[\rho \left(i_{\mathrm{hr}}\right){\mu }_0{E}_0\left(i_{\mathrm{hr}}\right)e^{-\tau \left(i_{\mathrm{hr}}\right)}\right]}_{i_{\mathrm{lr}}}+N\left(i_{\mathrm{lr}}\right)$$

(3)

where ${[·]}_{i_{\mathrm{lr}}}$ indicates the convolution with the slit function and the projection onto the spectral pixels $i_{\mathrm{lr}}$.

The spectral calibration is described by the CWs and SFPs, which change with pixel number. In the forward model, we calculate the CWs as:

$$\lambda \left(i_{\mathrm{lr}}\right)={\lambda }_0\left(i_{\mathrm{lr}}\right)+\Delta \left(i_{\mathrm{lr}}\right)$$

(4)

where ${\lambda }_0\left(i_{\mathrm{lr}}\right)$ are the reference CWs, for instance from the laboratory calibration, and $\Delta \left(i_{\mathrm{lr}}\right)$ is the wavelength shift from the reference. The slit function at each spectral pixel can be parametrized by a set of parameters, such as the FWHM and a shape factor, as suggested, for example, by Beirle et al. [31].

The forward model needs to include several absorbing gases. The optical depth of gas j can be written as:

$${\tau }_j={\int }_l{c}_j\left(l\right){\sigma }_j(\lambda ,l)dl\left(\lambda \right)$$

(5)

with optical path l, absorber number density $c_j$ and absorption cross-section ${\sigma }_j$. The absorption cross-section depends on temperature and pressure changing along the optical path. In the visible and near-infrared, important absorbing gases are oxygen (${\mathrm{O}}_2$), water (${\mathrm{H}}_2\mathrm{O}$), ozone (${\mathrm{O}}_3$), tetraoxygen (${\mathrm{O}}_4$) and nitrogen dioxide (${\mathrm{NO}}_2$).

In the near-infrared, absorption is mainly caused by ${\mathrm{O}}_2$ and ${\mathrm{H}}_2\mathrm{O}$ molecules. The absorption is caused by rotational and vibrational transitions, and the resulting absorption lines are narrow and depend strongly on temperature and pressure. For these two reasons, optical depth needs to be calculated as the sum of layer optical depths:

$${\tau }_j\left(\lambda \right)=\sum _k{x}_{jk}{A}_{jk}\left(\lambda \right){\tau }_{jk}^{\mathrm{ref}}\left(\lambda \right)$$

(6)

where each vertical layer k has constant temperature and pressure [32]. In each layer, a reference optical depth ${\tau }_{jk}^{\mathrm{ref}}$ is multiplied with a scaling factor $x_{jk}$ and layer air mass factor (AMF) $A_{jk}$. Since ${\mathrm{O}}_2$ and ${\mathrm{H}}_2\mathrm{O}$ absorb in the near infrared where atmospheric scattering is small, $A_k$ can be approximated by the geometric AMF, i.e.,:

$$A_k=\left\{\begin{array}{cc}{\mu }_0^{-1}+{\mu }^{-1}\hfill & \mathrm{if}{z}_k<H\hfill \\ {\mu }_0^{-1}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(7)

with ${\mu }_0$ and μ being the cosine of the solar and viewing zenith angle and H is the altitude of the instrument.

In the visible, atmospheric scattering is larger, and the optical path cannot be approximated by the geometric AMF. However, absorber optical depths are small, and their cross-sections depend only slightly on temperature and pressure. In this case, Equation (5) can be approximated as:

$${\tau }_j=S_j{\sigma }_j\left(\lambda \right)$$

(8)

where S is the slant column density (SCD) [6]. As a consequence, the surface reflectance $\rho \left(\lambda \right)$ in Equation (3) also includes the atmospheric albedo in the visible. We also assume that the slant column density is constant across the spectral window. The optical depths of ${\mathrm{NO}}_2$, ${\mathrm{O}}_3$ and ${\mathrm{O}}_4$ can reasonable well be fitted using this approximation.

Substituting Equations (6) and (8) with $n_{\mathrm{vis}}$ absorbers in the visible and $n_{\mathrm{nir}}$ absorbers in the near infrared in Equation (3), the final equation for calculating the spectral radiance $L\left(i_{\mathrm{lr}}\right)$ becomes:

$$L\left(i_{\mathrm{lr}}\right)={\left[{\mu }_0{E}_0\left(i_{\mathrm{hr}}\right)s\left(i_{\mathrm{hr}}\right)exp\left(-\sum _{j=1}^{n_{\mathrm{vis}}}{S}_j{\sigma }_j\left(i_{\mathrm{hr}}\right)-\sum _{j=1}^{n_{\mathrm{nir}}}\sum _{k=1}^{n_l}{x}_{jk}{A}_k{\tau }_{jk}^{\mathrm{ref}}\left(i_{\mathrm{hr}}\right)\right)\right]}_{i_{\mathrm{lr}}}+O\left(i_{\mathrm{lr}}\right).$$

(9)

The term $s\left(i\right)=\rho \left(i\right)e^{-\tau \left(i\right)}$ is the “albedo” term that combines surface reflectance and atmospheric extinction. The term $O\left(i\right)$ is a radiance offset that is the fraction of the error term $N\left(i\right)$ that varies smoothly with wavelength and can be caused, for example, by errors in the dark current calibration or stray light.

Spectral calibration algorithms based on DOAS are using low-order polynomials to approximate the wavelength dependency of the albedo, wavelength shifts and variations in slit function parameters across the wavelength window [6]. Figure 2c shows examples of spectral surface reflectance which need to be parametrized by the forward model. Most of these curves (e.g., albedo over green grass) cannot be parametrized accurately by a low-order polynomial. Parametrization errors in the wide window will become too large. The QDOAS software overcomes this limitation by dividing wide windows into smaller subwindows and performing the spectral calibration separately for each subwindow [16]. However, this approach can lead to issues with discontinuities at the boundaries between the subwindows. It is also not possible to use information from neighboring subwindows to constrain the retrieval in regions of the spectrum where clear spectral features are absent. Therefore, we propose to use splines instead of polynomials to approximate albedo, radiance offset, wavelength shifts and variations in slit function parameters in the forward model.

A spline consists of piecewise-defined polynomials and has a high degree of smoothness at the knots where the polynomials are connecting (e.g., [33]). Therefore, a spline is a continuous function suitable to approximate the spectral albedo or the slit function parameter over a wide spectral range. A spline $s\left(i\right)$ with coordinate i in the spectral direction is defined on $n+1$ knots:

$$\mathbf{t}=\{t_j:t_0<...<t_j<...<t_n\}$$

(10)

and can be written as the sum of piecewise polynomials:

$$s\left(i\right)=\sum _{j=0}^{n-1}{s}_j\left(i\right)$$

(11)

with:

$$s_j\left(i\right)=\left\{\begin{array}{cc}{P}_j\left(i\right)\hfill & \mathrm{for}{t}_j\le i<t_{j+1}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(12)

where $P_j\left(i\right)$ is a polynomial. There are various possibilities for defining $P_j\left(i\right)$. A simple, but for most applications, suitable variant is the Hermite form resulting in a cubic spline with the continuous first derivative, also called a C-spline. The polynomial can be written for knot vector $\mathbf{t}$ and control points $p_j$ as:

$$P_j\left(i\right)=(2{\xi }^3-3{\xi }^2+1)p_j+({\xi }^3-2{\xi }^2+\xi )m_j+(-2{\xi }^3+3{\xi }^2)p_{j+1}+({\xi }^3-{\xi }^2)m_{j+1}$$

(13)

with dimensionless number:

$$\xi =\frac{i-t_j}{t_{j+1}-t_j}$$

(14)

and three-point differences:

$$m_j=\frac{p_{j+1}-p_j}{t_{j+1}-t_j}+\frac{p_j-p_{j-1}}{t_j-t_{j-1}}$$

(15)

for $j=1,...,n-2$ and one-point differences at the end points. C-splines are defined locally, are continuous and have a continuous first derivative. An advantage of C-splines is that the control points $p_j$ are located on the spline curve and therefore directly describe the parameter of interest, for example a wavelength shift. This makes it easy to add a priori information. The number and positions of the knots can be selected based on the requirements of an accurate parametrization in the forward model.

3.1.2. Maximum A Posteriori Retrieval

For the retrieval, we define a measurement vector $\mathbf{y}$, whose elements are m spectral radiances computed by the forward model (Equation (9)) and write the model for all spectral bands as vector equation:

$$\mathbf{y}=\mathbf{F}(\mathbf{x},\mathbf{b})+\mathbf{e}$$

(16)

where $\mathbf{F}$ is the forward model that depends on the state vector $\mathbf{x}$ and the parameter vector $\mathbf{b}$. The state vector $\mathbf{x}$ contains parameters that are retrieved by the calibration algorithm. Its n elements are trace gas slant column densities and scaling factors, as well as the control points of the C-splines describing wavelength shift, slit function parameters, offset and albedo as a function of wavelength. The parameter vector contains parameters that are not retrieved by the algorithm, such as solar and viewing zenith angle or absorption cross-sections. The error vector $\mathbf{e}$ summarizes errors from the instrument and forward model.

To retrieve the spectral calibration for a given radiance spectrum, we find the maximum a posteriori (MAP) optimal estimate using the method described by Rodgers [34]. For a given measurement vector $\mathbf{y}$, the optimal state vector $\widehat{\mathbf{x}}$ minimizes of the following figure-of-merit function:

$${\chi }^2\left(\mathbf{x}\right)={(\mathbf{y}-\mathbf{F}\left(\mathbf{x}\right))}^T{\mathbf{S}}_{ϵ}^{-1}(\mathbf{y}-\mathbf{F}\left(\mathbf{x}\right))+{(\mathbf{x}-{\mathbf{x}}_a)}^T{\mathbf{S}}_a^{-1}(\mathbf{x}-{\mathbf{x}}_a)$$

(17)

where ${\mathbf{S}}_{ϵ}$ is the measurement error covariance matrix and ${\mathbf{x}}_a$ and ${\mathbf{S}}_a$ are the a priori state vector and error covariance matrix, respectively. Since the forward model is non-linear, we use the Gauss–Newton method to find $\widehat{\mathbf{x}}$ iteratively by:

$${\mathbf{x}}_{i+1}={\mathbf{x}}_i+{\widehat{\mathbf{S}}}_i\left({\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}(\mathbf{y}-\mathbf{F}({\mathbf{x}}_i,\mathbf{b}))-{\mathbf{S}}_a^{-1}({\mathbf{x}}_i-{\mathbf{x}}_a)\right)$$

(18)

with a posteriori error covariance matrix:

$${\widehat{\mathbf{S}}}_i={({\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}{\mathbf{K}}_i+{\mathbf{S}}_a^{-1})}^{-1}.$$

(19)

${\mathbf{K}}_i$ is the $m\times n$ Jacobian or weighting function matrix whose elements are the partial derivatives of the forward model with respect to the state vector. In other words, it describes the spectral response to variations in the state vector parameters. The Jacobian is calculated by finite-difference approximation in each iteration step. The iteration is stopped if:

$${({\mathbf{x}}_i-{\mathbf{x}}_{i+1})}^T{\widehat{\mathbf{S}}}_i^{-1}({\mathbf{x}}_i-{\mathbf{x}}_{i+1})<n$$

(20)

where n is the number of state vector elements.

The error in the retrieved state vector $\widehat{\mathbf{x}}$, i.e., its difference from the unobservable true state $\mathbf{x}$, can be expressed as:

$$\widehat{\mathbf{x}}-\mathbf{x}=(\mathbf{GK}-\mathbf{I})(\mathbf{x}-{\mathbf{x}}_a)+\mathbf{G}\mathbf{e}$$

(21)

with identity matrix $\mathbf{I}$ and gain matrix:

$$\mathbf{G}={({\mathbf{K}}^T{\mathbf{S}}_{\epsilon }^{-1}\mathbf{K}+{\mathbf{S}}_a^{-1})}^{-1}{\mathbf{K}}^T{\mathbf{S}}_{\epsilon }^{-1}.$$

(22)

The gain matrix is the sensitivity of the retrieval to the measurement and depends on measurement and a priori error covariance matrices. The first term on the right-hand side of Equation (21) is the smoothing error, and the second term is the retrieval error from instrument errors, model parameter errors and forward model errors [34].

The total retrieval error is estimated by the a posteriori error covariance matrix $\widehat{\mathbf{S}}$ as described by Equation (19) from measurement and the a priori error covariance matrix. It is the sum of smoothing error and the retrieval noise covariance matrix. However, the smoothing error is only described correctly if ${\mathbf{S}}_a$ is the covariance matrix of the true ensemble of state vectors (see Chapter 3 of Rodgers [34] for details). If the true ensemble is not known, $\widehat{\mathbf{S}}$ can be interpreted as the error of the smoothed true state.

If the measurement covariance matrix ${\mathbf{S}}_{\epsilon }$ is not known, it is possible to estimate it assuming a moderate quality of fit, i.e., ${\chi }^2=n+m$, and solving Equation (17) for ${\sigma }_{\epsilon }$ to obtain a single value for all spectral pixels [35].

For any retrieval algorithm, it is important to analyze the amount of information that can be retrieved from the measurement. This information is provided by the averaging kernel matrix $\mathbf{A}$, which is the sensitivity of the retrieval to the true state and is given by:

$$\mathbf{A}=\mathbf{G}\mathbf{K}.$$

(23)

The diagonal elements of the averaging kernel matrix are the degrees of freedom for signals $d_s$. The trace of the matrix is the total degrees of freedom [34].

3.2. APEX Spectral Calibration Algorithm

To demonstrate the application of the general framework presented above, we here describe the setup of the spectral calibration algorithm for an imaging spectrometer with the specifications of the APEX instrument. We estimate the expected instrument, model parameter and forward model errors and construct a suitable a priori state vector and covariance matrix. The algorithm was implemented in Python and consists of the forward model and retrieval module, as well as code for handling the APEX dataset and visualizing the results of the retrieval. It was built on top of the flexDOAS library, which is a Python library that we are currently developing for the flexible implementation of DOAS-like retrieval algorithms (available at [36]).

3.2.1. Forward Model

The forward model calculates spectral radiances with Equation (9) for a subset of 284 spectral pixels in the APEX VNIR channel between 385 and 870 nm (Band Numbers 32–316 in the unbinned mode). The bands have spectral sampling intervals between 0.47 and 6.11 nm. The choice of this spectral range was driven by the availability of high-resolution solar reference spectra and absorption cross-sections. Solar irradiance was taken from the Kitt Peak Flux Atlas 2005 with a 0.001-nm spectral resolution [20]. The high-resolution is required in the presence of strong absorbers with narrow absorption lines [32].

Atmospheric absorption of five gases was included in the forward model: ${\mathrm{H}}_2\mathrm{O}$, ${\mathrm{NO}}_2$, ${\mathrm{O}}_2$, ${\mathrm{O}}_3$ and ${\mathrm{O}}_4$. ${\mathrm{H}}_2\mathrm{O}$ and ${\mathrm{O}}_2$ reference optical depth profiles (${\tau }_{jk}^{\mathrm{ref}}$ in Equation (9)) were pre-computed with ARTS [22] using the HITRAN 2012 database [23] and mid-latitude summer profiles [21]. For the retrieval, we assumed a fixed profile shape; thus only a single scaling factor was used for each absorber. High-resolution cross-sections of ${\mathrm{NO}}_2$ (measured at a temperature of 220 K), ${\mathrm{O}}_3$ (223 K) and ${\mathrm{O}}_4$ (293 K) were taken from Vandaele et al. [24], Serdyuchenko et al. [25] and Thalman and Volkamer [26], respectively.

As described above, C-splines were used to parametrize wavelength shift, FWHMs, radiance offset and albedo. Except for albedo, these were parametrized with 29 equidistant control points and knots every five spectral pixels. Since the sampling interval increases with band number, the knot distance increases from 2.5 nm at 385 nm to 30 nm at 870 nm.

The knot spacing for the parametrization of surface reflectance and atmospheric extinction was chosen based on two opposing requirements: it needed to be small enough to avoid large parametrization errors and large enough to avoid tracing narrow absorption features not associated with these parameters. Figure 2c shows surface reflectances over “green grass”, “asphalt paving” and “metal roofing” from the ASTER database [29]. Based on these spectra, we selected a spline with 73 knots with a spacing of 10 nm between 385 and 480 nm, 5 nm between 480 and 750 nm and 15 nm between 750 and 860 nm.

The final state vector contains the control points of the splines and the five trace gas scaling factors. In total, the state vector has 165 elements.

3.2.2. Instrument and Forward Model Errors

The measurement covariance matrix includes instrument errors, model parameter errors and forward model errors that need to be estimated separately.

The APEX instrument error is a combination of sensor noise and radiometric calibration errors and has been estimated by the APEX Level 1 Processor [37]. The relative uncertainty is about 15% at 385 nm and decreases to about 5% for wavelengths larger than 450 nm. It is larger for small wavelengths, because the radiance signal is smaller for most surfaces. In contrast, the absolute error does not change strongly with wavelength and is about 1.31 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$.

Since the instrument error is partially random, it can be reduced by spatial binning. However, since the uncertainty has a random and systematic error component, the uncertainty of the binned spectra needs to be calculated as:

$${\sigma }_{\epsilon ,N}=\left(f_s+\frac{f_r}{\sqrt{N}}\right){\sigma }_{\epsilon }$$

(24)

where $f_s$ and $f_r$ are the fractions of systematic and random error ($f_s+f_r=1$), respectively, and N is the number of binned spectra. For example, binning 10 × 20 pixels reduces the error to 0.09 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ for pure random noise and to 0.46 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ for 70% random noise. The systematic error component in Equation (24) can be estimated by applying the spectral calibration to binned spectra where N is large enough such that the random component is negligible. The systematic error can then be estimated from the fitting residual assuming a moderate quality of fit (i.e., ${\chi }^2\approx m+n$, Equation (17)).

The forward model errors due to uncertainties in the solar reference spectrum and absorption cross-sections were estimated with a Monte Carlo approach, i.e., adding random noise to the reference spectra and calculating the measurement vector repeatedly using the a priori state vector. The uncertainty was then determined as the standard deviation of the ensemble. The uncertainty of the solar reference spectrum is between 0.1% and 1% [20], which corresponds to a forward model error between 0.004 and 0.040 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$. The uncertainty of the ${\mathrm{NO}}_2$, ${\mathrm{O}}_3$ and ${\mathrm{O}}_4$ absorption cross-sections is less than 5% [24,25,26] resulting in a model error of less than 0.006 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ for all wavelengths. The ${\mathrm{O}}_2$ and ${\mathrm{H}}_2\mathrm{O}$ reference optical depths are calculated by the ARTS model, which does not provide uncertainties. We simply assume an accuracy of 5%, which results in a forward model error of less than 0.030 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$.

The so-called Ring effect causes a filling-in of Fraunhofer lines and arises mainly due to rational Raman scattering (RRS) in the atmosphere [38]. Since the Ring effect is mainly important in the ultraviolet, the Ring effect was not implemented in our forward model. To determine the magnitude of the resulting forward model error, we calculated spectral radiances for the a priori state vector using the libRadtran radiative transfer model [18] with and without Raman scattering. The Ring spectrum, the difference between the two spectra, has a standard deviation of 0.350 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ at 400 nm and less than 0.050 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ above 450 nm. Since the error below 450 nm is quite large, its impact on the retrieved CWs and SFPs was also analyzed with synthetic spectra (see Section 4.3).

The C-splines cause a parametrization error that depends on the number and location of knots, as well as the parametrized curve. The knot spacing for wavelength shift and FWHM was chosen such that the spline can parametrize expected shifts and FWHM without any parametrization error. The offset spline fits the part of the error term that varies slowly with wavelength, and thus, its parametrization error is already described by the instrument error. For the albedo spline, knot spacing was chosen as a compromise between parametrization error and the number of knots using reference surface reflectances. For the three different surface reflectances, the parametrization errors were estimated between 0.038 and 0.174 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$.

Table 3. Estimates of expected instrument and forward model errors.

Error Source	σ (mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$)	Remark
APEX Instrument Error
single pixel	1.31	average instrument error
200 binned pixels	0.46	assuming 70% random noise
200 binned pixels	0.09	assuming 100% random noise
Forward Model Errors
Albedo parametrization (asphalt paving)	0.038	errors determined for the selected...
Albedo parametrization (green grass)	0.174	...knot spacing in the forward...
Albedo parametrization (metal roofing)	0.114	...model
Cross sections (${\mathrm{NO}}_2$, ${\mathrm{O}}_3$, ${\mathrm{O}}_4$)	<0.006	5% error on the cross section
Reference optical depth (${\mathrm{H}}_2\mathrm{O}$, ${\mathrm{O}}_2$)	<0.030	5% error on reference optical depth
Ring effect (<450 nm)	<0.350	standard deviation of the Ring spectrum...
Ring effect (>450 nm)	<0.050	...for a priori state vector
Solar reference spectrum	<0.004 and 0.040	0.1% and 1% uncertainty on data

summarizes measurement and forward model errors. The Ring error is the dominant component below 450 nm. Albedo parametrization errors are another important component that can become dominant for spatially-binned pixels.

3.2.3. A Priori State Vector Errors and Regularization

The MAP optimal estimate can be seen as a regularized maximum likelihood estimate. The optimal estimate is regularized by an a priori state vector and error covariance matrix in the second term on the right-hand side of Equation (17). Suitable a priori values need to be estimated for the five trace gases and the four splines describing wavelength shift, FWHM, offset and albedo.

A priori values of atmospheric absorbers were taken from the mid-latitude summer profiles. Their standard uncertainties were roughly estimated for their expected atmospheric variability to 10% for ${\mathrm{O}}_3$, 20% for ${\mathrm{NO}}_2$ and 30% for ${\mathrm{H}}_2\mathrm{O}$. ${\mathrm{O}}_2$ and ${\mathrm{O}}_4$ are very well known allowing for a small standard error of 3%.

For each C-spline, we estimated a priori mean and standard uncertainty of the control points. The albedo mean and standard uncertainty were calculated from the APEX surface reflectance product for measurements over Zurich [39]. A priori wavelength shifts were 0.0 nm with a standard uncertainty of 0.2 nm estimated from laboratory measurements [13]. For real data, laboratory measurements of the spectral smile effect (wavelength shift in the across-track direction) were additionally considered in the a priori shift. A priori FWHM is the laboratory calibration (FWHM${}_{\mathrm{lab}}$) with a standard uncertainty of 15%. The radiance offset was specified as 0.0 ± 5.0 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$. The standard uncertainty corresponds to about a third of the APEX instrument error (${\sigma }_{\epsilon }$ = 1.31 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$, see Table 3) or about 10% of the radiance over a dark surface (e.g., “green grass”). Overall, the offset standard deviation is quite large.

Some of the uncertainties of the a priori state vector elements are assumed to be correlated, which further regularizes the problem and makes the estimation of 165 individual state vector elements possible. Error covariances were considered for wavelength shift, FWHM and offset spline as:

$${\sigma }_{i,j}^2={\sigma }_i{\sigma }_{j}exp\left(-\frac{|t_i-t_j|}{L}\right)$$

(25)

where ${\sigma }_i$ and ${\sigma }_j$ are the standard deviations at knot $t_i$ and $t_j$, respectively, and L is the correlation length in units of spectral pixels. The correlation length for each spline was determined by minimizing the smoothing error (Equation (21)) for an ensemble of 100 state vectors. The motivation for this procedure comes from the fact that smoothing errors in CWs and FWHMs may introduce a systematic error in the retrieved variables and thus should be small. The ensemble of state vectors was created using three different surface reflectances, variable wavelength shifts due to pressure changes in the range ±15 hPa, FWHM scalings between 0.75 and 1.15 and with and without radiance offset described by Equation (1).

For each ensemble member, a measurement vector was computed using the state vector and forward model and adding measurement noise between 0.1 and 0.5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$. The retrieval was then done for each measurement vector and varying the correlation lengths between one and 10${}^6$ spectral pixels. We calculated standard deviation and mean bias to determine retrieval noise and smoothing error, respectively. The optimization was conducted for each spline separately once without retrieving the other splines and once with correlation lengths of the other splines set to the optimal value. Figure 4 shows retrieval noise and smoothing error assuming a measurement uncertainty of 0.2 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$. The retrieval noise is largest for small correlation lengths and decreases with larger correlation length, because with increasing correlation length, the solution is increasingly determined by the a priori information. On the other hand, the smoothing error is large for both small and large correlation lengths, because the retrieval under- and over-regularized in these cases, respectively. The correlation length has a minimum, which is used as the optimal value. The finally chosen values are listed in

Table 4. A priori state vector and covariance matrix.

Spline	${\mathit{x}}_{\mathit{a}}$	${\mathit{\sigma }}_{\mathit{a}}$	Correlation Length (L)
Wavelength shift	0.0 nm	0.2 nm	100 spectral pixels
Slit function	FWHM${}_{\mathrm{lab}}$	0.15 × FWHM${}_{\mathrm{lab}}$	100 spectral pixels
Offset	0 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$	5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$	1000 spectral pixels
Albedo	means and standard deviations of APEX surface reflectance dataset

.

For the albedo spline, we use only diagonal variances. In this case, the knot spacing determines the degree of regularization. For the other splines, the degree of regularization is mainly determined by the correlation length, as it is much larger than the distance between knots.

4. Results

In this section, we analyze the performance of the spectral calibration algorithm. First, we calculate the degrees of freedom, smoothing errors for the a priori state vector and error covariance matrix described in Section 3.1. Second, we analyze the performance with synthetic spectra. Finally, we apply the algorithm to the APEX dataset to identify potential variations of spectral calibration in the across- and along-track direction.

4.1. Degrees of Freedom

How well the state vector can be retrieved depends critically on the specified measurement and a priori covariance matrices. For the a priori state vector and covariance matrix defined in the previous section, averaging kernel matrices were calculated for both the case of a small measurement error of 0.1 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ and the case of a large error of 0.5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$. This range is based on the error estimates presented in Section 3.2.2 (Table 3). The state vector has 165 elements, and the total degrees of freedom for these errors are 139 and 112, respectively.

The degrees of freedom for trace gases are 0.68 and 0.36 for ${\mathrm{O}}_2$, 1.00 and 0.97 for ${\mathrm{H}}_2\mathrm{O}$, 0.06 and 0.00 for ${\mathrm{O}}_3$, 0.08 and 0.03 for ${\mathrm{O}}_4$ and 0.34 and 0.03 for ${\mathrm{NO}}_2$ for small and large errors, respectively. ${\mathrm{O}}_2$ and ${\mathrm{H}}_2\mathrm{O}$ can thus be retrieved well by the algorithm. The smaller degree of freedom for ${\mathrm{O}}_2$ is due to the a priori variance being very small. The other trace gases are not retrieved and mainly remain at their a priori values. An exception is ${\mathrm{NO}}_2$ for low noise, which has a relative high degree of freedom ($d_s$ = 0.34).

For the splines, the degrees of freedom depend on wavelength as shown in Figure 5. The albedo is mainly constrained by the absolute radiance level outside of the absorption features and can be successfully fitted across the whole spectral range. The other parameters are more difficult to fit. The retrieval works well below 500 nm, where several strong Fraunhofer lines provide useful information to constrain wavelength shift, slit function and radiance offset, and around the ${\mathrm{O}}_2$ A- and B-bands at 762 nm and 687 nm, respectively. Above 600 nm, Fraunhofer lines are less important, because they are less deep and the FWHM of the instrument larger. In between, the information content is lower such that wavelength shift and the slit function parameter can only be retrieved poorly, and the offset cannot be retrieved at all.

4.2. Smoothing Errors

The averaging kernel matrix $\mathbf{A}$ can also be used to analyze the smoothing errors of the retrieval. Figure 6 shows matrix $\mathbf{A}$ for a large measurement error (${\sigma }_{\epsilon }$ = 0.5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$) normalized by the a priori standard deviations. Each row of $\mathbf{A}$ can be regarded as the smoothing function where all non-zero elements in the row contribute to the retrieved state vector element (Equation (21)). Since the state vector is mainly composed of control points for the different splines, $\mathbf{A}$ can be divided into individual block matrices each representing a single spline. Figure 6 shows this division as vertical and horizontal lines. For an ideal retrieval method, $\mathbf{A}$ would be the identity matrix, and the smoothing error would be zero. In our case, the smoothing function for a control point of a spline is a more or less sharply peaked function with a maximum on the respective node. The width of the peak is related to the correlation length specified for the a priori information. Figure 6 shows that these functions have a sharp peak for wavelength shift and slit function. The radiance offset, in contrast, is smeared out over the neighboring nodes due to the large correlation length of the a priori covariance matrix. In addition, other splines can contribute to the control point if the corresponding parameters cannot be retrieved independently.

In the following, we estimate the smoothing error according to Equation (21) by creating a true state vector and adding or subtracting the a priori standard deviation from the a priori state vector (Table 4). We do this for each spline independently to identify any potential interactions between the different splines.

The wavelength shift has only a small smoothing error. Adding a shift of 0.2 nm results in an underestimation by less than 1.2% of spectral resolution increasing from 0.02 nm below 500 nm to 0.08 nm at 800 nm. The influence of the other splines and trace gases on the wavelength shift is negligible.

The smoothing error for the FWHM spline is also small. If the FWHM is increased by 15%, the retrieved FWHM is underestimated by less than 3%. A radiance offset of 5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ only increases the FWHM by less than 0.4%. A wavelength shift (±0.2 nm) has no significant effect on the FWHM. If ${\mathrm{O}}_2$ is increased by 3%, the FWHMs are underestimated by less than 0.5%. Similarly, if ${\mathrm{H}}_2\mathrm{O}$ is increased by 30%, the FWHMs are underestimated by less than 0.8%. The effect of other trace gases together is less than 0.3%.

The smoothing error of the radiance offset, in contrast, is large, which can already be seen in Figure 4c. Increasing the radiance offset by 5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ results in an underestimation of the retrieved offset between 0.25 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ at 400 nm and −1.0 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ at 750 nm. Increasing the FWHM has the opposite effect: a 15% larger FWHM results in an increase between 1.5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ at 400 nm and 2.5 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ at 750 nm. Changing the wavelength shift by +0.2 nm has no strong effect on the radiance offset (<0.1 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$). Increasing ${\mathrm{O}}_2$ by 3% or ${\mathrm{H}}_2\mathrm{O}$ by 30% increases the radiance offset, but the effect is small. The other trace gases increased simultaneously by 30% result in a less than 0.2 mW·nm${}^{-1}$·m${}^{-2}$·sr${}^{-1}$ increase of the offset.

In summary, smoothing errors are small for shift and FWHM, but important for the radiance offset. FWHMs and offset are positively correlated, but while the correlation causes a large smoothing error in the radiance offset, the error on the FWHMs is small.

4.3. Performance with Synthetic Spectra

The performance of the spectral calibration algorithm was analyzed with synthetic spectra constructed for a set of scenarios differing from the a priori assumptions used in the algorithm (see Section 2.1). Instrument errors were added to the synthetic spectra as random noise for small and large errors of 0.1 and 0.5 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$, respectively. To account for parametrization and forward model errors, an additional model error of 0.1 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ was considered for constructing the measurement covariance matrix.

Figure 7 presents the retrieval results for an ensemble of 20 spectra over “green grass” with a wavelength shift due to an instrument pressure change by 15 hPa and a FWHM scaling of 1.2. Figure 7a shows that the forward model fits the synthetic spectra on average extremely well and does not show any systematic errors due to, for example, the albedo parametrization. Note that the synthetic spectrum is hidden behind the averaged forward model spectrum (red line). Figure 7b–f shows fitted wavelength shift, FWHM and radiance offset. The blue line represents the true parameter.

Table 5. Estimated errors (Equation (19)) of retrieved wavelength shifts, FWHMs and radiance offsets for synthetic spectra with different surface reflectances. The wavelength shifts corresponds to a 15-hPa pressure change; the laboratory FWHMs were scaled by 1.2; and the offset term (Equation (1)) was added to the spectra. Measurement errors were 0.1 and 0.5 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$, and an additional error of 0.1 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ was considered in the calculation of the total retrieval error to account for forward model errors.

		Small Error: ${\mathit{\sigma }}_{\mathit{\epsilon }}$ = 0.1 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$			Large Error: ${\mathit{\sigma }}_{\mathit{\epsilon }}$ = 0.5 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$
	λ(nm)	Shift (nm)	FWHM (%)	Offset ($\frac{\mathbf{mW}}{{\mathbf{m}}^{\mathbf{2}}·\mathbf{nm}·\mathbf{sr}}$)	Shift (nm)	FWHM (%)	Offset ($\frac{\mathbf{mW}}{{\mathbf{m}}^{\mathbf{2}}·\mathbf{nm}·\mathbf{sr}}$)
	400	0.01	3.0	0.7	0.03	5.1	1.0
	450	0.03	4.3	0.9	0.06	6.8	1.4
(a)	500	0.06	4.9	1.3	0.08	7.2	1.7
	600	0.11	7.3	1.4	0.14	9.7	1.9
	760	0.08	4.3	0.9	0.13	7.1	1.7
	400	0.02	3.4	0.6	0.03	5.6	1.0
	450	0.04	5.0	0.9	0.06	7.4	1.4
(b)	500	0.06	5.6	1.3	0.09	7.9	1.8
	600	0.10	7.1	1.6	0.13	9.5	2.0
	760	0.02	1.3	1.5	0.05	2.4	2.1
	400	0.00	1.2	0.7	0.01	2.7	1.1
	450	0.01	1.4	1.0	0.02	3.3	1.5
(c)	500	0.02	1.8	1.3	0.04	3.7	1.9
	600	0.05	3.2	1.6	0.08	5.3	2.1
	760	0.02	1.0	1.5	0.05	2.2	2.2
surface reflectance: (a) “asphalt paving”, (b) “green grass” and (c) “metal roofing”

summarizes the retrieval errors according to Equation (19) for this case, as well as for “asphalt paving” and “metal roofing”. Since the true shift and FWHMs were known for these tests, we also calculated the retrieval errors from the ensemble. These errors agree well with the analytic errors calculated by the algorithm.

The wavelength shift is retrieved with good accuracy. The uncertainty of the shift is between 0.03 and 0.09 nm over a dark surface in the visible. In the case of large measurement errors, shifts above 700 nm are underestimated as expected from the smoothing error (Figure 7d). Adding the ring spectrum to the true spectrum has only a small impact on the retrieved shift (<0.004 nm).

The FWHM of the slit function is fitted with good accuracy even for the case with high measurement uncertainty. The a posteriori uncertainty (2%–10%) significantly reduced with respect to the a priori uncertainty (15%). For an extreme case where the true FWHM is twice the laboratory calibration, the retrieval would mostly underestimate the FWHM except near strong Fraunhofer lines and the ${\mathrm{O}}_2$ A-band. Furthermore, this would introduce substantial errors in the determination of the wavelength shift (<0.15 nm), because Fraunhofer lines and absorber lines are smoothed out by the wide slit function. Weak atmospheric absorbers, such as ${\mathrm{O}}_3$ and ${\mathrm{O}}_4$, have no impact on the FWHM retrieval. For example, setting these absorbers to zero or multiplying the a priori by five has no impact on the retrieved FWHM. Setting the radiance offset to zero, i.e., not fitting the offset, has a small effect on the FWHM (<0.4%). Adding a Ring spectrum to the measurement spectrum results only in a small misfit (<1% below 450 nm).

The radiance offset can only be retrieved well if the measurement errors are small. For large errors, the radiance offset has a large smoothing error. Ignoring the Ring effect overestimates the radiance offset by 1.0 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ at 400 nm, 0.5 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ at 500 nm and even smaller for larger wavelengths.

4.4. Performance with APEX Measurements

The spectral calibration algorithm was applied independently to each spectrum in the dataset to account for potential variations in both the along-track and across-track direction. Figure 8, Figure 9 and Figure 10 show the main results for the spectral calibration of the test dataset. Figure 8 presents the spectra, radiance offset, wavelength shift and FWHM for an individual line. Figure 9 shows the across- and along-track variability of shift, FWHM and radiance offset for different wavelength ranges, and Figure 10 shows the across- and along-track means over the full spectral range.

4.4.1. Spectral Residuals

The fraction of systematic error in the instrument error $f_s$ (Equation (24)) was estimated by applying the spectral calibration algorithm to spectra averaged in the along-track direction. These averaged spectra consist of about 15,000 individual spectra and, thus, have a negligible random component (<0.01 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$). The systematic error was estimated from the fitting residual to 0.40 ± 0.04 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ ($f_s$ = 31% of total error; see Section 3.2.2). This includes systematic errors due to not accounting for the Ring effect and due to parametrization errors (see Table 3), but is likely dominated by radiometric calibration errors of the Level-1 processor. For the test dataset, N is 339 ± 176 pixels, and thus, the average ${\sigma }_{\epsilon }$ is 0.45 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$.

For all spectra in the APEX dataset, the mean and standard deviation of the root mean square (RMS) of the residuals, i.e., the difference between the calculated and measured radiance spectrum, are 0.44 ± 0.07 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$. This RMS suggests a reasonable quality of the fit without “overfitting” and a realistic error estimation, because the RMS is very similar to the measurement error calculated by Equation (24). Figure 8a shows the averaged measured and fitted spectra, as well as the averaged residual of 50 spectra. The residual does not vanish by averaging many fits, suggesting a significant systematic error. Therefore, the estimate of the fitted state vector is also not improved by averaging. The RMS of the residuals changes with wavelength: it is of the order of 0.76 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ at 400 nm and 0.25 mW·m${}^{-2}$·nm${}^{-1}$·sr${}^{-1}$ at 500 and 760 nm, respectively.

4.4.2. Wavelength Shift

A wavelength-dependent spectral shift was observed (Figure 8c for all spectra in Line 10). The a priori shift is shown with cyan circles and accounts for the laboratory assessment of spectral smile. The black curve shows the averaged retrieved shifts and the median errors. The median errors are similar to the errors calculated for the synthetic spectra with large noise (Table 3).

Retrieved wavelength shifts change in the across- and along-track direction. Figure 9a,b shows the wavelength shift in the visible (410–520 nm) and near infra-red (720–815 nm). The shift does not depend strongly on surface type: we found only slightly higher values above forest, suggesting that errors due to the albedo parametrizations are negligible (compare Figure 9a with Figure 3). Near along-track Position 800, several pixels are missing or have a large error due to the issue of detector saturation over farmland surfaces mentioned earlier. The maps show the expected across-track smile and a previously unknown along-track dependency.

Figure 10a shows the spectral smile normalized by the laboratory calibrated FWHM for different wavelength bands. The spectral smile measured in the laboratory does not depend strongly on wavelength. The shapes of in-flight and laboratory measurements of the spectral smile agree well, but the in-flight values are shifted upwards, because the pressure sensitivity of the instrument affects different wavelengths differently [14].

Figure 10b shows the same, but in the along-track direction. The along-track drift is smaller, but is not negligible. For all wavelengths, the shift increases slowly in the along-track direction at a rate of 0.06 ± 0.02 normalized spectral pixels per 1000 along-track pixels. This corresponds to −30 ± 9 hPa per 1000 pixels, which is about twice the accuracy of the sensor regulating the instrument pressure. The housekeeping data suggest that the pressure was constant during the measurements.

4.4.3. Instrument Slit Function (FWHM) and Radiance Offset

Since FWHMs and radiance offsets depend on each other due to the smoothing error (cf. Section 4.2), we are analyzing the two parameters together in this section.

The wavelength-dependent FWHMs are similar to the laboratory calibration, but up to 50% larger, especially between 450 and 600 nm (Figure 8d). They also depend on surface type and are larger over bright surfaces (compare Figure 3 with Figure 9c for 410–520 nm). This behavior is unexpected and unrealistic. The retrieved radiance offsets do not clearly depend on the surface type (Figure 9d), but depend slightly on absolute levels of the measured radiance (Figure 11a). In the across-track direction, the radiance offset has a pronounced drop-off towards the edges of the detector (Figure 10d), although the measured radiance does not change strongly in the across-track direction. The wavelength-dependent radiance offset is retrieved with high accuracy below 450 nm, while it is mainly constant above 500 nm (Figure 8b).

Since FWHMs seem to depend on surface type, which is unlikely, we repeated the retrieval with a stronger constraint on the a priori standard deviation of 1% FWHM${}_{\mathrm{lab}}$. In this case, the retrieved FWHM remains close to the laboratory calibration while the radiance offset increases significantly and now depends on absolute levels of the measured radiance (Figure 11b) and surface type (Figure 11c). The radiance offsets are about 20% of the radiance signal.

5. Discussion

5.1. Usefulness of the Algorithm

A new algorithm has been presented for in-flight spectral calibration of spectrometers with wide spectral windows and moderately high spectral resolution. The algorithm has been carefully designed to be generally applicable to any spectrometer potentially benefiting from in-flight spectral calibration (see the Supplement for details) and has been exemplarily adapted to the APEX imaging spectrometer.

In-flight spectral calibration of airborne and satellite instruments is essential to identify deviations from the reference calibration obtained from on-ground or IFC facility measurements. Our algorithm significantly improved the accuracy of APEX’s spectral calibration and enabled a detailed characterization of the instrument’s across-tack spectral smile and a previously unknown along-track drift. This accurate spectral calibration is necessary for the retrievals of atmospheric and surface variables that rely on narrow spectral features, such as the retrievals of atmospheric trace gases [3,6] and chlorophyll fluorescence [5].

The presented algorithm is similar to algorithms based on DOAS that are frequently used for spectral calibration in the ultraviolet, visible and infra-red spectral range [6]. However, DOAS-based approaches are limited to small spectral windows, because they use low-order polynomials to approximate the wavelength dependency of albedo, wavelength shifts and variations of slit function parameters across the wavelength window. Our algorithm overcomes this limitation by using splines instead of low-order polynomials, allowing it to be applied to spectral windows of any width.

Laboratory and in-flight measurements of instrument behavior can be easily incorporated as a priori information in the MAP optical estimate. The a priori information can improve the quality of the fit and allows extending the calibration to spectral regions where Fraunhofer lines and atmospheric absorbers do not provide sufficient information for a successful fit (see Figure 5).

5.2. Potential Improvements of the Algorithm

The accuracy of the in-flight calibration of CWs and FWHMs is, in general, limited by instrument and forward model errors and is also determined by the spectral resolution of the instrument. In this study, we did not investigate the influence of spectral resolution, but we analyzed in detail the effects of instrument and forward model errors. For the APEX instrument, the increased RMS below 500 nm is likely caused in roughly equal parts by systematic errors in the radiometric calibration and the missing Ring effect in the forward model, while for larger wavelengths, the forward model error is smaller than the instrument error (see Table 3).

We found that systematic errors are critically limiting the quality of the retrieval of CWs and FWHMs. An apparently simple approach to overcome this issue would be to fit the mean residual as an additional term in the forward model. However, this approach is justified on the condition that the mean residual and the state vector do not influence each other. Unfortunately, this condition is not fulfilled in our case where the systematic structure of the residuals strongly depends on the spectral calibration described by the state vector. The mean residual can therefore not be added to the forward model. The systematic error can be reduced only by improving the radiometric calibration of the instrument.

The forward model was designed to be computationally efficient. It only requires a small set of inputs (solar irradiance spectrum, trace gas absorption cross-sections and optical depth profiles) and approximates parameters that vary smoothly with wavelength by cubic splines. The model is sufficiently accurate for the spectral resolution and noise level of an instrument similar to APEX, but it may have to be improved for instruments with higher resolution or better signal to noise performance. For example, the largest error in the forward model is caused by the missing Ring effect that can be included by fitting a Ring spectrum in Equation (8) as a pseudo absorber [40]. Additional potential improvements are: for strong atmospheric scattering, e.g., high aerosol load over dark surface, the geometric AMFs can be replaced by AMFs pre-computed for the a priori state vector using a radiative transfer model. To account for wavelength-dependent SCDs, wavelength-dependent AMFs can be included in the forward model.

5.3. Challenges and Issues Related to the Analysis of Real Data

Wavelength shifts could be retrieved for the APEX dataset with high accuracy with an uncertainty smaller than 0.1 nm (i.e., 0.05 spectral pixels). This is smaller than the uncertainty of the laboratory calibration of 0.2 nm. APEX’s across-track spectral smile could therefore be captured well by the algorithm and agreed reasonable well with the smile expected from the laboratory calibration. In addition, we identified an unexpected along-track drift of the spectral shift, which was attributed to a change in pressure within the spectrometer due to the limited accuracy of the pressure regulation.

Based on the analysis with synthetic spectra, FWHMs should be retrievable with high accuracy from the APEX spectra. However, FWHM and radiance offset turned out to be more difficult to retrieve, because both are very sensitive to the filling-in of Fraunhofer lines, which are a primary source of information for quantifying these parameters. Both the retrieved FWHMs and radiance offsets are larger than expected from the instrument characterization in the laboratory. The larger FWHMs were also reported by Popp et al. [3], who were not able to provide an explanation for this observation. The disagreement is likely too large to be solely explainable by a difference between laboratory and in-flight performance. In particular, it is unreasonable that FWHMs are higher over brighter scenes. A close look at the APEX spectra shows that the misfits are the result of Fraunhofer lines near 500 nm that are not as deep as expected from smaller FWHMs. This might be caused by filling-in of the lines by the Ring effect or an underestimation of ${\mathrm{O}}_3$ and ${\mathrm{O}}_4$ absorption. However, our tests with synthetic and real spectra show that all of these effects are too small to explain the misfit, making forward model errors an unlikely reason.

Spatial binning of radiance spectra with slightly varying CWs also causes a filling-in of sharp lines. The effect is largest in the across-track direction because of the spectral smile. We estimated the magnitude of the filling-in through spatial binning with synthetic spectra with a realistic smile and found that the FWHM is overestimated by less than 0.2% for a binning of forty spectra in the across-track direction. The effect is thus too small to explain the misfit in FWHM and consequently can be neglected here.

Radiometric calibration errors may also result in a similar filling-in of dark lines. The dichroic coating on the prisms can impact the radiometric calibration and is known to cause systematic errors in APEX [37]. However, the coating is not expected to absorb near 500 nm. Another possible explanation could be a not-fully-corrected radiance offset due to stray light or vignetting. A not-fully-corrected vignetting can cause a drop-off of radiance offsets towards the detector edge, as seen in Figure 10d. However, since the radiance offset would have to be has high as 20% of the radiance signal, i.e., much larger than the instrument error, vignetting and stray light are unlikely to be the sole explanation. Another effect is charge-coupled device (CCD) readout smear. It occurs when the radiance signal is read out by shifting charges between wells in the spectral direction without a shutter. This effect acts like a wider slit function and results in a spectral smoothing. The effect is corrected by a linear model in the Level 1 Processor, but our analysis suggests that the readout smear is non-linear and, thus, not fully corrected. Therefore, a combination of not-fully-corrected CCD readout smear, vignetting and stray light issues is likely responsible for the unreasonably high FWHMs and radiance offsets retrieved by the algorithm. This emphasizes the need for accurate radiometric calibration and characterization of stray light effects.

6. Conclusions

Imaging spectroscopy data products are increasingly used for advanced Earth science applications. Accurate spectral calibration is essential for the retrieval of various atmospheric and surface variables that rely on sharp spectral features, such as atmospheric trace gases and chlorophyll fluorescence. Since small changes of spectral calibration can occur during flights, e.g., due to temperature and pressure changes in the instrument, accurate in-flight calibration is essential.

We presented a new spectral calibration algorithm that retrieves continuous wavelength shifts and slit function parameters by fitting a high-resolution solar spectrum and atmospheric absorbers to in-flight spectra. The algorithm has three key features: (1) the forward model is designed to be computationally efficient and can be adapted easily to different instruments; (2) the use of splines allows conducting spectral calibrations across a wide spectral window with a single fit; (3) laboratory and in-flight measurements of instrument behavior can easily be incorporated as a priori knowledge in the maximum a posteriori optimal estimate.

The algorithm was comprehensively tested with synthetic spectra and applied to spectra from the APEX imaging spectrometer. Wavelength shifts were retrieved with high accuracy. APEX’s across-track spectral smile and a previously unknown along-track drift were captured well. The FWHM of the slit function could also be retrieved with high accuracy for synthetic spectra, while some misfits appear for the APEX spectra that are likely related to radiometric calibration issues, which emphasizes the necessity for comprehensive calibration and performance measurements of high resolution imaging spectrometers.

In conclusion, the algorithm can significantly improve the in-flight spectral calibration of APEX and similar instruments and can also help improve the retrievals of atmospheric and surface variables relying on sharp spectral features.