16
pages

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

Vous aimerez aussi

de rheinisch-westfalischen_technischen_hochschule_-rwth-_aachen

de rheinisch-westfalischen_technischen_hochschule_-rwth-_aachen

de technische_universitat_munchen

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

2π

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) 2a ∫ (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