一类基于回波信号包络拟合的气体超声波流量计信号处理方法与流程

本发明属于流量检测领域，为气体超声波流量计信号处理方法，特别是一类基于回波信号包络拟合的气体超声波流量计信号处理方法，具体包括基于回波信号上包络拟合、基于回波信号上、下包络拟合和改进的基于回波信号上、下包络拟合的气体超声波流量计信号处理方法，以实现气体流量的准确测量。

背景技术：

气体超声波流量计具有量程比大、测量精度高、无压损等诸多优点，特别是在大口径天然气流量测量方面，具备独特的优势。气体超声波流量计的测量原理分为传播速度差法和多普勒法等。传播时间差法又可以分为时差法、相差法和频差法，其中，时差法效果好，在气体超声波流量计中应用最为广泛。基于时差法测量原理的气体超声波流量计在测量气体流量时，首先依据逆压电效应，激励一个超声波换能器发射超声波信号；另一个超声波换能器接收到超声波信号，依据压电效应，转换为回波电信号；根据回波信号中某个稳定的特征点确定超声波信号顺、逆流传播时间，进而计算气体流量。但是，超声波信号在气体介质中传播能量衰减严重，回波信号存在幅值微弱、信噪比低和易受干扰等问题，随着流量增大，这些问题显得尤为严重。并且，直接利用寻找单个峰值点等简单方法确定特征点，容易产生较大波动，影响系统测量精度。

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

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

美国daniel公司采用检测能量突变的方法用于寻找回波信号特征点((williamfreund,winsorletton,jamesmc-clellan,baocangjia,anniwey,wenchang.methodandapparatusformeasuringthetimeofflightofasignal,uspatentno.5983730,nov.16,1999)。回波信号能量变化为一个从弱变强，再从强变弱的过程，相应的能量变化率的变化为先增大后减小，找到变化率变化的临界点，作为特征点，即可确定超声波信号的传播时间。具体步骤为：首先，求取回波信号各点幅值的平方，由于信号各点幅值的平方与信号能量呈线性关系，可表征信号能量大小。然后，利用滑动平均的方法，求取信号各点的平均能量，绘制信号前后两点平均能量之比变化曲线，即可得到能量变化曲线。最后，通过寻找能量变化率变化的临界点，定位回波信号，计算超声波信号的传播时间。但是，该专利中并未披露处理方法的关键参数。

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

从能量变化的角度出发，沈子文等也提出一种基于能量变化率的信号处理方法(沈子文,徐科军,方敏,等.基于能量变化率的气体超声波流量计信号处理方法[j].仪器仪表学报,2015,36(9):2138-2144.)。具体步骤为：首先，对回波信号进行滤波，对滤波后回波信号各峰值点幅值进行平方作为该峰值点的能量，通过前一个峰值点能量与后面峰值进行比较，得到各峰值点的能量变化率，进行包络拟合，得到不同流量下回波信号能量变化曲线；然后，根据能量变化曲线设定阈值，寻找特征波。此方法在实时测量时，需要实时求取回波信号的能量变化曲线，根据预先设定阈值寻找特征点，计算量较大，难以保证系统实时性。并且，需要在流量测量之前进行阈值参数选定，流程较复杂。

(3)基于回波上升段峰值拟合的方法

从轮廓包络的角度出发，田雷等提出了一种基于回波上升段峰值拟合的信号处理方法(田雷,徐科军,沈子文,等.基于回波峰值拟合的气体超声流量计信号处理[j].电子测量与仪器学报,2017,31(7):1107-1114.)。具体步骤为：首先，对回波信号上升段峰值进行包络拟合，得出回波信号上升段包络近似为一条直线；然后，确定处于一条直线的包络部分覆盖的峰值点幅值的范围；对幅值范围内的各峰值点进行最小二乘拟合，得到拟合直线；选取拟合直线与横轴采样点的交点作为特征点，定位回波信号，确定超声波信号顺、逆流传播时间，进而计算气体流量。此方法能快速、较准确地定位回波信号。但是，气体流量较小时，采样到的拟合峰值点存在较小的波动，拟合峰值点附近直线斜率会产生较小的变动。这些较小的变动在直线远离拟合峰值点与横轴采样点相交的部分被放大，导致小流量测量结果不太理想。

技术实现要素：

针对气体超声波流量计难以确定回波信号到达时刻、并且仪表对实时性要求较高的问题，本发明提出一类基于回波信号包络拟合的气体超声波流量计信号处理方法，具体包括基于回波信号上包络拟合的气体超声波流量计信号处理方法、基于回波信号上、下包络拟合的气体超声波流量计信号处理方法和改进的基于回波信号上、下包络拟合的气体超声波流量计信号处理方法，并在合肥工业大学研制的双声道气体超声波流量计硬件平台上(徐科军，方敏，汪伟，朱文姣，沈子文。一种基于fpga(现场可编程门阵列)和dsp(数字信号处理器)的气体超声流量计，申请中国发明专利，公开号cn104697593a，申请日2015.06.10)，实现这一类基于回波信号包络拟合的气体超声波流量计信号处理方法。

基于回波信号上包络拟合的气体超声波流量计信号处理方法，首先以上包络上升段定位点的幅值作为基准，对距离该点最近的三个峰值点进行直线拟合，得到一条特征直线；选取特征直线与x轴的交点作为特征点，确定超声波信号顺、逆流传播时间，计算气体流量。

基于回波信号上、下包络拟合的气体超声波流量计信号处理方法，首先以上包络上升段定位点的幅值作为基准，对距离该点最近的三个峰值点进行直线拟合，得到一条特征直线；然后，以下包络下降段定位点的幅值作为基准，对距离该点最近的三个峰值点进行直线拟合，得到第二条特征直线；选取两条特征直线的交点作为特征点，确定超声波信号顺、逆流传播时间，计算气体流量。

改进的基于回波信号上、下包络拟合的气体超声波流量计信号处理方法，首先以上包络上升段定位点的幅值作为基准，对距离该点最近的三个峰值点进行直线拟合，得到一条特征直线；然后，以下包络下降段定位点的幅值作为基准，对距离该点最近的三个峰值点进行直线拟合，得到第二条特征直线；将下包络拟合直线向上平移一定距离，选取平移过后两条特征直线的交点作为特征点，确定超声波信号顺、逆流传播时间，计算气体流量。

附图说明

图1是不同流量下的回波信号及其上、下包络。

图2是0m³/h流量下回波信号上包络上升段、下包络下降段以及两包络线的变化率曲线。图3是基于回波信号上包络拟合寻找特征点的信号处理方法示意图。

图4是0m³/h流量下50组回波信号上包络上升段拟合直线以及特征点的分布情况。

图5是基于回波信号上、下包络拟合寻找特征点的信号处理方法示意图。

图6是0m³/h流量下50组回波信号上、下包络拟合直线以及特征点的分布情况。

图7是改进的基于回波信号上、下包络拟合寻找特征点的信号处理方法示意图。

图8是0m³/h流量下50组回波信号上包络拟合直线、平移后下包络拟合直线以及特征点的分布情况。

图9是0m³/h流量下运用三种方法寻找的特征点的分布情况。

图10是不同流量下运用三种方法寻找的特征点的方差变化曲线。

图11是双声道气体超声波流量计系统的硬件框图。

图12是dsp软件流程图。

图13是改进的基于回波信号上、下包络拟合的信号处理方法的处理流程图。

具体实施方式

图1为不同流量下的回波信号及其上、下包络，可见回波信号上包络上升段和下包络的下降段中间部分的若干峰值点分别位于一条直线上。对0m³/h、200m³/h、400m³/h、600m³/h、800m³/h、1000m³/h流量下50组回波信号进行滤波、求和取均值、归一化处理，并用matlab对信号正、负峰值点分别进行曲线拟合，得到回波信号上、下包络，x轴为采样点与时间相对应。回波信号上包络下降段和下包络上升段随着流量增大，出现较大的畸变，而回波信号上包络上升段和下包络下降段的各峰值点受超声波换能器发射的激励信号直接作用，不存在反射波等因素的影响，峰值点幅值几乎不变。因此，回波信号上包络上升段和下包络下降段随着流量增大，轮廓几乎不变，近似于沿x轴平移。同时，回波信号上包络上升段和下包络下降段的中间部分的变化率基本相同，且随流量变化基本不变。该区间内的若干峰值点分布一条直线上。

图2为0m³/h流量下归一化后的回波信号上包络上升段、下包络下降段以及两包络的变化率曲线，通过该图可以确定上、下包络用于拟合直线的峰值点的范围。对0m³/h流量下50组回波信号进行滤波、求和取均值、归一化处理，运用matlab对信号正、负峰值点分别进行曲线拟合，得到效果较理想的6阶模型。

首先，对回波信号正峰值点进行拟合，得到回波信号上包络上升段包络。

对s1(x)进行求导，得到回波信号上包络上升段包络变化率曲线。

利用matlab中的max(最大值)函数求得包络变化率最大点的横坐标为877，带入式(1)得到包络线上变化率最大点的坐标为(877，0.52)。最大变化率90％范围内的包络线部分近似分布在一条直线上，可用一条直线作为特征直线来代替拟合曲线。用于拟合直线的若干峰值点为最大变化率90％范围内包络线部分覆盖的三个峰值点，是距离上包络上升段包络变化率最大点最近的三个峰值点。

同理，对回波信号负峰值点进行拟合，得到回波信号下包络下降段包络。

对s2(x)进行求导，得到回波信号下包络下降段包络变化率曲线。

利用matlab中的min(最小值)函数得到包络变化率最小点的横坐标为882，带入式(3)得到包络线上变化率最小点的坐标为(882,-0.56)。最小变化率90％范围内的包络线部分近似分布在一条直线上，可用一条直线作为特征直线来代替拟合曲线。用于拟合直线的若干峰值点为最小变化率90％范围内包络线部分覆盖的三个峰值点，是距离下包络下降段包络变化率最小点最近的三个峰值点。

本发明提出一类基于回波信号包络拟合寻找特征点的信号处理方法，利用两条拟合直线(特征直线)彼此相交且与x轴相交的特性，分别选取特征直线与x轴交点、两条特征直线的交点作为特征点，定位回波信号。

图3为一类基于回波信号包络拟合寻找特征点的信号处理方法中的方法一--基于回波信号上包络拟合寻找特征点的信号处理方法的示意图。利用拟合直线(特征直线)与x轴稳定相交的特性，选取交点作为特征点，可快速地定位回波信号。

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

(1)采用四阶巴特沃斯带通滤波器对采集到的一组回波信号进行滤波处理，通带范围为120khz～280khz，得到滤波后的回波信号。

(2)选取上包络变化率最大点归一化后的幅值0.52为基准，寻找距离该点最近的三个峰值点a、b、c，坐标记为(xn，yn)(n＝1，2，3)。对峰值点a、b、c进行直线拟合，可得特征直线y1＝k1x+b1，k1为直线的斜率，b1为直线的截距。

式中，为三个峰值点a、b、c横坐标的平均值，为三个峰值点a、b、c纵坐标的平均值。

(3)计算特征直线与x轴的交点即为特征点，进而确定超声波信号的传播时间，计算气体流量。

图4为0m³/h流量下50组回波信号上包络上升段分别进行直线拟合形成的直线簇，50条特征直线与x轴相交形成的交点集，分布较稀疏，特征点较为稳定。拟合峰值点附近直线簇的斜率稳定，分布密集，与x轴的交点(特征点)分布范围较小，但是，存在波动。这种寻找特征点的方法计算量较小，同时能较为准确的找到特征点。直线与x轴相交部分距离拟合直线的3个峰值点较远，峰值点产生细小的波动导致直线斜率出现较小偏差。这种偏差在直线离峰值点较远的部分被放大，直观上表现为直线簇与x轴的交点分布存在波动。

从图2可见，回波信号上包络上升段和下包络下降段的中间部分分别位于一条直线上，并且两直线相交。若取两直线交点作为特征点，则可以提取更多地流量信息。同时，正、负峰值点动态变化，产生相互补偿的作用，特征点变得更加稳定。

图5为一类基于回波信号包络拟合寻找特征点的信号处理方法中的方法二--基于回波信号上、下包络拟合寻找特征点的信号处理方法的示意图。利用两条拟合直线(特征直线)稳定相交的特性，选取交点作为特征点，可以较快速地定位回波信号。

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

(1)采用四阶巴特沃斯带通滤波器对采集到的一组回波信号进行滤波处理，通带范围为120khz～280khz，得到滤波后的回波信号。

(2)前一种方法已对上包络上升段距离上包络变化率最大点最近的三个峰值点a、b、c进行直线拟合，得到第一条特征直线y1＝k1x+b1，k1为直线的斜率，b1为直线的截距。同理，以下包络下降段包络变化率最小的点归一化后的幅值-0.56为基准，寻找距离该点最近的三个峰值点d、e、f，坐标记为(xn，yn)(n＝4，5，6)。对三个峰值点进行直线拟合，可得第二条特征直线y2＝k2x+b2，k2为直线的斜率，b2为直线的截距。

式中，为三个峰值点d、e、f横坐标的平均值，为三个峰值点d、e、f纵坐标的平均值。

(3)令y1＝y2，即k1x+b1＝k2x+b2，得到两条特征直线的交点m2即为特征点，进而确定超声波信号的传播时间。

图6为0m³/h流量下50组回波信号上、下包络分别进行直线拟合得到一对直线簇，每组回波信号的两条特征直线有一个交点共50个，形成一个交点集，分布较稀疏。该图中峰值点a、b、c、d、e、f附近直线分布密集，相交部分直线分布较为稀疏，直线交点(特征点)的分布范围较小，但是，同样存在波动。这种寻找特征点的方法包含了更多地流量信息，计算量较小，特征点较为稳定；两条特征直线相交部分距离拟合直线的峰值点较远，正、负峰值点波动的相互补偿作用微弱，特征点依然存在波动，相较于前一种方法提升不大。

由于拟合直线的峰值点a、b、c、d、e、f较为稳定，峰值点附近直线分布密集，距离峰值点较远的部分直线分布较为稀疏，图6中两直线簇的分布情况可证明此结论，所以，若两直线簇相交的部分为两直线簇分布较为密集部分，则可以有效减小直线交点的波动，增强特征点的稳定性。

图7为一类基于回波信号包络拟合寻找特征点的信号处理方法中的方法三--改进的基于回波信号上、下包络拟合寻找特征点的信号处理方法的示意图。利用拟合直线峰值点附近直线分布密集的特性，将下包络拟合直线向上平移一定距离，选取平移过后两条特征直线的交点作为特征点，可以较快速地定位回波信号。

确定特征点的具体步骤：

(1)采用四阶巴特沃斯带通滤波器对采集到的一组回波信号进行滤波处理，通带范围为120khz～280khz，得到滤波后的回波信号。

(2)前两种方法已对上包络上升段距离上包络变化率最大点最近的三个峰值点a、b、c进行直线拟合，得到一条特征直线y1＝k1x+b1，k1为直线的斜率，b1为直线的截距；对下包络下降段距离下包络变化率最小点最近的三个峰值点d、e、f进行直线拟合，得到第二条特征直线y2＝k2x+b2，k2为直线的斜率，b2为直线的截距。

(3)将下包络下降段拟合直线向上平移l2，取两条特征直线的交点m3作为特征点，能较为快速地定位回波信号，进而计算信号传播时间。

向上平移的距离l2决定两直线相交的部分，由于上、下包络两拟合直线近似关于x轴对称，保证两线簇的相交部分为两线簇的密集部分，最小平移距离约为两倍的峰值点d幅值的绝对值，记为l1；最大平移距离约为两倍的峰值点c幅值的绝对值，记为l3。为了保证两线簇的密集部分稳定相交，取l1、l3中间部分，即为平移距离l2，约为1.1，。

图8为0m³/h流量下50组回波信号上包络上升段拟合直线簇和平移l2距离后的下包络下降段拟合直线簇，每组回波信号的两条特征直线有一个交点共50个，形成一个交点集。该图中两线簇的相交部分在峰值点a、b、c、d、e、f附近，直线交点(特征点)的分布范围小，分布密集。这种寻找特征点的方法虽然增加了计算量，但是，峰值点细小的波动导致较小的直线斜率偏差，不会对特征点产生较大影响，增强了特征点的稳定性，同时也能保证系统具备较好的实时性。

图9为运用上述三种方法寻找的直线交点(特征点)的分布情况，可明显地看出运用方法三(改进的基于回波信号上、下包络拟合的信号处理方法)寻找的特征点分布密集。运用方法一(基于回波信号上包络拟合的信号处理方法)寻找的特征点，分布区间较大，约为4.5个采样间隔，分布较稀疏；运用方法二(基于回波信号上、下包络拟合的信号处理方法)寻找的特征点，分布区间与分布情况与方法一的基本相同；运用方法三寻找的特征点，分布区间小，约为1.6个采样间隔，分布密集。

图10为在不同流量下运用三种方法寻找的特征点的方差变化曲线，方法三效果最好。方差反映数据波动程度，方差数值越大代表数据波动越大。运用三种方法对0m³/h、30m³/h、60m³/h、100m³/h、200m³/h、300m³/h、400m³/h、500m³/h、600m³/h、700m³/h、800m³/h、900m³/h、1000m³/h流量下50组回波信号进行处理，寻找特征点，利用matlab中std(标准偏差)函数求取各流量下的特征点的方差，将方差随流量增大的变化情况绘制在一张图上。方法三相较于其他两种信号处理方法而言，在流量较小的情况下稳定性较高，优势明显；方法二与方法一相比，由于正、负峰值点动态补偿作用，特征点方差有一定改善；随着流量增大，三种方法寻找的特征点的方差都呈增长趋势，并且趋于相同。

图11为双声道气体超声波流量计系统的硬件框图，主要包括一次仪表、dsp最小系统、fpga最小系统、驱动信号生成与放大电路、发射/接收信号通道切换电路、回波信号调理和采集电路、电源管理模块。

所述一次仪表由四个换能器和流量管组成，四个换能器分别固定安装在流量管上；四个换能器分别为换能器1、换能器2、换能器3和换能器4，组建双声道直射式结构；换能器为收发两用型，通过发射/接收信号通道切换电路控制，既可作为发射换能器用于发射激励信号，又可作为接收换能器用于接收回波信号。

所述fpga最小系统由fpga芯片、时钟、复位电路和flash(闪存)组成，其中，fpga芯片型号为ep2c8q208c8n。其中，fpga内部包括rom(只读存储器)模块、dac(数模转换器)控制模块、分频模块和ram_2port模块。

所述dsp最小系统由dsp芯片、串口通信、外部看门狗、时钟、fram(铁电存储器)和人机接口组成，其中，dsp芯片型号为tms320f28335。

所述电源管理模块为系统提供所需幅值的数字电源和模拟电源。

图12为dsp软件流程图。系统上电后，dsp进行初始化，完成对gpio口的分配、定时器中断初始化、液晶显示初始化、流量计显示初始化等初始化任务。dsp向fpga发送激励起始信号，并切换发射、接收通道。系统检测定时器时间标志位是否为0，若为0，则进入等待转存信号状态；若为1，则更新液晶显示，完成与上位机通讯，并将定时器时间标志位复位。若dsp收到转存起始信号后，开始读取fpga暂存的回波信号，对回波信号进行处理，确定超声波信号顺、逆流传播时间，计算气体流量。

图13为方法三的处理流程图。首先，寻找距离回波信号上包络上升段最大变化率点最近的三个峰值点，对其进行直线拟合，得到一条特征直线；然后，再寻找距离回波信号下包络下降段最小变化率点最近的三个峰值点，对其进行直线拟合，得到第二条特征直线；将下包络拟合直线向上平移一定距离，求得平移后两条特征直线的交点，即为特征点。