介質離散方法對地震波場有限差分數值模擬準確性的影響

朱　強,姜蘆倩,張　偉

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

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

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

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

1　交錯網格有限差分方法

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

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

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

(3)

式中:h表示網格步長,αm表示差分系數。則交錯網格格式如下:

(4)

2　介質參數離散方法

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

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

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

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

2.1　直接取值法

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

2.2　格點平均法

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

2.3　體積分平均法

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

設界面位于x=0處,在一維情況下對波動方程在單位網格內做空間積分:

(7)

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

(8)

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

(9)

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

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

2.4　正交各向異性等效介質法

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

二維各向同性介質情況下:

(13)

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

(14)

其中,

(15)

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

(16)

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

(17)

(18)

設M=λ+2μ,則:

(19)

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

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

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

(20)

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

于Ricker子波,其PPW定義為:

(21)

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

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

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

3.1　透射波和反射波相對誤差隨網格的變化

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

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

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

3.2　P-S轉換波誤差隨網格的變化

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

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

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

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

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

3.3　多次反射波誤差隨網格的變化

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

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

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

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

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

個

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

4　結束語

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

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

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

在實際交錯網格計算中,首先應該確定計算精度要求,根據所采取的介質參數處理算法,參考界面誤差隨PPW的變化關系,以最小橫波速度為基準,確定網格的最大空間步長,再根據柯朗-弗里德里希斯-列維(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]裴正林.三維各向同性介質彈性波方程交錯網格高階有限差分法模擬[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]李勝軍,孫成禹,高建虎,等.地震波數值模擬中的頻散壓制方法分析[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]陳可洋.地震波數值模擬中差分近似的各向異性分析[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]陳東,梁文全,辛維.適用于聲波方程數值模擬的時間-空間域隱式有限差分算子優(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