Kirchhoff dynamic components

[mrava87]

This PR introduces dynamic components to the previously kinematically only Kirchhoff operator.

Features:

• Added wavfilter parameter that applies a frequency-domain shaping to the wavelet as described in Safron2018 eq2.37. If wavfilter=False the v1 behavior is preserved
• Added dynamic parameter that switches between kinematic-only vs dynamic Kirchhoff operator (all of the following additions are applied only when dynamic=True)
• Added aperture to limit the spatial aperture as function of the offset/depth ratio
• Added aperturetap to include a tapering at the edges of the spatial aperture
• Added angleaperture to limit the angle aperture between source and receiver apertures - the angle is computed directly from the spatial gradients of the traveltime maps.
• Added angleaperture to limit the angle aperture between source and receiver apertures - the angle is computed directly from the spatial gradients of the traveltime maps.
• Added angleref to allow users providing a map of angles between the normal of the reflectors and the vertical axis in degrees. If provided this will be used instead of the normal when choosing whether to include a contribution or not based on angleaperture. NOTE, this has not been implemented yet.
• Added snell to allow users choosing whether the Snell's law should be used to check if a pair of source-receivers should be included or not.

[mrava87]

@cako this is really still a WIP, but it would be nice to get some early feedback from you.

[cako]

Looks really good! Some comments:

• aperture as defined by H/Z is not great on heterogeneous velocity models. I think it's worth a warning.
• snell is being implemented as a filter on the scattering angle which is double the reflection angle. As per the name, it would make more sense to threshold by reflection angle. So passing snell = 45 means reflections above 45° are killed, not reflections above 22.5°.
• I would probably implement wavfilter as an option between 2D, 2.5D and 3D. Alternatively, we could keep a bool, but use the input data so figure out which filter:
  • 2D if srcs.ndim == 2 and recs.ndim == 2 and vel.ndim == 2
  • 2.5D if (srcs.ndim == 3 or recs.ndim == 3) and vel.ndim == 2
  • 3D if srcs.ndim == 3 and recs.ndim == 3 and vel.ndim == 2
• Parameters which are cutoff filters like aperture, angleaperture, anglerefl and snell are best implemented as tapers (e.g., applying snell without a taper will cause ugly artifacts). However, the exact form of the taper is not entirely important, a Hann taper between two values is generally good enough for 99% of cases. So we could implement them as a value or a tuple representing lowpass and highcut. For example, snell = (60, 90) will allow all angles below 60°, kill all angles above 90°, and taper in between. Submitting a single angle could work as a hard threshold.
• I'm not entirely sure what anglerefl is supposed to achieve. In high frequency, the migration vector ks + kr is anticollinear with the reflector normal n, assuming the imaging point is on a reflector and full aperture (see Bleistein 2001, Chapters 4 and 5, e.g., eq. (5.1.51)). This means that a filtering/scaling/etc on the reflector normal can be applied through a filtering/scaling/etc on the migration vector as given by the scattering angle.
  A common approach is to apply cos²(refl_angle) to the summation, weakening the far angles and returning a stack without backscattering (I believe this approach is the equivalent to β of equation (5.1.21) in Bleistein 2001, whereas not applying it would be equivalent to β₁ of equation (5.1.47)). I did notice the cosines in your code, but I couldn't understand if this was what it was doing. If so, we could substitute anglerefl for something like cospower: int defined as amp *= cos(0.5 * (ang_src + ang_rec))^cospower. Setting cospower = 1 would be equivalent to cos tapering, (IMO, the β), cospower = 0 would be equivalent to no cos tapering (IMO, the β₁) and cospower = -1 would be equivalent to cos boosting (IMO, the β₂).

[mrava87]

Alright, thanks a lot for the comments, things make sense at quick read :)

Let me slowly go through them one by one and I may have a few questions along the way

[mrava87]

Looks really good! Some comments:

I managed to go through these suggestions. Here are my replies. As you will see I have not implemented all of them and created a new Github issue with a list of remaining tasks to do for the Kirchhoff operator

• aperture as defined by H/Z is not great on heterogeneous velocity models. I think it's worth a warning.

I agree. However I am not sure what could be a more flexible strategy that does not require tracing rays and checking for example take-off angles. Any suggestion?

For the moment I have made the following changes:

• modified docstring to:

Maximum allowed aperture expressed as the ratio of offset over depth. If ``None``,
no aperture limitations are introduced.

• added following warning when user is providing an aperture different from None

Aperture is currently defined as ratio of offset over depth, 
and may be not ideal for highly heterogenous media

• snell is being implemented as a filter on the scattering angle which is double the reflection angle. As per the name, it would make more sense to threshold by reflection angle. So passing snell = 45 means reflections above 45° are killed, not reflections above 22.5°.

Mmmh I am not sure you understand what I tried to do here, and I must admit the Notes were not describing it. Here is my idea, perhaps is stupid as I never worked on Kirchhoff modelling (see figure). So snell should never exceed say 2-5° or should be instead left equal None. In the first case the modelling operator would become reflection based, in the second case scatting based. I have also worked on the Notes to add more context to the 3 different apertures introduced in this PR:

* ``aperture``: the maximum allowed aperture is expressed as the ratio of
      offset over depth. This aperture limitation avoid including grazing angles
      whose contributions can introduce aliasing effects. A taper is added at the
      edges of the aperture;
    * ``angleaperture``: the maximum allowed angle aperture is expressed as the
      difference between the incident or emerging angle at every image point and
      the vertical axis (or the normal to the reflector if ``anglerefl`` is provided.
      This aperture limitation also avoid including grazing angles whose contributions
      can introduce aliasing effects. Note that for a homogenous medium and slowly varying
      heterogenous medium the offset and angle aperture limits may work in the same way;
    * ``snell``: the maximum allowed snell's angle is expressed as the absolute value of
      the sum between incident and emerging angles defined as in the ``angleaperture`` case.
      This aperture limitation is introduced to turn a scattering-based Kirchhoff engine into
      a reflection-based Kirchhoff engine where each image point is not considered as scatter
      but as a local horizontal (or straight) reflector. 

What do you think?

• I would probably implement wavfilter as an option between 2D, 2.5D and 3D. Alternatively, we could keep a bool, but use the input data so figure out which filter:

  • 2D if srcs.ndim == 2 and recs.ndim == 2 and vel.ndim == 2
  • 2.5D if (srcs.ndim == 3 or recs.ndim == 3) and vel.ndim == 2
  • 3D if srcs.ndim == 3 and recs.ndim == 3 and vel.ndim == 2

I agree. I prepared the ground for this but honestly I cannot figure out what are the σ_S and σ_R in equation 2.3.7 (neither in the Bleistein and Gray, 2002 paper). If you do, we can easily add this.

• Parameters which are cutoff filters like aperture, angleaperture, anglerefl and snell are best implemented as tapers (e.g., applying snell without a taper will cause ugly artifacts). However, the exact form of the taper is not entirely important, a Hann taper between two values is generally good enough for 99% of cases. So we could implement them as a value or a tuple representing lowpass and highcut. For example, snell = (60, 90) will allow all angles below 60°, kill all angles above 90°, and taper in between. Submitting a single angle could work as a hard threshold.

Definitely, so far I sketched only one for aperture but I am going to add two more for angleaperture andsnell now. Good idea about scalar and tuple, let me do this.

• I'm not entirely sure what anglerefl is supposed to achieve. In high frequency, the migration vector ks + kr is anticollinear with the reflector normal n, assuming the imaging point is on a reflector and full aperture (see Bleistein 2001, Chapters 4 and 5, e.g., eq. (5.1.51)). This means that a filtering/scaling/etc on the reflector normal can be applied through a filtering/scaling/etc on the migration vector as given by the scattering angle.

See my comment above regarding anglerefl. Currently the angle scaling is based on incident and emerging angles computed against the vertical axis.anglerefl would allow for a more precise computation over the reflector normals, when those are known a-priori. Does it make sense? I cannot understand how without prior knowledge of reflector orientation (and normal) one can compute these scalings correctly?

A common approach is to apply cos²(refl_angle) to the summation, weakening the far angles and returning a stack without backscattering (I believe this approach is the equivalent to β of equation (5.1.21) in Bleistein 2001, whereas not applying it would be equivalent to β₁ of equation (5.1.47)). I did notice the cosines in your code, but I couldn't understand if this was what it was doing. If so, we could substitute anglerefl for something like cospower: int defined as amp *= cos(0.5 * (ang_src + ang_rec))^cospower. Setting cospower = 1 would be equivalent to cos tapering, (IMO, the β), cospower = 0 would be equivalent to no cos tapering (IMO, the β₁) and cospower = -1 would be equivalent to cos boosting (IMO, the β₂).

Hopefully by now you see what anglerefl is meant to do. But I can of course add an extra parameter cospower doing what you suggest. However I think my cosines do not really match what you say...

My cosines are from eq 2.36 in Safron, 2018 where instead of 2cos(θ) I thought it would make more sense to have cos(θ_S)+cos(θ_R) where θ_S is the incident angle and θ_R is the emerging angle. I dont see a square in his equations... but now you say to have cos(0.5 * (ang_src + ang_rec))^cospower which even with cospower=1 to me is different from cos(θ_S)+cos(θ_R)=2cos(θ). Only in the case of of θ=θ_S=θ_R then yours becomes cos(θ)^cospower and mine becomes 2 cos(θ) like in the Safron equation. As I think you are more familiar than me on this, I am happy to implement anything you suggest. In general, no matter the equation I can off

[mrava87]

@cako, have a look at my replies (and questions ;) ) as well as the implemented changes.

If you agree I would like to merge this and continue to finalize pylops v2. I made a new issue trying to summarize what needs to be done wrt Kirchhoff operator and our discussion, so that we can try to pick this up at any time without having to scratch our head trying to remember where we saw room for improvement. Feel free to comment here or to the issue directly :)

[cako]

I agree. However I am not sure what could be a more flexible strategy that does not require tracing rays and checking for example take-off angles. Any suggestion?

Yea, apart from tracing rays theres not much one can do about it (I think). I say sticking to the warning is probably fine. But I would put this in the docstring and not as a "logging".

Mmmh I am not sure you understand what I tried to do here, and I must admit the Notes were not describing it. Here is my idea, perhaps is stupid as I never worked on Kirchhoff modelling (see figure). So snell should never exceed say 2-5° or should be instead left equal None. In the first case the modelling operator would become reflection based, in the second case scatting based. I have also worked on the Notes to add more context to the 3 different apertures introduced in this PR:

Right, I think I got it. From the gradients of traveltimes, one can compute theta_S and theta_R, defined as the source/receiver ray angles with reference to the vertical (0, -1). At a reflector location, one could also possibly compute theta_S_n and theta_R_n, defined as the source/receive ray angles with reference to the reflector normal vector n. When n = (0, -1), theta_S = theta_S_n and theta_R = theta_R_n, which is the assumption underlying the filters currently implemented.

Onto the filters:

• angleaperture is a filter applied to each theta_S and theta_R separately. Basically it says if the angle is too high (ray is close to horizontal), kill it.

• snell is a filter applied to the sum of theta_S and theta_R, and says that they cannot deviate too much from each other.

So, I'm ok with these filters and they can remain as is in my opinion. What I thought you were implementing was actually a filter on the scattering angle theta_S - theta_R instead of theta_S + theta_R. My confusion stemmed from the definition of the reflection angle being -theta_R instead of theta_R. A filter on scattering angle (or reflection angle if implemented on 0.5 * (theta_S - theta_R) simply says that rays that intersect at an angle higher than some threshold will be killed. Commonly, one would apply something like refl_angle < 60 to keep small-angle specular reflections and kill large-angle refractions/diving waves. Alternatively, one can also apply cos²(refl_angle) to every image point. This second option is the equivalent of applying a laplacian operator to the RTM image (see Appendix A of Khalil & Hoeber).

But the question still remains:

I cannot understand how without prior knowledge of reflector orientation (and normal) one can compute these scalings correctly?

Well, we know that at reflector locations, the reflector normal should be anticollinear with ks + kr, that is, the migration vector. And the angle between ks + kr and the vertical is given exactly by theta_S - theta_R. So, under two assumptions (high-frequency and the image point is at a reflector), filtering on theta_S - theta_R should be exactly equivalent to a filtering on theta_S_n - theta_R_n even when n != (0, -1)!

I agree. I prepared the ground for this but honestly I cannot figure out what are the σ_S and σ_R in equation 2.3.7 (neither in the Bleistein and Gray, 2002 paper). If you do, we can easily add this.

Yea, these guys are tricky because they are defined through ray-tracing. In 2.5D one can separate in-plane geometrical spreading from out-of-plane geometrical spreading. In Bleistein & Gray 2002, the in-plane factor is given by 1 / A₂(x, xs)A₂(x, xr) and the out-of-plane factor can be computed analytically with sqrt(σs + σr). σ is simply the ray parameter. The formula could be rederived using time or arclength (s) as the running parameter, but then the out-of-plane geometrical spreading expression would change. In any case, it is not a super important factor and I wouldn't worry about it! For completeness, Bleistein 2001 has an expression for constant-velocity 2.5D migration which substituted the sigmas by distances:

Anyways, I'm going to go over the code and approve if everything is OK. We can revisit more features later!