Robust full waveform inversion of seismic data method and device

1. A method for calculating a velocity model for a subsurface of the earth, the method comprising:

receiving measured seismic data d;

calculating predicted seismic data p;

selecting a matching filter w such that a convolution operation between the matching filter w and one of the measured seismic data d or the predicted seismic data p is equal to the other one of the measured seismic data d or the predicted seismic data p;

selecting a misfit function J that calculates (1) a distance between the matching filter w and a Dirac Delta function or (2) a travel time shift associated with the measured seismic data;

calculating a new velocity model using the misfit function J, the measured seismic data d, and the predicted seismic data p; and

generating an image of the subsurface of the earth, based on the new velocity model, for identifying resources in the subsurface,

wherein the measured seismic data d includes wavefields generated by a seismic source and the wavefields propagate through the subsurface where they are attenuated and reflected, and the attenuated and reflected wavefields are recorded by plural seismic receivers.

2. The method of claim 1 , wherein the step of selecting the misfit function J comprises:

defining the misfit function J as a product of the matching filter w and a penalty function P that varies with a time t.

3. The method of claim 2 , wherein the penalty function P is equal to an absolute value of time t.

4. The method of claim 2 , further comprising:

applying a Rytov approximation to the matching filter w.

5. The method of claim 4 , wherein the Rytov approximation rewrites the matching filter w as a sum of shifted weighted Dirac Delta functions.

6. The method of claim 5 , further comprising:

calculating an adjoint source s, which is associated with a mismatch between the measured seismic data d and the predicted seismic data p, using the misfit function J.

7. The method of claim 6 , further comprising:

calculating an adjoint state wavefield u P using the adjoint source s and a wave equation operator L; and

calculating a gradient g of the wave equation operator L.

8. The method of claim 7 , further comprising:

calculating the new velocity model using the gradient g.

9. The method of claim 1 , further comprising:

conditioning the matching filter w to be a distribution whose integral over time is unity and each value of the distribution is equal to or larger than zero.

10. A computing unit for calculating a velocity model for a subsurface of the earth, the computing unit comprising:

an interface for receiving measured seismic data d; and

a processor connected to the interface and configured to:

calculate predicted seismic data p;

select a matching filter w such that a convolution operation between the matching filter w and the measured seismic data is equal to the predicted seismic data p;

select a misfit function J that that is indicative of (1) a distance between the matching filter w and a Dirac Delta function or (2) a travel time shift associated with the measured seismic data;

calculate a new velocity model using the misfit function J, the measured seismic data d, and the predicted seismic data p; and

generating an image of the subsurface of the earth, based on the new velocity model, for identifying resources in the subsurface,

wherein the measured seismic data d includes wavefields generated by a seismic source and the wavefields propagate through the subsurface where they are attenuated and reflected, and the attenuated and reflected wavefields are recorded by plural seismic receivers.

11. The computing unit of claim 10 , wherein the processor is further configured to:

calculate the misfit function J as a product of the matching filter w and a penalty function P that varies with a time t.