聲波有限差分?jǐn)?shù)值模擬

——基于完全匹配層吸收邊界條件聲波有限差分?jǐn)?shù)值模擬

——基于完全匹以下討論聲波模擬的有限差分方法地震資料的數(shù)值模擬作用：

1）模擬地震記錄，檢驗(yàn)處理結(jié)果的好壞、處理方法的有效性；

2）正演模擬可以作為反演研究的基礎(chǔ)。固體彈性介質(zhì)簡(jiǎn)化為聲學(xué)介質(zhì)：研究地震波傳播問(wèn)題及地震成像方法時(shí)，為了方便求解，只研究縱波的波場(chǎng)特征及成像方法。這種做法是對(duì)實(shí)際問(wèn)題的良好近似。因?yàn)榈乇砀浇嬖诘退賻?，地震反射記錄中橫波信息的能量非常微弱。以下討論聲波模擬的有限差分方法地震資料的數(shù)值模擬作用：固體彈內(nèi)容提要一、聲波方程二、聲波方程的差分近似三、吸收邊界條件內(nèi)容提要一、聲波方程

一、聲波方程聲波方程：通常，介質(zhì)密度相對(duì)于其速度變化很小，可以近似地將密度看作常數(shù)：即其中：為聲壓為密度為速度為震源為梯度算子一、聲波方程聲波方程：通常，介質(zhì)密度相對(duì)于其速度變化很小，二、聲波方程的差分近似三維常密度聲波方程：其中：為頻帶有限的震源子波，震源位于處。對(duì)微分方程進(jìn)行差分近似時(shí)，主要利用了Taylor級(jí)數(shù)展開(kāi)理論：（2-1）二、聲波方程的差分近似三維常密度聲波方程：其中：為頻帶有限二、聲波方程的差分近似設(shè)是一個(gè)多元函數(shù)，是它的一個(gè)自變量的差分近似時(shí)，先將記為為了表述方便，當(dāng)研究（2-2）（2-3）由（2-2）和（2-3）得（2-4）二、聲波方程的差分近似設(shè)是一個(gè)多元函數(shù)，是它的一個(gè)自變量的差二、聲波方程的差分近似（2-4）中分別取并做誤差截?cái)嗟茫毫睿?-5）二、聲波方程的差分近似（2-4）中分別取并做誤差截?cái)嗟茫毫睿ǘ?、聲波方程的差分近似?-5）可以寫(xiě)成（2-6）求解方程組（2-6）得（2-7）所以二、聲波方程的差分近似（2-5）可以寫(xiě)成（2-6）求解方程組的六階誤差截?cái)啾磉_(dá)式：用類似的方法可以求得的四階誤差截?cái)啾磉_(dá)式：（2-8）的二階誤差截?cái)啾磉_(dá)式：（2-9）（2-10）二、聲波方程的差分近似的六階誤差截?cái)啾磉_(dá)式：用類似的方法可以求得的四階誤差截?cái)啾磉_(dá)二、聲波方程的差分近似取x、y、z、t方向的差分網(wǎng)格間距分別為、、、令，在點(diǎn)考慮差分近似以時(shí)間二階空間四階誤差截?cái)酁槔校簩⑸厦媸阶哟肼暡ǚ匠蹋?-1）就可以得到對(duì)應(yīng)的差分方程：，震源項(xiàng)取為并將速度取為二、聲波方程的差分近似取x、y、z、t方向的差分網(wǎng)格間距二、聲波方程的差分近似（2-11）二維情況下，上式變?yōu)椋海?-12）其中，二、聲波方程的差分近似（2-11）二維情況下，上式變?yōu)椋海?整理（2-12）式，得整理（2-12）式，得在邊界處需要降階二、聲波方程的差分近似在邊界處需要降階二、聲波方程的差分近似三、吸收邊界條件（圖3-1）（圖3-2）在用有限差分法進(jìn)行聲波數(shù)值模擬時(shí)，若采用自然邊界條件，即令各個(gè)時(shí)間層邊界處的聲壓值為零，在波場(chǎng)沒(méi)有到達(dá)人為邊界時(shí)計(jì)算結(jié)果是正確的（圖3-1）；而當(dāng)波場(chǎng)超出人為邊界時(shí)，計(jì)算結(jié)果中混入了人為邊界的反射（圖3-2）。三、吸收邊界條件（圖3-1）（圖3-2）在用有限差分法進(jìn)行聲（圖3-2）（圖3-3）工程中的許多波動(dòng)問(wèn)題，其計(jì)算求解域往往很大，有時(shí)甚至是無(wú)界的。但實(shí)際情況下，我們總是在一個(gè)有限的區(qū)域內(nèi)進(jìn)行求解，因此，需要在所限定的區(qū)域的邊界上引入吸收邊界條件，從而最大限度地降低由于人為劃定的邊界而造成邊界反射。吸收邊界的建立對(duì)于波動(dòng)方程的數(shù)值模擬起到致關(guān)重要的作用，它將直接影響的計(jì)算結(jié)果的穩(wěn)定性和準(zhǔn)確性。圖3-2沒(méi)有加吸收邊界；圖3-3加入了吸收邊界。（圖3-2）（圖3-3）工程中的許多波動(dòng)問(wèn)題，其計(jì)算求解域往三、吸收邊界條件

一些用于消除人工邊界反射的方法包括：

Smith于1974年提出多次反射法；

Clayton及Engquist于1977年提出的旁軸近似法；

近年來(lái)，Berenger針對(duì)電磁波傳播情況，給出了一種高效的完全匹配層（PML）吸收邊界條件，并在理論上證明該方法可以完全吸收來(lái)自各個(gè)方向、各種頻率的電磁波，而不產(chǎn)生任何反射。目前已將完全匹配層的思想引入到了彈性介質(zhì)中，并建立了彈性介質(zhì)中的完全匹配層吸收邊界條件。以下利用聲波方程，導(dǎo)出聲波的完全匹配層的控制方程，并作為吸收邊界條件應(yīng)用于數(shù)值模擬。實(shí)際結(jié)果表明這是一種非常有效的方法。三、吸收邊界條件一些用于消除人工邊界反射的方法包聲波方程：其中（3-1）等價(jià)運(yùn)動(dòng)方程連續(xù)性方程（3-2）（3-3）質(zhì)點(diǎn)速度等價(jià)（3-6）（3-5）（3-4）聲波方程：其中（3-1）等價(jià)運(yùn)動(dòng)方程連續(xù)性方程（3-2）（3三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（3-9）（3-8）（3-7）完全匹配層邊界條件：當(dāng)時(shí)，上面方程就是聲波方程；可以證明，當(dāng)時(shí)，在方向上起到的衰減作用，以下給出證明。三、吸收邊界條件

——PerfectlyMat三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

將（3-7）、（3-8）、（3-9）變換到頻域（3-12）（3-11）（3-10）對(duì)x、y、z做下面坐標(biāo)拉伸變換：其中同理（3-13）（3-14）三、吸收邊界條件

——PerfectlyMat三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（3-13）、（3-14）代入（3-10）、（3-11）得（3-12）（3-11）（3-10）（3-10）、（3-11）、（3-12）變回到時(shí)域：（3-15）（3-14）（3-13）三、吸收邊界條件

——PerfectlyMat（3-9）（3-8）（3-7）（3-15）（3-14）（3-13）（3-6）（3-5）（3-4）兩組式子形式一樣坐標(biāo)變換因此，若聲波方程的解為PML方程的解應(yīng)該為（3-9）（3-8）（3-7）（3-15）（3-14）（3-三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

有了上面的推導(dǎo)作為基礎(chǔ)，下面以二維情況為例，說(shuō)明PML的衰減作用。設(shè)二維聲波方程的一個(gè)平面簡(jiǎn)諧波解為：若此簡(jiǎn)諧波在PML介質(zhì)中傳播，則對(duì)應(yīng)的聲壓解應(yīng)該為：坐標(biāo)變換兩個(gè)解的相位完全相同，且振幅比為：（3-16）三、吸收邊界條件

——PerfectlyMat三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（3-16）（3-16）可以解釋PMLABC效用和其參數(shù)的效用：當(dāng)時(shí)，，振幅不發(fā)生衰減。當(dāng)時(shí)，可以使，振幅發(fā)生衰減，且和分別控制著各自方向上的衰減。時(shí)，振幅無(wú)衰減，隨著它們的增大，振幅呈指數(shù)衰減。三、吸收邊界條件

——PerfectlyMat以下給出（3-7）（3-8）（3-9）交錯(cuò)網(wǎng)格離散表達(dá)式：（3-7）以為例，在進(jìn)行離散（3-8）(3-17)(3-18)以下給出（3-7）（3-8）（3-9）交錯(cuò)網(wǎng)格離散表達(dá)式：（對(duì)（3-17）和（3-18）做以下變量替換，替換時(shí)需要注意在（3-17）和（3-18）中，為了得到從到的聲壓波場(chǎng)的遞推公式，(3-17)(3-19)(3-18)(3-20)對(duì)（3-17）和（3-18）做以下變量替換，替換時(shí)需要注意在(3-9)(3-19)(3-20)交錯(cuò)網(wǎng)格時(shí)間層遞推(3-9)(3-19)(3-20)交錯(cuò)網(wǎng)格時(shí)間層遞推以下以二維情況為例說(shuō)明PML吸收邊界的實(shí)現(xiàn)方法：三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

（圖3-4）如圖3-4，中間的部分為進(jìn)行有限差分模擬的區(qū)域，為了加入PML吸收邊界條件，需要在它的周圍加入PML介質(zhì)，圖中，L代表左邊界，R代表右邊界，U代表上邊界，D代表下邊界，RD代表右下角邊界，LD代表左下角邊界，RU代表右上角邊界，LU代表左上角邊界。垂直方向設(shè)置為有衰減，這時(shí)可以讓從邊界處為零值開(kāi)始向外逐漸增大到某個(gè)值。水平方向設(shè)置為有衰減，這時(shí)可以讓從邊界處為零值開(kāi)始向外逐漸增大到某個(gè)值。在D和U處，水平方向設(shè)置為無(wú)衰減，即在L和R處，垂直方向設(shè)置為無(wú)衰減，即在RD、LD、RU、LU處，水平和垂直方向均設(shè)置為有衰減。

以下以二維情況為例說(shuō)明PML吸收邊界的實(shí)現(xiàn)方法：三、吸收邊界圖3-4中的細(xì)實(shí)線，代表每個(gè)數(shù)組的次邊界，粗實(shí)線代表每個(gè)數(shù)組的邊界，數(shù)值模擬時(shí)，每完成一次時(shí)間層遞推，一個(gè)數(shù)組里的邊界值要和相臨數(shù)組的次邊界值發(fā)生交換。三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

如圖3-5所示，下面的兩個(gè)“回”字矩形代表 時(shí)刻時(shí)間層上的數(shù)據(jù)，一個(gè)數(shù)組里的邊界值要和相臨數(shù)組的次邊界值發(fā)生交換，上面兩個(gè)矩形代表由 時(shí)刻的時(shí)間層數(shù)據(jù)遞推得到的 時(shí)刻的時(shí)間層數(shù)據(jù)，完成計(jì)算后，再進(jìn)行邊界和次邊界的數(shù)據(jù)交換，生成新的數(shù)組邊界數(shù)據(jù)，依此類推完成整個(gè)時(shí)間層的計(jì)算。（圖3-5）圖3-4中的細(xì)實(shí)線，代表每個(gè)數(shù)組的次邊界，粗實(shí)線代表每個(gè)數(shù)組三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

如圖3-6，為了程序?qū)崿F(xiàn)時(shí)便于操作數(shù)組，將D和U合并在一起存放的一個(gè)數(shù)組里，將L和R合并在一起存放的一個(gè)數(shù)組里，將RD、LD、RU、LU合并在一起存放的一個(gè)數(shù)組里。這種做法的另一好處就是，它將圖1-3-4中周圍的PML介質(zhì)四周也連接起來(lái)，這樣不需要在每次計(jì)算完一個(gè)時(shí)間層時(shí)對(duì)圖1-2-4的最外邊界賦零值，而是讓波場(chǎng)繼續(xù)傳播，多通過(guò)一倍的PML介質(zhì)，衰減更充分，卻沒(méi)有增加任何計(jì)算量。（圖3-6）（圖3-4）三、吸收邊界條件

——PerfectlyMat三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

下面兩個(gè)圖說(shuō)明了PML介質(zhì)中速度的一種選取方法。三、吸收邊界條件

——PerfectlyMat以下給出一些使用了PMLABC的二維聲波數(shù)值模擬的算例時(shí)間切片，震源坐標(biāo)為震源處使用的是一個(gè)峰值頻率為，峰值振幅為1，的Ricker子波。

介質(zhì)速度為延時(shí)為圖3-7、圖3-8、圖3-9、圖3-10，給出了四張時(shí)間分別為三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondition

以下給出一些使用了PMLABC的二維聲波數(shù)值模擬的算例時(shí)間（圖3-7）t=0.1s（圖3-8）t=0.25s（圖3-9）t=0.04s（圖3-10）t=0.55s（圖3-7）t=0.1s（圖3-8）t=0.25s（圖3-9三、吸收邊界條件

——PerfectlyMatchedLayerAbsorbingBoundaryCondi