基于射线理论的各向异性煤层地震响应

李 勤， 李庆春， 张 林

（1.西安科技大学地质与环境学院，陕西西安 710054；2.长安大学地质工程与测绘学院，陕西西安 710054 3.陕西省地质调查中心，陕西西安 710016）

基于射线理论的各向异性煤层地震响应

李 勤1， 李庆春2， 张 林3

（1.西安科技大学地质与环境学院，陕西西安 710054；2.长安大学地质工程与测绘学院，陕西西安 710054 3.陕西省地质调查中心，陕西西安 710016）

以VTI型煤层为研究对象，从介质中的地震波的射线方向、走时计算等多个方面展开讨论，分析VTI型煤层的试射法射线追踪技术，充分利用VC＋＋的面向对象类框架泛型编程的特点和Fortran数值计算方面的优势，通过混合编程实现VTI型煤层中的地震波试射法射线追踪；通过各向同性煤层和VTI型煤层的射线追踪试算对比，论证对VTI型煤层进行各向异性射线追踪的必要性，同时验证本文VTI型煤层射线追踪方法的正确性.最后，通过计算时间成本比较，表明该方法计算速度较快，能满足地震正演模拟的要求.该方法适用于各向异性煤层的多波射线追踪.

VTI型煤层；射线追踪；各向异性

0 引言

煤层裂隙是国内外煤矿地震勘探研究的重点和难点，裂隙发育复杂［1－3］，其中一种较为典型的裂隙系统可近似为平行于地表的裂隙系统，这种裂隙系统称为VTI型，而这种煤层就称之为VTI型煤层［4］.VTI型煤层作为各向异性煤层的一种常见模型，裂隙发育导致的各向异性会对地震波的记录产生较大影响［5－6］.但是，由于煤层相对地震子波来说较薄，各向异性煤层地震波正演模拟也相对复杂.目前，对于煤层进行各向异性的研究主要集中在实验室分析［7－8］，波动方程法地震波模拟［9－14］及AVO特征分析［15－17］.射线追踪法和波动方程法各有优缺点.波动方程法相对来说，所生成的地震记录波形丰富，信息量大，比较接近实际地震记录；而射线追踪法适合于追踪单一波形，所得模拟地震记录可直接用于地震波各向异性响应特征研究.射线追踪法地震波模拟由于其具有显示直观、高效计算的特点，在石油地震勘探方面取得了显著成效［18－21］，然而却一直未见用于各向异性煤层中.试射法是一种基于射线的算法，也是一种经典的射线追踪算法，有很大的发展空间［22］.目前的地震记录正演模拟软件很多，如Tesseral等，它们大多是对层厚较大的模型进行模拟，还没有对VTI型煤层的正演模拟模块.

本文以各向异性煤层为研究对象，采用试射法对煤层中地震波进行射线追踪，分析VTI型煤层中走时及射线路径的传播特征；并将VC＋＋与Fortran语言有机结合，实现混合编程，研发基于VTI型煤层的地震记录正演模拟模块，将其作为一个开放式模块，可加载到其它软件系统中，提高正演模拟的精度和效率.

1 各向异性煤层中试射法射线追踪

相对于一般VTI型介质来说，VTI型煤层的结构更为复杂，射线方向的计算精度要求更高.本文将VTI型煤层的试射法射线追踪分为以下几个步骤进行：给定群角进行射线试射，进行群相关系转换，寻找射线出射角度，计算走时.

1.1 群相关系

各向同性介质中，相速度与群速度传播方向相同，大小相等，而在各向异性介质中，波的相速度指波前传播的速度，沿着相角方向传播；群速度是指能量的传播速度，沿着群角方向传播，决定地震射线的分布.群、相速度关系式为

式中v（θ）为相速度，v（φ）为群速度.

1.2 VTI型煤层中的Snell定律

由惠更斯原理可知，介质中的反射与透射均遵循Snell定律.在VTI型煤层中，介质对称轴一般不与界面法线方向一致，因此，其反射与透射关系并不能用简单的Snell定律来描述.根据图1并结合Snell定律，将VTI型煤层中的Snell定律表述为

图1 VTI介质反射、透射示意图Fig.1 Reflection and transmission in VTImedia

式中α1，α2，α3分别为 P波的入射角、反射角、透射角，v（θ1），v（θ2），v（θ3）分别为入射、反射、透射时的相速度.

1.3 射线角度和方向计算

1）初始角度

试射法射线追踪首先需要给定射线的初始角度.因为群角是与射线走向相关的，射线追踪时初始群角是可以设定的，给定初始群角φ的范围为：arctan（xmin／2／z）～arctan（xmax／2／z），其中xmin，xmax分别为最小炮检距和最大炮检距，z为反射面的深度.

2）根据群角计算相角

根据上式作f（φ）的近似式，通过递推得到相角θ.

3）反射相角、透射相角计算

由式（2）可知，反射、透射角度不能直接通过Snell定律的计算得出，因此，首先计算射线参数P，并给定一个射线收敛的精度范围，然后通过多次修改初始射线参数进行多次试射，使出射点与接收点接近.本文采用非线性方程求解方法——对分法来求反射相角、透射相角.

根据Snell定律（式（2）），假设f（θ）＝P－sinα／v（θ）＝0，则对于非线性方程f（θ）＝0，设其搜索区间为［0，π／2］，从端点θ0＝0开始，以步长λ逐步向后进行搜索．对于在搜索区间范围内的［θi，θi＋1］，分f（θi）＝0，f（θi）f（θi＋1） ＝0，f（θi）f（θi－1） ＞0，f（θi）f（θi－1） ＜0四种情形搜索方程的实根，从区间的左端点开始，直到区间的右端点结束.采用对分法求取方程的实根时，步长的选择不应过大或过小，选择过大的步长，会导致一些存在的实根的丢失，而选择过小的步长，则在一定程度上增加计算量，因此以较大步长但不丢失实根为宜.

1.4 走时计算

对于给定的射线参数，先由Snell定律通过计算得出射线在第i层的相角θi及与相角相对应的相速度vi（θ），再求出与相角相对应的群速度vi（φ）以及群角φi，则可以得到射线经过n层介质的PP波的射线走时为

式中，hi为第i层的厚度，v（φ）i为φi对应的群速度.

2 模型试算

表1为模型参数，第2层和第4层为VTI型煤层，其中第2层的顶、底界面起伏，但厚度保持不变，第4层界面起伏，厚度变化.各层均表现各向异性.模型大小为1 000 m×700 m.震源位于地表.

表1 地质模型参数Table 1 Parameters in geologicalmodel

采用观测系统：单边放炮，每炮60道，道间距10 m，最小偏移距150 m.为了在图形显示时，射线能被清晰地分辨，对射线按5∶1的比例抽稀.通过对模型进行射线追踪，所得射线路径如图2所示，图中煤层的顶、底界面反射、透射清晰.如果将各层的各向异性参数ε，δ均设为0，则煤层为各向同性煤层，射线追踪结果如图3所示.将上述两图对比可知：各向同性煤层与VTI型煤层中的射线路径有明显区别；没有出现射线盲区.图4将各向同性模型与各向异性模型追踪的走时进行对比，可以看出：地震走时记录中有四条同相轴，能有效分辨薄层；各向同性煤层与各向异性煤层的走时有较大差别，更进一步说明在煤层地震波的研究中，考虑煤层的各向异性给地震波传播带来的影响是必要的.

图2 VTI型煤层中PP波射线路径Fig.2 Ray path of PPwave in VTI coal

图3 各向同性煤层中PP波射线路径Fig.3 Ray path of PPwave in isotropicmedia

在INTEL P4 2.8 GHz的CPU和1GB的内存配置下，采用Fortran和VC＋＋混合编程自主开发VTI型煤层的射线追踪软件，并比较了4个反射界面各向异性和各向同性射线追踪的计算效率.如图5所示，VTI型煤层中的计算时间略高于各向同性介质的射线追踪时间，反射界面多的模型射线追踪计算时间略高于简单模型（大约7.6% ～12%），但是总耗时不超过30 s，表明计算速度较快.

3 结语

分析了VTI型煤层的射线追踪，采用VC＋＋和Fortran语言混合编程技术，实现了VTI型煤层的PP波射线追踪，并自主研发出实用的地震记录正演模拟软件（已申请软件著作权保护），只需改变速度和角度的相应参数，该方法还可适用于PSV转换波的记录模拟，实现多种波形的射线追踪，并可将射线追踪法推广到HTI型煤层，以及更为复杂的且与实际煤层更相近的TTI型煤层中；通过模型试算，对各向同性煤层与VTI型煤层走时进行对比，论证了煤层中的各向异性对地震记录有明显影响，模型模拟结果与理论分析一致，为煤矿地震勘探中各向异性分析方法的验证提供了重要依据.

图4 各界面反射到地表的走时对比（实线表示煤层为VTI型，虚线表示煤层为各向同性）Fig.4 Travel time of reflection waves（Solid lines∶in VTImedia；Dotted lines∶in isotropic media）

图5 计算时间的比较（无填充柱体表示煤层为各向同性，填充柱体表示煤层为各向异性）Fig.5 Comparison of calculation time（Column∶in isotropic coal；Filling column∶in VTI coal）

［1］David G.Seismic anisotropy in coal beds［J］.SEG Expanded Abstracts，2005：142－146.

［2］Arjun T，Shilpi G，Archana S，et al.Seismic azimuthal anisotropy：An important tool for detection of coal bed methane，an unconventional source of energy—A review［J］.Geohorizons，2010，12：23－27.

［3］Liu E R，Crampin S，Queen J，et al.Velocity and attenuation anisotropy caused by microcracks and macrofractures in multiazimuth reverse VSP［J］.Canadian Journal of Exploration Geophysics，1993，29（1）：177－188.

［4］Chen T J，Cui R F，Liu E R.AVO forward modeling for VTI coal［J］.Journal of China Coal Society，2009，34（4）：438－442.

［5］Wang J，Chen D H，Zhang H L，etal.Studies on phase and group velocities from acoustic logging［J］.Applied Geophysics，2012，9（1）：108－113.

［6］Wei Z T，Fan Y R，Chen X L，et al.Application of shear wave anisotropy in fractures and in-situ stress analysis［J］. Progress in Geophys，2012，27（1）：217－224.

［7］Dong SH.Test on elastic anisotropic coefficients of gas coal［J］.Chinese JGeophys，2008，51930：947－952.

［8］Liu Y Y，Wang Y，Zhang MG.Discussion on“Teston elastic anisotropy coefficients ofgas coal”［J］.Progress in Geophys，2012，27（4）：1832－1836.

［9］Yang D Y，Peng S P，Wang Y，et al.Seismic record simulation for coal-seam containing fractures［J］.Coal Geology＆Exploration，2007，35（6）：61－66.

［10］Pei Z L，Dong Y S，Peng SP.Seismic record simulation azimuthally anisotropy for coal-seam containing fractures［J］.Oil Geophysical Prospecting，2007，42（6）：665－672.

［11］Li D H，Dong S H，Zhao X C，et al.Seismic wave simulation of anisotropy in coal-bed media［J］.Geophysical＆Geochemical Exploration，2010，34（5）：604－609.

［12］Li D H，Dong SH，Zhao X C，et al.Seismic wave simulation of bipbase EDA media in coal-bed media［J］.Progress in Geophys，2011，26（2）：654－663.

［13］Qian Jin，Cui Ruofei，Chen Tongjun.Anisotropic numerical simulation of coal-bearing strata with finite-difference［J］.Coal Geology＆Exploration，2010，38（2）：63－67.

［14］Sun ZQ，Sun JG，Han FX.Traveltime computation using fastmarchingmethod under complex topographical conditions［J］.Chinese Journal of Computational Physics，2010，27（2）：281－286.

［15］Chen S C，Ma Z T，Wu R S.Wave equation method for wavefield propagation path tracing［J］.Chinese Journal of Computational Physics，2008，25（2）：186－194.

［16］Peng SP，Gao Y F，Yang R S，et al.Theory and application of AVO for detection of coalbed methane-A case from the Huainan coalfield［J］.Chinese JGeophys，2005，48（6）：1475－1486.

［17］Shi Y，Wang Y，Lu J.Analysis of azimuthal anisotropy in coalmeasures by 3CVSP data［J］.Journal of China Coal Society，2007，32（8）：813－817.

［18］Gao EG，Xu GM，Jiang X Y，etal.Segment ray tracing in 3Dmedia［J］.Oil Geophysical Prospecting，2002，37（1）：11－16.

［19］Zhao A H，Ding Z F.New approximate expressions of seismic group velocities for weakly anisotropy media［J］.Progress in Geophys，2005，20（4）：916－919.

［20］Li L.Phase velocity and group velocity in TImedia［J］.Geophysical Prospecting for Petroleum，2008，47（4）：334－338.

［21］Li JG，Li Y P，Guo X L.Shooting ray tracing in VTImedia［J］.Oil Geophysical Prospecting，2010，45（1）：491－496.

［22］Li Q，Bai C Y.Review on seismic wavefront and raytracing in complex media［J］.Progress in Geophys，2012，27（1）：0092－0104.

［23］Thomsen L.Weak elastic anisotropy［J］.Geophysics，1986，51（10）：1954－1966.

［24］Li Juan，Zhang Qin.Mixed programming of Visual C＋＋ based on dynamic link library［J］.Computer Engineering and Design，2010，31（9）：2140－2144.

［25］Liu JX，JiW，Liu Y，et al.Simulation and visualization of 1D fCSEMbased on Matlab and Fortran mix-programming［J］. Computing Techniques for Geophysical and Geochemical Exploration，2012，34（2）：224－228.

Seism ic Response of Anisotropic Coal Based on Ray Theory

LIQin1，LIQingchun2，ZHANG Lin3
（1.Xi'an University ofScience and Technology，College ofGeological and Environmental Sciences，Xi'an 710054，China；2.Chang'an University，College ofGeology Engineering and Geomatics，Xi'an 710054，China；3.Shaanxi Center ofGeological Survey，Xi'an 710016，China）

We illustrate a technique of ray tracing by shooting method and analyze rules of ray in VTI coal.With mixed-language programming of VC＋＋and Fortran，ray tracing in VTI coal with a complex velocitymodel is realized.It demonstrates necessaries of anisotropic ray tracing in VTI coal by comparison of ray tracing in VTI coal and isotropic media.By contrast of calculation time，it indicates that thismethod is efficient tomeet needs of seismic records forward.The technique can be used inmulti-wave ray tracing in anisotropic coal as giving responding formulae and parameters.

VTI coal；ray tracing；anisotropy

date：2013－07－27；Revised date：2013－11－25

P631

A

1001-246X（2014）04-0444-05

2013－07－27；

2013－11－25

国家自然科学基金（41304105，41374145）；高等学校博士学科点专项科研基金（20120205130002）；陕西省教育厅专项科研计划（12JK0470）资助项目

李勤（1979－），女，博士，讲师，从事应用地球物理信息技术工作，E-mail：eriliqin＠126.com