CN115437006A - 一种煤岩界面折射波的多阶折射时间求取方法 - Google Patents

Abstract

本发明涉及一种煤岩界面折射波的多阶折射时间求取方法。当震源波长大于煤厚时，相邻阶的折射波相互混叠，常规反褶积方法难以对其稀疏化。本发明利用各阶折射波相似性高、阶数已知的特点，借鉴正交匹配追踪算法，能够从严重混叠的折射波信号中准确提取各阶折射时间。本发明考虑到了各阶折射波之间的相位差异，通过构造包含相位旋转角的字典矩阵，使各阶折射波能准确匹配对应子波，获得准确的脉冲信号，提高所得时间点的准确性。

Description

一种煤岩界面折射波的多阶折射时间求取方法

技术领域

本发明涉及一种多阶折射时间求取方法，属于煤层地质探测技术领域，具体是涉及一种煤岩界面折射波的多阶折射时间求取方法。

背景技术

煤层厚度是采煤过程的重要参数。在回采前准确探测出工作面内煤层厚度的分布情况，是采煤机动态调整位置和采高，从而实现煤矿工作面智能开采的前提和基础，对提高回采效率、减少资源浪费以及保障煤矿安全生产均具有重要意义。

震源在煤层内激发后，产生的地震波将在煤层与顶板和底板的两个分界面之间产生多次反射，其中仅有以临界角入射的地震波能够产生沿界面传播的折射波；因此当煤层的厚度稳定时，接收到的地震信号的折射波部分将具有明显的周期性，即在围岩速度对应的时间点上出现初至折射波之后，每隔固定的时间周期出现续至折射波；当工作面内存在薄煤带等煤厚异常区域时，续至折射波出现的时间点会发生移动，偏离正常的周期，因此求出各阶折射波出现的时间，就可以由此计算出异常煤层的厚度。

专利“一种基于折射波周期振幅衰减的薄煤带探测方法及系统”(公开号：202111434047.5)利用在煤层与围岩界面上传播的折射波的周期特性，探测采煤工作面内薄煤带等煤厚异常区的位置和边界。此专利方法仅能圈定出薄煤带的范围，不能准确估计出薄煤带内具体的煤厚及其分布情况。这是由于此专利中对折射波信号进行反褶积的方法较为简单，反褶积后的折射波信号不是稀疏的，从中无法准确求出多阶折射发生的具体时间，因此不能确定各个折射点对应的煤层厚度。为此，需要针对折射波特性采用更有效的算法对折射波信号稀疏化，并从中求取多阶折射时间，为准确估计煤厚分布创造条件。

发明内容

以下给出一个或多个方面的简要概述以提供对这些方面的基本理解。此概述不是所有构想到的方面的详尽综览，并且既非旨在指认出所有方面的关键性或决定性要素亦非试图界定任何或所有方面的范围。其唯一的目的是要以简化形式给出一个或多个方面的一些概念以为稍后给出的更加详细的描述之序。

本发明主要的目的是解决现有技术中所存在的技术问题，提供了一种煤岩界面折射波的多阶折射时间求取方法，该方法简单易用，适应于现有槽波探测工程的观测系统与施工条件，能够将基于折射波的煤厚探测方法的适用性条件由原来的巨厚煤层推广到厚煤层与中厚煤层中。

为解决上述问题，本发明的方案是：

一种煤岩界面折射波的多阶折射时间求取方法，包括：

计算折射波的折射阶数K，并将子波w_n(t)合成为字典矩阵W；

按以下迭代步骤计算得到稀疏信号r(t)：

迭代步骤1，初始化，令k＝1，A₀为空集，u₁＝s；

迭代步骤2，寻找索引i_k，使

即将W的所有列向量与u_k求内积，找到使内积最大的那一列对应的列序号；

迭代步骤3，基于下式求取迭代后的A_k：

A_k＝A_k-1∪a_ik；

迭代步骤4，求s＝A_kr_k的最小二乘解，即：

迭代步骤5，更新残差向量u_k+1＝s-A_kr_k；

迭代步骤6，令k＝k+1,如果k≤K则返回至迭代步骤2，否则将当前向量r_k作为稀疏信号r(t)输出；

迭代步骤7，拾取稀疏信号r(t)的K个脉冲峰值对应的时间，得到折射波的多阶折射时间；

其中，各迭代步骤中，k为当前迭代次数，r_k表示第k次迭代的稀疏向量，u_k表示第k次迭代的残差向量，i_k表示第k次迭代找到的索引，所述索引为W矩阵中对应的列序号，a_ik表示矩阵W中以i_k为序号的列向量，A_k为一个由所有挑选出的列向量a_ik所组成的集合矩阵。

优选的，上述的一种煤岩界面折射波的多阶折射时间求取方法，基于下式计算折射阶数K：

由已知的激发点位置(x_s,y_s)和接收点位置(x_r,y_r)，计算折射波的折射阶数K：

式中，地震波速度v_r、煤层的地震波速度v_c，T为折射波周期。

优选的，上述的一种煤岩界面折射波的多阶折射时间求取方法，基于下式计算折射波周期：

其中，围岩的地震波速度v_r、煤层的地震波速度v_c。

优选的，上述的一种煤岩界面折射波的多阶折射时间求取方法，基于下式产生子波：

其中t∈[0,NΔt]，周期相位旋转角Φ，折射波周期T，t表示时间，t∈[0,NΔt]，N为采集的折射波信号的长度，Δt是采样间隔。

优选的，上述的一种煤岩界面折射波的多阶折射时间求取方法，基于下式产生字典矩阵W：

其中符号∪表示集合并运算。

优选的，上述的一种煤岩界面折射波的多阶折射时间求取方法，基于下式求取使

达到最小时的

作为周期相位旋转角Φ：

其中，

主频为f_c的雷克子波w₀(τ)，τ表示时间，τ∈[-T,T]，T为折射波周期；

H(·)为希尔伯特变换，τ∈[-T,T]。

优选的，上述的一种煤岩界面折射波的多阶折射时间求取方法，基于下式构建主频为f_c的雷克子波w₀(τ)：

其中，τ表示时间，τ∈[-T,T]，T为折射波周期。

因此，相对于现有技术，本发明的优点是：

(1)本发明可实现震源波长大于煤厚条件下的折射波处理。当震源波长大于煤厚时，相邻阶的折射波相互混叠，常规反褶积方法难以对其稀疏化。本发明利用各阶折射波相似性高、阶数已知的特点，借鉴正交匹配追踪算法，能够从严重混叠的折射波信号中准确提取各阶折射时间。

(2)本发明考虑到了各阶折射波之间的相位差异，通过构造包含相位旋转角的字典矩阵，使各阶折射波能准确匹配对应子波，获得准确的脉冲信号，提高所得时间点的准确性。

附图说明

并入本文并形成说明书的一部分的附图例示了本发明的实施例，并且附图与说明书一起进一步用于解释本发明的原理以及使得所属领域技术人员能够制作和使用本公开。

图1为煤岩界面折射波的多阶折射时间求取方法的总流程图。

图2为工作面模型的示意图，其中，(a)为模型及观测系统平面图，(b)为模型剖面图；

图3为工作面模型不同主频震源的三维弹性波正演模拟结果，其中，(a)为200Hz主频震源产生的折射波信号，(b)为400Hz主频震源产生的折射波信号，(c)为R15接收到的两种主频的折射波信号的对比图；

图4为震源附近折射波信号及其频谱，其中，(a)为折射波信号，(b)为频谱；

图5为以200Hz为主频构建雷克子波；

图6为子波与折射波信号的误差函数

图7为R15道折射波信号s(t)与用Φ对子波w₀进行周期相位旋转得到的合成信号的对比图；

图8为字典矩阵W对应的图像；

图9为s(t)与得到的稀疏信号r(t)的示意图；

图10为对所有数据逐道处理，得到的稀疏信号；

将参照附图描述本发明的实施例。

具体实施方式

实施例

当震源波长大于煤厚时，相邻阶的折射波相互混叠，常规反褶积方法难以对其稀疏化，因而本实施例的方法利用各阶折射波相似性高、阶数已知的特点，借鉴正交匹配追踪算法，并且考虑到了各阶折射波之间的相位差异，通过构造包含相位旋转角的字典矩阵，使各阶折射波能准确匹配对应子波，获得准确的脉冲信号，提高所得时间点的准确性。

下面结合附图对本实实施例进行详细说明。如图1所示，本实施例的煤岩界面折射波的多阶折射时间求取方法包括：

步骤1，折射波周期计算步骤。由已知的围岩的地震波速度v_r、煤层的地震波速度v_c，煤层平均厚度h，计算出折射波周期T：

步骤2，由已知的激发点位置(x_s,y_s)和接收点位置(x_r,y_r)，计算折射波的折射阶数K：

步骤3，估计震源主频f_c。

选取激发点附近的折射波信号s₀(t)，其中t表示时间，t∈[0,NΔt]，N为采集的折射波信号的长度，Δt是采样间隔。用傅里叶变换求出其振幅谱，寻找到最大振幅对应的频率，就是震源主频f_c。

步骤4，构建主频为f_c的雷克子波w₀(τ)。

其中τ表示时间，τ∈[-T,T]，T为折射波周期。

步骤5，估计出周期相位旋转角Φ，即折射波信号每经过一个周期，其相位发生旋转的角度，具体包括：

步骤5.1，求出折射波的初至时间t₀:

步骤5.2，取

计算子波与折射波信号误差函数

其中：

H(·)为希尔伯特变换，τ∈[-T,T]。

步骤5.3，寻找使

达到最小时的

则

为估计出的周期相位旋转角Φ。

步骤6，产生一个N×N的字典矩阵W，N为信号长度，包括以下步骤：

步骤6.1令n＝0，w为空集；

步骤6.2，产生子波w_n(t)：

其中t∈[0,NΔt]。

步骤6.3，将子波w_n(t)记为N维向量w_n；

步骤6.4，令n＝n+1,如果n<N则返回至第6.2步；

步骤6.5，将所有w_n合并成一个矩阵W：

其中符号∪表示集合并运算。

步骤7，将折射阶数K作为输入参数，采用正交匹配追踪算法用字典矩阵W对折射信号s(t)进行稀疏反褶积，得到稀疏信号r(t)。具体实施为：

将长度为N的折射波信号s(t)记为N维向量s，设k为当前迭代次数，r_k表示第k次迭代的稀疏向量，u_k表示第k次迭代的残差向量，i_k表示第k次迭代找到的索引(W矩阵中对应的列序号)，a_ik表示矩阵W中以i_k为序号的列向量，A_k为一个由所有挑选出的列向量a_ik所组成的集合矩阵。符号∪表示集合并运算，符号<·，·>表示向量内积。按下述步骤迭代：

步骤7.1，初始化，令k＝1，A₀为空集，u₁＝s；

步骤7.2，寻找索引i_k，使

即将W的所有列向量与u_k求内积，找到使内积最大的那一列对应的列序号；

步骤7.3，A_k＝A_k-1∪a_ik；

步骤7.4，求s＝A_kr_k的最小二乘解：

步骤7.5，更新残差向量u_k+1＝s-A_kr_k；

步骤7.6，令k＝k+1,如果k≤K则返回至第(2)步，否则将当前向量r_k作为稀疏信号r(t)输出；

步骤8，稀疏信号r(t)具有K个脉冲，拾取这K个脉冲峰值对应的时间，得到折射波的多阶折射时间。

下面以正演模型数据为例说明本发明效果：

模型由顶板、煤层、底板三层组成，顶板与底板岩性相同。围岩纵波速度4000m/s，围岩横波速度2300m/s，密度2.56g/cm³；煤层纵波速度2000m/s，煤层横波速度1050m/s，密度1.4g/cm³。煤层厚度为10m，煤层中包含两条平行巷道，两巷间距200m。在煤层中存在一处薄煤区，大小为100m×50m，煤厚为5m。图2(a)为模型及观测系统平面图，图2(b)为模型剖面图。

激发点和接收点分别布置在两条巷道内，位置如图2(a)中所示。采用不同主频的雷克子波作为震源进行三维弹性波数值模拟，200Hz和400Hz主频震源产生的折射波信号分别如图3(a)和3(b)所示。图3(c)为R15接收到的两种频率的折射波信号的对比图，可以看出400Hz信号中各阶折射波彼此分隔，能直观的分辨出各阶折射波并从中找到最大振幅对应的各阶折射时间；而200Hz折射波则相互混叠，不能直观分辨出各阶折射波波形，因此无法找到各阶折射时间。

现依据本发明步骤对R15接收的折射波信号s(t)进行处理，求出各阶折射时间。

执行步骤1，由模型参数v_r＝4000m/s、v_c＝2000m/s，h＝10m，求出折射波周期T＝86.6ms。

执行步骤2，由激发点位置(x_s,y_s)＝(400,55)和接收点R15位置(x_r,y_r)＝(140,250)，求得K＝5。

执行步骤3，震源附近折射波信号s₀(t)，及其频谱分别如图4(a)和4(b)所示，得到主频f_c为200Hz。

执行步骤4，以200Hz为主频构建雷克子波w₀，如图5所示。

执行步骤5(1)，计算出折射波初至时间t₀＝71ms。

执行步骤5(2)，子波与折射波信号误差函数

如图6所示。

执行步骤5(3)，由图5估计出周期相位旋转角Φ＝1.24。图7为R15道折射波信号s(t)与用Φ对子波w₀进行周期相位旋转得到的合成信号的对比图。

执行步骤6，合成字典矩阵W。由于截取的折射波信号长度为1600个采样点，因此产生的W矩阵尺寸为1600×1600。将W中各元素作为图像中各点，则W矩阵可输出为一个1600×1600的图像，如图8所示，横坐标为矩阵的行，纵坐标为矩阵的列，图像中的值为矩阵中元素对应的颜色(例如灰度值)。

执行步骤7，对折射波信号s(t)进行如下处理：

执行步骤7(1)，完成初始化，令k＝1，将A₀设为空集，将信号s(t)中数据设为u₁；

执行步骤7(2)，将W中各列与u₁求内积，找到使内积最大的列为第1108列，则i_k＝1108；

执行步骤7(3)，将A₀与W中的第1108列合并成A₁；

执行步骤7(4)，计算出r₁；

执行步骤7(5)，计算出u₁；

执行步骤7(7)，由于k≤K，则令k＝k+1，返回步骤7(2)，直到k>K；s(t)与得到的稀疏信号r(t)如图9所示。

执行步骤8，从图9中的稀疏信号中，提取出各个最大值点对应的时间，分别为(76,849,938,102.4,1108,1174)ms。

可以看出虽然各阶折射波之间发生了混叠，但采用本发明获得的各阶折射时间相差基本为86ms，即一个折射波周期。因此，本发明能够准确获得的各阶折射时间。对所有数据逐道按上述步骤处理，得到的稀疏信号如图10所示，对比图3a，可以很明显看出右侧圈出区域内各阶折射时间存在异常。这种异常是由模型中的薄煤区引起的。

总体而言，本发明利用煤岩界面上传播的折射波的周期性对其进行稀疏反褶积，并从稀疏信号中获取多阶折射时间。准确的多阶折射时间是圈定煤厚异常区，并估计异常区内煤层厚度的先决条件，对保障煤矿工作面的智能化开采，提高回采效率和采出率，减少资源浪费具有重要意义。

本实施例中，尽管为使解释简单化将上述方法图示并描述为一系列动作，但是应理解并领会，这些方法不受动作的次序所限，因为根据一个或多个实施例，一些动作可按不同次序发生和/或与来自本文中图示和描述或本文中未图示和描述但本领域技术人员可以理解的其他动作并发地发生。

注意到，说明书中对“一个实施例”、“实施例”、“示例实施例”、“一些实施例”等的引用指示所描述的实施例可以包括特定特征、结构或特性，但是每个实施例可以不必包括所述特定特征、结构或特性。而且，这样的短语不必指代同一实施例。此外，当结合实施例描述特定特征、结构或特性时，无论是否明确描述，结合其他实施例来实现这样的特征、结构或特性将在所属领域的技术人员的知识范围内。

提供对本公开的先前描述是为使得本领域任何技术人员皆能够制作或使用本公开。对本公开的各种修改对本领域技术人员来说都将是显而易见的，且本文中所定义的普适原理可被应用到其他变体而不会脱离本公开的精神或范围。由此，本公开并非旨在被限定于本文中所描述的示例和设计，而是应被授予与本文中所公开的原理和新颖性特征相一致的最广范围。

Claims (7)

1.一种煤岩界面折射波的多阶折射时间求取方法，其特征在于，包括：

计算折射波的折射阶数K，并将子波w_n(t)合成为字典矩阵W；

按以下迭代步骤计算得到稀疏信号r(t)：

迭代步骤1，初始化，令k＝1，A₀为空集，u₁＝s；

迭代步骤2，寻找索引i_k，使

即将W的所有列向量与u_k求内积，找到使内积最大的那一列对应的列序号；

迭代步骤3，基于下式求取迭代后的A_k：

A_k＝A_k-1∪a_ik；

迭代步骤4，求s＝A_kr_k的最小二乘解，即：

迭代步骤5，更新残差向量u_k+1＝s-A_kr_k；

迭代步骤6，令k＝k+1,如果k≤K则返回至迭代步骤2，否则将当前向量r_k作为稀疏信号r(t)输出；

迭代步骤7，拾取稀疏信号r(t)的K个脉冲峰值对应的时间，得到折射波的多阶折射时间；

其中，各迭代步骤中，k为当前迭代次数，r_k表示第k次迭代的稀疏向量，u_k表示第k次迭代的残差向量，i_k表示第k次迭代找到的索引，所述索引为W矩阵中对应的列序号，a_ik表示矩阵W中以i_k为序号的列向量，A_k为一个由所有挑选出的列向量a_ik所组成的集合矩阵。

2.根据权利要求1所述的一种煤岩界面折射波的多阶折射时间求取方法，其特征在于，基于下式计算折射阶数K：

由已知的激发点位置(x_s,y_s)和接收点位置(x_r,y_r)，计算折射波的折射阶数K：

式中，地震波速度v_r、煤层的地震波速度v_c，T为折射波周期。

3.根据权利要求2所述的一种煤岩界面折射波的多阶折射时间求取方法，其特征在于，基于下式计算折射波周期：

其中，围岩的地震波速度v_r、煤层的地震波速度v_c。

4.根据权利要求1所述的一种煤岩界面折射波的多阶折射时间求取方法，其特征在于，基于下式产生子波：

其中t∈[0,NΔt]，周期相位旋转角Φ，折射波周期T，t表示时间，

t∈[0,NΔt]，N为采集的折射波信号的长度，Δt是采样间隔。

5.根据权利要求4所述的一种煤岩界面折射波的多阶折射时间求取方法，其特征在于，基于下式产生字典矩阵W：

其中符号∪表示集合并运算。

6.根据权利要求4所述的一种煤岩界面折射波的多阶折射时间求取方法，其特征在于，基于下式求取使

达到最小时的

作为周期相位旋转角Φ：

其中，

主频为f_c的雷克子波w₀(τ)，τ表示时间，τ∈[-T,T]，T为折射波周期；

H(·)为希尔伯特变换，τ∈[-T,T]。

7.根据权利要求6述的一种煤岩界面折射波的多阶折射时间求取方法，其特征在于，基于下式构建主频为f_c的雷克子波w₀(τ)：

其中，τ表示时间，τ∈[-T,T]，T为折射波周期。

Images