Prediction Model for Diffuser Induced Spectral Features in Imaging Spectrometers

Wide-field spectrometers for Earth Observation missions require inflight radiometric calibration, for which the sun can be used as a known reference. Therefor a diffuser is placed in front of the spectrometer in order to scatter the incoming light into the entrance slit and provide homogeneous illumination. The diffuser however, introduces interference patterns known as speckles into the system, yielding potentially significant intensity variations at the detector plane, called Spectral Features. There have been several approaches implemented to characterize the Spectral Features of a spectrometer, e.g. end-to-end 5 measurements with representative instruments. Additionally, in previous publications a measurement technique was proposed, which is based on the acquisition of monochromatic speckles in the entrance slit following a numerical propagation through the disperser to the detection plane. Based on this measurement technique we present a standalone prediction model for the magnitude of Spectral Features in imaging spectrometers, requiring only few input parameters and therefor mitigating the need for expensive measurement campaigns. 10


Introduction
Many current and future Earth Observation missions carry wide field spectrometer payloads such as the ENVISAT Medium Resolution Imaging Spectrometer (Olij et al. (1997)), the Sentinel-2 MultiSpectral Imager (Martimort et al. (2012)), the Sentinel-3a Ocean and Land Colour Imager (Nieke and Mavrocordatos (2017)), the Sentinel-4 UVN instrument (Clermont et al. (2019)), the Sentinel-5-UVNS instrument (Guehne et al. (2017)), or the GHGIS instrument of CO2M or former Car-15 bonSat (Fletcher et al. (2015)). These spaced based instruments require inflight radiometric calibration, for which the Sun can be used as a known reference. In order to ensure homogeneous illumination of the instrument a diffuser is used to scatter the incoming Sun light into the entrance slit. However, the diffuser introduces a statistical interference phenomenon known as speckles into the optical system. The speckles propagate through the disperser and are integrated at the detector plane, yielding intensity variations described as Spectral Features by van Brug et al. (2004). Since Spectral Features are of statistical nature 20 and cannot be mitigated by any post-processing steps, they may pose a significant contributor to the radiometric error budget. In particular for spectrometers with fine spectral resolution and strict demands to radiometric accuracy it is important to determine the severity of this error in an instrument even in early planing phases. The magnitude of this error is commonly described in terms of the Spectral Features Amplitude (SFA) first introduced by van Brug and Courrèges-Lacoste (2007).
The issue of diffuser induced Spectral Features in imaging spectrometers has first been documented by Richter et al. (2002) 25 and Wenig et al. (2004) in the context of the Global Ozone Monitoring Experiment (GOME). Ahlers et al. (2004) observed spectral oscillations caused by the onboard diffuser in the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) instrument, as well as van Brug et al. (2004). Several approaches for characterization and modeling have been proposed since, e.g. end-to-end measurements by van Brug and Courrèges-Lacoste (2007) as well as models for different speckle averaging effects derived by van Brug and Scalia (2012). However, a comprehensive and reliable model has 30 not been presented, yet. Isolating the Spectral Features by eliminating all other error sources in an representative end-to-end setup remains the main challenge to gain quantitative insights into the SFA dependence.
A different approach to quantify Spectral Features was taken by Burns et al. (2017) and improved by Richter et al. (2018).
It is based on the subsequent acquisition of monochromatic speckle patterns in the slit plane over a certain wavelength range, which are then propagated numerically through the disperser to the detection plane. Some simplifying assumptions are made 35 about the optical system which reduces the complexity and limits systematic error contributions. It is only limited by the SNR and the resolution achieved in the entrance slit and therefor capable of yielding comprehensive measurement data. Most important, it allows a step-by-step tracing of the speckle statistics from the slit to the detector plane.
Based on this SFA measurement technique we present a novel standalone SFA prediction model, which solely relies on mathematical descriptions of the speckle statistic and its SFA impact. It includes polarization effects of the diffuser, spatial and 40 spectral averaging as well as pixel averaging at the detector.
First we review the definition of the Spectral Features Amplitude (SFA) in Sect. 2. In Sect. 3 the revised SFA measurement technique used for our measurements is shown. We then present the standalone SFA prediction model in Sect. 4, which can be understood as a mathematical formulation of the SFA measurement technique. Finally, we compare the results of the prediction model to our measurement chain in detail in Sect. 5 to show its validity. In the last section we dicuss the applicability to a real 45 instrument.

Spectral Features Amplitude
The term Spectral Features Amplitude (SFA) was first proposed by van Brug and Courrèges-Lacoste (2007) as generic method to quantify diffuser induced "wiggles" in a spectrum measured by a space spectrometer instrument. They describe it as the magnitude of unwanted "features" that are left in the spectrum when subtracting all other features like emission lines from 50 the source and atmospheric absorption. The SFA value is then calculated as the standard deviation of the normalized signal over a certain spectral width, that includes multiple features. The SFA value holds information about the amplitude of features.
However, the data produced in this work, which is used to calculate the SFA, also allows for the estimation of the spectral extend of features. One usually may not draw conclusions about the absolute spectral positions of features with this approach.
We will show that for high performance volume diffuser, such as the one used in this work, lead to a spectral extend smaller 55 than the instrument detector pixel, which essentially allows for the treatment of the SFA as white noise at the detector level.
In this Sect. the used SFA measurement technique introduced by Burns et al. (2017) is presented in a revised state. The goal of this technique is the reduction of experimental complexity and therefore systematic error contributions during data acquisition.
First, the measurement principle is explained. Second, the used materials and the measurement procedures for the near infrared 60 (NIR) and the short wavelength infrared (SWIR) channel are presented. Figure 1 depicts the optical setup of an imaging spectrometer during solar calibration. The incoming Sun light is scattered by the diffuser. The scattering origin lies in the aperture plane with spatial coordinates g and h which is perpendicular to the light's direction of propagation. The angular field distribution at the aperture plane is imaged to the slit plane with coordinates x and 65 y. The light is collimated onto a dispersive element (e.g a diffraction grating), which is splitting it into its spectral components.

Principle
The spatial information in the y direction of the slit is transformed into spectral information at the detector with coordinate b by imaging the the diffracted beams of different wavelengths onto a 2D array detector. Beams of the same wavelength (within the spectral resolution) are assigned the same spectral detector coordinate b. The spatial information in x direction is retained in the detector coordinate a. We relate the coordinates via the simplified linear spectrometer equations where M x and M y are the respective magnification factors in x and y direction, k = db/dλ denotes the dispersion, and λ the wavelength. For these simplified equations to hold the magnification factors and the dispersion are assumed to be independent of the wavelength and the instrument point spread function (IPSF) is not accounted for. Additionally, a few properties of the 75 Sun's light and its detection need to be considered. It is assumed to be spatially coherent giving the distance from the Sun to the Earth and the limited acceptance angle of the spectrometer. Additionally, the temporal coherence is very short compared to the detector integration time, which is in the order of seconds. It follows that cross coherence terms of interfering fields of different wavelengths vanish and the net intensity distributions at the slit and the detector planes are well approximated by the superposition of monochromatic intensities. The Sun disk comprises of many incoherent point sources, which should be 80 considered for angular averaging contributions and is not part of this work, as it will only account for a single point source. For the purpose of the SFA measurement the sequence of optical components is subdivided into two parts. The first part ranging from the illuminated diffuser through the telescope to the entrance slit is represented by the optical setup in the lab. The second part comprises the rest of the optical system from the slit plane to the instrument detector plane. The data acquired in the first part is used as input for a numerical simulation of the optical setup after the slit plane. The setup layout is shown in Fig. 2.

85
The Sun is mimicked as a single field point with a tunable laser source, which is spectrally stabilized by a wavemeter, and illuminates the diffuser through a linear polarizer at normal incident with respect to the diffuser plane. The distance between the single mode fiber output and the diffuser is chosen such that the divergent beam illuminates the diffuser homogeneously over the size of the apertures. The second aperture blocks any unwanted angular contributions. A powermeter placed next to the diffuser records a fixed fraction of the emitted laser power. The telescope images the scattered light onto a 2D array detector 90 positioned in the focal plane. The focal plane of the telescope represents the slit plane in Fig. 1. The diffuser plane is tilted by 10°with respect to aperture and slit plane. This ensures that only scattered light contributes to the measurement. The telescope is aligned perpendicular to the aperture and slit plane.

Aperture plane (g, h)
Slit plane (x, y)  For a measurement, monochromatic speckle intensities are recorded subsequently over a wavelength range λ 1 ...λ N several times the spectral resolution λ res of the real spectrometer that is being mimicked. The spectral tuning step size ∆λ in between 95 images needs to be sufficiently small, in order to properly sample the change of the speckle patterns. The intermediate result is a three dimensional data set I slit (x, y, λ) consisting of a spectrum of monochromatic speckle images, where x and y are the spatial coordinates in the slit plane and λ is the wavelength. Every speckle image is mapped to a certain position (a, b) at the focal plane, where all images are summed up in intensity. The summation on intensity basis is justified as cross coherence terms involving interference of different wavelengths of actual Sun light will vanish for sufficiently long integration times. The

Aperture plane (g, h)
Slit plane (x, y)  summation procedure is detailed in Burns et al. (2017) and can be summarized as where slit coordinates are expressed in terms of the detector coordinates using Eq. (1) and (2) and the Heavyside function with Θ(y) = 0, y < 0 and Θ(y) = 1, y ≥ 0. The result of the sum is a two dimensional intensity distribution in the focal plane of the instrument I det (a, b). In a last step I det (a, b) is overlayed with the actual instruments detector pixel grid (ã,b)and intensities 105 belonging to the same pixel are summed. The SFA is calculated as standard deviation of the normalized detector pixel intensity distribution I det,binned (ã,b).

Materials and Procedure
Measurements are conducted in the NIR regime around 777 nm and in the SWIR regime around 1570 nm which represent wavelength bands with commonly monitored data products, such as water, clouds, CO2, Aerosols, or the O2 absorption which

Spectral Features Amplitude Prediction Model
The prediction model presented in the following is a mathematical formulation of the measurement method described in Sect. information about speckle averaging effects in the measurement chain lies in the intensity distributions. A single pattern I is sampled by a finite but sufficient amount of pixel, so that the individual pixel size is small compared to the speckle size. The magnitude of the speckle effect in I is described as the speckle contrast (Goodman, 2007, p. 28)

135
where σ I is the standard deviation and I is the mean value of I over all pixel. Under the general assumption, that the individual statistics of the underlying fields are circular complex Gaussian, a fully developed speckle pattern generated with linear polarized monochromatic light has a contrast of C = 1 (Goodman, 2007, p. 29). We adopt this assumption for this model.
The speckle contrast is reduced by several averaging effects introduced by the spectrometer instrument. A reduction of C is only achieved by the summation of intensity distributions showing a correlation smaller than unity. If the summation is on 140 amplitude basis (when the distributions can interfere), C is not reduced (Goodman, 2007, section 3.1.1). From this follows, that only distributions which can not interfere will impact C and are therefore subject of further discussions. Each one of the N independent averaging effects attributes to a certain amount of degrees of freedom M n or effectively uncorrelated intensity distributions, which can be combined to a total averaging factor M according to (Goodman, 2007, p. 186 145 The reduced speckle contrast will then calculate to In order to predict the contrast reduction we identified N = 3 contributors, which can be assigned to different steps of the SFA measurement chain: 1. Generation of monochromatic diffuse depolarized light in the aperture plane (g, h), 150 2. mapping intensities in the slit plane (x, y) to instrument detector positions (a, b),

Polarization Averaging
The laser source emits a single polarization state, which is ensured with the polarizer. The diffuse light leaving the volume diffuser can be treated as depolarized due to multi scattering (Lorenzo, 2012, p.85). This corresponds to two orthogonal 155 polarization configurations or two effective intensity distributions which can not interfere. Therefore step n = 1 introduces two degrees of freedom M polarization = 2 (Goodman, 2007, p. 49).

Spectral Averaging
Step n = 2 leads to spectral averaging at the detector. We recall the finding from Sect. 3.1, that the net intensities in the field planes (slit and detector) can be treated as superposition of monochromatic intensities for integration times greater than the 160 coherence time. Let us consider the acquired speckle intensities I n (x, y) and the underlying fields A n (x, y), which are related by I n = |A n | 2 . They are recorded at frequencies f n = c λn with a difference ∆f = f m − f n and c being the speed of light. The magnitude of the statistical change of subsequent speckle intensities I m and I n can be described in terms of the first order field correlation coefficient µ mn , with The field correlation is influenced by two effects, which in our case are both frequency dependent. The first effect is due to changing light paths through the diffuser medium. The second effect takes into account the spatial offset ∆b = k c ∆f at the detector plane between individual speckle patterns I n induced by the dispersion (see Eq. (2)). We start with the former contribution to the field correlation and follow the approach of Thompson et al. (1997a, b) and Webster et al. (2003), who used an analytic diffusion model to describe the light propagation in highly scattering, non-absorbing diffusers. The diffusion model yields the path length probability density function p(l) depending on the properties of the diffuser material, namely the scattering mean free path length l s , the refractive index of the material n s , and the thickness d. The characteristic function F l (∆f ) is the frequency dependent representation of p l (l) and is given by where c denotes the speed of light and α a constant factor taking into account the contribution of the tilted diffuser plane (e, f ) 175 with respect to the other planes and the specific geometry. The second contribution to the field correlation is due to changing spatial positions of speckle patterns which are distributed over the instrument detector in accordance with the spectral dispersion. This constitutes a spatial offset ∆b(∆f )) between the speckle intensities I n at the detector plane (a, b).The correlation of speckle fields, which are separated spatially by ∆b can be expressed as

180
where P (h) is the aperture function of the imaging system, h is the y-coordinate representation in the aperture plane, z is the distance between pupil and image plane, andλ the mean wavelength (Goodman, 2007, p. 169). Combining the two frequency dependent effects we can model the correlation between the speckle fields as µ mn (∆f ) = F l (∆f )Ψ(∆b(∆f )).
The accumulation of individual speckle patterns I n with field correlations µ mn at the detector can be interpreted as the sum-185 mation of partially correlated speckle intensities The amount of individual speckle intensities I n contributing to the sum at arbitrary detector coordinates (a, b) is equal to the ratio of the spectral resolution λ res with the step size ∆λ. This also applies to the mean intensities,
By diagonalization of J with a unitary linear transformation L 0 , the ensemble of correlated speckle fields is transformed to a basis with no correlation between them.
where † denotes the Hermitian transpose operation. The total mean intensity I det = n I n = n Ĩ n is conserved under this transformation but in general I n = Ĩ n . The complex coherence factor µ mn = |µ mn | exp(iΦ nm ) includes a phase Φ nm .
However, due to the specific construction of J, these phase terms can be omitted when calculating the eigenvalues (Dainty et al., 1975, section 4.7.2). Finally, for the spectral degrees of freedom we use the eigenvalues Ĩ n of the coherency matrix to Note that changing ∆λ to a smaller step size and therefor increasing N will not change the result of M spectral as long as ∆λ is sufficiently small to sample the covariance µ mn . The enabling property of the coherency matrix J is called Toeplitz, which implies an asymptotically behavior of its eigenvalues found by Grenander and Szegö (1958). Gray (2006) gives a simplified 205 prove in Corollary 2.1 and 2.2 that both, numerator and denominator in Eq. (16) converge for large N .

Detector Averaging
In step n = 3 an averaging due to the integration of the instrument detector pixel takes place. We already established, that the resultant intensity distribution at the detector I det (a, b) is given by the summation in Eq. (12). This effect impacts the speckle contrast if individual speckles are not sufficiently oversampled by the instrument detector pixel grid (ã,b). An analytical ex-210 pression for the degrees of freedom M detector introduced by stationary speckles in one detector pixel with relative coordinates (∆a, ∆b) is given by where A D is the area of a detector pixel, K D (∆a, ∆b) is the autocorrelation function of the detector pixel, and µ det (∆a, ∆b) is the field correlation at the detector plane (Goodman, 2007, p. 108). In order to accurately describe µ det one needs to account 215 for the evolution of the speckle size during the summation in Eq. (12). Let us consider a single speckle correlation area I 1 (S d /2 ≤ |x 1 | , S d /2 ≤ |y 1 | , f 1 ) with a spatial extend denoted by S d centered at (x 1 , y 1 ) in the slit. Its correlation relative to this position is described by 9 https://doi.org/10.5194/amt-2020-343 Preprint. Discussion started: 9 September 2020 c Author(s) 2020. CC BY 4.0 License.
with ∆x = x − x 1 and ∆y = y − y 1 being relative coordinates. The function F l (∆f ) introduced previously, characterizes how 220 the correlation area develops after n frequency steps at the same position, I n (S d /2 ≤ |x 1 | , S d /2 ≤ |y 1 | , f n ) with n > 1. In other words, it denotes the amount of spectral steps after which a single speckle seizes to exist at a fixed position. The initial position of the speckle at the detector is (a 1 , b 1 ) = (M x x 1 , M y y 1 + k c f1 ). The subsequent contributions relative to the initial position are shifted by k( c f1 − c fn ) and have a magnitude denoted by F l (f 1 − f n ). Therefore, the resultant speckle correlation function at the detector µ det (∆a, ∆b) is a convolution of |Ψ(∆a, ∆b(∆f ))| 2 with |F l (∆f )| 2 .

Predicted SFA
The predicted reduced speckle contrast at the instrument detector plane using Eq. (7) corresponds to the SFA and is

Results and Discussion
In the following we present and compare the SFA results from the measurement chain of Sect. 3 with the ones from the 230 prediction model of Sect. 4 in the NIR and SWIR regime. The values of relevant parameters used are depicted in Table 1. They were chosen to represent a proposed instrument for ESA's CO2M mission (Meijer et al. (2019)).
By calculating the frequency correlation of subsequent monochromatic speckle images I(x, y, c fn ) the fitted mean free path length of our diffuser sample is determined as l s = 53µm, which is close to the manufactures' specification of 55µm. This value may also reflect small instabilities of the speckle stationarity in the slit plane, which induce a smaller correlation between 235 subsequent speckle intensities and therefore a shorter effective mean free path length. Table 2 The results show a good agreement of prediction model and measurement. For the NIR regime the measured spectral averaging M spectral is 7% higher, than predicted. We assume, that detector noise is averaged in this step, which will yield a higher 250 effective averaging factor. The SWIR measurement spectral averaging factor shows a smaller deviation of 2%, which we also attribute to averaged detector noise. The measured detector averaging factor M detector for both wavelength regimes indicates no correlation of features between adjacent instrument detector pixel due to its high value (Goodman, 2007, p. 109). It is also dependent on detector noise, which explains the slightly higher averaging factors than predicted.
We demonstrated a comprehensive and numerical approach to quantify diffuser induced spectral features during solar calibration of space imaging spectrometers, which is based on established speckle theory concepts. We compared our prediction results with a current measuring method and observed a good agreement. The presented speckle averaging mechanisms are not a complete representation of the real in-orbit situation of an instrument. The effect of the Sun's disk, which consists of many incoherent points sources distributed over a 0.5 degree angle, needs to be taken into account as well as the averaging due to 260 the movement of the instrument relative to the Sun. Also, unlike the used laser point sources for the measurements, the Sun's light features an additional orthogonal polarization state, adding two polarization configurations to M polarization in the case of a highly scattering volume diffuser. The presented approach can be used for other diffuser types and optical geometries as well. It provides a solid starting point for future investigations into angular averaging mechanisms, which will complement the description of speckle reduction effects in imaging spectrometers of this type.

265
Data availability. The datasets generated and/or analyzed for this work are available from the corresponding author on reasonable request, subject to confirmation of Airbus Defence and Space GmbH. Competing interests. The authors declare, that they have no conflict of interest.