CN110988232A - 一种色谱基线降噪方法及其装置 - Google Patents

Abstract

本发明公开了一种色谱基线降噪方法，包括以下步骤：S1信号采集、S2信号放大、S3高斯噪声去除、S4脉冲噪声去除、S5标定拟合以及S6成分分析。本发明利用SG滤波器串联有形态学滤波器，综合两者的优点，摒除两者的缺点，先用SG最小二乘滤波算法去除高斯噪声，再用形态学滤波算法去除脉冲噪声。这样做的好处是，先用SG最小二乘滤波器去除大部分噪声，免除形态学滤波器对噪声的放大作用；再用形态学滤波器，可以滤除SG滤波器残留的异常点干扰问题。两个步骤顺序不可颠倒，否则会放大两种算法的缺点。最后，本发明中的两种滤波算法都是实时算法，可以实时对采集的信号进行滤波处理，延迟在0.5s以内。

Description

一种色谱基线降噪方法及其装置

技术领域

本发明涉及色谱基线降噪技术领域，具体为一种色谱基线降噪方法及其装置。

背景技术

随着我国工业化和城市化的快速发展以及能源消费的持续增长，以PM2.5 为特征的区域性复合型大气污染日益突出，区域内空气重污染现象大范围同时出现的频次日益增多，严重制约社会经济的可持续发展，威胁人民群众身体健康。而挥发性有机物(VOCs)是形成细颗粒物(PM2.5)、臭氧(O3)等二次污染物的重要前体物，会引发灰霾、光化学烟雾等大气环境问题，已成为现阶段我国大气环境领域的重点监测对象。

目前存在很多测定VOCs的技术，其中应用最广泛的是气相色谱技术(GC)。气相色谱法是利用气体作流动相的色层分离分析方法。但是在分离气体的过程中不免产生噪声污染，进而使得色谱图曲线不平滑，有毛刺，更严重的会影响对色谱成分的分析。因此，对色谱进行基线降噪非常必要。然而，现有的滤波方法都有其特定的适用场景及缺点：中值滤波，可以去除脉冲噪声，但是滤波窗口选择不当会造成失真；滑动平均滤波，可以去除高频噪声，但是会削平顶峰；高斯滤波，可以去除高斯白噪声，但是有明显的脉冲噪声残留，信号边缘受到严重模糊。并且，现实中的信号包含多种噪声，单一的方法必然会顾此失彼。

发明内容

本发明的目的在于提供一种色谱基线降噪方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明提供如下技术方案：一种色谱基线降噪方法，包括以下步骤：

S1：信号采集：通过采集设备，采集信号；

S2：信号放大：通过信号放大器，将步骤S1中采集的信号进行放大，并将放大后的信号通过USB传给PC端；

S3：高斯噪声去除：通过PC端，利用SG滤波器，以最小二乘滤波算法去除高斯噪声；

S4：脉冲噪声去除：利用形态学滤波器，对步骤S3去除高斯噪声的信号，以形态学滤波算法去除脉冲噪声；

S5：标定拟合：向采集设备中先通入已知浓度及成分的标准气体，通过观察谱图曲线，对每个峰的成分及浓度进行标定，对每种物质的峰高和浓度关系进行拟合，获得拟合曲线；

S6：成分分析：对步骤S4中去除脉冲噪音的信号，结合步骤S5中的拟合曲线，进行分析对比，得出结果。

优选的，所述采集设备为气相色谱仪。

优选的，所述标定拟合中拟合三次。

优选的，所述成分分析采用MATLAB自带的findpeaks寻峰方法，对色谱峰进行分析，进而获取峰高度、宽度、位置等信息。

其装置，包括气相色谱仪以及与气相色谱仪连接的放大器，所述放大器通过USB数据线连接有PC端，且PC端连接有SG滤波器，所述SG滤波器串联有形态学滤波器。

与现有技术相比，本发明的有益效果是：

1.保持峰高和峰型不变。因为峰高可以近似代表物质浓度，而峰型决定面积，面积也代表浓度。所以这个优点可以保证浓度测量的准确性。

2.计算简单、延迟小。SG滤波算法的时间复杂度为O(m*n)，其中m是滤波器长度，n是信号长度；形态学滤波算法的时间复杂度为O(m*n)，其中m， n同SG滤波器。且由于滤波过程是实时进行的，滤波窗口的大小是有限的，在二位数以内。所以整个算法的时间复杂度在常数级别，计算量小，从而计算速度快，延迟小。

3.可以实时滤波。在谱图曲线采集过程中及时滤波，可以提高用户体验。用户可以看到平滑的数据曲线，可以更直观的看到数据信息。

4.提高装置检测精度。曲线越平滑，越容易检测到低浓度的峰，最低检测浓度达0.005ppm。

附图说明

图1本发明的色谱基线降噪方法结构图；

图2本发明的简化版实时形态学滤波器的结构细节图；

图3本发明的色谱基线降噪装置整体结构图；

图4本发明的整体滤波效果图；

图5本发明滤波和只用SG滤波器的对比结果图；

图6本发明滤波和只用形态学滤波器和的对比结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1-6，本发明提供一种技术方案：

一种色谱基线降噪方法，包括以下步骤：

S1：信号采集：通过采集设备，采集设备为气相色谱仪，采集信号；

S2：信号放大：通过信号放大器，将步骤S1中采集的信号进行放大，并将放大后的信号通过USB传给PC端；

S3：高斯噪声去除：通过PC端，利用SG滤波器，以最小二乘滤波算法去除高斯噪声；

S4：脉冲噪声去除：利用形态学滤波器，对步骤S3去除高斯噪声的信号，以形态学滤波算法去除脉冲噪声；

S5：标定拟合：向采集设备中先通入已知浓度及成分的标准气体，采集设备为气相色谱仪，通过观察谱图曲线，对每个峰的成分及浓度进行标定，对每种物质的峰高和浓度关系进行拟合，获得拟合曲线，标定拟合中拟合三次；

S6：成分分析：对步骤S4中去除脉冲噪音的信号，结合步骤S5中的拟合曲线，进行分析对比，得出结果，成分分析采用MATLAB自带的findpeaks 寻峰方法，对色谱峰进行分析，进而获取峰高度、宽度、位置等信息。

其装置，包括气相色谱仪以及与气相色谱仪连接的放大器，放大器通过 USB数据线连接有PC端，且PC端连接有SG滤波器，SG滤波器串联有形态学滤波器。

Savitzky-Golay滤波器(简称SG滤波器)，是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变。其基本思想：基于多项式，通过移动窗口利用最小二乘法对数据进行最佳拟合，并且提供了一种利用多项式卷积的计算方法来计算滤波系数。

步骤：

1.构造多项式：

设一组信号数据为x[n]，n的取值为2m+1个连续的整值，即i＝-m,...,0,...,m。现构造一个n阶多项式(n≤2m+1)来拟合这一组数据：

2.构造点的误差平方和公式：

3.求解最优化问题

为了使得ε_N最小，可令其对各系数的偏导为0，即

化简后可得：

我们引入一个(2M+1)行(N+1)列的辅助矩阵A：

A＝{a_n,i}，其中a_n,i＝nⁱ,-M≤n≤M,0≤i≤N

再设一个辅助矩阵B，使得：

B＝A^TA

经计算可知：

定义：

则可以得到：

Ba＝A^TAa＝A^Tx

a＝(A^TA)^-1A^Tx，令H＝(A^TA)^-1A^T，则a＝Hx

4.分析确定滤波结果

定义滤波结果y[n]如下：

多项式阶数为N，对于每一个时间样本n，输入的信号数据集为x[n]，定义多项式p[n]在n＝0处的值为滤波结果y[n]。即

由此可见，只需要获得拟合多项式的常数项，等价于得到滤波结果。

经过计算可知，输入信号X已知，只需确定矩阵H的第一行行向量，即可确定a₀。

5.求解H矩阵，即滤波系数

由上述内容可知，只需求解矩阵H的第一行行向量。

令x'＝(0,0,...,1,...,0,0)^T

已知a'＝(A^TA)^-1A^Tx',

所以A^Tx'＝(1,0,...,0)^T。

那么(A^TA)^-1的第一列是(a'₀,a₁',...,a'_N)^T

又因为(A^TA)^-1是对称矩阵，所以(A^TA)^-1的第一行是(a'₀,a'₁,...,a'_N)

至此可以确定矩阵H的第一行为(p'(-M),p'(-M+1),...,p'(M-1),p'(M)),其中

6.用卷积求解滤波结果

由步骤4可知，y[n]＝a₀,并且从式11可知a₀是x[n]的线性组合，即y[n] 是x[n]的线性组合。令h[n]＝H[0,-n]作为单位脉冲，进而可以写成卷积的形式：

y[n]＝x[n]*h_N[n]

由此可以求得滤波结果。

并且，由于所选择的窗口是对称的，所以卷积计算可以简化为：

最终，SG滤波算法就变为用已知系数的窗口去扫描信号数据。

形态滤波的理论依据是数学形态学，在形状识别、边缘检测、纹理分析等方面应用广泛，其一维离散情况下的多值形态学变换，对抑制信号中的波峰(正脉冲)噪声、波谷(负脉冲)噪声，以及白噪声有很好的滤波效果。

该方法计算简单、速度快、对硬件要求较低，可广泛应用于各种实时信号***。

1、形态学基本运算定义：腐蚀、膨胀、形态开、形态闭

定义1 f(n)和g(m)分别为定义在F＝{0,1,...,N-1}和G＝{0,1,...,M-1}上的离散函数，且N＞＞M。其中f(n)为输入序列，g(m)为结构元素。

定义2 f(n)关于g(m)的腐蚀和膨胀分别为：

式中符号Θ和

分别表示腐蚀和膨胀运算。

定义3 f(n)关于g(m)的形态开和形态闭分别为：

式中符号

和·分别表示形态开和形态闭运算。

定义4 f(n)关于g(m)的形态开-闭(OC)和形态闭-开(CO)滤波器分别为：

2、开、闭运算的性质

形态开可以抑制信号中的波峰(正脉冲)噪声，具有反扩张性。而形态闭可以抑制信号中的波谷(负脉冲)噪声，且使邻近目标连接起来，具有扩张性。

3、结构元素的选取

结构元素的选取有如下两点依据：

1)结构元素横向尺寸(元素个数)要大于或等于要滤除的干扰脉冲横向宽度 (即干扰脉冲采样个数)。

根据待检测成分的半峰宽来确定结构元素的长度，可以在保留峰型的同时滤除脉冲噪声。

2)结构元纵向(幅度)的选取，主要是按要处理信号的细节形状来选取。

由于结构元素个数远小于采样分析的元素个数，可以认为结构元素的形状对波形没有影响。为简化运算，提高处理速度，可以将结构元幅值取为0。

综上，以待检测物质的半峰宽为5为例，本方法选取矩形结构元素为 {0,0,0,0,0}。注意，窗口大小一定要选为奇数，因为奇数大小的滤波器更容易确定其中心位置，而偶数大小的滤波器需要更加特定的指定一个位置为中心点。

4、滤波器的设计分析

由步骤2可知形态开可以用于过滤信号上方的峰值噪声，去除毛刺及小桥结构；而形态闭运算可用于平滑或抑制信号下方的波谷噪声，填平小沟结构。也就是说，形态开的输出信号总是位于原信号下方，而形态闭的输出信号总是位于原信号的上方。由此引发以下几个问题：

1)单独使用形态开运算或形态闭运算都会对原始信号造成更大的误差

2)在结构元素相同的情况下，交替使用形态开、形态闭运算。“先开后闭”或者“先闭后开”都会造成信号输出单向偏移的问题。

3)单纯地将形态开和形态闭运算的结果取平均，当结构元素较宽时，不能有效滤除脉冲噪声。

综上，结合形态学开闭运算的特点，并且针对上述几个问题，本方法提出交替混合形态学滤波方法(公式如下)：选择两种大小的结构元素，分别为 g1,g2,且g2>g1；同时进行a)形态开闭运算，及b)形态闭开运算，再将两个结果取平均。注意，对于a)方法，闭运算的结构元素要大于开运算的结构元素，用于滤除因开运算加大的负脉冲；对于b)方法，开运算的结构元素要大于闭运算的结构元素，用于弥补因闭运算减小的正脉冲。

具体步骤如下：

1、选取矩形结构元素{0,0,0,0,0}

2、执行形态开闭运算

a)收集数据：给定长度为g1的数组RawData[]，当输入数据不足g1个时，返回原始数据；当输入数据满足g1个时，执行步骤b)；

b)膨胀：计算数组RawData[]的最小值，存入长度为g1的数组ListOpenClose0[]，用于腐蚀计算；

c)腐蚀：当膨胀的计算结果不足g1个时，继续收集数据；当膨胀的计算结果满足g1个时，计算数组ListOpenClose0[]的最大值，存入长度为g2的数组 ListOpenClose1[]。

d)腐蚀：当腐蚀的计算结果不足g2个时，继续收集数据；当腐蚀的计算结果满足g2个时，计算数组ListOpenClose1[]的最大值，存入长度为g2的数组 ListOpenClose2[]。

e)膨胀：当腐蚀的计算结果不足g2个时，继续收集数据；当腐蚀的计算结果满足g2个时，计算数组ListOpenClose2[]的最小值，即为某时刻开闭运算的结果。

其中步骤a)，b)，c)统称开运算，步骤d)，e)统称为闭运算。

3、执行形态闭开运算，与步骤2同时进行

a)收集数据：给定长度为g1的数组RawData[]，当输入数据不足g1个时，返回原始数据；当输入数据满足g1个时，执行步骤b)；

b)腐蚀：计算数组RawData[]的最大值，存入长度为g1的数组ListOpenClose0[]，用于膨胀计算；

c)膨胀：当腐蚀的计算结果不足g1个时，继续收集数据；当腐蚀的计算结果满足g1个时，计算数组ListCloseOpen0[]的最小值，存入长度为g2 的数组ListCloseOpen1[]。

d)膨胀：当膨胀的计算结果不足g2个时，继续收集数据；当膨胀的计算结果满足g2个时，计算数组ListCloseOpen1[]的最小值，存入长度为g2 的数组ListCloseOpen2[]。

e)腐蚀：当膨胀的计算结果不足g2个时，继续收集数据；当膨胀的计算结果满足g2个时，计算数组ListCloseOpen2[]的最大值，即为某时刻闭开运算的结果。

其中步骤a)，b)，c)统称闭运算，步骤d)，e)统称为开运算。

4、将步骤2和3的结果取平均作为某时刻的滤波结果

5、重复执行步骤2,3,4

流程图参见附图2。

这里主要四个过程，数据采集、基线降噪、标定拟合以及成分分析。

1.数据采集：

即采集待检测物质的色谱数据信号。待检测物质通过氢火焰离子检测器，形成离子流被收集、输出，经阻抗转化，放大器(放大107～1010倍)便获得可测量的电信号，本装置获取该电信号，即完成数据采集。本装置可以显示谱图曲线。

2.基线降噪(SG\形态学滤波器)

首先确定滤波窗口大小：

SG滤波器的窗口宽度：19；

形态学滤波器的两个结构元素：g1＝{0 0 0 0 0}，g2＝{0 0 0 0 0 0 0}。

为了达到实时滤波的效果，本发明中SG滤波器和形态学滤波器两种方法依次且同时进行，每次只处理一个数据，两种窗口每次都只移动一格，只是窗口宽度不同。

第一步，采集数据流中的数据到一个宽度为19的数据容器，

若容器中的数据个数小于19，则返回原始值；

若容器中的数据个数等于19，则开始计算：利用最小二乘法进行多项式拟合，确定多项式系数，即为滤波窗口系数，将系数与窗口中的数据对应做卷积运算，得到滤波结果，将滤波结果赋值给窗口中数据的中心点。将经过滤波的数据传递给下一个数据容器用于形态学的计算。

接着移动窗口，也就是删除第一个数据，增加最后一个(新的)数据，重复滤波运算。

第二步，用于形态学计算的容器中存有经过SG滤波器滤波的数据，由于 SG滤波算法一直在进行，所以这也是以数据流的形式呈现。

同样的，若容器中的数据个数小于5，则直接返回容器中的数据；

若容器中的数据个数等于5，则依次进行形态学开闭运算，得到计算结果。

然后移动滤波窗口。

形态学方法分两路同时进行，分别是“先开后闭”和“先闭后开”过程，再将两个结果取平均，作为最终的滤波结果。

经过上述过程，获得一个数据的滤波结果。随着窗口的移动，实时获取相应的滤波结果。

本发明中的方法均是下一级收集上一级传来的数据，如果窗口中数据未填满，就返回原始值，如果填满就进行相应运算处理。移动窗口就是删除第一个数据，增加最后一个数据。而窗口扫描就是用窗口系数乘以窗口中对应的数据再相加，结果作为中心数据的滤波结果。

整体滤波效果见附图4。

SG滤波器和本方法的滤波效果对比见附图5。

形态学滤波器和本方法的滤波效果对比见附图6。

3.标定拟合

先通入已知浓度及成分的标准气体，通过观察谱图曲线，对每个峰的成分及浓度进行标定，标定两到三次，对每种物质的峰高和浓度关系进行拟合，获得拟合曲线Y＝a*X+b，其中Y是浓度，X是峰高。之后将峰高代入曲线方程，即可求得所测实时浓度。

两次实验，分别测了两种物质：

总烃：标定浓度17.3，标定高度27.171，

标定浓度17.6，标定高度27.642

标定关系曲线y＝0+0.6367*x

甲烷：标定浓度12.5，标定高度6.633，

标定浓度12.9，标定高度6.845

标定关系曲线y＝0+1.8844*x(注：小数位数有取舍)

4.成分分析

经过滤波之后，得到较为准确的真实信号，可以进行成分分析。所谓的成分分析就是根据谱图的出峰情况，分析待检测物质中的成分及浓度。

本模块使用MATLAB自带的findpeaks寻峰方法对色谱峰进行分析，进而获取峰高度、宽度、位置等信息。

本次实验分析了两种物质，分别是甲烷和总烃。测量的气体是实验室内的空气成分。测量过程在白天进行，时间为10小时，整体环境因素基本保持稳定。整个实验每1小时查看一次处理结果。下表为以总烃为例展示实验分析结果：

一般的，对于相同浓度的样本气体测量时连续进样保留时间的变化要满足小于±1％的要求，峰高、峰面积的变化要满足小于±3％的要求。由上表可知，测定10次数据后测得的保留时间、峰高和峰面积的相对标准差分别只有 0.73％，0.91％，和1.26％，所以以保留时间为标准的定性重复性分析满足±1％的要求，同时，以峰高及峰面积为标准的定量重复性分析满足小于±3％的要求。

注：相对标准偏差(RSD)的计算公式如下：

其中S为标准偏差(也可以表示为SD)，

为相应的平均值。

本发明利用SG滤波器串联有形态学滤波器，综合两者的优点，摒除两者的缺点，先用SG最小二乘滤波算法去除高斯噪声，再用形态学滤波算法去除脉冲噪声。这样做的好处是，先用SG最小二乘滤波器去除大部分噪声，免除形态学滤波器对噪声的放大作用；再用形态学滤波器，可以滤除SG滤波器残留的异常点干扰问题。两个步骤顺序不可颠倒，否则会放大两种算法的缺点。最后，本发明中的两种滤波算法都是实时算法，可以实时对采集的信号进行滤波处理，延迟在0.5s以内。

另外，本发明具有计算简单，延迟小，可以实时滤波等优点，且在保持峰高和峰型不变的情况下可以有效去除脉冲噪声和高斯噪声；并且使用本方法可以提高整个装置的检测精度，即色谱峰的识别准确度，可以检测浓度在最低检测浓度可达0.005ppm级别的物质。本发明在保持峰高和峰型不变的情况下去除脉冲噪声和高斯噪声，且计算简单，延迟小，可以实时滤波；并且提高了整个装置的检测精度，即色谱峰的识别准确度，可以检测浓度在0.01 级别的物质。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (5)

1.一种色谱基线降噪方法，其特征在于：包括以下步骤：

S1：信号采集：通过采集设备，采集信号；

S2：信号放大：通过信号放大器，将步骤S1中采集的信号进行放大，并将放大后的信号通过USB传给PC端；

S3：高斯噪声去除：通过PC端，利用SG滤波器，以最小二乘滤波算法去除高斯噪声；

S4：脉冲噪声去除：利用形态学滤波器，对步骤S3去除高斯噪声的信号，以形态学滤波算法去除脉冲噪声；

S5：标定拟合：向采集设备中先通入已知浓度及成分的标准气体，通过观察谱图曲线，对每个峰的成分及浓度进行标定，对每种物质的峰高和浓度关系进行拟合，获得拟合曲线；

S6：成分分析：对步骤S4中去除脉冲噪音的信号，结合步骤S5中的拟合曲线，进行分析对比，得出结果。

2.根据权利要求1所述的一种色谱基线降噪方法，其特征在于：所述采集设备为气相色谱仪。

3.根据权利要求2所述的一种色谱基线降噪方法，其特征在于：所述标定拟合中拟合三次。

4.根据权利要求3所述的一种色谱基线降噪方法，其特征在于：所述成分分析采用MATLAB自带的findpeaks寻峰方法，对色谱峰进行分析，进而获取峰高度、宽度、位置等信息。

5.一种用于权利要求4所述的色谱基线降噪方法的装置，其特征在于：包括气相色谱仪以及与气相色谱仪连接的放大器，所述放大器通过USB数据线连接有PC端，且PC端连接有SG滤波器，所述SG滤波器串联有形态学滤波器。

Images