賈 攀,王 元,張 洋
(西安交通大學(xué)流體機(jī)械及工程系,陜西 西安 710049)
作為一種非接觸式的測(cè)量技術(shù),PTV(Particle Tracking Velocimetry)已廣泛用于多相流[1]、生物流體力學(xué)[2]等研究領(lǐng)域,成為實(shí)驗(yàn)流體力學(xué)的有力工具。PTV 算法著眼于離散粒子及其拉格朗日描述[3],近年來(lái)取得了多樣化的發(fā)展。應(yīng)用PTV 算法分析高粒子濃度的PIV 流場(chǎng)圖像一直是流場(chǎng)測(cè)量的一個(gè)研究方向[4]。目前提出的算法中,ORX(Original Relaxation Method)算法為此提供了一個(gè)很好的選擇,它克服了一些常見(jiàn)算法的缺陷,可實(shí)現(xiàn)對(duì)剪切等復(fù)雜流場(chǎng)高濃度、高效率的測(cè)量[4]。Barnard 和Thompson[5]首先將ORX 用于標(biāo)記實(shí)景測(cè)速,此后Lee和Back[3]率先應(yīng)用ORX 分析湍流剪切流場(chǎng),證實(shí)了ORX 對(duì)復(fù)雜流場(chǎng)處理的高效性;Ohmi 和Huang[4]通過(guò)對(duì)ORX 無(wú)匹配概率的迭代方式進(jìn)行修正,提高了計(jì)算粒子團(tuán)無(wú)匹配的準(zhǔn)確率,提出了NRX(New Relaxation Method)算法;張洋等人[6]詳細(xì)分析了NRX 中“孤粒子無(wú)匹配問(wèn)題”的來(lái)源以及結(jié)果的偽邏輯性,提出了有效改進(jìn)該問(wèn)題的雙向法則。應(yīng)用NRX 為第1幀中的粒子在第2 幀中尋找匹配粒子時(shí),首先需要給定兩個(gè)固定閾值,分別用以確定臨近目標(biāo)粒子并與之保持相似運(yùn)動(dòng)形式的參選粒子集合和與目標(biāo)粒子預(yù)匹配的候選粒子集合,這兩個(gè)閾值的選取依賴(lài)于流場(chǎng)的參數(shù),同時(shí)也和流場(chǎng)示蹤粒子的分布情況有關(guān),選取合理與否將直接影響匹配計(jì)算的準(zhǔn)確率及效率。本文嘗試應(yīng)用DT(Delaunay Triangulation)方法確定候選和參選粒子集合來(lái)擺脫固定閾值以提高算法的獨(dú)立性,并通過(guò)模擬二維應(yīng)變流動(dòng)中的單極渦運(yùn)動(dòng)進(jìn)行了全面驗(yàn)證。
NRX 的原理如圖1所示。記錄連續(xù)運(yùn)動(dòng)的兩幀圖像分別記為X、Y,其粒子形心坐標(biāo)分別表示為坐標(biāo)矢量x、y。第1幀中第i個(gè)粒子的坐標(biāo)矢量為xi,第2幀中第j個(gè)粒子的坐標(biāo)矢量為yj;目標(biāo)是為第1幀中任一粒子xi在第2幀中尋找匹配粒子,步驟如下:
第一步,確定目標(biāo)粒子xi的參選和候選粒子集合。設(shè)候選粒子及參選粒子矢量的個(gè)數(shù)分別為Mi、Ni,候選粒子和參選粒子集合分別為Sc、Sr:
如圖1所示,Rs和Rn分別為用以圈選候選粒子矢量和參選粒子矢量的閾值,意義為:xi真正的匹配粒子位于以Rs為半徑的圓內(nèi);以Rn為半徑的圓內(nèi)的粒子真實(shí)位移矢量“擬平行”[3]。
圖1 NRX算法原理Fig.1 Algorithm for NRX
第二步,迭代概率初始化。設(shè)xi和yj匹配的概率為Pij;xi的無(wú)匹配概率為,由概率的歸一化準(zhǔn)則:
通過(guò)對(duì)計(jì)算效率和加權(quán)預(yù)判進(jìn)行折衷考慮,為了簡(jiǎn)單起見(jiàn),迭代初值和按如下格式確定:
如圖1(b)所示迭代初始值,箭頭表示匹配,數(shù)字為概率大小,n.m.表示“無(wú)匹配”。
第三步,對(duì)xi的所有Pij以及迭代運(yùn)算:
A=0.3,B=3,分別為松弛迭代的衰減和增益系數(shù)[4]。Θ 為參選位移矢量的集合:
其中dij=xi-yj,dkl=xk-xl,E、F均為常數(shù)。其形式是全局量E(由流場(chǎng)具體情況決定)與局部量(F=0.05)之和,從而能更好地適應(yīng)速度梯度較大的流場(chǎng)[4],圖1(a)中Rc=E+F|dij|。
其中,Qk和Zk表示xk的參選位移矢量個(gè)數(shù)和其所有參選粒子的候選粒子總數(shù),C、D為經(jīng)驗(yàn)常數(shù)。迭代得到和之后,歸一化:
至此,xi的所有匹配概率一次迭代完成。然后重復(fù)第三步,直至各匹配概率收斂。xi的匹配對(duì)象為代表最大概率的yi(或無(wú)匹配)。對(duì)第1幀圖像X中的每個(gè)粒子進(jìn)行第一至第三步計(jì)算,得到X到Y(jié)的匹配結(jié)果。然后根據(jù)雙向法則[6],按照上述步驟再次計(jì)算得到Y(jié)到X的匹配結(jié)果,合并兩份結(jié)果并統(tǒng)一至X到Y(jié)的模式,最后進(jìn)行概率擇優(yōu),得到最終的匹配結(jié)果。
如前文所述,應(yīng)用NRX 時(shí),需先確定兩個(gè)閾值Rs和Rn,分別用以確定參選和候選參選粒子集合。Rs可取值為兩幀時(shí)間間隔內(nèi)粒子可能的最大位移,一般來(lái)說(shuō),通過(guò)對(duì)所測(cè)流場(chǎng)的基本分析和半定量的估計(jì),可以得到比較合理的Rs,然而Rn并沒(méi)有明確的確定方法,一般都是通過(guò)考慮流場(chǎng)粒子示蹤的平均距離以及算法的計(jì)算效率經(jīng)驗(yàn)地給定一個(gè)估計(jì)值[3-4],這種經(jīng)驗(yàn)性的方法實(shí)際中不便應(yīng)用,選取的準(zhǔn)確與否對(duì)算法的效率等也會(huì)有明顯的影響。為了減少NRX對(duì)具體流場(chǎng)參數(shù)的依賴(lài)從而使之通用化,本文嘗試引入DT 方法確定參選和候選粒子集合。
DT 方法主要用于處理離散點(diǎn)的問(wèn)題,目前廣泛用于有限元分析和圖像的插值問(wèn)題[7]。DT 方法是將離散點(diǎn)有效連接為平面三角形或空間三棱錐,形成非結(jié)構(gòu)化網(wǎng)格[7]。目前DT 網(wǎng)格的生成方法很多,主要有Lawson算法和Bowyer-Watson算法[8]。此外,DT 方法生成的網(wǎng)格具有很好的“空?qǐng)A特性”,即DT網(wǎng)格中任一三角形的外接圓范圍內(nèi)不會(huì)有其它節(jié)點(diǎn)的存在[7]。
設(shè)目標(biāo)粒子為xi,確定參選粒子集合的方法為:以第1幀圖像X中的粒子為節(jié)點(diǎn)生成DT 網(wǎng)格,然后選取網(wǎng)格中所有直接和xi相連形成三角形的粒子組成參選粒子集合。相似地,可確定候選粒子集合:首先將xi置于第2幀圖像Y中形成新的離散點(diǎn)集并生成DT 網(wǎng)格,然后選取與xi直接相連形成三角形的粒子組成候選粒子集合。將xi置于第2幀圖像形成新的點(diǎn)集時(shí)需要判斷是否與第2幀中的粒子重合,若重合,則要作相應(yīng)的處理保證候選粒子集合中包括該重合粒子。
通過(guò)模擬二維應(yīng)變流動(dòng)中的單極渦流運(yùn)動(dòng)[9]驗(yàn)證DT 方法的有效性。根據(jù)運(yùn)動(dòng)方程,生成匹配關(guān)系明確且粒子分布均勻的兩幀自定義流場(chǎng)圖像[10]。選取典型流動(dòng)形式分別為剪切流動(dòng)、旋轉(zhuǎn)流動(dòng)以及雙曲流動(dòng)。設(shè)圖像大小為a,圖像粒子總數(shù)為Np,定義容積率Pv為圖像像素總數(shù)與粒子數(shù)目之比,即Pv=a2/Np,顯然a一定的情況下,Pv越小,粒子數(shù)Np越大。為保證粒子的數(shù)目,應(yīng)用無(wú)粘連粒子生成法生成第1幀圖像;然后根據(jù)給定的運(yùn)動(dòng)規(guī)律,生成第2幀圖像,同時(shí),為模擬圖像噪音,對(duì)第二幀圖像進(jìn)行隨機(jī)“添加”和“抹除”處理,定義“添加”和“抹除”的粒子數(shù)與第1幀粒子總數(shù)的比值分別為添加率μa和抹除率μd。在容積率Pv遞變和引入干擾的前提下,分別應(yīng)用DT 方法和傳統(tǒng)固定閾值方法確定參選和候選粒子集合進(jìn)行匹配計(jì)算,并與真實(shí)情況比較從而判斷優(yōu)劣。圖2所示為三種自定義流場(chǎng)圖像,其中圖像邊長(zhǎng)a=256像素、pv=21.8524。
圖2 自定義流場(chǎng)圖像Fig.2 Self-defintion flow field images
按照前文所述方法,3種流場(chǎng)分別采用10種容積率,同時(shí)取添加率μa=15%、抹除率μd=10%,生成30組圖像對(duì)。然后,分別應(yīng)用固定閾值Rn和DT兩種方法確定參選粒子集合,應(yīng)用式(1)確定候選粒子集合,對(duì)自定義流場(chǎng)圖像進(jìn)行雙向匹配計(jì)算[6],相關(guān)參數(shù)見(jiàn)表1。通過(guò)對(duì)3種流場(chǎng)共630 種工況進(jìn)行匹配計(jì)算,結(jié)果表明3種流場(chǎng)的計(jì)算結(jié)果都定性的一致,故此處僅針對(duì)旋轉(zhuǎn)流場(chǎng)的匹配結(jié)果進(jìn)行討論。
表1 圖像及算法相關(guān)參數(shù)(Ⅰ)Table 1 Parameters for images and algorithm(Ⅰ)
定義正確匹配矢量總數(shù)和第1幀粒子總數(shù)的比值為匹配準(zhǔn)確率η。對(duì)所有工況的計(jì)算準(zhǔn)確率進(jìn)行分析,結(jié)果表明:應(yīng)用固定閾值確定參選粒子集合時(shí),在Rn的變化范圍內(nèi),準(zhǔn)確率η的變化趨勢(shì)定性一致,并且都可保持在96%以上;應(yīng)用DT 方法確定參選粒子集合時(shí),準(zhǔn)確率η也可達(dá)到96%以上,也表現(xiàn)出和固定閾值情況下一致的變化趨勢(shì)。圖3所示為4種不同的Rn以及DT 方法確定參選粒子集合時(shí),雙向匹配計(jì)算[6]的準(zhǔn)確率η隨容積率Pv的變化趨勢(shì),圖示可知:5種情況下,準(zhǔn)確率η都保持在較高水平,并且變化趨勢(shì)都隨容積率的降低有所減小。因此,與合理選取的固定閾值Rn一樣,DT 方法也可在較大的容積率Pv范圍內(nèi)有效地確定參選粒子集合,計(jì)算結(jié)果保持高的匹配準(zhǔn)確率。
圖3 準(zhǔn)確率隨容積率的變化趨勢(shì)Fig.3 Tendency of matching accuracy
假設(shè)相同的條件下,N種工況的計(jì)算時(shí)間為ti,tmax=max{ti|i=1~N},通過(guò)對(duì)時(shí)間無(wú)量綱化定義計(jì)算效率ηT=1-ti/tmax,顯然,最大的計(jì)算時(shí)間對(duì)應(yīng)的計(jì)算效率為0。對(duì)所有工況的計(jì)算時(shí)間進(jìn)行分析,結(jié)果表明:一方面,應(yīng)用固定閾值選取參選粒子集合時(shí),各個(gè)固定閾值Rn下,ηT隨Pv的減小而減小,這是因?yàn)殡S著Pv的減小,粒子數(shù)目增加,數(shù)據(jù)的處理量增大;對(duì)于固定的容積率Pv,計(jì)算效率ηT隨著固定閾值Rn的增大有所降低,這是由于隨著Rn的增大圈選的參選粒子的數(shù)量增加,數(shù)據(jù)處理量增加,并且這種現(xiàn)象隨著Pv取值的減小逐漸顯著;然而當(dāng)固定閾值Rn取值過(guò)小(Rn<5像素)時(shí),計(jì)算效率會(huì)明顯降低,這是因?yàn)楣潭ㄩ撝颠^(guò)小時(shí),不能有效地確定參選粒子集合,雙向計(jì)算[6]概率擇優(yōu)模塊的耗時(shí)過(guò)大。另一方面,應(yīng)用DT 方法確定參選粒子集合時(shí),計(jì)算效率ηT也隨容積率Pv的變化趨勢(shì)與固定閾值情況下一致,并且對(duì)于給定的容積率Pv,應(yīng)用DT 方法確定參選粒子集合時(shí)計(jì)算效率達(dá)到最高。
圖4所示為4種不同的Rn以及DT 方法確定參選粒子集合時(shí),雙向匹配計(jì)算[6]的計(jì)算效率ηT隨容積率Pv的變化趨勢(shì),圖示可知:5種情況下,計(jì)算效率ηT都隨容積率Pv的減小而減?。籖n=2像素時(shí),由于取值過(guò)小,計(jì)算效率最低;DT 方法確定參選粒子集合時(shí),計(jì)算效率最高。
如圖5所示為應(yīng)用DT 方法和Rn=9像素確定參選粒子集合時(shí)計(jì)算效率ηT隨Pv的變化趨勢(shì),Δη表示DT 方法與固定閾值Rn相比計(jì)算效率的提高。圖示結(jié)果表明:計(jì)算效率的提高Δη隨容積率Pv的降低明顯上升。
圖4 計(jì)算效率的變化趨勢(shì)Fig.4 Tendency of matching efficiency
圖5 兩種方法計(jì)算效率比較Fig.5 Matching efficiency comparison
NRX 算法假設(shè)互相臨近的粒子的位移“擬平行”,這也是確定參選粒子集合的準(zhǔn)則。DT 方法以粒子為節(jié)點(diǎn)生成三角形網(wǎng)格,DT 網(wǎng)格的“空?qǐng)A特性”保證了與目標(biāo)粒子xi直接相連形成三角形的粒子也是幾何上臨近的粒子,這種“就近性”選取原則和NRX 的假設(shè)是一致的。因此,應(yīng)用DT 方法和固定閾值確定參選粒子集合的物理本質(zhì)是一致的,DT 方法可以有效地確定參選粒子集合,并且DT 方法可以自適應(yīng)容積率Pv的變化確定合理的參選粒子集合,顯示出對(duì)高粒子濃度粒子圖像分析的高效性。
首先,按照前文所述的方法,每種流場(chǎng)采用10種容積率,每種容積率采用8種粒子幀間最大位移,取添加率μa=15%、抹除率μd=10%,生成3種自定義流場(chǎng)共240對(duì)圖像。然后,分別應(yīng)用式(1)和DT 方法確定候選粒子集合,應(yīng)用式(2)確定參選粒子集合,對(duì)自定義流場(chǎng)圖像進(jìn)行雙向匹配計(jì)算[6],相關(guān)參數(shù)見(jiàn)表2。通過(guò)對(duì)3種流場(chǎng)共480組工況進(jìn)行匹配計(jì)算,結(jié)果表明,3種流場(chǎng)的計(jì)算結(jié)果都定性一致,故此處僅針對(duì)旋轉(zhuǎn)流場(chǎng)的匹配結(jié)果進(jìn)行討論。
表2 圖像和算法的相關(guān)參數(shù)(Ⅱ)Table 2 Parameters for images and algorithm(Ⅱ)
圖6所示為容積率Pv=16.384,粒子最大位移m=6像素時(shí),旋轉(zhuǎn)流場(chǎng)的匹配計(jì)算結(jié)果。圖6(a)為應(yīng)用固定閾值方法確定候選粒子集合時(shí)的計(jì)算結(jié)果,匹配準(zhǔn)確率η為97.45%,可以看出計(jì)算錯(cuò)誤矢量基本上是隨機(jī)分布在圖像中,數(shù)量較少;圖6(b)為應(yīng)用DT 方法確定候選粒子集合的情況下的計(jì)算結(jié)果,匹配準(zhǔn)確率 僅為86.26%,從圖6(b)中可以看出計(jì)算錯(cuò)誤矢量明顯多于圖6(a)的情況,且絕大部分的錯(cuò)誤匹配矢量在粒子位移較大的位置出現(xiàn),這是由于真正匹配粒子和目標(biāo)粒子xi之間的距離較大,DT 方法“就近性”確定的候選粒子集合中包含真正匹配粒子的概率減小。
圖6 匹配結(jié)果比較Fig.6 Comparison of the matching results
對(duì)計(jì)算結(jié)果分析,結(jié)果表明:一方面,應(yīng)用固定閾值確定候選粒子集合情況下,m在其變化范圍內(nèi)固定取值時(shí)η隨Pv的降低有所減?。籔v在變化范圍內(nèi)固定取值時(shí),η隨m的增大而基本保持不變。另一方面,應(yīng)用DT 方法確定候選粒子集合時(shí),η隨Pv的降低而減小,并且減小的幅度大于應(yīng)用固定閾值確定候選粒子的情況;并且η隨m的增大明顯減小。
圖7所示為m=6像素時(shí),旋轉(zhuǎn)流場(chǎng)匹配準(zhǔn)確率η隨Pv的變化趨勢(shì)。ηN和ηD分別為應(yīng)用固定閾值Rs和DT 方法確定候選粒子集合情況下,準(zhǔn)確率隨Pv的變化曲線。圖示可知:ηN在Pv變化范圍內(nèi)始終保持在95%以上,說(shuō)明只要根據(jù)粒子的最大位移m合理地選取固定閾Rs,就可有效的確定候選粒子集合;然而,ηD僅當(dāng)容積率Pv較高時(shí)保持高的值,隨著Pv進(jìn)一步降低,ηD則呈現(xiàn)明顯的下降趨勢(shì),當(dāng)Pv接近10時(shí),ηD低于85%,這是因?yàn)殡S著Pv的持續(xù)減小,粒子濃度增大,目標(biāo)粒子xi附近的非真正匹配粒子增多,DT 網(wǎng)格中目標(biāo)粒子和真正匹配粒子直接連接形成三角形的概率下降。
圖8所示為Pv=16.384時(shí),旋轉(zhuǎn)流場(chǎng)匹配準(zhǔn)確率η隨m的變化趨勢(shì)。ηN和ηD分別為應(yīng)用固定閾值Rs和DT 方法確定候選粒子集合情況下,η隨m的變化曲線。圖示可知:ηN在m變化范圍內(nèi)始終保持在95%以上,說(shuō)明只要對(duì)流場(chǎng)進(jìn)行基本分析,合理估計(jì)粒子最大位移m,給出合理的閾值Rs,就可得到合理的候選粒子集合;而然,ηD則隨m的增大呈現(xiàn)明顯的下降趨勢(shì),這是由于隨著m的增大,DT 網(wǎng)格中真正匹配粒子和目標(biāo)粒子的距離也增大,DT 方法“就近性”的確定候選粒子集合的缺陷也更加突出。
圖7 準(zhǔn)確率隨Pv 的變化Fig.7 Tendency of matching accuracy with Pv
圖8 準(zhǔn)確率隨m 的變化Fig.8 Tendency of matching accuracy with m
綜合以上分析可知:若要提高準(zhǔn)確率,一方面需提高粒子圖像的容積率,但是這一點(diǎn)和松弛迭代算法處理復(fù)雜流場(chǎng)高濃度粒子圖像的優(yōu)點(diǎn)相悖;另一方面是減小幀間粒子位移,對(duì)于特定的流場(chǎng),則需要提高粒子圖像的采樣頻率,這對(duì)實(shí)驗(yàn)平臺(tái)的圖像采集系統(tǒng)提出更加苛刻的要求,同時(shí)也會(huì)加大數(shù)據(jù)的處理量。
候選粒子集合中應(yīng)該包括目標(biāo)粒子的真實(shí)匹配粒子,實(shí)際應(yīng)用NRX 時(shí),若可對(duì)流場(chǎng)進(jìn)行基本分析,給定恰當(dāng)?shù)拈撝担涂傻玫胶侠淼暮蜻x粒子集合。DT 方法時(shí)僅從幾何位置出發(fā)確定候選粒子集合,由于DT 網(wǎng)格的“空?qǐng)A特性”,選取的粒子都在位置上臨近目標(biāo)粒子,然而對(duì)于一些復(fù)雜的流場(chǎng)情況,由于間粒子的最大位移m較大,真實(shí)的匹配粒子位于距目標(biāo)粒子較遠(yuǎn)的位置。因此,DT 方法這種“就近性”選取的特點(diǎn)與確定候選粒子的物理背景不統(tǒng)一,不能合理地確定候選粒子集合。
(1)由于DT 網(wǎng)格的“空?qǐng)A特性”所致的“就近性”選取參選粒子特征的與NRX 中確定參選粒子集合的物理本質(zhì)一致,DT 方法可以有效地確定參選粒子集合。
(2)DT 方法不能有效地確定候選粒子集合,由于DT 網(wǎng)格的“就近性”選取的特點(diǎn)與NRX 確定的候選粒子集合中必須包含真正的候選粒子的物理背景不一致,并且這種矛盾在容積率降低、幀間粒子位移增大的情況下更加突出。
(3)在NRX 的基礎(chǔ)上應(yīng)用DT 方法確定參選粒子集合提出一種新的DT-NRX 算法。DT 方法引入減少了一個(gè)固定閾值的確定,提高了算法的獨(dú)立性。DT-NRX 與NRX 算法比較表明:一方面,與NRX 類(lèi)似,DT-NRX對(duì)容積率變化有很好的適應(yīng)性保證了算法準(zhǔn)確率高的特點(diǎn),另一方面,DT 方法可以自適應(yīng)容積率的變化有效地確定參選粒子集合,表現(xiàn)出對(duì)高粒子濃度流場(chǎng)圖像的分析的高效性。
[1]TSUJI T,MIYAUCHI T,OH S,et al.Simultaneous measurement of particle motion and temperature in twodimensional fluidized bed with heat transfer[J].Kona PowderandParticleJournal,2010,28(28):167-179.
[2]EVANGELOS B,MICHELE G,UFUK O,et al.CFD and PTV steady flow investigation in an anatomically accurate abdominal aortic aneurysm[J].JournalofBiomechanicalEngineering,2009,131(1):1-15.
[3]BEAK S J,LEE S J.A new two-frame particle tracking algorithm using match probability[J].Experimentsin Fluids,1996,22(1):23-32.
[4]OHMI K,LI H Y.Particle-tracking velocimetry with new algorithms[J].MeasurementScienceandTechnology,2000,11(6):603-616.
[5]BARNARD S T,THOMPSON W B.Disparity analysis of images[J].IEEETrans.PatternAnalysisMachine Intelligence,1980,2(4):333-340.
[6]張洋,王元,李志強(qiáng).結(jié)合雙向法則的松弛迭代粒子追蹤測(cè)速法[J].空氣動(dòng)力學(xué)學(xué)報(bào),2010,28(3):250-254.
[7]馬逸塵,梅立泉,王阿霞.偏微分方程現(xiàn)代數(shù)值方法[M].北京:科學(xué)出版社,2006:25-29.
[8]WASTON D F.Computing the n dimensional DelaunayTessellation with application to Voronoi ploytopes[J].ComputerJournal,1981,24(2):167-172.
[9]TRELING R,BECKERS M,Van HEIJST.Dynamics of monopolar vortices in a strain flow[J].Journalof FluidMechanics,1997,345:165-201.
[10]OKAMOTO K,NISHIO S,SAGA T,et al.Standard images for particle-image velocimetry[J].Measurement ScienceandTechnology,2000,11(6):685-691.