上述实验指导书的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 change 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')

转载本文请联系原作者获取授权，同时请注明本文来自毛宁波科学网博客。
链接地址：https://blog.sciencenet.cn/blog-339326-422846.html

上一篇：毛宁波地震勘探原理实验参考---MATLAB命令大全
下一篇：毛宁波美国麻省理工学院访学见闻(1)---访问基本情况