地震学和地球结构（震源追踪）

2016年1月28日 3162点热度 1人点赞 0条评论

C-2. (a) 编写程序：震源在地表，选择震源一定范围的入射角，追踪穿过一个地球模型的射线。图形的输出显示震源、射线路径、地球表面、核幔边界以及内核-外核边界。

(b) 用PREM模型追踪震源地表的地震射线，表示出上地幔间断面和地核对射线的影响。

一、实验原理

运用地球的射线参数公式求取射线路径的最低点。

运用圆的几何角度关系，和角度公式计算每层的入射点

走时计算公式

读入prem模型的速度和深度数据，计算出一组数据，运用matlab绘图。

二、程序源代码

射线追踪程序

s=xlsread('perm.xlsx');

r=s(:,1);

v=s(:,2);

for i=1:26

zhong=r(i,1);r(i,1)=r(54-i,1);r(54-i,1)=zhong;

zhong=v(i,1);v(i,1)=v(54-i,1);v(54-i,1)=zhong;

end

for jiao=0:0.01:0.1

% jiao=input('请输入入射角');%入射角

p=r(1)*sin(jiao)./v(1);

N=size(r);

for i=2:N

k=v(i)*p./r(i);

if k>=1;

if k==1

M=i

Q=0;break;

else if r(i)==r(i-1)

M=i-1

R=r(M,1);

V=v(M,1);

y=R/V;

Q=p/(R*(y^2-p^2))

break;

else Xl=(v(i)-v(i-1))./(r(i)-r(i-1));

r(i)=(Xl*r(i-1)-v(i-1))./(Xl-1/p);

v(i)=r(i)/p;

M=i

Q=0

break;

end

else continue;

end

r(M,1)=int64(r(M,1));

t=0;

T=0;

zb1=zeros(1,6371);

zb2=zeros(1,6371);

for i=M:-1:2

K=r(i-1)-r(i);

for j=0:K

R=r(i)+j;

V=v(i)+j*(v(i-1)-v(i))./(r(i-1)-r(i));

y=R/V;

deta=p/(R*(y^2-p^2)^0.5);

T=deta+T;

zb1(j+t+1)=R*sin(T);

zb2(j+t+1)=R*cos(T);

end

t=K+1+t;

end

for i=1:6370

if zb1(i)==0 break

end

zb=zeros(2,i-1);

phi=i-1;

for j=1:i-1

zb(1,j)=zb1(j);

zb(2,j)=zb2(j);

end

zb3=zeros(2,2*phi);

for j=1:phi

zb3(1,phi+j)=zb(1,j);

zb3(1,phi+1-j)=-zb(1,j);

zb3(2,phi+j)=zb(2,j);

zb3(2,phi+1-j)=zb(2,j);

end

plot(zb3(1,:),zb3(2,:));hold on;

len=[6371 3630 1221.5];%边界

x=-6371:6371;

for i=1:3

y=sqrt(len(i).^2-x.^2);

plot(x,y);hold on;

end

走时曲线程序

s=xlsread('perm.xlsx');

r=s(:,1);

v=s(:,2);

Sj=zeros(2,50);

for i=1:26

zhong=r(i,1);r(i,1)=r(54-i,1);r(54-i,1)=zhong;

zhong=v(i,1);v(i,1)=v(54-i,1);v(54-i,1)=zhong;

end

xishur=r;

xishuv=v;

for g=1:50

r=xishur;

v=xishuv;

jiao=(51-g)*0.02;

p=6371*sin(jiao)/1450;

N=size(r);

for i=2:N

k=v(i)*p./r(i);

if k>=1;

if k==1

M=i;

Q=0;break;

else if r(i)==r(i-1)

M=i-1;

R=r(M,1);

V=v(M,1);

y=R/V;

Q=p/(R*(y^2-p^2));

break;

else Xl=(v(i)-v(i-1))./(r(i)-r(i-1));

r(i)=(Xl*r(i-1)-v(i-1))./(Xl-1/p);

v(i)=r(i)/p;

M=i;

Q=0;

break;

end

else continue;

end

r(M,1)=int64(r(M,1));

T=0;

D=0;

for i=M:-1:2

K=r(i-1)-r(i)

for j=0:K

R=r(i)+j;

V=v(i)+j*(v(i-1)-v(i))./(r(i-1)-r(i));

y=R/V;

deta=p/(R*(y^2-p^2)^0.5);

time=y^2/(R*(y^2-p^2)^0.5);

D=deta+D;

T=time+T;

end

Sj(1,g)=D/pi*180;

Sj(2,g)=T;

end

plot(Sj(1,:),Sj(2,:));

结果及分析

结果分析，此图在入射角为0.04的情况下得出。因为速度是连续变化的，所以绘制出图线为曲线。因为上地幔间断面为速度间断面，所以使射线弯折。入射角越小，地震波穿越地球越深。当入射角为0度的时候，地震波射线直接穿过地核，到达对称位置，入射角越大，地震波穿越深度越浅，深度很小时就到达最低点后返回。上地幔间断面和核幔边界因为速度突然减小，所以图线会下凹。核幔边界则会偏向地核。在内外核边界速度增大，本身图线应该下凹，但是下凹程度减小，近似为直线。走时曲线也为曲线，在弯折处为速度突变面。所以第一个弯折出为上地幔间断面。

地震学和地球结构（震源追踪）

文章评论