Level 0 to 1 processing of the imaging Fourier transform spectrometer GLORIA : generation of radiometrically and spectrally calibrated spectra

The Gimballed Limb Observer for Radiance Imaging of the Atmosphere (GLORIA) is an imaging Fourier transform spectrometer that is capable of operating on various high-altitude research aircraft. It measures the atmospheric emission in the thermal infrared spectral region in limb and nadir geometry. GLORIA consists of a classical Michelson interferometer combined with an infrared camera. The infrared detector has a usable area of 128× 128 pixels, measuring up to 16 384 interferograms simultaneously. Imaging Fourier transform spectrometers impose a number of challenges with respect to instrument calibration and algorithm development. The optical setup with extremely high optical throughput requires the development of new methods and algorithms for spectral and radiometric calibration. Due to the vast amount of data there is a high demand for scientifically intelligent optimisation of the data processing. This paper outlines the characterisation and processing steps required for the generation of radiometrically and spectrally calibrated spectra. Methods for performance optimisation of the processing algorithm are presented. The performance of the data processing and the quality of the calibrated spectra are demonstrated for measurements collected during the first deployments of GLORIA on aircraft.


Introduction
The upper troposphere/lower stratosphere (UTLS) is a region of particular importance for radiative forcing.Especially uncertainties of exchange processes between the tropical upper troposphere and the lowermost stratosphere are a major error source in the Earth radiation budget (Riese et al., 2012).In a more general sense, transport pathways of air between different compartments of the atmosphere affecting dynamics and chemistry in a changing climate need to be studied more thoroughly.In order to fully understand these processes, measurements are required which provide both high spatial resolution and sufficient coverage (Riese et al., 2014).
The Gimballed Limb Observer for Radiance Imaging of the Atmosphere (GLORIA) is an airborne imaging Fourier transform spectrometer (FTS) representing the first realisation of the infrared-emission limb-imaging technique (Riese et al., 2005;Friedl-Vallon et al., 2006).It is designed to measure two-and three-dimensional trace gas distributions in the UTLS region with high spatial resolution (Friedl-Vallon et al., 2014).For this purpose GLORIA combines a classical Michelson interferometer with a detector array providing a good spatial coverage and resolution together with a good spectral resolution and a sufficient signal-to-noise ratio.It is deployed in the belly pod of the German research aircraft Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Kleinert et al.: Level 0 to 1 processing of GLORIA HALO as well as on board the Russian high-altitude research plane M55 Geophysica.
The transition from a single detector FTS to a 2-D imaging FTS requires new methods in data evaluation and exploitation.The optical throughput of such a system is much higher than the throughput of a conventional FTS, allowing larger spatial coverage and higher sensitivity with a smaller instrument.As a challenge, such an instrument setup leads to different spectrometric properties for each pixel.In particular, the spectral axis is different for each pixel and new methods for an accurate spectral calibration are required in order to allow for a coherent scientific interpretation of different pixels.
In principle, each detector element together with the interferometer is an independent spectrometer.Properties like sensitivity, linearity and spectral cut-off of these several thousand independent spectrometers may vary randomly from pixel to pixel and/or systematically over the detector array.These properties have to be characterised and considered during calibration.Furthermore, an imaging FTS delivers much more data than a conventional spectrometer.This requires intelligent and efficient strategies for data processing.
This paper describes the level 0 and level 1 processing for GLORIA.Level 0 processing consists of resampling the raw interferograms on a space equidistant grid and the level 1 processing generates calibrated spectral radiances from the level 0 product.After a short instrument description (Sect.2) and an overview of the processing steps (Sect.3), the approaches for radiometric and spectral calibration are presented (Sects.4 and 5).The level 0 and level 1 processing steps are described in detail in Sect.6. Section 7 deals with the performance optimisation of the processor.The characterisation and processing results for scientific flights performed in 2012 are finally presented in Sect.8.

GLORIA instrument and data acquisition
The GLORIA spectrometer is a cooled imaging FTS with a large cryogenic HgCdTe detector array for the detection of infrared radiation in the spectral range of 780 to 1400 cm −1 .A detailed description of the instrument is given by Friedl-Vallon et al. (2014).
The heart of the FTS is a classical Michelson interferometer with a maximum optical path difference (MOPD) of ±9 cm.The interferogram length can be chosen freely within this range.
The infrared detector is an HgCdTe large focal plane array (LFPA) consisting of 256 × 256 pixels with a pitch of 40 µm and a stare-while-scan readout.Eight channels of 14 bit ADCs, synchronously clocked with 10 MHz on a separate front-end electronics, perform the digitalisation of the LFPA output.
For the GLORIA measurements, only a subset of the array is actually used.The optics is designed to cover a square of 128 × 128 pixels, but the typical flight configuration in the current setup uses only 48 pixels horizontally and 128 pixels vertically for an optimised interferogram sampling frequency.With this configuration, the sampling frequency or frame rate is 6281 Hz resulting in a detector raw data rate of 73.6 MiB s −1 (1 MiB = 1 mebibyte = 2 20 bytes) by using a 16-bit-wide integer value representation.
Each frame is marked with a unique time stamp by the interferometer electronics, and all frames of an interferogram are gathered into so-called data cuboids.The cuboids are transferred to the central computer where they are written onto a high-speed RAID-system.Figure 2 shows a schematic drawing of such a cuboid.
The optical path difference of the interferometer is measured with a laser reference system.A diode laser signal with a wavelength of about 646 nm is coupled into the interferometer, and the rising zero crossings of the laser interferogram are detected and marked with a time stamp.These time stamps are stored as laser data files.The time stamps for both the infrared and the laser signal are generated by a 80 MHz clock within the interferometer electronics.They are given in integer time ticks with one tick corresponding to 12.5 ns.
Each measurement consists of one cuboid containing the infrared signal sampled with constant frequency, and a laser file containing the time stamps of the rising zero crossings of the laser reference system.These two pieces of information make it possible to reconstruct the interferograms on an equidistant grid in space as a function of the optical path difference (see Sect. 6.1).
The spectrometer is mounted in a gimballed frame that allows agility in azimuthal, elevational and image rotation direction.The instrument is designed to measure in limb-or nadir-viewing geometry.In limb mode, the azimuthal viewing direction can be chosen between 45 and 132 • to the flight direction.In addition, the centre of the FOV can be directed to an elevation of about +10 • for the so-called deep space view and into two large area blackbodies for radiometric calibration.
Two principal measurement modes have been used for limb measurements: -The chemistry mode with focus on the atmospheric composition has an MOPD of ±8 cm and a spectral sampling of 0.0625 cm −1 .The azimuthal viewing direction is perpendicular to the flight direction.
-The dynamics mode with focus on small-scale dynamics of the atmosphere has an MOPD of ±0.8 cm and a spectral sampling of 0.625 cm −1 .The azimuth is scanning from 128 • to 44 • in 4 • steps measuring the same air mass from different viewing angles.
The optical interferometer velocity is set to 1.27 cm s −1 , leading to an interferogram acquisition time of about 12 s in chemistry mode and 1.2 s in dynamics mode, plus a turnaround time of 0.8 s in both modes.

Processing overview
The level 0 to 1 processing transforms the measured raw interferograms into radiometrically and spectrally calibrated spectra.The principal processing scheme is depicted in Fig. 1.Details of the individual processing steps are described in Sect.6.First of all, the cuboids are checked for missing frames and for spikes.In a next step, the raw interferograms are corrected for the nonlinearity of the detector system (see Sect. 4.2).This is a prerequisite for the radiometric calibration.Furthermore, all interferograms have to be aligned to the same ZOPD (see Sect. 6.1.5).
The core of the level 0 processing is the resampling of the interferograms on a space-equidistant grid using the information from the laser file.The resampling comprises the spectral calibration by taking the correct laser wavelength and the off-axis angle for each pixel into account.Laser wavelength and off-axis angle are determined from in-flight measurements.The determination is described in detail in Sect. 5.
The output of the level 0 processor consists of data cuboids containing nonlinearity corrected, spectrally calibrated interferograms that are sampled on a space-equidistant grid.
The level 1 processing comprises the Fourier transform converting the interferograms to complex spectra and the radiometric calibration.For the radiometric calibration, the atmospheric spectra are divided by the gain function and the offset is subtracted.Gain and offset spectra are determined from dedicated calibration measurements (see Sect. 4).The radiometric calibration is performed in the complex domain following the approach proposed by Revercomb et al. (1988).
The level 1 product consists of complex calibrated spectra.When the calibration has been performed correctly, all atmospheric contribution is contained in the real part of the spectrum while the imaginary part contains only noise.The imaginary part serves as quality control and is used to determine the noise equivalent spectral radiance (NESR).

Calibration approach
The radiometric calibration assigns absolute radiance units to the arbitrary intensity units of the measured spectra and the instrument self-emission contributing to the measured spectra is determined and subtracted from the signal.For a sound trace gas retrieval, the goal requirement for the radiometric gain accuracy is 1 % with a threshold of 2 % (Friedl-Vallon et al., 2014).This requirement is based on experiences with the MIPAS and CRISTA instruments and on studies performed for a satellite-borne imaging FTS (ESA, 2012).
The standard approach for radiometric calibration is to look at two blackbodies with different temperatures.Assuming a linear system, gain g and offset o can be derived from these two measurements (see e.g.Revercomb et al., 1988): S h and S c are the measured spectra of the hot and cold blackbody, respectively, and B h and B c are the corresponding radiance spectra which are calculated from the blackbody temperatures following Planck's law.The calibrated spectrum L atm is then calculated from the measured spectrum S atm using The gain is determined by the modulation efficiency and transmittance of the instrument and the responsivity of the detector element.The offset is dominated by the modulated part of the self-emission of the instrument.
The two-point calibration approach is only valid if the detected signal is linear with respect to the incoming radiance.Any nonlinearity of the detector system must be corrected before the calibration can be applied.The nonlinearity of the detection chain of GLORIA is addressed in Sect.4.2.
For the radiometric calibration of the data measured by GLORIA, two large-area blackbody radiation sources (126 mm × 126 mm) are flown on the aircraft.The blackbodies are mounted outside of the spectrometer placing all optical elements including the spectrometer entrance window within the light path during the calibration measurements.The blackbodies are temperature stabilised using thermoelectric coolers which allow likewise heating and cooling.Basic requirements are an emissivity of larger than 0.997 and a temperature homogeneity and knowledge of better than 0.1 K in order to achieve the radiometric accuracy goal for the calibrated spectra.Details of the blackbody design are given by Olschewski et al. (2013).
The blackbody measurements are complemented by socalled deep space measurements, where the instrument is looking upwards into space with an elevation angle of +10 • .Deep space is a commonly used cold calibration source for satellite measurements, because it directly gives the instrument self-emission.This is not entirely true for airborne measurements because of residual emission from trace gases above flight altitude.These residual atmospheric signatures have to be removed from the measured deep space signal in order to get the instrument self-emission.A method for the removal of atmospheric lines from deep space measurements is presented in Sect.4.3.So in fact three calibration points are available for GLORIA measurements.
A typical calibration sequence in flight consists of 20 hot and 20 cold blackbody measurements, followed by 10 deep space measurements.Blackbody measurements are performed with the dynamics mode spectral resolution, deep space measurements with the chemistry mode spectral resolution.For an optimised signal-to-noise ratio, the detector integration time is adjusted to the intensity of the source, i.e. the integration time for blackbody measurements is shorter than the integration time for atmospheric and deep space measurements.
In order to account for the thermal drift of the instrument, calibration measurements are performed about every 30 to 45 min.Gain and offset are interpolated linearly in time to the time of each atmospheric measurement.With a duration of about 5 min per calibration sequence, about 10 to 15 % of the available measurement time is currently used for calibration.
In principle, only two of the three calibration sources are needed to calculate gain and offset.The combination of any two of the three sources should give the same results.In reality, differences can be found due to 1. a non-perfect nonlinearity correction 2. errors in the temperature and emissivity of the blackbodies 3. a non-perfect removal of atmospheric contributions to the deep space measurements 4. measurement noise.
During the flights performed so far, the ambient temperature in the belly pod of the aircraft was typically 245 to 260 K and thus considerably above the outside air temperature, limiting the cooling of the cold blackbody.This led to a relatively small temperature difference of typically only 15 to 25 K instead of the target value of 40 K.The small temperature difference and thus the small difference in the measured blackbody spectra (S h −S c ) makes the blackbody-blackbody calibration very susceptible to measurement noise, and the extrapolation to zero input, which is needed to determine the instrument offset, is critical.Therefore the current baseline is to use primarily the cold blackbody and the deep space measurements for calibration.The additional data of the hot blackbody can be used to check the data for consistency and to enhance the data quality especially in the spectral regions where atmospheric contributions to the deep space measurements are strong and the determination of the instrument baseline is difficult.Improvements in the thermal design and the operation procedures are under development in order to enhance the quality of the blackbody measurements for future flights.
Using a blackbody-deep space calibration, Eqs. ( 1) and ( 2) are simplified to where bb and ds denote the (cold) blackbody and deep space measurements, respectively.

Nonlinearity correction
The detector system exhibits a certain nonlinearity, i.e. the output signal is not directly proportional to the input signal.It is assumed that the photovoltaic detector itself is basically linear, and the nonlinearity mainly stems from the readout electronics.The nonlinearity is characterised by varying the integration time and therewith the number of recorded electrons while looking at a constant radiation source (Hilbert, 2004;P. Giaccari, personal communication, 2010).The recorded DC values are mapped on a virtual detector system, which is linear, i.e. the output signal is directly proportional to the integration time.For the definition of the virtual detector, a linear fit was made to three measurement points from the middle region of the detector's dynamical range.Figure 3 shows the recorded signal as a function of integration time for a single pixel (black) together with the signal of the virtual linear detector system (red line).The relation between the measured and the linear DC values is expressed by a fourth-order polynomial fit.This polynomial is used to create a lookup table which assigns to each measured value of the interferogram the value that would have been recorded by the virtual linear detector.The characteristic of the nonlinearity is sufficiently similar for all nominally working pixels, such that a single lookup table can be applied to all pixels (Sha, 2013).
The quality of the nonlinearity correction is cross-checked by quantifying the out-of-band artefacts in measured blackbody spectra before and after nonlinearity correction (Kleinert, 2006).The artefacts arise from the distortion of the interferogram peak due to the nonlinearity.In first approximation, the nonlinearity can be described by a second-order polynomial: In the spectral domain, the quadratic term corresponds to a convolution of the spectrum with itself.If the original signal is limited to a spectral range [σ 1 , σ 2 ], the quadratic term leads to out-of-band artefacts in the spectral ranges [0, σ 2 −σ 1 ] and [2σ 1 , 2σ 2 ].Spectral artefacts from higher-order nonlinearity contributions are below noise level and not considered for quality control.Figure 4 shows the quadratic artefact before and after nonlinearity correction for a selected pixel for a blackbody measurement using two integration times.After nonlinearity correction, the artefact is clearly reduced, showing that the characterisation measurements that were made with variable integration time but constant photon flux can be transferred to in-flight data where the integration time is constant during one interferogram acquisition but the photon flux is variable.This has been shown to be valid for different integration times.
Figure 5 shows the intensity of the magnitude spectra in the spectral range of 20 to 400 cm −1 before and after nonlinearity correction.It can be seen that the artefact is considerably reduced by the nonlinearity correction for most of  the pixels.The residual value is mainly due to noise which gives a net positive contribution in the magnitude spectrum.
The pixels where the nonlinearity correction did not work have to be either discarded from processing or have to be treated separately.For the discarding of pixels, see Sect.8.3.

Removal of atmospheric contributions from deep space measurements
The determination of the gain and offset functions on basis of deep space and blackbody observations is complicated by the presence of atmospheric signatures in the deep space spectra.This atmospheric contamination depends on the atmospheric situation, the flight altitude and on the maximum upward viewing elevation angle.This maximum pointing elevation (10 • ) is limited by obstruction of the field of view for larger angles due to the mounting of the instrument in the belly pod of the aircraft.A typical deep space raw spectrum is shown as the black curve in Fig. 6 where various trace gas signatures are visible.For removal of these residual atmospheric signatures, a scheme first used for calibration of MIPAS-STR (Michelson Interferometer for Passive Atmospheric Sounding -STRatospheric aircraft) observations (Höpfner et al., 2001;Woiwode et al., 2012) was applied.This method relies on radiative transfer forward calculations which are applied to simulate as well as possible the actual atmospheric state using ECMWF analysis for temperature and water vapour together with climatological profiles of trace gases (Remedios et al., 2007).In an iterative process (i = 1, . .., i max ), the observed deep space spectrum is corrected as follows: S ds,i = S ds,0 − g i L (7) with where S ds,0 is the spectrally highly resolved deep space measurement, L is the forward calculated radiance, and S ds,i−1 is calculated from S ds,i−1 by spectral smoothing with a 4 cm −1 wide running average.Typically values for i max in the order of 5 have been shown to be adequate.Then, a linear fit is applied to obtain a first-order correction of the initial atmospheric parameters: As x we used altitude-constant scaling factors for the trace gas profiles exhibiting major signals in the spectrum and K are the Jacobians of the forward-calculated spectrum L with respect to the fit parameters x. x is subsequently used to update the simulated radiances: L cor is finally applied as corrected input for the iterative deep space correction (Eqs.7 and 8).
Resulting radiances are shown as red curve in Fig. 6.Obviously, many signatures of trace gases which were present in the initial measurement are strongly reduced -often below the measurement noise level.There are, however, some spectral regions where atmospheric residuals still show up after the correction in some cases, especially around the major ozone band at about 1040 cm −1 .Therefore the baseline between 995 and 1070 cm −1 is replaced by a parabola determined by a polynomial fit of the regions between 980 and 995 cm −1 and between 1070 and 1090 cm −1 .This additional correction makes microwindows up to 1012 cm −1 usable for the retrieval.The range between 1012 and 1070 cm −1 has been excluded from the subsequent retrievals of atmospheric parameters.

Spectral calibration
Spectral calibration is the assignment of an absolute spectral position to each sampling point.Experience has shown that for a good trace gas retrieval the spectral calibration accuracy should be better than 10 ppm (Friedl-Vallon et al., 2014).The main factors for a correct spectral calibration are the knowledge of the reference laser wavelength and the off-axis angle for each pixel.These two parameters are determined from measurement data and used for spectral calibration.An error in the knowledge of the laser wavelength leads to a corresponding error in the abscissa values of the spectrum.This error is pixel independent.The optical path difference (OPD) in the interferometer is dependent on the off-axis angle of the rays.The path through the interferometer is longer for off-axis rays as compared to the on-axis ray, but the optical path difference is actually shorter.The OPD for radiation falling on an off-axis pixel [i, j ] is given as where α i,j is the off-axis angle of the pixel [i, j ] and x 0 is the OPD for the optical axis (e.g.Davis et al., 2001).
If the same OPD is assumed for all pixels, the apparent position of a spectral line varies across the detector array according to with σ 0 being the line position for the on-axis position.
The off-axis angle for each pixel can be calculated geometrically from the distance r i,j of the pixel centre to the optical axis position on the array, and the image distance b, i.e. the distance between the focusing lens and the detector array.Ideally, when the detector is placed perfectly in the focus of the lens, the image distance is equal to the focal length.Figure 7 shows a schematic drawing of the off-axis angle.This simplified model neglects higher-order effects like chromatism.These effects are considered in the error budget (Sect.8.2).We calculate the off-axis angle for each pixel as In total, four parameters are needed for the spectral calibration of each pixel: -the wavelength of the reference laser; -the horizontal and vertical position of the optical axis on the detector array; -the image distance.
All these parameters may vary with the temperature of the instrument.It is mandatory to determine the spectral calibration parameters from in-flight measurements, because they are not stable enough in time to be determined from lab measurements on ground.The parameters are derived from atmospheric limb or deep space measurements taken with the chemistry mode spectral resolution (i.e. an MOPD of ±8 cm), using selected CO 2 lines with well-known spectral positions in the spectral range between 940 and 980 cm −1 (see Table 1).These lines are measured during flight with a relatively good signal-to-noise ratio, they are well isolated over a large altitude range and are globally available.
The spectral calibration parameters are deduced from a set of radiometrically calibrated spectra.For the generation of this data set, a first guess for the laser wavelength is used, and the off-axis effect is not taken into account.In this data set, the apparent position of a certain spectral line moves to smaller wavenumbers for pixels that are farther away from the optical axis (Eq.12).For each of the CO 2 lines listed in Table 1, the apparent line position is evaluated for each pixel.The interferogram is extended by zeros on both sides, such that the spectrum is strongly oversampled and reconstructed on a grid that is about 2 orders of magnitude finer than the original grid.The abscissa value of the line maximum of the oversampled spectrum is taken as line position.The line position across the detector array gives a two-dimensional bellshaped curve (Fig. 8, left).The position of the optical axis on the detector array is given by the maximum of this curve.A second-order 2-D polynomial, which is a good approximation to the cosine dependency of the values, is fitted to the data and the maximum of the fitted curve is taken as optical axis position.
Once the optical axis position is known, the distance r i,j of each pixel centre from the optical axis is calculated, and the line position found for each pixel is expressed as a function of this distance.Using Eqs. ( 12) and ( 13) and rearranging the terms, we get A curve fit of the form gives the fit results for b = √ a 0 and for σ 0 = √ a 1 .The actual laser wavelength is calculated from the deviation of σ 0 to the true line position of the respective line as given in the HITRAN database (Rothman et al., 2009): where λ laser is the laser wavelength, λ apriori is the first guess of the laser wavelength which was used to generate the spectra, σ 0 is the line position for the optical axis as derived from the measurement, and σ H is the true line position from the HITRAN database.The fit of the spectral calibration parameters is done separately for each of the lines listed in Table 1.The mean values over all spectral lines are used for spectral calibration.

Data processing
This section describes in more detail the processing steps which were briefly presented in Sect.3.

Level 0 processing
The level 0 processing resamples the interferograms on a space-equidistant grid including nonlinearity correction and spectral calibration.

Quality control
The time stamps within one cuboid are checked for equidistance.The number of time ticks between two consecutive frames should not differ by more than ±1.Any larger difference points towards lost frames.If a larger difference is detected, the file is discarded from further processing.

Spike detection and correction
The serial data link between IR detector front-end electronics and interferometer electronics is sporadically desynchronised, producing spikes in the interferograms.The immunity against the EMC disturbance depends on the gimbal yaw position of the instrument.In particular blackbody measurements, where the instrument looks approximately in the flight direction, are affected by spikes.These spikes need to be identified and corrected, as a single spike in the interferogram affects all data points of the spectrum.Two different methods are used to identify the spikes.One method is based on the identification of a specific pattern in affected frames.This method identifies spikes independently of the measurement type and the spike intensity.The second method is statistical and allows us to identify single large spikes that do not show the specific frame pattern.
If a frame is affected by spikes, the distribution of the affected pixels usually shows a specific pattern (see Fig. 9).In the case of a spike several pixels in a region show exactly the same value.In nominal measurements the probability that neighbouring pixels show the same value is very low because of measurement noise and the different sensitivity properties of each pixel.The strategy for searching spikes is a loop over all pixels in each frame.If at least one horizontal neighbour of a pixel has the same value as the pixel itself, a pixel is marked as spike candidate.A spike event is detected when at least for two rows the number of spike candidates per row passes an empirical threshold.
Manual inspection has shown that about 90 % of all spike events exhibit the pattern and are thus detected by this method.But if only one or few pixels are affected, this method fails.A second method, which does not use the characteristic pattern, is therefore applied to the data.For this method, all interferograms x i are normalised to the same scale with xi = where x i and σ i are the mean and standard deviation of the ith interferogram.The normalisation compensates for different sensitivities of the individual pixels.After normalisation each frame should have a similar internal variance.A spike which is much larger than typical read-out noise will increase the variance compared to non-affected frames.
Frames with a 9-times higher internal variance as compared to their neighbours are flagged as spike contaminated.Pixels with values deviating by more than three times the current standard deviation from the frame mean are then flagged as spiked.Close to the ZOPD the variation from frame to frame is generally high, therefore a region of ±0.06 cm around the peak is excluded from this method.
Figure 9 shows two examples of frames affected by spikes with the typical pattern which is identified by the first method.In the left picture, the affected pixels show the maximum value of 16 383 and nearly all columns are affected, whereas in the right picture, only few columns in the left part of the picture are affected, and the values of the affected pixels (light blue for the upper two and dark blue for the lower two rows of each affected block) are only slightly different from the values of the non-affected pixels, leading to spikes with small amplitudes.Figure 10 shows the normalised mean (top) and standard deviation (bottom) of a raw interferogram of an atmospheric measurement as calculated following the second method.The green curve shows three times the standard deviation.The red crosses mark identified spikes.
Detected spikes are corrected by replacing the spike with the mean value of the neighbouring frames.This correction method does not work for spikes close to the ZOPD (±0.02 cm).If a spike is detected in this range, the measurement is discarded.Outside the range close to the ZOPD, the effect of the corrected spike on the spectrum is well below noise level and the interferogram can be used.

Nonlinearity correction
The nonlinearity correction is performed by assigning to each interferogram value the value that would have been measured with a linear detector system (see Sect. 4.2).

Abscissa determination
In principle, the abscissa in the space domain can be chosen freely.Ideally, the sampling is chosen such that it corresponds to the measurement sampling in the time domain: x is the sampling in the space domain, v is the mean velocity and t is the sampling in the time domain.For the current flight configuration of GLORIA with a mean velocity of 1.27 cm s −1 and a sampling frequency of 6281 Hz, the natural sampling in the space domain as calculated from Eq. ( 18) would be 2.04 µm.A sampling denser than the natural sampling increases the number of points without adding new information, whereas a coarser sampling reduces the spectral bandwidth, leading to aliasing effects if out-of-band signals are not filtered.For an optimised processing (see Sect. 7) it is furthermore desired that the sampling distance is a multiple of 10 −4 nm and that the MOPD is hit by a sampling point.The sampling grid in the space domain is chosen as x = 2 µm, which is slightly smaller than the natural sampling of 2.04 µm.

Interferogram alignment
Using the complex calibration approach proposed by Revercomb et al. (1988), no explicit phase correction needs to be performed.It must be ensured, however, that the instrumental phase is stable and that the ZOPD position is the same for calibration and scene measurements, such that any instrument contribution cancels out during calibration.In real measurements, the ZOPD position may change from measurement to measurement due to a drift of the reference laser wavelength or due to fringe count errors (FCEs).FCEs occur when the laser fringes are not counted properly during turnaround, leading to a wrong starting point of the data acquisition.
A shift in the interferogram domain corresponds to a linear phase change in the spectral domain.Therefore the shift of the ZOPD position can be determined through the phase of the spectrum.Since the effects leading to a change of the ZOPD position are pixel independent, it is sufficient to determine the shift between interferograms for a single pixel.In order to enhance the signal-to-noise ratio, the interferograms of 11 × 11 pixels at the centre of the detector array are co-added.
In order to determine the linear phase change, the total phase has to be calculated and the linear phase change due to an interferogram shift has to be separated from other phase effects.In emission FTS, the calculated phase of the complex spectrum is strongly varying with the signal from the atmosphere, because instrument contributions, which are out of phase, are of the same order of magnitude as the atmospheric signal (Revercomb et al., 1988;Kleinert and Trieschmann, 2007).Figure 11 shows the real part of an uncalibrated but phase-corrected blackbody spectrum (top), the real part of a deep space spectrum (middle) and the imaginary part of a deep space spectrum (bottom).In the case of the blackbody spectrum, the instrument contributions are negligible for the phase calculation, because the signal coming from the source is one order of magnitude larger than the radiation emitted from the instrument.The calculated phase reflects the instrumental properties (e.g.dispersion of the beam splitter) which vary slowly with wavenumber (Fig. 12, black curve).In the case of the deep space spectrum, the situation is different.The instrument baseline is negative in a wide spectral range, and the beam-splitter emission, which is found in the imaginary part, gives a considerable contribution around 800 cm −1 .Therefore the phase angle changes rapidly with wavenumber, reflecting more the atmospheric spectrum than the instrumental phase (Fig. 12, red curve).In order to identify a linear phase change between different measurements, a spectral region has to be used where the contribution from the source is strong compared to the instrument self-emission for all types of measurement.In the case of GLORIA, the spectral range between 1025 and 1060 cm −1 is best suited for this purpose, because in this range the atmospheric signal is high due to the strong ozone signature, while instrument contributions are comparably low.As seen in Fig. 12, the phase in that range is the same as for the blackbody measurement.
The phase, which is calculated for a reference interferogram (usually a blackbody interferogram), is called the standard phase.For each interferogram, the difference to the standard phase is calculated, and a linear fit is applied in the range between 1025 and 1060 cm −1 .If the phase difference is purely caused by a shift in the interferogram domain, the intercept of the phase difference must be a multiple of 2 π .The calculated intercept is rounded to the next multiple of 2 π , and the final slope is calculated using this intercept and the mean phase value of the range between 1025 and 1060 cm −1 .The slope directly gives the shift in the interferogram domain.The shift is subtracted from the abscissa determined above.

Delay
The time stamps of the infrared data exhibit a certain delay with respect to the reference laser data.This is due to different signal travelling times in the electronics as well as the finite integration time.The delay t delay is given as f frame is the frame frequency, t int the integration time, t reset is the reset time for the integration capacitor of the detector, and t L is the runtime of the reference laser signal.
The frame frequency and the reset time are dependent on the array size used.The delay is calculated for the centre of the integration interval.
The delay is corrected by subtracting the delay time from the time stamps of the infrared signal.

Spectral calibration
The correct assignment of the spectral axis for a pixel on the optical axis is ensured through the use of the correct laser wavelength when building the space axis from the reference laser data.The laser wavelength is one of the spectral calibration parameters determined from selected atmospheric or deep space measurements.
In order to obtain the same abscissa for all pixels, the abscissa of each interferogram is divided by cos(α i,j ) (see Eq. 11) according to the off-axis angle of the respective pixel.The scaling factor for each pixel is stored in a lookup table, such that the cosine has to be calculated only once.

Interferogram resampling
Once the desired abscissa in space is defined, the interpolation points in time are calculated from the time-space relation given by the reference laser.Linear interpolation between adjacent points is sufficient to calculate the time for the desired point in space.The algorithm used to transfer the interferogram from the time domain into the spatial domain is an approximate Whittaker-Shannon interpolation (Shannon, 1949) as proposed by Brault (1996).The actual resampling is done by convolving the measured interferogram with an apodised sinc kernel.

Level 1 processing
The level 1 processing comprises the Fourier transform of the interferograms and the complex radiometric calibration.

Fourier transform
The Fourier transform calculates complex spectra from the resampled interferograms.In the standard processing, the interferograms are multiplied by the Norton-Beer strong apodisation function (Norton andBeer, 1976, 1977), since apodised spectra are used for the temperature and trace gas retrieval.

Radiometric calibration
Radiometrically calibrated spectra are obtained by dividing the atmospheric spectrum by the gain and subtracting the offset following Eq.(3).Gain and offset are drifting in time, therefore they are interpolated linearly between two calibration sequences to the time of the atmospheric measurement.
Ideally, the meaningful information (i.e. the measured atmospheric radiance) is contained in the real part of the complex calibrated spectrum while the imaginary part contains only noise.Any signature in the imaginary part suggests a measurement problem (e.g.pointing variations during interferogram acquisition or severe misalignment) or a processing problem (e.g. a phase shift between atmospheric and calibration measurements).
The real part of the calibrated spectrum serves as input to the temperature and trace gas retrieval.

Processor optimisation
The data processing as described in Sect.6 has to be performed for each individual pixel, i.e. several thousand times per measurement.Therefore especially time consuming operations like the interferogram resampling have to be optimised.
Since the processing is still under development, a twotrack strategy is followed: at first, the processing is programmed in IDL (Interactive Data Language), using the heritage of the processor that was developed for the balloonborne and airborne MIPAS instruments (Friedl-Vallon et al., 2004;Woiwode et al., 2012).The IDL processor with the character of a prototype allows for diagnostics, visualisation etc. and is easy to modify for testing purposes.It is used for quality control of the data, the definition of the abscissa and the interferogram shift determination.These processing steps need to be done only for one or few pixels per measurement.The computationally expensive resampling of the interferograms, which has to be done for each pixel individually, is encapsulated in a stand-alone C programme, which also performs the nonlinearity correction and spectral calibration and takes the externally computed delay and interferogram shift into account.
In view of a raw data rate of 73.6 MiB s −1 , the processor has to be optimised in order to approach or possibly even exceed the data acquisition speed.This is especially important for possible future long-duration balloon or satellite missions where it is not possible to store all the raw data.Then on-board processing may become necessary for data reduction or even for operating the experiment interactively from ground.

Description of the operational level 0 processor
Mathematically, the interpolation corresponds to a convolution of each pixel's interferogram data with a sinc function.This is implemented in the level 0 processor using a set of finite-width convolution kernels.For the sake of computing performance optimisation, each kernel consists of only 16 data points.The convolution kernel is apodised using a Kaiser filter in order to suppress the Gibbs oscillation.The different samplings of the convolution kernel are pretabulated to avoid computationally expensive runtime calculations of the sinc values.In order to reduce the size of the whole lookup table, the kernel values are stored in an integer (i.e.fixed point) representation.For each sample on the target axis, the correct kernel is chosen from the table depending on the position of the new sampling point.
In order to avoid cache misses (i.e. to ensure that the data needed for the convolution are more likely to be located in the cache), a prerequisite for the efficient convolution of the data is that the interferograms are contiguous in memory.Unfortunately, due to the nature of the data acquisition, the input is structured in a transposed form.The transposition of the input data is therefore the first processing step, followed by the interpolation setup and, finally, the convolution itself.
To avoid nonlinear memory access, a data structure has been specifically designed which maps each point of the predefined output abscissa to the correct input data and convolution kernel.
To ensure performance and compatibility, the level 0 processor is written in the C programming language.On top of the core processor, a Python interface has been added for more efficient interoperability with other processing steps.

Performance
The runtime of the level 0 processor is, apart from file input/output, dominated by two calculations.The first is the transposition of the measured image, the second is the convolution itself.
Initial optimisation targets were found by analysing both algorithms.Fortunately, both are by their nature very well suited for shared-memory parallelisation, which has been realised using the OpenMP protocol.Further potential for optimisation was discovered by analysing the assembly code generated by the C compiler.Using the SSE instruction set, the vector registers of modern CPUs could be used for optimised reimplementations of both the cuboid transposition and resampling.
Both approaches, parallelisation and vectorisation, were combined and resulted in a high-performance system able to match and exceed the data rate of the instrument.All following runtime measurements have been performed using a single GLORIA chemistry mode measurement.The cuboid file in question is about 945 MiB1 in size and took 12.8 s to record.The CPU of the benchmark computer was an AMD Opteron 6128 with eight cores clocked at 2000 MHz.In order to eliminate the impact of disk I/O, the relevant files have been read from and written to memory using tmpfs, a RAMbased file system for Linux.
Figures 13 and 14 show the runtimes and speedup factors for the convolution, the cuboid transposition and the total processing for the non-vectorised and the vectorised implementations, respectively.The non-vectorised version takes about 56.4 s to run without parallelisation, which is reduced to 12.4 s using eight concurrent threads.Note that this already slightly exceeds the real-time speed target.With the vectorisation (SSE) these timings reduce further to 36.9 s (single-thread) and 9.5 s (eight threads), respectively.
The most efficiently parallelised algorithm is by far the convolution, whose speedup scales almost linearly with the number of threads employed.The cuboid transposition involves a larger non-parallel overhead, making it scale less ideally.The use of SSE instructions mitigates the benefits of parallelisation only slightly, so that using both in conjunction is very effective.
The cuboid transposition benefits most from vectorisation with relative speedups larger than 4 in single-thread and still larger than 3 in eight-thread mode.The speedup is expected to become less prominent when more threads are used concurrently because memory bandwidth remains limited.In contrast, the speedup for the convolution is almost constant at 25 %, while the processor as a whole runs between 31 and 54 % faster with SSE enabled.
In total, the C code with all its optimisations is more than 2 orders of magnitude faster than our IDL prototype.

Outlook
Aside from x86-64, other computational architectures have been considered.Preliminary implementations of the level 0 processor on NVIDIA GPUs have been tested, but resulted in only a small overall performance gain of about 5 % compared with the CPU version.A more detailed study revealed that only 33 % of the GPU processing time was used for the computations, i.e. the transfer of the large data volumes between system and graphics memory resulted in a bottleneck.In the future, an implementation on a shared-memory CPU/GPU system such as the AMD Fusion platform may be investigated for better results.
The current level 1 processing uses the IDL routines developed for MIPAS with slight modifications.In the long run, this software is too inefficient for use with large volumes of GLORIA data and the degree of automation is limited.An integrated level 0 to 1 processing suite, which is implemented mainly in the Python programming language and includes optimisations for performance-critical sections, is under development.

Results
This section presents the quality of the calibrated spectra.Figure 15 shows as an example the real and imaginary part of a calibrated spectrum from a chemistry mode measurement after binning of one horizontal line (45 pixels used).The standard deviation of the imaginary part, which is a measure of the NESR, is about 15 nW cm −2 sr −1 cm.

Radiometric accuracy
An estimate for the radiometric accuracy can be achieved by analysing the calibrated spectra of the hot blackbody measurements.Since the data of the hot blackbody are not used for calibration in the current calibration approach, they provide independent measurements from a known source.
Figure 16, left, shows the mean relative radiance error in the spectral range from 900 to 950 cm −1 for each pixel, calculated from one calibration sequence of the HALO flight on 26 September 2012.Black and brown pixels exhibit an error of more than ±1.5 %.If a threshold of 1.5 % is applied to the data, less than 10 % of the pixels have to be discarded.This threshold provides a good compromise between the accuracy requirement of 1 to 2 % and the number of discarded pixels.Besides the outliers, the calibrated blackbody data show a positive bias of less than 1 % which is strongest in the middle of the detector.Such a bias can be evoked by temperature inhomogeneities of the blackbody surface.Currently the mean over five temperature sensors across the blackbody area is taken as blackbody temperature for all pixels, although the sensors show a temperature variation of approximately 0.3 K across the surface.This error is currently considered acceptable, but investigations are under way as to how the radiometric accuracy can be improved by taking the temperature inhomogeneity into account.Furthermore, the temperature homogeneity will be improved for future flights through improvements of the thermal control software and the operation procedures.

Drift of radiometric gain and offset
The drifts of the radiometric gain and offset are evaluated for the HALO flight on 13 September 2012.For the evaluation of the drift, a mean gain and offset function over the whole detector array is calculated for each calibration sequence.Then the mean gain and offset over all calibration sequences is calculated, and for each calibration sequence, the deviation of gain and offset from the mean is determined.Figure 17 shows the mean gain function over the flight in the upper panel, and in the lower panel it shows the relative deviation from the mean for each calibration sequence.The gain function varies by less than ±2 % over the nominal spectral range, showing that the calibration frequency is sufficiently high for the gain determination.Problems in the baseline determination around the strong ozone band at 1040 cm −1 (see Sect. 4.3) also affect the gain function in this range.
Figure 18 shows the offset, representing the instrument self-emission.The offset is negative, because instrument contributions having a phase of π with respect to the incoming radiance are stronger than instrument contributions which are in phase with the incoming radiance.The lower panel shows the difference from the mean spectrum for each calibration sequence.The drift above 900 cm −1 is governed by the temperature drift of the instrument.This drift can be captured well by the linear interpolation of the offset.The spectral feature below 900 cm −1 stems from the bulk emission of the germanium entrance window.The intensity of this feature is governed by the temperature of the entrance window, which changes faster than the instrument temperature.In this region a simple linear interpolation of the calibration measurements may be insufficient.Further characterisation work is planned in order to improve the radiometric calibration in this spectral range.

Spectral calibration accuracy
The accuracy of the spectral calibration is estimated by determining the line positions of selected lines after spectral calibration and comparing the values to the ones given in the HI-TRAN database.In order to enhance the signal-to-noise ratio, 48 × 8 pixels (horizontal × vertical) are co-added, such that the co-added data form an array of 1 × 16 super-pixels.For each of these super-pixels, the line positions of the 16 CO 2 lines listed in Table 1 as well as of two H 2 O lines at 1395.8 and 1399.2 cm −1 have been determined.The relative line position error expressed in ppm is shown in Fig. 19.The line indices 1 to 16 correspond to the ones given in Table 1, the line indices 17 and 18 correspond to the two H 2 O lines named above.For the CO 2 lines the line position error is mainly below 2 ppm, and no systematic effect relative to the optical axis is visible, showing that the simple cosine assumption worked sufficiently well.The H 2 O lines around 1400 cm −1 show a systematic deviation of about 5 ppm.This is still within the error budget of 10 ppm, but it shows that the simple scaling of the spectral axis only compensates first-order effects and that higher-order effects (e.g.chromatism) cause a systematic bias above noise level.

Stability of the spectral calibration parameters
In order to assess the stability of the spectral calibration parameters throughout a flight, the parameters have been determined from all deep space measurement sequences, separately for the forward and backward sweep direction (Sha, 2013).Figure 20 shows the variation of the spectral calibration parameters over time during the HALO flight on 13 September 2012.The parameters from top to bottom are laser wavelength, vertical and horizontal position of the optical axis (Opaxis_v and Opaxis_h), and image distance.The variation is mainly due to the temperature changes of the instrument.Because of the pressure-dependent refractive index of air, the apparent laser wavelength furthermore varies with pressure and thus with the flight altitude.The error bars denote the standard deviation of the results obtained from the 16 different spectral lines.The results for forward (index f) and backward (index b) do not show systematic differences.
When taking one parameter set for the whole flight, the maximum error for the corner pixel is about ±5 ppm.This error can be reduced by using time-dependent spectral calibration parameters.Currently this error is considered acceptable but an interpolation of the parameters may be included in the processing at a later stage.

Pixel binning
For the retrievals, all pixels of one row are co-added in the current approach to enhance the signal-to-noise ratio.Before co-addition, each pixel is characterised in terms of radiometric accuracy, in order to identify and discard bad pixels.Since the calibrated spectra are generated on a pixel basis, the pixel binning scheme can easily be modified for an optimised retrieval, applying various criteria or threshold values.When applying a threshold of 1.5 % for the radiometric accuracy as determined in Sect.8.1, about 10 % of the pixels are discarded.Furthermore, the first two columns, the last column and the first and last rows are discarded, because many pixels at the boundary of the array show a non-nominal behaviour.A possible binning mask is shown in Fig. 16, right panel.

Conclusions and outlook
This paper presents the level 0 to 1 processing of the data measured with the imaging FTS GLORIA.It is demonstrated that the radiometric and spectral requirements are generally fulfilled with the calibration concept and the processing approach for most of the pixels.About 15 % of the pixels are currently discarded due to an excess radiometric error or due to boundary effects of the detector array.
The radiometric properties for all pixels, namely the nonlinearity characteristics, shall be further improved by dedicated characterisation measurements.The better characterisation shall then reduce the number of discarded pixels.
Further improvements of the radiometric accuracy can be achieved by a better characterisation of the temperature inhomogeneity of the blackbodies and by optimising the thermal control parameters of the blackbodies for future flights.For the spectral region around 840 cm −1 improvements are expected through modelling of the emission from the germanium entrance window, taking its temperature into account.The combination of all three calibration points (hot and cold blackbody and deep space) should also allow for more precise gain and offset functions in spectral regions that are strongly affected by atmospheric signatures in the deep space spectrum, especially the range of the strong ozone band around 1040 cm −1 .
The spectral calibration accuracy is within the specification.The remaining spectral calibration error is considered uncritical.

Figure 3 .
Figure 3. Nonlinearity curve for a single pixel from laboratory measurements as a function of integration time.The red curve represents a virtual detector with linear properties.

Figure 4 .
Figure 4. Blackbody spectra of pixel[24,68]  measured during the flight on 13 September 2012 before (black) and after nonlinearity correction (red) for two different integration times: 54.9 µs in the upper panel and 112.8 µs in the lower panel.The blackbody temperature is 247 K.The maximum of each spectrum is scaled to 1.The reduction of the quadratic artefact can be seen in both cases showing that the nonlinearity characterisation curve determined from laboratory measurements can be applied to in-flight data.

Figure 5 .
Figure5.Mean signal of a normalised blackbody spectrum in the range between 20 and 400 cm −1 before (left) and after nonlinearity correction (right).The colour code denotes the relative intensity (the maximum of the spectrum is scaled to 1).The intensity of the artefact is clearly reduced after nonlinearity correction for most of the pixels.

Figure 6 .
Figure 6.GLORIA so-called deep space (10 • elevation) spectrum of a single detector pixel of 13 September 2012, 07:20 UT at 40.8 • S, 15.0 • E. For the spectrum six subsequent measurements have been co-added.Black: measured spectrum.Red: spectrum after removal of atmospheric contribution.Green: smoothed (4 cm −1 ) version of the red curve.

Figure 7 .
Figure 7.A ray diagram showing projections for the on-axis and an off-axis pixel.The [h, v] coordinates represent the on-axis position, [i, j ] represent the coordinates of an off-axis pixels, r is the distance from the optical axis, α is the off-axis angle, and b represents the image distance (sketch modified from Davis et al., 2001).

Figure 8 .
Figure 8. Apparent line position for the CO 2 line at 951.19 cm −1 across the detector array before (left) and after spectral calibration (right).The colours denote the shift with respect to the HITRAN position in cm −1 .
Figure 8 shows the determined line position for the CO 2 line at 951.19 cm −1 across the detector array before (left) and after spectral calibration (right) from a deep space measurement.While the bell shape is clearly visible in the left plot, the right plot does not show any systematic line position variation across the array.The pixel-to-pixel variation of the line position comes from the noise in the line position determination method.

Figure 9 .
Figure 9. Two examples of raw data frames affected by spikes.The colour code shows the intensity in ADC counts.For details, see text.

Figure 10 .
Figure 10.Normalised mean (top) and standard deviation (bottom) of a raw interferogram of an atmospheric measurement.The green curve shows three times the standard deviation.The red crosses mark identified spikes.

Figure 11 .
Figure 11.Raw spectra showing the real part of a blackbody spectrum (top), the real part of a deep space spectrum (middle) and the imaginary part of a deep space spectrum, which represents the beam-splitter emission (bottom).

Figure 12 .
Figure12.Phase calculated from a complex blackbody spectrum (black) and from a complex deep space spectrum (red).

Figure 13 .
Figure 13.Total runtimes and runtimes for the two most critical calculations with increasing level of OpenMP parallelisation, with and without the use of SSE vectorisation.Red: runtime for the convolution process, blue: runtime for the cuboid transposition, black: total runtime.Solid lines show the runtime with vectorisation, dashed lines without vectorisation.

Figure 14 .
Figure14.Speedup factor with increasing level of parallelisation for the total processor (black) and the most critical parts (red: convolution, blue: cuboid transposition).A speedup factor of α means that, at the indicated number of threads, the code in question runs α times as fast as the same code with only a single thread.Both the vectorised (solid lines) and the non-vectorised implementation (dashed lines) are shown for comparison.

Figure 15 .
Figure 15.Spectrum from a chemistry mode measurement on 28 September 2012, 10:45 UT at 29.67 • N, 43.41 • E for a tangent altitude of about 13.3 km after co-adding of one horizontal line (45 pixels used).The flight altitude was about 15.1 km.

Figure 16 .
Figure 16.Left: relative deviation of calibrated measurements of the hot blackbody (T = 256.6K) from the corresponding Planck function in the spectral range from 900 to 950 cm −1 .Right: bad pixel mask deduced from the radiometric error.Bad pixels are marked in black.

Figure 17 .
Figure 17.Mean gain function and relative deviation from the mean for the different calibration sequences during the HALO flight on 13 September 2012.

Figure 18 .
Figure 18.Mean instrument offset and deviation from the mean for the different calibration sequences during the HALO flight on 13 September 2012.

Figure 19 .
Figure 19.Relative line position error in ppm for selected spectral lines across the detector array, calculated for super-pixels consisting of 48 × 8 pixels.The line indices 1 to 16 correspond to the ones given in Table 1, the line indices 17 and 18 correspond to the two H 2 O lines given in the text.

Figure 20 .
Figure 20.Variation of the spectral calibration parameters over time during the HALO flight on 13 September 2012.

Table 1 .
The line position and intensity of the selected CO 2 lines from the HITRAN database used for the spectral calibration of the data measured during flight.