Method for Tomographically Laterally Varying Seismic Data Velocity Estimation

1. Field of the Invention.

The present invention relates to the processing of seismic data utilizing a tomographic process and, more particularly, but not by way of limitation, it relates to an improved method for velocity estimation of seismic data using the seismic tomography processing.

2. Brief Description of the Prior Art.

The prior art processes for ray-tracing evaluated the subsurface in terms of determinative cells of constant slowness values which were then reconstructed to approximate seismic ray paths, still other types employed finite difference ray-tracing methods. The prior types using constant slowness values for analyzing individual cell structures culminated in very severe reconstruction resolution limitations. The alternative finite difference tracing methods introduced uncertainties that make them unsuitable for iterative estimation schemes. None of the prior art types treat the field within a cell as continuous, and an assumption must be made that the cell either is constant in value or constant in the value of its gradient. This approach brings about severe limitation in terms of resolution and continuity of the reconstructed field. The prior art schemes also employ computationally expensive finite-difference ray-tracing schemes that are not always suitable for parameter adjustment algorithms and other attempts to generalize such diverse data.

In the present method, the seismic section or presentation medium is divided into cells of a preselected size and either slowness or velocity values are specified at the cell corners. Source and receiver positions are known relative to the medium and rays are plotted from the source at specific angular increments from a minimum to a maximum angle with subsequent iterative adjustment of tracing angles, i.e., source-reflector-receiver. The basic elements of the tomographic process for velocity estimation consists of the estimation of travel times of the rays through the medium by forward modeling, plus a following model parameter estimation. Ray tracing derives from individual cell analysis as the parameters of the medium are the velocities at the individual cell corners and the depths of reflector along vertical cell walls. The velocity field in a cell is reconstructed using the velocity values at the cell corners in a local two-dimensional neighborhood and ray-tracing is carried out utilizing derived expressions for position, orientation, and travel time of the ray in the cell as a function of cell dimensions, initial ray position, ray orientation, and the velocity values within a local two-dimensional neighborhood. Linearization is then used to set up a travel-time equation for parameter estimation, and derivations of travel time relative to velocities and reflector position enable construction of the model parameter estimation leading to finalized estimated data readout .

FIG. 1 is a flow diagram of a general tomographic velocity estimation computational flow;

FIG. 2 is a graphic depiction illustrating emanation of seismic rays from a source for different initial shooting angles;

FIG. 3 is a graphic depiction of bilinear interpolation of velocity values in a designated cell interior;

FIG. 4 is a graphic depiction of a ray patch within a cell;

FIG. 5 is a graphic depiction of a local neighborhood of a designated cell;

FIG. 6 is a graphic illustration of incident and reflected rays/from a reflector in the cell;

FIG. 7 is a graphic illustration of perturbation of the travel path due to change in the right-side reflector end point;

FIG. 8 is a graphic illustration of perturbation of the travel path due to change in the left-side end point of the reflector;

FIG. 9 is a graphic showing of travel times and points of intersection of the ray with the reflector for each of three reflector positions;

FIG. 10 illustrates graphically the point of intersection of a ray with the middle reflector position;

FIGS. 11A, 11B and 11C illustrate parameters of the medium in relation to an eigenray; and

FIG. 12 is a detailed flow diagram of tomographic parameter estimation process.

The general computational flow for tomographic velocity estimation is shown in FIG. 1. In relation thereto, reference will be made to the travel time measurements taken for a constant offset section; however, the concepts are equally applicable with configurations of data such as source gathers, common basement point data, etc.

The travel times between a source and receiver for specific reflectors are measured in constant offset sections of prestack data and these constitute the true travel times. Starting with an initial depth model, travel times are estimated by tracing rays, and subsequently the parameters of the model are updated by successively comparing the estimated travel times with the true travel times as will be described. The ray tracing and the parameter adjustment by iteration are continued until certain convergence criterion is satisfied. Thus, flow stage 10 of FIG. 1 inputs estimated travel times and designates sensitivity matrix elements by shooting the source-reflector-receiver rays with the initial model. The true field data travel times are available at flow stage 12 for input with data from flow stage 10 to flow stage 14 which functions to compute the ray-based parameter corrections.

In first parameterizing the depth model, the medium or data section is divided into a predetermined array of cells of certain designated size. The velocity file is parameterized in terms of velocities at the cell corners and reflectors are described in terms of depth of reflectors along vertical cell walls, as will be described. The model parameters then are the velocities of the medium at the cell corners and the depth of reflectors along cell walls.

FIG. 2 illustrates in generalized form a global coordinate system wherein the data medium is divided into cells of certain size where B(1), B(2), T(1) and T(2) define the medium expanse relative to the coordinate axes x and y. Thus, to trace a ray between a source position S and a receiver position R for reflector AA, rays are shot from S at different initial angles .theta. through .theta..sub.n for interpolated receiver positions R through R.sub.n by iterative establishment. The key computation of ray shooting from a source position at a certain initial angle is the tracing of a ray through a cell which depends upon the presence or absence of a reflector within the particular cell. Knowing the ray position and orientation as the point of entrance to a cell, and presence of a reflector, then enables computation of the ray position and orientation at the point of exit from the cell, as will be further described.

Interval velocities for the initial model are selected based upon available a priori information, e.g., well data or the like, or they may be derived from stacked velocities or by estimating interval velocities using stack optimization. For a particular basement point the equation that relates the interval velocities VI.sub.j, the RMS velocities VR.sub.i and zero offset travel times t.sub.i may be given as ##EQU1## Knowing the RMS velocities VR for the above equation can be inverted to yield the interval velocity VI.

Flow stage 16 (FIG. 1) then commences the operation of parameter adjustment as the stage functions to estimate parameter updating terms as a weighted combination of ray-based corrections as applied to data in flow stage 18. The ray tracing between a specific source S and receiver R for update of the initial model parameters by estimating the travel times and other necessary information is done in a manner as follows. Referring again to FIG. 2, rays are shot from the source S at specific angular increments, starting from a certain minimum angle and proceeding until a certain maximum angle is reached. The angular positions for which the rays and the receivers are obtained by iterative establishment of shooting angles. When the medium is divided into cells 20, the heart of the calculation involves tracing a ray through the cells in a designated manner.

After update of parameters at stage 18, a decision stage 22 performs the designated criteria tests which, if satisfied, passes flow through the final decision stage 24 to output. If criteria is not satisfied in stage 22, then flow recycles to repeat from stage 14 computing ray-based parameter corrections. If complete ray tracing must be repeated as detected in flow stage 24, then stage 26 performs flow stage 10 functions with an updated model for recycling for repetition through stage 14.

FIG. 3 illustrates the manner in which bilinear interpolation of velocity values is determined within a cell's interior. A coordinate system 28 is assumed to be positioned so that the origin coincides with the bottom lefthand corner of a cell, i.e., the coordinate (0,0). Given the velocity values at the corners of a cell, the bilinearly interpolated velocity value at a point r interior to the cell can be obtained as: ##EQU2## where r equals (x, y) .sup.T is the position vector of a point (x, y) in the cell, V.sub.o is the velocity value at the origin, V(r) is the interpolated velocity value of the point (x, y), and T denotes matrix transposition. Let V(x, y) be the velocity value at the point (x, y), and let a and b denote the dimensions of the cell. The quantities H.sub.1 and H.sub.2 are then defined as ##EQU3## One of the important features of the bilinear interpolation is that the interpolated velocity field is continuous across the cell boundaries from cell to cell.

It is next necessary to have a scheme for tracing the rays within the cell when the velocities are specified at the cell corners. Thus, given the position and the orientation of the ray entering the cell, the expressions may be derived for obtaining the position, orientation, and travel time of the ray as a function of path parameter. Referring to FIG. 4, the quantity r.sub.o is the position and d.sub.o is the orientation of the ray entering the cell. The equation which describes the path of a ray in an inhomogeneous medium can be written as ##EQU4## where ds is the differential path length along the ray, r is the position vector on the ray, and V(r) is the velocity value at position r. Integration of Equation (7) yields an expression for the position r as a function of path length s and the resulting expression is ##EQU5## Expanding r(s) in Taylor's series around s=o ##EQU6## Then from Equation (7) we obtain ##EQU7## and putting Equation (10) into Equation (9) yields ##EQU8## From Equation (2), after keeping only the terms of first order in H.sub.1 and H.sub.2, we get ##EQU9## Substitution then of Equation (12) and Equation (13) in Equation (9) yields ##EQU10##

The travel time in the cell as a function of path parameters is given as ##EQU11## Then, from Equations (14) and (16) we get ##EQU12## Equation (17) can be written as ##EQU13## where Equation groups 19, 20 and 21 are ##EQU14## and it then follows that ##EQU15##

FIG. 5 illustrates a local neighborhood of a cell in idealized reference form. A cell 30 is under consideration. The tracing of rays within the cell 30 is enabled by local neighborhood reconstruction of the velocity field. A two-dimensional local neighborhood 32 (.GAMMA..sub.o) of size M.times.N of the cell 30 is shown as being positioned at a designated place in the coordinate system of FIG. 5. The velocity field at a point in the cell 30 is reconstructed using the velocity values in its local neighborhood as ##EQU16## It is assumed that M and N are odd numbers. In terms of S's, the orthonormal basis wave forms, defined over the local neighborhood .GAMMA..sub.0, (24) becomes ##EQU17## where X is coordinate vector and the functions S's are defined in terms of polynomials orthogonal over one-dimensional neighborhood as ##EQU18## and X.sub.i is the domain of coordinate variable x.sub.i. The estimates of the coefficients a's are given as Use of Equations (26) and (27) simplify the denominators of equation (28) to ##EQU19##

In terms of matrix notation, equation (24) becomes ##EQU20## Let the ray be entering the cell at position X.sub.o and let V(X.sub.o) be the velocity field at X.sub.o, then from Equation (30), V(X) can be expressed as ##EQU21## Solving the Eikonel equation, expressions for the position X(s), orientation d(s), and the travel time t(s) of the ray in the cell can be written as ##EQU22## Using equations (31) and (35), the travel time expression can be rewritten as ##EQU23## Substitution of Equation (28) in Equation (36) yields ##EQU24## where X.sub.1 and X.sub.2 are domains of x.sub.1 and x.sub.2 and c(x.sub.1, x.sub.2) and e(x.sub.1, x.sub.2) are the weights of velocities at the cell corners of local two-dimensional neighborhood.

Next, the travel time must be expressed as a weighted combination of the velocities. For this, the Equation (38) can be written as ##EQU25## where C.sub.l and e.sub.l are the weights and v.sub.l is the velocity at the l cell corner, 1.ltoreq.l.ltoreq.MN. Let .beta..sub.l be defined as ##EQU26## Then from Equations (39) and (40) we derive ##EQU27## From Equation (41), after making use of equation (39), there is obtained ##EQU28## Differentiating Equation (41) with respect to v.sub.j yields ##EQU29## And using Equation (41), we get from Equation (43) ##EQU30## Equations (42) and (44) follow after noting that t and .beta.'s are homogeneous functions of velocities.

Ray tracing can also be carried out when slowness values are specified at the cell corners. At a point (x, y) in the cell, the bilinearly interpolated slowness value can be written as ##EQU31## The quantities G.sub.1 and G.sub.2 are defined as ##EQU32## where relationships are as follows: ##EQU33## Proceeding as in the preceding section relating to bilinear interpolation for the estimation of a velocity field (see FIG. 3), the expressions for the position, orientation, and travel time of the ray in the cell as a function of path parameter s can be obtained as ##EQU34## Equations (49), (50) and (51) can then be written as follows: ##EQU35## Where values derive as follows: ##EQU36## Equation (54) expresses the travel time in the cell as a linearly weighted combination of slowness values at the four corners of the cell.

The case is altered when there is a reflector in the cell. As shown in FIG. 6, there is depicted an incident ray 34, and a reflected ray 36, reflecting about a line 38 normal to the reflector point 40 in the cell. The position and orientation of the ray in the cell as a function of path parameter s can be written as ##EQU37## The reflector 40 in the cell is assumed to be described as a piece of cubic spline. That is ##EQU38## The coordinates (x.sub.2r, x.sub.1r), the point of intersection of the ray 34 with the reflector in the cell is computed from ##EQU39## and, where s.sub.r is a root of the equation ##EQU40## When the slope of the reflector is very small and the reflector is very close to a horizontal boundary, the coordinates of the point of intersection of the ray with the reflector can be approximately computed from ##EQU41## Let .nu. be the slope of the reflector at x.sub.1r and x.sub.2r. That is ##EQU42## The orientation of the reflected ray at the point (x.sub.1r, x.sub.2r) is given by (cos.theta.r, sin.theta.r), where ##EQU43##

Another function of flow stage 10 (FIG. 1) requires estimation of sensitivity matrix elements. The derivatives of travel time with respect to the velocities at the cell corners can be obtained from Equation (42). The estimation of derivatives of travel time with respect to the reflector position parameters is described below. In a first case for estimation without perturbing the reflector, reference is made to FIG. 7.

The change in the travel time of the reflector with respect to the change in position of the end point of the reflector can be obtained for both right-and left-side end points. For the right-side end point as shown in FIG. 7, let TT' be the tangent to the reflecting surface, and A(x.sub.r, y.sub.r) be the point of intersection of incident ray with the reflecting surface. Let .theta. define the angle of incidence and .beta. define the orientation of TT'. Let the right-side end point be changed from Q to Q', and let .DELTA.Y.sub.m define the change. Finally, let the ray path in the cell change from EAF to EBF. Assuming then that the changes in the reflector within the cell are small, we can write ##EQU44## As indicated in FIG. 7, we obtain ##EQU45## Let .DELTA.t be the travel time perturbation. That is ##EQU46## where V(x.sub.r, y.sub.r) is the value of the velocity field at the point (x.sub.r, y.sub.r). From Equations 70-72, the change of travel time of the ray in the cell with respect to the change in the right-side end point of the reflector can be obtained as ##EQU47##

Referring to FIG. 8, there is shown the change in the travel path of the ray for a change in left-side end point of the reflector. Proceeding as before, the change of travel time with respect to the change in the left-side end point of the reflector can be obtained as ##EQU48##

The derivatives of travel time with respect to the reflector position parameters can also be estimated by perturbing the reflector up and down by known amounts. Referring to FIG. 9, there is shown the rays between source S and receiver R during controlled perturbation up and down. The travel time of the ray and the coordinates of the point of intersection of the ray with the reflector are recorded for each reflector position. Let tu, tc, and td be the travel time of the rays with respect to the upper, middle, and lower reflector positions. Let {(x.sub.1u, x.sub.2u), (x.sub.1c, x.sub.2 c), (x.sub.1d, x.sub.2d)} be point of intersection of the ray with the upper, middle, and lower reflector positions, respectively.

Then let ##EQU49## The derivative of the travel time with respect to the x.sub.2 coordinate of the point of intersection of the ray with the middle reflection position can be estimated from ##EQU50## Letting the ray intersect the middle reflector in the (m-1) th horizontal cell position as depicted in FIG. 10, and letting .beta..sub.r be the angle of tangent of the reflector at the point of intersection of the ray with the middle reflector position, the derivatives of travel time with respect to the parameters of the reflector at the right-and left-side cell walls of the point (x.sub.1c, x.sub.2c) can be estimated from ##EQU51##

It is also necessary to derive the basic travel time equation. As the model parameters are updated (flow stage 18 of FIG. 1) from iteration k to iteration k+1, FIGS. 11A, 11B, and 11C pictorially depict in generalized form the medium parameters and the corresponding eigenray at various stages of parameter updating. Letting V.sup.k, p.sup.k, and j.sup.k be the medium velocity vector, reflection position parameter vector and the eigenray at the kth iteration, then j.sup.k is a function of V.sup.k, p.sup.k. For a particular source, receiver and reflector position, let t.sup.k be the estimated travel time of eigenray at the kth iteration. We can then write ##EQU52## where v.sup.k (s) the velocity along the ray j.sup.k and s is the path parameter, letting ##EQU53## The correction terms .DELTA.V.sup.k and .DELTA.p.sup.k are estimated such that the travel time t.sup.k+1 is made equivalent to true travel time t. Let ##EQU54## Define the variable I to denote the integrals ##EQU55## Expanding t.sup.k+1 around t.sub.1.sup.k+1, we get ##EQU56## Again using the linearization process, we can write ##EQU57## where the products of terms of first order smallness are assumed to be neglected. In both Equations (83) and (84), the change of I with respect to the change in the eigenray j is neglected to a first order approximation. Use of Equation (84) in (83) yields ##EQU58## where .xi..sub.l.sup.k are the derivatives of travel time with respect to the reflection position parameters. Letting {.beta. .sub.l.sup.k+1 } pertain to path j.sub.l.sup.k+1 and {.beta..sub.l.sup.k } pertain to path j.sup.k, then from Equation (42) we can write ##EQU59## Consider ##EQU60## Expanding .beta..sub.i.sup.k+1 around .beta..sub.i.sup.k there follows ##EQU61## From Equations (44) and (90), we obtain ##EQU62## Using Equation (91), (89) and (86) in Equation (85) yields ##EQU63## Letting t.sup.k+1 equal t, the true travel time, while ignoring the product of terms of first order smallness, the Equation (92) can be written as ##EQU64## For a particular source, receiver, and reflector position, Equation (93) relates the travel time errors to the model parameter correction terms under the linearization assumptions discussed above.

In the estimation of model parameters, the program process requires routines (a) for parameter updating when there are no sharp contrasts present in the velocity as well as the alternative mode of parameter updating (b) when there are sharp contrasts present in the velocity field. A ray based algorithm for the computation of parameter correction terms is presented below. The ray tracing and parameter correction parts are arranged as given in the flow chart of FIG. 1. From Equation (93), the travel time equation for the jth ray at the kth iteration can be written as ##EQU65## where .DELTA.v.sub.jl.sup.k and .DELTA.p.sub.jl.sup.k are based correction terms and are estimated as ##EQU66## The corrections to the velocity and the reflector position parameter vectors are made using the relationships ##EQU67## where n.sub.l and m.sub.l are the number of rays contributing to the th velocity parameter v.sub.l and l th reflector position parameter p.sub.l. The quantities .alpha..sub.v and .alpha..sub.p are estimated from ##EQU68## The parameters are adjusted iteratively until there is no appreciable reduction in the travel time error or until the travel time error begins to increase.

When sharp contrasts are present in the velocity field, the parameter estimation is arranged around swaths of rays. For a particular swath containing n rays, the estimation equation will be of the form

Where Z is a nxl travel time error vector, H is a known matrix, x is an unknown parameter vector (velocities at the cell corners and reflector position parameters) and v is a noise vector. Assuming the variables are gaussian, then a Bayesian solution for the iterative updating of unknown parameter vector x can be written as ##EQU69## The above scheme requires ray tracing and parameter updating for each swath of rays and the computations are stopped whenever certain convergence condition is satisfied.

FIG. 12 represents a detailed flow diagram of the tomographic parameter estimation process. Initial parameters will all have been written to file. That is, the .beta. parameters consisting of derivatives of travel time with respect to the velocities at the cell corners, and the .xi. parameters consisting of the derivatives of travel time with respect to the reflector position parameters, the number of rays visiting each cell corner IR(TI) array, the number of rays visiting each reflector position parameter IR(SI) array, and the travel time estimate of each ray.

In accordance with flow diagram 50, FIG. 12, flow stage 52 functions to read from the ray tracing files each of the sensitivity matrix elements: that is, the derivatives of travel time with respect to velocities at the cell corners--the .beta. parameters (Equation 41); the derivatives of travel time with respect to the reflector position parameters--.xi. parameters (Equations 77 and 78 or Equations 73 and 74); the number of rays visiting each cell corner--IR(T1) array; the number of rays visiting each reflector position parameter--IR(SI) array, and ray travel time estimates t.sub.p as derived by Equations (18), (36) or (51).

The flow stage 54 then initializes the iteration count and stage 56 initializes the sum of corrections for the velocities at cell corners, SOv (k), and the sum of the corrections for the reflector position parameters SOp(k). Flow stage 58 then computes the sum of the corrections for all rays for each velocity and reflector position parameter utilizing the relationships of Equations (96) and (97). Flow stage 60 then normalizes the correction terms and function proceeds to flow stage 62 for computation of alpha quantities using the Equations (102).

The next flow stage 64 functions to compute the estimated travel time correction terms while flow stage 66 applies or shoots the travel time corrections to the data with subsequent update of the travel time estimates and the model parameters as carried out at flow stage 68. Finally, convergence stage 70 functions to make a decision in response to predetermined convergence criteria and stops all computations when convergence is reached; that is, where there is no appreciable decrease in the travel time errors or the process has carried out a specified number of iterations. If decision stage 70 is affirmative then data is output, but if not then the process recycles through flow stage 72 to repeat iterations or begin again at flow stage 56 with initialization of the sum of correction terms through designated iterations. This cycling process continues until a convergence is satisfied to provide an output from decision stage 70.

The foregoing discloses a novel seismic processing technique that enables accurate determination of laterally varying velocity estimations over an extended assemblage of seismic section data. It is a subgoal of the objective to develop velocity estimation processes that will aid in depth migration and contribute to greater accuracy in seismic imaging. The present invention utilizes a tomographic technique in processing a line of seismic section data such that a move-along interval velocity estimation process is brought about.

Changes may be made in combination and arrangement of steps as heretofore set forth in the specification and shown in the drawings; it being understood that changes may be made in the embodiments disclosed without departing from the spirit and scope of the invention as defined in the following claims.

The embodiments of the invention in which an exclusive property or privilege is claimed are as follows: