GB2570330A

GB2570330A - Regularization of non-linear inversions of geophysical data

Info

Publication number: GB2570330A
Application number: GB1800980.3A
Authority: GB
Inventors: Causse Emmanuel; Iren Nordskag Janniche
Original assignee: Equinor Energy AS
Current assignee: Equinor Energy AS
Priority date: 2018-01-22
Filing date: 2018-01-22
Publication date: 2019-07-24
Anticipated expiration: 2038-01-22
Also published as: GB201800980D0; WO2019143253A1; GB2570330B

Abstract

A method of obtaining a 3D model of a sub-surface region comprises providing S1 a set of geologically realistic models, constructing a reference matrix using said models, and then describing the reference matrix using vectors. The most important vectors (e.g. singular vectors from a principal component analysis S2) are selected S3 to describe said reference matrix to an approximation. Upgoing wave energy is recorded S4 e.g. seismic, electromagnetic. An initial model of said sub-surface region is obtained S5. A change to the initial model is obtained S6 to provide an improved model. The change seeks to minimise a penalty defined as a function of (i) a difference between upgoing energy wave data simulated at locations using the improved model and the recorded upgoing energy wave data, and (ii) a measure of the difference between the improved model as defined approximately using said subset of vectors and a most probable model. This procedure is repeated iteratively S7 for each improved model until an appropriately accurate model is obtained S8.

Description

Regularization of Non-linear Inversions of Geophysical Data

Technical Field

The present invention relates to a method of regularization of 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 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 nonlinear 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:

m_k+1 = + C”¹) Ad + C^^l(m_k - m₀))

H_k g_k where and m^i 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 and the recorded data, and J = dAd I dm_k 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. Cjj¹ is the inverse of the covariance matrix of the data errors and Cm^x 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 C_m ¹ in the equations. Unfortunately, operator C_m ¹ is too complicated and costly to calculate in practice for EM inversion and regularization is done by other means. Operator Cd”¹ is also too complicated and costly to calculate and is typically replaced by a diagonal weight matrix. So the Tarantola approach is replaced by simpler approaches involving simplifications and approximations.

Summary

It is an object of the present invention to formulate non-linear inversion regularization in such a way that is numerically much less costly to implement than in Tarantola (1984), but still presents the same qualities that are important for obtaining good, geologically reasonable 3D images of the subsurface.

According to a first aspect of the present invention there is provided a method of obtaining a model of a sub-surface region, the model comprising a three-dimensional matrix of sub-surface material property values. The method comprises:

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) 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;

d) defining as a first approximation an initial model of said sub-surface region;

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 (i) a difference between upgoing energy wave data simulated at said locations using the improved model and the recorded upgoing energy wave data, and (ii) a measure of the difference between the improved model as defined approximately using said subset of vectors and a most probable model; and

f) iteratively repeating step e) for each improved model until an appropriately accurate model is obtained.

The analysis may be a principal component analysis with said vectors being singular vectors. Other singular vector approaches may alternatively be used.

The subset of singular vectors may comprise one hundred or fewer singular vectors, i.e. a number that significantly reduces computational expense. 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 electromagnetic energy waves, and said sub-surface material property values are electrical resistivity or conductivity values. Alternatively, said downgoing and upgoing energy waves may be seismic energy waves, and said sub-surface material property values are seismic velocity and density values.

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

At step e), determining a change may comprise solving a system of linear equations that link the difference between upgoing energy wave data simulated at said locations using the improved model and the recorded upgoing energy wave data and a search direction along which the model could optimally be updated.

The linear equations may be “normal equations”. Solving the system of linear equations may comprise using a conjugate gradient algorithm.

At step e), the step of determining a change may alternatively comprise calculating a gradient of (i) with respect to perturbations of the model.

The most probable model may be:

an average of said set of reference models; or a model that is perturbed to generate said set of reference models.

Step c) 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.

Step c) 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 present invention there is provided a method of selecting a region in which to drill in order to produce hydrocarbons and comprising applying the method of the above first aspect of the invention to obtain a model of a sub-surface region, and examining the model to identify locations within the subsurface likely to produce hydrocarbons.

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

Brief Description of the Drawings

Figure 1 illustrates a two dimensional vector space constructed by applying a SVD analysis to a set or reference models;

Figure 2 is a further view of the vector space of Figure 1, specifically indicating the relative importance of vector weights;

Figure 3 illustrates the vector space of Figure 1 specifically indicating an external difference between a model and its projection onto the vector space; and

Figure 4 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 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.

In summary, regularization is 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 approach of Tarantola, inversion is formulated differently. It uses a Bayesian formalism considering the statistics of the a priori information on the model. Instead of minimizing the data misfit only, one tries to find the model that maximizes the posterior probability. The posterior probability is a function of the prior probability of the model and a likelihood function that is high when the data misfit is small. This is a more physical approach in the sense that CSEM data comes simply as additional information modifying prior knowledge (by essence geologically reasonable) about the subsurface: the final model should be probable (consistent with prior information) and at the same time it should explain the recorded data by providing a sufficiently low data misfit. If such a model is found, it will be much better than an unphysical model. In this formulation, a cost function is reduced, which is a sum of the data misfit and a model penalty that is small for models in agreement with the available information we have on the subsurface, and high for models that are different from what we expect. Stabilization comes naturally: if there is an infinity of solutions, we prefer solutions that agree with what we already know about the subsurface. This approach leads naturally towards results that are stable and geologically reasonable. The model penalty term represents the distance to what we think is the most probable model, where the distance is defined with a covariance matrix. Unfortunately, this covariance matrix is too large to easily be calculated, inverted and used in practice.

The approach proposed here reformulates the model penalty term so that it is nearly equal to the expression using the covariance matrix, but uses a set of vectors instead. Matrix-vector operations are replaced by a limited set of vector-vector operations. This drastically reduces the number of operations and memory requirements involved and makes the inversion feasible.

Further detailed mathematical consideration of known approaches

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:

φα(ηί) = 1/2Δί/*Δί/, (1) where ^m is a column vector of size 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(w) - d(w) - d_{rec js g C0}|_{umn vec}t_or of size ^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 leastsquares sum over the different data points contained in .

When we want to find a model update Am that gives a smaller misfit, i.e. φα(ιη + Διη) <φ_ά(ιη), _{we must unc}|_erst_anc| how _a change in the model changes the simulated data. This information is contained in the Jacobian J -dd(m) I dm. γρθ first- and second-order derivatives of the cost function are the gradient and the Hessian 77 (w) and one can easily show that = J Ad. γ^θ _exp_ressj_onfor the Hessian is more complex, but it can be shown that a good approximation is H(m) = J J 7_ayior expansion of the cost function for model ^mk at iteration gives gk₊i = gk + H_khm_k (2)

When reaching a minimum for the cost function, its first-order derivatives must vanish. Setting “θ in the above equation gives

Am = -g , ₍₃₎ 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 (/7+/7,,,.,) Am =-g ₍₄₎ u

A typical value for '-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 ^4⁷⁷⁷) by

0(m) = φ_ά (m) + (m) = φ_ά (m) + (m - m₀ )^T H_reg (m-m₀)

If we do that, the system of equations becomes (H + H_reg)km = -(g + H_reg (m - m₀)) where ^m° is the initial model. The only difference between equations (4) and (6) is that ^m° is replaced by ^mk 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).

TT

With Tikhonov regularization, matrix ^reg can be written in the form

H_reg = R R (Y) 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 , is normally done by an iterative method itself (typically a conjugate gradient solver) which involves calculating a lot of matrixvector products of the form where ^v is a vector the same size as ^m . Because of the special diagonal structure of matrix R 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 π —

GZ reg — '—'m (θ) where is the covariance matrix of the model parameters, which explicitly defines a c^_1priori information about the subsurface. If the data are properly weighted and is known accurately, it should not be needed to adjust any parameter: the prior information is used before inversion is run. In practice, some small adjustments may c^_1 be done, e.g. by multiplying by an adjustable scalar parameter, but this represents a much simpler set of parameters than for Tikhonov regularization anyway. Unfortunately, we cannot simply update a CSEM inversion code that uses equation (7) by replacing equation (7) with equation (8). Some of the reasons are listed below:

• X is in practice not known analytically.

• In most practical cases has components far from the diagonals and is singular, meaning that its inverse does not exist.

• Even if an approximate inverse is calculated (e.g. by adding a constant diagonal to an estimate of C™ calculated from a sufficiently large population of models), this would normally have a lot of non-diagonal terms, so that products of the form could not be simplified and would result in numerically very costly matrixvector products.

The third point is probably the main reason why this approach is not used for 3D imaging of the subsurface. Instead, the cheaper Tikhonov regularization is used, despite its disadvantages.

Proposed solutions

In common with the approach proposed by Tarantola and discussed briefly above, it is proposed here to use a set of geologically realistic models (the “reference set”) to steer the inversion process. 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.

If we assume that each reference model comprises a set of cells, e.g. 1 to n, then the reference set can be represented by a matrix M where each column of the matrix corresponds to a given reference model. For example, if each model comprises 35 cells and there are 100 models, then M will be a 35 X 100 matrix. By applying a singular value-decomposition (SVD) method to the matrix M, the matrix can be represented using singular vectors (F_o) which represent the principal components of the reference set of realistic models, i.e. M = F^S^K^, where S_o are the singular values and K_o are the coefficients for each reference model.

The singular vectors (F_o) 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 (S_o) 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. 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. In other words, any proposed model can be projected onto the vectorial space spanned by the first few singular vectors and, if this proposed model is realistic, it can be expected that the approximation error (the residual) will be small. In contrast, if the residual is large it is likely that the model is not a probable one.

The singular values S_o and coefficients K_o can be combined into a single matrix W_owith components which represents the weights w, for the reference models. Again, only the first few coefficients for each model need to be considered and the remainder can be neglected. Figure 1 illustrates the case where only the first two weights, w/ and w₂, are considered. All of the models of the reference set can be represented by these two weights to a reasonably good approximation, with the discrete points in the figure indicating these approximations. The ellipse shown in the figure indicates the general spread of the approximations of the reference models. In this example the ellipse indicates that weight it’/ is of greater significance than weight w₂.

Of course any model can be projected onto the vectorial space. In Figure 1 a model m is shown projected onto the space as m. We can define a “regularisation penalty” Φ(ηί) in terms of the length of the residual = (m,m) and the distance to the most probable model d^ = i.e. Φ(ηι) is a function of d_(e) and d_(i).

The second of these terms may take into account the spread of the reference models as indicated by the ellipse. If the extent of the ellipse in the direction of axis w, is λ; and in axis w₂ is λ₂. then d_(i) may be defined as follows:

4) = (9) where N is the number of singular vectors and w, are the weights for the projection m. This equation biases the regularisation penalty to favour the most probable models and is advantageous, compared with other possible equations, as it keeps only meaningful (organized) terms and omits other (chaotic) terms, and it has relatively good numerical efficiency due to having fewer terms. Figure 2 illustrates this approach where the adjustment using the values λ means that, for example, the models m_lt m₂ and m₃ are all equally probable.

Considering the first regularisation penalty term d_(e), this is illustrated in Figure 3, where m = m + Δτη. The term d_(e) may be defined as:

d² _ey = αΔτη^τΔτη (10)

For a particular model m belonging to the set of reference models, m = m₀ + where/} is singular vector number / (column i of matrix F_o). In this case, = ft and Δτη^τΔτη = ^_=Nc+iwf

Comparing these equations, we can assume that a = (l/A_Wc+1)².

By way of example the regularisation penalty for a model m may be equal to the sum of the squares of the two distances (i.e. the internal and external distances), i.e.:

0^(^) = ^)(^) + ^)(^)

The first term corresponds to the internal error defined as

Nc

^)(^)=2(^/^)² , (12) where the coefficients ^Wj represent the projection ^m of model ^{m m}° on the basis formed by the first singular vectors, i.e.

m = Fw (13)

Because of the orthogonality of matrix F , these coefficients are easily formulated as a function of ^m :

w = F^T(m-m₀)

In the two last equations, ^w is a column vector formed by the F, coefficients ^W). Realizing that the values ^{1Z Wj} correspond to the terms on the diagonal of matrix we can apply (14) and rewrite equation (12) as dy₎(m) = (m-m_oy FS~²F^T(m-m₀)

This internal error is related to how probable or how “geological” model ^m° ^+m is.

The second term in equation (11) represents an external error related to the size of the residual vector +m) γ^θ fesjciugi vector represents what is unexpected in model ^m . It then follows that:

d² _e)(m) = (1/λ^₊ι) Sm^T Sm

We have

Sm = m-m₀-m = m-m₀- Fw = m-m₀ -FF^T(m-m₀) = (I -FF^T )(m-m₀ ).

Therefore, equation (16) can be rewritten as:

<7<² > (m) = (1 / /1¾+1) (m - m0 )^r (7 - FF^T)² (m-m0)^T

Because of the orthogonality of F , the squared matrix in the middle can be simplified: d^{ni) = (l/A^₊i) (m-m₀)^T (7-FF^T){m-m0)^T (17)

This internal error penalizes what is unexpected in ^m , i.e. the “surprises”. Therefore it is scaled with a larger number than the previous components in the internal error (

F is larger for 7-A_c+l _t|_ian f_or |_{Ower va}|_{ues O}f 7 because the eigenvalues quickly decrease). Thus, only small penalties are applied to reasonable models while the surprises are treated as with Tikhonov regularization. If we want to generalize equation (17) to the more general form of Tikhonov regularization we can use (m) = (1 / λ^₊₁ )^r(7 -FF^T)RR(/ -FF^T) (m -m₀)^T θ) by analogy with equation (7). This represents an alternative to equation (17).

Combining the internal and external error definitions in equations (15) and (17) we see

TT that our invention corresponds to a new definition of matrix ^reg , i.e.:

H_reg=FS-²F^T + (l/^_c+1)(7-FF^r)

Interpretation of equation (19)

In equation (19), the first term corresponds to the internal error. This treats what is organized and predictable in the model. The second term penalizes what is unpredictable and not organized in the model. It can be proven mathematically that equations (19) and (8) provide very similar model penalties provided a sufficient number of eigenvalues is used. In practice, since the eigenvalues decrease very quickly this number is expected to remain small for most applications.

The similarity of these penalties is confirmed by tests. In these tests, we have used an example where the covariance matrix is not singular and its inverse can therefore be calculated. This allows a comparison of the penalties for different models, measured with equations (8) and (19). This has been done for different values of F. When ^cis equal to or larger than 4, the differences between equations (8) and (19) is minimal. We can consider equation (19) as an approximation of the inverse covariance matrix ^m used in equation (8). This approximation is not perfect, but it is completely stable, so we can in principal use it also in the case where the covariance operator is singular.

With our approach we can also obtain approximations of the covariance operator itself (not only of its inverse) and we can confirm that those are very good when is equal to or larger than 4. In practice, the covariance operator is not known. For small problems like in our tests it can be estimated, but the estimate is probably not better than the approximations we get.

Numerical efficiency of equation (19)

The numerical efficiency of equation (19) is related to the fact that matrix F has a small number of columns : each time we need a matrix-vector product when using equation (8), we need instead vector-vector products when using equation (19). In terms of numerical operations this replaces numerical operations by operations. ^^,n is typically equal to several millions. should be equal to the effective number of degrees of freedom of the prior model space, i.e. to the number of independent parameters required to describe our a priori information. This is expected to be lower than 100 in most situations and significantly smaller values than 100 can probably be used in practice. So, while our approach gives similar results to the approach using covariance operators, it is numerically much more efficient (the problem is reduced by a factor with order of magnitude 10000). The fact that it is possible to increase numerical efficiency so drastically is related to the fact that the number of degrees of freedom of the prior information is limited (less than 100), while the model is described with millions of parameters. Embodiments of the present invention exploit this large scale difference.

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

1. A 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:

g) 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;

h) selecting the most important vectors from said set of vectors to obtain a subset of vectors that describe said reference matrix to an approximation;

i) 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;

j) defining as a first approximation an initial model of said sub-surface region;

k) 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, and (ii) a measure of the difference between the improved model as defined approximately using said subset of vectors and a most probable model; and

l) iteratively repeating step e) for each improved model until an appropriately accurate model is obtained.

2. A method according to claim 1, wherein said analysis is a principal component analysis and said vectors are singular vectors.

3. A method according to claim 2, wherein said subset of singular vectors comprises one hundred or fewer singular vectors.

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

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

6. A method according to any one of claims 1 to 4, wherein said downgoing and upgoing energy waves are seismic energy waves, and said sub-surface material property values are seismic velocity and density values.

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

8. A method according to any one of the preceding claims, wherein said penalty is defined as the sum of the squares of (i) and (ii).

9. A method according to any one of the preceding claims, wherein, at step e), determining a change comprises solving a system of linear equations that link the difference between upgoing energy wave data simulated at said locations using the improved model and the recorded upgoing energy wave data and a search direction along which the model could optimally be updated.

10. A method according to claim 9, wherein the linear equations are “normal equations”.

11. A method according to claim 10 and comprising solving the system of linear equations using a conjugate gradient algorithm.

12. A method according to any one of claims 1 to 8, wherein, at step e), determining a change comprises calculating a gradient of (i) with respect to perturbations of the model.

13. A method according to any one of the preceding claims, wherein said most probable model is:

is an average of said set of reference models; or is a model that is perturbed to generate said set of reference models.

14. A method according to any one of the preceding claims, wherein step c) 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.

15. A method according to any one of claims 1 to 13, wherein step c) 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

5 at said spaced apart locations to detect said upgoing energy waves, and recording the detected upgoing waves.

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

10 model of a sub-surface region, and examining the model to identify locations within the sub-surface likely to produce hydrocarbons.

17. A method of producing hydrocarbons and comprising selecting a region in which to drill using the method of claim 16, drilling one or more wells in that region, and

15 producing hydrocarbons from the well or wells.