GB2503640A - Quality Assurance in a Full Waveform Inversion Process - Google Patents

Abstract

A method for verifying the accuracy of a model of the subsurface of a portion of the earth generated by a full waveform inversion method. The method comprises receiving a measured seismic data set 100 obtained from a seismic measurement and starting modelled seismic data set 102 generated from a subsurface starting model. A time-domain window 106, 110 is then applied to both the measured and modelled seismic data sets prior to applying a Fourier transform 108, 112. The method further comprises determining the validity of said subsurface starting model for performing FWI on said measured seismic data set by calculating, from the phase component of said windowed measured values and the phase component of said windowed modelled values, the phase residual 114 at said one or more selected frequencies, identifying cycle-skipped data 116 within said phase residual data set; and evaluating the proportion of cycle-skipped data within said phase residual data set to determine the validity of the starting model for interpreting said measured seismic data set.

Description

Method of. and Apparatus for. Oualitv Assurance in a Full Waveform Inversion Process The present invention relates to a method of, and apparatus for, quality control of the full waveform inversion method for calculating a model of a portion of the Earth.

More particularly, the present invention relates to a method of, and apparatus for checking the validity of the full waveform inversion procedurn applied to mea ured seismic data when generating a subsurface model of at least a portion of the Earth.

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

A general schematic illustration of an experimental set up 10 for an undersea seismic survey is shown in Figure I. 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 capablc of crcating sufficicnt vibrational disturbanec. In most scismic survcy cxpcrirncnts, a singlc sourcc is uscd which is shot from multiplc locations.

A plurality of dctcctors 16 is providcd. Thc dctcctors 16 may cornprisc any suitablc vibrational dctcction apparatus. In practicc, dcviccs known as gcophoncs arc 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 scts of lines or in a grid (for 3D acquisition). The dctcctors 16 arc 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 usc, thc acoustic wavc 20 gcncratcd by thc sourcc 12 travcls downwards into the subsurface 22 of the Earth. The subsurface 22, in general, comprises one or more layers or strata 24, 26, 28 formed from rock or other materials. The acoustic wave 20 interacts with the medium, for example by being refracted through the layers and/or reflected off the interfaces between them, and a plurality of return signals 30 is detected by the detectors 16.

In general, the returning signals 30 comprise pressure or elastic waves having different polarisations. P-waves (also known as primary or pressure waves) are longitudinally polarised and comprises alternating rarefactions and compressions in the medium in which the wave is travelling. In other words, in an isotropic environment, the oscillations of a P-wave are parallel to the direction of propagation of the wave.

However, the polarisation of P-waves can vary depending upon the angle at which they strike boundaries. P-waves have the highest velocity and so are the first to be recorded. P-waves travel at a velocity VP in a particular medium. V is a scalar quantity and is, effectively, the speed of sound in a medium. Tt is this quantity V 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 j 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 survcy is composcd of a largc numbcr 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 thrne 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 ofthe acoustic wave 20 occurs at a boundary or interface between two dissimilar materials, or when the elastic properties of a material changes.

Traces are usually obtained at discrete intervals of the order of milliseconds. Each detected return signal 30 forming a seismic trace has a travel time from the source 12 which, for reflected waves, is a two-way travel time from the source 12 to the reflecting element (usually subsurface rock interface) and back up to the detectors 16.

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

Features of interest in prospecting include: faults, folds, anticlines, unconformities, salt domes, reefs.

Key to this process of modelling and imaging a portion of earth is the seismic velocity Vi,. Tn a portion of the volume of the Earth, V maybe estimated in various ways. One method is travel-time tomography. Other methods are velocity analysis from NMO or migration. In tomography, travel times of events are picked from the recorded data and thc vclocity that bcst prcdicts thcsc travcl timcs is found. Prcdictcd travcl tirncs arc calculatcd using ray-bascd approximations.

This approach is robust in thc scnsc that prcdictcd travcl-timcs do not vary significantly whcn thc modcl is varicd and produccs a smooth vclocity rnodcl on thc scaic of the Fresnel zone.

Howover, an alternative techniquc is fill 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 casc, an attcmpt is made to match the frill rccordcd wavcficld with predictcd data.

In other words, the technique involves generating a two or three dimensional model to represent the measured portion of the Earth and aftempts 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 stalling model. A sufficiently accurate starting model for FWI maybe provided by travel-time tomography.

FWI can extract many physical properties (V and V velocities, attenuation, density, anisotropy) of the modelled portion ofthe Earth. However, Vi,, the P-wave velocity, is a particularly important parameter which the subsequent construction ofthe other parameters depends heavily upon. The inversion of V 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 apparcnt to the skilled person.

The FWI technique seeks to obtain an accurate and high resolution V model of the subsurface which generates modelled data that matches the recorded data. Modelled data is calculated using the frill 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 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 thnction 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 ofthe 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 degrcc of accuracy or until sufficient convcrgcncc is obtained. Examples of this tcchniquc arc illustrated in "An overview offuli-waveform in version in exploration geophysics, S. 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 measurc of the mismatch between the recorded data and the modcllcd 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 ofthe 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. U. 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 noisc and the starting model is limited in accuracy.

Localised gradient-based methods arc used which iteratively update an existing model in a direction that derives from the objective function's direction of steepest dcsccnt. In its standard form it is thc L-2 norm of thc data rcsiduals that is minimiscd.

To dcscribc this proccdurc it is sufflcicnt to considcr a singic shot and singic frcqucncy in the objective function. The shot is located at points and w is the angular value of the frcqucncy bcing invcrtcd. Ultimatcly gradicnts from scparatc shots and frcqucncics will summcd ovcr whcn forming thc final modcl updatc.

1) E =+Lk(r,s)_dfr,s)2 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, usuaHy by flnite-diffcrcnce methods, by soMng the frequency domain (uniform-density) acoustic wave equation 2) (V2+oYm2(x))u(xs)=_S6(xs) here the subsurface model nix) is slowness, the inverse ofV (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 fmm the data. The solution is linear in S: 21) 3) uç,s)SGx,s).

where 0 is the green's function which we calculate by solving equation (2) with 8=1.

The gradient direction of objective function (1) is evaluated at the existing model.

This is the partial derivative of E for a point rnodcl perturbation at each subsurface position x'.

4) = Re (u(r, s) -d(r, s)) Thcn, from Equation 2): (v2 + 2m2 (x)) = -2w 2m(x')u(x',s)ö (x -x') This is the same equation as (2) but with a source term now located at model perturbation positionx'. The solution of equation 5)is a wavefield emitted by a source positioned atx', given by equation 6) au(x s) 6) 3m(x') = 2oym(x)u(x,s)G(x,x !0 Next, thc solution is evaluated at the receiver location, given by equation 7): 7) = 2w2Gfr,x')rn(x')u(xçs) Then, applying source-receiver reciprocity the wavefield is now emitted from the receiver location, as specified in equation 8): 8) = 2w2in(x')u(x', s)G(x r) Inserting equation 8) into equation 4) gives equation 9) for single receiver: 9) am(x') = 2w2m(x5Re(u(x,s)W(xP,r,s)) where: 10) w(x',r,s) = (1 (x',r) (u(r,s) -d(r,s)) The gradient of equation 9) is the product of two waveflelds: i) incident wavefield u emitted by a source at original source location s; and ii) a back-propagated wavefield w \vhich 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 11): 11) am') =2w2m(x9Re(u(xn,s)w*(xr,s)) whcrc: I 2) w(fs)= LG*(xP,r)(u(r,s)_dfr,s)) The above can still be calculated as a single wavefleld by taking a multi-point source. The sources are the wavefield residuals at each individual receiver location.

Finally, several frequencies and shots can be inverted simultaneously by summing the gradient contributions. As set out above, the gradient methods can be enhanced using an approximate form for the Hessian and conjugate directions for the model update.

Further a step-length is then calculated to scale the search direction and give the final model update.

Ideally, the above method will lead to a convergence to a model which is a correct representation of a portion of the subsurface of the Earth. However, there are some critical difficulties associated with obtaining correct convergence.

Firstly, FWT methodology for the objective function above relics 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(r4. 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 dfr,s) and ufr,s) is determined. The principal value (between -it and + it) of the phase residual is denoted as 4 and expressed in equation 13): 13) cji(r,s) phase u(r,s)d *(r,s) which is the wrapped difference of the individual phases of ufr4 and clfr,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 -it to + it. In other words, if the phase difference is it and increases, the value of the phase difference will become -it and move back towards it.

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 (asset 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 ofet receiver position towards the left hand side of Figure 4 c), it can be seen that the phase increases towards it before moving to -it (due to the wrapped phase) then increases towards it before moving to -it again (due to the wrapped phase). The apparent discontinuities, or +it transitions, can lead to convergence problems in numerical iterations.

At short offsets (r) is not cycle-skipped. However due to the +it transitions, at longer offsets cF (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 offiets 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 ij3 is moving upwards to convergence, whereas in fact it should be moving downwards through-n 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 FM/Il 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 modcl. Thcrcthrc, performing FWI starting with the starting model of Figurc 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 subsurfitce structures of real-world portions of the Earth. Figures 6 d) to 6 t) illustrate the operation ofFWl 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 frill recovered model after subsequent higher frequency iterations. As shown, there are clear inaccuracies in the recovered model when compared to the original model of Figure 6 a). Features are not resolved clearly, and are inaccurately located. This is a result of cycle skipped mis-convergence and, in practice, would give a sub-optimal final result which could lead to inaccurate identification of natural resources or other sites of interest.

A flirther example of this is given in Figures 7a) to 7fl. Figures 7a) and 7d) show different starting models -Figure 7a) shows a basic starting model and Figure 7d) shows a more accurate starting model.

Figures 7b) and 7e) show the respective models after a single iteration of an FWI process. The same measured data set is used in each case. Further, Figures 7c) and 7f) show the respective models after 15 iterations. As shown, there is a significant difference between the models after only 15 iterations. It can be seen that the less accurate starting model is converging to an inaccurate final model and the inaccuracies will only be amplified with each fbrther iteration.

The inaccuracy of final results is of concern to prospectois and data inteipreters alike and may potentially result in costly and inefficient misidentification of natural resources.

However, known FVTI techniques do not provide a robust indication that misconvcrgcncc has occurred, ict alonc providc tools as how to managc, mcasurc and corrcct deficiencies. Therefore, interpretation of the accuracy of the convergence is difficult using known methods. As a result, it is difficult to determine precisely the accuracy of a final model or to determine which final model (where multiple starting modcls are considered) represents thc most accurate depiction of the subsur&ce portion of the Earth being modelled.

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

According to a first aspect of the present invention, there is provided a method for checking the validity of a geophysical representation of a portion of the volume of the Earth generated from a seismic measurement, the method comprising: a) receiving a measured seismic data set obtained from said seismic measurement of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points; b) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points; c) applying a time-domain window to said modelled time-domain values of said modelled seismic data set to generate a windowed modelled seismic data set; d) performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowcd modcllcd values comprising amplitude and phasc components; c) 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; 0 performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components; wherein the method further comprises determining the validity of said subsurface starting model for performing ifill waveforn inversion on said measured seismic data set by: g) calculating, from the phase component of said windowed measured values and the phase component of said windowed modelled values, the phase residual at said one or more selected frequencies, the phase residual of a discrete data point comprising the mismatch between the phase of the windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point; h) generating a phase residual data set at said one or more selected frequencies and for said at least one physical parameter; i) identifying cycle-skipped data within said phase residual data set; and j) evaluating the proportion of cycle-skipped data within said phase residual data set to determine the validity of the starting model for interpreting said measured seismic data set.

In one embodiment, step j) comprises calculating the pioportion of cycle-skipped data as a function of the total data in the phase residual data set.

In one embodiment, the starting model is rejected if the proportion of cycle-skipped data exceeds a pre-determined threshold.

In one embodiment, the pre-determined threshold is 40%.

In one embodiment, the pre-determined threshold is 20%.

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

In one embodiment, step n) further comprises determining the change in proportion of cycle-skipped data between said phase residual set and said further phase residual set.

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

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

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

In one embodiment, if, in step t), the preportion of cycle-skipped data increases from said phase residual set to said further phase residual set, the method further comprises: u) rejecting said fbrther starting modeL In one embodiment, if, in step t), the preportion of cycle-skipped data decreases from said phase residual set to said further phase residual set, the method further comprises: v) rejecting said starting model.

In one embodiment, the discrete data points of said measured and modelled seismic data sets are a function of souzve and receiver locations.

In one embodiment, the method thiTher comprises the step of w) confirming the accuracy of an optimised subsurfitce model generated using said starting model.

In one embodiment, the method further comprises the step of: x) confirming or rejecting the validity of said sub surface starting model based upon stepj).

In one embodiment, if it is confirmed in step x) that the optimised subsurface model is sufficiently accurate, the method further comprises the step: y) utilising said optimised model for sub-surface exploration.

In one embodiment, said at least one physical parameter is pressure.

In one embodiment, said at least one physical parameter comprises particle velocity or displacement.

In one embodiment, the measured data set and the modelled data set comprise data taken at a plurality of discrctc frequcncics.

In one embodiment, the measured data set and the modelled data set comprise values of a plurality of physical parameters.

In one embodiment, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a full waveform inversion iterative process.

According to a second aspect of the present invention, there is provided a 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 third 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 thc basic starting model of Figurc 2 for an individual seismic shot; Figure 4 a) is a synthetic example of measured seismic frace 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 e) together with the phase residual after 4 iterations of a standard fill 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) -0 illustrates the effect of starting model on subsequent iterations; Figure 8 is a flow chart illustrating the method according to an embodiment of the present invention; Figures 9a) -d) are schematic illustrations of windowing modelled seismic trace data; Figure I Oa) -d) are schematic illustrations of windowing measured seismic trace data; Figures II a) to i) illustrate phase residual plots and resulting models for different starting models; and Figurc 12 is a flow chart illustrating thc method according to another cmbodimcnt of thc prcscnt invention; As dcscribcd above, FWI analysis can cxtract many physical propertics (attcnuation, density, anisotropy, elastic) from thc parametcrs of an accuratc Earth model.

However, the essential parameter is VI,, the P-wave velocity. The present invention is directed towards a method for verifying the accuracy of a V model of the subsurface of a portion of the Earth generated by a FWI method.

As described above, FWII is an effective procedure for generating a quantitative detailed model of the subsurface for use in the exploration process. However it is a localised optimisation and thereforc is strongly dcpendent on wheit thc optirnisation begins from.

The present method provides a tool to assess whether a given starting model is adequate for FWI on a particular seismic dataset i.e. whether it will converge successthlly or not. This is done by highlighting the cycle-skipped modelled data.

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

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

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

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

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

Step 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 I, 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 UPS devices. Additionally, the location of the source 12 is also well known by location tracking means.

The measured seismic data set may comprise multiple "shots" or source 12 emissions. An example of seismic trace data for a single shot is shown in Figure 4 a).

The gathered data is essentially a plot of pressure as a function of receiver position (on the x-axis) and time (on the y-axis). This is because, in general, a detector such as a hydrophone measures the scalar pressure at its location.

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

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

The actual gathering of the seismic data set is described here for clarity.

However, this is not to be taken as limiting and the gathering of the data may or may not form part of the present invention. The present invention simply requires a real-world measured seismic data set upon which analysis can be performed to fitcilitate subsurfhce exploration of a portion of the Earth. This may be supplied by a third party.

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

The method now proceeds to step 102.

Step 102: Obtain starting model At step 102, an initial starting model of the specified subsurface portion ofthc Earth is, optionally, provided, potentially by a third party. Thc model may bc providcd in either a two dimensional or a three dimensional form. Whilst the illustrated examples are of two dimensional form, the skilled person would be readily aware that the present invention is applicable to three dimensional approaches. Examples of typical starting models are shown in Figure 7.

Thc generated model consists of values of the physical parameter V, 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 pcrson.

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

The method then proceeds to step 104.

Step /04: Obtain starling modelled data set.

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

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

In othcr words, thc modcllcd data sct corrcsponds discrctc point to discrctc point to the measured dataset. The modelled data set may be generated in the time domain in ordcr to pcrform timc domain analysis with rcspcct to thc mcasurcd data sct obtaincd in stcp 100.

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. modcllcd pressure at the location of the measured data points as a ftinction of position and time.

Figure 9 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 9a) corresponds to a modelled estimate of the time evolution of pressure at a particular receiver number or region and corresponds to a time slice taken through a data set such as that shown in Figure 4a). As shown in Figure 9a), there is a peak at early times (early arrival pressure waves) with fUrther oscillations at longer times.

The inventor has found that the effect of applying a time-domain window to the modelled and measured seismic data sets is operable to improve the msults of the validation check.

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 Figurc 11 b). Thc Gaussian function 152 has various parameters which can be specified by a user such as thc tiinc (r) at which the function or window 152 is centred and the Full Width Half Maximum (FWHM) which spccifics thc spread of the function.

However, whilst a Gaussian operator has been specified as the function 152 in this embodiment, other suitable fbnctions 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 9c) shows the function 152 superimposed on the original modelled data prior to multiplication. Thc effect of the function 152 as shown will bc 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 9d) 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 appzuach 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 perfonued.

Thc diffcrcnt frcqucncy componcnts can bc cxtractcd through thc usc of a Fouricr transform if so desired. Since only a few frequency components may be required, a discrctc Fourier transform may bc computed from thc inputtcd seismic tracc data.

Alternatively, a conventional thst Fourier transform (FF1) approach may bc uscd. As a further alternative, partial or fill Fourier transforms may be used.

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

The modelled data set, a (r, s), comprises data at one or more frequencies. These frequencies correspond to that of the measured data set dfr, 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 ftmnction 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 tO a) shows a graph illustrating an aspect of the raw trace data 160. This figure shows a plot of pressure (in arbitrary units) at the source position as a function of time (in arbitrary units, which may be seconds). This graph corresponds to the time evolution of pressure at a particular receiver number or region and corresponds to a time slice taken through a data set such as that shown in Figure 4b). As shown in Figure I Oa), there is a pcak at early times (early arrival pressure waves) with a complicated series of oscillations cxtending out to longer times due to later arrivals.

Thc inventor has found that the effect of applying a time-domain window to the modelled and measured scismic data sets is operable to improve the itsults of thc validation check.

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 (i) at which the frmnction or window 162 is centred and the Full Width Half Maximum (FVsTHM) which specifies the spread of the function. The value of i 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. There is no need to pick first breaks in the measured data, which is a time-consuming and error-prone process.

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

In general, the same window (i.e. the same function) will be used as is applied to the modelled data in step 106.

Figure tOe) 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 oft may be selected from thc first arrival timc derived by ray tracing in othcr to centre the window over the first arrivals. However, the window ccntrc t of the function 152 may alternatively be selected to be at later times to fbcus the FWI on later arrivals if dcsiicd. Also r can simply be user-specified as an increasing function of offict.

The measured raw data 160 is then combined with the Gaussian function 162, in this example by multiplication. Figure lOd) 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 methodthen proceeds to step 112.

Step 112: Transfonn 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 lbr 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 ibrm one or more windowed measured data sets with the or each data set tepresenting the windowed seismic traces at a selected frequency.

Tn 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 associatcd cocfficicnts arc well known in thc art. Since only a few frequency components may bc required, a discrete Fouricr transform may bc 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: Determine phase residual As set out above, the measured seismic data is expressed as a function of receiver (r) and source (s) pair position. This is the fUnction dfr,s). The modelled seismic data is also expressed as a ifinction of receiver (ñ and source (s) pair position. This, as set out above, is the function uo(r,s) for the starting modelled data set. The modelled and measured data sets a 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 variable is obtained for the modelled data set, i.e. between the data sets d(rc) and uo(r,s).

The principal value (between -it and + it) of the phase residual is denoted as and expresscd in equation 14) for an iteration n: 14) cji(r4 phase (ufr,s)d*fr,s)) which is the wrapped difference of the individual phases of ufr,s) and dfr,s). The wrapped difference is the difference between the phases when the phase difference is restricted to the range -it to + it. In other words, if the phae difference is it and increases, the value of the phase difference will become -it and move back towards it. This can readily be visualised on an Argand diagram.

The phase residual 1(r,s) can be determined for the data points to generate a phase residual data set. In that case, each spatial data point now has associated therewith a value of the phase residual (r,s) for an iteration n. Note that the starting modelled data set, uofr,s) is generated prior to the first iteration (for which n =1).

In the ease of the starting modelled data set uo(r,s), do(r,s) can be calculated from equation 14).

The method then proceeds to step 116, Step 116: Identifj' Cycle-Skipped Data in Modelled Dataset Once the phase residual data set is generated in step 114, the cycle skipped data therein can be identified. Figure 11 shows this analysis. Figures 11 a) and 11 b) illustrate a suitable and a poor starting model respectively. Figures llc) and lid) illustrate phase residual data sets for the models of Figures 1 la) and 1 lb) respectively.

11(r,s) is plotted as a function of source and receiver in Figures 1 ic) and lid).

As shown in Figure 1 ic), 11fr,s) is plotted at 3.125Hz and varics smoothly in spacc.

Whcrc thc valucs of(rc) arc negative, the modelled data is early relative to the observed data and where positive, the values of 0(r,s), the modelled data is late relative to the observed data.

11(r,s) is close to zero at zero-offset and in the case of the starting model of Figure 11 a) remains well within the-n to + n range. Also of note is the symmetry in the zero-offset axis which occurs through source-receiver reciprocity.

In the case of the second starting model shown in Figure 1 lb), c11fr,s) as plotted in Figure lid) is similarly smoothly varying but in general with a larger magnitude than with the second starting model. At a certain offset its value decreases sufficiently to reach -n. Instead of passing through-n, it abruptly switches to +n before decreasing again.

This 2n transition in c(r,s) forms a continuous boundary through the data. It marks where the phase difference between the modelled and observed data exceeds half a cycle.

In this case it marks where the modelled data is more than half a cycle early relative to the observed data. This shows how the well-known cycle-skipping' relationship between predicted and observed data manifests itself in the low frequency phase residual.

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

The cycle skipped region can then be identified numerically or graphically for further computational processing and analysis. The cycle skipped region can be identified by the cycle-skip boundary, with data beyond the cycle skip boundary when moving away from zero-offset being cycle-skipped. In Fig lid, the cycle-skipped data is the data in the top right i.e. on the side with a + it'. Similarly, the cycle skipped data is the data beyond the boundary in the boftom left. In Fig lie there is no cycle-skipped data.

As set out above, the cycle-skipped boundaries can bc idcntificd numcrically and computationally and thc cycle-skippcd regions identified thcrefrom.

Thc method thcn procccds to stcp 118.

Step 118: Evaluate cycle-skipped data At step 118, the cycle-skipped data is evaluated. This may be done in a number of suitable ways. For example, the proportion of cycle-skipped data within the phase residual data set may be evaluated. This may be expressed as a percentage score wherein the lower the percentage score, the higher the accuracy of the starting model and the greater the probability of accurate convergence.

A parameter may be assigned. For example, a starting model may be considered to be suitable if less than 20% of the data is cycle-skipped at 3Hz for example. However, other examples may be used, e.g. 30% or 40% depending upon the application.

Other methods of numerical analysis may be used. The distribution of cycle-skipped data may be utilised.

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

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

A second embodiment of thc present invention will now be described with refrrence to Figure 12. The second embodiment includes all of the method steps 100 to above.

Steps 100 to 120 above enable identification, qualification and quantification of misconvergence which may occur from an inaccurate starting model. An analysis of the phase residual can then be performed to determine the accuracy of thc starting model.

In the second embodiment, this process is enhanced by carrying out a tbrther iteration ofthe model to see how the phase residual changes.

Step 122: Obtain further modelled data set The data obtained in step 104 comprises modelled seismic trace data in the time domain (i.e. modelled pressure at the location of the measured data points as a function of position and time).

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

In the case of an FWI method, a trial velocity is updated to minimise an objective function measuring the mismatch between modelled and recorded data. The objective function and the steepest descent direction used to update the model in conventional FWI ii is dcscribcd in cquations 1) to 12). This can bc donc in thc timc or frcqucncy domain, although time domain data is rcquircd for thc analysis.

A further modelled data sct in thc timc domain obtained from an itcration of thc mcthod dcscribcd above is thcn providcd.

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

Repeat Stepsl0& 108, 114-118 The method then proceeds, using the further modelled data set, through steps 106 and 108 to obtain uj fr,s).

Steps 110 and 112, having already been carried out, maybe omitted at this stage.

The method then proceeds to analyse the further modelled data set a1 (r,s) at steps 114, 116 and 118.

However, in this embodiment, the evaluation at step 118 comprises comparing the two phase residual sets.

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

For example, if the proportion of cycle-skipped data in the further modelled data set increases when compared to the starting modelled data set, the model will not result in an accurate convergence and should be rejected.

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

Thus, in this embodiment, in step 118 the change phase residual Ap is defined.

The change phase residual Ap is thus defined in equation 15) for the ease where the starting modelled data set uofr,s) and the first iteration (further) modelled data set uifr,s) are used: 15) A =phase ((u0d*) (u1d*)*)= pha.se(u0uj*,) or, in a more general form (with dependency on source and receiver locations shown) for an iteration n: 16) A1 (r4= phase (u,, fr,s) u, i jfr,s) *) Equations 15) and 16) describe the change in phase residual caused by a model update. Ac (r,s) is a wrapped quantity lying between -n and + 71. Further it is not explicitly dependent on d so is noise-free.

Alternatively, Ac, (r,s) can be calculated directly from u0(r,s) and Ui fr,s) and so it is not explicitly necessary to calculate or obtain the phase residual for iteration 1 I (r,s)) directly, although this may be done if required.

Ac]3 (r,s) may be obtained at one or many frequencies. These frequencies need not be those used in the inversion procedure.

An analysis of the change phase residual can provide further detailed information to vindicate or reject a particular starting model. This is an advantagc of thc present method which enables thc divcrgcncc or misconvergence of thc data to be obscrvcd as a function of source and receiver position. This is in contrast to any known methods which may simply sum the crrors over the whole dataset.

Step 124: Finish When, at step 124, it has been determined that the convcrgence 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 ehairnels 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 athantages over prior art methods. As described above, here we show how to validate a FM/Il approach enabling more robust and reliable convergence when using this technique. We can identif' cases where conventional FWI has or will be sufficient for a given starting model

which does not exist in the prior art.

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

Variations in the method will be apparent to the skilled person. In addition, alternative or additional steps may be implemented with respect to the present invention.

For example, whilst the present invention has been described with reference to a model which involves solving the acoustic wave equation, the present invention is equally applicable to the models which involve solving the elastic wave equation.

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

Further, whilst the example here has used the scalar parameter of pressure as its focus (i.e. P-waves), ft is also possible (with appropriate equipment such as geophones) to measure the particle motion in the receiver location and so determine S-wave parameters.

Data so-generated could then be utilised and processed in the present invention.

Embodiments of the present invention are operable to identii' the validity of a starting model without recourse to techniques such as comparing the change residual of each two or more iterations with the check residual variable of the nJb iteration.

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, ft should be understood, however, that the drawings and detailed description are not intended to limit the invention to the particular tbnn disclosed. Tt will be appreciated that variations and modifications may be made to the examples described within the scope of the present invention.

Claims (26)

1. CLAIMS1. A method for checking the validity of a geophysical representation of a portion of the volume of thc Earth generated from a seismic measurement, the method comprising: a) receiving a measured seismic data set obtained from said seismic measumment of said portion of the volume of the Earth, the measured seismic data set comprising measured time-domain values of at least one physical parameter measured at a plurality of discrete data points; b) receiving a starting modelled seismic data set generated from a subsurface starting model, said modelled data set comprising modelled time-domain values of at least one physical parameter at the same plurality of discrete data points; c) applying a time-domain window to said modelled time-domain values of said modelled seismic data set to generate a windowed modelled seismic data set; d) performing a Fourier transform on said windowed modelled seismic data set to obtain windowed modelled values taken at one or more selected frequencies, said windowed modelled values comprising amplitude and phase components; e) applying a time-domain window to said measured time-domain values of said measured seismic data set to generate a windowed measured seismic data set; 0 performing a Fourier transform on said windowed seismic measured data set to obtain windowed measured values taken at one or more selected frequencies, said windowed measured values comprising amplitude and phase components; wherein the method ftrther comprises determining the validity of said subsurface starting model for performing full wavefonii inversion on said measured seismic data set by: g) calculating, from the phase component of said windowed measured values and the phase component of said windowed modelled values, the phase residual at said one or more selected frequencies, the phase residual of a discrete data point comprising the mismatch between the phase of the windowed measured values and the windowed modelled values of the at least one physical parameter for that discrete data point; h) generating a phase residual data set at said one or more selected frequencies and for said at least one physical parameter; i) identifying cycle-skipped data within said phase residual data set; and j) evaluating the proportion of cycle-skipped data within said phase residual data set to determine the validity of the starting model for interpreting said measured seismic data set.
2. 2. A method according to claim 1, wherein step j) comprises calculating the proportion of cycle-skipped data as a function of the total data in the phase residual data set.
3. 3. A method according to claim 2, wherein the starting model is rejected if the proportion of cycle-skipped data exceeds a pre-determined threshold.
4. 4. A method according to claim 3, wherein the pie-determined threshold is 40%.
5. 5. A method according to claim 4, wherein the pie-determined threshold is 20%.
6. 6. A method according to any one of the preceding claims, further comprising: k) repeating steps b) to d) for a further modelled seismic data set generated from an nth iteration (where n is an integer) of a full waveform inversion process carried out on said modelled seismic data set to generate further windowed modelled values comprising amplitude and phase components; 1) repeating steps g) to h) utilising said further windowed modelled values to generate a further phase residual data set the phase residual at said one or more selected frequencies; m) identifying cycle-skipped data within said further phase residual data set; and wherein step j) further comprises: n) comparing the proportion of cycle skipped data within said further phase residual data set to said proportion of cycle-skipped data within said phase residual set to determine the validity of the starting model for interpreting said measured seismic data set.
7. 7. A method according to claim 6, wherein step n) further comprises determining the change in proportion of cycle-skipped data between said phase residual set and said further phase residual set.
8. 8. A method according to claim 7, wherein if; instep n), the proportion of cycle-skipped data increases from said phase residual set to said further phase residual set, the method further comprises: o) rejecting said starting model.
9. 9. A method according to claim 7, wherein if; in step n), the proportion of cycle-skipped data decreases from said phase residual set to said further phase residual set, the method further comprises: p) accepting said starting model or generating further iterations of said modelled data set.
10. 10. A method according to any one of claims ito 6, further comprising: q) repeating steps b) to d) for a further modelled seismic data set generated from a different starting model to generate further windowed modelled values comprising amplitude and phase components; r) repeating steps g) to h) utilising said further windowed modelled values to generate a further phase residual data set the phase residual at said one or more selected frequencies; s) identif'ing cycle-skipped data within said further phase residual data set; and wherein stepj) fUrther comprises: t) comparing the proportion of cycle skipped data within said further phase residual data set to said proportion of cycle-skipped data within said phase residual set to detennine which of said stalling model and said further starting model is most suitable for interpreting said measured seismic data set.
11. II. A method according to claim 10, wherein step t) further comprises determining thc change in proportion of cycle-skipped data between said phase residual sot and said further phasc residual sct.
12. 12. A mcthod according to claim 11, wherein if, in stcp t), the proportion of cycle-skipped data increases from said phasc residual sct to said further phasc residual set, thc method further comprises: u) rejecting said further starting model.
13. 13. A method according to claim 11, wherein if, in step t), the proportion of cycle-skipped data decreases from said phase residual set to said further phase residual set, the method further comprises: v) rejecting said starting model.
14. 14. A method according to any one of the preceding claims, wherein the discrete data points of said measured and modelled seismic data sets are a function of source and receiver locations.
15. 15. A method according to claim 1, further comprising the step of w) confirming the accuracy of an optimised subsurface model generated using said starting model.
16. 16. A method according to any one of the preceding claims, further comprising the step of: x) confirming or rejecting the validity of said sub surface starting model based upon step j).
17. I?. A method according to claim 16, wherein if it is confirmed in step x) that the optimised subsurface model is sufficiently accurate, the method further comprises the step: y) utilising said optimised model for sub-surface exploration.
18. 18. A mcthod according to any onc of thc prcccding claims, wherein said at least one physical parameter is pressure.
19. 19. A mcthod according to any onc of thc prcccding claims, whcrcin said at lcast onc physical paramctcr compriscs particlc velocity or displacemcnt.
20. 20. A method according to any one of the preceding claims, wherein the measured data set and the modelled data set comprise data taken at a plurality of discrete frequencies.
21. 21. A method according to any one of the preceding claims, wherein the measured data set and the modelled data set comprise values of a plurality of physical parameters.
22. 22. A method according to any one of the preceding claims, wherein, subsequent to step a), the method further comprises: receiving a subsurface starting model representing a portion of the volume of the Earth, the subsurface starting model being for use in a lull waveform inversion iterative process.
23. 23. A computer program product executable by a programmed or programmable processing apparatus, comprising one or more sofiware portions for performing the steps of any one of claims 1 to 22.
24. 24. A computer usable storage medium having a computer program product according to claim 23 stored thereon.
25. 25. A method substantially as shown in and/or described with reference to any one or more of Figures I to 12 of the accompanying drawings.
26. 26. A computer program product substantially as described with reference to any one or more of Figures I to 12 of the accompanying drawings.