Joint MT and CSEM data inversion using a multiplicative cost function approach

Joint MT and CSEM data inversion using a multiplicativecost function approach

A. Abubakar1, M. Li1, G. Pan1, J. Liu1, and T. M. Habashy1

ABSTRACT

We have developed an inversion algorithm for jointly

inverting controlled-source electromagnetic (CSEM) data

and magnetotelluric (MT) data. It is well known that

CSEM and MT data provide complementary information

about the subsurface resistivity distribution; hence, it is

useful to derive earth resistivity models that simultane-

ously and consistently fit both data sets. Because we are

dealing with a large-scale computational problem, one usu-

ally uses an iterative technique in which a predefined cost

function is optimized. One of the issues of this simultane-

ous joint inversion approach is how to assign the relative

weights on the CSEM and MT data in constructing the cost

function. We propose a multiplicative cost function instead

of the traditional additive one. This function does not

require an a priori choice of the relative weights between

these two data sets. It will adaptively put CSEM and MT

data on equal footing in the inversion process. The inver-

sion is accomplished with a regularized Gauss-Newton

minimization scheme where the model parameters are

forced to lie within their upper and lower bounds by a non-

linear transformation procedure. We use a line search

scheme to enforce a reduction of the cost function at each

iteration. We tested our joint inversion approach on syn-

thetic and field data.

INTRODUCTION

The controlled-source electromagnetic (CSEM) method has

the potential of providing useful information in applications

such as offshore hydrocarbon exploration. With a horizontal

electric dipole transmitter towed by a ship and multicomponent

electromagnetic receivers on the seafloor, this method has been

applied in several field surveys. The high contrast in resistivity

between saline-filled rocks and hydrocarbons (oil and gas)

makes this method well suited for detecting resistive targets

such as hydrocarbon reservoirs; see Constable et al. (1986),

Srnka (1986), Chave et al. (1991), MacGregor and Sinha

(2000), Eidesmo et al. (2002).

To make use of CSEM data one has to use nonlinear inver-

sion approaches, such as the Gauss-Newton method (Abubakar

et al., 2008; Abubakar et al., 2009), (preconditioned) nonlinear

conjugate gradient methods (Commer and Newman, 2008a) or

the limited memory quasi-Newton method (Plessix and Mulder,

2008). A history and brief technical overview of this CSEM

technique can be found in Constable (2010).

The magnetotelluric (MT) method has a long history of use in

diverse applications such as mining and geothermal explorations.

It also has been applied in offshore hydrocarbon explorations,

e.g., Constable et al. (1998), Hoversten et al. (1998), Constable

et al. (2009), Smirnov and Pedersen (2009). MT data largely are

insensitive to thin, high resistivity targets that are associated with

hydrocarbon deposits trapped in thin planar sedimentary layers.

This is because MT signals, which have characteristics of hori-

zontally polarized plane waves, are more sensitive to conductive

targets. However, it is known that MT has a larger depth of

investigation than CSEM. Hence, MT can be used effectively to

construct resistivity estimates of the subsurface (i.e., the back-

ground medium). MT recently has been used for improving seis-

mic depth imaging workflows (Colombo and De Stefano, 2007;

Virgilio et al., 2009). Similar to CSEM, interpretation of MT data

has to be done using full nonlinear inversion approaches such as

those described by Rodi and Mackie (2001), Kumar et al. (2007),

Abubakar et al. (2009).

MT and CSEM data can provide complementary information:

MT mainly provides information on background resistivity struc-

tures, whereas CSEM identifies thin resistive targets. Constable

and Weiss (2006) also point out the practical importance of MT

data: The same receivers can collect CSEM and MT data; hence

Manuscript received by the Editor 27 October 2010; revised manuscript received 00 0000; published online 5 May 2011.1Schlumberger-Doll Research, Cambridge, Massachusetts, U.S.A. E-mail: [email protected]; [email protected]; [email protected]; [email protected];

[email protected] 2011 Society of Exploration Geophysicists. All rights reserved.

F203

GEOPHYSICS, VOL. 76, NO. 3 (MAY-JUNE 2011); P. F203–F214, 13 FIGS., 1 TABLE.10.1190/1.3560898

Downloaded 06 May 2011 to 199.6.131.16. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

MT data come at a relatively low cost because they can be

recovered from time-series data when the source is turned off.

Usually, CSEM and MT data are inverted separately. In the

joint interpretation workflow, MT data first is inverted and then

the results are used as an initial model for CSEM data inversion.

Mackie et al. (2007) propose a joint inversion method where

one inverts both data sets simultaneously. By using this

approach, they can obtain reasonable results with a reduced

number of sources (for CSEM) and receivers. Commer and

Newman (2008b) also show that one can significantly improve

the inversion results by using this simultaneous joint inversion

approach. However, Mackie et al. (2007) and Commer and

Newman (2008b) point out that in this joint inversion approach

one has to carefully choose the relative weighting between MT

and CSEM data.

We propose to use a multiplicative cost function for the si-

multaneous joint inversion of MT and CSEM data. This multi-

plicative cost function will effectively put equal weighting for

MT and CSEM data. This approach has been used previously

for multifrequency electromagnetic inversions; see Abubakar

and van den Berg (2000). Unlike Mackie et al. (2007) and

Commer and Newman (2008b), we use the regularized Gauss-

Newton approach as our inversion scheme. The general descrip-

tion of the forward and inversion algorithms can be found in

Abubakar et al. (2008). The forward algorithm is implemented

based on the theory described in Druskin and Knizhnerman

(1994) and Ingerman et al. (2000).

The inversion algorithm uses the Gauss-Newton framework

described in Habashy and Abubakar (2004) augmented with the

multiplicative regularization technique introduced in van den

Berg et al. (1999). The inversion algorithm is equipped with

two regularization functions to produce either a smooth (using a

standard L2-norm function) or a blocky (using a weighted L2-

norm function) conductivity distribution (van den Berg and

Abubakar, 2001).

To enhance the robustness of the algorithm, we incorporated

a nonlinear transformation for constraining the minimum and

maximum values of the conductivity distribution and a line-

search procedure for enforcing the error reduction in the cost

function in the optimization process; see Habashy and Abubakar

(2004). To illustrate the performance of our inversion method,

we consider synthetic data inversion using isotropic and trans-

verse isotropic (TI) media.

We also will show the inversion results of field data collected

in the Norwegian shore. The survey is carried out between the

Møre basin and the Møre margin in the north Norwegian Sea.

The geological environment is quite complex, covering the mar-

gin between oceanic and continental crust (the Møre margin),

and the western part of the Møre basin, in which Eocene basalts

are present within the section; see Blystad et al. (1995). As

pointed out by Virgilio et al. (2009), one of the objectives of

this survey is to determine the thickness of prebasalt sediments,

and, if possible, identify the crystalline basement.

INVERSION METHOD

The nonlinear inverse problems are described by the follow-

ing operator equations:

dCSEM ¼ sCSEM mð Þ; (1)

dMT ¼ sMT mð Þ; (2)

where

dCSEM ¼ d CSEM rSi ; r

Rj ;xk

� �; i ¼ 1; 2;…; I;

h

j ¼ 1; 2;…; J; k ¼ 1; 2;…;K1�T

is the vector of CSEM data and

dMT ¼ dMT rRj ;xk

� �; j ¼ 1; 2;…; J; k ¼ 1; 2;…;K2

h iT

(4)

is the vector of MT impedance data, in which rSi , rR

j , and xk are

the source positions, the receiver positions, and the frequency of

operations, respectively. The superscript T denotes the transpose

of a vector. The vectors sCSEM mð Þ and sMT mð Þ represent simu-

lated data computed by solving Maxwell equations using a fi-

nite-difference method in the frequency domain.

In CSEM the source is an electric dipole oriented parallel to

the tow line, and the receiver can be any component of the mag-

netic or electric field vector. In MT the source is a plane wave

and the receivers measure the horizontal components of electric

and magnetic fields.

In this work we use the 2.5D approximation for CSEM data,

see Abubakar et al. (2008), and 2D approximation for MT data.

The forward algorithm uses the framework described in Druskin

and Knizhnerman (1994); however, we use a multifrontal LU

decomposition (Davis and Duff, 1997) to invert the stiffness ma-

trix. By using this direct matrix inversion technique, we can

simulate multisource experiments at nearly the cost of simulat-

ing only one single-source experiment. This feature is quite im-

portant for the inversion because it uses data from more than

one source position=orientation. The direct matrix inversion

technique is accomplished using the optimal grid technique

(Ingerman et al. 2000) to extend the boundaries of the computa-

tional domain to infinity, and a diagonal anisotropic material

averaging formula described in Keller (1964) to assign appropri-

ate conductivity values on the finite-difference grid nodes.

The unknown vector of model parameters is defined as:

m ¼ mh x‘; zq

� �;mv x‘; zq

� �; ‘ ¼ 1;…; L; q ¼ 1;…;Q

� �;

(5)

where x‘ and zq denote the center of the discretization cell.

The unknown model parameters mh rð Þ ¼ rh rð Þ=r0 rð Þ and

mv rð Þ ¼ rv rð Þ=r0 rð Þ are the normalized horizontal and vertical

conductivities where r0 rð Þ is the conductivity distribution of the

initial model used in the inversion. In this paper we assume that

the conductivity is TI-anisotropic (defined as rh ¼ rxx ¼ ryy,

rv ¼ rzz with all other components of the conductivity tensor

set to zero, where x, y, and z denote the Cartesian axes). Follow-

ing Abubakar and van den Berg (2000) and Abubakar et al.

(2008), we simultaneously invert CSEM and MT data using a

multiplicative cost function defined as

Un mð Þ ¼ /CSEM mð Þ � /MT mð Þ � /Rn mð Þ; (6)

F204 Abubakar et al.


where the CSEM data cost function is given by

/CSEM mð Þ ¼ 1

2K1

XK1

k¼1

PIi¼1

PJj¼1

wCSEMi;j;k dCSEM

i;j;k � sCSEMi;j;k mð Þ

h i�� 2

PIa¼1

PJb¼1

wCSEMa;b;k dCSEM

a;b;k

�� 2

¼ W CSEM dCSEM � sCSEM mð Þ� � 2

(7)

and the MT data cost function is given by

/MT mð Þ ¼ 1

2K2

XK2

k¼1

PJj¼1

wMTj;k dMT

j;k � sMTj;k mð Þ

h i�� 2

PJb¼1

wMTb;k dMT

b;k

�� 2

¼ WMT dMT � sMT mð Þ� � 2

; (8)

in which wCSEMi;j;k and wMT

j;k are real-valued data weighting matri-

ces. The normalized data weighting matrices WCSEM and WMT

are given by

W CSEMi;j;k ¼

wCSEMi;j;k

2K1

PIa¼1

PJb¼1

wCSEMa;b;k dCSEM

a;b;k

�� 2

26664

37775

12

(9)

and

WMTj;k ¼

wMTj;k

2K2

PJb¼1

wMTb;k dMT

b;k

�� 2

26664

37775

12

: (10)

In our implementation, wMTj;k is chosen to be an identity matrix

while wCSEMi;j;k is based on the so-called Jacobian weighting defined as

wCSEM ¼ diag JCSEM0 JCSEM

0

� �Hh in o�1

2; (11)

where JCSEM0 is the CSEM Jacobian matrix of the initial model.

The regularization cost function is given by

/Rn mð Þ ¼

Xj¼x;z

ðD

b2j;n x; zð Þ oj ln mh x; zð Þ½ �j j2

n

þ oj ln mv x; zð Þ½ �j j2 þ d2n

odxdz; (12)

in which ox and oz denote spatial differentiation with respect to

x and z. The weights bj;n x; zð Þ are given by

b2j;n x; zð Þ ¼ 1Ð

D oj ln mh;n x; zð Þ� �� 2 þ oj ln mv;n x; zð Þ

� �� 2 þ d2n

n odxdz

(13)

for the L2-norm regularizer and

b2j;n x; zð Þ ¼ 1

V

1

oj ln m h;n x; zð Þ� �� 2 þ oj ln m v;n x; zð Þ

� �� 2 þ d2n

;

V ¼ð

D

dxdz (14)

for the weighted L2-norm regularizer as introduced in van den

Berg and Abubakar (2001). In the above equations, D denotes

the inversion domain and mn denotes the unknown parameter at

iteration nth.

The L2-norm regularizer is known to favor smooth profiles,

while the weighted L2-norm regularizer is known for its ability

to preserve edges. Note that for the L2-norm regularizer, the

weight bj;n x; zð Þ is independent of the spatial position. This

weighted L2-norm regularization factor in equations 12 and 14

belongs to the same class as the well-known total variation reg-

ularization (Vogel and Oman, 1996; Charbonnier et al., 1997;

Farquharson and Oldenburg, 1998; Farquharson, 2008). The pos-

itive parameter dn is chosen to be

d2n ¼

/CSEM mnð Þ/MT mnð ÞDxDz

; (15)

where Dx and Dz are discretization cell widths. A discussion

about the choice of dn can be found in van den Berg et al.

(2003).

To minimize Un mð Þ in equation 6, we use a Gauss-Newton

minimization framework described in Habashy and Abubakar

(2004). At the nth iteration, we obtain a set of linear equations

for the search vector pn that determines the minimum of the

approximated quadratic cost function, namely,

Hnpn ¼ �gn; (16)

where the gradient vector gn is given by

gn ¼ /MT mnð Þr/CSEM mð Þ�m¼mn

þ/CSEM mnð Þr/MT mð Þ�m¼mn

þ/CSEM mnð Þ/MT mnð Þr/Rn mð Þ�m¼mn

¼ /MT mnð Þ JCSEMn

� �H WCSEM� �TWCSEM dCSEM � sCSEM mnð Þ

� �n o

þ/CSEM mnð Þ JMTn

� �H WMT� �TWMT dMT � sMT mnð Þ

� �n o

þ/CSEM mnð Þ/MT mnð ÞL 1ð Þ mnð Þmn; (17)

where H denotes the transpose conjugate of a matrix. In the

above equation we used the fact that /Rn mnð Þ ¼ 1. The matrices

JCSEMn and JMT

n are the Jacobian matrices and their explicit

expression can be found in Abubakar et al. (2008). The gradient

of the regularization cost function is given by

L 1ð Þ mnð Þvh i

‘;q¼ 1

mn;‘;q

Xj¼x;z

ojbj;n oj ln vð Þ½ � �

‘;q: (18)

The Hessian matrix is given by

Hn ¼ /MT mnð Þ JCSEMn

� �H WCSEM� �TWCSEMJCSEM

n

h iþ /CSEM mnð Þ JMT

n

� �H WMT� �TWMTJMT

n

h iþ /CSEM mnð Þ/MT mnð ÞL 2ð Þ mnð Þ: (19)

F205Joint MT and CSEM data inversion


In the above equation, to make the Hessian matrix nonnega-

tive definite, we neglect the second-order derivatives of the cost

function /CSEM mð Þ and /MT mð Þ as well as the nonsymmetric

terms. The second derivative of the regularization cost function

is given by

L 2ð Þ mnð Þ � vh i

‘;q¼ 1

mn;‘;q

Xj¼x;z

ojbj;n ojv

mn

� � �� ‘;q

:

(20)

For this multiplicative cost function the gradient is exact;

however, the Hessian matrix is an approximate one because we

neglected the second-order derivatives and the nonsymmetric

terms.

After the Gauss-Newton search vector pn is obtained by solv-

ing the linear system of equation 16 using a conjugate gradient

least-squares (CGLS) technique, Golub and Van Loan (1996),

the unknown model parameters are updated using a nonlinear

transformation and the line-search procedure described in Haba-

shy and Abubakar (2004).

This multiplicative cost function in equation 6 will make MT

and CSEM cost functions equally important. This is achieved by

using the MT data cost function at the previous iteration

/MT mnð Þ as the weight for the CSEM cost function, and vice

versa; see equations 17 and 19.

This strategy is reasonable when we do not have a priori in-

formation on which data set is more reliable or should be domi-

nant in the inversion process. When such a priori information is

available, a constant relative weighting that puts emphasis on

one of the cost functions preferably should be used.

As for the weight for the regularization factor, a large weight

at the beginning of the optimization process is used because the

values of /CSEM mð Þ and /MT mð Þ still are large. In this case, the

search direction is predominantly a steepest descent, which is

more appropriate to use in the initial steps of the iteration pro-

cess because it has a tendency to suppress large swings in the

search direction. As the iteration proceeds, the optimization pro-

cess gradually will reduce the error in the data misfit while the

regularization factor /Rn mð Þ remains at a nearly constant value

close to unity.

In this case, the search direction corresponds to a Newton

search direction, which is more appropriate to use as we get

closer to the minimum of the data misfit cost function of either

/CSEM mð Þ or /MT mð Þ where the quadratic model of the cost

function becomes more accurate.

If noise is present in the data, the data misfit cost functions

/CSEM mð Þ and /MT mð Þ will plateau to values determined by

the signal-to-noise ratio; hence, the weight on the regularization

factor is nonzero. In this way, the noise will, at all times, be

suppressed in the inversion process and the need for a larger

regularization when data contain noise will be automatically

fulfilled.

The optimization process is terminated if one of the following

stopping criteria is satisfied: The value of the data misfit U gets

below a prescribed error quantity, or the difference between the

data misfit U at two successive iterates is within a prescribed

error quantity, or the change in the model m at two successive

iterates is within a prescribed error quantity, or the total of itera-

tions exceeds a prescribed maximum.

NUMERICAL RESULTS AND DISCUSSIONS

Isotropic synthetic data inversion

We present a synthetic data inversion example in which two

thin reservoirs are embedded in a nonflat layered formation. The

horizontal distance between these reservoirs is about 1 km. The

widths of the reservoirs are about 5 km and the resistivities of

the reservoirs are 50 Ohm-m. The resistivity distribution of this

test model is shown in Figure 1a, and the color bar is given in

terms of the logarithm of the resistivity.

The seabed consists of three nonflat layers. Their resistivity

values from top to bottom are 1.67, 3.33, and 6.67 Ohm-m.

The seawater depth is 840 m and the water resistivity is 0.33

Ohm-m. For the survey, we used 21 receivers located at the

seafloor. These receivers are deployed uniformly from

�10 km to 10 km. They record CSEM and MT data and they

are indicated by circle signs in Figure 1a. For CSEM data,

we invert only the inline fields (horizontal electric fields in

the tow direction and the sources are horizontal electric

dipoles also oriented in the tow direction), while for MT data

we use transverse electric (TE) and transverse magnetic (TM)

polarization impedances. The CSEM source is towed from

�10 to 10 km at 50 m above the receivers. In the inversion

we only use CSEM sources located every 500 m; hence, we

have 41 CSEM sources.

We use three CSEM frequencies (0.25, 0.75, and 1.25 Hz)

and eight MT frequencies (0.005, 0.01, 0.025, 0.05, 0.1, 0.25,

0.5, and 1 Hz). The synthetic data sets are generated using a 3D

frequency-domain forward modeling code described in Zaslav-

sky et al. (2006). For each frequency, we also added random

white noise of 2% maximum amplitude of all data points to

CSEM and MT data. Hence, the noises are frequency-depend-

ent; however, they are not offset-dependent.

In the inversion we use a domain of 36 by 11.16 km with

grid size 200 by 60 m. The number of unknowns is 33,480. For

this case we assumed that the resistivity distribution is isotropic;

hence we have m ¼ mh ¼ mv.

Figure 1. (a) True and (b) initial models of the isotropic test case.The color bars are in terms of the logarithm of the resistivity(Ohm-m).



For the initial model we only use a model with a homoge-

neous seabed as shown in Figure 1b. The homogeneous

seabed resistivity is 6.67 Ohm-m. Figure 2a shows the recon-

structed image using only MT data. We observe that the

three seabed layers are properly reconstructed in shape and

conductivity values; however MT inversion failed to detect

the presence of the reservoir. This reflects the limitation of

MT data where the low-frequency electromagnetic plane

waves can have a large depth of investigation and good hori-

zontal resolution; however they lack sensitivity to resistive

layers.

Note that the regions x < �10 km and x > 10 km are not

well reconstructed because we do not have receivers in these

regions. In addition, the resistivity of the seabed for depths

larger than 8 to 9 km is not well reconstructed. The latter is

because we already reached the limit of the MT data depth of

investigation at the frequency of operations.

The inversion result of CSEM data is given in Figure 2b. In

this case we reconstruct a rough estimate model of the reser-

voir; however, the layered formation is not reconstructed prop-

erly. This mainly is because of the limited depth of investiga-

tion of CSEM data. Because the layered formation is not well

reconstructed, the CSEM inversion tries to compensate by

putting relatively shallow conductive and resistive objects

close to the interface between the water and seabed. In addi-

tion, the depth and shape of the reservoirs are not very well

reconstructed.

Figure 3 shows the joint inversion results using CSEM and

MT data. In Figure 3a, we first invert MT data and then use the

result as an initial model for the CSEM data inversion. By com-

paring the inversion result in Figure 3a with the result in Figure

2b, we observe quite significant improvements. In this case, the

thicknesses of the reservoirs are better reconstructed; however,

the depths of both reservoirs are more or less the same.

In Figure 3b, we present the result of inverting MT and

CSEM data simultaneously. We observe that the reservoir and

the layers both are well reconstructed. This is because MT data

are sensitive to the layering structure and CSEM data mainly

are sensitive to the resistive layers. At first glance it seems that

both joint inversion approaches (sequential and simultaneous)

produce similar results. However, close examination of the

results reveals that the depths of the reservoirs are better recon-

structed by using the simultaneous joint inversion approach (i.e.,

they are not located at the same depth).

In Figure 4 we plot the true model and inverted models

obtained by using joint inversion approaches only on the region

�10 km < x < 10 km and 1 km < z < 4 km. From Figure 4, we

also observe that the second layer and the third layer of the

seabed are reconstructed better using the simultaneous joint

inversion algorithm. This also is apparent when we compare the

normalized data misfits (the cost functions without including the

Figure 2. Inverted models of the isotropic test case using (a) MTdata and (b) CSEM data.

Figure 3. Inverted models of the isotropic test case using jointMT and CSEM data inversion with (a) the sequential approachand (b) the simultaneous approach.

Figure 4. Zoom-in of the isotropic test case (a) true model, (b)the sequential approach, and (c) the simultaneous approach.



data weighting matrices) for every frequency of the joint CSEM

and MT data inversion run; see Table 1.

In the table we show also the number of iterations nð Þ used

by each inversion run. From this table we observe that the

simultaneous joint CSEM and MT data produce the lowest data

misfits. It seems that this joint inversion helps the CSEM

data inversion more than the MT data inversion, because the

MT data misfits of the separate MT inversion do not reduce as

much as the CSEM data misfits. This exercise also shows that

without the use of a proper background medium (either derived

from an MT inversion or from other independent measurements)

the CSEM data interpretation will not produce reliable inversion

results.

FORCE field data inversion

We present the inversion results of the so-called FORCE

data set. The survey is carried out between the Møre basin

and the Møre margin in the north Norwegian Sea as shown

in Figure 5.

The geological environment is quite complex, covering the

margin between oceanic and continental crust (the Møre mar-

gin), and the western part of the Møre basin, in which Eocene

basalts are present within the section; see also Virgilio et al.

(2009). One of the objectives of this survey is to determine the

thickness of prebasalt sediments, and, if possible, identify the

crystalline basement. More information on the geology of

the area can be found in Blystad et al. (1995). The survey line

extends from northwest to southeast. Fourteen MT stations were

used. The frequency of the acquired MT data ranges from

0.0003 to 0.1 Hz, among which eight frequencies are used in

our inversion from 0.02 to 0.09 Hz.

For MT data we use TE and TM polarization impedances. In

the CSEM survey, there are 14 receivers and 324 sources. The

receiver locations are indicated in Figure 6a by small white rec-

tangles, and they are located between x ¼ �15 km and

x ¼ 25 km. We only use data at three frequencies (0.125 Hz,

0.25 Hz, and 0.375 Hz) for the CSEM data inversion. The num-

ber of measurement data points is 2,928.

For the inversion domain, we choose a rectangular domain

with size 70 by 8.8 km. The grid size is 250� 50 m. Hence, the

total of unknown resistivity cells is 42,450. The initial model

used in the inversion is a homogeneous seabed model with resis-

tivity 2.5 Ohm-m as shown in Figure 6a. The unknown medium

is assumed to be an isotropic one. Note that the interface

between the water and seabed is not flat.

Figure 6b shows the inversion results using only MT data.

The inversion process stops after 24 Gauss-Newton iterations. In

this figure, we observe a large resistive body on the left (its size

extends beyond 8 km in depth). These results are consistent

with the sensitivity of the MT data (we have deep penetration

depth; however, we obtain lower resolution because of the very

low operating frequency). Next, we invert the CSEM data only

and the results are shown in Figure 6c. In this case the top of

the salt is reconstructed clearly, but the image of the deep

region beyond 4 km cannot be trusted because of the lack of

sensitivity of the CSEM data. We also observe that the CSEM

data has a better resolution because it operates at higher fre-

quencies. For the CSEM inversion we use 31 iterations.

Next, we invert CSEM and MT data either sequentially or

simultaneously. The reconstructed images are shown in Figure

7a and 7b. The CSEM data inversion using the MT data inver-

sion results as its initial model takes 31 iterations, whereas the

simultaneous joint inversion takes 22 iterations. In both joint

inversion results the salt body is well reconstructed. However,

in the simultaneous joint inversion result the presence of a resis-

tive body (which probably is an artifact) located under the salt

layer is significantly weaker than in the sequential joint inver-

sion result.

When we compare these inversion results with the results

from the seismic interpretation as shown in Figure 8a and 8b, or

with the known geology of the area as shown in Figure 8c, we

Table 1. Normalized data misfit values (in %) for the iso-tropic test model inversion.

CSEM MT Sequential Simultaneous

CSEM, 0.25 Hz 35.169 3.296 1.822

CSEM, 0.75 Hz 68.436 4.473 2.270

CSEM, 1.25 Hz 97.403 5.273 2.163

MT, 0.005 Hz 5.006 4.476

MT, 0.01 Hz 3.064 2.513

MT, 0.02 Hz 2.416 2.343

MT, 0.05 Hz 1.853 1.871

MT, 0.1 Hz 1.752 1.592

MT, 0.25 Hz 1.865 1.759

MT, 0.5 Hz 1.680 1.642

MT, 1 Hz 1.860 1.792

n 4 11 31 21 Figure 5. Survey outline of the FORCE data set. This plot isobtained from Virgilio et al. (2009).



observe good correlations between them, especially for the si-

multaneous joint inversion result. This shows that the joint MT

and CSEM inversion scheme using the multiplicative cost func-

tion might greatly improve the data interpretation.

In Figure 9, we show the comparison between MT measured

data and simulated data for the initial and inverted models for

TM and TE polarization at frequency 0.05 Hz. The top plots

in Figure 9 show the real part of the impedances whereas the

bottom plots show the imaginary part.

In Figure 10, we also show the comparison between CSEM

measured data and simulated data for the initial and inverted

models for data at 0.125 Hz, 0.25 Hz, and 0.375 Hz and for the

receiver located close to the center of the inversion domain.

Figures 9 and 10 show excellent agreement between the meas-

ured and simulated data of the inverted model.

TI-anisotropic synthetic data inversion

By using this final example, we discuss pros and cons of car-

rying out the joint MT and CSEM data inversion for TI-aniso-

tropic media.

For the synthetic test example, we consider a thin reservoir

embedded in a seabed that consists of three nonflat layers as

shown in 11a and 11b. The width of the reservoir is about 6

km and its depth is about 1.5 km below the water-seabed

interface.

The horizontal and vertical resistivities of the reservoir are 20

Ohm-m and 50 Ohm-m, respectively. The horizontal resistivities

of the layers (from top to bottom) are 1, 2, and 4 Ohm-m, while

their vertical resistivities are 1.33, 3.33, and 6.67 Ohm-m. The

thickness of the first and second (from the top) layers are about

1.5 and 1 km. The seawater depth is 1 km and the water resis-

tivity is 0.33 Ohm-m.

For the survey, we used 21 receivers located at the seafloor.

These receivers are deployed uniformly from �10 to 10 km.

They record CSEM and MT data and are indicated by circles in

Figure 11.

For CSEM data we invert only the inline fields (horizontal

fields in the tow direction) while for MT data we use TE and

TM polarization impedances. The CSEM source is towed from

�10 to 10 km at 50 m above the receivers. In the inversion we

only use CSEM sources located every 500 m; hence, we have

41 CSEM sources.

We use two CSEM frequencies (0.25 and 0.75 Hz) and eight

MT frequencies (0.005, 0.01, 0.025, 0.05, 0.1, 0.25, 0.5, and 1

Hz). The synthetic data sets are generated using the 3D fre-

quency-domain forward modeling code described in Zaslavsky

et al. (2006). For each frequency, we added random white noise

of 2% maximum amplitude of all data points to CSEM and MT

data.

In the inversion we use a domain of 36 by 6 km, with grid

size 200 by 75 m. The number of unknowns is 14,400. For the

Figure 6. Initial model (a) and the FORCE datainversion results from (b) MT data and (c)CSEM data. The color bars are in terms of thelogarithm of the resistivity (Ohm-m), and theaxes of the figures are in meters.



Figure 7. The FORCE data inversion results from(a) joint sequential MT and CSEM data and (b)joint simultaneous MT and CSEM data. The colorbars are in terms of the logarithm of the resistivity(Ohm-m), and the axes of the figures are in meters.

Figure 8. Correlation of the joint MT and CSEM inversion results with seismic and geologic data (a) sequential inversion results, (b)simultaneous inversion results, and (c) plot of a geologic section from Virgilio et al. (2009).



Figure 9. MT data fit of the joint MT and CSEM inversion resultfor (a) TM impedance (Ohm) and (b) TE impedance (Ohm) at0.05 Hz. The top plots show the real-part of the impedances whilethe bottom plots show the imaginary-part.

Figure 10. CSEM data fit of the joint MT and CSEM inversionresults for CSEM frequencies at (a) 0.125 Hz, (b) 0.255 Hz,and (c) 0.375 Hz. The receiver used is located around the cen-ter of the inversion domain. The top plots show the amplitudeof the electric fields (in logarithmic scale) while the bottomplots show the phase of the electric fields. The electric fieldsare in V=m.



initial model we use a model with a homogeneous isotropic

seabed as shown in Figure 11c and 11d. The homogeneous

seabed resistivity is 6.67 Ohm-m.

Figure 12a and 12b show the reconstructed horizontal and

vertical resistivity images using only MT data. The inversion

uses 11 iterations.

We observe that the horizontal resistivity distribution of the

layered formation is reconstructed pretty well. However, the MT

data inversion fails to reconstruct the vertical resistivity distribu-

tion of the layered formation. In addition to that, similar to the

isotropic test case, the MT data inversion could not detect the

presence of the reservoir in either the horizontal or vertical re-

sistivity images. These observations are consistent with the sen-

sitivity of the MT measurements.

Next, we show the CSEM data inversion results in Figure 12c

and 12d. These results were obtained after 41 iterations. We

observe that the reservoir is present in the vertical resistivity

image. However, in the horizontal resistivity image, the CSEM

data inversion can produce only a vague image of the first and

second layers of the seabed. This is consistent with the 1D stud-

ies done in Everett and Constable (1999), Constable and Weiss

(2006), Ramananjaona et al. (2008), and the 2D study in Abuba-

kar et al. (2010), which show that the CSEM inline electric field

data are mainly sensitive to the vertical resistivity. Furthermore,

we notice that the CSEM data inversion only properly recon-

structs the top of the reservoir. The thickness and bottom shape

of the reservoir clearly are overestimated.

Next, we invert MT and CSEM data simultaneously. The

inversion results after 14 iterations are shown in Figure 12e and

12f. We observe some improvements in the reconstructed shape

of the reservoir and the layered formation of the seabed. More-

over the estimate of the vertical resistivity value of the reservoir

Figure 11. The (a, b) true and (c, d) initial TI-anisotropic models:(a, c) horizontal resistivity, (b, d) vertical resistivity. The colorbars are in terms of the logarithm of the resistivity (Ohm-m).

Figure 12. The TI-anisotropic inversion results: (a) MT, horizon-tal resistivity; (b) MT, vertical resistivity; (c) CSEM, horizontalresistivity; (d) CSEM, vertical resistivity; (e) joint, horizontal re-sistivity; (f) joint, vertical resistivity.

Figure 13. The isotropic inversion result of the simultaneous jointinversion algorithm.



is more accurate. However, we still could not detect the reser-

voir in the horizontal resistivity image. These results again are

consistent with the physics of the measurements: MT data

mainly are sensitive to the horizontal resistivity, while the

CSEM inline electric field data mainly are sensitive to the verti-

cal resistivity.

MT data also are not sensitive to the presence of thin resistive

targets. Hence, we observe only the presence of the resistive

reservoir in the vertical resistivity image. Note also that the arti-

facts in the left and right regions of the images are caused by

the lack of sources or receivers on top of these regions. In Fig-

ure 13, we also show the joint inversion result when we use the

isotropic assumption. Although we obtain an object that resem-

bles the reservoir, its location and shape are not accurate. Fur-

thermore, we also observe a few artifacts appearing as layers.

This shows the importance of using a correct resistivity anisot-

ropy model for CSEM and MT data inversion.

CONCLUSIONS

We presented an inversion algorithm to jointly invert the CSEM

and MT data. An important issue for this simultaneous joint inver-

sion approach is how to assign the relative weights between CSEM

and MT data. In this work we proposed to minimize the product of

CSEM and MT cost functions instead of the traditional additive

one. In this multiplicative cost function there is no need to select

the relative weighting between these two data sets. This cost func-

tion adaptively will put CSEM and MT data on equal footing in

the inversion process. The effectiveness of this approach was dem-

onstrated using synthetic and field data sets. The approach can be

extended straightforwardly for 3D geometries.

ACKNOWLEDGMENTS

We thank V. Druskin of Schlumberger-Doll Research (SDR)

and L. Knizhnerman of the Center for Geophysical Expedition for

providing the 2D finite-difference forward programs. We thank

WesternGeco-Electromagnetic (WG-EM) for providing the

FORCE data set, and L. Masnaghetti and R. Luetkemeier from

WG-EM for their help in interpreting the inversion results and pro-

viding the seismic plot. We also acknowledge valuable comments

from C.G. Farquharson, M. Jegen-Kulcsar, and two anonymous

reviewers who significantly improved the quality and readability of

this manuscript.

REFERENCES

Abubakar, A., T. M. Habashy, V. L. Druskin, L. Knizhnerman, and D.Alumbaugh, 2008, 2.5D forward and inverse modeling for interpretinglow-frequency electromagnetic measurements: Geophysics, 73, no. 4,F165–F177, doi:10.1190/1.2937466.

Abubakar, A., T. M. Habashy, M. Li, and J. Liu, 2009, Inversion algo-rithms for large-scale geophysical electromagnetic measurements:Inverse Problems, 25, 123012, doi:10.1088/0266-5611/25/12/123012.

Abubakar, A., J. Liu, M. Li, T. M. Habashy, and K. MacLennan, 2010,Sensitivity study of multi-sources receivers CSEM data for TI-anisot-ropy medium using 2.5D forward and inversion algorithm: 72nd EAGEConference and Exhibition.

Abubakar, A., and P. M. van den Berg, 2000, Iterative reconstructions ofelectrical conductivity from multiexperiment low-frequency electro-magnetic data: Radio Science, 35, 1293–1306, doi:10.1029/2000RS002349.

Blystad, P., H. Brekke, R. B. Færseth, B. T. Larsen, J. Skogseid, and B.Tørudbakken, 1995, Structural elements of the Norwegian continentalshelf. Part II: The Norwegian Sea region. Norwegian Petroleum Direc-torate bulletin, 8.

Charbonnier, P., L. Blanc-Feraud, G. Aubert, and M. Barlaud, 1997,Deterministic edge-preserving regularization in computed imaging:IEEE Transactions on Image Processing, 6, 298–311.

Chave, A. D., S. C. Constable, and R. N. Edwards, 1991, Electrical explo-ration methods for the seafloor, in Nabighian, M. N., ed., Electromag-netic methods in applied geophysics, vol. 2: SEG, 931–966.

Colombo, D., and M. De Stefano, 2007, Geophysical modeling via simul-taneous joint inversion of seismic, gravity, and electromagnetic data:Application to prestack depth imaging: The Leading Edge, 26, 326–331, doi:10.1190/1.2715057.

Commer, M., and G. A. Newman, 2008a, New advances in three-dimen-sional controlled-source electromagnetic inversion: Geophysical JournalInternational, 172, 513–535, doi:10.1111/j.1365-246X.2007.03663.x.

——–, 2008b, Optimal conductivity reconstruction using three-dimen-sional joint and model-based inversion for controlled-source and mag-netotelluric data: 78th Annual International Meeting, SEG, ExpandedAbstracts, 609–613.

Constable, S., 2010, Ten years of marine CSEM for hydrocarbon explora-tion: Geophysics, 75, 7567–7581.

Constable, S. C., C. S. Cox, and A. D. Chave, 1986, Offshore electromag-netic surveying techniques: 56th Annual International Meeting, SEG,Expanded Abstracts, 81–82.

Constable, S., K. Key, and L. Lewis, 2009, Mapping offshore sedimentarystructure using electromagnetic methods and terrain effects in marinemagnetotelluric data: Geophysical Journal International, 176, 431–442,doi:10.1111/j.1365-246X.2008.03975.x.

Constable, S., A. Orange, G. M. Hoversten, and H. F. Morrison, 1998,Marine magnetotellurics for petroleum exploration Part 1: A sea-floor equipment system: Geophysics, 63, 816–825, doi:10.1190/1.1444393.

Constable, S., and C. J. Weiss, 2006, Mapping thin resistors and hydrocar-bons with marine EM methods: Insights from 1D modeling: Geophy-sics, 71, no. 2, G43–G51, doi:10.1190/1.2187748.

Davis, T. A., and I. S. Duff, 1997, An unsymmetric pattern multifrontalmethod for sparse LU factorization: SIAM Journal on Matrix Analysisand Applications, 18, 140–158, doi:10.1137/S0895479894246905.

Druskin, V. L., and L. Knizhnerman, 1994, Spectral approach to solvingthree-dimensional Maxwell’s diffusion equations in the time andfrequency domains: Radio Science, 29, 937–953, doi:10.1029/94RS00747.

Eidesmo, T., S. Ellingsrud, L. M. MacGregor, S. Constable, M. C. Sinha,S. Johansen, F. N. Kong, and H. Westerdahl, 2002, Sea bed logging(SBL), a new method for remote and direct identification of hydrocar-bon filled layers in deepwater areas: First Break, 20, 144–152.

Everett, M. E., and S. Constable, 1999, Electric dipole fields over an aniso-tropic seafloor: theory and application to the structure of 40 Ma PacificOcean lithosphere: Geophysical Journal International, 136, 41–56,doi:10.1046/j.1365-246X.1999.00725.x.

Farquharson, C. G., 2008, Constructing piecewise-constant models in mul-tidimensional minimum-structure inversions: Geophysics, 73, no. 1,K1–K9, doi:10.1190/1.2816650.

Farquharson, C. G., and D. W. Oldenburg, 1998, Non-linear inversion usinggeneral measures of data misfit and model structure: Geophysical JournalInternational, 134, 213–227, doi:10.1046/j.1365-246x.1998.00555.x.

Golub, G. H., and C. F. Van Loan, 1996, Matrix computations, 3rd edition:Johns Hopkins University Press.

Habashyx, T. M., and A. Abubakar, 2004, A general framework for con-straint minimization for the inversion of electromagnetic measurements:Progress In Electromagnetics Research, 46, 265–312, doi:10.2528/PIER03100702.

Hoversten, G. M., H. F. Morrison, and S. C. Constable, 1998, Marinemagnetotellurics for petroleum exploration, Part 2: Numerical analysisof subsalt resolution: Geophysics, 63, 826–840, doi:10.1190/1.1444394.

Ingerman, D., V. Druskin, and L. Knizhnerman, 2000, Optimal finitedifference grids and rational approximations of the square root I. Ellip-tic problems: Communications on Pure and Applied Mathematics, 53,1039–1066, doi:10.1002/1097-0312(200008)53:8<1039::AID-CPA4>3.0.CO;2-I.

Keller, J. B., 1964, A theorem on the conductivity of a composite me-dium: Journal of Mathematical Physics, 5, 548–549, doi:10.1063/1.1704146.

Kumar, A., L. Wan, and M. S. Zhdanov, 2007, Regularization analysis ofthree-dimensional magnetotelluric inversion: 77th International Meet-ing, SEG, Expanded Abstracts, 26, 482–486.

MacGregor, L., and M. Sinha, 2000, Use of marine controlled-source elec-tromagnetic sounding for sub-basalt exploration: Geophysical Prospec-ting, 48, 1091–1106, doi:10.1046/j.1365–2478.2000.00227.x.



Mackie, R., M. D. Watts, and W. Rodi, 2007, Joint 3D inversion of marineCSEM and MT data: 77th International Meeting, SEG, ExpandedAbstracts, 26, 574–578.

Plessix, R-E., and W. A. Mulder, 2008, Resistivity imaging withcontrolled-source electromagnetic data: Depth and data we-ighting: Inverse Problems, 24, 034012, doi:10.1088/0266-5611/24/3/034012.

Ramananjaona, C. J., D. L. Andreis, and L. M. MacGregor, 2008, Charac-terisation of anisotropic resistivity from Marine CSEM data: 70thEAGE Conference and Exhibition.

Rodi, W., and R. L. Mackie, 2001, Nonlinear conjugate gradients algo-rithm for 2-D magnetotelluric inversion: Geophysics, 66, 174–187,doi:10.1190/1.1444893.

Smirnov, M. Y., and L. B. Pedersen, 2009, Magnetotelluric measurementsacross the Sorgenfrei-Tornquist Zone in southern Sweden and Denmark:Geophysical Journal International, 176, 443–456, doi:10.1111/j.1365-246X.2008.03987.x.

Srnka, L. J., 1986, Methods and apparatus for offshore electromagneticsounding utilizing wavelength effects to determine optimum source anddetector positions: U.S. patent 4,617,518.

van den Berg, P. M., and A. Abubakar, 2001, Contrast source inversionmethod: State of art: Progress in Electromagnetics Research, 34, 189–218, doi:10.2528/PIER01061103.

van den Berg, P. M., A. Abubakar, and J. T. Fokkema, 2003, Multiplica-tive regularization for contrast profile inversion: Radio Science, 38,8022–8031, doi:10.1029/2001RS002555.

van den Berg, P. M., A. L. van Broekhoven, and A. Abubakar, 1999,Extended contrast source inversion: Inverse Problems, 15, 1325–1344,doi:10.1088/0266-5611/15/5/315.

Virgilio, M., L. Masnaghetti, M. Clementi, and M. D. Watts, 2009, Simul-taneous joint inversion of MMT and seismic data for sub-basalt explora-tion on the Atlantic margin, Norway: 71st EAGE Conference andExhibition, P258.

Vogel, C. R., and M. E. Oman, 1996, Iterative methods for total variationdenoising: SIAM Journal on Scientific Computing, 17, 227–238,doi:10.1137/0917016.

Zaslavsk, M., S. Davydycheva, V. Druskin, A. Abubakar, T. Habashy, andL. Knizherman, 2006, Finite-difference solution of the three-dimen-sional electromagnetic problem using divergence-free preconditioners:76th Annual International Meeting, SEG, 775–778.



Date post:	15-Nov-2023
Category:	Documents
Upload:	slb
View:	1 times
Download:	0 times

Joint MT and CSEM data inversion using a multiplicative cost function approach

Documents