汪 勇 徐佑德 高 剛 桂志先 陳 英 王亞楠

(①油氣資源與勘探技術(shù)教育部重點實驗室(長江大學(xué))，湖北武漢 430100； ②長江大學(xué)地球物理與石油資源學(xué)院，湖北武漢 430100； ③中國石化勝利油田分公司勘探開發(fā)研究院西部分院，山東東營 257001；④東方地球物理公司裝備處測量中心，河北涿州 072751； ⑤東方地球物理公司研究院地質(zhì)研究中心，河北涿州 072751)

1 引言

目前，常用的地震波場數(shù)值模擬方法主要有射線追蹤法和波動方程法，其中波動方程法有偽譜法、有限元法、邊界元法、譜元法和有限差分法[1-5]。隨著數(shù)值模擬技術(shù)的發(fā)展和生產(chǎn)實踐的要求，圍繞著提高有限差分計算效率[6]、模擬精度[7]、算法穩(wěn)定性[8-9]、處理復(fù)雜介質(zhì)[10-11]、吸收邊界條件[12]和壓制數(shù)值頻散[13-15]等方面，已經(jīng)提出了許多方法,取得了大量有意義的研究成果。

提高差分格式數(shù)值計算精度最直接的方法就是在差分計算時增加網(wǎng)格節(jié)點數(shù)，但會大大增加計算量和存儲空間。緊致差分方法恰好能夠較好地解決這個矛盾，同時緊致差分是一種隱式差分格式，具有較好的穩(wěn)定性，這些優(yōu)勢也使之成為目前研究較多的有限差分方法之一。Dennis等[16]針對Navier-Stokes方程提出了空間四階的緊致差分格式，能夠使用粗網(wǎng)格獲得較高的精度；Lele[17]構(gòu)造了求解一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的緊致差分格式；Adams等[18]提出了緊致非震蕩(Essentially Non-oscillation，ENO)差分格式，用于求解激波湍流相互作用問題；Chu等[19]提出了三點組合型緊致差分格式，并將其用于求解對流擴(kuò)散方程。也有不少學(xué)者將緊致差分格式用于聲波、彈性波和復(fù)雜介質(zhì)等地震波場數(shù)值模擬，取得了許多研究成果[20-24]。

地層的黏彈性對地震波產(chǎn)生的吸收和衰減規(guī)律非常復(fù)雜，所以研究地震波在黏彈性介質(zhì)中的傳播規(guī)律對地震勘探有著重要的意義。李曉波等[25]模擬了地震波在斑狀飽和介質(zhì)中的傳播，分析了斑狀飽和介質(zhì)中孔隙度及含氣飽和度的變化對地震波傳播特征的影響；姚振岸等[26]進(jìn)行了黏彈各向異性介質(zhì)中的微地震波場模擬，并分析了微地震信號的傳播特性和偏振特性；汪勇等[27]利用近似解析離散化算子計算空間導(dǎo)數(shù)，對黏滯聲波方程進(jìn)行了數(shù)值模擬；馬靈偉等[28]、羅文山等[29]、姚振岸等[30]基于不同方法對不同的黏彈介質(zhì)模型進(jìn)行了數(shù)值模擬。

組合型緊致差分格式在黏滯聲波方程地震波場模擬中的應(yīng)用，尚未見到文獻(xiàn)報道。本文在黏滯聲波方程時間二階離散格式的基礎(chǔ)上，將組合型緊致格式應(yīng)用于位移場空間導(dǎo)數(shù)的求取，實現(xiàn)了二維黏滯聲波方程地震波場的數(shù)值模擬。

2 組合型緊致有限差分方法原理

早期的緊致差分格式是基于Hermite多項式構(gòu)造而來。Lele[17]對Hermite公式進(jìn)行了擴(kuò)展，構(gòu)造了求解一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的緊致差分格式(Compact Finite Difference，CD)，其特點是用相鄰節(jié)點的函數(shù)值和導(dǎo)數(shù)值計算待求節(jié)點上的導(dǎo)數(shù)值。而常規(guī)中心差分(Central Finite Difference，F(xiàn)D)僅利用函數(shù)值求解中心節(jié)點的導(dǎo)數(shù)值。一般情況下，在相同網(wǎng)格間距時，CD格式比常規(guī)FD格式具有更高的精度和更小的數(shù)值頻散[20]。

Chu[19]等構(gòu)造了精度更高的三點六階組合型緊致差分(Combined Compact Difference，CCD)格式

(1)

式中：f為一維函數(shù)；h為空間采樣間隔。上式CCD格式只需要相鄰的三個節(jié)點就可以同時求得一階和二階導(dǎo)數(shù)的六階精度近似值，比常規(guī)CD格式的節(jié)點數(shù)更少。并且式(1)中的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)是耦合的，既可以同時求出，又利于波形保真。

3 黏滯聲波方程組合型緊致有限差分法分析

經(jīng)常使用Kelvin-Voigt體模型描述聲波在非完全彈性介質(zhì)中的傳播，認(rèn)為地震波吸收系數(shù)與地震信號頻率的平方成正比，這與實際情況較為吻合。該模型假定介質(zhì)的應(yīng)力包括兩部分，一部分是彈性應(yīng)力，另一部分是黏滯應(yīng)力，其中彈性應(yīng)力與應(yīng)變成正比，黏滯應(yīng)力與應(yīng)變的時間變化率成正比。該模型的二維黏滯聲波方程可以表示為

(2)

利用時間二階導(dǎo)數(shù)的二階精度近似式

(3)

代入式(2)，可以得到位移場三層顯式差分格式

(4)

(5)

(6)

式中：A和C為式(1)左端的差分系數(shù)矩陣，階數(shù)分別為2M×2M和2N×2N；E和F為待求位移場空間一階和二階導(dǎo)數(shù)矩陣，階數(shù)分別為2M×N和M×2N；B和D為式(1)右端的差分系數(shù)矩陣，階數(shù)分別為2M×M和N×2N；U為位移場矩陣，階數(shù)為M×N。這些矩陣分別為

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

3.1 精度分析

不論是利用CCD格式，還是利用常規(guī)的七點六階FD格式或五點六階CD格式[17]，進(jìn)行黏滯聲波方程數(shù)值模擬時，它們在時間層推進(jìn)方式上是相同的，即都是利用差分格式式(5)，時間差分為二階精度。CD、FD和CCD不同之處在于它們分別利用不同的差分格式計算位移的空間導(dǎo)數(shù)。這三種格式計算的二階導(dǎo)數(shù)的截斷誤差主項系數(shù)分別為-17/(8!)、72/(8!)、-2/(8!)。雖然三種方法都能達(dá)到空間6階精度，但截斷誤差有較大的差別，F(xiàn)D和CD格式計算二階偏導(dǎo)數(shù)的截斷誤差約是CCD的36倍和8.5倍，說明CCD方法具有更小的截斷誤差和更高的差分精度。

應(yīng)用一維黏彈介質(zhì)模型比較CCD、CD和FD三種格式的數(shù)值模擬精度，其中模型長度為6000m、聲波速度為4000m/s、品質(zhì)因子Q=10。數(shù)值模擬中，Δt=1ms，h=40m。在模型的3000m處加載頻率為10Hz簡諧波震源，可得700ms時刻的模擬波形(圖1)。式(2)的一維黏滯聲波方程的解析解見文獻(xiàn)[27]。從局部放大圖中可以看出，CCD格式的數(shù)值模擬結(jié)果與精確解析解最接近，其模擬精度最高，CD和FD格式次之，這是由于在計算空間導(dǎo)數(shù)時截斷誤差不同。

3.2 頻散分析

頻散分析既是判斷數(shù)值模擬方法優(yōu)劣的重要指標(biāo)，也是確定空間網(wǎng)格大小的重要依據(jù)。通過對黏滯聲波方程的組合型緊致差分格式進(jìn)行數(shù)值頻散分析，進(jìn)一步分析方法的適用條件。

首先考慮x方向，令

(15)

圖1 三種差分格式的一維黏滯聲波方程數(shù)值模擬波形(700ms)(下)及其局部放大顯示(上)

代入差分格式式(1)，解方程組可得

(16)

同理可得z方向的解

(17)

將式(16)和式(17)代入黏滯聲波方程的差分格式(5)中，并令c0=vΔt/h為庫朗數(shù)，χ=kv′Δt(其中v′為數(shù)值波速，v是真波速)，利用歐拉公式，化簡可得CCD格式頻散關(guān)系

(18)

上式表明，黏滯聲波方程的頻散關(guān)系不僅與空間間距的取值有關(guān)，而且與介質(zhì)的品質(zhì)因子有關(guān)。通過上述頻散關(guān)系，確定φ后可以解得對應(yīng)的χ，且定義數(shù)值波速與真實速度的比值為

(19)

在理想情況下，如果不存在數(shù)值頻散則速度比γ恒等于1。γ偏離1越大，則說明該方法的數(shù)值頻散越嚴(yán)重，反之則說明該方法能更好地壓制數(shù)值頻散。取θ=π/4，計算CCD格式的速度比與φ的關(guān)系曲線，如圖2所示。

圖2 黏滯聲波方程CCD格式速度比曲線(a)Q=100，c0不同; (b)Q不同，c0 =0.2

取φ∈[0,π]作為橫坐標(biāo)，它是波數(shù)與空間間距的乘積，單位波長內(nèi)采樣點數(shù)N=2πφ，所以橫坐標(biāo)也可以看作N由∞逐漸減小至2。從圖中的速度比曲線可以看出：①隨著空間采樣點數(shù)的減少，速度比曲線逐漸偏離1，頻散現(xiàn)象逐步加劇。CCD方法的空間網(wǎng)格長度過大時，速度比曲線上翹，數(shù)值速度大于真速度。②當(dāng)庫朗數(shù)增加時，頻散加劇，說明即使空間網(wǎng)格間距選取合適時，加大時間采樣間距也會增加數(shù)值頻散。③品質(zhì)因子對速度比曲線影響不大，當(dāng)Q>5時，速度比曲線基本重合，具有相同的頻散特征。但這不能說明品質(zhì)因子大小與數(shù)值頻散無關(guān)，因為品質(zhì)因子影響著頻率衰減的快慢，進(jìn)而影響了傳播過程中地震子波的波長，在同樣網(wǎng)格間距情況下，每個波長內(nèi)的采樣點也在隨之改變，后文用波場快照做進(jìn)一步分析。④在不考慮Q值大小的情況下，當(dāng)c0分別為0.2、0.3和0.4時，為保證無數(shù)值頻散，每個波長內(nèi)需要6.7、10.0和13.3個樣點。

利用二維均勻模型驗證以上分析的結(jié)論。設(shè)置模型長度和深度均為4000m，聲波速度為4000m/s，震源子波為sin(2πf0t)exp(-π2f0t2/4)，其主頻為40Hz，Δt=0.75ms。圖3為CCD和常規(guī)FD格式不同品質(zhì)因子和網(wǎng)格間距時的波場快照，可以看出：①圖3a中品質(zhì)因子為2000，近似完全彈性介質(zhì)，傳播過程中近似無衰減。根據(jù)頻散分析結(jié)論，每個波長需采樣6.7個點，取網(wǎng)格間距為15m，340ms波場快照無明顯數(shù)值頻散。②當(dāng)網(wǎng)格間距為20m時，采樣略微不足，圖3b在對角線方向上頻散不明顯，而在x、z軸方向出現(xiàn)數(shù)值頻散，數(shù)值波速大于真實速度。③圖3c中的品質(zhì)因子僅為20，波場清楚，無數(shù)值頻散。這是由于地震波在黏彈介質(zhì)中傳播時，隨傳播距離的增加頻率降低衰減、波長增加，有利于壓制數(shù)值頻散。④由于地震子波頻率衰減，波長增加，圖3f所示結(jié)果平滑，無明顯數(shù)值頻散，而圖3d和圖3e均存在較嚴(yán)重的數(shù)值頻散，說明CCD格式比常規(guī)五點六階FD格式具有低數(shù)值頻散的優(yōu)勢，能夠適用于粗網(wǎng)格和大尺度模型的地震波場數(shù)值模擬。

圖3 均勻模型CCD和常規(guī)FD格式不同品質(zhì)因子和網(wǎng)格間距時的340ms時刻波場快照(a)CCD，Q=2000，Δx=15m; (b)CCD，Q=2000，Δx=20m; (c)CCD，Q=20，Δx=20m;(d)FD，Q=2000，Δx=15m; (e)FD，Q=2000，Δx=20m; (f)FD，Q=20，Δx=20m

3.3 穩(wěn)定性分析

穩(wěn)定性條件是有限差分?jǐn)?shù)值模擬中一個非常重要的問題，是影響差分方法計算效率的重要因素。采用與文獻(xiàn)[31]相同的Fourier方法對黏滯聲波方程的組合型緊致差分格式進(jìn)行穩(wěn)定性分析，可得穩(wěn)定性條件為

(20)

式中Ψ=19.2。可以看出，時間步長Δt的取值不僅與網(wǎng)格間距h和波速v有關(guān)，而且與介質(zhì)品質(zhì)因子Q和諧波頻率ω有關(guān)。當(dāng)Q趨向無窮大時，黏彈介質(zhì)可以看作完全彈性介質(zhì)，則時間二階、空間六階精度的CCD格式的穩(wěn)定性條件變?yōu)閏0=vmaxΔt/h<0.456，其中vmax為最大波速。

對式(20)進(jìn)行數(shù)值分析，依然定義c0=vΔt/h，并計算它隨速度、品質(zhì)因子、頻率的變化曲線，結(jié)果如圖4所示，可以看出：①當(dāng)品質(zhì)因子和頻率一定時，滿足穩(wěn)定性要求的庫朗數(shù)隨速度的增加而減小，即穩(wěn)定性變差，要求的時間步長越小。②當(dāng)?shù)貙铀俣纫欢?，品質(zhì)因子小于100時，庫朗數(shù)隨品質(zhì)因子的減小而急劇變小。當(dāng)品質(zhì)因子增加時，庫朗數(shù)隨之增加，穩(wěn)定性變好，滿足要求的時間步長越大。③當(dāng)品質(zhì)因子達(dá)到2000時，介質(zhì)可以近似為完全彈性介質(zhì)，庫朗數(shù)趨近于常數(shù)0.456，不隨速度的變化而改變。④穩(wěn)定性還和地震子波主頻有關(guān)，當(dāng)速度和品質(zhì)因子不變時，穩(wěn)定性隨頻率的降低而變差。由于黏彈介質(zhì)中地震波主頻會降低，降低的快慢與介質(zhì)的品質(zhì)因子有關(guān)，所以要求的穩(wěn)定性會隨傳播時間的增加而改變，即要求的時間步長應(yīng)該越來越小。若僅按照激發(fā)地震主頻得到的時間步長進(jìn)行模擬，在傳播一定距離后，差分格式可能就不穩(wěn)定了，這與完全彈性介質(zhì)波動方程數(shù)值模擬不同，因為完全彈性介質(zhì)中的地震波只會產(chǎn)生振幅衰減，而主頻不會降低。

圖4 差分格式式(5)的庫朗數(shù)與品質(zhì)因子、速度和頻率的關(guān)系曲線(a)速度—庫朗數(shù)曲線; (b)品質(zhì)因子—庫朗數(shù)曲線; (c)頻率—庫朗數(shù)曲線

3.4 差分格式優(yōu)化

為了進(jìn)一步提高組合型緊致差分格式的緊致性，在差分格式式(1)的基礎(chǔ)上建立四點組合型緊致差分格式

(21)

為了滿足2～6階精度泰勒公式截斷誤差的要求，上式中的差分系數(shù)ai,bi和ci必須滿足

(22)

為確定上述方程中的未知差分系數(shù)ai、bi和ci，按照Tam等[32]提出的頻散關(guān)系保持的思路，求取在最小數(shù)值頻散條件下的差分系數(shù)[33，34]。

與3.2相同，可以得到式(21)的頻散關(guān)系

(23)

式中：φ′=k′h；φ″=k″h。k、k′和k″分別表示真波數(shù)、一階和二階導(dǎo)數(shù)的數(shù)值波數(shù)(或稱修正波數(shù))。

優(yōu)化的目標(biāo)是在某個選定的波數(shù)范圍內(nèi)，確定式(21)中的8個未知差分系數(shù)，使得修正波數(shù)盡可能地接近真波數(shù)，所以定義誤差函數(shù)為

(24)

式中：R(φ′)表示φ′的實部；W(φ)為一個加權(quán)函數(shù)，目標(biāo)是使誤差函數(shù)解析可積，它是φ-R(φ′)的分母部分。選定c2和c3為變量，則確定Er的最小值為條件極值問題。采用拉格朗日乘數(shù)法進(jìn)行求解，可得

(25)

由式(22)和式(25)表示的8個方程，可以確定8個優(yōu)化后的差分系數(shù)為：a1=0.8873686;a3=0.0491178；b1=0.149532；b2=-0.2507682；b3=-0.0123598；c1=0.0163964；c2=-1.9692791；c3=1.9528828。將差分系數(shù)代入式(19)，得到優(yōu)化后的組合型緊致差分格式，簡稱OCCD。優(yōu)化前、后頻散關(guān)系曲線，如圖5所示。

圖5 CCD和OCCD格式頻散曲線(a)波數(shù)—修正波數(shù)曲線; (b)速度比曲線

從圖5可以看出，差分系數(shù)優(yōu)化后，不論是一階導(dǎo)數(shù)還是二階導(dǎo)數(shù)的修正波數(shù)都更接近真波數(shù)，速度比也更接近理想值1。優(yōu)化前，滿足一階和二階導(dǎo)數(shù)速度比等于1(或修正波數(shù)等于真波數(shù))的最小φ為0.1833，而優(yōu)化后滿足同樣條件的最小φ為0.2556，說明優(yōu)化后的OCCD格式提高了組合型緊致差分格式的緊致性，能更好地壓制數(shù)值頻散，從而能使用更粗的網(wǎng)格，進(jìn)一步提高了計算效率。

4 PML邊界條件

本文采用完全匹配層(Perfectly matched layer，PML)吸收邊界。按照王守東[35]和劉有山等[36]推導(dǎo)聲波PML控制方程的思路和方法，略去推導(dǎo)過程，直接給出二階黏滯聲波方程的PML邊界條件的控制方程

(26)

式中：u1、u2、A1和A2是引入的中間變量；d(x)和d(z)分別是x和z方向的衰減系數(shù)，分別衰減這兩個方向傳播的波。衰減函數(shù)采用高剛等[37]提出的余弦型吸收衰減函數(shù)，衰減幅度因子為500，吸收邊界厚度為20個網(wǎng)格間距。

5 模型試算

利用均勻介質(zhì)模型和Marmousi模型說明OCCD差分格式的黏滯聲波方程數(shù)值模擬效果。

5.1 均勻介質(zhì)模型

圖6為Q=20、50、200和完全彈性(Q=+∞)四種情況時模擬的220ms時刻波場快照，波場非常清晰，沒有數(shù)值頻散現(xiàn)象。

提取以上四個波場快照中的波形記錄(z=1250m)，如圖7所示。從中可以看出，波形曲線平滑，無數(shù)值頻散，且隨著品質(zhì)因子的減小，地震波振幅發(fā)生顯著衰減，頻率降低，波形展寬，波長增加。

為了說明本文差分方法在PML條件下對邊界反射的吸收效果，將圖6b的傳播時間延長至410ms，圖8a和圖8b分別為未使用和使用PML邊界條件時的波場快照，圖8c為二者之差，是只含邊界反射的波場快照?？梢钥闯?，本文差分格式在PML控制方程作用下，邊界反射得到了較好的吸收，有效波沒有明顯改變，邊界處理效果可靠。

圖6 均勻模型OCCD差分格式在Q=20(a)、50(b)、200(c)和+∞(d)時模擬的220ms時刻波場快照

圖7 均勻模型不同品質(zhì)因子時的波形曲線

5.2 Marmousi模型

用經(jīng)典的二維Marmousi縱波速度模型進(jìn)行數(shù)值模擬。速度范圍為1729～5500m/s，品質(zhì)因子根據(jù)經(jīng)驗公式Q=14v2.2計算，范圍為46.7～595.6(圖9)。模型網(wǎng)格數(shù)為501×501，網(wǎng)格間距為5m，時間步長為0.2ms，縱波震源位于(1250m，0)，激發(fā)40Hz的Ricker子波，采樣時間2s。

圖10為對Marmousi模型分別進(jìn)行黏滯聲波方程和聲波方程數(shù)值模擬得到的波場快照，可見波場快照平滑清晰，無數(shù)值頻散，且邊界吸收效果較好，無明顯邊界反射，這說明本文算法對復(fù)雜模型的適用性。圖11是兩種方程分別模擬得到的地面地震記錄。

圖8 未使用(a)和使用(b)PML邊界條件時的410ms時刻波場快照及二者之差(c)

圖9 Marmousi模型(a)速度; (b)品質(zhì)因子

圖10 Marmousi模型黏滯聲波方程(上)和聲波方程(下)模擬的波場快照(a)260ms; (b)460ms; (c)660ms; (d)860ms

圖11 Marmousi模型黏滯聲波方程(a)和聲波方程(b)模擬的地震記錄

由于波前擴(kuò)散造成振幅衰減，深部反射振幅明顯減弱，為了顯示清晰，對地震記錄進(jìn)行了瞬時自動增益控制(AGC)處理，時窗長度為500ms。對比二者可以看出，采用黏滯聲波方程數(shù)值模擬時，除了波前擴(kuò)散造成的振幅衰減以外，地震波還受到了黏彈介質(zhì)的吸收衰減作用，使得反射波振幅衰減幅度大于聲波方程模擬的結(jié)果，且地震波的頻率也會隨傳播的距離的增加而減小。

圖12是從模擬得到的地面地震記錄中提取并增益處理后的單道記錄波形圖(x=600m)。從圖中可以看出，在同樣的AGC處理時，黏滯聲波記錄的振幅小于聲波記錄。此外，黏滯聲波記錄的子波延續(xù)時間大于聲波記錄，同時反射波的個數(shù)也少于聲波記錄，并且隨傳播時間的增加，差異更加明顯，這也導(dǎo)致了黏滯聲波方程模擬結(jié)果的垂向分辨率小于聲波。

為了分析兩種模擬記錄頻譜上的差別，利用廣義S變換，對圖12單道地震記錄進(jìn)行時頻分析，結(jié)果如圖13所示。由時頻分析結(jié)果可見，聲波模擬記錄(圖13b)的振幅會隨時間的增加而減小，這是波前擴(kuò)散造成的。在振幅衰減的同時，反射波的主頻基本不變，與激發(fā)主頻一致，約為40Hz。而黏滯聲波記錄(圖13a)除了振幅衰減以外，其主頻會降低，在300ms時，主頻為30Hz。在1100ms時，主頻變?yōu)?0Hz左右，在1600ms～1800ms時，主頻進(jìn)一步降至15Hz左右。對比兩種記錄的時頻譜，在1400ms～1900ms時，聲波模擬結(jié)果中有4個明顯區(qū)分的振幅譜能量峰值區(qū)，對應(yīng)時間域的4個主要反射波，而黏滯記錄只有一個很難區(qū)分的峰值區(qū)，這也從頻率域說明了地震垂向分辨率的降低。

圖12 x=600m處的單道地震記錄對比

圖13 Marmousi模型黏滯聲波方程(上)和聲波方程(下)模擬的x=600m處單道地震記錄時頻譜

6 結(jié)論與建議

本文從組合型緊致差分格式出發(fā)，建立了時間二階、空間六階的黏滯聲波方程的差分格式，對該差分格式進(jìn)行了模擬精度、頻散關(guān)系、穩(wěn)定性和格式優(yōu)化分析。最后在黏滯聲波PML控制方程的基礎(chǔ)上，利用優(yōu)化組合型緊致差分格式，對均勻介質(zhì)和Marmousi模型進(jìn)行了數(shù)值模擬。數(shù)值模擬結(jié)果顯示地震波場特征清晰，邊界反射吸收效果較好，驗證了該方法的適用性，也說明該方法還可以進(jìn)一步推廣到二維或三維的各向異性介質(zhì)和雙相介質(zhì)等復(fù)雜介質(zhì)的聲波或彈性波數(shù)值模擬中。

文中建立的黏彈介質(zhì)聲波方程為時間二階精度，在今后的研究中，應(yīng)考慮提高時間差分精度，以提高數(shù)值模擬的計算精度。此外，可以通過其他有限差分系數(shù)優(yōu)化的方法[38，39]優(yōu)化該差分格式，提高計算精度和效率。