Method for estimating subsurface properties from geophysical survey data using physics-based inversion

Abstract

A hydrocarbon exploration method for determining subsurface properties from geophysical survey data. Rock physics trends are identified and for each trend a rock physics model is determined that relates the subsurface property to geophysical properties (103). The uncertainty in the rock physics trends is also estimated (104). A geophysical forward model is selected (105), and its uncertainty is estimated (106). These quantities are used in an optimization process (107) resulting in an estimate of the subsurface property and its uncertainty.

Claims

1. A computer-implemented method for estimating a subsurface property, and uncertainty associated with the estimate, from geophysical measurements, said method comprising: determining one or more rock type trends in the subsurface, and identifying a rock physics model for each trend that mathematically relates the subsurface property to one or more geophysical properties; determining a statistic describing uncertainty in each rock type trend and its rock physics model; selecting a geophysical forward model for predicting the geophysical measurements using current values of the subsurface property at a plurality of subsurface spatial locations; determining a statistic describing uncertainty in the geophysical forward model in predicting the geophysical measurements; forming an estimate of the subsurface property and its associated uncertainty, using a computer, which estimate minimizes misfit between the geophysical measurements and forward-modeled predictions of the geophysical measurements, using the one or more rock type trends and their rock physics models, the uncertainties in the rock type trends and in the forward model, and using the estimate of the subsurface property to assess commercial potential of a subsurface reservoir.

2. The method of claim 1 wherein the geophysical measurements are at least one of: seismic data traces, stacked seismic data, angle stacks of seismic data, pre-stack gathers of seismic data, seismic reflection amplitudes, resistivity survey traces, gravity survey traces, passively recorded seismic signals, and any combination thereof.

3. The method of claim 1, wherein the subsurface property is pore fluid density, acoustic velocity in the pore fluid, a fluid factor combining the density and velocity of the pore fluid, porosity of reservoir rock, volume fraction of shale in the reservoir rock, volume fraction of at least one pore type composing the reservoir rock, volume fraction of at least one mineral type composing the reservoir rock, permeability of the reservoir rock, or any combination thereof.

4. The method of claim 1, wherein the geophysical measurements are seismic data, and the rock type trends and the uncertainty in each rock type trend and its rock physics model are extrapolated guided by interpreted horizons in the seismic data.

5. The method of claim 1, further comprising: selecting an information metric that measures quality of the geophysical measurements; using the information metric to reduce amount of measured data with minimal loss of information relevant to the subsurface property.

6. The method of claim 5, wherein the subsurface property is a direct hydrocarbon indicator.

7. The method of claim 1, wherein the minimizing misfit uses a linearized optimization process.

8. The method of claim 7, wherein the linearized optimization process is iterated.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The present invention will be better understood by referring to the following detailed description and the attached drawings in which:

(2) FIG. 1 is a flow chart showing basic steps in one embodiment of the present inventive method;

(3) FIG. 2 is a flow chart showing basic steps in a second embodiment of the present inventive method;

(4) FIG. 3 shows an example of a 2D seismic section containing NMO corrected pre-stack CDP gathers;

(5) FIG. 4A shows an example of a common depth point gather or CDP gather near a well location for the full range of recorded travel time, and FIG. 4B shows the same data for a smaller time window;

(6) FIG. 5A shows an example of well log measurements of compressional wave velocity (Vp), shear wave velocity (Vs), and density (rho) versus depth; and FIG. 5B shows the clustering of the same well log measurements and how the relationship between well log measurements varies with rock type or lithofacies;

(7) FIG. 6 shows the measured reflection amplitude versus offset for the CDP point in FIG. 4 at a given reflector and the corresponding predicted reflection amplitude for a given choice of earth model properties;

(8) FIG. 7 shows the compressional wave velocity and density of the fluid at depth points within the time window shown in FIG. 4B as estimated using one embodiment of the inventive method; and

(9) FIG. 8 shows the probability that a given point in the subsurface contains oil (rather than water) as estimated using one embodiment of the present inventive method.

(10) The data shown in these images and used to generate the test examples shown herein are from Avseth, et al. 2010. The invention will be described in connection with example embodiments. To the extent that the following description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative only, and is not to be construed as limiting the scope of the invention. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the scope of the invention, as defined by the appended claims.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

(11) A first embodiment of the invention will now be described. With reference to FIG. 1, this embodiment involves procedures to estimate the rock and fluid properties of a subsurface reservoir from geophysical survey data. As illustrated in FIG. 1, geophysical survey data is obtained (step 101). The inventive method may be applied to geophysical survey data of a single type, for example to a seismic survey, a 2D case example of which is shown in FIG. 3, or of multiple types, for example combining at least two of: a seismic survey, a CSEM survey, an MT survey, a gravity survey, and a passive seismic survey, where the multiple survey data types sense properties of overlapping subsurface regions. Persons of ordinary skill in the art with the benefit of this disclosure will recognize that the present inventive method may be applied to these and other types of geophysical survey data singly or in combination.

(12) Step 102 is selecting an earth model which gives the subsurface properties of interest throughout the subsurface. Earth models that contain a single value for rock or fluid properties and those that provide information about the probability distributions of rock and fluid properties can both be used with the inventive method. Earth models of both fine and coarse gridding can also be used as well as models with adaptive gridding that adjusts grid size, for example, in response to geological complexity or scale. In one preferred embodiment of the invention, a novel earth model is used. This novel earth model is based on inventor recognition that within a reservoir-scale subsurface region, the earth can be approximated as consisting of a small number of rock types and that the individual layers or bodies of each rock type commonly cannot be separately distinguished in geophysical survey data. Thus the earth property sensed by each sample of geophysical survey data is an effective property reflecting an averaging (though not strictly an arithmetic average) of rock properties within a volume. Instead of specifying actual rock properties on a grid fine enough to represent the separate layers of each rock type or specifying effective properties on a coarse grid corresponding to geophysical survey resolution, the novel model specifies the properties of all rock types that contribute to each geophysical survey data point and the volume fraction of each rock type.

(13) Step 103 is determining a rock physics model relating the rock and fluid properties contained in the earth model to the geophysical properties sensed by the survey data. Various types of models could be employed in the inventive method to characterize the relationship between an effective material property of the sample and the volume fractions of the material phases that compose the rock. The use of any such model is within the scope of the invention. These models could be based on some simplifying physical approximations, as with theoretical models based on differential effective medium (DEM) theory or self-consistent (SC) theory. Alternatively, the functional forms used in empirical models could be employed. For example, the slope and intercept of a linear curve fit could be determined from laboratory or well-log measurements. For example, FIG. 5A shows a set of well log data containing Vp, Vs and density and FIG. 5B shows the corresponding scatter plots. Given that, lines could then be fit through the clusters of data points in FIG. 5B to determine a relationship among Vp, Vs, and rho. A different line would be fit for each rock type since a different trend clearly exists for each rock type. In one embodiment, horizons would be interpreted in a 2D or 3D seismic image which define boundaries between rock types according to methods familiar to those skilled in the art. Rock physics models would be determined within each of the intervals between horizons at the location where those horizons intersected a well. The rock physics model determined at the well would be applied throughout the 3D volume of the subsurface which was bounded above and below by the same horizons.

(14) Step 104 is determining a noise statistic characterizing the uncertainty in the rock physics model. Referring to FIG. 5B, the data for each trend may be fitted with a best straight line, but the scatter of data about that line is what is referred to as the noise associated with the selected rock physics model. In other words, the trends described by rock physics models have some uncertainty due to natural variability in structure and composition. This uncertainty can be found for a particular model by subtracting the best fit model predictions of bulk properties from measured bulk properties. Statistics describing the uncertainty may then be calculated, such as a variance or standard deviation of the error, an autocorrelation function describing the error correlation as a function of difference in depth between sample points, or a covariance matrix capturing the covariance between errors in model predictions of different bulk properties such compressional and shear wave velocities.

(15) Rock physics models are calibrated and their errors are characterized typically on core plug or well log measurements. Core plugs are samples of rock with dimensions on the order of 3 cm. Well logs sample volumes of rock on the order of 1 m in linear dimension. Geophysical surveys, however, sense rock properties at scales on the order of 50 m for seismic surveys, and on the order of 500 m for resistivity and gravity surveys. Statistics of rock physics model uncertainty derived at core plug or log scale should be compensated for this scale difference before being applied to rock property variability as sensed by geophysical surveys. This compensation should include adjusting both the trend line and the amount of scatter around that line.

(16) The constraints on the inversion provided by the rock physics model may be written as, f(,x,y,z)=n(x,y,z) where f is a deterministic, spatially varying vector function of the rock composition and bulk properties, , and n is the error vector representing the uncertainty in the rock physics models. In one embodiment, this relationship is linearized around a prior local estimate of the properties .sub.0(x,y,z). In this linearized form, the rock physics model constraints can be written as:
A(.sub.0)=n
where A is the matrix of derivatives of the error with respect to each property in , where the derivatives are evaluated at .sub.0. In one embodiment, n will be considered to be a zero-mean Gaussian random vector (one element per constraint) with covariance matrix P.

(17) Step 105 is selecting a geophysical forward model which relates the geophysical survey data to the sensed geophysical properties of the subsurface. For example, the geophysical forward model for a seismic survey would predict the recorded seismic data traces for each source-and-receiver pair as a function of the elastic properties and density throughout the subsurface. The geophysical forward model for a magnetotelluric survey would predict the measured components of the electrical and magnetic fields at each sensor position as a function of the electrical impedance tensor throughout the subsurface. The geophysical forward model for a CSEM survey would predict the measured components of the electrical and magnetic fields for each pair of source and receiver positions as a function of the electrical impedance tensor throughout the subsurface. The geophysical forward model for a gravity survey would predict the measured variation in the earth's gravitational field at each gravimeter position as a function of the density structure of the subsurface.

(18) With the benefit of this disclosure, persons of ordinary skill in the art will recognize various geophysical forward models that could be used in embodiments of the invention. A particular model for each data type in each application may be selected based on a tradeoff between geophysical forward model accuracy and computational speed.

(19) A novel geophysical forward model is now described which is used in one embodiment of the invention. This model approximates hydrocarbon reservoirs as being composed of two different rock types, a permeable rock which contains the hydrocarbons, such as a sandstone or a permeable limestone, and an impermeable seal rock which prevents the hydrocarbons from escaping from the reservoir, such as a shale or impermeable limestone. Several such layers may be contained in the 10 to 100-meter depth intervals resolvable by seismic data and numerous such layers are commonly contained in the 100 to 1000-meter depth intervals resolvable by resistivity or gravity surveys. Within the depth interval of seismic resolution, the properties of the permeable rock layers are often similar, and the properties of the impermeable rock layers are often similar. The seismic reflection amplitudes as a function of source-receiver offset from such an interbedded system of two distinct rock types is approximately that of a reflection from a single interface between the rock types as described by Zoeppritz equations or Aki-Richards equations or other equations that approximate the Zoeppritz equations, except for a scale factor that has a magnitude and sign that depends on the layer thicknesses within the interval of resolution. There is additional uncertainty in the scale factor due to imperfect knowledge of the seismic wave attenuation between the surface and reflector as well as uncertainty in the source amplitude, receiver sensitivity, coupling between the earth and the source and receivers, and spreading of the signal. Thus, in one embodiment, the geophysical forward model for seismic data traces d in a common depth point gather after moveout correction is,
d(t)=w(t)*[u(t)r((t),t]+v(t)
where * indicates convolution, r are the reflection coefficients from Zoeppritz equations for a reflection from an interface between two half-spaces of the rock types encountered at time t as a function of the incident angle at time t, w is the seismic wavelet, and u is the scale factor for each time point. The contributions to u include geometrical spreading, attenuation, uncertainties related to the data acquisition, and the relative abundance of the two rock types. The noise in this model is v, which includes both random environmental noises and systematic errors due to the approximations made to formulate the model. A novel feature of this model is that it does not predict the seismic data trace based on a specific value for the elastic properties or density of the subsurface at any time or depth point. Instead, it predicts characteristics of the seismic data trace based on the elastic properties of all the rock types in the vicinity of each time or depth point. Data that were thought to conform to this model would be prepared for inversion by performing a deconvolution to remove the effect of w, and a scaling to remove the effect of u, leaving the processed reflection amplitudes d(t)=r(,t)+m(t). An example of such noisy data d(t) is shown in FIG. 6 as the curve with square markers. The other curve with disc markers is the model prediction. The noise after this processing and deconvolution is m, which in one embodiment is assumed to be zero-mean Gaussian with covariance matrix R, and the shale-over-sand reflection coefficients r(,t) may now be analyzed to determine the shale and sand properties effective at the subsurface location interrogated by this CDP gather of data at time t. FIGS. 4A and 4B provide an example data set of seismic reflection at a sequence of CDP points with the same horizontal locations. The variations of reflection amplitude as function of offset are shown across a range of two-way travel time. The parameter vector contains the typical properties of the two rock types within the resolution interval. This method of specifying the subsurface properties is believed to be a novel feature of this disclosure. The geophysical forward model in one embodiment of the invention relates the survey data to the typical properties of two distinct rock types over the interval of data resolution without specifying which rock type exists (and thus without specifying the properties that exist) at any particular point in the subsurface.

(20) A second geophysical forward model is described which is used in one embodiment of the invention. In this embodiment, the geophysical forward model for seismic data traces in a common depth point gather after moveout correction may be represented as
d(t)=w(t)*r((t),t)+v(t)
where * indicates convolution, r are the reflection coefficients from Zoeppritz equations for a reflection from an interface between two half-spaces of the rock types whose interface occurs at time t as a function of the incident angle, and w is the seismic wavelet. The noise in this model is v, which includes both random environmental noises and systematic errors due to the approximations made to formulate the model. Data that were thought to conform to this model would be prepared for inversion by performing a deconvolution to remove the effect of w, leaving
d(t)=r(,t)+m(t)
where the noise after this deconvolution is m, which in one embodiment is assumed to be zero-mean Gaussian with covariance matrix R.

(21) In one embodiment of the invention, the relationship between the processed m reflection amplitudes and the rock properties is linearized around .sub.0 and d.sub.0=r(.sub.0) to give,
dd.sub.0=C(.sub.0)+m
where C is the matrix of derivatives of each reflection amplitude with respect to each rock property.

(22) Returning to FIG. 1, step 106 is determining a statistic describing the noise in the geophysical survey data. The geophysical survey data may be regarded as the sum of a component that is consistent with the geophysical forward model and a noise term. By definition, the noise term is the difference between the geophysical survey data and the geophysical forward model when applied to the correct geophysical properties for the subsurface. The noise term can be defined only with respect to a particular geophysical forward model and includes random or environmental noises as well as any systematic errors due to the approximations in the geophysical forward model. Statistics describing the noise can be calculated, for example, from a gather of geophysical survey data by finding the geophysical properties for the subsurface for which the geophysical forward model most nearly reproduces the gather of survey data, and subtracting the forward model result from the survey data. This difference is the noise in the survey data and statistics of this noise can be calculated. Examples of statistics of the noise include a variance or standard deviation of the noise, an autocorrelation function of the noise, cross-correlation functions between noises at different receiver locations or source-receiver offsets, and covariance matrices between noises at different receiver locations or source-receiver offsets.

(23) Step 107 is estimating a subsurface property or properties, for example in the case of seismic data, the fluid velocity and density or their linear combination. Let the property of interest be =h(). Linearizing around .sub.0 gives
=h(.sub.0)+H(.sub.0),
where H is a row vector with each element the derivative of with respect to each of the properties in . In one embodiment, the estimate of the subsurface property is a minimum mean squared error estimate of from d, call it .sub.est(d):
.sub.est(d)=h(.sub.0)+H[A.sup.TP.sup.1A+C.sup.TR.sup.1C].sup.1C.sup.TR.sup.1(dd.sub.0).

(24) In some cases, the system may not be sufficiently linear to obtain a good estimate of a with one iteration. In that case, in one embodiment, the system may be relinearized around each new estimate of
.sub.est=.sub.0+[A.sup.TP.sup.1A+C.sup.TR.sup.1C].sup.1C.sup.TR.sup.1(dd.sub.0)
until convergence is achieved. FIG. 6 shows one example of the measured data and the corresponding model prediction (smooth curve) obtained from a converged solution, at a specific CDP point.

(25) Step 108 is using the estimate of the subsurface property to assess the commercial potential of a subsurface reservoir. The commercial potential of a subsurface reservoir is determined by a variety of factors including the type or types of fluid in the pore space of the reservoir rock, the porosity of the reservoir rock, and the permeability of the reservoir rock. The subsurface property that is estimated may be one of these properties or may be another property that correlates with one of these properties or may be another factor that influences the commercial potential of the reservoir.

(26) A second embodiment of the invention will now be described. With reference to FIG. 2, this embodiment involves procedures to estimate the rock and fluid properties of a subsurface reservoir from geophysical survey data. As illustrated in FIG. 2, geophysical survey data is obtained (step 201). This step is the same as step 101 described previously. Steps 202-206 are also the same as steps 102-106, respectively.

(27) Step 207 selects one or more subsurface properties to be estimated instead of the full set of model parameters . These selected subsurface properties are typically quantities relevant to particular aspects of hydrocarbon commercial potential. For instance, this could be a certain direct hydrocarbon indicator (DHI) x as a function of the model parameters i.e. x=f(). The subsequent steps then estimate these properties from a set of reduced attributes that preserve all the information in the raw survey data about x.

(28) Step 209 is to apply a transformation to the geophysical data set d changing it into a much smaller data set of attributes, z. This transformation is selected to preserve substantially all of the information about the subsurface property or properties of interest that was present in d. In one embodiment, this transformation is linear, so that it is described by a matrix S which is an MN matrix with M<<N, and z=Sd. To determine S, an information metric is defined (step 208) that will measure the extent to which a given S produces a z that captures the information of interest. In one embodiment, this information metric is J(x|z=Sd), where
J(x|z=Sd)=tr{F[(SC).sup.t(SRS.sup.t).sup.1(SC)].sup.1F.sup.t}
where the rows of F consist of derivatives of f(). The matrices C, R are defined as previously. In one embodiment, we seek S to maximize the information metric,

(29) $S=\underset{S{R}^{MN}}{\arg \max }J\left(x|z=\mathrm{Sd}\right)$
In one embodiment, M is chosen to be the smallest integer value such that the value of J(x|z=Sy) is substantially the same as when S is equal to the identity matrix. In another embodiment, M is the largest integer that still produces a data set z of manageable size. The resulting data are of size M and much smaller than the original size N, however, they preserve substantially the same amount of information regarding x.

(30) Step 210 estimates the subsurface properties of interest, x, directly from the attributes z, instead of from the data d. In other words, the same optimization process described in connection with step 107 may be used in step 210, except that the reduced attribute data volume is used instead of the full data volume. Steps 209-210 may be iterated in nonlinear cases until the solution converges. As an example, FIG. 7 shows the scatter plot of the estimated fluid properties, including fluid velocity and density, throughout the spatial region defined in FIG. 4B.

(31) Step 211 uses the resulting estimate of x to assess the commercial hydrocarbon potential of the reservoir as described in step 108. For instance, the resulting estimates given in FIG. 7 allow fluid type classification and the subsequent calculation of the probability of the fluid being oil or brine at each spatial location, which is shown in the color heat map in FIG. 8 (presented here as a black-and-white reproduction due to patent law restrictions on the use of color). The circled region showing significantly increased oil probability over a continuous area is matched consistently, through the horizontal line, with the well log oil facie profile.

(32) Thus, a key advantage of the FIG. 2 embodiment of the present invention is that once a subset of parameters or DHI's has been identified (instead of using a fully parametric earth model that is fully populated over the computational grid), then all the input information that is needed is those portions of (or transformation of) the original data that contain the information relevant to these identified parameters or DHI's. That is, the reduced (in volume of data) attributes derived from the original data are specific to the chosen subset of parameters or the DHI. For example, the attributes that would be relevant for estimating rock porosity would be different than for estimating fluid velocity.

(33) Among the advantages of this approach is that it enables inverting directly for a DHI, which is not necessarily a parameter in the earth model to be inverted, or in current practice has to be interpreted after the inversion, e.g. computing Poisson's ratio from inverted Vp and Vs, which comes with the high cost of doing a full data-set inversion. Directly computing the DHI from data without multiple intermediate steps results in a smaller parameter model space, as well as a smaller attribute space to work with. Moreover the approach of FIG. 2 potentially leads to better inversion of the subset of interested parameters due to the fact that all the relevant information is preserved, which is not necessarily the case in multi-stage process. This may include improved convergence because the problem may be better conditioned due to the need to estimate only a smaller parameter set over a potentially patchy distribution. Where the other parameters are not of interest, this can be the preferred embodiment of the invention.

(34) The foregoing application is directed to particular embodiments of the present invention for the purpose of illustrating it. It will be apparent, however, to one skilled in the art, that many modifications and variations to the embodiments described herein are possible. With the benefit of the disclosure herein, persons of ordinary skill in the art will recognize that the inventive method may be applied to relate other effective properties besides seismic properties, permeability, and conductivity to structure and composition. For example, density, magnetic permeability, thermal conductivity, and nuclear magnetic resonance relaxation times can all be related using the inventive method to the structure and composition of rocks. All such modifications and variations are intended to be within the scope of the present invention, as defined in the appended claims. Persons skilled in the art will also readily recognize that in preferred embodiments of the invention, at least some of the steps in the present inventive method, for example steps 107 and 210, are performed on a computer, i.e. the invention is computer implemented. In such cases, the resulting model parameters may either be downloaded or saved to computer memory.

REFERENCES

(35) Avseth et al., Quantitative Seismic Interpretation: Applying Rock Physics Tools to Reduce Interpretation Risk, Cambridge University Press, particularly pages 336-338 (2010).