Bona, Andrej

Nadri, D.

Urosevic, Milovan

English

espace@Curtin

NONLINEAR ESTIMATION OF THOMSEN ANISOTROPY PARAMETERS IN TI MEDIA   Dariush Nadri1*, Anderj Bóna2, Milovan Urosevic3 CSIRO, Dariush.Nadri@csiro.au1, Curtin University of Technology,A.Bona@curtin.edu.au2  ,Curtin University of Technology, M.Urosevic@curtin.edu.au3  Key Words: Inverse Modelling, Traveltime Inversion, Thomsen Parameters, TI Anisotropy, Preconditioning Conjugate Gradient.  INTRODUCTION  Transverse isotropy with the horizontal axis of symmetry (HTI) is the simplest anisotropic model to characterize a geological formation with vertical fractures.  We formulate an azimuthally dependent parametric equation for P-wave traveltimes in a layered HTI medium. Also, assuming a known vertical velocity model, we estimate simultaneously the fracture orientation and Thomsen’s anisotropy parameters in a stack of horizontal layers using a modified preconditioning conjugate gradient algorithm.  In a numerical example, we show that in the context of weak anisotropy, Thomsen’s anisotropy parameters ( and  ) and fracture orientation can be efficiently estimated using the P-wave traveltimes without a prior information. .   TRAVELTIME EQUATIONS IN TI MEDIA  Ursin and Stovas (2006) proposed a parametric solution of the Christoffel equation for a transversely isotropic medium with vertical axis of symmetry (VTI) by the expression                     2 2 20v p H   ,                                                                                  (1) where v  is the phase velocity, 0  is the vertical  velocity, p is the ray parameter, and H  is an anisotropy parameter, which is a function of the ray parameter and Thomsen’s anisotropy parameters.  Equation (1) can be extended to an HTI medium using the equivalent VTI model formulated by Tsvankin (1997). This is accomplished by introducing the vertical Thomsen’s anisotropy parameters  v and  v , ( )1 2v   ,                                                                        (2)    ( ) 2 1 /1 2 1 2 /v ff      ,                                                              (3) where  0 01 /f     .                                                                   (4) In an HTI medium, P-wave and S-wave velocities depend on the azimuth of the seismic survey with respect to the axis of symmetry. Figure 1 shows the geometry of a ray propagating in a HTI medium with ray (group) angle  . The phase vector lies in the plane constructed from the symmetry axis and the ray vector.   is the phase angle,   is the azimuth, and  is the ray angle of incidence. zS RSymmetry axisCMP line  Figure 1:  A schematic of a HTI medium ( After Tsvankin (1997) with some modifications)  Using the relationship between the phase velocity and the phase angle with the ray velocity and the ray angle  given by Ursin and Hokstad (2003), together with the equation (1) and its  derivatives,  we express the offset  and traveltime for a given ray as:                       2, ,1 ,cosM VTI k HTI k kHTIk VTI kT XTX                                                                    (5)                      2,212 2 2 2 20, , , ,2,2 20, ,1 cos21Mk VTI kHTIkk VTI k k k k kVTI kk kz XXv T p z H pHXp v H             (6) where  0     / 2 , and / 2   defines the symmetry plane of the constant phase velocity. In this expression, , ,k kH dH dp   , kz is the thickness of the layer k, 0v  is the vertical P-wave velocity and M is the number of layers. VTIX  and VTIT  are offset and traveltime for a VTI medium given by  Ursin and Stovas (2006) as:                         , ,, 0, 2 20, ,1121 1k kVTI k k kk kH pHX pz vp v H   ,                                                         (7)                          3 20, ,, 2 20, 0, ,1121 1k kkVTI kk k kp v HzTv p v H .                                                               (8) Equations (7) and (8) have to be evaluated with the Thomsen’s anisotropy parameters in the vertical direction. To show the validity of equations (5) and (6), we have compared the resulting traveltimes with computed traveltimes for a HTI layer using the phase velocity given by Tsvankin (1997) with 1000 ,  1940 / ,  780 / ,  0.03,  and 0.1p sz m v m s v m s       . Figure 2 shows the comparison of the computed traveltimes between Tsvankin (1997) and this method for different azimuths. For longer offsets and small azimuths, Tsvankin’s equation results in smaller traveltimes.  0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000Offset(m)30002800260024002200200018001600140012001000Traveltime(ms)Azimuth:  Tsvankin (1997)    Parametric                   75                   60                   45                   30                     0            Figure 2: Traveltimes using Tsvankin’s equation and equation (5) for a HTI layer at different azimuths.  RAY PARAMETER ESTIMATION IN TI MEDIA  To estimate the ray parameter for a given ray path we minimize an objective function  2r m cf x x   which is the square of the difference between the computed cx  and a measured offset mx using the Newton method (Nocedal and Wright, 1999),  For a given ray parameter we compute the offset by summing the offsets from each ray element using the equations (6) and (7).  A random ray parameter is drawn from a uniform distribution in the range 1.e-05 – 1.e-4 and it undergoes a Newton method of minimization according to:    2 21n n r rp p f p f p        ,                                                           (9) where    2r m cf p x p x x        ,                                                           (10)     22 2 2 22 2r m cf p x p x x x p           .                                 (12) The first and second derivatives of offset with respect to ray parameter are computed analytically.  This method of ray parameter estimation is very accurate, fast and robust even for large values of Thomsen’s anisotropy parameters and for a large change in the vertical P-wave velocity at the interface of two layers.  PARAMETERS ESTIMATION   In contrast to the layer stripping methods in which the model parameters in each layer are minimized separately from the influence of the other layers, we consider the effect of all the rays contributing to a specific layer.  This, however, increase the nonlinearity of the objective function but results in a more robust parameter estimation.  We follow the Dahl and Ursin (1991) approach to compute the partial derivatives of traveltime with respect to model parameters in layered media analytically.  Since we express the traveltimes as a function of ray parameter, the derivative of traveltimes with respect to any model parameter has to consider the derivative of ray parameter with respect to any model parameter as well.  This will be accomplished by taking the derivative of offset with respect to each model parameters as:    .j j jp constdx dm x m x p p m         ,                                 (13) where jm  is the model parameter for a layer  j, x  is the source-receiver offset , and  p is the ray parameter.  Since the offset is given for a given ray, the derivative of the offset with respect to any model parameter is zero, 0jdx dm  .  This result in derivatives of the ray parameter with respect to a given model parameter as:    .j j p constp m x m x p       .                                             (14)  Given the equation (14), the derivative of the traveltime with respect to each model parameter can be expressed as the sum of the derivative of traveltime with respect to a model parameter for a given layer, at a constant ray parameter and a derivative term including the effect of the change in the ray parameter for a ray path as the consequence of change in a model parameter as:    .j j jp constdT dm T m T p p m        .                              (15)  A quadratic objective function, without imposing a specific model space structure given by Tarantola (2005), can be used to minimize the residual error as:    1 11 112TN Mobs syn ii obs synij ij D ij iji jf T T C T T                                        (16)  where TN  is the number of traces, obs  stands for measured traveltimes, and syn  for computed traveltimes.  We assume there is no correlation between the data, hence, the off-diagonal elements are zero.   The diagonal elements or the variances could simply be the errors in picking the traveltimes.  Using an iterative scheme we update the prior model vector m according to:  11Ti i i M D im m C fC e                                                                      (17)  In practice, we assume that the  data, as well as model parameters are independent, hence, the correlation between two different elements is zero.  So, the covariance matrices of model MC  and data DC  are diagonal and contain the variances.  Since, we assume that DC is proportional to identity matrix, we look at the product of gradient T f and data residual vector obs syne d d   1,1 1N MijTS k iji jk kTf e em m       ,                                                          (18)  where the indices i  , j , and k  repeat over the  number of offsets N ,   the number of layers 1M  , and the model parameters respectively.  To minimize the objective function (17) we use a nonlinear preconditioning conjugate gradient algorithm (Nocedal and Wright, 1999; Saad, 2000; Bonnans et al., 2003) with  modifications in search direction scheme by Kalkreuter and Simma (1996),      00 0 0 010 0 0 0 01Given a prior model m  and convergence tolerance >0,Evaluate , ;Set ,  , 0;whilecompute either exact or inexact line search ;;iii i i iif f m f f mh c H r r f irm m hr                1 1 111 11 1 1 111 1;          , , , ;     0          ;i ii ii i i i i i iii if mc H rr c r c r cifh r           1 1 1 1 1 1 1               , , ;      1;i i i i i i i i ielseh c h m m h m mi iend                     where H is the Hessian matrix, h  is the search direction, r  is the gradient  equation (19) ,   is the Polak-Ribiére  scalar,  is the machine precision, i  is the iteration number, and   is the step length and could be computed either by inexact line search criteria such as Wolfe conditions (Nocedal and Wright, 1999) or exact methods such as Brent algorithm (Press et al., 2002).  is the inner product operator.  We may use the objective function to stop the conjugate gradient algorithm.   NUMERICAL RESULTS  To evaluate the robustness of the inversion, we simulate noise free traveltimes for the model given in table 1 over 9 reflectors from a single shot in 96 geophone locations in 50 metres intervals. The nearest offset is at 200 metres, while the farthest is at 4950 metres.  The fracture directions are kept the same for all HTI layers.  We estimate the ray parameter for each ray path in the HTI and VTI layers using the Newton minimization of the offsets in equations (6) and (7), respectively.  In the following minimization we kept the thicknesses and vertical velocities constant and only optimized the Thomsen’s anisotropy parameters and the azimuth of axis of symmetry of HTI layers.  We draw a random prior model from a uniform distributions for the azimuth and the anisotropy parameters from the following intervals: 0 / 2   , 0.05 0.20   and 0 0.20  . We use a modified preconditioning conjugate gradient algorithm (Kalkreuter and Simma, 1996; Nocedal and Wright, 1999) to minimize the traveltimes.  Figure 3 shows the convergence of the objective function, which reached the solution after 60 iterations.  Table 2 shows the estimated model parameters, where all model parameters converged to the true solution.  Figures 4-8 show the convergence pattern of the azimuths for all HTI layers,  and   for layer 6 and 7. Azimuthal dependencies have significant impact on the value of the objective function and converge very fast, while Thomsen’s anisotropy parameters have smaller effect and converge slower.      Table 1:  True Layer properties which has        Table 2:  Prior and estimated model parameters  used for simulating observed  traveltime 0 20 40 60 80 100Iteration-0.1-0.08-0.06-0.04-0.020 20 40 60 80 100Iteration-0.2-0.18-0.16-0.14-0.12-0.1-0.080 20 40 60 80 100Iteration00.040.080.120.160.20 20 40 60 80 100Iteration0.040.060.080.10.120.140.160 2 4 6 8 10Iteration020406080AzimuthLayer 2 Layer 4 Layer 6 Layer 8 0 20 40 60 80 100Iteration1E-0051E-0031E-0011E+0011E+0031E+005Objective Function             Figure 3: Convergence of the objective function.         Figure 4: Convergence of the azimuths.               Figure 5: Convergence pattern of   for layer 7           Figure 6: Convergence pattern of   for (VTI).                                                                           layer 7 (VTI).                Figure 7: Cconvergence pattern of   for layer 6          Figure 8: Cconvergence pattern of   for  (HTI).                                                                           layer 6 (HTI).   CONCLUSION AND DISCUSSION  We have derived parametric traveltime and offset equations for ray tracing purposes for compressional waves in a horizontal transversely isotropic medium with horizontal axis of symmetry, where the elastic properties of the medium are laterally invariant. We have developed an unconstrained minimization algorithm using a preconditioning conjugate gradient method to estimate the fracture orientation and interval Thomsen’s anisotropic parameters in a succession of layers with VTI and HTI anisotropy from the inversion of  P-wave traveltimes. We have modelled  noise free traveltimes, more data from different azimuths will make the inversion robust in case of noisy traveltime.  ACKNOWLEDGMENTS  We would like to Acknowledge Curtin Reservoir Geophysics Consortium to support part of this project. We would like to thank Dr Bruce Hartley (Curtin University of Technology) and Professor Boris Gurevich (Curtin University of Technology) for supervision of one of us (DN).  REFERENCES  Bonnans, J. F., Gilbert, J. C., Lemaréchal, C. and Sagastizábal, C. A., 2003, Numerical optimization: Theoretical and practical  aspects: Springer-Verlag.  Dahl, T., Ursin, B. , 1991, Parameter estimation in one-dimensional anelastic medium: Journal of Geophysical Research, 96, 20217-20233.  Kalkreuter, T., Simma, H., 1996, An accelerated conjugate gradient algorithm to compute low-lying eigenvalues- a study for Dirac operator in SU(2) lattice QCD: Computer Physics Communications, 93, 33-47.  Nocedal, J.,  Wright, S. J., 1999, Numerical optimization: Springer-Verlag. Press, W. H., Teukolsky, S. A. , Vetterling, W. T., and Flannery, B. P., 2002, Numerical recipes in C++: the art of scientific computing: Cambridge University Press.  Saad, Y., 2000, Iterative Methods for Sparse Linear Systems.  Tarantola, A., 2005, Inverse problem theory and methods for model parameter estimation: SIAM.  Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, 1954-1966. Tsvankin, I., 1997, Reflection moveout and parameter estimation for horizontal transverse isotropy: Geophysics, 62, 614-629.  Ursin, B., Hokstad, K., 2003, Geometrical spreading in a layered transversely isotropic medium with vertical symmetry axis: Geophysics, 68, 2082-2091.  Ursin, B., Stovas, A., 2006, Traveltime approximations for a layered transversely isotropic medium: Geophysics, 71, D23-D33. 

Nonlinear estimation of Thomsen anisotropy parameters in TI media

https://espace.curtin.edu.au/bitstream/20.500.11937/25196/2/133217_12850_005ASEG09_DNadri.pdf

Nonlinear estimation of Thomsen anisotropy parameters in TI media

Abstract

Similar works

Full text

Available Versions

espace@Curtin