劉余丹，周楷文，劉應征，溫 新，*

(1.上海交通大學 機械與動力工程學院，上海 200240；2.空氣動力學國家重點實驗室，綿陽 621000)

0 引 言

飛行器實時控制是實現(xiàn)飛行器姿態(tài)調整、提高飛行器機動性能的有效方法之一，而飛行器環(huán)境參數(shù)的感知是實現(xiàn)其飛行控制的基礎。飛行器在飛行過程中往往面臨著復雜的非定常流動問題，但傳統(tǒng)的傳感器設置方案，如慣性測量裝置、磁力儀、風速儀等，受測量技術、安裝空間等限制，僅能感知到局部的、離散的參數(shù)，難以為飛行操縱系統(tǒng)及時提供飛行控制所需的全部信息。在此情況下，難以對飛行器的關鍵氣動參數(shù)進行精準識別，進而難以實現(xiàn)運行監(jiān)測與主動控制。因此，利用有限的觀測數(shù)據(jù)實時感知飛行器周圍流場的全貌具有重要意義。

許多學者研究了如何通過有限的測量數(shù)據(jù)重建整個流場的問題。其中，一部分學者的研究主要基于Everson和Sirovich[1]在1995年提出的一種恢復缺失數(shù)據(jù)的gappy POD方法展開。gappy POD方法可用于將風洞實驗所獲得的有限表面壓力測量數(shù)據(jù)與CFD計算數(shù)據(jù)融合，以感知流場的全貌[2]，但gappy POD所采用的求解流場重構系數(shù)的最小二乘法對噪聲較為敏感。還有一些學者基于機器學習方法進行了研究。Deng等[3]采用長短時間序列（LSTM）神經(jīng)網(wǎng)絡建立離散測量點數(shù)據(jù)與POD模態(tài)系數(shù)之間的關系，進而實現(xiàn)流場全重構。但機器學習的方法存在缺乏物理解釋性、過擬合等問題。

Candes等[4]2006年提出的壓縮感知理論指出，如果一個信號在某一特定基上具有稀疏表示，則該信號可以從相對較少的隨機采樣中得到恢復，該結論為從有限的測量數(shù)據(jù)實現(xiàn)流場實時重構指明了新的方向。此前已有部分學者嘗試將壓縮感知理論應用于流場重構。Huang[5]對壓縮感知在航空航天應用中的潛力進行了探討。Zhao等[6]提出了一種基于壓縮感知重構翼型表面壓力的精細方法，以獲得翼型表面完整的壓力分布信息以及精確的升力和力矩特性。Bright等[7]利用壓縮感知技術，通過圓柱體表面的有限離散壓力測量數(shù)據(jù)重構圓柱體周圍的流動特性。Wen等[8]開發(fā)了一種基于壓縮感知的數(shù)據(jù)挖掘方法，從低速圓柱尾流誘導的非定常壓力場的高噪聲測量中對流場數(shù)據(jù)進行重構。但基于壓縮感知的流場重構方法同樣不具魯棒性，易受噪聲影響。

本文所提出的基于壓縮感知與卡爾曼濾波器的實時流場重構方法，采用卡爾曼濾波器，針對基于壓縮感知的流場重構方法進行改進。該重構方法具體可分為線上和線下兩個部分：線下全場采樣，通過本征正交分解建立數(shù)據(jù)庫；線上稀疏采樣，通過壓縮感知優(yōu)化重構系數(shù)，然后采用基于內嵌動態(tài)模態(tài)分解的卡爾曼濾波器對所求解時間系數(shù)進行實時去噪，提高重構方法的魯棒性。

1 重構方法

基于壓縮感知與卡爾曼濾波器的實時流場重構方法流程圖如圖1所示。

1.1 本征正交分解

本征正交分解（POD）是提取離散數(shù)據(jù)主導特征的有效方法[9-10]。本文使用POD對帶有噪聲的全場測量數(shù)據(jù)Pfull進行處理，構建主導空間模態(tài) ψ 與對應時間系數(shù)afull組成的線下數(shù)據(jù)庫，該過程對應于圖1中紅色A所標記的過程。采用POD可將流場Pfull分解為以下形式：

圖1 重構方法流程Fig.1 Outline of reconstruction method

其中：k表示流場快照序號；M表示流場快照總數(shù)；表示流場快照均值； ψi表示POD主導空間模態(tài)；表示主導模態(tài)對應的POD時間系數(shù)。該分解過程可利用奇異值分解（SVD）實現(xiàn)：

其中：λ表示各POD模態(tài)能量貢獻。

以能量貢獻排序的前幾階POD模態(tài)捕獲了流場的主導流動結構，可用于構建原始流場的低階近似模型：

其中：m表示低階模型維度。構建低階模型可簡化對復雜流場的描述，同時，通過選擇主導空間模態(tài)并忽略其他不重要的模態(tài)，可在一定程度上起到濾除數(shù)據(jù)噪聲的效果。

1.2 壓縮感知

壓縮感知是一種尋找欠定線性系統(tǒng)的稀疏解的技術，已在各個領域得到廣泛應用[11-13]。壓縮感知理論指出，對于稀疏信號或可壓縮信號，可通過遠低于奈奎斯特采樣率的采樣數(shù)據(jù)實現(xiàn)其精確重構。

本文通過求解L1范數(shù)下的最優(yōu)化問題求解流場重構系數(shù)acs：

其中：Pscattered表示線上稀疏采樣數(shù)據(jù)；ψ表示線下所得POD主導模態(tài)；矩陣 ?存儲離散傳感器的位置信息。詳情參見文獻[14]。最終，可通過流場重構系數(shù)與POD主導模態(tài)的線性組合實現(xiàn)全場數(shù)據(jù)的重構：

其中：m表示用于實現(xiàn)流場重構的主導模態(tài)數(shù)量。

1.3 卡爾曼濾波器

卡爾曼濾波器提供了一種基于線性系統(tǒng)模型與測量數(shù)據(jù)對系統(tǒng)狀態(tài)進行最優(yōu)估計的方法[15-16]。本文引入卡爾曼濾波器對壓縮感知求解所得流場重構系數(shù)acs進行實時去噪，進一步提高流場重構方法的魯棒性?？柭鼮V波器中系統(tǒng)的線性模型可表示為：

其中：x表示系統(tǒng)狀態(tài)；A表示系統(tǒng)的線性狀態(tài)轉換矩陣；w表示協(xié)方差為Q的高斯分布的系統(tǒng)誤差。卡爾曼濾波器中的測量過程可表示為：

其中：y表示測量結果；H表示觀測矩陣；v表示協(xié)方差為R的高斯分布的測量誤差。

卡爾曼濾波器包含預估、矯正兩個步驟。預估步為：

矯正步為：

其中：上標 p表 示預估值；上標 c 表示校正值；F表示預估協(xié)方差矩陣；K表示卡爾曼增益。本文中，所要估測的系統(tǒng)狀態(tài)為POD時間系數(shù)a。

1.4 動態(tài)模態(tài)分解

采用動態(tài)模態(tài)分解構建卡爾曼濾波器的線性系統(tǒng)模型，對非線性流動問題進行線性化近似[17]：

本文中，

矩陣A的最佳形式可以表示為：

其中：上標+表示Moore-Penrose偽逆。

2 數(shù)據(jù)準備

采用商業(yè)軟件ANSYS模擬處于圓柱尾渦影響下的翼型表面壓力變化，CFD模型示意圖如圖2所示，選用k-ω SST作為湍流模型，來流速度設置為30 m/s，翼型表面隨機設置12個離散壓力傳感器，傳感器布置位置如圖中紅色·所示。

圖2 CFD模擬模型示意圖及離散傳感器布置位置Fig.2 Schematic diagram of CFD simulation model and deployment of discrete sensors

為模擬線下全場采樣中可能產(chǎn)生的噪聲問題，對所模擬翼型表面壓力數(shù)據(jù)加入標準差為流場參數(shù)均值10%的高斯噪聲，所得單周期無噪聲與含噪聲流場的演化過程對比如圖3所示。同時，為模擬線上稀疏采樣過程中可能產(chǎn)生的噪聲問題，對所合成離散數(shù)據(jù)加入標準差為流場參數(shù)均值10%的高斯噪聲，所得單周期無噪聲和有噪聲的離散數(shù)據(jù)隨時間變化的曲線圖如圖4所示。

3 方法驗證

采用POD對線下所獲得含有噪聲的全場數(shù)據(jù)Pfull進行處理，計算其主導空間模態(tài) ψ 與其對應的時間系數(shù)afull。前四階POD空間模態(tài)和時間系數(shù)如圖5所示。POD主導空間模態(tài)捕獲了流場的主要流動特征，與之對應的時間系數(shù)描述了其在時間上的動態(tài)演變。

圖5 前4階POD空間模態(tài)及其對應的時間系數(shù)Fig.5 The first four POD spatial modes and corresponding time coefficients

選取前十階POD模態(tài)作為壓縮感知的稀疏基，以含有噪聲的線上采樣稀疏數(shù)據(jù)Pscattered作為觀測值，在L1范數(shù)最小化的條件下求解最優(yōu)時間系數(shù)acs。壓縮感知計算所得前四階模態(tài)對應的時間系數(shù)如圖6所示，同樣計算所得時間系數(shù)嚴重扭曲，含有大量噪聲。直接通過這些含有大量噪聲的時間系數(shù)對流場進行重構，重構所得結果如圖7所示，與真實翼型表面壓力分布情況存在較大不符，無法滿足對翼型表面壓力進行全場感知的要求。

圖6 壓縮感知所得前4階模態(tài)時間系數(shù)Fig.6 The first four modal time coefficients obtained by compressed sensing

圖7 壓縮感知所重構翼型表面壓力Fig.7 Airfoil surface pressure reconstructed by compressed sensing

為了解決上述噪聲干擾問題，在壓縮感知流場重構方法的基礎上引入了卡爾曼濾波器。采用DMD對線下所得流場時間系數(shù)進行處理，以DMD所得線性近似模型x′=Ax作為卡爾曼濾波器的系統(tǒng)模型；以線上壓縮感知計算得到的含有噪聲的時間系數(shù)作為觀測值進行濾波?？柭鼮V波器中過程噪聲方差矩陣Q和 測量噪聲方差矩陣R應基于系統(tǒng)模型誤差與傳感器測量誤差的統(tǒng)計學信息設置，本文此處分別設置為0.0075Eeye(m)和 2Eeye(m)，其中Eeye(m)表示m階單位矩陣。經(jīng)卡爾曼濾波器實時去噪后所得的時間系數(shù)得到顯著平滑，如圖8所示。利用濾波后的時間系數(shù)進行重構，重構所得結果如圖9所示，與各時刻翼型表面真實壓力分布吻合良好，滿足對翼型表面壓力進行全場感知的要求。

圖8 經(jīng)卡爾曼濾波器處理后前4階模態(tài)時間系數(shù)Fig.8 The first four modal time coefficients processed by Kalman filter

圖9 經(jīng)卡爾曼濾波器所重構翼型表面壓力Fig.9 Airfoil surface pressure reconstructed by Kalman filter

基于翼型表面壓力數(shù)據(jù)積分可計算其升力l。經(jīng)壓縮感知所重構升力和卡爾曼濾波器處理所重構升力與真實升力變化曲線的對比如圖10所示?？梢钥闯觯?jīng)壓縮感知所重構升力在總體變化趨勢上與真實情況基本一致，但存在較嚴重的噪聲問題，無法滿足飛行器飛行過程中精準感知的要求。在采用卡爾曼濾波器進行處理之后，重構結果與真實值之間的差異顯著降低，重構結果基本滿足對升力進行準確感知的要求。

圖10 重構升力系數(shù)對比Fig.10 Comparison of lift coefficients reconstructed by compressed sensing and Kalman filter

為量化卡爾曼濾波器對升力重構準確性的提升效果，定義升力的重構誤差為：

其中：M表示流場快照總數(shù)；lrecon表示重構所得升力；ltruth表示真實升力。經(jīng)壓縮感知所重構流場與經(jīng)卡爾曼濾波器處理后所重構升力誤差分別為6.39%和18.15%，誤差相較下降64.79%

4 結 論

本文提出了一種魯棒的基于數(shù)據(jù)融合的實時流場重構方法，可以從有限的測量數(shù)據(jù)中實現(xiàn)全場數(shù)據(jù)的重構。通過壓縮感知和DMD線性模型內嵌的卡爾曼濾波器，顯著增強了其魯棒性，有效抑制了噪聲的影響。

以翼型模擬數(shù)據(jù)為例測試了該方法的性能，并對是否采用卡爾曼濾波器進行重構的結果進行了量化分析。對于翼型表面壓力數(shù)據(jù)，升力重構誤差從18.15%減小到6.39%，相較下降64.79%。

另外值得注意的是，卡爾曼濾波器的成功應用在很大程度上依賴于模型誤差方差矩陣及傳感器的誤差方差矩陣的準確性，而在實際應用中方差矩陣很難確定。為了對此進行改進，我們目前正在開發(fā)基于DMD線性模型內嵌的自適應卡爾曼濾波器，引入對噪聲統(tǒng)計特性的實時動態(tài)估測，解決實際應用中誤差方差矩陣難以準確估測的問題。

此外，本文目前所采用POD、DMD、卡爾曼濾波等流場分解與建模方法均為線性化方法，且DMD線性建模方法僅對較為簡單的周期性流動問題具有較好的適用性。后期考慮采用卷積神經(jīng)網(wǎng)絡、擴展DMD、集合卡爾曼濾波等非線性建模方法，進一步提高該方法的普適性。