一種全局優(yōu)化的隱式交錯(cuò)網(wǎng)格有限差分算法及其在彈性波數(shù)值模擬中的應(yīng)用

王洋， 劉洪， 張衡， 王之洋， 唐祥德

1 中國科學(xué)院地質(zhì)與地球物理研究所 中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 中國科學(xué)院大學(xué)， 北京 100049 3 國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室 廣州海洋地質(zhì)調(diào)查局， 廣州 510075

?

一種全局優(yōu)化的隱式交錯(cuò)網(wǎng)格有限差分算法及其在彈性波數(shù)值模擬中的應(yīng)用

王洋1, 2， 劉洪1*， 張衡3， 王之洋1, 2， 唐祥德1, 2

1 中國科學(xué)院地質(zhì)與地球物理研究所 中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 中國科學(xué)院大學(xué)， 北京 100049 3 國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室 廣州海洋地質(zhì)調(diào)查局， 廣州 510075

在數(shù)值模擬中，隱式有限差分具有較高的精度和穩(wěn)定性.然而，傳統(tǒng)隱式有限差分算法大多由于需要求解大型矩陣方程而存在計(jì)算效率偏低的局限性.本文針對(duì)一階速度-應(yīng)力彈性波方程，構(gòu)建了一種優(yōu)化隱式交錯(cuò)網(wǎng)格有限差分格式，然后將改進(jìn)格式由時(shí)間-空間域轉(zhuǎn)換為時(shí)間-波數(shù)域，利用二范數(shù)原理建立目標(biāo)函數(shù)，再利用模擬退火法求取優(yōu)化系數(shù).通過對(duì)均勻模型以及復(fù)雜介質(zhì)模型進(jìn)行一階速度-應(yīng)力彈性波方程數(shù)值模擬所得單炮記錄、波場快照分析表明：這種優(yōu)化隱式交錯(cuò)網(wǎng)格差分算法與傳統(tǒng)的幾種顯式和隱式交錯(cuò)網(wǎng)格有限差分算法相比不但降低了計(jì)算量，而且能有效的壓制網(wǎng)格頻散，使彈性波數(shù)值模擬的精度得到有效的提高.

隱式有限差分； 數(shù)值模擬； 交錯(cuò)網(wǎng)格； 彈性波方程 ； 模擬退火

1 引言

有限差分算法憑借其高計(jì)算效率以及易于編程實(shí)現(xiàn)等優(yōu)點(diǎn)，在地震波正演(Virieux, 1984, 1986; 裴正林，2005；嚴(yán)紅勇和劉洋，2012)、逆時(shí)偏移(Baysal et al., 1983; McMechan, 1983; Yan and Liu, 2013)和全波形反演(Virieux and Operto, 2009;Warner et al., 2013)等領(lǐng)域已經(jīng)得到了廣泛使用.尤其是交錯(cuò)網(wǎng)格有限差分，因具有更高的精度和更好的穩(wěn)定性而廣受歡迎.目前，通常使用的有限差分算法主要分為顯式和隱式兩種.傳統(tǒng)的波動(dòng)方程顯式有限差分系數(shù)通常由泰勒級(jí)數(shù)展開或優(yōu)化算法而得到(劉洋等，1998；董良國等，2000； Liu et al., 2014)，為了得到高精度的數(shù)值解，通常需要采用較長的差分算子.Kosloff 等(2010)認(rèn)為高階顯式差分算子存在“飽和效應(yīng)”，并且指出盡管高階顯式差分算法提高了精度，但精度提高的增量隨差分算子階數(shù)的增加而降低，經(jīng)驗(yàn)表明最經(jīng)濟(jì)的差分階數(shù)存在于四階到八階之間.相比之下隱式差分一般能以較短的差分算子長度得到較高的數(shù)值解精度，所以隱式差分備受關(guān)注.

與顯式差分不同，隱式差分算子是由有理展開得到的，它不但比級(jí)數(shù)展開收斂更快(Fornberg, 1998)，而且減小了高階差分算子的長度.盡管隱式差分算子已經(jīng)在許多領(lǐng)域中廣泛應(yīng)用(Lele, 1992; Ekaterinaris, 1999; Shang, 1999)， 但近年來才被引入解決地震波傳播問題(Liu and Sen, 2009a, 2009b; Kosloff et al., 2010; Chu and Stoffa, 2010，2012). Kosloff 等(2008)將隱式差分算法引入到針對(duì)常密度時(shí)間域聲波方程空間二階導(dǎo)數(shù)的求取中.并隨后在2010年，在聲波方程空間隱式二階導(dǎo)數(shù)的基礎(chǔ)上，推導(dǎo)出隱式一階導(dǎo)數(shù)交錯(cuò)網(wǎng)格有限差分算子，并將正演結(jié)果應(yīng)用于逆時(shí)偏移中，證明了該算子能夠有效的減少數(shù)值頻散.Chu 和Stoffa(2010)提出了頻率域二階聲波方程空間導(dǎo)數(shù)的隱式算法，之后在2012年給出了由二階隱式分裂時(shí)間積分和二階空間導(dǎo)數(shù)隱式差分相結(jié)合的聲波方程數(shù)值模擬方法.前面這些隱式算法，為了滿足高精度模擬的需求通常需要通過LU分解的方法完成五或更多對(duì)角矩陣求逆，極大的降低了算法的計(jì)算效率.Zhang 等 (2007)將空間導(dǎo)數(shù)算子分為x、y、z三個(gè)方向，推導(dǎo)了空間二階導(dǎo)數(shù)的分裂隱式差分格式，這種方法能夠給出精確的點(diǎn)源響應(yīng)并且能夠應(yīng)用于速度橫向變化的情況.在適當(dāng)參數(shù)的選擇下，這種算法是無條件穩(wěn)定的.Zhou 和 Zhang(2011)提出了預(yù)分解的策略，將隱式算法中原始的對(duì)角矩陣分解成兩個(gè)獨(dú)立的上三角矩陣和下三角矩陣，也就是使用L+U的預(yù)分解方法，替代了傳統(tǒng)的LU分解方法.并利用Taylor級(jí)數(shù)展開與最小二乘相結(jié)合的思想求取隱式算法優(yōu)化系數(shù)，極大的提高了模擬精度.然而，前述兩種方法只適用于二階聲波方程，無法將其推廣應(yīng)用于一階彈性波波動(dòng)方程求解中.Liu和Sen(2009a, b)在Claerbout (1985)基礎(chǔ)上，先后提出了針對(duì)一階、二階空間導(dǎo)數(shù)的規(guī)則網(wǎng)格和交錯(cuò)網(wǎng)格隱式有限差分新算法，并分別實(shí)現(xiàn)了該隱式算法來求解聲波方程和一階彈性波方程，從而進(jìn)行了數(shù)值模擬分析；之后，Liu (2014) 還在此隱式差分格式的基礎(chǔ)上，提出利用最小二乘方法來優(yōu)化該差分格式的差分系數(shù)，較大的提高了差分算子的計(jì)算精度和效率.

本文受Claerbout (1985)和Liu 和 Sen(2009b)思想的啟發(fā)，構(gòu)建出一種新的空間隱式差分格式，并采用非線性優(yōu)化法來計(jì)算交錯(cuò)網(wǎng)格有限差分的空間差分系數(shù).在差分系數(shù)優(yōu)化過程中，先使用二范數(shù)建立目標(biāo)函數(shù)，然后再用模擬退火法來求解優(yōu)化的系數(shù).最后，用此新隱式差分法來求解一階速度-應(yīng)力彈性波方程，并通過均勻介質(zhì)和復(fù)雜模型測試來驗(yàn)證新差分法的計(jì)算精度和計(jì)算效率，同時(shí)與傳統(tǒng)的顯式交錯(cuò)網(wǎng)格有限差分算法(Explicit Staggered-grid Finite-Difference Method, 簡稱ESFDM)和隱式交錯(cuò)網(wǎng)格有限差分算法(Implicit Staggered-grid Finite-Difference Method, 簡稱ISFDM)進(jìn)行了對(duì)比.

2 一階彈性波速度—應(yīng)力方程及高階交錯(cuò)網(wǎng)格有限差分

2.1 一階空間導(dǎo)數(shù)交錯(cuò)網(wǎng)格隱式有限差分格式

由平面波理論，令

(2)

(3)

將式(3)中的余弦與正弦項(xiàng)進(jìn)行泰勒級(jí)數(shù)展開，由左右系數(shù)對(duì)應(yīng)相等可以通過矩陣求逆的方法直接求取系數(shù)，見表1.

表1 傳統(tǒng)交錯(cuò)網(wǎng)格隱式有限差分系數(shù)Table 1 Coefficients of the conventional implicit staggered-grid finite-difference formulas

2.2 一階空間導(dǎo)數(shù)交錯(cuò)網(wǎng)格優(yōu)化隱式有限差分格式

(4)

即：

(5)

其中矩陣Bx，Ax分別由公式(1)中的am，bn，m=1,…M,n=1,2,…,N確定，矩陣大小為(NX-1)×(NX-1)或(NX-2)×(NX-2)，視交錯(cuò)網(wǎng)格中參數(shù)的定義情況而定，NX為橫向采樣點(diǎn)個(gè)數(shù).

為了達(dá)到高精度模擬的需求，公式(1)通常需要滿足N>1的條件.當(dāng)N=1時(shí)，有

(6)

此時(shí)Bx為三對(duì)角矩陣，可以通過追趕法或LU分解法進(jìn)行矩陣求逆.當(dāng)N=2時(shí)，有

(7)

此時(shí)Bx為五對(duì)角矩陣，追趕法已經(jīng)不適用，必須通過LU分解法進(jìn)行矩陣求逆.

由以上分析看出，當(dāng)N≥2時(shí)，求解每一個(gè)空間導(dǎo)數(shù)的過程中進(jìn)行五或更多對(duì)角矩陣求逆，此時(shí)已經(jīng)不能應(yīng)用追趕法快速求解，LU分解成為此時(shí)的最優(yōu)選擇.

表2 不同矩陣求逆方法平均CPU計(jì)算時(shí)間Table 2 Average CPU time of different algorithms for matrix inversion

隱式差分正演中，在求解每一個(gè)空間導(dǎo)數(shù)的過程中均需進(jìn)行一次對(duì)角矩陣求逆，因此，矩陣求逆的計(jì)算效率成為了整個(gè)算法效率至關(guān)重要影響因素.通過對(duì)追趕法和LU分解法計(jì)算時(shí)間對(duì)比測試可以看出，追趕法的求解速度明顯優(yōu)于LU分解.因此構(gòu)造一種新的可以由追趕法計(jì)算的隱式差分格式成為了提高隱式差分計(jì)算效率的有效途徑.

根據(jù)Claerbout(1985)以及Liu和Sen(2009b)的觀點(diǎn)，函數(shù)u(x)的二階有限差分算子可以通過增加一個(gè)高階項(xiàng)的方法來提高FD算子的精度，通常為了易于編程實(shí)現(xiàn)，它是以下面的形式出現(xiàn)的(Claerbout，1985) ，其中b為一個(gè)常量，表達(dá)式為

(8)

Liu 和 Sen (2009b)在此基礎(chǔ)上提出了交錯(cuò)網(wǎng)格隱式差分格式，通過增加三階導(dǎo)數(shù)項(xiàng)來提高精度，公式為

(9)

受他們的思想啟發(fā)，本文構(gòu)建出一種一階導(dǎo)數(shù)的優(yōu)化隱式交錯(cuò)網(wǎng)格有限差分算法(OptimizedImplicitStaggered-gridFinite-DifferenceMethod, 簡稱OISFDM)，表達(dá)式為

(10)

即：

(11)

其中a,c為可調(diào)常量，h為橫向網(wǎng)格間距.簡單起見，令

(12)

(13)

當(dāng)需要高階精度時(shí)可以將(12)寫成2M階差分形式：

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

(14)

將公式(10)寫成矩陣形式可以表示為

(15)

即：

(16)

其中

(17a)

(17b)

Cx和Ax矩陣大小均為(NX-1)×(NX-1)或(NX-2)×(NX-2)，視交錯(cuò)網(wǎng)格中參數(shù)的定義情況而定.同理，在z向時(shí)矩陣大小(NZ-1)×(NZ-1)或(NZ-2)×(NZ-2).NX和NZ分別表示數(shù)值模擬中橫向和縱向的網(wǎng)格點(diǎn)個(gè)數(shù).可以看出，公式(16)通過兩次三對(duì)角矩陣求逆來得到優(yōu)化ISFD算子，避免使用LU分解，提高了計(jì)算效率.

2.3 隱式算法中三對(duì)角矩陣的計(jì)算策略

由于在公式(16)中每求解一個(gè)空間導(dǎo)數(shù)都要進(jìn)行兩次三對(duì)角矩陣求逆.因此可以將一些重復(fù)的求逆過程提前計(jì)算出來，提高隱式算法的計(jì)算效率(Liu and Sen, 2009 a, b).通過對(duì)(17a)和(17b)三對(duì)角矩陣的觀察可以發(fā)現(xiàn)這兩個(gè)矩陣都只含有一個(gè)參數(shù)c或a.在每條對(duì)角線上，元素的值都是相等的，并且上下相對(duì)的兩條對(duì)角線的值也是相等的.在這種情況下，可以不需要真正的建立矩陣Cx和Ax，而是在追趕算法中將它們對(duì)應(yīng)的參數(shù)直接賦值，并且將追趕法求逆矩陣分解為LU兩個(gè)矩陣的過程在時(shí)間迭代前計(jì)算出來，這樣在節(jié)省內(nèi)存空間的同時(shí)降低了計(jì)算量.

3 模擬退火法求取交錯(cuò)網(wǎng)格優(yōu)化隱式有限差分系數(shù)

將u(x)=u0eikx，β=kh/2，(13)和(14)帶入公式(10)可得：

(18)

由此，本文中采用二范數(shù)原理建立2M階所對(duì)應(yīng)的目標(biāo)函數(shù)為

×ξ(β)dβ,

(19)

其中

(20)本文選擇目標(biāo)函數(shù)積分上下限分別為βl=0，βr=0.58π.

通過公式(11)可以看出，由于算子中存在系數(shù)相乘項(xiàng)，為非線性FD算子，不能通過常規(guī)的線性矩陣求逆的方法得到差分系數(shù).常用的非線性全局優(yōu)化方法有共軛梯度法、高斯牛頓法、模擬退火法等，模擬退火法因其能克服優(yōu)化過程陷入局部極小、克服初值依賴性、目標(biāo)函數(shù)靈活等優(yōu)點(diǎn)被廣泛認(rèn)可.因此，本文對(duì)公式(19)選用模擬退火法求解優(yōu)化系數(shù).

模擬退火算法最早的思想由Metropolis等(1953)提出，圖1展示了模擬退火法求解目標(biāo)函數(shù)的流程圖，該算法試圖在給定容差限Tolerance下尋找一組最佳的優(yōu)化系數(shù).圖中函數(shù)obj由公式(19)(20)定義，算法初值設(shè)定:MARKOV鏈長度N=1000, 初始溫度T=1000，終止溫度T0=100， 步長因子Stepfactor=0.02, 容差Tolerance=5×10-8,最大搜索值XMAX=[1；1], PreX=XMAX×rand, PreBestX=PreX, BestX=-XMAX×rand. 隨機(jī)選取下一點(diǎn)NextX=PreX+ Stepfactor ×XMAX×(rand-0.5)．

根據(jù)熱力學(xué)的原理，在溫度為T時(shí)，出現(xiàn)能量差為dE的降溫的概率為P(dE)，表示為

(21)

其中

dE=obj(NextX)-obj(PreX).

(22)

K是物理學(xué)中的波爾茲曼常數(shù)；exp表示自然指數(shù)，且dE>0.溫度越高，出現(xiàn)一次能量差為dE的降溫的概率就越大；溫度越低，則出現(xiàn)降溫的概率就越小.本質(zhì)上講，溫度T控制解的擾動(dòng)范圍.對(duì)于高溫度，更有可能使它達(dá)到一個(gè)更廣的范圍(Zhang and Yao, 2013).由于-dE/KT<0 ，所以P(dE)的函數(shù)取值范圍是(0,1).事實(shí)上，在給定的容差限Tolerance下將得到許多不同組解而不是像其他優(yōu)化方法所得到的唯一組解.因此可以進(jìn)一步通過一些精確的波數(shù)范圍和誤差之間權(quán)衡來選擇一組最優(yōu)的解(Zhang and Yao, 2013).

本文在離散化的假設(shè)下，利用優(yōu)化一階算子逼近真正的波數(shù)，進(jìn)而得到系數(shù)數(shù)組X=[a；c].通過以上計(jì)算流程，得到數(shù)組X=[ 0.09737101； -0.05866972]，即a=0.09737101,c= -0.05866972，此時(shí)M=1.當(dāng)選取高階差分時(shí)系數(shù)可見表3.事實(shí)上，由于模擬退火法求解的靈活性，也可以使用最大范數(shù)或者相對(duì)誤差等方式設(shè)定目標(biāo)函數(shù)(Zhang and Yao, 2013).

表3 優(yōu)化隱式交錯(cuò)網(wǎng)格有限差分系數(shù)Table 3 Coefficients of the optimized implicit staggered-grid finite-difference formulas

4 隱式顯式方法精度對(duì)比

傳統(tǒng)ESFDM精度分析表達(dá)式為(Liu and Sen, 2009b)

(23)

其中N為正整數(shù)，cn為相應(yīng)的差分系數(shù)，取值見Liu和Sen(2009b)中表1.

由公式(1)可以得到傳統(tǒng)ISFDM精度分析表達(dá)式為

(24)

其中M、N為正整數(shù)，am、bn為相應(yīng)的差分系數(shù)，取值見表1.

Liu和Sen(2009b)提出的ISFDM可表示為

(25)

其中，N為正整數(shù)，cn為相應(yīng)的差分系數(shù)，b為一個(gè)常量，取值見Liu 和 Sen(2009b)中表2.Liu(2014)采用ISFDM同公式(25)，在系數(shù)求取時(shí)采用最小二乘方法，取值見Liu(2014)中表9.

OISFDM精度分析表達(dá)式為

fOISFDM(β)≈

(26)

圖1 模擬退火法求解OISFD算子流程Fig.1 Flow chart of the ISFD scheme using the simulated annealing algorithm

其中a和c為常量，cm為相應(yīng)的差分系數(shù)，見表3.

由圖2—圖5中對(duì)比OISFD算子與其他差分算子的精度分析對(duì)比可以看出:(2+4)階OISFD算子的精度明顯高于傳統(tǒng)6階、8階ESFD算子甚至能達(dá)到10階ESFD算子精度.在與傳統(tǒng)ISFD算子精度的對(duì)比中，(2+4)階OISFD算子介于傳統(tǒng)ISFD(4+2)階與ISFD(6+2) 階算子之間.(2+4)階OISFD算子精度高于Liu和Sen(2009b)的(4+2) 階算子，低于Liu和Sen(2009b)的(6+2)階算子，與Liu(2014)的(4+2)階算子的精度接近.隨著差分階數(shù)的提高，OISFD算子的精度得到了明顯的改進(jìn).

5 數(shù)值模擬與分析

2D一階彈性波速度-應(yīng)力方程可以表示為(Virieux, 1986)

(27)

其中t為時(shí)間，x和z為坐標(biāo)方向，(vx,vz)為質(zhì)點(diǎn)的振動(dòng)速度，(τxx,τzz,τxz)為應(yīng)力的三個(gè)分量，λ和μ為拉梅常數(shù)，ρ為密度.

6 數(shù)值模擬實(shí)例

6.1 均勻介質(zhì)模型數(shù)值模擬

為了驗(yàn)證彈性介質(zhì)中OISFDM數(shù)值模擬的結(jié)果，本文首先設(shè)計(jì)了均勻模型，模型參數(shù)參見表4，分別采用傳統(tǒng)ESFDM、傳統(tǒng)ISFDM、Liu和Sen(2009b)的ISFDM和OISFDM進(jìn)行數(shù)值模擬.

圖2 一階導(dǎo)數(shù)傳統(tǒng)ESFD算子和OISFD算子精度對(duì)比Fig.2 Accuracy comparison between the conventional ESFD operators and OISFD operators for the first-order derivative

圖4 一階導(dǎo)數(shù)Liu 和 Sen(2009b)的 ISFD算子和 OISFD算子精度對(duì)比Fig.4 Accuracy comparison between Liu and Sen (2009b)′s ISFD operators and OISFD operators for the first-order derivative

圖5 一階導(dǎo)數(shù)Liu(2014)的 ISFD算子和 OISFD算子精度對(duì)比Fig.5 Accuracy comparison between Liu (2014)′sISFD operators and OISFD operators and for the first-order derivative

表4 均勻介質(zhì)模型正演模擬參數(shù)Table 4 Parameters of a homogeneous model for forward modeling

圖6為不同差分算法在計(jì)算時(shí)間為0.55 s處的波場快照，從圖中可以看到，各差分算法隨著差分階數(shù)的提高，頻散減弱，即模擬精度越高.通過快照對(duì)比可以發(fā)現(xiàn)傳統(tǒng)6階ESFDM會(huì)產(chǎn)生嚴(yán)重的波場畸變，在相同模擬參數(shù)條件下，由(2+4)階OISFDM引起的頻散則非常的微弱,可以得到與傳統(tǒng)10階ESFDM相當(dāng)或更高的模擬精度；(2+4)階OISFDM模擬結(jié)果與傳統(tǒng)(4+2)階 ISFDM相比波形保持得比較好，頻散相對(duì)較弱，與傳統(tǒng)(4+4)階ISFDM模擬精度相近，介于Liu 和 Sen(2009b)的(4+2)階ISFDM與Liu 和 Sen(2009b)的 (6+2)階ISFDM之間.

6.2 復(fù)雜介質(zhì)模型數(shù)值模擬

為了測試本文優(yōu)化隱式差分格式在復(fù)雜介質(zhì)波場模擬中的適應(yīng)能力，我們對(duì)修改的Sigsbee2A模型(如圖7)進(jìn)行數(shù)值模擬，模型參數(shù)見表5.

圖6 均勻模型0.55 s波場快照(左側(cè)為水平分量，右側(cè)為垂直分量) (a)(2+4)階OISFDM；(b)傳統(tǒng)6階 ESFDM；(c)傳統(tǒng)10階 ESFDM；(d)傳統(tǒng)(4+2)階ISFDM；(e)傳統(tǒng)(4+4)階ISFDM． (f)Liu 和 Sen(2009b)的(4+2)階ISFDM；(g)Liu 和 Sen(2009b)的(6+2)階ISFDM．Fig.6 Homogeneous model sanpshots at 0.55 s(Left:horizontal component. Right: vertical component) (a)(2+4) th-order OISFDM；(b)Conventional 6th-order ESFDM；(c)Conventional 10th-order ESFDM;(d)Conventional (4+2) th-order ISFDM；(e)Conventional (4+4) th-order ISFDM;(f) Liu and Sen (2009b)′s (4+2)th-order ISFDM；(g)Liu and Sen (2009b)′s (6+2)th- order ISFDM．

圖8分別給出了不同差分算法水平分量單炮記錄.圖10分別給出了不同差分算法垂直分量單炮記錄.為了對(duì)不同算法數(shù)值模擬的結(jié)果進(jìn)行比較，對(duì)圖8和圖10虛線框部分分別進(jìn)行局部放大得到圖9和圖11.從圖9和圖11中可以看到傳統(tǒng)6階ESFDM、傳統(tǒng)(4+2)階 ISFDM存在較為明顯的數(shù)值頻散.并且(2+4)階OISFDM可以得到與傳統(tǒng)(4+4)階ISFDM相當(dāng)有時(shí)甚至更佳的計(jì)算精度，介于Liu 和 Sen(2009b)的(4+2)階 ISFD算法(6+2)階 ISFD算法之間.

圖12分別給出了時(shí)間為0.92 s不同差分算法的波場快照.通過對(duì)比炮記錄和波場快照，可以看到(2+4)階OISFDM能夠較為有效的壓制數(shù)值頻散.將傳統(tǒng)10階ESFDM與(2+4)階OISFDM進(jìn)行對(duì)比時(shí)，它們的數(shù)值模擬精度幾乎是一樣的；傳統(tǒng)(4+2)階ISFDM會(huì)產(chǎn)生比(2+4)階OISFDM嚴(yán)重的數(shù)值頻散，當(dāng)將(2+4)階OISFDM與傳統(tǒng)(4+4)階ISFDM進(jìn)行對(duì)比時(shí)，(2+4)階OISFDM可以得到較高的精度；Liu 和 Sen(2009b)的(4+2)階 ISFDM會(huì)產(chǎn)生嚴(yán)重的數(shù)值頻散.然而，(2+4)階OISFDM精度介于Liu 和 Sen(2009b)的(4+2)階 ISFDM與(6+2)階 ISFDM之間.

圖7 復(fù)雜介質(zhì)速度模型Fig.7 Velocity model in complex medium

同時(shí)，通過圖13中4000個(gè)時(shí)間步長的計(jì)算時(shí)間對(duì)比，可以得出：(2+4)階OISFDM所需的CPU計(jì)算時(shí)間較少，具有更高的計(jì)算效率.此外，由公式(1)中N值更高的傳統(tǒng)ISFDM，在相同M值的情況下所帶來的精度提升空間有限，并不能顯著提升計(jì)算精度，反而因?yàn)榍竽婢仃噹捲黾訒?huì)明顯增加計(jì)算成本，降低計(jì)算效率.

圖8 復(fù)雜介質(zhì)水平分量單炮記錄 其他說明同圖6.Fig.8 Horizontal component in seismograms for complex medium Other the same as the Fig.6

圖9 復(fù)雜介質(zhì)水平分量炮記錄局部放大 其他說明同圖6.Fig.9 Zoom of horizontal component in seismograms for complex medium Other the same as the Fig.6

圖10 復(fù)雜介質(zhì)垂直分量單炮記錄 其他說明同圖6.Fig.10 Vertical component in seismograms for complex medium Other the same as the Fig.6

圖11 垂直分量單炮記錄局部放大 其他說明同圖6.Fig.11 Zoom of vertical component in seismograms for complex medium Other the same as the Fig.6.

7 結(jié)論

本文建立了一種空間OISFD格式，并在此基礎(chǔ)上應(yīng)用二范數(shù)原理建立絕對(duì)誤差的目標(biāo)函數(shù)；然后通過模擬退火法對(duì)目標(biāo)函數(shù)進(jìn)行多參數(shù)的優(yōu)化，以求取優(yōu)化系數(shù)；最后采用該優(yōu)化差分法對(duì)彈性波場進(jìn)行了數(shù)值模擬與分析，還與傳統(tǒng)的幾種交錯(cuò)網(wǎng)格差分方法進(jìn)行了對(duì)比.在均勻介質(zhì)以及復(fù)雜介質(zhì)中，對(duì)傳統(tǒng)ESFDM、傳統(tǒng)ISFDM、Liu 和 Sen(2009b)的ISFDM和OISFDM進(jìn)行數(shù)值模擬，通過單炮記錄、波場快照、計(jì)算效率等角度分析可以看出：在相同模擬參數(shù)的前提下，(2+4)階OISFDM可以比傳統(tǒng)6階 ESFDM更為有效的壓制數(shù)值頻散并且能達(dá)到傳統(tǒng)10階ESFDM的精度；在與傳統(tǒng)ISFDM的對(duì)比中可以看出，(2+4)階OISFDM較傳統(tǒng)(4+2)階ISFDM能夠更好的壓制頻散，可以得到比傳統(tǒng)(4+4)階ISFDM略高的精度；在與Liu 和Sen(2009b)的ISFDM的對(duì)比中，(2+4)階OISFDM可以達(dá)到介于Liu 和 Sen(2009b)的(4+2)階ISFDM與(6+2)階ISFDM之間的模擬精度.在高精度的數(shù)值模擬中可以應(yīng)用高階的OISFDM來達(dá)到更佳的模擬效果.因此，OISFDM具有更高的精度和計(jì)算效率，為彈性波正演模擬提供了一種新的數(shù)值算法.

此外，當(dāng)將OISFDM應(yīng)用于3D復(fù)雜介質(zhì)時(shí)，計(jì)算量和內(nèi)存需求也急劇增大，成為影響數(shù)值模擬計(jì)算效率的瓶頸.因此可以考慮引入GPU/CPU協(xié)同并行架構(gòu)，將模型在慢維方向分塊并行，在每塊中通過GPU實(shí)現(xiàn)每個(gè)時(shí)間步的波場數(shù)值計(jì)算，然而這一思想仍需要在下一步的研究中詳細(xì)討論實(shí)現(xiàn).

致謝 感謝嚴(yán)紅勇博士對(duì)本論文工作的建議和幫助.

Baysal E, Kosloff D, Sherwood W C. 1983. Reverse time migration.Geophysics, 48(11): 1514-1524.Chu C. 2009. Seismic modeling and imaging with the fourier method: numerical analyses and parallel implementation strategies[Ph. D. thesis]. Univ. of Texas.Chu C L, Stoffa P L. 2010. Frequency domain modeling using implicit spatial finite difference operators. // SEG Technical Program Expanded Abstracts, 3076-3080,

Chu C L, Stoffa P L. 2012. Implicit finite-difference simulations of seismic wave propagation.Geophysics, 77(2): T57-T67, doi: 10.1190/geo2011-0180. 1.

Claerbout J F. 1985. The craft of wavefield extrapolation. // Imaging the Earth′s Interior, 260-265, Oxford: Blackwell Scientific Publications.

Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid high-order difference method of one-order elastic wave equation.ChineseJ.Geophys. (in Chinese), 43(3): 411-419.

Ekaterinaris J A. 1999. Implicit, high-resolution, compact schemes for gas dynamics and aero acoustics.JournalofComputationalPhysics, 156(2): 272-299, doi: 10.1006/jcph.1999.6360. Fornberg B. 1998. Calculation of weights in finite difference formulas.SIAMReview, 40(3): 685-691, doi: 10.1137/S0036144596322507. Kosloff D, Pestana R, Tal-Ezer H. 2008. Numerical solution of the constant density acoustic wave equation by implicit spatial derivative operators. // 78th Annual International Meeting, SEG Expanded Abstracts, 2057-2061.

Kosloff D, Pestana R, Tal-Ezer H. 2010. Acoustic and elastic numerical wave simulations by recursive spatial derivative operators.Geophysics, 75(1): T167-T174, doi: 10.1190/1.348521.

Lele S K. 1992. Compact finite difference schemes with spectral-like resolution.JournalofComputationalPhysics, 103(1): 16-42, doi: 10.1016/0021-9991(92)90324-R.

Liu Y, Sen M K. 2009a. A practical implicit finite-difference method: examples from seismic modeling.JournalofGeophysicsandEngineering, 6(3): 231-249, doi: 10.1088/1742-2132/6/3/003.Liu Y, Sen M K. 2009b. An implicit staggered-grid finite-difference method for seismic modeling.GeophysicalJournalInternational, 179(1): 459-474, doi: 10.1111/j.1365-246X. 2009.04305.x.Liu Y. 2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling.GeophysicalJournalInternational, 197(2): 1033-1047, doi: 10. 1093/gji/ggu032.

Liu Y, Li C C, Mou Y G. 1998. Finite-difference numerical modeling of any even-order accuracy.OilGeophysicalProspecting, 33(1): 1-10.Liu Y, Yan H Y, Liu H. 2014. Least squares staggered-grid finite-difference for elastic wave modelling.ExplorationGeophysics, 45(4): 255-260. doi: 10.1071/EG13087.

McMechan G A. 1983. Migration by extrapolation of time-dependent boundary values.GeophysicalProspecting, 31(3): 413-420, doi: 10.1190/geo2012-0338.1.

Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. 1953. Equation of state calculations by fast computing machines.JournalofChemicalPhysics, 21(6): 1087-1092.

Pei Z L. 2005. Numerical simulation of elastic wave equation in 3-D isotropic media with staggered-grid high-order difference method.GeophysicalProspectingforPetroleum, 44(4): 308-315.Shang J S. 1999. High-order compact-difference schemes for time dependent Maxwell equations.JournalofComputationalPhysics, 153(2): 312-333, doi: 10.1006/jcph.1999. 6279.Virieux J. 1984. SH-wave propagation in heterogeneous media: velocity stress finite-difference method.Geophysics, 49(11): 1933-1957.

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

Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6): WCC1-WCC26.

Warner M, Ratclife A, Nangoo T, et al. 2013. Anisotropic 3D full-waveform inversion.Geophysics, 78(2): R59-R80, doi: 10.1190/geo2012-0338.1.

Yan H Y, Liu Y. 2012. Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media.ChineseJ.Geophys. (in Chinese), 55(4): 1354-1365, doi: 10.6038/j. issn. 0001-5733.2012.04.031.

Yan H Y, Liu Y. 2013. Visco-acoustic pre-stack reverse-time migration based on the time-space domain adaptive high-order finite-difference method.GeophysicalProspecting, 61(5): 941-954, doi: 10.1111/1365-2478.12046.

圖12 0.92s波場快照(左側(cè)為水平分量，右側(cè)為垂直分量) 其他說明同圖6.Fig.12 Snapshots at 0.92 s(Left:horizontal component. Right: vertical component) Other the same as the Fig.6.

圖13 測試時(shí)間對(duì)比 (a) OISFDM與傳統(tǒng)ESFDM對(duì)比； (b) OISFDM與傳統(tǒng)ISFDM對(duì)比; (c) OISFDM與Liu 和 Sen(2009b)的ISFDM對(duì)比．Fig.13 Comparison of computing costs (a) Comparison between OISFDM and conventional ESFDM; (b) Comparison between OISFDM and conventional ISFDM; (c) Comparison of between OISFDM and Liu and Sen (2009b)′s ISFDM．

Zhang H, Zhang Y, Sun J. 2007. Implicit splitting finite difference scheme for multi-dimensional wave simulation. // SEG Technical Program Expanded Abstracts, 2011-2015, doi: 10.1190/1.2792885.Zhang J H, Yao Z X. 2013. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm.Journalof Computational Physics, 250: 511-526,doi: 10.1016/j.jcp. 2013.04.029.

ZhouHB,ZhangGQ. 2011.Prefactoredoptimizedcompactfinitedifferenceschemesforsecondspatialderivatives. Geophysics, 76(5):WB87-WB95,doi: 10.1190/geo2011-0048.1.

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

董良國, 馬在田, 曹景忠等. 2000. 一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法. 地球物理學(xué)報(bào), 43(3): 411-419.

劉洋, 李承楚, 牟永光. 1998. 任意偶數(shù)階精度有限差分法數(shù)值模擬. 石油地球物理勘探, 33(1): 1-10.

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

嚴(yán)紅勇, 劉洋. 2012. 黏彈TTI介質(zhì)中旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬. 地球物理學(xué)報(bào), 55(4): 1354-1365, doi: 10.6038/j.issn.0001-5733.2012.04.031.

(本文編輯 張正峰)

A global optimized implicit staggered-grid finite-difference scheme for elastic wave modeling

WANG Yang1,2， LIU Hong1*， ZHANG Heng3， WANG Zhi-Yang1,2， TANG Xiang-De1,2

1KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China3MLRKeyLaboratoryofMarineMineralResources,GuangzhouMarineGeologicalSurvey,Guangzhou510075,China

The advanced algorithms, such as reverse time migration (RTM) and full waveform inversion (FWI), require iterative calculation of the wave field. Its wave propagation part is highly sensitive to computational accuracy and stability. Therefore, how to find a way to obtain a finite-difference scheme that can consider the both aspects is absolutely critical. Here we propose a global optimized implicit staggered-grid finite-difference scheme that inherits the merit of high stability from the implicit finite-difference scheme and improves the computational efficiency remarkably.The conventional implicit finite-difference operators generally have to obtain the high accuracy solutions from the inversion of either pentadiagonal or more diagonals matrices using the LU decomposition method, which greatly reduce the computation efficiency. Therefore, we first perform a global optimized implicit staggered-grid finite-difference (OISFD) scheme for first-order derivatives, which only needs to solve twice tridiagonal matrix inversion using the chasing method. In this way, it avoids the requirement of solving the large-scale matrix equation. Then we transform the new scheme from the time-space domain into the time-wavenumber domain, and develop the objective function by means of the two-norm theory. Different from traditional global optimized methods, we select the simulated annealing algorithm to get the optimized coefficients. And finally, after introducing the coefficients to the operator we establish the global optimized implicit staggered-grid finite-difference scheme.We use independent tests to show that for same size matrix inversion, the average CPU time of the chasing method occupies 0.0019～0.0058 times that of the LU decomposition method. And this comparison becomes the theoretical basis of improving the computational efficiency. From the accuracy curves of approximation between the conventional explicit staggered-grid finite-difference method (ESFDM) and implicit staggered-grid finite-difference method (ISFDM), we can indicate that the OISFD operators have wider frequency coverage and smaller accuracy error fluctuation than those of the conventional ESFD operators and ISFD operators when using the same order spatial difference operator. In order to demonstrate the effect in elastic wave modeling, we perform first-order velocity-stress elastic wave equation numerical simulations on a homogeneous model and the modified Sigsbee2A model. The resulting shot records and wave field snapshots indicate that the OISFDM is more effective in suppressing the numerical dispersion, which is fully consistent with the theoretical accuracy curve comparison. What′s more, the least CPU time cost in numerical simulation further confirms that the OISFD scheme reduces calculating amount while keeping high accuracy numerical simulation performance.We have presented an implicit optimized scheme with any order of accuracy for calculating first spatial derivatives on a staggered-grid. The operator is designed by fitting the response in the wavenumber domain using the simulated annealing algorithm. The method is based on applying optimized operators to each of the spatial coordinates in elastic wave equation. It could utilize smaller stencils to achieve more accurate results than the corresponding explicit operator without adding significantly higher computational cost to conventional implicit finite-difference methods. The dispersion relation analysis and the numerical results demonstrate that the OISFD operators can achieve higher accuracy compared with those using the conventional ESFD operators and ISFD operators. This means we can greatly save the memory and computational cost required when using our optimized OISFD methods.

Implicit finite-difference； Numerical simulations； Staggered-grid； Elastic wave equation； Simulated annealing algorithm

國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)項(xiàng)目(2012AA061202)，國家油氣重大專項(xiàng)(2011ZX05008-006-50)，國家油氣重大專項(xiàng)(2011ZX05003-003)，國家科技重大專項(xiàng)(2011ZX05023-005-002)聯(lián)合資助.

王洋，女，1989年生，中國科學(xué)院地質(zhì)與地球物理研究所在讀博士研究生，主要從事地震波成像方法的研究. E-mail: wangyang2010@mail.iggcas.ac.cn

*通訊作者 劉洪，男，1959年生，研究員，主要從事油儲(chǔ)地球物理研究．E-mail: liuhong@mail.iggcas.ac.cn

10.6038/cjg20150726

P631

2014-06-11，2015-05-16收修定稿

10.1190/1.3513485.

王洋， 劉洪， 張衡等.2015.一種全局優(yōu)化的隱式交錯(cuò)網(wǎng)格有限差分算法及其在彈性波數(shù)值模擬中的應(yīng)用.地球物理學(xué)報(bào)，58(7):2508-2524,doi:10.6038/cjg20150726.

Wang Y, Liu H, Zhang H, et al. 2015. A global optimized implicit staggered-grid finite-difference scheme for elastic wave modeling.ChineseJ.Geophys. (in Chinese),58(7):2508-2524,doi:10.6038/cjg20150726.