石油软件概论实验报告.docx
《石油软件概论实验报告.docx》由会员分享，可在线阅读，更多相关《石油软件概论实验报告.docx（30页珍藏版）》请在冰豆网上搜索。
石油软件概论实验报告
本科生实验报告
实验课程石油地球物理软件技术概论
学院名称地球物理学院
专业名称勘查技术与工程（石油物探）
学生姓名
学生学号
指导教师赵宪生
实验地点地球物理学院5417
实验成绩
二〇一五年4月-二〇一五年5月
目录
摘要1
第一章seg-y数据滤波2
1.1方法2
1.2源程序2
1.3结果分析8
1.3.1滤波因子8
1.3.2原始数据及结果显示与分析9
第二章自动统计剩余时差校正13
2.1方法13
2.2源程序14
2.3结果分析19
第三章交互集成软件21
3.1方法21
3.2源程序21
3.3结果显示22
参考文献23
心得及评语24
资料处理软件与交互软件设计
摘要
地震勘探采集的数据需要经过一定的处理，才能被解释使用，地震勘探采集数据多为seg_y格式，石油勘探采集处理越来越依赖于计算机技术的发展，使用计算机技术实现地球物理处理软件设计、交互软件设计和图形软件设计及可视化和集成化基本技术。
本实验报告中，主要针对seg_y格式数据进行滤波处理，比较带通滤波和低通滤波的结果差异；对给定信号进行剩余时差校正，并写出相应的C语言源程序；并在MFC框架程序下，设计出相应的交互集成软件。
关键字：
seg_y，带通滤波，低通滤波，剩余时差校正，交互集成软件
第一章seg-y数据滤波
1.1方法
分析已知seg-y数据格式的程序，补充该程序中计算滤波因子的子函数段，选择参数对数据文件366.sgy完成滤波处理，以及对某一道数据做谱分析，比较低通滤波和带通滤波及参数对结果的影响。
用Fimage.exe程序显示原数据和滤波处理数据，比较分析结果。
1.2源程序
#include
#include
#include
#include
#defineLEN3200
#defineNQ230000
#definePI3.1415926
//==========低通滤波因子=========//
floatfilter1（floath[],intn）
{
inti;
floatfc;
doublet,dt=0.002;
printf（"输入fc:
\n"）;
scanf（"%f",&fc）;
for（i=1;i<=n/2;i++）
{
t=i*dt;
h[i]=sin（2*PI*fc*t）/（PI*t）;
h[i+n/2]=h[i];
}
h[25]=2*fc;
for（i=1;i<=n/2;i++）
{
h[n/2-i]=h[n/2+i];
}
}
//============带通滤波因子=========//
floatfilter2（floath[],intn）
{
inti;
floatf1,f2,df,f0;
doublet，dt=0.002;
printf（"输入f1、f2:
\n"）;
scanf（"%f%f",&f1,&f2）;
f0=（f1+f2）/2;
df=abs（f2-f1）/2;
for（i=0;i<=n/2;i++）
{
t=i*dt;
h[i]=2*cos（2*PI*f0*t）*sin（2*PI*df*t）/（PI*t）;
h[i+n/2]=h[i];
}
h[25]=4*df;
for（i=0;i<=n/2;i++）
{
h[n/2-i]=h[n/2+i];
}
}
//==========求褶积函数==========//
intcvbx0（floatx[],floath[],floatr[],intkc,intnt）
{
intm,np,i,j;
floatsum;
m=kc/2;
for（np=0;npr[np]=0.0;
for（i=0;i{
sum=0.0;
for（j=0;j{
if（（i-m+j）<0||（i-m+j）>=nt）continue;
sum=sum+x[i-m+j]*h[j];
}
r[i]=sum;
}
}
//================傅里叶函数==============//
intIBIT（intj,intm）
{
inti,it,j1,j2;
it=0;
j1=j;
for（i=1;i<=m;i++）
{
j2=j1/2;
it=it*2+（j1-2*j2）;
j1=j2;
}
returnit;
}
voidfft（float*xr,float*xi,intnr,floatT）
{
inti,j,k,l,n,n2,nr1,i1,j1,k2,k1n2;
floatc,s,p,tr,ti,trc,tic,ars;
n=（int）（pow（2,nr））;
n2=n/2;
nr1=nr-1;
k=0;
for（l=1;l<=nr;l++）
{
loop1:
for（j=1;j<=n2;j++）
{
k2=k/（int）（pow（2,nr1））;
p=（float）（IBIT（k2,nr））;
ars=（float）（6.2831852*p/（float）（n））;
c=（float）（cos（ars））;
s=（float）（sin（ars））;
k1n2=k+n2;
tr=xr[k1n2]*c+xi[k1n2]*s*T;
ti=xi[k1n2]*c-xr[k1n2]*s*T;
xr[k1n2]=xr[k]-tr;
xi[k1n2]=xi[k]-ti;
xr[k]=xr[k]+tr;
xi[k]=xi[k]+ti;
k++;
}
k+=n2;
if（k{
gotoloop1;
}
else
{
k=0;
nr1=nr1-1;
n2/=2;
}
}
for（j=1;j<=n;j++）
{
i=IBIT（j-1,nr）+1;
if（i>j）
{
j1=j-1;
i1=i-1;
trc=xr[j1];
tic=xi[j1];
xr[j1]=xr[i1];
xi[j1]=xi[i1];
xr[i1]=trc;
xi[i1]=tic;
}
}
if（T<0.0）
{
for（j=0;j<=n;j++）
{
xr[j]/=n;
xi[j]/=n;
}
}
}
voidmain（）
{
floatki[NQ2],ko[NQ2],h[51],pr[1024],pi[1024];
intnflg,i,n,l,p;
intlen,len1,jhd;
charnamei[50],nameo[50];
FILE*fp1,*fp2,*fp3,*fp4,*fp5;
union
{
inthd[NQ2];
floatdata[NQ2];
charbuf[NQ2];
shortintshd[NQ2];
}dat;
//****************************/
printf（"inputfilename:
"）;
scanf（"%s",&namei）;
printf（"outputfilename:
"）;
scanf（"%s",&nameo）;
fp1=fopen（namei,"rb"）;
if（fp1==NULL）{printf（"can'topeninputsegyfile\n"）;exit（0）;}
fp2=fopen（nameo,"wb"）;
if（fp2==NULL）{printf（"can'topenoutputsegyfile\n"）;exit（0）;}
fp3=fopen（"滤波因子.csv","w"）;
if（fp3==NULL）{printf（"can'topenfilterfactorfile\n"）;exit（0）;}
fp4=fopen（"振幅谱.csv","w"）;
if（fp4==NULL）{printf（"can'topenamplitudefile\n"）;exit（0）;}
fp5=fopen（"相位谱.csv","w"）;
if（fp5==NULL）{printf（"can'topenphasefile\n"）;exit（0）;}
fread（&dat.buf[0],1,LEN,fp1）;
fwrite（&dat.buf[0],1,LEN,fp2）;
fread（&dat.buf[0],1,400,fp1）;
len1=dat.shd[10];
printf（"samplenumbers:
%d\n",len1）;
printf（"samplerate（ms）:
%u\n",dat.shd[8]/1000）;
printf（"seg_yflag:
%d\n",dat.shd[12]）;
/**********************************/
fwrite（&dat.buf[0],1,400,fp2）;
nflg=dat.shd[12];
if（nflg==1）len=len1*4+240;
if（nflg==2）len=len1*4+240;
if（nflg==3）len=len1*2+240;
printf（"\ninputrecordlength:
%d\n",len）;
n=0;
jhd=60;
//==========计算滤波因子===========//
printf（"请输入p:
\n"）;
scanf（"%d",&p）;
if（p==1）
{
filter1（h,51）;
for（i=0;i<51;i++）
fprintf（fp3,"%f\n",h[i]）;
}
if（p==2）
{
filter2（h,51）;
for（i=0;i<51;i++）
fprintf（fp3,"%f\n",h[i]）;
}
if（（p!
=1）&&（p!
=2））
printf（"输入错误！
请重新输入p:
\n"）;
//*******************************//
while（（l=fread（&dat.buf[0],1,len,fp1））>0）
{
for（i=0;iki[i]=dat.data[i+jhd];//dat.hd[],dat.shd[]
cvbx0（ki,h,ko,51,len1）;//对数据进行各种处理
for（i=0;idat.data[i+jhd]=ko[i];//dat.hd[],dat.shd[]
fwrite（&dat.buf[0],1,len,fp2）;
n++;
//=========对第100道做谱分析==========//
if（n==100）
{
for（i=0;i<1024;i++）
{
pr[i]=ki[i];
pi[i]=0;
}
fft（pr,pi,10,1.0）;
for（i=0;i<1024;i++）
{
fprintf（fp4,"%f\n",sqrt（pr[i]*pr[i]+pi[i]*pi[i]））;
fprintf（fp5,"%f\n",atan（pi[i]/pr[i]））;
}
}
printf（"tracenumber:
%d\n",n）;
}
fclose（fp1）;
fclose（fp2）;
fclose（fp3）;
fclose（fp4）;
fclose（fp5）;
}
1.3结果分析
1.3.1滤波因子
图1-1低通滤波因子
图1-2带通滤波因子
1.3.2原始数据及结果显示与分析
（a）、原始数据图（b）、第100道振幅谱图
图1-3原始数据及第100道振幅谱图
（a）、fc=45（b）、fc=50
（c）、fc=55d、fc=60
图1-4不同频率的低通滤波结果图
（a）、f1=25,f2=45（b）、f1=20,f2=50
（c）、f1=15,f2=55（d）、f1=10,f2=60
图1-5不同频率范围的带通滤波结果图
由图1-3（a）可知：
原始地震数据文件366.sgy所记录的地震记录能量分布不均匀，连续性差，噪声干扰大，信噪比低；
由图1-3（b）可知：
通过对第100道数据进行谱分析，可以发现其主频在50Hz左右；
由图1-4与图1-3（a）对比可知：
经过低通滤波后，所得的结果连续性有了一定好转，噪声减少，信噪比提高了；
由图1-4（a）（b）（c）（d）对比可知：
低通滤波过程中，截止频率fc对滤波结果影响较大，且越靠近主频，滤波效果越好；
由图1-5与图1-3（a）对比可知：
经过带通滤波后，所得的结果连续性有了一定好转，噪声减少，信噪比提高了；
由图1-5（a）（b）（c）（d）对比可知：
带通滤波过程中，低截止频率f1、高截止频率f2对滤波结果影响较大，若频率范围太小，则会导致有效波丢失（如a），若频率范围太大，则会导致噪声过多，信噪比低，滤波效果差（如d），达不到滤波目的；
由图1-4与图1-5对比可知：
低通滤波和带通滤波都可在一定程度上提高信噪比，使连续性变好；由于带通滤波和低通滤波的通频不同，带通滤波可将低频（低于f1）和高频（高于f2）干扰均可滤去，而低通滤波只能滤去高频干扰（高于fc），故而，带通滤波比低通滤波效果要好一些。
第二章自动统计剩余时差校正
2.1方法
反映地下界面的反射波时距曲线或同相轴一般是双曲线形状的。
其中只有在激发点处接收到的反射波时间（t0）代表界面的法线反射时间，故必须将各个观测点的时间值都变成相应各点的法线反射时间，时距曲线或同相轴才与地下界面的形态一致。
所以，必须从各观测点的时间值中减去一个相应的校正值。
当界面水平时，它等于观测时间减去法线反射时间。
即使对同一反射界面的相同深度，由于各接收点距激发点远近不同，校正量也不同；而对同一道来说，由浅层至深层的校正量亦不同，校正量是变化的，故称动校正。
其方法步骤如下：
1）读取原始文件数据；
2）求出模型道；
3）利用相关函数求出校正时差；
4）对道集进行时差校正；
5）迭代处理，一直到精度满足为止；
6）写文件，用fimage2012显示。
2.2源程序
#include
#include
#definePI3.1415926
//===========低通滤波因子===========//
floatfilter1（floath[],intn）
{
inti;
floatfc;
doublet,dt=0.002;
printf（"输入fc:
\n"）;
scanf（"%f",&fc）;
for（i=1;i<=n/2;i++）
{
t=i*dt;
h[i]=sin（2*PI*fc*t）/（PI*t）;
h[i+n/2]=h[i];
}
h[25]=2*fc;
for（i=1;i<=n/2;i++）
h[n/2-i]=h[n/2+i];
}
//========带通滤波因子========//
floatfilter2（floath[],intn）
{
inti;
floatf1,f2,df,f0;
doublet,dt=0.002;
printf（"输入f1、f2:
\n"）;
scanf（"%f%f",&f1,&f2）;
f0=（f1+f2）/2;
df=abs（f2-f1）/2;
for（i=0;i<=n/2;i++）
{
t=i*dt;
h[i]=2*cos（2*PI*f0*t）*sin（2*PI*df*t）/（PI*t）;
h[i+n/2]=h[i];
}
h[25]=4*df;
for（i=0;i<=n/2;i++）
h[n/2-i]=h[n/2+i];
}
//========求褶积函数======//
intcvbx0（floatx[],floath[],floatr[],intkc,intnt）
{
intm,np,i,j;
floatsum;
m=kc/2;
for（np=0;npfor（i=0;i{
sum=0.0;
for（j=0;j{
if（（i-m+j）<0||（i-m+j）>=nt）continue;
sum=sum+x[i-m+j]*h[j];
}
r[i]=sum;
}
}
//*******求模型道**********//
floatmodel（floatp[20][100],floaty[100],intn,intm）
{
inti,j;
floatsum,t;
for（j=0;j{
sum=0;
for（i=0;i{
t=p[i][j];
sum+=t;
}
y[j]=sum/n;
}
}
//***********求相关************//
floatxcorr（floatt[],floaty[],floatr[],intm）
{
inti,j,k;
floatsum;
for（k=-m+1;k{
sum=0;
for（i=0;i{
j=i+k;
if（（j>=0）&&（j<=m））
sum+=t[i]*y[j];
}
r[k+m-1]=sum;
}
}
//***********求时差***********//
floattime（floatxc[][199],intmax[]）
{
inti,j;
floatt;
for（i=0;i<20;i++）
{
t=xc[i][0];
for（j=1;j<199;j++）
{
if（xc[i][j]>=t）
{
t=xc[i][j];
max[i]=j;
}
}
}
}
//********动校正*************//
floatmove（floatp[][100],intmax[],floatym[][100]）
{
inti,j,k;
for（i=0;i<20;i++）
{
k=max[i]-99;
for（j=0;j<100;j++）
ym[i][j]=p[i][j-k];
}
}
voidmain（）
{
inti,j,k,m=0,J,max[20]={0};
floatx[100],y[100],h[51],p[20][100],ym[20][100];
floatt[100]={0},xc[20][199],r[199]={0},z[199]={0};
FILE*fp,*fp1,*fp2;
//====================================//
fp1=fopen（"nose20-100","rb"）;
if（fp1==NULL）{printf（"can'topeninputsegyfile\n"）;exit（0）;}
fp2=fopen（"nose20-100out","wb"）;
if（fp2==NULL）{printf（"can'topeninputsegyfile\n"）;exit（0）;}
fp=fopen（"滤波","wb"）;
printf（"选择低通
（1）或带通
（2）:
\n"）;
scanf（"%d",&J）;
if（J==1）
filter1（h,51）;
if（J==2）
filter2（h,51）;
for（i=0;i<20;i++）
{
fread（&x[0],4,100,fp1）;
for（j=0;j<100;j++）
printf（"%f\n",x[j]）;
cvbx0（x,h,y,51,100）;
for（j=0;j<100;j++）
p[i][j]=y[j];
fwrite（&y[0],4,100,fp）;
}
//========校正========//
lp:
model（p,y,20,100）;
for（i=0;i<20;i++）
{
for（j=0;j<100;j++）
t[j]=p[i][j];
xcorr（t,y,r,100）;
for（k=0;k<199;k++）
xc[i][k]=r[k];
}
time（xc,max）;
move（p,max,ym）;
m++;
if（m<5）
gotolp;
for（i=0;i<20;i++）
fwrite（&ym[i][0],4,100,fp2）;
fclose（fp）;
fclose（fp1）;
fclose（fp2）;
}
2.3结果分析
（a）、校正前（b）、校正后
图2-1无噪音数据时差校正图
（a）、校正前（b）、校正后
图2-2有噪声数据时差校正前后图
由图2-1（a）与图2-1（b）（或图2-1（a）与图2-1（b））对比可知：
进行时差校正后，同相轴对齐到同一直线上，便于水平叠加；
由图2-1（b）与图2-2（b）对比可知：
含有噪声的地震记录数据，进行相同的时差校正后，有噪声的数据记录不能完全校正到同一直线上。
第三章交互集成软件
3.1方法
1.建立一个单文档框架程序；
2.设计一个菜单；
3.设计一对话框，输入参数；
4.选菜单绘制图形;
5.对话框输入参数控制显示.
3.2源程序
voidCMy1111View:
:
OnMenuitem32771（）
{
//TODO:
Addyourcommandhandlercodehere
CClientDC*pdc=newCClientDC（this）;
CPenpen;
pen.CreatePen（PS_SOLID,5,RGB（250,0,0））;
CPen*oldpen=（CPen*）pdc->SelectObject（&pen）;
constdoublePI=3.14159;
inti,cx=10,x2=50,y2,m=50;
doublea[51];
intn=51;
floatfc=45;
doublet,dt=0.002;
for（i=0;i<=n/2;i++）
{
t=i*dt;
a[i]=sin（2*PI*fc*t）/（PI*t）;
a[i+n/2]=a[i];
}
a[25]=2*fc;
for（i=0;i<=n/2;
升级会员 