A problem with generating true seismograms for FWI/RTM

[NikoYura]

First, thank you for a great modeling/inversion package, you have done a great work in JUDI/Devito.

Now let me explain the problem I have noticed. Based on [05_custom_misfit.ipynb] I built my version with full Overthrust model.
I was pondering how to get a detailed image of the inhomogeneities located on the right side of the model, where there is a plane-parallel medium. I hoped to get full convergence avoiding cycle skipping by using Envelop FWI with a misfit defined in Hilbert space. While the left contrasting part turned out well, things were not so in the right plane-parallel. I was surprised to get an artifact in the right part (screenshot is attached). I started to deal with it and found the problem in the input data. I hope you can help me with this.

So, there is an original file with 97 shots for the full overtrust model. I also found a script in JUDI that generates the same data generate_data_overthrust_2D.jl in JUDI.jl/examples/software_paper/. Changing only the syntax to modern syntax (since the code doesn't run on the current version), I generated a file with the same 97 shots. Running the inversion on this file did not give me a better result. However, I was surprised to find that the original file and the generated file differed significantly. I attach screenshots of 25th shots for the generated and original file, their difference;

screenshots of the center trace (0-offset) of 25,50,75 shots superimposed on each other

and their difference in absolute values

. The graphs show that the discrepancies are strong, you can feel the displacements of the traces along the entire length. Thus, I am not sure that the inversion can eventually converge perfectly (as much as possible) in such a synthetic case. So I have several questions:

1. How is the forward modeling solved inside the inversion in JUDI (is it the same for all different functions or it may differ), for example in [05_custom_misfit.ipynb] ?
2. Does this implementation in 5th notebook of the forward modeling match with the implementation of the forward modeling inside "generate_data_overthrust_2D.jl" (d_obs = PrA_invadjoint(Ps)*q) ?
3. What is the correct way to generate shot gathers for inversion in my case? Can I use "generate_data_overthrust_2D.jl" or how to get the forward modeling version from [05_custom_misfit.ipynb]? Because I want to increase the inversion frequency and I need to generate gathers at higher frequency.
4. How sensitive is the inversion to the realization of the forward modeling if we take real data? I realize that the question is debatable, but maybe you can suggest what examples (notebooks, scripts) it is better to pay attention to ?

Thank you in advance!

[mloubout]

Hi

Thanks for reaching out. Lemme try to answer as much as possible and we can continue the discussions if some points need clarification.

So first, yes JUDI has gone through multipel updates since that dataset was generated. However you can clearly see that while there is differences due to changes in discretization and sligh modification of the formulation of the acoustic wave equation and boundary conditions, the difference are small amplitude change but do not impact the phase and traveltime. This means this is actually great for testing new algorithms since the data is very clean but you know you cannot have exact match (i.e not inverse crime).

So I want to make very clear that the error you get in the right part of the model cannot come from those small diffrerence in the data, these are actual cycle skipping issues coming from an initial model that is too far from the true model. These errors are most likely due to that very high velocity layer you get at the ocean bottom that should not be there. The notebook you pointed to shows that in most case the inversion work perfectly fine. And the fact you did not get better result creating new "inverse-crime" data shows the issue is in the algorithm.

Now on to the details:

1. The forward modeling is independent of the misfit function you choose, it only depends on the physics you choose (acoutic, density, tti)
2. Yes it does. We provide the abstractions for expresing those operation in Linear Algebra but the undrerlying modeling stays the same and defined by the Model and the Options (space_order, ...)
3. You can use that script as a base or any from the examples/script but I would first look at you inversion setup (misfit, optim algorithm, initial model) as this is what seems to introduce errors
4. This is a good question and has a lot of parameter. The discrete wave equations are tested so we know the result is correct within some margin for a given physics. In synthetic case you already know the physics is correct so you should always be able to get a good result even with discretization (different FD order, coeffs, different ABC) discrepencies.

[NikoYura]

First of all, I want to thank you for such a detailed response!
If I understood everything correctly, the difference in seismograms is caused by the updating of the forward modeling algorithm itself during these few years. In that case, the main problem is cycle skipping (I have an extremely bad starting model in the baseline scenario). Could you please suggest how you solve this problem if it is in the public domain? As I understand it, the main core of changes to the inversion process to adress this problem is the design of a "special" misfit function (with optimal transport, etc.)? Also, I've seen you use a new approach - TWRI? In any case I would be glad if you can guide me on this path, what to look out for, what examples to study?

[mloubout]

Well cycle skipping is still and "open problem". There is a lot of ongoing research and a lot of exsiting methods that mitigate it. Special misfit functions are some solutions that mitigates it, other solutions, like TWRI or extended source, also provides robustness in some cases but a lot of it will be spending time on a workflow with all the usual tricks (low-to-hogh freqency sweeps, data filtering for turning waves, adavanced imaging conditions).

So it really depends what your objective is.