張偉星，李建民，薛鵬飛，侯 文，王 高，劉爭光

(1. 中北大學信息與通信工程學院，山西太原030051；2. 山西省計量科學研究院，山西太原030032；3. 肇慶醫(yī)學高等專科學校，廣東肇慶526020)

0 引 言

與孔板法、渦輪法等傳統(tǒng)流速測量方法相比，超聲波流量計具有非侵入式測量、瞬態(tài)響應快、維護方便、適于大管徑測量等優(yōu)點，已成為管道流速測量的主要方法，廣泛應用于工業(yè)、醫(yī)學和民用等多個領域。隨著超聲波技術的迅速發(fā)展，管道流體的流速測量和監(jiān)測取得了長足發(fā)展。有多種超聲波流量計可用于測量和監(jiān)測液體或氣體流速[1-3]，其中傳輸時差流量計適用于各種尺寸管徑(9 mm～9 m)的流速測量。傳輸時差流量計的技術要點在于超聲波在流體中傳播時，順流和逆流產生的時間差的測量。針對小管徑、低流速的情形，采用高速數據采集電路，對采集的回波信號進行算法處理準確查找到目標波形以及計時終點方法[4]。文獻[4]中的流速測量方法對硬件電路性能要求較高。

電子技術的進步可以比較容易地處理 ns 級的時間間隔，但是確定超聲波信號到達換能器的時刻則比較困難[5-6]。文獻[7]對比了相關函數法與波包幅值法，結果顯示相關函數法精度較高，同時解決了信號信噪比低無法直接讀出波包峰值點的問題。文獻[8]表明互相關時延估計方法能夠在未進行降噪處理的情況下估計兩信號序列的時延。文獻[9]分析計算了管道中層流和湍流的流動狀態(tài)。文獻[10]提到近二十年來，國外的高性能微處理器技術的發(fā)展與成熟，使得超聲波流量計能夠實現高速和復雜的數學計算和邏輯控制，極大地簡化了電路，為超聲波流量計的高精度和智能化的發(fā)展提供了可靠的技術條件和多樣化的平臺以及廣闊的發(fā)展前景。國內對超聲波流量計的研究和應用已經取得較大進展，尤其在中等尺寸管徑的流量測量領域日趨成熟，適用于不同氣體、液體的超聲波流量計相繼問世，有較高的精度和穩(wěn)定性。但是目前，針對細管道中流速的超聲波測量的研究和應用仍不成熟，未能實現大范圍的應用，仍需進一步提高超聲波流量計的測量精度和穩(wěn)定性。文獻[11-13]分析了插值結合其他方法處理超聲波數據改善測量精度的原理。

以上幾種方法在大管徑、高流速的情形下取得了較為理想的測量結果。在實際應用中，對細管道中流速的測量需求越來越迫切，比如醫(yī)療中的精確注射、化工行業(yè)中的精確物料添加等。但是涉及相關技術的文獻較少。本文研究了在5 mm細管道中的流速的測量，利用多物理場有限元仿真軟件COMSOL Multiphysics建立三維管道流流動模型，設置背景流(水)的各種屬性及超聲波在管道流中的傳播狀態(tài)，以此作為實驗的基礎。細管道的特點是管徑小，超聲波傳播路程短，從而導致測量結果的穩(wěn)定性和精度較差，所以難點在于如何提高測量的穩(wěn)定性和精度。本文的解決方法為：(1) 采用相關函數法計算時間間隔，一定程度上減少噪聲的影響，提高穩(wěn)定性。(2) 由于細管道的流速測量對時間分辨率要求較高，達ns級別，否則測量結果的精度難以保證。這樣高要求的時間分辨率，硬件電路中難以直接測量，或者測量設備過于昂貴。采用插值后相關法對順、逆流波形進行處理，減小了時間軸間隔，提高了測量精度，節(jié)約了成本。傳統(tǒng)超聲波流速測量方法是使用硬件電路高頻脈沖和專用時間芯片獲取超聲波在管道中傳播的開始和終止時刻后，直接做差得到順、逆流時延。本文截取一個周期順、逆流超聲波波形圖，利用相關函數的算法使圖中所有數據點均參與運算，利用Matlab對超聲信號進行數據處理和分析，得到了較為精確的超聲波順、逆流時延，并且探索了在噪聲環(huán)境下測量結果的準確性。通過以上建立仿真實驗模型和數據分析，為管道流的研究和相關法處理數據提供了依據，并取得較為理想的測量結果。文中還探索了測量穩(wěn)定性與測量精度之間的關系。

1 有限元建模

COMSOL Multiphysics軟件建模流程可以概括為 (1) 聲學模型：調制超聲波的波形以及設置合適的頻率，要求超聲波波峰突出以便于后續(xù)的相關計算。由于管徑小，超聲波傳播路程短，衰減較弱，超聲波頻率適當取較高的頻率，穿透性好，接收波形較為理想。(2) 管道模型：主管道兩端的域設置為吸收層，邊界設置為硬聲場邊界。為了便于超聲波信號傳輸，加裝信號管，其兩端分別為超聲波發(fā)生器和超聲波接收器。(3) 流體模型：設置主管道中流體的各種屬性。(4) 流動模型：采用計算機的模擬流體的方法，模擬流體流動的特點。

由于后續(xù)主要計算超聲波沿順、逆背景流(水流)傳播產生的時間差，即兩個超聲波波形的時延，其本質上為差分運算，信號管對稱且參數完全相同，故信號管對背景流流速測量的影響極小。采用斜式換能器安裝方法可以減少管壁混響而產生的失真。

1.1 管道模型

為減少輸送管道的管壁對超聲波信號的衰減，建模時在輸送管道加裝用于傳播超聲信號的信號管道，輸送管道通徑d為 5 mm，信號管通徑為2 mm，與主管道夾角α=45°，管內壓強設為正常大氣壓，背景流流速均值為U0，信號管穿過主管的距離由L2表示，在任一側延伸的部分，軸向長度為L1(側分支)。因此，信號管的總長滿足L=2L1+L2。由于本文采用了三維對稱的幾何模型，為減少計算量，建立一半管道模型如圖1所示。

圖1 測試管道模型的示意圖Fig.1 Schematic diagram of the test pipeline model

建立仿真模型后，設置計算時間步長為：T0/12即時間步長為 0.033 3 μs，仿真過程的時間范圍為0～30T0，模擬30 MHz采樣頻率采樣，繪制點數為361。使用基于間斷伽遼金(Discontinuous Galerkin,DG)法的物理場接口時，由于過于小的網格使得計算時間過長，故采用自由四面體網格時設置單元質量優(yōu)化，以避免網格過小。網格化后的仿真模型如圖2所示。

圖2 測試管道模型網格化示意圖Fig.2 Schematic diagram of test pipeline model gridding

1.2 聲學模型

建立聲學模型：首先建立信號管兩端的模擬超聲波換能器，既能發(fā)射信號，又能接收信號。因此在實際設備中，這兩端既是發(fā)射器也是接收器。在當前模型中，先僅對向下游傳播的超聲波信號建模，向上游傳播的超聲波信號只需將接收器與發(fā)射器位置互換即可通過計算生成數據。超聲波信號采用信號管底面邊界添加法向加速度a來模擬。

調制超聲波波形解析式為

其中：A=0.1 mm，是信號振幅；頻率=f0=2.5 MHz；角速度ω0=2πf0；信號周期為T0=0.4 μs。此超聲波信號波峰明顯且單一，有利于后續(xù)的相關計算。

1.3 流體模型

為了更加接近真實情況，將背景流狀態(tài)設置為充分發(fā)展的湍流，其雷諾數Re=5×104。利用軟件中“計算流體動力學”(Computational Fluid Dynamics,CFD)模塊中的“湍流，k-ω”接口，對 1.2節(jié)中的模型進行計算。流體模型主要參數如表1所示。

對比兩組患者臨床療效和治療前后銀屑病嚴重程度以及血清干擾素-γ（IFN-γ）、白細胞介素-2（IL-2）水平。（1）臨床療效：無效：療效指數＜30%；一般：療效指數30%～59%；有效：療效指數＞60%[7]；（2）銀屑病嚴重程度：采用銀屑病面積和嚴重程度指數（PASI）評價，得分越低表示銀屑病嚴重程度越低[8]。

表1 流體模型主要屬性表Table 1 Main properties of the fluid model

1.4 流動模型

有限元方法在流體力學中應用時主要采用的就是伽遼金法，伽遼金法采用微分方程對應的弱形式，其原理為通過選取有限多項式函數(又稱基函數或形函數)，將它們疊加，再要求結果在求解域內及邊界上的加權積分(權函數為形函數本身)滿足原方程，便可以得到一組易于求解的線性代數方程，且能夠自動滿足自然邊界條件。超聲波在管道流中傳播的聲壓如圖3所示。

圖3 超聲波在管道流中傳播的聲壓云圖Fig.3 Nephogram of the sound pressure of ultrasonic propagation in pipe flow

2 仿真結果與分析

2.1 相關函數法

利用相關函數法計算兩信號的時延，其基本思想是將其中一個信號以采樣間隔為時間增量，沿著時間軸平移，將平移后的信號與另一個信號做相關運算，當相關函數達到最大值時，信號移動的時間增量即為兩個信號的時延。相關法計算兩信號的時延還具有抗噪聲干擾性質。

設一個信號x(t) 其中包含純凈信號X(t) 和噪聲信號n1(t)，經過延遲后的信號表示為y(t)，其中包含純凈信號Y(t)和噪聲信號n2(t)，表達式為

其中：τ為時延，Δt為采樣間隔，單位均為 s，N為采樣點數。

由于一個信號沿時間軸平移與另一個信號做相關運算時最小移動的時間為采樣間隔，故采樣間隔越短時延測得越精確，即通過互相關計算得到的時延分辨率由采樣頻率決定，但是采樣頻率提高空間有限，采樣頻率太高會導致處理數據時間增加，處理設備儲存量加大等問題[14]。本文將采用插值后相關法提高順逆流接收信號的采樣序列的精度，進而提高時延測量的精度。

本文方法利用采集到的大部分數據來求相關函數得出時延，再對比傳統(tǒng)門限法利用時間芯片獲取超聲波在管道中傳播的起始與終止時刻，相減得時延。由于只用到起始與終止時刻數據，如果起始和終止點數據受噪聲影響大，則結果誤差相對較大，而相關函數法在計算過程用到了大部分數據點，對采樣系統(tǒng)帶來的誤差進行了平均，提高了測量穩(wěn)定性。此外，相關函數法對噪聲還具有一定的抑制作用。

2.2 數據處理

設置背景流為液態(tài)水流，在管道中流動時為充分發(fā)展的湍流。當平均速度為 10 m·s-1時，平均流速云圖如圖4所示。主管道橫截面上的流速分布曲線如圖5所示。

圖4 管道內平均流速云圖Fig.4 Nephogram of the average flow velocity in pipeline

圖5 主管道橫截面上的流速分布曲線Fig.5 Flow velocity distribution curve on the cross-section of main pipeline

將模型上游信號管端設置為信號源，下游信號管端設置為接收端，接收到超聲波順水流傳播的波形，然后將信號源和接收端互換，接收到超聲波逆水流傳播的波形，最后將超聲波波形數據導入Matlab軟件中繪制超聲波沿順、逆流方向傳播的聲壓波形，結果如圖6所示。

圖6 超聲波沿順、逆流方向傳播的聲壓波形Fig.6 Sound pressure waveforms of ultrasonic wave propagation in the downstream and countercurrent directions

COMSOL Multiphysics軟件中導出的接收信號數據繪制的圖像點數為 361，文中模擬由硬件電路采集的原始數據采樣間隔為0.033 3 μs即33.3 ns。如果不采用插值處理，采用文中提及的傳統(tǒng)常規(guī)方法直接計算順、逆流時延，這個時延必然為33.3 ns的整數倍，顯然不滿足測量精度的要求，文中采用10倍插值和100倍插值后時間間隔變?yōu)?.33 ns和0.333 ns，則時延分別為3.33 ns和0.333 ns的整數倍，滿足超聲波在管道中傳播時間的測量需求，而這樣的時間分辨率在硬件電路中難以直接測量，或者測量設備過于昂貴，而對原信號分析處理可以達到所需時間精度，節(jié)約成本。這體現了信號分析的優(yōu)勢處。

根據輸出信號波形圖選用合適的插值方法。Matlab軟件中提供的插值函數interp1( )可以直接使用，并且提供了幾種插值方式：分段線性插值、臨近插值、球面插值、三次樣條插值等。為使圖像更加光滑，采用spline三次樣條插值方法，其插值函數以及其一階和二階導函數都連續(xù)，是最光滑的插值方法。

spline三次樣條插值后的圖像如圖7所示。插值后用相關法求時延時，由于相鄰采樣點間插值函數的方程不同，采樣點數又多，計算相關函數時如果使用連續(xù)函數乘積的積分，計算量就很大。為減小計算量，在完成三次樣條插值后的兩個采樣點之間平均取9個點作為插值點，每個插值點的縱坐標來自于本段的插值函數方程，于是將原采樣時間間隔平均劃分為 10份，則插值后時間軸時間間隔變?yōu)?.33 ns，即10倍插值。采用100倍插值時，原理同上。

流速設置為10 m·s-1時，采用spline三次樣條插值方法插值后的超聲波沿順、逆流方向傳播的聲壓波形如圖7所示，可以看到圖像變得比較光滑，有利于后續(xù)研究分析。

圖7 插值后的超聲波沿順、逆流方向傳播的聲壓波形Fig.7 Sound pressure waveforms of ultrasonic wave propagation in the downstream and countercurrent directions after interpolation

利用 Matlab軟件對發(fā)射和接收的信號做互相關計算：(1) 將兩離散信號點集做互相關運算。(2)檢測互相關函數最大值點，最大值點對應的橫坐標時間即為兩離散信號的延遲τ。

當流體平均流速為 10 m·s-1時，順、逆流超聲波接收信號之間相差139個插值間隔點(延遲點數)，插值間隔為 3.33×10-10s，則τ=139×3.33×10-10s=4.628 7×10-8s。設置流速值為 5、10、15、20 m·s-1時對應時延計算結果如表2所示。

表2 流速、插值間隔點數和時延表Table 2 List of flow velocity, numbers of interpolation intervals and time delay

流速測量值U0可由式(7)得到：

其中：c0為管道流中的聲速，c0=1 481 m·s-1。

當設置不同流速均值分別為5、10、15、20 m·s-1時，由式(7)計算得到的流速測量值如表3所示。

由表3可見，雖然流速不同，但相對誤差均小于3%，與一般工業(yè)儀器4級精度要求相當。

表3 流速測量結果與真值對比Table 3 Comparison between the measured flow velocity and its true value

3 噪聲對測量結果的影響

管道流體的流速測量過程中，外界環(huán)境復雜，基于超聲波非侵入式流量測量方法，易受環(huán)境噪聲影響。為探索噪聲對測量結果的影響，當流速為10 m·s-1時，在接收信號基礎上添加幅值不等的高斯白噪聲?？梢钥吹匠暡ㄐ斡捎谠肼暤募尤胱兊貌灰?guī)則，添加聲壓均值為5 Pa、方差為1 Pa2的高斯白噪聲后超聲波沿順、逆流方向傳播的聲壓波形如圖8所示。

圖8 加噪聲后的超聲波沿順、逆流方向傳播的聲壓波形圖Fig.8 Sound pressure waveforms of ultrasonic wave propagation in the downstream and countercurrent directions after adding noise

對含噪波形進行插值，噪聲對不同插值倍數下的時延點數的影響如表4～表6所示。

表4 流速為10 m·s-1時未插值的噪聲幅值、時延、流速測量值及相對誤差Table 4 List of noise amplitude, time delay, measured flow velocity and its relative error when the true flow velocity is 10 m·s-1 and without interpolation

表5 流速為10 m·s-1時10倍插值的噪聲幅值、時延、流速測量值及相對誤差Table 5 List of noise amplitude, time delay, measured flow velocity and its relative error when the true flow velocity is 10 m·s-1 and with 10 times interpolation

表6 流速為10 m·s-1時100倍插值時的噪聲幅值、時延、流速測量值及相對誤差Table 6 List of noise amplitude, time delay, measured flow velocity and its relative error when the true flow velocity is 10 m·s-1 and with 100 times interpolation

表4表明流速真值為10 m·s-1及未插值時的噪聲幅值、時延、流速測量值及流速相對誤差之間的關系。由于未插值時采樣間隔為0.033 3 μs即時間分辨率為 3.33×10-8s，順、逆流波形延遲點數為 1時相關函數值最大，即延遲時間為3.33×10-8s后計算流速測量值。由表4中可知，在噪聲幅值分別為0、1、5、10 Pa延遲點數均為1，流速相對誤差均為27.62%，未受噪聲影響。

表5表明流速真值為10 m·s-1及10倍插值時的噪聲幅值、時延、流速測量值及相對誤差之間的關系。插值后時間間隔為0.003 33 μs，即時間分辨率為3.33×10-9s。添加幅值為0、1、5 Pa的白噪聲后平均延遲點數為 15，流速相對誤差為 9.555%，添加幅值為10 Pa白噪聲后平均延遲點數變?yōu)?1.5，流速相對誤差為 16.007%，說明幅值為 10 Pa 的白噪聲對測量時延產生影響。

表6表明流速真值為10 m·s-1及100倍插值時的噪聲幅值、時延、流速測量值及流速相對誤差之間的關系。插值后時間間隔為 0.000 333 μs，即時間分辨率為3.33×10-10s。添加幅值為0、1 Pa的白噪聲后平均延遲點數為 139，流速相對誤差為1.522%，添加幅值為5 Pa白噪聲后平均延遲點數變?yōu)?39.2，流速相對誤差為1.668%，說明幅值5 Pa白噪聲對測量時間延時產生微弱影響。添加幅值為10 Pa白噪聲后平均延遲點數變?yōu)?02.2，流速相對誤差為 25.355%，說明幅值為 10 Pa 的白噪聲對測量時延產生較大影響。

為了排除單次實驗帶來的偏差，表中時延點數均為10次實驗的平均值。通過表4、表5和表6中的對比可以得出，利用相關數法計算時延對噪聲具有一定的抑制能力。但是隨著噪聲幅值的增加，插值的倍數越高，噪聲對測量結果的影響越大，測量結果對噪聲越敏感，噪聲對時延產生影響后，噪聲幅值越高會造成測量結果產生偏差越大。

以上實驗結果表明，相關法處理數據，可以對一定幅值的高斯白噪聲產生抑制效果，提高測量準確性。在工程實踐中，已經存在多種較為成熟的抑制白噪聲的方法，這些方法要比單一使用本文的相關法效果要好很多，在實際測量中可以先使用多次平均等手段降低白噪聲對波形的影響，再使用本文的方法進行測量效果會更好，由于本文主要研究相關函數法在細管道中流速測量的優(yōu)勢，對噪聲的抑制是其優(yōu)勢之一，故用其他聯合降噪方法此處不再討論。

4 結 論

本文在前人對管道流速研究的基礎上，針對細管道中流速的超聲測量中的諸多難點如：超聲波在管道中傳播路徑短，時間間隔短，對數據分析精度要求較高，易受噪聲影響等，提出了一種可行的測量方法。本文所做的主要工作是利用 COMSOL Multiphysics仿真軟件建立管道流仿真模型和聲學模型等，采用相關函數法對仿真得到的數據進行分析和計算，其相對于其他數據處理方法精度較高。為測試相關法的抗干擾能力，加入不同均值的高斯白噪聲。實驗表明，相關函數法具有一定抗干擾能力，可以在一定程度上降低噪聲干擾帶來的流量計算誤差。實驗和數據處理結果表明，對細管道中流速的測量時使用插值后的相關法具有穩(wěn)定性并能改善測量精度。值得注意的是，實際測量中，采集卡的時間基準電路自身也存在誤差，采用插值方法提高精度無法超越此時基電路誤差，所以可用的插值倍數是有限的，不能隨意插值獲取高精度。本文方法用于改善超聲波延時估計精度，為細管道中流速的測量提供一種技術參考。但是作為流速測量系統(tǒng)，還有許多其他的誤差來源有待進一步研究。