METHOD FOR OBTAINING ESTIMATES OF A MODEL PARAMETER SO AS TO CHARACTERISE THE EVOLUTION OF A SUBSURFACE VOLUME OVER A TIME PERIOD

Abstract

Abstract Disclosed is a method for characterising the evolution of a subsurface volume over time. The method comprises providing first and second surveys of the subsurface volume. Each survey comprises seismic data acquired by transmitting seismic signals into the subsurface volume and subsequently detecting some or all of the seismic signals after reflection within the subsurface. The first seismic data of the first survey corresponds to a first time and the second seismic data of the second survey corresponds to a second time. At least some of the first seismic data and the second seismic data is obtained with a non-zero offset. An inversion is performed to obtain estimates of changes having occurred between the first time and the second time in terms of at least one model parameter; wherein for the inversion: the first seismic data and the second seismic data is not processed to be equivalent to zero-offset data prior to the inversion; and it is assumed that the path taken by each received seismic signal between its transmission and reception is the same for the first survey and the second survey

Claims

1. A method for characterising the evolution of a subsurface volume over time, the method comprising: providing a first survey of the subsurface volume with first seismic data acquired by transmitting seismic signals into the subsurface volume and subsequently detecting some or all of the seismic signals after reflection within the subsurface, the first seismic data corresponding to a first time, at least some of the first seismic data being obtained with a non-zero offset; providing a second survey of the subsurface volume with second seismic data acquired by transmitting seismic signals into the subsurface volume and subsequently receiving some or all of the seismic signals after reflection within the subsurface volume, the second seismic data corresponding to a second time, at least some of the second seismic data being obtained with a non-zero offset; performing an inversion to obtain estimates of changes having occurred between the first time and the second time in terms of at least one model parameter; wherein for the inversion: the first seismic data and the second seismic data is not processed to be equivalent to zero-offset data prior to the inversion; and it is assumed that the path taken by each received seismic signal between its transmission and reception is the same for the first survey and the second survey.

2. The method as claimed in claim 1 wherein the inversion is performed using an objective function which is not dependent on the distance travelled by the received seismic signal.

3. The method as claimed in claim 2 comprising the step of performing a ray tracing operation to determine ray paths which approximate the paths taken by each received seismic signal between its transmission and reception.

4. The method as claimed in claim 3 wherein the subsurface volume is discretized into cells and the objective function comprises a forward modelling operator comprising the travel times per cell for each of the ray paths.

5. The method as claimed in claim 4 wherein the ray tracing operation is performed only once and the forward modelling operator is used for all iterations of the inversion.

6. The method as claimed in claim 1 wherein the at least one model parameter comprises the propagation velocity of the seismic signal.

7. The method as claimed in claim 1 wherein the inversion comprises finding the propagation velocity changes which describe observed time shifts between the first seismic data and the second seismic data.

8. The method as claimed in claim 1 wherein the at least one model parameter comprises the travel time of the seismic signal and/or time strain.

9. The method as claimed in claim 1 wherein the first seismic data and the second seismic data are obtained at a plurality of different non-zero offsets.

10. The method as claimed in claim 1 wherein the at least one model parameter is an anisotropic model parameter such that the inversion estimates the changes having occurred between the first time and the second time in terms of the at least one anisotropic model parameter.

11. The method as claimed in claim 10 comprising the step of decomposing the anisotropic model parameter as a function of the angle of the path taken by each received seismic signal.

12. The method as claimed in claim 10 wherein it is assumed that the anisotropic model parameter is elliptically anisotropic.

13. The method as claimed in claim 10, wherein the first seismic data and the second seismic data is multi-azimuth seismic data, and the inversion estimates azimuthally varying changes having occurred between the first time and the second time in terms of the at least one anisotropic model parameter.

14. The method as claimed in claim 1 wherein the first seismic data comprises a first set of seismic traces and the second seismic data comprises a second set of seismic traces.

15. The method as claimed in claim 14 further comprising: constructing an analytic filter operable to shift at least one seismic traces in dependence of the changes in the model parameter; and performing the inversion by using the analytic filter to align the first survey and the second survey.

16. The method as claimed in claim 15 wherein the analytic filter comprises a plane wave destruction filter, and is 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.

17. The method as claimed in claim 15 wherein the analytic filter is such that the model parameter is a function of the analytic filter instead of being a function of at least one of the second set of seismic traces.

18. The method as claimed in claim 1 further comprising the step of using the results of the method to aid hydrocarbon recovery from a reservoir.

19. A computer program comprising computer readable instructions which, when run on suitable computer apparatus, cause the computer apparatus to perform the method of claim 1.

20. A computer program carrier comprising the computer program of claim 19.

21. (canceled)

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

[0024] FIG. 1 shows schematic illustrations of (a) a base survey and (b) a monitor survey being performed;

[0025] FIG. 2 illustrates the effects of time-lapse changes on ray paths and travel times 1) for (a) near offsets and (b) far offsets; and

[0026] FIG. 3 is a flowchart describing a method according to an embodiment of the invention.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0027] Referring initially to FIGS. 1(a) and 1(b) there is illustrated a subsurface volume 10 comprising a reservoir containing hydrocarbons 12. A survey vessel 16 may comprise a transmitter 18, such as an acoustic transmitter, and an array of receivers 20. The survey vessel may perform a survey by travelling over the subsurface volume 10, transmitting seismic signals into the subsurface volume and subsequently receiving some or all of said seismic signals (the received seismic signals) after reflection of said seismic signal off reflectors 14 within the subsurface volume. A reflector 14 may comprise an interface between layers of contrasting acoustic, optical or electromagnetic properties. The first survey or initial survey, FIG. 1(a), may be referred to as a base survey and is typically performed in the exploration phase before production begins.

[0028] The base survey of the subsurface volume 10 provides first seismic data at a first time T, comprising a set of seismic traces. 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.

[0029] One or more wells 22 may be drilled in order to extract the hydrocarbons 12. As the subsurface volume 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 within the subsurface volume. Thus when a further survey is carried out, FIG. 1(b), 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.

[0030] Thus reservoir monitoring performs a second survey, which may be referred to as a monitor survey, of the subsurface volume 10. The monitor survey of the subsurface volume 10 provides second seismic data at a second time T+ΔT, comprising a set of seismic traces. In the simplest assumption, ΔT 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 ΔT 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.

[0031] 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.

[0032] Non-linear inversion, often referred to as warping, is a known technique for estimating time shifts from 4D base and monitor survey data. An iterative solver is used to find a minimum (with respect to time shifts) of an objective function of the form:

[00001]$E=\underset{j=1}{\overset{N}{.Math.}}.Math.{\left(b\left(t_j\right)-m\left(t_j+\Delta .Math..Math.t_j\right)\right)}^2$

is found using an iterative solver. The vector b(t.sub.i) contains base samples recorded at the discrete times of vector t.sub.j. The vector m(t.sub.j+Δt.sub.j) contains monitor samples recorded at the discrete times of vector t.sub.j plus the vector of time shifts, Δt.sub.j. The time sample index is j. To solve for Δt.sub.j involves linearizing the objective function around a Δt.sub.j estimate and solving for an incremental Δt.sub.j update. Linearization requires calculating the objective function derivative at each iteration.

[0033] Seismic data is acquired over a wide range of transmitter-receiver offsets. The offset is the horizontal distance between transmitter and receiver (e.g. the distance along the surface). In conventional practice, data over the full range of offsets are split into “data gathers”, or ranges. The number of gathers depends on many factors, for example: the acquisition parameters (in particular of the maximum offset that is acquired), the depth, and the seismic data quality. This data is often referred to as pre-stack data as it has not be processed using a technique known as stacking. Stacking is the summing of the seismic traces of this pre-stack data to improve the signal-to-noise ratio. Traces from different shotpoints with a common midpoint are migrated and stacked to form a single trace (post-stack data) during seismic processing. This stacking process gives the data the appearance of zero-offset data.

[0034] Estimating time shifts, and time strains, from time-lapse seismic data is usually done post-stack. This helps attenuate noise which may otherwise contaminate the estimate. However, it also discards all offset-dependent information and necessitates the use of invalid zero-offset ray path and zero-dip assumptions: that seismic waves have travelled on a vertical path through the subsurface, reflecting from horizontal reflectors. Such implementations rely on migration repositioning the recorded energy to its equivalent zero-offset location during processing. Stacking seismic data discards ray path information which is retained in the pre-stack data. The methods proposed herein intend to exploit this ray path information.

[0035] Time shift and time strain estimation methods which avoid the zero-offset assumption are in their infancy, but have already shown potential to yield extra information about time-lapse stress distribution and consequent velocity change anisotropy. Landro and Stammeijer, Geophysics, v. 69 949-957, (2004), Roste et al. Geophysics, v. 71 S219-S234 (2006), Herwanger et al. Annual Meeting SEG, Extended Abstracts (2007) and Herwanger and Horne, Geophysics, v. 74 W13-W33 (2009) have all used time shifts estimated from non-zero offsets to infer information about time-lapse stress distribution, compaction and/or anisotropic velocity changes. However, their methods are restricted to horizontally layered media. Fuck et al. Geophysics, v. 74 D25-D36 (2009) have developed a practical ray-based theory for modelling time shifts in heterogeneous, anisotropic media, but few have attempted to estimate time shifts and time strains from time-lapse seismic data in this way. One such example is given by Vesnaver et al. Geophysics, v. 68 815-823 (2003) in which a coupled tomography system is used to ensure there are no time shifts between time-lapse vintages away from the reservoir. In doing so they account for seasonal overburden variations unrelated to production through tomographic updates to vintage-dependent velocity models. Their methodology was successful in processing away unwanted time shifts, but was apparently not used to estimate true production related time-lapse effects.

[0036] Landro and Stammeijer (2004) showed that time shifts estimated from near and far offsets could be used to estimate simultaneously isotropic velocity changes and thickness changes. Their approach is derived from considering a single horizontal layer, which is statically and dynamically isotropic with regards to velocity. Roste et al. Geophysics, v. 71 S219-S234 (2006) extend this method to handle the full pre-stack data, but do not address the aforementioned assumptions. Landro and Stammeijer (2004) recognise that assuming dynamic velocity isotropy is a serious limitation of their method and Herwanger et al. Annual Meeting SEG, Extended Abstracts (2007) and Herwanger and Horne Geophysics, v. 74 W13-W33 (2009) have shown that very small dynamic velocity anisotropy invalidates the predictions of these methods. The latter authors suggest that ray tracing travel times is necessary to better understand the dependence of travel time on offset.

[0037] This technique disclosed herein falls into the category of non-linear inversion for estimating changes in a model parameter (e.g. estimation of time shifts, time strain or velocity change), but generalizes the inversion for use with non-zero offset pre-stack time-lapse seismic data. The inverse problem becomes that of finding the velocity changes which describe the pre-stack time shifts between time-lapse vintages. The link between velocity changes and pre-stack (offset-dependent) time shifts is:

[00002]$\begin{array}{cc}\Delta .Math..Math.t=\underset{i=1}{\overset{N}{.Math.}}.Math.\frac{d_{2,i}}{v_{2,i}}-\underset{j=1}{\overset{M}{.Math.}}.Math.\frac{d_{1,j}}{v_{1,j}}& \left(1\right)\end{array}$

Equation (1) states that an offset-dependent time shift is equal to the difference in summation of travel times along ray paths from image points common to the two vintages (denoted indices 1 and 2). The indices i and j relate to the index of each cell through which the rays travel within the two velocity fields (through a total of N cells and M cells respectively). The ratio of distance, d, travelled by a ray though a cell and velocity, v, within that cell and along that direction provides the travel time of the ray through the cell.

[0038] FIG. 2 illustrates the situation described by Equation (1) for both (a) near offsets and (b) far offsets which shows the effects of time-lapse changes on ray paths and travel times. The solid line d.sub.1,j is the base ray path and the broken line d.sub.2,i is the monitor ray path. Also shown is reflector R off which the base ray reflects as it travels from source to receiver, and same reflector R′ off which the monitor ray reflects as it travels from source to receiver. As can be seen in this example the position of reflector R, R′ has shifted between the base and monitor surveys, although this will not always be the case. The shaded box A represents a velocity anomaly which appears at some point in time between the base and monitor surveys. The differences between the base ray path d.sub.1,j and the monitor ray path d.sub.2,i are caused by the presence of the velocity anomaly as well as the change in position of the reflector R, R′.

[0039] The travel time of the base ray t.sub.1 and monitor ray t.sub.2 in the specific example illustrated in FIG. 2 (where M=4 and N=4) is:

[00003]$\begin{array}{cc}{t}_1=\frac{d_{1,1}}{v_{1,1}}+\frac{d_{1,2}}{v_{1,2}}+\frac{d_{1,3}}{v_{1,3}}+\frac{d_{1,4}}{v_{1,4}}=\underset{j=1}{\overset{4}{.Math.}}.Math.\frac{d_{1,j}}{v_{1,j}};.Math..Math.t_2=\frac{d_{2,1}}{v_{2,1}}+\frac{d_{2,2}}{v_{2,2}}+\frac{d_{2,3}}{v_{2,3}}+\frac{d_{2,4}}{v_{2,4}}=\underset{i=1}{\overset{4}{.Math.}}.Math.\frac{d_{2,i}}{v_{2,i}}& \left(2\right)\end{array}$

[0040] As can be seen, Equation (1) is dependent upon the distance each ray travels and therefore, to be solved properly, each ray has to be retraced for each iteration. This is computationally expensive and therefore impractical for most applications. To simplify the calculations, as previously mentioned, the time shifts are typically estimated post-stack using zero-offset and zero-dip assumptions.

[0041] Instead, the methodology proposed herein proposes to assume that the ray path geometries do not change between vintages (i=j and M=N):

[00004]$\begin{array}{cc}{t}_1=\underset{j=1}{\overset{M}{.Math.}}.Math.\frac{d_{1,j}}{v_{1,j}};.Math.t_2=\underset{i=1}{\overset{N}{.Math.}}.Math.\frac{d_{2,i}}{v_{2,i}}\approx \underset{j=1}{\overset{M}{.Math.}}.Math.\frac{d_{1,j}}{v_{2,j}}& \left(3\right)\\ \Delta .Math..Math.t=\underset{i=1}{\overset{N}{.Math.}}.Math.\frac{d_{2,i}}{v_{2,i}}-\underset{j=1}{\overset{M}{.Math.}}.Math.\frac{d_{1,j}}{v_{1,j}}\approx \underset{j=1}{\overset{M}{.Math.}}.Math.\left(\frac{d_{1,j}}{v_{2,j}}-\frac{d_{1,j}}{v_{1,j}}\right)& \left(4\right)\end{array}$

[0042] Dividing Equation (4) by

[00005]$t_{1,j}=\frac{d_{1,j}}{v_{1,j}}$

and rearranging in terms of velocity changes (Δv=v.sub.2,j−v.sub.1,j) yields:

[00006]$\begin{array}{cc}\frac{\Delta .Math..Math.t_j}{t_{1,j}}\approx \frac{d_{1,j}.Math.v_{1,j}}{v_{2,j}.Math.d_{1,j}}-\frac{d_{1,j}.Math.v_{1,j}}{v_{1,j}.Math.d_{1,j}}=\frac{v_{1,j}}{v_{2,j}}-1=\frac{-\Delta .Math..Math.v_j}{v_{2,j}}& \left(5\right)\\ \Delta .Math..Math.t\approx \underset{j=1}{\overset{M}{.Math.}}.Math.t_{1,j}.Math.\frac{-\Delta .Math..Math.v_j}{v_{2,j}}& \left(6\right)\end{array}$

[0043] It can now be seen that Equation (6) is no longer dependent on the distance travelled by the ray. The assumption that the ray path geometries do not change between vintages implies that velocity changes and geomechanical compaction/dilation must be small enough to leave the ray paths unchanged to within the accuracy that ray theory approximates seismic wave energy propagation. Seismic waves are extremely band limited and have very large wavelengths. While FIG. 2 illustrates the base and monitor ray paths as narrow; in reality the signals will propagate broadly over a much larger volume than illustrated and roughly conically. Consequently there is significant overlap between the base ray path d.sub.1,j and the monitor ray path d.sub.2,i which is not apparent in FIG. 2, and which makes the proposed approximation valid, to first order.

[0044] While Equation (6) describes only a single ray, it can be used to form a large linear tomographic system for all rays, which can then be inverted for the fractional change in velocity given a set of pre-stack time shifts. The forward modelling operator is a matrix of travel times per cell for every ray path, modelled by ray tracing through an initial velocity field (such as that used to migrate both datasets). In matrix notation this can be written as:

Δt=B−Δν/ν  (7)

[0045] Equation (7) describes a linear tomographic system linking offset-dependent time shifts, Δt, to model parameters, −Δν/ν, via the tomography matrix, B, containing the travel time per cell of each ray. Equation (7) can be used to solve for velocity changes given an existing set of offset-dependent time shifts, if available. Alternatively Equation (7) can be incorporated into the non-linear inversion for time-shifts disclosed in patent application GB1412400.2, as described below.

[0046] The time-lapse linear tomography scheme described by equation (7) has the objective function to be minimized with respect to time shifts:

ε=∥Δt−B−Δν/ν∥.sub.2  (8)

[0047] As can be seen from Equations (6) to (8), the tomographic method according to an embodiment described herein is distinguished from other such approaches by directly solving for velocity change, rather than solving for absolute velocity for each time-lapse vintage. This is accomplished by assuming the changes in ray path geometries between vintages have a minor effect on travel times compared to the change in velocity. Specifically, it is assumed that the distance travelled per cell is equal for rays common across vintages, so that this distance can be eliminated from the tomography (equations (3) to (7)). Since the ray path is assumed not to be a function of the model parameters, the new tomography scheme (equation (8)) is linear and the computationally intensive ray tracing needs to be performed only once.

[0048] The use of ray theory for approximating the paths travelled by seismic energy through the subsurface allows for a tomographic approach to estimate velocity changes. This eliminates the need for the zero-offset assumption, thereby generalizing the classical time shift inverse problem from stacked to pre-stack time-lapse seismic. Stacking seismic data discards ray path information which is retained in the pre-stack data. The extra information can be exploited by this tomographic method to better understand production induced subsurface changes, such as stress distribution, compaction, anisotropic velocity changes, etc. It also removes the requirement of zero-dip, so can be used in structurally complex regions.

[0049] In Equations (6) to (8) velocity changes are assumed to be isotropic. Such an assumption may limit the accuracy of the predictions obtained by this method. This assumption can optionally be removed by decomposing compressional velocity as a function of ray path angle, θ. From Thomsen (1986) anisotropic seismic velocity can be expressed as:

ν(θ)≈α(1+δ sin.sup.2θ cos.sup.2θ+ε sin.sup.4θ)  (9)

[0050] Assuming the velocity changes are elliptically anisotropic (δ=ε) and inserting Equation (9) into Equation (6) gives:

[00007]$\begin{array}{cc}\Delta .Math..Math.t\approx \underset{j=1}{\overset{M}{.Math.}}.Math.\frac{t_{1,j}}{1+{\mathrm{.Math.}}_{2,j}.Math.{\sin }^2.Math.{\theta }_j}.Math.\left(\frac{{\alpha }_{1,j}-{\alpha }_{2,j}}{{\alpha }_{2,j}}+\frac{{\alpha }_{1,j}.Math.{\mathrm{.Math.}}_{1,j}-{\alpha }_{2,j}.Math.{\mathrm{.Math.}}_{2,j}}{{\alpha }_{2,j}}.Math.{\sin }^2.Math.{\theta }_j\right)& \left(10\right)\end{array}$

[0051] Assuming elliptically anisotropic velocity changes allows the tomographic scheme to be altered to solve for the vertical and horizontal components of velocity change. No similar method currently exists, despite the importance placed on understanding such changes in geomechanical modelling.

[0052] Substituting in the definitions of the change in vertical velocity, Δν.sup.0=α.sub.2−α.sub.1, and the change in horizontal velocity, Δν.sup.90=α.sub.2ε.sub.2−α.sub.1ε.sub.1+Δν.sup.0, and assuming the changes are small, yields:

[00008]$\begin{array}{cc}\Delta .Math..Math.t\approx \underset{j=1}{\overset{M}{.Math.}}.Math.t_{1,j}.Math.{\cos }^2.Math.{\theta }_j.Math.\frac{-\Delta .Math..Math.v_j^0}{v_{2,j}^0}+t_{1,j}.Math.{\sin }^2.Math.{\theta }_j.Math.\frac{-\Delta .Math..Math.v_j^{90}}{v_{2,j}^{90}}& \left(11\right)\end{array}$

[0053] Equation (11), like Equation (6), describes only a single ray. In a similar manner to the isotropic example, Equation (11) can be used to form a large linear tomographic system for all rays, which in this anisotropic embodiment can be written as:

[00009]$\begin{array}{cc}\Delta .Math..Math.t=B\left[\begin{array}{c}\frac{-\Delta .Math..Math.v^0}{v}\\ \frac{-\Delta .Math..Math.v^{90}}{v}\end{array}\right]& \left(12\right)\end{array}$

[0054] Equation (12) describes a linear tomographic system linking offset-dependent time shifts, Δt, to elliptically anisotropic velocity changes

[00010]$\left(\mathrm{the}.Math..Math.\mathrm{model}.Math..Math.\mathrm{parameters}.Math..Math.\frac{-\Delta .Math..Math.v^0}{v},\frac{-\Delta .Math..Math.v^{90}}{v}\right)$

via the tomography matrix, B, containing the travel time per cell of each ray. This matrix contains the travel time per cell of each ray scaled by the squared sine or cosine of the ray angle in that cell. As with Equation (7), Equation (12) can be used to solve for velocity changes given an existing set of offset-dependent time shifts or can be incorporated into the non-linear inversion for time-shifts disclosed in patent application GB1412400.2, as described below.

[0055] The time-lapse anisotropic linear tomography scheme described by Equation (12) has the objective function to be minimized with respect to time shifts:

[00011]$\begin{array}{cc}\mathrm{.Math.}={.Math.\Delta .Math..Math.t-B\left[\begin{array}{c}\frac{-\Delta .Math..Math.v^0}{v}\\ \frac{-\Delta .Math..Math.v^{90}}{v}\end{array}\right].Math.}_2& \left(13\right)\end{array}$

[0056] If multi-azimuth time-lapse seismic data is available it is also possible to solve for azimuthally varying estimates of anisotropic velocity changes. This could be sensitive to overburden fracturing, which is important for health, safety and environmental (HSE) interests, integrity and well-planning. Given time-lapse volumes at n different azimuths, the tomographic system described by Equation (12) becomes:

[00012]$\begin{array}{cc}\left[\begin{array}{c}\Delta .Math..Math.t_1\\ \Delta .Math..Math.t_2\\ \mathrm{.Math.}\\ \Delta .Math..Math.t_n\end{array}\right]=\left[\begin{array}{cccc}{B}_1& 0& \mathrm{.Math.}& 0\\ 0& B_2& \ddots & \mathrm{.Math.}\\ \mathrm{.Math.}& \ddots & \ddots & 0\\ 0& \mathrm{.Math.}& 0& B_n\end{array}\right]\left[\begin{array}{c}\Delta .Math..Math.v_1\\ \Delta .Math..Math.v_2\\ \mathrm{.Math.}\\ \Delta .Math..Math.v_n\end{array}\right]& \left(14\right)\end{array}$

[0057] Equation (14) uses the notation

[00013]$\Delta .Math..Math.v_n=\left[\begin{array}{c}\frac{-\Delta .Math..Math.v_j^0}{v_j}\\ \frac{-\Delta .Math..Math.v_j^{90}}{v_j}\end{array}\right]$

for clarity. The index subscript n corresponds to a specific azimuth volume.

[0058] The time-lapse anisotropic linear tomography scheme described by this equation has the objective function to be minimized with respect to time shifts:

[00014]$\begin{array}{cc}\mathrm{.Math.}={.Math.\left[\begin{array}{c}\Delta .Math..Math.t_1\\ \Delta .Math..Math.t_2\\ \mathrm{.Math.}\\ \Delta .Math..Math.t_n\end{array}\right]-\left[\begin{array}{cccc}{B}_1& 0& \mathrm{.Math.}& 0\\ 0& B_2& \ddots & \mathrm{.Math.}\\ \mathrm{.Math.}& \ddots & \ddots & 0\\ 0& \mathrm{.Math.}& 0& B_n\end{array}\right]\left[\begin{array}{c}\Delta .Math..Math.v_1\\ \Delta .Math..Math.v_2\\ \mathrm{.Math.}\\ \Delta .Math..Math.v_n\end{array}\right].Math.}_2& \left(15\right)\end{array}$

[0059] A further advantage of this tomographic method is that it can be incorporated into the non-linear inversion for time shifts disclosed in patent application GB1412400.2. This allows inversion for a change in velocity directly from the pre-stack time-lapse data, without the need to first estimate pre-stack time shifts. This ensures a direct link between the recorded data and desired model parameter throughout the inversion, which is beneficial for understanding the model sensitivity to the data and for developing sensible regularization schemes.

[0060] Patent application GB1412400.2 describes a fast non-linear inversion method for estimating 4D time shifts and time strain. In this method, 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.

[0061] The method described in GB1412400.2 comprises:

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

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

[0064] 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

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

[0066] 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.

[0067] The method may comprise the steps of:

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

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

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

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

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

[0073] Plane wave destruction is taught by Sergey Fomel, in Geophysics, 67, 1946-1960 (2002)—herein incorporated by reference) From the plane wave description of a trace it is possible to use shifts σ, to parameterise a filter L(σ), that warps one trace to match another trace. The filter L(σ) destroys the energy difference between plane waves t.sub.1 and t.sub.2, and is consequently referred to as a “plane wave destruction (PWD) filter”. The system: t.sub.1−L(σ)t.sub.2=0, can be solved as a non-linear inversion for shifts, σ.

[0074] Modifying the objective function taught in GB1412400.2 for use with non-zero offset pre-stack time-lapse seismic data, the (unregularized) objective function is:

ε=∥L(−Δt)s.sub.1−L(Δt)s.sub.2∥.sub.2  (16)

[0075] The vectors s.sub.1 and s.sub.2 contain all pre-stack traces of two time-lapse vintages. The function L is a non-stationary phase-only filter parameterised by time shifts, Δt. Application of this filter shifts the target trace according to the array of time shifts. It should be noted that Equation (16) 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. The filter has an analytic expression derived from the plane wave model of a seismic trace, and has an analytic gradient. The filter may be a three term PWD filter, and the three-term filter coefficients b.sub.−1, b.sub.0, b.sub.1 of the filter may be described by:

[00015]$\begin{array}{cc}{b}_{-1}=\frac{\left(t-\Delta .Math..Math.t\right).Math.\left(2.Math.t-\Delta .Math..Math.t\right)}{12.Math.t^2}.Math..Math.b_0=\frac{\left(2.Math.t+\Delta .Math..Math.t\right).Math.\left(2.Math.t-\Delta .Math..Math.t\right)}{6.Math.t^2}.Math..Math.b_1=\frac{\left(t+\Delta .Math..Math.t\right).Math.\left(2.Math.t+\Delta .Math..Math.t\right)}{12.Math.t^2}& \left(17\right)\end{array}$

[0076] The coefficient t is the time sampling rate of the data and allows the local time shifts, Δt, to be expressed in real time units.

[0077] Knowing the filter coefficients b.sub.−1, b.sub.0, b.sub.1 allows a non-linear inversion for Δt to be constructed. The objective function given in Equation (16) implies finding the least-squares solution to the non-linear inverse problem:

[00016]$\begin{array}{cc}\left[\begin{array}{cc}L\left(-\Delta .Math..Math.t\right)& -L\left(\Delta .Math..Math.t\right)\end{array}\right]\left[\begin{array}{c}{s}_1\\ s_2\end{array}\right]=0& \left(18\right)\end{array}$

[0078] Making the substitutions:

[00017]$C_p=\left[\begin{array}{cc}L\left(-\Delta .Math..Math.t_p\right)& -L\left(\Delta .Math..Math.t_p\right)\end{array}\right],.Math.d_p=-C_p\left[\begin{array}{c}{s}_1\\ s_2\end{array}\right].Math..Math.\mathrm{and}$$G_p=\mathrm{diag}.Math.\left\{C_p^{\prime }\left[\begin{array}{c}{s}_1\\ s_2\end{array}\right]\right\}$

and linearizing Equation (18) gives:

G.sub.pδt.sub.p+1=d.sub.p  (19)

[0079] Equation (19) is the linearization of Equation (18) and is solved iteratively for an incremental time shift perturbation, δt.sub.p+1, given an initial model, Δt.sub.p, where p is the iteration number. The matrix C′ is the Jacobian of C, constructed from the partial derivatives of the filter coefficients with respect to Δt. Each iteration the inverted δt.sub.p+1 is added to the previous time shift estimate, Δt.sub.p, to give an updated set of time shifts, Δt.sub.p+1. This is used to rebuild G and d for the subsequent iteration.

[0080] Iterating on Equation (19) provides local time shifts which align the time-lapse vintages. For use with stacked data, it is possible to insert a z-direction integration matrix (weighted by the time sampling), A, into Equation (18), thereby changing the model space from time shifts to velocity change,

[00018]$\frac{-\Delta .Math..Math.v}{v}.Math.\text{:}$

[00019]$\begin{array}{cc}{G_p.Math.A.Math..Math.\delta .Math.\frac{-\Delta .Math..Math.v}{v}.Math.}_{p+1}=d_p& \left(20\right)\end{array}$

[0081] The matrix A simply integrates travel time along the vertical ray path, assumed valid for stacked data. Replacing matrix A with matrix B from Equation (7), Equation (12) or Equation (14) generalizes the method by changing this z-direction travel time integration with travel time integration along the traced ray paths:

[00020]$\begin{array}{cc}{G_p.Math.B.Math..Math.\delta .Math.\frac{-\Delta .Math..Math.v}{v}.Math.}_{p+1}=d_p& \left(21\right)\end{array}$

[0082] Iterating on equation (21) provides velocity changes which describe the pre-stack time shifts between time-lapse vintages. This system can be regularized using standard techniques in order to estimate velocity changes desirable for geomechanical modelling.

[0083] A similar treatment can be performed on the anisotropic embodiment. In this case, Equation (21) becomes:

[00021]$\begin{array}{cc}{G_p.Math.B.Math..Math.\delta \left[\begin{array}{c}\frac{-\Delta .Math..Math.v^0}{v}\\ \frac{-\Delta .Math..Math.v^{90}}{v}\end{array}\right].Math.}_{p+1}=d_p& \left(22\right)\end{array}$

[0084] A similar treatment can be performed on the azimuthally anisotropic embodiment described above. In this case, Equation (21) becomes:

[00022]$\begin{array}{cc}{G_p\left[\begin{array}{cccc}{B}_1& 0& \mathrm{.Math.}& 0\\ 0& B_2& \ddots & \mathrm{.Math.}\\ \mathrm{.Math.}& \ddots & \ddots & 0\\ 0& \mathrm{.Math.}& 0& B_n\end{array}\right]\left[\begin{array}{c}\delta .Math..Math.\Delta .Math..Math.v_1\\ \delta .Math..Math.\Delta .Math..Math.v_2\\ \mathrm{.Math.}\\ \delta .Math..Math.\Delta .Math..Math.v_n\end{array}\right].Math.}_{p+1}=d_p& \left(23\right)\end{array}$

[0085] Equation (23) uses the notation

[00023]$\Delta .Math..Math.v_n=\left[\begin{array}{c}\frac{-\Delta .Math..Math.v_j^0}{v_j}\\ \frac{-\Delta .Math..Math.v_j^{90}}{v_j}\end{array}\right]$

for clarity. The index subscript n corresponds to a specific azimuth volume.

[0086] FIG. 3 is a flowchart illustrating an embodiment of the techniques described herein. At step 300 a first survey, or base survey, is performed during which first seismic data 310 is acquired. At step 320 a second survey, or monitor survey, is performed during which second seismic data 330 is acquired. At step 340, a model of the initial velocity field is determined from standard 4D seismic data processing techniques to be used for seismic data migration. This may be optimized for migrating a particular 4D vintage (often the first seismic data), but could be any velocity model suitable for migrating all 4D vintages. Using the model of the velocity field, all sets of data are migrated using known techniques (for example, Kirchhoff Migration).

[0087] At step 350, a ray tracing operation is performed through the velocity field. In this example, the ray tracing is performed on the velocity field model used to migrate the first and second seismic data 310, 330, although a different velocity field model can be used. The ray tracing operation may use ray theory to approximate the paths travelled by seismic energy through the subsurface volume. The ray tracing needs only to be performed once as it is assumed that the approximated ray paths will be the same for both sets of seismic data. The ray tracing yields a model of travel times per cell for every ray path. This can be represented as a matrix. This matrix may be matrix B as used in the Equations above, specifically: Equation (7), Equation (12), Equation (14), Equation (21), Equation (22) or Equation (23).

[0088] At step 360, the inversion is performed to obtain estimates of changes having occurred between said base survey and monitor survey in terms of at least one model parameter. The matrix B determined at step 350 is used as the forward modelling operator in the inversion. The at least one model parameter may be (for example) propagation velocity or travel time of said seismic signal, and therefore the estimated changes may be in terms of time shift, time strain or velocity change. The inversion may be performed using any of the techniques disclosed herein, for instance using the objective function described by Equation (8), by Equation (13) or by Equation (15) or by iterating on Equation (21), Equation (22) or Equation (23). The inversion may be regularized as is well known in the art.

[0089] 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.

[0090] The concepts described herein find utility in all aspects (real time or otherwise) of surveillance, monitoring, optimisation and prediction of hydrocarbon reservoir and well systems, and may aid in, and form part of, methods for extracting hydrocarbons from such hydrocarbon reservoir and well systems.

[0091] 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.