Method for generating an image of a subsurface of an area of interest from seismic data

Abstract

The invention relates to a computer-implemented method for generating an image of a subsurface of an area of interest from seismic data. The method comprises providing seismic wavefields, providing a zero-offset seismic wavefield dataset, determining a seismic velocity parameter model w(x) comprising an initial model w.sub.0(x), a low frequency perturbation term δm.sub.b(x) and a high frequency perturbation term δm.sub.r(x), determining an optimal seismic velocity parameter model w.sub.opt(x) by computing a plurality of iterations, each iteration comprising calculating and optimizing a cost function, said cost function being dependent on the zero-offset seismic wavefield and on the low frequency perturbation term δm.sub.b(x) as a parameter in the optimization of the cost function, the high frequency perturbation term δm.sub.r(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 δm.sub.b(x).

Claims

1. A computer-implemented method for generating an image of a subsurface of an area of interest from seismic data, the method being carried out in a system for generating an image, the method comprising: 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 w.sub.0(x), a low frequency perturbation term δm.sub.b(x) and a high frequency perturbation term δm.sub.r(x), and determining an optimal seismic velocity parameter model w.sub.opt(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 δm.sub.b(x) as a parameter in the optimization of the cost function, the high frequency perturbation term δm.sub.r(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 δm.sub.b(x), wherein a linearized perturbation of cost function with respect to δm.sub.b comprises at least a first linearized perturbation term (F.sub.b(u.sub.0,δm.sub.r)) simulating the effect of the low frequency perturbation term on seismic data residuals and at least a second linearized perturbation term (F.sub.r(u.sub.0)F.sub.zo,r.sup.−1(u.sub.0) F.sub.zo,b(u.sub.0,δm.sub.r)) compensating the effect of the seismic data perturbation to keep the provided zero-offset seismic wavefield data invariant, the determination of the optimal seismic velocity parameter model w.sub.opt(x) comprising implementing a linear or a non-linear local optimization method, the implementation step comprising providing at least the gradient (J.sub.i.sup.TWIN) of the cost-function with respect to the low frequency perturbation term or/and at least one sensitivity (F.sub.j.sup.TWIN) of computed data derivatives with respect to the low frequency perturbation term, the computed data including a calculated refracted/transmitted wavefield (G.sup.cal(ω,g,s,w)) calculated from a seismic velocity model w obtained at a previous iteration and a computed pre-stack reflected wavefield (δG.sub.TWIN.sup.cal(ω,g,s,w,δG.sub.zo.sup.ob)) 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.

2. The method according to claim 1, wherein 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/v.sup.2.

3. The method according to claim 1, wherein the seismic velocity parameter model w(x) comprises a seismic anisotropy term or an elasticity parameter term.

4. The method according to claim 1, wherein the step for determining the optimal seismic velocity parameter model w.sub.opt(x) comprises expanding the low frequency perturbation term δm.sub.b(x) on basis functions.

5. The method according to claim 4, wherein the step for determining the optimal seismic velocity parameter model w.sub.opt(x) comprises determining associated expansion coefficients of the expansion of the low frequency perturbation term δm.sub.b(x) on the basis functions.

6. The method according to claim 4, wherein the basis functions are chosen among: uniform blocks, splines, polynomial or wavelet functions.

7. The method according to claim 1, comprising, at each iteration, updating the velocity parameter model and evaluating the high frequency perturbation term δm.sub.r(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.

8. The method according to claim 7, wherein the migration operator is F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,w.sub.0+δm.sub.b(x)), where F.sub.zo,r.sup.−1 is an inverse zero offset modeling operator associated with the high frequency perturbation term δm.sub.r(x), δG.sub.zo.sup.ob is an observed zero offset time domain data, w.sub.0+δm.sub.b is a perturbed smooth model taking into account the low frequency perturbation term.

9. The method according to claim 1, wherein the determining step comprises calculating at least the calculated wavefield G.sup.cal(ω,g,s,w) and calculating pre-stack forward modeled reflection data (δG.sub.TWIN.sup.cal(ω,g,s,w,δG.sub.zo.sup.ob)) on the basis of migrated zero-offset domain seismic and evaluating the misfit between the calculated wavefield and the recorded data (G.sup.ob).

10. The method according to claim 9, wherein the calculation of the pre-stack forward modeled reflection data (δG.sub.TWIN.sup.cal(ω,g,s,w,δG.sub.zo.sup.ob)) 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 (δG.sub.zo.sup.ob).

11. The method according to claim 9, comprising evaluating the misfit between the calculated wavefield, the recorded data (G.sup.ob) and a regularization term.

12. The method according to claim 1, wherein the implementation step comprises: summing a sensitivity contribution of a Full Wave Inversion and a sensitivity contribution of a Reflected Full Wave Inversion and subtracting a supplementary term resulting from a chain of operators, the chain of operators including a first operator linearly acting on the low frequency perturbation δm.sub.b(x) synthetizing a zero-offset wavefield by modeling linear scattering effects generated by the perturbation δm.sub.b(x), a second migration operator (F.sub.zo,r.sup.−1) 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.

13. The method according to claim 1, comprising providing seismic wavefields from several shots of seismic sources, and wherein providing at least the gradient (J.sub.i.sup.TWIN) 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 (J.sub.b), minus a second gradient term (J.sub.b.sup.zo), the first gradient term J.sub.b being representative of a sum of: a full wave inversion (FWI) gradient, computed based on the initial model (u.sub.0), and a reflected waveform inversion (RFWI) gradient, the second gradient term (J.sub.b.sup.zo), resulting from at least three operations: for each shot, a pre-stack depth migration of seismic data residuals (δR.sup.ob(ω,g,s,u.sub.0)) through the action of a first operator (F.sub.r*({circumflex over (ω)},g,{circumflex over (x)},s,u.sub.0)), 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 (F.sub.r*({circumflex over (ω)},g,{circumflex over (x)},s,u.sub.0)), demigration to zero-offset time domain through the action of a second operator ((F.sub.zo,r.sup.−1)*({circumflex over (ω)},r,{circumflex over (x)},r,u.sub.0)) and conversion 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 (F.sub.zo,b*({circumflex over (ω)},r,x,r,u.sub.0,F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0))).

14. A computer program product comprising instructions which, when the program is executed by a computer, cause the computer to carry out the method according to claim 1.

15. 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 w.sub.0(x), a low frequency perturbation term δm.sub.b(x) and a high frequency perturbation term δm.sub.r(x), an optimal seismic velocity parameter model module for determining an optimal seismic velocity parameter model w.sub.opt(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 δm.sub.b(x) as a parameter in the optimization of the cost function, the high frequency perturbation term δm.sub.r(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 δm.sub.b(x), a linearized perturbation of cost function with respect to δm.sub.b comprises at least a first linearized perturbation term (F.sub.b(u.sub.0,δm.sub.r)) simulating the effect of the low frequency perturbation term on seismic data residuals and at least a second linearized perturbation term (F.sub.r(u.sub.0)F.sub.zo,r.sup.−1(u.sub.0)F.sub.zo,b(u.sub.0,δm.sub.r)) compensating the effect of the seismic data perturbation to keep the provided zero-offset seismic wavefield data invariant, the optimal seismic velocity parameter model module being further configured to implement a linear or a non-linear local optimization method, by providing at least the gradient (J.sub.i.sup.TWIN) of the cost-function with respect to the low frequency perturbation term or/and at least one sensitivity (F.sub.j.sup.TWIN) of computed data derivatives with respect to the low frequency perturbation term, the computed data including a calculated refracted/transmitted wavefield (G.sup.cal(ω,g,s,w)) calculated from a seismic velocity model w obtained at a previous iteration and a computed pre-stack reflected wavefield (δG.sub.TWIN.sup.cal(ω,g,s,w,δG.sub.zo.sup.ob)) 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.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The invention will be better understood, upon reading of the following description, taken solely as an example, and made in reference to the following drawings, in which:

(2) FIG. 1 is a schematic cross-section of a subsurface of an area where a seismic survey acquisition is carried out,

(3) FIG. 2 is a schematic representation of system for generating an image of a subsurface of an area from seismic data, according to the invention,

(4) FIG. 3 is a flow chart of a method for generating an image of a subsurface of an area from seismic data, according to the invention,

(5) FIG. 4 is a schematic representation of a zero-offset seismic data,

(6) FIG. 5 is a schematic representation of the contributions of the different terms of the gradient of the modified cost function obtained with the method according to the invention,

(7) FIG. 6 is a graphical view of the convergence rate of a standard Reflection Full Waveform Inversion (RFWI) and the convergence rate of the method according to the invention, applied on the same 2D seismic dataset, and

(8) FIG. 7 and FIG. 8 are the inversion results obtained respectively with the standard Reflection Full Waveform Inversion (RFWI) and the method according to the invention, for a various number of iterations.

DETAILED DESCRIPTION OF THE DRAWINGS

(9) In the following description, the quantities in bold, such as x refer to vectors quantities.

(10) FIG. 1 show a schematic cross-section of a subsurface 10 of an area of interest 12, such as an oil and gas field, where a seismic survey acquisition is carried out using a seismic survey apparatus.

(11) The seismic survey acquisition is typically a land seismic acquisition or a marine seismic acquisition.

(12) The subsurface comprises a plurality of geological layers 14. Each geological layer 14 is characterized by seismic properties such as seismic velocities, density, anisotropy, etc.

(13) The seismic acquisition comprises generating seismic waves into the ground with seismic sources 16 of the seismic survey apparatus, such as explosive, and recording seismic wavefields at the surface with a plurality of seismic sensors 18 of the seismic survey apparatus, such as hydrophones or geophones.

(14) Typically, seismic wavefields from several shots of seismic sources are recorded.

(15) FIG. 1 shows examples of wavefields recorded by the seismic sensors resulting from primary reflections 20 on the interfaces of the geological layers, from transmissions 22 through the geological layer or from multiple reflections 24 in a geological layer 14.

(16) FIG. 2 shows a system 100 for generating an image of a subsurface of an area from seismic data, according to the invention.

(17) The system 100 comprises a calculator 110 for generating an image of a subsurface of an area from seismic data, a display unit 120 connected to the calculator 110 to display the results provided by the calculator and a man-machine interface 130.

(18) The calculator 110 comprises a database 140. The database 140 contains seismic wavefields recorded with a plurality of seismic sensors in the area of interest.

(19) The database 140 is moreover able to store the results provided by the calculator 110.

(20) In the example of FIG. 2, the database 140 is a local database comprised in the calculator 110. In a variant (not represented), the database 140 is a remote database connected to the calculator 110 by a network.

(21) The calculator 110 comprises a processor 150 and a memory 160 receiving software modules. The processor 150 is able to execute the software modules received in the memory 160 to carry out the method according to the invention.

(22) The memory contains a data recovery 170 module for providing seismic wavefields recorded with a plurality of seismic sensors over the area. Typically, the seismic wavefields are provided by loading the data from the database 140.

(23) The memory 160 comprises a zero-offset module 180 for providing a zero-offset seismic wavefield dataset.

(24) The memory 160 also comprises a seismic velocity parameter model module 190 for determining a seismic velocity parameter model w(x), said seismic velocity parameter model associating for each location x of the area 12 a seismic velocity parameter w, said seismic velocity parameter model comprising an initial model w.sub.0(x), a low frequency perturbation term δm.sub.b(x) and a high frequency perturbation term δm.sub.r(x). 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/v.sup.2.

(25) In the following description, the seismic velocity parameter model w(x) used as a reference is a squared slowness model u(x).

(26) In a variant, the seismic velocity parameter model w(x) comprises a seismic anisotropy term or an elasticity parameter term.

(27) Anisotropy terms are defined as terms accounting for the velocity of the waves in the different directions of the space. They are defined as combinations of the elements of the stiffness tensors. For instance, widely used “Thomsen” “weak” anisotropy parameters ε, δ, γ are defined in “Weak Elastic Anisotropy”, Thomsen, 1986, Geophysics, 51 (10): 1954-66.

(28) The memory further comprises an optimal seismic velocity parameter model module 200 for determining an optimal seismic velocity parameter model w.sub.opt(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 δm.sub.b(x) as a parameter in the optimization of the cost function, the high frequency perturbation term δm.sub.r(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 δm.sub.b(x).

(29) The cost function may include a classical regularization term for taking into account a priori information on the velocity parameter model. Examples of regularization terms and implementation may be found, for example, in “Parameter Estimation and Inverse Problems”, Aster et al., 2011, Volume 90 International Geophysics, Academic Press, or in “Inverse Problem Theory and Methods for Model Parameter Estimation”, Tarantola, 2005, SIAM.

(30) In a variant, the optimal seismic velocity parameter model module 200 comprises an optimization submodule 210 for computing the gradient and/or the Fréchet derivatives and performing a linear or nonlinear local optimization method.

(31) The display unit 120 is for example able to display the seismic images obtained with the method according to the invention. Moreover, the display unit 120 may display information relative to the progress of the method, such as the number of iterations, the lapsed time, etc.

(32) Typically, the display unit 120 is a standard computer screen.

(33) The man-machine interface 130 typically comprises a keyboard, a mouse and/or a touch screen to allow the user to activate the calculator 110 and the various software modules 170, 180, 190, 200 contained in the memory 160 to be processed by the processor 150.

(34) A flow chart of a method for generating an image of a subsurface of an area from seismic data, according to the invention, carried out with a system 100 as described above is shown in FIG. 3.

(35) Advantageously, the method may be implemented for generating an image of the subsurface of the area of interest 12 in two dimensions or in three dimensions.

(36) The method comprises a step 500 for providing seismic wavefields recorded with a plurality of seismic sensors over the area, such as the area presented in FIG. 1.

(37) The physical equation which governs the propagation of the seismic waves in the subsurface of the area is the wave equation.

(38) The scalar wave equation can be written in the frequency domain as follows:
[∇.sup.2+ω.sup.2u(x)]G(ω,x,s)=−δ(x−s)  (1)

(39) u(x) stands for the squared slowness at location x and it is related to heterogeneous wave propagation velocity model v(x) by the relation u(x)=1/v.sup.2(x).

(40) For the sake of clarity, in the following, u(x) is directly referred as the velocity model.

(41) G(ω,x,s) is the Green function solution of the wave equation (1) for temporal frequency ω, observed at position x, i.e. the position of the seismic receiver 18, and for a seismic source 16 at a spatial position s.

(42) ∇.sup.2 and δ are respectively the Laplacian operator and the Dirac function.

(43) A squared slowness model perturbation δu(x) to a reference or background model u.sub.0(x) is considered following the equation:
u(x)=u.sub.0(x)+δu(x)  (2)

(44) For the sake of clarity, in the following, δu(x) is directly referred as the velocity model perturbation.

(45) Similarly, an associated Green function perturbation associated to a position g of the seismic receiver 18, following the below mentioned relation is considered:
G(ω,g,s)=G.sub.0(ω,g,s)+δG(ω,g,s)  (3)

(46) G.sub.0(ω, x, s) is the solution of the wave equation for the reference velocity model u.sub.0(x).

(47) δG(ω,x, s) can be expressed as a solution of a non-linear integral equation, that takes the form:
δG(ω,g,s,u.sub.0,δu)=∫dx ω.sup.2G.sub.0(ω,g,x,u.sub.0)G(ω,x,s,u.sub.0+δu)δu(x)  (4)

(48) Equation (4) is known as the Lipmann-Schwinger equation. It is the exact integral equivalent form for the partial differential wave equation (1).

(49) Recursive expansion of equation (3) in the right-hand term of equation (4) leads to infinite scattering series. Truncating the series to the first two terms, it is obtained:

(50) $\begin{array}{cc}\delta G\left(\omega ,g,s,u_0,\delta u\right)=\int \mathrm{dx}{\omega }^2{G}_0\left(\omega ,g,x,u_0\right)G\left(\omega ,x,s,u_0\right)\delta u\left(x\right)+\int \mathrm{dx}\int d\hat{x}{\omega }^4{G}_0\left(\omega ,g,x,u_0\right)G_0\left(\omega ,x,\hat{x},u_0\right)G_0\left(\omega ,\hat{x},s,u_0\right)\delta u\left(x\right)\delta u\left(\hat{x}\right)+O\left(\delta u^3\right)& \left(5\right)\end{array}$

(51) where the vector variables x and {circumflex over (x)} stand for location of the perturbations.

(52) In a more compact form, it can be written:

(53) $\begin{array}{cc}\delta G\left(\omega ,g,s,u_0,\delta u\right)=\int {\mathrm{dxB}}_1\left(\omega ,g,x,s,u_0\right)\delta u\left(x\right)+\int \mathrm{dx}\int d\hat{x}{B}_2\left(\omega ,g,x,\hat{x},s,u_0\right)\delta u\left(x\right)\delta u\left(\hat{x}\right)+O\left(\delta u^3\right)& \left(6\right)\end{array}$

(54) This form makes clearer the contribution of the first and second order model perturbations and the associated first and second order Fréchet derivatives B.sub.1 and B.sub.2.

(55) B.sub.1 and B.sub.2 have the following forms:

(56) $\begin{array}{cc}\begin{array}{c}{B}_1\left(\omega ,g,x,s,u_0\right)=\frac{\partial G\left(\omega ,g,s,u_0,\delta u\right)}{\partial u\left(x\right)}\\ ={\omega }^2{G}_0\left(\omega ,g,x,u_0\right)G\left(\omega ,x,s,u_0\right),\end{array}& \left(7\right)\\ \begin{array}{c}{B}_2\left(\omega ,g,x,\hat{x},s,u_0\right)=\frac{{\partial }^{2}G\left(\omega ,g,s,u_0,\delta u\right)}{\partial u\left(x\right)\delta u\left(\hat{x}\right)}\\ ={\omega }^4{G}_0\left(\omega ,g,x,u_0\right)G_0\left(\omega ,x,\hat{x},u_0\right)G_0\left(\omega ,\hat{x},s,u_0\right)\end{array}& \left(8\right)\end{array}$

(57) B.sub.1 is classically known as the kernel of the 1.sup.st Born approximation.

(58) The method comprises modeling a scale separated model perturbation.

(59) To this aim, the velocity model perturbation δu(x) is considered as a sum of a smooth or low frequency perturbation δm.sub.b and a short scale or high frequency perturbation δm.sub.r.

(60) The subscripts b and r respectively refer to the background model perturbation and to the reflectivity model perturbation.

(61) The above-mentioned decomposition can be written as the following:

(62) $\begin{array}{cc}u\left(x\right)=u_0\left(x\right)+\delta u\left(x\right)& \mathrm{ }\left(9\right)\\ =u_0\left(x\right)+\delta m_r\left(x\right)+\delta m_b\left(x\right)& \left(10\right)\end{array}$

(63) Equation (10) is introduced into the second order Born series equation (6):

(64) $\begin{array}{cc}\delta G\left(\omega ,g,s,\delta m_r,\delta m_b\right)=\int dx{B}_1\left(\omega ,g,x,s,u_0\right)\delta m_r\left(x\right)+\int dx{B}_1\left(\omega ,g,x,s,u_0\right)\delta m_b\left(x\right)+\int dx\int d\hat{x}{B}_2\left(\omega ,g,x,\hat{x},s,u_0\right)\delta m_r\left(x\right)\delta m_b\left(\hat{x}\right)+\int dx\int d\hat{x}{B}_2\left(\omega ,g,x,\hat{x},s,u_0\right)\delta m_b\left(x\right)\delta m_r\left(\hat{x}\right)+\int dx\int d\hat{x}{B}_2\left(\omega ,g,x,\hat{x},s,u_0\right)\delta m_r\left(x\right)\delta m_r\left(\hat{x}\right)+\int dx\int d\hat{x}{B}_2\left(\omega ,g,x,\hat{x},s,u_0\right)\delta m_b\left(x\right)\delta m_b\left(\hat{x}\right)+O\left(\delta u^3\right)& \left(11\right)\end{array}$

(65) The first order effect of both high frequency perturbation δm.sub.r, and low frequency perturbation δm.sub.b are only considered.

(66) Then the last two terms of equation (11) are neglected, and equation (11) can be written as:

(67) $\begin{array}{cc}\delta G\left(\omega ,g,s,u_0,\delta m_r,\delta m_b\right)=\int {\mathrm{dxB}}_1\left(\omega ,g,x,s,u_0\right)\delta m_r\left(x\right)+\int {\mathrm{dxB}}_1\left(\omega ,g,x,s,u_0\right)\delta m_b\left(x\right)+\int \mathrm{dx}\int d\hat{x}\left\{B_2\left(\omega ,g,\hat{x},x,s,u_0\right)+B_2\left(\omega ,g,x,\hat{x},s,u_0\right)\right\}\delta m_r\left(\hat{x}\right)\delta m_b\left(x\right)+O\left({\delta }^2{m}_r,{\delta }^{2m_b}\right)& \left(12\right)\end{array}$

(68) The Frechet derivatives F.sub.r and F.sub.b associated respectively to δm.sub.r and δm.sub.b can be expressed as the following:

(69) $\begin{array}{cc}\begin{array}{c}{F}_r\left(\omega ,g,x,s,u_0\right)=\frac{\partial G\left(\omega ,g,s,u_0,\delta m_r,\delta m_b\right)}{\partial \delta m_r\left(x\right)}\\ =B_1\left(\omega ,g,x,s,u_0\right)\\ ={\omega }^2{G}_0\left(\omega ,g,x,u_0\right)G\left(\omega ,x,s,u_0\right)\end{array}\mathrm{and}& \left(13\right)\\ F_b\left(\omega ,g,x,s,\delta m_r\right)=\frac{\partial G\left(\omega ,g,s,u_0,\delta m_r,\delta m_b\right)}{\partial \delta m_b\left(x\right)}=B_1\left(\omega ,g,x,s,u_0\right)+\int d\hat{x}\left\{B_2\left(\omega ,g,\hat{x},x,s,u_0\right)+B_2\left(\omega ,g,x,\hat{x},s,u_0\right)\right\}\delta m_r\left(\hat{x}\right)={\omega }^2{G}_0\left(\omega ,g,x,u_0\right)G\left(\omega ,x,s,u_0\right)+{\omega }^4\int d\hat{x}{G}_0\left(\omega ,g,\hat{x},u_0\right)G_0\left(\omega ,\hat{x},x,u_0\right)G_0\left(\omega ,x,s,u_0\right)\delta m_r\left(\hat{x}\right)+{\omega }^4\int d\hat{x}{G}_0\left(\omega ,g,x,u_0\right)G_0\left(\omega ,x,\hat{x},u_0\right)G_0\left(\omega ,\hat{x},s,u_0\right)\delta m_r\left(\hat{x}\right)& \left(14\right)\end{array}$

(70) By introducing equations (13) and (14) in equation (12), the action of model perturbation on the wavefield can be written in the following compact form as:

(71) $\begin{array}{cc}\delta G\left(\omega ,g,s,u_0\delta m_r,\delta m_b\right)=\int {\mathrm{dxF}}_r\left(\omega ,g,x,s,u_0\right)\delta m_r\left(x\right)+\int {\mathrm{dxF}}_b\left(\omega ,g,x,s,u_0,\delta m_r\right)\delta m_b\left(x\right)& \left(15\right)\end{array}$

(72) The first term on the right hand side of equation (15) simulates the effect of the depth reflectivity or high frequency model perturbation keeping the smooth velocity model u.sub.0 fixed. In other words, it generates the reflected wavefield components. It corresponds to the classical first order Born modeling.

(73) The second term on the right hand side of equation (15) simulates the effect on the wavefield of a smooth velocity model perturbation keeping the depth reflectivity model fixed. It has a main effect of kinematic perturbation to both reflected wavefield components and to transmitted wavefield components.

(74) Previous equation (15) shows the high frequency model perturbation δm.sub.r as an argument of the Frechet derivative F.sub.b relative to the smooth perturbation δm.sub.b. This emphasizes clearly the explicit dependency of δm.sub.r in equation (14).

(75) If the depth reflectivity is set to zero, the second term of the right-hand side of equation (15) simplifies to the action of the kernel B.sub.1 and should generate perturbation only to direct arrival and diving wavefield (i.e. background wavefield).

(76) Compared to the linearized modeling given in the equation (6), equation (15) which represents the linearized modeling after scale separation splits explicitly the action of high frequency and smooth model components. It also designates the terms responsible respectively for reflected events generation (mainly amplitude effect), transmitted events kinematic perturbation and finally reflected events kinematic perturbation.

(77) The usual full wave inversion (FWI) method formulates the problem of retrieving the velocity model using the standard least-squares cost function C.sup.FWI:

(78) $\begin{array}{cc}{C}^{\mathrm{FWI}}\left(u\right)=\frac{1}{2}\int \mathrm{ds}\int \mathrm{dg}\int d\omega {.Math.G^{\mathrm{ob}}\left(\omega ,g,s\right)-G^{\mathrm{cal}}\left(\omega ,g,s,u\right).Math.}^2& \left(16\right)\end{array}$

(79) G.sup.ob and G.sup.cal are respectively the observed and simulated wavefields.

(80) It is assumed that the optimal model u occurs at the minimum value of the cost function C.sup.FWI which corresponds to the misfit or the discrepancies between the recorded seismic wavefields data and the modeled seismic wavefields data obtained using the seismic velocity parameter model w(x)

(81) Standard FWI solves this problem by iterative linearization around a background or reference model u.sub.0. Then, by introducing equation (2), it is obtained:

(82) 0$\begin{array}{cc}{C}^{\mathrm{FWI}}\left(u_0+\delta u\right)=\frac{1}{2}\int ds\int dg\int d\omega {.Math.\delta G^{ob}\left(\omega ,g,s,u_0\right)-\delta G^{\mathrm{ca}l}\left(\omega ,g,s,u_0,\delta u\right).Math.}^2& \left(17\right)\end{array}$

(83) wherein the notations of the observed and simulated wavefield residuals are introduced:
δG.sup.ob(ω,g,s,u.sub.0)=G.sup.ob(ω,g,s)−G.sup.cal(ω,g,s,u.sub.0)  (18)
δG.sup.cal(ω,g,s,u.sub.0,δu)=G.sup.cal(ω,g,s,u.sub.0+δu)−G.sup.cal(ω,g,s,u.sub.0)  (19)

(84) Assuming the linearization in equation (6) for a given linearized iteration around the velocity model u.sub.0 restricted to the first order (i.e. conventional Born approximation), the cost function in equation (17), can be approximated by

(85) $\begin{array}{cc}{C}^{\mathrm{FWI}}\left(\delta u+u_0\right)=\frac{1}{2}\int ds\int dg\int d\omega {.Math.\delta G^{ob}\left(\omega ,g,s,u_0\right)-\int dx{F}_{\delta u}\left(\omega ,g,x,s,u_0\right)\delta u\left(x\right).Math.}^2& \left(20\right)\end{array}$

(86) The notation of the kernel F.sub.δu=B.sub.1 is introduced to emphasize the meaning of Frechet derivative versus the unknown full model perturbation δu without any scale separation.

(87) The method according to the invention comprises introducing the scale separation of model perturbation (equation 2) in the linearized equation (15) of the wavefield perturbation δG.sup.cal in the FWI cost function C.sup.FWI.

(88) This leads to the standard Reflection Full Waveform Inversion (RFWI) cost function denoted C.sup.RFWI that takes the form:

(89) $\begin{array}{cc}{C}^{\mathrm{RFWI}}\left(\delta m_r,\delta m_b\right)=\frac{1}{2}\int ds\int dg\int d\omega .Math.\delta G^{ob}\left(\omega ,g,s,u_0\right)-\int dx{F}_r\left(\omega ,g,x,s,u_0\right){\delta m}_r(x)-\hspace{1em}{\hspace{1em}\int {\mathrm{dxF}}_b\left(\omega ,g,x,s,u_0,\delta m_r\right)\delta m_b\left(x\right).Math.}^2& \left(21\right)\end{array}$

(90) Equation (21) shows that assuming an estimated background velocity u.sub.0, the RFWI approach needs to estimate the two-coupled unknown parameters δm.sub.r and δm.sub.b.

(91) The method then comprises a step 510 for providing a zero-offset seismic wavefield dataset.

(92) FIG. 4 shows an example of zero-offset data 26 which correspond to data acquired with no horizontal distance between the source 16 and the receiver 18, compared to data 28 acquired with a horizontal distance between the source 16 and the receiver 18.

(93) In a variant, the zero-offset seismic wavefield dataset is obtained by time domain processing of the seismic wavefields data or by operating a zero-offset de-migration of an optimized seismic reflectivity image.

(94) Applying a move-out correction and the stacking seismic data acquired with separated sources 16 and receivers 18 gives the data the appearance of zero-offset data.

(95) Typically, zero-offset seismic data are available from direct measurements or from a reliable reconstructed version from seismic processing known from the skilled person, for example, normal moveout (NMO) stack, zero-offset common reflection surface (ZO-CRS), demigration of the pre-stack time migration (PSTM), pre-stack depth migration (PSDM), etc.

(96) In the following, the zero-offset time domain data is identified to the observed wavefield δG.sup.ob(ω,g=s,s).

(97) In the following, the observed and computed zero-offset time domain data at the coincident receiver-source position r are denoted by δG.sub.zo.sup.ob(ω,r) and δG.sub.zo.sup.cal(ω,r,u.sub.0,δm.sub.r,δm.sub.b), respectively.

(98) Equation (15) is now written for δG.sub.zo.sup.cal/(ω,r), the linearized calculated zero-offset time domain data:

(99) $\begin{array}{cc}\delta G_{zo}^{\mathrm{ca}l}\left(\omega ,r,u_0\delta m_r,\delta m_b\right)=\int dx{F}_{zo,r}\left(\omega ,r,x,r,u_0\right)\delta m_r\left(x\right)+\int dx{F}_{zo,b}\left(\omega ,r,x,r,u_0,\delta m_r\right)\delta m_b\left(x\right)& \left(22\right)\end{array}$

(100) The subscript zo for “zero-offset” is used here to avoid confusion with finite offset Fréchet derivative kernel.

(101) The method further comprises providing the inverse of zero-offset reflectivity modeling operator δm.sub.r(x)=F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0) generating a reflectivity model resulting from zero-offset data δG.sub.zo.sup.ob(ω,r) in a velocity model u.sub.0(x). The differential kernel associated to the operator δF.sub.zo,r.sup.−1 is denoted F.sub.zo,r.sup.−1.

(102) The operator F.sub.zo,r.sup.−1 may be approximated and provided by a true amplitude or least squares zero-offset migration implementation or the combination of both.

(103) The high frequency perturbation δm.sub.r, can be inverted by matching the simulated δG.sub.zo.sup.cal to δG.sub.zo.sup.ob, considering a given velocity model.

(104) It follows from this matching that the zero-offset data are independent of the background model u.sub.0+δm.sub.b, i.e. the computed and observed data zero-offset data are an invariant with respect to the velocity perturbation δm.sub.b.

(105) Then the reflectivity image F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0+δm.sub.b) consistent at the first order with the model u.sub.0+δm.sub.b and the zero-offset data δG.sub.zo.sup.ob is:
F.sub.zo,r.sup.−1(δG.sub.zo,r.sup.ob,u.sub.0+δm.sub.b)=∫dr∫dωF.sub.z0,r.sup.−1(ω,r,x,r,u.sub.0)  (23)
(δG.sub.z0.sup.ob(ω,r)−∫dyF.sub.z0,b(ω,r,y,r,u.sub.0,F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0))δm.sub.b(y))

(106) where y indicates the location of the low frequency perturbation δm.sub.b(y)

(107) Equation (23) is an explicit linearized expression which allows an estimation of the high frequency model perturbation or reflectivity model F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob(ω,r), u.sub.0+δm.sub.b) associated to the perturbed smooth model u.sub.0+δm.sub.b. This estimation is based on two input data.

(108) The first input is an initial reflectivity depth image that is kinematically consistent with the unperturbed u.sub.0 (i.e. produced by zero-offset migration using u.sub.0 as velocity model).

(109) The second input data is the zero-offset time domain data δG.sub.z0.sup.ob.

(110) The above-mentioned procedure can be viewed as a linearized zero-offset residual migration operator that takes into account the smooth model perturbation δm.sub.b effect and maps the initial zero-offset depth image δm.sub.r(u.sub.0)=F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0) to another zero-offset depth image δm.sub.r(u.sub.0+δm.sub.b)=F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob, u.sub.0+δm.sub.b).

(111) The method then comprises expanding the low frequency perturbation term δ.sub.m.sub.b on a family of basis functions β.sub.i.

(112) For example, basis functions are splines (A Practical Guide to Splines, Carl de Boor, Springer New York, 2001), uniform blocks or wavelet functions (A Wavelet Tour of Signal Processing: The Sparse Way, Stephane Mallat, third edition, Academic Press, 2008).

(113) It can be written:

(114) $\begin{array}{cc}\delta m_b\left(x\right)=\underset{i}{.Math.}{\beta }_i\left(x\right)\delta {\lambda }_i& \left(24\right)\end{array}$

(115) The set δΛ={δλ.sub.i} are the model expansion coefficients in the orthogonal basis function and then the coefficients δλ.sub.i become the unknown of the model building inverse problem.

(116) Using equations (23) and (24), the action of Frechet derivative F.sub.r(ω,g,x,s,u.sub.0) on the high frequency model perturbation F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0+δm.sub.b) is expressed by:

(117) $\begin{array}{cc}\int dx{F}_r\left(\omega ,g,x,s,u_0\right)F_{zo,r}^{-1}\left(\delta G_{zo}^{ob},u_0+\delta m_b\right)=\int \mathrm{dx}{F}_r\left(\omega ,g,x,s,u_0\right)\int dr\int d\hat{\omega }{F}_{zo,r}^{-1}\left(\hat{\omega },r,x,r,u_0\right)\delta G_{zo}^{ob}\left(\hat{\omega },r\right)-\underset{i}{.Math.}\delta {\lambda }_i\int dx{F}_r\left(\omega ,g,x,s,u_0\right)\int dr\int d\hat{\omega }{F}_{zo,r}^{-1}\left(\hat{\omega },r,x,r,u_0\right).& \left(25\right)\\ \int \mathrm{dy}{F}_{zo,b}\left(\hat{\omega },r,y,r,u_0,F_{zo,r}^{-1}\left(\delta G_{zo}^{ob},u_0\right)\right){\beta }_i\left(y\right)& \end{array}$

(118) where {circumflex over (ω)} indicates the frequency pulsation.

(119) A linear modeling operator F.sub.TWIN is introduced.

(120) $\begin{array}{cc}{F_{\mathrm{TWIN}}\left(\omega ,g,s,u_0\right)}_i=\int dx{F}_b\left(\omega ,g,x,s,u_0,\delta m_r\right){\beta }_i\left(x\right)-\int \mathrm{dx}{F}_r\left(\omega ,g,x,s,u_0\right)\int dr\int d\hat{\omega }{F}_{zo,r}^{-1}\left(\hat{\omega },r,x,r,u_0\right)\int \mathrm{dy}{F}_{zo,b}\left(\hat{\omega },r,y,r,u_0,\delta m_r\right){\beta }_i\left(y\right)& \left(26\right)\end{array}$

(121) This linear modeling operator simulates (up to first order) the effect of the low frequency background model perturbation on seismic data residuals δR.sup.ob which will be defined below.

(122) The first term F.sub.b(u.sub.0,δm.sub.r)=∫dx F.sub.b(ω,g,x,s,u.sub.0,δm.sub.r)β.sub.i(x) of the linear modelling operator F.sub.TWIN(ω,g,s,u.sub.0).sub.i models data perturbation for both transmitted waves components and reflected wavefields. The second term
F.sub.r(u.sub.0)F.sub.zo,r.sup.−1(u.sub.0)F.sub.zo,b(u.sub.0,δm.sub.r)=∫dxF.sub.r(ω,g,x,s,u.sub.0)∫dr∫d{circumflex over (ω)}F.sub.zo,r.sup.−1({circumflex over (ω)},r,x,r,u.sub.0)∫dyF.sub.zo,b({circumflex over (ω)},r,y,r,u.sub.0,δm.sub.r)β.sub.i(y)

(123) of the linear modelling operator F.sub.TWIN(ω,g,s,u.sub.0).sub.i compensates the partial effect of data perturbation to keep the zero offset data subset invariant.

(124) F.sub.TWIN(ω,g,s,u.sub.0).sub.i will be recognized later as the Frechet derivative of the modified cost function used by the method according to the invention.

(125) $\begin{array}{cc}\delta R^{\mathrm{ob}}\left(\omega ,g,s,u_0\right)=\delta G^{\mathrm{ob}}\left(\omega ,g,s,u_0\right)-\delta G_{\mathrm{TWIN}}^{cal}\left(\omega ,g,s,u_0\right)=G^{\mathrm{ob}}\left(\omega ,g,s\right)-\left\{G^{\mathrm{ca}l}\left(\omega ,g,s,u_0\right)+\delta G_{\mathrm{TWIN}}^{cal}\left(\omega ,g,s,u_0,\delta G_{zo}^{ob}\right)\right\}& \left(27\right)\end{array}$

(126) δR.sup.ob corresponds to the modified data residuals to be minimized.
δG.sub.TWIN.sup.cal(ω,g,s,u.sub.0,δG.sub.zo.sup.ob)=∫dxF.sub.r(ω,g,x,s,u.sub.0)
∫dx∫d{circumflex over (ω)}F.sub.zo,r.sup.−1({circumflex over (ω)},r,x,r,u.sub.0)δG.sub.zo.sup.ob({circumflex over (ω)},r)  (28)

(127) δG.sub.TWIN.sup.cal corresponds to the pre-stack forward modeled reflection data on the basis migrated zero-offset time domain seismic. In equation (28) the operator on the right-hand side represents then an operator acting on zero-offset data to generate data corresponding to a pre-stack acquisition geometry.

(128) Generally speaking, the computation of the computed transmitted data G.sup.cal(ω,g,s,u.sub.0) is done by solving a wave-equation or any approximation of the wave-equation in the initial velocity model u.sub.0(x). The computation of the computed scattered field δG.sub.TWIN.sup.cal(ω,g,s,u.sub.0,δG.sub.zo.sup.ob) is done by solving a scattering equation or any approximation of a scattering equation (as for instance the first order approximation known as the Born-approximation) in the velocity model u.sub.0(x) considering the high perturbation obtained from F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0).

(129) After introducing the last three terms in equation (21) relative to the cost function C.sup.RWI and by simple algebraic manipulation, the modified cost function C.sup.TWIN that is optimized by the method according to the invention is expressed by equation (29):

(130) $C^{\mathrm{TWIN}}\left(\delta {\lambda }_i\right)=\frac{1}{2}\int ds\int dg\int d\omega {.Math.\delta R^{ob}\left(\omega ,g,s,u_0\right)-\underset{i}{.Math.}{F_{\mathrm{TWIN}}\left(\omega ,g,s,u_0,F_{zo,r}^{-1}\left(\delta G_{zo}^{ob},u_0\right)\right)}_i\delta {\lambda }_i.Math.}^2$

(131) The proposed modified cost function C.sup.TWIN depends on a single kind of unknowns, the smooth model perturbation expressed by δλ.sub.i (or δm.sub.b after model expansion).

(132) Considering one linearized iteration, the minimization of C.sup.TWIN thus corresponds to solve a classical liner system y=Ax where the unknown x belong to a single class of parameters.

(133) Issues related to multi-parameters inversion between reflectivity and background velocity, such as weighting, request of the Hessian cross parameter terms and multi-modal behavior of the cost function, are largely relaxed thanks to the use of the modified cost function C.sup.TWIN.

(134) The residual term δR.sup.ob corresponds to the data residual not explained by the simulated data generated in the background medium u.sub.0.

(135) The first element G.sup.cal(ω,g,s,u.sub.0) of the simulated data corresponds to the calculated wavefield on the basis of the smooth model background u.sub.0, essentially, the direct arrival and diving waves.

(136) The second element of the modeled data δG.sub.TWIN.sup.cal(ω,g,s,u.sub.0) is the computed pre-stack reflected wavefield on the basis of the zero-offset time domain input data and of the smooth model background u.sub.0. The interpretation of this process is a mapping of the zero-offset time domain data to depth using the current background smooth model u.sub.0 (zero-offset depth migration in the true amplitude sense) followed by mapping to pre-stack time domain (using a pre-stack Born modeling).

(137) The Fréchet derivative F.sub.TWIN of the modified cost function is built from three components (see equation 26).

(138) The first two components are represented by the Frechet derivative kernel F.sub.b through the operators kernel B.sub.1 and B.sub.2. B.sub.1 simulates the effect of the smooth model perturbation on transmission or diving seismic events, whereas B.sub.2 simulates the effect of smooth model perturbation on reflected seismic events. The third component represented by the second term in the right hand side of equation (26) imposes the coupling constraint between smooth model perturbation (velocity) and high frequency model perturbation (reflectivity). That is to keep up to the first order the computed zero-offset data unchanged with respect to the perturbed background model.

(139) The method comprises estimating the smooth model perturbation by solving the normal equation associated to the modified cost function C.sup.TWIN:

(140) $\begin{array}{cc}\underset{i}{.Math.}{.Math.F_i^{\mathrm{TWIN}},F_j^{\mathrm{TWIN}}.Math.}_D\delta {\lambda }_j={.Math.F_j^{\mathrm{TWIN}},\delta R^{ob}.Math.}_D& \left(30\right)\end{array}$

(141) where .,..sub.D represents the scalar product of functions belonging to data space D.

(142) Direct inversion of the Hessian of the left hand side of equation (30) requires a very high computing effort.

(143) Advantageously, the method comprises using a preconditioning to diagonalize the Hessian to solve the normal equation associated to the modified cost function C.sup.TWIN.

(144) In a variant, the method may comprise implementing a linear or non-linear local optimization method, wherein the implementation step comprises providing at least the gradient J.sub.i.sup.TWIN of the cost-function with respect to the low frequency perturbation term. The gradient based optimization method is used to solve the normal equation associated to the modified cost function C.sup.TWIN.

(145) For example, the gradient based optimization method is a nonlinear conjugate gradient method (NLCG) or a limited memory Broyden-Fletcher-Goldfarb-Shanno method (I-BFGS).

(146) Without loss of generality, the method comprises for example using a simple steepest-descent optimization method described below.

(147) The model update at each linearized iteration will take the following classical form:
δλ.sub.i=δλ.sub.i.sup.0−α.sup.optJ.sub.i.sup.TWIN  (31)

(148) α.sup.opt is an optimal step-length of the updated model that is estimated later in the method, δλ.sub.i.sup.0 is the model expansion coefficients corresponding to the initial model

(149) J.sub.i.sup.TWIN=−F.sub.j.sup.TWIN,δR.sup.ob.sub.D is the gradient of the modified cost function that appears in the right hand side of the normal equation (30).

(150) An expanded expression of J.sub.i.sup.TWIN may be derived after a few algebraic manipulations of the previous relations:

(151) 0$\begin{array}{cc}{J}_i^{\mathrm{TWIN}}=-\int dx{\beta }_i\left(x\right)\int ds\int dg\int d\omega F_b^{*}\left(\omega ,g,x,s,u_0,F_{zo,r}^{-1}\left(\delta G_{zo}^{ob},u_0\right)\right)\delta R^{ob}\left(\omega ,g,s,u_0\right)+\int \mathrm{dx}{\beta }_i\left(x\right)\int \mathrm{dr}\int d\hat{\omega }{F}_{z0,b}^{*}\left(\hat{\omega },r,x,r,u_0,F_{zo,r}^{-1}\left(\delta G_{zo}^{ob},u_0\right)\right)\int d{\hat{x}\left(F_{z0,r}^{-1}\right)}^{*}\left(\hat{\omega },r,\hat{x},r,u_0\right)\int \mathrm{ds}\int dg\int d\omega F_r^{*}\left(\hat{\omega },g,\hat{x},s,u_0\right)\delta R^{ob}\left(\omega ,g,s,u_0\right)& \left(32\right)\end{array}$

(152) J.sub.i.sup.TWIN may be written in the symbolic form:
J.sub.i.sup.TWIN=∫dxβ.sub.i(x){J.sub.b(x)−J.sub.b.sup.zo(x)}  (33)

(153) By identification with the expanded expression of the modified cost function gradient J.sub.i.sup.TWIN, J.sub.b(x) and J.sub.b.sup.zo(x) are expressed as the following:
J.sub.b(x)=−∫ds∫dg∫dωF.sub.b*(ω,g,x,s,u.sub.0,F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0))δR.sup.ob(ω,g,s,u.sub.0)  (34)
and
J.sub.b.sup.zo(x)=−∫dr∫d{circumflex over (ω)}F.sub.zo,b*({circumflex over (ω)}r,x,r,u.sub.0,F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0))  (35)
∫d{circumflex over (x)}(F.sub.zo,r.sup.−1)*({circumflex over (ω)},r,{circumflex over (x)},r,u.sub.0)
∫ds∫dg∫dωF.sub.r*({circumflex over (ω)},g,{circumflex over (x)},s,u.sub.0)δR.sup.ob(ω,g,s,u.sub.0)

(154) The explicit expressions for the operators F.sub.r, F.sub.b and their zero-offset version F.sub.zo,r, F.sub.zo,b have already been established above. Their complex conjugates which appear in the the gradient of the modified cost function can therefore be expressed.

(155) F.sub.zo,r is the zero-offset forward Born modeling operator. Its inverse F.sub.zo,r.sup.−1 has to be understood as a reflectivity migration/inversion operator that may be approximated by a true-amplitude zero-offset depth migration or a least-squares zero-offset migration. In addition, the conjugate or the adjoint operator (F.sub.zo,r.sup.−1)* acts as a forward Born modeling operator, at least from the kinematic point of view.

(156) FIG. 5 summarizes the various contributions contained in the gradient of the modified cost function J.sub.i.sup.TWIN (equation 32).

(157) The first term J.sub.b(x) of the gradient J.sub.i.sup.TWIN has to be interpreted as the sum of the contribution of classical FWI gradient 600 (i.e. the diving sensitivity kernel or so-called “banana”) computed on the basis of smooth background velocity model and the contribution of standard RFWI gradient method, i.e. the RFWI sensitivity kernel or the so-called “two rabbit ears” as showed on FIG. 5, which connects the seismic source 16 and receiver 18 to the reflector.

(158) The second term J.sub.b.sup.zo(x) of the gradient J.sub.i.sup.TWIN has to be interpreted as the chain of the following sequence: a pre-stack depth migration of the residuals δR.sup.ob(ω,g,s,u.sub.0) for each shot through the action of ∫ ds∫ dg∫ dωF.sub.r*({circumflex over (ω)},g,{circumflex over (x)},s,u.sub.0), i.e. mapping from pre-stack time domain to depth domain which represents a first contribution 610 (FIG. 5), a demigration to zero-offset time domain through the action of ∫ d{circumflex over (x)}(F.sub.zo,r.sup.−1)*({circumflex over (ω)},r,{circumflex over (x)},r,u.sub.0), i.e. mapping from depth to zero-offset data, which represents a second contribution 620, and conversion and inversion of the obtained zero-offset data to get a gradient update along the zero-offset wave paths through the action of ∫ dr∫ d{circumflex over (ω)}F.sub.zo,b*({circumflex over (ω)},r,x,r,u.sub.0,F.sub.zo,r.sup.−1(δG.sub.zo.sup.ob,u.sub.0)) which represents a third contribution 630.

(159) For a single seismic reflection event that corresponds to a single couple of receiver-source position, the second term J.sub.b.sup.zo(x) adds a third area of contribution 640, i.e. a third “banana” or “rabbit ear” that follows the normal ray and connect the reflection point to the zero-offset position at the acquisition surface, as showed in FIG. 5.

(160) Therefore, the final outcome of the method according to the invention, is a sensitivity kernel made from four contributions 600, 610, 620, 630 (“bananas” or “rabbit ears”).

(161) As visible on FIG. 5, the contribution of a single reflected seismic event is then split along three wave-paths 640, 650, 660. Two of them 650, 660 connect the source 16 to a depth reflection point and back to the receiver 18 is spread or back-projected. The third wave-path 640 follows the normal ray from the reflection point up to the surface and back-projects along this wave-path a zero-offset residual that constrains the model update and depth positioning to be consistent.

(162) The adjoint-source term used in the computation of the gradient of the proposed approach in the previous passage corresponds to pre-stack data term depending on the type of misfit function. For a variety of misfit functions, expressions of this adjoint-source term may be obtained using a state of the art adjoint source derivation method based on Lagrange multipliers as described in “A review of the adjoint-state method for computing the gradient of a functional with geophysical applications”, Plessix, 2006, Geophysical Journal International, Volume 167, Issue 2, or automatic differentiation algorithms, described in “Automatic differentiation of algorithms: from simulation to optimization”, Corliss et al., 2002, Springer Science & Business Media. Other information may be found in “The adjoint method in seismology: I. Theory”, Fichtner et al., 2006, Physics of the Earth and Planetary Interiors, Vol. 157, Issues 1-2

(163) For instance, adjoint source term corresponding to the squared residual cost-function as defined above has for adjoint source term G.sup.ob(ω,g,s)−{G.sup.cal(ω,g,s,u)+δG.sub.TWIN.sup.cal(ω,g,s,u,δG.sub.zo.sup.ob))}

(164) Typically, the method then comprises determining the optimal step α.sup.opt.

(165) A first strategy comprises using the classical exact line search under the hypothesis that the modified cost function C.sup.TWIN may be accurately approximated by a quadratic function in the direction of the modified cost function gradient J.sub.i.sup.TWIN. Then, minimization of the optimization problem of the line-search along the gradient direction considered as a linear function versus the step length:

(166) $\begin{array}{cc}{C}^{\mathrm{TWIN}}\left(\alpha \right)=\frac{1}{2}\int ds\int dg\int d\omega {.Math.\delta R^{ob}\left(\omega ,g,s,u_0\right)-\alpha \delta R^{\mathrm{ca}l}\left(\omega ,g,s,u_0\right).Math.}^2& \left(36\right)\end{array}$
gives the classical closed formula for the optimal step-length:

(167) $\begin{array}{cc}{\alpha }^{opt}=\frac{{.Math.{\mathrm{\delta R}}^{ob}\left(\omega ,g,s,u_0\right)|\delta R^{cal}\left(\omega ,g,s,u_0\right).Math.}_D}{{.Math.{\mathrm{\delta R}}^{cal}\left(\omega ,g,s,u_0\right)|\delta R^{cal}\left(\omega ,g,s,u_0\right).Math.}_D},& \left(37\right)\end{array}$

(168) where

(169) $\begin{array}{cc}\delta R^{c\mathrm{al}}\left(\omega ,g,s,u_0\right)=-\underset{i}{.Math.}{F_{\mathrm{TWIN}}\left(\omega ,g,s,u_0,\delta m_r\right)}_i{J}_i^{\mathrm{TWIN}}& \left(38\right)\end{array}$

(170) A second strategy comprises proceeding by minimization of the modified cost function C.sup.TWIN considered as a non-linear function versus the step-length. Therefore, an inexact line search, such as proposed by Nocedal and Wright, in Numerical optimization, Springer Series in operations Research Financial Engineering, Berlin, Springer, 2006, has to be operated on the following form of the modified cost function:

(171) $\begin{array}{cc}{C}^{\mathrm{TWIN}}\left(\alpha \right)=\frac{1}{2}\int ds\int dg\int d\omega {.Math.\delta R^{ob}\left(\omega ,g,s,u_0-\alpha \underset{i}{.Math.}{\beta }_i\left(x\right)J_i^{\mathrm{TWIN}}\right).Math.}^2& \left(39\right)\end{array}$

(172) In a variant, the optimization comprises providing at least one sensitivity (F.sub.j.sup.TWIN) of a computed data derivative with respect to the low frequency perturbation term, taken as the Fréchet derivative of the computed data. The computed data includes a calculated wavefield G.sup.cal(ω,g,s,u.sub.0) calculated from the initial model u.sub.0 and a computed pre-stack reflected wavefield δG.sub.TWIN.sup.cal(ω,g,s,u.sub.0,δG.sub.zo.sup.ob) obtained on the basis of a zero-offset time domain input data and of the initial model u.sub.0.

(173) The implementation step then comprises adding a sensitivity contribution of a conventional Full Wave Inversion to a sensitivity contribution of a Reflected Full Wave Inversion and subtracting a supplementary term resulting from a chain of operators.

(174) The chain of operators includes a first operator linearly acting on the low frequency perturbation δm.sub.b(x) synthetizing a zero-offset wavefield by modeling linear scattering effects generated by the perturbation δm.sub.b(x), a second migration operator F.sub.zo,r.sup.−1 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.

(175) Finally, the method comprises generating an image of the subsurface of an area of interest from the optimal seismic velocity parameter model w.sub.opt(x) obtained after a plurality of optimization iterations. Methods for generating an image from the optimal seismic velocity parameter model w.sub.opt(x) are known in the art, for example from “Prestack Depth Migration and Velocity Model Building”, Jones, I. F., 2008, Geophysics reprint series, Society of Exploration Geophysicists.

(176) Said images may be two-dimension sections, three-dimension model block or slices extracted from the three-dimension model block.

(177) FIG. 6 is a graphical view of the convergence rate 710 of a standard Reflection Full Waveform Inversion (RFWI) and the convergence rate 720 of the method according to the invention, applied on the same synthetic 2D seismic dataset.

(178) The lowest usable frequency in the dataset is about 5 Hz and the maximal offset is 5 km.

(179) The only difference is that for the dataset inverted with the RFWI method, the effect of reflectivity-velocity coupling is not fully taken into account whereas with the method according to the invention, the effect is better taken into account.

(180) FIGS. 7 and 8 present an example of velocity vertical profile obtained respectively with the method according to the invention and with the RFWI method. The dashed line 800 corresponds to the true model and line referenced 801 corresponds to the initial velocity model.

(181) The lines referenced 805, 815, 825, 835 represent respectively the result obtained after 5, 15, 25 and 35 iterations.

(182) As it can be seen in FIG. 8, the standard RFWI method slowly updates the velocity model. With the method according to the invention (FIG. 7), the shallow part (depths from 1 km to 3 km) of the true model 800 are almost recovered whereas the RFWI method provides only a small update of the initial model 801.

(183) After 35 iterations (lines 835, FIGS. 7 and 8), the method according to the invention recovers most of the smooth features of the true model even in the deep part (4 km to 6 km). While, for the same number of iterations, the RFWI method only starts updating the deep part of the model.

(184) By comparing the misfit function values for both methods (FIG. 6), it can be noted that the misfit values obtained with the method according to the invention are about one order of magnitude smaller compared to those obtained with the standard RFWI method.

(185) Generally, the images obtained with the method according to the invention and optimal seismic velocity parameter models may be used directly to understand complex geological features and locate oil and gas reservoirs for example.

(186) In variant, said images and optimal seismic velocity parameter models are used as input for pore pressure prediction flows, time-lapse monitoring or reservoir characterization analysis.

(187) The drilling of at least a well can then be carried out based on a position determined on the image.

(188) The method according to the invention is computationally efficient.

(189) The method is not limited to seismic reflected events only. Both transmitted arrivals and reflected arrivals may be inverted simultaneously.

(190) The method is not limited to the full two-way equation solver but it remains valid for one way wave equation solver.