- Email: [email protected]
3D FDTD anisotropic and dispersive modeling for GPR using rotated staggered grid method
Journal Pre-proof 3D FDTD anisotropic and dispersive modeling for GPR using rotated staggered grid method Yongxu Lu, Suping Peng, Xiaoqin Cui, Dong li, Kang Wang PII:
S0098-3004(18)31218-4
DOI:
https://doi.org/10.1016/j.cageo.2019.104397
Reference:
CAGEO 104397
To appear in:
Computers and Geosciences
Received Date: 25 December 2018 Revised Date:
30 October 2019
Accepted Date: 3 December 2019
Please cite this article as: Lu, Y., Peng, S., Cui, X., li, D., Wang, K., 3D FDTD anisotropic and dispersive modeling for GPR using rotated staggered grid method, Computers and Geosciences (2020), doi: https://doi.org/10.1016/j.cageo.2019.104397. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
1
3D FDTD anisotropic and dispersive modeling for GPR using rotated
2
staggered grid method
3
Yongxu Lua,b,∗ , Suping Penga,b, Xiaoqin Cuib, Dong lia,b, Kang Wanga,b
4
a
5
Beijing 100083, China.
6
b
7
(Beijing), Beijing 100083, China.
College of Geoscience and Surveying Engineering, China University of Mining & Technology (Beijing),
State Key Laboratory of Coal Resources and Safe Mining, China University of Mining & Technology
8
Corresponding author:
9
Yongxu Lu, Ph.D. in Geophysics
10 11
Affiliation: College of Geoscience and Surveying Engineering, China University of Mining & Technology (Beijing).
12
Address: Ding No.11 Xueyuan Road, Haidian District, Beijing 100083 P. R.China
13
E-mail address: [email protected] (Y. Lu).
14
Tel.: +86 1565 071 6135
15
____________________________________
16
Yongxu Lu developed all of the theory, derived all the equations and wrote this paper and all the code.
17
Suping Peng and Xiaoqin Cui participated in developing the basic theory.
18
Dong Li and Kang Wang participated in debugging the code and revising this paper.
19 20 1
21
Abstract
22
The finite-difference time-domain (FDTD) method is widely used in numerical modeling of
23
ground penetrating radar (GPR). The well-known Yee’s method is the earliest and most popular
24
FDTD method. However, using Yee’s method to solve an anisotropic case requires considerable
25
interpolation of electric fields (assuming media to be non-magnetic). In this paper we present a
26
different approach called rotated staggered grid (RSG) method to simulate electromagnetic (EM)
27
wave propagation in anisotropic and dispersive media. RSG-FDTD places electric field
28
components on the same grid system and magnetic field components on the same grid system,
29
however, half a spacing away. The RSG-FDTD method was first used in elastic wave modeling. It
30
requires no interpolations of wave fields and is advantageous for modeling wave propagation in
31
anisotropic media. In this paper, the 3D anisotropic and dispersive EM wave equations are
32
presented, and the RSG-FDTD approximations are derived. The application of the new method is
33
tested using homogeneous media models. To validate the proposed method, results are compared
34
with those of GprMax, and good agreement is achieved. The common offset sections of complex
35
cases with anisotropy and dispersion are also tested and are compared with the corresponding
36
isotropic case.
37 38
Keywords: Finite-difference time-domain (FDTD), Rotated staggered grid (RSG) method, 3D
39
GPR anisotropic modeling, 3D GPR dispersive modeling
40
1. Introduction
41
Ground penetrating radar (GPR) is a well-known geophysical technique that is applied in 2
42
various applications across the fields of geology, archaeology, glaciology, and engineering. In
43
engineering, it is used to analyse asphalt, concrete, wood, and a few other materials
44
(Rodríguez-Abad et al., 2010; Kumlu and Erer, 2019). It transmits a high frequency (100 MHz to
45
4 GHz) electromagnetic (EM) signal into the earth or other target media, and receives the
46
reflections using a set of antennas. By processing and analysing the reflection signals, structures
47
and physical subsurface properties are detected.
48
Numerical modeling techniques can be used to understand the GPR response of certain models
49
or for testing new processing techniques (Giannakis et al., 2019). The finite-difference
50
time-domain (FDTD) method, which is straightforward and easy to be programmed, is well
51
established for EM wave simulation (Diamanti and Giannopoulos, 2009; Millington and Cassidy,
52
2010; Fang and Lin, 2012; Li et al., 2012; Ralchenko et al., 2015). Almost all FDTD methods
53
published are based on Yee’s algorithm (1966), which can be conveniently adopted. The basic idea
54
is using a staggered grid on which spatial derivatives are located between two original grid points.
55
Yee’s algorithm is adopted by many researchers to study high frequency EM waves propagating in
56
2D or 3D isotropic media (Wang and Tripp, 1996, Chen and Huang, 1998; Holliger and Bergman,
57
2002; Giannopoulos, 2005; Irving and Knight, 2006; Warren and Giannopoulos, 2011; Giannakis
58
et al., 2016; Warren et al., 2016; Warren et al., 2019).
59
As we mentioned, several papers discussing forward modeling of GPR focus on isotropic media.
60
However, natural or manmade materials are always anisotropic and/or dispersive. In homogenous
61
isotropic media, wave properties (like wave speed for example) are the same in all directions. In
62
anisotropic media, however, the wave properties are different depending on the propagation
3
63
direction. It is not appropriate to interpret data containing anisotropic information using isotropic
64
theory. EM anisotropy is the variation of electric conductivity, dielectric permittivity or magnetic
65
permeability with azimuth inside the considered material. Rodríguez-Abad et al. (2010) evaluated
66
moisture content in wood using GPR. They proposed analyzing the moisture content after
67
determining the wood anisotropic effects. Bradford et al. (2013) studied the EM wave velocity
68
anisotropy in temperate glaciers. The study of the propagation of EM waves in anisotropic media
69
may help us to learn more about the targets. Beker et al. (1989) and Strikel and Taflove (1990)
70
developed algorithms for anisotropic media based on Yee’s algorithm. However, the anisotropic
71
media adopted in their studies are vertical transverse isotropic (VTI), which means the electric
72
conductivity, dielectric permittivity or magnetic permeability tensors are diagonal. Schneider and
73
Hudson (1993) developed an algorithm for general anisotropic media. Their algorithm is also
74
based on Yee’s algorithm and requires interpolation of electrical and magnetic fields. Moreover,
75
they show results of their algorithm in just a one-dimension case.
76
In practice, several materials, such as water and soils, are dispersive (Teixeira et al., 1998).
77
Single and multi-pole Debye, Lorentz and Drude models are widely used to simulate dispersive
78
media (Bergmann et al., 1998). Several different FDTD based methods for EM wave simulation in
79
the dispersive model have been suggested (Alsunaidi and Al-Jabr, 2009; Giannakis and
80
Giannopoulos, 2014). The recursive convolution algorithms are efficient in handling the
81
dispersive case (Kelley and Luebbers, 1996; Fan and Liu, 2000; Giannakis and Giannopoulos,
82
2014).
83
In this study, we apply a rotated staggered grid (RSG) scheme first reported by Saenger
4
84
(Saenger and Gold, 2000; Saenger et al., 2004) for elastic waves to GPR applications. This grid
85
technique places components of the same wave fields on the same grids. For anisotropic media,
86
interpolation of EM wave fields is not required (Schneider and Hudson, 1993). In our simulation,
87
we place electric fields on a regular grid and magnetic fields on an “half-grid.” We also apply the
88
convolutional perfectly matched layer (C-PML) (Roden and Gedney, 2000; Irving and Knight,
89
2006; Li and Matar, 2009) as an absorbing boundary condition. C-PML, introduced by Roden and
90
Gedney (2000), dose not require splitting of components and has a better performance than a
91
perfectly matched layer (PML). First, we discuss the basic theory of GPR modeling in isotropic,
92
anisotropic and dispersive media, including Maxwell’s equation, and provide a brief introduction
93
of Yee’s FDTD method. Then we introduce the RSG method and derive an updated equation for
94
3D EM wave propagation in anisotropic and dispersive media. Furthermore, we give the
95
numerical stability and dispersion criteria for the RSG method. Finally, we present some examples
96
for testing the proposed GPR modeling method. A few of the results are compared with those
97
obtained with free software GprMax.
98
2. Basic theory of FDTD GPR modeling with anisotropy and dispersive
99
2.1. Maxwell’s equations for isotropic and anisotropic medium
100
GPR uses high frequency EM waves to detect media underground. The basic governing
101
equations for EM phenomena are the well-known Maxwell’s curl equations (Schneider and
102
Hudson, 1993), ∇ × E = −µ
∂H , ∂t
5
(1)
∇ × H = σE + ε
∂E , ∂t
(2)
103
where H and E are the magnetic field intensity and electric field strength vectors, respectively. σ ,
104
ε and µ are the electric conductivity tensor, dielectric permittivity tensor and magnetic
105
permeability tensor, respectively. In this paper we assume that the media are considered as
106
non-magnetic. Thus µ = µ0 , where µ0 = 1.256 × 10−6 H/m in vacuum. σ and ε can be written
107
as follows:
σ xx σ xy σ xz σ = σ yx σ yy σ yz , σ zx σ zy σ zz
ε xx ε = ε yx ε zx
ε xy ε xz ε yy ε yz . ε zy ε zz
(3)
108
Each element of σ and ε can be non-zero value. For example, a non-zero σ xy means Ey
109
produces a current in the x direction.
110
As we mentioned, the media in this article are considered to be non-magnetic. In this case we
111
only need to discuss σ and ε for isotropic and anisotropic media. For isotropic media, σ and
112
ε can be written as follows: 0 0 σ ISO σ = 0 σ ISO 0 , 0 0 σ ISO
ε ISO ε = 0 0
0
ε ISO 0
0 0 , ε ISO
(4)
113
where each tensor is diagonal and has three equal eigenvalues. Then σ and ε can be treated as
114
two scalar values, σ ISO and ε ISO , respectively. The Maxwell’s curl equations for isotropic media
115
can be written as
∂H , ∂t ∂E ∇ × H = σ ISO E + ε ISO , ∂t ∇ × E = − µ0
116
which are similar to Eqs. (1) and (2) but easier to solve.
6
(5) (6)
117
Anisotropy, which is the main topic of this study, is more complicated. There are different types
118
of anisotropy, such as transverse isotropy (TI) and orthorhombic anisotropy. In this paper, we only
119
show the results of media exhibiting TI (TI media). When a material exhibitsTI, its physical
120
properties are isotropic in certain planes but different along the direction that is orthogonal to
121
these planes. This direction is also called symmetry axis. For TI media, there is always a certain
122
coordinate system that makes both σ and ε diagonal. Among the three diagonal elements
123
(which are also the eigenvalues of the matrix) of each tensor, two and only two are equal. In other
124
words,
σ 1 0 0 σ = 0 σ 2 0 , 0 0 σ 3
ε1 0 ε = 0 ε 2 0 0
0 0 , ε 3
(7)
125
where σ 1 = σ 2 , ε1 = ε 2 . This means the plane with one and two axes defined is the isotropic
126
plane, and the third axis is the symmetry axis.
127
128 129 130 131 132
a
b
c
Fig. 1. a VTI model. b HTI model. c TTI model.
133
For further discussion, we define a Cartesian coordinate system whose x- and y- axes are in the
134
horizontal plane and the z-axis is vertical and downward. We assume the acquisition survey is
135
located on the surface (z = 0). According to the direction of the symmetry axis, there are three
7
136
types of TI in general: vertical transverse isotropy (VTI), horizontal transverse isotropy (HTI) and
137
tilted transverse isotropy (TTI), which correspond to the vertical, horizontal and tilted symmetry
138
axis respectively (Fig. 1). The σ and ε of VTI, HTI, and TTI are (Carcione and Schoenberg,
139
2000)
σVTI
σ HTI
σTTI
0 σ xx 0 = 0 σ xx 0 , 0 0 σ zz σ xx 0 = 0 σ zz 0 0
0 0 , σ zz
σ xx σ xy σ xz = σ yy σ yz , σ zz
εVTI
ε HTI
ε HTI
ε xx = 0 0
0
ε xx 0
0 0 , ε zz
ε xx = 0 0
ε zz
ε xx =
ε xy ε xz ε yy ε yz , ε zz
0 0
0 0 , ε zz
(8)
140
where σ and ε are diagonal for VTI and HTI and symmetric for TTI. It is not so
141
straightforward to obtain the σ and ε directly. Usually we first define a VTI medium and
142
rotate it in the 3D space to obtain a TTI medium. This rotation can be demonstrated using Euler's
143
rotation theorem: σTTI = R T σVTI R , εTTI = R T εVTI R ,
144
(9)
where R = R x R y R z and
0 1 R x = 0 cos α 0 − sin α
0 cos β sin α , R y = 0 sin β cos α
0 − sin β cos γ 1 0 , R z = − sin γ 0 0 cos β
sin γ cos γ 0
0 0 . 1
145
Angles α , β , and γ are the rotation angles along the x, y, and z axes, respectively.
146
2.2. Maxwell’s equations for isotropic, dispersive medium
147
For an isotropic, linear, non-magnetic, dispersive medium, the Maxwell’s equations are 8
(10)
∇ × E = − µ0
∂H , ∂t
N ∂E +ε 0 ∑ J s , ∂t s =1
∇ × H = σE + ε 0ε ∞
J s = χ s (t) ∗
∂E , ∂t
(11)
(12)
(13)
148
where ε 0 is the vacuum permittivity, which is equal to 8.854 ×10 −12 F/m, ε ∞ is the relative
149
electric permittivity for infinity frequency, P is the polarization density for each dispersive pole
150
and N is the total number of dispersive poles. As we only study the one pole Debye case, we set s
151
to 1 in the following equations and applications.
152
χ (t) is an inclusive function that can be written for Debye medium as χ (t) =
∆ε
τ0
e− t /τ 0 ,
(14)
153
where ∆ε = ε z − ε ∞ , and ε z is the relative permittivity for zero frequency, and τ 0 is the
154
relaxation time. Further detail can be found in Giannakis and Giannopoulos (2014).
155
2.3. Yee’s FDTD method for anisotropic media
156
The finite-difference technique is used for numerical analysis of differential equations.
157
Maxwell’s equations are time-dependent partial derivative equations that can be solved using
158
Yee’s FDTD method. Yee’s method is grid-based, and the grid system is called Yee cell. Fig. 2
159
shows the Yee cell for the 3D EM case. As shown in Fig. 2, the Yee cell staggers electric field
160
components, magnetic field components, and material properties on the grids in time and space.
161
Yee’s algorithm discretizes time and space partial derivatives using central-difference
162
approximations with respect to the Yee cell. To update the electric fields, we need the stored value
163
of the electric fields and the magnetic fields located on the staggered grids around them. The
9
164
update made to the magnetic fields is similar. To build a model for EM numerical modeling, we
165
usually set the electrical and magnetic properties on the regular grids. To calculate the EM field
166
components on the grid half a spacing away, interpolations of σ and ε are necessary.
167 168 169
Fig. 2. Wave field components distribution of 3D Yee grid.
170
As magnetic permeability is assumed to be a scalar value (the media is non-magnetic), Eq. (1)
171
for the anisotropic case is exactly the same as that for the isotropic case. For anisotropic media,
172
we will use a finite difference expression for Eq. (2). First, we solve the temporal derivatives in
173
Eq. (2). By defining E at time n∆t (where n is a integer value) and H at half a time step
174
away, the temporal derivatives in Eq. (2) can be estimated by
∂E ∂t 175
Furthermore, E at time n +
n+
1 2
=
(
)
1 n +1 n E −E . ∆t
1 in Eq. (2) can be estimated as (Schneider and Hudson, 1993) 2 1 1 n +1 n+ n E 2 = E +E . (16) 2
(
176
(15)
)
According to Eq. (11) and Eq. (12), Eq. (2) can be written as
10
E
177
n +1
= F E +G n
(∇ × H )
n +1 2
,
(17)
where -1
ε σ ε σ F=G − , G = + . ∆t 2 ∆t 2
(18)
178
In Eq. (18) ∆t is the time sampling interval. Notice that Eq. (18) is not unique to anisotropic
179
media.
180 181
Next, we allocate wave fields and media properties spatially according to Yee’ cell. It is evident that F and G have the same symmetry as σ and ε , and can be written as:
Fxx F =
Fxz Fyz , Fzz
Fxy Fyy
Gxx G =
Gxy G yy
Gxz Gyz . Gzz
(19)
182
It was demonstrated earlier that VTI and HTI media have similar diagonal σ and ε values as
183
isotropic media. This led to similar finite-difference equations for VTI, HTI, and isotropy
184
(Schneider and Hudson, 1993). It would be more meaningful to discuss the general anisotropic
185
case. We only give the detailed equation for E x in Eq. (15) for brevity. The equations for E y and
186
E z can be derived in a similar fashion. As shown in Fig. 3, we put E x at points ( i∆x, j ∆y, k ∆z ) ,
187
where i, j , k are indexes and ∆x, ∆ y, ∆ z are spatial intervals in x, y and z direction. Then,
188
E x (x, y, z) can be written as E x i , j ,k .
Ex
n +1 i , j ,k
= Fxx i , j ,k E x
n i , j ,k
+ Fxy
i , j ,k
Ey
n i , j ,k
n +1 2
+ Fxz
∂H z ∂H y + Gxx i , j ,k − + Gxy ∂z i , j ,k ∂y 189
i , j ,k
Ez
n i , j ,k
n +1 2
∂H x ∂H z − + Gxz i , j , k ∂z ∂x i , j ,k
n +1 2
∂H y ∂H x − i , j ,k ∂y i , j ,k ∂x
.
(20)
According to the Yee cell (Fig. 3), electric and magnetic field quantities are at staggered grids.
11
190
Therefore, E y
n i , j ,k
n
, Ez
i , j ,k
∂H x , ∂z
n +1 2
i , j ,k
∂H x , ∂y
n +1 2
, i , j ,k
∂H y
n +1 2
∂x
i , j ,k
n +1 2
∂H z , and ∂x
cannot be obtained i , j ,k
191
directly from the Yee cell. They should be interpolated from neighbouring quantities (Schneider
192
and Hudson, 1993). For Eq. (20), the approximations for the electric field quantities are as
193
follows:
Ey Ez
n
(
i , j ,k
n i , j ,k
)
n n n n 1 Ey + Ey + Ey + Ey , i +1 2, j +1 2,k i +1 2, j −1 2,k i −1 2, j +1 2,k i −1 2, j −1 2,k 4 1 n n n n = Ez + Ez + Ez + Ez . i +1 2, j ,k +1 2 i +1 2, j ,k −1 2 i −1 2, j ,k +1 2 i −1 2, j ,k −1 2 4
=
(21)
)
(
194
Spatial derivatives of magnetic fields estimated by second finite difference operator require
195
interpolations:
∂H x ∂z
∂H x ∂y
n +1 2
= i , j ,k
n +1 2
= i , j ,k
∂H y ∂x
∂H z ∂x
1 (H x 4∆z
n +1 2
− Hx
n +1 2
1 (H x 4∆z
n +1 2
− Hx
n +1 2
n +1 2
i +1 2, j +1 2,k −1 2
n +1 2
− Hx
n +1 2
+ Hx
n +1 2
− Hx
n +1 2
i +1 2, j +1 2,k +1 2
i +1 2, j −1 2,k +1 2
=
1 (H y 4∆x
n +1 2
=
1 (H z 4∆x
n +1 2
i , j ,k
n +1 2
i , j ,k
i +1 2, j +1 2,k +1 2
+ Hx
i +1 j ,k +1 2
i +1, j +1 2,k
i +1 2, j −1 2,k +1 2
i +1 2, j −1 2,k −1 2
i −1 2, j +1 2,k +1 2
i −1 2, j −1 2,k +1 2
+ Hy
+ Hz
n +1 2 i +1, j ,k −1 2
n +1 2 i +1, j −1 2,k
+ Hx
n +1 2
− Hx
n +1 2
i −1 2, j +1 2,k +1 2
i −1 2, j +1 2,k −1 2
+ Hx
n +1 2
− Hx
n +1 2
− Hy
− Hz
i +1 2, j +1 2,k −1 2
i +1 2, j −1 2,k −1 2
n +1 2 i −1, j ,k +1 2
n +1 2 i −1, j +1 2,k
+ Hx
n +1 2
− Hx
n +1 2
i −1 2, j −1 2,k +1 2
i −1 2, j −1 2,k −1 2
+ Hx
n +1 2
− Hx
n +1 2
− Hy
− Hz
, )
i −1 2, j +1 2,k −1 2
i −1 2, j −1 2,k −1 2
n +1 2
(22) )
),
i −1, j ,k −1 2
n +1 2 i −1, j −1 2,k
).
196
It is evident that when adopting higher order finite difference operators, more interpolations are
197
required.
198
12
199
3. 3D anisotropic and dispersive modeling using the RSG scheme
200
Yee’s algorithm has proven to be a robust, stable and accurate method in isotropic EM wave
201
simulation. However, as we mentioned, the anisotropic case requires more interpolations of wave
202
fields for Yee’s algorithm. We therefore introduced the RSG scheme, which is free of this
203
inconvenience. We also adopted the C-PML as a boundary condition.
204
3.1. RSG scheme and comparison with the Yee cell
205
In the Yee cell, all the components of the EM wave fields are put on different, staggered grids.
206
For the isotropic case, this discretization strategy will result in interpolations of model properties
207
2*ngrids times, where ngrids is the number of grids. For the anisotropic case, interpolations of the
208
wave fields are required. The RSG scheme puts all components of same type on the same grids,
209
and put other fields on the rotated staggered grids. Fig. 4 is the RSG for 3D case. As shown in Fig.
210
4, all the electric field components and σ and ε are put on the regular grid nodes, and all
211
magnetic field components and µ0 are put on grids in the centre of regular rectangular cell. The
212
spatial derivatives in Eq. (16) are rebuilt and will be demonstrated later.
213
214 215 216 217 218
Fig. 3. Wave field components distribution of the 3D RSG cell.
13
219 220 221 222 223 224 225 226
Fig. 4. Four directions in the 3D RSG cell.
3.2. Discretization of derivations using the RSG technique
r r r r First, four vectors l1 , l2 , l3 , and l4 are defined in the RSG cell, which are shown in Fig. 4. r r r Based on the vector addition, these directions can be expressed by x , y , and z :
r ∆x r ∆y r ∆z r l1 = x+ y+ z, ∆l ∆l ∆l r ∆x r ∆y r ∆z r l2 = x+ y− z, ∆l ∆l ∆l r ∆x r ∆y r ∆z r l3 = x− y+ z, ∆l ∆l ∆l r ∆x r ∆y r ∆z r l4 = x− y− z, ∆l ∆l ∆l
(23)
227
r r r where ∆x , ∆y , and ∆z are the spatial intervals in the x , y , and z directions, respectively;
228
∆l = ∆x 2 + ∆y 2 + ∆z 2 . Then, we have
r 1 ∆l r r r r x= l1 + l2 + l3 + l4 , 4 ∆x
(
)
r 1 ∆l r r r r y= l1 + l2 − l3 − l4 , 4 ∆y
(
)
(24)
r 1 ∆l r r r r z= l1 − l2 + l3 − l4 . 4 ∆z
(
)
229
Therefore the spatial derivatives of wave fields f ( x, y, z, t ) (can represent electrical or
230
r r r magnetic wave fields components) along x , y , and z directions can be expressed as a linear
231
r r r r combination of derivatives along the l1 , l2 , l3 , ad l4 directions (Saenger et al., 2000):
14
∂f 1 ∆l ∂f ∂f ∂f ∂f = + + + , ∂x 4 ∆x ∂l1 ∂l2 ∂l3 ∂l4 ∂f 1 ∆l ∂f ∂f ∂f ∂f = − − + , ∂y 4 ∆x ∂l1 ∂l2 ∂l3 ∂l4
(25)
∂f 1 ∆l ∂f ∂f ∂f ∂f = + − − . ∂z 4 ∆x ∂l1 ∂l2 ∂l3 ∂l4 232
The derivatives ∂f ∂l1 , ∂f ∂l2 , ∂f ∂l3 , and ∂f ∂l4 are estimated using: ∂f ∂l1 ∂f ∂l2 ∂f ∂l3 ∂f ∂l4
=
1 M (M) ∑ Cm f ∆l m =1
= i , j ,k
1 M (M) ∑ Cm f ∆l m =1
i , j ,k
1 M = ∑ Cm(M) f ∆l m =1
i , j ,k
i+
2 m −1 2 m −1 2 m −1 ,j+ ,k + 2 2 2
i+
2 m −1 2 m −1 2 m −1 ,j+ ,k − 2 2 2
−f
i−
2 m −1 2 m −1 2 m −1 ,j− ,k − 2 2 2
,
−f
i−
2 m −1 2 m −1 2 m −1 ,j− ,k + 2 2 2
, (26)
= i , j ,k
1 M (M) ∑ Cm f ∆l m =1
i+
2 m −1 2 m −1 2 m −1 ,j− ,k + 2 2 2
i+
2 m −1 2 m −1 2 m −1 ,j− ,k − 2 2 2
i−
2 m −1 2 m −1 2 m −1 ,j+ ,k − 2 2 2
,
−f
i−
2 m −1 2 m −1 2 m −1 ,j+ ,k + 2 2 2
,
( i∆x, j ∆y, k ∆z ) .
233
where f
234
order value; m = 1, 2,L , M ; Cm(M) are the differential coefficients, whose values can be obtained
235
by solving:
i, j ,k
represents f at position
−f
3 1 1 33 M M 2 M −1 1 3
M is half of the spatial difference
L 2 M − 1 C1(M) 1 L (2 M − 1)3 C2(M) 0 = . M M O M L (2 M − 1)2 M −1 CM(M) 0
(27)
236
Substituting Eqs. (25)–(27) into Eq.(20), we can obtain the updated equation for E x . The updated
237
equation for the other wave field components can be derived in a similar fashion. It is evident that
238
no interpolations for model properties are necessary, which saves some computational time
239
compared with Yee’s method. As all the electric fields are in the same grids, the interpolations that
15
240
Eq. (21) carries out are no longer necessary for Eq. (20). The flow chart of the presented method
241
is shown in Fig. 5.
242 243 244 245
Fig. 5. Flow chart of the RSG-FDTD method.
3.3. Numerical stability and dispersion criteria
246
To ensure the numerical simulation is stable, appropriate time and spatial discretization
247
intervals should be chosen (Taflove, 2000; Georgakopoulos et al., 2002). For Yee’s method
248
(second order in time and space) in 3D, the numerical stability constraint derived by Taflove and
249
Brodwin (1975) is
∆tvmax 1 ≤ , ∆h 3 250
where
vmax
251
∆h = ∆x = ∆y = ∆z is the spatial interval.
(28)
is the maximum EM wave phase velocity within the given model;
16
252 253
The numerical stability constraint for the RSG method for modeling the elastic wave is given by Saenger et al. (2000):
∆tv p ∆h
≤1
M
∑C m =1
(M) m
,
(29)
254
where Cm(M) are differential coefficients; v p denotes the P wave velocity. The numerical stability
255
constraint for RSG can be derived following the procedure that Taflove and Brodwin (1975)
256
proposed (see Appendix for details). The numerical stability constraint for the RSG method in EM
257
wave modeling is: ∆tvEM ≤1 ∆h
258
∑C m =1
(M) m
,
(30)
where vEM is the velocity of the EM wave, which can be calculated using vEM =
259
M
c
εr
=
0.3
m/ns,
εr
(31)
where c is the speed of light in vacuum. Note that Eq. (30) is only valid for low-loss materials.
260
Ideally, despite the limitation of spatial resolution, we tend to set ∆h and ∆t as large as
261
possible (under the numerical stability constraint in Eq.(31)) to accelerate the modeling process.
262
However, if ∆h is too large, the numerical dispersion of electric and magnetic fields will occur.
263
To control this, five field samples per minimum wavelength must be guaranteed for second order
264
in time and forth order in space strategy (Irving and Knight, 2006).
265
4. Numerical tests based on homogeneous media
266 267
To test the performance of the proposed method, the code was validated with the well-known open source software
GprMax. Two tests using one anisotropic medium and one dispersive
17
268
medium were performed. The media were homogeneous, linear, and non-magnetic. The EM
269
waves were triggered by electric currents on E z at the centre of both models. The wavelets of the
270
sources were rickers with central frequencies both set to 1 GHz.
271
4.1 Anisotropic case
272
We tested the proposed method first with a VTI medium. We used a 200 × 200 × 200 grid set
ε xx = ε yy = 25ε 0 , ε zz = 9ε 0 ,
273
with spacing of 0.005 m. The model properties were
274
σ xx = σ yy = 1 mS / m and σ zz = 0.5 mS / m . The time step we used was 9.62917 × 10−12 s. We
275
recorded EM waves using one receiver located at
276
comparison of the results obtained with GprMax and this method. The waveforms recorded on the
277
E x , E y , and E z components all show good agreement, which indicates that the proposed RSG
278
scheme is able to simulate EM waves propagating in anisotropic media.
( 0.5 m, 0.8 m, 0.5 m ) .
Fig. 6 shows the
279 280
Fig. 6.
281
Comparison between the RSG scheme (red circle) and GprMax (blue solid line) of the VTI 18
282
anisotropic case. The waveforms recorded on the E x , E y , and E z components are shown in a,
283
b, and c, respectively.
284
To test the validity, we then tested the method with an HTI medium. We used the same grid set
285
and, spatial and time step. The parameters of the HTI medium were as follows: ε xx = 9ε 0 ,
286
ε yy = ε zz = 25ε 0 , σ xx = 0.5 mS / m , and σ yy = σ zz = 1mS / m . The EM wave was recorded by a
287
receiver 0.3 m away from source in the y-direction. Fig. 7 shows a comparison of the results
288
obtained with GprMax and this method. The E x , E y , and E z components of the waveforms
289
show good agreement. The directed waves of the E x and E z components have different arrival
290
times, which may indicates the wavelet split into two EM waves with different velocities.
291 292
Fig. 7.
293
Comparison between the RSG scheme (red circle) and GprMax (blue solid line) of the HTI
19
294
anisotropic case. The E x , E y , and E z waveform components are shown in a, b, and c,
295
respectively.
296 297
4.2 Dispersive case
298
Then, we tested the RSG scheme in a dispersive case. The medium was an isotropic, one-pole
299
Debye case. The model properties were as follows: ε ∞ = 4.9 , ε z = 80.1 , σ = 0 mS / m , and
300
τ 0 = 9.231× 10−12 s . We use a 250 × 250 × 250 grid set with a spacing of 0.002 m. The time step
301
we used was 3.85167 × 10 −12 s. The EM waves were recorded using one receiver located at
302
( 0.25m,
303
agreement is observed, which validates the proposed method.
0.4m, 0.25m ) . Fig. 8 compares the results obtained with GprMax and this method; good
304 305
Fig. 8.
20
306
Comparison between the RSG scheme (red circle) and GprMax (blue solid line) of the dispersive
307
case. The E x , E y , and E z components of the waveform are shown in a, b, and c, respectively.
308
5. Numerical test with a complex case
309
We tested the RSG FDTD method with a complex case. Fig. 9 shows the model. Two cylinders
310
that simulate an air-filled pipe and a water-filled pipe were elongated along the x-direction. The
311
centre of the cross-section of the air-filled pipe and water-filled pipe was located at (y = 0.35 m, z
312
= 0.4 m) and (y = 0.65 m, z = 0.4 m), respectively. The radii were 0.05 m for both cylinders. The
313
relative dielectric permittivities of air and water were 1 and 81, respectively. The electric
314
conductivities of air and water were 0 and 0.01 mS/m. Sources were located on the surface every
315
0.05 m from 0.05 to 0.95 m along the y-direction. We use a 1 GHz ricker wavelet to excite the EM
316
wave. The time and spatial spacing were 5 × 10-12 s and 0.005 m, respectively. We set the
317
background media as homogeneous isotropic, anisotropic, and dispersive to simulate clay, finely
318
stratified sand, and moist sand, respectively. The model properties for clay were ε r = 9 and
319
σ = 2 mS / m . The properties for finely stratified sand were ε r xx = ε ryy = 15 , ε rzz = 9 ,
320
σ xx = σ yy = 1 mS/m, and σ zz = 2 mS/m. The properties for moist sand were ε z = 15ε 0 ,
321
ε ∞ = 9ε 0 , τ 0 = 15 × 10−9 s, and σ = 2 mS / m .
21
322 323
Fig. 9.
324
Complex model. An air-filled and a water-filled pipe are buried underground.
325 326
a
b
327
328 329
c
22
330
Fig. 10.
331
Common offset sections of isotropic (a), anisotropic (b), and dispersive (c) background media
332
models.
333
The common offset sections are shown in Fig. 10. The source-receiver offset was 0.05m. The
334
reflections events of two pipes were different in each case, which is caused by the different
335
dielectric permittivities and conductivities of air and water. Other than this, the travel time was
336
evidently different between anisotropic and isotropic cases. This illustrates that the EM wave
337
propagates with different velocities in anisotropic case. Compared with the isotropic case, the EM
338
wave energy was weaker in the dispersive case, which was expected. The above examples
339
illustrate that the proposed method can provide a reasonable solution for 3D anisotropic and
340
dispersive problems.
341
6. Discussion and conclusions
342
In this paper, we describe a new method known as the RSG-FDTD method for modeling EM
343
waves inside anisotropic and dispersive materials. This method defines new directions for
344
approximating spatial partial derivatives. This method may be superior to Yee’s method when
345
simulating EM waves propagating inside anisotropic media because it avoids the interpolation of
346
wave field components. Furthermore, it requires fewer interpolations of media properties than
347
Yee’s method.
348
Using this approach, we generated wave fields and records of electric fields components in
349
homogeneous anisotropic and dispersive media, as well as in a complex model. The results
350
showed good agreement with those obtained with GprMax, which highlights the validity of the 23
351
proposed method. The results of common offset sections reveal that anisotropy and dispersion can
352
apparently influence the interpretation of GPR data. The analysis of the arrival time is one of the
353
easiest and most intuitive ways to evaluate the anisotropy, and it can be used to obtain wave speed
354
expediently. The results reveal that the velocity of an EM wave travelling in anisotropic media
355
depends on the propagation direction.
356
We only studied the TI anisotropy (also known as uniaxial anisotropy) and assumed the
357
conductivities are relatively small. In the future, we plan to study the biaxial anisotropic model
358
( ε xx ≠ ε yy ≠ ε zz , σ xx ≠ σ yy ≠ σ zz ). Furthermore, we will study the effect of non-negligible
359
conductivities. The dispersion modeling will include Lorentz and Drude functions in the future as
360
well. In addition, a better absorbing boundary , such as a time-synchronized CPML (Giannakis
361
and Giannopoulos, 2015), will be adopted. These studies may help us acquire a deeper
362
understanding of underground media.
363
Acknowledgement
364
This research is supported by the National Key R&D Program of China (Grant No.
365
2018YFB0605503) and, the Open Fund of State Key Laboratory of Water Resource Protection
366
and Utilization in Coal Mining (Grant No. SHJT-17-42.16).
367
Computer Code Availability
368
The name of the code we developed for this method is GPR_3D_ANI_RSG_CPML.cpp. It was
369
developed by Dr. Yongxu Lu (one of the authors). His contact address is Room 410, Minzu
370
Building, Ding No.11 Xueyuan Road, Haidian District, Beijing 100083 P. R.China. His telephone 24
371
number is +86 15650716135, and his e-mail is [email protected]. The code was firstly made
372
available at 2018 and can be run on both Linux/Windows systems. A C++ development
373
environment such as G++ in Linux is needed. The code can be run on computers that support
374
Ubuntu 16.04 or Windows 7. However, the memory size limits the size of the model. The program
375
is
376
https://github.com/samhomes/GPR-RSG-ANI.git.
377
Appendix A: Derivation of the stability criterion for RSG
378 379
about
36
KB.
Furthermore,
the
code
382
∂ ( H + jE ) , ∂t
(A.1)
(A.2)
considering the following pair of eigenvalue problems:
V = ΛV ,
(A.3a)
numerical
(A.3b)
To ensure the algorithm’s stability, the imaginary part of Λ must fall within the range (Taflove and Brodwin, 1975) −
385
Github:
The stability of a particular numerical representation of Eq. A.2 can be examined simply by
j∇ numerical × V = ΛV .
384
on
or more simply as:
∂ ∂t
383
available
rewrite Maxwell’s equations as:
j∇ × V = ∂V ∂t , where V = H + jE . 381
is
We let µ = 1, ε = 1, σ = 0 and vEM = 1 as Taflove and Brodwin (1975) proposed. We can
j∇ × ( H + jE ) = 380
source
2 2 ≤ Im ( Λ ) ≤ . ∆t ∆t
We let
25
(A.4)
V I , J , K = V0 e
(
j k%x I ∆x + k% y J ∆y + k%z K ∆z
)
(A.5)
,
386
represent an arbitrary lattice spatial mode. Using the second order space differencing of RSG (Eqs.
387
25) – (26) to implement the derivatives of the curl operator, Eq. (A.3b) yields: xˆ % ∆x % ∆y % ∆z ∆x sin k x 2 cos k y 2 cos k z 2 + yˆ ∆x % ∆y % ∆z −2 cos k%x sin k y cos k z + × V I , J ,K = Λ V I ,J ,K , 2 2 2 ∆y zˆ ∆x % ∆y % ∆z cos cos k%x ky sin k z 2 2 2 ∆z
(A.6)
388
where xˆ , yˆ , and zˆ are unit vectors in the x-, y-, and z-coordinate directions. After performing
389
the cross product and writing out the x, y, and z component equations, we can obtain the
390
following:
1 ∆x 2 % ∆y 2 % sin 2 k%x cos k y cos k z 2 2 2 ( ∆x ) 1 ∆x 2 % ∆y 2 % Λ 2 = −4 cos 2 k%x sin k y cos k z 2 ( ∆y ) 2 2 1 ∆x 2 % ∆y 2 % cos 2 k%x cos k y sin k z 2 2 2 ( ∆z ) 391 392
394
( ∆h = ∆x = ∆y = ∆z ), it is clear that: 2 2 ≤ Im ( Λ ) ≤ . ∆h ∆h
396
(A.8)
To satisfy the stability condition (Eq. (A.4)) for the arbitrary lattice spatial mode, and for arbitrary value of vEM , we set:
∆tvEM ≤1. ∆h 395
(A.7)
Assuming the worst case for the trigonometric functions in Eq. A.7 and for equal grid spacing
− 393
∆z + 2 ∆z + . 2 ∆z 2
(A.9)
For arbitrary order of spatial difference, the stability condition can be obtained in the same manner. It yields: 26
∆tvEM ≤1 ∆h
M
∑C m =1
(M) m
.
(A.10)
397 398
References
399
Alsunaidi, M.A., Al-Jabr, A.A., 2009. A general ADE-FDTD algorithm for the simulation of dispersive
400
structures. IEEE photonics technology letters, 21(12), 817–819.
401
Beggs, J.H., Luebbers R.J., Yee K.S., Kunz K.S., 1992. Finite-difference time-domain implementation of
402
surface impedance boundary-conditions. IEEE Transactions on Antennas and Propagation 40, 49-56.
403
Beker, B., Umashankar, K.R., Taflove, A., 1986. Numerical analysis and validation of the combined field
404
surface integral equations for electromagnetic scattering by arbitrary shaped two-dimensional anisotropic
405
objects. IEEE Trans. Antennas Propagat. 37, 1573-1581.
406
Berenger, J. P., 1994. A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys.
407
144, 185-200.
408
Bergmann, T., Robertsson, J. O., Holliger, K., 1998. Finite-difference modeling of electromagnetic wave
409
propagation in dispersive and attenuating media. Geophysics, 63(3), 856-867.
410
Bradford, J. H., Nichols, J., Harper, J. T., and Meierbachtol, T., 2013. Compressional and EM wave velocity
411
anisotropy in a temperate glacier due to basal crevasses, and implications for water content estimation, Annals
412
of Glaciology 54(64), 168-178.
413
Carcione, J.M., 2007. Wave Propagation in Anisotropic, Anelastic, Porousand Electromagnetic Media, 2nd edn,
414
Elsevier Science & Technology Wave fields in real media, 538 pp.
415
Carcione J.M., Schoenberg M.A., 2000. 3-D ground-penetrating radar simulation and plane-wave theory in
416
anisotropic media. Geophysics, 65(5), 1527-1541. 27
417
Chen, H., Huang, T., 1998. Finite-difference time-domain simulation of GPR data. Journal of Applied
418
Geophysics 40, 139-163.
419
Diamanti, N., Giannopoulos A., 2009. Implementation of ADI-FDTD subgrids in ground penetrating radar
420
FDTD models. Journal of Applied Geophysics, 67, 309-317.
421
Fan, G.X., Liu, Q.H., 2000. An FDTD algorithm with perfectly matched layers for general dispersive media,”
422
IEEE transactions on antennas and propagation, 48(5), 637–646.
423
Fang, H., Lin, G., 2012. Symplectic partitioned Runge–Kutta methods for two-dimensional numerical model of
424
ground penetrating radar. Computers & Geosciences 49, 323-329.
425
Georgakopoulos, S.V., Birtcher C. R., Balanis C.A., et al., 2002. Higher-order finite-difference schemes for
426
electromagnetic radiation, scattering and penetration, part I, Theory. IEEE Antenna’s and Propagation Magazine,
427
44(1), 134-142.
428
Giannakis, I., Giannopoulos, A., 2014. A novel piecewise linear recursive convolution approach for dispersive
429
media using the finite-difference time-domain method. IEEE transactions on antennas and propagation, 62(5),
430
2669-2678.
431
Giannakis, I., Giannopoulos, A., 2015. Time-Synchronised convolutional perfectly matched layer for improved
432
absorbing performance in FDTD. IEEE Antennas and Wireless Propagation Letters, 14, 690-693.
433
Giannakis, I., Giannopoulos, A., Warren, C., 2016. A Realistic FDTD Numerical Modeling Framework of
434
Ground Penetrating Radar for Landmine Detection. IEEE Journal of Selected Topics in Applied Earth
435
Observations and Remote Sensing, 9(1), 37-51.
436
Giannakis, I., Giannopoulos, A., Warren, C., 2019. Realistic FDTD GPR antenna models optimized using a
437
novel linear/nonlinear full-waveform inversion. IEEE transactions on geoscience and remote sensing, 57(3),
28
438
1768-1778.
439
Giannopoulos, A., 2005. Modelling ground penetrating radar by GprMax. Construction and Building Materials
440
19, 755–762.
441
Holliger, K., Bergmann, T., 2002. Numerical modeling of borehole georadar data. Geophysics 67, 1249-1257.
442
Irving, J., Knight, R., 2006. Numerical modeling of ground-penetrating radar in 2-D using MATLAB.
443
Computers & Geosciences 32, 1247-1258.
444
Kelley, D.F., Luebbers, R.J., 1996. Piecewise linear recursive convolution for dispersive media using FDTD.
445
IEEE transactions on antennas and propagation, 44(6), 792–797.
446
Kumlu, D., and Erer, I., 2019. Detection of buried objects in ground penetrating radar data using incremental
447
nonnegative matrix factorization. Remote Sensing Letters 10(7), 649-658.
448
Li, J., Zeng, Z., Huang, L., et al., 2012. GPR simulation based on complex frequency shifted recursive
449
integration PML boundary of 3D high order FDTD. Computers & Geosciences 49, 121-130.
450
Li, Y. F., Matar, O. B., 2010. Convolutional perfectly matched layer for elastic second-order wave equation. J.
451
Acoust. Soc. Am. 127, 1318-1327.
452
Luebbers, R. J., Hunsberger, F., 1992. FDTD for nth-order dispersive media. IEEE transactions on antennas and
453
propagation 40, 1297-1301.
454
Millington T.M., Cassidy N.J., 2010. Optimising GPR modelling: A practical, multi-threaded approach to 3D
455
FDTD numerical modelling. Computers & Geosciences 36, 1135-1144.
456
Ralchenko, M., Svilans, M., Samson, C., et al., 2015. Finite-difference time-domain modelling of
457
through-the-Earth radio signal propagation. Computers & Geosciences 85, 184-195.
458
Rodríguez-Abad, I., Martínez Sala, R. M., García García, F., and Capuz Lladró, R., 2010. Non-destructive
29
459
methodologies for the evaluation of moisture content in sawn timber structures: ground-penetrating radar and
460
ultrasound techniques. Near Surface Geophysics 8(6), 475-482.
461
Ronden, J.A., Gedney, S.D., 2000. Convolutional PML (CPML): An efficient FDTD implementation of the
462
CFS-PML for arbitrary media. Microwave Opt. Technol. Lett. 27, 334-339.
463
Saenger, E.H., Gold, N., Shapiro, S.A., 2000. Modeling the propagation of elastic waves using a modified
464
finite-difference grid. Wave Motion 31, 77-92.
465
Saenger, E.H., Bohlen, T., 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation
466
using the rotated staggered grid. Geophysics 69, 583-591.
467
Schneider, J., Hudson, S., 1993. The finite-difference time-domain method applied to anisotropic material. IEEE
468
Trans. Antennas Propagat. 41, 994–999.
469
Strickel, M. A., Taflove, A., 1990. Time-domain synthesis of broadband absorptive coatings for two-dimensional
470
conducting targets.
471
Taflove, A., 2000. Computational electrodynamics: the finite-difference time-domain method. Boston, MA :
472
Artech House.
473
Taflove, A., Brodwin, M.E., 1975. Numerical solution of steady-state electromagnetic scattering problems using
474
the time-dependent Maxwell’s equations. IEEE transactions on microwave theory and techniques MTT-23,
475
623-630.
476
Teixeira, F. L., Chew, W. C., Straka, M., et al., 1998. Finite-difference time-domain simulation of ground
477
penetrating radar on dispersive, inhomogeneous, and conductive soils. IEEE transactions on geoscience and
478
remote sensing, 36(6), 1928-1937.
479
Ursin, B., 1983. Review of elastic and electromagnetic wave propagation in horizontally layered media,
IEEE Trans. Antennas Propagat. 38, 1084-1091.
30
480
Geophysics, 48(8), 1063–1081.
481
Wang, T., Tripp, A. C., 1996. FDTD simulation of EM wave propagation in 3-D media. Geophysics 61, 110–
482
120.
483
Warren, C., Giannopoulos, A., 2011. Creating finite-difference time-domain models of commercial
484
ground-penetrating radar antennas using Taguchi's optimization method. Geophysics, 76(2), G37-G47.
485
Warren, C., Giannopoulos, A., Giannakis, I., 2016. gprMax: Open source software to simulate electromagnetic
486
wave propagation for Ground Penetrating Radar. Computer Physics Communications, 209, 163-170.
487
Warren, C., Giannopoulos, A., Gray, A., et al., 2019. A CUDA-based GPU engine for gprMax: Open source
488
FDTD electromagnetic simulation software. Computer Physics Communications, 237, 208-218.
489
Yee, K., 1966. Numerical solution of initial boundary value problems involving Maxwell’s equations in
490
isotropic media. IEEE Transactions on Antennas and Propagation AP-14, 302–307.
31
Highlights: 1. Introduced a rotated staggered grid method to simulate the EM waves propagating. 2. The method can be used to simulate EM waves in anisotropic media and dispersive media. 3. Test the algorithm with complex models.
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
S0098-3004(18)31218-4
DOI:
https://doi.org/10.1016/j.cageo.2019.104397
Reference:
CAGEO 104397
To appear in:
Computers and Geosciences
Received Date: 25 December 2018 Revised Date:
30 October 2019
Accepted Date: 3 December 2019
Please cite this article as: Lu, Y., Peng, S., Cui, X., li, D., Wang, K., 3D FDTD anisotropic and dispersive modeling for GPR using rotated staggered grid method, Computers and Geosciences (2020), doi: https://doi.org/10.1016/j.cageo.2019.104397. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
1
3D FDTD anisotropic and dispersive modeling for GPR using rotated
2
staggered grid method
3
Yongxu Lua,b,∗ , Suping Penga,b, Xiaoqin Cuib, Dong lia,b, Kang Wanga,b
4
a
5
Beijing 100083, China.
6
b
7
(Beijing), Beijing 100083, China.
College of Geoscience and Surveying Engineering, China University of Mining & Technology (Beijing),
State Key Laboratory of Coal Resources and Safe Mining, China University of Mining & Technology
8
Corresponding author:
9
Yongxu Lu, Ph.D. in Geophysics
10 11
Affiliation: College of Geoscience and Surveying Engineering, China University of Mining & Technology (Beijing).
12
Address: Ding No.11 Xueyuan Road, Haidian District, Beijing 100083 P. R.China
13
E-mail address: [email protected] (Y. Lu).
14
Tel.: +86 1565 071 6135
15
____________________________________
16
Yongxu Lu developed all of the theory, derived all the equations and wrote this paper and all the code.
17
Suping Peng and Xiaoqin Cui participated in developing the basic theory.
18
Dong Li and Kang Wang participated in debugging the code and revising this paper.
19 20 1
21
Abstract
22
The finite-difference time-domain (FDTD) method is widely used in numerical modeling of
23
ground penetrating radar (GPR). The well-known Yee’s method is the earliest and most popular
24
FDTD method. However, using Yee’s method to solve an anisotropic case requires considerable
25
interpolation of electric fields (assuming media to be non-magnetic). In this paper we present a
26
different approach called rotated staggered grid (RSG) method to simulate electromagnetic (EM)
27
wave propagation in anisotropic and dispersive media. RSG-FDTD places electric field
28
components on the same grid system and magnetic field components on the same grid system,
29
however, half a spacing away. The RSG-FDTD method was first used in elastic wave modeling. It
30
requires no interpolations of wave fields and is advantageous for modeling wave propagation in
31
anisotropic media. In this paper, the 3D anisotropic and dispersive EM wave equations are
32
presented, and the RSG-FDTD approximations are derived. The application of the new method is
33
tested using homogeneous media models. To validate the proposed method, results are compared
34
with those of GprMax, and good agreement is achieved. The common offset sections of complex
35
cases with anisotropy and dispersion are also tested and are compared with the corresponding
36
isotropic case.
37 38
Keywords: Finite-difference time-domain (FDTD), Rotated staggered grid (RSG) method, 3D
39
GPR anisotropic modeling, 3D GPR dispersive modeling
40
1. Introduction
41
Ground penetrating radar (GPR) is a well-known geophysical technique that is applied in 2
42
various applications across the fields of geology, archaeology, glaciology, and engineering. In
43
engineering, it is used to analyse asphalt, concrete, wood, and a few other materials
44
(Rodríguez-Abad et al., 2010; Kumlu and Erer, 2019). It transmits a high frequency (100 MHz to
45
4 GHz) electromagnetic (EM) signal into the earth or other target media, and receives the
46
reflections using a set of antennas. By processing and analysing the reflection signals, structures
47
and physical subsurface properties are detected.
48
Numerical modeling techniques can be used to understand the GPR response of certain models
49
or for testing new processing techniques (Giannakis et al., 2019). The finite-difference
50
time-domain (FDTD) method, which is straightforward and easy to be programmed, is well
51
established for EM wave simulation (Diamanti and Giannopoulos, 2009; Millington and Cassidy,
52
2010; Fang and Lin, 2012; Li et al., 2012; Ralchenko et al., 2015). Almost all FDTD methods
53
published are based on Yee’s algorithm (1966), which can be conveniently adopted. The basic idea
54
is using a staggered grid on which spatial derivatives are located between two original grid points.
55
Yee’s algorithm is adopted by many researchers to study high frequency EM waves propagating in
56
2D or 3D isotropic media (Wang and Tripp, 1996, Chen and Huang, 1998; Holliger and Bergman,
57
2002; Giannopoulos, 2005; Irving and Knight, 2006; Warren and Giannopoulos, 2011; Giannakis
58
et al., 2016; Warren et al., 2016; Warren et al., 2019).
59
As we mentioned, several papers discussing forward modeling of GPR focus on isotropic media.
60
However, natural or manmade materials are always anisotropic and/or dispersive. In homogenous
61
isotropic media, wave properties (like wave speed for example) are the same in all directions. In
62
anisotropic media, however, the wave properties are different depending on the propagation
3
63
direction. It is not appropriate to interpret data containing anisotropic information using isotropic
64
theory. EM anisotropy is the variation of electric conductivity, dielectric permittivity or magnetic
65
permeability with azimuth inside the considered material. Rodríguez-Abad et al. (2010) evaluated
66
moisture content in wood using GPR. They proposed analyzing the moisture content after
67
determining the wood anisotropic effects. Bradford et al. (2013) studied the EM wave velocity
68
anisotropy in temperate glaciers. The study of the propagation of EM waves in anisotropic media
69
may help us to learn more about the targets. Beker et al. (1989) and Strikel and Taflove (1990)
70
developed algorithms for anisotropic media based on Yee’s algorithm. However, the anisotropic
71
media adopted in their studies are vertical transverse isotropic (VTI), which means the electric
72
conductivity, dielectric permittivity or magnetic permeability tensors are diagonal. Schneider and
73
Hudson (1993) developed an algorithm for general anisotropic media. Their algorithm is also
74
based on Yee’s algorithm and requires interpolation of electrical and magnetic fields. Moreover,
75
they show results of their algorithm in just a one-dimension case.
76
In practice, several materials, such as water and soils, are dispersive (Teixeira et al., 1998).
77
Single and multi-pole Debye, Lorentz and Drude models are widely used to simulate dispersive
78
media (Bergmann et al., 1998). Several different FDTD based methods for EM wave simulation in
79
the dispersive model have been suggested (Alsunaidi and Al-Jabr, 2009; Giannakis and
80
Giannopoulos, 2014). The recursive convolution algorithms are efficient in handling the
81
dispersive case (Kelley and Luebbers, 1996; Fan and Liu, 2000; Giannakis and Giannopoulos,
82
2014).
83
In this study, we apply a rotated staggered grid (RSG) scheme first reported by Saenger
4
84
(Saenger and Gold, 2000; Saenger et al., 2004) for elastic waves to GPR applications. This grid
85
technique places components of the same wave fields on the same grids. For anisotropic media,
86
interpolation of EM wave fields is not required (Schneider and Hudson, 1993). In our simulation,
87
we place electric fields on a regular grid and magnetic fields on an “half-grid.” We also apply the
88
convolutional perfectly matched layer (C-PML) (Roden and Gedney, 2000; Irving and Knight,
89
2006; Li and Matar, 2009) as an absorbing boundary condition. C-PML, introduced by Roden and
90
Gedney (2000), dose not require splitting of components and has a better performance than a
91
perfectly matched layer (PML). First, we discuss the basic theory of GPR modeling in isotropic,
92
anisotropic and dispersive media, including Maxwell’s equation, and provide a brief introduction
93
of Yee’s FDTD method. Then we introduce the RSG method and derive an updated equation for
94
3D EM wave propagation in anisotropic and dispersive media. Furthermore, we give the
95
numerical stability and dispersion criteria for the RSG method. Finally, we present some examples
96
for testing the proposed GPR modeling method. A few of the results are compared with those
97
obtained with free software GprMax.
98
2. Basic theory of FDTD GPR modeling with anisotropy and dispersive
99
2.1. Maxwell’s equations for isotropic and anisotropic medium
100
GPR uses high frequency EM waves to detect media underground. The basic governing
101
equations for EM phenomena are the well-known Maxwell’s curl equations (Schneider and
102
Hudson, 1993), ∇ × E = −µ
∂H , ∂t
5
(1)
∇ × H = σE + ε
∂E , ∂t
(2)
103
where H and E are the magnetic field intensity and electric field strength vectors, respectively. σ ,
104
ε and µ are the electric conductivity tensor, dielectric permittivity tensor and magnetic
105
permeability tensor, respectively. In this paper we assume that the media are considered as
106
non-magnetic. Thus µ = µ0 , where µ0 = 1.256 × 10−6 H/m in vacuum. σ and ε can be written
107
as follows:
σ xx σ xy σ xz σ = σ yx σ yy σ yz , σ zx σ zy σ zz
ε xx ε = ε yx ε zx
ε xy ε xz ε yy ε yz . ε zy ε zz
(3)
108
Each element of σ and ε can be non-zero value. For example, a non-zero σ xy means Ey
109
produces a current in the x direction.
110
As we mentioned, the media in this article are considered to be non-magnetic. In this case we
111
only need to discuss σ and ε for isotropic and anisotropic media. For isotropic media, σ and
112
ε can be written as follows: 0 0 σ ISO σ = 0 σ ISO 0 , 0 0 σ ISO
ε ISO ε = 0 0
0
ε ISO 0
0 0 , ε ISO
(4)
113
where each tensor is diagonal and has three equal eigenvalues. Then σ and ε can be treated as
114
two scalar values, σ ISO and ε ISO , respectively. The Maxwell’s curl equations for isotropic media
115
can be written as
∂H , ∂t ∂E ∇ × H = σ ISO E + ε ISO , ∂t ∇ × E = − µ0
116
which are similar to Eqs. (1) and (2) but easier to solve.
6
(5) (6)
117
Anisotropy, which is the main topic of this study, is more complicated. There are different types
118
of anisotropy, such as transverse isotropy (TI) and orthorhombic anisotropy. In this paper, we only
119
show the results of media exhibiting TI (TI media). When a material exhibitsTI, its physical
120
properties are isotropic in certain planes but different along the direction that is orthogonal to
121
these planes. This direction is also called symmetry axis. For TI media, there is always a certain
122
coordinate system that makes both σ and ε diagonal. Among the three diagonal elements
123
(which are also the eigenvalues of the matrix) of each tensor, two and only two are equal. In other
124
words,
σ 1 0 0 σ = 0 σ 2 0 , 0 0 σ 3
ε1 0 ε = 0 ε 2 0 0
0 0 , ε 3
(7)
125
where σ 1 = σ 2 , ε1 = ε 2 . This means the plane with one and two axes defined is the isotropic
126
plane, and the third axis is the symmetry axis.
127
128 129 130 131 132
a
b
c
Fig. 1. a VTI model. b HTI model. c TTI model.
133
For further discussion, we define a Cartesian coordinate system whose x- and y- axes are in the
134
horizontal plane and the z-axis is vertical and downward. We assume the acquisition survey is
135
located on the surface (z = 0). According to the direction of the symmetry axis, there are three
7
136
types of TI in general: vertical transverse isotropy (VTI), horizontal transverse isotropy (HTI) and
137
tilted transverse isotropy (TTI), which correspond to the vertical, horizontal and tilted symmetry
138
axis respectively (Fig. 1). The σ and ε of VTI, HTI, and TTI are (Carcione and Schoenberg,
139
2000)
σVTI
σ HTI
σTTI
0 σ xx 0 = 0 σ xx 0 , 0 0 σ zz σ xx 0 = 0 σ zz 0 0
0 0 , σ zz
σ xx σ xy σ xz = σ yy σ yz , σ zz
εVTI
ε HTI
ε HTI
ε xx = 0 0
0
ε xx 0
0 0 , ε zz
ε xx = 0 0
ε zz
ε xx =
ε xy ε xz ε yy ε yz , ε zz
0 0
0 0 , ε zz
(8)
140
where σ and ε are diagonal for VTI and HTI and symmetric for TTI. It is not so
141
straightforward to obtain the σ and ε directly. Usually we first define a VTI medium and
142
rotate it in the 3D space to obtain a TTI medium. This rotation can be demonstrated using Euler's
143
rotation theorem: σTTI = R T σVTI R , εTTI = R T εVTI R ,
144
(9)
where R = R x R y R z and
0 1 R x = 0 cos α 0 − sin α
0 cos β sin α , R y = 0 sin β cos α
0 − sin β cos γ 1 0 , R z = − sin γ 0 0 cos β
sin γ cos γ 0
0 0 . 1
145
Angles α , β , and γ are the rotation angles along the x, y, and z axes, respectively.
146
2.2. Maxwell’s equations for isotropic, dispersive medium
147
For an isotropic, linear, non-magnetic, dispersive medium, the Maxwell’s equations are 8
(10)
∇ × E = − µ0
∂H , ∂t
N ∂E +ε 0 ∑ J s , ∂t s =1
∇ × H = σE + ε 0ε ∞
J s = χ s (t) ∗
∂E , ∂t
(11)
(12)
(13)
148
where ε 0 is the vacuum permittivity, which is equal to 8.854 ×10 −12 F/m, ε ∞ is the relative
149
electric permittivity for infinity frequency, P is the polarization density for each dispersive pole
150
and N is the total number of dispersive poles. As we only study the one pole Debye case, we set s
151
to 1 in the following equations and applications.
152
χ (t) is an inclusive function that can be written for Debye medium as χ (t) =
∆ε
τ0
e− t /τ 0 ,
(14)
153
where ∆ε = ε z − ε ∞ , and ε z is the relative permittivity for zero frequency, and τ 0 is the
154
relaxation time. Further detail can be found in Giannakis and Giannopoulos (2014).
155
2.3. Yee’s FDTD method for anisotropic media
156
The finite-difference technique is used for numerical analysis of differential equations.
157
Maxwell’s equations are time-dependent partial derivative equations that can be solved using
158
Yee’s FDTD method. Yee’s method is grid-based, and the grid system is called Yee cell. Fig. 2
159
shows the Yee cell for the 3D EM case. As shown in Fig. 2, the Yee cell staggers electric field
160
components, magnetic field components, and material properties on the grids in time and space.
161
Yee’s algorithm discretizes time and space partial derivatives using central-difference
162
approximations with respect to the Yee cell. To update the electric fields, we need the stored value
163
of the electric fields and the magnetic fields located on the staggered grids around them. The
9
164
update made to the magnetic fields is similar. To build a model for EM numerical modeling, we
165
usually set the electrical and magnetic properties on the regular grids. To calculate the EM field
166
components on the grid half a spacing away, interpolations of σ and ε are necessary.
167 168 169
Fig. 2. Wave field components distribution of 3D Yee grid.
170
As magnetic permeability is assumed to be a scalar value (the media is non-magnetic), Eq. (1)
171
for the anisotropic case is exactly the same as that for the isotropic case. For anisotropic media,
172
we will use a finite difference expression for Eq. (2). First, we solve the temporal derivatives in
173
Eq. (2). By defining E at time n∆t (where n is a integer value) and H at half a time step
174
away, the temporal derivatives in Eq. (2) can be estimated by
∂E ∂t 175
Furthermore, E at time n +
n+
1 2
=
(
)
1 n +1 n E −E . ∆t
1 in Eq. (2) can be estimated as (Schneider and Hudson, 1993) 2 1 1 n +1 n+ n E 2 = E +E . (16) 2
(
176
(15)
)
According to Eq. (11) and Eq. (12), Eq. (2) can be written as
10
E
177
n +1
= F E +G n
(∇ × H )
n +1 2
,
(17)
where -1
ε σ ε σ F=G − , G = + . ∆t 2 ∆t 2
(18)
178
In Eq. (18) ∆t is the time sampling interval. Notice that Eq. (18) is not unique to anisotropic
179
media.
180 181
Next, we allocate wave fields and media properties spatially according to Yee’ cell. It is evident that F and G have the same symmetry as σ and ε , and can be written as:
Fxx F =
Fxz Fyz , Fzz
Fxy Fyy
Gxx G =
Gxy G yy
Gxz Gyz . Gzz
(19)
182
It was demonstrated earlier that VTI and HTI media have similar diagonal σ and ε values as
183
isotropic media. This led to similar finite-difference equations for VTI, HTI, and isotropy
184
(Schneider and Hudson, 1993). It would be more meaningful to discuss the general anisotropic
185
case. We only give the detailed equation for E x in Eq. (15) for brevity. The equations for E y and
186
E z can be derived in a similar fashion. As shown in Fig. 3, we put E x at points ( i∆x, j ∆y, k ∆z ) ,
187
where i, j , k are indexes and ∆x, ∆ y, ∆ z are spatial intervals in x, y and z direction. Then,
188
E x (x, y, z) can be written as E x i , j ,k .
Ex
n +1 i , j ,k
= Fxx i , j ,k E x
n i , j ,k
+ Fxy
i , j ,k
Ey
n i , j ,k
n +1 2
+ Fxz
∂H z ∂H y + Gxx i , j ,k − + Gxy ∂z i , j ,k ∂y 189
i , j ,k
Ez
n i , j ,k
n +1 2
∂H x ∂H z − + Gxz i , j , k ∂z ∂x i , j ,k
n +1 2
∂H y ∂H x − i , j ,k ∂y i , j ,k ∂x
.
(20)
According to the Yee cell (Fig. 3), electric and magnetic field quantities are at staggered grids.
11
190
Therefore, E y
n i , j ,k
n
, Ez
i , j ,k
∂H x , ∂z
n +1 2
i , j ,k
∂H x , ∂y
n +1 2
, i , j ,k
∂H y
n +1 2
∂x
i , j ,k
n +1 2
∂H z , and ∂x
cannot be obtained i , j ,k
191
directly from the Yee cell. They should be interpolated from neighbouring quantities (Schneider
192
and Hudson, 1993). For Eq. (20), the approximations for the electric field quantities are as
193
follows:
Ey Ez
n
(
i , j ,k
n i , j ,k
)
n n n n 1 Ey + Ey + Ey + Ey , i +1 2, j +1 2,k i +1 2, j −1 2,k i −1 2, j +1 2,k i −1 2, j −1 2,k 4 1 n n n n = Ez + Ez + Ez + Ez . i +1 2, j ,k +1 2 i +1 2, j ,k −1 2 i −1 2, j ,k +1 2 i −1 2, j ,k −1 2 4
=
(21)
)
(
194
Spatial derivatives of magnetic fields estimated by second finite difference operator require
195
interpolations:
∂H x ∂z
∂H x ∂y
n +1 2
= i , j ,k
n +1 2
= i , j ,k
∂H y ∂x
∂H z ∂x
1 (H x 4∆z
n +1 2
− Hx
n +1 2
1 (H x 4∆z
n +1 2
− Hx
n +1 2
n +1 2
i +1 2, j +1 2,k −1 2
n +1 2
− Hx
n +1 2
+ Hx
n +1 2
− Hx
n +1 2
i +1 2, j +1 2,k +1 2
i +1 2, j −1 2,k +1 2
=
1 (H y 4∆x
n +1 2
=
1 (H z 4∆x
n +1 2
i , j ,k
n +1 2
i , j ,k
i +1 2, j +1 2,k +1 2
+ Hx
i +1 j ,k +1 2
i +1, j +1 2,k
i +1 2, j −1 2,k +1 2
i +1 2, j −1 2,k −1 2
i −1 2, j +1 2,k +1 2
i −1 2, j −1 2,k +1 2
+ Hy
+ Hz
n +1 2 i +1, j ,k −1 2
n +1 2 i +1, j −1 2,k
+ Hx
n +1 2
− Hx
n +1 2
i −1 2, j +1 2,k +1 2
i −1 2, j +1 2,k −1 2
+ Hx
n +1 2
− Hx
n +1 2
− Hy
− Hz
i +1 2, j +1 2,k −1 2
i +1 2, j −1 2,k −1 2
n +1 2 i −1, j ,k +1 2
n +1 2 i −1, j +1 2,k
+ Hx
n +1 2
− Hx
n +1 2
i −1 2, j −1 2,k +1 2
i −1 2, j −1 2,k −1 2
+ Hx
n +1 2
− Hx
n +1 2
− Hy
− Hz
, )
i −1 2, j +1 2,k −1 2
i −1 2, j −1 2,k −1 2
n +1 2
(22) )
),
i −1, j ,k −1 2
n +1 2 i −1, j −1 2,k
).
196
It is evident that when adopting higher order finite difference operators, more interpolations are
197
required.
198
12
199
3. 3D anisotropic and dispersive modeling using the RSG scheme
200
Yee’s algorithm has proven to be a robust, stable and accurate method in isotropic EM wave
201
simulation. However, as we mentioned, the anisotropic case requires more interpolations of wave
202
fields for Yee’s algorithm. We therefore introduced the RSG scheme, which is free of this
203
inconvenience. We also adopted the C-PML as a boundary condition.
204
3.1. RSG scheme and comparison with the Yee cell
205
In the Yee cell, all the components of the EM wave fields are put on different, staggered grids.
206
For the isotropic case, this discretization strategy will result in interpolations of model properties
207
2*ngrids times, where ngrids is the number of grids. For the anisotropic case, interpolations of the
208
wave fields are required. The RSG scheme puts all components of same type on the same grids,
209
and put other fields on the rotated staggered grids. Fig. 4 is the RSG for 3D case. As shown in Fig.
210
4, all the electric field components and σ and ε are put on the regular grid nodes, and all
211
magnetic field components and µ0 are put on grids in the centre of regular rectangular cell. The
212
spatial derivatives in Eq. (16) are rebuilt and will be demonstrated later.
213
214 215 216 217 218
Fig. 3. Wave field components distribution of the 3D RSG cell.
13
219 220 221 222 223 224 225 226
Fig. 4. Four directions in the 3D RSG cell.
3.2. Discretization of derivations using the RSG technique
r r r r First, four vectors l1 , l2 , l3 , and l4 are defined in the RSG cell, which are shown in Fig. 4. r r r Based on the vector addition, these directions can be expressed by x , y , and z :
r ∆x r ∆y r ∆z r l1 = x+ y+ z, ∆l ∆l ∆l r ∆x r ∆y r ∆z r l2 = x+ y− z, ∆l ∆l ∆l r ∆x r ∆y r ∆z r l3 = x− y+ z, ∆l ∆l ∆l r ∆x r ∆y r ∆z r l4 = x− y− z, ∆l ∆l ∆l
(23)
227
r r r where ∆x , ∆y , and ∆z are the spatial intervals in the x , y , and z directions, respectively;
228
∆l = ∆x 2 + ∆y 2 + ∆z 2 . Then, we have
r 1 ∆l r r r r x= l1 + l2 + l3 + l4 , 4 ∆x
(
)
r 1 ∆l r r r r y= l1 + l2 − l3 − l4 , 4 ∆y
(
)
(24)
r 1 ∆l r r r r z= l1 − l2 + l3 − l4 . 4 ∆z
(
)
229
Therefore the spatial derivatives of wave fields f ( x, y, z, t ) (can represent electrical or
230
r r r magnetic wave fields components) along x , y , and z directions can be expressed as a linear
231
r r r r combination of derivatives along the l1 , l2 , l3 , ad l4 directions (Saenger et al., 2000):
14
∂f 1 ∆l ∂f ∂f ∂f ∂f = + + + , ∂x 4 ∆x ∂l1 ∂l2 ∂l3 ∂l4 ∂f 1 ∆l ∂f ∂f ∂f ∂f = − − + , ∂y 4 ∆x ∂l1 ∂l2 ∂l3 ∂l4
(25)
∂f 1 ∆l ∂f ∂f ∂f ∂f = + − − . ∂z 4 ∆x ∂l1 ∂l2 ∂l3 ∂l4 232
The derivatives ∂f ∂l1 , ∂f ∂l2 , ∂f ∂l3 , and ∂f ∂l4 are estimated using: ∂f ∂l1 ∂f ∂l2 ∂f ∂l3 ∂f ∂l4
=
1 M (M) ∑ Cm f ∆l m =1
= i , j ,k
1 M (M) ∑ Cm f ∆l m =1
i , j ,k
1 M = ∑ Cm(M) f ∆l m =1
i , j ,k
i+
2 m −1 2 m −1 2 m −1 ,j+ ,k + 2 2 2
i+
2 m −1 2 m −1 2 m −1 ,j+ ,k − 2 2 2
−f
i−
2 m −1 2 m −1 2 m −1 ,j− ,k − 2 2 2
,
−f
i−
2 m −1 2 m −1 2 m −1 ,j− ,k + 2 2 2
, (26)
= i , j ,k
1 M (M) ∑ Cm f ∆l m =1
i+
2 m −1 2 m −1 2 m −1 ,j− ,k + 2 2 2
i+
2 m −1 2 m −1 2 m −1 ,j− ,k − 2 2 2
i−
2 m −1 2 m −1 2 m −1 ,j+ ,k − 2 2 2
,
−f
i−
2 m −1 2 m −1 2 m −1 ,j+ ,k + 2 2 2
,
( i∆x, j ∆y, k ∆z ) .
233
where f
234
order value; m = 1, 2,L , M ; Cm(M) are the differential coefficients, whose values can be obtained
235
by solving:
i, j ,k
represents f at position
−f
3 1 1 33 M M 2 M −1 1 3
M is half of the spatial difference
L 2 M − 1 C1(M) 1 L (2 M − 1)3 C2(M) 0 = . M M O M L (2 M − 1)2 M −1 CM(M) 0
(27)
236
Substituting Eqs. (25)–(27) into Eq.(20), we can obtain the updated equation for E x . The updated
237
equation for the other wave field components can be derived in a similar fashion. It is evident that
238
no interpolations for model properties are necessary, which saves some computational time
239
compared with Yee’s method. As all the electric fields are in the same grids, the interpolations that
15
240
Eq. (21) carries out are no longer necessary for Eq. (20). The flow chart of the presented method
241
is shown in Fig. 5.
242 243 244 245
Fig. 5. Flow chart of the RSG-FDTD method.
3.3. Numerical stability and dispersion criteria
246
To ensure the numerical simulation is stable, appropriate time and spatial discretization
247
intervals should be chosen (Taflove, 2000; Georgakopoulos et al., 2002). For Yee’s method
248
(second order in time and space) in 3D, the numerical stability constraint derived by Taflove and
249
Brodwin (1975) is
∆tvmax 1 ≤ , ∆h 3 250
where
vmax
251
∆h = ∆x = ∆y = ∆z is the spatial interval.
(28)
is the maximum EM wave phase velocity within the given model;
16
252 253
The numerical stability constraint for the RSG method for modeling the elastic wave is given by Saenger et al. (2000):
∆tv p ∆h
≤1
M
∑C m =1
(M) m
,
(29)
254
where Cm(M) are differential coefficients; v p denotes the P wave velocity. The numerical stability
255
constraint for RSG can be derived following the procedure that Taflove and Brodwin (1975)
256
proposed (see Appendix for details). The numerical stability constraint for the RSG method in EM
257
wave modeling is: ∆tvEM ≤1 ∆h
258
∑C m =1
(M) m
,
(30)
where vEM is the velocity of the EM wave, which can be calculated using vEM =
259
M
c
εr
=
0.3
m/ns,
εr
(31)
where c is the speed of light in vacuum. Note that Eq. (30) is only valid for low-loss materials.
260
Ideally, despite the limitation of spatial resolution, we tend to set ∆h and ∆t as large as
261
possible (under the numerical stability constraint in Eq.(31)) to accelerate the modeling process.
262
However, if ∆h is too large, the numerical dispersion of electric and magnetic fields will occur.
263
To control this, five field samples per minimum wavelength must be guaranteed for second order
264
in time and forth order in space strategy (Irving and Knight, 2006).
265
4. Numerical tests based on homogeneous media
266 267
To test the performance of the proposed method, the code was validated with the well-known open source software
GprMax. Two tests using one anisotropic medium and one dispersive
17
268
medium were performed. The media were homogeneous, linear, and non-magnetic. The EM
269
waves were triggered by electric currents on E z at the centre of both models. The wavelets of the
270
sources were rickers with central frequencies both set to 1 GHz.
271
4.1 Anisotropic case
272
We tested the proposed method first with a VTI medium. We used a 200 × 200 × 200 grid set
ε xx = ε yy = 25ε 0 , ε zz = 9ε 0 ,
273
with spacing of 0.005 m. The model properties were
274
σ xx = σ yy = 1 mS / m and σ zz = 0.5 mS / m . The time step we used was 9.62917 × 10−12 s. We
275
recorded EM waves using one receiver located at
276
comparison of the results obtained with GprMax and this method. The waveforms recorded on the
277
E x , E y , and E z components all show good agreement, which indicates that the proposed RSG
278
scheme is able to simulate EM waves propagating in anisotropic media.
( 0.5 m, 0.8 m, 0.5 m ) .
Fig. 6 shows the
279 280
Fig. 6.
281
Comparison between the RSG scheme (red circle) and GprMax (blue solid line) of the VTI 18
282
anisotropic case. The waveforms recorded on the E x , E y , and E z components are shown in a,
283
b, and c, respectively.
284
To test the validity, we then tested the method with an HTI medium. We used the same grid set
285
and, spatial and time step. The parameters of the HTI medium were as follows: ε xx = 9ε 0 ,
286
ε yy = ε zz = 25ε 0 , σ xx = 0.5 mS / m , and σ yy = σ zz = 1mS / m . The EM wave was recorded by a
287
receiver 0.3 m away from source in the y-direction. Fig. 7 shows a comparison of the results
288
obtained with GprMax and this method. The E x , E y , and E z components of the waveforms
289
show good agreement. The directed waves of the E x and E z components have different arrival
290
times, which may indicates the wavelet split into two EM waves with different velocities.
291 292
Fig. 7.
293
Comparison between the RSG scheme (red circle) and GprMax (blue solid line) of the HTI
19
294
anisotropic case. The E x , E y , and E z waveform components are shown in a, b, and c,
295
respectively.
296 297
4.2 Dispersive case
298
Then, we tested the RSG scheme in a dispersive case. The medium was an isotropic, one-pole
299
Debye case. The model properties were as follows: ε ∞ = 4.9 , ε z = 80.1 , σ = 0 mS / m , and
300
τ 0 = 9.231× 10−12 s . We use a 250 × 250 × 250 grid set with a spacing of 0.002 m. The time step
301
we used was 3.85167 × 10 −12 s. The EM waves were recorded using one receiver located at
302
( 0.25m,
303
agreement is observed, which validates the proposed method.
0.4m, 0.25m ) . Fig. 8 compares the results obtained with GprMax and this method; good
304 305
Fig. 8.
20
306
Comparison between the RSG scheme (red circle) and GprMax (blue solid line) of the dispersive
307
case. The E x , E y , and E z components of the waveform are shown in a, b, and c, respectively.
308
5. Numerical test with a complex case
309
We tested the RSG FDTD method with a complex case. Fig. 9 shows the model. Two cylinders
310
that simulate an air-filled pipe and a water-filled pipe were elongated along the x-direction. The
311
centre of the cross-section of the air-filled pipe and water-filled pipe was located at (y = 0.35 m, z
312
= 0.4 m) and (y = 0.65 m, z = 0.4 m), respectively. The radii were 0.05 m for both cylinders. The
313
relative dielectric permittivities of air and water were 1 and 81, respectively. The electric
314
conductivities of air and water were 0 and 0.01 mS/m. Sources were located on the surface every
315
0.05 m from 0.05 to 0.95 m along the y-direction. We use a 1 GHz ricker wavelet to excite the EM
316
wave. The time and spatial spacing were 5 × 10-12 s and 0.005 m, respectively. We set the
317
background media as homogeneous isotropic, anisotropic, and dispersive to simulate clay, finely
318
stratified sand, and moist sand, respectively. The model properties for clay were ε r = 9 and
319
σ = 2 mS / m . The properties for finely stratified sand were ε r xx = ε ryy = 15 , ε rzz = 9 ,
320
σ xx = σ yy = 1 mS/m, and σ zz = 2 mS/m. The properties for moist sand were ε z = 15ε 0 ,
321
ε ∞ = 9ε 0 , τ 0 = 15 × 10−9 s, and σ = 2 mS / m .
21
322 323
Fig. 9.
324
Complex model. An air-filled and a water-filled pipe are buried underground.
325 326
a
b
327
328 329
c
22
330
Fig. 10.
331
Common offset sections of isotropic (a), anisotropic (b), and dispersive (c) background media
332
models.
333
The common offset sections are shown in Fig. 10. The source-receiver offset was 0.05m. The
334
reflections events of two pipes were different in each case, which is caused by the different
335
dielectric permittivities and conductivities of air and water. Other than this, the travel time was
336
evidently different between anisotropic and isotropic cases. This illustrates that the EM wave
337
propagates with different velocities in anisotropic case. Compared with the isotropic case, the EM
338
wave energy was weaker in the dispersive case, which was expected. The above examples
339
illustrate that the proposed method can provide a reasonable solution for 3D anisotropic and
340
dispersive problems.
341
6. Discussion and conclusions
342
In this paper, we describe a new method known as the RSG-FDTD method for modeling EM
343
waves inside anisotropic and dispersive materials. This method defines new directions for
344
approximating spatial partial derivatives. This method may be superior to Yee’s method when
345
simulating EM waves propagating inside anisotropic media because it avoids the interpolation of
346
wave field components. Furthermore, it requires fewer interpolations of media properties than
347
Yee’s method.
348
Using this approach, we generated wave fields and records of electric fields components in
349
homogeneous anisotropic and dispersive media, as well as in a complex model. The results
350
showed good agreement with those obtained with GprMax, which highlights the validity of the 23
351
proposed method. The results of common offset sections reveal that anisotropy and dispersion can
352
apparently influence the interpretation of GPR data. The analysis of the arrival time is one of the
353
easiest and most intuitive ways to evaluate the anisotropy, and it can be used to obtain wave speed
354
expediently. The results reveal that the velocity of an EM wave travelling in anisotropic media
355
depends on the propagation direction.
356
We only studied the TI anisotropy (also known as uniaxial anisotropy) and assumed the
357
conductivities are relatively small. In the future, we plan to study the biaxial anisotropic model
358
( ε xx ≠ ε yy ≠ ε zz , σ xx ≠ σ yy ≠ σ zz ). Furthermore, we will study the effect of non-negligible
359
conductivities. The dispersion modeling will include Lorentz and Drude functions in the future as
360
well. In addition, a better absorbing boundary , such as a time-synchronized CPML (Giannakis
361
and Giannopoulos, 2015), will be adopted. These studies may help us acquire a deeper
362
understanding of underground media.
363
Acknowledgement
364
This research is supported by the National Key R&D Program of China (Grant No.
365
2018YFB0605503) and, the Open Fund of State Key Laboratory of Water Resource Protection
366
and Utilization in Coal Mining (Grant No. SHJT-17-42.16).
367
Computer Code Availability
368
The name of the code we developed for this method is GPR_3D_ANI_RSG_CPML.cpp. It was
369
developed by Dr. Yongxu Lu (one of the authors). His contact address is Room 410, Minzu
370
Building, Ding No.11 Xueyuan Road, Haidian District, Beijing 100083 P. R.China. His telephone 24
371
number is +86 15650716135, and his e-mail is [email protected]. The code was firstly made
372
available at 2018 and can be run on both Linux/Windows systems. A C++ development
373
environment such as G++ in Linux is needed. The code can be run on computers that support
374
Ubuntu 16.04 or Windows 7. However, the memory size limits the size of the model. The program
375
is
376
https://github.com/samhomes/GPR-RSG-ANI.git.
377
Appendix A: Derivation of the stability criterion for RSG
378 379
about
36
KB.
Furthermore,
the
code
382
∂ ( H + jE ) , ∂t
(A.1)
(A.2)
considering the following pair of eigenvalue problems:
V = ΛV ,
(A.3a)
numerical
(A.3b)
To ensure the algorithm’s stability, the imaginary part of Λ must fall within the range (Taflove and Brodwin, 1975) −
385
Github:
The stability of a particular numerical representation of Eq. A.2 can be examined simply by
j∇ numerical × V = ΛV .
384
on
or more simply as:
∂ ∂t
383
available
rewrite Maxwell’s equations as:
j∇ × V = ∂V ∂t , where V = H + jE . 381
is
We let µ = 1, ε = 1, σ = 0 and vEM = 1 as Taflove and Brodwin (1975) proposed. We can
j∇ × ( H + jE ) = 380
source
2 2 ≤ Im ( Λ ) ≤ . ∆t ∆t
We let
25
(A.4)
V I , J , K = V0 e
(
j k%x I ∆x + k% y J ∆y + k%z K ∆z
)
(A.5)
,
386
represent an arbitrary lattice spatial mode. Using the second order space differencing of RSG (Eqs.
387
25) – (26) to implement the derivatives of the curl operator, Eq. (A.3b) yields: xˆ % ∆x % ∆y % ∆z ∆x sin k x 2 cos k y 2 cos k z 2 + yˆ ∆x % ∆y % ∆z −2 cos k%x sin k y cos k z + × V I , J ,K = Λ V I ,J ,K , 2 2 2 ∆y zˆ ∆x % ∆y % ∆z cos cos k%x ky sin k z 2 2 2 ∆z
(A.6)
388
where xˆ , yˆ , and zˆ are unit vectors in the x-, y-, and z-coordinate directions. After performing
389
the cross product and writing out the x, y, and z component equations, we can obtain the
390
following:
1 ∆x 2 % ∆y 2 % sin 2 k%x cos k y cos k z 2 2 2 ( ∆x ) 1 ∆x 2 % ∆y 2 % Λ 2 = −4 cos 2 k%x sin k y cos k z 2 ( ∆y ) 2 2 1 ∆x 2 % ∆y 2 % cos 2 k%x cos k y sin k z 2 2 2 ( ∆z ) 391 392
394
( ∆h = ∆x = ∆y = ∆z ), it is clear that: 2 2 ≤ Im ( Λ ) ≤ . ∆h ∆h
396
(A.8)
To satisfy the stability condition (Eq. (A.4)) for the arbitrary lattice spatial mode, and for arbitrary value of vEM , we set:
∆tvEM ≤1. ∆h 395
(A.7)
Assuming the worst case for the trigonometric functions in Eq. A.7 and for equal grid spacing
− 393
∆z + 2 ∆z + . 2 ∆z 2
(A.9)
For arbitrary order of spatial difference, the stability condition can be obtained in the same manner. It yields: 26
∆tvEM ≤1 ∆h
M
∑C m =1
(M) m
.
(A.10)
397 398
References
399
Alsunaidi, M.A., Al-Jabr, A.A., 2009. A general ADE-FDTD algorithm for the simulation of dispersive
400
structures. IEEE photonics technology letters, 21(12), 817–819.
401
Beggs, J.H., Luebbers R.J., Yee K.S., Kunz K.S., 1992. Finite-difference time-domain implementation of
402
surface impedance boundary-conditions. IEEE Transactions on Antennas and Propagation 40, 49-56.
403
Beker, B., Umashankar, K.R., Taflove, A., 1986. Numerical analysis and validation of the combined field
404
surface integral equations for electromagnetic scattering by arbitrary shaped two-dimensional anisotropic
405
objects. IEEE Trans. Antennas Propagat. 37, 1573-1581.
406
Berenger, J. P., 1994. A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys.
407
144, 185-200.
408
Bergmann, T., Robertsson, J. O., Holliger, K., 1998. Finite-difference modeling of electromagnetic wave
409
propagation in dispersive and attenuating media. Geophysics, 63(3), 856-867.
410
Bradford, J. H., Nichols, J., Harper, J. T., and Meierbachtol, T., 2013. Compressional and EM wave velocity
411
anisotropy in a temperate glacier due to basal crevasses, and implications for water content estimation, Annals
412
of Glaciology 54(64), 168-178.
413
Carcione, J.M., 2007. Wave Propagation in Anisotropic, Anelastic, Porousand Electromagnetic Media, 2nd edn,
414
Elsevier Science & Technology Wave fields in real media, 538 pp.
415
Carcione J.M., Schoenberg M.A., 2000. 3-D ground-penetrating radar simulation and plane-wave theory in
416
anisotropic media. Geophysics, 65(5), 1527-1541. 27
417
Chen, H., Huang, T., 1998. Finite-difference time-domain simulation of GPR data. Journal of Applied
418
Geophysics 40, 139-163.
419
Diamanti, N., Giannopoulos A., 2009. Implementation of ADI-FDTD subgrids in ground penetrating radar
420
FDTD models. Journal of Applied Geophysics, 67, 309-317.
421
Fan, G.X., Liu, Q.H., 2000. An FDTD algorithm with perfectly matched layers for general dispersive media,”
422
IEEE transactions on antennas and propagation, 48(5), 637–646.
423
Fang, H., Lin, G., 2012. Symplectic partitioned Runge–Kutta methods for two-dimensional numerical model of
424
ground penetrating radar. Computers & Geosciences 49, 323-329.
425
Georgakopoulos, S.V., Birtcher C. R., Balanis C.A., et al., 2002. Higher-order finite-difference schemes for
426
electromagnetic radiation, scattering and penetration, part I, Theory. IEEE Antenna’s and Propagation Magazine,
427
44(1), 134-142.
428
Giannakis, I., Giannopoulos, A., 2014. A novel piecewise linear recursive convolution approach for dispersive
429
media using the finite-difference time-domain method. IEEE transactions on antennas and propagation, 62(5),
430
2669-2678.
431
Giannakis, I., Giannopoulos, A., 2015. Time-Synchronised convolutional perfectly matched layer for improved
432
absorbing performance in FDTD. IEEE Antennas and Wireless Propagation Letters, 14, 690-693.
433
Giannakis, I., Giannopoulos, A., Warren, C., 2016. A Realistic FDTD Numerical Modeling Framework of
434
Ground Penetrating Radar for Landmine Detection. IEEE Journal of Selected Topics in Applied Earth
435
Observations and Remote Sensing, 9(1), 37-51.
436
Giannakis, I., Giannopoulos, A., Warren, C., 2019. Realistic FDTD GPR antenna models optimized using a
437
novel linear/nonlinear full-waveform inversion. IEEE transactions on geoscience and remote sensing, 57(3),
28
438
1768-1778.
439
Giannopoulos, A., 2005. Modelling ground penetrating radar by GprMax. Construction and Building Materials
440
19, 755–762.
441
Holliger, K., Bergmann, T., 2002. Numerical modeling of borehole georadar data. Geophysics 67, 1249-1257.
442
Irving, J., Knight, R., 2006. Numerical modeling of ground-penetrating radar in 2-D using MATLAB.
443
Computers & Geosciences 32, 1247-1258.
444
Kelley, D.F., Luebbers, R.J., 1996. Piecewise linear recursive convolution for dispersive media using FDTD.
445
IEEE transactions on antennas and propagation, 44(6), 792–797.
446
Kumlu, D., and Erer, I., 2019. Detection of buried objects in ground penetrating radar data using incremental
447
nonnegative matrix factorization. Remote Sensing Letters 10(7), 649-658.
448
Li, J., Zeng, Z., Huang, L., et al., 2012. GPR simulation based on complex frequency shifted recursive
449
integration PML boundary of 3D high order FDTD. Computers & Geosciences 49, 121-130.
450
Li, Y. F., Matar, O. B., 2010. Convolutional perfectly matched layer for elastic second-order wave equation. J.
451
Acoust. Soc. Am. 127, 1318-1327.
452
Luebbers, R. J., Hunsberger, F., 1992. FDTD for nth-order dispersive media. IEEE transactions on antennas and
453
propagation 40, 1297-1301.
454
Millington T.M., Cassidy N.J., 2010. Optimising GPR modelling: A practical, multi-threaded approach to 3D
455
FDTD numerical modelling. Computers & Geosciences 36, 1135-1144.
456
Ralchenko, M., Svilans, M., Samson, C., et al., 2015. Finite-difference time-domain modelling of
457
through-the-Earth radio signal propagation. Computers & Geosciences 85, 184-195.
458
Rodríguez-Abad, I., Martínez Sala, R. M., García García, F., and Capuz Lladró, R., 2010. Non-destructive
29
459
methodologies for the evaluation of moisture content in sawn timber structures: ground-penetrating radar and
460
ultrasound techniques. Near Surface Geophysics 8(6), 475-482.
461
Ronden, J.A., Gedney, S.D., 2000. Convolutional PML (CPML): An efficient FDTD implementation of the
462
CFS-PML for arbitrary media. Microwave Opt. Technol. Lett. 27, 334-339.
463
Saenger, E.H., Gold, N., Shapiro, S.A., 2000. Modeling the propagation of elastic waves using a modified
464
finite-difference grid. Wave Motion 31, 77-92.
465
Saenger, E.H., Bohlen, T., 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation
466
using the rotated staggered grid. Geophysics 69, 583-591.
467
Schneider, J., Hudson, S., 1993. The finite-difference time-domain method applied to anisotropic material. IEEE
468
Trans. Antennas Propagat. 41, 994–999.
469
Strickel, M. A., Taflove, A., 1990. Time-domain synthesis of broadband absorptive coatings for two-dimensional
470
conducting targets.
471
Taflove, A., 2000. Computational electrodynamics: the finite-difference time-domain method. Boston, MA :
472
Artech House.
473
Taflove, A., Brodwin, M.E., 1975. Numerical solution of steady-state electromagnetic scattering problems using
474
the time-dependent Maxwell’s equations. IEEE transactions on microwave theory and techniques MTT-23,
475
623-630.
476
Teixeira, F. L., Chew, W. C., Straka, M., et al., 1998. Finite-difference time-domain simulation of ground
477
penetrating radar on dispersive, inhomogeneous, and conductive soils. IEEE transactions on geoscience and
478
remote sensing, 36(6), 1928-1937.
479
Ursin, B., 1983. Review of elastic and electromagnetic wave propagation in horizontally layered media,
IEEE Trans. Antennas Propagat. 38, 1084-1091.
30
480
Geophysics, 48(8), 1063–1081.
481
Wang, T., Tripp, A. C., 1996. FDTD simulation of EM wave propagation in 3-D media. Geophysics 61, 110–
482
120.
483
Warren, C., Giannopoulos, A., 2011. Creating finite-difference time-domain models of commercial
484
ground-penetrating radar antennas using Taguchi's optimization method. Geophysics, 76(2), G37-G47.
485
Warren, C., Giannopoulos, A., Giannakis, I., 2016. gprMax: Open source software to simulate electromagnetic
486
wave propagation for Ground Penetrating Radar. Computer Physics Communications, 209, 163-170.
487
Warren, C., Giannopoulos, A., Gray, A., et al., 2019. A CUDA-based GPU engine for gprMax: Open source
488
FDTD electromagnetic simulation software. Computer Physics Communications, 237, 208-218.
489
Yee, K., 1966. Numerical solution of initial boundary value problems involving Maxwell’s equations in
490
isotropic media. IEEE Transactions on Antennas and Propagation AP-14, 302–307.
31
Highlights: 1. Introduced a rotated staggered grid method to simulate the EM waves propagating. 2. The method can be used to simulate EM waves in anisotropic media and dispersive media. 3. Test the algorithm with complex models.
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: