The interaction between particles in fusion plasmas leads to the emission of electromagnetic radiation ranging from the infrared to the SXR band. Most of the radiation is not absorbed by the plasma rather lost outside. By measuring the radiation in different spectral ranges, information about various plasma parameters can be obtained.

The measurement of SXR intensity is an important diagnostic method to get information about the magnetic field topology in tokamak devices.^[1] By the tomographic inversion of line-integrated emissivity measurements, the distribution of the SXR intensity in the poloidal cross section can be obtained. Equi-intensity surfaces can approximately reflect the configuration of the magnetic flux surface, which is essential for the study of the core MHD mode structure and impurity transport. The measurement of the fusion plasma radiation is also one of the primary methods to estimate the plasma energy loss, which affects the power balance and the plasma confinement. The AXUV photodiode detector system can measure the radiation power and radiation emissivity distribution in plasma discharges.^[2,3] With the radiation power and emissivity measured by the AXUV system, the 2D radiation distribution in the poloidal cross section can be reconstructed by tomographic inversion.

In tokamaks, the line-integrated signals are measured from relatively small apertures in the diagnostics systems. By combining a limited number of line-integrated signals,^[4] tomographic inversion can reconstruct the 2D distributions of physical quantities. Various tomographic techniques have been applied to fusion diagnostics, such as the Cormack method,^[5] maximum entropy method,^[6] and minimum Fisher information method.^[7] In recent studies, it is prevalent to use either the maximum entropy method or the Fisher minimal information method. However, they have no obvious advantages in uncertainty analysis.^[8] The random errors in experiments and the systematic errors in the modelling procedure are unavoidable, and these errors may even cause great uncertainties in the reconstruction results. Hence, a method with uncertainty estimation capability is demanded.

Typical tomography problems involve the existence, uniqueness and process stability of the solution. In recent years, a method that introduces a Bayesian probability theory-based GPT technique has been developed to reconstruct the SXR emissivity distribution. It has certain advantages in uncertainty analysis^[8] and overcomes typical tomographic problems by considering a probabilistic approach. This method has already been implemented in the SXR systems on multiple fusion devices, such as HL-2A^[7] and WEST.^[9] Here we employ it in the SXR and AXUV systems on EAST. This is the first application of this method in the AXUV diagnostic in fusion devices. Generally, it is difficult to reconstruct the 2D radiation distribution by only one direction detector arrays, but the GPT method is superior in this respect. From the GPT reconstructions, redistribution caused by MHD events and impurity transport can be clearly identified. By comparing with the previous reconstructions and combining with the analysis of other diagnostic signals, the reconstructions captured by GPT are validated. Additionally, the evolution of the core expulsion and the reformation details can be obtained. This is meaningful for the fine structure analysis of various MHD events in future works. A method of calling FORTRAN from MATLAB through a FORTRAN MEX file has been used to accelerate the tomography code in this work. Therefore, the GPT method has high potential for real-time application due to its fast computation.

The GPT based on Bayesian data analysis for the SXR and AXUV diagnostic systems on EAST is presented in this paper, which is organized as follows: the method is described in Section 2. The implementation of the GPT method in the EAST SXR and AXUV systems are described in Section 3 and Section 4, respectively. Section 5 summarizes this work.

In nuclear fusion devices, tomographic systems typically have relatively small apertures. Hence, the measurement of the j-th line of sight (LOS) d_j of a tomographic system can be approximated by the integration of the emissivity filtered by the detector along the corresponding LOS:

$$ \begin{eqnarray}{d}_{j}=\displaystyle {\int }_{{\rm{LOS}}}E{\rm{d}}{r}_{j}.\end{eqnarray}$$ (1)

The reconstruction region is the area within the last closed flux surface as shown in Fig. 1, which is determined based on both the actual extension of emitting zone and the coverage of available LOS. A square grid of 50 × 50 pixels discretizes the target region in space to solve the tomography problem. The emissivity within each pixel is treated as a constant. The light going into each specific LOS passes through many pixels, and the sum radiation of these pixels give the line-integrated radiation intensity. The advantage of this approach is that there is no limit to the shape of the magnetic flux surfaces, and the line-integrated emissivity can be converted into a linear system of equations by

$$ \begin{eqnarray}{{\boldsymbol{d}}}_{m}={{\boldsymbol{R}}}_{m\times n}\cdot {{\boldsymbol{E}}}_{n}+{\boldsymbol{\varepsilon }},\end{eqnarray}$$ (2)

Bayesian data analysis includes the measurements and prior information in the form of probability into the posterior probability distribution equation, from which unknown quantities can be inferred. This method contains all the model information provided by measurements and prior information of the unknown quantities before any measurement.^[11]

By applying the Bayesian approach to the tomography problem, the quantities both E_n and d_m can be expressed by probability distributions. The unknown random variable θ is introduced to summarize the prior information of the model. Then the relationship between E_n, d_m, and θ in a Bayesian framework is

$$ \begin{eqnarray}\begin{array}{lll}p({{\boldsymbol{E}}}_{n}|{{\boldsymbol{d}}}_{m},{\boldsymbol{\theta }}) & = & \displaystyle \frac{p({{\boldsymbol{d}}}_{m}|{{\boldsymbol{E}}}_{n},{\boldsymbol{\theta }})\cdot p({{\boldsymbol{E}}}_{n}|{\boldsymbol{\theta }})}{p({{\boldsymbol{d}}}_{m}|{\boldsymbol{\theta }})}\\ & \sim & p({{\boldsymbol{d}}}_{m}|{{\boldsymbol{E}}}_{n},{\boldsymbol{\theta }})\cdot p({{\boldsymbol{E}}}_{n}|{\boldsymbol{\theta }}),\end{array}\end{eqnarray}$$ (3)

The GPT method is based on Gaussian process modeling of the emissivity field, providing a posterior distribution in the context of bayesian inference. It is related to Gaussian process regression, which is a widely used nonparametric regression technique that can reduce the computational complexity compared with other algorithms.^[12] The distribution of a Gaussian process is the joint distribution of infinitely many normally distributed random variables and it is a distribution over a continuous domain of time or space. That is to say a Gaussian process (GP) is a generalization of the multivariate normal (MVN) distribution of the function space. It is described by a mean function μ and a covariance function Σ, where GP ∼ N(μ, Σ). The mean vector μ describes the expected value of the distribution, and each component of its distribution describes the mean of the corresponding pixel. Σ models the variance of each pixel. The prior emissivity distribution in n pixels with r_i coordinates is assumed to be multivariate Gaussian, its covariance matrix is defined by:

$$ \begin{eqnarray}{{\boldsymbol{\varSigma }}}_{{\rm{E}}}=\left[\begin{array}{ccc}k({r}_{1},{r}_{1}) & \cdots & k({r}_{1},{r}_{n})\\ \vdots & \ddots & \vdots \\ k({r}_{n},{r}_{1}) & \cdots & k({r}_{n},{r}_{n})\end{array}\right],\end{eqnarray}$$ (4)

$$ \begin{eqnarray}{k}_{{\rm{SE}}}={\sigma }_{{\rm{f}}}^{2}\exp \left(-\left(\displaystyle \frac{{{\boldsymbol{d}}}_{\perp }^{2}}{2{\sigma }_{{\rm{l}}\perp }^{2}}+\displaystyle \frac{{{\boldsymbol{d}}}_{\parallel }^{2}}{2{\sigma }_{{\rm{l}}\parallel }^{2}}\right)\right),\end{eqnarray}$$ (5)

$$ \begin{eqnarray}f(x)=\displaystyle \frac{1}{\sigma \sqrt{2\pi }}\exp \left[-\displaystyle \frac{{(x-\mu )}^{2}}{2{\sigma }^{2}}\right].\end{eqnarray}$$ (6)

$$ \begin{eqnarray}\begin{array}{lll} & & p({{\boldsymbol{E}}}_{n}|{\boldsymbol{\theta }})=\displaystyle \frac{1}{{(2\pi )}^{n/2}{\left|{{\boldsymbol{\varSigma }}}_{{\rm{E}}}^{{\rm{prior}}}\right|}^{1/2}}\exp \left[-\displaystyle \frac{1}{2}{({{\boldsymbol{E}}}_{n}-{{\boldsymbol{\mu }}}_{{\rm{E}}}^{{\rm{prior}}})}^{{\rm{T}}}{{\boldsymbol{\varSigma }}}_{{\rm{E}}}^{{\rm{prior}}-1}({{\boldsymbol{E}}}_{n}-{{\boldsymbol{\mu }}}_{{\rm{E}}}^{{\rm{prior}}})\right],\\ & & \end{array}\end{eqnarray}$$ (7)

$$ \begin{eqnarray}\begin{array}{lll} & & p({{\boldsymbol{d}}}_{m}|{{\boldsymbol{E}}}_{n},{\boldsymbol{\theta }})=\displaystyle \frac{1}{{(2\pi )}^{m/2}{|{{\boldsymbol{\varSigma }}}_{d}|}^{1/2}}\exp \left[-\displaystyle \frac{1}{2}({{\boldsymbol{R}}\cdot {{\boldsymbol{E}}}_{n}-{{\boldsymbol{d}}}_{m})}^{{\rm{T}}}{{\boldsymbol{\varSigma }}}_{d}^{-1}({\boldsymbol{R}}\cdot {{\boldsymbol{E}}}_{n}-{{\boldsymbol{d}}}_{m})\right],\\ & & \end{array}\end{eqnarray}$$ (8)

$$ \begin{eqnarray}\begin{array}{lll} & & {{\boldsymbol{\varSigma }}}_{d}=\left[\begin{array}{ccc}{(0.05\cdot {d}_{1})}^{2} & \cdots & 0\\ \vdots & \ddots & \vdots \\ 0 & \cdots & {(0.05\cdot {d}_{m})}^{2}\end{array}\right].\end{array}\end{eqnarray}$$ (9)

$$ \begin{eqnarray}p({{\boldsymbol{E}}}_{n}|{{\boldsymbol{d}}}_{m},{\boldsymbol{\theta }})\sim \exp \left[\begin{array}{c}-\displaystyle \frac{1}{2}{({\boldsymbol{R}}\cdot {{\boldsymbol{E}}}_{n}-{{\boldsymbol{d}}}_{m})}^{{\rm{T}}}{{\boldsymbol{\varSigma }}}_{d}^{-1}({\boldsymbol{R}}\cdot {{\boldsymbol{E}}}_{n}-{{\boldsymbol{d}}}_{m})\\ -\displaystyle \frac{1}{2}{({{\boldsymbol{E}}}_{n}-{{\boldsymbol{\mu }}}_{{\rm{E}}}^{{\rm{prior}}})}^{{\rm{T}}}{{\boldsymbol{\varSigma }}}_{{\rm{E}}}^{{\rm{prior}}-1}({{\boldsymbol{E}}}_{n}-{{\boldsymbol{\mu }}}_{{\rm{E}}}^{{\rm{prior}}})\end{array}\right],\end{eqnarray}$$ (10)

$$ \begin{eqnarray}\begin{array}{lll} & & p({{\boldsymbol{E}}}_{n}|{{\boldsymbol{d}}}_{m},{\boldsymbol{\theta }})\sim \exp \left[-\displaystyle \frac{1}{2}{({{\boldsymbol{E}}}_{n}-{{\boldsymbol{\mu }}}_{{\rm{E}}}^{{\rm{post}}})}^{{\rm{T}}}{{\boldsymbol{\varSigma }}}_{{\rm{E}}}^{{\rm{post}}-1}({{\boldsymbol{E}}}_{n}-{{\boldsymbol{\mu }}}_{{\rm{E}}}^{{\rm{post}}})\right].\end{array}\end{eqnarray}$$ (11)

In tokamaks, SXR measurements are performed to analyze the MHD instabilities. Using the line-integrated SXR measurements from multiple cameras to reconstruct the emissivity distribution can solve the problems such as determining the mode structure and rotational perturbation in the plasma core, which is of great significance to the study of the magnetic flux surface topology of the core plasma and the related MHD instabilities.

The SXR system on EAST before 2019 consists of 3 arrays, 2 up-down symmetric horizontal arrays and a vertical array. The vertical array has recently been developed into 2 arrays and the detailed layout of the detector arrays is shown in Fig. 1. With a high spatial resolution of approximately 2.5 cm and a temporal resolution of up to 10 μs, it is able to observe MHD instabilities on the poloidal cross section. The radiation reconstruction area consists of 50 × 50 pixels. Measurements from horizontal arrays and from both horizontal plus vertical arrays have been reconstructed. No clear difference between the reconstructions of these 2 cases have been found. Therefore, only the case reconstructed from the solely horizontal arrays is shown in this work. All reconstructed emissivity mentioned in this work have been normalized between 0 and 10 for the benefit of numerical stability and comparison.

SXR tomography is an effective tool for the studying of sawtooth crashes. The 2D images of the ramp-up phase of sawteeth reconstructed from SXR measurements during EAST discharge #70754 are shown in Fig. 2. As illustrated in the tomographic patterns, the nearly circular SXR emissivity core resembles the m/n = 1/1 oscillation and has ideal MHD nature before the sawtooth crash.

Figure 3 shows the 2D reconstructions of the internal kink mode patterns. Since the SXR original signal transmittance has an approximately linear relationship with the temperature, the SXR emissivity distribution can roughly reflect the heat distribution. The reconstructions show that during the SXR intensity ramp-up phase, the heat slowly accumulates in the core, as shown by pattern B to C and E to F, while the heat releases outwards during the SXR intensity decline phase. It is obvious that during the sawtooth crash phase, a considerable amount of heat rapidly flows from the core to the edge. According to the study of this target plasma, it is proved that the GPT method can provide a clear hot core shape by reconstructing the magnetic surface structure of the m = 1 mode.

The MHD mode structure has been analyzed using singular value decomposition (SVD). SVD considers both the time and spatial structure dependencies of different signals. SVD is a powerful mathematical tool, which is widely used in the study of the spatiotemporal behavior of magnetic islands, SXR emissivity and other fluctuations in fusion plasmas.^[15] The matrix-form of the SVD can be neatly expressed as An×p=Un×n⋅Sn×p⋅Vp×pT, where U_{n × n} and V_{p × p} are orthogonal. If each row of A_{n × p} contains a time series of a physical quantity at different locations, the columns of U_n × n represent spatial eigenvectors or topos, whereas the columns of V_{p × p} can be considered as temporal eigenvectors or chronos. S_{n × p} is diagonal and contains singular values. Due to the descending order of the singular values, the first topos/chronos pairs contain most of the information of A_{n × p}, whereas the components with the smallest singular values merely contain noise. In this work, the decomposed matrix consists of reconstructed emissivity vectors at p time slices. The first three topos corresponding to the m/n = 1/1 mode structure of the kink mode are shown in Fig. 4.

Figure 5 shows the reconstruction results for shot #70754 at 2.8284 s. The emissivity gradient shown in Fig. 5(b) has a hollow shape, which means the heat will transport from outside-in. This is consistent with the SXR intensity ramp-up phase. The diagonal elements of posterior covariance plotted in Fig. 5(c) is relatively low, indicating low uncertainty and relatively accurate reconstructed results. Moreover, a consistency check of that solution has been performed by checking whether the misfits between reconstructed and measured data are reasonable within the assumed measurement uncertainty. The so-called reconstructed signals are calculated from the reconstructed emission and the path lengths matrix according to Eq. (2). The reconstructed signals match well with the measured signals as evident from Fig. 5(d). These indicate that the GPT method can provide consistent and reliable reconstructions of the SXR measurements for EAST.

Sawtooth-like internal disruption is thought to appear due to interactions between two modes with identical toroidal mode numbers or double tearing modes.^[16,17] The toroidal MHD mode coupling due to the toroidal geometry of the tokamak and non-circular plasma cross-section is thought to be the dominant reason for sawtooth-like internal disruptions. The rapid growth of the m = 2 island or the coexistence of tearing modes with different helicities can lead to sawtooth-like internal disruptions. Different to normal sawteeth, sawtooth-like minor disruptions can cause loss of plasma confinement. This is shown in Fig. 6, where the electron density drops significantly when the sawtooth-like disruption occurs.

In this section, we present the reconstructed results of a sawtooth-like internal disruption. Figure 7 shows the tomographic reconstructions from SXR measurements during a sawtooth-like disruption event in EAST discharge #61450. It is anticipated that this mode erupts over a larger area than normal sawteeth. In addition, the conversion from an m = 2 islands to an m = 1 mode was obtained during the sawtooth-like internal disruption crash phase shown in Fig. 8, as suggested by previous studies.^[18] Figure 9 shows that the reconstructed signals match well with the measured signals.

The radiation power measurement system installed on EAST uses absolute extreme ultraviolet (AXUV) photodiodes. In 2015, four pinhole cameras with a total of 64 channels have been developed in the horizontal port to view the poloidal cross-section. As shown in Fig. 10(a), each camera contains a 16-channel array with a spatial resolution of 3 cm. Two additional vertical pinhole cameras with a total of 35 channels have been developed to view the divertor through the bulk plasma in 2019. This setup is shown in Fig. 10(b), in which vertical channels 21 to 25 are absent due to overlaps. The radiation power distribution measured by the AXUV system during the 2016 campaign are reconstructed using the GPT method. Only measurements of the horizontal cameras with 61 (3 bad LOS were eliminated) LOS have been employed in this work. The radiation reconstruction area consists of 50 × 50 pixels.

Figure 11 shows the GPT reconstructions from AXUV measurements obtained during EAST ELMy H mode discharge #63720. From the time slice A to E, the total radiation power in the plasma core gradually increases, which is due to the improved plasma confinement and the accumulation of molybdenum impurities in the plasma core. This is consistent with the evolution of molybdenum impurities in the plasma core as evident from Fig. 11(b). The concentration of molybdenum is measured by the Extreme UltraViolet spectrometer (EUV) diagnostics.^[19] From the radiation power contour plot of discharge # 63720 in Fig. 12, the radiation from the plasma core trails off after the L-mode onset. This is also related to the gradual decrease of core molybdenum concentration. These evidences support that the reconstructed results of AXUV measurements are reasonable and physically credible.

In general, it is difficult to reconstruct the radiation power of the boundary region by using only the measurements of AXUV horizontal arrays. The GPT method has certain advantages in this aspect. Experimental traces of shot #69042 in EAST demonstrate that the core radiation power could be mitigated by gas-puffing. As figure 13(e) shows, the core magnetic disturbance and the concentration of tungsten began to increase at t = 3 s, and the radiation power of the plasma core was also enhanced as shown in Fig. 13(d). At 4.5 s, the plasma confinement was degraded. With the injection of neon from the mid-plane at t = 5 s, the plasma density and energy storage increased, and the plasma confinement was improved. The radiation intensity of the tungsten impurity in the plasma core began to drop sharply, while the PXUV signals and the core Mirnov signal both showed that the amplitude of the plasma magnetic disturbance were significantly mitigated due to the injection of neon. The GPT reconstruction results of AXUV data of discharge #69042 are given by Fig. 15, which can clearly reflect the radiation power distribution before and after the injection of neon. This is in good agreement with the above analysis and the AXUV contour plot shown in Fig. 14(a).

The injection of light impurities such as nitrogen, neon, or argon can significantly reduce the thermal deposition at the first wall, i.e., reduce damage to the first wall and thus reduce impurities in the plasma core. Figure 16 gives the profiles of electron temperature and density measured by the Thomson scattering system,^[20] which shows that the electron temperature at the area outside ρ_t = 0.5 dropped about 500 eV and the density increased by about 30 % after the injection of neon. Here ρ_t is the normalized poloidal flux. The lower edge electron temperature means less thermal deposition at the first wall, which support the mechanism mentioned above. The SVD analysis in Fig. 17 shows that the ratio of the weight factor of m = 1 and m = 2 mode increased from 37 to 65 after the gas-puffing, which means that the further development of the m = 2 mode was suppressed or the m = 1 mode was promoted. The suppression of the m = 2 mode means fewer impurities are generated and transmitted to the plasma core, while the promotion of the m = 1 mode can help to exhaust the impurities from plasma core, both of which can reduce the impurities in the plasma core. There was more than one kind of mode during this period, and that led to the radiating shell-like structure in the reconstructions. From the above, injecting a certain gas into the mid-plane can modulate the MHD behavior, reduce the accumulation of impurities in the plasma core, and thus improve the confinement. This provides a useful reference for the radiation feedback control and plasma instability control experiments.

A new tomographic approach based on the GPT method has been developed, and reliable 2D reconstructions of the emissivity distribution are obtained. Examples of GPT implementation in the SXR and AXUV diagnostics on EAST show that Bayesian Gaussian process tomography is a feasible method that can realize robust reconstructions in fusion plasma applications. It is a remarkable fact that this method can provide credible 2D emissivity field by reconstructing from only horizontal arrays. It is the first time that this method is applied to the AXUV diagnostic, and satisfactory results have been achieved. Considering its fast computational speed, it is worth developing this method to provide real-time feedback on impurity transport and fast MHD events. Additionally, it is compelling for fine structure analysis of different MHD events in future works.

Acknowledgment. The authors would like to acknowledge the help provided by EAST Team.