CN106768109B - 基于回波上升段峰值拟合和基于回波能量点定位的气体超声波流量计信号处理方法 - Google Patents

Abstract

基于回波上升段峰值拟合的气体超声波流量计信号处理方法，首先根据拟合峰值范围，在回波信号的上升阶段，选取若干个峰值点；对各峰值点进行最小二乘拟合，得到特征直线；将该直线与横轴采样点的交点作为特征点，确定超声波信号的顺流、逆流传播时间，并计算气体流量。基于回波能量点定位的气体超声波流量计信号处理方法，首先以定位能量点的能量值作为基准，找到四个最邻近的峰值能量点进行直线拟合；再找出拟合直线上该能量值对应的能量点作为特征点，确定超声波信号的顺流、逆流传播时间，并计算气体流量。

Description

基于回波上升段峰值拟合和基于回波能量点定位的气体超声 波流量计信号处理方法

技术领域

本发明涉及流量测量领域，为气体超声波流量计信号处理方法，特别是一种基于回波上升峰值拟合和一种基于回波能量点定位的气体超声波流量计信号处理方法。采用本发明提出的两种气体超声波流量计信号处理方法处理超声回波信号，实现气体流量的测量。

背景技术

气体超声波流量计由于其精度高、无压损、量程比大等诸多优点，被广泛应用于大口径管道的气体流量测量。气体超声波流量计的测量方法(原理)主要包括传播时间差法、传播频率差法和多普勒法等。其中，传播时间差法是气体超声波流量计采用的最为广泛和有效的方法。传播时间差式气体超声流量计测量流量时，需要采用激励信号驱动发射换能器产生超声波信号，然后再由接收换能器将超声波信号转换为回波信号；通过测量激励信号与回波信号之间的时间差来计算出相应的顺流、逆流传播时间，进而计算气体流量。激励信号起始时刻是可以精确控制的，而回波信号到达时刻则要根据回波信号中的某一稳定的特征点来确定的。由于超声波信号在气体中传播时，其能量会衰减，并且随着气体流速的增加，超声波信号的传播路径会产生偏移，这种偏移将进一步导致能量的减弱，使得回波信号幅值变小，信噪比降低，容易受噪声等干扰。而直接通过最小峰值点、最大峰值点等简单方法确定的特征点，都容易产生波动，影响流量测量的精度。

国内外学者对气体超声流量计信号处理方法进行了研究，主要有以下几种方法。

(1)基于互相关和过零检测相结合的方法

德国西门子股份公司采用一种过零检测与互相关相结合的方法计算超声波顺流和逆流的绝对传播时间(Arthur Freund,Nils Kroemer.Method for measuring the timeof flight of electric,electromagnetic or acoustic signals,EP PatentNO.0797105A2,Mar.17,1997)，首先通过某种方法获得一段不含噪声的回波信号，假设其起始时刻为t₀，同时在该回波信号中选定某一个过零点，该过零点距离起始时刻t₀的时间间隔记为t_N。然后，将该段回波信号与换能器实际采集的回波信号做互相关运算。如果互相关运算结果的最大值表征的时间间隔为t_d，同时，假设实际回波信号的起始时刻相对于激励信号发射时刻的时间间隔为t₁，那么回波信号的传播时间为：

t＝t₁+t_d+t_N-t_korr (1)

式中，t_korr为修正时间，它包括换能器转换延时和电路传输延迟。这种方法能够可靠地从混杂有噪声信号的接收波形中检测出回波信号，并且由于过零点附近具有较高的幅值分辨率，所以，适合于嘈杂的工业现场。但是，该专利中没有披露如何获取不含噪声的回波信号；同时，还存在互相关运算量大等问题。

(2)基于互相关及其衍生的方法

由于在超声波流量计工作工程中，顺流和逆流时产生的回波信号之间具有相关性，而噪声与信号之间不具有相关性，因此，可以利用互相关及其衍生方法计算超声波的传播时间差。传统的互相关计算如式(2)所示，其中，假设顺流时的回波信号为x(t)，逆流时的回波为y(t)，那么，R_xy(τ)取得最大值时所对应的时间位移τ，即为顺流传播时间与逆流传播时间之差。传统的互相关方法运算量大，难以利用单片机的有限资源实时实现。为了减少互相关的运算量，日本东京计装(Keiso)株式会社采用了一种互相关衍生算法(Tokio Sugl,Tadao Sasaki.Ultrasonic flow meter,US patent NO.007299150B1,Nov.20,2007)，如式(3)所示，其中，m为移位步数，N为采样点数，那么，S_xy(m)取得最大值时所对应的移位步数m表征了顺流和逆流传播时间之差。通过对比公式(2)和(3)，不难发现，这种互相关衍生算法利用加法取模运算代替了原互相关算法中的乘法运算，减少了运算量。可是，同时也带了测量的误差，尤其当噪声干扰较大的时候，极有可能出现计算错误。

(3)基于可变阈值的过零检测方法

为了找到稳定的特征点以计算超声波传播时间，汪伟等提出一种基于可变阈值过零检测气体超声波流量信号处理方法(汪伟,徐科军,方敏,等.一种气体超声波流量计信号处理方法研究[J].电子测量与仪器学报,2015,29(9):1365-1373；徐科军，方敏，汪伟，朱文姣，沈子文。一种基于FPGA和DSP的气体超声流量计，申请中国发明专利，公开号CN104697593A，申请日2015.06.10.)，根据信号幅值的大小设定一个可变的阈值。这个阈值根据最大峰值实时调整，与最大峰值成一定比例。通过这个可变的阈值可以准确地找到回波信号的某个特征波，然后，进行过零检测，将过零点作为特征点，从而求得超声波信号顺流、逆流的传播时间，并计算出准确的气体流量。虽然该方法取得了较好的实验结果，但是，该方法未研究与之对应的效果更佳的激励信号，并且，该方法需要在测量流量之前进行阈值参数的选取，操作相对复杂。

(4)基于能量突变的方法

美国Daniel公司采用检测能量突变的方法查找特征点的方法(William Freund,Winsor Letton,James Mc-Clellan,Baocang Jia,Anni Wey,Wen Chang.Method andapparatus for measuring the time of flight of a signal,US patent NO.5983730,Nov.16,1999)。由于回波信号的能量经历了从弱到强，再从强到弱的过程，因此回波信号的能量变化率将先增大再减少，依据能量变化率的临界变化点，即可确定超声波传播时间。具体说来，此方法首先求取回波信号各点的幅值平方，幅值平方的大小表征了信号能量的大小。然后，利用滑动平均的方法，求取每一点的平均能量，再绘制出前、后两点平均能量之比的变化曲线，即能量变化率曲线。最后，通过能量变化率的临界点，确定超声波传播时间。但是，该专利中，没有披露实施过程中的关键技术，如滑动窗口的时间间隔，以及求取能量变化率时前、后两点之间的时间间隔等。

(5)基于能量变化率的方法

从能量变化的角度出发，沈子文等也提出了一种基于能量变化率的信号处理方法(沈子文,徐科军,方敏,等.基于能量变化率的气体超声波流量计信号处理方法[J].仪器仪表学报,2015,36(9):2138-2144.)，首先对采样的回波信号进行滤波，得出其各峰值点幅值，并对其求平方，作为该峰值点的能量，通过前一个峰值点能量与后面若干个峰值点能量做比较，得到每个峰值点的能量变化率，并对其进行包络线拟合，得到不同流量下的回波信号能量变化曲线，根据能量变化曲线设定能够区分特征波的阈值。实际测量时，实时求取回波信号的能量变化曲线，根据预先设定的阈值判别能量变化曲线对应的特征点。由于该方法需要在测量流量之前进行阈值参数选取，操作相对复杂。

发明内容

针对从气体超声波流量计回波信号中难以确定到达时刻，而仪表又要求具有较高实时性的问题，本发明提出了一种基于回波上升段峰值拟合和一种基于回波能量点定位的气体超声波流量计信号处理方法，并在合肥工业大学研制的双声道气体超声波流量计硬件平台上(徐科军，方敏，汪伟，朱文姣，沈子文。一种基于FPGA(现场可编程门阵列)和DSP(数字信号处理器)的气体超声流量计，申请中国发明专利，公开号CN104697593A，申请日2015.06.10)，实现基于回波上升段峰值拟合和基于回波能量点定位的气体超声波流量计信号处理系统。

基于回波上升段峰值拟合的气体超声波流量计信号处理方法的计算步骤为：首先根据拟合峰值范围，在回波信号的上升阶段，选取若干个峰值点；对各峰值点进行最小二乘拟合，得到特征直线；将该直线与横轴采样点(对应于时间)的交点作为特征点，确定超声波信号的顺流、逆流传播时间，并计算气体流量。

基于回波能量点定位的气体超声波流量计信号处理方法的计算步骤为：首先以定位能量点的能量值作为基准，找到四个最邻近的峰值能量点进行直线拟合；再找出拟合直线上该能量值对应的能量点作为特征点，确定超声波信号的顺流、逆流传播时间，并计算气体流量。

附图说明

图1是气体超声波流量计的硬件框图。

图2是DSP的软件功能框图。

图3是气体超声流量计系统主控芯片—FPGA和DSP的工作流程图。

图4是DSP的定时器中断服务程序。

图5是不同流量下的回波信号及其包络线。

图6是上升阶段包络线及其最大斜率切线。

图7是基于回波上升段峰值拟合的特征点确定方法的示意图。

图8是有流量情况下激励信号的起始时刻、回波信号起始时刻和特征点对应时刻之间的关系。

图9是零流量下激励信号的起始时刻、回波信号起始时刻和特征点对应时刻之间的关系。

图10是不同流量下的回波能量及其包络线。

图11是不同回波能量包络线的比较。

图12是回波信号上升阶段能量包络线及能量包络线变化率曲线。

图13是不同流量下回波能量峰值点比较。

图14是基于回波能量点定位的特征点确定方法的示意图。

具体实施方式

参见图1，本发明的硬件系统由换能器单元、发射/接收信号通道切换电路、激励信号产生与放大电路、回波信号调理与采样电路、DSP最小系统、FPGA最小系统和电源管理模块组成。

所述换能器单元由四个换能器组成，分别固定安装在输送气体流量的管道上；四个换能器分别为换能器1、换能器2、换能器3和换能器4，组建成为直射式双声道结构；每个换能器既作为发射换能器，又作为接收换能器。

所述DSP最小系统由DSP芯片、时钟模块、外部看门狗模块、FRAM(铁电存储器)模块、串口通讯模块和按键/液晶模块组成，其中DSP芯片的型号为TMS320F28335。

所述FPGA最小系统由FPGA芯片、时钟模块、复位电路模块和FLASH(闪存)模块组成，其中FPGA芯片的型号为EP2C8Q208C8N。其中，FPGA内部的模块包括RAM_2PORT数据存储模块、ROM数输出模块、DAC驱动控制模块、时钟分频模块和延时控制模块。

所述电源管理模块提供系统所需的模拟电源(+12V、-12V、+5V、-5V)，以及数字电源(+3.3V、+1.9V、+1.2V)。

图2所示为DSP软件功能框图，由主监控程序和各个子程序模块组成；其中，各个子程序模块包括初始化模块、发射/接收信号通道切换模块、与FPGA数据传输模块、中断模块、计算模块、FRAM读写模块、串口通讯模块、脉冲输出模块、按键输入模块以及液晶显示模块。主监控程序是整个系统的总调度程序，通过调动各个子程序模块实现系统的各项功能。

图3所示气体超声流量计系统主控芯片—FPGA和DSP的工作流程图。气体超声流量计系统上电后，FPGA进行初始化，建立相应的电路结构，并等待DSP发送激励起始信号；同时，DSP对GPIO口进行分配，完成定时器中断初始化、液晶显示初始化、从铁电存储器模块内读取累积计流量、初始化流量计各项参数等各项初始化任务。然后，DSP向FPGA发送激励起始信号；接着，DSP切换相应的激励、回波信号选通声道，并检测标志位是否为0。若为0，则直接进入等待转存起始信号状态；若为1，则更新液晶显示，完成上位机通讯，并将标志位置为0。而FPGA收到激励起始信号后，向延时控制模块发送延时起始信号，开始延时计数，并通过DAC驱动控制模块控制FPGA将内部ROM(只读存储器)数据输出模块中预存的激励信号波形数据输出至图1中所示的激励信号产生与放大电路。由该电路中的DAC(数模转换器)电路将数字信号转化为模拟信号，再经功率放大电路、电压放大电路后，通过发射/接收信号通道切换电路加载到选通的发射换能器上。发射换能器产生的超声波信号，在气体介质中传播一定时间后，由接收换能器接收，并转换为回波信号。回波信号传输至图1中所示的回波信号调理与采样电路，经放大、滤波后，再由FPGA内部的ADC采样控制模块，控制ADC以5MHz的采样频率对回波信号进行采样。采样的起始时刻由FPGA内部相应的延时控制模块决定，采样结果存放在FPGA内部的RAM_2PORT数据存储模块中。采样完成后，FPGA向DSP发送转存起始信号。DSP收到该信号后，开始转存FPGA中存储的采样数据，并采用本发明提出的基于回波上升段峰值拟合和基于回波能量点定位的信号处理方法找到特征点位置，确定超声波信号在顺流、逆流中的传播时间。由于双声道气体超声流量计计算一次瞬时流量需要确定每个声道的顺流、逆流传播时间，因此，循环完成四次上述过程后，可确定两个声道的顺流、逆流传播时间，完成一次瞬时流量的计算，再根据计算出的瞬时流量进一步计算出累积流量。

图4所示为定时器中断服务程序。定时器的时间计数为1s，当定时器计时达到1s时，DSP会进入该中断程序，读取瞬时流量，计算累积流量，并完成PWM输出；同时，将定时器的标志位置为1，此处所述标志位与图3中所述的标志位相同。

图5～图7是基于回波上升段峰值拟合的气体超声波流量计信号处理方法计算步骤的相关说明。

图5所示为不同流量下的回波信号及其包络线，通过该图可以确定拟合直线的大致区间。该图是将不同流量(0m³/h、200m³/h、400m³/h、600m³/h、800m³/h)下的各50组回波信号进行滤波、求平均和归一化预处理，利用MATLAB中曲线拟合工具对各流量下回波信号的所有峰值点进行包络线拟合得到的。由于回波信号采样点的位置与时间成对应关系，因此，拟合得到的连续曲线的横轴可看作时间。有研究表明，当被测流体流速超过一定范围时，回波信号会产生严重的畸变(SABATINI A M.A digital-signal-processing technique forultrasonic signal modeling and classification[J].IEEE Transactions onInstrumentation Measurement,2001,50(1):15-21.)。从该图中可以明显看出，随着流量的增大，回波信号的下降阶段确实发生畸变，但是，回波信号的上升阶段，由于受激励换能器的发射声波直接作用，基本不存在反射声波等影响。因此，上升阶段的波形轮廓基本不变。同时，在回波信号上升阶段的中间部分，由于其幅值变化率基本相同，该区间内有较多的峰值点分布在一条直线上。

图6所示为上升阶段包络线及其最大斜率切线，通过该图可以确定拟合直线的具体范围。该图是利用MATLAB中曲线拟合工具对零流量下50组滤波、归一化、求平均后的回波信号上升段的所有峰值点进行了包络曲线拟合得到的。比较不同阶次拟合得到曲线与实际数据的误差平方和(SSE)，得出采用5阶模型拟合时效果最佳。

选择靠近回波信号最大峰值附近若干个峰值点进行直线拟合时，由于最大峰值附近的幅值变化率较低，拟合直线的斜率较小，任一峰值点出现很小的波动，都会使拟合直线与横轴采样点(对应于时间)的交点(即特征点)发生较大波动，影响测量精度，因此，拟合直线的斜率越大，相对而言，任一峰值点出现波动时，拟合直线与横轴采样点(对应于时间)交点(即特征点)的波动会越小。选择分布在斜率最大直线上的峰值点进行直线拟合时，能有效地减少特征点波动，提高测量精度。

因此，根据采样信号的峰值点，求得上升段包络曲线S(x)。

S(x)＝-2.1924e^-13x⁵+5.7502e^-10x⁴-5.0412e^-7x³+1.5587e^-4x²

-6.6326e^-3x+0.07604) (4)

其中，由于回波信号采样点的位置与时间成对应关系，因此，拟合得到的连续曲线的横坐标x可看作时间。

为了得到上升阶段轮廓曲线的变化率，对曲线S(x)进行一次求导，得

S(x)′＝-1.0962e^-12x⁴+2.3001e^-9x³-1.5124e^-6x²+3.1174e^-4x-6.6326 (5)

为得到上升阶段轮廓曲线的最大变化率，即斜率最大的切线，通过MATLAB中的MAX函数，找出上升阶段轮廓曲线一次导数S(x)′的最大值，得到曲线S(x)对应的斜率最大的切线Y(x)。

Y(x)＝6.5987e^-3x-0.3248 (6)

为了得到分布在斜率最大的直线上的峰值点的范围，设定差值参数δ，根据式(4)、(6)、(7)得出对应x的范围。

|S(x)-Y(x)|≤δ (7)

具体的步骤为：

首先令S(x)-Y(x)＝δ，用MATLAB解方程，得出对应上升阶段轮廓曲线上的阈值下限横坐标X_min。然后，令S(x)-Y(x)＝-δ，用MATLAB解方程，得出对应上升阶段轮廓曲线上的阈值下限横坐标X_max。再将上、下限横坐标X_min、X_max代入公式(6)，可得出y的范围为Y_min～Y_max，即为拟合直线峰值点的阈值范围。

为了降低单个峰值波动对拟合直线带来的影响，增大确定特征点的稳定性，阈值的范围不宜过小。并且，回波信号峰值较小、峰值较大处相邻峰值差距较小，幅值变化率较小，容易受噪声影响，因此，阈值的范围应选择回波信号上升阶段的中间部分。综上所述，差值参数δ应折中选取。经过实验发现，将差值参数δ设置为0.01时，能得到较好的实验结果。

图7所示为基于回波上升段峰值拟合的特征点确定方法的示意图。根据该方法能确定超声波信号顺流、逆流的传播时间，进而计算气体流量。该方法计算简单，能保证系统具有较高的实时性；同时，能保证气体超声波流量计系统满足较高的测量精度，并且有效地扩大了气体超声波流量计的测量范围。具体的步骤为：

第一步：采用四阶巴特沃斯带通滤波器构成零相位滤波器，对采集的一组回波信号数据进行实时的带通滤波处理，通带范围为120kHz～280kHz。所述的零相位滤波器，即是将回波信号数据顺序、反序通过四阶巴特沃斯带通滤波器，得到滤波前后无相位差的回波信号。

第二步：根据确定的阈值范围Y_min～Y_max，找出回波信号上升阶段的N个峰值点。若设第k(k＝1,2,…,N)个峰值点坐标为(x_k,y_k)，则可计算出拟合直线斜率参数A。

再计算出拟合直线截距参数B。

式中，分别表示N个峰值点横坐标和纵坐标的平均值。

得到拟合直线y＝Ax+B。

第三步：得到该直线与横轴采样点(对应于时间)的交点该交点即为回波信号的特征点，根据该特征点即可计算超声波的顺流、逆流传播时间，以顺流传播时间T_d为例。

T_d＝T-t′ (12)

式中，T为有流量情况下激励信号起始时刻到特征点对应时刻的时间长度，t'为一个根据零流量下的传播时间求得的固定偏差值，其关系如图8所示。由于零流量下超声波的传播速度V_s不受流体的影响，是固定不变的；而超声波传播声道的长度L是已知的,因此，可直接计算出零流量下的顺流传播时间T_d|_Flow＝0。

再实测得到零流量下激励信号的起始时刻到特征点对应时刻的时间长度T|_Flow＝0，而零流量下激励信号的起始时刻、回波信号起始时刻和特征点对应时刻之间的关系，如图9所示，则t′为

t′＝T|_Flow＝0-T_d|_Flow＝0 (14)

重复步骤一至步骤三共4次，就完成了双声道顺流和逆流传播时间的测量，根据公式(15)，即可计算出两个声道的瞬时流量Q。

式中，为声道上的平均速度，S为管道的横截面积，α表示一个与声道数量、声道分布及雷诺数等因素有关的修正系数，L表示声道长度，θ表示换能器与气体管道之间夹角，T_d、T_u分别表示同一声道的超声波顺流、逆流传播时间。由于双声道是对称分布的，因此，两个声道计算的流量值权重各为0.5，则最终的瞬时流量Q′为

Q′＝0.5*Q₁+0.5*Q₂ (16)

式中，Q₁、Q₂分别表示两个不同声道，根据公式(15)计算出的瞬时流量。

在得到最终的瞬时流量后，根据如图4所示的DSP定时器的中断服务程序，每隔1秒钟计算一次累积流量，并完成相应显示、输出和通讯等功能。

表1所示为采用基于回波上升段峰值拟合信号处理方法的气体超声波流量计在安徽省计量科学研究院进行标定实验的结果。超声流量计检定规程(JJG1030—2007.超声流量计，中华人民共和国国家计量检定规程[S].国家质量监督检验检疫总局,2007.)要求1级精度气体超声波流量计流量范围大于3m/s(约85m³/h)时，平均实际脉冲系数误差小于±1％，重复性小于0.2％；流量范围小于3m/s(约85m³/h)时，平均实际脉冲系数误差小于±2％，重复性小于0.4％。从表1可知，最大平均实际脉冲系数误差为-0.349％，最大重复性为0.16％。实验结果表明，采用基于回波上升段峰值拟合信号处理方法的气体超声波流量计的测量精度满足超声流量计检定规程对1级精度气体超声波流量计的相关要求。同时，由于大流量时，回波信号衰减严重，容易受噪声干扰，使得气体超声波流量计可测的最大流量受限，处于国际领先地位的Daniel等公司的同类仪表的最大流量为850m³/h。而基于回波上升段峰值拟合的气体超声波流量计信号处理方法不仅计算简单，能保证气体超声波流量计系统具有较高的实时性，并且，其最大可测流量为1000m³/h，有效地提高了气体超声波流量计可测的最大流量。

表1标定结果

图10～图14是基于回波能量点定位的气体超声波流量计信号处理方法的计算步骤的相关说明。

图10所示为不同流量下的回波信号能量及其包络线。通过该图可以看出不同流量下回波信号上升阶段能量包络线是相同且平移的。该图是用不同流量(0m³/h、200m³/h、400m³/h、600m³/h、800m³/h、1000m³/h)下的各50组回波信号进行滤波、求平均、归一化预处理，然后，对回波信号所有采样点求平方，并找出了平方后的所有峰值能量点，再利用MATLAB中的envelope函数求其包络，得到的回波信号的能量包络线。由于信号的能量与信号的幅值平方呈比例关系，因此，为了简化计算，本发明使用回波信号幅值的平方表征回波信号的能量。同时，又由于回波信号采样点的位置与时间成对应关系，因此，拟合得到的连续曲线的横轴可看作时间。从图10中可以明显看出，不同流量下的回波信号上升阶段能量包络线基本相同，并且随着流量的增大，回波信号上升阶段能量包络线是等比例平移的。而回波信号能量包络线的下降阶段，由于受到反射波等影响，随着流量的增大逐渐产生变化。由于不同流量下回波信号上升阶段能量包络线是相同且平移的，所以，可以将回波信号上升阶段能量包络线上某一固定能量点作为特征点，用于计算超声波传播时间。

图11所示为不同回波能量包络线的比较，通过该图可以直观的说明易回波能量包络线靠近零和靠近峰值部分容易产生波动。该图是从0m³/h流量下的多组回波信号中，选择两组具有代表性的回波能量数据，并求其上升段能量包络线得到的。从图11中可以明显看出，回波能量包络线1和包络线2在靠近零和靠近峰值部分存在波动，靠近中间部分基本重合，而靠近零和靠近峰值部分能量包络线的变化率较小，靠近中间部分能量包络线的变化率较大。因此，说明选择能量包络线变化率大的能量点作为特征点，能有效地减少特征点的波动，提高特征点的稳定性。

图12所示为回波信号上升阶段能量包络线及能量包络线变化率曲线，通过该图可以找到回波信号上升阶段能量包络线上，能量包络线变化率大的能量点的范围，即为稳定性较高的能量点的范围。该图是利用MATLAB中曲线拟合工具对零流量下50组滤波、归一化、求平均后的回波信号上升段的所有平方后的峰值点进行了包络曲线拟合得到的。比较不同阶次拟合得到曲线与实际数据的误差平方和(SSE)，得出采用6阶模型拟合时效果最佳。拟合得到的回波信号上升阶段能量曲线H(x)为

H(x)＝-4.3125e^-15x⁶+3.5638e^-11x⁵-1.0899e^-7x⁴+1.6605e^-4x³-0.1360x²+57.4774x-9.8681 (17)其中，由于回波信号采样点的位置与时间成对应关系，因此，拟合得到的连续曲线的横坐标x可看作时间。

对该回波信号能量包络线进行一次求导，得到回波信号上升阶段能量包络线变化率曲线H(x)′为

H(x)′＝-2.5875e^-14x⁵+1.7819e^-10x⁴-4.3596e^-7x³+4.9816e^-4x²-0.2720x+57.4774(18)

可得能量包络变化率最大的能量点的能量值为0.53。选择最大能量变化率的90％作为标准，得到一个能量变化率较大的区间，即可得到能量包络变化率较大的能量值范围。选定变化率较大的能量值范围内的某一能量点作为特征点，通过该点即可计算出传播时间。

图13所示为不同流量下回波能量峰值点比较。该图说明峰值点位置随着流量的变化会产生偏移，直接根据某一能量值，找到不同流量下对应的特征点(图中用“*”号区别该特征点与实际处理方法得到的特征点)，存在较大偏差。该图是用不同流量(0m³/h、400m³/h、800m³/h、1000m³/h)下的各50组回波信号进行滤波、归一化和求平均后得到的。由于回波信号类似一个幅值受调制的连续正弦信号，存在上下振荡，从图13中可以看出，随着流量的增大，回波信号并不是完全平移变化的，其同一峰值点会发生偏移，导致其峰值能量点也会发生相应的偏移。但是，根据图10中对回波信号能量包络线的分析可知，随着流量的增大，虽然回波峰值能量点会随着流量发生偏移，但是，根据回波峰值能量点拟合的上升段能量包络线是相同且平移的。因此，可以将拟合的上升段能量包络线上的某一能量点作为特征点，计算传播时间。

图14所示为基于回波能量点定位的特征点确定方法的示意图。由于回波信号是通过某一基准能量点进行定位的，因此，在实际测量流量时，只用拟合该能量点附近的能量包络曲线即可。同时，气体超声流量计要求系统具有较高的实时性，采用曲线拟合会大大增加计算量。而通过分析回波信号能量包络曲线以及能量包络变化率曲线，不难发现，中间部分能量点的能量变化率基本相同，说明该区间内的能量点基本分布在一条直线上。

为了减少信号处理的计算量，进一步提高系统的实时性，同时，又避免单个峰值能量点的波动对拟合后能量点定位的影响，因此，采用能量变化率较大范围内的任一能量点的能量值作为基准，找到四个最邻近的峰值能量点进行直线拟合，再找出拟合直线上该能量值对应的能量点作为特征点，用于计算超声波传播时间，就能得到较为准确的结果。如图12所示，选择最大变化率能量点的能量值0.53为基准，由于采样值受采样率的和量化误差的影响，直接将采样点位置和平方后的采样值作为峰值能量点的坐标，往往不够准确的。而回波能量数据的峰值部分满足抛物线变化规律，因此，对平方后的采样峰值点和该点附近的两个平方后的采样点(共三个点)，进行二次拟合，得到二次曲线顶点坐标，作为峰值能量点的坐标，进一步保证拟合的精度。

确定特征点的具体步骤为：

第一步：采用四阶巴特沃斯带通滤波器构成零相位滤波器，对采集的一组回波信号数据进行实时的带通滤波处理，通带范围为120kHz～280kHz。所述的零相位滤波器，即是将回波信号数据顺序、反序通过四阶巴特沃斯带通滤波器，得到滤波前后无相位差的回波数据。再对回波信号求平方，得到其回波能量数据。

第二步：以最大变化率能量点的能量值0.53为基准，找到离该值最邻近的四个能量峰值点及每个峰值点附近的两个能量点。

第三步：对四个能量峰值点及其附近的两个能量点进行二次拟合，得到二次曲线顶点坐标，作为最终的峰值能量点坐标。以计算其中一个峰值能量点为例，具体的步骤为：

(1)将平方后的采样峰值点的坐标记为(x^k,y^k)，x表示采样点的位置，y表示平方后的采样值。峰值点附近的两个平方后的采样点的坐标，记为(x^k-1,y^k-1)、(x^k+1,y^k+1)(其中有x^k-x^k-1＝1，x^k+1-x^k＝1)。

(2)对这三个点进行二次函数拟合，得到二次函数曲线方程y′＝ax′²+bx′+c的相关参数。

b＝y^k-y^k-1-(2x^k-1+1)*a (20)

c＝y^k-a(x^k)²-bx^k (21)

(3)根据式(22)和式(23)得到二次拟合曲线的顶点坐标(x,y)，将该顶点作为最终的峰值能量点。

由于可直接通过公式(19)-(23)得到顶点坐标，并且，处理一组回波信号只用调用四次上述公式，因此，该步骤并不占用过多的计算量。

第四步：对第三步求得的四个峰值能量点进行直线拟合，找到最大变化率能量点的能量值0.53在拟合直线上对应的能量点作为特征点。根据该点的横坐标，确定传播时间。具体的步骤为：

(1)设最终的四个峰值能量点的坐标为(x_n,y_n)(n＝1,2,3,4)，可得拟合直线y＝Kx+B，K为拟合直线斜率参数，B为拟合直线截距参数。

式中，分别表示4个峰值能量点横坐标和纵坐标的平均值。

(2)令Kx+B＝0.53，即可得到拟合直线上能量值0.53对应的能量点的坐标将该点作为特征点。根据该特征点即可计算超声波的顺流、逆流传播时间。以顺流传播时间T_d为例。

T_d＝T-t′ (26)

式中，T为有流量情况下激励信号起始时刻到特征点对应时刻的时间长度，t'为一个根据零流量下的传播时间求得的固定偏差值，其关系如图8所示。由于零流量下超声波的传播速度V_s不受流体的影响，是固定不变的；而超声波传播声道的长度L是已知的,因此，可直接计算出零流量下的顺流传播时间T_d|_Flow＝0。

再实测得到零流量下激励信号的起始时刻到特征点对应时刻的时间长度T|_Flow＝0，而零流量下激励信号的起始时刻、回波信号起始时刻和特征点对应时刻之间的关系，如图9所示，则t′为

t′＝T|_Flow＝0-T_d|_Flow＝0 (28)

重复步骤一至步骤四共4次，就完成了双声道顺流和逆流传播时间的测量，根据公式(29)，即可计算出两个声道的瞬时流量Q。

式中，为声道上的平均速度，S为管道的横截面积，α表示一个与声道数量、声道分布及雷诺数等因素有关的修正系数，L表示声道长度，θ表示换能器与气体管道之间夹角，T_d、T_u分别表示同一声道的超声波顺流、逆流传播时间。由于双声道是对称分布的，因此，两个声道计算的流量值权重各为0.5，则最终的瞬时流量Q′为

Q′＝0.5*Q₁+0.5*Q₂ (30)

式中，Q₁、Q₂分别表示两个不同声道，根据公式(29)计算出的瞬时流量。

在得到最终的瞬时流量后，根据如图4所示的DSP定时器的中断服务程序，每隔1秒钟计算一次累积流量，并完成相应显示、输出、通讯等功能。

表2所示为采用基于回波能量点定位信号处理方法的气体超声流量计在安徽省计量科学研究院进行标定实验的结果。超声流量计检定规程(JJG1030—2007.超声流量计，中华人民共和国国家计量检定规程[S].国家质量监督检验检疫总局,2007.)要求1级精度气体超声流量计流量范围大于3m/s(约85m³/h)时，平均实际脉冲系数误差小于±1％，重复性小于0.2％；流量范围不大于3m/s(约85m³/h)时，平均实际脉冲系数误差小于±2％，重复性小于0.4％。从表1可知，流量范围大于3m/s(约85m³/h)时，最大平均实际脉冲系数误差为+0.772％，最大重复性为0.11％；流量范围不大于3m/s(约85m³/h)时，最大平均实际脉冲系数误差为+1.919％，最大重复性为0.21％。实验结果表明，采用基于回波能量点定位信号处理方法的气体超声流量计的测量精度满足超声流量计检定规程对1级精度气体超声流量计的相关要求。同时，由于大流量时，回波信号衰减严重，容易受噪声干扰，使得气体超声流量计可测的最大流量受限，处于国际领先地位的Daniel等公司的同类仪表的最大流量为850m³/h。而基于回波能量点定位的气体超声流量计信号处理方法不仅计算简单，能保证气体超声流量计系统具有较高的实时性，并且，其最大可测流量为1300m³/h，比同类仪表足足提高到约1.53倍。

表2标定结果

Claims (2)

1.基于回波上升段峰值拟合的气体超声波流量计信号处理方法，通过分析回波信号包络线的变化规律，建立回波信号上升段包络曲线的数学模型，得出拟合峰值点的范围，找到一个稳定的特征点，确定超声波信号顺流、逆流的传播时间，进而计算气体流量；其特征在于：

该方法的具体的步骤为：

第一步：采用四阶巴特沃斯带通滤波器构成零相位滤波器，对采集的一组回波信号数据进行实时的带通滤波处理，通带范围为120kHz～280kHz；所述的零相位滤波器，即是将回波信号数据顺序、反序通过四阶巴特沃斯带通滤波器，得到滤波前后无相位差的回波信号；

第二步：根据确定的阈值范围Y_min～Y_max，找出回波信号上升阶段的N个峰值点；若设第k(k＝1,2,…,N)个峰值点坐标为(x_k,y_k)，则可计算出拟合直线斜率参数A；

再计算出拟合直线截距参数B；

式中，分别表示N个峰值点横坐标和纵坐标的平均值；

得到拟合直线y＝Ax+B；

第三步：得到该直线与横轴采样点(对应于时间)的交点该交点即为回波信号的特征点，根据该特征点即可计算超声波的顺流、逆流传播时间，进而计算气体流量。

2.基于回波能量点定位的气体超声波流量计信号处理方法，通过分析回波信号能量包络线的变化规律，建立回波信号上升阶段能量曲线的数学模型，得出能量包络变化率较大的能量值范围，找到一个稳定的特征点，确定超声波信号顺流、逆流的传播时间，进而计算气体流量；其特征在于：

该方法的具体的步骤为：

第一步：采用四阶巴特沃斯带通滤波器构成零相位滤波器，对采集的一组回波信号数据进行实时的带通滤波处理，通带范围为120kHz～280kHz；所述的零相位滤波器，即是将回波信号数据顺序、反序通过四阶巴特沃斯带通滤波器，得到滤波前后无相位差的回波数据；再对回波信号求平方，得到其回波能量数据；

第二步：以最大变化率能量点的能量值0.53为基准，找到离该值最邻近的四个能量峰值点及每个峰值点附近的两个能量点；

第三步：对四个能量峰值点及其附近的两个能量点进行二次拟合，得到二次曲线顶点坐标，作为最终的峰值能量点坐标；以计算其中一个峰值能量点为例，具体的步骤为：

(1)将平方后的采样峰值点的坐标记为(x^k,y^k)，x表示采样点的位置，y表示平方后的采样值；峰值点附近的两个平方后的采样点的坐标，记为(x^k-1,y^k-1)、(x^k+1,y^k+1)(其中有x^k-x^{k -1}＝1，x^k+1-x^k＝1)；

(2)对这三个点进行二次函数拟合，得到二次函数曲线方程y′＝ax′²+bx′+c的相关参数；

b＝y^k-y^k-1-(2x^k-1+1)*a

c＝y^k-a(x^k)²-bx^k

(3)根据下式可得到二次拟合曲线的顶点坐标(x,y)，将该顶点作为最终的峰值能量点；

第四步：对第三步求得的四个峰值能量点进行直线拟合，找到最大变化率能量点的能量值0.53在拟合直线上对应的能量点作为特征点；根据该点的横坐标，确定传播时间；具体的步骤为：

(1)设最终的四个峰值能量点的坐标为(x_n,y_n)(n＝1,2,3,4)，可得拟合直线y＝Kx+B，K为拟合直线斜率参数，B为拟合直线截距参数；

式中，分别表示4个峰值能量点横坐标和纵坐标的平均值；

(2)令Kx+B＝0.53，即可得到拟合直线上能量值0.53对应的能量点的坐标将该点作为特征点；根据该特征点即可计算超声波的顺流、逆流传播时间，进而计算气体流量。