一种基于奇异值分解中值法的瓦斯浓度数据降噪方法与流程

本发明属于信号处理技术领域，具体涉及一种基于奇异值分解中值法的瓦斯浓度数据降噪滤波方法。

背景技术：

由于煤矿井下环境十分恶劣，布置在井下的瓦斯传感器经常受到各种干扰的影响，如烟尘、高温、水气等，并且还会受到电磁干扰的影响，致使采集到的瓦斯浓度数据普遍含有噪声。如果用含有噪声的瓦斯浓度数据直接进行分析处理，不仅不能准确预测瓦斯涌出量，及时预警危险，而且浪费时间，做大量无用工作。因此，对瓦斯浓度数据必须进行去噪处理还原其真实的发展变化趋势。

目前对瓦斯信号降噪的主要方法有小波变换降噪方法和支持向量回归机降噪方法，由于采集到的瓦斯信号数据往往具有混沌特征，其频谱散布于整个频率空间，这时，采用小波变换很难将有用信号和噪声频谱严格区分开来并且小波阈值和支持向量回归机的核函数无法准确选择。

技术实现要素：

针对现有技术中存在的上述技术问题，本发明提出了一种基于奇异值分解中值法的瓦斯浓度数据降噪方法，设计合理，克服了现有技术的不足，具有良好的效果。

为了实现上述目的，本发明采用如下技术方案：

一种基于奇异值分解中值法的瓦斯浓度数据降噪方法，包括以下步骤：

步骤1：导入含噪瓦斯浓度监测数据{xt,t＝1,2,…,n}；

步骤2：检测含噪瓦斯浓度数据中是否含有单个异常数据和缺失数据；

若含有单个异常数据和缺失数据，则通过移动平均线法处理单个异常数据，通过三次指数平滑法处理缺失数据；若不含有单个异常数据和缺失数据，则无需处理；

步骤3：将步骤2中的瓦斯浓度数据构造成hankel矩阵并进行svd(singularvaluedecomposition，奇异值分解)变换；

步骤4：基于奇异值中值滤波策略选出有效奇异值，选取λ1,λ2,…,λr/2作为有效奇异值；

步骤5：通过svd逆变换、重建hankel矩阵并进行信号重构，得到降噪后的信号。

优选地，在步骤2中，当有数值满足式(1)时，表示数据中含有单个异常数据；当有数值满足式(2)时，表示数据中含有缺失数据；

||xt-xt-1|-|xt-xt+1||>0.02％(1)；

式中：xt表示当前采样点的瓦斯浓度；xt-1表示当前采样点前一个点的瓦斯浓度；xt+1表示当前采样点后一个点的瓦斯浓度；0.02为判定是否有单个异常数据的阈值；％为导入的含噪瓦斯浓度监测数据单位；

xt＝null||xt＝？||xt＝*(2)；

式中：null表示为空；？和*表示特殊符号；

对单个异常数据使用移动平均线法进行处理：

若采样点t＝b时满足式(1)即出现单个瓦斯数据异常高或异常低的现象，则通过式(3)计算移动平均数值xb，单个异常数据用xb表示；

式中：k表示b点之前的瓦斯浓度数据的采样点；

对缺失数据使用三次指数平滑法进行处理：

监测数据{x(t),t＝1,2,…,n}中间含有缺失数据，先确定插入数据点数与平滑处理的步距l，以缺失数据之前的瓦斯浓度数据为基础数据，按照式(4)-式(8)进行平滑处理：

xt+l＝at+btl+ctl²/2(4)；

式中：xt+l为l步平滑值即得到的缺失数据；at,bt,ct为三次指数平滑法的预测参数；

at,bt,ct的计算公式如下：

at＝3s″t-3s″t+s″″t(5)；

式中：s′t为t点的一次指数平滑值；s″t为t点的二次指数平滑值；s″′t为t点的三次指数平滑值；α为权数；

s′t,s″t,s″′t平滑值的计算公式如下：

式中：xt为缺失数据之前的瓦斯浓度数据；s′t-1为t-1点的一次指数平滑值；s″t-1为t-1点的二次指数平滑值；s″′t-1为t-1点的三次指数平滑值。

优选地，取α＝0.5。

优选地，在步骤3中，基于相空间重构理论，对含噪瓦斯浓度数据构造如式(9)所示的p*q阶hankel矩阵：

式中：hpq为p*q阶矩阵，其中n为信号长度，n＝p+q-1并且p≥q；

对hpq进行svd变换如式(10)所示：

hpq＝u∑v^t(10)；

式中：u表示p*p阶的左奇异矩阵，v^t表示q*q阶的右奇异矩阵，∑表示p*q阶的对角奇异矩阵，其表达式如式(11)所示：

式中：λ1,λ2,…,λr为矩阵hpq的奇异值，且λ1≥λ2≥…≥λr≥0；

公式(10)的具体推导展开过程如下：

根据矩阵hpq在重构空间的特性将矩阵hpq表示为h＝d+w的形式；

其中，d表示纯净信号的p*q矩阵，w表示噪声干扰信号的p*q矩阵；

去噪的理想目标就是从矩阵hpq中恢复出d包含的信号即通过svd分解从矩阵hpq中恢复信号子空间；

假设d存在秩亏，即rank(d)＝r(r<q)，且具有如下svd分解：

式中：ux1为p*r阶矩阵，ux2为p*(p-r)阶矩阵，∑x1为r*r阶矩阵，vx1为r*q阶矩阵，vx2为(p-r)*q阶矩阵，r为矩阵hpq的秩，ux1张成的空间为d的列空间，即称为信号子空间；

根据将前面的带噪信号矩阵hpq重写如公式(14)所示：

式中：理想目标p1＝ux1，直接恢复d的信号子空间，但由于p1≠ux1，无法直接恢复d的信号子空间，就要寻找d的最佳逼近矩阵即选取有效的奇异值。

优选地，在步骤4中，对于无噪信号，对角矩阵s为满秩，即所有奇异值都是有效的；对于含噪信号，根据公式(12)以及奇异值分解理论和frobeious范数意义下矩阵最佳逼近定理得到：有效的信号包含在较大的奇异值中，噪声信号包含在较小的奇异值中，并且奇异值下降迅速，中值之前奇异值(λ1,λ2,…,λr/2)的和就占了全部奇异值(λ1,λ2,…,λr)之和的99％以上的比例，即

优选地，在步骤5中，保留大于中值的奇异值(λ1,λ2,…,λr/2)，将小于中值的奇异值设置为零，则源信号中的噪声被去除，即

再进行svd的逆变换即

得到矩阵相对于矩阵hpq而言少了一半的奇异值，不符合hankel矩阵的形式，将矩阵中的反对角线元素采用式(18)进行平均：

式中：i为矩阵hpq的行，j为矩阵hpq的列，m＝max(1，i-p+1)，n＝min(q，i)；由构成的即为降噪后的瓦斯信号。

本发明原理如下：

为了实现瓦斯浓度数据的有效降噪滤波，本发明针对瓦斯信号的随机性、非平稳性和混沌性的特点，对导入的含噪瓦斯信号数据进行检测，判断是否含有单个异常数据和缺失数据，然后构造hankel矩阵并进行svd变换，基于中值滤波策略选出有效奇异值，svd逆变换重建hankel矩阵即可得到降噪后的瓦斯信号。

本发明所带来的有益技术效果：

本发明基于中值滤波策略选出有效奇异值进行信号重构，在充分保留瓦斯信号随机性、非平稳特征的基础上，对瓦斯信号进行有效降噪滤波，能对瓦斯信号进行滤波，相比于小波变换和支持向量回归机降噪方法，有效解决了瓦斯信号数据的混沌特性，具有良好的稳定性、不变性和噪声鲁棒性，能够反映数据的内在属性，同时又能更好地保留信号细节；该方法简单易行、效果较为理想，具有很好的技术价值和应用前景。

附图说明

图1为本发明一种基于奇异值分解中值法的瓦斯浓度数据降噪方法的流程图。

图2为含噪瓦斯信号中含有单个异常数据的示意图。

图3为含噪瓦斯信号中缺失数据的示意图。

图4为对单个异常数据和缺失数据处理后的含噪瓦斯信号的示意图。

图5为使用本发明降噪后的瓦斯信号示意图。

图6为瓦斯信号数据所含噪声的示意图。

图7为几种降噪方法的对比图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

如图1所示，一种基于奇异值分解中值法的瓦斯浓度数据降噪方法，包括以下步骤：

步骤1：导入含噪瓦斯浓度监测数据{xt,t＝1,2,…,n}；

步骤2：检测含噪瓦斯浓度数据中是否含有单个异常数据和缺失数据；

若含有单个异常数据和缺失数据，则通过移动平均线法处理单个异常数据，通过三次指数平滑法处理缺失数据；若不含有单个异常数据和缺失数据，则无需处理；

当有数值满足式(1)时，表示数据中含有单个异常数据，如图2所示；当有数值满足式(2)时，表示数据中含有缺失数据，如图3所示；

||xt-xt-1|-|xt-xt+1||>0.02％(1)；

式中：xt表示当前采样点的瓦斯浓度；xt-1表示当前采样点前一个点的瓦斯浓度；xt+1表示当前采样点后一个点的瓦斯浓度；0.02为判定是否有单个异常数据的阈值；％为导入的含噪瓦斯浓度监测数据单位；

xt＝null||xt＝？||xt＝*(2)；

式中：null表示为空；？和*表示特殊符号；

对单个异常数据使用移动平均线法进行处理：

若采样点t＝b时满足式(1)即出现单个瓦斯数据异常高或异常低的现象，则通过式(3)计算移动平均数值xb，单个异常数据用xb表示；

式中：k表示b点之前的瓦斯浓度数据的采样点；

对缺失数据使用三次指数平滑法进行处理：

监测数据{xt,t＝1,2,…,n}中间含有缺失数据，先确定插入数据点数与平滑处理的步距l，以缺失数据之前的瓦斯浓度数据为基础数据，按照式(4)-式(8)进行平滑处理：

xt+l＝at+btl+ctl²/2(4)；

式中：xt+l为l步平滑值即得到的缺失数据；at,bt,ct为三次指数平滑法的预测参数。

at,bt,ct的计算公式如下：

at＝3s′t-3s″t+s″′t(5)；

式中：s′t为t点的一次指数平滑值；s″t为t点的二次指数平滑值；s″′t为t点的三次指数平滑值；α为权数取0.5；

s′t,s″t,s″′t平滑值得计算公式如下：

式中：xt为缺失数据之前的瓦斯浓度数据；s′t-1为t-1点的一次指数平滑值；s″t-1为t-1点的二次指数平滑值；s″′t-1为t-1点的三次指数平滑值。

如图4所示，单个异常数据和缺失数据处理之后的含噪瓦斯信号数据。

步骤3：将步骤2中的瓦斯浓度数据构造成hankel矩阵并进行svd变换；

基于相空间重构理论，对含噪瓦斯浓度数据构造如式(9)所示的p*q阶hankel矩阵(吸引子轨迹矩阵)：

式中：hpq为p*q阶矩阵，其中n为信号长度，n＝p+q-1并且p≥q；

对hpq进行svd变换如式(10)所示：

hpq＝u∑v^t(10)；

式中：u表示p*p阶的左奇异矩阵，v^t表示q*q阶的右奇异矩阵，∑表示p*q阶的对角奇异矩阵如式(11)、(12)所示：

式中：λ1,λ2,…,λr为矩阵hpq的奇异值，且λ1≥λ2≥…≥λr≥0；

公式(10)的具体推导展开过程如下：

根据矩阵hpq在重构空间的特性将矩阵hpq表示为h＝d+w的形式；

其中，d表示纯净信号的p*q矩阵，w表示噪声干扰信号的p*q矩阵；

去噪的理想目标就是从矩阵hpq中恢复出d包含的信号即通过奇异值分解从矩阵hpq中恢复信号子空间；

假设d存在秩亏，即rank(d)＝r(r<q)，且具有如下svd分解：

式中：ux1为p*r阶矩阵，ux2为p*(p-r)阶矩阵，∑x1为r*r阶矩阵，vx1为r*q阶矩阵，vx2为(p-r)*q阶矩阵，r为矩阵hpq的秩，ux1张成的空间为d的列空间，即称为信号子空间；

根据将前面的带噪信号矩阵hpq重写如公式(14)所示：

式中：和分别为其上一行括号中矩阵的svd，即和理想目标p1＝ux1，直接恢复d的信号子空间，但由于p1≠ux1，无法直接恢复d的信号子空间，就要寻找d的最佳逼近矩阵，就是如何选取有效的奇异值。

步骤4：基于奇异值中值滤波策略选出有效奇异值；

对于无噪信号，对角矩阵s为满秩，即所有奇异值都是有效的；对于含噪信号，根据公式(12)以及奇异值分解理论和frobeious范数意义下矩阵最佳逼近定理得到：有效的信号包含在较大的奇异值中，噪声信号包含在较小的奇异值中，并且奇异值下降迅速，中值之前奇异值(λ1,λ2,…,λr/2)的和就占了全部奇异值(λ1,λ2,…,λr)之和的99％以上的比例，即

所以本发明基于奇异值中值滤波策略来选取有效奇异值，选取λ1,λ2,…,λr/2作为有效奇异值。

步骤5：通过svd逆变换、重建hankel矩阵并进行信号重构，得到降噪后的信号。

保留大于中值的奇异值(λ1,λ2,…,λr/2)，将小于中值的奇异值设置为零，则源信号中的噪声被去除，即

再进行svd的逆变换即

得到矩阵相对于矩阵hpq而言少了一半的奇异值，不符合hankel矩阵的形式，将矩阵中的反对角线元素采用式(18)进行平均：

式中：i为矩阵hpq的行，j为矩阵hpq的列，m＝max(1，i-p+1)，n＝min(q，i)；由构成的即为降噪后的瓦斯信号。

如图5所示即为降噪后的瓦斯信号，图6为含噪瓦斯信号数据中所含的噪声示意图。

实验验证及分析

将svd中值法与小波阈值法、支持向量回归机降噪法和svd均值法进行对比。

从衡量指标信噪比(snr)和均方误差(mse)上对比以上方法性能的好坏，经过实验研究表明，如图7所示，本发明对含噪瓦斯数据降噪效果较好。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。