Electron Density Reconstruction of Solar Coronal Mass Ejections Based on a Genetic Algorithm: Method and Application

TL;DR

This paper introduces a genetic algorithm-based method for three-dimensional electron density reconstruction of solar coronal mass ejections, validated with synthetic and real observational data, improving understanding of CME structures.

Contribution

The paper presents the first application of a genetic algorithm for CME electron density reconstruction, demonstrating its effectiveness with synthetic and actual observational data.

Findings

GRM accurately reconstructs synthetic CMEs with various orientations.

The method performs well even with brightness incompleteness and noise.

Application to real data shows promising results for CME analysis.

Abstract

We present a new method to reconstruct the three-dimensional electron density of coronal mass ejections (CMEs) based on a genetic algorithm, namely the genetic reconstruction method (GRM). GRM is first applied to model CMEs with different orientations and shapes. A set of analytic model CMEs from Gibson and Low is employed to produce synthetic CME images for GRM reconstruction. Model CMEs with longitudes of 0, 45, 90, 135, 180 degrees and latitudes of 0, 15, 30, 45 degrees are used to test the performance of GRM. The model CMEs are obscured with a simulated occulter of a coronagraph to determine the influence of CME brightness incompleteness. We add random noise to some synthetic CME images to test the performance of GRM. The CME reconstructions are carried out using synthetic data from Solar Terrestrial Relations Observatory (STEREO) A and B with a separation angle of 90 degrees and from STEREO A and the Solar and Heliospheric Observatory (SOHO) with a separation angle of 73 degrees. The Pearson correlation coefficient and the mean relative absolute deviation are calculated to analyze the similarities in brightness and electron density between the model and reconstructed CMEs. Comparisons based on the similarity analysis under various conditions stated above give us valuable insights into the advantages and limitations of GRM reconstruction. The method is then applied to real coronagraph data from STEREO A and B, and SOHO on 2013 September 30.

Tables

Table 1. Table 1: M ​ R ​ A ​ D 𝑀 𝑅 𝐴 𝐷 MRAD and P ​ C ​ C 𝑃 𝐶 𝐶 PCC of electron density between model and reconstructed CMEs. Please click on the blue values of central longitude and latitude at the first column to go to corresponding FOS.

Modeled CME position	$MRAD$ of xyz space	$PCC$ of xyz space	$PCC$ of xy plane	$PCC$ of xz plane	$PCC$ of yz plane
(longitude, latitude)	1st , 2nd , 3rd	1st , 2nd , 3rd	1st , 2nd , 3rd	1st , 2nd , 3rd	1st , 2nd , 3rd
0, 0	0.612, 0.620, 0.592	0.784, 0.782, 0.780	0.848, 0.883, 0.899	0.634, 0.642, 0.666	0.957, 0.960, 0.949
0, 15	0.615, 0.608, 0.622	0.799, 0.798, 0.816	0.905, 0.891, 0.886	0.750, 0.751, 0.779	0.933, 0.949, 0.935
0, 30	0.652, 0.651, 0.666	0.808, 0.795, 0.800	0.892, 0.883, 0.864	0.853, 0.851, 0.851	0.965, 0.964, 0.959
0, 45	0.744, 0.741, 0.732	0.796, 0.800, 0.793	0.841, 0.849, 0.817	0.864, 0.871, 0.863	0.962, 0.964, 0.965
45, 0	0.569, 0.577, 0.569	0.845, 0.870, 0.862	0.923, 0.915, 0.925	0.762, 0.805, 0.796	0.992, 0.991, 0.991
45, 0aaThe CME is obscured by the modelled occulter of coronagraph.	0.640, 0.635, 0.632	0.737, 0.772, 0.759	0.703, 0.786, 0.739	0.502, 0.586, 0.562	0.988, 0.988, 0.986
45, 30	0.631, 0.630, 0.620	0.840, 0.823, 0.809	0.857, 0.850, 0.831	0.846, 0.803, 0.798	0.991, 0.990, 0.989
90, 0	0.608, 0.599, 0.612	0.810, 0.789, 0.799	0.938, 0.944, 0.924	0.665, 0.580, 0.652	0.946, 0.937, 0.927
90, 30	0.651, 0.668, 0.656	0.782, 0.798, 0.787	0.893, 0.890, 0.908	0.843, 0.802, 0.860	0.931, 0.926, 0.946
135, 0	0.592, 0.597, 0.581	0.838, 0.868, 0.855	0.980, 0.975, 0.977	0.685, 0.679, 0.726	0.923, 0.921, 0.900
135, 0aaThe CME is obscured by the modelled occulter of coronagraph.	0.578, 0.589, 0.579	0.821, 0.811, 0.817	0.790, 0.779, 0.782	0.649, 0.654, 0.653	0.927, 0.884, 0.932
135, 30	0.621, 0.618, 0.611	0.822, 0.821, 0.847	0.957, 0.953, 0.955	0.913, 0.899, 0.910	0.908, 0.898, 0.923
135, 30aaThe CME is obscured by the modelled occulter of coronagraph.	0.630, 0.643, 0.620	0.813, 0.835, 0.840	0.945, 0.953, 0.955	0.888, 0.872, 0.899	0.874, 0.886, 0.896
135, 30bbRandomly-distributed noise are added to the synthetic CME observations.	0.740, 0.727, 0.711	0.795, 0.801, 0.814	0.820, 0.881, 0.862	0.843, 0.838, 0.873	0.851, 0.826, 0.854
135, 30ccRandomly-distributed noise $\times$10 are added to the synthetic CME observations.	0.894, 0.897, 0.898	0.206, 0.398, 0.408	0.507, 0.284, 0.414	0.252, 0.499, 0.433	0.162, 0.381, 0.438
135, 30ddRandomly-distributed noise $\times$100 are added to the synthetic CME observations.	0.999, 0.999, 0.999	0.139, 0.184, 0.164	0.231, 0.139, 0.132	0.157, 0.218, 0.230	0.204, 0.112, 0.170
135, 30eeThe CME is reconstructed from modelled coronagraph images of STEREO A and SOHO as observation 1 and 2.	0.719, 0.710, 0.705	0.688, 0.731, 0.709	0.851, 0.857, 0.868	0.762, 0.773, 0.740	0.759, 0.763 0.744
135, 30ffThe CME is reconstructed from modelled coronagraph images of STEREO A, SOHO and B as observation 1, 2 and 3.	0.750, 0.761, 0.755	0.728, 0.716, 0.704	0.871, 0.869, 0.868	0.786, 0.822, 0.833	0.773, 0.773, 0.781
180, 0	0.596, 0.610, 0.608	0.789, 0.736, 0.788	0.871, 0.790, 0.891	0.653, 0.562, 0.649	0.958, 0.949, 0.939
180, 30	0.628, 0.655, 0.643	0.773, 0.792, 0.821	0.880, 0.856, 0.884	0.759, 0.775, 0.812	0.965, 0.953, 0.956
ABSggThe CME is reconstructed from real coronagraph images of STEREO A, B and SOHO as observation 1, 2 and 3.	0.000, 0.235, 0.252	1.000, 0.874, 0.836	1.000, 0.913, 0.927	1.000, 0.835, 0.791	1.000, 0.892, 0.780
ASBhhThe CME is reconstructed from real coronagraph images of STEREO A, SOHO and B as observation 1, 2 and 3.	0.000, 0.323, 0.299	1.000, 0.801, 0.837	1.000, 0.854, 0.936	1.000, 0.797, 0.847	1.000, 0.803, 0.898

Equations

$$B_{obs}={\int }_{{\theta }_1}^{{\theta }_2}{B}_e(\theta )N_e(\theta )d\theta +\eta .$$

$$cf=|\sum _{{\theta }_1}^{{\theta }_2}{B}_e(\theta )N_e(\theta )-B_{obs}{|}^p+regularization,$$

$$cf=|\sum _{{\theta }_1}^{{\theta }_2}{B}_e(\theta )N_e(\theta )-B_{obs}|.$$

$$N_e(i,k)=\frac{B_{obs}(i,k)}{B_e(\theta (i,k))}.$$

$$B_e(\theta (i,k))=\frac{\pi \sigma }{2z^2}[2((1-u)C+uD)-((1-u)A+uB){\cos }^2\theta (i,k)].$$

$$N_v(i,j,k)=\frac{L(i,j,k)}{{\sum }_{j=1}^{N_j}L(i,j,k)}{N}_e(i,k).$$

$$B_{av}(i,j,k)=B_e({\theta }_a(i,j,k))N_v(i,j,k),$$

$$B_e({\theta }_a(i,j,k))=\frac{\pi \sigma }{2z^2}[2((1-u)c+uD)-((1-u)A-uB){\cos }^2{\theta }_a(i,j,k)].$$

$$B_a(i,k)=\sum _{j=1}^{N_j}{B}_{av}(i,j,k).$$

$$B_b(j,k)=\sum _{i=1}^{N_i}{B}_{bv}(i,j,k)$$

$$cf(i,k)=|B_a(i,k)-B_{obs}(i,k)|$$

$$cf(j,k)=|B_b(j,k)-B_{obs}(j,k)|.$$

$$Fi{t}_q=\frac{B_{obs}}{c{f}_q}.$$

$$Fit[sort(Fit)]=findgen(N_q)+1,$$

$$Fi{t}_{qg}(i,k)=exp(Fi{t}_q(i,k)\frac{g}{N_g}0.03)$$

$$Fi{t}_{qg}(j,k)=exp(Fi{t}_q(j,k)\frac{g}{N_g}0.03),$$

$$L_m(i,j,k)\rightleftarrows L_n(i,j,k),$$

$$L_m(i,j,k)=L_m(i,j,k)+L_m(i,j,k)\times D_r,$$

$$\frac{1}{4\pi }(\nabla \times \mathbf{B})\times \mathbf{B}-\nabla p-\rho \mathbf{g}=0$$

$$\frac{1}{4\pi }(\nabla \times \mathbf{B})\times \mathbf{B}-\nabla p=0$$

$$\frac{1}{4\pi }(\nabla \times \mathbf{B})\times \mathbf{B}-\nabla p-\rho \mathbf{g}=\rho \left[\frac{\partial \mathbf{v}}{\partial t}+(\mathbf{v}\cdot \nabla )\mathbf{v}\right]$$

$$\rho =f(B_r,B_{\theta },B_{\varphi }).$$

$$MRAD=\frac{1}{N}\sum _{i=1}^N\frac{|x_i-y_i|}{x_i}$$

$$PCC=\frac{{\sum }_{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{{\sum }_{i=1}^N(x_i-\overline{x}{)}^2}\sqrt{{\sum }_{i=1}^N(y_i-\overline{y}{)}^2}},$$

Full text

Sun: corona — Sun: coronal mass ejections (CMEs) — techniques: genetic algorithm

††journal: ApJ

1 Introduction

Coronal mass ejections (CMEs) are usually observed by coronagraph (Lyot, 1939). The observation records radiation of Thomson scattering (Minnaert, 1930; van de Hulst, 1950; Billings, 1966; Vourlidas, & Howard, 2006; Howard & Tappin, 2009; Inhester, 2015) which is produced by interaction between radiation from the photosphere and free electrons inside CMEs. Information of the electron locations along line of sight (LOS) is hidden to the observer after the observation (Howard et al., 2008; Temmer et al., 2009).

Those hidden information can be restored by different methods of CME reconstruction (Mierla et al., 2010; Thernisien et al., 2011). We classify the various methods into three categories as follows.

(1) Methods to reconstruct point features of CMEs. The three dimensional (3D) coordinates of a CME feature can be calculated if the feature can be seen from different view points like Solar Terrestrial Relations Observatory (STEREO) A and B (Kaiser et al., 2008; Howard et al., 2008). For example, Mierla et al. (2008) calculated 3D position of the feature using height-time diagrams from COR1 onboard STEREO-A and B based on epipolar geometry (Inhester, 2006). Liu et al. (2010a, b, 2014) obtained the 3D position of CME feature from corona to helioshpere using time-elongation map from COR2, HI1 and HI2 onboard STEREO-A and B.

(2) Methods to reconstruct outline of CMEs. A global outline can be obtained from the CME boundary measurements in the coronagraph images. For example, Pizzo & Biesecker (2004) using the CME boundaries to construct a series of quadrilaterals to approximate the three dimensional CME outline based on the synthetic STEREO images. de Koning et al. (2009) applied this method to reconstruct the CME outline using coronagraph observations from STEREO. Thernisien et al. (2006, 2009) employed a forward model of Graduated Cylindrical Shell (GCS) to reconstruct flux rope-like CMEs using data from either single viewpoint or multiple viewpoints. Byrne et al. (2010) employed an elliptical tie-pointing technique to reconstruct a full CME front in 3D space. Feng et al. (2012, 2013) developed a 3D mask fitting reconstruction method using coronagraph images from three viewpoints to obtain the 3D morphology of a CME.

(3) Methods to reconstruct density inside CMEs. For example, Moran & Davila (2004) derived 3D position of electrons along LOS corresponding to every CME pixel through polarization analysis of single-view images. Dai et al. (2014) suggested a classification of ambiguity in this method to improve the reliability of the reconstruction. Antunes et al. (2009) reconstructed 3D CME electron density through combination of inversion using PIXON method (Puetter et al., 2005) and forward modeling (Thernisien et al., 2006, 2009). Frazin et al. (2009) developed a reconstruction method based on level-set (Jacob et al., 2006) algorithm using synthetic CME images produced from 2D slice of magnetohydrodynamic (MHD) simulation (Manchester et al., 2008). Hutton, & Morgan (2017) developed Automated CME Triangulation (ACT) to find the most probable location where the CME passes at height of 5 $R_{\odot}$. ACT is effective under different observations from three or two viewpoints over a sliding window of five hours.

Information of CME structure obtained from 3D reconstruction gradually increases by methods from category (1) to (3). For category (3), the reconstructed CME is more close to the real CME which is a collective of inhomogeneous magnetized plasmas. On the other hand, ambiguities of these methods in category (3) are also larger than those in category (1) and (2). In the current work, we develop a new method to reconstruct the 3D electron density of CME based on genetic algorithm (GA) (Holland, 1992; He et al., 2011), namely genetic reconstruction method (GRM). According to the classification mentioned above, GRM belongs to category (3).

The rest of this paper is organized as follows. We introduce GRM in Section 2, and present application of the method in Section 3.1 and 3.2 for model and observed CMEs. Finally, conclusions are given in Section 4.

2 Genetic Reconstruction Method

Before method description in this section, we briefly explain the approach to reconstruct CME electron density using GRM. Our approach is to randomly distribute test electrons through a 3D reconstruction space, with their populations within voxels used to calculate line of sight integrations that create synthetic observations. These can be compared to the true observations, giving a goodness of fit to the electron distribution. We can randomly distribute electrons in the reconstruction space to create many different test models, and their goodness of fit to the data sets a scoring to the models. This sets a basis for a genetic algorithm approach to improve the fitness of synthetic brightness to observation. According to the scoring of fitness, we select the better electron distribution to the next iteration. When the iteration process meet the termination condition, we obtain the final reconstruction of electron distribution. Process of the reconstruction with five steps is shown in Figure 1.

2.1 A general description for CME reconstruction

Observed brightness, $B_{obs}$, of each CME pixel is the corresponding LOS integral of Thomson Scattering including measurement noise $\eta$:

[TABLE]

Here, $N_{e}(\theta)$ is the electron density on LOS where the angle between the plane of sky (POS) and the line connecting the electron to the Sun center has value $\theta$. The analytic function $B_{e}(\theta)$ (Quémerais, & Lamy, 2002) is the value of Thomson scattering for a single electron at location $\theta$. $\theta_{1}$ and $\theta_{2}$ are LOS boundaries constrained by the CME shape.

CME reconstruction tries to obtain the vector $N_{e}(\theta)$ when $B_{obs}$ is observable. Methods in category (1) focus on the localized feature who contributes most of the brightness along the LOS. Methods in category (2) focus on the electron density who contributes the brightness of the observed CME boundary. Methods in category (3) focus on the electron density along the LOS between $\theta_{1}$ and $\theta_{2}$. We summarize part of the challenges and limitations of CME reconstruction (Frazin et al., 2009; Mierla et al., 2010; Thernisien et al., 2011) as follows:

(1) Equation 1 is highly underdetermined. There are multiple possible solutions of Equation 1 if constraints or a priori information are not enough. Possible distributions of the $N_{e}(\theta)$ may generate the same observed brightness $B_{obs}$.

(2) Mix of CME and background corona. The CME images always contain brightness of background corona because of the optically thin nature of Equation 1, $B_{obs}$ = $B_{cme}$ + $B_{bg}$. $B_{bg}$ can be removed by subtraction of pre-CME image from CME image (Vourlidas et al., 2010). Situation becomes complex when the background corona contains dynamic structures like streamers, jets or other CMEs. It may be difficult to do a clean subtraction. Some advanced techniques are developed to make more effective background subtraction by more sophisticated manners (Morgan et al., 2012; Morgan, 2015).

(3) Measurement noise in the coronagraph images. The noise, $\eta$, of any measurements is unavoidable. As pointed out by Thernisien et al. (2011), the presence of noise in the CME images may lead to unstable solutions of Equation 1.

(4) Incomplete observation caused by occulter of coronagraph. The CME may be obscured by occulter of coronagraph especially when the CME has not propagated far away from the Sun. And for different viewpoints, occulters may obscure different parts of the CME.

(5) The number of simultaneous observation viewpoints. Simultaneous observations from different viewpoints play important role to most of the reconstruction methods classified in Section 1. Inverse methods that try to restore $N_{e}(\theta)$ usually need data from many view directions, like tomography. With two or three views, like Solar and Heliospheric Observatory (SOHO) (Domingo et al., 1995; Brueckner et al., 1995) and STEREO, the inverse problem is extremely ill-posed. Solution of Equation 1 is not unique.

(6) Separation angle between view directions. For reconstructions using data from two views, separation angle of $0^{\circ}$ and $180^{\circ}$ is unsatisfactory while $90^{\circ}$ is ideal. STEREO A and B orbit the Sun with different speed and the angle between them varies about $44^{\circ}$ per year.

We should keep these limitations in mind through out the process of CME reconstruction.

CME reconstruction using inverse methods can be stated as an optimization problem:

[TABLE]

where $cf$ is called cost function. $p$ is a positive number which defines the error norm of the observations. Traditionally, p=2 and the regularization term also depends on $N_{e}^{2}$ so that minimizing Equation 2 results in a least-squared problem for which many solvers exist. According to the limitations of CME reconstruction that mentioned above, the optimization is highly ill-posed unless the regularization expression is chosen properly so that it yields different penalties for the many solutions in the null space of the data term.

As shown in Equation 2, an appropriate regularization needs to be added to obtain more stable solution to deal with the ill-posed problem. Since many solutions are possible to eliminate the data term in Equation 2, (in case that $\eta$=0, else the data $B_{obs}$ is probably inconsistent and there is no exact solution but there are many for which the delta error is of the order of $\eta$) the solution which is finally returned by any procedure depends on the error norm and the regularization we choose. We emphasize that this ambiguity lies in the heart of the CME reconstruction problem and is unavoidable because two views are not sufficient to determine a unique CME density distribution. Our goal is to test whether a certain combination of error norm p and regularization may yield more realistic solutions than others. Since conventional solvers only work for p=2 and $N_{e}^{2}$ dependent regularization, we decided to use GRM to find the minimum of Equation 2 so that we maintain a maximum of flexibility to choose p and the regularization expression.

In order to show the ability of GRM itself for CME reconstruction, the regularization has not been added in the current work and we choose p=1:

[TABLE]

The L1 norm is often favoured over least-squared solutions in image processing because it attributes less weight to measurement outliers. A systematic test of different regularization operators will be made in the next paper.

2.2 Calculations of electron number and white light brightness

The optimization is initiated by calculating electron number for each LOS and redistributing the electron density $N_{e}(\theta)$ in a discrete form along each LOS. Thomson scattering can then be calculated for the redistributed electrons and the cost function is ready for the optimization.

In order to explain the process of initiation, a cube of plasmas is used to mimic CME as shown in the middle panel of Figure 2. In the first and third panel, the simulated white light brightness is produced through Thomson scattering mechanism in the field of view (FOV) of STEREO A and B, respectively. There is no need to subtract the background corona since we only simulate the brightness of model CME itself. For brightness of real CME in section 3.2, we rescale the CME images to 128$\times$128 pixels and make a 3$\times$3 smoothing to minimize the influence of noise which is inherent in the measurements. A pre-CME background image is subtracted to get the excess brightness of CME (Vourlidas et al., 2010).

It is necessary to split the brightness for simplifying the reconstruction. As shown in Figure 2, the brightness is equally split up into ribbons with white lines for FOV A and B. For example, the corresponding part of Ribbon-A in FOV A is Ribbon-B in FOV B displayed by black lines. Pixels in this pair of ribbons belong to the same set of epipolar planes (Inhester, 2006) between two epipolar planes indicated in the middle panel of Figure 2. Process for finding ribbon B is carried out by projecting LOS of each pixel in ribbon A into FOV B. LOS A that is projected into FOV B should intersect with the CME boundary at least at two cross points unless the LOS is tangent to the CME boundary. Then the projection of LOS A and LOS A itself become finite line segments between those two cross points, not the infinite LOS from the observer to infinity. We assume that all the LOSs are parallel to each other in the FOV of a given telescope. Reconstruction can be carried out for each pair of the ribbons to obtain the electron density in the corresponding set of epipolar planes.

We use $i$ to index pixels along epipolar lines in FOV A, $j$ to index computation grid along LOS A, $k$ to index epipolar planes. In each epipolar plane, coordinate axes are along the view directions and usually non orthogonal. Pixels in FOV A can be indexed by $(i,k)$. Integration along LOS A is equivalent to summation of index $j$.

(1) First of all, electron number along LOS of corresponding pixel of the model CME brightness image in FOV A can be calculated as:

[TABLE]

where $B_{obs}(i,k)$ is brightness value of the pixel $(i,k)$. Thomson Scattering of a single electron, $B_{e}(\theta(i,k))$, can be calculated as

[TABLE]

$\theta$ is traditionally assumed to be zero (Vourlidas et al., 2010) when the electron location along the LOS is unknown. z is the distance from the observer to the scattering location along the LOS. Three dimensional reconstruction , like Moran & Davila (2004); Dai et al. (2015), can obtain the $\theta$ value of each pixel to get a more precise calculation of the electron number as done in this paper. $\sigma$ is the Thomson scattering cross section of an electron and $u$ is the limb darkening coefficient. $A$, $B$, $C$ and $D$ are the van de Hulst coefficients (van de Hulst, 1950).

In this work, electron number is always calculated from observation 1 as shown in Figure 2. In principle, observation 1 can be anyone of STEREO A, B or SOHO. Observation 2 is not restricted to STEREO-B but also possible for SOHO as shown in section 3. This enable us to carry out the reconstruction without data of STEREO B.

(2) Redistribution of the electrons of each pixel along the corresponding LOS. As shown in the top right panel of Figure 1, series of LOS are generated from Brightness-A and the extent of LOS is restricted by the length of Brightness-B. Then we get the LOS segment. Each LOS segment is uniformly separated into $N_{j}$=64 computational grids. The voxels on the computational grids are initially populated with a random distribution of electrons. An array of random decimal, $L(i,j,k)$, are used to represent the electron distribution.

We calculate electron number in the $j_{th}$ voxel on the grid point along the LOS segment by the following equation:

[TABLE]

(3) Calculation of the Thomson scattering for the redistributed electrons. Equation 7 converts the electron density $N_{v}(i,j,k)$ into white light emission $B_{av}(i,j,k)$:

[TABLE]

where

[TABLE]

Then we can do the summation for $j$ to get the total brightness of a CME pixel $(i,k)$ in FOV A:

[TABLE]

The total brightness of a CME pixel $(j,k)$ in FOV B can also be calculated by a similar summation for $i$:

[TABLE]

using Equation 7 and 8 when values of $\theta_{a}$ and the van de Hulst coefficients in FOV A are replaced by values of $\theta_{b}$ and the van de Hulst coefficients in FOV B. We emphasize that each $LOS\_B$ contains electrons from different $LOS\_A$. The values of $\theta$ and the van de Hulst coefficients should be calculated for Equation 7 and 8 according to the electron locations in the coordinate system of FOV A and FOV B, respectively.

2.3 CME reconstruction using genetic algorithm

Now, we can calculate the cost function for pixel $(i,k)$ in FOV A using Equation 3:

[TABLE]

and the cost function for pixel $(j,k)$ in FOV B:

[TABLE]

Once the cost function is ready, the process of CME reconstruction can be described in the following five steps which are typical for GA as shown in the top left panel of Figure 1.

(1) Population Initiation. The redistributed electrons along all of the LOS segments constitute one of the individual of GA. The same process of random redistribution is carried out for 200 times to produce the initial population $L_{q}(i,j,k)$, $q=1,2,3,...,N_{q},N_{q}=200$. This is the beginning of GA as shown in the top left panel of Figure 1.

(2) Fitness calculation. Fitness function used in GA is inversely proportional to the cost function:

[TABLE]

We show an example of the fitness of $N_{q}$ individuals in the top panel of Figure 3, from smallest to largest. In order to avoid local solution which may be produced by the steep distribution of fitness as shown in the top panel of Figure 3, we need to rescale the fitness function. Firstly, the values of fitness is changed to be linearly increasing by

[TABLE]

where $sort$ is the ranking function implemented in IDL. The sorted array of fitness $Fit[sort(Fit)]$ is incremental as shown in the top panel of Figure 3. The IDL function $findgen(N_{q})$ creates an arithmetic progression [0,1,2,…,$N_{q}-1$]. Then the values of fitness are transformed to [1,2,3,…,$N_{q}$] by Equation 14 as plotted in the middle panel of Figure 3. Finally, the IDL function $exp$ are used to transform the fitness values into natural exponential distribution by

[TABLE]

and

[TABLE]

where $g=1,2,3,...,N_{g},N_{g}=100$ is the total number of generation in the genetic evolution. As shown in the bottom panel of Figure 3, the fitness distribution becomes steeper from generation 40 to 60. It means that the algorithm can avoid local solution at the preliminary stage of genetic evolution and keep it more convergent at the later stage.

(3) Selection operator is employed to obtain the optimal electron distribution. In order to further simplify the reconstruction, we equally divide the brightness ribbon into 20 parts for FOV A and B respectively as shown in Figure 4. For each part of the brightness ribbon in FOV B, a group of electron distribution $L_{q}$ along the LOS segments inside a column which is marked by the purple lines for instance, contribute to the brightness of Thomson scattering. Mean value of fitness is already calculated for pixels inside each part of the brightness ribbon using Equation 16 in step (2). We randomly choose two candidates from $N_{q}$ candidates of $L_{q}$ group and compare their mean fitness. The candidate with larger fitness will be chosen and passed on to the next generation. Such tournament is repeated $N_{q}$ times for each part of brightness ribbon in FOV B to produce a new population. based on this new population, brightness in FOV A can be updated and the mean fitness of the whole brightness ribbon in FOV A can be calculated by Equation 15 for each of the $N_{q}$ individuals. Such mean fitness of the whole brightness ribbon is also calculated for FOV B based on the new population. Ten individuals with worst fitness in FOV B are replaced by ten individuals with best fitness in FOV A. This replacement optimizes the population for FOV A after the tournament selection for FOV B.

(4) In this work, we set the maximum number of generation to 100. If the genetic evolution meets this termination condition, the GA will stop and an optimal solution will be chosen to be the final reconstruction of CME. If the number of generation is less than 100, the algorithm will go to step (5).

(5) Crossover and mutation operators are used to update the population $L_{q}$. The basic concept of these operations is presented in the bottom panel of Figure 1. The crossover operator exchanges the electron distribution between two randomly selected individuals while the mutation operator changes the distribution inside one individual. Formula 17 and 18 describe these operators:

[TABLE]

where $L_{m}(i,j,k)\in CME_{m}$, $L_{n}(i,j,k)\in CME_{n}$ and $m,n\in[1,N_{q}],m\not=n$.

[TABLE]

where $D_{r}$ is a random decimal which can be positive or negative. In order to simplify the reconstruction, crossover and mutation operators are applied to all of the $L_{q}$ inside one of the 20$\times$20 areas with shape of parallelogram at each time.

After the genetic operation, we get the new CME population from parents to children. Probabilities of crossover and mutation decrease progressively to keep the optimization being global and convergence.

Process from (2) to (5) are repeated until the number of generation is equal to $N_{g}$.

In the following sections, GRM is firstly applied to a set of model CMEs with different directions and shapes. Then the method is employed to reconstruct a CME observed by coronagraphs of SOHO and STEREO.

3 Application

In order to test the GRM method in a more realistic way, we use a set of analytic GL98 (Gibson, & Low, 1998) model CMEs with various orientations and shapes, as well as occultation to produce synthetic coronagraph data instead of the cube of plasmas as used in Figure 2 at section 2.2. The GL98 model constructs a flux rope to create typical three-part CME structure. The flux rope is an analytical solution of magnetohydrostatic (MHS) equation

[TABLE]

and solenoidal condition of magnetic field, $\nabla\cdot\bm{B}=0$. The explicit solution of Equation 19 can be derived from solution of

[TABLE]

by applying a mathematical stretching transformation $r\rightarrow r-a$ to an axisymmetric sphere of magnetic flux (Chandrasekhar, 1956; Chandrasekhar, & Prendergast, 1956) with radius $r_{0}$. The sphere center is located at $r_{1}$ from solar center. After the transformation, the magnetic flux appears to be a tear drop. Then the magnetic field of GL98 model can be expressed by a Bessel function and a free parameter, $a_{1}$ (Prendergast, 1956; Lites et al., 1995). As pointed out by Gibson, & Low (1998), the solution of Equation 19 can even be used to solve MHD equation

[TABLE]

under self-similar theory (Low, 1984). GL98 model is implemented in Space Weather Modeling Framework (SWMF) (Tóth et al., 2005, 2012) and is successfully used in studies of CME simulation (Manchester et al., 2004a, b, 2014a, 2014b; Lugaz et al., 2005; Jin et al., 2016, 2017a). Recently, this model is applied in a user-friendly tool named Eruptive Event Generator using Gibson-Low configuration (EEGGL) (Jin et al., 2017b; Borovikov et al., 2017) which has been transitioned to the Community Coordinated Modeling Center (CCMC).

Finally, we obtain the distribution of plasmas density frozen into magnetic field:

[TABLE]

We change position of the sphere center, $r_{0}$, $r_{1}$, $a$ and $a_{1}$ to obtain model CMEs with different orientation, size, shape and strength of magnetic field. Since the distribution of plasmas density is known, the synthetic CME white light images can be produced through the Thomson scattering mechanism mentioned in section 2.2.

We apply GRM to these synthetic white light images to reconstruct the model CMEs. Comparisons of brightness and electron density between models and reconstructed CMEs are presented to show advantages and limitations of GRM in section 3.1. After systematic comparisons between model and reconstructed CME, GRM is applied to real observations from SOHO and STEREO in section 3.2.

3.1 Application for model CMEs

We simulate CME observations using the same parameters of STEREO/COR2 and SOHO/C3 including satellite position, view direction, size of FOV and spatial resolution on February 15th, 2013. Positions of the Earth (SOHO), STEREO A and B in our simulation are shown in Figure 5.

For simulations of two view observations from STEREO A and B as observation 1 and 2, the separation angle is set to $90^{\circ}$ to test interfering factors of reconstruction like longitude and latitude of CME, whether or not to use occultation, whether or not the CME is halo. For STEREO A and SOHO as observation 1 and 2, the separation angle is set to $135^{\circ}$. We add data of COR2 from STEREO B as the third view to improve the performance of GRM.

Central longitude of the model CME is set to $0^{\circ}$, $45^{\circ}$, $90^{\circ}$, $135^{\circ}$ and $180^{\circ}$ respectively as shown in Figure 5. Central latitude of the model CME is routinely set to $0^{\circ}$ and $30^{\circ}$ for each value of longitude. Latitude of $15^{\circ}$ and $45^{\circ}$ are also set to the model CME with longitude of $0^{\circ}$ to test latitude dependence of GRM more carefully. For convenience, we mark the CME with (longitude,latitude) for specific values of longitude and latitude. For example, CME with longitude of $0^{\circ}$ and latitude of $0^{\circ}$ can be called (0,0) CME.

In order to quantitatively evaluate the similarity between the model and reconstructed CME, we calculate mean relative absolute deviation ($MRAD$)

[TABLE]

and linear Pearson correlation coefficient ($PCC$) (Li et al., 2018)

[TABLE]

where $x$ and $y$ are values of brightness or electron density of the model and reconstructed CME, $N$ is the number of CME voxel. If y is similar to x, value of $MRAD$ should be close to 0 and value of $PCC$ should be close to 1. The upper limit for $MRAD$ is set to 1.0.

We show similarity between the model and reconstructed CMEs qualitatively with images and quantitatively with $PCC$ in Figure 7 and the following figures. For convenience, we name these figures “Figure of Similarity” (FOS). We run GRM three times for each model CME to prove consistency of the method as shown in FOS which is labelled with “1st”, “2nd” and “3rd”. Because of the ill-posed property of CME reconstruction and inherent randomness of GRM, reconstruction results including Thomson scattering brightness and electron density of the same model CME should be slightly different between each run.

For each run of GRM, we not only show the 3D distribution of electron density in the first column of last three rows of FOS but also the 2D distribution of density averaged along directions of x, y and z axis as shown in the last three columns of these rows. For example, the density averaged along direction of z forms the 2D distribution in xy plane as shown in the second column which is labelled with “xy” in the lower left corner. Performance of GRM is different for a specific CME model on xy, xz and yz plane. Average density in these planes enable us to make more comprehensive comparisons between the model and reconstructed CME.

For comparison of brightness, we display the original brightness of GL98 model CME in the first column of the first and second row of FOS. The brightness reconstructed by GRM at generation 100 is labelled with “1st”, “2nd” and “3rd” for different runs of GRM. We plot $PCC$ of brightness from generation 0 to 100 in the third row to show the process of genetic evolution. The plots of $PCC$ at “1st”, “2nd” and “3rd” runs for reconstruction 1 are shown in the first column with colour of red, green and blue. Plots of $PCC$ for reconstruction 2 are shown in the second column. $PCC$ becomes larger from generation 0 to 100, which means that the optimization is convergent.

For comparison of electron density, the electron density in 3D space as shown in FOS is the average inside cubes with width of 2.0 $R_{sun}$. The original density of GL98 model CME in 3D space, xy, xz and yz plane are firstly presented in the fourth row of FOS. Then reconstructed density of the “1st”, “2nd” and “3rd” runs of GRM are displayed in the last three rows of FOS. Heliocentric Earth Equatorial (HEEQ) (Thompson, 2006) coordinate system is used to illustrate the CME density in FOS. The corresponding heliographic coordinates is Stonyhurst. In Figure 7, longitude and latitude of the model CME are both equal to $0^{\circ}$. It means that direction of the (0,0) CME is along z axis. In order to simulate the brightness of (135,30) CME, longitude and latitude of the model CME are changed to $135^{\circ}$ and $30^{\circ}$ respectively. Based on the simulated brightness, we reconstruct the electron density of CME using GRM. Then we transform the model and reconstructed CME by $-135^{\circ}$ and $-30^{\circ}$ for longitude and latitude respectively. After this kind of transformation, the density of model and reconstructed CME in Figure 11 can be viewed under the same perspective as in Figure 7. Similar transformations are applied to the rest of FOS. $PCC$ of electron density distribution between model and reconstructed CME labelled in FOS are summarized in Table 1. In order to give a more complete picture of the GRM performance, values of $MRAD$ along with $PCC$ are shown in the left columns of Table 1.

We further discuss the interfering factors of GRM based on FOS and Table 1 as follows.

(1) As shown in Table 1 and FOS, values of $PCC$ in xz plane is usually lower then those in xy and yz plane when latitude is equal to $0^{\circ}$. Value of $PCC$ in xz plane is even less than 0.6 at the second run for (90,0) CME. An obviously improvement of $PCC$ in xz plane can be seen when latitude grows from $0^{\circ}$ to $30^{\circ}$ for CMEs with all values of longitude. Comparison between (0,0) CME in Figure 7 and (0,30) CME in Figure 8 validate such improvement as an example. For latitude from $30^{\circ}$ to $45^{\circ}$, the improvement of $PCC$ stops as shown for the $0^{\circ}$ longitude CME.

$PCC$ of the (135,30) CME are the best among all of the reconstructed CMEs. For this CME, $PCC$ in 3D space are larger than 0.8 while $PCC$ in the xy, xz and yz plane are even larger than 0.9.

(2) The (45,0) CME appears full halo in FOV B as a back side event and the (135,0) CME appears full halo in FOV A as a front side event. The (45,30) and (135,30) CME become partial halo when latitude is changed to $30^{\circ}$.

The halo CMEs are obscured by occulter of the coronagraph as shown in Figure 9, 10 and 12. For halo CMEs, influence of occulter is more obvious than limb CMEs. Central part of the halo CME is hidden behind the occulter. For GRM reconstruction, situations in Figure 9 and 10 is different, although these two CMEs are both full halo. The (45,0) CME is full halo for observation 2 as shown in Figure 9. On the other hand, the CME is almost not affected by the occulter since it is a limb event for observation 1. We obtain a relatively complete initial electron number of CME for GRM because the electron number is calculated from brightness in observation 1. In contrast, the initial electron number is obviously incomplete in Figure 10 for the (135,0) CME.

As a result, brightness of reconstruction 2 is obviously larger than that of observation 2 for the (45,0) CME as shown in Figure 9 because of the completeness of CME in observation 1 and incompleteness in observation 2. Situation is opposite for the (135,0) CME as shown in Figure 10.

(3) As shown in Figure 14, the (135,30) CME is reconstructed using simulation data from view of SOHO/C3 as observation 2. Separation angle between STEREO A and SOHO is $135^{\circ}$. Comparing to the orthogonal coordinates with $90^{\circ}$ separation angle, $135^{\circ}$ separation angle is not ideal for GRM reconstruction. In order to improve the performance of GRM under non orthogonal coordinate system, view of STEREO B is added into the reconstruction as observation 3 in Figure 15. Outline of the CME in observation 3 constrain the 3D electron distribution as shown in Figure 15 comparing to Figure 14. In the current work, we just use the outline of CME in observation 3 as a constraint. Such constraint makes the boundary of LOS more accurate. In the future, we may calculate the fitness function of brightness in observation 3 together with those of observation 1 and 2 to optimize the reconstruction of electron distribution.

(4) For the original synthetic observations of model CME, there is no noise. By contrary, noise is inevitable in real observations. We add noise to the synthetic images of (135,30) CME as shown in Figure 13. This is a critical test to evaluate the performance of GRM when the method is applied to real observations from coronagraphs. The noise image is produced by subtracting the background for a coronagraph image which is clear of CME. In a perfect subtraction of noiseless image, pixel values should be zero since there is no CME brightness. However, the brightness values usually deviate from zero for real observations as shown in the top panel of Figure 6. This deviation can approximate the normal noise level in real observations. In the second and third panel of Figure 6, average values of noise, $B_{n}$ and total brightness, $B_{t}$ are calculated within rings of 0.5 $R_{\odot}$ width. It is obvious that both of $B_{n}$ and $B_{t}$ decrease rapidly outward the Sun. However, the relative noise $B_{n}$/$B_{t}$ increase with height.

This normal noise level is added to the synthetic CME brightness as shown in Figure 13 and Table 1. Results of noise level multiplied by 10 and 100 are also shown in Table 1, respectively. We can see that difference of $PCC$ values between Figure 13 with normal noise and Figure 12 without noise is not obvious. In the case of normal noise level, values of $PCC$ are still larger than 0.8 and values of $MRAD$ are still less then 0.8. Difference becomes larger when the noise is multiplied by 10 and 100. An obvious decline of GRM performance can be seen in the case of noise multiplied by 100. In this case, values of $PCC$ are even less than 0.3 and values of $MRAD$ are close to 1.0. It means that the correlation between model and reconstructed CME is very low and GRM is no longer effective. On the other hand, GRM is still feasible under normal noise level which is common in the real CME observations.

3.2 Application for observed CME

In this section, we show the application of GRM for a real CME event observed by STEREO and SOHO on September 30th, 2013. The CME is observed by COR2 on STEREO A and B at 00:08 and by C3 on SOHO at 00:04. Difference of the observation time is four minutes between COR2 and C3. The general case is that difference of observation timings usually exist between COR2 and C3. What we can do is to choose the nearest observation time of C3 to match the COR2 observation. Coordinate system of this GRM reconstruction is non orthogonal because separation angle between STEREO A and B is $73.5^{\circ}$ while separation angle between STEREO A and SOHO is $147.1^{\circ}$.

In Figure 16, the CME is reconstructed from data of STEREO A, B and SOHO as observation 1, 2 and 3. In Figure 17, we reconstruct the CME using STEREO A, SOHO and B as observation 1, 2 and 3. For this real CME, $PCC$ of brightness can still be calculated between real observation and reconstruction. Plots of the brightness $PCC$ from generation 0 to 100 still keep convergence as shown in the FOS. However, $PCC$ of electron distribution can not be calculated between real and reconstructed CME because we do not know the 3D electron distribution of the real CME. Instead, $PCC$ of electron distribution can be calculated for the “2nd” and “3rd” reconstruction by comparing to the “1st” reconstruction. These values of $PCC$ are also summarized in Table 1 marked with “ABS” and “ASB” for observation 2 using data from STEREO B and SOHO, respectively. Because we use density distribution in the “1st” reconstruction as the reference (“model CME”), values of $PCC$ are equal to 1.0 and values of $MRAD$ are equal to 0.0 for the “1st” reconstruction as shown in Table 1. For the “2nd” and “3rd” reconstruction, values of $PCC$ are larger than 0.78 and values of $MRAD$ are less than 0.4. These results show us the stability of reconstructions in different runs of GRM.

GRM reconstructs this real CME using data from STEREO B and SOHO as observation 2 respectively. The CME is partial halo in these two FOVs. Thus we can not detect the complete brightness of CME core which is obscured by the occulter of coronagraph. On the other hand, the CME brightness is relatively complete in FOV A as observation 1 because the CME is a limb event in this FOV. A similar situation can be seen in Figure 9 for reconstruction of the model (45,0) CME with occulter. This model CME is a limb event in observation 1 and is a halo event in observation 2. Although the CME brightness is complete in observation 1, the reconstruction is still not correct because of the incomplete brightness in observation 2. The structure of reconstructed CME core is evidently destroyed as shown in Figure 9 comparing to the model CME. This may also be the reason that we can not find a typical dense core in the reconstruction of real CME as shown in Figure 16 and 17. Besides the influence of occulter, the CME itself may not be a typical three part event. The dense core may become indistinct during the outward expansion.

As pointed out by Morgan (2015), accuracy of reconstructed density depends on reliability of the input real data besides the reconstruction method. Calibration and background subtraction play important roles to obtain real CME brightness without contamination from vignetting, stray-light, F corona and other static structures in K corona like streamers and coronal holes, et al..

In this work, we download the level 0.5 fits files of coronagraph images and process them to level 1.0 using standard calibration routines from the Solarsoft library. For total brightness ($B_{t}$) observed by coronagraph C3 on board SOHO/LASCO, we apply the Solarsoft routine reduce_level_1.pro to make calibrations for dark current, flat field, stray light, distortion, vignetting, photometry, corrected time and position, et al.. For $B_{t}$ observed by coronagraph COR2 on board STEREO/SECCHI, we apply secchi_prep.pro to make calibrations for subtracting the CCD bias, multiplying by the calibration factor and vignetting function, and dividing by the exposure time, et al..

After the standard calibration, $B_{t}$ contains a combination of K corona brightness ($B_{k}$) and F corona brightness ($B_{f}$). Morgan (2015) developed a method to separate $B_{k}$ and $B_{f}$ for C2 observations with FOV between 2.2 and 6.0 $R_{\odot}$. In this method, an in-flight calibration using star brightness is employed to improve the standard calibration for polarized brightness ($B_{p}$) of C2 which is usually obtained once a day. Series of calibrated $B_{p}$ during half period of solar rotation is subtracted from the corresponding $B_{t}$ to produce $B_{f}$ approximation. Finally, $B_{k}=B_{t}-B_{f}$. For C3 observations with FOV between 3.7 and 30 $R_{\odot}$, we may use this method to obtain cleaner corona brightness by removing $B_{f}$ from $B_{t}$ in the future.

$B_{k}$ contains brightness of Thomson scattering from coronal structure like streamers, coronal holes and CME, et al.. For CME reconstruction, we need to remove the background including the relatively static components like streamers, coronal holes. Subtraction of a pre-CME image observed before the CME is a standard approach as what we have done in this work. After such subtraction, the excess brightness of CME is available. Morgan et al. (2012); Morgan (2015) developed a useful Dynamic Separation Technique (DST) to make the subtraction more effective. DST is valid under assumption that the quiescent structures like streamers are basically smooth in radial direction and slowly evolves in time while dynamic structures like CMEs are not radially smooth and evolve faster. Then DST is able to separate $B_{k}$ into quiescent and dynamic component using iterative deconvolution in the radial and time dimensions. In the future improvement, we will employ DST to subtract the background corona for CME reconstruction. Morgan (2015) presented a method to make cross calibration between C2 and COR2 which is necessary for electron density reconstruction using observations from different coronagraphs (Morgan, 2019). Such method will also be considered in the future work.

4 Conclusions

In this study, genetic reconstruction method (GRM) is used to reconstruct 3D distribution of CME electrons. We give a general description of CME reconstruction at first. Ill-posed property of the CME reconstruction is seriously pointed out. A set of analytic GL98 model CMEs with different orientations and shapes are employed to produce synthetic CME images without noise for the genetic reconstruction. Random noise is artificially added to some synthetic CME images to imitate the measurement noise which is unavoidable in the real CME observations. Since the electron distribution is known for the model CMEs, we make comparison based on $PCC$ for both of the electron distribution and its corresponding Thomson Scattering brightness between the model and reconstructed CME. $PCC$ of brightness from 0 to 100 generation presents the convergence of GRM for all of the model CMEs. $PCC$ of electron distribution in 3D space and 2D planes show us the ability of GRM to obtain stable and reasonable CME reconstruction. Values of $MRAD$ of electron distribution in 3D space validate the results of $PCC$.

Based on the comparison between model and reconstructed CME presented in FOS and Table 1, we understand more in depth about the advantages and limitations of GRM. Performance of GRM depends on the longitude and latitude of model CME as well as the completeness of observation and separation angle between points of view. A more reliable reconstruction can be obtained if: (1) The model CME in coronagraph is not much obscured by the occulter; (2) Central latitude of CME is about $30^{\circ}$; (3) Separation angle is about $90^{\circ}$. (4) Data of the third observation can be added into reconstruction.

Then the method is applied to real coronagraph data from STEREO A, B and SOHO. We compare the reconstructed brightness with observation to show the convergence of GRM. Comparisons of electron distribution between reconstructions from different runs tell us that the results of GRM reconstruction are stable. As pointed out in section 3.2, we should employ more effective techniques to calibrate the coronagraph data and make the background subtraction in the future work.

The purpose of this paper is to demonstrate how GRM could be used to find a solution. Because of the ill-posed nature of CME reconstruction using data from only two or three view points and the random nature of GRM, results of GRM for the same model CME are different between three runs. This just illustrates the range of solutions which are possible for the unregularized problem. Regularization as in Equation 2 is helpful to further stabilize GRM reconstruction and may also mitigate the ill-posed characteristic towards a unique solution. How realistic the solution is then depends on how reasonable the regularization operator is. The GRM may allow least restrictions in the choice of the regularization. In the future work, we may test different forms of regularization to find suitable constraints for the CME reconstruction.

We appreciate the anonymous reviewer for valuable comments which help us greatly improve the manuscript. We acknowledge SOHO and STEREO consortia for usage of the data. This work is jointly supported by the National Natural Science Foundation of China (NSFC) through grants 11603040, 11973058, 11473040, 11705210, 11873060 and U1831121.