Cet ouvrage fait partie de la bibliothèque YouScribe
Obtenez un accès à la bibliothèque pour le lire en ligne
En savoir plus

2-D Niblett-Bostick magnetotelluric inversion

De
16 pages
A simple and robust imaging technique for two-dimensional magnetotelluric interpretations has been developed following the well known Niblett-Bostick transformation for one-dimensional profiles. The algorithm processes series and parallel magnetotelluric impedances and their analytical influence functions using a regularized Hopfield artificial neural network. The adaptive, weighted average approximation preserves part of the nonlinearity of the original problem, yet no initial model in the usual sense is required for the recovery of the model
rather, the built-in relationship between model and data automatically and concurrently considers many half spaces whose electrical conductivities vary according to the data. The use of series and parallel impedances, a self-contained pair of invariants of the impedance tensor, avoids the need to decide on best angles of rotation for identifying TE and TM modes. Field data from a given profile can thus be fed directly into the algorithm without much processing. The solutions offered by the regularized Hopfield neural network correspond to spatial averages computed through rectangular windows that can be chosen at will. Applications of the algorithm to simple synthetic models and to the standard COPROD2 data set illustrate the performance of the approximation.
Voir plus Voir moins

Geologica Acta, Vol.8, Nº 1, March 2010, 15-31
DOI: 10.1344/105.000001513
Available online at www.geologica-acta.com
2-D Niblett-Bostick magnetotelluric inversion
*J. RoDRíguez F.J. espa Rza and e. gómez-TReviño
CiCese/Ciencias de la Tierra
Km 107 carr. Tijuana-Ensenada, Ensenada, B.C., 22860, México
* Corresponding author E-mail: fesparz@cicese.mx
ABSTRACT
A simple and robust imaging technique for two-dimensional magnetotelluric interpretations has been developed
following the well known Niblett-Bostick transformation for one-dimensional profles.
The algorithm processes series and parallel magnetotelluric impedances and their analytical infuence functions
using a regularized Hopfeld artifcial neural network. The adaptive, weighted average approximation preserves
part of the nonlinearity of the original problem, yet no initial model in the usual sense is required for the recovery
of the model; rather, the built-in relationship between model and data automatically and concurrently considers
many half spaces whose electrical conductivities vary according to the data. The use of series and parallel imped-
ances, a self-contained pair of invariants of the impedance tensor, avoids the need to decide on best angles of rota-
tion for identifying TE and TM modes. Field data from a given profle can thus be fed directly into the algorithm
without much processing. The solutions offered by the regularized Hopfeld neural network correspond to spatial
averages computed through rectangular windows that can be chosen at will. Applications of the algorithm to
simple synthetic models and to the standard COPROD2 data set illustrate the performance of the approximation.
KEYWORDS Niblett-Bostick inversion. Hopfeld Neural Network. magnetotelluric.
INTRODUCTION imposed to construct realistic solutions that ft the data to a
given degree. Most commonly, the norm of the solution or
Electromagnetic inverse problems in geophysics are of its frst or second derivative are minimized together with
nonlinear. Even relatively simple cases, like the one- the misft to the data. This technique, frst developed for
dimensional (1-D) magnetotelluric (MT) problem, require 1-D problems, can readily be applied in higher dimensions.
special treatment to fully handle nonlinearities (e.g., Applications to the 2-D MT inverse problem include those
Bailey, 1970; Weidelt, 1972; Parker, 1983). Methods based of Rodi (1989), de Groot-Hedlin and Constable (1990),
on linearization can be applied iteratively to handle the Smith and Booker (1988) and Rodi and Mackie (2001).
nonlinearity of the problem (e.g., Oldenburg, 1979; Smith Minimizing roughness avoids the appearance of sharp
and Booker, 1988). In practice, as well as not being linear, features in the solution models that are not strictly required
electromagnetic inverse problems are ill-posed and severely by the data. That is, the resulting models are as smooth and
underconstrained. Sensible external constraints are usually even as the data permit.
15J. RODRÍGUEZ et al. 2-D Niblett-Bostick
It is also possible to obtain useful and somewhat 2-D: a) it captures the basic physics in a simple integral form,
more general information by slightly shifting the focus b) it relates the data directly to an arbitrary conductivity
of attention; instead of looking for a single model that fts distribution and, c) it relies on the relatively simple theory
the data, one can ask for general properties of all possible of a homogeneous half-space, but it handles heterogeneous
models that ft the data. The method of Backus and Gilbert media by adapting its conductivity according to the data,
(e.g., Backus and Gilbert, 1968, 1970) allows for the without the need for an initial or reference model. These
computation of spatial averages by means of averaging features are exploited by the 1-D extensions mentioned
functions constructed as linear combinations of the Fréchet above in relation to controlled source methods. Extensions
derivatives of the data. The averaging functions are made to higher dimensions for special types of electromagnetic
to resemble box-car functions for the averages to have the measurements also proft from them (e.g. Pérez-Flores and
usual intuitive meaning. The results are average models for Gómez-Treviño, 1997; Pérez-Flores et al., 2001; Brunner
given window sizes. The models are not intended to produce et al., 2003; Friedel, 2003). In this paper, we explore the
responses that ft the data. In fact, they seldom do better possibilities of the same approach for the two-dimensional
in this respect than models designed specifcally for ftting (2-D) MT problem.
purposes. It is perhaps for this reason that average models
are not very popular among interpreters of feld data who
seem to prefer the assurance of a direct ft to the data. THE APPROXIMATION
Summarizing the above two paragraphs, we have on one In the magnetotelluric method, surface measurements
side an optimization process of ftting data with external of natural time-varying electric and magnetic felds are
constraints, and on the other side, an optimization process readily converted to four complex impedance values per
of ftting Fréchet derivatives to box-car functions. Although given angular frequency ω. In turn, these values are usually
the two processes are complementary, particularly for normalized to obtain apparent resistivities or, equivalently,
nonlinear problems, the latter is seldom applied, perhaps apparent conductivities, by referring the actual impedances
because it requires intensive computations, but most likely to those of a homogeneous half-space (Cagniard, 1953).
in view of the reason stated in the previous paragraph. In Here we use apparent conductivity s as derived from the a
this paper, we present an approach that combines features magnitude of a complex impedance Z which, for the mo-
from both methods. On one side, we keep the reassuring ment, represents any of the four elements of the impedance
feature of constructing models whose responses optimize tensor. The formula for apparent conductivity is simply
the ft to the data and, on the other, we maintain the concept given as
σ aof spatial averages. Spatial averages not only appeal to −2intuition, but when plotted against depth they resemble a σ ( x,ω) = ωµ |Z(x,ω) | , (1)(1) a 0σ aprofle of the property itself, fltered by the corresponding µ −20window. As shown in Gómez-Treviño (1996), a solution in where µ stands for the magnetic permeability of free-space σ (x,ω) = ωµ |Z(x,ω) | , (1)
0 a 0
σ (x,w) terms of averages represents a robust alternative to the well and x represents horizontal distance in an x-z coordinate a µ 0σ known 1-D Niblett-Bostick transformation. system whose z axis represents depth. s (x, w) represents a a2π
T = −2σ (x,w) the data at a given distance in a 2-D model with a fat a σ (x,ω) = ωµ |Z(x,ω) | , (1) a 0ωBetween linearization and nonlinear methods, there are topography and for a given angular ω. The data are usually 2πµ σ 0 a T = approximations simple enough to be handled analytically presented as individual sounding curves as a function of
ω 2 pσ (x,w) but that still keep some of the nonlinear features of the σ (x,w) period T = for different distances x, or in a pseudo-a a vσ aoriginal problem. Such is the case of the Niblett-Bostick section format contouring values of s as plotted over 2π a σ (,x z) T = approximate integral equation (e.g., Niblett and Sayn- (x,w) x-T coordinates. s (x, w) represents what is available; σ aω a Wittgenstein, 1960; Bostick, 1977; Jones, 1983) which what is required is s(x, z), the subsurface conductivity σ σ (,x z) a σ (x,T ) has inspired a number of generalizations. F a or instance: distribution.σ (x,w) a iteration of the corresponding exact equation (Esparza and σ (,x z)
σ (,x z) σ (x,T ) Gómez-Treviño, 1996) and, still within the approximation, A useful relationship between s (x, T) and s(x, z) for a a1 σ (x,T ) = F(x,xz', ',σσ,T ) (xz', ')dx'dz ', (2) a solution of the non-uniqueness problem by direct the problem at hand is (Gómez-Treviño, 1987a):σ (,x z) a ∫1−mσ (x,T ) computations of spatial averages (Gómez-Treviño, 1996). a 1d logσOther generalizations include the e σxtension of the basic (,x z) σ ( x,T ) = Fa(x,xz', ',σσ,T ) (xz', ')dx'dz ',(2)a m = ∫ , (3) 1−midea to 1-D inversion of controlled source electromagnetic d logT1
σ (x,T ) = F(x,xz', ',σσ,T ) (xz', ')dx'dz ', (2) d logσa ∫data (e.g., Boerner and Holladay, 1990; Smith et al., 1994; where: aF(xx, ',z ',σ ,T ) 1−m m = , (3)
Christensen, 1997). d logTd logσ aσ (x,T ) m = , (3) a (3)F(xx, ',z ',σ ,T ) d logT
The 1-D Niblett-Bostick approximation has the σ (,x z)
F(xx, ',z ',σ ,T ) σ (x,T ) afollowing basic features that are worthwhile exploiting in and F(x,x’,z’, s,T) represents the Fréchet derivative of σ (,x z)
σ (x,T ) σ (,x z) a
σ (x,T ) aσ (,x z) σ (,x z)
F(z ',σ ,)T = (1−mF) (z ',σ ,)T . (4) σ (,x z) haσ (x,T ) aFz( ',σ ,T ) σ (x,T ) 16Geologica Acta, 8(1), 15-30 (2010) a
F(z ',σ ,)T = (1−mF) (z ',σ ,)T . (4) DOI: 10.1344/105.000001513 haFz ( ',σ ,T ) F(z ',σ ,)T = (1−mF) (z ',σ ,)T . (4) ha haFz( ',σ ,T )
Fz( ',σ ,T ) F h
Fz ( ',σ ,T ) haFz ( ',σ ,T ) haFz ( ',σ ,T ) ha
F F hh ∞ 2 2'zz 2'−2z '/δa Fz (σσ',()σTe,=T ) cos( )+ sin( ) (z ')dz ', (5) Fz ( ',σ ,T ) a ha ha ∫ 0δ δδa aa
∞∞2 22'zz 2'2'zz 2'−2z '/δ −2z '/δa a σσ()Te= cos( )+ sin( ) (z ')dz ', (5) σσ()Te= cos( )+ sin( ) (z ')dz ', (5) a ∫ a δ = 503 T /δ 0 ∫δ δδ0a aa δ aa δδa aa
0.707δ δ = 503 T /δ aa aδ = 503 T /δ a a
−1
0.707 (0δ.707 σ ) a a0.707δ a−1(0.707σ ) . σ =Aσ (6) a −1 a(0.707σ ) a . σ =Aσ (6) σ aa . σ =Aσ (6) aσ a 2z 2z(−2z / δ ) (−2z / δ )j j+1j ai j +1 aiσ a =e cos( ) −e cos( ), (7) a 2z 2zij(−2z / δ ) (−2z / δ )j j+1j ai j +1 aiδ δ a =e cos( ) −e cos( ), (7) ai aiij 2z 2zδ (−2z / δ ) j δ (−2z / δ ) j+1j ai j +1 aiai ai a =e cos( ) −e cos( ), (7) z ijj
z δ δj ai ai
δ aiδ z ai j
δ ai


1 1
1
σ a
−2 σ (x,ω) = ωµ |Z(x,ω) | , (1) a 0
µ 0
σ (x,w) a

T =
ω σ a
σ a −2σ a σ (x,ω) = ωµ |Z(x,ω) | , (1) a 0σ (x,w) −2a σ (x,ω) = ωµ |Z(x,ω) | , (1) J. RODRÍGUEZ et al. a 0 2-D Niblett-Bostickµ 0σ (,x z)
µ 0σ (x,w) a
σ (x,w) 2π aσ (x,T ) a T =
2πω s (x, T) with respect to s(x, z). The integration is defned By analogy, making the same type of assumptions and σ (,x z) a T =
over the entire lower half-space. The recovery of s(x, z) from approximations in 2-D, equation (2) can be written asσ ωa 1
σ (x,T ) = F(x,xz', ',σσ,T ) (xz', ')dx'dz ', (2) s (x, T) is clearly a nonlinear problem since F depends on aa σ∫ σ (x,w) aa 1−ms. Otherwise, the integral equation could be readily solved (8) σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) a ha ∫σ (x,w) σ (,x z) d logσausing any of the available methods of linear analysis. a It m = , (3) σ F a hσ (,x z) is still possible to apply linear methods sequentially, as Analytical expressions for F are derived in Appendix A hd logT
−2 σ in traditional linearization, simply by updating a starting σ (x,ω) = ωµ |Z(x,ω) | , (1) for the traditional TE and TM modes, respectively. In turn, aσ (x,T ) a 0F(xx, ',z ', ,T ) a
model on the right-hand side of the equation. Esparza and they are used in Appendix C to derive the corresponding σ (x,T ) µσ (,x z) aσ (x,T ) 0a Gómez-Treviño (1996), working with the 1-D problem, expressions for series and parallel apparent conductivities.
σ σ (,x z) j1σ (x,w) σ (,x z) showed that reasonably good results can be obtained in a a σ (x,T ) = F(x,xz', ',σσ,T ) (xz', ')dx'dz ', (2) a ∫single iteration using an adaptive approximation. In 1-D, 1 T In 2-D, to construct equation (6) the half-space is 1−mσ (,x z) 2π ij σ (x,T ) = F(x,xz', ',σσ,T ) (xz', ')dx'dz ', (2) aT = ∫the approximation is divided into a large number of rectangular elements. The d logσ 1−ma σ , j = 1,...,N, σ (x,T ) ω ja m = , (3) integration over the elements can be performed analytically d logσd logT aσ a m = , (3) F(z ',σ ,)T = (1−mF) (z ',σ ,)T . ˆ (4) (4) as described in Appendix B. This is particularly useful for σ ,k =1,...,M ha ak
F(xx, ',z ',σ ,T ) d logTσ (x,w) handling the singularities at the points of measurement, Fz( ',σ ,T ) Na
F(xx, ',z ',σ ,T ) F(z’, s,T) on the left hand-side represents the true Fréchet and also for the fnal rectangles on the sides and bottom of σ (x,T ) ˆ σ = a σ ,k = 1,...,M. (9) σ (,xa z) ∑ak kj jFz ( ',σ ,T ) ha derivative for an arbitrary conductivity distribution, and the model. It is worth remarking that each of the elements j =1σ (x,T ) σ (,x z) a F h F (z’, s ,T) on the right stands for the much simpler of s is a weighted average of all the unknown conductivity h a aσ , σ (,x z) akσ (,x z) σ (,T ) Fréchet derivative of a homogeneous half-space, whose values and, that on virtue of equation (2), the elements Fz ( ',σ ,Ta) ha
σ (,x z) σ ˆ σ (x,T ) conductivity is known and equal to the measured apparent ak of matrix A are dimensionless. Furthermore, the sum of σ (,x z) a ∞ 2 2'zz 2'−2z '/δa M M Nvity. F is simply an attenuated cosine function the elements of any row of A is identically unity, which σσ()Te= cohs( )+ sin( ) (z ')dz ', (5) σ (x,T ) 1 1a  ∫ a F(1z ',σ ,)T = (1−mF) (z ',σ ,)T . (4) 2 20 haδ δδ C = ( σ − σ ˆ ) = [ σ − a σ ] . (10) athat gradually vσ (x,T ) = anishes with depth. aaF(x,xz', ',σσ,TThe f) (xzactor (1-', ')dx'dmz) drops ', (2) is a very useful property for checking the accuracy of the ∑ ak ak ∑ ak ∑ kj ja ∫ F(z ',σ ,)T = (1−mF) (z ',σ ,)T . (4) Fz( ',σ ,T ) 2 21−m ha k =1 k =1 j =1out of the approximation when substituting expression (4) computations involved. Notice that although equation (6) δ = 503 T /δ a a Fz( ',σ ,T ) d logσFz ( ',σ ,T ) in equation (2). Using the Fréchet deri a vative F (z’, s ,T) is a system of linear equations, the model it represents h aha m = , (3)
N N M N M M0.707δ for a homogeneous half-space (Gómez-Treviño, 1987b), is actually nonlinear, for A depends on the unknown Fz ( ',σ ,T ) a F d logT 1 1hah 2C = − [ − a a ] σ σ − [ − a σ + a σ ] σ +the approximation in 1-D can be written as distribution σ through the different values of s .−1 ∑ ∑ ∑ ki kj i j ∑ ∑ a ki i ∑ ki ak iF(xx, ',z ',σ ,T ) F (0.707σ )Fz ( ',σ ,T ) h 2 2 a ha i =1 j¹i =1 k =1 i =1 k =1 k =1
σ (x,T ) Fz ( ',σ ,T ) ha ∞ .a 2σ =Aσ 2'zz 2' (6) −2z '/δa a σσ()Te= cos( )+ sin( ) (z ')dz ', (5) a ∫σ (,x z) 0 ∞  HOPFIELD ARTIFICIAL NEURAL NETWORKS2 2'zz 2'δ δδ M−2z '/δσ a aaa 1a ()Te= cos( )+ sin( ) (z ')dz ', (5) 2a ∫ (5)σ (,x z) 0 ( σ ) . (11) δ δδ ∑ aka aaδ = 503 T /δ 2z 2z 2The application of artifcial neural networks to the a (−2z / δa ) j (−2z / δ ) j+1 k =1j ai j +1 ai σ (xa,T =) e cos( ) −e cos( ), (7) a ij where . In the original Niblett-Bostick δ = 503 T /δ inversion of MT data has been explored in various a aδ δ0.707δ ai aia integral equation, F(z ',σ ,)Tthe upper limit of inte= (1−mF) (z ',σ ,)gration T . is 0.707 d (4) directions. One way is to use the multi-layer feed forward a N N Nha 1−1 0.707δ z a -1j (0.707σ ) and the kernel is simply (0.707 s ) . neural netwE = − ork architecture (Rummelhart et al., 1986) T σ σ − I σ , (12) aFz( ',σ ,T )a ∑ ∑ ij i j ∑ i i
−1 2 i =1 j¹i =1 i =1which uses a set of responses and models presented to the (0.707σ ) δ ai a . σ =Aσ (6) Fz ( ',σ ,T ) aha To solve equation (5) numerically, we divide the input and output defned neurons, respectively. During a . σ =Aσ (6) F aσ h half-space into a large number of layers with uniform learning phase, the netwM ork bacMk-propagatesM (through its a 1 2 conducti σvities. The result is that the integral equation can T = − neurons and interconnection weights) errors due to misft a a and I = − a σ + a σ . (13) Fz ( ',σ ,T ) ∑ ∑ ∑a 2z 2z ij ki kj i ki i ki akha (−2z / δ ) j (−2z / δ ) j+1j ai j +1 ai 2k =1 k =1 k =1 be written as a matrix equationa =e cos( ) −e cos( ), (7) of the model and the obtained neural model. Learning from ij ∞  2z 2z2 2'zz 2'−2z '/δ (−2z / δ ) j (−2z / δ ) j+1δ j ai δ j +1 aia ai ai one response-model ‘pattern’ is achieved by updating the σσ()Te= cos( )+ sin( ) (z ')dz ', (5) a =e cos( ) −e cos( ), (7) 1 a ∫ ij0δ δδ= δ δs A s. a aa (6) inter neuron connection weights according to a gradient Nz a ai ai  j 1 1(t +1) (t)  descent minimization criteria. σ = + sgn TThe process is then applied σ +I . (14) z i ∑ ij j iδ = 503 T /δ  jδ a a 2 2ai j ≠i =1The vector s contains the data for the different to the complete response-model data set, thereby achieving a  
δ ai0.707δ periods. The vector s represents the unknown conductivity a learning epoch. a
distribution in its discrete form, and the matrix A contains (t +1)−1 σ (0.707σ ) ia the weights of the conductivity elements for all the Once the network is trained, it recovers a model in
(t) . σ =Aσ (6) available data. The elements of the matrix can be e σvaluated almost no time when provided with a sounding curve. The a j
1 analytically as distinctive feature of this learning approach is that there σ T a ij 1 is very little physics fed into the algorithms. In fact, as far
T2z 2z ()AA as the algorithms are concerned, the models and responses (−2z / δ ) j (−2z / δ ) j+1j ai j +1 ai ij a =e cos( ) −e cos( ), (7) ij used in the training sessions may or may not be related δ δ ai ai (7) through any physical link, the learning process is simply
z j the same. Hidalgo and Gómez-Treviño (1996) explored

where z is the top depth of the j-th layer, and d is the skin this approach for the 1-D problem with reasonably good j aiδ ai depth of the i-th measurement. results. However, extending the method to 2-D would be

2
1 17Geologica Acta, 8(1), 15-30 (2010)
DOI: 10.1344/105.000001513 σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) a ha σ (x,T ) = ∫ F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) a ha∫
F Fh h
σ σ a σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) aa ha∫
F h σ σ j jσ a
T T σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) ij a ha∫ ij
σ , j = 1,...,N, σ j, j = 1,...,N, σ F h jj
ˆσ ,k =1,...,M σ ˆT σ ak,k =1,...,M a akij
NN σ , j = 1,...,N, j ˆ σ = a σ ,k = 1,...,M. (9) ˆ ∑ σ ak = a kj σ j,k = 1,...,M. (9) ak ∑ kj jσ j j =1σ ˆ ,k =1,...,M j =1ak
T σ , Nij σ ak, ak
ˆ σ = a σ ,k = 1,...,M. (9) ak ∑ kj jσ , j = 1,..., σ ˆN, j σ ˆ ak akj =1
M M NM M Nσ ˆ ,k =1,...,M 1 12 2ak 1 1σ , 2 2 C = ( σ − σ ˆ ) = [ σ − a σ ] . (10) ak C = ∑ ( σ ak − σ ˆ ak) = ∑ [ σ ak − ∑ a kj σ j] . (10) ∑ ak ak ∑ ak ∑ kj jN2 2k =1 k =1 j =12 2ˆ k =1 k =1 j =1σ ak σ ˆ = a σ ,k = 1,...,M. (9) ak ∑ kj j M M N j =11 12 2N N M N M MN N M N M M C = ( σ − σ ˆ ) = [ σ − a σ ] . (10) 1 1∑ ak ak ∑ ak ∑ kj j 21 1σ , 2ak C = − [ − a a ] σ σ − [ − a σ + a σ ] σ +2 2k =1 C = −k =1 ∑ ∑ j[ = −1 ∑ a kia kj] σ σi j − ∑ [ − ∑ a ki σ i + ∑ a ki σ ak] σ i +∑ ∑ ∑ ki kj i j ∑ ∑ ki i ∑ ki ak i2 2 i =1 j¹i =1 k =1 i =1 k =1 k =12 2 ˆσ i =1 j¹i =1 k =1 i =1 k =1 k =1 ak
M M NN N M N M M1 1 1 12 2 2ˆ C = ( σ − σ ) = [ σ − a σ ] . (10) C = − [ − a a ] σ σ − [ − a σ + aM σ ] σ +∑ ∑ ∑ ∑ ak ∑ ak ∑ ∑ ak ∑∑ kj jki kj i j ki i Mki ak i1 22 2 12 2 k =1 k =1 j =1 2i =1 j¹i =1 k =1 i =1 k =1 k =1 ( σ ) . (11) ∑ ( σ ak) . ak2 k =1 2 k =1
J. RODRÍGUEZ et al. 2-D Niblett-Bostick
N N M N M M 1 1M 21C = − [ − a a ] σ σ − [ − a σ + a σ ] σ +2 N N N∑ ∑ ∑ ki kj i j ∑ ∑ ki i ∑ ki ak iN N N1 ( σ ) . (11) 1 ak2 2 i =1 j¹i =1 k =1E = − i =1T σ σ k −=1 I σ ,k =1 (12) 2 E = − ∑ ∑ Tij σ σi j − ∑ I σi i, k =1 ij i j i i22 i =1 j¹i =1 i =1i =1 j¹i =1 i =1 extremely cumbersome, for it is necessary to cover many where N N N M1 1possibilities during each epoch training session. 2M M MM M M E = − T σ σ − I σ , 1(12) ( σ ) . (11) ∑ ∑ ij i j ∑ i i ∑ 1 2ak 2 T = − a aand I = − a σ + a σ . (13) 2 ∑ 2 ∑ ∑i =1 j¹i = 1Tij = − ai =ki1a kj I i = − a ki σ i + a ki σ ak.k =1ij ki kj i ki i ki ak2Another approach to inversion based on neural networks k =1 k =1 k =12k =1 k =1 k =1
bypasses the learning sessions by directly providing the N N NM M M1 2algorithm with the relevant physics behind the particular N E = − T σ σ −  NI σ ,  (12) T = − a a and I a σ + a σ . (13) 1 1 (13)ij ∑ ki kj i ∑ ∑ki i ij ∑i j ki ∑ ak i i (t +1) (t)1 1(t +1) (t)application. Zhang and Paulson (1997) applied a simple   2 σ = + sgn T σ +I . (14) i =1 j¹i =1 i = 1 k =1 kσ i = k +=1 ∑ Tij σ j +I i .  i ij j i 2 2recursive regularization algorithm associated with a j ≠i =12 2  j ≠i =1  
Hopfeld artifcial neural network (HANN) in order to E is the Ising Hamiltonian (Ising, 1925) which Tank M M MN   11 1 2solve the 1-D MT inverse problem with excellent (t +results. 1) and Hopfeld (1986) sho (t) wed is a never increasing quantity (t +1) T = − a a and I = − a σ + a σ . (13) (t +1)   σ = sgn T σ + . (14) ij ∑ ki kj i ∑ ki i ∑ ki akσ i ∑ ij j iiσ  In the same paper the authors explored the application to i as the dynamics defned by the follo 2 wing equation evolve 2 2 kk =1 =1 k =1j ≠i =1 (t)(t)the 2-D full nonlinear inversion of MT data. Our interest from given starting values:σ σ j jin this algorithm was triggered by the fact that there are
N(t +1) T  1 1Tij σ (t +1) (t)no references in the literature attempting to continue their ij (14)i   σ = + sgn T σ +I . (14) ∑i ij j iT  T σ (x,T )research in this direction, e= F (x,x',z ',σσ,T ) (vxen though the results presented ',(t)z ')dx'dz '. (8) ()AA 2 2a ha j ≠i =1∫ ()AA ij  σ ijjare very promising. In this methodology there are no
F (t+l) (t)h T dedicated input-output defned units (as with the feed- s is the updated conductivity value of s which is either ij j j σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '(.t+1)a ha ∫ σ (x,T ) = F (x,x',z ',σσ,T ) (xσ ',z ')dx 'dz ' . (8) σ forward netw a ork architecture). Instead, ehaT very neuron i s a given starting value or the result of a previous iteration. ∫ i a ()AA ij F (t) h is interconnected with every other s neuron through a T The relevance of dynamic is that as iterations proceed j ij σ ( x,T ) =F F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) σ h ja ha∫ σ weight link as detailed below. Every neuron has a twofold E decreases at each step. Since E corresponds to C, the aσ σ j T F σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '. a input-output purpose. We assume M apparent conductiij vity dynamics ensure that the sum of the squared differences h a ha∫
TT measurements and N blocks whose conductivities are between data and model responses decrease at each step, σ σ ()AA ijF a j ijh unknown. According to equation (6), for gi ven conductivity thus leading to the model whose responses best ft the data σ σ , j = 1,...,N, j T σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) σ j 2 a ha ij∫ a values s , j = 1,...,N, the responses sˆ , k = 1,...,M, of the in a least square sense. It can be noted that with the HANN j ak, σ (x,T ) = F (x,x',z ',σσ,T ) (x',z ')dx'dz '. (8) a ha∫T σ σ , j = 1,...,N, σ ˆ ,k =1,...,M model can be expressed as the matrix elements T are known and provided by the j ij ijjF akh F Th physics of the problem, since they are the elements (A A) , ijNT σ , j = 1,...,N, σ σ ˆ ,k =1,...,M jσ j akija (9) which in turn come from the integral equation.ˆ 2 σ σ = a σ ,k = 1,...,M. (9) ak ∑ kj j NaT σ ˆ ,k =1,...,M σ , j = 1,...,N, j ak j =1ij ˆ σ = a σ ,k = 1,...,M. (9) ak ∑ kj j
N Strictly speaking, equation (14) applies only to binary j =1σ ˆσ , j σ = 1,,...,k =N 1,...,, σM, j j ak ak σ σ ˆ = a σ ,k = 1,...,M. (9) j The misft between s (the data) and sˆ (the model variables that can take values of zero and one, as in the HANN ak ∑ ak kj j ak,σ , ak NT ˆσ j =1σ ˆ ,k =1,...,M ij ak responses) can be written as architecture (Hopfeld, 1982). The corresponding equations for2 ak T σ ˆ = a σ ,k = 1,...,M. (9) ij ak ∑ kj jM M NN an assemblage of such variables, to make up for arbitrary values σ , σ , j = 1,...,N, 1ak 1j =1j 2 2M M Nσ , j = 1,...,N, C = ( σ 1− σ ˆ ) = [ σ 1 − a σ ] . (10) j σ ˆ = a σ ,k = 1,...,M. (9) 2 2 of conductivity, are very similar and are given in Appendix D.∑ ak ak ∑ ak ∑ kj jak ∑ kj j C = ( σ − σ ˆ ) = [ σ − a σ ] . (10) ˆσ σ , 2 ∑ ak 2 ak ∑ ak ∑ kj jσ ˆ ,k =1,...,M (10)akk =1 k =1 j =1ak j =1ak 2 2ˆ k =1 k =1 j =1σ ,k =1,...,M ak M M NN ˆ 1 1 The application of equation (14) is straightforward , σ σ 2 2ak ak N ˆ C = ( σ − σ ) = [ σ − a σ ] . (10) σ ˆ = a σ ,k = 1,...,M. (9) N N M ∑ ak akN ∑M ak ∑M kj jak ∑ kj j Squaring as indicated and rearranging terms, this when M>N, for the process reduces to the standard least N N M N M MM 1 M N 1ˆ σ =2 a σ ,k = 1,...,M2. 2 (9) 1 1ˆ 1 1 k =1 k =1 2j =1σ ak ∑ kj jj =1 2 2ak C = − [ − a a ] σ σ − [ − a σ + a σ ] σ +C = − [ − a a ] σ σ − [ − a σ + a σ ] σ +∑ ∑ expression corresponds to∑ ∑ ∑ ∑ squares problem. This corresponds to the over constrained ˆ ki kj i j ki i ki ak i C = ( σ − σ ) = ∑ ∑[ σ ∑− ki akj σi ]j . (10)∑ ∑ ki i ∑ ki ak ij =1∑ ∑ ∑ak ak ak kj j2 2 2 2 M M N i =1 j¹i =1 ik= =11 j¹i =1 k =1 i =1 i =k1 =1 k =1 k =1 k =1σ , 2 2 case when we have more data than unknown parameters. 1 k =1 1 k =1 j =1ak 2 2σ , ˆ N N M N M M C = ( σ − σ ) = [ σ − a σ ] . (10) ak∑ ak ak ∑ 1ak ∑ kj j 1 In general, however, M<N since there usually are fewer 2 ˆσ 2 2C = − [ − a a ] σ σ − [ − a σ + a σ ] σ +k =1 k1 j =1ak ∑ ∑ ∑ ki kj i j ∑ ∑ ki i ∑ ki ak iσ ˆ M data points than unknowns. We found, as did Zhang and N N M N M Mak2 1 2 2M M N 1 i =1 j¹i =1 1k =1 i =1 k =1 k =1 2 2 ( σ ) . (11) 1 1 Paulson (1997), that the algorithm still converges even for 2 C = − 2 [ − aMa ] σ σ − [ − M a σ ∑ + N aka σ ] σ + ( σ ) . (11) ∑ ∑ ∑ ki kj i j ∑ ∑ ki i ∑ ki ak iˆ 1 1 ak C = ( σ − σ ) = [ σ − a σ ] . (10) 22 2N N M N M Mk =1∑ ak ak ∑ ak ∑ kj j2 2 1 i =1 j¹i =1 Ck =1( σ − σ ˆ i =1) = k =1 [ σ − k =1a σ ] . (10) M<N, but the solution is hardly useful for any practical 22 2 ∑ ak ak ∑ ak ∑ kj jk =1 k =1 j =1C = − [ − a a ] σ σ − [ − a σ + a σ ] σ + ∑ ∑ ∑ ki kj2 i j ∑ ∑2 ki i ∑ ki ak iMk = 1 k =1 j =1 purpose because it is highly oscillatory and unrealistic. It is 12 2N N N 2 i =1 j¹i =1 k =1 i =1 k =1 k =11 ( σ ) . (11) N N N then necessary to make some extra assumptions about the ∑ ak E = − T σ σ − I σ , (12) 1N N M N M M M ∑ ∑ ij i j ∑ i i21 1 1 k =12 EN = −N M 2T σ σ 2 − NI σ , M M(12) solution. This issue is addressed in the next section.i =1 j¹i =1 i =1∑ ∑ ij i j ∑ i iC = − [ − a a ] σ σ − [ − a σ + a1 σ ] σ + ( σ ) . 1 (11) 2 (11)∑ ∑ ∑ ki kj i j ∑ ∑ ki i ∑ ki ak i ak2 i =1 j¹i =1 i =1C = − [M − a a ] σ σ − [ − a σ + a σ ] σ +2 2 ∑ ∑ ∑2 ki kj i j ∑ ∑ ki i ∑ ki ak i1i =1 j¹i =1 k =1 i =1 k =1 k =1 k =1 22 N N 2 N i =1 j¹i =1 k =1 i =1 k =1 k =1 M ( σ ) . M (11) M ∑ 1ak 1 2 2 E = − T σ σ − I σ , (12) T = − a a and I = − a σ + a σ . (13) k =1M It can be noted that this is an unnecessarily long, and ∑ M ∑ ij i M j ∑ i i SPATIAL AVERAGESij ∑ ki kj i ∑ ki i ki akN N N 12 221 k =1 i =1 j¹i =1 k =1 i =k1 =1M T = − a acertainly v and ery cumbersome w I = − aay of representing σ + a σ C. (13) , the sum 1 ij ∑ ki kj i ∑ ki i ∑ ki ak 2 E = − T σ σ − I σ , M (12) ∑ ∑ ij i j ∑ i i 2 1 ( σ ) . N N(11) Nk =1 k =1 k =21∑ of squares. However, this particular form helps to recognize True averagesak 21 i =1 j¹i =1 i =1 ( σ ) . (11) 2 ∑N akk =1 M M M E = − T σ σ − I σ ,  (12) 1 1 its correspondence to∑ ∑ ij i j ∑(t +1) i i 2 (t) 1k =1 2   σ = + sgn T σ +I . (14) 2 T = − i a a and ∑ I ij = −j i a σ + a σ . (13) i =1 j¹i =1 i =1 ij ∑ ki kj  i ∑  ki i ∑ ki akN As reviewed in the introductory section above, the 2 2M M  M  j ≠i =11 1  2 (t +1) (t)N N N k =1 1 k =1 k =12   σ = + sgn T σ +I . (14) methods of Backus and Gilbert allow for the computation of 1 T = − a a and i I = − N Na ∑σ +ij j a N σi . (13)  ij ∑ ki kj i ∑ ki i ∑ ki ak1 E = − T σ σ − I σ , (12) M 2M 2 M∑ ∑ ij i j ∑ i i j ≠i =1 2   spatial averages. This is done by means of averaging functions k =1 (12)1 k =1 k =1 (t +1) E = − T σ σ − I σ , (12) 22 ∑ ∑ ij i j ∑ i ii =1 j¹i =1 i =1 σ T = − a a and I = − a σ + a σ . (13) Niij ∑ ki kj i ∑ 2 ki i ∑ ki ak   constructed as linear combinations of the Fréchet derivatives 1 1 i =1 j¹i =1 i =1 (t +1) (t)2(t)  k =1 k =1 k =1 σ + sgn T σ +I . (14) ∑(t +1) σ i ij j ij   Nσ 2 2M M M  i 1 1 j ≠i =1  (t +1) (t)1 2 M   M M σ = + sgn T σ +I . (14) T i ij ∑ ij j i T = − a a and I = − a σ +(t) a σ . (13) 1  2ij ∑ ki kj i ∑ ki i ∑ ki ak N σ 2 2  T1 =1 − a a and I = − a σ + a σ . (13) j j ≠i =12 (t +1) T  (t) ij ∑ ki kj i ∑ ki i ∑ ki akk =1 k =1 k =1 ()AA  (t +1) σ = + sgn T σ +I . (14) ij 2i ∑ ij j ik =1 k =1 k =1σ  T i 18 Geologica Acta, 8(1), 15-30 (2010)2 2 ij j ≠i =1  (t) DOI: 10.1344/105.000001513(t +1) TN σ σ   j1 1 i ()AA (t +1) (t) ij N   1 1 σ = + sgn T σ +I . (14) (t +1) (t)i (∑t +1) (ijt) j i     T σ = + sgn T σ +I . (14) σ σ 2 2 ij i ∑ ij j ii jj ≠i =1    2 2 j ≠i =1 (t) T ()AA σ T ijj ij (t +1) Tσ (t +1) T ()AA i ij ij σ i(t) T 2 σ ()AA (t)j ij σ j T ij T ij T ()AA 2 Tij ()AA ij

2

2

2
2
2 J. RODRÍGUEZ et al. 2-D Niblett-Bostick
n
(t) (t) σ = w σ . (15) j ∑ l ( j +l)
l = −n
of the data in such a way that the averaging functions its performance in 1-D, and compare its results with
resemble boxcar functions for the averages to have the usual existing methods of average estimations. In the following N n 1 1(t +1) (t)meaning. The resulting models are not intended to produce σlines, we test the performance of equation (16) as it applies = + sgn T [ w σ ] +I . (16) i  ∑ ij ∑ l ( j +l) i 
2 2 j ≠i =1 l = −n responses that ft the data, but to address the non-uniqueness to the solution of equation (6). The elements of the matrix
(t +1)character of the inverse problem. Within the limitations of are computed using equation (7) which represents the 1-D σ i
linearization, they can be called true averages. More rigorous version of the proposed approximation. The synthetic data (t)σ japproaches for computing averages have been developed for for the tests, shown in Fig. 1A, correspond to apparent
1the 1-D magnetotelluric problem (e.g. Weidelt, 1992), but their conductivity responses of the model shown in Fig. 1B.w = ,k = −n,...,n kextension to higher dimensions are far from trivial. 2n +1
The tests are intended to demonstrate that the proposed
N nSimulated averages approach for the computation of averages has the same 1 1(t +1) (t) (t)σ = + sgn{ [(1 −b)T σ + β T w σ ] +I }. (17) i ∑ ij j ∑ ij l ( j +l) iproperties as the simple formula for averages in Gómez-2 2 j ≠i =1 l = −n
The shortcut that we take consists of computing models Treviño (1996). The formula is
that are averages of themselves everywhere with a given
σ (T ) δ (T ) − σ (T ) δ (T )a 2 a 2 a 1 a 1neighborhood, and that at the same time their responses ft (18) < σ (z ,z ) >= , (18) 1 2 δ (T ) − δ (T )the data at an optimal level. Consider that a 2 a 1
< σ (z ,z ) > n 1 2
(t) (t) (15) where < s(z ,z )> represents the average of conductivity σ = w σ . (15) 1 2∑ z = 0.707 δ (T ),i =1,2 j l ( j +l) i a in between z = 0.707 d (T ),i = 1,2. The averages are assigned l = −n i a i(t) (t) σ = w σ . (15) z = z z j ∑ l ( j +l) to the mean geometrical depth z = z z . Any two data 1 2 √ 1 2
l = −nThe dynamics of the HANN in (14) transforms to points can be used in the formula regardless of how far N n 1 1 (t +1) (t) apart they are in the sounding curve. If the points are σ = + sgn T [ w σ ] +I . (16) ANNEXES  ∑ ∑ i ij l ( j +l) iN n2 2  j ≠i1=1 1l = −n contiguous, the windows are narrow and the average (t +1)   (t) σ = + sgn T [ w σ ] +I . (16) (16) i ∑ ij ∑ l ( j +l) i(t +1) H (x,0) model, the average of the real earth, has the highest 2x2 2σ j ≠i =1 l = −nni   | σ | = ωµ | | , (A.1) a 0(t) (t) possible resolution. If, on the other hand, the chosen E (x,0)(t) (t +1) σ = w σ . (15) y∑j l ( j +l)σ σ ij Following Zhang and Paulson (1997), we will refer to points are wide apart, the windows are themselves wide, l = −n σ(x,z) (t) nthis as the regularized version of the HANN or RHANN. for z and z tend to separate. As the windows widen, the 1 σ 1 2 j (t) (t) δE (x,0)w = ,k = −n,...,n (t+l) δH (x,0) σ = w σ . (15) yk The updated conductivities s are obtained from averages models tend to fatten, as they should. x This intuitively jj ∑ l ( j +l) N n δ | σ | = −2 | σ | Re[ − ]. (A.2) 2n +1 1   a a1 1 (t)l = −n(t +1) (t)of the original values s . The choice of the appropriate appealing feature of spatial a E (x,0v)erages is bH (x,uilt into equation 0)jw = ,k = −n,...,n y x σ = + sgn T [ w σ ] +I . (16) k  i ∑ ij ∑ l ( j +l) i 2n +1 flter depends on the application. Zhang and Paulson (1997) (18) in spite of its simplicity. Figures 1C and 1D show 2 2 j ≠i =1 l = −nN n  1 1(t +1) (t) (t) N nfound binomial flters useful to stabilize their solution. this feature developing when considering averages taken σ = + (t +sgn{1) [(1 −b)T σ + β T w σ  ] +I }. (17) 1 E (x,z) i (t +1) ∑ ij j ∑(t) ij l ( j +l) i yσ N n σ =i + sgn T [ w σ ] +I . (16) 2 2 In our case, we use boxcar functions, for this leads to the from data values 1 and 5 periods apart, respectively, for 1 1 i j ≠i =1 ∑ ij ∑ l = −n( j +l) i(t +1) (t) (t)σ 2 =2 + sgn{ [(1 −b)T σ + β T w σ ] +I }. (17) (t) H (x,z) j ≠i =1 l = −ni ∑ ij j ∑ ij l ( j +l) isimplest and most intuiti ve type of av erages. The flter, with x the synthetic data presented in Fig. 1A. σ 2 2j j ≠i =1 l = −n-1(t +1) 2n+1 uniform weights w = (2n+1) , k = -n,..., n (square k E (x,z) σ yσ (T ) δ (T ) − σ (T ) δ (T )i a 2 a 2 a 1 a 1 1 < σ (z ,z ) >= , (18) window), can be of different widths. We introduce a further A corresponding sequence of full fltered models ( β =1) 1 2 2w = ,k = −n,...,n (t) ∇ E −i ωµ σ (x,z)E = 0, (A.3) k δ (T ) − δ (T )σ (T ) δ (T ) − σ (T ) δ (T ) y 0 yσ a 2 a 1parameter β to accommodate a traditional trade-ofa 2 a 2 a 1 a 1 f between is presented in Figs. 1E and 1F when the RHANN algorithm j 2n +1 < σ (z ,z ) >= , (18) 1 2 z< σ (z ,z ) > two contrasting features. We modifed the dynamics to is applied to the same set of data. It can be observ − 2i ed that δ (T ) − δ (T ) 1 21 a 2 a 1 δ E (x,z) = E(0)e , (A.4) yw = ,k = −n,...,n N n the behavior of the averages follows that of the averages kz = 0.707 δ (T ),i = 1,2 < σ (z ,z ) > 1 1i a i (t +1) (t) (t)2n +1 1 2 σ = + sgn{ [(1 −b)T σ + β T w σ ] +shoI }.wn in Figs. 1C and 1D. (17) 2 There are some differences 2i ∑ ij j ∑ ij l ( j +l) i δ = . (A.5) z = 0.707 δ (T ),i = 1,2 2z = z z j ≠i =1 l = −ni a i that can be traced back to basic differences between the 1 2 ωµ σ0N n1 1 two approaches. One is that equation (18) is based on the (t +1) (t) (t)z = z z 1 2 H (x,z) σ = + sgn{ [(1 −b)T σ + β T w σ ] +I }. (17) x (17)i ∑ ij j ∑ ij l ( j +l) i original Niblett-Bostick approximation that assumes a ANNEXES σ (T ) δ (T ) − σ (T ) δ (T )2 2 a 2 a 2 a 1 a 1 j ≠i =1 l = −n ∇ × E = −i ωµ H < σ (z ,z ) >= , (18) 0boxcar function as the kernel of equation (5). The RHANN 1 2
ANNEXES δ (T ) − δ (T ) a 2 a 1When β =1 we get the full fltered solution. On the other hand, algorithm uses the kernel as indicated in equation (5) which H (x,0) 2x | σ | = ωµ | | , (A.1) σ (T ) δ (T ) − σ (T ) δ (T )when a β =0 we return to the original unfltered solution0 given includes a small negative sidelobe (Gómez-Treviño, 1987b) < σ (z ,z ) >a 2 a 2 a 1 a 11 2 < σ (z ,z ) >= E (x,0) , (18) H (x,0)1 2 y 2xby equation (14). In this way, by varying this parameter we and extends to infnity. This is a somewhat less restrictive | σ | = ωµ | | , (A.1) δ (T ) − δ (T )z = 0.707 δ (T ),i =a1,2 0a a 1i a iσ(x,z) E (x,0)can gradually control the effect of the averaging process. approximation than the rather simpler boxcar function. The y
< σ (z ,z ) > 1 2 other difference originates in the fact that the windows in z = z z 3 δE (x,0)σ(x,z) δH (x,0)1 2 y x δ | σ | = −2 | σ | Re[ − ]. (A.2) each approach are necessarily different, because in the frst z = 0.707 δ (T ),i =a 1,2 ai a i E (x,0) HδE(x(,x0,)0) δH (x,0)y x y xAPPLICATIONS case they cannot be uniform, for they are determined by the δ | σ | = −2 | σ | Re[ − ]. (A.2) a aANNEXES z = z z 1 2 E (x,0) H (x,0) data, and in the second case they are uniform for all depths. y x
1-D averages The overall effect is a somewhat better performance for E (x,z)y H (x,0)x 2 the approximation given by equation (5) judging by the ANNEXES | σ | = ωµ | | , (A.1) H (x,z) E (x,z) a 0yx Before proceeding to the application in 2-D of the abo E (x,0) ve resemblance of the average models to the original model y
E (x,z) H (x,z) numerical averaging technique, it is convenient to illustrate from which the data were drawn.y x H (x,0)σ(x,z) x 2 | σ | = ωµ | | , (A.1) 2 a 0 E (x,z) ∇ E −i ωµ σ (x,z)E = 0, (A.3) y E (x,0)y 0 y δE (x,0)y δH (x,0)y x2 δ | σ | = −2 | σ | Re[ − ]. (A.2) za a ∇ E −i ωµ σ (x,z)E = 0, (A.3) − 2iσ(x,z) y 0 yE (x,0) H (x,0)δ y x E (x,z) = E(0)e , (A.4) y zδE (x,0) − 2i 19Geologica Acta, 8(1), 15-30 (2010) δH (x,0)y x δ δ | σ | = −2 | σ | Re[ E2 (x,z) =−E(0)e ]., (A.2) (A.4) DOI: 10.1344/105.0000015132a a y δ = . (A.5E (x,0) H (x,0)E (x,z) y xy ωµ σ 20 2 δ = . (A.5) H (x,z) H (x,z) x ωµ σx 0
E (x,z) y E (x,z) ∇ × E = −i ωµ H yH (x,z) 0 x
H (x,z) 2x ∇ × E = −i ωµ H ∇ E −i ωµ σ (x,z)E = 0, (A.3) 0 y 0 y
E (x,z) y z − 2i
δ2 E (x,z) = E(0)e , (A.4) ∇ E −i ωµ σ (x,z)E y = 0, (A.3) y 0 y
z 22− 2i 3 δ = . (A.5) δ E (x,z) = E(0)e , (A.4) y ωµ σ0 3
2H (x,z) 2x δ = . (A.5)
ωµ σ0∇ × E = −i ωµ H 0
H (x,z) x
∇ × E = −i ωµ H 0

3
3 J. RODRÍGUEZ et al. 2-D Niblett-Bostick
A B
2 3 4 5-4 -2 0 2 4 10 10 10 1010 10 10 10 10
Depth z (m)Period T(s)
60 60
C D
5050 T -TT -T i i+5i i+1
40 40
30 30
20 20
10 10
3 4 5 2 3 4 52 10 10 10 10 10 10 1010
Depth z (m)Depth z (m)
1000 1000
E F
800 800
600 600
400 400
200 200
Windowsize Windowsize
2 3 4 52 3 4 510 10 10 10 10 1010 10
Depth z (m) Depth z (m)
FIGURE 1 a) synthetic data for test model. B) Test model. C) Recovery of test model using the average approach given by equation (18). T and T are
1 2
chosen as follows: T -T . D) using T -T . Notice the smoothing effect on the models as the averaging windows widen. The windows, shown below the
i i+1 i i+5
model, are nonuniform because they depend on the data whose index is shown in the right axis. e) and F) show the recovery of full fltered test models
(β = 1) using a 500 layer RHaNN. in this case, the window size is chosen at will corresponding to a uniform window size in log depth (i.e., the same
averaging window size for the 500 layers). Notice the same overall behavior of the average models as compared with those in Figures 1C and 1D,
respectively.
More appealing from the practical point of view, is the profle. This is illustrated in Fig. 2, which presents the results
stabilizing property of the averages when applied to noisy data. of applying equation (18) to the noisy data shown in Fig. 2A. Figure1When equation (18) is applied to noisy data using contiguous The natural regularizing effect that the wider windows have
data points, which is equivalent to using the original Niblett- on the recovered profles can be observed in Figs. 2C and 2D.
Bostick transformation in terms of slopes, the estimation of The results of applying the RHANN algorithm to the same
the conductivity profle is very unstable. This is equivalent set of data are shown in Fig. 2E and 2F. The sequence of
to the use of a narrow window, which understandably leads models illustrates the evolution of the process as the windows
to a noisy recovery of a conductivity profle. Increasing the are widened, reaching a reasonably good recovery with the
width of the window improves the stability of the recovered widest window.
20Geologica Acta, 8(1), 15-30 (2010)
DOI: 10.1344/105.000001513
� (S/m)
� (S/m) �a(S/m)
-3 -2 -1 -3 -2 -1
-3 -2 -1
10 10
10 10 10 10
10 10 10
�a index number
neural model index number
� (S/m) � (S/m) �(S/m)
-3 -1
-2
-3 -2 -1 -3 -2 -1
10 10 10
10 10 10
10 10 10
�a index number
neural model index numberJ. RODRÍGUEZ et al. 2-D Niblett-Bostick
The above results show that the RHANN algorithm to noisy data. In both instances, the main beneft of the
mimics the performance of equation (18) when applied RHANN algorithm is that it can be readily applied to
to 1-D sounding data. First, we have that the simulated higher dimensions.
averages actually behave as true averages, with the added
beneft that the RHANN algorithm allows more freedom in Two-dimensional averages
choosing the size of the windows, for they can be chosen
at will and are independent of how the data were sampled. The problem at hand can be summarized as follows:
Second, we have the same type of regularization effect as evaluate the performance of equations (16) and (17) as they
provided by equation (18), when the algorithm is applied apply to the solution of equation (6), where the elements of
BA
2 3 4 5-4 -2 0 2 4 10 10 10 1010 10 10 10 10
Depth z (m)Period T(s)
60 60
C D
5050
T -Ti i+1 T -Ti i+5
40 40
30 30
20 20
10 10
2 3 4 5 2 3 4 510 10 10 10 10 10 10 10
Depth z (m) Depth z (m)
10001000
E F
800800
600600
400 400
200 200
Windowsize Windowsize
2 3 4 52 3 4 5 10 10 10 1010 10 10 10
Depth z (m) Depth z (m)
FIGURE 2 a) synthetic data with 5% noise added to test model. B) Test model. C) Recovery of test model using the average approach given by equation
(18). T and T are chosen as follows: T -T . D) using T -T . Notice the smoothing effect on the models as the averaging windows widen. The windows, 1 2 i i+1 i i+5
shown below the model, are non-uniform because they depend on the data whose index is shown in the right axis. Figures e) and F) show the recovery
of full fltered test models (β = 1) using a 500 layer RHaNN. in this case, the window size is chosen at will, corresponding to a uniform window size in
log depth (i.e., the same averaging window size for the 500 layers). Notice the same overall behavior of the average models as compared with those
in Figures 2C and 2D, respectively.
Figure2
21Geologica Acta, 8(1), 15-30 (2010)
DOI: 10.1344/105.000001513
� (S/m) ���(S/m) �a(S/m)
-3 -2 -1 -1 -3 -2 -1
-3 -2
10 10 10 10 10 10 10 10 10
�a index number
neural model index number
�(S/m)
� (S/m)
���(S/m)
-3 -2 -1
-3 -2 -1
-3 -2 -1
10 10 10 10 10 10
10 10 10
�a index number
neural model index numberJ. RODRÍGUEZ et al. 2-D Niblett-Bostick
the matrix are computed for 2-D on the basis of equation
(8), as developed in Appendix C. The evaluation is effected
frst using synthetic data sets and then employing the
standard COPROD2 (Jones, 1993) data set. For the case of
synthetic data, we generate a 2-D mesh array discretization
consisting of 21 lines and 60 columns in which we
embedded resistive and conductive blocks as shown in Fig.
3A. Figures 3B and 3C present average models derived
by convolving the original model with the two windows
indicated in each case. The object of the exercise is to
recover the average models from a given set of synthetic
data. Twenty-nine soundings at the indicated sites were
simulated in the period range of 0.01 to 100 s. The TM
and TE impedances were converted to corresponding series
and parallel quantities using equations (C7) and (C8).
Figure 4A shows the resultant model obtained when
using the RHANN algorithm. This result was achieved
after a ‘stable’ state of the dynamics was reached, when
a uniform half-space of resistivity ρ=10 Ω·m was used as
an initial condition for each neuron state. Figure 4B shows
the model obtained when a different uniform half-space FIGURE 4 a) Resultant model obtained when using the regularized
Hopfeld artifcial neural network (RHaNN). This result was achieved of ρ=1 Ω·m was used as an initial condition. Graphics 4A
after a stable state of the dynamics was reached, when uniform half-
ρspace of resistivity = 10 Ω·m was used as an initial condition for
each neuron state. B) model obtained when a different initial condition
ρis used. in this case a half-space of = 1 Ω·m was used. The above
graphics were obtained using both series and parallel data.
and 4B were obtained using both series and parallel data.
It can be observed that the models are practically identical,
as they should be, for the need of an initial state is an
internal requirement of the HANN algorithm and not of the
formulation of the inverse problem as given by equation
(6). It can also be noted that the conductive block is better
resolved than the resistive one, as could be expected from
the known features of the Niblett-Bostick approximation
as well as from the basic properties of electromagnetic
induction in general. We found that it is still possible to
emphasize either the resistor or the conductor by choosing
the data to be inverted. The series response is more
representative of resistors while conductors are better
refected in the parallel data (Romo et al., 2005). This is
illustrated in Fig. 5 for the case of the resistor and in Fig.
6 for the case of the conductor. In the preceding fgures,
and in the ones that follow, the rms values shown on top
of the models were computed as 100 times the square root
of the mean of the residuals squared. The residuals are
simply the difference between the data and the computed
response of the model given by equation (6), normalizing
each difference by the corresponding computed response.
We now turn to illustrate the performance of the
FIGURE 3 a) Two-dimensional test model. B) and C) average models
algorithm with feld data. The data corresponds to the obtained when the window shown at the right is applied to the original
model. These average models are the target of the RHaNN approach. original COPROD2 (Jones, 1993) data set, with no
22Geologica Acta, 8(1), 15-30 (2010)
DOI: 10.1344/105.000001513J. RODRÍGUEZ et al. 2-D Niblett-Bostick
FIGURE 7 model obtained by means of the Rodi and mackie (2001)
non-linear inverse method adapted for series and parallel data (Romo
et al., 2005) when applied to the CopRoD2 (Jones, 1993) standard data FIGURE 5 model obtained by means of the RHaNN, when using data
set.from the series mode only with a 3×3 window size and β = 1.0 . it
can be seen that the series-based inversion emphasizes the resistivity
nature of the anomalies.
each unknown, the 748x2,310 matrix took 20 min to
compute in an AMD Athlon 1.2 GHz processor, with 20
extra minutes for the minimization. Figure 8A shows the
result of jointly inverting the series and parallel data. It
can be observed that the approximation clearly recovers
FIGURE 6 model obtained by means of the RHaNN, when using data
from the parallel mode with only a 3×3 window size and β = 1.0. it can
be seen that the parallel-based inversion emphasizes the conductive
nature of the anomalies.
corrections for static effects. The full tensor was used to
compute the invariant series and parallel impedances using
equations (C3) and (C4). We used 34 sounding sites with
22 periods each, for a total of 748 apparent conductivity
values. The range of periods was from 1.3 s to 1365 s. This
same data set was recently inverted by Romo et al. (2005)
using Rodi and Mackie (2001) inverse code adapted for
series and parallel data. The model shown in Fig. 7 presents
two main conductive anomalies. The deeper anomaly to
the left of the model corresponds to the North American
Central Plains (NACP) conductivity anomaly, which has
been discussed in detail in the magnetotelluric literature
(Jones, 1993). The smaller lobular anomaly towards the
right end of the profle corresponds in turn to the Thomson
Nickel Belt (TOBE) anomaly.
FIGURE 8The challenge for the proposed 2-D Niblett-Bostick models obtained by the RHaNN algorithm when applied to
the CopRoD2 (Jones, 1993) standard data set. a) Result obtained using approximation is the recovery of these two anomalies. The
data from series and parallel modes. B) Result obtained using data from
model was discretized using a 36×70 mesh array for a total the series mode only. C) Result obtained using data from the parallel
of 2,310 conductivity unknowns. Using 22 bits to represent mode only.
23Geologica Acta, 8(1), 15-30 (2010)
DOI: 10.1344/105.000001513J. RODRÍGUEZ et al. 2-D Niblett-Bostick
only the NACP anomaly. The TOBE anomaly is somehow
diluted towards the right side of the model. The results are
somewhat similar when only series data are inverted, as
can be observed in Fig. 8B. The main difference is that in
this last case, the NACP anomaly decreases to about half
the size because of the removal of the parallel data which,
as mentioned above, tends to emphasize conductors. This
effect is explicitly shown in Fig. 8C, where only parallel
data were used in the inversion. Both conductive anomalies
are now clearly manifested, in reasonable agreement with
the model shown in Fig. 7, which does not depend on the
approximation. Notice also that, of the three models of Fig.
8, the one for the parallel data has the least rms misft. Still,
a level of 52% misft might seem too large. A corresponding
model with a 22% misft is shown in Fig. 9A. This was
obtained simply by lowering the trade-off parameter from
1 to ½. As expected, the higher-resolution model now
includes smaller scale anomalies, but the same two main
conductors still dominate the picture. Thus, summarizing,
when in search of conductors as in most MT surveys, we
should feed the approximation with parallel data, which
emphasize conductors.
FIGURE 10 zoom view of the frst 5 km for the models shown in Figure 9.
To test the robustness of the algorithm we perturbed
the data by multiplying each sounding curve by a random
factor of 60% to simulate the static shift effect produced by
small local variations in conductivity around the sounding model is the same high-resolution model described in the
sites. The results are shown in Fig. 9B. The reference previous paragraph, it is shown in Fig. 9A and has a misft
of 22%. When the soundings were perturbed by 60% the
best ft achieved was 40%, much better than an expected
FIGURE 9 models obtained by the RHaNN algorithm when random noi-
se was added to the individual sounding curves of the CopRoD2 stan-
dard data set to simulate static shift effects. The models were obtained
using data from the parallel mode. a) 0% of random noise added to
data. B) 60% of random noise added to data. FIGURE 11 same as for Figure 9 but using a wider window.
24Geologica Acta, 8(1), 15-30 (2010)
DOI: 10.1344/105.000001513