| 查看: 3142 | 回复: 6 | ||
| 【奖励】 本帖被评价5次，作者baobiao007增加金币 2.6 个 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[资源]
【分享】相移法偏移程序【原创】
|
||
|
#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 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 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 return p; } //释放申请的二维动态数组 void FreeMySpace(float ***p, int m) { int i; for(i=0; i 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 } 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 for(ix=0; ix v[ix][iz]=V[iz]; V[iz]=vt; } } for(iz=0; 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 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 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 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); } |
回复此楼
5楼2016-01-28 03:17:14
3楼2015-03-29 19:33:50
4楼2015-03-29 21:58:30
6楼2016-01-28 08:16:25
投票:
30