基于Morlet小波的时频分析方法在地震槽波精准反演煤厚中的应用

郭昌放，杨 真，武 祥，陈一鼎，顾 芗

(1.中国矿业大学人工智能研究院，江苏省徐州市，221116；2.中国矿业大学矿业工程学院，江苏省徐州市，221116；3.昆明煤炭科学研究所，云南省昆明市，650051；4.方圆标志认证集团，北京市海淀区，100048)

近些年来，我国袁亮院士[1]以及澳大利亚联邦科学与工业研究组织[2]等国内外众多学者与机构曾先后提出，工作面透明化作为煤矿智能精准开采实现的基础和前提，对推动煤炭资源的安全高效开采具有重要意义。而煤矿工作面透明化建设的核心是如何精准识别工作面内部的地质异常环境[3-4]。同时，在回采过程中，煤厚的突然变化以及工作面内部隐伏的地质构造经常会造成工作人员伤亡和巨大的经济损失[5-6]。因此，在工作面回采之前，煤厚及地质构造的精准预测对煤炭安全、智能、精准开采具有重要意义。

目前工作面内煤厚及地质构造探测方法已在煤矿得到广泛应用。其中，常见的钻探方法[7]对于单个钻孔对应的煤厚数据是准确的，然而受成本和时间限制，探测区域钻孔数量有限，存在探测盲区；三维地震勘探[8]是在地表全方位布设炮点和检波点，进行异常面积观测的一种地震勘探方法，然而由于探测距离较大而且易受上覆岩层和地形条件的影响，探测精度较低；无线电波坑透技术[9]是煤厚探测的常用物探方法，然而对于内部起伏较大的工作面探测效果较差，而且探测结果只能用于定性判断。

槽波地震勘探技术由于其探测距离大、精度高、抗电干扰能力强以及波形易于识别[10-11]等优点，成为目前国内外煤矿工作面内煤厚及地质构造探测和反演的主要方法[12-13]。其中，地震槽波走时的提取对煤厚及地质构造反演结果的精准性至关重要。笔者基于Morlet小波的高分辨性，利用时频分析方法将地震槽波在煤层中的传播以能量流动的形式进行分析，求取地震波信号的时间-能量谱，并拾取地震槽波的走时。基于Morlet小波的时频分析方法在邯郸矿区郭二庄煤矿22204工作面槽波地震勘探过程中，预测结果与回采揭露结果基本保持一致，得到了良好的验证。

1 槽波走时分析方法

地震波在煤层中激发产生之后，经过顶底界面的多次全反射以及相互干涉形成槽波。槽波在煤层中传播的过程中如果遇到断层、陷落柱以及煤厚变化等异常，槽波传播的速度、相位等参数就会发生变化。因此，槽波地震勘探方法就是利用槽波对异常体的高敏感性从而对煤层中的地质异常区域进行预测和解释[14]。

频散是地震槽波最基本的特性，槽波传播的速度主要取决于频率的大小。因此，槽波同时具有相速度和群速度，然而拾取单一相速度走时比较困难，而且误差较大。因此，通常选取槽波群速度的走时进行层析成像反演计算。在槽波走时分析过程中，通常是根据煤厚的探测范围，利用某一固定频率在频散曲线上拾取槽波的群速度和对应的走时[10,13]。然而，频率的选择并不具有唯一性，而且，只适用某一固定的煤厚探测范围。因此，对于较大范围的煤厚变化，必须选取不同的频率重复拾取走时，工作量较大。另一方面，频散曲线的不规整性，对槽波走时的精准拾取带来不利影响。为了解决上述问题，提高槽波走时拾取的准确性，采用时频分析方法获得透射槽波传播的时间-能量谱，以槽波能量包络中心点对应的时间作为槽波的走时，为地质异常的精准反演奠定了基础。

在地震勘探中，由于受到多次反射、地层吸收以及频散的影响，地震信号往往具有非平稳特征[15]。常规傅立叶变换无法表征任一时刻的频率成分，而时频分析作为分析非平稳时变信号的主要工具[16]，能够描述频谱随时间的变化情况，并可以使得地震波信号的能量密度分布在时域和频域中同时显示。其中，小波变换作为一种时间尺度域分析技术，广泛应用于地球物理勘探数据分析中[17]。通过变换能够对时间(空间)和频率进行局部化分析，通过伸缩平移运算对信号逐步进行多尺度细化，最终达到高频处时间细分、低频处频率细分，因此能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节，是一种多分辨率的分析方法。

小波分析得益于小波基函数的完备性、自相似性和多分辨性,母小波选择不当，应用效果将会受影响[18]。由于Morlet小波在时域和频域上具有较好的分辨率，它是一种复杂的单频正弦信号调制高斯波:

式中：ψσ(t)——Morlet小波随时间变化的波形振幅曲线；

t——时间；

Cσ——归一化常数；

σ——无量纲频率；

κσ——小波容许准则。

Morlet小波ψσ(t)对应的傅立叶变换表达式为：

(4)

w——频率。

(ωψ-σ)2-1=(ωψ2-1)e-σωψ

(5)

其中，参数σ决定了时域和频域之间的关系。一般而言，由于波形频率随时间变化较快，σ的取值往往大于5。因此κσ的取值也接近于零，而且中心频率ωψ≅σ。方程(1)也可以简化为：

(6)

由于σ为无量纲频率，归一化常数Cσ=1，则方程(4)可以简化为：

(7)

为了保证Morlet变换在每一个缩放因子s下具有可比性，需要对小波函数进行归一化处理，使其具有相同的单位能量：

ωk——第k个频率分量；

N——频率索引。

2 工程应用

2.1 地质概况

郭二庄煤矿22204工作面煤层结构简单，煤层厚度0.20～2.20 m，平均1.48 m，平均倾角23°，两巷落差近60 m，由于两巷落差较大，因此选择槽波地震勘探方法。在工作面两巷掘进过程中共揭露了10条断层，位于探测区域范围内的断层共有9条，均为正断层。22204工作面探测区域如图1所示。

图1 22204工作面探测区域

2.2 观测系统和数据采集

本次勘探以炸药爆破作为震源，通过单分量检波器采集地震波，炮孔深度为3 m，药量为300 g，发射点和接收点尽量布置在煤层中部。探测区域测点布置如图2所示。所有的发射点(S1～S31)布置在运输巷，间距10 m，要求炮孔直径为42 mm，尽量保持与煤层平行。所有的接收点(R1～R69)都固定在锚杆上，间距10 m。数据采集时，炮点从S31到S1逐一爆破，产生地震波并被检波器接收。在数据采集过程中，由于检波器R1、R2、R3和R69出现故障，因此共采集到2015道震波数据，图2中给出了3个炮点对应接收点的射线路径。

图2 探测观测系统图

在22204工作面透射槽波勘探中，采用以Morlet小波为母小波的小波变化时频分析方法来拾取槽波走时。震源点为S5时，在接收点R7处采集到的地震波数据时频分析结果如图3所示。从图中能够清楚看到槽波信号同时在时域和频域中的能量分布情况，并以槽波能量包络中心点对应的时间作为槽波的走时，对采集的2015道地震波数据进行时频分析，共拾取了1750个有效槽波走时数据，为进一步的层析成像反演奠定了基础。

图3 震源为S5时接收点R7采集到的地震波数据时频分析结果

2.3 层析成像反演

槽波层析成像(CT)技术已广泛应用于工作面内部的透射投影，利用槽波通过非均匀介质传播的时间，对其内部慢度分布的结果进行重构[19]。运用CT技术进行反演计算的常用算法有：最小平方法(LSQR)、代数重建法(ART)以及瞬时迭代法(SIRT)，由于SIRT对于稀疏、不规则以及低信噪比的探测数据具有较好的适用性，而且比较稳定[20]，因此选用SIRT算法用于层析成像反演。

(10)

对探测区域进行网格离散化，如图4所示，划分为N个单元，每个单元内的慢度值分别为S1,S2,S3,…,SN。公式(10)可转变为：

(11)

式中：tk——第k条射线走时；

dki——第k条射线经过第i面元的长度；

Si——第i面元的慢度。

若射线的条数为M，则所有射线路径走时方程可转化为：

T=D·S

(12)

式中：T——M维时间矩阵；

S——N维慢度矩阵；

D——(M*N)矩阵算子。

利用瞬时迭代重构法(SIRT)求解式(11)，从而可以得到慢度反演结果。

基于采集的地震槽波走时，通过SIRT算法进行层析成像迭代反演，由于网格的尺寸和数目都会影响层析成像的结果，因此在反演之前，需要对这些参数进行多次测试。最终网格尺寸定义为10.3 m × 10.3 m，共23行、61列。通过多次迭代计算，最终得到了探测区域槽波传播慢度分布结果，如图5所示。

图4 探测区域网格离散化分析

图5 探测区域慢度分布结果

2.4 慢度与煤厚之间的关系

槽波传播的慢速与传播的介质往往具有较好的对应关系，而且槽波在煤层中传播的慢度比在岩层中传播的值要高。图5显示了探测区域内部的槽波传播慢度分布结果，其中蓝色低慢度区域主要分布在400～680 m之间，初步判断该区域的煤厚值比周围的要低。然而，通过慢度分布只能对探测区域内煤厚的变化情况进行定性判断。

为了更清楚地了解探测区域煤厚的定量分布，针对22204工作面回风巷200～500 m以及运输巷400～750 m区段揭露的煤厚结果(0.2～2.2 m)和对应的慢度值进行数值拟合，由于试验数据近似于线性分布，因此选择一阶多项式拟合煤厚-慢度关系，如图6所示。

图6 煤厚-慢度拟合结果

3 煤厚分布结果

利用一阶多项式拟合的煤厚-慢度线性关系(图6)，将慢度分布结果(图5)转换成煤厚分布结果，如图7所示。22204工作面探测区域在700～900 m之间，由于受巷道布置限制，射线覆盖比较稀疏，会对结果的准确性造成影响。

由图7中可看出，探测区域煤厚中间薄、两端厚。其中两端平均煤厚都在1.5 m以上，而中部煤厚变化在0.5～1.5 m之间，由于该区域共揭露了6条正断层，判断由于受断层影响，顶底板岩层在张力的作用下延伸至煤层导致该区域煤层变薄。而且该薄煤区域延伸范围较大，对工作面回采影响较大。

22204工作面沿着回采方向生产过程中，在开采到706 m的位置时，由于煤厚太薄，无法继续回采，最终在470 m位置新开切眼，自右向左继续完成回采工作。因此经过生产验证该探测结果与实际相符，能够为工作面合理布置提供参考依据。

图7 探测区域煤厚分布结果

该探测区域回采完成之后，通过揭露的煤厚与预测煤厚之间的相关系数r进行对比验证，相关系数r定义为：

(12)

式中：D——揭露的煤厚结果，m；

P——预测的煤厚结果，m；

Cov(D,P)——揭露煤厚与预测煤厚之间的协方差；

Var[D]——揭露煤厚的方差；

Var[P]——预测煤厚的方差。

22204工作面回采后揭露的结果与预测结果进行对比分析，其结果如图8所示。揭露的煤厚与预测的煤厚之间的相关系数r为0.88，具有较好的吻合效果。

4 结论

基于Morlet小波的时频分析方法在郭二庄煤矿22204工作面煤厚探测工程应用中的预测结果与实际揭露结果基本保持一致，研究成果能够为工作面内地震槽波走时的精准提取提供借鉴意义，具体结论如下所述。

图8 预测煤厚与揭露煤厚对比结果

(1)选择基于Morlet小波的小波变换时频分析方法，进行槽波走时拾取时，不受固定频率和煤厚探测范围的限制，能够获得更加客观精准的槽波走时，能够有效提高煤厚及地质异常反演的精准性。

(2)由于22204工作面倾角较大，探测效果仍然比较理想，因此该方法对于大倾角工作面具有较好的适用性。

(3)由于地球物理方法往往存在多解性，而且探测结果的准确性还依赖于地球物理条件，因此应该结合钻探等其他的勘探手段，进一步提高预测结果的准确性。