二維黏滯聲波方程的優化組合型緊致有限差分數值模擬

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

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

1 引言

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

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

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

組合型緊致差分格式在黏滯聲波方程地震波場模擬中的應用，尚未見到文獻報道。本文在黏滯聲波方程時間二階離散格式的基礎上，將組合型緊致格式應用于位移場空間導數的求取，實現了二維黏滯聲波方程地震波場的數值模擬。

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

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

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

(1)

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

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

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

(2)

利用時間二階導數的二階精度近似式

(3)

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

(4)

(5)

(6)

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

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

3.1 精度分析

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

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

3.2 頻散分析

頻散分析既是判斷數值模擬方法優劣的重要指標，也是確定空間網格大小的重要依據。通過對黏滯聲波方程的組合型緊致差分格式進行數值頻散分析，進一步分析方法的適用條件。

首先考慮x方向，令

(15)

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

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

(16)

同理可得z方向的解

(17)

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

(18)

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

(19)

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

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

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

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

圖3 均勻模型CCD和常規FD格式不同品質因子和網格間距時的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 穩定性分析

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

(20)

式中Ψ=19.2?？梢钥闯觯瑫r間步長Δt的取值不僅與網格間距h和波速v有關，而且與介質品質因子Q和諧波頻率ω有關。當Q趨向無窮大時，黏彈介質可以看作完全彈性介質，則時間二階、空間六階精度的CCD格式的穩定性條件變為c0=vmaxΔt/h<0.456，其中vmax為最大波速。

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

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

3.4 差分格式優化

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

(21)

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

(22)

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

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

(23)

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

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

(24)

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

(25)

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

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

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

4 PML邊界條件

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

(26)

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

5 模型試算

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

5.1 均勻介質模型

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

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

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

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

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

5.2 Marmousi模型

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

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

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

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

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

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

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

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

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

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

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

6 結論與建議

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

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