子波穿过黏弹性节理的波形变化规律

胡世丽王星光王观石

（1.江西理工大学建筑与测绘工程学院，江西 赣州 341000；2.中国地质大学工程技术学院，北京 100193）

0 引 言

子波是应力波能量由震源通过复杂的地下路径传播到接收器所记录下来的质点运动远场时间域响应.子波在岩体的传播过程中，子波变化直接反映了波场的传播特征，是波场记录的基本单元，子波在地震资料的处理和解释中有着重要的作用，是波场正演和反演的前提，因而研究子波的波形变化规律具有重要的意义[1].

岩体存在大量不同尺度的结构面，导致岩体具有非均匀和非连续性等力学特性.针对这些形态各异和复杂多变的结构面，通常采用的方法是：在根据地质成因分类基础上，再根据岩体结构面的力学特性进行分类[2]，不同类型结构面，其力学特性相差非常大[3].应力波在岩体的传播过程中，不仅受到岩石本身的密度、孔隙率和各种微结构等内部因素的影响[4]，还受到温度、湿度、外荷载和地形等外部因素的影响，多因素共同作用影响应力波在介质中的传播特性的.Gaviglio[5]认为：在诸多的影响因素中，岩石中的结构面是最直接、最突出的影响因素之一，Mckenzie[6]的研究表明：应力波的衰减取决于裂隙的数量、宽度以及充填物的波阻抗，归纳起来，结构面对应力波传播主要有3个方面的影响：信号延迟、信号衰减和高频滤波[7-10].

目前，主要采用声时、波幅、频率3个参数进行超声检测和探测，已建立定量或半定量数学表达式，但波形还无法定量解释，其原因是声时、波幅和频率仅从某些方面揭示了应力波传播规律，而波形能够综合反映应力波传播规律[11-12].樊耀新[13]通过理论和试验研究表明子波波形变化比走时对断层的反应更为敏感，充分反映了波形变化规律在岩体力学参数测试和岩体结构探测方面的优势，因而开展应力波在岩体的传播过程中的波形变化规律研究具有广阔的应用前景.

本文将小尺度结构面抽象为黏弹性节理，采用Kelvin模型描述节理及其两侧岩体的力学特性，建立谐波在黏弹性节理的传播模型，分析Ricker子波穿过黏弹性节理后的波形变化规律.

1 谐波在黏弹性节理的传播模型

1.1 谐波在黏弹性节理的透射系数和反射系数

P波从岩体Ⅰ中以入射角θPI斜入射至黏弹性节理，进入岩体Ⅱ，P波在黏弹性节理处发生波形转换，采用Kelvin黏弹性模型描述岩体Ⅰ和岩体Ⅱ的力学特性，如图1所示.

图1 谐波在黏弹性节理的透射和反射

P波的位移位函数Φ为：

SV波的位移位函数Ψ为：

式（1）、（2）中下标m=1、2 时分别表示上行波和下行波，下标r、i分别表示复数的实部和虚部.B为P波的振幅，C为SV波的振幅，ω为角频率，t为时间，j为虚单位，dα为P波在z方向的波数，dβ为 SV波在z方向的波数.dα=kPcosθPI，dβ=kScosθSI.kP，kS分别为 P 波和SV 波的复波数，kP=ω/vP，kS=ω/vS，其中vP，vS分别为 P波的复波速和SV波的复波速.k为P波或SV波x方向波数，k=kPsinθPI.θPⅠ和 θPⅡ分别为 P 波的入射角和透射角，θSⅠ和θSⅡ分别SV波的反射角和透射角.

式（5a）、（5b）中 ν为介质的泊松比，为介质的复模量.

对于岩体Ⅰ和岩体Ⅱ均按Kelvin模型计算，可得复数形式的模量：

式（6）中为介质的复模量，E为介质的弹性模量，η0为介质的黏性系数.

在黏弹性节理处，满足广义Snell定律：

位移可根据质点的位移位与位移分量和应力分量关系，黏弹性节理处的位移分量和应力分量的计算公式为[14]：

式（8）至式（11）中u和w分别表示x和z方向的位移，σzz和τzx分别为节理处的法向应力和切向应力.

如图1所示，在岩体Ⅰ中，有下行P波、上行P波、SV波；在岩体Ⅱ中，有下行P波，下行SV波.由此得岩体Ⅰ中P波的位移位函数为：

岩体Ⅰ中SV波的位移位函数：

将式（12）和式（13）代入式（8）和式（9）可得岩体Ⅰ的位移，将式（12）和式（13）代入式（10）和式（11）可得岩体Ⅰ的应力.

岩体Ⅱ中P波的位移位函数为：

岩体Ⅱ中SV波的位移位函数为：

式（15）中kIi=kIIi，kIr=kIIr.式（12）至式（15）中下标Ⅰ和Ⅱ分别为岩体 I和Ⅱ.将式（14）和式（15）代入式（8）和式（9）可得岩体Ⅱ中的位移，将式（14）和式（15）代入式（10）和式（11）可得岩体Ⅱ中的应力.

黏弹性节理如图2所示，当节理厚度可以忽略不计，不考虑节理质量条件下，设谐波穿过黏弹性节理时，节理两侧应力连续和位移不连续，且位移不连续量等于节理两侧应力与刚度之比，即：

式（16a）至（16d）中Kx表示节理的切向等效刚度，Kz表示节理的法向等效刚度.

图2 黏弹性节理模型

将黏弹性节理的等效切向刚度和法向刚度分别定义为：

式（17）和式（18）中kn，kτ分别为黏弹性节理的法向刚度和切向刚度；ηn，ητ分别为节理法向黏性系数和切向的黏性系数.

将由式（8）至式（15）计算得到的位移和应力代入式（16），从而得到谐波斜入射黏弹性节理时透射和反射系数计算公式，将其整理为矩阵形式：

式（19）中A1、A2、A3、A4为 1×4 的矩阵；B为 4×1 的矩阵；C为4×1的矩阵.各矩阵的计算公式为：

式 （20a） 至 式 （20f） 中K=kIi-jkIr，PI=-dIαi1+jdIαr1，QI=-dIβi1+jdIβr1，PII=-dIIαi2+jdIIαr2，QII=-dIIβi2+jdIIβr2.

1.2 谐波在薄黏弹性节理的透、反射规律

设黏弹性节理两侧的岩体力学参数相同，密度ρ1=ρ2=2 500 kg/m3；泊松比 ν1=ν2=0.25；弹性模量E1=E2=30 GPa；黏性系数 η0=0.8 MPa·s，以下计算均采用这组岩体物理力学参数.入射谐波频率f=500 Hz.节理中各项参数取值为：法向黏性系数ηn=1.5 MPa·s；切向黏性系数ητ=1.3 MPa·s； 法向刚度kn=2.0 GPa/m；切向刚度kτ=2.0 GPa/m.

将上面的参数代入式（19）进行计算，得到P波斜入射黏弹性节理的透反射系数和入射角的关系如图3所示.从图3中可知，随着入射角的增加，反射P波的反射系数先减小后增加，在60°时有最小值；透射P波的透射系数呈现先增加后减小的趋势.而透射SV波的透射系数和反射SV波的反射系数随着入射角的增加，呈现先增加后减小的规律.

当入射角为30°时，由式（1）计算得到谐波的透射、反射系数与频率的关系，如图4所示.由图4可知，随着入射波频率的增加，透射P波的透射系数减小，而反射P波的反射系数逐渐增加.透射SV波的透射系数和反射SV波的反射系数随着入射波频率的增加，呈现逐渐增加的规律.可见，节理具有滤波特性，透射系数随着入射波频率的增大而减小，高频部分几乎不能通过节理.

图3 透、反射系数与入射角的关系

图4 透、反射系数与频率的关系

2 子波在薄黏弹性节理的波形变化

2.1 透射子波的计算

设入射波为振幅为1的Ricker子波，其时域表达式为：

式（21）中fM为峰值频率.

Ricker子波的频域表达式：

设由式（19）计算的透射系数为：

采用傅里叶逆变换得透射波的时域表达式：

式（24）中 ω=2πf.

将式（22）和式（23）代入式（24），计算得透射波的时域波形.设峰值频率fM=500 Hz的Ricker子波以30°入射至黏弹性节理，节理的力学参数与2.2节相同，由式（22）计算得到透射波的时域波形，如图5所示.由图5可知，Ricker子波穿过黏弹性节理后，透射波振幅减小，主瓣不突出，波形在时域上变得更平缓，即波形的分辨能力减弱.

图5 入射角为30°的入射波和透射波波形

为进一步描述Ricker子波穿过黏弹性节理后的波形变化，现采用波形相关系数和波形变化系数2个参数描述波形变化.

2.2 子波波形相关系数的变化规律

设入射波和透射波位移振幅分别为A0i和ATi，则两列波的波形相关系数为：

式（25）中和分别为 Ricker子波和透射波位移振幅的平均值.

当节理切向黏性系数ητ=1.3 MPa·s，法向刚度kn=2.0 GPa/m，切向刚度kτ=2.0 GPa/m时，将峰值频率fM=500 Hz的Ricker子波入射至黏弹性节理.改变法向黏性系数和入射角度，计算得到法向黏性系数和入射角度对波形相关系数的影响规律，如图6所示.图6表明，波形相关系数随着法向黏性系数增加近似成负指数减小，随着入射角的增加而增加.

当节理法向黏性系数ηn=1.5 MPa·s，法向刚度kn=2.0 GPa/m，切向刚度kτ=2.0 GPa/m时，将峰值频率fM=500 Hz的Ricker子波入射至黏弹性节理.改变切向黏性系数和入射角度，得到切向黏性系数和入射角度对波形相关系数的影响规律，如图7所示.图7表明，随着切向黏性系数的增加，波形相关系数基本不变，即切向黏性系数的改变对波形变化影响非常小.

图6 法向黏性系数与波形相关系数关系

图7 切向黏性系数与波形相关系数关系

当节理法向黏性系数ηn=1.5 MPa·s，切向黏性系数 ητ=1.3 MPa·s，切向刚度kτ=2.0 GPa/m，将 Ricker子波45°入射至黏弹性节理.改变节理法向刚度和入射波峰值频率，计算得到法向刚度和峰值频率对波形相关系数的影响规律，如图8所示.图8表明，随着峰值频率的增加，波形相关系数减小，波形相关系数随着法向刚度的增加而增加.

图8 法向刚度与波形相关系数的关系

当节理法向黏性系数ηn=1.5 MPa·s，切向黏性系数 ητ=1.3 MPa·s，法向刚度kτ=2.0 GPa/m，将 Ricker子波45°入射至黏弹性节理.改变节理切向刚度和入射波峰值频率，计算得到切向刚度和峰值频率对波形相关系数的影响规律，如图9所示.图9表明，随着峰值频率的增加，波形相关系数减小，随着切向刚度的增加，波形相关系数增加.

图9 切向刚度与波形相关系数的关系

2.3 子波波形变化系数的变化规律

为进一步研究透射波波形变化，现采用波形变化系数ξ来量化透射波形变化.定义波形变化系数：

图10采用和图6相同的节理参数，通过改变法向黏性系数及入射角，描述了它们对波形变化系数的影响规律.从图10中可知，随着法向黏性系数的增加，波形变化系数减小；随着入射角的增加，波形变化系数减小.

图10 法向黏性系数与波形变化系数的关系

图11中的节理参数同图7，通过改变切向黏性系数及入射角，描绘了它们对波形变化系数的影响规律.当入射角为0°时，随着切向黏性系数的增加，波形变化系数不变，即此时切向黏性系数的变化不改变透射波波形.当以其他角度入射时，随着切向黏性系数的增加，波形变化系数减小；随着入射角的增加，波形变化系数减小.

图11 切向黏性系数与波形变化系数的关系

图12采用和图8相同的节理参数，描述了改变法向刚度及峰值频率对波形变化系数的影响规律.随着法向刚度的增加，波形变化系数减小；随着峰值频率的增加，波形变化系数增加.

图12 法向刚度与波形变化系数的关系

图13中的节理参数同图9，描述了改变切向刚度和峰值频率对波形变化系数的影响规律.随着切向刚度的增加，波形变化系数减小；随着峰值频率的增加，波形变化系数增加.

图13 切向刚度与波形变化系数的关系

3 结 论

从波的位函数出发，基于应力连续和位移不连续假定，建立应力波在薄黏弹性节理中的传播模型.依据所建立的模型，获得了应力波在黏弹性节理中传播的透、反射系数解析解.将Ricker子波作为入射波，利用Ricker子波的频域表达式，运用傅里叶逆变换，得到透射波时域波形的计算表达式.采用Matlab计算得到透、反射系数及透射波波形，并用相关系数和波形变化系数2个参数描述子波穿过黏弹性节理的波形变化，讨论了子波峰值频率和入射角、黏弹性节理的法向和切向黏性系数、黏弹性节理的法向和切向刚度对子波波形变化的影响规律.

1）P波的反射系数，随入射角的增加，呈现先减小后增加的变化趋势；随入射波峰值频率的增加而增加.P波的透射系数，随入射角的增加先增加后减小；随入射波峰值频率的增加而减小.P波入射时形成反射SV波和透射SV波，反射系数和透射系数，随入射角的增加，都先增加后减小.

2）子波穿过黏弹性节理后，波形发生变化，波形相关系数和入射角及子波峰值频率有关，表现为随入射角的增加，波形相关系数增加；随峰值频率的增加，波形相关系数减小.同时波形相关系数还受到节理力学参数的影响，表现为：随节理黏性系数的增加，波形相关系数减小；随节理刚度的增加，波形相关系数增加；相对切向刚度和切向黏性系数，节理的法向刚度和法向黏性系数对子波波形变化影响更大，当子波垂直入射黏弹性节理分界面时，节理的切向黏性系数和切向刚度对波形相关系数的变化没有影响.

3）波形相关系数主要反映子波穿过黏弹性节理后发生的相位变化，而波形变化系数能综合反映子波的振幅和相位变化.波形变化系数随入射角的增加而逐渐减小；随入射波峰值频率的增加，波形变化系数增加.随节理黏性系数和刚度的增加，子波波形变化系数减小.

[1]Yang P J，Pan Y，Mu X，et al.Mixed-phase seismic wavelet extraction of SIMO systemby subspacemethod[J].Journal of China University of Petroleum，2010，34（1）：41-45.

[2]黄润秋.复杂岩体结构精细描述及其工程应用[M].北京：科学出版社，2002：4.

[3]赵莹.小波分析在松辽盆地北部高分辨率层序地层学中的应用[J].物探与化探，2013，37（2）：310-313.

[4]Shatilo A P，Sondergeld C S.Ultrasonic attenuation in Glenn Pool rocks，northeastern Oklahoma[J].Geophysics，1998，63： 465-478.

[5]Gaviglio P.Longitudinal waves propagation in a limestone：The relationship between velocities and density[J].Rock Mechanics and Rock Engineering，1989，22（4）：290-299.

[6]Mckenzie C K，Stacey G P，Gladwin M T.Sonic wave characteristics of rock mass[J].International Journal of Rock Mechanics and Mining Sciences，1982，19（1）：25-30.

[7]Cai J G，Zhao J.Effects of multiple parallel fractures on apparent attenuation of stress waves in rock masses[J].International Journal of Rock Mechanics and Mining Sciences，2000，37：661-682.

[8]Zhao J，Zhou Y X，Hefny A M，et al.Rock dynamic research related to cavern development for ammunition storage[J].Unnelling and Underground Space Technology，1999，14（4）：513-526.

[9]俞缙，钱七虎，林从谋，等.P波在改进的弹性非线性法向变形行为单节理处的传播特性研究[J].岩土工程学报，2009，31（8）：1156-1164.

[10]王卫华，李夕兵，左宇军.非线性法向变形节理对弹性P波传播的影响[J].岩石力学与工程学报，2006，25（6）：1218-1225.

[11]吴庆曾，李洪涛，杨进平.声波检测穿透信号波形变化规律和量化的探讨[J].物探与化探，2012，36（1）：144-148.

[12]谢忠球，肖宏彬，李珍玉，等.边坡软弱层评价中全波列声波测井的理论与实验研究[J].中南林业科技大学学报，2010，30（8）：108-111.

[13]樊耀新.断层对地震子波波形影响的研究[J].西北地震学报，1996，18（3）：54-59.

[14]王小杰，印兴耀，吴国忱.粘弹性介质地震波传播特征及反射特征研究[J].物探化探计算技术，2012，34（3）：258-266.