AU2016202792A1 - Automatic NMO correction and Full Common Depth Point NMO Velocity Field Estimation In Anisotropic and lateral heterogeneous Media - Google Patents

Abstract

Although several anisotropic traveltime equations have been written by different researchers to perform normal moveout (NMO) correction of common depth point (CDP) seismic gathers, none has been able to address the problems associated with the far offset data. In this paper, we present a computational process that, not only accurately flattens and NMO-correct the far offsets data to zero-offset times no matter how far, but simultaneously provide full NMO velocity (Vnmo) fields and corresponding anisotropy field. The method is based on doing a predefined number of NMO velocity iterations using linear vertical interpolation of different NMO velocities at each seismic trace individually. At each iteration, we measure the semblance between the zero offset trace and the next seismic trace, then after all the iterations are done, the one with the maximum semblance value is chosen, which is assumed to be the most suitable NMO velocity trace that accurately flattens seismic reflection events. The other traces follow the same process, and a final velocity field is then extracted. Isotropic, anisotropic and lateral heterogenous synthetic geological models were built based on a ray tracing algorithm to test the method. A range of synthetic background noise, ranging from 10 to 30 %, was applied to the models. In addition, the method was tested on Hess's VTI (vertical transverse isotropy) model and Alaska (gas field) pre-stack seismic CDP field example. The results from the presented examples show an excellent NMO correction and extracted a reasonably accurate NMO velocity field.

Description

1. Introduction

One of the major concerns in seismic data processing, is the seismic anisotropy. Seismic anisotropy is a seismological term that describes how Seismic wave velocity changes with distances or angles. A medium is called anisotropy if velocity at a fixed spatial location varies with the direction of propagation. In typical subsurface formations, velocity changes with both spatial position and propagation direction, which makes the medium heterogeneous and anisotropic (Tsvankin, 1998). In contrast to such mediums, (Alkhalifah, 1997) assumed the absence of the lateral change of velocity for the anisotropy calculation in VTI media, the same equations are used for HTI media (just a rotation of the VTI gives the HTI type of anisotropy). We present a new method that takes the lateral change of velocity in consideration as will be discussed later.

Anisotropy in sedimentary sequences is caused by the following main factors (Thomsen, 1986): • Intrinsic anisotropy due to preferred orientation of anisotropic mineral grains or the shapes of isotropic minerals. • Thin bedding of isotropic layers on a small scale compared to the wavelength (tilted or horizontal layers). • Vertical or dipping fractures or microcracks.

The combination of the three types is common to be produced by these factors. If transversely layers are dipping, the TI medium has a titled symmetry axis and becomes azimuthally anisotropic. (Grechka and Tsvankin, 1998), showed that the azimuthal variation of NMO velocity is an ellipse in the horizontal plane, irrespective of the type or degree of intrinsic anisotropy. For VTI anisotropy, the NMO ellipse will degenerate into a circle, with no azimuthal variation. For titled (TTI or HTI) the azimuthal variation in NMO velocity will be described by an ellipse. This type of azimuthal anisotropy has a profound effect on the seismic data (Williams and Jenner, 2002); (Lynn, 2004) and can be attributed to stress loading ( (Rasolo-fosaon, 1998)) or complex fracturing (Bakulin et al., 2000), see also (Tsvankin et al., 2010). A single set of vertically oriented fractures (open and fluid filled) can result in an azimuthally anisotropic medium in which the velocities are faster parallel to the fractures (V fast) than perpendicular to them (V slow).

Still, the anisotropic contribution to the NMO velocity seems to be the most pervasive reason for errors in converting time to depth in many exploration areas (Alkhalifah et al., 1996). In coal seam gas exploration, the problem of anisotropy is complex, where it exhibits VTI anisotropic type by the vertically aligned thin layers, HTI and TTI by the horizontal and titled aligned fractures with the coal seam beds

It is very difficult to know the presence of anisotropy in the subsurface geological formations only from P-wave seismic data and special analysis is required for this. A number of approaches to anisotropic velocity analysis and traveltime inversion have been developed (Tsvankin and Thomsen, 1995); (Alkhalifah and Tsvankin, 1995); (Baan and Kendall, 2002). Most of these approaches analyze higher-order terms in the NMO equation in order to improve images obtained in time processing. Such traditional non-hyperbolic NMO velocity analysis is often used for building initial velocity/anisotropy models for prestack depth migration. Then the interval velocities required by depth imaging are obtained from these initial models by a layer-stripping process using the Dix’s formula, (Dix, 1955). Dix inversion formula built a direct relation between the NMO velocity and the subsurface layers interval velocity.

The problem with the Dix’s formula is the error accumulation with time (Zhou et al., 2003), other methods were developed to overcome such a problem like tomography and full waveform inversion. In coal seam gas case, the target reservoir is very shallow compared to the shale unconventional reservoirs. The application of the Dix inversion formula in coal seam should give reasonable velocity information at the reservoir level due to the short time interval of the coal layers. NMO correction in isotropic media is very simple, and can be achieved using the hyperbolic travel time equation ( 1), (Sheriff, 1991), defined the normal moveout effect as the variation of reflection arrival time because of source-receiver offset.

In anisotropic case, anisotropy causes two major distortions to the normal moveout in seismic reflection data. First, in contrast to isotropic media, normal moveout differs from the vertical velocity, and the second is a substantial increase of deviations in hyperbolic moveout in an anisotropic layer, (Behera and Khare, 2010). Hence, with the help of conventional velocity analysis of short-spread moveout(stacking)velocities doesn’t provide enough information to determine the true vertical velocity in a transversely isotropic media. Therefore, it becomes essential to estimate the anisotropic effect from the long offset P-wave seismic data. The long offset P-wave seismic data, has observable non-hyperbolic moveout, which depends on the anisotropic effect and the normal moveout velocity (Vnmo).

To correct the pre-stack seismic reflection data due to the normal moveout effect (NMO) in anisotropic media, the anisotropic parameter (η) has to be calculated, according to Alkhalifah-Tsvankin equations for anisotropic media. The velocity information (e.g. horizontal and vertical velocity) that is required to calculate the anisotropic parameters, mostly, needs velocity logs data. Some other techniques estimate the horizontal velocity by scanning the Thomsens anisotropic parameter (η) with the NMO velocity, (Alkhalifah and Tsvankin, 1995).

Scanning the Thomsens anisotropic parameters with NMO velocity estimated from the semblance plot, requires high quality seismic data with long offsets, and offset-depth ratio (ODR) greater than 2.0. Otherwise the estimation of the horizontal velocity will not be stable, (Jenner, 2011).

Improvements have been done on the conventional normal moveout correction by many researchers. Taner and Koehler (1969),A1-Chalabi (1973) and M. Gidlow (1990) applied corrections using an order higher than 2. De-Bazlaire (1988) introduced the shifted hyperbola correction. Rupert and Chun (1975) presented the block move sum (BMS) concept, where a series of static shifts are applied to blocks of seismic reflection data followed by summation. Calderon-Macias C. (1988) proposed an automatic NMO correction by a feed forward neural network (FNN). Perroud H. and Tygel (2004) introduced an automatic non-stretch NMO. Abbad et al. (2008) proposed an automatic non-hyperbolic velocity analysis method, the method is based on efficient parameterization and automatic event picking to produce a zero-offset section and two attribute maps. Kazemi and Siahkoohi (2012) introduced a stretch free NMO correction method called local stretch zeroing (LSZ). It improves the NMO analysis through an optimized selection of the mute zone and avoiding interpolating samples.

In our work, we introduce a new method for NMO correction and non-hyperbolic velocity field estimation at each offset. At a given common depth point (CDP) the NMO velocity is constructed as a temporal interpolation of an offset dependent, NMO velocity and offset independent final velocity (this term will be discussed later). NMO velocities at individual offsets are scanned to match shifted traces at offsets. We call this new approach Full CDP Velocity (F-CDP-V) correction. The fit of the shifted traces against the zero offset trace was measured using a new semblance measure. In contrast to conventional semblance measures ours correlates the entire traces at each offset rather than inspecting correlation only at each sample time. The new semblance measure is more suitable for F-CDP-V correction as it allows to work on individual offsets. Our tests have also shown that the new semblance measure is suitable for a stable automatic NMO velocity picking as shown by the tested data sets. It is pointed out that our new approach does not require any prior information except for the final velocity, which for instance can be obtained using the conventional semblance analysis. This provides a big advantage in comparison to other conventional methods which typically require estimates of some of the anisotropy parameters (Like assuming the type of the anisotropy whether it is VTI, HTI or TTI before starting the processing workflow unlike our method which doesnt require a prior assumption of the anisotropy). In practical applications the F-CDP-V correction has provided results that are close to be stretching free at distant offsets and does not require to apply a mute for the far offset stretching.

In the following section 2 we will introduce our new NMO correction method. Section 3 will show results from synthetic data sets for an anisotropic media, an isotropic media with a thin layer and lateral heterogeneous media. An application to two field data sets is presented in section 4. First example, is a simple synthetic data set simulating a single layer in weak VTI media. The second data set, is also synthetic to simulate a lateral heterogenous media with lateral variation of P wave velocity and density. We also show the results using the well known Hess model representing a weak VTI anisotropy. Finally, we applied our method on a real data sets from Alaska, USA and one CDP example from a coal seam gas reservoir in Australia. In the final section we will draw some conclusions and discuss further work. 2. Methodology

Given is a set of seismic traces q for a common depth point (CDP) gathers. The value of the trace for sample k at offset x is denoted by q(k, x). To simplify the presentation we assume that the same sample interval s > 0 is used for all traces. As commonly done we assume that the primary signal amplitudes are higher than the background noise. The starting point of our new method is the hyperbolic travel time equation

(1) where tx is the arrival time at offset x for a signal observed at time to at offset zero. The stacking velocity is denoted vstack which in case of a layered media is approximately given by a weighted RMS velocity at least for small offsets (Dix, 1955). For the anisotropy the travel time equation is extended typically through higher order terms. Alkhalifah and Tsvankin (1995) suggested the following model which is widely used in practice:

(2) where vnmo is the vertical NMO velocity and η is an anitropy parameter. This equation can be written in the form 1 with

(3)

For shallow signals (t0vnmo « x) we have vstack « vnmo(l + 2η)~ζ. For long travel times (tovnmo >> x) t the stacking velocity vstack approaches vnmo. This observation motivates the generic velocity model

(4)

At each offset x the velocity is defined as a linear, temporal interpolation of the final value of velocity ^ set at time t — T and the NMO velocity v(x) which varies with offset x. v(x) has both initial and end value we called it ml and m2 to control the iteration process. Typically T is the time window length. The value for νχ is chosen constant. Using this velocity we can correct the arrival times of the signal q(k, x) at offset x. Following travel time equation 1 the new arrival time t(k, x) of sample k is given by the following equation

(5) where s is the sample interval. Using this formula we can identify for each sample k in the zero offset the corresponding sample k at offset x leading to a shifted and possibly stretched signal q as

(6) where k denotes the nearest integer round of It is pointed out that the sequence k is non-decreasing; however consecutive indices may differ by more than one which leads to a stretching in the shifted signal and can also be identical. The approximation of the sample number sometimes may shift between 1 or more, which causing a phase shift at some offsets but in practice this is not a problem.

The objective is now to find the best stacking velocity vstack which gives the best match of the time shifted traces q(.,x) at offsets x to the trace q(., 0) at offset zero. For this optimization process we need to define a quality measure (also called semblance measure). In our method we use the following measure: (7)

In fact is can be shown (e.g. see (Solomentsev, 2001)) that at any fixed offset x k(x) takes values between 0 and 1 where the value one is taken only if the traces at offset x and zero offset differ by a (offset x dependent) constant factor

Cx: (8)

It is our objective to chose the NMO velocity such that k(x) takes the value (or at least the largest possible value) at all offsets. This definition differs from conventional semblance (Taner and Koehler, 1969) which - in contrast to our semblance measure - is defined per sample k:

(9) where Nk gives the number of non-zero values across the offsets (Smoothing across samples is omitted). Using arguments similar to proof for κ one can show that a{k) is between 0 and 1 and for fixed sample k the value 1 is taking only if there is a factor (7¾ (varying with the sample) such that

(10) A comparison of this equation with equation (8) illustrates the difference objectives when maximizing the corresponding semblance function. While using the conventional semblance measure σ signal shifting is targeting to obtain the same value across all offsets for each sample while the new semblance measure is trying to align the entire trace with the zero offset trace and allows for re-scaling of the signal.

To find the best velocity profile the NMO velocity v(x) at given offset x is chosen between minimum and maximum values with an fixed increment. For each value of v{x) the trace is shifted using a temporal, linear interpolation between v(x) and the final velocity v<». Then the value for v{x) is chosen which provides the largest value for the semblance n(x) is calculated. Typically, values for v(x) are chosen between 100m/s and 2100m/s with an increment of 10m/s which requires 200 evaluations of the semblance measure. This scanning process is illustrated in Figure 1.

To quantify anisotropy we inspect the change of the NMO velocity with the arrival time and offset. Taking the derivative of (1) with respect to offset x we obtain (11)

As for isotropic media stacking velocity v is independent from offset, the second term is attributed to the anisotropy of the media. The magnitude of this term relative to the first term defines a measure for degree of anisotropy which is function of vertical travel time and offset:

(12)

The distribution of the anisotropy T at the offset Xi can be easily calculated from the covered velocity model 4 using

(13) (14)

For the two extreme cases of narrow (x — 0) and far (x — oo) the anisotropy indicator T takes the value zero. The maximum value Tmax(to) value across all offsets is taken at offset x — tovnmo with value

(15) (16) (17) (18)

An automatic linear upper muting can also be applied using the following equation: (19)

Where tstart is the start point that is defined as the initial time for muting. ten<i is the endpoint that is defined as the end time value. T is the total time window of the CDP. Trn is the trace number (n). 3. Synthetic Data 3.0.1. Isotropic media with multiples and noise: A synthetic CDP gather consistent of eighteen layer was created using Ricker wavelet with a band width of frequency from 5 to 50 Hz and the velocities range between 1500 to 2650 m/s. In order to test the robustness of the F-CDP-V method we have added 20% noise to the signal. Two layers were set up with a lower velocity than the layers below and above to simulate multiples. In Figure 2, these two layers can be noticed at the interval between 3 and 4 seconds. A linear automatic muting is applied using tstart — 0.2s and tend — 1.2s. The F-CDP-V correction is applied scanning a shallow velocity range from 900 to 2100 m/s, an increment 10 m/s and the final velocity is 4000 m/s. The resulting corrected traces are shown in right image in Figure 3.

Our automatic method assumes that the velocity is increasing with depth so automatically, it will not correct any layer with a lower velocity than the above layers. This can be considered both an advantage and a drawback of the method. According to Sheriff (1991), he mentioned that seismic signals associated with high velocities tend to represent primarily seismic reflections while those associated with lower velocities tend to represent multiples. So the method automatically ignores layers with lower velocities, which can be considered as the unwanted seismic multiples. In case of the lower velocity and the defined layers are not multiples; a solution can be done by specifying an interval so the method can run separately on this interval and merge the results back to the original matrix.

Figure 4a shows the new semblance measure k(x) as function offset x and values for NMO velocity. The blue dotted line indicates the selected NMO velocity as a function of offset with minimum semblance measure.

The recovered velocity field is shown in Figure 4b. The F-CDP-V method assumes that the velocity is increasing with depth so automatically, it will not correct any layer with a lower velocity than the layers above. This can be considered both an advantage and a drawback of the method. (Sheriff, 1991) mentioned that seismic signals associated with high velocities tend to represent primarily seismic reflections while those associated with lower velocities tend to represent multiples. So the method automatically ignores layers with lower velocities, which can be considered as unwanted seismic reflection events. 3.0.2. Anisotropic media:

To test and validate the F-CDP-V method in anisotropic media, we used a very simple one layer model using Tsvankin’s equation(2) to construct the synthetic non-hyperbolic seismic reflections, see Figure 5.

The following parameters were used:

Table 1: Synthetic anisotropic media: The parameters used to create the non-hyperbolic synthetic seismogram.

Again we use a Ricker wavelet with a band pass filter with these pa-rameters:0,5,50,60 Hz. The dominant frequency should be around 25 Hz. The vertical NMO velocity is set to vnmo — 1500m/s and horizonatal is choosen to Vh — 1800m/s giving an anisotropy parameter of η — --1) — 0.1.

The reflector is placed at 0.5.s. The F-CDP-V correction is applied by scanning NMO velocity range form 10 to 2500 m/s in increments 10m/s with the finial velocity of 2100to/s). The values of the semblance value n{x) as defined in equation (7) is shown in Figure 6. The vertical axis represents the scanned values for the NMO velocity v{x) while the horizontal axis represents offset. The actually picked NMO velocity with maximum semblance value is represented through the black, dashed line. At a given offset the values of the semblance for a concave shape a feature that makes it easy to pick a maximum value without the problem of local maximums. However the values near the maximum are typically form a plateau which makes picking the maximum values unstable. This can lead to unwanted oscillations in the selected NMO velocities indicated by the jagged black dashed line in Figure 6. To avoid this problem we apply simple smoothing to the semblance value before picking the maximum value. The smoothing has the effect that its creates distinguished maximum over the plateau. As a result the picked near surface velocities are smoother as function of the offset, see red line in Figure 6. The picked and interpolated velocity as function of offset and time is shown in Figure 7. The traces after NMO correction are shown in Figure 8. Notice that there is no stretching at the far offsets.

From the constructed velocity field we recover vnmo, the horizontal velocity Vh and the anisotropy parameter η. The results listed as follow in Table 2. The F-CDP-V correction method is able to reproduce the parameters of Tsvankin’s NMO equation (2) well Picking based on the smoothed semblance gives a slightly better match. Based on (Jenner, 2011)observation about how the increasing of the ODR can improve the stability and the accuracy of measuring NMO velocity, our semblance plot is in a perfect match with that observation, showing that at the far offsets the path of the automatic tracking becomes narrower and the uncertainty area becomes smaller with offset.

Table 2: Synthetic anisotropic media: The calculated anisotropic parameter η from the far offset velocity extracted from F-CDP-V velocity field. 3.0.3. lateral heterogeneous media:

The third example is testing the F-CDP-V method for a laterally het-erogeous medium. The example covers a lateral extend of 6km to depth of 2km and is inspired by a geological set-up as it typically found for coal seam gas (CSG) application (McConachie et al., 2015). The synthetic model that is built to test our method is to simulate the case of CSG layers, which are very thin and the gas is trapped within open fractures, see Figure 9. The media can be considered as anisotropic case where it exhibits VTI properties due to the vertical thin layers, HTI and TTI by the horizontal and titlted aligned fractures within the CSG beds. To simulate a similar case in a synthetic model, a thin layered model was built and filled with different velocity and density values along its length. The very thin layer near the center of the synthetic model, assumed to be a simulation to the thin CSG layer, and the variable velocities and densities are assumed to be the effect of the presented cleats and fluids, and it is shown as a slight variation in the color with in the layers in Figure 9. The control points that have been marked in Figure 9 by the asterisks are located at distances of 2, 4, and 6km from the left face of the domain. The slowness (4j·) and density values are laterally interpolated between the control points. The corresponding values are shown in Tables 3 and 4. The low-density values in layer 2 is representing the coal layer.

Table 3: Synthetic lateral heterogeneous model: The values representing the slowness2 at each controlling point.

Table 4: Synthetic lateral heterogeneous model: The values representing the density at each controlling point. A ray tracing algorithm is used based on Gaussian beam synthetic seismograms, The Gaussian beam method is an asymptotic method for the computation of wave fields in smoothly varying in-homogeneous media, and was proposed by (POPOV, 1981) and (POPOV, 1982) based on an earlier work of (BABICH and PANKRATOVA, 1973). The method was first applied by (POPOV and CERVENY, 1980), (KATCHALOV and POPOV, 1981) and (POPOV, 1981) and (CERVENY and PSENCIK, 1982) to describe high-frequency seismic wave fields by the summation of the par-axial Gaussian beams. One of the advantages of the method is that the individual Gaussian beams have no singularities either at caustics in the spatial domain or at pseudo-caustics in the wave number domain, Figure 10.

That method is used to create the synthetic seismic data, Figure 11. The total number of shots is 40 with a total of 60 traces at each CDP gather. The ray traces’ angle ranges from —75° to 75° with 0.004 s sampling interval. Figure 11 showing the whole seismic profile and how the amplitude changing along the model due to the anisotropic effect.

The extracted synthetic seismic profile in Figure 11 is a shot gather, and represents only one shot from the 40 shots. It is needed to be sorted into a CDP gather. Seismic Un*x (Stockwell, 2011) is used to sort the synthetic seismic profile. CDP number 1000 has been selected to be used (Figure 12). The selection criterion is based upon the number of traces and coverage area.

Conventional semblance analysis is tested to find the velocity information of the synthetic CDP, Figure 13. Where a range of Vnmo is applied from 1000 to 4000 m/s. Also the extracted velocity information (Vnmo) along with the time values (Tnmo) is listed in Table 5.

Table 5: Synthetic lateral heterogeneous model: Showing the picked arrival times and the corresponding NMO velocities using the conventional semblance velocity analysis.

As the semblance method only corrects the hyperbolic portion of the moveout curve. The corrected part (hyperbolic) most likely is the near offset, which hardly represent the true velocity of the formation. One can assume that the far offset traces can better represent the formation velocity as it cover more area than the near offset ones, so it can be used as a Vrms to simplify the interpolation or smoothness process between the CDPs.

The curved shape of the corrected seismic reflectors at Figure 13 (the plot at the right side) is because of the effect of the anisotropy that causes the non-hyperbolic shape.

Velocity field extraction and anisotropic effect calculation in lateral heterogeneous media:

The previously mentioned synthetic CDP 1000 is used to test the method. The values for m\ and m2 are 700 and 2100 m/s respectively, with 10 m/s as an increment step and the final velocity is 1100 m/s. The final output velocity field is shown by Figure 15. Where Figure 14, showing the automatic NMO velocity tracking and how our formula equation: 7 improves the resolution and the tracking of the NMO velocity at the far offsets, the black dashed hne showing the default tracking while the red one showing the smoothed/averaged tracking.

The velocity field from CDP 1000, showing an obvious NMO velocity variation with offset unlike the isotropic case where the NMO velocities were homogeneous. In this case, the anisotropic effect can be taken into consideration, and be calculated as it will be shown in the next section.

Figure 16, showing a combined picture of the original synthetic, the corrected CDP and the stacking output.

The stacking output at Figure 16, showing a very good result if it comes to the amplitude value and the accuracy of detecting the four layers compared to the conventional method. Figure 17, showing the stacking output from the conventional velocity analysis. 3.0-4- Hess model CDP example:

Hess model is a synthetic seismic model with a weak VTI media, built by the Colorado school of mines. The following figure showing a CDP number 500 after an automatic gain control is applied, Figure 18.

The extracted velocity variation with offset, Figure 19 shows that the Hess model has a very weak VTI media as the velocity seems homogeneous, and it is hard to tell that there is anisotropic effect. 4. Field data: 4-0.5. Alaska Pre-stack seismic CDP field data example:

Alaska marine 16-81 2D Pre-stack seismic line, G.Keys and J.Foster (1981) is used to test the method. Alaska seismic data is open for download, visiting U.S. Geological Survey. The data are sorted by CDP gather, and dip filter is applied.Figure 20, showing the pre-stack seismic data before and after the Automatic NMO correction. The following parameters were used for the correction: 1- Minimum initial velocity value = 5000 m/s. 2- Maximum initial velocity value = 15000 m/s. 3- Final velocity value = 20000 m/s.

Additionally, the conventional velocity analysis plot was used only to compare the range of the extracted NMO velocities from our method and the NMO linear range from the semblance plot, Figures 21, 22 showing that both the conventional manual method and our automatic method, have the same NMO velocity range which made the estimated full CDP velocity field more reliable. Figure 22 is the estimated NMO velocity field. 5. Conclusion:

This paper introduces a new method to automatically NMO correct Pre-stack seismic reflection events in a CDP gather in different medias (e.g. Isotropic, anisotropic and lateral heterogenous). The method can simply work without a need to predefine the anisotropic parameters. It provides a better far offset flattening with minimum stretching effect. Moreover, it provides a dense NMO velocity information, which is not available using the current techniques. The method automatically ignores the multiples, unlike the conventional method, where extra seismic processing steps need to be done to get rid of the unwanted multiple signals. Furthermore, It introduces a possibility to detect and interpret the lateral changes in NMO velocity (e.g. due to gas, fluids, fractures, faces change, etc.) as the extracted high dense velocity field covers all the possible azimuths and offsets. 6. Acknowledgement:

Thanks to US geological survey for providing the pre-stack Alaska seismic data online for free.

References

Abbad, B., Ursin, B., Rappin, D., 2008. Automatic non-hyperbolic velocity analysis. GEOPHYSICS 74, U1-U12.

Al-Chalabi, M., 1973. Series approximation in Velocity and travel time computations. Geophs. Prosp 21, 783-795.

Alkhalifah, T., 1997. Velocity analysis using non-hyperbolic moveout in transversely isotropic media. Geophysics .

Alkhalifah, T., Tsvankin, I., 1995. Velocity analysis for transversely isotropic media. Geophysics 60, 1550-1566.

Alkhalifah, T., Tsvankin, 1., Larner, K., Toldi, J., 1996. Velocity analysis and imagine in transversely isotropic media: Methodology and a case study. The leading edge 5, 371-378.

Baan, M., Kendall, J., 2002. Estimating anisotropy parameters and traveltimes in the -p domain. Geophysics 67, 1076-1086. BABICH, V.M., PANKRATOVA, T., 1973. On discontinuities of Green’s function of the wave equation with variable coefficient. Numer. Math. 6, 9-27.

Bakulin, A., Grechka, V., Tsvankin, I., 2000. Estimation of fracture parameters from reflection seismic data—Part I: HTI model due to a single fracture set. GEOPHYSICS 65, 1788-1802.

Behera, L., Khare, P., 2010. Improvement of stacking image by anisotropic velocity analysis using p-wave seismic data. National Geophysical Reasearch Institute (NGRI). 8th Biennial International Conference and Exposition of petroleum geophysics .

Calderon-Macias C., S.K.M.S.L.P., 1988. Automatic NMO correction and velocity estimation by feed-forward neural network. Geophysics 63, 1696-1707. CERVENY, V.P.M.M., PSENCIK, I., 1982. Computation of Wave Fields in Inhomogeneous Media- Gaussian Beam Approach. Geophysics 70, 109-128.

De-Bazlaire, E., 1988. Normal Moveout revisted Inhomogeneous media and Curved interfaces. Geophysics 53, 143-157.

Dix, C.H., 1955. Seismic velocities from surface measurements. Geophysics 20, 68-86. doi:10.1190/1.1438126. G.Keys, R., J.Foster, D., 1981. A Data set for Evaluating and comparing seismic inversion methods. Mobile Exploration and Producing Technical Center,

Dallas, Texas .

Grechka, V., Tsvankin, I., 1998. 3-d description of normal moveout in anisotropic media. Geophysics 63, 1079-1092.

Jenner, E., 2011. Combining VTI and HTI anisotropy in prestack time migration. Workflow and data examples. Leading edge . KATCHALOV, A.P., POPOV, M.M., 1981. Application of the Method of Summation of Gaussian Beams for Calculation of High-frequency Wave Fields. Sov. Phys 26, 604-606.

Kazemi, N., Siahkoohi, H.R., 2012. Local stretch zeroing NMO correction. Geophysical Journal International 188, 123-130.

Lynn, H., 2004. The winds of change: Anisotropic rocks their preferred direction of fluid flow and their associated seismic signatures. The fall 2004, SEG/AAPG Distinguished lecturer . M. Gidlow, J.F., 1990. Preserving far offset seismic data using non hyperbolic move-out corrections. SEG meeting, San Francisco, USA, Extended Abstracts 60, 1726-1729.

McConachie, B., Stanmore, P., Creech, M., Hodgson, L.M., Kushkarina, A., Lewis, E., 2015. Reserves estimation and influences on coal seam gas productivity in eastern australian basins. AAPG Geosciences Technology Workshop

Perroud H., Tygel, M., 2004. Non-stretch NMO. Geophysics 69, 599-607. POPOV, M.M., 1981. A New Method of Computing Wave Fields in the High-frequency Approximation. Zapiski Nauchnykh Seminarov Leningradskogo Ot-deleniya Matematicheskogo Instituta im. V. A.Steklova AN SSSR 104, 195-216. POPOV, M.M., 1982. A New Method of Computation of Wave Fields Using Gaussian Beams. Wave Motion 4, 85-97. POPOV, M.M.P.I., CERVENY, V., 1980. Uniform Ray Asymptotics for Seismic Wave Fields in Laterally Inhomogeneous Media (Abstract). Hungarian Geophysical Society 1, 143.

Rasolofosaon, P., 1998. STRESS-INDUCED SEISMIC ANISOTROPY REVISITED. REVUE DE L’lNSTITUT FRANQAIS DU PETROLE 53, 679-692.

Rupert, G.B., Chun, J.H., 1975. The block move sum normal move-out correction. Geophysics 1, 17-24.

Sheriff, R., 1991. Encyclopedic Dictionary of Exploration Geophysics. Society for Exploration Geophysicists, Tulsa, OK 1, 269-271.

Solomentsev, E., 2001. Cauchy inequality, in: Hazewinkel, M. (Ed.), Encyclopedia of Mathematics. Springer.

Stockwell, J.J.W., 2011. Seismic Un*x release No.8. Center for Wave Phenomena, Colorado School of Mines 8, 269-271.

Taner, M.T., Koehler, F., 1969. Velocity spectra-digital computer derivation and application of velocity functions. Geophysics 34, 859-881.

Thomsen, L., 1986. Weak elastic anisotropy. Geophysics 51, 1954-1966.

Tsvankin, I., 1998. Influence of seismic anisotropy on velocity analysis and depth imaging 10.

Tsvankin, I., Gaiser, J., Grechka, V., van der Baan, M., Thomsen, L., 2010. Seismic anisotropy in exploration and reservoir characterization: An overview. GEOPHYSICS 75, 75A15-75A29. doi:10.1190/1.3481775.

Tsvankin, I., Thomsen, L., 1995. Inversion of reflection traveltimes for transverse isotropy. Geophysics 60, 1095-1107.

Williams, M., Jenner, E., 2002. Interpreting seismic data in the presence of azimuthal anisotropy; or azimuthal anisotropy in the presence of the seismic interpretation. The leading Edge 21, 771-774.

Zhou, H., D., P., S., G., B., W., 2003. 3-d tomographic velocity analysis in transversely isotropic media. 73rd Annual, International Meeting, SEG, Expanded abstracts , 650-654doi:10.1190/1.1438126.

Claims (4)

1. Claims

   1. The new pre-stack seismic data processing method can automatically solve both the anisotropic and lateral heterogeneous problems in pre-stack seismic data without a need to predefine the anisotropic parameter to solve for the NMO (normal moveout) problem .
2. 2. The new method provides a full NMO velocity field that covers all the available offsets and arrival times.
3. 3. The new method also provides a calculation method of the anisotropic paramater.
4. 4. The method has been tested on Isotropic, anisotropic, lateral heterogenous media and real field data sets (Alaska pre-stack seismic data)