聲波方程數值模擬實驗報告學號：201107010230姓名：包靈指導老師：程冰潔老師完成時間：2014年11月26日星期三實驗要求：1、應用聲波方程作為正演模擬的波動方程；2、將所提供震源函數離散后繪圖；3、給定兩個二維速度-深度模型（一個小模型；一個大模型），繪出圖形來；4、對于小模型，整個區域的速度值可設為常數，即只有一種介質，將震源點放在模型中間，分別記錄兩個時刻的波前快照（即該時刻區域內所有網格點的波場值）。第一時刻為地震波還未傳播到邊界上的某時刻，第二時刻為地震波已經傳播到邊界上的某時刻，體會其人工邊界反射；5、對于大模型，定義為水平層狀速度模型（至少兩層）；做兩個實驗，一是將震源點放在區域表層任一點，記錄下某些時刻的波前快照，體會地震波在兩種介質的分界面上傳播規律；二是合成一個地震記錄，即記錄下與震源同一深度點的各點所有時刻的波場值，并指出記錄上的同向軸分別對應哪些波。實驗目的：1. 通過本次作業，加深對波動方程的理解，明白波動方程所代表的物理意義。2. 通過模擬地震波在介質中的傳播，理解實際勘探中地震波在地層中的傳播規律。3. 通過模擬水平層狀速度模型，體會地震波在兩種介質分界面的傳播規律，并能夠從地震記錄中識別出反射波，透射波，多次波，折射波和繞射波。4. 通過模擬人工合成的地震記錄，體會地震勘探基本原理和方法，驗證地震波傳播能量波形變化趨勢。需要的已知條件包括：1）震源函數2）地層速度（波速）3）邊界條件2彈性波方程：聲波方程的有限差分法數值模擬對于二維速度-深度模型，地下介質中地震波的傳播規律可以近似地用聲波方程描述： （4-1）是介質在點（x , z）處的縱波速度，為描述速度位或者壓力的波場，為震源函數。為求式（4-1）的數值解，必須將此式離散化，即用有限差分來逼近導數，用差商代替微商。為此，先把空間模型網格化（如圖4-1所示）。設x、z方向的網格間隔長度為，為時間采樣步長，則有： （i為正整數） (j為正整數) (n為正整數)表示在(i,j)點，k時刻的波場值。將在(i,j)點k時刻用Taylor展式展開： (4-2)將在(i,j)點k時刻用Taylor展式展開：(4-3)將上兩式相加，略去高階小量，整理得(i,j)點k時刻的二階時間微商為：(4-4)對于空間微分，采用四階精度差分格式，（以X方向為例）即將、分別在(i,j)點k時刻展開到四階小量，消除四階小量并解出二階微分得： (4-5)同理可得： (4-6)這就實現了用網格點波場值的差商代替了偏微分方程的微商，將上三個式子代入(4-1)式中得： (4-7)式中為介質速度的空間離散值，是空間離散步長，為時間離散步長，為震源函數，關于一般使用一個理論的雷克型子波代替，即：（1） （4-8）（2） （4-9）上式中，為時間， 為中心頻率，一般取為20-40HZ，為控制頻帶寬度的參數，一般取3-5。在實際計算過程中，需把此震源函數離散，參與波場計算。確定震源位置。穩定性條件：(4-10)這里表示的是地下介質的最大波速；若地下介質網格間隔、最大速度及時間采樣間隔不符合(4-8)式時，遞推求解（4-7）式，波場值會出現誤差（高階小量）累積，出現不穩定現象。頻散關系式： (4-11)式中為最小速度，為Nyquist（奈奎斯特）頻率。一般取震源子波中的主頻的2倍值參與計算，G為每個波長所占的網格點數，對于空間二階差分、時間二階差分取8，而對于空間為四階差分的情況則取4方能有效減少頻散。同理：對于空間微分，采用二階精度差分格式為 （4-12）考慮Reynolds邊界條件（詳細推到見文獻“基于Marmousi模型的聲波方程有限差分正演算法”），差分格式如下： （4-13）注：其中參數的設置：借用以前實驗的數據，當然可以根據限制條件4-10、4-11計算得到；至于震源放于101*5,101*5的矩形中心，根據速度與位移可以計算傳播到邊界時的時間，此處忽略。二實驗步驟1、 應用聲波方程作為正演模擬的波動方程，忽略轉換波的產生、傳播；2、 將所提供震源函數離散后繪圖；震源函數為雷克子波，離散繪圖如下：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、 對于小模型，整個區域的速度值可設為常數，即只有一種介質，將震源點放在模型中間，分別記錄兩個時刻的波前快照（即區域內所有網格點的波場值）。第一時刻為地震波還未傳播到邊界上的某時刻。由穩定條件，設v為2000，可以令DT=0.001,DH=5,此時精度較高，且滿足頻散關系程序，圖形如下：#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函數，當且僅當i,j同時為50時，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; /速度相同表示同一介質 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) /時間二階差分格式、空間四階差分格式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) /時間二階差分格式、空間二階差分格式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); /時間二階差分格式、空間四階差分格式、震源1function(0,1); /時間二階差分格式、空間二階差分格式、震源1function(1,0); /時間二階差分格式、空間四階差分格式、震源2function(1,1); /時間二階差分格式、空間二階差分格式、震源2（1）時間二階差分格式、空間四階差分格式、震源1（2）時間二階差分格式、空間二階差分格式、震源1（3）時間二階差分格式、空間四階差分格式、震源2（4）時間二階差分格式、空間二階差分格式、震源2第二時刻為地震波已經傳播到邊界上的某時刻，體會其人工邊界反射（此處用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函數，當且僅當i,j同時為50時，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; /速度相同表示同一介質 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) /時間二階差分格式、空間四階差分格式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) /時間二階差分格式、空間二階差分格式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; /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) /只需改動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); /時間二階差分格式、空間四階差分格式、震源1function(0,1); /時間二階差分格式、空間二階差分格式、震源1function(1,0); /時間二階差分格式、空間四階差分格式、震源2function(1,1); /時間二階差分格式、空間二階差分格式、震源2（1）時間二階差分格式、空間四階差分格式、震源1（2）時間二階差分格式、空間二階差分格式、震源1（3）時間二階差分格式、空間四階差分格式、震源2（4）時間二階差分格式、空間二階差分格式、震源24、 對于大模型，定義為水平層狀速度模型；做兩個實驗，一是將震源點放在區域表層任一點，記錄下某些時刻的波前快照，體會地震波在兩種介質的分界面上傳播規律，指出哪是反射波，哪是透射波；這時取小模型的常量，為減少頻散，速度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函數，在100，10為1 for(k=0;kKN;k+) if(0=flag2) /時間二階差分格式、空間四階差分格式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) /時間二階差分格式、空間二階差分格式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); /時間二階差分格式、空間四階差分格式、震源1function(0,1); /時間二階差分格式、空間二階差分格式、震源1function(1,0); /時間二階差分格式、空間四階差分格式、震源2function(1,1); /時間二階差分格式、空間二階差分格式、震源2（1）時間二階差分格式、空間四階差分格式、震源1（2）時間二階差分格式、空間二階差分格式、震源1（3）時間二階差分格式、空間四階差分格式、震源2（4）時間二階差分格式、空間二階差分格式、震源2二是合成一個地震記錄，即記錄下與震源同一深度點的各點所有時刻的波場值，并指出記錄上的同向軸分別對應哪些波。這時取小模型的常量，為減少頻散，速度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; /速度相同表示同一介質 f10010=1; /定義f函數，在100，10為1 for(k=0;kKN;k+) if(0=flag2) /時間二階差分格式、空間四階差分格式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) /時間二階差分格式、空間二階差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)