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

(1.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240；2.空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽 621000)

0 引 言

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

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

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

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

1 重構(gòu)方法

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

1.1 本征正交分解

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

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

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

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

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

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

1.2 壓縮感知

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

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

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

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

1.3 卡爾曼濾波器

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

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

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

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

矯正步為：

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

1.4 動(dòng)態(tài)模態(tài)分解

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

本文中，

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

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

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

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

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

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

3 方法驗(yàn)證

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

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

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

圖6 壓縮感知所得前4階模態(tài)時(shí)間系數(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

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

圖8 經(jīng)卡爾曼濾波器處理后前4階模態(tài)時(shí)間系數(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ù)積分可計(jì)算其升力l。經(jīng)壓縮感知所重構(gòu)升力和卡爾曼濾波器處理所重構(gòu)升力與真實(shí)升力變化曲線的對(duì)比如圖10所示。可以看出，經(jīng)壓縮感知所重構(gòu)升力在總體變化趨勢(shì)上與真實(shí)情況基本一致，但存在較嚴(yán)重的噪聲問題，無法滿足飛行器飛行過程中精準(zhǔn)感知的要求。在采用卡爾曼濾波器進(jìn)行處理之后，重構(gòu)結(jié)果與真實(shí)值之間的差異顯著降低，重構(gòu)結(jié)果基本滿足對(duì)升力進(jìn)行準(zhǔn)確感知的要求。

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

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

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

4 結(jié) 論

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

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

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

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