# Solution of the Inverse Problem of Subsurface Electric Exploration for Horizontally Stratified Medium

Authors: A. L. Karchevsky; Zh. O. Oralbekova; K. T. Iskakov
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

In this paper, we offer the generalization of the known technique of the construction of the gradient of the residual functional based on the statement of the conjugate problem for the case when the unknown function is complex valued. The notion of the reference frequency of the medium is introduced. Knowing the value of the reference frequency lets us judge the possibility of simultaneous definition of the dielectric permittivity and conductivity.


---

## Body

## 1. Introduction

The problem of definition of the dielectric permittivity and conductivity is the actual geophysical problem. Theoretical research of inverse problems of geoelectrics and some numerical examples of definition of these functions can be found in [1].There are many works devoted to the problems of reconstruction of conductivityσ and dielectric permittivity ε (or complex conductivity σ-iωε, or complex dielectric permittivity ε+iσ/ω), e.g., in [2–19]). According to the rules, the inverse problem is solved in frequency domain; that is, the external current source is harmonic. Despite the quite obvious idea that we can search for one complex function ε+iω0σ, the authors assume numerical algorithms for reconstruction either of two real functions σ and ε or of real and imaginary parts of correspondent complex function, considering that one of the functions σ or ε is known. It leads to the fact that the resulting formulas are too complicated and inconvenient for analysis and implement on a computer. Apparently, another reason such numerical algorithms are offered is that when solving the inverse problem of finding a complex function w by minimizing the functional J[w], we go through the following steps: first, the functional J:ℂ→ℝ and, secondly, by the definition of the gradient, increment of the functional should be presented in the form δJ[w]=〈J′[w],δw〉+o(∥δw∥); that is, there should be a scalar product for two complex values 〈·,·〉:ℂ×ℂ→ℝ.In this paper the authors propose a numerical method of simultaneous definition of the conductivityσ and dielectric permittivity ε considering the example of inverse problem of subsurface radiolocation. Instead of two unknown real functions, we consider one complex. The authors generalize the known technique of finding of the gradient of the residual functional, using the statement of the conjugate problem for the case when the unknown function is complex valued. The notion of the reference frequency of the medium is introduced, which helps to understand where these two functions can be determined simultaneously. Some test reconstructions for simulated data are offered.


## 2. Statement of the Inverse Problem

We consider the media—Nl-layered structure with interfaces zk (k=0,Nl¯), z0=0; m-layer is the interval [zm-1,zm], the last Nl+1 (underlying) layer is the half space [zNl,∞), and the air is the half space (-∞,0].Electromagnetic properties of each layer are defined by the permittivityε0ε, the conductivity σ, and by magnetic permeability μ0μ, ε0=8,854·10-12 F/m and μ0=4π·10-7 H/m, in most cases the relative permittivity ε belongs to the interval [1;80] and relative magnetic permeability μ=1. Since the medium is horizontally stratified, then ε and σ are piecewise-constant functions of the variable z (z∈(-∞,∞)).Let the source of external current be a cable disposed on the heightz* parallel to the y-axis.For the componentE2(t,x,z) from the Maxwell equations finally, we can obtain the differential equation of the second order (see, e.g., [1]):
(1)ε∂2E2∂t2+σ∂E2∂t=1μ0(∂2E2∂x2+∂2E2∂z2).
The Fourier transform with respect to the horizontal variable x and the time variable t gives the following equation:
(2)uzz-(λ2-ω2μ0ε0ε+iωμ0σ)u=0.
At points of discontinuity of the medium, we assume that the gluing conditions are as follows:
(3)[u]zk=0,[uz]zk=0,k=0,Nl¯.
Source concentrated at the point z*=0, which is equivalent to the gluing conditions at this point
(4)[u]z*=0,[uz]z*=-f(ω)μ0.
We assume that we have the conditions of damping in infinity
(5)u⟶0(z⟶±∞),
and, relative to the solution of the direct problem (2)–(5), the additional information is given
(6)u|z=0=g(ω,λ).Hereλ and ω are the Fourier parameters with respect to variables x and t, respectively, [·]z is the notation for gluing, that is, [w]z=w(z+0)-w(z-0), and the bar over the complex value will denote the complex conjunction.The inverse problem is to find the piecewise-constant functionsε and σ, if for solution of the direct problem (2)–(5), the additional information (6) is known.Introduce the notation(7)k2=λ2-ϰ,ϰ=ω2μ0ε0ε-iωμ0σ.Fix some values of the angular frequencyω0. In the inverse problem, we will recover the complex value ϰ0 (ϰ0=ω02μ0ε0ε-iω0μ0σ) which is a piecewise-constant function since the functions ε and σ are piecewise constant.As is easy to see, if we findϰ0 as a solution of the inverse problems (2)–(6), we immediately recover the functions ε and σ as follows:
(8)ε=Re[ϰ0]ω02μ0ε0,σ=-Im[ϰ0]ω0μ0.Inverse problem (2)–(6) may be solved by minimizing the residual functionals as follows:
(9)J[ϰ0]=∑ω‍hω|u(0,ω,λ0)-g(ω,λ0)|2,
(here hω are certain weight multipliers). Paying attention to that in the functional (9), we fix a value of the spatial frequency λ0 and prepare additional information (6) for different values of the angular frequency ω.For minimization, we will use the gradient method, since the rate of convergence of such method is higher than that of the method that uses only the values of the functional; therefore, we need to get the gradient of the residual functional (9), which in turn requires the definition of the scalar multiplication of two complex numbers.For two complex numbersw1 and w2, we introduce the following relation:
(10)〈w1,w2〉=12(w1w¯2+w¯1w2).
The relation (10) has all the properties of scalar multiplication. The proof of this is based on a geometric interpretation of complex numbers. Let w1=(α1,β1) and w2=(α2,β2); then, 〈w1,w2〉=α1α2+β1β2.For (10) note useful relation
(11)〈w1,w2〉=Re[w1w¯2]=Re[w¯1w2].


## 3. Gradient of the Residual Functional

First we note that(12)J[ϰ0]=J[ϰ01,…,ϰ0m,…,ϰ0Nl],
where ϰ0m are the values of a piecewise constant function ϰ0 in the segment [zm-1,zm].We obtain the expression of the gradient of the residual functional (9) by the statement of the conjugate problem
(13)ψzz-(λ2-ϰ)ψ=0,[ψ]zk=0,[ψz]zk=0,k=1,Nl¯,[ψ]z*=0,[ψz]z*=-2[u(0,ω0,λ)¯-g(ω0)¯],ψ⟶0(z⟶±∞).Let the valueϰ be incremented δϰ; then, the function u will get increment δu, which satisfies the following problem:
(14)δuzz-k2δu+uδϰ=0,[δu]zk=0,[δuz]zk=0,k=-1,Nl¯,u⟶0(z⟶±∞)
(here for, simplicity, we set z-1=z*).In this case, the increment of the residual functional up the second order can be obtained as follows:(15)δJ[ϰ0]=Re[∑λ‍2hλ(u¯-g¯)δu]=-Re(∑λ‍hλ[ψ]z0δu).Given statements of problems (14), (13), and the equality
(16)∫-∞∞‍ddz(δuψz)dz=-δu|z=z0[ψz]z0,
we get
(17)δJ[ϰ0]=Re[∑λ‍hλ∫z0zNl‍uψδϰdx].
Since in each layer [zm-1,zm], the function ϰ is constant and takes the value ϰm, we can write
(18)δJ[ϰ0]=∑m=1Nl‍Re[∑ω‍hω∫zm-1zm‍uψdz·δϰm].Let(19)Nm=∫zm-1zm‍uψdz.
Then
(20)Re[Nmδϰm]=Re[Nm]·ω2μ0ε0δεm+Im[Nm]·ωμ0δσm=Re[Nm]·ω2ω02·ω02μ0ε0δεm+Im[Nm]·ωω0·ω0μ0δσm=Re[N^mδϰ0m],
where
(21)N^m=Re[Nm]ω2ω02+iIm[Nm]ωω0.Therefore, the gradient of the residual functionalJ[ϰ0] will have the form
(22)J′[ϰ0]=(j¯1,…,j¯m,…,j¯Nl),jm=∑ω‍hωN^m.


## 4. Analytic Formulas for the Key Expressions

In order to solve the direct problem (2)–(5) and conjugate problem (13), we will reduce the differential equation of the second order to the Riccati equation. This method was successfully used, for example, in [20–23].For solving differential equation (2), we introduce the function s as follows:
(23)uz(z)=s(z)u(z),
which will satisfy the differential Riccati equation
(24)s′+s2=k2.
If we calculate from the right to the left, then the solution of (24) in each segment [zm-1,zm] will be as follows:
(25)s(z)=km(sm+km)e2km(z-zm)+(sm-km)(sm+km)e2km(z-zm)-(sm-km).
If we calculate from the left to the right, then
(26)s(z)=km(sm-1+km)+(sm-1-km)e2km(zm-1-z)(sm-1+km)-(sm-1-km)e2km(zm-1-z)
(here sm=s|z=zm, km is the value of the piecewise constant function k in the segment [zm-1,zm] and Re[km]≥0).For solving Riccati equation (24), we will do our recurrent calculations from the layer up to the layer moving to the point z*, where the source is disposed.The condition of damping in infinity (5) enables
(27)s(z)=-kNl+1,z∈[zNl,∞).
Due to conditions (6), we obtain the gluing conditions
(28)[s]zk=0.
Therefore, we may set
(29)sNl=-kNl+1
and begin the recurrent calculation of sm from the right to the left by formula
(30)sm-1=km(sm+km)e-2kmhm+(sm-km)(sm+km)e-2kmhm-(sm-km),hm=zm-zm-1.
Thus we get s*=s|z=z*+0.Analogously, taking into account damping in minus infinity (5) we get
(31)s(z)=k0,z∈(∞,z*].
Since the source is disposed in the half space (-∞,0], we may immediately take
(32)s|z=z*-0=k0.Gluing condition (28) in the point z=z* allows us to determine
(33)u|z=z*≡u*=-f(ω)μ0s*-k0.Integrating (23) in the interval [z*,0], we obtain
(34)u|z=0≡u0=u*2k0e-k0z*(s0+k0)e-2k0z*-(s0-k0).Further, integrating in each interval[zm-1,zm] (23), we derive the solution of the problem (2)–(5) in this interval
(35)u(z)=um-1l(z),l(z)=ekm(zm-1-z)(sm+km)e2km(z-zm)-(sm-km)(sm+km)e-2kmhm-(sm-km).The conjugate problem (13) is similar to direct problem (2)–(5); therefore, it can be solved similarly. Moreover, since the damping conditions in infinity are the same, the solution of the Riccati equation introduced for the function ψ will coincide with s(z) for all z∈[0,∞). Therefore, we have
(36)ψ|z=0≡ψ0=-2[u(0,ω,λ)¯-g(ω,λ)¯]s0-k0
and in each interval [zm-1,zm]. 						(37)ψ(z)=ψm-1l(z).Taking into account (35) and (37), we can obtain the formulas for the gradient components of Jω[ϰ0] (α=λ,ω); that is, calculate the integral
(38)∫zm-1zm‍uψdz=um-1ψm-1Rm[(sm+km)e-2kmhm-(sm-km)]2,
where
(39)Rm=1-e-2kmhm2km×[(sm+km)2e-2kmhm+(sm-km)2]-2hm(sm+km)(sm-km)e-2kmhm.


## 5. Numerical Experiment

### 5.1. Reference Frequency of the Mediumω0

First we note, that (see, e.g., formulas (25) and (35))
(40)u=u(k),
where
(41)k2=λ2-(ω2μ0ε0ε-iωμ0σ),k=[λ2-ω2μ0ε0ε]2+[ωμ0σ]24eiβ/2,β={arctgωμ0σλ2-ω2μ0ε0ε,λ2>ω2μ0ε0ε,π2,λ2=ω2μ0ε0ε,π+arctgωμ0σλ2-ω2μ0ε0ε,λ2<ω2μ0ε0ε.
That is, the variation of the solution of direct problem (2)–(5) depends directly on how the function k depends on variations of the functions ε and σ in the layers.Evidently, the greatest influence on change in the valuek is rendered by variations ε and σ, when
(42)ω2μ0ε0ε≈ωμ0σ,λ2~ω2μ0ε0ε.We putω2μ0ε0ε=ωμ0σ, whence we get
(43)ω0=σε0ε.
This means that we know for certain medium the mean values of dielectric permittivity and conductivity. Then we derive and fix the value of reference frequency of the medium (43).


### 5.2. Dependence of the Properties of the FunctionalJ[ϰ] from the Value ω

For numerical experiments below we choose model shown in Table1.Table 1
Model of the medium 1.
Number of the layer
1
2
3
4
5
6
ε
18.5
22.8
18.4
19.2
28.3
30.0
σ
0.017
0.024
0.016
0.017
0.022
0.024
z
k
0.11
0.21
0.39
0.58
0.82We set the mean values atε=20 and σ=0,020, and then the reference angular frequency is ω0=1.12·108. The capacity of the skin-layer, from here, is
(44)hs=2ω0μ0σ≈0.84(m).Fixing the valuesεm and σm (m=2,5¯), we will change the values of dielectric permittivity and conductivity of the first layer in the segment
(45)ε1∈[1.0,40.0],σ1∈[0.01,0.03],
for different values of ω=rω0, where r={1/10,1/7,1/5,1/2,1,2,5,7,10} and we will observe the behavior of the value
(46)Δω=|u(0,ω,λ0)-g(ω,λ0)|2,
which is one of the components of the residual functional (9). Results of numerical experiment are presented in Figure 1.Figure 1
The behavior ofΔω for various values ω.We note that first, ifω decreases with respect to the reference frequency of the medium ω0, then the residual functional loses sensitivity when ε1 varies. If ω increases, then the functional loses the sensitivity when σ1 varies. Secondly, the more the difference between ω and ω0 the less the sensitivity of the residual functional J[ϰ] to the change of desired unknown values, since it becomes more standing.Ifεm and σm (m=2,5¯) vary, the behavior of Δω will be similar; however, if the layer is lower, so is the sensitivity of the functional.Thus, the assumption that the value of the angular frequencyω must by such that
(47)ω2μ0ε0ε≈ωμ0σ
is verified (see (42)).In addition, we clarify the known geophysical condition of quasistationarity of the electromagnetic field:ωε0ε≪σ. The value of notation “≪” is lax. Proceeding from the numerical experiment, we may account the quasistationary approximation to be suitable if ωε0ε≤σ/10.


### 5.3. Dependence of the Properties of the Residual Functional from the Valueλ

For the model 1 (see Table1), the numerical experiment similar to those in the previous section was developed. The values of the dielectric permittivity and conductivity of the first layer were changing (in above mentioned segments) and the values of the residual functional J[ϰ0] were derived. The result is shown in Figure 2.Figure 2
The behavior of the residual functionalJ[ϰ0] for various values of λ. 450 values of ω were taken from the segment [ω0/10;10ω0] with constant step.The tendency is seen as the greaterλ the less is the sensitivity of residual functional to variations of the unknown parameters. The assumption that
(48)λ2~ω2μ0ε0ε
is verified.For model 1 we obtainω2μ0ε0ε≈2.84. It is seen that with increasing λ at first the residual functional J[ϰ0] loses the sensitivity to the variations of σ1 and then to variations of ε1. When λ2≥2ω2μ0ε0ε, the sensitivity of the residual functional rapidly decreases.


## 5.1. Reference Frequency of the Mediumω0

First we note, that (see, e.g., formulas (25) and (35))
(40)u=u(k),
where
(41)k2=λ2-(ω2μ0ε0ε-iωμ0σ),k=[λ2-ω2μ0ε0ε]2+[ωμ0σ]24eiβ/2,β={arctgωμ0σλ2-ω2μ0ε0ε,λ2>ω2μ0ε0ε,π2,λ2=ω2μ0ε0ε,π+arctgωμ0σλ2-ω2μ0ε0ε,λ2<ω2μ0ε0ε.
That is, the variation of the solution of direct problem (2)–(5) depends directly on how the function k depends on variations of the functions ε and σ in the layers.Evidently, the greatest influence on change in the valuek is rendered by variations ε and σ, when
(42)ω2μ0ε0ε≈ωμ0σ,λ2~ω2μ0ε0ε.We putω2μ0ε0ε=ωμ0σ, whence we get
(43)ω0=σε0ε.
This means that we know for certain medium the mean values of dielectric permittivity and conductivity. Then we derive and fix the value of reference frequency of the medium (43).


## 5.2. Dependence of the Properties of the FunctionalJ[ϰ] from the Value ω

For numerical experiments below we choose model shown in Table1.Table 1
Model of the medium 1.
Number of the layer
1
2
3
4
5
6
ε
18.5
22.8
18.4
19.2
28.3
30.0
σ
0.017
0.024
0.016
0.017
0.022
0.024
z
k
0.11
0.21
0.39
0.58
0.82We set the mean values atε=20 and σ=0,020, and then the reference angular frequency is ω0=1.12·108. The capacity of the skin-layer, from here, is
(44)hs=2ω0μ0σ≈0.84(m).Fixing the valuesεm and σm (m=2,5¯), we will change the values of dielectric permittivity and conductivity of the first layer in the segment
(45)ε1∈[1.0,40.0],σ1∈[0.01,0.03],
for different values of ω=rω0, where r={1/10,1/7,1/5,1/2,1,2,5,7,10} and we will observe the behavior of the value
(46)Δω=|u(0,ω,λ0)-g(ω,λ0)|2,
which is one of the components of the residual functional (9). Results of numerical experiment are presented in Figure 1.Figure 1
The behavior ofΔω for various values ω.We note that first, ifω decreases with respect to the reference frequency of the medium ω0, then the residual functional loses sensitivity when ε1 varies. If ω increases, then the functional loses the sensitivity when σ1 varies. Secondly, the more the difference between ω and ω0 the less the sensitivity of the residual functional J[ϰ] to the change of desired unknown values, since it becomes more standing.Ifεm and σm (m=2,5¯) vary, the behavior of Δω will be similar; however, if the layer is lower, so is the sensitivity of the functional.Thus, the assumption that the value of the angular frequencyω must by such that
(47)ω2μ0ε0ε≈ωμ0σ
is verified (see (42)).In addition, we clarify the known geophysical condition of quasistationarity of the electromagnetic field:ωε0ε≪σ. The value of notation “≪” is lax. Proceeding from the numerical experiment, we may account the quasistationary approximation to be suitable if ωε0ε≤σ/10.


## 5.3. Dependence of the Properties of the Residual Functional from the Valueλ

For the model 1 (see Table1), the numerical experiment similar to those in the previous section was developed. The values of the dielectric permittivity and conductivity of the first layer were changing (in above mentioned segments) and the values of the residual functional J[ϰ0] were derived. The result is shown in Figure 2.Figure 2
The behavior of the residual functionalJ[ϰ0] for various values of λ. 450 values of ω were taken from the segment [ω0/10;10ω0] with constant step.The tendency is seen as the greaterλ the less is the sensitivity of residual functional to variations of the unknown parameters. The assumption that
(48)λ2~ω2μ0ε0ε
is verified.For model 1 we obtainω2μ0ε0ε≈2.84. It is seen that with increasing λ at first the residual functional J[ϰ0] loses the sensitivity to the variations of σ1 and then to variations of ε1. When λ2≥2ω2μ0ε0ε, the sensitivity of the residual functional rapidly decreases.


## 6. Numerical Examples of Solution of the Inverse Problem

In order to test the operability of the proposed numerical algorithm, we carry out a number of reconstruction of electromagnetic properties of the medium, using simulated data. In order to obtain the additional information (6) we, first, solve the direct problem (2)–(5). Then we add the random value with this form
(49)g~(ω,λ0)=g(ω,λ0)(1+P100ξ),
where ξ is a random value from the unit circle and P is the percent of the introduced error. In Figure 3 we see the example of using the additional information with this error.Figure 3
The example of additional informationg(ω,λ) given with the error P=20% (ω∈[ω0/50,ω0/10], λ=5·10-1); the exact values of real and imaginary parts of g(ω,λ) are shown by the smooth line.In addition to the model of the medium 1 in Table1 were selected four more models (see Tables 2, 3, 4, and 5).Table 2
Model of the medium 2.
Number of the layer
1
2
3
4
5
6
ε
20.5
22.7
17.9
22.8
21.0
30.0
σ
0.021
0.023
0.020
0.021
0.020
0.025
z
k
0.11
0.22
0.39
0.62
0.82Table 3
Model of the medium 3.
Number of the layer
1
2
3
4
5
6
ε
2.05
2.27
2.05
2.08
1.93
2.50
σ
0.0021
0.0023
0.0020
0.0021
0.0020
0.0025
z
k
0.48
0.88
1.38
2.00
2.50Table 4
Model of the medium 4.
Number of the layer
1
2
3
4
5
6
7
8
9
10
11
ε
20.2
21.3
22.1
20.4
18.4
16.2
17.8
18.8
22.2
23.3
25.0
σ
0.0018
0.0021
0.0020
0.0021
0.0020
0.0018
0.0019
0.0021
0.0022
0.0023
0.0025
z
k
0.72
1.53
2.91
3.65
4.45
5.05
6.00
6.65
7.23
8.20Table 5
Model of the medium 5.
Number of the layer
1
2
3
4
5
6
7
8
ε
20.2
19.3
20.2
20.9
21.1
19.2
23.4
2.50
σ
0.021
0.020
0.021
0.019
0.018
0.019
0.022
0.025
z
k
0.11
0.22
0.33
0.44
0.59
0.72
0.83For each model of the medium the coordinate of the boundary of the last layerzNl coincides with the value of capacity of the skin-layer hs.The parameters used in the construction of the residual functionalJ[ϰ0] are collected in Table 6.Table 6
ω
0
Frequency interval
Number of frequencies
Model of the medium 1
1.12
·
1
0
8
[
ω
0
/
10,10
ω
0
]
2500
Model of the medium 2
1.12
·
1
0
8
[
ω
0
/
10,10
ω
0
]
2500
Model of the medium 3
1.12
·
1
0
8
[
ω
0
/
10,10
ω
0
]
2500
Model of the medium 4
1.12
·
1
0
7
[
ω
0
/
40,40
ω
0
]
6000
Model of the medium 5
1.12
·
1
0
8
[
ω
0
/
40,40
ω
0
]
6000In each layer, the initial approximations wereεNl and σNl.For minimizing of the residual functionalJ[ϰ0] the conjugate gradient method modified for complex values used
(50)ϰ0n+1=ϰ0n-αnpn,pn=J′[ϰ0n]-βnpn-1,p0=J′[ϰ00],βn=-∥J′[ϰ0n]∥2∥J′[ϰ0n-1]∥2,αn=argminα>0J[ϰ0n-αpn].The results of restoration of piecewise constant functionsε and σ are shown in Figure 4.The examples or recoveryε and σ: (a) model of the medium 1, (b) model of the medium 2, (c) model of the medium 3, (d) model of the medium 4, and (e) model of the medium 5. The exact solution is shown by the continuous line. The recovered solution is shown by the dotted line. For all the examples P=20%.
(a)
(b)
(c)
(d)
(e)Remember that we used the value of the reference angular frequencyω0 such that the condition ω02μ0ε0ε≈ω0μ0σ is satisfied. Recovery experiment shows that the less the value of ωμ0σ, the more difficult restoring of the functions ε and σ is, since small values of k2 in the differential equation (2) render little impact on changing values of the solution of the equation. The consequence of this is the large “flatness” of the residual functional and its low sensitivity to variations of ε and σ. Note: ω0μ0σ=2/hs2.


## 7. Conclusions

In this paper, the authors suggested the generalization of the known technique of constructing the gradient of the residual functional with the use of the statement of the conjugate problem, when the unknown function is complex.The numerical examples had shown that the conjugate gradient method in the case of complex valued gradient of the residual functional and the known function is applicable and we can find the minimum of the functional.Efficiency rate of the method is confirmed by examples of simultaneous reconstruction of the dielectric permittivity and conductivity on synthetic data with introduced random error.



---

*Source: 432121-2013-11-20.xml*