黏弹介质槽波超前探测数值模拟与波场分析

何帝秀，姬广忠，焦文杰，张亚伟，于 坤

(安徽理工大学 地球与环境学院，安徽 淮南 232001)

巷道掘进作为煤矿生产工作的重要环节之一，若煤巷掘进工作面(以下简称工作面)前方存在异常地质构造如断层、陷落柱、采空区等，很容易影响地下煤矿开采工作的安全性，严重时会引发矿井水害、顶板垮落、煤与瓦斯突出等地质灾害，给人员带来安全隐患，对煤矿生产造成经济损失。因此，采用物探方法进行工作面前方异常体超前探测工作尤其重要，常用于巷道超前探测的物探技术主要有反射地震波法、瞬变电磁法、地质雷达法等，而地震超前探测法中的槽波超前探测是利用仅在煤层中传播的槽波对工作面前方地质体进行超前探测，在探测断层、陷落柱、采空区等地质构造方面有优势。由于煤矿的全空间工作环境与工作面处的强反射波和绕射波干扰，造成了地震波的复杂性，因此，利用数值模拟方法研究含巷道煤层中反射槽波的传播规律，可以为实际巷道槽波超前探测提供一定的理论指导[1-5]。槽波数值模拟中通常将煤岩介质假设为完全弹性介质，而实际煤层更接近于黏弹性介质[6-7]。

近年来，诸多学者进行了槽波超前探测数值模拟研究，如丛皖平等[8]提出将Rayleigh 波运用于煤矿井下独头巷道超前探测；杨思通等[9]进行了Rayleigh 型槽波数值模拟并分析了工作面前方断层的波场特征，提出Rayleigh 槽波能量在垂直分量能量比水平分量更强；Yang Sitong 等[10]提出当震源位于巷道前方煤层中，且检波器沿着巷道中的煤壁布测时，使用反射Rayleigh 槽波超前探测煤矿巷道地质结构。蒋锦朋等[11]提出震源为垂向集中力源激发时，Rayleigh 槽波可应用于槽波超前探测；孙华超等[12]将Love 型槽波运用于小断层超前探测并分析其波场特征，并认为x分量的Love 槽波信噪比好；同年呼邦兵等[13]提出利用x和z分量的Rayleigh 型槽波对采空区进行超前探测；之后Sun Huachao 等[14]提出将x分量反射Love 槽波应用于陷落柱的超前探测；梁红波等[15]构建了4 个断层走向与巷道夹角不同的模型进行槽波超前探测模拟，提出夹角小于40°时偏移处理后的角度变化不大。

较多学者利用Kelvin-Voigt 模型研究黏弹性介质的地震波传播特征，例如苑春方等[16]研究了地震波在Kelvin-Voigt 均匀黏弹性介质中的传播速度和衰减系数；刘瑞珣等[17]发现采用Kelvin-Voigt 模型比其他黏弹模型更加符合地下介质是黏弹性质的实际情况；杨思通等[18]研究了Kelvin-Voigt 黏弹介质地震波衰减的影响因素，提出地震纵波振幅随波阻抗降低和震源子波频率增大衰减加快；严红勇等[19]研究了二维非均匀各向同性Kelvin-Voigt 黏弹性介质模拟，提出介质的黏滞性对转换波影响较大，且黏滞系数越大地震波能量衰减越大；姬广忠等[20]进行了Kelvin-Voigt 黏弹性TI 煤层介质Love 槽波三维模拟，同年张壹等[21]对Kelvin-Voigt 等典型黏弹性模型介质中的地震波衰减特征进行了归纳分析；焦文杰等[22]进行了Kelvin-Voigt黏弹性各向同性煤层介质槽波三维数值模拟。

目前黏弹介质下槽波超前探测的理论研究很少，为更好地模拟实际煤矿工作的槽波超前探测，笔者基于Kelvin-Voigt 黏弹性理论和黏弹三维一阶速度-应力方程，采用高阶有限差分方法进行独头巷道槽波超前探测的三维数值模拟，分析了煤层黏弹性及其变化、断层的落差和断层面与巷道的夹角等因素对反射槽波传播的影响，为实际槽波超前探测提供理论支持。

1 原理

1.1 槽波超前探测原理

槽波作为一种沿煤层传播的导波，传播距离远、能量强、具有明显的频散特性。根据不同的振动形式，通常将槽波分为两类，其一是Rayleigh 型槽波，由P 波与SV 波相互干涉形成，质点在与煤层面相互垂直、与传播方向平行的平面内振动；其二是Love 型槽波，由SH 波相互干涉形成的，质点在平行煤层面的平面内垂直于波传播方向上振动[23-24]。槽波大部分能量不向围岩辐射，携带大量地质构造信息，加之煤层的波导性，槽波被应用于超前探测[2]。槽波超前探测法为槽波反射法的一种特殊情况，两者原理类似。槽波超前探测原理为，激发震源后，槽波沿煤巷向工作面前方传播，若在超前方向遇断层面等反射界面，槽波会完全或不完全被阻断，产生反射槽波，由布置在巷道侧帮的检波器接收反射槽波信号，通过分析处理反射槽波记录，实现对工作面前方构造的探测[25-27]，如图1 所示。

图1 槽波超前探测原理Fig.1 Principle of in-seam wave advance detection

1.2 Kelvin-Voigt 一阶速度-应力方程

实际煤层具有黏弹特性，目前常用于地震正演模拟的黏弹性介质模型主要有Kelvin-Voigt 模型、Maxwell 模型、标准线性模型等。其中Kelvin-Voigt模型由一个弹簧元件和一个阻尼器元件并联构成[17]，如图2 所示。

图2 Kelvin-Voigt 模型Fig.2 Kelvin-Voigt model

Kelvin-Voigt 黏弹性介质本构关系为：

式中：σ为应力向量；σ1为 弹性单元体应力向量；σ2为黏弹性单元体应力向量；ε为应变向量；C为弹性矩阵；η为阻尼器黏滞矩阵；t为时间，s；Q为品质因子；ω为角频率。

用质点振动速度表示位移的一阶导数，来减少位移二阶导数的使用，简化整个计算过程。在无外力影响或外力影响消失后，根据Kelvin-Voigt 模型，获得三维黏弹一阶速度-应力方程[22]：

式中：σxx、σyy、σzz为正应力；τxy、τxz、τyz为剪应力；vx、vy、vz为质点振动速度；x、y和z分别为三维空间中的平面水平方向、平面竖直方向和纵向。

1.3 交错网格高阶有限差分法

应用交错网格高阶有限差分方法，将三维Kelvin-Voigt 一阶速度-应力方程进行离散化。高阶交错网格有限差分方法具有占用内存小、通用性强、高精度和高稳定性的特点。该方法是将一般的矩形差分网格替换为交错网格，借助交错网格的半程计算，来得到足够高阶的空间和时间精度的差分格式，以此来计算空间导数和时间导数。即在交错网格技术中，不仅要求空间网格的交错，在时间网格上也要求交错，变量的导数是在相应的变量网格点间的中点上计算[28-29]。

2 模型与参数

模型在x、y、z方向的尺寸分别为300 m×100 m×25 m，模型网格dx、dy、dz分别为1 m×1 m×0.25 m，煤层厚度为5 m。巷道在x、y、z方向的大小分别为150 m×4 m×4 m，位于x方向10～160 m、y方向48～52 m、z方向10.5～14.5 m，工作面位于x=160 m 处，断层位于x=240 m 处，如图3 所示。煤层顶底板岩性相同，巷道赋空值，模型各介质参数见表1，边界采用PML 吸收边界。

表1 模型介质参数Table 1 Medium parameters of the model

图3 含断层三维模型Fig.3 Schematic of the 3D model with faults

图4 为模型的观测系统，震源坐标为(110,54,12.5)，采用主频150 Hz 的雷克子波爆炸震源，延迟时间20 ms。检波器与震源同一水平线布置，过断层布测于x方向80～280 m，当落差大于1/2 倍煤厚时，断层面后的测线位于煤层顶板，道间距1 m，共201 个检波点。数值模拟采样间隔为0.05 ms。

图4 含断层xoy 切面图(z=12.5 m)Fig.4 xoy section with faults (z=12.5 m)

3 黏弹煤层介质对反射槽波的影响

3.1 煤层介质性质对反射槽波的影响

为研究煤层不同介质对反射槽波的影响，设计了2 个三维煤系含断层模型，分别是完全弹性介质模型和黏弹性介质模型(煤层QP=100，QS=50)，落差d=2.5 m(1/2 煤厚)，断面与巷道垂直(夹角α=90°)。

图5a、图5c 和图5e 分别是黏弹介质模型正演模拟合成地震记录的x分量、y分量和z分量，图5b、图5d和图5f 分别是完全弹性介质模型正演模拟合成地震记录的x分量、y分量和z分量，图中横轴为检波器编号，纵轴为传播时间。结合模型与时距关系分析可知，激发震源产生直达波，直达波沿巷道向前传播，在工作面处率先发生反射，形成工作面反射波，直达波继续向前传播遇断层产生断层反射波；直达波沿巷道向后传播，在巷道头处亦会发生反射，形成巷道头反射波。

图5 三分量合成地震记录Fig.5 Three components of synthetic seismograms

整体上看图5，各分量上相同波列在不同介质模型地震记录中的走时与分布相似，波振幅强度不同。图5a 和图5b 为分别是黏弹介质和完全弹性介质合成地震记录的x分量，依据时距曲线将易识别的波列按照被接收时间先后排序：直达纵波、直达横波、直达槽波、工作面反射波、反射横波、反射槽波和巷道头反射波。图5a、图5b 中，波组①由震源激发形成，出现时间最早，速度约为3 000 m/s，判断为直达纵波。波组②与直达纵波同相，振幅能量较强，速度约为1 700 m/s，判断为直达横波。波组③能量大，有频散现象，速度约为900 m/s，判断为直达槽波。波组④波列与直达波反相，出现较早，对其转换波组不进行分析，结合模型判断为工作面反射波。波组⑥由直达横波在断面处转换形成，速度约为1 700 m/s，判断为断层反射横波，简称为反射横波。波组⑦能量大，同相轴斜率比反射横波小，速度约为960 m/s，判断为断层反射槽波，为反射Rayleigh 槽波。在传播100 ms 后，出现与直达波同相的系列干扰波，结合模型判断其为巷道头反射波，对其转换波组不进行分析，统一标记为波组⑧。直达波遇断层面后，断面对不同地震波的阻断作用不同，透射纵波与透射横波能量衰减较小，槽波穿过断层后能量明显衰减；反射横波能量大，有清晰地震记录，而反射Rayleigh 槽波能量绝大部分衰减。完全弹性介质中各波列传播完整，断层反射波能量衰减小，且黏弹介质的反射横波波列延续时间比完全弹性的短。

图5c 和图5d 为分别是黏弹介质和完全弹性介质的y分量，依据时距曲线可以识别出：直达纵波、直达横波、直达槽波、工作面反射波、反射横波、反射槽波和巷道头反射波。波组⑥由直达横波在断面处转换形成，速度约为1 700 m/s，判断为反射横波。波组⑦波列呈现频散特征，波振幅大，速度约为900 m/s，判断为反射槽波，为反射Love 槽波。结合模型与时距关系，判断工作面反射波附近的双曲线干扰波为吸收边界反射波。直达波遇断层面后，透射横波与透射槽波能量衰减较小，透射纵波能量明显衰减；反射Love 槽波能量衰减较小，其他断层反射波能量相对衰减较大。完全弹性介质中各波列能量大于黏弹介质。

图5e 和图5f 为分别是黏弹介质和完全弹性介质的z分量，依据时距曲线可以识别出：直达纵波、直达横波、直达槽波、工作面反射波、P-S 波、反射横波、反射槽波和巷道头反射波。波组⑤由直达纵波在断面处转换形成，波列较宽，速度约为1 800 m/s，判断为直达纵波的转换横波，记为P-S 波。波组⑥与波组⑤相比斜率无明显变化，由直达横波在断面处转换形成，速度约为1 700 m/s，判断为反射横波。波组⑦波列呈现频散特征，能量大，速度约为900 m/s，判断为反射槽波，为反射Rayleigh 槽波。直达纵波遇断层产生反射时，另一部分能量转换为反射纵波，可能由于纵波波速较快导致反射纵波与工作面反射波混杂，未见明显同相轴。直达波遇断层面后，透射纵波能量增大，横波能量无明显变化，槽波能量衰减较大；在断面处产生的P-S 波、反射横波和反射Rayleigh 槽波能量一定程度上都减小。黏弹介质中反射Rayleigh 槽波能量明显随传播距离增大而衰减，完全弹性介质中无明显能量变化。

x分量上主要以反射横波为主，反射Rayleigh 槽波能量较小；y分量上反射横波波列范围大，反射Love 槽波能量强；z分量上P-S 波、反射横波和反射Rayleigh 槽波波形特征皆明显，信噪比好。与完全弹性介质模型相比，黏弹介质模型中地震波能量较小，地震波在黏弹介质中传播能量衰减速度比完全弹性介质快，导致部分波列在传播过程中衰减消失，没有完全弹性介质的波列完整。由此可知，黏弹介质模型更加符合实际煤层中地震波传播与衰减特征，更加贴近实际超前探测记录。

3.2 煤层品质因子变化对反射槽波的影响

为研究煤层不同品质因子对反射槽波的影响，设计了3 个不同煤层Q值的三维黏弹性煤系含断层模型，d=2.5 m，夹 角α=90°，QP、QS分别为40 与20、70 与35 和100 与50。

取测线第1 道至第81 道进行后续含断层模型的反射槽波对比分析。图6a、图6b 和图6c 分别为QP=40、QS=20 与QP=70、QS=35 及QP=100、QS=50 模型的x分量。图7、图8 分别为不同Q值模型正演模拟的y和z分量。根据时距关系3 个分量上皆可分辨出工作面反射波、反射横波、反射槽波和巷道头反射波，以及z分量的P-S 波。各分量上相同波列在不同Q值模型地震记录中的波形相似，波振幅强度不同。x分量上，反射Rayleigh 槽波和巷道头反射波能量随Q值增大明显增大，工作面反射波与反射横波能量随Q值增大而增大。当QP=40、QS=20 时，反射Rayleigh 槽波无明显同相轴。y分量上，各波列能量受Q值影响与x分量类似。随Q值增大，反射Love 槽波、工作面反射波、反射横波和巷道头反射波能量增大。当QP=40、QS=20 时，反射Love 槽波与反射横波能量小，波形模糊。z分量上，各波列能量受Q值影响与x分量类似。随Q值增大，反射Rayleigh 槽波和巷道头反射波能量明显增大；工作面反射波、反射横波能量增大；P-S 波能量变化较小。当QP=40、QS=20 时，反射Rayleigh 槽波与反射横波能量小。

图6 不同品质因子模型地震记录(x 分量，测线第1 道到第81 道)Fig.6 Seismograms of models with different quality factors (x component,the 1st to the 81th traces of the survey lines)

图7 不同品质因子模型地震记录(y 分量)Fig.7 Seismograms of models with different quality factors (y component)

图8 不同品质因子模型地震记录(z 分量)Fig.8 Seismograms of models with different quality factors (z component)

各分量上煤层Q值变化对不同地震波能量衰减影响不同，各波列的能量衰减都随Q值增大而减小。随Q值增大不同分量的工作面反射波和巷道头反射波能量增大；反射槽波在三分量上随Q值增大能量都明显增大；z分量的P-S 波受Q值影响小；反射横波在x分量受Q值影响小，在z分量受影响较大。QP=40、QS=20 时三分量上的反射槽波能量皆衰减大，传播时间越久能量越小，y分量上反射Love 槽波与反射横波能量都较小，因此，可利用x分量反射横波与z分量P-S 波进行超前探测。

4 断层变化对反射槽波的影响

4.1 断层落差变化对反射槽波的影响

为研究断层落差对反射槽波的影响，设计了3 个三维黏弹性煤系含断层模型，QP为100，QS为50，夹角α=90°，断层落差d分别为3.75、5 和7.5 m。

图9 为断层落差3.75、5 和7.5 m 模型的x分量，图10、图11 分别为不同落差模型正演模拟的y和z分量，根据时距关系3 个分量上皆可分辨出工作面反射波、反射横波、反射槽波和巷道头反射波，以及z分量的P-S 波。各分量上相同波列在不同落差模型地震记录中的走时相似，波振幅强度不同。x分量上，工作面反射波和巷道头反射波同相轴振幅随落差变化不明显，能量变化较小。随落差增大，断面对波的阻断作用增大，反射横波和反射Rayleigh 槽波能量增大。y分量上，各波列能量受落差影响与x分量类似。随落差增大，反射横波和反射Love 槽波能量增大。z分量上，各波列能量受落差影响与x、y分量不同。断层落差由2.5 m 增大至3.75 m 时，地震记录中的P-S 波和反射Rayleigh 槽波波形愈加明显，能量逐渐增大；落差由3.75 m 增大至7.5 m 时，直达波遇断面后转换波组能量迅速衰减，P-S 波和反射Rayleigh 槽波能量减小。落差2.5 m 增大至5 m，反射横波能量增大；落差5 m 增大至7.5 m，反射横波振幅变化较小，波列延续时间变短。

图9 不同落差模型地震记录(x 分量)Fig.9 Seismograms of models with different fault throws (x component)

根据上述分析可知，不同分量上地震波受落差变化的影响不同。工作面反射波和巷道头反射波在三分量上随落差增大能量变化都较小，说明断层落差对其影响小；在x和y分量上，断层落差越大断面对波的阻断作用越大，反射横波和反射槽波能量随落差增大而增大；z分量上，落差小于3/4 煤厚时，反射槽波、P-S 波和反射横波随落差增大能量增大；落差超过煤厚时，地震波的转换作用急剧减弱，由直达波转换形成的P-S 波和反射Rayleigh 槽波能量衰减大，反射横波能量变化较小。在断层上下盘完全错开后，z分量上反射Rayleigh 槽波能量小，可利用y分量的反射Love 槽波以及x和z分量的反射横波进行超前探测。

4.2 断层夹角变化对反射槽波的影响

为研究不同断层夹角对反射槽波的影响，设计了3 个三维黏弹性煤系含断层模型，QP均为100，QS均为50，d=2.5 m，夹角α分别为75°、60°和45°。

图12 为夹角75°、60°和45°模型的x分量，图13、图14 分别为不同夹角模型正演模拟的y和z分量，根据时距关系3 个分量上皆可分辨出工作面反射波、反射横波、反射槽波和巷道头反射波，以及z分量的P-S 波。各分量上相同波列在不同夹角模型地震记录中的走时不同，波振幅强度不同。x分量上，工作面反射波和巷道头反射波同相轴无明显变化，走时与能量变化较小；反射横波和反射Rayleigh 槽波同相轴斜率增大，波列间距减小，反射Rayleigh 槽波能量增大。y分量上，工作面反射波和巷道头反射波同相轴无明显变化，反射横波与反射Love 槽波被接收时间提前。断层夹角由90°减小至75°时，反射横波与反射Love 槽波能量增大，同相轴斜率增大；断层夹角由75°减小至45°时，反射横波与反射Love 槽波能量减小，波列间距明显减小。z分量上，各波列受夹角影响与y分量类似。P-S 波、反射横波和反射Rayleigh 槽波被接收到的时间提前，同相轴斜率随之增大，波能量呈现先增后减趋势。

图13 不同夹角模型地震记录(y 分量)Fig.13 Seismograms of models with different angles between a fault plane and a tunnel (y component)

图14 不同夹角模型地震记录(z 分量)Fig.14 Seismograms of models with different angles between a fault plane and a tunnel (z component)

图15 为夹角90°、75°、60°和45°模型的z分量波场快照在z=12.5 m 的切片图，结合模型与各波列速度可知，传播时间200 ms 切面上波列主要为巷道头干扰波、P-S 波、反射横波、反射Rayleigh 槽波、工作面绕射波和透射波。其中波组⑤的宽度最大，传播速度较快，判断为P-S 波；波组⑦波形最清晰，传播速度较慢，判断为反射Rayleigh 槽波；波组⑥传播速度介于波组⑤和波组⑦之间，波形较清晰，判断为反射横波；波组⑧与其他波列传播方向相反，产生于巷道头，判断为巷道头干扰波；波组⑨产生于工作面，波形混杂于断层反射波之中，判断为工作面绕射波。随夹角减小巷道头反射波波形无明显变化；P-S 波、反射横波和反射Rayleigh 槽波的反射角度随夹角减小而增大，测线位置的各波列清晰度在夹角小于75°后明显降低；工作面绕射波能量小于断层反射波，对波场影响较小，地震记录上难识别，对超前探测影响亦小。

图15 不同夹角模型200 ms 波场快照Fig.15 200 ms wave field snapshots of models with different angles between a fault plane and a tunnel

结合图14 地震记录与图15 波场快照分析可知，夹角角度越小，直达波离断层的相对距离越小，断层反射波产生的时间越早，反射角度越大。夹角由90°减小至75°时，检波器与断层面间距较大，反射槽波有足够的发育空间，各波列间隔较大，波形特征明显，且近似垂直反射，反射Rayleigh 槽波能量大；夹角由75°减小至45°时，测线与断层间距变小，各断层反射波出现时间明显提前，反射Rayleigh 槽波的反射角度增大，能量随之减小。P-S 波与反射横波同理。随夹角减小各断层反射波波列间隔减小，各波列难分离。

根据上述分析可知，不同分量上夹角对不同地震波的影响不同，且通过反射槽波的同相轴斜率可初步判断断层夹角大小。工作面反射波和巷道头反射波能量在3 个分量上变化较小，说明受断层夹角影响较小；随夹角角度减小，P-S 波、反射横波和反射槽波同相轴斜率增大。在x分量上，随断层夹角角度减小，反射横波与反射槽波被接收时间提前，能量增大；在y和z分量上，P-S 波、反射横波和反射槽波在夹角60°至90°模型中能量大，波形特征明显、易识别，槽波超前探测效果好；夹角小于60°时，反射槽波与其他波列难分离，探测分辨率随之降低，可利用x分量的反射Rayleigh 槽波进行超前探测。

5 结论

a.槽波在黏弹介质中能量衰减大，各分量上的反射槽波能量远小于直达槽波能量，z分量上反射Rayleigh 槽波随传播距离增大能量衰减速度明显大于完全弹性介质。说明黏弹介质中槽波的传播更加符合实际煤层槽波传播过程中能量的衰减特征。

b.在x、y和z分量上接收到的信号存在差异，x分量上主要以反射横波为主，反射Rayleigh 槽波能量较弱，不易识别；y分量上反射Love 槽波信噪比好；z分量上波列连续性较好，且P-S 波、反射横波和反射Rayleigh 槽波波形特征明显。

c.煤层Q值大小对反射槽波传播有较大的影响。随煤层Q值增大反射槽波能量增大，超前探测效果越好。QP=40、QS=20 时3 个分量的反射槽波能量皆小，同相轴都不明显，可利用x分量的反射横波与z分量的P-S 波进行超前探测。

d.断层落差与夹角在三分量上对反射槽波影响不同。断层落差小于煤厚时，y和z分量上反射槽波能量大，槽波超前探测效果都较好；落差大于煤厚时，z分量的反射Rayleigh 槽波能量小，可利用y分量的Love 槽波以及x和z分量的反射横波进行超前探测。断层夹角大于60°时，y和z分量上反射槽波能量大，槽波超前探测效果都较好；夹角小于60°时，x和y分量上反射槽波能量大，槽波超前探测都较好。

e.很多实际煤层Q值条件与QP=40、QS=20 断层模型相似，三分量上反射槽波能量都较小，槽波超前探测效果差，因此，可根据其他Q值断层模型数值模拟的分析结果，综合应用三分量上的反射槽波、反射横波和P-S 波，相互验证补充，提高煤巷槽波超前探测的精确度和可靠性。