有限差分縱波波場模擬

1、彈性波課程報告 有限差分縱波波場模擬姓名： 鄧 健 學號： 1201410215 課程名稱：彈性波理論 任課老師：顧 漢 明 專業(yè)：地球探測與信息技術 時間： 2015年5月 評 語對課程論文的評語:平時成績：課程論文成績：總 成 績：評閱人簽名：注：1、無評閱人簽名成績無效；2、必須用鋼筆或圓珠筆批閱，用鉛筆閱卷無效；3、如有平時成績，必須在上面評分表中標出，并計算入總成績。摘要地震波場數(shù)值模擬是勘探地震學的重要研究課題之一,也是研究波動現(xiàn)象的重要手段。地震波場數(shù)值模擬常用的方法有偽譜法、有限差分法、有限元法等。有限差分法具有計算速度快、占用內存小等優(yōu)點，該方法對于近遠場及復雜邊界都有廣泛的

2、適用性，能夠準確地模擬波在各種介質及復雜結構地層中的傳播規(guī)律，是勘探地震學中應用最廣泛的數(shù)值計算方法。本文以波動理論為理論基礎，利用泰勒級數(shù)展開式推導出波動方程的有限差分格式及其離散表達式。針對傳統(tǒng)二階精度差分方法模擬精度較低的問題，推導出時間域二階空間域四階精度的有限差分方程，并綜合分析了初始條件、震源項、算法穩(wěn)定性及數(shù)值頻散等因素在有限差分法數(shù)值模擬中的影響。有限差分數(shù)值模擬計算是在一個有限的區(qū)域內進行的，當波傳播到邊界時就會產(chǎn)生邊界反射，這是進行數(shù)值模擬計算時所不期望出現(xiàn)的。本報告通過調研大量文獻資料，對效果較好的透明邊界條件進行了詳細分析。在此基礎上,分別推導得到這種邊界條件方法的離散

3、表達式，并在數(shù)值模擬過程中進行驗證。通過建立均勻模型、層狀模型及高速透鏡體模型，基于Matlab語言編程，論文實現(xiàn)了波動方程有限差分算法的數(shù)值模擬。對不同地質理論模型進行模擬，從模擬結果的對比分析中可以看出，模型參數(shù)的選擇及各參數(shù)間的相互關系對模擬結果有著顯著影響。無論是模擬精度，還是模擬計算效率，有限差分算法都具有一定優(yōu)勢。通過綜合研究，論文認為波動方程有限差分算法具有算法簡單、計算效率高、模擬精度較高等特點，理論研究意義大，應用前景廣闊。關鍵詞：數(shù)值模擬，波動方程，有限差分，邊界條件目錄1 前言11.1 地震波場數(shù)值模擬概述11.2 有限差分法及其與幾種常用數(shù)值模擬方法的比較12 波動方程

4、有限差分法數(shù)值模擬22.1 有限差分原理22. 2 波動方程的建立33 二維聲波方程有限差分格式的建立63.1 聲波方程63.2 波動方程有限差分格式的建立73.3 高階差分方程94 波動方程有限差分的幾個問題124.1 初始條件124.2 震源函數(shù)124.3 差分方程的穩(wěn)定性及收斂性134.4 頻散問題135 有限差分算法中邊界條件的處理145.1 一維透明邊界條件155.2 二維透明邊界條件165.3 透明邊界條件的差分格式186 簡單模型試算206.1均勻模型206.2 水平層狀模型和高速透鏡體模型227結論與建議25參考文獻26附錄281 前言1.1 地震波場數(shù)值模擬概述地震波場數(shù)值模

5、擬簡單說來，就是已知地下介質構造及其參數(shù)，再用理論計算方法來研究地震波在地下介質中的傳播規(guī)律，合成地震記錄的一種技術方法。隨著地震勘探技術的發(fā)展，數(shù)值模擬方法已經(jīng)貫穿于地震數(shù)據(jù)的采集、處理和解釋的全過程，而且在確定觀測的合理性、檢驗處理和解釋的正確性等方面都有了廣泛的應用1-3。地震波數(shù)值模擬問題的研究主要包括三個方面：(1)地震波場數(shù)值模擬原理；(2)地震波場數(shù)值模擬算法；(3)數(shù)值模擬的計算機實現(xiàn)。地震波場數(shù)值模擬方法的理論基礎主要是波動理論和射線理論。射線理論方法是求出質點運動方程的漸進近似解，這種方法計算速度快，但是精度較低。而波動理論方法則是用數(shù)值計算方法直接求出波動方程的解，其模擬

6、結果中包含豐富的波動信息，模擬結果較為精確。數(shù)值模擬又可分為波動方程的解析解和數(shù)值解。對簡單的地質模型就可以得到其精確解，但對于復雜的地質模型，只有通過數(shù)值解來說明波在底下介質中的傳播。通常解析解只是作為數(shù)值解的一種檢驗4-6。1.2 有限差分法及其與幾種常用數(shù)值模擬方法的比較在地震勘探中，為了研究地震波在地下各種介質中的傳播規(guī)律，就需要對波場進行數(shù)值模擬。常用的數(shù)值模擬方法有偽譜法、有限差分法，有限元法等。有限元法依據(jù)變分原理，通過靈活的網(wǎng)格剖分，用簡單形態(tài)逼近實際的地質體，能處理多種介質和自然邊界條件。適用于非規(guī)則網(wǎng)格的計算，方便有效，模擬的精度高。但是有限元法的問題是不適用于大規(guī)模的模型

7、的計算，而且計算量大7。偽譜法在20世紀七十年代引入數(shù)值模擬計算領域的，它是利用傅立葉變換來計算波場的空間導數(shù)，用差分方法來計算時間的導數(shù)，可以看成是一種無限階的有限差分法，是傳統(tǒng)有限差分法的一個推廣。偽譜法在粗網(wǎng)格上也能實現(xiàn)高精度的計算，相對有限元法實現(xiàn)起來較容易，在非線性波動問題及氣象預測等領域有著廣泛的應用8-10。有限差分法是一種最常用的數(shù)值模擬方法，它是將波動方程中波場函數(shù)的空間導數(shù)和時間導數(shù)用相應的空間、時間的差分代替。有限差分法具有計算速度快、占用內存小等優(yōu)點，該方法對于近遠場及復雜邊界都有廣泛的適用性，能夠準確地模擬波在各種介質及復雜結構地層中的傳播規(guī)律。有限差分法方法簡單、高

8、效的優(yōu)點是其他方法難以比擬的，因此有限差分法目前仍然是勘探地震學中應用最廣泛的數(shù)值計算方法11。2 波動方程有限差分法數(shù)值模擬2.1 有限差分原理有限差分法就是把求解區(qū)域劃分為差分網(wǎng)格，然后用有限的網(wǎng)格節(jié)點代替連續(xù)的求解域，利用微商與差商的近似關系將描述介質傳播的微分方程轉化為差分方程進行求解。差分離散的方法有兩種，一種是將單變量的二階波動方程直接轉化為時間空間的二階中心差分進行離散求解；第二種方法是把用位移表示的二階波動方程轉化為用應力及質點速度表示的一階方程組，然后用應力和速度的交錯網(wǎng)格求解,8。構造差分的方法也有很多種，一般采用泰勒級數(shù)展開法。常見的差分格式有顯示差分、隱式差分和顯隱交替

9、格式；按照差分的精度又可以分為一階差分、二階差分和高階差分等，目前通常見到的差分格式主要是幾種差分格式的組合。我們知道，差分方程的建立首先選擇網(wǎng)格布局和差分形式，然后以有限差分代替無限微分，以差商代替微商，并以差分方程代替微分方程及其邊界條件，最后建立差分方程12。在建立差分方程的時候要注意到兩個方面。一是合理選擇網(wǎng)格布局與步長。我們將離散后各相鄰離散點之間的距離或者離散化單元的長度稱為步長。在所選定的區(qū)域內進行網(wǎng)格劃分是差分方程建立的第一步，其方法比較靈活，但是實際應用中往往遵循誤差最小原則。常用的典型的網(wǎng)格剖分方法有矩形網(wǎng)格、三角形網(wǎng)格等，網(wǎng)格樣式的選擇和區(qū)域的形狀有關。其次是將微分方程轉

10、化為差分方程。這個過程就是選擇一種差分格式代替微分形式的過程。構造差分的方法有多種形式，本文采用的是泰勒級數(shù)展開法2,5,13。地震波場數(shù)值模擬以地震波動理論為基礎，用有限差分法解波動方程時，對變量離散化，也即對連續(xù)的物理量只考慮其在離散的空間位置和離散時刻的值，然后把方程中的導數(shù)用這些離散的采樣值表示14。對于一個單變量的函數(shù)f(x)，將其離散化，那么在采樣點x =lx的采樣值就是f(lx),其中x表示步長，l為整數(shù)。則有限差分法中f(x)在采樣點x=lx的導數(shù)就可以近似表示為： 其中an是系數(shù)，N表示差分格式的長度，差分的格式是由這兩個系數(shù)來決定。在實際應用中常用的差分格式有向前差分、向后

11、差分以及中心差分：（1）向前差分：（2）向后差分：（3）中心差分：2. 2 波動方程的建立建立波動方程是在已知物體形狀、位置、彈性常數(shù)及外力分布等參數(shù)的情況下求出物體的位移、應力與應變的分布。這個問題具體包括以下內容：（1）應力部分的平衡方程（2）應變部分（3）應力與應變的關系將式（2-7）代入（2-6），推導可以得下式：方程（2-8）是物體在平衡狀態(tài)下的平衡方程。當物體處于不平衡狀態(tài)時，式（2-8）則變?yōu)椋涸诙S情況下，我們只需要考慮x，z方向的位移分量，這時式（2-9）就可轉換為關于x,z的方程：當fx=fz=0時，上面的這兩個方程就變?yōu)椋海?-11）就是二維非均勻介質彈性波方程。而在均勻

12、介質中，、和均為常數(shù)，則二維均勻介質彈性波方程就為：把（2-12）中的x和z分量合并就得到了二維均勻介質彈性波方程的矢量形式表達式在沒有彈性橫波只有彈性縱波存在時，對（2-13）兩邊取散度：利用旋度與散度的關系，交換上式的微分次序并化簡可得：其中表示縱波速度。根據(jù)赫姆霍茲在他著名的有關渦流運動的著作中證明了下屬定理：任何向量點函數(shù)，若它的散度和旋度具有位，則它可以表示為一個無旋部分和一個旋轉部分之和。亦即對于任何一種向量場，如果在其定義域內有散度和旋度，則該向量場可表示為一個標量位的梯度場和一個向量位的旋度場之和，而空間傳播的波動正是無旋運動和旋轉運動這兩種運動之和15-20 令代入式（2-8

13、）中，取其中的無旋部分則有：然后再消去式中的梯度,就得到了用位移位表示的縱波方程：3 二維聲波方程有限差分格式的建立3.1 聲波方程一般地，二維非均勻介質的聲波波動方程可表示為：式中 U = U ( x , z , t)表示聲壓，V表示聲波在介質中的傳播速度，是密度，它是隨空間各點而變的。其中 s ( x , z , t )為震源函數(shù)。對于均勻介質，密度為常數(shù)，則二維聲波波動方程就可以表示為：在實際問題中，我們總是用一個有限寬頻帶的時間函數(shù)來代替 s ( x , z , t )函數(shù)，以便能夠真實地反映地震波的傳播21。3.2 波動方程有限差分格式的建立上一節(jié)中討論了差分的原理和幾種常見的差分格

14、式，在這個基礎上，我們結合波動方程，對縱波波動方程的有限差分格式進行推導。令，其中x，z是空間間隔，t是時間間隔。用k表示時間方向的離散網(wǎng)格，m表示x方向的離散網(wǎng)格，n表示z方向的離散網(wǎng)格。利用泰勒級數(shù)展開式將 在 展開可得：同理，在處的展開式為：用（3-3）式減去（3-4）式，就可以推出關于t的一階中心差分：如果將（3-3）式與（3-4）式相加可得進一步推導就可以得到關于t的二階中心差分：同理也可以推導出關于x，z的中心差分格式：關于x的一階中心差分關于x的二階中心差分： 關于z的一階中心差分：關于z的二階中心差分：若令x =z=h，利用上面關于x,z,t的中心差分方程就可以得到二維波動方程

15、的有限差分方程：對上式繼續(xù)推導可得：其中，上式就是二階波動方程的有限差分格式。如果我們繼續(xù)利用泰勒級數(shù)展開式，則可以得到：這樣我們就能夠得到時間域二階空間域四階的波動方程有限差分格式：上式中k表示時間方向的離散網(wǎng)格，m表示x方向的離散網(wǎng)格，n表示z方向的離散網(wǎng)格，。需要注意的是，這里同樣也假設x=z=h，即網(wǎng)格步長相等5 ,22-24。3.3 高階差分方程將采用高階差分會有很多時間層，計算較復雜。為此時間上采用2階，空間上可采用高階。為了方便起見，本文空間步長在x，z軸上相等，則通過上述類似推倒，可以得到不同精度的聲波差分方程。時間2階，空間6階時間2階，空間8階時間2階，空間10階時間2階，

16、空間12階時間2階，空間14階圖3-1為均勻介質中各階精度波場快照對比圖。模型參數(shù)為網(wǎng)格500×500，震源位置（250，250），空間步長x和z方向均為1m,采樣時間間隔100us。所用的震源是雷克子波，子波頻率30Hz。從下圖可以看出，時間2階，空間2階精度波場出現(xiàn)明顯頻散，隨著空間精度的增高，頻散逐漸減弱。所以提高差分精度可以有效減輕頻散，但計算量也隨之增加。圖3-1均勻介質中各階精度波場快照圖4 波動方程有限差分的幾個問題4.1 初始條件在進行波動方程有限差分數(shù)值模擬時，初始時刻的波場值是不知道的，而在計算差分方程時，只有知道時間間隔前一時刻的值才能遞推計算出后面時刻的值。因

17、此在計算開始時就必須要考慮到初始條件。一般假設在震源發(fā)生震動前，各個質點均處于靜止狀態(tài)：初始時刻波場的速度也同樣為零，即4.2 震源函數(shù)震源函數(shù)的給出通常有兩種方式：一種是初值法，就是將理論結果當作初始值給出；另外一種是力源法，就是以力源的方式給出。初值法的震源除了自由表面或內界面附近，可以在模型的其他任意處進行定義；而力源法可以將震源定義在自由表面的附近，但必須是在網(wǎng)格點上。兩種方法各有其優(yōu)缺點1,26。在進行波動方程數(shù)值模擬計算的過程中，震源問題的處理對模擬的結果有著重要的影響。子波是描述震源的時間延續(xù)特征的時間函數(shù)，對于地震子波而言，子波延續(xù)的時間越短，頻帶越寬，子波的垂直分辨率也就越高

18、。由于做有限差分計算時會出現(xiàn)數(shù)值頻散，尤其當空間采樣不足時，子波的高頻成分頻散就會更嚴重。所以要根據(jù)模型的速度及網(wǎng)格間距合理選擇子波主頻。本文采用了雷克子波和高斯子波作為震源函數(shù)27。4.3 差分方程的穩(wěn)定性及收斂性一個有實際應用意義的數(shù)值模擬算法必須是穩(wěn)定的，也就是說對算法的數(shù)值解與解析解的差值是有限定的。對于有限差分法，一般來說如果差分方程的解的誤差不隨時間的推進而增加，那么就說該差分方程的解是穩(wěn)定的，這也就是說差分方程的解是有界的。差分的穩(wěn)定條件根據(jù)差分階數(shù)的不同而有所不同，差分的階數(shù)越高，穩(wěn)定性越差。對于不同的差分格式，參數(shù)的選取一般也是不同的5,28,29。對于均勻介質，差分方程的解

19、的穩(wěn)定性是確定的，穩(wěn)定性條件如下式所示：其中V是均勻介質中波的傳播速度，應當注意，這里的x是差分格點中x與z的最小值。而對于非均勻介質應滿足：其中Vmax是各種均勻介質中的最大波速值。4.4 頻散問題在進行地震波數(shù)值模擬時，我們希望盡可能同時達到較高的模擬精度和較快的計算速度。但是在用微分方程數(shù)值模擬的方法求解波動方程時，就會產(chǎn)生數(shù)值頻散或網(wǎng)格頻散。這種頻散是用微分方程求解波動方程時固有的本質特征，是避免不了的，只能盡量減弱這種頻散2。通過對有限差分法數(shù)值頻散的理論分析研究可以知道，影響頻散的因素有很多。在有限差分的計算過程中，不僅需要對空間進行網(wǎng)格化離散，同時還需要對時間進行離散。波動方程經(jīng)

20、過空間離散化之后就引入了差分算子誤差；而時間間隔的選擇也對數(shù)值模擬的精度和數(shù)值頻散有一定影響。當子波頻率越高，則一個波場內的網(wǎng)格點數(shù)就越少，則頻散現(xiàn)象越嚴重；對于一定頻率的子波，介質的速度越低，則一個波場內的網(wǎng)格點數(shù)也少，數(shù)值頻散現(xiàn)象也越嚴重。因此，從理論的角度來說，在進行數(shù)值模擬時，對于給定的子波頻率，提高時間和空間的差分階數(shù)，減小網(wǎng)格步長和時間間隔都可以提高數(shù)值計算的精度和改善數(shù)值頻散情況。但是過細的網(wǎng)格剖分會增大計算所需的存儲量，從而增加計算的時間，所以就需要選取合適的參數(shù)，在保證精度和計算速度的同時盡量減少頻散29-33。5 有限差分算法中邊界條件的處理我們在計算機上進行有限差分數(shù)值模

21、擬計算時，由于受到計算機內存和計算速度等因素的制約，只能在一個有限的區(qū)域內進行計算。這樣就會在兩側與底界上產(chǎn)生邊界，這個邊界就是人工邊界。當波傳播到人工邊界時就會產(chǎn)生反射，反射必然會造成一定的干擾，產(chǎn)生邊界效應。從圖5-1（左）中可以看出，在沒有加入邊界條件時，當波傳播到邊界時，在邊界處產(chǎn)生了很強的反射，這是我們在進行數(shù)值模擬計算時所不期望出現(xiàn)的。為了消除或減弱這種邊界反射效應，得到地質地層真實的反射信息，就需要對人工邊界進行處理，從而得到更接近于實際空間中波的傳播規(guī)律。在處理邊界反射時常用的邊界條件方法有高衰減區(qū)法、傍軸近似法、吸收邊界法和反周期擴展法1,2,25。高衰減區(qū)法是在靠近邊界處引

22、入一個高吸收區(qū)，通過耗散因子使入射波仔這個吸收區(qū)域內逐漸衰減，從而抑制邊界附近的人工反射。旁軸近似法是將邊界區(qū)的波動方程用單向的傍軸近似波動方程代替，只模擬波能量由差分網(wǎng)絡向外邊界的單向傳播，從而消除邊界反射的問題，不過這個方法只適用于小角度入射的情況。吸收邊界法是先導出吸收邊界條件方程，然后讓方程與計算區(qū)域內的波動方程聯(lián)立求解，從而使從人工邊界向計算區(qū)域反射回來的波能夠全部或部分被吸收。反周期擴展法，即利用正反周期函數(shù)極性相反的特點消除回繞波場，這個方法在理論上是能夠完美的消除回繞波場而不產(chǎn)生任何的反射，但是在實際應用時，由于這種方法需要進行反周期波場的擴展，大大的增加了計算量，而且為了得到

23、最后疊加的平均波場，正反周期波場都必須存儲，這樣同時也會占用很大內存空間25,34,35。下面將就本文將要用到的透明邊界條件進行討論。5.1 一維透明邊界條件對于一維均勻介質波動方程：式中 V 表示聲波的傳播速度，U 表示聲壓。我們引入一個向右行波（k =/V），將這個平面波的方程代入上面的均勻介質波動方程（5-1）中可以解得：其中r表示反射系數(shù)。在右邊界x=a上反射波與入射波的振幅相等，即|r|=1。當引入一個左行波時，左邊界x=a上同樣也會產(chǎn)生反射，且|r|=1。由此可以看出，如果邊界條件選取不當會使左右兩側的邊界上產(chǎn)生強反射，因此就需要在左右邊界處選取的邊界條件能使r=0。則右行波就滿足

24、方程左行波就滿足方程由此可以看出，如果在右行波上加上一個左行波項，則只有當r = 0時方程（5-5）才成立，即當選式（5-5）作為右邊界條件時，則在右邊界上不會產(chǎn)生反射，因此可以將式上式作為右邊界條件。同理可得左邊界條件為通常用右行波和左行波之和來表示波函數(shù)，因此也只有當 r = 0時，波函數(shù)才滿足式（5-4），所以在引入邊界條件以后左邊界和右邊界上都不會產(chǎn)生反射。5.2 二維透明邊界條件設在平面區(qū)域為，對于已知某初始條件下的二維波動方程如果給定一個邊界條件則平面波在邊界x = a、 x = a及 z = b處會產(chǎn)生反射?？紤]一個向右傳播的平面波其中表示入射角，即平面波前與 x 軸的夾角。那么

25、這時反射波為其中r為邊界產(chǎn)生的反射波的反射系數(shù)。則總波場可以表示為以x=a處為例，把式（5-11）代入U(a,z,t)=0可得|r|=1，即在右邊界上產(chǎn)了反射。由于對于波動方程（5-6）當波只在水平方向傳播時，即=0時在右邊界x = a，如果則由式（5-12）可得把式（5-11）代入式（5-15）可得r=0，這時在右邊界x=a上就不會產(chǎn)生反射。同理可得在左邊界x=a處，要使r=0，則取如果取并在邊界上令當x=a時取當x=a時取以上這兩種邊界條件都只有很小的邊界反射，但還需要考慮到有限差分法的穩(wěn)定性條件。如果把項的系數(shù)看成 當?shù)奶厥馇闆r下，的系數(shù)就為1/2。則（5-19）式可以寫為由波動方程式（

26、5-6）和上式聯(lián)立推導可得對式（5-22）進行整理化簡就可得，x = a時的右邊界條件為同理可得x=-a時的左邊界條件為當z=b時，底邊界條件為以上得到的式（5-23）、式（5-24）及式（5-25）就分別為二維均勻介質波動方程的右邊界條件、左邊界條件和底邊界條件。5.3 透明邊界條件的差分格式在推導出二維均勻介質波動方程的透明邊界條件基礎上，下面在進一步推導一下這些邊界條件的差分格式。對右邊界條件（5-24）兩邊同乘以h可得下式其中h=x=z，s=Vt/x。由于通過泰勒級數(shù)展開式有因此把式（5-28）和式（5-27）代入式（5-26）可得上式中的微分項可以表示為然后把式（5-30）各式代入式

27、（5-29）即得到右邊界條件的差分格式同理可以得到左邊界條件x=a的差分格式為底邊界條件z = b的差分格式為頂邊界條件z=0的差分格式為上面得到的式（5-34）、式（5-33）、式（5-32）及式（5-31）就是二維均勻介質波動方程透明邊界條件的差分格式。如下圖5-1（右）所示，可以看出在加入了邊界條件后，邊界反射明顯減少，說明透明邊界條件對邊界反射有很好的吸收效果。圖5-1 未加邊界條件（左）與加入邊界條件（右）波場快照對比圖（所加邊界條件為透明邊界，差分精度為時間2階空間2階）6 簡單模型試算6.1均勻模型模型速度為2000m/s，密度為常數(shù)，網(wǎng)格大小為500×500，網(wǎng)格空間

28、采樣步長為1m，采樣時間間隔為100s，震源坐標為（250m，20m），震源震動時間5ms，子波為雷克子波，頻率為30Hz，模型見圖6-1。圖6-2中包含了20-120ms的波場快照，差分精度為時間2階空間12階，邊界條件為透明邊界條件。從圖中可以看出，波以震源為中心，呈圓弧向外傳播，當波傳播至頂邊界時有微弱的反射波，說明利用透明邊界條件并不能無完全消除人工邊界產(chǎn)生的反射波。另外，從圖中還可以看出有輕微的頻散現(xiàn)象。 圖6-1均勻模型示意圖圖6-2 均勻模型波場快照圖（差分精度為時間2階，空間12階，邊界條件為透明邊界條件）6.2 水平層狀模型和高速透鏡體模型圖6-3（左）所示的水平層狀模型，模

29、擬網(wǎng)格大小為500×500，網(wǎng)格間距1m；表層速度為2000m/s，厚度為100m，第二層速度為3000m/s，厚度為200m，第三層速度為3500m/s，厚度為200m。采樣時間間隔為100s，震源坐標為（250m，20m），持續(xù)時間為5ms，子波選用雷克子波，子波頻率為30Hz。圖6-3（右）所示的高速透鏡體模型，模擬網(wǎng)格大小為500×500，網(wǎng)格間距1m。透鏡體為一個半長軸為100m，半短軸為50m的橢圓，中心坐標為：z=250m，x=250m；高速透鏡體的速度為3500m/s，周圍介質速度為2000m/s。采樣時間間隔為100s，震源坐標為（250m，20m），持續(xù)

30、時間為5ms，子波選用雷克子波，子波頻率為30Hz。圖6-4中T=40ms時，波在第一層與第二層的界面上開始傳播，從圖中T=60,80ms中可以清晰的看到第一個界面的反射波和透射波，T=100ms時波還沒有傳播到第二層與第三層的界面，T=120ms時波傳播第三層，從圖中可以看到產(chǎn)生了一個反射波。頂邊界反射波比較明顯，同時，隨著波向外傳播，頻散現(xiàn)象也逐漸明顯。圖6-5中T=80ms時，波時波傳播到透鏡體內，在透鏡體的頂界面產(chǎn)生了反射波，T=140ms波場快照中，出現(xiàn)了波從透鏡體內向外傳播時產(chǎn)生的反射波。T=100ms的波場快照中頻散現(xiàn)象較為明顯。圖6-3 層狀模型（左）和高速透鏡體模型（右）示意

31、圖圖6-4 水平層狀模型波場快照圖（差分精度為時間2階，空間12階，邊界條件為透明邊界條件）圖6-5 高速透鏡體模型波場快照圖（差分精度為時間2階，空間12階，邊界條件為透明邊界條件）7結論與建議地震波場的數(shù)值模擬是地震勘探學的重要研究內容之一，是進一步研究地震資料的采集、處理和解釋的有效輔助手段。論文以波動理論為基礎，對有限差分模擬的一些關鍵性問題進行了探討和研究，主要包括波動方程有限差分格式的建立、邊界條件、初始條件、震源問題、有限差分的穩(wěn)定性及頻散問題等。經(jīng)過本次課程報告，我對有限差分波場模擬有以下認識：（1）有限差分算法做波動方程數(shù)值模擬時，具有計算速度快、算法簡單等特征。（2）有限差

32、分波動方程是利用泰勒級數(shù)展開式展開推導得到的，差分的階數(shù)越高，有限差分的誤差就越小，而有限差分的解就越精確。（3）有限差分模擬時，網(wǎng)格空間步長的選擇不僅對模擬的精確度有影響，而且對邊界的處理效果也有著明顯的影響。同時由于震源項是以離散形式加入的，所以網(wǎng)格步長也必須滿足采樣定理，這樣才能取到完整的子波，以保證計算結果不會出現(xiàn)較大的誤差。因此，有限差分數(shù)值模擬計算時，必須合理地選擇網(wǎng)格步長。（4）數(shù)值模擬過程是在特定的區(qū)域內進行，因此就必然會產(chǎn)生人工邊界，波傳播到邊界就產(chǎn)生了反射，通過對幾種模型的模擬表明，本文所用的透明邊界條件和吸收邊界條件對邊界反射均有較好的吸收效果。本課程報告的工作在主要是在

33、前人的研究基礎上開展的，但是還存在一些問題需要進一步的完善和改進，建議主要從以下幾個方面進行：（1）網(wǎng)格剖分是影響模擬結果的主要因素之一，本文采用的仍然是常規(guī)的規(guī)則網(wǎng)格剖分，規(guī)則網(wǎng)格剖分方法在模擬計算中仍然存在一些缺陷，當模型較復雜或要求的精度較高時，往往不能達到理想的模擬效果。如果減小網(wǎng)格間距，增加網(wǎng)格數(shù)量，又會導致計算速度和計算效率降低，因此為了獲得較滿意的模擬結果，需要采用更優(yōu)的網(wǎng)格剖分方法，以克服存在的缺陷。（2）對于模擬中出現(xiàn)的數(shù)值頻散現(xiàn)象，還需要更深入的研究，以便能夠盡可能的減小數(shù)值頻散。（3）本報告的算法主要針對簡單地質模型，如果對于復雜的介質模型，算法還需要進一步的研究和完善。

34、參考文獻1裴正林，牟永光.地震波傳播數(shù)值模擬.地球物理學進展J，2004，19(4)：933941.2馬在田.計算地球物理學概論M.上海：同濟大學出版社，1997.3何樵登.地震波理論M .地質出版社，19884何樵登，熊維鋼.應用地球物理教程地震勘探M .地質出版社，19915吳清嶺，張平，施澤龍.波動方程正演模擬即應用J.大慶石油地質與開發(fā)，1998，17（3）6江玉樂，雷苑著.地球物理數(shù)據(jù)處理教程M.北京：地質大學出版社，2006，77姚姚著.地震波場與地震勘探.北京：地質出版社，2006，6.8張永剛.地震波數(shù)值模擬方法J.石油物探，2003，42（2）：143-147.9陳偉.起伏地

35、表條件下二維地震波場的數(shù)值模擬J.勘探地球物理進展，2005，28（1）10吳永剛，吳清嶺.有限差分法彈性波場數(shù)值模擬J.大慶石油地質與開發(fā)，1994，13（3）11Alterman. Propagation of elastic wave in layered media by finite difference methods. Bulletin of the Seismological Society of America,1968;58(1):36739812Alford, R. M. Kelly ,K .R. and Booer, D. M. Accuracy of finite d

36、ifference modeling of the acoustic wave equation.Geophysics,1974;39(6):83484213周家紀，賀振華.模擬地震波傳播的大網(wǎng)格快速差分算法J.地球物理學報，1994；37（增刊）：45045414董良國.一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究J.地球物理學報，2000，43（6）：85615邵秀民，藍志凌非均勻各向同性介質中地震波傳播的數(shù)值模擬J.地球物理學報，1935；38(1)：395516劉洋，李承楚，牟永光.任意偶數(shù)階精度有限差分數(shù)值模擬J.石油地球物理勘探，1998，53（1）17褚春雷，王修田.非規(guī)則三角網(wǎng)

37、格有限差分法地震正演模擬J.中國海洋大學學報，2005，35（1）：04304818董良國，李培明地震波傳播數(shù)值模擬中的頻散問題J.天然氣工業(yè)，2004，24(6)：53-5619奚先，姚姚.二維隨機介質及波動方程正演模擬J.石油地球物理勘探，2001，36（5）：54655220董良國等.一階彈性波方程交錯網(wǎng)格高階差分解法J.地球物理學報，2000，43（3）:41141921奚先，姚姚.二維彈性隨機介質中的波場特征J.石油地球物理勘探，2004，39（6）：67968522邵秀民，藍志凌非均勻各向同性介質中地震波傳播的數(shù)值模擬J.地球物理學報，1935；38(1)：395523孫若昧.地震

38、波傳播有限差分模擬的人工邊界問題J.地球物理學進展，1996，11（3）：5324 Cerjan C, Kosloff D, Kosloff R and ReshefM. A nonreflecting boundary condit -ion for discrete acoustic and elastic wave equation. Geophysics, 1985, 50 (4) : 70570825Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equation.B

39、ulletin of the Seismological Society of America, 1977;66(11):1529154026Reynolds. Boundary conditions for the numerical solution of wave propagation problem.Geophysics, 1978;43(6):1099111027崔興福，張關泉.地震波方程人工邊界的兩種處理方法J.石油物探，2003，42（4）：45245528張毅.聲波正演模擬中的人工邊界問題J.工程地球物理學報，2007，4（4）29李文杰，魏修成，劉洋.聲波正演中一種新的邊界

40、條件-雙重吸收邊界條件J.石油物探，2004，43（6）：52853130邢麗.地震波數(shù)值模擬中的吸收邊界條件J.上海第二工業(yè)大學學報，2006，23（4）31Dahlain M A. The application of high difference to the scalar waveequation. Geophysics, 1986;51(1): 546632 Higdon R L. Absorbing boundary condition for elastic waves. Geophysics,1991, 56 (2) :23124133孫成禹，張吉輝完全縱波方程有限差分數(shù)值模

41、擬J.石油地球物理勘探，2005；40(3)：28929434黃自萍，徐伶俐，周繼順.起伏地表聲波方程的數(shù)值模擬J.石油地球物理勘探，2006，41（3）：27528035宛書金，董敏煌等.橫向各向同性介質中的地震波傳播特性J.石油大學學報，2002,26(1)：2328附錄本課程報告所編matlab程序% 2015.5.20 縱波波場模擬%clcclear% 初始化mode=12; %選擇差分介數(shù)，最高支持16介Model=load('Model.mat');V_model=Model.UniformModel;%速度模型c=load('c.mat');c=c

42、.c; %系數(shù)矩陣cellnum=size(V_model); %網(wǎng)格數(shù)目dt=100*10(-6); %時間采樣間隔 100usdx=1; dz=1; %空間采樣間隔x0=250; z0=20; %震源坐標f0=30; %震源子波頻率S_t=2*10(-3); %震源持續(xù)時間 2ms T=140*10(-3); %波場快照時間 %U0=zeros(cellnum(1),cellnum(2); %波場初始值U1=zeros(cellnum(1),cellnum(2);U2=zeros(cellnum(1),cellnum(2);for k=1:T/dt; if k<S_t/dt U1(z

43、0+1,x0+1)=sin(2*pi*f0*k*dt)/(k*dt);%震源函數(shù)。 end % if mode=2 for m=mode/2+1:cellnum(2)-mode/2 % 2介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end % if mode=4 for m=mode/2+1:cellnum(2)-

44、mode/2 % 4介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end % % if mode=6 for m=mode/2+1:cellnum(2)-mode/2 % 6介差分方

45、程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. c(mode/2,3)*a2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end

46、% if mode=8 for m=mode/2+1:cellnum(2)-mode/2 %8介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. c(mode/2,3)*a2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1

47、(m,n+3)-4*U1(m,n)+. c(mode/2,4)*a2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end % if mode=10 for m=mode/2+1:cellnum(2)-mode/2 % 10介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. c(mode/2,3)*a2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n)+. c(mode/2,4)*a2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n)+. c(mode/2,5)*a2*(U1(m+5,n)+U1(m-5,n)+U1(