適于聲波方程數(shù)值模擬的時(shí)間-空間域有限差分算子改進(jìn)線性算法

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

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

?

適于聲波方程數(shù)值模擬的時(shí)間-空間域有限差分算子改進(jìn)線性算法

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

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

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

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

0 引言

數(shù)值模擬在實(shí)踐中得到了廣泛的應(yīng)用[1-2]。其中有限差分方法因?yàn)橛?jì)算效率高、所需內(nèi)存小、實(shí)現(xiàn)簡單而廣泛應(yīng)用于地震正演研究[3-10],同時(shí)構(gòu)成了逆時(shí)偏移成像和全波形反演的基礎(chǔ)[11-12]。

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

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

1 有限差分算子確定方法

對于二維的聲波方程,

(1)

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

(2)

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

(3)

對于時(shí)間偏導(dǎo)數(shù)，采用二階有限差分近似，得到

(4)

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

(5)

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

(6)

令

(7)

根據(jù)震源頻率、速度和網(wǎng)格間距確定需要滿足頻散關(guān)系的波數(shù)范圍占總波數(shù)的比例[16]，

(8)

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

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

(9)

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

(10)

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

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

由方程(6)可以得到，

(11)

令式(11)左邊

(12)

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

(13)

整理得到

(14)

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

(15)

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

(16)

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

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

2 頻散分析

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

(17)

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

有限差分方法在一個(gè)空間網(wǎng)格上傳播的時(shí)間誤差是[18]:

(18)

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

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

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

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

3 數(shù)值模擬

3.1 均勻介質(zhì)模型

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

(18)

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

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

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

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

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

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

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

5 結(jié)論

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

References)

[1] 連尉平,盧大偉,唐方頭,等.汶川MS8.0地震逆沖滑動(dòng)沿?cái)鄬由疃鹊姆植继卣骷捌鋽?shù)值模擬解釋[J].地震工程學(xué)報(bào),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] 董非非,溫燕林,鄧輝,等.地震海嘯在東海大陸架上傳播過程的數(shù)值模擬研究[J].地震工程學(xué)報(bào),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] 王德利,雍運(yùn)動(dòng),韓立國,等.三維黏彈介質(zhì)地震波場有限差分并行模擬[J].西北地震學(xué)報(bào),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] 嚴(yán)武建,王彥賓,石玉成.基于偽譜和有限差分混合方法的蘭州盆地強(qiáng)地面運(yùn)動(dòng)二維數(shù)值模擬[J].西北地震學(xué)報(bào),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] 黃翼堅(jiān).非均勻TI介質(zhì)P-SV波傳播交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬[J].西北地震學(xué)報(bào),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 提塔連科，等.地球物理數(shù)值反演問題[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)；福建省教育廳項(xiàng)目(JK2014050)；龍巖學(xué)院博士啟動(dòng)基金(LB2014010)

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

P631.4

A

1000-0844(2016)05-0815-07

10.3969/j.issn.1000-0844.2016.05.0815