朱浩然高 琪王洪平廖相巍趙 亮魏潤(rùn)杰王晉軍

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

0 引 言

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

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

技術(shù)利用實(shí)際粒子的特征形狀和大小來去除重構(gòu)強(qiáng)度場(chǎng)中的虛假粒子。抖盒子(Shake-The-Box,STB)方法(Schanz等)[9-10]使用已知軌跡來預(yù)測(cè)粒子分布。Wieneke[11]提出的體積粒子分布(Iterative Reconstruction of Volumetric Particle Distribution,IPR)迭代重構(gòu)來校正粒子位置。與MART相比,STB在精確度和粒子濃度方面都有很大的提高。對(duì)于時(shí)間分辨圖像采集,Lynch和Scarano[12]提出的序列運(yùn)動(dòng)跟蹤增強(qiáng)MART方法(Sequential Motion Tracking Enhancement MART,SMTE-MART)通過建立在前一時(shí)刻的對(duì)象進(jìn)行重構(gòu),也得到了一個(gè)基于增強(qiáng)猜測(cè)的時(shí)間推進(jìn)估計(jì)目標(biāo)強(qiáng)度場(chǎng)。與MART和MTE-MART相比,該方法具有更高的重構(gòu)質(zhì)量和更高的速度場(chǎng)測(cè)量精度(Novara等)[13]。對(duì)于三維粒子場(chǎng)重構(gòu),有學(xué)者提出了一些新的重構(gòu)方案。強(qiáng)度增強(qiáng)MART方法(Intensity-Enhanced MART,IntE-MART)使用基于直方圖的強(qiáng)度降低來抑制重影的強(qiáng)度(王洪平等)[14]。Gesemann等[15]使用基于約束最小二乘策略和L1正則化的優(yōu)化算法求解體積強(qiáng)度。葉志堅(jiān)等[16]提出了一種用于粒子重構(gòu)的雙基追蹤方法,該方法在二維模擬實(shí)驗(yàn)中與MART相比具有更高的重構(gòu)質(zhì)量。Bajpayee和Techet[17]提出了一種基于單應(yīng)匹配合成孔徑重聚焦的高效方法。Ben Salah等[18]提出了一種“面向?qū)ο蟆钡姆椒?稱為基于標(biāo)記點(diǎn)過程的迭代檢測(cè)對(duì)象體積重建方法(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ì)量。隨著機(jī)器學(xué)習(xí)在圖像處理領(lǐng)域的發(fā)展,設(shè)計(jì)一個(gè)基于機(jī)器學(xué)習(xí)的模型來處理與圖像相關(guān)的各種問題已成為一個(gè)熱門話題。近幾年來,神經(jīng)網(wǎng)絡(luò)已開始應(yīng)用于粒子圖像測(cè)速。機(jī)器學(xué)習(xí)已經(jīng)被用來代替?zhèn)鹘y(tǒng)的互相關(guān)算法進(jìn)行密集粒子運(yùn)動(dòng)的速度場(chǎng)估計(jì)(蔡聲澤等)[19-20]。在第13屆粒子圖像測(cè)速國(guó)際研討會(huì)上介紹了一系列工作。例如,Lagemann等[21]將卷積神經(jīng)網(wǎng)絡(luò)(CNN)應(yīng)用于PIV,取得了與傳統(tǒng)互相關(guān)算法相似的效果。然而,目前將機(jī)器學(xué)習(xí)應(yīng)用于PIV的研究大多還處于二維階段,將機(jī)器學(xué)習(xí)作為一種完全三維的應(yīng)用用于粒子重構(gòu)還缺乏研究。本文利用3D-CNN設(shè)計(jì)了一個(gè)機(jī)器學(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)的一種,作為計(jì)算機(jī)視覺技術(shù)的有力工具已經(jīng)得到了廣泛的應(yīng)用[22-23]。一般來說,一個(gè)深度神經(jīng)網(wǎng)絡(luò)由多層(包括輸入、輸出和隱藏層)組成,并且有許多可訓(xùn)練的參數(shù)。原始信息(如圖像)通過輸入層輸入到CNN,預(yù)測(cè)結(jié)果(如分類結(jié)果)通過輸出層給出。所有的信息通過網(wǎng)絡(luò)內(nèi)部的隱藏層進(jìn)行傳輸和處理。神經(jīng)網(wǎng)絡(luò)的應(yīng)用層一般包括卷積層、池化層、激活層和反卷積層等。

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

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

圖1 簡(jiǎn)單卷積層示意圖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)化問題相比),例如:小批量梯度下降、動(dòng)量梯度下降[24]、均方根支柱(RMSprop)[25]、自適應(yīng)矩估計(jì)(Adam)[26]等,這些算法已在卷積神經(jīng)網(wǎng)絡(luò)中得到了廣泛的應(yīng)用并取得了良好的效果。

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

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

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

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

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

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

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

3 結(jié)果與討論

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

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

圖4 粒子場(chǎng)截面圖Fig.4 Cross-sections of particle field

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

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

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

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

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

4 結(jié) 論

總的來說,基于機(jī)器學(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)確的粒子場(chǎng)分布。

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

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