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 0 (x), a low frequency perturbation term δm b (x) and a high frequency perturbation term δm r ( x ), - determining an optimal seismic velocity parameter model w 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 b ( x ) as a parameter in the optimization of the cost function, the high frequency perturbation term δm 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 b (x), a linearized perturbation of cost function with respect to δm b comprises at least a first linearized perturbation term (F b ( u 0 , δm r )) simulating the effect of the low frequency perturbation term on seismic data residuals and at least a second linearized perturbation term F r u 0 F zo , r − 1 u 0 F zo , b u 0 δm r compensating the effect of the seismic data perturbation to keep the provided zero-offset seismic wavefield data invariant. - 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 0 (x), a low frequency perturbation term δm b (x) and a high frequency perturbation term δm r ( x ), - determining an optimal seismic velocity parameter model w 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 b ( x ) as a parameter in the optimization of the cost function, the high frequency perturbation term δm 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 b (x), a linearized perturbation of cost function with respect to δm b comprises at least a first linearized perturbation term (F b ( u 0 , δm r )) simulating the effect of the low frequency perturbation term on seismic data residuals and at least a second linearized perturbation term F r u 0 F zo , r − 1 u 0 F zo , b u 0 δm r compensating the effect of the seismic data perturbation to keep the provided zero-offset seismic wavefield data invariant. <EOS>