基于数字滤波与自相关的涡街流量计抗强振动方法改进

徐科军， 王 沁， 方 敏， 邵春莉

（1.合肥工业大学电气与自动化工程学院，安徽合肥 230009；2.工业自动化安徽省工程技术研究中心，安徽合肥 230009）

基于数字滤波与自相关的涡街流量计抗强振动方法改进

徐科军1，2， 王 沁1， 方 敏1， 邵春莉1

（1.合肥工业大学电气与自动化工程学院，安徽合肥 230009；2.工业自动化安徽省工程技术研究中心，安徽合肥 230009）

涡街流量计容易受到管道振动的干扰，对测量精度造成较大影响，特别是当振动噪声的功率大于涡街信号的功率时，采用通常的数字信号处理方法也无法得到正确的涡街频率。为此，提出将数字滤波和自相关相结合的算法，运用于涡街流量传感器的信号处理，从含有强振动干扰的信号中识别出流量信息。针对该方法在实际实现中遇到的问题，对已有系统软件进行改进和完善。基于MSP430F5418单片机，使用汇编语言编制实数FFT算法实现了自相关函数的运算。

计量学；涡街流量计；数字滤波；自相关；强振动干扰；FFT；单片机

1 引 言

涡街流量计由于无机械可动部件，安装方便，量程比较宽，并且可以测量气体和液体介质等优点，被广泛应用于石油、化工等领域。但是，在实际应用中，由于涡街流量计普遍采用应力式传感方式，极易受到管道振动的干扰，无法保证测量精度，有时还出现错误的结果［1］。采用数字信号处理方法可以从含有噪声的信号中提取流量信息［2～6］。但是，当振动噪声的功率大于涡街信号的功率时，普通的数字信号处理方法也无法得到正确的流量信息。为了解决这个问题，文献［7］提出了基于数字滤波和自相关算法相结合的抗强振动干扰处理方法，并用超低功耗单片机实现，取得较好的实验结果。但是，要在实际应用中取得好的应用效果，还需要进一步改进和完善。为此，本文基于MSP430F5418单片机，编制汇编语言的实数FFT算法完成自相关函数的运算，提高计算速度；改进自相关衰减系数的计算方法，保证测量精度。

2 抗强振动干扰方法

2.1 基本原理

涡街流量信号在理想情况下为规则稳定的正弦信号，涡街体积流量与该正弦信号频率成正比。但是，在实际工作环境下，由于水泵或其它因素的影响均会引起管道内流场不稳定，因此实际的涡街流量信号并不是完全稳定的正弦波，其频率会因为干扰的影响存在晃动。而强振动信号一般由机泵或阀门开关等现场设备的机械振动引起，因此，其频率为现场设备的振动频率或其倍频，一般较固定。文献［8］通过对实验现场采集的大量数据进行分析和建模，以功率谱计算为基础，定性地分析了涡街信号和振动信号的带宽特征，得出涡街信号为宽带信号，而振动信号为窄带信号的结论。文献［7］分别对强振动信号和涡街信号进行自相关计算，得到不同的衰减特性。自相关计算可以反映数组序列不同时刻数据的相关性，因为强振动噪声信号为窄带信号，其频率值波动范围较小，而涡街信号为宽带信号，频率波动较大，因此强振动噪声的相关性要优于涡街信号，这一特征在自相关函数上表现为，涡街信号的自相关函数随着时间的延时，其幅值有较大程度的衰减，见图1；而强振动噪声则衰减较小，见图2。

图1 涡街信号自相关衰减函数图

图2 强振动信号自相关衰减函数图

为此，文献［7］提出了基于数字滤波和自相关相结合的涡街流量传感器抗强干扰信号处理方法，具体处理过程为：首先对采集到的时域信号进行频谱分析，按照一定的信噪比及幅值门限确定频率个数，可能为0，可能为1，也可能大于1。当频率个数为0时，系统输出0表示没有流量。当频率个数为1时，会有3种情况：（1）只有涡街流量；（2）只有噪声；（3）涡街流量信号与噪声信号同频。可以求原信号延迟一段时间后的自相关函数在半个周期内的幅值与在时间为0时的自相关函数的值的比值，称该比值为衰减系数，当衰减系数小于设定的阈值时，即只有涡街流量信号或涡街流量信号与噪声信号同频，输出涡街流量频率；当衰减系数大于等于阈值时，即只有噪声，输出0表示没有流量。当频率的个数大于1时，要先采用数字带通滤波器对信号进行滤波，将滤波以后得到的相应频率的信号进行自相关计算，然后计算0.5 s附近的正弦信号波峰值与自相关函数的第一点的比值，然后判断不同频率点信号所得到的衰减系数的大小，其中衰减系数最小的频率点对应的信号即为涡街信号。

2.2 几点局限

该方法根据振动信号和涡街信号自相关函数的衰减程度的不同，作为去除噪声信号的依据，在一定程度上解决了抗强振动干扰的问题。但是，由于为了满足低功耗、两线制过程控制仪表的要求，无法选用DSP去实现算法，只能选用低功耗单片机，而单片机内存单元和计算速度都很有限，所以，该方法在实时实现过程中还存在以下3方面的问题：

（1）该算法中采用了数字带通滤波器，对原始的涡街流量传感器输出信号进行滤波，得到相应频率点的信号。在Matlab仿真中，设定的带通滤波器带宽较窄，仅为［－5 Hz，5 Hz］。由于单片机的内存单元有限，只能进行2048点的滤波，而数字带通滤波器的带宽较窄，因此会有较长时间的收敛过程，如图3和图4所示。

图3 带通滤波后振动信号

当涡街信号频率和强振动信号频率值很靠近时，经过带宽很窄的数字带通滤波器滤波后，信号会出现较大幅度的波动；并且经过滤波后的振动信号以及涡街信号在前300个点均为收敛过程，同时，整个范围内幅值也存在失真的情况，这对下一步的自相关计算造成不良影响。

图4 带通滤波后涡街信号

（2）自相关函数可以描述随机过程中任意2个不同时刻之间相关性，反映函数在2个不同时刻的状态之间的统计关联程度。自相关估计一般需要点数较多，在Matlab中使用了16384个点来进行仿真，但是，在系统实时实现过程中单片机内存的限制，只能进行2048个点的自相关计算，点数减少，从而降低自相关计算的准确性。此外，用来计算自相关的序列是经过数字滤波得到的，而当数字滤波的结果不能真实反映相应频率点信号时，也会对自相关函数造成影响。图5和图6分别为振动信号和涡街信号的2048个点的自相关函数图。由图5可以看出，振动信号的自相关函数也存在一定的衰减，而由6可以看出，涡街信号的自相关函数的衰减曲线存在波动。

图5 振动信号自相关函数

（3）该方法采用无偏估计进行自相关计算：

图6 涡街信号自相关函数

式中：R为数据个数；r为延迟时间；x（n）为滤波后的原始信号；n＝0，1，2，…，R－r－1。

采用式（1）的方法乘法次数多，运算量较大。因此，为了节省计算时间，提高系统实时性，原来的系统中根据振动自相关函数和涡街信号自相关函数的特点进行了计算上的简化，因为在仿真中振动信号和涡街信号自相关函数的衰减规律都是单向衰减，并且振动信号自相关函数衰减较慢，因此我们一般选择0.5 s延时附近的峰值点与自相关函数的第一个峰值点进行比较，从而计算出振动信号和涡街信号的衰减系数，然后进行比较。即按照式（1）只计算了1个点的自相关系数。但是，经过分析，由图5和图6所示振动信号或者涡街信号的自相关函数并不完全是单向衰减的，其原因是滤波后信号的失真以及自相关点数不够多，因此原来软件中只计算0.5 s延时处的一点的自相关系数，会造成结果的误判。

3 方法的改进

针对上述局限性，对基于单传感器的抗强干扰系统软件进行了相应的改进。由于MSP430F5418单片机的内存为16 kB，为了能够在单片机上实现抗强干扰算法，自相关计算点数的上限为2048个点。因此不再仅仅计算0.5 s延时处的自相关系数，而是将2048个点的自相关系数都计算出来，然后在找出第一个最大衰减的点数，根据这一点与第1个峰值点的衰减比值来消除噪声。

3.1 采用实数FFT实现自相关计算

基于MSP430单片机用实数FFT方法实现自相关计算；自相关函数的无偏估计公式如式（1）所示，如果计算过程中R和r均较大，则计算过程非常繁琐，占用时间较长，因此很难在单片机中进行应用。自相关算法还有另一种估计方法即有偏估计方法：

式中~Rxx的点数为2R－1，对于规定的R值，当r远小于R时，~Rxx估计值才越接近真值。对~Rxx求傅里叶变换，则有

式中ω为圆周频率。

又由于2个R点的序列进行线性卷积，将得到一个长为2R－1点的序列，如果用DFT来实现线性卷积，则需要将序列补零延长至2R－1点，从而得到序列x2R（n）。由于x2R（n）的后R－1点均为0，因此，式（3）可以进行如下变换：

设m＝r＋n，且0≤m≤2R－1，则式（4）右边可以写成如下所示：

式中X2R（ejω）为序列x2R（n）的傅里叶变换。

由式（3）和式（5）可知：

式（6）说明，序列的有偏自相关估计和该序列的功率谱是一对傅里叶变换。而傅里叶变换可以通过FFT实现，因此，在单片机中通过FFT计算实现2048个点的自相关计算。

在以基于MSP430F5418单片机为核心的硬件系统上完成软件研制，因此算法的实时实现需要合理分配单片机资源。单片机内存为16kB，采样数组为2400个点，数据类型为整型，占用2 B，保存数组需要约4kB的内存，因此自相关计算最多只能选择2048个点。为了能通过FFT实现，首先需要将数据补零至4096个点，数据类型为整型，需要占用内存空间8 kB。

实现过程中利用实数FFT来进行自相关计算，而实数FFT是将序列按照奇偶分别存放于奇数序列和偶数序列。设存放原来序列的数组为x（n），长度为2N，奇数项h（r）＝x（2r＋1），其中r＝0，1，…，N－1，奇数项作为FFT计算的虚部；偶数项g（r）＝x（2r），其中r＝0，1，…，N，偶数项作为FFT计算的实部。具体实现过程如下：

（1）将经过数字滤波后的2048个点存放于x（n）中，该序列长度为2048点，然后通过补零的方法将序列x（n）补到4096个点，为了节省单片机内存，补出的0不是直接放于x（n）中，而是按照奇偶项直接分配到h（r）＝x（2r＋1）和g（r）＝x（2r）中，h（r）为2048点，后1024个点为0，同理g（r）也为2048个点，后1024个点也为0。

（2）h（r）和g（r）序列就绪后，进行实数FFT计算，由于h（r）和g（r）的总点数为4096个点，因此完成的实际是4096个点的FFT计算，设得到的傅里叶变换序列为X（k），因为X（k）是关于中心共轭对称的，因此只保存前2048个点的实部和虚部，后半部分根据其共轭对称的关系求出，其中序列h（r）和g（r）中分别存放着前2048个点的虚部和实部。同时求出X（k）的共轭序列X*（k）。

（3）根据第（2）歩计算X（k）和X*（k）的乘积，他们为共轭序列，因此乘积为实数，又因为X（k）后2048个点与前2048个点为对称共轭的关系，因此在这个过程中可以计算出4096个点的乘积值，按顺序存放于h（r）和g（r）序列中。

（4）对式（3）中求得的乘积进行IFFT计算，其计算公式为：

式中N＝2 048为自相关实际点数。其中利用FFT的子程序完成IFFT的计算。

3.2 优化自相关衰减系数计算方法

为了满足涡街流量计的低功耗要求，选择MSP430单片机作为处理器。但是，目前430系列单片机最大内存仅为18 kB，仍然不能实现多点数的数字滤波器和多点数的自相关计算，因此根据图（5）和图（6）的特性，在自相关衰减计算中做了相应的调整，以保证算法的有效性。虽然少点数的自相关特性有些变化，但是涡街信号的自相关函数仍然比振动信号衰减得更明显，由于滤波造成的失真，使得涡街信号自相关函数会先出现快速衰减的趋势，然后又再上升；而振动信号的自相关函数在越接近N的位置也会衰减比较快。因此不再选择0.5 s的位置来比较二者的衰减系数。对于涡街信号，通过逐点寻找，找出第一个衰减最大的正弦峰值点；对于振动信号仍采用此方法，如果自相关函数的衰减趋势为单向的，则仍选择1024个点附近的正弦峰值点。

4 实验结果

分别在重庆川仪自动化股份有限公司和合肥工业大学实验室的振动平台上做了抗强振动干扰的实验测试，验证在存在强振动干扰的情况下改进方法和系统的有效性。通过鼓风机吹风模拟气体流量，由振动实验平台模拟强振动干扰，并施加于涡街流量计上。分别选择膜盒式和悬臂梁式涡街传感器进行实验。

在重庆川仪自动化股份有限公司，针对DN50口径膜盒式涡街传感器进行抗强振动的实验。为了验证抗强干扰算法的有效性，在实验过程中使振动信号的功率大于涡街信号的功率。以振动信号频率为110 Hz，涡街信号为100 Hz时的功率谱图为例，从图7中可以看出振动信号的功率比涡街信号的大。该图是用数字存储示波器采集涡街流量传感器输出信号，用Matlab分析的结果。在实验过程中，选择了3个不同的流量点，分别大约为101 Hz，130 Hz以及148 Hz，在这3个流量点下，通过改变振动台的振动频率以及振动加速度来进行考核。当振动频率较低时，其振感很强，仪表液晶在振动的作用下其读数有些模糊。用基于MSP430F5418单片机的抗强振动型涡街流量计实时测量和处理，实验结果表明，改变振动频率和加速度，仪表均能显示正确的涡街频率，如表1所示。

图7 膜盒式涡街传感器输出信号功率谱

在合肥工业大学实验室搭建的振动实验平台上，采用DN50口径悬臂梁式涡街传感器进行抗强振动实验。实验过程中，当振动实验台加速度为0.8g及以上时，其振动强度能够将放于振动台表面的小工件慢慢振落下去。同样地，选择3个流量点分别大约为101 Hz、88 Hz和121 Hz，图8所示为用数字存储示波器采集涡街流量传感器输出信号，用Matlab分析的结果。用基于MSP430F5418单片机的抗强振动型涡街流量计实时测量和处理，实验结果如表2所示，可见，仪表在不同频率振动的影响下，均能正确显示涡街信号频率。

表1 DN50膜盒式涡街传感器抗强振动实验结果

图8 悬臂梁式涡街传感器输出信号功率谱

表2 DN50悬臂梁式涡街传感器抗强振动实验结果

5 结束语

根据现场实验，对数字滤波和自相关算法在单片机实现过程中存在的局限性进行分析，通过示波器采样数据在Matlab中的分析结果与单片机计算过程中得到的结果进行比较，分析得出：由于单片机内存较小，只能做2048点的滤波计算，因此滤波结果存在失真；单片机资源有限，选择了2048点自相关计算，无法达到Matlab的仿真效果。以现有单片机资源为考虑因素，进行了2方面的改进：一方面利用实数FFT实现自相关计算，将2048点的自相关计算快速计算出，节省计算时间。另一方面改进了自相关衰减系数计算方法，从而保证算法计算的准确性。

致谢：感谢重庆川仪自动化股份有限公司提供实验条件和涡街流量计一次仪表。

［1］ Miau J J，Hu C C，Chou J H.Response of a vortex flowmeter to impulsive vibrations［J］.FlowMeasurementandInstrumentation，2000，11：41－49.

［2］ Schlatter G L，Barrett W D，Waers J F，etal．Signal processing method and apparatus for flowmeters：WO，1990004230［P］.1990－04－19.

［3］ Clark D W.Designing phase-locked loops for instrumentation applications［J］.Measurement，2002，32：205－227.

［4］ Clark D W，Ghaoud T.A dual phase locked loop for vortex flow metering［J］.FlowMeasurementand Instrumentation，2003，14：1－11.

［5］ 曾宪俊，徐科军，周杨，等.低功耗、两线制涡街流量计数字信号处理系统［J］.仪表技术与传感器，2009，（1）：76－79.

［6］ 徐科军，刘三山，刘家祥，等.基于数字频谱与模拟带通滤波组的两线制涡街流量计［J］.计量学报，2011，32（1）：25－30.

［7］ Xu K J，Luo Q L，Wang G.Frequency-feature based anti-strong disturbance signal processing method and system for vortex flowmeterwith single sensor［J］.Review ofScientificInstuments，2010，7：075104－7.

［8］ 徐科军，罗清林，朱永强，等.涡街流量传感器信号的幅值和频率特征［J］.计量学报，2010，31（6）：504－508.

Im provements of Anti-strong-vibration Method Combing Digital Filter w ith Auto-correlation for Vortex Flowmeter

XU Ke-jun1，2， WANG Qin1， FANG Min1SHAO Chun-li1
（1.School of Electrical and Automation Engineering，Hefei University of Technology，Hefei，Anhui230009，China；
2.Engineering Technology Research Center of Industrial Automation，Anhui Province，Hefei，Anhui230009，China）

Vortex flowmeter is susceptible to the interference of pipe vibration，which affects on the measurement accuracy greatly.Adopting normal digital signal processingmethods cannotobtain the accurate vortex frequencieswhen the power of vibration noise is larger than that of the vortex signal.Therefore a digital signal processing technique combining digital filter with autocorrelation was proposed to apply the signal processing of vortex flow sensor for extracting the flow rate information from the signal containing strong vibration interference.Aiming at the implementing problems of thismethod in practice application，the existing system software is improved.The autocorrelation function is calculated by using real number FFT algorithm with assemble language based on MSP430F5418 MCU in order to speed up computation.

Metrology；Vortex flowmeter；Digital filter；Autocorrelation；Strong vibration interference；FFT；MCU

TB937

A

1000-1158（2014）05-0449-06

10．3969／j．issn．1000-1158．2014．05．09

2012-08-10；

2013-01-23

国家自然科学基金（61074164）

徐科军（1956－），男，江苏无锡人，合肥工业大学教授、博士生导师，主要从事传感器技术、自动化仪表和数字信号处理研究。dsplab＠hfut.edu.cn