黏聲波動方程變機制數(shù)有限差分正演

蔡瑞乾 孫成禹* 伍敦仕 李世中

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ②青島海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東青島 266071； ③中國石油勘探開發(fā)研究院西北分院，甘肅蘭州 730020)

0 引言

實際地球介質(zhì)普遍具有黏彈性，地震波在實際介質(zhì)中傳播會產(chǎn)生衰減。地球介質(zhì)的這種衰減特性可以用黏彈模型來描述[1]。

對于黏彈介質(zhì)的地震波場模擬，國內(nèi)外學(xué)者已經(jīng)做了大量研究。Liu等[2]建立了線性黏彈流變學(xué)的理論框架，可以很好地解釋實驗觀測到的地震波在地球介質(zhì)中傳播時的衰減現(xiàn)象；Day等[3]首次利用Padé近似方法將時間域的應(yīng)力—應(yīng)變褶積關(guān)系轉(zhuǎn)換成微分形式，在二維時間域進行黏彈性介質(zhì)數(shù)值模擬；Emmerich等[4]認(rèn)為Padé近似方法存在計算效率低和數(shù)值模擬效果較差的缺點，在流變學(xué)模型的基礎(chǔ)上提出了廣義標(biāo)準(zhǔn)線性固體(Generalized Standard Linear Solid，GSLS)模型，并將其應(yīng)用于標(biāo)量波動方程的二維有限差分正演；Robertsson等[5-6]基于GSLS模型，推導(dǎo)了交錯網(wǎng)格速度—應(yīng)力有限差分格式，并對二維和三維黏彈介質(zhì)進行了波場數(shù)值模擬與分析。Blanch等[7]針對時間域GSLS模型有限差分正演內(nèi)存需求大、計算時間長的問題提出了τ方法，顯著降低了存儲空間要求，提高了計算效率； Bohlen[8]提出三維并行黏彈有限差分正演方法，使三維地震波場模擬的效率和規(guī)模大大提高； Saenger等[9]將旋轉(zhuǎn)交錯網(wǎng)格有限差分法用于模擬黏彈性介質(zhì)中波的傳播。杜啟振等[10-11]先后采用有限元法、偽譜法和有限差分法對方位各向異性黏彈性介質(zhì)中的波場進行模擬，分析了速度頻散和衰減特征； 牛濱華[12]詳細(xì)介紹了黏彈性介質(zhì)的理論及地震波在黏彈性介質(zhì)中的傳播特征；嚴(yán)紅勇等[13]給出了適用于黏彈TTI介質(zhì)的二維三分量一階速度—應(yīng)力方程，并采用旋轉(zhuǎn)交錯網(wǎng)格高階有限差分實現(xiàn)了黏彈TTI介質(zhì)的地震波場數(shù)值模擬；何兵紅等[14]推導(dǎo)了基于傅里葉有限差分(FFD)算子的衰減介質(zhì)地震單程波波場傳播算子，通過空間—時間域與波數(shù)—頻率域之間的轉(zhuǎn)化實現(xiàn)了衰減介質(zhì)地震波場數(shù)值模擬；羅文山等[15]基于常Q模型的應(yīng)力—應(yīng)變關(guān)系推導(dǎo)了時間域近似常Q黏滯波動方程的一階速度—壓力形式，能較好地描述地震波在黏滯介質(zhì)中的衰減和頻散； 姚振岸等[16]對黏彈各向異性介質(zhì)中不同震源機制的微地震波場特征進行了研究。

關(guān)于黏彈介質(zhì)的研究還有很多[17-19]，但這些研究主要是針對更高維度或更復(fù)雜的介質(zhì)兩種情況，而對黏彈機制數(shù)的研究較少。對于GSLS模型有限差分正演，通常使用的黏彈單元(常稱為黏彈機制)數(shù)越多，精度越高，但相應(yīng)的求解難度也越高，計算效率也越低。一般情況下，當(dāng)Q值較小時需要使用較多的機制數(shù)才能得到較好的模擬效果，當(dāng)Q值較高時僅使用較少的機制數(shù)便能得到較好的模擬效果。Zhu等[20]對GSLS模型正演中單機制數(shù)和三機制數(shù)的精確度做了調(diào)研，結(jié)果表明介質(zhì)吸收衰減較弱時僅使用單機制數(shù)模擬就能達到較高的精度，介質(zhì)吸收衰減較強而波傳播距離很短時單機制數(shù)也能得到較好的模擬精度，因此一般情況下使用單機制模擬便可獲得較好的結(jié)果，無需使用更多的黏彈機制數(shù)。但隨著業(yè)界對精度的要求日益提高和近年來對近地表強吸收衰減影響的認(rèn)識逐步加深[21-23]，單機制的模擬精度已不能滿足需求。Cai等[24]對不同機制數(shù)的模擬精度做了調(diào)研，并分析了變機制數(shù)正演的可行性。

本文基于GSLS模型有限差分正演提出了一種黏聲波動方程變機制數(shù)有限差分正演方法，首先研究了不同機制數(shù)GSLS模型的地震波場正演模擬精度，然后對變機制數(shù)有限差分的適用性做了評估，最后實現(xiàn)地震波場交錯網(wǎng)格有限差分?jǐn)?shù)值模擬并對模擬結(jié)果進行了分析。

1 黏彈模型及其波動方程

1.1 GSLS衰減模型

廣義標(biāo)準(zhǔn)線性固體模型可以很好地模擬介質(zhì)的衰減特性，圖1為L個Maxwell體(一個彈性元和一個牛頓黏滯元串聯(lián)而成)與一彈性元并聯(lián)組成的GSLS模型示意圖。GSLS模型復(fù)模量可表示為[8]

(1)

式中：ω為角頻率；K0為與Maxwell體并聯(lián)的彈性元的彈性模量；Ki和ηi(l=1，2，…，L)分別為Maxwell體的彈性模量和牛頓黏滯系數(shù)；τσ,l和τε,l分別為第l個Maxwell體的應(yīng)力松弛時間和應(yīng)變松弛時間，可表示為[25]

(2)

(3)

圖1 GSLS衰減模型示意圖

品質(zhì)因子Q定義為[26]

(4)

式中RE(·)和IM(·)分別表示求復(fù)數(shù)的實部和虛部。由式(1)和式(4)可得

(5)

式中

(6)

當(dāng)GSLS模型所用機制數(shù)L已知時，參數(shù)τσ,l、τ的最優(yōu)值可由最小二乘反演得到

(7)

(8)

其中

(9)

1.2 二維黏聲波動方程

有關(guān)黏彈波動方程的推導(dǎo)Robertsson等[6]已做過詳細(xì)描述。二維黏聲微分方程為

(10)

式中：變量上面的點表示時間導(dǎo)數(shù)；P為壓力；MR為介質(zhì)松弛模量；rl為第l個Maxwell體的記憶變量；vx和vz分別為x方向和z方向的速度分量；ρ為密度；fx和fz分別為x方向和z方向施加的外力。

2 變機制數(shù)有限差分正演精度分析

2.1 數(shù)值解與解析解的對比

為評估不同機制數(shù)黏聲方程有限差分波場模擬的精度，本文首先計算二維均勻各向同性介質(zhì)中黏聲格林函數(shù)的解析解[27]，然后將不同機制數(shù)黏聲方程的模擬結(jié)果與解析解做對比。

對于GSLS模型而言，使用的機制數(shù)越多，其對黏彈介質(zhì)的描述就越準(zhǔn)確。為研究不同機制數(shù)的模擬精度，建立一個均勻模型，Q值設(shè)置為100，震源是主頻為30Hz為Ricker子波，頻帶為5～100Hz，參考頻率(主頻)處速度為3000m/s，密度為2.1g/cm3。采用時間二階、空間八階差分格式進行正演模擬，利用PML邊界條件吸收邊界反射。使用不同機制個數(shù)模擬的Q值曲線和相速度曲線如圖2所示。由圖2可見，單機制數(shù)(L=1)模擬的Q值只在參考頻率附近與理論值較為接近，在低頻和高頻端都遠(yuǎn)大于理論值，三機制數(shù)(L=3)模擬的Q值在給定頻帶上都與理論值相差不大，五機制數(shù)(L=5)模擬的Q值在給定頻帶上與理論值幾乎完全重合，表明GSLS模型使用的機制數(shù)越多其對黏彈介質(zhì)的模擬就越準(zhǔn)確。相速度曲線(圖2b)與Q值模擬曲線類似，即使用的機制數(shù)越多其模擬的相速度曲線越接近理論曲線。

圖2 Q=100時不同機制數(shù)的模擬結(jié)果

圖3為不同機制數(shù)黏聲方程模擬結(jié)果隨傳播距離的變化?？梢钥闯觯寒?dāng)傳播距離較近(500m)時，不同機制數(shù)的數(shù)值模擬結(jié)果都能與解析解相吻合； 中傳播距離(2500m)時單機制數(shù)和三機制數(shù)的模擬結(jié)果出現(xiàn)微小的誤差； 遠(yuǎn)傳播距離(4500m)時各機制數(shù)的模擬結(jié)果都出現(xiàn)誤差，但機制數(shù)越大誤差越小。為衡量不同機制數(shù)波場模擬的準(zhǔn)確性，計算模擬結(jié)果的均方根誤差(RMSE)

(11)

表1 不同機制數(shù)模擬結(jié)果均方根誤差(Q=100)

圖3 Q=100時不同機制數(shù)有限差分模擬結(jié)果隨傳播距離的變化

2.2 機制數(shù)的選取

實際介質(zhì)模擬中，模擬精度與機制數(shù)、Q值和波的傳播距離有關(guān)。由上文可知基于GSLS模型的黏聲波動方程有限差分正演使用的機制數(shù)越大，對地震波場的模擬就越準(zhǔn)確，但相應(yīng)的波動方程也越復(fù)雜，需要的存儲空間也越多，計算效率也越低。為達到計算效率與計算精度的統(tǒng)一，需要根據(jù)所需要的模擬精度、Q值和波在該Q值區(qū)域傳播的距離確定模擬時使用的機制數(shù)。具體研究思路為：

(1)采用的不同機制個數(shù)對不同Q值、不同傳播距離情況進行正演模擬；

(2)將不同機制數(shù)的模擬結(jié)果分別與相應(yīng)的解析解作對比，計算得出對應(yīng)的模擬精度；

(3)將不同機制數(shù)模擬結(jié)果的模擬精度與所需要的最低模擬精度進行比較，擬合得到不同機制數(shù)的適用范圍。

本文模擬精度設(shè)置為95%；Q值從2變化到200，涵蓋強吸收介質(zhì)和弱吸收介質(zhì)的情況； 波的傳播距離從100m變化到10000m。震源設(shè)置為Ricker子波，主頻為30Hz，頻帶為5～100Hz，參考頻率處(主頻)速度為3000m/s，密度為2.1g/cm3，采用時間二階、空間八階差分格式進行正演模擬，邊界設(shè)置為PML條件。圖4展示了L=1、3、5時的模擬精度。由圖可知三種機制數(shù)的模擬精度都隨著傳播距離減小而增大、隨Q值增大而增大。L=1時的模擬精度在傳播距離較大且Q值較小的時候出現(xiàn)負(fù)數(shù)情況，說明此時模擬的波場與實際波場相比出現(xiàn)了極大的波形畸變；L=3時的模擬精度在全區(qū)域大于50%并且大部分區(qū)域模擬精度在95%以上；L=5時的模擬精度在全區(qū)域大于70%并且絕大部分區(qū)域模擬精度都在95%以上。另外，模擬精度在傳播距離和Q值組成的平面上的投影坐標(biāo)值表示波傳播一定的距離時滿足特定精度所需的最小Q值，因此取不同機制數(shù)在模擬精度為95%時的投影點即可擬合得到不同機制數(shù)的模擬精度與Q值和波的傳播距離的關(guān)系，如圖5所示。L=1、3和5滿足模擬精度為95%時,最小Q值和波的傳播距離的關(guān)系分別為

Qmin=30.937lnD-157.2L=1

(12)

Qmin=7.6039lnD-40.364L=3

(13)

Qmin=5.7302lnD-37.40L=5

(14)

式中：Qmin為最小Q值；D為波的傳播距離。因此在實際介質(zhì)模擬中，控制精度為95%，在模型已知的情況下便可直接確定機制數(shù)。此外由圖5還可看出，當(dāng)傳播距離足夠遠(yuǎn)、Q值又足夠小時，L=5也不能滿足模擬精度要求，此時應(yīng)當(dāng)使用更大的機制數(shù)，但此種情況實際中極其少見。

圖4 模擬精度隨Q值和傳播距離變化

圖5 不同機制數(shù)模擬精度為95%時所需最小Q值與波的傳播距離間的關(guān)系

2.3 變機制數(shù)引起的數(shù)值干擾分析

黏聲波動方程變機制數(shù)有限差分正演在不同Q值區(qū)域會使用不同機制數(shù)，因此在Q值發(fā)生變化的邊界上，不同機制數(shù)黏聲方程差分格式的差異不可避免地會引起人為的數(shù)值干擾。因此需對這種數(shù)值干擾進行定量分析，從而更好地評價方法的適用性。

設(shè)計一個尺寸為1000m×1000m的均勻模型，空間采樣間隔Δx=Δz=5m，速度為2000m/s，密度為1.8g/cm3。為研究機制數(shù)變化引起的數(shù)值噪聲的大小，以z=500m為分界面，界面以上采用L=1或3、界面以下采用L=3或5黏聲方程進行波場模擬。震源采用Ricker子波，震源位于(500m，100m)，主頻為30Hz，時間步長為1ms，在z=400m的水平線上接收，間隔為5m；采用時間二階、空間八階差分格式及PML邊界條件。

圖6～圖8分別為Q=20、60、100時均勻模型的變機制數(shù)正演結(jié)果。由圖6～圖8可以看出，變機制數(shù)模擬確實會引起數(shù)值噪聲，但數(shù)值噪聲非常小，只有把增益設(shè)置很大的時候才能看到數(shù)值反射。同時可以看出，Q值越大，機制數(shù)變化引起的數(shù)值干擾越小，Q值不變時，機制數(shù)由1變?yōu)?、由1變?yōu)? 和由3變?yōu)?引入的數(shù)值干擾依次減小。

圖6 Q=20時變機制數(shù)正演記錄

圖7 Q=60時變機制數(shù)正演記錄

圖8 Q=100時變機制數(shù)正演記錄

為研究這些數(shù)值干擾對波場模擬的影響，本文將差分格式不同引起的數(shù)值干擾的大小與地下速度變化產(chǎn)生的反射波大小對應(yīng)起來。首先提取這些數(shù)值干擾的振幅，然后保持模型其他參數(shù)不變，將z=500m以下的速度適當(dāng)增大Δv后進行正演；提取正演記錄中的反射波振幅與數(shù)值反射的振幅作對比，然后根據(jù)對比結(jié)果再調(diào)整速度增量(Δv)大小進行正演，直至反射波振幅和數(shù)值反射振幅相等。此時的速度增量與原速度的比值(Δv/v)即和數(shù)值噪聲的影響相當(dāng)，速度變化的大小反映數(shù)值干擾對波場模擬的影響大小。表2為與數(shù)值干擾量級相當(dāng)?shù)姆瓷洳ㄋ鶎?yīng)的地層速度變化?？梢钥闯?，機制數(shù)變化相同時，Q值越小，數(shù)值干擾對波場模擬的影響越大，Q值相同時，機制數(shù)從1變到5引起的數(shù)值干擾最大，從1變到3次之，從3變到5最小。但不管是哪種情況，所引入的數(shù)值干擾在數(shù)量級上都比地下速度發(fā)生1%的擾動所產(chǎn)生的反射還要微弱，因此可以忽略不計，說明黏聲波動方程變機制數(shù)有限差分正演方法具有很高的精度，可以較為精確地對地下復(fù)雜構(gòu)造、細(xì)微巖性差異甚至流體變化引起的地震波場異常進行模擬。

表2 與變機制數(shù)引起數(shù)值干擾相當(dāng)?shù)牡貙铀俣茸兓?%)

3 模型測試

由上文可知，在實際介質(zhì)模擬中，L=5模擬結(jié)果在絕大多數(shù)情況下能保持很高的精度。因此，以L=5的模擬結(jié)果作為標(biāo)準(zhǔn)進行對比分析。

3.1 簡單模型

為了驗證變機制數(shù)有限差分正演的效果，設(shè)計了一個大小為1000m×500m的簡單層狀模型，空間采樣間隔Δx=Δz=1m，速度、密度和品質(zhì)因子參數(shù)如表3所示。其中各層品質(zhì)因子根據(jù)巖性信息給出，例如Q=10模擬沙層，Q=40模擬含水沙層，Q=60模擬砂質(zhì)膠泥層，Q=105模擬一般巖層。

表3 層狀模型參數(shù)

控制模擬精度為95%，震源采用Ricker子波，震源位于(500m，10m)，主頻為30Hz，時間步長為0.1ms，檢波器位于模型頂界面，間隔為1m。采用時間二階、空間八階差分格式和PML邊界條件對該模型做波場數(shù)值模擬。圖9是L=1、3、5和自適應(yīng)變機制數(shù)的正演結(jié)果，四者的增益相同。從圖9可以看出，L=1模擬記錄在小炮檢距時能保持較高的精度，但當(dāng)炮檢距較大時其記錄的振幅偏大，原因是L=1時模擬品質(zhì)因子整體高于理論值很多，不能將介質(zhì)的吸收特性描述準(zhǔn)確；而L=3、5和變機制數(shù)的正演結(jié)果的初至、反射和多次波等主要特征差別很小，說明三者具有相似的精度。圖10為圖9中x=900m處的單道波形對比，可以看出L=1模擬結(jié)果的振幅誤差超過實際振幅值的一倍，并且兩個鄰近波峰沒能區(qū)分開來，說明L=1模擬不滿足精度要求；L=3的模擬精度相較于L=1提高了很多，但波形還是出現(xiàn)了一些小的畸變；變機制數(shù)與L=5的模擬結(jié)果波形吻合較好。為進一步研究變機制數(shù)模擬結(jié)果的精度，將圖9c與圖9d數(shù)據(jù)相減，殘差的最大振幅僅為圖9c數(shù)據(jù)最大振幅的0.0055%，因此可以認(rèn)為變機制數(shù)和L=5具有相同的模擬精度。圖11為不同機制數(shù)模擬需要的計算時間，變機制數(shù)模擬效率比L=5提高了27.4%。

圖9 層狀模型L不同時正演記錄對比

圖10 層狀模型L不同時正演x=900m處單道對比(a)及其局部放大顯示(b)

圖11 層狀模型不同L模擬計算時間對比

3.2 SEAM模型

為更接近實際和檢驗算法的穩(wěn)定性，建立一個SEAM模型，模型有751×876個網(wǎng)格點，網(wǎng)格間距Δx=Δz=5m。模型中存在連續(xù)和不連續(xù)的參數(shù)變化，如圖12所示，其中密度和品質(zhì)因子模型根據(jù)速度模型給出。震源采用主頻為30 Hz的Ricker子波，震源位于(2190m，20m)，時間步長為0.5ms，檢波器位于模型表面，間隔為5m，采用時間二階、空間八階差分格式，PML邊界條件。圖13為L=1、3、5和變機制數(shù)的正演結(jié)果，四者的增益相同。L=1時的模擬結(jié)果在大炮檢距時出現(xiàn)較大誤差；L=3、5和變機制數(shù)模擬結(jié)果的初至、反射、折射、散射和多次波等主要特征差別很小，說明三者具有相似的精度。圖14為圖13中x=3000m的單道波形對比，L=1時對介質(zhì)的吸收衰減特性描述不準(zhǔn)確，其振幅明顯偏大；L=3時模擬精度較高，但波形出現(xiàn)了一些小的畸變； 而變機制數(shù)與L=5的道記錄波形依舊吻合較好。

圖12 SEAM模型

圖13 SEAM模型L不同時正演記錄對比

圖14 SEAM模型L不同時正演x=3000m單道對比(a)及其局部放大顯示(b)

將圖13c與圖13d數(shù)據(jù)相減，殘差的最大振幅僅為圖13c數(shù)據(jù)最大振幅的0.0038%，因此可以認(rèn)為變機制和5機制具有相同的精度。圖15為不同機制數(shù)模擬的計算耗時，變機制模擬效率比L=5提高了23.25%。

由以上結(jié)果可知，變機制數(shù)有限差分正演方法對復(fù)雜模型也能得到較精確的模擬結(jié)果，不僅具有與機制數(shù)L=5相同的模擬精度，并且計算效率約提高了四分之一。

圖15 SEAM模型不同L模擬計算時間對比

4 結(jié)論

隨著業(yè)界對近地表研究的日益關(guān)注和對精度要求的日益提高，單機制GSLS模型有限差分正演已逐漸不能滿足高精度地球物理研究的需要。本文基于GSLS模型有限差分正演，提出了一種變機制黏聲波動方程有限差分正演方法，研究結(jié)果表明：

(1)實際介質(zhì)模擬中，模擬精度與機制數(shù)、Q值和波的傳播距離有關(guān)。本文給出了機制數(shù)的選取方法。當(dāng)控制在一定精度時，Q值和波的傳播距離又是已知的情況下，利用本文給出的方法可直接確定機制數(shù)。

(2)在Q值發(fā)生變化的邊界上，不同機制數(shù)黏聲方程差分格式的差異會引起人為的數(shù)值干擾。但這些數(shù)值干擾對波場模擬的影響極小，說明本文提出的黏聲波動方程變機制數(shù)有限差分正演方法具有很高的精度，可以對地下復(fù)雜介質(zhì)的細(xì)微結(jié)構(gòu)引起的微弱波場異常進行精細(xì)模擬。

(3)小機制數(shù)GSLS模型有限差分在介質(zhì)吸收衰減較弱時能得到較高精度的模擬結(jié)果，而當(dāng)介質(zhì)吸收衰減較強時其正演結(jié)果會出現(xiàn)明顯的振幅誤差和波形畸變，無法滿足信號分析等地球物理研究的精度要求；大機制數(shù)GSLS模型有限差分能得到高精度的模擬結(jié)果，但模擬時所需的存儲空間大，計算效率低。變機制數(shù)GSLS模型有限差分針對不同模型自動選取機制數(shù)，因此對任意模型都能達到模擬精度與計算效率的統(tǒng)一。

尚需指出，本文僅研究了單機制數(shù)、三機制數(shù)和五機制數(shù)二維黏聲波動方程有限差分正演，將其推廣到三維黏彈介質(zhì)是今后的研究方向。