基于機器學(xué)習(xí)方法的三維粒子重構(gòu)技術(shù)

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

2 使用機器學(xué)習(xí)方法進(jìn)行粒子重構(gòu)

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

圖2 機器學(xué)習(xí)方法訓(xùn)練和預(yù)測示意圖Fig.2 Schematic diagram of machine learning method training and prediction

機器學(xué)習(xí)方法中所使用的卷積神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)如圖3所示,通過建立12層深度卷積神經(jīng)網(wǎng)絡(luò)來學(xué)習(xí)三維粒子場。輸入輸出層尺寸均為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)絡(luò)結(jié)構(gòu)示意圖Fig.3 Schematic diagram of neural network structure

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

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

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

3 結(jié)果與討論

由于粒子重構(gòu)技術(shù)很難用真實的實驗數(shù)據(jù)直接進(jìn)行驗證測試,為了得到準(zhǔn)確的比對結(jié)果,使用人工合成數(shù)據(jù)對AI-PR算法進(jìn)行了驗證測試。用與生成訓(xùn)練集相同的方式生成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次迭代的方法,與本文提出的機器學(xué)習(xí)方法進(jìn)行了比較。所有的訓(xùn)練和測試都是在Python編程的Tenserflow TMV1．13．1(Abadi等)[28]框架下進(jìn)行的(www．python．org),并且使用到了Matlab編程軟件(Math Works公司)。使用的計算機是一臺Intel x99工作站,包含一個CPU為E5-2696 V4,64GB DDR4內(nèi)存和一個RTX2080ti圖形處理單元。

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

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

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

圖5 幾種方法的質(zhì)量因子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)絡(luò)的訓(xùn)練耗時約為16 h,并且在測試中MART算法并未使用GPU加速。

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

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

4 結(jié) 論

總的來說,基于機器學(xué)習(xí)方法的三維粒子重構(gòu)技術(shù)在重構(gòu)質(zhì)量方面優(yōu)于傳統(tǒng)的基于MART的算法。

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

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

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