唐懷谷 何兵壽

（中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島266100）

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

唐懷谷 何兵壽*

（中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島266100）

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

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

1 引言

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

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

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

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

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

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

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

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

3 時(shí)間四階精度差分格式的偽譜法

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

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

由式（4）可得

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 時(shí)間四階精度偽譜法的頻散特性

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

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

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

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

3.4 簡(jiǎn)化的非分裂PML吸收邊界條件

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

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

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

4 模型算例及對(duì)比分析

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

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

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

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

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

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

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

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

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

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

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

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

表3 時(shí)間二階精度偽譜法和時(shí)間四階精度偽譜法的串行程序計(jì)算效率

5 結(jié)束語(yǔ)

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

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

感謝中國(guó)石油大學(xué)趙強(qiáng)博士對(duì)本文提出了建議，感謝中國(guó)海洋大學(xué)林琦博士為本文編制了圖件。

［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］ 董良國(guó)，馬在田，曹景忠等.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法.地球物理學(xué)報(bào)，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］ 杜增利，徐峰，高宏亮.虛譜法交錯(cuò)網(wǎng)格地震波場(chǎ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］ 董良國(guó)，李培明.地震波傳播數(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］ 吳國(guó)忱，王華忠.波場(chǎng)模擬中的數(shù)值頻散分析與校正策略.地球物理學(xué)進(jìn)展，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］ 牟永光，裴正林.三維復(fù)雜介質(zhì)地震數(shù)值模擬.北京：石油工業(yè)出版社，2005.

［10］ 秦臻，任培罡，姚姚等.彈性波正演模擬中PML吸收邊界條件的改進(jìn).地球科學(xué)：中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，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 二維傅氏變換的計(jì)算量

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

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

其計(jì)算量與傅氏變換的計(jì)算量相同。

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

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

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

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

對(duì)式（5）進(jìn)行空間傅氏變換，得到

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

則式（B-1）可寫成

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

A的3個(gè)特征值為

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

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

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

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

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

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

設(shè)一階聲波方程P分量簡(jiǎn)諧波形式的解析解為

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

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

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

該式等號(hào)兩端同時(shí)除以ω，進(jìn)一步整理得到

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

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

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

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

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

進(jìn)一步可寫為

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.011

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

1000-7210（2017）01-0071-10

*山東省青島市中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，266100。Email：hebinshou＠ouc.edu.cn

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

本項(xiàng)研究受國(guó)家自然科學(xué)基金項(xiàng)目（41174087）和“863”項(xiàng)目（2013AA064201）資助。

（本文編輯：朱漢東）

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