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

辛維， 閆子超, 梁文全， 陳雨紅， 楊長(zhǎng)春

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

?

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

辛維1,2， 閆子超1,2, 梁文全1， 陳雨紅1， 楊長(zhǎng)春1

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

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

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

1 引言

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

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

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

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

2.1 基本公式及已有算法

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

(1)

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

2.1.1 泰勒展開(kāi)法

一般來(lái)說(shuō)，空間導(dǎo)數(shù)的精度通過(guò)使用高階有限差分系數(shù)來(lái)提高，一階導(dǎo)數(shù)的2M階有限差分交錯(cuò)網(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)兩邊對(duì)波數(shù)k進(jìn)行泰勒展開(kāi)，然后匹配k的系數(shù)，就可以得到傳統(tǒng)的交錯(cuò)網(wǎng)格有限差分的系數(shù)，但是這種方法限定了波數(shù)(頻率)范圍，如果需要覆蓋較廣的波數(shù)域需要增加算子長(zhǎng)度增加計(jì)算量.

2.1.2 最小二乘算法

計(jì)算有限差分系數(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ù)的極值點(diǎn)在偏導(dǎo)數(shù)為零處，由此可以確定方程得到最小二乘交錯(cuò)網(wǎng)格有限差分系數(shù).但是這種方法不能精確的控制誤差在既定波數(shù)域內(nèi)的分布，此外，由于需要計(jì)算積分，此方法求得差分算子耗時(shí)較多.

2.2 線性化方法

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

(5)

(6)

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

(7)

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

2.3 模擬退火算法

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

(8)

并使得E→min.

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

3 頻散分析

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

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

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

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

(9)

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

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

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

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

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

4 數(shù)值模擬

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

4.1 均勻介質(zhì)模型

在簡(jiǎn)單的均勻介質(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ù)，時(shí)間步長(zhǎng)為1 ms.使用主頻為30 Hz的雷克子波，震源放在模型的中心位置.震源加載方式不同會(huì)影響波場(chǎng)生成，本文震源加在正應(yīng)力上，因而僅激發(fā)縱波.

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

為了進(jìn)一步比較頻散，抽取圖3(a—d)z=1900 m時(shí)的切片，如圖4所示.由圖4可以看到，泰勒展開(kāi)方法的頻散最為嚴(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.模擬試驗(yàn)中，使用了分裂的PML吸收邊界(Berenger, 1994).有限差分系數(shù)使用與圖2c相同的有限差分系數(shù)，空間步長(zhǎng)為10 m，時(shí)間步長(zhǎng)為1 ms.

圖6顯示了泰勒展開(kāi)法(TE)、模擬退火方法(SA)、最小二乘方法(LS)和線性化方法(LM)計(jì)算有限差分系數(shù)的地震記錄.圖6a是泰勒展開(kāi)法差分格式的地震記錄，數(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ù)值頻散.其中泰勒展開(kāi)方法有限差分系數(shù)的數(shù)值頻散最為明顯，其他3種方法的數(shù)值頻散都得到了壓制，同時(shí)從圖中可以看出線性方法和模擬退火方法有限差分系數(shù)的地震記錄幾乎重合.

5 結(jié)論

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

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

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

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

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

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ù)的波場(chǎng)快照?qǐng)D (a)泰勒展開(kāi)方法(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 波場(chǎng)快照的切片圖(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)泰勒展開(kāi)方法(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ù)地震道的對(duì)比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面波特性對(duì)比研究. 地球物理學(xué)報(bào), 55(1): 207-218, doi: 10.6038/j.issn.0001-5733.2012.01.020.

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

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

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

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

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

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

朱多林, 白超英. 2011. 基于波動(dòng)方程理論的地震波場(chǎng)數(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.

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

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

10.6038/cjg20150724

P631

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

辛維，閆子超, 梁文全等. 2015. 用于彈性波方程數(shù)值模擬的有限差分系數(shù)確定方法.地球物理學(xué)報(bào)，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.