有限差分縱波波場(chǎng)模擬

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

35、表?xiàng)l件下二維地震波場(chǎng)的數(shù)值模擬J.勘探地球物理進(jìn)展，2005，28（1）10吳永剛，吳清嶺.有限差分法彈性波場(chǎng)數(shù)值模擬J.大慶石油地質(zhì)與開發(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周家紀(jì)，賀振華.模擬地震波傳播的大網(wǎng)格快速差分算法J.地球物理學(xué)報(bào)，1994；37（增刊）：45045414董良國.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法穩(wěn)定性研究J.地球物理學(xué)報(bào)，2000，43（6）：85615邵秀民，藍(lán)志凌非均勻各向同性介質(zhì)中地震波傳播的數(shù)值模擬J.地球物理學(xué)報(bào)，1935；38(1)：395516劉洋，李承楚，牟永光.任意偶數(shù)階精度有限差分?jǐn)?shù)值模擬J.石油地球物理勘探，1998，53（1）17褚春雷，王修田.非規(guī)則三角網(wǎng)

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

38、波傳播有限差分模擬的人工邊界問題J.地球物理學(xué)進(jìn)展，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崔興福，張關(guān)泉.地震波方程人工邊界的兩種處理方法J.石油物探，2003，42（4）：45245528張毅.聲波正演模擬中的人工邊界問題J.工程地球物理學(xué)報(bào)，2007，4（4）29李文杰，魏修成，劉洋.聲波正演中一種新的邊界

40、條件-雙重吸收邊界條件J.石油物探，2004，43（6）：52853130邢麗.地震波數(shù)值模擬中的吸收邊界條件J.上海第二工業(yè)大學(xué)學(xué)報(bào)，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孫成禹，張吉輝完全縱波方程有限差分?jǐn)?shù)值模

41、擬J.石油地球物理勘探，2005；40(3)：28929434黃自萍，徐伶俐，周繼順.起伏地表聲波方程的數(shù)值模擬J.石油地球物理勘探，2006，41（3）：27528035宛書金，董敏煌等.橫向各向同性介質(zhì)中的地震波傳播特性J.石油大學(xué)學(xué)報(bào)，2002,26(1)：2328附錄本課程報(bào)告所編matlab程序% 2015.5.20 縱波波場(chǎng)模擬%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); %時(shí)間采樣間隔 100usdx=1; dz=1; %空間采樣間隔x0=250; z0=20; %震源坐標(biāo)f0=30; %震源子波頻率S_t=2*10(-3); %震源持續(xù)時(shí)間 2ms T=140*10(-3); %波場(chǎng)快照時(shí)間 %U0=zeros(cellnum(1),cellnum(2); %波場(chǎng)初始值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(