周誠堯,汪 勇,蔡偉祥,桂志先

(1.油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長江大學(xué)),湖北武漢430100;2.中石化重慶涪陵頁巖氣勘探開發(fā)有限公司,重慶408014)

地震正演模擬是探索地震波在介質(zhì)中傳播規(guī)律的重要手段,基于波動(dòng)方程的數(shù)值模擬方法是一種地震正演模擬方法[1]。目前波動(dòng)方程法主要包括偽譜法[2]、有限元法[3]、邊界元法[4]、譜元法[5]以及有限差分法。有限元法和有限差分方法相比,二者精度相同時(shí),有限元法耗時(shí)長[6]。有限差分法因其具有簡單靈活、計(jì)算效率高以及占用內(nèi)存小等優(yōu)勢被廣泛應(yīng)用于地震波場數(shù)值模擬[7]。近年來，隨著人們對數(shù)值模擬結(jié)果精度的要求不斷提升，傳統(tǒng)有限差分?jǐn)?shù)值模擬方法的不足逐步凸顯:一是數(shù)值計(jì)算時(shí)所需的網(wǎng)格點(diǎn)數(shù)多;二是具有嚴(yán)重的數(shù)值頻散[8]。相較傳統(tǒng)的差分格式,緊致差分格式在相同的計(jì)算網(wǎng)格時(shí),精度和分辨率高,且穩(wěn)定性高[9-10]。緊致差分格式的發(fā)展可以追溯到1989年DENNIS等針對Navier-Stokes方程首次提出了空間四階緊致格式[11];1992年LELE對Pade格式進(jìn)行了研究,提出了求解一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的對稱型緊致差分格式[12],該格式精度可達(dá)到譜方法的精度;1998年CHU等提出了三點(diǎn)六階超緊致差分(combined compact difference 6,CCD6)格式用于對流擴(kuò)散方程[13];林東等在三點(diǎn)CCD6格式的基礎(chǔ)上提出了組合型超緊致差分格式,其精度最高可達(dá)到十階[14]。在緊致差分格式中引入迎風(fēng)機(jī)制可有效抑制非物理振蕩,理論上可以提高緊致差分格式的穩(wěn)定性。1994年傅德薰在緊致差分格式中引入迎風(fēng)機(jī)制[15];1997年FU等又提出了更適合于多尺度復(fù)雜流場計(jì)算的五點(diǎn)五階精度的迎風(fēng)緊致格式[16];2002年王書強(qiáng)等在地震波場數(shù)值模擬時(shí),首次將常規(guī)緊致差分格式應(yīng)用于彈性波方程的數(shù)值模擬[17],截至目前，超緊致差分格式的地震波場數(shù)值模擬應(yīng)用未見文獻(xiàn)報(bào)道。

本文首先討論了五點(diǎn)八階超緊致差分(combined compact difference 8,CCD8)格式,然后通過引入迎風(fēng)機(jī)制,在五點(diǎn)八階超緊致差分格式的基礎(chǔ)上,構(gòu)造了五點(diǎn)七階迎風(fēng)超緊致差分(upwind combined compact difference 7,UCCD7)格式,并將這兩種格式應(yīng)用于位移場空間導(dǎo)數(shù)的求取,實(shí)現(xiàn)了二維均勻介質(zhì)模型、水平層狀介質(zhì)模型以及Marmousi模型的地震波場數(shù)值模擬。

1 超緊致差分格式及其迎風(fēng)型方法原理

早期的緊致差分格式是基于Hermite多項(xiàng)式構(gòu)造形成的,1992年LELE對Hermite公式進(jìn)行了擴(kuò)展,構(gòu)造了緊致差分格式[12]:

(1)

對上述差分格式進(jìn)行泰勒展開并求解差分系數(shù),可以得到一、二階導(dǎo)數(shù)的普通八階的緊致差分(compact difference 8,CD8)格式。

一階導(dǎo)數(shù)普通八階緊致差分格式:

(2)

二階導(dǎo)數(shù)普通八階緊致差分格式:

(3)

從公式(2)和公式(3)可知,普通緊致差分格式若要達(dá)到八階精度,需用到7個(gè)節(jié)點(diǎn)。此外,緊致差分格式也適用于交錯(cuò)網(wǎng)格,其中八階緊致交錯(cuò)網(wǎng)格差分(compact staggered difference 8,CSD8)格式可表示為:

(4)

式中:a1,b1,b2,b3是CSD8格式的差分系數(shù);xi,Δx分別為x方向上的第i個(gè)節(jié)點(diǎn)、x方向上的空間步長。為了達(dá)到2n+2階精度,只需2n個(gè)節(jié)點(diǎn)函數(shù)值。

超緊致有限差分格式在緊致差分格式基礎(chǔ)上發(fā)展而來,1998年KRISHNAN提出了五點(diǎn)CCD8格式[18]:

(5)

將迎風(fēng)機(jī)制引入五點(diǎn)CCD8格式,可得到如下的五點(diǎn)七階迎風(fēng)超緊致差分(UCCD7)格式,其精度可達(dá)七階[19]:

(6)

根據(jù)超緊致差分的原理,將超緊致差分方法用于求解聲波方程初值問題[20]。在完全彈性介質(zhì)中聲波方程可表示為:

(7)

式中:u(x,z)為二維介質(zhì)的質(zhì)點(diǎn)在x方向、z方向的位移;?2u/?t2為u(x,z)對時(shí)間t的二階導(dǎo)數(shù);?2u/?x2,?2u/?z2分別為u(x,z)在x方向和z方向的二階導(dǎo)數(shù);v為聲波速度。

對(7)式采用四階中心差分對時(shí)間導(dǎo)數(shù)離散可得:

(8)

式中:un+1,un,un-1分別為n+1,n和n-1時(shí)刻位移;Δt為時(shí)間步長。

以z方向?yàn)槔?說明利用五點(diǎn)CCD8格式求解公式(8)中u對z方向的空間偏導(dǎo)數(shù)的過程,x方向的求解過程不再贅述。假設(shè)U為位移場,其大小為m×n,則求?u/?z和?2u/?z2的過程可以用矩陣表示為AF=BU。其中A為公式(5)左端的差分系數(shù)矩陣,大小為2m×2m;B為公式(5)右端的差分系數(shù)矩陣,大小為2m×m;F為待求位移場空間一階和二階導(dǎo)數(shù)矩陣,大小為2m×n。矩陣分別為:

B=

(11)

根據(jù)AF=BU,可得F=A-1BU,其中F的奇數(shù)行矩陣元素為?u/?z,偶數(shù)行矩陣元素為?2u/?z2。同理,將位移場轉(zhuǎn)置,可求得?u/?x和?2u/?x2。利用UCCD7格式求公式(8)中的空間偏導(dǎo)數(shù)的方法與此相同,僅有差分系數(shù)矩陣不同,故不再贅述。

以聲波方程為例,說明了五點(diǎn)CCD8格式或五點(diǎn)UCCD7格式在地震聲波波場數(shù)值模擬中的應(yīng)用方法,該方法同樣可以用于彈性波方程的數(shù)值模擬。各向同性完全彈性介質(zhì)中的二維彈性波方程可以表示為:

(12)

式中:u(x,z)和w(x,z)分別表示x方向和z方向的質(zhì)點(diǎn)位移;vP和vS分別表示縱、橫波速度。u和w的時(shí)間四階精度差分格式為[21]:

(13)

(14)

由公式(13)和公式(14)可以看出,彈性波方程與聲波方程的差分格式相比,增加了許多二階和四階混合偏導(dǎo)數(shù)。對于這些混合偏導(dǎo)數(shù),均可以按不同方向進(jìn)行多次導(dǎo)數(shù)求取,求取順序?qū)Y(jié)果無影響,求取方法不再贅述。

2 數(shù)值頻散、穩(wěn)定性及精度分析

對波動(dòng)方程的時(shí)間和空間偏導(dǎo)數(shù)的離散化造成了數(shù)值頻散,它使相速度變成了網(wǎng)格間距的函數(shù)。數(shù)值模擬時(shí),如果空間網(wǎng)格長度過大,會(huì)造成較大的求解誤差,產(chǎn)生數(shù)值頻散,降低模擬精度[22-23]。壓制數(shù)值頻散是采用有限差分方法時(shí)面臨的關(guān)鍵問題之一[24]。頻散關(guān)系分析既是判斷數(shù)值模擬所用方法優(yōu)劣的重要依據(jù),也是確定空間網(wǎng)格大小的重要依據(jù)[25],五點(diǎn)CCD8的二階導(dǎo)數(shù)數(shù)值波數(shù)可以由公式(15)表示:

(15)

(16)

式中:v′為數(shù)值波速;v是真波速;vΔt/h稱為courant數(shù)。

數(shù)值波速與真波速的比值γ可表示為:

(17)

(18)

將多種有限差分方法的速度比進(jìn)行比較,結(jié)果如圖1所示。圖1a,圖1b和圖1c展示了不同的θ條件下4種常見網(wǎng)格差分格式的速度比曲線,其中α=vΔt/h;圖1d中增加了八階緊致交錯(cuò)網(wǎng)格格式的速度比曲線;圖1e增加了五點(diǎn)七階迎風(fēng)超緊致差分格式的速度比曲線。考慮到參與對比的方法的穩(wěn)定性,此處courant數(shù)均為0.25,橫坐標(biāo)φ∈[0,π],為波數(shù)與空間步長的乘積,縱坐標(biāo)γ為速度比。速度比為1表示數(shù)值波速與理論波速一致,說明該方法能很好地壓制數(shù)值頻散,反之,則說明該方法的數(shù)值頻散嚴(yán)重。

由圖1可以看出,五點(diǎn)CCD8格式的速度比曲線偏離1時(shí),其橫坐標(biāo)的數(shù)值最大,所以五點(diǎn)CCD8格式壓制數(shù)值頻散的能力強(qiáng)。圖1e說明將五點(diǎn)CCD8格式引入迎風(fēng)機(jī)制后的五點(diǎn)UCCD7格式,壓制頻散效果并無明顯降低。五點(diǎn)UCCD7,五點(diǎn)CCD8,三點(diǎn)CCD6,CD6,CD8和CSD8格式在偏離1時(shí)所對應(yīng)的橫坐標(biāo)數(shù)值分別為1.57,1.69,1.12,0.98,1.55和1.61,這6種方法在保證無數(shù)值頻散時(shí)所需的最低網(wǎng)格點(diǎn)數(shù)分別為4.00,3.70,5.60,6.40,4.05和3.90個(gè)。五點(diǎn)CCD8格式的數(shù)值頻散的壓制效果最好,CD6格式的數(shù)值頻散的壓制效果最差。

穩(wěn)定性條件是影響差分方法計(jì)算效率的重要因素。本文使用Fourier方法[26-27]進(jìn)行穩(wěn)定性分析,略去推導(dǎo)過程,直接給出聲波方程時(shí)間四階精度的五點(diǎn)CCD8和五點(diǎn)UCCD7格式的穩(wěn)定性條件αCCD8和αUCCD7分別為:

(19)

式中:vmax是指正演模擬時(shí),所建立模型的最大聲波速度,與五點(diǎn)CCD8格式相比,引入迎風(fēng)機(jī)制的五點(diǎn)UCCD7格式的穩(wěn)定性有了一定的提高,可使用略大的時(shí)間網(wǎng)格,在一定程度上減少計(jì)算時(shí)間,提高計(jì)算效率。

截?cái)嗾`差決定了有限差分格式和地震波場數(shù)值模擬精度,將文中所提到的幾種方法的二階導(dǎo)數(shù)差分格式截?cái)嗾`差進(jìn)行對比。由Taylor級數(shù)展開[28]可以計(jì)算得到CD6、三點(diǎn)CCD6、CD8和五點(diǎn)CCD8格式的截?cái)嗾`差分別為-4.2163×10-4,-4.9603×10-5,2.7000×10-5和2.2046×10-6,即五點(diǎn)CCD8格式的截?cái)嗾`差最小,數(shù)值模擬精度最大。

圖1 不同差分方法的速度比曲線a θ=0,α=0.25時(shí)4種常見差分格式; b θ=π/6,α=0.25時(shí)4種常見差分格式; c θ=π/3,α=0.25時(shí)4種常見差分格式; d θ=0,α=0.25時(shí)八階緊致交錯(cuò)網(wǎng)格格式與4種常見差分格式; e θ=π/3,α=0.25時(shí)五點(diǎn)七階迎風(fēng)超緊致差分格式與4種常見差分格式

3 PML邊界條件

數(shù)值模擬時(shí),由于模型大小的限制,必然會(huì)存在人工邊界,而人工邊界的處理直接影響到數(shù)值模擬的精度及計(jì)算效率,若不對其進(jìn)行處理則會(huì)干擾正常的地震波場。因此,本文在模型試算時(shí),采用完全匹配層(perfectly matched layer,簡稱PML)對人工邊界進(jìn)行處理,其基本思想是:在吸收邊界區(qū)域匹配與計(jì)算區(qū)域相同的波阻抗,使入射波無反射的進(jìn)入吸收層,吸收層為衰減介質(zhì),可使得入射波迅速衰減直至消除[29]。理論上,PML邊界能夠吸收任何入射角度和頻率的入射波,實(shí)踐也證明PML吸收邊界條件應(yīng)用效果良好,可應(yīng)用于聲波和彈性波的研究[30-31],下式為聲波方程的PML邊界條件的控制方程:

(20)

式中:u1,u2,J和L均為中間變量;d(x),d(z)分別為x,z方向的衰減系數(shù),用于衰減這兩個(gè)方向的波場。由文獻(xiàn)[32]可知,由于d(x)和d(z)的衰減函數(shù)和層數(shù)影響了邊界反射的衰減效果,因此余弦型吸收衰減函數(shù)效果較好,所以本文在模型試算中采用如下的余弦型吸收衰減函數(shù):

(21)

i=0,1,2,…,PML

式中:αi代表吸收衰減因子;λ為衰減幅度因子,λ=

500;PML為吸收邊界的層數(shù),PML=20。

4 模型試算

4.1 聲波方程均勻介質(zhì)模型波場模擬

均勻介質(zhì)模型的長度及深度均為4800m,聲波速度為4000m/s,激發(fā)位置位于模型中心,震源為20Hz雷克子波,設(shè)置縱、橫向空間步長均為48m,時(shí)間步長為1ms。

圖2為400ms時(shí)刻不同模擬方法得到的波場快照。

圖2 400ms時(shí)刻不同方法得到的波場快照a 普通六階緊致差分(CD6)格式波場快照; b 普通八階緊致差分(CD8)格式波場快照; c 三點(diǎn)六階超緊致差分(CCD6)格式波場快照; d 八階緊致交錯(cuò)網(wǎng)格差分(CSD8)格式波場快照; e 五點(diǎn)八階超緊致差分(CCD8)格式波場快照; f 迎風(fēng)五點(diǎn)七階超緊致差分(UCCD7)格式波場快照

可以看出,采用CD6和CD8格式模擬得到的波場快照存在嚴(yán)重的頻散;采用三點(diǎn)CCD6格式得到的結(jié)果仍有輕微的頻散;采用CSD8、五點(diǎn)CCD8和五點(diǎn)UCCD7格式模擬得到的波場清晰,無數(shù)值頻散。

4.2 彈性波方程均勻介質(zhì)模型波場模擬

采用五點(diǎn)CCD8格式及五點(diǎn)UCCD7格式進(jìn)行彈性波方程數(shù)值模擬。模型長度和深度均為4000m,縱橫向空間步長均為10m,介質(zhì)為均勻介質(zhì),縱波速度4000m/s,橫波速度2310m/s,模型中心處激發(fā),震源為20Hz雷克子波。參考五點(diǎn)CCD8格式與五點(diǎn)UCCD7格式的最大courant數(shù),設(shè)置時(shí)間步長為1ms。圖3為分別采用兩種格式在400ms時(shí)刻的波場快照,其中圖3a和圖3b為采用五點(diǎn)CCD8格式所得到的結(jié)果,圖3c和圖3d為采用五點(diǎn)UCCD7格式得到的結(jié)果。波場均清晰,無數(shù)值頻散,縱波和橫波位移在水平和垂向記錄上的特征明顯,說明采用五點(diǎn)CCD8格式與五點(diǎn)UCCD7格式均可有效而清晰地模擬彈性波在各向同性介質(zhì)模型中的傳播。

圖3 均勻介質(zhì)模型彈性波數(shù)值模擬400ms時(shí)刻的波場快照a 采用五點(diǎn)CCD8格式得到的x方向位移分量; b 采用五點(diǎn)CCD8格式得到的z方向位移分量; c 采用五點(diǎn)UCCD7格式得到的x方向位移分量; d 采用五點(diǎn)UCCD7格式得到的z方向位移分量

4.3 聲波方程水平層狀均勻介質(zhì)模型數(shù)值模擬

水平層狀均勻介質(zhì)模型包括兩層,各層的厚度均為1000m,長度與深度均為2000m。第1層和第2層的地震波速分別為4000m/s和4500m/s。在地表x=1000m處激發(fā),震源為20Hz雷克子波,時(shí)間步長Δt=0.5ms,空間步長Δx=Δz=20m,地表接收的單炮記錄長度為1000ms。圖4為采用五點(diǎn)CCD8及UCCD7格式對水平層狀介質(zhì)模型進(jìn)行聲波方程數(shù)值模擬得到的結(jié)果,地震記錄清晰,無數(shù)值頻散和不穩(wěn)定數(shù)值結(jié)果,地震記錄中的直達(dá)波和反射波顯示清楚,說明采用五點(diǎn)CCD8格式與五點(diǎn)UCCD7格式可清晰有效地模擬聲波在多層各向同性介質(zhì)模型中的傳播。

圖4 均勻水平層狀介質(zhì)模型聲波方程數(shù)值模擬地面地震記錄a 五點(diǎn)CCD8格式; b 五點(diǎn)UCCD7格式

4.4 聲波方程Marmousi模型數(shù)值模擬

基于二維Marmousi縱波速度模型進(jìn)行數(shù)值模擬。在地表中心處(1250m,0)激發(fā)震源為30Hz雷克子波,時(shí)間步長與空間步長分別為Δt=0.2ms,Δx=Δz=10m,采樣時(shí)間為2000ms。圖5為對Marmousi模型采用五點(diǎn)CCD8格式及五點(diǎn)UCCD7格式進(jìn)行聲波方程模擬得到的地面地震記錄,地震記錄平滑清晰,邊界吸收效果好,無明顯邊界反射和數(shù)值頻散,驗(yàn)證了五點(diǎn)CCD8格式和五點(diǎn)UCCD7格式對復(fù)雜模型的適用性。

圖5 Marmousi模型的完全彈性介質(zhì)聲波方程地面地震記錄a 五點(diǎn)CCD8格式; b 五點(diǎn)UCCD7格式

5 結(jié)論與認(rèn)識(shí)

本文重點(diǎn)研究了五點(diǎn)八階超緊致有限差分(CCD8)和五點(diǎn)七階迎風(fēng)超緊致有限差分(UCCD7)兩種格式,并將其應(yīng)用于聲波和彈性波方程的數(shù)值模擬,得到3點(diǎn)認(rèn)識(shí):

1) 與傳統(tǒng)的有限差分格式相比,五點(diǎn)CCD8格式僅需5個(gè)節(jié)點(diǎn)就能使二階空間偏導(dǎo)數(shù)的差分精度達(dá)到八階,該格式具有低數(shù)值頻散、高模擬精度及高穩(wěn)定性的特點(diǎn);

2) 引入迎風(fēng)機(jī)制,優(yōu)化五點(diǎn)CCD8格式,得到了五點(diǎn)UCCD7格式,提高了聲波方程差分格式的穩(wěn)定性,同時(shí)保持了較好的低數(shù)值頻散的優(yōu)勢;

3) 模型試算結(jié)果說明,五點(diǎn)CCD8格式與五點(diǎn)UCCD7格式不僅適用于聲波方程的波場模擬,還可以進(jìn)一步推廣應(yīng)用于二維或三維的各向異性介質(zhì)和粘彈介質(zhì)等復(fù)雜介質(zhì)的彈性波方程數(shù)值模擬。

本文迭代求取了空間四階導(dǎo)數(shù),這導(dǎo)致了四階導(dǎo)數(shù)的模擬精度降低。在今后的研究中,可建立一種新的組合型緊致差分格式,一次性求取四階導(dǎo)數(shù)及混合偏導(dǎo)數(shù),以提高數(shù)值模擬的精度。