基于Delaunay網(wǎng)格技術(shù)的松弛迭代粒子追蹤算法

賈 攀，王 元，張 洋

（西安交通大學(xué)流體機(jī)械及工程系，陜西 西安 710049）

0 引言

作為一種非接觸式的測量技術(shù)，PTV（Particle Tracking Velocimetry）已廣泛用于多相流［1］、生物流體力學(xué)［2］等研究領(lǐng)域，成為實(shí)驗(yàn)流體力學(xué)的有力工具。PTV 算法著眼于離散粒子及其拉格朗日描述［3］，近年來取得了多樣化的發(fā)展。應(yīng)用PTV 算法分析高粒子濃度的PIV 流場圖像一直是流場測量的一個研究方向［4］。目前提出的算法中，ORX（Original Relaxation Method）算法為此提供了一個很好的選擇，它克服了一些常見算法的缺陷，可實(shí)現(xiàn)對剪切等復(fù)雜流場高濃度、高效率的測量［4］。Barnard 和Thompson［5］首先將ORX 用于標(biāo)記實(shí)景測速，此后Lee和Back［3］率先應(yīng)用ORX 分析湍流剪切流場，證實(shí)了ORX 對復(fù)雜流場處理的高效性；Ohmi 和Huang［4］通過對ORX 無匹配概率的迭代方式進(jìn)行修正，提高了計(jì)算粒子團(tuán)無匹配的準(zhǔn)確率，提出了NRX（New Relaxation Method）算法；張洋等人［6］詳細(xì)分析了NRX 中“孤粒子無匹配問題”的來源以及結(jié)果的偽邏輯性，提出了有效改進(jìn)該問題的雙向法則。應(yīng)用NRX 為第1幀中的粒子在第2 幀中尋找匹配粒子時，首先需要給定兩個固定閾值，分別用以確定臨近目標(biāo)粒子并與之保持相似運(yùn)動形式的參選粒子集合和與目標(biāo)粒子預(yù)匹配的候選粒子集合，這兩個閾值的選取依賴于流場的參數(shù)，同時也和流場示蹤粒子的分布情況有關(guān)，選取合理與否將直接影響匹配計(jì)算的準(zhǔn)確率及效率。本文嘗試應(yīng)用DT（Delaunay Triangulation）方法確定候選和參選粒子集合來擺脫固定閾值以提高算法的獨(dú)立性，并通過模擬二維應(yīng)變流動中的單極渦運(yùn)動進(jìn)行了全面驗(yàn)證。

1 結(jié)合雙向法則的NRX的基本原理

NRX 的原理如圖1所示。記錄連續(xù)運(yùn)動的兩幀圖像分別記為X、Y，其粒子形心坐標(biāo)分別表示為坐標(biāo)矢量x、y。第1幀中第i個粒子的坐標(biāo)矢量為xi，第2幀中第j個粒子的坐標(biāo)矢量為yj；目標(biāo)是為第1幀中任一粒子xi在第2幀中尋找匹配粒子，步驟如下：

第一步，確定目標(biāo)粒子xi的參選和候選粒子集合。設(shè)候選粒子及參選粒子矢量的個數(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的無匹配概率為，由概率的歸一化準(zhǔn)則：

通過對計(jì)算效率和加權(quán)預(yù)判進(jìn)行折衷考慮，為了簡單起見，迭代初值和按如下格式確定：

如圖1（b）所示迭代初始值，箭頭表示匹配，數(shù)字為概率大小，n.m.表示“無匹配”。

第三步，對xi的所有Pij以及迭代運(yùn)算：

A＝0.3，B＝3，分別為松弛迭代的衰減和增益系數(shù)［4］。Θ 為參選位移矢量的集合：

其中dij＝xi－yj，dkl＝xk－xl，E、F均為常數(shù)。其形式是全局量E（由流場具體情況決定）與局部量（F＝0.05）之和，從而能更好地適應(yīng)速度梯度較大的流場［4］，圖1（a）中Rc＝E＋F｜dij｜。

其中，Qk和Zk表示xk的參選位移矢量個數(shù)和其所有參選粒子的候選粒子總數(shù)，C、D為經(jīng)驗(yàn)常數(shù)。迭代得到和之后，歸一化：

至此，xi的所有匹配概率一次迭代完成。然后重復(fù)第三步，直至各匹配概率收斂。xi的匹配對象為代表最大概率的yi（或無匹配）。對第1幀圖像X中的每個粒子進(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é)果。

2 應(yīng)用DT 方法確定參選和候選粒子集合

如前文所述，應(yīng)用NRX 時，需先確定兩個閾值Rs和Rn，分別用以確定參選和候選參選粒子集合。Rs可取值為兩幀時間間隔內(nèi)粒子可能的最大位移，一般來說，通過對所測流場的基本分析和半定量的估計(jì)，可以得到比較合理的Rs，然而Rn并沒有明確的確定方法，一般都是通過考慮流場粒子示蹤的平均距離以及算法的計(jì)算效率經(jīng)驗(yàn)地給定一個估計(jì)值［3－4］，這種經(jīng)驗(yàn)性的方法實(shí)際中不便應(yīng)用，選取的準(zhǔn)確與否對算法的效率等也會有明顯的影響。為了減少NRX對具體流場參數(shù)的依賴從而使之通用化，本文嘗試引入DT 方法確定參選和候選粒子集合。

DT 方法主要用于處理離散點(diǎn)的問題，目前廣泛用于有限元分析和圖像的插值問題［7］。DT 方法是將離散點(diǎn)有效連接為平面三角形或空間三棱錐，形成非結(jié)構(gòu)化網(wǎng)格［7］。目前DT 網(wǎng)格的生成方法很多，主要有Lawson算法和Bowyer－Watson算法［8］。此外，DT 方法生成的網(wǎng)格具有很好的“空圓特性”，即DT網(wǎng)格中任一三角形的外接圓范圍內(nèi)不會有其它節(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)集時需要判斷是否與第2幀中的粒子重合，若重合，則要作相應(yīng)的處理保證候選粒子集合中包括該重合粒子。

3 自定義圖像

通過模擬二維應(yīng)變流動中的單極渦流運(yùn)動［9］驗(yàn)證DT 方法的有效性。根據(jù)運(yùn)動方程，生成匹配關(guān)系明確且粒子分布均勻的兩幀自定義流場圖像［10］。選取典型流動形式分別為剪切流動、旋轉(zhuǎn)流動以及雙曲流動。設(shè)圖像大小為a，圖像粒子總數(shù)為Np，定義容積率Pv為圖像像素總數(shù)與粒子數(shù)目之比，即Pv＝a2／Np，顯然a一定的情況下，Pv越小，粒子數(shù)Np越大。為保證粒子的數(shù)目，應(yīng)用無粘連粒子生成法生成第1幀圖像；然后根據(jù)給定的運(yùn)動規(guī)律，生成第2幀圖像，同時，為模擬圖像噪音，對第二幀圖像進(jìn)行隨機(jī)“添加”和“抹除”處理，定義“添加”和“抹除”的粒子數(shù)與第1幀粒子總數(shù)的比值分別為添加率μa和抹除率μd。在容積率Pv遞變和引入干擾的前提下，分別應(yīng)用DT 方法和傳統(tǒng)固定閾值方法確定參選和候選粒子集合進(jìn)行匹配計(jì)算，并與真實(shí)情況比較從而判斷優(yōu)劣。圖2所示為三種自定義流場圖像，其中圖像邊長a＝256像素、pv＝21.8524。

圖2 自定義流場圖像Fig.2 Self－defintion flow field images

4 計(jì)算結(jié)果及分析

4.1 DT方法確定參選粒子集合的有效性驗(yàn)證

按照前文所述方法，3種流場分別采用10種容積率，同時取添加率μa＝15%、抹除率μd＝10%，生成30組圖像對。然后，分別應(yīng)用固定閾值Rn和DT兩種方法確定參選粒子集合，應(yīng)用式（1）確定候選粒子集合，對自定義流場圖像進(jìn)行雙向匹配計(jì)算［6］，相關(guān)參數(shù)見表1。通過對3種流場共630 種工況進(jìn)行匹配計(jì)算，結(jié)果表明3種流場的計(jì)算結(jié)果都定性的一致，故此處僅針對旋轉(zhuǎn)流場的匹配結(jié)果進(jìn)行討論。

表1 圖像及算法相關(guān)參數(shù)（Ⅰ）Table 1 Parameters for images and algorithm（Ⅰ）

定義正確匹配矢量總數(shù)和第1幀粒子總數(shù)的比值為匹配準(zhǔn)確率η。對所有工況的計(jì)算準(zhǔn)確率進(jìn)行分析，結(jié)果表明：應(yīng)用固定閾值確定參選粒子集合時，在Rn的變化范圍內(nèi)，準(zhǔn)確率η的變化趨勢定性一致，并且都可保持在96%以上；應(yīng)用DT 方法確定參選粒子集合時，準(zhǔn)確率η也可達(dá)到96%以上，也表現(xiàn)出和固定閾值情況下一致的變化趨勢。圖3所示為4種不同的Rn以及DT 方法確定參選粒子集合時，雙向匹配計(jì)算［6］的準(zhǔn)確率η隨容積率Pv的變化趨勢，圖示可知：5種情況下，準(zhǔn)確率η都保持在較高水平，并且變化趨勢都隨容積率的降低有所減小。因此，與合理選取的固定閾值Rn一樣，DT 方法也可在較大的容積率Pv范圍內(nèi)有效地確定參選粒子集合，計(jì)算結(jié)果保持高的匹配準(zhǔn)確率。

圖3 準(zhǔn)確率隨容積率的變化趨勢Fig.3 Tendency of matching accuracy

假設(shè)相同的條件下，N種工況的計(jì)算時間為ti，tmax＝max｛ti｜i＝1～N｝，通過對時間無量綱化定義計(jì)算效率ηT＝1－ti／tmax，顯然，最大的計(jì)算時間對應(yīng)的計(jì)算效率為0。對所有工況的計(jì)算時間進(jìn)行分析，結(jié)果表明：一方面，應(yīng)用固定閾值選取參選粒子集合時，各個固定閾值Rn下，ηT隨Pv的減小而減小，這是因?yàn)殡S著Pv的減小，粒子數(shù)目增加，數(shù)據(jù)的處理量增大；對于固定的容積率Pv，計(jì)算效率ηT隨著固定閾值Rn的增大有所降低，這是由于隨著Rn的增大圈選的參選粒子的數(shù)量增加，數(shù)據(jù)處理量增加，并且這種現(xiàn)象隨著Pv取值的減小逐漸顯著；然而當(dāng)固定閾值Rn取值過小（Rn＜5像素）時，計(jì)算效率會明顯降低，這是因?yàn)楣潭ㄩ撝颠^小時，不能有效地確定參選粒子集合，雙向計(jì)算［6］概率擇優(yōu)模塊的耗時過大。另一方面，應(yīng)用DT 方法確定參選粒子集合時，計(jì)算效率ηT也隨容積率Pv的變化趨勢與固定閾值情況下一致，并且對于給定的容積率Pv，應(yīng)用DT 方法確定參選粒子集合時計(jì)算效率達(dá)到最高。

圖4所示為4種不同的Rn以及DT 方法確定參選粒子集合時，雙向匹配計(jì)算［6］的計(jì)算效率ηT隨容積率Pv的變化趨勢，圖示可知：5種情況下，計(jì)算效率ηT都隨容積率Pv的減小而減?。籖n＝2像素時，由于取值過小，計(jì)算效率最低；DT 方法確定參選粒子集合時，計(jì)算效率最高。

如圖5所示為應(yīng)用DT 方法和Rn＝9像素確定參選粒子集合時計(jì)算效率ηT隨Pv的變化趨勢，Δη表示DT 方法與固定閾值Rn相比計(jì)算效率的提高。圖示結(jié)果表明：計(jì)算效率的提高Δη隨容積率Pv的降低明顯上升。

圖4 計(jì)算效率的變化趨勢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)格的“空圓特性”保證了與目標(biāo)粒子xi直接相連形成三角形的粒子也是幾何上臨近的粒子，這種“就近性”選取原則和NRX 的假設(shè)是一致的。因此，應(yīng)用DT 方法和固定閾值確定參選粒子集合的物理本質(zhì)是一致的，DT 方法可以有效地確定參選粒子集合，并且DT 方法可以自適應(yīng)容積率Pv的變化確定合理的參選粒子集合，顯示出對高粒子濃度粒子圖像分析的高效性。

4.2 DT方法確定候選粒子集合的有效性驗(yàn)證

首先，按照前文所述的方法，每種流場采用10種容積率，每種容積率采用8種粒子幀間最大位移，取添加率μa＝15%、抹除率μd＝10%，生成3種自定義流場共240對圖像。然后，分別應(yīng)用式（1）和DT 方法確定候選粒子集合，應(yīng)用式（2）確定參選粒子集合，對自定義流場圖像進(jìn)行雙向匹配計(jì)算［6］，相關(guān)參數(shù)見表2。通過對3種流場共480組工況進(jìn)行匹配計(jì)算，結(jié)果表明，3種流場的計(jì)算結(jié)果都定性一致，故此處僅針對旋轉(zhuǎn)流場的匹配結(jié)果進(jìn)行討論。

表2 圖像和算法的相關(guān)參數(shù)（Ⅱ）Table 2 Parameters for images and algorithm（Ⅱ）

圖6所示為容積率Pv＝16.384，粒子最大位移m＝6像素時，旋轉(zhuǎn)流場的匹配計(jì)算結(jié)果。圖6（a）為應(yīng)用固定閾值方法確定候選粒子集合時的計(jì)算結(jié)果，匹配準(zhǔn)確率η為97.45%，可以看出計(jì)算錯誤矢量基本上是隨機(jī)分布在圖像中，數(shù)量較少；圖6（b）為應(yīng)用DT 方法確定候選粒子集合的情況下的計(jì)算結(jié)果，匹配準(zhǔn)確率 僅為86.26%，從圖6（b）中可以看出計(jì)算錯誤矢量明顯多于圖6（a）的情況，且絕大部分的錯誤匹配矢量在粒子位移較大的位置出現(xiàn)，這是由于真正匹配粒子和目標(biāo)粒子xi之間的距離較大，DT 方法“就近性”確定的候選粒子集合中包含真正匹配粒子的概率減小。

圖6 匹配結(jié)果比較Fig.6 Comparison of the matching results

對計(jì)算結(jié)果分析，結(jié)果表明：一方面，應(yīng)用固定閾值確定候選粒子集合情況下，m在其變化范圍內(nèi)固定取值時η隨Pv的降低有所減??；Pv在變化范圍內(nèi)固定取值時，η隨m的增大而基本保持不變。另一方面，應(yīng)用DT 方法確定候選粒子集合時，η隨Pv的降低而減小，并且減小的幅度大于應(yīng)用固定閾值確定候選粒子的情況；并且η隨m的增大明顯減小。

圖7所示為m＝6像素時，旋轉(zhuǎn)流場匹配準(zhǔn)確率η隨Pv的變化趨勢。ηN和ηD分別為應(yīng)用固定閾值Rs和DT 方法確定候選粒子集合情況下，準(zhǔn)確率隨Pv的變化曲線。圖示可知：ηN在Pv變化范圍內(nèi)始終保持在95%以上，說明只要根據(jù)粒子的最大位移m合理地選取固定閾Rs，就可有效的確定候選粒子集合；然而，ηD僅當(dāng)容積率Pv較高時保持高的值，隨著Pv進(jìn)一步降低，ηD則呈現(xiàn)明顯的下降趨勢，當(dāng)Pv接近10時，ηD低于85%，這是因?yàn)殡S著Pv的持續(xù)減小，粒子濃度增大，目標(biāo)粒子xi附近的非真正匹配粒子增多，DT 網(wǎng)格中目標(biāo)粒子和真正匹配粒子直接連接形成三角形的概率下降。

圖8所示為Pv＝16.384時，旋轉(zhuǎn)流場匹配準(zhǔn)確率η隨m的變化趨勢。ηN和ηD分別為應(yīng)用固定閾值Rs和DT 方法確定候選粒子集合情況下，η隨m的變化曲線。圖示可知：ηN在m變化范圍內(nèi)始終保持在95%以上，說明只要對流場進(jìn)行基本分析，合理估計(jì)粒子最大位移m，給出合理的閾值Rs，就可得到合理的候選粒子集合；而然，ηD則隨m的增大呈現(xiàn)明顯的下降趨勢，這是由于隨著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ù)雜流場高濃度粒子圖像的優(yōu)點(diǎn)相悖；另一方面是減小幀間粒子位移，對于特定的流場，則需要提高粒子圖像的采樣頻率，這對實(shí)驗(yàn)平臺的圖像采集系統(tǒng)提出更加苛刻的要求，同時也會加大數(shù)據(jù)的處理量。

候選粒子集合中應(yīng)該包括目標(biāo)粒子的真實(shí)匹配粒子，實(shí)際應(yīng)用NRX 時，若可對流場進(jìn)行基本分析，給定恰當(dāng)?shù)拈撝担涂傻玫胶侠淼暮蜻x粒子集合。DT 方法時僅從幾何位置出發(fā)確定候選粒子集合，由于DT 網(wǎng)格的“空圓特性”，選取的粒子都在位置上臨近目標(biāo)粒子，然而對于一些復(fù)雜的流場情況，由于間粒子的最大位移m較大，真實(shí)的匹配粒子位于距目標(biāo)粒子較遠(yuǎn)的位置。因此，DT 方法這種“就近性”選取的特點(diǎn)與確定候選粒子的物理背景不統(tǒng)一，不能合理地確定候選粒子集合。

5 結(jié)論

（1）由于DT 網(wǎng)格的“空圓特性”所致的“就近性”選取參選粒子特征的與NRX 中確定參選粒子集合的物理本質(zhì)一致，DT 方法可以有效地確定參選粒子集合。

（2）DT 方法不能有效地確定候選粒子集合，由于DT 網(wǎng)格的“就近性”選取的特點(diǎn)與NRX 確定的候選粒子集合中必須包含真正的候選粒子的物理背景不一致，并且這種矛盾在容積率降低、幀間粒子位移增大的情況下更加突出。

（3）在NRX 的基礎(chǔ)上應(yīng)用DT 方法確定參選粒子集合提出一種新的DT－NRX 算法。DT 方法引入減少了一個固定閾值的確定，提高了算法的獨(dú)立性。DT－NRX 與NRX 算法比較表明：一方面，與NRX 類似，DT－NRX對容積率變化有很好的適應(yīng)性保證了算法準(zhǔn)確率高的特點(diǎn)，另一方面，DT 方法可以自適應(yīng)容積率的變化有效地確定參選粒子集合，表現(xiàn)出對高粒子濃度流場圖像的分析的高效性。

［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é)合雙向法則的松弛迭代粒子追蹤測速法［J］.空氣動力學(xué)學(xué)報，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.