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

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), a linearized perturbation of cost function with respect to δmb comprises at least a first linearized perturbation term (Fb(u0,δmr)) simulating the effect of the low frequency perturbation term on seismic data residuals and at least a second linearized perturbation term (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 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 (JiTWIN) of the cost-function with respect to the low frequency perturbation term or/and at least one sensitivity (FjTWIN) of computed data derivatives with respect to the low frequency perturbation term, the computed data including a calculated refracted/transmitted wavefield (Gcal(ω,g,s,w)) calculated from a seismic velocity model w obtained at a previous iteration and a computed pre-stack reflected wavefield (δGTWINcal(ω,g,s,w,δGzoob)) 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.