適于聲波方程數值模擬的時間-空間域有限差分算子改進線性算法

梁文全， 賴永明， 姜 琳， 吳超凡

(龍巖學院資源工程學院，福建 龍巖 364000)

?

適于聲波方程數值模擬的時間-空間域有限差分算子改進線性算法

梁文全， 賴永明， 姜 琳， 吳超凡

(龍巖學院資源工程學院，福建 龍巖 364000)

有限差分法廣泛應用于地震波場的數值延拓，確定合適的有限差分算子以減小數值頻散是有限差分法的一個重要研究內容。近年來為了進一步抑制數值頻散和增加時間步長，新的有限差分模板得到了應用，對于此，前人使用泰勒展開方法和最小二乘方法確定有限差分算子系數。本文在以前工作的基礎上，使用改進的線性方法確定新模板的有限差分系數，并與傳統模板線性方法進行對比；通過頻散分析和正演模擬驗證出新模板線性方法能夠更好地保持頻散關系，在相同的精度下效率提高了一倍，從而說明了改進的線性方法的有效性。

聲波正演; 時間-空間域; 有限差分格式; 線性方法

0 引言

數值模擬在實踐中得到了廣泛的應用[1-2]。其中有限差分方法因為計算效率高、所需內存小、實現簡單而廣泛應用于地震正演研究[3-10],同時構成了逆時偏移成像和全波形反演的基礎[11-12]。

數值頻散是有限差分研究的關鍵問題，直接影響著有限差分方法在波動方程中的應用。數值頻散是對波動方程的時間和空間偏導數的離散化造成的，使得相速度變成了網格間距的函數，從而導致地震波的數值相速度不等于地球介質的真實相速度，使得波場模擬的精度降低。一般來說，如果存在時間頻散，則高頻的相速度增大；如果存在空間頻散，則高頻的相速度減小[13]。

起初，波動方程空間偏導數的有限差分算子在空間域中確定。但地震波在時間域和空間域同時傳播，因此在時間-空間域確定空間偏導數的有限差分算子能同時抑制時間頻散和空間頻散。Liu等[14]通過平面波理論和泰勒展開方法在時間-空間域確定有限差分算子系數。Yang等[15]提出和改進了近似解析離散化方法(NADM 方法)。Zhang等[7]提出使用模擬退火方法確定有限差分系數并指出合適的誤差容限的重要性。筆者曾提出根據震源頻率范圍、波速和網格間距確定一個合適的波數范圍，在此圍內使用線性方法確定有限差分系數，并取得了良好效果[16]。為進一步壓制頻散和增加時間步長，Liu等[17]也提出新的有限差分模板并用最小二乘方法確定了相應的有限差分系數。本文使用線性方法確定新有限差分模板對應的系數，該方法比最小二乘方法計算系數的效率更高，可在保證有限差分算子精度的同時提高其計算效率。

1 有限差分算子確定方法

對于二維的聲波方程,

(1)

其中:p是壓力波場或位移波場，在(x,z)兩個方向上的二階空間偏導數，使用相同有限差分格式加以計算，并采用新的有限差分模板格式[17][圖1(b)]，可以得到

(2)

其中:cm和d表示有限差分系數[當d=0時新差分模板格式退化到圖1(a)所示的傳統差分模板格式];h表示網格長度(m);τ表示時間步長(s)，v表示波速(m/s);M表示差分格式的算子長度，其中

(3)

對于時間偏導數，采用二階有限差分近似，得到

(4)

將式(2)和式(4)帶入式(1)，可得

(5)

對式(5)進行傅里葉變換，可得

(6)

令

(7)

根據震源頻率、速度和網格間距確定需要滿足頻散關系的波數范圍占總波數的比例[16]，

(8)

其中:f是震源最高頻率;v表示地震波速度;h表示空間采樣間距。

在式(8)確定的波數范圍內，使用之前提出的線性方法[16]，假設M+2個均勻分布的波數點滿足式(6)所示的時間-空間域頻散關系,得到方程(9)所示的線性方程組：

(9)

由于矩陣的條件數較大，解方程(9)不能確定合適的有限差分算子系數，因此需要對線性方法進行改進，使之適應新有限差分模板。通過研究分析，假設M+1個均勻分布的波數點滿足式(6)所示的時間-空間域頻散關系, 并令d=r2/4，得到方程(10)所示的線性方程組：

(10)

方程(10)中，令有限差分算子系數d遍歷數組c中元素，并對頻散誤差進行比較，得到頻散誤差最小的一組有限差分算子系數。由于線性方法的計算效率極高，即使進行了3次計算，其效率仍遠高于最小二乘方法。在方程(9)和(10)中，k(i)(i=1,2,…,M+1)均勻分布于0和rπ/h之間，而r通過式(8)確定。ki的分量表示為ki,l，其中l=(x,z)。我們稱方程(10)為新模板改進線性方法，當d=0時，退化到傳統模板線性方法。

圖1 傳統及新的有限差分模板Fig.1 (a) The traditional and new FD stencil

由方程(6)可以得到，

(11)

令式(11)左邊

(12)

由式(11)可知g≤0。由于sin2(ωτ/2)≤1，要保證穩定性，則必須滿足條件

(13)

整理得到

(14)

一般來說，隨著波數(頻率)的增加頻散誤差增大，令

(15)

將式(12)和(15)代入不等式(14)，得到新模板線性方法穩定性條件

(16)

圖2是新模板改進線性方法不同有限差分算子長度M、不同r值時的穩定性分析曲線圖。從圖2中可以看到，M越大、r越大越容易造成不穩定，因此當r值較大時(波速較大時)，應該采用較小的M值，這和實際的模擬情況一致。

圖2 新模板線性方法的穩定性條件Fig.2 Stability condition of the new linear method

2 頻散分析

二維聲波方程數值頻散δ(θ)定義為[18]

(17)

1-δ(θ)表示數值頻散,其絕對值越小，數值頻散越小。

有限差分方法在一個空間網格上傳播的時間誤差是[18]:

(18)

表1比較了時間-空間域最小二乘方法, 傳統模板線性方法和新模板線性方法在不同速度時頻散誤差絕對值小于10-5所需的算子長度。由表1可見，假如時間步長τ=1ms，當地震波速分別為1 500、2 000 和2 500m/s時，最小二乘方法和傳統模板線性方法對應的算子長度M 分別為7、4 和3。由表1的最后一行可見，當時間步長從1ms增加到2.5ms時，新模板線性方法所需的算子長度并沒有改變(僅僅增加了d)，由此數值模擬需要的機時可以減少一半多。表2顯示了表1最后一行所對應的系數。

表 1 不同方法30 Hz頻率范圍內誤差小于10-5需要的算子長度對比

使用傳統模板和新模板的線性方法數值頻散如圖3所示。由此可見，使用新有限差分模板后，時間步長增加到τ=2.5 ms時仍然能夠在30 Hz頻率范圍內使得誤差小于10-5，而使用傳統的有限差分模板做不到這一點，由此證明了新有限差分模板和線性方法的有效性。

表 2 時間步長τ=2.5 ms時新有限差分模板對應的有限差分系數(d=r2/4)

3 數值模擬

3.1 均勻介質模型

首先考慮均勻介質模型，對于所有的有限差分算子υ=2 000 m/s,h=20 m,τ=2.5 ms,M=7。模型大小是700 m×700 m。震源在模型的中心位置，是高斯函數的二階導數[19]，

(18)

其中：f0=55 Hz；f0/π是主頻頻率;t0=4/f0;c是常數。

圖4顯示了5 000 ms時波場不同有限差分算法所對應快照的切片振幅圖(z/dz=120)，使用時間步長1 ms的傳統模板線性方法作為解析解。從圖4可見，使用新模板線性方法后時間步長增加到2.5 ms，數值頻散不明顯；而使用傳統模板線性方法時間步長增加到2.5 ms后，數值頻散明顯。

圖3 傳統模板與新模板線性方法的頻散誤差(h=20 m,τ=2.5 ms)Fig.3 Dispersion error of the former and new linear method (h=20 m，τ=2.5 ms)

4 鹽丘模型

圖5顯示了美國勘探地球物理學會的鹽丘模型，其速度從1 500 m/s變化到4 800 m/s[20]。震源位置用紅色菱形表示。震源函數使用式(18)，其中震源頻率等于45 Hz。對于所有的有限差分算子空間步長h=20 m。采用變算子長度方法(Liu and Sen, 2011)，當波速小于2 000 m/s時，算子長度M=7；當波速大于等于2 000 m/s小于2 500 m/s時，算子長度M=4；當波速大于等于2 500 m/s時，算子長度M=3。圖6顯示了不同有限差分算子在3 000 ms時的波場快照切片振幅圖 (z=400 m) ，其中時間-空間域泰勒展開方法時間步長τ=1 ms, 時間-空間域傳統模板線性方法時間步長τ=1 ms，新模板線性方法時間步長τ=2.5 ms。從圖6紅色橢圓中的部分對比可以看到，新模板線性方法在時間步長τ=2.5 ms時的空間頻散要小于時間-空間域泰勒展開方法在時間步長τ=1 ms時的空間頻散，而兩者的時間頻散相似(在波的傳播過程中，空間頻散使得高頻的相速度減小)。因此使用新模板線性方法既保持了精度又提高了計算效率。

圖4 不同有限差分算子對應的波場快照切片圖Fig.4 Slices of wave field snapshots using different methods

圖5 BP鹽丘速度模型Fig.5 BP salt velocity model

圖6波場快照的切片振幅對比圖。從上到下依次是時間-空間域泰勒展開法；時間-空間域傳統模板線性方法；時間-空間域新模板線性方法。

圖6 波場快照的切片振幅對比圖Fig.6 Comparision chart of slices of wave field snapshots

5 結論

基于前人提出的新有限差分模板，提出使用線性方法確定新有限差分模板對應的算子系數。其主要特點有二：其一是改進了線性方法，使之可以應用于新有限差分模板。使用新模板改進線性方法之后，在其他條件相同的情況下，時間步長可以從1 ms增加到2.5 ms，使得數值模擬的效率提高了一倍多；其二，線性方法計算新模板有限差分系數的效率要遠高于優化方法計算新模板有限差分系數的效率。鑒于新模板線性有限差分方法能夠有效去除頻散，在保持相同精度的情況下提高了計算效率，因此可以代替傳統的有限差分系數應用于聲波數值模中。

References)

[1] 連尉平,盧大偉,唐方頭,等.汶川MS8.0地震逆沖滑動沿斷層深度的分布特征及其數值模擬解釋[J].地震工程學報,2014,36(1):88-92.

LIAN Wei-ping,LU Da-Wei,TANG Fang-tou,et al.Distribution Characteristics of the Thrust Slip along Fault Depth in the WenchuanMS8.0 Earthquake and the Interpretation Based on Numerical Simulation[J].China Earthquake Engineering Journal,2014,36(1):88-92.(in Chinese)

[2] 董非非,溫燕林,鄧輝,等.地震海嘯在東海大陸架上傳播過程的數值模擬研究[J].地震工程學報,2013,35(1):154-159.

Dong Fei-fei,WEN Yan-lin,DENG Hui,et al.Numerical Simulation of Seismic Tsunami Dissemination at the Continental Shelf of the East China Sea[J].China Earthquake Engineering Journal,2013,35(1):154-159.(in Chinese)

[3] Marfurt K J.Accuracy of Finite-difference and Finite-element Modeling of the Scalar and Elastic Wave Equations[J].Geophysics,1984,49(5):533-549.

[4] Dablain M A.The Application of High-order Differencing to the Scalar Wave Equation[J].Geophysics,1986,51(1):54-66.

[5] Liu Y.Globally Optimal Finite-difference Schemes Based on Least Squares[J].Geophysics,2013,78(4):T113-T132.

[6] Chu C,Stoffa P L.Determination of Finite-difference Weights Using Scaled Binomial Windows[J].Geophysics,2012,77(3):17-26.

[7] Zhang J H,Yao Z X.Optimized Finite-difference Operator for Broadband Seismic Wave Modeling[J].Geophysics,2012,78(1):A13-A18.

[8] 王德利,雍運動,韓立國,等.三維黏彈介質地震波場有限差分并行模擬[J].西北地震學報,2007,29(1):30-34.

WANG De-li,YONG Yun-dong,HAN Li-guo,et al.Parallel Simulation of Finite Difference for Seismic Wavefield Modeling in 3-D Viscoelastic Media[J].Northwestern Seismological Journal,2007,29(1):30-34.(in Chinese)

[9] 嚴武建,王彥賓,石玉成.基于偽譜和有限差分混合方法的蘭州盆地強地面運動二維數值模擬[J].西北地震學報,2013,35(2):302-310.

YAN Wu-jian,WANG Yan-bin,SHI Yu-cheng.2D Simulation of the Strong Ground Motion in Lanzhou Basin with Hybrid PSM/FDM Method[J].China Earthquake Engineering Journal,2013,35(2):302-310.(in Chinese)

[10] 黃翼堅.非均勻TI介質P-SV波傳播交錯網格高階有限差分數值模擬[J].西北地震學報,2008,30(1):11-16.

HUANG Yi-jian.High-order Staggered-grid Finite Difference Numerical Modeling for P-SV Wave Propagation in Heterogeneous Transversely Isotropic Media[J].Northwestern Seismological Journal,2008,30(1):11-16.(in Chinese)

[11] Liu H W,Li B,Liu H,et al,The Issues of Prestack Reverse Time Migration and Solutions with Graphic Processing Unit Implementation[J].Geophysical Prospecting,2012,60(5):906-918.

[12] Yan H,Liu Y.Pre-stack Reverse-time Migration Based on the Time-space Domain Adaptive High-order Finite-difference Method in Acoustic VTI Medium[J].Journal of Geophysics and Engineering,2013,10(1):015010.

[13] Dablain M.The Application of High-order Differencing to the Scalar Wave Equation[J].Geophysics,1986,51(1):54-66.

[14] Liu Y,M K Sen.Acoustic VTI Modeling with a Time-space Domain Dispersion-relation-based Finite-difference Scheme[J].Geophysics,2010,75(3):7-13.

[15] Yang D,Lu M,Wu R S,et al.An Optimal Nearly-analytic Discrete Method for 2D Acoustic and Elastic Wave Equations[J].Bull Seismol Soc Am,2004,94(5):1982-1991.

[16] Liang W Q,Wang Y F,Yang C C.Comparison of Numerical Dispersion in Acoustic Finite-difference Algorithms[J].Exploration Geophysics,2014,46(2):206-212.

[17] Liu H,Dai N,Niu F,et al.An Explicit Time Evolution Method for Acoustic Wave Propagation[J].Geophysics, 2014,79(3):117-124.

[18] Liu Y,M K Sen.A New Time-space Domain High-order Finite-difference Method for the Acoustic Wave Equation[J].Journal of Computational Physics,2009,228:8779-8806.

[19] 王彥飛，I E 斯捷潘諾娃，V N 提塔連科，等.地球物理數值反演問題[M].北京：高等教育出版社，2011.

WANG Yan-fei,Stepanova I E,Titarenko V N, et al. Inverse Problems in Geophysicsand Solution Method[M].Beijing:Highter Education Press,2011.(in Chinse)

[20] Billette F J,S Brandsbery-Dahl. The 2004 BP Velocity Benchmark[C]//67thAnnual EAGE Meeting,EAGE, Expanded Abstracts,B305.

Acoustic Wave Equation Modeling with a New Time-Space Domain Finite Difference Stencil and an Improved Linear Algorithm

LIANG Wen-quan, LAI Yong-ming, JIANG Lin,WU Chao-fan

(CollegeofResourceEngineering,LongyanUniversity,Longyan364000,Fujian,China)

Numerical methods prove useful in exploration seismology. The most commonly used numerical methods are the finite difference, finite volume, and finite element methods. These methods constitute the basis for reverse time migration and full-waveform inversion. Finite difference methods are widely used in wave field extrapolation because of their higher computation efficiency, lower memory requirements, and easier implementation. During discretization of the time and the spatial derivatives in the wave equation, grid dispersion often occurs. Grid dispersion can result in artificial waves and inaccurate wave fields. Therefore, finding suitable finite difference operator coefficients to preserve the dispersion relationship of the wave equation, thus reducing grid dispersion, is one of the most important issues when using finite difference methods. To reduce grid dispersion, the traditional method uses high-order finite difference schemes in the spatial domain. However, waves are simultaneously propagated in time and space. Therefore, some researchers propose finite difference schemes based on the time-space domain dispersion relationship. Most commonly used are the high-order Taylor expansion method and the optimized method. However, the time step is relative small even in the optimized time-space domain method. Recently, a new finite difference stencil has been proposed to increase the time step while preserving the accuracy with the least-squares method. The time-space domain dispersion relationship of this new finite difference stencil is linear. Therefore, in this paper, we propose using this improved linear method, with the new finite difference stencil, to obtain the finite difference coefficients for the acoustic wave equation. We demonstrate, by dispersion analysis and numerical simulation, that with this new FD stencil and its improved linear solution, the wave equation simulation speed can be doubled compared with the previous linear method.

acoustic wave-equation modeling; time-space domain; finite difference scheme; linear method

2015-09-11

福建省科技廳青年基金(2016J05104)；福建省教育廳項目(JK2014050)；龍巖學院博士啟動基金(LB2014010)

梁文全(1983-)，男，博士，講師，主要從事地震波數值模擬與偏移研究。E-mail：onethink2002@foxmail.com。

P631.4

A

1000-0844(2016)05-0815-07

10.3969/j.issn.1000-0844.2016.05.0815