GB2481270A

GB2481270A - Method of, and apparatus for, full waveform inversion modelling

Info

Publication number: GB2481270A
Application number: GB1101345.5A
Authority: GB
Inventors: Nikhil Shah
Original assignee: Individual
Current assignee: Individual
Priority date: 2011-01-26
Filing date: 2011-01-26
Publication date: 2011-12-21
Anticipated expiration: 2031-01-26
Also published as: GB2481270B; GB201101345D0

Abstract

There is provided a method for subsurface exploration comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement. The method comprises: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 values of one or more frequency components of at least one physical parameter measured at a plurality of discrete data points and taken at one or more selected frequencies; b)generating, using a subsurface model of a portion of the Earth, at least one modelled seismic data set comprising modelled values of one or more frequency components of at least one physical parameter at the same plurality of discrete data points and at the same one or more selected frequencies; c) calculating, for a plurality of groups of two discrete data points, the principal value of the difference in phase mismatch 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 mismatch difference associated therewith, where the phase mismatch of a discrete data point comprises the mismatch between the phase of the measured and modelled values of the at least one physical parameter for that discrete data point;d) generating, from the phase mismatch difference of at least some of said groups, a phase mismatch difference data set at said one or more selected frequencies and for said at least one physical parameter; e) minimising the principal value of the phase mismatch difference data set to produce an updated subsurface model; and 0 utilising said updated model for subsurface exploration.

Description

Method of, and Apparatus for, Full Waveform Inversion Modelling The present invention relates to a method of, and apparatus for, full waveform inversion modelling of a portion of the Earth. More particularly, the present invention relates to a method of, and apparatus for full waveform inversion modelling of measured seismic data to generate a subsurface model of at least a portion of the Earth.

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

A general schematic illustration of an experimental set up 10 for an undersea seismic survey is shown in Figure 1. However, an equivalent experiment can be carried out on land. The present invention is applicable to subsurface exploration in any suitable environment, for example land or sea-based measurements of a portion of the subsurface of the Earth. However, the present invention is applicable in general to the subsurface exploration of any planetaiy 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 vaiy 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, in a particular medium. V, is a scalar quantity and is, effectively, the speed of sound in a medium. It is this quantity V which is of particular interest in seismic modelling.

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, each known as a seismic trace. For a two dimensional survey, the tens of thousands of individual traces may be taken. For the three dimensional case, this number may run into the millions.

A seismic trace comprises a measurement, by the multiplicity of detectors 16, of the returning reflected or refracted acoustic waves 30 originating from the source 12. In general, a partial reflection of the acoustic wave 20 occurs at a boundary or interface between two dissimilar materials, or when the elastic properties of a material changes.

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

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

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

Key to this process of modelling and imaging a portion of earth is seismic velocity Vi,. In of a portion of the volume of the Earth it 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 vaiy 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 aftempts to modif 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.

FWT can extract many physical properties (V and V velocities, aftenuation, density, anisotropy) of the modelled portion of the Earth. However, V, the P-wave velocity, is a particularly important parameter which the subsequent construction of the 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 apparent to the skilled person.

The FWT 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 full two-way wave equation. This is known as forward problem. This can be in the time domain or the frequency domain, elastic or acoustic, isotropic or anisotropic models can be used. 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 1 5 later) A modelled seismic single shot gather is shown in Figure 3. The modelled seismic trace has 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 ofJhll-wavejbrrn 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 fijnction through repeated steepest-descent direction calculation. An objective fijnction 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 flinction 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.

II 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 is makes it advantageous to extract the lowest useable frequencies in the recorded data and inverting these first before moving to successive higher frequencies. 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 fhnction. The shot is located at point s and w 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.

1) E=4Yu(r,s)-d(r,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, usually by finite-difference methods, by solving the frequency domain (uniform-density) acoustic wave equation 2) (V2 + oi2m2 (x)) u(x, s) = -Sö(x -s) here the subsurface model mx) is slowness, the inverse of V (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) ux, s)=rSGx, s).

where G is the green's function which we calculate by solving equation (2) with S1.

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) mx') = Re u(r,s) (u(r,s) -d(r,s))* Then, from Equation 2): (V2 + w2m2(x)) = -2w2m(x')u(x',s)6(x -x') 5) 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 wavefield emitted by a source positioned at x', given by equation 6) 6) = 2w2m(x')u(x',s)G(x,x') Next, the solution is evaluated at the receiver location, given by equation 7): =2ü2G(r,x')m(x')u(x',s) Then, applying source-receiver reciprocity the wavefield is now emitted from the receiver location, as specified in equation 8): 8) = 2w2rn(x')u(x',s)G(x',r) Inserting equation 8) into equation 4) gives equation 9) for single receiver: = 2ü2m(x') Re(u(x', s)w* (x', r, s)) where: 10) w(x',r,s) = G*(xf,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 11): 11) = 2oj2m(x)Re(u(x,s)w*(x,s)) where: 12) w(x',s) =G(x',r)(u(r,s)-d(r,s)) 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.

Ideally, the above method will lead to a convergence to a model which is a coffect 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 mismatch between d(r,s) and u(r,s) is determined. The principal value (between -it and + it) of the phase mismatch is denoted as and expressed in equation 13): 13) ep(r,s)= phase u(r,s)d*(r,s) 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 -it to + it. Tn 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 trace 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 ofet receiver location the phase mismatch 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 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 +ir transitions, can lead to convergence problems in numerical iterations.

At short offsets d2 (r) is not cycle-skipped. However due to the ±ir 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 4 is moving upwards to convergence, whereas in fact it should be moving downwards through -it 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 mis-converges to a cycle-skipped local minimum. However, for this condition to hold low frequency components or a particularly accurate starting model are required.

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. 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 stmctures of real-world portions of the Earth. Figures 6 d) to 6 f) illustrate the operation of FWJ 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 0 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.

However, alternative approaches have aftempted 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 at, Geophysical prospecting, 2007, 55, 449-464 and WO-A-2009/136708. In the above examples, the difference in the wavefield logarithms is minimised. This approach enables separation of amplitude and phase mismatch. II is well known that phase should be inverted first to match the kinematics of the recorded data before incorporating amplitudes into the procedure. Phase only inversion using this approach can minimise the unwrapped phase mismatch between the recorded and modelled data.

However, in practice estimating the number of cycle skips between the modelled and recorded data is problematic. This is particularly so in realistic situations when receiver and source locations will be a function of spatial dimensions (t,y,z) as opposed to a co-linear layout such as in Fig 4c. The difficulties arise in higher dimensions where the +1-it transitions will not in general uniquely determine the number of cycle-skips at each data point. Tn FWT the data-points are irregularly spaced making phase unwrapping especially difficult and prone to erroneous results.

As set out above, known FWI methods suffer from a technical problem of that cycle skipped mis-convergence can lead to significant enors in a final, generated model of a subsurface portion ofthe 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 values of at least one frequency component of at least one physical parameter measured at a plurality of discrete data points and taken at one or more selected frequencies; b) generating, using a subsurface model of a portion of the Earth, at least one modelled seismic data set comprising modelled values of at least one frequency component of at least one physical parameter at the same plurality of discrete data points and at the same one or more selected frequencies; c) calculating, for a plurality of groups of two discrete data points, the principal value of the difference in phase mismatch between the discrete data points in a group such that each group has a principal value of the phase mismatch difference associated therewith, where the phase mismatch of a discrete data point comprises the mismatch between the phase of the measured and modelled values of the at least one physical parameter for that discrete data point; d) generating, from the phase mismatch difference of at least some of said groups, a phase mismatch difference data set at said one or more selected frequencies and for said at least one physical parameter; e) minimising the principal value of the phase mismatch difference data set to produce an updated subsurface model; and futiIising said 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.

As a further advantage of the present method, cycle skipped data can be inverted without aftempting to determine the number of cycle-skips or phase unwrap the data as has been required in the prior art. The method of the present invention provides an automated and robust method for FWT in the frequency domain suitable for application to any common seismic data acquisition layout. This is because the method deals with differences in phase mismatch between points, not the phase mismatch itself For suitably closely-spaced receivers (or shots) this difference will on the whole remain in the principal value range between -it and + it whilst the actual phase mismatch becomes increasingly cycle-skipped with offset.

Further, working with differences between adjacent receivers (with a common source) means there is no need to estimate the source value. The difference quantity being minimised is independent of the source value. Therefore, the source value can be set to unity for forward modelling purposes.

In one embodiment, step a) comprises performing a Fourier transform on said seismic measurement data to obtain a plurality of frequency components of said data set, and selecting one or more frequencies therefrom. The frequencies at which the FT is taken with respect to may be real or complex. The measured data set contains frequency components of said measured seismic data at one or more frequencies.

In one embodiment, step b) wherein said modelled seismic data set at a selected frequency is generated directly from the subsurface model of a portion of the Earth.

In one embodiment, step b) comprises performing a Fourier transform on modelled seismic data generated from said subsurface model to obtain a plurality of frequency components of said data set, and selecting one or more frequencies therefrom.

In one embodiment, step c) comprises calculating the wrapped phase mismatch between the measured and modelled values at said plurality discrete data points such that each of said discrete data points has a wrapped phase mismatch value associated therewith, and calculating the phase mismatch difference therefrom; In one embodiment, step e) comprises:g) selecting / groups of discrete data points from said phase mismatch difference data set; h) defining an I dimensional variable b comprising said / groups; i) constructing an objective function from said variable b j) minimising said objective function by modifying at least one coefficient of said subsurface model of a portion of the Earth to provide an updated subsurface model of a portion of the Earth; and k) repeating steps b), c) and g) to j) for the updated subsurface model of a portion of the Earth until a convergence criterion is met.

In one embodiment, / comprises a subset of the total number of groups.

In one embodiment, / is of the order of the number of data points in the measured seismic data set.

In one embodiment, said plurality of discrete data points are identified by a source spatial location and a receiver spatial location and define a spatial separation between discrete data points in a group.

In one embodiment, the source spatial location comprises three coordinates and the receiver spatial location comprises three coordinates such that each discrete data point comprises six coordinate values.

In one embodiment, the I groups are selected on the basis that the spatial separation between discrete data points therein 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, the discrete data points within at least one group are spaced by no more than two discrete data points in any spatial direction.

In one embodiment, the discrete data points within at least one group are spatially adjacent one another. This may be within a non-uniform distribution of data points, or within a uniform distribution such as a grid.

In one embodiment, the discrete data points within at least one group are spatial nearest neighbours.

In one embodiment, each group comprises discrete data points having either a common source spatial location or a common receiver spatial location.

In one embodiment, if each group comprises discrete data points having a common source spatial location, the spatial separation between discrete data points in a group corresponds to the spatial separation between receivers.

In one embodiment, if each group comprises discrete data points having a common receiver spatial location, the spatial separation between discrete data points in a group corresponds to the spatial separation between sources.

In one embodiment, if all the groups in the phase mismatch difference data set fall within the spatial threshold, the 1 groups having the smallest values of wrapped phase difference are selected.

In one embodiment, the selected / groups have a wrapped phase difference less than a threshold value. In one embodiment, said threshold value is ir/2.

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

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

In one embodiment, said subsurface model of a portion of the Earth is updated by a localised gradient-based method.

In one embodiment, said at least one coefficient comprises the P-wave velocity, VP.

In one embodiment, said seismic measurement comprises a single source and said measured data points correspond to a plurality of receiver locations.

In one embodiment, said seismic measurement comprises a single receiver location and said measured data points correspond to a plurality of source locations.

In one embodiment, the measured data set, the modelled data set and the phase mismatch difference data set comprise data taken at a plurality of discrete frequencies.

In one embodiment, the measured data set, the modelled data set and the phase mismatch difference data set comprise values of a plurality ofphysical parameters.

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 the basic starting model of Figure 2for 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 mismatch between the measured data of Figure 4 a) and modelled data of Figure 4 b); Figure 5 shows the phase mismatch of Figure 4 c) together with the phase mismatch after 4 iterations of a standard full waveform inversion minimisation iteration; Figure 6 b) to 1) show examples of starting models and the resulting final models obtained from conventional full waveform inversion modelling; Figure 7 is a flow chart illustrating the method according to an embodiment of the present invention; Figure 8 is a schematic illustration of groups of two discrete data points according to an embodiment of the present invention; and Figure 9 illustrates two discrete data points and the phase mismatch difference therebetween.

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 being Vi,, the P-wave velocity. The present invention is directed towards a method for obtaining an accurate and high resolution V 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 7. Figure 7 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 (rj, 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: Transtbrin measured seismic data set The 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 necessaiy 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 measured data sets with the or each data set representing the seismic traces at a selected frequency.

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

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

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

Alternatively, these frequencies can be inverted simultaneously. d(r,s) is the notation we give for the measured data set when it has one or more measured physical parameters and one or more frequency components.

The method then proceeds to step 104.

Step 104: Provide starting model At step 104, 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, 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 106.

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

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

The isotropic acoustic constant-density version of this is equation 14) with shot location s: 14) (V2 +oj2 m2 (x))u(x,s) = -ö(x-s) where u (x, s) is the modelled data set at location x, m(x) is the subsurface model parameter, slowness (the inverse of Vu), and w is the angular value of the frequency being inverted.

As will be described later, the possible quantity being minimised (i.e. the difference in phase mismatch between receivers for a common source) is independent of the source value. This makes it sufficient to set S=1 in the acoustic wave equation.

15) ux,s)=Gx,s) where G is Green's function, x is the data point position and s is the source location.

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 data set u(r,s). The method now proceeds to step 108.

Step 108: Transtbrm modelled data set The 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.

However, alternatively, in common with step 104, 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 14) 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 106 is not needed.

The modelled data set, u (r, s), comprises data at one or more frequencies. These frequencies conespond to that of the measured data set d (r, s) such that meaningful analysis can be performed. Several frequencies can be inverted simultaneously.

The method now proceeds to step 110.

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

In step 110, the phase mismatch difference between d(r,s) and u(r,s) is determined. The principal value (between -it and + it) of the phase mismatch is denoted as and expressed in equation 16): 16) ep(r,s)= phase u(r,s)d*(r,s) 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 -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. This can readily be visualised on an Argand diagram.

The phase mismatch (r,s) can be determined for the data points. In that case, each spatial data point now has associated therewith a value of the phase mismatch (r,s).

However, in this embodiment, the phase mismatch difference can be calculated directly for data points (ri, Si) and (r2, S2) in accordance with expression 17): 17) phase u(r2,s2)d*(r2,s2)u *(ri, si)d(ri, Si) The phase mismatch 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 mismatch between each of the discrete data points in a group is calculated. A phase mismatch difference value is then associated with that group. The phase mismatch difference values then form a data set. The phase mismatch 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.

Spatial separation is related to both the source separation and receiver separation for the data points in the group. However, the present invention is particularly concerned with a set up where both data points within a group have some common source (or common receiver). In the case of a common source, the spatial separation of the group is the receiver spacing between the two data points.

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 8. In a set of data points 150, a number of groups 160 are shown. Any configuration of groups 160 may be envisaged, and one data point may form part of multiple groups as shown in Figure 8. As shown, the data points 150 are not uniformly arranged. This is because, in practice, the receiver locations are never evenly spaced due to experimental constraints. The phase mismatch difference data set is formed from the collection of groups each having a phase mismatch difference value.

In the present method, it is desirable that the spacing between data points 150 forming a group 160 does not exceed a threshold value. In the present embodiment, it is desirable that the spacing between two data points 150 in any group 160 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 150 in a group 160 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 160 which is arranged on a diagonal -i.e. the two data points arc 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 ftirther spaced apart. For example, a group may be formed from a discrete data point and a next-but-one data point. As a frirther alternative, any two data points may form a group.

By taking the phase mismatch difference between different discrete data points, we are effectively determining the gradient of a line between them. With reference to Figure 5, this could be interpreted to be the gradient of either trace shown therein.

However, in contrast to the simplified example shown in Figure 8, actual data obtained from seismic surveys is usually not usually regularly sampled due to experimental and practical limitations. Therefore, it is not always practical to take the gradient, and so the phase mismatch difference between adjacent discrete data points is used.

The reason for the use of the difference in phase mismatch between two discrete data points is illustrated with reference to Figure 9. Two data points are illustrated in Figure 9. One point has a phase mismatch of+170 degrees. The other has a phase mismatch of -170 degrees.

Under the new method the wrapped phase difference between the two points of 20 degrees is minimised. This enables us to avert cycle-skipped convergence without needing to consider the number of cycle-skips at each point.

In this example, only one component data (pressure) is used. Therefore, there is only a single value of ep(r,s) or of the phase mismatch difference between points.

However, in many situations, multiple variables may be used. For example, if pressure and particle velocities are measured, then there would be values of dp(r,s) for each direction of particle motion and one for pressure -leading to four separate 4j(r,s) values.

The phase mismatch difference for at least some of the groups 160 forms a phase mismatch data set as described above. The method now proceeds to step 112.

Step 112: Select groups of discrete data points Once all the value of the phase mismatch difference has been calculated for some or all of the groups of adjacent data points, and a phase mismatch difference data set generated, a subset of those groups is selected for frirther processing. Therefore, I groups (where 1 is a subset of the total number of groups) are selected.

The value of/is chosen to be of the magnitude of the total number of data points.

In a typical two dimensional survey, the number of data points is of the order of 100 x 100. In a three dimensional survey, the number of data points is of the order of 100 x 100 x 100. The two data points within a group, in one embodiment, either have a common source or a common receiver.

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 1 groups is primarily based on the spatial separation of the discrete data points. As discussed above, the criterion of twice the mimimum separation could be implemented, with the / 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 I. Additionally, should more than the desired / 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 mismatch difference.

The phase mismatch difference only provides meaningfil analysis and convergence if the phase mismatch values are less than a threshold value. In this embodiment, the threshold value is it/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 / of groups have the smallest values i.e. the values which are furthest away from the maximum possible values of +7-it.

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

Once the subset oft groups has been defined, the method proceeds to step 114.

Step 114: Define variables At step 114, the values of the phase mismatch difference in the / groups selected in step 112 are used to form an / dimensional vector b. The method then proceeds to step 116.

Step 116: Construct objective function In the present case, the present invention iteratively updates a trial velocity to minimise an objective function based on the norm of b such as the L-2 norm. Therefore an objective function is constructed based on equation 18) below: 18) E=--b12 Step 118. Minimise objective function and update model The model update derives from the gradient of(18) i.e. the partial derivative with respect to a point perturbation of the model iii(x) at each positionx'. In this example, groups are selected where both data points within the group have some common source.

It is sufficient to first consider one group labelled by index i with data points (ri,s) and (r2,s) at a single frequency.

The gradient for this term is set out in equation 19): _____ -((i,s) -�(i,s) Em(x') - 3m(x') 3m(x') ) i 19) To calculate this, we first note equation 20): 3(r,s) -j ( 1 öu(r,s) 20) m(x') -mU(TS) öm(x') Using equation 8) this gives 21) =2w2m(xP)Im('G(xP,r)) We then arrive at equation 22): 22) in() = 2co2m(x') Irn(u(x s)w*(x, r$, s)) where, as set out in equation 23): = G*(x,r1)) * (I) _G*(x,,) * (I) 23) u (r2,s) 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 residual source at receiver locations r1 and r2 The gradient for a whole shot gather comes from summing equation 21) over the all the pairs associated with that shot, as set out in equation 24) and 25): 24) " = 2o32m(xf)1114u(xf,s)w* (x', s)) 25) %t.(xs)=EG*(x,rhl) * (1) _G*(x,,hl) * (1) u (r2,s) u (i,s) The above is still calculated as a single wavefield with a multi-point source.

Hence the overall gradient calculation has the same computational structure as conventional FWJ.

The total gradient is then obtained by summing the gradient contributions from each individual shot and over a set of frequencies.

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

If the criteria 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 110 to 120 as discussed above.

Step 122: Finish When, at step 122, 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 havc 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 minimise the misfit or mismatch at each data point (r,s) individually. However, the present invention minimises the difference in misfit between groups or pairs of spatially adjacent (r,s) data points. In particular the wrapped-differences of the phase mismatch is minimised.

Therefore, in other words, the method of the present invention differs by not minimising the individual residuals themselves but the differences in their value, arising between pairs of spatially adjacent (r,s) data points. 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.

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.

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

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, 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 values of one or more frequency components of at least one physical parameter measured at a plurality of discrete data points and taken at one or more selected frequencies; b) generating, using a subsurface model of a portion of the Earth, at least one modelled seismic data set comprising modelled values of one or more frequency components of at least one physical parameter at the same plurality of discrete data points and at the same one or more selected frequencies; e) calculating, for a plurality of groups of two discrete data points, the principal value of the difference in phase mismatch 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 mismatch difference associated therewith, where the phase mismatch of a discrete data point comprises the mismatch between the phase of the measured and modelled values of the at least one physical parameter for that discrete data point; d) generating, from the phase mismatch difference of at least some of said groups, a phase mismatch difference data set at said one or more selected frequencies and for said at least one physical parameter; e) minimising the principal value of the phase mismatch difference data set to produce an updated subsurface model; and f) utilising said updated model for subsurface exploration.

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

3. A method according to claim I or 2, wherein step b) wherein said modelled seismic data set at a selected frequency is generated directly from the subsurface model of a portion of the Earth.

4. A method according to claim I or 2, wherein step b) comprises performing a Fourier transform on modelled seismic data generated from said subsurface model to obtain a plurality of frequency components of said data set, and selecting one or more frequencies therefrom.

5. A method according to claim 1, 2, 3 or 4, wherein step c) comprises calculating the principal value of the phase mismatch between the measured and modelled values at said plurality discrete data points such that each of said discrete data points has a principal value of the phase mismatch value associated therewith, and calculating the principal value of the phase mismatch difference therefrom;

6. A method according to any one of the preceding claims, wherein step e) comprises: g) selecting 1 groups of discrete data points from said phase mismatch difference data set; h) defining an 1 dimensional variable b comprising said 1 groups; i) constructing an objective function from said variable b j) minimising said objective function by modifying at least one coefficient of said subsurface model of a portion of the Earth to provide an updated subsurface model of a portion of the Earth; and k) repeating steps b), c) and g) to j) for the updated subsurface model of a portion of the Earth until a convergence criterion is met.

7. A method according to claim 6, wherein I comprises a subset of the total number of groups.

8. A method according to claim 7, wherein I is of the order of the number of data points in the measured seismic data set.

9. A method according to claim 6, 7 or 8, wherein said plurality of discrete data points are identified by a source spatial location and a receiver spatial location and define a spatial separation between discrete data points in a group.

10. A method according to claim 9, wherein the souite spatial location comprises three coordinates and the receiver spatial location comprises three coordinates such that each discrete data point comprises up to six coordinate values.

II. A method according to claim 9 or 10, wherein the 1 groups are selected on the basis that the spatial separation between discrete data points therein does not exceed a spatial threshold value.

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

13. A method according to any one of claims 9 to 12, wherein the discrete data points within at least one group are spaced by no more than two discrete data points in any spatial direction.

14. A method according to claim 13, wherein the discrete data points within at least one group are spatially adjacent one another.

15. A method according to claim 14, wherein the discrete data points within at least one group are spatial nearest neighbours.

16. A method according to any one of claims 9 to 15, wherein each group comprises discrete data points having either a common source spatial location or a common receiver spatial location.

17. A method according to claim 16, wherein, if each group comprises discrete data points having a common source spatial location, the spatial separation between discrete data points in a group corresponds to the spatial separation between receivers.

18. A method according to claim 16, wherein, if each group comprises discrete data points having a common receiver spatial location, the spatial separation between discrete data points in a group corresponds to the spatial separation between sources.

19. A method according to claim 11, wherein if all the groups in the phase mismatch difference data set fall within the spatial threshold, the 1 groups having the smallest values of principal value of the phase difference are selected.

20. A method according to claim 19, wherein the selected 1 groups have a principal value of the phase difference less than a threshold value.

21. A method according to claim 20, wherein said threshold value is ir/2.

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

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 said subsurface model of a portion of the Earth is updated by a localised gradient-based method.

25. A method according to any one of claims 6 to 24, wherein said at least one coefficient comprises the Pwave velocity, VP.

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

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

28. 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 I to 27.

29. A computer usable storage medium having a computer program product according to claim 28 stored thereon.

30. A method substantially as shown in and/or described with reference to any one or more of Figures 1 to 9 of the accompanying drawings.

31. A computer program product substantially as described with reference to any one or more of Figures ito 9 of the accompanying drawings.