田雪豐

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

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

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

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

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

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

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

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

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

(1)

其中，t為時間，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三個空間方向和時間方向的離散點(diǎn)序號，Δt為時間剖分步長。

式(2)中，分別在n+1時刻和n-1時刻對U做泰勒展開得：

(5)

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

(6)

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

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

(7)

(8)

(9)

(10)

(11)

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

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

(12)

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

1.2 邊界條件

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

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)中以φ值為自變量，可以求得對應(yīng)的w值，進(jìn)而求得相速度和真實(shí)速度的比值為：

(15)

(a)u分量 (b)ux分量 (c)uz分量圖1 截斷邊界效果驗(yàn)證圖Figure 1 Truncated boundary effect verification

圖2 NAD算法與常規(guī)差分算法頻散曲線對比圖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子波，震源置于模型中央，時間采樣間隔1ms。圖3為采用常規(guī)有限差分法和本文算法正演得到的300ms時的波場快照，常規(guī)有限差分的波場快照(圖3a)中存在嚴(yán)重的數(shù)值頻散現(xiàn)象，這種頻散的本質(zhì)是一種算法誤差而非介質(zhì)本身的聲學(xué)響應(yīng)。利用這種含有誤差的誤差模擬結(jié)果進(jìn)行波傳播機(jī)理的分析會帶來錯誤的分析結(jié)果，誤導(dǎo)生產(chǎn)實(shí)踐。同時，在一些基于波動方程求解的地震數(shù)據(jù)處理技術(shù)中，波場延拓過程中的數(shù)值頻散還會在處理結(jié)果中產(chǎn)生虛假同相軸，降低處理精度?；诒疚乃惴ǖ玫降牟▓隹煺諞]有數(shù)值頻散，模擬精度更高(圖3b)。

3.2 Marmousi模型

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

(a)常規(guī)有限差分法的波場模擬快照 (b)基于NAD算法的波場模擬快照圖3 常規(guī)有限差分法與基于NAD算法對均勻介質(zhì)模型的波場快照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ī)有限差分法波場快照 (b)基于NAD算法波場快照圖5 兩種算法的波場快照對比Figure 5 Comparison of wave field snapshots from two algorithms

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

4 結(jié)論

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