汪 勇,穆鵬飛,蔡文杰,王 鵬,桂志先

(1.油氣資源與勘探技術教育部重點實驗室(長江大學),湖北武漢430100;2.長江大學地球物理與石油資源學院,湖北武漢430100;3.中海石油(中國)有限公司天津分公司,天津300450)

地震波場有限差分數值模擬[1-9]是研究復雜地區(qū)地震資料采集、處理和解釋的有效輔助手段,也是逆時偏移和全波形反演的基礎。緊致有限差分格式是其中一種具有較高計算效率的方法,與常規(guī)中心差分格式相比,它能夠使用更少的計算節(jié)點達到相同的計算精度,減少了計算量和所需的存儲空間。此外,緊致差分格式是一種隱式差分格式,使用該格式計算空間導數時,不僅用到了待求網格點周圍的函數值,還用到了相鄰節(jié)點的導數值,具有更低的數值頻散誤差,因而能夠使用更粗的網格計算,提高數值模擬的計算效率。因此緊致差分格式被廣泛應用于聲波、彈性波和復雜介質等的地震波場數值模擬[10-17]。

在地震波場有限差分數值模擬中,需要高精度的空間離散差分算子。差分算子的精度取決于差分系數和差分階數,差分階數越高,差分算子越精確,差分精度就越高,數值頻散越小,但其計算量也隨之增大。優(yōu)化差分系數可使差分算子最大程度地逼近空間偏導算子。KIM等[18]利用頻散關系保持(dispersion-relation-preserving,DRP)的基本思想優(yōu)化了一階導數的緊致差分格式;LIU等[19]基于頻散關系提出了優(yōu)化的時空域有限差分系數,可在不增加計算成本的情況下顯著提高模擬精度;ZHANG等[20-21]使用最大范數的目標函數及模擬退火算法求解目標函數,對一階和二階常規(guī)中心差分的差分系數進行了優(yōu)化;LIU[22-23]使用最小平方法優(yōu)化了二階導數中心差分和一階導數交錯差分系數;YU等[24]基于DRP思想,優(yōu)化得到了5階精度的組合型緊致差分系數;REN等[25]利用優(yōu)化后的時空域差分格式進行了聲波和彈性波方程的數值模擬;YANG等[26]利用極小極大算法(minimax approximation)優(yōu)化了交錯網格差分系數,在保證中低波數范圍內差分格式分辨率的同時,提高了大波數情況下的模擬精度。

本文首先將二階導數的五對角緊致差分格式擴展到2N階差分精度,然后基于頻散關系保持的思想,根據最小平方法建立了波數誤差的目標函數,利用拉格朗日乘數法對目標函數進行求解,得到了優(yōu)化后的4～10階精度差分系數。與KIM等[18]方法相比,除了優(yōu)化目標是二階導數以外,主要區(qū)別在于使用的方程和網格節(jié)點數不同,KIM等[18]方法在優(yōu)化2～8階精度差分系數時均使用了7個節(jié)點,而本文根據差分精度的不同,使用不同的節(jié)點數,提高了計算效率。最后,分析和對比了優(yōu)化前后差分格式的模擬精度、數值頻散和聲波方程的穩(wěn)定性條件,并應用優(yōu)化后的緊致差分格式和基于輔助微分方程的完全匹配層邊界條件對二階聲波方程進行了數值模擬,驗證了方法的精度。

1 二維聲波方程的高階差分近似

根據彈性力學分析,二維情況下二階聲波方程(假設體力為零)可以表示為:

(1)

式中:u(x,z)為聲波位移場;v(x,z)為聲波速度場;x,z和t分別為空間和時間坐標。

1.1 時間2M階近似

利用截斷的泰勒公式表示n+1和n-1時刻的位移場可以得到:

式中:i,j和n分別為空間和時間網格坐標;Δt為時間步長。兩式相加,略去高次項,得到位移場時間2M階精度的差分格式:

(4)

當M=1時,根據聲波方程((1)式),將(4)式中位移對時間的導數轉化為位移對空間的導數,即可得時間二階精度的差分格式:

(5)

當M=2時,同理可得時間四階精度的差分格式:

(6)

公式(5)和公式(6)為位移場的三層顯式差分格式,利用它們就可以進行地震波場時間層的推進計算。公式中含有位移u對x和z的二階和四階導數,這些導數將利用緊致差分格式進行求取。

1.2 空間2N階緊致差分近似

1992年,LELE[27]對埃爾米特(Hermite)公式進行了擴展,構造了求解函數f(x)二階導數的緊致差分格式為:

(7)

對公式(7)進行擴展,可以構造二階導數的2N階精度緊致差分格式。由于公式(7)左邊有5項,系數矩陣為五對角矩陣,所以稱為五對角緊致有限差分格式(pentadiagonal compact finite difference,簡稱CFD5),表示為:

(8)

從公式(8)可以看出,CFD5格式只需要利用2N-3個節(jié)點即可達到2N階差分精度。利用泰勒級數展開和待定系數法以及下列方程組((9)式)可以求得CFD5格式2N(N≥3)階差分精度的差分系數,表1列出了6～12階精度CFD5格式差分系數。

(9)

利用公式(8)所表示的CFD5格式,可以求取公式(5)和公式(6)中的位移對空間的二階偏導數?2u/?x2和?2u/?z2,差分格式寫成矩陣形式為:

(10)

其中,

根據公式(10)即可求得空間二階導數?2U/?x2和?2U/?z2的2N(N≥3)階空間差分精度近似值,表示為:

(11)

此外,公式(6)中的四階偏導數和混合偏導數,可以利用公式(11)對二階偏導數再次求導進行求取。

表1 五對角緊致差分格式二階導數的差分系數

2 二階導數緊致格式的差分系數優(yōu)化及分析

地震波場數值模擬中,為了得到清晰的波場記錄,需要高精度的空間離散差分算子。精度越高的差分算子數值頻散越小,但算子長度越大,計算量也越大。優(yōu)化差分系數可使差分算子最大程度地逼近空間偏導算子,本節(jié)基于頻散關系保持的基本思想,對二階導數的緊致差分系數進行優(yōu)化,提高緊致差分格式的分辨率。

2.1 優(yōu)化差分系數

首先對CFD5格式(公式(8))進行數值頻散分析,令:

(12)

將公式(12)代入公式(8)中可以得到(2N階精度):

(13)

利用歐拉公式,并將B=-(k′)2A代入(13)式可得:

(14)

令ω=kh,ω′=k′h,則ω′可以表示為:

(15)

ω和ω′分別是真波數和數值波數與空間步長的乘積,在理想情況下,如果不存在數值頻散,則ω′=ω。它們的差別越大,則說明該方法的數值頻散越嚴重,反之則說明該方法能更好地壓制數值頻散。

下面以3點6階格式為例來說明差分系數的優(yōu)化方法。由公式(8)可以得到:

(16)

公式(16)中的差分系數a1,a2和b1為了滿足2階和4階精度泰勒公式截斷誤差的要求,必須滿足方程組:

(17)

按照TAM等[28]提出的DRP思路,在某個選定的波數范圍內,確定公式(17)中的3個未知差分系數a1,a2和b1,使得數值波數k′盡可能地接近真波數k。為了消除數值波數表達式中的根號,簡化計算,定義誤差函數為:

(18)

其中,W(ω)為一個加權函數,目標是使誤差函數解析可積,它是ω2-ω′2的分母部分。θ是積分上限,取值范圍是θ∈(0,π),取值大小與差分精度有關,一般低階取小值,高階取大值。它的取值不同,計算出來的差分系數也不同,這里θ=5π/8。

選定a1為變量,由方程(17)可得b1=(36a1+24)/23,a2=(1-10a1)/46。確定E的最小值為條件極值問題,采用拉格朗日乘數法進行求解,即根據?E/?a1=0,可以得到:

(19)

由公式(17)、公式(18)和(19)可以確定3個優(yōu)化后的差分系數,a1=0.127658,a2=-0.006013,b1=1.243290。將上述差分系數代入公式(8),則得到優(yōu)化后的4階精度五對角緊致差分格式(pentadiagonal optimized compact finite difference,以下簡稱OCFD5)。采用同樣的方法,對公式(8)表示的CFD5格式進行差分系數優(yōu)化,優(yōu)化后的4～10階差分系數見表2。為了方便,將二階導數的2N階精度優(yōu)化前后的緊致差分格式統一表示為:

(20)

CFD5:L=N-2,N≥3,a1,a2和bl見表1。OCFD5:L=N-1,N≥2,a1,a2和bl見表2。

若要達到2N階空間差分精度,優(yōu)化前的CFD5方法需要利用2N-3個節(jié)點,優(yōu)化后的OCFD5方法需要使用2N-1個節(jié)點。

表2 二階導數的五對角優(yōu)化緊致差分格式的差分系數

由于不同的積分上限可以得到不同的差分系數,這里補充說明積分上限選取的原則。通過選取θ=3π/8～7π/8,可以得到不同的差分系數,表3為4階精度差分格式優(yōu)化時,不同積分上限時得到差分系數。由不同的差分系數可以得到各自的波數比曲線,如圖1所示,其中波數比定義為數值波數與真波數之比,表示為R=k′/k=k′h/kh=ω′/ω。

由4階精度格式優(yōu)化結果可以看出(圖1a),當θ=5π/8,波數比曲線接近1的范圍最大。θ取值較小時,有效的波數范圍較小,達不到優(yōu)化的目的;θ取值過大時,曲線向上偏離1,也會造成數值頻散。定量來看,若給定頻散誤差標準E=|R-1|≤0.2%,當積分上限θ=3π/8,4π/8,5π/8,6π/8,7π/8時,波數比R滿足誤差標準的最大ω分別為1.55,1.69,1.92,1.12和0.96。也就是說,當θ=5π/8時的滿足誤差標準的波數范圍最寬,能在最大的范圍內壓制數值頻散。圖1b到圖1d是6,8,10階OCFD5格式在不同積分上限時得到的波數比曲線,與4階精度時確定積分上限和對應優(yōu)化系數的原則類似,確定6階和8階時的積分上限θ=6π/8,10階時的積分上限θ=7π/8。由于篇幅所限,6～10階不同積分上限時的差分系數不具體列出。

2.2 頻散分析

CFD5格式和OCFD5格式的數值波數與真波數之比R可以表示為:

圖1 不同積分上限θ優(yōu)化后格式的波數比曲線a 4階精度; b 6階精度; c 8階精度; d 10階精度

積分上限3π/84π/85π/86π/87π/8a10.1276580.1312350.1366110.1446050.156566a2-0.006013-0.006790-0.007959-0.009697-0.012297b11.2432901.2488901.2573041.2698171.288538

(21)

在理想情況下,即不存在數值頻散時,R恒等于1。R偏離1越大,說明該方法的數值頻散越嚴重,反之則說明該方法能更好地壓制數值頻散。使用表1和表2中不同精度的差分系數,可以計算得到CFD5和OCFD5格式的波數比曲線,如圖2所示。

由圖2波數比曲線可以看出:①在相同的差分精度情況下,優(yōu)化后的格式比優(yōu)化前的波數比曲線更接近1,所以每個波長內需要的最少采樣點數更少,這說明優(yōu)化過程起到了提高壓制數值頻散的作用,數值計算時能使用更大的空間網格,從而減少了計算內存,提高了計算效率;②優(yōu)化后的2N階精度格式,波數比曲線比2N+2階精度優(yōu)化前的頻散誤差更小,例如圖2a中的4階OCFD5要好于6階的CFD5。由圖2b到圖2d也可以得出相同結論。

2.3 差分精度分析

進行數值模擬時,在時間差分精度相同的條件下,不論是利用優(yōu)化前的CFD5格式,還是利用優(yōu)化后的OCFD5格式,它們在時間層遞推方式上是相同的,不同之處僅在于空間偏導數近似值求取時的差分格式不同。為了比較它們的近似精度,需要對其截斷誤差進行對比,結果見表4。從表4中可以看出,在兩種方法達到同樣空間近似精度的條件下,優(yōu)化后的格式具有更小的截斷誤差,例如6階精度時的CFD5截斷誤差約為OCFD5的3.69倍。

利用一維平面簡諧波初值問題的解析解來比較優(yōu)化前后格式的數值模擬精度。一維平面諧波初值問題可以表示為:

圖2 不同差分精度時CFD5和OCFD5格式的波數比曲線a 4階和6階精度比較; b 6階和8階精度比較; c 8階和10階精度比較; d 10階和12階精度比較

表4 優(yōu)化前后格式的二階導數截斷誤差主項系數比較

(22)

式中:v是平面波的波速;f0是平面簡諧波的頻率。

其初值問題的解析解為:

(23)

設置一維介質模型長度8km,平面波速v=2000m/s,f0=10Hz,數值模擬時的時間步長2ms,空間網格大小20m。數值模擬時,兩種差分格式均為時間4階、空間6階精度。通過數值模擬,可得圖3 所示的1s時刻的位移u的波場快照,圖中藍色實線是用公式(23)計算得到的精確解析解。從圖3可以看出,兩種差分格式模擬精度均很高,圖3a中的3條曲線基本重合,無明顯誤差。但從局部放大的圖3b 和圖3c來看,優(yōu)化后的OCFD5格式的數值模擬結果與精確的解析解更接近,說明優(yōu)化格式的模擬精度得到了提高。

為了避免邊界反射的影響,定量計算得到1s時刻3000～5000m區(qū)間的CFD5和OCFD5數值解與解析解之間的相對誤差分別為0.1005%和0.0167%。這是由于它們在計算空間導數近似值時存在不同截斷誤差大小造成的,這與表4的理論分析結果一致。相對誤差定義為:

(24)

圖3 一維聲波方程數值模擬快照a 3000～5000m范圍的波場快照; b 圖3a中A部分放大顯示; c 圖3a中B部分放大顯示

2.4 穩(wěn)定性條件

穩(wěn)定性條件是有限差分數值模擬中一個非常重要的問題,是影響差分方法計算效率的重要因素。我們采用Fourier方法對公式(5)表示的二維聲波方程的時間2階、空間2N精度差分格式進行穩(wěn)定性分析。

定義:

(25)

式中:kx和kz為x和z方向的視波數;i,j和n分別為空間和時間網格坐標;ξ和ζ1為振幅。代入公式(20)中可得:

(26)

令θ1=kxh,-π/2≤θ1≤π/2,并利用歐拉公式,將(26)式化簡可得:

(27)

解方程可得:

(28)

同理,對于z方向,定義:

(29)

代入公式(20)同樣可以得到:

(30)

公式(5)是一個三層顯式差分格式,為了分析其增長矩陣,將其改寫為:

(31)

(32)

若使差分格式(32)滿足穩(wěn)定性條件,則增長矩陣G最大特征值的模不大于1,由此可得穩(wěn)定性條件為:

(33)

由公式(28)和公式(30)可以得到:

(34)

所以二維聲波方程的時間2階、空間2N階精度差分格式的穩(wěn)定性條件可寫為:

(35)

式中:Δt為時間步長;h為空間網格大小;v為地震波速度;a1,a2和bl為表1和表2中的緊致差分格式差分系數。

定義α=vΔt/h,穩(wěn)定性條件寫作α≤C(C稱為Courant數),聲波方程時間2階精度的CFD5和OCFD5格式的穩(wěn)定性條件如表5所示。從表5來看,C隨時間差分精度的增加而增加,隨空間差分精度的增加而減小。在同樣的差分精度條件下,優(yōu)化后格式的穩(wěn)定性條件比優(yōu)化前要略微嚴格,也就是說在相同空間步長時,允許的時間步長要略小。聲波方程時間4階精度的穩(wěn)定性條件很難給出具體表達式,表5 給出的不同情況下的穩(wěn)定性條件是根據試驗獲得的。

2.5 邊界條件

在數值模擬時,由于模型的設置和計算模型大小的限制,必然會存在人工邊界。如果不對人工邊界進行處理,邊界反射會干擾正常的地震波場,所以人工邊界的處理直接影響到數值模擬的精度和效率。1994年,BERENGER[29]提出的完全匹配層(perfectly matched layer,簡稱PML)能對邊界反射起到很好的吸收作用。在此基礎上,DROSSAERT等[30-31]提出的復頻移完全匹配層(complex frequency shifted perfectly matched layer,簡稱CFS-PML)能對近掠入射的低頻波起到更好的衰減作用,而卷積PML[32](convolutional PML)和輔助微分方程PML[33](auxiliary differential equation PML,簡稱ADE-PML)的優(yōu)點是可以使用非分裂的波動方程。針對二階聲波方程,本文采用輔助微分方程結合復頻移完全匹配層[34](ADE-CFS-PML)進行邊界處理,略去推導過程,直接給出x方向控制方程為:

表5 二維聲波方程優(yōu)化前后緊致差分格式的穩(wěn)定性條件

(36)

式中:u1,u2和u3是引入的中間變量,可由(36)式中第一個方程到第三個方程依次計算u3,u2和u1;參數α′和d′分別是α(x)和d(x)對x的一階導數。α(x)和d(x)表示為:

(37)

式中:d(x)是x方向的衰減系數,起到衰減x方向波的作用;αmax=πf0,f0是震源信號主頻;δ是PML邊界厚度;x是節(jié)點至PML最內層邊界的距離。此外,d0=[3vmaxln(1/F)]/(2δ)[35],vmax是地震波傳播的最大速度,F是理論反射系數,表示為lgF=[1-lgN]/lg2-3,N是PML邊界層數,本文N=20。

公式(36)只表示了x方向的控制方程,同樣方法可以得到z方向的控制方程,文中不再贅述。

3 模型試算

3.1 均勻介質模型

由圖4可以看出:①當使用相對較粗網格計算時,圖4a所示的6階CFD5格式的波場快照頻散嚴重,與2.2節(jié)的理論分析結果一致;②增加網格點數和差分階數后,圖4b所示的8階CFD5格式的波場快照的頻散得到了壓制,但頻散還是較為明顯;③圖4c 所示的6階OCFD5格式的波場快照數值頻散明顯減弱,模擬結果不僅優(yōu)于優(yōu)化前同階精度的圖4a,也好于差分精度更高的圖4b;④圖4d所示的8階OCFD5格式的波場快照清晰,波形光滑,無明顯數值頻散。優(yōu)化前后波場快照的比較結果,說明了本文基于DRP思路進行優(yōu)化的結果達到了壓制數值頻散的目的。

圖4 不同差分格式計算的1000ms時刻波場快照a CFD5 (4,6); b CFD5 (4,8); c OCFD5(4,6); d OCFD5(4,8)

基于該模型進行計算效率分析,結果如表6所示。表6中所取空間網格大小滿足波場快照無數值頻散要求,同時為了簡化分析和避免時間步長過大產生的數值頻散,時間步長均為3ms。4種格式的數值模擬沒有本質的差別,計算中需要的變量(數組)個數一樣,只是差分系數和算子長度不同,它們均在相同的運行環(huán)境和程序下進行數值模擬。從表6中可以看出,優(yōu)化后的格式由于能使用更粗的網格計算,所以單個變量的網格數降低,既減少了內存的占用,也減少了計算時間,從而提高了數值模擬的計算效率。

表6 不同方法計算效率比較

3.2 Marmousi模型

為了驗證優(yōu)化后的緊致差分格式對復雜介質的適用性,用經典的二維Marmousi縱波速度模型進行數值模擬。模型如圖5所示,速度v的范圍是1729～5500m/s。模型大小為501×501個網格點,空間網格12m,時間步長1ms,震源位于(x=3000m,z=0)位置,激發(fā)震源采用30Hz的Ricker子波,采樣時間4s,利用ADE-CFS-PML邊界條件對人工邊界反射進行吸收處理。

圖5 Marmousi縱波速度模型

圖6給出了采用OCFD5(2,6)格式對Marmousi模型進行聲波方程數值模擬得到的不同時刻的波場快照。從圖6可以看出,Marmousi模型中地震波場復雜且清晰,能夠反映地震波的傳播特征。同時,邊界反射吸收效果較好,無明顯邊界反射,說明本文差分方法在ADE-CFS-PML邊界條件下,對邊界反射起到了很好的衰減作用。

圖7是分別使用優(yōu)化前后格式得到的地震記錄,精度均為時間2階、空間6階。為了顯示清晰,對地震記錄進行了瞬時自動增益控制(AGC)處理,時窗2s。

從圖7中方框所示范圍內的地震波場來看,優(yōu)化前CFD5格式的數值模擬結果出現了較嚴重的數值頻散,干擾了正常的地震波場,不利于波場特征分析及其它處理;而優(yōu)化后OCFD5格式的數值模擬結果,模擬精度高,無明顯數值頻散。同時,直達波、反射波及繞射波等波型清晰可見,且邊界吸收效果較好,無明顯邊界反射,驗證了本文算法對復雜模型的適用性。

圖6 Marmousi模型OCFD5(2,6)格式聲波模擬波場快照a 400ms; b 800ms; c 1200ms; d 1600ms

圖7 Marmousi模型CFD5(a)格式和OCFD5(b)格式模擬的地面地震記錄

4 結論與建議

本文基于頻散關系保持的思想,利用最小平方法和拉格朗日乘數法對二階導數的五對角緊致差分格式的差分系數進行了優(yōu)化,研究認為:

1) 同為2N階差分精度時,優(yōu)化后的差分格式具有更小的數值頻散和截斷誤差,能夠使用更粗的網格進行地震波場數值模擬,節(jié)省了計算內存,更加適用于粗網格下的大尺度的地震波場數值模擬。

2) 相同網格參數時,優(yōu)化后2N階OCFD5格式在壓制數值頻散方面不僅優(yōu)于2N階CFD5格式,也好于2N+2階CFD5格式,也就是可以使用更少的計算節(jié)點(即算子長度)的同時提高計算效率。

3) 優(yōu)化差分系數以后,聲波方程的穩(wěn)定性條件要比優(yōu)化前略微嚴格,但總體上差別不大,Courant數隨空間差分精度的增加而減小,隨時間差分精度的增加而增加。

4) 優(yōu)化后差分格式的數值模擬結果,地震波場特征清晰,驗證了該方法的適用性,為研究地震波傳播規(guī)律、逆時偏移、全波形反演等工作提供了一種有效的波場延拓方法。

在今后的研究中可以考慮使用不同的優(yōu)化方法,例如樣點逼近、模擬退火法和極小化極大算法等方法,進一步提高優(yōu)化差分系數的效果。