WO2013093468A2

WO2013093468A2 - Full waveform inversion quality control method

Info

Publication number: WO2013093468A2
Application number: PCT/GB2012/053198
Authority: WO
Inventors: Nikhil Koolesh SHAH; Mike Warner; Adrian UMPLEBY; Lluis GUASH; Tenice Nangoo; Sainath Lakshminarayanan
Original assignee: Shah Nikhil Koolesh; Mike Warner; Umpleby Adrian; Guash Lluis; Tenice Nangoo; Sainath Lakshminarayanan
Priority date: 2011-12-20
Filing date: 2012-12-20
Publication date: 2013-06-27
Also published as: GB201219828D0; GB201121932D0; GB2504158A; WO2013093468A3

Abstract

There is provided a method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement. The method comprises obtaining a measured phase data set derived from a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points as a function of source and receiver location, and the measured phase data set comprising the phase component of the measured values taken at one or more selected frequencies and for said at least one physical parameter. Next, a modelled phase data set is obtained generated from a subsurface starting model, the modelled phase data set comprising the phase component of modelled values taken at the same plurality of discrete data points, at said one or more selected frequencies and for said at least one physical parameter. For at least one source or receiver position of the measured phase data set, a measured phase distribution data set is generated comprising the phase distribution for a single source or receiver location as a function of spatial location. Then, for said at least one source or receiver position of the modelled phase data set, a modelled phase distribution data set is generated comprising the phase distribution for a single source or receiver location as a function of spatial location. Once these sets have been generated, the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set can be analysed and a determination, based on said correlation, of the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set can be made.

Description

Full Waveform Inversion Quality Control Method

The present invention relates to a method of, and apparatus for, quality control of the full waveform inversion method for calculating a model of a portion of the Earth. More particularly, the present invention relates to a method of, and apparatus for checking the validity of the full waveform inversion procedure applied to measured seismic data when generating a subsurface model of at least a portion of the Earth.

There is considerable interest in surveying sections of the Earth to detect natural mineral resources or other sites of geological interest. One approach to detection of such natural features is seismic surveying. In a seismic survey, vibrational energy generated by a seismic source is directed into the surface of the Earth. Returning signals are then detected and analysed. The returning signals comprise reduced amplitude reflections or refractions of the original vibrational wave which have been reflected or refracted from boundaries between different rock layers or strata under the surface. These signals can be used to map the subsurface of the Earth.

A general schematic illustration of an experimental set up 10 for an undersea seismic survey is shown in Figure 1. However, an equivalent experiment can be carried out on land. The present invention is applicable to subsurface exploration in any suitable environment, for example land or marine measurements of a portion of the subsurface of the Earth. However, the present invention is applicable in general to the subsurface exploration of any planetary body or moon. The skilled person would be readily aware of the suitable environments in which data could be gathered for analysis and exploration purposes as set out in the present disclosure.

In this example, the experimental set up 10 comprises a source 12 (in this case located on a ship 14). The source 12 generates an acoustic wave having sufficient vibrational energy to penetrate the subsurface of the Earth and generate sufficient return signals to aid useful detection. The source 12 may comprise, for example, an explosive device (for example, comprising dynamite), an air gun or other mechanical device capable of creating sufficient vibrational disturbance. In most seismic survey experiments, a single source is used which is shot from multiple locations.

A plurality of detectors 16 is provided. The detectors 16 may comprise any suitable vibrational detection apparatus. In practice, devices known as geophones are used which detect particle motion and hydrophones which detect pressure variations. Commonly, a large number of detectors 16 are laid out in long lines (in the case of 2D data acquisition) or in sets of lines or in a grid (for 3D acquisition). The detectors 16 are connected to trace acquisition apparatus such as a computer or other electronic storage device. In this example, the acquisition apparatus is located on a further ship 18.

In use, the acoustic wave 20 generated by the source 12 travels downwards into the subsurface 22 of the Earth. The subsurface 22, in general, comprises one or more layers or strata 24, 26, 28 formed from rock or other materials. The acoustic wave 20 interacts with the medium, for example by being refracted through the layers and/or reflected off the interfaces between them, and a plurality of return signals 30 is detected by the detectors 16.

In general, the returning signals 30 comprise pressure or elastic waves having different polarisations. P-waves (also known as primary or pressure waves) are longitudinally polarised and comprises alternating rarefactions and compressions in the medium in which the wave is travelling. In other words, in an isotropic environment, the oscillations of a P-wave are parallel to the direction of propagation of the wave. However, the polarisation of P-waves can vary depending upon the angle at which they strike boundaries. P-waves have the highest velocity and so are the first to be recorded. P-waves travel at a velocity V_p in a particular medium. V_p is a scalar quantity and is, effectively, the speed of sound in a medium. It is this quantity V_p which is of particular interest in seismic inversion

S-waves (also known as shear or secondary waves) are also generated. S-waves have a transverse polarisation (i.e. perpendicular to the direction of travel). S-waves travel slower than P- waves in materials such as rock. Whilst S-wave analysis is possible and falls within the scope of the present invention, the following description will focus on the analysis of P-waves.

A seismic survey is composed of a large number of individual source excitation events. The earth's response to these events is recorded at each receiver location, as a seismic trace for each source-receiver pair. For a two dimensional survey, the tens of thousands of individual traces may be taken. For the three dimensional case, this number may run into the millions.

A seismic trace comprises a measurement, by the multiplicity of detectors 16, of the returning reflected or refracted acoustic waves 30 originating from the source 12. In general, a partial reflection of the acoustic wave 20 occurs at a boundary or interface between two dissimilar materials, or when the elastic properties of a material changes. Traces are usually obtained at discrete intervals of the order of milliseconds. Each detected return signal 30 forming a seismic trace has a travel time from the source 12 which, for reflected waves, is a two-way travel time from the source 12 to the reflecting element (usually subsurface rock interface) and back up to the detectors 16.

Seismic surveys at the surface or seabed can be used to extract rock properties and construct reflectivity images of the subsurface. Such surveys can, with the correct interpretation, provide an accurate picture of the subsurface structure of the portion of the Earth being surveyed. This may include subsurface features associated with cavities or pockets of mineral resources such as hydrocarbons (for example, oil and natural gas). Features of interest in prospecting include: faults; folds; anticlines; unconformities; salt domes, reefs.

Key to this process of modelling and imaging a portion of earth is the seismic velocity V_p. In a portion of the volume of the Earth, V_p may be estimated in various ways. One method is travel- time tomography. Other methods are velocity analysis from NMO or migration. In tomography, travel times of events are picked from the recorded data and the velocity that best predicts these travel times is found. Predicted travel times are calculated using ray-based approximations. This approach is robust in the sense that predicted travel-times do not vary significantly when the model is varied and produces a smooth velocity model on the scale of the Fresnel zone.

However, an alternative technique is full waveform inversion (FWI). FWI is a technique that seeks to extract the properties of subsurface rocks from a given seismic dataset recorded at the surface or seabed. In FWI, the aim is to go further and produce a more detailed velocity estimate with variations on the scale of a seismic wavelength. In this case, an attempt is made to match the full recorded wavefield with predicted data.

In other words, the technique involves generating a two or three dimensional model to represent the measured portion of the Earth and attempts to modify the properties and parameters of the Earth model to generate modelled data that matches the experimentally obtained seismic trace data.

The predicted data is calculated from the subsurface model using the full two-way wave equation. The final model has potentially far higher resolution and accuracy however the method can fail due to the sensitivity of the predicted waveforms to the model. FWI is an iterative process requiring a starting model. A sufficiently accurate starting model for FWI may be provided by travel- time tomography.

FWI can extract many physical properties (V_p and V_s velocities, attenuation, density, anisotropy) of the modelled portion of the Earth. However, V_p, the P-wave velocity, is a particularly important parameter which the subsequent construction of the other parameters depends heavily upon. The inversion of V_p is focussed on here. However, other parameters may be used with the present invention, either alone or in combination. The nature and number of parameters used in a model of a portion of the Earth will be readily apparent to the skilled person. The FWI technique seeks to obtain an accurate and high resolution V_p model of the subsurface which generates modelled data that matches the recorded data. Modelled data is calculated using the full two-way wave equation. This is known as forward problem. This equation can be in the time domain or the frequency domain, elastic or acoustic, isotropic or anisotropic. In most cases FWI proceeds using the acoustic approximation with a single component modelled wavefield which in the marine case is pressure.

An example of a basic starting model is shown in Figure 2. The model shows a simple estimation of the subsurface of a portion of the Earth. The source of acoustic waves is shown as a star and a plurality of receivers shown as circles. Both the source and the receivers are located at or above the seabed. As shown, the basic model shows a gradually increasing V_p with increasing depth without any of the detailed structure present in a Marmousi true earth model (which is shown in Figure 6a) and described later).

A modelled seismic gather is shown in Figure 3 for one shot and one line of receivers. The modelled seismic traces in the gather have been generated using the basic model shown in Figure 2. This is done by applying the isotropic acoustic wave equation to the model of Figure 2 and then modelling the reflected and refracted signals as they would be detected. The modelled seismic shot gather is made up of individual traces at surface receiver position showing pressure recorded as a function of time.

In general, the parameters of the model are estimated at a plurality of points set out in a grid or volume. The model is used to generate a modelled representation of the seismic data set. The modelled seismic data set is then compared to the real-world experimentally obtained seismic data set. Then, through use of a convergent numerical iteration, the parameters of the model are modified until the modelled seismic data set generated from the Earth model matches the actual measured seismic data to a sufficient degree of accuracy or until sufficient convergence is obtained. Examples of this technique are illustrated in "An overview of full-waveform inversion in exploration geophysics ", J. Virieux and S. Operto, Geophysics Vol. 74 No. 6 and US-A-7,725,266.

A general method to update a model will now be described. FWI operates on the principle of iteratively updating the starting model to minimise an objective function through repeated steepest- descent direction calculation. An objective function represents some measure of the mismatch between the recorded data and the modelled data.

Due to the non-linearity in the relationship between the model and the data, the objective function used in FWI will oscillate with changes in the model. This makes it necessary to have a sufficiently accurate starting model for global minimum convergence. It can be formulated in either the frequency domain or time domain. The choice of domain allows the use of pre-conditioning on either the data or the model update direction that could improve convergence or the linearity of the inverse problem.

Frequency domain inversion is equivalent to time domain inversion if all the frequencies are inverted simultaneously. However, the global minimum broadens at lower frequencies reducing how accurate the starting model needs to be for localised inversion to be successful. This makes it advantageous to extract the lowest useable frequencies in the recorded data and invert these first before moving to successive higher frequencies. This is illustrated in "Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies " Sirgue, L., and R. G. Pratt, Geophysics, 69, 231-248 (2004).

This multi-scale approach aims to first recover the long wavelength structures in Vp which is necessary for global minimum convergence and accurate subsequent reflectivity imaging. However local minimum convergence can occur in using existing techniques in practical situations where low frequencies are corrupted by noise and the starting model is limited in accuracy.

Localised gradient-based methods are used which iteratively update an existing model in a direction that derives from the objective function's direction of steepest descent. In its standard form it is the L-2 norm of the data residuals that is minimised. To describe this procedure it is sufficient to consider a single shot and single frequency in the objective function. The shot is located at point s and ω is the angular value of the frequency being inverted. Ultimately gradients from separate shots and frequencies will summed over when forming the final model update.

where u and d are the complex frequency components of the modelled and measured data respectively at receiver-source pair position (r,s) for the given frequency being inverted, u may be computed, usually by finite-difference methods, by solving the frequency domain (uniform-density) acoustic wave equation

2) ( ² +co²m²(x))u(x,s) = -SS(x-s) here the subsurface model m(x) is slowness, the inverse of V_p (which we seek to recover) and x denotes subsurface position. S is the complex point-source which in general is also unknown and needs to estimate from the data. The solution is linear in S:

3) u(x, s)=SG(x, s). where G is the green's function which we calculate by solving equation (2) with S=l.

The gradient direction of objective function (1) is evaluated at the existing model. This is the partial derivative of E for a point model perturbation at each subsurface position x

4) a¾ = R ∑^(»(^)-rf(r,*))^*

r

Then, from Equation 2):

(V² + ²m x))^e^=-2(o²m{x')u{x',s)5{x-x')

This is the same equation as (2) but with a source term now located at model perturbation position x'. The solution of equation 5) is a wave field emitted by a source positioned at*', given by equation 6)

6) ¾^=26)²ffl(x ( ',s)G(x, ') Next, the solution is evaluated at the receiver location, given by equation 7):

7) =

2(D²G(r, x')m(x')u(x', s)

Then, applying source-receiver reciprocity the wavefield is now emitted from the receiver location, as specified in equation 8):

8) dm(x') v \ > \ > Inserting equation 8) into equation 4) gives equation 9) for single receiver:

where:

10) w(x', r, s) = G (x' ,r) (u(r, s) - d (r, s))

The gradient of equation 9) is the product of two wavefields:

i) incident wavefield u emitted by a source at original source location s; and

ii) a back-propagated wavefield w which is emitted by a residual source at receiver location r. This expression is the gradient for a single residual. The gradient for a whole shot gather comes from summing equation 10) over all the receivers to give equation 1 1):

where:

w(x',s) = ) (u(r,s) - d(r, s))

12) y G^* (x' , r

The above can still be calculated as a single wavefield by taking a multi-point source. The sources are the wavefield residuals at each individual receiver location.

Finally, several frequencies and shots can be inverted simultaneously by summing the gradient contributions. As set out above, the 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. Ideally, the above method will lead to a convergence to a model which is a correct representation of a portion of the subsurface of the Earth. However, there are some critical difficulties associated with obtaining correct convergence.

Firstly, FWI methodology for the objective function above relies upon the Bom

approximation of the inverse problem. This requires that the starting model matches the observed travel times to within half a cycle. If the starting model data does not match major events in the recorded data to within half a wave cycle, FWI mis-converges to a cycle-skipped local minimum.

This is illustrated using the basic example above. In this illustration, the 5Hz frequency component of both the (synthetic) measured (i.e. real world) and modelled seismic data is used. The measured seismic data is expressed as a function of receiver (r) and source (s) pair position. This is the function d(r,s). The modelled seismic data is also expressed as a function of receiver (r) and source (s) pair position. This is the function u(r,s).

Next, the phase residual between d(r,s) and u(r,s) is determined. The principal value (between -π and + π) of the phase residual is denoted as φ and expressed in equation 13):

which is the wrapped difference of the individual phases of u(r,s) and d(r,s). The wrapped difference, or principal value of the difference, is the difference between the phases when the phase difference is restricted to the range -π to + π. In other words, if the phase difference is π and increases, the value of the phase difference will become -π and move back towards π.

Measured time-domain data for a single shot is shown in Figure 4a) together with the corresponding shot's modelled data in Figure 4b). Additionally, Figure 4 c) shows φ (as set out in equation 13)) as a function of receiver position.

As can be seen in Figure 4 c), at the zero offset receiver location the phase residual is zero, i.e. the phase is matched. However, moving away from the zero offset receiver position towards the left hand side of Figure 4 c), it can be seen that the phase increases towards π before moving to - π (due to the wrapped phase) then increases towards π before moving to -π again (due to the wrapped phase). The apparent discontinuities, or ±π transitions, can lead to convergence problems in numerical iterations.

At short offsets φ (r) is not cycle-skipped. However due to the ±π transitions, at longer offsets φ (r) is cycle-skipped, by first one and then two cycles moving towards the left hand side of Figure 4c). This is the number of cycles that the modelled and recorded data are misaligned by at 5Hz. When the number of cycle-skips is zero, the modelled data lies within plus or minus half a wave cycle of the recorded data.

Now, suppose that the model is updated. The initial φ (data set 1) is shown together with an updated φ (data set 2) after 4 iterations of the model. In general, updating the model leads to φ shifting toward zero. However, this has the implication that the modelled waveform is shifted to the nearest cycle of the recorded waveform. At short offsets this is acceptable. However, at longer offsets where φ is cycle-skipped by one or more cycles, the iteration converges to an incorrect, cycle- skipped zero.

Therefore, taking the position of the arrow in Figure 5, for correct convergence φ is moving upwards to convergence, whereas in fact it should be moving downwards through -π to a non-cycle skipped zero.

Consequently, as illustrated above, for successful convergence the modelled data must lie within half a wave cycle of the recorded data at the lowest useable frequency. In other words, the starting model must be sufficiently accurate to match the actual data to within half a wave cycle, otherwise FWI may mis-converge to a cycle-skipped local minimum.

For this condition to hold low frequency components or a particularly accurate starting model are required. However, in practice, low frequency data is too noisy to be used in meaningful calculations and the starting model is of limited accuracy. Consequently, cycle skipping is a problem that leads to errors in the convergence procedure and to inaccurate final models.

An example of this is shown in Figures 6 a) to 6 f). Figure 6 a) shows the synthetic Marmousi model generated to represent the subsurface structure of a portion of the Earth. Figure 6 b) shows an accurate starting model which, in general shape, follows that of the actual Marmousi model.

Therefore, performing FWI starting with the starting model of Figure 6 b) leads to accurate convergence as shown in Figure 6 c). Figure 6 c) shows the final, converged model which is a good representation of the original Marmousi model.

However, in practice, such an accurate starting model is unlikely given the lack of information about subsurface structures of real-world portions of the Earth. Figures 6 d) to 6 f) illustrate the operation of FWI when provided with a poor (inaccurate) starting model as shown in Figure 6 d).

Figure 6 e) shows the recovered model after convergence using only a single 5Hz frequency component. Figure 6 f) shows the full recovered model after subsequent higher frequency iterations. As shown, there are clear inaccuracies in the recovered model when compared to the original model of Figure 6 a). Features are not resolved clearly, and are inaccurately located. This is a result of cycle skipped mis -convergence and, in practice, would give a sub-optimal final result which could lead to inaccurate identification of natural resources or other sites of interest. A further example of this is given in Figures 7a) to 7f). Figures 7a) and 7d) show different starting models - Figure 7a) shows a basic starting model and Figure 7d) shows a more accurate starting model.

Figures 7b) and 7e) show the respective models after a single iteration of an FWI process. The same measured data set is used in each case. Further, Figures 7c) and 7f) show the respective models after 15 iterations. As shown, there is a significant difference between the models after only 15 iterations. It can be seen that the less accurate starting model is converging to an inaccurate final model and the inaccuracies will only be amplified with each further iteration. The inaccuracy of final results is of concern to prospectors and data interpreters alike and may potentially result in costly and inefficient misidentification of natural resources.

However, known FWI techniques do not provide a robust indication that mis-convergence has occurred, let alone provide tools as how to manage, measure and correct deficiencies. Therefore, interpretation of the accuracy of the convergence is difficult using known methods. As a result, it is difficult to determine precisely the accuracy of a final model or to determine which final model (where multiple starting models are considered) represents the most accurate depiction of the subsurface portion of the Earth being modelled.

Therefore, to date, known FWI methods suffer from a technical problem of that cycle skipped mis-convergence can lead to significant errors in a final, generated model of a subsurface portion of the Earth and that is difficult to determine when this has occurred. Consequently, a need in the art exists to provide a method that is operable to provide an indication regarding the accuracy of a final model obtained through the FWI method. The present invention seeks to address these issues.

According to a first aspect of the present invention, there is provided a method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising: a) obtaining a measured phase data set derived from a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points as a function of source and receiver location, and the measured phase data set comprising the phase component of the measured values taken at one or more selected frequencies and for said at least one physical parameter; b) obtaining a modelled phase data set generated from a subsurface starting model, the modelled phase data set comprising the phase component of modelled values taken at the same plurality of discrete data points, at said one or more selected frequencies and for said at least one physical parameter; c) generating, for at least one source or receiver position of the measured phase data set, a measured phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location; d) generating, for said at least one source or receiver position of the modelled phase data set, a modelled phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location; e) analysing the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set; f) determining, based on said correlation, the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set.

In one embodiment, step f) further comprises: g) comparing the spatial location of phase distribution features within the measured and modelled phase data sets with respect to the source or receiver position to determine the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set.

In one embodiment, said phase distribution features comprise a plurality of concentric ringlike features, and step g) further comprises comparing the radial spacing from the source or receiver location of a phase distribution feature within a selected ring-like feature in said modelled phase data set with respect to the radial spacing from the source or receiver location of a corresponding phase distribution feature in a corresponding ring-like feature in said measured phase data set.

In one embodiment, if, in a radial direction, a phase distribution feature within an n^th ring-like feature in said modelled phase data set is closest to a corresponding phase distribution feature in an m^th (where m≠n) ring-like feature in said measured phase data set, the method further comprises: h) rejecting the subsurface starting model.

In one embodiment, if, in a radial direction, a phase distribution feature within an n^th ring-like feature in said modelled phase data set is closest to a corresponding phase distribution feature in an n^th ring-like feature in said measured phase data set, the method further comprises: i) proceeding further with said subsurface starting model. In one embodiment, if, in a radial direction, a phase distribution feature within an n^th ring-like feature is directly adjacent a corresponding phase distribution feature in an n^th ring-like feature in said measured phase data set, step i) further comprises: j) validating said subsurface starting model.

In one embodiment, step i) further comprises: k) obtaining a further modelled phase data set generated from an «^Λ iteration (where n is an integer) of a full waveform inversion process carried out on said subsurface starting model, the further modelled phase data set comprising the phase component of further modelled values taken at the same plurality of discrete data points, at said one or more selected frequencies and for said at least one physical parameter; 1) generating, for said at least one source or receiver position of the further modelled phase data set, a further modelled phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location; m) analysing the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set and the further modelled phase data set; n) determining, based on said correlation, the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set. In one embodiment, step m) further comprises: o) analysing the change in phase distribution between the modelled phase distribution data set and the further modelled phase distribution data set to identify convergence or divergence from the measured phase distribution data set. In one embodiment, step o) further comprises: p) monitoring the location, distribution and/or movement of phase distribution rings within said phase data sets. In one embodiment, if, in step o) the phase distribution change from said modelled phase distribution data set to said further modelled phase distribution data set is divergent with respect to said measured phase distribution data set, the method further comprises: q) rejecting said starting model or modifying said starting model.

In one embodiment, if, in step o), the change in phase distribution from said modelled phase distribution data set to said further modelled phase distribution data set is convergent with respect to said measured phase distribution data set, the method further comprises :r) accepting said starting model or generating further iterations of said modelled data set.

In one embodiment, step a) further comprises: s) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points; t) applying a time-domain window to said measured time-domain values of said measured seismic data set to generate a windowed measured seismic data set; u) performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components; v) selecting said phase components for said measured phase data set; and wherein step b) comprises: w) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points; x) applying a time-domain window to said modelled time-domain values of said modelled seismic data set to generate a windowed modelled seismic data set; y) performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowed modelled values comprising amplitude and phase components; and z) selecting said phase components for said modelled phase data set.

In one embodiment, step t) and x) further comprises applying a shift to the measured or modelled data. In one embodiment, prior to step e), the method further comprises: aa) determining whether the spatial consistency of the measured and modelled phase distribution data sets meets a predetermined criterion.

In one embodiment, if, in step aa) said predetermined criterion is not met, the method further comprises: bb) selecting a new time-domain window parameter; and cc) repeating steps s) to z) with said new time-domain window parameter.

In one embodiment, the discrete data points of said measured and modelled seismic data sets are a function of source and receiver locations. In one embodiment, if it is confirmed in step m) that the optimised subsurface model is sufficiently accurate, the method further comprises the step: dd) utilising said optimised model for sub-surface exploration.

In one embodiment, said at least one physical parameter is pressure. In one embodiment, said at least one physical parameter comprises particle velocity or displacement. In one embodiment, the measured data set and the modelled data set comprise data taken at a plurality of discrete frequencies. In one embodiment, the measured data set and the modelled data set comprise values of a plurality of physical parameters. In one embodiment, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a full waveform inversion iterative process.

According to a second aspect of the present invention, there is provided a method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising: a) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points; b) applying a time-domain window to said measured time-domain values of said measured seismic data set to generate a windowed measured seismic data set; c) performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components;

d) generating, for at least one source or receiver position, a measured phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location; e) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points; f) applying a time-domain window to said modelled time- domain values of said modelled seismic data set to generate a windowed modelled seismic data set; g) performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowed modelled values comprising amplitude and phase components; and h) generating, for said at least one source or receiver position, a modelled phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location; i) determining whether the spatial consistency of the measured and modelled phase data sets meets a predetermined criterion; and, if predetermined criterion is not met, the method further comprises: j)selecting a new time-domain window parameter; and k) repeating steps b) to d) and f) to i) with said new time-domain window parameter until said predetermined criterion is met; and 1) determining the validity of said subsurface starting model for interpreting said measured seismic data set.

In one embodiment, step 1) comprises: m) receiving at least one further modelled seismic data set generated from an «^Λ iteration (where n is an integer) of a full waveform inversion process carried out on said modelled seismic data set, said at least one further modelled seismic data set comprising modelled values of at least one physical parameter at the same plurality of discrete data points; n) determining, for a plurality of discrete data points, a starting check residual variable from the measured values and the modelled values, the check residual variable of a discrete data point comprising the mismatch between a variable of the measured values and said variable of the modelled values of the at least one physical parameter for that discrete data point; o) calculating a change residual value for each of a plurality of discrete data points, the change residual for a discrete data point comprising the difference between the starting check residual variable and the check residual variable for the n^th iteration; and p) comparing the change residual of each of a plurality of discrete data points with the check residual variable of the n^th iteration for the same discrete data points to determine the validity of the starting model for interpreting said measured seismic data set.

In one embodiment, step 1) comprises: q) calculating, from the phase component of said windowed measured values and the phase component of said windowed modelled values, the phase residual at said one or more selected frequencies, the phase residual of a discrete data point comprising the mismatch between the phase of the windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point;

r) generating a phase residual data set at said one or more selected frequencies and for said at least one physical parameter; s) identifying cycle-skipped data within said phase residual data set; and t) evaluating the proportion of cycle-skipped data within said phase residual data set to determine the validity of the starting model for interpreting said measured seismic data set. In one embodiment, step 1) comprises: u) analysing the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set; and f) determining, based on said correlation, the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set.

According to a third aspect of the present invention, there is provided a method for subsurface exploration, the method comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement, and comprising the steps of:

a) providing, from said seismic measurement, a measured seismic data set obtained from said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain or frequency-domain values of at least one physical parameter measured at a plurality of discrete data points; b) generating, using a subsurface model of a portion of the Earth, at least one modelled seismic data set comprising modelled time-domain values or frequency-domain values of at least one physical parameter at the same plurality of discrete data points;

c) if, in step b) the data is provided in the time-domain, performing a Fourier transform on said modelled seismic data set to obtain modelled values taken at one or more selected frequencies, said modelled values comprising amplitude and phase components; d) if, in step a) the data is provided in the time-domain, performing a Fourier transform on said seismic measured data set to obtain measured values taken at one or more selected frequencies, said measured values comprising amplitude and phase components; e) calculating, from the phase component of said measured values and the phase component of said modelled values, the unwrapped phase residual at said one or more selected frequencies, the unwrapped phase residual of a discrete data point comprising the mismatch between the phase of the measured values and the modelled values of the at least one physical parameter for that discrete data point; f) generating a phase residual data set from the measured and modelled dataset; g) applying a weight to the discrete data points within said phase residual dataset in dependence upon the magnitude of the unwrapped phase residual for one or more discrete data points; h) utilising an FWI process to minimise the weighted phase residual dataset to produce an updated subsurface model; and i) utilising said updated model for subsurface exploration.

In one embodiment, step g) comprises weighting a particular discrete data point with a value between 0 and 1 in dependence upon the magnitude of the unwrapped phase residual for that discrete data point. In one embodiment, in step g) a discrete data point is given a value of 0 if said data point has a value of 180° (π) or more. In one embodiment, in step g) a discrete data point is given a value of 0 if said data point has a value of 160° (0.88π) or more.

In one embodiment, if a data point has a value of 0, said data point is ignored or deleted from the measured and the updated modelled seismic data sets. In one embodiment, when, in step g), a discrete data point has a non-zero value between 0 and 1, the data point is given a weighting to the minimisation in step h) based on said non-zero value, with a value of 1 being given full weight.

In one embodiment, step g) comprises weighting a particular discrete data point with a value of 0 or 1 in dependence upon the magnitude of the unwrapped phase residual for that discrete data point, data with a value of 0 being deleted from the measured and the updated modelled seismic data sets or ignored in step h).

According to a fourth aspect of the present invention, there is provided a method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising: a) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured values of at least one physical parameter measured at a plurality of discrete data points; b) receiving a starting modelled seismic data set generated from a subsurface starting model, said starting modelled seismic data set comprising starting modelled values of at least one physical parameter at the same plurality of discrete data points; c) receiving at least one further modelled seismic data set generated from an iteration (where n is an integer) of a full waveform inversion process carried out on said starting modelled seismic data set, said at least one further modelled seismic data set comprising modelled values of at least one physical parameter at the same plurality of discrete data points; wherein the method further comprises determining, for at least one iteration, the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set by: d) determining, for a plurality of discrete data points, a starting check residual variable from the measured values and the starting modelled values, the check residual variable of a discrete data point comprising the mismatch between a variable of the measured values and said variable of the modelled values of the at least one physical parameter for that discrete data point: e) calculating a change residual value for each of a plurality of discrete data points, the change residual for a discrete data point comprising the difference between the starting check residual variable and the check residual variable for the n^th iteration; and f) comparing the change residual of each of a plurality of discrete data points with the check residual variable of the n^th iteration for the same discrete data points to determine the validity of the starting model for interpreting said measured seismic data set.

In one embodiment, step d) further comprises: determining, for a plurality of discrete data points, a check residual variable from the measured values and the modelled values from said n^th iteration. In one embodiment, step e) further comprises: determining, for a plurality of discrete data points, a change residual variable from the starting modelled values and the modelled values from said n^th iteration; In one embodiment, the residual variable is a phase residual and step d) comprises calculating a phase residual for a plurality of discrete data points, the phase residual of a discrete data point comprising the mismatch between the phase of the measured values and the modelled values of the at least one physical parameter for that discrete data point. In one embodiment, the residual variable is an unwrapped phase residual and step d) comprises calculating an unwrapped phase residual for a plurality of discrete data points, the unwrapped phase residual of a discrete data point comprising the mismatch between the unwrapped phase of the measured values and the unwrapped phase of the modelled values of the at least one physical parameter for that discrete data point.

In one embodiment, said step of calculating an unwrapped phase residual comprises: g) calculating, for a plurality of groups of two discrete data points, the principal value of the difference in phase residual at said one or more frequencies between the discrete data points in a group such that each group has a principal value of the phase residual difference associated therewith; and h) determining, from phase residual difference of at least some of said plurality of groups, the unwrapped phase residual for at least some of said plurality of data points.

In one embodiment, step g) further comprises selecting a subset of groups of discrete data points from the total number of groups based on a selection criterion. In one embodiment, said selection criterion comprises selecting the groups having the smallest difference between the phase residuals of their constituent points.

In one embodiment, said selection criterion comprises selecting groups on the basis that the spatial separation between discrete data points within a group does not exceed a spatial threshold value. In one embodiment, said spatial threshold value is twice the minimum spatial separation of the discrete data points. In one embodiment, step g) further comprises applying spatial averaging to reduce the effects of noise. In one embodiment, step d) comprises updating said subsurface model of a portion of the Earth by a localised gradient-based method. In one embodiment, step d) further comprises: i) constructing an objective function; and j) minimising said objective function by modifying at least one coefficient of said subsurface model of a portion of the Earth to provide an updated subsurface model of a portion of the Earth.

In one embodiment, said objective function comprises a logarithmic objective function consisting of the unwrapped phase residual summed over a plurality of data points. In one

embodiment, step a) comprises performing a Fourier transform on said seismic measurement data to obtain a plurality of frequency components of said data set, and selecting one or more frequencies therefrom.

In one embodiment, wherein in step b) and/or c) said modelled seismic data set at a selected frequency is generated directly from the subsurface model of a portion of the Earth. In one embodiment, steps b) and/or c) comprises performing a Fourier transform on modelled seismic data generated from said subsurface model to obtain a plurality of frequency components of said data set, and selecting one or more frequencies therefrom.

In one embodiment, the check residual variable is a time-offset residual and step d) comprises calculating, for a plurality of discrete data points, the cross-correlation in the time-domain between modelled and measured data at various time offsets between the datasets, with the time-offset residual for a discrete data point comprises being the time offset that maximises the cross-correlation.

In one embodiment, step f) comprises comparing the sign of the change residual with the sign of the residual variable of the n^th iteration for at least some of said plurality of data points to determine the validity of the starting model for interpreting said measured seismic data set.

In one embodiment, step f) further comprises determining a check ratio for at least some of said plurality of data points, where the check ratio is the ratio of the number of data points for which the sign of the change residual is the same as the sign of the residual variable of the n^th iteration to the number of data points for which the sign of the change residual is different from the sign of the residual variable of the n^th iteration. In one embodiment, the method further comprises the step of deleting the data points from the measured data set for which sign of the change residual is different from the sign of the residual variable of the n^th iteration to generate a modified measured data set. In one embodiment, the method further comprises utilising said modified measured data set in an FWI minimisation process to generate a new starting model. In one embodiment, the method further comprises repeating steps a) to f) utilising the new starting model.

In one embodiment, the validity of the starting model for interpreting said measured seismic data set is determined at least in part based upon said check ratio. In one embodiment, the starting model is considered valid if the check ratio is 50% or greater. In one embodiment, the starting model is considered valid if the check ratio is 70% or greater. In one embodiment, step f) comprises calculating the ratio of the change residual to the residual variable of the n^th iteration for at least some of said plurality of data points. In one embodiment, the discrete data points of said measured and modelled seismic data sets are a function of source and receiver locations.

In one embodiment, step f) further comprises comparing the change residual with the residual variable of the n^th iteration be performed as a function of receiver and/or source locations. In one embodiment, based upon the comparison in step f), the method may comprise one of: switching to a different full waveform inversion method; repeating the full waveform inversion again starting from the starting model; or modifying the starting model.

In one embodiment, the method further comprises the step of k) confirming the accuracy of an optimised subsurface model generated using said starting model. In one embodiment, the method further comprises the step of: 1) confirming or rejecting the validity of said sub surface starting model based upon step f). In one embodiment, if it is confirmed in step 1) that the optimised subsurface model is sufficiently accurate, the method further comprises the step:

m) utilising said optimised model for sub-surface exploration. In one embodiment, said at least one physical parameter is pressure. In one embodiment, said at least one physical parameter comprises particle velocity or displacement. In one embodiment, the measured data set and the modelled data set comprise data taken at a plurality of discrete frequencies. In one embodiment, the measured data set and the modelled data set comprise values of a plurality of physical parameters. In one embodiment, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a full waveform inversion iterative process.

According to a fifth aspect of the present invention, there is provided a method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising: a) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points; b) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points; c) applying a time-domain window to said modelled time-domain values of said modelled seismic data set to generate a windowed modelled seismic data set; d)performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowed modelled values comprising amplitude and phase components; e) applying a time-domain window to said measured time-domain values of said measured seismic data set to generate a windowed measured seismic data set; f) performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components; wherein the method further comprises determining the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set by: g) calculating, from the phase component of said windowed measured values and the phase component of said windowed modelled values, the phase residual at said one or more selected frequencies, the phase residual of a discrete data point comprising the mismatch between the phase of the windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point; h) generating a phase residual data set at said one or more selected frequencies and for said at least one physical parameter; i) identifying outlying data within said phase residual data set.

In one embodiment, said outlying data comprises cycle-skipped data and the method further comprises j) evaluating the proportion of cycle-skipped data within said phase residual data set to determine the validity of the starting model for interpreting said measured seismic data set. In one embodiment, step j) comprises calculating the proportion of cycle-skipped data as a function of the total data in the phase residual data set. In one embodiment, the starting model is rejected if the proportion of cycle-skipped data exceeds a pre-determined threshold. In one embodiment, the pre-determined threshold is 40%. In one embodiment, the pre-determined threshold is 20%.

In one embodiment, the method further comprises: k) repeating steps b) to d) for a further modelled seismic data set generated from an «^Λ iteration (where n is an integer) of a full waveform inversion process carried out on said modelled seismic data set to generate further windowed modelled values comprising amplitude and phase components; 1) repeating steps g) to h) utilising said further windowed modelled values to generate a further phase residual data set the phase residual at said one or more selected frequencies; m) identifying cycle-skipped data within said further phase residual data set; and wherein step j) further comprises: n) comparing the proportion of cycle skipped data within said further phase residual data set to said proportion of cycle-skipped data within said phase residual set to determine the validity of the starting model for interpreting said measured seismic data set.

In one embodiment, step n) further comprises determining the change in proportion of cycle- skipped data between said phase residual set and said further phase residual set. In one embodiment, if, in step n), the proportion of cycle-skipped data increases from said phase residual set to said further phase residual set, the method further comprises: o) rejecting said starting model.

In one embodiment, if, in step n), the proportion of cycle-skipped data decreases from said phase residual set to said further phase residual set, the method further comprises: p) accepting said starting model or generating further iterations of said modelled data set.

In one embodiment, the method further comprises: q)repeating steps b) to d) for a further modelled seismic data set generated from a different starting model to generate further windowed modelled values comprising amplitude and phase components; r) repeating steps g) to h) utilising said further windowed modelled values to generate a further phase residual data set the phase residual at said one or more selected frequencies; s)identifying cycle-skipped data within said further phase residual data set; and wherein step j) further comprises: t) comparing the proportion of cycle skipped data within said further phase residual data set to said proportion of cycle-skipped data within said phase residual set to determine which of said starting model and said further starting model is most suitable for interpreting said measured seismic data set. In one embodiment, step t) further comprises determining the change in proportion of cycle- skipped data between said phase residual set and said further phase residual set. In one embodiment, if, in step t), the proportion of cycle-skipped data increases from said phase residual set to said further phase residual set, the method further comprises: u) rejecting said further starting model. In one embodiment, if, in step t), the proportion of cycle-skipped data decreases from said phase residual set to said further phase residual set, the method further comprises: v) rejecting said starting model.

In one embodiment, the discrete data points of said measured and modelled seismic data sets are a function of source and receiver locations. In one embodiment, the method further comprises the step of: w) confirming the accuracy of an optimised subsurface model generated using said starting model. In one embodiment, the method further comprises the step of: x) confirming or rejecting the validity of said sub surface starting model based upon step j).

In one embodiment, if it is confirmed in step x) that the optimised subsurface model is sufficiently accurate, the method further comprises the step: y) utilising said optimised model for subsurface exploration. In one embodiment, said at least one physical parameter is pressure. In one embodiment, said at least one physical parameter comprises particle velocity or displacement. In one embodiment, the measured data set and the modelled data set comprise data taken at a plurality of discrete frequencies.

In one embodiment, wherein the measured data set and the modelled data set comprise values of a plurality of physical parameters. In one embodiment, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a full waveform inversion iterative process.

According to a sixth aspect of the present invention, there is provided a computer program product executable by a programmed or programmable processing apparatus, comprising one or more software portions for performing the steps of any one of the first to fifth aspects.

In one embodiment, the method further comprises the step of, subsequent to step h), unwrapping the phase residual data set to generate an unwrapped phase residual data set. In one embodiment, the method further comprises the step of deleting the data points associated with outlying data in step i) from the measured seismic data set to generate a modified measured data set. In one embodiment, the method further comprises the step of deleting the data points from the measured data set for which sign of the change residual is different from the sign of the residual variable of the n^th iteration to generate a modified measured data set.

In one embodiment, the method further comprises utilising said modified measured data set in an FWI minimisation process to generate a new starting model. In one embodiment, the method further comprises repeating steps a) to f) utilising the new starting model.

According to a seventh aspect of the present invention, there is provided a computer usable storage medium having a computer program product according to the sixth aspect stored thereon.

Embodiments of the present invention will now be described in detail with reference to the accompanying drawings, in which:

Figure 1 is a schematic illustration of a typical seismic survey experiment in which seismic traces are obtained from an undersea portion of the Earth;

Figure 2 is a schematic illustration of a basic starting model for full waveform inversion modelling;

Figure 3 is a schematic illustration of modelled seismic trace data generated from the basic starting model of Figure 2 for an individual seismic shot;

Figure 4 a) is a synthetic example of measured seismic trace data taken from a seismic survey for an individual seismic shot;

Figure 4 b) is the seismic trace data of Figure 3 generated from the basic starting model of Figure 2 for an individual seismic shot;

Figure 4 c) is the phase residual between the measured data of Figure 4 a) and modelled data of Figure 4 b);

Figure 5 shows the phase residual of Figure 4 c) together with the phase residual after 4 iterations of a standard full waveform inversion minimisation iteration; Figure 6 a) to f) show examples of starting models and the resulting final models obtained from conventional full waveform inversion modelling;

Figure 7 a) - f) illustrates the effect of starting model on subsequent iterations;

Figure 8 is a flow chart illustrating the method according to a first embodiment of the present invention;

Figure 9a) illustrates a modelled seismic data set;

Figure 9 b) illustrates a measured seismic data set;

Figure 9 c) illustrates the phase residual between the data of Figures 7 a) and 7 b); Figure 10 a) illustrates a phase residual for a starting model plotted as a function of source and receiver location;

Figure 10 b) illustrates a change residual plot for the data of Figure 10 a) Figure 11 is a flow chart illustrating the method according to a second embodiment of the present invention;

Figure 12 a) illustrates the phase residual of Figure 9 c) on a larger scale; Figure 12 b) illustrates phase unwrapping of the data of Figure 12 a);

Figure 13 illustrates wrapped phase as a function receiver and source location showing terminating transitions in the prior art; Figure 14 is similar to figure 13 but shows windowed wrapped phase and data points and groups comprised thereof;

Figure 15 illustrates unwrapped data; Figure 16 a) illustrates a phase residual for a starting model plotted as a function of source and receiver location similar to Figure 10 a); and

Figure 16 b) illustrates a change residual plot for the data of Figure 16 a);

Figure 17 is a flow chart illustrating the method according to a third embodiment of the present invention;

Figures 18a) - d) are schematic illustrations of windowing modelled seismic trace data;

Figure 19a) - d) are schematic illustrations of windowing measured seismic trace data;

Figures 20a) to i) illustrate phase residual plots and resulting models for different starting models; and

Figure 21 is a flow chart illustrating the method according to a fourth embodiment of the present invention;

Figure 22 is a schematic of a source-receiver configuration for the following examples;

Figures 23 a) to f) illustrate two different starting models and resultant phase residual plots;

Figure 24 is a flow chart illustrating the method according to a fifth embodiment of the present invention;

Figure 25 is a flow chart illustrating the method according to a sixth embodiment of the present invention;

Figures 26 a) to g) illustrate a measured phase data plot, two different starting models and resultant modelled phase data plots;

Figures 27 a) to d) illustrate a measured phase data plot and resultant modelled phase data plots for a good starting model; Figures 28 a) to d) illustrate a measured phase data plot and resultant modelled phase data plots for a poor starting model;

Figure 29 is a flow chart illustrating the method according to a seventh embodiment of the present invention;

Figures 30 a) and 30 b) illustrate measured and measured windowed time domain data; and

Figures 31 a) and 31 b) illustrate phase plots for unwindowed and windowed data respectively.

As described above, FWI analysis can extract many physical properties (attenuation, density, anisotropy, elastic) from the parameters of an accurate Earth model. However, the essential parameter is V_p, the P-wave velocity. The present invention is directed towards a method for verifying the accuracy of a V_p model of the subsurface of a portion of the Earth generated by a FWI method.

As described above, FWI is an effective procedure for generating a quantitative detailed model of the subsurface for use in the exploration process. However it is a localised optimisation and therefore is strongly dependent on where the optimisation begins from.

The present method provides a tool to assess whether a given starting model is adequate for FWI on a particular seismic dataset i.e. whether it will converge successfully or not. This is done by highlighting modelled data shifting in the wrong direction. In other words, data shifting away from the recorded waveform at any given iteration.

A method according to a first embodiment of the present invention will now be described with reference to Figure 8. Figure 8 shows a flow diagram of an embodiment of the present invention. The first embodiment will be described with reference to a verification of the suitability of a subsurface starting model for use in interpreting a particular measured data set. Therefore, the general principle of the invention could be applied in a number of situations.

For example, the starting subsurface model, initial modelled data and a first iteration of the modelled data could be provided along with the measured data in order to provide guidance or verification of the likelihood of success of the model reaching an accurate convergence using that particular starting model following conventional FWI.

Alternatively, an already-run convergence procedure could be critically analysed by reviewing the residuals at each iteration of the FWI method to determine the validity of the final result. As a further alternative, the method could be run alongside or within an FWI method to provide an indication as to the accuracy of the FWI process as the iterations are computed.

Step 100: Obtain measured seismic data set

The starting point for the subsurface exploration analysis is a set of experimentally gathered seismic trace data. This may be gathered by an experimental arrangement such as the experimental set up shown and described with reference to Figure 1. As shown in Figure 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 in a line (for a two dimensional measurement) or in a grid (for a three dimensional measurement). The position of the detectors 16 is well known from location tracking devices such as GPS devices. Additionally, the location of the source 12 is also well known by location tracking means.

The measured seismic data set may comprise multiple "shots" or source 12 emissions. An example of seismic trace data for a single shot is shown in Figure 4 a). The gathered data is essentially a plot of pressure as a function of receiver position (on the x-axis) and time (on the y- axis). This is because, in general, a detector such as a hydrophone measures the scalar pressure at its location.

The seismic trace data comprises a plurality of measured data points. Each measured discrete data point has six associated location values - three spatial directions (x, y and z) for receiver (or detector) position (r), and three spatial dimensions (x, y, z) for source location (s), together with pressure magnitude data. The spatial coordinates for each discrete data point define its location in space. It also comprises one or more measurement parameters which denote the physical property being measured. In this embodiment, a single measurement parameter, pressure is utilised.

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

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 invention. The present invention simply requires a real-world measured seismic data set upon which analysis can be performed to facilitate subsurface exploration of a portion of the Earth. This may be supplied by a third party.

The measured seismic trace data comprises data in the time domain, i.e. measured pressure as a function of position and time. The data can be analysed and FWI performed on time-domain data, and this is intended to comprise an aspect of the present invention.

Alternatively, FWI can be performed in the frequency domain. The measured data will, however, contain multiple frequencies and, in order for data to be analysed against modelled data, it is necessary to extract data taken at discrete frequencies. Therefore, if frequency domain analysis is to be performed, it is necessary to resolve the different frequency components within the seismic trace data to form one or more measured data sets with the or each data set representing the seismic traces at a selected frequency. In other words, a measured seismic data set is constructed which comprises the frequency components of the trace data at discrete data points. The measured data set may contain one or more measurement parameters and one or more frequency components.

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.

One or more measured data sets at one or more frequencies may be extracted. For example, a plurality of measured 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. Alternatively, these frequencies can be inverted simultaneously. d(r,s) is the notation we give for the measured data set when it has one or more measured physical parameters and one or more frequency components.

By way of an example, Figure 9a) shows actual, measured time-domain data. In both cases, the data represents a single source shooting into a line of (461) receivers. Fourier transforming the time-domain data gives complex-valued wavefield d(r,s). The method now proceeds to step 102.

Step 102: Obtain starting model

At step 102, an initial starting model of the specified subsurface portion of the Earth is, optionally, provided, potentially by a third party. 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 invention is applicable to three dimensional approaches. Examples of typical starting models are shown in Figure 7.

The generated model consists of values of the physical parameter V_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.

It is noted that this step is optional and, instead, a third party may supply merely the starting modelled data set as set out in step 104. The method then proceeds to step 104.

Step 104: Obtain starting modelled data set

In step 104, a starting data set is obtained. The actual generation of the data in this step is, as described above, not material to the present invention. The data may be merely provided or obtained from another source. Instead, the present invention relates to confirmation of the validity of the generated data. However, for a complete understanding of the approach used, the method to generate the starting modelled data set will be described.

In order to proceed with the FWI method in order to model the subsurface region under investigation, it is necessary to utilise modelled data which corresponds to the same source-receiver location data position as the actual measured trace data so that the modelled and measured data can be compared.

In other words, the modelled data set corresponds discrete point to discrete point to the measured dataset. The modelled data set may be generated in the time domain in order to perform time domain analysis with respect to the measured data set obtained in step 100. Alternatively, the modelled data is generated for the same measurement parameter(s) at the same frequency or frequencies in order to perform frequency-domain analysis.

The modelled data set is generated using the full two-way wave equation. Two methods may be employed here. In one approach, the modelled trace data may be generated using the time domain full two-way wave equation and then performing a Fourier transform on that data in common with that done for the measured data in step 102.

Alternatively, the required frequency components may be calculated directly using the frequency domain version of the full two wave equation being considered.

The isotropic acoustic constant-density version of this is equation 14) with shot location s :

14) ( ¹ + <O²m² (x)} u(x, s) = -S (x - s) where u (x, s) is the modelled data set at location x , m(x) is the subsurface model parameter, slowness (the inverse of V_p), and ω is the angular value of the frequency being modelled.

From the above equation, modelled seismic data can be generated for one or more physical parameters and at one or more selected frequencies. This forms the starting modelled seismic data set ii₀(r,s).

An example is illustrated in Figure 8b) which shows modelled time-domain data. When Fourier transformed, this becomes the starting modelled seismic data set udr,s). The method now proceeds to step 106. Step 106: Obtain further modelled data set

The data obtained in step 104 comprises modelled seismic trace data in the time domain (i.e. modelled pressure at the location of the measured data points as a function of position and time) or in the frequency domain as discussed.

The further modelled data set required in step 106 is obtained after an iteration of FWI has been carried out on the modelled data of step 106. As set out in relation to step 104, the actual generation of the data in this step is, as described above, not material to the present invention. The data may be merely provided or obtained from another source. Instead, the present invention relates to confirmation of the validity of the generated data. However, for a complete understanding of the approach used, the method to generate the further modelled data set will be described in the context of a single iteration. In the case of an FWI method, a trial velocity is updated to minimise an objective function measuring the mismatch between modelled and recorded data. The objective function and the steepest descent direction used to update the model in conventional FWI is described in equations 1) to 12). This can be done in the time or frequency domain. A further modelled data set obtained from an iteration of the method described above is then provided. This forms the data set ui(r,s). The method then proceeds to step 108.

Step 108: Determine check residual As set out above, the measured seismic data is expressed as a function of receiver (r) and source (s) pair position. This is the function d(r,s). The modelled seismic data is also expressed as a function of receiver (r) and source (s) pair position. This, as set out above, is the function udr,s) for the starting modelled data set and ui(r,s) for the modelled data set generated from the first iteration of the procedure as obtained in step 106. The modelled and measured data sets u and d may comprise values of one or more physical parameters (e.g. pressure) taken at one or more frequencies.

In step 108, a residual variable is obtained for each of the modelled data sets, i.e. between the data sets d(r,s) and uo(r,s) and between d(r,s) and ui(r,s). The principal value (between -π and + π) of the phase residual is denoted as φ and expressed in equation 15) for an iteration n :

which is the wrapped difference of the individual phases of u_n(r,s) and d(r,s). The wrapped difference is the difference between the phases when the phase difference is restricted to the range -π to + π. In other words, if the phase difference is π and increases, the value of the phase difference will become -π and move back towards π. This can readily be visualised on an Argand diagram.

The phase residual φ_η(*",*) can be determined for the data points. In that case, each spatial data point now has associated therewith a value of the phase residual φη(*",«) for an iteration n. Note that the starting modelled data set, udr,s) is generated prior to the first iteration (for which n =1). In the case of the starting modelled data set uo(r,s), φο(/*,«) can be calculated from equation

15). The phase residual for one or more subsequent n^th iterations, φ_η(ι%8) can also be calculated in this manner for subsequent steps of the method. However, this is not necessary and the change phase residual can be indirectly obtained for further stages of the calculation by utilising the mathematical relationship between uo(r,s) and u_n(r,s) as will be described later. The method then proceeds to step 1 10.

Step 110: Unwrap phase?

The phase residual (or phase mismatch) for the data shown in Figures 9a) and 9b) is shown in Figure 9c). The phase residual is calculated in accordance with equation 15) as set out above. In Figure 9c), the phase residual data is shown at 3.125Hz and is plotted as a function of receiver number.

As shown in Figure 9c), the phase residual is zero at the source location (or zero-offset point) and oscillates within +π to -π until the receiver number increases beyond about 320. At this point, the phase residual increases beyond +π and wraps to -π. The phase residual increases beyond +π again at around receiver number 400 and wraps around to -π. This data set exhibits what is known as cycle skipping. The effect of cycle skipping can be to distort or render ineffective the determination of the validity of the starting model. In step 1 10, it is determined, based on the phase residual determined in step 108, whether the phase residual comprises cycle-skipped data and whether further processing to address this is required. If further processing is required, then the method proceeds to step 112. Otherwise, the method proceeds to step 114.

As an alternative, if it is desired to carry out the analysis with time domain data, the cross- correlation of the modelled and measured data at each receiver location for various time-offsets would be calculated. Then the time offset, T, would be selected corresponding to the value of maximum cross -correlation as our check residual variable.

Step 112: Unwrap phase residual

In dataset where we have such occurrences of 2 π wrap around we phase unwrap by adding multiples of 2 π to the phase residual prior to analysis (Figure 12). Whilst straightforward to do in one dimension we describe how to do this in a multi-dimensional way and attain a result consistent across all dimensions.

In real-world data, φ(τ,«) is a function of three spatial dimensions for both receiver and source. Therefore, φ(τ,«) will be a function of r_x, r_y, r_z, s_x, s_y, s_z and frequency / However, for clarity, consider a two dimensional case where Φ(Γ,«) is a function of just one spatial dimensions for each of the receiver and source. Therefore, φ(τ,«) is a function of r_x and s_x. Equation 16) illustrates the use of differential phase g to calculate the unwrapping: 16) g= phase w(r₂,s₂)i/*(r₂,s₂)w *(r₁, s_\)d(ri, s_\) which enables calculation of the phase residual difference between data points (r_\, s_\) and (r₂,

¾). Now the phase residual difference defined in equation 16) can be utilised for multidimensional unwrapping of the phase. The phase residual difference is taken between adjacent discrete data points. In other words, adjacent discrete data points are grouped into groups of two associated discrete data points. Then, the difference in phase residual between each of the discrete data points in a group is calculated. A phase residual difference value is then associated with that group.

The phase residual difference g values then form a data set which can be utilised to solve for unwrapped phase. The phase residual difference data set is for one or more physical parameters and one or more frequencies, as is the case for the measured and modelled data sets. The phase is thus unwrapped by solving equation 17):

17) L <Mr,_S)= g where:

18) φ (r,s)= unwrapped phase (u(r,s)d*(r,s))

And where g is the phase residual difference (or differential phase residual) and L is a matrix operator defining the groups (or pairs) of points selected in g. We select at least one pair for each of the data points where phi' is to be calculated.

The matrix operator L is arranged such that each column corresponds to a particular data point and each row corresponds to a particular pair. Each row takes values of 1 and -1 in the columns corresponding to the two points belonging to the pair, and 0 elsewhere.

Equation 17) is then solved as a system of equations. For an over-determined system a least- squares solution may be used. One effective such solution would involve applying a QR factorization to the matrix L.

Further the system can be weighted for example according to the ranks given to the selected groups. For instance those groups with a greater value of g or greater space distance between data points would be de-weighted and contribute less to the unwrapped phase solution.

As an additional step, it is possible to apply a spatial averaging to g by utilising a low-pass filter to reduce the effects of noise prior to solving the computation.

As set out above, selection of groups and points to enable unwrapping are important. Spatial separation is related to both the source separation and receiver separation for the data points in the group.

In the context of this example, the two discrete data points in each group are spatially nearest neighbours. This is shown schematically for a two dimensional example in Figure 14 which shows data points 170 superimposed over the phase residual plot of Figure 13. In a set of data points 170, a number of groups 180 are shown. Any configuration of groups 180 may be envisaged, and one data point may form part of multiple groups as shown in Figure 14. As shown, the data points 170 are not uniformly arranged. This is because, in practice, the receiver locations are never evenly spaced due to experimental constraints. The phase residual difference data set is formed from the collection of groups each having a phase residual difference value.

In one embodiment, it is desired that the spacing between data points 170 forming a group 180 does not exceed a threshold value. In the present embodiment, it is desirable that the spacing between two data points 170 in any group 180 does not exceed twice the minimum spacing. To explain, since the receiver spacing is defined by the measured data set (where each point corresponds to a point in space where a hydrophone or geophone is located), the spacing between points can also be represented as a physical distance. Therefore, for example, if the minimum spacing between receivers in the measured data set is 20 m, then the threshold for spacing between two data points 170 in a group 180 should, in this embodiment, not exceed 40 m.

Additionally, it is envisaged that the term adjacent data points may include groups such as one group 180 which is arranged on a diagonal - i.e. the two data points are adjacent one another but not necessarily nearest neighbours. In addition to the example given above, a group may be formed of data points which are further spaced apart. For example, a group may be formed from a discrete data point and a next-but-one data point. As a further alternative, any two data points may form a group.

Once all the value of the phase residual difference has been calculated for some or all of the groups of adjacent data points, and a phase residual difference data set generated, a subset of those groups can be selected for further processing. Therefore, n_g groups (where n_g is a subset of the total number of groups) are selected. As set out above, each discrete data point has three coordinates for source spatial location and three coordinates for receiver spatial location. This defines the spatial location of the discrete data points and, concomitantly, their separation within a group.

The selection criteria for the n_g groups are primarily based on the spatial separation of the discrete data points. In one embodiment, the pairs are ranked in order of the pairs having the smallest difference between the phase residuals of their constituent points. Alternatively, the criterion of twice the minimum separation could be implemented, with the n_g groups having a point separation of twice the minimum group spacing or less, with the remaining groups being unselected.

As an alternative, adjacent or nearest neighbour discrete data points could be selected. As a further alternative, the groups comprising the closest spaced discrete data points could be selected for the subset n_g. Additionally, should more than the desired n_g groups fall within the spatial separation threshold, then the groups chosen to form part of the subset are selected on the basis of their value of phase residual difference.

The phase residual difference only provides meaningful analysis and convergence if the phase residual values are less than a threshold value. In this embodiment, the threshold value is π/2 (or 90 degrees). However, other threshold values may be used. In other words, the groups which are chosen to form part of the subset n_g of groups have the smallest values i.e. the values which are furthest away from the maximum possible values of +/- π.

As set out above, the groups may be selected on the magnitude of the wrapped phase residual difference as long as their spatial separation falls below a threshold value. Therefore, the n_g groups having the smallest values of the wrapped phase residual difference could be selected for further processing.

In this example as described, only one component data (pressure) is used. Therefore, there is only a single value of Φ(Γ,«) or of the phase residual difference between points. However, in many situations, multiple variables may be used. For example, if pressure and particle displacements are measured, then there would be values of Φ(Γ,«) for each direction of particle motion and one for pressure - leading to four separate Φ(Γ,«) values.

Once the phase residual has been unwrapped, the exemplary data plot of Figure 13 then becomes the exemplary plot shown in Figure 15. As shown in Figure 15, the phase unwrapping removes the cycle skip boundaries in the data plot.

In this example, φ (r,s) is defined as the check residual variable. Additionally or alternatively, the differential phase residual, g (equation 16) could also be used directly as our check residual. Note that steps 110 and 112 described above are optional and the phase residual may simply be calculated without unwrapping. Further other residuals may be utilised with the present invention. The skilled person would be readily aware of suitable residual variables for use with the present invention. The method now proceeds to step 114.

Step 114: Calculate change residual from the check residuals

As described above, the present method provides a tool to assess whether a given starting model is adequate for FWI on a particular seismic dataset i.e. whether it will converge successfully or not under conventional FWI. This is done by highlighting modelled data shifting in the wrong direction i.e. away from the recorded waveform, at any given iteration.

In order to determine how the data is changing with each iteration, the change in the residual variable between iterations (or, as described here, between the starting data set and the first iteration) is determined and analysed.

Thus, in step 114 the change residual Δφ is defined. The residuals discussed here focus on frequency domain analysis, although time domain analysis is to be considered within the scope of the present invention. The change residual Δφ is thus defined in equation 19) for the case where the starting modelled data set uo(r,s) and the first iteration modelled data set u i(r,s) are used:

19) Δφ = phase ((u₀d*) (uid*)*)= phase(w₀Wi *) or, in a more general form (with dependency on source and receiver locations shown) for an iteration «. ·

20) Δφ_η (Γ,«)= phase (u_n (r,s) u_n+1(r,s)*)

Equations 19) and 20) describe the change in phase residual caused by a model update. Mis- convergence can then be detected by comparing Δφ_η with φ_α We note that Δφ_η (r,s) is retained as a wrapped quantity lying between - π and + π. Further it is not explicitly dependent on d so is noise- free. Further, as described above, Δφ_η (r,s) can be calculated directly from u₀(r,s) and u i (r,s) and so it is not explicitly necessary to calculate or obtain the check residual for iteration 1 (φ1(Γ,8)) directly, although this may be done if required. φ_η (r,s) and Δφ_η (r,s) may be compared at one or many frequencies. These frequencies need not be those used in the inversion procedure. In addition, the change residual could be calculated for g as set out in equation 21 ):

21 ) Ag_n (r,s)= phase w„(r₂,s₂)w„*(r₂,s₂)w,_!+; *(r₁, si)u_n+1(ri, si)

If using the time-offset variable T, we would here calculate the change variable delta T. Step 116: Compare change residual with check residual variable The mis -convergence or otherwise of the data after an iteration can be gleaned through a comparison of Δφ_η with (f)_nas calculated in the previous step of the method for frequency domain analysis. Figure 9a) shows a graphical example of φ_η plotted as a function of receiver and source position. In Figure 9a), φ_η is taken at a frequency of 3.125Hz plotted as a function of s_x and r_x. Figure 9b) shows a corresponding plot of Δφ_η plotted as a function of the same receiver and source positions after 1 iteration of conventional FWI.

As shown by a comparison between these figures, in the circled area it can be seen that φηίβ greater than zero whereas Δφ_η is less than zero. Therefore, the change residual is moving in the opposite direction to the phase residual and so the data is shifting in the wrong direction. Therefore, from a general, qualitative, analysis of the data, it can be stated that the starting model for the iterations shown is inadequate and further iterations will only cause the model to converge to an inaccurate local minimum. Consequently, the method concludes that the given ID starting model is inadequate for conventional FWI. Therefore, as a general check method, where φ_η and Δφ_η oppose each other (i.e. have a different sign from one another), the modelled data shifting away from the recorded data. We could similarly compare g_n and Ag_n A further similar approach can be applied when processing data in the time domain. In this case, compare residual T_n with ΔΤ_η A more complex analysis of the change residual can provide further detailed information to vindicate or reject a particular starting model. This is an advantage of the present method which enables the divergence or mis-convergence of the data to be observed as a function of source and receiver position. This is in contrast to any known methods which may simply sum the errors over the whole dataset. In this example we see the long offsets are shifting in the wrong direction. This validates what we expect in a cycle-skipped example such as shown in Figure 9c).

A more complex analysis may be to determine a check ratio for at least some of said plurality of data points. The check ratio is the ratio of the number of data points for which the sign of the change residual Δφ_η is the same as the sign of the residual variable of the n^th iteration φ_η to the number of data points for which the sign of the change residual Δφ_η is different from the sign of the residual variable of the n^th iteration φ_η.

This ratio may be expressed as a percentage and the validity of the starting model for interpreting said measured seismic data set is determined at least in part based upon the check ratio. As examples, the starting model may be considered valid if the check ratio is 50% or greater.

Alternatively, the starting model may be considered valid if the check ratio is 70% or greater.

Step 118: Complete analysis

At step 118, the analysis can then be completed and a final score or rating given to the starting model. This may be in terms of a percentage score or regions of the starting model identified as inaccurate or problematic. A full report could then be generated. Variations on the above method may be used. For example, whilst the above embodiment has been described with reference to utilising the starting modelled data set u₀ (r, s) and the modelled data set after the first iteration of the FWI ui (r, s), other iterations could be used. For example, the modelled data sets after the first and second iterations (i.e. ui (r, s) and 112 (r, s)) could be used.

Alternatively, the next but one iteration could be used, e.g. the starting modelled data set and the second iteration could be used to amplify the divergence of the data.

Additionally, further iterations could be used and Δφ_η calculated for a plurality of iterations n. This could be used to highlight problem points in the data as the iterations proceed. Furthermore, the data sets provided in steps 104 and 106 need not be the modelled data sets and, instead, values of the residual variable could be provided instead.

A second embodiment of the present invention will now be described with reference to Figure 11. The second embodiment includes all of the method steps 100 to 118 above. Steps 100 to 118 above enable identification, qualification and quantification of mis- convergence which may occur from an inaccurate starting model. However, the second embodiment, once a starting model is identified as inaccurate for a particular FWI method, provides for an alternative technique to be implemented which may have an improved chance of accurate convergence.

Once it has been demonstrated that the starting model for a particular measured data set is inadequate, an alternative iterative method can be attempted. This is known as phase unwrapped FWI. This avoids the difficult-to -attain requirement of modelled data being within half a cycle of measured data for successful convergence.

If, in step 110, it was determined that further processing was necessary on the data, and, at step 112, phase unwrapping was carried out, the phase unwrapped residual data set can be used in in the inversion process itself to generate the model update. This model update can then be quality checked by itself and in addition to the original model used.

Step 120: Construct objective function

In the present case, the present invention iteratively updates a trial velocity to minimise an objective function based on some norm of the unwrapped phase such as the L-2 norm. Therefore an objective function is constructed based on equation 23) below:

22) This is equivalent to the phase component of the logarithmic approach could be used as disclosed in "Comparison of waveform inversion, part 1: Conventional Wavefield vs Logarithmic WavefieW Shin et al, Geophysical prospecting, 2007, 55, 449-464 and WO-A-2009/136708.

Step 122: Minimise objective function and update model

The model update direction derives from the gradient of equation 23) i.e. the partial derivative with respect to a point perturbation of the model m(x) at each position x

To describe this procedure it is sufficient to consider a single shot, located at point s, and single frequency in the objective function. Ultimately gradients from separate shots and freq will summed when forming the final model update.

£ = |∑φ'² (r,s)

r

The gradient of this objective function with respect to the model is: dE

dm(x')

To calculate this we note:

Using equation 25) gives: dq>'<

dm ₇(¾x = 2co½(x')Im(¾#G(x',,))

26)

We then arrive at:

₂₇^ g = 2co½( ')Im(M( ',s)w^*( ',s)) where: w(x',s) = 2^G (x',r)—— - 28) r ^u (^r>^s)

As set out above in relation to the conventional FWI method, this is the product of two wavefields:

i) an incident wavefield u emitted by a source at original source location s; and a back-propagated wavefield w which is emitted by a (multi-point) source located at receiver positions. 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. Step 124: Verify data?

At step 124 it is determined whether convergence criteria have been met. In this case, it is possible to perform a check using the check residual steps outlined in steps 114 to 118 to check whether the unwrapping of the phase has improved the convergence of the results. If it has not improved substantially, we may retune the parameters of the unwrapped step 112.

As an example of this, Figure 16a) shows the same data for (f)_n as for Figure 10a). However, Figure 16b) shows Δφ_η for a phase unwrapped iteration as described above. As shown (and particularly in comparison with Figure 10b), it can be seen from the circled area that both (f)_nand Δφ_η are the same sign, indicating that the data is converging in the correct direction. Therefore, it can be inferred that the starting model used is inadequate for accurate convergence utilising a conventional FWI method, but is adequate for use with a phase unwrapped method.

If the criteria have been met, then the method proceeds back to repeat steps 104 to 122 as discussed above to continue further iterations until a convergence criterion is met. Once the convergence criterion is met, the method proceeds to step 126.

Step 126: Finish When, at step 126, 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 involves 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.

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.

The method according to the present invention has numerous advantages over prior art methods. As described above, here we show how to validate a FWI approach enabling more robust and reliable convergence when using this technique. We can identify cases where conventional FWI has or will be sufficient for a given starting model which does not exist in the prior art.

Variations in the method of the above embodiments will be apparent to the skilled person. In addition, alternative or additional steps may be implemented with respect to the above embodiments. For example, if time domain data is provided, it is possible to window the measured and modelled data to improve the accuracy of the final result. An example of windowing is given below with respect to modelled data.

The data obtained in steps 104 and 106 may comprise modelled seismic trace data in the time domain, i.e. modelled pressure at the location of the measured data points as a function of position and time. In this case, a time window can be applied to the modelled data set in steps 104 and 106.

An example of this is shown in the set of figures comprising Figure 18. Figure 18 a) shows a graph illustrating an aspect of the modelled seismic trace data 150. This figure shows a plot of pressure (in arbitrary units) at the source position as a function of time (in arbitrary units, which may be seconds). The source position corresponds to that of the measured data set, such that each spatial data point in the measured data set (representing a discrete source position) has a corresponding data point in the modelled data set. The graph of Figure 18a) corresponds to a modelled estimate of the time evolution of pressure at a particular receiver number or region and corresponds to a time slice taken through a data set such as that shown in Figure 9b). As shown in Figure 12a), there is a peak at early times (early arrival pressure waves) with further oscillations at longer times. The inventor has found that the effect of applying a time-domain window to the modelled and measured seismic data sets prior to calculating the phase residual is operable to improve the results of the validation check.

A time-domain window is then to the modelled seismic trace data 150. This is done by applying a function 152 to the modelled data. A suitable function may be a Gaussian function 162 as shown in Figure 18b). The Gaussian function 152 has various parameters which can be specified by a user such as the time (τ) at which the function or window 152 is centered and the Full Width Half Maximum (FWHM) which specifies the spread of the function.

However, whilst a Gaussian operator has been specified as the function 152 in this embodiment, other suitable functions may be used with the present invention. A non-exhaustive list may comprise: exponential functions, step functions, delta functions, triangular functions, quadratic or linear functions, Hann, Hamming or cosine windows. In general, the same window (i.e. the same function) will be applied to the measured data obtained in step 100 for this embodiment.

Figure 18c) shows the function 152 superimposed on the original modelled data 150 prior to multiplication. In common with step 102, the effect of the function 152 as shown will be to enhance relatively the amplitude of early arrivals with respect to later arrivals.

The modelled data 150 is then combined with the Gaussian function 152, in this example by multiplication. Figure 18d) shows the resulting windowed modelled seismic data set 164 obtained from the multiplication. Whilst this step has been illustrated with respect to a single receiver position, it is to be understood that the windowing approach will, in general, be carried out on the whole, or a significant portion thereof, of the modelled seismic data set to generate a windowed modelled seismic data set 154 for all traces under consideration. A method according to a third embodiment of the present invention will now be described with reference to Figure 17. Figure 17 shows a flow diagram of an embodiment of the present invention. The third embodiment will be described with reference to a verification of the suitability of a subsurface starting model for use in interpreting a particular measured data set. Therefore, the general principle of the invention could be applied in a number of situations.

For example, the starting subsurface model and initial modelled data could be provided along with the measured data in order to provide guidance or verification of the likelihood of success of the model reaching an accurate convergence using that particular starting model following conventional FWI. Alternatively, the starting subsurface model, initial modelled data and a first iteration of the modelled data could be provided along with the measured data in order to provide guidance or verification of the likelihood of success of the model reaching an accurate convergence using that particular starting model following conventional FWI.

As a further alternative, an already-run convergence procedure could be critically analysed by reviewing the residuals at each iteration of the FWI method to determine the validity of the final result. As a further alternative, the method could be run alongside or within an FWI method to provide an indication as to the accuracy of the FWI process as the iterations are computed.

As a yet further alternative, two different starting models could be compared using the following analysis. This could be done by generating an initial modelled data set for each starting model and then comparing the results of the two to determine which starting model is most likely to give rise to accurate convergence.

Step 200: Obtain measured seismic data set

The starting point for the subsurface exploration analysis is a set of experimentally gathered seismic trace data. This may be gathered by an experimental arrangement such as the experimental set up shown and described with reference to Figure 1.

As shown in Figure 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 in a line (for a two dimensional measurement) or in a grid (for a three dimensional measurement). The position of the detectors 16 is well known from location tracking devices such as GPS devices. Additionally, the location of the source 12 is also well known by location tracking means.

The measured seismic data set may comprise multiple "shots" or source 12 emissions. An example of seismic trace data for a single shot is shown in Figure 4 a). The gathered data is essentially a plot of pressure as a function of receiver position (on the x-axis) and time (on the y- axis). This is because, in general, a detector such as a hydrophone measures the scalar pressure at its location. The seismic trace data comprises a plurality of measured data points. Each measured discrete data point has six associated location values - three spatial directions (x, y and z) for receiver (or detector) position (r), and three spatial dimensions (x, y, z) for source location (s), together with pressure magnitude data. The spatial coordinates for each discrete data point define its location in space. It also comprises one or more measurement parameters which denote the physical property being measured. In this embodiment, a single measurement parameter, pressure is utilised.

The measured seismic trace data comprises data in the time domain, i.e. measured pressure as a function of position and time. The data can be analysed and FWI performed on time-domain data, and this is intended to comprise an aspect of the present invention. The method now proceeds to step 202.

Step 202: Obtain starting model

At step 202, an initial starting model of the specified subsurface portion of the Earth is, optionally, provided, potentially by a third party. 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 invention is applicable to three dimensional approaches. Examples of typical starting models are shown in Figure 7.

The generated model consists of values of the physical parameter V_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. It is noted that this step is optional and, instead, a third party may supply merely the starting modelled data set as set out in step 204. The method then proceeds to step 204.

Step 204: Obtain starting modelled data set

In step 204, a starting data set is obtained. The actual generation of the data in this step is, as described above, not material to the present invention. The data may be merely provided or obtained from another source. Instead, the present invention relates to confirmation of the validity of the generated data. However, for a complete understanding of the approach used, the method to generate the starting modelled data set will be described.

In other words, the modelled data set corresponds discrete point to discrete point to the measured dataset. The modelled data set may be generated in the time domain in order to perform time domain analysis with respect to the measured data set obtained in step 200. The method now proceeds to step 206.

Step 206: Time Window Modelled Seismic Data Set

The data obtained in step 204 comprises modelled seismic trace data in the time domain, i.e. modelled pressure at the location of the measured data points as a function of position and time. Figure 18 a) shows a graph illustrating an aspect of the modelled seismic trace data 150. This Figure shows a plot of pressure (in arbitrary units) at the source position as a function of time (in arbitrary units, which may be seconds). The source position corresponds to that of the measured data set, such that each spatial data point in the measured data set (representing a discrete source position) has a corresponding data point in the modelled data set.

The graph of Figure 18a) corresponds to a modelled estimate of the time evolution of pressure at a particular receiver number or region and corresponds to a time slice taken through a data set such as that shown in Figure 4a). As shown in Figure 18a), there is a peak at early times (early arrival pressure waves) with further oscillations at longer times. The inventor has found that the effect of applying a time-domain window to the modelled and measured seismic data sets is operable to improve the results of the validation check.

Step 206 is operable to apply a time-domain window to the modelled seismic trace data 150. This is done by applying a function 152 to the modelled data. A suitable function may be a Gaussian function 152 as shown in Figure 18b). The Gaussian function 152 has various parameters which can be specified by a user such as the time (τ) at which the function or window 152 is centred and the Full Width Half Maximum (FWHM) which specifies the spread of the function.

However, whilst a Gaussian operator has been specified as the function 152 in this embodiment, other suitable functions may be used with the present invention. A non-exhaustive list may comprise: exponential functions, step functions, delta functions, triangular functions, quadratic, linear functions, Hann, Hamming or Cosine functions.

Figure 18c) shows the function 152 superimposed on the original modelled data 150 prior to multiplication. The effect of the function 152 as shown will be to enhance relatively the amplitude of early arrivals with respect to later arrivals.

The modelled data 150 is then combined with the Gaussian function 152, in this example by multiplication. Figure 18d) shows the resulting windowed modelled seismic data set 154 obtained from the multiplication. Whilst this step has been illustrated with respect to a single receiver position, it is to be understood that the windowing approach will, in general, be carried out on the whole, or a significant portion thereof, of the modelled seismic data set to generate a windowed modelled seismic data set 154 for all data points under consideration. The method then proceeds to step 208. Step 208: Transform modelled data set

The windowed modelled seismic trace data generated in step 206 may, as set out in step 206, include multiple frequencies. These must be resolved into single or groups of discrete frequencies so that direct analysis between the modelled and measured data can be performed. The different frequency components can be extracted through the use of a Fourier transform if so desired. 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.

One or more modelled data sets at one or more frequencies may be extracted. For example, a plurality of modelled data sets may be provided having frequency components of, for example, 4.5 Hz, 5 Hz and 5.5 Hz. These can be inverted individually or simultaneously.

The modelled data set, u (r, s), comprises data at one or more frequencies. These frequencies correspond to that of the measured data set d(r, s) (as will be described later) such that meaningful analysis can be performed. Several frequencies can be inverted simultaneously. The method now proceeds to step 210.

Step 210: Time Window measured seismic data set

The data obtained in step 200 comprises measured seismic trace data in the time domain, i.e. measured pressure as a function of position and time. A time window is then applied to the modelled data set in step 210. This is generally similar to step 206 above but is applied to the measured data set.

Figure 19 a) shows a graph illustrating an aspect of the raw trace data 160. This figure shows a plot of pressure (in arbitrary units) at the source position as a function of time (in arbitrary units, which may be seconds). This graph corresponds to the time evolution of pressure at a particular receiver number or region and corresponds to a time slice taken through a data set such as that shown in Figure 4b). As shown in Figure 19a), there is a peak at early times (early arrival pressure waves) with a complicated series of oscillations extending out to longer times due to later arrivals.

The inventor has found that the effect of applying a time-domain window to the modelled and measured seismic data sets is operable to improve the results of the validation check. Step 210 is operable to apply a time-domain window to the measured seismic trace data 160. This is done by applying a function 162 to the measured raw data. A suitable function may be a Gaussian function 162 as shown in Figure 19b). The Gaussian function 164 has various parameters which can be specified by a user such as the time (τ) at which the function or window 162 is centred and the Full Width Half Maximum (FWHM) which specifies the spread of the function. The value of τ is, as set out for step 206, calculated from the modelled data first arrival time or from starting model i.e. by ray tracing to compute first arrival travel times. The same window is then applied to the measured data. There is no need to pick first breaks in the measured data, which is a time-consuming and error-prone process.

However, whilst a Gaussian operator has been specified as the function 162 in this embodiment, other suitable functions may be used with the present invention. A non-exhaustive list may comprise: exponential functions, step functions, delta functions, triangular functions, quadratic, linear functions, Hann, Hamming or Cosine functions. In general, the same window (i.e. the same function) will be used as is applied to the modelled data in step 206. Figure 19c) shows the function 162 superimposed on the raw data 160 prior to multiplication.

As shown, the effect of the function 162 as shown will be to enhance relatively the amplitude of early arrivals with respect to later arrivals. Therefore, the application of the window will ultimately focus the method on the early arrival data and reduce proportionally the effect of later arrivals on the iterative solving of the model. The value of τ may be selected from the first arrival time derived by ray tracing in order to centre the window over the first arrivals. However, the window centre τ of the function 152 may alternatively be selected to be at later times to focus the FWI on later arrivals if desired. Also τ can simply be user-specified as an increasing function of offset.

The measured raw data 160 is then combined with the Gaussian function 162, in this example by multiplication. Figure 20d) shows the resulting windowed measured seismic data set 164 obtained from a multiplication of the raw measured data 160 and the Gaussian function 162. As previously described, the effect of the time windowing procedure is to enhance the early arrival times and dampen the later arrival times. Whilst this step has been illustrated with respect to a single receiver position, it is to be understood that the windowing approach will, in general, be carried out on the whole, or a significant portion thereof, of the measured seismic data set to generate a windowed measured seismic data set 164 over all traces. The method then proceeds to step 212. Step 212: Transform measured seismic data set

The windowed measured seismic trace data 164 comprises data in the time domain, i.e. measured pressure as a function of position and time. This data will contain multiple frequencies and, in order for data to be analysed against modelled data, it is necessary to extract data taken at discrete frequencies. Therefore, it is necessary to resolve the different frequency components within the seismic trace data to form one or more windowed measured data sets with the or each data set representing the windowed seismic traces at a selected frequency. In other words, a windowed measured data set is constructed which comprises the frequency components of the windowed trace data at discrete data points. The measured data set may contain one or more measurement parameters and one or more frequency components.

One or more measured data sets at one or more frequencies may be extracted. For example, a plurality of measured 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.

Alternatively, these frequencies can be inverted simultaneously. d(r,s) is the notation we give for the measured data set when it has one or more measured physical parameters and one or more frequency components. The method then proceeds to step 214. Step 214: Determine phase residual

As set out above, the measured seismic data is expressed as a function of receiver (r) and source (s) pair position. This is the function d(r,s). The modelled seismic data is also expressed as a function of receiver (r) and source (s) pair position. This, as set out above, is the function udr,s) for the starting modelled data set. The modelled and measured data sets u and d may comprise values of one or more physical parameters (e.g. pressure) taken at one or more frequencies.

In step 214, the phase residual variable is obtained for the modelled data set, i.e. between the data sets d(r,s) and uo(r,s). The principal value (between -π and + π) of the phase residual is denoted as φ and expressed in equation 29) for an iteration n:

The phase residual φη(*",«) can be determined for the data points to generate a phase residual data set. In that case, each spatial data point now has associated therewith a value of the phase residual φ_η(*",*) for an iteration n. Note that the starting modelled data set, udr,s) is generated prior to the first iteration (for which n =1). In the case of the starting modelled data set udr,s), φο(/*,«) can be calculated from equation 29). The method then proceeds to step 216.

Step 216: Identify Cycle-Skipped Data in Modelled Dataset

Once the phase residual data set is generated in step 214, the cycle skipped data therein can be identified. Figure 20 shows this analysis. Figures 20 a) and 20 b) illustrate a suitable and a poor starting model respectively. Figures 20 c) and 20 d) illustrate phase residual data sets for the models of Figures 20 a) and 20 b) respectively.

Φ„(Γ,«) is plotted as a function of source and receiver in Figures 20 c) and 20 d). As shown in Figure 20 c), φ_η(*",*) is plotted at 3.125Hz and varies smoothly in space. Where the values of φ_η(*",*) are negative, the modelled data is early relative to the observed data and where positive, the values of Φ„(Γ,«), the modelled data is late relative to the observed data. φ_η(/·,«) is close to zero at zero-offset and in the case of the starting model of Figure 20 a) remains well within the -π to + π range. Also of note is the symmetry in the zero-offset axis which occurs through source-receiver reciprocity.

In the case of the second starting model shown in Figure 20 b), φ_η(*",*) as plotted in Figure 20 d) is similarly smoothly varying but in general with a larger magnitude than with the second starting model. At a certain offset its value decreases sufficiently to reach -π. Instead of passing through -π, it abruptly switches to +π before decreasing again. This 2π transition in φ_η(Γ,«) forms a continuous boundary through the data. It marks where the phase difference between the modelled and observed data exceeds half a cycle. In this case it marks where the modelled data is more than half a cycle early relative to the observed data. This shows how the well-known 'cycle-skipping' relationship between predicted and observed data manifests itself in the low frequency phase residual.

The presence of this cycle-skipped section, beyond the 2π transition when moving away from zero-offset, provides a technical indicator from the data mismatch that the starting model of Figure 20 b) is less accurate than the starting model of Figure 20 a).

The cycle skipped region can then be identified numerically or graphically for further computational processing and analysis. The cycle skipped region can be identified by the cycle-skip boundary, with data beyond the cycle skip boundary when moving away from zero-offset being cycle-skipped. In Fig 20 d), the cycle-skipped data is the data in the top right i.e. on the side with a '+ π'. Similarly, the cycle skipped data is the data beyond the boundary in the bottom left. In Fig 20 c) there is no cycle-skipped data. As set out above, the cycle-skipped boundaries can be identified numerically and

computationally and the cycle-skipped regions identified therefrom. The method then proceeds to step 218.

Step 218: Evaluate cycle-skipped data

At step 218, the cycle-skipped data is evaluated. This may be done in a number of suitable ways. For example, the proportion of cycle-skipped data within the phase residual data set may be evaluated. This may be expressed as a percentage score wherein the lower the percentage score, the higher the accuracy of the starting model and the greater the probability of accurate convergence. A parameter may be assigned. For example, a starting model may be considered to be suitable if less than 20% of the data is cycle-skipped at 3 Hz for example. However, other examples may be used, e.g. 30% or 40% depending upon the application. Other methods of numerical analysis may be used. The distribution of cycle-skipped data may be utilised.

Step 220: Complete analysis

At step 220, the analysis can then be completed and a final score or rating given to the starting model. This may be in terms of a percentage score or regions of the starting model identified as inaccurate or problematic. A full report could then be generated.

Variations on the above method may be used. For example, whilst the above embodiment has been described with reference to utilising the starting modelled data set u₀ (r, s), other iterations could be used. For example, the modelled data sets after the first and second iterations (i.e. ui (r, s) and 112 (r, s)) could be used. Alternatively, the next but one iteration could be used, e.g. the second iteration could be used.

A fourth embodiment of the present invention will now be described with reference to Figure 21. The fourth embodiment includes all of the method steps 200 to 220 above. Steps 200 to 220 above enable identification, qualification and quantification of mis-convergence which may occur from an inaccurate starting model. An analysis of the phase residual can then be performed to determine the accuracy of the starting model. In the fourth embodiment, this process is enhanced by carrying out a further iteration of the model to see how the phase residual changes.

Step 222: Obtain farther modelled data set

The data obtained in step 204 comprises modelled seismic trace data in the time domain (i.e. modelled pressure at the location of the measured data points as a function of position and time).

The further modelled data set required in step 222 is obtained after one or more iterations of FWI has been carried out on the modelled data of step 204. As set out in relation to step 204, the actual generation of the data in this step is, as described above, not material to the present invention. The data may be merely provided or obtained from another source. Instead, the present invention relates to confirmation of the validity of the generated data. However, for a complete understanding of the approach used, the method to generate the further modelled data set will be described in the context of a single iteration.

In the case of an FWI method, a trial velocity is updated to minimise an objective function measuring the mismatch between modelled and recorded data. The objective function and the steepest descent direction used to update the model in conventional FWI is described in equations 1) to 12). This can be done in the time or frequency domain, although time domain data is required for the analysis. A further modelled data set in the time domain obtained from an iteration of the method described above is then provided.

Instead of the example given above, the further modelled data set may comprise a modelled data set from starting, or subsequent, iteration of a calculation using a different model. This would then allow the models to be compared.

Repeat Steps 206, 208, 214 - 218

The method then proceeds, using the further modelled data set, through steps 206 and 208 to obtain ui (r,s). Steps 210 and 212, having already been carried out, may be omitted at this stage. The method then proceeds to analyse the further modelled data set ui (r,s) at steps 214, 216 and 218. However, in this embodiment, the evaluation at step 218 comprises comparing the two phase residual sets.

Step 218: Evaluate cycle-skip change

As described above, the present method provides a tool to assess whether a given starting model is adequate for FWI on a particular seismic dataset i.e. whether it will converge successfully or not under conventional FWI. This is done by identifying cycle-skipped data in the data set and analysing the change between iterations.

For example, if the proportion of cycle-skipped data in the further modelled data set increases when compared to the starting modelled data set, the model will not result in an accurate convergence and should be rejected. Examples of this are shown in Figures 22 and 23. Figure 22 shows a schematic diagram illustrating the configuration of Figures 23 a) to f). In Figure 22, a plurality of lines of receivers is shown in a field of sources. The shooting direction for the sources is shown as left to right. A subterranean gas cloud is highlighted at approximately the centre of the figure. Finally, the plots shown in Figures 23 c) to f) relate to the receiver position shown in Figure 22.

Two starting models for the configuration shown in Figure 22 are shown in Figures 23 a) and 23 b). Figure 23 a) shows a good starting model and Figure 23 b) shows a poorer starting model.

Figures 23 c) and d) illustrate the phase residual for the starting models of Figures 23 a) and b) respectively. As shown, the cycle skipped region is greater in Figure 23 d). Next, an iteration of FWI method outlined above is performed on the data, and steps 206 to 208 and 214 to 216 are carried out on the first iteration of the starting model. This process generates the phase residual plots shown in Figures 23 e) and f).

By comparing the phase residual plots of Figure 23 c) and 23 e), it can be seen that, for the better starting model of Figure 23 a), between the starting model and the first iteration, the cycle skipped data region has shrunk. However, with reference to Figure 23 d) and 23 f), it can be seen that, for the poorer starting model of Figure 23 b), the cycle skipped region has grown. Therefore, it can be concluded that further iterations utilising the starting model of Figure 23 a) are likely to converge to global minimum and produce a successful minimisation. However, the model of Figure 23 b) is likely to converge to a cycle-skipped local minimum, providing an erroneous result.

Alternatively, to determine how the data is changing with each iteration, the change in the phase residual variable between iterations (or, as described here, between the starting data set and the first iteration) can be determined and analysed.

Thus, in step 218 the change phase residual Δφ can be defined. The change phase residual Δφ is thus defined in equation 30) for the case where the starting modelled data set uo(r,s) and the first iteration (further) modelled data set ui(r,s) are used:

30) Δφ = phase ((u₀d*) (uid*)*)= p ase(u₀ui *) or, in a more general form (with dependency on source and receiver locations shown) for an iteration «. ·

31 ) Δφ_η (r,s)= phase (u_n (r,s) u_n+ ](r,s) *) Equations 30) and 31 ) describe the change in phase residual caused by a model update. Δφ_η

(r,s) is a wrapped quantity lying between - π and + π. Further it is not explicitly dependent on d so is noise-free. Alternatively, Δφ_η (r,s) can be calculated directly from u₀(r,s) and ui (r,s) and so it is not explicitly necessary to calculate or obtain the phase residual for iteration 1 ^l(r,s)) directly, although this may be done if required. Δφ_η (r,s) may be obtained at one or many frequencies. These frequencies need not be those used in the inversion procedure.

An analysis of the change phase residual can provide further detailed information to vindicate or reject a particular starting model. This is an advantage of the present method which enables the divergence or mis-convergence of the data to be observed as a function of source and receiver position. This is in contrast to any known methods which may simply sum the errors over the whole dataset.

Step 224: Finish When, at step 224, 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 involves 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.

In addition, the above first to fourth embodiments enable identification of outlying (i.e. cycle- skipped data, or close to cycle-skipped data) within the data set. This could be implemented in an alternative method for improving the accuracy of an FWI process as follows. A fifth embodiment of the present invention is illustrated with respect to Figure 24. Method steps 200 to 216 of the fourth embodiment are identical to steps 200 to 216 described previously. Step 226: Unwrap phase residual

The phase residual (or phase mismatch) is calculated in step 214 as set out above. The value calculated is the wrapped phase residual. As shown in Figure 9c), the phase residual is zero at the source location (or zero-offset point) and oscillates within +π to -π until the receiver number increases beyond about 320. At this point, the phase residual increases beyond +π and wraps to -π. The phase residual increases beyond +π again at around receiver number 400 and wraps around to -π.

However, additional advantages can be achieved by unwrapping the phase residual to obtain further information about outlying data. Therefore, step 226 may optionally be performed to provide further information for analysis.

In dataset where we have such occurrences of 2 π wrap around we phase unwrap by adding multiples of 2 π to the phase residual prior to analysis (Figure 12 above. In real-world data, Φ(Γ,«) is a function of three spatial dimensions for both receiver and source. Therefore, Φ(Γ,«) will be a function of r_x, r_y, r_z, s_x, s_y, s_z and frequency / However, for clarity, consider a two dimensional case where Φ(Γ,«) is a function of just one spatial dimensions for each of the receiver and source. Therefore, Φ(Γ,«) is a function of r_x and s_x. The phase residual can then be unwrapped in accordance with equations 16) to 18) outlined above and as described with reference to steps 110 and 112. Alternatively, the unwrapping may be described in terms of degrees, in which case for each multiple of π, there would be a multiple of 180 degrees. The method then proceeds to step 228.

Step 228: Identify Outlying Data in Modelled Dataset

Once the phase residual data set is generated in step 214, the outlying data therein can be identified. This may take two forms in this embodiment. If step 226 is not carried out, the method proceeds in accordance with step 216 and cycle skipped data is identified.

However, if step 226 is implemented, the unwrapped phase residual is obtained. This can then be analysed. In this example, it is possible to identify not only cycle-skipped data, but data which may enter a cycle-skipped region should an iteration proceed. This is known as outlying data and this term is intended to cover both cycle-skipped data and data in a region close to cycle skipping. This is possible in this embodiment because unwrapped phase provides the magnitude of the data as opposed to just identifying cycle skipped data.

The outlying data region can then be identified numerically or graphically for further computational processing and analysis. The outlying data region can be identified, in this

embodiment, by any data having a magnitude in excess of π when unwrapped. This is cycle skipped data. Alternatively, a higher threshold could be implemented and, for example, any data having a magnitude greater than 0.88 π when unwrapped may be considered outlying data. We may set the threshold above pi, say 1. lpi. Alternatively, they may be expressed in degrees, so 160 degrees or greater, or 200 degrees or greater, for example. As set out above, the outlying regions can be identified numerically and computationally and the cycle-skipped regions identified therefrom. The method then proceeds to step 230.

Step 230: Weight outlying data

As previously described, at step 228 the outlying (and, in certain embodiments, cycle- skipped) data can be explicitly identified. The outlying data in the phase diagrams has a spatial coordinate associated therewith. In other words, each data point has a receiver location r and a source location s. Therefore, this coordinate location can be used to identify the outlying data within the measured data set. It can then be deleted from the measured data set.

The outlying data of a particular data point may be given a value between 0 and 1 in dependence upon the magnitude of the unwrapped phase residual for that discrete data point.

For example, a discrete data point may be given a value of 0 if said data point has a value of 180° (π) or more. A more stringent criterion would be to identify data close to being cycle skipped but not actually cycle skipped; for example, a discrete data point may be given a value of 0 if said data point has a value of 160° (0.88π) or more.

In one embodiment, if a data point has been assigned a weight of 0, the data point is ignored in subsequent steps. Alternatively, the data point could be deleted from the measured and the updated modelled seismic data sets. The method proceeds to step 232. Step 232: Generate new starting model

Once the measured data set has been amended to weight the outlying data, then the modified measured data set can be used as a basis for a further iteration of FWI to generate a new starting model. This is advantageous since there is no requirement to modify the original starting model in any way. By simply weighting the outlying data and the re-running an initial iteration of an FWI process to generate a new starting model, a subsequently run FWI process may more reliably converge to a global minimum without the cycle skipped or outlying data points in the data set. Where a weighting of 0 has been given to the data point, this data point is to be ignored.

Alternatively it may be deleted. Where a discrete data point has a non-zero value between 0 and 1 , the data point is given a weighting to the minimisation based on the non-zero value, with a value of 1 being given full weight. The FWI process involves minimising an objective function defined to measure the difference between the measured and modelled dataset. The objective function for frequency domain inversion is given above in equation (1). It will have the outlying data down -weighted to simply deleted.

Step 234: Utilise new starting model

In step 234, the new starting model generated in step 232 is utilised. This may involve repeating the entire quality control process, starting at step 200 with the original measured data set and the new starting model, or the new starting model may simply be used as the basis for a complete FWI minimisation process.

Whilst the fifth embodiment has been described in the context of the third and fourth embodiments above, it is to be understood that the fifth embodiment is also applicable to the first and second embodiments above and can be used with any method that is operable to generate a phase residual data set to enable identification of cycle skipped data.

A method according to a sixth embodiment of the present invention will now be described with reference to Figure 24. Figure 24 shows a flow diagram of the sixth embodiment of the present invention. The sixth embodiment will be described with reference to a verification of the suitability of a subsurface starting model for use in interpreting a particular measured data set. Therefore, the general principle of the invention could be applied in a number of situations.

As a further alternative, an already-run convergence procedure could be critically analysed by reviewing the residuals at each iteration of the FWI method to determine the validity of the final result. As a further alternative, the method could be run alongside or within an FWI method to provide an indication as to the accuracy of the FWI process as the iterations are computed. As a yet further alternative, two different starting models could be compared using the following analysis. This could be done by generating an initial modelled data set for each starting model and then comparing the results of the two to determine which starting model is most likely to give rise to accurate convergence. Step 300: Obtain measured seismic data set

The starting point for the subsurface exploration analysis is a set of experimentally gathered seismic trace data. This may be gathered by an experimental arrangement such as the experimental set up shown and described with reference to Figure 1. In this regard, this step corresponds to step 200 of the third and fourth embodiments.

However, it is possible to measure other parameters using appropriate technology; for example, to measure the velocity or displacements of particles in three spatial directions in addition to the pressure. The present invention is applicable to the measurement of these additional variables. 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 invention. The present invention simply requires a real-world measured seismic data set upon which analysis can be performed to facilitate subsurface exploration of a portion of the Earth. This may be supplied by a third party.

The measured seismic trace data comprises data in the time domain, i.e. measured pressure as a function of position and time. The data can be analysed and FWI performed on time-domain data, and this is intended to comprise an aspect of the present invention. The method now proceeds to step 302.

Step 302: Obtain starting model

At step 302, an initial starting model of the specified subsurface portion of the Earth is, optionally, provided, potentially by a third party. 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 invention is applicable to three dimensional approaches. Examples of typical starting models are shown in Figure 7.

The generated model consists of values of the physical parameter V_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. It is noted that this step is optional and, instead, a third party may supply merely the starting modelled data set as set out in step 304. The method then proceeds to step 304.

Step 304: Obtain starting modelled data set

In step 304, a starting data set is obtained. The actual generation of the data in this step is, as described above, not material to the present invention. The data may be merely provided or obtained from another source. Instead, the present invention relates to confirmation of the validity of the generated data. However, for a complete understanding of the approach used, the method to generate the starting modelled data set will be described. In order to proceed with the FWI method in order to model the subsurface region under investigation, it is necessary to utilise modelled data which corresponds to the same source-receiver location data position as the actual measured trace data so that the modelled and measured data can be compared. In other words, the modelled data set corresponds discrete point to discrete point to the measured dataset. The modelled data set may be generated in the time domain in order to perform time domain analysis with respect to the measured data set obtained in step 300. The method now proceeds to step 306. Step 306: Time Window Modelled Seismic Data Set

The data obtained in step 304 comprises modelled seismic trace data in the time domain, i.e. modelled pressure at the location of the measured data points as a function of position and time. This data may, in this embodiment, be directly Fourier transformed and so step 306 is optional. In which case, the method proceeds directly from step 304 to step 308.

Referring again to Figure 18 a), this figure shows a graph illustrating an aspect of the modelled seismic trace data 150. This Figure shows a plot of pressure (in arbitrary units) at the source position as a function of time (in arbitrary units, which may be seconds). The source position corresponds to that of the measured data set, such that each spatial data point in the measured data set (representing a discrete source position) has a corresponding data point in the modelled data set.

Step 306 is operable to apply a time-domain window to the modelled seismic trace data 150. This is done by applying a function 152 to the modelled data in the manner of that disclosed in the third and fourth embodiments and as set out in step 206. A suitable function may be a Gaussian function 152 as shown in Figure 18b). The Gaussian function 152 has various parameters which can be specified by a user such as the time (τ) at which the function or window 152 is centred and the Full Width Half Maximum (FWHM) which specifies the spread of the function.

Figure 18c) shows the function 152 superimposed on the original modelled data 150 prior to multiplication. The effect of the function 152 as shown will be to enhance relatively the amplitude of early arrivals with respect to later arrivals. The modelled data 150 is then combined with the Gaussian function 152, in this example by multiplication. Figure 18d) shows the resulting windowed modelled seismic data set 154 obtained from the multiplication.

Whilst this step has been illustrated with respect to a single receiver position, it is to be understood that the windowing approach will, in general, be carried out on the whole, or some portion thereof, of the modelled seismic data set to generate a windowed modelled seismic data set 154 for all data points under consideration.

For example, in the dataset shown in Figure 22, the method as described would be carried out for a handful of the total 1440 receivers. In addition, either before or after windowing, the time- domain modelled dataset might be shifted in time. The shift parameter may be a pre-determined increasing function of offset. The method then proceeds to step 308. Step 308: Transform modelled data set

The modelled seismic trace data generated in step 304, or the windowed data generated in step 306, may include multiple frequencies. These must be resolved into single or groups of discrete frequencies so that direct analysis between the modelled and measured data can be performed.

The different frequency components can be extracted through the use of a Fourier transform if so desired. 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. One or more modelled data sets at one or more frequencies may be extracted. For example, a plurality of modelled data sets may be provided having frequency components of, for example, 4.5 Hz, 5 Hz and 5.5 Hz. These can be inverted individually or simultaneously.

The modelled data set, u (r, s), comprises data at one or more frequencies. These frequencies correspond to that of the measured data set d(r, s) (as will be described later) such that meaningful analysis can be performed. Several frequencies can be inverted simultaneously. The method now proceeds to step 310.

Step 310: Time Window measured seismic data set

The data obtained in step 300 comprises measured seismic trace data in the time domain, i.e. measured pressure as a function of position and time. A time window may then applied to the modelled data set in step 310. This is generally similar to step 306 above but is applied to the measured data set. Again, step 310 is optional in this embodiment. Step 310 will not be carried out if step 306 is omitted as set out above. Either both the modelled and measured data sets are windowed, or neither modelled nor measured data sets are windowed.

Figure 19 a) shows a graph illustrating an aspect of the raw trace data 160. This figure shows a plot of pressure (in arbitrary units) at the source position as a function of time (in arbitrary units, which may be seconds). This graph corresponds to the time evolution of pressure at a particular receiver number or region and corresponds to a time slice taken through a data set such as that shown in Figure 4b). As shown in Figure 19a), there is a peak at early times (early arrival pressure waves) with a complicated series of oscillations extending out to longer times due to later arrivals. The inventor has found that the effect of applying a time-domain window to the modelled and measured seismic data sets is operable to improve the results of the validation check.

Step 310 is operable to apply a time-domain window to the measured seismic trace data 160. This is done by applying a function 162 to the measured raw data. A suitable function may be a Gaussian function 162 as shown in Figure 19b). The Gaussian function 164 has various parameters which can be specified by a user such as the time (τ) at which the function or window 162 is centred and the Full Width Half Maximum (FWHM) which specifies the spread of the function. The value of τ is, as set out for step 306, calculated from the modelled data first arrival time or from starting model i.e. by ray tracing to compute first arrival travel times. The same window is then applied to the measured data. There is no need to pick first breaks in the measured data, which is a time- consuming and error-prone process.

However, whilst a Gaussian operator has been specified as the function 162 in this embodiment, other suitable functions may be used with the present invention. A non-exhaustive list may comprise: exponential functions, step functions, delta functions, triangular functions, quadratic, linear functions, Hann, Hamming or Cosine functions. In general, the same window (i.e. the same function) will be used as is applied to the modelled data in step 306.

Figure 19c) shows the function 162 superimposed on the raw data 160 prior to multiplication. As shown, the effect of the function 162 as shown will be to enhance relatively the amplitude of early arrivals with respect to later arrivals. Therefore, the application of the window will ultimately focus the method on the early arrival data and reduce proportionally the effect of later arrivals on the iterative solving of the model. The value of τ may be selected from the first arrival time derived by ray tracing in order to centre the window over the first arrivals. However, the window centre τ of the function 152 may alternatively be selected to be at later times to focus the FWI on later arrivals if desired. Also τ can simply be user-specified as an increasing function of offset.

The measured raw data 160 is then combined with the Gaussian function 162, in this example by multiplication. Figure 20d) shows the resulting windowed measured seismic data set 164 obtained from a multiplication of the raw measured data 160 and the Gaussian function 162. As previously described, the effect of the time windowing procedure is to enhance the early arrival times and dampen the later arrival times. In addition, either before or after windowing, the time-domain modelled dataset might be shifted in time. The shift parameter may be a pre-determined increasing function of offset. The shift of the measured data will be equal to that applied to the modelled data, if any.

Whilst this step has been illustrated with respect to a single receiver position, it is to be understood that the windowing approach will, in general, be carried out on the whole, or a significant portion thereof, of the measured seismic data set to generate a windowed measured seismic data set 164 over all traces. The method then proceeds to step 312.

Step 312: Transform measured seismic data set

The windowed measured seismic trace data 164 comprises data in the time domain, i.e. measured pressure as a function of position and time. This data will contain multiple frequencies and, in order for data to be analysed against modelled data, it is necessary to extract data taken at discrete frequencies. Therefore, it is necessary to resolve the different frequency components within the seismic trace data to form one or more windowed measured data sets with the or each data set representing the windowed seismic traces at a selected frequency.

In other words, a windowed measured data set is constructed which comprises the frequency components of the windowed trace data at discrete data points. The measured data set may contain one or more measurement parameters and one or more frequency components.

One or more measured data sets at one or more frequencies may be extracted. For example, a plurality of measured 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. Alternatively, these frequencies can be inverted simultaneously. d(r,s) is the notation we give for the measured data set when it has one or more measured physical parameters and one or more frequency components. The method then proceeds to step 314.

Step 314: Obtain farther modelled data set

The data obtained in step 304 comprises modelled seismic trace data in the time domain (i.e. modelled pressure at the location of the measured data points as a function of position and time). It is, in one embodiment, possible to perform meaningful analysis on a single modelled data set and a single measured data set as will be described below.

However, in one embodiment, further modelled data sets may be provided in order to determine the convergence of an iterative procedure based on a given starting model. Alternatively or additionally, a modelled data set using a different starting model may be generated at this stage. In this example, the further modelled data set required in step 314 is obtained after one or more iterations of FWI has been carried out on the modelled data of step 304. As set out in relation to step 304, the actual generation of the data in this step is, as described above, not material to the present invention. The data may be merely provided or obtained from another source. Instead, the present invention relates to confirmation of the validity of the generated data. However, for a complete understanding of the approach used, the method to generate the further modelled data set will be described in the context of a single iteration.

In the case of an FWI method, a trial velocity is updated to minimise an objective function measuring the mismatch between modelled and recorded data. The objective function and the steepest descent direction used to update the model in conventional FWI is described in equations 1) to 12). This can be done in the time or frequency domain, although time domain data is required for the analysis. A further modelled data set in the time domain obtained from an iteration of the method described above is then provided. Instead of the example given above, the further modelled data set may comprise a modelled data set from starting, or subsequent, iteration of a calculation using a different model. This would then allow the models to be compared.

Repeat at least Step 308

The method then proceeds, using the further modelled data set, through step 308, or through steps 306 and 308 (depending upon whether windowing is utilised) to obtain, in this example, ui (r,s). Steps 310 and 312, having already been carried out, may be omitted at this stage.

Step 316: Plot phase of measured data set

As set out above, the measured seismic data is expressed as a function of receiver (r) and source (s) pair position. This is the function d(r,s). Subsequent to the transforming of the measured seismic data in step 312 (or the transforming of the windowed measured seismic data if step 310 is carried out), the phase of the function d(r,s) for particular frequencies can be plotted. For the example of the system of Figure 22 described earlier, such a plot of the phase of the transformed (or windowed and transformed) measured data d(r,s) at a frequency of 3 Hz is shown in Figure 26 a). Figure 26 a) is plotted as a function of source x- and y-position for the receiver shown in Figure 22. Note that the phase as plotted is the wrapped, or principal value (i.e. between -π and + π) of the phase. As shown, the -π /+ π transitions in the phase form rings which encircle the receiver position. The method proceeds to step 318.

Step 318: Plot phase of modelled data set(s)

As set out above, the modelled seismic data is expressed as a function of receiver (r) and source (s) pair position. For the starting modelled data set, this is the function udr,s). However, whilst this is given as an example, any iteration of the modelled data set may be used, where u_n (r,s) denotes the data set and n is the iteration number. Nevertheless, the following examples describe the use of at least the starting model data set ¾o(f,s).

Subsequent to the transforming of the modelled seismic data in step 312 (or the transforming of the windowed measured seismic data if step 310 is carried out), the phase of the function uo(r,s) for particular frequencies can be plotted. For the example of the system of Figure 22 described earlier, such a plot of the phase of the transformed (or windowed and transformed) starting modelled data uo (r,s) at a frequency of 3 Hz is shown in Figure 26 c) and d) for the starting models shown in Figure 26 b) and 26 c) respectively. Note that Figures 26 b) and 26 c) show the same data as Figures 23 a) and 23 b).

Figure 26 d) and e) are plotted as a function of source x- and y-position for the receiver shown in Figure 22. The same scale and coordinate system is used as for the measured data plotted in Figure 26 a). Note that again the phase as plotted is the wrapped, or principal value (i.e. between -π and + π) of the phase. As shown, the -π /+ π transitions in the phase form rings which encircle the receiver position. The method then proceeds to step 320.

Step 320: Analyse data sets

Once the measured and modelled data sets are plotted, they can be analysed. The following example focusses on the effect of a respective iteration of the starting modelled data set or sets to determine whether convergence is occurring to a correct minimum. However, the principles of the analysis are applicable to data sets originating from different starting models or different iterations of the same starting model. Numerous techniques could be utilised to analyse the data. In this example, the relative positioning and movement of the rings of phase variation within the measured and modelled data sets is analysed to provide an indication of the likely convergence of the FWI process to a global minimum. Figures 27 a) to d) illustrate measured and modelled data sets for the starting model shown in Figures 23 a) and 25 b). Figures 28 a) to d) illustrate measured and modelled data sets for the starting model shown in Figures 23 b) and 26 c).

With respect to the measured data set, increased movement and/or misalignment of the phase rings of the modelled data with increasing iterations of the modelled data sets indicates that the modelled data is converging to an incorrect minimum. In other words, movement, divergence or spreading of the ring patterns indicates that convergence is occurring an incorrect, cycle skipped minimum. This is the case where the modelled data is more than half a cycle early relative to the observed data.

Consider first Figures 27 a) to d). Figure 27 a) shows the phase of the measured data set d(r,s) as shown in Figure 26 a). Figure 27 b) shows the starting modelled data set ¾o(f,s). Figure 27 c) shows the modelled data set W (r,s) after one iteration. Finally, Figure 27 d) shows the modelled data set ¾ os(r,s) taken after 108 iterations.

In Figures 27 a) to d), the fourth and sixth phase rings (as counted radiating out from the receiver position) are highlighted. By a comparison of Figures 27 b) and c), it can be seen that, after only one iteration, the rings are converging to the same position, spacing and configuration as the rings of the measured data set shown in Figure 27 a). This follows for the next 107 iterations, leading to the plot shown in Figure 27 d). Consequently, the analysis shows that the starting model of Figure 26 b) is an accurate starting model which is leading to an accurate convergence and that cycle skipping is avoided or minimised.

Consider, in contrast, Figures 28 a) to d). Figure 28 a) shows the phase of the measured data set d(r,s) as shown in Figures 26 a) and 27 a). Figure 28 b) shows the starting modelled data set uo(r,s) of the model of Figure 26 c). Figure 28 c) shows the modelled data set Wi(r,s) after one iteration. Finally, Figure 28 d) shows the modelled data set w,/os(r,s) taken after 108 iterations.

In Figures 28 a) to d), the fourth and sixth phase rings (as counted radiating out from the receiver position) are highlighted. By a comparison of Figures 28 b) and c), it can be seen that, after only one iteration, the rings are diverging in spacing and configuration when compared to the rings of the measured data set shown in Figure 28 a). This follows for the next 107 iterations, leading to the plot shown in Figure 28d). In the plot of Figure 28 d), it can be seen that an entire ring has disappeared and that the data is significantly different from the measured plot.

Consequently, the analysis shows that the starting model of Figure 26 c) is a poor starting model which is leading to mis-convergence of the iterative procedure to a cycle skipped local minimum. In the method of the present invention, it is possible to determine this from one iteration of FWI which saves the need to proceed with the erroneous inversion. In the examples given, the inversion was run to completion and the final model validates this finding.

Numerous approaches may be taken to analyse the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set. In this embodiment, the spatial location of phase distribution features within the measured and modelled phase data sets with respect to the source or receiver position is compared to determine the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set.

The radial spacing from the source or receiver location of a phase distribution feature within a selected ring-like feature in said modelled phase data set is compared with the radial spacing from the source or receiver location of a corresponding phase distribution feature in a corresponding ring-like feature in said measured phase data set.

By way of a qualitative assessment, if a phase distribution feature within an n^th ring-like feature in said modelled phase data set is closest to a corresponding phase distribution feature in an m^th (where m≠n) ring-like feature in said measured phase data set, then it is likely that the data, in subsequent iterations, will converge to a cycle-skipped local minimum. The data should, therefore, be rejected.

However, if, in a radial direction, a phase distribution feature within an n^th ring-like feature in said modelled phase data set is closest to a corresponding phase distribution feature in an n^th ring-like feature in said measured phase data set, then it is likely that the data will converge towards a non- cycle skipped minimum and further analysis of the model could be done. If, in a radial direction, a phase distribution feature within an n^th ring-like feature is directly adjacent a corresponding phase distribution feature in an n^th ring-like feature in said measured phase data set, the data could be directly validated as an accurate subsurface starting model. Whilst, as described above, the phase distribution features comprise a plurality of concentric ring-like features, this need not be so. Other phase distribution features may be used to interpret the data.

Step 322: Complete analysis

At step 322, the analysis can then be completed and a final score or rating given to the starting model. This may be in terms of a percentage score or regions of the starting model identified as inaccurate or problematic. A full report could then be generated. Step 324: Finish

When, at step 324, 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 involves 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. 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. Variations on the above method may be used. For example, whilst the above embodiment has been described with reference to utilising the starting modelled data set u₀ (r, s), other iterations could be used. For example, the modelled data sets after the first and second iterations (i.e. ui (r, s) and 112 (r, s)) could be used. Alternatively, the next but one iteration could be used, e.g. the second iteration could be used.

A seventh embodiment of the present invention will now be described with reference to Figure 29. The seventh embodiment includes at least some of the method steps 300 - 324 described above and these will not be described here further. Steps 300 to 324 above enable identification, qualification and quantification of mis-convergence which may occur from an inaccurate starting model. An analysis of the phase residual can then be performed to determine the accuracy of the starting model. In the sixth embodiment, this process is enhanced by optimising the windowing technique (optional in the above embodiment) to optimise the spatial consistency of the data sets. Steps in common with the above embodiment will not be described again here and share the same reference numerals.

Step 326: Spatial consistency acceptable?

Consider the example of Figures 30 a) and b). Figure 30 a) shows measured seismic data in the time domain. Figure 30 b) shows windowed measured seismic data in the time domain windowed in step 310 with a Gaussian window having a full-width-half-maximum of 0.8s. This window is applied to both the measured and modelled data sets in steps 310 and 306 respectively.

The windowed measured data sets and windowed modelled data sets are then Fourier transformed at 3Hz in steps 312 and 308 respectively. Once the data has been windowed and Fourier transformed, it can then be plotted in steps 316 and 318.

Next, the data is analysed for spatial consistency. The 360 degree transitions in the phase rings encircle the receiver position. Spatial consistency is, in one embodiment, characterised by the closing of the phase rings. This is illustrated by a comparison of Figure 31 a) with Figure 31 b). Figure 31 a) shows unwindowed modelled data taken at 3 Hz. Figure 31 b) shows the same data, but windowed with the Gaussian window described above.

As shown, in the windowed data of Figure 31 b) the rings are closed, signifying spatial consistency. Put another way, if the data is spatially consistent, then any path points A and B in Figure 31 b) will cross exactly one 360 degree transition.

If, in step 326 it is determined that the modelled data and corresponding measured data is spatially consistent, then the method proceeds to step 320 for further analysis of the divergence or convergence of the modelled data. If, however, the spatial consistency is insufficient, then the method proceeds to step 328.

Step 328: Modify window parameter

In step 328, the window parameter is modified to attempt to achieve greater spatial consistency. This may be done by, for example, narrowing the width of the Gaussian window from the example of 0.8 s given above or shifting the window.

However, any appropriate change to the window used may be performed. For example, the window may be narrowed or widened, shifted, or a different window may be used. The method then proceeds back to step 306 where the measured and modelled data sets are then windowed with the new window parameters and Fourier transformed prior to further spatial consistency analysis in step 326. The method according to the present invention has numerous advantages over prior art methods. As described above, here we show how to validate a FWI approach enabling more robust and reliable convergence when using this technique. We can identify cases where conventional FWI has or will be sufficient for a given starting model which does not exist in the prior art.

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 invention. In addition, alternative or additional steps may be implemented with respect to the present invention. For example, whilst the present invention has been described with reference to a model which involves solving the acoustic wave equation, the present invention is equally applicable to the models which involve solving the elastic wave equation.

Further, different windowing techniques may be used. The modelled data, the field data, and the further iterations or variations of the modelled data are all windowed in the above-described embodiments. However, these parameters may, in alternative arrangements, use the same window or different windows, where the windows differ in width, centre or function.

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 utilised and processed in the present invention.

Embodiments of the present invention 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 invention 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 invention.

Claims

1. A method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising:

a) obtaining a measured phase data set derived from a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points as a function of source and receiver location, and the measured phase data set comprising the phase component of the measured values taken at one or more selected frequencies and for said at least one physical parameter;

b) obtaining a modelled phase data set generated from a subsurface starting model, the modelled phase data set comprising the phase component of modelled values taken at the same plurality of discrete data points, at said one or more selected frequencies and for said at least one physical parameter;

c) generating, for at least one source or receiver position of the measured phase data set, a measured phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location;

d) generating, for said at least one source or receiver position of the modelled phase data set, a modelled phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location;

e) analysing the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set;

f) determining, based on said correlation, the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set.

2. A method according to claim 1, wherein step f) further comprises:

g) comparing the spatial location of phase distribution features within the measured and modelled phase data sets with respect to the source or receiver position to determine the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set.

3. A method according to claim 2, wherein said phase distribution features comprise a plurality of concentric ring-like features, and step g) further comprises comparing the radial spacing from the source or receiver location of a phase distribution feature within a selected ring-like feature in said modelled phase data set with respect to the radial spacing from the source or receiver location of a corresponding phase distribution feature in a corresponding ring-like feature in said measured phase data set.

4. A method according to claim 3, wherein if, in a radial direction, a phase distribution feature within an n^th ring-like feature in said modelled phase data set is closest to a corresponding phase distribution feature in an m^th (where m≠n) ring-like feature in said measured phase data set, the method further comprises:

h) rejecting the subsurface starting model.

5. A method according to claim 3, wherein if, in a radial direction, a phase distribution feature within an n^th ring-like feature in said modelled phase data set is closest to a corresponding phase distribution feature in an n^th ring-like feature in said measured phase data set, the method further comprises:

i) proceeding further with said subsurface starting model.

6. A method according to claim 5, wherein if, in a radial direction, a phase distribution feature within an n^th ring-like feature is directly adjacent a corresponding phase distribution feature in an n^th ring-like feature in said measured phase data set, step i) further comprises:

j) validating said subsurface starting model.

7. A method according to claim 5, wherein step i) further comprises:

k) obtaining a further modelled phase data set generated from an «^Λ iteration (where n is an integer) of a full waveform inversion process carried out on said subsurface starting model, the further modelled phase data set comprising the phase component of further modelled values taken at the same plurality of discrete data points, at said one or more selected frequencies and for said at least one physical parameter;

1) generating, for said at least one source or receiver position of the further modelled phase data set, a further modelled phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location;

m) analysing the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set and the further modelled phase data set;

n) determining, based on said correlation, the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set.

8. A method according to claim 7, wherein step m) further comprises:

o) analysing the change in phase distribution between the modelled phase distribution data set and the further modelled phase distribution data set to identify convergence or divergence from the measured phase distribution data set.

9. A method according to claim 8, wherein step o) further comprises:

p) monitoring the location, distribution and/or movement of phase distribution rings within said phase data sets.

10. A method according to claim 8 or 9, wherein if, in step o) the phase distribution change from said modelled phase distribution data set to said further modelled phase distribution data set is divergent with respect to said measured phase distribution data set, the method further comprises: q) rejecting said starting model or modifying said starting model.

11. A method according to claim 8 or 9, wherein if, in step o), the change in phase distribution from said modelled phase distribution data set to said further modelled phase distribution data set is convergent with respect to said measured phase distribution data set, the method further comprises: r) accepting said starting model or generating further iterations of said modelled data set.

12. A method according to any one of claims 1 to 11, wherein step a) further comprises:

s) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points;

t) applying a time-domain window to said measured time-domain values of said measured seismic data set to generate a windowed measured seismic data set;

u) performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components;

v) selecting said phase components for said measured phase data set; and

wherein step b) comprises: w) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points;

x) applying a time-domain window to said modelled time-domain values of said modelled seismic data set to generate a windowed modelled seismic data set;

y) performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowed modelled values comprising amplitude and phase components; and

z) selecting said phase components for said modelled phase data set.

13. A method according to claim 12, wherein step t) and x) further comprises applying a shift to the measured or modelled data.

14. A method according to claim 12 or 13, wherein, prior to step e), the method further comprises:

aa) determining whether the spatial consistency of the measured and modelled phase distribution data sets meets a predetermined criterion.

15. A method according to claim 14, wherein if, in step aa) said predetermined criterion is not met, the method further comprises:

bb) selecting a new time-domain window parameter; and

cc) repeating steps s) to z) with said new time-domain window parameter.

16. A method according to any one of the preceding claims, wherein the discrete data points of said measured and modelled seismic data sets are a function of source and receiver locations.

17. A method according to claim 6 or 11 , wherein if it is confirmed in step m) that the optimised subsurface model is sufficiently accurate, the method further comprises the step:

dd) utilising said optimised model for sub-surface exploration.

18. A method according to any one of the preceding claims, wherein said at least one physical parameter is pressure.

19. A method according to any one of the preceding claims, wherein said at least one physical parameter comprises particle velocity or displacement.

20. A method according to any one of the preceding claims, wherein the measured data set and the modelled data set comprise data taken at a plurality of discrete frequencies.

21. A method according to any one of the preceding claims, wherein the measured data set and the modelled data set comprise values of a plurality of physical parameters.

22. A method according to any one of the preceding claims, wherein, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a full waveform inversion iterative process.

23. A computer program product executable by a programmed or programmable processing apparatus, comprising one or more software portions for performing the steps of any one of claims 1 to 22.

24. A computer usable storage medium having a computer program product according to claim 23 stored thereon.

25. A method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising:

a) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points;

b) applying a time-domain window to said measured time-domain values of said measured seismic data set to generate a windowed measured seismic data set;

c) performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components;

d) generating, for at least one source or receiver position, a measured phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location;

e) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points;

f) applying a time-domain window to said modelled time-domain values of said modelled seismic data set to generate a windowed modelled seismic data set;

g) performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowed modelled values comprising amplitude and phase components; and

h) generating, for said at least one source or receiver position, a modelled phase distribution data set comprising the phase distribution for a single source or receiver location as a function of spatial location;

i) determining whether the spatial consistency of the measured and modelled phase data sets meets a predetermined criterion; and, if predetermined criterion is not met, the method further comprises:

j) selecting a new time-domain window parameter; and

k) repeating steps b) to d) and f) to i) with said new time-domain window parameter until said predetermined criterion is met; and

1) determining the validity of said subsurface starting model for interpreting said measured seismic data set.

26. A method according to claim 25, wherein step 1) comprises:

m) receiving at least one further modelled seismic data set generated from an n^th iteration (where n is an integer) of a full waveform inversion process carried out on said modelled seismic data set, said at least one further modelled seismic data set comprising modelled values of at least one physical parameter at the same plurality of discrete data points;

n) determining, for a plurality of discrete data points, a starting check residual variable from the measured values and the modelled values, the check residual variable of a discrete data point comprising the mismatch between a variable of the measured values and said variable of the modelled values of the at least one physical parameter for that discrete data point;

o) calculating a change residual value for each of a plurality of discrete data points, the change residual for a discrete data point comprising the difference between the starting check residual variable and the check residual variable for the n^th iteration; and

p) comparing the change residual of each of a plurality of discrete data points with the check residual variable of the n^th iteration for the same discrete data points to determine the validity of the starting model for interpreting said measured seismic data set.

27. A method according to claim 25, wherein step 1) comprises:

q) calculating, from the phase component of said windowed measured values and the phase component of said windowed modelled values, the phase residual at said one or more selected frequencies, the phase residual of a discrete data point comprising the mismatch between the phase of the windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point;

r) generating a phase residual data set at said one or more selected frequencies and for said at least one physical parameter;

s) identifying cycle-skipped data within said phase residual data set; and

t) evaluating the proportion of cycle-skipped data within said phase residual data set to determine the validity of the starting model for interpreting said measured seismic data set.

28. A method according to claim 25, wherein step 1) comprises:

u) analysing the correlation between the spatial location of one or more phase distribution features within the measured phase distribution data set and the spatial location of one or more corresponding phase distribution features within the modelled phase data set; and

29. A method for subsurface exploration, the method comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement, and comprising the steps of:

a) providing, from said seismic measurement, a measured seismic data set obtained from said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain or frequency-domain values of at least one physical parameter measured at a plurality of discrete data points;

b) generating, using a subsurface model of a portion of the Earth, at least one modelled seismic data set comprising modelled time-domain values or frequency-domain values of at least one physical parameter at the same plurality of discrete data points;

c) if, in step b) the data is provided in the time-domain, performing a Fourier transform on said modelled seismic data set to obtain modelled values taken at one or more selected frequencies, said modelled values comprising amplitude and phase components;

d) if, in step a) the data is provided in the time-domain, performing a Fourier transform on said seismic measured data set to obtain measured values taken at one or more selected frequencies, said measured values comprising amplitude and phase components;

e) calculating, from the phase component of said measured values and the phase component of said modelled values, the unwrapped phase residual at said one or more selected frequencies, the unwrapped phase residual of a discrete data point comprising the mismatch between the phase of the measured values and the modelled values of the at least one physical parameter for that discrete data point;

f) generating a phase residual data set from the measured and modelled dataset;

g) applying a weight to the discrete data points within said phase residual dataset in dependence upon the magnitude of the unwrapped phase residual for one or more discrete data points;

h) utilising an FWI process to minimise the weighted phase residual dataset to produce an updated subsurface model; and

i) utilising said updated model for subsurface exploration.

30. A method according to claim 29, wherein step g) comprises weighting a particular discrete data point with a value between 0 and 1 in dependence upon the magnitude of the unwrapped phase residual for that discrete data point.

31. A method according to claim 30, wherein in step g) a discrete data point is given a value of 0 if said data point has a value of 180° (π) or more.

32. A method according to claim 31, wherein in step g) a discrete data point is given a value of 0 if said data point has a value of 160° (0.88π) or more.

33. A method according to claim 30, 31 or 32, wherein if a data point has a value of 0, said data point is ignored or deleted from the measured and the updated modelled seismic data sets.

34. A method according to claim 30, wherein when, in step g), a discrete data point has a nonzero value between 0 and 1 , the data point is given a weighting to the minimisation in step h) based on said non-zero value, with a value of 1 being given full weight.

35. A method according to claim 30, wherein step g) comprises weighting a particular discrete data point with a value of 0 or 1 in dependence upon the magnitude of the unwrapped phase residual for that discrete data point, data with a value of 0 being deleted from the measured and the updated modelled seismic data sets or ignored in step h).

36. A method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising:

a) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured values of at least one physical parameter measured at a plurality of discrete data points;

b) receiving a starting modelled seismic data set generated from a subsurface starting model, said starting modelled seismic data set comprising starting modelled values of at least one physical parameter at the same plurality of discrete data points;

c) receiving at least one further modelled seismic data set generated from an n^th iteration (where n is an integer) of a full waveform inversion process carried out on said starting modelled seismic data set, said at least one further modelled seismic data set comprising modelled values of at least one physical parameter at the same plurality of discrete data points;

wherein the method further comprises determining, for at least one iteration, the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set by:

d) determining, for a plurality of discrete data points, a starting check residual variable from the measured values and the starting modelled values, the check residual variable of a discrete data point comprising the mismatch between a variable of the measured values and said variable of the modelled values of the at least one physical parameter for that discrete data point;

e) calculating a change residual value for each of a plurality of discrete data points, the change residual for a discrete data point comprising the difference between the starting check residual variable and the check residual variable for the n^th iteration; and

f) comparing the change residual of each of a plurality of discrete data points with the check residual variable of the n^th iteration for the same discrete data points to determine the validity of the starting model for interpreting said measured seismic data set.

37. A method according to claim 36, wherein step d) further comprises:

determining, for a plurality of discrete data points, a check residual variable from the measured values and the modelled values from said n^th iteration;

38. A method according to claim 36, wherein step e) further comprises:

determining, for a plurality of discrete data points, a change residual variable from the starting modelled values and the modelled values from said n^th iteration;

39. A method according to any one of claims 36 to 38, wherein the residual variable is a phase residual and step d) comprises calculating a phase residual for a plurality of discrete data points, the phase residual of a discrete data point comprising the mismatch between the phase of the measured values and the modelled values of the at least one physical parameter for that discrete data point.

40. A method according to any one of claims 36 to 38, wherein the residual variable is an unwrapped phase residual and step d) comprises calculating an unwrapped phase residual for a plurality of discrete data points, the unwrapped phase residual of a discrete data point comprising the mismatch between the unwrapped phase of the measured values and the unwrapped phase of the modelled values of the at least one physical parameter for that discrete data point.

41. A method according to claim 40, wherein said step of calculating an unwrapped phase residual comprises:

g) calculating, for a plurality of groups of two discrete data points, the principal value of the difference in phase residual at said one or more frequencies between the discrete data points in a group such that each group has a principal value of the phase residual difference associated therewith; and

h) determining, from phase residual difference of at least some of said plurality of groups, the unwrapped phase residual for at least some of said plurality of data points.

42. A method according to claim 41, wherein step g) further comprises selecting a subset of groups of discrete data points from the total number of groups based on a selection criterion.

43. A method according to claim 42, wherein said selection criterion comprises selecting the groups having the smallest difference between the phase residuals of their constituent points.

44. A method according to claim 42, wherein said selection criterion comprises selecting groups on the basis that the spatial separation between discrete data points within a group does not exceed a spatial threshold value.

45. A method according to claim 44, wherein said spatial threshold value is twice the minimum spatial separation of the discrete data points.

46. A method according to any one of claims 40 to 44, wherein step g) further comprises applying spatial averaging to reduce the effects of noise.

47. A method according to any one of claims 36 to 46, wherein step c) comprises updating said subsurface model of a portion of the Earth by a localised gradient-based method.

48. A method according to any one of claims 36 to 47, wherein step c) further comprises:

i) constructing an objective function; and

j) minimising said objective function by modifying at least one coefficient of said subsurface model of a portion of the Earth to provide an updated subsurface model of a portion of the Earth.

49. A method according to claim 48 when dependent upon claim 41, wherein said objective function comprises a logarithmic objective function consisting of the unwrapped phase residual summed over a plurality of data points.

50. A method according to any one of claims 36 to 49, wherein step a) comprises performing a Fourier transform on said seismic measurement data to obtain a plurality of frequency components of said data set, and selecting one or more frequencies therefrom.

51. A method according to any one of claims 36 to 50, wherein in step b) and/or c) said modelled seismic data set is generated at a selected frequency directly from the subsurface model of a portion of the Earth.

52. A method according any one of claims 36 to 51 , wherein step b) and/or c) comprises performing a Fourier transform on modelled seismic data generated from said subsurface model to obtain a plurality of frequency components of said data set, and selecting one or more frequencies therefrom.

53. A method according to claim 36, 37 or 38, wherein the check residual variable is a time-offset residual and step d) comprises calculating, for a plurality of discrete data points, the cross-correlation in the time-domain between modelled and measured data at various time offsets between the datasets, with the time-offset residual for a discrete data point comprises being the time offset that maximises the cross-correlation.

54. A method according to any one of the preceding claims, wherein step f) comprises comparing the sign of the change residual with the sign of the residual variable of the n^th iteration for at least some of said plurality of data points to determine the validity of the starting model for interpreting said measured seismic data set.

55. A method according to claim 54, wherein step f) further comprises determining a check ratio for at least some of said plurality of data points, where the check ratio is the ratio of the number of data points for which the sign of the change residual is the same as the sign of the residual variable of the n^th iteration to the number of data points for which the sign of the change residual is different from the sign of the residual variable of the n^th iteration.

56. A method according to claim 54, further comprising the step of deleting the data points from the measured data set for which sign of the change residual is different from the sign of the residual variable of the n^th iteration to generate a modified measured data set.

57. A method according to claim 56, further comprising utilising said modified measured data set in an FWI minimisation process to generate a new starting model.

58. A method according to claim 57, further comprising repeating steps a) to f) utilising the new starting model.

59. A method according to claim 55, wherein the validity of the starting model for interpreting said measured seismic data set is determined at least in part based upon said check ratio.

60. A method according to claim 59, wherein the starting model is considered valid if the check ratio is 50% or greater.

61. A method according to claim 60, wherein the starting model is considered valid if the check ratio is 70% or greater.

62. A method according to any one of claims 36 to 54, wherein step f) comprises calculating the ratio of the change residual to the residual variable of the n^th iteration for at least some of said plurality of data points.

63. A method according to any one of claims 36 to 62, wherein the discrete data points of said measured and modelled seismic data sets are a function of source and receiver locations.

64. A method according to claim 63, wherein step f) further comprises comparing the change residual with the residual variable of the n^th iteration be performed as a function of receiver and/or source locations.

65. A method according to any one of claims 36 to 64, wherein based upon the comparison in step f), the method may comprise one of: switching to a different full waveform inversion method;

repeating the full waveform inversion again starting from the starting model; or modifying the starting model.

66. A method according to claim 65, further comprising the step of

k) confirming the accuracy of an optimised subsurface model generated using said starting model.

67. A method according to any one of claims 36 to 66, further comprising the step of:

1) confirming or rejecting the validity of said sub surface starting model based upon step f).

68. A method according to claim 67, wherein if it is confirmed in step 1) that the optimised

subsurface model is sufficiently accurate, the method further comprises the step:

m) utilising said optimised model for sub-surface exploration.

69. A method according to any one of claims 36 to 68, wherein said at least one physical parameter is pressure.

70. A method according to any one of claims 36 to 69, wherein said at least one physical parameter comprises particle velocity or displacement.

71. A method according to any one of claims 36 to 70, wherein the measured data set and the modelled data set comprise data taken at a plurality of discrete frequencies.

72. A method according to any one of claims 36 to 71 , wherein the measured data set and the modelled data set comprise values of a plurality of physical parameters.

73. A method according to any one of any one of claims 36 to 72, wherein, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a full waveform inversion iterative process.

74. A method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising:

b) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points;

c) applying a time-domain window to said modelled time-domain values of said modelled seismic data set to generate a windowed modelled seismic data set;

d) performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowed modelled values comprising amplitude and phase components;

e) applying a time-domain window to said measured time-domain values of said measured seismic data set to generate a windowed measured seismic data set;

f) performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components;

wherein the method further comprises determining the validity of said subsurface starting model for performing full waveform inversion on said measured seismic data set by:

g) calculating, from the phase component of said windowed measured values and the phase component of said windowed modelled values, the phase residual at said one or more selected frequencies, the phase residual of a discrete data point comprising the mismatch between the phase of the windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point;

h) generating a phase residual data set at said one or more selected frequencies and for said at least one physical parameter; i) identifying outlying data within said phase residual data set.

75. A method according to claim 74, wherein said outlying data comprises cycle-skipped data and the method further comprises j) evaluating the proportion of cycle-skipped data within said phase residual data set to determine the validity of the starting model for interpreting said measured seismic data set.

76. A method according to claim 75, wherein step j) comprises calculating the proportion of cycle-skipped data as a function of the total data in the phase residual data set.

77. A method according to claim 76, wherein the starting model is rejected if the proportion of cycle-skipped data exceeds a pre-determined threshold.

78. A method according to claim 77, wherein the pre-determined threshold is 40%.

79. A method according to claim 78, wherein the pre-determined threshold is 20%.

80. A method according to any one of claims 75 to 79, further comprising:

k) repeating steps b) to d) for a further modelled seismic data set generated from an n^th iteration (where n is an integer) of a full waveform inversion process carried out on said modelled seismic data set to generate further windowed modelled values comprising amplitude and phase components;

1) repeating steps g) to h) utilising said further windowed modelled values to generate a further phase residual data set the phase residual at said one or more selected frequencies;

m) identifying cycle-skipped data within said further phase residual data set; and wherein step j) further comprises:

n) comparing the proportion of cycle skipped data within said further phase residual data set to said proportion of cycle-skipped data within said phase residual set to determine the validity of the starting model for interpreting said measured seismic data set.

81. A method according to claim 80, wherein step n) further comprises determining the change in proportion of cycle-skipped data between said phase residual set and said further phase residual set.

82. A method according to claim 81, wherein if, in step n), the proportion of cycle-skipped data increases from said phase residual set to said further phase residual set, the method further comprises: o) rejecting said starting model.

83. A method according to claim 81, wherein if, in step n), the proportion of cycle-skipped data decreases from said phase residual set to said further phase residual set, the method further comprises:

p) accepting said starting model or generating further iterations of said modelled data set.

84. A method according to any one of claims 75 to 83, further comprising:

q) repeating steps b) to d) for a further modelled seismic data set generated from a different starting model to generate further windowed modelled values comprising amplitude and phase components;

r) repeating steps g) to h) utilising said further windowed modelled values to generate a further phase residual data set the phase residual at said one or more selected frequencies;

s) identifying cycle-skipped data within said further phase residual data set; and wherein step j) further comprises:

t) comparing the proportion of cycle skipped data within said further phase residual data set to said proportion of cycle-skipped data within said phase residual set to determine which of said starting model and said further starting model is most suitable for interpreting said measured seismic data set.

85. A method according to claim 84, wherein step t) further comprises determining the change in proportion of cycle-skipped data between said phase residual set and said further phase residual set.

86. A method according to claim 85, wherein if, in step t), the proportion of cycle-skipped data increases from said phase residual set to said further phase residual set, the method further comprises: u) rejecting said further starting model.

87. A method according to claim 85, wherein if, in step t), the proportion of cycle-skipped data decreases from said phase residual set to said further phase residual set, the method further comprises:

v) rejecting said starting model.

88. A method according to any one of claims 75 to 87, wherein the discrete data points of said measured and modelled seismic data sets are a function of source and receiver locations.

89. A method according to claim 75, further comprising the step of

w) confirming the accuracy of an optimised subsurface model generated using said starting model.

90. A method according to any one claims 75 to 89, further comprising the step of:

x) confirming or rejecting the validity of said sub surface starting model based upon step j).

91. A method according to claim 90, wherein if it is confirmed in step x) that the optimised

y) utilising said optimised model for sub-surface exploration.

92. A method according to any one of claims 75 to 91, wherein said at least one physical parameter is pressure.

93. A method according to any one of of claims 75 to 92, wherein said at least one physical parameter comprises particle velocity or displacement.

94. A method according to any one of claims 75 to 93, wherein the measured data set and the modelled data set comprise data taken at a plurality of discrete frequencies.

95. A method according to any one of claims 75 to 94, wherein the measured data set and the modelled data set comprise values of a plurality of physical parameters.

96. A method according to any one of of claims 75 to 95, wherein, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a full waveform inversion iterative process.

97. A computer program product executable by a programmed or programmable processing apparatus, comprising one or more software portions for performing the steps of any one of claims 74 to 96.

98. A computer usable storage medium having a computer program product according to claim 97 stored thereon.

99. A method according to any one of claims 74 to 96, further comprising the step of, subsequent to step h), unwrapping the phase residual data set to generate an unwrapped phase residual data set.

100. A method according to any one of claims 74 to 96, further comprising the step of deleting the data points associated with outlying data in step i) from the measured seismic data set to generate a modified measured data set.

101. A method according to claim 100, further comprising the step of deleting the data points from the measured data set for which sign of the change residual is different from the sign of the residual variable of the n^th iteration to generate a modified measured data set.

102. A method according to claim 101, further comprising utilising said modified measured data set in an FWI minimisation process to generate a new starting model.

103. A method according to claim 102, further comprising repeating steps a) to f) utilising the new starting model.