liufeng2317 / ADsurf

An automatic differentiation based (AD-based) multimodal surface wave inversion tools
MIT License
30 stars 3 forks source link

Questions regarding multimodal inversion #3

Open Geodude-93 opened 6 months ago

Geodude-93 commented 6 months ago

Hi, first of all, thank you for sharing your code. I have run the examples, have gone a bit through the code and came up with a few questions:

  1. Is the multi-mode inversion implemented yet? The "inv_param" class get a "mode" argument which default is 0 and which is apparently never used in "cpu_iter_inversion". Or do I have to invert the modes separately while specifying the argument (if I have multiple modes) and then merge the respective velocity profiles afterwards somehow ?

  2. I have attached a figure showing exemplary modes that I need to invert. They are largely overlapping in the periods. Will this be a problem for the inversion as I have to specify a "fundamental range" for the initialization of the model ? modes4github

  3. Is it possible to retrieve the predicted dispersion curve that is used to compute the loss ?

  4. This is probably a bug. When I use none of the initialization methods and just define my own starting model parameters (thickness, vp,vs,rho) I will end up in a NameError ("The initialize method can only select from [Brocher,Constant]") in the cpu_iter_inversion.

  5. Lastly, is there a paper with a more detailed description of this method and which paper do I cite when I use you code ?

Thank you!

liufeng2317 commented 6 months ago

Thansk for your question. (1) the misfit function we use is named "determinant misfit", which not need to seperate different mode of dispersion curve. (More information about this misfit function can be found at doi: 10.1111/j.1365-246X.2010.04703.x). The "mode" option just used for initialize your velocity mode (Usually empirical formulas are only suitable for dispersion curves of the fundamental order mode)

(2) I apologize for this problem, but I've been working on it recently. The correct input to the program should contain three columns (frequency, phase velocity, and mode order), but right now it doesn't do a very good job of initializing multi-order inputs. I'll be working on an updated version for easier initialization of the model.

(3) It is possible to retrieve the dispersion curves corresponding to the inversion results, just by adding a bisection search module. One of my recent work also combines a traditional loss function with a deterministic loss function, which also requires extracting the dispersion curve of the inversion result.

(4) Once you have used your own defined model parameters, you need to specify the initialization method as "Constant" to determine how vp/rho is converted during the update process (These two parameters are shown to be insensitive to changes in the dispersion curve).

(5) Our article is being submitted for review by Geophysical Journal International and I will be submitting an axriv version soon, thanks for your comments!

Geodude-93 commented 6 months ago

Thank you for the quick response and the reference. I bypassed the initializing problem feeding a second "tvs_obs" to the function only containing the fundamental mode. The multimodal inversion seems to work fine.

regarding (3): when I compute the determinant matrix for the inversion result I retrieve the following image: modes4github2 I assume the predicted phase velocities for each frequency are at the zero-crossings such that I can extract them looking for the index where the amplitude is (close to) zero ?

I would also like to use the MonteCarlo sampling approach. Can I here use the "cpu_multi_inversion" (with zero horizontal damping?) to invert from many starting models simultaneously for a single location? Could you explain how the input data for this case should look like in terms of dimension? It says " # observed dispersion data [[t,v],[t,v]] 【station,pvs_num,[t,v]]] " but I guess its not a list of lists that is requested here. Maybe you could post an example. Also it appears that the init_model_MonteCarlo instance can not be passed to the inversion directly as it does not contain a init_model instance with all the MonteCarlo samples, not sure if this is a bug or I do something wrong here. Thanks!

liufeng2317 commented 5 months ago

Thanks for your question, I have recently tweaked the program: (1) the input data supports three dimensions, [tlist,phase velocity, dipsersion order], but note that the empirical method of initializing our model only supports the 0-order dispersion curves (2) the CPU and GPU runs of the iterative inversion process are fused together, which you can determine by simply indicating the DEVICE. (3) I added the forward module, which can help you directly compute the surface wave dispersion curves for a given model. Some examples are added into the folder. (4) 0 horizontal damping is mainly used for 2D linear table array inversion to constrain the velocity model at different locations in the lateral direction, and only vertical damping is mainly used for 1D inversion. (5) The initial model of Montecarlo method is 2D data: the model is [initial number, Vp/Vs/rho/h], and the observation data is 3D data [initial number, disp_t, disp_v].

Geodude-93 commented 2 months ago

Hi, thanks for the update and congrats for the nice paper! I am back working on the project where I would like to use your code and have more questions:

The data I am using is offshore so I need to include a water layer (80m). I tried to just add a water layer as the uppermost layer to the initial models (thickness=0.08, vs=0.0, vp=1.5 and rho=1.0) When I run the inversion now the "llw-flag is activated" and I am getting dimension mismatches in "_cps/_surf96_vectorAll_gpu.py" in the dltar4_vector function in 386 "xka = omega / alpha[0]". If I debug this one I get another one at 392 at " p = ra * dpth" and afterwards one in the "var_vector" function 393/394. Hence there seems to be a systematic dimension mismatch. It seems like those functions are not adapted to handle the multiple starting models. Without adding the water layer, and hence with llw=-1, the inversion works (multi-inversion on GPU) , however the the results seem not well determined.