Domain freezing in joint inversion

Abstract

Method for estimating geological properties in a subsurface region using multiple types of geophysical data (21). An initial physical properties model 22 is constructed. Some parameters in the model are frozen (23) and optionally portions of the model wave number and spatial domains (24) and the data frequency and data time domains (25), are also frozen. Then, a joint inversion (26) of the multiple data types is performed to calculate an update to the model only for the portions that are not frozen. The converged model (27) for this inversion is used as a new starting model, and the process is repeated (28), possibly several times, unfreezing more parameters and data each time until the desired spatial and parameter resolution (29) has been achieved.

Description

FIELD OF THE INVENTION

This invention relates generally to the field of geophysical prospecting and, more particularly, to processing of geophysical data. Specifically, the invention is a method for performing joint inversion of two or more different geophysical data types.

BACKGROUND OF THE INVENTION

This invention pertains to joint inversion of remote geophysical data to infer geological properties of the subsurface. Remote geophysical data are likely to include active seismic reflection data; electromagnetic data (either controlled source (“CSEM”) or magneto-telluric (“MT”); and/or gravity measurements; however, any type of data that can be used to remotely infer the properties of subsurface rocks in the region of interest may be included. When multiple data types (e.g. reflection seismic and electromagnetic data) are inverted simultaneously, it is known as a joint inversion. During inversion the aim is to minimize the difference between the measured data and the data predicted by the inversion model. By combining multiple different types of geophysical data in a joint inversion, one often aims to invert for a model with multiple different types of model parameters (e.g. porosity and fluid type) rather than just a single parameter (e.g. p-wave impedance).

Due to the large number of model parameters and the often large computational cost of the forward calculation (“synthesizing” the data from a test model), one is often limited to linearized, local optimization techniques for inversion. These involve starting at an initial model and updating it by moving along a path in model-parameter space that decreases the misfit between measured data and synthesized data (known as the objective function). Geophysical inversion in general, and joint inversion in particular, often has a highly non-linear objective function which can result in poor convergence properties due to the solution becoming stuck in a local minimum of the objective function.

Lack of convergence due to strong non-linearity of the inversion problem often arises in geophysical inversion due to the nature of the seismic reflection data—specifically, the relative lack of low-frequency content in the data. This problem can be mitigated to some degree by first inverting a low-pass filtered version of the data to find a long-spatial-wavelength model. Using this model as a starting model for subsequent inversions of higher-frequency data can serve to stabilize the inversion process. The technique of first inverting low-frequency portion of the seismic reflection data can be combined with inverting only the earliest portions of the recorded data first, i.e. the earliest arrivals at the detectors. By limiting the time window during the inversion, the more complicated deeper reflections, which are overprinted by multiples, can be excluded to obtain a good shallow model.

Bunks et al. describe a multiscale approach to full waveform seismic inversion. (Bunks, C., Saleck, F. M., Zaleski, S., and Chavent, G., “Multiscale seismic waveform inversion,” Geophysics 50, 5, pp 1457-1473 (1995)) They propose to low-pass filter the seismic data and increase the model grid size in order to avoid many of the local minima normally encountered when inverting full waveform reflection seismic data. At each step, they add more frequencies to the data and reduce the grid size to realize the full resolution available in the data set. This method, however, does not describe how to stabilize a joint inversion of multiple data and parameter types.

Hu et al (2009) perform a joint inversion of electromagnetic and seismic data. (Hu, W., Abubakar, A., and Habashy, T. M., “Joint electromagnetic and seismic inversion using structural constraints,” Geophysics 74, 6, pp R99-R109 (2009)) In order to prevent high-frequency data from dominating the inversion and thus becoming trapped in local minima, they apply a weight to the data such that lower frequency portions of the data are emphasized. The data weighting does not change during the course of the inversion. This technique does not allow one to increase the influence of higher frequency data or to alter which parameters the inversion is solving for as the solution approaches the global minimum.

SUMMARY OF THE INVENTION

In one of its aspects, the invention, with reference to the flowchart of FIG. 2, is a method for estimating geological properties in a subsurface region using multiple types of geophysical data (21). An initial physical properties model 22 is constructed. Some parameters in the model are frozen (23) and optionally portions of the model wave number and spatial domains (24) and the data frequency and data time domains (25), are also frozen. Then, a joint inversion (26) of the multiple data types is performed to calculate an update to the model only for the portions that are not frozen. The converged model (27) for this inversion is used as a new starting model, and the process is repeated (28), possibly several times, unfreezing more parameters and data each time until the desired spatial and parameter resolution (29) has been achieved. In all practical applications, joint inversion is carried out on a suitably programmed digital computer.

In one embodiment, the invention is a method for estimating a physical properties model in a subsurface region using data comprising two or more geophysical data types, said model comprising numbers in a model parameter domain representing one or more physical properties, said method comprising jointly inverting (using a computer) the data in a plurality of sequential phases, wherein in each phase until a last phase, only a portion of the data is inverted in order to infer a subset of the model's parameter domain.

BRIEF DESCRIPTION OF THE DRAWINGS

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 illustrates the relative frequency content of data that may be used in a joint geophysical inversion, in this case reflection seismic data and CSEM data;

FIG. 2 is a flow diagram showing how the general domain freezing process of the present invention may fit into a typical joint inversion scheme;

FIG. 3 is a flow diagram showing basic steps in one embodiment of the present inventive method; and

FIG. 4 illustrates a test example of the present invention.

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

The present inventive method applies specifically to joint inversion. In joint inversion, the different data types often encode information about the subsurface at different wavelengths as illustrated in FIG. 1, which shows the relative frequency content of data that may be used in a joint geophysical inversion of seismic and electromagnetic data. The solid line is the frequency content of reflection seismic data and the dashed line is the frequency content of CSEM (controlled-source electromagnetic) data. This difference in frequency content can be used effectively when attempting to stabilize the inversion by using low-frequency data during the initial stages of the inversion. Also, the present inventive method goes further than simply using different parts of the data during different stages of the inversion. Particular model parameters are fixed during the inversion to improve inversion convergence by using knowledge of the physics of the model parameters and their likely effect on the non-linearity of the inverse problem.

The invention allows for robust joint inversion of geophysical data via knowledge of the frequency content of the different types of data, and the corresponding wavelengths of the different model parameters. The invention allows for a more robust inversion process by approaching the true model at longer wavelengths first before solving for model parameters expected to vary at short wavelengths.

This is accomplished by “freezing” various data and model domains of the inversion problem. In this document, freezing is defined as fixing, damping, down-weighting, or removing particular parts of the objective function pertaining to data or model parameters that might contribute to poor convergence properties during an inversion. As the solution to the inversion problem evolves towards the global minimum of the objective function, the frozen domains are gradually unfrozen; that is, they contribute more strongly to the inversion process.

By first inverting for the longer wavelength portions of the model domain, the inversion becomes more robust and less likely to become stuck in a local minimum. Once a long wavelength estimate of the model has been determined, it is effectively used as a starting model for a subsequent phase of the inversion. This next phase of the inversion (or set of iterations) solves for finer wavelength parameters in the model. The processes can continue, adding higher frequencies and solving for smaller wavelength features until the highest frequency data available has contributed to the inversion solution. At each phase, the previous solution is used as a new starting model. This is the un-freezing process.

The primary freezing domain is the model parameter domain. In many practical applications, the inversion model domain contains more than one physical parameter. For example, a geophysical joint inversion may invert for both conductivity and seismic velocity, or perhaps porosity and water saturation. If one or more of the model parameters are believed to vary on a length scale smaller than the other parameters (i.e., they vary rapidly), one can use that fact to avoid local minima by freezing the rapidly varying parameters for early stages of the inversion process.

Four other example domains are presented here as candidates to be frozen to improve joint inversion convergence. This list is not intended to be exhaustive, and further, domains may be used alone or in combination. The particular problem of interest will determine the domain(s) used, as well as the order of domain freezing in the inversion.

Data Frequency Domain: In a joint inversion each data set is likely to have a different frequency content. For example, electromagnetic data, such as controlled source electromagnetic (CSEM) or magnetotelluric (MT) data has a lower frequency content than seismic reflection data. The high frequency nature of the seismic reflection data (and specifically the lack of intermediate and low frequencies) often cause the inversion to suffer from local minima problems. By low-pass filtering the higher frequency data sets so that all data contribute to model resolution at similar wavelengths, we can avoid local minima of the inversion objective function. After solving for a long wavelength model, we can gradually thaw (add back in) higher frequency portions of the filtered data sets.

Data Time Domain: Multiples and converted waves in the inversion of seismic reflection data make the problem very nonlinear. To avoid these late-arriving phases in the initial steps of the inversion, one can limit the time window to the early portion of the seismograms. This restricts the inversion to the shallow part of the model and it is better behaved due to the lack of multiples and reduced parameter space. Having a more accurate solution for the shallow part of the model will aide the inversion of the deeper portions of the model, which correspond to times in the data that include multiples and converted waves from the shallow part of the model.

Wave Number Domain: Local minima in the geophysical joint inversion processes are directly related to short-wavelength features of the model being mislocated in the physical model space. Because of this, the freezing of short-wavelength features in the model domain can improve inversion stability. This can be achieved by increasing the grid size of the physical model, i.e. completely disallowing small-wavelength features in the inversion model during early steps. This can also be achieved by using a spatial smoothing constraint on the inversion that penalizes small-wavelength fluctuations in the model.

Spatial Domain: Spatial masks can be applied to the model domain. For example, one could freeze deeper portions of the model during early stages of the inversion. After finding a solution for the shallow portion of the model, the deeper portions are then solved. Alternatively, regions of the model that vary rapidly (for example, near faults or major lithologic boundaries) may be frozen until large-scale structure is obtained.

In general, the invention may be applied successfully in any generic inversion workflow, as illustrated in the flowchart of FIG. 2. Further, different benefits may be obtained using different solvers. For example, while these techniques mitigate local minima in gradient-based solvers, they also can drastically reduce computational expense in stochastic or genetic algorithms by substantially reducing the scope of model parameter space searched in each phase of the inversion.

FIG. 3 illustrates an embodiment of the invention that jointly inverts multiple geophysical data types in multiple phases by domain freezing. In this example, it is assumed that rock physics relationships are used to link the different data types. However, the present invention is not concerned with the type of coupling between data types that is used to make the inversion joint.

In the phase labeled (a) in FIG. 3, geophysical data are collected or otherwise obtained. The data might include at least two of active seismic reflection and/or refraction, controlled source electromagnetic, MT, and gravity data, each acquired over the same region of interest.

In phase (b), only the low frequency information of seismic data 301, as obtained by low pass filtering or any other mathematical transformation that eliminates high frequencies, is inverted together with naturally low frequency data 302, such as CSEM, MT or gravity (frequency domain freezing). Because only low frequencies are used, the inversion is performed on a coarse grid 303 (wavenumber domain freezing). Furthermore, any of several other types of domain freezing 304 can be used in phase (b). Simplified rock physics relationships—possibly implemented by fixing or strongly damping one or more inversion partameters—can be used in this phase (model parameter domain freezing). For example, if water saturation were an inversion parameter, one could fix the water saturation in this phase because large scale fluctuations of water saturation are not expected. Further reduction of free model parameters that are inverted for can be achieved by fixing or strongly damping portions of the model space (spatial domain freezing), or the data time domain by applying a taper to the seismic data.

Phase (b) can consist of multiple sub-phases, for example starting by inverting only the early arrivals and in the following phases include information arriving later in time. (In the claims appended hereto, the prefix “sub” is omitted and each one of multiple sub-phases is referred to as a “phase.” Each “phase” is defined by completion of a joint inversion process, which means in the case of iterative inversion, satisfying a convergence criterion.) The final result of phase (b) is a coarsely sampled model space containing information only at low wavenumbers.

In phase (c), the frequency content of the seismic data is increased 305 and the grid size decreased 306 over several stages 307 until frequency content is used together with the desired grid in the final stage 310. The model obtained in each previous sub-phase is used as initial model 308. Furthermore the model may also be used to add an additional damping constraint for regularization. The low frequency data are still inverted jointly with the seismic data, but their contribution to the objective function is decreased 309 as the grid size decreases (data frequency domain freezing). This weighting scheme enables the low frequency data to inherit the resolution of the high frequency data while stabilizing the inversion. The damping at each sub-phase 308 can be increased to ensure that only information at smaller scale length is changed. Furthermore, the complexity of the link between the data types is increased. For example by increasing the number of rock physics parameters to invert for or by decreasing the damping of portions of the rock physics domain. Again, spatial masks, which can change as a function of iteration, can be applied 311 to the model space to reduce the model parameters that are inverted for in each sub-phase.

The workflow is stopped at phase (d) once no more information can be added and the highest resolution is achieved. Reservoir properties 312 (porosity, lithology, fluid type) are the result of the inversion. Other products may be Vp, Vs, density and resistivity, all of which inherited the high spatial resolution of the seismic data.

For hydrocarbon detection, an additional phase can be added at the end. In this phase all parameters except the water saturation are fixed. A very strong damping may be used to the model resulting at 312, and the weights of the data containing information about resistivity (CSEM, MT) and thus water saturation are increased strongly compared to the seismic data. This utilizes the information about water saturation contained in these data types and only allows changes to the model if required by these data. The strong damping prevents a potential loss of resolution due to the large weight of the low frequency data (like CSEM and MT).

For hydrocarbon exploration, a re-gridding of the water saturation can be performed. All layers identified as sand are gridded into one layer. Alternatively a constraint can be added allowing each sand layer to have only one value for the water saturation (thus the same for each layer inside the sand).

EXAMPLES

FIG. 4 illustrates an embodiment of the invention such as that illustrated in FIG. 3 applied to an image problem of a shallow body that has fast seismic velocities and large resistivity. The data used are seismic reflection (amplitude and travel times), CSEM, MT, and gravity. The data are coupled using a rock physics model described by porosity, shale volume fraction (“Vshale”), and water saturation. On the right side of the drawing, 1D models of the three geological properties are shown at the completion of each of the three phases (a)-(c) from FIG. 3.

In the first phase (a), as stated at 401, the seismic data are low-pass filtered with a corner frequency of 10 Hz and the inversion is performed on a coarse grid. The low frequency data are added using a large weight and the water saturation is fixed as indicated by the arrow 501.

In the next phase (b), as stated at 402, more frequencies are added to the seismic data and the inversion is performed on finer grids. The solution from phase (a) is used as a starting model (interpolated on the finer grid) and a damping term is added in the objective function (502). Because the reflections from the shallow, fast, and reflective body dominate the misfit, that portion of the model (503) is kept fixed in subsequent phases. In a first sub-phase the water saturation is fixed, in the second it is allowed to change via the inversion process.

In phase (c), as stated at 403, the inversion is performed using the full resolution of the seismic data and the model from phase (b) as starting model and model to damp against (504). The portion of the model corresponding to the fast and resistive layer (505) is kept fixed and the weight of the low frequency data is decreased.

The progressive convergence of the inverted parameters (gray lines) to the “true” parameters (dark black lines) through the three phases (a) to (c) are shown as the 1D models are viewed from top to bottom.

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.

Claims

1. A computer-implemented method for estimating a physical properties model in a subsurface region using data comprising two or more geophysical data types, said model comprising numbers in a model parameter domain representing one or more physical properties, said method comprising: using a computer to jointly invert the data in a plurality of sequential phases, wherein in each phase until a last phase, only a portion of the data is inverted in order to infer a subset of the model's parameter domain, said portion being determined at least partly based on relative frequency content of the two or more data types,wherein in a first phase of the sequential phases, the portion of the data inverted is longer wavelength data from the data's frequency domain, and model parameters expected to vary at shorter wavelengths are frozen, further comprising, between said first phase and said last phase, one or more intermediate phases, each intermediate phase and the last phase progressively using more of the data frequency domain and inferring more of the model parameter domain, andwherein a last of the plurality of sequential phases yields a final subsurface physical property model that provides reservoir properties for the subsurface region; andprospecting for hydrocarbons with the reservoir properties of final subsurface physical property model.

2. The method of claim 1, wherein a converged model from one phase is used as a starting model, or a damping model, or both, for a next phase.

3. The method of claim 1, wherein the intermediate phases is topped when a pre-selected spatial and parameter resolution is reached.

4. The method of claim 1, wherein the joint inversion takes place in an inversion domain comprising, in addition to the model parameter domain and the data frequency domain, one or more of a data time domain, model wave number domain, and a model spatial domain.

5. The method of claim 4, further comprising, between said first phase and said last phase, one or more intermediate phases, each intermediate phase and the last phase progressively unfreezing more of: the model parameter domain, the model spatial domain, or the model wave number domain; andthe data frequency domain or the data time domain.

6. The method of claim 5, wherein as the sequential phases progress, one or more of the following occurs: the model wave number domain moves from a coarser inversion grid to a finer inversion grid;the data time domain moves from only earlier arrivals to including later arrivals;the model spatial domain moves from freezing deeper portions of the model to solving for the deeper portions, or from freezing regions of the model expected to vary rapidly to solving for such rapidly varying regions;the model parameter domain moves from freezing model parameters for a physical property expected to vary at shorter wavelengths to solving for model parameters for the physical property expected to vary at shorter wavelengths and four model parameters and for model parameters for any other of said one or more physical properties; andthe data frequency domain progresses from lower frequency data to include higher frequency data, which is accomplished by using a filter or by selecting from among the two or more data types by frequency content.

7. The method of claim 6, wherein the selecting from among the two or more data types by frequency content is accomplished by adjusting weights in an objective function used in the joint inversion to measure misfit between measured data and model-predicted data, said objective function having a term for each of the two or more data types, each term containing an adjustable weight.

8. The method of claim 7, wherein the two or more geophysical data types are selected from a group consisting of active seismic reflection data, a high-frequency data type; electromagnetic data, either controlled source or magneto-telluric, both low-frequency data types; and gravity measurements, also a low-frequency data type.

9. The method of claim 8, wherein the weight for a low-frequency data type is progressively decreased at each successive phase of the joint inversion through said last phase.

10. The method of claim 9, wherein the two or more geophysical data types are seismic data and either electromagnetic or magneto-telluric data, the physical properties are rock porosity, lithology, and water saturation, lithology being parameterized as shale volume fraction, and model parameters for water saturation are frozen through said last phase of the joint inversion, further comprising an additional phase after said last phase during which all model parameters except for water saturation are fixed, a strong damping is applied constraining the model to agree with a starting model obtained from the last phase of the joint inversion, and a large weight is given to the electromagnetic or magneto-telluric data compared to the seismic data.

11. The method of claim 1, wherein the one or more physical properties are rock porosity, lithology, and water saturation, lithology being parameterized as shale volume fraction, and each sand layer or interval identified during the joint inversion is assigned a single value for water saturation.

12. The method of claim 1, wherein the one or more physical properties include water saturation, and the model parameters for water saturation are frozen for at least a first phase of said plurality of sequential phases.

13. The method of claim 1, wherein in each phase of the plurality of sequential phases, the portion of the data inverted corresponds to a selected frequency band.