Revamped LSM module

[mrava87]

Motivation

PyLops has support for least-squares migration (pylops.waveeqprocessing.lsm.LSM) using a simplified Kirchhooff modelling operator (pylops.waveeqprocessing.lsm.Demigration). This PR is intended to revamp this module in a number of ways:

• Improve the overall Kirchhoff demigration/migration operators
• Introduce new demigration/migration operators implementing Born modelling of the acoustic two-way wave equation

Changes/Addition

• The pylops.waveeqprocessing.lsm.Demigration operator has been renamed into
  pylops.waveeqprocessing.kirchhoff.Kirchhoff. Its internal working has been modified taking into account the
  geometrical spreading of propagation (see the Green's function equation in Notes). To maintain the previous behaviour one can simply fill the distance table dist with ones. Moreover, to increase the flexibility of this operator and allow for the inclusion of additional weighting and aperture filters, the matvec and rmatvec are now natively implemented (instead of calling Spread
• A new acoustic wave equation engine that uses Devito's Born modelling has been added under pylops.waveeqprocessing.twoway.AcousticWave2D.
• The optional parameter engine has been added to pylops.waveeqprocessing.lsm.LSM to allow users to choose
  between the original pylops.waveeqprocessing.kirchhoff.Kirchhoff modelling operator and the new
  pylops.waveeqprocessing.twoway.AcousticWave2D modelling operator.
• Tests have been added for pylops.waveeqprocessing.twoway.AcousticWave2D, but we are currently lacking a tutorial.

Examples

• Kirchhoff LSM with data modeled using FD: https://github.com/PyLops/pylops_notebooks/blob/master/developement/LeastSquaresMigration-fdmod.ipynb
• LSRTM (small example): https://github.com/PyLops/pylops_notebooks/blob/master/developement/LSTM-PyLops_small.ipynb
• LSRTM: https://github.com/PyLops/pylops_notebooks/blob/master/developement/LSTM-PyLops.ipynb

[mrava87]

@cako, whenever you feel please take a look at this and provide any input you may have (including adding/modifying codes).

The Kirchoff propagator could be further improved, it would be good to devise a roadmap (even if we dont do it right now).

For the Devito based ones, this is just a start which could serve as a template to add more (elastic, anisotropic, etc.). Feel free to play with it and see if you like it and suggest anything that may improve it. One obvious thing I haven't tried yet is to couple it with VStack with multiprocessing to speed up computations... but I think so far I am more interested about usability and making nice examples than speed (that is easy to add later).

[cako]

Really good stuff. I think the overall structure is great and I completely agree with the name change to Kirchhoff. On thing I would suggest is to use "amp" or "weight" instead of "dist" and use that as a multiplicative factor instead of a division. The name change is related to the fact that the geometrical spreading is only dist in constant velocity. Second, the amplitude factor need not only represent geometrical spreading, but it can also incorporate other corrections. I think a generic "traveltime table & amplitude table" is more consistent. Then the change to multiplication from division makes more sense, as we rely on the user to sanitize division by zero instead of offloading that to the demigration class.

One thing to keep in mind as well, regarding nomenclature. Kirchhoff inversion is a linearization on reflectivity, Born inversion is a linearization on velocity. For Kirchhoff modeling, the formulas are, e.g.,

• Via integral (Santos et al., 2000):
  
• Via WKBJ (Jaramillo & Bleistein, 1999):
  
• Or via wave-equation (Yang & Zhang, 2019):
  
  

On the other hand, Born modeling is given by:

• Via WKBJ (Jaramillo & Bleistein , 1999):
  
• Via wave-equation (Yang & Zhang, 2019):
• 
  

The connection between them is the angle- and frequency-dependent (Yang & Zhang, 2019 - https://library.seg.org/doi/10.1190/geo2018-0438.1):

In practice what this means is that to model Kirchhoff, one needs to basically add a derivative for the result (in 3D), or alternatively, apply the appropriate rho filter in 2D or 2.5D. It also means that if one wants to do Kirchhoff inversion (instead of Born) using wave equation, we need to remove one derivative from the traditional right-hand side (-dtt u(x,t) dm(x)).

It is important to not that none of these weights are true-amplitude, they just have the correct/appropriate phase. True amplitude weights is a whole other beast we can try to implement at a later point - they should reduce the amount of required iterations for the inversions. This thesis has a lot of useful pseudocode for implementing Kirchhoff modeling via weighted sums: Safron, 2018.

[mrava87]

Thanks, this is great feedback :)

I have been recently trying to get my head around modelling/inversion with reflectivity and velocity perturbation - this sums it very well!

Although I think it is good to have an overall LSM wrapper I will make sure to modify the documentation to remark this point (I’ll ask you to read it once ready).

And I much prefer your idea of using amp and asking the user to handle 1/0 cases, I’ll change that right away. This way we can later on modify such amplitude weights with more than ‘geometrical spreading’ without having to affect the function signature.

I will be nice to have true-amplitude modelling. Let me take a look at that thesis, and this is something we can ask new users so they learn the inner working of pylops whilst doing something meaningful at the same time (could be good for some of my students doing imaging courses)

[cako]

Feel free to update and merge!