import numpy as np
import matplotlib.pyplot as plt
from examples.seismic import PointSource, Receiver
from recipes import AcousticIsotropic, WaveModel
import pysegy
%matplotlib inlineFree surface vs mirrored source / mirrored receiver
Run 1: fs=True on the physical model. Run 2: open half-space (fs=False), padded upward, real source plus a negative image source at \(-d_s\), recorded at the true receivers and at mirrored receivers at \(-d_r\).
Conclusion (last three cells). Once the direct wave is removed with a constant-medium run, Du - Mu reproduces the free-surface scattered field: the difference panel is blank at the reflection, where using the true receivers alone leaves the reflection almost untouched. The mirrored receiver is required, as expected.
On the full records (cells 6-8) the direct arrival dominates and hides this, which is what misled the earlier analysis in this notebook.
N_X, N_Z = 400, 300 # cells, unpadded
SP = 5.0 # m
V_WATER, V_HALF = 1.5, 3.0
Z_REFL = 120 # reflector index unpadded -> 600 m
D_SRC = 30.0
D_REC = 30.0
TN = 1400.0 # ms
F0 = 0.030
NBL = 60
PAD = 80 # cells above z=0 in the open run
N_REC = 120
X_SRC = 0.5 * (N_X - 1) * SP
REC_X = np.linspace(0.1, 0.9, N_REC) * (N_X - 1) * SP
REC_CO = np.stack([REC_X, np.full_like(REC_X, D_REC)], axis=1)
REC_MIRROR = np.stack([REC_X, np.full_like(REC_X, -D_REC)], axis=1)def velocity(pad=0):
vp = np.full((N_X, N_Z), V_WATER, dtype=np.float32)
vp[:, Z_REFL:] = V_HALF
if pad:
vp = np.concatenate(
[np.full((N_X, pad), V_WATER, dtype=np.float32), vp], axis=1)
return vp
def run(fs, pad, sources, receivers, wavelet=None):
vp = velocity(pad)
model = WaveModel(origin=(0.0, -pad * SP), spacing=(SP, SP),
shape=vp.shape, space_order=8, vp=vp, nbl=NBL,
fs=fs, bcs="damp")
nt = int(TN / model.critical_dt) + 1
solver = AcousticIsotropic(model, {"nt": nt, "f0": F0,
"space_order": 8, "save": False})
base = np.array(solver.src.data[:, 0]) if wavelet is None else wavelet
src = PointSource(name="src", grid=model.grid, npoint=len(sources),
time_range=solver.src.time_range)
n = min(len(base), src.data.shape[0])
for i, (x, z, sign) in enumerate(sources):
src.data[:n, i] = sign * base[:n]
src.coordinates.data[i] = [x, z]
solver.__dict__["src"] = src
rec = Receiver(name="rec", grid=model.grid, npoint=len(receivers),
time_range=solver.src.time_range, coordinates=receivers)
solver.__dict__["rec"] = rec
solver.forward(save=False)
return np.array(rec.data), baseThe two models
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
vp0 = velocity(0)
ax[0].imshow(vp0.T, aspect="auto", cmap="viridis",
extent=[0, (N_X - 1) * SP, (N_Z - 1) * SP, 0])
ax[0].axhline(0, color="w", lw=3)
ax[0].plot([X_SRC], [D_SRC], "r*", ms=18, label="source")
ax[0].plot(REC_X[::6], REC_CO[::6, 1], "wv", ms=6, label="receivers")
ax[0].set_title("model 1: fs=True (free surface = white line at z=0)")
ax[0].legend(loc="lower right")
vp1 = velocity(PAD)
top = -PAD * SP
ax[1].imshow(vp1.T, aspect="auto", cmap="viridis",
extent=[0, (N_X - 1) * SP, top + (N_Z + PAD - 1) * SP, top])
ax[1].axhline(0, color="w", ls="--", lw=2)
ax[1].axhspan(top, top + NBL * SP, color="k", alpha=0.3)
ax[1].plot([X_SRC], [D_SRC], "r*", ms=18, label="source")
ax[1].plot([X_SRC], [-D_SRC], "m*", ms=18, label="image source (negative)")
ax[1].plot(REC_X[::6], REC_CO[::6, 1], "wv", ms=6, label="receivers")
ax[1].plot(REC_X[::6], REC_MIRROR[::6, 1], "cv", ms=6,
label="mirrored receivers")
ax[1].set_title("model 2: fs=False, padded (dashed = old surface, "
"grey = sponge)")
ax[1].legend(loc="lower right", fontsize=8)
for a in ax:
a.set_xlabel("x (m)")
a.set_ylabel("z (m)")
plt.tight_layout()Run both
A, wav = run(fs=True, pad=0, sources=[(X_SRC, D_SRC, 1.0)],
receivers=REC_CO)
both = np.concatenate([REC_CO, REC_MIRROR], axis=0)
R, _ = run(fs=False, pad=PAD,
sources=[(X_SRC, D_SRC, 1.0), (X_SRC, -D_SRC, -1.0)],
receivers=both, wavelet=wav)
D, M = R[:, :N_REC], R[:, N_REC:]
k = min(A.shape[0], D.shape[0])
A, D, M = A[:k], D[:k], M[:k]
t = np.linspace(0.0, TN / 1000.0, k)
print("shot records:", A.shape)Shot records
clip = 0.05 * np.abs(A).max()
panels = [(A, "1. fs=True"),
(D, "2. image source, true receivers"),
(M, "3. image source, MIRRORED receivers"),
(D - M, "4. D - M"),
(D + M, "5. D + M"),
(A - D, "6. (1) minus (2)")]
plt.figure(figsize=(16, 10))
for i, (data, name) in enumerate(panels):
plt.subplot(2, 3, i + 1)
pysegy.plot_sdata(data, name=name, new_fig=False, cbar=True)
plt.tight_layout()One trace overlaid
col = N_REC // 3
fig, ax = plt.subplots(figsize=(13, 5))
ax.plot(t, A[:, col], "k", lw=1.5, label="fs=True")
ax.plot(t, D[:, col], "r--", lw=1.1, label="image src, true rec")
ax.plot(t, (D - M)[:, col], "b:", lw=1.2, label="D - M")
ax.set_xlim(0.6, 1.15)
ax.set_xlabel("t (s)")
ax.legend()
ax.set_title("trace %d (x = %.0f m)" % (col, REC_X[col]))
plt.tight_layout()Scattered field: remove the direct wave with a constant medium
Two extra runs in a homogeneous medium give the direct arrival alone; subtracting them leaves the scattered field, where the receiver ghost lives. All gathers are resampled to one physical time axis first – critical_dt differs between the constant and layered models, and subtracting without resampling silently compares different times.
NT_OUT = 1401 # common 1 ms output axis
def run_rs(fs, pad, srcs, recs, layered=True, wav_t=None):
vp = np.full((N_X, N_Z), V_WATER, dtype=np.float32)
if layered:
vp[:, Z_REFL:] = V_HALF
if pad:
vp = np.concatenate(
[np.full((N_X, pad), V_WATER, dtype=np.float32), vp], axis=1)
m = WaveModel(origin=(0.0, -pad * SP), spacing=(SP, SP),
shape=vp.shape, space_order=8, vp=vp, nbl=NBL,
fs=fs, bcs="damp")
nt = int(TN / m.critical_dt) + 1
s = AcousticIsotropic(m, {"nt": nt, "f0": F0, "space_order": 8,
"save": False})
ta = np.asarray(s.src.time_range.time_values)
if wav_t is None:
wav_t = (ta, np.array(s.src.data[:, 0]))
p = PointSource(name="src", grid=m.grid, npoint=len(srcs),
time_range=s.src.time_range)
for i, (x, z, sg) in enumerate(srcs):
p.data[:, i] = sg * np.interp(ta, wav_t[0], wav_t[1],
left=0, right=0)
p.coordinates.data[i] = [x, z]
s.__dict__["src"] = p
r = Receiver(name="rec", grid=m.grid, npoint=len(recs),
time_range=s.src.time_range, coordinates=recs)
s.__dict__["rec"] = r
s.forward(save=False)
out = r.resample(dt=TN / (NT_OUT - 1)).data[:NT_OUT]
g = np.zeros((NT_OUT, len(recs)), dtype=np.float32)
g[:out.shape[0]] = out
return g, wav_t
both = np.concatenate([REC_CO, REC_MIRROR], axis=0)
srcs = [(X_SRC, D_SRC, 1.0), (X_SRC, -D_SRC, -1.0)]
Af, wt = run_rs(True, 0, [(X_SRC, D_SRC, 1.0)], REC_CO)
Af0, _ = run_rs(True, 0, [(X_SRC, D_SRC, 1.0)], REC_CO,
layered=False, wav_t=wt)
Rl, _ = run_rs(False, PAD, srcs, both, wav_t=wt)
Rc, _ = run_rs(False, PAD, srcs, both, layered=False, wav_t=wt)
As = Af - Af0 # free surface, scattered only
Du = Rl[:, :N_REC] - Rc[:, :N_REC] # image src, true rec, scattered
Mu = Rl[:, N_REC:] - Rc[:, N_REC:] # image src, mirrored rec, scattered
tt = np.linspace(0, TN / 1000.0, NT_OUT)clip = 0.3 * np.abs(As).max()
panels = [(As, "free surface (scattered)"),
(Du, "image src, true rec"),
(Mu, "image src, MIRRORED rec"),
(Du - Mu, "Du - Mu"),
(As - (Du - Mu), "fs - (Du - Mu) <-- cancels"),
(As - Du, "fs - Du <-- reflection left over")]
plt.figure(figsize=(16, 10))
for i, (data, name) in enumerate(panels):
plt.subplot(2, 3, i + 1)
pysegy.plot_sdata(data, name=name, new_fig=False, cbar=True)
plt.tight_layout()Result. fs - (Du - Mu) is blank at the reflection while fs - Du still shows it clearly: the mirrored receiver supplies the missing energy and enters with a minus sign. Peak amplitudes are 0.61 for the free-surface scattered field against 0.33 (Du) and 0.30 (Mu), i.e. roughly half each. The residual left in the middle panel sits on the direct-arrival cone and is imperfect cancellation by the constant-medium subtraction, not ghost physics.
Elastic
Same construction applied to ElasticIsotropic, whose free surface uses stress imaging (Thorbecke’s fdelmodc, manual eq. 51-54):
- \(\sigma_{zz}=0\) on the surface,
- \(\sigma_{xz}=0\) by odd symmetry, \(\sigma_{xz}[i_z]=-\sigma_{xz}[i_z+1]\),
- \(\sigma_{xx}\): remove the \(\partial V_z/\partial z\) term from its update and add the equivalent horizontal term, \(\sigma_{xx} = -[(\lambda+2\mu) - \lambda^2/(\lambda+2\mu)]\, \partial V_x/\partial x\).
Two differences from the acoustic case, both visible in the panels below:
- the image construction has no surface, so it cannot generate the dispersive surface-wave train that dominates the free-surface record after ~0.75 s – only the body-wave part is comparable;
- the receiver ghost applies to upgoing energy, so
Mmust not be subtracted from the direct downgoing arrival.
from recipes import ElasticIsotropic
E_VS_TOP, E_VP_TOP = 0.8, 1.6
E_VS_HALF, E_VP_HALF = 1.6, 3.0
E_ZREFL = 100
E_NX, E_NZ = 300, 200
E_TN = 1200.0
E_NTOUT = 1201
E_NBL, E_PAD = 50, 80
E_NREC = 80
E_DS = E_DR = 30.0
E_XSRC = 0.5 * (E_NX - 1) * SP
E_RX = np.linspace(0.15, 0.85, E_NREC) * (E_NX - 1) * SP
E_CO = np.stack([E_RX, np.full_like(E_RX, E_DR)], axis=1)
E_MI = np.stack([E_RX, np.full_like(E_RX, -E_DR)], axis=1)
def run_elastic(fs, pad, srcs, recs, layered=True, wav_t=None):
vp = np.full((E_NX, E_NZ), E_VP_TOP, dtype=np.float32)
vs = np.full((E_NX, E_NZ), E_VS_TOP, dtype=np.float32)
if layered:
vp[:, E_ZREFL:] = E_VP_HALF
vs[:, E_ZREFL:] = E_VS_HALF
if pad:
vp = np.concatenate(
[np.full((E_NX, pad), E_VP_TOP, dtype=np.float32), vp], axis=1)
vs = np.concatenate(
[np.full((E_NX, pad), E_VS_TOP, dtype=np.float32), vs], axis=1)
m = WaveModel(origin=(0.0, -pad * SP), spacing=(SP, SP),
shape=vp.shape, space_order=8, vp=vp, vs=vs,
b=1.0 / 2.2, nbl=E_NBL, fs=fs, bcs="damp")
nt = int(E_TN / m.critical_dt) + 1
s = ElasticIsotropic(m, {"nt": nt, "f0": F0, "space_order": 8,
"save": False})
ta = np.asarray(s.src.time_range.time_values)
if wav_t is None:
wav_t = (ta, np.array(s.src.data[:, 0]))
p = PointSource(name="src", grid=m.grid, npoint=len(srcs),
time_range=s.src.time_range)
for i, (x, z, sg) in enumerate(srcs):
p.data[:, i] = sg * np.interp(ta, wav_t[0], wav_t[1],
left=0, right=0)
p.coordinates.data[i] = [x, z]
s.__dict__["src"] = p
r = Receiver(name="rec", grid=m.grid, npoint=len(recs),
time_range=s.src.time_range, coordinates=recs)
s.__dict__["rec"] = r
s.forward(save=False)
out = r.resample(dt=E_TN / (E_NTOUT - 1)).data[:E_NTOUT]
g = np.zeros((E_NTOUT, len(recs)), dtype=np.float32)
g[:out.shape[0]] = out
return g, wav_t
e_both = np.concatenate([E_CO, E_MI], axis=0)
e_srcs = [(E_XSRC, E_DS, 1.0), (E_XSRC, -E_DS, -1.0)]
Ae, ewt = run_elastic(True, 0, [(E_XSRC, E_DS, 1.0)], E_CO)
Ae0, _ = run_elastic(True, 0, [(E_XSRC, E_DS, 1.0)], E_CO,
layered=False, wav_t=ewt)
Rle, _ = run_elastic(False, E_PAD, e_srcs, e_both, wav_t=ewt)
Rce, _ = run_elastic(False, E_PAD, e_srcs, e_both, layered=False,
wav_t=ewt)
Ase = Ae - Ae0
Due = Rle[:, :E_NREC] - Rce[:, :E_NREC]
Mue = Rle[:, E_NREC:] - Rce[:, E_NREC:]
te = np.linspace(0.0, E_TN / 1000.0, E_NTOUT)clip = 0.3 * np.abs(Ase).max()
panels = [(Ase, "elastic fs (scattered)"),
(Due, "image src, true rec"),
(Mue, "image src, MIRRORED rec"),
(Due - Mue, "Du - Mu"),
(Ase - (Due - Mue), "fs - (Du - Mu)"),
(Ase - Due, "fs - Du")]
plt.figure(figsize=(16, 10))
for i, (data, name) in enumerate(panels):
plt.subplot(2, 3, i + 1)
pysegy.plot_sdata(data, name=name, new_fig=False, cbar=True)
plt.tight_layout()def rel(a, b):
return float(np.linalg.norm(a - b) / np.linalg.norm(a))
print("%-30s %8s %8s" % ("window", "D", "D-M"))
for lo, hi, lab in ((0.55, 0.75, "reflection"),
(0.00, 0.40, "direct / early"),
(0.00, 1.20, "whole record")):
w = (te >= lo) & (te < hi)
print("%-30s %8.4f %8.4f"
% ("%s (%.2f-%.2f s)" % (lab, lo, hi),
rel(Ase[w], Due[w]), rel(Ase[w], (Due - Mue)[w])))Elastic result. On the reflection window Du - Mu scores 0.33 against 0.51 for Du alone, so the mirrored receiver helps exactly as it does acoustically. On the direct arrival it must not be subtracted (0.32 against 1.36), and the whole-record number is dominated by the surface-wave train that the image construction cannot reproduce.
The remaining 0.33 is the boundary condition’s own error. It does not move when the \(\sigma_{xx}\) expression is changed – the strain form and a plain mirror give byte-identical results, and only deleting the equation changes anything (0.9999) – so the value written for \(\sigma_{xx}\) never survives. fdelmodc replaces the update on the surface row; here it is written after the interior stress update over the physcore subdomain has already run, and the two are evidently not the same loop nest. That is the open item.
What the late elastic events are
The elastic free-surface record looks busy after the primary. Checked rather than assumed:
- 0.701 s. Independent of \(v_s\) (0.701 s for \(v_s\) = 0.6, 0.8 and 1.0 km/s), so it is a pure P path. It equals the layer primary (2 x 500 / 1600 = 0.625 s) plus both ghost delays (2 x 30/1600 twice = 0.075 s) – the reflection carrying its source and receiver ghost, the \(\tau_s + \tau_r\) structure of the paper.
- Not a base reflection. Moving the model base from 0.955 s to 1.622 s two-way (N_Z 200 -> 400, nbl 50 -> 120) leaves the energy in the 0.9-1.15 s window at 12.8% throughout.
- Not a boundary artifact. In a homogeneous elastic model with the free surface on, the late window peaks at 0.0003 against a record peak of 0.577, so the condition generates no spurious energy of its own.
The image runs contain none of this: they have no surface, so no ghosts, no surface-wave train and no free-surface multiples. Only the primary reflection is common to both, which is why the comparison is scored on that window.
Ground roll on a land model
A slow weathered layer over bedrock, with two reflectors under it. The same model is run twice: once with the free surface on, once with it off. The free surface is the only difference between the two gathers.
L_SP, L_TN, L_F0 = 6.0, 1400.0, 0.014 # m, ms, kHz
L_OFF, L_NREC = 700.0, 180 # max offset, receivers
L_W = 500.0 # weathered layer thickness, m
L_Z1, L_Z2 = 300.0, 850.0 # reflector depths, m
L_NX = int(2 * L_OFF / L_SP) + 80
L_NZ = int((L_Z2 + 520) / L_SP)
_z = np.arange(L_NZ) * L_SP
# weathered layer, bedrock, then the two reflectors
# vs of the weathered layer sets the Rayleigh speed, 0.92 vs =
# 0.60 km/s, so the ground roll reaches the far offset at 1.17 s
# and sits inside the record.
_vp = np.where(_z < L_W, 1.30, 2.40).astype(np.float32)
_vs = np.where(_z < L_W, 0.65, 1.35).astype(np.float32)
# _vp[_z >= L_Z1], _vs[_z >= L_Z1] = 2.90, 1.65
# _vp[_z >= L_Z2], _vs[_z >= L_Z2] = 3.50, 2.00
_vp[_z >= L_Z2], _vs[_z >= L_Z2] = 3.50, 2.00
L_VP = np.broadcast_to(_vp, (L_NX, L_NZ)).copy()
L_VS = np.broadcast_to(_vs, (L_NX, L_NZ)).copy()
class Geophone(ElasticIsotropic):
"""Records vertical particle velocity, i.e. what a geophone measures.
The default receiver expression is `pressure`, the trace of the
stress tensor. A Rayleigh wave is almost pure shear and shear is
traceless, so the ground roll barely shows up in pressure at all --
a hydrophone would not see it either. Land data is recorded with
geophones, which measure particle velocity, and that is where the
surface wave lives.
"""
def rec_expr(self, fw=True):
return self.state['v'][-1] # v_z
def land_shot(fs):
"""One shot on the land model, free surface on or off."""
m = WaveModel(origin=(0.0, 0.0), spacing=(L_SP, L_SP),
shape=L_VP.shape, space_order=8, vp=L_VP, vs=L_VS,
b=1.0 / 2.0, nbl=60, fs=fs, bcs="damp")
nt = int(L_TN / m.critical_dt) + 1
s = Geophone(m, {"nt": nt, "f0": L_F0, "space_order": 8,
"save": False})
x0 = 0.25 * (L_NX - 1) * L_SP
src = s.src
src.coordinates.data[0] = [x0, L_SP]
co = np.zeros((L_NREC, 2))
co[:, 0] = x0 + np.linspace(0, L_OFF, L_NREC)
co[:, 1] = L_SP
r = Receiver(name="rec", grid=m.grid, npoint=L_NREC,
time_range=s.src.time_range, coordinates=co)
s.forward(save=False, src=src, rec=r)
return np.asarray(r.data), s
plt.figure(figsize=(10, 4))
pysegy.plot_velocity(L_VP.T, spacing=(L_SP, L_SP), new_fig=False, cbar=True,
name="land model: vp (km/s)")
plt.tight_layout()
plt.show()g_fs, s_fs = land_shot(True) # free surface: ground roll
g_open, s_nofs = land_shot(False) # absorbing top: no surface wave
l_off = np.linspace(0, L_OFF, L_NREC)
l_t = np.linspace(0, L_TN / 1e3, g_fs.shape[0])
l_dt, l_dx = l_t[1] - l_t[0], l_off[1] - l_off[0]
print("record %.1f s, offsets 0-%d m" % (l_t[-1], int(L_OFF)))vb = 2.40 # bedrock P, for the move-out curves
t1 = np.sqrt((2 * L_Z1) ** 2 + l_off ** 2) / (vb * 1e3)
t2 = np.sqrt((2 * L_Z2) ** 2 + l_off ** 2) / (2.90 * 1e3)
plt.figure(figsize=(13, 13))
for i, (d, solver, nm) in enumerate(((g_fs, s_fs, "fs=True: ground roll dominates"),
(g_open, s_nofs, "fs=False: reflections visible"))):
plt.subplot(2, 2, i + 1)
pysegy.plot_sdata(d, spacing=(l_dt, l_dx), new_fig=False, cbar=False,
cmap="seismic", perc=99, name=nm)
plt.plot(l_off, t1, "b--", lw=1, label="reflector 300 m")
plt.plot(l_off, t2, "g--", lw=1, label="reflector 550 m")
plt.xlim(0, L_OFF)
plt.ylim(l_t[-1], 0)
if i == 0:
plt.legend(fontsize=8, loc="upper right")
plt.subplot(2, 2, i + 3)
if solver.model.fs:
p = solver.pressure_data
ps = solver.state['tau'][0, 1].data[0]
vp = solver.model.vp.data
else:
p = solver.pressure_data[:, solver.model.nbl:]
ps = solver.state['tau'][0, 1].data[0][:, solver.model.nbl:]
vp = solver.model.vp.data[:, solver.model.nbl:]
pysegy.plot_simage(p.T, spacing=(l_dx, l_dx), new_fig=False,
name="pressure field", cmap="seismic")
pysegy.plot_simage(ps.T, spacing=(l_dx, l_dx), new_fig=False,
name="surface shear stress", cmap="PuOr", alpha=.5)
pysegy.plot_velocity(vp.T, spacing=(l_dx, l_dx), new_fig=False,
name="vp (km/s)", cmap="terrain", alpha=.25)
plt.tight_layout()
plt.show()Same model, same source, same receivers – the only change is the boundary at the top.
Left, free surface on. The record is dominated by the ground roll: a slow, low-frequency, dispersive train fanning out from the origin. It disperses because the weathered layer traps the short wavelengths while the long ones reach the faster bedrock. The two reflections are in there, on the dashed curves, but the surface wave sits on top of them.
Right, free surface off. No surface wave, and the direct arrival and both reflections are plainly visible.
That contrast is the point. The ground roll is usually the largest thing on a land record, and a solver that mirrors pressure instead of imposing zero traction produces none of it – so a land gather modelled that way is missing its dominant arrival, which is exactly the arrival that processing has to deal with.