CN105785440A - 一种矿井槽波双分量地震信号频散曲线提取方法 - Google Patents

Abstract

本发明涉及一种矿井槽波双分量地震信号频散曲线提取方法，本发明首先从槽波双分量地震信号中获取勒夫型槽波信号；然后对获取的勒夫型槽波信号进行广义S变换，以得到勒夫型槽波群速度与频率的关系图；将所得到群速度与频率的关系图中的振幅极大值进行连线，所得曲线即为所求的频散曲线。本发明利用线性变化的调节参数改变S变换中高斯窗口宽度，提高了槽波信号频散曲线的精度。此外，本发明在进行广义S变换之前依据采集信号中纵波在双分量检波器的能量分布准则，通过计算确定地震波传播方向及偏振角度，进行坐标旋转，以获得在垂直波传播方向水平振动的勒夫型槽波信号，再对该信号进行频散计算，获取的频散曲线频散特征更加精确。

Description

一种矿井槽波双分量地震信号频散曲线提取方法

技术领域

本发明涉及一种矿井槽波双分量地震信号频散曲线提取方法，属于地震信号处理技术领域。

背景技术

随着2010年德国DMT公司槽波地震仪在国内的推广，矿井槽波地震勘探技术迅速发展，国内多家煤矿积极开展了相关应用研究。频散是槽波地震信号最明显的特点，即槽波的传播速度是频率的函数，这是在矿井中特有岩性组合“顶板-煤层-底板”条件下产生。槽波信号频散特征中包含有围岩、煤层速度和结构信息。依据不同厚度煤层槽波频散特征，义马煤业集团地质研究所等科研单位利用槽波勘探技术在工作面煤层厚度探测方面取得了显著成果。程建远、姬广忠等学者利用正演模拟获得煤层槽波理论频散曲线，并对其特征进行了详细地分析，但对于实际槽波数据频散曲线计算方法的讨论相对较少。

槽波频散曲线计算方法是研究槽波频散特征的基础，计算方法的准确性决定了频散特征分析的可靠性。目前槽波地震信号频散曲线计算方法主要借鉴地面工程面波勘探领域中频散曲线提取方法，主要分为两类：一是基于单道信号处理，如杨真、冯涛等在《薄煤层SH型槽波频散特征及波形模式》中公开了包括多次滤波法和时频分析法，两者虽然分别从频率域和时间域进行滤波，但处理性质相同；二是基于多道信号处理，例如f-k变换法(二维傅立叶变换方法)、τ-p变换法和相移法等。f-k变换法要求空间和时间域的数据采样间隔相等，如果采样数据有坏道还会对结果产生较大的误差，τ-p变换法对高阶模态频散曲线的计算效果较好，但对基阶模态低频段频散曲线的计算效果较差。虽然多道信号频散曲线提取方法能够从多道数据中提取频散信息，获得稳定的频散曲线，但在煤矿槽波层析成像实际应用中，需要分析每道槽波信号的频散特征。因此，研究更精确地单道信号频散计算方法对槽波信号更具有实际意义。

发明内容

本发明的目的是提供一种矿井槽波双分量地震信号频散曲线提取方法，提高频散曲线的精度。

本发明为解决上述技术问题提供了一种矿井槽波双分量地震信号频散曲线提取方法，该提取方法的步骤如下：

1)从矿井采集的双分量地震信号中获取勒夫型槽波信号；

2)对获取的勒夫型槽波信号进行广义S变换，以得到勒夫型槽波群速度与频率的关系图；

3)将所得到群速度与频率的关系图中的振幅极大值进行连线，所得曲线即为所求的频散曲线。

该方法还包括对步骤1)获取的勒夫型槽波信号进行坐标转换的步骤，该坐标转换依据纵波信号在双分量检波器的能量分布准则，通过计算地震波传播方向及偏振角度实现坐标旋转，以获得在垂直波传播方向水平振动的勒夫型槽波信号。

所述的槽波信号群速度-频率广义S变换为：

$ST(\frac{x}{u},f)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}s\left(\mathrm{\τ}\right)\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-(\mathrm{\α}+\mathrm{\β}f)e\frac{f^2{(\mathrm{\τ}-\frac{x}{u})}^2}{2}}{e}^{-j2\mathrm{\π}f\mathrm{\τ}}d\mathrm{\τ}$

x为震源到检波器的距离，u为槽波传播群速度，α和β为常数，用来控制S变换的窗口函数宽度，f为频率，τ表示窗函数中心。

所述的坐标旋转角度θ为：

$\mathrm{\θ}=\frac{1}{2}arctan\left(\frac{2\underset{i=1}{\overset{n}{\Σ}}{x}_i{y}_i}{\underset{i=1}{\overset{n}{\Σ}}(x_i^2-y_i^2)}\right)$

其中x_i、y_i为原始采集双分量信号数据点。

本发明的有益效果是：本发明首先从槽波双分量地震信号中获取勒夫型槽波信号；然后对获取的勒夫型槽波信号进行广义S变换，以得到勒夫型槽波群速度与频率的关系图；将所得到群速度与频率的关系图中的振幅极大值进行连线，所得曲线即为所求的频散曲线。本发明利用线性变化的调节参数改变S变换中高斯窗口宽度，提高了槽波信号频散曲线的精度。此外，本发明在进行广义S变换之前依据采集信号中纵波在双分量检波器的能量分布准则，通过计算确定地震波传播方向及偏振角度，进行坐标旋转，以获得在垂直波传播方向水平振动的勒夫型槽波信号，再对该信号进行频散计算，获取的频散曲线频散特征更加精确。

附图说明

图1是不同厚度煤层对应理论频散曲线示意图；

图2是S变换原理图；

图3是不同参数对应调节系数示意图；

图4是不同参数对高斯函数窗口影响示意图；

图5是地震频散正演模型图；

图6-a是本发明实施例中α取0.8且β取0.01时的频散图；

图6-b是本发明实施例中α取0.8且β取0.5时的频散图；

图6-c是本发明实施例中α取10且β取0.1时的频散图；

图7是地震双分量检波器记录的原始信号图；

图8是经过旋转后得到的双分量信号示意图；

图9是纵波质点振动矢量图；

图10-a是槽波信号原始X分量频散图；

图10-b是槽波信号原始Y分量频散图；

图10-c是槽波信号旋转后水平垂直分量频散图。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步的说明。

槽波可以看作是由频率不同的简谐波在煤层中传播时相互干涉形成的一种复杂的合成振动信号，由于频散作用，合成信号振幅以一种独立的速度传播，振幅极大值沿着煤层传播的速度即为槽波传播的群速度，信号同一相位传播的速度为相速度，槽波频散分析实质上就是从采集的槽波地震数据中提取群速度或者相速度曲线。由于群速度和相速度可以相互转换，两种速度只要能求出一个，就能换算成另一种速度，转换公式如下：

$u=\frac{d\mathrm{\ω}}{d(\mathrm{\ω}/c)}=\frac{d\left(ck\right)}{dk}=c+k\frac{dc}{dk}---\left(1\right)$

式中c是相速度，u是群速度，k是圆波数，该公式表示了群速度与相速度间的关系。由于地震波能量与振幅的平方成正比，而地震波的能量主要集中在振幅极大值附近，因此群速度也就是波的能量传播速度，通过时频分析方法获得的为时频能量谱，从中可以获得槽波信号的能量分布情况，进而可以得到槽波传播的群速度。

实施例1

本实施例首先从槽波双分量地震信号中获取勒夫型槽波信号；然后对获取的勒夫型槽波信号进行广义S变换，以得到勒夫型槽波群速度与频率的关系图；将所得到群速度与频率的关系图中的振幅极大值进行连线，所得曲线即为所求的频散曲线。该方法的具体实现过程如下。

1.从槽波双分量地震信号中获取勒夫型槽波信号。

槽波分为勒夫型槽波和瑞雷型槽波两种，煤矿井下使用的水平双分量检波器接收到的槽波地震数据以勒夫型槽波为主，瑞雷型槽波不太发育，本发明针对的是勒夫型槽波频散曲线的提取。对于简单的顶班-煤层-底板三层理论模型，在煤层中心激发，由于煤层中地震波速度明显低于顶板、底板传播速度，将会在煤层中产生多次反射而干涉形成的槽波，其中勒夫型槽波形成条件相对容易满足，质点振动轨迹简单，在垂直波传播平面内水平方向上振动。

$\frac{\mathrm{\ω}d}{c_L}\sqrt{\frac{c_L^2}{{\mathrm{\υ}}_{s2}^2}-1}=\arctan \[\frac{{\mathrm{\μ}}_1}{{\mathrm{\μ}}_2}\frac{\sqrt{1-\frac{c_L^2}{{\mathrm{\υ}}_{s1}^2}}}{\sqrt{\frac{c_L^2}{{\mathrm{\υ}}_{s2}^2}-1}}\]+n\mathrm{\π}n=0,1,2,\mathrm{....}---\left(2\right)$

其中d为煤层厚度的1/2；υ_s1、υ_s2为围岩与煤层横波速度；μ₁、μ₂为围岩与煤层的剪切模量；ω为圆频率；C_L勒夫型槽波相速度，且υ_s2≤c_L≤υ_s1。n为频散阶数，当n＝0时，频散为基阶模式，薄煤层槽波以基阶模式频散为。从勒夫型槽波传播方程可知，煤层中槽波传播相速度受煤层厚度和震源激发频率影响，不同频率分量以不同速度传播，从而产生频散，随着传播距离的增大，波形长度逐渐变宽。依据对煤层实测岩石物理数据，如表1所示，利用公式(1)和(2)，建立煤层厚度分别为1m、3m、5m和7m的勒夫型槽波基阶模式群速度理论频散曲线，如图1所示。

表1

从图1中可以看出，煤层厚度对频散影响显著，煤层越厚，勒夫型槽波群速度随频率变化越快，煤层越薄，群速度随频率变化越慢，并且群速度随着频率升高逐渐降低，当频率高到一定程度时，波速将趋于稳定，接近煤层中横波速度。

2.对获取的槽波信号进行广义S变换

S变换是利用多尺度的高斯窗口进行滤波，如图2所示，在槽波信号低频区域，采用小尺度窗口进行滤波，提高了频率域分辨率；在信号高频区域，增加频率域窗口宽度以减小时间域窗口尺度，进而提高时间域分辨率，因此，S变换能够根据信号分辨率特征，自适应调整滤波窗口，进而提高槽波信号时频分析的精度。在S变换中，高斯窗口函数由频率控制，导致窗口宽度不能自由调整。而广义S变换是在S变换基础上加入调节系数得到，因此，其窗口宽度可以自由调整。

假定矿井槽波信号定义为s(t)，对应S变换的一般公式为：

$ST(t,f)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}s\left(\mathrm{\τ}\right)w(\mathrm{\τ}-t)e^{-j2\mathrm{\π}f\mathrm{\τ}}d\mathrm{\τ}---\left(3\right)$

式中，ST(t,f)表示信号s(t)的S变换，t和f分别代表时间和频率，τ表示窗函数中心，w(t)为窗函数，定义为随着频率变化的高斯函数，公式为：

$w\left(t\right)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{f^2{i}^2}{2}}---\left(4\right)$

此时S变换公式为：

$ST(t,f)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}s\left(\mathrm{\τ}\right)\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{f^2{(\mathrm{\τ}-t)}^2}{2}}{e}^{-j2\mathrm{\π}f\mathrm{\τ}}d\mathrm{\τ}---\left(5\right)$

S变换中高斯窗口函数受频率控制，窗口宽度不能自由调整，若要实现窗口的调整，需在窗口函数中加入调节系数，而窗口变换类型较多，可统称为广义S变换，本发明基于广义S变换的原理，通过增加调节系数来实现对窗口宽度的控制，其采用的窗口函数为：

$w\left(t\right)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\mathrm{\λ}\left(f\right)\frac{f^2{t}^2}{2}}---\left(6\right)$

λ(f)＝α+βf(7)

其中λ(f)为调节系数函数，为频率的线性函数，参数α和β为常数，用来控制S变换的窗口函数宽度。随着α和β的增加，高斯窗口调节系数呈线性增加，如图3所示。更改调节系数函数中参数α和β的数值，可自由调整高斯窗口的宽度，其中α用于调整高斯窗口宽度的范围，数值越大，窗口越窄；β用于调整高斯窗口边界的斜率，数值越高，窗口边界越陡，如图4所示。

槽波信号进行S变换后得到时间-频率能量分布图，因频散分析通常分析速度与频率的关系，可采用公式(8)将时间转换为速度，根据前述槽波速度概念，该速度为群速度。

u(f)＝x/t(8)

公式中x为震源到检波器的距离，将公式(6)、(7)和(8)代入公式(3)，可到槽波信号速度-频率S变换公式：

$ST(\frac{x}{u},f)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}s\left(\mathrm{\τ}\right)\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-(\mathrm{\α}+\mathrm{\β}f)\frac{f^2{(\mathrm{\τ}-\frac{x}{u})}^2}{2}}{e}^{-j2\mathrm{\π}f\mathrm{\τ}}{d}_{\mathrm{\τ}}---\left(9\right)$

其中u为槽波传播群速度，其它变量定义同前。将单道槽波信号所有数据点利用上式进行变换后，即得到槽波群速度与频率的关系图。需要注意将时间转换为速度后，原来等间隔的时间采样点，在速度域中数据点不再是等间隔分布，因此需要进行重新插值网格化后再成图，

3.将图中振幅极大值进行连线，所得结果即为对应的频散曲线。

下面对本发明的效果进行分析验证，本实施例根据图2中煤层厚度3m对应的频散曲线数据建立槽波正演模型，设定已知频散曲线为v(f)。假定炮点激发地震子波用以下函数表示。

w(t,f_i)＝ae^-btsin(2πf_it)(10)

其中a和b为常数，分别代表振幅系数和阻尼系数，用来控制地震子波的振幅大小和衰减快慢，速度v_i＝v(f_i)是频率的函数，随着频率的改变而变化。当给定炮点到检波器距离x，即可计算不同频率f_i时，地震子波的延迟时间。

${\mathrm{\Δt}}_i=\frac{x}{v_i}---\left(11\right)$

将所有频率的地震子波进行求和计算即可合成频散地震波形，合成的地震记录为：

$s(t,x)=\underset{i=1}{\overset{n}{\Σ}}w(t+\frac{x}{v\left(f\right)})=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{ae}}^{-b(t+\frac{x}{v\left(f_i\right)})}sin\left(2{\mathrm{\πf}}_i\right(t+\frac{x}{v\left(f_i\right)}\left)\right)---\left(12\right)$

其中n为已知频散曲线中所有数据点个数，x为炮点到检波器的距离，本实施例利用频散地震记录合成公式得到的正演模型如图5所示，其中炮点在左边激发，地震道间距为10m。以图5中的第15道数据为例进行频散计算，调节系数可采用不同数值的α和β，当α取0.8且β取0.01时，其对应的频散计算结果如图6-a所示，该图中能量集中，极大值连线与煤厚3m对应理论频散曲线一致，当α和β取其他数值，得到的频散计算结果如图6-b和图6-c所示，能量发散，不集中。因此利用调节参数可改变高斯窗口宽度的频散，从而能够控制信号频散的计算精度。

实施例二

由于原始采集信号中质点振动方位与波传播方向关系不明确，采集信号中包含多种类型地震波，导致频散图中频散规律不明显，能量极大值连线发散，不连续。为此，本发明在实施例一的基础上，进行了进一步改进，在对槽波信号进行广义S变换之前，将获取的双分量槽波信号分别转换到水平平行波传播方向和水平垂直波传播方向，从而获得质点在垂直波传播方向水平振动的勒夫型槽波信号，再进行频散计算，以得到勒夫型槽波频散曲线。下面就槽波信号进行坐标转换的过程进行详细说明。

当前槽波地震井下仪器多采用水平双分量检波器(X分量和Y分量)，其中X分量平行煤壁，Y分量垂直煤壁放置，由于槽波从震源传到各检波器方位角不同，采集点各检波器记录的质点振动轨迹与波传播方向关系不明确，导致频散分析时，难以确定是对X分量处理，还是对Y分量进行处理。为统一双分量检波器接收信号的传播方向，可利用极化滤波方法进行波场分离，该方法精度高但实现复杂。经过旋转后，将槽波水平传播方向定为X’，槽波水平传播垂直方向为Y’，坐标旋转公式为：

$\left\{\begin{array}{c}{x}^{\′}=x\cos \mathrm{\θ}+y\sin \mathrm{\θ}\\ y^{\′}=y\cos \mathrm{\θ}-x\sin \mathrm{\θ}\end{array}\right.---\left(13\right)$

式中θ即为坐标旋转的角度，x、y为原始采集双分量信号，x′分量为槽波水平传播方向接收到的信号，y′分量为槽波水平传播垂直波传播方向接收到的信号。

为获得坐标旋转角度θ，利用纵波传播能量准则法进行计算。煤层中检波器接收到地震波信号，除了槽波以外，最先到达的是纵波。图7为某矿实际采集的双分量槽波信号，该检波器接收到的X分量和Y分量中，最先接收到振动(50-70ms)记录即为纵波信号，在时间185-225ms为槽波信号中最明显的埃里相。由于质点振动方向与波传播方向关系不明确，水平垂直波传播方向振动的勒夫型槽波在X和Y分量中都存在。

由纵波传播规律可知，纵波质点振动方向与波传播方向一致，也就是说纵波在波传播方向能量最大，而在波传播垂直方向能量最小。而地震波能量与振幅平方呈正比，因此，可以定义纵波在波传播方向能量为P_x(θ)，纵波在波传播垂直方向上能量为P_y(θ)，公式如下：

$\left\{\begin{array}{c}{P}_x\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{n}{\Σ}}{\left(x_i^{\′}\right)}^2=\underset{i=1}{\overset{n}{\Σ}}{(x_i\cos \mathrm{\θ}+y_i\sin \mathrm{\θ})}^2\\ P_y\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{n}{\Σ}}{\left(y_i^{\′}\right)}^2=\underset{i=1}{\overset{n}{\Σ}}{(y_i\cos \mathrm{\θ}-x_i\sin \mathrm{\θ})}^2\end{array}\right.---\left(14\right)$

式中n为数据点个数，当P_x(θ)取最大值，同时P_y(θ)为最小值时的方位角即为地震波需要进行坐标旋转的角度θ，设定目标函数M(θ)：

$M\left(\mathrm{\θ}\right)=\frac{\underset{i=1}{\overset{n}{\Σ}}{(y_{i}cos\mathrm{\θ}-x_i\sin \mathrm{\θ})}^2}{\underset{i=1}{\overset{n}{\Σ}}{(x_{i}cos\mathrm{\θ}+y_i\sin \mathrm{\θ})}^2}---\left(15\right)$

当目标函数M(θ)取极小值时，对应θ为所求角度。其取极小值的必要条件为经过计算可以求出旋转角度θ：

$\mathrm{\θ}=\frac{1}{2}arctan\left(\frac{2\underset{i=1}{\overset{n}{\Σ}}{x}_i{y}_i}{\underset{i=1}{\overset{n}{\Σ}}(x_i^2-y_i^2)}\right)---\left(16\right)$

根据公式(16)利用纵波数据进行计算，获得该检波器需要旋转的角度θ为105°。旋转后的信号如图8所示，从图中可以看出，在波传播水平平行分量上，因纵波振动方向与传播方向一致，纵波能量最强，而在水平垂直分量上，纵波能量最弱。同样，质点振动方向与波传播方向垂直的勒夫型槽波，在旋转后的水平垂直分量上能量最强，特征最明显。

将旋转前后X分量和Y分量的纵波信号数据点，在直角坐标系中按照采样点数值大小和符号标出，并连成曲线(如图9所示)。在理想无干扰信号情况下，纵波振动方向与波传播方向一致，运动轨迹应该为一条直线，但由于随机噪音的干扰，表现为极扁的椭圆。这个椭圆的长轴方向与x轴之间的夹角，就是纵波质点振动方向与x轴之间的夹角θ。可以看出，旋转前，椭圆长轴与X轴角度约为105°，经过坐标旋转后，纵波质点振动方向与X轴方向(波传播方向)一致。经过坐标旋转后，根据勒夫型槽波质点振动方向垂直波传播方向的特征，可认为旋转后的Y’分量记录的信号为勒夫型槽波信号。

根据上述坐标旋转后，根据勒夫型槽波质点振动方向垂直波传播方向特征，旋转后的分量记录的信号Y’即为勒夫型槽波信号，对比坐标旋转前后结果，同时对原始X分量和Y分量信号分别进行频散计算，结果如图10-a和图10-b所示，频散图中频散规律不明显，能量极大值连线发散，不连续，经过坐标旋转后，在水平垂直分量上主要为勒夫型槽波，频散特征明显，规律清晰，如图10-c所示，与理论频散曲线(如图1所示)特征一致。

本发明利用线性变化的调节参数改变S变换中高斯窗口宽度，提高了槽波信号频散曲线计算的精度，依据纵波信号在双分量检波器的能量分布准则，计算地震波传播方向及偏振角度，进行坐标旋转，以获得在水平垂直分量上振动的勒夫型槽波信号。经过本发明提取的勒夫型槽波地震信号频散特征明显，槽波速度受煤层厚度及频率影响大，当频率增高时，槽波速度快速减小。

Claims (4)

1.一种矿井槽波双分量地震信号频散曲线提取方法，其特征在于，该提取方法的步骤如下：

1)从矿井采集的双分量地震信号中获取勒夫型槽波信号；

2)对获取的勒夫型槽波信号进行广义S变换，以得到勒夫型槽波群速度与频率的关系图；

3)将所得到群速度与频率的关系图中的振幅极大值进行连线，所得曲线即为所求的频散曲线。

2.根据权利要求1所述的矿井槽波双分量地震信号频散曲线提取方法，其特征在于，该方法还包括对步骤1)获取的勒夫型槽波信号进行坐标转换的步骤，该坐标转换依据纵波信号在双分量检波器的能量分布准则，通过计算地震波传播方向及偏振角度实现坐标旋转，以获得在垂直波传播方向水平振动的勒夫型槽波信号。

3.根据权利要求1或2所述的矿井槽波双分量地震信号频散曲线提取方法，其特征在于，所述的槽波信号群速度-频率广义S变换为：

$ST(\frac{x}{u},f)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}s\left(\mathrm{\τ}\right)\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-(\mathrm{\α}+\mathrm{\β}f)\frac{f^2{(\mathrm{\τ}-\frac{x}{u})}^2}{2}}{e}^{-j2\mathrm{\π}f\mathrm{\τ}}d\mathrm{\τ}$

x为震源到检波器的距离，u为槽波传播群速度，α和β为常数，用来控制S变换的窗口函数宽度，f为频率，τ表示窗函数中心。

4.根据权利要求2所述的矿井槽波双分量地震信号频散曲线提取方法，其特征在于，所述的坐标旋转角度θ为：

$\mathrm{\θ}=\frac{1}{2}arctan\left(\frac{2\underset{i=1}{\overset{n}{\Σ}}{x}_i{y}_i}{\underset{i=1}{\overset{n}{\Σ}}(x_i^2-y_i^2)}\right)$

其中x_i、y_i为原始采集双分量信号数据点。