疊加地震記錄的相移波動(dòng)方程正演模擬

1、本科生實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)課程 數(shù)值模型模擬 學(xué)院名稱(chēng) 地球物理學(xué)院 專(zhuān)業(yè)名稱(chēng) 勘測(cè)技術(shù)與工程（石油物探） 學(xué)生姓名 學(xué)生學(xué)號(hào) 指導(dǎo)教師 熊高君 實(shí)驗(yàn)地點(diǎn) 5417 實(shí)驗(yàn)成績(jī) 2015年5月 成都理工大學(xué)地震數(shù)值模擬實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)時(shí)間2015年 5月開(kāi)課單位地球物理學(xué)院指導(dǎo)教師熊高君實(shí)驗(yàn)題目： 疊加地震記錄的相移波動(dòng)方程正演模擬姓名學(xué)號(hào)班級(jí)專(zhuān)業(yè)勘測(cè)技術(shù)與工程（石油物探）院（系）地球物理學(xué)院 地球探測(cè)與信息技術(shù)系單項(xiàng)成績(jī)內(nèi)容理解寫(xiě)作結(jié)構(gòu)程序設(shè)計(jì)模型設(shè)計(jì)計(jì)算結(jié)果結(jié)果分析總成績(jī)實(shí)驗(yàn)報(bào)告一、 實(shí)驗(yàn)題目疊加地震記錄的相移波動(dòng)方程正演模擬二、實(shí)驗(yàn)?zāi)康恼莆崭飨蛲越橘|(zhì)任意構(gòu)造、水平層狀速度結(jié)構(gòu)地質(zhì)模型的相移波動(dòng)方程正

2、演模擬基本理論、實(shí)現(xiàn)方法與程序編制，由正演記錄初步分析地震信號(hào)的分辨率。三、原理公式1、地震波傳播的波動(dòng)方程設(shè)（x,z）為空間坐標(biāo)，t為時(shí)間，地震波傳播速度為v(x,z)，則二維介質(zhì)中任意位置、任意時(shí)刻的地震波場(chǎng)為p(z,x,t)：壓縮波縱波。則二維各向同性均勻介質(zhì)中地震波傳播遵循的聲波方程：2p(x,z,t)x2+2p(x,z,t)z2=1v2(x,z)2p(x,z,t)t2 (1)2、傅里葉變換的微分性質(zhì)p(t)與其傅里葉變換的P()的關(guān)系：P=-pte-itdt 正傅里葉變換 pt=12-Peitdt 逆傅里葉變換 2則有時(shí)間微分性質(zhì)：(i)P=-dptdte-itdt 一階微分 (i)

3、2P=-d2ptdt2e-itdt 二階微分 3為頻率，=2T，T為周期。同理有空間微分性質(zhì)：(ik)Pk=-dpxdxe-ikxdx 一階微分 (ik)2Pk=-d2pxdx2e-ikxdx 二階微分 (4)k為波數(shù)，k=2， 為波長(zhǎng)3、地震波傳播的相移外推公式令速度v 不隨x 變化，只隨z 變化，則利用傅里葉變換微分性質(zhì)（3）和（4）式，把波動(dòng)方程(1)式變換到頻率波數(shù)域，得：ik2Pk,zi,+2P(k,z,)z2=i2vz2P(k,z,)或：2P(k,z,)z2=-2vz2-k2Pk,z, (5)令：kz2=2vz2-k2則(5)式的解為： Pk,z,=c1e-ikzz+c2eikzz

4、 (6)包括上行波和下行波兩項(xiàng)，正演模擬取上行波：Pk,z,=c1e-ikzz (7)若Zj和Zj+1間隔為z，速度v(z)為在此間隔內(nèi)不隨Z變的常數(shù)，(7)式實(shí)現(xiàn)波場(chǎng)從Zj+1到Zj的延拓，即：Pk,zj,=c1e-ikzz在深度Zj+1開(kāi)始向上延拓到Zj，若延拓深度為零，即：Z=Zj+1-Zj，則 Pk,zj=zj+1,=c1e-ikz(zj+1-zj)=ce-ikz×0=c (8) 對(duì)于任意深度Zj+1到Zj的延拓，可得正演模擬中地震波的傳播方程（延拓公式Pk,zj,=Pk,zj+1,c1e-ikz(zj+1-zj) (9)4、初始條件和邊界條件按照爆炸界面理論，反射界面震源在

5、 t=0 時(shí)刻同時(shí)起爆，此時(shí)刻的波場(chǎng)就是震源。根據(jù)不同情況，可直接使用反射系數(shù)脈沖或子波作震源。如果直接使用反射系數(shù)作震源脈沖，則初始條件可表示為： p0x,z,t=rx,z t=0 0 t=其他 (10) p0x,z,t對(duì)時(shí)間t和空間x做二維傅立葉變換，則得頻率-波數(shù)域的初始波場(chǎng) p0k,z,。邊界條件： p0x,z,t=rx,z t=0,xmin<x<xmax,zmin<z<zmax 0 t=其他,x=其他,z=其他 (11) 其他參數(shù)都是在xmin<x<xmax,zmin<z<zmax范圍內(nèi)定義的。5、邊界處理（1）邊界反射問(wèn)題把實(shí)際無(wú)窮空

6、間區(qū)域中求解波場(chǎng)的問(wèn)題化為有窮區(qū)域求解時(shí)，左右兩邊使用零邊界條件。物理上假設(shè)探區(qū)距Xmin與Xmax兩個(gè)端點(diǎn)很遠(yuǎn)，在兩個(gè)端點(diǎn)上收到的反射波很弱。但是，上述條件在實(shí)際中不能成立，造成零邊界條件反而成為絕對(duì)阻止波通過(guò)的強(qiáng)反射面。在正演模擬的剖面上出現(xiàn)了邊界假反射干涉正常界面的反射。（2）邊界強(qiáng)反射的處理鑲邊法、削波法、吸收邊界都能有效消除邊界強(qiáng)反射。削波法就是在波場(chǎng)延拓過(guò)程中，每延拓一次，在其兩側(cè)均勻衰減到零，從而消除邊界強(qiáng)反射的影響。假設(shè)橫向總長(zhǎng)度為NX，以兩邊Lx道吸波為例，有以下吸波公式：AbsNx-Ix=AbsIx=sin(2×LxLx-1) 0<=Ix<Lx Abs

7、Ix=1. Ix=其它 (12) 6、數(shù)字化根據(jù)數(shù)字信號(hào)處理的采樣定理，把連續(xù)的信號(hào)變?yōu)橛?jì)算機(jī)能處理的數(shù)字信號(hào)，使相移法正演模擬得以實(shí)現(xiàn)。頻域抽樣定理：一個(gè)頻譜受限信號(hào) f(t)，如果時(shí)間只占據(jù) -tm-tm的范圍，若在頻域以不大于 1/2tm 頻率間隔 f1/2tm 對(duì)信號(hào)f(t)的頻譜F()采樣，則抽樣到的離散信號(hào) F1 可以唯一表示原信號(hào)。時(shí)域抽樣定理：一個(gè)時(shí)間受限信號(hào) f(t)，如果頻譜只占據(jù)-m-m的范圍，則信號(hào) f(t)可以用等間隔的抽樣值唯一表示出來(lái)，而時(shí)間 t 抽樣間隔必須不大于1/2fm，m=2fm，t=1/2fm。四、實(shí)驗(yàn)內(nèi)容削波法相位移正演模擬（1）點(diǎn)繞射構(gòu)造和水平層狀速

8、度模型（參數(shù)如圖1所示）的正演數(shù)值模擬； 1）削波的正演；2）無(wú)削波的正演；（2）計(jì)算中點(diǎn)和兩個(gè)邊界的 信號(hào)位置，分析實(shí)驗(yàn)結(jié)果的正確性；（3）做同樣模型的褶積模型 數(shù)值模擬，對(duì)比分析兩者的異同；圖1 點(diǎn)繞射構(gòu)造與水平速度模型（4）改變繞射點(diǎn)位置、速度，再做正演模擬；五、 方法路線1、參數(shù)初始化；2、形成邊界削波數(shù)據(jù)；3、波場(chǎng)初始化；3、Zmax層波場(chǎng)延拓到深度Zmax-1；5、Zi+1層波場(chǎng)延拓到深度Zi；6、重復(fù)5，從Iz=Nz-1開(kāi)始，直到Iz=1，得測(cè)線上的頻率空間域波場(chǎng)；7、頻率-空間域波場(chǎng)對(duì)頻率做反傅里葉變換，得時(shí)間-空間波場(chǎng)；8、使用Fimage軟件顯示所得結(jié)果。六、實(shí)驗(yàn)結(jié)果1、結(jié)

9、果顯示及分析（點(diǎn)繞射結(jié)構(gòu)）圖2 削波（左）正演模擬結(jié)果與未削波正演模擬結(jié)果（右）由圖2分析可知：削波前，由于在兩個(gè)端點(diǎn)上收到的反射波很弱，造成零邊界條件反而成為絕對(duì)阻止波通過(guò)的強(qiáng)反射面，在正演模擬的剖面上出現(xiàn)了邊界假反射干涉正常界面的反射。削波后，邊界假反射的影響消除了，且曲線變得更光滑了。a、速度V1=5000、V2=5500，深度h=2000時(shí)，改變繞射點(diǎn)x的位置；圖3 x=Nx/4-1時(shí)，反射系數(shù)（左）與削波正演模擬結(jié)果（右）圖4 x=Nx/2-1時(shí)，反射系數(shù)（左）與削波正演模擬結(jié)果（右）由圖3與圖4對(duì)比分析可知：當(dāng)速度，深度一定時(shí)，改變繞射點(diǎn)的位置，曲線隨繞射點(diǎn)位置的改變而左右移動(dòng)，且

10、移動(dòng)趨勢(shì)與繞射點(diǎn)位置的移動(dòng)一致。b、繞射點(diǎn)x=Nx/2-1，深度h=2000 ，速度V2=5500時(shí)，改變速度V1；V2=5500V1=4000圖5 V1=4000時(shí)，速度模型（左）與削波正演模擬結(jié)果（右）V1=6000V2=5500圖6 V1=6000時(shí)，速度模型（左）與削波正演模擬結(jié)果（右）由圖5與圖6對(duì)比分析可知：當(dāng)繞射點(diǎn)位置固定，深度一定時(shí)，隨著介質(zhì)速度增大，曲線變得越平緩。這是由于速度增大時(shí)，波的傳播時(shí)間減小所致。c、速度V1=5000、V2=5500，繞射點(diǎn)x=Nx/2-1，改變深度h；圖7 h=1000時(shí)，反射系數(shù)（左）與削波正演模擬結(jié)果（右）圖8 h=3000時(shí)，反射系數(shù)（左）

11、與削波正演模擬結(jié)果（右）由圖7與圖8對(duì)比分析可知：當(dāng)介質(zhì)速度一定，繞射點(diǎn)位置不變時(shí)，隨著繞射點(diǎn)深度的增加，曲線變得越來(lái)越平滑。顯然是由于深度變大，波的傳播時(shí)間變小所致。2、結(jié)果分析（水平層速度模型）圖9 未削波（左）正演模擬結(jié)果與削波正演模擬結(jié)果（右）由圖9分析可知：削波前，由于在兩個(gè)端點(diǎn)上收到的反射波很弱，造成零邊界條件反而成為絕對(duì)阻止波通過(guò)的強(qiáng)反射面，在正演模擬的剖面上出現(xiàn)了邊界假反射干涉正常界面的反射。削波后，邊界假反射的影響消除了，且曲線變得更光滑了。a、 深度h=2000 ，速度V2=5500時(shí)，改變速度V1；V1=4000V2=5500圖10 V1=4000時(shí)，速度模型（左）與削波

12、正演模擬結(jié)果（右）V2=5500V1=6000圖11 V1=6000時(shí)，速度模型（左）與削波正演模擬結(jié)果（右）由圖10與圖11對(duì)比分析可知：水平層深度不變時(shí)，隨著速度的增加，結(jié)果顯示的水平層向上移動(dòng)。這是由于速度增大時(shí)，波的傳播時(shí)間減小的。b、 速度V1=5000、V2=5500，改變深度h；圖12 h=1000時(shí)，反射系數(shù)（左）與削波正演模擬結(jié)果（右）圖13 h=3000時(shí)，反射系數(shù)（左）與削波正演模擬結(jié)果（右）由圖12與圖13對(duì)比分析可知：當(dāng)波的傳播速度不改變時(shí)，隨著水平層深度的增加，結(jié)果顯示的水平層向下移動(dòng)。這是由于深度增加，波的傳播時(shí)間增大。七、討論建議1、實(shí)驗(yàn)收獲通過(guò)本次實(shí)驗(yàn)，對(duì)疊后

13、相移波動(dòng)方程正演模擬的過(guò)程有了更深入的理解，以及對(duì)其公式和意義都有了進(jìn)一步的認(rèn)識(shí)。編程能力又有了一定提高。2、存在問(wèn)題圖14 x=3*Nx/4-1時(shí)，反射系數(shù)（左）與削波正演模擬結(jié)果（右）如圖14，當(dāng)速度，深度一定時(shí)，繞射點(diǎn)位置x=3*Nx/4-1，其反射系數(shù)不能用Fimage正確顯示，由顯示的削波正演模擬結(jié)果可知，使用削波函數(shù)并未起到削波的作用。3、心得體會(huì)本次實(shí)驗(yàn)，并不容易。在實(shí)驗(yàn)中，必須要理解理論知識(shí)，才能夠正確的編寫(xiě)程序。在實(shí)驗(yàn)過(guò)程中，遇到許多問(wèn)題，通過(guò)和同學(xué)討論，向老師求解，最終才得出實(shí)驗(yàn)結(jié)果。在以后的學(xué)習(xí)工作中，當(dāng)一個(gè)人的力量不夠時(shí)，要積極和其他人求教，共同學(xué)習(xí)，共同進(jìn)步。附程序代

14、碼：#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.

15、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 I

16、x;float AbsNx;if(fp_Abs=fopen("Absb.dat","wb")=NULL)printf("Connot open file ""Absb""");for(Ix=0;Ix<Nx;Ix+)AbsIx=1.0;/1.Absorb Boundary Initializing ？for(Ix=0;Ix<Labs;Ix+)AbsIx=sqrt(sin(Pai/2)*Ix/(Labs-1);AbsNx-Ix-1=AbsIx;/2.Absorb Boundary Com

17、pute？for(Ix=0;Ix<Nx;Ix+)fwrite(&AbsIx,sizeof(AbsIx),Nx,fp_Abs);/3.Byte Number offclose(fp_Abs);return(1);/function03: Form Reflect Structure Modelint Rflct( )FILE *fp_Rflct;int Ix,Iz;float RflctNz;if(fp_Rflct=fopen("Rflct.dat","wb")=NULL)printf("Connot open file "

18、;"Reflection""");for(Ix=0;Ix<Nx;Ix+)for(Iz=0;Iz<Nz;Iz+)RflctIz=0.;if(Ix=Nx/4-1&&Iz=Nz/2-1) RflctIz=1;fwrite(&RflctIz,sizeof(RflctIz),1,fp_Rflct);/4.fclose(fp_Rflct);return(1);/function04: Form Velocity Modelint Vlcty( )FILE *fp_Vlcty;int Ix,Iz;float VlctyNz;if(f

19、p_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+)VlctyIz=5000.;fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);/5. for(Iz=(int)(3*Nz/4);Iz<Nz;Iz+)VlctyIz=5500.;fwrite(

20、&VlctyIz,sizeof(VlctyIz),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 fr4096,fi4096;int k=0;long Ln=0;for(k=0;n-Ln>0;k+)Ln=(long)pow(2,k);k=k-1;for (

21、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;frit = pris;fiit = piis;pr0 = 1.0;pi0 = 0.0;p = 6.283185306/(1.0*n);pr1 = (float) cos(p);pi1 = -(float)sin(p);if (l!=0)pi1=-pi1;for (i=2; i<=n-1; i+)p = pri-1*pr1;q = pii-1*pi1;s = (pri-1+pii-1)*(pr1+pi1);p

22、ri = p-q;pii = s-p-q;for (it=0; it<=n-2; it=it+2)vr = frit;vi = fiit;frit = vr+frit+1;fiit = vi+fiit+1;frit+1 = vr-frit+1;fiit+1 = vi-fiit+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 = prm*j*frit+j+nv/2;q = pim*j

23、*fiit+j+nv/2;s = prm*j+pim*j;s = s*(frit+j+nv/2+fiit+j+nv/2);poddr = p-q;poddi = s-p-q;frit+j+nv/2 = frit+j-poddr;fiit+j+nv/2 = fiit+j-poddi;frit+j = frit+j+poddr;fiit+j = fiit+j+poddi;if(l!=0)for(i=0; i<=n-1; i+)fri = fri/(1.0*n);fii = fii/(1.0*n);if(il!=0)for(i=0; i<=n-1; i+)pri = sqrt(fri*f

24、ri+fii*fii);if(fabs(fri)<0.000001*fabs(fii)if (fii*fri)>0)pii = 90.0;elsepii = -90.0;elsepii = atan(fii/fri)*360.0/6.283185306;for(i=0;i<n;i+)pri=fri;pii=fii;return ( 1 );/function05: Form Inition Wave Field: Reflect Coefficiency /int WvFld0()FILE *fp_Rflct,*fp_Wfld0r,*fp_Wfld0i;int Ix,Iz,I

25、t;float RflctNz;float Wfld0rNt,Wfld0iNt;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", &

26、quot;rb")=NULL)printf("Connot open Rflct.dat");for(Ix=0;Ix<Nx;Ix+)for(Iz=0;Iz<Nz;Iz+)fread(&RflctIz,sizeof(RflctIz),1,fp_Rflct);for(It=0;It<Nt;It+)Wfld0rIt=0.0;/7Wfld0iIt=0.0;/8.Initial Wave field Initializing: Wfld0r0=RflctIz;/9.Initial Wave field Initializingif( kkfft(W

27、fld0r,Wfld0i,Nt,0 ) !=1 )/10.FFT or Inver FFT?printf("FFT is error");exit(0);for(It=0;It<Nt/2+1;It+)/fwrite(&Wfld0rIt,sizeof(Wfld0rIt),1,fp_Wfld0r);/11fwrite(&Wfld0iIt,sizeof(Wfld0iIt),1,fp_Wfld0i);/12 fclose(fp_Rflct);fclose(fp_Wfld0r);fclose(fp_Wfld0i);return(1);/ Function06.1

28、.1: Read In Velocity Data and Absorb Boundary Dataint 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(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty

29、);/19.fclose(fp_Vlcty);if(fp_Absb = fopen("Absb.dat", "rb") =NULL)printf("Connot open Absb.dat" );for(Ix=0;Ix<Nx;Ix+)fread(&AbsbIx,sizeof(AbsbIx),1,fp_Absb);/20.fclose(fp_Absb);remove("Absb.dat");return(1);/ 06.1.2.1: Read Data From (Ix,Iz,Iw)Oder to (I

30、w,Iz,Ix)Oderint ReadIxIzIwToIwIzIx(FILE *fp_Wfld0r,FILE *fp_Wfld0i,float Wfld0r,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(Wfld0rIx)*Mgrtn,0);/26fread(&Wfld0rIx,sizeof(Wfld0rIx),1,fp_Wfld0r); /27fseek(fp_Wfl

31、d0i, sizeof(Wfld0iIx)*Mgrtn,0);/28fread(&Wfld0iIx,sizeof(Wfld0iIx),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 W

32、fld0i,int Iz,int Iw);int Ix,Nw=Nt;float Wfld0rNx,Wfld0iNx;if(ReadIxIzIwToIwIzIx(fp_Wfld0r,fp_Wfld0i,Wfld0r,Wfld0i,Iz,Iw)!=1)printf("exp_ikzDz is error");exit(0);for(Ix=0;Ix<Nx;Ix+)WfldrIx+=Wfld0rIx;/21.WfldiIx+=Wfld0iIx;/22. WfldrIx=WfldrIx*AbsbIx;/23. WfldiIx=WfldiIx*AbsbIx;/24.return(

33、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;eikzdz0 = 0.;eikzdz1 = 0.;kz =sqrt(pow(Iw*Dw/Vc,2)-pow(Ix*Dkx,2);/38. if(kz >0.0)/39.eikzdz0 = (float)cos(kz*Dz);/40. eikzdz1 = (float)sin(-1)*kz*Dz);/41.return(1);/ 06.1.3: Ext

34、rapolate One Depth Stepint MoveOneDz(float Wfldr,float Wfldi,float Vz,float Dkx,float Dw,int Iw)int Ikx,Nkx=Nx;float kz2,Wfld_r,Wfld_i;if(kkfft(Wfldr,Wfldi,Nx,0)!=1 ) /30printf("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("e

35、xp_ikzDz is error");exit(0);Wfld_r =kz0*WfldrIkx-kz1*WfldiIkx;/32Wfld_i =kz0*WfldiIkx+kz1*WfldrIkx;/33 WfldrIkx = Wfld_r;WfldiIkx = Wfld_i;if(Ikx!=0)&&(Ikx!=Nkx/2)/34Wfld_r =kz0*WfldrNkx-Ikx-kz1*WfldiNkx-Ikx;/35.Wfld_i =kz0*WfldiNkx-Ikx+kz1*WfldrNkx-Ikx;/36. WfldrNkx-Ikx = Wfld_r;WfldiN

36、kx-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,flo

37、at,int);FILE *fp_Wfldr,*fp_Wfldi,*fp_Wfld0r,*fp_Wfld0i;int Ix,Iz,Iw,Nw=Nt;float VlctyNz,AbsbNx;float WfldrNx,WfldiNx;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.da

38、t","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_W

39、fldi = fopen("Wfldi.dat", "wb") =NULL) printf("Connot open Wfldi.dat" );for(Iw=0;Iw<Nw/2+1;Iw+)printf("Iw=%dn",Iw);for(Ix=0;Ix<Nx;Ix+)WfldrIx=0.0;/15. WfldiIx=0.0;/16for(Iz=Nz-1;Iz>0;Iz-)if( FrmNewWfld(fp_Wfld0r,fp_Wfld0i,Wfldr,Wfldi,Absb,Iz,Iw)!=1)pr

40、intf("FrmNewWfld Error");exit(0);if(MoveOneDz(Wfldr,Wfldi,(float)(VlctyIz/2.),Dkx,Dw,Iw)!=1)printf("WfldMoveOneStep Error");exit(0);for(Ix=0;Ix<Nx;Ix+)fwrite(&WfldrIx,sizeof(WfldrIx),1,fp_Wfldr);/17 fwrite(&WfldiIx,sizeof(WfldiIx),1,fp_Wfldi);/18fclose(fp_Wfld0r); remo

41、ve("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 WfldtrNt,WfldtiN

42、t;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;I