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

2017-10-23 22:36:12唐懷谷何兵壽

石油地球物理勘探 2017年1期

唐懷谷何兵壽

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

唐懷谷何兵壽*

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

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

一階聲波方程數值模擬偽譜法時間高階差分格式數值頻散穩定性條件

1 引言

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

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

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

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

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

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

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

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

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

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

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

由式（4）可得

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

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

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

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

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

而在均勻介質中，時間二階精度差分格式的交錯網格偽譜法求解格式的穩定性條件為

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

表1 偽譜法在二維均勻介質中滿足穩定性條件的最大時間網格步長

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

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

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

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

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

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

時間四階精度偽譜法的頻散關系可表示為

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

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

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

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

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

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

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

4 模型算例及對比分析

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

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

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

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

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

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

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

圖4 不同參數條件下偽譜法求取空間導數產生的吉布斯效應比較

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

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

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

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

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

5 結束語

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

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

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

［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］董良國，馬在田，曹景忠等.一階彈性波方程交錯網格高階差分解法.地球物理學報，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］杜增利，徐峰，高宏亮.虛譜法交錯網格地震波場數值模擬.石油物探，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］李慶忠.走向精確勘探的道路.北京：石油工業出版社，1993.

［7］董良國，李培明.地震波傳播數值模擬中的頻散問題.天然氣工業，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］吳國忱，王華忠.波場模擬中的數值頻散分析與校正策略.地球物理學進展，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］牟永光，裴正林.三維復雜介質地震數值模擬.北京：石油工業出版社，2005.

［10］秦臻，任培罡，姚姚等.彈性波正演模擬中PML吸收邊界條件的改進.地球科學：中國地質大學學報，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的數組x n=（x1，x2，…，x N）變換后的結果為

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

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

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

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

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

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

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