In [1]:
# gdapy16_00 clears all variables and import various modules
# you must run this cell first, before running any of others
# note that these examples use inverse theory techniques that
# are only covered in later chapters of the book
# History:
# 2023/04/17 - Created by W. Menke
# 2024/05/23 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patched dot product to return scalar value
# changed sum() to np.sum in gdapy16_05
# A reset deletes any previously-created variables from memory
%reset -f
# import various packages
import numpy as np # matrices & etc
from matplotlib import pyplot as plt # general plotting
from matplotlib import patches # plotting shapes
from datetime import date # date and time arithmetic
from math import exp, pi, sin, cos, tan, sqrt, floor, log10, nan, atan2 # math functions
from scipy.special import erf
import scipy.linalg as la # linear algebra functions
import scipy.signal as sg # signal processing
import os # operating system
import matplotlib.cm as cm
from matplotlib.colors import ListedColormap
import scipy.signal as sg
import scipy.stats as st
import scipy.sparse.linalg as las
from scipy import sparse
import matplotlib
# function to make a numpy N-by-1 column vector
# c=eda_cvec(v1, v2, ...) from a list of several
# array-like entities v1, v2, including a number
# a list of numbers, a tuple of numbers, an N-by-0 np array
# and a N-by-1 np array. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it
# It promotes integers to floats, and integers
# and floats to complex, by context.
# This version concatenates many argments,
# whereas c=eda_cvec1(v1) takes just one argiment.
# I recommend always using eda_cvec(v1, v2, ...)
def gda_cvec(*argv):
t = int;
Nt = 0;
for a in argv:
v = gda_cvec1(a);
N,M = np.shape(v);
Nt = Nt + N;
if( N==0 ):
continue; # skip vector of zero length
if (t==int) and isinstance(v[0,0],float):
t=float;
elif isinstance(v[0,0],complex):
t=complex;
w = np.zeros((Nt,1),dtype=t);
Nt = 0;
for a in argv:
v = gda_cvec1(a);
N,M = np.shape(v);
#w[Nt:Nt+N,0] = v[0:N,0];
w[Nt:Nt+N,0:1] = v; # patch 20230418 was #w[Nt:Nt+N,0] = v[0:N,0];
Nt = Nt + N;
return w;
# function to make a numpy N-by-1 column vector
# c=gda_cvec1(v) from entity v that is array-like,
# including a number, a list of numbers, a tuple
# of numbers, an N-by-0 np array and a N-by1 np array.
# It promotes integers to floats, and integers
# and floats to complex, by context. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it.
# This version takes just one input argment.
# whereas c=gda_cvec(v1,v2,...) concatenates
# many argiments.
def gda_cvec1(v):
if isinstance(v, int) or isinstance(v, np.int32):
w = np.zeros((1,1),dtype=int);
w[0,0] = v;
return w;
elif isinstance(v, float):
w = np.zeros((1,1),dtype=float);
w[0,0] = v;
return w;
elif isinstance(v, complex):
w = np.zeros((1,1),dtype=complex);
w[0,0] = v;
return w;
elif isinstance(v, np.ndarray):
s = np.shape(v);
if len(s) == 1:
return np.copy(np.reshape(v,(s[0],1)));
else:
[r,c]=s;
if( c==1 ):
return(np.copy(v));
elif(r==1):
return(np.copy(v.T));
else:
raise TypeError("gda_cvec: %d by %d ndarray not allowed" % (r, c));
elif isinstance(v, list):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("gda_cvec: list contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = np.array(v); # patch v -> np.array(v)
return w;
elif isinstance(v, tuple):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("gda_cvec: tuple contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = np.array(list(v)); # patch v -> np.array(list(v));
return w;
else:
raise TypeError("gda_cvec: %s not supported" % type(v));
# gda_draw function makes a "pictorial matrix equation"
# arguments are vectors, matrices and strings
# which are plotted in the order that the appear
# except that strings starting with 'title ' are plotted
# under the subseqeunt matrix or vector
# always returns a status of 1
def gda_draw(*argv):
DOCOLOR=True;
if( DOCOLOR ):
bwcmap = matplotlib.colormaps['jet'];
else:
bw = np.zeros((256,4));
v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
bw[:,0] = v;
bw[:,1] = v;
bw[:,2] = v;
bw[:,3] = np.ones(256);
bwcmap = ListedColormap(bw);
# size of plot
W = 16;
H = 4;
fig1 = plt.figure(1);
# figsize width and height in inches
fig1.set_size_inches(W,H);
ax1 = plt.subplot(1,1,1);
plt.axis([0, W, -H/2, H/2]);
plt.axis('off');
LM = W/6; # matrix width and heoght
LV = W/40; # vector width
FS = 0.12; # character width
TO = 0.4; # title vertical offset
SP = 0.2; # space between objects
LS = 0.2; # leading space
p = LS; # starting x-position
istitle=0; # flags presence of a title
for a in argv:
if isinstance(a,np.ndarray):
sh = np.shape(a);
if len(sh) == 1: # conversion to nx1 array
n = sh[0];
m = 1;
ap = a;
a = np.zeros((n,1));
a[:,0] = ap;
else:
n = sh[0];
m = sh[1];
if m==1:
pold=p;
left=p;
right=p+LV;
bottom=-LM/2;
top=LM/2;
plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
p = p+LV;
pm = (p+pold)/2;
if istitle:
plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
istitle=0;
p = p+SP;
else:
pold=p;
left=p;
right=p+LM;
bottom=-LM/2;
top=LM/2;
plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
p = p+LM;
pm = (p+pold)/2;
if istitle:
plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
istitle=0;
p = p+SP;
elif isinstance(a,str):
ns = len(a);
istitle=0;
if( ns>=6 ):
if 'title ' in a[0:6]:
istitle=1;
titlestr=a[6:];
if( istitle != 1):
plt.text(p,0,a);
p = p + ns*FS + SP;
plt.show();
return 1;
# bandpass filter, used in seismological example, but hand
# in a variety of settings involving time series
def gda_chebyshevfilt(d, Dt, flow, fhigh):
# (dout,u,v)=gda_chebyshevfilt(d, Dt, flow, fhigh);
# chebyshev IIR bandpass filter
# d - input array of data
# Dt - sampling interval
# flow - low pass frequency, Hz
# fhigh - high pass frequency, Hz
# dout - output array of data
# u - the numerator filter
# v - the denominator filter
# these filters can be used again using dout=filter(u,v,din);
# make sure input timeseries is a column vector
s = np.shape(d);
N = s[0];
if(N==1):
dd = np.zeros((N,1));
dd[:,0] = d;
else:
dd=d;
# sampling rate
rate=1/Dt;
# ripple parameter, set to ten percent
ripple=0.1;
# normalise frequency
fl=2.0*flow/rate;
fh=2.0*fhigh/rate;
# center frequency
cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );
# bandwidth
bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );
# ripple parameter factor
rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));
a = 0.5*(rpf-1.0/rpf);
b = 0.5*(rpf+1.0/rpf);
u=np.zeros((5,1));
v=np.zeros((5,1));
theta = 3*pi/4;
sr = a * cos(theta);
si = b * sin(theta);
es = sqrt(sr*sr+si*si);
tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;
v[0,0] = 1.0;
v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;
v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;
v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;
v[4,0] = (16.0 + 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw + 4.0*bw*cf*sr + cf*cf)/tmp;
tmp = 4.0*es*es*bw*bw/tmp;
u[0,0] = tmp;
u[1,0] = 0.0;
u[2,0] = -2.0*tmp;
u[3,0] = 0.0;
u[4,0] = tmp;
dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());
return (gda_cvec(dout),gda_cvec(u),gda_cvec(v));
# function to perform the multiplication FT (F v)
def gda_FTFmul(v):
# this function is used by bicongugate gradient to
# solve the generalized least squares problem Fm=F.
# Note that "F" must be called "gdaFsparse".
global gdaFsparse;
s = np.shape(v);
if(len(s)==1):
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
# weirdly, "*" multiplies sparse matrices just fine
temp = gdaFsparse*vv;
return gdaFsparse.transpose()*temp;
# function to perform the multiplication FT f
def gda_FTFrhs(f):
# this function is used by bicongugate gradient to
# solve the generalized least squares problem Fm=F.
# Note that "F" must be called "gdaFsparse".
global gdaFsparse;
return gdaFsparse.transpose()*f;
def gda_matrixmin(E):
E1 = np.min(E,axis=0);
k = np.argmin(E,axis=0);
Emin = np.min(E1);
ic = np.argmin(E1);
ir = k[ic];
return ir, ic, Emin;
def gda_Esurface( xc, yr, E, sigmad2 ):
# gda_Esurface: analysis of 2D error surface
# input:
# E: matrix of errors xc increases with columns, yr with rows
# xc, yr: corresponding vectors of (xc, yr) values; must be evenly-spaced
# sigmad2: variance of data that went into E=e'*e with e=dobs-dpre
# output:
# x0, y0: location of minimum in E
# E0: value of E at minimum
# covxy: covariance matrix
# status: 1 on success, 0 on fail
# method: quadratic approximation
# default return values
x0=0.0;
y0=0.0;
E0=0.0;
covxy=np.zeros((2,2));
status=0;
Nx,i = np.shape(xc);
Ny,i = np.shape(yr);
Dx = xc[1,0]-xc[0,0];
Dy = yr[1,0]-yr[0,0];
ir,ic,Emin = gda_matrixmin(E)
if( (ir<1) or (ir>(Ny-2)) or (ic<1) or (ic>(Nx-2)) ):
# minimum not in interior
return x0, y0, E0, covxy, status;
# quadratic model with 9 data
# E = m1 + m2*x + m3*y + m4*x*x/2 + m5*x*y + m6*y*y/2
# dE/dx = m2 + m4*x + m5*y
# dE/dy = m3 + m5*x + m6*y
# at minimum
# 0 = d1 + D2 [x,y]' so [x,y] = -inv(D2)*d1;
x = Dx*gda_cvec(-1, 0, 1, -1, 0, 1, -1, 0, 1);
y = Dy*gda_cvec(-1, -1, -1, 0, 0, 0, 1, 1, 1);
d = gda_cvec( E[ir-1,ic-1], E[ir-1,ic], E[ir-1,ic+1],
E[ir,ic-1], E[ir,ic], E[ir,ic+1],
E[ir+1,ic-1], E[ir+1,ic], E[ir+1,ic+1]);
G = np.concatenate((np.ones((9,1)),x,y,
0.5*np.multiply(x,x),np.multiply(x,y),0.5*np.multiply(y,y)),axis=1);
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,d);
m = la.solve(GTG,GTd);
d1 = gda_cvec(m[1,0], m[2,0] );
D2 = np.array( [[m[3,0], m[4,0]], [m[4,0], m[5,0]]] );
invD2 = la.inv(D2);
v = -np.matmul(invD2,d1);
x0 = xc[ic,0]+v[0,0];
y0 = yr[ir,0]+v[1,0];
#E0= m(1) + m(2)*v(1) + m(3)*v(2) + m(4)*v(1)*v(1)/2 + m(5)*v(1)*v(2) + m(6)*v(2)*v(2)/2;
E0 = m[0,0] + m[1,0]*v[0,0] + m[2,0]*v[1,0] + m[3,0]*v[0,0]*v[0,0]/2.0 + m[4,0]*v[0,0]*v[1,0] + m[5,0]*v[1,0]*v[1,0]/2.0;
covxy = (sigmad2)*2.0*invD2;
status=1;
return x0, y0, E0, covxy, status;
def gda_delay( u, Dt, t0 ):
N,i = np.shape(u);
fmax=1/(2.0*Dt);
df=fmax/(N/2);
Nf=int(N/2+1);
dw=2*pi*df;
wp = gda_cvec(dw*np.linspace(0,Nf-1,Nf));
wn = gda_cvec(-np.flipud(wp[1:Nf-1,0]));
w = gda_cvec(wp,wn);
u=np.fft.ifft(np.multiply(np.exp(complex(0,-1)*w*t0),np.fft.fft(u,axis=0)),axis=0);
u=np.real(u);
return u;
def gda_timedelay( uA, uB, Dt ):
# time delay between two similar timeseries, by
# cross-correlation followed by parabilic refinment.
# delay is positive when pulse on B is later than pulse on A
# input:
# uA, uB, two signals to be compared, as N-by-1 np.arrays
# should have zero mean
# if of unequal length, the shorter is zero padded
# Dt, sampling of both signals
# output
# tBmA_est, estimated time delay, positive when dB is delated with respect to uA
# Cmax_est, cross-correlation at estimated delay
# cross correlate
c = gda_cvec(np.correlate(uA.ravel(), uB.ravel(), mode='full'));
Nc,i = np.shape(c);
# rough estimate of lag
NA, i = np.shape(uA); # length of u
NB, i = np.shape(uB); # length of v
cmax = np.max(c);
icmax = np.argmax(c);
tBmArough = -Dt * (icmax-NB+1);
# parabolic refinement
# (c-cmax) = m(1) + m(2) dt + m(3) dt^2
# d(c-cmax)/ddt = 0 = m(2) + 2 m(3) dt
# dt = -0.5*m(2)/m(3)
# c = cmax + m(1) + m(2) dt + m(3) dt^2
Dt2 = Dt**2;
g1 = gda_cvec([1.0, 1.0, 1.0]);
g2 = Dt*gda_cvec([-1.0, 0.0, 1.0]);
g3 = Dt2*gda_cvec([1.0, 0.0, 1.0]);
G = np.concatenate( (g1,g2,g3), axis=1);
v1 = c[icmax-1,0]-cmax;
v2 = c[icmax,0]-cmax;
v3 = c[icmax+1,0]-cmax;
d = gda_cvec( [v1,v2,v3] );
m = la.solve( np.matmul(G.transpose(),G), np.matmul(G.transpose(),d) );
dt = -0.5*m[1,0]/m[2,0];
Cmax_est = cmax + m[0,0] + m[1,0]*dt + m[2,0]*Dt2;
tBmA_est = tBmArough - dt;
return tBmA_est, Cmax_est;
In [3]:
# gdama16_01
# apply Pavlis and Booker (1980) variable partitioning
# method to a simple inverse problem, randomly generated
# The script shows that the solution to the partitioned
# problem is the same as to the unpartitioned problem.
print("gdapy16_01");
N=100;
M=10;
G=np.random.uniform(low=0.0,high=1.0,size=(N,M)); # random data kernel
mtrue=np.random.uniform(low=-1.0,high=1.0,size=(M,1)); # random true model parameters
dtrue=np.matmul(G,mtrue); # true data
dobs=dtrue; # no noise in this problem, so error should be zero
# solve by least squares
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mest=la.solve(GTG,GTd);
# some solution statistics
dpre=np.matmul(G,mest);
e=dobs-dpre;
E=np.matmul(e.T,e); E=E[0,0];
print("Prediction error of full inverse problem: %f" % (E) );
l = mest-mtrue;
L = np.matmul(l.T,l); L=L[0,0];
print("Distance between full and true solutions: %f" % (L) );
# partitioned problem
# divide model parameters into two groups m=[m1',m2']'
M1=5;
M2=M-M1;
m1true=mtrue[0:M1,0:1];
m2true=mtrue[M1+1:M,0:1];
G1=G[0:N,0:M1];
G2=G[0:N,M1:M];
# svd of first data kernel
U,lam,VT = la.svd(G1);
V = VT.T;
LAM=np.diag(lam);
Up = U[0:N,0:M1];
U0 = U[0:N,M1:M];
LAMp=LAM[0:M1,0:M1];
Vp=V[0:M1,0:M1];
# first partitioned problem, Gp*m2=dp
dp = np.matmul(U0.T,dobs);
Gp = np.matmul(U0.T,G2);
# solve by least squares
GTG = np.matmul(Gp.T,Gp);
GTd = np.matmul(Gp.T,dp);
m2est = la.solve(GTG,GTd);
# second partitioned problem, Gpp*m1=dpp
dpp = np.matmul(Up.T,dobs-np.matmul(G2,m2est));
Gpp = np.matmul(LAMp,Vp.T);
# solve by least squares
GTG = np.matmul(Gpp.T,Gpp);
GTd = np.matmul(Gpp.T,dpp);
m1est = la.solve(GTG,GTd)
# some solution statistics
d12pre=np.matmul(G1,m1est)+np.matmul(G2,m2est);
e12=dobs-d12pre;
E12=np.matmul(e12.T,e12); E12=E12[0,0];
print("Prediction error of partitioned inverse problem: %f" % (E12) );
m12est = np.concatenate((m1est,m2est),axis=0);
l12 = m12est-mtrue;
L12 = np.matmul(l12.T,l12); L12=L12[0,0];
print("Distance between partitioned and true solutions: %f" % L12 );
gdapy16_01 Prediction error of full inverse problem: 0.000000 Distance between full and true solutions: 0.000000 Prediction error of partitioned inverse problem: 0.000000 Distance between partitioned and true solutions: 0.000000
In [4]:
# gdapy16_02
# Partial derivative of wavefield error for the
# acoustic case using adjoint methods.
print("gdapy16_02");
# mycase controls time of error, 1, 2, 3 allowed
# must be run for each case to produce sequence of
# figures in the book
mycase = 1;
# background slowness
s0 = 1;
m = 1;
# independent variable x
Nx = 101;
Dx = 1;
x = gda_cvec(Dx*np.linspace(0,Nx-1,Nx));
xmin = x[0,0];
xmax = x[Nx-1,0];
# independent variable y
Ny = 101;
Dy = 1;
y = gda_cvec(Dy*np.linspace(0,Ny-1,Ny));
ymin = y[0,0];
ymax = y[Ny-1,0];
# independent variable t
Nt = 2*Nx;
Dt = 1.0;
t = gda_cvec(Dt*np.linspace(0,Nt-1,Nt));
tmin=t[0,0];
tmax=t[Nt-1,0];
# source time
t0 = tmin + (tmax-tmin)/10;
# source time function is Gaussian curve
fS = np.exp( -0.1*np.power(t-t0,2) );
# source location
xS = 80.0;
yS = 51.0;
# receiver location
xR = 20.0;
yR = 51.0;
# heterogeneity location
xH = 51.0;
yH = 79.0;
# distance and traveltime from Source to Receiver
RSR = sqrt( (xS-xR)**2 + (yS-yR)**2 );
TSR = s0*RSR;
nSR = int(floor(TSR/Dt));
# distance and traveltime from Source to Heterogenety
RSH = sqrt( (xS-xH)**2 + (yS-yH)**2 );
TSH = s0*RSH;
nSH = int(floor(TSH/Dt));
# distance and traveltime from Heterogenety to Receiver
RHR = sqrt( (xH-xR)**2 + (yH-yR)**2 );
THR = s0*RHR;
nHR = int(floor(THR/Dt));
# reference field at the receiver
uR = np.zeros((Nt,1));
uR[nSR:Nt] = fS[0:Nt-nSR,0:1]/(4*pi*RSR);
# waveform error, a gaussian pulse at time te
e = np.zeros((Nt,1));
Ae = 0.1;
if( mycase == 1 ):
te = TSR+tmax/8.0;
elif( mycase == 2 ):
te = TSR+tmax/5.0;
elif( mycase == 3 ):
te = TSR+tmax/3.5;
e = e + Ae*np.exp( -0.4*np.power(t-te,2) );
# second derivative of error
edd = np.zeros((Nt,1));
edd[1:Nt-1,0:1] = (e[0:Nt-2,0:1]-2*e[1:Nt-1,0:1]+e[2:Nt,0:1])/(Dt**2);
NTPLOTS = 5;
# plot source time function
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(4,1,1);
plt.axis( [tmin, tmax, -1.25, 1.25] );
plt.plot( t, fS/np.max(np.abs(fS)), "k-", lw=2 );
plt.ylabel("fS");
# plot reference wavefield at heterogeneity
fH = np.zeros((Nt,1));
fH[nSH:Nt,0:1] = fS[0:Nt-nSH,0:1]/(4*pi*RSH);
ax1=plt.subplot(4,1,2);
plt.axis( [tmin, tmax, -1.25, 1.25] );
plt.plot( t, fH/np.max(np.abs(fH)), "k-", lw=2 );
plt.ylabel("fH");
# plot error
ax1=plt.subplot(4,1,3);
plt.axis( [tmin, tmax, -1.25, 1.25] );
plt.plot( t, e/np.max(np.abs(e)), "k-", lw=2 );
plt.ylabel("eR");
# plot second derivarive of error at heterogeneity
eddH = np.zeros((Nt,1));
eddH[0:Nt-nHR] = edd[nHR:Nt,0:1];
# plot reference wavefield and error at heterogeneity
ax1=plt.subplot(4,1,4);
a = fH/np.max(np.abs(fH));
b = eddH/np.max(np.abs(eddH));
plt.axis( [tmin, tmax, -1.25, 1.25] );
plt.plot( t, a, "k-", lw=2 );
plt.plot( t, b, "r-", lw=2 );
plt.ylabel("fH & eddH");
plt.show();
print("Caption: (Top) Source time function fS. (2nd) Reference wavefield");
print("at heterogeneity fH. (3rd) Wavefield error at receiver eR. (Bottom)");
print("Reference wavefield (fH, black) and second derivative of wavefield error");
print("at heterogeneity (eddH, red)");
# derivative dEdm
d = np.zeros((Nx,Ny));
# tabulate dEdm on 2D grid
for ix in range(Nx):
for iy in range(Ny):
# heterogeneity location
xh = x[ix,0];
yh = y[iy,0];
#source to heterogeneity parameters
RSh = sqrt( (xS-xh)**2 + (yS-yh)**2 + Dx**2);
TSh = s0*RSh;
NSh = int(floor(TSh/Dt));
# heterogeneity to receiver parameters
RhR = sqrt( (xR-xh)**2 + (yR-yh)**2 + Dx**2);
ThR = s0*RhR;
NhR = int(floor(ThR/Dt));
# reference wavefield
a = np.zeros((Nt,1));
a[NSh:Nt,0:1] = fS[0:Nt-NSh,0:1]/(4*pi*RSh);
# adjoint wavefield
b = np.zeros((Nt,1));
b[0:Nt-NhR,0:1] = 4*s0*edd[NhR:Nt,0:1]/(4*pi*RhR);
# correlate the two wavefields
d[ix,iy] = np.sum(np.multiply(a,b));
# plot 2D inage of dEdm with source, heterogeneity,
# receiver positions superimposed
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [xmin, xmax, ymin, ymax] );
dmax = np.max(np.abs(d));
left=xmin; right=xmax;
top=ymax; bottom=ymin;
plt.imshow(d, cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.plot( yS, xS, "ko", lw=4 );
plt.plot( yR, xR, "ko", lw=4 );
plt.plot( yH, xH, "ro", lw=4 );
plt.xlabel("x");
plt.ylabel("y");
plt.show();
print("Caption: The derivtive dEdm as a function of position (x,y),");
print("with source and receiver (blavk circles) and heterogenity (red circle)");
print("superimposed");
gdapy16_02
Caption: (Top) Source time function fS. (2nd) Reference wavefield at heterogeneity fH. (3rd) Wavefield error at receiver eR. (Bottom) Reference wavefield (fH, black) and second derivative of wavefield error at heterogeneity (eddH, red)
Caption: The derivtive dEdm as a function of position (x,y), with source and receiver (blavk circles) and heterogenity (red circle) superimposed
In [6]:
# gdapy16_03
# seismic migration in a medium with a uniform
# reference slowness s0, and with a reference
# (unperturbed) wavefield that is a downgoing
# plane wave.
# Note: This script uses a function gda_delay()
# that is not desribed in the text (sorry about
# that) that delays a time series by a given amount.
# It is roughly equivalent to the Numpy method
# roll(), except that it allows delays that
# are non-integral multiples of the sampling interval
# Dt. It uses a Fourier technique.
print("gdapy16_03");
# reference slowness
s0 = 1;
# axes setup. Note: Although the calulations
# are for a plane wave that starts at z=0 and
# advances in the positive z direction, I plot
# eveything upside down so that it looks "downgoing"
# (Sorry about that).
# independent variable x
xmin = -0.5;
xmax = 0.5;
Nx = 101;
Dx = (xmax-xmin)/(Nx-1);
x = gda_cvec( xmin + Dx*np.linspace(0,Nx-1,Nx) );
# independent variable z
zmin = 0.0;
zmax = 1.0;
Nz = 101;
Dz = (zmax-zmin)/(Nz-1);
z = gda_cvec( zmin + Dz*np.linspace(0,Nz-1,Nz) );
# independent variable t
# sampling Dt depends on sampling Dz
Dt = s0*Dz;
Nt = 512;
tmin = 0.0;
t = gda_cvec( tmin + Dt*np.linspace(0,Nt-1,Nt) );
tmax = t[Nt-1,0];
# in order to prevent singularities, geometrical
# spreading 1/R is replaced with 1/(R+epsi)
epsi = 4*Dx;
# source time function is a Gaussian pulse
sigmat = 2.0*Dt;
t0 = 2.0*sigmat;
p = np.exp( -np.power(t-t0,2)/(2.0*sigmat*sigmat) );
pdd = gda_cvec( 0.0, np.diff(p.ravel(),n=2), 0.0); # second derivative
pdd = (0.1/(4*pi))*pdd/np.max(pdd); # normalize
# unperturbed plane wave pulse on the grid
u0 = np.zeros((Nz,Nx,Nt));
for iz in range(Nz):
tau0 = s0*(z[iz,0]-zmin); # phase delay for a plane wave
d=gda_delay(p,Dt,tau0);
for ix in range(Nx):
# note that pwave wave experiences no geometerical spreading
u0[iz,ix,0:Nt]=d.ravel();
# array of stations (receivers), all at zmin
istas=np.arange(0,Nx,5).astype(int);
Nstas=np.shape(istas);
# point scatterer
xs = 0.00
zs = 0.25;
# scattered wavefield
tau0 = s0*(zs-zmin); # delay of planewave down to scatterer
du = np.zeros((Nz,Nx,Nt));
for iz in range(Nz):
for ix in range(Nx):
R=sqrt(((x[ix,0]-xs)**2)+((z[iz,0]-zs)**2));
Ri = 1/(R+epsi); # fudged geometrical spreading
dtau = s0*R; # delay of scattered wave
tau1 = dtau+tau0; # total delay
d=gda_delay(pdd,Dt,tau1);
du[iz,ix,0:Nt]=d.ravel()*(4*pi*Ri);
# plot wavefield
itf = 20;
fig1 = plt.figure(1,figsize=(12,8));
ax1=plt.subplot(2,2,1);
mycmap=matplotlib.colormaps['jet'];
it = itf;
plt.axis( [xmin, xmax, zmin, zmax] );
left=xmin; right=xmax;
top=zmax; bottom=zmin;
D=u0[0:Nz,0:Nx,it]+du[0:Nz,0:Nx,it]
dmax = np.max(np.abs(D)) ;
plt.imshow(D, cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.plot(xs,zmax-zs,"wo",lw=3);
plt.plot(x[istas], (zmax-0.03)*np.ones((Nstas)), "k^", lw=2 );
plt.colorbar();
it = itf*2;
ax1=plt.subplot(2,2,2);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [xmin, xmax, zmin, zmax] );
left=xmin; right=xmax;
top=zmax; bottom=zmin;
D=u0[0:Nz,0:Nx,it]+du[0:Nz,0:Nx,it]
dmax = np.max(np.abs(D)) ;
plt.imshow(D, cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.plot(xs,zmax-zs,"wo",lw=3);
plt.plot(x[istas], (zmax-0.03)*np.ones((Nstas)), "k^", lw=2 );
plt.colorbar();
it = itf*3;
ax1=plt.subplot(2,2,3);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [xmin, xmax, zmin, zmax] );
left=xmin; right=xmax;
top=zmax; bottom=zmin;
D=u0[0:Nz,0:Nx,it]+du[0:Nz,0:Nx,it]
dmax = np.max(np.abs(D)) ;
plt.imshow(D, cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.plot(xs,zmax-zs,"wo",lw=3);
plt.plot(x[istas], (zmax-0.03)*np.ones((Nstas)), "k^", lw=2 );
plt.colorbar();
it = itf*4;
ax1=plt.subplot(2,2,4);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [xmin, xmax, zmin, zmax] );
left=xmin; right=xmax;
top=zmax; bottom=zmin;
D=u0[0:Nz,0:Nx,it]+du[0:Nz,0:Nx,it]
dmax = np.max(np.abs(D));
plt.imshow(D, cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.plot(xs,zmax-zs,"wo",lw=3);
plt.plot(x[istas], (zmax-0.03)*np.ones((Nstas)), "k^", lw=2 );
plt.colorbar();
plt.xlabel("x");
plt.ylabel("z");
plt.show();
print("Caption: Plot of wavefield");
# plot record section on surface
fig2 = plt.figure(2,figsize=(12,8));
ax1=plt.subplot(1,1,1);
itp = 101;
tp = Dt*(itp-1);
plt.axis( [xmin, xmax, 0, tp] );
mycmap=matplotlib.colormaps['jet'];
left=xmin; right=xmax;
top=tp; bottom=0;
D=u0[0,0:Nx,0:itp]+du[0,0:Nx,0:itp];
dmax = np.max(np.abs(D));
plt.imshow(np.flipud(D.T), cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.xlabel("x");
plt.ylabel("t");
plt.show();
print("Caption: record section on z=0 surface");
# back-propagate du
dubT = np.zeros((Nz,Nx,Nt));
for ixr in istas.tolist(): # for each receiver on surface
xr = x[ixr,0]; # receiver position
zr = z[0,0];
dub = np.zeros((Nz,Nx,Nt)); # wavefield at receiver
a = gda_cvec( du[0,ixr,0:Nt] ); # reduce to vector
for iz in range(Nz): # back-propagate to each gridpoint
for ix in range(Nx):
R=sqrt(((x[ix,0]-xr)**2)+((z[iz,0]-zr)**2)); # distance
Ri = 1.0/(R+epsi); # fudged geometrical spreading
taub = -s0*R; # time advance
d=gda_delay(a,Dt,taub);
dub[iz,ix,0:Nt]=d.ravel()*(4*pi*Ri); # for this observation
dubT=dubT+dub; # stack all observations
# correlate
c = np.zeros((Nz,Nx));
for iz in range(Nz): # each gridpoint
for ix in range(Nx):
a = gda_cvec(u0[iz,ix,0:Nt]); # reduce to vector
add = gda_cvec(0.0,np.diff(a.ravel(),n=2),0.0); # 2nd derivative
b = gda_cvec( dubT[iz,ix,0:Nt] ); # reduce to vector
myc = np.matmul(add.T,b); myc=myc[0,0];
c[iz,ix]=myc; # correlate (= dot product)
cmax = np.max(np.abs(c));
# plot correlation
fig3 = plt.figure(3,figsize=(12,8));
ax1=plt.subplot(1,1,1);
plt.axis( [xmin, xmax, zmin, zmax] );
mycmap=matplotlib.colormaps['jet'];
left=xmin; right=xmax;
top=zmax; bottom=zmin;
plt.imshow(c/cmax, cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.plot(x[istas], (zmax-0.03)*np.ones((Nstas)), "k^", lw=2 );
plt.xlabel("x");
plt.ylabel("t");
plt.show();
print("Caption: Plot of correlation (migrated record section).");
gdapy16_03
Caption: Plot of wavefield
Caption: record section on z=0 surface
Caption: Plot of correlation (migrated record section).
In [8]:
# gdapy16_04
# How much does a small secondary pulse near a main pulse
# affect differential traveltime as determined by cross
# correlation? Both timeseries are bandpass filtered before
# cross correlation. Calculation is performed both by brute
# force and Marquering et al (1999) formula
# Note: This script uses a function gda_timedelay() that is
# not described in the text (sorry about that) that does a
# standard cross-correlation to find the delay between two
# signals, and then refines that delay by fitting a parabola
# to the peak of the cross-correlation and finding the
# peak in the parabola. This process allows the delay to
# be estimated to a resolution much better than 1Dt (often
# several orders of magnitude better).
# Note: This script uses a function gda_chebyshevfilt() that is
# not described in the text (sorry about that) to bandpass
# filter a pulse.
print("gdapy16_04");
# narrow bandpass frequencies
flow = 0.2;
fhigh = 0.3;
pulsesize = 0.2; # secondary pulse size (main pulse is 1)
# time setup
Dt = 0.01;
N=4096;
Na = int(floor(5*N/16)); # left limit of secondary pulse position
Nb = int(floor(11*N/16)); # left limit of secondary pulse position
No2 = int(floor(N/2));
t = gda_cvec(Dt*np.linspace(0,N-1,N));
t0 = t[No2,0];
ta = t[Na,0];
tb = t[Nb,0];
# reference signal u0 has just one main pulse
u0 = np.zeros((N,1));
u0[No2,0] = 1.0;
u0f,uu,vv = gda_chebyshevfilt(u0, Dt, flow, fhigh); # bandpass filter
ilist = np.arange(Na,Nb,10).astype(int); # times at which delay is calculated
ilist2 = np.arange(Na,Nb,100).astype(int); # times at which time series is plotted
Nlist = np.shape(ilist);
Nlist2 = np.shape(ilist2);
Nlist = Nlist[0]; # result above is a tupple with one element
Nlist2 = Nlist2[0];
tpulse = np.zeros((Nlist,1)); # save calculation in array
tau = np.zeros((Nlist,1));
tau2 = np.zeros((Nlist,1));
# plot of time series
fig1 = plt.figure(1,figsize=(12,8));
ax1=plt.subplot(1,1,1);
plt.axis( [10.0, 35.0, -1.0, 10.0] );
plt.title("u(t)");
plt.xlabel("time t (s)");
j=(-1); # counts iterations
for i in ilist.tolist(): # loop over pulse positions
# pulse position
j=j+1;
tpulse[j,0]=Dt*(i-No2);
# another signal u has two pulses
du = np.zeros((N,1));
du[i,0] = pulsesize;
duf,uu,vv = gda_chebyshevfilt(du, Dt, flow, fhigh); # band pass filtered perturbation
u = u0+du;
uf = u0f + duf; # bandpass filtered signal
if( i in ilist2.tolist() ):
plt.plot( t, 100*uf+0.06*j, "k-", lw=2 );
# brute force calculation of the delay by
# cross-correlation followed by parabolic refinement
# see the comments in gda_timedelay() for details.
tBmA_est, Cmax_est = gda_timedelay( u0f, uf, Dt );
tau[j,0] = tBmA_est;
cmax = Cmax_est;
# formula from Marquering et al, GJI 137, 805-815, 1999
u0fd = gda_cvec(np.diff(u0f.ravel(),n=1), 0.0)/Dt; # 1st derivative
u0fdd = gda_cvec(0.0, np.diff(u0f.ravel(),n=2), 0.0)/(Dt**2); # 2nd derivative
top = np.matmul(u0fd.T,duf); top=top[0,0];
bot = np.matmul(u0fdd.T,u0f); bot=bot[0,0]
tau2[j,0] = (Dt*top)/(Dt*bot);
plt.show();
print("Caption: Series of time series with a reference pulse (large) all with the same start");
print("time and smaller perturbation pulses with different start times.");
print(" ");
print("amplitude %.3f lowpass %.3f highpass %.3f" % (pulsesize, flow, fhigh) );
# plot of delays
scale = 80.0;
fig2 = plt.figure(2,figsize=(12,8));
ax1=plt.subplot(1,1,1);
plt.axis( [ta-t0, tb-t0, (ta-t0)/scale, (tb-t0)/scale] );
# reference lines
plt.plot( gda_cvec(0.0,0.0), gda_cvec((ta-t0)/scale, (tb-t0)/scale), "b:", lw=2 );
plt.plot( gda_cvec(ta-t0, tb-t0), gda_cvec(0.0,0.0), "b:", lw=2 );
plt.xlabel("secondary pulse position (s)");
plt.ylabel("perturbation in arrival time");
plt.plot( tpulse, tau, "k-", lw=3);
plt.plot( tpulse, tau2, "r-", lw=2);
plt.show();
print("Caption: Delay associated with small perturbation pulse, calculated using");
print("cross-correlation (black) and pertubative formula (red).");
gdapy16_04
Caption: Series of time series with a reference pulse (large) all with the same start time and smaller perturbation pulses with different start times. amplitude 0.200 lowpass 0.200 highpass 0.300
Caption: Delay associated with small perturbation pulse, calculated using cross-correlation (black) and pertubative formula (red).
In [11]:
# gdapy16_05
# bannana-doughnut kernel for the acoustic case4
print("gdapy16_05");
# reference slowness
s0 = 1;
m = 1;
# independent variable x
Nx = 101;
Dx = 1.0;
x = gda_cvec(Dx*np.linspace(0,Nx-1,Nx));
xmin = x[0,0];
xmax = x[Nx-1,0];
# independent variable y
Ny = 101;
Dy = 1;
y = gda_cvec(Dy*np.linspace(0,Ny-1,Ny));
ymin = y[0,0];
ymax = y[Ny-1,0];
# independent variable t
SSF = 10;
Nt = 2*Nx*SSF; # number of times depends on number of x's
Dt = 1.0/SSF; # time sampling depends on spatial sampling
t = gda_cvec(Dt*np.linspace(0,Nt-1,Nt));
tmin=t[0,0];
tmax=t[Nt-1,0];
# source time
t0 = tmin + (tmax-tmin)/2;
# define three frequency bands, one for each case
wny = (2*pi)/(2*Dt);
w0 = gda_cvec( 0.005*wny, 0.025*wny, 0.1*wny) ;
# loop over three cases of finite frequencies
myplt=0;
for iw in range(3):
# a band-limited source time function
fS = np.sinc(2*w0[iw,0]*(t-t0)) - np.sinc(0.5*w0[iw,0]*(t-t0));
# d/dt of source time function
fSd = np.zeros((Nt,1));
fSd[0:Nt-1,0:1] = (fS[1:Nt,0:1]-fS[0:Nt-1])/(Dt);
# d2/dt2 of source time function
fSdd = np.zeros((Nt,1));
fSdd[1:Nt-1,0:1] = (fS[0:Nt-2,0:1]-2.0*fS[1:Nt-1,0:1]+fS[2:Nt,0:1])/(Dt**2);
# source location
xS = 80.0;
yS = 51.0;
# receiver location
xR = 20.0;
yR = 51.0;
# Source to Receiver distance and time
RSR = sqrt( (xS-xR)**2 + (yS-yR)**2 );
TSR = s0*RSR;
nSR = int(floor(TSR/Dt));
# a constant used later
D = - (1/(4*pi*RSR)**2) * Dt * np.sum(np.power(fSd,2));
d = np.zeros((Nx,Ny)); # the kernel on the (x,y) plane
for ix in range (Nx): # loop over x
for iy in range(Ny): # loop over y
# location of heterogeneity
xh = x[ix,0];
yh = y[iy,0];
# Source to heterogenity distance and time
RSh = sqrt( (xS-xh)**2 + (yS-yh)**2 );
TSh = s0*RSh;
# heterogenity to receiver distance and time
RhR = sqrt( (xR-xh)**2 + (yR-yh)**2 );
ThR = s0*RhR;
# some constants
NShR = int(floor((TSh+ThR-TSR)/Dt));
c = 4.0*pi*RSR * 4.0*pi*(RSh+Dx) * 4.0*pi*(RhR+Dx);
# note, the Dx is a kludge to prevent singularies
# adjoint field
a = np.zeros((Nt,1));
a[0:Nt-NShR,0:1] = fSd[NShR:Nt,0:1];
# second derivative of reference field
b = fSdd;
# kernel is an inner product of a and b
d[ix,iy] = (2.0*s0*D)*Dt*np.sum(np.multiply(a,b))/c;
# plot kernel
myplt=myplt+1;
fig1 = plt.figure(myplt,figsize=(12,8));
ax1=plt.subplot(1,1,1);
plt.axis( [xmin, xmax, ymin, ymax] );
mycmap=matplotlib.colormaps['jet'];
left=xmin; right=xmax;
top=ymax; bottom=ymin;
dmax = np.max(np.abs(d))/3.0;
plt.imshow(d, cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );
plt.plot( yS, xS, 'wo', 'LineWidth', 4 );
plt.plot( yR, xR, 'wo', 'LineWidth', 4 );
plt.xlabel("y");
plt.ylabel("x");
plt.show();
print("Caption: Bannana-doughnut kernel for a frequency %d" % (iw) );
# plot time series
myplt=myplt+1;
fig2 = plt.figure(myplt,figsize=(12,4));
ax1=plt.subplot(1,1,1);
plt.axis( [tmin, tmax, -1.1*np.max(np.abs(fS)), 1.1*np.max(np.abs(fS))] );
plt.plot( t, fS, "k-", lw=2 );
plt.xlabel("time, t");
plt.ylabel("fS");
print("Caption: Source wavelet for frequency %d" % (iw) );
gdapy16_05
Caption: Bannana-doughnut kernel for a frequency 0 Caption: Source wavelet for frequency 0
Caption: Bannana-doughnut kernel for a frequency 1 Caption: Source wavelet for frequency 1
Caption: Bannana-doughnut kernel for a frequency 2 Caption: Source wavelet for frequency 2
In [ ]: