# Seismic Waveform Inversion Using the Finite-Difference Contrast Source Inversion Method

Authors: Bo Han; Qinglong He; Yong Chen; Yixin Dou
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper extends the finite-difference contrast source inversion method
to reconstruct the mass density for two-dimensional elastic wave inversion in the
framework of the full-waveform inversion. The contrast source inversion method
is a nonlinear iterative method that alternatively reconstructs contrast sources and
contrast function. One of the most outstanding advantages of this inversion method
is the highly computational efficiency, since it does not need to simulate a full forward
problem for each inversion iteration. Another attractive feature of the inversion
method is that it is of strong capability in dealing with nonlinear inverse problems
in an inhomogeneous background medium, because a finite-difference operator
is used to represent the differential operator governing the two-dimensional elastic
wave propagation. Additionally, the techniques of a multiplicative regularization
and a sequential multifrequency inversion are employed to enhance the quality of
reconstructions for this inversion method. Numerical reconstruction results show
that the inversion method has an excellent performance for reconstructing the objects
embedded inside a homogeneous or an inhomogeneous background medium.


---

## Body

## 1. Introduction

Full-waveform inversion (FWI) is a powerful tool to get high-resolution parameter models in the context of seismic tomography. This feature of FWI makes it very suitable for estimating subsurface parameters, which is the ultimate goal in exploration seismology. In recent years, with the development of computer technology and the requirement of higher quality of seismic imaging, FWI becomes a very promising technique for exploration seismology as acquisition improves in size and density. The process of FWI is to find an optimal parameter model by minimizing an objective functional that measures the discrepancy between the observed seismograms and the synthetic seismograms. In other words, the specific implementation of FWI is an iterative procedure that is usually based on gradient-based optimization, either in time domain or frequency domain (see [1–9]). In this paper, we only focus on frequency domain case. Simulation in frequency domain has its inherent advantages: it can obtain a relatively good reconstruction result with only a few frequencies, which means that it can significantly reduce the data redundancy; in addition, a parallel computation programming can be easily implemented for different frequency inversions, which results from the independence of each frequency. Most methods in the framework of FWI are based on gradient or Newton method. A common feature of these gradient-based optimization methods is the highly expensive forwarding calculation for estimating the gradient of the object functional with respect to the interested parameter. One special method to estimate the gradient for a minimization problem based on gradient method is the adjoint method (see [10]). Additionally, seismic inverse problems usually are nonlinear and illposed, which make the seismic inversion difficult to handle. In short, FWI incorporated both traveltime and amplitude information is a highly computational process, no matter it is implemented in time domain or frequency domain, and this feature greatly prevents from its widespread applications in exploration seismology.An alternative method is proposed by van den Berg and Kleinman (see [11]), extended by van den Berg and Abubakar to effectively reconstruct the interested parameters in the field of the microwave imaging (see [12]). The method is referred to as the so-called contrast source inversion (CSI) method, which has two versions: IE-CSI based on integral equation formulation (see [11, 12]) and FD-CSI based on the finite-difference operator (see [13]). The IE-CSI method is very efficient for the homogeneous or simple layered background medium inverse problems, since the IE-CSI method needs to calculate Green function of the background medium, while the FD-CSI method has a good performance for solving high complex problems in an inhomogeneous background medium. Another highly efficient scattering inversion method called subspace-based optimization method (SOM) shares this good property of CSI method, which also has two variants adopted in homogeneous and inhomogeneous background medium, respectively (see [14, 15]). The comparisons in detail between SOM and CSI are covered in [14], and this reference also concludes that the CSI is close in spirit to a special case of SOM and the SOM could be convergent for a few special problems. The interesting advantage of CSI method is that it does not require solving any full forward problem in each inversion iteration; therefore, it can significantly reduce the computational cost, and it may be very efficient for large scale or high dimension inverse problems (see [16]). Both the IE-CSI and FD-CSI methods incorporating with multiplicative regularization constraint have been successfully applied in electromagnetic wave imaging and acoustic wave imaging (see [17–21]).In this paper, we extend the FD-CSI method to reconstruct the mass density based on a two-dimensional elastic wave equation (in fact, all theories in this paper are suitable for three-dimensional problems without any changes) in an isotropic medium, either homogeneous or inhomogeneous background medium, and its IE-version is presented in [22], which is successfully employed in elastodynamics when the background medium is homogeneous. For our forward modeling, a perfectly matched layer (PML) absorbing boundary condition is employed to reduce the boundary reflections (see [23, 24]), and, considering the fact that the stiffness matrix usually is large and highly sparse, the techniques of matrix compression store and LU decomposition of a sparse matrix are employed to decrease the storage space and increase the computation efficiency (see [25]). For the inversion processing, since all the finite-difference operators involving the FD-CSI algorithm are only dependent on the background medium and angular frequency that both of them are invariant throughout the inversion process, we use a direct solver such as an LU decomposition method to solve these finite-difference responses to improve the efficiency of the forward modeling. We invert these finite-difference operators only once at the start of the inversion processing, and the results can be reused for multiple source positions and successive iterations of the inversion.In order to get a better reconstruction profile and some geometrical information of an abnormality, a multiplicative regularization factor called FD-MRCSI is incorporated, which is first introduced by van den Berg et al. (see [26]). The multiplicative regularization factor has the same function as the total variation (TV) regularization. The most promising advantage of the multiplicative regularization is that the regularization parameter is automatically calculated during each inversion. Furthermore, in order to mitigate the nonlinearity and enhance the reconstruction results, the sequential multifrequency inversion approach is employed; that is, the seismic inversion starts from low frequency and moves upward to high frequency. To illustrate the performance of the finite-difference contrast source inversion method in seismic inversion problems, we present several numerical experiments, either in a homogeneous or an inhomogeneous background medium.


## 2. Forward Modeling

In this section, we present the two-dimensional forward modeling using finite-difference method, which is of great importance for the inversion method. The forward modeling simulates the propagation process of the elastic wave in an isotropic inhomogeneous medium and drives the inversion algorithm. In time domain, the governing PDE system, called elastodynamic system, is written as the following differential equation system:(1)
ρ
∂
t
2
u
=
∇
·
σ
+
f
,
where u is the displacement vector, ρ is the mass density, σ denotes the stress tensor of the elastic medium, ∇
· represents the divergence operator with respect to spatial variables, and the external force is denoted by the term f. In an isotropic medium, the stress tensor σ is given by the following formulas:
(2)
σ
i
j
=
λ
δ
i
j
ε
k
k
+
2
μ
ε
i
j
,
(3)
ε
i
j
=
1
2
(
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
)
,
where the indices i and j can be 1 or 2, δ
i
j is the well-known Kronecker delta symbol, ε
i
j denotes the strain tensor, and λ and μ are the Lame parameters of the isotropic medium.For simplicity, all physical variables in frequency domain are also denoted by the same symbols as in time domain unless it makes an ambiguity. For example,u also represents the displacement vector in frequency domain. By substituting (2) and (3) into (1) and employing the Fourier transformation on both sides of (1), we get the Navier equations as
(4)
ω
2
ρ
(
r
)
u
(
r
,
ω
,
r
s
)
+
∇
[
λ
(
r
)
∇
·
u
(
r
,
ω
,
r
s
)
]
+
∇
·
[
μ
(
r
)
{
∇
u
(
r
,
ω
,
r
s
)
+
[
∇
u
(
r
,
ω
,
r
s
)
]
T
}
]
=
f
(
r
,
ω
,
r
s
)
,
where r denotes the spatial variable, r
s represents the source position, ω is the angular frequency, the source term f is a Ricker wavelet represented in frequency domain throughout this paper, T is a matrix transpose, and ∇ is a gradient operator with respect to the spatial variable r.System (1) and (4) defines the wave propagation in an infinite medium in time and frequency domains, respectively. However, the wave field is numerically modeled in a finite domain. Therefore, waves will reflect from the computation edges of the domain, and the artificial waves will be recorded. In time domain, these reflections can be eliminated by the time windowing; however, such reflections cannot be eliminated in frequency domain. Therefore, we should attenuate or suppress the unwanted waves to simulate a modeling in a finite domain, and one type of absorbing boundary conditions should be applied. For our forward modeling, a stabilized unsplit convolutional perfectly matched layer (PML) is employed (see [24]), and the source termf we added is in term of p-wave incident excitation.The parameters and wave fields are uniformly discretized on a regular two-dimensional domain at intervals ofΔ
x and Δ
z by using the second-order finite-differences described in [27, 28]. Note that the spatial sampling intervals Δ
x and Δ
z must be small enough so that the numerical solutions of the discretization linear equation system of (4) are sufficiently accurate. Generally, this five-point stencil requires about 10 grid points per shear wave wavelength in the lowest velocity zone for an accurate modeling (see Pratt [28]). After the unknowns at each grid point are ordered by row-major rule, we obtain a linear equation system as
(5)
A
(
ω
)
U
(
ω
)
=
F
(
ω
)
,
where A is the finite-difference matrix, or stiffness matrix, which is asymmetric and extremely sparse, U is the unknowns vector, and F denotes the discrete vector of the source term. To solve the linear equation system (5), we choose a direct solver method based on the sparse matrix LU decomposition rather than an iterative solver. The main reason we employ a direct solver method is that the stiffness matrix A is only dependent on the background medium and an angular frequency ω. Hence, we can get all solutions for all source excitations simultaneously by factoring the finite-difference matrix A only once. However, a new matrix decomposition is required for another frequency.


## 3. Finite-Difference Contrast Source Inversion Algorithm

### 3.1. Problem Formulation

We consider the two-dimensional elastic wave inversion problem in a homogeneous or an inhomogeneous background medium. Throughout the paper, the inversion domain or object domain in which the scatterers are embedded is denoted by symbolD, and the data domain S is the surface where the sources and receivers are located. We assume that both of the domains are contained within the total domain D
1 called finite-difference computation domain. By this configuration of this inverse problem, we can slightly reduce the inversion domain. Assume that there are N
s sources successively illuminating the inversion domain and the index for the jth source is denoted by subscript j. For each source j, the total field vector u
j
(
r
,
ω
) and incident field vector u
j
inc
(
r
,
ω
) satisfy the following Navier equations, respectively:
(6)
N
[
u
j
(
r
,
ω
)
]
≡
ω
2
ρ
(
r
)
u
j
(
r
,
ω
,
r
s
)
+
∇
[
λ
(
r
)
∇
·
u
j
(
r
,
ω
,
r
s
)
]
+
∇
·
[
μ
(
r
)
{
∇
u
j
(
r
,
ω
,
r
s
)
+
[
∇
u
j
(
r
,
ω
,
r
s
)
]
T
}
]
=
f
(
r
,
ω
,
r
s
)
,
r
∈
D
1
,
(7)
N
b
[
u
j
inc
(
r
,
ω
)
]
≡
ω
2
ρ
b
(
r
)
u
j
inc
(
r
,
ω
,
r
s
)
+
∇
[
λ
(
r
)
∇
·
u
j
inc
(
r
,
ω
,
r
s
)
]
+
∇
·
[
μ
(
r
)
{
∇
u
j
inc
(
r
,
ω
,
r
s
)
+
[
∇
u
j
inc
(
r
,
ω
,
r
s
)
]
T
}
]
=
f
(
r
,
ω
,
r
s
)
,
r
∈
D
1
.
In (7), the parameter involving the subscript b denotes the relative parameters when the equation is employed to describe the seismic wave propagation in a known background medium. Subtracting (7) from (6) and using the relationship between the total filed u
j, incident filed u
j
inc, and the scattered filed u
j
sct, the scattered fields satisfy
(8)
N
b
[
u
j
sct
(
r
,
ω
)
]
=
-
ω
2
ρ
b
(
r
)
w
j
(
r
,
ω
)
,
r
∈
D
1
,
where w
j
(
r
,
ω
) are the contrast source vectors, defined as
(9)
w
j
(
r
,
ω
)
=
χ
(
r
)
u
j
(
r
,
ω
)
,
r
∈
D
,
in which the contrast function χ
(
r
) is given by
(10)
χ
(
r
)
=
ρ
(
r
)
-
ρ
b
(
r
)
ρ
b
(
r
)
=
ρ
(
r
)
ρ
b
(
r
)
-
1
,
r
∈
D
.According to (7), the differential operator N
b with the known Lame parameters λ and μ only depends on the background mass density ρ
b
(
r
). The inverse of the operator N
b is formally given by the linear operator L
b
[
·
]
=
N
b
-
1
[
-
ω
2
ρ
b
(
r
)
(
·
)
]. Therefore, the solution vectors of (8) can be symbolically denoted by the following equation:
(11)
u
j
sct
(
r
,
ω
)
=
N
b
-
1
[
-
ω
2
ρ
b
(
r
)
w
j
(
r
,
ω
)
]
=
L
b
[
w
j
(
r
,
ω
)
]
,
r
∈
D
1
.
By using (11), the measured scattered field data vectors at the receiver r
R location can be formally presented by the following equation:
(12)
d
j
(
r
R
,
ω
)
=
M
S
L
b
{
[
w
j
(
r
,
ω
)
]
}
,
r
R
∈
S
,
r
∈
D
1
,
and the object equation or state equation satisfies
(13)
w
j
(
r
,
ω
)
=
χ
(
r
)
u
j
inc
(
r
,
ω
)
+
χ
(
r
)
M
D
{
L
b
[
w
j
(
r
′
,
ω
)
]
}
,
r
∈
D
,
r
′
∈
D
1
,
where M
S is an operator selecting the scattered fields on the measurement surface S among the scattered fields in the total computation domain, and the operator M
D selects the fields inside the inversion domain D.Note that (12) and (13) are the fundamental equations for the CSI method, which are used to establish the cost functional. From (12) and (13), the data error r
j
S
(
r
,
ω
) and the object (or state) error r
j
D
(
r
,
ω
) are denoted by the following equations, respectively:
(14)
r
j
S
(
r
,
ω
)
=
d
j
(
r
R
,
ω
)
-
M
S
L
b
{
[
w
j
(
r
,
ω
)
]
}
,
r
R
∈
S
,
r
∈
D
1
,
(15)
r
j
D
(
r
,
ω
)
=
χ
(
r
)
u
j
inc
(
r
,
ω
)
+
χ
(
r
)
M
D
{
L
b
[
w
j
(
r
′
,
ω
)
]
}
-
w
j
(
r
,
ω
)
,
r
∈
D
,
r
′
∈
D
1
.To establish the objective functional, we should firstly give the specific expressions of the inner product andL
2-norm as
(16)
〈
u
,
v
〉
Ω
=
∫
Ω
‍
u
(
r
)
·
v
(
r
)
¯
d
r
,
(17)
∥
r
∥
Ω
2
=
∫
Ω
‍
u
(
r
)
·
u
(
r
)
¯
d
r
,
where the overbar denotes the complex conjugate for each component of a complex vector and symbol · represents the product of two vectors. Ω is the integral domain, which can be either S or D.For simplicity, the explicit dependence of the field quantities onr and r
′ is dropped in the remainder of this paper.


### 3.2.  Finite-Difference Contrast Source Inversion Method

In the framework of the contrast source inversion method (see [11]), it does not eliminate the contrast source quantities w
j in (12) and (13). On the contrary, it alternatively reconstructs contrast sources w
j and contrast function χ by minimizing a cost functional. In order to avoid minimizing a nonquadratic objective function caused by the presence of contrast function χ in both numerators and denominators, we borrow the idea from [14] and slightly modify the cost functional of CSI method. Therefore, the modified cost functional to be minimized at the nth step is defined as the following form:
(18)
F
n
(
χ
,
w
j
)
=
∑
j
∥
r
j
S
∥
S
2
‍
∑
j
∥
d
j
S
∥
S
2
‍
+
∑
j
∥
r
j
D
∥
D
2
‍
∑
j
∥
χ
n
-
1
u
j
inc
∥
D
2
‍
=
F
S
(
w
j
)
+
F
n
D
(
χ
,
w
j
)
,
where F
S
(
w
j
) represents the data equation error and F
n
D
(
χ
,
w
j
) represents the object equation error. The scattered field data is obtained by subtracting the simulated background field from the simulated model field. From (14), (15), and (18), we note that the background medium does not change throughout the inversion process; therefore, the computation cost of the inversion of the finite-difference operator N
b is relatively cheap for each frequency, which results from the LU decomposition required only once.The main idea of the CSI method is to reconstruct two interlaced sequences,w
j
,
n and χ
n, by minimizing the objective functional (18). The first step of FD-CSI method is to update contrast source vectors w
j
,
n by the conjugate-gradient method with Polak-Ribiere search directions. While updating the contrast source, the contrast function χ is assumed to be χ
n
-
1. At the nth step, the contrast sources w
j
,
n are updated by the following formula:
(19)
w
j
,
n
=
w
j
,
n
-
1
+
α
j
,
n
ω
v
j
,
n
,
where α
j
,
n
ω is an update-step size and v
j
,
n is the Polak-Ribiere search directions (see [29]) given by
(20)
v
j
,
0
=
0
,
v
j
,
n
=
g
j
,
n
ω
+
∑
k
〈
g
j
,
n
ω
,
g
j
,
n
ω
-
g
j
,
n
-
1
ω
〉
∑
k
∥
g
j
,
n
-
1
ω
∥
D
2
v
j
,
n
-
1
,
in which g
j
,
n
ω is the gradient of the cost function with respect to w
j
,
n evaluated at the (
n
-
1
)th iteration. The details of this type of the conjugate-gradient can be found in [30]. After updating the gradient g
j
,
n
ω, the update-steps α
j
,
n
ω can be obtained by minimizing the cost functional F
n; that is,
(21)
α
j
,
n
ω
=
arg
min
⁡
α
{
F
n
(
χ
n
-
1
,
w
j
,
n
-
1
+
α
v
j
,
n
)
}
,
and the explicit expression of the minimizer α
j
,
n
ω of (21) can be founded by setting the derivative with respect to α being equal to zero; that is,
(22)
α
j
,
n
ω
=
(
-
〈
g
j
,
n
ω
,
v
j
,
n
〉
D
)
×
(
η
S
∥
M
S
{
L
b
[
v
j
,
n
]
}
∥
S
2
000
+
η
n
D
∥
v
j
,
n
-
χ
n
-
1
M
D
{
L
b
[
v
j
,
n
]
}
∥
D
2
)
-
1
,
where the normalization factors η
S and η
n
D are
(23)
η
S
=
(
∑
j
∥
d
j
∥
S
2
)
-
1
,
η
n
D
=
(
∑
j
∥
χ
n
-
1
u
j
inc
∥
D
2
)
-
1
.The second step of the CSI method is to update contrast functionχ. Once we update the contrast source vectors w
j with conjugate gradient method, the contrast function χ is updated by finding the minimizer of the functional F
n
D
(
χ
,
w
j
). We assume the contrast source vectors w
j are constant during this process. It can be shown that the updating formula has the closed form as follows:
(24)
χ
n
=
∑
j
Re
(
w
j
,
n
·
u
¯
j
,
n
)
∑
j
|
u
j
,
n
|
2
,
where
(25)
u
j
,
n
=
u
j
,
n
inc
+
M
D
{
L
b
[
w
j
,
n
]
}
.
In (24), the symbol Re indicates getting the real part of a complex number. The details of finding the minimizer with respect to χ of the objective functional (18) can be found in [11].To start the FD-CSI algorithm, Considering the nonlinearity and illposedness of the inverse problem, and the local convergence of the FD-CSI method, a good initial guess is very important to get a better reconstruction result. We initialize the algorithm by the following way:(26)
w
j
,
0
=
∥
L
b
*
{
M
S
*
[
d
j
]
}
∥
D
2
∥
M
S
(
L
b
[
L
b
*
{
M
S
*
[
d
j
]
}
]
)
∥
S
2
L
b
*
{
M
S
*
[
d
j
]
}
,
and the contrast function is calculated by using (24) and (25). In (26), L
b
* is the adjoint operator of the operator L
b, which is defined as
(27)
L
b
*
[
·
]
=
(
N
b
*
)
-
1
[
-
ω
2
ρ
b
(
r
)
¯
(
·
)
]
,
where
(28)
N
b
*
[
·
]
≡
ω
2
ρ
(
r
)
¯
(
·
)
+
∇
[
λ
(
r
)
¯
∇
·
(
·
)
]
+
∇
·
[
μ
(
r
)
¯
{
∇
(
·
)
+
[
∇
(
·
)
]
T
}
]
,
and the symbol M
S
* has the same meaning with L
b
*, which maps the field values on the measurement surface into the total computation domain. This initialization method is called back-propagation method, and it has been proved that it can provide a good initial guess (see [8]).To improve the quality of reconstruction results, the multiplicative regularization technique introduced by van den Berg et al. (see [26]) is applied to this inverse problem, which has the same effect as the total variational regularized method. After incorporating the multiplicative regularization factor, the cost functional is rewritten as
(29)
C
n
(
χ
,
w
j
)
=
[
F
S
(
w
j
)
+
F
n
D
(
χ
,
w
j
)
]
F
n
R
(
χ
)
,
where the regularization cost function F
n
R
(
χ
) for the weighted L
2-norm regularization is defined as
(30)
F
n
R
(
χ
)
=
∫
D
‍
b
n
2
(
|
∇
χ
(
r
)
|
2
+
δ
n
2
)
d
r
,
where
(31)
b
n
2
(
r
)
=
1
A
(
|
∇
χ
n
-
1
(
r
)
|
2
+
δ
n
2
)
,
δ
n
2
=
F
n
D
(
χ
n
-
1
,
w
j
,
n
-
1
)
Δ
x
Δ
z
,
where A is the area of the inversion domain and Δ
x
Δ
z is the cell area. More details on finite-difference multiplicative regularization contrast source inversion (FD-MRCSI) method can be found in Appendix A of [8].It is very similar to minimizing functional (29) using conjugate-gradient method. There is no change while updating the contrast source vectors w
j, since the regularization factor F
n
R
(
χ
) is only dependent on contrast function χ and F
n
R
(
χ
n
-
1
)
≡
1. Nevertheless, the contrast function χ is updated by using the conjugate-gradient method with Polak-Ribiere search directions.Compared with an additive regularization method, it is not necessary to determine the artificial weighting regularization parameter because it is automatically calculated at each iteration. As we know that selecting the regularization parameter is a time-consuming and difficult work, especially there is no prior information about the measurement data, multiplicative regularization method is very suitable to invert experimental data as shown in [20].


### 3.3.  Computational Complexity Analysis

The highest computational cost part of FD-MRCSI is to calculate these operatorsL
b and L
b
* to estimate the conjugate g
j and update-steps α
j
ω. Hence, an efficient method to compute the action of these operators is important to increase the efficiency of the inversion algorithm. Considering the fact that these operators only depend on the background medium, which does not change throughout the inversion process for a given frequency ω, we choose an efficient direct method, that is, the LU decomposition, and the decomposition for each of the two discretization operators L
b and L
b
* is done only once at the beginning of the inversion process and then reused at each subsequent step. The decomposition process takes O
(
n
1.5
) operations, where n is the number of the unknowns of the discrete system (n
=
2
N, N is the number of the grid cells) (see [25]). The back-substitution needs O
(
n
log
⁡
n
) operations for each source position. In the framework of FD-MRCSI, back-substitution is required twice for each of the two operators L
b and L
b
* per iteration. Hence, the computational complexity of FD-MRCSI for each frequency is approximately given by
(32)
T
(
n
,
N
iter
,
N
s
)
≈
O
(
2
n
1.5
)
+
N
iter
O
(
2
N
s
n
log
⁡
n
)
,
where N
iter and N
s are the number of inversion iteration and the number of sources, respectively.For usual inversion methods (UIM), the objective functional involving the scattered observed datad
j and scattered synthetic data u
j
sct is given in the following formula:
(33)
F
(
χ
)
=
∑
j
‍
∥
d
j
-
u
j
sct
(
χ
)
∥
S
2
+
α
R
(
χ
)
.
The first term of (33) ensures that the minimizer χ
α of F will be an optimal parameter model, while the second is the regularized term, which guarantees the stabilization of the problem and forces the minimizer χ
α of (33) to satisfy certain regularity properties, and α denotes the regularization parameter that balances the data error and regularized term. If functional (33) is also minimized by using a nonlinear conjugate gradient method, then the computational complexity of this algorithm for this inverse problem can be approximately estimated in analogy to FD-MRCSI method as follows:
(34)
T
1
(
n
,
N
iter
,
N
s
)
≈
N
iter
{
O
(
2
n
1.5
)
+
O
[
(
2
N
S
)
n
log
⁡
n
]
}
,
in which O
(
n
1.5
) represents the LU decomposition computational cost for forward modeling and O
(
n
log
⁡
n
) denotes the back-substitution computation for each source during each iteration process. Table 1 indicates the details of computational complexity for FD-MRCSI and usual inversion methods.Table 1
Computational complexity for FD-MRCSI and UIM.
Method
LU decomposition
Back-substitution
Total cost
UIM
N
iter
O
(
2
n
1.5
)
N
iter
O
(
2
N
S
n
log
⁡
n
)
N
iter
{
O
(
2
n
1.5
)
+
O
(
2
N
S
n
log
⁡
n
)
}
FD-MRCSI
O
(
2
n
1.5
)
N
iter
O
(
2
N
S
n
log
⁡
n
)
O
(
2
n
1.5
)
+
N
iter
O
(
2
N
S
n
log
⁡
n
)
N
iter is the iteration number, N
S presents the number of source, and n denotes the number of the grid cells.From Table1, we observe that the FD-MRCSI method does not need to do an LU decomposition for each iteration, while usual methods require it, and this is the highest computational cost part. Therefore, the FD-MRCSI method is a greatly efficient method.


## 3.1. Problem Formulation

We consider the two-dimensional elastic wave inversion problem in a homogeneous or an inhomogeneous background medium. Throughout the paper, the inversion domain or object domain in which the scatterers are embedded is denoted by symbolD, and the data domain S is the surface where the sources and receivers are located. We assume that both of the domains are contained within the total domain D
1 called finite-difference computation domain. By this configuration of this inverse problem, we can slightly reduce the inversion domain. Assume that there are N
s sources successively illuminating the inversion domain and the index for the jth source is denoted by subscript j. For each source j, the total field vector u
j
(
r
,
ω
) and incident field vector u
j
inc
(
r
,
ω
) satisfy the following Navier equations, respectively:
(6)
N
[
u
j
(
r
,
ω
)
]
≡
ω
2
ρ
(
r
)
u
j
(
r
,
ω
,
r
s
)
+
∇
[
λ
(
r
)
∇
·
u
j
(
r
,
ω
,
r
s
)
]
+
∇
·
[
μ
(
r
)
{
∇
u
j
(
r
,
ω
,
r
s
)
+
[
∇
u
j
(
r
,
ω
,
r
s
)
]
T
}
]
=
f
(
r
,
ω
,
r
s
)
,
r
∈
D
1
,
(7)
N
b
[
u
j
inc
(
r
,
ω
)
]
≡
ω
2
ρ
b
(
r
)
u
j
inc
(
r
,
ω
,
r
s
)
+
∇
[
λ
(
r
)
∇
·
u
j
inc
(
r
,
ω
,
r
s
)
]
+
∇
·
[
μ
(
r
)
{
∇
u
j
inc
(
r
,
ω
,
r
s
)
+
[
∇
u
j
inc
(
r
,
ω
,
r
s
)
]
T
}
]
=
f
(
r
,
ω
,
r
s
)
,
r
∈
D
1
.
In (7), the parameter involving the subscript b denotes the relative parameters when the equation is employed to describe the seismic wave propagation in a known background medium. Subtracting (7) from (6) and using the relationship between the total filed u
j, incident filed u
j
inc, and the scattered filed u
j
sct, the scattered fields satisfy
(8)
N
b
[
u
j
sct
(
r
,
ω
)
]
=
-
ω
2
ρ
b
(
r
)
w
j
(
r
,
ω
)
,
r
∈
D
1
,
where w
j
(
r
,
ω
) are the contrast source vectors, defined as
(9)
w
j
(
r
,
ω
)
=
χ
(
r
)
u
j
(
r
,
ω
)
,
r
∈
D
,
in which the contrast function χ
(
r
) is given by
(10)
χ
(
r
)
=
ρ
(
r
)
-
ρ
b
(
r
)
ρ
b
(
r
)
=
ρ
(
r
)
ρ
b
(
r
)
-
1
,
r
∈
D
.According to (7), the differential operator N
b with the known Lame parameters λ and μ only depends on the background mass density ρ
b
(
r
). The inverse of the operator N
b is formally given by the linear operator L
b
[
·
]
=
N
b
-
1
[
-
ω
2
ρ
b
(
r
)
(
·
)
]. Therefore, the solution vectors of (8) can be symbolically denoted by the following equation:
(11)
u
j
sct
(
r
,
ω
)
=
N
b
-
1
[
-
ω
2
ρ
b
(
r
)
w
j
(
r
,
ω
)
]
=
L
b
[
w
j
(
r
,
ω
)
]
,
r
∈
D
1
.
By using (11), the measured scattered field data vectors at the receiver r
R location can be formally presented by the following equation:
(12)
d
j
(
r
R
,
ω
)
=
M
S
L
b
{
[
w
j
(
r
,
ω
)
]
}
,
r
R
∈
S
,
r
∈
D
1
,
and the object equation or state equation satisfies
(13)
w
j
(
r
,
ω
)
=
χ
(
r
)
u
j
inc
(
r
,
ω
)
+
χ
(
r
)
M
D
{
L
b
[
w
j
(
r
′
,
ω
)
]
}
,
r
∈
D
,
r
′
∈
D
1
,
where M
S is an operator selecting the scattered fields on the measurement surface S among the scattered fields in the total computation domain, and the operator M
D selects the fields inside the inversion domain D.Note that (12) and (13) are the fundamental equations for the CSI method, which are used to establish the cost functional. From (12) and (13), the data error r
j
S
(
r
,
ω
) and the object (or state) error r
j
D
(
r
,
ω
) are denoted by the following equations, respectively:
(14)
r
j
S
(
r
,
ω
)
=
d
j
(
r
R
,
ω
)
-
M
S
L
b
{
[
w
j
(
r
,
ω
)
]
}
,
r
R
∈
S
,
r
∈
D
1
,
(15)
r
j
D
(
r
,
ω
)
=
χ
(
r
)
u
j
inc
(
r
,
ω
)
+
χ
(
r
)
M
D
{
L
b
[
w
j
(
r
′
,
ω
)
]
}
-
w
j
(
r
,
ω
)
,
r
∈
D
,
r
′
∈
D
1
.To establish the objective functional, we should firstly give the specific expressions of the inner product andL
2-norm as
(16)
〈
u
,
v
〉
Ω
=
∫
Ω
‍
u
(
r
)
·
v
(
r
)
¯
d
r
,
(17)
∥
r
∥
Ω
2
=
∫
Ω
‍
u
(
r
)
·
u
(
r
)
¯
d
r
,
where the overbar denotes the complex conjugate for each component of a complex vector and symbol · represents the product of two vectors. Ω is the integral domain, which can be either S or D.For simplicity, the explicit dependence of the field quantities onr and r
′ is dropped in the remainder of this paper.


## 3.2.  Finite-Difference Contrast Source Inversion Method

In the framework of the contrast source inversion method (see [11]), it does not eliminate the contrast source quantities w
j in (12) and (13). On the contrary, it alternatively reconstructs contrast sources w
j and contrast function χ by minimizing a cost functional. In order to avoid minimizing a nonquadratic objective function caused by the presence of contrast function χ in both numerators and denominators, we borrow the idea from [14] and slightly modify the cost functional of CSI method. Therefore, the modified cost functional to be minimized at the nth step is defined as the following form:
(18)
F
n
(
χ
,
w
j
)
=
∑
j
∥
r
j
S
∥
S
2
‍
∑
j
∥
d
j
S
∥
S
2
‍
+
∑
j
∥
r
j
D
∥
D
2
‍
∑
j
∥
χ
n
-
1
u
j
inc
∥
D
2
‍
=
F
S
(
w
j
)
+
F
n
D
(
χ
,
w
j
)
,
where F
S
(
w
j
) represents the data equation error and F
n
D
(
χ
,
w
j
) represents the object equation error. The scattered field data is obtained by subtracting the simulated background field from the simulated model field. From (14), (15), and (18), we note that the background medium does not change throughout the inversion process; therefore, the computation cost of the inversion of the finite-difference operator N
b is relatively cheap for each frequency, which results from the LU decomposition required only once.The main idea of the CSI method is to reconstruct two interlaced sequences,w
j
,
n and χ
n, by minimizing the objective functional (18). The first step of FD-CSI method is to update contrast source vectors w
j
,
n by the conjugate-gradient method with Polak-Ribiere search directions. While updating the contrast source, the contrast function χ is assumed to be χ
n
-
1. At the nth step, the contrast sources w
j
,
n are updated by the following formula:
(19)
w
j
,
n
=
w
j
,
n
-
1
+
α
j
,
n
ω
v
j
,
n
,
where α
j
,
n
ω is an update-step size and v
j
,
n is the Polak-Ribiere search directions (see [29]) given by
(20)
v
j
,
0
=
0
,
v
j
,
n
=
g
j
,
n
ω
+
∑
k
〈
g
j
,
n
ω
,
g
j
,
n
ω
-
g
j
,
n
-
1
ω
〉
∑
k
∥
g
j
,
n
-
1
ω
∥
D
2
v
j
,
n
-
1
,
in which g
j
,
n
ω is the gradient of the cost function with respect to w
j
,
n evaluated at the (
n
-
1
)th iteration. The details of this type of the conjugate-gradient can be found in [30]. After updating the gradient g
j
,
n
ω, the update-steps α
j
,
n
ω can be obtained by minimizing the cost functional F
n; that is,
(21)
α
j
,
n
ω
=
arg
min
⁡
α
{
F
n
(
χ
n
-
1
,
w
j
,
n
-
1
+
α
v
j
,
n
)
}
,
and the explicit expression of the minimizer α
j
,
n
ω of (21) can be founded by setting the derivative with respect to α being equal to zero; that is,
(22)
α
j
,
n
ω
=
(
-
〈
g
j
,
n
ω
,
v
j
,
n
〉
D
)
×
(
η
S
∥
M
S
{
L
b
[
v
j
,
n
]
}
∥
S
2
000
+
η
n
D
∥
v
j
,
n
-
χ
n
-
1
M
D
{
L
b
[
v
j
,
n
]
}
∥
D
2
)
-
1
,
where the normalization factors η
S and η
n
D are
(23)
η
S
=
(
∑
j
∥
d
j
∥
S
2
)
-
1
,
η
n
D
=
(
∑
j
∥
χ
n
-
1
u
j
inc
∥
D
2
)
-
1
.The second step of the CSI method is to update contrast functionχ. Once we update the contrast source vectors w
j with conjugate gradient method, the contrast function χ is updated by finding the minimizer of the functional F
n
D
(
χ
,
w
j
). We assume the contrast source vectors w
j are constant during this process. It can be shown that the updating formula has the closed form as follows:
(24)
χ
n
=
∑
j
Re
(
w
j
,
n
·
u
¯
j
,
n
)
∑
j
|
u
j
,
n
|
2
,
where
(25)
u
j
,
n
=
u
j
,
n
inc
+
M
D
{
L
b
[
w
j
,
n
]
}
.
In (24), the symbol Re indicates getting the real part of a complex number. The details of finding the minimizer with respect to χ of the objective functional (18) can be found in [11].To start the FD-CSI algorithm, Considering the nonlinearity and illposedness of the inverse problem, and the local convergence of the FD-CSI method, a good initial guess is very important to get a better reconstruction result. We initialize the algorithm by the following way:(26)
w
j
,
0
=
∥
L
b
*
{
M
S
*
[
d
j
]
}
∥
D
2
∥
M
S
(
L
b
[
L
b
*
{
M
S
*
[
d
j
]
}
]
)
∥
S
2
L
b
*
{
M
S
*
[
d
j
]
}
,
and the contrast function is calculated by using (24) and (25). In (26), L
b
* is the adjoint operator of the operator L
b, which is defined as
(27)
L
b
*
[
·
]
=
(
N
b
*
)
-
1
[
-
ω
2
ρ
b
(
r
)
¯
(
·
)
]
,
where
(28)
N
b
*
[
·
]
≡
ω
2
ρ
(
r
)
¯
(
·
)
+
∇
[
λ
(
r
)
¯
∇
·
(
·
)
]
+
∇
·
[
μ
(
r
)
¯
{
∇
(
·
)
+
[
∇
(
·
)
]
T
}
]
,
and the symbol M
S
* has the same meaning with L
b
*, which maps the field values on the measurement surface into the total computation domain. This initialization method is called back-propagation method, and it has been proved that it can provide a good initial guess (see [8]).To improve the quality of reconstruction results, the multiplicative regularization technique introduced by van den Berg et al. (see [26]) is applied to this inverse problem, which has the same effect as the total variational regularized method. After incorporating the multiplicative regularization factor, the cost functional is rewritten as
(29)
C
n
(
χ
,
w
j
)
=
[
F
S
(
w
j
)
+
F
n
D
(
χ
,
w
j
)
]
F
n
R
(
χ
)
,
where the regularization cost function F
n
R
(
χ
) for the weighted L
2-norm regularization is defined as
(30)
F
n
R
(
χ
)
=
∫
D
‍
b
n
2
(
|
∇
χ
(
r
)
|
2
+
δ
n
2
)
d
r
,
where
(31)
b
n
2
(
r
)
=
1
A
(
|
∇
χ
n
-
1
(
r
)
|
2
+
δ
n
2
)
,
δ
n
2
=
F
n
D
(
χ
n
-
1
,
w
j
,
n
-
1
)
Δ
x
Δ
z
,
where A is the area of the inversion domain and Δ
x
Δ
z is the cell area. More details on finite-difference multiplicative regularization contrast source inversion (FD-MRCSI) method can be found in Appendix A of [8].It is very similar to minimizing functional (29) using conjugate-gradient method. There is no change while updating the contrast source vectors w
j, since the regularization factor F
n
R
(
χ
) is only dependent on contrast function χ and F
n
R
(
χ
n
-
1
)
≡
1. Nevertheless, the contrast function χ is updated by using the conjugate-gradient method with Polak-Ribiere search directions.Compared with an additive regularization method, it is not necessary to determine the artificial weighting regularization parameter because it is automatically calculated at each iteration. As we know that selecting the regularization parameter is a time-consuming and difficult work, especially there is no prior information about the measurement data, multiplicative regularization method is very suitable to invert experimental data as shown in [20].


## 3.3.  Computational Complexity Analysis

The highest computational cost part of FD-MRCSI is to calculate these operatorsL
b and L
b
* to estimate the conjugate g
j and update-steps α
j
ω. Hence, an efficient method to compute the action of these operators is important to increase the efficiency of the inversion algorithm. Considering the fact that these operators only depend on the background medium, which does not change throughout the inversion process for a given frequency ω, we choose an efficient direct method, that is, the LU decomposition, and the decomposition for each of the two discretization operators L
b and L
b
* is done only once at the beginning of the inversion process and then reused at each subsequent step. The decomposition process takes O
(
n
1.5
) operations, where n is the number of the unknowns of the discrete system (n
=
2
N, N is the number of the grid cells) (see [25]). The back-substitution needs O
(
n
log
⁡
n
) operations for each source position. In the framework of FD-MRCSI, back-substitution is required twice for each of the two operators L
b and L
b
* per iteration. Hence, the computational complexity of FD-MRCSI for each frequency is approximately given by
(32)
T
(
n
,
N
iter
,
N
s
)
≈
O
(
2
n
1.5
)
+
N
iter
O
(
2
N
s
n
log
⁡
n
)
,
where N
iter and N
s are the number of inversion iteration and the number of sources, respectively.For usual inversion methods (UIM), the objective functional involving the scattered observed datad
j and scattered synthetic data u
j
sct is given in the following formula:
(33)
F
(
χ
)
=
∑
j
‍
∥
d
j
-
u
j
sct
(
χ
)
∥
S
2
+
α
R
(
χ
)
.
The first term of (33) ensures that the minimizer χ
α of F will be an optimal parameter model, while the second is the regularized term, which guarantees the stabilization of the problem and forces the minimizer χ
α of (33) to satisfy certain regularity properties, and α denotes the regularization parameter that balances the data error and regularized term. If functional (33) is also minimized by using a nonlinear conjugate gradient method, then the computational complexity of this algorithm for this inverse problem can be approximately estimated in analogy to FD-MRCSI method as follows:
(34)
T
1
(
n
,
N
iter
,
N
s
)
≈
N
iter
{
O
(
2
n
1.5
)
+
O
[
(
2
N
S
)
n
log
⁡
n
]
}
,
in which O
(
n
1.5
) represents the LU decomposition computational cost for forward modeling and O
(
n
log
⁡
n
) denotes the back-substitution computation for each source during each iteration process. Table 1 indicates the details of computational complexity for FD-MRCSI and usual inversion methods.Table 1
Computational complexity for FD-MRCSI and UIM.
Method
LU decomposition
Back-substitution
Total cost
UIM
N
iter
O
(
2
n
1.5
)
N
iter
O
(
2
N
S
n
log
⁡
n
)
N
iter
{
O
(
2
n
1.5
)
+
O
(
2
N
S
n
log
⁡
n
)
}
FD-MRCSI
O
(
2
n
1.5
)
N
iter
O
(
2
N
S
n
log
⁡
n
)
O
(
2
n
1.5
)
+
N
iter
O
(
2
N
S
n
log
⁡
n
)
N
iter is the iteration number, N
S presents the number of source, and n denotes the number of the grid cells.From Table1, we observe that the FD-MRCSI method does not need to do an LU decomposition for each iteration, while usual methods require it, and this is the highest computational cost part. Therefore, the FD-MRCSI method is a greatly efficient method.


## 4. Numerical Results

In this section, we present some numerical examples based on the FD-MRCSI algorithm. All inversion domains of these tests are inside a square with side lengthL. In order to measure the discrepancy between the exact profile and reconstructed profile, the relative error with respect to contrast quantity χ is introduced, which is given by the following formula:
(35)
Err
=
∥
χ
n
(
r
)
-
χ
exact
(
r
)
∥
D
∥
χ
exact
(
r
)
∥
D
×
100
%
,
where χ
n is the reconstructed result after n iterations and χ
exact is the exact solution in which we are interested.Since the “inverse crime” phenomenon may happen, the synthetic data are generated by using a finite-difference time-domain forward method. A Ricker waveletf with a dominant frequency 30 Hz is employed for generating the synthetic data. However, the inversion algorithm is carried out in frequency domain, and hence the datasets generated by using time-domain forward method are transformed into frequency domain using a fast Fourier transform (FFT) (see [31]). After that, we select some frequencies for the FD-MRCSI algorithm. Once we get the synthetic data, a random noise is added to each frequency synthetic data using this formula as follows:
(36)
d
j
noise
=
d
j
+
max
⁡
j
(
|
d
j
|
)
η
2
(
ξ
Re
I
+
i
ξ
Im
⁡
I
)
,
where d
j
noise and d
j are the field data vectors with noise and the noiseless data, respectively, ξ
Re and ξ
Im
⁡ are uniformly distributed random numbers on the interval (
-
1,1
), I is a two- dimensional vector whose components are equal to one, and η is the desired fraction of the noise.


### 4.1. Reconstruction in a Homogeneous Background Medium

The first example is concerned with recovering a square scatterer from noise observed datad
j
noise in a homogeneous background medium, and the exact model is given by Figure 1(a). The inversion domain D is set as a square with side length 6 km, which does not include PML thickness. For this square scatterer, the second Lame parameter μ is equal to 2.0 GPa, and a Poisson’s ratio is 0.25. The square with the side length 1.5 km is located in the inversion domain, and its centre is at (3.3) km; the mass density of the scatterer is 1.8 g/cm3, 1.5 g/cm3 outside the scatterer. Hence, the contrast χ for the scatterer is 0.2. During the whole inversion procedure, the angular frequency ω is 8 Hz.Applying the FD-MRCSI method to one abnormity inversion. (a) True model. (b) Reconstruction result after 1000 inversion iterations. (c) Discrepancy between true model and reconstruction result. (d) Slices at horizontal distances of 3 km for the true model and reconstruction result.
(a)
True model
(b)
Reconstruction result
(c)
Discrepancy
(d)
Model sliceThe inversion domainD is divided into N
×
N uniform grids using the second-order scheme, adding the n
1 grids of PML; so, the size of the total computation domain is (
N
+
2
n
1
)
×
(
N
+
2
n
1
). For this experiment, N
=
61 and n
1
=
15; therefore, the grid dimension of the total computation domain is 91
×
91. Both the sources and receivers are assigned to 16, which are equally distributed on the four sides of the inversion domain, and all sources are added in the form of P-wave incident excitation. To show the denoising performance of the FD-MRCSI method, the numerical inversion results are presented when the measurement datasets are added 5% random noise. Figures 1(b) and 1(c) are the reconstruction profile of the experiment and discrepancy between the true model and recovering model after 1000 iterations. Figure 1(d) shows the slices along the depth at horizontal distances of 3 km for the true model and inversion model. The relative error is 12.26%.In this type of numerical tests of a homogeneous background medium, we will test that the FD-MRCSI method possesses the ability of reconstructing multiple scatterers with different geometrical information and different physical property parameters. For this purpose, we employ a model consisting of two rectangular scatterers as shown in Figure2(a). The configuration consists of two rectangular objects, which have different dimensions. One is 1.2 km × 1.2 km, and another is 0.8 km × 1 km. We also assume that the inversion domain is a square domain with the dimension of 6 km × 6 km and has a second Lame coefficient μ of 5.2 GPa, a Poisson’s ratio of 0.25, and a mass density of 1.5 g/cm3. The large obstacle has a mass density of 2.7 g/cm3, and the small one has a mass density of 2.1 g/cm3. Hence, the contrasts for the two abnormities are 0.8 and 0.4, respectively.Applying the FD-CSI method to two abnormities inversion. (a) True model. (b) Reconstruction result after 1000 inversion iterations. (c) Discrepancy between true model and reconstruction result. (d) Variation of the data term error with respect to iterations.
(a)
True model
(b)
Reconstruction result
(c)
Discrepancy
(d)
Variation of the data equation errorIn the inversion process, we discretize the total computation domain into151
×
151 uniform grids, including 15 PML grids, and 24 sources and 24 receivers are used for this example. These sources and receivers are distributed along the inversion domain, six sources or receivers for each side. After the measurement data is generated using the same way as the previous example, we select them at frequency 8 Hz to drive the inversion algorithm, and 5% noise is added. After 1000 iterations, the reconstruction profile is given in Figure 2(b). Figure 2(c) denotes the discrepancy between the true model and recovering model after 1000 iterations, and Figure 2(d) shows the variation of the data equation term error with respect to iterations. The relative error is 16.05%.At first glance, we should solve the forward modeling for many times for these 1000 inversion iterations. Since the contrast sourcew and contrast function χ are incorporated into this algorithm, which makes the forward stiffness matrix only involve the background medium for each given frequency, the LU decomposition is required only once throughout the whole inversion procedure; therefore, the forward modeling computation cost is very cheap. However, the price that we have to pay is the storing of the LU decomposition arrays of the stiffness matrix of the background medium.From these numerical results, we observe that the FD-MRCSI method can effectively reconstruct the scatterers including their locations and shapes using the measurement scattered data, and the method is of strong denoising ability. Additionally, Figure2(d) demonstrates that the FD-MRCSI method is convergent.


### 4.2.  Reconstruction in an Inhomogeneous Background Medium

In this section, we examine the performance of the FD-CSI using an inhomogeneous medium as a background medium. The configuration of the background consists of a compressional wave speed of 2500 m/s, a shear wave speed of 1443 m/s, and a four-layer distribution of a mass density of 1.8, 1.6, 1.5, and 1.8 g/cm3 from top to bottom. Figure 3(a) shows the mass density distribution of the background medium. The inversion model consisting of a circle and an ellipse is presented as shown in Figure 3(b). The circular abnormity with a contrast of 0.8 has a radius of 0.550 km, which is centered at (3.4, 3.6) km. For the elliptical scatterer, it has a semimajor axis of 0.7 km. and a semiminor axis of 0.36 km, and its centre is (2.3, 2) km. We also assume that the inversion domain is a square with dimension of L
×
L. In our numerical experiments, L is assigned to 6 km.Applying the FD-CSI method to two abnormities inversion in an inhomogeneous background medium. (a) Background medium. (b) True model. (c) Reconstruction result after 1000 inversion iterations. (d) Variation of the data term with respect to iterations.
(a)
Background medium
(b)
True model
(c)
Reconstruction result
(d)
Variation of the data equation errorIn the inversion, the contrast functionχ, contrast sources w
, and state variables u, u
inc, and u
sct in the inversion domain are uniformly discretized on a regular two-dimensional grid at intervals of Δ
x and Δ
z. Both Δ
x and Δ
z are equal to 0.05 km; in other words, the grid size is 121 × 121. The measurement domain surface S where 24 receivers and 24 sources are distributed is a square around the inversion domain. After the synthetic datasets are generated by using FFT on the finite-difference time-domain (FD-DT) data at frequency 10 Hz, 5% noise is added by using (35). After running 1000 iterations, the reconstruction results are given in Figure 3(c). The relative errors are 18.83%.As shown in Figure3, the FD-MRCSI method can successfully reconstruct the objects including their locations and shapes as those present in a homogeneous background medium. The feature of the inversion method results from the flexibility of the finite-difference operator of the governing differential equations.


### 4.3.  Reconstruction Improvement Using Multiple Frequencies

In these examples, we use the previous lower frequency inversion reconstruction results as an initial guess for another frequency inversion, which is the so-called sequential multifrequency inversion approach, that one usually employs this technique in the Gauss-Newton or nonlinear conjugate gradient method. This technique applied to improve the inversion profiles is guaranteed by the following two observations: (1) each frequency is independent of FD-MRCSI algorithm; (2) the inversion algorithm is locally convergent, and hence a good initial guess will result in a better inversion results. Additionally, multiple frequencies inversion technique can migrate the nonlinearity of the elastic inverse problem.To drive the inversion algorithm, we should use another initialization method, which was introduced by Abubakar et al. (see [8]). Firstly, use the previous inversion result as χ
0; then, solve the forward problem with ρ
=
(
1
+
χ
0
)
ρ
b to get the total field vectors u
j
,
0; finally, the initial contrast sources w
j
,
0 are obtained with the equation w
j
,
0
=
χ
0
u
j
,
0.We employ the strategy to rebuild the two test examples given in Section4.1, and all relative parameters are the same as those present in Section 4.1. For the one scatterer model case, we successively select three frequencies (4, 6, and 8 Hz). For the two scatterers model case, the sequential frequency applied during the multiple frequencies inversion is 5 Hz, 7 Hz, and 10 Hz. Figures 4(a) and 4(b) show the corresponding reconstruction results. The relative errors are 8.13% and 12.60%, respectively, which are smaller than those in Section 4.1. Table 2 indicates the reconstruction relative errors of the two models using FD-MRCSI and FD-MRCSI with multiple frequencies technique.Table 2
Relative errors using FD-CSI and FD-CSI with multiple frequencies technique.
Model
FD-CSI
FD-CSI with multifrequency
One scatterer model
12.26%
8.13%
Two scatterers model
16.05%
12.60%Applying sequential multifrequency technique to improve reconstruction results. (a) Reconstruction result of one square scatterer model. (b) Reconstruction result of two rectangle scatterers model.
(a)
Multifrequency inversion result
(b)
Multifrequency inversion resultFrom Figure4 and Table 2, we can observe that the use of the sequential multifrequency approach can improve reconstructed model to some extent. This can be explained by the fact that a good initial guess can result in a better computation result for a locally convergent method.


## 4.1. Reconstruction in a Homogeneous Background Medium

The first example is concerned with recovering a square scatterer from noise observed datad
j
noise in a homogeneous background medium, and the exact model is given by Figure 1(a). The inversion domain D is set as a square with side length 6 km, which does not include PML thickness. For this square scatterer, the second Lame parameter μ is equal to 2.0 GPa, and a Poisson’s ratio is 0.25. The square with the side length 1.5 km is located in the inversion domain, and its centre is at (3.3) km; the mass density of the scatterer is 1.8 g/cm3, 1.5 g/cm3 outside the scatterer. Hence, the contrast χ for the scatterer is 0.2. During the whole inversion procedure, the angular frequency ω is 8 Hz.Applying the FD-MRCSI method to one abnormity inversion. (a) True model. (b) Reconstruction result after 1000 inversion iterations. (c) Discrepancy between true model and reconstruction result. (d) Slices at horizontal distances of 3 km for the true model and reconstruction result.
(a)
True model
(b)
Reconstruction result
(c)
Discrepancy
(d)
Model sliceThe inversion domainD is divided into N
×
N uniform grids using the second-order scheme, adding the n
1 grids of PML; so, the size of the total computation domain is (
N
+
2
n
1
)
×
(
N
+
2
n
1
). For this experiment, N
=
61 and n
1
=
15; therefore, the grid dimension of the total computation domain is 91
×
91. Both the sources and receivers are assigned to 16, which are equally distributed on the four sides of the inversion domain, and all sources are added in the form of P-wave incident excitation. To show the denoising performance of the FD-MRCSI method, the numerical inversion results are presented when the measurement datasets are added 5% random noise. Figures 1(b) and 1(c) are the reconstruction profile of the experiment and discrepancy between the true model and recovering model after 1000 iterations. Figure 1(d) shows the slices along the depth at horizontal distances of 3 km for the true model and inversion model. The relative error is 12.26%.In this type of numerical tests of a homogeneous background medium, we will test that the FD-MRCSI method possesses the ability of reconstructing multiple scatterers with different geometrical information and different physical property parameters. For this purpose, we employ a model consisting of two rectangular scatterers as shown in Figure2(a). The configuration consists of two rectangular objects, which have different dimensions. One is 1.2 km × 1.2 km, and another is 0.8 km × 1 km. We also assume that the inversion domain is a square domain with the dimension of 6 km × 6 km and has a second Lame coefficient μ of 5.2 GPa, a Poisson’s ratio of 0.25, and a mass density of 1.5 g/cm3. The large obstacle has a mass density of 2.7 g/cm3, and the small one has a mass density of 2.1 g/cm3. Hence, the contrasts for the two abnormities are 0.8 and 0.4, respectively.Applying the FD-CSI method to two abnormities inversion. (a) True model. (b) Reconstruction result after 1000 inversion iterations. (c) Discrepancy between true model and reconstruction result. (d) Variation of the data term error with respect to iterations.
(a)
True model
(b)
Reconstruction result
(c)
Discrepancy
(d)
Variation of the data equation errorIn the inversion process, we discretize the total computation domain into151
×
151 uniform grids, including 15 PML grids, and 24 sources and 24 receivers are used for this example. These sources and receivers are distributed along the inversion domain, six sources or receivers for each side. After the measurement data is generated using the same way as the previous example, we select them at frequency 8 Hz to drive the inversion algorithm, and 5% noise is added. After 1000 iterations, the reconstruction profile is given in Figure 2(b). Figure 2(c) denotes the discrepancy between the true model and recovering model after 1000 iterations, and Figure 2(d) shows the variation of the data equation term error with respect to iterations. The relative error is 16.05%.At first glance, we should solve the forward modeling for many times for these 1000 inversion iterations. Since the contrast sourcew and contrast function χ are incorporated into this algorithm, which makes the forward stiffness matrix only involve the background medium for each given frequency, the LU decomposition is required only once throughout the whole inversion procedure; therefore, the forward modeling computation cost is very cheap. However, the price that we have to pay is the storing of the LU decomposition arrays of the stiffness matrix of the background medium.From these numerical results, we observe that the FD-MRCSI method can effectively reconstruct the scatterers including their locations and shapes using the measurement scattered data, and the method is of strong denoising ability. Additionally, Figure2(d) demonstrates that the FD-MRCSI method is convergent.


## 4.2.  Reconstruction in an Inhomogeneous Background Medium

In this section, we examine the performance of the FD-CSI using an inhomogeneous medium as a background medium. The configuration of the background consists of a compressional wave speed of 2500 m/s, a shear wave speed of 1443 m/s, and a four-layer distribution of a mass density of 1.8, 1.6, 1.5, and 1.8 g/cm3 from top to bottom. Figure 3(a) shows the mass density distribution of the background medium. The inversion model consisting of a circle and an ellipse is presented as shown in Figure 3(b). The circular abnormity with a contrast of 0.8 has a radius of 0.550 km, which is centered at (3.4, 3.6) km. For the elliptical scatterer, it has a semimajor axis of 0.7 km. and a semiminor axis of 0.36 km, and its centre is (2.3, 2) km. We also assume that the inversion domain is a square with dimension of L
×
L. In our numerical experiments, L is assigned to 6 km.Applying the FD-CSI method to two abnormities inversion in an inhomogeneous background medium. (a) Background medium. (b) True model. (c) Reconstruction result after 1000 inversion iterations. (d) Variation of the data term with respect to iterations.
(a)
Background medium
(b)
True model
(c)
Reconstruction result
(d)
Variation of the data equation errorIn the inversion, the contrast functionχ, contrast sources w
, and state variables u, u
inc, and u
sct in the inversion domain are uniformly discretized on a regular two-dimensional grid at intervals of Δ
x and Δ
z. Both Δ
x and Δ
z are equal to 0.05 km; in other words, the grid size is 121 × 121. The measurement domain surface S where 24 receivers and 24 sources are distributed is a square around the inversion domain. After the synthetic datasets are generated by using FFT on the finite-difference time-domain (FD-DT) data at frequency 10 Hz, 5% noise is added by using (35). After running 1000 iterations, the reconstruction results are given in Figure 3(c). The relative errors are 18.83%.As shown in Figure3, the FD-MRCSI method can successfully reconstruct the objects including their locations and shapes as those present in a homogeneous background medium. The feature of the inversion method results from the flexibility of the finite-difference operator of the governing differential equations.


## 4.3.  Reconstruction Improvement Using Multiple Frequencies

In these examples, we use the previous lower frequency inversion reconstruction results as an initial guess for another frequency inversion, which is the so-called sequential multifrequency inversion approach, that one usually employs this technique in the Gauss-Newton or nonlinear conjugate gradient method. This technique applied to improve the inversion profiles is guaranteed by the following two observations: (1) each frequency is independent of FD-MRCSI algorithm; (2) the inversion algorithm is locally convergent, and hence a good initial guess will result in a better inversion results. Additionally, multiple frequencies inversion technique can migrate the nonlinearity of the elastic inverse problem.To drive the inversion algorithm, we should use another initialization method, which was introduced by Abubakar et al. (see [8]). Firstly, use the previous inversion result as χ
0; then, solve the forward problem with ρ
=
(
1
+
χ
0
)
ρ
b to get the total field vectors u
j
,
0; finally, the initial contrast sources w
j
,
0 are obtained with the equation w
j
,
0
=
χ
0
u
j
,
0.We employ the strategy to rebuild the two test examples given in Section4.1, and all relative parameters are the same as those present in Section 4.1. For the one scatterer model case, we successively select three frequencies (4, 6, and 8 Hz). For the two scatterers model case, the sequential frequency applied during the multiple frequencies inversion is 5 Hz, 7 Hz, and 10 Hz. Figures 4(a) and 4(b) show the corresponding reconstruction results. The relative errors are 8.13% and 12.60%, respectively, which are smaller than those in Section 4.1. Table 2 indicates the reconstruction relative errors of the two models using FD-MRCSI and FD-MRCSI with multiple frequencies technique.Table 2
Relative errors using FD-CSI and FD-CSI with multiple frequencies technique.
Model
FD-CSI
FD-CSI with multifrequency
One scatterer model
12.26%
8.13%
Two scatterers model
16.05%
12.60%Applying sequential multifrequency technique to improve reconstruction results. (a) Reconstruction result of one square scatterer model. (b) Reconstruction result of two rectangle scatterers model.
(a)
Multifrequency inversion result
(b)
Multifrequency inversion resultFrom Figure4 and Table 2, we can observe that the use of the sequential multifrequency approach can improve reconstructed model to some extent. This can be explained by the fact that a good initial guess can result in a better computation result for a locally convergent method.


## 5. Conclusions

In this paper, the finite-difference contrast source inversion method incorporating a multiplicative regularization factor is successfully employed in seismic full-waveform inversion to reconstruct the mass density. The method is very efficient because it does not require solving a full forward problem in each inversion iteration, and the highly efficient calculation of the forward modeling makes it have a potential for dealing with large scale computational inverse problems.From the results of numerical experiments, the algorithm can well recognize the shapes of the scatterers and precisely identify their locations, no matter the background mediums are homogeneous or inhomogeneous medium. The algorithm is very robust when the measurement data is polluted by random noise. Furthermore, the FD-MRCSI using sequential multifrequency technique can successfully improve the reconstruction results, which is very hopeful for getting higher resolution in exploration seismology. Estimating other subsurface parameters using this method, such as Lame coefficients and velocity, is our future work.



---

*Source: 532159-2014-08-26.xml*