流量传感器(2)—超声流量传感器，相位差和相关性原理

之前我们简单介绍了基于差压的流量传感器，这次我们要开始基于超声的流量计。超声流量计是一种测量液体或气体流量的仪器，其工作原理基于超声波在介质中传播的特性。这种流量计可以在无需接触或干扰流体本身的情况下进行测量，因此被广泛应用于各种工业环境，尤其是在需要无接触或无损测量的场合。

超声流量计的应用领域非常广泛，从水处理和污水处理设施，到化工和石油行业，再到食品和饮料生产线。这些设备还在航空、电力和热力行业中发挥着关键作用，帮助这些行业的专业人士监控系统的效率，评估设备的性能，以及确保工艺流程的连续性。

而今，随着科技的不断发展，超声流量计的准确性和可靠性得到了大幅提升，让我们能够在不同的场景和环境中，对流体流动有着更精确的控制和理解。

2. 相关性原理

• 超声波发射：超声波流量计中的超声换能器发出的超声波通常会以一定的角度，被发射到流体中。精准记录发射、接收时间是进行相关性分析的一个基本要素。

• 接收信号：流量计中的两个超声接收器会分别接收顺流和逆流的超声波信号。由于流速的原因，这两个信号从各自的发射到接收之间，会有一定的时间差。如果不考虑多普勒频移，接收到的两组信号可以认为是同频及波形相似的。

• 对收到的信号进行预处理：包括但不限于降噪、滤波等。

• 计算相关性：将接收到的两个信号进行相关性计算，得到相关函数。相关函数的峰值表示两个信号最相关（也就是最相似）的位置以及所表示的相位差。

• 计算相位差：根据计算相关性序列中得到的对应时序，得到所代表的两信号间在波形上的相位差。而经过给定的参数条件，即可获取两组信号在接收时的时间差。

• 计算流速：通过已知的声速，路径长度，以及由相位差得到的时延，本质上也是通过计算TOF的差异，再得到流速。一般可以假定流体为无损耗的流体，流体的湍流程度对测量结果影响可以根据实际情况调整算法。一定要注意：温度对于声速在流体中的影响是明显的，实际应用中一定要处理好温度对于流量计算的补偿。

• 计算流量：最后，利用流速和管道断面积的信息，可以得到流体的体积流量，即通过管道的流体体积与时间的比值：Q=Av （A为管道横截面积，v为流速），本文在计算得到流速后，统一不再计算流量。

• c = 4500     # 流体中的声速 m/s

• d = 0.5       # 管径(m)

• θ = math.pi / 6       # 超声换能器的安装角,30 deg.

• L = d / math.cos(theta) # 理想情况下的声音传播路径长(m)

• v = 4           # 流体流速(m/s)-用于生成模拟信号，也是待比较值

• f = 1e6        # 超声波的驱动频率 1MHz

• T = 1 / f      # 驱动频率的周期时间(s)

• Fs = 125e6  # 采样频率为125MHz

• Ts = 1/Fs     # 采样间隔

图-4    无噪声正弦信号的相关性检测

• Length of t lags: 2497
• Length of t serial: 1249
• Ts Delayed cycles: 14
• Delayed time:  1.1208967173738991e-07s 
• shift degree: 40.352281825460366
• Sample time interval: 8.00640512409928e-09s
• Estimated velocity:  3.9314356314212877(m/s)

图-5    幅值不同的两个正弦信号的相关性检测

图-6    带加性噪声的两个正弦信号的相关性检测

• Length of t lags: 2497

• Length of t serial: 1249

• Ts Delayed cycles: 17

• Delayed time:  1.3610888710968775e-07 

• degree: 48.9991993594876

• Sample time interval: 8.00640512409928e-09

• Estimated velocity:  4.773885690505651(m/s)

========================

import numpy as np

import matplotlib.pyplot as plt

from scipy.signal import correlate

import math

def cor_demo():

    try:

        # Parameters

        c = 4500    # Sound speed in the fluid in m/s

        d = 0.5     # Pipe diameter in meters

        theta = math.pi / 6     # Angle of the emitted ultrasound wave, in radians. 30 degrees here.

        L = d / math.cos(theta) # The sound path length in meters parallel to the flow direction

        v = 4                   # Fluid velocity in m/s

        f = 1e6     # Frequency of the signal in 1MHz

        T = 1 / f   # Period of the signal in s

        Fs = 125e6 # 采样频率为125MHz

        Ts = 1/Fs  # 采样间隔

        # Calculate time delay caused by flow velocity

        time_up = L / (c - v * math.sin(theta)) # Time of flight upstream

        time_down = L / (c + v * math.sin(theta)) # Time of flight downstream

        #t_delay = L / c * (1 / np.sqrt(1 - (v / c * np.sin(theta)) ** 2) - 1)

        print("T1_0:",time_up-time_down)

        t_delay = 2 * L * v * np.sin(theta) / (c**2 - v**2 * np.sin(theta)**2)

        # Time array

        # 返回一个从0到period*T范围内的有num_samples个元素的一维数组，数组中的数值是等间距分布的

        period = 10

        num_samples = period*T/Ts

        t = np.linspace(0, period*T, int(num_samples))

        # Signals

        s0 = np.sin(2 * np.pi * f * t)

        s1 = np.sin(2 * np.pi * f * (t - t_delay))  #*2.0 # Signal 1

        #=====mask this seg to remove noise from signals=========

        noise0 = np.random.normal(0, 0.5, s0.shape)

        noise1 = np.random.normal(0, 0.5, s1.shape)

        s0 = s0 + noise0

        s1 = s1 + noise1

        #=====mask the seg above to remove noise from signals=========

        # 计算Cross-correlation，并找到最大值对应的位置

        correlation = correlate(s1, s0, method='direct', mode='full')

        lags = np.arange(-len(s1) + 1, len(s1)) # Lags array

        print("Length of t lags:", len(lags))

        # Calculate flow speed using the estimated time delay

        # Find the peak of the cross-correlation corresponds to the time delay

        print("Length of t serial:", len(t))

        delay = lags[np.argmax(correlation)]            # 相位差所对应的信号序列值

        sample_time = (period*T) / len(t)               # 采样间隔时间

        time_delay = delay * sample_time                # 两个信号间的延迟时间

        phase_shift = (time_delay / T) * 2 * np.pi      # 由时间延迟换算成的两个信号序列的相位差

        phase_shift_deg = phase_shift * (180 / np.pi)   # 由相位差换成的两个信号序列的角度差

        print("Ts Delayed cycles:", delay, ", Delayed time: ", time_delay, "degree:", phase_shift_deg) # 将延迟转换为相位差

        print("Sample time interval:",sample_time)

        # 计算所得的流速

        v_estimated = (math.sqrt((L**2)+(time_delay**2)*(c**2))-L)/(time_delay * math.sin(theta))

        print('Estimated velocity: ', v_estimated)

        # 画图

        #fig, (ax_origin, ax_corr, ax_shift) = plt.subplots(2, 1, figsize=(12, 8))

        fig, (ax_origin, ax_corr) = plt.subplots(2, 1, figsize=(12, 8))

        # 原始信号图

        ax_origin.plot(t, s0, label='Signal 1')

        ax_origin.plot(t, s1, label='Signal 2')

        ax_origin.set_title('Original signals')

        ax_origin.legend()

        # 互相关图

        ax_corr.plot(correlation)

        ax_corr.axvline(x = len(correlation)//2 + delay, color = 'r', linestyle = '--', label = "Max correlation at delay")

        ax_corr.set_title('Cross-correlation between signal 1 and signal 2')

        ax_corr.legend()

        plt.tight_layout()

        plt.show()

        return

    except Exception as e:

        print("Error:",e)

if __name__=='__main__':

    cor_demo()  

========================

需要了解的是，在信号采集过程中，足够快的采样率是保证相位差有足够分辨率的保证。对应的是模拟代码中“Ts Delayed cycles”。容易设想，如果我们可以在处理器足够快的前提下，还可以得到更多的“Ts Delayed cycles”，就意味着相位差的分辨率就可以足够高。

在上面的模拟过程中，DAQ的采样频率是125MHz，对应一些高速ADC器件。反之呢？如果实际使用的ADC器件的采样率因为某些原因速度并不是那么高时，我们的模拟结果会出现什么情况？

下面的波形是将ADC的采样频率调整为12.5MHz时的输出。

图-7    降低信号采集频率后的相关性检测

      Estimated velocity:  2.828550434230693(m/s)——插值前

图-8    插值法提高相关性检测精度的输出

图-9    带噪声未滤波情况下的插值后相关性检测输出

插值运算相位差和流速。下面是模拟代码，仅供参考。以下代码中，噪声部分是注释掉的。

========================

import numpy as np

import matplotlib.pyplot as plt

from scipy.signal import correlate

import math

from scipy.interpolate import interp1d

def cor_demo():

    try:

        # Parameters

        c = 4500    # Sound speed in the fluid in m/s

        d = 0.5     # Pipe diameter in meters

        theta = math.pi / 6     # Angle of the emitted ultrasound wave, in radians. 30 degrees here.

        L = d / math.cos(theta) # The sound path length in meters parallel to the flow direction

        v = 4                   # Fluid velocity in m/s

        f = 1e6     # Frequency of the signal in 1MHz

        T = 1 / f   # Period of the signal in s

        Fs = 125e5 # 采样频率为12.5MHz

        Ts = 1/Fs  # 采样间隔

        # Calculate time delay caused by flow velocity

        time_up = L / (c - v * math.sin(theta)) # Time of flight upstream

        time_down = L / (c + v * math.sin(theta)) # Time of flight downstream

        #t_delay = L / c * (1 / np.sqrt(1 - (v / c * np.sin(theta)) ** 2) - 1)

        print("T1_0:",time_up-time_down)

        t_delay = 2 * L * v * np.sin(theta) / (c**2 - v**2 * np.sin(theta)**2)

        # Time array

        # 返回一个从0到period*T范围内的有num_samples个元素的一维数组，数组中的数值是等间距分布的

        period = 10

        num_samples = period*T/Ts

        t = np.linspace(0, period*T, int(num_samples))

        # Signals

        s0 = np.sin(2 * np.pi * f * t)

        s1 = np.sin(2 * np.pi * f * (t - t_delay))  #*2.0 # Signal 1

        """

        noise0 = np.random.normal(0, 0.5, s0.shape)

        noise1 = np.random.normal(0, 0.5, s1.shape)

        s0 = s0 + noise0

        s1 = s1 + noise1

        """

        # Define interpolation factor

        interp_factor = 10

        # New time vector after interpolation

        t_new = np.linspace(t.min(), t.max(), t.size * interp_factor)

        # Create a function based on the original signals, which can be used to generate the interpolated signals

        interp_func_s0 = interp1d(t, s0, kind='cubic')

        interp_func_s1 = interp1d(t, s1, kind='cubic')

        # Generate the interpolated signals

        s0_new = interp_func_s0(t_new)

        s1_new = interp_func_s1(t_new)

        # 计算Cross-correlation，并找到最大值对应的位置

        #correlation = correlate(s1, s0, method='direct', mode='full')  # old

        correlation = correlate(s1_new, s0_new, method='direct', mode='full')

        lags = np.arange(-len(s1_new) + 1, len(s1_new)) # Lags array

        print("Length of t lags:", len(lags))

        # Calculate flow speed using the estimated time delay

        # Find the peak of the cross-correlation corresponds to the time delay

        #print("Length of t serial:", len(t))

        delay = lags[np.argmax(correlation)]            # 相位差所对应的信号序列值

        #sample_time = (period*T) / len(t)               # 采样间隔时间  old

        sample_time = (period*T) / len(t_new)

        time_delay = delay * sample_time                # 两个信号间的延迟时间

        phase_shift = (time_delay / T) * 2 * np.pi      # 由时间延迟换算成的两个信号序列的相位差

        phase_shift_deg = phase_shift * (180 / np.pi)   # 由相位差换成的两个信号序列的角度差

        print("Ts Delayed cycles:", delay, ", Delayed time: ", time_delay, "degree:", phase_shift_deg) # 将延迟转换为相位差

        print("Sample time interval:",sample_time)

        # 计算所得的流速

        v_estimated = (math.sqrt((L**2)+(time_delay**2)*(c**2))-L)/(time_delay * math.sin(theta))

        print('Estimated velocity: ', v_estimated)

        # 画图

        #fig, (ax_origin, ax_corr, ax_shift) = plt.subplots(2, 1, figsize=(12, 8))

        fig, (ax_origin, ax_interpo, ax_corr) = plt.subplots(3, 1, figsize=(12, 8))

        # 原始信号图

        ax_origin.plot(t, s0, label='Signal 1')

        ax_origin.plot(t, s1, label='Signal 2')

        ax_origin.set_title('Original signals')

        ax_origin.legend()

        # 插值后的信号图

        ax_interpo.plot(t_new, s0_new, label='Signal 1')

        ax_interpo.plot(t_new, s1_new, label='Signal 2')

        ax_interpo.set_title('interpolated signals')

        ax_interpo.legend()

        # 互相关图

        ax_corr.plot(correlation)

        ax_corr.axvline(x = len(correlation)//2 + delay, color = 'r', linestyle = '--', label = "Max correlation at delay")

        ax_corr.set_title('Cross-correlation between signal 1 and signal 2')

        ax_corr.legend()

        plt.tight_layout()

        plt.show()

        return

    except Exception as e:

        print("Error:",e)

if __name__=='__main__':

    cor_demo()  

========================

另外一个大家可以尝试比较的是：根据奈奎斯特采样定理，信号的采集频率至少应该是被采集信号带宽的2倍以上，才可以从采集信号中回复原信号。是不是还可以再降低一些采样频率然后通过插值的方式来提高测量精度呢？

5. 小结

本文中提高的模拟信号及分析方式仅供理论参考，实际的超声接收信号不会是图示中的那样都是完美的正余弦波形。

超声流量计通过利用超声波在流体中传播的特性来测量流量，其中使用相位差方式是一种常见的方法。它的特点主要包括：

• 非接触测量：通过在管道外壁上安装超声探头，可以避免对流体产生任何影响，适合处理易激发、高温或腐蚀性的流体。

• 较高精度：相位差法可以实现较高的测量精度，因为它直接测量的是流体的流动速度，而非其他的，可能受到流体性质变化影响的参数。

• 广泛适用：适用于各种类型的流体，包括气体、液体和蒸汽，并且管道大小范围广泛。

• 流体的温度、压力和流速状态变化：这些因素可能影响超声波在流体中的传播速度，从而影响测量结果，尤其是温度对于超声波在液体中的传播速度的影响。

• 超声探头的安装位置和角度：如果探头的安装位置或角度不准确，会导致超声波不能正确地穿过流体，影响测量结果。

• 管道和流体的不同：如果流体的性质（如粘度、密度）或者管道的材料和状况发生变化，可能需要重新校准设备。
• 流体中的颗粒或气泡：这些可能会吸收、散射或者反射超声波，从而影响信号质量。
• 以上并未考虑由于流速引起的多普勒频移影响。

在超声换能器领域，安费诺传感器（Amphenol Sensors）有Piezo Technologies和Airmar两个自有知名品牌。需要温度和压力补偿时，我们当然也是品牌产品多有备选。

[1] 《数字信号处理—原理、算法与应用》，第四版，普罗克斯著；方艳梅等译， 电子工业出版社，2014.8

[本公众号演示代码版权所有 (c) [2024] Amphenol Sensors]

特此免费授予任何获得此软件及相关文档文件（“软件”）副本的人无限制处理软件的权利，包括但不限于使用、复制、修改、合并、发布、分发、再许可和/或销售软件副本的权利，并允许向其提供软件的人员行使这些权利，但需满足以下条件：

上述版权声明和此权限声明应包含在所有副本或实质部分的软件中。

本软件按“原样”提供，无任何形式的明示或暗示保证，包括但不限于适销性、特定目的适用性和非侵权性的保证。在任何情况下，作者或版权所有者均不承担任何索赔、损害或其他责任，无论是在合同行为、侵权行为或其他行为中，还是与软件或软件的使用或其他交易有关。