Introduction

Sufficient spatial resolution is required in experimental studies of high-temperature plasma physics including the fusion research and development. However, determining the spatial behaviour is not straightforward due to the very poor accessibility of fusion relevant plasmas. Active diagnostic systems which determine plasma properties from its interaction with a physical probe (e.g. laser light, particle beam) offer an excellent spatial resolution in a very specific region of plasma only. Majority of the passive diagnostic systems, which can diagnose plasma from outside its border, allow for spatial resolution of the integral plasma emission, i.e. of the plasma projections. The aim of plasma tomography is to determine local plasma properties from the measured projections.

Spatial resolution of the reconstructed plasma image is, however, rather poor due to the sparse data caused by the limited accessibility, with detector positions typically confined to ports of the vacuum vessel. Implementation of additional constraints (e.g. non-negativity of the emission, zero border emission) and a-priori information (e.g. increased smoothess along field lines) can improve the resolution considerably. On the other hand, temporal resolution of the plasma tomography can be sufficiently high for studies of rapid emissivity evolution [1,2,3].

Plasma tomography has been applied on several different diagnostic systems. It is rather widespread with Soft X-ray (SXR) diagnostics, often using linear pinhole cameras, see e.g. [2, 4, 5]. Due to the high temporal resolution, the SXR tomography can contribute to the MHD analyses [2] and to the impurity transport studies [6]. At several facilities, data from bolometric diagnostic systems are regularly analysed using tomography, e.g. at the Joint European Torus JET [1]. In the perspective of fusion reactors the key application of tomography is linked to analyses of data from the neutron cameras, see [7, 8], where it is expected to contribute substantially to regular monitoring of the fusion power.

In this paper, our recent experience with the plasma tomography development for fusion research (in particular at JET) is summarised. In “Current Methods of Plasma Tomography” section the main methods currently applied in plasma tomography are outlined, with a focus on the Tikhonov regularisation. In “Progress in the Minimum Fisher Regularisation” section, the Minimum Fisher Regularisation method is briefly explained and its further development into current applications, fostered by our team, is presented. In particular, a novel approach to the combined neutron tomography with spectral unfolding is proposed. The key recommendations for tomography applications and development in fusion reseach are recapitulated in the conclusions.

Current Methods of Plasma Tomography

The aim of a typical two-dimensional tomographic reconstruction is to retrieve spatial distribution of a scalar field in a plane (e.g. emissivity cross-section image) from its line integrated values (e.g. projections of the emissivity, provided the object is optically thin). Tomographic reconstruction is a special case of an inversion task. Even with the complete knowledge of the line integrals (for all angles and distances) the inversion task is ill-posed. As a consequence, a minor error in the data may result in substantial errors in the reconstructed image. In order to mitigate this issue, different regularisation methods are usually applied. The regularisation process typically introduces some a-priori knowledge, e.g. some restriction on the result complexity (commonly enforcing smoothness of the result).

In plasma tomography, however, the inversion task is further complicated by the fact that the experimental constraints usually prevent measurements in arbitrary directions. Therefore, many plasma projections are missed, i.e. data are sparse, making the tomographic inversion underdetermined for any practical spatial resolution, see Fig. 1. Semi-analytical inversion methods, e.g. based on filtered backprojection (FBP) [1], encounter diverse difficulties under these circumstances. Therefore, dedicated methods with reliable and robust regularisation principles had to be developed for plasma tomography.

Fig. 1
figure 1

The Soft X-ray tomography setup at tokamak TCV [3]. While number of the measured line integrated projections is rather high (left, with R the radius and z the vertical coordinate) the actual coverage in the projection space is still quite sparse (right, with \(\theta\) the angle and and p distance of lines from the vessel centre at \(R=0.88\, {\text {m}}, z=0 \, {\text {m}}\))

Full size image

In order to process digitised projection data, FBP or Fourier transform based approaches to tomography inversion implement digitisation of the analytical inversion, i.e. of the inverse Radon transform. However, given sparse data, digitisation of the unknown image prior to the tomography inversion proves to be better adapted to implementation of a-priori knowledge [1]. Therefore, in a large majority of current computer codes for plasma tomography, the unknown emissivity distribution \(g(\mathbf{r })\) is decomposed into basis functions \(b_j(\mathbf{r })\). The basis functions can be global, each covering the whole size of the reconstructed image (e.g. special functions) or local, decomposing the image into regions (e.g. square pixels), see [1]. Consequently, the plasma tomography task is reduced to finding the amplitudes \(g_j\):

$$\begin{aligned} g(\mathbf{r }) = \sum _j^\infty g_j b_j(\mathbf{r })\, , \end{aligned}$$
(1)

While with infinite series the image can be (in principle) reconstructed to full perfection, in practice the sum is truncated. In other words, the choice of a limited amount of basis functions \(g_j \, , \, j=1, \ldots ,N\) contribute to the regularisation of the task. In this respect, implementing any other additional information, in particular the magnetic field geometry directly into the shape of the basis functions is discouraged. Indeed, knowledge of the magnetic field geometry is burdened with errors and, consequently, any adapted shape of the basis functions enhances artefacts. Notice that this holds in particular for the notorius Abel inversion (see e.g. [9]), where the unfolding procedure stems from the unprecise knowledge of the magnetic axis position. In the current plasma tomography, it is recommended to apply standard 2D tomography with preferential smoothness along magnetic flux surfaces (see “Progress in the Minimum Fisher Regularisation” section, Eq. 11) instead of the 1D Abel inversion.

In the vast majority of diagnostic methods, the projections are determined for a finite number of chords (lines of sight, views) \(f_i \, , \, i=1, \ldots ,P\) which will be further referred to as the line integrated measurements, although in general the chords need not be strictly linear in the projection plane. Actually, the real width of the chords is often evaluated within the tomographic algorithms. A simple relation between the finite number of amplitudes \(g_j\) and the finite number of projection measurements \(f_i\) can then be proposed:

$$\begin{aligned} f_i = \sum _j^N T_{ij}g_j + \xi _i ,\quad i \in {1,\ldots ,P} , \end{aligned}$$
(2)

where \(T_{ij}\) stands for elements of the contribution matrix (sometimes referred to as the geometric matrix) and \(\xi _i\) for misfits which compensate, in a real experimental setup, both systematic errors and noise.

Different regularisation methods are applied in algorithms of tomographic codes in order to find unknown \(g_j\) from sparse data \(f_i\). A reliable performance of the algorithm, with sufficient precision, robustness against artefacts caused by \(\xi _i\), low demand on CPU etc. is usually requested from any applicable tomography method. Some of the methods will actually find a reconstruction matrix \({\mathbf {M}}\):

$$\begin{aligned} g_j = \sum _i^P M_{ji}f_i \, , \end{aligned}$$
(3)

Tikhonov regularisation is probably the most widespread method where the reconstruction matrix \({\mathbf {M}}\) is derived. In the Tikhonov regularisation

$$\begin{aligned} M_{ji} = \sum _k^N \left( \sum _l^P T^T_{kl}T_{lj}+ \lambda H_{kj}\right) ^{-1} T^{T}_{ki} \end{aligned}$$
(4)

where \(\lambda\) is the regularisation factor which sets the relative strength of the regularisation, and the regularisation matrix \({\mathbf {H}}\) is given, in the simplest case, by

$$\begin{aligned} H_{kj} = \sum _m^N B^T_{km}B_{mj} \, , \end{aligned}$$
(5)

where the smoothing matrix \({\mathbf {B}}\) corresponds to numerical differentiation, for details see e.g. [2]. The regularisation factor \(\lambda\) can be determined iteratively with the Pearson’s \(\chi ^2\) test:

$$\begin{aligned} \chi ^2 = \frac{1}{P} \sum _i^P \frac{\xi _i^2}{\sigma _i^2} \rightarrow 1 \, , \end{aligned}$$
(6)

where \(\xi _i\) is given by Eq. 2 and \(\sigma _i\) represents expected data errors (standard deviations of the projection data \(f_i\)), see e.g. [2]. Figure 2 gives a flowchart of the Tikhonov regularisation as implemented in practice. Another efficient and reliable implementation of the Tikhonov regularisation can be made using the Generalised Singular Value Decomposition (GSVD), see e.g. [10].

Fig. 2
figure 2

Flowchart of algorithm solving the Tikhonov regularisation, including the implemented value of the convergence threshold of the Pearson’s test and the typical maximum number of iteration loops

Full size image

Besides the Tikhonov regularisation, considerable efforts in plasma tomography have been invested into the development and application of reconstruction methods based on Bayesian approach, see e.g. [5, 11]. At present, the Maximum likelihood (ML) tomography represents application of a probabilistic method where systematic expertise in experimental data analyses has been collected, see [12]. It is a nonlinear iterative algorithm that attempts to find the estimate of the emissivity distribution that is most consistent with the measured tomographic projections in the sense of maximizing the likelihood (conditional probability of the data given the parameters). The emission is considered to be a Poisson process and therefore \(f_i\) presents a sample from a Poisson distribution, whose expected value is f. Consequently, the probability of obtaining the measurement \(f={f_i |i=1,\ldots ,P }\) from the emissivity \(g={g_j |j=1,\ldots ,N }\) is given by the likelihood function:

$$\begin{aligned} L(f/g) = \prod _{i} \frac{1}{f_i !} (f)^{f_i} \times \exp (f) \, \end{aligned}$$
(7)

The ML estimate is obtained by maximizing the above expression:

$$\begin{aligned} g_{\mathrm{ML}} = {\text {argmax}}_g L(f/g) \end{aligned}$$
(8)

The mathematical basis for a broadly applicable algorithm has been first applied to images by Richardson [11] and Lucy [13] but the method has been started to be used extensively in tomography only after the introduction of an iterative solution for finding the ML estimate by Shepp and Vardi [14] and, independently, by Lange and Carson [15]:

$$\begin{aligned} g_j^{(k+1)} = \frac{g_j^{(k)}}{s_j} \sum _{i} \frac{f_i}{\sum _{m} T_{im}g_m^{(k)}} \, \, T_{ij} \, \end{aligned}$$
(9)

where k indexes the iterations and \(s_j = \sum _{i} T_{ij}\) is the probability that emission originating in pixel n will be recorded in a projection bin.

As already mentioned the tomography problem is a highly undetermined inversion leading to an ill-posed mathematical problem. In order to obtain realistic and robust solutions, it is therefore strongly recommended also for the ML method to introduce a regularisation procedure which consists of imposing smoothness along the magnetic surfaces, given by the plasma equilibrium. The ML method uses 1-D average filtering on a sliding window, which moves along the magnetic contour lines. The iterative reconstruction formula 9 has a particular form, which is advantageous with respect to modelling of the projection noise propagation. In other words, it is possible to retrieve a variance image which accompanies the reconstructions. Accurate modelling of the projection noise propagation is important for both qualitative interpretation and quantitative analysis of the reconstructed images. Following the ideas first introduced by Barret [16] two approximations can be introduced for obtaining the variance image: to consider that the noise is small compared to the mean value of the reconstruction and to assume that the convergence of the ML algorithm is fast enough, so that the projection of the current estimate is close to the noise-free projection. Details of the implementation for JET tomography are given in [17].

In JET, the ML method has been applied for gamma and neutron tomography (see e.g. [12, 18] for representative examples). More recently it has been implemented also for bolometry [19], where the evaluation of the uncertainties accompanying the reconstruction is particularly important, allowing the estimation of the confidence intervals for the radiated power (Fig. 3).

Fig. 3
figure 3

Illustration of the ML bolometric tomography at JET: reconstruction (top-left), image variance (top-right) and radiation profile versus the normalised \(\varPsi\) coordinate with the estimate of the uncertainties in the emitted power (bottom)

Full size image

The tomographic inversion process can also be performed with neural networks. A first attempt at doing this involved the use of a neural network to find the parameters of a two-dimensional gaussian distribution that would best fit the measurements of a horizontal and a vertical soft X-ray camera [20]. This assumed that the plasma profile could be approximated by a 2D gaussian shape, and the neural network would learn to predict the amplitude, the horizontal/vertical position, and the horizontal/vertical width of the gaussian distribution.

Later, the assumption that the plasma profile had a particular shape was relaxed to allow for models with more degrees of freedom. Namely, a model was developed to take into account—among other parameters—ellipticity, triangularity and the Shafranov shift [21]. A neural network was used to learn the parameters of this model from measurements of a horizontal and a vertical neutron camera [22]. In this case, the neural network could learn up to 16 model parameters.

With the advent of deep learning [23], it became possible to train neural networks with many more layers and with a much larger number of parameters. One of the most successful applications of deep learning was image classification using convolutional neural networks (CNNs), where a 2D input image is transformed into a 1D output vector of class probabilities. By reversing the architecture of a CNN, it is possible to devise a deconvolutional neural network for plasma tomography, where a 1D input vector of bolometer measurements is transformed into a 2D output reconstruction of the plasma profile [24].

In this case, the neural network was trained to reproduce each single pixel of the tomographic reconstruction. For reconstructions with a resolution of \(200\times 120\) pixels, the network had 24,000 outputs, and was able to achieve a similarity score above 90% on previously unseen data. The main advantage is that, once trained, such network can compute hundreds or even thousands of reconstructions per second, making it possible to visualise the plasma profile over the course of an entire discharge [25] and, potentially, in real-time applications.

To conclude, the neural network tomography is the fastest among the three methods presented above, however, it is critically dependent on the quality and range of the training data. The Tikhonov regularisation is highly versatile, robust and independent on previous knowledge. Although it requires more computation time than the network tomography, it is considerably faster than any proper implementation of the Bayesian approach, including ML. The ML tomography represents the most sophisticated method which in various tests results in higher accuracy of tomographic reconstructions than the Tikhonov regularisation [26] and allows for semi-analytic calculation of the projection noise propagation. In this context it should be reminded that any tomographic method solves an ill-posed task, so that its accuracy can be significantly damaged by minor systematic errors including detector misalignments or varying sensitivities of individual detectors.

Progress in the Minimum Fisher Regularisation

The Minimum Fisher Regularisation (MFR) presents a widespread method of tomography in current tokamak research. It relies on Tikhonov regularisation with an iterative optimisation of the results so that the Fisher information of the reconstructed image is minimised:

$$\begin{aligned} I_F = \int \sum _{i,j=1}^2 \left| \frac{\partial g(x_1,x_2)}{\partial x_i}\frac{\partial g(x_1,x_2)}{\partial x_j}\frac{1}{g(x_1,x_2)}\right| \, {\mathrm {d}}x_1 \, {\mathrm {d}}x_2 \,\, . \end{aligned}$$
(10)

In plasma tomography on tokamaks, \(x_1\) and \(x_2\) correspond to the R and z coordinates (Figs. 1, 3). Optimisation according to Eq. 10 allows—in simple terms—for better spatial resolution (less smoothness) in regions with high levels of emissivity. The MFR method was first developed and applied at TCV [2] and at present it is implemented and updated, among others, at JET [9], COMPASS [4] and Tore Supra [27].

The first major extension of the MFR aimed at a possibility of rapid analyses of large amounts of data via temporal averaging of the smoothess matrix, see [3]. Nowadays the temporal averaging is hardly ever applied in practice due to the substantial increase in the CPU performance of the computers. With robust and rapid performance of MFR, further efforts have been invested into the development of real-time relevant version of MFR. As a result the idea of a rolling iteration was introduced, see [28]. In the rolling iteration, the time index of analysed data is increased by one with every new round of the iterative process. This allows for sufficient precision in the case of smooth, slowly evolving data from line integrated measurements. Importantly, in [28] it is demonstrated that the artefacts linked to sudden changes in data have also rather low and rapidly decreasing amplitude in MFR, i.e. the rolling iteration is stable. It can be concluded that a real-time version of the MFR is foreseeable.

Furthermore, preferential smoothness of the reconstructed image along the flux surfaces was introduced at JET due to the low number of the lines of view. As a result, the reconstructed image features slowly changing emissivity along the flux surfaces, while the emissivity gradient in plasma radius may be steep. This corresponds to the expected plasma emissivity behaviour. In oder to enforce this smoothness anisotropy, a new recipe for the regularisation matrix was implemented instead of Eq. 5:

$$\begin{aligned} H_{kj}=S(\eta )\sum _m^N B^T_{\Vert km}w_m B_{\Vert mj} + S(-\eta )\sum _m^N B^T_{\bot km}w_m B_{\bot mj} \, , \end{aligned}$$
(11)

In this set of equations, \({\mathbf {B_{\Vert }}}\) is the smoothing matrix in the direction parallel to the magnetic field, \({\mathbf {B_{\bot }}}\) the smoothing matrix in the direction perpendicular to the magnetic field, the weights \(w_m\) allow e.g. for implementation of the Minimum Fisher Information according to Eq. 10 (see [2] for details), and the function \(S(\eta )\) controls the anisotropy amplitude. In practice, the logistic sigmoid function \(S(\eta )=(1+e^{-\eta })^{-1}\) is applied. This amendment of the matrix \({\mathbf {H}}\) proved so reliable and beneficial [9] that it is nowadays applied as a routine feature in MFR.

Due to the non-linear character of the Minimum Fisher Information and to the sophisticated smoothing procedure it is not possible to determine analytically the error transmision from line integrated data to the emissivity distribution like in the ML method, see “Current Methods of Plasma Tomography” section. Instead, the Monte Carlo method was introduced, which tests statistically the MFR reconstruction response to random errors in data. The method was first used in extensive studies of the MFR performance at JET [29]. In this work it was shown that the MFR is stable against artefacts and that a Gaussian noise in data transmits with a good precision to a Gaussian noise in the reconstruction. Recently MFR contributed to detailed evaluation of accuracy and precision of the ITER neutron profile reconstruction from the simulated Radial Neutron Camera (RNC) data. The studies proved high robustness of the MFR method and acceptable level of precision of the neutron profiles with a sufficient temporal resolution, in particular in the high performance discharges [8].

Besides plasma tomography, MFR was successfully applied in unfolding of neutron spectra \(\varPhi _j\) from the pulse height data \(\varPsi _i\) measured by neutron scintillation detector with detailed knowledge of its response function with matrix elements \(R_{ij}\), for details see [30] :

$$\begin{aligned} \varPsi _i = \sum _{j} R_{ij}\varPhi _j , \end{aligned}$$
(12)

Unlike in plasma tomography, in the case of unfolding the task need not be underdetermined so that the L-curve principle based on data statistics can be employed in search of the value of the regularisation factor \(\lambda\) instead of the \(\chi ^2\) Pearson’s method, see [30, 31]. Notice that the unfolding is still an ill-posed problem with a tendency to create artefacts, therefore a reliable calibration of the response matrix \({\mathbf {R}}\) was indispensable.

Based on the observed robustness of the MFR method as well as on experience with spectral unfolding, a new approach to analyses of spatially resolved pulse height data from the RNC scintillation detectors at ITER can be proposed. In contrast to the step-by-step approach presented in [32], where tomography analyses is run for separate energy bins, it is recommended to directly combine the contribution (i.e. geometric) matrix and the response matrix of the energy calibrated detectors into a single inversion problem as follows:

  • Denote \(f_{ik}\) elements of a matrix of data where the rows correspond to \(i=1,\ldots ,P\) line integrated measurements from the P scintillation spectrometers and the columns to \(k=1,\ldots ,C\) pulse height measurements from each spectrometer.

  • Seek elements of a matrix of spectrally resolved emissivities \(g_{jl}\), where the rows \(j=1,\ldots ,N\) correspond to the pixel index in the spatially resolved mesh of pixels, and the columns \(l=1,\ldots ,B\) to discrete bins of the unfolded neutron energies.

  • The double inversion problem is described—according to Eqs. 2 and 12—by the following set of equations

    $$\begin{aligned} f_{ik} = \sum _l^B R_{kl}\sum _j^N T_{ij}g_{jl} + \xi _{ik} , \end{aligned}$$
    (13)

  • This set of equations can be re-indexed so that a standard regularisation procedure (e.g. the MFR code) can be applied

    $$\begin{aligned} f_{\alpha } = \sum _\beta ^{B\times N} S_{\alpha \beta }g_{\beta } + \xi _{\beta } \;\;\;\; \rightarrow \;\;\;\; g_{\beta } = \sum _\alpha ^{C\times P}G_{\beta \alpha }f_{\alpha } \end{aligned}$$
    (14)
Fig. 4
figure 4

Preliminary results of the combined tomography and unfolding MFR code, with phantom functions in left column and combined 2D tomography and unfolding results in the right column. In the top row, the spatial distribution of the DT neutron emissivity is shown, while the bottom row presents the model temperature profile and the unfolded temperature distribution, based on spectral width of the DT neutron energy. The simulation was run for the JET neutron emissivity profile monitor geometry [12]

Full size image

The main advantage of the proposed procedure is that the combined tomography and spectral unfolding problem is treated as a single ill-posed task, so that the amplification of artefacts is prevented.

The method was tested on phantom functions with random noise. The contribution matrix \({\mathbf {T}}\) was based on geometry of the JET neutron emissivity profile monitor [12], and the response matrix \({\mathbf {R}}\) had values of the well calibrated liquid scintillation detector that was successfully tested at JET as a neutron spectrometer [33]. One of the test functions and the resulting reconstruction is exemplified in Fig. 4. In this example, a simple gaussian profile of neutron emissivity was used as a phantom function on the grid of \(10\times 15\) pixels, with 100 energy bins in each. The phantom neutron spectrum was based on the DT neutron emission (14.1 MeV) with Doppler broadening corresponding to a parabolic plasma temperature profile. As a result, the combined procedure resulted in considerably improved precision of reconstruction (by approx. 15%) compared to the step-by-step unfolding and tomography. It can be concluded that the preliminary results are promising.

Since its introduction more than 20 years ago, the MFR also went through several code remakes and numerical optimisations. The most important extension happened in 2012 when the original MatLab MFR package was ported to the Python platform with new numerical options included [9]. The Python version was significantly refactored again in 2017, with a new streamlined hierarchic module structure and Mercurial version control. In the version currently maintained for the JET tomography, a new implementation of the anisotropic smoothness matrix based on gradient of the poloidal magnetic flux (instead of the flux surface interpolation) was introduced, see Fig. 5. Besides, the same mesh of pixels was pre-defined for the three JET tomography diagnostic systems: the neutron, the soft X-ray and the bolometry cameras, and a new simple access to JET data is under development. A new graphical user interface (GUI) was introduced, which allows, among others, to evaluate evolution of radiation in between regions of interest pre-defined by the user, profiting from data of all the three diagnostic systems.

Fig. 5
figure 5

Left: Two-dimensional plot of the normalised magnetic flux function \(\varPsi _N(R,z)\) as reconstructed by EFIT at JET, including Last Closed Flux Surface in black (the value of the flux is normalised to the value at this contour). The red arrows mark the gradient of the \(\varPsi\) function and the white arrows the perpendicular direction. Right: A similar plot of inverse tangent of the direction perpendicular to the gradient of \(\varPsi\). This function is then used in order to find proper weights for the derivative (smoothing) matrix with preferential smoothing along the magnetic flux contours

Full size image

The MFR method was also applied in reconstruction of data from fast 2D matrix cameras with tangential view of the plasma in the visible region (based on the plasma axial symmetry) on the COMPASS tokamak, see e.g. [34], and tentatively to unfolding in analyses of the data from the activation probes at JET, see [35]. Its potential to reconstruct the plasma current density from the magnetic data is under discussion. Last, but not least a possible merger of the MFR with the neural networks can be considered in future, as proposed in [36]. According to this scheme, MFR regularisation parameters in Eqs. 4 and 11 would be determined by a trained neural network, combining robustness of MFR with the real-time relevance of the neural networks.

Conclusions and Outlook

In this contribution, research and development efforts in plasma tomography were presented. Particular focus was given to the large number of applications that a reliable method resolving the inversion problem may offer in fusion data analyses. Indeed, the potential of the tomography development for fusion is still to be exploited. At present, several methods with their advantages and disadvantages compete, in which a rich set of constraints may (but need not) be applied. Importantly, the JET contributors involved in tomography analyses have been currently working on quantitative comparison of performance of the three methods presented in the second part of this contribution. The conclusions are yet to be drawn, however, the preliminary results demonstrate that (1) with sparse data, the room for improvement is limited so that in this respect, augmented data precision is rather called for, and (2) the correct attitude is to maintain a few different inversion methods in order to be able to critically compare their results in analyses of important data events.

For future fusion reactors, it will be instrumental to develop a real-time tomography algorithm with low susceptibility of developing major artefacts. This task is of a particular importance in determination of the fusion neutron emissivity distribution. In this respect, several works and conceptual studies based on Minimum Fisher Regularisation method were also presented in this contribution. Obviously, similar ideas can be also anticipated in other plasma tomography methods.