有限差分法不同邊界條件下的數(shù)值模擬

摘 要：文章介紹了地震數(shù)據(jù)處理中所使用的數(shù)值模擬法，對(duì)采用有限差分法所使用不同邊界條件處理方式進(jìn)行了數(shù)值模擬，通過(guò)波場(chǎng)快照直觀的得出了不同的邊界吸收條件的吸收效果，對(duì)結(jié)果進(jìn)行了對(duì)比，分析總結(jié)了各種方法的優(yōu)缺點(diǎn)。

關(guān)鍵詞：數(shù)值模擬；有限差分；邊界條件

隨著近年來(lái)國(guó)家宏觀經(jīng)濟(jì)調(diào)控，經(jīng)濟(jì)增長(zhǎng)的速度逐步減緩，能源行業(yè)受此影響最為嚴(yán)重，許多煤礦是在虧損的情況下生產(chǎn)，直接導(dǎo)致了地質(zhì)行業(yè)投入的減少。物探行業(yè)壓力也越來(lái)越大，物探行業(yè)應(yīng)該抓緊發(fā)展先進(jìn)技術(shù)，提高能源勘探的效率。

在物探行業(yè)中，地震勘探作為一個(gè)重要的手段，發(fā)揮著巨大的作用。數(shù)據(jù)處理作為地震勘探的一大重要環(huán)節(jié)，所采用的各種方法和技術(shù)手段也一直在更新和進(jìn)步。在地震勘探處理方法研究中，地震數(shù)值模擬技術(shù)可以在室內(nèi)完成地震數(shù)據(jù)模型的建立，并對(duì)其地震數(shù)據(jù)進(jìn)行各種方法的處理，查看處理方法的效果和數(shù)據(jù)的好壞，另一方面，地震數(shù)值模擬進(jìn)行正演獲得的數(shù)據(jù)也可以作為反演的基礎(chǔ)進(jìn)行比對(duì)。

在地震數(shù)據(jù)處理的過(guò)程中，如何模擬地震波的傳播便是需要解決的問(wèn)題。在二十世紀(jì)70年代開(kāi)始采用顯示差分格式來(lái)模擬地震波的傳播。由于有限差分法適用條件廣，計(jì)算速度比較快，占用計(jì)算機(jī)內(nèi)存少，編程比較容易實(shí)現(xiàn)，模擬精度相對(duì)較高而得到廣泛應(yīng)用。

但是有限差分法模擬地震波場(chǎng)時(shí)，由于計(jì)算機(jī)運(yùn)算核心的限制，有限差分方法只能得到有限的數(shù)據(jù)點(diǎn)，地震波動(dòng)方程只能是在有限差分方程中求得近似解，這時(shí)就考慮到人工邊界問(wèn)題，如果不對(duì)邊界進(jìn)行處理，波在通過(guò)邊界時(shí)會(huì)產(chǎn)生反射，因此我們希望對(duì)添加的邊界進(jìn)行處理來(lái)消除這些反射。

在20世紀(jì)70年代，地球物理學(xué)界陸續(xù)采用了不同的邊界條件來(lái)實(shí)現(xiàn)削弱地震波在通過(guò)邊界時(shí)的反射，比如reynold邊界、clayton邊界、cerjan邊界，以及后來(lái)提出的PML層邊界條件，每種邊界條件都在不同程度上實(shí)現(xiàn)了地震波通過(guò)邊界時(shí)的衰減。

為了驗(yàn)證以上邊界條件在數(shù)值模型的效果，在文章中，我們?cè)O(shè)計(jì)了一些簡(jiǎn)單的數(shù)值模型，給出了不同的邊界條件，通過(guò)波形在通過(guò)不同邊界條件時(shí)反射進(jìn)行比較，觀察每種方法衰減反射的效果。

1 reynold邊界條件

文章首先采用的邊界條件為reynold邊界條件。這是ALBERT C.REYNOLDS于1978年提出的一種邊界條件的處理方式。地下介質(zhì)中地震波的傳播規(guī)律可以近似地用聲波方程描述，對(duì)于二維模型，首先我們給出二維情況下聲波方程，

在此方程中S（t）為震源函數(shù)，使用了雷克子波，其中時(shí)間域采樣間隔為0.001s，采樣點(diǎn)數(shù)為100點(diǎn)，總的時(shí)間長(zhǎng)度為0.1s，則頻率域采樣間隔為10hz。v（x，z）為介質(zhì)在點(diǎn)（x，z）處的縱波速度，u為描述速度位或者壓力的波場(chǎng)。震源函數(shù)波形如圖1：

建立的數(shù)值模型是一個(gè)均勻速度模型，其中網(wǎng)格數(shù)為300x300，速度vp=3000，其中上邊界采用自由反射邊界，其余三個(gè)邊界采用reynold邊界，可以由圖2得到的波場(chǎng)快照看出波在通過(guò)邊界時(shí)反射得到了有效衰減。

2 clayton and engquist邊界條件

采用同樣的模型參數(shù)，當(dāng)下邊界使用自由反射界面，其余三面采用clayton and engquist邊界條件，得到的波場(chǎng)快照如圖3所示。

兩種邊界條件得到的單炮記錄如圖4所示。

由所得到的波場(chǎng)快照可以直觀的看出，兩種邊界處理方法都有明顯效果，當(dāng)?shù)卣鸩ㄍㄟ^(guò)未經(jīng)過(guò)處理的邊界時(shí)，地震波的反射現(xiàn)象十分嚴(yán)重，產(chǎn)生一系列反射波并多次反射，對(duì)數(shù)據(jù)形成干擾，但是經(jīng)過(guò)處理的邊界時(shí)反射得到了明顯的削弱。

3 Cerjan吸收邊界條件

下面所使用的數(shù)值模型是時(shí)間二階、空間二階精度二維彈性波模擬，采用的震源仍然是雷克子波，其時(shí)間域采樣間隔為0.001s，主頻為30。模型采用300x300的網(wǎng)格均勻常速度介質(zhì)，其中P波速度Vp=2000，S波速度Vs=1400，邊界條件為上邊界采用自由反射邊界，其余邊界采用Cerjan邊界條件。圖5為全波震源得到的波場(chǎng)快照，上圖為x分量，下圖為z分量。

由波場(chǎng)快照可以看出，地震波在傳播時(shí)，由于縱波速度比橫波速度快，所產(chǎn)生的波形為兩個(gè)同心圓，外層為直達(dá)縱波，內(nèi)層為直達(dá)橫波。在經(jīng)過(guò)邊界時(shí)，由于上邊界采用了自由邊界條件，可以看出地震波的反射現(xiàn)象也是比較明顯，其他邊界由于添加了衰減吸收函數(shù)，所產(chǎn)生的反射比較微弱。同時(shí)由于該吸收采用的衰減函數(shù)比較簡(jiǎn)單，計(jì)算時(shí)所需要的時(shí)候也比較少，計(jì)算起來(lái)速度比較快。

4 結(jié)束語(yǔ)

由以上圖可以看出，在進(jìn)行地震模型數(shù)值模擬時(shí)，所采用的邊界條件都達(dá)到了一定了效果，但是也有一些不足，采用reynold邊界條件和Clayton條件時(shí)，其邊界吸收對(duì)比自由邊界效果比較明顯，但是仍然具有一定的反射，采用cerjan邊界時(shí)，波在邊界反射后產(chǎn)生的頻散也比較嚴(yán)重，除了以上方法之外，PML邊界條件是近年來(lái)提出的吸收效果比較好的一種邊界處理?xiàng)l件，其原理是在構(gòu)建的地質(zhì)模型外圍添加一層吸收邊界，這種邊界能迅速吸收邊界中的反射，效果比較好，所需要的計(jì)算量也比較少，使用也比較普遍。由于采用這種邊界的模型文獻(xiàn)較多，在此不再進(jìn)行敘述。

參考文獻(xiàn)

[1]賈躍瑋.基于MATLAB的地震正演模型實(shí)現(xiàn)[J].石油工業(yè)計(jì)算機(jī)應(yīng)用，2009，2：10-13+65.

[2]左瑩.基于高階交錯(cuò)網(wǎng)格的有限差分地震波場(chǎng)數(shù)值模擬[D].長(zhǎng)安大學(xué)，2009.

[3]王開(kāi)燕，周妍，劉丹，等.有限差分波動(dòng)方程正演模擬中的吸收邊界條件[J].當(dāng)代化工，2014，5：752-755.

[4]李桂花，朱光明，馬德堂.井間地震數(shù)據(jù)處理方法——以大慶某采油廠井間地震勘探為例[J].石油地球物理勘探，2011，1：76-82+164+171.

[5]朱光明，李桂花，程建遠(yuǎn).煤礦巷道內(nèi)地震勘探的數(shù)值模擬[J].煤炭學(xué)報(bào)，2008，11：1263-1267.

作者簡(jiǎn)介：韓浩（1992-），男，山東濟(jì)寧人，在讀研究生，專業(yè)為地球探測(cè)與信息技術(shù)，主要從事地球物理勘探數(shù)據(jù)處理方法的研究。