Frequency domain electromagnetic data inversion Caterina Fenu Joint - - PowerPoint PPT Presentation

Frequency domain electromagnetic data inversion

Caterina Fenu

Joint work with Gian Piero Deidda and Giuseppe Rodriguez (University of Cagliari) University of Pisa, Italy

PING Workshop, Florence, April 6, 2016

• C. Fenu

Frequency domain electromagnetic data inversion

Electromagnetic induction techniques

Electromagnetic induction (EMI) techniques are often used for non-destructive investigation of soil properties, as they are affected by electromagnetic properties of the subsurface layers, namely the electrical conductivity σ and the magnetic permeability µ. Knowing such parameters allows one to ascertain the presence of particular substances, with many important applications: hydrological characterizations hazardous waste studies archaeological surveys precision-agriculture unexploded ordnance detection

• C. Fenu

Frequency domain electromagnetic data inversion

Ground conductivity meter

A ground conductivity meter is the basic instrument for EMI. A (GCM) contains two coils (a transmitter and a receiver) placed at a fixed distance. An alternating sinusoidal current in the transmitter produces a primary magnetic field HP, which induces small eddy currents in the subsurface. These currents produce a secondary magnetic field HS, which is sensed by the receiver. The ratio of the secondary to the primary magnetic fields is then used, along with the instrumental parameters, to estimate electrical properties of the subsurface.

• C. Fenu

Frequency domain electromagnetic data inversion

Ground conductivity meter

The coils axes can be aligned either vertically or horizontally with respect to the ground surface, producing different measures. vertical & horizontal alignments

• C. Fenu

Frequency domain electromagnetic data inversion

Ground conductivity meter

Under some assumptions: instrument at ground level (h = 0), in vertical orientation, soil with uniform magnetic permeability µ0 = 4π10−7 H/m, soil with uniform electrical conductivity σ, small induction number B = r

• 1

2µ0ωσ ≪ 1, where r is the inter-coil distance (∼ 1 m), ω = 2πf , f operating frequency (∼ 10 kHz).

• C. Fenu

Frequency domain electromagnetic data inversion

Ground conductivity meter

Under some assumptions: instrument at ground level (h = 0), in vertical orientation, soil with uniform magnetic permeability µ0 = 4π10−7 H/m, soil with uniform electrical conductivity σ, small induction number B = r

• 1

2µ0ωσ ≪ 1, the instruments measures the apparent conductivity m = 4 Im(HS/HP) µ0ωr2 , which, under the above assumptions, coincides with σ.

• C. Fenu

Frequency domain electromagnetic data inversion

The inversion problem

In real applications the assumption of uniform soil conductivity is not realistic. On the contrary, geophysicists are particularly interested in non homogeneous soil. Moreover, apparent conductivity gives no information on the depth localization of inhomogeneities. To face the problem of data inversion multiple measures are needed to recover the distribution of conductivity with respect to depth.

• C. Fenu

Frequency domain electromagnetic data inversion

Parameters which influence the GCM measures

• h

z r ω

• rientation (vert/horiz)

height over the ground h inter-coil distance r angular frequency ω they can be varied in order to generate multiple measures for each geographical position, and realize data inversion

• C. Fenu

Frequency domain electromagnetic data inversion

Parameters which influence the GCM measures

• h

z r ω

• rientation (vert/horiz)

height over the ground h inter-coil distance r angular frequency ω measurements are taken for h = hk = (k − 1)∆h for k = 1, . . . , 20, ∆h = 0.1 m

• C. Fenu

Frequency domain electromagnetic data inversion

The linear model

In 1980, McNeill (Geonics) described a linear model, based on the response curves in the vertical and horizontal positions, which relates the apparent conductivity to the height over the ground mV (h) = ∞ φV (h + z)σ(z) dz mH(h) = ∞ φH(h + z)σ(z) dz where σ(z) is the real conductivity, φV (z) = 4z (4z2 + 1)3/2 , φH(z) = 2 − 4z (4z2 + 1)1/2 , and z is the ratio between the depth and the inter-coil distance r.

• C. Fenu

Frequency domain electromagnetic data inversion

The linear model

In 1980, McNeill (Geonics) described a linear model, based on the response curves in the vertical and horizontal positions, which relates the apparent conductivity to the height over the ground mV (h) = ∞ φV (h + z)σ(z) dz mH(h) = ∞ φH(h + z)σ(z) dz The linear model is valid for: uniform magnetic permeability µ0 = 4π10−7 H/m, small induction number, moderate conductivity (σ 100 mS/m).

• C. Fenu

Frequency domain electromagnetic data inversion

The nonlinear model [Wait 1982], [Hendrickx et al 2002]

R T Depth (z) Height (h) d1 d2 dn-1 dn z1=h1 z2 z3 zn-1 zn Ground surface hi hm σ1 σ2 σn-1 σn µ1 µ2 µn-1 µn Halfspace r h2

• C. Fenu

Frequency domain electromagnetic data inversion

The nonlinear model

The electric conductivity σ(z) and magnetic permeability µ(z) are assumed to be piecewise constant (i.e., layered soil) σ = σ0, σ1, . . . , σm, µ = µ0, µ1, . . . , µm, in each layer of width dk (σ0, µ0 in free space). The model is derived from Maxwell’s equations, keeping into account the cylindrical symmetry of the problem, (H at the receiver is independent of the vertical rotation of the GCM). Let the characteristic admittance inside the kth layer be Nk(λ) = uk(λ) iµkω , uk(λ) =

• λ2 + iσkµkω,

for k = 0, 1, . . . , m, where λ has no direct physical meaning.

• C. Fenu

Frequency domain electromagnetic data inversion

The nonlinear model

The surface admittance Yk at the top of the k-th layer can be

• btained by the following recursion

Yk(λ) = Nk(λ)Yk+1(λ) + Nk(λ) tanh(dkuk(λ)) Nk(λ) + Yk+1(λ) tanh(dkuk(λ)), setting Ym(λ) = Nm(λ) at the lowest layer and letting k = m − 1, m − 2, . . . , 1. Let the reflection factor be R0(λ) = N0(λ) − Y1(λ) N0(λ) + Y1(λ).

• C. Fenu

Frequency domain electromagnetic data inversion

The nonlinear model

Then, the (vert/horiz) apparent conductivity is given by mV (σ, µ; h) = 4r µ0ωH0

• −λe−2hλ Im(R0(λ))
• (r)

mH(σ, µ; h) = 4 µ0ωH1

• −e−2hλ Im(R0(λ))
• (r)

where h is the height above the ground. We denote by Hν[f ](r) = ∞ f (λ)Jν(rλ)λ dλ, ν = 1, 2, the Hankel transform, while J0, J1 are first kind Bessel functions.

• C. Fenu

Frequency domain electromagnetic data inversion

Inversion of the nonlinear model

We consider the residual vector r(σ) = b − m(σ), as a function of the conductivities σ = (σ1, . . . , σn)T, where ri(σ) =

• bV

i − mV (σ, hi),

i = 1, . . . , m, bH

i−m − mH(σ, hi−m),

i = m + 1, . . . , 2m, bV

i , bH i

(i = 1, . . . , m) are the V/H measurements and mV /H(σ, h) are the V/H apparent conductivities at height h. In our analysis, we let the magnetic permeability take the same value µ0 in the n layers. This assumption is approximately met if the ground does not contain ferromagnetic materials.

• C. Fenu

Frequency domain electromagnetic data inversion

Inversion by Newton’s method

The problem of data inversion consists of computing the conductivity σi of each layer which determine a given data set b ∈ R2m, i.e., min

σ∈Rn f (σ),

with f (σ) = 1 2r(σ)2, Newton’s method consists of the iteration σk+1 = σk + sk, where sk solves the linear system f′′(σk)s = −f′(σk). While f′(σ) is available, computing f′′(σ) is hard.

• C. Fenu

Frequency domain electromagnetic data inversion

Inversion by linearization

The Gauss–Newton method is based on the least squares minimization of a linear approximation of r(σk+1) = r(σk + s) min

s∈Rn r(σk) + Jks,

where Jk = J(σk) =

• ∂ri(σk)

∂σj

• , or some approximation.
• C. Fenu

Frequency domain electromagnetic data inversion

Armijo–Goldstein principle

The iteration becomes σk+1 = σk + αksk = σk − αkJ†

k r(σk),

where J†

k is the Moore–Penrose pseudoinverse of Jk and αk is a

damping parameter determined by Armijo–Goldstein principle, that is, αk as the largest number in the sequence 2−i, i = 0, 1, . . . , for which the following inequality holds r(σk)2 − r(σk + αksk)2 ≥ 1 2αkJksk2. This choice of αk ensures convergence of the method, provided that σk is not a critical point and allows us to obtain positive solutions.

• C. Fenu

Frequency domain electromagnetic data inversion

Gauss–Newton vs. Newton

Newton f′′(σk)s = −f′(σk) Gauss–Newton JT

k Jks = −JT k r(σk)

We have f′(σ)= J(σ)Tr(σ) f′′(σ)= J(σ)TJ(σ)+

2m

• i=1

ri(σ)Hi(σ) where [Hi(σ)] is the Hessian of the ith residual ri(σ). The additional term is neglectable if each ri(σ) is either small (compatible problem) or mildly nonlinear at σk.

• C. Fenu

Frequency domain electromagnetic data inversion

Computation of the Jacobian matrix - 1

To obtain derivatives Y ′

kj = ∂Yk ∂σj (k, j = 1, . . . , n) we start from

Y ′

nn =

1 2un , Y ′

nj = 0,

j = 1, . . . , n − 1, and proceeding recursively for k = n − 1, n − 2, . . . , 1 by Y ′

kj= N2 kbkY ′ k+1,j,

j = n, n − 1, . . . , k + 1, Y ′

kk= ak

2uk + bk 2

• N2

kdk − Yk+1

• dkYk+1 +

1 iµkω

• ,

Y ′

kj= 0,

j = k − 1, k − 2, . . . , 1. All functions of λ are complex valued, and

ak = Yk+1 + Nk tanh(dkuk) Nk + Yk+1 tanh(dkuk), bk = 1 [Nk + Yk+1 tanh(dkuk)]2 cosh2(dkuk).

• C. Fenu

Frequency domain electromagnetic data inversion

Computation of the Jacobian matrix - 1

To obtain derivatives Y ′

kj = ∂Yk ∂σj (k, j = 1, . . . , n) we start from

Y ′

nn =

1 2un , Y ′

nj = 0,

j = 1, . . . , n − 1, and proceeding recursively for k = n − 1, n − 2, . . . , 1 by Y ′

kj= N2 kbkY ′ k+1,j,

j = n, n − 1, . . . , k + 1, Y ′

kk= ak

2uk + bk 2

• N2

kdk − Yk+1

• dkYk+1 +

1 iµkω

• ,

Y ′

kj= 0,

j = k − 1, k − 2, . . . , 1. Numerical implementation needs care: cancellation may occur in some formulae arrangement; NaNs appear for large λ’s, because of cosh(√z) with z ∈ C.

• C. Fenu

Frequency domain electromagnetic data inversion

Computation of the Jacobian matrix - 2

Finally, ∂mV (h) ∂σj = 4r µ0ωH0

• λe−2hλ Im

∂R0(λ) ∂σj

• (r),

∂mH(h) ∂σj = 4 µ0ωH1

• e−2hλ Im

∂R0(λ) ∂σj

• (r),

where ∂R0(λ) ∂σj = ∂ ∂σj N0(λ) − Y1(λ) N0(λ) + Y1(λ) = −2iµ0ωλ (λ + iµ0ωY1(λ))2 · ∂Y1 ∂σj . The complexity is lower than using the finite differences approach.

• C. Fenu

Frequency domain electromagnetic data inversion

Approximation of the Jacobian matrix

Finite difference approximation ∂ri σj (σ) = ri(σ + δj) − ri(σ) δ , where δj = δej = (0, . . . , 0, δ, 0, . . . , 0)T, for i = 1, . . . , 2m and j = 1, . . . , n; 2m is the size of the data set and n the number of layers. This is the approach used in the existing package for nonlinear inversion (em38nonlin). Computing the Jacobian is about 5 times faster.

• C. Fenu

Frequency domain electromagnetic data inversion

Update of the Jacobian matrix

Broyden update Starting with J0 = J(σ0), we set Jk = Jk−1 + (yk − Jk−1sk)sT

k

sT

k sk

, where sk = σk − σk−1 and yk = r(σk) − r(σk−1). This formula makes the linearization r(σk) + Jk(σ − σk) exact in σk−1, and guarantees the least change in Jk − Jk−1F. The update makes the computation 5 times faster.

• C. Fenu

Frequency domain electromagnetic data inversion

The problem is ill-conditioned, regularization is needed.

Let J(σk) = UΣV T, σk = 100*rand(n,1), k = 1, . . . , 1000.

5 10 15 20 10

−20

10

−15

10

−10

10

−5

10 5 10 15 20 25 30 35 40 10

−20

10

−15

10

−10

10

−5

10 n=10 n=20 n=30 n=40

Left-side: average singular values and errors (n = 20). Right-side: average singular values for n = 10, 20, 30, 40. The conditioning increases with dimension and does not change much during iteration.

• C. Fenu

Frequency domain electromagnetic data inversion

Nonlinear regularization method: TSVD approach

We replace the ill-conditioned Jacobian J = Jk with the best rank ℓ approximation according to the Euclidean norm Aℓ = min

rank(A)=ℓ J − A2.

Let J = UΓV T and p = rank(J). Then, the regularized Gauss-Newton step vector s(ℓ) can be expressed as s(ℓ) = −A†

ℓr(σk) = − ℓ

• i=1

uT

i r(σk)

γi vi, where ℓ = 1, . . . , p is the regularization parameter.

• C. Fenu

Frequency domain electromagnetic data inversion

Using a regularization matrix: TGSVD

Let us introduce a regularization matrix L ∈ Rt×n (t ≤ n). Common choices are L = D1, D2. We call minimal L-norm solution the vector which solves min

s∈S Ls,

S = {s ∈ Rn | JTJs = −JTr(σk)}, under the assumption N(J) ∩ N(L) = {0}. By the generalized singular value decomposition (GSVD) of (J, L) J = UΣJZ −1, L = V ΣLZ −1, it is possible to define the truncated GSVD (TGSVD) solution sℓ, where ℓ = 0, 1, . . . , p is the regularization parameter and p depends upon p and t.

• C. Fenu

Frequency domain electromagnetic data inversion

Inversion algorithm

We apply the regularized damped Gauss–Newton method σ(ℓ)

k+1 = σ(ℓ) k

+ αks(ℓ)

k

= σ(ℓ)

k

− αkA†

ℓ r(σ(ℓ) k ),

with ℓ fixed, αk determined by Armijo–Goldstein/positivity. We iterate until σ(ℓ)

k

− σ(ℓ)

k−1 < τσ(ℓ) k

• r

k > 100

• r

αk < 10−5, for a given τ. The solution at convergence is denoted by σ(ℓ). The choice of the regularization parameter ℓ is crucial.

• C. Fenu

Frequency domain electromagnetic data inversion

Choosing ℓ: discrepancy

Not all the methods for determining a regularization parameter available for linear problems can be adapted to the nonlinear case. The optimal choice would be the discrepancy principle b − m(σ(ℓdiscrepancy)) ≤ κe, κ > 1, but it can seldom be applied to EMI techniques because in applications the noise on the data is not necessarily equally distributed, an accurate estimate of e is often unknown.

• C. Fenu

Frequency domain electromagnetic data inversion

Choosing ℓ: L-curve

Among the heuristic methods, which do not require an estimate of e, the L-curve principle [Hansen, O’Leary 1993] can be adapted quite naturally to the nonlinear case. It chooses the value of ℓ which identifies the corner of the curve connecting the points

• log r(σ(ℓ)), log Lσ(ℓ)
• ,

ℓ = 1, . . . , p. The curve is L-shaped in many discrete ill-posed problems. To detect the corner of the L-curve we used two algorithms: the L-corner method [Hansen, Jensen, Rodriguez 2007]; the ResReg method [Regi´ nska 1996], [Reichel, Rodriguez 2013].

• C. Fenu

Frequency domain electromagnetic data inversion

Experimental results with synthetic data

exact L−corner discrepancy

m(σk) - b

2 10-2 10-3

D2 σk 2

100 10-1 0.0 0.5 1.0 1.5 2.0 2.5 Depth (m) 1.0 0.8 0.6 0.4 0.2 Conductivity (S/m)

piecewise linear conductivity distribution 10+10 measurements, 40 layers, noise 10−3, L = D1 ℓoptimal = 4, ℓL-corner = 3, ℓdiscrepancy = 2

• C. Fenu

Frequency domain electromagnetic data inversion

Experimental results with synthetic data

exact L−corner discrepancy

m(σk) - b

2 10-2

D2 σk 2

10-1 0.0 0.5 1.0 1.5 2.0 2.5 Depth (m) 1.0 0.8 0.6 0.4 0.2 Conductivity (S/m)

smooth conductivity distribution 10+10 measurements, 40 layers, noise 10−2, L = D2 ℓoptimal = 2, ℓL-corner = 2, ℓdiscrepancy = 1

• C. Fenu

Frequency domain electromagnetic data inversion

Experimental results with synthetic data

1.0 0.8 0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0 0.0 1.0 2.0 0.0 1.0 2.0 0.0 1.0 2.0 Depth (m) Depth (m) Depth (m) Conductivity (S/m) Conductivity (S/m)

identification of thin conductive layers fξ(z) = 1 for z ∈ [0.5, 0.5 + ξ] ξ = 1.0, 0.5, 0.2 10+10 measurements, 40 layers, noise 10−2, L = D1

• C. Fenu

Frequency domain electromagnetic data inversion

Experiments on field data: the Airport data set

• C. Fenu

Frequency domain electromagnetic data inversion

Experiments on field data: the Airport data set

0.5 1 1.5 2 50 100 150 200 250 300 350 Height (m) Apparent conductivity (mS/m)

20 measurements for each hi = 0, 0.1, . . . , 1.9 m error bars are standard deviations (multiplied by 10)

• C. Fenu

Frequency domain electromagnetic data inversion

Experiments on field data: the Airport data set

Depth (m) Depth (m)

1 0.4 0.8 1.2 1.6 2 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 5 10 15 20

Distance (m) Conductivity (S/m) Conductivity (S/m)

(a) (b)

Electrical resistivity tomography (ERT)

Iris Instruments Syscal Pro Switch 48 resistivity meter 48 electrodes set up with an inter-electrode spacing of 0.5 m There is a conductivity maximum at depth ∼ 1.6 m

EM data measured at x = 16 m

GF Instruments CMD-1 conductivity meter

• perating frequency 10 kHz, coil separation 0.98 m
• C. Fenu

Frequency domain electromagnetic data inversion

Experiments on field data: the Airport data set

2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 0.0 1.0 2.0 0.0 1.0 2.0 0.0 1.0 2.0 Depth (m) Depth (m) Depth (m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m)

σ2 σ3 σ4 σ5 σ6 σ7

Gauss–Newton solutions - ERT conductivity profile

20+20 measurements, 40 layers, L = I ℓL-corner = 5, ℓResReg = 4

• C. Fenu

Frequency domain electromagnetic data inversion

Experiments on field data: the Airport data set

2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 0.0 1.0 2.0 0.0 1.0 2.0 0.0 1.0 2.0 Depth (m) Depth (m) Depth (m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m)

σ1 σ2 σ3 σ4 σ5 σ6

Gauss–Newton solutions - ERT conductivity profile

20+20 measurements, 40 layers, L = D1 ℓL-corner = 2, ℓResReg = 1

• C. Fenu

Frequency domain electromagnetic data inversion

Experiments on field data: the Airport data set

2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 0.0 1.0 2.0 0.0 1.0 2.0 0.0 1.0 2.0 Depth (m) Depth (m) Depth (m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m)

σ1 σ2 σ3 σ4 σ5 σ6 σ1 σ2 σ3 σ4 σ5 σ6

Gauss–Newton solutions - ERT conductivity profile

20+20 measurements, 40 layers, L = D2 ℓL-corner = 2, ℓResReg = 1

• C. Fenu

Frequency domain electromagnetic data inversion

Experiments on field data: the Airport data set

2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 2.0 1.5 1.0 0.5 0.0 0.0 1.0 2.0 0.0 1.0 2.0 0.0 1.0 2.0 Depth (m) Depth (m) Depth (m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m) Conductivity (S/m)

σ1 σ2 σ3 σ1 σ2 σ3

m = 10 m = 10 m = 10 m = 5 m = 5 m = 5

Gauss–Newton solutions - ERT conductivity profile

10+10 and 5+5 measurements, 40 layers, L = D2, ℓL-corner = 2

• C. Fenu

Frequency domain electromagnetic data inversion

Conclusions

Numerical experiments showed that: a small number of solutions σℓ has to be computed the L-corner method produces acceptable solutions meaningful results can be obtained with m = 5 Broyden update is extremely effective computation is very fast even on a notebook L = D1, D2 better than L = I In the future: work with multifrequency/multi intercoil distance GCM measurements consider also the variation of the magnetic permeability consider also the “in-phase component” of the signal

• C. Fenu

Frequency domain electromagnetic data inversion

• C. Fenu

Frequency domain electromagnetic data inversion