吳長松,白旭峰,張洋,周立飛

(1.西安交通大學能源與動力工程學院,710049,西安;2.西安陜鼓動力股份有限公司,710075,西安)

隨著圖像識別技術(shù)的提高和精細測量的需要,粒子追蹤測速(PTV)算法在近年來得到了不斷的進步[1]。在需要觀測粒子運動軌跡的情況下[2],PTV則直接跟蹤流場中的粒子運動,從而避免了PIV方法的平均效應[3]。在眾多PTV算法中,有一類算法基于粒子叢思想,即以目標粒子與其周圍粒子的相對位置關(guān)系為判斷基準來實現(xiàn)前后幀粒子匹配。粒子叢算法預設(shè)參數(shù)少,因此運行速度較快,而在構(gòu)建了合理的粒子叢特征后,此類算法在粒子丟失或圖像噪音等問題程度較輕時能夠獲得可觀的準確率。典型的粒子叢類算法有Okamoto等提出的彈簧模型算法[4],以及Ishikawa等提出的速度梯度張量算法[5],這兩種算法的問題在于粒子叢的大小需要提前規(guī)定。因此,當臨近粒子所組成粒子叢的大小隨流動發(fā)生改變時,算法的準確率驟降,而采用適當?shù)木W(wǎng)格劃分技術(shù)則可以避免此問題。Song等提出的基于德勞內(nèi)劃分的算法(DT-PTV)將圖像中離散粒子分布轉(zhuǎn)化成為三角形非結(jié)構(gòu)化網(wǎng)格,通過匹配前后兩幀的三角形單元來實現(xiàn)粒子的匹配[6],張洋等通過引入邏輯判斷模塊實現(xiàn)了DT-PTV運算結(jié)果中錯誤矢量的有效剔除[7],賈攀等將這種方法推廣到了三維情況,即利用四面體單元的信息進行粒子匹配[8]。然而,上述算法存在兩個問題:首先,由于三角形或四面體結(jié)構(gòu)簡單,在粒子濃度較高時會出現(xiàn)形狀趨同,從而導致匹配難度增加;其次,無論對于二維還是三維問題,利用德勞內(nèi)劃分進行粒子匹配時,需要對三角形或四面體進行表層匹配,然后再通過三角形或四面體的頂點配對來實現(xiàn)里層的粒子匹配,這不僅造成算法結(jié)構(gòu)的復雜化,而且在粒子分布相對均勻時將增加里層匹配的邏輯甄別難度。本文以沃羅諾伊劃分為基礎(chǔ),提出了更為簡潔高效的粒子叢類PTV算法,嘗試構(gòu)建以原始粒子為中心的匹配單元,從而僅通過一層匹配完成運算,同時引入人工模擬流場和實際復雜流場來檢驗算法的準確性和實用性。

1 算法準備

1.1 確定候選粒子

粒子叢類算法利用目標粒子與其周圍粒子的相對位置關(guān)系構(gòu)建匹配單元。第1幀中的每個目標粒子需要從第2幀挑選出自己的候選粒子,原則如下

|Xn-Ym|

(1)

式中:Xn、Ym分別為第1、2幀粒子的矢量坐標;Rs為查問域半徑,其取值需要大于兩幀間粒子最大位移,而后者可以通過實際流場的主流速度進行估算。構(gòu)建與粒子一一對應的匹配單元,然后將目標粒子的匹配單元與候選粒子的匹配單元逐一比較,最匹配單元所對應的候選粒子則被定義為匹配粒子。

1.2 兩種對偶劃分

對于二維平面中的離散點,沃羅諾伊劃分是將平面區(qū)域離散點群落轉(zhuǎn)化為互不重疊的多邊形,而德勞內(nèi)劃分是將平面區(qū)域離散點連接成為成互不重疊的三角形,如圖1所示。在構(gòu)建沃羅諾伊網(wǎng)格時,常先將離散點構(gòu)成德勞內(nèi)三角網(wǎng)格,再將德勞內(nèi)三角形各邊的中垂線段進行連接可得到沃羅諾伊劃分,因此沃羅諾伊網(wǎng)格和德勞內(nèi)網(wǎng)格互為對偶,而原始粒子是沃羅諾伊多邊形的中心[9]。在三維空間中,沃羅諾伊劃分(多面體)和三維德勞內(nèi)劃分(四面體)同樣互為對偶。兩種結(jié)構(gòu)在運用時各有好處,可根據(jù)實際情況選擇合適的結(jié)構(gòu),德勞內(nèi)網(wǎng)格曾被用于改進松弛迭代類PTV算法中用來圈選匹配粒子的模塊[10]。

(a)二維平面粒子點 (b)沃羅諾伊多邊形圖1 粒子群的沃羅諾伊劃分

2 算法原理

2.1 二維算法

圖1中原始粒子用“+”符號表示,經(jīng)劃分之后得到的封閉區(qū)域被稱為沃羅諾伊多邊形,它們緊密排列且絕對沒有交疊,并與原始粒子呈一一對應的關(guān)系,即每一個沃羅諾伊多邊形都擁有一個核。對參與匹配的前后兩幀粒子群分別作沃羅諾伊劃分,基于兩幀流動相關(guān)性假設(shè),同一粒子在前后兩幀中對應的沃羅諾伊多邊形變形有限。利用這一特性,將沃羅諾伊多邊形定義為基本匹配單元。

(1)提取多邊形即匹配單元的特征。二維算法中的匹配單元如圖2所示,以原始粒子(x0,y0)為中心建立笛卡爾坐標系,它所對應的匹配單元由若干個三角形Ⅰ,Ⅱ,…,Ⅴ組成,三角形的極半徑r與極角α的關(guān)系為

(2)

最后,此多邊形就被轉(zhuǎn)換成特征曲線,結(jié)果如圖3所示。

圖2 二維算法中的匹配單元

圖3 沃羅諾伊多邊形特征曲線

(3)

式中:Γ為協(xié)方差運算;D為方差運算?？紤]到多邊形跨幀過程中有可能進行旋轉(zhuǎn),需要比較兩條曲線在0°～360°不同相位差下的相似性,本文選擇的相位變化步長為2°。根據(jù)式(3),相比較的兩條曲線長度相同,因此當相位差不為零時,僅取兩條曲線相位重疊的部分進行比較,任意一對特征曲線會計算出多個相關(guān)系數(shù),選擇最大值作為最終相關(guān)系數(shù)。最后,所有候選者中對應最大相關(guān)系數(shù)的粒子被認定為匹配粒子。

考慮到候選粒子中真實匹配粒子已經(jīng)丟失的情況,例如粒子在第2幀中逃逸至圖像之外,只有足夠相似的匹配才會被認為有效,若某一粒子與其所有候選粒子的相關(guān)系數(shù)均小于0.5,則認為此粒子在第2幀中沒有匹配。由于實際流場變形方式及強度均不相同,0.5實際上對應均勻概率,即兩個多邊形完全相同與完全不同兩種極端情況之間的中間點。對于變形較強的流動,可適當放寬匹配條件即相關(guān)性閾值小于0.5,并基于匹配結(jié)果的人眼判斷返回修正這一閾值,直至達到滿意結(jié)果,反之亦然。

(3)遍歷第1幀所有粒子,完成整場粒子匹配,基于匹配關(guān)系和時間間隔可得向前差分格式的瞬時速度場,即實現(xiàn)追蹤測速,而二維算法僅有Rs這一個獨立預設(shè)參數(shù)。

2.2 三維算法

(a)空間中的粒子n (b)粒子n的參考四面體圖4 三維算法中的匹配單元

對三維空間離散的粒子點,采用三維德勞內(nèi)劃分將其轉(zhuǎn)化為空間四面體群,此時名義上的匹配單元為所有包含同一粒子n的四面體所聚合而成的多面體,或者此多面體可以被劃分成若干個以粒子n為共同頂點的四面體,如圖4中所示。但是,實際操作中的基本匹配單元則是這些四面體,如此拆分的原因在于多面體的參數(shù)化遠比二維多邊形的要困難。此處引入一種全新的四面體投票機制,所涉及的多面體匹配與三維DT-PTV中的表層匹配有根本不同,這些作為基本匹配單元的四面體被稱為粒子的參考四面體。如果一個四面體單元可在空間中做任意旋轉(zhuǎn)且具有一個基點,那么它一共有6個自由度。因此,選取6個線性無關(guān)量作為參考四面體的特征量,包括從粒子n指向其余3個頂點的向量的長度,以及代表這3個向量間夾角的點積。將上述特征量匯總為如下兩個四面體特征向量

(4)

參考四面體之間的差異dij,計算式為

dij=‖lni-lmj‖1+‖ani-amj‖1

(5)

式中:下標n和m分別為第1和第2幀中的粒子;下標i和j分別為第1幀和第2幀中粒子的參考四面體;‖‖1為向量的1范數(shù)。dim為dij的最小值,取得最小值的參考四面體為j′,對應的四面體特征向量為lmj′和amj′。對dim進行無量綱化處理,即

(‖lni‖1+‖ani‖1+‖lmj′‖1+‖amj′‖1)-1

(6)

為了避免多面體局部大變形或者粒子缺失所帶來的影響,引入?yún)⒖妓拿骟w投票機制來確定匹配粒子。粒子n的全部參考四面體即為選民,粒子n的候選粒子即為候選人。選民根據(jù)如下公式來決定投票與否,獲得最多認可的候選人當選,即

(7)

vm/Ni>β

(8)

式中:β∈[0,1]為投票系數(shù)。如果式(8)成立,對粒子n的所有候選粒子進行相同操作,否則認為本次投票結(jié)果無效,此時vm為0。vm值最大的候選粒子即為粒子n在第2幀中的匹配粒子,若vm的最大值為0,則認為粒子n在第2幀粒子中沒有匹配結(jié)果。

α代表兩個參考四面體間的差異,在算法中決定是否能投票,β代表兩個粒子之間相似性的預判,在算法中決定投票最終是否有效。如果α增大,根據(jù)式(7)可知,投票的可能性增大,即總票數(shù)增加,此時β也相應增大,從而提高了投票的有效性,使算法在邏輯上更加合理,反之亦然。由式(7)和vm的定義可知,α增加一個較小量就可能會導致總票數(shù)有較大增長,即β有較大增長,β的增速應大于α的。本文通過人工流場測試證明α∈(0,0.25]和β=2α是較為有效的選擇,三維算法僅有Rs和α兩個獨立的預設(shè)參數(shù)。

3 人工流場測試

(a)人工流場速度矢量圖

(b)兩種算法的準確率隨查問域半徑的變化圖5 人工流場及測試結(jié)果

以二維算法為例進行人工流場測試。通過疊加流動解析解生成復雜人工流場。首先在第1幀隨機布置粒子,然后賦予一定的流動規(guī)則,從而在第2幀形成新的粒子位置,并使得前后兩幀聯(lián)合描述一個瞬時流動。流動規(guī)則分多步實施:首先形成具有一定速度梯度的剪切流,然后設(shè)置若干旋渦和偶極子,最后增大流場速度和粒子濃度,并添加相對誤差為5%的位置擾動,以及隨機抹除第2幀數(shù)量比例為5%的粒子。流場大小為256×256像素,旋渦區(qū)大小為128×128像素,偶極子的距離為80像素,兩幀之間最大位移為10像素,由于人工流場真實匹配結(jié)果已知,任一種PTV的匹配結(jié)果與其比較之后即可得到計算準確率,人工流場及測試結(jié)果如圖5所示。隨著Rs的增大,新二維算法的準確率更加穩(wěn)定,這是因為所有的三角形在拓撲結(jié)構(gòu)上是相同的,Rs值的增大將包納更多形狀相似、但實際上毫無關(guān)系的候選匹配單元,這將對DT-PTV的表層匹配產(chǎn)生較大干擾,而多邊形的個體差異遠大于三角形之間的差異,所以能夠更好抵抗Rs的變化所帶來的影響。因此,采用一個足夠大的固定Rs值可使算法滿足一個圖像序列中所有瞬時或者一個瞬時中所有粒子的匹配需要,這兩種情況分別對應流動的非定常性和非均勻性。

4 實驗測試

4.1 二維算法測試

采用脈沖流檢驗算法,流動通過模擬腹主動脈瘤中血液流動實現(xiàn),流動循環(huán)系統(tǒng)如圖6所示。

圖6 人工腹主動脈瘤中脈沖流動循環(huán)系統(tǒng)

系統(tǒng)主要包括位于測試段的腹主動脈瘤模型、一個離心泵以及兩個電磁閥。其中離心泵用來給流體提供動力,兩個電磁閥用來模擬周期約為1 s的人體血流脈沖波形,電磁流量計用來測試主流速度。腹主動脈瘤模型由透明聚甲基丙烯酸甲酯制成,人造血液為45%蒸餾水加55%甘油混合而成,并添加了一定濃度的熒光示蹤粒子。通過高速激光器和狹縫裝置在測試段的流向中心縱剖面上形成片光,然后采用高速相機以最高100 Hz的頻率在側(cè)面記錄粒子圖像。實驗共記錄20個脈沖周期的圖像,選擇其中一個周期測試二維算法,每個周期包含20對圖像,每對圖像對應一個瞬時流場。記人造血液的流動方向為x,鉛垂方向為z。用動態(tài)閾值算法[11]識別圖像中的示蹤粒子并將其轉(zhuǎn)化為單像素點。模型總長L=80 mm,選擇固定查問域半徑為L/5,這已超過主流速度所對應的兩幀間最大位移,選取2個時刻進行分析,結(jié)果與文獻[12]的分析結(jié)果一致,二維算法對人工脈沖流動的重構(gòu)如圖7所示。新算法可以在選取較大Rs值的條件下應對流動序列的非定常性以及瞬時流場的非均勻性。

(a)模擬脈沖流動波形

(b)時刻1的流場重構(gòu)

(c)時刻2的流場重構(gòu)圖7 二維算法對人工脈沖流動的重構(gòu)

4.2 三維算法測試

本實驗在循環(huán)水洞中進行。在水中播撒直徑為10 μm的熒光粒子作為示蹤粒子,在測試段上游布置電磁流量計,測得主流速度約為0.16 m/s。尾跡流實驗段示意圖如圖8所示,圖中透明測試段長約200 cm,截面尺寸為30×30 cm2。V3V系統(tǒng)[13]安裝在測試段的頂部,通過3個鏡頭重構(gòu)示蹤粒子的瞬時三維坐標,拍攝頻率為15 Hz。V3V的測量區(qū)域為一長方體,其厚度為7 cm,截面積為20 cm×20 cm,距離水洞底部10 cm。拍攝區(qū)域上游10 cm處放置一個高為30 cm、截面積為6 cm×6 cm的長方體障礙物以生成尾跡流。

圖8 尾跡流實驗段示意圖

(a)流場全局圖

(b)流場俯視圖

將V3V系統(tǒng)計算得到的粒子三維坐標導入新發(fā)展的三維算法中,預設(shè)參數(shù)α=β/2=0.05,綜合主流速度和采集頻率信息設(shè)置查問域半徑為1.2 cm,通過三維算法重構(gòu)得到的瞬時流場如圖9所示,圖中x方向為主流方向。由圖9可知:算法成功重構(gòu)了速度梯度較大的區(qū)域,即勢流-尾跡流結(jié)合區(qū)域以及尾跡內(nèi)部的旋渦結(jié)構(gòu);圖9c所示的少量錯誤匹配可通過基于邏輯判斷的后處理方法移除[7,14-15]。目前,V3V系統(tǒng)的后處理模塊采用基于松弛迭代的粒子追蹤算法進行粒子匹配,其預設(shè)參數(shù)較多,用戶友好性一般[11],當前的三維流場測試結(jié)果表明,預設(shè)參數(shù)少、結(jié)構(gòu)簡潔的新三維算法亦可用來與V3V等三維成像硬件系統(tǒng)搭配使用。

5 結(jié) 論

本文提出了一種基于沃羅諾伊劃分的粒子叢類PTV算法,并針對匹配單元的構(gòu)建原則提出了相應的二維及三維算法。二維算法充分利用了多邊形形狀的特異性,是一種結(jié)構(gòu)簡單、穩(wěn)定性與準確率兼得的算法;三維算法則采用沃羅諾伊的對偶劃分來構(gòu)建最小匹配單元,從而規(guī)避了三維多面體自由度過多、難以參數(shù)化的問題。人工脈沖流場和尾跡流流場的測試結(jié)果表明,新發(fā)展的算法能很好處理速度梯度較大的復雜流場,具有良好的實用價值。