Theory of interval traveltime parameter estimation in layered anisotropic media

Moveout approximations for reflection traveltimes are typically based on a truncated Taylor expansion of traveltime squared around zero offset. The fourthorder Taylor expansion involves NMO velocities and quartic coefficients. We derive general expressions for layer-stripping both secondand fourth-order parameters in horizontally-layered anisotropic strata and specify them for two important cases: horizontally stacked aligned orthorhombic layers and azimuthally rotated orthorhombic layers. In the first of these cases, the formula involving the out-of-symmetry-plane quartic coefficients has a simple functional form and possesses some similarity to the previously known formulas corresponding to the 2D in-symmetry-plane counterparts in VTI media. The error of approximating effective parameters by using approximate VTI formulas can be significant in comparison with the exact formulas derived in this paper. We propose a framework for deriving Dix-type inversion formulas for interval parameter estimation from traveltime expansion coefficients both in the general case and in the specific case of aligned orthorhombic layers. The averaging formulas for calculation of effective parameters and the layer-stripping formulas for interval parameter estimation are readily applicable to 3D seismic reflection processing in layered anisotropic media.


INTRODUCTION
Moveout approximations play an important role in conventional seismic data processing (Yilmaz, 2001). Bolshykh (1956) and Taner and Koehler (1969) laid the groundwork for studies on moveout approximations by proposing to employ the Taylor expansion of reflection traveltimes around zero offset. This approach has led to many developments in traveltime approximations for isotropic and anisotropic media (Malovichko, 1978;Hake et al., 1984;Sena, 1991;Tsvankin and Thomsen, 1994;Alkhalifah and Tsvankin, 1995;Alkhalifah, 1998;Grechka and Tsvankin, 1998;Taner et al., 2005;Ursin and Stovas, 2006;Blias, 2009;Fomel and Stovas, 2010;Aleixo and Schleicher, 2010;Golikov and Stovas, 2012;Stovas, 2015). Some of the early developments are summarized by Castle (1994). In the small range of offsets, the reflection traveltime has the well-known Interval parameter estimation hyperbolic form and its processing involves only one controlling parameter, namely the NMO (normal moveout) velocity. In the case of horizontally stacked isotropic layers, the effective NMO velocity can be related to the interval velocity through Dix inversion (Dix, 1955). Its 3D counterpart was described as generalized Dix equation Tsvankin, 1998, 1999;Tsvankin and Grechka, 2011) In the case of long-offset seismic data and more complex media, both the hyperbolic moveout approximation and the Dix inversion formula need to be modified due to the moveout nonhyperbolicity (Fomel and Grechka, 2001). Hake et al. (1984) studied this problem in VTI (vertically transversely isotropic) media and proposed a 2D averaging formula for the quartic coefficient of the traveltime-squared expansion. Tsvankin and Thomsen (1994) introduced a different functional form for nonhyperbolic moveout approximation with correct asymptotic behavior at large offsets and propose a Dix-type formula for inversion of the quartic coefficient for interval anisotropic parameters. The moveout approximations by Tsvankin and Thomsen (1994) and its 3D extensions have led to subsequent developments on this topic (Alkhalifah and Tsvankin, 1995;Pech et al., 2003;Pech and Tsvankin, 2004;Xu et al., 2005). Various alternative nonhyperbolic moveout approximations have been investigated in the literature. Several of them can be related to the generalized moveout approximation (Fomel and Stovas, 2010;Sripanich and Fomel, 2015a;Sripanich et al., 2016).
In the case of 3D orthorhombic media, several well-known moveout approximations make use of the rational approximation  in combination with the weak anisotropy assumption (Pech and Tsvankin, 2004;Xu et al., 2005;Vasconcelos and Tsvankin, 2006;Grechka and Pech, 2006;Farra et al., 2016). This kind of approximation is readily applicable for a single homogeneous orthorhombic layer. For a horizontally layered model, the rational approximation suggests averaging azimuthally dependent interval quartic coefficients using expressions for VTI media (Hake et al., 1984;Tsvankin and Thomsen, 1994;Tsvankin, 2012). However, this is justifiable only when the azimuthal anisotropy is mild Vasconcelos and Tsvankin, 2006). In this paper, we derive exact expressions for averaging interval quartic coefficients in a 3D horizontal stack of general anisotropic layers. Next, we specify these expressions explicitly for two particular settings: a horizontal stack of aligned orthorhombic layers and a horizontal stack of azimuthally rotated orthorhombic layers. These expressions lead to exact Dix-like layer-stripping formulas for interval parameter estimation in layered anisotropic media.

TRAVELTIME EXPANSION
Assuming the Einstein repeated-indices summation convention, we can expand the one-way traveltime t into a Taylor series of half offset h i (i = 1 or 2 in 3D) around zero offset as follows: where t 0 is one-way vertical traveltime, t ij = ∂ 2 t ∂h i ∂h j and t ijkl = ∂ 4 t ∂h i ∂h j ∂h k ∂h l are secondand fourth-order derivative tensors, respectively. Both tensors are symmetric thanks to the symmetry of mixed derivatives. Analogously, we can also derive, for the negative half offset −h i , Assuming pure-mode reflections with source-receiver reciprocity, we can sum the two expansions (equations 1 and 2) for the two legs of rays to derive the expansion of the two-way traveltime as follows : Equation 3 can be additionally transformed into the series of the squared two-way traveltime in terms of the full offset x i = 2h i as follows: where In consideration of the symmetry of the time derivative tensors, the quadratic and quartic terms in equation 4 reduce to the following known expressions : In the derivation of the general formulas for moveout coefficients in the next section, we keep the tensor notation, which simplifies the use of tensor operations. We also use the fact that, in the case of horizontally stacked layers, the half-offset h i and reflection traveltime 2t can be expressed in terms of horizontal slownesses (ray parameters) p 1 and p 2 in h 1 and h 2 directions as follows: 2t(p 1 , p 2 ) = 2 where D (n) and Q (n) (p 1 , p 2 ) denote the thickness and the vertical slowness of the n-th layer. The derivation of equations 9 and 10 is included in the appendix. The general dependence Q (n) (p 1 , p 2 ) follows directly from the Christoffel equation. Throughout the text, we use the subscript index in parentheses to indicate the corresponding layer. The upper-case and lower-case letters denote interval and effective parameters respectively.

GENERAL FORMULAS FOR TRAVELTIME DERIVATIVE TENSORS
Using equations 9 and 10 and applying the chain rule, we can differentiate the one-way traveltime t with respect to half offset h i to derive the following equations: where the derivatives with respect to p 1 and p 2 are represented by comma (e.g, ∂t ∂p i corresponds to t ,i ), δ ij denotes the Kronecker delta, g ij denotes ∂p i ∂h j , andî,ĵ,k,l represent dummy indices. Equations 11-14 can be used to compute t ij and t ijkl terms needed by equations 5 and 6 using explicit relationships for h(p 1 , p 2 ) and t(p 1 , p 2 ).
According to the chain rule and the symmetry of the time derivative tensors, the second-derivative tensor g ij and its derivatives in equations 13 and 14 can be related to the derivatives of half offset h as follows: g ij,kl = 2g im h mm,k gm n h nn,l gn j − g im h mm,kl gm j , where m,m, n,n are dummy indices. The matrix h ji = ∂h j ∂p i is the inverse of the matrix g ij (Grechka and Tsvankin, 1998). Substituting equations 15 and 16 into equations 13 and 14, we subsequently arrive at expressions which only involve derivatives of explicitly defined functions. Subsequently, we have at zero offset (p i =0):

Interval parameter estimation
We can now substitute expressions for time derivative tensors from equations 21-24 into equations 5 and 6 and rewrite them as follows For convenience in subsequent derivation, we let elements of the matrix inverse of a ij be denoted as Equations 25 and 26 represent the most general forms of the traveltime derivatives for pure-mode reflections in arbitrary anisotropic layered media. Throughout the rest of this text, we continue to use the upper-case letters (A, B, ...) to denote interval parameters. The lower-case letters (a, b, ...) are used to indicate effective values corresponding to the same quantities. In the case of effective values, the subscript refers to the effective value from the surface down to the bottom of the corresponding layer. Therefore, in this notation, we can denote the interval expression in equation 25 for the n-th layer and the effective counterpart down to the n-th layer as A ij(n) and a ij(n) respectively.
Let us first consider the second-order term in equation 27. In the case of multiple layers, there are direct accumulations of traveltime and offset as We can deduce where T 0(n) and B ji(n) denote the vertical one-way traveltime and the inverse of interval matrix A ij(n) in the n-th layer, and t 0(N ) and b ji(N ) denote the effective values of the same two parameters at the bottom of the N -th layer. Therefore, we can find the interval B ji(N ) of the N -th layer from Equation 30 is referred to as the generalized Dix equation (Tsvankin and Grechka, 2011). Applying matrix inversion to its result produces the second-order interval coefficients (A ij(N ) ), which are related to 3D NMO ellipse parameters. In the isotropic case, equation 30 reduces to classic Dix inversion (Dix, 1955). Interval parameter estimation We can now turn to equation 26 for the quartic coefficients and follow an analogous procedure. In this expression, only h mm,kl term needs to be considered because other terms can be simply related to b ji . Equation 26 leads to By substituting equations 25 and 27, we arrive at which can be used to find the interval H mm,kl(N ) in the N -th layer by a simple subtraction: We can then compute the interval quartic coefficient A ijkl(N ) in the N -th layer from the interval form of equation 26 given the interval value of equation 25 obtained from the generalized Dix equation (equation 30), as follows: where T 0(N ) = t 0(N ) −t 0(N −1) denotes the vertical one-way traveltime in the N -th layer. Thus, the second-and fourth-order interval coeffcients for the traveltime expansion can be found from equations 30 and 34 respectively, The exact expressions for the traveltime expansion coefficients in two particular cases of a horizontal stack of aligned orthorhombic layers and a horizontal stack of azimuthally rotated orthorhombic layers are detailed in the subsequent sections.

COEFFICIENTS OF TRAVELTIME EXPANSION FOR A STACK OF HORIZONTALLY ALIGNED ORTHORHOMBIC LAYERS
In the special case of horizontally stacked homogeneous orthorhombic layers with aligned symmetry planes, we simplify the equations derived in the previous section and compare them with some known expansion coefficients.

General formulas for coefficients
Considering equations 9-14 in the case of an orthorhombic medium, only certain derivatives of half-offset h i and t with respect to ray parameters p 1 and p 2 are nonzero at zero offset. The non-zero terms are listed in Table 1.
As a consequence, equation 7 reduces to the expressions given below and derived previously by . Considering the case where the [x 1 ,x 3 ] plane corresponds to one of the vertical symmetry planes of an orthorhombic medium,  show that, for pure modes, where v 2 2 and v 2 1 denote NMO velocities squared associated with x 1 and x 2 directions, respectively. The quartic term can be expressed as where, in the notation of the previous section, the exact expressions of these coefficients are given by and Equations 37-41 appear somewhat cumbersome, we can simplify them further by relating the derivatives of t to derivatives of h i according to equation 10. Therefore, we can conveniently express both derivatives in terms of the derivatives of vertical slowness Q (equation 9) thanks to implicit differentiation. By using the notation Interval parameter estimation we rewrite equations 37-41 in a simpler form as follows: Alternatively, equations 43-47 can also be derived directly from equations 25-26. They provide an easy way to adding layers by a simple accumulation in the ψ i,j parameter.  Vasconcelos and Tsvankin, 2006), is no longer needed.

Formulas for interval NMO velocities
Equations 30, 43 and 44 lead to the following Dix-type inversion formula where a ii(N ) and a ii(N −1) with i = 1, 2 denote the coefficients for the reflections from the top and bottom of the N -th layer respectively, and A ij(N ) denotes the interval coefficient in the N -th layer. t 0 is vertical one-way traveltime with the similar notation rules. We can convert equation 48 to a familiar form (Dix, 1955) by expressing it in terms of interval moveout velocity V 2 1 = 1/A 22 and V 2 2 = 1/A 11 as follows: Equation 49 is a scalar form of the tensor equation 30 for the special case of aligned orthorhombic layers. Interval parameter estimation

Formulas for interval quartic coefficients
Similarly to the derivation for A ii , we have first, which can be calculated for interval Ψ 4,0(N ) in the N -th layer using, Subsequently, interval A 1111(N ) in the N -th layer can be computed from, where T 0(N ) = t 0(N ) − t 0(N −1) . Equation 52 is similar to the Dix-type formula proposed by Tsvankin and Thomsen (1994) for the VTI case. Similar expressions can be derived for A 2222 by considering Ψ 0,4 and A 22 instead of Ψ 4,0 and A 11 . To derive the corresponding expression for A 1122 , we follow an analogous procedure ,which leads to Using Ψ 2,2(N ) = ψ 2,2(N ) − ψ 2,2(N −1) , A 1122(N ) in the N -th layer can be computed as Comparison with known expressions for quartic Taylor series coefficients  derived the following expressions for a 1111 and a 2222 : where f 1 = 1 − v 2 S1 /v 2 P 0 and f 2 = 1 − v 2 S0 /v 2 P 0 , v P 0 is the vertical qP-wave velocity, v S0 is the qS-wave velocity polarized in [x 1 , x 3 ], and v S1 is the qS-wave velocity Interval parameter estimation polarized in [x 2 , x 3 ]. The subscript in anisotropic parameters denotes the index of the normal direction to the associated plane. Note that equation 56 is similar to the VTI expression with corresponding changes in the anisotropic parameters for different planes. This result is due to the similarity of the Christoffel equations in the case of VTI media and in-symmetry-plane orthothombic media (Tsvankin, 2012). The expression for a 1122 is more complicated and is given in . These expressions are equivalent to the equations 45-47.
If the pseudoacoustic assumption (Alkhalifah, 2003) is applied or ,equivalently, f 1 = f 2 = 1, the expressions for the three coefficients simplify as follows: where η 1 , η 2 , and η 3 are the time-processing parameters suggested by Alkhalifah (2003), v i denote the NMO velocities in the plane defined by x i -axis, and denoting an anelliptic parameter responsible for the fourth-order cross term suggested by Stovas (2015). To obtain effective values of the moveout coefficients from interval value, an averaging formula can be used. In the particular case of horizontally stacked VTI layers, the commonly used averaging formula can be expressed in this paper's notation as (Hake et al., 1984;Tsvankin and Thomsen, 1994), where V (n) , A (n) , and T 0(n) denote the interval NMO velocity, VTI quartic coefficient, and one-way vertical traveltime in the n-th layer specified in a given stack of VTI layers. Several schemes have been proposed to obtain estimates of effective quartic coefficients from this averaging formula (equation 62) in the setting of stacked orthorhombic layers. Because of their approximate nature, the application of this expression in orthorombic media is suggested to be valid only if the azimuthal anisotropy is not severe Vasconcelos and Tsvankin, 2006).
To test the accuracy of the previously suggested approximations for computing effective coefficients based on equation 62 against the exact expressions given in equations 45-47, we consider two three-layered models with parameters given in Table 2. The respective layer thicknesses are 0.25, 0.45, and 0.3 km. We compute the quartic term (equation 36) at the offset radius of 1 km and different azimuths. Interval parameter estimation Model 1 (Weak)  Table 2: Parameters of sample orthorhombic models under the notation of Tsvankin (1997) with the subscript denoting the normal direction to the corresponding plane: Model 1 (Weak anisotropy) parameters are extracted and modified from a layered model by (Ursin and Stovas, 2006;Stovas and Fomel, 2012) and Model 2 (Moderate-Strong anisotropy) parameters are modified from Schoenberg and Helbig (1997) and Tsvankin (1997).

Comparison with the method of Al-Dajani et al. (1998)
Al-  suggest to compute the effective quartic term using the VTI formula given by equation 62. The independent parameters in each layer are substituted by their azimuthally variant counterparts as follows: A(α) = A 1111 cos 4 α + A 1122 cos 2 α sin 2 α + A 2222 sin 4 α , where α denotes the azimuth angle measured from [x 1 ,x 3 ] plane, A 1111 , A 1122 , and A 2222 are given in equations 58-60. The effective a (N ) computed based on equation 62 is then multiplied by the source-receiver distance along the CMP line (offset in polar coordinates) and summed with lower-order terms for traveltime computation in polar coordinates. The relative error in the quartic term computation is shown in Figure 1a.
Comparison with the method of Xu et al. (2005) Under the weak anisotropy assumption, Xu et al. (2005) show that the azimuthally dependent quartic coefficients in each layer can be approximated by where V denotes the interval azimuthally varying NMO velocity (equation 63) and η i denotes the interval η value. Vasconcelos and Tsvankin (2006) suggest that in the case of a mild azimuthal anisotropy, one can also use equation 62 to approximate the effective value of A. Therefore, the independent parameters in each layer are substituted by their azimuthally variant counterparts (equations 63 and 65). Analogous to the case of equation 64, the effective A is then multiplied by the source-receiver distance along the CMP line in order to evaluate traveltime in polar coordinates. Interval parameter estimation Figure 1 shows the resultant azimuthal error plots of the quartic term a ijkl x i x j x k x l /r 4 where r = x 2 1 + x 2 2 denotes the offset along the CMP line. Two models of threelayered aligned orthorhombic models with different degree of anisotropy are considered and the model parameters are given in Table 2. qP reflection traveltimes and errors from including traveltime approximation up to hyperbolic and quartic terms are shown in Figure 2 and 3. The VTI based approximation produces small errors when applied in weak anisotropic media. However, the errors increase noticeably with the strength of anisotropy of the media. Note that the observed errors result from the cumulative effect of the use of the pseudoacoustic approximation (equations 58-60) and of the VTI averaging formula (equation 62). The separate effects from pseudoacoustic approximation alone are shown in the same figures and are denoted with large-dashed blue lines. It can be seen from the results that the error from pseudoacoustic approximation dominates the error from the VTI averaging formula in media with a higher degree of anisotropy. The total error on the quartic term a ijkl x i x j x k x l can be computed from the azimuthal error amplified by r 4 and grows with larger distance r along the CMP line.

Comparison of interval parameter estimation
Next, we test the accuracy of the proposed interval parameter estimation formulas (equations 30 and 34) in Model 2 (Table 2) Table 2, b) the difference between the exact reflection traveltime and the hyperbolic moveout, and c) the difference between the exact reflection traveltime and the quartic moveout. Since the reflector depth at the bottom of the third layer is equal to 1 in this case, the nonhyperbolicity effects should have noticeable when x 1 and x 2 are greater than 1. Taking into account the quartic term improves the accuracy of traveltime approximation at large offsets. Interval parameter estimation  Table 2, b) the difference between the exact reflection traveltime and the hyperbolic moveout, and c) the difference between the exact reflection traveltime and the quartic moveout. Since the reflector depth at the bottom of the third layer is equal to 1 in this case, the nonhyperbolicity effects should have noticeable when x 1 and x 2 are greater than 1. Taking into account the quartic term improves the accuracy of traveltime approximation at large offsets.
interval parameters from known effective parameters are shown in Table 3

COEFFICIENTS OF TRAVELTIME EXPANSION FOR A STACK OF AZIMUTHALLY ROTATED ORTHORHOMBIC LAYERS
In this section, we specify the exact expressions for the coefficients of the traveltime expansion in the case of a stack of azimuthally rotated orthorhombic layers.
The following formulas are a specification of the general coefficient formulas for such coefficients given in equations 25 and 26. The quadratic and quartic terms in the traveltime expansion become a ij x i x j = a 11 x 2 1 + a 12 x 1 x 2 + a 22 x 2 2 , and, a ijkl x i x j x k x l = a 1111 x 4 1 + a 1112 x 3 1 x 2 + a 1122 x 2 1 x 2 2 + a 1222 x 1 x 3 2 + a 2222 x 4 2 .

DISCUSSION
In this study, we have focused only on pure-mode reflections where the source-receiver reciprocity holds and the moveout approximation around zero-offset is an even function (equations 3 and 4). This symmetry is generally absent in the case of converted waves.
In addition, we assume that there is no relection dispersal, which is equivalent to assuming that the one-way traveltime can be expressed in terms of half-offset only t(h). This assumption is appropriate when the two legs of the ray path are symmetric with respect to zero-offset and is related to the case of laterally homogeneous, horizontal anisotropic layers with a horizontal symmetry plane. When this assumption does not hold, for example, where tilted anisotropic media without horizontal symmetry planes or arbitrarily shaped interfaces are considered, one-way traveltime is no longer a function of merely half-offset h but also a function of reflection point y(h). Grechka and Tsvankin (1998) emphasized that the effect of reflection dispersal can be neglected in consideration of NMO velocity (Hubral and Krey, 1980). However, reflection dispersal becomes important when higher-order coefficients are considered. The general form for the quartic and higher-order coefficients that honor reflection dispersal effects was first studied by Fomel (1994) based on the extension of Normal-Incident-Point (NIP) theorem (Gritsenko, 1984) to the higher-order. Similar derivation for the quartic coefficientis given by Pech et al. (2003).
In general, with the assumption of laterally homogeneous and horizontal layers, one can conveniently express one-way traveltime and half-offset in each sublayer as functions of horizontal slownesses (p i ) as shown in equations 9 and 10. Their derivatives in the most general form are given in equations 17-20. In consideration of pure modes, the zero-offset corresponds to p i = 0 and, which results in the simplified expressions in equations 21-24. In the case of converted waves or lateral variations or anisotropic media without horzontal symmetry plane, the location of p i = 0 no longer corresponds to the zero-offset but instead to the location of the minimum traveltime. The traveltime derivatives in equations 17-20 can still be used with proper specifications of which position for the derivatives to be evaluated at. On the other hand, if the interfaces are arbitrarily shaped, the appropriate traveltime derivatives can be computing from ray tracing with the knowledge of approximate group velocity in each sublayer Fomel, 2014, 2015b). The proposed framework has a straightforward extension from quartic to higher-order coefficients as long as all the mentioned assumptions are satisfied.

CONCLUSIONS
We have derived novel exact formulas for averaging and inverting the interval quartic coefficients of the reflection traveltime expansion in a general layered anisotropic medium. Expressions for the specific case of a stack of aligned orthorhombic layers are compared with the previous approximations of using VTI averaging formulas. The Interval parameter estimation specific expressions in the case of a stack of azimuthally rotated orthorhombic layers are also provided. The proposed formulas for both averaging interval coefficients and their direct inversion are readily applicable to 3D seismic processing in layered anisotropic media.

ACKNOWLEDGMENTS
We are grateful to A. Stovas, M. Vyas, and an anonymous reviewer for helpful comments and suggestions. We thank the sponsors of the Texas Consortium for Computational Seismology (TCCS) for financial support. The first author was also supported by the Statoil Fellows Program at the University of Texas at Austin.

APPENDIX: TRAVELTIME AND OFFSET AS FUNCTIONS OF RAY PARAMETERS
In this appendix, we show the derivation of the offset (equation 9) and traveltime (equation 10) functions in terms of two ray parameters (p 1 and p 2 ) in 3D. The total offset is constituted of offset increments from each individual layer and can be expressed as where the derivative dh i dz represent the change in the offset h i direction with respect to the vertical direction z, and D (n) denote the thickness of the n-th layer. According to the ray theory (Červený, 2001), this derivative can be related to the derivative of the Christoffel equation with respect to the ray slownesses as follows where F (n) (p 1 , p 2 , Q (n) ) = 0 and Q (n) denote the Christoffel equation and the vertical ray slowness of the interested wave mode in the n-th layer. Equation A-2 can be simplified further due to implicit differentiation in the Christofel equation as which after substitution in equation A-1 results in the function of h(p 1 , p 2 ) given in equation 9. Analogously, we can follow the same line of reasoning and derive an expression of traveltime. We start from the total traveltime expression given by where the derivatives ∂t ∂h 1 = p 1 and ∂t ∂h 2 = p 2 represent the ray parameters in the two directions of the local coordinates. Since the ray parameters are conserved in the sequence of horizontal layers due to Snell's law, we can further transform equation A-4 into equation 10 as follows: where dt dz = Q (n) denotes the vertical slowness in the n-th layer.