Migration velocity analysis in factorized VTI media Debashish Sarkar and Ilya Tsvankin Department of Geophysics, Colorado School of Mines Summary One of the main challenges in anisotropic velocity analysis and imaging is reliable estimation of both velocity gradients and anisotropic parameters from reflection data. Approximating the subsurface by a factorized VTI (transversely isotropic with a vertical symmetry axis) medium provides a convenient way of building vertically and laterally heterogeneous anisotropic models for prestack depth migration. The algorithm for P-wave migration velocity analysis (MVA) introduced here is designed for 2-D models composed of factorized VTI layers or blocks with constant vertical and lateral gradients in the vertical velocity VP 0 . The MVA is implemented as an iterative two-step procedure that includes prestack depth migration (imaging step) followed by an update of the medium parameters (velocity-analysis step). The residual moveout of the migrated events, which is minimized during the velocity analysis, is described by a nonhyperbolic equation governed by two moveout parameters determined from semblance analysis. To build a model for prestack depth migration based on the MVA results, the vertical velocity VP 0 has to be either specified at a single point in each block or assumed to be continuous across layer boundaries. A synthetic test for a layered model with a syncline structure confirms the accuracy of our MVA algorithm and illustrates the influence of errors in the vertical velocity on the image quality. Currently this method is being applied to offshore data from West Africa.
Introduction Most existing velocity-analysis methods for VTI media approximate the subsurface with homogeneous or vertically heterogeneous layers or blocks (e.g., Alkhalifah and Tsvankin, 1995; Le Stunff and Jeannot, 1998; Tsvankin, 2001; Grechka et al., 2002). Anisotropic layers, however, are often characterized by non-negligible lateral velocity gradients that may distort the shape of underlying reflectors and cause errors in the anisotropic parameters. Here, we consider a model composed of factorized VTI blocks or layers separated by plane or irregular interfaces. The vertical P-wave velocity VP 0 is assumed to vary linearly within each block, so the vertical (kz ) and lateral (kx ) gradients in VP 0 are constant. The kinematics of Pwave propagation in each block can be described by five parameters: the velocity VP 0 defined at a certain spatial location, the gradients kz and kx , and the anisotropic parameters and δ. Although it is possible to introduce jumps in velocity across the boundaries of the blocks, this
model can be conveniently used to generate smooth velocity fields required by many migration algorithms (in particular, those based on ray tracing). The parameter-estimation methodology employed here is based on the results of Sarkar and Tsvankin (2002), hereafter referred to as Paper I. Although depth imaging of P-wave data requires knowledge of the five parameters listed above (VP 0 , kz , kx , , and δ) in each block, Paper I shows that the moveout of events in image gathers is governed by the following four combinations of these parameters: (1) the NMO velocity at a certain point on√the surface of the factorized layer or block: Vnmo ≡ VP 0 1 + 2δ; (2) the vertical velocity gradient kz ; (3) the lateral velocity gradient combined with the pa√ ˆx = kx 1 + 2δ; rameter δ: k (4) the anellipticity parameter η ≡ ( − δ)/(1 + 2δ). If prestack migration is performed with the correct values ˆx , and η, the image gathers for reflections of Vnmo , kz , k from both horizontal and dipping interfaces are flat. Below we describe our MVA methodology designed to estimate the parameters of factorized VTI media from Pwave reflection data and demonstrate its robustness in reconstructing the spatially varying anisotropic velocity field.
Parameter estimation in a factorized layer In a factorized VTI layer, the moveout of a horizontal event in an image gather is governed by the effective values of the NMO velocity and the parameter η (Paper I): e2kz t0 − 1 , 2kz t0 1 (1 + 8η) (e2kz t0 − 1) kz t0 ηˆ(x, t0 ) = +1 ; 8 2(ekz t0 − 1)2 2 (x, t0 ) = V 2 (x)(1 + 2δ) vnmo
(1) (2)
V (x) ≡ VP 0 + kx x is the vertical P-wave velocity at the surface, and t0 ≡ t0 (x, z) is the zero-offset time at location x from a horizontal reflector at depth z. If long-offset data needed to constrain ηˆ have been acquired, moveout analysis of a single event can yield estimates of both vnmo (x, t0 ) and ηˆ(x, t0 ). The ratio of the NMO velocities [equation (1)] from two horizontal reflectors sufficiently separated in depth (vnmo,1 and vnmo,2 ) can be used to find kz : 2 vnmo,1 (x, t0,1 ) t0,2 (e2kz t0,1 − 1) , = 2 vnmo,2 (x, t0,2 ) t0,1 (e2kz t0,2 − 1)
(3)
where t0,1 and t0,2 are the zero-offset times for the two events.
Migration velocity analysis and anisotropy
Algorithm for migration velocity analysis We apply anisotropic prestack depth migration (the migration algorithm is described in detail in Paper I) and tomographic velocity update to P-wave data acquired over the subsurface composed of factorized v(x, z) VTI blocks. The iterations are stopped when the residual moveout for at least two reflectors in each factorized block is close to zero. The overall organization of our MVA algorithm is similar to that developed by Liu (1997) for isotropic media, but the VTI model is characterized, for P-waves, by two additional parameters – and δ. The tomographic update of the medium parameters is based entirely on the residual moveout of events in image gathers. Although it is difficult to express the migrated depth zM in laterally heterogeneous media analytically, the residual moveout equation can be cast in a form similar to that obtained in Paper I for homogeneous models: 2 2 zM (h) ≈ zM (0) + A h2 + B
h2
2h4 , 2 (0) + zM
(4)
where h is the half-offset, and A and B are dimensionless constants that describe the hyperbolic and nonhyperbolic portions of the moveout curve, respectively. Numerical tests (see below) confirm that equation (4) with fitted coefficients A and B provides a good approximation for P-wave moveout in long-spread image gathers. To apply equation (4) in velocity analysis, we first pick an approximate value of the zero-offset reflector depth zM (0) on the migrated stacked section. The parameters A and B are obtained by a 2-D semblance scan on image gathers at each migrated zero-offset depth point. The best-fit combination of A and B that maximizes the semblance value is substituted into equation (4) to describe the residual moveout. It should be emphasized that the coefficients A and B in our algorithm are not directly inverted for the model parameters; their only role is in providing a close approximation for the residual moveout. After estimating the residual moveout in image gathers, we update the N -element parameter vector λ by solving the system of linear equations, AT A ∆λ = AT b .
(5)
2000 0
depth (m) (m) depth
Knowledge of kz and the zero-offset time t0 is sufficient for obtaining the anellipticity parameter η from equation (2) applied to one or both reflection √ events. The remaining√two key quantities, Vnmo = VP 0 1 + 2δ and ˆx = kx 1 + 2δ, can then be computed from equak tion (1), if the effective NMO velocities are determined at two or more locations x. We conclude that the moveout of horizontal events at two different depths and two image locations may provide enough information to estiˆx , and η. To decouple mate the parameters Vnmo , kz , k the horizontal gradient kx from the coefficient δ and determine the anisotropic coefficient , the velocity VP 0 has to be known at a certain point within the factorized block.
surface coordinate cmp coordinate (m)(m) 4000
A a
1000
B b
2000
(a)
(b)
Fig. 1: (a) Image of two reflectors embedded in a factorized VTI medium obtained using a homogeneous isotropic velocity field with VP 0 = 2600 m/s. (b) Semblance contour plot computed from equation (4) for the shallow reflector at the surface location 3 km.
Here A is a matrix with M · P rows (P is the total number of image gathers used in the velocity analysis and M is the number of offsets) and N columns that includes the derivatives of the migrated depth with respect to the medium parameters. The superscript T denotes the transpose, and the vector b contains the migrated depths that define the residual moveout.
Example with a single factorized layer First, we consider two irregular reflectors embedded in a factorized v(x, z) VTI medium (Figure 1a) with kz > kx > 0 and a positive value of η typical for shale formations [VP 0 (x = 3000 m, z = 0 m) = 2600 m/s, kz = 0.6 s−1 , kx = 0.2 s−1 , = 0.1, δ = −0.1]. A homogeneous, isotropic medium with the velocity VP 0 = 2600 m/s is chosen as the initial model for prestack depth migration. We start the velocity-updating process by manually picking along both imaged reflectors to outline their shapes. The velocity analysis operates with the image gathers at 12 equally spaced surface locations between 3 km and 4.2 km. The maximum offset-to-depth ratio for the selected image gathers at the shallow reflector is close to two, which is marginally suitable for estimating the parameter η. Tighter constraints on η are provided by the NMO velocities of reflections from the dipping segments of the shallow reflector (the dips exceed 30◦ in the middle of the section). Equation (4) is used to compute two-parameter (A and B) semblance scans for each reflector and evaluate the residual moveout in the image gathers (Figure 1b). Although a certain degree of trade-off exists between A and B, any pair of values inside the innermost semblance contour gives almost the same (quite small) variance of the migrated depths. Note that the interplay between A and B is similar to that between the NMO velocity and parameter η in the inversion of P-wave nonhyperbolic reflection moveout (Grechka and Tsvankin, 1998; Tsvankin, 2001). After the residual moveout has been estimated, we fix the vertical velocity VP 0 (x = 3000 m, z = 0 m) = 2600 m/s at the correct value and update the parameters kz , kx , , and δ using equation (5). The stacked images after
Migration velocity analysis and anisotropy 2000 0
surface coordinate cmp coordinate (m) (m) 4000
surface coordinate cmp coordinate (m)(m) surface coordinate (m)
2000 0
4000
1000
2000
1000
depth (m)
depth(m) (m) depth
depth(m) (m) depth
0
2000
(a)
surface coordinate (km) 5 10
2000
(b)
Fig. 2: Stacked image for the model from Figure 1 after (a) four iterations; (b) eight iterations.
four (Figure 2a) and eight (Figure 2b) iterations illustrate the improvements in the focusing and positioning of both reflectors during the velocity update. The inverted model parameters are close to the correct values: kz = 0.58 ± 0.02 s−1 , kx = 0.2 ± 0.0 s−1 , = 0.12 ± 0.01, and δ = −0.09 ± 0.01.
Test for a multilayered model Next, the algorithm is applied to a three-layer model shown in Figure 3. Each layer contains two reflecting interfaces, as required in our method, with every second reflector serving as the boundary between layers. The first and third layers are vertically heterogeneous [v(z)] and isotropic, while the second layer is a factorized, laterally heterogeneous [v(x, z)] VTI medium. All interfaces are quasi-horizontal, with the largest dips (at the flanks of the syncline) of 10◦ or less. The model is designed to represent a typical depositional environment in the Gulf of Mexico, where anisotropic shale layers are often embedded between isotropic sands. For the velocity analysis we use image gathers located along the left flank of the syncline with the surface coordinates ranging from 4400 m to 5600 m; the maximum offset-to-depth ratio for the image gathers is close to two. The parameter estimation is carried out in the layerstripping mode starting at the surface. For the first (top) layer, the vertical velocity is assumed to be known at a single surface location [VP 0 (x = 4000 m, z = 0 m) = 1500 m/s], which is sufficient for obtaining accurate values of kz , kx , , and δ (Figure 4). If we also use the correct vertical velocity (e.g., found from borehole data) at the top of the second and third layers, the inverted interval parameters for the whole model are close to the true values. The shapes and depths of the reflectors imaged for the reconstructed velocity field (Figure 5) are practically indistinguishable from those on the true section (Figure 3). If no borehole information is available, the velocity VP 0 can be assumed to be a continuous function of depth at a certain horizontal coordinate. To identify this point of continuity at the boundary between the first and second layers, we examine the moveout along the third and fourth reflectors (only for offsets smaller than 1000 m) after migration with an isotropic homogeneous velocity field in the second layer. The migration velocity was cho-
Fig. 3: True image of a three-layer factorized medium. Every second reflector (indicated here with arrows) represents the bottom of a layer. The parameters of the first subsurface layer are VP 0 (x = 4000 m, z = 0 m) = 1500 m/s, kz = 1.0 s−1 , and kx = = δ = 0; for the second layer, VP 0 (x = 4000 m, z = 800 m) = 2300 m/s, kz = 0.6 s−1 , kx = 0.1 s−1 , = 0.1, and δ = −0.1; for the third layer, VP 0 (x = 4000 m, z = 1162 m) = 2718 m/s, kz = 0.3 s−1 , and kx = = δ = 0. 1 0.5 0
kx
kz
ε
δ
Fig. 4: Estimated (◦) and true (?) parameters of the first layer obtained using the correct VP 0 (x = 4000 m, z = 0 m) = 1500 m/s on the surface.
sen to be equal to the estimated velocity at the bottom of the first layer (i.e., at the second reflector). To select the point of continuity, we pick the surface coordinate with the smallest residual moveout on the image gathers at the third and fourth reflectors. This criterion yielded x = 3900 m, which is sufficiently close to the true point of continuity for the second reflector (x = 4000 m, Figure 8). Using the estimated vertical velocity at x = 3900 m [VP 0 (x = 3900 m, z = 800 m) = 2316 m/s], we obtain the parameters of the second layer with high accuracy. (Figure 6). Applying the criterion of minimum residual moveout, we find the point of continuity between the second and third layers at (x = 5937 m, z = 1483 m). Although the estimated continuity point is shifted by 1000 m from its true location, the results of the velocity analysis for the third layer (Figure 7) are quite satisfactory, and the imaged reflectors are well focused and properly positioned (Figure 8).
Discussion and conclusions We developed an iterative migration velocity analysis (MVA) algorithm to estimate the parameters of VTI media composed of factorized blocks. At each iteration step, the residual moveout of P-wave reflection events in image gathers is minimized by solving a system of linear equations for the parameter updates. Since the parame-
Migration velocity analysis and anisotropy
surface coordinate (km) 5 10
0 depth (m)
ter estimation is performed in the post-migrated domain, the algorithm is robust in the presence of random noise, significant lateral heterogeneity, and dipping structures.
2000
Fig. 5: Stacked image obtained after prestack depth migration with estimated medium parameters. The vertical velocity VP 0 at the top of each layer was known.
0.6 0.4 0.2 0 −0.2
k
x
k
z
ε
δ
Fig. 6: Estimated (◦) and true (?) parameters of the second layer obtained assuming that VP 0 is continuous between the first and second layers at the point (x = 3900 m, z = 800 m).
The residual moveout of both horizontal and dipping events in each factorized block is governed by the NMO velocity Vnmo , the vertical gradient kz , the effective latˆx , and the anisotropic coefficient η. Stable eral gradient k recovery of all four parameters requires reflection moveout from at least two interfaces within each block sufficiently separated in depth. However, the vertical velocity VP 0 , needed to reconstruct the velocity field in the depth domain, is generally unconstrained by P-wave reflection moveout. The velocity VP 0 can often be estimated from borehole data using either check shots or sonic logs. If no borehole information is available, a suitable model for depth imaging can be constructed by assuming that VP 0 is continuous across layer boundaries. Then, given the value of the vertical velocity at a single point on the surface, the entire velocity model in depth can be estimated from the residual moveout of P-wave reflection events. It should be mentioned that if the overburden contains interfaces with significant dip or curvature, P-wave reflection moveout and, therefore, residual moveout on image gathers become dependent on the vertical velocity and the parameters and δ (Grechka et al., 2002). For models of this type, the layer-stripping approach adopted in our MVA algorithm is not always adequate because the parameters of a given layer may remain unconstrained in the absence of reflection data from deeper interfaces (Le Stunff et al., 2001).
0.6
References
0.4 0.2 0 −0.2
k
x
k
z
ε
δ
Fig. 7: Estimated (◦) and true (?) parameters of the third layer obtained assuming that VP 0 is continuous between the first and second layers at the point (x = 3900 m, z = 800 m) and between the second and third layers at the point (x = 5937 m, z = 1483 m).
depth (m)
0
surface coordinate (km) 5 10
2000
Fig. 8: Stacked image obtained by assuming the vertical velocity to be continuous at the points marked by (◦). The true points of continuity are marked by (•).
Alkhalifah, T. and Tsvankin, I., 1995, Velocity analysis for transversely isotropic media: Geophysics, 60,1550–1566. Grechka, V., and Tsvankin, I., 1998, Feasibility of nonhyperbolic moveout inversion in transversely isotropic media: Geophysics, 63, 957–969. Grechka, V., Pech, A., and Tsvankin, I., 2002, P-wave stacking-velocity tomography for VTI media: Geophys. Prosp., 50, 151–168. Liu, Z., 1997, An analytical approach to migration velocity analysis: Geophysics, 62, 1238–1249. Le Stunff, Y., and Jeannot, J.P., 1998, Pre-stack anisotropic depth imaging: 60th EAGE Conference, Extended Abstracts. Le Stunff, Y., Grechka, V., and Tsvankin, I., 2001, Depth-domain velocity analysis in VTI media using surface P-wave data: Is it feasible?: Geophysics, 66, 897–903. Sarkar, D., and Tsvankin, I., 2002, Analysis of image gathers in factorized VTI media: 72nd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 2305–2308. Tsvankin, I., 2001, Seismic signatures and analysis of reflection data in anisotropic media: Elsevier Science Publ. Co., Inc.