CN107704831A - 一种基于奇异值分解中值法的瓦斯浓度数据降噪方法 - Google Patents
一种基于奇异值分解中值法的瓦斯浓度数据降噪方法 Download PDFInfo
- Publication number
- CN107704831A CN107704831A CN201710939526.XA CN201710939526A CN107704831A CN 107704831 A CN107704831 A CN 107704831A CN 201710939526 A CN201710939526 A CN 201710939526A CN 107704831 A CN107704831 A CN 107704831A
- Authority
- CN
- China
- Prior art keywords
- matrix
- data
- formula
- gas concentration
- signal
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
Abstract
本发明公开了一种基于奇异值分解中值法的瓦斯浓度数据降噪方法，属于信号处理技术领域，包括如下步骤：导入含噪瓦斯浓度数据X；检测含噪瓦斯浓度数据中是否含有单个异常数据和缺失数据，若有则使用移动平均线法处理单个异常数据，使用三次指数平滑法处理缺失数据，若没有则无需处理；对含噪瓦斯浓度数据构造Hankel矩阵；对Hankel矩阵进行奇异值分解(Singular Value Decomposition，简称svd)变换；基于奇异值中值滤波策略选出有效奇异值；svd逆变换重建Hankel矩阵，得到降噪后的瓦斯浓度数据。本发明提出的瓦斯信号降噪方法，通过实测瓦斯信号进行了降噪实验，结果表明该方法具有较好的降噪性能，能有效提高瓦斯信号分析精度。
Description
技术领域
本发明属于信号处理技术领域，具体涉及一种基于奇异值分解中值法的瓦斯浓度数据降噪滤波方法。
背景技术
由于煤矿井下环境十分恶劣，布置在井下的瓦斯传感器经常受到各种干扰的影响，如烟尘、高温、水气等，并且还会受到电磁干扰的影响，致使采集到的瓦斯浓度数据普遍含有噪声。如果用含有噪声的瓦斯浓度数据直接进行分析处理，不仅不能准确预测瓦斯涌出量，及时预警危险，而且浪费时间，做大量无用工作。因此，对瓦斯浓度数据必须进行去噪处理还原其真实的发展变化趋势。
目前对瓦斯信号降噪的主要方法有小波变换降噪方法和支持向量回归机降噪方法，由于采集到的瓦斯信号数据往往具有混沌特征，其频谱散布于整个频率空间，这时，采用小波变换很难将有用信号和噪声频谱严格区分开来并且小波阈值和支持向量回归机的核函数无法准确选择。
发明内容
针对现有技术中存在的上述技术问题，本发明提出了一种基于奇异值分解中值法的瓦斯浓度数据降噪方法，设计合理，克服了现有技术的不足，具有良好的效果。
为了实现上述目的，本发明采用如下技术方案：
一种基于奇异值分解中值法的瓦斯浓度数据降噪方法，包括以下步骤：
步骤1：导入含噪瓦斯浓度监测数据{Xt,t＝1,2,…,N}；
步骤2：检测含噪瓦斯浓度数据中是否含有单个异常数据和缺失数据；
若含有单个异常数据和缺失数据，则通过移动平均线法处理单个异常数据，通过三次指数平滑法处理缺失数据；若不含有单个异常数据和缺失数据，则无需处理；
步骤3：将步骤2中的瓦斯浓度数据构造成Hankel矩阵并进行svd(Singular ValueDecomposition，奇异值分解)变换；
步骤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+ctL2/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∑VT (10)；
式中：U表示p*p阶的左奇异矩阵，VT表示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+ctL2/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∑VT (10)；
式中：U表示p*p阶的左奇异矩阵，VT表示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所示，本发明对含噪瓦斯数据降噪效果较好。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。
Claims (6)
1.一种基于奇异值分解中值法的瓦斯浓度数据降噪方法，其特征在于，包括以下步骤：
步骤1：导入含噪瓦斯浓度监测数据{Xt,t＝1,2,…,N}；
步骤2：检测含噪瓦斯浓度数据中是否含有单个异常数据和缺失数据；
若含有单个异常数据和缺失数据，则通过移动平均线法处理单个异常数据，通过三次指数平滑法处理缺失数据；若不含有单个异常数据和缺失数据，则无需处理；
步骤3：将步骤2中的瓦斯浓度数据构造成Hankel矩阵并进行svd变换；
步骤4：基于奇异值中值滤波策略选出有效奇异值，选取λ1,λ2,…,λr/2作为有效奇异值；
步骤5：通过svd逆变换、重建Hankel矩阵并进行信号重构，得到降噪后的信号。
2.根据权利要求1所述基于奇异值分解中值法的瓦斯浓度数据降噪方法，其特征在于，在步骤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表示；
<mrow>
<msub>
<mi>x</mi>
<mi>b</mi>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<msubsup>
<mo>&Sigma;</mo>
<mrow>
<mi>t</mi>
<mo>=</mo>
<mi>k</mi>
</mrow>
<mrow>
<mi>b</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<msub>
<mi>x</mi>
<mi>t</mi>
</msub>
</mrow>
<mrow>
<mi>b</mi>
<mo>-</mo>
<mi>K</mi>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
式中：K表示b点之前的瓦斯浓度数据的采样点；
对缺失数据使用三次指数平滑法进行处理：
监测数据{X(t),t＝1,2,…,N}中间含有缺失数据，先确定插入数据点数与平滑处理的步距L，以缺失数据之前的瓦斯浓度数据为基础数据，按照式(4)-式(8)进行平滑处理：
xt+L＝at+btL+ctL2/2 (4)；
式中：Xt+L为L步平滑值即得到的缺失数据；at,bt,ct为三次指数平滑法的预测参数；
at,bt,ct的计算公式如下：
at＝3St′-3St″+St″′ (5)；
<mrow>
<msub>
<mi>b</mi>
<mi>t</mi>
</msub>
<mo>=</mo>
<mfrac>
<mi>&alpha;</mi>
<mrow>
<mn>2</mn>
<msup>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&lsqb;</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>-</mo>
<mn>5</mn>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mo>&prime;</mo>
</msubsup>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>10</mn>
<mo>-</mo>
<mn>8</mn>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
<mo>+</mo>
<mrow>
<mo>(</mo>
<mn>4</mn>
<mo>-</mo>
<mn>3</mn>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
<mrow>
<msub>
<mi>c</mi>
<mi>t</mi>
</msub>
<mo>=</mo>
<mfrac>
<msup>
<mi>&alpha;</mi>
<mn>2</mn>
</msup>
<msup>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mfrac>
<mo>&lsqb;</mo>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mo>&prime;</mo>
</msubsup>
<mo>-</mo>
<mn>2</mn>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>7</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
式中：St′为t点的一次指数平滑值；St″为t点的二次指数平滑值；St″′为t点的三次指数平滑值；α为权数；
St′,St″,St″′平滑值的计算公式如下：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mo>&prime;</mo>
</msubsup>
<mo>=</mo>
<msub>
<mi>&alpha;x</mi>
<mi>t</mi>
</msub>
<mo>+</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>S</mi>
<mrow>
<mi>t</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mo>&prime;</mo>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>&alpha;S</mi>
<mi>t</mi>
<mo>&prime;</mo>
</msubsup>
<mo>+</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>S</mi>
<mrow>
<mi>t</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>S</mi>
<mi>t</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>&alpha;S</mi>
<mi>t</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
<mo>+</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>&alpha;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>S</mi>
<mrow>
<mi>t</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msubsup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>8</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
式中：xt为缺失数据之前的瓦斯浓度数据；S′t-1为t-1点的一次指数平滑值；S″t-1为t-1点的二次指数平滑值；S″′t-1为t-1点的三次指数平滑值。
3.根据权利要求2所述的基于奇异值分解中值法的瓦斯浓度数据降噪方法，其特征在于，取α＝0.5。
4.根据权利要求1所述基于奇异值分解中值法的瓦斯浓度数据降噪方法，其特征在于，在步骤3中，基于相空间重构理论，对含噪瓦斯浓度数据构造如式(9)所示的p*q阶Hankel矩阵：
式中：Hpq为p*q阶矩阵，其中N为信号长度，N＝p+q-1并且p≥q；
对Hpq进行svd变换如式(10)所示：
Hpq＝U∑VT (10)；
式中：U表示p*p阶的左奇异矩阵，VT表示q*q阶的右奇异矩阵，∑表示p*q阶的对角奇异矩阵，其表达式如式(11)所示：
<mrow>
<mo>&Sigma;</mo>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>S</mi>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>11</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
式中：λ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分解：
<mrow>
<mi>D</mi>
<mo>=</mo>
<msub>
<mi>U</mi>
<mi>x</mi>
</msub>
<msub>
<mi>&Sigma;</mi>
<mi>x</mi>
</msub>
<msub>
<mi>V</mi>
<mi>x</mi>
</msub>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
</msub>
</mtd>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>&Sigma;</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>13</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
式中：Ux1为p*r阶矩阵，Ux2为p*(p-r)阶矩阵，∑x1为r*r阶矩阵，Vx1为r*q阶矩阵，Vx2为(p-r)*q阶矩阵，r为矩阵Hpq的秩，Ux1张成的空间为D的列空间，即称为信号子空间；
根据将前面的带噪信号矩阵Hpq重写如公式(14)所示：
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>H</mi>
<mo>=</mo>
<mi>D</mi>
<mo>+</mo>
<mi>W</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mi>D</mi>
<mo>+</mo>
<mi>W</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
</msub>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
<mi>T</mi>
</msubsup>
<mo>+</mo>
<msub>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
</msub>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
<mi>T</mi>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>DV</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
</msub>
<mo>+</mo>
<msub>
<mi>WV</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
<mi>T</mi>
</msubsup>
<mo>+</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>WV</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
<mi>T</mi>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mn>1</mn>
</msub>
<msub>
<mi>S</mi>
<mn>1</mn>
</msub>
<msubsup>
<mi>Q</mi>
<mn>1</mn>
<mi>T</mi>
</msubsup>
<mo>)</mo>
</mrow>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
<mi>T</mi>
</msubsup>
<mo>+</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mn>2</mn>
</msub>
<msub>
<mi>S</mi>
<mn>2</mn>
</msub>
<msubsup>
<mi>Q</mi>
<mn>2</mn>
<mi>T</mi>
</msubsup>
<mo>)</mo>
</mrow>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
<mi>T</mi>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<msub>
<mi>P</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>P</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<msub>
<mi>S</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mi>S</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>Q</mi>
<mn>1</mn>
<mi>T</mi>
</msubsup>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>1</mn>
</mrow>
<mi>T</mi>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>Q</mi>
<mn>2</mn>
<mi>T</mi>
</msubsup>
<msubsup>
<mi>V</mi>
<mrow>
<mi>x</mi>
<mn>2</mn>
</mrow>
<mi>T</mi>
</msubsup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>14</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
式中：理想目标P1＝Ux1，直接恢复D的信号子空间，但由于P1≠Ux1，无法直接恢复D的信号子空间，就要寻找D的最佳逼近矩阵即选取有效的奇异值。
5.根据权利要求1所述的基于奇异值分解中值法的瓦斯浓度数据降噪方法，其特征在于，在步骤4中，对于无噪信号，对角矩阵S为满秩，即所有奇异值都是有效的；对于含噪信号，根据公式(12)以及奇异值分解理论和Frobeious范数意义下矩阵最佳逼近定理得到：有效的信号包含在较大的奇异值中，噪声信号包含在较小的奇异值中，并且奇异值下降迅速，中值之前奇异值(λ1,λ2,…,λr/2)的和就占了全部奇异值(λ1,λ2,…,λr)之和的99％以上的比例，即
<mrow>
<msubsup>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>r</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msubsup>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>&ap;</mo>
<msubsup>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</msubsup>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>15</mn>
<mo>)</mo>
</mrow>
<mo>.</mo>
</mrow>
6.根据权利要求1所述基于奇异值分解中值法的瓦斯浓度数据降噪方法，其特征在于，在步骤5中，保留大于中值的奇异值(λ1,λ2,…,λr/2)，将小于中值的奇异值设置为零，则源信号中的噪声被去除，即
再进行svd的逆变换即
<mrow>
<mover>
<mi>H</mi>
<mo>^</mo>
</mover>
<mo>=</mo>
<msubsup>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>r</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msubsup>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<msub>
<mi>u</mi>
<mi>k</mi>
</msub>
<msubsup>
<mi>v</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>17</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
得到矩阵 相对于矩阵Hpq而言少了一半的奇异值，不符合Hankel矩阵的形式，将矩阵中的反对角线元素采用式(18)进行平均：
<mrow>
<msub>
<mover>
<mi>x</mi>
<mo>&OverBar;</mo>
</mover>
<mi>i</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mi>m</mi>
<mo>-</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</mfrac>
<msubsup>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</msubsup>
<msub>
<mi>H</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mi>j</mi>
<mo>+</mo>
<mn>1</mn>
<mo>,</mo>
<mi>j</mi>
</mrow>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>18</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
式中：i为矩阵Hpq的行，j为矩阵Hpq的列，m＝max(1，i-p+1)，n＝min(q，i)；由构成的即为降噪后的瓦斯信号。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710939526.XA CN107704831A (zh) | 2017-10-11 | 2017-10-11 | 一种基于奇异值分解中值法的瓦斯浓度数据降噪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710939526.XA CN107704831A (zh) | 2017-10-11 | 2017-10-11 | 一种基于奇异值分解中值法的瓦斯浓度数据降噪方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107704831A true CN107704831A (zh) | 2018-02-16 |
Family
ID=61183542
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710939526.XA Pending CN107704831A (zh) | 2017-10-11 | 2017-10-11 | 一种基于奇异值分解中值法的瓦斯浓度数据降噪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107704831A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108399385A (zh) * | 2018-02-23 | 2018-08-14 | 中国石油大学(华东) | 一种风力发电机组振动监测信号降噪方法 |
| CN108645920A (zh) * | 2018-04-09 | 2018-10-12 | 华南理工大学 | 一种基于去噪和对齐的钢轨超声探伤的直达波抑制方法 |
| US11431976B2 (en) | 2019-01-28 | 2022-08-30 | Kla Corporation | System and method for inspection using tensor decomposition and singular value decomposition |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103810394A (zh) * | 2014-02-28 | 2014-05-21 | 东北电力大学 | 旋转设备故障信号奇异值分解降噪的设计方法 |
| CN105761223A (zh) * | 2016-02-16 | 2016-07-13 | 四川用联信息技术有限公司 | 一种基于图像低秩性的迭代降噪方法 |
| CN106446829A (zh) * | 2016-09-22 | 2017-02-22 | 三峡大学 | 一种基于svd与vmd模态自相关分析的水电机组振动信号降噪方法 |
-
2017
- 2017-10-11 CN CN201710939526.XA patent/CN107704831A/zh active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103810394A (zh) * | 2014-02-28 | 2014-05-21 | 东北电力大学 | 旋转设备故障信号奇异值分解降噪的设计方法 |
| CN105761223A (zh) * | 2016-02-16 | 2016-07-13 | 四川用联信息技术有限公司 | 一种基于图像低秩性的迭代降噪方法 |
| CN106446829A (zh) * | 2016-09-22 | 2017-02-22 | 三峡大学 | 一种基于svd与vmd模态自相关分析的水电机组振动信号降噪方法 |
Non-Patent Citations (1)
| Title |
|---|
| 董丁稳: "基于安全监控系统实测数据的瓦斯浓度预测预警研究", 《中国博士学位论文全文数据库_工程科技Ⅰ辑》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108399385A (zh) * | 2018-02-23 | 2018-08-14 | 中国石油大学(华东) | 一种风力发电机组振动监测信号降噪方法 |
| CN108399385B (zh) * | 2018-02-23 | 2021-10-15 | 中国石油大学(华东) | 一种风力发电机组振动监测信号降噪方法 |
| CN108645920A (zh) * | 2018-04-09 | 2018-10-12 | 华南理工大学 | 一种基于去噪和对齐的钢轨超声探伤的直达波抑制方法 |
| CN108645920B (zh) * | 2018-04-09 | 2020-12-22 | 华南理工大学 | 一种基于去噪和对齐的钢轨超声探伤的直达波抑制方法 |
| US11431976B2 (en) | 2019-01-28 | 2022-08-30 | Kla Corporation | System and method for inspection using tensor decomposition and singular value decomposition |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gao et al. | Denoising nonlinear time series by adaptive filtering and wavelet shrinkage: a comparison | |
| CN107704831A (zh) | 一种基于奇异值分解中值法的瓦斯浓度数据降噪方法 | |
| CN103576060A (zh) | 基于小波自适应阈值的局部放电信号去噪方法 | |
| CN111679328A (zh) | 一种基于变分模态分解的瞬变电磁探测信号的降噪方法 | |
| Bing et al. | DeepCEDNet: an efficient deep convolutional encoder-decoder networks for ECG signal enhancement | |
| CN116088320B (zh) | 稀疏双时空非凸罚自适应Chirp模态交叉混叠分解方法 | |
| CN113887398A (zh) | 一种基于变分模态分解和奇异谱分析的gpr信号去噪方法 | |
| CN108847828A (zh) | 一种具有随机建模误差的非线性事件触发滤波方法 | |
| CN113568058A (zh) | 一种基于多分辨率奇异值分解的大地电磁信噪分离方法及系统 | |
| CN103839239A (zh) | 一种电缆瓷套终端红外图像自适应去噪方法 | |
| CN116127287A (zh) | 一种电阻率法勘探信号的降噪方法 | |
| CN112817056B (zh) | 一种大地电磁信号去噪方法及系统 | |
| CN112327260B (zh) | 一种sar回波数据中脉冲式干扰信号的抑制方法和装置 | |
| CN113009564B (zh) | 地震数据处理方法和装置 | |
| CN106845448B (zh) | 一种基于非负约束2d变分模态分解的红外弱小目标检测方法 | |
| CN115034978A (zh) | 基于子空间表示结合图嵌入约束的高光谱图像去噪方法 | |
| CN110703089B (zh) | 一种用于低频振荡Prony分析的小波阈值去噪方法 | |
| CN114355440A (zh) | 一种时频峰值滤波微震数据随机噪声的压制方法 | |
| Thoufeer et al. | A Deep Learning Scheme for Efficient Power System Waveform De-noising and Reconstruction | |
| CN112230200A (zh) | 一种基于激光雷达回波信号的改进型组合降噪方法 | |
| CN114325846B (zh) | 一种利用时间相干性抑噪的磁异常检测方法 | |
| Ting et al. | Eeg signal processing based on proper orthogonal decomposition | |
| Su et al. | A nonparametric derivative-based method for R wave detection in ECG | |
| CN113887362B (zh) | 一种局部放电信号的特征提取方法 | |
| Zhang et al. | Vibration signal filtering algorithm based on singular value subspace decomposition |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20180216 |