基于Biot機制的孔隙介質地震波傳播特征研究

李倩宇，李紅星，黃 宇，章晨望，徐文斌，樊嘉偉，梅振繁

（東華理工大學地球物理與測控技術學院，江西南昌330013）

為了滿足日益增長的勘探需求，以及研究地震波傳播過程中的衰減及頻散等問題，國外早在20世紀中葉開展有關雙相孔隙介質的相關科研工作。Gassmann于1951年提出了Gassmann理論[1]，將孔隙度這個概念引入；1956年Biot建立了Biot理論，并首次預言慢縱波的存在[2-3]；Panneton等人對其進行改進，實現了Biot模型的三維波場正演模擬[4]；Stoll等人對Biot模型進行改進，建立了Biot-Stoll模型[5]；Dvorkin等人將Biot流動機制和噴射流動機制綜合考慮，即結合固流相相互作用的力學機制，提出了BISQ（Biot-squirt）模型[6]，其能夠較好反映孔隙介質中彈性波的傳播規律；Nicholas Chotiros在Biot和BISQ模型的基礎上建立了BICSQS模型[7]。與國外相比，國內對雙相孔隙介質的研究開展較晚，門福錄于1965年最早開始相關研究[8]；隨后唐應吾、劉銀斌、牟永光等人也相繼展開研究，這一系列孔隙介質理論為油氣等資源的勘探提供了極其重要的理論支持[9]。本文基于Biot機制對地震波的傳播進行數值模擬，取得較好效果，之后對其波場特征進行分析。并且基于Biot模型探討波速、衰減與頻率的相關性。

1 Biot模型地震波傳播數值模擬

1.1 交錯網格差分方法

平衡方程（應力部分）：

幾何關系（應變部分）：

應力與應變的關系：

二維下，彈性介質中質點在x、z軸方向上的速度分量如果用Vx、Vz表示的話，那么，在外力為零時的介質部分：

將式（2）和Vx、Vz代入式（3）并對時間t求導可得一階微分方程組：

將方程組（4）和（5）聯合，便可得到質點速度和應力表示的一階彈性波波方程：

在時間上，使用二階中心差分對式（6）進行差分離散；在空間上，采用任意高階差分：

其中：

其余各式同理，其中為交錯網格差分系數。

從交錯網格差分格式可看出，Vx、Vz在 (k+1/2)Δt時取得式(k-1)Δt時刻 σxx、σzz、σzx的值，而應力張量Vx、Vz則在 kΔt時取得式在 (k+1/2)Δt的值。在程序設計時，這點顯得尤為重要，理解交錯網格的意義不僅僅是空間上的交錯，同時還有時間上的交錯，這才是交錯網格最特別的地方。

1.2 完全匹配邊界條件

在利用計算機進行波動方程數值模擬時，有個不可避免的問題，那就是邊界問題，因為需要在空間域選取一個有限范圍。學者們通常會在有限區域的外圍增加邊界吸收條件，其目的是盡可能利用有限計算區域，實現更為真實的無限空間的波動特征，那么如何有效設置人為邊界對數值模擬的可靠程度及精確度至關重要。

在理論意義上，可以完全吸收任意角度、任意頻率入射的波動則稱為完全匹配層。由JP.Berenger在計算電磁場有限域波動問題中首先引入。迄今為止，比較系統的基于二階線性雙曲型偏微分彈性動力方程是仍然缺少的。在這給出Collino的二維各向同性介質PML吸收邊界條件與交錯網格差分格式，這很重要，它能夠在各向異性介質彈性波場模擬中運用，得到的效果比較好。

下列方程組為波動控制方程的一種形式，可以根據它得到不同形式的波動模型：

式中：A、B——m×m維矩陣；

ν——m維向量。

引入一個新變量，能夠將待求解問題表達式改寫成如下格式：

其中，‖為只保留與交界面方向平行的偏導數分量。如對y軸求偏導（而⊥則為只保留x軸的偏導分量）。

設下軸方向上的阻尼因子為d(x)，另外有個新變量，它對與交界面成法向的波動分量上施加阻尼，使左半空間的解答滿足系統在右半空間新波動的解答，這樣一個波動變量u同時被定義[10]。

于是，使u滿足下列控制方程是非常重要的：

和：

在頻域內，由（11）式可以得到：

由（12）式可以得到：

對自變量作復變換如下：

現在利用平面諧波解來研究該模型特性，設：

其中，iν0w-iAν0kx-iBν0ky=0 ，用相同方法，能夠找到式（12）的一組特解：

式（12）的解，u‖和u⊥需滿足條件如下：

取a‖+a⊥=ν0，即：

平面波動系統（11）的解可以寫成如下形式：

并且滿足下列條件：

（1）左半空間內，u≡v x≤0，可以保證在界面處無反射；

（2）右半空間內，u因為受到阻尼影響，發生衰減；

（3）完全匹配吸收層內，阻尼系數如下所示：

1.3 交錯網格算法流程圖

對有限差分法的技術路線（見圖1）和實現步驟進行研究后，能發現有限差分法的實現相對簡單并且模擬出的效果也比較好，因此本文在對其深入研究的基礎上，得到了令人滿意的效果。

圖1 交錯網格差分算法流程圖[10]

1.4 Biot地震波傳播數值模擬

根據Biot模型建立相關地質模型，數值模擬結果見圖2。

從圖2中的模擬結果可以清晰看出快縱波、慢縱波及橫波。對固流兩相的波場快照進行分析，可以得到兩者快縱波的波形相位具有一致性，與單相介質中的縱波類似。而固流兩相中慢縱波的波形呈現出相位相反的特征。

2 Biot模型波速和衰減特征分析

2.1 孔隙介質相速度和衰減

Biot機制基于孔隙介質的雙相特性，提出了流體的粘滯性對孔隙流體相對運動產生影響，并造成彈性波在孔隙介質內傳播過程中逐漸衰減。Biot理論模型利用變分原理得到的波動方程來自流體孔隙介質中波的傳播，之后對地震波在流體孔隙介質中的傳播進行理論預測，發現存在3種體波：快縱波、慢縱波及橫波，三者中快縱波頻散相對較小，而慢縱波在傳播過程中衰減與耗散嚴重。

Stoll用勢Φs、Φf、Ψs、Ψf定義固體骨架的位移向量u及液體的位移向量U：

圖2 孔隙介質地震波場快照

?為沉積物孔隙度。Biot關于標量勢的方程（用平面波解exp[i(k·x-ωt)]）：

式中：k——波數

α——構造常數或曲度因子；

κ——滲透率；

η——孔隙液體粘滯系數；

ρs——沉積物顆粒密度；

ρf——孔隙液體密度。

并有：

式中：f——頻率；

c——波速；

ρ——沉積物總體密度；

F——隨頻率增高的泊肅葉流偏移；

Kf——孔隙液體體積模量；

Ks——沉積物顆粒體積模量；

μ——框架剪切模量。

假定孔隙是圓筒狀，Biot可得F的表達式：

式中：a——孔隙尺寸參數。

由上述方程可得：

其中，O=MH-C2，P=2BC-EM-AH ，Q=

由方程（23）可得：

其中：

則縱波相速度和衰減可表示為：

在Kb=0，μ=0的特殊情況時（海底沉積物的這兩項值比 Ks小得多，大多情況下可以忽略），C=M=H ，形成EDF模型：

2.2 波速和衰減數值分析

通過對孔隙度、粘滯系數、孔隙氣體含量等參數依次進行改變，分析其縱波衰減和相速度隨頻率發生的變化，可以就上述參數對地震波衰減特征的影響進行探討，并且在理論上對地震波衰減特征進行分析。

起初對Biot模型中孔隙度、孔隙氣體含量、粘滯系數對速度的影響進行探討。數值分析中所用的頻率范圍為1～108Hz，孔隙度依次選擇20%、40%及60%，分別選取0.1Pa·s、0.01Pa·s和0.001Pa·s作為粘滯系數，氣體含量則取0、0.1及0.2，為了研究它們對地震波衰減的影響，因此對這些參數值進行改變。

從圖3～圖5中分析可知，Biot模型速度在低頻段上基本保持穩定，不過也存在著不足之處，那就是頻散相對來說較明顯；同時，Biot模型速度在頻率增大的同時也在增大，這能從速度與頻率的關系中分析得到。頻段相同的條件下，孔隙度、孔隙氣體和粘滯系數都與速度成反比。不同于速度的是，孔隙度、孔隙氣體含量和粘滯系數對衰減幾乎沒有影響，并且衰減與頻率成正比，103Hz作為頻率的一個分界點，衰減在不超過這個頻率的范圍內穩定增大，超過這個頻率后，衰減增長的速度變緩慢。

3 結論

（1）本文利用較易實現的有限差分法數值模擬Bi?ot模型地震波場，快縱波、慢縱波及橫波都可以從模擬結果中清晰看出。通過比較流相與固相的波場快照，可以看出兩者快縱波的波形相位一致，而慢縱波在固

圖3 不同粘滯系數下頻率與衰減和速度關系

圖4 不同孔隙度下頻率與速度和衰減關系

圖5 不同孔隙氣體含量下頻率與速度和衰減關系

相與流相中的波形則表現出相位相反的特征。在固相波場中，橫波波形最明顯，能量最強；在流相波場中，慢縱波比快縱波及橫波更清晰，這是不言而喻的。

（2）分析Biot模型速度與衰減的關系，顯而易見，Biot模型速度在低頻段上基本保持穩定，不過也存在著不足之處，那就是頻散相對來說較明顯；同時，Biot模型速度在頻率增大的同時也在增大，這能從速度與頻率的關系中分析得到。頻段相同的條件下，孔隙度、孔隙氣體和粘滯系數都與速度成反比。不同于速度的是，孔隙度、孔隙氣體含量和粘滯系數對衰減幾乎沒有影響，并且衰減與頻率成正比，103Hz作為頻率的一個分界點，衰減在不超過這個頻率的范圍內穩定增大，超過這個頻率后，衰減增長的速度變緩慢。

[1]Gassmann F. über die Elastizit?t Por?ser Medien[J].Vier?teljschr.naturforsch.ges.zürich,1951,96(96):1-23.

[2]Biot M A.Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid.I.Low-Frequency Range[J].Journal of the Acoustical Society of America,1956,28(2):168-178.

[3]Biot M A.Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid.II.Higher Frequency Range[J].Journal of the Acoustical Society of America,1956,28(2):179-191.

[4]吳從輝,李紅星,李勤武,等.基于波場分離的雙相各向同性介質正演模擬[J].東華理工大學學報：自然科學版,2017(1):58-65.

[5]Stoll R D,Bryan G M.Wave Attenuation in Saturated Sedi?ments[J].Journal of the Acoustical Society of America,1970,47:1440-1447.

[6]Dvorkin J,Nur A.Dynamic poroelasticity:a Unified Model with the Squirt and the Biot Mechanisms[J].Geophysics,1993,58(4):524-533.

[7]Chotiros N P.An Inversion for Biot Parameters in Water-satu?rated Sand[J].Journal of the Acoustical Society of America,2002,112(5):1853-1868.

[8]門福錄.波在飽水孔隙彈性介質中的傳播[J].地球物理學報,1965:107-114.

[9]李紅星.杭州灣沉積物原位聲學特性分析及淺表低速層研究[D].吉林大學,2008.

[10]梅振繁.雙相介質地震波場模擬與特征分析[D].東華理工大學,2014.