JOINT INVERSION WITH UNKNOWN LITHOLOGY

CROSS-REFERENCE TO RELATED APPLICATION

[0001] This application claims the benefit of U.S. Provisional Patent Application 61/492,624, filed June 2, 2011, entitled JOINT INVERSION WITH UNKNOWN LITHOLOGY, the entirety of which is incorporated by reference herein.

FIELD OF THE INVENTION

[0002] This invention relates generally to the field of geophysical prospecting and, more particularly, to joint inversion of two or more different types of geophysical data to infer physical property models of the subsurface. Specifically, the invention is a method for inferring lithology as well as geologic properties in the joint inversion.

BACKGROUND OF THE INVENTION

[0003] This invention pertains to using geophysical data in a joint inversion to infer geological properties of the subsurface. During an inversion, the aim is to minimize the difference between the measured data and the data predicted by the inversion model. In order to perform a predicted data calculation, geophysical parameters such as seismic velocity (or elastic coefficients), or electrical conductivity must be known. When multiple data types (e.g. reflection seismic and electromagnetic data) are inverted simultaneously, it is known as a joint inversion. Geophysical data are likely to include active seismic reflection data; active seismic refraction data; electromagnetic data (either controlled source or magneto-telluric); and/or gravity measurements; however, it may in addition include any other type of data that can be used to infer the properties of subsurface rocks in the region of interest.

[0004] Geophysical properties, such as elastic coefficients, density, and electrical conductivity, can be converted to the geological properties of interest in hydrocarbon exploration (e.g., porosity and fluid type) via rock physics relationships (obtained empirically or theoretically). In this way the different geophysical data types are linked in the joint inversion. These rock physics relationships can be embedded in a joint inversion of geophysical data. They are used to calculate the needed geophysical parameters (elastic coefficients, electrical conductivities, and density) that are necessary for forward calculating the predicted data. Using the forward predicted data, a misfit between the predicted and observed data is computed. The model is then iteratively updated using some optimization scheme to minimize the difference between predicted and measured data. Attorney Docket No. 2011EM148

[0005] In general, to perform joint inversions of this type, one must assume a priori a particular rock physics relationship between the geophysical parameters (for example sonic velocity, shear wave velocity, density or conductivity) that predict the data and the geological parameters (for example porosity or water saturation) of interest. By assuming a rock physics relationship we are assuming a lithology and depositional environment present in the subsurface. A lithologic class is a rock type that is considered to possess unifying rock physics behavior for the purposes of the inversion; e.g., elastics and carbonates might be considered two distinct lithologic classes, each with their own rock physics relationship. However, the lithology in the subsurface of a particular region of interest is often not known beforehand, and further, a single physical volume sampled by the data may contain more than one lithology with an unknown spatial distribution of those lithologies.

[0006] One way to jointly invert multiple geophysical data for geophysical properties is to assume structural coupling (e.g., Haber and Oldenburg, 1997) where anomalies in one of the geophysical properties (e.g., velocity) are required to occur in the same location as anomalies in one or more of the other geophysical properties (e.g., resistivity). The problem with this approach is that it is highly nonlinear, when data of very different resolutions are being inverted. This makes it practically challenging to invert, for example, high frequency seismic data together with low frequency CSEM data.

[0007] A joint inversion for geophysical properties can also be performed by assuming explicit or implicit relationships between the parameters. For example two parameters can be assumed to be correlated (see for example Farquharson et al, 2010). The problem with this approach is that these relationships have to be known beforehand and must be adequate for the subsurface area of interest, or the inversion will fail.

[0008] In order to infer geological properties from an inversion for geophysical properties, a rock physics model can be used to convert the inverted geophysical properties into geological properties. Even though this approach allows one to infer geological properties, it relies on the inverted geophysical properties. The conversion does not rely on the measured data and thus does not allow for feedback between measured data and geological properties.

[0009] Doetsch et al. (2010) perform a joint inversion using structural coupling.

Following the inversion, they analyze the inverted geophysical properties for patterns of similar properties, such as zones that are fast and resistive versus zones that are slow and conductive. These zones of similar geophysical properties are then treated as one model cell Attorney Docket No. 2011EM148

and they invert for average geophysical properties for each zone. In the next step they use the average properties to do an after-the-fact conversion to the average geological properties, using rock-physics relationships. The problem with this approach is that it relies on the data being of similar resolution due to the structural coupling. Furthermore, as described before, there is no feedback between the inferred average geological properties and the measured data. In general applications, therefore, this method may not succeed.

[0010] Another approach to joint inversion uses statistical methods. A lithology in these methods is simply defined as a class of rocks that can be assigned a probability density function ("pdf ') of continuous parameters (e.g. seismic p-wave velocity, or porosity) - no explicit rock physics equations are necessary. The use of the statistical method is, for example, demonstrated by Guillen et al (2004) who use gravity and magnetic data to invert for lithology of the subsurface; Buland et al. (2008) use a similar technique to invert for seismic reflection data. Unfortunately, this approach assumes that the pdf for each lithology is known beforehand, which is rarely the case for most exploration settings.

[0011] In another approach for joint inversion, the geologic environment is assumed known and the corresponding rock physics model is applied (see Abubakar et al, 2010; Jing et al, 2010; Hoversten, 2010). For example, if the lithology is assumed to be clastic, a clastic rock physics model is used to relate the velocity, density and resistivity to the rock properties (e.g., Vshale, porosity, Water saturation). This has the advantage that the geological parameters are inverted for directly, i.e., this approach allows for feedback between geologic parameters and the measured data. But unless the lithology is known beforehand, the assumption of a specific rock physics model can strongly bias the inverted result. For example, in the event the lithology is in fact a volcanic, such as basalt, and not clastic, the estimates of velocity and density will be incorrect and the characterization of rock in terms of Vshale will not make sense. That is, unless the lithology is known beforehand, the assumption of a specific rock physics model can strongly bias the inverted result, leading to incorrect results.

[0012] DiCaprio et al. (2010) present an approach that allows jointly inverting geophysical data for subsurface properties in cases where the lithology class is not known beforehand. Instead of assigning a rock type a priori, their invention prescribes using the lithology classes as a discrete inversion parameter to be found during the inversion. At each step in the inversion, the appropriate rock physics relationship is used on the resolution cells depending on what lithology they are currently assigned. The lithology parameter is allowed Attorney Docket No. 2011EM148

to vary both as the inversion evolves and as a function of space (allowing for mixed lithologies in a single physical volume). A drawback of this approach is that specialized optimization schemes must be used, and the model space is greatly expanded. This may not be computationally practical for all applications.

[0013] The invention presented here is an alternative approach to DiCaprio et al.

(2010). It also allows for the inversion of geological properties in cases where the lithology is unknown. There is no restriction to use data of similar resolution. Instead of one inversion with additional parameters to be inverted for (DiCaprio et al, 2010), it uses different data coupling strategies in different stages of the inversion to arrive at a model of geological properties and lithologies.

SUMMARY OF THE INVENTION

[0014] In one embodiment, the invention is a computer-implemented method for joint inversion of two or more sets of geophysical data of different types, measured in surveys of a subsurface region, to obtain a model of at least one geological property for the subsurface region, said method comprising:

using a computer to jointly invert the sets of geophysical data, using structural coupling between the different data types, to obtain models of geophysical properties corresponding to the sets of geophysical data, said structural coupling being determined from the geophysical data or from a priori knowledge of the subsurface region;

partitioning the subsurface region into sub-regions based on similar combinations of geophysical parameters in the geophysical property models that correlate to particular lithologies, thereby defining lithology sub-regions;

determining mathematical rock physics relationships appropriate for each lithology, said rock physics relationships relating geological properties to geophysical properties; and using a computer to jointly invert the sets of geophysical data, using the rock physics relationships according to lithology sub-region, to obtain a model of one or more of the geological properties for the subsurface region.

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:

Fig. 1 is a flowchart showing basic steps in one embodiment of the present inventive method; Attorney Docket No. 2011EM148

Fig. 2 is a schematic diagram illustrating how the method of Fig. 1 works for a case of mixed lithologies (clastic and basalt);

Fig. 3 is a table of lithologies commonly encountered in hydrocarbon exploration;

Fig. 4 is a flowchart showing basic steps in one embodiment of the present inventive method; Fig. 5 is a flowchart showing basic steps in an embodiment of the present inventive method for determining average geological properties for lithologic units;

Fig. 6 is a flowchart showing basic steps in the embodiment of Fig. 5 modified with an additional step of geometric coupling;

Fig. 7 is a flowchart showing basic steps in the embodiment of Fig. 6 modified to assess uncertainties of the locations of the lithologic boundaries; and

Fig. 8 is a flowchart showing how the method of Fig. 4 may be modified to add an additional intermediate step of geometric coupling.

[0016] The invention will be described in connection with example embodiments.

However, to the extent that the following detailed 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. Persons skilled in the technical field will readily recognize that in practical applications of the present inventive method, it must be performed on a computer, typically a suitably programmed digital computer.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

[0017] This invention is a method for a rock physics based joint inversion of geophysical data for subsurface properties in cases where the lithology class - also referred to as rock type or facies - is not known beforehand. This is done in a workflow (figures 1 and 2), in which the coupling strategies may change at each step. Figure 1 is a flowchart showing basic steps in one embodiment of the present inventive method.

[0018] In step 1, a structural coupling method, for example Haber and Oldenburg

(1997), is used to jointly invert the geophysical data of at least two data types. In structural coupling the sole assumption is that the inverted geophysical properties, for example sonic velocity and conductivity, have structural similarity, i.e. boundaries and changes in the properties are co-located ("structural constraints"). The result of a joint inversion using Attorney Docket No. 2011EM148

structural coupling is geophysical properties, e.g. sonic velocity, density, conductivity. The term structural coupling is used herein to mean include all forms of coupling that would make the inversion joint and that could be used without knowledge of lithology. It also includes additional constraints that might be interpreted from the data as well as come from previous knowledge, such as known boundaries that can be enforced in the inverted geophysical parameters.

[0019] In step 2, the lithology is determined by identifying zones that have similar geophysical properties, for example zones that have a large sonic velocity and a small conductivity as opposed to small sonic velocity and high conductivity. Once these zones are identified, the lithology of each zone ("lithology volume") can be determined and an appropriate rock physics relation is assigned.

[0020] In step 3, a joint inversion using rock physics coupling is performed. In this approach the rock physics relations are embedded in the joint inversion of the different data types for a common set of geological properties or parameters (see, for example, Jing et al. 2010).

[0021] Figure 2 is a schematic diagram illustrating how the method of Fig. 1 works for a case of mixed lithologies (clastic and basalt).

[0022] Note, that the first two steps of the workflow are similar to Doetsch et al.

(2010) with the exception that the present invention uses the resulting lithology zones to assign a rock physics model, while Doetsch et al. invert for average geophysical properties of each zone. Doetsch et al. (2010) use the average properties to do an after-the-fact conversion to the average geological properties, while the current invention performs a rigorous joint inversion using rock physics to couple the different geophysical data. Another key difference is that Doetsch et al (2010) determine only average geological properties for each zone, while the present joint inversion method yields geological properties at the employed grid spacing.

[0023] In frontier hydrocarbon exploration the broad lithology classes of the rocks are often not known; in this case, possible lithologies for the inversion may be clastic, carbonate, salt, and basalt. The particular possible lithologies chosen will be based on prior geological knowledge of the area. The data sets collected must be able to distinguish among the chosen lithologies. Figure 3 demonstrates how a combination of different data types can allow distinguishing among several possible lithologies in step 2 of Fig. 1. Attorney Docket No. 2011EM148

[0024] Figure 4 is a flowchart showing, in more detail, steps in a particular embodiment of the present inventive method.

[0025] In step 41, geophysical data are collected. The data set might include collecting two or more of active seismic reflection and/or active seismic refraction, controlled source electromagnetic, MT, and gravity data over a region of interest.

[0026] Because structural coupling becomes highly nonlinear if data types of very different resolution are combined, a low pass filter may be applied to the high frequency data types (step 42), for example seismic reflection data. This decreases the resolution of these data types and gives more stable structural coupling. Furthermore, the filtering may mitigate the local minima problem associated with the inversion of seismic reflection data.

[0027] In step 43, the subsurface is discretized into inversion grid cells. Because only low frequencies are used in this step, it is possible to use a coarser grid than if all frequencies would be considered. This gives an additional advantage of reduced computation time for this step. Furthermore, an initial guess for the geophysical properties in each grid cell is created using any prior information available.

[0028] In step 44, a cross gradient approach, for example, may be used to couple the different data in a joint inversion. The so called cross gradient constraint, which is a form of structural coupling, is added to the objective function (e.g. Gallardo and Meju (2003, 2004)). This constraint tries to co-locate changes in each of the model parameters independent of their magnitude. The result of the inversion will be geophysical property volumes 45.

[0029] In step 46, clusters/zones with similar combinations of geophysical parameters, such as fast sonic velocity, high density, and high resistivity are identified. Some type of cluster analysis (e.g. MacQueen (1967), Kaufman and Rousseeuw (2005), Doetsch et al. (2010)) or pattern recognition (e.g. Specht (1990), Sarle (1994)) may be used to identify these zones with similar geophysical properties.

[0030] Next in step 46, a finite set of lithologic classes (e.g. basalt, salt, clastic, and carbonate) is picked. Each class has a corresponding set of equations that relate the geophysical properties (e.g., velocity, density, conductivity) to the geological properties (e.g., porosity, lithology, fluid type).

[0031] Finally in step 46, each zone of similar geophysical properties is assigned a lithology - a basalt, for example; rock physics equations associated with each lithology are chosen. This can be done as an interpretation step, using a data base (e.g. Fig. 3), or using Attorney Docket No. 2011EM148

some more rigorous mathematical approach that assesses which rock physics relation predicts the geophysical properties best.

[0032] Once the rock physics relations are assigned, the subsurface is discretized in step 47 into inversion grid cells at the desired scale length appropriate for an inversion utilizing all available frequencies. An initial guess for the geological properties is constructed also in step 47 using the geophysical properties inverted in step 44 together with the rock physics equations assigned in step 46.

[0033] In step 48, a joint inversion using rock physics is performed (see, for example,

Jing et al, 2010). Each zone gets assigned the rock physics relationship of the lithology determined in step 46. The data sets are now coupled using rock physics relationship and for each zone the inversion is performed for a common set of geological parameters. The inversion may be performed using the full resolution of all data and the desired grid spacing.

[0034] During the joint inversion, the aim is to minimize the difference between the measured data and the data predicted by the inversion model. In order to perform a predicted data calculation, geophysical parameters such as seismic velocity (or elastic coefficients), and electrical conductivity must be known. Using the geologic parameters from the current geological properties model as input quantities, the rock physics models assigned to each zone are used to calculate the needed geophysical parameters such as elastic coefficients, electrical conductivities, and densities that are necessary for forward calculating the predicted data. Using the forward predicted data, a misfit between the predicted and observed data is computed. The model is iteratively updated using any one of a number of optimization schemes. In determining the update, damping terms and lithologic/rock physics constraints may be included in the objective function. This total objective function (misfit and additional terms) is to be minimized in the inversion, resulting in best estimates of the geological properties. The constraints are non-deterministic parts of the rock physics model that help to restrict solutions to geologically realistic combinations of parameters. For example, in a cell with a clastic lithology, porosity and Vclay might be expected to be inversely correlated, so a term could be added to the objective function that penalizes solutions that do not follow this trend.

[0035] When the optimization process has converged on a model, a possible solution

49 to the inverse problem in the geological property space has been found. Attorney Docket No. 2011EM148

[0036] There may be circumstances in which the values of any continuous inversion parameters are not of interest and only lithologic packages and their average geological parameters are to be found. In this case, illustrated by the flowchart of Fig. 5, the joint inversion using structural coupling and the cluster analysis may be performed as described. But instead of performing the inversion on a fine grid, each cluster is treated as one inversion cell at step 57, and the inversion is performed at step 58 for just the average geological properties of each cluster 59. This is similar to the approach by Doetsch et al. (2010) with the main difference that Doetsch et al. (2010) invert for average geophysical properties and then do a conversion to the geological properties, while the present invention assigns a rock physics model to each zone and performs a rigorous inversion of all geophysical data types for the average geological properties, i.e., the present inventive method allows for feedback between data and geological parameters.

[0037] The embodiment of Fig. 5 may be used with an additional step before the inversion of the average geological parameters, as shown in the flowchart of Fig. 6. The purpose of this intermediate step is to provide a better starting model for the iterative joint inversion using rock physics coupling (58, 68) by allowing the cluster boundaries to change to account for the fact that each cell gets assigned only average geophysical properties. In this intermediate step 66, after the clusters are identified in step 46, the average geophysical parameters are inverted for each cluster, and the coupling of the data is accomplished by performing that inversion for a common geometry of all geophysical parameters, for example the thickness of each cluster. In other words, the thickness of each cluster is an additional unknown, in addition to the average geophysical parameters for each cluster, that is inferred in the inversion of step 66. Since the cluster thickness applies to each data type, this provides the coupling that makes the inversion "joint" rather than sequential. This is similar to, and serves the same "jointness" purpose as, the structural coupling in step 44; however, the various types of geometrical coupling would not be known at the outset, i.e. at the stage of step 44. Using a common geometry to couple different data is described in Mataracioglu and Asci (2010), who jointly invert different geophysical data for the location of an ore body. The resulting new geometry 67 is then used in step 68 to invert the average geological parameters 69. In other words, the average geophysical parameters resulting from step 66 are substituted into the appropriate rock physics relationships to yield average geological properties for a starting model for the iterative joint inversion in step 68, and the better defined clusters (locations of the cluster boundaries are better known) that are also inferred in step 66 define the grid cells for the inversion of step 68. Attorney Docket No. 2011EM148

[0038] Once the average geological properties for each cluster are determined using one of the embodiments of the invention such as those illustrated in Figs. 5 and 6, they can be perturbed at step 71 in Fig. 7, and an inversion for the common geometry can be performed at step 72 for different perturbations. This can be used to infer the possible range of thicknesses of lithologic units. This does not provide estimates of uncertainty of the geological properties, as for example obtained by a stochastic inversion (Hoversten, 2010), but instead provides uncertainties related to the geometry of a lithological unit 73.

[0039] The embodiment of the invention illustrated in Fig. 4, with the additional intermediate step outlined earlier (step 66 in Fig. 6), during which the data are coupled by inverting for a common geometry of the identified zones/cluster, is outlined in Fig. 8.

[0040] The foregoing patent 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. All such modifications and variations are intended to be within the scope of the present invention, as defined in the appended claims.

References

Abubakar, A., Gao, G., Habashy, T.,Liu, J., "Joint electromagnetic and seiosmic data inversion algorithm for geophysical applications," American Geophysical Union, Fall Meeting 2010, Abstract #NS43A-02 (2010).

Buland, A., Kolbjornsen, O., Hauge, R., Skjaeveland, O., and Duffaut, K., "Bayesian lithology and fluid prediction from seismic prestack data "Geophysics 73(3), pp. C13-C21 (2008).

DiCaprio, C, Saltzer, R., Mullur, A., "Inverting geophysical data for geological parameters or lithology," U.S. Provisional Patent Serial No. 61/368,035.

Doetsch J., Linde, N., Coscia, I., Greenhalgh, S.A., Green, A.G., "Zonation for 3D aquifer characterization based on joint inversions of multi method crosshole geophysical data," Geophysics 75(6), 53-64 (2010).

Farquharson, C.G., Leviere, P.G., Hurich, C.A., "Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exportation," AGU Fall Meeting 2010, NS43A-05 (2010). Attorney Docket No. 2011EM148

Gallardo, L.A and M.A. Meju, "Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data," Geophysical Research Letters 30, 1658 (2003).

Gallardo, L.A and M.A. Meju, "Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradient constraints," Journal of Geophysical Research 109, B0331 1 (2004).

Guillen, A., Courrioux, G., Calcagno, P., Lane, R., Lees, T., Mclnerney, P., "Constrained gravity 3D litho-inversion applied to Broken Hill," ASEG 17 ^th Geophsisical Conference and Exhibition, Sydney (2004).

Haber, E and Oldenburg D, "Joint inversion: a structural approach," Inverse Problems 13, 63- 77 (1997).

Hoversten G.M., "Stochastic inversion ofgeophysical data for estimatingearth model parameters," U.S. Patent Application Publication US2010-0185422 (2010).

Jing, C, Carazzone, J.J., Rumpfhuber, E.M., Saltzer, R.L., Dickens, T.A., Mullur, A.A., "Hydrocarbon detection with passive seismic data," PCT Patent Publication No. WO2010- 080366 (2010).

MacQueen, J., "Some Methods for classification and Analysis of Multivariate Observations," Proceedings of 5th Berkeley Symposium on Mathematical Statistics and Probability, Berkeley, University of California Press, 1 :281-297 (1967).

Mataracioglu M., and M. Asci, "Effectivness of joint inversion for mine survey," 2010 AGU Fall Meeting, NS41B-1509 (2010).

Sarle, W.S., "Neural networks and statistical models," Proceedings of the Nineteenth Annual SAS Users Group International Conference, April, 1-13 (1994).

Specht, D.F., "Probabilistic neural networks," Neural Networks 3, 109-118 (1990).