CROSS REFERENCE TO RELATED APPLICATION

This is a U.S. National Phase Application under 35 U.S.C. § 371 of International Patent Application No. PCT/162018/000757, filed Jun. 8, 2018. The entire contents of which are hereby incorporated by reference.

FIELD OF THE INVENTION

The present invention concerns a computer-implemented method for generating an image of a subsurface of an area of interest from seismic data.

BACKGROUND

In order to determine the structure of ground, the geological nature of the subsurface or the location of oil and gas reservoirs, seismic data may be used.

Indeed, when a seismic wave propagates through the ground from a seismic source to a seismic receiver, the recorded seismic wave contains information relative to the geophysical properties of the underground. Such properties may include seismic wave velocity, density and anisotropy properties.

In exploration geophysics, an image of the subsurface may be determined by recording seismic waves at the surface with a plurality of seismic sensors such as hydrophones or geophones, using various seismic processing methods.

For example, the time between transmission of the seismic wave and reception of the reflected seismic waves may be used to determine the distance to a geological interface. More advanced techniques, such as reflection velocity tomography, are able to use the kinematic information contained in the seismic data to determine a detailed seismic image of the underground.

Another method known as “full waveform inversion” or “full wavefield inversion” (FWI) uses full waveform information provided by the recorded seismic waves, i.e. both phase and amplitude data, both transmission seismic data and reflection seismic data, both primaries and multiples, contained in the recorded seismic dataset, to build a high resolution image of the underground geological formation.

Due to the high computational cost associated with the method, classical implementations utilize local optimization techniques to estimate optimal model parameters. A variety of local optimization methods are generally used. However, the classical FWI inverse problem is highly non-linear and the FWI method generally generates smooth velocity model updates mainly located in the shallow parts of the model. Therefore, as long as the starting model is not sufficiently close to the true model, the local optimization may not converge toward the right solution and may remain trapped in a wrong local solution.

Other approaches based on a model scale separation may also be used. However, these approaches have similar drawbacks and lead to a high computational cost and a slow convergence of the method.

SUMMARY

One aim of the invention is to provide a fast and efficient method for generating a very accurate image of a subsurface of an area of interest from seismic data.

To this aim, the subject-matter of the invention comprises the following steps:providing seismic wavefields recorded with a plurality of seismic receivers over the area,providing a zero-offset seismic wavefield dataset,determining a seismic velocity parameter model w(x), said seismic velocity parameter model associating for each location x of the area a seismic velocity parameter w, said seismic velocity parameter model comprising an initial model w0(x), a low frequency perturbation term δmb(x) and a high frequency perturbation term δmr(x),determining an optimal seismic velocity parameter model wopt(x) by computing a plurality of iterations, each iteration comprising calculating and optimizing a cost function, said cost function being a measure of discrepancies between the recorded seismic wavefields data and the modeled seismic wavefields data obtained using the seismic velocity parameter model w(x), said cost function being dependent on the zero-offset seismic wavefield and on the low frequency perturbation term δmb(x) as a parameter in the optimization of the cost function, the high frequency perturbation term δmr(x) being related to the velocity parameter model w(x) to keep the provided zero-offset seismic wavefield data invariant with respect to the low frequency perturbation term δmb(x).

The method according to the invention may comprise one or more of the following features, taken into consideration in isolation, or according to any one of any technically feasible combination:a linearized perturbation of cost function with respect to δmbcomprises at least a first linearized perturbation term, advantageously Fb(u0,δmr) simulating the effect of the low frequency perturbation term on seismic data residuals and at least a second linearized perturbation term advantageously Fr(u0)Fzo,r−1(u0)Fzo,b(u0,δmr) compensating the effect of the seismic data perturbation to keep the provided zero-offset seismic wavefield data invariant,the seismic velocity parameter model w(x) comprises a compression waves velocity model v or a squared slowness model u(x), the squared seismic slowness u being related to the seismic velocity v by the relation u=1/v2,the seismic velocity parameter model w(x) comprises a seismic anisotropy term or an elasticity parameter term,the step for determining the optimal seismic velocity parameter model wopt(x) comprises expanding the low frequency perturbation term δmb(x) on basis functions, the basis functions being advantageously chosen among: uniform blocks, splines, polynomial or wavelet functions,the step for determining the optimal seismic velocity parameter model wopt(x) comprises determining associated expansion coefficients of the expansion of the low frequency perturbation term δmb(x) on the basis functions,the method comprises, at each iteration, updating the velocity parameter model and evaluating the high frequency perturbation term δmr(x) from the zero-offset seismic dataset and the updated seismic velocity parameter model using a migration operator, the migration operator being an approximate or an exact iterative or non-iterative migration operator,the migration operator is Fzo,r−1(δGzoob,w0+δmb(x)), where Fzo,r−1is an inverse zero offset modeling operator associated with the high frequency perturbation term δmr(x), δGzoobis an observed zero offset time domain data, w0+δmbis a perturbed smooth model taking into account the low frequency perturbation term,the determining step comprises calculating at least the calculated wavefield and calculating pre-stack forward modeled reflection data on the basis of the migrated zero-offset domain seismic and evaluating the misfit between the calculated wavefield and the recorded data and advantageously a regularization term,the calculation of the pre-stack forward modeled reflection data comprises the evaluation of a zero-offset model using a migration operator, based on the velocity model w and the zero offset time domain observed data,the determination of the optimal seismic velocity parameter model wopt(x) comprises implementing a linear or non-linear local optimization method,the implementation step comprises providing at least the gradient of the cost-function with respect to the low frequency perturbation term or/and at least one sensitivity of computed data derivatives with respect to the low frequency perturbation term, the computed data including a calculated refracted/transmitted wavefield calculated from a seismic velocity model w obtained at a previous iteration and a computed pre-stack reflected wavefield obtained from a Born modeling based on the migrated zero-offset time domain data and the seismic velocity model obtained w from the previous iteration,the implementation step comprises:summing a sensitivity contribution of a conventional Full Wave Inversion and a sensitivity contribution of a Reflected Full Wave Inversion andsubtracting a supplementary term resulting from a chain of operators,the chain of operators including a first operator linearly acting on the low frequency perturbation δmb(x) synthetizing a zero-offset wavefield by modeling linear scattering effects generated by the perturbation δmb(x), a second migration operator transforming the synthetized zero-offset wavefield in a reflectivity image using the velocity model of the current iteration, and a third Born modeling operator generating from the reflectivity image a pre-stack modeled data at the same locations than a pre-stack observed data.the method comprises providing seismic wavefields from several shots of seismic sources, and wherein providing at least the gradient of the cost-function with respect to the low frequency perturbation term comprises generating a gradient expression corresponding to a difference of a first gradient term minus a second gradient term, the first gradient term being representative of a sum of:a full wave inversion gradient, computed based on the initial model, anda reflected waveform inversion gradient, the second gradient term, resulting from at least three operations:for each shot, a pre-stack depth migration of seismic data residuals through the action of a first operator, advantageously Fr*({circumflex over (ω)},g,{circumflex over (x)},s,u0), or a pre-stack depth migration of an adjoint source depending on the misfit function between the calculated wavefields and the recorded data through the action of a first operator, advantageously Fr*({circumflex over (ω)},g,{circumflex over (x)},s,u0),demigration to zero-offset time domain through the action of a second operator, advantageously (Fzo,r−1)*({circumflex over (ω)},r,{circumflex over (x)},r,u0), andconversion and inversion of the obtained zero-offset data to get a gradient update along zero-offset wave paths through the action of a third operator, advantageously Fzo,b*({circumflex over (ω)},r,{circumflex over (x)},r,u0,Fzo,r−1(δGzoob,u0)).the method may comprise obtaining an image of the area of interest by extracting data from the optimal seismic velocity parameter model, advantageously by interpolation;the recorded and modeled seismic wavefields comprise transmitted and reflected wavefields;the method is implemented in a mathematical domain chosen among time, Fourier, Laplace or wavelets.

The invention also relates to a process for exploring and/or exploiting an area of interest comprising:applying the method as defined above to obtain an image of a subsurface of the area of interest,drilling the area of interest based on said obtained image.

The invention also deals with a system for generating an image of a subsurface of an area of interest from seismic data, said system comprising:a data recovery module for providing seismic wavefields recorded with a plurality of seismic receivers over the area,a zero-offset module for providing a zero-offset seismic wavefield dataset,a seismic velocity parameter model module for determining a seismic velocity parameter model w(x), said seismic velocity parameter model associating for each location x of the area a seismic velocity parameter w, said seismic velocity parameter model comprising an initial model w0(x), a low frequency perturbation term δmb(x) and a high frequency perturbation term δmr(x),an optimal seismic velocity parameter model module for determining an optimal seismic velocity parameter model wopt(x) by computing a plurality of iterations, each iteration comprising calculating and optimizing a cost function, said cost function being a measure of discrepancies between the recorded seismic wavefields data and the modeled seismic wavefields data obtained using the seismic velocity parameter model w(x), said cost function being dependent on the zero-offset seismic wavefield and on the low frequency perturbation term δmb(x) as a parameter in the optimization of the cost function, the high frequency perturbation term δmr(x) being related to the velocity parameter model w(x) to keep the provided zero-offset seismic wavefield data invariant with respect to the low frequency perturbation term δmb(x).

The invention also deals with a computer program product comprising instructions which, when the program is executed by a computer, cause the computer to carry out the method described above.

