叠加地震记录的相移波动方程正演模拟

本科生实验报告实验课程数值模型模拟学院名称地球物理学院专业名称勘测技术与工程（石油物探）学生姓名学生学号指导教师熊高君实验地点5417实验成绩2015年5月成都理工大学《地震数值模拟》实验报告实验报告一、实验题目叠加地震记录的相移波动方程正演模拟二、实验目的掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论、实现方法与程序编制，由正演记录初步分析地震信号的分辨率。

三、原理公式1、地震波传播的波动方程设（x,z）为空间坐标，t为时间，地震波传播速度为v(x,z)，则二维介质中任意位置、任意时刻的地震波场为p(z,x,t)：压缩波——纵波。

则二维各向同性均匀介质中地震波传播遵循的声波方程：ð2p(x,z,t)ðx2+ð2p(x,z,t)ðz2=1v2(x,z)ð2p(x,z,t)ðt2(1)2、傅里叶变换的微分性质p(t)与其傅里叶变换的P(ω)的关系：{P(ω)=∫p(t)e−iωt dt∞−∞正傅里叶变换p(t)=12π∫P(ω)e iωt dt∞−∞逆傅里叶变换(2)则有时间微分性质：{ (iω)P (ω)=∫dp (t )dt e −iωt dt ∞−∞ 一阶微分 (iω)2P (ω)=∫d 2p (t )dt 2e −iωt dt ∞−∞二阶微分 (3) ω为频率，ω=2πT ⁄，T 为周期。

同理有空间微分性质：{ (ik)P (k )=∫dp (x )dx e −ikx dx ∞−∞ 一阶微分 (ik)2P (k )=∫d 2p (x )dx 2e −ikx dx ∞−∞二阶微分 (4) k 为波数，k=2πλ⁄， λ为波长 3、地震波传播的相移外推公式令速度v 不随x 变化，只随z 变化，则利用傅里叶变换微分性质（3）和（4）式，把波动方程(1)式变换到频率‐波数域，得：(ik )2P (k,z i ,ω)+ð2P(k,z,ω)ðz 2=(iω)2v (z )2P(k,z,ω) 或：ð2P(k,z,ω)ðz 2=−(ω2v (z )2−k 2)P (k,z,ω) (5) 令：k z 2=ω2v (z )2−k 2 则(5)式的解为： P (k,z,ω)=c 1e −ik z z +c 2e ik z z (6)包括上行波和下行波两项，正演模拟取上行波：P (k,z,ω)=c 1e −ik z z (7)若Z j和Z j+1间隔为∆z，速度v(z)为在此间隔内不随Z变的常数，(7)式实现波场从Z j+1到Z j的延拓，即：P(k,z j,ω)=c1e−ik z∆z在深度Z j+1开始向上延拓到Z j，若延拓深度为零，即：∆Z=Z j+1−Z j，则P(k,z j=z j+1,ω)=c1e−ik z(z j+1−z j)=ce−ik z×0=c (8)对于任意深度Z j+1到Z j的延拓，可得正演模拟中地震波的传播方程（延拓公式P(k,z j,ω)=P(k,z j+1,ω)c1e−ik z(z j+1−z j) (9)4、初始条件和边界条件按照爆炸界面理论，反射界面震源在t=0 时刻同时起爆，此时刻的波场就是震源。

根据不同情况，可直接使用反射系数脉冲或子波作震源。

如果直接使用反射系数作震源脉冲，则初始条件可表示为：r(x,z) t=0p0(x,z,t)={(10)0 t=其他p0(x,z,t)对时间t和空间x做二维傅立叶变换，则得频率-波数域的初始波场p0(k,z,ω)。

边界条件：p0(x,z,t)=r(x,z) t=0,x min<x<x max,z min<z<z max{(11)0 t=其他,x=其他,z=其他其他参数都是在x min<x<x max,z min<z<z max范围内定义的。

5、边界处理（1）边界反射问题把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。

物理上假设探区距X min与X max两个端点很远，在两个端点上收到的反射波很弱。

但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。

在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。

（2）边界强反射的处理镶边法、削波法、吸收边界都能有效消除边界强反射。

削波法就是在波场延拓过程中，每延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。

假设横向总长度为NX ，以两边Lx 道吸波为例，有以下吸波公式：{Abs (Nx −Ix )=Abs (Ix )=√sin (π2×Lx Lx−1) 0<=Ix <Lx Abs (Ix )=1. Ix =其它 (12) 6、数字化根据数字信号处理的采样定理，把连续的信号变为计算机能处理的数字信号，使相移法正演模拟得以实现。

频域抽样定理：一个频谱受限信号 f(t)，如果时间只占据 −t m −t m 的范围，若在频域以不大于 1/2t m 频率间隔 ∆f ≤1/2t m 对信号f(t)的频谱F(ω)采样，则抽样到的离散信号 F 1(ω) 可以唯一表示原信号。

时域抽样定理：一个时间受限信号 f(t)，如果频谱只占据−ωm −ωm 的范围，则信号 f(t)可以用等间隔的抽样值唯一表示出来，而时间 ∆t 抽样间隔必须不大于1/2f m ，ωm =2πf m ，∆t =1/2f m 。

四、实验内容削波法相位移正演模拟（1）点绕射构造和水平层状速度模型（参数如图1所示）的正演数值模拟；1）削波的正演；2）无削波的正演；（2）计算中点和两个边界的信号位置，分析实验结果的正确性；（3）做同样模型的褶积模型数值模拟，对比分析两者的异同；（4）改变绕射点位置、速度，图1 点绕射构造与水平速度模型再做正演模拟；五、方法路线1、参数初始化；2、形成边界削波数据；3、波场初始化；3、Zmax层波场延拓到深度Zmax-1；5、Zi+1层波场延拓到深度Zi；6、重复5，从Iz=Nz-1开始，直到Iz=1，得测线上的频率—空间域波场；7、频率-空间域波场对频率做反傅里叶变换，得时间-空间波场；8、使用Fimage软件显示所得结果。

六、实验结果1、结果显示及分析（点绕射结构）图2 削波（左）正演模拟结果与未削波正演模拟结果（右）由图2分析可知：削波前，由于在两个端点上收到的反射波很弱，造成零边界条件反而成为绝对阻止波通过的强反射面，在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。

削波后，边界假反射的影响消除了，且曲线变得更光滑了。

a、速度V1=5000、V2=5500，深度h=2000时，改变绕射点x的位置；图3 x=Nx/4-1时，反射系数（左）与削波正演模拟结果（右）图4 x=Nx/2-1时，反射系数（左）与削波正演模拟结果（右）由图3与图4对比分析可知：当速度，深度一定时，改变绕射点的位置，曲线随绕射点位置的改变而左右移动，且移动趋势与绕射点位置的移动一致。

b 、绕射点x=Nx/2-1，深度h=2000 ，速度V 2=5500时，改变速度V 1；图5 V 1=4000时，速度模型（左）与削波正演模拟结果（右）V 1=4000V 2=5500图6 V1=6000时，速度模型（左）与削波正演模拟结果（右）由图5与图6对比分析可知：当绕射点位置固定，深度一定时，随着介质速度增大，曲线变得越平缓。

这是由于速度增大时，波的传播时间减小所致。

c、速度V1=5000、V2=5500，绕射点x=Nx/2-1，改变深度h；图7 h=1000时，反射系数（左）与削波正演模拟结果（右）V2=5500 V1=6000图8 h=3000时，反射系数（左）与削波正演模拟结果（右）由图7与图8对比分析可知：当介质速度一定，绕射点位置不变时，随着绕射点深度的增加，曲线变得越来越平滑。

显然是由于深度变大，波的传播时间变小所致。

2、结果分析（水平层速度模型）图9 未削波（左）正演模拟结果与削波正演模拟结果（右）由图9分析可知：削波前，由于在两个端点上收到的反射波很弱，造成零边界条件反而成为绝对阻止波通过的强反射面，在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。

削波后，边界假反射的影响消除了，且曲线变得更光滑了。

a 、深度h=2000 ，速度V 2=5500时，改变速度V 1；图10 V 1=4000时，速度模型（左）与削波正演模拟结果（右）图11 V 1=6000时，速度模型（左）与削波正演模拟结果（右）由图10与图11对比分析可知：水平层深度不变时，随着速度的增加，结果显示的水平层向上移动。

这是由于速度增大时，波的传播时间减小的。

V 1=4000V 2=5500V 2=5500V 1=6000b、速度V1=5000、V2=5500，改变深度h；图12 h=1000时，反射系数（左）与削波正演模拟结果（右）图13 h=3000时，反射系数（左）与削波正演模拟结果（右）由图12与图13对比分析可知：当波的传播速度不改变时，随着水平层深度的增加，结果显示的水平层向下移动。

这是由于深度增加，波的传播时间增大。

七、讨论建议1、实验收获通过本次实验，对叠后相移波动方程正演模拟的过程有了更深入的理解，以及对其公式和意义都有了进一步的认识。

编程能力又有了一定提高。

2、存在问题图14 x=3*Nx/4-1时，反射系数（左）与削波正演模拟结果（右）如图14，当速度，深度一定时，绕射点位置x=3*Nx/4-1，其反射系数不能用Fimage正确显示，由显示的削波正演模拟结果可知，使用削波函数并未起到削波的作用。

3、心得体会本次实验，并不容易。

在实验中，必须要理解理论知识，才能够正确的编写程序。

在实验过程中，遇到许多问题，通过和同学讨论，向老师求解，最终才得出实验结果。

在以后的学习工作中，当一个人的力量不够时，要积极和其他人求教，共同学习，共同进步。

附程序代码：#include <stdlib.h>#include <conio.h>#include <math.h>#include <stdio.h>#include <string.h>#define Nx 128 // Trace Number#define Nt 256 // Record Number#define Nz 100 // Depth Number#define Labs 10 // Length Of Boundary Absorbing#define Dx 20 // Trace Interval#define Dt 0.004 // Record Interval#define Dz 20 // Depth Interval#define Pai 3.1415926//function01: Judge the power of 2int Po2Judge(int N){int k=0;long Ln=0;for(k=0;N-Ln>0;k++)Ln=(long)pow(2,k);Ln=(long)pow(2,k-1);if(fabs(Ln-N)>=1)return(0);return(1);}//function02: Absorb Bounderyint Absorb( ){FILE *fp_Abs;int Ix;float Abs[Nx];if((fp_Abs=fopen("Absb.dat","wb"))==NULL)printf("Connot open file ""Absb""");for(Ix=0;Ix<Nx;Ix++)Abs[Ix]=1.0;//1.Absorb Boundary Initializing ？for(Ix=0;Ix<Labs;Ix++){Abs[Ix]=sqrt(sin((Pai/2)*Ix/(Labs-1)));Abs[Nx-Ix-1]=Abs[Ix];//2.Absorb Boundary Compute？}for(Ix=0;Ix<Nx;Ix++)fwrite(&Abs[Ix],sizeof(Abs[Ix]),Nx,fp_Abs);//3.Byte Number of fclose(fp_Abs);return(1);}//function03: Form Reflect Structure Modelint Rflct( ){FILE *fp_Rflct;int Ix,Iz;float Rflct[Nz];if((fp_Rflct=fopen("Rflct.dat","wb"))==NULL)printf("Connot open file ""Reflection""");for(Ix=0;Ix<Nx;Ix++){for(Iz=0;Iz<Nz;Iz++){Rflct[Iz]=0.;if(Ix==Nx/4-1&&Iz==Nz/2-1) Rflct[Iz]=1;fwrite(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct);//4.}}fclose(fp_Rflct);return(1);}//function04: Form Velocity Modelint Vlcty( ){FILE *fp_Vlcty;int Ix,Iz;float Vlcty[Nz];if((fp_Vlcty=fopen("Vlcty.dat","wb"))==NULL)printf("Connot open file ""Vlcty""");for(Ix=0;Ix<Nx;Ix++){for(Iz=0;Iz<(int)(3*Nz/4);Iz++){Vlcty[Iz]=5000.;fwrite(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);//5.}for(Iz=(int)(3*Nz/4);Iz<Nz;Iz++){Vlcty[Iz]=5500.;fwrite(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);//6.}}fclose(fp_Vlcty);return(1);}//Fourier Transform: FFT(i=0) or IFFT(l=1)int kkfft(float pr[], float pi[], int n, int l){int it,m,is,i,j,nv,l0,il=0;float p,q,s,vr,vi,poddr,poddi;float fr[4096],fi[4096];int k=0;long Ln=0;for(k=0;n-Ln>0;k++){Ln=(long)pow(2,k);}k=k-1;for (it=0; it<=n-1; it++){m = it;is = 0;for(i=0; i<=k-1; i++){j = m/2;is = 2*is+(m-2*j);m = j;}fr[it] = pr[is];fi[it] = pi[is];}pr[0] = 1.0;pi[0] = 0.0;p = 6.283185306/(1.0*n);pr[1] = (float) cos(p);pi[1] = -(float)sin(p);if (l!=0)pi[1]=-pi[1];for (i=2; i<=n-1; i++){p = pr[i-1]*pr[1];q = pi[i-1]*pi[1];s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1]);pr[i] = p-q;pi[i] = s-p-q;}for (it=0; it<=n-2; it=it+2){vr = fr[it];vi = fi[it];fr[it] = vr+fr[it+1];fi[it] = vi+fi[it+1];fr[it+1] = vr-fr[it+1];fi[it+1] = vi-fi[it+1];}m = n/2;nv = 2;for (l0=k-2; l0>=0; l0--){m = m/2;nv = 2*nv;for(it=0; it<=(m-1)*nv; it=it+nv)for (j=0; j<=(nv/2)-1; j++){p = pr[m*j]*fr[it+j+nv/2];q = pi[m*j]*fi[it+j+nv/2];s = pr[m*j]+pi[m*j];s = s*(fr[it+j+nv/2]+fi[it+j+nv/2]);poddr = p-q;poddi = s-p-q;fr[it+j+nv/2] = fr[it+j]-poddr;fi[it+j+nv/2] = fi[it+j]-poddi;fr[it+j] = fr[it+j]+poddr;fi[it+j] = fi[it+j]+poddi;}}if(l!=0){for(i=0; i<=n-1; i++){fr[i] = fr[i]/(1.0*n);fi[i] = fi[i]/(1.0*n);}}if(il!=0){for(i=0; i<=n-1; i++){pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);if(fabs(fr[i])<0.000001*fabs(fi[i])){if ((fi[i]*fr[i])>0)pi[i] = 90.0;elsepi[i] = -90.0;}elsepi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;}}for(i=0;i<n;i++){pr[i]=fr[i];pi[i]=fi[i];}return ( 1 );}//function05: Form Inition Wave Field: Reflect Coefficiency //int WvFld0(){FILE *fp_Rflct,*fp_Wfld0r,*fp_Wfld0i;int Ix,Iz,It;float Rflct[Nz];float Wfld0r[Nt],Wfld0i[Nt];if((fp_Wfld0r = fopen("Wfld0r.dat","wb"))==NULL)printf("Connot open Wfld0r.dat");if((fp_Wfld0i = fopen("Wfld0i.dat","wb"))==NULL)printf("Connot open Wfld0i.dat");if((fp_Rflct = fopen("Rflct.dat", "rb"))==NULL)printf("Connot open Rflct.dat");for(Ix=0;Ix<Nx;Ix++){for(Iz=0;Iz<Nz;Iz++){fread(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct);for(It=0;It<Nt;It++){Wfld0r[It]=0.0;//7Wfld0i[It]=0.0;//8.Initial Wave field Initializing: }Wfld0r[0]=Rflct[Iz];//9.Initial Wave field Initializingif( kkfft(Wfld0r,Wfld0i,Nt,0 ) !=1 )//10.FFT or Inver FFT?{printf("FFT is error");exit(0);}for(It=0;It<Nt/2+1;It++)//{fwrite(&Wfld0r[It],sizeof(Wfld0r[It]),1,fp_Wfld0r);//11fwrite(&Wfld0i[It],sizeof(Wfld0i[It]),1,fp_Wfld0i);//12}}}fclose(fp_Rflct);fclose(fp_Wfld0r);fclose(fp_Wfld0i);return(1);}// Function06.1.1: Read In Velocity Data and Absorb Boundary Data int ReadVlctyAbsb(float Vlcty[],float Absb[]){FILE *fp_Vlcty,*fp_Absb;int Iz,Ix;if((fp_Vlcty = fopen("Vlcty.dat", "rb")) ==NULL)printf("Connot open Rflct.dat" );for(Iz=0;Iz<Nz;Iz++){fread(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);//19.}fclose(fp_Vlcty);if((fp_Absb = fopen("Absb.dat", "rb")) ==NULL)printf("Connot open Absb.dat" );for(Ix=0;Ix<Nx;Ix++){fread(&Absb[Ix],sizeof(Absb[Ix]),1,fp_Absb);//20.}fclose(fp_Absb);remove("Absb.dat");return(1);}// 06.1.2.1: Read Data From (Ix,Iz,Iw)Oder to (Iw,Iz,Ix)Oderint ReadIxIzIwToIwIzIx(FILE *fp_Wfld0r,FILE *fp_Wfld0i,floatWfld0r[],float Wfld0i[],int Iz,int Iw){int Ix,Nw=Nt;long Mgrtn;for(Ix=0;Ix<Nx;Ix++){Mgrtn=(Ix*Nz+1+Iz)*(Nt/2+1)+Iw;//25.fseek(fp_Wfld0r, sizeof(Wfld0r[Ix])*Mgrtn,0);//26fread(&Wfld0r[Ix],sizeof(Wfld0r[Ix]),1,fp_Wfld0r); //27fseek(fp_Wfld0i, sizeof(Wfld0i[Ix])*Mgrtn,0);//28fread(&Wfld0i[Ix],sizeof(Wfld0i[Ix]),1,fp_Wfld0i); //29. }return(1);}// 06.1.2: Form New Wave Fieldint FrmNewWfld(FILE *fp_Wfld0r,FILE *fp_Wfld0i,float Wfldr[],float Wfldi[],float Absb[],int Iz,int Iw){int ReadIxIzIwToIwIzIx(FILE *fp_Wfld0r,FILE *fp_Wfld0i,float Wfld0r[],float Wfld0i[],int Iz,int Iw);int Ix,Nw=Nt;float Wfld0r[Nx],Wfld0i[Nx];if(ReadIxIzIwToIwIzIx(fp_Wfld0r,fp_Wfld0i,Wfld0r,Wfld0i,Iz,Iw)!= 1){printf("exp_ikzDz is error");exit(0);}for(Ix=0;Ix<Nx;Ix++){Wfldr[Ix]+=Wfld0r[Ix];//21.Wfldi[Ix]+=Wfld0i[Ix];//22.Wfldr[Ix]=Wfldr[Ix]*Absb[Ix];//23.Wfldi[Ix]=Wfldi[Ix]*Absb[Ix];//24.}return(1);}// 06.1.3.1: Compute out PhaseShift Dataint exp_ikzDz(float eikzdz[],int Ix,float Vc,int Iw,float Dw,float Dkx) {float kz;eikzdz[0] = 0.;eikzdz[1] = 0.;kz =sqrt(pow(Iw*Dw/Vc,2)-pow(Ix*Dkx,2));//38.if(kz >0.0)//39.{eikzdz[0] = (float)cos(kz*Dz);//40.eikzdz[1] = (float)sin((-1)*kz*Dz);//41.}return(1);}// 06.1.3: Extrapolate One Depth Stepint MoveOneDz(float Wfldr[],float Wfldi[],float Vz,float Dkx,float Dw,int Iw){int Ikx,Nkx=Nx;float kz[2],Wfld_r,Wfld_i;if(kkfft(Wfldr,Wfldi,Nx,0)!=1 ) //30{printf("FFT is error");exit(0);}for(Ikx=0;Ikx<Nkx/2+1;Ikx++) //31.{if( exp_ikzDz(kz,Ikx,Vz,Iw,Dw,Dkx) !=1 ){printf("exp_ikzDz is error");exit(0);}Wfld_r =kz[0]*Wfldr[Ikx]-kz[1]*Wfldi[Ikx];//32Wfld_i =kz[0]*Wfldi[Ikx]+kz[1]*Wfldr[Ikx];//33Wfldr[Ikx] = Wfld_r;Wfldi[Ikx] = Wfld_i;if((Ikx!=0)&&(Ikx!=Nkx/2))//34{Wfld_r =kz[0]*Wfldr[Nkx-Ikx]-kz[1]*Wfldi[Nkx-Ikx];//35.Wfld_i =kz[0]*Wfldi[Nkx-Ikx]+kz[1]*Wfldr[Nkx-Ikx];//36.Wfldr[Nkx-Ikx] = Wfld_r;Wfldi[Nkx-Ikx] = Wfld_i;}}if(kkfft(Wfldr, Wfldi,Nkx,1) !=1 )//37.{printf("FFT is error");exit(0);}return(1);}//function06.1: PhaseShift calculateint PhaseShift( ){int ReadVlctyAbsb(float *,float *);int FrmNewWfld(FILE *,FILE *,float *,float *,float *,int,int);int MoveOneDz (float *,float *,float,float,float,int);FILE *fp_Wfldr,*fp_Wfldi,*fp_Wfld0r,*fp_Wfld0i;int Ix,Iz,Iw,Nw=Nt;float Vlcty[Nz],Absb[Nx];float Wfldr[Nx],Wfldi[Nx];float Dkx,Dw;Dkx=(float)(Pai/Dx)/Nx;//13.Dw =(float)(Pai/Dt)/Nt;//14.if(ReadVlctyAbsb(Vlcty,Absb)!=1){printf("ReadVlctyAbsb Error");exit(0);}if((fp_Wfld0r = fopen("Wfld0r.dat","rb")) ==NULL)printf("Connot open Wfld0r.dat");if((fp_Wfld0i = fopen("Wfld0i.dat","rb")) ==NULL)printf("Connot open Wfld0i.dat");if((fp_Wfldr = fopen("Wfldr.dat", "wb")) ==NULL)printf("Connot open Wfldr.dat" );if((fp_Wfldi = fopen("Wfldi.dat", "wb")) ==NULL)printf("Connot open Wfldi.dat" );for(Iw=0;Iw<Nw/2+1;Iw++){printf("Iw=%d\n",Iw);for(Ix=0;Ix<Nx;Ix++){Wfldr[Ix]=0.0;//15.Wfldi[Ix]=0.0;//16}for(Iz=Nz-1;Iz>0;Iz--){if( FrmNewWfld(fp_Wfld0r,fp_Wfld0i,Wfldr,Wfldi,Absb,Iz,Iw)!=1) {printf("FrmNewWfld Error");exit(0);}if(MoveOneDz(Wfldr,Wfldi,(float)(Vlcty[Iz]/2.),Dkx,Dw,Iw)!=1) {printf("WfldMoveOneStep Error");exit(0);}}for(Ix=0;Ix<Nx;Ix++){fwrite(&Wfldr[Ix],sizeof(Wfldr[Ix]),1,fp_Wfldr);//17 fwrite(&Wfldi[Ix],sizeof(Wfldi[Ix]),1,fp_Wfldi);//18 }}fclose(fp_Wfld0r); remove("Wfld0r.dat");fclose(fp_Wfld0i); remove("Wfld0i.dat");fclose(fp_Wfldr); fclose(fp_Wfldi);return(1);}// Function03: Wave Field Fourier Transform From Frequency Domain to Time Domainint Frqcy2Time( ){FILE *fp_Wfldr,*fp_Wfldi;FILE *fp_Record;int Ix,It,Iw,Nw=Nt;float Wfldtr[Nt],Wfldti[Nt];long AddFrmStrt;if((fp_Wfldr = fopen("Wfldr.dat", "rb"))==NULL){printf("Connot open Wfldr.dat" );exit(0);}if((fp_Wfldi = fopen("Wfldi.dat", "rb"))==NULL){printf("Connot open Wfldi.dat" );exit(0);}if((fp_Record = fopen("Record.dat", "wb"))==NULL){printf("Connot open Record.dat");exit(0);}for(Ix=0;Ix<Nx;Ix++){for(Iw=0;Iw<Nw/2+1;Iw++)//42.{AddFrmStrt=Iw*Nx+Ix;//43fseek(fp_Wfldr, sizeof(Wfldtr[Iw])*AddFrmStrt,0);//44fread(&Wfldtr[Iw],sizeof(Wfldtr[Iw]), 1,fp_Wfldr); //45fseek(fp_Wfldi, sizeof(Wfldti[Iw])*AddFrmStrt,0);//46fread(&Wfldti[Iw],sizeof(Wfldti[Iw]), 1,fp_Wfldi); //47if((Iw!=0)&&(Iw!=Nw/2)) //48.{Wfldtr[Nw-Iw]=Wfldtr[Iw]; //49.Wfldti[Nw-Iw]=-Wfldti[Iw]; //50.}}if(kkfft(Wfldtr,Wfldti,Nw,1)!=1) //51.{printf("FFT is error");exit(0);}for(It=0;It<Nt;It++){fwrite(&Wfldtr[It],sizeof(Wfldtr[It]),1,fp_Record);//52.}}fclose(fp_Wfldr); remove("Wfldr.dat" );fclose(fp_Wfldi); remove("Wfldi.dat" );fclose(fp_Record);return(1);}//function06: PhaseShift Forward Modelingint PsFrwd( ){int PhaseShift( );int Frqcy2Time( );if( PhaseShift( ) !=1 ){ printf("PhaseShift is error\n"); exit(0); } if( Frqcy2Time( ) !=1 ){ printf("Frqcy2Time is error\n"); exit(0); } return(1);}void main( ){if( Po2Judge(Nt) !=1 ){ printf("Nt=%d is not the Power of 2\n",Nt); exit(0); } if( Po2Judge(Nx) !=1 ){ printf("Nx=%d is not the Power of 2\n",Nx); exit(0); }if( Absorb( ) !=1 ) { printf("Absorb is error\n"); exit(0); } if( Rflct( ) !=1 ) { printf("Rflction is error\n"); exit(0); } if( Vlcty( ) !=1 ) { printf("Vlcty is error\n"); exit(0); }if( WvFld0( ) !=1 ) { printf("WvFld is error\n"); exit(0); } if( PsFrwd( ) !=1 ) { printf("PsFwrd is error\n"); exit(0); } }。