SYSTEMATIC DE-SATURATION OF IMAGES FROM THE ATMOSPHERIC IMAGING ASSEMBLY IN THE SOLAR DYNAMICS OBSERVATORY

Since its launch on 2010 February 11, the Atmospheric Imaging Assembly (AIA) on board the Solar Dynamics Observatory (SDO; Lemen et al. 2012) has been contributing to our knowledge of solar flare emission, permitting full-disk imaging in 10 wavebands (7 of which are centered at EUV wavelengths), at a spatial resolution of about 1 arcsec and with a temporal cadence of 12 s in each of the EUV bands. However, the scientific potential of AIA/SDO images is limited by the presence of significant saturation, which flattens and muddles the brightest, physically most interesting core of such images for even relatively modest flares.

Saturation refers to the condition where, in a CCD-based imaging system, pixels lose their ability to accommodate additional charge (primary saturation), causing charge to spill into their neighbors. This secondary effect, named blooming, also induces bright rays typically along the north–south (for AIA) axis in the image.

We now introduce a computational method that utilizes the diffraction component of the point-spread function (PSF) of AIA/SDO (Boerner et al. 2012; Poduval et al. 2013; Grigis et al. 2011) to recover the signal in the primary saturation region. The PSF of each AIA telescope models the system response to a point source and describes blurring, dispersion, and diffraction. The various AIA PSFs depend on the wavelength of the incoming photons and present a two-component structure, made of a diffraction component and a diffusion component, that is similar to those of the TRACE passbands (Gburek et al. 2006). The diffraction component arises from the two square wire grid meshes that support thin EUV filters and that are offset from each other by an angle of 10°, giving rise to four diffraction patterns that appear as eight rays of fringes that cross at the source of the emission. During the observation of a flaring event, the diffraction component of the AIA PSF creates copies of the bright flare image that would be otherwise obscured by saturation at its point of origin in the CCD. The goal of this Letter is to illustrate the use of automatic image processing tools to exploit the diffraction pattern showing up outside the core of the saturated image to recover the true image inside the saturated region. The feasibility of this approach was demonstrated by previous analyses of TRACE and AIA diffraction patterns (Raftery et al. 2011; Lin et al. 2001; Grigis et al. 2011) where, however, the inverse diffraction problem is addressed by means of non-automatic, laborious, and semi-heuristic procedures. Our approach is based on three image-processing steps. In the first step (correlation), we identify the region of primary saturation, distinguishing it from blooming on the basis of a data-dependent correlation procedure. In the second step (convolution), the diffraction fringes in the image are localized by convolving the diffraction PSF with the pixels in the region of primary saturation. The third step (back-projection) is used to reconstruct the true image in the primary saturation region using the information contained in the diffracted images. Details about these steps are given in Section 2. Section 3 illustrates the performances of the method in the case of AIA/SDO data recorded during the 2013 July 8 event. Our conclusions are offered in Section 4.

The general model equation that mimicks the signal formation process in an optical system is

$$\begin{equation} g = K * f + b, \end{equation} \tag{ 1 }$$

where g is the measured data, K is the global PSF, f is the true flux distribution that has to be reconstructed from the data, and b is the background. This equation applies to all points in the data image. However, in the case of the AIA/SDO image reconstruction problem, the use of this equation should account for the fact that:

• 1.  
  the unknown true image can be recovered only in a subdomain of the field of view (FOV) corresponding to the primary saturation region;
• 2.  
  since information on the primary saturation region is contained in the diffraction fringes, the usable data, g, are just those pixels containing the diffraction pattern;
• 3.  
  only the diffraction component of the PSF carries the information that can be used to reconstruct the true image, and therefore we are free to neglect the remaining components of the PSF.

In fact, the PSF associated with each AIA/SDO wavelength can be realistically modeled as

$$\begin{equation} K(x,y) = K_C(x,y) + K_D(x,y), \end{equation} \tag{ 2 }$$

where K_C(x, y) is the core, Gaussian-like component mimicking diffusion, and K_D(x, y) is the diffraction component wherein about 20% of the total flux is coherently scattered due to the two square wire grid meshes supporting the thin EUV filters. Figure 1 shows the PSFs corresponding to the 131 Å and 304 Å passbands. In order to make them more visible in the same figure, the 304 Å PSF has been rotated by 45° and a threshold has been applied to better show the diffraction patterns over the full range of the FOV. The figure points out the presence of eight arms of diffraction for each passband, how diffraction scales with wavelength, and how dispersion in the diffraction peaks scales with respect to dispersion in the core.

Zoom In Zoom Out Reset image size

The diffraction component of the PSF plays a crucial role in the de-saturation process described in this section. The workflow of such a process is illustrated in Figure 2 for a flare occurring on 2013 July 8, measured at 131 Å (more details on this event will be given in Section 3). The computational steps at the basis of its realization are as follows.

• Correlation. This step separates the overall saturation domain into the primary saturation region, where the true image can be reconstructed, and the blooming region, which is extraneous to the reconstruction process. Let us denote with S the overall saturation domain in the AIA image in the top left panel of Figure 2. S is the sum of the primary saturation region PS and the blooming region B and all its pixels have intensities of around 16, 000 DN pixel⁻¹. Such pixels are masked away to obtain the image g' represented in the top right panel of Figure 2. The back-projection matrix $K_D^T$ is the transpose of the diffraction matrix K_D and has a size equal to the number of unsaturated pixels times the number of saturated pixels. We used a slightly modified version of the routine "aia_calc_psf.pro" in Solar SoftWare (SSW) for the generation of the PSF (Grigis et al. 2011) and extracted the diffraction component by masking out the central peak. We then computed the correlation map

  $$\begin{equation} C = K_D^T g^{\prime } \end{equation} \tag{ 3 }$$

  for all pixels in S, and show it in the middle left panel of Figure 2. Since the primary saturation is identified by those pixels in S for which C is high, a thresholding procedure allowed the identification of PS in the correlation map.
• Convolution. This step computationally finds in the AIA frame the location of the diffraction pattern that can be used to infer information on the true image. To do this, in the middle right panel of Figure 2, we highlighted with a white boundary the domain PS in the observed saturated image where information can be recovered by means of image reconstruction. We computed the convolution between PS and K_D to localize the diffraction fringe in g' and utilized just that information, indicated with g in the reconstruction procedure. In fact, g contains pixels where the diffraction effect is most significant.
• Back-projection. This step implements the recovery of information f in the primary saturation region using the diffraction fringes g as input data. Since we do not have at our disposal the inverse matrix of K_D, we use the back-projection matrix $K^T_D$ within the expectation maximization (EM) iterative scheme (Dempster et al. 1977),

  $$\begin{equation} f^{k+1} = \frac{f^{k}}{K_D^T I} K^T_D \left(\frac{g}{K_Df^k + b}\right), \end{equation} \tag{ 4 }$$

  to find a value of f that is consistent with g. Here, f^k denote successive approximations of the unknown spatial distribution of the signal in the primary saturation region PS, $K_D^T$ has size equal to the number of pixels in the diffraction fringes times the number of pixels in the primary saturation region PS, I is the unit vector, and b is an estimate of the background. EM (also known as the Lucy-Richardson method in this context; Lucy 1974) implements the maximum likelihood in the case of Poisson data and imposes positivity on the solution space (CCD data are more correctly quasi-Poisson, but EM effectiveness relies more on positivity than on the way the data-prediction discrepancy is computed). In our implementation, the stopping rule that avoids numerical instability is based on the concept of a constrained residual (Benvenuto & Piana 2014) and imposes that the empirical value of this random variable coincides with its expectation value (Benvenuto et al. 2013).

The result of the back-projection step is the de-saturated image in the bottom left panel of Figure 2 where the blooming effects are still present, the diffraction pattern is clearly visible, and the recovery of information in the region originally affected by primary saturation is evident. Such a recovery of information is even more evident in the bottom right panel of the figure, a zoom on the de-saturated domain. We note that the use of the log scale better shows the presence of the diffraction fringes with respect to the reconstructed signal.

Zoom In Zoom Out Reset image size

On 2013 July 8, the Sun experienced a C9.7 GOES class flare in the time range 01:13–01:23 UT with peak at around 01:22 UT. This event was very well observed by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI; Lin et al. 2002) but, despite the moderate peak flux, the corresponding AIA images saturated at most wavelengths and times. In principle, simultaneous solar flare observations with SDO and RHESSI provide spatially resolved information on hot plasma and energetic particles during flares, and the unavailability of EUV maps due to saturation clearly limits this kind of analysis. As a test of our de-saturation approach, we considered this event and, specifically, we applied the method to a number of saturated images collected at three different wavelengths (131 Å, 211 Å, and 304 Å), two times for each wavelength. We downloaded AIA Level 1.5 data from http://www.lmsal.com/get_aia_data and again used "aia_calc_psf.pro" to generate the PSF. The results of our analysis are given in Table 1, which reports the Cash statistic (Cash 1979) values predicted by the de-saturated images. These values, also named C-statistic values, measure the distance between the observation g and the data predicted by the reconstructed signal f throughout the diffraction PSF K_D (the way C-statistic is computed accounts for the fact that the experimental observations are affected by Poisson noise). We have computed the C-statistic of a set made of 100 random samples of the fringes, repeated this process 10 times (every time changing the set of random samples), and reported the average values over these 10 times. We point out that for values ⩾10, the C-statistic is analogous to χ², while for smaller values it denotes a reliable fit. In particular, the values close to 1 corresponding to stronger fringes indicate a good agreement between the data and the model. Finally, Figure 3 visually compares saturated and de-saturated maps corresponding to one of the two observation times for each wavelength. We point out that the linear scale softens the effect of diffraction fringes and blooming (blooming cannot be eliminated by means of reconstruction, but its visual effect could be further reduced by means of interpolation).

Zoom In Zoom Out Reset image size

The potential of this method for a more extensive use of EUV data is illustrated in the example of Figure 4. To construct it, we first considered a RHESSI image of hard X-ray flux reconstructed by Pixon (Puetter & Gosnell 2006) in the time interval 01:22:14–01:22:26 UT for the energy channel 12–25 keV. The contours of this X-ray image are superimposed on the de-saturated 131 Å image corresponding to the AIA/SDO observation that occurred at 01:22:20 UT (an offset of around 2 arcsec is visible between the RHESSI and AIA images).

Zoom In Zoom Out Reset image size

The EUV analysis of flaring events by means of AIA/SDO is notably limited by saturation and blooming, occurring at all wavelengths and times even in the case of moderate peak flux. This Letter illustrates a method that addresses, for the first time in a systematic way, the problem of eliminating primary saturation from these maps, utilizing standard image processing techniques like correlation, convolution, and statistical deconvolution (the problems of recovering information in the blooming region and removing the multiple diffraction components will be addressed in future research). This method strongly relies on the characteristics of the applied PSF; here, we used a synthetic form generated by a routine in SSW, but it would be interesting to run the reconstruction algorithm also in the case of the semi-empirical PSF described by Poduval et al. (2013). We applied our approach to data recorded during the 2013 July 8 flare at three different EUV wavelengths and showed its effectiveness to reduce saturation in a variety of artifact configurations and intensities. Once fully automated, these de-saturation procedures will be included in Solar SoftWare. Further, we will be extending this technique to include the wavelength-dependent dispersion of the PSF to integrate the method of Raftery et al. (2011) into this schema.

The impact of this method on the exploitation of EUV information in solar flares is notable, for example, in imaging spectroscopy studies of the dynamic response of the solar corona and chromosphere to solar flares. Another broad and significant applicability domain involves multispectral studies. For example, combined analysis of AIA/SDO and RHESSI data may provide information on crucial plasma parameters like differential emission measure, temperature, and density, which would notably contribute to unveiling still mysterious processes concerned with energy transport mechanisms in flares.

R.S. was supported by NASA grant NNX14AG06G. G.T. was supported during his stay at Goddard by the Catholic University of America in Washington, DC, using RHESSI funds under grant number NNX11AB37G. Anna Maria Massone and Federico Benvenuto are kindly acknowledged for fruitful discussions.