有限差分縱波波場模擬

,.——有限差分縱波波場模擬波理論與信息技術(shù)任課老師：顧漢明評(píng)語對(duì)課程論文的評(píng)語:對(duì)課程論文的評(píng)語:課程課程論文成績：評(píng)閱人簽名：平時(shí)成績：總成績：注：1、無評(píng)閱人簽名成績無效；2、必須用鋼筆或圓珠筆批閱，用鉛筆閱卷無效；3、如有平時(shí)成績，必須在上面評(píng)分表中標(biāo)出，并計(jì)算入總成績。摘要地震波場數(shù)值模擬是勘探地震學(xué)的重要研究課題之一,也是研究波動(dòng)現(xiàn)象的重要手段。地震波存小等優(yōu)點(diǎn)，該方法對(duì)于近遠(yuǎn)場及復(fù)雜邊界都有廣泛的適用性，能夠準(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í)間域二階空間域四階精度的模擬中的影響。有限差分?jǐn)?shù)值模擬計(jì)算是在一個(gè)有限的區(qū)域內(nèi)進(jìn)行的，當(dāng)波傳播到邊界時(shí)就會(huì)產(chǎn)好的透明邊界條件進(jìn)行了詳細(xì)分析。在此基礎(chǔ)上,分別推導(dǎo)得到這種邊界條件方法的離散表達(dá)式，并在數(shù)值模擬過程中進(jìn)行驗(yàn)證。模型參數(shù)的選擇及各參數(shù)間的相互關(guān)系對(duì)模擬結(jié)果有著顯著影響。無論是模擬精度，還是模擬計(jì)簡單、計(jì)算效率高、模擬精度較高等特點(diǎn)，理論研究意義大，應(yīng)用前景廣闊。關(guān)鍵詞：數(shù)值模擬，波動(dòng)方程，有限差分，邊界條件1前言 11.1地震波場數(shù)值模擬概述 11.2有限差分法及其與幾種常用數(shù)值模擬方法的比較 12波動(dòng)方程有限差分法數(shù)值模擬 22.1有限差分原理 22.2波動(dòng)方程的建立 43二維聲波方程有限差分格式的建立 73.1聲波方程 73.2波動(dòng)方程有限差分格式的建立 83.3高階差分方程 114波動(dòng)方程有限差分的幾個(gè)問題 144.1初始條件 144.2震源函數(shù) 144.3差分方程的穩(wěn)定性及收斂性 154.4頻散問題 155有限差分算法中邊界條件的處理 165.1一維透明邊界條件 175.2二維透明邊界條件 185.3透明邊界條件的差分格式 216簡單模型試算 23 6.2水平層狀模型和高速透鏡體模型 267結(jié)論與建議 30參考文獻(xiàn) 32 1前言1.1地震波場數(shù)值模擬概述地震波場數(shù)值模擬簡單說來，就是已知地下介質(zhì)構(gòu)造及其參數(shù)，再用理論計(jì)算方法來研究地解釋的正確性等方面都有了廣泛的應(yīng)用[1-3]。地震波數(shù)值模擬問題的研究主要包括三個(gè)方面：(1)地震波場數(shù)值模擬原理；(2)地震波場數(shù)值模擬算法；地震波場數(shù)值模擬方法的理論基礎(chǔ)主要是波動(dòng)理論和射線理論。射線理論方法是求出質(zhì)點(diǎn)運(yùn)直接求出波動(dòng)方程的解，其模擬結(jié)果中包含豐富的波動(dòng)信息，模擬結(jié)果較為精確。數(shù)值模擬又可分為波動(dòng)方程的解析解和數(shù)值解。對(duì)簡單的地質(zhì)模型就可以得到其精確解，但對(duì)于復(fù)雜的地質(zhì)模型，只有通過數(shù)值解來說明波在底下介質(zhì)中的傳播。通常解析解只是作為數(shù)值解的一種檢驗(yàn)[4-6]。1.2有限差分法及其與幾種常用數(shù)值模擬方法的比較大[7]。偽譜法在20世紀(jì)七十年代引入數(shù)值模擬計(jì)算領(lǐng)域的，它是利用傅立葉變換來計(jì)算波場的空間導(dǎo)數(shù)，用差分方法來計(jì)算時(shí)間的導(dǎo)數(shù)，可以看成是一種無限階的有限差分法，是傳統(tǒng)有限差分法的動(dòng)問題及氣象預(yù)測等領(lǐng)域有著廣泛的應(yīng)用[8-10]。有限差分法是一種最常用的數(shù)值模擬方法，它是將波動(dòng)方程中波場函數(shù)的空間導(dǎo)數(shù)和時(shí)間導(dǎo)近遠(yuǎn)場及復(fù)雜邊界都有廣泛的適用性，能夠準(zhǔn)確地模擬波在各種介質(zhì)及復(fù)雜結(jié)構(gòu)地層中的傳播規(guī)震學(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]。的差分格式主要是幾種差分格式的組合。我們知道，差分方程的建立首先選擇網(wǎng)格布局和差分形式，然后以有限差分代替無限微分，以差商代替微商，并以差分方程代替微分方程及其邊界條件，最后建立差分方程[12]。在建立差分方程的時(shí)候要注意到兩個(gè)方面。一是合理選擇網(wǎng)格布局與步長。我們將離散后各相鄰離散點(diǎn)之間的距離或者離散化單元的長度稱為步長。在所選定的區(qū)域內(nèi)進(jìn)行網(wǎng)格劃分是差分本文采用的是泰勒級(jí)數(shù)展開法[2,5,13]。即對(duì)連續(xù)的物理量只考慮其在離散的空間位置和離散時(shí)刻的值，然后把方程中的導(dǎo)數(shù)用這些離散的lx其中an是系數(shù)，N表示差分格式的長度，差分的格式是由這兩個(gè)系數(shù)來決定。在實(shí)際應(yīng)用中常用的差分格式有向前差分、向后差分以及中心差分： (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的方程： (2-11)就是二維非均勻介質(zhì)彈性波方程。而在均勻介質(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è)一個(gè)標(biāo)量位的梯度場和一個(gè)向量位的旋度場之和，而空間傳播的波動(dòng)正是無旋運(yùn)動(dòng)和旋轉(zhuǎn)運(yùn)動(dòng)這兩種運(yùn)動(dòng)之和[15-20]令令代入式(2-8)中，取其中的無旋部分則有：然后再消去式中的梯度,就得到了用位移位表示的縱波方程：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é)中討論了差分的原理和幾種常見的差分格式，在這個(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的二階中心差分：關(guān)于x的一階中心差分關(guān)于x的二階中心差分：關(guān)于z的一階中心差分：關(guān)于z的二階中心差分：若令Δx=Δz=h，利用上面關(guān)于x,z,t的中心差分方程就可以得到二維波動(dòng)方程的有限差分方程：式。對(duì)上式繼續(xù)推導(dǎo)可得：如果我們繼續(xù)利用泰勒級(jí)數(shù)展開式，則可以得到：這樣我們就能夠得到時(shí)間域二階空間域四階的波動(dòng)方程有限差分格式：10階上式中k表示時(shí)間方向的離散網(wǎng)格，m表示x方向的離散網(wǎng)格，n表示z方向的離散網(wǎng)格，。需要注意的是，這里同樣也假設(shè)Δ。需要注意的是，這里同樣也假設(shè)Δx=Δz=h，即網(wǎng)格步長相等[5,22-24]。3.3高階差分方程高階。為了方便起見，本文空間步長在x，z軸上相等，則a=vt==vt通過上述類似推xz倒，可以得到不同精度的聲波差分方程。階階12階14階250)，空間步長x和z方向均為1m,采樣時(shí)間間隔100us。所用的震源是雷克子波，子波頻率30Hz。從下圖可以看出，時(shí)間2階，空間2階精度波場出現(xiàn)明顯頻散，隨著空間精度的增高，頻散逐漸減弱。所以提高差分精度可以有效減輕頻散，但計(jì)算量也隨之增加。4波動(dòng)方程有限差分的幾個(gè)問題4.1初始條件在進(jìn)行波動(dòng)方程有限差分?jǐn)?shù)值模擬時(shí)，初始時(shí)刻的波場值是不知道的，而在計(jì)算差分方程時(shí)，只有知道時(shí)間間隔前一時(shí)刻的值才能遞推計(jì)算出后面時(shí)刻的值。因此在計(jì)算開始時(shí)就必須要考慮一般假設(shè)在震源發(fā)生震動(dòng)前，各個(gè)質(zhì)點(diǎn)均處于靜止?fàn)顟B(tài)：初始時(shí)刻波場的速度也同樣為零，即4.2震源函數(shù)震源函數(shù)的給出通常有兩種方式：一種是初值法，就是將理論結(jié)果當(dāng)作初始值給出；另外一種有其優(yōu)缺點(diǎn)[1,26]。在進(jìn)行波動(dòng)方程數(shù)值模擬計(jì)算的過程中，震源問題的處理對(duì)模擬的結(jié)果有著重要的影響。子波的高頻成分頻散就會(huì)更嚴(yán)重。所以要根據(jù)模型的速度及網(wǎng)格間距合理選擇子波主頻。本文采用了雷克子波和高斯子波作為震源函數(shù)[27]。4.3差分方程的穩(wěn)定性及收斂性一個(gè)有實(shí)際應(yīng)用意義的數(shù)值模擬算法必須是穩(wěn)定的，也就是說對(duì)算法的數(shù)值解與解析解的差就說該差分方程的解是穩(wěn)定的，這也就是說差分方程的解是有界的。差分的穩(wěn)定條件根據(jù)差分階不同的[5,28,29]。對(duì)于均勻介質(zhì)，差分方程的解的穩(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ù)值頻散的理論分析研究可以知道，影響頻散的因素有很多。在有限差分離散化之后就引入了差分算子誤差；而時(shí)間間隔的選擇也對(duì)數(shù)值模擬的精度和數(shù)值頻散有一定影說，在進(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ù)值模擬計(jì)算時(shí)，由于受到計(jì)算機(jī)內(nèi)存和計(jì)算速度等因素的制界。當(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]。式中式中V表示聲波的傳播速度，U表示聲壓。我們引入一個(gè)向右行波(k=ω/V)，高衰減區(qū)法是在靠近邊界處引入一個(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)消除回繞波場，這個(gè)方法在理論上是能夠完美的消除回繞波場而不產(chǎn)生任何的反射，但是在實(shí)際應(yīng)用時(shí)，由于這種方法需要進(jìn)行反周期波場樣同時(shí)也會(huì)占用很大內(nèi)存空間[25,34,35]。下面將就本文將要用到的透明邊界條件進(jìn)行討論。5.1一維透明邊界條件對(duì)于一維均勻介質(zhì)波動(dòng)方程：將這個(gè)平面波的方程代入上面的均勻介質(zhì)波動(dòng)方程(5-1)中可以解得：rxar取不當(dāng)會(huì)使左右兩側(cè)的邊界上產(chǎn)生強(qiáng)反射，因此就需要在左右邊界處選取的邊界條件能使r=0。由此可以看出，如果在右行波上加上一個(gè)左行波項(xiàng)由此可以看出，如果在右行波上加上一個(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)榫S波動(dòng)方程如果給定一個(gè)邊界條件xx=a及z=b處會(huì)產(chǎn)生反射。考慮一個(gè)向右傳播的平面波其中θ表示入射角，即平面波前與x軸的夾角。那么這時(shí)反射波為其中r為邊界產(chǎn)生的反射波的反射系數(shù)。則總波場可以表示為以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í)則由式(5-12)可得把式(5-11)代入式(5-15)可得r=0，這時(shí)在右邊界x=a上就不會(huì)產(chǎn)生反射。同理可得在左邊如果取,.并在邊界上令以上這兩種邊界條件都只有很小的邊界反射，但還需差條分法的穩(wěn)定性條分法的穩(wěn)定性系數(shù)看成當(dāng)?shù)奶厥馇闆r下，的系數(shù)就為1/2。則(5-19)式可以寫為由波動(dòng)方程式(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=VΔt/Δx。由于通過泰勒級(jí)數(shù)展開式有因此把式(5-28)和式(5-27)代入式(5-26)可得上式中的微分項(xiàng)可以表示為然后把式(5-30)各式代入式(5-29)即得到右邊界條件的差分格式同理可以得到左邊界條件x=a的差分格式為底邊界條件z=b的差分格式為頂邊界條件z=0的差分格式為上面得到的式(5-34)、式(5-33)、式(5-32)及式(5-31)就是二維均勻介質(zhì)波動(dòng)方程透明邊界條件的差分格式。如下圖5-1(右)所示，可以看出在加入了邊界條件,.后，邊界反射明顯減少，說明透明邊界條件對(duì)邊界反射有很好的吸收效果。圖5-1未加邊界條件(左)與加入邊界條件(右)波場快照對(duì)比圖(所加邊界條件為透明邊界，差分精度為時(shí)間2階空間2階)6簡單模型試算6.1均勻模型震源震動(dòng)時(shí)間5ms，子波為雷克子波，頻率為分精度為時(shí)間2階空間12階，邊界條件為透明邊界條件。從圖中可以看出，波以震源為中心，呈圓弧向外傳播，當(dāng)波傳播至頂邊界時(shí)有微弱的反射波，說明利用透明邊界條件并不能無完全消除人工邊界產(chǎn)生的反射波。另外，從圖中還可以看出有輕微的頻散現(xiàn)象。圖6-2均勻模型波場快照?qǐng)D(差分精度為時(shí)間2階，空間12階，邊界條件為透明邊界條件)6.2水平層狀模型和高速透鏡體模型圖6-3(左)所示的水平層狀模型，模擬網(wǎng)格大小為500×500，網(wǎng)格間距1m；表層速度為雷克子波，子波頻率為30Hz。圖6-3(右)所示的高速透鏡體模型，模擬網(wǎng)格大小為500×500，網(wǎng)格間距1m。透鏡體為速度為3500m/s，周圍介質(zhì)速度為2000m/s。采樣時(shí)間間隔為100μs，震源坐標(biāo)為(250m，mms波選用雷克子波，子波頻率為30Hz。圖6-4中T=40ms時(shí)，波在第一層與第二層的界面上開始傳播，從圖中T=60,80ms中可以清晰的看到第一個(gè)界面的反射波和透射波，T=100ms時(shí)波還沒有傳播到第二層與第三層的界面，T=120ms時(shí)波傳播第三層，從圖中可以看到產(chǎn)生了一個(gè)反射波。頂邊界反射波比較明顯，同時(shí)，隨著波向外傳播，頻散現(xiàn)象也逐漸明顯。波場快照中，出現(xiàn)了波從透鏡體內(nèi)向外傳播時(shí)產(chǎn)生的反射波。T=100ms的波場快照中頻散現(xiàn)象較圖6-3層狀模型(左)和高速透鏡體模型(右)示意圖圖6-4水平層狀模型波場快照?qǐng)D(差分精度為時(shí)間2階，空間12階，邊界條件為透明邊界條件)圖6-5高速透鏡體模型波場快照?qǐng)D(差分精度為時(shí)間2階，空間12階，邊界條件為透明邊界條件)7結(jié)論與建議地震波場的數(shù)值模擬是地震勘探學(xué)的重要研究內(nèi)容之一，是進(jìn)一步研究地震資料的采集、處理和解釋的有效輔助手段。論文以波動(dòng)理論為基礎(chǔ)，對(duì)有限差分模擬的一些關(guān)鍵性問題進(jìn)行了探穩(wěn)定性及頻散問題等。經(jīng)過本次課程報(bào)告，我對(duì)有限差分波場模擬有以下認(rèn)識(shí)： (1)有限差分算法做波動(dòng)方程數(shù)值模擬時(shí)，具有計(jì)算速度快、算法簡單等特征。 (2)有限差分波動(dòng)方程是利用泰勒級(jí)數(shù)展開式展開推導(dǎo)得到的，差分的階數(shù)越高，有限差分的誤差就越小，而有限差分的解就越精確。 (3)有限差分模擬時(shí)，網(wǎng)格空間步長的選擇不僅對(duì)模擬的精確度有影響，而且對(duì)邊界的處理效果的，所以網(wǎng)格步長也必須滿足采樣定理，這樣才能取到完整的子波，以保證計(jì)算結(jié)果不會(huì)出現(xiàn)較大 (4)數(shù)值模擬過程是在特定的區(qū)域內(nèi)進(jìn)行，因此就必然會(huì)產(chǎn)生人工邊界，波傳播到邊界就產(chǎn)生了反射，通過對(duì)幾種模型的模擬表明，本文所用的透明邊界條件和吸收邊界條件對(duì)邊界反射均有較好的吸收效果。本課程報(bào)告的工作在主要是在前人的研究基礎(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á)到了獲得較滿意的模擬結(jié)果，需要采用更優(yōu)的網(wǎng)格剖分方法，以克服存在的缺陷。 (2)對(duì)于模擬中出現(xiàn)的數(shù)值頻散現(xiàn)象，還需要更深入的研究，以便能夠盡可能的減小數(shù)值頻散。 (3)本報(bào)告的算法主要針對(duì)簡單地質(zhì)模型，如果對(duì)于復(fù)雜的介質(zhì)模型，算法還需要進(jìn)一步的研究參考文獻(xiàn)[2]馬在田.計(jì)算地球物理學(xué)概論[M].上海：同濟(jì)大學(xué)出版社，1997.M1988[4]何樵登，熊維鋼.應(yīng)用地球物理教程——地震勘探[M].地質(zhì)出版社，1991[5]吳清嶺，張平，施澤龍.波動(dòng)方程正演模擬即應(yīng)用[J].大慶石油地質(zhì)與開發(fā)，1998，17(3)M學(xué)出版社，2006，7.起伏地表?xiàng)l件下二維地震波場的數(shù)值模擬[J].勘探地球物理進(jìn)展，2005，28(1)[10]吳永剛，吳清嶺.有限差分法彈性波場數(shù)值模擬[J].大慶石油地質(zhì)與開發(fā)，1994，13(3)[11]Alterman.Propagationofelasticwaveinlayeredmediabyfinitedifferencemethods.BulletinoftheSeismologicalSocietyofAmerica,1968;58(1):367~398]Alford,R.M.Kelly,K.R.andBooer,D.M.Accuracyoffinitedifferencemodelingoftheacousticwaveequation.Geophysics,1974;39(6):834~842模擬地震波傳播的大網(wǎng)格快速差分算法[J].地球物理學(xué)報(bào)，1994；37(增刊)：450~454[17]褚春雷，王修田.非規(guī)則三角網(wǎng)格有限差分法地震正演模擬[J].中國海洋大學(xué)學(xué)報(bào)，2005，35 J業(yè)，2004，24(6)：53-56552[20]董良國等.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法[J].地球物理學(xué)報(bào)，2000，43(3):411~419J53[24]CerjanC,KosloffD,KosloffRandReshefM.Anonreflectingboundarycondit-ionfordiscreteacousticandelasticwaveequation.Geophysics,1985,50(4):705~708[25]ClaytonR,EngquistB.Absorbingboundaryconditionsforacousticandelasticwaveequation.BulletinoftheSeismologicalSocietyofAmerica,1977;66(11):1529~1540[26]Reynolds.Boundaryconditionsforthenumericalsolutionofwavepropagationproblem.Geophysics,1978;43(6):1099~1110J物理學(xué)報(bào)，2007，4(4)[30]邢麗.地震波數(shù)值模擬中的吸收邊界條件[J].上海第二工業(yè)大學(xué)學(xué)報(bào)，2006，23(4)[31]DahlainMA.Theapplicationofhighdifferencetothescalarwaveequation.Geophysics,1986;51(1):54~66[32]HigdonRL.Absorbingboundaryconditionforelasticwaves.Geophysics,1991,56[35]宛書金，董敏煌等.橫向各向同性介質(zhì)中的地震波傳播特性[J].石油大學(xué)學(xué)報(bào)，2002,26(1)：3~28本課程報(bào)告所編matlab程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2015.5.20縱波波場模擬%%%%%%%%%%%%%%%%%%%%%%%%%%%%clcclear%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化mode=12;%選擇差分介數(shù)，最高支持16介Model=load('Model.mat');V_model=Model.UniformModel;%速度模型c=load('c.mat');c=c.c;%系數(shù)矩陣cellnum=size(V_model);%網(wǎng)格數(shù)目dt=100*10^(-6);%時(shí)間采樣間隔100usdx=1;dz=1;%空間采樣間隔x0=250;z0=20;f0=30;%震源坐標(biāo)%震源子波頻率S_t=2*10^(-3);T=140*10^(-3);%震源持續(xù)時(shí)間2ms%波場快照時(shí)間%%%%%%%%%%%%%%%%%%%%%%%%%%%%U0=zeros(cellnum(1),cellnum(2));%波場初始值U1=zeros(cellnum(1),cellnum(2));U2=zeros(cellnum(1),cellnum(2));fork=1:T/dt;ifk<S_t/dtU1(z0+1,x0+1)=sin(2*pi*f0*k*dt)/(k*dt);%震源函數(shù)。end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ifmode==2form=mode/2+1:cellnum(2)-mode/2forn=mode/2+1:cellnum(1)-mode/2%2介差分方程a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/2,1)*a^2*(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);endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ifmode==4form=mode/2+1:cellnum(2)-mode/2forn=mode/2+1:cellnum(1)-mode/2%4介差分方程a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+...c(mode/2,2)*a^2*(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);endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ifmode==6form=mode/2+1:cellnum(2)-mode/2forn=mode/2+1:cellnum(1)-mode/2%6介差分方程a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+...c(mode/2,2)*a^2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n))+...c(mode/2,3)*a^2*(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);endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ifmode==8form=mode/2+1:cellnum(2)-mode/2forn=mode/2+1:cellnum(1)-mode/2%8介差分方程a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+...c(mode/2,2)*a^2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n))+...c(mode/2,3)*a^2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n))+...c(mode/2,4)*a^2*(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);endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%form=mode/2+1:cellnum(2)-mode/2ifmode==10%10介差分方程forn=mode/2+1:cellnum(1)-mode/2a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+...c(mode/2,2)*a^2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n))+...c(mode/2,3)*a^2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n))+...c(mode/2,4)*a^2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n))+...c(mode/2,5)*a^2*(U1(m+5,n)+U1(m-5,n)+U1(m,n-5)+U1(m,n+5)-4*U1(m,n))+...2*U1(m,n)-U0(m,n);endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%form=mode/2+1:cellnum(2)-mode/2forn=mode/2+1:cellnum(1)-mode/2ifmode==12%12介差分方程a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+...c(mode/2,2)*a^2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n))+...c(mode/2,3)*a^2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n))+...c(mode/2,4)*a^2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n))+...c(mode/2,5)*a^2*(U1(m+5,n)+U1(m-5,n)+U1(m,n-5)+U1(m,n+5)-4*U1(m,n))+...c(mode/2,6)*a^2*(U1(m+6,n)+U1(m-6,n)+U1(m,n-6)+U1(m,n+6)-4*U1(m,n))+...2*U1(m,n)-U0(m,n);endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%form=mode/2+1:cellnum(2)-mode/2forn=mode/2+1:cellnum(1)-mode/2ifmode==14%14介差分方程a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+...c(mode/2,2)*a^2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n))+...c(mode/2,3)*a^2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n))+...c(mode/2,4)*a^2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n))+...c(mode/2,5)*a^2*(U1(m+5,n)+U1(m-5,n)+U1(m,n-5)+U1(m,n+5)-4*U1(m,n))+...c(mode/2,6)*a^2*(U1(m+6,n)+U1(m-6,n)+U1(m,n-6)+U1(m,n+6)-4*U1(m,n))+...c(mode/2,7)*a^2*(U1(m+7,n)+U1(m-7,n)+U1(m,n-7)+U1(m,n+7)-4*U1(m,n))+...2*U1(m,n)-U0(m,n);endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ifmode==16form=mode/2+1:cellnum(2)-mode/2forn=mode/2+1:cellnum(1)-mode/2%16介差分方程a=V_model(m,n)*dt/dx;U2(m,n)=c(mode/