上述实验指导书PDF文件下载：地震勘探原理 实验1
以下是P_and_S_waves.m的matlab程序
function P_and_S_waves
clear all
close all
f = 2;
T = 1/f;
omega = 2*pi*f;
vel = 3000;
%xg =[-1.5,1.5]
%zg =[0,0]
% Set up array of depths
x1 = [0:150:5000];
x2 = x1;
nx = length(x1);
k = 2*pi*f/vel;
%Start time step loop
n=40;
M=moviein(n);
for it=1:n
t=T*(it)/n;
for ix=1:nx
A1(ix) = 0.2*cos(k*x1(ix)-omega*t);
A2(ix) = 0.0;
x2(ix) = x1(ix) 100*cos(k*x1(ix)-omega*t);
end
% Plot S-wave
plot(x1*0.001,A1 1.0,'o');
hold on
plot(x1*0.001,A1 1.3,'o');
plot(x1*0.001,A1 1.15,'o');
plot(x1*0.001,A1 0.85,'o');
plot(x1*0.001,A1 0.7,'o');
plot(x1(20)*0.001,A1(20) 1.30,'ro')
plot(x1(20)*0.001,A1(20) 1.15,'ro')
plot(x1(20)*0.001,A1(20) 1,'ro')
plot(x1(20)*0.001,A1(20) 0.85,'ro')
plot(x1(20)*0.001,A1(20) 0.70,'ro')
plot(x1(1)*0.001,A1(1) 1.30,'ro')
plot(x1(1)*0.001,A1(1) 1.15,'ro')
plot(x1(1)*0.001,A1(1) 1,'ro')
plot(x1(1)*0.001,A1(1) 0.85,'ro')
plot(x1(1)*0.001,A1(1) 0.70,'ro')
% Plot P-wave
plot(x2*0.001,A2-1.3,'o');
plot(x2*0.001,A2-1.15,'o');
plot(x2*0.001,A2-1.0,'o');
plot(x2*0.001,A2-0.85,'o');
plot(x2*0.001,A2-0.7,'o');
plot(x2(20)*0.001,A2(20)-1.30,'ro');
plot(x2(20)*0.001,A2(20)-1.15,'ro');
plot(x2(20)*0.001,A2(20)-1,'ro');
plot(x2(20)*0.001,A2(20)-0.85,'ro');
plot(x2(20)*0.001,A2(20)-0.70,'ro');
plot(x2(1)*0.001,A2(1)-1.30,'ro');
plot(x2(1)*0.001,A2(1)-1.15,'ro');
plot(x2(1)*0.001,A2(1)-1,'ro');
plot(x2(1)*0.001,A2(1)-0.85,'ro');
plot(x2(1)*0.001,A2(1)-0.7,'ro');
hold off
axis([-0.5,5,-2.2,2.2])
title('P-waves and S-waves')
ylabel ('A')
xlabel ('(km)')
M(:,it)=getframe;
end
movie(M,10)
print -dpsc P-wave-and-S-wave.ps
以下是rayleigh_waves.m的matlab程序
function rayleigh_waves
clear all
close all
f = 0.5;
T = 1/f;
omega = 2*pi*f;
vel = 3000;
ncol = 7;
% Set up array of depths
x = [0:300:15000];
z = [0:-100:-2000];
Ax = 200;
Az = 200;
nx = length(x);
nz = length(z);
k = 2*pi*f/vel;
lamda = 750;
%Start time step loop
n=40;
M=moviein(n);
for it=1:n
t=T*(it)/n;
for ix=1:nx
dx(ix) = 0.001*(x(ix) Ax*exp(z(1)/lamda)*sin(k*x(ix)-omega*t));
dz(ix) = 0.001*(z(1) Az*exp(z(1)/lamda)*cos(k*x(ix)-omega*t));
end
plot(dx,dz,'o-')
hold on
plot(dx(20),dz(20),'or')
axis([-0.2,15.2,-2.05,1])
title('Rayleigh wave f=0.5 Hz')
ylabel ('depth(km)')
xlabel ('(km)')
for iz=2:ncol-1
for ix=1:nx
dx(ix) = 0.001*(x(ix) Ax*exp(z(iz)/lamda)*sin(k*x(ix)-omega*t));
dz(ix) = 0.001*(z(iz) Az*exp(z(iz)/lamda)*cos(k*x(ix)-omega*t));
end
plot(dx,dz,'o')
end
for ix=1:nx
dx(ix) = 0.001*(x(ix) Ax*exp(z(ncol)/lamda)*sin(k*x(ix)-omega*t));
dz(ix) = 0.001*(z(ncol) Az*exp(z(ncol)/lamda)*cos(k*x(ix)-omega*t));
end
plot(dx,dz,'ro')
for iz=ncol 1:nz
for ix=1:nx
dx(ix) = 0.001*(x(ix) Ax*exp(z(iz)/lamda)*sin(k*x(ix)-omega*t));
dz(ix) = 0.001*(z(iz) Az*exp(z(iz)/lamda)*cos(k*x(ix)-omega*t));
end
plot(dx,dz,'o')
end
hold off
M(:,it)=getframe;
end
movie(M,10)
print -dpsc rayleigh-wave.ps
以下是synseismic.m的matlab程序
dt=0.001; % sampling rate in seconds
freq=15; % peak frequency of wave in Hertz
%%Do not cange the following
tracelength=1000; % length of trace
lengthwave=round(1.5/freq/dt); % length of wave
temp=[0:1:lengthwave]; % temporary variable
timeaxis=temp*dt; % time axis
% wave=-sin(temp*dt*2*pi*freq).*hanning(lengthwave+1)'; % this is the wave
wave=-sin(temp*dt*2*pi*freq)
tracelength=max(lengthwave,tracelength); % error checking
tracewave=zeros(1,tracelength);
tracewave(1,1:lengthwave+1)=wave;
freqaxis=[0:length(tracewave)-1]/length(tracewave)/dt; % frequency axis
fftwave=abs(fft(tracewave)); % fourier transform of wave
%% The following plots the wave and the fourier transform of the wave
figure
subplot(2,1,1)
plot(timeaxis,wave)
title('Input wave')
xlabel('Time (sec)')
ylabel('Amplitude')
grid
subplot(2,1,2)
plot(freqaxis(1:round(length(tracewave)/2)),fftwave(1:round(length(tracewave)/2)))
title('Fourier transform of input wave')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
grid
%% The following calculates the reflection series and seismic trace
%% Do not change the following
refserlength=1001; % length of reflection series
refsertime=[0:1:refserlength-1]*dt; % time axis of reflection series
refser=zeros(1,refserlength); % reflection series
refser(1,150)=0.8;
refser(1,175)=0.5;
refser(1,300)=0.25;
refser(1,333)=-0.75;
refser(1,500)=0.69;
refser(1,600)=0.2;
refser(1,750)=-0.33;
refser(1,833)=0.45;
refserfft=abs(fft(refser)); % fourier transform of reflection series
freqaxisrefser=[0:refserlength-1]/refserlength/dt; % frequency axis
seismictracetemp=conv(wave,refser); % convolution of wave with reflection series
seismictrace=seismictracetemp((round(lengthwave/2)+1):(round(lengthwave/2)+refserlength));
seismictracefft=abs(fft(seismictrace)); % fourier transform of seismic trace
figure
subplot(2,1,1)
plot(refsertime,refser)
title('Reflection series')
xlabel('Time (sec)')
ylabel('Amplitude')
subplot(2,1,2)
plot(freqaxisrefser(1:round(refserlength/2)),refserfft(1:round(refserlength/2)))
title('Frequency spectrum of reflection series')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
figure
subplot(2,1,1)
plot(refsertime,seismictrace)
title('Seismic trace')
xlabel('Time (sec)')
ylabel('Amplitude')
subplot(2,1,2)
plot(freqaxisrefser(1:round(refserlength/2)),seismictracefft(1:round(refserlength/2)))
title('Frequency spectrum of seismic trace')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
转载本文请联系原作者获取授权，同时请注明本文来自毛宁波科学网博客。
链接地址：http://blog.sciencenet.cn/blog-339326-422846.html
上一篇：毛宁波地震勘探原理实验参考---MATLAB命令大全
下一篇：毛宁波美国麻省理工学院访学见闻(1)---访问基本情况