周誠堯， 汪 勇， 桂志先， 于曉東

(油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長江大學(xué))，武漢 430100)

地震正演模擬研究是探索地震波在介質(zhì)中傳播規(guī)律的重要手段,基于波動(dòng)方程進(jìn)行的數(shù)值模擬方法是地震正演模擬研究中的一種方式[1]。在目前常用的方法中，有限差分解法是求解波動(dòng)方程的一種重要數(shù)值方法。隨著地震正演模擬的發(fā)展，人們對(duì)于數(shù)值結(jié)果的精度要求越來越高[2]，傳統(tǒng)的有限差分?jǐn)?shù)值模擬方法有著兩個(gè)較為明顯的缺點(diǎn)：①需要用來進(jìn)行數(shù)值模擬的網(wǎng)格點(diǎn)數(shù)較多；②傳統(tǒng)的顯式差分格式具有較大的數(shù)值頻散[3]。緊致差分格式比傳統(tǒng)的差分格式有更高的精度和分辨率，且具有較高的穩(wěn)定性[4]。因此，緊致差分格式為近年來應(yīng)用和研究得較多的有限差分方法。1989年，Dennis和Hundson針對(duì)Navier-Stokes方程提出了空間四階的緊致格式[5]，LeLe于1992年對(duì)pade格式進(jìn)行了研究與總結(jié)，提出了求解一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的對(duì)稱型緊致差分格式，該格式精度可達(dá)到譜方法的精度[6]。為進(jìn)一步提高精度加強(qiáng)計(jì)算效率，Chu等在2000年時(shí)提出了三點(diǎn)六階超緊致差分格式用于對(duì)流擴(kuò)散方程[7]。隨后針對(duì)不同的問題，人們提出了許多種不同的緊致差分格式。如林東等在Chu的三點(diǎn)六階超緊致差分的基礎(chǔ)上提出了組合型超緊致差分，用于對(duì)淺水方程的求解[8]，其格式精度最高可達(dá)到十階。所應(yīng)用的五點(diǎn)八階超緊致差分格式是1998年Krishnan提出的一族具有高光譜分辨率的高階緊致差分格式的一種[9]。在將緊致差分運(yùn)用于地震正演模擬方面，王書強(qiáng)等于2002年將LeLe的理論運(yùn)用于對(duì)彈性波的導(dǎo)數(shù)的空間離散上[10]。然而在討論地震波傳播理論時(shí)，絕大部分情況下都是把大地傳播介質(zhì)看做完全彈性介質(zhì)[11]。實(shí)際上地下介質(zhì)屬于具有黏滯性的黏彈性介質(zhì)，將地震波當(dāng)成彈性波來研究其傳播規(guī)律往往不能有效反映其特征[12]。

將通過五點(diǎn)八階超緊致有限差分格式(five point eight order combined compact difference scheme， CCD8)應(yīng)用到位移場(chǎng)空間偏導(dǎo)數(shù)的求取中，并對(duì)該格式進(jìn)行穩(wěn)定性、頻散及精度分析。最后將其應(yīng)用到黏彈介質(zhì)的聲波方程正演數(shù)值模擬中來驗(yàn)證該格式的可行性。

1 五點(diǎn)八階超緊致有限差分方法原理及應(yīng)用

關(guān)于五點(diǎn)超緊致差分格式的產(chǎn)生，早期的緊致差分格式是基于Hermite多項(xiàng)式構(gòu)造而來的。LeLe在1992年對(duì)Hermite公式進(jìn)行了擴(kuò)展，構(gòu)造了如下的緊致差分格式：

(1)

式(1)中：h為網(wǎng)格間距；a、b、c、α、β均為差分系數(shù)。

對(duì)上述LeLe差分格式進(jìn)行泰勒公式展開和求解差分系數(shù)，可以得到二階導(dǎo)數(shù)的七點(diǎn)八階(seven point eight order compact difference 8， CD8)緊致差分格式：

(2)

從式(2)可以看出普通緊致差分格式若需要其精度為八階的時(shí)候，需用到七個(gè)節(jié)點(diǎn)。

然而超緊致有限差分格式只需要五個(gè)節(jié)點(diǎn)就可以達(dá)到八階精度， 1998年Krishnan提出了如下緊致差分格式：

(3)

式(3)中：a1、a0、a2、b1、b0、b2、c1、c2、c0、c3、c4為常數(shù)，可以根據(jù)不同精度要求進(jìn)行確定。式(3)可以達(dá)到六階以及八階精度，相比一般的緊致差分格式所用到的網(wǎng)格點(diǎn)數(shù)更少。利用其五點(diǎn)八階的超緊致格式(CCD8)：

(4)

根據(jù)上述超緊致差分的原理，將超緊致差分方法運(yùn)用于求解黏滯聲波方程初值問題中來[13]。黏彈介質(zhì)中的二維聲波方程可以表示為

(5)

(6)

利用式(6)就可以進(jìn)行地震波場(chǎng)時(shí)間層的推進(jìn)。采用五點(diǎn)八階超緊致差分格式對(duì)式(6)中u和?u/?t對(duì)x、z的導(dǎo)數(shù)進(jìn)行離散。為方便計(jì)算，令p=?u/?t，則式(6)改寫為

(7)

式(7)中：un+1為u在n+1時(shí)刻位移；un-1為u在n-1時(shí)刻位移；un為u在n時(shí)刻位移，pn為p=?u/?t在n時(shí)刻位移。

超緊致差分格式的應(yīng)用關(guān)鍵點(diǎn)在于對(duì)其空間導(dǎo)數(shù)的計(jì)算。利用五點(diǎn)八階超緊致差分格式對(duì)u和p的導(dǎo)數(shù)進(jìn)行逼近，式(7)左側(cè)可以得到如下形式的系數(shù)矩陣：

(8)

系數(shù)矩陣B則為五點(diǎn)八階超緊致差分公式右側(cè)的系數(shù)：

(9)

記矩陣F與u分別為

(10)

式(10)中：u為大小為m×n的位移場(chǎng)矩陣，F(xiàn)為大小為2m×n的位移場(chǎng)空間一階和二階導(dǎo)數(shù)矩陣。

(11)

(12)

同理p的導(dǎo)數(shù)部分亦可由上述步驟為模式進(jìn)行修改后得出，式(7)變?yōu)槿缦滦问剑?/p>

(13)

2 頻散及穩(wěn)定性分析

關(guān)于x方向上的修正波數(shù)二次方如下:

cos(kxh)+2(-17/216+23/216)cos(2kxh)]-1-[2(-1/108+374/243-1 391/1 944)×

cos(2kxh)+4(-17/3 888+23/1 944)×

cos(kxh)cos(2kxh)][1-17/108-23/108+2(-1/6+17/36)cos(kxh)+2(-17/216+23/216)cos(2kxh)]-1

(14)

同理z方向上的修正波數(shù)二次方為

(15)

式中：kxh=khcosθ，kzh=khsinθ，θ即波的傳播方向與x軸的夾角。

(16)

根據(jù)上述方法所得到速度比，將CCD8與CD8得到的速度比進(jìn)行一個(gè)比較。在理想情況下，如果不存在數(shù)值頻散則速度比恒等于1。若偏離1越大，則說明該方法的數(shù)值頻散越嚴(yán)重。反之，則說明該方法能更好地壓制數(shù)值頻散。圖1顯示了以上幾種方法在不同的品質(zhì)因子Q與θ時(shí)候的速度比曲線。

圖1 α=0.25時(shí)兩種差分方法在不同角度及品質(zhì)因子下的速度比曲線Fig.1 Velocity ratio curves of the numerical simulation in different methods and different angle with α=0.25

如圖1所示，這里所用的courant數(shù)α=0.25，橫坐標(biāo)為波數(shù)與空間步長的乘積(φ)，縱坐標(biāo)為速度比(y)，其中y=1表示修正波速與理論波速一致。可以看出：①在不同的角度下，兩種方法的的偏移趨勢(shì)沒有明顯的變化，即兩種方法在不同角度下其頻散壓制效果均無有明顯改變。②從速度比曲線的信息得知CCD8在偏離1時(shí)所對(duì)應(yīng)的橫坐標(biāo)數(shù)值為1.05，CD8所對(duì)應(yīng)的橫坐標(biāo)數(shù)值為0.8，因此CCD8與CD8在保證壓制頻散的時(shí)候所需的最低網(wǎng)格點(diǎn)數(shù)分別是6.0和7.85。③不同品質(zhì)因子的條件下，從圖1得知品質(zhì)因子的改變對(duì)速度比曲線影響不大，但是品質(zhì)因子減小會(huì)影響頻率衰減的速率，從而使傳播過程中地震子波的波長加長。這一影響表現(xiàn)為在網(wǎng)格長度不變時(shí)，每個(gè)波長內(nèi)的采樣點(diǎn)數(shù)會(huì)增加，進(jìn)而有利于壓制數(shù)值頻散。然而品質(zhì)因子并不是越小越好，這就要考慮到差分格式的穩(wěn)定性問題。分析可得：CCD8對(duì)于頻散的壓制有較好的效果。

關(guān)于CCD8的穩(wěn)定性分析這里使用Fourier分析[19-20]將公式ui, j=ξneIh(kx i + kz j)代入位移場(chǎng)及位移場(chǎng)空間一階和二階導(dǎo)數(shù)，可得：

(17)

(18)

將式(18)代入式(4)得到：

(19)

對(duì)其進(jìn)行求解可以得到：

(20)

在得到所需要的振幅關(guān)系式后，對(duì)于聲波方程來說，就是求在放大因子滿足穩(wěn)定性下的條件，再與現(xiàn)有條件做對(duì)比，而聲波方程的放大因子是下一時(shí)刻的振幅與當(dāng)前時(shí)刻的振幅的比值[21]。將式(6)變形為

(21)

(22)

因此放大因子為

(23)

αCCD8=(vp)maxΔt/h<0.5

(24)

式(24)中：αCCD8代表CCD8在黏彈介質(zhì)中的最大courant數(shù)。

3 截?cái)嗾`差對(duì)比

截?cái)嗾`差決定了有限差分方法的精度，表1示出CCD8與CD8的二階導(dǎo)數(shù)的截?cái)嗾`差，其數(shù)值是

表1 二階導(dǎo)數(shù)不同方法截?cái)嗾`差Table 1 Second derivative truncation error by different methods

通過Taylor截?cái)嗾`差分析所得。通過表1得到以下結(jié)論： CCD8的截?cái)嗾`差高于CD8。

4 PML邊界條件

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

(25)

式(25)中：u1、u2、L和J是引入的中間變量，d(x)和d(z)是x和z方向的衰減系數(shù)，分別起到衰減這兩個(gè)方向波的作用。在模型試算中采用余弦型吸收衰減函數(shù)如下所示：

(26)

式(26)中：αi代表吸收衰減因子；λ為衰減幅度因子，設(shè)為500；PML為吸收邊界的層數(shù)，取為20。

5 模型試算

采用五點(diǎn)八階超緊致格式對(duì)一些典型的介質(zhì)進(jìn)行黏滯聲波方程數(shù)值模擬。

5.1 均勻介質(zhì)模型

同樣是具有八階精度的有限差分格式，相比CD8，CCD8不僅在完全彈性介質(zhì)的數(shù)值模擬結(jié)果對(duì)比中具有頻散低的特點(diǎn)，而且在黏彈介質(zhì)聲波方程數(shù)值正演模擬上亦具有良好的壓制頻散的特性。且觀察正演模擬過程中可知使用CCD8進(jìn)行正演模擬運(yùn)算速度要優(yōu)于CD8。因此可以在實(shí)際的正演數(shù)值模擬應(yīng)用中考慮使用CCD8，以減少運(yùn)算時(shí)間。

5.2 水平層狀介質(zhì)模型

模型參數(shù)如表2所示。

表2 多層層狀介質(zhì)模型參數(shù)

模型為四層層狀介質(zhì)，各層參數(shù)如表2所示，層狀介質(zhì)的長度和深度均為1 000 m，為了比較黏彈和完全彈性介質(zhì)中地震波場(chǎng)的區(qū)別，第一層和第四層設(shè)置的為完全彈性介質(zhì)，因此品質(zhì)因子取為∞，而第二、三層的品質(zhì)因子則取Q=50和Q=100。以地表x=500 m處進(jìn)行激發(fā)，所用的能量震源為30 Hz雷克子波震源，時(shí)間步長與空間步長分別取Δt=0.5 ms與Δx=Δz=10 m。地表接收的單炮記錄長度為1 s，如圖3所示，邊界條件為PML邊界條件。

圖2 0.2 s時(shí)刻波場(chǎng)快照 Fig.2 Snapshots at 0.2 s

由圖3可得知四層水平層狀介質(zhì)模型數(shù)值模擬得到的地震記錄非常清晰，沒有數(shù)值頻散和邊界反射，地震記錄中的直達(dá)波和反射波顯示清晰，因此證明了CCD8可以有效地模擬黏滯聲波在多層介質(zhì)模型中的傳播。

5.3 Marmousi模型

為了驗(yàn)證本文算法對(duì)復(fù)雜介質(zhì)的適用性，用經(jīng)典的二維Marmousi縱波速度模型進(jìn)行數(shù)值模擬，模型的長度和深度均為2 500 m，其他的參數(shù)如下，縱波速度范圍：1 729～5 500 m/s，品質(zhì)因子根據(jù)經(jīng)驗(yàn)公式Q=14×Vp2.2得到的范圍為46.7～595.6。以地表中心處(1 250 m，0 m)激發(fā)，所用的能量震源為30 Hz雷克子波震源。時(shí)間步長與空間步長分別取Δt=0.2 ms與Δx=Δz=10 m，采樣時(shí)間為2 s。

圖3 CCD8水平層狀模型的黏彈介質(zhì)聲波方程地面地震記錄Fig.3 CCD8 is used to the ground seismic record of the acoustic wave equation of the viscoelastic medium of the horizontal layered dielectric model

圖4為Marmousi模型采用CCD8在黏滯聲波方程中得到的地面地震記錄。

圖4 CCD8 marmousi模型的黏彈介質(zhì)聲波方程地面地震記錄Fig.4 CCD8 is used to the ground seismic record of the acoustic wave equation of the viscoelastic medium of the marmousi model

由圖4可知使用CCD8進(jìn)行正演數(shù)值模擬所得到的地面地震記錄平滑清晰，邊界吸收效果較好，無明顯邊界反射且該地震記錄無數(shù)值頻散。從而證明了CCD8可以應(yīng)用于復(fù)雜模型的數(shù)值正演模擬。

6 結(jié)論

重點(diǎn)研究了五點(diǎn)八階超緊致有限差分格式(CCD8)并將此用于黏滯聲波方程的數(shù)值模擬中。通過研究與分析，得到如下幾點(diǎn)認(rèn)識(shí)。

(1)與七點(diǎn)八階緊致有限差分格式(CD8)相比，CCD8僅需使用五個(gè)節(jié)點(diǎn)，就能使二階空間偏導(dǎo)數(shù)的差分精度達(dá)到八階，且證明了CCD8具有低數(shù)值頻散，較高穩(wěn)定性和模擬精度的特點(diǎn)，適用于大尺度和長時(shí)間的地震波場(chǎng)數(shù)值模擬的特點(diǎn)。

(2)模型試算的結(jié)果說明，CCD8不僅可以用于均勻介質(zhì)的黏滯聲波方程的波場(chǎng)模擬，還可以進(jìn)一步推廣到二維或三維的各向異性介質(zhì)以及復(fù)雜介質(zhì)的數(shù)值模擬中。

(3)在進(jìn)行黏滯聲波方程數(shù)值模擬時(shí)使用的是均勻網(wǎng)格，這對(duì)于地震正演模擬來說還是不夠，今后的研究中不僅可以研究建立一個(gè)穩(wěn)定性更高的差分格式，也應(yīng)考慮與其他非均勻的網(wǎng)格相結(jié)合，例如交錯(cuò)網(wǎng)格及變步長網(wǎng)格。