VSP data inversion for vertical velocity gradient and elliptical anisotropy model

Inversion of velocity parameters for the walkaway VSP data in a multilayered medium can be impeded by velocity gradients and anisotropy in some layers. A problem occurs if we compare velocities obtained from borehole seismic profiling which are equal to their vertical components with the velocities calculated with paths coming from far offsets where the horizontal component plays an important role, especially when the vertical gradient exists and the ray paths are curve-shaped. In this contribution we present the results of velocity model inversion for VSP data considering velocity gradient and elliptical anisotropy. The algorithm consists of two steps, optimization of velocity parameters and optimization of ray paths for the given model. Both procedures use the Nelder-Mead simplex method which finds local minima. Due to the character of optimization we performed also multistart analysis which can provide information about possible equivalences between parameters. Analysis was conducted for different parameterizations, in some cases allowing introduction of additional parameters: vertical gradient and elliptical anisotropy coefficient. The optimal model for a specific set of data is chosen with the help of Bayesian Information Criterion to balance complexity of model with quality of approximation of traveltimes.


Introduction
Correct inversion of the velocity model is one of the most important problems in seismic data analysis. In the simplest version, geological model is composed of several parallel layers characterized with one parameter: velocity. Commonly, the layers has not one homogeneous velocity but it varies on depth or direction of wave propagation. In practical investigations, simplifications of real geological conditions are used. One of them is vertical velocity gradient which is especially useful in deposit rocks where velocity increases with depth due to compaction. Another factor influencing obtained seismic data is elastic anisotropy. Changes of velocity with direction of wave propagation are common, but not always plays significant role. Especially important is model of vertical transverse isotropy where velocity is changing only with respect to vertical angle of propagation. One of the useful simplification to describe this phenomenon is model of elliptical anisotropy utilized in many recent researches, e.g. [1].
In this paper we used the model of elliptical anisotropy combined with vertical velocity gradient to inverse velocity model from walk-away VSP data from Newfoundland Shelf. Since the model was composed of three layers, ray tracing based on Fermat's principle was used. Optimal number of parameters was chosen with Bayesian Information Criterion to achieve balance between data fitting and model complexity.

Elliptical anisotropy
We consider a medium described with 3 parameters: vertical velocity, its gradient and elliptical anisotropy coefficient. The last of them is commonly used to simplify the phenomenon of anisotropy and formulate it as a one value. It is especially useful in horizontally layered models where velocity of wave propagation in horizontal direction (along the laminas) is slightly bigger than in vertical direction. If the layers are dipped with small angles, the same approximation can be used. Elliptical anisotropy coefficient connects the horizontal ( ℎ ) and vertical ( ) velocities with the formula [1]: This model is simplification of Thomsen's model with two anisotropy coefficients δ and ε fulfilling the condition that both has the same value [2]. In the medium characterized as VTI, this approximation of the real situation is acceptable.
The value of the P-wave velocity is dependent on the vertical angle θ of wave propagation. With the vertical velocity 0 (equivalent of (0)), the formula for the actual wave velocity is given by Thomsen [3]: where δ is coefficient connecting elements of elastic stiffness tensor : The Rogister and Slawinski [4] formulated analytical formula for the traveltime in such a medium, with the source in the point (0,0) and receiver in the point (x, z): The mentioned formulae are valid for one layer only. If one has to calculate traveltime in multi-layered medium, the ray path optimization have to be performed (for two layer model  analytical solution is given in [4]). The algorithm of finding correct ray path is described in the next section.

Inversion algorithm
Since the ray paths in medium characterized by vertical velocity gradient and elliptical anisotropy are not straight lines, we decided to use ray tracing method not basing on the Snell's law. The real wave path through the boundaries can be retrieved with the Fermat's principle used directly. It means that for every calculation of the traveltime one has to find the path giving the least time between source and receiver located in specified coordinates.
The base software algorithm consisted of two-step optimization. At first, the starting velocity model is prepared. For the current model, optimization of ray path is performed with the target function equal to traveltime and X-coordinates of crossing through boundaries being subject of optimization. When the traveltimes for all pairs of sources and receivers are correctly calculated, the target function value for velocity model optimization is found as a sum of squared differences between modelled and measured traveltimes. Then, velocity parameters are changed according to the target function value and new ray paths, for new velocity model are computed.
For both these optimizations, Nelder-Mead simplex method is used. The method has local character and is not based on gradient so it can be used also for undifferentiable target functions. In course of optimization, a polyhedron called simplex is created in the space of solutions. Through the modifications of its vertices the simplex converges to the point with minimum of target function. For the detailed description see [5].

Bayesian Information Criterion
Since the geophysical data are always burdened with measurement errors one should not use the model describing the reality in too much detailed way. Excessive simplification of model is also not desired. In order to balance these two restrictions, some objective criterion is necessary. Commonly used tool is Bayesian Information Criterion (BIC) [6]. The BIC value grows either with number of parameters or growth of target function so it ensures optimal number of model parameters (it is more critical than other information criterion e. g. Akaike). The original formula for BIC value is: where L is maximized likelihood, M is the number of observations and k is the number of parameters in the model. According to [7] we can change searching for maximum probability in (6) to finding minimum of following formula: where 2 is error variance. The method is used also in other similar researches, like [8].

Data description
In this paper we used data from VSP measurements on the Newfoundland Shelf presented in [9]. The core of the dataset is walk-away VSP data recorded with 5 receivers inside the borehole at the depth between 1980 m and 2020 m. The airgun was source of a wave, shots were performed with interval about 25 meters and maximal offsets reaching 1000 meters at "short side" and 4000 meters at "long side".
The zero-offset VSP data (Fig. 1) were used to create a starting model for the inversion. According to smoothed interval velocity curve interpreted from the traveltimes, the velocity model should be divided into three layers (in its part significant for walk-away VSP data inversion). The boundaries are visible at the depth of about 1300 meters and 1750 meters and in starting model of inversion these values are used as the fixed depths.

The results of inversion
The optimization of velocity model was performed for many starting models which differ on number of parameters, starting values and starting steps of optimization (as an element of Nelder-Mead algorithm). In the course of tests, the details of the software have changed, e.g. in order to avoid problems with vertical rays the near-offset data were removed from the dataset used to calculate target function. For the near-vertical segments of seismic ray, the numerical errors occurs since the precision of argument of atanh function in equation (4) is limited. Similar conclusions were reached in [9], thus author suggests to not include nearoffset data in the calculations. Excluding the cases disturbed by those numerical errors, all remaining optimizations, despite of many different starting models lead to the traveltimes consistent with measured ones (Fig. 2). The seismic ray paths obtained as a result of optimization based on Fermat's principle seems also correct -the angles of wave propagation for neighbouring rays are consistent and follow the rules connecting velocities and ray angles (Fig. 3).

The multistart analysis
Since the optimization has a local character, the multistart analysis was performed to ensure that obtained results are exactly the best models. The whole procedure of velocity model inversion was repeated 1000 times with differing starting values of parameters for each type of model regarding number of parameters. Depending on the starting model, different values were obtained, however some dependencies between parameters values can be noticed. In the Fig. 4 the values of parameters obtained as a result of multistart analysis for 7 parameter model are presented as crossplots. The data exhibits strong correlation between a and b parameters of the first layer (correlation coefficient about -0.966). Negative correlation means that lower velocity at the depth of zero is compensated with increased vertical gradient. This dependence is not related to number of model parameters and value of elliptical anisotropy coefficient if it is considered. Similarly, in the second layer correlation between a and χ parameters has coefficient about -0.679. In this case gradient does not compensate changes of velocity at the top of the layer.
The comparison of the target function value achieved with different values of specific parameters also provide interesting remarks (Fig. 5). The values of parameters a0 and χ1 in 8 parameter model plotted against target function value create curve-shaped optimization fronts. Some of the results demonstrate that with given value of some parameter the target function cannot be lower than specific number even with changes of another parameters.

BIC analysis
For all the best results for each parameterization, BIC value was calculated to choose the most suitable model with equation (7) (see section 2.3). The obtained values shows that 7 parameter model is the best one as BIC has the lowest value. For both more and less complicated models the BIC value grows with the growth of the target function or model complexity. It means that only in the middle layer anisotropy is significant enough to allow one to introduce elliptical anisotropy coefficient. In other layers the impact of possible anisotropy on the traveltimes can be either not present or hidden under measurement errors. Dependence between level of measurements errors and optimal number of parameters chosen with BIC criterion is analysed in [10] with example of synthetic datasets.

Conclusions
The analysis of results of the velocity model optimization with elliptical anisotropy assumption was conducted. The used algorithm worked correctly finding local minima which allow measured and modelled traveltimes to differ not more than measurement precision. Using Nelder-Mead optimization method brings good results but because of its local character it is important to provide a starting model close to a real one. The created tool can be successfully used to find velocity parameters in multi-layered medium with a starting model created on the base of zero-offset VSP or other a priori data.
The multistart analysis simulates using of global optimization in finding velocity parameters simultaneously showing the regions where local minima occur. Comparing target function value with values of specific parameters can lead to conclusion about effectiveness of optimization algorithm in heading to global minimum. Comparison of velocity parameter values in pairs or groups can suggest existence of correlations between parameters, e.g. a0 -b0, a1 -χ1 (indices refer to number of the layer as described in Fig.4). It may result in difficulties with finding the proper parameter value as they can compensate traveltimes in pairs or greater groups.
The BIC analysis demonstrates that in this case the 7 parameter model is the best one. The optimal values are presented in Table 1. The existence of significant anisotropy described with elliptical anisotropy coefficient helps to obtain model fitting better to measured data simultaneously without complicating the model too much. The presented algorithm composed of finding minima of target function for different parameterization and then choosing the optimal parameterization with BIC criterion provides user with good capabilities to solve the problem of velocity model inversion in multi-layered anisotropic medium.