地球物理学报  2012, Vol. 55 Issue (2): 645-654   PDF    
煤矿井下槽波三维数值模拟及频散分析
姬广忠1 , 程建远1 , 朱培民2 , 李辉2     
1. 中煤科工集团西安研究院, 西安 710077;
2. 中国地质大学地球物理与空间信息学院, 武汉 430074
摘要: 采用交错网格高阶有限差分法编制了地震波场三维正演模拟软件,设计了基于镜像法原理处理煤矿井下近水平和起伏巷道特殊空间的算法;模拟了煤矿井下含巷道和不含巷道情况下煤层中传播的地震波场,并分析其频散特征.结果发现:由于巷道的影响,巷道壁上产生很强的巷道振型槽波,煤层中则出现了以Love型为主的槽波,据此分析了实际槽波记录的形成机理,研究结果对今后煤矿井下巷道地震超前探测和工作面弹性波透视等具有重要的理论意义和实际价值.
关键词: 煤矿      槽波      三维数值模拟      巷道      交错网格高阶有限差分      频散     
3-D numerical simulation and dispersion analysis of in-seam wave in underground coal mine
JI Guang-Zhong1, CHENG Jian-Yuan1, ZHU Pei-Min2, LI Hui2     
1. Xi'an Branch of China Coal Research Institute, Xi'an 710077, China;
2. Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China
Abstract: The staggered-grid high-order finite difference method is used in developing of a 3-D seismic wavefield forward modeling software. Based on the principle of image method, the algorithm is designed to process the special spaces of subhorizontal and undulate roadway in underground coal mine. The seismic wave fields of 3-D coal seam models with and without roadway are simulated and dispersion characteristics are analyzed. It is found that strong roadway-type waves appear in the roadway wall and in-seam waves mainly of Love type appear in the coal seam. According to the simulation, we have analyzed the forming mechanism of practical in-seam wave records. The result has important theoretical significance and actual value for future seismic roadway probing ahead and coal working face perspective of elastic wave in underground coal mine..
Key words: Coal mine      In-seam wave      3-D numerical simulation      Roadway      Staggered-grid high-order finite difference      Dispersion     
1 引言

在煤炭开采过程中,常常因不能准确预测采掘前方的小构造及陷落柱等地质异常体而发生灾害事故.在地面勘探阶段,要发现这些小构造比较困难,因此需要发展煤矿井下地球物理探测技术.

目前煤矿井下巷道超前探测的地震方法有:槽波、瑞雷波、隧道地震勘探TSP等.煤矿井下槽波地震勘探方法的研究在20世纪80年代出现了一个高峰期,刘天放等(1994)编写的《槽波地震勘探》[1]一书,系统总结了国内外关于槽波的研究现状、槽波特性、模拟技术、处理方法和软件等.

对于煤层槽波的数值模拟,许多学者采用了多种方法模拟槽波,包括有限差分法[2]、有限元法[3]和解析法[4]等.Krey[5](1963)应用Haskell(1953)提出的方法,以Baumgarte和Krey[6](1961)的工作为基础,实现了理论频散曲线计算.1991年,Liu等[7]研究了煤层及围岩中裂隙对Love型槽波的影响.徐果明等[8]进行了瑞雷型槽波本征方程的研究.2001年,杨文强模拟了二维三层非对称煤层P-SV波的传播[9].杨真[10](2009)系统分析了薄煤层槽波频散特征,采用数值分析法求得基阶、高阶模式的频散曲线.

上述研究工作主要局限于槽波特性的二维数值模拟,对三维含巷道煤层的数值模拟还没有人做过.实际上,只有含巷道的三维地震地质模型才符合实际,所以需要研究三维含巷道煤层模型的数值模拟方法.

交错网格高阶有限差分地震模拟方法具有精度高、计算速度快的优点,但在处理起伏自由界面(如巷道)时较为麻烦,三维情况下做得较少.Robertsson[11]提出把不规则起伏表面分解成许多水平、竖直、顶角等类型网格单元的思路,把用于水平地表的镜像法推广到不规则界面,进行了二维起伏界面数值模拟.Wang等[12]改进了Robertsson的方法,采用不同的交错网格类型,简化了计算.汪利民[13]也按照这种思想,将AEA(Acoustic-ElasticboundaryApproach)[14]方法用于二维和三维的起伏自由表面模拟.其他处理方法有旋转交错网格法、映射网格法和坐标变换法等[15-16],但这几种方法不太适合煤矿井下三维巷道的处理.

本文采用交错网格高阶有限差分法[17]来模拟三维煤层槽波,边界吸收采用完全匹配层(PML)方法,着重研究含巷道特殊地质模型在数值模拟过程中的处理方法,通过分析比较现有自由界面数值模拟处理方法的优缺点,改进了镜像法,设计了适用于水平、倾斜和起伏不太剧烈的巷道三维镜像编程算法,并利用程序模拟含巷道和不含巷道三维煤层模型,来分析槽波特性.

2 三维巷道界面处理方法

巷道界面类似于地表自由界面,由于煤矿井下巷道有多个面,处理起来更加复杂.在有限元方法中,自由界面可以直接满足,但在交错网格有限差分方法中则较为麻烦.

2.1 自由界面处理方法

理论上,在自由界面处z=0的水平界面,自由边界条件为

(1)

σzzz方向正应力,τxzxz面上剪应力.

常用的处理自由界面的方法有三类:以真空形式处理自由界面(Heterogeneous Approach, 简称HA 法)、镜像法[12]和利用横向各向同性介质近似代替自由界面[18](由Mittet提出,简称Mittet法).Xu[14]对最后一种方法进行了改进,提出了AEA 方法.对于水平界面的处理,镜像法、AEA 法效果最好.

2.2 起伏界面及三维巷道处理方法

巷道壁是自由界面,可应用镜像法,对倾斜和起伏巷道借用Robertsson[11]的处理思想.但Robertsson方法将网格类型分类过多,三维情况下难于编程实现,本文将之改进为除去顶角网格等类型,只处理在自由界面(巷道壁)上的网格.有的网格可能有几个面是自由界面,对每个面分别设置镜像法.交错网格选择Wang[12]的类型.

按照以直代曲、变换交错网格类型、每个自由网格面逐个设置镜像法的思路,本文将传统的镜像法加以改进,来模拟三维起伏巷道.具体实现如下:

采用将正应力放置在网格中心的交错网格,正应力不在自由界面上,处理起来简单(图 1表 1).

图 1 三维交错网格节点分布 Fig. 1 Nodes distribution of 3-D staggered-grid
表 1 交错网格各节点处的弹性参数 Table 1 Elastic parameters at nodes of staggered-grid

表 1中,σ 为正应力,τ 为剪应力,V为质点振动速度,λμ 为拉梅常数,ρ 为密度,以下符号含义同此.

以水平地表自由界面为例,设水平自由表面位于z=kΔz,空间阶数为4,三维镜像法的公式是

(2)

式(2)是相对地表而言,巷道底板与此相同,但巷道顶板与地表朝向相反,对离散网格来说地表界面设置在上表面,对巷道顶板则自由界面设在网格的下表面,镜像方程发生变化,设巷道顶板面位于z=kΔz,(2)式应修改为

(3)

与巷道底板不同,质点振动速度Vz(ijk+1)因在自由界面上(图 1 节点4)需要计算,但此离散网格的物性值为0,计算时密度需用上面网格的密度.巷道侧壁处理方法同上.

对每个自由网格面采用上述镜像法,巷道内物性为0的网格可用于设置镜像应力值(公式(2)、(3)左边变量).

2.3 三维巷道模拟编程算法

巷道处理思路简单,关键在于编程.巷道壁上网格有多种类型,镜像公式不同,增加了编程难度,需要设计周密、精巧的算法,才能保证巷道界面的正确处理,且避免增加过多计算量.

以在xoz面的三维网格(图 2)为例,斜线标注的是需要镜像法计算的网格,阴影是巷道壁自由界面.以网格A 为例,设空间阶数为4,计算质点振动速度需用到周围各个方向两个网格的应力值,A 所用网格xz方向各有一网格在巷道内,需设置镜像条件.巷道网格应力值,设置为与自由界面相对称的介质网格应力的负值(公式(2)、(3)左边变量),具体过程如下:

图 2 巷道处理示意图 Fig. 2 Schematic diagram of roadway processing

(1) 根据物性参数判断网格类型,包括介质、巷道、巷道壁类型,做上标记.

(2) 根据网格与巷道自由界面的距离,判断需用镜像法的网格,做上标记.

(3) 在循环计算中,对正常介质网格采用正常算法,巷道网格不计算,对需用镜像法的网格在计算质点振动速度之前,对需用到的巷道网格设置应力镜像值:

(a) 在每个方向上进行搜索,找出(1)中做了标记的巷道壁网格,即可确定自由界面位置,然后设置镜像条件(正负方向镜像公式不同).这样解决了网格有多个面是自由界面的问题.

(b) 对有些巷道壁网格,以图 2 中网格A 右边的网格为例,网格右面是自由界面,振动速度Vy节点位置在面中心(图 1中节点3),但Vy对应的是旁边巷道网格密度值0,而Vy是有值的,所以用为0的密度值计算不正确,计算时需用巷道壁网格的密度.在对此自由界面设置镜像条件时,可将Vy一同算出.其他类型同上.

(c) 设置完后,仍按原来算法顺序计算A 的质点振动速度.

此算法只增加了很小的内存和计算量,保证了程序的高效率,对起伏地表和隧道同样适用.需注意的是:网格自由界面前方须有L个(空间阶数2L)巷道网格,保证镜像条件设置;高阶容易发散,低阶比较稳定,对规则巷道空间阶数可选用4或6为宜,对倾斜或起伏巷道应降低阶数,可以选用2或4阶.

3 煤矿井下三维地质模型地震波场数值模拟及频散分析 3.1 模型1:三层对称煤层模型(煤层厚5m)

模型1 中xyz方向的大小为200 m×200 m×50m, z方向高度50m, 中间为煤层,煤厚5m, 两边是岩性相同的围岩,xyz方向网格大小1m×1m×0.5m, 煤层纵波速度1700m/s, 横波速度1000m/s, 密度1300kg/m3,围岩纵波速度3000 m/s, 横波速度2000m/s, 密度2200kg/m3.时间采样间隔dt=0.1ms, 炮点位置(x=4,y=100,z=25)(图 3),纵波震源雷克子波主频120Hz.

图 3 三层对称煤层模型(模型1) Fig. 3 A geological coal model with three symmetrical layers (Model 1)

地震波场快照和合成记录如图 4所示.槽波在煤层中干涉形成,能量集中,占总能量的大部分比重,最前面是折射波,看不到直达波,质点振动速度Vz分量在围岩中能量泄露较多,Vx则很少.

图 4 地震波场快照(模型1) Fig. 4 Seismic wavefield snapshots (Model 1)

选取过炮点位置沿y方向在z=25m, y=100m处测线(图 3所示)合成记录(图 5)做频散分析(图 6).

图 5 测线(z=25m,y=lOOm)合成记录 Fig. 5 Synthetic traces of survey line (z=25 m, y=lOO m)
图 6 测线V-f域功率谱 Fig. 6 Power spectrum in V-f domain of survey line

合成记录中,槽波为一串波列,波场外形呈斜的窄条带状,近似一个斜率,绝大部分能量分布在这一区域,变化区间小,速度稳定,对实际勘探有较大应用价值.

应用FK 法提取频散曲线(图 6).5m 煤层槽波频散曲线能量Vx分量主要集中在100~200Hz, Vz集中于150Hz附近,相速度约1500~1600m/s, 这是优势频带,其他范围能量很弱,所以从图中不能看到完整的频散曲线形状;频散曲线呈稍微倾斜的直线状,槽波能量团传播速度十分稳定,为利用槽波埃里相勘探提供了理论依据.由于是纵波震源,模拟的槽波属于瑞雷型槽波.

3.2 空巷模拟算法检验

为检验含巷道条件下的算法效果,设计了含水平长方体空巷和倾斜空巷的均匀介质模型,空巷物性值为0.含水平空巷模型大小100 m×100 m×100 m, 网格大小1 m×1 m×1 m, 时间采样间隔dt=0.1ms, 纵波震源主频120 Hz, 空间阶数为6,介质为煤,参数与模型1的煤层参数相同.含倾斜空巷模型大小100m×100m×100m, 在xz面上倾斜18.5°,y方向水平,网格大小x=0.5 m、y=1m、z=0.5m, 时间采样间隔dt=0.1ms, 纵波震源主频100Hz, 空间阶数为2.震源位置都在空巷壁附近.模拟实验结果如图 78所示.

图 7 水平空巷40 ms波场快照 Fig. 7 Wavefield snapshots of horizontal and empty cuboid at 40 ms
图 8 倾斜空巷40 ms波场快照 Fig. 8 Wavefield snapshots of dip and empty cuboid at 40 ms

水平空巷模型网格的剖分稀疏,最小波长占8~9个网格,接近差分模拟最低要求,模拟结果(图 7)符合瑞雷面波的特征.空巷两侧壁波场传播正常,说明在网格不同位置的自由面,设置不同镜像条件的模拟方法是正确的.

倾斜空巷模型网格剖分较密,xz方向最小波长占16~17个网格,y方向最小波长占8~9 个网格,模拟结果较好,波形轮廓符合瑞雷面波特征(图 8).

3.3 模型2:含巷道煤层模型(煤厚5m)

模型xyz方向大小200 m×200 m×50 m, 在模型1的基础上,增加两条y方向平行巷道,巷道截面4m×4m, 与边界距离20m, y方向延伸长度160m, 与竖直切片对称分布,炮点在左边巷道壁旁煤层里,位置为(x=26,y=100,z=25),如图 9 中圆圈所示.绿线表示测线,测线1(x=24,z=25)在左边的巷道壁上,测线2(y=100,z=25)为过炮点与测线1垂直的绿线(图 9).

图 9 含巷道煤层模型(模型2) Fig. 9 A geological coal model including roadway (Model 2)

网格大小x=1 m、y=1 m、z=0.5 m, 时间采样间隔dt=0.1ms, 纵波震源雷克子波主频120Hz, 巷道模拟采用上述镜像法,空间阶数为4.

3.3.1 60ms波场快照

选取60ms的波场快照(图 10).

图 10 60 ms波场快照(模型2) Fig. 10 Wavefield snapshots at 60 ms (Model 2)

含巷道模型波场(图 10)与无巷道波场(图 5)差别很大.巷道壁是自由界面,产生能量很强的面波及转换横波,面波在顶底板间相互干涉,形成一种更加复杂的巷道振型槽波.煤层内则出现很强的由横波及其转换波在煤层内干涉形成的槽波.波场能量主要集中在煤层内,围岩中能量小没有显示出来.煤层中整个波场形状类似二维瑞雷面波,前部分是纵波及转换波干涉形成的瑞雷型槽波,后面是横波及转换波形成的槽波,巷道壁上是面波形成的槽波,而无巷道模型波场只有前部分.

由波场快照可说明槽波的形成过程为:震源在巷道壁附近激发,产生瑞雷面波和纵横体波,在薄煤层内形成槽波,由于煤层薄,速度、密度低,大部分能量被挤压在煤层内,所以整个波场外形类似二维瑞雷面波.

震源在巷道壁产生了面波和转换SV 波,但在此三维坐标系中,炮点所在巷道壁是竖直方向,大部分SV 波在平行煤层的水平方向上传播(尤其当槽波传播一段距离后),在坐标系中大部分“SV"波变成了“SH"波,则横波型槽波以SH 波为主,属Love型槽波,与实际接收到的槽波记录相符;巷道振型槽波由P波和“SH"波(相对于巷道壁是SV 波)为主的面波在顶底板间干涉形成,并包含有少量SV 波,它不同于Love型和Rayleigh 型槽波,性质更加复杂,目前尚无理论论述.

通过本次三维模拟,我们论证了巷道在SH 型槽波生成的作用,能够较合理地解释实际槽波记录的形成原因.另外,如果震源位于煤质松散、不均匀性强的位置,产生的震源波场成分则十分复杂,纵波能量不一定占大部分比例,SV 波可能也很强,形成的波场可能有较大不同,需要继续深入研究.

3.3.2 合成记录及频散分析

取测线1、2 各分量合成记录(图 11、12),并做频散分析(图 13、14)

图 11 测线l(x=24,z=25)各分量合成记录 Fig. 11 Various component synthetic traces of survey line No.1 (x=24,z=25)
图 12 测线2(z= 25,y=100)各分量合成记录 Fig. 12 Various component synthetic traces of survey line No.2 (z= 25,y = 100)
图 13 测线l(x = 24,z=25)各分量V-f域功率谱 Fig. 13 Various component power spectrum in V-f domain of survey lineNo.l (x= 24,z= 25)
图 14 测线2(z=25,y=100)各分量V-f域功率谱 Fig. 14 Various component power spectrum in V-f domain of survey line No.2 (z=25 m, y = 100 m)

测线1合成记录中,槽波能量团呈条带状,传播速度稳定,Vx分量(图 11a)巷道振型槽波占绝大部分;Vy分量(图 11b)分为两部分,巷道振型槽波能量强,能量团群速度约800 m/s, 瑞雷型稍弱,群速度约1500m/s, 低于煤层纵波速度1700m/s;Vz分量(图 11c)瑞雷型槽波能量比巷道振型的强.频散图(图 13)与合成记录图(图 11)相对应,Vx分量(图 13a)基本上是巷道振型槽波频散曲线,Vy分量(图 13b)以巷道振型槽波频散曲线为主,瑞雷型频散曲线与无巷道模型(图 6a)形状相同,Vz分量(图 13c)瑞雷槽波能量最强,但能量团较散,没形成清晰的曲线形状,巷道振型相对较弱.巷道振型槽波频散曲线趋向的截止速度低于煤层横波速度1000m/s, 符合地表瑞雷面波特征.

测线2Vz合成记录(图 12c)明显有两种模式瑞雷型槽波,上面槽波以纵波能量为主,下面以SV 波为主.Vx分量(图 12a)纵波型槽波能量很强、SV 型很弱,Vy(图 12b)主要是SH 型槽波.频散图(图 14)与之对应,值得注意的是:Vz分量(图 14c)SV 型槽波频散曲线形状与Vy分量(图 14b)相似,但位置向后推移了一段距离,与纵波型频散曲线最大能量团相连.合成记录中还包括有其他类型的波,如直达波、折射波、多次反射波等(直达波能量小,湮灭在其他波中),波场复杂,频散严重,在频散曲线周围(图 14)出现了阴影区域.

4 结论及建议

本文研究了含巷道特殊地质模型在数值模拟中的处理方法,设计了适用于水平、倾斜和起伏不太剧烈巷道的三维镜像编程算法,首次进行了含巷道煤层模型三维槽波数值模拟.通过对含巷道和不含巷道煤层模型的模拟,得到以下认识:

(1) 含巷道和不含巷道煤层波场差别很大,含巷道模拟结果符合实际情况.巷道壁上产生很强的巷道振型槽波,波场成分特殊,相对于巷道壁以P、SV波为主,在坐标系中以P、SH 波为主,目前理论描述很少.

(2) 通过分析含巷道煤层模型的模拟结果,论证了实际数据中Love型槽波形成机理.由于巷道的存在,纵波震源在巷道壁上产生瑞雷面波及转换横波,转换SV 横波在坐标系中为水平方向,变成SH波,所以实际记录中出现了能量很强的Love型槽波.

(3) 巷道振型槽波和地表瑞雷面波有明显差异,目前研究很少,在基础理论上还需深入研究.

文中不足之处包括:实际煤层是弱各向异性,各向同性只是近似情况,模拟结果和实际会有一些出入;巷道内有空气介质,对波场稍有影响,文中近似当作真空处理;对实际槽波形成机制的解释限于数值模拟分析,可能产生认识偏差.煤层槽波十分特殊,很多地面勘探方法和处理技术不适用于槽波,迫切需要研究槽波性质和传播规律,为实际勘探提供理论基础.

参考文献
[1] 刘天放, 潘冬明, 李德春, 等. 槽波地震勘探. 徐州: 中国矿业大学出版社, 1994 . Liu T F, Pan D M, Li D C. In-Seam Seismic Exploration (in Chinese). Xuzhou: China University of Mining and Technology Press, 1994 .
[2] Korn M, Stckl H. Reflection and transmission of Love channel waves at coal seam discontinuities computed with a finite difference method. J. Geophys. , 1982, 50: 171-176.
[3] Edwards S A, Asten M W, Drake L A. P-SV wave scattering by coal-seam inhomogeneities. Geophysics , 1985, 50(2): 214-223. DOI:10.1190/1.1441911
[4] Buchanan D J. The scattering of SH-channel waves by a fault in a coal seam. Geophysical Prospecting , 1986, 34(3): 343-365. DOI:10.1111/gpr.1986.34.issue-3
[5] Krey T C. Channel waves as a tool of applied geophysics in coal mining. Geophysics , 1963, 28(6): 701-714.
[6] Baumgarte J, Krey T. Reflexion und Brechung Beim Schragen Durchgang Ebener Seismisgher Wellen durch N Planparallele Medien. Geophysical Prospecting , 1961, 9(2): 242-260. DOI:10.1111/gpr.1961.9.issue-2
[7] Liu E R, Crampin S, Roth B. Modelling channel waves with synthetic seismograms in an anisotropic in-seam seismic survey. Geophysical Prospecting , 1991, 40(5): 513-540.
[8] 徐果明, 倪四道, 王汉标. 瑞利型槽波的本征方程及其应用. 煤炭学报 , 1998, 23(2): 124–129. Xu G M, Ni S D, Wang H B. Eigen-equation of Rayleigh guide waves and its application. Journal of China Coal Society (in Chinese) , 1998, 23(2): 124-129.
[9] 杨文强. 槽波地震勘探的数学物理模型研究. 地质与勘探 , 2001, 37(3): 58–60. Yang W Q. Modeling research of channel wave seismic exploration. Geology and Prospecting (in Chinese) , 2001, 37(3): 58-60.
[10] 杨真. 基于ISS的薄煤层采空边界探测理论与试验研究. 中国矿业大学, 2009. Yang Z. Theory of using an In-Seam Seismic (ISS) method for uncertain boundary detection in a thin coal seam (in Chinese). China University of Mining and Technology, 2009.
[11] Robertsson J O A. A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography. Geophysics , 1996, 61(6): 1921-1934. DOI:10.1190/1.1444107
[12] Wang X M, Zhang H L. Modeling of elastic wave propagation on a curved free surface using an improved finite-difference algorithm. Sciences in China Ser. G Physics, Mechanics & Astronomy , 2004, 47(5): 633-648.
[13] 汪利民. 三维带地形瑞雷面波交错网格有限差分法正演技术研究. 武汉: 中国地质大学, 2009. Wang L M. Modeling of Rayleigh-wave propagation in three-dimensional media with topography using staggered-grid finite difference scheme (in Chinese). Wuhan: China University of Geosciences, 2009.
[14] Xu Y X, Xia J H, Miller R D. Numerical investigation of implementation of air-earth boundary by acoustic-elastic boundary approach. Geophysics , 2007, 72(5): 147-153. DOI:10.1190/1.2753831
[15] Bohlen T, Saenger E H. Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves. Geophysics , 2006, 71(4): 109-115. DOI:10.1190/1.2213051
[16] Hestholm S, Ruud B. 2D finite-difference elastic wave modelling including surface topography. Geophysics Prospecting , 1994, 42(5): 371-390. DOI:10.1111/gpr.1994.42.issue-5
[17] 牟永光, 裴正林. 三维复杂介质地震数值模拟. 北京: 石油工业出版社, 2005 . Mou Y G, Pei Z L. Seismic Numerical Modeling for 3-D Complex Media (in Chinese). Beijing: Petroleum Industry Press, 2005 .
[18] Mittet R. Free-surface boundary conditions for elastic staggered-grid modeling schemes. Geophysics , 2002, 67(5): 1616-1623. DOI:10.1190/1.1512752