Prediction model for diffuser-induced spectral features in imaging spectrometers

Wide-field spectrometers for Earth observation missions require in-flight radiometric calibration for which the Sun can be used as a known reference. Therefore, 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 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 stand-alone prediction model for the magnitude of spectral features in imaging spectrometers, requiring only few input parameters and, therefore, mitigating the need for expensive measurement campaigns.

1 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), and the Greenhouse Gas Information System (GHGIS) instrument of CO2M or the former CarbonSat (Fletcher et al., 2015). These space-based instruments require in-flight 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 sunlight into the entrance slit. However, the diffuser introduces an interference pattern 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). When the solar spectrum is used for the calibration of the top of atmosphere measured spectra, the spectral features create a radiometric error as they alter the calibration function. Speckle effects are commonly known from applications involving highly coherent laser light, such as holographic imaging (Bianco et al., 2018) or laser speckle contrast imaging (Heeman et al., 2019). In general, diffuse sunlight does not yield a significant net speckle pattern when incident on a detector due to its broad spectrum. However, for spectrometers with a fine spectral resolution, the spectral width per channel is narrow enough to give rise to an amplitude of the spectral features that can be as large as other radiometric errors, e.g., due to straylight and polarization. Spectral features depend on various geometric conditions, implying that their exact position at the detector plane cannot be reliably predicted. This also renders any mitigating post-processing steps ineffective. Widely known speckle suppression techniques, such as the rotation or tilting of elements in the optical system (Goodman, 2007, Sect. 5), are only viable for on-ground calibration with a static setup. For space applications, additional moving parts are typically not implemented because they pose a supplementary risk of failure. In early planning phases, each contribution to the radiometric error needs to be quantitatively estimated to allow a global optimization of the instrument. This work presents a novel method to predict the radiometric error due to spectral features. As in most applications, the radiometric error is driven by the amplitude of the spectral features, and van Brug and Courrèges-Lacoste (2007) introduced the spectral features amplitude (SFA) as a standardized measure.

The issue of diffuser-induced spectral features in imaging spectrometers was first documented by Richter et al. (2002) and Wenig et al. (2004) in the context of the Global Ozone Monitoring Experiment (GOME). Ahlers et al. (2004) and van Brug et al. (2004) observed spectral oscillations caused by the onboard diffuser in the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) instrument. Several approaches to characterization and modeling have been proposed since. van Brug and Courrèges-Lacoste (2007) presented an end-to-end measurement setup, featuring a complete instrument, including a source, telescope optics, disperser, and detector. However, since spectral features need to be suppressed by the design of certain instrument parameters, a quantitative analysis with a representative spectrometer is usually difficult to achieve. van Brug and Scalia (2012) introduced models for different speckle averaging effects. However, a comprehensive and reliable model has not been presented yet. Isolating the spectral features by eliminating all other error sources in a representative end-to-end setup remains the main challenge with respect to gaining quantitative insights into the SFA dependence.

A different approach to quantifying 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 several spectral channels, which are then propagated numerically through the disperser to the detection plane. Some simplifying assumptions are made about the optical system, which reduces the complexity and limits systematic error contributions. It is only limited by the signal-to-noise ratio (SNR) and the resolution achieved in the entrance slit and is, therefore, capable of yielding comprehensive measurement data for most instrument designs. In particular, it allows for 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, stand-alone SFA prediction model, which solely relies on mathematical descriptions of the speckle statistic and its SFA impact. It includes the polarization effects of the diffuser, spatial and spectral averaging, and 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 stand-alone 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 discuss the applicability to a real instrument.

2 Spectral features amplitude

The term spectral features amplitude (SFA) was first proposed by van Brug and Courrèges-Lacoste (2007) as standardized way to quantify diffuser-induced wiggles in a spectrum measured by a space spectrometer instrument. They describe it as the magnitude of the features in a spectrum that are solely caused by the diffuser altering the solar reference spectrum. The SFA value is calculated as the standard deviation of the mean normalized detector signal over multiple spectral channels (see details in Sect. 3.1). The SFA value holds information about the amplitude of features. However, the data produced in this work, which are used to calculate the SFA, also allow 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 the instrument parameters used in this work lead to a spectral speckle extent smaller than the instrument detector pixel. This essentially allows for the treatment of the SFA as white noise at the detector level.

3 SFA measurement chain

In this section, the SFA measurement technique, introduced by Burns et al. (2017), is presented in a revised state. The goal of this technique is the reduction in 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 (NIR) and the short wavelength infrared (SWIR) channel are presented.

3.1 Principle

Figure 1 depicts the optical setup of an imaging spectrometer during solar calibration. The incoming sunlight 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 the coordinates x and y. The light is collimated onto a dispersive element (e.g., a diffraction grating), which splits it into its spectral components. The spatial information in the y direction of the slit is transformed into spectral information at the detector with coordinate b by imaging 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 the x direction is retained in the detector coordinate a. We relate the coordinates via the simplified linear spectrometer equations as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }a=M_{x}x,\text{(2)}& {\displaystyle }b=M_{y}y+k\mathit{\lambda },\end{array}$$

where M_x and M_y are the respective magnification factors in the x and y direction, $k=\mathrm{d}b/\mathrm{d}\mathit{\lambda }$ denotes the dispersion, and λ is the wavelength. For these simplified equations to hold the magnification factors and the dispersion, they are assumed to be independent of the wavelength and the spatial position (x,y). Also, the instrument point spread function (IPSF) is not accounted for. Sunlight is spatially coherent (Agarwal et al., 2004), and we can assume collimated light is illuminating the diffuser. Additionally, the sunlight's temporal coherence is on the order of femtoseconds (Hecht and Lippert, 2018) and much smaller than the typical detector integration time of several hundred milliseconds of an optical instrument. As a consequence, cross-coherence terms of interfering fields of different wavelengths average out, and the net intensity distributions at the slit and the detector planes are well approximated by the superposition of monochromatic intensities. The Sun disk is comprised of many incoherent point sources, which should be considered for angular averaging contributions and are not part of this work as they 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 are used as input for a numerical simulation of the optical setup after the slit plane. The setup layout is shown in Fig. 3. The Sun is mimicked as a single field point with a tunable laser source, which is spectrally stabilized by a wavemeter, and it 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 power meter 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 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.

Figure 1 Optical setup of an imaging spectrometer during solar calibration. The sequence of optical components is subdivided into two parts. The first part is covered by the experimental setup in the lab, starting at the illuminated diffuser and ending at the slit in the telescope focal plane. The second part numerically propagates the images recorded in the slit plane to the instrument focal plane.

Figure 2 Speckle patterns in the NIR band (a–c) and SWIR band (d–f) at different stages in the measurement chain. Panels (a) and (d) are an example of a monochromatic speckle pattern in the slit plane; panels (b) and (e) are the speckle pattern integrated at the detector plane using Eq. (3) and normalized to its mean, where the horizontal lines denote the instrument detector pixel grid $(\stackrel{\mathrm{̃}}{a},\stackrel{\mathrm{̃}}{b})$; and panels (c) and (f) are the final normalized integrated detector signal. The standard deviation taken over the pixel rows is the SFA.

Figure 3 Layout of the experimental setup for measuring diffuser-induced monochromatic speckle patterns in the slit plane. Single-mode fibers are used to connect the laser source with the wavemeter and the output via a fiber tab and are indicated by curved lines.

For a measurement, monochromatic speckle intensities are recorded subsequently over a wavelength range λ₁…λ_N several times the spectral resolution λ_res of the real spectrometer that is being mimicked. An example of such an monochromatic speckle pattern is shown in Fig. 2a and d in the slit plane. The spectral tuning step size Δλ in between images needs to be sufficiently small, in order to properly sample the change of the speckle patterns. The intermediate result is a 3D data set $I_{\mathrm{slit}}(x,y,\mathit{\lambda })$ 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 the intensity basis is justified as cross-coherence terms involving the interference of different wavelengths of actual sunlight, which will vanish for sufficiently long integration times. The summation procedure is detailed in Burns et al. (2017) and can be summarized as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{I}_{\det }(a,b)& ={\displaystyle \frac{\mathrm{\Delta }\mathit{\lambda }}{{\mathit{\lambda }}_{\mathrm{res}}}}\sum _{\mathit{\lambda }={\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_N}{I}_{\mathrm{slit}}\left({\displaystyle \frac{a}{M_x}},{\displaystyle \frac{b-k\mathit{\lambda }}{M_y}},\mathit{\lambda }\right)\\ & \times \mathrm{\Theta }\left(b-k\mathit{\lambda }\right)\mathrm{\Theta }\left(k{\mathit{\lambda }}_{\mathrm{res}}-b+k\mathit{\lambda }\right),\end{array}\end{array}$$

where slit coordinates are expressed in terms of the detector coordinates using Eqs. (1) and (2) and the Heaviside function with $\mathrm{\Theta }\left(y\right)=\mathrm{0},y<\mathrm{0}$ and $\mathrm{\Theta }\left(y\right)=\mathrm{1},y\ge \mathrm{0}$. The result of the sum is a 2D intensity distribution in the focal plane of the instrument I_det(a,b), which is depicted in Fig. 2b and e. In a last step, I_det(a,b) is overlain with the actual instrument's detector pixel grid $(\stackrel{\mathrm{̃}}{a},\stackrel{\mathrm{̃}}{b})$, and intensities belonging to the same pixel are summed, which are shown in Fig. 2c and f. The SFA is calculated as the standard deviation of the mean normalized detector pixel intensity distribution $I_{\det ,\mathrm{binned}}(\stackrel{\mathrm{̃}}{a},\stackrel{\mathrm{̃}}{b})$.

3.2 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 vapor, clouds, CO₂, aerosols, or the O₂ absorption, which are commonly used to calculate the effective path length and the air mass factor (see Irizar et al., 2019, Meijer et al., 2019, or Voors et al., 2017). The experimental setup is shown in Sect. 3.1. As light sources, they serve tunable monochromatic external cavity diode lasers with single-mode output and an integrated optical isolator. They are stabilized via a proportional–integral–derivative (PID) loop with feedback data from a wavemeter, which uses a Fizeau interferometer. Connecting fibers between the laser, wavemeter, and output serve single-mode (SM) fibers since the spectral tuning range will be narrow. Also, SM fibers will introduce no additional speckle in contrast to multi-mode fibers. A linear polarizer ensures polarization stability. The round diffuser plate has a diameter of 70 mm and a thickness of 3 mm. It is made of high-scattering fused silica HOD^®-500 material featuring inhomogeneities of 20 µm or less. The spatial coherence length of the Sun's light is around 60 µm in the NIR and 120 µm in the SWIR band, according to Divitt and Novotny (2015). Hence, the material, which was selected for the Sentinel-5/UVNS according to Irizar et al. (2019), is suited to the wavelengths in question and deemed a good choice for this study. The data collected with the power meter are used to normalize the acquired images. The round apertures are used to control the size of individual speckle correlation areas. The laser beam's uniformity at the diffuser plane was measured to be around 3 % in the NIR band and 6 % in the SWIR band over the size of the aperture. There were no additional speckle contribution by the fibers detected. The telescope has a focal length of f_tel=1100 mm. For the NIR, the laser source has a center wavelength of 780 nm and a nominal linewidth of $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{7}}\mathrm{nm}$. The Thorlabs 780HP SM fiber (SMF) type was used for all fiber connections in the NIR range. The charge-coupled device (CCD) detector features a 12.5 mm × 10.0 mm active area with 2750 × 2200 pixels of size 4.54 µm × 4.54 µm. The expected average 1D width of the speckle correlation area S_d in the slit plane is calculated using the following:

$$\begin{array}{}\text{(4)}& S_{\mathrm{d}}={\displaystyle \frac{\mathrm{2}\mathit{\lambda }{f}_{\mathrm{tel}}}{D\sqrt{\mathit{\pi }}}},\end{array}$$

which was derived by Goodman (2007). The aperture's diameter is set to D=10 mm, and, thus, the speckles are sampled by 19 pixels in one dimension, which is deemed sufficiently fine to ensure that the sampling with the detector does not alter the speckle pattern significantly. The laser wavelength was tuned over the range of 776.4–777.7 nm, with a step size of Δλ=1 pm. The step size is chosen so that there is a non-zero correlation between subsequent speckle images. For the SWIR measurement, the laser source and the detector and the fiber splitter are replaced, and the Thorlabs SMF-28 was used for the fiber connections. The SWIR laser source center wavelength is 1550 nm, with a single-mode output of nominal $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}\mathrm{nm}$ linewidth. The detector is a 640 × 512 pixel InGaAs camera with a pixel size of 15.5 µm × 15.5 µm. The tuning range was 1571–1577.5 nm, with a step size of Δλ=3.1 pm. The aperture's diameter is set to D=13 mm, which means a sampling of 10 pixels per speckle in one dimension. This is smaller than for the NIR configuration but is a trade off between sampling and SNR.

4 Spectral features amplitude prediction model

The prediction model presented in the following is a mathematical formulation of the measurement method described in Sect. 3. It relies on the determination of the speckle statistics at different steps of the measurement chain. The relevant physical 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 pixels, 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) as follows:

$$\begin{array}{}\text{(5)}& C={\displaystyle \frac{{\mathit{\sigma }}_I}{\langle I\rangle }},\end{array}$$

where σ_I is the standard deviation and 〈I〉 is the mean value of I over all pixels. Under the general assumption that the individual statistics of the underlying fields are a 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 in C is only achieved by the summation of intensity distributions showing a correlation smaller than unity. If the summation is on an amplitude basis (when the distributions can interfere), C is not reduced (Goodman, 2007, Sect. 3.1.1). From this, it follows that only distributions which cannot interfere will impact C and are therefore the subject of further discussions. Each one of the N independent averaging effects is attributed 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 (Goodman, 2007, p. 186) by the following:

$$\begin{array}{}\text{(6)}& M=\prod _n^N{M}_n.\end{array}$$

The reduced speckle contrast will then be calculated as follows:

$$\begin{array}{}\text{(7)}& C_{\mathrm{reduced}}={\displaystyle \frac{\mathrm{1}}{\sqrt{M}}}={\left(\sqrt{\prod _n^N{M}_n}\right)}^{-\mathrm{1}}.\end{array}$$

In order to predict the contrast reduction, we identified N=3 contributors, which can be assigned to different steps of the SFA measurement chain as follows:

1. generation of monochromatic diffuse depolarized light in the aperture plane (g,h) yielding a factor M_polarization;

2. mapping intensities in the slit plane (x,y) to instrument-detector positions (a,b) contributing a factor M_spectral; and

3. integrating the instrument-detector pixels with a factor denoted by M_detector.

The predicted reduced speckle contrast at the instrument-detector plane, using Eq. (7), corresponds to the SFA and is as follows:

$$\begin{array}{}\text{(8)}& \mathrm{SFA}=C_{\mathrm{reduced}}={\displaystyle \frac{\mathrm{1}}{\sqrt{M_{\mathrm{polarization}}{M}_{\mathrm{spectral}}{M}_{\mathrm{detector}}}}}.\end{array}$$

4.1 Polarization averaging

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

4.2 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 the superposition of monochromatic intensities for integration times much greater than the 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={\left|A_n\right|}^{\mathrm{2}}$. They are recorded at various wavelengths λ_n. The magnitude of the statistical change of subsequent speckle intensities I_m and I_n, with a wavelength difference $\mathrm{\Delta }{\mathit{\lambda }}_{nm}=\left|{\mathit{\lambda }}_n-{\mathit{\lambda }}_m\right|$, can be described in terms of the first-order field correlation coefficient μ_mn, with the following:

$$\begin{array}{}\text{(9)}& {\mathit{\mu }}_{mn}\left({\mathit{\lambda }}_n,{\mathit{\lambda }}_m\right)={\displaystyle \frac{\langle {\mathit{A}}_m{\mathit{A}}_n^{\ast }\rangle }{\sqrt{\langle I_m\rangle \langle I_n\rangle }}},\end{array}$$

where ∗ denotes the complex conjugate. 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Δλ_nm 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 Zhu et al. (1991), who presented an analytic equation for the angular correlation function of a slab geometry of a scattering media of thickness d, which can also be used for wavelength correlations as follows:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& F\left({\mathit{\lambda }}_n,{\mathit{\lambda }}_m\right)=\\ & {\displaystyle \frac{\left(d+\mathrm{2}B\right)\left[\sinh \left(z_{\mathrm{0}}\sqrt{Q}\right)+B\sqrt{Q}\cosh \left(z_{\mathrm{0}}\sqrt{Q}\right)\right]}{\left(z_{\mathrm{0}}+B\right)\left[\left(\mathrm{1}+B^{\mathrm{2}}Q\right)\sinh \left(d\sqrt{Q}\right)+\mathrm{2}B\sqrt{Q}\cosh \left(d\sqrt{Q}\right)\right]}},\end{array}\end{array}$$

where we set $Q=q^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}$ and $q=\sqrt{i\mathrm{6}\mathit{\pi }\left|\frac{\mathrm{1}}{{\mathit{\lambda }}_n}-\frac{\mathrm{1}}{{\mathit{\lambda }}_m}\right|\mathit{\beta }{n}_{\mathrm{s}}/l_{\mathrm{t}}}$, with n_s denoting the refractive index of the scattering material, l_t being the transport mean free path used for anisotropic multi-scattering systems, and β being a constant factor taking into account the contribution of the tilted diffuser plane (e,f) with respect to the other planes and the specific geometry. The value z₀ describes the average penetration depth after which the light is scattered for the first time. It does not have great influences in transmission geometry and is approximated with the transport mean free path l_t. We set α=0, thereby ignoring absorption. The parameter for the boundary condition is given by $B=l_{\mathrm{t}}\frac{\mathrm{2}(\mathrm{1}+R)}{\mathrm{3}(\mathrm{1}-R)}$, where R is the reflection coefficient which is calculated using the Fresnel equations. It accounts for internal reflection due to the index of refraction mismatch at the boundaries.

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 between the speckle intensities I_n and I_m at the detector plane (a,b). Keeping in mind Eq. (1) and (2) for the transformation from the slit to the detector plane, the correlation of speckle fields, which are separated spatially by Δb, can be expressed as follows:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\left.\mathrm{\Psi }(\mathrm{\Delta }a,\mathrm{\Delta }b)\right|}_{\mathrm{\Delta }a=\mathrm{0}}=\\ & {\left.{\displaystyle \frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\left|P(g,h)\right|}^{\mathrm{2}}{e}^{-i\frac{\mathrm{2}\mathit{\pi }}{\stackrel{\mathrm{̃}}{\mathit{\lambda }}z}(g\mathrm{\Delta }a+h\mathrm{\Delta }b)}\mathrm{d}g\mathrm{d}h}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\left|P(g,h)\right|}^{\mathrm{2}}\mathrm{d}g\mathrm{d}h}}\right|}_{\mathrm{\Delta }a=\mathrm{0}},\end{array}\end{array}$$

where P(g,h) is the aperture function of the imaging system and may be defined, e.g., for a circular aperture of diameter D, as $P(g,h)=\mathrm{circ}\left(\frac{\mathrm{2}\sqrt{g^{\mathrm{2}}+h^{\mathrm{2}}}}{D}\right)$, which is 1 inside the aperture and 0 otherwise. g and h are the a- and b-coordinate representations in the aperture plane, z is the distance between aperture and slit plane, and $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ is the mean wavelength (Goodman, 2007, p. 169). We set Δa=0 since we are only interested in the offset in the spectral direction. Combining the two effects, we can model the correlation between the speckle fields as follows:

$$\begin{array}{}\text{(12)}& {\mathit{\mu }}_{mn}\left({\mathit{\lambda }}_n,{\mathit{\lambda }}_m\right)=F\left({\mathit{\lambda }}_n,{\mathit{\lambda }}_m\right)\mathrm{\Psi }\left(\mathrm{\Delta }b\right).\end{array}$$

The accumulation of individual speckle patterns I_n with field correlations μ_mn at the detector can be interpreted as the summation of partially correlated speckle intensities as follows:

$$\begin{array}{}\text{(13)}& I_{\det }(a,b)=\sum _{n=\mathrm{1}}^{N={\mathit{\lambda }}_{\mathrm{res}}/\mathrm{\Delta }\mathit{\lambda }}{I}_n\left({\displaystyle \frac{a}{M_x}},{\displaystyle \frac{b-k\mathit{\lambda }}{M_y}},{\mathit{\lambda }}_n\right).\end{array}$$

The number 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, as follows:

$$\begin{array}{}\text{(14)}& \langle I_{\det }(a,b)\rangle =\sum _{n=\mathrm{1}}^{N={\mathit{\lambda }}_{\mathrm{res}}/\mathrm{\Delta }\mathit{\lambda }}〈I_n\left({\displaystyle \frac{a}{M_x}},{\displaystyle \frac{b-k\mathit{\lambda }}{M_y}},{\mathit{\lambda }}_n\right)〉.\end{array}$$

Using an established method by Bevan (2009) and Goodman (2007), we define a coherency matrix with entries $J_{nm}=\langle {\mathit{A}}_m{\mathit{A}}_n^{\ast }\rangle$ and use Eq. (9) to obtain the following:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& \mathbf{J}=\\ & \left[\begin{array}{cccc}\langle I_{\mathrm{1}}\rangle & \sqrt{\langle I_{\mathrm{1}}\rangle \langle I_{\mathrm{2}}\rangle }{\mathit{\mu }}_{\mathrm{1},\mathrm{2}}& \mathrm{\cdots }& \sqrt{\langle I_{\mathrm{1}}\rangle \langle I_N\rangle }{\mathit{\mu }}_{\mathrm{1},N}\\ \sqrt{\langle I_{\mathrm{1}}\rangle \langle I_{\mathrm{2}}\rangle }{\mathit{\mu }}_{\mathrm{1},\mathrm{2}}^{\ast }& \langle I_{\mathrm{2}}\rangle & \mathrm{\cdots }& \sqrt{\langle I_{\mathrm{2}}\rangle \langle I_N\rangle }{\mathit{\mu }}_{\mathrm{2},N}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \sqrt{\langle I_{\mathrm{1}}\rangle \langle I_N\rangle }{\mathit{\mu }}_{\mathrm{1},N}^{\ast }& \sqrt{\langle I_{\mathrm{2}}\rangle \langle I_N\rangle }{\mathit{\mu }}_{\mathrm{2},N}^{\ast }& \mathrm{\cdots }& \langle I_N\rangle \end{array}\right].\end{array}\end{array}$$

By the diagonalization of J with a unitary linear transformation L₀, the ensemble of correlated speckle fields is transformed to a basis with no correlation between them.

$$\begin{array}{}\text{(16)}& {\mathbf{J}}^{\prime }={\mathbf{L}}_{\mathrm{0}}{\mathbf{JL}}_{\mathrm{0}}^{†}=\left[\begin{array}{cccc}\langle \stackrel{\mathrm{̃}}{I_{\mathrm{1}}}\rangle & \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \langle \stackrel{\mathrm{̃}}{I_{\mathrm{2}}}\rangle & \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \langle \stackrel{\mathrm{̃}}{I_N}\rangle \end{array}\right],\end{array}$$

where † denotes the Hermitian transpose operation. The total mean intensity $\langle I_{\det }\rangle ={\sum }_n\langle I_n\rangle ={\sum }_n\langle \stackrel{\mathrm{̃}}{I_n}\rangle$ is conserved under this transformation but, in general, $\langle I_n\rangle \ne \langle \stackrel{\mathrm{̃}}{I_n}\rangle$. The complex coherence factor ${\mathit{\mu }}_{mn}=\left|{\mathit{\mu }}_{mn}\right|\exp \left(i{\mathrm{\Phi }}_{nm}\right)$ 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, Sect. 4.7.2). Finally, for the spectral degrees of freedom, we use the eigenvalues $\langle \stackrel{\mathrm{̃}}{I_n}\rangle$ of the coherency matrix to obtain the following:

$$\begin{array}{}\text{(17)}& M_{\mathrm{spectral}}={\left({\displaystyle \frac{\langle I_{\det }\rangle }{{\mathit{\sigma }}_{\det }}}\right)}^{\mathrm{2}}={\displaystyle \frac{{\left({\sum }_n\langle \stackrel{\mathrm{̃}}{I_n}\rangle \right)}^{\mathrm{2}}}{{\sum }_n\langle \stackrel{\mathrm{̃}}{I_n}{\rangle }^{\mathrm{2}}}}.\end{array}$$

Note that changing Δλ and, thus, changing 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 asymptotic behavior of its eigenvalues found by Grenander and Szegö (1958). Gray (2006) gives a simplified proof in corollary 2.1 and 2.2 that both the numerator and denominator in Eq. (17) converge for large N.

4.3 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. (13). This effect impacts the speckle contrast if individual speckles are not sufficiently oversampled by the instrument detector pixel grid $(\stackrel{\mathrm{̃}}{a},\stackrel{\mathrm{̃}}{b})$. An analytical expression for the degrees of freedom M_detector introduced by stationary speckles in one detector pixel with relative coordinates (Δa,Δb) is given by the following:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& M_{\mathrm{detector}}=\\ & {\left[{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{D}}^{\mathrm{2}}}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iint }}{K}_{\mathrm{D}}(\mathrm{\Delta }a,\mathrm{\Delta }b){\left|{\mathit{\mu }}_{\det }(\mathrm{\Delta }a,\mathrm{\Delta }b)\right|}^{\mathrm{2}}\mathrm{d}\mathrm{\Delta }a\mathrm{d}\mathrm{\Delta }b\right]}^{-\mathrm{1}},\end{array}\end{array}$$

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 for the evolution of the speckle size during the summation in Eq. (13). Let us consider a single speckle correlation area $I_{\mathrm{1}}(S_{\mathrm{d}}/\mathrm{2}\le \left|x-x_{\mathrm{1}}\right|,S_{\mathrm{d}}/\mathrm{2}\le \left|y-y_{\mathrm{1}}\right|,{\mathit{\lambda }}_{\mathrm{1}})$ with a spatial extent denoted by S_d centered at (x₁,y₁) in the slit. Its correlation relative to this position is described by the following:

$$\begin{array}{}\text{(19)}& \mathrm{\Psi }(\mathrm{\Delta }x,\mathrm{\Delta }y)={\displaystyle \frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\left|P(g,h)\right|}^{\mathrm{2}}{e}^{-i\frac{\mathrm{2}\mathit{\pi }}{\stackrel{\mathrm{̃}}{\mathit{\lambda }}z}(g\mathrm{\Delta }x+h\mathrm{\Delta }y)}\mathrm{d}g\mathrm{d}h}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\left|P(g,h)\right|}^{\mathrm{2}}\mathrm{d}g\mathrm{d}h}},\end{array}$$

with $\mathrm{\Delta }x=x-x_{\mathrm{1}}$ and $\mathrm{\Delta }y=y-y_{\mathrm{1}}$ being relative coordinates. The function F(λ₁,λ_n), introduced previously, characterizes how the correlation area develops after n wavelength steps at the same position, $I_n(S_{\mathrm{d}}/\mathrm{2}\le \left|\mathrm{\Delta }x\right|,S_{\mathrm{d}}/\mathrm{2}\le \left|\mathrm{\Delta }y\right|,{\mathit{\lambda }}_n)$, with n>1. In other words, it denotes the number of spectral steps after which a single speckle ceases to exist at a fixed position in the slit. The initial position of the speckle at the detector is $(a_{\mathrm{1}},b_{\mathrm{1}})=(M_x{x}_{\mathrm{1}},M_y{y}_{\mathrm{1}}+k{\mathit{\lambda }}_{\mathrm{1}})$. The subsequent contributions relative to the initial position are shifted in the b direction by k|λ₁−λ_n| and have a magnitude denoted by F(λ₁,λ_n). Therefore, the resultant speckle correlation function at the detector ${\left|{\mathit{\mu }}_{\det }(\mathrm{\Delta }a,\mathrm{\Delta }b)\right|}^{\mathrm{2}}$ is a convolution of |Ψ(Δa,Δb)|² with ${\left|F\left({\mathit{\lambda }}_n,{\mathit{\lambda }}_m\right)\right|}^{\mathrm{2}}$ as follows:

$$\begin{array}{}\text{(20)}& {\left|{\mathit{\mu }}_{\det }(\mathrm{\Delta }a,\mathrm{\Delta }b)\right|}^{\mathrm{2}}={\left|\mathrm{\Psi }(\mathrm{\Delta }a,\mathrm{\Delta }b)\right|}^{\mathrm{2}}\textcolor{red}{}{\left|F\left({\mathit{\lambda }}_n,{\mathit{\lambda }}_m\right)\right|}^{\mathrm{2}}.\end{array}$$

The symbol denotes the convolution.

5 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 prediction model of Sect. 4 in the NIR and SWIR regime. M_spectral can be interpreted as the average number of speckle patterns generated by the diffuser per spectral channel at the given wavelengths. M_detector can be understood as the average number of speckle correlation areas influencing the measurement in a detector pixel. The values of relevant parameters used are given in Table 1. They were chosen to represent a proposed instrument for the European Space Agency's (ESA) CO2M mission (Meijer et al., 2019). In Fig. 4, the measured Pearson correlations between speckle patterns I_n(λ_n) and I_m(λ_m), with respect to their relative spectral distance, are shown as blue stars for the NIR band at 776 nm and the SWIR band at 1571 nm, respectively. All wavelength shift combinations up to 0.1 nm are averaged for 120 images. Error bars are omitted, because the standard error of the mean is too small to be displayed. The red line denotes a fit to Eq. (10), leading to the values for l_t given in Table 1. They are in agreement with the supplier's given value of l_t=56 µm at 500 nm wavelength.

Table 1 Sample spectrometer parameters used for the measurement and prediction. They were chosen to represent a proposed instrument for the European Space Agency's CO2M mission (Meijer et al., 2019).

Figure 4 Measurement of the correlation function ${\left|F\left({\mathit{\lambda }}_n,{\mathit{\lambda }}_m\right)\right|}^{\mathrm{2}}$ in (a) the NIR (776 nm) and (b) the SWIR band (1571 nm). Blue stars denote the measured Pearson correlations between speckle patterns I_n(λ_n) and I_m(λ_m). All wavelength shift combinations up to 0.1 nm are averaged for 120 images. Error bars are omitted because the standard error of the mean value is too small to be displayed. The red graph denotes the fit of the measured data points to Eq. (10), with (a) $l_{\mathrm{t}}=(\mathrm{59.3}\pm \mathrm{0.4})$ µm and (b) $l_{\mathrm{t}}=(\mathrm{67.8}\pm \mathrm{0.5})$ µm.

Table 2 Comparison SFA results of the measurement chain with the prediction model. Measurement uncertainties are given in the 1σ interval.

Table 2 shows the SFA values of the measurements and predictions with their corresponding intermediate averaging factors M_polarization, M_spectral, and M_detector, as introduced in Sect. 4. Their counterparts from the measurement are deducted by calculating the speckle contrast at intermediate steps in the measurement chain. To verify the factor M_polarization=2, a linear polarizer is placed after the diffuser, and the measured speckle contrast, compared to the nominal case, rises by a factor of $\sqrt{\mathrm{2}}$. Additionally, the polarization axis is rotated to different random positions without changing the result, which confirms the assumption made in Sect. 4 that the light exiting the diffuser is depolarized. With an ideal measurement chain the speckle contrast expected in the slit plane for monochromatic polarized speckles is $C_{\mathrm{slit},\mathrm{ideal}}=\mathrm{1}$. This is the numerator in Eq. (8). The measured contrast in the slit is smaller, probably due to detector noise, as suggested by Postnov et al. (2019). Webster et al. (2003) attributed the reduced measured contrast to straylight from multiple reflections in their setup. The experimental factors M_{spectral,measured} and M_{detector,measured} are calculated in similar way, using the following relations:

$$\begin{array}{}\text{(21)}& {\displaystyle }{C}_{\mathrm{spectral},\mathrm{measured}}={\displaystyle \frac{\langle C_{\mathrm{slit},\mathrm{measured}}\rangle }{\sqrt{M_{\mathrm{spectral},\mathrm{measured}}}}},\text{(22)}& {\displaystyle }{C}_{\mathrm{detector},\mathrm{measured}}={\displaystyle \frac{C_{\mathrm{spectral},\mathrm{measured}}}{\sqrt{M_{\mathrm{detector},\mathrm{measured}}}}}.\end{array}$$

The results show a good agreement between the prediction model and measurement and are well within the estimated error margins, which are given in the 1σ interval. While the spectral extent of a speckle correlation area at the detector can be derived from the width of μ_det, one can also see from the values of M_detector that the extent must be small compared to a pixel. This allows for the treatment of the SFA as white noise on the detector level, which can be a goal of the instrument design. The uncertainties ΔM_{spectral,measured} and ΔM_{detector,measured} are estimated by considering fluctuations in the two contributors, F and Ψ, to the field correlation μ_mn in Eq. (12) and μ_det in Eq. (20). For F, the measured standard deviation of each data point in Fig. 4 is taken as a measure, which amounts to 1 %–2 % in the NIR and 2 %–3 % in the SWIR band. The average fluctuations of Ψ were determined by measuring its width over 120 speckle patterns with autocorrelation to be 1.3 % and 2.7 % in the NIR and SWIR, respectively. The impact on the M_{spectral,measured} and M_{detector,measured} is estimated by a Monte Carlo propagation of both uncertainties. It yields that those two factors are the major contributors to the uncertainty ΔM_spectral. This approach also accounts for detector noise, if it caused the fluctuations. The impact on M_detector, however, is small. It was deduced from Ahn and Fessler (2003), that ΔM_detector is primarily caused by the uncertainty of the standard deviation estimator for a small number of samples. In our case, the sample number equals the number of detector pixels (see Fig. 2), which is n=30 for the NIR and n=48 for the SWIR measurement. The uncertainties presented also reflect the fact that the SFA is not constant over a wavelength range of a few nanometers, as illustrated in Fig. 5. It shows plots of the SFA scaling, with the wavelength of a CO2I-like instrument in the NIR and SWIR band, calculated with the prediction model. The linear dependency only holds under the assumption of a constant dispersion, which is usually not the case for a realistic instrument. However, this result gives a clear indication of the general scaling of the speckle effect.

Figure 5 Scaling of the SFA with the wavelength of a CO2I-like instrument in the (a) NIR and (b) SWIR band, using the prediction model.

6 Conclusions

We demonstrated a comprehensive and numerical approach to quantify diffuser-induced spectral features during the solar calibration of space imaging spectrometers, which is based on established speckle theory concepts. We compared our prediction results with our current measuring method and observed a good agreement. We also gave an indication regarding the wavelength dependency of the effect. 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 point sources distributed over a 0.5^∘ angle, needs to be taken into account together with the averaging due to 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 high-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.