% gdama16
% W. Menke (2024), Geophysical Data Analysis and Inverse Theory
% with MATLAB(R) or Python, 5th Edition
% Code History:
% 2022-10-01, Imported from 4th Edition by W. Menke
% 2023-04-08, Edit and integration with text
% 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.
clear all;
fprintf('gdama16_01\n');
N=100;
M=10;
G=random('uniform',0,1,N,M); % random data kernel
mtrue=random('uniform',-1,1,M,1); % random true model parameters
dtrue=G*mtrue; % true data
dobs=dtrue; % no noise in this problem, so error should be zero
% solve by least squares
mest=(G'*G)\(G'*dobs);
% some solution statistics
dpre=G*mest;
e=dobs-dpre;
E=e'*e;
fprintf('Prediction error of full inverse problem: %f\n', E );
l = mest-mtrue;
L = l'*l;
fprintf('Distance between full and true solutions: %f\n', L );
% partitioned problem
% divide model parameters into two groups m=[m1',m2']'
M1=5;
M2=M-M1;
m1true=mtrue(1:M1);
m2true=mtrue(M1+1:end);
G1=G(:,1:M1);
G2=G(:,M1+1:end);
% svd of first data kernel
[U,LAMBDA,V] = svd(G1);
lambda=diag(LAMBDA);
Up = U(:,1:M1);
U0 = U(:,M1+1:end);
LAMBDAp=LAMBDA(1:M1,1:M1);
Vp=V(1:M1,1:M1);
% first partitioned problem, Gp*m2=dp
dp = U0'*dobs;
Gp = U0'*G2;
% solve by least squares
m2est = (Gp'*Gp)\(Gp'*dp);
% second partitioned problem, Gpp*m1=dpp
dpp = Up'*(dobs-G2*m2est);
Gpp = LAMBDAp*Vp';
% solve by least squares
m1est = (Gpp'*Gpp)\(Gpp'*dpp);
% some solution statistics
d12pre=G1*m1est+G2*m2est;
e12=dobs-d12pre;
E12=e12'*e12;
fprintf('Prediction error of partitioned inverse problem: %f\n', E12 );
m12est = [m1est', m2est']';
l12 = m12est-mtrue;
L12 = l12'*l12;
fprintf('Distance between partitioned and true solutions: %f\n', L12 );
% gdama16_02
% Partial derivative of wavefield error for the
% acoustic case using adjoint methods.
clear all;
fprintf('gdama16_02\n');
% 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 = Dx*[0:Nx-1]';
xmin = x(1);
xmax = x(end);
% independent variable y
Ny = 101;
Dy = 1;
y = Dy*[0:Ny-1]';
ymin = y(1);
ymax = y(end);
% independent variable t
Nt = 2*Nx;
Dt = 1;
t = Dt*[0:Nt-1]';
tmin=t(1);
tmax=t(end);
% source time
t0 = tmin + (tmax-tmin)/10;
% source time function is Gaussian curve
fS = exp( -0.1*(t-t0).^2 );
% source location
xS = 80;
yS = 51;
% receiver location
xR = 20;
yR = 51;
% heterogeneity location
xH = 51;
yH = 79;
% distance and traveltime from Source to Receiver
RSR = sqrt( (xS-xR)^2 + (yS-yR)^2 );
TSR = s0*RSR;
nSR = floor( TSR/Dt );
% distance and traveltime from Source to Heterogenety
RSH = sqrt( (xS-xH)^2 + (yS-yH)^2 );
TSH = s0*RSH;
nSH = floor( TSH/Dt );
% distance and traveltime from Heterogenety to Receiver
RHR = sqrt( (xH-xR)^2 + (yH-yR)^2 );
THR = s0*RHR;
nHR = floor( THR/Dt );
% reference field at the receiver
uR = zeros(Nt,1);
uR(nSR+1:Nt) = fS(1:Nt-nSR)/(4*pi*RSR);
% waveform error, a gaussian pulse at time te
e = zeros(Nt,1);
Ae = 0.1;
if( mycase == 1 )
te = TSR+tmax/8;
elseif( mycase == 2 )
te = TSR+tmax/5;
elseif( mycase == 3 )
te = TSR+tmax/3.5;
end
e = e + Ae*exp( -0.4*(t-te).^2 );
% second derivative of error
edd = zeros(Nt,1);
edd(2:Nt-1) = (e(1:Nt-2)-2*e(2:Nt-1)+e(3:Nt))/(Dt^2);
NTPLOTS = 5;
% plot source time function
figure();
clf;
subplot(4,1,1);
hold on;
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -1.25, 1.25] );
plot( t, fS/max(abs(fS)), 'k-', 'LineWidth', 2 );
ylabel('fS');
% plot reference wavefield at heterogeneity
fH = zeros(Nt,1);
fH(nSH+1:Nt) = fS(1:Nt-nSH)/(4*pi*RSH);
subplot(4,1,2);
hold on;
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -1.25, 1.25] );
plot( t, fH/max(abs(fH)), 'k-', 'LineWidth', 2 );
ylabel('fH');
% plot error
subplot(4,1,3);
hold on;
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -1.25, 1.25] );
plot( t, e/max(abs(e)), 'k-', 'LineWidth', 2 );
ylabel('eR');
% plot second derivarive of error at heterogeneity
eddH = zeros(Nt,1);
eddH(1:Nt-nHR) = edd(nHR+1:Nt);
% plot reference wavefield and error at heterogeneity
subplot(4,1,4);
hold on;
set(gca,'LineWidth',2);
hold on;
a = fH/max(abs(fH));
b = eddH/max(abs(eddH));
axis( [tmin, tmax, -1.25, 1.25] );
plot( t, a, 'k-', 'LineWidth', 2 );
plot( t, b, 'r-', 'LineWidth', 2 );
ylabel('fH & eddH');
fprintf('Caption: (Top) Source time function fS. (2nd) Reference wavefield\n');
fprintf('at heterogeneity fH. (3rd) Wavefield error at receiver eR. (Bottom)\n');
fprintf('Reference wavefield (fH, black) and second derivative of wavefield error\n');
fprintf('at heterogeneity (eddH, red)\n');
% derivative dEdm
d = zeros(Nx,Ny);
% tabulate dEdm on 2D grid
for ix = [1:Nx]
for iy = [1:Ny]
% heterogeneity location
xh = x(ix);
yh = y(iy);
% source to heterogeneity parameters
RSh = sqrt( (xS-xh)^2 + (yS-yh)^2 + Dx^2);
TSh = s0*RSh;
NSh = floor(TSh/Dt);
% heterogeneity to receiver parameters
RhR = sqrt( (xR-xh)^2 + (yR-yh)^2 + Dx^2);
ThR = s0*RhR;
NhR = floor(ThR/Dt);
% reference wavefield
a = zeros(Nt,1);
a(NSh+1:Nt) = fS(1:Nt-NSh)/(4*pi*RSh);
% adjoint wavefield
b = zeros(Nt,1);
b(1:Nt-NhR) = 4*s0*edd(NhR+1:Nt)/(4*pi*RhR);
% correlate the two wavefields
d(ix,iy) = sum( a .* b );
end
end
% plot 2D inage of dEdm with source, heterogeneity,
% receiver positions superimposed
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
colormap('jet');
hold on
axis ij
axis( [xmin, xmax, ymin, ymax] );
dmax = max(abs(d(:)));
imagesc(d,[-dmax,dmax]);
plot( yS, xS, 'ko', 'LineWIdth', 4 );
plot( yR, xR, 'ko', 'LineWIdth', 4 );
plot( yH, xH, 'ro', 'LineWIdth', 4 );
xlabel('x');
ylabel('y');
fprintf('Caption: The derivtive dEdm as a function of position (x,y),\n');
fprintf('with source and receiver (blavk circles) and heterogenity (red circle)\n');
fprintf('superimposed.\n');
% gdama16_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 MatLab function
% circshift(), except that it allows delays that
% are non-integral multiples of the sampling interval
% Dt. It uses a Fourier technique.
clear all
fprintf('gdama16_03\n');
% 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 = xmin + Dx*[0:Nx-1]';
% independent variable z
zmin = 0;
zmax = 1;
Nz = 101;
Dz = (zmax-zmin)/(Nz-1);
z = zmin + Dz*[0:Nz-1]';
% independent variable t
% sampling Dt depends on sampling Dz
Dt = s0*Dz;
Nt = 512;
tmin = 0;
t = tmin + Dt*[0:Nt-1]';
tmax = t(end);
% 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*Dt;
t0 = 2*sigmat;
p = exp( -((t-t0).^2)/(2*sigmat*sigmat) );
pdd = [0;diff(p,2);0]; % first derivative
pdd = (0.1/(4*pi))*pdd/max(pdd); % second derivative
% unperturbed plane wave pulse on the grid
u0 = zeros(Nz,Nx,Nt);
for iz = (1:Nz)
tau0 = s0*(z(iz)-zmin); % phase delay for a plane wave
d=gda_delay(p,Dt,tau0);
for ix = (1:Nx)
% note that pwave wave experiences no geometerical spreading
u0(iz,ix,:)=d;
end
end
% array of stations (receivers), all at zmin
istas=(1:5:Nx)';
% point scatterer
xs = 0;
zs = 0.25;
% scattered wavefield
tau0 = s0*(zs-zmin); % delay of planewave down to scatterer
du = zeros(Nz,Nx,Nt);
for iz = (1:Nz)
for ix = (1:Nx)
R=sqrt(((x(ix)-xs)^2)+((z(iz)-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,:)=d*(4*pi*Ri);
end
end
% plot wavefield
figure();
clf;
subplot(2,2,1);
hold on;
itf = 20;
it = itf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [xmin, xmax, zmin, zmax] );
colormap('jet');
imagesc([xmin, xmax], [zmax, zmin], u0(:,:,it)+du(:,:,it), [-1, 1]);
plot(xs,zmax-zs,'wo','LineWidth',3);
plot(x(istas), zmax-0, 'k^', 'LineWidth', 2 );
subplot(2,2,2);
it = itf*2;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
axis( [xmin, xmax, zmin, zmax] );
colormap('jet');
imagesc([xmin, xmax], [zmax, zmin], u0(:,:,it)+du(:,:,it), [-1, 1]);
plot(xs,zmax-zs,'wo','LineWidth',3);
plot(x(istas), zmax-0, 'k^', 'LineWidth', 2 );
subplot(2,2,3);
it = itf*3;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
axis( [xmin, xmax, zmin, zmax] );
colormap('jet');
imagesc([xmin, xmax], [zmax, zmin], u0(:,:,it)+du(:,:,it), [-1, 1]);
plot(xs,zmax-zs,'wo','LineWidth',3);
plot(x(istas), zmax-0, 'k^', 'LineWidth', 2 );
subplot(2,2,4);
it = itf*4;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
axis( [xmin, xmax, zmin, zmax] );
colormap('jet');
imagesc([xmin, xmax], [zmax, zmin], u0(:,:,it)+du(:,:,it), [-1, 1]);
plot(xs,zmax-zs,'wo','LineWidth',3);
plot(x(istas), zmax-0, 'k^', 'LineWidth', 2 );
fprintf('Caption: Plot of wavefield\n');
% plot record section on surface
figure();
clf;
hold on;
itp = 101;
tp = Dt*(itp-1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [xmin, xmax, 0, tp] );
colormap('jet');
imagesc([xmin, xmax], [0, tp], squeeze(u0(1,:,1:itp)+du(1,:,1:itp))', [-1, 1]);
xlabel('x');
ylabel('t');
fprintf('Caption: record section on z=0 surface\n');
% back-propagate du
dubT = zeros(Nz,Nx,Nt);
for ixr = istas' % for each receiver on surface
xr = x(ixr); % receiver position
zr = z(1);
dub = zeros(Nz,Nx,Nt); % wavefield at receiver
a = squeeze(du(1,ixr,:)); % reduce to vector
for iz = (1:Nz) % back-propagate to each gridpoint
for ix = (1:Nx)
R=sqrt(((x(ix)-xr)^2)+((z(iz)-zr)^2)); % distance
Ri = 1/(R+epsi); % fudged geometrical spreading
taub = -s0*R; % time advance
d=gda_delay(a,Dt,taub);
dub(iz,ix,:)=d*(4*pi*Ri); % for this observation
end
end
dubT=dubT+dub; % stack all observations
end
% correlate
c = zeros(Nz,Nx);
for iz = (1:Nz) % each gridpoint
for ix = (1:Nx)
a = squeeze(u0(iz,ix,:)); % reduce to vector
add = [0;diff(a,2);0]; % 2nd derivative
b = squeeze(dubT(iz,ix,:)); % reduce to vector
c(iz,ix)=add'*b; % correlate (= dot product)
end
end
cmax = max(max(c));
% plot correlation
figure();
clf;
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
axis( [xmin, xmax, zmin, zmax] );
colormap('jet');
imagesc([xmin, xmax], [zmax, zmin], c/cmax, [-1, 1]);
plot(x(istas), zmax-0, 'k^', 'LineWidth', 2 );
xlabel('x');
ylabel('t');
fprintf('Caption: Plot of correlation (migrated record section)\n');
% gdama16_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.
clear all;
fprintf('gdama16_04\n');
% 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 = floor(5*N/16); % left limit of secondary pulse position
Nb = floor(11*N/16); % left limit of secondary pulse position
No2 = floor(N/2);
t = Dt*[0:N-1]';
t0 = t(No2);
ta = t(Na);
tb = t(Nb);
% reference signal u0 has just one main pulse
u0 = zeros(N,1);
u0(No2) = 1.0;
u0f = gda_chebyshevfilt(u0, Dt, flow, fhigh); % bandpass filter
ilist = [Na:10:Nb]'; % times at which delay is calculated
ilist2 = [Na:100:Nb]'; % times at which time series is plotted
Nlist = length(ilist);
tpulse = zeros(Nlist,1); % save calculation in array
tau = zeros(Nlist,1);
tau2 = zeros(Nlist,1);
% plot of time series
figure()
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [10, 35, -1, 10] );
title('u(t)');
xlabel('time t (s)');
j=0; % counts iterations
for i=ilist' % loop over pulse positions
% pulse position
j=j+1;
tpulse(j)=Dt*(i-No2);
% another signal u has two pulses
du = zeros(N,1);
du(i) = pulsesize;
duf = gda_chebyshevfilt(du, Dt, flow, fhigh); % band pass filtered perturbation
u = u0+du;
uf = u0f + duf; % bandpass filtered signal
if( ~isempty(find(ilist2==i,1)) )
plot( t, 100*uf+0.06*j, 'k-', 'LineWidth', 2 );
end
% 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) = tBmA_est;
cmax = Cmax_est;
% formula from Marquering et al, GJI 137, 805-815, 1999
u0fd = [diff(u0f); 0]/Dt; % 1st derivative
u0fdd = [0; diff(u0f,2); 0]/(Dt^2); % 2nd derivative
tau2(j) = (Dt*u0fd'*duf)/(Dt*u0fdd'*u0f);
end
fprintf('Caption: Series of time series with a reference pulse (large) all with the same start\n');
fprintf('time and smaller perturbation pulses with different start times\n');
% plot of delays
figure();
fprintf('amplitude %.3f lowpass %.3f highpass %.3f\n', pulsesize, flow, fhigh);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
hold on;
scale = 80;
axis( [ta-t0, tb-t0, (ta-t0)/scale, (tb-t0)/scale] );
% reference lines
plot( [0,0]', [(ta-t0)/scale, (tb-t0)/scale]', 'b:', 'LineWidth', 2 );
plot( [ta-t0, tb-t0]', [0,0]', 'b:', 'LineWidth', 2 );
xlabel('secondary pulse position (s)');
ylabel('perturbation in arrival time');
plot( tpulse, tau, 'k-', 'LineWidth', 3);
plot( tpulse, tau2, 'r-', 'LineWidth', 2);
fprintf('Caption: Delay associated with small perturbation pulse, calculated using\n');
fprintf('cross-correlation (black) and pertubative formula (red)\n');
% gdama16_05
% bannana-doughnut kernel for the acoustic case4
clear all;
fprintf('gdama16_05\n');
% images of kernel
figure(1);
clf;
% time series
figure(2);
clf;
% reference slowness
s0 = 1;
m = 1;
% independent variable x
Nx = 101;
Dx = 1;
x = Dx*[0:Nx-1]';
xmin = x(1);
xmax = x(end);
% independent variable y
Ny = 101;
Dy = 1;
y = Dy*[0:Ny-1]';
ymin = y(1);
ymax = y(end);
% independent variable t
SSF = 10;
Nt = 2*Nx*SSF; % number of times depends on number of x's
Dt = 1/SSF; % time sampling depends on spatial sampling
t = Dt*[0:Nt-1]';
tmin=t(1);
tmax=t(end);
% source time
t0 = tmin + (tmax-tmin)/2;
% define three frequency bands, one for each case
wny = (2*pi)/(2*Dt);
w0 = [0.005*wny; 0.025*wny; 0.1*wny];
% loop over three cases of finite frequencies
for iw=(1:3)
fprintf('Frequenvcy band %d: w0: %f\n', iw, w0(iw) );
% a band-limited source time function
fS = sinc(2*w0(iw)*(t-t0)) - sinc(0.5*w0(iw)*(t-t0));
% d/dt of source time function
fSd = zeros(Nt,1);
fSd(1:Nt-1) = (fS(2:Nt)-fS(1:Nt-1))/(Dt);
% d2/dt2 of source time function
fSdd = zeros(Nt,1);
fSdd(2:Nt-1) = (fS(1:Nt-2)-2*fS(2:Nt-1)+fS(3:Nt))/(Dt^2);
% source location
xS = 80;
yS = 51;
% receiver location
xR = 20;
yR = 51;
% Source to Receiver distance and time
RSR = sqrt( (xS-xR)^2 + (yS-yR)^2 );
TSR = s0*RSR;
nSR = floor( TSR/Dt );
% a constant used later
D = - (1/(4*pi*RSR)^2) * Dt * sum( fSd .^ 2 );
d = zeros(Nx,Ny); % the kernel on the (x,y) plane
for ix = (1:Nx) % loop over x
for iy = (1:Ny) % loop over y
% location of heterogeneity
xh = x(ix);
yh = y(iy);
% 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 = floor((TSh+ThR-TSR)/Dt);
c = 4*pi*RSR * 4*pi*(RSh+Dx) * 4*pi*(RhR+Dx);
% note, the Dx is a kludge to prevent singularies
% adjoint field
a = zeros(Nt,1);
a(1:Nt-NShR) = fSd(NShR+1:Nt);
% second derivative of reference field
b = fSdd;
% kernel is an inner product of a and b
d(ix,iy) = (2*s0*D)*Dt*sum( a .* b )/c;
end
end
% plot kernel
figure();
clf;
hold on;
colormap('jet');
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on
axis ij
axis( [xmin, xmax, ymin, ymax] );
dabsmax = max(abs(d(:)))/3;
imagesc(d, [-dabsmax, dabsmax]);
plot( yS, xS, 'ko', 'LineWidth', 4 );
plot( yR, xR, 'ko', 'LineWidth', 4 );
xlabel('y');
ylabel('x');
fprintf('Caption: Bannana-doughnut kernel for a frequency %d\n',iw);
% plot time series
figure();
clf;
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [tmin, tmax, -1.1*max(abs(fS)), 1.1*max(abs(fS))] );
plot( t, fS, 'k-', 'LineWidth', 2 );
xlabel('time, t');
ylabel('fS');
fprintf('Caption: Source wavelet for frequency %d\n',iw);
fprintf('\n');
fprintf('\n');
fprintf('\n');
end
