畢曉君, 何明潔, 于長東, 范毅偉
(1.中央民族大學 信息與工程學院, 北京 100081; 2.哈爾濱工程大學 信息與通信工程學院, 黑龍江 哈爾濱 150001; 3.哈爾濱工程大學 船舶工程學院, 黑龍江 哈爾濱 150001)
粒子圖像測速(particle image velocimetry,PIV)是一種非接觸的、全局的、定量的流體可視化技術[1]。通過從連續粒子圖像中獲取流體的速度場信息,可以幫助研究者更深入地了解復雜流體的物理特性。從粒子圖像中計算速度場的軟件算法主要有互相關法[2]和變分光流法[3]2種。互相關法通過提取圖像對應的2個查詢窗口來執行互相關計算,搜索查詢窗口的互相關最大值來確定位移。最具有代表性的一種互相關算法是多重網格迭代算法(window deformation iterative multi-grid,WIDIM)[4]。由于互相關算法計算的是窗口內的統計平均位移,仍無法實現稠密(像素級別)的速度場估計。另一常用方法是變分光流方法,其優點是易于嵌入先驗物理知識和特定幾何約束知識[5-6],但是光流法在進行變分優化的過程會花費大量時間和計算量,而且對噪聲非常敏感。
此外,在實驗流體力學中,發展兩相流、多相流場景下的PIV算法已成為一個有前景的研究方向。例如,在船舶與海洋工程領域,研究物體入水的兩相流PIV具有廣泛的應用價值[7-8]。針對具有復雜邊緣的物體入水兩相流PIV圖像,對液相區域的速度場計算具有很大的挑戰性。研究者通常手動完成非計算區域的掩模工作,而后期仍通過互相關算法WIDIM計算速度場,但基于查詢窗口準則的互相關法,在不同相位邊緣位置很難獲得準確的稠密速度場[9]。隨著深度學習的發展和成熟,部分研究工作開始將深度學習技術用于單相流PIV估計問題的研究[10-13]。在計算機視覺領域,深度學習光流模型[14-15]在運動估計任務中取得了杰出的性能。受此啟發,Cai等[12]提出改進的PIV-LiteFlowNet-en模型將其成功用于單相流PIV估計,該模型在速度場估計精度和計算效率上具有一定的優勢。然而,目前的深度學習算法在單相流PIV速度場計算中取得了突破性的進展,但是并未在兩相流PIV實驗領域獲得廣泛應用。
本文提出一種同時用于單相流和物體入水兩相流速度場估計的深度學習運動估計器。本文對先進的光流模型RAFT[16]進行改進,是對文獻[17]的進一步擴展。去掉了模型的冗余部分,旨在不影響估計精度的同時,減少模型參數量。為了更好地訓練和優化模型參數,制作了相應的單相流和物體入水兩相流PIV數據集,并對訓練后的模型進行多種測試和評估,結果表明提出的算法不僅能用于單相流速度場估計,也可用于物體入水的兩相流PIV圖像中的液相速度場計算。
本文利用的網絡基礎架構是一種用于光流估計任務的卷積神經網絡循環全對場變換(recurrent all-pairs field transforms,RAFT)的深度學習框架[16]。該網絡結構如圖1所示,主要由3部分組成。

圖1 光流模型RAFT網絡結構Fig.1 The optical flow model RAFT network structure
1)特征提取模塊。
該模塊利用卷積神經網絡來提取圖像的特征。特征編碼器Feature encoder功能為提取輸入第1幀、第2幀圖像中的特征F1和F2,將輸入圖像映射為1/8分辨率的特征圖。Feature encoder包含6個殘差塊,每經過2個殘差塊后特征圖的分辨率變為原來一半,對應的特征通道數分別為64、128、192和256。同時存在的上下文編碼器Context encoder與特征提取模塊的結構基本類似,作用為提取第1幀圖像的上下文信息特征。在光流估計任務中通常會存在遮擋和大位移等問題,Context encoder可以為下一階段的特征信息融合提供更多有用的上下文信息。
2)相關計算層。
扭曲層[18](Warp)和代價體層[19](cost volume)是光流法在計算光流過程中非常重要和關鍵的部分。Warp層的作用是減少圖像特征和的距離,從而解決大位移的光流估計問題。而Cost volume的作用是構建和已經扭曲的像素匹配損失,從而得到粗略的速度矢量估計。而RAFT在結構設計中不再使用Warp操作層,創新性的直接構建4D Cost volume代價體模塊來進行全局的特征匹配搜索。該模塊通過對所有向量進行對應內積來構建4D(W×H×W×H)的相關層。相關計算層由4層不同分辨率的相關金字塔結構構成{C1,C2,C3,C4},從而可以獲得圖像的全局和局部的特征信息。然后相關性查詢操作L被用來在不同分辨率的相關層上進行查詢,生成的特征圖用于后面的迭代式光流計算。
3)更新迭代模塊。
更新迭代模塊利用卷積GRU循環網絡迭代計算出光流。模塊輸入由4部分組成:Context編碼器輸出,Cost volume輸出,上層隱藏狀態信息以及上層迭代過程輸出的光流。每次計算出相對于上次迭代過程輸出的殘差光流Δf,然后與上次迭代輸出的光流進行相加用于目前的流估計:fk+1=Δf+fk。RAFT無論在訓練還是測試的過程中,都可以選擇更新迭代次數,直到估計出最合適的估計光流。
該模型的損失函數為L1范數:

(1)
式中:γ取值0.8;fgt表示真值場。
實際上,基于深度學習的光流估計和PIV估計任務存在一定的區別,其區別主要反映在應用場景和數據方面。光流估計任務所描述的通常是生活場景中具有穩定突出特征的剛性運動或準剛性運動。因此,光流數據集是一種復雜而具有挑戰性的數據集,圖像帶有RGB彩色信息并含有豐富的特征信息。光流數據集中通常包括遮擋、真實光照、大位移等其他經典問題;另一方面,粒子圖片通常是灰度圖像,具有低分辨率、特征信息少的特點。具體的比較如圖2所示。

圖2 粒子圖像和光流數據集圖像Fig.2 Particle image and optical flow dataset image
經過以上分析,本文對RAFT進行如下改進:
1)RAFT特征提取部分的上下文編碼器Context encoder可以捕獲額外的上下文信息,這對于光流估計非常重要。然而, 由于灰度粒子圖像的簡單性,所以上下文網絡提取的信息對PIV的估計幾乎不起作用。因此,本文移除RAFT的Context encoder模塊。同時,將第2幀圖像的提取特征直接饋入到4D Cost volume和更新迭代模塊。根據實驗驗證,移除Context encoder對PIV估計的準確性幾乎沒有影響,同時刪除Context encoder也大大減少了參數量;
2)光流數據集的圖像通常具有高分辨率的特點,相反,構建的PIV數據集中的粒子圖像通常是具有較低分辨率(256×256像素)。因此,粒子圖像經過特征編碼器Feature encoder的8倍下采樣后,特征信息太少,對后面的流推理階段無法起到應有的作用。本文對特征編碼器Feature encoder進行了改進,刪除相應的2個殘差模塊,將其輸出分辨率提高到原始圖像的1/4。在提高分辨率和獲得更多圖像特征信息的同時,也進一步減少了模型的參數量。圖3對比了改進前后特征編碼器的結構,當輸入粒子圖像為256×256大小時,編碼器最后的輸出分辨率為64×64:

圖3 改進前后的特征編碼結構對比Fig.3 Comparison of feature encoder structure before and after improvement
g(Ii)=EF(Ii):R256×256×3→R64×64×256
(2)
式中EF(·)表示特征編碼器函數。
改進的模型如圖4所示,模型的輸入既可以是只包含液相的單相流PIV圖像對,也可以是對非計算區域掩膜后的兩相流圖像對,以重點對液相區域進行速度場估計。

圖4 改進的RAFT網絡結構Fig.4 Improved RAFT network structure diagram
提出的流動估計深度學習模型采用監督學習策略,所以需要具有真實值的數據來訓練和優化模型參數。按照實驗流體力學的方法[11],首先生成粒子圖像以及流動速度場,然后通過速度場對稱地移動粒子的位置以獲得粒子圖像對。注意,本文在粒子圖像中疊加了少量的高斯噪聲以模擬真實的流動條件。將具有不同參數的粒子圖像和流動速度場隨機組合,即可以組成單相流PIV數據集。最后,本文生成的PIV數據集的數量和流場種類分布如表1所示,關于數據集的具體介紹可參考文獻[17]。制作的PIV數據集合包括據集包括超過29 000項的粒子圖像對及速度場真值,流場種類達13種。龐大的數據量和數據的多樣性可以幫助更好的訓練模型,增加模型的泛化能力。

表1 PIV數據集的描述Table 1 Description of the PIV dataset
此外,為充分模擬兩相流物體入水的場景,本文制作了兩相流PIV數據集。首先生成掩膜數據。使用Matlab軟件中的 PIVLab 工具包[23], 對物體入水圖像進行人工精確標定,對非計算區域進行掩膜處理,標定結果以二值化圖片形式儲存,非液相區域像素值為0, 含有粒子的液相區域像素值為1,從而獲取物體入水的掩膜數據。本文研究的入水物體目標是楔形體,其手動掩膜具體標準可參考文獻[24],經人工標定的部分結果示例如圖5所示。

圖5 人工標定的楔形體入水的兩相流PIV圖像示例Fig.5 Samples of manually calibrated two-phase flow PIV image of a wedge entering water
其次,將單相流數據集和生成的掩膜數據相乘。為了使模型更好的學習到兩相流圖像中物體的邊緣特征,本文將掩膜數據和單相流的粒子圖像對及真實速度場隨機地相乘,從而獲得兩相流PIV數據集。掩膜相乘過程為:
I′(i,j)=I(i,j)?M(i,j)
(3)
V′(i,j)=V(i,j)?M(i,j)
(4)
式中:M(i,j)表示掩膜數據;?表示點乘操作;I(i,j)和V(i,j)表示原來的圖像和真實速度場;I′(i,j)和V′(i,j)分別表示掩膜后的圖像和速度場。通過在此數據集訓練中,模型可以更好地學習兩相流PIV圖像的邊緣特征信息,提高對此的泛化能力。部分數據集示例如圖6所示。
本次實驗中采用的硬件環境配置為Intel(R) Core(TM) i7-9700K CPU@3.00 GHz,32 G內存,并采用RTX 1080Ti GPU進行運算加速,操作系統為64位Ubuntu 18.04采用基于Python的深度學習框架PyTorch來完成程序編程。
在訓練過程中,首先將改進的RAFT在表1的單相流PIV數據集進行訓練,使其具有估計多種典型速度場的能力,訓練集和測試集按照9比1的比例進隨機劃分,在訓練過程中,批尺寸大小為6,初始學習率設為λ=10-4,采用AdamW梯度優化器,其衰減權重為0.000 05。經過1.18×106迭代后,模型參數達到收斂;其次,將訓練好的模型在本文構建的兩相流PIV數據集上再進行訓練微調,使得模型進一步有效學習圖像中的相位邊緣特征信息,最后訓練而成的模型稱為PIV-RAFT-2P。
當在具有真實值的數據集上評估時,評估粒子圖像測速算法最 常 采 用 的 是 均 方 根 誤 差(root mean square error,RMSE):
(5)

本文選擇具有代表性的3種圖像測速算法進行對比,一是基于相關分析法的窗口變形迭代網格方法[4](window deformation iterative multi-grid, WIDIM);多尺度金字塔的HS變分光流算法[3]。WIDIM算法的參數設置為三通道多重判讀窗口(判讀窗口大小分別為64、32、16)。HS變分光流法則采用3層金字塔迭代結構,光滑系數在實驗中進行適當調整;3)深度學習基準PIV算法PIV-LiteFlowNet-en[12]。
本文首先在單相流PIV數據集上對包含不同種類的流場的圖像進行了測試。為了增加對比性,這里本文選擇與WIDIM相關、HS光流法和PIV-LiteFlowNet-en進行了精度的比較,如表2所示。

表2 單相流PIV數據集測試誤差Table 2 Single-phase flow PIV dataset test error
如表2所示,本文改進的光流神經網絡PIV-RAFT-2P在不同種類的流場都取得了最高精度的計算結果。除了在復雜流場DNS-turbulence取得0.115的測試誤差外,在其他流場上計算精度都提高了一個量級。尤其在一些簡單的流場測試中,如Sink流場,更是取得了0.017的低誤差結果。這證明了本文提出的方法具有很好的泛化能力和魯棒性。接著通過不同方法計算的速度場來比較算法獲取細節渦結構的能力。二維湍流流場DNS-turbulence的粒子圖像是國際公認的 PIV 測試基準之一, 對其進行了測試和比較。
如圖7所示,WIDIM和HS算法所計算的渦心位置不夠連續,與真實值相差較多,而本文的模型所計算的速度場,基本與真實值一致,誤差只有0.076。尤其是在小尺度渦結構的估計上(圈處),神經網絡模型PIV-RAFT的優勢更明顯,這是由于該方法能夠提供稠密速度場。
接著在兩相流PIV數據集上對模型性能進行了評估。圖8展示了不同算法在Back-step, DNS-turbulence和Cylinder流場的估計速度場。由于HS算法和PIV-LiteFlowNet-en在兩相流PIV估計上表現較差,誤將非計算區域當作可測區域,這里對其結果不進行展示。可以看到,本文方法的速度場計算結果與真實流場基本一致。在掩模區域與計算區域交界面位置處,可以對邊界矢量進行有效估計。對于復雜流場DNS-turbulence,PIV-RAFT-2P可以獲得更稠密的速度場(圈處)。這是因為相關法WIDIM是基于窗口查詢的方式,獲得的速度場是稀疏的形式,而本文的方法則可以獲得更多的小尺度渦結構。

圖8 不同方法在兩相流圖像估計的速度場Fig.8 The velocity field estimated by different methods in two-phase flow images
為了定量的分析和比較,本文對包含4種經典的流場的兩相流測試集進行測試并對結果取平均。如表3所示,WIDIM計算的誤差較大,而深度學習估計器PIV-RAFT-2P在不同流場上都取得了最低的誤差,同時本文的方法可以獲得像素級別的高分辨率速度場。

表3 兩相流PIV圖像測試RMSE結果Table 3 Test error of two-phase flow PIV images
本文首先實際流場單相流圖像進行測試,所選取的實際粒子圖像對是渦流對[25]。該圖像成像質量較高(粒子分布均勻、濃度合適、粒直徑合適),是非常適合 PIV 分析的粒子圖像,如圖9、10所示。

圖9 不同方法在渦流對圖像估計的速度場Fig.9 Velocity fields estimated from laboratory PIV images

圖10 不同方法在兩相流圖像估計的速度場Fig.10 The velocity field estimated by different methods in two-phase flow images
可以看到,由于受到噪聲影響,HS光流法估計的速度場不夠光滑,結果較差。而WIDIM算法因為其發展多年,取得比較好的估計效果,但是其因為計算得到低分辨率的速度場,難以獲得渦心位置的細節流動。本文的方法也可以取得和WIDIM方法相一致的效果,獲得了光滑的速度場和更多的細節信息,尤其在渦心位置和在下面的噴嘴位置也可以得到有效估計。
然后在PIV實驗系統平臺中拍攝獲得實際的楔形不同時刻入水的粒子圖像對(t=5 ms, 10 ms, 20 ms),對其進行測試。可以看到2種算法都能展現楔形體堆積區射流根部的高速射流現象,但是PIV-RAFT-2P在堆積區的計算的速度場分布更為平滑,而基于窗口查詢方式計算的WIIDM算法因為獲得的速度場是稀疏的,所以在邊界位置的速度場呈現鋸齒形狀的分布。在楔形體入水不同時刻算法都展現出了相同的趨勢。
本文分別利用大小為256×256的PIV圖像對不同算法進行了計算效率的評估。如表4所示,PIV-RAFT-2P在計算效率上占有很大的優勢。在計算時間上,PIV-RAFT-2P遠遠超過了傳統算法,快于PIV-LiteFlowNet-en。相比于WIDIM,本文的方法可以計算出更多的速度矢量個數。另外,與深度學習方法相比,本文的改進后的模型PIV-RAFT-2P參數量只有3.73 MB,僅為PIV-LiteFlowNet-en的59.8%。

表4 計算效率對比Table 4 Comparison of computational efficiency
1)本文提出的改進的光流卷積神經網絡PIV-RAFT-2P方法,可同時用于單相流和物體入水兩相流的液相區域速度場估計。對光流模型RAFT的上下文編碼器進行移以及將特征編碼器部分的輸出分辨率提高到原來的1/4以獲取更高分辨率的粒子圖像特征;相應的單相流PIV數據集和兩相流PIV數據集用于訓練和優化模型參數。
2)本文提出的模型可以同時用于單相流和物體入水兩相流PIV估計并取得高精度的估計效果,表明模型具有較強的泛化能力。此外,模型還具有參數量少(3.73 MB)、計算速度快的輕量化優勢。
深度學習光流模型在PIV估計任務中具有很大的潛力和優勢。然而,本文兩相流中的液相區域的提取是利用人工掩膜方式,比較耗時和消耗人力,本文將進一步深入研究,提出一種級聯的深度學習框架自動進行非液相區域的掩膜及液相區域的速度場計算。