WO2013093467A1 - Method of, and apparatus for, full waveform inversion - Google Patents

Abstract

There is provided a method for subsurface exploration. The method comprises generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement. The method comprises: 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 values of at least one physical parameter measured at a plurality of discrete data points; generating, using a subsurface model of a portion of the Earth, at least one modelled seismic data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points; 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; 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; 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; 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; calculating, from the phase component of said windowed measured values and the phase component of said windowed 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 windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point; generating an unwrapped phase residual data set at said one or more selected frequencies and for said at least one physical parameter; minimising the value of the unwrapped phase residual data set to produce an updated subsurface model; and providing an updated model for subsurface exploration.

Description

Method of, and Apparatus for. Full Waveform Inversion

The present invention relates to a method of, and apparatus for, windowed phase-iinwrapped full waveform inversion in order to solve for a best-fit model a portion of the Earth. More particularly, the present invention relates to a method of, and apparatus for full waveform inversion of measured seismic data to generate a subsurface model of at least a portion of the Earth. The model is used to migrate the seismic data to generate an image of the subsurface.

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 is 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 retum 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. Here 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 (Sirgue, L., and R. G. Pratt, 2004, Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies: Geophysics, 69, 231-248).

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) (v² +toW(x))t.(x, s) = -S<5(x- s) where 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) l¾ = Re∑!¾ (_M(r, _S) - </(r_{) S}))^*

Then, from Equation 2):

(V ² + ω ²m ² (x)) ¾g = -2ω ²m (x (x s)S (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 atx', given by equation 6)

Next, the solution is evaluated at the receiver location, given by equation 7): du(r ,s)

dm(x')

V)

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

8) ^"(*') 2co²m ..-(x~ ') ^/u--{x~ s)G( vx" ^> r- )^/ Inserting equation

8) into equation 4) gives equation 9) for single receiver: dE

dtn(x') = 2<D²m(x^r) Re (u(x^r, s)w* (x^r, r, s)^

9) 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:

12) w(x', s) =∑ G^* (χ', r) (u(r, s) - d(r, s))

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 difficulties associated with obtaining correct convergence. Firstly, FWI methodology for the objective function above relies upon the Born 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.

Alternative approaches have been suggested to mitigate the issues discussed above. One such approach is shown and described in "Comparison of waveform inversion, part 1: Conventional Wavefield vs Logarithmic Wavefield' Shin et al, Geophysical prospecting, 2007, 55, 449-464 and WO-A-2009/136708. In this disclosure, a logarithmic wavefield residual, log w-log d is used instead of a conventional wavefield residual, u-d as the measure of mismatch between modelled and measured data. This enables a separation of phase and amplitude in the residual. As a result the inversion can be geared towards phase at early iterations as opposed to amplitude to increase the chance of global minimum convergence. Further the phase component of the logarithmic wavefield residual can be unwrapped (adjusted to account for cycle-skips). This leads to the phase-unwrapped objective function which consists of some norm of unwrapped phi summed over source and receiver positions and selected frequencies. Again the model is updated in the direction of steepest descent.

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 frequencies will summed when forming the final model update.

14) The gradient of this objective function with respect to the model φ' (r,s)

dm (x') <P'(r,s)

To calculate this we note:

dg>'(r,s) _ j_m 1 du (r,s) \

Using equation X gives: I V) ¾^ = 2co ½( ') Im (^G( ',r))

We then arrive at:

₁^ = 2cd ²m(x ') Im (u(x', s)w^* (x ', s)⁾)

18) where:

19) , 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 ii) 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.

An advantage of using unwrapped phase is that it no longer requires the difficult to attain assumption that the modelled wavefield lies within half a cycle of the measured wavefield. This approach has the potential to reduces greatly convergence constraints and widen the applicability of the FWI technique. However, to date, phase unwrapping has rarely been attempted in FWI. This is because unwrapping raw seismic data is problematic and prone to failure even at low end of the measurable frequency range where less of the data is cycle-skipped. In addition, whilst phase unwrapping is technically unchallenging for a single dimension, in multi-dimensional data it is not always possible to unwrap the phase consistently across the entire data set.

An example is illustrated in Figures 7 to 9. Figure 7 a) shows modelled time-domain data and Figure 7b) 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 wavefields u (for the modelled data) and d (for the measured data).

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

As shown in Figure 7c), 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 -π. At this point, the phase residual is cycle-skipped by one cycle. The phase residual increases beyond +π again at around receiver number 400 and wraps around to -π. Therefore, beyond receiver number 400, the phase residual data is cycle-skipped by two cycles.

Figure 8a) shows the data plotted in Figure 7c) on a larger scale where the same data (i.e. phase residual at 3.125 Hz) is plotted on a scale of +/-4 π, with +/- π bounds shown by the dotted line. The cycle-skipped data can clearly be seen.

It is possible to avoid cycle-skipped misconvergence by unwrapping the phase residual in order to obtain the unwrapped phase residual. This may be done by unwrapping the phase of each of the modelled and measured data individually, or by unwrapping their difference (i.e. unwrapping the phase residual).

Figure 8b) shows the phase residual data of Figure 8a) unwrapped. The wrapped cycle- skipped data is also shown for clarity. The unwrapped phase residual data is obtained by adding 2π to the data for each cycle skip as shown. Therefore, 2π is added for a single cycle skip, 4π for two cycle skips etc.

This approach is straightforward when considering a one dimensional phase residual, i.e. Φ(Γ,«) as a function of r_x alone. However, 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.

Figure 9 shows a plot of φ(τ,«) as a function of r_x and s_x The zero-offset point occurs along diagonal of symmetry. The shading represents the phase residual values as shown on the scale to the right of the main figure. As shown see the + π /- π transitions abruptly terminate. Therefore, it is not possible to undo the wrap around in a spatially consistent manner to provide a spatially-consistent unwrapped phase residual. Consequently, in known arrangements, the phase unwrapping method cannot be used to calculate with sufficient accuracy many-dimensional results utilising real-world data.

In summary, estimating the number of cycle skips (to unwrap the phase) between the modelled and recorded data is problematic in known arrangements. This is particularly so in realistic situations when receiver and source locations will be a function of multiple spatial dimensions (x, ,z) as opposed to a co-linear layout such as in Fig 4c. Furthermore, in FWI the data-points are irregularly spaced making phase unwrapping especially difficult and prone to erroneous results.

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. Known methods to address this problem have complications which prevent their widespread use. The present invention seeks to address these issues.

According to a first 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, the method comprising: 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 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 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; g) calculating, from the phase component of said windowed measured values and the phase component of said windowed 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 windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point; h) generating an unwrapped phase residual data set at said one or more selected frequencies and for said at least one physical parameter; i) minimising the value of the unwrapped phase residual data set to produce an updated subsurface model; and j) providing an updated model for subsurface exploration. By providing such a method, cycle skipped mis -convergence is avoided, enabling a truer convergence towards a more accurate model of the subsurface of the portion of the Earth under investigation. Additionally, the present method enables the use of a less accurate starting model or data taken at higher frequencies than prior art methods allow, whilst still achieving excellent accuracy of convergence.

In one embodiment, step g) comprises: k) 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 1) 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 k) 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 k) further comprises applying spatial averaging to reduce the effects of noise. In one embodiment, step i) further comprises: m) constructing an objective function which consists of the sum over some or all the data points of the unwrapped phase residual; n) 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; and o) repeating steps b) to h) and m) to o) for the updated subsurface model of a portion of the Earth until a convergence criterion is met.

In one embodiment, said objective function comprises a logarithmic wavefield objective function consisting of some norm of the unwrapped phase residual summed over a plurality of data points. In one embodiment, the method further comprises the step of: p) utilising a conventional objective function at later iterations of steps m) to o). In one embodiment, the step of applying a time domain window in step c) and step e) comprises multiplying the measured or modelled seismic data set by a predetermined function. In one embodiment, the function is selected from the group of:

Gaussian; cosine, Hann, Hamming, linear; exponential; triangular; quadratic; non-linear; step; and delta.

In one embodiment, the method further comprises, after step c): q) applying a time-domain offset to said time-domain window to shift said time-domain window by a predetermined time period; r) applying said shifted time-domain window to said modelled time-domain values of said modelled seismic data set to generate a further windowed modelled seismic data set; and s) generating, from said windowed modelled data sets, a collated windowed modelled data set, said collated windowed modelled data set comprising data values which are a function of said time-domain window offset; and after step e): s) applying a time-domain offset to said time-domain window to shift said time- domain window by a predetermined time period; t) applying said shifted time-domain window to said measured time-domain values of said measured seismic data set to generate a further windowed measured seismic data set; generating, from said windowed measured data sets, a collated windowed measured data set, said collated windowed measured data set comprising data values which are a function of said time-domain window offset.

In one embodiment, the convergence criterion in step o) comprises: v) 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 nth iteration; and w) comparing the change residual of each of a plurality of discrete data points with the check residual variable of the nth 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 j) further comprises: utilising said updated model for sub-surface exploration. In one embodiment, said at least one physical parameter is pressure.

According to a second 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, the method comprising: 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 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 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) applying a time-domain offset to said time- domain window to shift said time-domain window by a predetermined time period; e) applying said shifted time-domain window to said modelled time-domain values of said modelled seismic data set to generate a further windowed modelled seismic data set; f) generating, from said windowed modelled data sets, a collated windowed modelled data set, said collated windowed modelled data set comprising data values which are a function of said time-domain window offset; g) 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; h) applying a time-domain offset to said time-domain window to shift said time-domain window by a predetermined time period; i) applying said shifted time- domain window to said measured time-domain values of said measured seismic data set to generate a further windowed measured seismic data set;j) generating, from said windowed measured data sets, a collated windowed measured data set, said collated windowed measured data set comprising data values which are a function of said time-domain window offset; k) calculating, from the windowed measured values and windowed modelled values, the residual for a plurality of discrete data points to generate a residual data set, where the residual of a discrete data point comprises the mismatch between the windowed measured values and windowed modelled values of the at least one physical parameter for that discrete data point; 1) minimising the residual data set to produce an updated subsurface model; and m) providing an updated model for subsurface exploration.

In one embodiment, step 1) further comprises: n) constructing an objective function; o) 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; and p) repeating steps b) to k), n) and o) for the updated subsurface model of a portion of the Earth until a convergence criterion is met.

In one embodiment, the step of applying a time domain window in one or more of steps c), e), g) and i) comprises multiplying the measured and/or modelled seismic data set by a predetermined function. In one embodiment, the function is selected from the group of: Gaussian; cosine, Hann, Hamming, linear; exponential; triangular; quadratic; non-linear; step; and delta. In one embodiment, steps d), e), i) and j) are performed a plurality of times to generate multiple windowed measured data sets and multiple windowed modelled data sets. In one embodiment, step o) further comprises:

utilising said updated model for sub-surface exploration.

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.

According to a third 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 the first to third aspects.

According to a fourth aspect of the present invention, there is provided a computer usable storage medium having a computer program product according to the second 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) illustrates a modelled seismic data set;

Figure 7 b) illustrates a measured seismic data set;

Figure 7 c) illustrates the phase residual between the data of Figures 7 a) and 7 b); Figure 8 a) illustrates the phase residual of Figure 7 c) on a larger scale; Figure 8 b) illustrates phase unwrapping of the data of Figure 8 a);

Figure 9 illustrates pressure as a function receiver and source location showing terminating transitions in the prior art;

Figure 10 is a flow chart illustrating the method according to an embodiment of the present invention;

Figures 1 la) - d) are schematic illustrations of windowing modelled seismic trace data;

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

Figure 13 illustrates pressure as a function receiver and source location similar to Figure 9 showing terminating transitions in the prior art;

Figure 14 illustrates calculation of non-zero curl to determine terminating transition points;

Figure 15 illustrates windowed phase residual data;

Figure 16 is similar to figure 15 but shows data points and groups comprised thereof;

Figure 17 illustrates windowed and unwrapped data;

Figure 18a) and b) illustrate detection of convergence of data;

Figure 19 is a flow chart illustrating the method according to another embodiment of the present invention;

Figure 20a) to c) illustrates the sliding window approach of the third embodiment; and Figure 21a) to c) illustrates the sliding window approach of the second embodiment. 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 obtaining an accurate and high resolution V_p model of the subsurface of a portion of the Earth which generates modelled data that matches the measured seismic data.

A method according to the present invention will now be described with reference to Figure 10. Figure 10 shows a flow diagram of an embodiment of the present invention. 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 dimensions (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 measured. The measured data set is defined as d (r,s) and may be taken at a single frequency or may comprise data taken at plurality of frequencies. However, it is possible to measure other parameters using appropriate technology; for example, to measure the velocity or displacements of particles in three spatial dimensions in addition to the pressure. The present 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. The method now proceeds to step 102.

Step 102: Provide starting model

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

The generated model consists of values of the coefficient 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. Whilst, as illustrated above, it is normally necessary to provide a starting model of sufficient accuracy to ensure non-cycle skipped convergence, in the case of the present invention a less accurate initial model could be used, saving on time and resources whilst enabling accurate convergence. Also not requiring such an accurate starting model would extend the applicability to a wider range of possible applications. The method then proceeds to step 104.

Step 104: Generate modelled data set

In order to model accurately the subsurface region under investigation, it is necessary to generate 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 is generated for the same measurement parameter(s) at the same frequency or frequencies.

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.

From the above analysis, modelled seismic data can be generated for one or more physical parameters and at one or more selected frequencies. This forms the modelled seismic data set u(r,s). The method now proceeds to step 106.

Step 106: Time Window Modelled Seismic 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.

Figure 11 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 11a) 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 7a). As shown in Figure 1 la), 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 eliminate or reduce the problematic points of non-zero curl in g as described previously.

Step 106 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 l ib). 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 11c) shows the function 152 superimposed on the original modelled data 160 prior to multiplication. The effect of the function 162 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 l id) 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 108.

Step 108: Transform modelled data set

The windowed modelled seismic trace data generated in step 106 may, as set out in step 106, 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.

However, as discussed in step 104, a modelled data set at a predetermined or preselected single frequency of interest can be generated directly. Using equation 15) generates single-frequency data. Therefore, the model of the subsurface of a portion of the Earth can be used to readily generate data sets at one or more frequencies. Consequently, in this case, step 112 is not needed.

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

Step 110: Time Window measured seismic data set

The data obtained in step 100 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 110. This is generally similar to step 106 above but is applied to the measured data set.

Figure 12 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 7b). As shown in Figure 12a), 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 measured seismic data set prior to calculating the phase residual is operable to eliminate or reduce the points of non-zero curl in g as described previously.

Step 110 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 12b). 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 106, 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.

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

Figure 12c) 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 12d) 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 112.

Step 112: 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.

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. The method then proceeds to step 114.

Step 114: Calculate wrapped phase residual

As set out above, the windowed measured seismic data is expressed as a function of receiver (r) and source (s) pair position. This is the function d(r,s). The windowed 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 u(r,s). The windowed 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 114, the phase residual between the windowed data sets 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 20):

which is the wrapped difference of the individual phases of u(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 Φ(Γ,«).

By utilising the present method where the original data is windowed in the time domain before being Fourier transformed and the phase residual determined, the abrupt +2π/-2π terminations described previously and associated with known arrangements can be avoided.

Figure 13 shows a graph similar to that of Figure 9. Figure 13 shows a plot of phase residual as a function of source location and receiver location for a single source/receiver dimension (e.g. s_x, r_x) for a conventional FWI approach which does not use the windowing steps 106 and 110.

In Figure 13, the abrupt termination points (for which the phase cannot be unambiguously unwrapped) are shown as circles (open circles for -2π points, closed circles for +2π points). These points inhibit reliable and accurate unwrapping of phase and lead to inaccuracies in the eventual data model.

This condition can be defined precisely by identifying the 'terminations' utilising the differential phase g (i.e. the difference in phase between adjacent or closely related data points) as shown in equation 21 ) below: 21 ) g= phase w(r₂,s₂)i/*(r₂,s₂)w *(r₁, si)d(r_u s ) which enables calculation of the phase residual difference between data points r , «i) and (r₂,

¾).

The abrupt terminations are identified by a non-zero curl of g. This is illustrated in Figure 14. Figure 14 shows a central point X (which may be an abrupt transition) surrounded by other data points. Data points (r_\, s_\) and (r₂, s₂) are illustrated as examples which have a phase residual difference of +80 between them. However, to calculate the curl, it is necessary to calculate the phase residual difference between all of the points in a closed loop as shown in Figure 14.

In other words, when g is summed around a closed loop the sum is non-zero. This means when g is summed along a contour (in r_x,r_y,s_x,S_y and / space) enclosing one of the abrupt transition points it returns a non-zero result. So when points of non zero curl exist in g, the number of cycle- skips at each point cannot be uniquely specified. Therefore, it is not possible to unwrap accurately the phase residual at these points.

However, when the measured and modelled seismic data sets are windowed in accordance with steps 106 and 110 above, these abrupt transitions are reduced or eliminated. Figure 15 shows a graph similar to that of Figure 13. Figure 15 shows a plot of phase residual as a function of source location and receiver location for a single source/receiver dimension (e.g. s_x, r_x) for the present invention utilising the windowing steps 106 and 1 10.

The effect of applying a time-domain window in steps 106 and 1 10 prior to calculation of the phase residual difference eliminates or reduces the points of non-zero curl in g. As shown, the cycle skip boundaries are generally connected and absent abrupt terminations. This makes phase unwrapping possible and at higher frequencies with noisier data. Taking the example of Figure 14, the +2π/-2π transitions now uniquely defines the cycle-skipped portion of the data, enabling unwrapping to be carried out in step 1 16. The method now proceeds to step 1 16.

Step 116: Unwrap phase residual

At step 1 16, the phase residual Φ(Γ,«) calculated in step 114 is unwrapped. For a one dimensional system (such as that illustrated in Figures 8a) and 8b)) it is relatively straightforward to add 2π for each cycle skip where required to unwrap the phase. However, in practice, seismic data comprises at least seven dimensions (three dimensions of space for each of source and receiver, and one of frequency). Therefore, an alternative procedure for unwrapping phase is required.

In this regard, the phase residual difference defined in equation 18) can be utilised to unwrap the phase. As described, 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 show where +2π/-2π transition contours occur and to which appropriate multiples of 2π can be added to unwrap the 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 and can be used to unwrap the phase of the data.

In practice, when working in multi-dimensions and with noisy data where a reduced number of points of non zero curl may still exist after windowing, a robust approach is required. The phase is thus unwrapped by solving equation 22):

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

23) φ (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 the data points where φ (r,s) 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 22) 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 optional step, it is possible to apply spatial averaging to g (low-pass filtering) to reduce the effects of noise prior to subsequent steps. An appropriate filter may be:

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 16 which shows data points 170 superimposed over the phase residual plot of Figure 15. 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 16. 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 is 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 mimimum 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 16 then becomes the exemplary plot shown in Figure 17. As shown in Figure 17, the phase unwrapping removes the cycle skip boundaries in the data plot. Therefore, the objective function can subsequently be minimised without risk of cycle skipped misconvergence occurring. Consequently, the unwrapped phase data set is more likely to be able to converge to an accurate model than a similar wrapped phase data set. The method now proceeds to step 118.

Step 118: 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 24) below:

24) a

Step 120: Minimise objective function and update model The model update derives from the gradient of equation 24) i.e. the partial derivative with respect to a point perturbation of the model m(x) at each position*'. 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 frequencies will summed when forming the final model update.

The gradient of this o ect to the model is:

26) r

To calculate this we note:

Using equation X gives:

2_& ²m(x') lm (^G(x',r))

28)

We then arrive at:

2co ²m(x') lm ( (x'₎ s)w^* (x' , s

29)

where:

As with computational structure of 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 ii) a back-propagated wavefield w which is emitted by a (multi-point) source located the 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 122: Convergence criteria met?

At step 122 it is determined whether convergence criteria have been met. For example, when the method is deemed to have reached convergence when the difference between the data sets reaches a threshold percentage. Alternatively and optionally, the convergence criteria could be analysed more completely as follows.

This involves highlighting modelled data shifting in the wrong direction i.e. away from the recorded waveform, at any given iteration within the procedure. In order to determine how the data is changing with each iteration, the change in the unwrapped phase residual between iterations is determined and analysed. The change residual Δφ is thus defined in equation 31 ) for the case where the starting modelled data set udr,s) and the first iteration modelled data set ui(r,s) are used:

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

32) Δφ_η (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.

Misconvergence 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_n(r,s) and u_n+i (r,s) and so it is not explicitly necessary to calculate or obtain the phase residual for iteration n+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 33): 33) Ag_n (r,s)= phase w„(r₂,s₂)w„*(r₂,s₂)w,_!+; *(r₁, si)u„+i(n, si)

The misconvergence or otherwise of the data of an iteration can be gleaned through a comparison of Δφ_η with ^_nas calculated. Figure 18a) shows a graphical example of φ_η plotted as a function of receiver and source position. In Figure 18a), φ_η is taken at a frequency of 3.125Hz plotted as a function of s_x and r_x. Figure 18b) shows a corresponding plot of Δφ_η plotted as a function of the same receiver and source positions after 1 iteration of phase unwrapped FWI. As shown by a comparison between these figures, in the circled area it can be seen that both φ_nand Δφ_η are greater than zero, indicating that the data is converging in the correct direction.

Therefore, from a general, qualitative, analysis of the data, it can be stated that the starting model for the iterations shown is adequate and converge to an accurate minimum has occurred.

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 and where they agree, the modelled data is shifting in the correct direction. Alternatively, g_n and Ag_n could be compared.

A more complex analysis of the change residual can provide further detailed information to vindicate or reject a particular iteration or final result. This is an advantage of the present method which enables the divergence or misconvergence 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. 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 iteration for interpreting said measured seismic data set is determined at least in part based upon the check ratio. As examples, an iteration may be considered valid if the check ratio is 50% or greater. Alternatively, the iteration may be considered valid if the check ratio is 70% or greater. The analysis can then be completed and a final score or rating given to the iteration or the final result. 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 if required. Also parameters chosen for the inversion may be re-tuned to improve the result of the final result.

If the criteria as set out above have been met, then the method proceeds to step 124 and finishes with the resultant Earth model generated. If the criteria have not been met, then the method proceeds back to repeat steps 104 to 120 as discussed above.

Step 124: Finish

When, at step 124, 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, prior art approaches are unable to unwrap the phase residual in real- world situations where the cycle skips are not uniquely defined. Therefore, the present invention provides a reliable method for unwrapping the phase and providing more accurate convergence. Thus, this minimisation approach is operable to avert cycle-skipped local minimum convergence giving a more accurate final model of a subsurface region of the Earth.

An alternative embodiment of the invention is illustrated in Figure 19. In the alternative embodiment of Figure 19, method steps 208 and 214 are also intended as alternative or additional approaches to the method steps of the earlier embodiment and can be utilised with this embodiment. Step 200: Obtain measured seismic data set

Step 200 corresponds substantially to method step 100 of the previous embodiment.

Therefore, this will not be described again here. Only the steps which are new to this embodiment of the method of the present invention will be described. The method now proceeds to step 202.

Step 202: Provide starting model

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

The generated model consists of values of the coefficient 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. Whilst, as illustrated above, it is normally necessary to provide a starting model of sufficient accuracy to ensure non-cycle skipped convergence, in the case of the present invention a less accurate initial model could be used, saving on time and resources whilst enabling accurate convergence . Also not requiring such an accurate starting model would extend the applicability to a wider range of possible applications. The method then proceeds to step 204.

Step 204: Generate modelled data set

In order to model accurately the subsurface region under investigation, it is necessary to generate 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 is generated for the same measurement parameter(s) at the same frequency or frequencies. 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. From the above analysis, modelled seismic data can be generated for one or more physical parameters and at one or more selected frequencies. This forms the modelled seismic data set u(r,s). The method now proceeds to step 206.

Step 206: Time Window Modelled Seismic Data Set

The data obtained in step 206 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. A time window is then applied to the modelled data set in step 206. This is generally similar to steps 106

Figure 20 a) shows a graph illustrating an aspect of the modelled data 350. 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 7a).

The inventor has found that the effect of applying a time-domain window to the measured seismic data set prior to calculating the phase residual is operable to eliminate or reduce the points of non-zero curl in g as described previously.

Step 206 is operable to apply a time-domain window to the modelled seismic trace data 250. This is done by applying a function 352 to the modelled data. A suitable function may be a Gaussian function 252 as shown in Figure 20b). The Gaussian function 254 has various parameters which can be specified by a user such as the time (τ) at which the function or window 252 is centred and the Full Width Half Maximum (FWHM) which specifies the spread of the function, τ may be determined based upon the earliest arrival times as determined by ray tracing.

However, whilst a Gaussian operator has been specified as the function 252 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. The resulting multiplication corresponds substantially to the data illustrated in Figures 11c) and l id). 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 254 for all data points. The method then proceeds to step 208. Step 208: Repeat Time Windowing of modelled seismic data set

At step 208, the time windowing procedure is repeated with the same general function 252 and same data set 250. However, this time, τ is varied and a new data set generated. Figure 20c) shows the function 252 with τ shifted towards longer times. A new data set is generated therefrom using the modelled data set. This step is repeated as many times as is necessary to generate a new variable in the modelled data set in which the data is dependent upon τ. The method then proceeds to step 210.

Step 210: Transform modelled data set

The windowed modelled seismic trace data generated in step 208 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. However, 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 time-domain 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. However, as discussed in step 306, a modelled data set at a predetermined or preselected single frequency of interest can be generated directly. Using equation 15) generates single-frequency data. Therefore, the model of the subsurface of a portion of the Earth can be used to readily generate data sets at one or more frequencies. Consequently, in this case, step 210 is not needed.

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

Step 212: 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. Figure 21 a) shows a graph illustrating an aspect of the raw trace data 260. 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 7b). As shown in Figure 20a), there is a peak at early times (early arrival pressure waves) with a complicated series of oscillations extending out to longer times. There is also considerable noise in the data 350 due to experimental conditions and other factors in the measuring environment.

The inventor has found that the effect of applying a time-domain window to the measured seismic data set prior to calculating the phase residual is operable to eliminate or reduce the points of non-zero curl in g as described previously.

Step 212 is operable to apply a time-domain window to the measured seismic trace data 260. This is done by applying a function 262 to the measured raw data. A suitable function may be a Gaussian function 262 as shown in Figure 21b). The Gaussian function 354 has various parameters which can be specified by a user such as the time (τ) at which the function or window 262 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 262 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. The resulting multiplication corresponds substantially to the data illustrated in Figures 12c) and 12d)

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 264 for the whole trace. The method then proceeds to step 214.

Step 214: Repeat Time Windowing of measured seismic data set

At step 214, the time windowing procedure is repeated with the same general function 262 and same data set 260. However, this time, τ is varied and a new data set generated. Figure 21c) shows the function 262 with τ shifted towards longer times. A new data set is generated therefrom.

This step is repeated as many times as is necessary to generate a new variable in the data set in which the data is dependent upon τ. The method proceeds to step 216.

Step 216: Transform measured seismic data set

The windowed measured seismic trace data 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.

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. The method then proceeds to step 218.

Step 218: Calculate phase residual

As set out above, the windowed measured seismic data is expressed as a function of receiver (r) and source (s) pair position. This is the function i (r,sx). The windowed modelled seismic data is also expressed as a function of receiver (r) and source (s) pair positionx. This, as set out above, is the function u(r,sr). The windowed 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 218, the phase residual between the windowed data sets 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 34):

which is the wrapped difference of the individual phases of u(r,s,r) and d r,s,r). 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 φ(/·,«,τ). Optionally, it may be desired to unwrap the phase at this point. If so, this is done in accordance with step 116 above. The method now proceeds to step 220.

Step 220: 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 35) below:

35) Ε = ±∑φ

The summation is over source, receiver, τ and selected frequencies

Step 222: Minimise objective function and update model

The model update derives from the gradient of equation 35) 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 frequencies and separate values of r will summed when forming the final model update.

E = ∑<p'² (r, _S)

36)

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

To calculate this we note:

Using equation X gives: ₃₉₎ ¾¾ = 2co ½( ') Im (¾#G( ',,))

We then arrive at: dE

x') 2ω ²m(x') lm (u(x', s)w^* (χ'

40) ^dm(

where:

As with computational structure of 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 ii) a back-propagated wavefield w which is emitted by a (multi-point) source located the receiver positions.

Step 224: Convergence criteria met?

At step 224 it is determined whether convergence criteria have been met. For example, when the method is deemed to have reached convergence when the difference between the data sets reaches a threshold percentage.

Alternatively, the convergence may be checked using the check residual described with reference to the first embodiment.

If the criteria have been met, then the method proceeds to step 226 and finishes with the resultant Earth model generated. If the criteria have not been met, then the method proceeds back to repeat steps 204 to 222 as discussed above. Step 226: Finish

When, at step 226, 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 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. 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, 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. Additionally, whilst the above embodiments have been described with reference to unwrapping the phase residual, this need not be so. Instead, it is possible to unwrap u and d independently prior to calculating the phase residual. The present invention is intended to cover such variations. 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 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 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 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;

g) calculating, from the phase component of said windowed measured values and the phase component of said windowed 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 windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point;

h) generating an unwrapped phase residual data set at said one or more selected frequencies and for said at least one physical parameter;

i) minimising the value of the unwrapped phase residual data set to produce an updated subsurface model; and

j) providing an updated model for subsurface exploration.

2. A method according to any one of the preceding claims, wherein step g) comprises: k) 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

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

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

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

5. A method according to claim 3, 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.

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

7. A method according to any one of claims 2 to 6, wherein step k) further comprises applying spatial averaging to reduce the effects of noise.

8. A method according to any one of the preceding claims, wherein step i) further comprises:

m) constructing an objective function which consists of the sum over some or all the data points of the unwrapped phase residual;

n) 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; and o) repeating steps b) to h) and m) to o) for the updated subsurface model of a portion of the Earth until a convergence criterion is met.

9. A method according to claim 8, wherein said objective function comprises a logarithmic wavefield objective function consisting of some norm of the unwrapped phase residual summed over a plurality of data points.

10. A method according to claim 9, further comprising the step of:

p) utilising a conventional objective function at later iterations of steps m) to o) 11. A method according to any one of the preceding claims, wherein the step of applying a time domain window in step c) and step e) comprises multiplying the measured or modelled seismic data set by a predetermined function.

12. A method according to claim 11, wherein the function is selected from the group of: Gaussian; cosine, Hann, Hamming, linear; exponential; triangular; quadratic; non-linear; step; and delta.

13. A method according to any one of the preceding claims, further comprising, after step c):

q) applying a time-domain offset to said time-domain window to shift said time-domain window by a predetermined time period;

r) applying said shifted time-domain window to said modelled time-domain values of said modelled seismic data set to generate a further windowed modelled seismic data set; and

s) generating, from said windowed modelled data sets, a collated windowed modelled data set, said collated windowed modelled data set comprising data values which are a function of said time- domain window offset; and

after step e):

s) applying a time-domain offset to said time-domain window to shift said time-domain window by a predetermined time period;

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

u) generating, from said windowed measured data sets, a collated windowed measured data set, said collated windowed measured data set comprising data values which are a function of said time- domain window offset.

14. A method according to claim 8, 9 or 10, wherein the convergence criterion in step o) comprises: v) 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 nth iteration; and w) 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.

15. A method according to any one of the preceding claims, wherein step j) further comprises: utilising said updated model for sub-surface exploration.

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

17. 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, the method 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 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 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) applying a time-domain offset to said time-domain window to shift said time-domain window by a predetermined time period;

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

f) generating, from said windowed modelled data sets, a collated windowed modelled data set, said collated windowed modelled data set comprising data values which are a function of said time- domain window offset;

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

h) applying a time-domain offset to said time-domain window to shift said time-domain window by a predetermined time period;

i) applying said shifted time-domain window to said measured time-domain values of said measured seismic data set to generate a further windowed measured seismic data set; j) generating, from said windowed measured data sets, a collated windowed measured data set, said collated windowed measured data set comprising data values which are a function of said time- domain window offset;

k) calculating, from the windowed measured values and windowed modelled values, the residual for a plurality of discrete data points to generate a residual data set, where the residual of a discrete data point comprises the mismatch between the windowed measured values and windowed modelled values of the at least one physical parameter for that discrete data point;

1) minimising the residual data set to produce an updated subsurface model; and

m) providing an updated model for subsurface exploration.

18. A method according to claim 17, wherein step 1) further comprises:

n) constructing an objective function;

o) 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; and p) repeating steps b) to k), n) and o) for the updated subsurface model of a portion of the Earth until a convergence criterion is met.

19. A method according to claim 17 or 18, wherein the step of applying a time domain window in one or more of steps c), e), g) and i) comprises multiplying the measured and/or modelled seismic data set by a predetermined function.

20. A method according to claim 19, wherein the function is selected from the group of: Gaussian; cosine, Hann, Hamming, linear; exponential; triangular; quadratic; non-linear; step; and delta.

21. A method according to any one of claims 17 to 20, wherein steps d), e), h) and i) are performed a plurality of times to generate multiple windowed measured data sets and multiple windowed modelled data sets.

22. A method according to any one of claims 17 to 21 , wherein step o) further comprises: utilising said updated model for sub-surface exploration.

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

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

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

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

27. A computer usable storage medium having a computer program product according to claim 26 stored thereon.