一階聲波方程時間四階精度差分格式的偽譜法求解

唐懷谷 何兵壽

（中國海洋大學海底科學與探測技術(shù)教育部重點實驗室，山東青島266100）

一階聲波方程時間四階精度差分格式的偽譜法求解

唐懷谷 何兵壽*

（中國海洋大學海底科學與探測技術(shù)教育部重點實驗室，山東青島266100）

在地震波場數(shù)值模擬中，偽譜法不產(chǎn)生由空間網(wǎng)格離散引起的數(shù)值頻散，而常規(guī)偽譜法用于求解時間二階精度差分格式時，則會受到時間差分精度較低的影響而產(chǎn)生數(shù)值頻散。本文基于一階聲波方程，提出將差分格式的時間差分精度增至四階，并利用偽譜法求解，從而在避免由空間網(wǎng)格離散引起的數(shù)值頻散的同時，降低由時間網(wǎng)格離散引起的數(shù)值頻散。此外，與時間二階精度差分格式偽譜法相比，時間四階精度差分格式偽譜法的穩(wěn)定性條件更為寬松，進而可通過增大時間網(wǎng)格步長提高計算效率。

一階聲波方程 數(shù)值模擬 偽譜法 時間高階差分格式 數(shù)值頻散 穩(wěn)定性條件

1 引言

在地震波場數(shù)值模擬中，從聲波方程或彈性波方程中求取空間導數(shù)的方法主要有偽譜法［1，2］和有限差分法［3］等。與有限差分法相比，偽譜法不會產(chǎn)生對空間導數(shù)項的差分引起的數(shù)值頻散，因此對空間導數(shù)的求取精度更高［4］。尤其是在子波頻率較高及模型中存在低速層的情況下，偽譜法在數(shù)值模擬的精度上表現(xiàn)出明顯優(yōu)勢。

偽譜法雖然不產(chǎn)生空間差分引起的數(shù)值頻散，但其數(shù)值模擬精度仍受兩方面制約：一是由空間網(wǎng)格離散和傅里葉變換引起的吉布斯效應造成的波前畸變；另一方面是對時間導數(shù)項進行差分引起的數(shù)值頻散。目前，偽譜法交錯網(wǎng)格求解技術(shù)［5］可有效壓制吉布斯效應造成的波前畸變。因此，進一步提高偽譜法數(shù)值模擬精度的關(guān)鍵在于降低由時間差分引起的數(shù)值頻散。

地震分辨率決定于子波的頻帶寬度和主頻［6］。然而子波主頻越高，由時間差分引起的數(shù)值頻散越嚴重，也就意味著信噪比越低。因此，在高分辨率地震勘探中，必須降低由時間差分引起的數(shù)值頻散。減小時間網(wǎng)格步長是壓制頻散的方式之一，然而減小時間網(wǎng)格步長意味著正演到相同時間情況下增加了迭代次數(shù)，導致計算效率降低。壓制由時間差分引起的數(shù)值頻散的另一途徑是采用高階時間精度的差分格式進行求解［7，8］。

在一階聲波方程和一階彈性波方程中，可將時間高階導數(shù)項的求取轉(zhuǎn)化為對空間高階導數(shù)項的求取，從而得到任意階時間精度的差分格式［3，9］。然而利用交錯網(wǎng)格有限差分法對高階時間精度差分格式進行求解，一方面對空間高階導數(shù)的計算量非常大，另一方面其空間差分引起的數(shù)值頻散相比時間差分引起的數(shù)值頻散要嚴重得多。本文通過偽譜法求取空間導數(shù)項，在簡化高階空間導數(shù)求取的同時避免了空間差分引起的數(shù)值頻散。本文基于二維情況下的一階聲波方程給出了時間四階精度差分格式的偽譜法求解方法，該方法可以拓展到三維情況以及一階彈性波方程。此外，本文證明了時間四階精度差分格式偽譜法比時間二階精度差分格式偽譜法的穩(wěn)定性條件寬松，因此可以通過增大時間步長來提高運算效率。

2 一階聲波方程及其時間2M階差分近似

在不考慮體力的情況下，各向同性介質(zhì)一階聲波速度—應力方程可寫成

式中：D x、Dz分別為x和z方向的空間微分算子，分別表示為P為應力分量；v x和v z分別為x和z方向的速度分量。將其表示為矩陣形式

應用交錯網(wǎng)格中心差分的思想求解一階速度應力方程，根據(jù)一元函數(shù)的Taylor展式，得到聲波方程的時間2M階差分近似［9］

3 時間四階精度差分格式的偽譜法

3.1 時間四階精度差分格式及交錯網(wǎng)格偽譜法求解

令式（5）中的M為2，則

由式（4）可得

以方程形式表示式（6），得到時間四階精度差分格式如下

利用交錯網(wǎng)格偽譜法［5］求解式（8）中的空間導數(shù)，其計算公式為

式（9）中：q指代速度分量或應變分量，q F為q的空間二維傅氏變換；k x和k z分別為x和z方向的波數(shù)；Δx和Δz分別為x和z方向的空間網(wǎng)格步長；±表示交錯網(wǎng)格中不同分量之間的位置關(guān)系。將式（8）中的空間導數(shù)利用式（9）求解，得到時間四階差分精度交錯網(wǎng)格偽譜法計算公式為

3.2 時間四階精度偽譜法的穩(wěn)定性條件和計算效率

在均勻介質(zhì)中，式（10）所示的時間四階精度差分格式的交錯網(wǎng)格偽譜法求解格式的穩(wěn)定性條件（附錄B）是

而在均勻介質(zhì)中，時間二階精度差分格式的交錯網(wǎng)格偽譜法求解格式的穩(wěn)定性條件為

根據(jù)式（11）和式（12）的計算結(jié)果，表1列舉了幾種參數(shù)條件下，時間二階精度偽譜法和時間四階精度偽譜法滿足穩(wěn)定性條件的最大時間網(wǎng)格步長Δtmax。

表1 偽譜法在二維均勻介質(zhì)中滿足穩(wěn)定性條件的最大時間網(wǎng)格步長

在非均勻介質(zhì)中，穩(wěn)定性條件要苛刻一些。從表1可看出，與時間二階精度偽譜法相比，時間四階精度偽譜法對Δt的取值范圍有顯著的放寬。

在內(nèi)存消耗方面，與時間二階精度偽譜法相比，時間四階精度偽譜法需要對三階空間微分項進行計算。由于式（8）中對不同空間微分項的計算相互獨立，且每個微分項只需計算一次，在計算空間微分項的過程中可以利用相同的內(nèi)存空間進行存儲。因此時間四階精度偽譜法相比時間二階精度偽譜法并不增加內(nèi)存的消耗。

由于對規(guī)模為N x×N z的二維波場進行一次傅氏變換或反傅氏變換需要進行4（N x+N z）N x N z次乘法運算和4（N x+N z）N x N z次加法運算（附錄A），因此偽譜法數(shù)值模擬的計算量主要集中在對傅氏變換和反傅氏變換的計算上，可以將正、反傅氏變換的計算次數(shù)作為評定偽譜法計算量的指標。

表2 二維情況下時間二階差分精度偽譜法和時間四階差分精度偽譜法進行一次迭代的傅氏變換次數(shù)

數(shù)值模擬總的計算量是迭代一次的計算量和迭代次數(shù)的乘積。從表2可見，在二維情況下，時間四階精度偽譜法迭代一次的計算量約為時間二階精度偽譜法的2.2倍；通過式（11）和式（12）的比較，若時間步長均取理想條件下的最大值，則時間四階精度差分格式的時間步長可以取到時間二階差分格式的2.84倍。因此可以通過增大時間網(wǎng)格步長的方式減少迭代次數(shù)，進而提高時間四階精度差分格式偽譜法的計算效率。

3.3 時間四階精度偽譜法的頻散特性

時間四階精度偽譜法的頻散關(guān)系可表示為

圖1 時間二階精度和時間四階精度偽譜法頻散曲線

從圖1上看，時間二階精度偽譜法的相速度超前于實際地震波速度，1／G越大則相速度與實際地震波速度之比越大；時間四階精度偽譜法的相速度落后于實際地震波速度，1／G越大則相速度與實際地震波速度之比越小。

1／G越大則意味著一次振動包含的時間網(wǎng)格點數(shù)越少。對于子波中的低頻成分，由于周期較長，一次振動包含的網(wǎng)格數(shù)較多，1／G較小，從圖1中可看出，在1／G較小的情況下，時間二階精度偽譜法和時間四階精度偽譜法的低頻成分的相速度與實際地震波速度差異均比較小；反之，對于高頻成分，則1／G較大。從圖1可見，在1／G較大的情況下，時間二階精度偽譜法與時間四階精度偽譜法相比，前者的相速度與實際地震波速度的差異遠大于后者。因此，在子波主頻較高的情況下，時間四階精度偽譜法的頻散特性優(yōu)于時間二階精度偽譜法。

3.4 簡化的非分裂PML吸收邊界條件

考慮到時間高階差分格式存在混合偏導數(shù)，不易將波場按傳播方向進行拆分，因此采用非分裂PML吸收邊界條件［10］。本文對非分裂PML吸收邊界條件做進一步簡化，采用了一種衰減系數(shù)不隨時間變化的非分裂PML邊界條件。加載了邊界條件的波動方程表示為

式中Ω為衰減系數(shù)，根據(jù)波場傳播方向和邊界位置之間的關(guān)系，將Ω分為四部分：非衰減區(qū)域、x方向邊界區(qū)域（左邊界區(qū)域和右邊界區(qū)域）、z方向邊界區(qū)域（上邊界區(qū)域和下邊界區(qū)域）、角部邊界區(qū)域

這種簡化的PML吸收邊界的衰減系數(shù)不隨時間變化，只與位置和PML邊界的厚度有關(guān)，因此衰減系數(shù)的計算簡單且計算量小。加載了邊界條件的波動方程不需要對波場按傳播方向進行劃分，而針對整體波場進行衰減，因此內(nèi)存占用小。

4 模型算例及對比分析

在圖2a所示的模型中，空間網(wǎng)格點數(shù)為200×200，空間網(wǎng)格步長Δx=Δz=5m。令時間網(wǎng)格步長Δt=0.8ms，震源位于（100，50），震源子波函數(shù)為

令子波主頻f=80 Hz，鑲嵌30層吸收邊界，分別利用一階聲波方程時間二階精度，空間十二階精度的有限差分法、時間二階差分精度的偽譜法以及時間四階差分精度的偽譜法進行正演模擬，記錄第400ms的波前快照。

從圖2b、圖2c、圖2d三者的對比來看，由于波速較低，子波頻率較高，圖2b所示的有限差分法正演模擬結(jié)果中存在明顯的數(shù)值頻散，尤其是由空間差分引起的數(shù)值頻散非常嚴重，圖2c和圖2d所示的偽譜法模擬結(jié)果中不含空間差分引起的數(shù)值頻散。從圖2c和圖2d的對比來看，在子波頻率較高的情況下，時間二階精度差分格式的偽譜法正演模擬結(jié)果中仍然存在明顯的由時間差分引起的數(shù)值頻散，將時間差分精度提高至四階能夠有效減弱頻散效應。

抽取圖2b、圖2c、圖2d中橫坐標650m、縱坐標900～1250m位置處的振動圖（圖3）。進一步觀察波前形態(tài)，并與震源子波形態(tài)對比。可見有限差分法受空間離散影響造成的數(shù)值頻散在波形上滯后于波前，二階時間精度偽譜法受時間差分精度低的影響造成的數(shù)值頻散在波形上超前于波前，四階時間精度的偽譜法的數(shù)值頻散效應顯著降低，總體上數(shù)值頻散的波形略滯后于波前。這與本文3.3節(jié)所述時間二階精度偽譜法和時間四階精度偽譜法的頻散特征相符。

圖2 含有一個反射層的層狀模型及三種不同正演方法得到的波前快照

由圖3可見，盡管時間四階精度偽譜法相比時間二階精度偽譜法的頻散特性有了改善，但在子波頻率較高的情況下，受時間差分截斷誤差的影響，時間四階精度偽譜法的波前振動圖仍然呈現(xiàn)一定的滯后特性。此外，子波旁瓣兩側(cè)仍有擾動，這是由吉布斯效應引起的誤差。雖然使用交錯網(wǎng)格技術(shù)可壓制由吉布斯效應導致的波前畸變［5］，但由于偽譜法只能利用離散頻譜求取空間導數(shù)的本質(zhì)特性，將離散的波數(shù)域波場反變換回空間域必然會產(chǎn)生吉布斯效應。如圖4所示，在子波頻率較高、波場速度較低以及空間網(wǎng)格步長較大的情況下，偽譜法求取空間導數(shù)的精度受吉布斯效應影響較為嚴重。換言之，可通過降低子波主頻、提高地震波速度以及減小空間網(wǎng)格步長三個方面減弱吉布斯效應。

圖3 抽取圖2b、圖2c、圖2d的橫坐標650m、縱坐標900～1250m位置處的振動圖及震源雷克子波波形

圖4 不同參數(shù)條件下偽譜法求取空間導數(shù)產(chǎn)生的吉布斯效應比較

為進一步驗證時間四階精度偽譜法相對于時間二階偽譜法在精度上的優(yōu)勢，對圖5e所示水平層狀模型進行正演模擬，正演模擬的參數(shù)為Δt=0.5ms，Δx=Δz=3.5，炮點位置（1250，0），f=80Hz，記錄長度為700ms。通過對地震記錄的比較，可觀察到在子波頻率較高時，時間二階精度偽譜法的時間頻散明顯強于時間四階精度偽譜法，利用時間四階精度偽譜法得到的地震記錄分辨率更高。

在Marmousi（局部）模型中分別應用上述三種方法進行正演模擬，通過圖6所示地震記錄的對比可見，時間四階精度差分格式的偽譜法能夠有效克服頻散效應，得到的地震記錄具有較高的分辨率。時間二階精度偽譜法和時間四階精度偽譜法的計算效率如表3所示。

圖5 有限差分法、時間二階精度和時間四階精度偽譜法水平層狀模型地震記錄

圖6 利用不同正演方法對Marmousi（局部）模型得到的地震記錄

表3 時間二階精度偽譜法和時間四階精度偽譜法的串行程序計算效率

5 結(jié)束語

通過將一階聲波方程的時間差分精度提高至四階，在參數(shù)相同的情況下，偽譜法的數(shù)值模擬精度比時間差分精度為二階時有顯著提升，有效降低了由時間差分引起的數(shù)值頻散。

時間四階精度差分格式偽譜法迭代一次的計算量大于時間二階精度差分格式，但時間四階精度差分格式的偽譜法具有穩(wěn)定性條件寬松的優(yōu)勢，因此可采用更大的時間步長以降低迭代次數(shù)，在一定程度上提高其計算效率。

感謝中國石油大學趙強博士對本文提出了建議，感謝中國海洋大學林琦博士為本文編制了圖件。

［1］ Kreiss H O，Oliger J.Comparison of accurate methods for the integration of hyperbolic equations.Tellus，1972，24（3）：199-215.

［2］ Kosloff D D，Baysal E.Forward modeling by a Fourier method.Geophysics，1982，47（10）：1402-1412.

［3］ 董良國，馬在田，曹景忠等.一階彈性波方程交錯網(wǎng)格高階差分解法.地球物理學報，2000，43（3）：411-419.Dong Liangguo，Ma Zaitian，Cao Jingzhong et al.A staggered-grid high-order difference method of oneorder elastic wave equation.Chinese Journal of Geophysics，2000，43（3）：411-419.

［4］ Fornberg B.The pseudospectral method：Comparisons with finite differences for the elastic wave equation.Geophysics，1987，52（4）：483-501.

［5］ 杜增利，徐峰，高宏亮.虛譜法交錯網(wǎng)格地震波場數(shù)值模擬.石油物探，2010，49（5）：430-437.Du Zengli，Xu Feng and Gao Hongliang.Pseudo spectral seismic wavefield simulation with staggered grid.GPP，2010，49（5）：430-437.

［6］ 李慶忠.走向精確勘探的道路.北京：石油工業(yè)出版社，1993.

［7］ 董良國，李培明.地震波傳播數(shù)值模擬中的頻散問題.天然氣工業(yè)，2004，24（6）：53-56.Dong Liangguo and Li Peiming.Dispersive problem in seismic wave propagation numerical modeling.Natural Gas Industry，2004，24（6）：53-56.

［8］ 吳國忱，王華忠.波場模擬中的數(shù)值頻散分析與校正策略.地球物理學進展，2005，20（1）：58-65.Wu Guochen and Wang Huazhong.Analysis of numerical dispersion in wavefield simulation.Progress in Geophysics，2005，20（1）：56-65.

［9］ 牟永光，裴正林.三維復雜介質(zhì)地震數(shù)值模擬.北京：石油工業(yè)出版社，2005.

［10］ 秦臻，任培罡，姚姚等.彈性波正演模擬中PML吸收邊界條件的改進.地球科學：中國地質(zhì)大學學報，2009，34（4）：658-664.Qin Zhen，Ren Peigang，Yao Yao et al.Improvement of PML absorbing boundary conditions in elastic wave forward modeling.Earth Science—Journal of China University of Geosciences，2009，34（4）：658-664.

［11］ Marchuk G I.Methods of Numerical Mathematics.New York：Springer-Verlag，1975，22-30.

［12］ Irving R S.Integers，Polynomials and Rings：A Course in Algebra.New York：Springer Science ＆ Business Media，2003.

附錄A 二維傅氏變換的計算量

在一維傅氏變換中，長度為N的數(shù)組x n=（x1，x2，…，x N）變換后的結(jié)果為

由式（A-1）計算一個Xm需要N次復數(shù)乘法運算和N次復數(shù)加法運算，因此長度為N的數(shù)組進行傅氏變換需要N2次復數(shù)乘法運算和N2次復數(shù)加法運算。反傅氏變換公式為

其計算量與傅氏變換的計算量相同。

對一個規(guī)模為N x×N z的矩陣進行二維傅氏變換，其過程由兩次一維傅氏變換組成，即對矩陣按行（列）進行一維傅氏變換，所有的行（列）做一維傅氏變換后仍然按行（列）構(gòu)成新矩陣，再對該矩陣按列（行）進行一維傅氏變換。

按行進行一維傅氏變換，即對N x個長度為N z的一維數(shù)組進行傅氏變換，計算量為次復數(shù)乘法運算和復數(shù)加法運算；同理，按列進行一維傅氏變換，即對N z個長度為N x的一維數(shù)組進行傅氏變換，計算量為次復數(shù)乘法運算和復數(shù)加法運算。因此對一個規(guī)模為N x×N z的矩陣進行二維傅氏變換或反變換的總計算量為次復數(shù)乘法運算和復數(shù)加法運算。

計算機完成一次復數(shù)乘法運算需要進行4次實數(shù)乘法運算和2次實數(shù)加法運算，完成一次復數(shù)加法運算需要2次實數(shù)加法運算。因此計算機完成一個規(guī)模為N x×N z的矩陣進行二維傅氏變換或反變換需要進行4（N x+N z）N x N z次實數(shù)加法運算和4（N x+N z）N x N z次實數(shù)乘法運算。

附錄B 偽譜法求解一階聲波方程時間2M階差分格式的穩(wěn)定性條件

對式（5）進行空間傅氏變換，得到

式中G（kx，kz）是Q（x，z）的傅氏變換，將U（t，kx，kz）的系數(shù)矩陣記為A，即

則式（B-1）可寫成

記系數(shù)矩陣A的特征值為λ，則根據(jù)傅里葉分析方法［11］可得到式（B-3）穩(wěn)定的條件是A的最大模長的特征值的模小于2。

A的3個特征值為

式中Dkx、D kz為空間微分算子的差分格式D x和D z的空間傅里葉變換。利用偽譜法求解式（B-3），則

γ的模在奈奎斯特波數(shù)處取到最大值。即當k x=時，就有

因此，利用偽譜法求解式（B-3）穩(wěn)定的條件是

當M=1時，得到時間二階精度的穩(wěn)定性條件，即式（11）；當M=2時，得到時間四階精度的穩(wěn)定性條件，即式（12）。

附錄C 一階聲波方程偽譜法頻散特性

偽譜法求解時間二階精度聲波方程的差分格式為

設一階聲波方程P分量簡諧波形式的解析解為

式中：θ為P的傳播方向x與z方向的夾角；ω為角頻率（ω=2πf），f為頻率；k為波數(shù)λ為波長，則有為相速度。

將式（C-2-1）分別代入式（C-1-2）和式（C-1-3）中，得到v x和v z的解析解為

將式（C-2）代入式（C-1-1）中，得到

該式等號兩端同時除以ω，進一步整理得到

將式（C-2）代入式（C-1-2）中，得到

同理，將式（C-2）代入式（C-1-3）中，得到

在式（C-6）和式（C-7）的等號兩端同時除以ω，做進一步整理，對應地可得到式（C-4）和式（C-5）。因此，時間二階精度差分格式偽譜法滿足式（C-5）所描述的頻散關(guān)系。

求取時間四階精度差分格式的頻散特性，將式（C-2）代入式（8-1）中，得到

經(jīng)過化簡，得到

進一步可寫為

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.011

唐懷谷，何兵壽.一階聲波方程時間四階精度差分格式的偽譜法求解.石油地球物理勘探，2017，52（1）：71-80.

1000-7210（2017）01-0071-10

*山東省青島市中國海洋大學海底科學與探測技術(shù)教育部重點實驗室，266100。Email：hebinshou＠ouc.edu.cn

本文于2015年12月23日收到，最終修改稿于2016年12月18日收到。

本項研究受國家自然科學基金項目（41174087）和“863”項目（2013AA064201）資助。

（本文編輯：朱漢東）

唐懷谷 碩士研究生，1992年生；2014年本科畢業(yè)于中國海洋大學勘查技術(shù)與工程專業(yè)，獲工學學士學位；現(xiàn)在該校海洋地球科學學院攻讀地球探測與信息技術(shù)專業(yè)碩士學位。