0 引言

煤层中裂隙通常比较发育，其存在必然导致地震波传播时出现各向异性现象.同时，由于煤层的厚度相对地震子波波长来说较小，而且裂缝形态复杂，在地震勘探中如何通过各向异性分析确定地下裂隙的发育情况是一个难度较大的问题.煤层中存在一种常见的裂隙为近似垂向裂隙，对此通常可以用HTI(Transverse isotropy with a horizontal axis of symmetry)介质来描述, 将存在垂向裂隙的煤层近似视为HTI型煤层.常用的裂隙各向异性模型主要包括Hudson模型、Thomsen模型和Schoenberg模型(刘财等, 2013, 2016).Hudson(1981)使用刚度矩阵建立Hudson模型，首次分析了裂缝参数对速度的影响，提出了介质含定向裂缝时有效弹性模量的计算方法；Thomsen(1986)及Tsvankin(1997)提出各向异性参数，分别明确横向各向异性介质中各向异性参数的物理意义；Schoenberg和Douma(1988)借助柔度矩阵，应用线性滑移(LSD)理论进行裂隙介质模拟；Bakulin等(2000)将Thomsen各向异性参数与法向弱度和切向弱度关联，分析了弱各向异性介质中饱和水和干裂隙中的弱度变化规律.然而，这些裂隙理论中的参数不能直接与裂隙相关联，引入与裂隙相关的参数进行数值模拟，对于分析不同类型裂隙引起的煤层地震波响应尤为重要.针对不同类型的裂隙介质，许多学者曾借助不同的数值模拟方法，得到地震多波合成记录(张文生和何樵登，1998；刘洋和李承楚，2000；轩义华等，2006；杜启振等，2009；刘财等, 2013, 2016).然而，地震波各向异性参数的处理是基于单一波场记录进行的，通过数值模拟得到的分量合成记录必须经过波场分离才能使用，波场分离必然产生一定误差，且由于误差的传递，在参数反演结果中产生更大误差，影响成像结果.因此，如何生成适于HTI煤层的单一波场记录是煤层各向异性分析的基础.

地震波各向异性响应是地震波各向异性研究的一个重要内容，尹志恒等(2012)和陈楠等(2017)建立物理模型或从理论方面分析HTI介质的方位各向异性.Li(1997)、Zhang和McMechan(2010)、裴正林等(2007)及Ekanem等(2009)通过P波或多波各向同性数据走时、衰减等参数进行了裂隙探测研究.这些研究成果相继被引入煤田地震勘探中进行裂隙含水或含气的识别，如：王赟等(2008)和杨德义等(2011)通过正演模拟分析弱度参数与煤层裂隙的联系，论证了可以通过法向弱度与切向弱度比来鉴别裂隙含瓦斯的情况；董守华(2008)和李东会(2012)通过测定不同方向的速度参数、弹性参数，计算煤样的各向异性参数；陈同俊等(2010)、彭苏萍等(2013)等从AVO方面分析HTI型煤层的振幅变化，指出其变化规律与裂隙密度及充填物有关.描述地震波各向异性响应的参数主要包括波场、速度、振幅等相关参数，而地震记录中难以直观体现影响地震波各向异性响应的相关参数，需加强参数与响应间的理论模拟与分析，为各向异性处理时选用相应地震数据，为提取参数信息提供依据，并进一步应用于裂隙煤层各向异性的研究.

本文针对HTI型煤层，推导不同裂隙类型的各向异性参数公式，实现了基于波场分解的地震波数值模拟，分析各向异性介质中地震波传播产生的各种波形的不同响应，并对HTI型煤层地震数据各向异性响应给出总结，为各向异性参数分析时如何选取适合的地震数据提出建议.

1 HTI介质各向异性参数及波场分解理论 1.1 HTI介质各向异性参数

当横向各向同性介质的对称轴与水平方向相重合时，称为HTI介质，用于近似地表示由于构造应力产生空间排列垂直裂隙群体而引起的横向各向同性.HTI介质中弹性矩阵为

(1)

其中, .

公式(1)一直被用于传统弹性参数的分析中，Tsvankin(1997)为了让各参数有更明确的物理意义，提出HTI介质各向异性参数，公式为

(2)

式中，ε、δ、γ为HTI介质各向异性参数，无量纲.Bakulin等(2000)引入裂隙弱度，用法向弱度Δ_N和切向弱度Δ_T描述各向异性参数与裂缝纵横比，裂缝密度相互关系为

(3)

其中，g≡v_s²/v_p²，k和μ为背景介质的体积模量和剪切模量，k'和μ'是裂隙中充填物的体积模量和剪切模量，c/a是裂隙纵横比，e为裂隙密度.

根据地震波弹性参数特征可知：干裂隙中k'=μ'=0，饱和水裂隙中μ'=0.本文为了对不同充填物的裂隙进行地震波波场响应分析，结合裂隙弱度与各向异性弹性矩阵，经推导得：

(4)

式中，上标为d、f的参数分别表示干裂隙、饱和水裂隙对应的参数.该公式适用于常温常压.

由公式(4)可以看出，当裂隙为干裂隙或饱和水裂隙时，各向异性参数主要与裂隙密度及纵横波速度比有关.

1.2 波场分解及模拟理论

波场分解不同于波场分离，波场分离是用散度或旋度等方法将纵、横波场分离，而波场分解是将原始波场中的纵、横波场进行分解.波场分离得到的波场与原始波场有区别，地震波各向异性分析对于数据精度的要求较高，采用波场分离得到的波场数据通常会因为误差的传递，影响各向异性数据分析的精度；而波场分解得到的波场仍为矢量波场，与原始输入矢量波场具有相同的振幅、相位等信息，不产生分离波场的误差，为后期做各向异性处理奠定基础.

Dellinger(1991)给出各向异性介质中的Christoffel方程为

(5)

式中，Γ为Christoffel矩阵，Γ_i，j=c_ijkln_jn_l，c_ijkl为弹性张量，n_j、n_l为归一化后的j、l方向的波数(i，j，k，l=1, 2, 3)，V为相速度,

将Christoffel方程展开，得到：

(6)

求解Christoffel方程得到A^P和A^SV为(A^P和A^SV分别表示P波和S波的极化向量)：

(7)

对于各向异性介质，通过式(8)在波数域实现实现P波和S波分解(Zhang and McMechan, 2010)，公式为

(8)

式中，为质点速度分量在波数域中的表示，对、作傅立叶反变换得到分解后的P波、SV波场.

各向异性介质中, 波场分解针对P波和S波的极化向量，利用其相互垂直的特性来实现纵横波场分离.波场分解法在波数域实现波场分解，当各向异性参数在空间上出现变化时，极化向量可能会发生变化.本文在数值模拟时，将真空法(Zeng et al., 2012)引入交错网格高阶有限差分中，将在交错网格单元网格的节点上分别定义应力和弹性常量，不需要采用弹性常量平均法来求取应力，只需通过密度平均法求取密度的平均值来获取速度.这样能尽可能保留模型既定参数，减少数值模拟中可能产生的振幅畸变，以适应煤层顶底界面物性参数突变时的数值模拟.

2 数值模拟与分析 2.1 不同充填物煤层地质模型

根据煤层特性，本文建立模型如图 1所示，煤层的顶、底界面为各向同性介质，煤层为常温常压HTI型介质(模型参数见表 1)，将交错网格高阶有限差分法与真空法相结合，采用400×800网格，网格大小Δx=Δz=2 m，震源位于网格点(400, 25)，时间步长为0.2 ms，纵波震源，Ricker子波主频70 Hz，参数满足收敛性条件，实现地震波数值模拟及波场分解，针对煤层中裂隙饱和水和干裂隙两种情况分别展开讨论.

以模拟结果的x分量及x分量的分解为例，图 2为煤层中裂隙密度为0.11且裂隙分别为干裂隙或饱和水裂隙时，通过文中所述的地震波数值模拟法得到的地震波记录.图 2a1、图 2b1分别为干裂隙和饱和水裂隙模型中地震波x分量的合成地震记录，图 2a2、图 2b2分别为干裂隙和饱和水裂隙模型中从x分量分解的P波记录，图 2a3、图 2b3分别为干裂隙和饱和水裂隙模型中从x分量分解的SV波记录.每一个合成记录(如图 2a1、图 2b1)都可以被完整分解到P波和SV波记录(如图 2a2和图 2a3，图 2b2和图 2b3)中.由于煤层较薄，煤层顶、底界面相互干涉，形成P波或SV波的复合反射波，来自煤层的反射P波或SV波均表现为单个界面的反射特征，调谐振幅较大.图中x分量的SV波反射能量较强；x分量的SV波干涉现象更明显，SV波各向异性特性随着偏移距的增加更加显著.这一特征与薄层复合反射系数变化规律相符.

由于地震单炮记录中能量以及走时等的变化差异较微弱，不便观察.从图 2中各记录抽取第150道，并按照各模型进行组合，如图 3所示.图 3a从左到右的各道依次为干裂隙模型的混合波场、P波波场、SV波波场的第150道记录；图 3b从左到右的各道依次为饱和水裂隙模型的混合波场、P波波场、SV波波场的第150道记录；图 3c为饱和水裂隙模型(图 3a)与干裂隙模型(图 3b)的差值，为了使差值对比更明显，该图作了增益调整.从图中可以看出，分解后的P波与混合波场中的P波部分形态、振幅、相位一致，SV波与混合波场中的SV波部分形态、振幅、相位一致，地震记录得到分解；SV波能量强于P波，且SV波差值明显大于P波差值，也就是说，SV波对于不同裂隙表现分异性能力更强.

结合图 2、图 3可知，由于裂隙的存在导致地震波各向异性的出现，使得各波形的同相轴在大偏移距上的双曲线形状有较明显影响，所以在进行地震波各向异性分析时，选取大偏移距数据是比较有利的；在双曲线对于各向异性的敏感程度上，SV波的非双曲线特性更显著.

因为对于不同裂隙密度的模型来说，直达波没有变化，为了使P波、SV波对比更明显，图 4将不同裂隙密度地震波混合波场记录中的直达波切除后，对比不同裂隙密度的P波、SV波波场单道记录.图 4a、图 4b分别为干裂隙模型和含水裂隙模型的第150道记录，各图中的第一道、第二道、第三道数据分别对应裂隙密度为0.01、0.11、0.33时的地震记录.从图中可以看出：干裂隙模型中，P波振幅随着裂隙密度的增大而明显减小，SV波变化不明显；饱和水裂隙模型中，随着裂隙密度的增加P波振幅减小，而SV波随着裂隙密度的增加，其振幅值先减小而后增大，也就是说当裂隙密度较大时，相对于P波来说，SV波具有较强的反射能量，且随着裂隙密度的增加，不同裂隙充填物类型的SV波差值进一步扩大.

薄层模型中波形记录相对复杂，结合波场快照，可以更好地识别各类波形，分析波场变化规律.图 5为煤层中裂隙密度为0.11且裂隙为干裂隙或饱和水裂隙时，地震波经过HTI型煤层(t=300 ms)时的波场快照，两图中的a1、b1分别为干裂隙和饱和水裂隙中地震波在x分量上的波场快照，图a2、b2分别为干裂隙和饱和水裂隙中x分量的P波波场快照，图a3、b3分别为干裂隙和饱和水裂隙中x分量的SV波波场快照，每一个分量(如图中的a1，b1)的波场快照都可以被完整地分解到P波和SV波波场快照(如图中的a2和a3, b2和b3)中.

图 5的波场快照与图 2的地震记录相对应，在介质分界面处，由于煤层较薄，上下界面反射耦合，形成复合反射波；裂隙充填物的差异在波场快照中得到表现，干裂隙介质在薄层中的波干涉现象比饱和水裂隙的反射波干涉现象明显，饱和水裂隙模型中分离出来的SV波能量较强.对比图 5中x分量的不同波形波场快照可知：P波或SV波经反射或透射后，沿x轴或z轴方向均发生相位反转，反射P波能量比透射P波能量弱，反射SV波能量强于透射SV波能量.在HTI型煤层中，相对于P波波前来说，SV波波前形状受各向异性变化更显著，椭圆扁率大，波前面能量发散.

根据Tsvankin(1997)对于HTI介质各向异性参数的定义可知，HTI介质可以看作是将VTI介质坐标轴旋转90°而得到，文中借用Thomsen(1986)推导的VTI介质的群、相速度公式，通过旋转生成HTI型煤层的群、相速度关系图，如图 6所示为表 1所示模型具有不同裂隙密度时地震波群、相速度随群、相角变化关系图，图中紫色、蓝色、红色曲线分别表示裂隙密度e为0.01、0.11、0.33时的速度-角度关系图.从图 6可以看出，在速度各向异性响应方面，当煤层裂隙密度较小(如e=0.01)时，无论是饱和水裂隙还是干裂隙，介质的各向异性表现不明显，基本接近各向同性，群、相速度关系图几乎为圆形；当裂隙密度逐渐增大到一定程度时，由于裂隙充填物不同，速度响应方面出现明显差异.同时，速度随着裂隙密度变化明显：裂隙为饱和水时，P波的群、相速度受裂隙密度影响相对显著，SV波的群、相速度受裂隙密度影响明显，相速度波前形状随着裂隙密度的增大由近似圆形变化为近似方形，群速度随着裂隙密度的增大出现倾向三叉现象，表现出明显的各向异性；而当裂隙为干裂隙时, P波的群、相速度明显受裂隙密度影响，波前形状变化剧烈，呈现强各向异性，虽然P波对于干裂隙变化敏感，但是由于其变化趋势为随着裂隙密度的增大呈现迅速减小趋势，也就是说，当裂隙密度较大时，会因为速度过小而不适合于做各向异性参数分析，SV波的群、相速度接近各向同性介质，波前形状为近似圆形.由此可见，在上述模型中，P波对干裂隙比较敏感，SV波对含水裂隙更敏感，但是无论是对于干裂隙或含水裂隙，当裂隙密度较大时，地震波速度会表现较强各向异性，此时，采用SV波数据进行各向异性速度分析仍是较优的数据选择.

2.2 不同围压煤层地质模型

当煤层处于压力状态下，可以通过实验室进行弹性参数测定(董守华, 2008; 李东会, 2012).对于煤层压力状态下的介质进行模拟时，本文仍采用如图 1所示的地质模型，其中：上、下层介质为各向同性层，纵、横波波速分别为2700 m·s^-1、1550 m·s^-1，密度为2.3 g·cm^-3；模型中间层为各向异性层，借用李东会(2012)对不同围压下(轴压为4 MPa)的煤样测试结果(如表 2所示)，并对弹性参数进行计算.通过本文的波场模拟及分解方法对含有HTI型煤层的地质模型进行正演模拟，得到不同围压下的地震记录图，如图 7所示.

从图 7可以看出, 由于煤层为低速薄层，对于反射波有明显的调谐效应, 地震波遇煤层产生反射时，煤层作为一个薄层的整体，煤层反射波为复合反射波，同时受煤层顶、底界面的综合影响；PP波反射振幅随着偏移距的增大而减小，且不同围压的PP波反射差异极微弱；大偏移距SV波上能分辨出不同围压而导致的地震记录的较小差异.作者在文献(Li et al., 2018)中研究了薄层介质复合反射系数，当煤层较薄时，其反射系数不再是单一介质的反射系数，而是可以用复合反射系数来表示.复合反射系数受单一界面反射影响，而且可能大于单一界面的反射系数，其反射振幅也可能大于顶、底单一界面的反射振幅.本文借用该文中的复合反射系数理论，并将其用于表 2所示的模型中，以地震记录上差异并不明显的PP反射波为例，得到薄层复合反射系数随着模型及观测参数变化的规律(如图 8所示)，分析不同围压下的煤层在地震反射振幅方面的响应.

图 8a为不同子波频率(煤层厚度7 m，围压4 MPa)时的煤层复合反射系数与入射角的关系，煤层复合反射系数随着地震子波的频率呈现非线性变化规律；图 8b为不同围压(子波频率70 Hz，煤层厚度7 m)时，煤层复合反射系数与入射角的关系，当煤层围压变化时，煤层复合反射系数出现明显差异；图 8c以围压为4 MPa的煤层为例，采用地震子波频率70 Hz，选取5个煤层厚度的参数进行试算，分析不同厚度对煤层PP波复合反射系数的影响.从图中可以看出，当煤层较厚，如17 m(接近λ接近)时，在小入射角范围内，复合反射系数趋于零，当入射角达到一定程度时，反射系数随着入射角的增加而增加.煤厚为9 m(≈λλ在)时，复合反射系数最大；当煤层较薄(如7 m)，复合反射系数小于调谐厚度时的复合反射系数，反射振幅也较小，在一定入射角度范围内，复合反射系数随着入射角的增大而减小，这一规律与图 7中PP波反射地震记录相符.综合图 8中的三组图可知，煤层复合反射系数受薄层的厚度影响最大，煤层厚度的调谐作用对反射波振幅有明显的影响；即使地震合成记录对于不同围压的响应差异不显著，在分解后的PP波复合反射系数图中，仍能看到明显差异，因此，不同围压导致的煤层各向异性的差异，在反射振幅方面的响应是可分辨的.

从以上模型试算结果可知，通过波场分解后的波场、速度、振幅响应特性，能明确分辨裂隙的不同状态，验证了通过地震波各向异性分析预测煤层裂隙的可行性.

3 结论

(1) 饱和水和干裂隙状态对于研究煤层裂隙、预测煤层裂隙型灾害有较大意义，而目前缺乏直接将裂隙与各向异性响应关联起来的公式.论文从HTI弹性矩阵及Tsvankin各向异性参数出发，结合Bakulin关于弱度的表示，推导出能表征干裂隙及饱和水裂隙模型的各向异性参数表达式.

(2) 由于目前针对裂隙进行各向异性反演时缺乏理想数据，论文实现了在波动方程有限差分数值模拟计算中从混合波场分解出独立的P波及S波，为各向异性参数分析、偏移成像或参数提取提供了数据.

(3) 论文通过建立不同类型模型，从地震波波场、速度、振幅等多方面，较系统地分析了不同裂隙煤层的地震波响应特征，论证对煤层裂隙进行各向异性分析是可行的，并为地震波各向异性分析提供数据选择建议.

(4) 本文为了明确分析薄煤层的各向异性响应特征，选取水平薄层模型进行试算，该算法同样适用于复杂模型，但是由于复杂模型中波场太复杂，不利于对变化规律做普遍性分析，对复杂模型如何总结其响应变化规律，有待进一步研究.