斯奈尔定律和Zoeppritz方程

1、斯奈尔定律和Zoeppritz方程姓名： 学号： 专业：地球物理勘察技术2012级一、实验目的1. 利用密度、上下界面的纵横波速度通过斯奈尔定律求出该界面的0°90°入射角的反射角度和透射角度。2. 利用Zoeppritz方程绘制反射系数和透射系数曲线。二、实验步骤1、模型：Vp1=3300 m/s Vs1=1585 m/s 1 =2.4 g/cm3 Vp1=3100 m/s Vs1=1989 m/s 1 =2.24 g/cm3 计算从090度的反射透射系数曲线 如图1-1-4三、实验结果P波入射的反射透射角度：图1P波入射的反射透射系数曲线：图2SV波入射的反射透射角度正

2、弦值：图3SV波入射的反射透射系数曲线：图4SH波入射的反射透射角度正弦值：图5SH波入射的反射透射系数曲线：图6四、实验分析因为各波反射和透射波振幅系数与其能量成正比，由此可以看出其能量的变化以及入射波能量的分配。当上层介质为密介质，下层介质为疏介质时：由图2知，P波从上层密介质入射到界面时：随着入射角增大，P波反射系数、P波透射系数在减小，即随着入射角的增大，反射P波、透射P波的能量在减小，而S波反射系数、S波透射系数在增大，即随着入射角的增大，反射S波、透射S波的能量在增大。说明当入射角发生变化时，入射波的能量分配在改变。由图4可知，SV波从上层密介质入射到界面时：入射角小于29度时，随

3、着入射角增大，SV波反射系数在减小，P波反射系数、SV波透射系数、P波透射系数在增大，SV波透射系数最大，即反射SV波能量在减小，反射P波、透射SV波、透射P波能量在增大，透射SV波的能量最大；入射角大于等于29度且小于31度时，SV波反射系数、SV波透射系数减小，P波反射系数、P波透射系数在增大，SV波透射系数最大，即随着入射角的增大，反射SV波、透射SV波能量在减小，反射P波、透射P波能量在增大，透射SV波能量最大；入射角大于等于31度且小于53度时，SV波透射系数增大，SV波反射系数、P波反射系数、P波透射系数先减小后增大，SV波透射系数仍最大，即随着入射角的增大，透射SV波能量增大，反

4、射SV波、反射P波、透射P波能量先减小后增大，透射SV波能量最大；入射角大于53度时， SV波反射系数为1，SV波透射系数、P波反射系数、P波透射系数减小，即随着入射角的增大，反射SV波能量不变，透射SV波、反射P波、透射P波能量减小。由图6可知，SH波从上层密介质入射到界面时：不产生转换波，入射角小于临界角时，SH波的反射系数、透射系数均随着入射角的增大而增大，即反射、透射SH波能量增大；入射角大于等于临界角时，随着入射角的增大，SH波的反射系数几乎不变，透射系数减小，即反射SH波的能量减小，透射SH波的能量几乎不变。五、附：源程序代码 P波入射时，程序：#include<stdio.

5、h>#include<math.h>#include"6GAUS.C"#define PI 3.1415926void main()FILE *fp1,*fp2;int i,n91;double ipp,x1,x2,y1,y2,pr,sr,pt,st,a195,a295,b195,b295,vp1=3300, vp2=3100,vs1=1585,vs2=1989,den1=2.4,den2=2.24,k=den2/den1;static double a44=0.0,b4=0.0;fp1=fopen("snell.csv","

6、w");fp2=fopen("P波入射反射透射系数.csv","w");for(i=0;i<=90;i+)ni=i;ipp=i*PI/180;x1=ipp;a1i=x1*180/PI;x2=sin(ipp)*vs1/vp1;a2i=asin(x2)*180/PI;y1=sin(ipp)*vp2/vp1;b1i=asin(y1)*180/PI;y2=sin(ipp)*vs2/vp1;b2i=asin(y2)*180/PI;fprintf(fp1,"%d,%f,%f,%f,%fn",ni,a1i,a2i,b1i,b2i);

7、/*输出P波入射反射投射角度*/ pr=a1i*PI/180;sr=a2i*PI/180;pt=b1i*PI/180;st=b2i*PI/180;a00=sin(pr);a01=cos(sr);a02=-sin(pt);a03=-cos(st); a10=cos(pr);a11=-sin(sr);a12=cos(pt);a13=-sin(st);a20=sin(2*pr);a21=cos(2*sr)*vp1/vs1;a22=sin(2*pt)*k*vp1*pow(vs2/vs1,2)/vp2;a23=cos(2*st)*k*vp1*vs2/pow(vs1,2);a30=cos(2*sr);a3

8、1=-sin(2*pr)*vs1/vp1;a32=-cos(2*st)*k*vp2/vp1;a33=sin(2*pt)*k*vs2/vp1;/*输入系数矩阵*/b0=-sin(pr);b1=cos(pr);b2=sin(2*pr);b3=-cos(2*sr);/*输入常数矩阵*/if(gaus(a,b,4)!=0)fprintf(fp2,"%d,%f,%f,%f,%fn",ni,b0,b1,b2,b3);fclose(fp1);fclose(fp2);SV波入射时，程序：#include<stdio.h>#include<math.h>#includ

9、e"4CINV.C"#include"4TCMUL.C"#define PI 3.1415926void main()FILE *fp1,*fp2;int i,n91;double in,ipp1,ipp2,ipp3, x1,x2,y1,y2, pr,sr,pt,st,p, rsp,rss,tsp,tss, a191,a291,b191,b291, vp1=3300,vp2=3100,vs1=1585,vs2=1989, den1=2.4,den2=2.24,k=den2/den1; static double ar44,ai44,br4,bi4,cr4

10、1,ci41; fp1=fopen("snell.csv","w");fp2=fopen("SV波入射反射透射系数.csv","w");ipp1=asin(vs1/vp1);ipp2=asin(vs1/vs2);ipp3=asin(vs1/vp2);/*临界角*/for(i=0;i<=90;i+)ni=i;in=i*PI/180;p=sin(in)/vs1;x1=p*vs1;x2=p*vp1;y1=p*vs2;y2=p*vp2;/*snell定律a1i=x1;a2i=x2;b1i=y1;b2i=y2;fpri

11、ntf(fp1,"%d,%f,%f,%f,%fn",ni,a1i,a2i,b1i,b2i);/*输出S波入射反射透射角度正弦值*/ if(in<ipp1) sr=asin(x1); pr=asin(x2); st=asin(y1); pt=asin(y2); ar00=-sin(pr); ar01=cos(sr); ar02=sin(pt); ar03=cos(st); ar10=-cos(pr); ar11=-sin(sr); ar12=-cos(pt); ar13=sin(st); ar20=2*cos(pr)*sin(pr)*vs1/vp1; ar21=-(1-

12、2*sin(sr)*sin(sr); ar22=2*cos(pt)*sin(pt)*k*vs2*vs2/(vs1*vp1); ar23=(1-2*sin(st)*sin(st)*k*vs2/vs1; ar30=(1-2*sin(sr)*sin(sr)*vp1/vs1; ar31=2*cos(sr)*sin(sr); ar32=-(1-2*sin(st)*sin(st)*k*vp2/vp1; ar33=2*cos(st)*sin(st)*k*vs2/vp1; /*输入系数矩阵 br0=cos(sr); br1=sin(sr); br2=1-2*sin(sr)*sin(sr); br3=2*cos

13、(sr)*sin(sr); if(cinv(ar,ai,4)!=0)/*求系数矩阵的逆矩阵*/tcmul(ar,ai,br,bi,4,4,1,cr,ci);rsp=sqrt(cr00*cr00+ci00*ci00); rss=sqrt(cr10*cr10+ci10*ci10);tsp=sqrt(cr20*cr20+ci20*ci20);tss=sqrt(cr30*cr30+ci30*ci30);/*求S波入射反射透射系数*/fprintf(fp2,"%d,%f,%f,%f,%fn",ni,rsp,rss,tsp,tss);/*输出S波入射反射透射系数*/else if(in

14、>=ipp1&&in<ipp3) sr=asin(x1); st=asin(y1); pt=asin(y2); ar00=-p*vp1,ai00=0; ar01=cos(sr),ai01=0; ar02=sin(pt),ai02=0; ar03=cos(st),ai03=0; ar10=0,ai10=-sqrt(p*p*vp1*vp1-1); ar11=-sin(sr),ai11=0; ar12=-cos(pt),ai12=0; ar13=sin(st),ai13=0; ar20=0,ai20=2*p*vp1*sqrt(p*p*vp1*vp1-1)*vs1/vp1;

15、 ar21=-(1-2*sin(sr)*sin(sr),ai21=0; ar22=2*cos(pt)*sin(pt)*k*vs2*vs2/(vs1*vp1),ai22=0; ar23=(1-2*sin(st)*sin(st)*k*vs2/vs1,ai23=0; ar30=(1-2*sin(sr)*sin(sr)*vp1/vs1,ai30=0; ar31=2*cos(sr)*sin(sr),ai31=0; ar32=-(1-2*sin(st)*sin(st)*k*vp2/vp1,ai32=0; ar33=2*cos(st)*sin(st)*k*vs2/vp1,ai33=0; /*输入系数矩阵 b

16、r0=cos(sr),bi0=0; br1=sin(sr),bi1=0; br2=1-2*sin(sr)*sin(sr),bi2=0; br3=2*cos(sr)*sin(sr),bi3=0; /*输入常数矩阵if(cinv(ar,ai,4)!=0)/*求系数矩阵的逆矩阵*/tcmul(ar,ai,br,bi,4,4,1,cr,ci);rsp=sqrt(cr00*cr00+ci00*ci00); rss=sqrt(cr10*cr10+ci10*ci10);tsp=sqrt(cr20*cr20+ci20*ci20);tss=sqrt(cr30*cr30+ci30*ci30);/*求S波入射反射透

17、射系数*/fprintf(fp2,"%d,%f,%f,%f,%fn",ni,rsp,rss,tsp,tss);/*输出S波入射反射透射系数*/else if(in>=ipp3&&in<ipp2) sr=asin(x1); st=asin(y1); ar00=-p*vp1,ai00=0; ar01=cos(sr),ai01=0; ar02=p*vp2,ai02=0; ar03=cos(st),ai03=0; ar10=0,ai10=-sqrt(p*p*vp1*vp1-1); ar11=-sin(sr),ai11=0; ar12=0,ai12=-sq

18、rt(p*p*vp2*vp2-1); ar13=sin(st),ai13=0; ar20=0,ai20=2*p*vp1*sqrt(p*p*vp1*vp1-1)*vs1/vp1; ar21=-(1-2*sin(sr)*sin(sr),ai21=0; ar22=0,ai22=2*p*vp2*sqrt(p*p*vp2*vp2-1)*k*vs2*vs2/(vs1*vp1); ar23=(1-2*sin(st)*sin(st)*k*vs2/vs1,ai23=0; ar30=(1-2*sin(sr)*sin(sr)*vp1/vs1,ai30=0; ar31=2*cos(sr)*sin(sr),ai31=0

19、; ar32=-(1-2*sin(st)*sin(st)*k*vp2/vp1,ai32=0; ar33=2*cos(st)*sin(st)*k*vs2/vp1,ai33=0; /*输入系数矩阵 br0=cos(sr),bi0=0; br1=sin(sr),bi1=0; br2=1-2*sin(sr)*sin(sr),bi2=0; br3=2*cos(sr)*sin(sr),bi3=0; /*输入常数矩阵if(cinv(ar,ai,4)!=0)/*求系数矩阵的逆矩阵tcmul(ar,ai,br,bi,4,4,1,cr,ci);rsp=sqrt(cr00*cr00+ci00*ci00); rss=

20、sqrt(cr10*cr10+ci10*ci10);tsp=sqrt(cr20*cr20+ci20*ci20);tss=sqrt(cr30*cr30+ci30*ci30);/*求S波入射反射透射系数*/fprintf(fp2,"%d,%f,%f,%f,%fn",ni,rsp,rss,tsp,tss);/*输出S波入射反射透射系数*/else if(in>=ipp2) sr=asin(x1); ar00=-p*vp1,ai00=0; ar01=cos(sr),ai01=0; ar02=p*vp2,ai02=0; ar03=0,ai03=sqrt(p*p*vs2*vs2-

21、1); ar10=0,ai10=-sqrt(p*p*vp1*vp1-1); ar11=-sin(sr),ai11=0; ar12=0,ai12=-sqrt(p*p*vp2*vp2-1); ar13=p*vp2,ai13=0; ar20=0,ai20=2*p*vp1*sqrt(p*p*vp1*vp1-1)*vs1/vp1; ar21=-(1-2*sin(sr)*sin(sr),ai21=0; ar22=0,ai22=2*p*vp2*sqrt(p*p*vp2*vp2-1)*k*vs2*vs2/(vs1*vp1); ar23=(1-2*p*p*vs2*vs2)*k*vs2/vs1,ai23=0; a

22、r30=(1-2*sin(sr)*sin(sr)*vp1/vs1,ai30=0; ar31=2*cos(sr)*sin(sr),ai31=0; ar32=-(1-2*p*p*vs2*vs2)*k*vp2/vp1,ai32=0; ar33=0,ai33=2*p*vs2*sqrt(p*p*vs2*vs2-1)*k*vs2/vp1; /*输入系数矩阵 br0=cos(sr),bi0=0; br1=sin(sr),bi1=0; br2=1-2*sin(sr)*sin(sr),bi2=0; br3=2*cos(sr)*sin(sr),bi3=0; /*输入常数矩阵if(cinv(ar,ai,4)!=0)

23、/*求系数矩阵的逆矩阵tcmul(ar,ai,br,bi,4,4,1,cr,ci);rsp=sqrt(cr00*cr00+ci00*ci00); rss=sqrt(cr10*cr10+ci10*ci10);tsp=sqrt(cr20*cr20+ci20*ci20);tss=sqrt(cr30*cr30+ci30*ci30);/*求S波入射反射透射系数*/fprintf(fp2,"%d,%f,%f,%f,%fn",ni,rsp,rss,tsp,tss);/*输出S波入射反射透射系数*/fclose(fp1);fclose(fp2);SH波入射时，程序：#include<stdio.h>#include<math.h>#define PI 3.1415926void main()FILE *fp1,*fp2;int i,n91;double in,ipp,x,y,st,rsr,rsi,tsr,tsi,rsh,tsh,p,q,r91,t91,vs1=1585,v