Mathematical modeling to reconstruct elastic and geoelectrical parameters

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.


Introduction
Diffraction tomography is an imaging technique 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 resolution 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 development of the diffraction tomography is introduced 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).
The applicability of the first order Born approximation in diffraction tomography is shown (Slaney et al., 1984) by numerical simulation 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 sourcereceiver pairs.We also implement multiparametric 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 multiparametric reconstruction in the elastic case (reconstruction of the Lame parameters, mass density and as a corollary -shear and longitudinal 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 observation 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 electromagnetic cases.

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 (2.2) the wave field u(x, t) 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 − u rf 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) 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 where d u are the samples of the scattered field.The final version of this system after introducing regularizing terms reads as 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.

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) q p s j r j q s q j q p s j

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 The reference medium (rf ) is supposed to be known (σ rf ,ε rf ,µ rf ) and electrical field E rf satisfies to the equation As in the case of scattering by elastic inhomogeneities 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 − E rf 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: where in a two-dimensional case the sources from (2.4) are used.
Introducing the tomography functionals (3.8) ( ) 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.˜, ,  , ,   , , ,  ˜, ,  , ,   , , , ( and 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 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).

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 The values (from (3.13)) describe the space characteristics of the scattered fields, which are produced by the elementary inhomogeneities ( , , µ) 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 conductivity.
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 ).
) Graphic representation of directivity diagrams for the elastic and electromagnetic cases can be found in Wu and Aki (1985) and Saintenoy and Tarantola (2001).The results of reconstruction of p are represented 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 reconstruction 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 handside 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.

Electromagnetic case
The model of inhomogeneity of electrical conductivity together with the results of reconstruction 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 .value of electrical conductivity is 10 3 S/m.An apparent frequency of the sounding signal is 5 × 10 6 Hz.We use the observation scheme which is similar to the observation scheme represented 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 dependence on the values of the regularizing coefficients in comparison with the elastic case.This difference can be explained by greater distortion of a shape of the scattered signal in the case of reconstruction of electrical conductivity.
of n = 0 the values , and are represented by the amplitude factors of the zero approximation A 0i , A and eikonals ~, (from ( Elastic caseAs 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.

Fig. 1 .
Fig. 1.Model of medium and observation scheme.S is the region for reconstruction; location of inhomogeneityblack color; 1, 2, 3 are the source and the observation points locations; s ; s / p = 1/ .
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.

Fig. 4a -
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).