THE FIRST STRAY LIGHT CORRECTED EXTREME-ULTRAVIOLET IMAGES OF SOLAR CORONAL HOLES

One of the longest standing puzzles in astrophysics concerns the processes that energize the solar wind. The fast solar wind, which has a speed of ∼800 km s⁻¹ near the Earth, emanates from extreme-ultraviolet (EUV) and X-ray faint regions in the Sun's atmosphere called coronal holes (Krieger et al. 1973). Physics-based modeling efforts to identify the processes that heat and accelerate the solar wind have met with very limited success (Cranmer 2010). It is likely that much more observational input will be required to adequately constrain the numerical models used to investigate this question. NASA's EUV imaging instruments (SOHO/EIT, TRACE, STEREO/EUVI, and SDO/AIA) provide a comprehensive data set covering 1(1/2) solar cycles that can provide powerful constraints on temperatures and densities in coronal holes (Frazin et al. 2009). However, coronal holes and other faint structures are severely contaminated with stray light, which must be corrected before these vast data sets may be utilized for this purpose.

All EUV imaging instruments are afflicted to some degree by long-range scattering, which distributes a haze of stray light over the whole imaging plane. This haze is not noticeable in the brighter areas of the image, but can completely overwhelm the emissions of faint regions. Stray light arises from two sources: non-specular reflection from microrough mirror surfaces, and diffraction due to the pupil and a mesh obstructing the pupil (Howard et al. 2008). Here we describe a method for, and results of, correction of stray light in STEREO-B/EUVI, henceforth abbreviated as EUVI-B.

Stray light correction requires determination of the instrument's point-spread function (PSF). The four EUVI channels (171, 195, 284, and 304 Å) have different optical paths and therefore different PSFs. Except for pixel-scale variations due to optical aberration (whereas the scattering of interest here is significant over ranges of hundreds of pixels), the EUVI PSFs are spatially invariant (Howard et al. 2008). Thus, stray light contamination is modeled by convolution with the PSF, and the correction process is deconvolution (Starck et al. 2002).

The instrument PSFs are difficult to characterize experimentally due to the lack of a sufficiently strong EUV source, so they must be determined primarily from in-flight observations. Solar flare images can provide significant information about the entrance aperture diffraction (Gburek et al. 2006), as diffraction orders are easily visible around the flare. Unfortunately, most of the scatter is too diffuse to be observed clearly around flares. The best information about diffuse scatter is obtained from transiting bodies that do not emit in the EUV, so any apparent emission is instrumental in origin. DeForest et al. (2009) used a Venus transit to estimate the mirror scattering in the TRACE instrument by fitting a truncated Lorentzian. Here, we present the first self-consistent determination of the mirror scattering via blind deconvolution in which both the true solar image and various PSF parameters are taken to be unknown simultaneously. The information making this effort successful comes from calibration rolls and the 2007 February 25 STEREO-B lunar transit, each exposure of which provided about 50,000 lunar disk pixels illuminated only by instrumental effects.

We characterize scattering in a given EUVI filter band with a three-component PSF. The first two components, h^p and h^g, account for diffraction through the pupil and pupil mesh, and were determined analytically using Fraunhofer diffraction theory (Goodman 1996). The third component, h^m, accounts for the remaining diffuse scatter, much of which derives from the mirror microroughness. The formula for this component contains free parameters φ which we determine empirically by blind deconvolution, and we write h^m_φ to indicate the dependence. The total PSF h_φ is the convolution of these three components:

$$\begin{equation} h_\varphi = h^g * h^p * h_\varphi ^m, \end{equation} \tag{ 1 }$$

where

$$\begin{equation} (u * v)(x) = \sum _{x^{\prime } \in I} u(x - x^{\prime }) v(x^{\prime }) \end{equation} \tag{ 2 }$$

denotes the convolution of two discrete functions u and v over the index set I of the 2048 × 2048 array of CCD pixels. (We set u(x − x') = 0 at pixels x − x'∉I.) This model assumes that each component is independent, neglecting phase correlations that arise as light propagates from the entrance aperture to the primary mirror. To determine the proportionality constants of the three PSFs, we assume that each sums to unity.

We used the EUV mirror literature to help us choose an appropriate parametric model for the empirical PSF component h^m_φ. In a typical EUV scatter model, a significant fraction of the light is not scattered, while the rest is scattered by broad wings (Krautschik et al. 2002). Up to scaling and normalization constants, these wings are described by the power spectral density (PSD) of the mirror surface height function, and this PSD has been directly measured for mirrors similar to EUVI's (Martínez-Galarce et al. 2010). A log–log plot of the measured PSD versus spatial frequency is roughly piecewise linear, implying that h^m_φ is a piecewise power law whose exponent depends on pixel distance r from the origin.

Accordingly, a parametric formula for a family of piecewise power laws was used to describe h^m_φ. A series of breakpoints 1 = r₀ < r₁ < r₂ < ⋅⋅⋅ < r_b < ∞ was chosen, and the number of breakpoints, b = 8, was selected according to a χ² goodness of fit criterion (we performed several fits using different values of b, and the fit did not improve significantly for b > 8). On each subinterval [r_{i − 1}, r_i), the formula is given by $p_{\alpha,\boldsymbol{\beta }}(r) \propto r^{-\beta _i}$, where β_i ⩾ 0, and $\boldsymbol{\beta }= (\beta _1,\ldots,\beta _b)$. The parameter α represents the fraction of non-scattered light: $p_{\alpha,\boldsymbol{\beta }}(0) = \alpha$. This radially symmetric profile did not adequately describe the anisotropic scatter observed in the transit and calibration roll images, even after removal of the scatter due to h^p and h^g. We therefore generalized the formula to allow h^m_φ to have an elliptical cross section. Let M_{s, θ} denote the 2 × 2 matrix that dilates the plane by a factor of s along a line rotated θ radians counterclockwise from the horizontal axis, then

$$\begin{equation} h^m_\varphi (x) = p_{\alpha,\boldsymbol{\beta }}(\Vert M_{s,\theta } \cdot x \Vert _2+\delta), \end{equation} \tag{ 3 }$$

where the constant δ = 1 was added to avoid the power law's singularity at the origin. The free parameters of the PSF model h_φ are then $\varphi = (\boldsymbol{\beta },\alpha,\theta,s) \in \mathbb {R}^{p}$, where p = b + 3 = 11.

To determine the PSF from the lunar transit data, we described the observation process by a statistical image formation model and sought a maximum likelihood estimate of the PSF under this model. Let u^true denote the ideal image, that is, u^true(x) is the expected solar photon count that would be measured at pixel x by an ideal instrument. Due to scatter and the Poisson photon arrival process, the actual number of photon arrivals is a Poisson random variable with expected value h^true*u^true. The difference between the expected and the actual photon count is the photon noise n_phot. The observed photon count f deviates from h^true*u^true + n_phot due to CCD dark current and read noise n_ccd. Based on histograms of dark images, we find n_ccd is reasonably modeled as a Gaussian white noise process with standard deviation σ_ccd ≈ 1 digital number (DN). The variance must be divided by the photometric gain factor (the recorded DN per incident photon) to obtain the CCD noise level in units of photons. Combining n_phot and n_ccd into a single variable n, we obtain the statistical image formation model

$$\begin{equation} f = h^{{\rm true}} * u^{{\rm true}} + n. \end{equation} \tag{ 4 }$$

We would like to find a maximum likelihood estimate of the PSF under the model of Equation (4), but the Poisson–Gaussian distribution of n leads to a difficult large-scale nonlinear optimization problem. We therefore apply a variance stabilizing transform (VST) to make n approximately standard normal (mean zero and variance one), which leads to an easier least-squares problem. A VST for a Poisson–Gaussian random variable where the Gaussian has variance σ² is derived in Murtagh et al. (1995). Their full formula allows for CCD bias and non-unity camera gain; however, we correct for these in pre-processing, so the formula simplifies to ${\rm vst}(I) = 2\sqrt{I + {3}/{8} + \sigma ^2}.$ If the VST is applied to both sides of Equation (4) with σ² = σ²_ccd, we obtain

$$\begin{equation} {\rm vst}(f) \approx {\rm vst}(h^{{\rm true}} * u^{{\rm true}}) + n_{{\rm vst}}, \end{equation} \tag{ 5 }$$

where, for each x, n_vst(x) is roughly standard normal (Murtagh et al. 1995). The approximation begins to break down when I ≲ 5 photons, but this only occurs far from the solar disk in the lunar transit images. In these low-intensity regions we replace the observed image pixel value with a local average of sufficient size that the total photon counts exceed 5 in the averaged neighborhood.

We found an estimate h_φ of the EUVI-B PSF h^true by solving a blind deconvolution problem on a series of eight images f₁, ..., f₈ from the lunar transit series. We assume the scatter-free images u^true_i are positive everywhere and zero on the lunar disk pixels Z_i, leading to the constraints

$$\begin{equation} u_i^{{\rm true}} \ge 0, \quad u_i^{{\rm true}}(Z_i) = 0, \quad i=1,\ldots,8. \end{equation} \tag{ 6 }$$

The lunar disk was identified by detecting its edge pixels with gradient thresholding, then fitting a circle to the detected edge pixels. An approximate maximum likelihood estimate for the PSF h_φ and images {u_i} is obtained by solving

$$\begin{eqnarray} \displaystyle\mathop{\rm minimize \quad}_{\varphi = (\boldsymbol{\beta},\alpha,s,\theta), \{ u_i\}} &&\sum _{i=1}^8 \Vert {\rm vst}(h_\varphi * u_i) - {\rm vst}(f_i)\Vert ^2 \nonumber\\ {\rm subject\ to} &&\quad u_i \ge 0, \, \, u_i(Z_i) = 0 \, \, {\rm for\ all}\ i,\nonumber \\ &&\boldsymbol{\beta }\ge 0, \, 0 \le \alpha \le 1. \end{eqnarray} \tag{ 7 }$$

Since each u_i is a 2048 × 2048 image, this problem has over 32 million variables and is intractable by general-purpose numerical methods. We solved this problem with a customized variant of the variable projection method (Golub & Pereyra 2003), which is efficient for problems of this kind.

Figure 1 shows a log–log plot of the 171 Å profile function obtained by blind deconvolution, along with views of the total PSF and its core. Note that the profile is constant for 200 < r < 1200 pixels, implying considerably more long-range scatter than a constant power-law exponent would predict. The total PSF is dominated by the elliptical power-law decay, although the diffraction orders from the pupil and mesh PSFs are visible near the PSF core.

Zoom In Zoom Out Reset image size

Once we determined the PSF from blind deconvolution, we were able to correct any observed EUVI-B image f via more standard deconvolution methods. Let $\mathcal {C}(h) : \mathbb {R}^I \rightarrow \mathbb {R}^I$ denote the linear operator that convolves an input image u with a fixed PSF h: $\mathcal {C}(h)u = h*u$. The stray light correction of a given image f is obtained by applying the inverse of $\mathcal {C}(h)$ to f:

$$\begin{equation} u = \mathcal {C}(h)^{-1}f. \end{equation} \tag{ 8 }$$

The inverse operator does not exist for arbitrary h. However, our PSFs have an origin value h(0) > 1/2, which makes $\mathcal {C}(h)$ diagonally dominant and invertible. The inversion of Equation (8) can be performed with the conjugate gradient method in less than a minute on a laptop.

A first test of our stray light correction was performed by applying cross-validation (CV) to the lunar transit series. For each of the eight images f_i, we determined a PSF h_i by solving Equation (7) with f_i removed from the data set. We then calculated the stray light corrected CV image $u_i^\star = \mathcal {C}(h_i)^{-1} f_i$. Note that both h_i and u^⋆_i are calculated without any assumption of a zero lunar disk in image i, so the values of u^⋆_i on the lunar disk represent an independent check on the correction's effectiveness. In Figure 2, we compare the lunar disks of f_i and u^⋆_i in 171 Å. We find that the lunar disk values after deconvolution are very strongly clustered near zero, as seen in a histogram of the lunar disk values (bottom right).

Zoom In Zoom Out Reset image size

Calibration roll images provide direct evidence of anisotropic scatter and our deconvolution's ability to correct it. We will present the analysis for 171 Å here, leaving the other bands to future publications. On 2011 November 8, STEREO-B executed a 360° calibration roll, and in each band, EUVI-B acquired nine solar images at roll angles of 0°, 60°, 90°, 120°, 180°, 240°, 270°, 300°, and 360° relative to the pre-roll position. These images are useful because the direction of anisotropic scatter rotates with STEREO-B, introducing discrepancies between the images. To pick out these discrepancies, an 8 × 8 boxcar was applied to reduce noise and the images were rotated into a common solar coordinate system. This procedure was repeated on stray light corrected roll images, resulting in series {f_ρ}³⁶⁰_{ρ = 0} and {u_ρ}³⁶⁰_{ρ = 0} indexed by roll angle ρ. A haze of scattered light can be observed off the limb in the image f₀ (Figure 4, top left), which is much reduced in the corrected image u₀ (middle left).

Anisotropic scatter rotation was tracked using the difference images Δf_ρ = f_ρ − f₀ and Δu_ρ = u_ρ − u₀. To avoid analyzing regions of the Sun that changed significantly over the course of the roll, we examined the difference Δf₃₆₀ between the pre- and post-roll images, and all pixels where |Δf₃₆₀| > 1  DN were masked out of the difference images. In Δf₉₀ (top right), we observe a "dark axis" and a ≈90°-rotated "light axis" corresponding to the preferential scattering directions for f₀ and f₉₀, respectively. (These regions are very diffuse, so the separation angle can only be estimated up to about ±10°.) These axes are eliminated in Δu₉₀ because the stray light correction greatly reduces the anisotropic haze. Similar reductions are observed at other ρ values.

To estimate the error in the corrected images, we decompose the total error ε = u − u^true into components due to noise and PSF error, and estimate their contributions separately. To obtain the decomposition we substitute Equation (4) into Equation (8) and introduce the PSF error variable δh = h − h^true, obtaining

$$\begin{equation} \varepsilon = u-u^{{\rm true}} =\mathcal {C}(h)^{-1} \mathcal {C}(\delta h) u^{{\rm true}} + \mathcal {C}(h)^{-1} n. \end{equation} \tag{ 9 }$$

The first term on the right-hand side gives the error due to PSF misestimation, which we name ε_psf. The second, ε_noise, is the noise in the corrected image.

To estimate ε_psf, we first apply an 8 × 8 moving average filter to u and f to reduce noise, so that u − u^true ≈ ε_psf. (Averaging also reduces ε_psf on small scales, potentially leading to underestimation of the error. However, we are only interested in the large-scale effect of stray light, and on large scales, ε_psf is unaffected by averaging.) Next, we assume that ε_psf can be bounded by a quantity σ_psf which is proportional to the magnitude of the stray light correction |u − f|:

$$\begin{equation} |u - u^{{\rm true}}| \le B|u - f| = \sigma _{{\rm psf}}. \end{equation} \tag{ 10 }$$

We use the bounding term B|u − f| because it shares two distinctive properties with ε_psf: both are approximately invariant if u^true is changed by scaling or addition of a constant. These invariance properties help ensure that an empirical bound based on lunar transit images alone will continue to be valid for other solar images.

We determined the constant B by requiring that Equation (10) holds true over almost all of the lunar disk pixels of f_i and u^⋆_i. Setting u = u^⋆_i, f = f_i, dividing both sides by |u^⋆_i − f_i|, and noting that u^true_i = 0 on the lunar disk, we find |u^⋆_i|/|u_i^⋆ − f_i| ⩽ B is required on the lunar disk pixels. We call the left-hand side ratio b^⋆_i and collect the values of b^⋆_i on the eight lunar disks into a histogram. A normalized histogram for 171 Å is given in Figure 3, right, and values at the 68th, 95th, and 99.7th percentiles, corresponding to the first three standard deviations of a Gaussian, are 0.08, 0.13, and 0.16. We set B equal to the 95th percentile of the histogram, which gives a conservative 2σ error bound robust to outliers. To estimate the contribution of noise to the error, we analytically calculate its covariance matrix using the known variance of n. The diagonal of this covariance matrix, σ²_noise, is our estimate of the squared error due to noise. We add the two independent errors in quadrature to obtain the final error estimate σ: σ² = σ²_psf + σ_noise².

Zoom In Zoom Out Reset image size

Here we note just a few effects of stray light correction on faint regions. Figure 5 shows a coronal hole before and after correction. The stray light corrected coronal hole is significantly dimmer: the percent change (u − f)/f ranges from −40% to −70% over most of the coronal hole, implying that 40%–70% of the apparent coronal hole emissions are stray light. The plot of intensity versus position (right) contains the error bars σ. Note that they are small relative to the size of the correction, which is representative of typical results. Stray light correction also removes the haze seen above the limb (Figure 4, top and center left), often reducing the intensity by 75% or more (bottom left). A filament just above the solar center is also reduced in intensity by 50%–70% after stray light correction (bottom left).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These large downward corrections to faint regions have major impacts on the plasma diagnostics available from EUV images, which in turn are related to the electron temperature T_e and density n_e. Assuming electron impact excitation, the intensity of the emission in images from the kth spectral band is I_k∝∫_LOSdl n²_e(l)Q_k(T_e(l)), where Q_k(T_e) comes from an optically thin plasma emission model and ∫_LOSdl represents integration along the line of sight (Brown et al. 1991). This relationship may be utilized to create three-dimensional tomographic maps of n_e and T_e (Frazin et al. 2009). A downward correction of the observed intensity causes a proportional reduction of the estimated value of n²_e; indeed, analysis of the 171, 195, and 284 Å intensities shows a downward revision of the coronal hole column density [∫_LOSdl n²_e(l)]^1/2 of ∼40%. The effect on the estimated T_e is more complex. The removal of stray light from the off-limb causes a dramatic steepening of the profile function n_e(h) (where h is the height above the photosphere). The specific impact of stray light correction, including constraints on solar wind models from the corrected profiles n_e(h), T_e(h), is currently under investigation. Obvious consequences include reduction of the plasma β, electron–ion collision rates, and the mass of solar wind plasma requiring acceleration.

We have obtained PSFs for all four bands of STEREO-B/EUVI using lunar transit data, which enable us to correct all EUVI-B images for stray light and provide uncertainties. Similar methods may be applied to treat the stray light problems in the other solar EUV imaging instruments (SOHO/EIT, TRACE, STEREO-A/EUVI, SDO/AIA). This work and its heliophysical implications will be reported in more detail in upcoming publications.

We thank Jean-Pierre Wuelser for continuous assistance with EUVI technical issues; Frederic Auchere and Raymond Mercier, for useful discussions of PSF modeling; and Simon Plunkett, for granting our request for EUVI calibration roll data.