Method of constraining seismic inversion

Abstract

Disclosed is a method a seismic inversion for petrophysical properties of a subsurface volume comprising the steps of: obtaining petrophysical data relating to valid geological and/or dynamical scenarios, converting this data into valid combinations of elastic parameters; projecting the valid combinations of elastic parameters onto a spherical plot; and determining a penalty term from the distances between each cell of the spherical plot and the nearest valid combination of elastic parameters within the subsurface volume. Valid geological and/or dynamical scenarios comprise those which are petrophysically possible. The penalty term is then used to constrain an inversion minimizing a cost function associated with seismic mismatch between two or more seismic surveys.

Claims

1. A method of performing a geometric inversion of seismic data comprising the steps of: obtaining petrophysical data relating to valid geological and/or dynamical scenarios within a subsurface volume comprising a hydrocarbon reservoir or region thereof, wherein valid geological and/or dynamical scenarios comprise those which are petrophysically possible; converting said petrophysical data into valid combinations of elastic parameters; projecting said valid combinations of elastic parameters onto a spherical plot; determining a penalty term from the distances between each cell of the spherical plot and the nearest valid combination of elastic parameters; performing an inversion which inverts for changes in dynamic properties of the hydrocarbon reservoir by minimizing a cost function associated with seismic mismatch between two or more seismic surveys and at least a synthetic dataset computed from said elastic parameters; wherein said penalty term is used to constrain said inversion; and using the result of said inversion in predicting performance of said hydrocarbon reservoir.

2. A method as claimed in claim 1 wherein the projection of said valid combinations of elastic parameters onto a spherical plot comprises projection onto the surface of a plotted sphere, and the determined distance between each cell of the spherical plot and the nearest valid combination of elastic parameters is the great circle distance.

3. A method as claimed in claim 1 wherein, when the inversion is a 3D inversion, one of the surveys is the synthetic dataset computed from the inverted parameters.

4. A method as claimed in claim 1 comprising determining the penalty term for only a region of said subsurface volume, wherein said valid geological and/or dynamical scenarios and therefore said valid combinations of elastic parameters are those which are valid for that particular region.

5. A method as claimed in claim 4 comprising determining different penalty terms for different regions within the subsurface volume.

6. A method as claimed in claim 4 wherein each region comprises a different layer of the subsurface volume.

7. A method as claimed in claim 4 comprising the step of using knowledge of the predominant fluid effect in a region in determining valid geological and/or dynamical scenarios for that region.

8. A method as claimed in claim 1 comprising the step of using a simulation to determine said valid geological and/or dynamical scenarios within the subsurface volume.

9. A method as claimed in claim 8 wherein said simulation is a MonteCarlo simulation.

10. A method as claimed in claim 1 comprising the step of using a rock physics/petro-elastic model to convert said petrophysical data into valid combinations of elastic parameters.

11. A method as claimed in claim 1 comprising the step of determining valid combinations of elastic parameters directly from observation.

12. A method as claimed in claim 1 wherein said method comprises 4D modelling of said subsurface volume and the spherical plot is a plot of changes over time in any three elastic parameters, one on each axis.

13. A method as claimed in claim 1 wherein said method comprises 3D modelling of said subsurface volume and the spherical plot is a plot of any three elastic parameters, with one elastic parameter on each axis.

14. A method as claimed in claim 13 comprising determining an origin of the spherical plot of elastic values wherein said determination comprises: grouping said petrophysical data according to their facies or seismic-facies characteristics; calculating average values for the elastic parameters being plotted, from the values of each of these parameters for each group of facies characteristic; and locating said origin at the point defined by said calculated average values for each parameter.

15. A method as claimed in claim 12 wherein said elastic parameters comprise any three of: p-wave velocity or impedance, s-wave velocity or impedance, density and any Lam parameter.

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

17. 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.

18. A computer program carrier comprising the computer program of claim 17.

19. Apparatus specifically adapted to carry out the steps of the method as claimed in claim 1.

20. A method as claimed claim 1 comprising using the result of said inversion in determining future well placement.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

(2) FIG. 1 is a graph of reflectivity (y-axis) against angle (x-axis) showing AVO curves for an increase in V, for both the best determined and worst determined combinations of parameters;

(3) FIG. 2 is a flowchart illustrating a method according to an embodiment of the invention;

(4) FIG. 3 is a plot of valid combinations of parameters in spherical space;

(5) FIG. 4 is a plot of a cost function resultant from the plotted valid combinations of parameters in FIG. 3; and

(6) FIG. 5 shows spherical plots of valid elastic parameter combinations for a number of scenarios: (a) a random scenario; (b) water replacing oil; (c) gas replacing water; and (d) pressure changes only.

DETAILED DESCRIPTION OF THE EMBODIMENTS

(7) Amplitude-versus-offset (AVO) inversion is an ill-posed problem and it is commonly admitted that only two parameters can be extracted from the AVO. More recently, it has been shown that not only is 4D AVO an ill-posed problem, but there is also a large null space in the inversion. This means that, with the same elastic data, a large ensemble of rock physics solutions are equivalent. The first common way of by-passing the problem is to perform the inversion in the rock physics domain. This approach relies on an accurate Petro Elastic Model (PEM) and therefore on prior knowledge of certain petrophysical data (e.g. porosity) throughout the reservoir. Another method is to use a Bayesian inversion, but the weighting of the prior can be extremely cumbersome.

(8) To illustrate the problem, consider the pre-stack inversion AVO term R(). This is a very ill-posed and badly determined problem (considering here the 3D situation):

(9) $R () = ((\frac{1}{2} + \frac{1}{2} \tan^{2}) \frac{V_{P}}{V_{P}} + (\frac{1}{2} - 2^{2} \sin^{2}) \frac{}{} - 4^{2} \sin^{2} \frac{V_{s}}{V_{s}})$
and therefore:

(10) $R () ((\frac{1}{2} + \frac{1}{2} \tan^{2}) \frac{(V_{P})}{{\hat{V}}_{P}} + (\frac{1}{2} - 2^{2} \sin^{2}) \frac{()}{\hat{}} - 4^{2} \sin^{2} \frac{(V_{s})}{{\hat{V}}_{s}})$
where is the seismic angle of incidence, V.sub.p/V.sub.p is the relative change of the p-wave velocity at the interface between two layers (e.g. sand/shale), V.sub.s/V.sub.s, is the relative change of the s-wave velocity at the interface between two layers, / is the relative change of the density at the interface between two layers, is V.sub.s/V.sub.p (i.e. the mean value) and {circumflex over (V)}.sub.p, {circumflex over ()} and {circumflex over (V)}.sub.s are the mean values of V.sub.p, V.sub.s and respectively (as relative changes are being estimated).

(11) The skilled person will appreciate that other equations are possible which are not functions of , V.sub.p, V.sub.s, but other combinations of elastic parameters e.g. Impedances I.sub.p I.sub.s, or , , (Lam parameters). The skilled person will further appreciate that it is possible to compute the constraint in essentially the same way for all sets of parameters.

(12) Solving the system of equations in a linear manner would mean solving for:

(13) $[\begin{matrix} R (_{1}) \\ R (_{2}) \\ .Math. \\ R (_{n}) \end{matrix}] = [\begin{matrix} \frac{1}{2} + \frac{1}{2} \tan^{2}_{1} & \frac{1}{2} - 2^{2} \sin^{2}_{1} & 4^{2} \sin^{2}_{1} \\ \frac{1}{2} + \frac{1}{2} \tan^{2}_{2} & \frac{1}{2} - 2^{2} \sin^{2}_{2} & 4^{2} \sin^{2}_{2} \\ .Math. & .Math. & .Math. \\ \frac{1}{2} + \frac{1}{2} \tan^{2}_{n} & \frac{1}{2} - 2^{2} \sin^{2}_{n} & 4^{2} \sin^{2}_{n} \end{matrix}] [\begin{matrix} \frac{V_{P}}{V_{P}} \\ \frac{}{} \\ \frac{V_{S}}{V_{S}} \end{matrix}]$
for any number of angles (equal or greater than 3), assuming a constant . This matrix equation can be represented in shorthand as: d=Am.

(14) To analyse properties of the system, a covariance matrix can be calculated:
C.sub.m=(A.sup.TC.sub.d.sup.1A).sup.1

(15) By way of example, say that a data uncertainty or data error is introduced, for example 0.0001 (reflection coefficient units). Considering then an example (best case) of angles of 0-45 degrees with the uncertainty at each angle of 0.0001, for =0.5:

(16) $C_{m} = (\begin{matrix} 0.0033 & - 0.0031 & 0.0045 \\ - 0.0031 & 0.0030 & - 0.0043 \\ 0.0045 & - 0.0043 & 0.0063 \end{matrix})$

(17) The principal axis and vectors can then be determined to find the best determined and worst determined combinations of parameters:

(18) $_{m} = (\begin{matrix} 3.5 10^{- 6} & 0 & 0 \\ 0 & 5.8 10^{- 5} & 0 \\ 0 & 0 & 1.25 10^{- 2} \end{matrix})$ $_{m} = (\begin{matrix} 0.83 & - 0.22 & 0.51 \\ 0.51 & 0.72 & - 0.49 \\ - 0.26 & 0.66 & 0.71 \end{matrix})$
where .sub.m are the eigenvalues of the covariance matrix C.sub.m and .sub.m are the associated eigenvectors. These in turn reveal the worst determined and best determined combinations of parameters.

(19) From the first column of each of these examples, it can be seen that the best determined combinations of parameters are when:

(20) $\frac{}{} = \frac{0.51}{0.83} \frac{V_{p}}{V_{p}} \frac{V_{s}}{V_{s}} = \frac{- 0.26}{0.83} \frac{V_{p}}{V_{p}}$
and from the last column of each of these examples, it can be seen that the worst determined combinations of parameters are when:

(21) $\frac{}{} = \frac{- 0.49}{0.51} \frac{V_{p}}{V_{p}} \frac{V_{s}}{V_{s}} = \frac{- 0.71}{0.51} \frac{V_{p}}{V_{p}}$

(22) FIG. 1 shows calculated reflectivity curves for two values of V.sub.p:V.sub.p=0.1 (thin lines) and V.sub.p=0.2 (thick lines), representing 9 and 10 percent increases respectively. The lines within rectangle 10 comprise plots for the best determined combinations of parameters, while the lines within rectangle 12 comprise plots for the worst determined combinations of parameters. The Figure highlights how solutions of the worst determined combinations of parameters for V.sub.p=0.1 and V.sub.p=0.2 will be virtually indistinguishable from each other (there are two curves shown within rectangle 12 which overlap for most angles).

(23) V.sub.p is highly constrained by the time shift t, near angle and far angle. is determined through near angle and mid angle amplitudes. V.sub.s is determined through mid angle offsets. As the majority of the cost is reduced through the time shift: V.sub.p fits the time shift then compensates amplitudes at near offsets. V.sub.s then compensates for at mid offsets.

(24) As a consequence, unconstrained pre-stack seismic inversions will generally not be adequate to provide realistic combinations of elastic parameters for dynamic inversion.

(25) It is proposed to constrain AVO inversions in a way that is simple, efficient and respects the information contained in the data (not a hard constraint). This approach introduces an additional cost into the inversion process that is a function of the combination of elastic parameters only. The magnitude of the elastic parameters has no impact on the additional cost. By using this domain, loose prior information is provided to the elastic inversion which enables selection of model solutions that are consistent with prior geological and dynamic considerations.

(26) FIG. 2 is a flowchart of a method according to a first embodiment. The approach entails a first main step 200 of creating distance penalties for combinations of elastic parameters. This may be done on a region-by-region basis for a particular reservoir, such that separate penalties are created for each reservoir layer, for example.

(27) Creation of distance penalties 200 may be performed using geological and dynamic information 205 for the reservoir, or region thereof. This information may be obtained from a reservoir model, geological model and/or the well. Alternatively or in addition, the geological and dynamic information can also come from prior knowledge of expected geology and/or expected changes in dynamic properties (a priori assumptions).

(28) The geological and dynamic information is used to create all possible (i.e. valid) situations relative to a geological/dynamical context 210. This may be done using a MonteCarlo simulation, for example. Possible situations are those which could occur in actuality (i.e. those situations that are petrophysically sensible) as opposed to mathematical solutions which are nonsensical in practice. The geological and dynamic information may be applicable to a particular region (e.g. layer) of the reservoir, with the possible situations being determined for that region. This is particularly the case where the predominant fluid interaction for a region is known. Examples of predominant fluid interaction which may be known for a particular region include: water replacing oil, oil replacing water, gas replacing water, water replacing gas, pressure changes only.

(29) These possible/valid solutions are converted that into elastic parameters 215. This may be done using a rock physics model or petro-elastic model (PEM), for example. Using such a model, all expected changes in V.sub.p, V.sub.s and (for example, other elastic parameters may be chosen) expected within a region of the reservoir should be forward modelled. This step may comprise normalising all of the combinations of V.sub.p, V.sub.s and that are bigger than a certain threshold (it would not be expected to detect very small changes in seismic properties (e.g less than 0.5%)).

(30) The combinations of elastic parameters are then projected onto spherical space (step 220) comprised of a 3D spherical plot of changes in V.sub.p against changes in V.sub.s against changes in density p. FIG. 3 shows an example plot 300, comprising calculated valid combinations of elastic parameters 310 expected around a water injector, including a variety of 3D geological properties.

(31) At step 225 a distance penalty is computed by measuring the minimum great circle distance or orthodromic distance (i.e. the distance along the surface of the sphere) between each cell of the sphere (360 azimuths*180 inclinations) to the nearest simulated combination (or valid solution). This may such that, for example, a cell representing a valid solution is attributed a distance of 0, while the maximum distance on the sphere (between two antipodal points) is attributed a distance of 1. A normalisation can be applied to scale the distance between chosen values; this should not simply be a linear transformation if the constraint is to be made sharper. In the 3D example (described below), isolated points may be removed. Other distance computation methods are possible. For example, rather than computing the minimum distance, an average distance to all valid points or an average distance to the n nearest points can be used should the data be noisy.

(32) FIG. 4 shows a resultant cost function 400 produced using the data of FIG. 3. Darker shading indicates lower cost, with conversely the unshaded region being the highest cost.

(33) Output is a distance penalty 230 for the region or layer in question (or for the whole reservoir, if the method is not being performed on a region-by region basis). Distance penalties for other layers can then be calculated, such that a different constraint sphere is assigned to different parts of the model (usually each layer) via an associated penalty.

(34) Once all penalties have been calculated, the elastic penalty can be used to constrain the inversion. During the inversion 240 the penalty term is applied to the cost function associated with seismic mismatch , which is always between 0 and 1 in this embodiment (although it may take other values in different embodiments).

(35) The computation of the constraint may be done in the same process as the inversion. During the inversion, at each evaluation of the cost function, the current inversion parameters (V.sub.p, V.sub.s, ) can be projected onto the sphere, with the penalty term corresponding to the value of the constraint at this position.

(36) FIGS. 5(a)-(d) show spherical plots of valid elastic parameter combinations for a number of scenarios. FIG. 5(a) shows a random scenario (as you may expect to see if performing the methods herein throughout whole reservoir without regionalisation). FIG. 5(b) shows a plot for a region/layer where the predominant effect is water replacing oil. FIG. 5(c) shows a similar plot where the predominant effect is gas replacing water, and FIG. 5(d) shows a plot where there have only been pressure changes. In each case, the dark regions represent valid combinations. These figures illustrate the advantage of determining constraints on a regional or layer-by-layer basis. It is likely that the predominant effect of any given layer is known, and therefore only valid combinations for such a scenario needs to be plotted As a result, FIGS. 5(b) to (d) show far fewer and more localised valid solutions than FIG. 5(a), therefore resulting in a better constrained inversion.

(37) The cost function may take the form of:

(38) $C = {.Math.}_{j = 1}^{M} {.Math.}_{i = 1}^{N} (b (t_{i},_{j}) - m (t_{i} - {.Math.}_{1}^{i} \frac{V_{Pk}}{V_{Pk}},_{j}) - * ((\frac{1}{2} + \frac{1}{2} \tan^{2}_{j}) \frac{V_{Pk}}{V_{Pk}} + (\frac{1}{2} - 2^{2} \sin^{2}_{j}) \frac{_{k}}{_{k}} - 4^{2} \sin^{2} \frac{V_{Sk}}{V_{Sk}})) +^{2} D (\frac{V_{Pk}}{V_{Pk}}, \frac{_{k}}{_{k}}, \frac{V_{Sk}}{V_{Sk}})$

(39) Wherein the first two terms are the base and monitor traces respectively, the third term is the AVO term and the fourth term is the constraint, with D being the calculated distance penalty.

(40) While the examples above relate to 4D seismic inversions, the concepts herein are also applicable to 3D seismic inversions. The main difference in the latter case is the need to locate an origin for the spherical plot as an initial step. With 4D seismic (measuring changes of elastic parameters over time), there is a natural origin representing no change, i.e. the initial conditions at zero time. With 3D seismic data sets there is no such natural origin.

(41) To find the origin, the data (elastic values of the logs) is initially grouped according to facies log, i.e. the data is grouped according to whether it predominately comprises sand/shale/carbon etc. This grouping can be made according to true geological facies or seismically distinguishable facies. The origin for the V.sub.p axis is calculated from the average (midpoint) of the values of V.sub.p for the different groups. The origin for the V.sub.s and axes are similarly found from averages of the values of V.sub.s and for each group. For example, where the facies log shows only sand and shale, the origin will be calculated thusly:

(42) $V_{Porigin} = \frac{V_{Psand} + V_{Pshale}}{2}; V_{Sorigin} = \frac{V_{Ssand} + V_{Sshale}}{2};_{origin} = \frac{_{sand} +_{shale}}{2}$

(43) Once the origin is found, the rest of the method is much the same as that shown in FIG. 2. However, a rock physics model is unnecessary when petro-elastic information can be obtained from the well-log (elastic values have been measured at the well). If it is not the case, synthetic elastic logs can be computed from petro-physics logs using a rock physics model. In this embodiment, valid combinations of parameters V.sub.p, V.sub.s and (not changes in these parameters as with the 4D plots) are projected onto spherical space centred on the calculated origin.

(44) The techniques described herein provide a model driven prestack inversion workflow that inverts for changes in the dynamic properties of the reservoir. The (3D or 4D) prestack elastic inversion is constrained using a combination of rock physics and reservoir engineering information. The workflow also allows for the simultaneous inversion of amplitudes and time-shifts.

(45) The techniques described herein stabilize the elastic inversion, with results showing almost no additional residual energy in the seismic data compared to an unconstrained inversion. This means these constraints provide a cost-equivalent solution which is better suited to the dynamic changes that are expected to be seen in the reservoir. This is extremely beneficial since it is not the elastic parameters that matter in reservoir characterization but rather petrophysical parameters (in 3D) or dynamic parameters (in 4D). These techniques are particularly beneficial for using seismic data in building geomodels or using time lapse seismic in quantitative assisted history matching.

(46) The constraint can be implemented throughout the 3D space or on a layer by layer basis. The latter approach is particularly suited where the predominant fluid effect in each layer is known, i.e. which fluid is being replaced by which fluid in each layer.

(47) By computing more realistic but cost-equivalent inversion results, the inversion for pressure and saturation yields more useful data, in that it is more consistent with the expected changes in reservoir properties. This allows for the interpretation of 4D seismic directly in the dynamic domain (e.g. changes in pressure and water saturation), potentially enabling the detection of by-passed oil, or reservoir compartmentalisation. A better understanding of the dynamic reservoir properties helps improve future well placement and prediction of reservoir performance.

(48) 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.

(49) The concepts described herein find utility in all aspects 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.

(50) 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.

Method of constraining seismic inversion

Assignee

Inventors

Cpc classification

Classification Explorer

G01V2210/6122

PHYSICS

Classification Explorer

G01V1/306

PHYSICS

Classification Explorer

G01V2210/6169

PHYSICS

Classification Explorer

G06F17/18

PHYSICS

Classification Explorer

G01V2210/6242

PHYSICS

Classification Explorer

G01V1/308

PHYSICS

Classification Explorer

G01V20/00

PHYSICS

Classification Explorer

G01V2210/614

PHYSICS

International classification

Classification Explorer

G01V99/00

PHYSICS

Classification Explorer

G01V1/30

PHYSICS

Classification Explorer

G06F17/18

PHYSICS

Abstract

Claims

Description