Tomographic imaging of non-local media based on space-fractional diffusion models Free

Tomography, or tomographic imaging, is a particular class of inverse problems that enables reconstructing an image of the internal structure of a body (or of specific material properties) without requiring intrusive access to it. Depending on the specific method, different types of penetrating waves or field transport phenomena can be selected to probe the body. Typical examples include electrostatic and electromagnetic waves (employed in Electrical Impedance Tomography-EIT and Magnetic Resonance Imaging-MRI, respectively),¹ acoustic waves (acoustic tomography),² and thermal fields (thermography).³ In the past few decades, tomographic techniques have found many engineering applications including biomedics, atmospheric science, geology, oceanography, plasma physics, materials science, astrophysics, and acoustics.^4–6 More recently, tomography was also successfully applied to structural health monitoring and non-destructive testing for the inspection and characterization of structures and materials.^3,7–12

As many inverse problems, tomographic methods are heavily dependent on the accuracy of the underlying mathematical model used to predict the response of the medium. Typically, inverse methodologies fit, in a least-squares sense, a mathematical model to experimentally measured data by iteratively tuning the model parameters.^1,13,14 It is evident that the ability to reconstruct accurate images and to improve the resolution is strongly dependent on the ability of the model to capture and properly simulate the underlying physical response of the system. Typical applications of tomographic methods deal with domains of analysis that are highly inhomogeneous (e.g., including geometric or material inhomogeneities), have complex material properties (e.g., viscoelastic materials), and are eventually characterized by hybrid and non-local transport behavior. If not properly accounted for, these properties limit both the sensitivity and resolution of the methodology and can even prevent entirely the possibility to reconstruct a meaningful image.

To date, several numerical and experimental studies have demonstrated that field transport processes in complex inhomogeneous media (e.g., random and/or porous materials) can occur according to hybrid modalities^15–20 and anomalous behavior.^21–24 The hybrid transport typically involves a combination of propagating and diffusive mechanisms.^25–29 The diffusion processes can be either standard that is characterized by a typical Gaussian distribution of the field quantities or anomalous where instead an α-stable distribution emerges. In general, diffusive processes are characterized based on the functional relationship that describes the evolution of the variance of the particle displacement with time. This relationship is expressed as $⟨x⟩2∼tγ$ where the normal diffusion is characterized by a linear scaling γ = 1, while the anomalous diffusion occurs for $γ≠1$⁠. In particular, processes described by γ > 1 correspond to diffusion phenomena that evolve in time faster than the standard (Gaussian) process and are therefore denoted as superdiffusive. On the contrary, processes that are characterized by 0 < γ < 1 are denoted as subdiffusive, because the variance increases more slowly than what predicted by the classical (Gaussian) diffusion model. Subdiffusive processes can be described at a microscopic level by continuous time random walk (CTRW) models employing a probability density function for the waiting time of a particle jump characterized by heavy tailed distribution.^30,31 On the other hand, superdiffusive phenomena can be described by Lévy flights or Lévy walks where the probability density function of the step length is a heavy-tailed distribution (that emerges as a consequence of the long-range interactions among particles).

These unconventional transport mechanisms are typically driven by non-local effects, such as long-range interactions^32,33 and memory and hereditary mechanisms,^34,35 and have been shown to be well described by fractional order models. At macroscopic scales, the CTRW subdiffusive processes can be mapped to a time-fractional diffusion equation, while the superdiffusive processes map into a space-fractional diffusion equation. It is worth emphasizing that the first type of process is linked to a power-law distribution of the waiting time between successive particle jumps, while the second is due to power-law distributions of the jump length.

Examples of these unusual transport behaviors can be observed in different media and fields of application including, but not limited to, porous soils,^36–40 heterogeneous aquifers,^21,41 and underground fluid flow.^24,42

Meerschaert et al.^42,43 performed an extensive research campaign focused on the transport of contaminants in heterogeneous soils and porous materials and showed that the underlying advective-diffusive behavior exhibits indeed a fractional nature. Other studies^{28,31,44–47} have highlighted the occurrence of regular and anomalous diffusion mechanisms in the propagation of electromagnetic and acoustic waves in random and periodic materials. Propagation of light through fog⁴⁸ or murky media,⁴⁹ of sound through a forest or a urban environment,^50–52 and of seismic waves through the ground^53,54 are all examples of practical applications where hybrid and non-local transport mechanisms can take place. In all these situations, fractional calculus was shown to be a powerful tool to properly capture and simulate complex transport mechanisms in continuum media.

From the discussion above, it appears that the integration of fractional order models within tomographic methods could improve the imaging performance and even allow capturing details and inhomogeneities that would otherwise go undetected when using conventional imaging approaches based on integer order models.

In this paper, we present a numerical investigation into the application of fractional order continuum models to tomographic imaging techniques. A few applications of fractional calculus to inverse problems have already been explored in the literature. Cheng et al.⁵⁵ successfully solved a one-dimensional fractional diffusion inverse problem to determine the order of the temporal fractional derivative and the spatially varying diffusion coefficients. Miller and Yamamoto⁵⁶ formulated a parameter identification problem based on the fractional diffusion equation. Battaglia et al.⁵⁷ proposed the solution of an inverse heat conduction problem for parameter estimation based on non-integer forward models. The precision of the inverse procedure proposed in their method heavily depended on the accuracy of the identified model that had to be determined by experiments. Murio^58,59 presented the numerical solution of an ill-posed problem consisting in determining analytical functions of the boundary temperatures and the heat fluxes based on transient temperature measurements at some internal point of a 1D conductor. Jin and Rundell⁶⁰ investigated the degree of ill-posedness of a series of theoretical and numerical inverse problems based on fractional differential equations involving the Caputo definition in both time and space. They found out that the fractional character of the operator can either improve or worsen the conditioning of the inverse problem based on the type of input data and quantities to be reconstructed. Kirane et al.⁶¹ solved a 2D source reconstruction problem for a time fractional diffusion equation using biorthogonal sets of functions.

In this study, we explore the capabilities and performance of fractional tomography applied to two-dimensional domains. More specifically, we refer to fractional tomography as a tomographic approach relying on fractional-order continuum models for the solution of the forward problem. The use of these models within a tomographic imaging framework is expected to impact the reconstruction performance at multiple levels:

1. It was previously discussed that anomalous diffusion processes, such as those arising in imaging of highly aberrating and scattering media, are associated with non-Gaussian distributions of the field quantities having noticeable heavy-tails. Note that a distribution of a random variable is said to be heavy-tailed if the moment generating function of the distribution function always diverges.⁶² Typically, the content of the tails is discarded because it is not considered as a primary source of information for the imaging process. However, while this could be a reasonable assumption when the domain is dominated by classical diffusion (in this case the tails of the Gaussian distribution have a low content of information), in the presence of anomalous diffusion the tails contain a non-negligible amount of information about the interior structure of the domain to be imaged. Hence, neglecting the heavy tails results in considerable and irremediable loss of information.

2. Recognizing and exploiting the fractional nature of the host medium and the information contained in the heavy tails can represent a turning point in the development of accurate imaging technologies capable of sensing through highly inhomogeneous and scattering media.

3. In order to extract meaningful information from the non-Gaussian distribution of the field quantities, the underlying mathematical model should be able to simulate these hybrid and anomalous physical mechanisms and to relate them to the physical parameters of the medium.

In order to clearly illustrate the methodology and without limiting the generality of our approach, we focus the following study on diffusive transport processes. An example of a classical tomographic method that falls under this category is Thermal Tomography (TT),^3,8,63,64 where diffusive heat transfer is used as a probing mechanism to sense the medium. Note that the results presented below could be easily generalized to other approaches whose transport mechanism can be effectively described by fractional models.

Tomographic imaging techniques are usually formulated as iterative inverse problems in which measured data are fit, in a least-squares sense, to a mathematical model that simulates the response of the host medium.¹ Although many inversion methodologies have been presented over the years, the most common approach to solve complex geometries is based on iterative numerical techniques. These methods typically consist of two main steps: the forward problem and the inverse problem. The forward problem simulates the response of the system assuming that the excitation conditions and the system parameters (e.g., material properties, etc.) are known. The inverse problem relies on an optimization approach designed to minimize the least-squares error between the measured and the predicted response at selected measurement points. The system parameters, or a subset of them, are typically selected as design variables.

In this work, we use fractional order models in the forward problem in order to describe the heat transport process at the basis of the thermographic technique. We will show that the introduction of fractional models has two main effects on the general formulation of the tomographic problem: (1) if the behavior of the real physical system is fractional, the model will capture information that would otherwise go undetected when using an integer order model; however, (2) if the behavior is not fractional but the system is simply inhomogeneous, the fractional formulation will provide an enlarged parameter space (which now includes the order of the differential operator) in which the reconstruction can be performed.

From a general perspective, thermal tomography reconstructs the internal (thermal) properties of a solid medium based on the measurements of boundary (or surface) temperature fields. The data acquisition procedure is performed by heating the sample at a specific location and measuring the resulting surface temperature at multiple locations. The same procedure is repeated for several source locations in order to improve the conditioning of the inverse problem and increase the accuracy of the reconstruction. In the numerical analyses performed in this study, the stationary heat sources were placed near the boundaries as shown in Fig. 1. The response of the system was estimated in terms of absolute temperature measurements performed at the sensor locations (Fig. 1). The thermographic technique is typically formulated as an optimization process based on iterative solution of the forward and inverse problem. The thermal properties, such as the conductivity or the heat capacity, are typically used as design variables.

In the forward problem, the thermal properties are assumed known (either due to an initial guess or to values available from the previous iteration) so that the mathematical model allows estimating the thermal field in the entire domain. In the inverse problem, the measured temperature values at selected locations are compared with numerical predictions in order to find an updated distribution of thermal properties that minimizes the least-squares error. The reconstruction of different thermal properties is an indicator either of material inhomogeneities or of possible defects, depending on the specific application of the tomographic technique.

The fractional heat diffusion process in a 2D domain is described by the following fractional partial differential equation:⁶⁵ 

where T(x, y, t) is the temperature field at time t in the finite domain a < x < b, c < y < d, the terms $c1±(x,y)$⁠, $c2±(x,y)$ are spatially dependent functions that conceptually play the same role as thermal diffusivity in integer-order models. These coefficients have dimension $Lαt−γ$⁠.⁴⁰ In the following, we will refer to them as thermal diffusivities but with the understanding that this terminology should be intended in the context of the fractional equation. Note also that these coefficients could potentially be functions of the temperature, in an analogous way to the dependence observed in the thermal conductivity and heat capacity of the integer-order model. However, such dependence would make the equation nonlinear and therefore it would require a different type of approach.

The sign subscript is associated with the left- and right-handed fractional operators. In the following simulations, we will assume that $c1+(x,y)=c1−(x,y)=c2+(x,y)=c2−(x,y)=c(x,y)$ which means that we consider a symmetric non-local diffusion process. The notation $∂±$ in Eq. (1) refers to the left- and right-handed Riemann-Liouville fractional derivative defined as

where Γ is the Euler gamma function and n is an integer such that $n−1≤α≤n$⁠.

Equation (1) represents a general form of the fractional wave-diffusion equation capable of representing several forms of transport mechanisms.⁶⁶ When $0<γ≤2$ and $α=β=2$⁠, the equation describes a time-fractional process having diffusion-like characteristics when $0<γ≤1$ and wave-like characteristics when $1<γ≤2$⁠. When γ = 1 and $α=β=2$⁠, we recover the classical diffusion equation. For γ   1 and $0<[α,β]≤2$⁠, Eq. (1) describes a space-fractional diffusion process whose solutions belong to the superdiffusion regime.⁶⁷ The superdiffusive behavior is the results of dynamics dominated by long-range interactions (the so-called Lévy flights), where the step length distribution decays asymptotically as $x−(α+1)$⁠. This is also the interval of interest for the present study in which we want to explore the regime characterized by heavy-tailed distributions of the field quantities due to Lévy flight dominated dynamics. It is also worth highlighting that this macroscopic diffusion model is very general in nature and applicable to heat transfer problems independently of the microscopic mechanism dominating the transport process itself (e.g., either phonon- or electron-driven transport).

In the following numerical analysis, we will consider the case of γ = 1 and α = β. Once complemented with proper boundary conditions, this space-fractional diffusion equation (1) can be solved numerically.

The thermographic problem considered below is solved using the following initial conditions:

and boundary conditions

where Ω and $∂Ω$ indicate the 2D domain and its boundary [Fig. 1(b)], respectively. We assume a square domain having a side length of l = 1 m.

The initial value problem defined by the set of equations (1, 3, 4) can be solved numerically using a shifted Grunwald-Letnikov (GL) finite difference (FD) scheme.⁶⁷ More specifically, the two-sided shifted Grunwald-Letnikov scheme⁶⁵ was selected due to its ability to respect the local symmetry of the solution at every node of the domain. Figure 1(b) shows a uniform discretization grid used for the numerical solution of equation Eq. (1). The grid has a uniform spacing with $Δx=(b−a)/I$ and $Δy=(d−c)/J$⁠, where I and J are the total number of discretization nodes in the x and y directions, respectively. The finite difference approximation to $T(xi,yj,tn)$ will be denoted as $Tijn$ where $xi=iΔx, yj=jΔy$ and $tn=nΔt$⁠.

Following the above definition, the general forms of the shifted Grunwald-Letnikov formulas are⁶⁷ 

where $M+, M−$ are positive integers, $h+=(x−a)/M+, h−=(b−x)/M−$ and $k=1,2,3…M(±)$⁠. The coefficients g_k are the Grunwald weights expressed as

Using the shifted GL relations (5) to discretize Eq. (1), we obtain⁶⁵ 

which can be explicitly solved for $Tijn+1$⁠.

The stability condition of this scheme for the case $1<α<2$ is given by⁶⁵ 

The solution of the inverse problem is performed according to a classical iterative minimization approach following the Levenberg-Marquardt (LM) method. According to the LM method

where P is the vector of design variables consisting in the system parameters (e.g., the coefficients c or the fractional order α), T is the vector of numerically estimated temperatures, Y is the vector of measured temperatures, μ^k is a regularization parameter,¹ $Ωk$ is a regularization matrix chosen in our simulations as the identity matrix, J is the Jacobian matrix, and the superscript $()T$ indicates the transpose. The term $μkΩk$ is a regularization term used to damp the instabilities introduced by the ill-conditioned Hessian matrix $JTJ$⁠.¹³ The generic element of the Jacobian, or sensitivity matrix, is defined as

with $l=1,2,…,L$ and $m=1,2,…,M$⁠, where M is the total number of unknown parameters, L is the total number of measurements (⁠$L≥M$⁠), T_l is the lth estimated temperature, and P_m is the mth unknown parameter. The iterative problem described by Eq. (9) is repeated until a prescribed error threshold is satisfied. In our simulations, the exit conditions were chosen to be a threshold value of the residual error $|Y−T(Pk)|<10−5$⁠.

In this section, we apply the methodology discussed above to reconstruct the spatial distribution of the thermal diffusivity of a 2D domain. The main difference, compared to traditional thermography, is that we let the domain have a non-homogeneous distribution of the order of the underlying differential operator. In other terms, we consider media in which the transport process is non-homogeneous and can result from a combination of standard and anomalous diffusion. The following three scenarios will be considered: (1) the reference domain exhibits a fractional-order behavior with localized integer-order sub-domains, (2) the reference domain exhibits integer-order behavior with localized fractional-order sub-domains, and (3) the reference domain exhibits integer-order behavior with local sub-domains having either different fractional order or diffusivity coefficients with respect to the background. The selection of these cases of studies has been made to illustrate some specific characteristics of the fractional operator within the context of an inverse problem.

This case addresses a situation in which the thermal field transport within the physical domain is described by a space-fractional diffusion equation. The domain includes two localized inhomogeneities consisting of areas where the transport is of integer order. The reference map showing the distribution of the order of the spatial differential operator is shown in Fig. 2(a). Note that this condition can be mapped to a physical situation where the domain of interest is characterized by long-range interactions (i.e., an anomalous diffusion process) other than for small localized areas that are controlled instead by short-range interaction (hence resulting in a standard diffusion process).

In order to illustrate the need for fractional order models in inverse problems involving anomalous transport behavior, we compared the results of the reconstruction obtained using both a fractional- and an integer-order forward model. In the case of the integer order model, the only unknown parameter that can be reconstructed is the thermal diffusivity. The reconstructed spatial distribution of the thermal diffusivity c(x, y) is shown in Fig. 2(b). A direct visual inspection of the results indicates that the integer-order model is completely unable to capture the inhomogeneities in the domain due to the inability to simulate the underlying anomalous transport process. Conversely when a fractional-order model is used, the process successfully identifies (both in amplitude and location) the inhomogeneities present in the domain [Fig. 2(c)]. In this case, the order of the fractional operator was selected as the unknown parameter describing the characteristics of the domain.

The results above are likely the most significant among the three cases and clearly highlight the necessity of using a fractional-order model when the imaging is sought in domains characterized by anomalous field transport. They also highlight another important feature of the fractional inverse problem. The fractional order can be used as an additional parameter to track the presence, the magnitude, and the location of the inhomogeneities. In other terms, the fractional order expands the parameter space that can be used to characterize inhomogeneous domains (which could be relevant in applications such as non-intrusive sensing and non-destructive evaluation).

The numerical experiments presented in the next two cases are intended to address and characterize the role of the fractional order as an additional indicator of inhomogeneities that is beyond the traditional parameters such as the transport coefficients c(x, y).

This case explores the ability to image a domain characterized by local fractional-order inhomogeneities distributed in an otherwise integer-order background. Such a situation physically corresponds to a domain in which the heat transfer in the background occurs according to a standard (Fourier) diffusion process, but it experiences localized areas in which the transport is non-local.

The reference case to be reconstructed is shown in Fig. 3(a) where the background of the 2D domain has order α = 2 other than for a small area in the center where $α=1.8$⁠. The reconstruction was performed by using a fractional-order model where the thermal diffusivity coefficients and the fractional order were both treated as unknown variables.

To solve the inverse problem, the iterative approach was initialized using the true values of both c and α in the background. The inverse procedure exited with a residual error of $5.44$ × 10⁻¹⁰ well below the set threshold. Figures 3(b) and 3(c) show the results of the reconstructed maps. In general, both parameter distributions show clear traces of the location of the inhomogeneity. However, a marked difference is visible in their corresponding amplitude. More specifically, the order of the differential operator is well captured by returning a value of $α=1.835$ which is within 2% error with respect to the reference [Fig. 4(b)].

This was expected because the original inhomogeneity was set in terms of an anomaly on the order of the operator. The thermal diffusivity distribution [Fig. 3(b)] instead shows a variation of 16% [Fig. 4(a)] with respect to the background in correspondence to the inhomogeneity despite the original configuration being uniform in terms of c(x, y). This result was not unexpected because previous studies⁶⁸ had shown that there is a direct correspondence between spatially inhomogeneous coefficients of an integer order differential equation and the corresponding order of a matched fractional equation. We note that, although the location and extent of the anomaly is still well captured, the absolute value of the reconstructed coefficient is not necessarily meaningful and does not allow quantifying of the inhomogeneity.

In this scenario, we explore the possibility to simultaneously reconstruct two distinct types of inhomogeneities affecting either the fractional order or the thermal diffusivity. This case is intended to underline the fundamental difference in the role played by both the thermo-physical properties (captured in the coefficients of the differential equation) and by the physical transport mechanisms (captured by the fractional operator). In fact, variations of the thermal diffusivity do not affect the functional form of the solution of the heat conduction equation, while perturbations on the order of the operator do affect the transport process by introducing non-local effects. The underlying idea is that a change in the nonlocal parameter α indicates the occurrence of non-local heat transfer which ultimately is connected to changes in the properties of the host medium.

Figure 5 summarizes the results of the numerical reconstruction. The background values of c and α are the same as the previous case, and the two inhomogeneities are represented by an anomaly in the thermal diffusivity of intensity c = 0.02 $[Lα/t]$ [Fig. 5(b)] and an anomaly in the fractional order distribution equal to $α=1.2$ [Fig. 5(a)]. The iterative inverse procedure was initialized using as initial guess the true values of the background distributions of c and α. The inverse procedure exited with a residual error of $1$ × 10⁻⁵.

The reconstructed distributions could correctly identify the location and intensity of both the inhomogeneities in terms of fractional order [Fig. 5(c)] and thermal diffusivity [Fig. 5(d)]. A marked inhomogeneity emerged in the reconstructed map of the diffusion coefficients. This result should have been expected given that, as shown in case 2, there is a strict correspondence between the spatially inhomogeneous coefficients of an integer order differential equation and the corresponding order of a matched fractional equation. This link between the system parameters also tends to reduce the accuracy of the identification of the fractional order, as seen in Fig. 5(c). The indeterminacy that occurs when both parameters are perturbed cannot be resolved mathematically unless more information on the host system was provided. As an example, in cases where a priori information on the inhomogeneity was available, then the inverse problem could be solved using constraints on the optimization problem so to facilitate resolving indeterminacy. Such situations might correspond to cases in which, as an example, knowledge of the materials involved or of possible structural defects was available from other measurements. Figures 6(a) and 6(b) show the percentage error maps for the fractional order and the thermal diffusivity.

The results of the numerical experiments suggest that the augmented parameter space resulting from the use of fractional order models has the potential to improve significantly the accuracy of the reconstruction procedure when compared to the conventional tomographic techniques.

In fact, the nonlocal parameter introduced by the fractional formulation (i.e., the order of the operator) serves as a new design variable able to capture the underlying nature of the heat conduction process as well as of describing the anomalous behavior observed in complex and inhomogeneous media. In addition, fractional tomography opens the way to exploit hybrid field transport (i.e., a combination of propagating and diffusive mechanisms) as a probing mechanism, whereas classical tomography is constrained to the use of a single transport process.

Overall, the above results suggest that the enrichment of the parameter space due to the integration of fractional order models in tomographic methods may have important effects on the sensitivity of the reconstruction process with respect to classical tomography.

This study presented the general theoretical and mathematical framework for the implementation of fractional tomographic imaging. This general class of techniques was conceived to achieve non-intrusive and non-destructive imaging in a variety of applications where the host medium is characterized by anomalous and hybrid field transport. Among these applications, we can include imaging of viscous, nonlocal, and highly scattering media.

The proposed fractional tomographic method was tested via numerical experiments. In this study, the fractional nature of the model was limited to the spatial derivative in order to target specifically nonlocal and anomalous diffusion processes. Numerical results showed that the order of the operator can be used as a sensitive and effective parameter to probe and image the interior of a medium. This approach provides, not only a tomographic framework with an augmented parameter space, but also a powerful method for sensing and imaging in nonlocal complex media that traditionally pose extremely challenging conditions. Even though further research is required to establish constitutive relations able to connect the inhomogeneity of the system parameters to the fractional order, the present study shows clear evidence of the potential of the fractional tomographic framework and its ability to considerably enhance sensitivity and resolution in remote sensing applications.

The authors gratefully acknowledge the financial support from the National Science Foundation under the Grant CMMI CAREER $#1621909$⁠.