吳 凡, 周 騖, 蔡小舒

(上海理工大學(xué) 顆粒與兩相流測量研究所,上海 200093)

0 引 言

流動廣泛存在于化工和能源等行業(yè)，復(fù)雜流動現(xiàn)象中流動參數(shù)的獲取通常十分重要[1- 2]。 目前有很多測量顆粒速度的方法，其中非侵入式的光學(xué)測量按照原理有兩大類，即基于多普勒原理的測量方法和基于TOF(Time of Flight)的測量方法[3]。 前者利用激光照射顆粒產(chǎn)生的多普勒頻移測速，后者則根據(jù)單位時間內(nèi)顆粒移動的距離來計算速度[4]。

典型的基于TOF的測量方法包括PIV(Particle Image Velocimetry)和PTV(Particle Tracking Veloci- metry)。PIV是20世紀(jì)70年代末發(fā)展起來的一種非接觸式流體力學(xué)測速方法，是目前應(yīng)用最為廣泛的全場無干擾測速手段[5]。PTV技術(shù)本質(zhì)上是PIV技術(shù)的衍生，兩者的不同之處在于:PIV獲取的速度來自于觀測區(qū)域內(nèi)多個粒子位移的統(tǒng)計平均值，而PTV則是通過追蹤單個粒子的位移來獲得速度。PTV由于其追蹤顆粒運(yùn)動的特點(diǎn)，所得結(jié)果具有直觀性和準(zhǔn)確性，不過同時也具有粒子濃度不宜過高的局限性[6]。早期PTV技術(shù)由于CCD攝像系統(tǒng)水平的

限制，其時間和空間分辨率均不高，精度也不如PIV[7]。隨著計算機(jī)和CCD相機(jī)等硬件水平的不斷提高，在輸移現(xiàn)象和傳送現(xiàn)象等尤其關(guān)注顆粒的拉格朗日運(yùn)動過程的情況下，PTV的作用逐漸顯現(xiàn)出來，因而也得到了廣泛的應(yīng)用[8]。

根據(jù)圖像記錄的方式不同，PTV可以分為3種類型：單幀單脈沖、單幀多脈沖和多幀單脈沖[9]。近年來PTV領(lǐng)域的研究大多集中在多幀單脈沖形式[10]。多幀單脈沖形式具有更高的精度且能很好地解決速度二義性問題，因而得到廣泛應(yīng)用。不過當(dāng)涉及湍流這種十分復(fù)雜、瞬息萬變的流動時， 通常需要頻率非常高的激光脈沖， 才能實(shí)現(xiàn)對復(fù)雜流動顆粒的精確追蹤， 但同時也會造成成本的急劇增加。

單幀單曝光圖像法，即PTV的單幀單脈沖形式，具有系統(tǒng)簡潔、圖像直觀等優(yōu)點(diǎn)。由于只有一次曝光，因此不需要功率非常大的激光器即可完成對示蹤粒子的追蹤。不過單幀單曝光圖像法的圖像處理過程相對復(fù)雜，不同于多幀單脈沖形式的PTV，單幀單曝光雖然沒有粒子匹配的問題，但由于其處理對象是軌跡而非單個點(diǎn)，復(fù)雜度也有所提升。

近年來對PTV的研究多集中在多幀單脈沖形式，對單幀單曝光圖像法的處理算法研究相對較少。其中，劉春嶸等[11]通過形心主軸的概念來計算運(yùn)動軌跡的速度； 陳晶麗等[12]提出針對顆粒運(yùn)動軌跡采用外接矩形的方法來提取顆粒速度與粒度；Jing等[13]則提出了通過檢測軌跡邊緣8個特征點(diǎn)來計算速度與粒度的方法。上述3種方法各有優(yōu)勢和特點(diǎn)，不過三者有一個共同點(diǎn)，那就是針對一條軌跡處理出一個速度矢量，導(dǎo)致信息的利用率較低。同時，在處理彎曲軌跡時，這些方法的精度均有所下降，而復(fù)雜流場中的彎曲軌跡是無法避免的。因此，本文提出一種方法，不僅能同時提取直線段軌跡和彎曲段軌跡的速度矢量，同時也提高軌跡上運(yùn)動信息的利用率。

1 單幀單曝光圖像法測速原理

在利用單幀單曝光圖像法測量示蹤粒子的速度時，一般CCD相機(jī)固定不動，通過視窗記錄測量區(qū)的圖像，此時，離散相顆粒(如水滴、煤粉、氣泡等)與CCD相機(jī)發(fā)生相對位移， 其成像在曝光時間內(nèi)的疊加就形成了模糊的“拖影”，如圖1所示[14]。

單幀單曝光運(yùn)動模糊圖像包含了顆粒與相機(jī)的相對位移信息。因?yàn)槠毓鈺r間已知，在物體勻速直線運(yùn)動假設(shè)的前提下，可以計算出運(yùn)動速度，如式(1)所示。

圖1 運(yùn)動模糊圖像及其成像模型

(1)

式中，T為曝光時間且已知，故速度v的計算關(guān)鍵就在于運(yùn)動長度L的計算?，F(xiàn)有的很多算法都將軌跡近似視為一段直線，且最終得出一個速度矢量。以外接矩形法為例，該方法求出軌跡的最小外接矩形，提取矩形的長和寬，長度減去寬度便得到了顆粒的運(yùn)動距離。其他幾種算法各有不同，但都具有一定的相似性。上述算法均需滿足顆粒運(yùn)動軌跡接近直線的前提，在顆粒軌跡為曲線的情況下會產(chǎn)生一定誤差。所以需要提出另一種計算速度的方法，能同時適用于直線軌跡和彎曲軌跡的情況。

2 基于骨架提取的顆粒軌跡處理算法

骨架提取，也被稱為二值圖像細(xì)化。這種算法能將一個連通區(qū)域細(xì)化成一個像素的寬度，同時保持圖像的形態(tài)學(xué)特征，主要用于特征提取和目標(biāo)拓?fù)浔硎尽膱D2可以很直觀地看出骨架提取的特點(diǎn)。

圖2 骨架提取示意圖

結(jié)合骨架提取的特點(diǎn)可以將顆粒的軌跡進(jìn)行細(xì)化，得到單像素寬度的軌跡圖，也就是像素精度的跡線。在此基礎(chǔ)上進(jìn)行軌跡長度計算要準(zhǔn)確得多，且對于彎曲軌跡也更能區(qū)分各像素點(diǎn)處速度方向的細(xì)微區(qū)別。為得到效果更好的軌跡骨架圖，需要對原圖依次進(jìn)行去噪、自適應(yīng)閾值分割、閉運(yùn)算、去除小顆粒及骨架提取等操作，主要的處理流程如圖3所示[15]。

2.1 基于軌跡骨架圖的速度提取

在原理介紹部分已經(jīng)提到，速度計算的前提是計算運(yùn)動長度。在獲得軌跡的單像素寬度跡線之后，通過連通區(qū)域分割函數(shù)可以識別出每條軌跡上的所有像素點(diǎn)，傳統(tǒng)算法采用跡線起點(diǎn)與終點(diǎn)之間的距離作為長度，這樣在處理彎曲軌跡時就會產(chǎn)生較大誤差。本文提出類似于求線積分的求長度算法，首先將每2個像素點(diǎn)之間的距離求出來，然后迭加得到總的軌跡長度，2種算法的對比如圖4所示。

圖3 圖像處理主要流程[15]

圖4 新舊算法求軌跡長度的對比

Fig.4Thedifferenceontrajectorylengthsbetweenthetraditionalandthenewalgorithm

對于速度的方向，以往的計算方式是采取軌跡首尾相連的一條直線作為方向向量，如圖4中的舊算法。但這種方法在處理彎曲軌跡的時候會產(chǎn)生很大誤差，且一條軌跡上所有的點(diǎn)共享一個速度矢量，十分籠統(tǒng)。因此，提出了求跡線切線的方法，將像素點(diǎn)前后若干個像素點(diǎn)之間的割線近似視為切線,并作為軌跡上該點(diǎn)的速度方向，不僅對于彎曲軌跡能很好地提高精度，而且同一條軌跡上各個像素點(diǎn)的速度矢量也有區(qū)分，能很好地體現(xiàn)顆粒的實(shí)際運(yùn)動過程。圖5顯示了新舊算法求速度方向的區(qū)別，新算法(左側(cè))各個像素點(diǎn)處的速度由其切線決定，能更好地顯示粒子的流動過程，舊算法(右側(cè))則過于籠統(tǒng)，所有像素點(diǎn)的速度矢量完全相同，不能體現(xiàn)實(shí)際的流動狀態(tài)。此外，若將軌跡上各個點(diǎn)處的速度矢量逐幀播放，也能更清晰生動地顯示顆粒的動態(tài)流動過程。

圖5 新舊算法求速度方向的對比

Fig.5Thedifferenceonvelocitydirectionbetweenthetraditionalandthenewalgorithm

圖6為新舊算法在計算不同彎曲程度軌跡時的精度對比。從仿真圓形軌跡上截取不同弧長的軌跡，然后分別采用2種算法求長，并與弧長的理論值進(jìn)行對比。圖中，橫坐標(biāo)為軌跡的弧長，也代表了軌跡的彎曲程度，縱坐標(biāo)為軌跡長度。紅色圓點(diǎn)為舊算法計算的長度，藍(lán)色三角點(diǎn)為新算法計算的長度，黑色方點(diǎn)為實(shí)際弧長。由圖6可見，在弧長較低，即軌跡接近直線的情況下，2種算法的結(jié)果均十分接近真值；隨著弧長的增加(即軌跡彎曲程度增加)，新算法則逐漸體現(xiàn)出了優(yōu)勢。

圖6 2種速度提取算法精度比較

Fig.6Accuracycontrastbetweenthetraditionalandthenewalgorithm

2.2 速度二義性的判斷

在前文的敘述中，速度的計算忽略了一個很重要的問題，就是速度二義性的判斷。對一條軌跡來說，一般有2個端點(diǎn)，那么軌跡是從端點(diǎn)A運(yùn)動到B還是從B運(yùn)動到A是2種完全不同的情況，所以需要對該問題進(jìn)行判斷。下面提供2種方法：一種是利用流體的角速度進(jìn)行判斷，另一種則利用多幀圖像之間的運(yùn)動關(guān)系進(jìn)行判斷。

2.2.1 角速度判別法

角速度判別法是一種利用右手螺旋定則來確定軌跡運(yùn)動方向的方法。如圖7所示，有一組圍繞同一中心點(diǎn)進(jìn)行旋轉(zhuǎn)的4條軌跡，此時執(zhí)行前述的速度提取流程后，上下2條軌跡的速度方向會保持一致(水平向右)。但實(shí)際上，由于這是一組旋轉(zhuǎn)軌跡，上下2條軌跡的方向應(yīng)該相反。為使同一旋轉(zhuǎn)中心的一組軌跡旋轉(zhuǎn)方向保持一致，使用角速度判別法進(jìn)行判斷。圖7中，通過右手螺旋定則可以知道上方軌跡的角速度垂直紙面朝里，而下方的角速度垂直紙面向外，二者互相矛盾，通過改變其中一條軌跡的速度矢量即可解決速度二義性問題。在實(shí)際執(zhí)行過程中可預(yù)先定義軌跡的旋轉(zhuǎn)方向(順時針或逆時針)，從而統(tǒng)一各個軌跡的速度方向，消除二義性。

圖7 角速度判據(jù)

對于直線運(yùn)動狀態(tài)，可以將旋轉(zhuǎn)中心點(diǎn)選在垂直于軌跡的無窮遠(yuǎn)處。此外，對于復(fù)雜的流動，可以通過選取多個旋轉(zhuǎn)中心分別進(jìn)行角速度判別，由于過程比較復(fù)雜，故方法具有一定的局限性。

2.2.2 多幀連續(xù)曝光圖像匹配

單幀圖像很難解決二義性問題，而多幀圖像則可以。因此，在相機(jī)進(jìn)行第一段曝光之后，立刻進(jìn)行第二段曝光，這樣得到2幀在時間序列上幾乎銜接的圖片，由于2幀圖片的時間先后順序已知，就能很方便地求出顆粒的速度方向。

在多幀圖像匹配方法中，關(guān)鍵在于2幀圖片中同一條軌跡的匹配識別[16]。由于2次曝光之間的時間間隔非常短暫，所以2段軌跡理論上是能首尾相連的。實(shí)際處理中，對于第一幀圖片里的所有軌跡，于2個端點(diǎn)附近搜索其在第二幀圖片中是否有軌跡，搜索半徑由2段曝光之間的間距時間以及流場的速度決定。為了保證在搜索半徑中只有該軌跡的下一幀軌跡出現(xiàn)，因此顆粒密度不宜太大。此外，考慮到邊緣部分的軌跡可能由于進(jìn)出視場而導(dǎo)致無法匹配，故應(yīng)提前將邊緣部分的軌跡剔除。圖8展示了2幀連續(xù)曝光圖片(曝光時間為7ms，間隔時間為0.077ms)匹配并疊加的過程，其軌跡匹配率為84.5%。

圖8 2幀圖片的匹配與疊加過程

在完成2幀圖片之間的軌跡匹配之后，速度二義性就解決了。舉例來說：假設(shè)AB是第一幀圖片的某條軌跡，端點(diǎn)分別是A和B，而CD是該軌跡在第二幀圖片上的位置，端點(diǎn)分別是C和D；若B和C是2條軌跡的連接點(diǎn)，則軌跡的方向就依次是A—B—C—D。

多幀圖像匹配之后，不僅解決了速度二義性問題，同時，由于得到了不同時間序列上的速度值，也能因此獲得加速度。

2.3 二值化與骨架提取的誤差分析

基于骨架提取的軌跡處理算法可以分成2部分：圖像處理部分和速度計算部分。速度計算部分的誤差已經(jīng)在2.1節(jié)進(jìn)行了描述，此處對二值化和骨架提取的誤差進(jìn)行分析。

2.3.1 二值化過程誤差分析

單閾值分割通常采用Ostu (最大類間方差)閾值，其通過圖像灰度將圖像分為背景和目標(biāo)，通過最大化背景和目標(biāo)之間的類間方差來找尋最優(yōu)閾值?？紤]到單閾值分割在灰度分布不均的情況下并不適用，故使用了自適應(yīng)閾值分割。這是一種多閾值分割方法，通過分析像素周圍局部范圍的灰度特性來決定該像素點(diǎn)的閾值，能夠使閾值分割后的圖片保留更多信息。

二值化過程的誤差分析實(shí)驗(yàn)采用了仿真的運(yùn)動模糊圖片，如圖9所示。根據(jù)實(shí)際射流實(shí)驗(yàn)圖片，將仿真顆粒圖像的背景灰度設(shè)為40，目標(biāo)灰度設(shè)為200，粒徑為5pixel，將顆粒圖像進(jìn)行卷積并添加方差為0.008的高斯加性噪聲即得到仿真運(yùn)動模糊圖像，卷積核則由運(yùn)動長度決定。

圖9 多種運(yùn)動的仿真運(yùn)動模糊圖片

仿真實(shí)驗(yàn)的結(jié)果如圖10所示，可見隨著運(yùn)動長度的增加，2種閾值分割的誤差均顯著下降，且自適應(yīng)閾值分割的誤差要明顯低于Ostu 閾值分割。

圖10 不同二值化過程的誤差

2.3.2 骨架提取過程誤差分析

骨架提取過程將二值圖轉(zhuǎn)變成骨架圖，故誤差分析過程采用仿真的軌跡二值圖(由一系列圓形疊加形成的類似于橢圓的圖形)進(jìn)行一系列仿真實(shí)驗(yàn)。初步認(rèn)為骨架提取的結(jié)果受軌跡大小、軌跡長短軸比、軌跡方向(位置)三方面的影響。實(shí)驗(yàn)過程中分別計算二值圖與骨架圖中軌跡對應(yīng)的運(yùn)動長度。在二值圖中，運(yùn)動長度為軌跡的長度減去寬度；在骨架圖中，軌跡長度即為像素距離迭代求和的長度。仿真實(shí)驗(yàn)結(jié)果如圖11～13所示。

從圖11可以看出骨架提取結(jié)果與軌跡的方向無關(guān)；從圖12和13看出，骨架提取操作對運(yùn)動長度的影響很小。進(jìn)一步考察原始數(shù)據(jù)發(fā)現(xiàn)，在所有情況下，骨架圖中的運(yùn)動長度都比二值圖的運(yùn)動長度大一

個固定值1，這是由于長度計量的起點(diǎn)不同所導(dǎo)致的。舉例來說，假設(shè)有一個僅包含一個像素點(diǎn)的軌跡，其長度應(yīng)該為1，但長度減去寬度卻等于0，從而導(dǎo)致了二者相差1的現(xiàn)象，這種情況在實(shí)際應(yīng)用過程中統(tǒng)一即可。

圖11 骨架提取結(jié)果與軌跡方向的關(guān)系

圖12 不同軌跡尺寸下二值圖與骨架圖的運(yùn)動長度

Fig.12Motionlengthofbinaryandskeletonimageondifferenttrajectorysizes

圖13 不同長寬比(φ)下二值圖與骨架圖的運(yùn)動長度

Fig.13Motionlengthofbinaryandskeletonimagewithdifferentwidth/lengthrate(φ)

3 流場重構(gòu)算法

在執(zhí)行了上述軌跡處理算法之后，可以得到每條軌跡的速度信息。然而多數(shù)情況下實(shí)際需要的是整個流場的速度信息，因此需要進(jìn)行流場重構(gòu)工作。

3.1 RBF插值法重構(gòu)流體速度場

從隨機(jī)的數(shù)據(jù)擬合出整個平面場的信息屬于數(shù)學(xué)上插值的范疇。選用徑向基函數(shù)(Radial Basis Function,RBF)插值的方法得出流場的速度場。該插值方法的基本假設(shè)是插值點(diǎn)的函數(shù)值受插值源函數(shù)值影響，且影響程度與插值點(diǎn)和插值源之間的距離有關(guān)，其數(shù)學(xué)表達(dá)式如下：

(2)

根據(jù)不同應(yīng)用實(shí)例，可以選取不同的距離系數(shù)函數(shù)Φ以達(dá)到最佳的插值效果。比較多種函數(shù)的插值效果，從中選取名為“multiquadric” 的距離系數(shù)，其距離系數(shù)函數(shù)如下：

(3)

考慮到流體速度場的特點(diǎn)，任意一點(diǎn)的速度受來流方向速度值的影響大于下游方向速度值的影響。需要在RBF插值的基礎(chǔ)上進(jìn)行類似于“迎風(fēng)插值”的優(yōu)化。優(yōu)化的基本思想是在插值點(diǎn)的上游增加對插值點(diǎn)的影響；在插值點(diǎn)的下游減弱對插值點(diǎn)的影響。因此，可以構(gòu)建迎風(fēng)修正系數(shù)函數(shù)ψ，ψ的值由插值點(diǎn)與已知點(diǎn)之間速度向量的夾角決定，當(dāng)夾角為銳角時，ψ大于0，對結(jié)果進(jìn)行增強(qiáng)；當(dāng)夾角為鈍角時，ψ小于0，對結(jié)果削弱。其表達(dá)式如下:

(4)

修正后的插值表達(dá)式如下:

(5)

3.2 RBF插值的誤差分析

對插值的誤差進(jìn)行分析，采用仿真輸出的流場結(jié)果，隨機(jī)選取一定的數(shù)據(jù)點(diǎn)之后進(jìn)行插值，分析插值結(jié)果與原流場的差別即可。原流場與插值流場的速度云圖如圖14所示。

為了量化2幅圖之間的相似度，以便更直觀地獲知不同插值方法的優(yōu)劣，定義2幅圖像之間的平均距離δ，如式(6)所示。δ越大，表示兩者之間的差別越大，δ越小，表示兩者越為接近。

(6)

表1顯示了選用不同插值函數(shù)并采用“迎風(fēng)插值”優(yōu)化之后的插值結(jié)果的平均距離。從中可以看出基于 “multiquadric” 并采用“迎風(fēng)插值”優(yōu)化之后得到插值效果是最好的，與肉眼直觀相符。

圖14 仿真流場(a)與插值流場(b)對比

Fig.14Contrastbetweensimulationflowfield(a)andinterpolationflowfield(b)

表1 插值結(jié)果誤差分析Table 1 Error analysis on interpolation

4 處理結(jié)果

為驗(yàn)證本文算法在實(shí)際情況下的適用性，采用射流實(shí)驗(yàn)中拍攝的圖片作為處理對象，如圖15所示。實(shí)驗(yàn)中射流出口流速為0.15m/s，拍攝部位為射流口上部的卷吸區(qū)域。射流口位于圖片右下方，視場大小為2.88mm×3.86mm。圖16展示了軌跡骨架圖，圖17和18分別是最終速度場的云圖和矢量圖。

圖15 原始圖像

在速度云圖中，每個像素點(diǎn)處的顏色代表了該處的速度大小，因而能直觀地看出速度場在空間上的能量分布；在速度矢量圖中，箭頭方向?yàn)樗俣确较颍^長度表征了速度大小，從而可以看出流場整體運(yùn)動趨勢。

圖16 軌跡骨架圖

圖17 速度云圖

圖18 速度矢量圖

為驗(yàn)證實(shí)驗(yàn)結(jié)果的可信度，利用插值后的速度場來計算速度的散度場。考慮到此處速度場并非連續(xù)，因此采用變分方式來計算散度，同時假定流場為二維流動，散度計算公式如下：

(7)

采用云圖方式展示計算的散度場，如圖19所示，由于常溫常壓下水流可視為不可壓流動，速度的散度應(yīng)處處為零，圖19中大多數(shù)散度值接近為零，少數(shù)幾處出現(xiàn)稍大的散度值(正值或負(fù)值)，不過總體上結(jié)果符合事實(shí)。

圖19 速度的散度云圖

5 結(jié) 論

為獲取顆粒的速度分布，針對單幀單曝光圖像的處理提出一種顆粒軌跡處理算法。通過骨架提取獲得軌跡的跡線，然后對跡線的像素坐標(biāo)進(jìn)行迭代求和來計算軌跡長度，同時通過求跡線的切線來計算速度方向。相比傳統(tǒng)算法，該算法在速度大小和方向上均有更高精度，且不同于傳統(tǒng)算法一條軌跡一個速度矢量的模式，該算法的一條軌跡可以獲得軌跡上任意點(diǎn)處的速度矢量，大大提高了信息利用率。

在得到原始的離散速度場之后，利用迎風(fēng)優(yōu)化后的RBF 插值算法，可以重構(gòu)一個相對連續(xù)且平滑的全分布流場，可在此基礎(chǔ)上計算加速度、散度等參數(shù)，方便后續(xù)研究；也可以很方便地與CFD 等仿真模擬結(jié)果進(jìn)行比較。

本算法處理結(jié)果的精度上限取決于閾值分割過程，因而需要在拍攝過程中調(diào)整曝光時間、光源、顆粒濃度等因素，盡可能提高所拍攝圖像中軌跡的清晰度，以期達(dá)到好的閾值分割效果，減少圖片中的干擾軌跡，從而獲得更高精度的結(jié)果。