SuNeRF: 3D Reconstruction of the Solar EUV Corona Using Neural Radiance Fields

To validate our method, we leverage a magnetohydrodynamic simulation, developed by Predictive Science Inc., of the solar corona used to forecast the state of the solar atmosphere prior to the 2019 July 2 total solar eclipse (see the website ¹⁰ and Boe et al. 2021, for more details). The data set consists of a single snapshot from which we extract full-disk images from 256 viewpoints at a resolution of 1024 × 1024 pixels (263 pixel^–1) in three different wavelengths (171, 193, and 211 Å). We preprocess the simulated data analogously to the observational data introduced in Section 2.2.

From 2010 to 2014, three satellites capable of imaging in similar EUV channels were orbiting the Sun, namely, the SDO equipped with the AIA instrument (Lemen et al. 2011) and the twin STEREO spacecraft, each equipped with an EUVI (Wülser et al. 2004). All three satellites orbit in the ecliptic plane, with SDO in geosynchronous orbit and the STEREO satellites progressively advancing (STEREO-A) and trailing (STEREO-B) Earth. When the separation angle between SDO and each of the STEREO spacecraft exceeded 90°, it allowed a view of the full Sun at one instant for the first time.

An important issue in the context of this project is the design, performance, and operational differences among the various instruments. We apply Instrument-To-Instrument translation (ITI; Jarolim et al. 2024) to translate full-disk observations from STEREO to the SDO domain. With ITI, we can account for differences in resolution, calibration, and filter bands between the individual instruments. The AIA and EUVI instruments have common observing wavelengths (171 and 304 Å). For the wavelengths 195 and 284 Å, we use ITI to approximate the 193 and 211 Å channels of AIA (see Jarolim et al. 2024). Although labeled as single wavelength channels, the SDO/AIA and STEREO/EUVI spectral channels actually result from spectral integration as specified by the spectral response of the instruments and wavelength filters (Szenicer et al. 2019). The multichannel translation of ITI is a prerequirement to have an appropriate approximation to combine homogeneous observations of this channel. The resulting homogenized data set serves as approximation of multiviewpoint SDO/AIA observations in four spectral bands (171, 193, 211, and 304 Å).

We resample all data to 12 pixel^–1. The focal length f, as defined by our ray tracing method (discussed below), is then accordingly computed individually for each observation,

$$\begin{eqnarray}&&f=\displaystyle \frac{W/2}{\arctan (1.2\cdot W/2)},\end{eqnarray} \tag{ 1 }$$

where W refers to the width of the image. Note that this focal length is computed from the satellite position and observed image and is not directly related to the actual focal length of the instrument due to the adjustment of the image size and resolution. Furthermore, we only consider square images. In this context, we analyze the propagation of diverging rays from the observer, with the parameter f specifying the angle of view. This approach departs from the assumption of small angles, where rays would be considered parallel. For each observer, we specify the time of the observation, the distance from the Sun, and the Carrington longitude and latitude.

We adapt the approach from Mildenhall et al. (2021). We use a neural network to represent the 3D scene and ray tracing methods for learning the representation. Information on the entire scene is stored in the weights of the neural network, rather than having an explicit grid representation. To be precise, we train a neural network for each individual scene, in contrast to classical deep-learning applications, which learn to perform a general task.

For each image pixel, we determine the light ray path using known information from the observer position and viewing angle. Along this path, we sample points (x, y, z) within the simulation volume and use classical ray tracing methods to compute the integrated pixel value. In the NeRF approach proposed by Mildenhall et al. (2021), each point is associated with a color and density. The pixel value is computed by the weighted sum of the colors of the sampled points along the ray, where the weights are the cumulative sum of the density values. In other words, this approach computes the likelihood of a ray being absorbed along the line of sight. The training images, together with the observer position, serve as the target for model training. By iteratively updating the neural network to match the pixel values of the training images, we obtain the underlying scene representation. This approach has shown the ability to obtain photorealistic images at interpolated positions (Mildenhall et al. 2021; Li et al. 2022).

In contrast to a classical NeRF, where the volume is composed of colors and densities, the solar atmosphere consists of emitting and absorbing plasma. In addition, NeRFs require dozens of images in order to learn the scene representation, while only three simultaneous observations are available for observations of the solar EUV corona. Therefore, we apply four primary changes to the NeRF approach. (1) We replace density and color with wavelength channel–specific emission and absorption. These capture how plasma with temperatures observed by a given wavelength channel emits and absorbs light. Each pixel value refers to the total intensity that is observed along the line of sight for a given wavelength channel and is expressed in units of data number per second (DN s⁻¹). (2) We utilize the temporal evolution and solar rotation to obtain more viewpoints (see Frazin et al. 2005). From this approach, we obtain a full scan of the Sun per observing instrument within one solar rotation (∼27 days). For the definition of viewpoints, we consider all observations in heliographic Carrington coordinates, which are solarcentric. We explicitly include time in the query coordinates of the NeRF to account for a consistent treatment of dynamic changes in the solar atmosphere. (3) We account for the known geometry of the Sun by only sampling points up to the solar surface. With this, we avoid an invalid placement of emitting plasma from below the solar surface. The change to monochromatic intensity values, instead of colors, further implies that there is no need for a background (see Mildenhall et al. 2021). (4) An optional component for NeRFs is a viewing-angle dependence of the pixel color. For our implementation, we do not consider any view dependence of the observed intensities.

In Figure 1, we provide an overview of our approach. For our 3D representation, we use a fully connected sinusoidal representation network (Sitzmann et al. 2020) with eight layers, where each contains 512 neurons, resulting in about 2M free parameters. Our model takes four input coordinates (x, y, z, t) and outputs the estimated emission (x, y, z, t) and absorption κ(x, y, z, t) coefficient for a given wavelength. We use an exponential and ReLU activation (setting negative values to 0) for the emission and absorption coefficients, respectively. With this, we obtain positive values for both coefficients. We then compute the total emission E at each point by multiplying with the spacing of the sampled points ds,

$$\begin{eqnarray}&&E=\epsilon \cdot {ds}.\end{eqnarray} \tag{ 2 }$$

Each emitted intensity value propagates through the line of sight to the observer, where absorption can reduce the intensity value at points between the observer and the point of emission. We model absorption values between [0, 1], where 1 refers to full transmission and 0 to total absorption. We compute the absorption A for each sampled point based on the spacing,

$$\begin{eqnarray}&&A=\exp (-\kappa \cdot {ds}).\end{eqnarray} \tag{ 3 }$$

The total intensity I is then computed as the sum of the absorption reduced emissions,

$$\begin{eqnarray}&&I=\displaystyle \sum _{k}^{N}{E}_{k}\displaystyle \prod _{i}^{k-1}{A}_{i},\end{eqnarray} \tag{ 4 }$$

where i and k refer to the points along the sampled ray and N refers to the total number of sampled points along the ray. Therefore, each emission E is reduced by all absorption values A in the range from the emitting point to the observer. We note that the first point is always transmitted without absorption (A₁ = 1), and the last sampled absorption point has no effect.

Zoom In Zoom Out Reset image size

For all computations, the spatial scales are normalized to 1 R_⊙. We use positional encoding, which implies periodic boundary conditions at 2π. Therefore, we scale the coordinates to units of 2 R_⊙ prior to the encoding (simulation volume of [−2π, 2π] R_⊙) to avoid reflections at the boundary.

For training, we adjust the pixel-value range based on the average maximum value (see Jarolim et al. 2024) and use an asinh stretch to mitigate the contributions by bright pixels. To account for this scaling, we use the same stretch function when computing the final pixel values (I_scaled) of our rendered images,

$$\begin{eqnarray}&&{I}_{\mathrm{scaled}}=\displaystyle \frac{\mathrm{asinh}(\widehat{I}/a)}{\mathrm{asinh}(1/a)},\end{eqnarray} \tag{ 5 }$$

with a = 0.05.

We apply two sampling strategies to prioritize points from the most important regions. (1) We only sample from −1.3 to +1.3 R_⊙ with respect to the solar center, such that the default sampling range spans 2.6 R_⊙. In addition, we stop the sampling where light rays intersect the solar surface (1 R_⊙), such that a ray that is pointing at the center of the solar disk is sampling a range of 0.3 R_⊙. Therefore, the total line of sight is dependent on the viewing direction. We note that this sampling is independent of the observer distance. (2) Similar to Mildenhall et al. (2021), we use a coarse and a fine network. For the coarse network, we sample 64 points uniformly distributed along the ray. The emission values of the coarse network are then used as a probability distribution to sample 128 additional nonuniformly distributed sampling points. These points and the 64 uniformly sampled points are then used as input in the fine network, which serves as the primary training objective. For all our evaluations, we use the outputs of the fine network.

For our model training, we select random rays from the full set of image pixels, where the number of rays per batch is adjusted to the available computational resources. For our reconstructions, we consider each wavelength separately. Therefore, our model is trained with observations from a single wavelength channel.

We initiate the model training with a centered crop of the image, where we use pixels within [−1000'', 1000''] in helioprojective coordinates. We train on this data set for one epoch (∼100,000 iterations) to prevent divergence at the beginning of the model training and emphasize the reconstruction of the solar disk in the initial phase of training. Model training is then continued with the full field of view of the instruments.

When training for smaller fields of view (subframes), we use a model trained on the full field of view as an initial starting point. This limits divergences and already incorporates the global geometry of the Sun. In the second training step, we crop the subframes at a fixed coordinate point from the full-disk observations, which serves as our new training set. Therefore, pixels outside of the region of interest are neglected and no longer suitable for further analysis after fitting the model to the subregion. The reduced spatial information that the model needs to fit enables the increase of temporal and spatial resolution (Section 3.4).

During model training, we separate one image of the training set to validate performance. We monitor the peak signal-to-noise ratio (S/N), mean squared error (MSE), and structural similarity index (SSIM; Wang et al. 2004). For all our model runs, we found a monotonic performance increase and a similar performance independent of the initialization. We note that starting with a field of view that is too large can result in irreversible divergence within the first epoch.

During each training step, we optimize both the coarse and fine models to match the target pixels. We use MSE as a loss function.

We introduce an additional regularization of the modeled absorption that reduces nonphysical results above the poles. Viewpoints from the ecliptic provide no observations of the solar surface and therefore only model low-intensity values. This allows for arbitrary absorption values at these regions. We add an additional loss factor ${{ \mathcal L }}_{\mathrm{regularization}}$ to our optimization that suppresses absorption values κ above 1.2 R_⊙ to more strongly penalize large absorption at greater heights,

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{regularization}}=\max \{d-1.2,0\}\times (1-\kappa ),\end{eqnarray} \tag{ 6 }$$

where d refers to the distance from the solar center in R_⊙. When training with full field-of-view images, we add the regularization to the MSE loss, whereas we set the regularization to 0 when training with subframes (i.e., no polar regions).

To evaluate our method, we compute the peak S/N and SSIM between our reconstruction and the ground-truth image. Here, we compare scaled images (see Equation (5)) to avoid bright pixels dominating our evaluation.

To provide an uncertainty estimation, we train five independent models with randomly initialized weights. To generate uncertainty maps, we render the same viewpoint per model and compute the standard deviation of the resulting set of images. Therefore, for regions that are less constrained by observations (e.g., single viewpoint), different solutions can be found. This degree of freedom is then reflected by our uncertainty metric.

Calculating height from imaging data is challenging, since there exists no exact spatial point that can be associated with an observed intensity value. We determine height information by computing the average distance from the solar center for all sample points along the ray using channel-specific emission as a weight. In other words, for each ray, we compute the height of the average emission value. This implies that dark regions (e.g., coronal holes, filament channels) appear as elevated structures, since the primary emission is coming from higher layers, and should not be considered for mapping the spatial location of the solar feature.

For the visualization of 2D slices of the solar atmosphere, we sample radial points along a given longitude. We pass each point through our network to obtain the respective channel-specific emission and absorption coefficients. For the emission, we use an asinh stretch, analogous to Equation (5).

As a baseline approach for rendering novel viewpoints, we use a simple reprojection. More specifically, we combine observations at a given time from STEREO-A and -B and SDO into a synchronic (Carrington) map. In the case of overlapping regions, we compute the mean value of the intersecting pixels. For missing regions (i.e., polar regions that are not covered by any instrument), we use the mean value over the entire image as a filling value. To synthesize novel viewpoints, we project the map onto a unit sphere and rotate the perspective to the new viewpoint.

This approach has several shortcomings. By projecting to a 2D map, we lose the off-disk information; consequently, this method cannot account for the faint emission that is observed over the solar limb. This further approximates the atmosphere as flat, neglecting the height information when projecting to new viewpoints. This is especially problematic for extended structures (e.g., coronal loops), where projection effects become especially prominent closer to the limb. Finally, combining overlapping regions is problematic. Due to the projection effects, the observations of the individual instruments do not align pixelwise in the synchronic map. Therefore, computing the mean of the overlapping pixels leads to blurred features.

Here, we only consider a basic approach as our baseline but note that there are more advanced data homogenization methods (e.g., Hamada et al. 2020) and image stitching methods (e.g., Caplan et al. 2016) available.

To validate our method, we use the simulated data introduced in Section 2.1. For our model training, we consider each wavelength separately and use observations viewed from ∣latitude∣ ≤ 7° from the ecliptic (i.e., 32 images). This constraint matches the observing capabilities of the SDO and STEREO satellites, which are solely located on the ecliptic plane. Our data set corresponds to a single snapshot; therefore, we consider a static atmosphere where we set the temporal coordinate to 0. For our evaluation, we focus on the comparison of ground-truth and SuNeRF reconstructed images. With this, we assess the global plasma distribution and can estimate the model performance with standard image metrics (i.e., peak S/N, SSIM).

In Table 1, we summarize the quantitative evaluation of a baseline spherical reprojection and SuNeRF reconstructions. We compute the peak S/N, SSIM, relative mean absolute error (MAE), and relative mean error (ME) based on images from nonecliptic viewpoints that are not seen by the model during training (i.e., a test set of 224 images) and the synthesized SuNeRF images. As a baseline, we use a standard reprojection of synchronic maps, where each point is mapped to a spherical surface (Section 2.4). We note that this method cannot account for the off-limb emission of the Sun, while our SuNeRF method can take the extent of the 3D corona into account. The evaluation shows that our method provides a clear improvement over the baseline approach. Moreover, the metrics averaged over all considered wavelength channels indicate a strong similarity to the ground truth, with an MAE of 0.5%, an SSIM of 0.97, and a peak S/N of 40.4 dB, which throughout outperforms the baseline with an MAE of 3.2%, an SSIM of 0.73, and a peak S/N of 22.6 dB. The ME indicates that there is no systematic over- or underestimation (maximum ME: −0.3%).

In Figure 2, we compare ground-truth and reconstructed EUV images at different latitudes. Our method results in almost identical images to the ground truth. The difference maps indicate that errors primarily originate from off-limb features and pixelwise shifts. In addition, there is only a small quality decrease with increasing latitudes (e.g., peak S/N from 46 to 41 dB for 193 Å), which suggests a valid approximation of the 3D geometry of the solar corona. This is also visible from the pixelwise error maps. Our SuNeRF model achieves the highest similarity in quiet Sun regions, as well as the interior of coronal holes. On the other hand, features such as active regions are more complex, in addition to spanning a larger range of emission values. This results in larger differences when compared to the ground truth. Larger errors are expected close to the limb as the model integrates over a larger range of uncertainties. Therefore, independent of the viewing angle, inferences on the limb will have an increased uncertainty.

Zoom In Zoom Out Reset image size

To further estimate model errors, we render images from an ensemble (Lakshminarayanan et al. 2016) of five SuNeRFs with different random weight initializations and compute the standard deviation per pixel (Figure 2). The resulting uncertainty maps show a good spatial overlap with the pixelwise difference maps. We perform a pixelwise comparison of the uncertainty estimates and the MAE over the full test set (see also Appendix A for pixelwise scatter plots), where we find a Pearson correlation of 0.63 and a Spearman correlation of 0.69. This indicates that the ensemble approach reflects the model errors and can serve as an error estimate in the absence of a ground truth.

We evaluate our method using multiviewpoint EUV observations from two different space missions, SDO and STEREO, with the latter mission consisting of twin satellites, STEREO-A and STEREO-B. All three spacecraft orbit in the ecliptic plane, with SDO in geosynchronous orbit around Earth and the STEREO spacecraft leading (STEREO-A) and trailing (STEREO-B) Earth. Even for constellations where we have full 360° coverage of the Sun, we find that three simultaneous images are insufficient for our model training. Therefore, we make use of the solar rotation that allows us to have a full scan of the Sun in ∼27 days, even for a single instrument. For our model training, we consider a sequence of 14 days at 1 hr cadence (from 2012 August 24 00:00 to 2012 September 7 00:00 UT; a total of 681 observations). We explicitly encode time as a coordinate for the sampling along light rays (Figure 1), which allows SuNeRFs to model temporal changes. We randomly sample 32,768 rays within the volume that we then use in our training set. Our model is trained for ∼250,000 iterations (three epochs) until the model converges (about 2 days on 8 × A100 GPUs). For our analysis, we primarily focus on the 193 Å channel, but we also provide reconstructions of the 171, 211, and 304 Å channels. Before training, we adjust instrumental differences, including differences in the wavelength bands, using ITI (Jarolim et al. 2024), allowing us to generate homogeneous observations of SDO and STEREO.

In Figure 3, we show snapshots from a global reconstruction of the solar atmosphere. The arrows in Figure 3(a) indicate the viewpoints used for the reconstruction, with latitudes ranging from −4° to 8°. The central image shows the reconstruction of the solar south pole at 2012 August 30 00:00:00 (UT). This directly enables us to further investigate the polar coronal hole, which we segment using a deep-learning tool for automatic coronal hole detection (CHRONNOS; Jarolim et al. 2021). We then compare two approaches in identifying the coronal hole boundaries. (1) We use the individual satellite observations (i.e., full-disk images) to identify coronal holes and then reproject the combined binary maps to the polar viewpoint (blue contours). (2) We use SuNeRF to reconstruct the polar observation and then apply the detection method (red contours). The comparison shows that the strong projection effects only allow for a partial identification of the coronal hole boundary from the ecliptic plane. Coronal holes are also frequently obscured by loop systems in adjacent active regions, which makes detection and consequent space-weather predictions challenging. Our method can account for these projection effects by mapping the 3D geometry (see the online animation of Figure 3). From our comparison to simulated coronal holes, the SuNeRF reconstruction provides a more consistent reference to estimate the coronal hole boundary, specifically in the case of strong reprojection effects, where errors are expected to be in the range of the uncertainty maps (Figure 3(b)).

In panel (b), we compare viewpoints from different latitudes against the corresponding baseline reprojections. SuNeRF clearly recovers the polar coronal hole by combining the multi-instrument data. Artifacts occur only for the off-limb region that is not sampled by the ray tracing (i.e., bottom left at −90°).

In the bottom row of Figure 3(b), we provide uncertainty estimates derived via an ensemble of five individual runs performed with random initializations. The uncertainty maps serve as a proxy for the model error in the absence of a ground truth (Lakshminarayanan et al. 2016). The uncertainty is generally higher for observations than for the model data in Section 3.1. This uncertainty likely arises from the increased complexity of the observations (i.e., the high-frequency features are more difficult to match) compared to the smoother atmosphere in the simulations. Note the higher number of small regions (e.g., bright points) with errors >10%. These could be associated with a misalignment in temporal evolution, where the 1 hr temporal resolution may be too coarse to capture small brightness variations. Overall, large-scale features (i.e., active regions, coronal holes, filaments) are consistently reconstructed by our method, as can be gauged by the uncertainty maps, and the main errors are primarily related to pixelwise shifts (e.g., the coronal hole boundary).

Using SuNeRFs, we can now determine the height of the solar features from the 3D reconstructions. In Figure 4(a), we extract a meridional slice to analyze the emission profile from 1.0 to 1.3 R_⊙. From north to south, we observe the emission profile of an active region, a coronal hole (blue), a filament channel (orange), and the polar coronal hole. The 2D coronal slice in Figure 4(b) shows regions of lower emission close to the surface. A possible explanation for the modeled emission profile is that our model accounts for the solar chromosphere and transition region, which are not visible in EUV due to temperatures being below the million K range, to generate EUV emission. In contrast, bright active region loops are rooted in the solar photosphere. The coronal hole can be clearly identified from the decreased emission of the entire column. The solar filament channel can be identified from the dark cavity that appears as a distinct and sharp drop in intensity. Lower layers show a similar brightness as the quiet Sun corona, where the reduced emission primarily originates from the dark filament channel. The most distinct emission decrease occurs at the solar south pole, where we observe the polar coronal hole (Section 3.2). From the emission, we can estimate the relative plasma density by assuming that the emission is proportional to the squared density. This provides a direct estimate of the coronal hole boundaries, which appear as distinct jumps in the integrated emission profiles for the polar and low-latitude coronal holes (Figure 4(c)).

Zoom In Zoom Out Reset image size

To investigate solar transient events, we increase the temporal cadence of observations in our training set and focus on a subframe instead of the full field of view. Here, we apply our method to the solar eruption that started at 2012 August 31 19:30 (UT), and we use observations from STEREO-B and SDO in the 304 Å channel. We use the parameters learned from the 304 Å model in Section 3.2 as a starting point, which already provides the general geometry and structure of the considered time period. We focus the training on 1024 × 1024 pixel subframes centered at the erupting solar filament (Carrington longitude of 90° and latitude −20°) at a 1 minute cadence and train our model for an additional 100,000 iterations (∼four epochs).

In Figure 5(a), we give an overview of the filament eruption event and show SuNeRF-rendered images during and after the eruption (top and bottom row, respectively). We provide height maps that give the distance from the solar center for the observed emission (Section 3.3). The reconstruction shows that the prominence rises slowly (yellow arrow), followed by a rapid ejection of the filament plasma into interplanetary space (coronal mass ejection). Afterward, the hot flare ribbons, which are clearly mapped close to the solar surface, are formed, in agreement with our physical understanding of eruptions (blue arrow). In addition, we show the total absorption along the line of sight, where the filament structure can be clearly identified. The ejection can be related to the depletion of absorbing plasma.

Triangulation is a useful approach for height estimations of dense plasma but can lead to ambiguous results for an optically thin medium (Aschwanden 2005, 2011). With the full mapping of the solar atmosphere, we can better analyze the underlying 3D structure and study emission and absorption profiles along their radial extent. In Figure 5(b), we show the time evolution of an extracted slice of the modeled solar atmosphere (blue line in Figure 5(a)). The top row shows the emission profile. After the eruption (∼19:40), postflare loops and their related flare ribbon footpoints are formed (blue arrow). The loop height can be directly estimated from the profiles (∼0.7 R_⊙ or 500 Mm above the solar surface). The absorption profiles show the evolution of the filament plasma and allow one to determine the direction and velocity of the erupting structure.