一種應(yīng)用于水流的紋影特性光流測速算法

黃天立，王 倩

上海交通大學(xué) 機(jī)械與動力工程學(xué)院，上海 200240

0 引言

作為流動類型判斷、壓力計算等的基礎(chǔ)數(shù)據(jù)，流速是流體力學(xué)研究中最為重要的參數(shù)之一。流速測量技術(shù)是開展流體相關(guān)實(shí)驗(yàn)研究和工程應(yīng)用的重要技術(shù)手段，對于研究基本物理現(xiàn)象和深入解釋流動機(jī)理具有重要意義。在過去幾十年里，流速測量技術(shù)不斷發(fā)展，經(jīng)歷了從介入式到非介入式、從單點(diǎn)到平面再到全場、從低速測量到超聲速測量的過程[1]。傳統(tǒng)的流速測量主要為介入式測量，包括皮托管（Pitot tube）、熱線熱膜風(fēng)速儀（HWFA,Hot Wire Film Anemometer）等，不足之處是會對流場產(chǎn)生較大的擾動。近年來發(fā)展了基于圖像處理技術(shù)的流速測量方法，包括粒子圖像測速（PIV,Particle Image Velocimetry）、基于紋影成像的速度測量等。PIV 是目前最常用的全場測速技術(shù)，通過在流場中布撒示蹤粒子，以脈沖激光片光源入射所測流場區(qū)域，采用相機(jī)連續(xù)兩次或多次采集粒子圖像，通過互相關(guān)算法求得流場速度分布。在示蹤粒子濃度較低時，也可對單個粒子進(jìn)行測量，即粒子跟蹤測速技術(shù)（Particle Tracking Velocimetry,PTV）[2]。紋影成像基于折射率與密度梯度的關(guān)系將透明工質(zhì)流動可視化，是一種常用的流場觀察方法。近年來，有學(xué)者結(jié)合互相關(guān)算法和光流算法，探討了利用紋影圖像測速的可能性。與PIV 相比，紋影圖像測速以普通光源代替激光光源，節(jié)省了成本，且無需添加示蹤粒子，可以實(shí)現(xiàn)無擾動測量[3]。對于可壓縮湍流邊界層內(nèi)的流動，因難以添加示蹤粒子無法使用PIV 方法，采用紋影圖像測速則可獲得較好結(jié)果，已被應(yīng)用于激波、邊界層、帶電流場及腐蝕性流場等的定性和定量研究中[4-6]。此外，紋影成像測速具有高幀頻采樣的獨(dú)特優(yōu)勢，僅需普通光源（如LED、氙燈等）即可獲得微秒級曝光時間，可采用高速相機(jī)的最高幀頻。因此，研究紋影圖像測速方法有助于推動流場測速技術(shù)的發(fā)展和測速方法的推廣應(yīng)用[7]。

基于紋影圖像的測速方法主要分為互相關(guān)法和光流法兩大類?；ハ嚓P(guān)法將紋影圖像中的顯著運(yùn)動結(jié)構(gòu)（如渦團(tuán)、邊界）作為追蹤對象，運(yùn)用互相關(guān)算法計算速度場。Kegerise 等[8]進(jìn)行了點(diǎn)熱源的軸對稱湍流羽流實(shí)驗(yàn)，通過PIV 方法測得中心面速度，并對互相關(guān)算法的積分效應(yīng)進(jìn)行了評估。結(jié)果表明，PIV 方法和互相關(guān)法的歸一化平均流速剖面一致性較好?；ハ嚓P(guān)法的偏差誤差有兩項(xiàng)，一是由窗口的平均效應(yīng)引起，二是由紋影系統(tǒng)無法觀測到過小的羽流邊緣溫度梯度引起，最終可能導(dǎo)致高達(dá)中心線速度10%的偏差誤差。隨機(jī)誤差來自時間間隔的不確定性和互相關(guān)算法，局部流速的平均隨機(jī)誤差為7%～8%。根據(jù)偏差誤差和隨機(jī)誤差得到羽流中心線上速度估計的總不確定度為8%，羽流邊緣附近最大測量不確定度為中心線速度的10%～12%。

光流法基于物體在視網(wǎng)膜上的成像原理提出，能夠?qū)崿F(xiàn)運(yùn)動目標(biāo)的檢測、跟蹤等功能[9]，可以基于單個像素獲得速度矢量，具有較高的空間分辨率和測量精度。目前應(yīng)用最廣泛的光流算法基于Horn 等[10]提出的剛體運(yùn)動假設(shè)發(fā)展而來。光流算法假設(shè)運(yùn)動物體在連續(xù)兩幀圖像上的亮度保持恒定，從而得到一個包含亮度與速度u、v的約束條件，即為物理約束條件。Horn 等[10]提出了鄰域內(nèi)像素點(diǎn)速度不能相差太大的假設(shè)，由此得出空間平滑約束條件，從而將速度求解問題轉(zhuǎn)化為一個包含了物理約束和空間平滑約束的最優(yōu)化問題，可以通過變分法對能量方程進(jìn)行最小化得到速度場。由于光流法的約束條件是針對剛體性質(zhì)提出的，Martínez-González 等[11]將傳統(tǒng)光流算法應(yīng)用于自然對流流場紋影圖像，計算結(jié)果僅能顯示主流運(yùn)動趨勢，準(zhǔn)確度和精度均較低。

為得到更適用于流體性質(zhì)的測速方法，許多學(xué)者對光流法的約束條件進(jìn)行了改進(jìn)。Suter[12]引入了二階散度–旋度空間平滑約束形式，可根據(jù)運(yùn)動特征調(diào)節(jié)散度和旋度項(xiàng)的相對比重，能較好地保留流體特性。Liu 等[13]針對紋影/陰影成像以及PIV 等流場可視化手段，推導(dǎo)了通用形式的投影運(yùn)動方程，構(gòu)建了具有流體特性的物理約束條件，誤差分析表明，光流法的不確定度與圖像亮度梯度成反比，與兩幀連續(xù)圖像之間的時間間隔成正比。Arnaud 等[14]結(jié)合紋影亮度方程與連續(xù)性方程推導(dǎo)了新的亮度守恒約束，采用二階形式的空間平滑約束，提出了符合流體運(yùn)動特點(diǎn)的紋影特性光流測速算法。Wang 等[15]在紋影特性光流測速算法中加入金字塔算法和魯棒優(yōu)化等算法，優(yōu)化后的算法能夠得到連續(xù)性較好的速度場，并能捕捉更多精細(xì)的流動結(jié)構(gòu)。梅笑寒等[16]分析了紋影特性光流測速算法因權(quán)重系數(shù)取值不當(dāng)而產(chǎn)生的兩種錯誤，并給出了權(quán)重系數(shù)的合理取值范圍。

上述研究大都針對氣體流場，氣體折射率與密度的關(guān)系滿足Gladstone-Dale 定律[3]。目前針對液體流場的紋影測速方法還鮮見于文獻(xiàn)。與氣體不同，液體折射率與密度的關(guān)系更為復(fù)雜，由Lorentz-Lorenz公式聯(lián)系，一般通過測量給出[17]。在液體中，水是流體研究中應(yīng)用最為廣泛的工質(zhì)，關(guān)于水的密度的研究已經(jīng)比較充分。因此，本文從水的折射率–密度關(guān)系出發(fā)，基于光流法推導(dǎo)適用于水的紋影亮度方程，構(gòu)建新的物理約束條件，采用Wang 等[15]的優(yōu)化算法求解速度場，對浮力羽流紋影圖像進(jìn)行分析，并與互相關(guān)算法、光流算法的計算結(jié)果進(jìn)行對比。

1 算法原理

1.1 互相關(guān)算法

對于時間間隔Δt的圖像序列，當(dāng)Δt適當(dāng)小時，第二幀圖像相對于第一幀圖像有微小位移和變形，此時兩幀圖像存在相關(guān)性。用互相關(guān)算法估計圖像位移場，除以時間間隔Δt，即可得到速度場估計，這種測速方法稱為“圖像相關(guān)測速法”[8]，其算法原理如圖1所示。

圖1 圖像相關(guān)測速法原理[8]Fig.1 Schematic of image-correlation velocimetry algorithm[8]

互相關(guān)算法是在兩幀圖像中劃分若干窗口，通過計算窗口間的互相關(guān)函數(shù)（Cross-correlation function）最小值來得到窗口區(qū)域的最佳匹配關(guān)系，從而估計圖像位移場?；ハ嚓P(guān)函數(shù)的表達(dá)式如式（1）所示。式中，M和N表示在第一幀圖像中選取的窗口區(qū)域大小為M×N個像素；（i,j）表示區(qū)域內(nèi)某點(diǎn)坐標(biāo)，g（i,j）表示該點(diǎn)的亮度值，m和n分別表示該點(diǎn)可能的水平位移和垂直位移；D（m,n）取最小值時，對應(yīng)的m和n即為窗口區(qū)域最佳的位移估計。將整個圖像劃分為若干窗口，即可通過式（1）求得全場位移估計。

1.2 光流算法

光流算法假設(shè)圖像亮度值I（x,y,t）是關(guān)于坐標(biāo)（x,y）和時間t的連續(xù)函數(shù)，某一物體在t時刻和（t+dt）時刻的亮度相等：

將等式右邊進(jìn)行泰勒展開，消去I（x,y,t），再除以dt，即可得到亮度守恒約束方程：

Horn 等[10]提出了空間平滑約束假設(shè)：相鄰點(diǎn)具有相近的速度，u和v應(yīng)滿足方程：

從而將速度場求解轉(zhuǎn)化為一個最優(yōu)化問題，構(gòu)建能量泛函：

式中，α為正則項(xiàng)的權(quán)重，用于平衡正則項(xiàng)在能量泛函中的影響。將式（5）代入變分法中的歐拉–拉格朗日方程，得到兩個等式：

將等式離散化，再進(jìn)行迭代求解，即可得到全場速度u和v。

1.3 液體紋影特性光流測速算法

在擴(kuò)展光源紋影系統(tǒng)中，從光源發(fā)出的光束通過透鏡后，形成平行光穿過測試區(qū)域，再由第二個透鏡匯聚形成光源像，光源像經(jīng)刀口遮擋部分光線后，在像平面上形成紋影圖像。當(dāng)測試區(qū)域的折射率均勻時，像平面的背景亮度為IK，其值與刀口的位置有關(guān)：

式中，h為刀口平面垂直于刀口方向上光源像的寬度，a為未被刀口遮擋的光源像寬度，I0為無刀口遮擋時的背景亮度，如圖2所示。a越小，像平面的背景亮度越暗，紋影圖像就越清晰、靈敏度越高。通常設(shè)置a/h= 0.5，此時刀口平分光源像。

圖2 擴(kuò)展光源紋影系統(tǒng)示意圖Fig.2 Diagram of schlieren system with an extended light source

測試區(qū)域中某點(diǎn)的密度變化會改變該點(diǎn)附近的折射率，經(jīng)過該點(diǎn)的光線將發(fā)生偏折（偏折角為εy），光源像在垂直于刀口的方向上產(chǎn)生偏移Δa，被刀口遮擋的光線數(shù)量改變，使得像平面上該點(diǎn)對應(yīng)位置的亮度發(fā)生變化[7]。亮度變化與折射率滿足以下關(guān)系：

式中，ΔI表示折射率變化引起的像平面上的亮度變化，Z為測試區(qū)域距第二個透鏡的距離，n為測試區(qū)域點(diǎn)（x,y,z）的折射率，n0為周圍介質(zhì)的折射率，z為刀口平面的法向，即光軸方向。

式（9）為Weiss 等[18]給出的水的密度ρ與波長λ、折射率n的關(guān)系式，相關(guān)具體參數(shù)如表1所示。在100 kPa 下、水溫0 ～100 ℃時，水的密度如表2所示[19]。

表1 與水的密度相關(guān)的參數(shù)Table 1 Coefficient related to the density of water

表2 水溫0～100 ℃時水的密度Table 2 Density of water in the range 0-100 ℃

根據(jù)表2，當(dāng)波長在可見光范圍時，通過對數(shù)據(jù)點(diǎn)的擬合，可將式（9）近似化簡為：

由此可以得到亮度變化與密度的關(guān)系：

當(dāng)測試區(qū)域?yàn)橐粋€薄層時，光軸方向上的密度變化可以忽略，上式簡化為：

式中，下標(biāo)0 和1 分別代表第一幀圖像和第二幀圖像，數(shù)據(jù)守恒項(xiàng)為：

空間平滑約束建立了測試區(qū)域內(nèi)相鄰點(diǎn)的速度聯(lián)系，光流算法中的式（4）可對低散度和低渦量速度場進(jìn)行估計，但不適用于流體速度場。Suter 等[12]提出的二階散度–旋度平滑約束被證實(shí)能夠較好地保存流體特性，表達(dá)式如下：

式中，F(xiàn)2為正則項(xiàng)，S1和S2為引入的兩個標(biāo)量場，作為散度、旋度的計算值以簡化能量泛函。β為標(biāo)量場的權(quán)重，β越大，速度場計算值越平滑；β越小，平滑約束越弱，可能造成速度場的劇烈變化。

根據(jù)數(shù)據(jù)守恒約束和空間平滑約束構(gòu)建能量泛函：

式中，正則項(xiàng)權(quán)重α用于平衡正則項(xiàng)F2在能量泛函中的影響。α取值過小時，正則項(xiàng)在能量泛函中比重較小，數(shù)據(jù)項(xiàng)易受噪聲干擾，較弱的平滑約束可能引起計算發(fā)散。與傳統(tǒng)光流法類似，應(yīng)用變分法對式（16）進(jìn)行最小化即可得到速度場，詳見文獻(xiàn)[20]。

2 實(shí)驗(yàn)與計算

2.1 實(shí)驗(yàn)裝置

浮力羽流發(fā)生裝置參考了Gono 等[21]的PIV實(shí)驗(yàn)臺布置，主要由透明水槽、粒子混合腔、流量計、水泵和恒溫水槽等組成，如圖3所示。透明水槽盛裝冷水，底部開口（邊長2 cm）通過上升管道連接粒子混合腔，其后依次與流量計、水泵和恒溫水槽串聯(lián)。

圖3 浮力羽流發(fā)生裝置示意圖Fig.3 Diagram of buoyant plume generator

拍攝紋影圖像時，混合腔內(nèi)無需添加粒子，保持恒溫水槽內(nèi)熱水溫度恒定，通過水泵和流量計將熱水以一定流量泵入粒子混合腔，待腔內(nèi)熱水溫度穩(wěn)定后，打開上升管道，熱水從透明水槽底部開口緩慢流入，在浮力作用下形成羽流緩慢上升。實(shí)驗(yàn)工況為：混合腔內(nèi)熱水溫度30 ℃，透明水槽中冷水溫度15 ℃，熱水流量保持為192 mL/min，透明水槽底部開口處熱水平均速度為8 mm/s。

紋影光路采用了“Z”形雙反射鏡紋影系統(tǒng)，主要由點(diǎn)光源、聚光透鏡、（刀口）狹縫、拋物面鏡、刀口和高速相機(jī)組成，布置方案如圖4所示。

圖4 “Z”形紋影裝置示意圖Fig.4 Diagram of Z-type Schlieren arrangement

在紋影系統(tǒng)中，拋物面鏡之間的距離約為其焦距的2 倍，測試區(qū)域布置于平行光路中央。點(diǎn)光源發(fā)出的光線被聚光透鏡匯聚后通過狹縫，經(jīng)拋物面鏡反射形成平行光穿過測試區(qū)域；受流場影響，平行光發(fā)生偏折，經(jīng)拋物面鏡匯聚后被刀口遮擋，最終在高速相機(jī)成像平面上呈現(xiàn)紋影圖像。

以鹵素?zé)糇鳛閷?shí)驗(yàn)光源，光源功率300 W。兩個拋物面鏡直徑0.3 m、焦距3 m。相機(jī)為Photron FASTCAM SA-Z，圖像尺寸為1024 pixel × 1024 pixel，采樣頻率500 幀/s。標(biāo)定點(diǎn)距離100.0 mm（以千分尺測量），圖像中兩點(diǎn)距離802 pixel，由此計算得到圖像空間分辨率為0.1247 mm/pixel。

2.2 計算方法

采集實(shí)驗(yàn)的紋影圖像數(shù)據(jù)后，分別用互相關(guān)算法、光流算法和液體紋影特性光流測速算法（后文簡稱“液體紋影測速算法”）對圖像進(jìn)行計算。

互相關(guān)算法的計算由PIVlab 完成。PIVlab 是一個基于互相關(guān)算法的MATLAB 工具包，可用于數(shù)字粒子圖像測速（DPIV）分析，也可用于其他環(huán)境（如細(xì)胞內(nèi)流動、沙礫變形等）的分析[22]。PIVlab 的分析過程如下：首先對圖像進(jìn)行前處理，采用自適應(yīng)直方圖均衡化對圖像進(jìn)行增強(qiáng)，可將檢測有效向量的效率提高4.7% ± 3.2%[23]；在計算環(huán)節(jié)，依次將窗口邊長設(shè)置為64、32 和16 pixel，進(jìn)行3 輪離散傅里葉變換，每輪結(jié)果在隨后輪次中用于補(bǔ)償窗口區(qū)域，降低由粒子位移引起的信息損失，同時提升分辨率[24]；最后手動設(shè)置速度閾值以過濾異常值。

光流算法采用Sun 等[25]的優(yōu)化光流算法，綜合平均端點(diǎn)誤差為0.319 pixel/幀[20]。光流算法參數(shù)設(shè)置為：時間分辨率50 幀/s，空間分辨率0.1247 mm/pixel。

液體紋影測速算法的參數(shù)設(shè)置為：時間分辨率50 幀/s，空間分辨率0.1247 mm/pixel。多組計算結(jié)果表明，權(quán)重系數(shù)α= 300、β= 100 時能得到收斂的計算結(jié)果。

3 結(jié)果與分析

3.1 實(shí)驗(yàn)結(jié)果

高速相機(jī)在10 s 內(nèi)記錄了5 000 幀紋影圖像，通過觀察可知羽流呈現(xiàn)周期性變化，周期為1 s。分別進(jìn)行瞬時速度場和平均速度場分析：分析瞬時速度場時，選取的兩幀紋影圖像間隔20 ms（如圖5所示），分別以互相關(guān)算法、光流算法和液體紋影測速算法對兩幀圖像進(jìn)行計算，導(dǎo)出速度矢量圖、速度云圖和渦量圖；分析平均速度場時，選取100 對間隔為20 ms 的圖像序列進(jìn)行計算，對垂直方向的速度場v進(jìn)行平均，導(dǎo)出平均速度矢量和云圖。

圖5 原始紋影圖像Fig.5 Original schlieren image

3.2 瞬時速度場分析

互相關(guān)算法的計算結(jié)果如圖6（a）所示。從速度矢量圖可以看出，在浮力作用下，熱流體在冷流體中向上運(yùn)動，形成湍流結(jié)構(gòu)。從速度云圖可以看出，羽流速度由下至上逐漸增大，出現(xiàn)最大值后又逐漸減小，速度最大值約為0.04 m/s。為了清晰對比3 種方法的差異，將速度云圖局部放大，如圖7所示。從圖7（a）可以看出，互相關(guān)算法的空間分辨率不高，邊界層呈鋸齒狀，出現(xiàn)了一些速度與主流差距較大的塊狀區(qū)域。從渦量圖則可以看出，局部出現(xiàn)了一些較大的渦量。

圖6 3 種方法計算的流場速度結(jié)果Fig.6 Velocity results of three methods

光流算法計算結(jié)果如圖6（b）所示。從速度矢量圖可以看出，速度矢量基本位于主流范圍內(nèi)，在邊緣處速度急劇減小，主流外部區(qū)域的速度接近于零。從圖7（b）可以看出，速度等值區(qū)域呈塊狀分布，即在主流內(nèi)部速度分布相對統(tǒng)一，而在邊緣處變化明顯，表現(xiàn)出較大的速度梯度。從渦量圖可以看出，渦量沿流動結(jié)構(gòu)邊緣呈明顯的線狀分布，在邊緣處有較大的渦量值，而在主流內(nèi)部渦量較小，在高速區(qū)和低速區(qū)之間沒有平緩的過渡。

液體紋影測速算法的計算結(jié)果如圖6（c）所示。從速度矢量圖可以看出，速度矢量在邊緣處逐漸減小，即便是微弱的流動也能顯示出速度矢量。在圖7（c）中，沒有出現(xiàn)過大的速度梯度，體現(xiàn)了流體連續(xù)變化的特性，而且流動結(jié)構(gòu)的邊界較厚，符合邊界層的特點(diǎn)。從渦量圖可以看出，渦量結(jié)果均勻分布于左右兩側(cè)邊界，體現(xiàn)了流體的剪切和旋轉(zhuǎn)特性。

圖7 3 種方法計算流場速度的局部展示Fig.7 Velocity contour of the selected part using three methods

提取垂直速度分量v沿圖6 白色水平虛線的分布，如圖8所示。通過比較速度分布曲線發(fā)現(xiàn)：3 種方法得出的速度值在整體趨勢上較為接近，但在局部位置存在一定差異。互相關(guān)算法和光流算法均在20～30 mm 范圍內(nèi)表現(xiàn)出較大的速度梯度，在30～50 mm范圍內(nèi)出現(xiàn)了一些離群值；而液體紋影測速算法得到的速度場連續(xù)性更好，沒有出現(xiàn)離群值。

圖8 垂直速度分量沿直線的分布Fig.8 Distribution of vertical velocity component along a line

對比3 種算法，由于窗口尺寸限制，互相關(guān)算法的精度和空間分辨率受到影響，當(dāng)追蹤對象超過窗口范圍時易產(chǎn)生錯誤矢量；光流算法的空間平滑約束條件會使速度相近的相鄰點(diǎn)的速度越來越相近，最后整體保持較為統(tǒng)一的速度，如同剛體一般，沒有考慮流體特性；液體紋影測速算法采用二階散度–旋度正則化導(dǎo)出的空間平滑約束條件，結(jié)合紋影亮度方程和流體連續(xù)性方程導(dǎo)出了適用于液體的物理約束條件，可以辨別流體內(nèi)部的細(xì)微差異。通過對比可以看出，采用了具有物理含義特性方程的液體紋影測速算法具有更好的連續(xù)性，能捕捉更多流場細(xì)節(jié)，沒有出現(xiàn)明顯錯誤，整體性能優(yōu)于另外兩種算法。

3.3 平均速度場分析

采用3 種算法對100 對時間間隔為20 ms 的圖像序列進(jìn)行計算。測試流場的速度主要為垂直方向，因此導(dǎo)出y方向速度分量v的平均速度云圖和矢量圖進(jìn)行對比，如圖9所示。圖中v0為熱水在透明水槽底部開口處的平均速度，d為底部開口邊長。

圖9 垂直速度分量的平均速度云圖和矢量圖Fig.9 Mean velocity contours and vectors of vertical velocity component

從圖9 可以看出：垂直方向上，垂直速度分量v隨著與底部開口距離y的增大而增大，在y/d= 2 時達(dá)到最大值，之后隨著y的增大而減?。凰椒较蛏?，羽流的寬度隨著y的增大而增大，體現(xiàn)了垂直于中心線的動量擴(kuò)散。

互相關(guān)算法得到的平均場如圖9（a）所示，白色虛線框內(nèi)部區(qū)域顯示計算出現(xiàn)了離群值，此外，在y/d< 1 的區(qū)域，速度值與其他兩種方法相比偏低；光流算法得到的平均場如圖9（b）所示，羽流邊緣區(qū)域速度梯度仍然較大，過渡不平緩，主流區(qū)域速度略高于其他兩種方法；液體紋影測速算法得到的平均場如圖9（c）所示，速度場近似于軸對稱分布，沒有出現(xiàn)離群值，羽流邊緣區(qū)域速度梯度較光流算法更小，速度緩慢減小。上述對比結(jié)果表明，液體紋影測速算法的計算穩(wěn)定性以及對流體運(yùn)動特性的體現(xiàn)優(yōu)于其他兩種方法。

4 結(jié)論

本文提出了一種基于高速紋影系統(tǒng)的水流流場測速算法。應(yīng)用該算法計算了浮力羽流速度場，與互相關(guān)算法、光流算法的計算結(jié)果進(jìn)行了對比。該算法結(jié)合紋影亮度方程和流體連續(xù)性方程導(dǎo)出了適用于水流的物理約束條件，采用二階散度–旋度正則化作為空間平滑約束條件，計算結(jié)果體現(xiàn)了流體內(nèi)部的細(xì)微差異，得到了邊界層速度。與光流算法相比，該算法能捕捉到微弱的流動結(jié)構(gòu)，保留更多細(xì)節(jié)，避免了過大的速度梯度；與互相關(guān)算法相比，該算法精度和空間分辨率更高，更加穩(wěn)定，不易產(chǎn)生錯誤矢量。