黃天立，王 倩

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

0 引言

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

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

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

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

上述研究大都針對氣體流場，氣體折射率與密度的關(guān)系滿足Gladstone-Dale 定律[3]。目前針對液體流場的紋影測速方法還鮮見于文獻。與氣體不同，液體折射率與密度的關(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)化算法求解速度場，對浮力羽流紋影圖像進行分析，并與互相關(guān)算法、光流算法的計算結(jié)果進行對比。

1 算法原理

1.1 互相關(guān)算法

對于時間間隔Δt的圖像序列，當Δt適當小時，第二幀圖像相對于第一幀圖像有微小位移和變形，此時兩幀圖像存在相關(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)系，從而估計圖像位移場。互相關(guān)函數(shù)的表達式如式（1）所示。式中，M和N表示在第一幀圖像中選取的窗口區(qū)域大小為M×N個像素；（i,j）表示區(qū)域內(nèi)某點坐標，g（i,j）表示該點的亮度值，m和n分別表示該點可能的水平位移和垂直位移；D（m,n）取最小值時，對應(yīng)的m和n即為窗口區(qū)域最佳的位移估計。將整個圖像劃分為若干窗口，即可通過式（1）求得全場位移估計。

1.2 光流算法

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

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

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

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

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

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

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

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

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

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

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

式中，ΔI表示折射率變化引起的像平面上的亮度變化，Z為測試區(qū)域距第二個透鏡的距離，n為測試區(qū)域點（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，當波長在可見光范圍時，通過對數(shù)據(jù)點的擬合，可將式（9）近似化簡為：

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

當測試區(qū)域為一個薄層時，光軸方向上的密度變化可以忽略，上式簡化為：

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

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

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

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

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

2 實驗與計算

2.1 實驗裝置

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

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

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

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

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

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

以鹵素燈作為實驗光源，光源功率300 W。兩個拋物面鏡直徑0.3 m、焦距3 m。相機為Photron FASTCAM SA-Z，圖像尺寸為1024 pixel × 1024 pixel，采樣頻率500 幀/s。標定點距離100.0 mm（以千分尺測量），圖像中兩點距離802 pixel，由此計算得到圖像空間分辨率為0.1247 mm/pixel。

2.2 計算方法

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

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

光流算法采用Sun 等[25]的優(yōu)化光流算法，綜合平均端點誤差為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 實驗結(jié)果

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

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

3.2 瞬時速度場分析

互相關(guān)算法的計算結(jié)果如圖6（a）所示。從速度矢量圖可以看出，在浮力作用下，熱流體在冷流體中向上運動，形成湍流結(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)的邊界較厚，符合邊界層的特點。從渦量圖可以看出，渦量結(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)算法的精度和空間分辨率受到影響，當追蹤對象超過窗口范圍時易產(chǎn)生錯誤矢量；光流算法的空間平滑約束條件會使速度相近的相鄰點的速度越來越相近，最后整體保持較為統(tǒng)一的速度，如同剛體一般，沒有考慮流體特性；液體紋影測速算法采用二階散度–旋度正則化導(dǎo)出的空間平滑約束條件，結(jié)合紋影亮度方程和流體連續(xù)性方程導(dǎo)出了適用于液體的物理約束條件，可以辨別流體內(nèi)部的細微差異。通過對比可以看出，采用了具有物理含義特性方程的液體紋影測速算法具有更好的連續(xù)性，能捕捉更多流場細節(jié)，沒有出現(xiàn)明顯錯誤，整體性能優(yōu)于另外兩種算法。

3.3 平均速度場分析

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

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

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

互相關(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)定性以及對流體運動特性的體現(xiàn)優(yōu)于其他兩種方法。

4 結(jié)論

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