基于NAD算法的聲波方程時(shí)間四階差分解法

田雪豐

(中國(guó)煤炭地質(zhì)總局地球物理勘探研究院，河北 涿州 072750)

數(shù)值模擬是研究復(fù)雜介質(zhì)中地震波傳播機(jī)理的重要工具[1-3]，高精度的正演模擬結(jié)果也可用于指導(dǎo)地震資料的采集方案設(shè)計(jì)和處理、反演流程設(shè)置[4-8]?；诟唠A差分離散技術(shù)的有限差分算法[9-11]在以往的科學(xué)研究和工業(yè)生產(chǎn)中發(fā)揮了巨大作用。但隨著業(yè)界對(duì)模擬精度要求的不斷提高和地震勘探所面臨的問(wèn)題的復(fù)雜程度的不斷增加，這種方法由于其固有缺陷的限制而越來(lái)越難以滿(mǎn)足學(xué)術(shù)界和工業(yè)界的需求，其缺陷主要表現(xiàn)在：當(dāng)波動(dòng)方程差分階數(shù)過(guò)小或地下介質(zhì)中的地震波傳播速度過(guò)低或地震波的頻率過(guò)高時(shí)，常規(guī)有限差分算法往往會(huì)產(chǎn)生嚴(yán)重的數(shù)值頻散問(wèn)題，這種數(shù)值頻散不是由介質(zhì)本身的彈性性質(zhì)引起的，而是由模擬算法本身引起的，其本質(zhì)是一種誤差，該誤差會(huì)使地震波的振幅、頻率和相位產(chǎn)生畸變，嚴(yán)重時(shí)還會(huì)產(chǎn)生虛假同相軸，因此帶有數(shù)值頻散的波動(dòng)方程正演模擬結(jié)果是無(wú)法指導(dǎo)生產(chǎn)實(shí)踐的，也無(wú)法用于地震波傳播機(jī)理的分析，必須研究并采用適當(dāng)?shù)募夹g(shù)消除這種誤差，確保模擬結(jié)果的精度。

目前，業(yè)界解決數(shù)值頻散壓制問(wèn)題的主要思路有五：①減小時(shí)間和空間網(wǎng)格的剖分步長(zhǎng)[12]；②提高波動(dòng)方程導(dǎo)數(shù)項(xiàng)的差分離散階數(shù)[13-15]；③優(yōu)化偏導(dǎo)數(shù)項(xiàng)的差分系數(shù)[16-17]；④在數(shù)值頻散產(chǎn)生之后采用適當(dāng)算法對(duì)其進(jìn)行壓制[18-21]；⑤研究并采用新技術(shù)提高差分離散精度，減小頻散誤差[22-25]。

在常規(guī)速度模型和低頻條件下，上述第一和第二種思路對(duì)數(shù)值頻散的壓制非常有效，實(shí)現(xiàn)起來(lái)也最為簡(jiǎn)單，但其缺陷也非常明顯，表現(xiàn)為：空間剖分步長(zhǎng)的減小意味著網(wǎng)格點(diǎn)的指數(shù)級(jí)增加，時(shí)間剖分步長(zhǎng)的減小意味著成倍增加計(jì)算時(shí)間步數(shù)以獲取相同記錄長(zhǎng)度的合成記錄，兩者會(huì)導(dǎo)致計(jì)算量和內(nèi)存需求的指數(shù)級(jí)增大。提高波動(dòng)方程導(dǎo)數(shù)項(xiàng)的差分離散階數(shù)也會(huì)成倍增加計(jì)算量，并且當(dāng)離散階數(shù)增加到一定程度時(shí)，差分階數(shù)的提高對(duì)于數(shù)值頻散壓制效果的提升會(huì)變得越來(lái)越緩慢。第三種思路一般不會(huì)增加計(jì)算量和內(nèi)存需求，但對(duì)頻散壓制效果的提升有限。第四種思路以通量校正傳輸(FCT)技術(shù)[19-21]為代表，它基于先產(chǎn)生后壓制的思路在波場(chǎng)延拓的過(guò)程中對(duì)波動(dòng)方程的各個(gè)分量進(jìn)行校正，當(dāng)參數(shù)設(shè)置合理時(shí)，F(xiàn)CT技術(shù)具有明顯的數(shù)值頻散壓制效果，但FCT技術(shù)的缺陷表現(xiàn)在：漫射通量和反漫射通量的求取以及反漫射通量的校正均需要消耗巨大的計(jì)算量，同時(shí)漫射因子和反漫射因子的設(shè)置沒(méi)有明確的標(biāo)準(zhǔn)，難以合理設(shè)置，此外FCT技術(shù)在壓制數(shù)值頻散的同時(shí)還可能對(duì)有效信號(hào)造成損害，降低波場(chǎng)延拓的精度。

本文方法屬于第五種思路，即在保持空間和時(shí)間網(wǎng)格大小、差分階數(shù)等參數(shù)不變的前提下，依據(jù)數(shù)值頻散的產(chǎn)生機(jī)理研究新方法來(lái)避免波場(chǎng)延拓過(guò)程中的數(shù)值頻散，獲取更高精度的模擬結(jié)果。具體為：利用NAD算子[23-28]求解各網(wǎng)格點(diǎn)處的高階空間微分項(xiàng)，實(shí)現(xiàn)聲波方程的差分離散，得到聲波方程的水平梯度場(chǎng)和垂直梯度場(chǎng)，在此基礎(chǔ)上推導(dǎo)基于NAD算子的時(shí)間四階精度的差分格式，實(shí)現(xiàn)數(shù)值頻散的壓制。數(shù)值算例表明，本文基于NAD算子的時(shí)間四階精度的波動(dòng)方程差分格式對(duì)波場(chǎng)中的空間頻散和時(shí)間頻散均具有明顯壓制作用。

本文方法和第四種思路的區(qū)別在于：第四種思路本質(zhì)上是先產(chǎn)生后壓制，而本文方法是在事前避免其產(chǎn)生。本文方法和前四種思路的區(qū)別在于：前四種思路只能壓制波動(dòng)方程數(shù)值求解過(guò)程中產(chǎn)生的空間頻散，不能壓制時(shí)間頻散，而本文算法可以實(shí)現(xiàn)這兩種頻散的同時(shí)壓制。

1 基于NAD算子的聲波方程時(shí)間四階差分精度解法

1.1 差分格式推導(dǎo)

各向同性介質(zhì)中的聲波方程為：

(1)

其中，t為時(shí)間，x、y、z為直角坐標(biāo)，v為縱波傳播速度，u為質(zhì)點(diǎn)位移。

在方程(1)中引入新變量U、W、P：

(2)

顯然可得到以下關(guān)系式：

(3)

(4)

其中i、j、k、n分別為x、y、z三個(gè)空間方向和時(shí)間方向的離散點(diǎn)序號(hào)，Δt為時(shí)間剖分步長(zhǎng)。

式(2)中，分別在n+1時(shí)刻和n-1時(shí)刻對(duì)U做泰勒展開(kāi)得：

(5)

將式(5)中的兩個(gè)方程求和并整理可得：

(6)

式(3)、式(4)、式(6)表明下一時(shí)刻的波場(chǎng)值和其梯度值可通過(guò)求解前一時(shí)刻波場(chǎng)的空間偏導(dǎo)數(shù)得到，因此聲波方程的求解問(wèn)題變成求解分量高階空間偏導(dǎo)數(shù)問(wèn)題。

將式(2)代入(6)式，可得：

(7)

(8)

(9)

(10)

(11)

將式(8—11)代入式(7)即可得到基于NAD算子的三維聲波方程時(shí)間4階精度差分格式。

經(jīng)推導(dǎo)，式(7)的穩(wěn)定性條件為：

(12)

其中vmax為計(jì)算空間中的最大地震波速度。

1.2 邊界條件

邊界條件是波動(dòng)方程正演的重要內(nèi)容，PML邊界條件在許多波動(dòng)方程正演中得到了成功應(yīng)用[29-31]，本文也采用PML邊界條件吸收入射到計(jì)算邊界的外行波。圖1為PML吸收邊界條件的效果驗(yàn)證圖，圖中入射到邊界的外行波被邊界完全吸收，說(shuō)明PML邊界條件能夠有效解決本文差分算法的邊界偽反射問(wèn)題。

2 數(shù)值頻散分析

二維情況下可假設(shè)式(1)的平面諧波解為：

u(x,z,t)=u0exp(i(wt-kxx-kzz))

(13)

將上式代入式(7)中，整理可得式(7)的頻散關(guān)系為：

(14)

式中，v為數(shù)值模擬得到的地震波傳播速度，v0為模型速度，k為波數(shù)，w=kvΔt，φ=kΔx，φx=φcosθ，φz=φsinθ，θ為波的空間傳播方向與x方向的夾角。

式(14)中以φ值為自變量，可以求得對(duì)應(yīng)的w值，進(jìn)而求得相速度和真實(shí)速度的比值為：

(15)

(a)u分量 (b)ux分量 (c)uz分量圖1 截?cái)噙吔缧Ч?yàn)證圖Figure 1 Truncated boundary effect verification

圖2 NAD算法與常規(guī)差分算法頻散曲線(xiàn)對(duì)比圖Figure 2 Comparison of frequency dispersion curves from NAD algorithm and conventional difference algorithm

3 模型算例

3.1 均應(yīng)介質(zhì)模型

二維均勻介質(zhì)模型的縱波速度為2 000m/s，模型大小2 000m×2 000m，按照10m×10m的空間網(wǎng)格大小差分離散該模型并進(jìn)行正演模擬，模擬所用的震源為主頻40Hz的Ricker子波，震源置于模型中央，時(shí)間采樣間隔1ms。圖3為采用常規(guī)有限差分法和本文算法正演得到的300ms時(shí)的波場(chǎng)快照，常規(guī)有限差分的波場(chǎng)快照(圖3a)中存在嚴(yán)重的數(shù)值頻散現(xiàn)象，這種頻散的本質(zhì)是一種算法誤差而非介質(zhì)本身的聲學(xué)響應(yīng)。利用這種含有誤差的誤差模擬結(jié)果進(jìn)行波傳播機(jī)理的分析會(huì)帶來(lái)錯(cuò)誤的分析結(jié)果，誤導(dǎo)生產(chǎn)實(shí)踐。同時(shí)，在一些基于波動(dòng)方程求解的地震數(shù)據(jù)處理技術(shù)中，波場(chǎng)延拓過(guò)程中的數(shù)值頻散還會(huì)在處理結(jié)果中產(chǎn)生虛假同相軸，降低處理精度?；诒疚乃惴ǖ玫降牟▓?chǎng)快照沒(méi)有數(shù)值頻散，模擬精度更高(圖3b)。

3.2 Marmousi模型

Marmousi模型如圖4所示，分別采用常規(guī)有限差分法和本文算法進(jìn)行正演模擬，模擬所用的參數(shù)為：震源為主頻40Hz的Ricker子波，震源位于地表中間位置處，空間網(wǎng)格大小10m×10m，時(shí)間采樣間隔1ms，全排列接收，接收點(diǎn)位于地表，道間距10m。圖5為上述兩種算法得到的t=1 400ms時(shí)的波場(chǎng)快照，圖6為兩種算法得到的合成地震記錄，對(duì)比表明在上述時(shí)間和空間網(wǎng)格尺寸條件下，常規(guī)有限差分算法會(huì)產(chǎn)生嚴(yán)重的數(shù)值頻散問(wèn)題，模擬精度底，而本文算法得到的波場(chǎng)清晰干脆，沒(méi)有數(shù)值頻散。

(a)常規(guī)有限差分法的波場(chǎng)模擬快照 (b)基于NAD算法的波場(chǎng)模擬快照?qǐng)D3 常規(guī)有限差分法與基于NAD算法對(duì)均勻介質(zhì)模型的波場(chǎng)快照Figure 3 Homogeneous medium model wave field snapshots from conventional finite difference method and NAD algorithm

圖4 Marmousi速度模型圖Figure 4 Marmousi velocity model

(a)常規(guī)有限差分法波場(chǎng)快照 (b)基于NAD算法波場(chǎng)快照?qǐng)D5 兩種算法的波場(chǎng)快照對(duì)比Figure 5 Comparison of wave field snapshots from two algorithms

(a)常規(guī)有限差分法合成記錄 (b)基于NAD算法合成記錄圖6 兩種算法的合成記錄對(duì)比Figure 6 Comparison of synthetic seismic records from two algorithms

4 結(jié)論

數(shù)值頻散是有限差分求解波動(dòng)方程時(shí)產(chǎn)生的一種誤差，這種誤差會(huì)使地震波的頻帶和波形發(fā)生畸變，還會(huì)產(chǎn)生虛假波場(chǎng)。常規(guī)有限差分算法中，當(dāng)?shù)卣鸩l率較高或差分網(wǎng)格較大時(shí)，這種頻散尤其嚴(yán)重。本文針對(duì)這一問(wèn)題，推導(dǎo)了三維各向同性介質(zhì)聲波方程近似解析離散化數(shù)值模擬方法，相較于常規(guī)差分算法，本文算法在求解空間微分時(shí)同時(shí)利用了更多的點(diǎn)的波場(chǎng)信息，得到的空間微分更為準(zhǔn)確，有效地壓制了數(shù)值頻散誤差，得到了更高精度的波場(chǎng)快照和模擬結(jié)果。同時(shí)，由于本文算法在大網(wǎng)格尺寸條件下能得到無(wú)頻散的波場(chǎng)模擬結(jié)果，因此相同模型條件下，可采用大網(wǎng)格對(duì)模擬空間進(jìn)行差分離散，減少網(wǎng)格數(shù)，這在客觀(guān)上降低了模擬所需的內(nèi)存需求，也提高了計(jì)算效率。