劉 洋 ， 唐 湘 蓉 （ 成 都 理 工 大 學 油 氣 藏 地 質(zhì) 及 開 發(fā) 工 程 國 家 重 點 實 驗 室 成都理工大學地球物院，四川 成都610059）

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

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

1 頻散的理論分析

在聲學中，三維波動方程式如下：

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

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

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

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

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

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

式中，O（ ）為高階無窮小。

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

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

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

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

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

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

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

由于有式 （6）的影響，在二維波動方程式正演下，當網(wǎng)格變小，有限差分總計算量將會成3次方增長；在二維波動方程式正演下，有限差分總計算量將會成4次方增長。圖1是使用標準波動方程式在各向同性介質(zhì)中250ms處的波場快照，震源使用的均是雷克子波，其中圖1（a）、（b）、（c）使用的計算網(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）需要的計算量大約是圖1（a）的8倍。如果能夠在2階精度的有限差分法的基礎上壓制頻散，那么就能夠使用比較低的計算量而得到更好效果的結(jié)果。

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語

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

圖3 修改的波動方程式2階精度有限差分法的頻譜圖 （a）與標準波動方程式2階精度有限差分法的頻譜圖 （b）對比

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

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

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

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

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

[6]肖云飛 .面向波場分析的波動式有限差分正演方法 [D].北京：中國石油大學，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.