基于壓縮感知的高頻響流場(chǎng)重構(gòu)方法及其應(yīng)用

陳子玉，周楷文，劉應(yīng)征，溫 新

(上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240)

0 引 言

流場(chǎng)的時(shí)空變化包含了重要的流動(dòng)信息，但實(shí)際測(cè)量卻難以同時(shí)得到較高時(shí)間分辨率和空間分辨率的流場(chǎng)。粒子圖像測(cè)速（PIV）[1-4]誕生于20世紀(jì)80年代，是一種可獲得高空間分辨率的光學(xué)流場(chǎng)測(cè)試技術(shù)，但其時(shí)間分辨率通常限于15 Hz以下。盡管利用高速相機(jī)和激光器可以獲得相對(duì)較高時(shí)間分辨率的PIV流場(chǎng)信息，但價(jià)格昂貴，且有時(shí)受實(shí)驗(yàn)條件限制，而實(shí)際應(yīng)用中高頻流場(chǎng)非常常見，因此有必要發(fā)展高頻PIV流場(chǎng)重構(gòu)方法。

為了提高PIV測(cè)量流場(chǎng)的時(shí)間分辨率，國(guó)內(nèi)外諸多學(xué)者進(jìn)行了相關(guān)研究。Druault等[5]對(duì)所測(cè)得的PIV數(shù)據(jù)進(jìn)行本征正交分解，然后直接對(duì)模態(tài)系數(shù)進(jìn)行樣條插值，從而得到高時(shí)間分辨率的流場(chǎng)。這種方法魯棒性好且易于實(shí)施，但這種方法本身依賴于高時(shí)間分辨率采樣，頻率提高有限。之后，許多研究者借用具有高時(shí)間分辨率的局部點(diǎn)測(cè)量技術(shù)，如熱線、麥克風(fēng)等，來提高PIV的時(shí)間分辨率。He等[6]在使用PIV 測(cè)量技術(shù)的同時(shí)，附加了34個(gè)局部點(diǎn)測(cè)量的高頻采樣傳感器，結(jié)合線性隨機(jī)估計(jì)和數(shù)據(jù)同化，重構(gòu)出時(shí)間分辨率的湍流流場(chǎng)，但這種方法需要解控制方程，計(jì)算量較大，且在某些工況時(shí)，無法安裝大量的高頻采樣傳感器。Jin[7]等使用雙向循環(huán)神經(jīng)網(wǎng)絡(luò)，并借助于高頻局部點(diǎn)測(cè)量和PIV測(cè)量，恢復(fù)出時(shí)間分辨率的圓柱繞流流場(chǎng)，但這種方法仍需要消耗較大計(jì)算量，且加入的局部測(cè)量點(diǎn)同樣會(huì)擾亂流場(chǎng)。黃路[8]針對(duì)非定常流場(chǎng)，結(jié)合互相關(guān)算法和卡爾曼濾波，利用仿真驗(yàn)證了時(shí)間分辨率實(shí)時(shí)自適應(yīng)調(diào)整算法的有效性，但這種方法依賴于納秒級(jí)分幅成像的硬件裝置。因此，目前針對(duì)該問題的研究存在許多缺陷，如僅使用PIV測(cè)量的時(shí)間分辨率重構(gòu)依賴于自身的高時(shí)間分辨率測(cè)量或依賴于昂貴的硬件，而結(jié)合侵入式點(diǎn)測(cè)量的重構(gòu)在某些工況下甚至無法安裝點(diǎn)測(cè)量?jī)x器或會(huì)擾亂流場(chǎng)。

考慮上述問題，本文提出了一種基于壓縮感知（CS）從亞采樣的PIV數(shù)據(jù)中重構(gòu)高時(shí)間分辨率流場(chǎng)的方法。該方法無需借助侵入式的高頻點(diǎn)測(cè)量，使用壓縮感知對(duì)亞采樣的POD系數(shù)進(jìn)行高分辨率重構(gòu)，可實(shí)現(xiàn)12倍的頻率放大，且不依賴于高時(shí)間分辨率的PIV測(cè)量，利用常用的PIV測(cè)量裝置即可實(shí)現(xiàn)。該方法將隨機(jī)采樣得到的PIV數(shù)據(jù)進(jìn)行本征正交分解，將所得到的時(shí)間系數(shù)利用壓縮感知進(jìn)行重構(gòu)，繼而得到高時(shí)間分辨率的系數(shù)，結(jié)合分解得到的空間模態(tài)，重構(gòu)出時(shí)間分辨率的流場(chǎng)。本方法首先在較為簡(jiǎn)單的周期性實(shí)驗(yàn)數(shù)據(jù)中進(jìn)行驗(yàn)證，并且為探索該方法的適用邊界，將該方法進(jìn)一步應(yīng)用至非周期的復(fù)雜流場(chǎng)（不同直徑的雙圓柱尾渦流場(chǎng)）中。研究結(jié)果表明，該方法可以對(duì)周期性流場(chǎng)實(shí)現(xiàn)較好的重構(gòu)，而對(duì)于非周期流場(chǎng)則表現(xiàn)較差。

1 流場(chǎng)重構(gòu)原理

1.1 壓縮感知概述

壓縮感知理論[9-10]由Candès、陶哲軒及Donoho等提出。該理論指出，只要信號(hào)在某個(gè)變換域是稀疏的，那么就可以用一個(gè)與變換基不相關(guān)的觀測(cè)矩陣將變換所得的高維信號(hào)投影到一個(gè)低維空間上，然后通過求解一個(gè)優(yōu)化問題就可以從這些少量的投影中以高概率重構(gòu)出原信號(hào)[11]。由于壓縮感知方法基于信號(hào)的隨機(jī)采樣，觀測(cè)矩陣 ?定義了在高時(shí)間分辨率全流場(chǎng)F上的隨機(jī)采樣規(guī)則，基于這樣的隨機(jī)采樣規(guī)則，亞采樣的PIV數(shù)據(jù)G得以產(chǎn)生，可以定義為：

如式（3）所示，對(duì)隨機(jī)采樣的PIV數(shù)據(jù)G進(jìn)行本征正交分解，可以獲得空間模態(tài) π和相應(yīng)于數(shù)據(jù)G的亞采樣時(shí)間系數(shù)Asa。在本文中，由于時(shí)間系數(shù)在離散余弦變換的變換空間中是稀疏的[12]，因此可將離散余弦變換（DCT）作為所選取的變換域。通過壓縮感知，可以從亞采樣時(shí)間系數(shù)Asa中重構(gòu)出系數(shù)A，重構(gòu)過程求解最優(yōu)化問題：

其中，ψ為離散余弦變換，S為時(shí)間分辨率系數(shù)A在離散余弦變換中的系數(shù)。但是，求解上述問題是一個(gè)NP難問題[13]。Chen、Donoho和Saunders[14]指出，求解一個(gè)更加簡(jiǎn)單的L1優(yōu)化問題會(huì)產(chǎn)生同等的解，使得原問題簡(jiǎn)化成一個(gè)線性規(guī)劃問題。本文中，選擇典型的基追蹤算法，使用Matlab中的linprog函數(shù)進(jìn)行求解。結(jié)合所得到的空間模態(tài)，時(shí)間分辨率的流場(chǎng)重建完成。更多細(xì)節(jié)可見Baraniuk等[13]、Zhao等[15]、Stankovic等[16]對(duì)壓縮感知理論的論述。

1.2 本征正交分解概述

本征正交分解[17-18]是一種對(duì)流場(chǎng)降維，獲取流場(chǎng)中主要流動(dòng)成分的方法，已經(jīng)得到了廣泛應(yīng)用。記g為瞬態(tài)流場(chǎng)數(shù)據(jù)集G在某時(shí)刻的數(shù)據(jù)，g可以表示為：

其中：表 示平均流場(chǎng)， πi表示空間模態(tài)，w表示模態(tài)總數(shù)，aisa表示相應(yīng)于某一模態(tài)的時(shí)間系數(shù)。由于在該方法中，數(shù)據(jù)G為PIV采樣數(shù)據(jù)的集合，因此aisa所記錄的時(shí)刻點(diǎn)相應(yīng)于數(shù)據(jù)G的采樣時(shí)間。

時(shí)間系數(shù)asa可通過奇異值分解計(jì)算得到：

其中，C表示對(duì)流動(dòng)數(shù)據(jù)內(nèi)積 (G,G′)進(jìn)行空間互相關(guān)計(jì)算得到的矩陣。將數(shù)據(jù)G映射到式（4）計(jì)算所得的時(shí)間系數(shù)asa上并進(jìn)行歸一化即可得到空間模態(tài)π 。特征值 λ 是逐漸遞減的序列，表示模態(tài) π所占有的能量。模態(tài)所占有的能量比例通常被記為

1.3 流場(chǎng)重構(gòu)流程

如圖1所示，首先對(duì)低時(shí)間分辨率隨機(jī)采樣的PIV數(shù)據(jù)G進(jìn)行本征正交分解，得到空間模態(tài)及相應(yīng)的時(shí)間系數(shù)，由于該時(shí)間系數(shù)與PIV數(shù)據(jù)G遵循相同的采樣規(guī)則，因此式（1）中的采樣矩陣即為時(shí)間系數(shù)的觀測(cè)矩陣，選取離散余弦變換作為稀疏基，從而重構(gòu)得到高時(shí)間分辨率的模態(tài)系數(shù)，結(jié)合前面所得到的空間模態(tài)，即可完成流場(chǎng)重構(gòu)。

圖1 流場(chǎng)重構(gòu)流程Fig. 1 Flowchat of the flow field reconstruction method

2 算例驗(yàn)證與分析

2.1 振蕩器流場(chǎng)的驗(yàn)證與分析

振蕩器無需任何移動(dòng)部件即可產(chǎn)生周期性的射流。在流動(dòng)控制、換熱等應(yīng)用領(lǐng)域中，這種振蕩器射流具有良好的魯棒性[19]，但其振蕩頻率遠(yuǎn)高于PIV采樣頻率。針對(duì)這個(gè)問題，該時(shí)間分辨率流場(chǎng)重構(gòu)方法首先應(yīng)用于振蕩器實(shí)驗(yàn)數(shù)據(jù)，來驗(yàn)證該方法的可靠性。PIV拍攝區(qū)域如圖2所示。為驗(yàn)證重構(gòu)準(zhǔn)確性，使用TR-PIV獲取真實(shí)流場(chǎng)信息，有關(guān)該實(shí)驗(yàn)的更多詳細(xì)信息可見Wen等[20]的論文。在本文中，PIV采樣數(shù)據(jù)為間隔采樣時(shí)間超過0.444 s的隨機(jī)采樣，振蕩器頻率為2.687 Hz，而重構(gòu)所得的高時(shí)間分辨率全流場(chǎng)數(shù)據(jù)采樣頻率為27 Hz。

圖2 振蕩器PIV拍攝區(qū)域示意圖Fig. 2 Schematic of the PIV measurement region of the oscillator flow field

為了選擇適當(dāng)數(shù)量的模態(tài)，對(duì)PIV數(shù)據(jù)的POD模態(tài)累積能量譜進(jìn)行分析，由圖3所示，可以發(fā)現(xiàn)，在一共7208個(gè)模態(tài)中，隨著模態(tài)數(shù)的增加，每一模態(tài)所占據(jù)的能量比例逐漸降低，前7個(gè)模態(tài)占據(jù)了振蕩器射流約72%的能量，而前100個(gè)模態(tài)占據(jù)了約88%的能量，因此，在下面的分析中，首先選取前7個(gè)模態(tài)用于重構(gòu)，繼而選擇前100個(gè)模態(tài)用于對(duì)比。

圖3 振蕩器流場(chǎng)POD模態(tài)能量譜Fig. 3 Energy spectrum of the POD modes of the oscillator flow field

圖4對(duì)比了該模態(tài)系數(shù)在時(shí)間域上的真實(shí)值及重構(gòu)值。低頻采樣數(shù)據(jù)G所得的模態(tài)系數(shù)用黑色菱形標(biāo)出，可以看出該重構(gòu)得到的模態(tài)系數(shù)與真實(shí)值有著較好的擬合。盡管相應(yīng)于其他模態(tài)的系數(shù)這里未給出，但重構(gòu)效果基本類似，不再贅述。

基于圖4重構(gòu)得到的模態(tài)系數(shù)，結(jié)合POD分解所得空間模態(tài)，完成高時(shí)間分辨率流場(chǎng)重構(gòu)。隨機(jī)選取某周期中未被采樣的三張快照觀察其重構(gòu)效果（圖5）?？梢姡?階模態(tài)重構(gòu)流場(chǎng)準(zhǔn)確獲取了真實(shí)流場(chǎng)的主模態(tài)相關(guān)信息，但由于所選模態(tài)較少，部分細(xì)節(jié)信息有所丟失。針對(duì)該問題，進(jìn)一步增加重構(gòu)階數(shù)至100階。當(dāng)使用100階模態(tài)進(jìn)行重構(gòu)時(shí)，重構(gòu)流場(chǎng)的細(xì)節(jié)信息更為豐富。

圖4 振蕩器流場(chǎng)時(shí)間域上重構(gòu)模態(tài)系數(shù)與真實(shí)值對(duì)比Fig. 4 Comparison of reconstructed POD coefficients with ground truth data in the temporal domain for the oscillator flow field

圖5 振蕩器流場(chǎng)非采樣點(diǎn)重構(gòu)與真實(shí)值對(duì)比Fig. 5 Comparison of the reconstructed flow field with the ground truth data at a non-sampling point for the flow oscillator

為量化重構(gòu)準(zhǔn)確率，誤差由Ruscher等[21]的分析定義為：

其中，uriecon表示在一張PIV快照中空間i位置的重構(gòu)值，utiruth表示空間i位置的真實(shí)參考值，和分別表示真實(shí)PIV數(shù)據(jù)F在時(shí)間域和空間域的最大值和最小值。最終每張快照的誤差計(jì)算值對(duì)上式所有所得ei進(jìn)行平均，繼而排除采樣點(diǎn)位置處的快照數(shù)據(jù)，對(duì)所得每張快照的誤差進(jìn)行平均。所得誤差曲線如圖6所示。當(dāng)選取模態(tài)階數(shù)由3階升至20階時(shí)，重構(gòu)誤差從2.8%降至2.7%，而隨著重構(gòu)階數(shù)進(jìn)一步增多至100階時(shí)，誤差稍微增加。經(jīng)過分析，認(rèn)為該誤差的增加是由高階模態(tài)包含更多噪聲所導(dǎo)致的。

圖6 振蕩器流場(chǎng)重構(gòu)誤差隨選取模態(tài)數(shù)量的變化Fig. 6 Error incurred in data reconstruction vs the number of POD modes employed of oscillator flow field

2.2 雙圓柱流場(chǎng)的驗(yàn)證與分析

為了探索提出的方法在非周期流場(chǎng)中的重構(gòu)效果，對(duì)不同大小的雙圓柱尾渦流場(chǎng)原始PIV數(shù)據(jù)進(jìn)行了處理分析（圖7）。實(shí)驗(yàn)在循環(huán)水槽[22]中進(jìn)行，實(shí)驗(yàn)的具體細(xì)節(jié)參見文獻(xiàn)[23]。圓柱采樣頻率為250 Hz，因此本文用于重構(gòu)的數(shù)據(jù)間隔采樣大于0.048 s。

圖7 雙圓柱PIV拍攝區(qū)域示意圖Fig. 7 Schematic of the PIV measurement region of double-cylinder wake flow

圖8為 POD模態(tài)累積能量譜，可以發(fā)現(xiàn)，第一個(gè)模態(tài)占據(jù)了超過80%的能量，前7個(gè)模態(tài)占據(jù)了約92%的能量，因此，在本案例中，仍然選用前7個(gè)模態(tài)進(jìn)行重構(gòu)。

圖8 雙圓柱流場(chǎng)POD模態(tài)能量譜Fig. 8 Energy spectrum of the POD modes of the double-cylinder wake flow

相應(yīng)于第一模態(tài)，圖9對(duì)比了該模態(tài)系數(shù)的真實(shí)值及重構(gòu)值，如圖9（a）所示，由低頻采樣數(shù)據(jù)G所得的模態(tài)系數(shù)用藍(lán)色五角星標(biāo)出。在非采樣點(diǎn)，重構(gòu)得到的模態(tài)系數(shù)偏離于真實(shí)參考值，存在較多的高頻噪聲，而結(jié)合圖9（b）所示，盡管在非采樣點(diǎn)重構(gòu)值偏離于真實(shí)值，但也是在真實(shí)值附近波動(dòng)。這是由于L1范數(shù)無法區(qū)分稀疏系數(shù)尺度的位置，所以存在低尺度能量搬移到了高尺度的現(xiàn)象，從而容易出現(xiàn)一些人工效應(yīng)。即，如圖9（a）所示的一維信號(hào)高頻振蕩現(xiàn)象。整體上來說，重構(gòu)信號(hào)在歐氏距離上仍然逼近原信號(hào)。相應(yīng)于其他模態(tài)的系數(shù)未在圖中示出，但重構(gòu)效果基本類似，所以不再贅述。

圖9 雙圓柱流場(chǎng)時(shí)間域上重構(gòu)模態(tài)系數(shù)與真實(shí)值對(duì)比Fig. 9 Comparison of the reconstructed POD coefficients with the ground truth data in the temporal domain for the double-cylinder wake flow

隨機(jī)選取1張?jiān)诜遣蓸狱c(diǎn)的重構(gòu)快照，見圖10。可以發(fā)現(xiàn)雖然壓縮感知方法可以重構(gòu)流場(chǎng)主要特征，但是難以捕捉高頻流場(chǎng)細(xì)節(jié)。相比于周期性的振蕩器射流流場(chǎng)，誤差顯著增大。這說明，該方法對(duì)非周期性流場(chǎng)重構(gòu)適應(yīng)性較低。

圖10 雙圓柱流場(chǎng)非采樣點(diǎn)7階模態(tài)重構(gòu)流場(chǎng)與真實(shí)流場(chǎng)對(duì)比Fig. 10 Comparison of the reconstrcucted flow field using leading 7 POD modes with the ground truth data at a nonsampling point for the double-cylinder wake flow

3 結(jié) 論

研究提出了一種基于壓縮感知和POD的PIV流場(chǎng)高時(shí)間分辨率重構(gòu)方法，對(duì)低頻隨機(jī)采樣數(shù)據(jù)進(jìn)行本征正交分解，然后對(duì)時(shí)間系數(shù)進(jìn)行壓縮感知重構(gòu)，結(jié)合分解所得的空間模態(tài)，重構(gòu)出時(shí)間分辨率的流場(chǎng)。利用該方法處理了周期性的亞采樣振蕩器射流實(shí)驗(yàn)流場(chǎng)和非周期性的雙圓柱實(shí)驗(yàn)流場(chǎng)。結(jié)果表明，對(duì)于周期性流場(chǎng)，該重構(gòu)方法的重構(gòu)誤差在3%以下，能夠準(zhǔn)確恢復(fù)時(shí)間分辨率的流場(chǎng)。而對(duì)于非周期性復(fù)雜流場(chǎng)，則出現(xiàn)較大的高頻噪聲。

該方法為一種純數(shù)據(jù)驅(qū)動(dòng)的方法，可以實(shí)現(xiàn)任意選擇的高階或低階重構(gòu)，而且無需任何接觸式的傳感器，不會(huì)擾亂流場(chǎng)，且計(jì)算量較低。該方法不依賴于高時(shí)間分辨率的PIV測(cè)量，利用常用的PIV測(cè)量裝置即可實(shí)現(xiàn)，且可實(shí)現(xiàn)12倍的頻率放大。因此，該方法對(duì)于周期性流場(chǎng)的時(shí)間分辨率重構(gòu)是一種行之有效的方案。但目前該方法對(duì)非周期流場(chǎng)的重構(gòu)適應(yīng)性較低，需要做進(jìn)一步研究。