三維流場測量中的微風量多傳感器數據融合

肖欣招，劉建旭，伍國靖，付東翔

(1.上海理工大學 信息化辦公室,上海 200093;2.上海理工大學 光電信息與計算機工程學院,上海 200093)

三維空間的流場測量方法包括接觸式和非接觸式。接觸式測量通過風速傳感器如熱線風速儀、多孔方向壓力探針等常規傳感器來測量空間內的各點風速，通過數據擬合得到整個空間的流場分布。非接觸式測量包括基于光學方法的激光多普勒測速(Laser Doppler Velocimetry，LDV)法和粒子成像測速(Particle Image Velocimetry，PIV)法[1-3]，以及基于超聲波的超聲多普勒測速(Ultrasound Doppler Velocimetry，UDV)測量方法。LDV應用多普勒效應，利用激光的高相干性和高能量測量流體或固體的流速[4]。PIV是一種基于粒子示蹤和粒子圖像處理算法的流場測量技術[5-7]。UDV基于脈沖超聲回波技術，通過檢測流體中分散的微小示蹤粒子相位變化來實現流體流場的測量[8]。上述非接觸式測量方法通常只能應用于在局部小范圍空間進行三維流場測量，例如發動機內部某個區域、局部管道或者閥門內部等[9]。因此，基于接觸式測量的三維流場測量仍然是一種常用的測量方法。

本文針對三維流場測量中的多傳感器數據融合進行研究。該測量對象為有人工作環境的封閉艙室，微風速流場分布測量可為其通風系統的風速、風向及風機功率等參數提供設計依據。較大室內空間尺寸(例如10 m×2 m×2 m)由于測量空間較大，故采用多只傳感器組成的傳感器陣列進行風速測量。在這種測量方式中，雖然每個傳感器經過了標準數據的標定，但實際使用中仍然會出現無效數據。微風量風速傳感器的量程較低，基于單純最小二乘法建立的線性模型也無法完全反映傳感器的實際測量特性[10]。此外，多只傳感器數據需要融合，同一測量點微風速條件下的多個傳感器測量值跳變范圍大，無效數據時有出現。因此排除無效數據、融合有效數據是保證微風速流場測量值準確的關鍵。

本文提出一種基于相關性函數-卡爾曼濾波算法的多傳感器數據融合算法。在上述三維微風速流場測量中，對多個風量傳感器進行數據測量，采用相關性判定方法，剔除測量異常數據和無效數據。基于卡爾曼濾波方法[11]，將傳感器標定的輸出數據和方差作為初始預估值和初始方差估計，進行多個風量傳感器數據的融合。與普通標定的傳感器測量數據相比，新方法提高了風量數據的測量精確。測量實驗驗證了該方法的有效性，說明使用該方法可確保三維空間微風速流場測量數據的精確性和可靠性。

1 測量原理

根據測量要求，測量系統是對一個封閉的類長方體空間進行各個點的風速測量，采集空間內的各點風速數據，經過數據處理即可得到各個截面的二維風速流場以及整個三維空間的流場分布，為該封閉艙室的通風系統設計提供依據。因此測量方案需在測量空間中劃分出一系列測量點。如圖1所示，這些點上限至天花板附近，下限至地板附近，左、右限接近墻壁。通過分別測量這些網格點風速，可測定出測量空間流場分布。

圖1(a)是一個網格測量點截面，整個測量空間被劃分成很多個截面，需對多個位置的點進行測量。如果每點都布置傳感器，一方面需要大量傳感器，增加了系統成本；另一方面密布的大量傳感器及安裝支架會對空間流場產生較大影響，使得測量結果無法反映真實的測量空間流場分布。因此，測量方案采用傳感器支架，搭載一定數量的風速傳感器和三維移動平臺，通過移動掃描測量方式來采集各個空間點的風速。

圖1 測量空間、測量點及傳感器分布圖(a)測量空間測量點網格劃分 (b)H型支架 (c)風速傳感器角度Figure 1.Measurement space,points and sensor distribution map(a)Meshing of measurement points (b)H style stand (c)Wind speed sensor angle

測量空間網格點劃分如圖1(a)所示，圖中圓點為測量位置點。傳感器支架如圖1(b)所示，呈H型，共8個端點(見圖中A～H位置)。每個端點上呈垂直角度布置2個熱球式風速儀，如圖1(c)所示。這里需要考慮風向的因素，根據測試要求，風速傳感器選擇TSI 的全向熱球式風速傳感器。雖然是全向傳感器，但實際上熱球式也有一定的方向性。由于測試環境內各個點的風向是不確定的，傳感器的布置和測量數據處理方法對保證風速測量的準確性至關重要。

測量以單個截面為單位，依次對每個截面進行掃描測量。圖2為單個測量截面中傳感器選擇測量點的示意圖。如圖2(a)所示，傳感器節點A～D在該截面內。測量根據就近原則，單個傳感器節點選擇離它最近的位置點測量，例如傳感器頂點D測量第1個測量點。本文規定先在垂直方向測量，下一個測量位置是測量點2，對該點的測量可以由傳感器節點A或D來進行。如圖2(b)所示，節點A距離該點最近。此外，還需要判斷傳感器支架向下運行是否有限位，例如在臨近地板的測量點。如果支架往下移動，可能會觸碰到垂直方向下方的限位擋板。因此，如果傳感器節點A不受限位條件限制，則測量點2由節點A測量，否則由節點D測量，如圖2(c)所示。

圖2 測量截面網格點及測量過程示意圖(a)測量位置1 (b)測量位置2 (c)測量位置3Figure 2.Schematic diagram of measurement section grid points and measurement process(a)Positon 1 (b)Positon 2 (c)Positon 3

本研究通過機電控制系統和上位機軟件系統采集測量點的移動和數據。機電系統的主要功能是通過控制電機驅動該傳感器支架平臺做三維移動至各個測量點進行測量。16只風速傳感器的測量值由上位機軟件系統讀取和處理。

2 微風速多傳感器數據融合方法

2.1 相關函數

對多個風速傳感器測量同一個坐標測點的風速數據，記第i個傳感器和第j個傳感器測得的數據分別為Yi和Yj，均服從高斯分布。以它們的概率分布函數作為傳感器的特征函數，記成pi(y)，pj(y)。yi、yj分別為Yi和Yj的一次觀測值。為評估yi和yj之間偏差程度，計算兩者之間的置信距離測度dij[12]

(1)

同理，yj和yi之間的偏差程度置信距離測度dij為

(2)

根據測量誤差概率積分函數Φ(t)計算上述置信測度dij和dji[9]。

(3)

(4)

n個風速傳感器測量同一個坐標測點的風速數據，置信距離測度dij(i,j=1,2,…,n)構成多個風速傳感器數據的置信距離矩陣Dn，如式(5)所示。

(5)

根據上述定義，dij值為概率密度函數是正態分布的兩個測量值[yi,yj]范圍內的積分，0≤dij≤1。該值越小說明第i個風速傳感器與第j個風速傳感器測量輸出越接近，第j個傳感器對第i個傳感器的支持程度越高。為消除傳感器的無效數據和多個傳感器數據的融合，根據模糊理論中的相關性函數定義，可令

f(i|j)=1-dij，i,j=1,2,…,n

(6)

相關性函數f(i|j)的大小表示傳感器i被傳感器j支持的程度,相關性函數定義為

(7)

構造由f(i|j)值組成的矩陣J，其為方陣且秩為n，i,j=1,2…,n。

若f(i|j)=1，則認為第i個風速傳感器與第j個風速傳感器相融性好，稱第i個風速傳感器支持第j個風速傳感器。若f(i|j)=f(j|i)=1，稱第i個風速傳感器與第j個風速傳感器互相支持。如果一個風速傳感器被一組風速傳感器所支持，則這個傳感器的讀數是有效的。若一個風速傳感器不被其它風速傳感器所支持，或只被少數的風速傳感器所支持，則這個風速傳感器的讀數是無效的，進行數據融合時，需刪掉該數據，保證數據有效性的客觀判斷。

為確定每個風速傳感器被其它傳感器支持的程度，令

(8)

傳感器支持度矩陣Ji給出了第i個風速傳感器被其它風速傳感器支持的程度。通過計算處理排除掉支持度為0或支持程度低的傳感器數據，得到最終的多個傳感器的有效測量數據。

2.2 基于Kalman濾波算法的多傳感器數據融合

設有m個風速傳感器，通過風洞標定的每個傳感器風速數據Yk=(y1,y2,…,yk)，對m個風速傳感器風速數據進行融合。當前時刻m個風速傳感器測量數據為Yk,m=(yk,1,yk,2，…，yk,m)。隨機噪聲和測量噪聲取高斯噪聲，隨機噪聲和測量噪聲協方差分別為Qk,m和Rk,m。用Xk,n|j表示根據j時刻及j以前時刻的傳感器測量值對k時刻狀態Xk,n的估計值，得出n個估計輸出，即可得到風速傳感器數據融合的Kalman濾波方程[13]，具體過程如下：

首先估計狀態的更新，Φk,n|k-1為狀態轉移矩陣。

(9)

多次更新之后得到最終的估計數據Xk,n

(10)

式中，εk,n為每次狀態更新產生的信息[14-15]；Kk,n為增益矩陣；Hk,nm為n×m矩陣，其元素為風速測量數據。

(11)

(12)

最終得到估計方差矩陣Pk,n。Kk,n|k-1,n為預測方差矩陣[16-18]

Pk,n=(I-Kk,nHk,n)Pk,n|k-1,n

(13)

(14)

經過上述Kalman濾波方程融合之后得到新的預測融合傳感器數據輸出

(15)

融合后，最終估計方差Pk,n=(Pk,1,Pk,n,Pk,2,…，Pk,n)。

3 測量實驗與結果

3.1 風洞數據測試

為驗證上述數據相關性判斷和融合方法，選取3只相同型號的風速傳感器在風洞中進行實驗。本文所使用的傳感器型號為TSI 公司的熱球式8475型傳感器。該傳感器主要特點為：(1)全向探頭尖端；(2)精確測量風速范圍0.05～2 m·s-1；(3)適合未知或不同的流向，輸出為4～20 mA信號。

設定風洞理論風速0.6 m·s-1作為標準風速數據。每個風速傳感器在該標準風速下同時采集10組，每組10個數據。3個傳感器得到30組測點的風速數據。以30組測點風速數據作為樣本，先對風速數據進行相關性判定，剔除部分無效數據點和支持程度低的傳感器數據點，再對余下的有效數據點運用上述Kalman濾波方法進行數據融合。融合后數據的風速對比和數據均方根誤差(Root Mean Square Error，RMSE)對比結果如圖3(a)和圖3(b)所示。

(a)

風洞實驗數據表明，數據經過相關性判斷和融合后，給出的測量值更接近真實值。圖3表明融合后的測量值誤差更小。

3.2 三維空間流場測量實驗

上述傳感器風洞測量數據融合驗證了多傳感器測量數據融合的有效性，下文將進行三維空間的流場測量實驗。本研究中封閉空間大小(x×y×z)為5 000 mm×1 700 mm×1 800 mm，測量系統結構示意圖如圖4所示。該系統由便攜式工業PC、風速傳感器陣列、三維移動測量平臺及其H型支架、導軌和信號與電源箱及拖鏈電纜等主要部分組成。16只風速傳感器陣列布置在H型支架上，搭載于三維移動測量平臺上，隨平臺三維移動。

圖4 測量系統示意圖Figure 4.Schematic diagram of measurement system

圖4中的每個傳感器節點由上述兩個呈垂直角度布置的兩個熱球式風速8475型傳感器組成，風速傳感器如圖5所示。電源箱內布放兩個24 VDC電源，分別為所有傳感器和3個步進電機供電。信號線纜和電源電纜置于拖鏈內，當測量平臺移動時，拖鏈保證了線纜的拉伸不受損壞。測量平臺每次移動到測量點，為了消除平臺移動導致的氣流擾動，在每個測量點平臺都要停留若干秒鐘，該值可以在軟件系統中調整。

圖5 8475型低風速傳感器Figure 5.The 8475 type low wind speed sensor

測量過程中，風速傳感器位置移動如前所述，傳感器陣列共16只傳感器，分別布置在A～H共8個位置上，為8個節點。如圖1(b)所示，每個節點上有兩個傳感器，傳感器編號S={s1,…，si，…,s16}。

測量步驟如下：(1)對每個測量網格點p，距離最近的傳感器i移動至點p，讀取傳感器節點i的兩個測量值si和si+1,1≤i≤16，該點的傳感器數據支持矩陣J由si和si+1決定；(2)根據傳感器支持度判斷該兩個傳感器是否互相支持，如果si和si+1相互支持，則由si和si+1融合得到該點的風速值；如果兩個傳感器相互不支持，則由其它7個傳感器節點中距離點p最近且滿足運動至該點條件的一個傳感器節點運動至點p。設該節點為j,其測量數據為si和si+1,1≤i≤16，這里j≠i。該點的傳感器數據支持矩陣J由si、si+1、sj、sj+1決定。如果si、si+1、sj、sj+1這4個傳感器值相互不被支持的傳感器超過2個，則再由其它6個傳感器節點中最近的而且滿足運動條件的一個傳感器節點加入測量。這樣該點的傳感器數據支持矩陣J由6個傳感器數據決定。以此類推，同一個測量點的測量至少被測量2次，最多可以16次，即16只傳感器全部加入測量；(3)根據傳感器數據支持矩陣J，接下來對測量數據進行判斷，排除無效數據，然后對有效數據進行融合，給出測量點的最終風速值。

測量系統對測量空間進行了掃描測量，以X方向間隔為200 mm，測得10個yoz截面的風速數據，如圖6(a)所示。x=450 mm，x=850 mm，x=1 250 mm處3個yoz截面風速流場如圖6(b)～圖6(d)所示。

(a)

由于測量空間大，空間網格測量點數量較多，如果每個網格點都采用8個傳感器節點和16個傳感器進行測量，會大幅度延長總測量時間。在實際測量中，大部分測量點的測量次數為2～4次。這種測量方式還可以較快發現某個傳感器的故障，及時將其屏蔽，使其不參與測量。

4 結束語

為滿足某特殊裝備的艙內低風速流場測量精度的要求，本文研究了基于接觸式測量的多傳感器測量數據融合方法。針對風速傳感器陣列的測量數據處理，提出一種基于相關性函數-卡爾曼濾波算法的多傳感器數據融合算法，對多個傳感器的測量數據，運用相關性判定方法，排除測量數據中的異常數據和無效數據；把傳感器標定輸出數據和方差作為卡爾曼濾波的初始預估值和初始方差估計，進行多傳感器數據的融合。與普通傳感器標定的數據相比，該方法融合得到的風速數據測量誤差更小。測量數據融合實驗結果表明，該方法能有效提高測量精度，基于該數據處理方法的三維空間流場測量結果準確可靠。