Short note: Polarization and slowness component inversion for TI anisotropy

Abstract

Seismic anisotropy can be measured using walkaway VSPs. One popular method is through slowness curve construction from cross-plots of traveltime derivatives calculated with respect to offset and depth. This method works where the overburden is relatively simple. However, in areas where the geological structure cannot be considered as approximately one-dimensional, the method fails. For these situations a more general method is required which relies on local homogeneity around the receiver array and which is independent of structural complexities in the overburden. Local wavefield measurements include slowness components resolved along the receiver array direction and polarizations, both of which can be measured in walkaway VSP experiments ( Fig. 1). We invert qP and qSV polarization and slowness component measurements obtained from a walkaway VSP experiment using a global optimization method. Schematic illustration of vertical slowness component and polarization measurements recorded across a VSP receiver array. For qP waves in anisotropic media the propagation direction is generally not the same as the polarization direction . We apply our inversion to a walkaway VSP experiment shot in the Java Sea region as described by Leaney (1994). Data were collected on a five-level non-gimballed three-component receiver tool locked at depths between 1630 m and 1690 m in a near-vertical well. Maximum source offsets were 2.5 km with a 25 m approximate intershot spacing. The geology is essentially horizontal which allows us to construct slowness curves from which we can derive anisotropy estimates using existing inversion procedures ( Miller and Spencer 1994). This provides us with a comparison for our polarization inversion results. After the waveform data are rotated to compensate for horizontal tool twist, little energy is observed on the cross-line geophone component, indicating propagation in a vertical symmetry plane. Given that the geology local to the receiver array consisted of a homogeneous and horizontally layered shale unit, it is reasonable to assume that the local anisotropy can be well approximated by transverse isotropy with a vertical symmetry axis (TIV). Converted shear-arrivals are present within the data and are generated by a 12 m thick carbonate stringer located above the shale unit. Decomposition into qP and qSV wavefields is achieved using a parametric inversion ( Esmersoy 1990) applied to data from three geophone levels centred at the second deepest geophone level. In our implementation of the parametric inversion we include both upgoing and downgoing waves simultaneously ( Leaney 1990) and treat the slowness components and polarizations independently, so that there is no ‘hidden’ assumption of isotropic wavefield behaviour. Conveniently the parametric inversion also yields polarization and slowness component measurements for the separated qP and qSV arrivals. Each shot in the walkaway is processed to give slowness component and polarization measurements for downgoing transmitted qP, upgoing reflected qP, downgoing converted qSV and upgoing reflected qSV converted events. These are shown in Fig. 2 as plots of polarization angle and slowness component (with 90° removed from the shear-wave polarizations for the purposes of display). Although the data show some scatter, particularly in the case of the qSV arrivals, there is a clear trend which does not conform to isotropic behaviour. These are the data we invert to obtain local anisotropy estimates. Vertical slowness components and polarizations for qP and qSV arrivals. Downgoing arrivals have polarizations greater than 90° and negative vertical slowness components. The qSV arrivals have 90° removed for the purposes of display. For a given anisotropic model we use a numerical procedure to find the propagation directions which correspond to the measured qP and qSV polarizations. Once we have found these propagation directions, we calculate the slowness vectors and project these on to the receiver-array direction. This gives a set of modelled slowness components which are compared with the observed slowness components using an appropriate misfit function. The aim of the inversion is to minimize this quantity as a function of the model parameters. The essential part of our approach is finding a propagation direction for a given qP polarization direction. We do this with a simple and apparently robust algorithm which relies on the qP polarization generally being close to the propagation direction. This algorithm is illustrated in Fig. 3 and proceeds as follows: Schematic illustration of the iterative polarization search algorithm. p●0 is the required polarization direction to be found and i is the polarization corresponding to the propagation direction n●i. The difference between p●0 and p●i is Δi. The next propagation direction is calculated by adding Δi to n●i and normalizing to obtain n●i+1. 1 Solve the Kelvin–Christoffel equation ( Musgrave 1970) for the propagation direction n^i which is initially equal to the measured qP polarization direction p^0. This gives the calculated qP polarization direction p^i. 2 Calculate the difference between the calculated qP polarization direction and the required qP polarization direction to give Δi. Thus 3 Add this difference to the propagation direction, so that 4 Renormalize to obtain a new propagation direction n^i+1, where 5 Solve the Kelvin–Christoffel equation for this new propagation direction giving a new calculated qP polarization direction. 6 If this new qP polarization is not close enough to the required polarization direction then iterate from step 2. In numerical tests for TI, orthorhombic and triclinic symmetry systems, the algorithm converged. In the case of TI symmetries where the anisotropy is not too severe, i.e. the qP polarizations are approximately equal to the propagation direction, this algorithm typically converges within five iterations. Convergence is slower for symmetries lower than TI and in more strongly anisotropic models where more than 50 iterations were required. The extension to symmetries lower than TI and our use of non-linear optimization methods distinguishes our work from previous polarization inversions such as those of de Parscau (1991) and Hsu and Schoenberg (1991). It is also possible to improve this algorithm using perturbation methods for the polarization directions ( Pšencik and Gajewski 1998). So far we have focused on the search process for qP polarizations. Unfortunately we cannot use the same approach for qS waves in arbitrary anisotropic media. This is because the two qS waves are orientated in a plane orthogonal to the qP polarization, so that a single qS polarization could correspond to a whole range of directions. Furthermore, for symmetries lower than TI, the qS polarizations can show complicated behaviour due to the presence of point singularities where the shear-wave surfaces touch and there is a rapid change in the polarizations. However, for TIV symmetries, one of the shear waves is polarized in the vertical plane containing the propagation direction (qSV) and the second is polarized perpendicular to this plane and is a true shear wave (SH). We can exploit this information so that, if we have observations of qSV polarizations, we can use our existing qP search algorithm to find propagation directions for a qP polarization which is perpendicular to the required qSV polarization direction. Our inversion uses the four-parameter model, The first two elements define the qP and qSV anisotropy for TIV media ( Carrion et al. 1992 ), and The parameter εP is similar to the Thomsen parameter ε and essentially measures the difference between the horizontal and vertical qP velocities ( Thomsen 1986). The parameter εA defines the anellipticity of the material with εA=0 defining an elliptically anisotropic medium. If εA is positive, then the qP slowness curve bulges outwards and the qSV slowness curve bulges inwards from a circle near 45°. We chose this parametrization scheme because, unlike the Thomsen parameter δ, the anellipticity parameter can be readily interpreted in terms of the geometrical nature of the slowness curves. and are the ratios of the final to the initial estimates of the vertical Vp and Vs velocities which are and respectively. Our inversion method can also include inversion parameters describing a non-vertical symmetry axis, borehole deviation and gain mismatch that may occur between horizontal and vertical geophone components. We chose a least-absolute-values misfit in preference to the more commonly used least-squares criterion because of the apparent presence of outliers in the data. The misfit function to be minimized in the inversion is where Smod are the modelled slowness components and Sobs are the observed slowness components which have estimated errors of ΔS as indicated in Table 1. The estimated slowness component errors are 0.02 s/km and 0.04 s/km for the downgoing and upgoing qP events, respectively. The increased errors for the upgoing events reflect the greater degree of scatter compared with the downgoing data. For similar reasons the estimated qSV slowness component errors are 0.04 s/km and 0.08 s/km for the downgoing and upgoing events, respectively (these error bars are shown in Fig. 5). N is the number of observations. After editing outliers from the data there are 851 observations consisting of 408 and 443 qP and qS measurements, respectively. Comparison of observed and modelled data for the isotropic model and the two anisotropic models from the slowness curve inversion and the polarization and slowness component inversion. Adaptive simulated annealing is used to minimize the misfit function as a function of the model parameters. This technique has the advantage of being an efficient, non-linear and global search method which is not critically sensitive to control parameters ( Ingber 1993). Each run of the inversion samples 2000 models and we run the inversion five times under different initial conditions to improve sampling of the model space. In total we sampled 10 000 models. In each of our five inversion runs the algorithm converged to essentially the same model. This is because the problem is essentially unimodal for this particular inversion as can be seen in the εP versus εA cross-plot shown in Fig. 4. In this case our choice of such a powerful optimization method as simulated annealing is probably excessive. Nonetheless, for other polarization inversions we have found the problem to be multimodal ( Horne et al. 1998 ). Cross-plot of εP versus εA for models sampled by the algorithm with a misfit of less than one. On the scatter plots we have marked the best-fitting isotropic model, the model found by our inversion and the model derived from the slowness curve inversion. The isotropic slowness curve, the slowness curve inversion and the polarization inversion results are compared in Fig. 5 with the corresponding model parameters given in Table 2. The isotropic model clearly represents a poorer fit to the observed data than the anisotropic models. The similarity of the two anisotropic models and misfits suggests that the polarization and slowness component data are consistent with the slowness curve data. Furthermore, if we compare the modelled data for the two anisotropic models there appears to be very little difference in the qP modelled data and only a relatively small difference in the qS modelled data. This relative insensitivity of the qP data to the anisotropy can be shown using cross-sections through the model space about the best model as shown in Fig. 6. Here we have plotted the contributions to the misfit function from the qP and qS observations as we have varied εP and also εA. It can be seen that the qP misfit appears much flatter than the corresponding qS misfit. This implies that the qP data are not as sensitive to the anisotropy parameters as the qS data. Inclusion of parameters describing gain mismatch and non-vertical symmetry axis orientation did not produce significantly better results. Cross-sections about the best model showing the qP, qSV and total misfits. The arrows indicate the best-fitting model parameters. The shaded area indicates the region of acceptable models which have a misfit of less than one. Our method has also been tested in the inversion of qP polarizations from multi-azimuthal walkaway VSPs. Here the objective was to determine fracture-induced anisotropy estimates ( Horne et al. 1998 ). In this case the inverted fracture orientation is in good agreement with a priori geological and geophysical results. This provides us with further confidence that polarizations and slowness components can be inverted to give realistic anisotropy estimates. Polarization and slowness component measurements can be inverted to yield realistic local anisotropy estimates. Such a technique will be most useful where conventional slowness curves cannot be constructed because of structural complexities in the overburden. Although the slowness component along the receiver array and the polarizations are local wavefield measurements, they will require reliable measurement which may be difficult in such complex heterogeneous regions. Even when the overburden is sufficiently simple to allow both horizontal and vertical slowness components to be calculated, it will be advantageous also to invert for polarizations. We are grateful to P.T. Maxus Southeast Sumatra for the release of this data and to Schlumberger for allowing us to publish this work.