石仁剛 魏周拓 葛新民 鄧少貴 祁磊 姚志廣

摘要：天然氣水合物儲層可以看作由骨架、水合物和孔隙水構(gòu)成的三相孔隙介質(zhì)。以三相Biot一階速度應(yīng)力波動方程作為聲波在水合物儲層傳播的控制方程，用有限差分法（FDTD）進(jìn)行數(shù)值模擬。波場快照可以清楚展現(xiàn)3個縱波和兩個橫波，與Biot理論相符。根據(jù)實測速度優(yōu)化控制方程參數(shù)，使數(shù)值速度與實測速度相符。結(jié)果表明：孔隙度對水合物剪應(yīng)力分量影響最顯著，對水合物、孔隙水和骨架速度分量的影響次之，所以儲層中孔隙度最適合用水合物剪應(yīng)力分量來分析。骨架速度分量和應(yīng)力分量受飽和度的影響較小，水合物剪應(yīng)力分量受飽和度的影響最大，所以分析儲層中水合物飽和度的最佳選擇是水合物剪應(yīng)力分量。

關(guān)鍵詞：水合物儲層可視化； 三相Biot波動方程； 有限差分法;孔隙度； 水合物飽和度； 完全匹配層

中圖分類號：P 631.4?? 文獻(xiàn)標(biāo)志碼：A

文章編號：1673-5005（2024）03-0057-08?? doi：10.3969/j.issn.1673-5005.2024.03.006

Numerical simulation of hydrate reservoir based on three-phase Biot velocity-stress wave equaiton

SHI Rengang1，2， WEI Zhoutuo2，3， GE Xinmin2，3， DENG Shaogui2，3， QI Lei4， YAO Zhiguang4

（1.School of Science in China University of Petroleum（East China）， Qingdao 266580， China;2.State Key Laboratory of Deep Oil and Gas， China University of Petroleum（East China）， Qingdao 266580， China;3.School of Geosciences in China University of Petroleum（East China）， Qingdao 266580， China;4.CPOE Research Institute of Engineering Technology， Tianjin 300451， China）

Abstract：Natural gas hydrate reservoirs are considered as three-phase porous media， comprising skeleton， hydrate， and pore water. Utilizing the three-phase Biot first-order velocity stress wave equation as the governing equation for sound wave propagation in gas hydrate reservoirs， numerical simulations are conducted employing the finite difference time domain （FDTD） method. The wavefield snapshots vividly illustrate three compressional waves and two shear waves， consistent with Biots theory. Through parameter optimization of the control equation based on measured velocities， numerical velocities are aligned with measured velocities. The results highlight porosity as the most influential factor on the shear stress component of gas hydrates， followed by impacts on hydrates， pore water， and skeleton velocities. Consequently， porosity is best analyzed using the shear stress component of gas hydrates in reservoirs. The skeleton velocity and stress components are less affected by saturation， while the shear stress component of gas hydrates is most affected. Therefore， the optimal choice for analyzing gas hydrate saturation in reservoirs is the shear stress component of gas hydrates.

Keywords：hydrate reservoir visualization; three-phase Biot velocity-stress wave equation; finite difference time domain （FDTD）; porosity; hydrate saturation; perfectly matched layer

海域天然氣水合物是目前最具潛力的能源之一。在其開發(fā)過程中，最重要的工作是確定水合物儲層的位置和水合物飽和度，聲波探測［1］是完成這項工作的有效技術(shù)手段之一。為了研究聲波在儲層介質(zhì)中的傳播規(guī)律，科學(xué)家們提出了各種聲波模型，例如有效介質(zhì)理論［2］、Biot-Gassmann 模型［3］、三相Biot理論［4］等等。由于水合物儲層被視為由骨架、水合物和孔隙水3種介質(zhì)構(gòu)成的多孔介質(zhì)［5］，因此用三相介質(zhì)模型研究聲波在水合物儲層的速度與衰減規(guī)律，才能分析出三相介質(zhì)的相關(guān)物理參數(shù)對聲波的影響。Liu等［6］通過三相介質(zhì)Biot二階控制方程的特征方程，分析聲波在三相介質(zhì)中的速度與衰減，得出水合物儲層的物理參數(shù)對聲波速度與衰減的影響規(guī)律。Gei等［7］把三相介質(zhì)Biot模型推廣到四相介質(zhì)，把孔隙液體中包含的氣泡單獨看做一相，考慮氣液混合流體對聲波的影響。上述模型能夠在給定的條件下從理論上分析水合物儲層中的聲波傳輸特征。另外一種分析儲層物理參數(shù)對聲波影響的方法是數(shù)值模擬。波場響應(yīng)數(shù)值模擬建立在水合物儲層聲波控制方程的數(shù)值解之上，目前這方面的科研工作多集中在兩相介質(zhì)（巖石骨架和孔隙液體）。Yan等［8］分析了兩相介質(zhì)中薄裂縫帶對反射斯通利波的影響，發(fā)現(xiàn)反射斯通利波的振幅隨裂縫的寬度和密度增大而增大，對裂縫的傾角并不敏感。Ou等［9］通過雙相介質(zhì)波動方程數(shù)值分析發(fā)現(xiàn)滲透率對斯通利波的影響規(guī)律。Carcione等［10］構(gòu)建三相介質(zhì)Biot波動方程，并通過波場快照展示了介質(zhì)中的3種縱波和兩種橫波。但Caricone給出的三相介質(zhì)Biot波動方程，對應(yīng)的波速約為5 km/s，遠(yuǎn)高于水合物儲層所測數(shù)據(jù)。因此需要對模型加以改進(jìn)，以符合水合物儲層的實測波速。波速的主控因素是體積模量、剪切模量和密度。在密度一定的情況下，要校正波速只有通過調(diào)節(jié)體積模量和剪切模量來實現(xiàn)。傳統(tǒng)三相介質(zhì)Biot模型，由固結(jié)系數(shù)反映孔隙固體對儲層骨架體積模量和剪切模量的影響。但是通過Liu等［6］的數(shù)值研究發(fā)現(xiàn)，固結(jié)系數(shù)對波速的影響非常微弱，難以描述水合物儲層中速度的復(fù)雜變化規(guī)律。Carione等［10］使用一種基于實驗數(shù)據(jù)擬合得到的模量公式，把骨架模量表示成儲層等效壓力的函數(shù)，能夠切實反映儲層速度變化規(guī)律。數(shù)值模擬是分析參數(shù)對本構(gòu)方程影響規(guī)律的有效技術(shù)手段之一［11-13］。筆者基于三相Biot一階速度應(yīng)力波動方程對水合物儲層進(jìn)行數(shù)值模擬，使用Carione模量公式表示儲層骨架模量，用遺傳算法對模量參數(shù)進(jìn)行優(yōu)化，使數(shù)值速度與實測速度誤差最小。數(shù)值模擬采用的是時空交錯網(wǎng)格上的有限差分（FDTD）算法，時間2階精度，空間4階精度，計算區(qū)域使用以微分方程為輔助函數(shù)的完全匹配層［14］ （ADEPML）進(jìn)行截斷。

1 三相介質(zhì)Biot速度應(yīng)力波動方程

根據(jù)三相Biot理論，水合物儲層的一階速度應(yīng)力波動方程可以表示為

ρ11si+ρ12wi+ρ13hi+b12（vsi-vwi）+b13（vsi-vhi）=σsij，j，（1）

ρ21si+ρ22wi+ρ23hi+b12（vwi-vsi）+b23（vwi-vhi）=Si，（2）

ρ31si+ρ32wi+ρ33hi+b13（vhi-vsi）+b23（vhi-vwi）=σhij，j，（3）

sij=（K1θs+C12θw+C13θh）δij+2μ1dsij+μ13dhij， （4）

=-φwPw=C12θs+K2θw+C23θh，（5）

hij=（C13θs+C23θw+K3θh）δij+2μ3dhij+μ13dsij.（6）

其中θs=div（vs），θw=div（vw），θh=div（vh）為體應(yīng)變隨時間變化率，偏應(yīng)變變化率如下：

dsij=（vsi，j+vsj，i）2-div（vs）δij/3， （7）

dwij=（vwi，j+vwj，i）2-div（vw）δij/3， （8）

dhij=（vhi，j+vhj，i）2-div（vh）δij/3. （9）

式中，vs、vw和vh分別為骨架、孔隙水和水合物的速度向量；

σs和σh分別為骨架和水合物的應(yīng)力張量;S為與孔隙水壓力成比例的量，本文中在不引起混淆的地方統(tǒng)稱應(yīng)力。上標(biāo)s、w和h分別表示關(guān)于骨架、孔隙水和水合物的相關(guān)變量。下標(biāo)i，j（i，j=1，2，3，i≠j）表示變量在某個坐標(biāo)方向上的分量，例如vs1表示骨架在x方向上的速度分量，下標(biāo)*，i表示函數(shù)關(guān)于某個坐標(biāo)求方向?qū)?shù)，例如vs1，1表示函數(shù)

vs1關(guān)于x的導(dǎo)數(shù)。

水合物儲層的彈性常數(shù)表示為

K1=［（1-c1）φs］2Kav+Ksm，K2=φ2wKav， （10）

K3=［（1-c3）φh］2Kav+Khm，（11）

C12=（1-c1）φsφwKav， C13=（1-c1）（1-c3）φsφhKav，（12）

C23=（1-c3）φhφwKav， μ1=［（1-g1）φs］2μav+μsm， （13）

μ3=［（1-g3）φh］2μav+μhm，（14）

μ13=（1-g1）（1-g3）φsφhμav+μsh0（φsφh）2.（15）

式中，Ksm、μsm、Khm 、μhm分別為骨架和水合物的體積模量和剪切模量；Kav和μav為三相介質(zhì)的平均模量。

Kav=（1-c1）φsKs+φwKw+（1-c3）φhKh-1， μav=0.（16）

c1=KsmφsKs， g1=μsmφsKs， c3=KhmφhKh， g3=μhmφhμh .（17）

式中，c1、g1、c3、g3為固結(jié)系數(shù)。Lee用固結(jié)參數(shù)α把骨架和水合物模量表示為

Ksm=（1-φw-εφh）Ks1+α（φw+εφh）， μsm=（1-φw-εφh）μs1+αγ（φw+εφh） . （18）

Khm=φhKh1+α（1-φh）， μhm=φhμh1+αγ（1-φh） .（19）

其中

γ=1+2α1+α .

式中，φh和φw分別為水合物和孔隙水的體積占比，并滿足φ=φw+φh。本文中使用儲層的有效壓力來表示儲層骨架的模量

Ksm（pe）=r（1.73+0.23pe-1.7exp（-pe/0.17）），（20）

μsm（pe）=r（0.54+0.06pe-0.54exp（-pe/0.12））.（21）

式中，pe為儲層的等效壓力；r為速度調(diào)節(jié)系數(shù)。

儲層的等效密度表示為

ρ11=α13φsρs+（α12-1）φwρw+（α31-1）φhρh，（22）

ρ22=（α12+α23-1）φwρw，（23）

ρ33=（α13-1）φsρs+（α23-1）φwρw+α31φhρh.（24）

不同相之間的耦合密度表示為

ρ12=-（α12-1）φwρw，ρ23=-（α23-1）φwρw，（25）

ρ13=-（α13-1）φsρs-（α31-1）φhρh.（26）

這里αij表示孔隙空間拓?fù)浣Y(jié)構(gòu)的彎曲度，

α12=r12φs（φwρw+φhρh）φwρw（φw+φh）+1，（27）

α23=r23φh（φwρw+φsρs）φwρw（φw+φs）+1，（28）

α13=r13φh（φsρs+φhρh）φsρs（φh+φs）+1，（29）

α31=r31φs（φsρs+φhρh）φhρh（φh+φs）+1.（30）

式中，rij為孔隙固體形狀的參數(shù)，當(dāng)孔隙固體是球形時rij=0.5。黏滯系數(shù)bij表示為

b12=ηφ2wks， b23=ηφ2wkh，? b13=（1-ε）b013（φhφs）2.（31）

式中，η為孔隙流體黏滯系數(shù)；ks和kh分別為骨架和水合物的有效滲透率。

ks=ks0S3r， kh=kh0（φ/φh）2（φw/φs）3.（32）

其中Sr=φw/φ是流體飽和度。

2 數(shù)值模擬

要實施FDTD方法，需要把

s、w和h從方程（1）～（3）求解出來，這里給出水合物的第二個速度分量方程，其他兩相的速度分量可以用類似的方法獲得。

h2=［ρ11ρ22（-b13（vh2-vs2）-b23（vh2-vw2）+σh21，1+σh22，2）-ρ12ρ21（-b13（vh2-vs2）-b23（vh2-vw2）+σh21，1+σh22，2）-ρ11ρ32（-b12（vw2-vs2）-b23（vw2-vh2）+S2）+ρ12ρ31（-b12（vw2-vs2）-b23（vw2-vh2）+S2）+ρ21ρ32（-b12（vs2-vw2）-b13（vs2-vh2）+σs21，1+σs22，2）-ρ22ρ31（-b12（vs2-vw2）-b13（vs2-vh2）+σs21，1+σs22，2）］/det（ρ） .（33）

其中det（ρ）為密度矩陣（ρij）的行列式。物理計算區(qū)域由300×300個單元構(gòu)成，并且在外側(cè)覆蓋10個單元厚度的PML層。網(wǎng)格單元的長度為0.058 m，時間步長為11.254 μs，總迭代步數(shù)為 M=3000，正方形總計算區(qū)域的邊長是1.736 m?？傆嬎銋^(qū)域及PML層分布情況見圖1，兩個接收器與激勵源在平行于x軸的同一條直線上，兩個接收器分別離激勵源為50個單元和100個單元。

各相速度分量和應(yīng)力分量在網(wǎng)格上的位置見圖2。其中Qzσs12、Qxσh11、Qzvs1、Qxvs2、Qxσs12、Qxσh12、Qzvs2、Qzvh2、Qxvs1和Qxvh1是相應(yīng)分量的PML輔助函數(shù)［14］，例如Qzσs12整體是一個變量，表示σs12在z方向上PML層中的輔助函數(shù)，其他PML輔助函數(shù)可以類似表示。這些輔助函數(shù)只定義在PML層中，其他位置為零［15］。

2.1 差分格式

給出差分格式，其中時間方向為二階精度，空間方向為四階精度。為了便于展示，省略了差分格式中的空間坐標(biāo)，各網(wǎng)格函數(shù)的空間坐標(biāo)可以從圖2對應(yīng)位置處獲得。

第一步：速度分量更新。

Lτ1=-b12（vs2-vw2）|n-b13（vs2-vh2）|n+（Dxσs12/kx+Qxσs12+Dzσs22/kz+Qzσs22））|n+1/2，（34）

Lτ2=-b12（vw2-vs2）|n-b23（vw2-vh2）|n+（DzS/kz+QzS））|n+1/2，（35）

Lτ3=-b13（vh2-vs2）|n-b23（vh2-vw2）|n+（Dxσh12/kx+Qxσh12+Dzσh22/kz+Qzσh22））|n+1/2，（36）

vs，n+12=vs，n2+dt（D11Lτ1+D12Lτ2+D13Lτ3），（37）

vw，n+12=vw，n2+dt（D21Lτ1+D22Lτ2+D23Lτ3），（38）

vh，n+12=vh，n2+dt（D31Lτ1+D32Lτ2+D33Lτ3）.（39）

式中，|n為網(wǎng)格函數(shù)在第n時間層上的數(shù)值，表示相關(guān)變量的空間坐標(biāo);dt為離散時間步長;系數(shù)D31=（ρ21ρ32-ρ22ρ31）/det（ρ）從式（33）～（36）得到。Dxσh11為函數(shù)σh11沿x方向的四階中心差分格式；Qxσs11為函數(shù)σs11在x方向的PML輔助函數(shù)。Qxσs11迭代公式如下：

Qxσs，n+3211=-dtdxDxσs，n+1211k2x+（1-dt（dxkx+ax））Qxσs，n+1211.（40）

式中，dx、kx和ax為PML層參數(shù)。其他參數(shù)定義和變量迭代公式可以類似得到。

第二步：應(yīng)力分量更新。

θn+1s=Dxvs1kx+Qxvs1+Dzvs2kz+Qzvs2，

θn+1h=Dxvh1kx+Qxvh1+Dzvh2kz+Qzvh2，

θn+1w=Dxvw1kx+Qxvw1+Dzvw2kz+Qzvw2，

dn+111，s=Dxvs1kx+Qxvs1-θs/3，

dn+122，s=Dzvs2kz+Qzvs2-θs/3，

dn+111，h=Dxvh1kx+Qxvh1-θh/3，

dn+122，h=Dzvh2kz+Qzvh2-θh/3，

dn+112，h=（Dzvh1kz+Qzvh1+Dxvh2kx+Qxvh2）/2，

σs，n+3211=σs，n+1211+dt（K1θs+C12θw+C13θh+2μ1d11，s+μ13d11，h），

σs，n+3222=σs，n+1222+dt（K1θs+C12θw+C13θh+2μ1d22，s+μ13d22，h），

σs，n+3212=σs，n+1212+dt（2μ1d12，s+μ13d12，h），

σh，n+3211=σh，n+1211+dt（C13θs+C23θw+K3θh+2μ3d11，h+μ13d11，s），

σh，n+3222=σh，n+1222+dt（C13θs+C23θw+K3θh+2μ3d22，h+μ13d22，s），

σh，n+3212=σh，n+1212+dt（2μ3d12，h+μ13d12，s），

Sn+32=Sn+12+dt（C12θs+K2θw+C23θh）.

格式中的PML輔助函數(shù)Qxvh1用如下公式進(jìn)行迭代更新：

Qxvs，n+11=-dtdxDxvs，n1k2+（1-dt（dxkx+ax））Qxvs，n1.

其他PML輔助函數(shù)的迭代公式可以類似得到。

雷克子波源放置在中心單元各主應(yīng)力分量上

f（t）={1-2［πf0（t-τ0）］2}exp{-［πf0（t-τ0）］2}.

其中τ0=1/f0，主頻f0=1000 Hz，為了壓制數(shù)值頻散，本文中采用FCT［16］磨光技術(shù)。計算過程中三相Biot波動方程的參數(shù)見表1

2.2 波場可視化

二維三相Biot波動方程總共有13個分量，各分量的波場快照見圖3～4，其中所有子圖采用相同色標(biāo)（色標(biāo)見圖4）。

根據(jù)三相Biot理論，三相Biot波動方程有5種波，3個縱波（每相有一個）和兩個橫波（兩種固相各有一個）。從圖3中可以看出有3個縱波，可以通過波速差異識別出來，在骨架、水合物和孔隙水中分別表現(xiàn)明顯，其中第一縱波在儲層骨架中比較明顯，第二縱波在水合物中表現(xiàn)突出，第三縱波在孔隙水中表現(xiàn)明顯。

兩個橫波在圖4中剪應(yīng)力分量表現(xiàn)明顯，可以通過不同的波速分別出來，其中第一橫波在骨架中表現(xiàn)明顯，第二橫波在水合物中表現(xiàn)突出。圖3和圖4清楚展現(xiàn)這5種波場，仿真結(jié)果符合Biot理論，從而驗證了數(shù)值解的有效性。

在圖3和圖4中，看到從上到下發(fā)現(xiàn)波形越來越弱，振幅越來越小。也就是說，在骨架中波幅最大，水合物中波幅次之，孔隙水中波幅最小。從而也就解釋了在孔隙介質(zhì)聲波測速巖石物理實驗中，難以檢測到孔隙水中的縱波的原因，波幅較小，很容易受到噪聲干擾。

2.3 速度調(diào)節(jié)參數(shù)優(yōu)化

對速度調(diào)節(jié)參數(shù)r，使用遺傳算法對數(shù)值速度進(jìn)行優(yōu)化，使之與實測速度差最小。

優(yōu)化方案如下：

（1）取計算區(qū)域上兩點作為接收器，與波源在同一條平行于x軸的直線上，見圖1;

（2）計算每個接收器上首波到達(dá)時間;

（3）根據(jù)兩個接收器上首波到達(dá)時間差和兩點間距，計算數(shù)值解速度vP，n（r）;

（4）構(gòu)造目標(biāo)函數(shù)

J（r）=minvP，n（r）-vP，d.

這里vP，d是實測數(shù)據(jù)［17］;

（5）用遺傳算法求解速度調(diào)節(jié)參數(shù)r。

通過遺傳算法求解，得出當(dāng)r=1004.2時，數(shù)值速度達(dá)到vP，d=3215.2（m/s），與實際測速較為接近。在后面的計算中采用本節(jié)優(yōu)化結(jié)果。

2.4 孔隙度分析

數(shù)值仿真計算可以用來分析水合物儲層相關(guān)物理參數(shù)對聲波響應(yīng)的影響。用數(shù)值分析的方法研究孔隙度對聲波響應(yīng)的影響規(guī)律。

數(shù)值分析所采用的參數(shù)與前節(jié)相同。選取距離激勵源為50個單元的接收器采集數(shù)據(jù)。各相的速度分量和應(yīng)力分量隨孔隙度的變化規(guī)律見圖5～7。

圖5為三相介質(zhì)中速度分量隨孔隙度變化情況。從圖中可以看出水合物和孔隙水中的速度分量對孔隙度的變化更加敏感，而儲層骨架中的速度分量對孔隙度的變化不太敏感。而且隨孔隙度變大，三相介質(zhì)中的速度變慢，但是變慢的幅度逐漸變小。

圖6是三相介質(zhì)中正應(yīng)力分量隨孔隙度變化情況。從圖中可以看出儲層骨架和水合物中的正應(yīng)力分量對孔隙度的變化更加敏感，而孔隙水中的壓力分量對孔隙度的變化不太敏感。

從圖6中前兩個子圖可以看出，儲層骨架中的敏感區(qū)間來得更早一些，相對于水合物中的情況來說。而且隨孔隙度變大，三相介質(zhì)中的正應(yīng)力變化相速度變慢，但是變慢的幅度逐漸變小。

圖7是儲層骨架和水合物中剪應(yīng)力分量隨孔隙度變化情況。從圖中可以看出兩相中的剪應(yīng)力分量對孔隙度的變化都很敏感，比較而言，水合物中敏感區(qū)間持續(xù)時間長些。

綜合圖5～7可以得出，水合物中剪應(yīng)力分量對孔隙的變化最為敏感，最適合用來判斷儲層中孔隙度變化情況，其次是水合物儲層中的速度分量。

2.5 飽和度分析

數(shù)值分析所采用的參數(shù)與前節(jié)相同。選取距離激勵源為50個單元的接收器采集數(shù)據(jù)。各相的速度分量和應(yīng)力分量隨飽和度的變化規(guī)律見圖8～10。

圖8為三相介質(zhì)中速度分量隨水合物飽和度變化情況。從圖中可以看出水合物和孔隙水的速度分量對飽和度的變化更加敏感，而骨架的速度分量對飽和度的變化不太敏感。而且隨飽和度變大，三相介質(zhì)中的速度快慢變化很小，不容易分辨。

圖9為三相介質(zhì)中正應(yīng)力分量隨水合物飽和度變化情況。從圖中可以看出水合物和孔隙水的正應(yīng)力分量對飽和度的變化更加敏感，而骨架的正應(yīng)力分量對飽和度的變化不太敏感。而且隨飽和度變大，三相介質(zhì)中的速度快慢變化不大。從圖9可以看出，水合物飽和度對水合物的正應(yīng)力分量和孔隙水的壓力影響最大。

圖10為儲層骨架和水合物中剪應(yīng)力分量隨水合物飽和度變化情況。從圖中可以看出水合物的剪應(yīng)力分量對飽和度的變化十分敏感，飽和度為27%時波形幾乎是線性，而飽和度為47%時，波形變化十分劇烈。儲層骨架中的剪應(yīng)力分量對飽和度的變化，相對于儲層骨架中的速度分量來說，敏感度大些，當(dāng)飽和度變大時，相速度變大，振幅變大。

圖8～10表明水合物剪應(yīng)力分量對水合物飽的變化最為敏感，最適合用來定量分析水合物飽和度，其次是水合物儲層中的速度分量。

3 結(jié) 論

（1） 數(shù)值仿真技術(shù)能夠清楚展示三相Biot波動方程的3個縱波和兩個橫波，從而說明數(shù)值解的可靠性。

（2）水合物速度分量和剪應(yīng)力分量對孔隙度的響應(yīng)更加敏感。

（3）水合物速度分量和剪應(yīng)力分量對水合物飽和度的響應(yīng)更加敏感。水合物剪應(yīng)力分量是用來定量分析儲層孔隙度和水合物飽和度的最佳技術(shù)手段。

參考文獻(xiàn)：

［1］ WANG X H， XU Q， HE Y N， et al. The acoustic properties of sandy and clayey hydrate-bearing sediments［J］. Energies，2019，12（10）：1-11.

［2］ HELGERUD M B， DVORKIN J， NUR A， et al. Elastic-wave velocity in marine sediments with gas hydrates： effective medium modeling［J］. Geophysical Research Letters， 1999，26（13）：2021-2024.

［3］ LEE M W. Modified Biot-Gassmann theory for calculating elastic velocities for unconsolidated and consolidated sediments［J］. Marine Geophysical Researches， 2002，23（5）：403-412.

［4］ LEE M W， COLLETT T S. Elastic properties of gas hydrate-bearing sediments［J］. Geophysics， 2001，66（3）：763-771.

［5］ LECLAIRE P H， COHEN-TNOUDJI F，AGUIRRE-PUENTE J. Extension of Biots theory of wave propagation to frozen porous media［J］. The Journal of the Acoustical Society of America， 1994，96（6）：3753-3768.

［6］ LIU L， ZHANG X， WANG X. Wave propagation characteristics in gas hydrate-bearing sediments and estimation of hydrate saturation［J］. Energies， 2021，14（4）：804.

［7］ GEI D， CARCIONE J M， PICOTTI S. Seismic rock physics of gas-hydrate bearing sediments

［M］//MIENERT J， BERDT C， TREHO A M， et al. World Atlas of Submarine Gas Hydrates in Continental Margins， Cham： Springer，? 2022：55-63.

［8］ YAN S G， XIE F L， GONG D， et al. Borehole acoustic fields in porous formation with tilted thin fracture［J］. Chinese Journal of Geophysics， 2015，58（1）：307-317.

［9］ OU Weiming， WANG Zhuwen. Simulation of Stoneley wave reflection from porous formation in borehole using FDTD method［J］. Geophysical Journal International， 2019，217（3）：2081-2096.

［10］ CARCIONE J M， SERIANI G. Wave simulation in Frozen porous media［J］. Journal of Computational Physics， 2001，170（2）：676-695.

［11］ 范翔宇，么勃衛(wèi)，張千貴，等.基于聲波信號預(yù)測Kaiser應(yīng)力點的水平地應(yīng)力［J］.西南石油大學(xué)學(xué)報（自然科學(xué)版），2022，44（2）：40-48.

FAN Xiangyu， YAO Bowei， ZHANG Qiangui， et al. Study on the calculated method of horizontal in-situ stress based on the Kaiser stress point predicted by acoustic signal［J］. Journal of Southwest Petroleum University （Science & Technology Edition）， 2022，44（2）：40-48.

［12］ 鄭多明，汪家洪，肖又軍，等. 基于地震數(shù)值模擬的溶洞型儲層地震特征分析［J］. 西南石油大學(xué)學(xué)報（自然科學(xué)版），2023，45（6）：57-68.

ZHENG Duoming， WANG Jiahong， XIAO Youjun， et al. Seismic characteristics analysis of karst cavity reservoirs based on seismic numerical simulation［J］. Journal of Southwest Petroleum University （Science & Technology Edition）， 2023，45（6）：57-68.

［13］ 張濤，孫天禮，陳偉華，等.不同孔-縫-洞組合碳酸鹽巖儲層氣水兩相孔隙尺度流動模擬［J］.西南石油大學(xué)學(xué)報（自然科學(xué)版），2023，45（5）：88-96.

ZHANG Tao， SUN Tianli， CHEN Weihua， et al. Pore-scale simulation of gas-water two phase flow in carbonate reservoir with different combinations of pore， network and hole［J］. Journal of Southwest Petroleum University （Science & Technology Edition）， 2023，45（5）：88-96.

［14］ MARTIN R， KOMATITSCH D， GEDNEY S D， et al. A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations（ADE-PML） ［J］. Computer Modeling in Engineering and Sciences， 2010，56（1）：17-40.

［15］ 楊書博，趙琪琪，王磊，等.基于柱坐標(biāo)系復(fù)頻移完全匹配層的隨鉆套管井聲場有限差分模擬［J］.中國石油大學(xué)學(xué)報（自然科學(xué)版），2022，46（1）：62-71.

YANG Shubo， ZHAO Qiqi， WANG Lei， et al. Finite difference modeling of acoustic fields in cased holes under LWD conditions based on CFS-NPML in a cylindrical coordinate system［J］. Journal of China University of Petroleum（ Edition of Natural Science），2022，46（1）：62-71.

［16］ HUANG M X， XUE T X， WANG Y X. Seismic wave field simulation on the symplectic and FCT algorithm［J］. Progress in Geophysics， 2015，30（2）：565-570.

［17］ LIU J， ZHANG J， MA F，? et al. Estimation of seismic velocities and gas hydrate concentrations： a case study from the Shenhu area， northern South China Sea［J］. Marine and Petroleum Geology， 2017，88：225-234.

（編輯 修榮榮）

基金項目：CNPC重大專項（ZD2019-184-001）; 中央高?；究蒲袠I(yè)務(wù)費專項（20CX05011A）;國家自然科學(xué)基金項目（42174142，42374152，41404091）;山東省自然科學(xué)基金項目（ZR2023YQ034，ZR2020MD050）;中國石油科技創(chuàng)新基金項目（2021DQ02-0402）

第一作者及通信作者：石仁剛（1978-），男，講師，博士，研究方向為波場模擬與成像。E-mail：shirg@upc.edu.cn。

引用格式：石仁剛，魏周拓，葛新民，等.基于三相Biot速度應(yīng)力波動方程的水合物儲層波場數(shù)值模擬［J］.中國石油大學(xué)學(xué)報（自然科學(xué)版），2024，48（3）：57-64.

SHI Rengang， WEI Zhoutuo， GE Xinmin， et al. Numerical simulation of hydrate reservoir based on three-phase Biot velocity-stress wave equaiton ［J］. Journal of China University of Petroleum （ Edition of Natural Science ）， 2024，48（3）：57-64.