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

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

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

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

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

1 測量原理

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

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

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

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

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

圖2 測量截面網(wǎng)格點及測量過程示意圖(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

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

2 微風(fēng)速多傳感器數(shù)據(jù)融合方法

2.1 相關(guān)函數(shù)

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

(1)

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

(2)

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

(3)

(4)

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

(5)

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

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

(6)

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

(7)

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

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

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

(8)

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

2.2 基于Kalman濾波算法的多傳感器數(shù)據(jù)融合

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

首先估計狀態(tài)的更新，Φk,n|k-1為狀態(tài)轉(zhuǎn)移矩陣。

(9)

多次更新之后得到最終的估計數(shù)據(jù)Xk,n

(10)

式中，εk,n為每次狀態(tài)更新產(chǎn)生的信息[14-15]；Kk,n為增益矩陣；Hk,nm為n×m矩陣，其元素為風(fēng)速測量數(shù)據(jù)。

(11)

(12)

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

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

(13)

(14)

經(jīng)過上述Kalman濾波方程融合之后得到新的預(yù)測融合傳感器數(shù)據(jù)輸出

(15)

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

3 測量實驗與結(jié)果

3.1 風(fēng)洞數(shù)據(jù)測試

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

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

(a)

風(fēng)洞實驗數(shù)據(jù)表明，數(shù)據(jù)經(jīng)過相關(guān)性判斷和融合后，給出的測量值更接近真實值。圖3表明融合后的測量值誤差更小。

3.2 三維空間流場測量實驗

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

圖4 測量系統(tǒng)示意圖Figure 4.Schematic diagram of measurement system

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

圖5 8475型低風(fēng)速傳感器Figure 5.The 8475 type low wind speed sensor

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

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

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

(a)

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

4 結(jié)束語

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