1. 引言

之前我們簡單介紹了基于差壓的流量傳感器，這次我們要開始基于超聲的流量計。超聲流量計是一種測量液體或氣體流量的儀器，其工作原理基于超聲波在介質(zhì)中傳播的特性。這種流量計可以在無需接觸或干擾流體本身的情況下進行測量，因此被廣泛應(yīng)用于各種工業(yè)環(huán)境，尤其是在需要無接觸或無損測量的場合。

 

超聲流量計的應(yīng)用領(lǐng)域非常廣泛，從水處理和污水處理設(shè)施，到化工和石油行業(yè)，再到食品和飲料生產(chǎn)線。這些設(shè)備還在航空、電力和熱力行業(yè)中發(fā)揮著關(guān)鍵作用，幫助這些行業(yè)的專業(yè)人士監(jiān)控系統(tǒng)的效率，評估設(shè)備的性能，以及確保工藝流程的連續(xù)性。

 

此外，超聲流量計還在醫(yī)療領(lǐng)域有著重要應(yīng)用。例如，醫(yī)護人員可以利用超聲流量計監(jiān)測血液或藥物的輸送速度，對心血管疾病的診斷和治療中起到關(guān)鍵作用。

 

而今，隨著科技的不斷發(fā)展，超聲流量計的準確性和可靠性得到了大幅提升，讓我們能夠在不同的場景和環(huán)境中，對流體流動有著更精確的控制和理解。

2. 相關(guān)性原理

相關(guān)性是一個統(tǒng)計學(xué)術(shù)語，用來描述兩個或多個隨機變量之間的統(tǒng)計關(guān)系。在信號處理中，原則上這個關(guān)系也可以不必非要是線性的，相關(guān)性分析也可以衡量兩個信號的非線性關(guān)系。相關(guān)性的計算通常通過相關(guān)系數(shù)（Pearson相關(guān)系數(shù)為最常用的一種）來實現(xiàn)。

 

相關(guān)系數(shù)的值范圍在-1到1之間，0表示兩個變量獨立無關(guān)，-1表示兩個變量完全負相關(guān)，1表示兩個變量完全正相關(guān)。

 

在信號處理中，通常會用到互相關(guān)（Cross-Correlation）函數(shù)，它是對兩個信號的相似性進行度量。這在很多應(yīng)用中非常有用，比如在無線通信中的信號識別，圖像處理中的模式匹配，以及地球物理勘探中地震波的分析等。計算互相關(guān)函數(shù)的公式為：

其中，x和y分別為兩個信號，Rxy(t)為在延時 t 時的互相關(guān)值。在離散時間序列中，可以用求和代替積分進行計算。

其中，l>=0時，i=l,k=0，并且l<0時，i=0，k=l。

以上公式中，t和l表示時間的移位；下標xy的順序，表示一個變量或序列相對于留一個變量或序列的移位方向：上面兩個等式中，變量x或者序列x(n)未移動，而y或者y(n)在時間上移動了t或l個單位：向右移動則t或l為正；向左移動則t或l為負。數(shù)據(jù)相關(guān)性分析對于有加性噪聲的信號分析有其獨到的特點。舉現(xiàn)成的栗子[1]： 

這里x(n)（紅色部分正弦線條）是一個周期性的序列，周期N未知，ω(n)表示隨機噪聲，藍灰色信號線為y(n)=x(n)+ω(n)。我們怎么去了解該帶噪聲信號y(n)的周期？[1](這里省掉了對于周期信號相關(guān)性計算特點的陳述) 

盡管小編也很想花些篇幅來說明相關(guān)性應(yīng)用的特點，限于篇幅，我們對于上面的周期為N的信號自相關(guān)的結(jié)果中，Rxx在l=0，N和2N等處有相對較大的峰值；信號x(n)和加性噪聲ω(n)的互相關(guān)Rxω和Rωx相對很??；而隨機噪聲信號的自相關(guān)Rωω序列，除了在l=0處有峰值外（每個對應(yīng)點自相乘后累加），其余點將很快衰減到0。

這里忘了說明：周期信號的自相關(guān)序列也是一個同周期的序列。也就是如下公式所示的周期為N的序列x(n)自相關(guān)[1]：

在這種情況下，我們計算y(n)的自相關(guān)序列后，可以得到噪聲淹沒下的周期，如下的圖示的每2個紅點之間的時間所代表的周期結(jié)果：噪聲下的周期信號的相關(guān)性周期檢查演示代碼如下所示：

import numpy as np

 

import matplotlib.pyplot as plt

 

from scipy.signal import find_peaks

 

import numpy as np

 

import matplotlib.pyplot as plt

 

from scipy.signal import find_peaks

 

 

# 創(chuàng)建一個正弦信號

 

T = 1/50 # 50Hz信號，周期為1/50，這里實際設(shè)置為0.02秒

 

fs = 5000 # 抽樣頻率

 

t = np.arange(0, 1, 1/fs) # 持續(xù)時間為1秒

 

signal = np.sin(2 * np.pi * 50 * t)

  

 

# 添加噪聲

 

noise = np.random.normal(0, 2, signal.shape)

 

noisy_signal = signal + noise

 

 

# 計算信噪比SNR，以dB表示

 

SNR = 10 * np.log10(np.mean(signal**2) / np.mean(noise**2))

 

print('SNR: {:.2f} dB'.format(SNR))

 

 

# 計算并顯示帶噪聲的正弦波和清潔的正弦波

 

plt.figure(figsize=(12, 8))

 

plt.subplot(2, 1, 1)

 

plt.plot(t, noisy_signal, label='Noisy signal')

 

plt.plot(t, signal, 'r', label='Original signal')

 

plt.legend(loc='upper right')

 

plt.xlabel('Time [s]')

 

plt.ylabel('Amplitude')

 

plt.title('Noisy and Original Signal')

 

 

# 計算自相關(guān)函數(shù)

 

auto_corr = np.correlate(noisy_signal, noisy_signal, mode='full')

 

# 找到自相關(guān)函數(shù)的峰值檢查閾值[max()/12]來估計周期

 

auto_corr = auto_corr[auto_corr.size // 2:]     #取相關(guān)信號結(jié)果的后半截數(shù)據(jù)

 

peaks, _ = find_peaks(auto_corr, height=auto_corr.max()/12, distance=fs/60) # 最小距離設(shè)置為一個頻率60Hz的周期

 

estimated_period = np.diff(peaks) / fs

 

print('Estimated period: {:.2f} s'.format(estimated_period.mean()))

 

 

# 繪制自相關(guān)函數(shù)并顯示峰值

 

plt.subplot(2, 1, 2)

 

plt.plot(t, auto_corr / fs, 'g') # 除以fs以便時間以秒為單位

 

plt.plot(peaks / fs, auto_corr[peaks] / fs, 'ro')

 

plt.xlabel('Time delay [s]')

 

plt.ylabel('Autocorrelation function')

 

plt.title('Autocorrelation of Noisy Signal')

 

plt.tight_layout()

 

plt.show()

 

運行測試可以得到正確的結(jié)果：Estimated period: 0.02 s

圖-2    周期信號通過相關(guān)性序列的周期檢測

 

計算互相關(guān)函數(shù)需要注意以下幾點：

相關(guān)并不表示因果關(guān)系。即使兩變量存在高度相關(guān)，也只能表明它們之間存在某種關(guān)系，但并不能確定是哪個變量影響了另一個變量。

為了準確計算互相關(guān)函數(shù)，關(guān)鍵在于靈敏地選擇合適的時延值。合適的時延值取決于所考慮系統(tǒng)的特性及其應(yīng)用。

對于非平穩(wěn)信號（信號的統(tǒng)計特性會隨著時間變化），可能需要使用互相關(guān)函數(shù)的改進版，如互協(xié)方差函數(shù)等。

需要注意，互相關(guān)函數(shù)對噪聲非常敏感，信號的噪聲會直接影響到相關(guān)運算的結(jié)果，因此在進行相關(guān)運算之前，通常需要對信號進行去噪處理。

計算復(fù)雜度高，尤其是數(shù)據(jù)量大的時候。對于大數(shù)據(jù)集，可采用快速傅里葉變換（FFT）進行降維計算。FFT算法的時間復(fù)雜度為O(NlogN)，而直接在時域計算互相關(guān)函數(shù)的時間復(fù)雜度則是O(N^2)。

 

3. 超聲波流量計中的相關(guān)性原理的應(yīng)用

圖-3    超聲流量計的配置圖

 

要開啟對超聲流量計工作原理的討論，關(guān)鍵在于理解超聲波在流體中傳播的特性以及兩者之間如何互動。

 

超聲流量計采用了所謂的“TOF-時間飛逝”方法進行測量。超聲收發(fā)的配置會存在差異。假設(shè)在這種設(shè)備中，有兩對超聲換能器A和B，按照上圖所示被置于管道的兩側(cè)，都是向下發(fā)射，下面的接收。當沒有流體流過的時候，從一個傳感器發(fā)出的超聲波到達另一個傳感器所需的時間應(yīng)該是恒定的。然后，當流體開始流過時，這個時間將會發(fā)生變化。

 

超聲波流量計的構(gòu)造有很多種類，簡單以上面的圖-3所示。具體來說，當流體如圖示的方向流動時，從A發(fā)出的超聲波會被流體帶著一起移動，所以它到達接收端的時間t2會少于沒有流體流動的情況。反之，從B發(fā)出的超聲波在到達接收端的過程中需要逆流而上，所以耗費的時間t1會更多。

 

在流量計內(nèi)部，這兩個時間被用來計算出一種稱為“時間差”的值delay，這就是流體流速的直接度量參數(shù)。根據(jù)這個delay所對應(yīng)的時間來計算流速的公式見3.2段。

 

然后，通過對管道的直徑和流速的測量，便可以計算出流量。有關(guān)超聲波在流體中的傳播和接收時間差的形成的推導(dǎo)部分，請參考我們之前我們的博文：《如何消除超聲波在流量測量及液位測量中溫度的影響》，其中對超聲檢測流量稍微有一些說明，不過是基于單純TOF時差的方式，在實際應(yīng)用中，超聲信號可能會受到各種類型噪聲的干擾。如果信號受到過多的干擾，將無法正確測定時差。為了獲取準確的流速讀數(shù)，我們需要盡量消除這些噪聲帶來的影響。

 

既然我們看到了相關(guān)性對于噪聲影響的抵消能力，我們就此引出利用信號的相關(guān)性進行流速測量的思路方法。

 

本文不是重點于如何設(shè)計，而是提供一種解決問題的思路。

 

3.1 使用相關(guān)性檢測相位差法的超聲流量計的處理大致步驟

超聲波發(fā)射：超聲波流量計中的超聲換能器發(fā)出的超聲波通常會以一定的角度，被發(fā)射到流體中。精準記錄發(fā)射、接收時間是進行相關(guān)性分析的一個基本要素。

接收信號：流量計中的兩個超聲接收器會分別接收順流和逆流的超聲波信號。由于流速的原因，這兩個信號從各自的發(fā)射到接收之間，會有一定的時間差。如果不考慮多普勒頻移，接收到的兩組信號可以認為是同頻及波形相似的。

對收到的信號進行預(yù)處理：包括但不限于降噪、濾波等。

計算相關(guān)性：將接收到的兩個信號進行相關(guān)性計算，得到相關(guān)函數(shù)。相關(guān)函數(shù)的峰值表示兩個信號最相關(guān)（也就是最相似）的位置以及所表示的相位差。

計算相位差：根據(jù)計算相關(guān)性序列中得到的對應(yīng)時序，得到所代表的兩信號間在波形上的相位差。而經(jīng)過給定的參數(shù)條件，即可獲取兩組信號在接收時的時間差。

計算流速：通過已知的聲速，路徑長度，以及由相位差得到的時延，本質(zhì)上也是通過計算TOF的差異，再得到流速。一般可以假定流體為無損耗的流體，流體的湍流程度對測量結(jié)果影響可以根據(jù)實際情況調(diào)整算法。一定要注意：溫度對于聲速在流體中的影響是明顯的，實際應(yīng)用中一定要處理好溫度對于流量計算的補償。

計算流量：最后，利用流速和管道斷面積的信息，可以得到流體的體積流量，即通過管道的流體體積與時間的比值：Q=Av （A為管道橫截面積，v為流速），本文在計算得到流速后，統(tǒng)一不再計算流量。

 

以上就是使用相關(guān)性檢測相位差法的超聲流量計的基本處理步驟。需要注意的是，這只是大致步驟，每一步可能會根據(jù)具體實現(xiàn)，使用的硬件設(shè)備以及信號處理方法的不同而有所差異。

 

3.2 使用相關(guān)性檢測相位差法的超聲流量計的仿真

依舊參考前面的構(gòu)造圖，并且設(shè)定一些已知條件（比如管徑，流體的聲速—溫度暫不考慮，流速，超聲換能器的安裝角，超聲波的驅(qū)動頻率，信號的采樣頻率等）來看模擬的結(jié)果和給定的已知參數(shù)是否相近或者相等。這里我們以給定的流速和模擬測量結(jié)果中的流速為比較。

 

在以下的模擬過程中，兩路信號各自的生成采集起始時間是需要盡可能準確記錄，這樣，我們以其中的一路為信號為相關(guān)性處理的相對基準序列時，在存在流速的情況下，另外一路和該路信號中間就存在相位差。

在計算所得的相關(guān)性序列數(shù)據(jù)中，找到相關(guān)數(shù)值Rxy最大值所對應(yīng)的時間“l(fā)”——在波形上是相位——就找到了兩組信號的相位差，亦即時間差。以上接收信號序列中的每個標識，意味著一次數(shù)據(jù)采集的結(jié)果——后文會提到：數(shù)據(jù)的采樣頻率，也會影響計算的精度。在這里，兩路接收信號之間的相位差，就是我們要相關(guān)性處理測量計算的中間變量。整個模擬過程中，我們按照之前的分析，暫時忽略由于流速導(dǎo)致的多普勒頻偏的影響。根據(jù)信號的相位差得到時差后，換算成流速的公式如下：

無噪聲情況下的模擬
由以下條件生成的模擬信號，然后看看通過信號的處理是否可以還原流體的流速。

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        # 超聲波的驅(qū)動頻率 1MHz

T = 1 / f      # 驅(qū)動頻率的周期時間(s)

Fs = 125e6  # 采樣頻率為125MHz

Ts = 1/Fs     # 采樣間隔

圖-4    無噪聲正弦信號的相關(guān)性檢測

 

仿真的結(jié)果如下：

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)

可以看到，通過模擬信號進行相關(guān)性處理之后得到的流體速度計算結(jié)果，和模擬信號中設(shè)定的“實際”流速很接近，僅1.71%的誤差。
如果接收信號的幅值不一樣會有什么影響？我們來模擬一下：兩個超聲換能器接收到的信號乃至放大之后采集到的信號仍然存在不一致時，對測量結(jié)果是否會有影響？下圖中信號#2的幅值已經(jīng)是信號#1幅值的2倍了。

圖-5    幅值不同的兩個正弦信號的相關(guān)性檢測

 

我們看到，接收信號的幅值不一樣的情況下，測量所得的流速和前面信號同幅值時是一樣的。可見理想情況下，幅值對測量值影響不是首位的，但是實際應(yīng)用中肯定存在噪聲，而有噪聲我們就需要強調(diào)信噪比，因此雖然模擬結(jié)果尚可，但是信號的有效幅值我們還是希望更高些。

Estimated velocity:  3.9314356314212877(m/s)

 

有噪聲情況下的模擬

信噪比是隨便選的，但是我們可以看到在有噪聲的影響下，單純使用時差檢測的方法將無法再檢測時，通過信號序列的相關(guān)性處理，仍然可以有相應(yīng)的輸出。我們將兩路信號加上各自幅值的0.5倍的隨機噪聲。

圖-6    帶加性噪聲的兩個正弦信號的相關(guān)性檢測

 

雖然可以檢測，但是此時得到的結(jié)果是4.77m/s，和理想模擬值流速=4m/s有了較大差異。

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)

可見，噪聲將影響測量的精度。如果我們在設(shè)計硬件過程中，已經(jīng)充分考慮了信號的收發(fā)部分的濾波處理，也將非常有利于最后測量值的精度和穩(wěn)定性。關(guān)于信號濾波這部分，在稍后也需要用到。
以下是上面模擬有、無噪聲信號處理所用的python代碼，僅供參考。代碼中在生成的信號中添加了隨機噪聲，不需要時請注釋掉。

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范圍內(nèi)的有num_samples個元素的一維數(shù)組，數(shù)組中的數(shù)值是等間距分布的

 

        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，并找到最大值對應(yīng)的位置

 

        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)]            # 相位差所對應(yīng)的信號序列值

 

       

        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) # 將延遲轉(zhuǎn)換為相位差

 

        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()

 

 

        # 互相關(guān)圖

 

        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()

 

4. 超聲波流量計中的信號相關(guān)性處理中的插值處理

 

需要了解的是，在信號采集過程中，足夠快的采樣率是保證相位差有足夠分辨率的保證。對應(yīng)的是模擬代碼中“Ts Delayed cycles”。容易設(shè)想，如果我們可以在處理器足夠快的前提下，還可以得到更多的“Ts Delayed cycles”，就意味著相位差的分辨率就可以足夠高。

在上面的模擬過程中，DAQ的采樣頻率是125MHz，對應(yīng)一些高速ADC器件。反之呢？如果實際使用的ADC器件的采樣率因為某些原因速度并不是那么高時，我們的模擬結(jié)果會出現(xiàn)什么情況？

下面的波形是將ADC的采樣頻率調(diào)整為12.5MHz時的輸出。
 

圖-7    降低信號采集頻率后的相關(guān)性檢測
 

采集到的波形來看，正弦波已經(jīng)有點異常了。從模擬結(jié)果來看，即使是其他都是理想的情況下，流速結(jié)果和4m/s的設(shè)定相差就很多了。
Estimated velocity:  2.828550434230693(m/s)——插值前

在這種情況下，有無其他的方式提高最后輸出結(jié)果的精度？答案是肯定的，那就是插值——我們已經(jīng)設(shè)定這個波形是近正余弦的，在python中，我們看到軟件包scipy有3階樣條插值的實現(xiàn)方式函數(shù)：scipy.interpolate.interp1d。以下的仿真結(jié)果就是基于該插值方式實現(xiàn)的。
插值的過程，就是在已有的時間序列中，每個間隔之間均勻地插入設(shè)定的時間點數(shù)，然后通過插值函數(shù)，在已采集的信號（這里是模擬的采集信號）中，對應(yīng)插入的時間點上，使用3階樣條插值，添加額外的插值。
 

圖-8    插值法提高相關(guān)性檢測精度的輸出

上圖是通過插值方法得到的輸出圖。第一個圖，是沒有插值時的2路信號輸出圖，該圖和前一個圖中的信號波形是一致的，都是基于12.5MHz的采樣率生成的。第二個圖和第三個圖，就是在12.5MHz采樣率的采集信號的前提下插值得到的。這里每兩個源數(shù)據(jù)之間插入點數(shù)有10個，實際應(yīng)用過程中，應(yīng)該需要充分考慮處理器的處理速度，流量測量的實時性要求。
Estimated velocity:  3.959970232994945m/s——插值后得到的模擬流速測量結(jié)果也是很接近設(shè)定的模擬流速（4m/s)的。并且，其中的數(shù)據(jù)點，有原先的247個，經(jīng)過插值后，達到2479個，相對于增加了分辨率。需要特別提出，信號的完整性在這里是重要的，濾波！否則插值的結(jié)果也將是亂七八糟的。如下圖，直接在噪聲疊加的信號上進行插值。模擬結(jié)果也測試了，不盡如人意，和設(shè)定值相差較大。大家可以用本文附帶的代碼進行驗證測試一下。

圖-9    帶噪聲未濾波情況下的插值后相關(guān)性檢測輸出
 

插值運算相位差和流速。下面是模擬代碼，僅供參考。以下代碼中，噪聲部分是注釋掉的。

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范圍內(nèi)的有num_samples個元素的一維數(shù)組，數(shù)組中的數(shù)值是等間距分布的

 

        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，并找到最大值對應(yīng)的位置

 

        #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)]            # 相位差所對應(yīng)的信號序列值

    

 

        #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) # 將延遲轉(zhuǎn)換為相位差

 

        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()

 

        # 互相關(guān)圖

 

        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()

 

另外一個大家可以嘗試比較的是：根據(jù)奈奎斯特采樣定理，信號的采集頻率至少應(yīng)該是被采集信號帶寬的2倍以上，才可以從采集信號中回復(fù)原信號。是不是還可以再降低一些采樣頻率然后通過插值的方式來提高測量精度呢？

 

5. 小結(jié)

 

本文中提高的模擬信號及分析方式僅供理論參考，實際的超聲接收信號不會是圖示中的那樣都是完美的正余弦波形。

超聲流量計通過利用超聲波在流體中傳播的特性來測量流量，其中使用相位差方式是一種常見的方法。它的特點主要包括：

非接觸測量：通過在管道外壁上安裝超聲探頭，可以避免對流體產(chǎn)生任何影響，適合處理易激發(fā)、高溫或腐蝕性的流體。

較高精度：相位差法可以實現(xiàn)較高的測量精度，因為它直接測量的是流體的流動速度，而非其他的，可能受到流體性質(zhì)變化影響的參數(shù)。

廣泛適用：適用于各種類型的流體，包括氣體、液體和蒸汽，并且管道大小范圍廣泛。

然而，在使用超聲流量計時，我們需要注意以下可能影響測量精度的問題：

流體的溫度、壓力和流速狀態(tài)變化：這些因素可能影響超聲波在流體中的傳播速度，從而影響測量結(jié)果，尤其是溫度對于超聲波在液體中的傳播速度的影響。

超聲探頭的安裝位置和角度：如果探頭的安裝位置或角度不準確，會導(dǎo)致超聲波不能正確地穿過流體，影響測量結(jié)果。

管道和流體的不同：如果流體的性質(zhì)（如粘度、密度）或者管道的材料和狀況發(fā)生變化，可能需要重新校準設(shè)備。

流體中的顆?；驓馀荩?/strong>這些可能會吸收、散射或者反射超聲波，從而影響信號質(zhì)量。

以上并未考慮由于流速引起的多普勒頻移影響。

以上都是在使用超聲流量計時應(yīng)注意的問題，以確保測量的準確性。
 

6. 安費諾傳感器時間

在超聲換能器領(lǐng)域，安費諾傳感器（Amphenol Sensors）有Piezo Technologies和Airmar兩個自有知名品牌。需要溫度和壓力補償時，我們當然也是品牌產(chǎn)品多有備選。

 

Amphenol的Piezo和Airmar品牌的超聲傳感器以其優(yōu)越的性能和廣泛的應(yīng)用領(lǐng)域，贏得了行業(yè)內(nèi)的高度認可。

 

Piezo品牌的超聲傳感器以其高靈敏度，穩(wěn)定性和耐用性突出，適用于廣泛的密度、粘度和顏色的液壓和液體應(yīng)用等。特別適合需要實現(xiàn)精確控制和測量的場合，例如在醫(yī)療、工業(yè)自動化和汽車等領(lǐng)域的流量控制、液位測量和距離測量等。

 

Airmar品牌的超聲傳感器同樣擁有強大的性能，主要用于航海、船舶、海洋科學(xué)研究和漁業(yè)等領(lǐng)域。它們以其卓越的防水、防腐蝕特性以及在苛刻環(huán)境中的可靠性備受推崇。其超聲傳感器常用于下述場合：深海探測、魚群定位、水深和水溫測量，以及其它海洋環(huán)境參數(shù)的監(jiān)測。

 

無論是需求突出的Piezo，還是適用于海洋環(huán)境的Airmar，Amphenol Sensors都能提供一流的解決方案，滿足各種行業(yè)的需求。