24小时热门版块排行榜

返回列表

【奖励】本帖被评价5次，作者baobiao007增加金币 2.6 个

baobiao007

木虫 (职业作家)

应助: 201 (大学生)
金币: 6482.7
帖子: 3050
在线: 1009.9小时
虫号: 505962

[资源] 【分享】相移法偏移程序【原创】

#include
#include
#include
#include"FFT.h"

//将文件filename中的内容读入二维数组a
void MyReadFile(char *filename,float **a,int m,int n)
{
   int i,j;
   FILE *fp;
   fp=fopen(filename,"rb"

;
   for(i=0; i    for(j=0; j    fread(&a[j],sizeof(float),1,fp);
   fclose(fp);
}
//将二维数组a中的数据写入文件filename中
void MyWriteFile(char *filename,float **a,int m,int n)
{
   int i,j;
   FILE *fp;
   fp=fopen(filename,"wb";
   for(i=0; i    for(j=0; j    fwrite(&a[j],sizeof(float),1,fp);
   fclose(fp);
}
//申请二维动态数组的函数
float **MySpace(int m, int n)
{
int i;
float **p;
p=(float **)calloc(m,sizeof(float *));
if(p==NULL)
{  printf("申请空间失败\n";
   exit(1);
}
for(i=0; i p=(float *)calloc(n,sizeof(float));
return p;
}
//释放申请的二维动态数组
void FreeMySpace(float ***p, int m)
{
int i;
for(i=0; i free((*p));
free(*p);
}
//矩阵转置a[m][n]->b[n][m]
void Zhuan(float ***a,float ***b,int m,int n)
{
int i,j;
for(i=0; i for(j=0; j *(*(*b+i)+j)=*(*(*a+j)+i);
}

void main()
{
/////////////////////////////////设置初始参数////////////////////////////////////
const int DNUM=256;//零偏移距剖面接收道数
const int TNUM=512;//时间采样点数
const int ZNUM=256;//深度采样点数
const float DT=0.002;//时间采样间隔2ms
const float DX=12.5;//道间距
const float DZ=5.0;//深度延拓步长
const float pi=3.1415926;
//float V=2000.0;
float V[ZNUM];
char fname1[]="Datapost-8-8-256-512.dat";//原始地震记录文件名
char fnamev[]="Modelpost-8-8-256-256.dat";//速度文件
char fname2[]="after.dat";//偏移后的剖面文件名
////其他参数计算区//////////////////////////////////
float **recordr, **recordi, **poum;
float **midr, **midi, *tempr, *tempi, *temp2r, *temp2i;
float p, **v, vt;
float df;//频率采样间隔
float dkx;//x方向波数采样间隔
float w;//圆频率
float z;//深度
float kz;//z方向圆波数
float kx;
int iw,iz,ix,k1,k2;
df=1.0/(DT*TNUM); dkx=1.0/(DX*DNUM);
///相移偏移过程/////////////////////////////////////
//申请空间
tempr=(float *)calloc(DNUM,sizeof(float));
tempi=(float *)calloc(DNUM,sizeof(float));
temp2r=(float *)calloc(DNUM,sizeof(float));
temp2i=(float *)calloc(DNUM,sizeof(float));
recordr=MySpace(DNUM,TNUM);
recordi=MySpace(DNUM,TNUM);
poum=MySpace(DNUM,ZNUM);
midr=MySpace(TNUM,DNUM);
midi=MySpace(TNUM,DNUM);
v=MySpace(DNUM,ZNUM);
//读入叠加剖面
MyReadFile(fname1,recordr,DNUM,TNUM);
MyReadFile(fnamev,v,DNUM,ZNUM);
for(iz=0; iz       V[iz]=v[0][iz];
for(ix=0; ix if(v[ix][iz] vt=v[ix][iz];
v[ix][iz]=V[iz];
V[iz]=vt;
}
}
for(iz=0; iz printf("%f\n",V[iz]);
//recordr[127][127]=1.0;
//recordr[50][50]=1.0;
//将输入剖面对时间做一维傅里叶变换
k1=log(DNUM)/log(2.0);
k2=log(TNUM)/log(2.0);
if(DNUM>pow(2,k1)) k1=k1+1;
if(TNUM>pow(2,k2)) k2=k2+1;
for(ix=0; ix       fft(recordr[ix],recordi[ix],k2,1);
Zhuan(&recordr,&midr,DNUM,TNUM);
Zhuan(&recordi,&midi,DNUM,TNUM);
//圆频率循环与深度延拓
for(iw=0; iw<=TNUM/2; iw++)//频率循环
{
      w=2.0*pi*iw*df;
if(w==0) continue;
for(ix=0; ix { tempr[ix]=midr[iw][ix];
tempi[ix]=midi[iw][ix];
}
for(iz=0; iz {
         fft(tempr,tempi,k1,1);//对x进行一维傅里叶变换
for(ix=0; ix { //计算相移因子并相乘
kx=2.0*pi*ix*dkx;
if(ix>DNUM/2) kx= 2.0*pi*(DNUM-ix-1)*dkx;
p=1.0-pow(0.5*V[iz]*kx/w,2);
if(p>=0.0){
         kz=w*sqrt(p)/(0.5*V[iz]);
      temp2r[ix]=tempr[ix]*cos(kz*DZ)-tempi[ix]*sin(kz*DZ);
      temp2i[ix]=tempi[ix]*cos(kz*DZ)+tempr[ix]*sin(kz*DZ);
}
else{
temp2r[ix]=0.0;
temp2i[ix]=0.0;
}
}
//kx做一维傅里叶反变换
fft(temp2r,temp2i,k1,-1);
for(ix=0; ix { tempr[ix]=temp2r[ix];
tempi[ix]=temp2i[ix];
            //成像
poum[ix][iz]+=2.0*temp2r[ix];
}
}
}
//输出偏移结果文件
MyWriteFile(fname2,poum,DNUM,ZNUM);
//释放空间
FreeMySpace(&recordr,DNUM);
FreeMySpace(&recordi,DNUM);
FreeMySpace(&midr,TNUM);
FreeMySpace(&midi,TNUM);
FreeMySpace(&poum,DNUM);
free(tempr);free(tempi);
free(temp2r);free(temp2i);
FreeMySpace(&v,DNUM);
}