Acta Polytechnica CTU Proceedings


doi:10.14311/APP.2017.13.0043
Acta Polytechnica CTU Proceedings 13:43–48, 2017 © Czech Technical University in Prague, 2017

available online at http://ojs.cvut.cz/ojs/index.php/app

NON-INVASIVE INVERSE PROBLEM IN CIVIL ENGINEERING

Jan Havelka∗, Jan Sýkora, Anna Kučerová

Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 2077/7, 166 29 Prague 6, Czech
Republic

∗ corresponding author: jan.havelka.1@fsv.cvut.cz

Abstract. 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.

Keywords: Calderón problem, Finite Element Method, Boundary inverse method, Neumann-to-
Dirichlet map, Heat transfer.

1. Introduction
Non-invasive methods are gaining increasing interest
in various disciplines spanning from medicine to geo-
physics. 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, of-
ten 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 de-
scription 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 Tomog-
raphy (EIT). The mathematical procedure was first
rigorously described by an Argentinian mathemati-
cian 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 elec-
tric 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 distri-
bution.
In this contribution we will employ similar tech-

niques of Calderón inverse problem paradigm to
a parabolic partial differential equation with two in-
dependent 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.

2. Inverse problem
Parameters recovery of a given model is often regarded
to as ill-posed problem [4, 5]. For such task it is typi-
cal 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 insuffi-
ciency 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 informa-
tion, 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 regular-
isation 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

θ̃ = arg min
θ
‖F (θ)− um‖p +G(θ), (1)

43

http://dx.doi.org/10.14311/APP.2017.13.0043
http://ojs.cvut.cz/ojs/index.php/app


J. Havelka, J. Sýkora, A. Kučerová Acta Polytechnica CTU Proceedings

where F (θ) is output of the forward operator, um 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]

G(θ) = λ2‖L(θ − θr)‖2`2
, (2)

where the hyper-parameter, λ controls the trade-off
between solution stability, given constrains and a dis-
tance from the actual solution. Additional parame-
ters θ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 cho-

sen to be `2-norm due to the convenience for com-
putational purposes and preference of smooth solu-
tions [7]. In such case, by combining the minimi-
sation 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]

θk+1 = θk + δθk, (3)

δθk =
(
JTk Jk + λ2R

)−1 (JTk δu− λ2R (θk − θr)
)
,
(4)

where Jk is the Jacobian evaluated at θk, δu is a vec-
tor containing the difference between the model out-
put and real measurements, θr is a conductivity field
containing its prior assumptions and R = LTL is reg-
ularisation operator composed from a pre-calculated
Laplacian L of a piece-wise constant functions on
a finite element mesh.

δu = ur − F (θk) , (5)

where F (θ) ∈ Rvw 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 ur ∈ Rvw. 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

λk = max
(
max

(
JTk Jk

))
. (6)

Jacobian Jk is updated in each iteration and was
calculated numerically in a following way

J
(jkl)
i = ∂ujk

∂θli
, (7)

where Ji 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 mul-
tiple 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 supple-
mentary constrains which provide the solver robust-
ness. 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. Tp in a following way

Tp(θ(x)) =
{
θ(x) > 0→ x+, θ(x+) = θ(x+),
θ(x) < 0→ x−, θ(x−) = min(θ(x+)).

(8)
Since the problem is inherently non-linear and the

iteration scheme proceeds with non-restricted linear
steps we introduce a cap operator Tc limiting the
resulting step size to a chosen value k in a following
way

Tc(θ(x)) =
{
θ(x) > k, θ(x) = sign(θ(x)) · k,
θ(x) < k, θ(x) = θ(x).

(9)
The final modified GN iteration scheme after apply-

ing (8) and (9) operators to (3) then takes following
form

θk+1 = Tp (θk + Tc (δθk)) . (10)

3. 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. Numeri-
cal solution is obtained using Finite Element Method
space-time discretisation [9–12].
The main idea of employing the Calderón’s tech-

nique is created upon foundational research conducted
by K. S. Chen [13] in 1989, who proposed and experi-
mentally validated so called Complete Electrode Model
(CEM) which allows to consider electrodes actual phys-
ical size and introduces a contact impedance layer

44



vol. 13/2017 Non-invasive Inverse Problem in Civil Engineering

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 ap-
plication 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 sus-
tain 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

ρs(x)cp(x)∂u
∂t

(x, t)−∇ · (λs(x)∇u(x, t)) = f̃(x, t),

n(x) · (λs(x)∇u(x, t)) = f̃N (x, t),

α(i)(u(x, t)− u0(x, t)) = f̃T (x, t),

u(x, t) = f̃D(x, t),

t > 0, ∂Ω(i) = ∂Ω(i)
N ∪ ∂Ω(i)

T ∪ ∂Ω(i)
D , i = 1, . . . , I,

(11)
where ρs is volumetric mass density, cp is specific heat
capacity and ∂Ω(i)

(N,T,D) are disjoint subsets of bound-
ary ∂Ω(i) in i-th loading condition with corresponding
environmental factors u0, α(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

Ψ =
(
∂Ω(i)

N ∪ ∂Ω(i)
T

)
: (Ψ ∩ Γ) /∈ ∅, ∀i,

must hold, i.e. the boundary subjected to measure-
ments 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, cp
and transfer coefficients α(i) satisfy following

(i) λs ∈ L∞(Ω;R), inf
x∈Ω

λs(x) = λ+
s > 0,

(ii) 0 < ρ−s ≤ ρs ≤ ρ+
s <∞,

(iii) 0 < c−p ≤ cp ≤ c+p <∞,

(iv) 0 < α
(i)
− ≤ α(i) ≤ α(i)

+ <∞, ∀i.

In a definition of this model (11), one can notice
that there is no mention of electrodes, thermocou-
ples 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 simi-
lar 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 indepen-

dence on excitation device, i.e. stimulation electrodes
or thermocouples, and therefore has the least require-
ments 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
Nt at once, giving rise to a following formulation

F (λs, cv) = um, um ∈ RM ·Nt , (12)

where the inputs are heat conductivity λs and volu-
metric capacity cv = ρs · cp. The measurement matrix
um ∈ RM ·Nt , where M is number of measurements
within single time step and Nt is number of time steps,
is obtained as a subset of complete solution vectors ui,
i = 1, . . . Nt from a following set of linear equations

KCui+1 =
(

KC −∆tK̃
)

ui−1 + . . .

∆t
(

(1− τ)f̃n,i−1 + τ f̃n,i
)

+ . . .

∆t
(

(1− τ)f̃t,i−1 + τ f̃t,i
)

+ . . .

∆t
(

(1− τ)f̃d,i−1 + τ f̃d,i
)
,

(13)

K̃ = K + Kt, (14)

KC = C + τ∆tK̃, (15)

where ∆t = ti+1− ti 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

Kkl =
∫

Ω
λs∇φk · ∇φl dA k, l = 1 . . . N, (16)

1The foregoing computations were calculated using Crank-
Nicolson integration method, i.e. τ = 0.5

45



J. Havelka, J. Sýkora, A. Kučerová Acta Polytechnica CTU Proceedings

Ckl =
∫

Ω
ρscpφkφl dA k, l = 1 . . . N, (17)

Kt,kl =
∫
∂ΩT

αφkφl dS k, l = 1 . . . N, (18)

f̃t,li =
∫
∂ΩT

αu0,iφl dS i = 1 . . . Nt, l = 1 . . . N,

(19)

f̃n,li = −
∫
∂ΩN

f̃N (ti)φl dS i = 1 . . . Nt, l = 1 . . . N,

(20)

f̃d,i = − (K + Kt) f̃D(ti) i = 1 . . . Nt. (21)

For detailed information and derivation of individ-
ual terms, the interested reader is referred to [9–12]
and literature therein.

4. Results

We consider two finite element meshes for each physi-
cal domain, i.e. one is used to simulate real measure-
ments and the second is used for the reconstruction
algorithm. Both finite element meshes are shown on
Figure 1. Further improvements were made by inte-
grating measurement errors in the calculation, specifi-
cally we assumed only the direct measured quantity
is subjected to errors in a following way

ume = um + ε, (22)

where ume is some measured quantity with the error
included, um is its true value and ε ∼ N

(
µ, σ2) is the

error. We assume the sensor is calibrated and its devi-
ation is known to fall within 3·σ region, i.e. µ = 0 and
deviation dm = 0.25 oC given by manufacturer implies
σ = dm/3. All the sensor errors are considered to be
mutually independent and identically distributed ran-
dom variables, although using thermal camera would
rather impose spatially correlated random variables
forming a random field.

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.

We assume changes in temperature on the observed
boundary over several days can provide sufficient
amount of data for the algorithm to recover mate-
rial fields θ = [λs; cv], i.e thermal conductivity and
volumetric capacity respectively. In case the observed
structure is located in an environment without suf-
ficient changes in loading conditions, one may use
a device that excites a different boundary conditions,
e.g. a heater or air conditioning system. The mea-
sure of reconstruction algorithm performance is repre-
sented by an error on the reconstructed material field
in comparison to the original field in a following way

εθi
= ‖θi,r − θi,o‖2

‖θi,o‖2
, (23)

where θi,r is the i-th recovered material field ob-
tained from (10) and θi,o is the i-th original material
field.
In order to simulate a natural environment con-

ditions and capture temperature fluctuations in the
interior and exterior, we use a realistic measurements
shown on Figure 2. The environment was in total mon-
itored for eleven consecutive days at minute intervals.
For the calculation purposes data were further sparsi-
fied 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

46



vol. 13/2017 Non-invasive Inverse Problem in Civil Engineering

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

λ
∂u

∂n
(x, t)|∂Ω1 = 10 · (u(x, t)− T1(t)) ,

λ
∂u

∂n
(x, t)|∂Ω2 = 10 · (u(x, t)− T2(t)) .

(24)

The magnitude of reference fields were set according
to regular building materials, i.e.

λr = 1 [ W m−1 K−1 ],
cv,r = 1× 106 [ J m−3 K−1 ].

(25)

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 measure-
ments of boundary Γ122 =

⋃2
i=1 ∂Ωi. For purpose

of this contribution a set of 500 realisations of mea-
sured temperatures was given to the reconstruction
algorithm. Resulting mean and standard deviation of
reconstructed fields are displayed on following figures

5. Conclusions
In oppose to EIT where the crucial part of successful
material field recovery is a precise placement of mea-
surement electrodes together with boundary shape and
knowledge of loading conditions2, in our application
the problems might arise from insufficient environ-
mental 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 parame-
ters of such system in the first place. The difference
between Figures 3a and 4a shows that the variation
of temperature in time was not sufficient to fully re-
construct the volumetric capacity field cv resulting
into higher reconstruction error of E[εcv

] = 0.2641 in
oppose to E[ελs ] = 0.1414. Looking closer at devia-
tions of errors std[εcv

] = 0.0014 and std[ελs
] = 0.0017

2Especially in EIT, where the electrodes are placed on a hu-
man body, the boundary shape and electrodes are in continuous
motion as the patient breath and move.

(a) . E[ελs ] = 0.1414 [-]

(b) . std[ελs ] = 0.0017 [-]

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

one can notice the insensitivity of reconstruction al-
gorithm to induced error on measurements which can
be also seen by comparing the magnitudes of plots

47



J. Havelka, J. Sýkora, A. Kučerová Acta Polytechnica CTU Proceedings

(a) . E[εcv ] = 0.2641 [-]

(b) . std[εcv ] = 0.0014 [-]

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

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 ambi-
ent temperature, transfer coefficient or a geometrical
level inaccuracies.

Acknowledgements
The financial support of this research by the GA15-07299S
and SGS17/042/OHK1/1T/11 is gratefully acknowledged.

References
[1] A. P. Calderón. On an inverse boundary value
problem. Computational & Applied Mathematics 25:133
– 138, 2006.

[2] L. B. Slichter. The interpretation of the resistivity
prospecting method for horizontal structures. Physics
4(9):307–322, 1933. doi:10.1063/1.1745198.

[3] R. E. Langer. An inverse problem in differential
equations. Bulletin of the American Mathematical
Society 39(10):814–820, 1933.

[4] M. Vauhkonen. Electrical impedance tomography and
prior information. Phd thesis, Kuopio University, 1997.

[5] A. Kirsch. An Introduction to the Mathematical
Theory of Inverse Problems. Applied Mathematical
Sciences 120. Springer-Verlag New York, 2nd edn., 2011.

[6] D. S. Holder. Electrical Impedance Tomography:
Methods, History and Applications. Series in Medical
Physics and Biomedical Engineering. Taylor & Francis,
1st edn., 2004.

[7] Y. Mamatjan, A. Borsic, D. Grsoy, A. Adler.
Experimental/clinical evaluation of EIT image
reconstruction with `1 data and image norms. Journal
of Physics: Conference Series 434(1):012078, 2013.

[8] H. Niu, P. Guo, L. Ji, et al. Improving image quality
of diffuse optical tomography with a
projection-error-based adaptive regularization method.
Optics express 16(17):12423–12434, 2008.

[9] K. J. Bathe. Finite element procedures. Klaus-Jurgen
Bathe, 2006.

[10] T. J. R. Hughes. The finite element method: linear
static and dynamic finite element analysis. Courier
Corporation, 2012.

[11] G. Strang, G. J. Fix. An analysis of the finite
element method, vol. 212. Prentice-hall Englewood
Cliffs, NJ, 1973.

[12] O. C. Zienkiewicz, R. L. Taylor. The finite element
method, vol. 3. McGraw-hill London, 1977.

[13] K. S. Cheng, D. Isaacson, J. C. Newell, D. G. Gisser.
Electrode models for electric current computed
tomography. IEEE Transactions on Biomedical
Engineering 36(9):918–924, 1989. doi:10.1109/10.35300.

[14] E. Somersalo, M. Cheney, D. Isaacson. Existence and
uniqueness for electrode models for electric current
computed tomography. SIAM Journal on Applied
Mathematics 52(4):1023–1040, 1992.
doi:10.1137/0152060.

[15] R. M. Brown, G. A. Uhlmann. Uniqueness in the
inverse conductivity problem for nonsmooth
conductivities in two dimensions. Communications in
Partial Differential Equations 22(5-6):1009–1027, 1997.
doi:10.1080/03605309708821292.

[16] J. Sylvester, G. Uhlmann. A global uniqueness
theorem for an inverse boundary value problem. Annals
of Mathematics 125(1):153–169, 1987.

[17] J. Sylvester, G. Uhlmann. A uniqueness theorem for
an inverse boundary value problem in electrical
prospection. Communications on Pure and Applied
Mathematics 39(1):91–112, 1986.
doi:10.1002/cpa.3160390106.

[18] A. I. Nachman. Global uniqueness for a
two-dimensional inverse boundary value problem.
Annals of Mathematics 143(1):71–96, 1996.

48

http://dx.doi.org/10.1063/1.1745198
http://dx.doi.org/10.1109/10.35300
http://dx.doi.org/10.1137/0152060
http://dx.doi.org/10.1080/03605309708821292
http://dx.doi.org/10.1002/cpa.3160390106

	Acta Polytechnica CTU Proceedings 13:44–49, 2017
	1 Introduction
	2 Inverse problem
	3 Transient model
	4 Results
	5 Conclusions
	Acknowledgements
	References