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

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

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

0 引言

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

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

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

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

1 黏彈模型及其波動方程

1.1 GSLS衰減模型

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

(1)

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

(2)

(3)

圖1 GSLS衰減模型示意圖

品質因子Q定義為[26]

(4)

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

(5)

式中

(6)

當GSLS模型所用機制數L已知時，參數τσ,l、τ的最優值可由最小二乘反演得到

(7)

(8)

其中

(9)

1.2 二維黏聲波動方程

有關黏彈波動方程的推導Robertsson等[6]已做過詳細描述。二維黏聲微分方程為

(10)

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

2 變機制數有限差分正演精度分析

2.1 數值解與解析解的對比

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

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

圖2 Q=100時不同機制數的模擬結果

圖3為不同機制數黏聲方程模擬結果隨傳播距離的變化。可以看出：當傳播距離較近(500m)時，不同機制數的數值模擬結果都能與解析解相吻合； 中傳播距離(2500m)時單機制數和三機制數的模擬結果出現微小的誤差； 遠傳播距離(4500m)時各機制數的模擬結果都出現誤差，但機制數越大誤差越小。為衡量不同機制數波場模擬的準確性，計算模擬結果的均方根誤差(RMSE)

(11)

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

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

2.2 機制數的選取

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

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

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

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

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

Qmin=30.937lnD-157.2L=1

(12)

Qmin=7.6039lnD-40.364L=3

(13)

Qmin=5.7302lnD-37.40L=5

(14)

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

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

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

2.3 變機制數引起的數值干擾分析

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

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

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

圖6 Q=20時變機制數正演記錄

圖7 Q=60時變機制數正演記錄

圖8 Q=100時變機制數正演記錄

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

表2 與變機制數引起數值干擾相當的地層速度變化(%)

3 模型測試

由上文可知，在實際介質模擬中，L=5模擬結果在絕大多數情況下能保持很高的精度。因此，以L=5的模擬結果作為標準進行對比分析。

3.1 簡單模型

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

表3 層狀模型參數

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

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

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

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

3.2 SEAM模型

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

圖12 SEAM模型

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

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

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

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

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

4 結論

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

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

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

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

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