聲波方程數(shù)值模擬實(shí)驗(yàn)報(bào)告.doc

聲波方程數(shù)值模擬實(shí)驗(yàn)報(bào)告學(xué)號(hào)：201107010230姓名：包靈指導(dǎo)老師：程冰潔老師完成時(shí)間：2014年11月26日星期三實(shí)驗(yàn)要求：1、應(yīng)用聲波方程作為正演模擬的波動(dòng)方程；2、將所提供震源函數(shù)離散后繪圖；3、給定兩個(gè)二維速度-深度模型（一個(gè)小模型；一個(gè)大模型），繪出圖形來(lái)；4、對(duì)于小模型，整個(gè)區(qū)域的速度值可設(shè)為常數(shù)，即只有一種介質(zhì)，將震源點(diǎn)放在模型中間，分別記錄兩個(gè)時(shí)刻的波前快照（即該時(shí)刻區(qū)域內(nèi)所有網(wǎng)格點(diǎn)的波場(chǎng)值）。第一時(shí)刻為地震波還未傳播到邊界上的某時(shí)刻，第二時(shí)刻為地震波已經(jīng)傳播到邊界上的某時(shí)刻，體會(huì)其人工邊界反射；5、對(duì)于大模型，定義為水平層狀速度模型（至少兩層）；做兩個(gè)實(shí)驗(yàn)，一是將震源點(diǎn)放在區(qū)域表層任一點(diǎn)，記錄下某些時(shí)刻的波前快照，體會(huì)地震波在兩種介質(zhì)的分界面上傳播規(guī)律；二是合成一個(gè)地震記錄，即記錄下與震源同一深度點(diǎn)的各點(diǎn)所有時(shí)刻的波場(chǎng)值，并指出記錄上的同向軸分別對(duì)應(yīng)哪些波。實(shí)驗(yàn)?zāi)康模?. 通過(guò)本次作業(yè)，加深對(duì)波動(dòng)方程的理解，明白波動(dòng)方程所代表的物理意義。2. 通過(guò)模擬地震波在介質(zhì)中的傳播，理解實(shí)際勘探中地震波在地層中的傳播規(guī)律。3. 通過(guò)模擬水平層狀速度模型，體會(huì)地震波在兩種介質(zhì)分界面的傳播規(guī)律，并能夠從地震記錄中識(shí)別出反射波，透射波，多次波，折射波和繞射波。4. 通過(guò)模擬人工合成的地震記錄，體會(huì)地震勘探基本原理和方法，驗(yàn)證地震波傳播能量波形變化趨勢(shì)。需要的已知條件包括：1）震源函數(shù)2）地層速度（波速）3）邊界條件2彈性波方程：聲波方程的有限差分法數(shù)值模擬對(duì)于二維速度-深度模型，地下介質(zhì)中地震波的傳播規(guī)律可以近似地用聲波方程描述： （4-1）是介質(zhì)在點(diǎn)（x , z）處的縱波速度，為描述速度位或者壓力的波場(chǎng)，為震源函數(shù)。為求式（4-1）的數(shù)值解，必須將此式離散化，即用有限差分來(lái)逼近導(dǎo)數(shù)，用差商代替微商。為此，先把空間模型網(wǎng)格化（如圖4-1所示）。設(shè)x、z方向的網(wǎng)格間隔長(zhǎng)度為，為時(shí)間采樣步長(zhǎng)，則有： （i為正整數(shù)） (j為正整數(shù)) (n為正整數(shù))表示在(i,j)點(diǎn)，k時(shí)刻的波場(chǎng)值。將在(i,j)點(diǎn)k時(shí)刻用Taylor展式展開： (4-2)將在(i,j)點(diǎn)k時(shí)刻用Taylor展式展開：(4-3)將上兩式相加，略去高階小量，整理得(i,j)點(diǎn)k時(shí)刻的二階時(shí)間微商為：(4-4)對(duì)于空間微分，采用四階精度差分格式，（以X方向?yàn)槔┘磳?、分別在(i,j)點(diǎn)k時(shí)刻展開到四階小量，消除四階小量并解出二階微分得： (4-5)同理可得： (4-6)這就實(shí)現(xiàn)了用網(wǎng)格點(diǎn)波場(chǎng)值的差商代替了偏微分方程的微商，將上三個(gè)式子代入(4-1)式中得： (4-7)式中為介質(zhì)速度的空間離散值，是空間離散步長(zhǎng)，為時(shí)間離散步長(zhǎng)，為震源函數(shù)，關(guān)于一般使用一個(gè)理論的雷克型子波代替，即：（1） （4-8）（2） （4-9）上式中，為時(shí)間， 為中心頻率，一般取為20-40HZ，為控制頻帶寬度的參數(shù)，一般取3-5。在實(shí)際計(jì)算過(guò)程中，需把此震源函數(shù)離散，參與波場(chǎng)計(jì)算。確定震源位置。穩(wěn)定性條件：(4-10)這里表示的是地下介質(zhì)的最大波速；若地下介質(zhì)網(wǎng)格間隔、最大速度及時(shí)間采樣間隔不符合(4-8)式時(shí)，遞推求解（4-7）式，波場(chǎng)值會(huì)出現(xiàn)誤差（高階小量）累積，出現(xiàn)不穩(wěn)定現(xiàn)象。頻散關(guān)系式： (4-11)式中為最小速度，為Nyquist（奈奎斯特）頻率。一般取震源子波中的主頻的2倍值參與計(jì)算，G為每個(gè)波長(zhǎng)所占的網(wǎng)格點(diǎn)數(shù)，對(duì)于空間二階差分、時(shí)間二階差分取8，而對(duì)于空間為四階差分的情況則取4方能有效減少頻散。同理：對(duì)于空間微分，采用二階精度差分格式為 （4-12）考慮Reynolds邊界條件（詳細(xì)推到見(jiàn)文獻(xiàn)“基于Marmousi模型的聲波方程有限差分正演算法”），差分格式如下： （4-13）注：其中參數(shù)的設(shè)置：借用以前實(shí)驗(yàn)的數(shù)據(jù)，當(dāng)然可以根據(jù)限制條件4-10、4-11計(jì)算得到；至于震源放于101*5,101*5的矩形中心，根據(jù)速度與位移可以計(jì)算傳播到邊界時(shí)的時(shí)間，此處忽略。二實(shí)驗(yàn)步驟1、 應(yīng)用聲波方程作為正演模擬的波動(dòng)方程，忽略轉(zhuǎn)換波的產(chǎn)生、傳播；2、 將所提供震源函數(shù)離散后繪圖；震源函數(shù)為雷克子波，離散繪圖如下：fm=30;r=3;t=0.002;for n=1:50w(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n);endsubplot(211)plot(w)t1=0.002;for n=1:50w1(n)=(1-2*(pi*fm*(t1*n-1/fm)2)*exp(-(pi*fm*(t1*n-1/fm)2);endsubplot(212)plot(w1)3、 對(duì)于小模型，整個(gè)區(qū)域的速度值可設(shè)為常數(shù)，即只有一種介質(zhì)，將震源點(diǎn)放在模型中間，分別記錄兩個(gè)時(shí)刻的波前快照（即區(qū)域內(nèi)所有網(wǎng)格點(diǎn)的波場(chǎng)值）。第一時(shí)刻為地震波還未傳播到邊界上的某時(shí)刻。由穩(wěn)定條件，設(shè)v為2000，可以令DT=0.001,DH=5,此時(shí)精度較高，且滿足頻散關(guān)系程序，圖形如下：#include #include #include #define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /不能直接初值為0 float vXNZN,wKN,uu0,uu1,uu2;for(i=0;iXN;i+) /定義f函數(shù)，當(dāng)且僅當(dāng)i,j同時(shí)為50時(shí)，f為1，其余為0for(j=0;jZN;j+)if(i=50&j=50)fij=1;else fij=0; for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0;u2ij=0.0;u3ij=0.0;u4ij=0.0;vij=2000; /速度相同表示同一介質(zhì) if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(k=0;kKN;k+) if(0=flag2) /時(shí)間二階差分格式、空間四階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; else if(1=flag2) /時(shí)間二階差分格式、空間二階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=u2i+1j-2*u2ij+u2i-1j; uu2=u2ij+1-2*u2ij+u2ij-1; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(m=0;mXN;m+) for(n=0;nZN;n+) u1mn=u2mn; u2mn=u3mn; if(k=100) for(m=0;mXN;m+) for(n=0;nZN;n+) u4mn=u3mn;/記錄波前快照,中間點(diǎn) if(0=flag1 & 0=flag2)fp=fopen(wavefront14.dat,w);else if(1=flag1 & 0=flag2)fp=fopen(wavefront24.dat,w);else if(0=flag1 & 1=flag2)fp=fopen(wavefront12.dat,w);else if(1=flag1 & 1=flag2)fp=fopen(wavefront22.dat,w); if(fp!=NULL) for(i=0;iXN;i+)for(j=0;jZN;j+)fprintf(fp,%d %d %fn,i,j,u4ij);fclose(fp); void main()function(0,0); /時(shí)間二階差分格式、空間四階差分格式、震源1function(0,1); /時(shí)間二階差分格式、空間二階差分格式、震源1function(1,0); /時(shí)間二階差分格式、空間四階差分格式、震源2function(1,1); /時(shí)間二階差分格式、空間二階差分格式、震源2（1）時(shí)間二階差分格式、空間四階差分格式、震源1（2）時(shí)間二階差分格式、空間二階差分格式、震源1（3）時(shí)間二階差分格式、空間四階差分格式、震源2（4）時(shí)間二階差分格式、空間二階差分格式、震源2第二時(shí)刻為地震波已經(jīng)傳播到邊界上的某時(shí)刻，體會(huì)其人工邊界反射（此處用Reynolds邊界條件）；程序圖形如下：#include #include #include #define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /不能直接初值為0 float vXNZN,wKN,uu0,uu1,uu2,uu3; for(i=0;iXN;i+) /定義f函數(shù)，當(dāng)且僅當(dāng)i,j同時(shí)為50時(shí)，f為1，其余為0for(j=0;jZN;j+)if(i=50&j=50)fij=1;else fij=0; for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0;u2ij=0.0;u3ij=0.0;u4ij=0.0;vij=2000; /速度相同表示同一介質(zhì) if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(k=0;kKN;k+)if(0=flag2) /時(shí)間二階差分格式、空間四階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; else if(1=flag2) /時(shí)間二階差分格式、空間二階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=u2i+1j-2*u2ij+u2i-1j; uu2=u2ij+1-2*u2ij+u2ij-1; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(m=0;mXN;m+) /波向前傳播記錄連續(xù)三個(gè)時(shí)刻的值for(n=0;nZN;n+)u1mn=u2mn; u2mn=u3mn; /Reynolds邊界條件uu3=vij*(DT/DH);u30j=u20j+u21j-u11j+uu3*(u21j-u20j-u12j+u11j);/左邊界u3XNj=u2XNj+u2XN-1j-u1XN-1j+uu3*(u2XN-1j-u2XNj-u1XN-2j+u1XN-1j);/右邊界u3i0=u2i0+u2i1-u1i1+uu3*(u2i1-u2i0-u1i2+u1i1);/下邊界u3iZN=u2iZN+u2iZN-1-u1iZN-1+uu3*(u2iZN-1-u2iZN-u1iZN-2+u1iZN-1);/上邊界if(k=180) /只需改動(dòng)K值，即顯示邊界反射for(m=0;mXN;m+)for(n=0;nZN;n+)u4mn=u3mn;/記錄波前快照，邊界 if(0=flag1 & 0=flag2)fp=fopen(wavefront14.dat,w);else if(1=flag1 & 0=flag2)fp=fopen(wavefront24.dat,w);else if(0=flag1 & 1=flag2)fp=fopen(wavefront12.dat,w);else if(1=flag1 & 1=flag2)fp=fopen(wavefront22.dat,w); if(fp!=NULL) for(i=0;iXN;i+)for(j=0;jZN;j+)fprintf(fp,%d %d %fn,i,j,u4ij);fclose(fp); void main()function(0,0); /時(shí)間二階差分格式、空間四階差分格式、震源1function(0,1); /時(shí)間二階差分格式、空間二階差分格式、震源1function(1,0); /時(shí)間二階差分格式、空間四階差分格式、震源2function(1,1); /時(shí)間二階差分格式、空間二階差分格式、震源2（1）時(shí)間二階差分格式、空間四階差分格式、震源1（2）時(shí)間二階差分格式、空間二階差分格式、震源1（3）時(shí)間二階差分格式、空間四階差分格式、震源2（4）時(shí)間二階差分格式、空間二階差分格式、震源24、 對(duì)于大模型，定義為水平層狀速度模型；做兩個(gè)實(shí)驗(yàn)，一是將震源點(diǎn)放在區(qū)域表層任一點(diǎn)，記錄下某些時(shí)刻的波前快照，體會(huì)地震波在兩種介質(zhì)的分界面上傳播規(guī)律，指出哪是反射波，哪是透射波；這時(shí)取小模型的常量，為減少頻散，速度v至少為2400程序圖形如下：#include #include #include #define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 200#define ZN 100#define DH 5#define DT 0.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN; /不能直接初值為0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2;if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0; u2ij=0.0; u3ij=0.0; u4ij=0.0; fij=0.0; if (j=30) vij=2400; /第一層速度為2400 else vij=3000; /第二層速度為3000 f10010=1; /定義f函數(shù)，在100，10為1 for(k=0;kKN;k+) if(0=flag2) /時(shí)間二階差分格式、空間四階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; else if(1=flag2) /時(shí)間二階差分格式、空間二階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=u2i+1j-2*u2ij+u2i-1j; uu2=u2ij+1-2*u2ij+u2ij-1; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(m=0;mXN;m+) for(n=0;nZN;n+) u1mn=u2mn; u2mn=u3mn; if(k=100)for(m=0;mXN;m+)for(n=0;nZN;n+)u4mn=u3mn;/記錄波前快照 if(0=flag1 & 0=flag2)fp=fopen(wavefront14.dat,w);else if(1=flag1 & 0=flag2)fp=fopen(wavefront24.dat,w);else if(0=flag1 & 1=flag2)fp=fopen(wavefront12.dat,w);else if(1=flag1 & 1=flag2)fp=fopen(wavefront22.dat,w); if(fp!=NULL) for(i=0;iXN;i+)for(j=0;jZN;j+)fprintf(fp,%d %d %fn,i,j,u4ij);fclose(fp); void main()function(0,0); /時(shí)間二階差分格式、空間四階差分格式、震源1function(0,1); /時(shí)間二階差分格式、空間二階差分格式、震源1function(1,0); /時(shí)間二階差分格式、空間四階差分格式、震源2function(1,1); /時(shí)間二階差分格式、空間二階差分格式、震源2（1）時(shí)間二階差分格式、空間四階差分格式、震源1（2）時(shí)間二階差分格式、空間二階差分格式、震源1（3）時(shí)間二階差分格式、空間四階差分格式、震源2（4）時(shí)間二階差分格式、空間二階差分格式、震源2二是合成一個(gè)地震記錄，即記錄下與震源同一深度點(diǎn)的各點(diǎn)所有時(shí)刻的波場(chǎng)值，并指出記錄上的同向軸分別對(duì)應(yīng)哪些波。這時(shí)取小模型的常量，為減少頻散，速度v至少為2400程序圖像如下：#include #include #include #define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 200#define ZN 180#define DH 5#define DT 0.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNKN; /不能直接初值為0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2;if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0; u2ij=0.0; u3ij=0.0; u4ij=0.0; fij=0.0; if (j=30) vij=2400; /第一層速度為2400 else vij=3000; /第二層速度為3000 for(i=0;iXN;i+) for(k=0;kKN;k+) u5ik=0.0; /速度相同表示同一介質(zhì) f10010=1; /定義f函數(shù)，在100，10為1 for(k=0;kKN;k+) if(0=flag2) /時(shí)間二階差分格式、空間四階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; u5ik=u3i10; else if(1=flag2) /時(shí)間二階差分格式、空間二階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)