由于水和瓦斯在煤层中的运移通道只能是裂隙、微裂隙，所以煤层中水和瓦斯的运移方向、含量与煤层裂隙的密度、走向、填充物类型和连通性等参数密切相关[1]。如果能从地震数据中提取裂隙信息以掌握煤层裂隙的发育情况，对于煤田安全开采，降低透水、瓦斯突出等灾害性事故的发生具有重要意义[2-3]。煤层中广泛存在定向的彼此平行的裂隙系统，因此，为了接近煤系地层的实际情况，将煤层看作TTI等效介质[4]。对于TTI煤层来说，由于裂隙的作用，使得地震波传播时存在各向异性的差异，同时随着裂隙密度和裂隙填充物的不同，这种差异也随之变化[5]。因此可以利用这种各向异性来预测煤层裂隙的发育。
目前对于各向异性介质的研究多采用交错网格的方法，董良国提出了一阶弹性波动方程交错网格的高阶有限差分解法，同时给出了相应的稳定性条件[6];裴正林采用高阶交错网格有限差分算法对三维各向异性介质进行数值模拟，并推导出了三维各向异性介质PML法吸收层边界条件公式和相应的交错网格差分格式[7];李东会、董守华利用交错网格有限差分法对煤储层双相各向异性介质进行数值模拟，并将其应用于煤层气富集区的研究[8]。但是在交错网格中每个变量的不同分量都是交错定义的，对于没有定义的点需要进行变量插值，这降低了模拟精度。因此笔者在前人的基础上，给出了TTI介质的二维三分量应力-速度弹性波方程的Lebedev网格有限差分解法，利用Lebedev网格对TTI介质煤层进行正演模拟，同时结合AVO正演模拟的方法，分析煤层裂隙参数对正演结果的影响。从而为裂隙型煤层的地震数据处理、各向异性参数反演及解释方法的研究提供可靠依据。
1 基本原理
1.1 TTI介质二维三分量应力-速度方程
在弹性波正演模拟中，采用一阶应力-速度弹性波方程，结合运动微分方程、本构方程和几何方程可得到三维各向异性介质的波动方程(无体力作用)，并令
可得一阶应力-速度二维三分量弹性波
方程组:
(1)
式中，v，σ分别代表速度和应力;Cij为TTI介质弹性张量系数矩阵的元素。
1.2 TTI介质的Lebedev网格有限差分格式
在使用常规的交错网格进行数值模拟时，由于需要复杂的插值计算对未定义在交错位置上的变量进行插值，极大地降低了运行效率，同时引入了数值误差，降低了运算精度，因此Lisitsa等基于Lebedev的思想提出了Lebedev网格，该网格的特点是将应力、速度分量定义在同一网格上，但是应力、速度两个分量在网格点上交错分布，这相当于把需要插值的分量直接在网格点上定义出来，可以准确计算每个变量在各个方向的数值，而无需进行插值[9-11]，由此得到Lebedev网格机制下TTI介质波动方程的有限差分格式为
(2)
其中，N为空间差分精度;Cn为Lebedev网格N阶有限差分系数。同理可以推出vy,vz,σzz,σyz,σxz,σxy的差分格式。
1.3 反射系数
当煤层为TTI介质时，煤层顶底板界面为各向异性介质分界面。对于各向同性介质分界面上的P波反射系数，可以通过解Zoeppritz方程获得。当介质分界面两侧的地层为各向异性介质时，Zoeppritz方程不再适用。
首先研究双层TTI介质的反射系数公式，如图1所示，在xoz界面中P波入射时产生4种散射波，在上界面产生反射P波和反射SV波的，下层产生透射P波和透射SV波的。在上述坐标系下，5种波的波函数可以表示为
Uk=(fk+gk)Akexp{iω(t-pkx-rkz)}
(3)
式中，fk,gk为偏振方向;Ak为平面波振幅;pk,rk分别为平面波的水平慢度和垂直慢度;k=1,2，3,4，5分别表示入射P波、反射P波、反射SV波、透射P波、透射SV波;ω为圆频率。
根据弹性介质理论，这5种波在弹性界面上应满足位移连续和应力连续的边界条件，同时结合TTI介质的Hooke定律[12]，可以得到应力与应变的关系
(4)
其中，
(5)
图1 双层TTI介质模型
Fig.1 Double layer TTI medium model
将波函数方程代入边界条件，可以得到反射系数表达式[2]为
(RPP RPSH RPSV TPP TPSH TPSV)T=M-1N
(6)
式中，R和T分别代表反射系数和透射系数，下标代表波的类型;M和N为包含已知参数，如弹性系数、偏振方向和慢度值的矩阵[2]。同理还可以求出SV波入射时的反射、透射系数。
同理，可以分别求得横波向下入射、纵波向上入射、横波向上入射时的反射系数以及透射系数。当介质模型为3层时，反射系数可以通过双层介质模型的反射系数求得。
2 TTI煤层正演模拟
2.1 单层裂隙介质的波场特征
笔者采用Lebedev网格对TTI介质煤层进行正演模拟，分析了裂隙密度、裂隙填充物、裂隙开度等因素对波场传播的影响。设定模型大小为300 m×300 m，空间采样间隔dx=dz=1 m，时间采样间隔为0.1 ms，采用爆炸震源激发，震源位于中心，主频为100 Hz，其他物性参数见表1。
图2为不同裂隙密度下含饱和流体裂隙的TTI煤层二维三分量数值模拟波场快照。从图2可以看出，在裂隙密度一定的情况下，不同方向上的波场能量有所差异。当裂隙密度为0时，TTI介质退化为各向同性介质，3个方向上只存在纵波和横波，并且各个方向上速度相同。随着裂隙密度的增大，各向异性程度越来越高，具体表现为快慢横波分离程度逐渐增大，横波波前面由圆形变为椭圆形，长短轴的比值增大，而纵波波前面的变化较小，因此横波的速度方位各向异性对含水裂隙较为敏感。同时可以看出在含水裂隙中，3个方向上的横波能量较强，在y方向上纵波能量很弱，波场不明显。
表1 TTI介质模型参数
Table 1 TTI media model parameters
图2 不同裂隙密度下含流体裂隙煤层介质弹性波场快照
Fig.2 Elastic wave field of fluid fractured seam under different fracture density
图3为不同裂隙密度下含干裂隙的TTI煤层二维三分量数值模拟波场快照。由图3可以看出，在干裂隙的情况下地震波场传播与裂隙含流体时存在较大差异，3个方向上的能量差异较小，在y方向上纵波波形较裂隙含水时明显，当裂隙密度为0时，TTI介质退化为各向同性介质，波场传播与含流体时相同，各个方向上速度相同，但是随着裂隙密度的增大，纵波波前面由圆形变为椭圆形，长短轴的比值增大，各向异性程度变高，而横波的分裂程度较含流体裂隙小，因此纵波的速度方位各向异性对干裂隙较为敏感。
图3 不同裂隙密度下干裂隙煤层介质弹性波场快照
Fig.3 Elastic wave field of dry fractured coal seam under different fracture density
造成以上现象的原因如下，通过对Thomsen的各向异性参数的研究发现，纵波各向异性系数ε(v)主要受弹性常数C11和C33的影响，而横波各向异性系数γ(v)主要受弹性常数C44和C66的影响[13]。根据沈金松等的研究[14-15]，当煤层裂隙为干裂隙时，C11和C33对裂隙密度的变化更为敏感，而横波的速度各向异性变化明显较小，因此在波场记录中显示为纵波的速度方位各向异性程度高，横波的分裂程度较小。当煤层裂隙为含流体裂隙时，弹性常数C11和C33基本不受裂隙密度的影响，而C66则较为敏感，因此在波场记录中显示为纵波变化较小，而横波分裂现象明显。这一结论在下文AVO分析中亦能体现，裂隙密度对干裂隙煤层顶板的纵波反射系数的影响较大，而对含水裂隙煤层顶板的纵波反射系数基本无影响。
2.2 3层裂隙模型弹性波传播规律
为研究实际地层中地震波在煤层中的传播规律，设计3层介质模型，模型物性参数见表2，上下两层围岩为砂岩，中间为TTI煤层介质厚度为10 m。模型大小为500 m×500 m，空间采样间隔dx=dz=1 m，时间采样间隔为0.1 ms，采用爆炸震源激发，震源位于中心，主频为100 Hz。3层裂隙介质弹性波场数值模拟结果如图4所示。图4(a)为干裂隙情况下，裂隙密度为0.2时不同时刻x分量上的波场快照，其中各种类型的波区分较为明显，图中1表示煤层顶板的反射纵波，2表示煤层顶板的转换横波，3表示煤层底板的透射纵波，4表示煤层底板的反射纵波，5表示煤层顶板反射横波，6表示煤层底板的透射横波，图4(b)同。图4(b)为含流体情况下，裂隙密度为0.2时不同时刻x分量上的波场快照，其中各种类型波的划分与干裂隙时相同。此外，可以看出在煤层中存在多次波。根据黄饶、陈小宏等的研究发现多次波会引起地震波振幅的变化，对煤层的AVO分析会产生较大的影响，包括波形特征、频谱特征和AVO响应规律等[16]。因此在实际的处理中，会通过衰减多次波以降低对地震波振幅的影响[17]。
表2 3层裂隙介质模型物性参数
Table 2 Physical parameters of fractured medium model with three layers
图4 干裂隙下和含流体裂隙下裂隙密度为0.2时不同时刻波长快照
Fig.4 Snapshot of the wavelength at different wavelengths under a fluid fracture and dry fracture with a crack density of 0.2
2.3 TTI裂隙介质AVO正演分析
为考察煤层裂隙密度、裂隙填充物以及煤层顶底板等因素对煤层反射系数的影响，笔者设计了不同的裂隙密度并选用了不同的裂隙填充物。模型如图5所示，采用上下两层围岩厚度为100 m、宽200 m，中间煤层厚度为10 m的3层介质模型，在计算煤层弹性常数时采用CHENG模型[18]，其参数设置见表3。地震观测系统采用端点放炮的二维观测系统，取道间距2 m，共51道接收，最小炮间距20 m，震源采用频率100 Hz的Ricker子波。
通过表3建立TTI煤层3层介质模型，煤层中所含裂隙的纵横比为0.001，裂隙密度x分别设定为0,0.05,0.10,0.15,0.20，所裂隙填充物分比为干裂隙和含饱和流体裂隙。煤层的顶底板岩性设为泥岩，具体参数见表3。其中煤层厚度为10 m，属于薄层，所获取的煤层顶底板记录以复合波的形式存在。由于煤层顶板反射波极性为负，煤层底板反射波极性为正，笔者分别拾取目标层的正、负极性的最大振幅来进行研究。其中通过拾取负相位最大振幅绘制了图6～7所示的AVO曲线，拾取正相位最大振幅绘制了图8所示的AVO曲线。
图5 3层介质模型
Fig.5 Three layer dielectric model
图6(a)反映了裂隙密度对干裂隙煤层反射系数的影响，通过更改煤层介质的裂隙密度，得到不同煤层裂隙的AVO曲线。由图6(a)可以看出，当煤层裂隙为干裂隙时，裂隙密度的变化对于煤层纵波反射系数存在较大的影响，随着裂隙密度的增大纵波反射系数曲线的斜率基本保持不变，而截距随之减小。在入射角相同的情况下，裂隙密度越大其反射系数值越小。这是由于随着裂隙密度的增大，煤层与围岩之间的物性差异也不断变大，因此AVO曲线截距的绝对值在图中表示出随裂隙密度增大而不断变大的趋势。
表3 煤层模型参数
Table 3 Parameters of coal seam model
图6 裂隙密度对干裂隙煤层和含水裂隙煤层反射系数的影响
Fig.6 Effect of fracture density on reflection coefficient of dry fractured seam and water cut seam roof
图7 裂隙密度较小和较大时裂隙填充物对反射系数的影响
Fig.7 Effect of fracture filling material on reflection coefficient when fracture density is small and larger
由图6(b)中可以看出，当煤层裂隙为含饱和水裂隙时，裂隙密度的变化对于煤层纵波反射系数影响相对较小，随着裂隙密度的增大反射系数曲线的斜率和截距基本保持不变，这一结论与论文2.1节的数值模拟结果保持一致。造成这种现象的主要原因是盐水与煤之间的弹性模量差异较小，在裂隙开度很小的情况下，对纵波反射系数而言可以忽略裂隙密度的影响。
图7分别反映了裂隙填充物对小裂隙密度和大裂隙密度煤层反射系数的影响。所用裂隙密度分别为0.05和0.20，裂隙填充物为干裂隙和含水裂隙。从图7(a)可以看出，在相同裂隙密度的情况下，含水裂隙与干裂隙的反射系数斜率基本相同，含水裂隙的截距较大，曲线位于干裂隙反射系数曲线的上方。而图7(b)反映了在较大裂隙密度的情况下裂隙填充物对反射系数曲线的影响，可以看出含水裂隙的反射系数曲线截距依然大于干裂隙曲线，但是相对于图7(a)而言两种介质下的曲线差距变大。对比图7(a),(b)可以明显看出随着裂隙密度的增大，裂隙填充物对于反射系数的影响也随之增大，在裂隙密度为0.05时两种介质中的反射系数曲线差异要远远小于裂隙密度为0.2的情况。造成这一现象的原因是由于纵波各向异性系数ε(v)对干裂隙和含水裂隙的敏感程度不同，裂隙密度越大，2者之间的差异也越大。
图8 裂隙密度对干裂隙煤层反射系数的影响
Fig.8 Effect of fracture density on reflection coefficient of dry fractured seam
图8反映了干裂隙情况下煤层的反射系数曲线，与图6(a)对比可以看出正负相位提取的煤层反射曲线显著不同，正相位的煤层反射曲线其斜率为负值，随着入射角的增大其反射系数值不断减小，曲线的截距为正值，随着裂隙密度的增大其截距也随之增大，但是增大的幅度不断减小，这与负相位的反射曲线情况正好相反。
4 结 论
(1)煤层介质为干裂隙时，纵波随裂隙密度变化明显，而横波变化相对较小，纵波的速度各向异性对裂隙反应敏感，当煤层介质为含流体裂隙时，横波分裂现象明显，纵波变化较小，因此横波的速度各向异性对裂隙反应敏感。
(2)煤层顶板的反射系数为负值，曲线斜率为正值，随着入射角的增大，反射系数也随之增大。煤层底板的反射系数为正值，曲线斜率为负值，随入射角度的减小，反射系数减小，并且煤层反射系数的变化受多种因素的影响。
(3)煤层的裂隙密度、介质的各向异性程度、裂隙中的填充物性质都会对煤层反射系数产生较大影响。在煤层裂隙介质为干裂隙时，随着裂隙密度的增大，煤层反射系数曲线的斜率保持不变，截距的绝对值随之增加。在煤层裂隙介质为含饱和流体时，裂隙密度对于煤层反射系数基本无影响。对于裂隙填充物对煤层反射系数的影响受裂隙密度的影响，随着裂隙密度的增大，裂隙填充物对于反射系数的影响也相对增大。
参考文献(References):
[1] 朱培民,王家映,於文辉,等.用纵波AVO数据反演储层裂隙密度参数[J].石油物探,2001,40(2):1-12.
ZHU Peimin,WANG Jiaying,YU Wenhui,et al.Inverting reservoir fracture density using P-wave AVO data[J].Geophysical Prospecting for Petroleum,2001,40(2):1-12.
[2] 王宏伟,彭苏萍,杜文凤,等.HTI构造煤方位AVO正演[J].石油地球物理勘探,2014,49(6):1122-1130.
WANG Hongwei,PENG Suping,DU Wenfeng,et al.Azimuthal AVO in HTI tectonic coal[J].Journal of Geophysical Research Atmospheres,2014,49(6):1122-1130.
[3] 刘振峰,曲寿利,孙建国,等.地震裂缝预测技术研究进展[J].石油物探,2012,51(2):191-198.
LIU Zhenfeng,QU Shouli,SUN Jianguo,et al.Progress of seismic fracture characterization technology[J].Geophysical Prospecting for Petroleum,2012,51(2):191-198.
[4] 邓小娟,彭苏萍,林庆西,等.基于各向异性的薄煤层AVO正演方法[J].煤炭学报,2010,35(12):2053-2058.
DENG Xiaojuan,PENG Suping,LIN Qingxi,et al.AVO forward method for thin seam based on anisotropy[J].Journal of China Coal Society,2010,35(12):2053-2058.
[5] 陈同俊,王新,崔若飞.基于方位AVO正演的HTI构造煤裂隙可探测性分析[J].煤炭学报,2010,35(4):640-644.
CHEN Tongjun,WANG Xin,CUI Ruofei.The detectability analysis on HTI tectonic coal cracks by azimuthal AVO forward modeling[J].Journal of China Coal Society,2010,35(4):640-644.
[6] 董良国,马在田,曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究[J].地球物理学报,2000,43(6):411-419.
DONG Liangguo,MA Zaitian,CAO Jingzhong.A study on stability of the staggered-grid high-order diference method of first-order elastic wave equation[J].Chinese Journal of Geophysics,2000,43(6):411-419.
[7] 裴正林.三维各向异性介质中弹性波方程交错网格高阶有限差分法数值模拟[J].中国石油大学学报:自然科学版,2004,28(5):23-29.
PEI Zhenglin.Numerical simulation of high order finite difference method for elastic wave equation in three-dimensional anisotropic medium[J].Journal of the University of Petroleum：Natural Science，2004,28(5):23-29.
[8] 李东会,董守华,赵小翠,等.煤储层双相EDA介质的地震波场模拟[J].地球物理学进展,2011,26(2):654-663.
LI Donghui,DONG Shouhua,ZHAO Xiaocui,et al.Seismic wave simulation of biphase EDA medium in coal-bed media[J].Progress in Geophys,2011,26(2):654-663.
[9] LISITSA V,VISHNEVSKIY D.Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity [J].Geophysical Prospecting,2010,58(4):619-635.
[10] BERNTH H,CHAPMAN C.A comparison of the dispersion relations for anisotropic elastodynamic finite-difference grids[J].Geophysics,2011,76(3):WA43-WA50.
[11] RIMANO′CZY A,VATHY I.P-wave reflection coefficients for transversely isotropic media with vertical and horizontal axis of symmetry[J].Seg Technical Program Expanded Abstracts,1995,62(14):1566.
[12] 吴国忱.各向异性介质地震波传播与成像[M].东营:中国石油大学出版社,2006:43-15.
[13] THOMSEN L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954-1966.
[14] 沈金松,詹林森,马超.裂缝等效介质模型对裂缝结构和充填介质参数的适应性[J].吉林大学学报:地球科学版,2013,43(3):993-1003.
SHEN Jinsong,ZHAN Linsen,MA Chao.The applicability of fracture effective model to diffeerent parameters of fracture geometric strutures and media filled in fracture pores[J].Journal of Jinlin University:Earth Science Edition,2013,43(3):993-1003.
[15] 镇晶晶,刘洋,李敏.几种裂缝模型的波场传播特征比较[J].石油地球物理勘探,2012,47(6):908-917.
ZHEN Jingjing,LIU Yang,LI Min.Comparison of wave field characteristics among the fracture rock models[J].Oil Geophysical Prospecting,2012,47(6):908-917.
[16] 黄饶,陈小宏,李景叶.层间多次波模拟与影响分析[J].地球物理学进展,2009,24(3):981-987.
HUANG Rao,CHEN Xiaohong,LI Jingye.Simulation and effect analysis of interbed multiples[J].Progress in Geophys,2009,24(3):981-987.
[17] 张军华,缪彦舒,郑旭刚,等.预测反褶积去多次波几个理论问题探讨[J].物探化探计算技术,2009,31(1):6-10.
ZHANG Junhua,MIU Yanshu,ZHENG Xugang,et al.Several theoretical problems of deconvolution of prediction deconvolution[J].Computing Techniques for Geophysical and Geochemical Exploration,2009,31(1):6-10.
[18] CHENG C H.Crack models for a transversely anisotropic medium[J].Journal of Geophysical Research Atmospheres,1993,98(S1):675-684.
[19] 全红娟,朱光明,李桂花.地震波横波分裂现象的数值模拟特征[J].西北大学学报:自然科学版,2014,44(5):745-750.
QUAN Hongjuan,ZHU Guangming,LI Guihua.Numerical simulation characteristics of Swave splitting in seismic waves[J].Journal of Northwest University(Natural Science Edition),2014,44(5):745-750.