弹性波在岩石中的衰减包含扩散衰减、吸收衰减和散射衰减.通常所说的岩石衰减指吸收衰减和散射衰减，受岩石的裂纹、孔隙结构、流体的相互作用等(冯若等，1999；李振春和王清振，2007； 陈颙等，2009)影响.研究岩石的弹性波衰减特性，对于储层流体预测、裂缝识别、岩石内部结构预测等具有重要意义.

国内外学者针对多孔岩石的弹性波衰减机理与规律开展了大量的研究.在理论研究方面，主要围绕岩石内孔隙流体的影响展开，如Biot流理论(Biot，1956)、喷射流模型(Dvorkin et al., 1994)、BISQ模型(Dvorkin et al., 1995; 杨顶辉和张中杰，2000)、斑块饱和模型(Müller et al., 2008)等；魏修成等(2008)基于Biot理论研究了衰减随流体性质、渗透率、入射角的变化规律；巴晶等(2012)分析了非饱和岩石中的纵波频散与衰减.刘炯等(2010，2011)基于Biot理论和随机过程理论分析了随机孔隙介质中散射和流体流动引起的地震波衰减.张显文等(2010)分析了三维双相正交裂隙介质的衰减及频散方位特性.在实验研究方面，Klimentos and McCann(1990)，徐果明等(2002)分析了砂岩的声波衰减规律；史謌和邓继新(2005)分析了地层条件下泥页岩衰减的各向异性；赵群和郝守玲(2006)对煤样超声速度和衰减的各向异性进行了测试；刘向君等(2012)分析了含气饱和度对碳酸盐岩纵、横波衰减的影响.此外，国内外学者基于相似性原理构建地震物理模型，对简单缝、洞结构的弹性波衰减(Gómez et al., 2000；李琼等，2007；尹志恒等，2012)进行了实验研究.但是除孔隙流体的影响外，“干”岩石本身的非均质性、孔隙结构、骨架颗粒等引起的衰减还有待进一步的研究.对于孔隙结构复杂的储层，如碳酸盐岩，研究孔隙结构对干岩心声波衰减的影响对于提高储层流体预测的可靠性及裂缝识别的定量化具有重要作用.

通过真实岩心实验研究孔隙结构对干岩心声波衰减的影响，除了要测量超声波在天然岩心中的衰减，还要进行岩心CT分析，实验程序较繁琐，对设备、经费需求较高.与实验方法相比，数值模拟方法更易实现单因素控制，可定量考察各因素的影响，并能开展大量复杂孔隙结构的模拟实验，降低研究成本和周期，是对实验研究的重要辅助手段.国内学者对岩石的电阻率(Liu et al., 2009)、破坏(Ju et al., 2010)、渗流特性(Zhao et al., 2010)和裂缝介质声衰减(陈乔等，2012)分别开展了微观数值模拟研究，这表明岩石物理性质的微观模拟实际可行.本文采用数值模拟的方法对岩心尺度弹性波衰减进行了研究.首先根据弹性波动力学及相似性理论，构建二维多孔介质的弹性波传播数值模型，在此基础上，计算弹性波在不同孔隙结构“干”模拟岩心中的传播，分析孔隙结构对干岩心弹性波衰减的影响，为揭示复杂孔隙结构天然岩心的弹性波衰减规律提供一种有效途径.

对于沿z方向传播的平面波，不需要考虑扩散衰减，其振幅随传播距离的变化遵循指数衰减规律(陈颙等，2009)：

式中，α为衰减系数；z为波传播距离.根据式(1)，如果能记录下岩石入射端的激励信号和另一端的接收信号，则岩石的衰减系数α可表示为

式中，A₀为激励信号的幅值；A为接收信号的幅值；l 为岩心的长度.

天然岩心的孔隙形状复杂，孔隙尺寸平均在10~100 um量级.如果按照1 1的比例构建几何模型，网格尺寸将比孔隙尺寸更小，导致网格质量变差、计算量剧增、瞬态问题计算不收敛性等问题.此外，不规则孔隙形状使网格划分困难，数值计算出现奇点.本文在构建多孔介质二维几何模型时进行如下简化：

(1)采用椭圆孔代替不同形状的孔隙，通过椭圆的半长轴、纵横比控制孔隙的几何形状与尺寸；

(2)将真实岩心和孔隙按照适当比例进行放大.为了使数值模拟得到的结果与实际介质中波的动力学特性(振幅及衰减等)相似，需遵守几何相似原理和动力学相似性(牟永光，2003).

基于Auto CAD二次开发，随机生成每个椭圆孔的中心坐标和倾角，控制每个椭圆孔的半长轴和孔隙纵横比，构建几何模型.Auto CAD生成的几何模型可导入大多数有限元软件或网格剖分软件(如ANSYS，ABAQUS，Hypermesh等)中完成网格划分.

在弹性介质中，波动方程为

式中，K为介质的体积模量；μ为介质的剪切模量；ρ为介质的密度； u 为位移矢量.

由式(3)可知，为了使模拟介质与实际介质中的波动方程相同，需满足

式中，下标mod表示与数值模拟相关的量.当满足式(4)时，弹性波在模拟介质中的传播速度v_mod与在实际岩心中的传播速度v相同，根据几何相似原理应满足(牟永光，2003)：

式中，f，λ 为室内岩心超声波测试的频率和波长；l 为岩心的长度；d为孔隙直径.

对于同一块岩心，测量频率不同，测得的衰减系数也不同；一般随着频率升高，衰减系数增大(陈颙等，2009).声学理论证明，介质的吸收衰减系数(αa)与频率一次方成正比，即

在实验室超声脉冲测试条件下，一般波长远大于孔隙尺寸，其散射衰减满足长波长瑞利衰减，散射衰减系数(αs)可表示为

如果模拟岩心骨架与实际岩石骨架具有相同的材料属性，则由式(5)-(7)可得：

介质的衰减系数通常表示为吸收衰减系数和散射衰减系数之和，即：

由式(8)和(9)可得：

由式(1)，(5)和(10)得：

综上分析可知，当模拟岩心的骨架性质(弹性参数、密度)与实际岩心相同，且满足式(5)时，模拟得到的衰减系数与实际岩心的衰减系数之比等于几何相似比；波在模拟介质中传播的振幅衰减与在实际岩心中的振幅衰减(A/A₀)相同.

模拟岩心可以分为固体骨架和孔隙两部分(如图 1所示).骨架为均匀各向同性弹性材料，与实际岩心骨架具有相同的弹性模量和密度，该区域的运动方程为式(3).孔隙内饱和空气，由于空气的密度、弹性参数与固体骨架相比近似为零，在构建几何模型时将孔隙做挖空处理，内部孔隙边界设置为自由边界.将模拟岩心的一端作为发射端，施加激励信号；在发射端对侧边界中间部分等间距的设置信号接收点.

图 1

Fig. 1

图 1 数值模型构建示意图 Fig. 1 Numerical model of elastic wave propagation in porous media

在信号发射端y方向施加频率为f_mod的脉冲激励信号，产生纵波脉冲：

为了降低边界反射的影响，将信号发射端以外的三个外边界设置为低反射边界.在低反射边界上，根据邻近单元格的数据，施加一个合适的阻抗，使弹性波从模拟区域的边界传播出去而基本无反射(Cohen and Jennings，1983)：

式中，n，τ 为边界的单位法向量和切向量，σ 为应力张量，cp，cs为质点沿平行和垂直波传播方向的振动速度.

数值模拟的初始条件为

在几何模型基础上，利用有限元软件求解式(3)及(12)~(14)描述的定解问题，可以获得不同时刻的波场信息，并在不同接收点提取接收信号，根据式(2)可计算模拟岩心的衰减系数α_mod及对应实际岩心的衰减系数α.

由于碳酸盐岩孔隙结构复杂，研究孔隙结构对其干岩心声波衰减的影响具有代表性，本文数值模拟参数取值以碳酸盐岩为例.根据碳酸盐岩岩心数字图像统计的孔隙结构数据(Weger et al., 2009)，几何模型构建时，控制主要孔隙尺寸变化范围为100~1000 um，孔隙纵横比变化范围为0.05~1.根据相似性原理确定的模拟参数与比例因子如表 1所示.该比例因子条件下，数值模拟得到的衰减系数为实际岩心衰减系数的1/10.为了保证数值模拟的收敛和准确性，时间步长取脉冲信号周期的1/(25~50)，有限元网格的最大尺寸为脉冲信号波长的1/(5~8).

表 1(Table 1)

表 1 碳酸盐岩声衰减模拟的比例因子 Table 1 Scale factors used in numerical modeling of acoustic attenuation in carbonates 参数 室内岩心 超声波测试 比例因子 数值模拟 骨架属性 (95%方解石+5%粘土) K=71.9 GPa μ=29.6 GPa ρ=2.70 g/cm3 1 1 Kmod=71.9 GPa μmod=29.6 GPa ρmod=2.70 g/cm3 频率 f=0.5~2 MHz 10 1 fmod =0.05~0.2 MHz 波长 λ=2~10 mm 1 10 λmod=20~100 mm 岩心尺寸 l=2.5~8 cm 1 10 lmod=25~50 cm 孔隙尺寸 d=10~103 um 1 10 dmod=0.1~10 mm 孔隙纵横比 AR=0.01~1 1 1 ARmod=0.01~1 波长/孔隙尺寸 一般λ/d>10 1 1 λmod/dmod≈20 衰减系数(待求) α 10 1 αmod

表 1 碳酸盐岩声衰减模拟的比例因子 Table 1 Scale factors used in numerical modeling of acoustic attenuation in carbonates

对于模拟岩心中的每一个椭圆孔，其中心点坐标在模拟区域(125 mm×250 mm)内随机均匀分布，其长轴与波传播方向的夹角在0°~180°范围内随机均匀分布.在孔隙尺寸、模拟频率不变的条件下，分别在孔隙度为4%，10%，16%时重复进行10次随机模拟，结果如图 2所示.

图 2

Fig. 2

图 2 不同孔隙度时重复随机模拟的结果 Fig. 2 Results of repeated r and om modeling at different porosity

从图 2可以看出，当孔隙随机分布参数相同时，10次随机模拟的衰减系数结果基本一致，偏差小于5%.这表明，在孔隙随机均匀分布条件下，岩石衰减系数的数值模拟结果无论在单点还是在整体变化规律上都具有很好的可重复性.这为分析孔隙结构参数对岩石弹性波衰减的影响奠定了基础.

保持模拟频率、孔隙尺寸、孔隙数目不变(fmod=0.1 MHz， dmod=3 mm，n=800)，模拟孔隙规则均匀分布、随机均匀分布两种条件下弹性波在模拟岩心中的传播.图 3为孔隙纵横比AR=1时不同时刻的波场快照，图 4为孔隙纵横比AR=0.1时不同时刻的波场快照.

纵横比AR=1的孔为圆形孔，对应于铸模孔、溶孔等孔隙类型；而纵横比AR=0.1的孔隙呈扁平状，对应于粒间孔、吼道、张开的微裂隙等孔隙类型，在二维平面内具有各向异性.相同条件下，孔隙随机均匀分布(图 3d~f，图 4g~i)比规则均匀分布(图 3a~c，图 4a~f)引起的弹性波散射更强，形成复杂、不规则的散射，导致更多的散射能量损失.对于扁平孔，从图 4a~f可以看出，孔隙长轴与波传播方向夹角对弹性波传播的速度、衰减均会产生影响.对比90°夹角(图 4a~c)与0°夹角(图 4d~f)相同时刻对应的波场快照，90°夹角时弹性波传播速度更慢，质点振动速度的幅值更小，衰减更大；而孔隙夹角随机取值时(图 4g~i)，弹性波在模拟岩心中的传播速度、衰减介于0°和90°对应的速度、衰减之间.孔隙的空间分布即影响岩石的吸收衰减，也影响其散射衰减.

图 3

Fig. 3

图 3 AR=1时不同时刻质点Y方向速度的波场快照(a)~(c)孔隙规则均匀分布，(d)~(f)孔隙随机均匀分布. Fig. 3 Snapshot of vibration velocity in Y direction with AR=1(a)~(c)Pores are regular distributed，(d)~(f)Pores are uniformly r and om distributed.

图 4

Fig. 4

图 4 AR=0.1时不同时刻质点Y方向速度的波场快照(a)~(c)孔隙规则均匀分布，孔隙长轴与波传播方向夹角为90°；(d)~(f)孔隙规则均匀分布，孔隙长轴与波传播方向夹角为0°；(g)~(i)孔隙随机均匀分布，孔隙长轴与波传播方向夹角随机取值. Fig. 4 Snapshot of vibration velocity in Y direction with AR=0.1(a)~(c)Pores are regular distributed， and the angle between the pore major axis and the wave propagation direction is 90°;(d)~(f)Pores are regular distributed， and the angle is 0°;(g)~(i)Pores are uniformly r and om distributed， and the angle is r and omly valued.

在包含椭圆孔的平面区域S₀内，孔隙度可表示为

式中，φ为孔隙度；N为孔隙数目.类比裂隙密度的概念，记ε=Nd²/S₀为孔隙密度.已发表的实验结果(Klimentos and McCann， 1990; 徐果明等，2002)表明，随着孔隙度的增大，弹性波在岩石中传播时的衰减增大.从式(15)可知，孔隙数目增多、孔隙尺寸增大、孔隙纵横比增大都可使岩石的孔隙度增大.下面将根据控制变量法的原理，进一步分析孔隙随机均匀分布、不同孔隙度变化方式下，衰减系数的变化规律.

保持孔隙尺寸、孔隙数目不变，则孔隙密度不变，孔隙度随纵横比的增大线性增大.取孔隙尺寸d=300 um，在孔隙数目为300，800，1200时，分别模拟不同孔隙纵横比情况下弹性波的传播并计算衰减系数，结果如图 5所示.

图 5

Fig. 5

图 5 孔隙纵横比对衰减系数的影响 Fig. 5 Effects of aspect ratio on attenuation coefficient

图 5中，实心圆点为通过数值模拟计算的数据点；实线为根据对应的数值模拟数据点采用最小二乘法得到的拟合曲线；虚线为在孔隙尺寸不变条件下，等效无量纲柔量与孔隙纵横比的关系(Sevostianov and Kachanov， 2011; Vernik and Kachanov， 2012).从图 7可以看出，在孔隙尺寸不变的条件下，衰减系数随纵横比的增大呈指数规律增大，其关系可表示为

图 7

Fig. 7

图 7 传播相同距离后质点Y方向速度的波场快照，孔隙数目增多，散射场相近，衰减增大主要为吸收衰减(a)n=300，dmod=3 mm;(b)n=500，dmod=3 mm;(c)n=900，dmod=3mm Fig. 7 Snapshot of vibration velocity in Y direction after propagating the same distance， and cores with more pores show similar scattering but stronger absorbing and results in larger attenuation

在孔隙数目、孔隙尺寸保持不变时，衰减系数与等效无量纲柔量随纵横比的变化趋势基本吻合：当纵横比小于0.2时，衰减系数与等效无量纲柔量均对纵横比的变化不敏感；随着孔隙纵横比增大，特别是纵横比大于0.5后，岩石孔隙度增大，其等效柔量显著增大(或等效弹性模量减小)，衰减系数迅速增大.这种相关性表明孔隙纵横比影响岩石的等效柔量是影响弹性波在岩石中传播衰减的内在机制之一.对比图 3d~f和图 4g~i可以看出，在孔隙数目、孔隙尺寸相同时，不同纵横比孔隙造成的散射场基本相同，因此该条件下岩石的弹性波衰减随纵横比增大主要为吸收衰减增大.

孔隙密度受孔隙尺寸和孔隙数目的控制.保持孔隙纵横比(AR=0.55)不变，波长与孔隙尺寸之比近似不变，分别模拟不同孔隙尺寸(d=100~800 um)和孔隙数目(100~1500)时弹性波在模拟岩心中的传播并计算对应的衰减系数，结果如图 6所示.

图 6

Fig. 6

图 6 AR=0.55时衰减系数随孔隙密度的变化规律 Fig. 6 Relationship between attenuation coefficient and pore density at AR=0.55

如图 6所示，不论是通过改变孔隙尺寸还是孔隙数目来改变孔隙密度，在纵横比一定的条件下，岩石的衰减系数都随孔隙密度的增大而呈幂函数规律增大，即

当孔隙尺寸不变、孔隙数目增多时，岩石的总衰减系数增大，但散射场能量基本不变，衰减增大主要为吸收衰减增大(如图 7所示).当孔隙数目不变、孔隙尺寸增大时，岩石的总衰减系数也增大，这是由于岩石内岩石-空气界面增大，阻碍弹性波传播路径的范围增大，散射场能量增大，衰减增大主要为散射衰减增大(如图 8所示).

图 8

Fig. 8

图 8 传播相同距离后质点Y方向速度的波场快照，孔隙尺寸增大，散射场增强，衰减增大主要为散射衰减(a)n=300，dmod=1 mm;(b)n=300，dmod=2 mm ;(c)n=300，dmod=3 mm Fig. 8 Snapshot of vibration velocity in Y direction after propagating the same distance， and cores with larger pores show stronger scattering and results in larger attenuation