When using sample.py to perform multi-mode dispersion curve inversion, the higher mode data is not fitting well. Why?

[jhsa26]

When using sample.py to perform multi-mode dispersion curve inversion, the higher mode data is not fitting well. Why?

Full code

import numpy as np

import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalize
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

from disba import PhaseDispersion, depthplot

from evodcinv import EarthModel, Layer, Curve


# Generate synthetic data
# Example taken from <https://www.geopsy.org/wiki/index.php/Dispersion_curve_inversion>
velocity_model = np.array([
    [7.5, 500.0, 200.0, 1700.0],
    [25.0, 1350.0, 210.0, 1900.0],
    [0.0, 2000.0, 1000.0, 2500.0],
]) * 1.0e-3

pd = PhaseDispersion(*velocity_model.T)
f = np.linspace(2.0, 20.0, 50)
t = 1.0 / f[::-1]
cp_0 = pd(t, mode=0, wave="rayleigh")
cp_1 = pd(t, mode=1, wave="rayleigh")
# print(cp_1)

# 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]))

# Constant density (=2 g/cm3)
# First argument is P-wawe velocity [km/s]
density = lambda vp: 2.0

# Configure model
model.configure(
    optimizer="pso",
    misfit="rmse",
    density=density,
    optimizer_args={
        "popsize": 10,
        "maxiter": 100,
        "workers": -1,
        "seed": 0,
    },
)

# Define dispersion curves to invert
curves = [Curve(cp_0.period, cp_0.velocity, 0, "rayleigh", "phase"),
          Curve(cp_1.period, cp_1.velocity, 1, "rayleigh", "phase"),
         
         ]

# Run inversion
# See stochopy's documentation for optimizer options <https://keurfonluu.github.io/stochopy/>
res = model.invert(curves)

res = res.threshold(5 * np.min(res.misfits))
print(res.misfits)
# Plot results
fig, ax = plt.subplots(1, 3, figsize=(15, 6))

for a in ax:
    a.grid(True, linestyle=":")

zmax = 0.04
cmap = "viridis_r"

# Velocity model
res.plot_model(
    "vs",
    zmax=zmax,
    show="all",
    ax=ax[0],
    plot_args={"cmap": cmap},
)
depthplot(
    velocity_model[:, 0],
    velocity_model[:, 2],
    zmax=zmax,
    ax=ax[0],
    plot_args={
        "color": "black",
        "linewidth": 2,
        "label": "True",
    },
)
res.plot_model(
    "vs",
    zmax=zmax,
    show="best",
    ax=ax[0],
    plot_args={
        "color": "red",
        "linestyle": "--",
        "label": "Best",
    },
)
ax[0].legend(loc=1, frameon=False)

# Dispersion curve
# res.plot_curve(
#     t, 0, "rayleigh", "phase",
#     show="all",
#     ax=ax[1],
#     plot_args={
#         "type": "semilogx",
#         "xaxis": "frequency",
#         "cmap": cmap,
#     },
# )

ax[1].semilogx(
    1.0 / cp_0.period, cp_0.velocity,
    color="black",
    linewidth=2,
    label="True_mode0",
)

ax[1].semilogx(
    1.0 / cp_1.period, cp_1.velocity,'--',
    color="black",
    linewidth=2,
    label="True_mode1",

)

res.plot_curve(
    t, 0, "rayleigh", "phase",
    show="best",
    ax=ax[1],
    plot_args={
        "type": "semilogx",
        "xaxis": "frequency",
        "color": "red",
        "linestyle": "--",
        "label": "Best--mode0",
    },
)


res.plot_curve(
    t, 1, "rayleigh", "phase",
    show="best",
    ax=ax[1],
    plot_args={
        "type": "semilogx",
        "xaxis": "frequency",
        "color": "red",
        "linestyle": "--",
        "label": "Best-mode1",
    },
)


ax[1].set_xlim(2.0, 20.0)
ax[1].xaxis.set_major_formatter(ScalarFormatter())
ax[1].xaxis.set_minor_formatter(ScalarFormatter())
ax[1].legend(loc=1, frameon=False)

# Misfit
res.plot_misfit(ax=ax[2])

# Colorbar
norm = Normalize(vmin=res.misfits.min(), vmax=res.misfits.max())
smap = ScalarMappable(norm=norm, cmap=cmap)
axins = inset_axes(
    ax[1],
    width="150%",
    height="6%",
    loc="lower center",
    borderpad=-6.0,
)

cb = plt.colorbar(smap, cax=axins, orientation="horizontal")
cb.set_label("Misfit value")

plt.show()

[jhsa26]

Hi, keurf
You have done a great job for opening evodcinv library. It is a very good library for surface wave dispersion inversion task. When I try evodcinv to perform multiple-mode dispersion inversion with synthetic dataset. The higher mode is not fitting well. I have little experience for this task. It is challenge to search many parameters such like thickness, vp, vs, density, in parameter space?