基于卡爾曼濾波的翼型表面壓力實時重構(gòu)方法

劉余丹，周楷文，劉應(yīng)征，溫 新，*

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

0 引 言

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

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

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

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

1 重構(gòu)方法

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

1.1 本征正交分解

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

圖1 重構(gòu)方法流程Fig.1 Outline of reconstruction method

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

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

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

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

1.2 壓縮感知

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

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

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

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

1.3 卡爾曼濾波器

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

其中：x表示系統(tǒng)狀態(tài)；A表示系統(tǒng)的線性狀態(tài)轉(zhuǎn)換矩陣；w表示協(xié)方差為Q的高斯分布的系統(tǒng)誤差?？柭鼮V波器中的測量過程可表示為：

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

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

矯正步為：

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

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

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

本文中，

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

其中：上標(biāo)+表示Moore-Penrose偽逆。

2 數(shù)據(jù)準(zhǔn)備

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

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

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

3 方法驗證

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

圖5 前4階POD空間模態(tài)及其對應(yīng)的時間系數(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)對應(yīng)的時間系數(shù)如圖6所示，同樣計算所得時間系數(shù)嚴(yán)重扭曲，含有大量噪聲。直接通過這些含有大量噪聲的時間系數(shù)對流場進(jìn)行重構(gòu)，重構(gòu)所得結(jié)果如圖7所示，與真實翼型表面壓力分布情況存在較大不符，無法滿足對翼型表面壓力進(jìn)行全場感知的要求。

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

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

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

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

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

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

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

為量化卡爾曼濾波器對升力重構(gòu)準(zhǔn)確性的提升效果，定義升力的重構(gòu)誤差為：

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

4 結(jié) 論

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

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

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

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