CN105092023A - 基于白噪声统计特征的锚杆振动信号矫正方法 - Google Patents

Abstract

本发明公开了一种基于白噪声统计特征的锚杆振动信号矫正方法，对布设于被测试锚杆上的锚杆振动检测单元所检测振动信号进行矫正，包括步骤：一、振动信号经验模态分解：对锚杆振动检测单元所检测振动信号进行经验模态分解，获得n个本征模态函数分量和一个趋势项；二、本征模态函数分量有效性识别：调用本征模态函数分量有效性识别模块且根据预先建立的白噪声近似熵分布区间模板，对分解出的n个本征模态函数分量分别进行有效性识别；三、信号重构：利用识别出的振动信号经验模态分解后的所有有效本征模态函数分量重构振动信号。本发明方法步骤简单、设计合理且实现方便、使用效果好，能对锚杆振动信号进行有效矫正，获得锚杆的真实振动状态。

Description

基于白噪声统计特征的锚杆振动信号矫正方法

技术领域

本发明涉及一种信号矫正方法，尤其是涉及一种基于白噪声统计特征的锚杆振动信号矫正方法。

背景技术

根据煤矿企业安全生产规程(即《煤矿安全规程》)规定，生产前要采用钻探、物探及巷探等技术手段，探测综采工作面(简称综采面)安全距离内的地质构造、水文地质情况、煤、岩层位及其它相关地质情况，防止瓦斯突出和突水事故的发生，确保安全生产。但是，随着近十年来我国煤矿开采深度的加大和开采规模的提高，导致开采条件更趋复杂；而依据传统的超前探测方法无法动态地、高分辨地跟踪综采面近场煤岩裂隙场变化。其中，综采面近场煤岩指的是综采工作面(简称综采面)安全距离内的煤岩。因此，迫切需要找到能够实时跟踪综采面近场采动煤岩的结构破坏与裂隙演化的动态可视化方法，从而对即将发生事故的概率进行科学评价和预测。另外，目前常用的综采面近场煤岩探测方法有：直流电法探测、音频电透法、瞬变电磁探测。但是上述现有的探测方法必须在综采面停产/准停产的条件下进行专门探测施工。对处于工作状态下的综采面近场煤岩结构破坏与裂隙演化无法进行动态测量，从而无法为煤矿安全的动态预测预警提供有效的帮助。

如能在采煤过程中对综采面近场煤岩的振动状态进行探测，从而掌握处于工作状态下的综采面近场煤岩结构破坏与裂隙演化规律，具有重大的实际意义。采煤过程中对综采面近场煤岩的振动状态进行探测时，将振动检测单元布设在综采面两侧回采巷道内的锚杆上，将测试得出的锚杆振动信号作为综采面近场煤岩该位置处的振动信号，具有安装布设方便、投入成本低、检测效果好等优点。为提高综采面近场煤岩振动状态的探测精度，需对检测出的锚杆振动信号进行矫正，主要目的在于去除锚杆振动信号中含有的噪声信号。经验模态分解(EMD)是一种常用的信号去噪方法，EMD是将信号自适应地分解为一系列IMF分量(即本征模态函数分量)，但分解出的部分IMF分量并没有从物理上解释信号的特征，该分量被称之为虚假分量(falseintrinsicmodefunctions，FIMF)。由于EMD利用信号的局部极值信息筛选IMF分量是导致虚假模态产生的一个原因，另外由于所检测的振动信号中混入噪声干扰，使得真实信号中出现虚假极值点，同样导致了一些高频噪声和低频分量等虚假模态。因而，需设计一种方法步骤简单、设计合理且实现方便、使用效果好的基于白噪声统计特征的锚杆振动信号矫正方法，能对锚杆振动信号进行有效矫正，获得锚杆的真实振动状态。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于白噪声统计特征的锚杆振动信号矫正方法，其方法步骤简单、设计合理且实现方便、使用效果好，能对锚杆振动信号进行有效矫正，获得锚杆的真实振动状态。

为解决上述技术问题，本发明采用的技术方案是：一种基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：采用数据处理设备对布设于被测试锚杆上的锚杆振动检测单元所检测振动信号x(t)进行矫正，包括以下步骤：

步骤一、振动信号经验模态分解：调用经验模态分解模块，对锚杆振动检测单元所检测振动信号x(t)进行经验模态分解，获得n个本征模态函数分量和一个趋势项；n个本征模态函数分量和一个趋势项均为离散函数；

其中，振动信号x(t)为锚杆的振动位移或振动加速度随时间变化的函数；振动信号x(t)的经验模态分解结果为：式(1)中c_i(t)为振动信号x(t)经验模态分解后获得的第i个本征模态函数分量，i为正整数且其本征模态函数分量的编号，i＝1、2、…、n；r_n(t)为振动信号x(t)经验模态分解后获得的趋势项；

步骤二、本征模态函数分量有效性识别：调用本征模态函数分量有效性识别模块，且根据预先建立的白噪声近似熵分布区间模板，对步骤一中分解出的振动信号x(t)的n个本征模态函数分量分别进行有效性识别；n个本征模态函数分量的有效性识别方法均相同；

对所述白噪声近似熵分布区间模板进行建立时，采用所述数据处理设备进行建立，过程如下：

步骤Ⅰ、样本获取：产生多个随机白噪声信号作为样本，多个所述白噪声信号的数据长度均相同；

其中，产生的第j个随机白噪声信号，记作ζ_j(t)；j为随机白噪声信号的编号且其为正整数，j＝1、2、…、N；N为本步骤中所产生随机白噪声信号的总数量；

步骤Ⅱ、经验模态分解：调用所述经验模态分解模块，对步骤Ⅰ中各白噪声信号分别进行经验模态分解，获得各白噪声信号的n个本征模态函数分量和一个趋势项；n个本征模态函数分量和一个趋势项均为离散函数；

其中，ζ_j(t)的经验模态分解结果为：式(1)中c_i,j(t)为ζ_j(t)经验模态分解后获得的第i个本征模态函数分量，r_n,j(t)为ζ_j(t)经验模态分解后获得的趋势项；

步骤Ⅲ、近似熵与平均周期自然对数计算：调用近似熵计算模块，对步骤Ⅱ中各白噪声信号的n个本征模态函数分量的近似熵分别进行计算；同时，调用平均周期计算模块，对各白噪声信号的n个本征模态函数分量的平均周期分别进行计算，并求解出各白噪声信号的n个本征模态函数分量的平均周期对数；其中，平均周期对数为平均周期的自然对数；

步骤Ⅳ、二维联合概率密度函数确定：根据步骤Ⅲ中计算得出的各白噪声信号的n个本征模态函数分量的近似熵和平均周期对数，调用二维联合概率密度函数确定模块，获取白噪声信号的n个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数；

其中，根据步骤Ⅲ中计算得出的N个所述白噪声信号的第i个本征模态函数分量的近似熵和平均周期对数，且调用所述二维联合概率密度函数确定模块，得出白噪声信号的第i个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数，记作f(p_i,t_i)；其中，p_i和t_i分别表示白噪声信号的第i个本征模态函数分量的近似熵与平均周期对数；

步骤Ⅴ、白噪声近似熵分布区间模板获取：根据步骤Ⅳ中获得的白噪声信号的n个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数，获得白噪声信号的n个本征模态函数分量的置信区间；所述置信区间的置信水平为c％，其中c％＝92％～99％；

其中，白噪声信号的n个本征模态函数分量的置信区间组成所述白噪声近似熵分布区间模板；

对步骤一中分解出的第i个本征模态函数分量c_i(t)进行有效性识别时，包括以下步骤：

步骤201、近似熵与平均周期自然对数计算：调用步骤Ⅲ中所述近似熵计算模块，计算得出c_i(t)的近似熵；同时，调用步骤Ⅲ中所述平均周期计算模块，计算得出c_i(t)的平均周期对数；

步骤202、有效性识别：根据步骤201中计算得出的c_i(t)的近似熵与平均周期自然对数，判断c_i(t)是否落入步骤Ⅴ中获得的白噪声信号的第i个本征模态函数分量的置信区间内：当落入步骤Ⅴ中获得的白噪声信号的第i个本征模态函数分量的置信区间内时，说明c_i(t)为有效本征模态函数分量；否则，说明c_i(t)为虚假本征模态函数分量；

步骤203、多次重复步骤201至步骤202，直至完成步骤一中分解出的n个本征模态函数分量的有效性识别过程，获得振动信号x(t)经验模态分解后的所有有效本征模态函数分量；

步骤三、信号重构：利用步骤203中获得的振动信号x(t)经验模态分解后的所有有效本征模态函数分量和步骤一中分解出的趋势项，重构出振动信号x'(t)；x'(t)为对振动信号x(t)进行矫正后获得的信号。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：步骤一中进行振动信号经验模态分解之前，先调用归一化处理模块，对振动信号x(t)进行归一化处理；再调用所述经验模态分解模块，对归一化处理后的振动信号x(t)进行经验模态分解；

步骤Ⅱ中进行经验模态分解之前，先调用所述归一化处理模块，对步骤Ⅰ中各白噪声信号分别进行归一化处理；再调用所述经验模态分解模块，对归一化处理后的各白噪声信号分别进行经验模态分解。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：步骤Ⅲ中调用所述近似熵计算模块对各白噪声信号的n个本征模态函数分量的近似熵进行计算时，所有本征模态函数分量的近似熵的计算方法均相同；

其中，c_i,j(t)的近似熵，记作IMF-ApEn_i,j；对c_i,j(t)的近似熵进行计算时，过程如下：

步骤Ⅲ-1、排列组成m维向量O(k)：按照c_i,j(t)中的各数据先后顺序，将c_i,j(t)中的m个数据排列组成M₀-m+1个m维矢量；

其中，c_i,j(t)中包含M₀个数据且M₀个数据分别为c_i,j(1)、c_i,j(2)、…、c_i,j(M₀)，M₀为c_i,j(t)的数据长度且其为正整数；m为预先设定的嵌入维数，m为正整数且m≥2；

M₀-m+1个m维矢量中第k个m维矢量，记作O(k)；

并且，O(k)＝[c_i,j(k),c_i,j(k+1),…,c_i,j(k+m-1)](2-0)，其中k为正整数且k＝1、2、…、M₀-m+1；

步骤Ⅲ-2、Hausdorff空间距离计算：调用Hausdorff空间距离计算模块，计算得出M₀-m+1个所述m维矢量中任一个m维矢量与其它M₀-m个m维矢量之间的Hausdorff空间距离；

其中，m维矢量O(k)与m维矢量O(g)之间的Hausdorff空间距离，记作d[O(k),O(g)]；m维矢量O(g)为步骤Ⅲ-1中M₀-m+1个m维矢量中第g个m维矢量，g为正整数且g＝1、2、…、M₀-m+1，g≠k；

步骤Ⅲ-3、m维矢量间相似度计算：根据预先限定的相似容限r，计算得出M₀-m+1个所述m维矢量中任一个m维矢量与其它M₀-m个m维矢量之间的相似度；

其中，m维矢量O(k)与其它M₀-m个m维矢量之间的相似度，记作根据公式 $C_k^m\left(r\right)=\underset{g=1}{\overset{N-m+1}{\Σ}}\mathrm{\Θ}\{r-d\[O\left(k\right),O\left(g\right)\]\}/(M_0-m+1)---(2-2),$ 进行计算；r＞0；

式(2-2)中， $\mathrm{\&Theta;}\{r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\}=\left\{\begin{array}{cc}1,& r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\&GreaterEqual;0\\ 0,& r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;<0\end{array}\right.---(2-3);$

步骤Ⅲ-4、Φ^m(r)计算：根据步骤Ⅲ-3中计算得出的M₀-m+1个所述m维矢量中任一个m维矢量与其它M₀-m个m维矢量之间的相似度，并按照公式 ${\mathrm{\&Phi;}}^m\left(r\right)=\frac{1}{M-m-1}\underset{k=1}{\overset{N-m+1}{\&Sigma;}}{\mathrm{lnC}}_k^m\left(r\right)---(2-4),$ 计算得出Φ^m(r)；

步骤Ⅲ-5、Φ^m+1(r)计算：按照步骤Ⅲ-1至步骤Ⅲ-4中所述的方法，计算得出Φ^m+1(r)；

步骤Ⅲ-6、近似熵计算：根据公式IMF-ApEn_i,j＝Φ^m(r)-Φ^m+1(r)(2-5)，计算得出IMF-ApEn_i,j。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：所述的m＝2，r＝c₀×S_td[c_i,j(t)]；其中，c₀＝0.1～0.25，S_td[c_i,j(t)]为c_i,j(t)的标准差。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：步骤Ⅲ中调用所述平均周期计算模块，对各白噪声信号的n个本征模态函数分量的平均周期进行计算时，所有本征模态函数分量的平均周期的计算方法均相同；

其中，对c_i,j(t)的平均周期进行计算时，采用快速傅里叶变换方法且按照公式进行计算；其中，F[c_i,j(t)]表示c_i,j(t)的快速傅里叶变换。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：所述的 $F\&lsqb;c_{i,j}\left(t\right)\&rsqb;={\max }_{h=1}^{M}F\&lsqb;c_{i,j}\left(h\right)\&rsqb;---(3-2);$

c_i,j(t)中包含M个数据，M为c_i,j(t)的数据长度且其为正整数；M个数据中第h个数据，记作c_i,j(h)；其中，h为正整数且h＝1、2、…、M。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：步骤Ⅰ中多个所述随机白噪声信号的数据长度均与振动信号x(t)的数据长度相等。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：步骤Ⅰ中所述的N≥200。

上述基于白噪声统计特征的锚杆振动信号矫正方法，其特征是：所述锚杆布设在煤矿开采过程中综采面两侧的工作面运输巷或工作面回风巷内，所述锚杆振动检测单元安装在锚杆的外端。

本发明与现有技术相比具有以下优点：

1、方法步骤简单、设计合理且实现方便、使用效果好，能在短时间内(具体是几分钟内)完成信号矫正过程。

2、利用白噪声近似熵分布区间模板作为“筛子”提取振动信号经EMD分解后的有效信息分量，由于经验模态(EMD)分解能保留原信号的非线性、非平稳性的特征，因而先对待矫正信号进行EMD分解；由于EMD分解由于噪声干扰会引起的过分解现象，因而利用预先获取的白噪声近似熵分布区间模板进识别出有效本征模态函数分量，这样能够在待矫正信号的先验特征未知的条件下提取信号的有效信号分量，因而信号矫正结果准确、可靠，并且实现简便，数据处理速度快。

3、两个锚杆振动测量阵列的各锚杆振动检测单元与近场煤岩层析成像工作站之间均以无线方式进行通信，并且各锚杆振动检测单元分别布设在工作面运输巷和工作面回风巷内的锚杆上，因而形成以无线通信方式构建的煤岩近场层析成像监测网络，方便地实现了综采面煤岩振动工况监测信号的大范围覆盖，解决了***布线对采煤工作的干扰问题。

4、通过综采面两侧回采巷道内的锚杆振动检测单元对综采面近场煤岩对应监测位置的振动状况进行监测，各锚杆振动检测单元所处位置为振动监测位置，并且振动监测位置固定且位置信息易于确定，因而能对综采面近场煤岩不同监测位置的振动状况进行有效监测。由于两个锚杆振动测量阵列分别布设在回采巷道中的工作面运输巷和工作面回风巷内，因而两个锚杆振动测量阵列能对综采面近场围岩的振动状况进行有效监测；并且，锚杆振动测量阵列中包含多个锚杆振动检测单元，各锚杆振动检测单元均布设在巷道壁上所布设锚杆的外端，因而锚杆振动检测单元布设简便、实现方便且监测数据可靠。采用基于白噪声统计特征的振动有效量提取方法对锚杆振动检测单元的振动监测结果进行矫正，能进一步保证振动监测结果的准确性和可靠性，提高动态层析成像效果。

综上所述，本发明方法步骤简单、设计合理且实现方便、使用效果好，能对锚杆振动信号进行有效矫正，获得锚杆的真实振动状态。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明的方法流程框图。

图1-1为本发明动态层析成像***的电路原理框图。

图2为本发明锚杆振动检测单元的电路原理框图。

图2-1为本发明无线汇集节点的电路原理框图。

图3为本发明工作面运输巷内锚杆振动检测单元的布设位置示意图。

图4为本发明进行动态层析成像时的方法流程框图。

图5为本发明所获取的综采面近场煤岩层析图像。

图5-1为本发明综采面前方第一个振动检测断面上两个测点的布设位置示意图。

图5-2为本发明综采面前方测点的布设位置示意图。

图6为本发明所构建需矫正信号及其EMD分解结果的示意图。

图7-1为本发明所构建模拟信号经EMD分解后的IMF1的近似熵和平均周期对数分布示意图。

图7-2为本发明所构建模拟信号经EMD分解后的IMF2的近似熵和平均周期对数分布示意图。

图7-3为本发明所构建模拟信号经EMD分解后的IMF3的近似熵和平均周期对数分布示意图。

图7-4为本发明所构建模拟信号经EMD分解后的IMF4的近似熵和平均周期对数分布示意图。

图7-5为本发明所构建模拟信号经EMD分解后的IMF5的近似熵和平均周期对数分布示意图。

图7-6为本发明所构建模拟信号经EMD分解后的IMF6的近似熵和平均周期对数分布示意图。

附图标记说明:

1-1—开采位置检测单元；1-2—锚杆振动测量阵列；

2—近场煤岩层析成像工作站；2-1—上位机；

2-2—第二无线通信模块；3—采煤机；4—工作面运输巷；

5—工作面回风巷；6—锚杆；7—综采面；

8—锚杆振动检测单元；8-1—振动传感器；8-2—信号调理电路；

8-3—第一处理器；8-4—无线信号发射模块；

8-5—第一数据存储模块；8-6—信号采集电路；

9—锚杆测量阵列汇集节点；9-1—第二处理器；

9-2—无线信号接收模块；9-3—第一无线通信模块；

9-4—第二数据存储模块。

具体实施方式

如图1所示的一种基于白噪声统计特征的锚杆振动信号矫正方法，采用数据处理设备对布设于被测试锚杆6上的锚杆振动检测单元8所检测振动信号x(t)进行矫正，包括以下步骤：

步骤一、振动信号经验模态分解：调用经验模态分解模块，对锚杆振动检测单元8所检测振动信号x(t)进行经验模态分解，获得n个本征模态函数分量和一个趋势项；n个本征模态函数分量和一个趋势项均为离散函数；

其中，振动信号x(t)为锚杆6的振动位移或振动加速度随时间变化的函数；振动信号x(t)的经验模态分解结果为：式(1)中c_i(t)为振动信号x(t)经验模态分解后获得的第i个本征模态函数分量，i为正整数且其本征模态函数分量的编号，i＝1、2、…、n；r_n(t)为振动信号x(t)经验模态分解后获得的趋势项；

步骤二、本征模态函数分量有效性识别：调用本征模态函数分量有效性识别模块，且根据预先建立的白噪声近似熵分布区间模板，对步骤一中分解出的振动信号x(t)的n个本征模态函数分量分别进行有效性识别；n个本征模态函数分量的有效性识别方法均相同；

对所述白噪声近似熵分布区间模板进行建立时，采用所述数据处理设备进行建立，过程如下：

步骤Ⅰ、样本获取：产生多个随机白噪声信号作为样本，多个所述白噪声信号的数据长度均相同；

其中，产生的第j个随机白噪声信号，记作ζ_j(t)；j为随机白噪声信号的编号且其为正整数，j＝1、2、…、N；N为本步骤中所产生随机白噪声信号的总数量；

步骤Ⅱ、经验模态分解：调用所述经验模态分解模块，对步骤Ⅰ中各白噪声信号分别进行经验模态分解，获得各白噪声信号的n个本征模态函数分量和一个趋势项；n个本征模态函数分量和一个趋势项均为离散函数；

其中，ζ_j(t)的经验模态分解结果为：式(1)中c_i,j(t)为ζ_j(t)经验模态分解后获得的第i个本征模态函数分量，r_n,j(t)为ζ_j(t)经验模态分解后获得的趋势项；

步骤Ⅲ、近似熵与平均周期自然对数计算：调用近似熵计算模块，对步骤Ⅱ中各白噪声信号的n个本征模态函数分量的近似熵分别进行计算；同时，调用平均周期计算模块，对各白噪声信号的n个本征模态函数分量的平均周期分别进行计算，并求解出各白噪声信号的n个本征模态函数分量的平均周期对数；其中，平均周期对数为平均周期的自然对数；

步骤Ⅳ、二维联合概率密度函数确定：根据步骤Ⅲ中计算得出的各白噪声信号的n个本征模态函数分量的近似熵和平均周期对数，调用二维联合概率密度函数确定模块，获取白噪声信号的n个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数；

其中，根据步骤Ⅲ中计算得出的N个所述白噪声信号的第i个本征模态函数分量的近似熵和平均周期对数，且调用所述二维联合概率密度函数确定模块，得出白噪声信号的第i个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数，记作f(p_i,t_i)；其中，p_i和t_i分别表示白噪声信号的第i个本征模态函数分量的近似熵与平均周期对数；

步骤Ⅴ、白噪声近似熵分布区间模板获取：根据步骤Ⅳ中获得的白噪声信号的n个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数，获得白噪声信号的n个本征模态函数分量的置信区间；所述置信区间的置信水平为c％，其中c％＝92％～99％；

其中，白噪声信号的n个本征模态函数分量的置信区间组成所述白噪声近似熵分布区间模板；

对步骤一中分解出的第i个本征模态函数分量c_i(t)进行有效性识别时，包括以下步骤：

步骤201、近似熵与平均周期自然对数计算：调用步骤Ⅲ中所述近似熵计算模块，计算得出c_i(t)的近似熵；同时，调用步骤Ⅲ中所述平均周期计算模块，计算得出c_i(t)的平均周期对数；

步骤202、有效性识别：根据步骤201中计算得出的c_i(t)的近似熵与平均周期自然对数，判断c_i(t)是否落入步骤Ⅴ中获得的白噪声信号的第i个本征模态函数分量的置信区间内：当落入步骤Ⅴ中获得的白噪声信号的第i个本征模态函数分量的置信区间内时，说明c_i(t)为有效本征模态函数分量；否则，说明c_i(t)为虚假本征模态函数分量；

步骤203、多次重复步骤201至步骤202，直至完成步骤一中分解出的n个本征模态函数分量的有效性识别过程，获得振动信号x(t)经验模态分解后的所有有效本征模态函数分量；

步骤三、信号重构：利用步骤203中获得的振动信号x(t)经验模态分解后的所有有效本征模态函数分量和步骤一中分解出的趋势项，重构出振动信号x'(t)；x'(t)为对振动信号x(t)进行矫正后获得的信号。

本实施例中，c％＝95％。

实际使用时，可根据具体需要，对c％的取值大小进行相应调整。

经验模态分解后，获得的n个本征模态函数分量(即IMF分量)的频率分辨率是不同的，这点不同于小波分析中时间与频率相互影响特性。因而，经验模态(EMD)分解保留了原信号的非线性、非平稳性的特征。

本实施例中，步骤一中进行振动信号经验模态分解之前，先调用归一化处理模块，对振动信号x(t)进行归一化处理；再调用所述经验模态分解模块，对归一化处理后的振动信号x(t)进行经验模态分解；

步骤Ⅱ中进行经验模态分解之前，先调用所述归一化处理模块，对步骤Ⅰ中各白噪声信号分别进行归一化处理；再调用所述经验模态分解模块，对归一化处理后的各白噪声信号分别进行经验模态分解。

本实施例中，步骤201中对c_i(t)的近似熵和平均周期对数进行计算时，所采用的近似熵和平均周期对数的计算方法与步骤Ⅲ中采用的计算方法相同。

其中，近似熵(ApproximateEntropy，ApEn)是由StevePincus提出的一种从统计学角度反映随机时间序列复杂度的指标。

本实施例中，步骤Ⅲ中调用所述近似熵计算模块对各白噪声信号的n个本征模态函数分量的近似熵进行计算时，所有本征模态函数分量的近似熵的计算方法均相同；

其中，c_i,j(t)的近似熵，记作IMF-ApEn_i,j；对c_i,j(t)的近似熵进行计算时，过程如下：

步骤Ⅲ-1、排列组成m维向量O(k)：按照c_i,j(t)中的各数据先后顺序，将c_i,j(t)中的m个数据排列组成M₀-m+1个m维矢量；

其中，c_i,j(t)中包含M₀个数据且M₀个数据分别为c_i,j(1)、c_i,j(2)、…、c_i,j(M₀)，M₀为c_i,j(t)的数据长度且其为正整数；m为预先设定的嵌入维数，m为正整数且m≥2；

M₀-m+1个m维矢量中第k个m维矢量，记作O(k)；

并且，O(k)＝[c_i,j(k),c_i,j(k+1),…,c_i,j(k+m-1)](2-0)，其中k为正整数且k＝1、2、…、M₀-m+1；

步骤Ⅲ-2、Hausdorff空间距离计算：调用Hausdorff空间距离计算模块，计算得出M-m+1个所述m维矢量中任一个m维矢量与其它M-m个m维矢量之间的Hausdorff空间距离；

其中，m维矢量O(k)与m维矢量O(g)之间的Hausdorff空间距离，记作d[O(k),O(g)]；m维矢量O(g)为步骤Ⅲ-1中M-m+1个m维矢量中第g个m维矢量，g为正整数且g＝1、2、…、M-m+1，g≠k；

步骤Ⅲ-3、m维矢量间相似度计算：根据预先限定的相似容限r，计算得出M-m+1个所述m维矢量中任一个m维矢量与其它M-m个m维矢量之间的相似度；

其中，m维矢量O(k)与其它M-m个m维矢量之间的相似度，记作 根据公式 $C_k^m\left(r\right)=\underset{g=1}{\overset{N-m+1}{\&Sigma;}}\mathrm{\&Theta;}\{r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\}/(M-m+1)---(2-2),$ 进行计算；r＞0；

式(2-2)中， $\mathrm{\&Theta;}\{r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\}=\left\{\begin{array}{cc}1,& r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\&GreaterEqual;0\\ 0,& r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;<0\end{array}\right.---(2-3);$

步骤Ⅲ-4、Φ^m(r)计算：根据步骤Ⅲ-3中计算得出的M-m+1个所述m维矢量中任一个m维矢量与其它M-m个m维矢量之间的相似度，并按照公式 ${\mathrm{\&Phi;}}^m\left(r\right)=\frac{1}{M-m-1}\underset{k=1}{\overset{N-m+1}{\&Sigma;}}{\mathrm{lnC}}_k^m\left(r\right)---(2-4),$ 计算得出Φ^m(r)；Φ^m(r)为平均值；

步骤Ⅲ-5、Φ^m+1(r)计算：按照步骤Ⅲ-1至步骤Ⅲ-4中所述的方法，计算得出Φ^m+1(r)；

步骤Ⅲ-6、近似熵计算：根据公式IMF-ApEn_i,j＝Φ^m(r)-Φ^m+1(r)(2-5)，计算得出IMF-ApEn_i,j。

其中，步骤Ⅲ-2中所述的Hausdorff空间距离简称“Hausdorff距离”，即

其中，步骤Ⅲ-2中所述的Hausdorff空间距离，简称“Hausdorff距离”，即豪斯多夫距离，是豪斯多夫距离量度度量空间中真子集之间的距离。豪斯多夫距离量度度量空间中真子集之间的距离。

由上述近似熵的计算过程可知，近似熵的计算值与嵌入维度m和相似容限r有关，若相似容限r减少，则满足相似条件的模式就会减少，其近似熵的估计准确度就会变差；相反，如果相似容限r增加，其相似模式过多，时间序列的细节信息又会丢失。因而，嵌入维度m和相似容限r的取值至关重要。

本实施例中，所述的m＝2。

并且，r＝c₀×S_td[c_i,j(t)]；其中，c₀＝0.1～0.25，S_td[c_i,j(t)]为c_i,j(t)的标准差。c_i,j(t)为一个时间序列。本实施例中，c₀＝0.2。

实际使用时，可根据具体需要，对m和c₀的取值进行相应调整。

由于白噪声或分形高斯噪声的经验模态分解(即EMD分解)相当于一组二分滤波器组，其各IMF分量的傅里叶变化谱图沿着周期或频率的对数近似为单一的形状。考虑到尺度丢失效应，白噪声的IMF分量个数应该小于N₀为噪声序列的长度。IMF分量(即本征模态函数分量)被定义为由局部极大值和极小值形成的包络函数，其极值数目与过零点相同。IMF分量的平均周期从其定义上可以看出是由其局部峰值的个数决定。此处，通过快速傅里叶变换(FastFourierTransform，FFT)方法对平均周期进行计算。其中，平均周期也称为平均时间周期。

本实施例中，步骤Ⅲ中调用所述平均周期计算模块，对各白噪声信号的n个本征模态函数分量的平均周期进行计算时，所有本征模态函数分量的平均周期的计算方法均相同；

其中，对c_i,j(t)的平均周期进行计算时，采用快速傅里叶变换方法且按照公式进行计算；其中，F[c_i,j(t)]表示c_i,j(t)的快速傅里叶变换。

本实施例中，所述的 $F\&lsqb;c_{i,j}\left(t\right)\&rsqb;={\max }_{h=1}^{M_0}F\&lsqb;c_{i,j}\left(h\right)\&rsqb;---(3-2);$

c_i,j(t)中包含M个数据，M为c_i,j(t)的数据长度且其为正整数；M个数据中第h个数据，记作c_i,j(h)；其中，h为正整数且h＝1、2、…、M。

本实施例中，步骤Ⅰ中多个所述随机白噪声信号的数据长度均与振动信号x(t)的数据长度相等。其中，宽信号的时宽也称数据长度。

并且，步骤Ⅰ中所述的N≥200。实际使用时，N的取值越大，步骤Ⅳ中所确定的二维联合概率密度函数的精度越高。

步骤Ⅳ中进行二维联合概率密度函数确定时，所确定的二维联合概率密度函数为二维随机变量的联合概率密度函数，具体是近似熵和平均周期对数的联合概率密度函数。

步骤Ⅳ中进行二维联合概率密度函数确定时，样本为根据步骤Ⅲ中计算得出的各白噪声信号的n个本征模态函数分量的近似熵和平均周期对数，并且需对计算得出的各白噪声信号的n个本征模态函数分量的近似熵和平均周期对数分别进行编号；

并且，根据公式进行确定，其中n₀为样本数量，h₀为核估计窗宽，K()表示二维核函数，x和y分别表示近似熵，为样本中第个近似熵，为第个平均周期对数。其中，i₀为正整数且i₀＝1、2、…、n₀。

由于EMD分解将信号自适应地分解为一系列IMF分量，但是其中有一部分IMF分量并没有从物理上解释信号的特征，该分量被称之为虚假分量。EMD分解利用信号的局部极值信息筛选IMF分量是导致虚假模态产生的一个原因，另外，由于信号中混入噪声干扰，使得真实信号中出现虚假极值点，同样导致了一些高频噪声和低频分量等虚假模式。

由上述内容可知，采用基于白噪声统计特征的振动有效量提取方法，对锚杆振动检测单元8检测到的振动信号x(t)进行矫正时，待步骤Ⅳ中白噪声信号的n个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数确定后，获得白噪声信号的n个本征模态函数分量的置信区间组成所述白噪声近似熵分布区间模板，并利用白噪声近似熵分布区间模板作为“筛子”提取振动信号x(t)经EMD分解后的有效信息分量，从而完成振动信号x(t)的矫正过程。

并且，利用白噪声近似熵分布区间模板作为“筛子”提取振动信号x(t)经EMD分解后的有效信息分量时，先用对振动信号x(t)进行归一化处理，再对振动信号x(t)进行EMD分解，得到不同时间尺度下的信号模态分量(即IMF分量)，再计算各IMF分量的近似熵和平均周期对数，其中近似熵记作ApEn，平均周期对数记作Log-M-P；然后，再根据所获得的白噪声近似熵分布区间模板，比较振动信号x(t)经EMD分解后的各IMF分量的近似熵和平均周期对数，以判别振动信号x(t)中各IMF分量中是否包含有虚假模式分量，其中落在白噪声近似熵分布区间模板内的IMF分量为高频噪声IMF分量与低频虚假IMF分量，落在白噪声近似熵分布区间模板外的IMF分量为振动信号x(t)的有效信号分量(effectiveintrinsicmodefunctions，EIMF)。

此处，构建一个模拟信号y(n₁)对本发明采用的基于白噪声统计特征的振动有效量提取方法进行验证，信号y(n₁)由60Hz的正弦信号、15Hz的正弦衰减信号和白噪声信号组成，y(n₁)＝s(n₁)+3·sin(2π·60·n₁·T₁)+v(n₁)(1-1)；其中， $s\left(n_1\right)=\left\{\begin{array}{cc}\sin (2\mathrm{\&pi;}\&CenterDot;15\&CenterDot;n_1\&CenterDot;T_1)e^{-a(n_1\&CenterDot;T_1-b)}& n_1{T}_1<b\\ 0& n_1{T}_1\&GreaterEqual;b\end{array}\right.---(1-2),$ 式(1-2)中T₁为采样时间间隔且T₁＝0.001s，v(n₁)为均值为0且方差为1的白噪声信号。如图6所示，该模拟信号y(n₁)的EMD分解结果中，第2个IMF分量(即IMF2)和第3个IMF分量(即IMF3)代表正弦信号，第4个IMF分量(即IMF4)代表正弦衰减信号，IMF2、IMF3和IMF4均为有效信号分量(即EIMF)。而第1个IMF分量(即IMF1)能从信号的物理意义上进行解释，并不是信号的有效构成模式，是由于噪声干扰导致的为虚假的模态分量；第5个IMF分量(即IMF5)为EMD过分解而产生的虚假低频模式。图6中，signal表示模拟信号y(n₁)，Res表示模拟信号y(n₁)经EMD分解后的趋势项。

结合图7-1、图7-2、图7-3、图7-4、图7-5和图7-6中模拟信号y(n₁)的6个IMF分量的近似熵和平均周期对数分布能看出：模拟信号y(n₁)的IMF1、IMF5和IMF6均落在白噪声近似熵分布区间模板(具体是白噪声信号的对应本征模态函数分量的置信区间)内，因而IMF1、IMF5和IMF6分别为高频噪声分量和低频虚假分量；模拟信号y(n₁)的IMF2、IMF3和IMF4落在白噪声近似熵分布区间模板(具体是白噪声信号的对应本征模态函数分量的置信区间)外，属于有效本征模态函数分量(即EIMF)。

由于EMD分解是数据自驱动的，其特性受数据自身特性的影响，EMD分解由于噪声干扰会引起的过分解现象(如上述IMF2和IMF3)。IMF2和IMF3具有几乎相同的平均周期，因而二者组合在一起能够更精确的代表60Hz的正弦信号。由此说明，本发明采用的基于白噪声统计特征的振动有效量提取方法能够在待测信号的先验特征未知的条件下提取信号的有效信号分量。

本实施例中，所述锚杆6布设在煤矿开采过程中综采面7两侧的工作面运输巷4或工作面回风巷5内，所述锚杆振动检测单元8安装在锚杆6的外端。如图1-1所示，所述工作面运输巷4和工作面回风巷5内均布设有多个所述锚杆振动检测单元8，并且工作面运输巷4和工作面回风巷5内的多个所述锚杆振动检测单元8均形成一个锚杆振动测量阵列1-2。同时，在对综采面7进行开采的采煤机3上布设开采位置检测单元1-1，所述开采位置检测单元1-1为采煤过程中对振动激励源的位置进行实时检测的位置检测单元。所述综采面7为回采工作面，所述振动激励源为对综采面7进行开采的采煤机3，所述工作面运输巷4和工作面回风巷5均为回采巷道。所述开采位置检测单元1-1和两个所述锚杆振动测量阵列1-2组成采煤过程中对所述振动激励源作用下综采面近场煤岩的振动状态进行实时检测的阵列信息同步采集装置，所述阵列信息同步采集装置与近场煤岩层析成像工作站2组成煤机采动激励下综采面近场煤岩动态层析成像***，且所述阵列信息同步采集装置与近场煤岩层析成像工作站2相接。所述工作面运输巷4和工作面回风巷5的巷道支护结构均包括多个布设在巷道壁上的锚杆6。所述综采面近场煤岩为位于综采面7前方的探测区域内的煤岩体，所述阵列信息同步采集装置位于综采面7前方。所述近场煤岩层析成像工作站2包括上位机2-1。

每个所述锚杆振动测量阵列1-2的多个所述锚杆振动检测单元8均布设在工作面运输巷4或工作面回风巷5内的多个所述锚杆6外端。两个所述锚杆振动测量阵列1-2中的所有锚杆振动检测单元8均与上位机2-1相接。所述开采位置检测单元1-1布设在采煤机3上且其与上位机2-1相接。

根据《煤矿安全规程》的规定，综采工作面前方安全距离指的是直线距离不小于100米；拐1个直角弯时不小于70米。因而，综采面7前方的探测区域为综采工作面前方的安全距离范围内的区域。煤矿开采过程中，回采巷道是指形成采煤工作面及为其服务的巷道，如开切眼、工作面运输巷、工作面回风巷等。

层析成像技术是借鉴医学CT，根据射线扫描，对所得到的信息进行反演计算，重建被测范围内岩体弹性波和电磁波参数分布规律的图像，从而达到圈定地质异常体的一种物探反演解释方法。

本实施例中，所述开采位置检测单元1-1为对采煤机3的截割滚筒位置进行实时检测的位置检测单元。

实际布设安装时，所述开采位置检测单元1-1安装在采煤机3的截割滚筒上，具体是安装在采煤机3的截割滚筒的滚筒轴。

实际使用时，两个所述锚杆振动测量阵列1-2中的所有锚杆振动检测单元8沿综采面7的推进方向由前至后分为多个振动检测组，多个所述振动检测组分别布设在多个振动检测断面上，所述振动检测断面与综采面7呈平行布设，多个所述振动检测断面沿综采面7的推进方向由前至后进行布设；每个所述振动检测组均包括布设在同一振动检测断面上的多个所述锚杆振动检测单元8，每个所述振动检测组中布设于工作面运输巷4内的锚杆振动检测单元8的数量为一个或多个，每个所述振动检测组中布设于工作面回风巷5内的锚杆振动检测单元8的数量为一个或多个。

本实施例中，每个所述振动检测组中布设于工作面运输巷4内的锚杆振动检测单元8的数量为多个，每个所述振动检测组中布设于工作面回风巷5内的锚杆振动检测单元8的数量为多个。并且，每个所述振动检测组中布设于工作面运输巷4内的锚杆振动检测单元8的数量与布设于工作面回风巷5内的锚杆振动检测单元8的数量相同且二者呈对称布设。

实际使用时，可根据具体需要，对每个所述振动检测组中布设于工作面运输巷4内的锚杆振动检测单元8的数量和布设于工作面回风巷5内的锚杆振动检测单元8的数量进行相应调整。

本实施例中，所述工作面运输巷4和工作面回风巷5内的多个所述锚杆6均由前至后分为多个锚杆组，多个所述锚杆组分别布设在多个所述振动检测断面上。每个所述锚杆组均包括布设在同一振动检测断面上的多个所述锚杆6，每个所述锚杆组中的多个所述锚杆6均包括多个由左至右布设在巷道顶板上的顶板锚杆、多个由上至下布设在内侧巷道帮部的内侧巷帮锚杆和多个由上至下布设在外侧巷道帮部的外侧巷帮锚杆，所述工作面运输巷4和工作面回风巷5的内侧巷道帮部均为靠近综采面7一侧的巷道帮部。两个所述锚杆振动测量阵列1-2中的所有锚杆振动检测单元8均布设在所述内侧巷帮锚杆上，每个所述内侧巷帮锚杆上所布设锚杆振动检测单元8的数量均为一个。两个所述锚杆振动测量阵列1-2组成对采煤机3采煤过程中产生的应力波信号进行实时检测的应力波检测装置。

本实施例中，两个所述锚杆振动测量阵列1-2中的所有锚杆振动检测单元8与上位机2-1之间均以无线方式进行通信。

如图2所示，所述锚杆振动检测单元8为无线振动检测单元且其包括振动传感器8-1、第一处理器8-3、与振动传感器8-1相接的信号调理电路8-2以及分别与第一处理器8-3相接的无线信号发射模块8-4和第一数据存储模块8-5，所述信号调理电路8-2与第一处理器8-3相接。

本实施例中，所述振动传感器8-1包括对所布设锚杆6径向上的振动状态进行实时检测的径向振动传感器和对所布设锚杆6轴向上的振动状态进行实时检测的轴向振动传感器，所述径向振动传感器和轴向振动传感器均与信号调理电路8-2相接。

本实施例中，所述径向振动传感器和轴向振动传感器均为振动位移传感器。

实际使用时，所述径向振动传感器和轴向振动传感器也可以采用振动速度传感器或振动加速度传感器。

本实施例中，所述无线振动检测单元还包括信号采集电路8-6，所述信号调理电路8-2通过信号采集电路8-6与第一处理器8-3相接，所述信号采集电路8-6为A/D转换模块。

并且，所述无线振动检测单元还包括为各用电单元进行供电的电源模块，所述振动传感器8-1、信号调理电路8-2、第一处理器8-3、无线发射模块8-4、数据存储模块8-5和信号采集电路8-6均与所述电源模块相接。

本实施例中，本发明所述的煤机采动激励下综采面近场煤岩动态层析成像***，还包括两个锚杆测量阵列汇集节点9，两个所述锚杆振动测量阵列1-2分别通过两个所述锚杆测量阵列汇集节点9与上位机2-1进行通信。两个所述锚杆测量阵列汇集节点9分别为布设在工作面运输巷4和工作面回风巷5内的第一汇集节点和第二汇集节点，所述第一汇集节点和第二汇集节点均为无线汇集节点。如图2-1所示，所述无线汇集节点包括第二处理器9-1以及分别与第二处理器9-1相接的无线信号接收模块9-2、第一无线通信模块9-3和第二数据存储模块9-4，所述第一处理器8-3通过无线信号发射模块8-4和无线信号接收模块9-2与第二处理器9-1进行通信；所述近场煤岩层析成像工作站2还包括与上位机2-1相接的第二无线通信模块2-2，所述第二处理器9-1通过第一无线通信模块9-3和第二无线通信模块2-2与上位机2-1进行通信。

实际使用时，所述锚杆测量阵列汇集节点9将与其相接的锚杆测量阵列汇集节点9的检测信号传送至上位机2-1之前，需先采用第二处理器9-1对锚杆测量阵列汇集节点9中多个所述锚杆振动检测单元8分别进行地址编码。

本实施例中，所述近场煤岩层析成像工作站2安装在工作面运输巷4中。实际进行安装布设时，可根据具体需要，对近场煤岩层析成像工作站2的安装位置进行相应调整。

本实施例中，所述锚杆振动检测单元8中的第一处理器8-3、信号调理电路8-2、无线信号发射模块8-4、第一数据存储模块8-5和信号采集电路8-6均布设在一个防爆壳体内，且该防爆壳体通过紧固螺钉固定在所述内侧巷帮锚杆上，详见图3。

本实施例中，两个所述锚杆测量阵列汇集节点9分别布设在工作面运输巷4和工作面回风巷5内。

实际使用时，通过两个所述锚杆振动测量阵列1-2中的多个锚杆振动检测单元8对采煤机3在不同开采位置进行采煤时产生的应力波引起的锚杆6的振动信号进行实时检测；两个所述锚杆测量阵列汇集节点9分别对两个所述锚杆振动测量阵列1-2中多个锚杆振动检测单元8所检测信号进行无线中继传送。本实施例中，所述无线信号发射模块8-4和无线信号接收模块9-2为基于TCP/IP网络协议的无线通信模块。

实际使用时，锚杆振动测量阵列1-2中的多个所述锚杆振动检测单元8对采煤机3在一个开采位置进行采煤时多个所述锚杆振动检测单元8所检测信号(即采煤机3在一个开采位置的振动检测结果)进行同步采集。并且，所述锚杆测量阵列汇集节点9对锚杆振动测量阵列1-2中多个所述锚杆振动检测单元8所检测信号进行传送时，将采煤机3在同一开采位置进行采煤时多个所述锚杆振动检测单元8所检测信号进行同步传送。

如图4所示，实际进行煤机采动激励下综采面近场煤岩动态层析成像时，包括以下步骤：

步骤A、锚杆振动检测单元及开采位置检测单元布设：对需开采煤层进行回采之前，先在工作面运输巷4和工作面回风巷5内由前至后分别布设多个所述锚杆振动检测单元8，并在采煤机3上布设开采位置检测单元1-1；每个所述锚杆振动检测单元8的布设位置均为一个测点；

步骤B、煤机采动激励下综采面近场煤岩动态层析成像：沿综采面7的推进方向，采用采煤机3由后向前对需开采煤层进行回采；回采过程中，对需开采煤层中沿综采面7推进方向上各位置的综采面近场煤岩分别进行层析成像，且沿综采面7推进方向上各位置的综采面近场煤岩的层析成像方法均相同；其中，沿综采面7推进方向上的各位置，指的是沿综采面7推进方向上的各综采面推进位置；

对需开采煤层中沿综采面7推进方向上任一位置的综采面近场煤岩进行层析成像时，过程如下：

步骤B1、综采面开采及振动同步监测：采用采煤机3且沿综采面7的长度方向由后向前对当前综采面进行开采，开采过程中通过开采位置检测单元1-1对采煤机3的开采位置进行实时检测，并将所检测信号同步传送至近场煤岩层析成像工作站2；同时，通过位于当前综采面前方的两个所述锚杆振动测量阵列1-2对振动激励源作用下综采面近场煤岩的振动状态进行实时检测，并将所检测信号同步传送至近场煤岩层析成像工作站2，获得当前综采面开采过程中采煤机3不同开采位置的振动监测结果；

所述当前综采面为当前状态下进行开采的综采面7；

所述采煤机3各开采位置的振动监测结果，均包括开采位置检测单元1-1检测到的该开采位置的位置信号和采煤机3在该开采位置进行采煤时两个所述锚杆振动测量阵列1-2中所有振动检测单元8检测到的振动信号；

步骤B2、综采面近场煤岩层析成像：所述近场煤岩层析成像工作站2接收到步骤B1中采煤机3不同开采位置的振动监测结果后，所述上位机2-1根据接收到采煤机3不同开采位置的振动监测结果，并结合两个所述锚杆振动测量阵列1-2中各锚杆振动检测单元8的布设位置，采用上位机2-1且按照速度层析成像方法，获得沿综采面7推进方向上当前位置的综采面近场煤岩层析图像，完成综采面7当前位置的综采面近场煤岩层析成像过程；

步骤B3、综采面向前推进至下一位置的综采面近场煤岩层析成像：按照步骤B1至步骤B2中所述的方法，完成综采面7向前推进至下一位置的综采面近场煤岩层析成像过程；

步骤B4、一次或多次重复步骤B3，直至完成需开采煤层中沿综采面7推进方向上各位置的综采面近场煤岩层析成像过程。

需说明的是：步骤B2中所述的沿综采面7推进方向上当前位置为综采面7的当前推进位置，也是所述当前综采面所处的位置。因而，步骤B2中沿综采面7推进方向上当前位置的综采面近场煤岩层析图像，为当前综采面开采过程中综采面近场煤岩的层析图像。相应地，步骤B3中综采面7向前推进至下一位置为沿综采面7推进方向将综采面7向前推进后的下一个推进位置。步骤B4中需开采煤层中沿综采面7推进方向上各位置的综采面近场煤岩层析成像过程，指的是需开采煤层中综采面7向前动态推进过程中各推进位置的综采面近场煤岩层析成像过程，完成综采面7向前动态推进过程中的综采面近场煤岩动态层析成像过程。

本实施例中，需开采煤层中一个综采面7的综采面近场煤岩层析成像，详见图5。图5中，Length表示综采面7的长度，Width表示综采面7的宽度，velocity表示速度。

本实施例中，步骤B1中获得的采煤机3不同开采位置的振动监测结果为综采面7当前推进位置开采过程中采煤机3在不同开采位置进行采煤时产生的应力波的检测结果，即所述当前综采面开采过程中采煤机3在不同开采位置进行采煤时产生的应力波的检测结果；

步骤B2中所述近场煤岩层析成像工作站2接收到步骤B1中采煤机3不同开采位置的振动监测结果后，所述上位机2-1先对接收到的采煤机3不同开采位置的振动监测结果进行同步存储；

步骤B2中进行层析成像之前，先采用上位机2-1选取综采面7当前推进位置开采过程中采煤机3的m个开采位置作为有效振动监测位置，并对m个所述振动监测位置的位置信息进行记录；之后，上位机2-1根据m个所述有效振动监测位置的位置信息，从步骤B1中获得的采煤机3不同开采位置的振动监测结果中提取出m个所述有效振动监测位置的振动监测结果；其中，m为正整数且m≥3；

步骤B2中进行层析成像时，所述上位机2-1根据提取出的m个所述有效振动监测位置的振动监测结果，并结合两个所述锚杆振动测量阵列1-2中各锚杆振动检测单元8的布设位置，且按照速度层析成像方法，获得沿综采面7推进方向上当前位置的综采面近场煤岩层析图像；

步骤B4中完成需开采煤层中沿综采面7推进方向上各位置的综采面近场煤岩层析成像过程后，获得综采面7各推进位置开采过程中综采面近场煤岩的层析图像；所获得的综采面7各推进位置开采过程中综采面近场煤岩的层析图像，按照沿综采面7推进方向上的开采先后顺序组成需开采煤层的综采面近场煤岩的动态层析图像。

本实施例中，步骤B2中进行层析成像时，包括以下步骤：

步骤B21、速度模型建立：建立成像区域层析成像用的速度模型，所述速度模型为初始模型且其为网格模型；所述成像区域为当前综采面前方的探测区域且其为当前综采面的综采面近场煤岩所处区域；

步骤B22、振动监测结果分析处理：对提取出的m个所述有效振动监测位置的振动监测结果分别进行分析处理，m个所述有效振动监测位置的振动监测结果的分析处理方法均相同；其中，对任一个所述有效振动监测位置的振动监测结果进行分析处理时，包括以下步骤：

步骤B221、任意两测点间应力波实际走时差计算：根据当前开采位置的振动监测结果，计算得出工作面运输巷4内的各测点与工作面回风巷5内的各测点之间的应力波实际走时差；

其中，工作面运输巷4内的测点p与工作面回风巷5内的测点q之间的应力波实际走时差，记作ΔT^(p,q)；ΔT^(p,q)为当前分析应力波到测点p的初至时间与到达测点q的初至时间的时间差；其中，所述当前分析应力波为采煤机3在当前开采位置进行采煤时产生的应力波；所述当前开采位置为当前分析处理的有效振动监测位置；

其中，步骤一中两个所述锚杆振动测量阵列1-2中所包括锚杆振动检测单元8的总数量为M个，所述综采面7前方探测区域内的测点总数量为M个，M个所述测点分别为测点1、测点2、…、测点M；p为正整数且p＝1、2、…、m1，q为正整数且q＝m1+1、m1+2、…、M，m1为正整数且其为工作面运输巷4内的测点总数量，m2为正整数且其为工作面回风巷5内的测点总数量，m1+m2＝M；

步骤B222、任意两测点间应力波理论走时差计算：先按照射线追踪方法，求解出当前分析应力波穿越成像区域的M条射线的轨迹和旅行时；之后，根据计算得出的M条射线的旅行时，计算得出工作面运输巷4内的各测点与工作面回风巷5内的各测点之间的应力波理论走时差，完成当前开采位置的振动监测结果分析处理过程；

其中，M条射线为当前开采位置分别与M个所述测点之间的射线，M条射线中任一条射线的旅行时均为根据预先建立的所述成像区域的速度模型计算得出的该射线的正演理论走时；M条射线分别为射线1、射线2、…、射线M；

测点p与测点q之间的应力波理论走时差，记作为射线p的旅行时与射线q的旅行时的时间差；其中，射线p为当前分析应力波穿越成像区域到测点p的射线且其为当前开采位置与测点p之间的射线，射线q为当前分析应力波穿越成像区域到测点q的射线且其为当前开采位置与测点q之间的射线；

步骤B221中计算得出的工作面运输巷4内各测点与工作面回风巷5内各测点之间的应力波实际走时差和步骤B222中计算得出的工作面运输巷4内各测点与工作面回风巷5内各测点之间的应力波理论走时差，组成当前开采位置的一组反演求解用模型修正数据；

步骤B223、多次重复步骤B221至步骤B222，直至完成m个所述有效振动监测位置的振动监测结果分析处理过程，获得m组反演求解用模型修正数据；

步骤B23、反演求解：调用SIRT算法，对步骤B21中所述速度模型进行修正，并计算得出修正后的速度模型中各网格的慢度，使得步骤B223中获得的每组反演求解用模型修正数据中测点p与测点q之间的应力波理论走时差与应力波实际走时差之间的偏差均满足预先设计的误差要求，获得采煤机3当前开采位置的综采面近场煤岩的层析图像。

本实施例中，两个所述锚杆振动测量阵列1-2中共包括M个所述振动检测单元8，M个所述振动检测单元8沿综采面7的推进方向由前至后分为m1个振动检测组，m1个所述振动检测组分别布设在m1个振动检测断面上，每个所述振动检测断面上均包含两个振动检测单元8，两个所述振动检测单元8包括一个布设在工作面运输巷4内的振动检测单元8和一个布设在工作面回风巷5内的振动检测单元8。如图5-1所示，位于综采面7前方的第一个振动检测断面上两个测点的布设位置示意图，即测点1与测点m1+1的布设位置示意图。

实际使用过程中，为操作简便，对工作面运输巷4内布设的多个振动检测单元8的布设位置，按照位置先后顺序由后向前进行编号，即对工作面运输巷4和工作面回风巷5内所有测点由前至后进行编号。本实施例中，工作面运输巷4内所布设振动检测单元8的数量为m1个，其中m1为正整数且m1≥3；工作面回风巷5内所布设振动检测单元8的数量为m2个，其中m2为正整数且m2≥3；相应地，工作面运输巷4内所布设测点的总数量为m1个，工作面回风巷5内所布设测点的总数量为m2个。本实施例中，m1＝m2，并且工作面运输巷4内的m1个所述测点与工作面回风巷5内的m2个所述测点呈对称布设，详见图5-2。所述工作面运输巷4内的m1个所述测点，沿工作面推进方向由后向前分别为测点1、测点2、…、测点m1，所述工作面回风巷5内的m2个所述测点，沿工作面推进方向由后向前分别为测点m1+1、测点m1+2、…、测点M。

本实施例中，M个所述测点的布设高度均相同

本实施例中，步骤B221中对M个测点中任意两个测点的应力波实际走时差进行计算时，采用互相关运算方法进行计算。由于两个测点的应力波为同一激励源，因而，能简便计算得出M个测点中任意两个测点的应力波实际走时差。

本实施例中，步骤B222中按照射线追踪方法对当前分析应力波穿越成像区域的M条射线的轨迹和旅行时进行求解时，采用常规的射线追踪方法进行求解。

并且，步骤B2中对提取出的m个所述有效振动监测位置的振动监测结果进行分析处理时，按照射线追踪方法，先建立反映所述速度模型中各速度节点与各条射线的关系的雅可比矩阵A_(M1,N1)，雅可比矩阵A_(M1,N1)为一个M1行N1列矩阵，其中M1为m个所述有效振动监测位置对应的射线总数量且M1＝m×M，N1为正整数且其为所述速度模型中网格的数量。雅可比矩阵A_(M1,N1)中第i1行第j1列的数据，记作A_(i1,j1)；A_(i1,j1)为第i1条射线在所述速度模型的第j1网格中的长度。实际使用时，当最短射线路径确定后，A_(i1,j1)即可确定，其中i1为正整数且i1＝1、2、…、M1，j1为正整数且j1＝1、2、…、N1。

其中，对M1条射线的旅行时进行求解时，待所建立的速度模型进行网格离散化后，结合预先建立的雅可比矩阵A_(M1,N1)，并根据每条射线经过所有网格的时间之和计算得出各条射线的旅行时(即旅行时理论值)。其中，第i1条射线的旅行时，根据公式进行计算，其中S_(j1)为修正前所建立的速度模型中第j1个网格的慢度。

其中，SIRT算法为同步迭代算法，也称为联合迭代重建法或瞬时迭代重建法。

本实施例中，步骤C中进行层析成像之前，先对采煤机3不同开采位置的振动监测结果进行矫正；对采煤机3任一开采位置的振动监测结果进行矫正时，对采煤机3在该开采位置进行采煤时两个所述锚杆振动测量阵列1-2中所有锚杆振动检测单元8检测到的振动信号分别进行矫正，所有锚杆振动检测单元8检测到的振动信号的矫正方法均相同且均按照如图1所示的基于白噪声统计特征的振动有效量提取方法进行矫正。

本实施例中，对采煤机3不同开采位置的振动监测结果进行矫正时，采用上位机2-1进行矫正。

实际使用时，也可以采用两个分别与两个所述锚杆振动测量阵列1-2的振动测量信号矫正单元进行矫正。所述振动测量信号矫正单元包括第三处理器，所述第三处理器与锚杆振动测量阵列1-2中各锚杆振动检测单元8的第一处理器8-3相接。所述振动测量信号矫正单元同样装在一个防爆壳体内。

因而，所述数据处理设备为上位机2-1或所述振动测量信号矫正单元。

实际使用时，只需在工作面运输巷和工作面回风巷内分别布设在一个锚杆振动测量阵列，并在振动激励源(即采煤机)上位置一个开采位置检测单元即可，以煤机采动作为振动激励源，无需在采煤工作面上开设炮孔并布放***，因而能大幅度减少炮孔的钻孔与装药所需的人力、物力与原料成本，并能节约施工成本和施工时间，省工省时，同时能避免因炮孔***带来的采煤工作面开采过程存在的安全隐患，并且不会对巷道支护体系造成不良影响，避免煤矿井下的环境污染问题，因而本发明以采煤机截割响应作为分析对象，并通过综采面两侧回采巷道上布设的振动检测单元进行振动监测，能有效适用于采动过程的煤岩动态成像，具有非常重要的经济效益和社会效益。同时，由于进行动态层析成像之前，均先按照如图1所示的基于白噪声统计特征的振动有效量提取方法进行矫正，因而能有效保证动态层析成像质量，成像效果好。

以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。

Claims (9)

1.一种基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：采用数据处理设备对布设于被测试锚杆(6)上的锚杆振动检测单元(8)所检测振动信号x(t)进行矫正，包括以下步骤：

步骤一、振动信号经验模态分解：调用经验模态分解模块，对锚杆振动检测单元(8)所检测振动信号x(t)进行经验模态分解，获得n个本征模态函数分量和一个趋势项；n个本征模态函数分量和一个趋势项均为离散函数；

其中，振动信号x(t)为锚杆(6)的振动位移或振动加速度随时间变化的函数；振动信号x(t)的经验模态分解结果为：式(1)中c_i(t)为振动信号x(t)经验模态分解后获得的第i个本征模态函数分量，i为正整数且其本征模态函数分量的编号，i＝1、2、…、n；r_n(t)为振动信号x(t)经验模态分解后获得的趋势项；

步骤二、本征模态函数分量有效性识别：调用本征模态函数分量有效性识别模块，且根据预先建立的白噪声近似熵分布区间模板，对步骤一中分解出的振动信号x(t)的n个本征模态函数分量分别进行有效性识别；n个本征模态函数分量的有效性识别方法均相同；

对所述白噪声近似熵分布区间模板进行建立时，采用所述数据处理设备进行建立，过程如下：

步骤Ⅰ、样本获取：产生多个随机白噪声信号作为样本，多个所述白噪声信号的数据长度均相同；

其中，产生的第j个随机白噪声信号，记作ζ_j(t)；j为随机白噪声信号的编号且其为正整数，j＝1、2、…、N；N为本步骤中所产生随机白噪声信号的总数量；

步骤Ⅱ、经验模态分解：调用所述经验模态分解模块，对步骤Ⅰ中各白噪声信号分别进行经验模态分解，获得各白噪声信号的n个本征模态函数分量和一个趋势项；n个本征模态函数分量和一个趋势项均为离散函数；

其中，ζ_j(t)的经验模态分解结果为：式(1)中c_i,j(t)为ζ_j(t)经验模态分解后获得的第i个本征模态函数分量，r_n,j(t)为ζ_j(t)经验模态分解后获得的趋势项；

步骤Ⅲ、近似熵与平均周期自然对数计算：调用近似熵计算模块，对步骤Ⅱ中各白噪声信号的n个本征模态函数分量的近似熵分别进行计算；同时，调用平均周期计算模块，对各白噪声信号的n个本征模态函数分量的平均周期分别进行计算，并求解出各白噪声信号的n个本征模态函数分量的平均周期对数；其中，平均周期对数为平均周期的自然对数；

步骤Ⅳ、二维联合概率密度函数确定：根据步骤Ⅲ中计算得出的各白噪声信号的n个本征模态函数分量的近似熵和平均周期对数，调用二维联合概率密度函数确定模块，获取白噪声信号的n个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数；

其中，根据步骤Ⅲ中计算得出的N个所述白噪声信号的第i个本征模态函数分量的近似熵和平均周期对数，且调用所述二维联合概率密度函数确定模块，得出白噪声信号的第i个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数，记作f(p_i,t_i)；其中，p_i和t_i分别表示白噪声信号的第i个本征模态函数分量的近似熵与平均周期对数；

步骤Ⅴ、白噪声近似熵分布区间模板获取：根据步骤Ⅳ中获得的白噪声信号的n个本征模态函数分量的近似熵与平均周期对数的二维联合概率密度函数，获得白噪声信号的n个本征模态函数分量的置信区间；所述置信区间的置信水平为c％，其中c％＝92％～99％；

其中，白噪声信号的n个本征模态函数分量的置信区间组成所述白噪声近似熵分布区间模板；

对步骤一中分解出的第i个本征模态函数分量c_i(t)进行有效性识别时，包括以下步骤：

步骤201、近似熵与平均周期自然对数计算：调用步骤Ⅲ中所述近似熵计算模块，计算得出c_i(t)的近似熵；同时，调用步骤Ⅲ中所述平均周期计算模块，计算得出c_i(t)的平均周期对数；

步骤202、有效性识别：根据步骤201中计算得出的c_i(t)的近似熵与平均周期自然对数，判断c_i(t)是否落入步骤Ⅴ中获得的白噪声信号的第i个本征模态函数分量的置信区间内：当落入步骤Ⅴ中获得的白噪声信号的第i个本征模态函数分量的置信区间内时，说明c_i(t)为有效本征模态函数分量；否则，说明c_i(t)为虚假本征模态函数分量；

步骤203、多次重复步骤201至步骤202，直至完成步骤一中分解出的n个本征模态函数分量的有效性识别过程，获得振动信号x(t)经验模态分解后的所有有效本征模态函数分量；

步骤三、信号重构：利用步骤203中获得的振动信号x(t)经验模态分解后的所有有效本征模态函数分量和步骤一中分解出的趋势项，重构出振动信号x'(t)；x'(t)为对振动信号x(t)进行矫正后获得的信号。

2.按照权利要求1所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：步骤一中进行振动信号经验模态分解之前，先调用归一化处理模块，对振动信号x(t)进行归一化处理；再调用所述经验模态分解模块，对归一化处理后的振动信号x(t)进行经验模态分解；

步骤Ⅱ中进行经验模态分解之前，先调用所述归一化处理模块，对步骤Ⅰ中各白噪声信号分别进行归一化处理；再调用所述经验模态分解模块，对归一化处理后的各白噪声信号分别进行经验模态分解。

3.按照权利要求1或2所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：步骤Ⅲ中调用所述近似熵计算模块对各白噪声信号的n个本征模态函数分量的近似熵进行计算时，所有本征模态函数分量的近似熵的计算方法均相同；

其中，c_i,j(t)的近似熵，记作IMF-ApEn_i,j；对c_i,j(t)的近似熵进行计算时，过程如下：

步骤Ⅲ-1、排列组成m维向量O(k)：按照c_i,j(t)中的各数据先后顺序，将c_i,j(t)中的m个数据排列组成M₀-m+1个m维矢量；

其中，c_i,j(t)中包含M₀个数据且M₀个数据分别为c_i,j(1)、c_i,j(2)、…、c_i,j(M₀)，M₀为c_i,j(t)的数据长度且其为正整数；m为预先设定的嵌入维数，m为正整数且m≥2；

M₀-m+1个m维矢量中第k个m维矢量，记作O(k)；

并且，O(k)＝[c_i,j(k),c_i,j(k+1),…,c_i,j(k+m-1)](2-0)，其中k为正整数且k＝1、2、…、M₀-m+1；

步骤Ⅲ-2、Hausdorff空间距离计算：调用Hausdorff空间距离计算模块，计算得出M₀-m+1个所述m维矢量中任一个m维矢量与其它M₀-m个m维矢量之间的Hausdorff空间距离；

其中，m维矢量O(k)与m维矢量O(g)之间的Hausdorff空间距离，记作d[O(k),O(g)]；m维矢量O(g)为步骤Ⅲ-1中M₀-m+1个m维矢量中第g个m维矢量，g为正整数且g＝1、2、…、M₀-m+1，g≠k；

步骤Ⅲ-3、m维矢量间相似度计算：根据预先限定的相似容限r，计算得出M₀-m+1个所述m维矢量中任一个m维矢量与其它M₀-m个m维矢量之间的相似度；

其中，m维矢量O(k)与其它M₀-m个m维矢量之间的相似度，记作 根据公式 $C_k^m\left(r\right)=\underset{g=1}{\overset{N-m+1}{\&Sigma;}}\mathrm{\&Theta;}\{r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\}/(M_0-m+1)---(2-2),$ 进行计算；r＞0；

式(2-2)中， $\mathrm{\&Theta;}\{r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\}=\left\{\begin{array}{cc}1,& r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;\&GreaterEqual;0\\ 0,& r-d\&lsqb;O\left(k\right),O\left(g\right)\&rsqb;<0\end{array}\right.---(2-3);$

步骤Ⅲ-4、Φ^m(r)计算：根据步骤Ⅲ-3中计算得出的M₀-m+1个所述m维矢量中任一个m维矢量与其它M₀-m个m维矢量之间的相似度，并按照公式 ${\mathrm{\&Phi;}}^m\left(r\right)=\frac{1}{M-m-1}\underset{k=1}{\overset{N-m+1}{\&Sigma;}}{\mathrm{lnC}}_k^m\left(r\right)---(2-4),$ 计算得出Φ^m(r)；

步骤Ⅲ-5、Φ^m+1(r)计算：按照步骤Ⅲ-1至步骤Ⅲ-4中所述的方法，计算得出Φ^m+1(r)；

步骤Ⅲ-6、近似熵计算：根据公式IMF-ApEn_i,j＝Φ^m(r)-Φ^m+1(r)(2-5)，计算得出IMF-ApEn_i,j。

4.按照权利要求3所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：所述的m＝2，r＝c₀×S_td[c_i,j(t)]；其中，c₀＝0.1～0.25，S_td[c_i,j(t)]为c_i,j(t)的标准差。

5.按照权利要求1或2所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：步骤Ⅲ中调用所述平均周期计算模块，对各白噪声信号的n个本征模态函数分量的平均周期进行计算时，所有本征模态函数分量的平均周期的计算方法均相同；

其中，对c_i,j(t)的平均周期进行计算时，采用快速傅里叶变换方法且按照公式进行计算；其中，F[c_i,j(t)]表示c_i,j(t)的快速傅里叶变换。

6.按照权利要求5所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：所述的

c_i,j(t)中包含M₀个数据，M₀为c_i,j(t)的数据长度且其为正整数；M₀个数据中第h个数据，记作c_i,j(h)；其中，h为正整数且h＝1、2、…、M₀。

7.按照权利要求1或2所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：步骤Ⅰ中多个所述随机白噪声信号的数据长度均与振动信号x(t)的数据长度相等。

8.按照权利要求1或2所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：步骤Ⅰ中所述的N≥200。

9.按照权利要求1或2所述的基于白噪声统计特征的锚杆振动信号矫正方法，其特征在于：所述锚杆(6)布设在煤矿开采过程中综采面(7)两侧的工作面运输巷(4)或工作面回风巷(5)内，所述锚杆振动检测单元(8)安装在锚杆(6)的外端。