Finite Difference Migration Imaging of ... - Semantic Scholar

JOURNAL OF COMPUTERS, VOL. 8, NO. 10, OCTOBER 2013

2553

Finite Difference Migration Imaging of Magnetotellurics Runlin Du China University of Petroleum (East China) Qingdao 266555 China Qingdao Institute of Marine Geology Qingdao 266071 China Email: [email protected]

Zhan Liu* China University of Petroleum (East China) Qingdao 266555 China Email: [email protected]

Abstract—we put forward a new migration imaging technique of Magnetotellurics (MT) data based on improved finite difference method, which increased the accuracy of difference equation and imaging resolution greatly. We also discussed the determination of background resistivity and reimaging. The processing results of theoretical model and case study indicated that this method was a more practical and effective for MT imaging. Finally the characteristics of finite difference migration imaging were summarized and the factors which can affect the migration imaging were analyzed. Index Terms—Magnetotellurics, Migration imaging, Finite difference

I. INTRODUCTION With the deepening of oil and gas exploration, the exploration blocks have become complex increasingly. MT sounding method with the restrictions of other geophysical exploration methods has become an increasingly important exploration method [1-3]. The range of MT method is much wider than seismic exploration and the probing depth also deeper than seismic and other methods. Therefore, in deep and complex lithology region, MT has certain advantages, especially for deep earth structure and the distribution of igneous rocks. It can also benefit to searching for a high impedance strata of deposited depression with oil and gas structures in the overthrust fault zone and volcanic rocks [4-5]. MT imaging techniques have been improved considerably in recent years on the basis of propagation similarity between the EM wave in the conductive media and the seismic wave in the elastic media. The PSI (Pseudo Seismic Interpretation) method can provide the average velocity of the EM wave depending on the similar process of reflective seismic method. The phase-shift migration method realized the downward continuation of the separated MT and used the U/D theory to image the subsurface conductivity boundaries by importing the advanced seismic wave equation theory such as phase-shift migration and PSPI [6-7]. Take the seismic reverse time concept as reference, Zhdanov put forward a new method of resistivity imaging based on © 2013 ACADEMY PUBLISHER doi:10.4304/jcp.8.10.2553-2557

frequency-domain electromagnetic migration [8-10]. Instead of transforming the diffusive field into wave field, this method transferred the principles of wave field analysis to interpretation of the EM field and the diffused time-reversed fields were called migrated fields. This method could be applied for determining both the position and resistivity imaging of the anomalous structures or interfaces [11-13]. The methods have been significantly improved in the recent years. However, there are still no fully developed approaches to image the resistivity property itself, which is the key problem in the inversion of EM data [14-15]. We present a new method of the resistivity imaging based on frequency-domain electromagnetic migration. II. MIGRATION IMAGING OF MT Zhdanov put forward the idea of electromagnetic field offset. The imaging problems of electromagnetic field offset and elastic wave seismic exploration offset is similar . It puts the analysis principle of seismic wave field method into the interpretation of the electromagnetic field under the conditions of the diffusion equation, instead of putting the electromagnetic wave meeting the conditions of the diffusion field into a wave field. The elastic wave field satisfies the general wave equation, while the electromagnetic field in our study accords with the diffusion equation. There is a similar mechanism for electromagnetic field and the elastic wave field. The imaging function named map can characterize the interface of different mediums [16], that is

map( x, z ) 

u( x, z, tu ) d ( x, z, td )

where u and d denote the up going and down going wave functions respectively , x, z indicate the horizontal and vertical coordinates respectively, and tu , td represent the travel time of up going and down going waves. The same phase of incident wave and reflected wave indicates that the time for the initial incidence is equivalent to the beginning of reflection. The electromagnetic field received on the earth surface is the superimposed result of

2554

JOURNAL OF COMPUTERS, VOL. 8, NO. 10, OCTOBER 2013

the up going and down going fields. The up going field shows a kind of back propagation that goes from the observation surface into the earth. In other words, it is a superimposed field of various kinds of source fields with different mediums under the ground. The migration imaging method can relocate the up going field to each source field successively, i.e., the interface of different mediums in a reverse time transmission[17-19]. A. TE polarization mode Suppose y indicates the landscape orientation of underground target in 2D (with variable property of electricity). When the variation of medium resistivity is not obvious, we can express the component y of electromagnetic field approximately as follows

E y ( x, z,  )  GEd ( x, z,  ) exp(ikn z ) GEu ( x, z,  ) exp(ikn z )

(1)

Combining equation (5) and (6) can get

where G Ed , G Eu are unknown coefficients relating to depth,

kn  i0 n ( x, z) ,  n ( x, z ) is back-

ground conductivity. Down going wavefield and up going wavefield of electric field are represented as

H yd ( x, z ,  )  GMd ( x, z,  ) exp(ikn z )

(2)

H yu ( x, z ,  )  GMu ( x, z ,  ) exp( ikn z ) (3) The up going and down going fields satisfy Helmholtz equation,

 2  2  u ,d   x 2 z 2  H y ( x, z,  )   2  kn ( x, z,  ) H yu ,d ( x, z,  )  0

kn  2GEd  3GEd  3GEd   (2 iz ) x 2 z z 3 x xz k k 2 G d (2i n  2iz n ) E x xz x k  2GEd (2ikn  2iz n ) x z 2 (6) kn kn2 GEd (4ikn  2iz 2 ) z z z kn kn d (2kn  2 z ) GE z z k 2 k GEd (2kn z  2i )( n2 GEd  n )0 z z z

(4)

Integrate equation (2) and (4)

kn GEd  2GEd  2GEd   (2 iz ) x 2 z 2 x x k G d (5) (2iz n  2ikn ) E z z k (2kn z  2i ) n GEd  0 z The above equation to z for partial derivative can get

kn  2GEd  3GEd  (2 ik  2 iz ) n x 2 z x x 2 G d [(2ikn ) 2  6(kn z  i )] E z 2 d k  GE (2iz n ) x xz k k k G d [(4kn z  2i ) n  n n 4 z 2 ] E x x z x 2 k k (4ikn2 z n  2kn n2 )GEd z z

(7)

The above equation ignores G Ed to z three derivative terms, and

k n to x, z second order derivative and the

first order derivative of the squared term. Repeating the process of equation (5) - (7), similar results are available for the up going wavefield of the electric field.

kn  2GEu  3GEu   (2 ik  2 iz ) n x 2 z x x 2 G u [(2ikn ) 2  6(kn z  i )] E z 2 u k  GE (2iz n ) x xz k k k G u [(4kn z  2i ) n  n n 4 z 2 ] E x x z x 2 k k (4ikn2 z n  2kn n2 )GEu z z

(8)

Based on function (7) and (8), there is the following equation

© 2013 ACADEMY PUBLISHER

JOURNAL OF COMPUTERS, VOL. 8, NO. 10, OCTOBER 2013

kn  2GEd ,u  3GEd ,u   (2 ik  2 iz ) n x 2 z x x 2 G d ,u [(2ikn ) 2  6(kn z  i )] E z 2 d ,u k  GE  (2iz n ) x xz k k k G d ,u [(4kn z  2i ) n  n n 4 z 2 ] E x x z x 2 k k (4ikn2 z n  2kn n2 )GEd ,u z z

 3GEpu, M

(9)

Where ’+’ indicates the down going wavefield and ’-’ indicates up going wavefield. B. TM (Transverse- Electric) polarization mode As for the TM mode, repeating the process of equation (5) - (7) can get similar results

kn  2GMd ,u  3GMd ,u   (2 ik  2 iz ) n x 2 z x x 2 G d ,u [(2ikn ) 2  6(kn z  i )] M z 2 d ,u k  GM  (2iz n ) x xz k k k G d ,u [(4kn z  2i ) n  n n 4 z 2 ] M x x z x 2 k k (4ikn2 z n  2kn n2 )GMd ,u z z

(10)

H ypu ( x, z,  )  GMp ( x, z, ) exp(kn z ) (12) *

*

and magnetic field

of migration up going fields satisfy

E u ( x, z ,  ) E d ( x, z ,  )

electrical

and

   pu*  x 2  z 2  E y ( x, z ,  )  

 Ea ( x, z,  )  Ea ( x, z,  )  ( x, z,  ) n  Ma ( x, z,  )  Ma  Ma ( x, z,  )

(13)

kn2 ( x, z ,  ) E ypu ( x, z ,  )  0 *

 2  2  pu*  x 2  z 2  H y ( x, z ,  )  

(17)

reflection (18)

(19)

 En,Ma ( x, z,  ) 

1 N n   E ,Ma ( x, z,  j ) N j 1

(20)

Realize the resistivity imaging based on the normalized

 En,Ma ( x, z,  ) and n E，Ma

the

( x, z,  ) . For

inhomogeneous background resistivity n ( x, z ) , we can

m ( x, z) as

calculate the migrated apparent resistivity 2

1   En, Ma ( x, z ) En, Ma ( x, z )    n ( x, z ) 1   En, Ma ( x, z ) En, Ma ( x, z ) 

 m ( x, z )   (14)

 kn2 ( x, z ,  ) H ypu ( x, z ,  )  0 *

For equation (11) - (14), doing the derivation and simplification similar to equation (5) - (10) can obtain © 2013 ACADEMY PUBLISHER

(16)

The superimposed equation is

apparent reflectivity function 

2

magnetic

n  Ea ( x, z,  ) 

migrated reflectivity function

Helmholtz equation 2

 Ma ( x, z,  ) 

Normalizing coefficients

*

pu*

E u ( x, z ,  ) E d ( x, z ,  )

G u ( x, z ,  )  Md exp(2ikm z L ) GM ( x, z ,  )

Eypu ( x, z,  )  GEp ( x, z, ) exp(kn z) (11)

component H y

D. Joint imaging of TE and TM polarization modes For the TE and TM polarization modes, the apparent electric reflection coefficient is the ratio of up going wave and down going wave according to the U/D imaging principle

G u ( x, z ,  )  Ed exp(2ikm z L ) GE ( x, z ,  )

Go on with of the offset field using the following equations

pu

2 pu kn  GE , M ) x 2 z x x 2 GEpu, M 2 [(2kn )  6(kn z  1)] z 2 pu k  GE , M (15) (2 z n ) x xz G pu k k k [(4kn z  2) n  n n 4 z 2 ] E , M x x z x 2 k k (4kn2 z n  2kn n2 )GEpu, M z z

 (2kn  2 z

 Ea ( x, z,  ) 

C. Migration electromagnetic field

The electric field component E y

2555

(21) Furthermore, another improvement also has been made by joining the TE and TM polarized mode together using

2556

JOURNAL OF COMPUTERS, VOL. 8, NO. 10, OCTOBER 2013

weighted coefficient ( 0  1、2

 1, 1  2  1 ) to

achieve a joint imaging result 2

n n n n 1  (1  Ea ( x, z )  2  Ma ( x, z ))(1 Ea ( x, z )  2  Ma ( x, z ))  m ( x, z )     (1 En ( x, z )  2  Mn ( x, z )) n n n n 1  (1  Ea ( x, z )  2  Ma ( x, z ))(1 Ea ( x, z )  2  Ma ( x, z )) 

Two modes joint imaging gives prominence to the consistency of the same geoelectric model and reinforces the single information. III. MODEL TEST AND REAL OBSERVED MT DATA IMAGING Devise a geoelectric model as shown in Figure 1: the resistivity of first layer is 100  and thickness is 1000m. Oblique bottom is at 1400m. The resistivity of second layer is 10  and thickness is 1000m.The resistivity of the third layer is 100  . Figure 2 gives the result of TE, TM joint imaging. We can see that the migrated result is better and clearly reflects the syncline structure.

Figure 1. Schematic Diagram of Syncline

(22) From the second figure with the third one, the result of two modes of TE, TM joint imaging is more clear and closer to the real geological model, indicating the method of two modes of TE, TM joint imaging better than continuous rapid inversion. However, if using conventional inversion to obtain the preliminary estimate of the underground structure firstly, the resistivity values are similar to the actual resistivity distribution, with the result as the background conductivity offset imaging effect will be better. So the two methods complement each other in practice in order to serve better for MT interpretation. IV. CONCLUSION This paper put forward a new multi-parameter MT migration imaging technique using the improved finite difference method, which had greatly increased the imaging accuracy and resolution. It is a multi-parameter, multi-mode imaging method and can simultaneously obtain information on the interface and resistivity. In comparison with conventional inversion, two modes of TE, TM joint imaging can get better vertical resolution with less false information. ACKNOWLEDGMENT

Figure 2. The Output of the Two Modes of TE, TM Joint Imaging

We use the method to deal with the MT data of NE region of the Junggar Basin, China. The gobi and desert are the main terrain in this region.

This study was supported by The joint inversion of CSEM and seismic(CX2013001) and 1:250000 Qingdao site marine areas Geological Survey "(pilot) project (GZH200900501) and Innovative Scientific Research Plan of China University of Petroleum (12CX06001A) and High-precision gravity, magnetic and electrochemical exploration of oil and gas information detection and comprehensive evaluation (CN-2012-JS-88). REFERENCES [1]

[2]

[3]

[4]

[5]

a: Geological profile interpretation chart. b:The result of continuous rapid inversion. c:The result of Two Modes of TE,TM joint imaging. Figure 3 © 2013 ACADEMY PUBLISHER

[6]

M.S.Zhdanov, P.Traynin, J.R.Booker, “Underground imaging by frequency domain electromagnetic migration,”Geophysics, vol. 61, no. 3, pp. 666-682, 1996. J.T.Smithand, J.R.Booker, “Rapid Inversion of Two and Three-Dimensional Magnetotelluric. Data,” JGR, vol. 96, no. 3, pp. 3902-3922, 1991. H.DeGroot, “Constable S. Occam inversion to generate smooth two-dimensiona1 models from magnetotelluric data,” Geophysics, vol. 55, no. 12, pp. 1613-1624, 1990. H.Mauriello, “Principles of Probability Tomography for Natural Source electromagnetic Induction Fields,” Geophyslcs, vol. 64, no. 5, pp. 1403-1417, 1999. W.Rodl, R.L.Mackie, “Nonlinear conjugate gradients algorithm for 2-D Magnetotelluric inversion,” Geophysics , vol. 66, no. 1, pp. 55-67, 2001. L.Lee, G.A.McMechan, C.L.Aiken, “Phase field imaging: The electromagnetic equivalent of seismic migration,” Geophysics, vol. 52, no. 5, pp. 687-693, 1987.

JOURNAL OF COMPUTERS, VOL. 8, NO. 10, OCTOBER 2013

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

L.Lee, J. R.Booker, “Underground imaging by electromagnetic migration,” Mtg. ExPI. GeoPhys. Expanded Abstractst ,1993, pp. 355-357. M.S.Zhdanov, P.Traynin, “Portniaguine.,Migration and analytic continuation in geoelectric imaging,” 64th SEG EMI. 1994, pp. 1325-1342. M.S.Zhdanov, “Joint iterative migration of electric and magnetic-field data,” SEG Expanded Abstracts 26, 2007, pp. 259-268. M.S.Zhdanov, “Electromagnetic migration of marine CSEM data in areas with rough bathymetry,” SEG Expanded Abstracts 28, 2009, pp. 149-158. M.S.Zhdanov, S.Fang, “Electromagnetic inversion using quasi-linear approximation,” Geophysics, vol. 65, no. 5, pp. 1501-1513, 2000. T.Sirlpunvaraporn,G.Weerachai, “An efficient datasubspace inversion method for 2-D magnetotelluric data,” Geophysics, vol. 65, no. 3, pp. 791-803, 2000. A.Kumar, “Regularization analysis of three- dimensional magnetotelluric inversion,” SEG Expanded Abstracts 26, 2007, pp.439-443. F.Toshiko, “Two-dimensional time domain electromagnetic migration using integral transformation,” SEG Expanded Abstracts 26, 2007, pp. 1724-1739. W.Le, “Frequency domain 3-D electromagnetic migration and imaging of sea-bottom geoelectrical structures,” SEG Expanded Abstracts 26, 2007, pp. 1241-1258.

© 2013 ACADEMY PUBLISHER

2557

[16] W.Le, “Rapid seabed imaging by frequency domain electromagnetic migration,” SEG Expanded Abstracts 24, 2005, pp. 321-337. [17] Yanqin Zhang, Xibei Yag, “Intuitionistic fuzzy dominance-based rough set approach: model and attribute reductions.” Journal of Software, 2012, 7(3), pp.551-563. [18] Qianjin Wei, Gu Tianlog, “Symbolic Representation for Rough Set Attribute Reduction Using Ordered Binary Decision Diagrams.” Journal of Software, 2011, 6(6), pp.977-984. [19] Kaiji Liao, Huihui Xiong, Donghai Ye, Chunyang Zhang. “A Method of Emergency Management based on Knowledge Element Theory.” Journal of Software, 2012, 7(1) pp.41-48.

Runlin Du was born in Shandong province of China on 29 February, 1987. He received the B.S and M.S degree in China University of Petroleum. His research interests focus on gravity and magnetic data processing and interpretation. Zhan Liu was born in Sichuan province of China on 27 April, 1957. He received Ph.D. degree in December 1999. His major field of study is joint exploration.