NON-INVASIVE INVERSE PROBLEM IN CIVIL ENGINEERING

In this contribution we focus on recovery of spatial distribution of material parameters utilizing only non-invasive boundary measurements. Such methods has gained its importance as imaging techniques in medicine, geophysics or archaeology. We apply similar principles for non-stationary heat transfer in civil engineering. In oppose to standard technique which rely on external loading devices, we assume the natural fluctuation of temperature throughout day and night can provide sufficient information to recover the underlying material parameters. The inverse problem was solved by a modified regularised Gauss-Newton iterative scheme and the underlying forward problem is solved with a finite element space-time discretisation. We show a successful reconstruction of material parameters on a synthetic example with real measurements. The virtual experiment also reveals the insensitivity to practical precision of sensor measurements.


Introduction
Non-invasive methods are gaining increasing interest in various disciplines spanning from medicine to geophysics.The main advantage of such methods is that they preserve the observed sample intact while getting information about its inner properties.Such tasks can therefore play fundamental role in analysis of existing structures, improving future design, minimising the gap between model prediction and real observations etc. Moreover today trend is in favour of complex models described by multiple parameters which have to be properly identified from real measurements, often requiring a specially designed set destructive tests with a high precision laboratory equipment.Such process often requires a set of samples which do not have to be available and features a rather limited description of the studied material.In this contribution we propose a non-invasive parameter identification technique based merely on boundary observations of some structure.This idea was inspired by a medical imaging technique named Electric Impedance Tomography (EIT).The mathematical procedure was first rigorously described by an Argentinian mathematician Alberto Calderón in his foundational paper [1] in 1980, although the idea was studied earlier in 1930 in geophysics [2,3].The goal of EIT is to recover an electric conductivity field inside some domain using only boundary measurements.Basic idea of this method lies in the difference of surface measurements due to variations in the subsurface conductivity distribution.Meanwhile a single set of surface measurements for a given loading conditions might result in a number of possible conductivity fields, Calderón surpassed this problem by sequentially implying multiple loading conditions each followed by a measurement, which in essence has the potential to contain the complete information about the underlying conductivity distribution.
In this contribution we will employ similar techniques of Calderón inverse problem paradigm to a parabolic partial differential equation with two independent parameters to describe thermal properties of some structure.Our method is designed to be external loading-free, i.e. it completely relies on the environmental factors as a source of changes which are necessary in order to identify the material inside a domain.

Inverse problem
Parameters recovery of a given model is often regarded to as ill-posed problem [4,5].For such task it is typical that even a small perturbation in the input data can result into a large variation in the solution.This can be due to several reasons e.g.uniqueness issues stemming from information shortage, i.e. data insufficiency or it can be an intrinsic property of the system itself.For this reason needs a stabilization provided by additional assumptions, e.g.some prior information, enforcing smoothness, preferring solution with the smallest norm, provide bounds to the unknown entity, etc. Procedures for determining the constrains for solving such problems are generally called regularisation methods and within this contribution we shall employ a deterministic approach.An inverse problem is usually an optimization task solved by minimising some functional which is often stated as follows, with the first term representing the cost function and the second term being the regularisation penalty where F (θ) is output of the forward operator, u m is the vector of measured quantity, p indicates p -norm, G(θ) is the regularisation functional introducing the additional constrains to the solution and θ are the parameter fields being reconstructed.
A typical choice of the regularisation penalty term in (1) is of following form [4] where the hyper-parameter, λ controls the trade-off between solution stability, given constrains and a distance from the actual solution.Additional parameters θ r represent some a priori known and possibly non-smooth behaviour of the unknown parameters [6].
A straightforward way of explaining the regularisation operator L is such that it draws the solution towards its null space, i.e. ker(L).In our case it takes a form of a discrete approximation to the Laplacian.In this contribution, both terms in (1) were chosen to be 2 -norm due to the convenience for computational purposes and preference of smooth solutions [7].In such case, by combining the minimisation scheme (1) for p = 2, penalty term (2) and a Gauss-Newton approximation of Newton-Raphson multi-variable method, one can obtain a following iterative formula [6] ) where J k is the Jacobian evaluated at θ k , δu is a vector containing the difference between the model output and real measurements, θ r is a conductivity field containing its prior assumptions and R = L T L is regularisation operator composed from a pre-calculated Laplacian L of a piece-wise constant functions on a finite element mesh.
where F (θ) ∈ R vw represents a discrete Neumann-to-Dirichlet operator of a forward model with v being the number of measurement points, i.e. a number of nodes of FE mesh on a subset Γ v of boundary ∂Ω that is being observed and w is the number of experiments.The a priori measured quantity is stored in vector u r ∈ R vw .In all cases the reference field θ r (x) = θ 0 (x) = 1 for x ∈ Ω.From our experience, the most stable and simple choice of hyper-parameter λ k is the one utilized by Levenberg-Marquardt regularisation (LMR) [8], which is gradually decreasing the parameter during iteration and takes following form Jacobian J k is updated in each iteration and was calculated numerically in a following way where J i is the third-order tensor in i-th iteration, indexes jk are representing measurement nodes in a finite element mesh and individual measurements respectively.Index l denotes a conductivity change on l-th FE element.The tensor is calculated numerically and matricised along indexes jk for convenience of calculation purposes.In general θ can represent multiple parameter fields and in such case equations ( 3) to (7) are evaluated for each parameter individually.

Additional numerical treatment
The iteration scheme (3) does not include any supplementary constrains which provide the solver robustness.Therefore one have to adopt further numerical treatment to maintain the assumptions 1 and 2 valid within the iteration.Introduction of the following two operators increased the stability and robustness of results substantially.
The constrain for strict positivity of conductivity field during the iteration is treated with an ad-hoc algorithm which is not allowing the occurrence of negative numbers in the solution and gives rise to a positivity operator, i.e.T p in a following way ) Since the problem is inherently non-linear and the iteration scheme proceeds with non-restricted linear steps we introduce a cap operator T c limiting the resulting step size to a chosen value k in a following way The final modified GN iteration scheme after applying ( 8) and ( 9) operators to (3) then takes following form

Transient model
Since the numerical model is repeatedly used in each iteration of (10) it plays a fundamental role in the entire inverse procedure.In our case it will to some extent also substitute an actual experiment.Numerical solution is obtained using Finite Element Method space-time discretisation [9][10][11][12].
The main idea of employing the Calderón's technique is created upon foundational research conducted by K. S. Chen [13] in 1989, who proposed and experimentally validated so called Complete Electrode Model (CEM) which allows to consider electrodes actual physical size and introduces a contact impedance layer between the skin and the electrode.Further proof for CEM solution uniqueness was given by E. Somersalo, et al. in 1992 [14].Significant contributions to the theory were papers [15,16] proving the solution is uniquely determined in dimension n ≥ 3 for complete Dirichlet-to-Neumann or Neumann-to-Dirichlet data.Unlike in n ≥ 3, in dimension n = 2 the problem is undetermined and the solution was proved to be unique up to a change of coordinates [17,18].
Specific choice of this model is motivated by real application on historical and other structures, where the standard procedure of taking samples is too difficult or the structure should not be weakened or damaged.Also in real conditions it is not an easy task to sustain a stable and steady state conditions.Not only the surrounding temperature will fluctuate even in a laboratory environment, but for standard building materials like bricks, concrete, wood, etc. the steady state, after changing the loading conditions, can be reached after several hours or days depending on the volumetric capacity, heat conductivity and material thickness.Therefore, we intend to apply the identical principles used in a Calderón problem for transient models.
To capture a time dependent heat transfer, one can adopt following set of equations where ρ s is volumetric mass density, c p is specific heat capacity and ∂Ω (i) (N,T,D) are disjoint subsets of boundary ∂Ω (i) in i-th loading condition with corresponding environmental factors u 0 , α (i) and f N,T,D .
In order to maintain Neumann-to-Dirichlet sensing and uniqueness of the inverse problem, the set of equations in ( 11) is subjected to following constrain Assumption 1.Let Γ m be a subset of boundary ∂Ω that is subjected to measurements.Then must hold, i.e. the boundary subjected to measurements must contain at least some Neumann conditions.
Another set of constrains similar to the ones in [14] is represented by following assumptions Assumption 2. The material parameters λ s , ρ s , c p and transfer coefficients α (i) satisfy following In a definition of this model (11), one can notice that there is no mention of electrodes, thermocouples or measurements indicating our intention not to consciously intervene in the system itself, but only to rely on external influences and natural fluctuation of temperature throughout the day and night.For systems where their boundaries ∂Ω are not exposed to different external influences, one can adapt a similar techniques to Complete Electrode Model, i.e. to attach a device that can excite a different boundary conditions, or control the ambient temperature in the second or third condition in equation (11).
A special feature of this model lies in its independence on excitation device, i.e. stimulation electrodes or thermocouples, and therefore has the least requirements on equipment and the full measurement setup consists of arrays of thermometers and/or thermal cameras.

Numerical solution of the forward model
The forward operator F consists of two independent parameters and solves the system for all time steps N t at once, giving rise to a following formulation where the inputs are heat conductivity λ s and volumetric capacity , where M is number of measurements within single time step and N t is number of time steps, is obtained as a subset of complete solution vectors u i , i = 1, . . .N t from a following set of linear equations where ∆t = t i+1 − t i is a time step, τ ∈ 0; 1 is a time integration parameter1 and the system matrices and right hand side vectors are computed in a following way For detailed information and derivation of individual terms, the interested reader is referred to [9][10][11][12] and literature therein.

Results
We consider two finite element meshes for each physical domain, i.e. one is used to simulate real measurements and the second is used for the reconstruction algorithm.Both finite element meshes are shown on Figure 1.Further improvements were made by integrating measurement errors in the calculation, specifically we assumed only the direct measured quantity is subjected to errors in a following way where u me is some measured quantity with the error included, u m is its true value and ∼ N µ, σ 2 is the error.We assume the sensor is calibrated and its deviation is known to fall within 3•σ region, i.e. µ = 0 and deviation d m = 0.25 o C given by manufacturer implies σ = d m /3.All the sensor errors are considered to be mutually independent and identically distributed random variables, although using thermal camera would rather impose spatially correlated random variables forming a random field.We assume changes in temperature on the observed boundary over several days can provide sufficient amount of data for the algorithm to recover material fields θ = [λ s ; c v ], i.e thermal conductivity and volumetric capacity respectively.In case the observed structure is located in an environment without sufficient changes in loading conditions, one may use a device that excites a different boundary conditions, e.g. a heater or air conditioning system.The measure of reconstruction algorithm performance is represented by an error on the reconstructed material field in comparison to the original field in a following way where θ i,r is the i-th recovered material field obtained from (10) and θ i,o is the i-th original material field.
In order to simulate a natural environment conditions and capture temperature fluctuations in the interior and exterior, we use a realistic measurements shown on Figure 2. The environment was in total monitored for eleven consecutive days at minute intervals.For the calculation purposes data were further sparsified to one hour intervals which satisfy both sufficient precision in describing the temperature curve and reasonable computational load.Boundary conditions were set accordingly to Figure 2  (

as follows
The magnitude of reference fields were set according to regular building materials, i.e.
In Figure 2, the measurement starts from day two and the observation period lasts 9 day.Measurement is taken every hour which yields into 216 measurements of boundary Γ 122 = 2 i=1 ∂Ω i .For purpose of this contribution a set of 500 realisations of measured temperatures was given to the reconstruction algorithm.Resulting mean and standard deviation of reconstructed fields are displayed on following figures

Conclusions
In oppose to EIT where the crucial part of successful material field recovery is a precise placement of measurement electrodes together with boundary shape and knowledge of loading conditions2 , in our application the problems might arise from insufficient environmental factors variability leading to imperfections and unwanted artefacts in the reconstructed fields.This can be identified although it can be remedied by a use of an external devices providing different boundary conditions.
In a simple case shown on Figures 3 and 4 one can observe the possibility of recovering material parameters of such system in the first place.The difference between Figures 3a and 4a     in Figures 3a and 3b or 4a and 4b, i.e. the standard deviation is in both cases approximately two orders of magnitude lower than the mean.Our intention is then to extend the existing model with imperfections or errors on all the other levels, i.e. uncertainty in ambient temperature, transfer coefficient or a geometrical level inaccuracies.

Figure 1 .
Figure 1.Finite elements meshes: coarse (top) and fine (bottom).Observed boundary Γm depicted in red colour, where m = 122 nodes of finite element mesh being measured.

Figure 2 .
Figure 2. Interior and exterior temperatures.Preceding temperature in grey, actual observation period in colour.

Figure 3 .
Figure 3. Mean (A) and standard deviation (B) of reconstructed conductivity field.In grey: original field, in colour: reconstructed field.

Figure 4 .
Figure 4. Mean (A) and standard deviation (B) of reconstructed capacity field.In grey: original field, in colour: reconstructed field.