朱浩然高 琪王洪平廖相巍趙 亮魏潤杰王晉軍

1．浙江大學 航空航天學院,杭州 310027;2．中國科學院力學研究所,北京 100190;3．鞍鋼集團鋼鐵研究院 冶金工藝研究所,遼寧 鞍山 114009;4．海洋裝備用金屬材料及其應用國家重點實驗室,遼寧 鞍山 114009;5．北京立方天地科技發(fā)展有限責任公司,北京 100083;6．流體力學教育部重點實驗室 北京航空航天大學,北京 100191

0 引 言

粒子圖像測速技術(Particle Image Velocity,PIV)是一種應用廣泛的速度測量技術,對于復雜流體流動,可以根據(jù)體PIV中三維三分量(3D3C)流動結構進行研究。在所有的3D3C測量方法中,Elsinga等[1]提出的層析PIV技術(Tomographic PIV,Tomo-PIV)已經(jīng)被證明能夠在相當高的粒子播種密度(每像素0．05個粒子)下進行精確的空間分辨率測量。粒子重構(Particle Reconstruction,PR)是Tomo-PIV的關鍵步驟,它是解決從二維粒子圖像到三維粒子分布的逆投影問題的過程。在Elsinga等[1]的Tomo-PIV原始文章中,提出了基于倍增代數(shù)重構技術(Multiplicative Algebraic Reconstruction Technique,MART)的三維粒子重構方法。從那時起,許多先進的技術被陸續(xù)開發(fā)出來用來優(yōu)化三維粒子重構,以提高重構精度或效率,這已經(jīng)被Scarano[2]和高琪等[3]充分驗證。大多數(shù)可用的粒子重構技術都是基于MART算法的,例如空間濾波MART(Spatial Filtering MART,SFMART),它在每次MART迭代后對重構的粒子強度場進行空間濾波(Discetti等)[4]。SF-MART算法比傳統(tǒng)的MART算法具有更好的重構質量,該算法將在本文中用來與新技術進行測試和比較。

對于粒子重構問題,隨著粒子播種密度的增加,在視線交叉處會產(chǎn)生虛假粒子,重構質量將急劇下降。在高粒子濃度情況下,為了提高重構質量,目前已經(jīng)有很多算法。Worth和Nickels[5]將乘性第一猜測(Multiplicative First Guess,MFG)作為標準MART方法的先驗條件,它可以提供一個相當精確的解作為MART迭代的初始條件,并加快收斂速度。Atkinson和Soria[6]進一步提出了一種乘性LOS(Multiplicative LOS,MLOS)估計來確定可能的粒子位置,而不需要加權矩陣作為MFG。MLOS除了具有良好的初始化水平外,還可以有效去除虛假粒子從而顯著提高重構質量。通過分析峰值強度和軌跡長度的聯(lián)合分布,可以在某些情況下成功分離虛假粒子和實際粒子(Elsinga和Toggoz)[7]。de Silva等[8]提出的基于模擬匹配的重構增強(Simulacrum Matching-based Reconstruction Enhancement,SMRE)

技術利用實際粒子的特征形狀和大小來去除重構強度場中的虛假粒子。抖盒子(Shake-The-Box,STB)方法(Schanz等)[9-10]使用已知軌跡來預測粒子分布。Wieneke[11]提出的體積粒子分布(Iterative Reconstruction of Volumetric Particle Distribution,IPR)迭代重構來校正粒子位置。與MART相比,STB在精確度和粒子濃度方面都有很大的提高。對于時間分辨圖像采集,Lynch和Scarano[12]提出的序列運動跟蹤增強MART方法(Sequential Motion Tracking Enhancement MART,SMTE-MART)通過建立在前一時刻的對象進行重構,也得到了一個基于增強猜測的時間推進估計目標強度場。與MART和MTE-MART相比,該方法具有更高的重構質量和更高的速度場測量精度(Novara等)[13]。對于三維粒子場重構,有學者提出了一些新的重構方案。強度增強MART方法(Intensity-Enhanced MART,IntE-MART)使用基于直方圖的強度降低來抑制重影的強度(王洪平等)[14]。Gesemann等[15]使用基于約束最小二乘策略和L1正則化的優(yōu)化算法求解體積強度。葉志堅等[16]提出了一種用于粒子重構的雙基追蹤方法,該方法在二維模擬實驗中與MART相比具有更高的重構質量。Bajpayee和Techet[17]提出了一種基于單應匹配合成孔徑重聚焦的高效方法。Ben Salah等[18]提出了一種“面向對象”的方法,稱為基于標記點過程的迭代檢測對象體積重建方法(Iterative Object Detection-Object Volume Reconstruction based on Marked Point Process,IOD-OVRMPP),用于粒子群總體的三維重構,用這種方法可以直接得到粒子的位置。

當前一些傳統(tǒng)的粒子重構方法雖然已經(jīng)取得比較理想的效果,但是在粒子重構過程中還存在一定的虛假粒子,影響粒子重構質量。隨著機器學習在圖像處理領域的發(fā)展,設計一個基于機器學習的模型來處理與圖像相關的各種問題已成為一個熱門話題。近幾年來,神經(jīng)網(wǎng)絡已開始應用于粒子圖像測速。機器學習已經(jīng)被用來代替?zhèn)鹘y(tǒng)的互相關算法進行密集粒子運動的速度場估計(蔡聲澤等)[19-20]。在第13屆粒子圖像測速國際研討會上介紹了一系列工作。例如,Lagemann等[21]將卷積神經(jīng)網(wǎng)絡(CNN)應用于PIV,取得了與傳統(tǒng)互相關算法相似的效果。然而,目前將機器學習應用于PIV的研究大多還處于二維階段,將機器學習作為一種完全三維的應用用于粒子重構還缺乏研究。本文利用3D-CNN設計了一個機器學習算法來解決三維粒子重構問題,并嘗試進一步提高粒子重構質量,該算法稱為AI-PR(Artificial Intelligence Particle Reconstruction)。

1 卷積神經(jīng)網(wǎng)絡的概念

卷積神經(jīng)網(wǎng)絡(CNN)是人工神經(jīng)網(wǎng)絡(Artificial Neural Network,ANN)的一種,作為計算機視覺技術的有力工具已經(jīng)得到了廣泛的應用[22-23]。一般來說,一個深度神經(jīng)網(wǎng)絡由多層(包括輸入、輸出和隱藏層)組成,并且有許多可訓練的參數(shù)。原始信息(如圖像)通過輸入層輸入到CNN,預測結果(如分類結果)通過輸出層給出。所有的信息通過網(wǎng)絡內部的隱藏層進行傳輸和處理。神經(jīng)網(wǎng)絡的應用層一般包括卷積層、池化層、激活層和反卷積層等。

其中,卷積層是卷積神經(jīng)網(wǎng)絡結構的核心模塊,其中輸出的每個神經(jīng)元(yj)僅連接到輸入的一小部分神經(jīng)元(xj),如圖1(a)所示。相反,在全連接層中,任何輸出神經(jīng)元(yj－1)與每個輸入神經(jīng)元之間存在連接,如圖1(a)中的黃色虛線箭頭所示。因此,卷積層的正向計算通常更有效,涉及的可訓練參數(shù)更少。輸入和輸出之間的關系可以被表示為Y＝WX＋b,輸入向量X＝(x1,x2,．．．,xi)T,輸出向量Y＝(y1,y2,．．．,yj)T,對于卷積層,權重矩陣W是稀疏的。

圖1(a)中的局部連接性本質上是由卷積運算確定的,圖1(b)中用一個簡單的例子說明了卷積計算的原理(b＝0)。灰色立方體為過濾器或卷積核,其大小為3×3×4,其深度始終與輸入層相同(輸入層為5×5×4的藍色立方體)。卷積核以一定的步幅(這里取1)在水平和垂直方向上移動,并將相應的輸入元素乘以3×3×4卷積核內相應的系數(shù)。然后,將36組系數(shù)的乘積相加,所得的總和作為相應位置的輸出值(青色正方形)。輸出立方體的深度等于卷積核數(shù),為便于說明,示例中卷積核數(shù)設為1,也可以使用多個卷積核來從輸入中提取更多的特征。在卷積神經(jīng)網(wǎng)絡中,卷積核是需要學習的未知量。

圖1 簡單卷積層示意圖Fig.1 Sketch of the Simple convolutional layer

一旦卷積神經(jīng)網(wǎng)絡結構建立起來,就可以通過基于梯度下降的技術,通過最小化損失來訓練網(wǎng)絡的參數(shù)。由于訓練數(shù)據(jù)量往往很大,為了提高訓練效率,需要采用特殊的優(yōu)化算法(與經(jīng)典的參數(shù)優(yōu)化問題相比),例如:小批量梯度下降、動量梯度下降[24]、均方根支柱(RMSprop)[25]、自適應矩估計(Adam)[26]等,這些算法已在卷積神經(jīng)網(wǎng)絡中得到了廣泛的應用并取得了良好的效果。

2 使用機器學習方法進行粒子重構

神經(jīng)網(wǎng)絡的訓練主要使用合成的三維粒子場作為“真實粒子場”和將其投影成多個二維圖像后再使用MLOS方法重構出的三維粒子場作為神經(jīng)網(wǎng)絡的訓練集。其中,MLOS方法被認為是一種非常高效的粒子場初始猜測方法,將其作為神經(jīng)網(wǎng)絡的輸入可以提高計算的效率和成功率。通過多次迭代訓練,獲取MLOS粒子場和真實粒子場二者之間的虛擬映射關系。通過輸入MLOS重構的初始三維粒子場,利用機器學習方法在訓練過程中學習到的映射關系,得到接近“真實粒子場”的結果(如圖2所示)。即本文提出的方法主要通過兩個步驟來獲得最終的粒子場。第一步通過相機成像從多個二維粒子圖像中計算出初始三維粒子場(EMLOS),這與傳統(tǒng)的PR算法相同。第二步用EMLOS作為輸入,使用訓練好的粒子重構的卷積神經(jīng)網(wǎng)絡計算得到最終的粒子三維分布。

圖2 機器學習方法訓練和預測示意圖Fig.2 Schematic diagram of machine learning method training and prediction

機器學習方法中所使用的卷積神經(jīng)網(wǎng)絡的結構如圖3所示,通過建立12層深度卷積神經(jīng)網(wǎng)絡來學習三維粒子場。輸入輸出層尺寸均為64×64×32×1。最后一層的激活函數(shù)為Sigmoid函數(shù),其他層為修正線性單元(Re LU)函數(shù)。除輸入層和輸出層外,其余層的尺寸均為64×64×32×16,并均采用批量歸一化法[27]。輸入層的卷積核大小為3×3×3×1,其他層的卷積核大小為3×3×3×16。

圖3 神經(jīng)網(wǎng)絡結構示意圖Fig.3 Schematic diagram of neural network structure

在實驗中,通常無法通過測量獲得粒子的精確位置和強度分布。因此,采用合成粒子場作為訓練和測試數(shù)據(jù)。合成粒子場及其圖像的生成使用了PR算法測試中常用的典型方法,詳見王洪平[14]和葉志堅[16]等的研究。主要根據(jù)給定的映射函數(shù)計算粒子場的4個投影,模擬攝像機成像。然后計算初始MLOS字段,并將其作為上述3D-CNN的輸入。

為了提高卷積神經(jīng)網(wǎng)絡的魯棒性,我們在用于訓練的20%的粒子圖像中加入了高斯噪聲。并在4幅粒子圖像中加入不同程度的高斯噪聲。按照添加噪聲的典型方法[14,19],圖像噪聲的標準差σ用PR測試的nσ值來計算,其中噪聲強度n的取值區(qū)間為0～0．20,間隔為0．05。實驗結果表明:當訓練樣本數(shù)大于500時,新算法具有足夠的穩(wěn)定性和精確度。訓練網(wǎng)絡的損失函數(shù)定義為:

式中:F是真值目標矩陣,Fnn是網(wǎng)絡輸出矩陣,皆表示粒子空間灰度的矩陣;ε是一個小值,以防止分母為零,文章中ε取10－3;運算符“·”表示矩陣的內積,即矩陣中對應元素的乘積;sum()表示矩陣所有元素的灰度值總和;i是訓練樣本的索引序號;M是訓練樣本的數(shù)目。

3 結果與討論

由于粒子重構技術很難用真實的實驗數(shù)據(jù)直接進行驗證測試,為了得到準確的比對結果,使用人工合成數(shù)據(jù)對AI-PR算法進行了驗證測試。用與生成訓練集相同的方式生成780×780×140測試粒子場。粒子播種密度ppp取值為0．05～0．30,間隔為0．05。噪聲強度n取值為0．05～0．30,間隔為0．05。使用3種傳統(tǒng)的PR方法:MLOS、MART-5次迭代和MART-10次迭代的方法,與本文提出的機器學習方法進行了比較。所有的訓練和測試都是在Python編程的Tenserflow TMV1．13．1(Abadi等)[28]框架下進行的(www．python．org),并且使用到了Matlab編程軟件(Math Works公司)。使用的計算機是一臺Intel x99工作站,包含一個CPU為E5-2696 V4,64GB DDR4內存和一個RTX2080ti圖形處理單元。

圖4提供了ppp＝0．15的重構粒子場的中心截面圖。很明顯,MLOS只給出了一個非常粗略的粒子位置和強度分布的初始猜測,而AI-PR和MART可以更好的重構粒子場。進一步比較AI-PR和MART方法,可以發(fā)現(xiàn)MART比AI-PR產(chǎn)生更多的虛假粒子,并且強度分布更差。如果仔細觀察粒子形狀,可以發(fā)現(xiàn)MART重構的粒子具有更多的橢球形狀,而AI-PR則可以更好地恢復成球形形狀。

圖4 粒子場截面圖Fig.4 Cross-sections of particle field

如圖5所示,在重構質量方面,AI-PR方法要優(yōu)于MART方法。用質量因子Q來評價新技術的準確性和穩(wěn)定性,即合成場與重構場之間的相關系數(shù)。圖5(a)為無成像噪聲情況下測得所有方法質量因子Q隨粒子播種密度的變化。結果表明,AI-PR可以有效地從MLOS場中恢復粒子。它比MART方法有更好的Q值。當粒子播種密度為0．25時,AI-PR的Q值保持在0．7左右,而MART-10次迭代則降至0．6以下。在圖5(b)粒子播種密度為0．15的樣本中加入了噪聲的影響,從圖中可以發(fā)現(xiàn)對于所有方法,Q隨著噪聲強度的增加而減小,但是AI-PR對偏差的穩(wěn)定性最好。

圖5 幾種方法的質量因子Q變化圖Fig.5 Variation of quality factor Q of several methods

就計算效率而言,MLOS、MART-5次迭代、MART-10次迭代和AI-PR算法的計算時間分別為512．5、5333．5、9881．5和524．5 s。由 于AI-PR處 理包含了MLOS和CNN的計算時間,實際計算時間僅為12 s左右,但值得注意的是,前期卷積神經(jīng)網(wǎng)絡的訓練耗時約為16 h,并且在測試中MART算法并未使用GPU加速。

除了評估新方法統(tǒng)計的品質因子之外,還針對粒子重構中經(jīng)常出現(xiàn)的虛假粒子(Ghost particle)問題進行了測試。雖然虛假粒子也會在AI-PR中出現(xiàn),其出現(xiàn)概率會隨著粒子濃度增加而增大,但是其出現(xiàn)率非常低,以至于影響粒子重構質量因子的主要原因是灰度值分布的偏差。因為本文方法是以MLOS方法的初始場作為輸入進行粒子場重構優(yōu)化,其對粒子場是一個非常初步的猜測,只是為了實現(xiàn)粒子空間稀疏性的猜測,而最終粒子都將分布在MLOS初始場的灰度非零體素上,幾乎不會出現(xiàn)粒子丟失現(xiàn)象。而機器學習有別于傳統(tǒng)Tomo-PIV標定獲得固定的映射函數(shù),其能通過大量樣本實現(xiàn)未知映射約束關系的學習,因此其求解稀疏欠定問題的能力更強,甚至可以求解1或0范數(shù)的優(yōu)化問題?；谶@些原因,機器學習能很好地避免虛假粒子的出現(xiàn),尤其是在大樣本訓練集的前提下。

此外,有一點值得注意的是:目前測試數(shù)據(jù)為人工合成數(shù)據(jù)。最理想的方法是通過真實實驗來驗證重構結果,但由于在真實實驗中難以得到真實粒子空間分布的信息,所以難以通過實驗來進行算法驗證。另一方面,粒子重構質量對流場統(tǒng)計特性的影響效果過于復雜且未知,很難評估具體哪個環(huán)節(jié)造成誤差,因此也沒有采用通過類似湍流邊界層統(tǒng)計特性的評估來驗證粒子重構算法。新方法在真實實驗中的應用效果將在后續(xù)探索中進一步深化。

4 結 論

總的來說,基于機器學習方法的三維粒子重構技術在重構質量方面優(yōu)于傳統(tǒng)的基于MART的算法。

1)AI-PR技術相較于基于MART的傳統(tǒng)算法,可以進一步剔除虛假粒子,并從二維粒子圖像中有效的恢復粒子位置和形狀,得到更加準確的粒子場分布。

2)AI-PR技術相較于基于MART的傳統(tǒng)算法,具有更高的運算處理效率。

基于AI-PR技術的優(yōu)勢,該技術很有希望應用于更真實的實驗。同時需要注意的是,對于不同的成像映射函數(shù)的實驗,AI-PR需要在整個訓練過程中來學習投影規(guī)則。今后的研究可以集中在體PIV的標定與AI-PR訓練相結合,或者直接從AI-PR中進行粒子重構,而不需要標定以及針對不同實驗情況來進行附加網(wǎng)絡訓練。