彈性高超聲速滑翔飛行器的狀態/參數聯合估計

陳爾康，荊武興，高長生

哈爾濱工業大學 航天學院，哈爾濱 150001

高超聲速滑翔飛行器一般指在大氣層內以馬赫數大于5的高速[1]作無動力滑翔的飛行器。為在大氣層內高速飛行，高超聲速滑翔飛行器一般采用細長體外形[2]并使用碳纖維等輕質復合材料，結構質量較小，因此彈性振動模態固有頻率較低，彈性振動與剛體運動間存在不可忽略的耦合效應[3-4]。此外，在現代飛行器上廣泛應用的高增益數字控制系統具有較大帶寬，與高超聲速滑翔飛行器的彈性振動模態固有頻率接近，又會引起彈性振動和控制系統間相互耦合的氣動伺服彈性問題[5]：導彈的結構彈性振動信號與剛體運動信號通過傳感器送至控制系統，經處理后驅動舵面偏轉，而舵面偏轉所產生的力又會激勵結構的彈性振動。文獻[6]研究了彈性導彈的氣動伺服彈性問題，發現彈性振動對慣導的測量信號有很大影響，慣導位置選取不當可能導致失穩。綜上，在高超聲速滑翔飛行器的總體設計和控制系統設計中必須考慮其結構彈性[7]。

但是傳統的控制律都是基于剛體模型設計的。為抑制彈性振動模態對控制系統的不利影響，可在控制回路中加入陷波器，或通過反饋撓度等測量信號進行主動振動抑制控制[8]。這些方法都需要精確已知彈性振動的模態信息，但這些信息很難準確計算或直接測量。在過去的幾十年中，自適應控制[9]等控制方法被用于彈性高超聲速飛行器的控制器設計，有效提高了控制性能。這些方法將彈性模態作為擾動處理，雖然能夠保證穩定性，但對控制性能有不利影響。因此，為實現對彈性高超聲速飛行器的高精度控制，需要對彈性狀態進行估計，并將其用于模型預測控制[10]、滑模控制[11]等先進控制方法。除狀態以外，彈性振動模態的固有頻率和振型等參數對飛行器的動態特性有很大影響，在狀態估計、控制器設計和故障診斷等應用中十分重要。但這些參數很難在地面準確測得，且受飛行器狀態和飛行環境影響而并非常值[12]，因此對其進行實時估計同樣是十分必要的。針對這一問題，文獻[13-14]對參數進行估計，并將估計結果用于高超聲速飛行器的自適應控制。但這類方法只對參數進行估計，無法實現狀態和參數的聯合估計。因此，開展對彈性高超聲速滑翔飛行器狀態/參數聯合估計方法[15]的研究是十分必要的。

作為一種基于模型的遞歸濾波器，擴展卡爾曼濾波(Extended Kalman Filter,EKF)利用線性化技術對卡爾曼濾波進行擴展以用于非線性系統，現已廣泛用于狀態/參數估計、故障診斷和組合導航[16]等領域。但是EKF難以處理約束和參數不確定性，而這正是滾動時域估計(Moving Horizon Estimation,MHE)的優勢。

MHE利用一段時間內的測量值，將估計問題轉化為優化問題進行求解，因此能夠處理帶有約束的狀態與參數聯合估計問題，且具有更好的魯棒性和估計精度[17-18]。MHE的優化目標函數由過程噪聲、量測噪聲和到達代價組成[19]。到達代價描述了歷史數據對估計的影響，是保證MHE穩定性的關鍵[20]。傳統的到達代價計算方法從概率和統計的角度出發，采用估計協方差矩陣計算到達代價[21]，但這類方法并未充分利用滾動時域估計的結果，在非線性較強的情況下存在估計精度較差和結果發散的問題。

事實上，估計精度不僅與估計方法有關，還與傳感器布置有關。對于高超聲速飛行器常用的速率陀螺，其量測信號中同時包含剛體分量和彈性分量，可以在不增加額外傳感器的情況下用來估計彈性狀態和參數。但陀螺位置會影響量測信號中彈性分量的大小，進而影響系統可觀性和估計精度[22]。因此還需要對速率陀螺的布置方案展開研究。

為了實現彈性高超聲速滑翔飛行器的狀態和參數聯合估計，本文提出了一種傳感器布置策略和一種利用正交三角(QR)分解更新到達代價的滾動時域估計算法(Moving Horizon Estimation with arrival cost updated by QR decomposition,MHE-QR)。通過對陀螺位置的優化提高了系統可觀性，并在最優布置方案下利用MHE-QR算法對狀態和參數進行聯合估計。MHE-QR算法將到達代價的計算問題轉化為最小二乘問題并利用QR分解進行求解，能夠保證算法穩定性，有效提高了估計精度和計算速度。

1 彈性高超聲速滑翔飛行器的數學模型

本文研究的彈性高超聲速滑翔飛行器具有如圖1所示的雙錐體外形，其幾何參數如表1所示。

由于高超聲速滑翔飛行器結構彈性的主要來源是縱向彎曲，因此本文使用的縱向動力學模型為[5]

(1)

式中：m為飛行器質量；V為飛行器速度大小；θ為彈道傾角；L為升力；g為重力加速度；?為俯仰角；q為俯仰角速度；M為俯仰力矩；I為轉動慣量；α為迎角；ηi、ζi、ωi和Ni分別為第i階彈性模態的廣義坐標、阻尼比、固有頻率和廣義力。

圖1 彈性高超聲速滑翔飛行器幾何外形Fig.1 Flexible hypersonic glide vehicle geometry

表1 幾何參數Table 1 Geometry parameters

為估計飛行器狀態、彈性模態固有頻率和振型信息，彈性高速飛行器的傳感器輸出為4個速率陀螺的角速度信號qm,j[12]：

(2)

(3)

式中：w(t)為過程噪聲；v(t)為量測噪聲；函數f描述了系統動力學；函數h描述了系統輸出。

由于MHE處理的是離散系統，因此利用多重打靶法將系統(3)離散:

(4)

2 傳感器布置策略

傳感器布置可視為一個優化問題，其中優化指標的選擇較為關鍵。由陀螺的量測模型(2)可知，陀螺測得的角速度由剛體分量和彈性分量組成。當陀螺位于振型斜率為零處時，其測得的角速度中不含彈性分量，這對彈性狀態和參數的估計是非常不利的。較好的傳感器布置方案應使得測量信號中包含較多的彈性分量以減小估計誤差。因此應選擇能夠反映系統可觀測度的優化指標。傳統的可觀性判據只能判斷狀態是否可觀，不適合作為優化指標。在這方面，Popov-Belevitch-Hautus(PBH)判據的改進形式能夠反映系統可觀測度，可選為性能指標[22-23]。該性能指標由PBH判據導出，PBH判據的具體內容如下所述。

考慮式(5)所示的系統:

(5)

若對于Ae的每一個特征值λi，式(6)均成立，則該系統是可觀的。

rank(PBH(λi,Ae,Ce))=n

(6)

式中：n為Ae的維數，PBH(λi,Ae,Ce)為式(7)所示的PBH矩陣：

(7)

式中：In為n維單位矩陣。

PBH矩陣的每個條件數能夠描述其對應模態的可觀測度。每個條件數越大，其對應的模態就越接近不可觀[21]。因此，優化指標J選為所有PBH矩陣條件數中的最大值，即

J=max(cond(PBH(λi,Ae,Ce)))i=1,2,…,n

(8)

式中：Ae和Ce可由式(1)～式(3)經線性化得到。

在僅考慮一階彈性振動模態的情況下，傳感器位置與性能指標間的關系如圖2所示。當陀螺位于飛行器中部時，性能指標最大，可觀測度最差。此時陀螺測量信號中幾乎不含彈性分量，對狀態/參數的估計非常不利。而布置在飛行器頭部或尾部時，陀螺測量信號中彈性分量較大，對狀態/參數估計較為有利。因此上述優化指標的選擇是合理的。

在考慮三階彈性模態的情況下，單傳感器位置與性能指標間的關系如圖3所示。受二階和三階彈性模態的影響，一些位置的可觀測度變差，傳感器位置與性能指標間的關系變的復雜。

圖2 單傳感器位置與性能指標間的關系(一階)Fig.2 Relationship between single sensor position and performance index (one order)

圖3 單傳感器位置與性能指標間的關系(三階)Fig.3 Relationship between single sensor position and performance index(three orders)

因此，多傳感器的布置方案需要通過優化確定。最終，多傳感器布置問題被轉化為如下優化問題：

s.t.xmin

xm,i+1-xm,i>lmin

(9)

式中：xmin和xmax分別為陀螺到飛行器頭部的最小和最大距離；lmin為相鄰兩個陀螺間的最小距離。本文中xmin、xmax和lmin的取值如表2所示。

利用序列二次規劃算法求解優化問題(9)即可得到最優傳感器布置方案，如表3所示。

表2 傳感器布置參數Table 2 Sensor placement parameters

表3 傳感器布置方案Table 3 Sensor placement scheme

3 滾動時域估計

3.1 問題描述

在本文中，為表示方便，采用如下記法：

(10)

(11)

式中：a為m維向量；W為m×m維矩陣。

記滾動窗口長度為N。MHE利用當前時刻之前N個時間點的數據，將狀態/參數估計問題轉化為優化問題：

(12)

s.t.

(13)

式中：wk為過程噪聲；Q為過程噪聲協方差矩陣，vk為量測噪聲；R為量測噪聲協方差矩陣。ΘT-N(xT-N)為到達代價，到達代價概括了滾動窗口之前的數據，對MHE的性能和穩定性有很大的影響。

3.2 到達代價更新算法

根據MHE問題的定義和前向動態規劃原理，到達代價的定義為[20]：

(14)

s.t.

xk+1=F(xk,uk)+wkk=0,1,…,T-N-1

(15)

由式(14)和式(15)可知，到達代價的計算較為復雜，對于約束非線性系統更是無法得到解析表達式[20]。因此一般采用二次函數估算到達代價，即

(16)

PT-N=Q+GPT-N-1GT-

GPT-N-1C(R+CPT-N-1CT)CPT-N-1GT

(17)

(18)

在本文中，利用式(17)更新到達代價的MHE算法稱為利用EKF更新到達代價的滾動時域估計方法(Moving Horizon Estimation with arrival cost updated by EKF,MHE-EKF)。MHE-EKF算法未能充分利用MHE的估計結果，難以進一步提高MHE的估計精度和計算速度。對此，本文在確定性框架下給出一種利用QR分解更新到達代價的算法。

由式(14)可知，T+1時刻的到達代價為

(19)

根據前向動態規劃原理和式(14)，式(19)可改寫為

(20)

s.t.

(21)

利用近似到達代價代替到達代價，將式(16)代入式(19)可得

(22)

為表示方便，對權重矩陣作Cholesky分解可得

(23)

(24)

(25)

式中：

(26)

將式(25)代入式(24)可得

(27)

至此，到達代價的更新問題被轉化為如下最小二乘問題：

(28)

式中：

(29)

為保證數值穩定性和提高計算速度，使用QR分解求解這一最小二乘問題。對A作QR分解可得

(30)

式中：E為正交矩陣；F1和F2為上三角矩陣。設

(31)

由式(29)～式(31)可得

(32)

顯然可得此最小二乘問題的解為

(33)

將式(33)代入式(28)可得

(34)

由式(34)可得到達代價更新的遞推方程為

(35)

基于QR分解的到達代價更新算法根據上一時刻MHE的估計結果計算到達代價，且能夠保證數值穩定性，能夠提高MHE的估計精度和計算速度。此外，該到達代價更新算法還具有如下性質。

(36)

(37)

注2該到達代價更新算法給出的到達代價是有界的，說明歷史數據對估計結果的影響是有限度的，不會無限增大。即對于任意向量v，有

(38)

式中：λmax為矩陣Q-1最大的特征值。

由式(29)可得

(39)

由式(30)可得

(40)

比較式(39)和式(40)等號右側矩陣第2行第2列元素可得

(41)

由式(35)和式(41)可知，對于任意向量v，有

(42)

注3該到達代價更新算法能夠保證MHE的穩定性。

根據式(20)，該到達代價算法具有如下性質：

(43)

反復應用式(43)可得

(44)

式(44)表明該達到代價更新算法滿足文獻[20] 中的條件C2：

(45)

因此根據文獻[20]中的命題3.4，該到達代價算法能夠保證MHE的穩定性。

3.3 MHE的數值解法

在3.2節提出的到達代價更新算法的基礎上，利用序列二次規劃算法求解式(46)所示的優化問題即可實現對狀態和參數的估計，形成MHE-QR算法。

(46)

s.t.

(47)

MHE算法的流程如圖4所示。

圖4 滾動時域估計算法流程圖Fig.4 Flow chart of MHE

4 仿真分析

在仿真中，窗口長度設為N=15，采樣周期為Δt=0.05 s，量測噪聲方差矩陣為R=diag(2.5×10-5,2.5×10-5,2.5×10-5,2.5×10-5)，系統噪聲方差矩陣為Q=diag(1×10-6,1× 10-2,1×10-6,2.5×10-3)。對一階彈性模態固有頻率和振型斜率加入了正弦擾動信號。為了驗證不同輸入下算法的性能，系統輸入信號如圖5 所示。

為驗證傳感器布置策略，并分析MHE-QR、MHE-EKF和EKF 3種算法的性能，分別在表4所示的仿真場景下進行仿真。

在相同條件下分別使用MHE-QR(場景1)、MHE-EKF(場景2)和EKF(場景3)3種算法進行蒙特卡羅仿真，均方根誤差(Root Mean Square Error,RMSE)如圖6和圖7所示。

圖5 輸入信號Fig.5 Input signal

表4 仿真場景Table 4 Simulation cases

圖6 不同算法下?、q和η1的均方根誤差Fig.6 RMSE of ?、qandη1 with different algorithms

圖7 不同算法下和的均方根誤差 different algorithms

同樣使用MHE-QR算法，分別在傳感器最優布置(場景1)、傳感器非最優布置(場景4)和不估計振型斜率(場景5)的情況下進行蒙特卡羅仿真，均方根誤差如圖8和圖9所示。

由圖8和圖9可知，由于傳感器并非最優布置，場景4的RMSE存在較多的突然增大的情況，估計誤差大于采用最優傳感器最優布置的場景1。而在15～20 s的時間段內，由于系統狀態已經收斂至穩態值，雖然場景1的均方根誤差大于場景4和場景5，但場景4和場景5均存在均方根誤差突然增大的情況，這對實際應用非常不利。因此，仿真結果表明傳感器非最優布置對估計精度有不利影響。此外，受不對振型斜率進行估計的影響，場景5中彈性狀態和彈性模態固有頻率的RMSE均顯著大于場景1和場景4，說明振型斜率的偏差對估計結果影響較大。因此，傳感器布置策略和包含振型斜率的狀態/參數聯合估計能夠有效提高估計精度。

圖8 不同傳感器布置下?、q和η1的均方根誤差Fig.8 RMSE of ?、q and η1 with different sensor placements

圖9 不同傳感器布置下的均方根誤差Fig.9 RMSE of with different sensor placements

對于估計算法的實際應用，計算耗時是非常重要的。表5給出了不同場景下的計算耗時，仿真在Windows 10系統(CPU為i5-7400，主頻為3.00 GHz)中MATLAB R2017a環境下完成。其中EKF(場景3)的耗時明顯短于其他算法，而MHE-QR(場景1、4、5)的平均耗時短于MHE-EKF。雖然同樣使用MHE-QR算法，但傳感器非最優布置的場景4的最大耗時是所有場景中最大的，說明傳感器最優布置能夠縮短計算耗時。由于估計狀態數量最小，場景5的平均耗時和最大耗時都是使用MHE-QR算法的所有場景中最小的，但其估計精度較差，難以實際應用。而場景1則在估計精度和計算耗時之間做到了較好的平衡。得益于傳感器最優布置和基于QR分解的到達代價更新策略，場景1的估計精度是最高的，但其平均耗時和最大耗時都僅長于使用EKF的場景2和不估計振型斜率的場景5，且其最大耗時并未超出采樣周期，表明其具備實時應用潛力。因此，本文提出的傳感器布置策略和到達代價更新策略不僅能夠提高估計精度，也能夠提高計算速度。

表5 計算耗時Table 5 CPU time consumption

5 結 論

1)為提高彈性高超聲速滑翔飛行器的狀態/參數估計精度，本文提出了一種傳感器布置優化方法。該優化方法選擇PBH矩陣條件數的最大值為優化指標。仿真結果表明，傳感器布置優化方法能夠有效提高狀態/參數的估計精度。

2)針對彈性高超聲速滑翔飛行器的狀態/參數聯合估計問題，本文提出一種新型的滾動時域估計算法——MHE-QR算法。該算法將到達代價的估計問題轉化為最小二乘問題，并給出一種新型的利用QR分解更新到達代價的算法。MHE-QR算法限制了到達代價過度增長，能夠保證滾動時域估計算法的穩定性。仿真結果表明，MHE-QR算法的估計精度高于EKF和MHE-EKF算法。此外，MHE-QR的計算耗時小于采樣周期，具備實時應用的潛力。

后續需要研究滾動時域估計算法的數值求解方法，以進一步提高計算速度，滿足實時應用的要求。