朱　強,姜蘆倩,張　偉

(1.中國科學技術(shù)大學地球和空間科學學院,安徽合肥230026;2.南方科技大學地球與空間科學系,廣東深圳518055)

地震波場數(shù)值模擬是研究地震波傳播規(guī)律的主要手段,也是基于波動方程進行地震反演和成像的基礎。有限差分方法具有實現(xiàn)簡單、應用范圍廣、內(nèi)存需求小、易于并行的優(yōu)點,被廣泛應用于理論地震學和勘探地震研究中。目前的有限差分算法可以根據(jù)變量在計算網(wǎng)格上的分布不同而分成同位網(wǎng)格和交錯網(wǎng)格兩種網(wǎng)格系統(tǒng)。交錯網(wǎng)格布局又可以根據(jù)波場分量在計算網(wǎng)格上的相對位置關(guān)系被大致分為三種格式:標準交錯網(wǎng)格[1-2]、完全交錯網(wǎng)格[3-4]和旋轉(zhuǎn)交錯網(wǎng)格[5]。交錯網(wǎng)格有限差分方法通常采用中心差分的方式求解一階速度-應力方程組,可以避免同位網(wǎng)格簡單中心差分格式的奇偶失聯(lián)問題[6],是目前求解波動方程數(shù)值解的主流方法之一。

有限差分波動方程數(shù)值解存在不同類型的誤差,其中與地震波頻率有關(guān)的色散/耗散誤差可以通過使用高階差分算子或者優(yōu)化格式來壓制[7-10]。而包含內(nèi)部速度間斷面的復雜速度模型采用有限差分網(wǎng)格離散時,由于網(wǎng)格與速度界面不重疊,同一物理模型采用不同的介質(zhì)離散方法,因此會導致不同的網(wǎng)格離散模型,進而引起正演結(jié)果的不同。ALTERMAN等[11]采用虛擬網(wǎng)格點和顯式實施內(nèi)部邊界條件的方式處理二階位移方程中的水平間斷面問題。VIRIEUX[1,12]發(fā)展了基于一階速度-應力方程的交錯網(wǎng)格有限差分方法,通過網(wǎng)格點離散取值表征速度模型,并且展示了這種方法在大泊松比情況下的穩(wěn)定性。GRAVES[13]給出了交錯網(wǎng)格有限差分算法處理三維彈性介質(zhì)參數(shù)的格點平均方法。MUIR等[14]首次將等效介質(zhì)參數(shù)理論[15]應用到有限差分方法中。MOCZO等[16]提出了體積分平均法,用來處理交錯網(wǎng)格有限差分算法中的介質(zhì)參數(shù)間斷面,解決了物理界面在網(wǎng)格格點之間變化時的表征問題。MOCZO等[17]、KRISTEK等[18]將離散后界面附近的介質(zhì)近似為正交各向異性介質(zhì)來實施間斷面的邊界條件,從而有效地壓制了界面誤差。LISITSA等[19]通過平面波方程分析了交錯網(wǎng)格有限差分方法處理界面的理論誤差。VISHNEVSKY等[20]對標準交錯網(wǎng)格、完全交錯網(wǎng)格和旋轉(zhuǎn)交錯網(wǎng)格在層狀界面下的誤差收斂情況進行了總結(jié)。但是,上述這些研究都沒有定量評估和給出在任意復雜模型中準確計算(滿足給定誤差要求)反射/透射波的網(wǎng)格大小要求,尤其是沒有考慮準確計算P-S轉(zhuǎn)換波所要求的網(wǎng)格大小,所以目前進行有限差分地震模擬通常僅考慮格式本身的色散和耗散誤差對網(wǎng)格大小的限制,而沒有考慮介質(zhì)離散方法不同導致的與介質(zhì)界面有關(guān)的反射/透射波計算準確性的影響。

本文首先介紹了交錯網(wǎng)格有限差分方法原理,然后總結(jié)了目前有限差分地震模擬中所使用的介質(zhì)離散方法,最后定量分析了不同介質(zhì)離散方法對不同地震波場模擬準確性的影響。

1　交錯網(wǎng)格有限差分方法

在二維笛卡爾坐標系下,波動方程的速度-應力形式表示為:

式中:vx與vz代表速度分量;τxx,τzz與τxz代表應力分量;ρ表示介質(zhì)密度;t表示時間;Cij(i=1,2,3;j=1,2,3)為彈性剛度矩陣,各向同性介質(zhì)情況下,Cij由兩個獨立的分量構(gòu)成,即拉梅常數(shù)λ與μ。

為了分析問題方便,假定計算網(wǎng)格是矩形,顯式M(M=2,4,6,8,…)階交錯網(wǎng)格差分算子Dx定義為:

(3)

式中:h表示網(wǎng)格步長,αm表示差分系數(shù)。則交錯網(wǎng)格格式如下:

(4)

2　介質(zhì)參數(shù)離散方法

對于介質(zhì)參數(shù)不連續(xù)的內(nèi)部界面,其邊界條件應該滿足位移連續(xù)與法向牽引力連續(xù)的條件。目前較為常用的介質(zhì)參數(shù)離散方法有直接取值法、格點平均法[13]、體積分平均法[16]與正交各向異性等效介質(zhì)法[17-18]。本文以一個傾斜界面模型為例,說明不同介質(zhì)離散方法的差異(圖1)。如圖1a所示,不同介質(zhì)的分界面傾斜穿過網(wǎng)格,傾角為30°,上層介質(zhì)參數(shù)為vP1=2000m/s,vS1=1410m/s,ρ1=1160kg/m3;下層介質(zhì)參數(shù)為vP2=4000m/s,vS2=2310m/s,ρ2=2050kg/m3。圖1b 是采用廣義反射-透射系數(shù)方法[21-23]按圖1a中的觀測系統(tǒng)計算的接收波形(通過水平層計算后旋轉(zhuǎn)得到),震源為炸藥震源。其中第一個事件(第一個檢波器0～1.0s之間)是直達P波,第二個事件(第一個檢波器1.5～2.0s之間)是反射P波,第三個事件是P-S轉(zhuǎn)換波(第一個檢波器2.3～2.6s之間)。以下介紹不同的介質(zhì)離散方法、圖1a所示模型用不同介質(zhì)離散方法得到的不同離散模型,以及與圖1b對應的不同事件計算結(jié)果的差異(圖2,圖3)。

圖1　斜層模型a 模型示意; b 道集

圖2　修正介質(zhì)參數(shù)示意a 直接取值法; b 格點平均法; c 體積分平均法; d 正交各向異性等效介質(zhì)法

圖3　斜層模型局部道集展示(圖1b中2～3s波形)a 直接取值法; b 格點平均法; c 體積分平均法; d 正交各向異性等效介質(zhì)法

2.1　直接取值法

直接取值法即直接將物理模型根據(jù)計算域劃分的網(wǎng)格予以離散,每個格點取其與物理模型相對應點的介質(zhì)參數(shù)。界面傾斜穿過網(wǎng)格時,該方法會形成階梯狀的介質(zhì)參數(shù)間斷界面,見圖2a。根據(jù)惠更斯原理,此時階梯上的每一個點相當于一個點震源,從而導致界面上產(chǎn)生人工的虛假波。圖3a展示了采用直接取值法離散介質(zhì)所產(chǎn)生的散射波形,可以看到在每道地震記錄的P-S轉(zhuǎn)換波后都伴隨著明顯的波形扭曲(第一個檢波器2.6s之后)。因此直接取值法不能精確模擬界面的實際位置,誤差較為明顯。

2.2　格點平均法

格點平均法對界面附近的介質(zhì)做了簡單的平均化處理,避免了介質(zhì)參數(shù)的劇烈間斷(圖2b),但這種方法本身無法對單位網(wǎng)格之內(nèi)的界面變化進行細致的描述,因此當界面位置與網(wǎng)格位置沒有耦合的時候,其模擬精度較差(圖3b)。

2.3　體積分平均法

當界面在單位網(wǎng)格間變化的時候,簡單的格點平均法無法準確反映界面位置。MOCZO等[16]提出的體積分平均法可以更好地描述界面位置,提高反射/透射波的模擬精度。在體積分平均法中,某個點在計算域的介質(zhì)參數(shù)不僅僅由這個點本身所在位置的介質(zhì)參數(shù)決定,還受到以它為中心的單位網(wǎng)格區(qū)域內(nèi)的介質(zhì)參數(shù)的影響。

設界面位于x=0處,在一維情況下對波動方程在單位網(wǎng)格內(nèi)做空間積分:

(7)

對(7)式應用均值定理,得到:

(8)

則可以定義密度分量ρ的積分平均為:

(9)

ρA反映了這個被積分單元的平均響應。將其推廣到二維情況下,同理有:

如公式(10)至公式(12)所示,對于網(wǎng)格上的某一點,應用體積分平均法對其周圍單位空間的介質(zhì)參數(shù)進行積分然后取平均,即可得到修正后的彈性介質(zhì)參數(shù)。圖2c展示了體積分平均法處理后的介質(zhì)參數(shù),與直接取值法和格點平均法相比,這種方法能夠反映界面在單位網(wǎng)格間的具體變化。同時,體積分平均法對傾斜界面做了合適的處理,與直接取值法相比有效壓制了因傾斜界面介質(zhì)離散不準確而產(chǎn)生的散射波(圖3c)。

2.4　正交各向異性等效介質(zhì)法

層狀介質(zhì)研究表明,由五個彈性系數(shù)描述的均勻橫向各向同性介質(zhì)是層狀各向同性介質(zhì)疊加物的長波等效[15,24]。其等效介質(zhì)的彈性參數(shù)可以通過對層狀介質(zhì)拉梅常數(shù)的合理平均來近似獲得。正交各向異性等效介質(zhì)法基于這一等效理論,通過將介質(zhì)界面附近的彈性參數(shù)等效為一種正交各向異性介質(zhì)彈性參數(shù),模擬界面在單位網(wǎng)格內(nèi)的不同位置。

二維各向同性介質(zhì)情況下:

(13)

式中:εxx,εxz與εzz表示應變分量。將(13)式應力-應變關(guān)系應用于一個水平介質(zhì)界面,τzz,τxz與εxx在界面兩端連續(xù),τxx,εxz與εzz在界面兩端不連續(xù),用上標“+”表示在界面兩端具有連續(xù)性的分量,上標“-”表示在界面兩端間斷的分量,則(13)式可以改寫為:

(14)

其中,

(15)

(14)式中彈性剛度矩陣被橫線與豎線分為4個部分,將每一個部分看作一個整體變量,解出界面兩端間斷的分量τ-與ε-,同時對等式兩端取平均可得:

(16)

式中,“〈〉”表示對單位網(wǎng)格求積分平均。(16)式反映了位于界面上的微小單元的應力-應變關(guān)系,由于〈τ+〉=τi,j,〈ε+〉=εi,j,因此可改寫為:

(17)

(18)

設M=λ+2μ,則:

(19)

式中:〈〉H表示調(diào)和平均,通過這樣的平均化處理將界面附近的介質(zhì)視為等效的正交各向異性介質(zhì),將其彈性剛度矩陣代入應力-應變關(guān)系以模擬界面上的邊界條件。與體積分平均法一樣,正交各向異性等效介質(zhì)法能夠準確地描述界面在單位網(wǎng)格之間的實際位置(圖2d),從而能夠獲得更加準確的地震記錄(圖3d)。

3　不同介質(zhì)離散方法對不同事件模擬準確性的影響

本文通過改變計算域網(wǎng)格的每波長網(wǎng)格數(shù)(points per wavelength,PPW)來觀察界面誤差的變化情況,分析了當每波長網(wǎng)格數(shù)逐漸增大時,由界面計算產(chǎn)生的不同反射/透射波的相對誤差。為了定量評估交錯網(wǎng)格有限差分算法對界面計算的準確性,首先定義相對誤差為:

(20)

式中:‖u‖2表示u的二范數(shù),ufd表示有限差分解,uref表示參考解。本文通過截取相應事件窗口的波形計算相對誤差,使用PPW來評估網(wǎng)格的分辨率。對

于Ricker子波,其PPW定義為:

(21)

式中:fc為Ricker子波中心頻率,Δh為空間網(wǎng)格步長,vmin表示P波或S波的最小速度。

在具體探索各種介質(zhì)參數(shù)處理方法在不同波形下的相對誤差PPW約束之前,首先觀察不同反射/透射角度對相對誤差的影響。如圖4所示設置一個兩層模型,上層介質(zhì)參數(shù)為vP1=2000m/s,vS1=1000m/s,ρ1=1200kg/m3,下層介質(zhì)參數(shù)為vP2=4000m/s,vS2=2000m/s,ρ2=2100kg/m3,震源位于(1000m,1500m)處。在距離界面上方400m與下方300m深度,分別沿著圖4中的虛線布置兩排檢波器,相對誤差隨不同反射/透射角的變化如圖5所示。由于某些位置反射/透射波與P-S轉(zhuǎn)換波發(fā)生了耦合,因此在展示相對誤差隨反射/透射角度的變化時,使用整體波形計算其相對誤差。

圖4　不同PPW相對誤差測試模型

3.1　透射波和反射波相對誤差隨網(wǎng)格的變化

在圖5中相對誤差較大的位置附近,同時在時間上能夠區(qū)分反射/透射波與P-S轉(zhuǎn)換波的前提下,將反射波接收點置于(1930m,2060m)處,透射波接收點置于(2260m,3160m)處。通過保持震源、檢波器與界面相對位置不變,逐漸縮小網(wǎng)格空間步長以得到不同的PPW。為了使獲得的結(jié)論更加一般化,將實際物理界面位置置于兩個網(wǎng)格之間(不與任何網(wǎng)格格點發(fā)生耦合),使用主頻為10Hz的Ricker子波作為震源子波,分別對比不同反射/透射波在使用不同介質(zhì)離散方法后的相對誤差。為了避免色散誤差的影響,本文采用基于泰勒展開的10階交錯網(wǎng)格格式,該格式滿足格式的色散誤差要求的PPW為5左右,而下面測試的PPW最小是6,因此所有測試都滿足格式的色散誤差要求。

圖5　相對誤差隨偏移距的變化a 反射波; b 透射波

圖6分別展示了反射波與透射波相對誤差隨PPW的變化,可以看到不同介質(zhì)離散方法對P波影響很大。直接取值法具有極大的不穩(wěn)定性,其對于波形的誤差顯著大于其它方法。就不同波形而言,相同情況下透射波的準確性普遍高于反射波。圖7展示了反射波與透射波分別在相對誤差約為10%與5%時的波形對比。

3.2　P-S轉(zhuǎn)換波誤差隨網(wǎng)格的變化

圖8展示了P-S轉(zhuǎn)換波相對誤差隨PPW的變化。由于界面上發(fā)生了P-S波的轉(zhuǎn)換,因此轉(zhuǎn)換波(反射)的精度顯著低于其它波形。由于直接取值法與格點平均法不能精確描述界面在網(wǎng)格間的具體位置,因此它們在粗網(wǎng)格情況下誤差較大,而通過增大PPW提高模擬精度的計算代價高昂。P-S轉(zhuǎn)換波在相對誤差約為10%與5%時的波形對比見圖9。

圖6　相對誤差隨PPW的變化a 反射波; b 透射波

圖7　相對誤差波形對比(10%與5%)a 反射波波形; b 透射波波形

圖8　P-S轉(zhuǎn)換波相對誤差隨PPW的變化a 反射波; b 透射波

圖9　P-S轉(zhuǎn)換波相對誤差波形對比(10%和5%)a 反射波波形; b 透射波波形

3.3　多次反射波誤差隨網(wǎng)格的變化

為了觀察多次反射波相對誤差隨網(wǎng)格的變化,設置一個平層模型,模型大小為nx×nz=2500m×1500m,包括3層介質(zhì):頂層與底層縱波速度為4000m/s,橫波速度為2000m/s,密度2050kg/m3;中間低速層縱波速度為2000m/s,橫波速度為1000m/s,密度1410kg/m3。震源(700m,800m)與檢波點(1330m,770m)均位于模型中間的低速層。圖10展示了多次反射波在不同PPW情況下,使用不同介質(zhì)離散方法得到的相對誤差。圖11和圖12展示了誤差約10%和5%時的波形對比(實際誤差分別為9.61%與4.79%),可以看到在相對誤差約為5%時,波形的擬合已經(jīng)相對較好。

圖10　多次反射波相對誤差隨PPW的變化

圖11　多次反射波相對誤差波形對比

表1與表2分別展示了4種介質(zhì)離散方法達到10%與5%的相對誤差時,不同類型反射/透射波各自的PPW要求。

圖12　多次反射波相對誤差波形對比局部(圖11中1.0～1.6s波形)

個

表2　不同反射/透射波PPW約束(相對誤差5%)　個

4　結(jié)束語

本文針對介質(zhì)的彈性系數(shù)處理對于交錯網(wǎng)格波場模擬準確性造成的影響進行了分析,給出了不同介質(zhì)離散方式下各個反射/透射波相對誤差的PPW約束參考,得到以下結(jié)論：

1) 在使用交錯網(wǎng)格有限差分進行正演計算時,不僅要考慮格式本身的色散誤差,還需要考慮介質(zhì)離散方式對模擬精度的影響,特別是在復雜模型中,即使在滿足色散誤差的情況下,仍然存在較為嚴重的界面誤差,這類誤差對整體波形有著不可忽視的影響。

2) 界面誤差對PPW的要求普遍大于4階格式的色散誤差,因此,在選取交錯網(wǎng)格的網(wǎng)格步長時,為了保證誤差在一定范圍之內(nèi),應該以橫波最小速度對應的PPW為約束,在可以接受的誤差范圍內(nèi)首先確定網(wǎng)格的最大空間步長。

在實際交錯網(wǎng)格計算中,首先應該確定計算精度要求,根據(jù)所采取的介質(zhì)參數(shù)處理算法,參考界面誤差隨PPW的變化關(guān)系,以最小橫波速度為基準,確定網(wǎng)格的最大空間步長,再根據(jù)柯朗-弗里德里希斯-列維(CFL)穩(wěn)定性條件確定時間步長,從而在確保計算準確性的前提下,減少計算資源需求。

[1]VIRIEUX J.P-SV wave propagation in heterogeneous media:velocity-stress finite-difference method[J].Geophysics,1986,51(4):889-901

[2]LEVANDER A R.Fourth-order finite-difference P-SV seismograms[J].Geophysics,1988,53(11):1425-1436

[3]MCGARRY R,PASALIC D,ONG C.Anisotropic elastic modeling on a Lebedev grid:dispersion reduction and grid decoupling[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2829-2833

[4]LISITSA V,VISHNEVSKIY D.Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity[J].Geophysical Prospecting,2010,58(4):619-635

[5]SAENGER E H,GOLD N,SHAPIRO S A.Modeling the propagation of elastic waves using a modified finite-difference grid[J].Wave Motion,2000,31(1):77-92

[6]裴正林.三維各向同性介質(zhì)彈性波方程交錯網(wǎng)格高階有限差分法模擬[J].石油物探,2005,44(4):308-315

PEI Z L.Numerical simulation of elastic wave equation in 3-D isotropic media with staggered-grid high-order difference method[J].Geophysical Prospecting for Petroleum,2005,44(4):308-315

[7]李勝軍,孫成禹,高建虎,等.地震波數(shù)值模擬中的頻散壓制方法分析[J].石油物探,2008,47(5):444-448

LI S J,SUN C Y,GAO J H,et al.Analysis of dispersion suppression in wave equation numerical simulation[J].Geophysical Prospecting for Petroleum,2008,47(5):444-448

[8]陳可洋.地震波數(shù)值模擬中差分近似的各向異性分析[J].石油物探,2010,49(1):19-22

CHEN K Y.Anisotropic analysis of difference approximation in seismic wave numerical modeling[J].Geophysical Prospecting for Petroleum,2010,49(1):19-22

[9]梁文全,王彥飛,楊長春.基于優(yōu)化方法的時間-空間域隱格式有限差分算子確定方法[J].石油物探,2015,54(3):254-259

LIANG W Q,WANG Y F,YANG C C.Determination on the implicit finite-difference operator based on optimization method in time-space domain[J].Geophysical Prospecting for Petroleum,2015,54(3):254-259

[10]陳東,梁文全,辛維.適用于聲波方程數(shù)值模擬的時間-空間域隱式有限差分算子優(yōu)化方法[J].地球物理學報,2016,59(4):1491-1497

CHEN D,LIANG W Q,XIN W.Acoustic wave equation modeling based on implicit finite difference operators in the time-space domain[J].Chinese Journal of Geophysics,59(4):1491-1497

[11]ALTERMAN Z,KARAL F C.Propagation of elastic waves in layered media by finite difference methods[J].Bulletin of the Seismological Society of America,1968,58(1):367-398

[12]VIRIEUX J.SH-wave propagation in heterogeneous media:velocity-stress finite-difference method[J].Geophysics,1984,49(11):1933-1942

[13]GRAVES R W.Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[J].Bulletin of the Seismological Society of America,1996,86(4):1091-1106

[14]MUIR F,DELLINGER J,ETGEN J,et al.Modeling elastic fields across irregular boundaries[J].Geophysics,1992,57(9):1189-1193

[15]BACKUS G E.Long-wave elastic anisotropy produced by horizontal layering[J].Journal of Geophysical Research,1962,67(11):4427-4440

[17]MOCZO P,KRISTEK J,GLIS M.The finite-difference modelling of earthquake motions:waves and ruptures[M].Cambridge:Cambridge University Press,2014:199-226

[18]KRISTEK J,MOCZO P,CHALJUB E,et al.An orthorhombic representation of a heterogeneous medium for the finite-difference modelling of seismic wave propagation[J].Geophysical Journal International,2017,208(2):1250-1264

[19]LISITSA V,PODGORNOVA O,TCHEVERDA V.On the interface error analysis for finite difference wave simulation[J].Computational Geosciences,2010,14(4):769-778

[20]VISHNEVSKY D,LISITSA V,TCHEVERDA V,et al.Numerical study of the interface errors of finite-difference simulations of seismic waves[J].Geophysics,2014,79(4):T219-T232

[21]CHEN X.Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method I,theory of two-dimensional SH case[J].Bulletin of the Seismological Society of America,1990,80(6A):1696-1724

[22]CHEN X.Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method Ⅱ,applications for 2D SH case[J].Bulletin of the Seismological Society of America,1995,85(4):1094-1106

[23]CHEN X.Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method III,theory of 2D P-SV case[J].Bulletin of the Seismological Society of America,1996,86(2):389-405

[24]SCHOENBERG M,MUIR F.A calculus for finely layered anisotropic media[J].Geophysics,1989,54(5):581-589