Estimating interval velocities

Abstract

A method of estimating a velocity of a geological layer includes a. providing a first, initial model including an interval velocity associated with a subsurface location and an uncertainty associated with the interval velocity; b. providing data including an actual or approximated root-mean-square (RMS) velocity associated with a subsurface location and an uncertainty associated with the RMS velocity; and c. estimating a second model including an interval velocity associated with a subsurface location and an uncertainty associated with the interval velocity, based on the interval velocity and the uncertainty of the first model, and the RMS velocity and the uncertainty of the data.

Claims

1. A method of imaging a velocity structure of a subsurface region of the Earth, the method comprising: providing data comprising a data interval velocity of the subsurface and a data uncertainty associated therewith; determining an output grid onto which velocity profiles are to be estimated, and calculating a separation between each pair of desired grid locations; using said separation between locations to evaluate a correlation function, and using correlation values to form a covariance matrix that provides a model of a covariance between a velocity to be predicted in a prediction location on said output grid and said data interval velocity; providing a first, initial model based on said data comprising a first model interval velocity of the subsurface and a first model uncertainty associated therewith; performing a spatial prediction of the data interval velocity to form a second model comprising a predicted second model interval velocity of the subsurface and a second model uncertainty associated therewith, the predicted second model interval velocity being associated with the prediction location in the subsurface region, wherein the prediction is performed based on the first model interval velocity and the first model uncertainty, and based on the data interval velocity and the data uncertainty, and based on said covariance matrix; displaying a 2D or 3D spatial image of the velocity structure of the subsurface region from the second model.

2. The method as claimed in claim 1, wherein the first model interval velocity is associated with a first model location in the subsurface region and the data interval velocity is associated with a data location in the subsurface region, and one or more of the first model location, the data location and/or the prediction location is a different subsurface location.

3. The method as claimed in claim 1, wherein performing the prediction includes predicting an interval velocity of the data from a first, data location, to a second, prediction location, to thereby determine an interval velocity for the second model for the second location.

4. The method as claimed in claim 1, wherein performing the prediction includes predicting a plurality of data interval velocities associated with a set of subsurface data locations from the subsurface data locations to a set of subsurface prediction locations, to thereby determine a plurality of second model interval velocities for the set of prediction locations.

5. The method as claimed in claim 1, wherein performing the prediction includes performing a Bayesian prediction.

6. The method as claimed in claim 1, wherein performing the prediction includes performing a constrained least-squares prediction.

7. The method as claimed in claim 1, wherein the interval velocities of the data are obtained from performing a Dix inversion of actual or approximated RMS velocities.

8. The method as claimed in claim 1, wherein the first model, the data and/or the second model are represented by Gaussian distributions or multi-Gaussian distributions, where the data and the first and second models include a plurality of interval velocities.

9. The method as claimed in claim 1, wherein the first model interval velocity is associated with a first model location in the subsurface region and the data interval velocity is associated with a data location in the subsurface region, and one or more of the first model location, the data location and/or the prediction location is the same subsurface location.

10. The method as claimed in claim 1, further comprising the steps of: building a depth or time-migrated image of the Earth's subsurface by using seismic data provided based on the predicted second model interval velocity; and detecting or exploring for hydrocarbons based on the depth or time-migrated image of the Earth's subsurface.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) There will now be described, by way of example only, embodiments of the invention with reference to the accompanying drawings, in which:

(2) FIG. 1 is schematic representation of a region of the earth containing geological layers;

(3) FIG. 2 is a flow chart relating to a method for estimating interval velocities;

(4) FIG. 3 is a graph showing results from an inversion using a first parameterisation;

(5) FIG. 4 is a graph showing results from an inversion using a second parameterisation;

(6) FIG. 5 is a graph showing results from an inversion using a third parameterisation; and

(7) FIG. 6 is a graph showing results following a lateral prediction procedure applied to interval velocity data at a depth of 1500 m.

DETAILED DESCRIPTION OF THE INVENTION

(8) In this example, an estimate of the Dix interval velocity is obtained from RMS velocity data using the relations of Equations 1 to 3 as outlined above. The estimate in this example is determined by a Bayesian inversion (although in other embodiments a stabilised least squares inversion technique could be used). The Bayesian inversion technique uses Bayes' Rule to provide a solution in the form of a conditional, posterior probability function p(m|d) (the probability of m given d), expressed as:
p(m|d)p(d|m)p(m)(Equation 4)

(9) where

(10) m is a vector of model parameters

(11) d is a vector of data

(12) p(d|m) is the likelihood function (the probability of d given m)

(13) p(m) is the prior distribution probability function

(14) In the present case, the data d comprises squared RMS velocity values for different travel times. RMS velocity and travel time pairs may be derived during a processing flow of time-series seismic reflection amplitude data for detecting subsurface structure, for example, at the stage of velocity analysis as known in the art of seismic data processing. In practice, the RMS velocities may be approximations to the RMS velocities. The RMS velocities may be approximated by for example stacking velocities or any other velocity that may be taken to represent an RMS velocity.

(15) The prior distribution p(m) in this case provides an initial or prior model (before taking account of the data and its uncertainty) comprising interval velocities and an uncertainty associated with each interval velocity. This initial model may be formed based on knowledge about the subsurface, such as its rock properties, which could for instance be obtained from geological well log samples and tests. The uncertainties may be determined accordingly.

(16) The data d and the prior model p(m) are provided in advance in order determine the posterior distribution p(m|d) and represent inputs in the present implementation, although they are not the only inputs, as described further below.

(17) The likelihood function or model p(d|m) in this case incorporates a forward model operator linking the RMS velocity data with the interval velocities of the model solution and an error vector containing uncertainties associated with each RMS velocity datum. This provides a generally non-deterministic link between the model and the data. The uncertainties in the data may be estimated during a seismic data processing step. For example, the processing step of semblance analysis could be modified to allow variances and temporal covariances to be determined. The data uncertainties are a further input to the inversion. The operator linking the squared RMS velocity data d and the squared Dix interval velocity parameters m, are based on the relationships of Equations 1 to 3. Mathematically, the forward model can be expressed as:
d=Gm+e(Equation 5)

(18) where

(19) d is the data vector

(20) m is the model parameter vector

(21) G is the forward operator matrix

(22) e is an error vector

(23) In this expression, G is the forward operator in the form of an integral/summation matrix, representing the link between the model squared Dix interval velocity parameters m and the data of squared RMS velocities d. The error vector e represents the associated uncertainties of the data. The error vector e constitutes a data error model which may be represented by means of a zero expectation and a set of covariances, and may be multi-Gaussian.

(24) The posterior distribution therefore takes into account pre-existing information, including the prior model of interval velocities and its uncertainties, and the RMS velocity data and its uncertainties, and provides a best estimate model including Dix interval velocities and associated uncertainties. Details of the numerical formulation for performing the inversion and determining the posterior distribution for this example are described in Appendix A.

(25) Implementation

(26) As shown in Appendix A, the posterior distribution is characterised by the posterior expectation vector .sub.m|d and the associated covariance matrix .sub.m|d according to the following equations:
.sub.m|d=.sub.m+.sub.mG.sup.T(G.sub.mG.sup.T+.sub.e).sup.1(dG.sub.m)(Equation 6)
.sub.md=.sub.m.sub.mG.sup.T(G.sub.mG.sup.T+.sub.e).sup.1G.sub.m(Equation 7)

(27) where

(28) d is the data vector

(29) .sub.m is the prior model expectation vector

(30) .sub.m is the prior model covariance matrix

(31) G is the forward operator linking the model to the data

(32) .sub.e is the covariance matrix of the error vector of the data

(33) G.sub.mG.sup.T+.sub.e.sub.d is the covariance of the prior model for the data vector

(34) The posterior distribution is determined by evaluating these equations using a computer program. These equations produce squared Dix interval velocities and corresponding covariances. The posterior expectation vector .sub.m|d defines the best estimate of (squared) Dix interval velocities for different depths. The associated covariance matrix .sub.m|d provides the corresponding variance for the squared Dix interval velocity in each relevant interval as well as the covariance between each pair of squared Dix interval velocity distributions. These are evaluated from input data and/or user-specified parameters, as detailed below.

(35) With reference firstly to FIG. 1, RMS velocity data are initially obtained in various lateral x,y positions 3 in a geological region of interest 1 where seismic time series reflection data have been acquired. The RMS velocity data typically comprises RMS velocities at several TWT points defining various intervals or subsurface layers 5. The TWT represents the travel time of a seismic pulse. This will include signal paths from the Earth's surface, into the subsurface where it is reflected and back again. In typical data processing procedures, seismic reflection amplitudes are plotted against two way time at different x,y locations. Primary reflections from deeper subsurface layers will tend to arrive at late times in the time series so that different TWTs can be representative of different depths and enable a depth image for a region of the Earth's subsurface to be produced. Generally, it is not determinable from such plots the true geometric structure of the subsurface, for example the true thickness between layers, without further estimating interval velocities of the medium of successive layers of the subsurface, in order to convert raw travel times to migrated time, or to convert the travel time to depth. The velocities may be used to reposition the events in time (i.e. a corrected travel time by a process of time-migration), and may be iterated with several steps of velocity analysis where RMS velocities are estimated.

(36) The RMS velocity can be approximated in the seismic data processing procedure, for example by seismic data processors selecting a stacking velocity. Each RMS velocity sample has a standard deviation associated with it, which may be derived when an RMS velocity is estimated from the seismic data. In addition, covariances within the RMS velocities can also be included. The standard deviation and the covariance could be determined upon selection of the stacking velocity.

(37) A prior model is also derived for the geological region 1. This model represents the best knowledge of a user of what the interval velocities should be, without knowledge of the sampled RMS velocity trend. Typically, this may take account of geological knowledge of the region and may include data regarding seismic or elastic rock properties from nearby wells or the like.

(38) In FIG. 2, an example workflow 10 for inverting, in 1D (or globally in 2 or 3 dimensions), the RMS velocity data using Bayes' Rule or equivalently a constrained least squares inversion is shown. The workflow 10 includes an initial stage of reading in various inputs 12, which are then operated on by an inversion routine 20 to evaluate Equations 6 and 7, and thereby generate an output 22a containing Dix interval velocities and associated uncertainties. In the case of a 1D inversion, only the top two options of the output 22a may be produced. As can be seen, the inputs 12 include an RMS velocity data set 14 containing RMS velocities obtained as described above. In practice, the RMS velocity data set 14 comprises collections of TWT-RMS velocity pairs for one or more horizontal (x,y) locations 3, and a standard deviation assigned to each RMS velocity. It will be appreciated that the velocity data locations do not need to be provided on a regular spatial grid.

(39) The inputs 12 include a further set of desired output TWTs 16 which are defined at the same x,y locations as the TWT-RMS velocity pairs of data set 14. The TWTs 16 are used to define time intervals for which Dix interval velocities are to be estimated in the output 22a. The TWTs 16 typically correspond to a top and bottom of the interval for which a Dix interval velocity estimate is desired, and may differ from those of the TWT-RMS velocity pairs of data set 14. Temporal covariances, or in the case of a global inversion additionally spatial covariances, can also be assigned in a suitable covariance matrix for the RMS data.

(40) In addition, the inputs 12 include a prior model 18 comprising expected or initial guess interval velocities as a function of x,y and of TWT and a standard deviation associated with each interval velocity. The covariance may also be specified. The expected interval velocities of the prior model 18 could be discrete or continuous with TWT provided they can be evaluated for the set of TWTs 16 and x,y locations in which interval velocities estimates are required in the outputs 22a, 22b.

(41) From the inputs, the matrices and/or vectors .sub.m, .sub.m, G, .sub.e, and d are formed and the posterior distribution p(m|d) determined by the inversion routine 20 in line with Equations 6 and 7. For a 1D inversion as described here, the inversion routine 20 performs the following operations repeated for each lateral (x,y) location. Firstly, the covariance matrix .sub.e of the error vector is formed from the standard deviations of the RMS velocities of the data set 14. This is done by squaring the standard deviations to form variances and placing the variances of the RMS velocities on the leading diagonal of the covariance matrix. Temporal covariances between errors on the RMS velocities, if known or estimated on input, are provided in the off-diagonal elements of the covariance matrix.

(42) The prior model expectation vector (.sub.m) is then formed comprising the squared expected interval velocities of the prior model 18 (the best guess model) for the desired TWTs 16 for the output. Thereafter, the prior model covariance matrix (.sub.m) is formed from the standard deviations associated with the interval velocities of the prior model 18. In a similar manner to the covariance matrix .sub.e, covariance matrix .sub.m for the prior model is formed by squaring the standard deviations to form variances and placing the variances on the leading diagonal of the matrix for each of the desired output TWTs 16. Covariances between parameters in the prior model are included by means of the off-diagonal elements of the covariance matrix .sub.m.

(43) The operator G of the likelihood model is then formed. This is a discrete matrix representation of the integral of Equation 2. It acts as a summation/integral operator. It is populated from a vector (i.e. a 1N matrix) of TWT times (t.sub.i) of the TWT-RMS velocity pairs of data set 14, and from the desired output TWTs (t.sub.j) 16, as follows:

(44) $\begin{array}{cc}G=\left[\begin{array}{ccccccc}{t}_1/t_1& 0& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& & & & & & \\ t_1/t_i& \mathrm{.Math.}& \mathrm{.Math.}& t_j/t_i& 0& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& & & \mathrm{.Math.}& & & \\ t_1/t_N& \mathrm{.Math.}& \mathrm{.Math.}& t_1/t_i& \mathrm{.Math.}& \mathrm{.Math.}& t_M/t_N\end{array}\right]& \left(\mathrm{Equation}8\right)\end{array}$

(45) This is a discrete version of the forward operator that takes a vector of squared interval velocities defined for layers, j, with two-way-time thicknesses t.sub.j and transforms it into a vector of squared RMS velocities defined at times t.sub.i. If t.sub.i lies within the interval j, then t.sub.j refers to the two-way time from the top of that interval down to t.sub.i only, and not the whole TWT thickness of that interval.

(46) The covariance matrix for the prior model for the data vector .sub.d (i.e. G.sub.mG.sup.T+.sub.e) is evaluated and then matrix inverted. The posterior distribution comprising posterior expectation vector .sub.m|d and posterior covariance matrix .sub.m|d are then readily determined using previously formed vectors and matrices.

(47) Having evaluated the posterior expectation vector .sub.m|d and the posterior covariance matrix .sub.m|d, the non-squared Dix interval velocity and standard deviations are derived. The posterior expectation vector .sub.m|d comprises squared Dix interval velocities for different two-way time intervals j, and the posterior covariance matrix .sub.m|d provides the variance on the squared Dix interval velocities. However, the non-squared Dix interval velocities V are typically of most interest physically. The expectation value, , of the Dix interval velocity V can be estimated from evaluating the equation:
=[E(V.sup.2).sup.2Var(V.sup.2)].sup.1/4(Equation 9)

(48) The corresponding standard deviation .sup.2 of the (non-squared) Dix interval velocity vector can be estimated by evaluating the following equation:
.sup.2=E(V.sup.2).sup.2(Equation 10)

(49) From the standard deviations, confidence intervals for the Dix interval velocities V can be determined. The derivation of the expectation and standard deviations is described further in Appendix A. The Dix interval velocities and standard deviations may then be output in output 22a of the workflow 10.

(50) In a different implementation, all available TWT-RMS velocity pairs and covariances could be input and used to form global matrices of the data and covariances, which can then be inverted to provide an output of Dix interval velocities. The sequence above may then only need to be carried out once.

(51) FIG. 3 shows an example of the output 22a in the form of a graph 30 of velocity against travel time, showing an interval velocity curve 34 derived from the inversion routine 20 alongside a prior art velocity curve 32, at a single x,y location. The interval velocity curve 34 comprises a best-fit velocity curve 34a and associated uncertainties 34b.

(52) The process above is typically repeated for each x,y location where RMS velocity-TWT pairs of data set 14 have been defined. It will also be appreciated that the sequence of formation of the various matrices may be different in other examples. It is sometimes desirable that the Dix interval velocity estimates are converted to define depth intervals. This may be user-selectable as an input option, or may be done automatically. For example, the program could be designed to output both if desirable.

(53) Selectable Input/Output Parameterisation

(54) In a variation of the 1D inversion method described above, the TWT points of the set of TWTs 16 defining the intervals at which the output is to be determined, and/or the TWTs of the prior model 18, are selectable between three different pre-defined parameterisations, to produce a model output: 1) At the same points in TWT as the input RMS velocities for each x,y location; 2) At an arbitrarily chosen set of points in time; or 3) as a densely sampled or pseudo-continuous function of TWT. The parameterisation required can be selected as an extra selection input 19.

(55) These parameterisation options of the selection input 19 seek to provide complete generality of input and output. Any number of data samples and prior model parameters for any TWT position can be used. Similarly, interval velocities can be evaluated at any set of points and/or intervals in terms of TWT or depth, not necessarily corresponding to the same TWT points specified in the input data set 14.

(56) The modelling results shown in FIG. 3 are based on the first parameterisation where N=M, i.e. the number of interval velocities is the same as the number of RMS velocities, and each depth/TWT point at which the RMS velocities are provided defines a corresponding top/bottom of the interval over which the interval velocity is determined. FIG. 4 shows similar results to those of FIG. 3, but demonstrates the second parameterisation in which the number of interval velocities is fewer than the number of RMS velocities. Features in FIG. 4 are denoted with the same reference numerals as corresponding or similar features in FIG. 3, but incremented by 10.

(57) It should be appreciated that in other embodiments, a global 2D or 3D inversion may be performed instead of the 1D method above. In this case, the method will be performed in all x,y locations together. In such embodiments, given sufficient computational power, all available RMS velocities could be inverted simultaneously in a single, global, Dix inversion, to predict Dix interval velocities and covariances as a function of arbitrary locations in x,y,t or x,y,z. Alternatively, a 2D or 3D inversion may be performed in a limited range of x,y locations as a sub-set of the whole 3D volume. For example, a cube of the full inversion volume may be defined and the inversion may be performed in respect of that cube, and then the cube may be moved sequentially to new lateral positions until estimates for the full volume have been obtained. This may be applied to areas of influence around lateral positions where correlations tend to zero at long distances.

(58) Smoothed Output in TWT/Depth

(59) In a further variation, where a pseudo-continuous output may be specified, the number of specified interval velocities is much larger than the number of picked RMS velocities such that the output of interval velocities is initially represented as a relatively finely sampled discrete series of points with TWT or depth.

(60) When such pseudo-continuous output is specified this is ensured by adding a requirement of smoothness within the inversion itself. This is done by means of a correlation function. This allows the output of stable Dix interval velocities from the inversion with a finer temporal sampling than would be possible without the imposition of such a correlation function (see smooth output option of output 22a).

(61) Many different correlation functions are suitable to achieve this. For example, the correlation function could be an exponential correlation function:

(62) $\begin{array}{cc}C\left(\right)=\exp \left[-3{\left(\frac{}{R}\right)}^v\right]& \left(\mathrm{Equation}11\right)\end{array}$

(63) where

(64) is the time difference between two output points

(65) R is correlation range parameter

(66) v is order of the correlation function between 1 and 2

(67) In this case, the correlation range and correlation order need to be specified. These parameters are user-specified as part of the selection input 19 for the pseudo-continuous parameterisation. The correlation function is imposed within the inversion by modifying the model covariance matrix .sub.m. This results in a model covariance matrix .sub.m, which is similar to .sub.m, except it now includes non-zero off-diagonal elements which define the covariance of parameters.

(68) FIG. 5 shows interval velocity results for a continuous output using the function above with correlation range R=300 ms and v=1.8. Features in FIG. 5 are denoted with the same reference numerals as similar features in FIG. 3, but incremented by 20.

(69) Prediction Procedure

(70) The method is provided with further flexibility to estimate interval velocities and the associated standard deviations in any x,y,z or x,y,t location based on data from arbitrary, and sometimes sparsely available, x,y,t locations. This allows the interval velocity estimates to be regularised, and mapped onto a grid at a number of finely spaced locations, not only in terms of depth or TWT as described so far, but also in x and y directions, thereby forming a prediction model of the subsurface region.

(71) To do this, a separate kriging or prediction procedure 24 is performed using estimated interval velocities and the associated uncertainties at different x,y,t or x,y,z points, for example, as previously determined via the Dix inversion procedure described above. However, it will be appreciated that interval velocity estimates could be derived from other sources, where each velocity has been assigned an uncertainty, for use in the prediction procedure 24. An appropriate prior model would also need to be determined in this case.

(72) The prediction procedure 24 is also carried out in this example by an implementation of Bayes' Rule, in a similar way to the inversion routine 20 as described above but it will be appreciated that a constrained least squares inversion could be used as an alternative. In this case, interval velocity profiles (i.e. interval velocities with TWT/depth for a particular x,y location) and associated uncertainties from the output 22a of the inversion 20 now act as input data in the form of initial velocity profiles for the prediction procedure 24. The forward model in this procedure can be formulated as:
d=v+e(Equation 12)

(73) where

(74) d is a vector of imperfect observations of the velocity

(75) v is a model vector of interval velocities

(76) e is an error vector

(77) The observations d are the interval velocities estimated from the Dix inversion procedure and v is a vector which contains in each element the desired velocities to be estimated at each lateral output location. The prior distribution of interval velocities v is Gaussian in a similar way to the velocity inversion procedure described further above. The model matrix of interval velocities thus has a Gaussian prior distribution that has a mean vector .sub.v and covariance matrix .sub.v. The distribution of the error e is also Gaussian with zero expectation and a covariance matrix .sub.e.

(78) A correlation function used in constructing the prior covariance matrix constrains the smoothness of a final prediction model solution output 22b containing Dix interval velocity profiles at any x,y,t or x,y,z location. In the presently described example, the correlation function allows the initial profiles to be interpolated smoothly to any x,y,t or x,y,z location, to form new estimates in specified x,y,t or x,y,z locations along with uncertainties. One possible correlation function is:

(79) $\begin{array}{cc}C\left(r,z\right)=\exp \left[-3{\left(\frac{r}{R_r}\right)}^{v_r}\right]\exp \left[-3{\left(\frac{z}{R_z}\right)}^{v_z}\right]& \left(\mathrm{Equation}13\right)\end{array}$

(80) where

(81) .sub.r and .sub.z are the lateral and vertical distance or separation between two locations

(82) R.sub.r and R.sub.z are lateral and vertical correlation range parameters v.sub.r and v.sub.z are orders of the lateral and vertical correlation functions

(83) This function provides a correlation structure for the prior velocity model. It also requires correlation order and length to be specified. These parameters are typically specified as an extra input in the selection input 19. The model results may be mapped into a regular x,y,z or x,y,t grid if desired.

(84) The solution is obtained by evaluating a posterior distribution in terms of the posterior expectation .sub.0|d for the velocity in the new arbitrary position x.sub.0 and its variance .sub.0|d.sup.2 according to the following equations:
.sub.0|d=.sub.0|.sub.0u(.sub.r|.sub.e).sup.1(d.sub.v),(Equation 14)
and
.sub.0|d.sup.2=.sub.0.sup.2.sub.0u(.sub.u+.sub.e).sup.1.sub.0v.sup.T.(Equation 15)

(85) where

(86) .sub.0 is an expected interval velocity parameter of the original prior model

(87) d is estimated interval velocities from the inversion

(88) .sub.u is the expected interval velocity of the prior model for the prediction

(89) .sub.v is the covariance matrix for the prior model for the prediction

(90) .sub.0v is the matrix transpose of .sub.v0 a modified covariance matrix formed from the covariance matrix of the original prior model .sub.0 multiplied by the correlation matrix for the correlation function C

(91) .sub.0 is a standard deviation of the original prior model

(92) In the present example of lateral prediction, evaluation of these equations 14 and 15 is carried out in the following steps. The initial interval velocity profiles (output from inversion 20) are contained in the vector d. In order then to perform the prediction to arbitrary x,y locations, the desired output grid onto which the velocity profiles are to be estimated is first determined, and then the separation between each pair of desired x,y grid locations is calculated. Correlation length and order parameters for the correlation function are read in from the selection input 19. Using the correlation range and order and the separations between locations, the correlation function above is evaluated, and the correlation values are used to form the modified covariance matrix .sub.v0 and its transpose .sub.0v. This matrix provides a model of the covariance between a velocity V.sub.0 to be predicted in a new arbitrary position x.sub.0 and the observed velocities of d. Specifically, it incorporates the lateral spatial correlation according to the function of Equation 13.

(93) The prior model for the prediction procedure is then formed. This prior model is based on a combination of the original prior model 18 as used in the inversion 20, and the estimated model of interval velocities (including uncertainties) from the inversion output 22a. The parameter .sub.0 represents an interval velocity value of the original prior model 18 at a particular depth or TWT location. The prior model for the prediction procedure is then defined by an expectation .sub.v, which is the value of the interval velocity of the original prior model 18 evaluated at desired output positions, and by an uncertainty associated with each interval velocity which is contained in the covariance matrix for the prior model .sub.v.

(94) The equations above can then be evaluated, providing as a result the values of interval velocity .sub.0|d in the predicted positions with an associated variance .sub.0|d.sup.2.

(95) Finally, a correction is applied to the values of interval velocity and variance to take account of using the (original) prior model twice (first in constraining the inversion and then in the prediction procedure). The corrected values are given respectively by these relationships:

(96) $\begin{array}{cc}d={}_v+\frac{\left({}_v^2+{}_e^2\right)}{{}_v^2}\left({}_{v|d}-{}_v\right)& \left(\mathrm{Equation}16\right)\\ {}_e^2=\frac{\left({}_v^2{}_{v|d}^2\right)}{{}_v^2-{}_{v|d}^2}& \left(\mathrm{Equation}17\right)\end{array}$

(97) The results from Equations 16 and 17 represent the final output 22b from the prediction procedure 24. An example set of such results is shown in a graph 100 of FIG. 6 for a depth slice at 1500 m. The y-axis represents velocity whilst the x-axis represents lateral distance. The graph 100 shows a smooth curve 110 comprising a best-fit velocity 110a and its uncertainty 110b having applied the prediction procedure. For comparison, it also shows a curve 112 containing the raw interval velocities 112a and uncertainties 112b from the inversion.

(98) Although FIG. 6 is in the form of a 2D graph, typically the interval velocity and uncertainty distribution is provided in grid form, for example in x,y,z or x,y,TWT positions, suitable for providing a 2D or 3D spatial image of the velocity structure of the earth's subsurface.

(99) The derivation of the numerical implementation of the prediction procedure is described in further detail in Appendix B.

(100) In general it is preferable to implement the present method to provide maximum flexibility for a user in specifying and handling the inputs 12 and outputs 22a, 22b. For example, it may be desirable for a user to be able to easily incorporate prior information from different sources and it may be desirable to be able to freely specify prior model interval velocity values, uncertainties attached to one or more of those values, the desired x,y,t and optionally x,y,z output locations for the inversion and prediction procedure. Further, it may be designed to accommodate RMS velocity data and uncertainties in various formats as may be produced from different industrial software packages. It may also allow manual or automatic input of some or all of such data. In addition, it may allow selection of the smooth or discrete output, and control of correlation parameters for smoothing vertically or in TWT in relation to the velocity inversion and in the prediction procedure.

(101) It may also be desirable to select the nature of calculation for example to perform a global velocity inversion of RMS velocities for a large range of x,y locations simultaneously in parallel, rather than perform successive inversions in sequence for different x,y locations as described in examples above.

(102) The present invention provides significant advantages, and the flexibility of input and output is particularly beneficial. It allows inclusion of information about the expected distribution of interval velocity with TWT from any source, and is not limited to a set of particular functional forms. Output of velocities onto an arbitrary grid makes it easier to obtain structure-consistent interval velocities even when input RMS velocities have not been picked along the same time horizons that define the structure. Another advantage is the flexibility of time sampling provided by the different output parameterisations from the inversion. Moreover, specifying uncertainties and prior knowledge as inputs provides a simpler, clearer and more intuitive way to constrain the inversion.

(103) A particularly important advantage is that the method includes a priori information and estimated uncertainty in the inversion in a statistically correct and flexible way. It is therefore possible to use the measures of uncertainty provided in the inputs and the outputs to better understand the reliability of estimated interval velocities. This applies both to x,y,TWT or x,y,z images of velocity distribution.

(104) As a user can specify a prior model and its uncertainty and estimate uncertainties in the input data, the effects of the prior information and velocity uncertainty on the model estimation can be assessed. In particular, a user can explore such effects on seismic imaging algorithms that use interval velocities (in TWT or depth). Further, where the RMS velocities are obtained from seismic reflection sections, the results of the present method can help to interpret the sources of error contributing to a mis-tie of well data with seismic reflection data.

(105) In the time domain (e.g. TWT output), the inversion results can be used to facilitate interpretation of features, for example to help assess the likelihood of particular structural features being present in a seismic reflection section such as whether there is structural closure above a reservoir.

(106) Since there is an uncertainty attached to the velocity, the uncertainty in the geometrical structure seen in a depth converted image can also be assessed. In turn, it is possible therefore to make best- and worse-case calculations of reservoir volumes, based on the structures seen in the depth images, such that the commercial risks can be calculated.

(107) In a similar manner, the results can advantageously be used to estimate probability distributions of other physical properties (e.g., resistivity) of the subsurface which may have a theoretical or empirically-derived relationship with interval or instantaneous velocity. Further, it may also be possible to perform a quantitative assessment of uncertainties on other geological estimates that are based on interval velocities, such as geological uplift estimates.

(108) The Dix interval velocity distribution from the present method may be used in other contexts, for example incorporated in stochastic seismic modelling or stochastic seismic inversion for producing a model of reflectivity or more typically impedance characteristics of the subsurface.

(109) Various modifications may be made without departing from the scope of the invention herein described. For example, it will be appreciated that the method could be applied to form a prior model including previously known surface seismic interval velocities and well constrained vertical velocities from well data in a statistically correct way for depth conversion. Further, the velocity trend of the prior model could be interpolated between well locations with an uncertainty that depends on the distance from the nearest well data location.

(110) It will be appreciated that although a Bayesian inversion technique is described above, the interval velocities may be estimated by a constrained least-squares inversion of the data. In such an embodiment, the expectation vectors and associated covariances used in the examples above are used as inputs to the least-squares inversion. The constrained least-squares inversion is equivalent to the Bayesian inversion only when the prior model, the error and the likelihood model are represented by Gaussian or multi-Gaussian distributions.

Appendix A

(111) This appendix describes the underlying derivation of the posterior distribution from which the estimate of Dix interval velocities with associated uncertainties is obtained.

(112) The simple linear relation between the squared interval velocities and the squared RMS velocities is

(113) $U_n^2=\frac{1}{t_n}\underset{i=1}{\overset{n}{.Math.}}{V}_i^2{t}_i$
and the following defines the relation between squared RMS velocities and the instantaneous velocity

(114) $U^2\left(t\right)=\frac{1}{t}{}_0^t{V}^2\left(t\right)\mathrm{dt}.$

(115) These equations propose a practical choice of data and model parameters. Let the data vector d contain the squared RMS velocities picked from velocity analysis
d[U.sub.1.sup.2,U.sub.2.sup.2, . . . ,U.sub.N.sup.2].sup.T
and let m be the unknown squared interval velocities
m[V.sub.1.sup.2,V.sub.2.sup.2, . . . ,V.sub.M.sup.2].sup.T
where N is the number of picked RMS velocities, and M is the number of interval velocities to be predicted. In classical Dix inversion as defined in

(116) $V_n^2=\frac{U_n^2{t}_n-U_{n-1}^2{t}_{n-1}}{t_n}$
N=M, but we will also allow for estimation of M interval velocities where M is smaller or larger than the number of picked RMS velocities N.

(117) The relation between the squared picked RMS velocities in d and the unknown squared interval velocities in m can be written
dGm

(118) where G is an integral or summation operator with dimension NM,

(119) 0$G=\left[\begin{array}{ccccccc}{t}_1/t_1& 0& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& & & & & & \\ t_1/t_i& \mathrm{.Math.}& \mathrm{.Math.}& t_j/t_i& 0& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& & & \mathrm{.Math.}& & & \\ t_1/t_N& \mathrm{.Math.}& \mathrm{.Math.}& t_j/t_N& \mathrm{.Math.}& \mathrm{.Math.}& t_M/t_N\end{array}\right]$
or simply defined as

(120) $G_{\mathrm{ij}}=\{\begin{array}{cc}\frac{t_j}{t_i}\text{:}& \underset{k=1}{\overset{j}{.Math.}}{t}_k{t}_i\\ 0\text{:}& \underset{k=1}{\overset{j}{.Math.}}{t}_k>t_i\end{array}$
where t.sub.i is the two-way travel time corresponding to RMS velocity U.sub.i, and t.sub.j is the two way travel-time thickness of interval or layer j. If t.sub.i is within interval j, then t.sub.j represents the thickness in two-way time from top of interval j to t.sub.i.

(121) The errors in the picked RMS velocities are assumed to be additive
d=Gm+e
where the error vector e has dimension N1. In general, the RMS velocity uncertainty can be estimated from semblance analysis and gather flattening and stack sensitivity. Usually the uncertainties increase with depth, but the uncertainty may also increase at shallow depths due to the limited offset range available where far offset data may be muted to avoid other, undesirable effects on the data. This form does not explicitly include error in the time t.sub.i of a picked RMS velocity U.sub.i but this can indirectly be accounted for as an error in the RMS velocity. The statistical properties of the error can generally be described by a probability density function (PDF) denoted p(e). The forward model together with the statistical error model define the likelihood model p(d|m).

(122) Geophysical inversion problems are often solved using the least squares approach. In a statistical setting, this is closely related to the Gaussian assumption. If we assume that the error vector e is multi-Gaussian with zero expectation and covariance matrix .sub.e, compactly denoted
eN(0,.sub.e)
then the likelihood function p(d|m) is Gaussian with expectation vector Gm and covariance matrix .sub.e,
d|mN(Gm,.sub.e)
where means distributed as. Note that the covariance matrix allows for varying uncertainty for the different picked RMS velocities, and may also model covariances between the RMS velocity errors. The special case of uncorrelated (white) errors defines a diagonal covariance matrix.

(123) Assume further the prior model p(m) can be represented by a multi-Gaussian distribution with expectation vector .sub.m and covariance matrix .sub.m, denoted
mN(.sub.m,.sub.m)
where the expectation vector represents the velocity depth trend for the M intervals, and the MM covariance matrix contains the prior variances on the diagonal representing the prior trend uncertainty, and the off-diagonal elements represent spatial covariance.

(124) In this linear Gaussian case, an elegant way of finding the posterior distribution is to define the Gaussian joint distribution for m and d and then use a standard result for the conditional distribution for m given d, see text books on multivariate statistics (e.g. Anderson, T. W., 1984, An introduction to multivariate statistical analysis: John Wiley & Sons Inc). In general, a prior distribution for the data vector (before the data are available) can be derived from the likelihood model and the prior model by the integral:
p(d)=p(d|m)p(m)dm

(125) In the linear Gaussian case, the prior model for the data vector is also Gaussian
dN(G.sub.m,G.sub.mG.sup.T+.sub.e)

(126) If we assume that the error vector e and the model parameter vector m are stochastic independent, the covariances are
Cov(d,m)=Cov(Gm+e,m)=G.sub.m
and
Cov(m,d)=Cov(m,Gm+e)=.sub.mG.sup.T

(127) The joint distribution of the model parameter vector m and the data vector d can now be written:

(128) $\left[\begin{array}{c}m\\ d\end{array}\right]N\left(\left[\begin{array}{c}{}_m\\ G{}_m\end{array}\right],\left[\begin{array}{cc}{}_m& {}_m{G}^T\\ G{}_m& G{}_m{G}^T+{}_e\end{array}\right]\right)$
and the posterior distribution is defined by the conditional Gaussian distribution for m given d:
m|dN(.sub.m|d,.sub.m|d)
where the posterior expectation vector is
.sub.m|d=.sub.m+.sub.mG.sup.T(G.sub.mG.sup.T+.sub.e).sup.1(dG.sub.m)
and the posterior covariance matrix is
.sub.m|d=.sub.m.sub.mG.sup.T(G.sub.mG.sup.T+.sub.e).sup.1G.sub.m.

(129) The complete solution for the squared interval velocities is defined by a multi-Gaussian distribution with explicit expressions for the posterior expectation vector and posterior covariance matrix. The posterior expectation represents a stabilized optimal solution in the least squares sense, and the covariance matrix contains the uncertainties and parameter covariances. Since the vector of squared interval velocities is multi-Gaussian, the marginal distribution for the squared interval velocity is a specific layer is also Gaussian. This means that the marginal distribution for the squared interval velocity in layer i, p(V.sub.i.sup.2), is completely described by element i in the posterior vector .sub.m|d and the corresponding variance element on the diagonal of the posterior covariance matrix .sub.m|d.

(130) The forward model and the Dix inversion above is defined for a system of squared RMS velocities and squared interval velocities, but usually we are interested in the interval velocity itself. The distribution of the interval velocity V and the squared velocity V.sup.2 cannot both be Gaussian, but if the standard deviation is small relative to the expectation, this is a good approximation. Assume now that V has a Gaussian distribution with expectation and standard deviation ,
VN(,.sup.2)
and assume also that >0 and >> such that we can neglect the probability that V<0.

(131) A standard normal form is obtained by subtracting the expectation value and dividing by the standard deviation:

(132) $Z=\frac{V-}{}N\left(0,1\right)$$\mathrm{and}$$V=+ZN\left(,{}^2\right)$

(133) The squared velocity is now
V.sup.2=.sup.2+2Z+.sup.2Z.sup.2.

(134) The squared standard normal random variable is chi-squared distributed with one degree of freedom, Z.sup.2.sub.1.sup.2, with expectation 1 and variance 2. We then get that the expectation of V.sup.2 is
E(V.sup.2)=.sup.2+.sup.2
and the variance is
Var(V.sup.2)=4.sup.2.sup.2+2.sup.4.

(135) Given that we know the expectation and variance of V.sup.2, the inverse relations for the expectation and variance of V are
=[E(V.sup.2).sup.2Var(V.sup.2)].sup.1/4,
and
.sup.2=E(V.sup.2).sup.2.

(136) A total vertical velocity profile can be represented by a Gaussian vector v with posterior expectation vector and covariance matrix .

Appendix B

(137) This appendix describes the mathematical derivation underlying the prediction procedure used to interpolate the inversion results onto an arbitrary grid in x,y,t or x,y,z.

(138) We assume that the velocity can be modelled as a Gaussian random field. The prior velocity model is defined by the expectation function (x) and the variance function .sup.2(x), where x is a spatial location coordinate, and a spatial correlation function C(x,y) between the velocities in two positions x and y. The spatial correlation is typically defined by the distance between the two positions, and such that the correlation range is long laterally along the stratigraphy, and shorter vertically. A simple spatial correlation function of an exponential form can be written

(139) $C\left(r,z\right)=\exp \left[-3{\left(\frac{r}{R_r}\right)}^{v_r}\right]\exp \left[-3{\left(\frac{z}{R_z}\right)}^{v_z}\right]$

(140) where r and z are the lateral and vertical distance between two locations, R.sub.r and R.sub.z are the lateral and vertical correlation range parameters, and v.sub.r and v.sub.z are orders of the lateral and vertical correlation functions, respectively. The variance and correlation functions define the spatial covariance function (x,y).

(141) Instead of performing a 3D Dix inversion using all RMS velocity data simultaneously, we here propose a practical two step procedure: First independent 1D Bayesian Dix inversion (as per Appendix A) in all positions with RMS velocity data, and then spatial prediction of velocities in all positions in an arbitrarily sampled cube.

(142) Let us now consider the separate problem of spatial prediction of velocities in arbitrary locations x.sub.0 from some sparse observations of the velocity. Assume that we have imperfect observations d of the velocities v in the positions x.sub.1, x.sub.2, . . . , x.sub.n,
d=v+e,

(143) Where d=[d(x.sub.1),d(x.sub.2), . . . , d(x.sub.n)].sup.T is the imperfect observations of the velocities, v=[V(x.sub.1),V(x.sub.2), . . . , V(x.sub.n)].sup.T is the velocities, and e is an error term. From the assumptions above, the velocities have a Gaussian prior distribution,
vN(.sub.v,.sub.v),
where .sub.v and .sub.v are the prior mean and covariance. The error vector is assumed to have a Gaussian distribution
eN(0,.sub.e).

(144) We now want to predict the velocity V.sub.0 in a new position x.sub.0 with prior expectation p.sub.0 and variance .sub.0.sup.2. Under the Gaussian assumptions above we can now define the joint distribution

(145) $\left[\begin{array}{c}{V}_0\\ d\end{array}\right]N\left(\left[\begin{array}{c}{}_0\\ {}_v\end{array}\right],\left[\begin{array}{cc}{}_0^2& {}_{0v}\\ {}_{0v}^T& {}_v+{}_e\end{array}\right]\right).$
where .sub.0v is the covariance between V.sub.0 and the observed velocities in vector d, and .sub.0v=Cov(V.sub.0, v), where we assume independence between v and e. According to standard rules of multivariate analysis, we now end up with a Gaussian posterior distribution for the predicted velocity V.sub.0,
V.sub.0dN(.sub.0|d,.sub.0|d.sup.2),
with expectation
.sub.0|d=.sub.0+.sub.0v(.sub.v+.sub.e).sup.1(d.sub.v),
and variance as
.sub.0|d.sup.2=.sub.0.sup.2.sub.0v(.sub.v+.sub.e).sup.1.sub.0v.sup.T.

(146) Based on the observed velocities in positions x.sub.1, x.sub.2, . . . , x.sub.n (contained in vector d), the spatial prediction algorithm can now be used repeatedly to predict the velocity in any arbitrary spatial location x.sub.0. A typical use will be to fill a densely sampled 3D model, where each cell is defined by the predicted velocity and the associated uncertainty.

(147) Using the Dix velocities as the observed velocities in the spatial prediction algorithm now proposes itself, where the posterior mean and variance from the Bayesian Dix inversion appear as natural candidates for the velocity observations and the error. However, the use of the Bayesian Dix velocities directly will cause double use of parts of the prior information, first in the Dix inversion and next in the spatial prediction. In the spatial prediction, the prior model used in the Bayesian Dix inversion is extended to include lateral correlations to impose lateral smoothness, compare the correlation function used in Appendix B

(148) $C\left(r,z\right)=\exp \left[-3{\left(\frac{r}{R_r}\right)}^{v_r}\right]\exp \left[-3{\left(\frac{z}{R_z}\right)}^{v_z}\right]$
with the correlation function used in the main example of the velocity inversion

(149) $C\left(\right)=\exp \left[-3{\left(\frac{}{R}\right)}^v\right].$

(150) It is however possible to wash away the effect of the prior information used in both steps.

(151) For a large data set, the vector d containing the velocity data may be large. However, since the spatial correlations may go to zero at large separations, it may be a good approximation to calculate the velocity V.sub.0 using only Dix velocities with a short lateral separation and a significant correlation. Since the Dix velocities are vertically correlated, the value of including values above and below the depth plan defined by x.sub.0 is low and thus not used. Let the prior and posterior distribution in a location x.sub.i be N(.sub.v, .sub.v.sup.2) and N(.sub.v|d, .sub.v|d.sup.2) respectively, where index i is suppressed to simplify the notation. We then compute the Dix velocity in location x.sub.i without the effect of the prior

(152) $d={}_v+\frac{\left({}_v^2+{}_e^2\right)}{{}_v^2}\left({}_{v|d}-{}_v\right)$
and the corresponding Dix uncertainty

(153) ${}_e^2=\frac{\left({}_v^2{}_{v|d}^2\right)}{{}_v^2-{}_{v|d}^2}.$

(154) The values computed in these expressions for d and .sub.e.sup.2 can be inserted into d and .sub.e for location i in the expressions for the expectation .sub.0|d and .sub.0|d.sup.2 above. This defines a two step procedure of Bayesian Dix and spatial prediction where the repeated use of the prior information is accounted for.

(155) It will be appreciated that this prediction procedure could equally and alternatively be applied to results from the 1D inversion in terms of two-way travel time rather than depth. In the case of applying the technique to velocity depth profiles, a conversion is carried out, for example as part of the inversion routine, to transform two-way times to depths. Data may be converted from two-way seismic travel time t to depth z by the transform
z(t)=.sub.0.sup.eV(f)dt,
where V is the estimated Dix interval velocity (posterior expectation).