基于三相Biot速度应力波动方程的水合物储层波场数值模拟

石仁刚 魏周拓 葛新民 邓少贵 祁磊 姚志广

摘要：天然气水合物储层可以看作由骨架、水合物和孔隙水构成的三相孔隙介质。以三相Biot一阶速度应力波动方程作为声波在水合物储层传播的控制方程，用有限差分法（FDTD）进行数值模拟。波场快照可以清楚展现3个纵波和两个横波，与Biot理论相符。根据实测速度优化控制方程参数，使数值速度与实测速度相符。结果表明：孔隙度对水合物剪应力分量影响最显著，对水合物、孔隙水和骨架速度分量的影响次之，所以储层中孔隙度最适合用水合物剪应力分量来分析。骨架速度分量和应力分量受饱和度的影响较小，水合物剪应力分量受饱和度的影响最大，所以分析储层中水合物饱和度的最佳选择是水合物剪应力分量。

关键词：水合物储层可视化； 三相Biot波动方程； 有限差分法;孔隙度； 水合物饱和度； 完全匹配层

中图分类号：P 631.4   文献标志码：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

海域天然气水合物是目前最具潜力的能源之一。在其开发过程中，最重要的工作是确定水合物储层的位置和水合物饱和度，声波探测［1］是完成这项工作的有效技术手段之一。为了研究声波在储层介质中的传播规律，科学家们提出了各种声波模型，例如有效介质理论［2］、Biot-Gassmann 模型［3］、三相Biot理论［4］等等。由于水合物储层被视为由骨架、水合物和孔隙水3种介质构成的多孔介质［5］，因此用三相介质模型研究声波在水合物储层的速度与衰减规律，才能分析出三相介质的相关物理参数对声波的影响。Liu等［6］通过三相介质Biot二阶控制方程的特征方程，分析声波在三相介质中的速度与衰减，得出水合物储层的物理参数对声波速度与衰减的影响规律。Gei等［7］把三相介质Biot模型推广到四相介质，把孔隙液体中包含的气泡单独看做一相，考虑气液混合流体对声波的影响。上述模型能够在给定的条件下从理论上分析水合物储层中的声波传输特征。另外一种分析储层物理参数对声波影响的方法是数值模拟。波场响应数值模拟建立在水合物储层声波控制方程的数值解之上，目前这方面的科研工作多集中在两相介质（岩石骨架和孔隙液体）。Yan等［8］分析了两相介质中薄裂缝带对反射斯通利波的影响，发现反射斯通利波的振幅随裂缝的宽度和密度增大而增大，对裂缝的倾角并不敏感。Ou等［9］通过双相介质波动方程数值分析发现渗透率对斯通利波的影响规律。Carcione等［10］构建三相介质Biot波动方程，并通过波场快照展示了介质中的3种纵波和两种横波。但Caricone给出的三相介质Biot波动方程，对应的波速约为5 km/s，远高于水合物储层所测数据。因此需要对模型加以改进，以符合水合物储层的实测波速。波速的主控因素是体积模量、剪切模量和密度。在密度一定的情况下，要校正波速只有通过调节体积模量和剪切模量来实现。传统三相介质Biot模型，由固结系数反映孔隙固体对储层骨架体积模量和剪切模量的影响。但是通过Liu等［6］的数值研究发现，固结系数对波速的影响非常微弱，难以描述水合物储层中速度的复杂变化规律。Carione等［10］使用一种基于实验数据拟合得到的模量公式，把骨架模量表示成储层等效压力的函数，能够切实反映储层速度变化规律。数值模拟是分析参数对本构方程影响规律的有效技术手段之一［11-13］。笔者基于三相Biot一阶速度应力波动方程对水合物储层进行数值模拟，使用Carione模量公式表示储层骨架模量，用遗传算法对模量参数进行优化，使数值速度与实测速度误差最小。数值模拟采用的是时空交错网格上的有限差分（FDTD）算法，时间2阶精度，空间4阶精度，计算区域使用以微分方程为辅助函数的完全匹配层［14］ （ADEPML）进行截断。

1 三相介质Biot速度应力波动方程

根据三相Biot理论，水合物储层的一阶速度应力波动方程可以表示为

ρ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）为体应变随时间变化率，偏应变变化率如下：

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分别为骨架和水合物的应力张量;S为与孔隙水压力成比例的量，本文中在不引起混淆的地方统称应力。上标s、w和h分别表示关于骨架、孔隙水和水合物的相关变量。下标i，j（i，j=1，2，3，i≠j）表示变量在某个坐标方向上的分量，例如vs1表示骨架在x方向上的速度分量，下标*，i表示函数关于某个坐标求方向导数，例如vs1，1表示函数

vs1关于x的导数。

水合物储层的弹性常数表示为

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为三相介质的平均模量。

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为固结系数。Lee用固结参数α把骨架和水合物模量表示为

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为速度调节系数。

储层的等效密度表示为

ρ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表示孔隙空间拓扑结构的弯曲度，

α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为孔隙固体形状的参数，当孔隙固体是球形时rij=0.5。黏滞系数bij表示为

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

式中，η为孔隙流体黏滞系数；ks和kh分别为骨架和水合物的有效渗透率。

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

其中Sr=φw/φ是流体饱和度。

2 数值模拟

要实施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）的行列式。物理计算区域由300×300个单元构成，并且在外侧覆盖10个单元厚度的PML层。网格单元的长度为0.058 m，时间步长为11.254 μs，总迭代步数为 M=3000，正方形总计算区域的边长是1.736 m。总计算区域及PML层分布情况见图1，两个接收器与激励源在平行于x轴的同一条直线上，两个接收器分别离激励源为50个单元和100个单元。

各相速度分量和应力分量在网格上的位置见图2。其中Qzσs12、Qxσh11、Qzvs1、Qxvs2、Qxσs12、Qxσh12、Qzvs2、Qzvh2、Qxvs1和Qxvh1是相应分量的PML辅助函数［14］，例如Qzσs12整体是一个变量，表示σs12在z方向上PML层中的辅助函数，其他PML辅助函数可以类似表示。这些辅助函数只定义在PML层中，其他位置为零［15］。

2.1 差分格式

给出差分格式，其中时间方向为二阶精度，空间方向为四阶精度。为了便于展示，省略了差分格式中的空间坐标，各网格函数的空间坐标可以从图2对应位置处获得。

第一步：速度分量更新。

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为网格函数在第n时间层上的数值，表示相关变量的空间坐标;dt为离散时间步长;系数D31=（ρ21ρ32-ρ22ρ31）/det（ρ）从式（33）～（36）得到。Dxσh11为函数σh11沿x方向的四阶中心差分格式；Qxσs11为函数σs11在x方向的PML辅助函数。Qxσs11迭代公式如下：

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

式中，dx、kx和ax为PML层参数。其他参数定义和变量迭代公式可以类似得到。

第二步：应力分量更新。

θ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辅助函数Qxvh1用如下公式进行迭代更新：

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

其他PML辅助函数的迭代公式可以类似得到。

雷克子波源放置在中心单元各主应力分量上

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

其中τ0=1/f0，主频f0=1000 Hz，为了压制数值频散，本文中采用FCT［16］磨光技术。计算过程中三相Biot波动方程的参数见表1

2.2 波场可视化

二维三相Biot波动方程总共有13个分量，各分量的波场快照见图3～4，其中所有子图采用相同色标（色标见图4）。

根据三相Biot理论，三相Biot波动方程有5种波，3个纵波（每相有一个）和两个横波（两种固相各有一个）。从图3中可以看出有3个纵波，可以通过波速差异识别出来，在骨架、水合物和孔隙水中分别表现明显，其中第一纵波在储层骨架中比较明显，第二纵波在水合物中表现突出，第三纵波在孔隙水中表现明显。

两个横波在图4中剪应力分量表现明显，可以通过不同的波速分别出来，其中第一横波在骨架中表现明显，第二横波在水合物中表现突出。图3和图4清楚展现这5种波场，仿真结果符合Biot理论，从而验证了数值解的有效性。

在图3和图4中，看到从上到下发现波形越来越弱，振幅越来越小。也就是说，在骨架中波幅最大，水合物中波幅次之，孔隙水中波幅最小。从而也就解释了在孔隙介质声波测速岩石物理实验中，难以检测到孔隙水中的纵波的原因，波幅较小，很容易受到噪声干扰。

2.3 速度调节参数优化

对速度调节参数r，使用遗传算法对数值速度进行优化，使之与实测速度差最小。

优化方案如下：

（1）取计算区域上两点作为接收器，与波源在同一条平行于x轴的直线上，见图1;

（2）计算每个接收器上首波到达时间;

（3）根据两个接收器上首波到达时间差和两点间距，计算数值解速度vP，n（r）;

（4）构造目标函数

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

这里vP，d是实测数据［17］;

（5）用遗传算法求解速度调节参数r。

通过遗传算法求解，得出当r=1004.2时，数值速度达到vP，d=3215.2（m/s），与实际测速较为接近。在后面的计算中采用本节优化结果。

2.4 孔隙度分析

数值仿真计算可以用来分析水合物储层相关物理参数对声波响应的影响。用数值分析的方法研究孔隙度对声波响应的影响规律。

数值分析所采用的参数与前节相同。选取距离激励源为50个单元的接收器采集数据。各相的速度分量和应力分量随孔隙度的变化规律见图5～7。

图5为三相介质中速度分量随孔隙度变化情况。从图中可以看出水合物和孔隙水中的速度分量对孔隙度的变化更加敏感，而储层骨架中的速度分量对孔隙度的变化不太敏感。而且随孔隙度变大，三相介质中的速度变慢，但是变慢的幅度逐渐变小。

图6是三相介质中正应力分量随孔隙度变化情况。从图中可以看出储层骨架和水合物中的正应力分量对孔隙度的变化更加敏感，而孔隙水中的压力分量对孔隙度的变化不太敏感。

从图6中前两个子图可以看出，储层骨架中的敏感区间来得更早一些，相对于水合物中的情况来说。而且随孔隙度变大，三相介质中的正应力变化相速度变慢，但是变慢的幅度逐渐变小。

图7是储层骨架和水合物中剪应力分量随孔隙度变化情况。从图中可以看出两相中的剪应力分量对孔隙度的变化都很敏感，比较而言，水合物中敏感区间持续时间长些。

综合图5～7可以得出，水合物中剪应力分量对孔隙的变化最为敏感，最适合用来判断储层中孔隙度变化情况，其次是水合物储层中的速度分量。

2.5 饱和度分析

数值分析所采用的参数与前节相同。选取距离激励源为50个单元的接收器采集数据。各相的速度分量和应力分量随饱和度的变化规律见图8～10。

图8为三相介质中速度分量随水合物饱和度变化情况。从图中可以看出水合物和孔隙水的速度分量对饱和度的变化更加敏感，而骨架的速度分量对饱和度的变化不太敏感。而且随饱和度变大，三相介质中的速度快慢变化很小，不容易分辨。

图9为三相介质中正应力分量随水合物饱和度变化情况。从图中可以看出水合物和孔隙水的正应力分量对饱和度的变化更加敏感，而骨架的正应力分量对饱和度的变化不太敏感。而且随饱和度变大，三相介质中的速度快慢变化不大。从图9可以看出，水合物饱和度对水合物的正应力分量和孔隙水的压力影响最大。

图10为储层骨架和水合物中剪应力分量随水合物饱和度变化情况。从图中可以看出水合物的剪应力分量对饱和度的变化十分敏感，饱和度为27%时波形几乎是线性，而饱和度为47%时，波形变化十分剧烈。储层骨架中的剪应力分量对饱和度的变化，相对于储层骨架中的速度分量来说，敏感度大些，当饱和度变大时，相速度变大，振幅变大。

图8～10表明水合物剪应力分量对水合物饱的变化最为敏感，最适合用来定量分析水合物饱和度，其次是水合物储层中的速度分量。

3 结 论

（1） 数值仿真技术能够清楚展示三相Biot波动方程的3个纵波和两个横波，从而说明数值解的可靠性。

（2）水合物速度分量和剪应力分量对孔隙度的响应更加敏感。

（3）水合物速度分量和剪应力分量对水合物饱和度的响应更加敏感。水合物剪应力分量是用来定量分析储层孔隙度和水合物饱和度的最佳技术手段。

参考文献：

［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］ 范翔宇，么勃卫，张千贵，等.基于声波信号预测Kaiser应力点的水平地应力［J］.西南石油大学学报（自然科学版），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］ 郑多明，汪家洪，肖又军，等. 基于地震数值模拟的溶洞型储层地震特征分析［J］. 西南石油大学学报（自然科学版），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］.西南石油大学学报（自然科学版），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］ 杨书博，赵琪琪，王磊，等.基于柱坐标系复频移完全匹配层的随钻套管井声场有限差分模拟［J］.中国石油大学学报（自然科学版），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）; 中央高校基本科研业务费专项（20CX05011A）;国家自然科学基金项目（42174142，42374152，41404091）;山东省自然科学基金项目（ZR2023YQ034，ZR2020MD050）;中国石油科技创新基金项目（2021DQ02-0402）

第一作者及通信作者：石仁刚（1978-），男，讲师，博士，研究方向为波场模拟与成像。E-mail：shirg@upc.edu.cn。

引用格式：石仁刚，魏周拓，葛新民，等.基于三相Biot速度应力波动方程的水合物储层波场数值模拟［J］.中国石油大学学报（自然科学版），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.