波動(dòng)方程有限差分法頻散衰減方法研究

劉 洋 ， 唐 湘 蓉 （ 成 都 理 工 大 學(xué) 油 氣 藏 地 質(zhì) 及 開(kāi) 發(fā) 工 程 國(guó) 家 重 點(diǎn) 實(shí) 驗(yàn) 室 成都理工大學(xué)地球物院，四川 成都610059）

張璽華 （成都理工大學(xué)沉積地質(zhì)研究院，四川 成都610059）

近年來(lái)，相對(duì)于其他方法，有限差分法顯示了其數(shù)值方法的簡(jiǎn)潔性和生命力。這里通過(guò)修改波動(dòng)方程式來(lái)壓制有限差分中的頻散，新的解法是利用原始波動(dòng)方程式與一個(gè)衰減因子相乘。因此，數(shù)值頻散在大波數(shù)情況下成指數(shù)方式衰減，使得有限差分法的效果像高階一樣，但是計(jì)算成本要小得多。對(duì)模型和真實(shí)數(shù)據(jù)的試驗(yàn)表明了該方法的有效性。由于有限差分法的簡(jiǎn)單性和可擴(kuò)展性，有限差分法已經(jīng)廣泛運(yùn)用于正演和偏移中波場(chǎng)的數(shù)值外推，但是必須特別注意的是一種有限差分法固有的誤差，即數(shù)值頻散。為此，筆者提出了一種簡(jiǎn)單和低計(jì)算成本的改進(jìn)的波動(dòng)方程式用來(lái)壓制數(shù)值頻散。

1 頻散的理論分析

在聲學(xué)中，三維波動(dòng)方程式如下：

式中，v＝v（x，y，z）為空間速度，m/s；x、y、z為空間坐標(biāo)；u為位移，m；t為時(shí)間，s。

在有限差分法中，數(shù)值頻散是一個(gè)不希望出現(xiàn)的現(xiàn)象。先討論一維上行波動(dòng)方程式：

對(duì)空間Z方向進(jìn)行離散，并且用1階中心差分代替微分：

式中，Δz為Z方向的空間采樣間隔，m。

以u(píng)＝exp [i （kz＋ωt） ]（其中，kz為Z方向的波數(shù)；ω為角頻率，rad/s）代入式（3）得：

空間離散引起的相速度與群速度之比為：

式中，O（ ）為高階無(wú)窮小。

由式 （5）可見(jiàn)，空間離散使得相速度低于群速度，并且波數(shù)越大的分量滯后得越多。而且當(dāng)時(shí)間和空間采樣間隔相對(duì)于所使用的有限差分方法過(guò)于粗糙，數(shù)值頻散就會(huì)發(fā)生。從物理意義上看，這種現(xiàn)象的產(chǎn)生是由于不同的頻率具有不同的相速度產(chǎn)生的。而體現(xiàn)在以式 （1）為基礎(chǔ)的有限差分法正演模擬上，數(shù)值頻散主要是由下面的原因造成的：①網(wǎng)格太粗糙；②差分算法不精確。

圖1（a）顯示了一個(gè)由式 （1）推導(dǎo)出的有限差分方法正演結(jié)果，圖1（a）使用的是時(shí)間2階精度差分和空間4階精度差分。雖然在圖1（a）中可以看到比較清晰的圓形波場(chǎng)，但是頻散幾乎覆蓋了環(huán)形波場(chǎng)內(nèi)部的全部區(qū)域，其主要原因是地震波的算法采用的是低階 （2階）精度差分方法和太粗糙的差分計(jì)算網(wǎng)格。

圖1 在各向同性介質(zhì)中250ms處的波場(chǎng)快照

克服上述問(wèn)題的一個(gè)通常的做法是使用較長(zhǎng)的差分算子 （高階有限差分法），另一個(gè)常用的方法是減小計(jì)算網(wǎng)格的尺寸。由于要達(dá)到一定的精確度，每波長(zhǎng)上網(wǎng)格點(diǎn)的最小數(shù)須足以表示一個(gè)波形，算子越長(zhǎng)，需要的網(wǎng)格點(diǎn)數(shù)越少，反之亦然。但必須注意，每波長(zhǎng)的網(wǎng)格點(diǎn)數(shù)不能少于2個(gè)。

圖1顯示的是各種不同空間精度有限差分法合成的波場(chǎng)快照，模型是各向同性介質(zhì)，速度是3500m/s，時(shí)間步長(zhǎng)是0.001s，X方向采樣間隔和Z方向采樣間隔均是10m，有限差分計(jì)算網(wǎng)格是201×201，震源函數(shù)使用的是雷克子波。在圖1（b）與圖1（a）相比，圖1（b）較圖1（a）的頻散現(xiàn)象已經(jīng)明顯改善。當(dāng)更進(jìn)一步，使用8階精度有限差分法計(jì)算同樣的模型，其他參數(shù)保持不變，結(jié)果如圖1（c），頻散現(xiàn)象較圖1（a）和圖1（b）進(jìn)一步改善。圖1（d）使用了2階精度有限差分算法，但是計(jì)算所使用的網(wǎng)格大小是圖1（a）～ （c）所使用的1/4（即網(wǎng)格點(diǎn)數(shù)圖1（d）是圖1（a）～ （c）的4倍），其他參數(shù)保持不變。綜合比較以上4種效果圖，圖1（b）～ （d）展示了比圖1（a）好得多的效果。

在三維波動(dòng)方程式正演中，每一時(shí)間計(jì)算間隔，減小網(wǎng)格會(huì)使計(jì)算量呈3次方增長(zhǎng)；另外，增加有限差分算子的長(zhǎng)度，即用更高階精度的有限差分算法，將使計(jì)算量線(xiàn)性增長(zhǎng)。對(duì)比這兩個(gè)方法，后者具有明顯的優(yōu)勢(shì)。但是，有限差分算子越長(zhǎng)有限差分式的穩(wěn)定性越差。再討論有限差分法的穩(wěn)定條件，根據(jù)柯朗－弗里德里希斯－列維條件，在三維正演模擬中，每個(gè)時(shí)間步長(zhǎng)必須滿(mǎn)足下式：

式中，Δt為時(shí)間采樣間隔，s；vmax為最大模型速度，m/s；Δx、Δy分別為X、Y 方向的空間采樣間隔，m。

由于有式 （6）的影響，在二維波動(dòng)方程式正演下，當(dāng)網(wǎng)格變小，有限差分總計(jì)算量將會(huì)成3次方增長(zhǎng)；在二維波動(dòng)方程式正演下，有限差分總計(jì)算量將會(huì)成4次方增長(zhǎng)。圖1是使用標(biāo)準(zhǔn)波動(dòng)方程式在各向同性介質(zhì)中250ms處的波場(chǎng)快照，震源使用的均是雷克子波，其中圖1（a）、（b）、（c）使用的計(jì)算網(wǎng)格是201×201。圖1（d）使用的是401×401，但是圖1（d）在X方向和Z方向的采樣間隔是圖1（a）～ （c）的1/2，網(wǎng)格大小是圖1（a）～ （c）所使用的1/4。因此，盡管圖1（d）比圖1（a）的效果要好得多，但是圖1（d）需要的計(jì)算量大約是圖1（a）的8倍。如果能夠在2階精度的有限差分法的基礎(chǔ)上壓制頻散，那么就能夠使用比較低的計(jì)算量而得到更好效果的結(jié)果。

2 聲波波動(dòng)方程式優(yōu)化

筆者將式 （1）的解用下式表示：

式中，r（x，y，z）為隨空間變化的非負(fù)實(shí)數(shù)。表面上看，式（9）比式（1）要復(fù)雜得多；然而，式（9）等號(hào)右邊多出的部分只是對(duì)式（1）簡(jiǎn)單的修改。式（9）能夠有效地利用低階（2階或者4階）有限差分法。

圖2驗(yàn)證了不同空間精度有限差分法計(jì)算的單炮正演模型，圖2（a）與圖2（b）均使用的是空間2階有限差分法，區(qū)別在于圖2（a）使用的是式 （1），而圖2（b）使用的是式 （9），圖2（b）中的頻散比圖2（a）中的頻散現(xiàn)象明顯少得多。圖2（c）和圖2（d）均使用式 （1）計(jì)算，但是前者使用的是空間8階精度有限差分，而后者使用的是空間2階精度有限差分，且后者使用的網(wǎng)格在X、Z方向的采樣間隔均是前者的1/2。圖2（b）的效果幾乎可以達(dá)到圖2（c）、圖2（d）的效果，在單炮正演記錄上基本觀測(cè)不到頻散的影響；而圖2（b）的計(jì)算量比圖2（d）的計(jì)算量要小得多。

為了保證式 （9）的解有意義，r（x，y，z）必須滿(mǎn)足下式：

還有一個(gè)簡(jiǎn)單的處理r（x，y，z）值的方法是把r（x，y，z）設(shè)置為一個(gè)常數(shù)R，R的值與差分方法的階數(shù)和差分法計(jì)算所用的網(wǎng)格大小有關(guān)，在該次正演模擬中，R值為0.5。

圖3給出了修改后的波動(dòng)方程式空間2階精度有限差分 （圖2（b））的正演頻譜圖 （圖3（a））和標(biāo)準(zhǔn)波動(dòng)方程式空間2階精度有限差分 （圖2（d））的正演頻譜圖 （圖3（b））。比較這2幅圖，能夠證明修改后的波動(dòng)方程式能量在頻率的分布上幾乎沒(méi)有變化。

式中，X ＝ （x，y，z）表示空間位置；p為式（1）的一個(gè)解；K ＝ （kx，ky，kz）表示空間波數(shù)；kx、ky、kz分別為X、Y、Z方向的波數(shù)。

上面已經(jīng)證明了在用有限差分法模擬波場(chǎng)時(shí)，大波數(shù)的相速度與小波速的相速度不同，而且波數(shù)越大的，相比滯后群速度得越多。因此，給出下式：

式中，g（K）為阻尼因子。相比式（7），式（8）多出一個(gè)阻尼因子g（K），其是平面波解的阻尼函數(shù)，g（K）將在大波數(shù)時(shí)壓制頻散噪聲。下面給出具有式（1）形式的解：

圖2 對(duì)由3層不同介質(zhì) （速度從上到下分別是2700、3100、3500m/s）組成的模型制作的單炮正演記錄

3 結(jié) 語(yǔ)

在波動(dòng)方程式有限差分解法中，數(shù)值頻散是固有的一種非物理數(shù)值噪聲。數(shù)值頻散可以簡(jiǎn)單地通過(guò)大幅度的計(jì)算量來(lái)壓制。筆者從分析空間差分引起的有限差分法數(shù)值頻散出發(fā)，提出了一個(gè)修改的波動(dòng)方程式，在這個(gè)新式中添加了一個(gè)阻尼因子，使得新式的解能夠在大波數(shù)時(shí)高速衰減，從而實(shí)現(xiàn)低階精度有限差分法和較粗糙的網(wǎng)格計(jì)算得到高階精度有限差分算法和精細(xì)網(wǎng)格才能得到的效果。

圖3 修改的波動(dòng)方程式2階精度有限差分法的頻譜圖 （a）與標(biāo)準(zhǔn)波動(dòng)方程式2階精度有限差分法的頻譜圖 （b）對(duì)比

[1]賀振華 .反射地震資料偏移處理與反演方法 [M].重慶：重慶大學(xué)出版社，1989.1～2.

[2]胡建偉，湯懷民 .微分式數(shù)值解法 [M].第2版 .北京：科學(xué)出版社，2007.1～2.

[3]高剛，賀振華，黃德濟(jì)，等 .完全匹配層人工邊界條件中的衰減因子分析 [J].石油物探，2011，50（5）：430～433.

[4]孫成禹，宮同舉，張玉亮，等 .波動(dòng)式有限差分法中的頻散與假頻分析 [J].石油地球物理勘探，2009，44（1）：43～48.

[5]李賓 .橫向各向同性介質(zhì)有限差分法波場(chǎng)模擬方法研究 [D].北京：中國(guó)石油大學(xué)，2009.

[6]肖云飛 .面向波場(chǎng)分析的波動(dòng)式有限差分正演方法 [D].北京：中國(guó)石油大學(xué)，2010.

[6]Li Zhiming.Compensating finite－difference errors in 3－D migration and modeling [J].Geophysics，1991，56 （10）：1650～1660.

[7]Cerjan C，Kosloff D.A none reflection boundary condition for discrete acoustic and elastic wave equation [J].Geophysics，1985，50（4）：705～708.

[8]Tessmer E.Seismic finite－difference modeling with spatially varying time steps [J].Geophysics，2000，65 （4）：1290～1293.