keurfonluu / evodcinv Goto Github PK

Inversion of dispersion curves using Evolutionary Algorithms

License: BSD 3-Clause "New" or "Revised" License

Python 100.00%

dispersion python inversion seismology geosciences surface-wave

evodcinv's Introduction

evodcinv

evodcinv is a Python library to invert surface wave dispersion data (e.g., phase velocity dispersion curves) for an isotropic layered velocity model using Evolutionary Algorithms. It relies on stochopy for the evolutionary optimizers while forward modeling is heavy-lifted by disba.

Inversion of phase velocity dispersion curve (fundamental mode).

Features

Invertible data curves:

Love-wave phase and/or group velocity dispersion curves,
Rayleigh-wave phase and/or group velocity dispersion curves,
Rayleigh-wave ellipticity (experimental).

Installation

The recommended way to install evodcinv and all its dependencies is through the Python Package Index:

pip install evodcinv --user

Otherwise, clone and extract the package, then run from the package location:

pip install . --user

To test the integrity of the installed package, check out this repository and run:

pytest

Documentation

Refer to the online documentation for detailed description of the API and examples.

Alternatively, the documentation can be built using Sphinx:

pip install -r doc/requirements.txt
sphinx-build -b html doc/source doc/build

Usage

The following example inverts a Rayleigh-wave phase velocity dispersion curve (fundamental mode).

from evodcinv import EarthModel, Layer, Curve

# Initialize model
model = EarthModel()

# Build model search boundaries from top to bottom
# First argument is the bounds of layer's thickness [km]
# Second argument is the bounds of layer's S-wave velocity [km/s]
model.add(Layer([0.001, 0.1], [0.1, 3.0]))
model.add(Layer([0.001, 0.1], [0.1, 3.0]))

# Configure model
model.configure(
    optimizer="cpso",  # Evolutionary algorithm
    misfit="rmse",  # Misfit function type
    optimizer_args={
        "popsize": 10,  # Population size
        "maxiter": 100,  # Number of iterations
        "workers": -1,  # Number of cores
        "seed": 0,
    },
)

# Define the dispersion curves to invert
# period and velocity are assumed to be data arrays
curves = [Curve(period, velocity, 0, "rayleigh", "phase")]

# Run inversion
res = model.invert(curves)
print(res)

Expected output:

--------------------------------------------------------------------------------
Best model out of 501 models (1 run)

Velocity model                                    Model parameters
----------------------------------------          ------------------------------
         d        vp        vs       rho                   d        vs        nu
      [km]    [km/s]    [km/s]   [g/cm3]                [km]    [km/s]       [-]
----------------------------------------          ------------------------------
    0.0298    0.5033    0.2055    2.0000              0.0298    0.2055    0.4000
    1.0000    2.0586    0.9935    2.0000                   -    0.9935    0.3482
----------------------------------------          ------------------------------

Number of layers: 2
Number of parameters: 5
Best model misfit: 0.0038
--------------------------------------------------------------------------------

Contributing

Please refer to the Contributing Guidelines to see how you can help. This project is released with a Code of Conduct which you agree to abide by when contributing.

evodcinv's People

Contributors

Stargazers

Watchers

evodcinv's Issues

FileNotFoundError

model.invert: concatenation of run results

That's a cool code! Thanks for posting. I found it very useful but also found some bugs.

BUG 1.
When running model.invert for maxrun>1 with split_results=False (default value) TypeError is raised.

my example:

period=np.array([ 0.66666667, 1.        , 2.        ])
velocity=np.array([ 3.49215148, 3.5623525 , 2.74965395])
curves = [Curve(period, velocity, 0, "rayleigh", "phase")]
res = model.invert(curves, maxrun=2,split_results=False)

Full error message

Run 1  |████████████████████| 100% [0:00:01 / 0:00:00] - Misfit: 0.0529
Run 2  |████████████████████| 100% [0:00:01 / 0:00:00] - Misfit: 0.0528f
�[?25hTraceback (most recent call last):

  File "/tmp/ipykernel_131271/491983976.py", line 1, in <module>
    res = model.invert(curves, maxrun=2,split_results=False)

  File "/home/wgajek/.local/lib/python3.9/site-packages/evodcinv/_model.py", line 304, in invert
    out += result

  File "/home/wgajek/.local/lib/python3.9/site-packages/evodcinv/_result.py", line 88, in __add__
    raise TypeError()

TypeError

BUG 2
What's more, prior to setting maxrun>1 running, i.e. using:
res = model.invert(curves, split_results=False)
with default maxrun, the code worked. Then I tried maxrun=2 (bug 1 reported above) and returned to maxrun=1. However, the code didn't work anymore, raising this error:

Run 1  |████████████████████| 100% [0:00:01 / 0:00:00] - Misfit: 0.0529
�[?25hTraceback (most recent call last):

  File "/tmp/ipykernel_131271/824901905.py", line 1, in <module>
    res = model.invert(curves, split_results=False)

  File "/home/wgajek/.local/lib/python3.9/site-packages/evodcinv/_model.py", line 304, in invert
    out += result

  File "/home/wgajek/.local/lib/python3.9/site-packages/evodcinv/_result.py", line 88, in __add__
    raise TypeError()

TypeError

I tried redefining the model (model = EarthModel() etc as in tutorial) but it did not solve the problem. Restarting the kernel (I'm using spyder) solved the issue (till trying maxrun>1 again).
Still, it seems solving bug 1 will also solve bug 2

System information

Operating system: Linux
Python installation (Anaconda, system, ETS): Anaconda
Version of Python: Python 3.9.7
Version of this package: v2.01

Side bug:
when trying to plot the results afterwards by using res.plot_model('rho',show='all') (for 1 run) the colors of all the lines are same so I am not able to reproduce plots from the readme. Not sure if it is another problem or you need more than 1 run for it (I guess 1 run is enough to plot density of solutions correctly). The figure is plotted then, but all lines are yellow and following error message is given:

/home/wgajek/anaconda3/lib/python3.9/site-packages/matplotlib/colors.py:1254: RuntimeWarning: invalid value encountered in true_divide
  resdat /= (vmax - vmin)
/home/wgajek/anaconda3/lib/python3.9/site-packages/matplotlib/colors.py:1254: RuntimeWarning: invalid value encountered in true_divide
  resdat /= (vmax - vmin)
/home/wgajek/anaconda3/lib/python3.9/site-packages/matplotlib/colors.py:1254: RuntimeWarning: invalid value encountered in true_divide
  resdat /= (vmax - vmin)
/home/wgajek/anaconda3/lib/python3.9/site-packages/matplotlib/colors.py:1254: RuntimeWarning: invalid value encountered in true_divide
  resdat /= (vmax - vmin)
/home/wgajek/anaconda3/lib/python3.9/site-packages/matplotlib/colors.py:1254: RuntimeWarning: invalid value encountered in true_divide
  resdat /= (vmax - vmin)

Poisson's ratio

Hello I am a postdoc trying to use this software, thanks very much for making it available.
Is there a way to prescribe a range of Poisson's ratio (for example, 0.25 to 0.35) in the configuration options ?
Is there a way to tie Vp to Vs based on some empirical relationship ?
Thanks.

Controlling size of inversion resuls

When running large inversions (for example - 2000 iterations, 10000 particles, 20 layer model), the size of the output becomes difficult to manage on small machines, especially when running through Jupyter Notebook.

Is it possible to add a flag for a "diminshed" output, something like:

Only best model per iteration
Random decimation (with a factor of choice)
Memory dump every N iterations

Thanks!

Support for Scholte wave inversion

Hi Keurfon,

More and more (including myself) people are looking into Scholte wave inversion, especially for fiber cables in the ocean.
Do you think it would be possible to extend the forward modeling/inversion to include cases where the first layer is water (S-wave = 0)?

If I just try setting Vs=0,Vp=1.5km/s,rho=1 gr/cc for the test model you provide, I get an error that seems to be coming from the Herrmann library.

Thanks!!!
Ariel

Adding a constraint for velocity increasing with depth

Geologically (especially at large scales), we often want the velocity models to increase with depth.
The current model parametrization doesn't really allow to do that (unless layers are mutually exclusive). Could it be possible to add a flag for this constraint, so that for the layer number n+1, the velocity will be between [velocity of layer n, max allowed velocity of layer n+1]?

Thanks again!

Example of plot results

Could you provide me with the script to plot the results of the inversion curve? Thanks

Multi-mode inversion yields models that do not predict the data

I seem to be getting a non-physical result when I run a 6-mode inversion of Rayleigh waves. When I lower the number of modes, the code behaves as predicted. You can see the notebook reproducing the error here:
https://www.dropbox.com/sh/9gi37zgoqg7evx0/AADD4UzR8sbPqntI9COWWaDNa?dl=0

I also join some screenshots. Sorry if this is a mistake on my side.
Picking file attached is non-standard (although a text file), there is a reader function in the jupyter notebook in the link.

Thanks!
Ariel
6_modes.txt

2): Symbol not found: _PyBytes_FromString

Hi Dr. Luu，

I got an error when i tried to run example_dc.py. Do you have any idea how did this happen?
I am using a macOS platform.

ImportError: dlopen(/Users/Mars/anaconda/lib/python3.6/site-packages/evodcinv-1.0.0-py3.6-macosx-10.7-x86_64.egg/evodcinv/_dispcurve.cpython-36m-darwin.so, 2): Symbol not found: _PyBytes_FromString Referenced from: /Users/Mars/anaconda/lib/python3.6/site-packages/evodcinv-1.0.0-py3.6-macosx-10.7-x86_64.egg/evodcinv/_dispcurve.cpython-36m-darwin.so Expected in: flat namespace in /Users/Mars/anaconda/lib/python3.6/site-packages/evodcinv-1.0.0-py3.6-macosx-10.7-x86_64.egg/evodcinv/_dispcurve.cpython-36m-darwin.so

Example not working

Description of the problem

I tried to run your example on the github page by directly copying and running. There is no "period" variable and hence the call to Curve() fails. Can you please provide a working example so I can see how your code is running? I would like to try your method with my own data, but there is no working example on your site.

Full code that generated the error

from evodcinv import EarthModel, Layer, Curve

# Initialize model
model = EarthModel()

# Build model search boundaries from top to bottom
# First argument is the bounds of layer's thickness [km]
# Second argument is the bounds of layer's S-wave velocity [km/s]
model.add(Layer([0.001, 0.1], [0.1, 3.0]))
model.add(Layer([0.001, 0.1], [0.1, 3.0]))

# Configure model
model.configure(
    optimizer="cpso",  # Evolutionary algorithm
    misfit="rmse",  # Misfit function type
    optimizer_args={
        "popsize": 10,  # Population size
        "maxiter": 100,  # Number of iterations
        "workers": -1,  # Number of cores
        "seed": 0,
    },
)

# Define the dispersion curves to invert
# period and velocity are assumed to be data arrays
curves = [Curve(period, velocity, 0, "rayleigh", "phase")]

# Run inversion
res = model.invert(curves)
print(res)

Full error message

Traceback (most recent call last):
  File "/home/dmi/PROJECT/Kumamoto/evo_example.py", line 26, in <module>
    curves = [Curve(period, velocity, 0, "rayleigh", "phase")]
NameError: name 'period' is not defined

System information

Operating system: Linux
Python installation (Anaconda, system, ETS): Conda python, installed evodcinv via pip
Version of Python: Python 3.8.12
Version of this package: 2.0.1

Output models

Hi Keurfon,
How do you output all models + misfit by iteration? I'm looking to create a figure like on the homepage with all models color-coded by misfit and Fig 7 in your 2020 paper in Prospects. Also, what's a good starting model in the event that there is no information? I am looking at using your code for near-surface investigations; < ~100 m depth

Also, how do you estimate Vp and Density of layers in the inversion? Are you using a scaling relationship between these parameters and Vs?

Has your code, both disba and evocdinv been benchmarked?