Method for Obtaining Estimates of a Model Parameter so as to Characterise the Evolution of a Subsurface Volume The present invention relates generally to the field of geosciences and more particularly to seismic data processing. Specifically the invention relates to a method for obtaining estimates of a model parameter (such as time shift, velocity change or time strain) so as to characterise the evolution of a subsurface volume over a time period.

In the oil and gas industry, seismic surveys are carried out in order to provide subsurface images so that accumulations of hydrocarbons or other fluids might be identified. In a seismic survey, one or several sources emit elastic waves in the form of pressure or ground motion modulation from specific locations (wavefield), at or below the land or sea surface or in a borehole. This wavefield propagates away from the source(s) through the subsurface. Along with this propagation, a fraction of the incident wavefield is reflected from the heterogeneities in the elastic material properties of the subsurface (such as acoustic impedance). This excitation by the incident wavefield generates a reflected wavefield from the heterogeneities, which manifests as pressure, particle motion or some derived quantities and can be detected and recorded at the surface or in a borehole at a number of receiver locations.

Processing of the measurements is undertaken so as to construct a 3D image of the subsurface. Repeated surveys at selected time intervals (days, months, years) allow observation of the changes in, over or under a given reservoir across the time interval - e.g. before oil or gas production starts and after some period of production or injection and to compare the results of measurements. This is called 4D seismic and involves comparing 2D or 3D seismic surveys carried out at different time instances. The aim is to observe changes in the state of the formations and fluids consequent upon production of hydrocarbons from or the injection of fluids into a reservoir. Proper detection of the changes and proper identification of the effects, factors and processes requires specialised acquisition techniques and data processing steps.

The data within the seismic data sets may be first processed to compensate for variations in acquisition (or non-repeatability of seismic surveys) and changes in velocity in the sub-surface.

In EP 1 865 340 to the Applicant, and incorporated herein by reference, the evolution of an oil reservoir in the process of producing is carried out by jointly inverting for the changes in the propagation times and seismic amplitudes of a seismic wavelet along propagation paths in the ground. Inverting allows to back filter, in effect, deriving the original from the solution. A base survey of the reservoir is provided, with a set of seismic traces at a first time T associated to a first velocity field Vb ; a monitor survey of the reservoir is provided, the monitor survey being taken at a second time T + ΔΤ, with a set of seismic traces associated to the same positions as in the base survey; the monitor survey is associated to a second velocity field I For a set of samples ; ^' in the base survey, one computes over the samples of the set the sum S of a norm of the difference between:

- the amplitude bi of the seismic trace in the base survey at each sample ; ^' and - the sum of the amplitude mr of the seismic trace at a time-corresponding ; ' in the monitor survey and the amplitude due to the reflectivity change local to the time- corresponding sample ; ' induced by the difference between the first velocity field Vb and the second velocity field V _m ; the time-corresponding sample /' being shifted in time by a time-shift derived from the velocity changes along the propagation path from the surface to time-corresponding sample ;'. This sum is minimised to derive the velocity changes from the base survey to the monitor survey and thus characterise the evolution of the reservoir.

This analysis is based on the fact that changes in the reservoir, due to exploitation, will cause changes to the petrophysical properties of the rock and therefore to the seismic velocity field. Practically, oil will be substituted by gas or water and/or the fluid pressure will change, modifying saturation, porosity, permeability and pressure, and consequently in velocity. Changes within the reservoir may also perturb the stress and strain state of the surrounding rocks, further altering their velocities. These changes to velocity will produce time shifts in the seismic response of underlying reflectors and associated changes in reflectivity, causing an alteration of the local wavefield. By using an inversion technique, for every point in the 3D volume, an estimate of the 4D changes having occurred in the time lapse between collection of the base and monitor surveys is provided. It is therefore possible to deduce a field of 4D velocity changes without having to proceed with cross correlation of the traces.

Although the 4D inversion problem appears relatively easy to formulate as the minimisation of a difference between base and monitor seismic data, it is an ill- posed problem that has multiple solutions: for instance, any smooth zero-mean velocity changes map into zero time-shift and does not generate any 4D amplitude difference. Moreover the inversion becomes even more highly non-linear for fields that induce subsidence and have potentially large time shift.

In EP 1 865 340, the crucial step is in minimising the difference between base and monitor seismics. Essentially this is an optimisation problem which requires minimising of the objective function or cost function over all choices of variables i.e. velocity changes that satisfy the modelled constraints. Usually the cost function is computed over all the available time-samples but it can be also calculated for decimated time samples or the sample number can be increased by interpolation to improve the accuracy of the solution.

However, as in almost any inverse problem, this cost function does not go identically to zero. In fact the forward model used for this inversion, is just an approximation which implies some assumptions, and therefore a residual still exists. To partially overcome this problem a constraint is added to the cost function. However, it is difficult to choose an optimal weight for this constraint. Several strategies can be used and a lot of interpretation is required, meaning that the solution is not unique. The problem of estimating time shifts in 4D seismic data in the overburden (where a signal originates from the geomechanical effects and no amplitude change exists) has been studied for over 10 years. Non-linear inversion (warping) techniques such as described above, comprise one approach that has been used to address this problem. Such techniques can struggle to accurately estimate large shifts due to cycle skipping. Other techniques include local cross-correlation and image processing warping, but these techniques also have their drawbacks.

It is desirable to provide a non-linear inversion technique for estimating time shifts which address at least one of the drawbacks of the prior art.

SUMMARY OF INVENTION

In a first aspect of the invention there is provided a method for characterising the evolution of a subsurface volume, said method comprising:

providing a first survey of the reservoir with a first set of seismic traces corresponding to a first time;

providing a second survey of the reservoir with a second set of seismic traces corresponding to a second time;

constructing an analytic filter operable to shift one or more seismic traces in dependence of a model parameter;

performing an inversion to obtain estimates of said model parameter such that said analytic filter aligns said first survey and said second survey. Other aspects of the invention comprise a computer program comprising computer readable instructions which, when run on suitable computer apparatus, cause the computer apparatus to perform the method of the first aspect; and an apparatus specifically adapted to carry out all the steps of any of the method of the first aspect. Other non-essential features of the invention are as claimed in the appended dependent claims. BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the invention will now be described, by way of example only, by reference to the accompanying drawings, in which:

Figures la and lb are schematic illustrations of respectively, a base survey and a monitor survey being performed;

Figures 2a and 2b are seismic trace diagrams illustrating the principle of plane wave destruction filters; and

Figure 3 is a flowchart describing a method according to an embodiment of the invention.

DETAILED DESCRIPTION OF THE EMBODIMENTS

Referring initially to Figures la and lb there is illustrated a reservoir, generally indicated by reference numeral 10, containing hydrocarbons 12 in the sub-surface 14. A survey vessel 16, upon which is located a sonar transmitter 18, being an acoustic source, and an array of receivers 20, performs a survey by travelling over the reservoir 10. The first or initial survey, Figure la, may be referred to as a base survey and is typically performed in the exploration phase before production begins.

The base survey of the reservoir 10 provides a set of seismic traces at a first time T. For a given trace, the base survey provides amplitudes that are a function of time. With digital recording and processing the trace is sampled at a set of values; typical trace lengths correspond to around 1000 samples. The trace is then handled as a set of values.

One or more wells 22 may be drilled in order to extract the hydrocarbons 12. As the reservoir 10 is produced, hydrocarbons will be substituted by other fluids and the fluid pressure will change. Additionally, enhanced oil recovery techniques may be applied wherein a fluid is injected into the reservoir at one or more locations giving changes in fluid pressure and saturation. Changes within the reservoir may also change the stress and strain state of the surrounding rocks. Thus when a further survey is carried out, Figure lb, these changes will be seen due to a consequential change in the velocity field. These changes to velocity will produce time shifts in the seismic expression of underlying reflectors and associated changes in reflectivity, causing a change in the local wavefield.

Thus reservoir monitoring performs a monitor survey of the reservoir 10, taken at a second time T + ΔΤ, with a set of seismic traces. In the simplest assumption, ΔΤ is a positive quantity, and the monitor survey is taken at a time later than the base survey; however, the order in which the surveys are taken is irrelevant to the operation of the process of the invention and, in principle, the time lapse ΔΤ could as well be negative - which amounts to comparing the earlier survey to the later one. As for the base survey, a sampled trace in the monitor survey is represented as a set of values. Ideally, the traces in the monitor survey are associated to the same positions as in the base survey. This is carried out by using, inasmuch as possible, the same equipment, acquisition geometry and processes for running the base and monitor surveys. Techniques such as interpolation may be used where traces in the monitor survey and in the base survey do not fulfil this condition.

Non-linear inversion, often referred to as warping, is a known technique for estimating time shifts from 4D base and monitor survey data in the overburden. An iterative solver is used to find a minimum (with respect to time shifts) of an objective function of the form E = \\b— /(At) ll _p - The vector b contains the base traces, the vector At is the time shift estimate and the non-linear function f(At) applies At to the monitor traces. To solve for At involves linearizing f(At) around a At estimate and solving for an incremental At update. Linearization of the function f(At) requires calculating its time derivative at each iteration. It is proposed herein to utilise the theory of plane wave destruction filters, as described by Sergey Fomel, in Geophysics, 67, 1946-1960 (2002)- herein incorporated by reference) to minimise the difference between base and monitor surveys. The concepts described herein use the plane wave filter construct described by Fomel, but alters the filter parameterization from local dip to 4D time shift. Doing so generates a new set of filters suitable for aligning 4D base and monitor traces when parameterized with the correct 4D time shifts. This leads to a non-linear inverse problem for 4D time shifts which has three advantages over existing non-linear inversions for 4D time shifts: 1) Within the inverse problem the 4D time shifts are now a function of the filters, not the monitor trace. For each non-linear iteration it is necessary to only rebuild and apply the filters, which is a quick approximation to warping the monitor.

2) The filters are analytic, as are their gradients, making the Jacobian fast to calculate. Other inversions for 4D time shift must warp and numerically differentiate the monitor traces to rebuild the Jacobian at each non-linear iteration.

3) Due to the efficiency of this method, spatial regularization can be applied which reduces the 4D noise while maintain the 4D signal, without significant additional computing cost. Equivalent spatial regularization is not normally performed in 4D warping as it would slow the computation down due to the comparatively large number of iterations required to obtain a solution. Instead only regulation in the z-direction is normally performed.

In a preferred, but optional, embodiment, rather than inverting for the time shifts which align the monitor to the base, it is proposed to invert for the time shifts which move the base and monitor an equal amount in opposite directions so they align at a mid-point. This brings two further benefits:

4) The time shifts required to align the base and monitor are halved, which makes interpolation errors smaller, increases the range of time shifts for which the inversion is accurate, and reduces the number of non-linear iterations required to converge to a solution.

5) Interpolation errors are equal in both base and monitor traces, so have no contribution to the cost function.

Plane Wave Destruction

From the plane wave description of a trace it is possible to use shifts σ, to parameterise a filter L(o), that warps one trace to match another trace. This is illustrated in Figures 2a and 2b. Figure 2a shows a first trace ti and a second trace t∑. Also shown is the difference between the two traces t2-ti (the residual) without application of a filter. Figure 2b again shows first trace ti and a second trace t∑. Next to this is the filtered second trace L(o)t2. This filtered second trace L(o)t2 \s such that trace ti - L(o)t2 is zero, as shown by the rightmost trace in this Figure (a straight line). In other words, the filter L(o) destroys the energy difference between plane waves ti and ¾ and is consequently referred to as a "plane wave destruction (PWD) filter". The system: ti - L(o)t2 = 0, can be solved as a non-linear inversion for shifts, σ.

A PWD filter, L(o), maps a trace, ¾ to a trace, ti, through shifts, σ:

ί _! - (σ)ί ₂ = 0 W

Representing traces as plane waves, the frequency domain ti trace, Ti, can b written as a function of its amplitude, Ai, and phase, e ^ie :

For a constant time shift the frequency domain trace, T2, can be written as:

J ₌ 4 gi{kx-o)(t+o)) _ gi(kx-<at) g-ίωσ [3]

If it is assumed that there is no energy change between Ti and T2, Ai equals A2 and:

A constant shift of a plane wave is a frequency domain multiplication with phas shift operator, ε ^ίωσ .

For time varying shifts, one option would be to perform short window Fourier transforms and invert for the optimum σ per window. Instead it is proposed to find the time domain convolution filter analogue to the phase shift operator ε ^ίωσ .

The condition of no energy change between Ti and T2 is met by an all pass (phase only) digital filter, which in Z-transform notation is:

T ₂ e ^laa « T Z _t ) ^B{ ^ _r ^{(5) 2 2} B(z ¹ )

For a 3-term filter: b^Z •; -

5 ₃ (Z _t ) 1 ¹ + b ₀ + b _x Z _t b_ _x e ^m + b ₀ +

(6)

5 ₃ (Z ¹ ) b_ _x Z _t + b ₀ + b _x Z b_ _x e ^ia + b ₀ + b _x e

The filter coefficients, b _n , can be found by fitting its response at low frequencies to that of the phase shift operator ε ^ίωσ . Using a Taylor expansion around zero frequency, the phase shift operator and the all pass filter can be approximated as:

2 2 3 3

e »1 + ιωσ 1 + ... (7)

2 6

B(Z _t ) _|

Equating like-terms provides a system of non-linear equations relating σ to the filter coefficients b _n :

σ_ - ) ₌ jb -b _σ 3 _2(^ -hjbl -4b ₀ {b +fc,X¾ -4ft A + ¾ ₁ ² )j (8) b^+b.+b, ^~ {b _i+ b ₀ +b ² 3(^+¾+^) ³

As this system does not uniquely define the coefficients, b _n , an additional constraint is required. An additional constraint ensures the filter E>3(Z _t ) does not alter the zero frequency component:

_Bj (\) = b_ _l +b ₀ +b _l =\ (9)

This also allows the system to be solved uniquely: , ₌ (ΐ-2σ)(2-2σ) _{b =} (2 + 2σ)(2-2σ) _& _ (l + 2σ)(2 + 2σ) ₍₁₀ )

^_1 12 ⁰ 6 ¹ 12

These are the coefficients for the centred 3-term PWD filter, which is the time domain analogue of the frequency domain phase shift operator: Te^ * ±b _H {a)t,_ _H =L{a)t _fnl

=-1 ^ where s is the time sample index (running from 1 to the number of samples) of trace t. The time domain filter can be parameterised with a different σ at every time sample.

The PWD filter derivation of equations (5) through (10) follows precisely the method of Fomel (2002) wherein σ equals the local dip from one trace to its neighbour. Local dip can be defined as the ratio of time shift to horizontal shift required to align neighbouring traces. Rather than predicting one trace from its neighbour, the 4D time shift estimation problem can be thought of as predicting a base trace from a collocated monitor trace by shifting it along local time shifts. The filter that performs such shifting is analogous to that defined by equation (10), with the local slope term, σ, replaced by the local time shift, At:

The coefficient t is the time sampling rate of the data and allows the local time shifts, At, to be expressed in real time units. Note that the filter described by coefficients in Equation (12) is that required to shift the monitor to the base, or indeed the base to the monitor. The preferred embodiment is to invert for the time shifts which move the base and monitor an equal amount in opposite directions so they align at a midpoint, thereby exploiting the aforementioned advantages of doing so. In this case the filter coefficients become:

(13)

Knowing the coefficients b _n (equation 13) allows us to construct a non-linear inversion for At:

-At)t _b - L(At)t _m = [L(-At) - L(At)] C(At)a = 0 (14)

The inverse problem can be solved by iteratively updating an initial time shift estimate, Ato, using the linear system for incremental time shift, 6t:

diag[C(At ₀ )a]5t + C(At ₀ )a = 0 C ¹⁵⁾ For one pair of traces, the matrix C(At) is:

The PWD filter, L _s (At), at each sample index, s, is:

CYAtJ is the differentiation of C(At) with respect to At:

Where L _s '(At) is the differentiation of L _s (At) with respect to At

Substituting G = diag[C'(At)a], d = -C(At)a, gives the formal least squares solution of equation 15:

5t = (G ^T G) ^l G ^T d ™

However, the inversion is unstable. G ^T d and G ^T G are diagonal and there is noise in the system. To stabilise the inversion a regularizing operator, R can be used to force correlation between model parameters:

(21)

The regularizing operator, R, may comprise a Tikhonov regularization operator which acts on all elements in St. Its impact is controlled by the parameter λ. This can be used to force 6t to have some property, such as low energy or smoothness.

It is possible to recast the inversion to solve for time strain, u, by recognising that the time shift, At(x), at time x = T is:

At(T) = t∑u(x) (22) Substituting this into the inversion yields:

(23)

Where S _z maps 5u to St by performing the running sum of 5u in the z-direction before multiplying by the time sampling rate, t.

Practical application Figure 3 is a flowchart describing a practical embodiment of the techniques described above. The method uses a Gauss-Newton scheme as follows:

Step 200: Make an initial guess of the time shifts (usually zero). Step 210: Construct the time shift filters and their derivatives from the current time shift model, At, (where i | i £ M, 0 < i < N and indicates to the current nonlinear iteration number). The filter coefficients are those of Equation (13). While these are described in relation to single traces i¾, t _m , the objective function is cast in terms of vectors b and m, which comprise all base and monitor traces respectively. The unregularized objective function to be minimised may take the form:

0 = \\L(-M)b - L(Ai)m\\ ₂ (24) where L is a non-stationary, phase-only filter parameterised by vector At representing all the time shifts. It should be noted that Equation (24) shifts both sets of base and monitor traces by equal amounts in opposite directions, bringing them to common midpoints in time. This is optional, and shifting the monitor to the base, or indeed the base to the monitor, is envisaged within the scope of this disclosure. As already described, the filter has an analytic expression derived from the plane wave model of a seismic trace, and has an analytic gradient. Here, the three-term filter coefficients are described by Equation (13) and the coefficients of the derivative of the three term filter are described by Equation (19). Step 220: Generate the data residuals for the current iteration, rf,, by applying the time shift filters to the base and monitor data in the following way:

b

[L(-At _t ) -L(M _i )] r \ = d _i (25)

m

Step 230: Build the Jacobian matrix, Gi, by applying the filter derivatives to the base and monitor data and diagonalizing the result.

Step 240: Solve the linearized inverse problem for an incremental time shift update, 6tr.

GiStt = -di (27)

Iterating through this process provides local time shifts which align the base and monitor data. However, the time shift at a particular subsurface location is an accumulation of all time shifts down to that point, and not just a local anomaly. To understand the physical cause of an unobserved time shift, it should preferably be an integrated measurement. Changing the inversion model parameters from time shifts, At, to time strains, Au, provides a more physically understandable output which is better suited to geomechanical modelling. This can be done by substituting the relationship between time shift and time strain expressed by Equation (22) into Equation (27). Additionally, one or more regularisation terms may be added to stabilise the inversion in the presence of noise.

Recasting Equation (27) in terms of time strain and adding Tikhonov regularization yields:

Where the operator V produces the gradient in the direction specified by the subscripts x, y and z. Naturally, the inverse of this operator performs integration. The weighting terms λ control the severity of the model damping in each of the subscripted directions. Iterating using Equation (28) provides spatially smooth time strains and vertically smooth accumulated time shifts which align the base and monitor data. Equations (23) and (28) are equivalent when the regularization operator R is chosen to be the 3D gradient, V _xyz . The result is well suited for building geomechanical models. Using this formulation, all variants of the existing warping methods can be implemented, e.g. multi-monitor warping as disclosed in WO2010/139998, or layered based approaches such as Propagation 4D as disclosed in WO2012/038647. Both of these documents are hereby incorporated by reference.

Step 250: Differentiate the time strain estimate to give an incremental time shift update (where the inversion model parameters were changed to time strains) for adding to the current time shift model. At each iteration, the inverted 5t, is added to the current time shift estimate Δ£,, to give an updated set of time shifts Δ£,+ί.

Step 260: Repeat steps 200 to 250 until convergence.

In an optional embodiment, it is proposed to stabilise the inversion described herein using a novel regularization method. In this regularization method, it is proposed to apply constraints per region of subsurface volume, the regions defined according to a high order discrete property (e.g. geomechanical facies) attributed thereto in accordance with the region's geologic characteristics. Regions attributed with the same high order discrete property can be assumed to have similar geomechanical properties.

One way of applying constraints per region defined according to a high order discrete property is to define velocity changes (or time stain or time shifts depending on the chosen parameterization) as a regional function defined by a limited number of parameters at each position per region (e.g. layer or group of layers) used in geomechanical modelling. Geological surfaces can be used to divide the model space into layers (or strata) of geomechanical facies (rocks with similar geomechanical properties). Coefficients of a regional function are then inverted for at each spatial location within a layer. The regional function describes how an elastic property, such as time strain, varies within each region (e.g. layer). The inversion forward modelling operator includes the regional function, which expands the coefficients at each spatial location into the desired property, such as time strain, at every seismic sample point within the region/layer. Determination of layer positions (through geological interpretation) does not need to be very accurate for measurements within the overburden/underburden as the wavelength of the geomechanical property variation are very long.

A simple example of a regional function is constant time strain per region or layer, in which case only a single coefficient is required at each spatial location within a region/layer. However other regional functions may also be appropriate depending on the subsurface geology. Therefore the constraint enables inversion within each region or layer using any function. The example below explains this method of geological constraints applied to an inversion for time strain from noisy 4D seismic, but the regularization method is general to any 4D seismic inversion (for geomechanical or elastic properties). Introducing geological layer constraints into equation (28), without the Tikhonov regularization, gives:

Gi V^HS i = di (29)

The new model parameters, ίϊ, are layer-dependent coefficients of a function, H, which is evaluated at every sample point. Geological surfaces are used to define the layers of geomechanical facies; rocks that can be assumed to have similar geomechanical properties. The simplest expression of this constraint is that time strain is constant within each layer, but any layer-dependent rock mechanics relationships could be expressed by Hu. For the case of constant time strain per layer, ίϊ is a single time strain value at each spatial location per layer. The mapping operator, H, for each trace location is built of N matrix blocks, each with N columns, where N is the number of layers. The number of rows per block varies according to the number of z-samples in the corresponding layer. Each block has a single column of ones, which is found in the first column for the first block, the second column for the second block, and so on. All other elements are zero. The final matrix, H, is a vertical concatenation of the blocks. Under these circumstances the action of V ^{~ 1} H can be expressed by:

St(T) = SuQ) Su(j)x(j) j \ j £ n ≤ j (30)

Where 6t(T) is the incremental time shift at two way time T, Su(j) is the time strain in layer j, x(j) is the time thickness of layer j, and the index j runs from one up to the layer / in which time T exists. Equation (29) can be further regularized using classical Tikhonov methods to give:

It is important to note that the lateral smoothing action of _X V _X and _y V _y in equation (31) takes place within geological layers by construction, thereby following layer dip. The vertical smoothing governed by A _; -V _j occurs across layer boundaries and guards against the estimated model parameters swinging (successive over- and under-estimation) from one layer to the next. This constraint is not applied for reservoir layers.

One or more steps of the methods and concepts described herein may be embodied in the form of computer readable instructions for running on suitable computer apparatus, or in the form of a computer system comprising at least a storage means for storing program instructions embodying the concepts described herein and a processing unit for performing the instructions. As is conventional, the storage means may comprise a computer memory (of any sort), and/or disk drive, optical drive or similar. Such a computer system may also comprise a display unit and one or more input/output devices.

The concepts described herein find utility in all aspects of surveillance, monitoring, optimisation and prediction of a subsurface volume comprising a hydrocarbon reservoir and well systems, and may aid in, and form part of, methods for extracting hydrocarbons from such hydrocarbon reservoir and well systems.

It should be appreciated that the above description is for illustration only and other embodiments and variations may be envisaged without departing from the spirit and scope of the invention.