疊前地震記錄的相移波動(dòng)方程正演模擬實(shí)驗(yàn)

1、本科生實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)課程 油氣勘探新方法 學(xué)院名稱(chēng) 地球物理學(xué)院 專(zhuān)業(yè)名稱(chēng) 勘查技術(shù)與工程（石油物探） 學(xué)生姓名 學(xué)生學(xué)號(hào) 指導(dǎo)教師 熊高君 實(shí)驗(yàn)地點(diǎn) 5417 實(shí)驗(yàn)成績(jī) 2015年12月 成都理工大學(xué)油氣勘探新方法實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)時(shí)間2015年 12月開(kāi)課單位地球物理學(xué)院指導(dǎo)教師熊高君實(shí)驗(yàn)題目：疊前地震記錄的相移波動(dòng)方程正演模擬實(shí)驗(yàn)姓名學(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)二、實(shí)驗(yàn)?zāi)康恼莆崭飨蛲越橘|(zhì)任意構(gòu)造、水平層狀速度結(jié)構(gòu)地質(zhì)模型的

2、疊前地震記錄相移波動(dòng)方程正演模擬基本理論、實(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=-dptd

3、te-itdt 一階微分 (i)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=i2vz2Pk,z, （5）或： 2P(k,z,)z2=-2vz2-k2Pk,z, (6)令：kz2=2vz2-k2則(6)式的解為： Pk

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

5、1)4、初始條件和邊界條件按照爆炸界面理論，反射界面震源在 t=0 時(shí)刻同時(shí)起爆，此時(shí)刻的波場(chǎng)就是震源。根據(jù)不同情況，可直接使用反射系數(shù)脈沖或子波作震源。如果直接使用反射系數(shù)作震源脈沖，則初始條件可表示為： p0x,z,t=rx,z t=0 0 t=其他 (12)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=其他 (13) 其他參數(shù)都是在xmin<x<xmax,zmin<z<zmax范

6、圍內(nèi)定義的。5、邊界處理（1）邊界反射問(wèn)題把實(shí)際無(wú)窮空間區(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，以?xún)蛇匧x道吸波為例，有以下吸波公式： AbsNx-Ix=AbsIx=sin(2

7、5;LxLx-1) 0<=Ix<Lx AbsIx=1. Ix=其它 (14)6、定位原理某一震源在某一固定檢波點(diǎn)產(chǎn)生的記錄， 就等于該固定檢波點(diǎn)主動(dòng)到所有可能反射點(diǎn)處接收的由這一震源激發(fā)的反射波的總和。形象地稱(chēng)為檢波點(diǎn)下延記錄原理。7、反射記錄的形成 1）待激發(fā)的二次源的描述：介質(zhì)的離散化把整個(gè)介質(zhì)空間離散化，每個(gè)離散點(diǎn)都可以看成一個(gè)繞射點(diǎn)，某深度的反射界面與這一深度線的交點(diǎn)為這一反射界面的分布點(diǎn)。 每個(gè)離散點(diǎn)就是待激發(fā)的二次源 2）反射界面的分布某一深度層的反射界面分布在這一深度的離散點(diǎn)上，即反射界面與這一深度線的交點(diǎn)為這一反射界面的分布點(diǎn)。 3）反射界面二次源的形成 ：繞射波的

8、產(chǎn)生震源從炸點(diǎn)開(kāi)始向下傳播，遍歷介質(zhì)的每一個(gè)離散點(diǎn)繞射點(diǎn)。每遇到一個(gè)繞射點(diǎn)，都會(huì)激發(fā)該點(diǎn)的繞射波。反射界面與某深度離散點(diǎn)的 交叉點(diǎn)的反射系數(shù)不為零， 產(chǎn)生的繞射波也不為零， 而其他沒(méi)有反射界面的地方，反 反射系數(shù)為零，是均勻介質(zhì)， 產(chǎn)生的繞射波為零。 4）反射記錄的形成把檢波器的檢波機(jī)制看成一個(gè)數(shù)學(xué)算子，或檢波因子，可用脈沖或雷克子波。在延拓震源的同時(shí)，把檢波器即檢波因子也做完全相同的向下延拓，同樣遍歷每個(gè)離散點(diǎn)。震源和檢波因子每延拓一個(gè)深度步長(zhǎng)，震源、檢波因子及反射系數(shù)三者做相關(guān)運(yùn)算，則檢波點(diǎn)就記錄到了該深度每個(gè)繞射點(diǎn)的繞射波。所有二次源地震波都傳播到測(cè)線某檢波點(diǎn)，疊加而成該點(diǎn)的記錄炮集中的

9、某一記錄道。測(cè)線排列多少道檢波器，則都可記錄道所在位置的地震模擬記錄。8、數(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)容削波法相位移

10、正演模擬（1）點(diǎn)繞射構(gòu)造和水平層狀速度模型（參數(shù)如圖1所示）的正演數(shù)值模擬； 1）削波的正演；2）無(wú)削波的正演；（2）計(jì)算中點(diǎn)和兩個(gè)邊界的信號(hào)位置，分析實(shí)驗(yàn)結(jié)果的正確性；（3）做同樣模型的褶積模型數(shù)值模擬，對(duì)比分析兩者的異同；（4）改變繞射點(diǎn)位置、速度，再做正演模擬；圖1 點(diǎn)繞射構(gòu)造與水平速度模型五、 方法路線1、參數(shù)初始化；2、二次源的形成：1）震源波場(chǎng)定義：Scr(x,z,t) Scrx,z,t=1 x=x0,z=z0,t=00 其它 （15）2）界面待激發(fā)爆炸源定義：Rsht(x,z,t) Rshtx,z,t=Rflctx,z t=00 其它 （16）3）形成削波數(shù)據(jù)：Absb( x )

11、在所有深度都用同樣的削波數(shù)據(jù)： AbsNx-Ix=AbsIx=sin(2×LxLx-1) 0<=Ix<Lx AbsIx=1. Ix=其它 （17）4）從測(cè)線深度開(kāi)始向下做波場(chǎng)延拓 在所有深度都用同樣的削波數(shù)據(jù)震源波場(chǎng)做削波處理 Scrx,zj,t=Scrx,zj,t×Absbx （18）震源做 x 和t 的傅里葉變換到頻率-波數(shù)域： Scrx,zj,tFFTScrx,zj,t （19）計(jì)算相移延拓算子： eikzdzkx,zj,=e-i2vzj+12-kx2(zj+1-zj) （20）震源波場(chǎng)延拓計(jì)算 Scrkx,zj+1,= Scrkx,zj,×ei

12、kzdzkx,zj, （21） 延拓后的震源波場(chǎng)變換到時(shí)間-空間域： Scrkx,zj+1,FFT-1 Scrx,zj,t （22）激發(fā)當(dāng)前深度的二次源：震源波場(chǎng)與待激發(fā)波場(chǎng)相關(guān) Rscndx,zj+1,t=-Scrx,zj+1,Rshtx,zj+1,t+d （23）將二次源變換到頻率域存儲(chǔ)： Rscndx,zj+1,tFFTRscndx,zj+1, 24從目前深度開(kāi)始，重到步，震源波場(chǎng)繼續(xù)向下延拓,直至介質(zhì)最深處，生成介質(zhì)所有點(diǎn)的繞射波- 整個(gè)介質(zhì)的二次源形成。3、正演記錄產(chǎn)生的數(shù)學(xué)過(guò)程 1）定義頻率域延拓波場(chǎng):wfld(x,z,)；2）調(diào)入削波數(shù)據(jù)：Absb(x)；3）在頻率域延拓波場(chǎng)：

13、每給一個(gè)頻率j，延拓波場(chǎng)初始化：wfldx,z,j=0；從最深處開(kāi)始，波場(chǎng)向上延拓： a、取出當(dāng)前頻率和深度的二次源：Rscnd(x,zj,j)；b、形成新延拓波場(chǎng)： wfldx,zj,j=wfldx,zj,j+Rscndx,zj,j （25）c、新延拓波場(chǎng)做削波處理 wfldx,zj,j=wfldx,zj,j×Absbx （26）d、新延拓波場(chǎng)變換到波數(shù)域： wfldx,zj,jFFTwfldkx,zj,j （27）e、計(jì)算相移延拓算子： eikzdzkx,zj,j=e-ij2vzj2-kx2(zj-zj-1) （28）f、波場(chǎng)延拓計(jì)算 wfldkx,zj-1,j=wfldkx,z

14、j,j×eikzdzkx,zj,j （29）g、波場(chǎng)變換到空間域： wlfdkx,zj-1,jFFT-1wlfdx,zj-1,j （30）h、重復(fù)a到g的計(jì)算，直至波場(chǎng)延拓到測(cè)線深度，ZJ=Z1，得到測(cè)線上的頻率域波場(chǎng)：wlfdx,z1,j，存儲(chǔ)測(cè)線上的波場(chǎng):wlfdx,z1,j;重復(fù)到的計(jì)算，直到完成所有頻率的循環(huán)；4)把波場(chǎng)變換到時(shí)間域，得一炮的疊前正演模擬結(jié)果： wlfdx,z1,FFT-1wlfdx,z1,t （31）4、頻率-空間域波場(chǎng)對(duì)頻率做反傅里葉變換，得時(shí)間-空間波場(chǎng)；5、使用Fimage軟件顯示所得結(jié)果。六、實(shí)驗(yàn)結(jié)果1、結(jié)果顯示及分析（點(diǎn)繞射結(jié)構(gòu)）圖2 削波（左）正

15、演模擬結(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)，且移動(dòng)趨勢(shì)與繞射點(diǎn)位置的移動(dòng)一致。b、繞射點(diǎn)x

16、=Nx/2-1，深度h=2000 ，速度V2=5500時(shí)，改變速度V1；V1=4000V2=5500圖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ù)（左）與削波正演模擬結(jié)果（右）由圖7與圖8對(duì)比分析

17、可知：當(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í)，速度模型（左）與削波正演模擬結(jié)果（右）V2=5500V1=600

18、0圖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ì)疊后相移波動(dòng)方程正演模擬的過(guò)程有了更深入的理解，

19、以及對(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é)果。附程序代碼：/Include #include <stdlib.h> #include <conio.h> #

20、include <math.h> #include <stdio.h> #include <string.h> /Parameters set #define Nx 128 /Trace Number #define Nt 256 /Record Number #define Nz 100 /Depth Number #define Labs 15 /Length Of Boundary Absorbing #define Dx 20 /Trace Interval #define Dt 0.004 /Record Interval #define Dz 2

21、0 /Depth Interval #define Pai 3.1415926 /function:Judge the power of 2 int 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); /function:Form Absorb Bounderyint Absorb( ) FILE *fp_Abs; int Ix; double AbsNx; if(fp

22、_Abs=fopen("Absb.dat","wb")=NULL)printf("Connot open file ""Absb"""); for(Ix=0;Ix<Nx;Ix+) AbsIx=1.; for(Ix=0;Ix<=Labs-1;Ix+) AbsIx=sqrt(sin(Pai/2)*(Ix/(Labs-1);AbsNx-Ix-1=AbsIx; for(Ix=0;Ix<Nx;Ix+) fwrite(&AbsIx,sizeof(AbsIx),1,fp_Abs);

23、fclose(fp_Abs); return(1); /function:Form Reflect Structure Model int Rflct( ) FILE *fp_Rflct; int Ix,Iz; float RflctNz; 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+)

24、 RflctIz=0.; if(Iz=Nz/5&&Ix<Nx/2) RflctIz=1.; if(Iz=2*Nz/5&&Ix>=Nx/2) RflctIz=1.; /if(Iz=3*Nz/5&&fmod(Ix+1,16)=0.) RflctIz=5.; /if(Iz=4*Nz/5&&fmod(Ix+1,32)=0.) RflctIz=7.; fwrite(&RflctIz,sizeof(RflctIz),1,fp_Rflct); fclose(fp_Rflct); return(1); /function:Rs

25、ht(x,z,t)int Rsht()int Ix,Iz,It;float RshtNt,RflctNz;FILE *fp_Rflct=fopen("Rflct.dat","rb");FILE *fp_Rsht=fopen("Rsht.dat","wb");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+)RshtIt=0.;Rsht0=Rf

26、lctIz;for(It=0;It<Nt;It+)fwrite(&RshtIt,sizeof(RshtIt),1,fp_Rsht);fclose(fp_Rflct);fclose(fp_Rsht);return(1); /function:Form Velocity Modelint Vlcty( ) FILE *fp_Vlcty; int Ix,Iz; float VlctyNz; if(fp_Vlcty=fopen("Vlcty.dat","wb")=NULL) printf("Connot open file "&

27、quot;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); for(Iz=(int)(3*Nz/4);Iz<Nz;Iz+) VlctyIz=5500.; fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty); fclose(fp_Vlcty); return(1); /function：Fourier Tr

28、ansform: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 (it=0; it<=n-1; it+)m = it; is = 0; for(i=0; i<=k-1; i+)j = m/2; is = 2

29、*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); pri = p-q; pii = s-p-q;for (it=0; it<=n-2; it=it+2)vr = fr

30、it; 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*fiit+j+nv/2; s = prm*j+pim*j; s = s*(frit+j+nv/

31、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*fri+fii*fii); if(fabs(fri)<0.000001

32、*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 );/function:Form Inition Wave Field:ReflectCoefficiency int Scr0()/Function05.1:kkfft(Wfld0r,Wfld0i,Nt,0) FILE *fp_Scr0r,*fp_Scr0i; int Ix,Iz,It; float Scr

33、0rNt,Scr0iNt; if(fp_Scr0r=fopen("Scr0r.dat","wb")=NULL)printf("Connot open Wfld0r.dat"); if(fp_Scr0i=fopen("Scr0i.dat","wb")=NULL)printf("Connot open Wfld0i.dat"); /if(fp_Rflct=fopen("Rflct.dat", "rb")=NULL) printf("

34、;Connot open Rflct.dat"); for(Ix=0;Ix<Nx;Ix+) printf("Scr0_FFT:Ix=%dn",Ix); for(Iz=0;Iz<Nz;Iz+) for(It=0;It<Nt;It+) Scr0rIt=0; Scr0iIt=0; if(Ix=64&&Iz=0)Scr0r0=1.; if(kkfft(Scr0r,Scr0i,Nt,0)!=1) printf("FFT is error"); exit(0); for(It=0;It<Nt/2+1;It+) fwr

35、ite(&Scr0rIt,sizeof(Scr0iIt),1,fp_Scr0r); fwrite(&Scr0iIt,sizeof(Scr0iIt),1,fp_Scr0i);fclose(fp_Scr0r); fclose(fp_Scr0i); 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(

36、"Vlcty.dat","rb")=NULL)printf("Connot open Rflct.dat" );for(Iz=0;Iz<Nz;Iz+)fread(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);fclose(fp_Vlcty); if(fp_Absb=fopen("Absb.dat","rb")=NULL)printf("Connot open Absb.dat"); for(Ix=0;Ix<Nx;Ix+)frea

37、d(&AbsbIx,sizeof(AbsbIx),1,fp_Absb);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_Scr0r,FILE *fp_Scr0i,float Scr0r,float Scr0i,int Iz,int Iw)int Ix,Nw=Nt; long AddfrmStrt; for(Ix=0;Ix<Nx;Ix+)/Da

38、ta Numbers From start Position To Current Position AddfrmStrt=(Ix*Nz+Iz)*(Nt/2+1)+Iw; /Byte Numbers From start Position To Current Position fseek(fp_Scr0r, sizeof(Scr0rIx)*AddfrmStrt,0); fread(&Scr0rIx,sizeof(Scr0rIx),1,fp_Scr0r); fseek(fp_Scr0i, sizeof(Scr0iIx)*AddfrmStrt,0); fread(&Scr0iIx

39、,sizeof(Scr0iIx),1,fp_Scr0i);return(1);/function06.1.2:Form New Wave Field int FrmNewWfld(FILE *fp_Scr0r,FILE *fp_Scr0i,float Scrr,float Scri,float Absb,int Iz,int Iw) /function06.1.2.1 ReadIxIzIwToIwIzIx(.) int ReadIxIzIwToIwIzIx(FILE *fp_Scr0r,FILE *fp_Scr0i,float Scr0r,float Scr0i,int Iz,int Iw);

40、 int Ix,Nw=Nt; float Scr0rNx,Scr0iNx; /Take out Initial Wave Field Data With Individual Depth if(ReadIxIzIwToIwIzIx(fp_Scr0r,fp_Scr0i,Scr0r,Scr0i,Iz,Iw)!=1)printf("exp_ikzDz is error");exit(0); /Form New Wave Field for(Ix=0;Ix<Nx;Ix+) /Compute Current New Wave Field ScrrIx=ScrrIx+Scr0rI

41、x; ScriIx=ScriIx+Scr0iIx; /Boundary Obsorb of New Wave Fied ScrrIx=ScrrIx*AbsbIx; ScriIx=ScriIx*AbsbIx; return(1);/function06.1.3.1:Compute out PhaseShift Data int exp_ikzDz(float eikzdz,int Ix,float Vc,int Iw,float Dw,float Dkx) float kz=0;eikzdz0=0.;eikzdz1=0.; kz=(Iw*Dw*Iw*Dw)/(Vc*Vc)-Ix*Ix*Dkx*D

42、kx; if(kz>0) kz=sqrt(kz); eikzdz0=(float)cos(-kz*Dz); eikzdz1=(float)sin(-kz*Dz); return(1); /function06.1.3: Extrapolate One Depth Stepint MoveOneDz(float Scrr,float Scri,float Vz,float Dkx,float Dw,int Iw)/06.1.3.1:exp_ikzDz(kz,Ix,Vz,Iw,Dw,Dkx) /06.1.3.2:kkfft(Wfldr, Wfldi,Nx,j) int Ikx,Nkx=Nx;

43、 float kz2,Scr_r,Scr_i; if(kkfft(Scrr,Scri,Nx,0)!=1) printf("FFT is error");exit(0); for(Ikx=0;Ikx<Nx/2+1;Ikx+) /4.2.3.1 Computing Phaseshift Function if(exp_ikzDz(kz,Ikx,Vz,Iw,Dw,Dkx)!=1)printf("exp_ikzDz is error");exit(0); /4.2.3.2 WaveField multiply Phaseshift Function /Co

44、mpute WaveField Phaseshift Scr_r=ScrrIkx*kz0-ScriIkx*kz1; Scr_i=ScriIkx*kz0+ScrrIkx*kz1; ScrrIkx=Scr_r; ScriIkx=Scr_i; if(Ikx!=0&&Ikx!=Nkx/2) Scr_r=kz0*ScrrNkx-Ikx-kz1*ScriNkx-Ikx; Scr_i=kz1*ScrrNkx-Ikx+kz0*ScriNkx-Ikx; ScrrNkx-Ikx=Scr_r; ScriNkx-Ikx=Scr_i; /4.2.4 IFFT of New Wave field From

45、 WaveNumber to Space if(kkfft(Scrr,Scri,Nkx,1)!=1)/37.FFT or Inverse FFT printf("FFT is error");exit(0); return(1);/function06.1:PhaseShift calculate int PhaseShift( ) int ReadVlctyAbsb(float *,float *); int FrmNewWfld(FILE *,FILE *,float *,float *,float *,int,int); int MoveOneDz (float *,float *,float,float,float,int); /1.2 Define Varibles FILE *fp_Scrr,*fp_Scri,*fp_Scr0r,*fp_Scr0i; int Ix,Iz,I