Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume

1. A method for characterising evolution of a subsurface volume comprising at least one of a hydrocarbon reservoir and well system, said method comprising:

obtaining, via at least one acoustic source and an array of acoustic receivers, a first survey of the subsurface volume with a first set of seismic traces corresponding to a first time;

obtaining, via the at least one acoustic source and the array of acoustic receivers, a second survey of the subsurface volume with a second set of seismic traces corresponding to a second time;

constructing an analytic filter operable to shift at least one seismic trace in dependence of time shift, wherein said analytic filter comprises an all-pass plane wave destruction filter parameterized by time shifts 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;

performing an inversion to obtain time shift estimates for a model parameter such that said analytic filter aligns said first survey and said second survey; and

wherein said time shift estimates characterize the evolution of the subsurface volume comprising at least one of a hydrocarbon reservoir and well system over a time period between the first time and the second time.

2. The method as claimed in claim 1 , wherein said analytic filter is such that said model parameter is a function of said analytic filter and is not a function of at least one of said second set of seismic traces.

3. The method as claimed in claim 1 , wherein said analytic filter comprises a plane wave destruction filter.

4. A method as claimed in claim 1 , wherein said analytic filter 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.

5. The method as claimed in claim 1 , wherein filter coefficients of said analytic filter are non-stationary with respect to time and space coordinates.

6. The method as claimed in claim 1 , wherein said model parameter is at least one of: time shift, time strain or velocity change.

7. The method as claimed in claim 1 , wherein said method comprises the steps of:

constructing said analytical filter and their derivatives from a current model parameter model; and

iteratively updating an initial model parameter estimate with an incremental model parameter update until convergence on a solution.

8. The method as claimed in claim 7 wherein the step of iteratively updating an initial model parameter estimate comprises:

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

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

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

9. The method as claimed in claim 7 , wherein:

said model parameter comprises time shift such that said analytic filter is parameterized by time shifts and the model parameter model is a time shift model, and

a conversion to time strains is performed prior to the inversion, a resultant time strain estimate being converted back to an incremental time shift update and used to update the time shift model.

10. The method as claimed in claim 1 , wherein the inversion is regularized spatially in three directions.

11. The method as claimed in claim 1 , wherein the inversion comprises shifting each of said first survey and said second survey to a common midpoint.

12. The method as claimed in claim 1 , wherein said subsurface volume is divided into regions according to a high order discrete classification attributed in accordance with a particular region's geologic characteristics, and said inversion is constrained per region of subsurface volume.

13. The method as claimed in claim 12 , wherein said regions comprise geological layers and said high order discrete classification comprises facies.

14. The method as claimed in claim 12 , wherein said inversion comprises inverting for coefficients of a regional function at each spatial location within a region, wherein said regional function describes how a property varies within a region.

15. The method as claimed in claim 14 , wherein said coefficients of the regional function are region-dependent coefficients, said region-dependent coefficients being dependent on the region of subsurface volume they correspond to.

16. The method as claimed in claim 14 , wherein said regional function comprises a forward modelling operator imposing a region-dependent constraint, said region-dependent constraint being dependent on the region of subsurface volume.

17. The method as claimed in claim 16 , wherein said forward modelling operator expands the coefficients of the regional function at each spatial location to every seismic sample point within the region.

18. The method as claimed in claim 16 , wherein the region-dependent constraint is that the model parameter is constant within each region.

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

20. The method as claimed in claim 1 , further comprising the step of using the estimates of said model parameter in prediction of future evolution of said subsurface volume.

21. A computer-program product comprising a non-transitory computer-usable medium having computer-readable program code embodied therein, the computer-readable program code adapted to be executed to implement a method comprising the steps of claim 1 .