WO2011039149A1

WO2011039149A1 - Improved estimation of time shift based on multi-vintage seismic data

Info

Publication number: WO2011039149A1
Application number: PCT/EP2010/064288
Authority: WO
Inventors: Espen Oen Lie
Original assignee: Statoil Asa
Priority date: 2009-09-30
Filing date: 2010-09-27
Publication date: 2011-04-07
Also published as: GB0917141D0; GB2474022A; GB2474022B

Abstract

The invention is a method of calculating seismic time shifts Δt(t), comprising : - conducting a number (n) of three or more seismic surveys with large time differences such as months or years, producing a number (n) of three or more vintages of seismic data d₁(t), d₂(t),.., d_n(t)= d (t), each said vintage of seismic data dx(t) comprising one or more seismic traces, - calculating partial time derivatives from one or more base trace ∂/∂t d_b(t) multiplied by the sought time shift Δt(t) and setting the at product like ∂/∂t d_b (t) Δ t(t) = Δd, simplified to LΔt(t) = Δd wherein L comprises the partial time derivatives, and forming a set of (n-1) equations ĽΔt =Ad which is overdetermined in the series of consecutive time shifts Δt (t), - and solving for the consecutive time shifts Δt (t).

Description

IMPROVED ESTIMATION OF TIME SHIFT BASED ON MULTI-VINTAGE SEISMIC DATA

Introduction

The present invention relates to seismic data processing.

More specifically, the invention relates to a method of

calculating seismic time shifts between vintages of seismic data. The purpose of the method of the invention is to obtain an

improved measurement of time shift in the seismic data in order to interpret such time shift to relate to reservoir parameters which may change during production of petroleum fluids. The time shifts may be used during petroleum production activities. The time

shifts may be used for monitoring production from subsurface

geological reservoirs. The time shifts may be used, for example along with other petroleum reservoir data, as input data for

controlling further production from said subsurface reservoirs, such as controlling fluid injection rates, gas production rates, petroleum production rates, changing production positions, or setting new production or injection wells.

For finding petroleum fluids in geological strata the seismic method is the prime method. Seismic signals are generated at the surface and propagate downward and are partly reflected by every seismic impedance contrast. Seismic impedance is the product of seismic acoustic velocity and density. The seismic signals are acquired by a series of seismic sensors after having been

reflected, and the time series collected at a seismic sensor for each seismic transmission from the seismic source is called a seismic trace. For monitoring or controlling the development of the fluid content of the geological strata during petroleum fluid production so-called time-lapse seismic data are conducted during the life of the petroleum field. Material changes within the

geological strata may cause changes of the local seismic impedance and may be seen as a time shift between seismic data acquired at different times during petroleum production. Knowing parameters about the material changes of the geological strata may provide key information to how to control the petroleum fluid production such as adjusting the production rate of gas or petroleum, adjusting the depth of which petroleum fluids are produced, or determining

injection rates of gases or fluids in order to support the

petroleum fluid production.

For obtaining the seismic data for the method according to the present invention one may conduct a number of three or more seismic surveys of said subsurface geological reservoirs with large time differences such as months or years, producing a number of three or more vintages of seismic data, each said vintage of seismic data comprising a number of seismic traces of registrations of seismic data acquired each 4 milliseconds or so. Each seismic survey

usually comprises millions of seismic samples and it requires huge processing capacity to analyse the seismic data.

The present invention relates to estimation of time shifts in seismic data acquired in different vintages. In other words, the present invention uses seismic data measured over a geological structure or geographical in which the seismic data acquisition is repeated in three or more vintages, each vintage separated in time by, say, one or two years.

Long time spans may separate the vintages of seismic data. In one example, the initial or "base" 3-D seismic data survey may have been conducted for mapping in three-dimensional detail a petroleum prospect as early as in 1989, then a seismic survey may have been conducted in generally the same manner in the same place eight years later in 1997 before production drilling has started in the area, then repeated seismic surveys have been conducted during gas or petroleum production each second year in 2001, 2003, 2005, 2007, 2009.

Two-way propagation reflection time (twt) seismic data are generally collected as in ordinary 2-D or 3-D seismic surveys, such as in marine or land surveys. In the seismic data collected during the years, the geological structures as such may generally not have changed significantly during the time span in question. (An exception to this assumption may be petroleum reservoirs experiencing

significant lithological compaction due to decreased liquid content such as after petroleum production in the Ekofisk field.) However, due to changes in gas or liquid content, particularly in reservoir layers of porous or permeable geological formations, the local seismic velocity, such as seismic p-wave or s- ave velocity, in one or more of the geological layers affected may change due to fluid migration such as production of fluids such oil or gas, or injection of fluids such as gas or water or C0₂.

A change in the seismic velocity in one layer of significant thickness or a consecutive sequence of layers may give rise to a detectable time shift of seismic peaks and troughs as a function of the two way time t in the layer itself. The local time shift will be proportional to the local change in seismic velocity and also

proportional to the thickness of the layer in question. Please notice that a calculated time shift may be far less than the seismic sampling time resolution corresponding to 4 ms normal sampling or 8 ms for sparse sampling. A time shift occurring in one layer will shift all layers below the actual layer the same amount of time shift. Time shift in two or more layers is a cumulative function. If several consecutive layers have a generally even changed seismic velocity they may contribute to a generally increasing time shift as counted from the top to the bottom of the layers in question. Moreover, such a stack of layers will contribute to a constant time shift for the geological layers below. This is similar to the effect of the one single layer inducing a time shift. In general, the time shift in the two-way time trace is cumulative from the upper layer displaying time shift and downwards, increasing or decreasing cumulatively for each new layer having its own time shift.

A time shift At (t) may be calculated on a trace-by-trace basis by subtracting a single, first seismic trace di (t) acquired at a specific position at a first time ti from a single, second seismic trace d₂ (t) acquired at the same specific position at a second time t₂. The dimensions in a broader aspect may be extended so as for a three- dimensional time shift cube or volume of time shifts At(t) may be calculated by subtracting a first set of traces, stacked or not, from a first set of seismic traces di(t) acquired from a specifically positioned cube or volume or region of the earth at a first time ti from a generally equally sampled second cube or volume of seismic traces d₂(t) acquired at the same specific position at a second time t₂. The sampled Earth volume may be represented by e.g. 400 time samples x 250 traces in-line times x 250 traces cross-line totalling 25000000 seismic samples, all samples comprising some noise.

One significant area of use for estimating time shifts is for petroleum fluid production from a given part of the geological

subsurface formations and monitoring the time shifts in that

geological underground as the production process progresses.

Significant changes in the observed or calculated time shifts may be interpreted as changes in the fluid masses, pressure changes in the gas masses, and / or changes of the fluid interface positions in the geological strata or geological traps. The changes in the fluid masses may thus be used for controlling the production processes including adjusting production rates or determining injection rates or even halting injection. Further, changes in the fluid masses may be used in the decision process of determining the position of auxiliary

production drilling.

In Fig. 1 is shown several vintages of seismic data acquired as di in 1989, d₂ in 1997, d₃ in 2001, d₄ in 2003, d₅ in 2005, d₆ in 2007, and d₇ in 2009. Consecutive time shifts At_±(t) may be calculated by using a subtraction method producing e.g. the closest consecutive time shifts Ati-₂(t) between 1989 and 1997, At₂-₃(t) between 1997 and 2001, At₃-₄(t) between 2001 and 2003, up to At₆-₇(t) between 2007 and 2009, all in all (n-1) time shifts out of (n) consecutive data sets. Two simple facts arise: The first fact is that a cumulative time shift calculated by adding all the consecutive time shifts from Ati-₂(t) to At₆-₇(t), should ideally be the same as the total time shift Ati,₇(t) between the first data set di from 1989 and the later data set d₇(t) from 2009. The second fact is that all the seismic measurements contain noise due to general seismic noise during the acquisition and instrumental noise. Thus Ati , At₂ , . . . At_n contain noise and do not as such represent true, intrinsic time shifts At_±.

Brief summary of the invention

The above problems may be remedied by the invention which is a method of calculating seismic time shifts At (t) between vintages of seismic data d(t) for monitoring production from subsurface geological reservoirs, comprising: - conducting a number (n) of three or more seismic surveys of said subsurface geological reservoirs with large time differences such as months or years, producing a number (n) of three or more vintages of seismic data di(t), d₂(t),.., d_n(t)= d (t) , each said vintage of seismic data di(t) comprising one or more seismic trace,

- formulating an operator between up to n(n-l)/2 possible pairs of said seismic data sets d (t) and their corresponding up to n(n-l)/2 time shifts At (t) ;

- solving for (n-1) time shifts At (t) between (n-1) selected pairs of vintages of said seismic data di(t), d₂(t),.., d_n(t) as an

overdetermined system of equations,

- applying one or more of said calculated time shifts Ati_₂(t), At 2 -

3 (t) , .. , At_n-i,n(t) as input data for controlling further production from said subsurface reservoirs .

According to an advantageous embodiment of the method, in the step of selecting the time shifts to be solved, it is preferred to select the (n-1) time consecutive shifts Ati_₂(t), At₂_₃(t),.., At_n-i,_n(t) between consecutive vintages of said seismic data di(t), d₂(t),.., d_n(t) . Other time spans between vintages to select time shifts to be solved for may as well be used without contributing anything new.

In a preferred embodiment of the invention, the method comprises, in the step of formulating said operator between the n(n-l)/2 possible pairs of said seismic data sets d (t) , to

- calculate partial time derivatives from at least one base trace

- formulate the operator between said n(n-l)/2 possible pairs of seismic data sets d (t) and their corresponding n(n-l)/2 time shifts

At (t) as —a— d_b(t)At{t) = Ad , simplified to LAt{t) = Ad wherein L comprises at

the partial time derivatives . In a preferred embodiment of the invention, the method of further comprises, in the step of formulating said operator,

- forming a set of n(n-l)/2 equations LAt = Ad which is overdetermined in said series of time shifts, preferably consecutive, such as Ati_

2 (t) , At₂-3 (t) , .. , At_n-i,n(t) .

According to a preferred embodiment of the invention, each said vintage of seismic data di(t) comprises seismic traces representing a three-dimensional volume of part of the underground geological formations, said resulting series of consecutive time shifts At (t) representing time shifts in said volume part of the underground geological formations.

According to a preferred embodiment of the invention, the method further comprises calculating partial time derivatives — d_b (t) of the

d t

traces by estimating coefficients of basis functions.

According to a preferred embodiment of the invention, the method of estimating coefficients comprises estimating spline function coefficients. The order may be low, such as at least order one or higher, up to four or even five.

If using Laplacian basis functions, coefficients may be very low or rather low, such as four or five.

According to a first aspect of the present invention there is provided a method of determining time shifts At (t) between vintages of seismic data d(t) relating to a subsurface geological formation, the method comprising: providing a number n of three or more vintages of seismic data di(t), d₂(t), .., d_n(t) obtained from a corresponding number n of seismic surveys of the subsurface geological formation carried out at least a plurality of months or years apart, each vintage of seismic data di(t) comprising at least one seismic trace; formulating an operator between a first number of the n(n-l)/2 possible amplitude differences Ad_b,_m(t) between pairs d_b(t), d_m(t) of the seismic data sets and a second number of the corresponding n(n-l)/2 possible time shifts At_b,m(t) , where the second number is arranged to be less than the first number by using a known inter-relationship between at least some of the n(n-l)/2 possible time shifts At_b,_m(t) , thereby forming an overdetermined set of equations; and solving the overdetermined set of equations for the second number of time shifts At_b,_m(t) to determine the time shifts At (t) .

The method may comprise using one or more of the determined time shifts At (t) as input for monitoring or controlling production or further production from a subsurface reservoir, or at least arranging for or causing such use.

Each amplitude difference Ad_b,_m(t) may derive from a difference between a seismic trace in data set d_m(t) and a seismic trace in data set

The first number may be n(n-l)/2, such that the amplitude differences Ad_b,m(t) between all n(n-l)/2 possible pairs d_b(t), d_m(t) of the seismic data sets are used.

The second number may be (n-1), with the (n-1) time shifts At_b,_m(t) being the (n-1) consecutive time shifts Ati,₂ (t) , At₂,3 (t) , .. , At_n-i,_n(t) between consecutive vintages of said seismic data.

The step of formulating the operator may comprise calculating, for each of the first number of amplitude differences Ad_b,_m(t), a partial time derivative of a trace from the pair of seismic data sets d_b(t), d_m(t) forming that amplitude difference Ad_b,_m(t), the partial time derivative being represented as L_b;m.

The partial time derivative may be calculated by estimating

coefficients of basis functions.

The basis functions may comprise spline functions or Laplacian functions .

The order of spline function coefficients may be low such as one to four . The number of Laplacian basis function coefficients may be between two and five.

The partial time derivative for Ad_b,_m(t) may be — A d_b (i) , where d_b(t) is dt

the earlier seismic data set of the pair of seismic data sets d_b(t), d_m(t) .

The set of equations, before use of the known inter-relationship, may be represented as L_b;mAt_b,_m=Ad_b,_m or in matrix form as LAt = Ad where L is a diagonal matrix and At and Ad are vectors. Use of the known interrelationship may involve reducing a dimension of the vector At and a corresponding dimension of the matrix L, and duplicating and/or rearranging at least one of the elements of the matrix L to form matrix L , thereby creating the overdetermined set of equations

L At = Ad .

The method may comprise first determining A = L^TL and B = L^T Ad and then solving for At based on At = A^~lB , which is equivalent to

At = (L^T L)^~l L^T Ad , or alternatively solving for At directly based on

At = {L^TL)-^lL^TAd .

Each said vintage of seismic data di(t) may comprise a seismic trace representing a three-dimensional volume of part of the subsurface geological formation, said resulting time shifts At (t) representing time shifts in said volume part of the subsurface geological

formations .

The known inter-relationship may be that the sum of a plurality of individual consecutive time shifts is equal to the time span from the time span from the start of the first time shift of the plurality to the end of the final time shift of the plurality, for example

^Δί1,3 = ^Δί1,2 + ^Δί2,3 ° ^{r t}X,Y = ^tX,X₊l + ^tX₊l,X₊2 +'" + ^Y-2,Y-l + ^Y-I ^■ The method may comprise acquiring the seismic data sets.

The method may comprise using the determined time shifts during petroleum production activities, or at least arranging for or causing such use.

According to a second aspect of the present invention there is provided a method of determining or detecting a change in a

geophysical property of a subterranean region of the earth, such as material changes in geological strata, comprising performing the method according to the first aspect of the present invention, wherein the seismic data sets are acquired from the region, and further comprising: using the determined time shifts to determine a change in the geophysical property of the region.

The method may comprise using the determined change in the geophysical property during petroleum production activities, or at least arranging for or causing such use.

According to a third aspect of the present invention there is provided a method in which the time shifts determined according to the method of the first aspect of the present invention or the change in the geophysical property determined according to the second aspect of the present invention is/are used during petroleum production activities.

The above-mentioned petroleum production activities may comprise at least one of: monitoring or controlling petroleum fluid production; adjusting the production rate of gas or petroleum; adjusting the depth of which petroleum fluids are produced; and determining injection rates of gases or fluids in order to support the petroleum fluid production .

The subsurface geological formation may comprise a subsurface

geological reservoir. According to a fourth aspect of the present invention there is provided an apparatus for performing a method according to any of the first to third aspects of the present invention.

According to a fifth aspect of the present invention there is provided a program for controlling an apparatus to perform a method according to any of the first to third aspects of the present invention or which, when loaded into an apparatus, causes the apparatus to become an apparatus according to the fourth aspect of the present invention. The program may be carried on a carrier medium. The carrier medium may be a storage medium. The carrier medium may be a transmission medium.

According to a sixth aspect of the present invention there is provided an apparatus programmed by a program according to the fifth aspect of the present invention.

According to a seventh aspect of the present invention there is provided a storage medium containing a program according to the fifth aspect of the present invention.

Short Figure captions

The invention is illustrated in the following drawings, in which:

Fig. 1 shows background art in which time shift series are calculated based on pairs of seismic data vintages such as 1989 and 1997, 1997 and 2001, and so on.

Fig. 2 shows a similar set of seismic data vintages as

illustrated in Fig. 1, and a number of potential time shift data which may be calculated between all of the seismic vintage data sets.

However, calculating all such data sets will be a laborious task increasing with the square of the number n of data sets .

Fig. 3 illustrates in (a) a set of two vintages of one seismic trace, in (b) a difference between the two vintages of the seismic trace, in (c) a partial time derivative of the base seismic trace, and in (d) an estimated time shift, here constant for the actual interval.

Fig. 4 illustrates a seismic vessel at the sea surface conducting a seismic reflection survey. The seismic data collected may represent a volume of the subsurface and two vintages of one seismic trace is illustrated, one of them exhibiting a time shift down in the volume. Such a time shift will usually have to be laterally continuous to neighbouring traces in order to be significant, but here we illustrate only one such partially time shifted trace.

Fig. 5 illustrates results from time shifts using single combinations of vintage seismic data using the current method of estimated basis function coefficients. To the upper right, a

horizontal section is shown illustrating time shifts calculated from single combinations, here between 1999 and 2003, and in the lower right, a vertical section of 400 time samples of the same. To the left is shown a vertical section is shown calculated as the difference between 1985 and 1996 seismic vintages, also using the current method of estimated basis function coefficients. The amplitude scale to the right of each section is in units of samples.

Fig. 6 illustrates corresponding results from time shifts calculated according to the present invention using multi-vintage seismic data for calculating seismic time shifts. To the upper right, a horizontal section is shown illustrating time shifts calculated according to the present invention, also here between 1999 and 2003, and in the lower right, a vertical section of 400 time samples of the same. To the left is a vertical section calculated according to the present invention as the difference between 1985 and 1996 seismic vintages .

Fig. 7 shows a comparison of time shift calculations represented by summing background art time shifts as calculated for consecutive vintages of seismic data, cumulatively from 1985 to 2003 to the lower left, and in one difference calculated by the difference from 1985 to 2003 directly, to the right. In the top right diagram is shown the inconsistency between the two methods, which ideally should be zero, but which is not. This inconsistency of the background art is one of the problems sought relieved.

Fig. 8 shows time shift calculations over a long time span as calculated according to the present invention, in which the time shift data over a long time span from 1985 to 2003 is consistent with the accumulated consecutive time shifts from 1985 to 2003. Embodiments of the invention

Contrary to calculating time shifts based on pairs of seismic vintage data sets, usually consecutive seismic data sets, several vintages of seismically acquired data may produce several possible combinations of time shifts. Having three consecutive data sets di di d₃ may give rise to three time shift sets Ati-₂, At₂-3, Ati-₃ of which

At₂-3 are directly estimated and

Combinatorially, the number of combinations of possible time shift volumes increase as (n-l)n/2:

2 data sets may give 1 combination.

3 data sets may give 3 combinations .

4 data sets may give 6 combinations .

5 data sets may give 10 combinations.

6 data sets may give 15 combinations, and so on.

Thus 6 data sets provide 15 measurements of (6-1) = 5 actual physical time shift parameters. Notice that the time shift parameters are multidimensional and may be of significant size. This provides a redundancy which may be utilized according to the present invention to reduce noise and more precisely calculate the intrinsic time shifts. We may estimate (n-l)n/2 time shifts between all combinations between n data sets. If solving for the (n-1) intrinsic physical time shifts from (n-l)n/2 estimated time shifts, instead of attempting to

estimating all the time shifts, we may reduce the number of volume cubes while utilizing the potential of increasing the signal to noise ratio S/N increasing with the square root of n provided by the large number of measurements di . . d_n .

A method currently used in time shift estimation by the inventor may be illustrated using the trace illustrations of Fig. 3. Fig. 3a describes a base trace di here named d_b , the index "b" indicating base. Further, a later acquired monitor trace d₂ , d₃ , .. di here named d_m , the index "m" for "monitor". Fig. 3d indicates a time shift, here shown constant, of which d_m has been displaced relative to d_b. The amplitude difference Ad between the monitor trace d_m and the base trace d_b is formed simply by subtracting d_b from d_m as shown in Fig. 3b. A time derivative of the base trace

—a

a d_b(t)

t

is shown in Fig. 3c. Now we have the elements for setting up a basic equation for the amplitude difference Ad between the time shifted traces :

d_b(t)At(t) + ε(t) + d_4D(t)

or

Ad =—a d_b At + ε + d_4D

at that is: the amplitude difference Ad(t) between the monitor trace and the base trace may be expressed as the partial time derivative of the base trace multiplied with the time shift At (t) , plus random noise, plus possible new events d_4D(t) at one or more depths, new events which were not in the base trace d_b.

According to the background art, the time shift At (t) has been

calculated discretely as time shift vectors for every point of d_m(t) . This is a highly laborious task which requires vast amounts of

processing time for seismic data. Given the noise of both d_m(t) and d_b(t) the calculated time shift At(t) calculated according to the background art may be highly discontinuous and requires some smoothing techniques to be applied. Such discrete time shifts have also been attempted to be correlated laterally to provide better consistency. A priori, it is not likely that the time shift should be discontinuous; rather, it should be a quite smooth function of two way seismic time t, having only sharp gradients on gas/liquid surfaces which have been displaced.

In a preferred embodiment of the invention, the time shift function At (t) may, instead of being calculated as a time shift vector for every 4 millisecond or 8 millisecond sample, rather be calculated far more efficiently by expressing a time shift function by calculating coefficients of basis functions. Instead of calculating 400 time shift vectors for a 400 sample trace, and conducting smoothing operations, one may calculate a few first coefficients of a selected basis function and obtain an estimated time shift function At_e(t) . The inventor has conducted several experiments on real data set vintages by calculating a few coefficients of basis functions such as Laplacian functions, spline functions, polynomials, etc. Using Laplacian functions a specific set of Laplacian coefficients will have an impact on the entire part of the seismic trace selected. The spline function coefficients of low order may be preferred because a low order spline function need not have a global influence on the entire two-way time range selected, but will have limited time (depth) range. This is a preferred method for providing a smooth estimated time shift function in the method of the invention described below. Other methods for finding the time shift function, such as brute force vector

calculation methods of the background art may be used. It is believed that the use of basis function coefficients significantly reduces the computational load for finding a time shift between traces of

different vintages.

The above method was developed for finding the time shift function between two vintages. The above formula

Ad(t) =—a d_b (t) A t(t) + ε( + d_4D (t)

a t may be rewritten as

Ad = LAt or

LAt = Ad in which L contains the partial time derivatives of the base trace, i.e. — d_b(t) , but not the noise and the possible 4D contribution.

d t The material dimensions of the above two equations need not only comprise one single trace. The material dimensions may be broadened to relate to a seismic section profile or a seismic volume with length along a set of generally parallel seismic survey lines, and width along a cross-line direction, and the time samples along the vertical direction. Please notice, as above, that depth in this context is considered as constant (i.e. depth shift is not allowed) while time shift is allowed.

So L may be a 3-D volume matrix n_t x ni x n_x of 400 samples x 1000 traces inline x 100 lines, for example, please see Fig. 4. In practice only a relevant depth part of the traces is used, not the entire depth of all traces, for finding time shift near a petroleum production relevant depth of the geological layers.

We may express two differences between three consecutive vintages of seismic volume data Li, L₂ expressed in matrix form in an analogue way as the above, as follows:

Here, the left matrix is a diagonal matrix of matrixes L, the time shifts Ati_₂(t) and At₂-3<t) are "observed", i.e. the only two

consecutive time shifts as estimated or calculated between three consecutive seismic data volume vintages, noise inclusive.

We now illustrate the method using three vintages of seismic data d(t) . Due to the fact that we actually have three vintages in the example above, we may also use the difference between the first and third data set:

Written out in full, this would be:

L_{l 3}At_{l 3} = Ad_{l 3}

Physically it can be stated that:

Δ*ι_>3 =Δί_1>2+Δί_2>3 so that the final one of the above three equations can be re-written as :

>3(Δίι_>2+Δί_2>3) = Δί_{1 >}3

Thus we may rewrite the above matrix formulation as

(or LAt = Ad )

which is an overdetermined system with one degree of freedom less. In other words, prior to the substitution of At_{l 3} = At_{l 2} + At_{2 3} , there were n(n-l)/2 equations and the same number n(n-l)/2 of unknowns ( At ) , which in the above example is three equations and three unknowns.

After the substitution, there were a first number n(n-l)/2 of

equations with only a second number (n-1) of unknowns (the n-1 time shifts), which in the above example is three equations and two unknowns; the second number (number of unknowns) is smaller than the first number (number of equations) . This is an overdetermined set of equations . The known or predetermined inter-relationship between the n(n-l)/2 time shifts ( At_{l 3} = At_{l 2} + At_{2 3} ) has been used to create an overdetermined set of equations having n(n-l)/2 equations with only (n-1) unknowns (the n-1 time shifts) . It is noted that not all the n(n-l)/2 possible amplitude differences Ad_{b m} need be selected for use in the formulation; there is still a benefit in using the

predetermined inter-relationship between the n(n-l)/2 time shifts to ensure that the second number (number of unknowns) is smaller than the first number (number of equations) so that an overdetermined set of equations is created. In the above expression LAt = Ad , L can be considered to be an operator between Ad and At , or in other words an operator between a first number of the n(n-l)/2 possible amplitude differences Ad_b,_m(t) between pairs d_b(t), d_m(t) of the seismic data sets and a second number of the corresponding n(n-l)/2 possible time shifts At_b,_m(t) .

In the above expression Li_/2 and L_2;3 contain only the derivatives, and Ati,2 (t) and At_2;3(t) are the intrinsic, noise-free time shift functions we seek, while the Adi,₂(t) , Ad_2;3(t) and Adi,₃(t) are our measured differences between the three vintages of seismic data.

There are some subtle differences which arise from using multiple vintages for calculating time shifts according to the present

invention. In the previous method using only the difference between two vintages, one would have to adapt the time shift estimates by allowing noise, while according to the present invention one may adapt the seismic data by allowing noise in the seismic data. We assume that it is easier to match time shifts assuming that there is noise in the seismic data than allowing systematic errors in the time shift estimates .

There are two main advantages to be obtained using the present invention :

a) An improved signal-to-noise ratio is obtained for the calculated time shifts.

b) The geophysicist needs only to deal with (n-1) calculated time shift volumes (or so-called time shift "cubes") which are consistent, i.e. the calculated time shift volumes satisfy all the seismic data optimally .

The above method of finding all the derivatives in the equations:

by calculating all the local time shifts would be a very laborious task requiring huge amounts of processing time. If, however, using the above-mentioned method according to a preferred embodiment of the invention, calculating coefficients of basis functions for obtaining the derivatives, the time shifts for large multi-vintage seismic data sets may calculated using far fewer processor steps.

The expressions above may express multi-vintage consecutive intrinsic time shift while keeping the above general structure:

L At = Ad

wherein we keep the same constraints as in the two vintage method.

We then calculate

A = L^TL

and

B = L^TAd

We then solve:

At = A^~lB (which is equivalent to At = (L^T L)^~l L^T Ad ) .

This is a well-known procedure to make an overdetermined matrix become square and honouring all input data in a least squares sense.

In the above example, having three vintages of seismic data, we may use something analogous to Pascal's triangle to find the combinations:

At_{l 2} t

At₃ in the given example only three time shifts arise as expected.

When we now build our big L matrix:

we have different sums of matrix products

buildin our new data vector is easier

From these we can calculate At = A B .

The example above used three vintages of seismic data,

For a series of four vintages, the data set is:

and the number of potential time shifts is

In a bigger example, six consecutive vintages may be written as d—

_{? ?} and the matrix of derivatives "big L" may be written as

for being used in the overdetermined set of equations .

These equations needs to be solved in a sensible way for n seismic vintages. At the moment, we use a geometric scheme for building the matrix without exploiting the symmetry of the matrix, which includes some duplicate storage of data and some extra calculations. The calculation time grows quadratically with the number of vintages, and has the same order as conducting calculation of all possible

combinations of time shifts.

Examples of results of the background art and of the invention

horizontal section is shown illustrating time shifts calculated from single combinations, here between 1999 and 2003, and in the lower right, a vertical section of 400 time samples of the same. To the left is shown a vertical section is shown calculated as the difference between 1985 and 1996 seismic vintages, also using the current method of estimated basis function coefficients.

Fig. 6 illustrates corresponding results from time shifts calculated according to the present invention using multi-vintage seismic data for calculating seismic time shifts. To the upper right, a horizontal section is shown illustrating time shifts calculated according to the present invention, also here between 1999 and 2003, and in the lower right, a vertical section of 400 time samples of the same. To the left is a vertical section calculated according to the present invention as the difference between 1985 and 1996 seismic vintages. It is clear that the horizontal section shown in Fig. 6 is smoother than the corresponding horizontal section shown in Fig. 5. But the feature in the horizontal section between 130 to 180 metres in one direction and between 50 and 150 metres in the other direction in Fig. 6 is consistent with the corresponding feature shown in Fig. 5. Further, the features shown in Fig. 6 in the vertical sections also correspond with what is shown in Fig. 5, but Fig. 6 shows generally less variation than what is shown in Fig. 5. The smoother results according to the method shown in the images of Fig. 6 is attributed to the general noise-reducing ability of the method of the invention.

Fig. 8 shows time shift calculations over a long time span as calculated according to the present invention, in which the time shift data over a long time span from 1985 to 2003 is consistent with the accumulated consecutive time shifts from 1985 to 2003.

Comparing the right diagram in the lower left part of Fig. 7, which is the long term time shift from 1985 to 2003 calculated according to the background art using single combinations, with the corresponding right diagram of Fig. 8 which is the long term time shift as calculated according to the present invention, they are in good correspondence, i.e. their differences are rather small, as would be expected. The left diagram of Fig. 8 which may represent the accumulated data is of course consistent with the long term time shift shown to the right .

An advantage of the invention is that our time shift calculation method is expanded for use with multi-vintage data and is still rather fast. The method is significantly more memory demanding, but improves resolution and provides consistency to the results as shown above. The method may also improve the signal to noise ratio.

On lateral constraints

We may expect that in order for time shifts to be significant, they should show some lateral continuity, except for across geological faults. In an advantageous embodiment of the invention one may emphasize lateral continuity by, if considering one time trace in one geographical location, adding a weighted sum of the near neighbour traces. In a preferred embodiment, the four nearest neighbour traces would be used. The preferred embodiment works even if though faults are present, because if one of the neighbour traces is actually across a fault, the contribution may be only 1/4 and it will not dominate.

On nodes and grids

Calculating a few coefficients of basis functions may be based on as Laplacian functions, spline functions, polynomials, etc. The spline function coefficients such as bsplines, which at the moment are preferred, requires a set of nodes to be defined. These nodes specify the points where time shifts are discontinuous in some selected derivative. We may use a regular grid for the nodes which is bound to regular sample intervals of the seismic data, but in an embodiment of the invention we may arrange some of the nodes so as for following known interfaces of the subsurface reservoir or geological interfaces, that is, interfaces of which time shift changes typically should change behaviour.

By knowing the materially caused time shifts and their depths, the change of seismic impedance may be calculated. The change in seismic impedance may be either due to a change in seismic velocity or density, or both, of the geological layers in question. Knowing such parameters about the material changes of the geological strata may provide key information input to how to control the petroleum fluid production such as adjusting the production rate of gas or petroleum, adjusting the depth of which petroleum fluids are produced, or determining injection rates of gases or fluids in order to support the petroleum fluid production.

It will be appreciated that an apparatus can provided to perform a method as described above. Such an apparatus can be controlled to perform a method as described above by way of a computer program.

Such a computer program can be stored on a computer-readable medium, or could, for example, be embodied in a signal such as a downloadable data signal provided from an Internet website. The appended claims are to be interpreted as covering an operating or computer program by itself, or as a record on a carrier, or as a signal, or in any other form .

Embodiments of the present invention have been described with

particular reference to the examples illustrated. However, it will be appreciated that variations and modifications may be made to the examples described within the scope of the present invention.

Claims

1. A method of determining time shifts At (t) between vintages of seismic data d(t) relating to a subsurface geological formation, the method comprising:

providing a number n of three or more vintages of seismic data di(t) , d₂(t) , .., d_n(t) obtained from a corresponding number n of seismic surveys of the subsurface geological formation carried out at least a plurality of months or years apart, each vintage of seismic data di(t) comprising at least one seismic trace;

formulating an operator between a first number of the n(n-l)/2 possible amplitude differences Ad_b,_m(t) between pairs d_b(t), d_m(t) of the seismic data sets and a second number of the corresponding n(n- l)/2 possible time shifts At_b,_m(t) , where the second number is arranged to be less than the first number by using a known inter-relationship between at least some of the n(n-l)/2 possible time shifts At_b,_m(t) , thereby forming an overdetermined set of equations; and

solving the overdetermined set of equations for the second number of time shifts At_b,_m(t) to determine the time shifts At (t) .

2. A method as claimed in claim 1, comprising using one or more of the determined time shifts At (t) as input for monitoring or

controlling production or further production from a subsurface reservoir, or at least arranging for or causing such use.

3. A method as claimed in claim 1 or 2, wherein each amplitude difference Ad_b,_m(t) derives from a difference between a seismic trace in data set d_m(t) and a seismic trace in data set d_b(t) .

4. A method as claimed in claim 1, 2 or 3, wherein the first number is n(n-l)/2, thereby using the amplitude differences Ad_b,_m(t) between all n(n-l)/2 possible pairs d_b(t), d_m(t) of the seismic data sets.

5. A method as claimed in any preceding claim, wherein the second number is (n-1) , with the (n-1) time shifts At_b,_m(t) being the (n-1) consecutive time shifts Ati,₂ (t) , At₂ , 3 (t) , .. , At_n-i,_n(t) between consecutive vintages of said seismic data.

6. A method as claimed in any preceding claim, wherein the step of formulating the operator comprises calculating, for each of the first number of amplitude differences Ad_b,_m(t), a partial time derivative of a trace from the pair of seismic data sets d_b(t), d_m(t) forming that amplitude difference Ad_b,_m(t), the partial time derivative being represented as L_b,_m.

7. A method as claimed in claim 6, wherein the partial time

derivative is calculated by estimating coefficients of basis

functions .

8. A method as claimed in claim 7, wherein the basis functions comprise spline functions or Laplacian functions.

9. A method as claimed in claim 8, wherein the order of spline function coefficients is low such as one to four.

10. A method as claimed in claim 8, wherein the number of Laplacian basis function coefficients is between two and five.

11. A method as claimed in any one of claims 6 to 10, wherein the partial time derivative for Ad_b,_m(t) is —— d_b (t) , where d_b(t) is the

a t

earlier seismic data set of the pair of seismic data sets d_b(t), d_m(t) .

12. A method as claimed in any one of claims 6 to 11, wherein the set of equations, before use of the known inter-relationship, is

represented as L_b;mAt_b,_m=Ad_b,_m or in matrix form as LAt = Ad where L is a diagonal matrix and At and Ad are vectors, and wherein use of the known inter-relationship involves reducing a dimension of the vector At and a corresponding dimension of the matrix L, and duplicating and/or rearranging at least one of the elements of the matrix L to form matrix L , thereby creating the overdetermined set of equations L At = Ad .

13. A method as claimed in claim 12, comprising first determining

A = L^TL and B = L^T Ad and then solving for At based on At = A^~lB , which is equivalent to At = (L^T L)^~l L^T Ad , or alternatively solving for At directly based on At = (L^T L)^~l L^T Ad .

14. A method as claimed in any preceding claim, wherein each said vintage of seismic data di(t) comprises a seismic trace representing a three-dimensional volume of part of the subsurface geological

formation, said resulting time shifts At (t) representing time shifts in said volume part of the subsurface geological formations.

15. A method as claimed in any preceding claim, wherein the known inter-relationship is that the sum of a plurality of individual consecutive time shifts is equal to the time span from the time span from the start of the first time shift of the plurality to the end of the final time shift of the plurality, for example At_{l 3} = At_{l 2} + At_{2 3} or ^χ,γ ⁼ ^χ,χ+ι ^~*^~ ί_{χ+ι χ+2} ^"I I^" At_Y__{2 Y}_y + At_Y__{l Y} .

16. A method as claimed in any preceding claim, comprising acquiring the seismic data sets.

17. A method of determining or detecting a change in a geophysical property of a subterranean region of the earth, such as material changes in geological strata, comprising performing the method according to any preceding claim, wherein the seismic data sets are acquired from the region, and further comprising: using the determined time shifts to determine a change in the geophysical property of the region .

18. A method as claimed in any preceding claim, comprising using the determined time shifts of claim 1 or the determined change in the geophysical property of claim 17 during petroleum production activities, or at least arranging for or causing such use.

19. A method in which the time shifts determined according to the method of any one of claims 1 to 16 or the change in the geophysical property determined according to the method of claim 17 is/are used during petroleum production activities.

20. A method as claimed in claim 18 or 19, wherein the petroleum production activities comprise at least one of: monitoring or

controlling petroleum fluid production; adjusting the production rate of gas or petroleum; adjusting the depth of which petroleum fluids are produced; and determining injection rates of gases or fluids in order to support the petroleum fluid production.

21. A method as claimed in any preceding claim, wherein the

subsurface geological formation comprises a subsurface geological reservoir .

22. An apparatus for performing a method according to any preceding claim .

23. A program for controlling an apparatus to perform a method as claimed in any one of claims 1 to 21, optionally being carried on a carrier medium such as a storage medium or a transmission medium.

24. A storage medium containing a program as claimed in claim 23.

25. A method of calculating seismic time shifts At (t) between vintages of seismic data d(t) for monitoring production from subsurface geological reservoirs, comprising:

- conducting a number (n) of three or more seismic surveys of said subsurface geological reservoirs with large time differences such as months or years, producing a number (n) of three or more vintages of seismic data di(t), d₂(t),.., d_n(t)= d (t) , each said vintage of seismic data di(t) comprising one or more seismic trace, - formulating an operator between up to n(n-l)/2 possible pairs of said seismic data sets d (t) and their corresponding up to n(n-l)/2 time shifts At (t) ;

- solving for (n-1) time shifts At (t) between (n-1) selected pairs of vintages of said seismic data di(t), d₂(t) , .., d_n(t) as an

overdetermined system of equations,

- applying one or more of said calculated (n-1) time shifts At (t) as input data for controlling further production from said subsurface reservoirs .

26. The method of claim 25, - while selecting the time shifts to be solved, selecting the (n-1) time consecutive shifts Ati-₂ (t) , At₂- ₃(t) , .., At_n-i,n(t) between consecutive vintages of said seismic data di (t) , d₂ (t) , .. , d_n(t) .

27. The method of claim 25, in the step of formulating said operator between the n(n-l)/2 possible pairs of said seismic data sets d (t) by

- calculating partial time derivatives from at least one base trace

- formulating said operator between said n(n-l)/2 possible pairs of seismic data sets d (t) and their corresponding n(n-l)/2 time shifts

At (t) as —a— d_b (t) A t{t) = Ad , simplified to L A t{t) = Ad wherein L comprises

a t

the partial time derivatives .

28. The method of claim 27, further in the step of formulating said operator,

- forming a set of n(n-l)/2 equations LAt = Ad which is overdetermined in said series of consecutive time shifts Ati_₂(t) , At₂_₃(t) , .., At_n- l,n(t) .

29. The method of claim 25, each said vintage of seismic data di(t) comprising seismic traces representing a three-dimensional volume of part of the subsurface geological formations, said resulting series of consecutive time shifts At (t) representing time shifts in said volume part of the subsurface geological formations .

30. The method of calculating seismic time shifts of claim 25, comprising calculating partial time derivatives — d_b (t) of the traces d t

by estimating coefficients of basis functions.

31. The method of calculating seismic time shifts of claim 30, said basis functions being spline functions or Laplacian functions.

32. The method of calculating seismic time shifts of claim 31, the order of spline function coefficients being low such as one to four.

33. The method of calculating seismic time shifts of claim 31, the number of Laplacian basis function coefficients being between two and five .