GB2585216A - Improved inversions of geophysical data - Google Patents

Abstract

A method of obtaining a three-dimensional model of a sub-surface region, comprises providing S1 a set of geologically realistic models, constructing a reference matrix using said models, and describing the reference matrix using vectors e.g. by principal component analysis. A subset of the most important vectors is obtained S2 which describe the reference matrix to an approximation. An initial model is defined S3. Upgoing energy wave data e.g. electromagnetic or seismic is recorded S4 at a plurality of spaced apart locations. A change to said initial model is obtained S5 to provide an improved model, the change minimising a penalty defined as a function of a difference between upgoing energy wave data simulated at said locations using the improved model and the recorded upgoing energy wave data, wherein the improved model comprises a combination of said most important vectors. Improved models are iteratively generated S6 until an appropriately accurate model is obtained.

Description

Improved Inversions of Geophysical Data

Technical Field

The present invention relates to a method of estimating properties of a geological structure and more specifically to non-linear inversions of geophysical data including electromagnetic and seismic data.

Background

Seismic imaging has long been used to determine the geological structure of a subsurface region, for example as an aid to the identification of hydrocarbon bearing regions. The basis of all seismic imaging techniques is the generation at the surface, or in a body of water above the surface, of acoustic sound waves which propagate through the surface and the detection of upgoing waves at the surface (or again within the body of water) which result from reflection of the downgoing waves from subsurface reflectors. The term upgoing wave' refers generally to a wave which has been reflected at a sub-surface location, but should not be interpreted narrowly as travelling only in vertical direction because the wave will propagate in all directions. The process of generating an image of the sub-surface from the seismic data is known as "inversion".

More recently, imaging processes that make use of electromagnetic (EM) techniques have been developed. These processes generally make use of a source of electromagnetic energy at the surface (or possibly above the surface or downhole) and a plurality of EM detectors at or close to the surface. One particular EM technique used for sub-sea imaging is known as Controlled Source EM (CSEM) and uses an active source towed behind a boat [Eidesmo, T., Ellingsrud, S., MacGregor, L. M., Constable, S.,. Sinha, M. C, Johansen, S., Kong, F. N. and Westerdahl, H., 2002, Sea Bed Logging (SBL), a new method for remote and direct identification of hydrocarbon filled layers in deepwater areas, First Break (20), pp 144-152.]. Again, the process of generating a sub-surface image (in effect a resistivity profile of the sub-surface) from the EM data is known as "inversion".

Imaging is often enormously computationally expensive, meaning that it is both financially expensive and time consuming to perform. If an assumption turns out to be incorrect, it might be several days before that assumption can be corrected and an analysis re-run. Taking shortcuts can lead to poor imaging quality leading, in a worst-case scenario, to drilling wells in a poor location or making incorrect assumptions about existing hydrocarbon reserves.

In both seismic and EM (specifically CSEM) approaches one may attempt to build a geological model of the sub-surface which, when the known seismic or EM energy is injected into it, generates simulated detected data that closely matches the actual detected data. Such approaches start with a basic geophysical model and iteratively adjust the model until the simulated and detected results converge at which point one can assume that the model closely approximates the real sub-surface geology. A finite difference or finite element method may be used to generate the simulated data using the model and the input data. This process is an example of non-linear inversion.

Taking an EM approach as an example, the model consists of a 3-dimensional representation of the sub-surface made up of a large number of cells. Each cell is assigned a resistivity (or its inverse, conductivity) value. At each iterative step, the simulated results are compared with the recorded data to generate some difference measure and the model adjusted to reduce the difference. The process of selecting the adjustment to make can be very computationally expensive as it typically requires adjusting the value of each cell in turn to determine which adjustment generates the best reduction in the determined difference.

Even if the principle is very simple there is a large difference between existing non-linear inversion methods. The most advanced methods (i.e. the most advanced formulas applied to update the model) are basically too computer intensive to be applied to the large-scale imaging problems that are met in petroleum exploration. The algorithms are therefore simplified at the cost of less efficient model updates. Many types of simplifications can be applied, and they can have large implications on the final quality of the images. There is a general trade-off between inversion quality and numerical efficiency or "cost" (i.e. how much disk memory is required and how much time it will take to get the inversion result).

Several formulas exist to calculate the adjustment. A very general one is the equation given by Tarantola (1984). With some slight adjustments, we can rewrite this equation as follows: mk-Ft = mk-UsCiril + Cm-1) 1(fCd-1 + (Ink -mo)) Hk gk where Mk and Mki-/ are vectors defining the current model and the subsequent improved model (in terms of resistivity values within each cell of the model), mo is the initial model, Ad is the difference between the simulated data in model ink and the recorded data, and J =CAd / anifr is the jacobian matrix of data derivatives with respect to the model parameters and * denotes a conjugate transpose. The information contained by J is important: it shows how the data difference is changed by changing the resistivity independently in each cell of the model. Cd-1 is the inverse of the covariance matrix of the data errors and Cirt is the inverse of the covariance matrix of model parameters describing a priori information on the model. Hr will be explained later.

If the 3-dimensional representation consists for example of one million cells, the inversion has to find one million unknowns. A first problem is the very large numerical cost this generates. A second problem is that the recorded CSEM data are not sufficient to properly constrain the problem: basically there are too few equations for the number of unknowns, so that an infinity of solutions exists. In other words, the inversion problem is ill-posed and this generates mathematical instability that can result in errors and artefacts in the obtained resistivity models. To reduce instability, one can use a mathematical approach called "regularization".

In the formulation of Tarantola explained above, the regularization is naturally obtained by using Cffil in the equations. Unfortunately, operator Cm*1 is too complicated and costly to calculate in practice for EM inversion and regularization is done by other means. Operator Ccii is also too complicated and costly to calculate and is typically replaced by a diagonal weight matrix.

Summary

According to a first aspect of the invention there is provided a computer implemented method of obtaining a model of a sub-surface region, the model comprising a three-dimensional matrix of sub-surface material property values, and comprising: a) providing a set of geologically realistic models, constructing a reference matrix using said models, and applying an analysis to that reference matrix to describe the reference matrix using vectors; b) selecting the most important vectors from said set of vectors to obtain a subset of vectors that describe said reference matrix to an approximation; c) defining as a first approximation an initial model; d) obtaining upgoing energy wave data recorded at a plurality of spaced apart locations within, at or above the surface, the upgoing energy waves depending on downgoing energy waves and the sub-surface material properties; e) determining a change to said initial model to provide an improved model, said change being one that seeks to minimise a penalty defined as a function of a difference between upgoing energy wave data simulated at said locations using the improved model and the recorded upgoing energy wave data, wherein the improved model comprises a combination of said most important vectors; and f) iteratively repeating step e) for each improved model until an appropriately accurate model is obtained.

The initial model may comprise a combination of said most important vectors.

Optionally, the analysis is a principal component analysis and said vectors are singular vectors. The singular vectors may comprise one hundred or fewer singular vectors. The number of vectors in said subset of vectors may be less than the number of models in said set of geologically realistic models.

The downgoing and upgoing energy waves may be electromagnetic energy waves, and said sub-surface material property values may be electrical resistivity or conductivity values. Alternatively, said downgoing and upgoing energy waves may be seismic energy waves, and said sub-surface material property values may be seismic velocity and density values.

The set of geologically realistic models may comprise five hundred or fewer geologically realistic models.

Step d) may comprise operating one or more energy sources at, above, or within a surface of said sub-surface region to produce said downgoing energy waves and operating detectors at said spaced apart locations to detect said upgoing energy waves, and recording the detected upgoing waves.

Alternatively, step d) may comprise operating one or more energy sources within a body of water above a surface of said sub-surface region to produce said downgoing energy waves and operating detectors at said spaced apart locations to detect said upgoing energy waves, and recording the detected upgoing waves.

According to a second aspect of the invention there is provided a method of selecting a region in which to drill in order to produce hydrocarbons and comprising applying the method according to the first aspect to obtain a model of a sub-surface region, and examining the model to identify locations within the sub-surface likely to produce hydrocarbons.

According to a third aspect of the invention there is provided a method of producing hydrocarbons and comprising selecting a region in which to drill using the method of the second aspect, drilling one or more wells in that region, and producing hydrocarbons from the well or wells.

According to a fourth aspect there is provided a computer program product comprising program instructions operative to be executed by a processor to cause the processor to perform a method according to the first aspect.

Brief Description of the Drawings

Figure 1 illustrates a set or reference models in a two dimensional vector space; Figure 2 also illustrates a set or reference models in a two dimensional vector space; and Figure 3 is a flow diagram illustrating a method of obtaining a model of a sub-surface region.

Detailed Description

It will be appreciated that the approaches described below will be executed primarily on computers comprising processors and memory. However, the size and capacity of both may be significantly reduced as compared with the implementation of prior art imaging approaches, and or imaging quality may be improved. It will also be appreciated that the approaches may be configured to provide images in renderable form for display on a computer display. Indeed, in some approaches this will be required in order to allow relevant experts to inspect the imaging results.

Non-linear inversion of CSEM data is a process by which an initial model of the subsurface is iteratively improved by reducing the misfit between data simulated in that model and the data recorded at seabed during an acquisition campaign. The model is typically represented as a coloured 3D image with small cubic pixels, where the colour represents the value of the resistivity (or its inverse, conductivity). At each iteration, the value of resistivity (or conductivity) in each pixel is updated to provide a better model.

What we mean here by "better" is a model that will result in simulated data that are more like the recorded data than at previous iteration, i.e. a model that reduces the data misfit.

To improve the current model, a small model update is added to it. To calculate the update, a large system of equations needs to be solved. The unknowns of the system are the value of the update for each pixel, and the equations are provided by the data measurements. Unfortunately, this system presents mathematical instability for several reasons. This can occur if we have too few data points, i.e. not enough equations, or because the equations are not sufficiently independent. One can also easily understand that it would not be possible in practice to estimate the value of resistivity in a small pixel that is very far away from the CSEM source and receivers because the influence of that pixel on the propagation of the EM wave through the earth is far too small.

When a system of equations has too many unknowns for the number of equations, there are an infinity of solutions. In the present case, this means that the model update provided by the formula can be very unstable: it may for instance provide a physically unrealistic model (i.e. the colours are wrong) or a model with a lot of artefacts (the image is very noisy). The instability can be treated in a mathematical way. For instance, adding a small constant to the diagonal of the system matrix will stabilize the system. This is a particular form of what is generally called Tikhonov regularization. Basically, regularization enforces a form of smoothness in the solution of the system. Here, it means that two neighbouring pixels will tend to have similar values rather than very different values. This reduces artefacts in the image.

Currently CSEM inversion methods typically use Tikhonov regularization or some variant of it. When using regularization, the inversion tries not to reduce the data misfit, but rather tries to reduce a cost function which is the sum of the data misfit and a second term which can be called "model penalty", and which is large for models that present a lot of contrasts. In practice, when running an inversion, the user has to adjust many parameters to tune regularization, and those can have a large impact on the final inversion result. For CSEM one tries for instance to choose parameters that favour vertical contrasts rather than horizontal contrasts to increase the chance to detect hydrocarbon reservoirs, which are typically horizontally extending resistive objects.

Regularization may be introduced as a way to stabilize inversion, which otherwise provides unstable results. Regularization introduces a second term representing a model penalty, in addition to the data misfit that inversion tries to reduce. This second term penalizes variations in the 3D image. In practice, the user chooses regularization parameters, runs the inversion (this can take days), checks the image, and tunes the regularization parameters to improve the image. This can be a difficult and time-consuming process.

In the classical approach to regularization (Tikhonov regularization), the model penalty term is high for rough models and small for smooth models. This naturally guides the inversion towards a solution which is a smooth model: regularization stabilizes the inversion by telling which types of models we prefer, which reduces the size of the solution space (most of the solutions among the infinity of solutions previously mentioned are not smooth at all). In the approach of Tarantola, the types of models that are preferred are models described by a prior model and a model covariance matrix Cm. This is in theory a better approach because it naturally guides the inversion towards realistic models rather than smooth models (in petroleum exploration, where hydrocarbons accumulations are present the hydrocarbon filled reservoir typically creates a strong discontinuity, i.e. the model is not smooth). The approach of Tarantola is unfortunately numerically too costly for many practical applications. The inventors have realised that the use of the covariance matrix can be replaced by the use of a set of vectors to have the benefit of the Tarantola approach at an affordable numerical cost. However, with that approach there may still be in the order of one million unknowns, with the related high computational cost.

The inventors have realised that the numerical cost can be reduced significantly by forcing the inversion to provide realistic models and at the same time to drastically reduce the number of unknowns, from one million to typically less than 100. When doing this, the number of numerical operations involved for calculating the model update is dramatically reduced.

Inversion seeks to find a model that minimizes a cost function. The model update done at each iteration depends on how this cost function is defined. The simplest approach is to define the cost function as the data misfit only: 0,1(m) =1 / 2Ad*Ad, (1) where 711 is a column vector of size Nm containing the value of the model (e.g. the resistivity or a simple function of resistivity, e.g. the conductivity) in each of the small 3D pixels, Ad -Ad (111) -d(ll) dr" is a column vector of size NJ representing the difference between simulated and recorded data, normally weighted by the inverse of the data uncertainty, and * represents complex conjugation. Equation (1) is a least-squares sum over the different data points contained in Ad.

When we want to find a model update Ant that gives a smaller misfit, i.e. OdOnT <41(m), we must understand how a change in the model changes the J = ed(m)/ simulated data. This information is contained in the Jacobian The first-and second-order derivatives of the cost function Od(n)are the gradient g(m) and the Hessian HOn) and one can easily show that g(m) = /*Ad The expression for the Hessian is more complex, but it can be shown that a good approximation is H (m)= J* . Taylor expansion of the cost function for model ink at iteration k gives gk = g,, + &An (2) When reaching a minimum for the cost function, its first-order derivatives must vanish.

Setting gk-1T ° in the above equation gives H Am = -g (3) where we have dropped the iteration index k for simplicity. This system of "normal equations" can be solved using, for example, a Conjugate Gradient approach, at each iteration, to get the model update Am. As explained above, however, this system is unstable. To stabilize the system, we replace it by (H + H,,),Arn = -g (4) A typical value for Hirg is a constant times the identity matrix. If the constant is sufficiently large, the system becomes stable enough.

Note that applying equation (4) rather than equation (3) is similar to replacing the cost function Od (/1) by 0(377) = Od(m) + (m) = Od (m)+ (in -mo)T Iir, (in -mo) (5) If we do that, the system of equations becomes (H +11,0Ani = -(g + mo)) (6) where mo is the initial model. The only difference between equations (4) and (6) is that m° is replaced by 377,, in equation (4), making the second term in the gradient vanish.

Both equations (4) and (6) normally give much more stable model updates than equation (3).

With Tikhonov regularization, matrix Mg can be written in the form H, = R (7) where T denotes transpose and where R is a square matrix typically containing a diagonal only (usually R is a simple multiple of the identity matrix, as introduced above) or a diagonal and a couple of subdiagonals and/or superdiagonals. Solving the system of equations (3), for a given iteration k, is normally done by an iterative method itself (typically a conjugate gradient solver) which involves calculating a lot of matrix-vector products of the form Hregv, where V is a vector the same size as m. Because of the special diagonal structure of matrix 1? these matrix-vector products have a low numerical cost, the same as a couple of vector-vector products.

Conventional (Tikhonov) regularization applied in inversion of CSEM data is therefore relatively computationally cheap to implement. However, it has the disadvantages explained above: the user needs to choose a number of parameters, must run inversions (often waiting several days), compare the results to what the user thinks a reasonable result would look like, and then adjust the parameters. In this workflow, the user uses a priori knowledge about resistivity after inversion is completed.

If we follow a Bayesian approach, like Tarantola (1984), we obtain equations very similar to equation (6), but with = (8) where Cm is the covariance matrix of the model parameters, which explicitly defines a priori information about the subsurface. As previously explained this approach is however numerically too costly for most practical applications.

In a specific embodiment, a set of geologically realistic models (the "reference set") is used to represent the prior information available about the subsurface properties. This is discussed below in the context of electromagnetic (specifically, CSEM) based methods where the model consists of a matrix of resistivity values. However, the approach is equally applicable to seismic approaches in which case the model may consist of a matrix of acoustic impedance values.

Non-linear inversion algorithm leads typically to one single solution but the interpreter must be aware that other solutions may be valid. In practice, one will normally run several inversions with slightly different stabilization parameters leading to different 3D images. The best images can be picked for more detailed interpretation. Running many inversion is however costly.

Because the stabilization often is based on simplistic assumptions that do not reflect the complexity of the earth (e.g. in reality spatial resistivity variations show abrupt variations at several locations in the subsurface) the solutions found by inversion can be quite unrealistic, i.e. can badly fit the geological understanding of the subsurface obtained from other type of knowledge, e.g. based on seismic data, on measurements of rock properties in wells or general knowledge about the geological setting. For instance, high resistivity may be expected at the location of a hydrocarbon-filled reservoir. The location of this potential reservoir may be known quite accurately from interpretation of seismic data and measurements in neighboring wells. However, the inversion may provide a solution where high resistivity in the 3D cube is erroneously imaged too shallow or too deep, at a depth where it is highly unlikely to find any hydrocarbons. In other words, the 3D image is geologically speaking unrealistic. The interpreter must consider the possibility of such errors, which complicates the interpretation process and reduces the value of the inversion result. New inversion runs with new sets of stabilization parameters may be required to cope with that problem, adding more computer costs.

The proposed solution is to describe the space of possible, realistic geological models beforehand to reduce the possibility of such errors. The inversion will be forced to look for models that are by construction respecting the available information available on the subsurface's geology. Such information is referred to as "a-priori" information, in the sense that it is available prior to any use of CSEM data.

The prior information is translated into a large number of possible 3D resistivity models, for example 250 resistivity models that we may call "prior models". These N realistic 3D prior resistivity models describe a linear space that we can call the prior space.

Extracting the principal components of these N 3D prior models, we can describe the prior space almost exactly with Nc singular vectors, where Nc is smaller than N, i.e. each of the N 3D models can be represented almost exactly by a linear combination of the N6 singular vectors. N6 represents the number of degrees of freedom in the space of possible 3D resistivity models built from the available prior information. Nc may typically be lower than 50.

If we assume that each reference model comprises a set of cells, e.g. 1 to Nm, then the reference set can be represented by a matrix M where each column of the matrix corresponds to a given reference model, or to the difference between a given reference model and the average of all reference models (the most probable model). For example, if each model comprises 1000000 cells and there are 250 models, then M will be a 1000000 X 250 matrix. By applying a singular value-decomposition (SVD) method to the matrix M, the matrix can be represented using singular vectors (Fo) which represent the principal components of the reference set of realistic models, i.e. M = FoSoKo, where So are the singular values and Ko are the coefficients for each reference model.

The singular vectors (F0) are ordered so that the first vectors are relatively more organised whilst the later ones are less organised, i.e. more chaotic. The singular values (So) are sorted in decreasing order and only the first few are retained. Whilst the matrix M now no longer exactly describes the reference models, the approximation can still be very good, i.e. M FSK, whereFconsists of the N6 first columns of Fo, S corresponds to the N6 first rows and columns of So and K corresponds to the Nc first rows of Ko.

It can also be expected that the approximation will be very good for any realistic model, even if not strictly belonging to the set of reference models.

The singular values S and coefficients K can be combined into a single matrix W with components which represents the weights w, for each of the reference models.

Once the singular vectors are calculated, the inversion is done by optimizing for N, coefficients, one coefficient for each singular vector. The obtained model belongs naturally to the prior space. For instance any subsurface model, consisting of resistivity or conductivity values in cells 1 to Nm, can now be described as m = map + FSc,or map + Fiv, where map is the average of the reference models (the most probable model) and c and w are column vectors with N, values, where N, is normally smaller than 50 while Nm is typically of the order of a million. The known inversion approaches explained above for m can be applied to w or c instead. Since the number of unknowns is strongly reduced, the computational burden is now much smaller. As soon as the unknown values in vectors Iv or c are found by inversion, the values in vector m are easily obtained by the formulae above.

Figure 1 illustrates the case where N, = 2 and vector w has only two weights, w1 and w2. All of the models of the reference set can be represented by these two weights to a reasonably good approximation, and are shown by the circles in the figure. The confidence ellipses correspond to one and two standard deviations and indicate the general spread of the reference models, i.e. the prior uncertainty on model properties. We see that weight w1 is of greater significance than weight w2 and has a larger prior uncertainty.

Figure 2 shows how the same set of reference models is represented by the coefficients c1 and c2 of vector c instead. The ellipses are now circles, indicating that the prior covariance matrix (or "the covariance matrix of model parameters") for vector c is the identity matrix. By analogy with equation (8), we see that matrix Hreg is in this case equal to the identity matrix.

Figure 3 illustrates a method, comprising the steps of: 51, providing a set of geologically realistic models, constructing a reference matrix using said models, and applying an analysis to that reference matrix to describe the reference matrix using vectors; S2 selecting the most important vectors from said set of vectors to obtain a subset of vectors that describe said reference matrix to an approximation; S3 defining as a first approximation an initial model; S4 obtaining upgoing energy wave data recorded at a plurality of spaced apart locations within, at or above the surface, the upgoing energy waves depending on downgoing energy waves and the sub-surface material properties; S5 determining a change to said initial model to provide an improved model, said change being one that seeks to minimise a penalty defined as a function of (i) a difference between upgoing energy wave data simulated at said locations using the improved model and the recorded upgoing energy wave data, wherein the improved model comprises a combination of said most important vectors; and S6 iteratively repeating step e) for each improved model until an appropriately accurate model is obtained.

A technical effect of the claimed method is a great reduction in computational cost. The computer runtime for a typical inversion process may be in the order of weeks or even months, with a corresponding amount of electrical energy consumed by the computer. In specific examples illustrated above, the number of unknowns may be reduced from one million to less than 100, and the computational runtime will be reduced correspondingly.

A further effect is an improved estimation accuracy because the models are restricted to geologically realistic models, without allowing the algorithm to arrive at solutions which produce a good match to the upgoing energy wave data while providing a model of the subsurface which geologically cannot be true. This improved estimation avoids situations like drilling wells in a poor location or making incorrect assumptions about existing hydrocarbon reserves. For example, hydrocarbons may be located in a trap formed by a suitably shaped medium which inhibits upward migration of hydrocarbons through permeable rock formation. If the shape of the medium is incorrectly determined, the location of the trapped hydrocarbons is also incorrectly determined and an incorrect drilling operation may be carried out.

It will be readily appreciated that the embodiments described above may be incorporated into an exploration process to identify regions likely to produce hydrocarbons. As part of that process data may be obtained by performing EM, seismic or other surveys in the area of interest. Of course, historical data may be reanalysed to provide new insight into a region. The result of the process may be the drilling of one or more wells within the location.

It will be appreciated by the person of skill in the art that various modifications may be made to the above described embodiments without departing from the scope of the present invention as defined in the claims. For example, whilst reference has been made above to EM and seismic imaging, the approach described may be applied to other forms of imaging, e.g. gravitation-magnetic (gravmag) imaging. In another modification, rather than using SVD, a factor analysis may be used.

Claims (13)

1. Claims 1. A computer implemented method of obtaining a model of a sub-surface region, the model comprising a three-dimensional matrix of sub-surface material property values, and comprising: a) providing a set of geologically realistic models, constructing a reference matrix using said models, and applying an analysis to that reference matrix to describe the reference matrix using vectors; b) selecting the most important vectors from said set of vectors to obtain a subset of vectors that describe said reference matrix to an approximation; c) defining as a first approximation an initial model; d) obtaining upgoing energy wave data recorded at a plurality of spaced apart locations within, at or above the surface, the upgoing energy waves depending on downgoing energy waves and the sub-surface material properties; e) determining a change to said initial model to provide an improved model, said change being one that seeks to minimise a penalty defined as a function of a difference between upgoing energy wave data simulated at said locations using the improved model and the recorded upgoing energy wave data, wherein the improved model comprises a combination of said most important vectors; and f) iteratively repeating step e) for each improved model until an appropriately accurate model is obtained.
2. 2. The method of claim 1, wherein said initial model comprises a combination of said most important vectors.
3. 3. A method according to claim 1 or 2, wherein said analysis is a principal component analysis and said vectors are singular vectors.
4. 4. A method according to claim 3, wherein said singular vectors comprises one hundred or fewer singular vectors.
5. 5. A method according to claim 1, wherein the number of vectors in said subset of vectors is less than the number of models in said set of geologically realistic models.
6. 6. A method according to any one of the preceding claims, wherein said downgoing and upgoing energy waves are electromagnetic energy waves, and said sub-surface material property values are electrical resistivity or conductivity values.
7. 7. A method according to any one of claims 1 to 5, wherein said downgoing and upgoing energy waves are seismic energy waves, and said sub-surface material property values are seismic velocity and density values.
8. 8. A method according to any one of the preceding claims, wherein said set of geologically realistic models comprises five hundred or fewer geologically realistic models.
9. 9. A method according to any one of the preceding claims, wherein step d) comprises operating one or more energy sources at, above, or within a surface of said sub-surface region to produce said downgoing energy waves and operating detectors at said spaced apart locations to detect said upgoing energy waves, and recording the detected upgoing waves.
10. 10. A method according to any one of claims 1 to 8, wherein step d) comprises operating one or more energy sources within a body of water above a surface of said sub-surface region to produce said downgoing energy waves and operating detectors at said spaced apart locations to detect said upgoing energy waves, and recording the detected upgoing waves.
11. 11. A method of selecting a region in which to drill in order to produce hydrocarbons and comprising applying the method of any one of the preceding claims to obtain a model of a sub-surface region, and examining the model to identify locations within the sub-surface likely to produce hydrocarbons.
12. 12. A method of producing hydrocarbons and comprising selecting a region in which to drill using the method of claim 11, drilling one or more wells in that region, and producing hydrocarbons from the well or wells.
13. 13. A computer program product comprising program instructions operative to be executed by a processor to cause the processor to perform a method according to any of claims 1 to 12.