李世中,孫成禹,彭鵬鵬

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

有限差分方法廣泛應(yīng)用于地震數(shù)值模擬中。20世紀(jì)70年代,有限差分法首次應(yīng)用于波動(dòng)方程正演模擬[1]。之后,VIRIEUX[2]提出了穩(wěn)定的二階交錯(cuò)網(wǎng)格彈性波有限差分格式。當(dāng)用固定步長的有限差分格式對地下低速帶或者復(fù)雜構(gòu)造進(jìn)行數(shù)值模擬時(shí),必須選擇較小的網(wǎng)格步長來滿足計(jì)算精度要求或者避免漏掉重要的波場信息,因而不可避免地增加了模擬的計(jì)算成本。1989年,MOCZO[3]首次提出采用可變網(wǎng)格思想來提高有限差分?jǐn)?shù)值模擬精度。隨后,JASTRAM等[4-5]提出了網(wǎng)格步長可取任意倍數(shù)的可變網(wǎng)格算法。李勝軍等[6]比較了傳統(tǒng)網(wǎng)格與變網(wǎng)格有限差分算法的內(nèi)存需求、計(jì)算效率等,總結(jié)了變網(wǎng)格算法的優(yōu)點(diǎn)。李振春等[7]從理論上推導(dǎo)了變網(wǎng)格界面處虛假反射誤差的數(shù)學(xué)表達(dá)式,并通過引入Lanczos濾波算子,實(shí)現(xiàn)了一種基于交錯(cuò)網(wǎng)格的穩(wěn)定且高精度的雙變網(wǎng)格正演模擬算法。孫成禹等[8]證實(shí)了變網(wǎng)格有限差分算法能提高模擬結(jié)果的分辨率,同時(shí)降低內(nèi)存需求量,減少計(jì)算時(shí)間。孟凡順等[9]將可變網(wǎng)格有限差分方法引入到交錯(cuò)網(wǎng)格高階差分?jǐn)?shù)值模擬中,并將可變網(wǎng)格步長倍數(shù)變化范圍推廣至任意奇數(shù)倍。KANG等[10]基于一階速度應(yīng)力方程,進(jìn)一步改進(jìn)有限差分的變時(shí)間步長非連續(xù)空間交錯(cuò)網(wǎng)格方法。LI等[11]實(shí)現(xiàn)了一種靈活的非均勻時(shí)間步長不連續(xù)曲線網(wǎng)格有限差分方法,用于地震波數(shù)值模擬。孫林潔等[12]實(shí)現(xiàn)了空間與時(shí)間步長任意整數(shù)倍變化的雙相介質(zhì)彈性波可變網(wǎng)格正演模擬。馬光克等[13]將可變網(wǎng)格有限差分算法引入到逆時(shí)偏移中,并分析了可變網(wǎng)格對有限差分算法數(shù)值頻散的影響。

采用有限差分法進(jìn)行數(shù)值模擬時(shí),由于利用離散網(wǎng)格點(diǎn)的差分來代替空間偏導(dǎo)數(shù),不可避免地產(chǎn)生數(shù)值頻散現(xiàn)象。許多學(xué)者對此做了大量研究,旨在提高模擬的精度。ALFORD等[14]在研究聲波方程有限差分模擬精度時(shí)指出,為消除數(shù)值頻散,對于時(shí)間2階、空間4階精度的有限差分法,在源子波的Nyquist頻率所對應(yīng)的一個(gè)波長內(nèi)不能少于5.5個(gè)網(wǎng)格節(jié)點(diǎn)。董良國等[15]指出影響數(shù)值頻散的因素為地震波傳播方向、差分精度和一個(gè)波長內(nèi)離散點(diǎn)數(shù)。為減少數(shù)值頻散,可以通過增加有限差分算子的長度去計(jì)算空間偏導(dǎo)數(shù)[16],但會(huì)增加計(jì)算成本。TAM等[17]給出了顯式頻散關(guān)系保存方案。LIU等[18]通過使用Taylor級(jí)數(shù)展開法來最小化時(shí)空域頻散關(guān)系的誤差,推導(dǎo)出有限差分系數(shù),以此來減少數(shù)值頻散。杜啟振等[19]通過引入強(qiáng)約束條件和弱約束條件,構(gòu)造了不同的Lagrange函數(shù),然后通過求取條件極值得到優(yōu)化差分算子。梁文全等[20]采用前人提出的新的有限差分模板(在保持相同精度的情況下增大了時(shí)間步長),基于非線性優(yōu)化確定時(shí)間-空間域隱格式有限差分系數(shù)。YANG等[21]基于頻散關(guān)系和最小二乘理論推導(dǎo)出一階空間導(dǎo)數(shù)的差分系數(shù)。YAN等[22]提出一種基于最小二乘理論的優(yōu)化差分算子進(jìn)行聲波方程正演和逆時(shí)偏移。IGEL等[23]使用高斯窗截?cái)嗟玫接邢薏罘炙阕酉禂?shù)。雍鵬等[24]給出一種顯式時(shí)間遞推格式,并采用共軛梯度法得到精確時(shí)間遞推匹配系數(shù),實(shí)現(xiàn)時(shí)空差分算子的同時(shí)優(yōu)化。梁文全等[25]提出使用線性方法壓制聲波方程矩形網(wǎng)格有限差分算子的數(shù)值頻散。此外,還有一些關(guān)于有限差分算子的優(yōu)化方法,如牛頓法[26]、模擬退火法[27]、樣點(diǎn)逼近法[28]等,都能不同程度地降低較大波數(shù)范圍內(nèi)的數(shù)值頻散,提高地震波數(shù)值模擬精度。

以上差分系數(shù)優(yōu)化都是針對固定步長的網(wǎng)格,關(guān)于變網(wǎng)格差分系數(shù)優(yōu)化的研究相對較少。本文在前人研究的基礎(chǔ)上,提出一種基于最小二乘理論計(jì)算可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)算法。通過頻散分析和正演模擬算例表明,與Taylor系數(shù)展開所獲得的差分系數(shù)相比,采用優(yōu)化后的可變交錯(cuò)網(wǎng)格差分系數(shù),能更有效地壓制頻散現(xiàn)象。

1 方法原理

二維各向同性介質(zhì)中的一階應(yīng)力-速度聲波方程可表示為:

(1)

式中:vx,vz分別為質(zhì)點(diǎn)在x和z方向上的速度分量;t為時(shí)間;P為應(yīng)力;vP為縱波速度;ρ為介質(zhì)密度。

1.1 變網(wǎng)格差分算法

常規(guī)網(wǎng)格采用固定的網(wǎng)格步長對模擬區(qū)域進(jìn)行剖分,而變網(wǎng)格一般采用不同空間步長來對模擬區(qū)域進(jìn)行離散,如圖1所示。圖1中紅色線框區(qū)域?yàn)樾〔介L網(wǎng)格區(qū)域,黑色線框區(qū)域?yàn)榇蟛介L網(wǎng)格區(qū)域。當(dāng)對大小步長網(wǎng)格邊界點(diǎn)做插值處理時(shí),一方面會(huì)降低數(shù)值模擬效率,另一方面也會(huì)積累計(jì)算誤差,影響模擬精度。在大小網(wǎng)格邊界處引入過渡帶,能避免因插值計(jì)算引起的誤差,使得數(shù)值模擬時(shí),網(wǎng)格長度能自然過渡,增加穩(wěn)定性。

為了保證過渡帶的所有半節(jié)點(diǎn)導(dǎo)數(shù)不用插值就可以求出高階差分?jǐn)?shù)值,大小步長網(wǎng)格的縱向步長比需為奇整數(shù)倍。以縱向上大小網(wǎng)格步長比3∶1為例,圖2給出了可變交錯(cuò)網(wǎng)格在過渡帶的差分處理方式[13]。

圖1 差分網(wǎng)格分布

圖2 過渡帶網(wǎng)格變化及波場計(jì)算

以一階導(dǎo)數(shù)?P/?z的可變交錯(cuò)網(wǎng)格差分形式為例,大網(wǎng)格區(qū)域的高階差分計(jì)算公式為:

(2)

小網(wǎng)格區(qū)域差分格式:

(3)

網(wǎng)格過渡帶差分格式:

(4)

1.2 優(yōu)化差分算子

對于可變交錯(cuò)網(wǎng)格聲波波動(dòng)方程,一般采用時(shí)間二階、空間高階有限差分算子來提高數(shù)值模擬的精度。對于大小網(wǎng)格過渡帶的波場f(x)來說,其空間一階導(dǎo)數(shù)交錯(cuò)網(wǎng)格差分格式可近似地表示為:

(5)

假設(shè)f(x)為沿x方向傳播的單頻波,有:

(6)

(7)

其中β=kh/2,且β∈(0,0.5π)。運(yùn)用Taylor公式將(7)式的三角函數(shù)展開為多項(xiàng)式,可得:

(8)

比較β的系數(shù),得:

(9)

首先對(7)式構(gòu)建基于頻散關(guān)系的平方誤差函數(shù),有:

(10)

式中:b是積分上限,其取值與K有關(guān)。在給定的積分上限取值范圍[d,u]內(nèi),均勻分布著N個(gè)采樣點(diǎn),即b1,b2,…,bN。當(dāng)K確定時(shí),積分上限b取對應(yīng)的采樣點(diǎn)值bK。

(10)式表示可變交錯(cuò)網(wǎng)格有限差分法求取空間導(dǎo)數(shù)時(shí)在給定間隔[0,b]內(nèi)引入的平方誤差,由最小二乘理論,求取誤差的最小值。此外,加上近零波數(shù)約束條件以提高可變交錯(cuò)網(wǎng)格有限差分法的數(shù)值模擬精度[21]。對于可變交錯(cuò)網(wǎng)格,有:

(11)

令

(12)

由(10)式和(12)式可得極值的必要條件為:

(13)

式中:λ為Lagrange乘數(shù)。將(10)式、(12)式代入(13)式中,最終得到:

(14)

其中,

(15)

1.3 頻散分析

根據(jù)公式(7),定義以下公式描述可變交錯(cuò)網(wǎng)格有限差分的頻散關(guān)系:

(16)

當(dāng)δ(β)的值越接近0時(shí),數(shù)值頻散就越小;當(dāng)δ(β)的值越遠(yuǎn)離0時(shí),數(shù)值頻散就越大。為了對比可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)和Taylor級(jí)數(shù)展開所得系數(shù)的頻散關(guān)系,以N=5為例,取積分上限取值范圍u=0.4π,將可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)和Taylor級(jí)數(shù)展開所得系數(shù)分別代入(16)式中,得到不同大小網(wǎng)格步長比l下的δ(β)值隨β變化的曲線,如圖3所示。附錄A給出了N=5的可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)和Taylor級(jí)數(shù)展開的差分系數(shù)。

從圖3可以看出,相同的可變交錯(cuò)網(wǎng)格差分算子長度下,不同的大小網(wǎng)格步長比l對應(yīng)的積分上限取值范圍d值也不同,l越大,d越小。隨著過渡帶中小網(wǎng)格內(nèi)的點(diǎn)逐漸遠(yuǎn)離大小網(wǎng)格交界處,優(yōu)化的效果越好。因此,可變交錯(cuò)網(wǎng)格優(yōu)化差分算子相比基于Taylor級(jí)數(shù)展開的差分算子能夠減少差分近似帶來的誤差,更好地壓制數(shù)值頻散。下面將對大小網(wǎng)格步長比l=5時(shí)的可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)和基于Taylor級(jí)數(shù)展開的差分系數(shù)分別進(jìn)行數(shù)值模擬和對比驗(yàn)證。

圖3 大小網(wǎng)格邊界區(qū)域差分精度分析及對比a l=3,d=0.18π; b l=5,d=0.12π; c l=9,d=0.07π; d l=15,d=0.05π

2 數(shù)值模擬

2.1 低速夾層模型

如圖4所示的一個(gè)低速夾層模型寬為4500m,高為4500m。由地表往下各層的速度分別為3800,2500,4200m/s。第2層為低速層,上、下界面深度分別為2355m和2400m。對模型分別利用可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法和Taylor級(jí)數(shù)展開法進(jìn)行聲波方程正演模擬,模型虛線框外為大網(wǎng)格區(qū)域,步長15m×15m,虛線框內(nèi)為小網(wǎng)格區(qū)域,步長Δx=15m,Δz=3m。再用常規(guī)交錯(cuò)網(wǎng)格(15m×15m,3m×3m)的Taylor級(jí)數(shù)展開法對低速夾層模型進(jìn)行聲波方程正演模擬。介質(zhì)密度從上往下依次為1700,1500,1800kg/m3,時(shí)間采樣間隔均為0.1ms,采用時(shí)間2階、空間10階有限差分算子進(jìn)行數(shù)值模擬。震源點(diǎn)坐標(biāo)為(2250m,1800m),子波為雷克子波,主頻為40Hz。

圖4 低速夾層模型

圖5是低速夾層模型正演模擬在0.36s的波場快照,可以看出,當(dāng)采用固定步長的大網(wǎng)格模擬時(shí),由于低速層的影響,與圖5a,圖5b和圖5d相比,圖5c出現(xiàn)非常明顯的數(shù)值頻散。再對比圖5a,圖5b和圖5d 可以看出,即使采用了基于Taylor級(jí)數(shù)展開的變網(wǎng)格有限差分?jǐn)?shù)值模擬,也出現(xiàn)了較嚴(yán)重的數(shù)值頻散,干擾了正常波場。而采用可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法進(jìn)行數(shù)值模擬時(shí),精度更高,壓制了頻散。

圖5 低速夾層模型聲波方程波場快照a 可變交錯(cuò)網(wǎng)格Taylor級(jí)數(shù)展開法(3m/15m); b 可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法(3m/15m); c 常規(guī)大網(wǎng)格Taylor級(jí)數(shù)展開法(15m×15m); d 常規(guī)小網(wǎng)格Taylor級(jí)數(shù)展開法(3m×3m)

圖6是低速夾層模型正演模擬在檢波點(diǎn)(2250m,1800m)的部分接收記錄,以常規(guī)小網(wǎng)格的波形記錄為準(zhǔn)確值,可以算出可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法波形記錄與準(zhǔn)確值的均方根誤差為0.002289;而可變交錯(cuò)網(wǎng)格Taylor級(jí)數(shù)展開法波形記錄與準(zhǔn)確值的均方根誤差為0.008089。均方根誤差定義為:

圖6 低速夾層模型正演模擬在檢波點(diǎn)(2250m,1800m)的部分接收記錄

(17)

2.2 起伏界面模型

如圖8所示的起伏界面模型寬、高均為4000m,分為4層,第1層速度為1500m/s,第2層速度為2500m/s,第3層速度為2800m/s,第4層速度為3400m/s。第1層為低速層,且第1層和第2層的分界面是起伏的。因此,進(jìn)行可變交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬時(shí),令深度600m以上為小網(wǎng)格區(qū)域,其步長Δx=10m,Δz=2m,深度600m以下為大網(wǎng)格區(qū)域,步長10m×10m。再用常規(guī)網(wǎng)格(10m×10m,2m×2m)的Taylor級(jí)數(shù)展開法對起伏界面模型進(jìn)行聲波方程正演模擬。介質(zhì)密度設(shè)為常數(shù)2000kg/m3,時(shí)間采樣間隔均為0.1ms。震源點(diǎn)坐標(biāo)為(2000m,10m),子波為雷克子波,主頻為24Hz。

圖7 低速夾層模型聲波方程地震記錄a 可變交錯(cuò)網(wǎng)格Taylor級(jí)數(shù)展開法(3m/15m); b 可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法(3m/15m); c 常規(guī)大網(wǎng)格Taylor級(jí)數(shù)展開法(15m×15m); d 常規(guī)小網(wǎng)格Taylor級(jí)數(shù)展開法(3m×3m)

圖8 起伏界面模型

圖9是空間10階有限差分算子起伏界面模型的正演模擬地震記錄,可以明顯看出,采用常規(guī)大網(wǎng)格進(jìn)行正演模擬得到的地震記錄數(shù)值頻散最嚴(yán)重,采用可變交錯(cuò)網(wǎng)格Taylor級(jí)數(shù)展開法得到的地震記錄也出現(xiàn)了相當(dāng)程度的數(shù)值頻散,不利于實(shí)際的波場分析;而在同樣條件下采用可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法得到的地震記錄出現(xiàn)的數(shù)值頻散較少。因此,在采用相同的空間差分算子長度下,可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法相比Taylor級(jí)數(shù)展開法能有效降低數(shù)值模擬的頻散誤差。此外,穩(wěn)定性條件是有限差分?jǐn)?shù)值模擬中十分重要的問題。對于可變交錯(cuò)網(wǎng)格,若是網(wǎng)格步長變化梯度太劇烈,不僅會(huì)引起人為的數(shù)值干擾,造成計(jì)算誤差,甚至?xí)钩绦蛴?jì)算溢出,無法運(yùn)行。故空間網(wǎng)格出現(xiàn)太大的步長變化時(shí),往往需要平緩的、經(jīng)過幾次步長變化的處理來保證長時(shí)程計(jì)算的穩(wěn)定性。這里長時(shí)程的地震記錄也間接證明了可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法計(jì)算的穩(wěn)定性,適用于大區(qū)域數(shù)值模擬。

圖9 空間10階差分算子起伏界面模型正演模擬地震記錄a 可變交錯(cuò)網(wǎng)格Taylor級(jí)數(shù)展開法(2m/10m); b 可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法(2m/10m); c 常規(guī)大網(wǎng)格Taylor級(jí)數(shù)展開法(10m×10m); d 常規(guī)小網(wǎng)格Taylor級(jí)數(shù)展開法(2m×2m)

由于過渡帶的優(yōu)化差分算子與空間偏導(dǎo)數(shù)算子之間仍存在誤差,因此,如果需要進(jìn)一步降低頻散,可以通過提高空間差分階數(shù)來實(shí)現(xiàn),但提高差分階數(shù)將會(huì)增大數(shù)值模擬計(jì)算量,從而增加計(jì)算成本。圖10是空間18階差分算子起伏界面模型可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法地震記錄,與圖9b相比,數(shù)值頻散得到了進(jìn)一步壓制。圖11是起伏界面模型數(shù)值模擬計(jì)算成本對比結(jié)果,可以看出,當(dāng)空間差分階數(shù)相同時(shí),可變交錯(cuò)網(wǎng)格的優(yōu)化差分系數(shù)法與可變交錯(cuò)網(wǎng)格Taylor級(jí)數(shù)展開法計(jì)算時(shí)間幾乎一樣,表明用優(yōu)化差分算子解波動(dòng)方程不會(huì)增加額外的計(jì)算成本;而隨著空間差分階數(shù)的增大,計(jì)算成本也會(huì)增加;采用常規(guī)小網(wǎng)格做數(shù)值模擬,計(jì)算成本將急劇增加。

圖10 空間18階差分算子起伏界面模型正演模擬地震記錄(積分上限取值范圍d=0.12π,u=0.42π)

圖11 起伏界面模型數(shù)值模擬計(jì)算成本對比

3 結(jié)束語

在常規(guī)交錯(cuò)網(wǎng)格差分系數(shù)優(yōu)化的基礎(chǔ)上,推導(dǎo)出可變交錯(cuò)網(wǎng)格的優(yōu)化差分算子系數(shù),既保證了變網(wǎng)格的計(jì)算效率,又能進(jìn)一步壓制數(shù)值頻散,便于波場特征分析。通過對變網(wǎng)格過渡區(qū)域優(yōu)化差分系數(shù)的求取和采用優(yōu)化系數(shù)進(jìn)行數(shù)值模擬,并將模擬結(jié)果與Taylor級(jí)數(shù)展開法的模擬結(jié)果對比,可以得到以下認(rèn)識(shí)。

1) 可變交錯(cuò)網(wǎng)格差分算子的精度取決于差分系數(shù)和差分階數(shù),提高差分階數(shù)能壓制數(shù)值頻散,但也導(dǎo)致計(jì)算成本增加。與基于Taylor級(jí)數(shù)展開的可變交錯(cuò)網(wǎng)格差分算子相比,基于最小二乘理論的可變交錯(cuò)網(wǎng)格優(yōu)化差分算子在不增加計(jì)算量的情況下能有效地壓制數(shù)值頻散。

2) 可變交錯(cuò)網(wǎng)格優(yōu)化差分算子的長度不同時(shí),對應(yīng)的積分上限取值范圍u值也不同。相同算子長度下,不同的大小網(wǎng)格步長比l對應(yīng)的積分上限取值范圍d值也不同,一般來說,當(dāng)差分算子長度相同時(shí),l越大,d越小。

此外,由于時(shí)間步長未采用可變算法,為進(jìn)一步提高可變交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)法數(shù)值模擬的穩(wěn)定性,可以采用變空間步長和變時(shí)間步長相結(jié)合的方法。

參 考 文 獻(xiàn)

[1] ALTERMAN Z S,LOEWENTHAL D.Seismic wave in a quarter and three quarter plane[J].Geophysical Journal of the Royal Astronomical Society,1970,20(1):101-126

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

[3] MOCZO P.Finite-difference technique for SH-waves in 2-D media using irregular grids:application to the seismic response problem[J].Geophysical Journal International,1989,99(2):321-329

[4] JASTRAM C,BEHLE A.Acoustic modeling on a grid of vertically varying spacing[J].Geophysical Prospecting,1992,40(1):157-169

[5] JASTRAM C,TESSMER E.Elastic modeling on a grid with vertically varying spacing[J].Geophysical Prospecting,1994,42(2):357-370

[6] 李勝軍,孫成禹,倪長寬,等.聲波方程有限差分?jǐn)?shù)值模擬的變網(wǎng)格步長算法[J].工程地球物理學(xué)報(bào),2007,4(3):207-212

LI S J,SUN C Y,NI C K,et al.Acoustic equation numerical modeling on a grid of varying spacing[J].Chinese Journal of Engineering Geophysics,2007,4(3):207-212

[7] 李振春,李慶洋,黃建平,等.一種穩(wěn)定的高精度雙變網(wǎng)格正演模擬與逆時(shí)偏移方法[J].石油物探,2014,53(2):127-136

LI Z C,LI Q Y,HUANG J P,et al.A stable and high-precision dual-variable grid forward modeling and reverse time migration method[J].Geophysical Prospecting for Petroleum,2014,53(2):127-136

[8] 孫成禹,李勝軍,倪長寬,等.波動(dòng)方程變網(wǎng)格步長有限差分?jǐn)?shù)值模擬[J].石油物探,2008,47(2):123-128

SUN C Y,LI S J,NI C K,et al.Wave equation numerical modeling by finite difference method with varying grid spacing[J].Geophysical Prospecting for Petroleum,2008,47(2):123-128

[9] 孟凡順,李洋森,王璐,等.基于可變交錯(cuò)網(wǎng)格的地震波數(shù)值模擬[J].地球物理學(xué)進(jìn)展,2013,28(6):2936-2943

MENG F S,LI Y S,WANG L,et al.Based on variable and staggered grids seismic numerical simulation[J].Progress in Geophysics,2013,28(6):2936-2943

[10] KANG T S.Finite-difference seismic simulation combining discontinuous grids with locally variable timesteps[J].Bulletin of the Seismological Society of America,2004,94(1):207-219

[11] LI H,ZHANG W,ZHANG Z,et al.Elastic wave finite-difference simulation using discontinuous curvilinear grid with non-uniform time step:two-dimensional case[J].Geophysical Journal International,2015,202(1):102-118

[12] 孫林潔,劉春成,張世鑫.PML邊界條件下孔隙介質(zhì)彈性波可變網(wǎng)格正演模擬方法研究[J].石油物探,2015,54(6):652-664

SUN L J,LIU C C,ZHANG S X.A variable grid finite-difference scheme with PML boundary condition in elastic wave of porous media[J].Geophysical Prospecting for Petroleum,2015,54(6):652-664

[13] 馬光克,李洋森,孫萬元,等.可變網(wǎng)格高階有限差分法逆時(shí)偏移研究[J].石油物探,2016,55(5):719-736

MA G K,LI Y S,SUN W Y,et al.Acoustic pre-stack reverse time migration using variable grid finite-difference method[J].Geophysical Prospecting for Petroleum,2016,55(5):719-736

[14] ALFORD R M,KEILY K R,BORE D M.Accuracy of finite difference modeling of the acoustic wave equation[J].Geophysics,1974,39(6):834-842

[15] 董良國,李培明.地震波傳播數(shù)值模擬中的頻散問題[J].天然氣工業(yè),2004,24(6):53-56

DONG L G,LI P M.Dispersion problem in seismic wave propagation numerical modeling[J].Natural Gas Industry,2004,24(6):53-56

[16] DU Q Z,LI B,HOU B.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered-grid finite difference scheme[J].Applied Geophysics,2009,6(1):42-49

[17] TAM C K W,WEBB J C.Dispersion-relation-preserving finite difference schemes for computational acoustics[J].Journal of Computational Physics,1993,107(2):262-281

[18] LIU Y,SEN M K.Scalar wave equation modeling with time-space domain dispersion-relation-based staggered-grid finite-difference schemes[J].Bulletin of the Seismological Society of America,2011,101(101):141-159

[19] 杜啟振,白清云,李賓.橫向各向同性介質(zhì)優(yōu)化差分系數(shù)法地震波場數(shù)值模擬[J].石油地球物理勘探,2010,45(2):170-176

DU Q Z,BAI Q Y,LI B.Optimized difference coefficient method seismic wave field numerical simulation in vertical transverse isotropic medium[J].Oil Geophysical Prospecting,2010,45(2):170-176

[20] 梁文全,王彥飛,楊長春.基于優(yōu)化方法的時(shí)間-空間域隱格式有限差分算子確定方法[J].石油物探,2015,54(3):254-259

LIANG W Q,WANG Y F,YANG C C.Determination on the implicit finite-difference operator based on optimization method in time-space domain[J].Geophysical Prospecting for Petroleum,2015,54(3):254-259

[21] YANG L,YAN H,LIU H.Least squares staggered-grid finite-difference for elastic wave modelling[J].Exploration Geophysics,2014,45(4):255-260

[22] YAN H Y,YANG L,LIU H.Acoustic reverse-time migration using optimal staggered-grid finite-difference operator based on least squares[J].Acta Geophysica,2015,63(3):715-734

[23] IGEL H,MORA P,RIOLLET B.Anisotropic wave propagation through finite-difference grids[J].Geophysics,1995,60(4):1203-1216

[24] 雍鵬,黃建平,李振春,等.一種時(shí)空域優(yōu)化的高精度交錯(cuò)網(wǎng)格差分算子與正演模擬[J].地球物理學(xué)報(bào),2016,59(11):4223-4233

YONG P,HUANG J P,LI Z C,et al.Optimized staggered-grid finite-difference method in time-space domain based on exact time evolution schemes[J].Chinese Journal of Geophysics,2016,59(11):4223-4233

[25] 梁文全,王彥飛,楊長春.聲波方程數(shù)值模擬矩形網(wǎng)格有限差分系數(shù)確定法[J].石油地球物理勘探,2017,52(1):56-62

LIANG W Q,WANG Y F,YANG C C.Acoustic wave equation modeling with rectangle grid finite difference operator and its linear time space domain solution[J].Oil Geophysical Prospecting,2017,52(1):56-62

[26] KINDELAN M,KAMEL A,SGUAZZERO P.On the construction and efficiency of staggered numerical differentiators for the wave equation[J].Geophysics,1990,55(1):107-110

[27] ZHANG J,YAO Z.Optimized explicit finite-difference schemes for spatial derivatives using maximum norm[J].Journal of Computational Physics,2013,250(10):511-526

[28] YANG L,YAN H,LIU H.Optimal implicit staggered-grid finite-difference schemes based on the sampling approximation method for seismic modeling[J].Geophysical Prospecting,2016,64(3):595-610

附錄A 可變交錯(cuò)網(wǎng)格過渡帶差分系數(shù)

表A-2 大小網(wǎng)格步長比l=5的可變交錯(cuò)網(wǎng)格差分系數(shù)(d=0.12π,u=0.4π)

表A-3 大小網(wǎng)格步長比l=9的可變交錯(cuò)網(wǎng)格差分系數(shù)(d=0.07π,u=0.4π)

表A-4 大小網(wǎng)格步長比l=15的可變交錯(cuò)網(wǎng)格差分系數(shù)(d=0.05π,u=0.4π)