1996 C. Jordi et al.
inversions on a variety of scales and targets, such as in hydrological
near-surface studies (e.g. Linde et al. 2006;Doetschet al. 2010b),
targets in mineral exploration (e.g. Leli
`
evre et al. 2012) or hydro-
carbon exploration (e.g. Colombo et al. 2012). The majority of the
2-D applications (e.g. Gallardo 2007;Huet al. 2009; Colombo et al.
2012) as well as 3-D examples (Linde et al. 2006; Fregoso & Gal-
lardo 2009;Doetschet al. 2010b; Moorkamp et al. 2011;Shiet al.
2017) use regular grid discretizations. Only Leli
`
evre et al. (2012)
work on irregular tetrahedral (triangular in 2-D) grids. They show
an application of structural joint inversion using cross-gradients on
irregular meshes in a 2-D synthetic example.
While irregular meshes are superior in representing complicated
subsurface structures, topography and complex acquisition geome-
tries (R
¨
ucker et al. 2006), the design of gradient and regulariza-
tion operators is not trivial for irregular discretizations (Leli
`
evre &
Farquharson 2013). A suitable option for regularization on irreg-
ular meshes are geostatistical regularization operators (Jordi et al.
2018). Geostatistical regularization is based on a correlation model
and allows to incorporate prior knowledge such as length scales over
which geological properties are expected to be correlated (correla-
tion lengths). In particular, for geological settings with a preferential
direction (e.g. layering), geostatistical regularization can improve
the inversion results compared to traditional cell-to-cell smooth-
ness constraints (e.g. Hermans et al. 2012; Jordi et al. 2018).
Concerning the model gradients that are needed to calculate
the cross-gradients for structural coupling, Leli
`
evre & Farquharson
(2013) propose to calculate the components of the model gradient
at the centre of a particular cell by fitting a linear trend through the
model values of that cell and its direct neighbours. For each cell,
a small inverse problem is solved, resulting in the sought model
gradients. We combine the ideas of Leli
`
evre & Farquharson (2013)
for gradient calculation and those of Jordi et al. (2018) for design-
ing the operator size based on a geostatistical model, in order to
develop a flexible and robust way of calculating cross-gradients on
unstructured meshes for structural coupling. Instead of only using
the nearest neighbours of a cell, we use a correlation-based opera-
tor design that includes a larger neighbourhood around a particular
cell. The main advantage of the correlation-based operators is their
reduced dependency on a particular mesh discretization. Due to the
underlying correlation model, they act on a physical length-scale
rather than on the length-scale that is prescribed by the cell sizes
of a mesh. In particular, for irregular meshes, where cell sizes and
distances between neighbouring cells can vary, decoupling of the
operators from the mesh is important.
Using these operators, we present a structural joint inversion al-
gorithm that is suitable for 2-D and 3-D joint inversions on irregular
meshes. We tested our joint inversion on a 3-D synthetic cross-well
example, typical for hydrological near-surface studies. In a second
example, joint inversion was applied to 2-D seismic refraction and
ERT data that were recorded in a complex karstified area.
2 METHODS
2.1 Joint inversion scheme
The objective function of the joint inverse problem has the following
form:
=
d
+ λ
m
+ λ
cg
cg
(1)
where λ and λ
cg
are the geostatistical regularization and cross-
gradient weights, respectively.
d
is the combined data misfit term
for all data types
d
=
K
k=1
C
−0.5
d,k
(d
k
− F
k
(m
k
)
2
2
, (2)
where K is the number of data types in the joint inverse problem. We
consider here only the case K=2, but extension to more data types
is possible (Gallardo 2007;Doetschet al. 2010b). C
d,k
is the data
covariance matrix (assumed to be diagonal), d
k
are the observed
data, F
k
(m
k
) is the calculated forward response of the model m
k
.
Theregularizationtermisdefinedas
m
=
K
k=1
C
−0.5
m,k
m
k
2
2
. (3)
C
−0.5
m,k
is the regularization operator for model m
k
and m
k
is the
difference between the model m
k
and a reference model m
ref
k
.The
cross-gradients term for two data types (K = 2) and the correspond-
ing models m
1
and m
2
is
cg
=t(m
1
,m
2
)
2
2
, (4)
where t(m
1
,m
2
) is the cross-product calculated between m
1
and m
2
, that is the model updates with respect to the reference
model (see below).
For two different data types d
1
(e.g. apparent resistivity) and d
2
(e.g. traveltime), the following system of equations is solved in a
least-squares sense (adapted from Linde et al. 2006):
⎡
⎢
⎢
⎢
⎢
⎣
C
−0.5
d
J
p
C
−0.5
m
λ
cg
B
p
x
λ
cg
B
p
y
λ
cg
B
p
z
⎤
⎥
⎥
⎥
⎥
⎦
m
p+1
=
⎡
⎢
⎢
⎢
⎢
⎣
C
−0.5
d
ˆ
d
p
C
−0.5
m
m
p
λ
cg
(B
p
x
m
p
− t
p
x
)
λ
cg
(B
p
y
m
p
− t
p
y
)
λ
cg
(B
p
z
m
p
− t
p
z
)
⎤
⎥
⎥
⎥
⎥
⎦
, (5)
where superscript p indicates the iteration number,
m
p+1
=
m
p+1
1
m
p+1
2
, J
p
=
J
p
1
J
p
2
,
ˆ
d
p
=
d
1
− F
1
(m
p
1
) + J
p
1
m
p
1
d
2
− F
2
(m
p
2
) + J
p
2
m
p
2
,
C
−0.5
d
=
w
1
C
−0.5
d1
0
0 w
2
C
−0.5
d2
,
C
−0.5
m
=
λ
p
w
1
C
−0.5
m1
0
0 λ
p
w
2
C
−0.5
m2
.
m
p+1
1
(e.g. resistivity) and m
p+1
2
(e.g. slowness) are the sought
model updates that are used to iteratively update the model vectors
with respect to the reference models m
ref
1,2
(Doetsch et al. 2010a):
m
p+1
1
m
p+1
2
=
m
ref
1
m
ref
2
+
m
p+1
1
m
p+1
2
. (6)
J
p
1
and J
p
2
are the Jacobian matrices containing the sensitivities of
a change in data d
1,2
related to a change in model m
1,2
, w
1
and w
2
are weighting factors (see below), λ
p
is the overall regularization
weight that is allowed to change for each iteration and λ
cg
is the
cross-gradient weight. B
p
x,y,z
are the Jacobians of the cross-gradient
functions t
p
x,y,z
with respect to changes in the respective models and
respective directions. These sensitivities, for example, B
p
x
,areused
to estimate the cross-gradient function t
p+1
x
at the present iteration
by a first order Taylor expansion around the cross-gradient function
t
p
x
of the previous iteration (Linde et al. 2008):
t
p+1
x
∼
=
t
p
x
+ B
p
x
m
p
. (7)
Downloaded from https://academic.oup.com/gji/article/220/3/1995/5658692 by China University of Geosciences user on 12 December 2022