Automatic estimation of traveltime parameters in VTI media using similarity-weighted clustering

Compared with hyperbolic velocity estimation methods, nonhyperbolic methods (such as shifted hyperbola) are better choices for large offsets or vertical transverse isotropy (VTI) media. Since local seismic event slope contains subsurface information, they can be used to estimate zero-offset two-way traveltime and normal moveout velocity. The traditional velocity estimation methods require a great deal of manual work and are also prone to human error. In order to estimate the traveltime parameters for VTI media automatically, in this paper, we propose to use predictive painting and similarity-weighted clustering to obtain traveltime parameters. The predictive painting is used to estimate zero-offset two-way traveltime, and the shifted-hyperbola traveltime equation is used to obtain velocity and anisotropy attributes. We first map local slopes to zero-offset two-way traveltime and moveout-parameters domain and then use similarity-weighted k-means clustering to find the maximum likelihood anisotropy parameters of the main subsurface structures. In order to demonstrate that, we apply the similarity-weighted clustering method to synthetic and field data examples and the results are of higher accuracy when compared to the ones obtained using multiparameter semblance-based method. From estimation error section, it can be seen that the estimation error of multiparameter semblance-based method is about 3–5 times that of the proposed method.


Introduction
Building a subsurface velocity model is one of the most important issues in exploration geophysics. An accurate velocity model is significant for migration and inversion, which can also directly be used as a lithology indicator. Conventional velocity analysis involves several procedures of scanning effective moveout velocity from seismic gathers, computing semblance and creating velocity spectra for later velocity picking (Taner and Koehler 1969;Yilmaz 2001;Luo and Hale 2012;Chen 2018c). These procedures intend to extract the normal moveout (NMO) velocity as a function of the two-way zero-offset traveltime at the selected common midpoint (CMP) locations along the seismic line. However, these algorithms use search methods along trajectories described by the hyperbolic equation (Dix 1955), whose validity is restricted to layered isotropic media with offsetreflector depth ratio smaller than one (Thore et al. 1994). Using the hyperbolic traveltime approximations reduces the velocity estimation accuracy for anisotropy media (Fowler 2003). The nonhyperbolic traveltime approximation with anisotropy parameter as reported in Koshnavaz et al. (2016a, b) can describe the reflection moveouts with a high accuracy for large offsets or vertical transverse isotropy (VTI) media. Several nonhyperbolic traveltime approximations have been proposed for the description of the behaviors of reflection curves in stratified isotropic and anisotropic media, and such traveltimes are based on the shifted hyperbola (Malovichko 1978;de Bazelaire 1988;Castle 1994;Siliqi and Bousquié Edited by Jie Hao 2000), rational moveout approximation (Alkhalifah and Tsvankin 1995;Ursin and Stovas 2006) and the commonreflection-surface (CRS) traveltime for weakly VTI media (Bloot et al. 2018).
The conventional semblance-based techniques for moveout-parameters estimation are time-consuming, and local-slope-based techniques are comparatively faster than the conventional methods. Local slopes of seismic events represent the apparent ray parameters or apparent slownesses in the seismic record. They can be used to estimate seismic velocity automatically because they contain subsurface velocity information. Ottolini (1983) uses the idea decomposing seismic data into a range of local slopes to propose the velocity-independent, time-domain seismic imaging method. This idea is generalized by Fomel (2007a), who develops analytical velocity-independent techniques to transform seismic data from prestack domain to the image domain for different time-domain imaging operators. Velocity-independent time imaging can be extended to anisotropic medium (Burnett and Fomel 2009a, b), since the local slopes also depend on the anisotropic properties of the media. Stovas and Fomel (2016) proposed the use of local slopes as well as curvatures, which is the derivative of two-way traveltime with respect to offset, to obtain kinematic attributes, including normal moveout (NMO) velocity, anisotropy anellipticity coefficient and zero-offset two-way traveltime t 0 for different nonhyperbolic approximations. Khoshnavaz et al. (2016a) proposed to use predictive painting (Fomel 2010) of local slopes to estimate the kinematic attributes, which substitutes the need for curvatures by using the result of predictive painting. The local slopes are also applied to other steps of seismic data processing and imaging (Cooke et al. 2009;Bóna 2011;Khoshnavaz et al. 2016b;Khoshnavaz 2017;Decker et al. 2017;Casasanta and Fomel 2011). The local slopes can be obtained by so-called plane-wave destruction (PWD) (Claerbout 1992;Fomel 2002;Fomel and Stovas 2010).
The k-means clustering originally proposed by Forgy (1965) and MacQueen (1967) is a common algorithm in the field of machine learning (Drineas et al. 2004;Vattani 2011, Little andJones 2011). Clustering is a classic unsupervised learning method but is relatively new in the seismic community. In recent years, clustering algorithm has been applied to microseismic data analysis (Chen 2018a, b;Chen et al. 2019). The clustering algorithm is also used in seismic velocity estimation. Recently, two automatic velocity estimation methods in isotropic media have been proposed under the use of the concepts of accelerated clustering (Zhang and Lu 2016) and similarity-weighted clustering (Liu et al. 2018). In this paper, we combine the automatic velocity estimation method for VTI media described in Khoshnavaz et al. (2016a) with the similarity-weighted k-means clustering described in Liu et al. (2018). As proposed in Zhang and Lu (2016) and Liu et al. (2018), estimations of velocity and effective anisotropic parameters can be estimated as k-means clustering centers of suitable mixed distributions of traveltimes and slopes. Such approach replaces the conventional one based on semblance.
In this paper, we first explain the procedure to obtain velocity and anisotropic attributes derived from the shiftedhyperbola traveltime approximation, which uses the estimated local slopes and zero-offset two-way traveltime. Next, we review the theory of zero-offset two-way traveltime estimation using the predictive painting approach. Then, we describe the method of automatic velocity estimation using k-means clustering algorithm, which calculates the cluster centroids of the velocity attributes and anisotropic parameter. Finally, we use synthetic and field data examples to evaluate the accuracy and efficiency of the proposed method.

Moveout attributes from local slopes by shifted-hyperbola approximation
We consider the hyperbolic normal moveout (NMO) equation commonly used in the industry today, which is a small-offset approximation for a horizontally layered-earth model (Dix 1955). In Eq. (1), t is the two-way traveltime from the source to the reflector and back to the receiver, t 0 is the zero-offset two-way traveltime, and x is the offset from the shot to the receiver.
For VTI media or at large offset, the traveltime is approximated by nonhyperbolic equations. Fomel and Stovas (2010) used five parameters ( a , b , c , and t 0 ) to describe the moveout behavior and proposed the following general form of the moveout approximation: The number of parameters in Eq.
(2) can be further reduced to three by the shifted hyperbola (Malovichko 1978;de Bazelaire 1988;Castle 1994), which is written as, where t , t 0 and v denote two-way traveltime, zero-offset twoway traveltime and the NMO velocity, respectively. Parameter s is related to heterogeneity and anisotropy of seismic velocities. Khoshnavaz et al. (2016a) investigate estimation of kinematic attributes using several approximation of nonhyperbolic moveout. In this paper, we use shifted hyperbola, which is one of the most popular moveout approximations for explaining nonhyperbolic moveouts, to study the estimation of moveout parameters with weighted clustering.
Differentiation of (3) with respect to offset x gives the local slope p: The moveout parameters v and s can be obtained by (Khoshnavaz et al. 2016a): and Equations (5) and (6) show that the moveout parameters v and s are dependent on local slope p and zero-offset traveltime t 0 . Note that the local slope and zero-offset traveltime are functions of x and t , namely p(t, x) and t 0 (t, x) , respectively. If we have p(t, x) and t 0 (t, x) , we can define a mapping from p(t, x), t 0 (t, x) into t 0 , v and t 0 , s by Eqs. (5) and (6). The local slope p(t, x) can be obtained, e.g., applying a plane-wave destruction (PWD) algorithm on each CMP gather. The zero-offset traveltime t 0 (t, x) can be obtained by different methods. Fomel (2007a) proposed a method to obtain t 0 (t, x) by using the second derivative of traveltime.
Predictive painting method is also used to estimate t 0 (t, x) with higher accuracy by Fomel (2010). In this paper, we Velocity, m/s Eta Using Eqs. (4-7), a single group of data points for reflection can be represented as so-called mixture distribution model (Zhang and Lu 2016) in t 0 , v and t 0 , domain. Zhang and Lu (2016) assume the error of one reflection caused by the traveltime assumption and slope estimation algorithms follows a Gaussian distribution. For the case of more than one single reflector, a simple linear superposition of Gaussian components (so-called Gaussian mixture model) can be introduced. Once we have these points t 0 , v and t 0 , , we can use the similarity-weighted clustering (Liu et al. 2018) to obtain the maximum likelihood velocity v and anisotropic parameter of the subsurface structures.

Predictive paintings
The main concept of predictive painting (Fomel 2010) is that each seismic event can be predicted from its neighbors by shifting it along the local event slopes. We assume that local slopes p have been previously estimated, for example, by the plane-wave destructor (PWD) method (Claerbout 1992;Fomel 2002). Then, the prediction of a trace from a distant neighbor can be accomplished by a simple recursion (Fomel 2010): where i,j describes prediction of trace j from trace i . Prediction of a trace consists of shifting the original trace along dominant event slopes. Predictive painting can spread information from a given trace to its neighbor recursively by following the local structure of seismic events. Khoshnavaz et al. (2016a) give an illustrative figure regarding the predictive painting of a seismic event using local slopes (see Fig. 1 in Khoshnavaz's paper), which vividly illustrates the basic principle of predictive painting. For each trace number i , the prediction painting method estimates its corresponding zero-offset traveltime t 0,i following Eq. (8). Note that, along the traces of a same event, the estimated zero-offset traveltimes will be clustered around the correct zero-offset traveltime. The method to predict the zerooffset traveltime can be used for hyperbolic and nonhyperbolic event. We use the shifted-hyperbola traveltime to estimate the zero-offset traveltime t 0 (t, x) for all samples of each trace location x.

Moveout-parameters estimation by similarity-weighted clustering
After the local slope p(t, x) and zero-offset traveltime t 0 (t, x) are, respectively, estimated by PWD and predictive painting, one can obtain the points t 0 , v and t 0 , s from p(t, x), t 0 (t, x) by Eqs. (5) and (6). The distribution density of these points t 0 , v or t 0 , represents the uncertainty of the estimated moveout parameters. The greater the local density, the greater the accuracy of the parameters (Zhang and Lu 2016). The similarity-weighted clustering algorithm (Liu et al. 2018) is used to find the centers of these clustering points t 0 , v or t 0 , . In centroids update step of k-means clustering, the local similarity is used by (Liu et al. 2018) where (iter+1) i is the update mean of the cluster Q iter i , x p is the point that belongs to cluster Q iter i ,w p is the weight of each  data point, p is the local similarity measurement between prestack and reference traces, and T is a parameter of thresholding. Local similarity as a weight implies that those data samples with high signal-to-noise ratio are more reliable for the estimation of slope p(t, x) and zero-offset traveltime t 0 (t, x) . Thus, the reliability of data points t 0 , v or t 0 , mapped by p(t, x), t 0 (t, x) is also related to local similarity. Weighted clustering considers another kind of uncertainty on the basis of clustering density. The clustering center represents the place where the clustering density is high, so the moveout parameters corresponding to clustering centers are more accurate than those of other places. Figure 1 shows the flowchart of automatic estimation of moveout parameters in VTI media using similarity-weighted clustering. We first estimate local slopes by means of the PWD algorithm applied on the CMP gathers. Using the obtained local slopes, we employ the predictive painting method to estimate the zero-offset traveltime of each trace for all samples at that trace. Then, we use Eqs. (4-7) to map the local slopes p(t, x) and zero-offset traveltime t 0 (t, x) from each CMP gather onto the t 0 , v and t 0 , domains. To select the samples that belong to the actual events, we calculate the local similarity (Fomel 2007b;Liu et al. 2018) between the prestack seismic trace and zero-offset trace. A user-selected threshold value is selected to cut the low similarity and irrelevant data points to construct the t 0 , v and t 0 , distributions. Finally, the similarity-weighted k-means clustering algorithm is used to automatically estimate the velocity and anisotropic parameters according to the v and maps. The weighted k-means clustering can get the clustering centers of these points t 0 , v or t 0 , , and the corresponding velocity or anisotropic parameter of these clustering centers is the estimated one. The final velocity or anisotropic-parameter section can be obtained by interpolating (e.g., cubic spline interpolation) along these clustering centers.
Moveout-parameters estimation can be implemented either by semblance-based method (Yilmaz 2001;Luo and Hale 2012) or by oriented local-slope-based method (Fomel 2007a;Khoshnavaz et al. 2016a). In these methods, a moveout-parameter map (semblance-scanned map or oriented-parameter map) is first obtained, and then automatic or manual picking is done to get the final velocity section. The oriented-parameter map has higher resolution than semblance-scanned map because only the true event slopes are identified by the slope estimation algorithm (Fomel 2007a;Khoshnavaz et al. 2016a). On the other hand, the semblance map contains the obvious uncertainty of moveout parameters. The greater the energy in the semblance-scanned spectrum, the more accurate the moveout parameter is. Because the oriented-parameter map from local slopes is obtained by directly mapping the amplitude information (e.g., instantaneous amplitude) of seismic data (Fomel 2007a), the uncertainty in parameter-oriented map is not very clear; that is, the physical meaning of energy on map is not as intuitive as with the semblance-scanned map. In the estimation based on clustering of oriented moveout parameters (Zhang and Lu 2016;Liu et al. 2018), the distribution density of the mapped points t 0 , v or t 0 , provides the uncertainty information of the estimated parameters. In this sense, the proposed method takes the uncertain information into account, and the uncertain information can fully play their role in estimating traveltime parameters by weighted clustering.
In the proposed method, the local slopes and clustering are used together. The local slope establishes the mapping relationship between the curve bending degree of event and the velocity, and the weighted clustering realizes the automatic selection of the most possible velocity. The traditional traveltime semblance method is to pick up the most possible velocity points manually. Therefore, the proposed method is equivalent to the two steps of the semblance-based estimation method for calculating the velocity or spectrum and picking.

Synthetic example
In this section, we use a synthetic layered seismic data (Fig. 2a) to show how to utilize the proposed method to estimate velocity v(t 0 ) and anisotropic parameter (t 0 ) . The CMP gather in this example contains 13 shifted-hyperbola events with the maximum amplitude 0.035. The red and blue curves in Fig. 2a indicate hyperbolic traveltime and shiftedhyperbola traveltime, respectively, which shows an obvious better performance of the latter in the difference in far offset. The blue curves in Fig. 3a and b are the exact velocity v(t 0 ) and anisotropic parameter (t 0 ) used for modeling. Figure 2b shows the local slopes estimated by applying PWD to the CMP gather. The zero-offset traveltimes (Fig. 2c) are estimated by predictive painting with the local slopes (Fig. 2b). Figure 2d shows the flattened CMP gather using predictive painting. We clearly see that the CMP gather is well flattened by using the zero-offset traveltimes estimated by predictive painting. Figure 2e shows the local similarity between prestack and zero-offset trace, which can be used to define a threshold value and weights to be used in the clustering algorithm to select the data point with the highest SNR. The estimated local slopes and zero-offset traveltimes are used to measure the moveout parameter using the shifted-hyperbola approximations. Figure 3 shows the comparison of velocity t 0 , v and anisotropic parameter t 0 , estimated with the proposed and the multiparameter semblance methods, respectively. The black dots are the mapped points t 0 , v and t 0 , from local slopes p(t, x) and zero-offset traveltime t 0 (t, x) , and the red dots indicate the cluster centroids after weighted k-means clustering of the black-dot points. The result of the proposed method (the red curves) is obtained by interpolating the red dots. The yellow lines indicate the automatic picking velocity and anisotropy parameter obtained by the multiparameter semblance. The blue curves indicate the modeled velocity v and anisotropy parameter . One can find that the results from the weighted clustering method (red curves) are closer to the ones of the exact model (blue curves) than the corresponding ones obtained by the multiparameter semblance method. The velocity points t 0 , v (Fig. 3a) are very focused, while the points of anisotropy domain, which are generated by using shifted-hyperbola approximation for Fig. 7a. The yellow curves indicate the automatically picking velocity and anisotropy parameter from multiparameter semblance. The blue curves indicate the exact effective velocity v and anisotropy parameter used for modeling. The black dots are mapped t 0 , v and t 0 , , and the red dots indicate the centroids of clusters after weighted k-means clustering for black dots. The red curve is obtained by interpolating the red dots parameter t 0 , (Fig. 3b) are less focused than that of velocity. The estimation of is more sensitive difficult than velocity. Fortunately, there is little difference between the centroid (red dots) of these black points and the true parameter values (blue curve).
The black dots are mapped t 0 , v and t 0 , , and the red dots indicate the centroids of clusters after weighted k-means clustering for black dots. The red curve is obtained by interpolating the red dots.
To illustrate the uncertainty of moveout-parameter estimation, we show the local density (Fig. 4) of the points t 0 , v and t 0 , indicated by black dots in Fig. 3. Here, we use 2D histogram function as the approximation of local density. White color in Fig. 4 denotes that the points (black dots) are dense. One can see that the trend is similar to the real parameter values. To test the effect of noise on the algorithm, we add Gaussian random noise on CMP gather (Fig. 5a). The variance of noise is 10 −7 . The noise has little effect on velocity estimation, but the points for t 0 , are more dispersed (Fig. 4d) than the noise-free case (Fig. 3b). From the local density of ( Fig. 5f), we can find that the uncertainty of has worsened. The weighted k-means clustering method gives a more acceptable result than the multiparameter semblance method. Figure 6 shows a prestack dataset from the Viking Graben dataset. In order to test our method, we use conventional semblance-based velocity analysis (Taner and Koehler 1969) to flatten the events in CMP gathers.

Field example
In the same way as in Koshnavaz et al. (2016a, b), we inject an anisotropy into the original gather to generate new CMP gathers which contain a known VTI anisotropic effect characterized by the anisotropy parameter = 0.1+0.07 * t 0 + 0.01 * sin 4 * t 0 + 0.01 * sin(0.03 * x) , where t 0 is zero-offset two-way traveltime and x is CMP number. Figure 7a shows a generated CMP gather obtained  Fig. 10 a Exact effective velocity section used for forward modeling. b Estimated velocity section using weighted k-means clustering. c Estimated velocity section using multiparameter semblance analysis. d Is the difference between b and a, and e is the difference between c and a by the inverse NMO-corrected (anisotropic) CMP gather using the previous velocity v and the new anisotropy parameter . The estimated local slopes are shown in Fig. 7b. Figure 8a shows the estimated zero-offset traveltimes computed with local slopes and predictive painting. The flattened CMP gather is shown in Fig. 8b. Figure 9 shows the estimated anisotropy parameter and velocity v using local slopes and zero-offset traveltimes with the shifted-hyperbola approximation. In Fig. 9, the yellow curves indicate the velocity and anisotropy parameter obtained by multiparameter semblance. The blue curves indicate the modeled velocity v and anisotropy parameter . For this real dataset, the focused velocity points t 0 , v are still better, but the focused anisotropy parameter points t 0 , are more dispersed. Compared with multiparameter semblance, the clustering method is closer to the exact model parameters. We obtain the velocity and anisotropy parameters of all CMP gathers by interpolating clustering centers (red dots in Fig. 9). Figures 10 and 11 show the results of velocity and anisotropy parameters, respectively. Similarly, the weighted k-means clustering method gives a more accurate result than the multiparameter semblance analysis. From estimation error sections (Figs. 10d and e, 11d and e), it can be seen that the estimation error of multiparameter semblance-based method is about 3-5 times that of the proposed method.

Conclusions
We applied a method for automatic velocity and anisotropy parameter estimation by predictive painting and similarityweighted clustering. By applying the predictive painting method to CMP gathers, the zero-offset traveltimes are estimated with local slopes. The velocity and anisotropy parameters are derived by means of the shifted-hyperbola Fig. 11 a Exact anisotropy parameter section used for forward modeling. b Estimated anisotropy parameter section using weighted k-means clustering. c Estimated anisotropy parameter section using multiparameter semblance analysis. d Is the difference between b and a, and e is the difference between c and a traveltime approximation, which can describe traveltime curves for VTI media. Using the similarity-weighted k-means clustering algorithm, the proposed approach can simultaneously and automatically build both the velocity and the anisotropy parameter models. The synthetic and field examples demonstrate a better performance of our proposed approach as compared with the multiparameter semblance method.
The computational cost of this method is slightly larger than that of the multiparameter semblance-based method, but the latter method needs manual picking and more manual work. In this paper, we did not consider the case of threedimensional CMP gather. If we want to extend it to threedimensional case, we need further investigation because more involving parameters for 3D gather may make the mapping from local slopes to moveout parameters and the clustering more complex.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.