adg vol5 n02 Troyn 339_349.pdf ANNALS OF GEOPHYSICS, VOL. 45, N. 2, April 2002 339 Mathematical modeling to reconstruct elastic and geoelectrical parameters Vladimir N. Troyan and Yurii V. Kiselev Institute of Physics, St. Petersburg State University, St. Petersburg, Russia Abstract The monitoring of the underground medium requires estimation of the accuracy of the methods used. Numerical simulation of the solution of 2D inverse problem on the reconstruction of seismic and electrical parameters of local (comparable in size with the wavelength) inhomogeneities by the diffraction tomography method based upon the first order Born approximation is considered. The direct problems for the Lame and Maxwell equations are solved by the finite difference method that allows us to take correctly into account the diffraction phenomenon produced by the target inhomogeneities with simple and complex geometry. For reconstruction of the local inhomogeneities the algebraic methods and the optimizing procedures are used. The investigation includes a parametric representation of inhomogeneities by the simple and complex functions. The results of estimation of the accuracy of the reconstruction of elastic inhomogeneities and inhomogeneities of electrical conductivity by the diffraction tomography method are represented. 1. Introduction Diffraction tomography is an imaging tech- nique that makes use of a large volume of input data (recorded traces) to produce the image of underground medium parameters with high spatial resolution. In contrast to ray tomography (travel time tomography), for which the re- solution is connected with the Fresnel zone and the large number of the source-receiver pairs is required, diffraction tomography (Devaney, 1984; Devaney and Zhang, 1991; Zhou et al., 1993; Ryzhikov and Troyan, 1994; Alumbaugh and Morrison, 1995; Kiselev and Troyan, 1997) provides information on the medium parameters with subwavelength resolution. In our study (a complete review of dev- elopment of the diffraction tomography is in- troduced in (Devaney and Zhang, 1991)) most attention is given to an estimation of the accuracy of multiparametric reconstruction of elastic parameters and reconstruction of an electrical conductivity with the use of sounding by elastic wave and electromagnetic wave correspondingly. As a tool for our investigation, we use numerical simulation in 2D space domain. The direct problem is solved by the finite difference method that allows us correctly to take into account the diffraction phenomenon produced by the target inhomogeneities. The linearized inverse problem is solved by the diffraction tomography method with the use of the first order Born approximation (Keller, 1969). Mailing address: Dr. Vladimir N. Troyan, Univer- sitetskaya nab. 7/9, St. Petersburg, 199034 Russia; e-mail: troyan@hq.pu.ru Key words diffraction tomography – mathematical modeling – seismic tomography – electromagnetic tomography 340 Vladimir N. Troyan and Yurii V. Kiselev The applicability of the first order Born approximation in diffraction tomography is shown (Slaney et al., 1984) by numerical sim- ulation conducted on a single cylinder using analytical expressions for the exact scattering field. Beylkin and Burridge (1990) describe the multiparametric inversion in the time domain in the cases of acoustic and elasticity on the basis of the generalized Radon transform. This approach requires a large number of the source- receiver pairs. We also implement multi- parametric reconstruction in the time domain, but our numerical simulation is realized with not more than three sources and three receivers (nine source-receiver pairs). For the multi- parametric reconstruction in the elastic case (reconstruction of the Lame parameters, mass density and as a corollary – shear and lon- gitudinal velocities) the optimizing procedures are used. The multiparametric reconstruction makes use of amplitude information connected to the scattering characteristic (Wu and Aki, 1985; Beylkin and Burridge, 1990) of the elementary disturbances of the reconstructed parameters. Under tomography experiment, these scattering characteristics should be taken into account by the use of the relevant ob- servation schemes. Scattering by the elementary disturbances of the parameters can be described with the use of tomography functionals (Ryzhikov and Troyan, 1994; Troyan and Ryzhikov, 1994), which in a form of the ray series are represented for elastic and electro- magnetic cases. 2. Basic equations and algorithms for elastic case The numerical simulation to reconstruct the parameters of local inhomogeneities with a smooth change of the elastic parameters , µ and mass density is carried out for two-dimensional model of the elastic medium, containing the free surface and plane-parallel welded interfaces. The source f f (x, t), located in the point (x = x S , z = z S ) of the Cartesian system of coordinates (x, y ,z; e 1 , e 2 , e 3 ), produces the wave field u u (x, z, t) u (x, t) which satisfies the equation (2.1) Boundary condition at the free surface (z = 0) is and at welded interfaces (z = z i ) we assume that (2.3) The field u will be produced by the sources f , f 1, f 3 (2.4) with the time dependences f (t) and (t), which are located at the source point x = (x S , z S ) and at the observation point x = (x r , z r ). The velocities of longitudinal ( p ) and shear ( s ) waves, expressed through the quantities , µ, , read as (2.5) We introduce the differences = rf , µ = µ µ rf and = rf of the values , µ, , for the unknown medium, which are connected with µ µ2u u u+ + × + ( ) . 0 0tz z == f f x z t x x z zs s s s( ) = ( ) ( )ˆ ˆ , , ff t( )e3 x z t x x z z tr r r r( ) = ( ) ( ) ( )1 1 1ˆ ˆ , ,f f e x z t x x z z tr r r r( ) = ( ) ( ) ( )3 3 3ˆ ˆ , ,f f e µ µ p s= + = 2 , . (2.2) µ µ 2 2 u f u u t = + +( ) + +ˆ u u z z z zi i= = +=0 0 , t tz z z z z zi i= = +=0 0 . 341 Mathematical modeling to reconstruct elastic and geoelectrical parameters the wave field u(x, t) (2.6) and the values λrf (x), µrf (x), ρrf (x) for the known, reference (rf) medium for which the wave field is u rf (x,t) (2.7) Assuming δλ, δµ, and δρ are small we can write (2.8) where δu = u − urf is the difference field. The right hand side of (2.8) (2.9) can be considered as a source of this field. We shall represent the components of the difference field δu i from (2.8) at the observation point of x = x r , as following: (2.10) where S is the region of reconstruction; and , are the solutions of the eqs. (with sources from (2.4)) (2.11) and (2.12) respectively. After introducing the tomography functionals (Troyan and Ryzhikov, 1994) (2.13) the components of the difference field δu i (2.10) can be written down as (2.14) Using the linear relations between δλ (x), δρ (x) and δµ (x) δλ (x) = cλδµ (x), δρ (x) = cρδµ (x), (2.15) (cλ = const, cρ = const) λ µ µL L= − ≡ +( ) ∇∇ ⋅ + +ˆu f, u u u∆ λ µ µL ≡ +( )∇∇ ⋅ + +rf rf rf rf rf rf rfu u u∆ u urf rf≈ −δ δL L u u erf rf rf= ∇ × ∇ ×( ) + ∇ ⋅ ∇   + = ∑δ δµ δµ 1 3 2L j j j , u δu ti s rx x, ,( ) = tru u x, x ,˜ ˜≡ ( )1 1 tru u x, x ,˜ ˜≡ ( )3 3 tsu u x, x ,≡ ( ) L u f̂= −rf L i iu f˜ ˆ= −rf p ti r s ρ x x x, , ,( ) = t di r sτ ∂ ∂τ τ τu x x u x x˜ , , , ,= − −( ) ⋅ ( ) ∞ ∫ 2 20 p ti s r λ x x x, , ,( ) = t t di r sτ τ τ u x x u x x˜ , , , ,= − ∇ ⋅ −( )∇ ⋅ −( ) ∞ ∫0 p ti s r µ x x x, , ,( ) = = + ∇ × ˜̃ , , , ,u x x u x xi r st −( ) ⋅ ∇ × ( ) −     ∞ ∫ τ τ0 δ δλ λu t p ti s r i s r S x x x x x x, , , , ,( ) = ( ) ( ) +[∫ δµ δρµ ρp t p t di s r i s rx x x x x x x x x, , , , , ,+ ( ) ( ) + ( ) ( )] . (i = 1, 3) L = − ˆrf rfu f , ˜ , , , , , x x x xij r j j st d) − ∇ −( ) ⋅ ∇ ( )     = ∑ τ τ τ2 1 3 u u u u uu u t + ∇ ∇ ⋅ + ∇ × + ∇ ⋅ ∇( ) −λ µ µ ρ ∂ ∂ 2 2 2 t rf rf rf rf rf rf rf u u u u + ∇ ∇ ⋅ + ∇ × + ∇ ⋅ ∇( ) −λ µ µ ρ ∂ ∂ rf 2 2 2 . δλ δρ ∂ ∂t u u rf rf+ ∇ ∇ ⋅( ) − 2 2 τ δ τ τt L d di s r ru x x u x x x˜ , , , ,) = −( ) ⋅ ( ) ∞ ∫∫ 0 342 Vladimir N. Troyan and Yurii V. Kiselev the eq. (2.14) can be rewritten as (2.16) After the digitization of the eq. (2.16), the system of equations for determination of δµ (vector dµ ), cλ and cρ can be written as P(cλ ,cρ )dµ = du (2.17) where d u are the samples of the scattered field. The final version of this system after introducing regularizing terms reads as (2.18) where α1, α2, α3 are the regularizing coefficients; matrices B x and B z are the finite difference images of second partial derivatives with respect to x and z correspondingly; C and D are penalty matrices for non-zero values of δµ at boundary and near boundary points of the reconstructed region S. We find δµ, cλ and cρ by using an iterative procedure. At the first step the system of linear eqs. (2.18) is solved with some initial values cλ (0) and cρ (0). By minimizing the sum of squared differences of the left-hand side and the right-hand side of (2.17), we find c λ (1) and cρ (1), which are the corrected values of cλ (0) and cρ (0). At the second step, the values c λ (1) and cρ (1) are used for solution of the system of linear eqs. (2.18). The convergence of this procedure is based on the distinctions of the scattering diagrams (Wu and Aki, 1985) created by the elementary disturbances of λ, µ, ρ. Similar distinctions can be studied, for example, using formulas (2.21) given below. 2.1. Ray representation of the tomography functionals The components of the wave field δu i (x S ,x r ,t) (i = 1, 3) scattered by local inhomogeneity (δλ, δµ, δρ) for incident wave field u (x, x S , t) are given by relations (2.14), (2.13). Using the ray method the fields u (x, x S ,t), u~ i (x,x r ,t) (2.11)-(2.13) can be represented as following: (2.19) where τ~ q ≡ τ~ q (x,x r ) (τ q ≡ τ q (x,x S )) is the time of propagation of a wave from point x r , (x S ) to a point x. Substituting of the relations (2.19) in (2.13) we write down the quantities in the form of the ray series (2.20) thus the coefficients read as (2.21) x x z zc c P c c B B B B, , ' '( ) ( ) + +( ) +[ λ ρ λ ρ αP' 1 uC C D D P c c' ' ' ,+ + ] = ( )µ λ ρα α d d2 3 f t f tj j( ) ,= ′ ( )+1 ˜ ( ) ˜f t f tj j= ′ ( )+1 p p pi i i, , ρ λ µ +( ) = ′ ( )j jt t1φ φ , ∫ ∞ ( ) = −( ) ( )t d dt f t f d0 2 2 00 0φ τ τ τ˜ p p piqq j iqq j iqq j′ ′ ′ λ µ ρ, , p p p iqq iq q iqq iq q q q iqq q iq q q ′ ′ ′ ′ ′ ′ ′ ′ = − ⋅( ) = − ⋅ ∇( ) ⋅ ∇( ) = ∇ ×[ ]⋅ ∇ ×[ ] − ρ λ µ τ τ τ τ ˜ ˜ ˜ ˜ ˜ 0 0 0 0 0 0 0 0 0 A A A A A A iq q′− ⋅( )˜2 0 0A A ∇∇ ⋅ ∇( )′˜ .τ τq q j =( )0 p p ti iqq j j q q q p sq p sj ρ ρ φ τ τ= − −( )′ ′ ′ === ∞ ∑∑∑ ˜ ,,0 p p ti iqq j j q q q p sq p sj λ λ φ τ τ= − −( )′ ′ ′=== ∞ ∑∑∑ ˜ ,,0 p p ti iqq j j q q q p sq p sj µ µ φ τ τ= − −( )′ ′ ′=== ∞ ∑∑∑ ˜ ,,0 δ λ λ u t c p ti s r i s rx x x x x, , , , ,( ) ≈ ( ) + S [∫ µpi sx x, ,+ xx x x x x xr i s rt c p t d, , , ,( ) + ( )] ( )ρ ρ δµ . ˜ , , ˜ , , ˜ ˜ , , , , , , u x x A x x u x x A x x i r iqj q p sj r j q s qj q p sj s j q t t f t t t f t ( ) = ( ) −( ) ( ) = ( ) −( ) == ∞ == ∞ ∑∑ ∑∑ 0 0 τ τ 343 Mathematical modeling to reconstruct elastic and geoelectrical parameters 3. Basic equations and algorithms for electromagnetic case Numerical simulation to reconstruct the local inhomogeneities of electrical conductivity σ = σ (x), located in the uniform space (electrical conductivity σ = const, electrical permittivity ε′ = εε0 = const, magnetic permittivity µ′ = µµ0 = const) is implemented for the 2D problem. Electrical (E = E(x,t)) and magnetic (H = H(x, t)) fields excited by a current density j ex = j ex (x, t) satisfy to the Maxwell equations (3.1) Electrical field E = E(x,t) in the medium, containing the local inhomogeneity (σ = σ (x), ε′ = ε′(x), µ′ = µ′(x)) (1) is given by a solution of the equation (3.2) The reference medium (rf ) is supposed to be known (σrf , εrf , µrf ) and electrical field Erf satisfies to the equation (3.3) As in the case of scattering by elastic inhomo- geneities discussed earlier, we assume that magnitudes of the values δσ = σ − σ rf , δε = ε′ − ε rf and δµ = µ′ − µ rf make it possible to write an approximate equality (3.4) where δE = E − Erf is the difference field. Thus, the value (3.5) can be considered as a source of this field. The components of the difference field δE i is possible to write down as (3.6) that coinciding to within notations with (2.10). Wave fields E(x, x s , τ ) and E ~ i (x, x r , t − τ ) satisfy, respectively, the following equations: (3.7) where in a two-dimensional case the sources from (2.4) are used. Introducing the tomography functionals (3.8) ∇ × = + + ∇ ⋅ = = ′ t ex ∂ ∂ ρ µ, ,H j j D D B H. = −L t ex ∂ ∂ E j = ′ ∇ × ∇ × + ′ + −L t tµ ε ∂ ∂ σ ∂ ∂ 1 2 2 E E E E − ′ ∇ ′ × ∇ ×[ ] µ µ . 1 2 E L t exrf rfE j= − ∂ ∂ L t rf= ∇ × ∇ × + + µ ε ∂ ∂ rf rf rf rf rfE E E1 2 2 t + − ∇ × ∇ ×[ ]σ ∂ ∂ µ µ 1 2 rf rf rf rf rf E E . L L≈ −δ δrf rfE E δ δε ∂ ∂ δσ ∂ ∂ δµ µ L t t E E E Erf rf rf rf rf= + − ∇ × ∇ × 2 2 2 ( ) δ ti s r( ) =E x x, , L = − ˆrf E f L i i= −˜ ˆrf E f ( = 1, 3)i (1) The main relations considered below will be written down in general case, i.e. we will assume an existance of anomalies of electrical permittivity and magnetic permittivity. ∇ × = − ∇ ⋅ = = ′ =E B B D E j E ∂ ∂ ε σ t , , ,0 τ δ τ τE x x E x x xi S r st L d d ˜ , , , ,= −( ) ⋅ ( ) ∞ ∫∫ 0 p t t d p t t d p t i r s i r s i r s i r s i r s i ε σ µ τ ∂ ∂τ τ τ τ ∂ ∂τ τ τ x x x E x x E x x x x x E x x E x x x x x E , , , ˜ , , , , , , , ˜ , , , , , , , ˜ ( ) = −( ) ⋅ ( ) ( ) = −( ) ⋅ ( ) ( ) = − ∇ × ∞ ∞ ∫ ∫ 0 2 2 0 00 ∞ ∫ −( ) ⋅ ∇ × ( )x x E x x, , , ,r st dτ τ τ and 344 Vladimir N. Troyan and Yurii V. Kiselev the components of the difference field can be written as (3.9) We will consider the numerical experiments on reconstruction of the anomalies of electrical conductivity ( = 0, µ = 0). In this case, after the digitization, the integral eq. (3.9) can write as the system of linear equations Pd = d E (3.10) with respect to vector d , which is sought for the value (x), where d E is the time samples of the components of the wave field scattered by inhomogeneity. The final version of these equations, after introducing the regularization terms, is coincident with the system of linear equations similar to (2.18). 3.1 Ray representation of the tomography functionals In Section 2.1 the expressions for the tomography functionals (2.13) were written out in the case of the ray description of the wave fields u and u~ i , which are included into the integrands of the right-hand sides of (2.13). By assuming that the value of electrical conductivity for the reference medium is equal to zero, we can represent the wave fields E and E ~ i from (3.8) as (3.11) and can write the tomography functionals (3.8) in a form of the ray series (3.12) In a case of n = 0 the values , and are re- presented by the amplitude factors of the zero approximation A 0i , A and eikonals ~, (from (3.11)) (3.13) The values (from (3.13)) describe the space characteristics of the scattered fields, which are produced by the elementary inhomoge- neities ( , , µ) at the far-field region. The directivity diagrams in the cases of scattering by perturbation of electrical conductivity and electrical permittivity coincide. The wave field scattered by perturbation of electrical permittivity contains higher frequencies than the wave field scattered by perturbation of electrical con- ductivity. It should be noted that formula (3.13) for the values and coincides with the formula for scattering the elastic wave field by the elementary perturbation of mass density (2). f ti ni n n ˜ ˜ ˜ ˜ ,= ( ) = 0 E A f tn n n ,= ( ) = 0 E A f t f tn n ,( ) = ( )+1 ˜ ˜f t f tn n( ) = ( )+1 0 2 2 0 0 0t d dt f t f d( ) = ( ) ( )˜ , 1t tn n( ) = ( )+ . p i0 p pi i i0 0 0 0= = ( )à A p i0 µ, ,0 0p pi i p i0 p i0 p i0 p i0 µ p i i0 0 0 µ = ×[ ] ×[ ]( )˜ ˜ .A A (2) Graphic representation of directivity diagrams for the elastic and electromagnetic cases can be found in Wu and Aki (1985) and Saintenoy and Tarantola (2001). ~ p p ti ni n n= ( ) =0 ˜ p p ti ni n n= ( ) = + 0 1 ˜ p p ti ni n n µ µ= ( ) =0 ˜ E t Ei s r ix x, ,( ) = µ µµp p p di i i S = + +( )rf .2 x 345 Mathematical modeling to reconstruct elastic and geoelectrical parameters 4. Numerical simulation 4.1. Elastic case As the result of numerical simulation, we obtain the values , µ and . With formula from (2.5) we get value of p . The observation scheme and the location of inhomogeneity are represented in fig. 1. The sources and observation points (9 pairs) are located at the free surface. The results of reconstruction of p are rep- resented in figs. 2a-d and 3a-d. Figure 2a and fig. 3a show the models of two inhomogeneities with 20% contrast relatively to the reference medium. These inhomogeneities are located inside the layer with intermediate velocity. Figure 2b,c and fig. 3b show the results of reconstruction obtained by the solution of the system of eqs. (2.18) with the different values of the regularizing coefficients. The relative error of the recon- struction in these cases is 20-25% (with the use of two components of the wave field). For the cases of more contrast inhomogeneities (20-40%) the accuracy of reconstruction can be ~ 50%. Figure 2d and fig. 3c,d show the results of reconstruction using one (fig. 2d and fig. 3d) or three (fig. 3c) parametric functions. In these cases, we solve the system of the eqs. (2.17) by minimization of the sum of squared differences between the left-hand side and the right hand- side of (2.17). The reconstruction with the use of just one simple parametric function is very stable. From numerical simulation, we conclude that the realization of the diffraction tomography method offered in Ryzhikov and Troyan (1994), under appropriate observation scheme, allows satisfactory accuracy for velocity parameter reconstruction in the case of not very contrast inhomogeneities of size of ~ p at small number of source-receiver pairs. 4.2. Electromagnetic case The model of inhomogeneity of electrical conductivity together with the results of re- construction are represented in fig. 4a-c. The inhomogeneity, comparable in size with the wavelength of the sounding signal, is located inside the uniform space (reference medium) with parameters: = 10, µ = 1, = 0. The maximum . Fig. 1. Model of medium and observation scheme. S is the region for reconstruction; location of inhomogeneity - black color; 1, 2, 3 are the source and the observation points locations; s ; s / p = 1/ .3 346 Vladimir N. Troyan and Yurii V. Kiselev a b c d Fig. 2a-d. Reconstruction of p for symmetric inhomogeneity. a) The model; b,d) results of reconstruction; b,c) solution of the system (2.18) with 3 = 0 and 3 0 respectively; d) reconstruction in case of representation of reconstructed inhomogeneity by simple parametric function. value of electrical conductivity is 10 3 S/m. An apparent frequency of the sounding signal is 5 × 106 Hz. We use the observation scheme which is similar to the observation scheme represented in fig. 1. Distance between the observation line and the center of inhomogeneity is 10 ( = 20 m). The system of the linear eqs. (3.10) is solved under introducing the minimum magnitudes of 347 Mathematical modeling to reconstruct elastic and geoelectrical parameters a b c d Fig. 3a-d. Reconstruction of p for asymmetric inhomogeneity. a) The model; b,d) results of reconstruction; b) solution of the system of eq. (2.18) with 3 = 0; c,d) representation of reconstructed inhomogeneity by three and one simple parametric function correspondingly. the regularizing coefficients 1 , 2 , 3 (2.18). The errors of reconstruction in considered numerical example are 50% (fig. 4b) and 30% (fig. 4c). The accuracy of reconstruction of inhomogeneities of electrical conductivity has a stronger de- pendence on the values of the regularizing coefficients in comparison with the elastic case. This difference can be explained by greater 348 Vladimir N. Troyan and Yurii V. Kiselev Fig. 4a-c. Reconstruction of electrical conductivity. a) The model; b,c) results of reconstruction; in the case (b) the regularizing coefficients are ten times greater than in the case (c). a b c 349 Mathematical modeling to reconstruct elastic and geoelectrical parameters distortion of a shape of the scattered signal in the case of reconstruction of electrical con- ductivity. Acknowledgements This investigation was supported by INTAS- 99-1102, Russian Foundation of Basic Research (Grant 02-05-65081), Ministry of Education (Grant E00-9.0-82) and «Intergeofizika». REFERENCES ALUMBAUGH, D.L. and H.F. MORRISON (1995): Theoretical and practical consideration for crosswell elec- tromagnetic tomography assuming a cylindrical geometry, Geophysics, 60, 846-870. BEYLKIN, G. and R. BURRIDGE (1990): Linearized inverse scattering problems in acoustics and elasticity, Wave Motion, 12 (1), 15-52. DEVANEY, A.J. (1984): Geophysical diffraction tomography, IEEE Trans. Geosci. Remote Sensing, GE-22 (1), 3-13. DEVANEY, A.J. and D. ZHANG (1991): Geophysical diffraction tomography in a layered background, Wave Motion, 14 (3), 243-265. KELLER, J.B. (1969): Accuracy and validity of the Born and Rytov approximations, J. Opt. Soc. Am., 59, 1003-1004. KISELEV, YU.V. and V.N. TROYAN (1997): Estimation of elastic parameters in diffraction tomography, in Proceedings of the International Seminar «Day on Diffraction-97», St. Petersburg, Russia, 3-5 June 1997, 135-139. RYZHIKOV, G.A. and V.N. TROYAN (1994): Tomography and Inverse Sounding Problems (St. Petersburg University, St. Petersburg), pp. 220. SAINTENOY, A.C. and A. TARANTOLA (2001): Ground- penetrating radar: analysis of point diffractors for modeling and inversion, Geophysics, 66, 540-550. SLANEY, M., A.C. KAK and L.E. LARSEN (1984): Limitation of imaging with first-order diffraction tomography, IEEE Trans. Microwave Theory Tech., MTT-32 (8), 860-874. TROYAN, V.N. and G.A. RYZHIKOV (1994): Three- dimensional diffraction tomography functionals, in Proceedings of 2nd Geophysics Congress, Thessaloniki, 340-348. WU, R. and K. AKI (1985): Scattering characteristics of elastic waves by an elastic heterogeneity, Geophysics, 50, 585-592. ZHOU, Q., A. BECKER and H.F. MORRISON (1993): Audio- frequency electromagnetic tomography, Geophysics, 58, 482-495.