Tomographic iterative reconstruction of a passive scalar in a 3D turbulent flow

Turbulence in the stable planetary boundary layers often encountered in high latitudes influences the exchange fluxes of heat, momentum, water vapor and greenhouse gases between the Earth’s surface and the atmosphere. In climate and meteorological models, such effects of turbulence need to be parameterized, ultimately based on experimental data and/or high-resolution models (Figure 1). A novel experimental approach is being developed within the COMTESSA project in order to study turbulence statistics with high resolution. Using controlled tracer releases, high-resolution camera images and estimates of the background radiation, different tomographic algorithms can be applied in order to obtain time series of 3D representations of the scalar dispersion.

Figure 1: A turbulent flow at the forest-atmosphere interface as simulated by a meter scale Large Eddy Simulation and colored by the temperature.

The synthetic data was generated from large eddy simulations with the PALM model (PArallelized Large-Eddy Simulation Model for Atmospheric and Oceanic Flows) and used for the tomography observing system simulation experiments. The maximum available resolution is currently 1025x258x325. Figure 2a shows the sample output of the model in low resolution used for initial testing.

 

Figure 2: a) Sample output of the model in low resolution used for initial testing. b)  sample reconstruction with 50 iterations of the Kaczmartz method and 20 cameras symmetrically distributed around the vertical axis.

Tomography with iterative algorithms: Tomography was introduced by Radon in 1917. Classical medical tomography based on Radon transform based methods is not applicable to the current setting with sparse measurements. We focus on different iterative solvers instead, both sequential (Kaczmartz aka. ART and variants) and simultaneous (Landweber, Cimmino, SART, SIRT). We also used CGLS for comparison only in an idealized setting. We use a custom version of the AIRTtools and the ASTRA packages. Figure 2b shows a sample reconstruction at the same resolution of the phantom with 50 iterations of the Kaczmartz method and 20 cameras symmetrically distributed around the vertical axis.

Cloud computing with GPUs: We performed most simulations in a dedicated instance of the Amazon Web Services cloud.  The instance consists of 2500 GPU cores providing substantial acceleration with respect to the standard CPU compiled algorithms.

Preliminary results

1. Semi convergence and stopping rules: The quality of the reconstruction of Ax=b depends on the number of iterations, and for noisy data the over-fitting error can blow up. Stopping rules terminate the iterations near “semi-convergence”. Stopping rules include

• Fit to Noise Level (FN): smallest k for which

• Unbiased predictive risk estimation (UPRE): smallest k that minimizes

• Generalized cross validation (GCV): smallest k that minimizes

(ρ is the residual, η the noise level, m is the dimension of b). They require estimates of t_k , the trace of an analogue of AA^T. Figure 3a shows the stopping functions and 3b illustrates the evolution of the error in an OSSE for three stopping rules and two estimates of the trace t_k.

 

Figure 3: Above: Stopping functions for FN (left), UPRE (center) and GCV (right) with two trace estimates. Below: Evolution of the error in an OSSE for three stopping rules and two estimates of the trace t_k.

2. Discretisation resolution: The “inverse crime” occurs when the same ingredients are employed to synthesize as well as to invert data. In order to prevent trivial results while calibrating the methodology, different resolutions for the phantom and the sinogram can be used. Figure 4 illustrates the effect using a slice of a 3D synthetic flow data as a phantom.

Figure 4: Assessing the dependence on reconstruction resolution. 1^st row: phantom. 2^nd row: sinogram. 3^rd row: reconstruction.  1^st column: high-resolution phantom, high-resolution reconstruction. 2^nd column: high-resolution phantom, low-resolution reconstruction. 3^rd column: low-resolution phantom, low-resolution reconstruction.

3. Geometry of the camera positions: In contrast with medical tomography where many images are obtained by rotating an object, remote sensing with individual cameras yields significant gaps between views and the camera locations have a large influence on the quality of the reconstruction. Using simulation experiments, the theoretical optimal distribution can be estimated for a given set of restrictions.  Figure 5 a-c illustrates 6 steps in the approximation process using simulated annealing for an ideal 2D case.

Figure 5: Assessment of the optimal camera distribution using simulated annealing for an ideal 2D case showing successive approximation steps.

 

Outlook: We preliminary described some settings that can affect the quality of the tomographic reconstruction for a synthetic tracer plume modeled in large eddy simulations. When the target resolution of the LES model and reconstructions exceeds the available observation, dimensionality reduction methods are necessary. The first measurement campaign took place this summer yielding high frequency time series of UV camera images.