辛維， 閆子超, 梁文全， 陳雨紅， 楊長春

1 中國科學(xué)院地質(zhì)與地球物理所 中國科學(xué)院油氣資源研究重點實驗室， 北京 100029 2 中國科學(xué)院大學(xué)， 北京 100049

?

用于彈性波方程數(shù)值模擬的有限差分系數(shù)確定方法

辛維1,2， 閆子超1,2, 梁文全1， 陳雨紅1， 楊長春1

1 中國科學(xué)院地質(zhì)與地球物理所 中國科學(xué)院油氣資源研究重點實驗室， 北京 100029 2 中國科學(xué)院大學(xué)， 北京 100049

有限差分方法是波場數(shù)值模擬的一個重要方法，交錯網(wǎng)格差分格式比規(guī)則網(wǎng)格差分格式穩(wěn)定性更好，但方法本身都存在因網(wǎng)格化而形成的數(shù)值頻散效應(yīng)，這會降低波場模擬的精度與分辨率.為了緩解有限差分算子的數(shù)值頻散效應(yīng)，精確求解空間偏導(dǎo)數(shù)，本文把求解波動方程的線性化方法推廣到用于求解彈性波方程交錯網(wǎng)格有限差分系數(shù)；同時應(yīng)用最大最小準(zhǔn)則作為模擬退火(SA)優(yōu)化算法求解差分系數(shù)的數(shù)值頻散誤差判定標(biāo)準(zhǔn)來求解有限差分系數(shù).通過上述兩種方法，分別利用均勻各向同性介質(zhì)和復(fù)雜構(gòu)造模型進(jìn)行了數(shù)值正演模擬和數(shù)值頻散分析，并與傳統(tǒng)泰勒展開算法、最小二乘算法進(jìn)行比較，驗證了線性化方法和模擬退火方法都能有效壓制數(shù)值頻散，并比較了各個算法的特點.

彈性波方程； 數(shù)值模擬； 線性化方法； 模擬退火算法

1 引言

采用有限差分法對彈性波進(jìn)行數(shù)值模擬是一種常用方法(王秀明等，2004；王東等，2006；周竹生等，2007；高靜懷等，2012).交錯網(wǎng)格有限差分(FD)方法因為穩(wěn)定性好(Igel et al., 1992; Dong et al., 2000)、計算效率高、所需內(nèi)存小、實現(xiàn)簡單，而廣泛應(yīng)用于地震正演研究(Luo and Schuster, 1990; Graves, 1996; Saenger and Bohlen, 2004)，同時是地震逆時偏移成像技術(shù)得以快速發(fā)展的基礎(chǔ)(Kelly et al., 1976；Dablain, 1986；Levander, 1988；Liu et al., 2012；Yan and Liu, 2013).

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

為了緩解或者壓制網(wǎng)格數(shù)值頻散效應(yīng)，一般采用兩種措施：其一是采用低階差分格式，要求時間步長和空間間距非常小，這會造成所需內(nèi)存和計算量過大而無法應(yīng)用于三維的情況；其二是采用高階差分格式(裴正林，2005；李振春等，2007)，一般情況下，時間方向采用二階差分格式，空間方向采用高階差分格式.高階差分格式主要是利用等波紋準(zhǔn)則或最小誤差準(zhǔn)則優(yōu)化設(shè)計差分系數(shù)實現(xiàn)對空間偏導(dǎo)數(shù)的近似.偽譜法(朱多林和白超英，2012；唐小平等，2012)則是通過正、反傅里葉變換來實現(xiàn)空間偏導(dǎo)數(shù)的精確求解.兩種方法相比，高階空間差分格式因為具有精度適中、計算效率高的特征而得到廣泛應(yīng)用，特別是用于地震逆時偏移成像.Zhang和Yao(2013a，b)通過模擬退火方法確定聲波方程有限差分系數(shù)，并且給出了萬分之一的頻散誤差容限；Liu(2013)以及Yang等(2013)通過最小二乘方法確定波動方程的有限差分系數(shù).這些方法都減小了數(shù)值頻散，提高了數(shù)值模擬的計算效率.在利用有限差分法求解聲波方程時，可以使用線性方法確定交錯網(wǎng)格有限差分系數(shù)以得到較好模擬效果(Liang et al., 2013; Liang et al., 2014a,b).本文主要把聲波方程的線性化方法推廣到彈性波方程用于求解有限差分系數(shù)，同時提出在既定波數(shù)域內(nèi)以最大最小準(zhǔn)則為判定依據(jù)并用模擬退火算法對彈性波方程交錯網(wǎng)格有限差分系數(shù)進(jìn)行求解.然后將這兩種方法與泰勒展開法、最小二乘方法進(jìn)行了頻散分析比較和一階速度應(yīng)力交錯網(wǎng)格彈性波對比數(shù)值模擬，數(shù)值模擬對比表明三種方法都較泰勒展開法能有效地壓制數(shù)值頻散.同時試驗表明：利用線性化方法直接求解有限差分系數(shù)，其計算效率要遠(yuǎn)遠(yuǎn)高于最小二乘方法和模擬退火算法.

2 空間域交錯網(wǎng)格有限差分系數(shù)的計算方法

2.1 基本公式及已有算法

非均質(zhì)的速度-應(yīng)力彈性波動方程可以為(Virieux, 1986)

(1)

其中，t為時間，x和z為空間變量，(vx,vz)是速度向量，(τxx,τzz,τxz)是應(yīng)力向量，λ(x,z)和μ(x,z)是Lamé系數(shù)，b=b(x,z)是密度的倒數(shù).求解方程(1)，通常時間方向采用二階差分格式，空間方向采用高階差分格式.

2.1.1 泰勒展開法

一般來說，空間導(dǎo)數(shù)的精度通過使用高階有限差分系數(shù)來提高，一階導(dǎo)數(shù)的2M階有限差分交錯網(wǎng)格定義為

-u(x-mh+0.5h)]，

(2)

其中am是有限差分系數(shù)，h是空間網(wǎng)格間距.把平面波u(x)=u0eikx帶入上面公式(2)，其中u0為常數(shù)，i為虛數(shù)單位，k為波數(shù)，并令β=kh/2得到(Yangetal., 2013)：

(3)

公式(3)兩邊對波數(shù)k進(jìn)行泰勒展開，然后匹配k的系數(shù)，就可以得到傳統(tǒng)的交錯網(wǎng)格有限差分的系數(shù)，但是這種方法限定了波數(shù)(頻率)范圍，如果需要覆蓋較廣的波數(shù)域需要增加算子長度增加計算量.

2.1.2 最小二乘算法

計算有限差分系數(shù)的目的是用差分方程逼近微分方程，在更廣的頻率域范圍之內(nèi)使得公式(3)得以滿足.為此，建立目標(biāo)函數(shù)為 (Yang et al.,2013)

(4)

其中I是有限差分系數(shù)的非線性函數(shù)，并使I→min，fm(β)=sin((2m-1)β)，b為積分上限且小于等于π/2.目標(biāo)函數(shù)的極值點在偏導(dǎo)數(shù)為零處，由此可以確定方程得到最小二乘交錯網(wǎng)格有限差分系數(shù).但是這種方法不能精確的控制誤差在既定波數(shù)域內(nèi)的分布，此外，由于需要計算積分，此方法求得差分算子耗時較多.

2.2 線性化方法

從公式(3)出發(fā)，可以導(dǎo)出求取彈性波方程交錯網(wǎng)格有限差分系數(shù)的線性化方程.假設(shè)有M個均勻分布的波數(shù)點k(i)滿足公式(3)所表示的頻散關(guān)系，則由(3)立即得到線性方程組Ax=b，其中x為有限差分系數(shù)，系數(shù)矩陣A和右端項b分別為

(5)

(6)

其中β(i)(i=1,2,…,M)均勻分布于rπ/2M與rπ/2之間，其中r是用來確定波數(shù)范圍的比值，使得在特定的波數(shù)范圍內(nèi)頻散誤差小于萬分之一.本文定義r為

(7)

其中ku為波數(shù)的上限，f為最大頻率，v為波速.解線性方程組Ax=b，就可以得到彈性波方程交錯網(wǎng)格有限差分系數(shù).

2.3 模擬退火算法

求解有限差分系數(shù)是函數(shù)逼近問題，這里可以選擇最大最小準(zhǔn)則為逼近準(zhǔn)則：通過優(yōu)化系數(shù)使得在同樣波數(shù)覆蓋范圍內(nèi)，數(shù)值頻散帶來的最大誤差E最小，即可以通過方程(3)建立目標(biāo)函數(shù)(5).在此準(zhǔn)則可以使得誤差可控性增強(qiáng)，在既定的波數(shù)域范圍內(nèi)最大頻散誤差最小，實現(xiàn)最大風(fēng)險最小化，使整個波數(shù)域范圍內(nèi)的頻散誤差都在可以接收的范圍內(nèi)，其他準(zhǔn)則達(dá)不到同樣效果.該最大最小準(zhǔn)則下的問題為非線性問題，可以利用全局尋優(yōu)算法在一定準(zhǔn)則下對系數(shù)進(jìn)行遍歷求解：

(8)

并使得E→min.

模擬退火算法的程序框圖如圖1所示.

3 頻散分析

為了對上述各種方法進(jìn)行比較，本節(jié)首先分別利用泰勒展開算法(TE)、模擬退火算法(SA)、最小二乘算法(LS)和線性化方法(Linearmethod(LM))求解彈性波方程的有限差分系數(shù)，然后對其進(jìn)行頻散分析.所求解的有限差分算子長度在M=4、6、8時的差分系數(shù)如表1所示，其中最小二乘系數(shù)來自文獻(xiàn)(Yang et al.,2014).

表1 M=4、6、8時不同方法求得的有限差分系數(shù)Table 1 Finite difference coefficients in various methods with M=4, 6, 8

圖1 模擬退火算法程序流程Fig.1 Flowchart of simulated annealing algorithm

對于有限差分算子，可以使用公式(9)衡量頻散誤差，表達(dá)式為

(9)

其中E(β)越接近1，誤差越小，數(shù)值頻散對比見圖2，從圖中可得到：

(1)隨著算子長度的增加，各種方法都能夠在更大的波數(shù)范圍內(nèi)保持頻散關(guān)系.

(2)在一定的頻散誤差范圍內(nèi)，泰勒展開方法覆蓋的波數(shù)范圍最小，模擬退火和最小二乘方法覆蓋的波數(shù)范圍最大，線性化方法覆蓋的波數(shù)范圍遠(yuǎn)大于泰勒展開所覆蓋的波數(shù)范圍而接近于優(yōu)化方法覆蓋的波數(shù)范圍.

(3)雖然線性化方法覆蓋的波數(shù)范圍略小于優(yōu)化方法，但是該方法在低波數(shù)(低頻)時的頻散誤差小于優(yōu)化方法.同時考慮到震源的頻譜，線性化方法得到的有限差分系數(shù)能夠有效減少數(shù)值頻散.

(4)利用線性化方法直接求解有限差分系數(shù)，其計算效率要遠(yuǎn)遠(yuǎn)高于最小二乘方法和模擬退火算法.以算子長度M=8為例，線性化方法的CPU計算時間在1秒鐘以內(nèi)，而模擬退火算法則需要2個小時左右.當(dāng)算子長度長時，線性化方法的計算優(yōu)勢更為明顯.

4 數(shù)值模擬

在本次正演模擬實驗中，分別選取了簡單的均勻介質(zhì)模型和較為復(fù)雜的Marmousi模型進(jìn)行對比計算.

4.1 均勻介質(zhì)模型

在簡單的均勻介質(zhì)模型中，設(shè)定縱波速度為2000 m·s-1， 橫波速度為1000 m·s-1，密度為1800 kg·m-3.網(wǎng)格大小為350×350，網(wǎng)格間距h=10 m，有限差分系數(shù)使用與圖1b相同的有限差分系數(shù)，時間步長為1 ms.使用主頻為30 Hz的雷克子波，震源放在模型的中心位置.震源加載方式不同會影響波場生成，本文震源加在正應(yīng)力上，因而僅激發(fā)縱波.

不同方法得到的彈性波有限差分系數(shù)在800 ms時的波場快照如圖3所示，其中x為縱波水平分量，z為縱波垂直分量.由圖3 可以看出，泰勒展開方法有限差分系數(shù)的數(shù)值頻散最大，模擬退火方法、最小二乘方法和線性方法的有限差分系數(shù)數(shù)值頻散效應(yīng)都得到了有效的壓制.

為了進(jìn)一步比較頻散，抽取圖3(a—d)z=1900 m時的切片，如圖4所示.由圖4可以看到，泰勒展開方法的頻散最為嚴(yán)重，其他3種方法都有效減小了數(shù)值頻散.

4.2 Marmousi模型

圖 5顯示了Marmousi II速度和密度模型，其縱波速度從1500 m·s-3變化到4800 m·s-3，橫波速度從305 m·s-3變化到2800 m·s-3.震源函數(shù)使用主頻為30 Hz的雷克子波，震源在(10 m, 3000 m)處，檢波器記錄深度為10 m.模擬試驗中，使用了分裂的PML吸收邊界(Berenger, 1994).有限差分系數(shù)使用與圖2c相同的有限差分系數(shù)，空間步長為10 m，時間步長為1 ms.

圖6顯示了泰勒展開法(TE)、模擬退火方法(SA)、最小二乘方法(LS)和線性化方法(LM)計算有限差分系數(shù)的地震記錄.圖6a是泰勒展開法差分格式的地震記錄，數(shù)值頻散明顯；圖6(b—d)分別是用優(yōu)化方法和線性化方法有限差分系數(shù)得到的x方向的地震記錄.可以看出，數(shù)值頻散得到了明顯的壓制.

圖2 不同方法求得的有限差分系數(shù)的頻散誤差 (a)M=4；(b)M=6；(c)M=8.Fig.2 Dispersion errors of different FD coefficients determination methods

圖7抽取了x=3000 m處的地震記錄以進(jìn)一步比較地震記錄的數(shù)值頻散.其中泰勒展開方法有限差分系數(shù)的數(shù)值頻散最為明顯，其他3種方法的數(shù)值頻散都得到了壓制，同時從圖中可以看出線性方法和模擬退火方法有限差分系數(shù)的地震記錄幾乎重合.

5 結(jié)論

本文研究了確定彈性波方程有限差分系數(shù)的模擬退火算法和線性化方法.通過對泰勒展開方法、最小二乘方法、模擬退火方法和線性化方法交錯網(wǎng)格有限差分系數(shù)數(shù)值頻散的對比分析，可以得到以下幾點認(rèn)識：

(1)線性化方法、模擬退火方法和最小二乘等三種方法的數(shù)值頻散均小于泰勒展開方法.

(2)在低頻范圍，線性化方法的數(shù)值頻散要小于其他方法；同時注意到震源頻譜的主能量集中在中低頻，因此線性化方法能夠有效去除數(shù)值頻散，其數(shù)值效果與最小二乘方法以及模擬退火方法相當(dāng).

(3)利用線性化方法直接求解有限差分系數(shù)，其計算效率要遠(yuǎn)遠(yuǎn)高于最小二乘方法和模擬退火算法.

通過均勻介質(zhì)模型和復(fù)雜構(gòu)造模型的數(shù)值模擬及其結(jié)果的對比分析認(rèn)為，線性化方法、最小二乘方法和模擬退火方法都顯著地壓制了數(shù)值頻散(它們的地震記錄幾乎重合)，因而線性化方法確定交錯網(wǎng)格有限差分系數(shù)可以用于速度-應(yīng)力交錯網(wǎng)格彈性波數(shù)值模擬和逆時偏移成像，以提高計算的精度和效率.致謝 感謝審稿專家提出的寶貴修改意見和編輯部的大力支持！

Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves.JournalofComputationalPhysics, 114(2): 185-200.

Dablain M A. 1986. The application of high-order differencing to the scalar wave equation.Geophysics, 51(1): 54-66.

圖3 不同方法有限差分系數(shù)的波場快照圖 (a)泰勒展開方法(TE); (b)模擬退火算法(SA); (c)最小二乘方法(LS); (d)線性方法(LM)，左為x分量，右為z分量.Fig.3 Wave field snapshots of FD coefficients determined by different methods (a) Taylor method (TE); (b) Simulated annealing algorithm (SA); (c) Least squares method (LS); (d) The linear method (LM); left for the x component, right for the z component.

圖4 波場快照的切片圖(x分量)Fig.4 Slices of snapshots of the wave field (x component)

圖5 Marmousi速度模型 (a) P波；(b) S波；(c) 密度.Fig.5 Marmousi velocity model (a) P waves;(b) S waves;(c) Density.

圖6 不同有限差分系數(shù)算法模擬得到的Marmousi模型的地震記錄 (a)泰勒展開方法(TE);(b)模擬退火算法(SA);(c)最小二乘方法(LS); (d)線性方法(LM).Fig.6 Seismic records of the Marmousi model using different FD coefficient determination methods (a) Taylor method (TE); (b) Simulated annealing algorithm (SA); (c) Least squares method (LS); (d) The linear method (LM).

圖7 不同方法有限差分系數(shù)地震道的對比Fig.7 Comparison of seismic traces obtained using different FD coefficient determination methods

Dong L G, Ma Z T, Cao J Z. 2000. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation.ChineseJournalofGeophysics, 43(6): 856-864.

Gao J H, He Y Y, Ma Y C. 2012. Comparison of the Rayleigh wave in elastic and viscoelastic media.ChineseJournalGeophysics(in Chinese), 55(1): 207-218, doi: 10.6038/j.issn.0001-5733.2012.01.020. Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences.BulletinoftheSeismologicalSocietyofAmerica, 86(4): 1091-1106.

Igel H, Riollet B, Mora P. 1992. Accuracy of staggered 3-D finite-difference grids for anisotropic wave propagation. //62nd Annual International Meeting, SEG, Expanded Abstracts, 1244-1246. Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms: a finite-difference approach.Geophysics, 41(1): 2-27. Levander A R. 1988. Fourth-order finite-difference P-SV seismograms.Geophysics, 53(11): 1425-1436. Li Z C, Zhang H, Liu Q M, et al. 2007. Numeric simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm.OilGeophysicalProspecting(in Chinese), 42(5): 510-515. Liang W Q, Wang Y F, Yang C C, et al. 2013. Acoustic wave equation modeling with new time-space domain finite difference operators.ChineseJournalofGeophysics, 56(6): 840-850, doi: 10.1002/cjg2.20076.

Liang W Q, Wang Y F, Yang C C. 2014a. Comparison of numerical dispersion in acoustic finite-difference algorithms.ExplorationGeophysics, doi: 10.1071/EG13072.

Liang W Q, Wang Y F, Yang C C. 2014b. Determining finite difference weights for the acoustic wave equation by a new dispersion-relationship-preserving method.GeophysicalProspecting, 63(1): 11-22. Liu H W, Li B, Liu H, et al. 2012. The issues of prestack reverse time migration and solutions with Graphic Processing Unit implementation.GeophysicalProspecting, 60(5): 906-918.

Liu Y. 2013. Globally optimal finite-difference schemes based on least squares.Geophysics, 78(4): T113-T132.

Luo Y, Schuster G. 1990. Parsimonious staggered grid finite-differencing of the wave equation.GeophysicalResearchLetters, 17(2): 155-158.Pei Z L. 2005. Numerical simulation of elastic wave equation in 3-D isotropic media with staggered-grid high-order difference method.GeophysicalProspectingforPetroleum(in Chinese), 44(4): 308-315.

Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid.Geophysics, 69(2): 583-591.

Tang X P, Bai C Y, Liu K H. 2012. Elastic wavefield simulation using separated equations through pseudo-spectral method.OilGeophysicalProspecting(in Chinese), 47(1): 19-26.

Virieux J. 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method.Geophysics, 51(4): 889-901.

Wang D, Zhang H L, Wang X M. 2006. A numerical study of acoustic wave propagation in partially saturated poroelastic rock.ChineseJournalofGeophysics(in Chinese), 49(2): 524-532, doi: 10.3321/j.issn:0001-5733.2006.02.027.

Wang X M, Zhang H L, Wang D. 2003. Modelling of seismic wave propagation in heterogeneous poroelastic media using a high order staggered finite difference method.ChineseJournalofGeophysics(in Chinese), 46(6): 842-849.

Yan H Y, Liu Y. 2013. Acoustic VTI modeling and pre-stack reverse-time migration based on the time-space domain staggered-grid finite-difference method.JournalofAppliedGeophysics, 90: 41-52. Yang L, Yan H Y, Liu H. 2013. Least squares staggered-grid finite-difference for elastic wave modelling.ExplorationGeophysics, doi: 10.1071/EG13087. Zhang J H, Yao Z X. 2013a. Optimized finite-difference operator for

broadband seismic wave modeling.Geophysics, 78(1): A13-A18.

Zhang J H, Yao Z X. 2013b. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm.JournalofComputationalPhysics, 250: 511-526.

Zhou Z S, Liu X L, Xiong X Y. 2007. Finite-difference modelling of Rayleigh surface wave in elastic media.ChineseJournalofGeophysics(in Chinese), 50(2): 567-573, doi: 10.3321/j.issn:0001-5733.2007.02.030.

Zhu D L, Bai C Y. 2011. Review on the seismic wavefield forward modelling.ProgressinGeophysics(in Chinese), 26(5): 1588-1599, doi: 10.3969/j.issn.1004-2903.2011.05.011.

附中文參考文獻(xiàn)

高靜懷, 何洋洋, 馬逸塵. 2012. 黏彈性與彈性介質(zhì)中Rayleigh面波特性對比研究. 地球物理學(xué)報, 55(1): 207-218, doi: 10.6038/j.issn.0001-5733.2012.01.020.

李振春, 張華, 劉慶敏等. 2007. 彈性波交錯網(wǎng)格高階有限差分法波場分離數(shù)值模擬. 石油地球物理勘探, 42(5): 510-515.

裴正林. 2005. 三維各向同性介質(zhì)彈性波方程交錯網(wǎng)格高階有限差分法模擬. 石油物探, 44(4): 308-315.

唐小平, 白超英, 劉寬厚. 2012. 偽譜法分離波動方程彈性波模擬. 石油地球物理勘探, 47(1): 19-26.

王東, 張海瀾, 王秀明. 2006. 部分飽和孔隙巖石中聲波傳播數(shù)值研究. 地球物理學(xué)報, 49(2): 524-532, doi: 10.3321/j.issn:0001-5733.2006.02.027.

王秀明, 張海瀾, 王東. 2004. 利用高階交錯網(wǎng)格有限差分法模擬地震波在非均勻孔隙介質(zhì)中的傳播. 地球物理學(xué)報, 46(6): 842-849.

周竹生, 劉喜亮, 熊孝雨. 2007. 彈性介質(zhì)中瑞雷面波有限差分法正演模擬. 地球物理學(xué)報, 50(2): 567-573, doi: 10.3321/j.issn:0001-5733.2007.02.030.

朱多林, 白超英. 2011. 基于波動方程理論的地震波場數(shù)值模擬方法綜述. 地球物理學(xué)進(jìn)展, 26(5): 1588-1599, doi: 10.3969/j.issn.1004-2903.2011.05.011.

(本文編輯 張正峰)

Methods to determine the finite difference coefficients for elastic wave equation modeling

XIN Wei1,2, YAN Zi-Chao1,2, LIANG Wen-Quan1,CHEN Yu-Hong1, YANG Chang-Chun1

1KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China

Numerical simulation of the elastic wave equation with staggered-grid finite-difference algorithms is widely used to synthesize seismograms theoretically, and is also the basis of the reverse time migration. With some stability conditions, grid dispersion often appears because of the discretization of the time and the spatial derivatives in the wave equation. How to suppress the grid dispersion is therefore a key problem for finite-difference approaches. Different methods have been proposed to address this issue. The most commonly used methods are the high order Taylor expansion (TE) methods. In this paper, we extend the linear method for solving the acoustic wave equation to the staggered grid finite difference method for solving the elastic wave equation. We also apply the maximum-minimum criterion to measure the dispersion error when performing simulated annealing (SA) algorithm. Dispersion analysis and numerical simulation demonstrate that a linear method without iteration is nearly equal to the SA method and the least squares (LS) method in the space domain, and is better than the TE methods.For the finite difference coefficients obtained by the two methods, using homogeneous isotropic and complex structural model, we performed a numerical forward modeling and numerical dispersion analysis firstly, then compared it with the traditional Taylor expansion (TE) method and least squares(LS) method.Dispersion analysis and numerical simulation demonstrate the following conclusions: (1) With the increase of the length of the operator, various methods are able to maintain the dispersion relation in a larger wave number range. (2) The coefficients obtained by the TE method covers the minimal wave number range, coefficients from SA and LS method cover the maximal wave number range, the wave number range of linearization method is much larger than that of TE method, and is very close to that of the optimization method. (3)Although the wave number range of the linearization method is slightly less than the optimization method, but the dispersion error of this method is smaller in lower wave number (low frequency). Taking the spectrum of seismic source into account, the coefficient of linearization finite difference method is able to effectively reduce the numerical dispersion. (4) The linearization method can be used to solve finite difference coefficients directly. Its computational efficiency is much higher than the LS method and SA method.Comparison of the numerical simulation results from the uniform medium model and the complex structure model indicates that the linearization method, LS method and SA method can significantly suppress the numerical dispersion (seismograms almost coincide), thus the coefficients confirmed by the linearization finite difference method can be used to speed-stress staggered grid elastic wave modeling and time-reverse migration to improve the accuracy and efficiency of the calculation.

Elastic wave equation; Numerical simulation; Linear method; Simulated annealing

10.6038/cjg20150724.

國家自然科學(xué)基金項目(41325016和11271349)資助．

辛維，男，1986年生，中國科學(xué)院地質(zhì)與地球物理研究所博士畢業(yè)，主要從事地震波正演研究和地球物理儀器研發(fā)． E-mail：332630765@qq.com

10.6038/cjg20150724

P631

2014-05-14，2015-04-23收修定稿

辛維，閆子超, 梁文全等. 2015. 用于彈性波方程數(shù)值模擬的有限差分系數(shù)確定方法.地球物理學(xué)報，58(7):2486-2495,

Xin W, Yan Z C,Liang W Q, et al. 2015. Methods to determine the finite difference coefficients for elastic wave equation modeling.ChineseJ.Geophys. (in Chinese),58(7):2486-2495,doi:10.6038/cjg20150724.