吳 凡, 周 騖, 蔡小舒

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

0 引 言

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

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

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

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

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

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

1 單幀單曝光圖像法測(cè)速原理

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

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

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

(1)

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

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

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

圖2 骨架提取示意圖

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

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

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

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

圖4 新舊算法求軌跡長(zhǎng)度的對(duì)比

Fig.4Thedifferenceontrajectorylengthsbetweenthetraditionalandthenewalgorithm

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

圖5 新舊算法求速度方向的對(duì)比

Fig.5Thedifferenceonvelocitydirectionbetweenthetraditionalandthenewalgorithm

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

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

Fig.6Accuracycontrastbetweenthetraditionalandthenewalgorithm

2.2 速度二義性的判斷

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

2.2.1 角速度判別法

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

圖7 角速度判據(jù)

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

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

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

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

圖8 2幀圖片的匹配與疊加過(guò)程

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

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

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

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

2.3.1 二值化過(guò)程誤差分析

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

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

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

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

圖10 不同二值化過(guò)程的誤差

2.3.2 骨架提取過(guò)程誤差分析

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

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

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

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

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

Fig.12Motionlengthofbinaryandskeletonimageondifferenttrajectorysizes

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

Fig.13Motionlengthofbinaryandskeletonimagewithdifferentwidth/lengthrate(φ)

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

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

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

從隨機(jī)的數(shù)據(jù)擬合出整個(gè)平面場(chǎng)的信息屬于數(shù)學(xué)上插值的范疇。選用徑向基函數(shù)(Radial Basis Function,RBF)插值的方法得出流場(chǎng)的速度場(chǎng)。該插值方法的基本假設(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)

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

(4)

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

(5)

3.2 RBF插值的誤差分析

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

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

(6)

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

圖14 仿真流場(chǎng)(a)與插值流場(chǎng)(b)對(duì)比

Fig.14Contrastbetweensimulationflowfield(a)andinterpolationflowfield(b)

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

4 處理結(jié)果

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

圖15 原始圖像

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

圖16 軌跡骨架圖

圖17 速度云圖

圖18 速度矢量圖

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

(7)

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

圖19 速度的散度云圖

5 結(jié) 論

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

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

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