2 RADAR TOMOGRAPHY WITH TIME DOMAIN SIGNALS
2.2 The forward model
The forward model is formulated based on the hyperbolic wave equation 2.13. The formulation is shown here for the two-dimensional case where the transmitted sig-nal produces a transverse electric field perpendicular to the direction of propagation.
The reason for choosing to present the two-dimensional case is that the higher-order Born approximation developed during this thesis (Publication III) was formulated for the two-dimensional case. The three-dimensional formulation of the forward problem including the polarization terms resulting from electromagnetic wave in-teraction with medium interfaces is given in [74]. To match the notations in the original publications introducing the forward model[64, 74, 75], and Publications I, III, V, and VI of this thesis, the electric fieldE⃗ is denoted by the scalar potential u. The antenna current densityJ⃗s acting as the source function for the radar pulse signal transmitted from the antenna is given by f(t). The forward model can hence be expressed as
ϵµ∂2u
∂t2 +σµ∂u
∂t −∆u=µ∂ f
∂t. (2.30)
The equation 2.30 is sought in a form emphasizing the recovery of the real part of the relative electric permittivityϵr. By expandingϵ=ϵrϵ0(Eq. 2.24), and assuming that the magnetic permeability of the mediumµ=µ0, it is written as
ϵrϵ0µ0∂2u
∂t2 +σµ0∂u
∂t −∆u=µ0∂ f
∂t. (2.31)
Table 2.1 Spatial scaling of SI-units to the unitless computational domain. The spatial scaling factor iss and the wave velocityc=1/⎷ϵ0µ0=1whenϵr =1. The physical constantsϵ0= 8.85·10−12F/m is the electric permittivity, andµ0=4π·10−7H/m the magnetic permeability of vacuum.
Quantity Unitless SI-units
Time t (ϵ0µ0)1/2s t
Position ⃗x s x⃗
Velocity c=1 (ϵ0µ0)−1/2c Frequency f (ϵ0µ0)−1/2s−1f Electric permittivity ϵr ϵr
Electric conductivity σ (ϵ0/µ0)1/2s−1σ
To simplify the numerical parameters, the computational problem is solved with a unitless set of parameters of time t, frequency f, space x⃗, electric permittivity ϵr, conductivityσ and velocityc by choosing that the wave velocity in free space (ϵr =1) isc= (ϵ0µ0)−1/2=1.
By multiplying both sides of the equation 2.31 by(ϵ0µ0)−1 = 1 to obtain the unitless system, and by simplifying the resulting equation, the forward model takes the form
ϵr∂2u
∂t2 +σ∂u
∂t −∆u=∂ f
∂t for all(t,x⃗)∈[0,T]×Ω (2.32) with the initial conditions u|t=0 = u0 and (∂u/∂t)|t=0 = u1. The signal source
∂ f/∂t is given by
∂ f(t,x⃗)
∂t =δp⃗(x⃗)∂ f˜(t)
∂t , (2.33)
where the notation f˜(t) refers to the time-dependent part of f, andδp⃗(⃗)x is the Dirac’s delta function with respect to the signal transmission point ⃗p.
The quantities which contain spatial dependency are scaled by the spatial scaling factor sobtained by the ratio of the real SI-unit size of the target object divided by the unitless size of the target domainD. The parameters in the equation 2.32 can be scaled to SI-units according to Table 2.1.
By defining an auxiliary function⃗gas g⃗=
∫︂ t
0
∇u(τ,x⃗)dτ, (2.34) the hyperbolic wave equation 2.32 can be expressed as the first-order differential sys-tem
ϵr∂u
∂t +σu− ∇ ·g⃗= f (2.35)
∂⃗g
∂t − ∇u=0 (2.36)
in the domainΩ×[0,T]. To obtain the weak form of this system, the equation 2.35 is multiplied by the test functionv:[0,T]→H1(Ω), and the equation 2.36 by the test functionw⃗ :[0,T]→L2(Ω)followed by integration by parts yields the weak formulation of the system:
∂
∂t
∫︂
ΩϵruvdΩ+
∫︂
ΩσuvdΩ+
∫︂
Ωg⃗· ∇vdΩ=
⎧
⎨
⎩
f˜(t) if⃗x=p⃗
0 otherwise (2.37)
∂
∂t
∫︂
Ωg⃗·w⃗dΩ−
∫︂
Ωw⃗· ∇udΩ=0. (2.38) This weak form has the unique solutionu :[0,T]→H1(Ω)when the domain and the parameters are regular enough[27].
2.2.1 Evaluation of wave propagation in a domain
The problem of full wave propagation in the domain Ω can be associated with a situation where a signal ˜f(t)is transmitted from the point p⃗1 and received at the point p⃗2. The received signald˜(t)at the receiver point ⃗p2can be expressed as the linear convolution
d˜(t) =G⃗p1,⃗p2∗f˜(t) =
∫︂∞
−∞
Gp⃗1,p⃗2(t−τ)˜f(τ)dτ, (2.39)
whereG⃗p1,⃗p2 is the Green’s functionG =G(ϵ,ε), a functional of the permittivityϵ and a nuisance parameterε, and representing the impulse response to an infinitely short monopolar pulse transmitted from p⃗1 and received in p⃗2. The nuisance pa-rameter ε is related to the uncertainties in the numerical model, including signal attenuation due to, for example, conduction losses, reflections and refractions, and also numerical modelling errors. These uncertainties can be modelled with a an ad-ditive error termεresulting in a simplified expression for the Green’s function as
G(ϵ,ε) =G(ϵ) +ε. (2.40) The equation 2.39 expresses the wavefield in a domain with constant permittivity distributionϵ. In realistic scenarios, the domain includes electromagnetic scatterers, perturbationsρto the overall permittivity distributionϵ. To investigate the effect of such small scattering obstacles, assume one at pointr⃗causing a local perturbation ρin the overall permittivity distributionϵ. The permittivity at that point is then defined as
ϵr⃗=ϵ+ρ. (2.41)
An electromagnetic wavefield signal˜ transmitted from the pointf ⃗p1and travel-ling through the point⃗r which contains the permittivity perturbation, is altered at that point. The point⃗rcan therefore be considered to act as a new point source trans-mitting an altered signal h˜. Because the wavefield is altered only at ⃗, the Green’sr function equalsG(ϵ)elsewhere. The perturbed wavefieldd˜ received at the point⃗p2 is hence
d˜=G(ϵ)⃗p1,p⃗2∗˜f+c(ρ)G(ϵ)⃗r,p⃗2∗h˜, (2.42) wherec(ρ)is a constant contrast factor whose numerical dependence on the pertur-bationρis determined later in the section 2.3.2, and the wavefield
h˜=G(ϵ⃗r)⃗p1,r⃗∗˜.f (2.43) Due to the uncertainties included in the error termε(Eq. 2.40), the exact form of the Green’s function is not known. However, it is approximated numerically by using a regularised deconvolution process (Figure 2.1) which proceeds according to
Figure 2.1 A schematic presentation of regularised deconvolution which is applied in computing wave propagation between points⃗p1and⃗p2where there is a scattering detail altering the wave-field at point⃗r. Adapted from Publication III. Reprinted with permission.
the following steps:
1. Propagate a single wave from both p⃗1 and p⃗2 and evaluate the terms h˜ = G(ϵ)⃗p1,r⃗∗˜ andf p˜=G(ϵ)p⃗2,r⃗∗˜ for each scattering pointf ⃗.r
2. Use the Tikhonov-regularised deconvolution with a suitably chosen regulari-sation parameterνto estimate the Green’s function
g˜=G(ϵ)r⃗,⃗p2 =G(ϵ)⃗p2,r⃗. (2.44) 3. Estimate the wavefield originating fromr⃗at⃗p2byd˜=g˜∗˜.f
The principle of reciprocity of wave propagation ensures that the equation 2.44 holds. The entries of the vectors˜,f h˜,p˜,g˜, andd˜ contain the pointwise time evolu-tion of the corresponding variables. Backscattering data is obtained when p⃗1=⃗p2. 2.2.2 Born approximation and higher-order scattering
In a complex, bounded geometry, a propagating wavefield may experience multiple scattering events related to the same perturbation. Such events can be modelled with the Born series. The first-order term of the series is called Born approximation, where the total field is replaced by the incident field and this incident field is assumed to be the driving field at each point in the scatterer (Figure 2.2). Therefore, the Born approximation (BA), or the first-order term of the Born series, can be obtained by substitutingG(ϵr⃗)in the equation 2.43 with the corresponding Green’s function of the incident fieldG(ϵ):
hBA,1˜ =G(ϵ)p1⃗,r⃗∗˜.f (2.45)
Figure 2.2 Schematic representation of the order BA (left) and second-order BA (right). The first-order BA accounts for the scattered wavefields (black, dashed arrows) which originate from a single point and its interaction withr⃗with respect to the details in the computational geometry.
The second-order BA accounts for the wavefields which have scattered twice from⃗r (solid grey arrows). Adapted from Publication III. Reprinted with permission.
The higher-order terms of the series modelling the secondary, tertiary, and higher order scattering events from the same scatterer are derived analogously to the equa-tion 2.42 from the recursive equaequa-tion
h˜BA,n=G(ϵ)p1⃗,r⃗∗˜f+c(ρ)G(ϵ)r⃗,r⃗∗h˜BA,n−1, (2.46) wheren gives the order of the scattering event. As the first-order Born approxima-tion is based onG alone, the nonlinear wave propagation effects are only recovered by the higher order terms where the total field is replaced by the incident field as given by the equation 2.43.
2.2.3 Discretisation of the forward problem
The forward problem 2.32 and its weak form are solved numerically by the finite element time domain (FETD) method, where the spatial d-dimensional (d ={2, 3}) domainΩis discretised into a d-simplex finite element meshT consisting of m ele-ments{T1,T2, . . . ,Tm}. Each of these elementsTi,i=1, . . . ,mis associated with an element indicator functionχi∈L2(Ω). Here, the discretisation is presented for the two-dimensional case. The meshT consists of a set ofnmesh nodes{r⃗1,r⃗2, . . . ,⃗rn} identified with piecewise linear nodal basis functions{ϕ1,ϕ2, . . . ,ϕn} ∈H1(Ω). The potential and gradient fields in equations 2.35 and 2.36 are approximated in each dimension of the meshT as finite sums
u=
n
∑︂
j=1
pjϕj and g⃗=
d
∑︂
k=1
g(k)e⃗k, (2.47)
where g(k) = ∑︁m
i=1qi(k)χi, the weighted sum of the element indicator functions {χ1,χ2, . . . ,χm} ∈L2(Ω).
By defining the test functionsv :[0,T]→(V)⊂ b1(Ω)and w⃗ :[0,T]→ W ⊂ L2(Ω)withV =span{ϕ1,ϕ2, . . . ,ϕn}andW =span{χ1,χ2, . . . ,χm}, the weak form (Eqs 2.35-2.36) is expressed in the discretised Ritz-Galerkin form[15]:
∂
∂tCp+Rp+Sp+∑︂d
k=1
B(k)Tq(k)=f (2.48)
∂
∂tAq(k)−B(k)p+T(k)q(k)=0, (2.49) wherep= (p1,p2, . . . ,pn)andq(k)= (q1(k),q2(k), . . . ,qm(k))are the coordinate vectors associated with the finite sums in the equations 2.47 which contain the scalar poten-tial fieldu, and its gradientsg⃗. The elements in the matricesA∈Rm×m,B∈Rm×n, C∈Rn×n,S∈Rn×n, andT∈Rm×mare defined in the two-dimensional case as:
Ci,j =
∫︂
ΩϵrϕiϕjdΩ (2.50)
Ri,j =
∫︂
ΩσϕiϕjdΩ (2.51)
Bi,j =
∫︂
Ti
⃗ek· ∇ϕjdΩ (2.52)
fi=
∫︂
Ω
fϕidΩ (2.53)
Ai,i=
∫︂
Ti
dΩ Ai,j =0 ifi̸= j (2.54)
Si,j =
∫︂
Ωξ ϕiϕjdΩ (2.55)
Ti,i(k)=
∫︂
Ti
ζ(k)dΩ Ti,(k)j 0 ifi̸= j. (2.56)
The matricesSand T(k) are additional matrices associated with the split-field per-fectly mached layer which is defined for the outermost part of the computational domainΩto eliminate reflections from the boundary∂Ωback to the inner part of
Ω. The outermost part is defined as{x⃗∈Ω|ϱ1≤maxk|xk| ≤ϱ2}. The parametersξ andζ(k)are absorption parameters of the formξ(x⃗) =ς, whenϱ1≤maxk|xk| ≤ϱ2, andζ(k)(x⃗) =ς, whenϱ1≤ |xk| ≤ϱ2. In other parts of the modelξ(x⃗) =ζ(k)(x⃗) = 0.
The time domain[0,T]is discretized by the standard finite difference approach into∆t-spaced regular grid ofN time points. The temporal discretization of the spatially discretized weak formulation equations 2.48 and 2.49 yields the leapfrog time integration system:
pℓ+1=pℓ+∆tC−1
fℓ−Rpℓ−Spℓ−
∑︂d k=1
B(k)Tq(k)
ℓ+12
(2.57) q(k)
ℓ+12 =q(k)
ℓ−12 +∆tA−1
B(k)pℓ−T(k)q(k)
ℓ−12
, (2.58)
in whichℓ=1, 2, . . . ,N are the time points used to simulate the signal propagation in the spatial domainΩ, andq(k)
ℓ−12 is the gradient ofpℓintegrated over time.
By defining the auxiliary time-evolution vectors
a(k)
ℓ−12 =A−1
B(k)pℓ−T(k)q(k)
ℓ−12
(2.59) bℓ+1
2 =−Rpℓ−Spℓ−
∑︂d k=1
B(k)Tq(k)
ℓ+12, (2.60)
the equations 2.57 and 2.58 can be expressed in a more concise form, emphasising the permittivity-containing mass matrixC. The electromagnetic field is hence solved in the Cartesian components by the leapfrog iteration
pℓ+1=pℓ+∆tC−1
fℓ+bℓ+1 2
(2.61) q(k)
ℓ+12 =q(k)
ℓ−12+∆ta(k)
ℓ−12. (2.62)
2.2.4 Formulation of the higher-order Born approximation
The matrix C in the leapfrog formulation of the time-evolution of the wavefield propagation (Eq. 2.61) is a permittivity-weighted positive definite mass matrix which entries are given by the equation 2.50. To formulate this matrix in the similar way the point-wise permittivity is formulated as the sum of the overall permittivity distribu-tion and perturbadistribu-tion (Eq. 2.41), the matrixCcan be decomposed as the difference
C=C1−C2, (2.63)
whereC1andC2correspond toϵandρ, respectively. When the perturbationρis small enough so that||C−11 C2||<1, the inverse ofCcan be expanded as a geometric series for(I−C−11 C2)−1as
C−1= (C1−C2)−1= (I−C−11 C2)−1C−11
= (I+C−11 C2+C−21 C22+. . .)C−11 . (2.64) Substituting the first-degree polynomial approximationC−1≈(I+C−11 C2)C−11 of the equation 2.64 in the leapfrog equation 2.61, the numerical leapfrog iteration for-mulae (Eqs 2.61 and 2.62) are formulated as
pℓ+1=pℓ+∆tC−11
hℓ+fℓ+bℓ+1
2
(2.65) q(k)
ℓ+12 =q(k)
ℓ−12 +∆ta(k)
ℓ−12, (2.66)
where the termhℓ=C2C−11 (fℓ+bℓ+1
2)is an auxiliary source vector. It can be inter-preted as the source vector originating from the perturbation point and hence as the first-order BA of the system. Therefore, it is labelled as
hBA,1ℓ =C2C−11 (fℓ+bℓ+1
2). (2.67)
The higher order Born approximations are obtained by applying equation 2.67 recursively. Physically, this corresponds to replacing the total field with the incident field at each step at the point of perturbation. The auxiliary source term for the
higher-order BA is thus given by
hBA,nℓ =C2Pn(C−11 C2)C−11 (fℓ+bℓ+1
2), (2.68)
where Pn(C−11 C2) =I+C−11 C2+· · ·+C−n1 +1C2n−1is a matrix-valued polynomial.
The detailed derivation of the equation 2.68 and its convergence is shown in Publi-cation III.
2.2.5 Multiresolution approach for the forward and inverse solvers To speed up the forward model computation, a multiresolution approach[15, 64, 75] is adopted. The fine mesh T (section 2.2.3) with n nodes is used as basis to recover a coarser, nested d-simplex meshT′withM elements{T1′,T2′, . . . ,TM′}and N nodes{r⃗1,r⃗2, . . . ,⃗rN}, which are shared by the fine mesh (N < n). While the resolution of the fine meshT is determined by the geometrical constraints of the forward simulation such as the domain structure and the applied wavelength, the resolution of the coarse mesh T′ is chosen based on the desired precision of the reconstruction.
The electric permittivity distribution of the domain is sought in the form
ϵ=ϵb g+ϵρ, (2.69)
where ϵb g is a constant, fixed, background distribution and ϵρ a variable pertur-bation. The permittivity distribution ϵwithin the domain is defined with respect to the elements of the coarse mesh T′and is given by ϵ = ∑︁M
j=1cjχj′. Assuming a piecewise constant permittivity distribution in the fine meshT, the permittivity distribution of the fine mesh is given byϵ=∑︁m
j=1cjχj. Such a nested mesh structure leads to a natural conclusion that the maximal theoretical reconstruction resolution is lower than that which is used in propagating a wave through the domain.