[go: up one dir, main page]

WO2016193180A1 - Procédé amélioré de modélisation d'inversion - Google Patents

Procédé amélioré de modélisation d'inversion Download PDF

Info

Publication number
WO2016193180A1
WO2016193180A1 PCT/EP2016/062092 EP2016062092W WO2016193180A1 WO 2016193180 A1 WO2016193180 A1 WO 2016193180A1 EP 2016062092 W EP2016062092 W EP 2016062092W WO 2016193180 A1 WO2016193180 A1 WO 2016193180A1
Authority
WO
WIPO (PCT)
Prior art keywords
data set
objective function
model
observed data
data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/EP2016/062092
Other languages
English (en)
Inventor
Mike Warner
Lluis Guasch
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
SUB SALT SOLUTIONS Ltd
Original Assignee
SUB SALT SOLUTIONS Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by SUB SALT SOLUTIONS Ltd filed Critical SUB SALT SOLUTIONS Ltd
Publication of WO2016193180A1 publication Critical patent/WO2016193180A1/fr
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/66—Subsurface modeling

Definitions

  • the present invention relates to a method of, and apparatus for, inversion modelling. More particularly, the present invention relates to an improved methodology for inverse problems which enables a more accurate determination of step length. Additionally, the present invention relates to an improved method of, and apparatus for, inversion modelling which requires fewer computational resources than known methods.
  • Seismic surveys are the principal means by which the petroleum industry can explore the subsurface of the Earth for oil and gas reserves. Typically, seismic survey data is acquired and analysed with regard to identifying locations suitable for direct investigation of the sub-surface by drilling. Seismic surveying also has applications within the mining industry and within other industrial sectors that have an interest in details of the subsurface of the Earth.
  • one or more natural or artificial seismic sources are arranged to generate vibrational energy which is directed into the subsurface of the Earth. Reflected, refracted and other signals returned from subsurface features are then detected and analysed. These signals can be used to map the subsurface of the Earth.
  • FIG. 1 A schematic illustration of an experimental set up 10 for an undersea seismic survey is shown in Figure 1.
  • this example is intended to be non-limiting and 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.
  • the present invention may be applicable to identification of numerous subsurface resources, and is intended to include oil exploration and gas prospecting.
  • the experimental set up 10 comprises a source 12.
  • the source 12 is located on a ship 14, although this need not be the case and the source may be located on land, or within the sub-surface, or on any other suitable vessel or vehicle.
  • the source 12 generates acoustic and/or elastic waves 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, or alternatively an air gun or other mechanical device capable of creating sufficient vibrational disturbance. Commonly, for many seismic survey experiments a single source is used which is shot from multiple locations. Naturally occurring sources may also be employed.
  • a plurality of detectors 16 is provided.
  • the detectors 16 may comprise any suitable vibrational detection apparatus. Commonly, two types of device are used. Geophones which detect particle motion, and hydrophones which detect pressure variations. Commonly, a large number of detectors 16 are laid out in lines for 2D data acquisition. Alternatively, the detectors 16 can be arranged in sets of lines or in a grid for 3D data acquisition. Detectors 16 may also be located within the subsurface, for example down boreholes.
  • 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. However, this need not be the case and other arrangements are possible.
  • elastic waves 20 generated by the source 12 propagate 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 elastic waves 20 are transmitted and refracted through the layers and/or reflected off the interfaces between them and/or scattered from other heterogeneities in the sub-surface and a plurality of return signals 30 is detected by the detectors 16.
  • the returning signals 30 comprise elastic waves having different polarisations.
  • Primary or pressure waves (known as P-waves) are approximately longitudinally polarised and comprise alternating rarefactions and compressions in the medium in which the wave is travelling.
  • P-waves typically have the highest velocity and so are typically the first to be recorded.
  • P-waves travel at a velocity V p in a particular medium.
  • V p may vary with position, with direction of travel, with frequency, and with other parameters, and is, effectively, the speed of sound in a medium. It is this quantity V p which is most commonly of particular interest in seismic inversion.
  • Shear or secondary waves may also be generated.
  • S-waves have an approximately transverse polarisation. In other words, in an isotropic environment, the polarisation is perpendicular to the direction of propagation.
  • S-waves are in general, more slowly moving 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 typically composed of a large number of individual source excitation events. The Earth's response to these events is recorded at each receiver location, as a seismic trace for each source-receiver pair. For a two dimensional survey, the tens of thousands of individual traces may be taken. For the three dimensional case, this number may run into the millions.
  • a seismic trace comprises a sequence of measurements in time made by one or more of the multiplicity of detectors 16, of the returning reflected, refracted and/or scattered acoustic and/or elastic waves 30 originating from the source 12.
  • 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 sampled in time at discrete intervals of the order of milliseconds.
  • 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 mineral resources such as hydrocarbons (for example, oil and natural gas).
  • Mineral resources such as hydrocarbons (for example, oil and natural gas).
  • Features of interest in prospecting include: faults, folds, anticlines, unconformities, salt domes, reefs.
  • V p the seismic velocity
  • V p may be estimated in various ways.
  • Full-wavefield inversion is one known method for analysing seismic data.
  • FWI is able to produce models of physical properties such as V p in a subsurface region that have high fidelity and that are well resolved spatially.
  • FWI seeks to extract the properties of subsurface rocks from a given seismic dataset recorded at the surface or seabed.
  • a detailed velocity estimate can be produced using an accurate model with variations on the scale of a seismic wavelength.
  • the FWI technique involves generating a two or three dimensional model to represent the measured portion of the Earth and attempting to modify the properties and parameters of the Earth model to generate predicted data that matches the experimentally obtained seismic trace data.
  • the predicted data is calculated from the subsurface model typically 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.
  • V p and V s velocities, attenuation, density, anisotropy are particularly important parameter which the subsequent construction of the other parameters depends heavily upon.
  • 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.
  • 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.
  • the basic model shows a gradually increasing V p with increasing depth without any detailed structure.
  • 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 positions showing pressure recorded as a function of time.
  • the parameters of the model are estimated at a plurality of points set out in a grid or volume, but they may be estimated from any suitable parameterisation.
  • the model is used to generate a modelled representation of the seismic data set.
  • the modelled data set is then compared to the real-world experimentally obtained seismic data set.
  • the parameters of the model are modified until the modelled data set generated from the Earth model matches the actual observed seismic data to a sufficient degree of accuracy or until sufficient convergence is obtained. Examples of this technique are illustrated in "An overview of full-waveform inversion in exploration geophysics", J. Virieux and S. Operto, Geophysics Vol. 74 No. 6 and US-A-7, 725,266.
  • FWI typically operates on the principle of iteratively updating the starting model to minimise or maximise an objective function through repeated steepest-descent direction calculation, or an analogous technique.
  • An objective function represents some measure of the mismatch or some measure of similarity between the recorded data and the predicted data.
  • a measure of mismatch obtained for example by subtracting two traces, should be minimised; whereas a measure of similarity, obtained for example by cross-correlating two traces, should be maximised.
  • 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.
  • the objective function can be formulated in the frequency domain, the time domain or other suitable domain. The choice of domain allows the use of pre-conditioning on either the data or the model update direction that could improve convergence or the linearity of the inverse problem.
  • Frequency domain inversion is equivalent to time domain inversion if all the frequencies are inverted simultaneously. However, the global minimum broadens at lower frequencies reducing how accurate the starting model needs to be for localised inversion to be successful.
  • the residual expresses the misfit between the two datasets as a single number. This parameter is known as the objective function, although it often takes other names such as the objective function, the cost function or the functional.
  • the objective function f is a real, positive, scalar quantity, and it is a function of the model m.
  • d(m) is the data where d(m) is the data modelled using a model m and d 0 is the observed data.
  • FWI is a local iterative inversion scheme.
  • a starting model m 0 that is assumed to be sufficiently close to the true, ideal model is prepared.
  • a typical, generalised FWI-type method involves the following steps:
  • step length a controls the scaling of an update to be applied to the model. This is, clearly, a critical aspect of the inversion process.
  • m k+i m k + sSm ⁇
  • 5m k+ is the normalised update for the k th iteration of the model m k .
  • the step length a can be calculated in a number of different ways. Two possible approaches are as follows:
  • the value of the functional f is plotted as a function of a.
  • the thicker line shows the actual distribution of f as a function of a, and the thinner line the estimated parabola.
  • the parabola is estimated based on the two obtained a values, a gue ssi and a gue ss2-
  • the optimal step length is the value of a where f is at a minimum. As shown, this method can result in errors when only a few points are used to generate an estimated parabola. As shown, the estimated optimal step-length is shifted from the actual (real) optimal step-length.
  • the second approach also has disadvantages.
  • the residuals are 5d.
  • the Born approximation may not hold in all cases, leading to an inaccurate step length determination.
  • model-dependent Wiener filters w(m) are defined which enable the modelled data to be fitted to the measured data.
  • the method of AWI is described in United Kingdom Patent GB2509223B.
  • a method of subsurface exploration comprising generating a geophysical representation of a portion of the volume of the Earth from a seismic measurement of at least one physical parameter, the method comprising the steps of: providing an observed data set, the observed data set comprising data values derived from seismic measured values of said portion of the volume of the Earth; generating, using a subsurface model of a portion of the Earth comprising a plurality of model coefficients, a modelled data set comprising a plurality of modelled data values; updating the subsurface model by: generating an objective function operable to measure the mismatch or similarity between the observed data set and the predicted data set; determining, using a first proportion of the total number of data values of the observed data set, the gradient of the objective function; determining, using a reduced objective function utilising a second proportion of the total number of data values of the observed data set, the step length for a subsurface model update, wherein the second proportion comprises 40% or less of the total number of data values of
  • the second proportion is 30% or less of the total number of data values of the first proportion of the observed data set. In one embodiment, the second proportion is 15% or less of the total number of data values of the first proportion of the observed data set. In one embodiment, the second proportion is 10% or less of the total number of data values of the first proportion of the observed data set. In one embodiment, the second proportion is 5% or less of the total number of data values of the first proportion of the observed data set. In one embodiment, the first proportion is 100% of the total number of data values in the observed data set.
  • the observed data set comprises a plurality of shots, each shot comprising a plurality of seismic traces grouped by source or receiver location and comprising data values derived from seismic measured values of said portion of the volume of the Earth; and the modelled data set comprises a corresponding plurality of shots comprising a plurality of modelled traces grouped by source or receiver location and comprising modelled data values.
  • the first proportion of the observed data set comprises a first number of shots and the second proportion of the observed data set comprises a second number of shots, the second number being 40% or less than the first number.
  • the second number is 30% or less of the first number.
  • the second number is 15% or less of the first number.
  • the second number is 10% or less of the first number.
  • the second number is 5% or less of the first number.
  • step f) comprises determining the step length for a subsurface model update by obtaining a plurality of trial values of the step length and identifying the minimum step length therefrom.
  • the number of trial values of the step length is 5 or more. In one embodiment, the number of trial values of the step length is 10 or more. In one embodiment, the number of trial values of the step length is 50 or more.
  • step d) further comprises: generating at least one convolutional filter, which when convolved with the observed data, increases the similarity between the predicted data and the convolved observed data; generating at least one convolutional reference filter, which when convolved with the observed data, leaves the observed data unaltered to some level of precision; and generating an objective function operable to measure the similarity or mismatch between the filter coefficients for the or each convolutional filter and the filter coefficients for the or each convolutional reference filter.
  • step d) further comprises: generating at least one convolutional filter, which when convolved with the predicted data, increases the similarity between the observed data and the predicted data; generating at least one convolutional reference filter, which when convolved with the predicted data, leaves the predicted data unaltered to some level of precision; and generating an objective function operable to measure the similarity or mismatch between the filter coefficients for the or each convolutional filter and the filter coefficients for the or each convolutional reference filter.
  • the method comprises: generating at least one non-trivial convolutional filter, the or each filter comprising three or more non-zero filter coefficients; generating a convolved observed data set by convolving the or each convolutional filter with said observed seismic data set; generating one or more filter objective functions operable to measure the similarity and/or mismatch between said convolved observed dataset and said predicted dataset; maximising and/or minimising at least one of said filter objective functions by modifying at least one filter coefficient of the or each convolutional filter; and generating one or more pre-determined reference filters comprising at least three reference coefficients.
  • said at least one convolutional filter is operable to transform at least a portion of said observed data set to render the observed data set and predicted data set approximations of one another.
  • one or more of the convolutional filters comprise Wiener filters.
  • a plurality of observed data sets and/or a plurality of modelled data sets is provided.
  • the objective function and/or the reduced objective function comprise different functions. In one embodiment, the objective function and/or the reduced objective function comprises a norm misfit objective function.
  • the objective function and/or the reduced objective function comprises an l_i-norm misfit objective function.
  • the objective function and/or the reduced objective function comprises a least-squares misfit objective function.
  • at least a portion of the observed data set is numerically propagated to a subsurface region of the model that is spatially removed from the location at which they were originally recorded.
  • step g further comprising, subsequent to step g), repeating steps b) to f) for the updated subsurface model of a portion of the Earth until a convergence criterion is met.
  • step h) further comprises: utilising said updated model for sub-surface exploration.
  • said at least one physical parameter comprises pressure, particle velocity or displacement.
  • the observed data set and the predicted data set comprise values of a plurality of physical parameters.
  • 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 aspect.
  • Figure 1 is a schematic illustration of a typical seismic survey experiment in which seismic traces are obtained from an undersea portion of the Earth;
  • Figure 2 is a schematic illustration of a basic starting model for full waveform inversion modelling
  • Figure 3 is a schematic illustration of modelled seismic trace data generated from the basic starting model of Figure 2 for an individual seismic shot;
  • Figure 4 shows a graph of an objective function as a function of the step length a for a typical functional. Guess values of a are shown to enable the illustrated parabolic fit;
  • Figure 5 shows a seismic trace data shot for observed and modelled data sets as a function of time;
  • Figure 6 shows the misfit of the observed and modelled data sets in Figure 5 as a function of time for conventional FWI;
  • Figure 7 shows a method according to a first embodiment of the present invention
  • Figure 8 shows a method according to a second embodiment of the present invention
  • Figure 9 shows a graph of an objective function as a function of step length a for a full observed data set of 91 shots and for a reduced observed data set of 3 shots;
  • Figure 10 shows a trial model for testing the accuracy of FWI;
  • Figure 1 1 shows the convergence achieved for conventional FWI relative to said trial model
  • Figure 12 shows the convergence achieved for AWI relative to said trial model
  • Figure 13 shows a schematic representation of the variation of the objective function with V p and V s with different step sizes.
  • the present invention relates to an improved methodology for inverse problems which enables a more accurate determination of step length. Therefore, more accurate and reliable convergence can be obtained using the present invention when compared to known approaches. More particularly, the present invention provides a robust method to calculate the step-length in Full Waveform Inversion (FWI)-like methods where the residuals to be minimised do not have a linear relation with the model parameters.
  • FWI Full Waveform Inversion
  • the present invention is further operable to improve the convergence rate of conventional FWI by computing more accurate step-lengths. This, therefore, assists in eliminating the linear assumption commonly adopted by the scientific community that leads to step-length calculations relying on a parabolic fit of relatively a small number of data points.
  • the present invention provides a method whereby only a small subset of the original data or datasets is used to perform an explicit line search of the functional values in the update (preconditioned gradient) direction. When the number of points computed in the line search times the number of shots used to evaluate the functional value is equal to the total number of shots used in the inversion (or in any particular iteration) then the cost of both the conventional and our new method is exactly the same.
  • the following embodiments illustrate the application of the present invention in practice.
  • the first embodiment outlines the general approach of the present invention.
  • the second embodiment relates to a specific implementation of the present invention utilising AWI methods using Wiener filters.
  • Figure 7 shows a flow diagram of a first embodiment of the present invention.
  • Step 100 Obtain observed data set Initially, it is necessary to obtain a set of experimentally gathered data in order to initiate subsurface exploration. This may be gathered by an experimental arrangement such as the set up shown and described with reference to Figure 1.
  • the gathered seismic data may be optionally pre-processed in various ways including by propagating numerically to regions of the surface or subsurface where experimental data have not been acquired directly.
  • the skilled person would readily be able to design and undertake such pre-processing as might be necessary or desirable.
  • the resultant seismic dataset representing experimentally-gathered data is known as an "observed data set".
  • a large number of receivers or detectors 16 are positioned at well known positions on the surface of the portion of the Earth to be explored.
  • the detectors 16 may be arranged in a two dimensional (such as a line) or a three dimensional (such as a grid or plurality of lines) arrangement.
  • the physical location of the detectors 16 is known from, for example, location tracking devices such as GPS devices. Additionally, the location of the source 12 is also well known by similar location tracking means.
  • the observed data set generally multiple source 12 emission traces known in the art as "shots".
  • a typical observed data set may comprise hundreds of individual traces (or shots).
  • a shot is a group of observed data traces belonging to the same source or receiver.
  • a shot would be a group of observed data trace taken at different receiver locations for the same source.
  • a shot identified by individual receiver a shot would be a group of observed data trace taken at the same receiver location for different source locations.
  • Multiple shots comprise the "observed data set”.
  • Figure 5 shows an example of an observed data shot which forms part of the observed data set do.
  • Figure 5 also shows a predicted seismic data trace forming part of a predicted data set d(m). The predicted data set will be described in the next steps. Both traces are shown as a function of time (on the vertical axis).
  • the data comprises pressure as a function of receiver position (on the x-axis) with respect to time (on the y-axis). This is because, in general, a detector such as a hydrophone measures the scalar pressure at its location. However, other arrangements may be used.
  • the seismic trace data comprises a plurality of observed data points. In one example, each measured discrete data point has a minimum of seven associated location values - three spatial dimensions (x, y and z) for receiver (or detector) position (r), three spatial dimensions (x, y, z) for source location (s), and one temporal dimension measuring the time of observation relative to the time of source initiation, together with pressure magnitude data.
  • the seven coordinates for each discrete data point define its location in space and time.
  • the seismic trace data also comprises one or more measurement parameters which denote the physical property being measured.
  • a single measurement parameter, pressure is measured.
  • the observed data set is defined as d 0 and, in this embodiment, is in the time domain.
  • 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 observed data set upon which analysis can be performed to facilitate subsurface exploration of a portion of the Earth. The method now proceeds to step 102.
  • Step 102 Provide starting model At step 102, an initial starting model of the specified subsurface portion of the Earth is provided.
  • the model may be provided in either a two dimensional or a three dimensional form. Whilst the illustrated examples are of a two-dimensional form, the skilled person would be readily aware that the present invention is applicable to three dimensional approaches.
  • the form of the model is not material to the present invention, and may take any suitable form. However, a general example will be described.
  • a commonly generated model generally consists of values of the coefficient V p and, possibly, other physical values or coefficients, typically defined over a discrete grid representing the subsurface.
  • Such starting models are routinely generated and represent the general trends of the major features within the subsurface region to be modelled and could be readily generated by the skilled person.
  • Predicted seismic data may be generated based on an analysis of the acoustic isotropic two- way wave equation as set out below in equation 4):
  • the wave equation 4) represents a linear relationship between a wave field p and the source s that generates the wave field. After discretisation (with, for example, finite differences) we can therefore write equation 4) as a matrix equation 5):
  • equation 4) can be rewritten as equation 7):
  • G(m) p
  • m is a column vector that contains the model parameters. Commonly these will be the values of c (and p if density is an independent parameter) at every point in the model, but they may be any set of parameters that is sufficient to describe the model, for example slowness l i e , acoustic modulus c 2 p , or impedance cp .
  • G is not a matrix. Instead it is a non-linear Green's function that describes how to calculate a wavefield p given a model m.
  • Step 104 Generate predicted data set
  • a modelled data set is generated.
  • the predicted data is required to correspond to the same source-receiver location data positions as the actual measured trace data so that the modelled and observed data can be compared.
  • the predicted data set corresponds discrete point to discrete point to the observed dataset.
  • the predicted data set is generated for the same measurement parameter(s) at the same frequency or frequencies. From the above analysis, predicted seismic data can be generated for one or more physical parameters in the time domain. If done in the frequency domain it could be done for one or more selected frequencies. This forms the modelled data set d(m).
  • Figure 5 shows an predicted seismic data trace forming part of a predicted data set d(m). The method now proceeds to step 106.
  • Step 106 Construct objective function
  • an objective (or misfit) function is configured.
  • the objective function is configured to measure the dis-similarity between the modelled and observed data sets (or analogues thereof).
  • an objective function may be configured that measures similarity; in this case, step 108 will operable to maximise rather than minimise the objective function.
  • the objective function f is a real, positive, scalar quantity, and it is a function of the model m. In practice, a factor of a half is often included in the definition of the objective function f because it makes later results simpler.
  • d(m) is the data modelled using model m and d 0 is the observed data.
  • the difference 5d also known as the residual
  • a model is sought that minimises the L 2 -norm of the data residuals.
  • the L 2 -norm expresses the misfit between the two datasets as a single number.
  • This parameter is known as the objective function, although it often takes other names such as the misfit function, the cost function or the functional.
  • the objective function f is a real, positive, scalar quantity, and it is a function of the model m.
  • n s , n r and n t are the number of sources, receivers and time samples in the observed data set.
  • Figure 6 shows an exemplary objective function f corresponding to the data of Figure 5.
  • the value of the objective function, or misfit is oscillatory as a function of time as shown.
  • step 108 The skilled person would readily understand how to design such objective functions and how to minimise or maximise them.
  • the method then proceeds to step 108.
  • Step 108 Calculate gradient
  • FWI localised gradient-based methods
  • These methods iteratively update an existing model in a direction that derives from the objective function's direction of steepest descent.
  • the central purpose of FWI is to find a model of the subsurface that minimises the difference between an observed seismic dataset and a predicted seismic dataset generated by the model for the same real-world spatial data points as the observed seismic dataset.
  • FWI or AWI are local iterative inversion schemes. Improvements to the starting model are made which successively reduces the objective function towards zero. Therefore, across an iterative step of the calculation, the objective function needs to be considered for a starting model m 0 and a new model m, at iteration 1.
  • 5m k is the normalised update for the k* 1 iteration of the model m k which minimises the objective function, and is the a step length.
  • the first stage is to determine the normalised update by calculating the gradient. This may be done by any suitable method. However, a generic example of gradient of the objective function is obtained in accordance with known conventional FWI.
  • An expression can be derived which expresses the update to the model 5m k :
  • V m / is the gradient of the objective function f with respect to the model parameters
  • H is the Hessian matrix of second differentials, both evaluated at m k .
  • the gradient is a column vector of length M and the Hessian is an M x M symmetric matrix.
  • Equation 1 1 Equation 1 1
  • the forward wavefield p is calculated, the numerical operator A is differentiated with respect to the model parameters and the final term of equation 13) is calculated, which represents a back-propagated residual wavefield. These terms are then multiplied together for all times and all sources, and summed together to give a value corresponding to each parameter within the model, typically to give one value of the gradient at each grid point within the model.
  • a c5p ⁇ Equation 14) simply describes a wavefield p that is generated by a (virtual) source 5d, and that is propagated by the operator A T which is the adjoint of the operator in the original wave equation. So the term that we need to compute in equation 11 ) is just the solution of a modified wave equation with the data residuals used as a source.
  • the gradient for each of the seismic trace groups (or shots) forming the observed data set can then be obtained.
  • the gradient is determined for the total number (or a large subset thereof) of shots in the observed data set. This may, typically, be in the region of 100 shots for a single parameter analysis. However, if more than one parameter is used (e.g. Vp and Vs are calculated), then the number of shots increases. For a 3D computation, the number of shots may be around 10000.
  • Step 110 Determine step length
  • Step 1 10 the step length a is determined.
  • Step 108 enables determination of the gradient for each shot, summed across each shot.
  • the gradient identifies in which direction the model should be changed. However, it does not specify by how much. To determine this, the step length a is computed.
  • the value of the functional f is plotted as a function of a.
  • the thicker line shows the actual distribution of f as a function of a, and the thinner line the estimated parabola.
  • the parabola is estimated based on the two obtained a values, a gue ssi and a gue ss2-
  • the optimal step length is the value of a where f is at a minimum.
  • This approach is advantageous since it only requires one step of forward modelling to determine a.
  • the present invention is operable to determine the step length more accurately than for known methods.
  • each forward model run to determine the step length is carried using a reduced objective function utilising only a sub-set of the full observed data set.
  • a small sub-set of the full observed data set is used. For example, 40% or less of the data of the observed data set used for calculation of the gradient in the previous step may be used in the determination of step length in step 1 10. In embodiments, 15% or less of the data or shots) of the data of the observed data set used for calculation of the gradient in the previous step may be used in the determination of step length.
  • the sub-set of the observed data set used for determination of step length can be reduced below 10% to, for example, 5% or even lower.
  • the inventors have found that useful results can be obtained using only a very small sub-set of the observed data set used in the objective function for determining the gradient in step 108. For example, 2-4% of the observed data set used in the calculation of the gradient may be used to determine in the or each calculation of the step length.
  • the sub-set of the observed data set used to determine step length may be differentiated in terms of numbers of "shots". As described above, a shot is a group of observed data traces belonging to the same source or receiver. For the same source, a shot would be a group of observed data traces taken at different receiver locations for the same source. Conversely, for a shot identified by individual receiver, a shot would be a group of observed data trace taken at the same receiver location for different source locations.
  • the size of the data set is, at least in part, determined by the size of the observed data set. This is because, as described above, a predicted data set is generated based on the individual data points (i.e. source and receiver locations) in the observed data set. The differences between the two are then minimised. Therefore, the number of data points and number of sampled parameters in the observed data set defines the overall size of the data set for processing in the objective function.
  • composite shot data sets may be used with the method of the present invention.
  • the form and nature of the data set is not material to the present invention, provided that the data set used to determine the step length is reduced in size when compared to the data set used to determine the gradient in step 108. Indeed, the two data sets need not overlap.
  • a sub-set of the original observed data set to be used to calculate the gradient in step 108. This may be done either through appropriate selection of data or through processing such as composition of data shots or other pre-calculation processing.
  • a different sub-set of the original data may then be used for calculation of the step length, provided the absolute size of the data set used is smaller than the data set used for determination of the gradient.
  • Figure 9 shows an example of step length determination for an objective function using the AWI method.
  • An objective function computed using all shots i.e. the full observed data set
  • the objective function was computed using an observed data set comprising 91 shots.
  • a reduced objective function computed using only 3 shots is also shown. Note that the reduced objective function is identical to the objective function save for the use of a reduced number of shots in its calculation.
  • the two graphs vary in form, as may be expected.
  • the reduced objective function calculated using only three shots has the same global minimum in a as the full objective function calculated using 91 shots. This is the key parameter to determine the optimum step length a, and so this example shows that the reduced objective function calculation can provide an accurate determination of a from a reduced sub-set of shots.
  • the calculation of the reduced objective function is significantly less expensive computationally, and so the calculation can be completed with a reduced load on computational resources. It has been determined above that the computation of the reduced objective function can provide accurate information on the minimum value of a, albeit whilst requiring significantly fewer computational resources.
  • a further advantage of this approach is that a number of different computations of a can be carried out for the same or similar computational resource cost as a single conventional objective function determination of step length. Therefore, a can be sampled at a greater number of points along the objective function curve to obtain a more accurate value of the step length. This is particularly important for methods such as AWI, where the Born approximation may not hold under all circumstances. However, it also has significant implications for FWI-type methods where a local minimum in a may be found using conventional approaches, leading to convergence of the model to an inaccurate local minimum.
  • the number of sampled points is five or more. However, 50 to 100 samples could be used to sample the curve of the objective functional at multiple points to provide a highly accurate estimate for the step length. Such a large number of samples would be computationally unfeasible using conventional methods.
  • the model is updated using the gradient obtained in step 108 and the step length obtained in step 1 10.
  • the model update derives from the gradient i.e. the partial derivative with respect to a point perturbation of the model m at each position. Ultimately gradients from separate shots will be summed when forming the final model update.
  • this is the result of two wavefields: an incident wavefield emitted by a source at the source location, and a back- propagated wavefield which is emitted by a (multi-point) source located the receiver positions.
  • Step 114 Convergence criteria met?
  • step 1 14 it is determined whether convergence criteria have been met. For example, when the method is deemed to have reached convergence when the difference between the data sets reaches a threshold percentage or other value. If the criteria as set out above have been met, then the method proceeds to step 1 16 and finishes with the resultant Earth model generated. If the criteria have not been met, then the method proceeds back to repeat steps 104 to 112 as discussed above.
  • Step 116 Finish
  • the method finishes and the modelled subsurface portion of the Earth is deemed to be sufficiently accurate to be used for subsurface exploration.
  • This may involve the direct interpretation of the recovered model, and/or involve the process of depth-migration to generate a subsurface reflectivity image to be used for the identification of subsurface features such as cavities or channels which may contain natural resources such as hydrocarbons. Examples of such hydrocarbons are oil and natural gas.
  • FIG. 8 A second embodiment of the invention is illustrated in Figure 8. The second embodiment focuses on a specific application of the present invention to an AWI-type method using Wiener filters. Step 200: Obtain observed data set
  • Step 200 corresponds substantially to method step 100 of the previous embodiment. Therefore, this will not be described again here. Only the steps which are new to this embodiment of the method of the present invention will be described. The method now proceeds to step 202.
  • Step 202 Provide starting model
  • 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 model is generated in this step in accordance with step 102 above.
  • the generated model consists of values of the coefficient V p and, possibly, other physical values or coefficients over a discrete grid representing the subsurface.
  • Such starting models are routinely generated and represent the general trends of the major features within the subsurface region to be modelled and could be readily generated by the skilled person.
  • the method then proceeds to step 204.
  • Step 204 Generate predicted data set
  • the predicted data set corresponds discrete point to discrete point to the observed dataset.
  • the predicted data set is generated for the same measurement parameter(s) at the same frequency or frequencies.
  • the predicted data set is generated using the full two-way wave equation.
  • the modelled trace data may be generated using the time domain full two-way wave equation as set out above. From the above analysis, predicted seismic data can be generated for one or more physical parameters and at one or more selected frequencies. This forms the modelled data set d P reci(r,s). The method now proceeds to step 206. Step 206: Scale data
  • the data is scaled so that the predicted d (r,s) and observed d 0 (r,s) data sets comprise matching root mean square amplitudes. This may be accomplished by any known approach. Either or both of the predicted d(r,s) and observed d 0 (r,s) data sets may be scaled as appropriate.
  • Step 208 Design Wiener filter
  • a Wiener filter is a convolutional filter of finite length that is operable to convert an input wavelet into a desired output wavelet in a least-squares manner. In other words, if b is the input wavelet, and c is the desired output wavelet, then the Wiener filter w minimises expression 21 ).
  • B is a rectangular matrix with M columns and M + N -1 rows that contains the elements of b arranged as:
  • B T B represents a matrix containing the autocorrelation of b in each column with the zero lag on the main diagonal
  • B T c is the cross-correlation of b and c.
  • the approach is to find the auto-correlation B T B of the input trace b, the cross- correlation B T c of the input trace with the desired output trace c, and deconvolve the latter using the former.
  • fast algorithms such as the Levinson algorithm
  • equation 25 would be solved:
  • Equation 26) and 27) denote the approach of this embodiment.
  • Step 210 Generate reference filter
  • an analogous convolutional filter is designed such that, when applied to an input dataset, this reference filter will generate an output dataset that provides an equivalent approximation to all or parts of the input dataset. In other words, the reference filter does not modify the input data in essence.
  • Step 212 Construct objective function
  • an objective function is configured.
  • the objective function (or objective function) is configured to measure the dis-similarity between the actual filter coefficients and reference filter coefficients.
  • an objective function may be configured that measures similarity; in this case, step 214 (described later) will operable to maximise rather than minimise the objective function.
  • An example of an objective function configured to measure the dis-similarity between a simple one-dimensional convolutional filter in time and a reference function that consist of only one non-zero value at zero lag, would be to weight the convolution filter coefficients by the modulus of their temporal lag. The objective function would then consist of some norm of these weighted coefficients divided by the same norm of the unweighted coefficients. If the L 2 norm is used here, then this objective function will provide the least-squares solution, but other norms are possible and potentially desirable.
  • the norm of the weighted coefficients must be normalised by the norm of the unweighted coefficients in this example otherwise the objective function could be simply minimised by driving the predicted data to large values, and hence driving the filter coefficients to small values.
  • the coefficients generated for each source receiver pair r,s in step 208 are weighted as a function of the modulus of the temporal lag.
  • the coefficients are weighted based on the data position in time for a time domain analysis.
  • weighting could be used.
  • more complicated functions of the temporal lag such as weighting with a Gaussian function of lag centred on zero lag.
  • weighting with a Gaussian function of lag centred on zero lag.
  • two types of weighting are desirable; those that increase monotonically away from zero lag, such as the modulus, and those that decrease monotonically away from zero lag, such as a Gaussian weighting.
  • the former type of weighting will lead to objective functions that should be minimised and the latter type will lead to objective functions that should be maximised. Combinations of these two types are also possible. The skilled person would readily understand how to design such objective functions and how to minimise or maximise them.
  • a model is then sought that makes the Wiener filters as close as possible to being just a spike at zero time lag.
  • the predicted and modelled data matches apart from a scale factor. Therefore, a model m is desired that minimises or maximises:
  • Step 214 Minimise objective function (determine gradient)
  • the gradient method is used to determine the minimum of the objective function.
  • the gradient is straightforward to derive from equation 28):
  • the gradient can be obtained based on equation 29) as follows.
  • a set of Wiener filters is found in step 208, one filter per data trace.
  • the filter coefficients are normalised by their inner product, trace-by-trace.
  • the normalised coefficients are weighted by a function of temporal lag.
  • the resultant sequence is deconvolved by the auto-correlation of the observed data, and convolved with the observed data. This forms an adjoint source for each source- receiver pair. These adjoint sources are then back-propagated and combined in the normal way with the forward wavefield to produce the gradient.
  • the gradient calculation is performed using a "full" objective function. In other words, the gradient calculation is performed using the entire observed data set, pre-processed derivative thereof, or pre-defined sub-set.
  • Step 216 Minimise objective function (determine step length)
  • the step length a is determined.
  • the gradient identifies in which direction the model should be changed. However, it does not specify by how much.
  • the step length a is computed. Aspects of the method are as described in step 1 10. Subject matter in common will not be repeated here.
  • the use of the present invention is particularly advantageous in the case of AWI. This is because the residuals and the model parameters are not linearly related, and so the second method of step 1 10, i.e. the Born approximation method for calculating step length, is not valid. Additionally, the shape of the objective function in the update direction is a combination of concave and convex curves. Therefore, the method of fitting a parabola is often inaccurate. Consequently, the present invention provides a robust method for determining step length in such cases.
  • step 1 model update m k+1 which will have residuals 5m k .
  • each forward model run to determine the step length is carried using a reduced objective function utilising only a sub-set of the full observed data set.
  • a small sub-set of the full observed data set is used. For example, less than 40% of the data of the observed data set used for calculation of the gradient in the previous step may be used in the determination of step length . In embodiments, less than 15% of the data or shots) of the data of the observed data set used for calculation of the gradient in the previous step may be used in the determination of step length.
  • the sub-set of the observed data set used for determination of step length can be reduced below 10%, or even further.
  • the inventors have found that useful results can be obtained using only a very small sub-set of the observed data set used in the objective function for determining the gradient in step 108. For example, 2-4% of the observed data set used in the calculation of the gradient may be used to determine in the or each calculation of the step length.
  • the sub-set of the observed data set used to determine step length may be differentiated in terms of numbers of "shots". As described above, a shot is a group of observed data traces belonging to the same source or receiver. For the same source, a shot would be a group of observed data trace taken at different receiver locations for the same source. Conversely, for a shot identified by individual receiver, a shot would be a group of observed data trace taken at the same receiver location for different source locations.
  • the size of the data set is, at least in part, determined by the size of the observed data set. This is because, as described above, a predicted data set is generated based on the individual data points (i.e. source and receiver locations) in the observed data set. The differences between the two are then minimised. Therefore, the number of data points and number of sampled parameters in the observed data set defines the overall size of the data set for processing in the objective function. It is noted that different types of data set may be used. For example, composite shot data sets may be used with the method of the present invention.
  • the form and nature of the data set is not material to the present invention, provided that the data set used to determine the step length is reduced in size when compared to the data set used to determine the gradient in step 212. Indeed, the two data sets need not overlap entirely (although in practice this will be necessary). For example, it is to be considered to be within the scope of the present invention for a sub-set of the original observed data set to be used to calculate the gradient in step 212. This may be done either through appropriate selection of data or through processing such as composition of data shots or other precalculation processing. A different sub-set of the original data may then be used for calculation of the step length, provided the absolute size of the data set used is smaller than the data set used for determination of the gradient.
  • Figure 9 shows an example of step length determination for an objective function using the AWI method.
  • An objective function computed using all shots i.e. the full observed data set
  • the objective function was computed using an observed data set comprising 91 shots.
  • a reduced objective function computed using only 3 shots is also shown. Note that the reduced objective function is identical to the objective function save for the use of a reduced number of shots in its calculation.
  • the two graphs vary in form, as may be expected.
  • the reduced objective function calculated using only three shots has the same global minimum in a as the full objective function calculated using 91 shots. This is the key parameter to determine the optimum step length a, and so this example shows that the reduced objective function calculation can provide an accurate determination of a from a reduced sub-set of shots.
  • the calculation of the reduced objective function is significantly less expensive computationally, and so the calculation can be completed with a reduced load on computational resources.
  • the number of sampled points is five or more. However, 50 to 100 samples could be used to sample the curve of the objective functional at multiple points to provide a highly accurate estimate for the step length. Such a large number of samples would be computationally unfeasible using conventional methods.
  • the model is updated using the gradient obtained in step 214 and the step length in step 216.
  • the model update derives from the gradient of equation 1 1 ) i.e. the partial derivative with respect to a point perturbation of the model m at each position. Ultimately gradients from separate shots will summed when forming the final model update.
  • this is the product of two wavefields: an incident wavefield emitted by a source at source location back-propagated wavefield which is emitted by a (multi-point) source located the receiver positions.
  • gradient methods can be enhanced using an approximate form for the Hessian and conjugate directions for the model update. Further a step-length is then calculated to scale the search direction and give the final model update.
  • step 220 The method now proceeds to step 220.
  • Step 220 Convergence criteria met?
  • step 220 it is determined whether convergence criteria have been met. For example, when the method is deemed to have reached convergence when the difference between the data sets reaches a threshold percentage. If the criteria as set out above have been met, then the method proceeds to step 222 and finishes with the resultant Earth model generated. If the criteria have not been met, then the method proceeds back to repeat steps 204 to 218 as discussed above. Step 222: Finish
  • the method finishes and the modelled subsurface portion of the Earth is deemed to be sufficiently accurate to be used for subsurface exploration.
  • This may involve the direct interpretation of the recovered model, and/or involve the process of depth-migration to generate a subsurface reflectivity image to be used for the identification of subsurface features such as cavities or channels which may contain natural resources such as hydrocarbons. Examples of such hydrocarbons are oil and natural gas.
  • 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 effectiveness of the AWI method (with or without the additional modification described in the present invention) is illustrated with reference to Figures 10 to 12.
  • Figure 10 shows a relatively extreme subsurface model including a strong anomaly to resolve.
  • Figure 1 1 shows a model derived from conventional FWI. As shown, the model fails to converge correctly due to cycle skipping. Consequently, the minimised model bears little resemblance to the actual model as set out in Figure 10.
  • Figure 12 shows the model recovered using the method of the present invention. As shown, the model is fully recovered subject only to the limitations of geometry and bandwidth. Clearly, the present invention enables accurate modelling of complex subsurface features where conventional FWI fails.
  • Sources and other sub-sets of the data may be selected, modified and/or combined in various ways, and these used in the inversion in preference to using the original physical sources. Many of these enhancements are discussed and described in "An overview of full-waveform inversion in exploration geophysics", J. Virieux and S. Operto, Geophysics Vol. 74, WCC1- WCC26, doi: 10.1 190/1.3238367.
  • the scalar parameter of pressure as its focus (i.e. P-waves)
  • the present invention has been described in relation to single parameter inversions, the present invention can also be used in multi-parameter inversions. Indeed, the use of a reduced data set for step length calculations could provide improvements where multi- parameter initial data sets are very large.
  • the target is the global minimum in a surface with as many dimensions as different inverted parameters.
  • the number of parameters will be two or three.
  • V p , and V s may be inverted in the case of elastic inversion.
  • V p and e could be inverted for anisotropic inversion.
  • Figure 13 shows a schematic representation of an objective function for a multi-parameter inversion. It is represented as a 3D surface as a function of different V p , and V s step sizes.
  • a limitation when applying the method to multi-parameter inversion is that the density of points mapped along each parameter step size is decimated by a factor equal to the number of different parameters to invert. For example, if an acoustic inversion maps the reduced functional with 49 points, an elastic equivalent could generate a 2D map with a resolution of 7x7 points as shown in Figure 13.
  • the present invention may be an alternative to subspace methods (or other techniques that aim at balancing the different parameter updates at each iteration).
  • the proportion of data selected for either the gradient calculation (i.e. the first proportion of the original data set) or the proportion of the data used for the gradient calculation to be used in the step length calculation may be determined in a number of organised ways.
  • the proportion of data to be used may be done by selecting a fraction of the receivers in the original observed data set.
  • only a fraction of the sources in the observed data set may be selected.
  • different post-processing techniques could be used to remove less significant shots or data elements from the observed data set to reduce the data load.
  • the method may start coarser and focus on the interesting regions in a 'second pass'.
  • the step length calculated in the previous iteration of the model update could be used to inform the next iteration of the model update.
  • the method is applicable to composite shots without modification. Problems arising from cross-talk in composite shot methods are unlikely to be problematic because only the objective function is evaluated. Cross-talk problems arise from the fact that forward- and back-propagated energy from different shots correlate. The use of composite shots is likely to be restricted to acquisition geometries where different shots have receivers in common.
  • multidimensional filters could also be used.
  • the filter for each source-receiver pair is designed only using data from that source-receiver pair.
  • the scheme is also implemented in time - the two sequences that are being matched represent data that varies in time, and the Wiener filter has its coefficients arranged by temporal lag.
  • the predicted data can be calculated from that source not only at the receiver also in a three dimensional volume around it.
  • a one-dimensional observed dataset and a four-dimensional predicted dataset is now available. That is, three space dimensions plus time.
  • a convolutional filter can be designed that takes all or some 1 D, 2D or 3D subset of this 4D predicted dataset as input and generates the 1 D observed dataset as output. This can be done in step 208.
  • this convolutional filter is a multi-dimensional Wiener convolutional filter
  • the corresponding reference filter will be unity at the coefficient that corresponds to zero lag in time and zero lag in all of the space dimensions - that is at the position of the receiver - and it will be zero everywhere else.
  • the method remains unchanged but now a 1 D, 2D, 3D or 4D convolutional filter can be designed, and this can be combined with an appropriate function of all the temporal and spatial lags involved. The associated objective functional can then be minimised.
  • Equivalent schemes may also be generated that match the observed data to the predicted data.
  • These multi-dimensional filters could be extended by applying them to virtual sources and receivers located anywhere within the sub-surface including at the source position.
  • frequency domain analysis could also be used.
  • the different frequency components can be extracted through the use of a Fourier transform. Fourier transforms and associated coefficients are well known in the art. Since only a few frequency components may be required, a discrete Fourier transform may be computed from the inputted seismic trace data. Alternatively, a conventional fast Fourier transform (FFT) approach may be used.
  • FFT fast Fourier transform
  • 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 observed data sets at one or more frequencies may be extracted.
  • a plurality of observed data sets may be provided having frequency components of, for example, 4.5 Hz, 5 Hz and 5.5 Hz. These 3 frequencies can be inverted individually with the output of one becoming the input of the next.
  • the process could be applied using the reverse filter, that using a filter that takes the observed data as input and that generates an approximation to the predicted data as output.
  • the data could be pre-processed in any desired manner the observed and/or field data prior to the formation of the Wiener filter coefficients.
  • the auto and/or cross correlation functions could be pre-processed in order to form the Wiener filter.
  • the Wiener filter coefficients and/or the weighted Wiener coefficients could be processed before forming the quantity to be minimised or maximised.
  • post-processing could be performed on the quantity to be minimised or maximised.
  • the scheme could be applied to data, to models, and/or to filter coefficients that have been transformed into some other domain or domains including but not limited to the single or multi-dimensional Fourier domain, the Laplace domain, the Radon domain, the wavelet domain, the curvelet domain.
  • the observed and/or predicted data could be transformed into a sequence of complex numbers, for example by forming an imaginary part using the Hilbert transform, such that the Wiener coefficients are themselves complex numbers, and optionally combine this with the minimisation of the imaginary portion of the Wiener coefficients. Additional constraints could be added to the Wiener coefficients including but not limited to constraints generated by the observed data and/or the predicted data and/or the model and/or the other Wiener coefficients and/or the accuracy of the convolutional filter.
  • Wiener filters could be used that vary continuously or discontinuously in time, in space, with receiver position, with source position, with source-receiver azimuth, with source-receiver offset, with source-receiver midpoint, and/or with other characteristics of the data, the model, the model updates, the source, the receiver, the filer coefficients.
  • the Wiener filter need not be one dimensional and may encompass multiple dimensions.
  • Non-limiting examples may include: Kalman filters; least mean squares filters; normalised least mean squares filters; and recursive least squares filters.
  • 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 visco-acoustic, elastic, visco-elastic or poro-elastic wave equation.

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

L'invention concerne un procédé d'exploration souterraine, le procédé consistant à générer une représentation géophysique d'une partie du volume de la Terre à partir d'une mesure sismique d'au moins un paramètre physique, le procédé comprenant les étapes consistant : à fournir un ensemble de données observées, l'ensemble de données observées comprenant des valeurs de données dérivées de valeurs de mesure sismique de ladite partie du volume de la Terre ; à générer, à l'aide d'un modèle souterrain d'une partie de la Terre comprenant une pluralité de coefficients de modèle, un ensemble de données modélisées comprenant une pluralité de valeurs de données modélisées ; à mettre à jour le modèle souterrain : en générant une fonction objective pouvant être utilisé pour mesurer la non-correspondance ou la similarité entre l'ensemble de données observées et l'ensemble de données prédites ; à déterminer, à l'aide d'une première proportion du nombre total de valeurs de données de l'ensemble de données observées, le gradient de la fonction objective ; à déterminer, à l'aide d'une fonction objective réduite utilisant une seconde proportion du nombre total de valeurs de données de l'ensemble de données observées, la longueur d'étape pour une mise à jour de modèle souterrain, la seconde proportion comprenant 40 % ou moins du nombre total de valeurs de données de la première proportion de l'ensemble de données observées ; et en mettant à jour le modèle souterrain d'une partie de la Terre à l'aide du gradient déterminé et de la longueur d'étape ; et en fournissant un modèle souterrain mis à jour d'une partie de la Terre pour l'exploration souterraine.
PCT/EP2016/062092 2015-05-29 2016-05-27 Procédé amélioré de modélisation d'inversion Ceased WO2016193180A1 (fr)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
GB1509332.1 2015-05-29
GB1509332.1A GB2538804A (en) 2015-05-29 2015-05-29 Improved method for inversion modelling

Publications (1)

Publication Number Publication Date
WO2016193180A1 true WO2016193180A1 (fr) 2016-12-08

Family

ID=53677477

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2016/062092 Ceased WO2016193180A1 (fr) 2015-05-29 2016-05-27 Procédé amélioré de modélisation d'inversion

Country Status (2)

Country Link
GB (1) GB2538804A (fr)
WO (1) WO2016193180A1 (fr)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR102026063B1 (ko) * 2018-05-14 2019-09-27 서울대학교 산학협력단 반복적 파형 역산을 이용한 강건한 지하 구조 영상화 방법
CN110914718A (zh) * 2017-05-22 2020-03-24 沙特阿拉伯石油公司 计算频域中地震速度反演的与振幅无关的梯度
US10908305B2 (en) 2017-06-08 2021-02-02 Total Sa Method for evaluating a geophysical survey acquisition geometry over a region of interest, related process, system and computer program product
US11397273B2 (en) 2019-12-20 2022-07-26 Saudi Arabian Oil Company Full waveform inversion in the midpoint-offset domain
US11768303B2 (en) 2021-04-22 2023-09-26 Saudi Arabian Oil Company Automatic data enhancement for full waveform inversion in the midpoint-offset domain
CN116992780A (zh) * 2023-09-25 2023-11-03 中南大学 数字化电解槽的温度传感器布置方法
WO2024051834A1 (fr) * 2022-09-09 2024-03-14 中国石油化工股份有限公司 Procédé et dispositif d'inversion de forme d'onde complète et support de stockage

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2509223A (en) * 2013-10-29 2014-06-25 Imp Innovations Ltd Reducing or eliminating cycle-skipping in full waveform inversion

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130311149A1 (en) * 2012-05-17 2013-11-21 Yaxun Tang Tomographically Enhanced Full Wavefield Inversion

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2509223A (en) * 2013-10-29 2014-06-25 Imp Innovations Ltd Reducing or eliminating cycle-skipping in full waveform inversion

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
MIKE WARNER* ET AL: "Adaptive waveform inversion: Theory", SEG TECHNICAL PROGRAM EXPANDED ABSTRACTS 2014, 5 August 2014 (2014-08-05), pages 1089 - 1093, XP055296840, DOI: 10.1190/segam2014-0371.1 *
Q. LIU ET AL: "Seismic imaging: From classical to adjoint tomography", TECTONOPHYSICS., vol. 566-567, 22 July 2012 (2012-07-22), NL, pages 31 - 66, XP055254424, ISSN: 0040-1951, DOI: 10.1016/j.tecto.2012.07.006 *
S F REKER ET AL: "We ELI1 09 Applications of Random Dynamic Shot Decimation in Full Waveform Inversion", SHELL GLOBAL SOLUTIONS INTERNATIONAL B.V.) 76 TH EAGE CONFERENCE & EXHIBITION 2014 AMSTERDAM RAI, 16 June 2014 (2014-06-16), pages 1 - 5, XP055216896, DOI: 10.3997/2214-4609.20141119 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110914718A (zh) * 2017-05-22 2020-03-24 沙特阿拉伯石油公司 计算频域中地震速度反演的与振幅无关的梯度
US11269097B2 (en) 2017-05-22 2022-03-08 Saudi Arabian Oil Company Computing amplitude independent gradient for seismic velocity inversion in a frequency domain
CN110914718B (zh) * 2017-05-22 2022-07-15 沙特阿拉伯石油公司 计算频域中地震速度反演的与振幅无关的梯度
US10908305B2 (en) 2017-06-08 2021-02-02 Total Sa Method for evaluating a geophysical survey acquisition geometry over a region of interest, related process, system and computer program product
KR102026063B1 (ko) * 2018-05-14 2019-09-27 서울대학교 산학협력단 반복적 파형 역산을 이용한 강건한 지하 구조 영상화 방법
US11397273B2 (en) 2019-12-20 2022-07-26 Saudi Arabian Oil Company Full waveform inversion in the midpoint-offset domain
US11768303B2 (en) 2021-04-22 2023-09-26 Saudi Arabian Oil Company Automatic data enhancement for full waveform inversion in the midpoint-offset domain
WO2024051834A1 (fr) * 2022-09-09 2024-03-14 中国石油化工股份有限公司 Procédé et dispositif d'inversion de forme d'onde complète et support de stockage
CN116992780A (zh) * 2023-09-25 2023-11-03 中南大学 数字化电解槽的温度传感器布置方法
CN116992780B (zh) * 2023-09-25 2024-01-02 中南大学 数字化电解槽的温度传感器布置方法

Also Published As

Publication number Publication date
GB201509332D0 (en) 2015-07-15
GB2538804A (en) 2016-11-30

Similar Documents

Publication Publication Date Title
EP3063562B1 (fr) Procédés pour l'exploration souterraine, produit programme d'ordinateur et support d'enregistrement lisible par ordinateur
RU2693495C1 (ru) Полная инверсия волнового поля с компенсацией показателя качества
RU2582480C2 (ru) Кодирование одновременных источников и разделение источников в качестве практического решения по инверсии полного волнового поля
US11029432B2 (en) De-aliased source separation method
WO2016193180A1 (fr) Procédé amélioré de modélisation d'inversion
EP1616202A2 (fr) Imagerie sismique par l'intermediaire de la migration d'ondes utilisant une expansion spatiale de krylov de l'operateur d'exposant de la racine carree
EP2730949A2 (fr) Procédé d'atténuation de bruit d'interférence et appareil
WO2013080128A1 (fr) Séparation de données de source simultanée
US20180164453A1 (en) Method for Improved Geophysical Investigation
WO2013093468A2 (fr) Procédé de contrôle de qualité d'inversion de forme d'onde totale
WO2013093467A1 (fr) Procédé et appareil d'inversion de forme d'onde totale
RU2570827C2 (ru) Гибридный способ для полноволновой инверсии с использованием способа одновременных и последовательных источников
GB2481270A (en) Method of, and apparatus for, full waveform inversion modelling
EP4423539A1 (fr) Séparation de sources sismiques par interpolation et démélange conjoints
EP3938810B1 (fr) Récupération de signaux pendant un démixage et une séparation de sources simultanées
EP3776012B1 (fr) Propagateur de champ d'ondes de milieux orthorhombiques inclinés
US20240159930A1 (en) Method and apparatus for implementing full waveform inversion using angle gathers
EP4088144B1 (fr) Procédé et appareil de réalisation de modélisation efficace de sources sismiques en mouvement à durée étendue
Przebindowska Acoustic full waveform inversion of marine reflection seismic data
WO2024147948A1 (fr) Procédé et appareil d'estimation de l'incertitude d'un modèle de vitesse d'une région souterraine
Sergio-Alberto et al. Optimal coding of blended seismic sources for 2D full waveform inversion in time
EA044564B1 (ru) Построение модели скорости

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 16725539

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 16725539

Country of ref document: EP

Kind code of ref document: A1