斯奈尔定律和Zoeppritz方程.docx
《斯奈尔定律和Zoeppritz方程.docx》由会员分享，可在线阅读，更多相关《斯奈尔定律和Zoeppritz方程.docx（19页珍藏版）》请在冰豆网上搜索。
斯奈尔定律和Zoeppritz方程
斯奈尔定律和Zoeppritz方程
姓名：
学号：
专业：
地球物理勘察技术2012级
一、实验目的
1.利用密度、上下界面的纵横波速度通过斯奈尔定律求出该界面的0°~90°入射角的反射角度和透射角度。
2.利用Zoeppritz方程绘制反射系数和透射系数曲线。
二、实验步骤
1、模型：
Vp1=3300m/sVs1=1585m/sρ1=2.4g/cm3
Vp1=3100m/sVs1=1989m/sρ1=2.24g/cm3
计算从0～90度的反射透射系数曲线如图1-1-4
三、实验结果
P波入射的反射透射角度：
图1
P波入射的反射透射系数曲线：
图2
SV波入射的反射透射角度正弦值：
图3
SV波入射的反射透射系数曲线：
图4
SH波入射的反射透射角度正弦值：
图5
SH波入射的反射透射系数曲线：
图6
四、实验分析
因为各波反射和透射波振幅系数与其能量成正比，由此可以看出其能量的变化以及入射波能量的分配。
当上层介质为密介质，下层介质为疏介质时：
由图2知，P波从上层密介质入射到界面时：
随着入射角增大，P波反射系数、P波透射系数在减小，即随着入射角的增大，反射P波、透射P波的能量在减小，而S波反射系数、S波透射系数在增大，即随着入射角的增大，反射S波、透射S波的能量在增大。
说明当入射角发生变化时，入射波的能量分配在改变。
由图4可知，SV波从上层密介质入射到界面时：
入射角小于29度时，随着入射角增大，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波能量增大，反射SV波、反射P波、透射P波能量先减小后增大，透射SV波能量最大；
入射角大于53度时，SV波反射系数为1，SV波透射系数、P波反射系数、P波透射系数减小，即随着入射角的增大，反射SV波能量不变，透射SV波、反射P波、透射P波能量减小。
由图6可知，SH波从上层密介质入射到界面时：
不产生转换波，入射角小于临界角时，SH波的反射系数、透射系数均随着入射角的增大而增大，即反射、透射SH波能量增大；入射角大于等于临界角时，随着入射角的增大，SH波的反射系数几乎不变，透射系数减小，即反射SH波的能量减小，透射SH波的能量几乎不变。
五、附：
源程序代码
P波入射时，程序：
#include
#include
#include"6GAUS.C"
#definePI3.1415926
voidmain（）
{
FILE*fp1,*fp2;
inti,n[91];
doubleipp,x1,x2,y1,y2,pr,sr,pt,st,a1[95],a2[95],b1[95],b2[95],vp1=3300,vp2=3100,vs1=1585,vs2=1989,den1=2.4,den2=2.24,k=den2/den1;
staticdoublea[4][4]={0.0},b[4]={0.0};
fp1=fopen（"snell.csv","w"）;
fp2=fopen（"P波入射反射透射系数.csv","w"）;
for（i=0;i<=90;i++）
{
n[i]=i;
ipp=i*PI/180;
x1=ipp;
a1[i]=x1*180/PI;
x2=sin（ipp）*vs1/vp1;
a2[i]=asin（x2）*180/PI;
y1=sin（ipp）*vp2/vp1;
b1[i]=asin（y1）*180/PI;
y2=sin（ipp）*vs2/vp1;
b2[i]=asin（y2）*180/PI;
fprintf（fp1,"%d,%f,%f,%f,%f\n",n[i],a1[i],a2[i],b1[i],b2[i]）;
//*输出P波入射反射投射角度*//
pr=a1[i]*PI/180;
sr=a2[i]*PI/180;
pt=b1[i]*PI/180;
st=b2[i]*PI/180;
a[0][0]=sin（pr）;
a[0][1]=cos（sr）;
a[0][2]=-sin（pt）;
a[0][3]=-cos（st）;
a[1][0]=cos（pr）;
a[1][1]=-sin（sr）;
a[1][2]=cos（pt）;
a[1][3]=-sin（st）;
a[2][0]=sin（2*pr）;
a[2][1]=cos（2*sr）*vp1/vs1;
a[2][2]=sin（2*pt）*k*vp1*pow（vs2/vs1,2）/vp2;
a[2][3]=cos（2*st）*k*vp1*vs2/pow（vs1,2）;
a[3][0]=cos（2*sr）;
a[3][1]=-sin（2*pr）*vs1/vp1;
a[3][2]=-cos（2*st）*k*vp2/vp1;
a[3][3]=sin（2*pt）*k*vs2/vp1;
//*输入系数矩阵*//
b[0]=-sin（pr）;
b[1]=cos（pr）;
b[2]=sin（2*pr）;
b[3]=-cos（2*sr）;
//*输入常数矩阵*//
if（gaus（a,b,4）!
=0）
fprintf（fp2,"%d,%f,%f,%f,%f\n",n[i],b[0],b[1],b[2],b[3]）;
}
fclose（fp1）;
fclose（fp2）;
}
SV波入射时，程序：
#include
#include
#include"4CINV.C"
#include"4TCMUL.C"
#definePI3.1415926
voidmain（）
{
FILE*fp1,*fp2;
inti,n[91];
doublein,ipp1,ipp2,ipp3,
x1,x2,y1,y2,
pr,sr,pt,st,p,
rsp,rss,tsp,tss,
a1[91],a2[91],b1[91],b2[91],
vp1=3300,vp2=3100,vs1=1585,vs2=1989,
den1=2.4,den2=2.24,k=den2/den1;
staticdoublear[4][4],ai[4][4],br[4],bi[4],cr[4][1],ci[4][1];
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++）
{
n[i]=i;
in=i*PI/180;
p=sin（in）/vs1;
x1=p*vs1;
x2=p*vp1;
y1=p*vs2;
y2=p*vp2;//*snell定律
a1[i]=x1;
a2[i]=x2;
b1[i]=y1;
b2[i]=y2;
fprintf（fp1,"%d,%f,%f,%f,%f\n",n[i],a1[i],a2[i],b1[i],b2[i]）;
//*输出S波入射反射透射角度正弦值*//
if（in{
sr=asin（x1）;
pr=asin（x2）;
st=asin（y1）;
pt=asin（y2）;
ar[0][0]=-sin（pr）;
ar[0][1]=cos（sr）;
ar[0][2]=sin（pt）;
ar[0][3]=cos（st）;
ar[1][0]=-cos（pr）;
ar[1][1]=-sin（sr）;
ar[1][2]=-cos（pt）;
ar[1][3]=sin（st）;
ar[2][0]=2*cos（pr）*sin（pr）*vs1/vp1;
ar[2][1]=-（1-2*sin（sr）*sin（sr））;
ar[2][2]=2*cos（pt）*sin（pt）*k*vs2*vs2/（vs1*vp1）;
ar[2][3]=（1-2*sin（st）*sin（st））*k*vs2/vs1;
ar[3][0]=（1-2*sin（sr）*sin（sr））*vp1/vs1;
ar[3][1]=2*cos（sr）*sin（sr）;
ar[3][2]=-（1-2*sin（st）*sin（st））*k*vp2/vp1;
ar[3][3]=2*cos（st）*sin（st）*k*vs2/vp1;
//*输入系数矩阵
br[0]=cos（sr）;
br[1]=sin（sr）;
br[2]=1-2*sin（sr）*sin（sr）;
br[3]=2*cos（sr）*sin（sr）;
if（cinv（ar,ai,4）!
=0）//*求系数矩阵的逆矩阵*//
{
tcmul（ar,ai,br,bi,4,4,1,cr,ci）;
rsp=sqrt（cr[0][0]*cr[0][0]+ci[0][0]*ci[0][0]）;
rss=sqrt（cr[1][0]*cr[1][0]+ci[1][0]*ci[1][0]）;
tsp=sqrt（cr[2][0]*cr[2][0]+ci[2][0]*ci[2][0]）;
tss=sqrt（cr[3][0]*cr[3][0]+ci[3][0]*ci[3][0]）;//*求S波入射反射透射系数*//
fprintf（fp2,"%d,%f,%f,%f,%f\n",n[i],rsp,rss,tsp,tss）;//*输出S波入射反射透射系数*//
}
}
elseif（in>=ipp1&&in{
sr=asin（x1）;
st=asin（y1）;
pt=asin（y2）;
ar[0][0]=-p*vp1,ai[0][0]=0;
ar[0][1]=cos（sr）,ai[0][1]=0;
ar[0][2]=sin（pt）,ai[0][2]=0;
ar[0][3]=cos（st）,ai[0][3]=0;
ar[1][0]=0,ai[1][0]=-sqrt（p*p*vp1*vp1-1）;
ar[1][1]=-sin（sr）,ai[1][1]=0;
ar[1][2]=-cos（pt）,ai[1][2]=0;
ar[1][3]=sin（st）,ai[1][3]=0;
ar[2][0]=0,ai[2][0]=2*p*vp1*sqrt（p*p*vp1*vp1-1）*vs1/vp1;
ar[2][1]=-（1-2*sin（sr）*sin（sr））,ai[2][1]=0;
ar[2][2]=2*cos（pt）*sin（pt）*k*vs2*vs2/（vs1*vp1）,ai[2][2]=0;
ar[2][3]=（1-2*sin（st）*sin（st））*k*vs2/vs1,ai[2][3]=0;
ar[3][0]=（1-2*sin（sr）*sin（sr））*vp1/vs1,ai[3][0]=0;
ar[3][1]=2*cos（sr）*sin（sr）,ai[3][1]=0;
ar[3][2]=-（1-2*sin（st）*sin（st））*k*vp2/vp1,ai[3][2]=0;
ar[3][3]=2*cos（st）*sin（st）*k*vs2/vp1,ai[3][3]=0;
//*输入系数矩阵
br[0]=cos（sr）,bi[0]=0;
br[1]=sin（sr）,bi[1]=0;
br[2]=1-2*sin（sr）*sin（sr）,bi[2]=0;
br[3]=2*cos（sr）*sin（sr）,bi[3]=0;
//*输入常数矩阵
if（cinv（ar,ai,4）!
=0）//*求系数矩阵的逆矩阵*//
{
tcmul（ar,ai,br,bi,4,4,1,cr,ci）;
rsp=sqrt（cr[0][0]*cr[0][0]+ci[0][0]*ci[0][0]）;
rss=sqrt（cr[1][0]*cr[1][0]+ci[1][0]*ci[1][0]）;
tsp=sqrt（cr[2][0]*cr[2][0]+ci[2][0]*ci[2][0]）;
tss=sqrt（cr[3][0]*cr[3][0]+ci[3][0]*ci[3][0]）;//*求S波入射反射透射系数*//
fprintf（fp2,"%d,%f,%f,%f,%f\n",n[i],rsp,rss,tsp,tss）;//*输出S波入射反射透射系数*//
}
}
elseif（in>=ipp3&&in{
sr=asin（x1）;
st=asin（y1）;
ar[0][0]=-p*vp1,ai[0][0]=0;
ar[0][1]=cos（sr）,ai[0][1]=0;
ar[0][2]=p*vp2,ai[0][2]=0;
ar[0][3]=cos（st）,ai[0][3]=0;
ar[1][0]=0,ai[1][0]=-sqrt（p*p*vp1*vp1-1）;
ar[1][1]=-sin（sr）,ai[1][1]=0;
ar[1][2]=0,ai[1][2]=-sqrt（p*p*vp2*vp2-1）;
ar[1][3]=sin（st）,ai[1][3]=0;
ar[2][0]=0,ai[2][0]=2*p*vp1*sqrt（p*p*vp1*vp1-1）*vs1/vp1;
ar[2][1]=-（1-2*sin（sr）*sin（sr））,ai[2][1]=0;
ar[2][2]=0,ai[2][2]=2*p*vp2*sqrt（p*p*vp2*vp2-1）*k*vs2*vs2/（vs1*vp1）;
ar[2][3]=（1-2*sin（st）*sin（st））*k*vs2/vs1,ai[2][3]=0;
ar[3][0]=（1-2*sin（sr）*sin（sr））*vp1/vs1,ai[3][0]=0;
ar[3][1]=2*cos（sr）*sin（sr）,ai[3][1]=0;
ar[3][2]=-（1-2*sin（st）*sin（st））*k*vp2/vp1,ai[3][2]=0;
ar[3][3]=2*cos（st）*sin（st）*k*vs2/vp1,ai[3][3]=0;
//*输入系数矩阵
br[0]=cos（sr）,bi[0]=0;
br[1]=sin（sr）,bi[1]=0;
br[2]=1-2*sin（sr）*sin（sr）,bi[2]=0;
br[3]=2*cos（sr）*sin（sr）,bi[3]=0;
//*输入常数矩阵
if（cinv（ar,ai,4）!
=0）//*求系数矩阵的逆矩阵
{
tcmul（ar,ai,br,bi,4,4,1,cr,ci）;
rsp=sqrt（cr[0][0]*cr[0][0]+ci[0][0]*ci[0][0]）;
rss=sqrt（cr[1][0]*cr[1][0]+ci[1][0]*ci[1][0]）;
tsp=sqrt（cr[2][0]*cr[2][0]+ci[2][0]*ci[2][0]）;
tss=sqrt（cr[3][0]*cr[3][0]+ci[3][0]*ci[3][0]）;//*求S波入射反射透射系数*//
fprintf（fp2,"%d,%f,%f,%f,%f\n",n[i],rsp,rss,tsp,tss）;//*输出S波入射反射透射系数*//
}
}
elseif（in>=ipp2）
{
sr=asin（x1）;
ar[0][0]=-p*vp1,ai[0][0]=0;
ar[0][1]=cos（sr）,ai[0][1]=0;
ar[0][2]=p*vp2,ai[0][2]=0;
ar[0][3]=0,ai[0][3]=sqrt（p*p*vs2*vs2-1）;
ar[1][0]=0,ai[1][0]=-sqrt（p*p*vp1*vp1-1）;
ar[1][1]=-sin（sr）,ai[1][1]=0;
ar[1][2]=0,ai[1][2]=-sqrt（p*p*vp2*vp2-1）;
ar[1][3]=p*vp2,ai[1][3]=0;
ar[2][0]=0,ai[2][0]=2*p*vp1*sqrt（p*p*vp1*vp1-1）*vs1/vp1;
ar[2][1]=-（1-2*sin（sr）*sin（sr））,ai[2][1]=0;
ar[2][2]=0,ai[2][2]=2*p*vp2*sqrt（p*p*vp2*vp2-1）*k*vs2*vs2/（vs1*vp1）;
ar[2][3]=（1-2*p*p*vs2*vs2）*k*vs2/vs1,ai[2][3]=0;
ar[3][0]=（1-2*sin（sr）*sin（sr））*vp1/vs1,ai[3][0]=0;
ar[3][1]=2*cos（sr）*sin（sr）,ai[3][1]=0;
ar[3][2]=-（1-2*p*p*vs2*vs2）*k*vp2/vp1,ai[3][2]=0;
ar[3][3]=0,ai[3][3]=2*p*vs2*sqrt（p*p*vs2*vs2-1）*k*vs2/vp1;
//*输入系数矩阵
br[0]=cos（sr）,bi[0]=0;
br[1]=sin（sr）,bi[1]=0;
br[2]=1-2*sin（sr）*sin（sr）,bi[2]=0;
br[3]=2*cos（sr）*sin（sr）,bi[3]=0;
//*输入常数矩阵
if（cinv（ar,ai,4）!
=0）//*求系数矩阵的逆矩阵
{
tcmul（ar,ai,br,bi,4,4,1,cr,ci）;
rsp=sqrt（cr[0][0]*cr[0][0]+ci[0][0]*ci[0][0]）;
rss=sqrt（cr[1][0]*cr[1][0]+ci[1][0]*ci[1][0]）;
tsp=sqrt（cr[2][0]*cr[2][0]+ci[2][0]*ci[2][0]）;
tss=sqrt（cr[3][0]*cr[3][0]+ci[3][0]*ci[3][0]）;//*求S波入射反射透射系数*//
fprintf（fp2,"%d,%f,%f,%f,%f\n",n[i],rsp,rss,tsp,tss）;//*输出S波入射反射透射系数*//
}
}
}
fclose（fp1）;
fclose（fp2）;
}
SH波入射时，程序：
#include
#include
#definePI3.1415926
voidmain（）
{
FILE*fp1,*fp2;
inti,n[91];
doublein,ipp,x,y,st,rsr,rsi,tsr,tsi,rsh,tsh,p,q,r[91],t[91],vs1=1585,vs2=1989,den1=2.4,den2=2.24,k=den2/den1;
staticdoublear[2][2],ai[2][2],br[2],bi[2];
fp1=fopen（"snell.csv","w"）;
fp2=fopen（"SH波入射反射透射系数.csv","w"）;
ipp=asin（vs1/vs2）;//*临界角
for（i=0;i<=90;i++）
{
n[i]=i;
in=i*PI/180;
x=sin（in）;
y=sin（in）*vs2/vs1;//*snell定律
r[i]=x;
t[i]=y;
fprintf（fp1,"%d,%f,%f\n",n[i],r[i],t[i]）;
//*输出SH波入射反射透射角度正弦值
if（in{
st=asin（y）;
p=vs1*cos（in）;
q=vs2*cos（st）;
rsh=（p-k*q）/（p+q）;
tsh=（2*p+（1-k）*q）/（p+q）;
fprintf（fp2,"%d,%f,%f\n",n[i],rsh,tsh）;//*输出SH波入射反射透射系数
}
else
{
p=vs1*cos（in）;
q=vs2*sqrt（y*y-1）;
rsr=（p*p-k*q*q）/（p*p+q*q）;
rsi=-（1+k）*p*q/（p*p+q*q）;
tsr=（2*p*p+（1-k）*q*q）/（p*p+q*q）;
tsi=-（1+k）*p*q/（p*p+q*q）;
rsh=sqrt（rsr*rsr+rsi*rsi）;
tsh=sqrt（tsr*tsr+tsi*tsi）;
fprintf（fp2,"%d,%f,%f\n",n[i],rsh,tsh）;//*输出SH波入射反射透射系数
}
}
fclose（fp1）;
fclose（fp2）;
}
升级会员 