利用沃羅諾伊劃分的粒子追蹤測(cè)速算法及其在大速度梯度流場(chǎng)中的應(yīng)用

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

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

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

1 算法準(zhǔn)備

1.1 確定候選粒子

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

|Xn-Ym|

(1)

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

1.2 兩種對(duì)偶劃分

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

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

2 算法原理

2.1 二維算法

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

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

(2)

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

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

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

(3)

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

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

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

2.2 三維算法

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

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

(4)

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

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

(5)

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

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

(6)

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

(7)

vm/Ni>β

(8)

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

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

3 人工流場(chǎng)測(cè)試

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

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

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

4 實(shí)驗(yàn)測(cè)試

4.1 二維算法測(cè)試

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

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

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

(a)模擬脈沖流動(dòng)波形

(b)時(shí)刻1的流場(chǎng)重構(gòu)

(c)時(shí)刻2的流場(chǎng)重構(gòu)圖7 二維算法對(duì)人工脈沖流動(dòng)的重構(gòu)

4.2 三維算法測(cè)試

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

圖8 尾跡流實(shí)驗(yàn)段示意圖

(a)流場(chǎng)全局圖

(b)流場(chǎng)俯視圖

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

5 結(jié) 論

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