Method of, and apparatus for, full waveform inversion

Abstract

A method of subsurface exploration includes generating a representation of a portional volume of the Earth from a seismic measurement of a physical parameter. The method includes providing observed seismic dataset having three distinct nonzero data values derived from three distinct nonzero seismic measured values of said portional volume of the Earth, generating a predicted seismic dataset having three distinct nonzero data values, generating a nontrivial convolutional filter including three nonzero filter coefficients, generating a convolved observed dataset by convolving the convolutional filter with said observed seismic dataset, generating primary objective functions to measure the similarity between said convolved observed dataset and said predicted dataset, maximizing and/or minimizing said primary objective functions by modifying at least one filter coefficient of the convolutional filter, generating predetermined reference filters having at least three reference coefficients generating secondary objective functions to measure the similarity between filter coefficients for the nontrivial filter and reference coefficients for the predetermined reference filters, and minimizing and/or maximizing said secondary objective functions by modifying a model coefficient of a subsurface model of a portion of the Earth to produce an updated subsurface model of a portion of the Earth.

Claims

1. A method of subsurface exploration, the method comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement of at least one physical parameter, and comprising the steps of: providing an observed seismic data set comprising at least three distinct non-zero data values derived from at least three distinct non-zero seismic measured values of said portion of the volume of the Earth; generating, using a subsurface model of a portion of the Earth comprising a plurality of model coefficients, a predicted seismic data set comprising at least three distinct non-zero data values; utilizing one or more non-trivial convolutional filters, each filter comprising three or more non-zero filter coefficients and configured to, when convolved with the observed seismic data set, maximize the similarity or minimize the mismatch between the convolved observed data set and the predicted data set; minimizing and/or maximizing one or more objective functions operable to measure the similarity and/or mismatch between the filter coefficients for the or more non-trivial convolutional filters and at least three reference coefficients of one or more reference filters by modifying at least one model coefficient of said subsurface model of a portion of the Earth to produce an updated subsurface model of a portion of the Earth; and providing an updated subsurface model of a portion of the Earth for subsurface exploration.

2. The method according to claim 1, wherein said at least one non-trivial convolutional filter is operable to transform at least a portion of said observed seismic data set to render the observed seismic data set and predicted data set approximations of one another.

3. The method according to claim 1, wherein the method further comprises: utilizing at least one further non-trivial convolutional filter, the at least one further non-trivial convolutional filter comprising three or more non-zero filter coefficients and configured to, when convolved with the predicted seismic data set, maximize the similarity or minimize the mismatch between the convolved observed data set and the convolved predicted data set.

4. The method according to claim 3, wherein the least one non-trivial convolution filter is operable to transform at least a portion of said observed seismic data set and the at least one further nontrivial convolutional filter is operable to transform said predicted data set to render the observed seismic data set and predicted data set approximations of one another.

5. The method according to claim 3, wherein a plurality of further non-trivial convolutional filters is provided.

6. The method according to claim 1, further comprising generating one or more normalized objective functions operable to measure the similarity and/or mismatch between the filter coefficients for the one or more non-trivial convolutional filters and the reference coefficients for the one or more reference filters, wherein the one or more normalized objective function is insensitive or has reduced sensitivity to the relative scaling of the filter coefficients for the one or more non-trivial convolutional filters with respect to the reference coefficients for the one or more reference filters.

7. The method according to claim 6, wherein at least one of said objective functions is normalized by a norm of the filter coefficients of the at least one non-trivial convolutional filter.

8. A method according to claim 7, wherein at least one of said objective functions is normalized by the L.sub.1-norm or the L.sub.2-norm of the filter coefficients of the at least one non-trivial convolutional filter.

9. The method according to claim 1, wherein a plurality of non-trivial convolutional filters is provided.

10. The method according to claim 1, wherein a plurality of observed seismic data sets and/or a plurality of predicted seismic data sets is provided.

11. The method according to claim 1, wherein one or more of the non-trivial convolutional filters involve convolution in time or space.

12. The method according to claim 1, wherein one or more of the non-trivial convolutional filters involve convolution in more than one dimension.

13. The method according to claim 1, wherein one or more of the non-trivial convolutional filters are selected from the group consisting of: Wiener filters, Kalman filters, least mean square filters, normalized least mean square filters, and recursive least squares filters.

14. The method according to claim 1, wherein one or more of the non-trivial convolutional filters are calculated in the time domain, the temporal-frequency domain, the space domain, or the wave-number domain.

15. The method according to claim 1, wherein the filter coefficients are regularised and/or interpolated between said data sets and/or subsets of said data sets, and/or interpolated between filters.

16. The method according to claim 1, wherein the objective function involves the application of a weighting to the filter coefficients of said at least one non-trivial convolutional filter, and wherein the weighting is dependant upon the temporal lag of said at least one non-trivial convolutional filter.

17. The method according to claim 16, wherein the objective function consists of one or more of: a norm of the weighted convolutional filter coefficients in combination with a norm of the unweighted filter coefficients; a norm of the weighted convolutional filter coefficients divided by a norm of the unweighted filter coefficients; or a norm of the unweighted convolutional filter coefficients divided by a norm of the weighted filter coefficients.

18. The method according to claim 1, wherein at least one of said objective functions comprises one or more of a norm misfit objective function, an L.sub.1-norm misfit objective function, or a least-squares misfit objective function.

19. The method according to claim 1, wherein the step of minimizing the one or more objective functions with respect to said model parameters comprises utilizing the gradient of the objective.

20. The method according to claim 1, wherein at least one of the reference filters is a convolutional filter which leaves the dataset to which it is applied unaltered, or approximately unaltered, other than by one or more of: multiplication by a constant, by a shift in time or space, and by a rotation of phase.

21. The method according to claim 1, wherein at least a portion of the observed seismic dataset is numerically propagated to a subsurface region of the model that is spatially removed from the location at which they were originally recorded.

22. The method according to claim 1, further comprising repeating the steps of utilizing and minimizing and/or maximizing for the updated subsurface model of a portion of the Earth until a convergence criterion is met.

23. The method according to claim 1 wherein, further to said step of providing, the method comprises utilising said updated model for sub-surface exploration.

24. The method according to claim 1, wherein said at least one physical parameter comprises pressure, particle velocity or displacement.

25. The method according to claim 1, wherein the observed data set and the predicted data set comprise values of a plurality of physical parameters.

26. The method according to claim 1, wherein said at least one convolutional filter is operable to render the cross-correlation of at least a portion of the observed seismic dataset with at least a portion of the predicted seismic dataset and the cross-correlation of at least a portion of the predicted seismic dataset with at least a portion of the observed seismic dataset approximations of one another.

27. A method of subsurface exploration, the method comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement of at least one physical parameter, and comprising the steps of: providing an observed seismic data set comprising at least three distinct non-zero data values derived from at least three distinct non-zero seismic measured values of said portion of the volume of the Earth; generating, using a subsurface model of a portion of the Earth comprising a plurality of model coefficients, a predicted seismic data set comprising at least three distinct non-zero data values; utilizing one or more non-trivial convolutional filters, the one or more non-trivial convolutional filters comprising three or more non-zero filter coefficients and configured to, when convolved with the predicted seismic data set, maximize the similarity or minimize the mismatch between the convolved predicted data set and the observed data set; minimizing and/or maximizing one or more objective functions operable to measure the similarity and/or mismatch between the filter coefficients for the one or more non-trivial convolutional filters and at least three reference coefficients for one or more reference filters by modifying at least one model coefficient of said subsurface model of a portion of the Earth to produce an updated subsurface model of a portion of the Earth; and providing an updated subsurface model of a portion of the Earth for subsurface exploration.

28. The method according to claim 27, wherein said at least one non-trivial convolutional filter is operable to transform at least a portion of said predicted seismic data set to render the observed seismic data set and predicted data set approximations of one another.

29. The method according to claim 27, further comprising: utilizing at least one further non-trivial convolutional filter, the at least one further non-trivial convolutional filter comprising three or more non-zero filter coefficients and configured to, when convolved with the observed seismic data set, maximize the similarity or minimize the mismatch between the convolved observed data set and the convolved predicted data set.

30. The method according to claim 29, wherein the or more non-trivial convolutional filters is operable to transform at least a portion of said predicted seismic data set and the at least one further non-trivial convolutional filter is operable to transform said observed data set to render the observed seismic data set and predicted data set approximations of one another.

31. A non-transitory computer usable storage medium having a computer program product executable by a programmed or programmable processing apparatus, comprising one or more software portions for performing a method of subsurface exploration, the method comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement of at least one physical parameter, and comprising the steps of: providing an observed seismic data set comprising at least three distinct non-zero data values derived from at least three distinct non-zero seismic measured values of said portion of the volume of the Earth; generating, using a subsurface model of a portion of the Earth comprising a plurality of model coefficients, a predicted seismic data set comprising at least three distinct non-zero data values; utilizing one or more non-trivial convolutional filters, each convolutional filter comprising three or more non-zero filter coefficients and operable, when convolved with the observed seismic data set, to maximize the similarity or minimize the mismatch between the convolved observed data set and the predicted data set minimizing and/or maximizing one or more objective functions operable to measure the similarity and/or mismatch between the filter coefficients for the one or more non-trivial convolutional filters and at least three reference coefficients of one or more reference filters by modifying at least one model coefficient of said subsurface model of a portion of the Earth to produce an updated subsurface model of a portion of the Earth; and providing an updated subsurface model of a portion of the Earth for subsurface exploration.

32. A non-transitory computer usable storage medium having a computer program product executable by a programmed or programmable processing apparatus, comprising one or more software portions for performing a method of subsurface exploration, the method comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement of at least one physical parameter, and comprising the steps of: providing an observed seismic data set comprising at least three distinct non-zero data values derived from at least three distinct non-zero seismic measured values of said portion of the volume of the Earth; generating, using a subsurface model of a portion of the Earth comprising a plurality of model coefficients, a predicted seismic data set comprising at least three distinct non-zero data values; utilizing one or more non-trivial convolutional filters, each convolutional filter comprising three or more non-zero filter coefficients and operable, when convolved with the predicted data set, to maximize the similarity or minimize the mismatch between the convolved predicted data set and the observed seismic data set; minimizing and/or maximizing one or more objective functions operable to measure the similarity and/or mismatch between the filter coefficients for the one or more non-trivial convolutional filters and at least three reference coefficients of one or more reference filters by modifying at least one model coefficient of said subsurface model of a portion of the Earth to produce an updated subsurface model of a portion of the Earth; and providing an updated subsurface model of a portion of the Earth for subsurface exploration.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Embodiments of the present disclosure will now be described in detail with reference to the accompanying drawings, in which:

(2) FIG. 1 is a schematic illustration of a typical seismic survey experiment in which seismic traces are obtained from an undersea portion of the Earth;

(3) FIG. 2 is a schematic illustration of a basic starting model for full waveform inversion modelling;

(4) FIG. 3 is a schematic illustration of modelled seismic trace data generated from the basic starting model of FIG. 2 for an individual seismic shot;

(5) FIG. 4 shows seismic trace data for observed and predicted seismic data sets as a function of time;

(6) FIG. 5 shows the misfit of the observed and predicted seismic data sets in FIG. 4 as a function of time for conventional FWI;

(7) FIG. 6 a) to f) show examples of starting models and the resulting final models obtained from conventional full waveform inversion modelling;

(8) FIG. 7 shows a method according to a first embodiment of the present disclosure;

(9) FIG. 8 shows a method according to a second embodiment of the present disclosure;

(10) FIG. 9 shows the misfit of the observed and predicted seismic data sets in FIG. 4 as a function of time using the method of the present disclosure;

(11) FIG. 10 shows a method according to a third embodiment of the present disclosure;

(12) FIG. 11 shows a trial model for testing the accuracy of FWI;

(13) FIG. 12 shows the convergence achieved for conventional FWI relative to said trial model;

(14) FIG. 13 shows the convergence achieved for the method of the present disclosure relative to said trial model;

(15) FIG. 14 shows a further, more complex trial model for testing the accuracy of FWI;

(16) FIG. 15 shows an inaccurate starting model with which to begin the inversion;

(17) FIG. 16 shows the convergence achieved for conventional FWI relative to the trial model of FIG. 14; and

(18) FIG. 17 shows the convergence achieved for the method of the present disclosure relative to said trial model of FIG. 14.

DETAILED DESCRIPTION

(19) The present disclosure provides an improved methodology for FWI which avoids the problem of cycle skipping. Therefore, more accurate and reliable convergence can be obtained using the present disclosure when compared to known arrangements.

(20) Conventional FWI seeks to minimize the data misfit between the modelled and observed data sets directly or, alternatively, to maximize the similarity between the modelled and observed data sets.

(21) The present disclosure presents an alternative method which, for the first time, reduces the likelihood of cycle skipping. The inversion problem is divided into two stages. The first stage involves the solution of a relatively straightforward inversion problem that can be solved without the detrimental effects of cycle skipped mis-convergence. The output from the first stage provides a parametric description of the main inversion problem to be solved and of the misfit function that is to be minimized.

(22) The second stage is then to find the sub-surface model that minimizes this new parametric misfit function. If the parameterization is chosen appropriately, then cycle skipping in the major inversion step at the second stage is mitigated or eliminated.

(23) The following embodiments illustrate the application of the present disclosure in practice. The first embodiment outlines the general approach of the present disclosure. The second and third embodiments detail a specific implementation of the present disclosure utilizing Wiener filters.

(24) The following embodiments are described with regard to time domain analysis. However, the described methods are equally applicable to other domains, for example, frequency domain analysis. This may be achieved by performing a Fourier analysis on the time domain data and extracting particular frequencies for analysis. These aspects are to be considered to form part of the present disclosure and the skilled person would be readily aware of implementations of this.

(25) A method according to the present disclosure will now be described with reference to FIG. 7. FIG. 7 shows a flow diagram of a first embodiment of the present disclosure.

(26) Step 100: Obtain Observed Seismic Data Set

(27) Initially, it is necessary to obtain a set of experimentally gathered seismic data in order to initiate subsurface exploration. This may be gathered by an experimental arrangement such as the set up shown and described with reference to FIG. 1.

(28) The gathered seismic data may be optionally pre-processed in various ways including by propagating numerically to regions of the surface or subsurface where experimental data have not been acquired directly. The skilled person would readily be able to design and undertake such pre-processing as might be necessary or desirable. With or without such pre-processing, the resultant seismic dataset representing experimentally-gathered data is known as an observed seismic data set.

(29) As shown in FIG. 1, a large number of receivers or detectors 16 are positioned at well-known positions on the surface of the portion of the Earth to be explored. The detectors 16 may be arranged in a two dimensional (such as a line) or a three dimensional (such as a grid or plurality of lines) arrangement. The physical location of the detectors 16 is known from, for example, location tracking devices such as GPS devices. Additionally, the location of the source 12 is also well known by similar location tracking means.

(30) The observed seismic data set may comprise multiple source 12 emissions known in the art as shots. The data comprises pressure as a function of receiver position (on the x-axis) with respect to time (on the y-axis). This is because, in general, a detector such as a hydrophone measures the scalar pressure at its location. However, other arrangements may be used.

(31) The seismic trace data comprises a plurality of observed data points. Each measured discrete data point has a minimum of seven associated location valuesthree spatial dimensions (x, y and z) for receiver (or detector) position (r), three spatial dimensions (x, y, z) for source location (s), and one temporal dimension measuring the time of observation relative to the time of source initiation, together with pressure magnitude data. The seven coordinates for each discrete data point define its location in space and time.

(32) The seismic trace data also comprises one or more measurement parameters which denote the physical property being measured. In this embodiment, a single measurement parameter, pressure is measured. The observed data set is defined as d.sub.obs(r,s,t) and, in this embodiment, is in the time domain. For clarity, the following discussion considers a single source-receiver pair and so r, s are not needed.

(33) However, it is possible to measure other parameters using appropriate technology; for example, to measure the velocity or displacements of particles in three spatial dimensions in addition to the pressure. The present disclosure is applicable to the measurement of these additional variables.

(34) The actual gathering of the seismic data set is described here for clarity. However, this is not to be taken as limiting and the gathering of the data may or may not form part of the present disclosure. The present disclosure simply requires a real-world observed seismic data set upon which analysis can be performed to facilitate subsurface exploration of a portion of the Earth. The method now proceeds to step 102.

(35) Step 102: Provide Starting Model

(36) At step 102, an initial starting model of the specified subsurface portion of the Earth is provided. The model may be provided in either a two dimensional or a three dimensional form. Whilst the illustrated examples are of a two-dimensional form, the skilled person would be readily aware that the present disclosure is applicable to three dimensional approaches.

(37) The generated model consists of values of the coefficient V.sub.p and, possibly, other physical values or coefficients, typically defined over a discrete grid representing the subsurface. Such starting models are routinely generated and represent the general trends of the major features within the subsurface region to be modelled and could be readily generated by the skilled person. As discussed above, it is normally necessary to provide a starting model of high accuracy to ensure non-cycle skipped convergence. However, in the case of the present disclosure a less accurate initial model could be used. This may reduce the required time and resources to generate the starting model whilst still enabling accurate convergence.

(38) Predicted seismic data may be generated based on an analysis of the acoustic isotropic two-way wave equation as set out below in equation 17):

(39) $\begin{array}{cc}\frac{1}{c^2}\frac{{}^{2}p}{t^2}-.Math.\left(\frac{1}{}p\right)=s& 17)\end{array}$
where the acoustic pressure p and driving source s vary in both space and time, and the acoustic velocity c and density vary in space. This equation applies to small-amplitude pressure waves propagating within an inhomogeneous, isotropic, non-attenuating, non-dispersive, stationary, fluid medium. It is relatively straightforward to add elastic effects, attenuation and anisotropy to the wave equation. Introduction of these parameters changes the detailed equations and numerical complexity, but not the general approach.

(40) The wave equation 17) represents a linear relationship between a wave field p and the source s that generates the wave field. After discretization (with, for example, finite differences) we can therefore write equation 17) as a matrix equation 18):
Ap=s18)
where p and s are column vectors that represent the source and wavefield at discrete points in space and time, and A is a matrix that represents the discrete numerical implementation of the operator set out in equation 3):

(41) $\begin{array}{cc}\frac{1}{c^2}\frac{{}^2}{t^2}-.Math.\left(\frac{1}{}\right)& 19)\end{array}$

(42) Although the wave equation represents a linear relationship between p and s, it also represents a non-linear relationship between a model m and wavefield p. Thus equation 17) can be rewritten as equation 20):
G(m)=p20)
where m is a column vector that contains the model parameters. Commonly these will be the values of c (and if density is an independent parameter) at every point in the model, but they may be any set of parameters that is sufficient to describe the model, for example slowness 1/c, acoustic modulus c.sup.2 p, or impedance cp.

(43) In equation 20), G is not a matrix. Instead it is a non-linear Green's function that describes how to calculate a wavefield p given a model m. Once the model has been generated, the method then proceeds to step 104.

(44) Step 104: Generate Predicted Data Set

(45) At step 104, a predicted seismic data set is generated. The predicted data is required to correspond to the same source-receiver location data positions as the actual measured trace data so that the modelled and observed data can be compared. In other words, the predicted data set corresponds discrete point to discrete point to the observed dataset. The predicted data set is generated for the same measurement parameter(s) at the same frequency or frequencies.

(46) From the above analysis, predicted seismic data can be generated for one or more physical parameters in the time domain. If done in the frequency domain it could be done for one or more selected frequencies. This forms the predicted seismic data set d.sub.pred. The method now proceeds to step 106.

(47) Step 106: Design a Convolutional Filter

(48) In step 106, a convolutional filter, or filters is designed. The aim of the convolutional filter is to take as an input all or part of the predicted data d.sub.pred(r,s), and from that generate as an output an approximation to all or part of the observed data d.sub.obs(r,s). A filter is designed for each source-receiver pair r,s to generate a set of coefficients specific to particular r,s values. More than one filter may be designed in this step if required.

(49) It is to be understood that the term convolutional filter in the present disclosure may have broad applicability. For example, the convolution may be non-causal, it may be applied in one or more dimensions, including spatial, temporal or other dimensions, the number of dimensions and/or number of data points on input need not the same as the number of dimensions or data points on output, the filter may vary in space, in time and/or in other dimensions, and it may involve post-processing following convolution.

(50) The convolutional filter may comprise any generalized convolutional operation that can be described by a finite set of parameters that depend upon both the predicted d.sub.pred(r,s) and observed d.sub.obs(r,s) data, such that when the convolution and associated operations are applied to all or part of the predicted data d.sub.pred(r,s) an accurate or generally approximate model of the observed data is generated.

(51) The skilled person would be readily aware of how to design and apply such convolutional filters in order to match together two datasets. For example, an objective function could be generated and this could be maximized or minimized to achieve the necessary matching. The method proceeds to step 108.

(52) Step 108: Generate Reference Filter

(53) Once the convolutional filter is designed in step 106, an analogous filter is designed such that, when it is applied to an input dataset, the filter will generate an output dataset that provides an equivalent approximation to all or parts of the input dataset. In other words, this filter does not significantly modify the input data other than optionally in simple, predetermined and well-defined ways, for example by scaling in amplitude.

(54) In one embodiment, the reference filter comprises a sequence of zeros and a non-zero value at the zero lag time. This is, in essence, similar to an impulse function centered upon time t=0. This corresponds to an identity filter in that it does not transform the form of the data, although the magnitude of the output data may not necessarily be the same as the input data.

(55) However, other types of reference filter will be apparent to the skilled person. Given a particular form for a convolutional filter, the skilled person would readily know how to design such an identity filter. The method then proceeds to step 110.

(56) Step 110: Construct Misfit Function

(57) At step 110, a misfit (or objective) function is configured. In one example, the misfit function (or objective function) is configured to measure the dis-similarity between the actual filter coefficients and reference filter coefficients. Alternatively, an objective function may be configured that measures similarity; in this case, step 112 will operable to maximize rather than minimize the objective function.

(58) An example of an objective function configured to measure the dis-similarity between a simple one-dimensional convolutional filter in time and a reference function that consist of only one non-zero value at zero lag, would be to weight the convolution filter coefficients by the modulus of their temporal lag. The objective function would then consist of some norm of these weighted coefficients divided by the same norm of the unweighted coefficients. If the L.sub.2 norm is used here, then this objective function will provide the least-squares solution, but other norms (for example, the L.sub.1 norm) are potentially utilizable.

(59) The norm of the weighted coefficients must be normalized by the norm of the unweighted coefficients in this example otherwise the objective function could be simply minimized by driving the predicted data to large values, and hence driving the filter coefficients to small values.

(60) In this example, the coefficients generated for each source receiver pair r,s in step 106 are weighted as a function of the modulus of the temporal lag. In other words, the coefficients are weighted based on the data position in time for a time domain analysis.

(61) However, it is to be understood that other types of weighting could be used. For example, more complicated functions of the temporal lag are possible such as weighting with a Gaussian function of lag centered on zero lag.

(62) In general two types of weighting are desirable; those that increase monotonically away from zero lag, such as the modulus, and those that decrease monotonically away from zero lag, such as a Gaussian weighting. The former type of weighting will lead to objective functions that should be minimized and the latter type will lead to objective functions that should be maximized. Combinations of these two types are also possible. The skilled person would readily understand how to design such objective functions and how to minimize or maximize them. The method then proceeds to step 112.

(63) Step 112: Minimize or Maximize the Misfit Function

(64) This may be done by any suitable method. In this embodiment, the gradient method is used.

(65) However, in contrast to conventional FWI, the present disclosure seeks to minimize the misfit between the convolutional and reference filters, rather than between the data sets themselves. Therefore, because the convolutional filters are not periodic, cycle skipping does not occur.

(66) In one embodiment, the gradient of the misfit functional is obtained in accordance with conventional FWI, albeit minimizing the filter coefficients. The method then proceeds to step 114.

(67) Step 114: Update Model

(68) At step 114, the model is updated using the gradient obtained in step 116. The model update derives from the gradient of equation 11) i.e. the partial derivative with respect to a point perturbation of the model m at each position. Ultimately gradients from separate shots will summed when forming the final model update.

(69) As with the computational structure of conventional FWI method, this is the result of two wavefields: an incident wavefield emitted by a source at the source location, and a back-propagated wavefield which is emitted by a (multi-point) source located the receiver positions.

(70) As noted above, gradient methods can be enhanced using an approximate form for the Hessian and conjugate directions for the model update. Further a step-length is then calculated to scale the search direction and give the final model update.

(71) For any useful convolutional filter, and for any useful measure of misfit or similarity, as the convolutional filter moves towards the reference filter, the predicted seismic data set will move towards the observed seismic data set. Thus, the starting model will move towards the true model.

(72) With an appropriate choice of convolutional filter and measure of misfit/similarity, this scheme is unaffected or less affected by cycle skipping than is conventional FWI. It can therefore be applied more readily to observed data that lack low frequencies and/or in the absence of a good starting model. The method now proceeds to step 116.

(73) Step 116: Convergence Criteria Met?

(74) At step 116 it is determined whether convergence criteria have been met. For example, when the method is deemed to have reached convergence when the difference between the data sets reaches a threshold percentage or other value. If the criteria as set out above have been met, then the method proceeds to step 118 and finishes with the resultant Earth model generated. If the criteria have not been met, then the method proceeds back to repeat steps 104 to 114 as discussed above.

(75) Step 118: Finish

(76) When, at step 118, it has been determined that the convergence criteria has been met, the method finishes and the modelled subsurface portion of the Earth is deemed to be sufficiently accurate to be used for subsurface exploration. This may involve the direct interpretation of the recovered model, and/or involve the process of depth-migration to generate a subsurface reflectivity image to be used for the identification of subsurface features such as cavities or channels which may contain natural resources such as hydrocarbons. Examples of such hydrocarbons are oil and natural gas.

(77) Once these features have been identified in the subsurface model and/or reflectivity image, then measures can be taken to investigate these resources. For example, survey vessels or vehicles may be dispatched to drill pilot holes to determine whether natural resources are indeed present in these areas.

(78) The method according to the present disclosure has numerous advantages over prior art methods. In experimental tests, the method has proven to be robust against cycle skipping with starting models that have errors that are several times larger than those that are fatal to conventional FWI.

(79) It is also possible to use the above method in conjunction with conventional FWI. For example, the above method may be used initially to refine the starting model. Once the starting model is sufficiently accurate, conventional FWI could be used.

(80) A second embodiment of the disclosure is illustrated in FIG. 8. The second embodiment focuses on a specific application of Wiener filters to the present disclosure.

(81) Step 200: Obtain Observed Seismic Data Set

(82) Step 200 corresponds substantially to method step 100 of the previous embodiment. Therefore, this will not be described again here. Only the steps which are new to this embodiment of the method of the present disclosure will be described. The method now proceeds to step 202.

(83) Step 202: Provide Starting Model

(84) At step 202, an initial starting model of the specified subsurface portion of the Earth is provided. The model may be provided in either a two dimensional or a three dimensional form. Whilst the illustrated examples are of two-dimensional form, the skilled person would be readily aware that the present disclosure is applicable to three dimensional approaches. The model is generated in this step in accordance with step 102 above.

(85) The generated model consists of values of the coefficient V.sub.p and, possibly, other physical values or coefficients over a discrete grid representing the subsurface. Such starting models are routinely generated and represent the general trends of the major features within the subsurface region to be modelled and could be readily generated by the skilled person. The method then proceeds to step 204.

(86) Step 204: Generate Predicted Data Set

(87) In order to model accurately the subsurface region under investigation, it is necessary to generate predicted data which corresponds to the same source-receiver location data position as the actual measured trace data so that the modelled and observed data can be compared. In other words, the predicted data set corresponds discrete point to discrete point to the observed dataset. The predicted data set is generated for the same measurement parameter(s) at the same frequency or frequencies.

(88) The predicted data set is generated using the full two-way wave equation. In one approach, the modelled trace data may be generated using the time domain full two-way wave equation as set out above.

(89) From the above analysis, predicted seismic data can be generated for one or more physical parameters and at one or more selected frequencies. This forms the predicted seismic data set d.sub.pred(r,s). The method now proceeds to step 206.

(90) Step 206: Scale Data

(91) At step 206, the data is scaled so that the predicted d.sub.pred(r,s) and observed d.sub.obs(r,s) data sets comprise matching root mean square amplitudes. This may be accomplished by any known approach. Either or both of the predicted d.sub.pred(r,s) and observed d.sub.obs(r,s) data sets may be scaled as appropriate.

(92) Additionally, this step may be optional if it is possible to generate predicted data which intrinsically matches the observed data. The method proceeds to step 208.

(93) Step 208: Design Wiener Filter

(94) A Wiener filter is a convolutional filter of finite length that is operable to convert an input wavelet into a desired output wavelet in a least-squares manner. In other words, if b is the input wavelet, and c is the desired output wavelet, then the Wiener filter w minimizes expression 21).
b*wc.sup.221)
Where expression 21) represents, inter alia, the convolution of input wavelet b with filter w. However, any convolution can be expressed as a matrix operation, and so b*w=c can be expressed as equation 22):
Bw=c 22)

(95) If w is a column vector of length M, and c is a column vector of length M+N1, then B is a rectangular matrix with M columns and M+N1 rows that contains the elements of b arranged as:

(96) $\begin{array}{cc}B=\left[\begin{array}{ccccc}{b}_0& 0& \mathrm{.Math.}& \mathrm{.Math.}& 0\\ b_1& b_0& & & 0\\ b_2& b_1& & & \mathrm{.Math.}\\ b_3& b_2& & & 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & & b_0\\ b_{N-1}& b_{N-2}& & & b_1\\ 0& b_{N-1}& & & \mathrm{.Math.}\\ 0& 0& & & \mathrm{.Math.}\\ \mathrm{.Math.}& \mathrm{.Math.}& & & b_{N-2}\\ 0& 0& \mathrm{.Math.}& \mathrm{.Math.}& b_{N-1}\end{array}\right]& 23)\end{array}$
Based on equation 23), B.sup.TB represents a matrix containing the autocorrelation of b in each column with the zero lag on the main diagonal, and B.sup.Tc is the cross-correlation of b and c. The solution to the least squares problem in expression 21) is now:
w=(B.sup.TB).sup.1B.sup.Tc24)

(97) Therefore, the approach is to find the auto-correlation B.sup.TB of the input trace b, the cross-correlation B.sup.Tc of the input trace with the desired output trace c, and deconvolve the latter using the former. Because of the Toeplitz structure of B, fast algorithms (such as the Levinson algorithm) to find the filter coefficients if B is square which can be ensured by suitably padding b and c with zeros. In practice, equation 25) would be solved:
w=(B.sup.TB+I).sup.1B.sup.Tc25)

(98) This approach can be applied to FWI. At each receiver, for each source, a Wiener filter w can be designed that operates on the observed seismic data set dobs to convert the observed seismic data set d.sub.obs into the predicted data set d.sub.pred. In this embodiment, the filter w operates to convert the observed data into the predicted data. The third embodiment discloses an alternative approach. Equation 26) and 27) denote the approach of this embodiment.
d.sub.obs*w=d.sub.pred26)
w=(B.sup.TB).sup.1B.sup.Td.sub.pred27)
Here B is the matrix representation of convolution by d.sub.obs, and B.sup.T is the matrix representation of cross-correlation with d.sub.obs. The method proceeds to step 210.

(99) Step 210: Generate Reference Filter

(100) Once the convolutional filter is designed in step 208, an analogous convolutional filter is designed such that, when applied to an input dataset, this reference filter will generate an output dataset that provides an equivalent approximation to all or parts of the input dataset. In other words, the reference filter does not modify the input data in essence.

(101) In one embodiment, the reference filter comprises a sequence of zeros and a value of one at the zero lag time. This is, in essence, similar to an impulse function centered upon time t=0. The method then proceeds to step 212.

(102) Step 212: Construct Misfit Function

(103) At step 212, a misfit function is configured. The misfit function (or objective function) is configured to measure the dis-similarity between the actual filter coefficients and reference filter coefficients. Alternatively, an objective function may be configured that measures similarity; in this case, step 214 (described later) will operable to maximize rather than minimize the objective function.

(104) An example of an objective function configured to measure the dis-similarity between a simple one-dimensional convolutional filter in time and a reference function that consist of only one non-zero value at zero lag, would be to weight the convolution filter coefficients by the modulus of their temporal lag. The objective function would then consist of some norm of these weighted coefficients divided by the same norm of the unweighted coefficients. If the L.sub.2 norm is used here, then this objective function will provide the least-squares solution, but other norms are possible and potentially desirable.

(105) The norm of the weighted coefficients must be normalized by the norm of the unweighted coefficients in this example otherwise the objective function could be simply minimized by driving the predicted data to large values, and hence driving the filter coefficients to small values.

(106) In this example, the coefficients generated for each source receiver pair r,s in step 208 are weighted as a function of the modulus of the temporal lag. In other words, the coefficients are weighted based on the data position in time for a time domain analysis.

(107) However, it is to be understood that other types of weighting could be used. For example, more complicated functions of the temporal lag are possible such as weighting with a Gaussian function of lag centered on zero lag. In general two types of weighting are desirable; those that increase monotonically away from zero lag, such as the modulus, and those that decrease monotonically away from zero lag, such as a Gaussian weighting. The former type of weighting will lead to objective functions that should be minimized and the latter type will lead to objective functions that should be maximized. Combinations of these two types are also possible. The skilled person would readily understand how to design such objective functions and how to minimize or maximize them.

(108) A model is then sought that makes the Wiener filters as close as possible to being just a spike at zero time lag. In other words, the predicted and modelled data matches apart from a scale factor. Therefore, a model m is desired that minimizes or maximizes:

(109) $\begin{array}{cc}f=\frac{\frac{1}{2}{.Math.\mathrm{Tw}.Math.}^2}{{.Math.w.Math.}^2}& 28)\end{array}$
where w is a column vector containing all the Wiener filter coefficients, and T is the weighting function, for example the modulus of temporal lag. The filters vary with source and receiver position, and with the model. The Wiener filters are obtained by solving equation (27). The method proceeds to step 214.

(110) Step 214: Minimize Misfit Function

(111) This may be done by any suitable method. In this embodiment, the gradient method is used. In this embodiment, the gradient is straightforward to derive from equation 28):

(112) $\begin{array}{cc}\frac{f}{m}=-{p^T\left(\frac{A}{m}\right)}^T{A}^{-T}\left\{{B\left(B^{T}B\right)}^{-1}\left(\frac{T^2-2\mathrm{fI}}{w^{T}w}\right)w\right\}& 29)\end{array}$
This is the usual expression for FWI, but now the adjoint source that must be used in the back propagation is no longer simply the residual data but is instead the entity represented by:

(113) $\left\{{B\left(B^{T}B\right)}^{-1}\left(\frac{T^2-2jT}{w^{T}w}\right)w\right\}.$

(114) As set out above, in contrast to conventional FWI, the present disclosure seeks to minimize the difference between the convolutional and reference filters, rather than between the data sets themselves. Therefore, because the convolutional filters are not periodic, cycle skipping does not occur.

(115) The gradient can be obtained based on equation 29) as follows. A set of Wiener filters is found in step 208, one filter per data trace. The filter coefficients are normalized by their inner product, trace-by-trace. The normalized coefficients are weighted by a function of temporal lag. The resultant sequence is deconvolved by the auto-correlation of the observed data, and convolved with the observed data. This forms an adjoint source for each source-receiver pair. These adjoint sources are then back-propagated and combined in the normal way with the forward wavefield to produce the gradient. The method then proceeds to step 216.

(116) Step 216: Update Model

(117) At step 216, the model is updated using the gradient obtained in step 214. The model update derives from the gradient of equation 11) i.e. the partial derivative with respect to a point perturbation of the model m at each position. Ultimately gradients from separate shots will summed when forming the final model update.

(118) As with the computational structure of conventional FWI method, this is the product of two wavefields: an incident wavefield emitted by a source at source location back-propagated wavefield which is emitted by a (multi-point) source located the receiver positions.

(119) As noted above, gradient methods can be enhanced using an approximate form for the Hessian and conjugate directions for the model update. Further a step-length is then calculated to scale the search direction and give the final model update.

(120) For any useful convolutional filter, and for any useful measure of difference or similarity, as the convolutional filter moves towards the reference filter, the predicted seismic data set will move towards the observed seismic data set. Thus, the starting model will move towards the true model.

(121) With an appropriate choice of convolutional filter and measure of difference/similarity, this scheme is unaffected or less affected by cycle skipping than is conventional FWI. It can therefore be applied more readily to observed data that lack low frequencies and/or in the absence of a good starting model. The method now proceeds to step 218.

(122) Step 218: Convergence Criteria Met?

(123) At step 218 it is determined whether convergence criteria have been met. For example, when the method is deemed to have reached convergence when the difference between the data sets reaches a threshold percentage. If the criteria as set out above have been met, then the method proceeds to step 220 and finishes with the resultant Earth model generated. If the criteria have not been met, then the method proceeds back to repeat steps 204 to 216 as discussed above.

(124) Step 220: Finish

(125) When, at step 220, it has been determined that the convergence criteria has been met, the method finishes and the modelled subsurface portion of the Earth is deemed to be sufficiently accurate to be used for subsurface exploration. This may involve the direct interpretation of the recovered model, and/or involve the process of depth-migration to generate a subsurface reflectivity image to be used for the identification of subsurface features such as cavities or channels which may contain natural resources such as hydrocarbons. Examples of such hydrocarbons are oil and natural gas.

(126) Once these features have been identified in the subsurface model and/or reflectivity image, then measures can be taken to investigate these resources. For example, survey vessels or vehicles may be dispatched to drill pilot holes to determine whether natural resources are indeed present in these areas.

(127) FIG. 9 shows an example of a misfit function as a function of time, similar to FIG. 5 for conventional FWI. As can be seen, the misfit is not cyclic and so cycle skipping is eliminated.

(128) A third embodiment of the disclosure is illustrated in FIG. 10. The third embodiment focuses on an alternative application of Wiener filters to the present disclosure.

(129) Step 300: Obtain Observed Seismic Data Set

(130) Step 300 corresponds substantially to method steps 100 and 200 of the previous embodiments. Therefore, this will not be described again here. Only the steps which are new to this embodiment of the method of the present disclosure will be described. The method then proceeds to step 302.

(131) Step 302: Provide Starting Model

(132) At step 302, an initial starting model of the specified subsurface portion of the Earth is provided. The model may be provided in either a two dimensional or a three dimensional form. Whilst the illustrated examples are of two dimensional form, the skilled person would be readily aware that the present disclosure is applicable to three dimensional approaches. The model is generated in this step in accordance with step 102 above.

(133) The generated model consists of values of the coefficient V.sub.p and, possibly, other physical values or coefficients over a discrete grid representing the subsurface. Such starting models are routinely generated and represent the general trends of the major features within the subsurface region to be modelled and could be readily generated by the skilled person. The method then proceeds to step 304.

(134) Step 304: Generate Predicted Data Set

(135) In order to model accurately the subsurface region under investigation, it is necessary to generate predicted data which corresponds to the same source-receiver location data position as the actual measured trace data so that the modelled and observed data can be compared. In other words, the predicted data set corresponds discrete point to discrete point to the observed dataset. The predicted data set is generated for the same measurement parameter(s) at the same frequency or frequencies.

(136) The predicted data set is generated using the full two-way wave equation. In one approach, the modelled trace data may be generated using the time domain full two-way wave equation as set out above.

(137) From the above analysis, predicted seismic data can be generated for one or more physical parameters and at one or more selected frequencies. This forms the predicted seismic data set d.sub.pred(r,s). The method now proceeds to step 306.

(138) Step 306: Scale Data

(139) At step 306, the data is scaled so that the predicted d.sub.pred(r,s) and observed d.sub.obs(r,s) data sets comprise matching root mean square amplitudes. This may be accomplished by any known approach. Either or both of the predicted d.sub.pred(r,s) and observed d.sub.obs(r,s) data sets may be scaled as appropriate.

(140) Additionally, this step may be optional if it is possible to generate predicted data which intrinsically matches the observed data. The method proceeds to step 308.

(141) Step 308: Design Wiener Filter

(142) In step 308, the Wiener filter is designed. In this embodiment, it is desired to design a filter w which converts the predicted data set d.sub.pred into the observed data d.sub.obs. In other words, equations 30) and 31) set out these parameters:
d.sub.pred*x=d.sub.obs30)
x=(C.sup.TC).sup.1C.sup.Td.sub.obs31)
where C is the matrix representation of convolution by d.sub.pred. A Wiener filter can be designed on this basis. The method proceeds to step 310.

(143) Step 310: Weight Coefficients

(144) At step 310, the coefficients generated for each source receiver pair r,s in step 308 are weighted as a function of the modulus of the temporal lag. As set out in equation 32), T is a diagonal matrix arranged to scale the Wiener filter coefficients as a function of the modulus of the time lag.

(145) Therefore, if the value of this function increases as the lag increases, then it is desired to minimize functional g. Conversely, if the function decreases as the lag increases, then g should be maximized. The former is more likely to provide better resolution, whilst the latter approach is likely to be more resistant to noise.

(146) Step 312: Generate Reference Filter

(147) Once the convolutional filter is designed in step 308, an analogous convolutional filter is designed such that, when applied to an input dataset, the filter will generate an output dataset that provides an equivalent approximation to all or parts of the input dataset.

(148) In other words, this filter does not modify the input data in essence and so is denoted the reference filter, which comprises the auto parameter set.

(149) In one embodiment, the reference filter comprises a sequence of zeros and a value of one at the zero lag time (i.e. at the receiver). This is, in essence, similar to an impulse function centered upon time t=0. The method then proceeds to step 314.

(150) Step 314: Construct Misfit Function

(151) At step 314, a misfit function is configured. The misfit function (or objective function) is configured to measure the difference between the auto and cross parameter sets. Alternatively, the misfit function may be operable to measure similarity and seek to maximize it.

(152) A model is then sought that makes the Wiener filters as close as possible to being just a spike at zero time lag. In other words, the predicted and modelled data matches apart from a scale factor. Therefore, a model m is desired that minimizes or maximizes:

(153) $\begin{array}{cc}g=\frac{\frac{1}{2}{.Math.\mathrm{Tx}.Math.}^2}{{.Math.x.Math.}^2}& 32)\end{array}$
where x is a column vector containing all the Wiener filter coefficients, and T is a weighting function applied in step 210. The filters vary with source and receiver position, and with the model. The Wiener filters are obtained by solving equation (31). Equation 32) is a normalised misfit (or objective) function due the normalisation by a norm of the column vector x. This eliminates, or reduces, the sensitivity of the function to the scale of the reference filter coefficients. The method proceeds to step 316.

(154) Step 316: Minimize Misfit Function

(155) This may be done by any suitable method. In this embodiment, the gradient method is used. In this embodiment, the gradient is straightforward to derive from equation 32):

(156) $\begin{array}{cc}\frac{g}{m}=-{p^T\left(\frac{A}{m}\right)}^T{A}^{-T}\left\{X^T{C\left(C^{T}C\right)}^{-1}\left(\frac{2\mathrm{gI}-T^2}{x^{T}x}\right)x\right\}& 33)\end{array}$
where X.sup.T represents cross-correlation with the Wiener filter generated in step 308.

(157) Therefore, in summary, to find the gradient in this approach, the Wiener filter is found in step 308, it is normalized, weighted as a function of lag, deconvolved using the auto-correlation of the predicted data, convolved with the predicted data, and cross-correlated with the Wiener filter. This forms an adjoint source for each source-receiver pair. These adjoint sources are then back-propagated as normal, and combined with the forward wavefield to generate the gradient. The method then proceeds to step 318.

(158) Step 318: Update Model

(159) At step 318, the model is updated using the gradient obtained in step 316. The model update derives from the gradient of equation 33) i.e. the partial derivative with respect to a point perturbation of the model m at each position. Ultimately gradients from separate shots will summed when forming the final model update.

(160) For any useful convolutional filter, and for any useful measure of difference or similarity, as the convolutional filter moves towards the reference filter, the predicted seismic data set will move towards the observed seismic data set. Thus, the starting model will move towards the true model.

(161) With an appropriate choice of convolutional filter and measure of difference/similarity, this scheme is unaffected or less affected by cycle skipping than is conventional FWI. It can therefore be applied more readily to observed data that lack low frequencies and/or in the absence of a good starting model. The method now proceeds to step 320.

(162) Step 320: Convergence Criteria Met?

(163) At step 320 it is determined whether convergence criteria have been met. For example, when the method is deemed to have reached convergence when the difference between the data sets reaches a threshold percentage. If the criteria as set out above have been met, then the method proceeds to step 322 and finishes with the resultant Earth model generated. If the criteria have not been met, then the method proceeds back to repeat steps 304 to 318 as discussed above.

(164) Step 322: Finish

(165) When, at step 322, it has been determined that the convergence criteria has been met, the method finishes and the modelled subsurface portion of the Earth is deemed to be sufficiently accurate to be used for subsurface exploration. This may involve the direct interpretation of the recovered model, and/or involve the process of depth-migration to generate a subsurface reflectivity image to be used for the identification of subsurface features such as cavities or channels which may contain natural resources such as hydrocarbons. Examples of such hydrocarbons are oil and natural gas.

(166) Once these features have been identified in the subsurface model and/or reflectivity image, then measures can be taken to investigate these resources. For example, survey vessels or vehicles may be dispatched to drill pilot holes to determine whether natural resources are indeed present in these areas.

(167) The effectiveness of the above-described embodiments is illustrated with reference to FIGS. 11 to 13. FIG. 11 shows a relatively extreme subsurface model including a strong anomaly to resolve.

(168) FIG. 12 shows a model derived from conventional FWI. As shown, the model fails to converge correctly due to cycle skipping. Consequently, the minimized model bears little resemblance to the actual model as set out in FIG. 11.

(169) In contrast, FIG. 13 shows the model recovered using the method of the present disclosure. As shown, the model is fully recovered subject only to the limitations of geometry and bandwidth. Clearly, the present disclosure enables accurate modelling of complex subsurface features where conventional FWI fails.

(170) A further example of the effectiveness of the above-described embodiments is illustrated with reference to FIGS. 14 to 17.

(171) FIG. 14 shows a Marmousi model used to test the method of the present disclosure.

(172) FIG. 15 shows a simple and inaccurate starting model used to test the present disclosure.

(173) FIG. 16 shows a model derived from conventional FWI. The data used here had a dominant frequency of 10 Hz. As shown, the resulting model bears little resemblance to the correct model, which is significantly cycle skipped with respect to the inaccurate starting model. In other words, the model fails to converge correctly due to cycle skipping. Consequently, the minimized model bears little resemblance to the actual model as set out in FIG. 14.

(174) In contrast, FIG. 17 shows the model recovered using the method of the present disclosure. As shown, the model is fully recovered subject only to the limitations of geometry and bandwidth. Clearly, the present disclosure enables accurate modelling of complex subsurface features where conventional FWI fails.

(175) Variations of the above embodiments will be apparent to the skilled person. The precise configuration of hardware and software components may differ and still fall within the scope of the present disclosure.

(176) Firstly, whilst the above embodiments are illustrated with regard to one dimensional Wiener filters, multidimensional filters could also be used. In the simple one-dimensional scheme above, the filter for each source-receiver pair is designed only using data from that source-receiver pair. The scheme is also implemented in timethe two sequences that are being matched represent data that varies in time, and the Wiener filter has its coefficients arranged by temporal lag.

(177) However, it is possible to implement an analogous scheme in multiple dimensions in space, or in a mixture of time and one or more space dimensions, or in any of many transformed spaces, for example in intercept-slowness (tau-p) space, in frequency-wavenumber (f-k) space, or in parabolic Radon space.

(178) Considering a single source and a single receiver, then, using a source estimate, the predicted data can be calculated from that source not only at the receiver also in a three dimensional volume around it. For that receiver, a one-dimensional observed dataset and a four-dimensional predicted dataset is now available. That is, three space dimensions plus time. Now a convolutional filter can be designed that takes all or some 1D, 2D or 3D subset of this 4D predicted dataset as input and generates the 1D observed dataset as output. This can be done in any of steps 106, 208 or 308.

(179) If this convolutional filter is a multi-dimensional Wiener convolutional filter, then the corresponding reference filter will be unity at the coefficient that corresponds to zero lag in time and zero lag in all of the space dimensionsthat is at the position of the receiverand it will be zero everywhere else.

(180) In this context, the method remains unchanged but now a 1D, 2D, 3D or 4D convolutional filter can be designed, and this can be combined with an appropriate function of all the temporal and spatial lags involved. The associated misfit functional can then be minimized.

(181) In practice, it is unlikely to be necessary to use all four dimensions together because the computational effort would be large. However, different subsets of dimensions on different iterations may be required.

(182) Further, it is possible to generate data at or around the receiver, for multiple source locations, and for different source initiation times. It takes four dimensions to describe the source time and location. Therefore, up to eight dimensions of predicted data are available. A filter can be designed using any combination of these as input though in a practical scheme we would be unlikely to want to use all of them simultaneously. There is also redundancy between the two time dimensions unless at least one space dimension is also involved in designing the filter.

(183) The dimensionality of the output can also be increased. Thus, rather than design a filter that seeks to match the observed data in just one dimension, the observed data could be matched in two or more dimensions, for example with the observed data predicted in time over a two-dimensional array of surface receivers, giving a 3D output from the convolutional filter.

(184) Equivalent schemes may also be generated that match the observed data to the predicted data. In this case, we would use observed data from multiple receivers and/or multiple sources as input, generating predicted data over the same multi-dimensional volume or some subset of it.

(185) These multi-dimensional filters could be extended by applying them to virtual sources and receivers located anywhere within the sub-surface including at the source position.

(186) Further, whilst the above embodiments have been illustrated in the time domain, frequency domain analysis could also be used. The different frequency components can be extracted through the use of a Fourier transform. Fourier transforms and associated coefficients are well known in the art. Since only a few frequency components may be required, a discrete Fourier transform may be computed from the inputted seismic trace data. Alternatively, a conventional fast Fourier transform (FFT) approach may be used. As a further alternative, partial or full Fourier transforms may be used. The Fourier transform may be taken with respect to complex or real valued frequencies.

(187) One or more observed data sets at one or more frequencies may be extracted. For example, a plurality of observed data sets may be provided having frequency components of, for example, 4.5 Hz, 5 Hz and 5.5 Hz. These 3 frequencies can be inverted individually with the output of one becoming the input of the next.

(188) As can be understood, further variations to this method are possible. For example, the process could be applied using the reverse filter, that using a filter that takes the observed data as input and that generates an approximation to the predicted data as output;

(189) The data could be pre-processed in any desired manner the observed and/or field data prior to the formation of the Wiener filter coefficients.

(190) The auto and/or cross correlation functions could be pre-processed in order to form the Wiener filter.

(191) The Wiener filter coefficients and/or the weighted Wiener coefficients could be processed before forming the quantity to be minimized or maximized. Alternatively, post-processing could be performed on the quantity to be minimized or maximized.

(192) The scheme could be applied to data, to models, and/or to filter coefficients that have been transformed into some other domain or domains including but not limited to the single or multi-dimensional Fourier domain, the Laplace domain, the Radon domain, the wavelet domain, the curvelet domain.

(193) The observed and/or predicted data could be transformed into a sequence of complex numbers, for example by forming an imaginary part using the Hilbert transform, such that the Wiener coefficients are themselves complex numbers, and optionally combine this with the minimization of the imaginary portion of the Wiener coefficients.

(194) Additional constraints could be added to the Wiener coefficients including but not limited to constraints generated by the observed data and/or the predicted data and/or the model and/or the other Wiener coefficients and/or the accuracy of the convolutional filter.

(195) Wiener filters could be used that vary continuously or discontinuously in time, in space, with receiver position, with source position, with source-receiver azimuth, with source-receiver offset, with source-receiver midpoint, and/or with other characteristics of the data, the model, the model updates, the source, the receiver, the filer coefficients.

(196) In other words, the Wiener filter need not be one dimensional and may encompass multiple dimensions.

(197) In addition, whilst embodiments have illustrated the use of Wiener filters, the present disclosure is not limited to such filters. The skilled person would readily be aware of other filter types suitable for use with the present disclosure. Non-limiting examples may include: Kalman filters; least mean squares filters; normalized least mean squares filters; and recursive least squares filters.

(198) Further, whilst the present disclosure has been described with reference to a model which involves solving the acoustic wave equation, the present disclosure is equally applicable to the models which involve solving the visco-acoustic, elastic, visco-elastic or poro-elastic wave equation.

(199) Further, whilst the example here has used the scalar parameter of pressure as its focus (i.e. P-waves), it is also possible (with appropriate equipment such as geophones) to measure the particle motion in the receiver location and so determine S-wave parameters. Data so-generated could then be utilized and processed in the present disclosure.

(200) Embodiments of the present disclosure have been described with particular reference to the examples illustrated. While specific examples are shown in the drawings and are herein described in detail, it should be understood, however, that the drawings and detailed description are not intended to limit the disclosure to the particular form disclosed. It will be appreciated that variations and modifications may be made to the examples described within the scope of the present disclosure.