一种水声声速的时频测量综合算法研究

鲍晶晶， 蒋志迪， 刘尉悦

(1.宁波大学 信息科学与工程学院，浙江 宁波 315201；2.宁波大学 科学技术学院，浙江 宁波 315300)

1 引 言

随着超声波技术的不断发展，超声探伤和超声测距等技术已越来越成熟，被广泛应用到各工业领域[1]。在环境恶劣的工业生产环境，尤其是一些具有强氧化性或易挥发性溶液场合，如何实时监测环境参数，如液体浓度、液体流速与流量，成为高端自动化生产的必要条件。通过实时测量能够及时进行反馈操作，提高生产效率，减少生产成本等，较传统的测量方式更具优势。

目前采用超声波测量液体浓度、流速与流量的方法已有多种，如波束偏移法、相关法、相位差法、时差法和多普勒法等，其中相位差的测量应用最为广泛[2,3]。液体浓度测量一般侧重于静态环境测量，而流速与流量测量则偏重于动态环境测量。时差法超声波测量方法因其测量精度和结果可重复性等突出的优点受到越来越多研究者的青睐。

超声时差法用于水表计量，文献[4]应用高速驱动及高速切换电路，减少时间测量的误差以提高时差法的精度；而文献[5]通过高精度计时芯片TDC-GP22为核心获得传播时间差，从而提高水表的测量精度。对于单频信号的延迟时间差测量，也可采用文献[6]四参数拟合测量方法进行，但仅给出了解决方案。在时差法测量过程中，时差误差主要由信号传递时延不同、环境因素影响、噪声及超声波接收器器件等引起。由环境因素如温度等引起的时差校正已在很多文献中进行了讨论，如文献[7]，也有文献[8、9]利用超声回波的波速变化获取温度场分布情况。

诸多学者通过构建测量系统以期获得超声飞行时间测量的高精度。文献[10]研制了双反射脉冲回波技术测量液相及超临界流体声速的实验系统，获得较高的声速测量精度，但系统结构相对复杂；而文献[11]应用差分飞行时间法，设计研制了双超声腔体结构用于高压液体中声速的精密测量；文献[12]更是从换能器安装、换能器模型等多角度分析，结合阈值法、互相关法等提出一种温度自适应的超声波渡越时间测量方法，但仅从理论上进行分析并通过软件进行仿真来验证方法的可行性，并未给出软硬件系统实现结果。

为提高超声信号的时间测量精度，精确的时间单元或时间基准单元是首要条件，文献[13]搭建硬件系统，利用FPGA设计固定时延单元、运用延迟线法提高测时精度。文献[14]则构建粗延时单元及精延时单元获得回波信号的延时高精度测量。文献[15]以数字增益补偿方法减少时间基准误差，可提高飞行时间的测量精度。

实际上，超声信号波形易受环境干扰影响，测量波形的传输时间，其难度在于信号起点的确定。文献[16]提出了一种基于模型的超声波渡越时间测量的新方法，通过超声波接收探头获取声波信号，然后建立超声波接收信号从起振到稳定的数学模型，最后通过模型参数拟合得到渡越时间参数最优值。当信号起振点受污染时，信号的傅里叶变换相位信息并不会受到影响，文献[17]则利用此信息提出频率估计算法，获得较好的精度与稳定性。

针对传统阈值确定比较困难，文献[18]提出一种统计阈值方法，运用小波分析方法完成噪声映射，并通过概率建模获得最佳阈值等方法测得渡越时间。文献[3]则提出双阈值比较法，对接收信号进行细致分析，获得较为正确的超声到达时间。

文献[19]提出单周期相关滤波法测量超声波渡越时间，结构简单，减少了传统相关算法的计算量，但相关峰值点与信号采样频率相关。

本文根据超声传播原理，针对飞行时间测量中接收起振点难以确定的难题，提出了一种时频测量综合算法，该方法将时域互相关法和频域相位测量方法相结合，提高飞行时间的测量精度，从而提高超声声速的测量精度。

2 基本原理

超声波测量的基本原理如图1所示，发射单频正弦信号激励超声波换能器的发射端发射声波，经过一定时间传播，换能器接收端接收到超声信号并将其转换为电信号。根据测量原理，超声波的速度如式(1)所示：

图1 声速测量模型

v=L/T

(1)

式中：L为超声波测量装置中两个换能器之间的距离；T为超声波的飞行时间。

对飞行时间直接测量，提高测量精度存在以下几个问题较难解决：一是特定频率下进行采样时不能获得信号到达时间的同步信息，即信号到达时刻与采样时刻不同步；二是超声波接收器在接收信号时刻存在起振阶段，导致接收信号的开始部分不稳定，其波形如图2所示，时域内受噪声等因素干扰较大，起振点Q时刻无法精确检测。

图2 接收信号示意图

在不考虑高精度测量的背景下，往往采用阈值法计算超声波的飞行时间，即对超声波接收器的起振时刻进行了时域滤波。从图2可知，接收信号的起点不仅与接收换能器的起振有关，而且受噪声等因素影响较大，从而制约了此类时域方法的测量精度。

为提高飞行时间测量精度，解决超声波接收器起振时刻测不准问题，本文提出了一种采用时域互相关算法与频域相位测量相结合的综合算法对飞行时间进行精确测量，即将飞行时间分成2部分，如式(2)所示：

T=T0+T1

(2)

式中：T0是指通过粗尺度的方法测量，即采用时域互相关的方法测量飞行时间中信号的整周期部分时间，在这个阶段，不用考虑测量的精度；T1则指通过细尺度的方法测量，即采用FFT的方法在频域比较发送信号与接收信号的相位，两者间的相位差即飞行时间的相位差，其具体示意图如图3所示。

图3 时间测量模型

3 时频测量综合算法

超声波速度等参数的测量，主要是先求超声波的飞行时间。在本文中，先通过互相关法粗略地估算出超声波飞行时间的整周期数；然后通过FFT处理，在频域计算出激励信号与输出信号之间的相位差，避免了由于换能器起振点的因素对测量结果造成的影响；最后通过两信号之间的相位差校准整周期数，从而得到高精度的测量结果，其算法如图4所示。

图4 时频测量综合算法

3.1 基于互相关运算的飞行时间正周期测量

为了抑制系统噪声和环境噪声，通常采用时域互相关算法来检测特定的高性能信号的到达时刻。超声波测量中，接收信号参数，包括频率，与原始发送信号的参数基本相同。运用互相关运算检测相关峰时刻，获得超声到达时刻，因存在信号起振点易受测量系统各种因素影响，该测量值存在较大的不确定性。该阶段仅需提取飞行时间的整周期部分。

假设激励信号为频率f0的正弦信号，则x1(t)、x2(t)分别为超声波测量系统中的发送和接收连续信号。通过双路AD采样模块以一定的同步采样率fs对其进行采样，形成离散序列x1(n)和x2(n)，n=0,1,2,…,N-1。它们的互相关函数如式(3)所示：

r=0,1,2，…，N-l

(3)

式中：l为进行互相关运算的信号长度，取4到8个超声激励周期即可。互相关函数R(r)的峰值表示为R(q)，如式(4)所示：

R(q)=rmax, s.t.

R(rmax)=max(R(r)),r=0,1,…,N-l

(4)

式中：q点位于振动起点的相邻区域内，max()函数是选择互相关函数R(r)的最大值。基于以上原因，q点不是直接用来计算超声波的飞行时间，而是用来计算传播延时的整数波形，波数M的计算公式如式(5)所示：

M=rounding[qf0/fs]

(5)

式中：rounding()函数是向下取整函数，这样传输延时的整周期数延时如式(6)所示：

T0=M*1/f0

(6)

式(6)中可以获得超声波延时的整数周期延时，其受系统波动、噪声与外部噪声的影响较小。

3.2 基于频域分析的相位差测量

根据超声波的传播原理，正弦信号在传播过程中频率不发生偏移，同时信号的相位信息随着传播时间增加而不断变化。在超声波测量过程中，存在两种振动起振过程。一个由发射信号激励形成振动起振，另一个由接收换能器刚接收到信号而振动的起振现象。因此，超声波测量时对单频信号测量仅考虑相位时，可剔除这些振动起振阶段，从而提高信号的相位测量精确性。

考虑到实际情况，将前2～3个周期的发射信号作为振动起始部分，经1～2周期后，其后的M周期的接收信号作为接收部分。将采样序列x1(n)、x2(n)转化为序列z1(r)和z2(r)，其表达式如式(7)、式(8)所示：

r=0,1,…,N

(7)

r=0,1,…,N

(8)

序列x1(n)和x2(n)之间的相位差与新序列z1(r)和z2(r)之间的相位差相同。此外经过序列变换后，信号中包含的系统干扰更少了。z1(r)和z2(r)的频谱可以通过FFT计算，分别标记为z1(k)和z2(k)。由于发送信号为单频信号，频率为f0，因此0到n-1范围内的对应点m计算如式(9)所示：

(9)

式中：n是FFT变换的数据长度；f0是发射频率；fs是AD转换器的采样频率。

数据序列z1(r)和z2(r)经过FFT变换后，可以通过求峰值获得频率点m，然后可以通过计算三角函数的方法得到其相位信息，其计算如式(10)、式(11)所示。式中Re[Zi(m)]和Im[Zi(m)]分别是Zi(m)的实部和虚部。

(10)

(11)

利用正切三角运算可以求得信号z1(r)和z2(r)的相位差正切值，如式(12)所示：

(12)

通过对式(12)求反正切三角函数运算求得两个信号的相位差Δφ，如式(13)所示。在求得过程中，将式(10)和式(11)代入式(12)，然后求反切就可。

(13)

通过式(13)求得两信号的相位差Δφ的范围为[-π/2，π/2]，实际上两信号的相位差Δφ为[-2π，0]，为了求得两信号真实的相位差，还需通过判断两个信号在频域的虚实部信息得到其相位差的真实信息，如表1所示。

表1 相位修正

3.3 飞行时间的高精度测量修正

当测量得到的飞行时间接近或等于若干个信号周期时，由于系统噪声如量化噪声等影响，极有可能使相位差落入[-2π，0]的边缘。当相位差Δφ位于边缘点0时，将在前一个波形周期引入一个负误差。当相位差Δφ位于一个边缘点-2π时，将在后一个周期引入一个正误差。

为了测量正确的飞行时间，需利用修正后的相位差Δφ来判断整周期初测量M值是否存在误差，也就是根据修正后的相位差Δφ反馈到时域的整周期测量。具体修正方法如下：

(1)判断起振点q是否在整周期倍点的附近[q%20≤2‖q%20≥18]，其中符号%是求余运算，如果不在，M值不变，否则进入下一步；

(2)若q满足q%20≤2，检测的起振点q在整周期倍点的右附近，这时利用相位差来对初测的整周期值M进行验证及修正如式(14)所示：

(14)

(3)若q满足q%20≥18，检测的起振点q在整周期倍点的左附近，这时利用相位差来对初测的整周期值M进行验证及修正如式(15)所示：

(15)

因此最终声速值v的计算如式(16)所示：

(16)

利用频域反馈到时域的分析方法，确保了整周期时间测量的精度，很好地解决了由接收信号点不准而造成的测量误差。

4 实验结果

实验测量平台如图5所示，采用的是XILINX系列XC3S500E-4VQG100C，具有高性能、低功耗的特点，能够灵活进行软件配置和数据存储器配置。AD转换电路选择的是10位双通道模数转换器AD9218，最高采样率可达到105 Msps，能够与TTL和LVCMOS兼容。主控处理器采用的是C8051F340，主要用来计算处理数据，同时将处理结果通过液晶显示，将测量数据通过串口传送到PC机上。

图5 系统测量平台

实验中，测量平台将收发换能器置于水中，发射换能器采用单频信号1 MHz连续正弦波激励。实测中换能器之间的距离L将影响接收信号的数据：当距离L很小时，回波信号会对接收信号产生回波干扰；而距离L过大时需要消耗更大的数据存储单元，不利于数据采集保存与处理。通过实时测量验证，换能器之间的距离在7 mm到25 mm之间可获得很好的接收信号进行分析。

确定本测量平台，换能器之间距离L=12 mm，采用8 k的存储单元，图6为连续测量100组数据的飞行时间测量结果，图7为100组数据对应的飞行时间测量绝对误差，对100组数据进行可信度分析可以得到当可信度为99.73%时，飞行时间的测量误差在±0.5 ns。实验表明时频测量综合算法能够提高飞行时间的测量精度。

图6 100组超声波延时测量

图7 100组超声波延时绝对误差

飞行时间测量结果：

t=(8.344 49±0.15)ns, 置信概率Pa=68.27%;

t=(8.344 49±0.47)ns, 置信概率Pa=99.73%。

系统每3 s进行一次声速测量，并通过串口上传数据。图8为连续测量100组声速值的数据，图9为声速测量的绝对误差，对声速测量值进行可信度分析结果可以看出当可信度为99.73%时，声速的测量误差为±0.008 m/s。

图8 100组超声波声速测量

图9 100组声速测量的绝对误差

声速测量结果：

v=(1 452.279±0.002 5)m/s, 置信概率Pa=68.27%；

v=(1 452.279±0.007 5)m/s, 置信概率Pa=99.73%。

5 结 论

本文提出了一种基于时频的综合测量算法以获得高精度超声波声速。首先，通过时域互相关法求取飞行时间的整周期值，而后，通过相位差的高精度测量获得声波传输的群延时的测量。最后，由相位差的实时测值来对整周期延时时间进行修正，确保声波传输延时时间的准确性。因此，互相关与相位测量相结合的综合算法可提高飞行时间的测量精度，从而实现高精度的超声波声速测量。