Compressed sensing and Sequential Monte Carlo for solar hard X-ray imaging

TL;DR

This paper introduces two inversion methods combining compressed sensing and Sequential Monte Carlo techniques for reconstructing solar hard X-ray images, validated on real and synthetic data from RHESSI and STIX instruments.

Contribution

The paper presents novel inversion algorithms that integrate compressed sensing with Sequential Monte Carlo methods for improved solar X-ray imaging.

Findings

Methods successfully reconstruct solar X-ray images from RHESSI data.

Algorithms perform well on synthetic STIX data.

Enhanced image quality compared to traditional techniques.

Abstract

We describe two inversion methods for the reconstruction of hard X-ray solar images. The methods are tested against experimental visibilities recorded by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) and synthetic visibilities based on the design of the Spectrometer/Telescope for Imaging X-rays (STIX).

Figures

Equations

$$HFx=V,$$

$$\hat{x}=\underset{x}{\min }\{\Vert HFx-V{\Vert }_2^2+\lambda \Vert Wx{\Vert }_1\},$$

$$p(x|y)=\frac{p(y|x)p(x)}{p(y)},$$

Full text

1 Introduction

The NASA Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) [8] and the ESA Spectrometer/Telescope for Imaging X-rays (STIX) [3] are two space telescopes for imaging hard X-rays that rely on rather similar imaging technologies. RHESSI has been decommissioned on August 16 2018 after more than $16$ years of successful operations, while STIX is going to fly in the next two years. Both hardwares allow the modulation of the X-ray flux coming from the Sun, providing as a result sparse samples of its Fourier transform, named visibilities, picked up at specific $(u,v)$ points of the Fourier plane. Therefore, for both RHESSI and STIX, image reconstruction is needed to determine the actual spatial photon flux distribution from the few Fourier components acquired by the hard X-ray collimators [1] [2] [5] [7] [9] [10]. In Section 2 of this paper we briefly overview a reconstruction method based on compressed sensing [6]. In Section 3 we provide more insights on a Monte Carlo method for the Bayesian estimation of several imaging parameters [11] [12]. Our conclusions are offered in Section 4.

2 Compressed sensing for hard X-ray image reconstruction

Figure 1 shows how RHESSI and STIX grids sample the $(u,v)$ frequency domain. From this design, it follows that the mathematical model for data formation in the framework of these two instruments is, in a matrix form:

[TABLE]

where $x$ is the photon flux image to reconstruct, $V$ are the experimental visibilities, $F$ is the discretized Fourier transform, $H$ is a mask that realizes the sampling in the $(u,v)$ plane. The reconstruction of $x$ from $V$ is an ill-posed problem and therefore regularization is required to mitigate the numerical instabilities induced by the observation noise. A possible approach is to apply an $l_{1}$ penalty term in some transformation domain. This can be realized by solving the minimum problem [6]

[TABLE]

where the regularization term $\|Wx\|_{1}$ is designed to penalize reconstructions that would not exhibit the sparsity property with respect to the Finite Isotropic Wavelet Transform [6]. Figure 2 compares the reconstructions provided by this compressed sensing algorithm to the ones obtained by using other four visibility-based imaging methods currently implemented in the RHESSI pipeline [5] [7] [9].

3 Sequential Monte Carlo for hard X-ray image reconstruction

Sequential Monte Carlo (SMC) samplers are computational methods aiming at sampling target distributions of interest, and are often applied to sample the posterior distribution $p(x|y)$ as given by Bayes’ theorem

[TABLE]

where $x$ is the unknown, $y$ is the observation, $p(x)$ is the prior probability encoding all a priori information, $p(y|x)$ is the likelihood encoding the image formation model (1) and the noise model, and the marginal likelihood $p(y)$ is a normalization factor. In the case of RHESSI and STIX imaging, $x$ is the image to reconstruct and $y$ denotes the set of recorded visibilities. We modeled $x$ as $x(N,T_{1:N},\Theta_{1:N})$ where $N$ is the number of sources in the image, $T_{1:N}=(T_{1},\ldots,N)$ represents the source types (Gaussian, elliptical, loop-like) and $\Theta_{1:N}=(\theta_{1},\ldots,\theta_{N})$ contains the parameters characterizing each source. We chose a prior distribution factorized as the product of a Poisson distribution for $N$, uniform distributions for the source types and uniform distributions for the source parameters [11] [12]. Sequential Monte Carlo [4] computes the posterior distribution iteratively, by constructing a sequence of converging approximate distributions. Once the posterior is determined, it can be used to compute the solution image and all image parameters. Figures 3 and 4 show results provided by this approach using simulated STIX visibilities and experimental RHESSI visibilities, respectively.

4 Conclusions

This paper shows the performances of two image reconstruction methods formulated for hard X-ray solar visibilities. The implementation of the corresponding tools within Solar SoftWare (SSW) is under construction.