Synthetic Velocity Forward Model Invertion #814
Description
Problem Description
I am trying to compute travel time with pygimli's "tt.TravelTimeManager", and then invert travel time to velocity from a forward model of synthetic velocity. I keep getting an error message.
My incorrect code will output a travel time model that has some travel time in seconds as negative values.
Environment
OS : Windows (10 10.0.22631 SP0 Multiprocessor Free)
CPU(s) : 32
Machine : AMD64
Architecture : 64bit
RAM : 31.7 GiB
Environment : Jupyter
Python 3.11.10 | packaged by conda-forge | (main, Oct 16 2024, 01:17:14)
[MSC v.1941 64 bit (AMD64)]
pygimli : 1.5.3
pgcore : 1.5.0
numpy : 1.26.4
matplotlib : 3.9.2
scipy : 1.14.1
tqdm : 4.67.0
IPython : 8.29.0
pyvista : 0.44.1
Steps To Reproduce
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.traveltime as tt
import matplotlib.pyplot as plt
import scipy.io
import scipy.interpolate as interp
from scipy.interpolate import griddata
vertices = [[1, 2000], [70, 2007], [70, 1968], [1, 1961]]
polygon = mt.createPolygon(vertices, isClosed=True)
mesh = mt.createMesh(polygon, quality=2, area=10, smooth=[1, 10])
xcenter=[]
ycenter=[]
for c in mesh.cells():
tempx=c.center().x()
xcenter.append(tempx)
tempy=c.center().y()
ycenter.append(tempy)
mat_data = scipy.io.loadmat('Model3.mat')
x = mat_data['x']
z = mat_data['z']
perm = mat_data['perm']
phi = mat_data['phi']
sw = mat_data['sw']
sa = 1 - sw
phic = 0.4
BulkSolid = 30
BulkWater = 2.25
BulkAir = 0.01
BulkM = BulkSolid * (1 - phi / phic)
Bulkfluidmix = sw * BulkWater + sa * BulkAir
BulkMfluid = BulkM + (1 - BulkM / BulkSolid)**2 / (phi / Bulkfluidmix + (1 - phi) / BulkSolid - (BulkM / BulkSolid)**2)
ShearSolid = 20
Pa = 0.001
Pw = 1
pm = 2.5
ShearM = ShearSolid * (1 - phi / phic)
ShearMfluid = ShearM
DensitySR = (1 - phi) * pm + phi * (sa * Pa + sw * Pw)
Vp = np.sqrt((ShearMfluid + 4/3 * ShearMfluid) / DensitySR) * 1000
VpDG = interp.griddata((x.flatten(), z.flatten()), Vp.flatten(), (xcenter,ycenter))
numberGeophones = 70
sensors = np.linspace(1., 70, numberGeophones)
scheme = tt.createRAData(sensors)
slope = (2007 - 2000) / 70
pos = np.array(scheme.sensors())
for x in pos[:, 0]:
i = np.where(pos[:, 0] == x)
new_y = x * slope + 2000
pos[i, 1] = new_y
scheme.setSensors(pos)
print(mesh.cellSizes())
ax1, _ =pg.show(mesh, VpDG, cMin=1000, cMax=4500, logScale=False, label='v in m/s')
pg.viewer.mpl.drawSensors(ax1, scheme.sensors(), diam=1.0,
facecolor='white', edgecolor='black')
data = tt.simulate(slowness=1.0 / VpDG, scheme=scheme, mesh=mesh,
noiseLevel=0.001, noiseAbs=0.001, verbose=True)
tt.show(data)
mgr = tt.TravelTimeManager(data)
vest = mgr.invert(secNodes=2, paraMaxCellSize=20,
maxIter=15, verbose=True)
np.testing.assert_array_less(mgr.inv.inv.chi2(), 1.1)
However all of this uses data from matlab
Perm: 40 * 70 double
Phi: 40 * 70 double
sw: 40 * 70 double
x: 40 * 70 double
z: 40 * 70 double
These were randomly generated matrices