Support for Scholte wave inversion

[ariellellouch]

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

[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