Method of constraining an inversion in the characterisation of 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 of constraining an inversion in the characterisation of the evolution of a subsurface volume. 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, the overburden and underburden.

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.

Even greater difficulty is encountered when estimating time shifts in 4D seismic data in the overburden or underburden (that is in the regions above or below the reservoir). In these regions, a signal originates from the geomechanical effects and very weak amplitude change exists. These signals are much weaker, over much larger depth than signals from within the reservoir. Non-linear inversion (warping) techniques such as described above, can be used to address this problem, but struggle in handling this type of signal. Other techniques include local cross- correlation and image processing warping, but these techniques also have their drawbacks.

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. Moreover, noise in the data is quite high and sometimes of the same order of magnitude as the 4D signal. To partially overcome this problem a regularisation is applied to the inversion. However, it is difficult to choose an optimal regularization without having a detrimental impact on the magnitude and the localisation of the result in geomechanical modelling. The regularization may comprise a combination of model energy damping, model smoothing and/or model sparsifi cation, and is controlled by unintuitive weighting parameters. Several strategies can be used and a lot of interpretation is required, meaning that the solution is not unique. Particularly, in the case of overburden or underburden analysis, the nature of the signal, the ambient noise, and the very large number of unknowns in the inversion, means that the inversion is unstable without a very high level of regularisation.

It would therefore be desirable to improve regularisation when performing a 4D inversion.

SUMMARY OF INVENTION

In a first aspect of the invention there is provided a method for characterising the evolution of a subsurface volume over time, 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;

performing an inversion to obtain estimates of changes having occurred between said first time and said second time in terms of at least one model parameter; wherein said subsurface volume is divided into regions according to a high order discrete classification attributed in accordance with the region's geologic characteristics, and said inversion is constrained per region of subsurface volume. 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;

Figure 2a shows two traces from a classical inversion; and

Figure 2b shows two traces from an inversion constrained according to a method 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 (especially pressure changes) will change the stress state inside the reservoir and 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.

4D time shifts in the over-burden and sub-reservoir can be caused by production related geomechanical changes. Estimating such time shifts from 4D seismic data is important for calibrating geomechanical models. This calibration can be achieved by deriving the fractional change in travel-time ('time strain') from 4D time shifts, and assuming a relationship between time strain, velocity change and true vertical strain. The estimation of time strain is crucial, but complicated by relatively weak geomechanical signals in the over-burden and sub-reservoir (compared to the 4D signal in a reservoir due to changes of pressure and saturation). It is therefore very difficult to correctly estimate its magnitude and locality.

Time shifts induced by geomechanical effects accumulate through the overburden, making the time shift at a particular point equal to the sum of the travel-time difference due to velocity and depth changes at all points above. Methods reliant on local cross-correlation are either overly sensitive to noise (if the correlation window is small) or provide poorly localized results (if the correlation window is large) and are therefore ill-suited for calibrating geomechanical models. Dynamic time warping methods estimate local time shifts sample-by-sample and are extremely sensitive to noise. They too are poorly suited for calibrating geomechanical models.

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. The inversion should be regularized in order to converge to a stable solution.

Non-linear inversion for time shifts can be directly parameterized in terms of relative interval velocity change or time strain, thereby implicitly providing cumulative time shift solutions. However, solving for time strain is less stable and more noise prone than solving for local time shifts. Significant inversion regularization is required which can have a detrimental impact on the magnitude and the localisation of the result in geomechanical modelling. This regularization is typically a combination of model energy damping, model smoothing and/or model sparsification, and is controlled by unintuitive weighting parameters.

In order to regularize inversion of 4D differences induced by geomechanical changes, 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 a (e.g. elastic) model parameter, 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).

A non-linear inversion for time shifts can be constructed using an objective function of the form E = \\b— /(At) || _p . The vector b contains the base traces, the vector At is the time shift estimate and the non-linear function f(M) applies At to the monitor traces. The system of equations solved at every non-linear iteration can be written as:

GCvi _n = /(At _n ) - b (--Q The matrix G is the Jacobian of the function f(M _n ) at non-linear iteration n. The matrix C performs integration (scaled by the time sampling) of the time strain μ„ (to convert to time shifts). It is because μ„ is defined at every sample point that this inversion is unstable. Introducing geological layer constraints into equation (1) gives: GCLy^ = f(M _n ) - b (2)

The new model parameters μ _η , are layer-dependent or region-dependent coefficients of a function L, which is evaluated at every sample point. Geological surfaces are used to define the regions/layers of geomechanical facies; rocks that can be assumed to have similar geomechanical properties. For this specific example the constraint is that time strain is constant within each layer, but any region- or layer- dependent function could be expressed by Ιμ _η .

For the case of constant time strain per layer, model parameters μ _η comprise a single time strain value at each spatial location per layer. The mapping operator L, 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, L, is a vertical concatenation of the blocks.

A significant advantage of using geological constraints when inverting 4D seismic is the ability to stabilize the inversion without using physically unrealistic regularization methods with unintuitive weighting parameters. Using geological surfaces to segregate the model space into layers of geomechanical facies allows inversion for a small number of coefficients per layer at each spatial location which define the model parameter variation. Layer-dependent rock mechanics relationships are incorporated directly into the inversion. This is in contrast to existing methods which invert for properties at every sample point, and stabilise the process by smoothing, sparsifying or limiting the energy of the model space; all of which create a trade-off between an unstable or over-regularized result. Unlike existing regularization methods, geological constraints are directly related to prior knowledge of subsurface geomechanics. Using geological surfaces to constrain the inversion also helps to impose discontinuities where they are expected; at the boundaries of geomechanical facies. Using geological constraints does not preclude the use of the traditional regularization methods. Inverting noisy 4D seismic may require additional model energy damping or smoothing. Combining such regularization with geological constraints ensures its effect is consistent with geological layering. This is particularly advantageous when regularizing towards a spatially smooth solution; without geological layering smoothness is enforced horizontally (often conflicting with structural dip), while geological layering ensures the smoothness constraint follows layer dip by construction. Similarly, such an inversion will naturally take into account any fault, dip or similar anomaly. Therefore the throw of two stratigraphically neighbouring points on either side of a fault will be honoured.

A practical advantage of using geological constraints is the decreased computational cost of solving the inverse problem due to the vastly reduced model space. Without using layers to constrain the inversion, the inverted model size is the product of the number of inlines, crosslines and z-samples. Using layers the model size is the product of the number of inlines, crosslines, layers and model coefficients per layer. Noting that there are many z-samples per layer, this represents a very large saving. This saving is illustrated in Figures 2a and 2b. Figure 2a shows two traces following inversion by a classical method. The inversion needs to be performed for each sample point, of which there can be thousands per trace. While the layers are shown in Figure 2a, they are not used in classical inversion. Figure 2b shows the inverted traces with the inversion constrained according to the teaching disclosed herein. In this example, the constraint is that the parameter is constant per layer. Consequently there is only one parameter to invert for per layer, per trace. It should be appreciated that the result of the inversion for each layer is not simply an average of the values within a layer following classical inversion, as the classical inversion results will have been skewed by the regularization necessary to stabilise the inversion. A number of specific examples of the basic concept of regularising the inversion according to layer or facies region will now be provided. The first of these is related to Propagation 4D techniques as disclosed in WO2012/038647, which is hereby incorporated by reference. This document describes methods for the estimation of (e.g. elastic) model parameters (such as density, pressure and/or shear wave propagation velocity) in pervious layers located in a dense underground area along an array of horizontal positions, by inverting 4D seismic data. This inversion is layer based, whereby one or several seeds are provided (generally using Cal4D as described in WO2012/010790 herein incorporated by reference) with a set of layers and an estimation of the changes of model parameters in each layer. By using the seeds as an initial solution, the Propagation 4D method inverts for the position of the layers and the values of the changes of each model parameters at the neighbouring traces around the seeds. Propagation 4D is typically performed to measure parameters within a reservoir, rather than the overburden. It typically uses an inversion to find layer positions. It uses knowledge from a previous trace to invert for changes in the next trace (spatially). Therefore, beginning with known positions of layers of a first trace, an inversion is performed to find the layer positions of the next trace (which should not have changed much). Then using this next trace, the layer positions of the next trace is inverted for, and so on.

The Propagation 4D method for estimating model parameters of a region of the subsoil meshed according to a network of horizontal positions comprises the following steps, according to an embodiment:

(a) measuring, in a first phase, base seismic traces associated with said horizontal positions;

(b) measuring, in a second phase, monitor seismic traces associated with said horizontal positions and corresponding to the base seismic traces; (c) obtaining an estimate of the variations of model parameters between the first and the second phases in permeable layers of the subsurface in at least one starting position of the network, the model parameters including a density and a speed of propagation of a pressure wave in each of said permeable layers;

(d) determining a boundary consisting of the positions of the network which are adjacent to a starting position;

(e) estimating, at each position of the boundary, the variations of the model parameters in the permeable layers and/or the depths and/or thicknesses of said permeable layers by minimizing a cost function derived from a seismic wave propagation model accounting for the evolution between the base and monitor seismic traces associated with said position of the boundary, the minimization taking into account the variations estimated at the adjacent starting position, and marking each position of the boundary as having been processed;

(f) removing from the boundary any position where the minimized cost function exceeds a limit value;

(g) selecting at least one propagation position in the boundary, where the minimized cost function has, within the boundary, a minimal value, then removing the selected propagation position from the boundary;

(h) adding to the boundary each position of the network not marked as processed and adjacent to a propagation position selected in the preceding step (g);

(i) estimating, at each position added to the boundary in the preceding step (h), the variations of the model parameters in the permeable layers and/or the depths and/or thicknesses of said permeable layers by minimizing the cost function derived from the seismic wave propagation model accounting for the evolution between the base and monitor seismic traces associated with said position added to the boundary, the minimization taking into account the variations estimated at the selected propagation position, and marking each position added to the boundary in the preceding step (h) as having been processed;

fj) removing from the boundary any position added in the preceding step (h) where the minimized cost function exceeds the limit value; and

(k) repeating from step (g) if the boundary still contains at least one position adjacent to a position of the network not marked as processed.

The Propagation 4D algorithm described above makes it possible to evaluate, step by step, the model parameters and/or the geometry of the layers. The propagation is performed in such a way as to retain the best possible results for the inversion with limited perturbations of the model. The propagation stops when the assumption of relative continuity of the values of the parameters is no longer in agreement with the measurements.

The Propagation4D and Cal4D algorithms, as disclosed in the aforementioned documents, require strong changes of (e.g. elastic) model parameters between layers, which correspond to high frequency amplitude variation of the seismic signal at the layer boundaries. These methods are therefore presently not suited for estimating time strains in the overburden, which are weak and essentially continuous.

It is proposed to use geological layer or region constraints, consistent with layers or subsurface regions as used in geomechanical modelling, to define and fix the layer positions within the Propagation4D method. As already mentioned, geological surfaces divide the model space into layers of geomechanical facies; rocks that can be assumed to have similar geomechanical properties. Coefficients of a regional function are inverted for at each spatial location within a layer. The inversion forward modelling operator expands the coefficients at each spatial location per layer (or region) into some desired property (such as time strain or relative velocity change) at every seismic sample point within the layer (or region), in accordance with the regional function. The example below explains this method of geological constraints applied to the Propagation4D method for inverting time strain from noisy 4D seismic, but the regularization method is general to any 4D seismic inversion (for geomechanical properties).

A non-linear inversion for time shifts can be constructed using an objective function of the form:

J ∑n(B _n - MB _n )J

Λ N

j (3)

∑n(Bn- M _n )i

Λ N where JBM is the coherency function between the base B and the monitor M computed over the number of samples N. MB is the monitor warped into the base by the non linear function f defined in equation (2) which, in this case, can be expressed as:

The vector M contains the monitor traces, the vector— is the velocity change

v _{V l} y &

(related to time strain) estimate in each layer 1. TI _Q ^ts represents a box function with values of 1 between [0, nt _s ] and 0 outside, with n being the sample number and t _s being the time sampling interval. Compared to classical Propagation 4D, the most important point to notice here is that the limits of the layers T ₁₍ T^ are no longer parameters of the equation, but are prior information corresponding to the geological constraint.

It should be noted that equation (4) shows the inversion being constrained to a constant value function for variation of velocity per layer. However, as before, it is possible to introduce any function other than this.

In some cases it is possible to mix this approach (where layer positions are given) with some uncertainty on the position of certain layers (e.g. in the reservoir), such that the amplitude term drives the accurate estimation of the layer positioning.

Another method disclosed in a patent application having agent's reference P166099.GB.01, filed the same day as this application (to the same applicant) discloses a fast non-linear inversion method for estimating 4D time shifts and time strain. This method is fundamentally different in construction to those described above, as the inversion is cast as a phase-only plane wave filter design problem. These filters are described analytically, making inversion for their coefficients fast and efficient. The method discloses:

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 time varying parameter (such as time shift, time strain or velocity change); and

performing an inversion to obtain estimates of the time varying parameter such that the analytic filter aligns the first survey and the second survey.

The analytic filter may comprise a plane wave destruction filter, and may be a time domain analogue of a frequency domain phase shift operator operable to impose a constant time shift to a plane wave model of a seismic trace. The analytic filter may be such that the time varying parameter is a function of the analytic filter instead of being a function of one or more of the second set of seismic traces.

The method may comprise the steps of:

constructing the analytical filter and their derivatives from a current time varying parameter model; and

iteratively updating an initial time varying parameter estimate with an incremental time varying parameter update until convergence on a solution by:

applying the analytical filter to at least one of the first survey and the second survey to generate residuals;

building a Jacobian matrix by applying the filter derivatives to the first survey and the second survey and diagonalizing the result; and

iteratively solving a linear inverse function constructed from the residuals and the Jacobian matrix.

The implementation of time strain inversion using plane wave filters described in the aforementioned document estimates a single time strain value at each seismic sample point, with a weak 3D correlation enforced between sample points. The stability of this method relies on signal variation between samples being significantly larger than the noise. This is often not the case in the overburden where time strains are weak and approximately continuous. Consequently it is proposed to constrain this inversion using the techniques described herein, that is to use geological layer or region constraints, consistent with layers or subsurface regions as used in geomechanical modelling, to divide the model space into layers/regions of geomechanical facies.

The unregularized objective function to be minimised may take the form:

E = \\L(-M)b - L(Ai)m\\ ₂ (5) where L is a non-stationary, phase-only filter parameterised by vector At representing all the time shifts. It should be noted that Equation (5) 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.

Making the substitutions: C _p = Δί _ρ ) — L(At _p )], d _p =—C _p G _p — diag c _p [ ^ j jand u _p = V _z At _p equation (1) can be linearized to give:

G _p tV _z ¹ Su _p = d _p (6)

Equation (6) is solved iteratively for an incremental time strain perturbation, 6u _p , where p is the iteration number. The matrix C is constructed from the partial derivatives of the filter coefficients with respect to At. The operator, V, produces the gradient in the direction specified by the accompanying subscript. Naturally, the inverse of this operator performs integration. The scalar, t, is the time sampling interval of the 4D data.

Introducing geological layer constraints into equation (2) gives:

G _p tV _z ¹ HSu _p = d _p

(7) 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 (8)

Where 5t(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 (7) 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 (9) 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. Equation (9) can be used to estimate layer-dependent overburden and sub-reservoir geomechanical properties (such as time strain) which vary smoothly in 3D, while honouring the a priori layering scheme.

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.