基于粒子圖像分割的混合PIV-PTV算法

李 拓，張清福，潘 翀，2，陳 爽，申俊琦，王宏偉，李曉輝，黃 湛，王晉軍

(1. 北京航空航天大學 流體力學教育部重點實驗室，北京 100191；2. 北京航空航天大學寧波創(chuàng)新研究院 先進飛行器與空天動力創(chuàng)新研究中心, 寧波 315800；3. 中國空氣動力研究與發(fā)展中心 設(shè)備設(shè)計與測試技術(shù)研究所, 綿陽 621000；4. 中國航天空氣動力技術(shù)研究院, 北京 100074)

0 引言

粒子圖像測速（particle image velocimetry, PIV）和粒子追蹤測速（particle tracking velocimetry，PTV）等基于數(shù)字圖像處理的速度場測量技術(shù)被廣泛應(yīng)用于空氣動力學測試領(lǐng)域[1]，是科研工作者探索流體運動規(guī)律的有效工具[2]。

PIV技術(shù)[3]的原理是向流動中添加示蹤粒子，示蹤粒子隨流體一同運動，通過分析相機曝光得到的粒子圖像上局部光流信息（一般由查詢窗口內(nèi)的一組示蹤粒子提供）隨時間的變化得到當?shù)氐乃俣刃畔ⅰTV技術(shù)通過追蹤每個示蹤粒子在雙幀曝光圖像中的位置變化實現(xiàn)速度場的測量。不同的測速原理使得兩種方法適用于不同的場景：PIV技術(shù)要求每個查詢區(qū)至少有三個及以上粒子；而PTV技術(shù)則要求測量區(qū)域中示蹤粒子稀疏，易于識別和匹配。

PIV能夠適用于示蹤粒子濃度較高的場景，但需對查詢窗口內(nèi)的光流信息進行追蹤或匹配，導致其空間分辨率受限；而PTV能夠追蹤單個粒子的跨幀位移，但難以處理示蹤粒子間距小、跨幀位移大的場景。近年來，PIV和PTV技術(shù)手段迅猛發(fā)展[4-11]，極大地提高了測量性能，但眾多學者仍希望能通過混合使用PIV和PTV以更進一步地提高測量性能[12-16]。Cowen[13]開發(fā)了DPTV技術(shù)， 通過使用PIV的估計值以限制在第二幅圖像中粒子的搜索窗口，使得PTV可以在較高的粒子密度下進行匹配。Kim等[14]將PIV結(jié)果用于確定PTV過程中的匹配參數(shù)和結(jié)果驗證，提高了粒子匹配效率。Benkovic等[15]將基于匹配概率的松弛算法集成到基于視覺的特征關(guān)聯(lián)概念中, 通過仿真證明了結(jié)合PIV預(yù)分析的混合方法有助于減少錯誤向量的數(shù)量。此外，Zhu等[16]對風沙兩相流中離散相顆粒和連續(xù)相示蹤粒子進行分相處理，并分別使用PIV和PTV技術(shù)，實現(xiàn)了風沙兩相速度場的同步測量。這些混合算法有效地拓寬了粒子圖像測速能力，提高了測量的精度。

另外，在復(fù)雜流動中，示蹤粒子在流場空間難以均勻分布的問題十分常見。例如，在邊界層流動測量實驗中，示蹤粒子難以進入邊界層內(nèi)的高剪切區(qū)，致使示蹤粒子在邊界層內(nèi)的分布過于稀疏，無法進行有效的互相關(guān)計算，導致PIV得到的速度場在近壁區(qū)域存在顯著誤差。針對這一問題，通常的解決手段是提高示蹤粒子的整體濃度。然而，該方法對流動分離區(qū)、旋渦內(nèi)部和激波后方流區(qū)的示蹤粒子濃度改善效果有限。Brooks等[17]通過壁面粒子注射裝置直接從壁面處將示蹤粒子注入到邊界層內(nèi)的近壁流區(qū)，提高近壁區(qū)域的粒子密度。這種方法在邊界層測量中取得了一定的效果，但不具有泛化能力。

針對上述問題，本文提出了一種基于粒子圖像分割的混合PIV-PTV算法。該算法使用維諾多邊形圖像分割方法，對粒子圖像按粒子分布濃度劃分出粒子稀疏區(qū)和稠密區(qū)，對稀疏區(qū)和稠密區(qū)分別調(diào)用PTV和PIV算法，實現(xiàn)示蹤粒子分布不均情況下的高精度流場測量。

1 算法介紹

1.1 算法流程

本文所發(fā)展的基于粒子圖像分割的混合PIVPTV算法主要分為以下4個步驟：示蹤粒子識別、應(yīng)用維諾多邊形計算示蹤粒子密度、設(shè)定閾值并使用支持向量機（support vector machine, SVM）進行粒子圖像劃分、分別調(diào)用PIV和PTV計算速度場。其算法概述如圖1所示，本節(jié)將對實現(xiàn)該算法的幾個關(guān)鍵步驟進行介紹。

圖1 混合PIV-PTV算法技術(shù)路線Fig. 1 Hybrid PIV-PTV algorithm flow chart

1.2 粒子識別

Ohmi等[18]首先提出了基于動態(tài)閾值的示蹤粒子識別方法。根據(jù)粒子圖像上每個粒子的平均亮度水平自適應(yīng)地調(diào)整閾值，消除了傳統(tǒng)方法基于單一或多個閾值解帶來的人工確定閾值的不確定性。Mikheev等[19]對此方法進行了改進，其基本思想是將每個粒子圖像的分割閾值設(shè)置為峰值的比例值，而不是原始算法中的平均值。近幾年基于機器學習等方法的粒子識別技術(shù)也得到迅猛發(fā)展[20]，并展現(xiàn)出了優(yōu)異的性能。綜合考量現(xiàn)有方法的易用性和魯棒性，本文選用經(jīng)典的改進動態(tài)閾值法[19]，其總體方案如圖2所示。

圖2 改進動態(tài)閾值法實現(xiàn)流程[19]Fig. 2 Process of improved dynamic threshold method[19]

首先使用一定的閾值進行濾波，以排除噪聲背景的干擾；尋找濾波后粒子圖像的所有局部極大值，即亮度峰值點，如圖2（b）。然后將每個峰值像素點擴展一個像素大小的邊界，并檢查每個新像素。邊界像素強度與峰值像素強度之比應(yīng)大于預(yù)設(shè)的對比度閾值，即：

式中：Ii,j為新添加邊界像素的灰度值，Imax為峰值點的灰度值，Cthr為對比度閾值。如果滿足此條件，則認為該像素與峰值像素點屬于同一個粒子；否則，則認為該像素不屬于該粒子。這個擴展過程持續(xù)向外進行，直到不存在滿足要求的像素點。擴展過程中若遇到灰度值反常增大，則認為出現(xiàn)拐點，沿此方向的擴展停止。這一粒子識別算法使用了更合理的閾值準則，并能通過拐點檢測在擴展過程中自動分離出重疊的粒子，具有很高的魯棒性。

1.3 粒子濃度定義

由粒子圖像測速原理可知，示蹤粒子的濃度是決定PIV-PTV測速精度的重要因素。本文使用一種基于維諾多邊形劃分的粒子局部濃度評估準則。維諾分析廣泛應(yīng)用于氣象學、細胞學和建筑結(jié)構(gòu)領(lǐng)域，是一種特殊的空間劃分方法。Monchaux等[21]首先將維諾分析引入到顆粒湍流領(lǐng)域，其后得到廣泛應(yīng)用[22-25]。對于二維平面而言，維諾劃分可將包含N個點的平面唯一劃分成N個維諾多邊形，每個多邊形包含且僅包含一個點。對于某一特定的粒子空間分布，維諾劃分具有客觀性和唯一性。基于這一特性，可以用維諾單元面積的倒數(shù)衡量局部粒子濃度。

維諾多邊形生成流程如圖3所示。通過粒子識別得到示蹤粒子的中心點坐標，對點集生成Delaunay三角形網(wǎng)，Delaunay三角形的外接圓圓心是與三角形相關(guān)的維諾多邊形的一個頂點，連接各頂點即可得到維諾圖。

圖3 維諾多邊形生成流程Fig. 3 Voronoi polygon generation process

用包含單個示蹤粒子的維諾多邊形面積的倒數(shù)來表征當?shù)氐牧Ｗ訚舛龋Q為DBV（density based on voronoi）密度， 其定義為：

其中，S(pi)為第i個維諾多邊形的面積，單位為pixel2。可以看到，DBV密度度量與PPP（particles per pixel）密度度量本質(zhì)相同。基于每個粒子當?shù)氐腄BV密度和設(shè)定的閾值，可對每個粒子當?shù)貪舛冗M行稀疏或稠密的判定，進一步在粒子稀疏區(qū)域應(yīng)用PTV，在粒子稠密區(qū)域使用PIV。然而，對一個非常小的區(qū)域應(yīng)用PTV是沒有意義的。若將稀疏或稠密看作是每個粒子的特征，我們希望找到一條稀疏區(qū)和稠密區(qū)之間的光滑分割線。

1.4 基于高斯核的支持向量機

SVM是一種二分類算法，其基本模型是定義在特征空間上的間隔最大的線性分類器[26]。數(shù)據(jù)集根據(jù)自身特征可分為如圖4所示的線性可分、軟間隔可分和線性不可分3種情況[27]。

圖4 3種數(shù)據(jù)集類型示意圖Fig. 4 Three types of data diagram

以線性可分數(shù)據(jù)集為例，給定一些數(shù)據(jù)點，它們分別屬于兩個不同的類，用X表示數(shù)據(jù)點，用Y表示類別（Y取1或者-1，分別代表兩個不同的類），SVM學習目標便是要在空間中通過間隔最大化策略找到一個最優(yōu)的分離超平面（hyper plane）。

線性可分問題的處理手段已經(jīng)非常成熟，而對于非線性可分問題，可通過映射函數(shù)將原空間映射至線性可分新空間。映射函數(shù)由核函數(shù)確定。設(shè)H為特征空間，如果存在一個從X到H的映射函數(shù)?(x):X→H，使得對所有x1,x2∈X, 函數(shù)K(x1,x2)滿足條件K(x1,x2)=?(x1)·?(x2)，則稱K(x1,x2)為核函數(shù)。一個常用的核函數(shù)是高斯核：

其中σ為高斯核函數(shù)的超參數(shù)。SVM中還存在另一個超參數(shù)C，C越大則代表模型對誤分類的懲罰越大。σ和懲罰參數(shù)C的設(shè)置會顯著影響劃分結(jié)果，本文選用的超參數(shù)基于網(wǎng)格法進行尋優(yōu)得到，分別為σ= 0.2，懲罰參數(shù)C= 0.0023。在實際應(yīng)用中需要根據(jù)具體情況調(diào)整，以使得最終的劃分符合預(yù)期。

2 粒子圖像分割

2.1 仿真粒子圖像合成

由于真實測速實驗中無法獲得先驗的參考速度場，在仿真測試階段，為了驗證算法的有效性，按照驗證PIV-PTV算法的一般慣例，使用人工生成仿真粒子圖像進行算法測試。相比以往的仿真粒子圖像[28-30]，本文合成的粒子圖像具有顯著的粒子分布不均的特征。其合成步驟如下：

1）由二維高斯函數(shù)分布生成單個虛擬粒子的灰度分布，以制備粒子集備用。粒子灰度分布函數(shù)為：

其中，(x0,y0)為粒子的中心位置，粒子灰度等級的峰值強度為I0∈(150,250)，粒子直徑dτ∈(2,5)。

2）在大小為512 pixel × 512 pixel的灰度矩陣（第一幀粒子圖像）中隨機劃定部分區(qū)域為稀疏區(qū)域，其余部分設(shè)定為稠密區(qū)域。

3）在稀疏區(qū)域和稠密區(qū)域分別以PPP ＜0.002的密度和PPP = 0.004的密度隨機放置粒子，并記錄粒子坐標，產(chǎn)生第一幀粒子圖像。

4）建立一個512 × 512的標簽矩陣，在該標簽矩陣中由步驟3定義的稠密區(qū)的標簽記為0、稀疏區(qū)的標簽記為1。

5）給定速度場，將步驟3中第一幀圖像中的每個粒子按照當?shù)厮俣冗M行跨幀移動，產(chǎn)生第二幀PIV圖像。

6）在產(chǎn)生的兩幀粒子圖像中添加噪聲，使得生成的粒子圖像更接近真實實驗圖像。最終輸出具有粒子稀疏分布的PIV圖像對和其對應(yīng)的標簽矩陣。

2.2 模型評價指標IoU

能否正確地劃分稀疏區(qū)和稠密區(qū)，是決定混合PIV-PTV算法表現(xiàn)的關(guān)鍵，為此定義稀疏區(qū)交并比（intersection over union, IoU），用以定量評估SVM算法對稀疏區(qū)劃分的合理性。IoU系數(shù)是一種集合相似度度量函數(shù)，通常用于計算兩個樣本的相似度，取值范圍在[0,1]。其公式如下：

A∩B是A和B之間的交集，A∪B是A和B的并集。在本文研究的問題中，A代表預(yù)設(shè)的稀疏區(qū)域，B表示SVM劃分出的稀疏區(qū)域。IoU數(shù)值越高表明A區(qū)域與B區(qū)域重合程度越高，SVM的劃分越準確。IoU具有尺度不變性和非負性等特性，在后文，我們將以IoU為評價指標，探究SVM劃分稀疏區(qū)的表現(xiàn)。

2.3 閾值的選擇

閾值的選擇會深刻影響劃分結(jié)果，應(yīng)充分考慮現(xiàn)有PIV和PTV算法的處理能力。若現(xiàn)有PIV算法要求計算區(qū)域密度至少為ρ1才能進行有效的互相關(guān)計算，而PTV方法能處理的最高密度為ρ2，則此時粒子濃度的閾值ρ的可選擇范圍應(yīng)為ρ1＜ρ＜ρ2。通常情況下，PIV算法要求在一個查詢區(qū)（32 pixel × 32 pixel）內(nèi)示蹤粒子數(shù)目應(yīng)至少大于3；反之，若粒子數(shù)目低于3，則認為使用PTV技術(shù)得到的速度場更有參考價值。基于此，本文工作中以DBV密度ρ= 0.003作為臨界閾值。

2.4 仿真PIV圖像分割

使用2.1節(jié)仿真粒子圖像合成方法生成的一幀粒子圖像如圖5所示，其速度呈蘭金渦結(jié)構(gòu)。蘭金渦的渦核區(qū)域為預(yù)設(shè)的粒子稀疏區(qū)，其邊界如圖6中橙色線條所示。渦心位置為（295, 293），半徑a為92，速度分布在極坐標系下表示為：

圖5 粒子圖像劃分Fig. 5 Illustration of particle image partitioning

圖6 粒子圖像的維諾圖Fig. 6 The Voronoi diagram of the tracer particles

其中，Γ為渦環(huán)量，r為任一點到渦心的距離。

對生成的仿真粒子圖像進行粒子識別，并生成維諾多邊形，如圖6所示。根據(jù)維諾多邊形計算DBV密度，并以密度ρ= 0.003為閾值對粒子進行分類。圖6中，紅色點代表該示蹤粒子的DBV密度小于設(shè)定閾值，即為當?shù)叵∈瑁缓谏c代表該示蹤粒子的DBV密度大于設(shè)定閾值，即為當?shù)爻砻堋？梢钥吹剑t色點主要分布在圖像中心區(qū)域，即中心區(qū)域存在一個較大的稀疏區(qū)，在圖像邊緣處也分布著許多小的稀疏區(qū)。

基于高斯核函數(shù)的SVM尋找兩類粒子的光滑分割邊界。圖6中，紫色曲線為SVM劃分得到的稀疏區(qū)與稠密區(qū)邊界，橙色曲線為預(yù)先劃定的粒子稀疏區(qū)域的邊界。使用公式（5）計算預(yù)設(shè)稀疏區(qū)與預(yù)測稀疏區(qū)交并比IoU為0.96。基于SVM結(jié)果，得到適宜使用PIV和PTV的區(qū)域如圖5所示。

3 混合PIV-PTV測速仿真分析

基于圖5給出的劃分結(jié)果，對SVM劃分出的粒子稠密區(qū)使用基于多重迭代的光流PIV算法MILK[31]計算速度場，查詢區(qū)大小設(shè)為32 pixel × 32 pixel，步長設(shè)為1 pixel。應(yīng)用PIV算法計算速度場的結(jié)果如圖7中藍色箭頭所示。使用基于蟻群優(yōu)化的混合粒子匹配算法[32]進行PTV計算，該方法可以看作是傳統(tǒng)蟻群PTV算法的一個更新，其核心思想是通過改進的蟻群算法來尋求混合目標函數(shù)在全局中的最小化解，該方法在跨幀位移較大的場景下表現(xiàn)出良好的魯棒性[32]。對劃分出的稀疏區(qū)域應(yīng)用該粒子追蹤算法，其計算結(jié)果如圖7紅色箭頭所示。

圖7 混合PIV-PTV結(jié)果展示Fig. 7 Illustration of Hybrid PIV-PTV results

圖7中可以看到，PIV算法計算結(jié)果與理論速度總體分布具有一致性。圖8給出了PIV算法和PTV算法得到的跨幀位移切向分量的誤差沿徑向r的分布。可以看到，在蘭金渦渦核內(nèi)部，PIV結(jié)果的平均誤差在1 pixel附近，顯著高于PTV的誤差（約0.1～0.5 pixel）。在距渦心75～100 pixel半徑區(qū)域內(nèi)，PIV算法誤差較大，究其原因有兩點：一是該區(qū)域示蹤粒子稀疏，難以進行有效的互相關(guān)計算；二是因蘭金渦的存在，該區(qū)域具有速度梯度較大的剪切流動，平均效應(yīng)的存在導致PIV計算產(chǎn)生較大誤差。而PTV方法適用于粒子稀疏區(qū)域，且對速度梯度不敏感，故表現(xiàn)出較低的計算誤差。由此可見，分區(qū)域分別應(yīng)用PIV、PTV進行速度場計算是必要的。

圖8 PIV-PTV誤差分布Fig. 8 Error distribution of PIV-PTV

4 混合PIV-PTV測速應(yīng)用實例

在高超聲速平板邊界層流動中，強剪切流使得示蹤粒子難以進入邊界層近壁區(qū)[33]，導致PIV計算速度存在偏差。在本節(jié)中，混合PIV-PTV測速方法被應(yīng)用于來流馬赫數(shù)Ma= 6的高超聲速平板邊界層速度測量中，實驗在航天空氣動力技術(shù)研究院的FD-03高速風洞中進行，實驗示意圖見圖9。如圖9所示，模型主要由玻璃平板和轉(zhuǎn)捩帶組成。其中，玻璃平板長度240 mm，平板前緣距轉(zhuǎn)捩帶10 mm、距模型最前緣60 mm；采用Nd:YAG脈沖式雙曝光激光器，該激光器輸出能量達600 mJ；使用跨幀CCD相機進行流場拍攝，相機分辨率為2048 pixel ×2048 pixel。實驗參數(shù)如表1所示。

表1 實驗參數(shù)與實驗工況Table 1 Experimental conditions

圖9 Ma = 6超聲速PIV實驗示意圖Fig. 9 PIV experimental sketch for flat-plate flow at Ma = 6

為了便于算法驗證，截取實驗圖像部分區(qū)域進行分析，如圖10所示。截取粒子圖像底部，靠近壁面且與壁面平行，大小為1024 pixel × 1024 pixel；對應(yīng)物理空間為11.2 mm × 11.2 mm，成像區(qū)域最上游距模型最前沿235 mm，該位置邊界層厚度為4.4 mm。

圖10 Ma = 6超聲速PIV實驗結(jié)果Fig. 10 PIV experimental image for flat-plate flow at Ma = 6

使用改進的動態(tài)閾值法進行粒子識別并進行維諾多邊形劃分，得到如圖11（a）所示的維諾圖。計算DBV密度并以密度閾值ρ= 0.003進行二分類，使用基于高斯核函數(shù)的SVM進行稀疏區(qū)和稠密區(qū)劃分。考慮到對小區(qū)域應(yīng)用PTV沒有意義，故剔除劃分結(jié)果中面積較小的稀疏區(qū)域，得到最終的劃分邊界如圖11（a）中橙色線條所示。

圖11 超速平板邊界層混合PIV-PTV結(jié)果Fig. 11 Hybrid PIV-PTV results for flat-plate flow at Ma = 6

超聲速平板邊界層混合PIV-PTV測量的瞬時速度場結(jié)果如圖11（b）所示，其中藍色箭頭為PIV測速結(jié)果，紅色箭頭為PTV測速結(jié)果。可以看到，使用基于粒子圖像分割的混合算法，可以實現(xiàn)對粒子稀疏區(qū)域的分割，進而對近壁區(qū)域使用誤差更小的PTV算法，從而緩解近壁粒子濃度低帶來的PIV測量不確定度高、空間分辨率低的問題。

5 結(jié)論

本文發(fā)展了一種基于粒子圖像分割的混合PIVPTV測速方法，用以解決PIV測速實驗中局部區(qū)域示蹤粒子分布不均勻的問題。算法的基本流程是：首先進行單個粒子識別，對每個粒子計算維諾多邊形和DBV密度；設(shè)定閾值對粒子打標簽，使用基于高斯核的SVM、依據(jù)DBV密度對粒子進行分類，從而實現(xiàn)對粒子圖像稀疏區(qū)和稠密區(qū)的分割；對分割出的粒子稀疏區(qū)和稠密區(qū)分別應(yīng)用PTV和PIV算法計算速度場，實現(xiàn)混合PIV-PTV測速。仿真結(jié)果和高速風洞實驗表明，SVM方法可以實現(xiàn)對粒子圖像基于粒子分布密度的分割，混合PIV-PTV測速方法可以顯著提高粒子分布不均場景下的速度測量精度。

需要指出的是，對粒子濃度分布進行簡單的二分類不能總是滿足工程的需求，在后續(xù)研究中可以發(fā)展基于粒子濃度的分級策略，對不同濃度的粒子區(qū)域在應(yīng)用PIV時自適應(yīng)地使用不同的查詢窗口，從而進一步挖掘PIV的后處理潛力。