Support for Scholte wave inversion about evodcinv

[keurfonluu]

HI @ariellellouch,

I tested disba's example with a water layer, and the forward modeling seems to be working:

import numpy as np
from disba import PhaseDispersion
import matplotlib.pyplot as plt

velocity_model = np.array([
   [5.0, 1.5, 0.0, 1.00],
   [10.0, 7.00, 3.50, 2.00],
   [10.0, 6.80, 3.40, 2.00],
   [10.0, 7.00, 3.50, 2.00],
   [10.0, 7.60, 3.80, 2.00],
   [10.0, 8.40, 4.20, 2.00],
   [10.0, 9.00, 4.50, 2.00],
   [10.0, 9.40, 4.70, 2.00],
   [10.0, 9.60, 4.80, 2.00],
   [10.0, 9.50, 4.75, 2.00],
])
t = np.logspace(0.0, 3.0, 100)

pd = PhaseDispersion(*velocity_model.T)
cp = pd(t, mode=0, wave="rayleigh")

plt.semilogx(cp.period, cp.velocity)

But indeed, water layers are not supported yet in evodcinv (it's not possible to set Vp). For the water layer, I guess we are only looking for its thickness (Vp and density are known I guess).

[ariellellouch]

It also works for me in disba, which is already super useful!
I think this is a reasonable assumption - changes in water velocity/density are very small and would require multiple layers.

Maybe adding a flag of "first_layer_water" or something like that would be a good wraparound without opening the option to change Vp for the inversion?

Ariel

[keurfonluu]

Hi @ariellellouch,

Please update to v2.2.2. The fix was pretty straightforward in the end, I simply added a condition when transforming the model parameter vector to a velocity model. The first layer's Vp is determined by its Vs as follows:

• if Vs > 0, Vp determined using Poisson's ratio (no water layer)
• if Vs = 0, Vp = 1.5 km/s (default value)
• if Vs < 0, Vp = -Vs km/s

import numpy as np
import matplotlib.pyplot as plt

from disba import PhaseDispersion
from evodcinv import EarthModel, Layer, Curve, factory

# Generate synthetic data
velocity_model = np.array([
   [5.0, 1.5, 0.0, 1.0],
   [30.0, 4.6, 2.6, 2.0],
])
f = np.linspace(0.05, 1.0, 100)
t = 1.0 / f[::-1]
pd = PhaseDispersion(*velocity_model.T)
cp = pd(t, mode=0, wave="rayleigh")

# Set up model
model = EarthModel()
model.add(Layer([3.0, 7.0], [-2.0, 1.0]))  # Invert for water layer's Vp or first layer's Vs
# model.add(Layer([3.0, 7.0], 0.0))  # Fix Vp to default value (1.5 km/s)
# model.add(Layer([3.0, 7.0], -1.45))  # Fix Vp to 1.45 km/s
model.add(Layer(30.0, [1.5, 3.5]))

model.configure(
    optimizer="cpso",
    misfit="rmse",
    density=lambda vp: 2.0,
    optimizer_args={
        "popsize": 10,
        "maxiter": 200,
        "workers": -1,
        "seed": 0,
    },
)
curves = [Curve(cp.period, cp.velocity, 0, "rayleigh", "phase")]

# Invert
res = model.invert(curves)
print(res)

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

Velocity model                                    Model parameters
----------------------------------------          ------------------------------
         d        vp        vs       rho                   d        vs        nu
      [km]    [km/s]    [km/s]   [g/cm3]                [km]    [km/s]       [-]
----------------------------------------          ------------------------------
    5.0274    1.4991    0.0000    1.0000              5.0274   -1.4991    0.3243
    1.0000    4.7551    2.5827    2.0000                   -    2.5827    0.2908
----------------------------------------          ------------------------------

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

[ariellellouch]

That's a great solution!
Thanks