唐懷谷 何兵壽
(中國(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)定性條件
在地震波場(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)算效率。
在不考慮體力的情況下,各向同性介質(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]
令式(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ì)算公式為
在均勻介質(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ì)算效率。
時(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í)間二階精度偽譜法。
考慮到時(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)存占用小。
在圖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ì)算效率
通過將一階聲波方程的時(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é)位。