viabilidade da aplicação do método quase-Newton

Prévia do material em texto

UNIVERSIDADE FEDERAL DO PARÁ
INSTITUTO DE GEOCIÊNCIAS
PROGRAMA DE PÓS-GRADUAÇÃO EM GEOFÍSICA
DISSERTAÇÃO DE MESTRADO
Analysis of quasi-Newton inversion
of magnetoteluric data
ALEX FABRICIO DE ANDRADE SILVA
Belém – Pará
2020
ALEX FABRICIO DE ANDRADE SILVA
Analysis of quasi-Newton inversion
of magnetoteluric data
Dissertação apresentada ao Programa de Pós-Graduação
em Geofísica do Instituto de Geociências da Universi-
dade Federal do Pará para obtenção do título de Mestre
em Geofísica.
Área de concentração:
Modelagem e inversão de dados geofísicos.
Linha de pesquisa:
Geofísica aplicada à exploração de hidrocarbonetos
Orientador: Cícero Roberto Teixeira Régis
Belém – Pará
2020
Dados Internacionais de Catalogação na Publicação (CIP) de acordo com ISBD
Sistema de Bibliotecas da Universidade Federal do Pará
Gerada automaticamente pelo módulo Ficat, mediante os dados fornecidos pelo(a) autor(a)
S586a Silva, Alex Fabricio de Andrade.
 Analysis of quasi-Newton inversion of magnetoteluric data /
Alex Fabricio de Andrade Silva. — 2020.
 57 f. : il. color.
 Orientador(a): Prof. Dr. Cícero Roberto Teixeira Régis
 Dissertação (Mestrado) - Universidade Federal do Pará,
Instituto de Geociências, Programa de Pós-Graduação em
Geofísica, Belém, 2020.
 1. Inversão. 2. Método quase-Newton. 3. Método
Magnetotelúrico. I. Título.
CDD 550
Powered by TCPDF (www.tcpdf.org)
ALEX FABRICIO DE ANDRADE SILVA
Analysis of quasi-Newton inversion
of magnetoteluric data
Dissertação apresentada ao Programa de Pós-Graduação
em Geofísica do Instituto de Geociências da Universi-
dade Federal do Pará para obtenção do título de Mestre
em Geofísica.
Data de aprovação: 23 de outubro de 2020
Banca Examinadora:
Cícero Roberto Teixeira Régis (Orientador)
Universidade Federal do Pará
Prof. Dr. Paulo de Tarso Luiz Menezes
Universidade do Estado do Rio de Janeiro
Prof. Dr. Marcos Welby Correa Silva
Universidade Federal do Pará
À minha família,
que muito me incentiva e apoia.
AGRADECIMENTOS
Ao professor Cícero Roberto Teixeira Régis pelos ensinamentos e paciência nas orien-
tações, bem como pela sugestão do tema.
Ao professor Paulo de Tarso Luiz Menezes pelas correções e sugestões para o apri-
moramento do presente trabalho.
Ao professor Marcos Welby Correa Silva pela participação com recomendações e
ajustes que colaboraram para o aperfeiçoamento deste trabalho.
Aos meus colegas do PROEM - PPGf/UFPa, em particular ao Carlos Mateus pelo
esclarecimento de questões sobre métodos geofísicos e pelo auxílio na implementação do
PARADISO.
Aos meus amigos Jéssica Larissa, Karolina Almeida, René Michel e Tarciana Valéria,
que contribuíram com boas conversas e tornaram este período mais agradável.
Ao convênio Petrobrás/UFPa/Fadesp pelo apoio financeiro com bolsa de mestrado do
projeto de pesquisa 4600572512.
RESUMO
O presente trabalho analisa a viabilidade da aplicação do método quase-Newton de
Broyden-Fletcher-Goldfarb-Shanno (BFGS) de memória limitada na inversão de dados
magnetotelúricos. O BFGS de memória limitada (ou simplesmente LM-BFGS) consiste
em um algoritmo de otimização da classe de métodos quase-Newton que usa uma quanti-
dade menor de memória por não armazenar matrizes Jacobianas e Hessianas. A aplicação
consiste em inverter dados sintéticos magnetotelúricos de modelos bidimensionais, com
variações topográficas no relevo. Uma análise comparativa mostra o ganho de memória
alcançado pelo método LM-BFGS em relação ao método de Gauss-Newton.
Palavras-chave: Inversão. Método quase-Newton. Método Magnetotelúrico.
ABSTRACT
This work explores the application of the limited memory Broyden-Fletcher-Goldfarb-
Shanno (BFGS) quasi-Newton method in the inversion of magnetotelluric data. The
limited memory BFGS (or LM-BFGS) consists of an optimization algorithm of the quasi-
Newton methods category that uses a smaller amount of computer memory because it
does not need storing Jacobian and Hessian matrices. The application is to invert magne-
totelluric synthetic data from two-dimensional models, including variations in topographic
relief. A comparative analysis shows the gain in memory requirements achieved by the
LM-BFGS method quantitatively to the traditional Gauss-Newton method.
Keywords: Inversion. quasi-Newton method. Magnetotelluric method.
LIST OF FIGURES
2.1 Scheme of a geoelectric model typically structured in geophysical electro-
magnetic investigations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Schemes illustrating the directions of associations between the parameters
in the two-dimensional case . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Topographic surface model, a ramp. . . . . . . . . . . . . . . . . . . . . . . 20
3.1 Example of iteration, time and memory spent by both methods. . . . . . . 21
3.2 Values stored and iterations by both methods for 2D inversion. . . . . . . . 22
3.3 Model 1 - Ramp and rectangular body. . . . . . . . . . . . . . . . . . . . . 23
3.4 The absolute values of successive φ differences in the inversion of model 1. 24
3.5 Result of the model 1 inversion with Gauss-Newton for a grid with 2720
cells. The true target body position is shown by the white line. . . . . . . . 25
3.6 Results of the model 1 inversion with the LM-BFGS quasi-Newton method
for a 2720 cells grid. The true target body position is shown by the white
line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.7 Model 1 inversions in contour lines . . . . . . . . . . . . . . . . . . . . . . 26
3.8 Observed and estimated data for the first station in model 1. . . . . . . . . 27
3.9 Observed and estimated data for the last station in model 1. . . . . . . . . 27
3.10 Apparent resistivity pseudosection (Model 1 - Observation). . . . . . . . . 28
3.11 Apparent resistivity pseudosection (Model 1 - Gauss-Newton). . . . . . . . 28
3.12 Apparent resistivity pseudosection (Model 1 - LM-BFGS). . . . . . . . . . 28
3.13 Phase angle pseudosection (Model 1 - Observation). . . . . . . . . . . . . . 29
3.14 Phase angle pseudosection (Model 1 - Gauss-Newton). . . . . . . . . . . . . 29
3.15 Phase angle pseudosection (Model 1 - LM-BFGS). . . . . . . . . . . . . . . 29
3.16 Model 2 - Medium with three target bodies and complex topographic relief. 30
3.17 The absolute values of successive φ differences in the inversion of model 2. 31
3.18 Model 2 inversion with Gauss-Newton for a grid with 2112 cells . . . . . . 32
3.19 Model 2 inversion with LM-BFGS for a 2112 cells grid . . . . . . . . . . . 32
3.20 Model 2 inversions in contour lines . . . . . . . . . . . . . . . . . . . . . . 33
3.21 Observation and calculated data for the first station (model 2) . . . . . . . 34
3.22 Observation and calculated data for the last station (model 2) . . . . . . . 34
3.23 Apparent resistivity pseudosection (Model 2 - Observation). . . . . . . . . 35
3.24 Apparent resistivity pseudosection (Model 2 - Gauss-Newton). . . . . . . . 35
3.25 Apparent resistivity pseudosection (Model 2 - LM-BFGS). . . . . . . . . . 35
3.26 Phase angle pseudosection (Model 2 - Observation). . . . . . . . . . . . . . 36
3.27 Phase angle pseudosection (Model 2 - Gauss-Newton). . . . . . . . . . . . . 36
3.28 Phase angle pseudosection (Model 2 - LM-BFGS). . . . . . . . . . . . . . . 36
3.29 Model 3 - Simplified model based on Reyes-Wagner et al. (2017) of a re-
gional line in South Volcanic Zone (SVZ) of the Chilean Andes. . . . . . . 37
3.30 The absolute values of successive φ differences in the inversion of model 3. 38
3.31 Result of the model 3 inversion with Gauss-Newton for a grid with 3500
cells. The true target bodies positions are shown by the white lines. . . . . 39
3.32 Results of the model 3 inversion with the LM-BFGS quasi-Newton method
for a 3500 cells grid. The true target bodies positions are shown by the
white lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.33 Model 3 inversions in contour lines . . . . . . .. . . . . . . . . . . . . . . 40
3.34 Observed and estimated data for the first station in model 3. . . . . . . . . 41
3.35 Observed and estimated data for the last station in model 3. . . . . . . . . 41
3.36 Apparent resistivity pseudosection (Model 3 - Observation). . . . . . . . . 42
3.37 Apparent resistivity pseudosection (Model 3 - Gauss-Newton). . . . . . . . 42
3.38 Apparent resistivity pseudosection (Model 3 - LM-BFGS). . . . . . . . . . 42
3.39 Phase angle pseudosection (Model 3 - Observation). . . . . . . . . . . . . . 43
3.40 Phase angle pseudosection (Model 3 - Gauss-Newton). . . . . . . . . . . . . 43
3.41 Phase angle pseudosection (Model 3 - LM-BFGS). . . . . . . . . . . . . . . 43
CONTENTS
1 INTRODUCTION 1
2 METHODOLOGY 3
2.1 THE MAGNETOTELLURIC METHOD . . . . . . . . . . . . . . . . . . 5
2.2 INVERSE PROBLEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.1 Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.2 Gauss-Newton method . . . . . . . . . . . . . . . . . . . . . . 10
2.2.3 Quasi-Newton method . . . . . . . . . . . . . . . . . . . . . . 11
2.2.4 Limited-memory BFGS (LM-BFGS) . . . . . . . . . . . . . 13
2.3 CALCULATION OF SENSITIVITIES . . . . . . . . . . . . . . . . . . . . 18
2.4 SPACIAL DERIVATIVES . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 RESULTS - APPLICATION TO THE INVERSION OF MAGNE-
TOTELLURIC DATA 21
3.1 MODEL 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 MODEL 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.3 MODEL 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4 CONCLUSIONS 44
REFERENCES 45
1 INTRODUCTION
Quasi-Newton (QN) methods to invert geophysical data work by creating an iterative
approximation to the Hessian matrix that appears in the Gauss-Newton method, or to its
inverse. In a given iteration of the inversion, these methods estimate the Hessian by using
the gradient information from some or all previous iterations. The matrix is updated and
adjusted on each iteration and it can be produced in several different ways.
Broyden (1970a,b) introduced an important quasi-Newton method and presented gen-
eral considerations about the class of methods known as the Broyden class. Fletcher (1970)
reached a formula from an algebraic derivation of the first known quasi-Newton method,
the Davidon-Fletcher-Powell (DFP) formula. Goldfarb (1970) developed the same for-
mula from the Greenstadt (1970) paper. Shanno (1970) developed this strategy to solve
a conditioning problem of the DFP method. These authors worked independently and
developed what is now known as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-
Newton method.
The method approximates the Hessian (or its inverse) from a formula that uses all
previous iterations. The difference with the limited memory BFGS (L-BFGS or LM-
BFGS) is in the storage of this information for the approximation. The limited memory
methodology uses only the information from the most recent iterations.
The formulation of the limited memory strategy for quasi-Newton methods, originated
in the works of Perry (1977) and Shanno (1978a), were later developed and analyzed
by Buckley (1978), Nazareth (1979), Nocedal (1980), Shanno (1978b), Gill and Murray
(1979), Buckley and Lenir (1983) and it is based on the storage of pairs of vectors,
gradients and parameters, of the most recent iterations during the inversion process to
be used in approximations of the products of the inverse of the Hessian matrix with the
gradient vector. The LM-BFGS performs the estimates for the product of the Hessian
inverse matrix with the gradient vector using an approximation algorithm.
The LM-BFGS has the advantage of not saving matrices (Jacobian and Hessian) in
memory during the inversion iterations. The number of the recent vectors stored depends
on the amount of memory used and influences the approximations made, reducing or
increasing the total number of iterations of the inversion.
The main advantage of implementing the LM-BFGS method is the lower demand for
computational memory in the inverse problem than other often used inversion methods,
such as the Gauss-Newton method. In addition, the LM-BFGS method is particularly ap-
plicable for optimization problems with a large number of variables, such as the inversion
of 2D and 3D geophysical data.
1
2
Over the past few decades, many studies have shown that the use of quasi-Newton
methods, such as the LM-BFGS algorithm, requires more iterations when compared to the
Newton method or the Gauss-Newton method, but with a total memory demand smaller,
as presented by Nocedal and Wright (1999).
The LM-BFGS method has been used for many large-scale problems. In inverse prob-
lems of electromagnetic geophysical data, quasi-Newton methods have become popular
tools based on the work of Newman and Boggs (2004) and Haber (2005). Avdeeva and
Avdeev (2006) present the quasi-Newton LM-BFGS applied to the 1D magnetotelluric
inverse problem, as a first step for the larger 3D inverse problems, based on the opti-
mization work of Ni and Yuan (1997) and the work to implement the Wolfe conditions
according to Byrd et al. (1995).
This work evaluates the feasibility of applying the LM-BFGS quasi-Newton method
comparing the inversion of magnetotelluric data with the Gauss-Newton method, in quali-
tative (data adjustment and inversion) and quantitative (spent computer memory) terms.
The analysis’s motivation is the search for an inversion alternative that requires a smaller
amount of computational memory, due to the limitations of maintenance and obtaining
robust machines for processing problems that require a large volume of memory. In an-
alyzing the possibility of applying the method, the amount of time required to achieve a
satisfactory inversion result was also observed.
2 METHODOLOGY
Geophysical investigations are based on the indirect measurement of physical variables,
without the need for direct contact with the site to be studied in a subsurface. The
geophysical exploration of the Earth’s interior consists of making measurements on or
near the surface, where the internal distribution of physical properties influences these
measurements. These measurements make up sets of observations, which is sometimes
the only one information we have, contaminated with noise. The observations can be
established in terms of fields components, such as magnetic or electrical.
The subsurface properties are then analyzed by measuring, processing and interpreting
the data. In the analysis and handling of data, we make the modeling and inversion
processes. The modeling or forward problem, in short, consists of establishing simulations
through mathematical models. In the inverse problem, the values and distribution of a
specific physical property of the subsurface are estimated that best fit the observations to
the modeling.
The observation data and the known information of the region of interest offer us
a set of hypotheses to set up mathematically a function that correlates the hypotheses
and the observations. The function is characterized by the parameters and based on
these it is possible to determine the data set by means of which it is possible to check,
compared with the observation data, if the hypotheses are sufficient in the explanation of
the observations.
The forward problem is based on the function between the parameters, which configure
the hypotheses, and the observation data. The forward problem in the present work
are based on the development of electromagnetic theory in the geophysical perspective,
governed by Maxwell’s equations.
To stipulate the geoelectric features on the subsurface, the electromagnetic geophysical
methods are based on the physical considerations of the fields generated from the trans-
mitters (Tx ), the configuration of the receivers (Rx ), the formations and the geological
environment .
In the frequency domain, Maxwell’s equationsfor conductive media are
∇ · �E = 0, (2.1)
∇×H − (σ + iω�)E = Jext, (2.2)
∇× E + iωµH = 0, (2.3)
∇ · µH = 0, (2.4)
where the electric current of the transmitter is represented by Jext.
3
4
Figure 2.1: Scheme of a geoelectric model typically structured in geophysical electromagnetic
investigations.
A geoelectric model typically idealized in electromagnetic geophysical studies is shown
in the figure 2.1. In the scheme, we can differentiate two particular regions: one referring
to the embedding medium, that is, the stratified medium without the heterogeneities and
the other representing the heterogeneities.
Taking z = iωµ as the impeditivity of the medium and η = σ + iω� its admittivity, in
an electrically homogeneous medium Maxwell’s equations condense into
∇×H − ηE = Jext, (2.5)
∇× E + zH = 0. (2.6)
In order to have these equations related to the fields induced by the heterogeneities,
the magnetic and electric fields are separated in primary and secondary terms:
E = Ep + Es, (2.7)
H = Hp +Hs. (2.8)
Thus, in terms of primary and secondary fields, Maxwell’s equations are written as
∇×Hs − ηEs = ∆ηEp, (2.9)
∇× Es + zHs = −∆zHp, (2.10)
where ∆η = ∆σ + iω∆� and ∆z = iω∆µ represent the variations in electrical properties.
5
The components of the electric and magnetic primary fields are obtained analytically.
In the case of the secondary fields Es and Hs, these can be obtained from the Finite
element method which consists of a numerical technique for calculating approximate so-
lutions of certain differential equations. The method is based on subdividing the problem
domain into smaller elements (or finite elements). The problem conceived in the secondary
fields is subject to homogeneous Dirichlet boundary condition across the mesh border.
2.1 THE MAGNETOTELLURIC METHOD
The magnetotelluric (MT) method is a well-known study procedure among electro-
magnetic methods as a source of information on the geoelectric properties of subsurface
structures. These surveys give knowledge about wide regions and without damage to the
environment. Among the electromagnetic methods in the frequency domain, the MT uses
natural sources, encompassing a frequency range ranging from 10.000 Hz to 0.0001 Hz.
As a starting point for the forward problem, the magnetotelluric method can be ad-
mitted to the one-dimensional model of variation of geoelectric properties, one in which
the earth is assumed as a medium consisting of several horizontal layers (1D).
For a stratified layered environment, the appreciation of the orthogonal components of
the electric field and the magnetic field on the surface, at a certain frequency, allows the
deduction of the apparent resistivity of the subsurface, which corresponds to the resistivity
of a proportional half-space to the stratified medium. This apparent resistivity ρa can be
mathematically measured from the expression
ρa =
1
ωµ
∣∣∣∣EH
∣∣∣∣2
em z=0
, (2.11)
where µ is the magnetic permeability, generally assumed to be the same as the vacuum
(µ = µ0), and ω the angular frequency (ω = 2πf), then this function is dependent on
both the resistivities of the layers in depth and the defined frequency.
The expression 2.11 above can be rewritten and the apparent resistivity modeled as
ρa =
1
ωµ0
∣∣∣Ẑ1∣∣∣2, (2.12)
where Ẑ1 is the apparent impedance on the surface of the layered medium. The apparent
impedance Ẑ1 can be computed based on the following recursive algorithm
Ẑj = Zj
Ẑj+1 + Zj tanh(ikjhj)
Zj + Ẑj+1 tanh(ikjhj)
, (2.13)
where hj corresponds to the thickness of the j-th layer and kj is the wave number.
6
2.2 INVERSE PROBLEM
The inverse problem intends to estimate the values and distribution of a particular
physical property of the subsurface, the parameters of the problem, that correlate the
model data with the observations.
The observation vector y0 has n observations. Each value refers to the observed data
(apparent resistivity and phase angle, respectively, the real and imaginary parts of the
impedance tensor) from a given station at a given frequency. Similarly, the parameter
vector p is composed of m parameters, which are the resistivities of each homogeneous
cell in a grid encompassing the area of interest in the subsurface.
The problem defines a synthetic data vector y, each element of which is associated
with each observation. The functional relationship between parameters and the synthetic
data is then defined as
y = f(x, ω,p), (2.14)
where x refers to the coordinates of the measurement stations and ω represents the angular
frequencies used in the survey.
The inverse problem is based on the search for a parameters vector p̂ that generates
a computed data set that approximates the observation vector. The `2 norm (also known
as the Euclidean norm) is used regularly to measure this proximity condition between the
computed and observed data. This defines the objective function φ(p) as
φ(p) = ||y0 − y||2, (2.15)
φ(p) =
n∑
1
(
y0 − y
)2
, (2.16)
or
φ(p) =
(
y0 − y
)T (
y0 − y
)
. (2.17)
From this, the inversion process is based on the calculation of a vector p̂ to minimize
this functional φ(p). This means looking for a vector p̂ so that the gradient of the
objective function in this vector is null, this being an extreme point of the objective
function. Convergence is evaluated based on the variation of the φ function in consecutive
inversion cycles, that is,
|φk(p)− φk−1(p)| < �, (2.18)
where � refers to the value for interrupting the process.
There are several methods in the literature and used to search for a p̂ vector that
minimizes a function. The often methodology develops as an iterative investigation that
starts from an initial estimate p0 and proceeds with the determination of increments δp,
the steps. These increments are applied to the estimates successively and the process
7
ends when a vector p̂ is found such that the desired approximation is satisfied, that is,
the stopping criterion is satisfied.
Among the most used and known methods, such as the Newton and Gauss-Newton
method, there is the class of methods that start from the same principle as these, the group
of quasi-Newton methods. These procedures start from the premise of the computation
of increments δp from the Taylor series expansion of the function φ(p) to be minimized.
This expansion is assumed until the second order around the first approximation p0, with
p = p0 + δp,
φ(p0 + δp) ≈ φ(p0) + δpT∇φ(p0) + 1
2
δpT∇
[
∇φ(p0)
]
δp, (2.19)
where the gradient of φ, considering the equation 2.17, can be expressed as
∇pφ(p) =
∂
∂pj
[(
y0 − f
)T (
y0 − f
)]
, (2.20)
∇pφ(p) =
∂φ
∂pj
=
[
−∂f
T
∂pj
(
y0 − f
)]
+
[
−
(
y0 − f
)T ∂f
∂pj
]
, (2.21)
or better,
∂φ
∂pj
= −2∂f
T
∂pj
(
y0 − f
)
, (2.22)
where
∂f
∂pj
is a n×m dimension matrix expressed as
∂f
∂pj
=

∂f1
∂p1
∂f1
∂p2
· · · ∂f1
∂pm−1
∂f1
∂pm
∂f2
∂p1
∂f2
∂p2
· · · ∂f2
∂pm−1
∂f2
∂pm
...
... . . .
...
...
∂fn−1
∂p1
∂fn−1
∂p2
· · · ∂fn−1
∂pm−1
∂fn−1
∂pm
∂fn
∂p1
∂fn
∂p2
· · · ∂fn
∂pm−1
∂fn
∂pm

= An×m. (2.23)
This Jacobian matrix of f(x, ω,p) is called sensitivity matrix, since each element of a
line represents its sensitivity of this element in relation to the variation of the parameters.
We can then rewrite the equation 2.22 as
∇pφ(p) = −2AT
(
y0 − f
)
= −2ATc, for c = y0 − f . (2.24)
In the last part of the equation 2.19 we have the term ∇{∇φ(p0)} referring to the
8
Hessian matrix which can be developed as
∂φ
∂pi
[
∂φ
∂pj
]
=
∂φ
∂pi
[
−2∂f
T
∂pj
(
y0 − f
)]
, (2.25)
∂φ
∂pi
[
∂φ
∂pj
]
=
[
−2 ∂
2fT
∂pi∂pj
(
y0 − f
)]
+
[
2
∂fT
∂pj
∂f
∂pi
]
= Hm×m. (2.26)
The Gauss-Newton method and the class of quasi-Newton methods have as their
essence the determination of an approximation for this Hessian matrix (or directly for
its inverse, depending on the procedure).
Evaluated the sensitivity matrix and the Hessian approximation for the function φ(p)
in 2.19 as the second order approximation of φ(p) around p0,we intend to find an ap-
proximate vector p̂ that minimizes the function, as before said, such that its gradient is
equal to zero. Considering, then, the gradient of φ,
∂φ
∂δp
= −2ATc + Hδp = 0. (2.27)
Thus, the system to be solved in the iterative computation of the steps δp to approx-
imate p̂ takes the form
−2ATc + Hδp = 0, (2.28)
δp = H−1
(
2ATc
)
= −H−1∇pφ(p). (2.29)
The new estimate of the parameter vector in an k + 1 iteration is given by
pk+1 = pk + δpk+1 (2.30)
or
pk+1 = pk −H−1∇pφ(p). (2.31)
2.2.1 Regularization
Inverse problems in geophysics usually present themselves as ill-posed, mainly due to
instability and/or lack of uniqueness in solutions. This can be caused both by the presence
of noise in the observed data set and by the essence of the inverse geophysical problem
that, in general, is based on the determination of a set of parameters from a limited
amount of information. To reduce the obstacles in the inverse problem and achieve a
stable solution by imposing a priori constraints, the regularization process is used.
The inverse problem regularized seeks to estimate the parameters to both adjust the
observed data and satisfy the imposed constraints. The constraints imposed introduce
initial information that may be of a purely mathematical or geophysical nature. This
information is added to the problem from functions dependent on certain variables. With
9
this addition of information, the objective function φ takes a new form, according to the
work of Tikhonov (1963), as
φ(p, µ) = φd(p) + µφv(p), (2.32)
where φd(p) is defined by 2.15, µ is the regularization parameter and φv(p) refers to the
regularization relationship imposed to the solution.
The µ regularization parameter can be defined in several ways, depending on the order
of the data set. The contribution of the regularization parameter µ is based on controlling
the relative relevance of the φv(p) function. It is necessary to choose a value for µ that
favors a satisfactory adjustment of the data and make the solution stable.
The most common practice is to use smoothing constraints that impose proximity
between adjacent parameters, that is, a global smoothing function. The application of
this function is appropriate when the computation of the parameters must be tenuous.
The use of this function is based on the understanding that contiguous parameters must
be computed with close values, that is, they must not have values changing abruptly
or that the contrast between neighboring parameters must be as small as possible. The
vector of the differences between the close parameters can be expressed as
v =

p2 − p1
p3 − p2
p4 − p3
...
pn − pn−1
 =

−1 1 0 · · · 0 0 0
0 −1 1 · · · 0 0 0
0 0 −1 · · · 0 0 0
...
...
... . . .
...
...
...
0 0 0 · · · −1 1 0
0 0 0 · · · 0 −1 1


p1
p2
p3
...
pn−1
pn

= Mp. (2.33)
Therefore, the φv(p) function is expressed as
φv(p) = ||v||2 = ||Mp||2 = pTMTMp. (2.34)
In addition, the gradient φ′v(p) of this function and the Hessian φ′′v(p) can be written,
respectively, as
φ′v(p) = 2M
TMp (2.35)
and
φ′′v(p) = M
TMp. (2.36)
The relationship between adjacent parameters in the application of the association can
be established in more than one direction in the two-dimensional case. For these models,
the bond between neighboring parameters (or reference cells) can be vertical, horizontal
or vertical and horizontal, as illustrated in the schemes below
10
Figure 2.2: Schemes illustrating the directions of associations between the parameters in the
two-dimensional case
In addition, when evaluating the number of associations established between adjacent
cells in one or another orientation, with Cv being the number of cells vertically and Ch
horizontally, we have Vv = (Cv − 1)Ch e Vh = (Ch − 1)Cv, with Vv referring to the total
of vertical relations and Vh horizontally.
2.2.2 Gauss-Newton method
The Gauss-Newton method is structured to minimize the objective function based on
the calculation of the Jacobian and Hessian matrices. For the method, the term involving
the second derivatives of f(x, ω,p) in the expression 2.26 is ignored, marking the difference
of the Newton with this one, remaining the following approximation for the hessian
∂φ
∂pi
[
∂φ
∂pj
]
≈ 2∂f
T
∂pj
∂f
∂pi
= 2ATA ≈ Hm×m. (2.37)
The approximation in the hessian calculation considers the significant contribution of
the term 2
∂fT
∂pj
∂f
∂pi
and avoids costly computation of the cross derivatives in−2 ∂
2fT
∂pi∂pj
(
y0 − f
)
,
which has smaller and smaller values the closer it gets to adjusting the function and has
considerably less influence.
The equation for solving the inverse problem using the Gauss-Newton method, in the
computation of the steps δp can be written as
−2ATc + 2ATAδp = 0, (2.38)
or better,
2ATAδp = 2ATc, (2.39)
with the inverse problem regularized equation being
[
2ATA + µ
(
2MTM
)]
δp = 2ATc− µ
(
2MTMp
)
, (2.40)
[
ATA + µ
(
MTM
)]
δp = ATc− µ
(
MTM
)
p, (2.41)
11
soon,
δp =
[
ATA + µ
(
MTM
)]−1 [
ATc− µ
(
MTMp
)]
. (2.42)
The Gauss-Newton method may diverge, however, if the quadratic approximation
for φ from δpk+1 is not adequate, where the evaluated Hessian matrix it may not be
positive definite. The strategy used to avoid the problem and ensure convergence is the
Levenberg–Marquardt method. For this, the Levenberg–Marquardt method introduces
an addition to the Gauss-Newton method.
To guarantee an acceptable direction for approximations δpk+1 beyond the validity
zone of the quadratic convergence, the Levenberg–Marquardt strategy changes the inverse
of the hessian so that it is a positive definite matrix, not singular, making it strongly
diagonal-dominant and with positive diagonal elements. For this, a term is introduced as
follows to the diagonal of the hessian
δp =
[
ATA + µ
(
MTM
)
+ λI
]−1 [
ATc− µ
(
MTMp
)]
. (2.43)
The value for the λ parameter can increase or decrease in each iteration of the inversion
process, depending on the approximation for the Hessian matrix. The value of λ thus
moderates the step size towards the evaluated gradient.
2.2.3 Quasi-Newton method
The first method known as quasi-Newton appeared in the mid 1950s. The algorithm
emerged as a solution to a problem of computational performance and was initially devel-
oped by the physics and mathematics professor William Cooper Davidon, at the Argonne
National Laboratory.
Fletcher and Powell later showed that the new technique, now known as the quasi-
Newton DFP method, or Davidon–Fletcher–Powell formula, was more reliable and faster
than the existing algorithms at the time, which made the idea a milestone in development
of nonlinear optimization problems.
Quasi-Newton methods only need to calculate the gradient of the objective function
in each iteration of the inversion process. According to Nocedal and Wright (1999), quasi-
Newton methods can be more efficient than Newton’s methods and now they are presented
in a plurality of algorithms for application to several large-scale optimization problems.
Among the quasi-Newton methods, the most popular algorithm is the BFGS method,
named by the authors Broyden, Fletcher, Goldfarb and Shanno. In the quasi-Newton
method BFGS, the expression for the inverse of the approximate Hessian matrix is deter-
mined so that there is no need to solve a system for calculating δpk+1, in the equation
2.29, as in the case of the Gauss-Newton method shown, thus reducing the overall compu-
tational cost, in the case of the total number of operations in each iteration. The iterative
12
search expressed in the equation 2.31 is similar to the Gauss-Newton method. The main
difference is the approximate form for the Hessian matrix, or better, for its inverse.
In an k+1 inversion iteration, the expression in the BFGS method for an approximation
of the Hessian inverse is as follows
H−1k+1 = (I− ρksky
T
k )H
−1
k (I− ρkyks
T
k ) + ρksks
T
k , (2.44)
in which
ρk =
1
yTk sk
, (2.45)sk = pk+1 − pk, (2.46)
yk = ∇φk+1 −∇φk. (2.47)
The expression 2.44 is based on the so-called secant equation (Nocedal and Wright,
1999), which can be determined as follows
H−1k+1 = (H
−1
k − ρkH
−1
k sky
T
k )(I− ρkyksTk ) + ρksksTk , (2.48)
H−1k+1 = H
−1
k − ρkH
−1
k yks
T
k − ρkH−1k sky
T
k + ρ
2
kH
−1
k sky
T
k yks
T
k + ρksks
T
k , (2.49)
H−1k+1yk = H
−1
k yk−ρkH
−1
k yks
T
k yk−ρkH−1k sky
T
k yk+ρ
2
kH
−1
k sky
T
k yks
T
k yk+ρksks
T
k yk, (2.50)
H−1k+1yk = H
−1
k yk − ρkH
−1
k yk
1
ρk
− ρkH−1k sky
T
k yk + ρ
2
kH
−1
k sky
T
k yk
1
ρk
+ ρksk
1
ρk
, (2.51)
H−1k+1yk = H
−1
k yk −H
−1
k yk − ρkH
−1
k sky
T
k yk + ρkH
−1
k sky
T
k yk + sk, (2.52)
H−1k+1yk = sk, (2.53)
or also,
Hk+1sk = yk. (2.54)
The secant equation requires that the matrix Hk+1 be positive definite, possible when
a certain variation vector sk and the vector of difference of the gradients yk satisfy the
curvature condition
sTk yk > 0. (2.55)
This condition, however, will not always be guaranteed for non-convex functions and
in these situations it is necessary to apply certain restrictions to impose 2.55, such as the
Wolfe conditions. From the imposition of these conditions over the line search strategy, a
step size can be determined, that is, a choice for α (Nocedal and Wright, 1999). Thus, a
new estimate of the parameters vector can be expressed as
pk+1 = pk + αkδpk+1 (2.56)
13
or also
pk+1 = pk − αkH−1∇pφ(p). (2.57)
In determining α, we want a step that provides a decrease in function, that is,
φ(pk+1) = φ(pk + αkδpk+1) < φ(pk). However, we may experience very small reduc-
tions in the value of φ relative to the size of this step. Therefore, you must include a
condition that we have not only φ(pk + αkδpk+1) < φ(pk), but also that the value of the
function decreases at least a predetermined fraction from the previous one. To prevent
φ from having a very small decrease from its previous value, we can demand that the
reduction in φ be proportional to the size of αk and the term ∇φ(pk)T δpk+1. So, we have
the first Wolfe condition, also called Armijo condition or condition of sufficient decrease,
φ(pk + αkδpk+1) ≤ φ(pk) + c1αk∇φ(pk)T δpk+1, (2.58)
with 0 < c1 < 1. However, this condition is not enough to guarantee the efficiency of the
strategy, as we can see that it is satisfied for all sufficiently small values of αk. To avoid
this problem, with extremely small steps, the following condition is added that evaluates
the derivative ∇φ(pk + αkδpk+1)T δpk+1 = ∇φ(pk+1)T δpk+1,
∇φ(pk+1)T δpk+1 ≥ c2∇φ(pk)T δpk+1. (2.59)
This condition is due to Wolfe, according to Nocedal and Wright (1999), with c2 ∈ ]c1, 1[.
The guarantee of 2.55 based on Wolfe conditions can be observed from 2.59, 2.46 and
2.47, since (
yTk +∇φ(pk)T
)
sk ≥ c2∇φ(pk)T sk, (2.60)
yTk sk ≥ c2∇φ(pk)T sk −∇φ(pk)T sk, (2.61)
yTk sk ≥ (c2 − 1)αk∇φ(pk)T δpk+1. (2.62)
Therefore, the term on the right-hand side of the expression 2.62 will always be positive,
as long as c1 < c2 < 1 and δpk+1 is the descent direction, ensuring thus the curvature
condition.
2.2.4 Limited-memory BFGS (LM-BFGS)
Quasi-Newton methods should not be directly applicable to large-scale problems, ac-
cording to Nocedal and Wright (1999), as their approximations to the Hessian matrix (or
inverse of the Hessian) are, in general, dense, demanding computation machines with large
requirements. However, it is possible to change and extend the definitions and strategies
of quasi-Newton methods in several ways to make them suitable for large-scale inverse
problems.
14
In the work of Nocedal (1980), a limited memory approach is presented to get approx-
imations of the Hessian matrices stored compactly, from some vectors imposed depending
on the computational capacity. Despite relatively modest memory requirements, this
treatment generally produces an acceptable convergence.
The main idea of the limited memory methodology, according to Nocedal (1980),
is to use curvature information only from the most recent iterations to construct the
approximations for the Hessian inverse. Curvature information from older iterations,
which is less likely to be relevant to the actual behavior of the inverse hessian matrix in
the current iteration, is discarded for saving storage.
Each step of the LM-BFGS method takes the form following the expression 2.57 for
the complete BFGS, where αk+1 is the length of the step, determined in the line search
procedure following strong Wolfe conditions. As already shown, the inverse of the Hessian
matrix H−1k+1 is updated according to the formula
H−1k+1 = V
T
kH
−1
k Vk + ρksks
T
k , (2.63)
where ρk, sk and yk are defined, respectively, by the expressions 2.45, 2.46 and 2.47. Vk
is expressed as
Vk = I− ρkyksTk , (2.64)
where the matrix H−1k+1 is obtained by updating H
−1
k with the pairs {sk,yk}. To get
around the problem of the cost of storing and handling a dense hessian inverse, a number
of the vector pairs {si,yi} are stored for use in a modified version of expression for H−1k .
The product defined by H−1k ∇φk, in particular, can be determined from a sequence of dot
products and vectors addition involving the vector pairs {si,yi}. After a new iteration of
the inversion, the oldest vector pair of the pairs {si,yi} is forgotten and replaced by a new
pair {sk,yk}. To show this strategy, initially establishing an estimate for the parameters
set x0, we have determined g0 and H−10 . It is important to note that for the application
of the LM-BFGS algorithm the choice of an first approximation for the inverse hessian
matrix H−10 must be defined. A possible solution to be used for the demand for the
first Hessian inverse H−10 , presented by Nocedal and Wright (1999), is to start with the
approximation by the expression H−10 = β||g0||−1I, with g0 being the gradient ∇pφ(p0)
and β a multiple defined a priori.
For k = 0, we have calculated
x1 = x0 − α0H−10 g0, g1, (2.65)
s0 = x1 − x0, (2.66)
y0 = g1 − g0, (2.67)
15
H−11 = (I− ρ0s0yT0 )H−10 (I− ρ0y0sT0 ) + ρ0s0sT0 , (2.68)
or better,
H−11 = V
T
0 H
−1
0 V0 + ρ0s0s
T
0 . (2.69)
Then, for k = 1,
x2 = x1 − α1H−11 g1, g2, (2.70)
s1 = x2 − x1, (2.71)
y1 = g2 − g1, (2.72)
H−12 = V
T
1 H
−1
1 V1 + ρ1s1s
T
1 , (2.73)
in other words,
H−12 = V
T
1
(
VT0 H
−1
0 V0 + ρ0s0s
T
0
)
V1 + ρ1s1s
T
1 , (2.74)
H−12 = V
T
1 V
T
0 H
−1
0 V0V1 + V
T
1 ρ0s0s
T
0V1 + ρ1s1s
T
1 . (2.75)
The term H−10 can be updated in the iterations considering the curvature information.
A methodology defined and demonstrated by Nocedal and Wright (1999) for calculating
H−1k0 is given by H
−1
k0
= γkI, where
γk =
sTk−1yk−1
yTk−1yk−1
. (2.76)
This term would be a scaling factor related to the computation of the Hessian matrix
along the most recent descent direction. Thus succeeding for k = 2,
x3 = x2 − α2H−12 g2, g3, (2.77)
s2 = x3 − x2, (2.78)
y2 = g3 − g2, (2.79)
H−13 = V
T
2 H
−1
2 V2 + ρ2s2s
T
2 , (2.80)
H−13 = V
T
2 V
T
1 V
T
0 H
−1
k0
V0V1V2 + V
T
2 V
T
1 ρ0s0s
T
0V1V2 + V
T
2 ρ1s1s
T
1V2 + ρ2s2s
T
2 . (2.81)
Then, following the strategy of keeping the pairs {sk,yk} of the inversion, we make
ρ0 = s0 = y0 = 0, then V0 = (I− ρ0s0yT0 ) = I and
H−13 = V
T
2 V
T
1 H
−1
k0
V1V2 + V
T
2 ρ1s1s
T
1V2 + ρ2s2s
T
2 , (2.82)
as well as,
H−14 = V
T
3 V
T
2 H
−1
k0
V2V3 + V
T
3 ρ2s2s
T
2V3 + ρ3s3s
T
3 , (2.83)
H−15 = V
T
4 V
T
3 H
−1
k0
V3V4 + V
T
4 ρ3s3s
T
3V4 + ρ4s4s
T
4 , (2.84)
16
and so on. From these deductions we can generalize the LM-BFGS to k + 1 ≤ m
H−1k+1 =
(
VTkV
T
k−1...V
T
0
)
H−10 (V0...Vk−1Vk) (2.85)
+ VTk ...V
T
1 ρ0s0s
T
0V1...Vk
+ VTk ...V
T
2 ρ1s1s
T
1V2...Vk
+ ... +
+ VTkV
T
k−1ρk−2sk−2s
T
k−2Vk−1Vk
+ VTk ρk−1sk−1s
T
k−1Vk
+ ρksks
T
k ,
and for k + 1 > m,
H−1k+1 =
(
VTk ...V
T
k−m+1
)
H−10 (Vk−m+1...Vk) (2.86)
+ VTk ...V
T
k−m+2ρk−m+1sk−m+1s
T
k−m+1Vk−m+2...Vk
+ VTk ...V
T
k−m+3ρk−m+2sk−m+2s
T
k−m+2Vk−m+3...Vk
+ ... +
+ VTk ρk−1sk−1s
T
k−1Vk
+ ρksks
T
k
For m = k + 1 the expressions are analogous. These expressedmatrices are also called
special BFGS matrices.
The product H−1k ∇φk can be recursively computed, following the algorithm presented
by Nocedal and Wright (1999). To deduce the algorithm, suppose k = 3, where we already
have H−14 (equation 2.83) and m = 2, which represents the number of pairs {sk,yk} of
the most recent stored iterations (or “memory size”) for use. Therefore, in particular,
(H−14 g4) = V
T
3 V
T
2 H
−1
0 V2V3g4 + V
T
3 ρ2s2s
T
2V3g4 + ρ3s3s
T
3 g4, (2.87)
(H−14 g4) = (I− ρ3s3yT3 )(I− ρ2s2yT2 )H−10 (I− ρ2y2sT2 )(I− ρ3y3sT3 )g4 (2.88)
+ (I− ρ3s3yT3 )ρ2s2sT2 (I− ρ3y3sT3 )g4 + ρ3s3sT3 g4,
(H−14 g4) = (I− ρ3s3yT3 )(I− ρ2s2yT2 )H−10 (I− ρ2y2sT2 )(g4 − y3ρ3sT3 g4) (2.89)
+ (I− ρ3s3yT3 )ρ2s2sT2 (g4 − y3ρ3sT3 g4) + s3ρ3sT3 g4.
Making γ3 = ρ3sT3 g4 and q3 = g4 − y3γ3,
(H−14 g4) = (I−ρ3s3yT3 )(I−ρ2s2yT2 )H−10 (I−ρ2y2sT2 )q3+(I−ρ3s3yT3 )ρ2s2sT2 q3+s3γ3, (2.90)
17
(H−14 g4) = (I− ρ3s3yT3 )(I− ρ2s2yT2 )H−10 (q3 − ρ2y2sT2 q3) + (I− ρ3s3yT3 )s2ρ2sT2 q3 + s3γ3.
(2.91)
Now assuming γ2 = ρ2sT2 q3 and q2 = q3 − y2γ2,
(H−14 g4) = (I− ρ3s3yT3 )(I− ρ2s2yT2 )H−10 q2 + (I− ρ3s3yT3 )s2γ2 + s3γ3. (2.92)
Computed all γi and qi, for i = 3, 2, r2 = H−10 q2, so
(H−14 g4) = (I− s3ρ3yT3 )(I− s2ρ2yT2 )r2 + (I− s3ρ3yT3 )s2γ2 + s3γ3, (2.93)
(H−14 g4) = (I− s3ρ3yT3 )
[
(I− s2ρ2yT2 )r2 + s2γ2
]
+ s3γ3, (2.94)
(H−14 g4) = (I− s3ρ3yT3 )
(
r2 − s2ρ2yT2 r2 + s2γ2
)
+ s3γ3. (2.95)
Then taking β2 = ρ2yT2 r2,
(H−14 g4) = (I− s3ρ3yT3 ) (r2 − s2β2 + s2γ2) + s3γ3, (2.96)
(H−14 g4) = (I− s3ρ3yT3 ) [r2 + s2(γ2 − β2)] + s3γ3. (2.97)
With r3 = r2 + s2(γ2 − β2), it follows that
(H−14 g4) = (I− s3ρ3yT3 )r3 + s3γ3, (2.98)
and
(H−14 g4) = r3 − s3ρ3yT3 r3 + s3γ3. (2.99)
Lastly, with β3 = ρ3yT3 r3, we have
(H−14 g4) = r3 − s3β3 + s3γ3, (2.100)
(H−14 g4) = r3 + s3(γ3 − β3). (2.101)
After calculating all rj and βj, with j = 2, 3, we finally have the product (H−14 g4)
desired from dot products only. We can then generalize the product (H−1k+1gk+1) with a
recursive scheme starting with the computation of the values γi and the vector q (with
i = k, ..., k −m+ 1), where
γi = ρis
T
i q (2.102)
and
q = q− yiγi, (2.103)
with the initial value for q = gk+1. After that, the calculation for βj and r follows, with
18
the first r = H−10 q, where j = k −m+ 1, ..., k, to
βj = ρjy
T
j r (2.104)
and
r = r + sj(γj − βj), (2.105)
where H−1k+1gk+1 = rlast.
Considering the limited memory logic, which is based only on vectors determined
from the number of parameters, we can say that the greater the number of parameters
considered, the greater the difference in memory required in the inverse problem, as is
demonstrated by the analysis of the models in the results chapter.
2.3 CALCULATION OF SENSITIVITIES
With the global matrix and the source vector for the solution of the system and
definition of the vector referring to the secondary field, it is also possible to find the
sensitivities of the field, that is, the derivatives related to the parameters and associated
with the positions of the receivers (stations). Particularly in the case of TE mode, G
being the global matrix and f the source vector, built from the contributions of each
element and its nodes, the system solution is the secondary electric field at each node and
the system can be written
GEsy = f . (2.106)
To calculate the sensitivity of a cell (that represents a pk parameter) in relation to the
receivers positions, we start by deriving (Meunier, 2008) the matrix equation 2.106
∂
∂pk
[
GEsy = f
]
, (2.107)
∂G
∂pk
Esy + G
∂Esy
∂pk
=
∂f
∂pk
, (2.108)
G
∂Esy
∂pk
=
∂f
∂pk
−G
∂Esy
∂pk
, (2.109)
then,
∂Esy
∂pk
= G−1
(
∂f
∂pk
−G
∂Esy
∂pk
)
. (2.110)
The equation 2.110 returns the sensitivities in all Nn nodes of the mesh. In the interest
of obtaining the sensitivities only in the Ne receiver positions, we can use a matrix Q,
with dimensions of Ne rows and Nn columns , where each line has a value of 1 in the place
19
that corresponds to the node of interest and the others are admitted as 0, so
Q
∂Esy
∂pk
= QG−1
(
∂f
∂pk
−G
∂Esy
∂pk
)
. (2.111)
In the electromagnetic geophysical problems it is frequent that the number of param-
eters Np is greater than the number of stations Ne. Aiming to reduce the number of
operations associated with the solutions of the developed systems, we can work with the
transposed in 2.111 and get(
Q
∂Esy
∂pk
)T
=
(
∂f
∂pk
−G
∂Esy
∂pk
)T
G−1QT , (2.112)
thus, we perform only Ne system solutions for each factoring, in place of Np solutions.
With the derivatives of the fields, the sensitivities of the apparent resistivities ρa and
phases φ are determined following the relations presented by Farquharson and Oldenburg
(1996)
∂ρa
∂ρk
= 2ρa
[
<
(
1
E
∂E
∂ρk
)
−<
(
1
H
∂H
∂ρk
)]
, (2.113)
∂φ
∂ρk
= =
(
1
E
∂E
∂ρk
)
−=
(
1
H
∂H
∂ρk
)
. (2.114)
In the inversion process of the parameters values, in particular, the log scale was used,
the natural logarithm of the properties and apparent resistivity with the main purpose
being to reduce the disparities between their values. That said, in computation of the
sensibilities of ρa and phases φ we have
∂`n(ρa)
∂`n(ρk)
=
1
ρa
ρk
∂ρa
∂ρk
, (2.115)
∂φ
∂`n(ρk)
= ρk
∂φ
∂ρk
, (2.116)
therefore
∂`n(ρa)
∂`n(ρk)
= 2ρk
[
<
(
1
E
∂E
∂ρk
)
−<
(
1
H
∂H
∂ρk
)]
, (2.117)
∂φ
∂`n(ρk)
= ρk
[
=
(
1
E
∂E
∂ρk
)
−=
(
1
H
∂H
∂ρk
)]
. (2.118)
2.4 SPACIAL DERIVATIVES
To simulate the response of the magnetotelluric propagation modes to two-dimensional
earth topography, taking Wannamaker et al. (1986) as a reference, care must be taken as
to the nodal values in the calculation of spatial derivatives. To calculate the solutions of
the nodal fields parallel to the strike (Ex for the TE mode and Hx for the TM mode), the
20
auxiliary fields transversal to the strike are obtained through a nodal values numerical
differentiation. The nodal values of fields can be taken within the uniform earth media
in two ways. Only in one direction, as under a flat surface, horizontally (gray dots
- alternative A1 in the figure 2.3) or parallel to the topographic surface (black dots -
alternative A2). As in Wannamaker et al. (1986), placing the sensors horizontally reduces
topographic anomalies compared to locating sensors parallel to the slope.
Figure 2.3: Topographic surface model, a ramp.
3 RESULTS - APPLICATION TO THE INVERSION OF
MAGNETOTELLURIC DATA
The LM-BFGS method was applied to the inversion of magnetoteluric synthetic data
in two-dimensional media, with the primary goal of analyzing its feasibility and efficiency.
The results are compared with those obtained with the Gauss-Newton method.
MT data were simulated using finite elements in 2D models, with the addition of
uniformly distributed zero mean noise. All the code is in FORTRAN.
Using the smoothing constraint, a value was accepted for the parameter µ = 10−3,
with tolerances for the inversion process � = 10−4 (model 1) and � = 10−5 (model 2) for the
variation of the objective function between two successive iterations. The initial resistivity
values (first step) for the inversions were automatically defined as an average value of the
set of apparent input resistivities, with the initial model, then, being a homogeneous
medium.
In applying the LM-BFGS method, practical experiences have shown that few stored
pairs {sk,yk} (between 3 and 20 pairs) are often enough to produce satisfactory results
(Nocedal, 1980). With this information and based on the evaluations carried out, the
number of 5 {sk,yk} pairs stored from the most recent iterations in the performed inver-
sions was defined.
As an example of the investigations to determine the sufficient number 5 of stored
pairs, figure 3.1 shows results from the inversion of data from a one-dimensional layered
model. In this case the amount of memory used in the forward problem is a small fraction
of the total and the measured memory used can be considered exclusively that of the
inverse problem. The plots in the figure show the numberof iterations, processing time
and memory as functions of the number of pairs stored to carry out the inversion.
This problem presents little demand for memory and takes very little time, even in the
3 4 5 6 7
Storage
120
130
140
150
160
170
180
I
te
r
a
ti
o
n
s
LM −BFGS
3 4 5 6 7
Storage
4.5
5
5.5
6
6.5
7
T
im
e
(s
e
c)
LM −BFGS
3 4 5 6 7
Storage
110
111
112
113
114
115
M
e
m
o
r
y
(k
B
)
LM −BFGS
Gauss−Newton : 23 ite
Gauss−Newton : 245.8 kB
Gauss−Newton : 0.9 sec
Figure 3.1: Example of iteration, time and memory spent by both methods.
21
22
traditional Gauss-Newton approach. However, it provides a sufficient case to show that
with 5 pairs stored there is a certain stabilization in the necessary number of iterations
to reach the desired accuracy of the solution and there is not a considerable variation
in memory for the LM-BFGS. The information of memory, time and iterations of the
Gauss-Newton method are given as a reference for comparison.
For the analysis of the total elements stored in the inverse problem with the increase in
the number of stored pairs, supposing now the inversion of an idealized two-dimensional
model with 19 frequencies and 41 stations, from apparent resistivity and phase, with a
grid of 2000 parameters (80 by 25), the Gauss-Newton method, in general, stores elements
in the hessian and in the sensitivity matrix, this number being dependent on both the
number of parameters and the number of observations of the problem. In the LM-BFGS
method, this number depends only on the number of parameters of the problem and on the
number of stored pairs of the most recent iterations that is desired, as a linear dependency,
as shown in the following figure. This figure also shows the number of iterations with the
use of around 5 stored pairs.
3 4 5 6
Storage
1.2
1.4
1.6
1.8
2
2.2
2.4
N
u
m
be
r
o
f
e
le
m
e
n
ts
×10
4 LM − BFGS
3 4 5 6
Storage
58
58.5
59
59.5
60
I
te
r
a
ti
o
n
s
LM − BFGS
Gauss−Newton : 7116000
Gauss−Newton : 26
Figure 3.2: Values stored and iterations by both methods for 2D inversion.
All results obtained and the evaluations were carried out in a Workstation Silix®12
cores coffee lake intel®core, I7-8700 3.2 GHZ 12 MB, 64GB DDR4 2400, SSD 480GB, 1
TB SATA3, GF GTX 1060 6GB 1280 CUDA.
3.1 MODEL 1
The model consists of a square body with 1 km in side, with resistivity ρ1 = 10 Ωm, 3
km deep (taking as reference the highest point of the surface), imbeded in a homogeneous
medium with ρ2 = 100 Ωm. The topography consists of a slope with 3 km horizontal
length and 1.2 km in the vertical direction. Along this surface, a survey line was modeled
with 28 stations in the interval from -5 km to 5 km. The data was simulated at 37
frequencies, ranging from 10−3 to 103 Hz.
23
Figure 3.3: Model 1 - Ramp and rectangular body.
The inversion grid was defined following the topography variation. The cells, normally
rectangular for plane surfaces, are deformed to accommodate the relief. At the greatest
depths they adjust to the lower plane limit of the grid. There are 2720 cells (80 × 34)
125 m in width, with thickness of 100 m at the top and increasing with depth at a rate
of 5%, with a 10 km limit.
The LM-BFGS inversion converged to a solution with 78 iterations in more than 149
minutes. For the same data set, the inversion with Gauss-Newton reached the same
solution in 28 iterations and over 49 minutes. Figure 3.4 shows the absolute values of the
variation of the function φ in successive iterations.
The LM-BFGS achieved the same results as the Gauss-Newton method, as can be
seen in figures 3.5 and 3.6. The same data are shown in figure 3.20 with contour lines.
The Gauss-Newton code needs to store the Hessian matrix (size 2720 × 2720) as well
as the sensitivity matrix (2072 × 2720) for a total of 13 034 240 values. Together with
the requirements of the forward problem, this results in a total of 221.04 MB of memory
used.
In an iteration of the LM-BFGS program, only the 5 {sk,yk} pairs of vectors from
the previous 5 iterations are stored, so that each iteration stores 5 × 2 × 2720, or 27200
values, only 0.2% of the amount stored in the Gauss-Newton method. The total memory
24
0 10 20 30 40 50 60 70 80
Iterations
10
-6
10
-4
10
-2
10
0
10
2
|φ
k
+
1
−
φ
k
|
Convergence
LM − BFGS
Gauss−Newton
Figure 3.4: The absolute values of successive φ differences in the inversion of model 1.
Tabela 3.1: The main differences between the inversions of the model 1
Method Stored units Number of itera-tions Processing time
Gauss-Newton 13034240 28 49 min
LM-BFGS quasi-
Newton 27200 78 149 min
used by the LM-BFGS program reaches a maximum of 58.4 MB, which is about 26,5% of
the total memory requirement of the Gauss-Newton program.
Table 3.1 summarizes these differences between the two methodologies, in terms of
stored units, number of iterations and processing time.
The data fit reached by both methods is shown figures 3.8 and 3.9 for the observations
from the first and last stations on the survey line.
In this comparison, the differences in memory usage and processing time are all in the
inverse problem. The total memory used in the forward problem is the same for both
methods.
Apparent resistivity and phase pseudosections estimated by Gauss-Newton and LM-
BFGS inversions are shown in the figures 3.11, 3.12, 3.14 and 3.15. The observation
pseudosections of apparent resistivity and phase are the figures 3.10 and 3.13.
25
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
(Gauss−Newton) Inverted data from ρa and φ − TE mode − MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
Figure 3.5: Result of the model 1 inversion with Gauss-Newton for a grid with 2720 cells. The
true target body position is shown by the white line.
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
(LM −BFGS) Inverted data from ρa and φ − TE mode − MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
Figure 3.6: Results of the model 1 inversion with the LM-BFGS quasi-Newton method for a 2720
cells grid. The true target body position is shown by the white line.
26
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
(Gauss−Newton) Inverted data from ρa and φ − TE MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
(a) Model 1 Gauss-Newton inversion
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
(LM −BFGS) Inverted data from ρa and φ − TE MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
(b) Model 1 LM-BFGS inversion
Figure 3.7: Model 1 inversions in contour lines
27
80
85
90
95
100
105
110
ρ
a
(Ω
·
m
)
Observation
LM −BFGS
Observation
Gauss−Newton
0.001 0.01 0.1 1 10 100 
T (s)
40
45
50
φ
(◦
)
0.01 0.1 1 10 100 1000
T (s)
Figure 3.8: Observed and estimated data for the first station in model 1.
90
95
100
105
110
ρ
a
(Ω
·
m
)
Observation
LM −BFGS
Observation
Gauss−Newton
0.001 0.01 0.1 1 10 100 
T (s)
40
45
50
φ
(◦
)
0.01 0.1 1 10 100 1000
T (s)
Figure 3.9: Observed and estimated data for the last station in model 1.
28
(Observation) Apparent Resistivity
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
Distance (m)
80
85
90
95
100
105
110
115
120
ρa
Figure 3.10: Apparent resistivity pseudosection (Model 1 - Observation).
(Gauss−Newton) Apparent Resistivity
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
Distance (m)
80
85
90
95
100
105
110
115
120
ρa
Figure 3.11: Apparent resistivity pseudosection (Model 1 - Gauss-Newton).
(LM −BFGS) Apparent Resistivity1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
Distance (m)
80
85
90
95
100
105
110
115
120
ρa
Figure 3.12: Apparent resistivity pseudosection (Model 1 - LM-BFGS).
29
(Observation) Phase
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
Distance (m)
43
44
45
46
47
48
49
50
φ
Figure 3.13: Phase angle pseudosection (Model 1 - Observation).
(Gauss−Newton) Phase
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
Distance (m)
43
44
45
46
47
48
49
φ
Figure 3.14: Phase angle pseudosection (Model 1 - Gauss-Newton).
(LM −BFGS) Phase
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000
Distance (m)
43
44
45
46
47
48
49
φ
Figure 3.15: Phase angle pseudosection (Model 1 - LM-BFGS).
30
3.2 MODEL 2
To verify the consistency of the comparison results, a second model was defined with
greater variations in topography and three target bodies, as shown in the figure 3.16.
The three bodies have dimensions of 1 km by 1 km and resistivities ρ1 = 10 Ωm,
positioned directly under the topographic features, which are typified as a slope, a hill
and a valley. The background medium has resistivity ρ2 = 100 Ωm. As shown in figure
3.16, the bodies are located at a depth of 3 km below the highest point of the topography.
Along this topographic surface, 79 measurement stations were simulated in the interval
from -12 km to 12 km. 31 frequencies ranging from 10−3 to 102 Hz were used in each
station.
Figure 3.16: Model 2 - Medium with three target bodies and complex topographic relief.
The inversion grid was defined with 2112 cells (96 × 22) 250 m wide, with initial
thickness of 200 m and increasing in depth with a growth rate of 5% and a 10 km depth
limit. The successive φ variations between iterations are shown in figure 3.17.
Table 3.2 summarizes the quantitative differences between the two inversions. The
Gauss-Newton inversion code (217.68 MB) converged with 30 iterations in 96 minutes,
and required the storage of 2112 × 2112 values from the Hessian matrix and 2112 × 4898
from the Jacobian matrix, for a total of 14805120 values. The quasi-Newton LM-BFGS
program (52.95 MB total) met the stopping criterion in 415 minutes and 110 iterations,
31
0 20 40 60 80 100 120
Iterations
10
-6
10
-4
10
-2
10
0
10
2
|φ
k
+
1
−
φ
k
|
Convergence
LM − BFGS
Gauss−Newton
Figure 3.17: The absolute values of successive φ differences in the inversion of model 2.
Tabela 3.2: The main differences between the inversions of the model 2
Method Stored units Number of itera-tions Processing time
Gauss-Newton 14805120 30 96.8 min
LM-BFGS quasi-
Newton 21120 110 415.7 min
with 5× 2× 2112 = 21120 stored values per iteration. The LM-BFGS code requires only
about 0.14% of the amount of storage for the inverse problem and about 24% of the total
memory needed by the Gauss-Newton program.
Figure 3.17 presents the values of |φk(p)− φk−1(p)| along the iterations.
Figures 3.18 and 3.19 show the resulting resistivity distribution achieved by both
methods. Figure 3.20 shows the same results using contour lines.
Figures 3.21 and 3.22 present the data fit computed by the methods, for apparent
resistivities and phases.
For this model, the apparent resistivity and phase observation pseudosections are
shown in figures 3.23 and 3.26. The pseudosections estimated by the inversions are the
figures 3.24, 3.25, 3.27 and 3.28.
32
-1 -0.5 0 0.5 1
×10
4
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
(Gauss−Newton) Inverted data from ρa and φ − TE mode − MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
Figure 3.18: Model 2 inversion with Gauss-Newton for a grid with 2112 cells
-1 -0.5 0 0.5 1
×10
4
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
(LM −BFGS) Inverted data from ρa and φ − TE mode − MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
Figure 3.19: Model 2 inversion with LM-BFGS for a 2112 cells grid
33
-1 -0.5 0 0.5 1
×10
4
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
(Gauss−Newton) Inverted data from ρa and φ − TE MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
(a) Model 2 Gauss-Newton inversion
-1 -0.5 0 0.5 1
×10
4
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
(LM −BFGS) Inverted data from ρa and φ − TE MT
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
log10ρ
(b) Model 2 LM-BFGS inversion
Figure 3.20: Model 2 inversions in contour lines
34
85
90
95
100
105
110
ρ
a
(Ω
·
m
)
Observation
LM − BFGS
Observation
Gauss−Newton
0.01 0.1 1 10 100 
T (s)
40
45
50
φ
(◦
)
0.01 0.1 1 10 100 1000
T (s)
Figure 3.21: Observation and calculated data for the first station (model 2)
85
90
95
100
105
110
115
ρ
a
(Ω
·
m
)
Observation
LM − BFGS
Observation
Gauss−Newton
0.01 0.1 1 10 100 
T (s)
40
45
50
φ
(◦
)
0.01 0.1 1 10 100 1000
T (s)
Figure 3.22: Observation and calculated data for the last station (model 2)
35
(Observation) Apparent Resistivity
1000
100 
10 
1 
0.1 
0.01
P
e
r
io
d
(s
)
-1 -0.5 0 0.5 1
Distance (m) ×104
70
80
90
100
110
120
130
ρa
Figure 3.23: Apparent resistivity pseudosection (Model 2 - Observation).
(Gauss−Newton) Apparent Resistivity
1000
100 
10 
1 
0.1 
0.01
P
e
r
io
d
(s
)
-1 -0.5 0 0.5 1
Distance (m) ×104
70
80
90
100
110
120
130
ρa
Figure 3.24: Apparent resistivity pseudosection (Model 2 - Gauss-Newton).
(LM −BFGS) Apparent Resistivity
1000
100 
10 
1 
0.1 
0.01
P
e
r
io
d
(s
)
-1 -0.5 0 0.5 1
Distance (m) ×104
70
80
90
100
110
120
130
ρa
Figure 3.25: Apparent resistivity pseudosection (Model 2 - LM-BFGS).
36
(Observation) Phase
1000
100 
10 
1 
0.1 
0.01
P
e
r
io
d
(s
)
-1 -0.5 0 0.5 1
Distance (m) ×104
41
42
43
44
45
46
47
48
49
50
51
φ
Figure 3.26: Phase angle pseudosection (Model 2 - Observation).
(Gauss−Newton) Phase
1000
100 
10 
1 
0.1 
0.01
P
e
r
io
d
(s
)
-1 -0.5 0 0.5 1
Distance (m) ×104
41
42
43
44
45
46
47
48
49
50
φ
Figure 3.27: Phase angle pseudosection (Model 2 - Gauss-Newton).
(LM −BFGS) Phase
1000
100 
10 
1 
0.1 
0.01
P
e
r
io
d
(s
)
-1 -0.5 0 0.5 1
Distance (m) ×104
41
42
43
44
45
46
47
48
49
50
φ
Figure 3.28: Phase angle pseudosection (Model 2 - LM-BFGS).
37
3.3 MODEL 3
Based on Reyes-Wagner et al. (2017), a regional line example was included as a third
model. The model consists of a profile of stations in the South Volcanic Zone (SVZ) of
the Chilean Andes to simply represent the subduction zone and its relationship with the
volcanic arc at this latitude.
The Andean SVZ extends between 33° and 46° S and was formed as the Nazca Plate is
obliquely subducted below the South American Plate and is characterized by many active
volcanoes in the volcanic arc. As in other subduction zones, the fluids released from the
fusion of the subducted plate facilitate partial fusion in the upper mantle, reducing the
melting point of the rocks. The electrical resistivity of rocks in the earth’s crust and
the upper mantle is sensitive to the presence of aqueous fluids and partial melts (Reyes-
Wagner et al., 2017). The profile cuts across the ZVS, extending for 150 km with 25
stations covering a period from 0.001 to 1000s, with 25 periods. The Andean subduction
zone, at this latitude, is composed of three morphostructural units (Fig. 3.29), from west
to east: Coastal Cordillera (CC), Central Valley (CV) and Principal Cordillera (PC). A
variation in relief was considered in the synthetic model to represent the small elevation
in the region altitude.
Figure 3.29: Model 3 - Simplified model based on Reyes-Wagner et al. (2017) of a regional line
inSouth Volcanic Zone (SVZ) of the Chilean Andes.
Below the Cordillera Principal we consider three bodies with ρ1 = 10Ω m that char-
acterize the region of volcanic activity, with partial melting. Two structures were also
characterized as a product of the region’s tectonic activity, associated with sedimentary
rocks and weathered basement rock, with resistivities of ρ2 = 100Ω m, below the Central
Valley and the western Main Range. The value for the medium is ρ3 = 1000Ω m. The
inversion grid was defined with 3500 cells (100 × 35) 1.5 km in width, with thickness of
750 m at the top and increasing with depth at a rate of 5%, with a 60 km limit.
38
0 50 100 150 200
Iterations
10
-4
10
-2
10
0
10
2
10
4
|φ
k
+
1
−
φ
k
|
Convergence
LM − BFGS
Gauss−Newton
Figure 3.30: The absolute values of successive φ differences in the inversion of model 3.
The LM-BFGS inversion converged to a solution with 186 iterations in more than
347 minutes. For the same data set, the inversion with Gauss-Newton reached the same
solution in 22 iterations and over 90 minutes. Figure 3.30 shows the absolute values of
the variation of the function φ in successive iterations.
The LM-BFGS achieved similar results as the Gauss-Newton method, as can be seen
in figures 3.31 and 3.32. The same data are shown in figure 3.33 with contour lines.
The Gauss-Newton code needs to store the Hessian matrix (size 3500 × 2720) as well
as the sensitivity matrix (1250 × 3500) for a total of 16 625 000 values. Together with
the requirements of the forward problem, this results in a total of 312.7 MB of memory
used.
In an iteration of the LM-BFGS program, only the 5 {sk,yk} pairs of vectors from
the previous 5 iterations are stored, so that each iteration stores 5 × 2 × 3500, or 35000
values, only 0.2% of the amount stored in the Gauss-Newton method. The total memory
used by the LM-BFGS program reaches a maximum of 81.2 MB, which is about 25.9 %
of the total memory requirement of the Gauss-Newton program.
Table 3.3 summarizes these differences between the two methodologies, in terms of
stored units, number of iterations and processing time.
The data fit reached by both methods is shown figures 3.34 and 3.35 for the observa-
tions from the first and last stations on the survey line.
In this comparison, the differences in memory usage and processing time are all in the
inverse problem. The total memory used in the forward problem is the same for both
methods.
Apparent resistivity and phase pseudosections estimated by Gauss-Newton and LM-
BFGS inversions are shown in the figures 3.37, 3.38, 3.40 and 3.41. The observation
pseudosections of apparent resistivity and phase are the figures 3.36 and 3.39.
39
-6 -4 -2 0 2 4 6
×10
4
0
1
2
3
4
5
×10
4 (Gauss−Newton) Inverted data from ρa and φ − TE mode − MT
0.5
1
1.5
2
2.5
3
log10ρ
Figure 3.31: Result of the model 3 inversion with Gauss-Newton for a grid with 3500 cells. The
true target bodies positions are shown by the white lines.
-6 -4 -2 0 2 4 6
×10
4
0
1
2
3
4
5
×10
4 (LM −BFGS) Inverted data from ρa and φ − TE mode − MT
0.5
1
1.5
2
2.5
3
log10ρ
Figure 3.32: Results of the model 3 inversion with the LM-BFGS quasi-Newton method for a
3500 cells grid. The true target bodies positions are shown by the white lines.
Tabela 3.3: The main differences between the inversions of the model 3
Method Stored units Number of itera-tions Processing time
Gauss-Newton 16625000 22 90.1 min
LM-BFGS quasi-
Newton 35000 186 347.2 min
40
-6 -4 -2 0 2 4 6
×10
4
0
1
2
3
4
5
6
×10
4 (Gauss−Newton) Inverted data from ρa and φ − TE MT
0.5
1
1.5
2
2.5
3
log10ρ
(a) Model 3 Gauss-Newton inversion
-6 -4 -2 0 2 4 6
×10
4
0
1
2
3
4
5
6
×10
4 (LM −BFGS) Inverted data from ρa and φ − TE MT
0.5
1
1.5
2
2.5
3
log10ρ
(b) Model 3 LM-BFGS inversion
Figure 3.33: Model 3 inversions in contour lines
41
700
800
900
1000
1100
ρ
a
(Ω
·
m
)
Observation
LM − BFGS
Observation
Gauss−Newton
0.001 0.01 0.1 1 10 100 
T (s)
40
45
50
φ
(◦
)
0.01 0.1 1 10 100 1000
T (s)
Figure 3.34: Observed and estimated data for the first station in model 3.
500
600
700
800
900
1000
1100
ρ
a
(Ω
·
m
)
Observation
LM − BFGS
Observation
Gauss−Newton
0.001 0.01 0.1 1 10 100 
T (s)
35
40
45
50
55
φ
(◦
)
0.01 0.1 1 10 100 1000
T (s)
Figure 3.35: Observed and estimated data for the last station in model 3.
42
(Observation) Apparent Resistivity
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-6 -4 -2 0 2 4 6
Distance (m) ×104
200
300
400
500
600
700
800
900
1000
1100
1200
ρa
Figure 3.36: Apparent resistivity pseudosection (Model 3 - Observation).
(Gauss−Newton) Apparent Resistivity
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-6 -4 -2 0 2 4 6
Distance (m) ×104
200
300
400
500
600
700
800
900
1000
1100
1200
ρa
Figure 3.37: Apparent resistivity pseudosection (Model 3 - Gauss-Newton).
(LM −BFGS) Apparent Resistivity
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-6 -4 -2 0 2 4 6
Distance (m) ×104
200
300
400
500
600
700
800
900
1000
1100
ρa
Figure 3.38: Apparent resistivity pseudosection (Model 3 - LM-BFGS).
43
(Observation) Phase
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-6 -4 -2 0 2 4 6
Distance (m) ×104
30
35
40
45
50
55
60
65
φ
Figure 3.39: Phase angle pseudosection (Model 3 - Observation).
(Gauss−Newton) Phase
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-6 -4 -2 0 2 4 6
Distance (m) ×104
30
35
40
45
50
55
60
65
φ
Figure 3.40: Phase angle pseudosection (Model 3 - Gauss-Newton).
(LM −BFGS) Phase
1000 
100 
10 
1 
0.1 
0.01 
0.001
P
e
r
io
d
(s
)
-6 -4 -2 0 2 4 6
Distance (m) ×104
30
35
40
45
50
55
60
65
φ
Figure 3.41: Phase angle pseudosection (Model 3 - LM-BFGS).
4 CONCLUSIONS
Because the LM-BFGS methodology works with an approximation of the inverse of
the Hessian matrix, it has a lower convergence rate and needs more iterations to achieve
the same solutions as the Gauss-Newton method. In this implementation, both methods
require the calculation of sensitivities, which is the most time consuming task in the
inversion. The LM-BFGS program only needs a few columns of the sensitivity matrix
at the same time. Still, the number of calculations ends up being the same because all
sensitivities are eventually required in each iteration.
Apart from the calculation of sensitivities, an iteration of the Gauss-Newton method
needs to perform the operations to solve equations with the Hessian matrix, whereas the
LM-BFGS approach only needs matrix-vector multiplications, therefore an iteration in
the former takes longer than one in the latter. However, the number of iterations required
by the LM-BFGS program is much higher, so that the total computation time is longer
than that of the Gauss-Newton method.
The real advantage of the LM-BFGS approach is obviously the reduction in the amount
of memory required to perform the inversion. Because there’s no need to store the entire
Hessian and Jacobian matrices, this gain becomes bigger as the number of parameters
and observations increase.
The next steps in this research are: (1) to search for methods to avoid the calculation
of sensitivities by calculating the gradients directly from the equations of the forward
problem; (2) to implement applications to 3D problems, where the number of parameters
and observations are much greater; and (3) to implement applications to the inversion of
data from electromagnetic methods other than magnetotellurics.
44
REFERENCES
Avdeeva, A., and D. Avdeev, 2006, A limited-memory quasi-newton inversion for 1d mag-
netotellurics: Geophysics, 71, G191–G196.
Broyden, C. G., 1970a, The Convergence of a Class of Double-rank Minimization Algo-
rithms 1. General Considerations: IMA Journal of Applied Mathematics, 6, 76–90.
Broyden, C. G. 1970b,The Convergence of a Class of Double-rank Minimization Algo-
rithms: 2. The New Algorithm: IMA Journal of Applied Mathematics, 6, 222–231.
Buckley, A., and A. Lenir, 1983, Qn-like variable storage conjugate gradients: Mathemat-
ical Programming, 27, 155–175.
Buckley, A. G., 1978, A combined conjugate-gradient quasi-newton minimization algo-
rithm: Mathematical Programming, 15, 200–210.
Byrd, R., P. Lu, J. Nocedal, and C. Zhu, 1995, A limited memory algorithm for bound
constrained optimization: SIAM Journal on Scientific Computing, 16, 1190–1208.
Farquharson, C. G., and D. W. Oldenburg, 1996, Approximate sensitivities for the elec-
tromagnetic inverse problem: Geophysical Journal International, 126, 235–252.
Fletcher, R., 1970, A new approach to variable metric algorithms: The Computer Journal,
13, 317–322.
Gill, P. E., and W. Murray, 1979, "conjugate-gradient methods for large-scale nonlin-
ear optimization.": Technical Report SOL 79-15, Department of Operations Research,
Stanford University.
Goldfarb, D., 1970, A family of variable-metric methods derived by variational means:
Mathematics of Computation, 24, 23–26.
Greenstadt, J., 1970, Variations on variable-metric methods: Mathematics of Computa-
tion, 24, 1.
Haber, E., 2005, Quasi-newton methods for large-scale electromagnetic inverse problems:
Inverse Problems, 21, 305–323.
Meunier, G., 2008, The finite element method for electromagnetic modeling.
Nazareth, L., 1979, A relationship between the bfgs and conjugate gradient algorithms
and its implications for new algorithms: 16.
Newman, G., and P. Boggs, 2004, Solution accelerators for large scale 3d electromagnetic
inverse problems: Inverse Problems, 20, S151.
Ni, Q., and Y. Yuan, 1997, A subspace limited memory quasi-newton algorithm for large-
scale nonlinear bound constrained optimization: Math. Comput., 66, 1509–1520.
Nocedal, J., 1980, Updating quasi newton matrices with limited storage: Mathematics of
Computation, 35, 951–958.
Nocedal, J., and S. J. Wright, 1999, Numerical optimization: Springer.
Perry, A., 1977, A class of conjugate gradient algorithms with a two-step variable metric
memory: Discussion Papers 269, Northwestern University, Center for Mathematical
Studies in Economics and Management Science.
Reyes-Wagner, V., D. Diaz, D. Cordell, and M. Unsworth, 2017, Regional electrical struc-
ture of the andean subduction zone in central chile (35°–36°s) using magnetotellurics:
Earth, Planets and Space, 69.
Shanno, D. F., 1970, Conditioning of quasi-newton methods for function minimization:
Mathematics of Computation, 24, 647–656.
Shanno, D. F. 1978a, Conjugate gradient methods with inexact searches: Mathematics of
Operations Research, 3, 244–256.
45
46
Shanno, D. F. 1978b, On the convergence of a new conjugate gradient algorithm: SIAM
Journal on Numerical Analysis, 15, 1247–1257.
Tikhonov, A. N., 1963, Regularization of incorrectly posed problems: Soviet Mathematics
Doklady, 4, 1624–162.
Wannamaker, P., J. Stodt, and L. Rijo, 1986, Two-dimensional topographic responses in
magnetotellurics using finite elements: Geophysics, 51.