基于機器視覺的流體速度量化研究

陳金鵬 李貴山 李雙科 馬宏鋒

(蘭州工業(yè)學(xué)院電氣工程系，甘肅 蘭州 730050)

0 引言

不管是普通的流體視覺化，還是由計算機輔助的流體視覺化，它們在流體動力學(xué)中對于觀察和監(jiān)控流體結(jié)構(gòu)都是非常有用的工具。流體視覺化通過在流體中引入染色物或煙氣形式的污染物(即被動示蹤體)，使流體的結(jié)構(gòu)變得可視。

根據(jù)不同流體的密度，只有選擇合適的污染物，才能確保很好地跟蹤流體的流動，視覺化才能集中于被動示蹤體以及它在時間和空間上的濃度的變化，從而識別和確認(rèn)特定流體的特性。

對于用數(shù)碼相機采集的流體/被動示蹤體流動的圖像，可以進一步對其進行處理，以用于不同目的的研究，如渦流分離分析［1］、回流區(qū)位置確認(rèn)［2］、監(jiān)控氣蝕結(jié)構(gòu)［3］、流體流動對沉積的影響［4］等。在這個階段出現(xiàn)的問題是能否借助計算機視覺獲取流體的定量特性，也就是說，能否采用被動示蹤體將速度場這樣的流體運動量從一組流體圖像中抽取出來。

目前，流體速度場測量主要采用粒子圖像速度測量儀。該測量儀通過在流體中引入固體微粒，采用機器視覺的方法進行測量。借助機器視覺對流體進行定性和定量的研究，至少應(yīng)同時使用兩種方法或系統(tǒng)［5－6］，但這樣的組合系統(tǒng)往往很復(fù)雜。此外，大部分普通的測量方法都不能測量非平穩(wěn)或周期性流體結(jié)構(gòu)的流體速度。

本文僅通過采用視覺方法對加入示蹤體的流體的速度場進行測量，并在實際使用中檢驗這種方法的可行性。該方法是利用一種已知的物理關(guān)系，把污染物濃度和流體速度結(jié)合起來的方法。

1 研究方法

在所描述的研究方法中，最重要的步驟是把速度向量場與濃度的標(biāo)量場聯(lián)系起來，后者用灰度圖像表示。通過濃度N(即單位體積流體中污染體的分子數(shù))表征污染體引入流體的情況。如果采用如下速度供給污染體分子，則可以使x方向上N的梯度保持不變。

式中:Φ為X方向的擴散流量;D為分子擴散系數(shù)。

式(1)為菲克擴散定律。如果考慮非穩(wěn)態(tài)擴散過程，基本體積中濃度N的變化率可以用流出和流入該體積流體邊界的污染體的微分表示，則式(1)可以改寫為［7］:

當(dāng)有流體運動(即對流)時，用總的導(dǎo)數(shù)代替時間偏導(dǎo)數(shù)，則式(2)可進一步展開為:

式(3)即為著名的對流擴散方程，它表示示蹤體/污染體與流體動力學(xué)的基本關(guān)系。其中，等式左邊第二項包括流體對流的影響。這樣，如果濃度N的值(即N的導(dǎo)數(shù)已知)和流體中污染體的擴散系數(shù)D已知，則唯一的未知量就是速度v，即只有xi方向上的分量vi為未知量。

分子擴散系數(shù)D是物質(zhì)的屬性。在一定壓力和溫度下，D與另一種物質(zhì)有關(guān)，污染體濃度N在不斷變化，因此，應(yīng)該即時測量。在這種情況下，應(yīng)該將它的值從流體灰度圖像中抽取出來。

引入跟蹤物的流體灰度圖像可形成具有不同灰度的像素的矩陣。一般來說，灰度的整數(shù)值從0(黑)到255(白)，它可以在0～1區(qū)間被歸一化。如果引入流體的污染物是用外光源照明的，那么它的灰度值一般比流體的灰度值高。流體中污染體的濃度越高，在灰度圖像中污染體的灰度值越高。因此，有:

式中:A為灰度圖像中選定像素窗口的平均灰度值。文獻［8］和文獻［9］對此已有論述，式(4)中灰度值和污染體濃度的比值是線性的。選定的像素窗口的平均灰度可以用如下關(guān)系［2］得到:

式中:E(i，j)為l×m像素窗口中實測的(i，j)點單個像素的灰度強度。

濃度的時間導(dǎo)數(shù)?N/?t可以利用如圖1所示的連續(xù)圖像中固定窗口的平均灰度計算得到。

圖1 流體連續(xù)圖像窗口Fig.1 Continous image windows of the fluid

利用已知的兩連續(xù)圖像的時間間隔Δt(通過相機的圖像采集頻率得到)，則?N/?t項可以用式(6)估算:

式中:ΔA為固定窗口中相距Δt的兩連續(xù)圖像之間平均灰度的差值。污染體濃度在空間上的偏導(dǎo)數(shù)?N/?xi可以通過計算流體的單個圖像獲得。為了得到平均灰度A在空間上的差值，將現(xiàn)在的窗口向選定的xi方向移動，其示意圖如圖2所示。

圖2 在窗口中移動單個圖像示意圖Fig.2 Diagram of moving single image in the window

空間導(dǎo)數(shù)可以用標(biāo)準(zhǔn)數(shù)字技術(shù)完成。式(7)為忽略截斷誤差的中心差分法［10］:

式(3)在兩維系統(tǒng)中表示兩個線性微分方程，它有兩個未知數(shù)，即兩個速度分量。求解這樣一組方程需要明確的初始條件，非常繁瑣和困難，因為流體的運動特性一般都未知。但是，在速度分量和其空間導(dǎo)數(shù)仍未知的情況下，這組微分方程可以轉(zhuǎn)換為線性方程組。這組方程即使是超定的(即方程的個數(shù)多于未知數(shù)的個數(shù))，也可以通過用數(shù)字的方法很容易地解出。與微分方程組相比，線性方程組的未知數(shù)的個數(shù)增加了一倍，這就意味著對兩維系統(tǒng)至少需要四個一般方程。每一個方程需要兩個流體連續(xù)的圖像(因為要計算時間差)，即至少需要五個連續(xù)的流體圖像才能得到流體的兩個速度分量的信息。

上述研究方法是基于已知的物理關(guān)系，反映了借助機器視覺確定流體動力學(xué)特性的基本原理。它的有效性可以用傳統(tǒng)的測量方法驗證。為了達到這個目的，采用熱絲風(fēng)速計進行試驗驗證。

2 試驗結(jié)構(gòu)

測量來自垂直管空氣射流速度場的測量裝置如圖3所示。

圖3 測試部分的試驗裝置Fig.3 Experimental device of the test section

流體的雷諾數(shù)在管道的出口處進行測定，記為Re=1300。由于管道中的流體是層流，可以采用熱線風(fēng)速計法(hot-wire anemometer，HWA)和計算機輔助視覺法(computer-aided visualization，CAV)完成測量。在選定的測量點，用兩分支熱絲風(fēng)速計測定速度的分布。兩分支熱絲風(fēng)速計測量流向(y)和側(cè)向(x)瞬時速度波的步驟與 Bruun［11］和 Jφrgenson［12］的研究步驟一致。在9×17網(wǎng)孔中選定153個測量點，每兩個點之間的距離用定位表設(shè)置為5 mm。

為了用熱絲風(fēng)速計測量流體速度，本文采用DANTEC 55P62恒溫?zé)峤z風(fēng)速計傳感器系統(tǒng)。該傳感器用鍍鉑鎢絲制成，長度和直徑分別為1.25 mm、φ5 μm。風(fēng)速計放大器的截止頻率設(shè)置為10 kHz，風(fēng)速計的工作溫度為250℃。風(fēng)道中的兩分支熱絲風(fēng)速計傳感器的位置要求能從兩個方向測量瞬時速度。因此，風(fēng)速計傳感器采用PC控制的精確定位裝置進行定位。熱絲風(fēng)速計的信號采用NI公司的16位數(shù)據(jù)采集板獲得。數(shù)據(jù)采集時間為6 s，采集頻率為50 kHz。

風(fēng)速計的輸出在傳送到具有SCXI模塊的A/D轉(zhuǎn)換器之前，先用頻率為10 kHz的低通四階貝塞爾濾波器予以濾波。系統(tǒng)采用LabVIEW軟件采集并存儲數(shù)據(jù)。熱絲風(fēng)速計的校準(zhǔn)是在專門校準(zhǔn)風(fēng)速計的測量臺上完成的。為了與國王定律［11］保持一致，對第一和第二個熱絲的常數(shù)設(shè)定如下:

此外，在校準(zhǔn)和測量過程中還對溫度進行了補償。

溫度的測量采用A級四線電阻溫度計Pt100和Agilent 34970A數(shù)據(jù)采集器。根據(jù)國王定律方程式［11］和Jφrgenson方程［12］，從熱絲風(fēng)速計的輸出計算得到實際的流向流速vy和側(cè)向流速vx。測量也是用同樣的設(shè)備完成的，且在測量之前進行了校準(zhǔn)［13］。測量的總誤差［14］是由工作點的選擇和波動的不確定性、定位起始點的選擇、風(fēng)扇的轉(zhuǎn)速、定位表對熱絲溫度傳感器的定位、線性化、A/D的分辨率、數(shù)據(jù)采集、溫度補償、濕度以及熱絲風(fēng)速計的頻率響應(yīng)有限等因素引起的。風(fēng)扇由一臺分離變壓器供電，便于與測量設(shè)備的供電電源分開。瞬時速度的整體測量誤差估計為測量值的2.8%。

在CAV測量中，為了使流體可視化，在充分發(fā)展流中引入含有氣化石蠟油的被動示蹤煙氣［15］。測量點與上述熱絲風(fēng)速計的測量點完全相同。通過流體入口將被動示蹤體引入管道。相機放置于與風(fēng)道的測量段垂直并且距離測量段表面1 m的位置。相機的圖像采集頻率為400 Hz，拍攝圖像的快門速度為0.6 ms，所記錄的圖像具有8位灰度值、600×250像素的分辨率。本文采用具有直線導(dǎo)光的不間斷無閃爍光源(Vega Velum DC 150 W)進行照明;利用LabVIEW軟件包實現(xiàn)對相機的設(shè)置、圖像的采集和保存。在每一個測量點采集1000個連續(xù)的圖像，圖像的個數(shù)受被動示蹤體煙氣的發(fā)生時間限制。被動示蹤體瞬時濃度的總體測量誤差估計為測量值的 3.8%［15］。

3 試驗結(jié)果

試驗中采集的圖像序列用Matlab程序包進行處理。兩個速度分量的測定與試驗測量(熱絲法)的情況一樣，都是在同樣的測量點、在x-y平面利用機器視覺測定的。根據(jù)文獻［16］，空氣(T=23℃、p=100 kPa)中的癸烷(C10H22)蒸汽 D=7.01 ×10－12m2/s。用 30 個連續(xù)圖像序列來確定x-y平面上的流速場，序列中圖像的數(shù)量保證了方程式(3)即使在被動示蹤體稀少(測量平面x-y外沿)的情況下，對所有的測量點都有效。另一方面，當(dāng)圖像序列很大時，求取速度分量的平均值需要非常多的計算時間，有可能導(dǎo)致給定采集頻率上的信息丟失。為了減小導(dǎo)數(shù)數(shù)值計算的截斷誤差，式(7)中步長Δxi應(yīng)選得盡可能小。因此，步長接近于1個像素，它的大小取決于圖像分辨率。

本文采用HWA和CAV兩種方法測量流速場。測量結(jié)果顯示，在大部分測量點，HWA法和CAV法的測量結(jié)果一致。

兩種方法的測量結(jié)果表明，離開節(jié)流孔的射流流速減小，射流外邊沿的流體沿徑向朝外側(cè)運動(射流沿軸向和徑向散開)。盡管如此，兩種方法所測得的速度矢量的方向和大小都有差別，尤其在射流的外邊沿。為了正常工作，序列中的每個圖像都要有足夠的被動示蹤體濃度值，而在這些位置，由于石蠟油蒸汽很少，影響了CAV法的效果。

對兩種方法測得的速度分量分別從流向方向(y)和正常方向(x)進行定量分析。用x方向和y方向的歸一化值來表示結(jié)果，這里 x*、y*∈［0，1］。x*=0 和 x*=1分別對應(yīng)于圖3最左邊和最右邊的測量點;類似地，y*=0和y*=1分別對應(yīng)于距離管道節(jié)流孔最遠(yuǎn)和最近的測量點。

經(jīng)歸一化流向，當(dāng)觀測點y*的值分別為0.125、0.375、0.625 和 0.875 時，兩種測量方法所測得的流速分量vy的值如圖4所示。

圖4 不同流向觀測點流速分量vyFig.4 Velocity component vyin different directional observing points

通過比較兩種測量結(jié)果可知，不管采用何種測量方法，y方向的流速具有相似的外形。但是，采用CAV法時vy的波動比HWA法的測量結(jié)果要大得多，在遠(yuǎn)離管道節(jié)流閥處表現(xiàn)得更為突出。這種現(xiàn)象很可能是流體的漩渦和特殊測量方法的數(shù)據(jù)采集頻率共同作用的結(jié)果。當(dāng)選定雷諾數(shù)，且觀測點y*＞0.3時，管道中的流體形成了強烈的漩渦。

HWA法的數(shù)據(jù)采集頻率大約為50 kHz，在每個點的測量歷時約2 s。因此，HWA法獲得的結(jié)果是這段時間的平均值，這個時間足夠長，足以保證大漩渦的旋轉(zhuǎn)對流體速度的平均值影響不大。對于同樣的時間段，一方面CAV法大約需要700個圖像序列，這將極大地影響處理/計算時間;另一方面，這種情況下CAV法的采集頻率僅為HWA法采集頻率的1/140。為了獲得更精確的結(jié)果，應(yīng)該增加圖像序列的數(shù)量，但該方法只適用于流體平均速度(波動)偏差較大、持續(xù)時間足夠長且被CAV法檢測到的情況。即與HWA法相比，CAV法需要相機具有較高的采集頻率，兩者才有更好的可比性。

根據(jù)HWA法，圖4中每一個流向觀測點的vy(沿側(cè)向距離x的計算值)平均值的相對差如表1所示。

表1 vy平均值的相對差Tab.1 Relative difference of vyaverage values

表1表明四個選定的流向觀測點的平均值相對差大約為10%。因此，可以得出結(jié)論:與HWA法相比，即使在圖像序列很短的情況下，CAV法也能很好地預(yù)測vy的值。

根據(jù)CAV法和HWA法，vx平均值的相對差比較如表2所示。

表2 vx平均值的相對差Tab.2 Relative difference of vxaverage values

表2中，vx平均值相對差比表1中vy平均值相對差稍大一些。其主要原因是:與vy的值相比，vx值更小(接近于0，最多到0.04 m/s)。由于相對差再大也不會超過15%，所以，對基本CAV法而言，這可以認(rèn)為是較好的結(jié)果。

4 結(jié)束語

本文闡述了借助于機器視覺對流體動力學(xué)性能量化分析的基本原則。通過流體中污染體(被動示蹤體)的連續(xù)圖像序列，測量了流體的兩維速度場。這是一種通過對流擴散方程把速度場和濃度場聯(lián)系起來的方法。所需參數(shù)的計算是用圖像處理和普通的數(shù)字技術(shù)來實現(xiàn)的。該方法本身是用傳統(tǒng)的熱絲風(fēng)速計法與機器視覺法同時進行而予以驗證的。

兩種測試方法的結(jié)果對比表明:視覺法(即使最基本的形式)能夠相當(dāng)精確地計算出流體的大小和方向;兩種方法在不同流向觀測點的速度分量外形是相似的，兩種結(jié)果的差別非常小。由此可以得出利用機器視覺測量流體速度場是可行的結(jié)論。

從現(xiàn)階段來看，視覺法的精度劣于傳統(tǒng)方法如HWA、PIV和LDA。因此，未來的研究不應(yīng)該僅僅專注于數(shù)字計算領(lǐng)域(如步長大小、不同導(dǎo)數(shù)計算方法的應(yīng)用對最終結(jié)果的影響等)，而是更應(yīng)該專注于測量準(zhǔn)備工作領(lǐng)域(如照明、三維流體結(jié)構(gòu)、非定常性、序列中圖像的個數(shù)等對結(jié)果誤差的影響)。根據(jù)現(xiàn)有的算法和設(shè)備要求(合適的圖像分辨率和采集頻率)，視覺方法僅可應(yīng)用于有很高雷諾數(shù)的湍流，具有一定的有效性和局限性。

［1］Sirok B，Potocar E，Novak M.Analysis of the flow kinematics behind a pulsating adaptive airfoil using computer-aided visualization［J］.Strojniski Vestnik-Journal of Mechanical Engineering，2000，46(6):330－341.

［2］Sirok B，Bajcar T，Dular M.Reverse flow phenomenon in a rotating diffuser［J］.Journal of Flow Visualization and Image Processing，2002，9(2＆3):193 －210.

［3］Sirok B，Blagojevic B，Bajcar T，et al.Simultaneous study of pressure pulsation and structural fluctuations of a cavitated vortex core in the draft tube of a Francis turbine［J］.Journal of Hydraulic Research，2003，41(4):541 －548.

［4］Bajcar T，Blagojevic B，Sirok B，et al.Influence of flow properties on a structure of a mineral wool primary layer［J］.Experimental Thermal and Fluid Science，2007，32(2):440 －449.

［5］Hu H，Saga T，Kobayashi T，et al.Analysis of a turbulent jet mixing flow by using a PIV-PLIF combined system［J］.Journal of Visualization，2004，7(1):33 －42.

［6］Reungoat D，Riviere N，F(xiàn)aure J P.3C PIV and PLIF measurement in turbulent mixing［J］.Journal of Visualization，2007，10(1):99 －110.

［7］McComb W D.The physics of fluid turbulence［M］.Oxford:Clarendon Press，1996.

［8］Aider J L，Westfreid J E.Visualizations and PDF of the fluctuations of a passive scalar in a turbulent Goertler flow［J］.Experimental and Numerical Flow Visualization，1995，218:123 －130.

［9］Simoens S，Ayrault M.Concentration flux measurements of a scalar quantity in turbulent flows［J］.Experiments in Fluids，1994，16(3 －4):273－281.

［10］Mathews J H，Kurtis D F.Numerical methods using Matlab［M］.4th Edition.London:Prentice-Hall，2004.

［11］Bruun H H.Hot-wire anemometry-principles and signal analysis［M］.New York:Oxford University Press，1995.

［12］Dantec Dynamics A/S.How to measure turbulence with hot-wire anemometers-a practical guide［EB/OL］.［2002 －02 －01］［2012 －03 － 02］.Http://www.dantecdynamics.com/Admin/Public/Download.aspx?file=files%2fsupport_and_download%2fresearch_and_education%2fpracticalguide.pdf.

［13］Coleman H W，Steele W G.Experimentation and Uncertainty Analysis for Engineers［M］.New York:John Wiley ＆ Sons Inc.，1999.

［14］Eberlinc M，Sirok B，Hocevar M，et al.Numerical and experimental investigation of axial fan with trailing edge self-induced blowing［J］.Journal of Engineering Research，2009，73(3):129 －138.

［15］Eberlinc M，Dular M，Sirok B，et al.Influence of blade deformation on integral characteristic of axial flow fan［J］.Journal of Mechnical Engineering，2008，54(3):159 －169.

［16］Won D，Shaw C Y.Determining coefficients for mass-transfer models for volatile organic compound emissions from architectural coatings［C］//Proceedings of the 7th International Conference of Healthy Buildings，2003:1－6.