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

（東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌330013）

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

1 Biot模型地震波傳播數(shù)值模擬

1.1 交錯(cuò)網(wǎng)格差分方法

平衡方程（應(yīng)力部分）：

幾何關(guān)系（應(yīng)變部分）：

應(yīng)力與應(yīng)變的關(guān)系：

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

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

將方程組（4）和（5）聯(lián)合，便可得到質(zhì)點(diǎn)速度和應(yīng)力表示的一階彈性波波方程：

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

其中：

其余各式同理，其中為交錯(cuò)網(wǎng)格差分系數(shù)。

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

1.2 完全匹配邊界條件

在利用計(jì)算機(jī)進(jìn)行波動(dòng)方程數(shù)值模擬時(shí)，有個(gè)不可避免的問題，那就是邊界問題，因?yàn)樾枰诳臻g域選取一個(gè)有限范圍。學(xué)者們通常會在有限區(qū)域的外圍增加邊界吸收條件，其目的是盡可能利用有限計(jì)算區(qū)域，實(shí)現(xiàn)更為真實(shí)的無限空間的波動(dòng)特征，那么如何有效設(shè)置人為邊界對數(shù)值模擬的可靠程度及精確度至關(guān)重要。

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

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

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

ν——m維向量。

引入一個(gè)新變量，能夠?qū)⒋蠼鈫栴}表達(dá)式改寫成如下格式：

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

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

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

和：

在頻域內(nèi)，由（11）式可以得到：

由（12）式可以得到：

對自變量作復(fù)變換如下：

現(xiàn)在利用平面諧波解來研究該模型特性，設(shè)：

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

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

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

平面波動(dòng)系統(tǒng)（11）的解可以寫成如下形式：

并且滿足下列條件：

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

（2）右半空間內(nèi)，u因?yàn)槭艿阶枘嵊绊?，發(fā)生衰減；

（3）完全匹配吸收層內(nèi)，阻尼系數(shù)如下所示：

1.3 交錯(cuò)網(wǎng)格算法流程圖

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

圖1 交錯(cuò)網(wǎng)格差分算法流程圖[10]

1.4 Biot地震波傳播數(shù)值模擬

根據(jù)Biot模型建立相關(guān)地質(zhì)模型，數(shù)值模擬結(jié)果見圖2。

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

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

2.1 孔隙介質(zhì)相速度和衰減

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

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

圖2 孔隙介質(zhì)地震波場快照

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

式中：k——波數(shù)

α——構(gòu)造常數(shù)或曲度因子；

κ——滲透率；

η——孔隙液體粘滯系數(shù)；

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

ρf——孔隙液體密度。

并有：

式中：f——頻率；

c——波速；

ρ——沉積物總體密度；

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

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

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

μ——框架剪切模量。

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

式中：a——孔隙尺寸參數(shù)。

由上述方程可得：

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

由方程（23）可得：

其中：

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

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

2.2 波速和衰減數(shù)值分析

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

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

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

3 結(jié)論

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

圖3 不同粘滯系數(shù)下頻率與衰減和速度關(guān)系

圖4 不同孔隙度下頻率與速度和衰減關(guān)系

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

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

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

[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]吳從輝,李紅星,李勤武,等.基于波場分離的雙相各向同性介質(zhì)正演模擬[J].東華理工大學(xué)學(xué)報(bào)：自然科學(xué)版,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]門福錄.波在飽水孔隙彈性介質(zhì)中的傳播[J].地球物理學(xué)報(bào),1965:107-114.

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

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