譚文卓 吳幫玉 李 博 雷 軍

(①西安交通大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，陜西西安 710049； ②中國石化石油物探技術(shù)研究院，江蘇南京 211103； ③中國石油長慶油田分公司第六采氣廠，陜西西安 710018)

0 引言

地震波場數(shù)值模擬是復(fù)雜地區(qū)地震勘探常用手段[1]，尤其在逆時偏移成像[2]、全波形反演[3]等領(lǐng)域。這些模擬方法可分為波動方程法和幾何射線法兩種主要類型。相較于幾何射線法，波動方程法模擬的地震波場包含了更多的地震波傳播信息，能為地震波傳播機理研究及后期地質(zhì)綜合解釋提供更多基礎(chǔ)數(shù)據(jù)。因此，波動方程法在地震波場數(shù)值模擬中得到越來越多的應(yīng)用。

波動方程數(shù)值解法主要有有限差分法[4-5]、偽譜法[6-7]、有限元法[8-10]等。其中，基于快速傅里葉變換的偽譜法對于帶限波場的正演可達到無限階差分算子精度[11]，被廣泛應(yīng)用于復(fù)雜介質(zhì)的波場正演模擬與偏移成像。傳統(tǒng)的波動方程求解方法往往采用規(guī)則的均勻網(wǎng)格對模型進行離散，以模型的最小速度決定網(wǎng)格尺寸，在高速區(qū)易產(chǎn)生“過采樣”現(xiàn)象。而變網(wǎng)格方法能根據(jù)精度要求，在模型速度較低區(qū)域采用細網(wǎng)格、在速度較高區(qū)域采用粗網(wǎng)格，以便兼顧波場模擬的效率與精度。

Chen等[12]提出基于坐標(biāo)變換的梯形網(wǎng)格方法，使用梯形網(wǎng)格對介質(zhì)進行剖分，并通過坐標(biāo)變換在新坐標(biāo)系下采用偽譜法對彈性波波場做數(shù)值模擬。該方法計算效率高，同時避免了一般變網(wǎng)格方法易造成的虛假反射。

高靜懷等[13]最近提出基于深度域均勻采樣的梯形網(wǎng)格有限差分模擬法，應(yīng)用橫向采樣間隔隨深度增加而逐漸線性增大、縱向采樣間隔恒定的梯形網(wǎng)格剖分模型，并采用有限差分法求解梯形坐標(biāo)系波動方程，模擬的地震波場與實際情況相符。但由于在深度域進行均勻采樣，對于含有高速夾層的模型，仍會出現(xiàn)一定程度過采樣。

Wu等[14]對上述梯形網(wǎng)格剖分方法[13]做了改進，在深度方向進行非均勻采樣，進一步減少了網(wǎng)格數(shù)量。通過Marmousi模型試算顯示，梯形網(wǎng)格所占內(nèi)存只有常規(guī)網(wǎng)格的25%，同時計算時間減少了66%，大大提高了計算效率。

本文基于Wu等[14]的研究，提出一種梯形網(wǎng)格偽譜法地震波模擬方法。內(nèi)容包括：定義與速度相關(guān)的梯形坐標(biāo)變換； 推導(dǎo)梯形坐標(biāo)系下波動方程表達式，同時為消除人工截斷邊界反射波的影響，定義了相關(guān)的吸收邊界； 研究時間差分的穩(wěn)定性條件； 利用典型模型測試數(shù)值結(jié)果，并與常規(guī)網(wǎng)格偽譜法和有限差分法結(jié)果進行對比； 最后給出所得結(jié)論。

1 基本理論

1.1 梯形坐標(biāo)變換

根據(jù)Wu等[14]的理論，可定義如下直角坐標(biāo)系與梯形坐標(biāo)系之間的變換關(guān)系

(1)

式中： (x,z)為梯形坐標(biāo)系坐標(biāo)， (x0,z0)為對應(yīng)的直角坐標(biāo)系坐標(biāo)；α為橫向位置參數(shù)，一般與震源橫向位置相同；γ為形態(tài)參數(shù)，與模型速度隨深度變化的程度相關(guān)；g(z)為深度域映射函數(shù)。

由式(1)可看出，梯形坐標(biāo)系中在矩形網(wǎng)格的橫向采樣間隔確定為Δx時，對應(yīng)直角坐標(biāo)系中梯形網(wǎng)格在深度z0處橫向采樣間隔為

Δx0(z0)=(1+γz0)Δx

(2)

得到Δx后，設(shè)梯形坐標(biāo)系矩形網(wǎng)格的縱向采樣間隔Δz=Δx，則對于深度采樣函數(shù)g(z)，可用以下迭代方法得到

(3)

式中NGz為一個波長內(nèi)的最小縱向采樣點數(shù)。

通過三次樣條插值得到g(z)的一階、二階導(dǎo)數(shù)g′(z)與g″(z)。最后，可得到形態(tài)參數(shù)

(4)

式中z0 max為給定模型的最大深度。

圖1為梯形網(wǎng)格剖分示意圖。

圖1 梯形網(wǎng)格剖分示意圖

1.2 帶PML吸收邊界的聲波方程的梯形坐標(biāo)系表征

由式(1)通過鏈?zhǔn)椒▌t可求得兩個坐標(biāo)系之間一階、二階偏導(dǎo)的關(guān)系

(5)

(6)

(7)

(8)

(9)

將式(9)代入空間混合偏導(dǎo)項，即得

(10)

根據(jù)完全匹配層吸收邊界條件相關(guān)理論[16-19]，可將波場u分裂為三項u1、u2、u3，在計算區(qū)域外部引入完全匹配層。直角坐標(biāo)系中帶PML吸收邊界的二階聲波方程分以下三種情形表示。

在內(nèi)部計算區(qū)域

(11)

在左側(cè)和右側(cè)PML區(qū)域

(12)

在上側(cè)和下側(cè)PML區(qū)域

(13)

式中： (x0s,z0s)為震源在直角坐標(biāo)系中位置；衰減函數(shù)d(i)及其一階導(dǎo)數(shù)d′(i)可表示為

其中：i=x0或z0；L為PML吸收層的厚度；si為PML層內(nèi)節(jié)點到計算區(qū)域邊界的距離；R為衰減系數(shù)，一般取10-3。

將式(5)～式(8)代入式(11)～式(13)，即得梯形坐標(biāo)系中帶PML吸收邊界的三種情形的二階聲波方程表達式。

在內(nèi)部計算區(qū)域

(14)

在左側(cè)和右側(cè)PML區(qū)域

(15)

在上側(cè)和下側(cè)PML區(qū)域

(16)

式中(xs,zs)為震源在梯形坐標(biāo)系中位置。

此處的衰減函數(shù)d(i)及其一階導(dǎo)數(shù)d′(i)的表達式在形式上與直角坐標(biāo)系時相同，且i=x或z。

1.3 傅里葉偽譜法原理

由傅里葉變換的相關(guān)理論[20]，若函數(shù)f(x)在Hilbert空間中可積，則有

(17)

式中F(k)為f(x)的傅里葉變換。

以均勻網(wǎng)格離散f(x)，則式(17)變?yōu)?/p>

n=0,…,N-1

(18)

通過上述方法就能得到f(x)在各個離散點的導(dǎo)數(shù)值。對于波場函數(shù)，只需改變相關(guān)參數(shù)即可計算其偏導(dǎo)值。在數(shù)值模擬中常采用快速傅里葉變換加速計算。

1.4 穩(wěn)定性條件

與z軸夾角為θ的平面波的解析式為

(19)

不帶PML邊界的梯形坐標(biāo)系偽譜法求解聲波方程的差分格式為

(20)

將式(19)代入式(20)，可得

(21)

(22)

(23)

(24)

2 數(shù)值算例

利用不同速度模型脈沖響應(yīng)驗證本文方法的正確性，并與常規(guī)網(wǎng)格偽譜法及常規(guī)網(wǎng)格有限差分法結(jié)果進行對比。測試均基于Intel(R) Core(TM) i5-6300HQ CPU及MATLAB編程。其中常規(guī)網(wǎng)格偽譜法與常規(guī)網(wǎng)格有限差分法均采用中心差分近似時間偏導(dǎo)，常規(guī)網(wǎng)格有限差分法采用八階差分算子近似空間偏導(dǎo)[21-22]。由于梯形坐標(biāo)系偽譜法可在大尺度時間步長下保持良好穩(wěn)定性，因此采用Erik等[23-24]提出的時域離散變換消除大尺度時間步長帶來的時間頻散。

2.1 均勻介質(zhì)模型

在一個速度為3000m/s、深度為2000m、寬度為1000m的梯形坐標(biāo)系均勻介質(zhì)模型中，設(shè)置橫向、縱向采樣間隔均為25m，形態(tài)參數(shù)γ=5×10-4。在模型中心放置主頻為10Hz的Ricker子波作為震源，震源時延為0.075s。使用時間步長為1ms的梯形坐標(biāo)系偽譜法做波場模擬。得到t=0.4s時的波場快照(圖2a)。對該波場快照做坐標(biāo)變換即可得到直角坐標(biāo)系波場快照。但由于梯形網(wǎng)格剖分的空間采樣點數(shù)較少，因而得到的直角坐標(biāo)系波場快照會缺失部分波場信息，可通過對直角坐標(biāo)系波場快照插值得到更多波場信息(圖2b)。圖3顯示由梯形坐標(biāo)系中位于(500m，0m)處的檢波器接收到的信號插值得到的直角坐標(biāo)系(1000m，0)處的波場時間曲線，可見梯形網(wǎng)格偽譜法求得的數(shù)值解與解析解[25]高度吻合，驗證了本文方法的正確性。

圖2 均勻介質(zhì)模型的梯形坐標(biāo)系(a)和直角坐標(biāo)系(b)波場快照

圖3 均勻介質(zhì)模型(1000m，0)處波場時間曲線

2.2 三層速度模型

測試選用圖4所示的2000m×2000m的三層速度模型。各層的速度、厚度(從上到下)依次為2000、4000、3000m/s，500、500、1000m。梯形網(wǎng)格偽譜法的實際模擬范圍是圖中紅線包圍的(梯形)區(qū)域。設(shè)置直角坐標(biāo)系梯形網(wǎng)格的橫向采樣間隔由最頂層的10.5470m增至最底層的18.5185m，深度域映射函數(shù)如圖5所示，形態(tài)參數(shù)為γ=3.697×10-4。

圖4 三層速度模型

在模型頂層正中央放置主頻為20Hz的Ricker子波作為震源，使用時間步長為1ms的梯形坐標(biāo)系偽譜法做波場模擬，得到t1=0.18s、t2=0.42s、t3=0.81s三個時刻梯形坐標(biāo)系波場快照(圖6a～圖6c)。經(jīng)過坐標(biāo)變換和插值，得到相應(yīng)的直角坐標(biāo)系波場快照(圖6d～圖6f)。同時，在相同的三層速度模型中使用時間步長為1ms，橫向、縱向采樣間隔均為10m的常規(guī)網(wǎng)格偽譜法和常規(guī)網(wǎng)格有限差分法得到同樣三個時刻的波場快照(圖6g～圖6i、圖6j～圖6l)?？梢娞菪尉W(wǎng)格因未覆蓋左右兩個角形區(qū)域，使用梯形網(wǎng)格偽譜法求得的波場快照會缺失部分淺層大角度信息，但在其他區(qū)域與常規(guī)網(wǎng)格方法結(jié)果一致。

2.3 Marmousi模型

選取的Marmousi模型(圖7a)中紅色線圈定范圍是梯形網(wǎng)格偽譜法的實際模擬區(qū)域。設(shè)置梯形網(wǎng)格的橫向采樣間隔由最淺層的8.7184m增至最底層的17.3493m，深度域映射函數(shù)如圖7c所示，可算出形態(tài)參數(shù)γ=3.3003×10-4。在模型頂層正中間放置主頻為15Hz的Ricker子波震源，震源時延為0.125s。使用時間步長為0.943ms的梯形網(wǎng)格偽譜法做波場模擬，得到t1=1s、t2=2s時的梯形坐標(biāo)系波場快照(圖8a、圖9a)。

圖5 梯形坐標(biāo)系深度域映射函數(shù)

圖6 三層速度模型波場快照

經(jīng)過坐標(biāo)變換和插值，得到t1=1s、t2=2s時的直角坐標(biāo)系波場快照(圖8b、圖9b)。在Marmousi原始速度模型中使用時間步長為0.436ms，橫向采樣間隔為12.5m，縱向采樣間隔為4m的常規(guī)網(wǎng)格偽譜法做波場模擬，得到t1=1s、t2=2s時的波場快照(圖8c、圖9c)。使用時間步長為0.537ms，橫向采樣間隔為12.5m，縱向采樣間隔為4m的常規(guī)網(wǎng)格有限差分法進行波場模擬，得到的t1=1s、t2=2s時的波場快照(圖8d、圖9d)。可見三種方法在Marmousi模型中正演得到的同一時刻波場快照并無差別。為了對結(jié)果做更詳細對比，抽取了三種方法得到的t1=1s時波場快照的第269、第369及第469道的深度域波形(圖10)，可見三者波形高度一致，充分體現(xiàn)了梯形網(wǎng)格偽譜法的高精度。

表1列出使用三種方法進行地震波數(shù)值模擬時各項參數(shù)?？煽吹较噍^于常規(guī)網(wǎng)格剖分，梯形網(wǎng)格剖分的網(wǎng)格數(shù)減少了69.26%； 相較于常規(guī)網(wǎng)格偽譜法，梯形網(wǎng)格偽譜法的計算時間減少了58.16%； 相較于常規(guī)網(wǎng)格有限差分法，梯形網(wǎng)格偽譜法的計算時間減少了60.33%。考慮到梯形網(wǎng)格偽譜法的實際模擬區(qū)域比常規(guī)方法減少約25%，因此對于同一模擬區(qū)域，梯形網(wǎng)格偽譜法的計算時間比常規(guī)方法只減少約45%。同時，三種方法的精度相當(dāng)。這充分展現(xiàn)了深度非均勻采樣梯形網(wǎng)格偽譜法進行地震波模擬時的高效率和高精度。

圖7 Marmousi模型與深度域映射函數(shù)

(a)Marmousi原始模型，紅線框為梯形網(wǎng)格偽譜法實際模擬區(qū)域； (b)梯形坐標(biāo)系Marmousi模型； (c)Marmousi模型深度域映射函數(shù)

圖8 Marmousi模型t=1s時的波場快照

圖9 Marmousi模型t=2s時的波場快照

圖10 Marmousi模型第269道(a)、第369道(b)及第469道(c)波形對比圖

表1 三種方法Marmousi模型地震模擬參數(shù)

3 結(jié)論

本文提出了一種梯形網(wǎng)格偽譜法地震波模擬方法。相較于常規(guī)均勻網(wǎng)格方法，梯形網(wǎng)格更契合地震波傳播速度由淺入深逐漸增大的物理現(xiàn)象，在低速區(qū)域采用細網(wǎng)格、高速區(qū)域采用粗網(wǎng)格，同時在深度域采用非均勻采樣，這樣可大幅度減少內(nèi)存需求，提高計算效率。對均勻介質(zhì)模型的數(shù)值測試表明，本文所提方法符合解析解，同時數(shù)值頻散較小。三層速度模型的測試結(jié)果表明，本文方法正演得到的地震波場與常規(guī)方法正演得到的地震波場相比，除了部分淺層大角度信息缺失，無其他區(qū)別。針對Marmousi模型的測試表明，梯形網(wǎng)格偽譜法在復(fù)雜模型中同樣可得很好模擬效果。相較于常規(guī)網(wǎng)格剖分，梯形網(wǎng)格剖分能減少約69%的采樣點數(shù)； 相較于常規(guī)網(wǎng)格偽譜法，梯形網(wǎng)格偽譜法計算時間能減少約58%； 相較于常規(guī)網(wǎng)格有限差分法，梯形網(wǎng)格偽譜法的計算時間能減少約60%。上述測試結(jié)果表明梯形網(wǎng)格偽譜法地震波模擬方法是一種高效、高精度地震波模擬方法，具有較高應(yīng)用價值。

當(dāng)然，本文方法也存在一定的局限性： ①梯形網(wǎng)格并未覆蓋左右兩個角形區(qū)域，因此使用梯形網(wǎng)格偽譜法做波場模擬會缺失此部分波場信息； ②由于采用梯形網(wǎng)格剖分只擬合地震波傳播速度由淺入深逐漸增加的整體趨勢，因此對于淺層高速或深層低速模型仍會出現(xiàn)“過采樣”。