基于期望最大化算法和求容積卡爾曼平滑器的氣動參數(shù)辨識算法

劉藝蕾, 張煒, 劉明航, 王宏偉, 左軍毅

?

基于期望最大化算法和求容積卡爾曼平滑器的氣動參數(shù)辨識算法

劉藝蕾1,2, 張煒1,3, 劉明航1,2, 王宏偉1,3, 左軍毅1,3

摘要：針對初值和噪聲統(tǒng)計特性未知情形下的飛行器系統(tǒng)辨識的問題，提出了基于期望最大化(expectation maximization, EM)和求容積卡爾曼平滑器(cubature Kalman smoother, CKS)的辨識算法。該算法用期望最大化算法對初值和噪聲的統(tǒng)計特性進行估計；用求容積卡爾曼平滑器估計狀態(tài)向量和未知參數(shù)。在期望最大化算法的求期望步驟中，所求的期望值通過求容積規(guī)則獲得，用較少的采樣點保證了估計精度；在期望最大化算法的最大化步驟中，未知量的最優(yōu)值以解析解形式給出，減小了計算量。仿真結(jié)果說明，該算法在飛行器氣動參數(shù)辨識問題中，能給出較好的辨識結(jié)果。與其他方法的對比驗證說明新算法具有辨識精度高、收斂速度快等優(yōu)點。

關鍵詞：參數(shù)辨識；期望最大化；求容積規(guī)則；求容積平滑器；飛機參數(shù)辨識

飛行器氣動力模型對于飛行器設計的很多任務環(huán)節(jié)如飛行控制系統(tǒng)的設計等都是必不可少的，其處理結(jié)果常被作為新研飛機和改型飛機最終定型的依據(jù)[1-3]。

目前，飛行器氣動力參數(shù)辨識領域已經(jīng)存在一些較為經(jīng)典的算法[2-4]。傳統(tǒng)的參數(shù)辨識算法[5-8]往往假設飛行試驗中描述飛行器運動的系統(tǒng)方程不受系統(tǒng)噪聲的干擾，然而在實際的飛行試驗中系統(tǒng)噪聲由于各種原因是不可避免的，強行忽略會產(chǎn)生較大估計誤差。針對這種情況，可以采用基于狀態(tài)增廣的濾波方法，如卡爾曼濾波(KF)[9]，擴展卡爾曼濾波(EKF)[10]，無跡卡爾曼濾波(UKF)[11]和求容積卡爾曼濾波(CKF)[13]等。相對于濾波，平滑方法因為利用更多的量測信息而獲得更高精度的估計，例如擴展卡爾曼平滑器(EKS)，正反向無跡平滑器[14](URTSS)和求容積卡爾曼平滑器(CKS)[15]等。其中，求容積卡爾曼平滑器估計精度高、計算量相對較小且適用于高維狀態(tài)量的估計問題，具有一定的優(yōu)勢。

然而，上述濾波器及平滑器在進行參數(shù)估計時需要噪聲的方差統(tǒng)計特性，而不準確的甚至錯誤的噪聲假設往往對估計結(jié)果的精度影響很大，進而影響飛行器氣動力參數(shù)的辨識。因此，需要另外的方法對噪聲協(xié)方差及狀態(tài)的初值等統(tǒng)計量進行估計。期望最大化(EM)算法是由Dempster等[16]提出的極大似然估計的一種經(jīng)典算法，其數(shù)值穩(wěn)定性好[17]并且應用廣泛。因此本文提出一種基于EM和CKS的聯(lián)合估計算法，其中CKS用于估計狀態(tài)向量和未知參數(shù)，EM算法估計狀態(tài)均值初值、狀態(tài)協(xié)方差初值和噪聲協(xié)方差等統(tǒng)計量。近期，運用EM算法結(jié)合估計算法的參數(shù)辨識技術(shù)也取得了一些進展。Bavdekar等[18]將EM算法與EKF相結(jié)合，Sch?n等[19]將粒子濾波與EM算法相結(jié)合，但仍存在精度與計算量等的問題。Yokoyama[20]在2011年提出了將URTSS與EM算法相結(jié)合的參數(shù)估計技術(shù)，并且推導了EM在此類問題中的解析形式，提高了計算精度并減小了計算量；然而URTSS解決高維問題時效果不佳且響應速度慢[20]。因此本文將CKS與EM算法相結(jié)合，有望使算法具有性能穩(wěn)定、精度較高、收斂速度快以及計算量相對較小等特點。

1問題描述

考慮以下非線性離散系統(tǒng)狀態(tài)空間模型中，對未知參數(shù)θ的估計問題：

(1)

(2)

為了對未知參數(shù)進行估計,需要對狀態(tài)向量進行擴維。狀態(tài)擴維后的狀態(tài)向量表示為xk]T,原系統(tǒng)方程可重新構(gòu)造為

(3)

2基于EM的辨識算法

將EM算法應用于上述非線性離散系統(tǒng)狀態(tài)空間模型中,將擴維狀態(tài)向量XN作為EM中的不可觀測數(shù)據(jù),量測向量ZN作為EM中的可觀測數(shù)據(jù),λ為未知統(tǒng)計量。由EM算法的原理可知,其代價函數(shù)為

(4)

在給定參數(shù)λ的條件下,系統(tǒng)狀態(tài)和觀測值得對數(shù)聯(lián)合概率分布為

(5)

(6)

式中

內(nèi)蒙古錫盟～江蘇泰州±800kV特高壓直流輸電工程天津段路徑總長196.153km，線路經(jīng)過天津市薊縣、寶坻區(qū)、武清區(qū)、西青區(qū)、靜海區(qū)、河北省廊坊市6個地市級行政區(qū)。全線以平地為主，沿線樹木較多。線路天津段跨越河道18條，蓄滯洪區(qū)5個。依據(jù)《中華人民共和國防洪法》規(guī)定：在洪泛區(qū)、蓄滯洪區(qū)內(nèi)建設非防洪建設項目，應當就洪水對建設項目可能產(chǎn)生的影響和建設項目對防洪可能產(chǎn)生的影響作出評價，編制洪水影響評價報告，提出防御措施。

(7)

(8)

(9)

(10)

(11)

(12)

(13)

式中

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

為了更加簡潔明了,本文提出的算法步驟歸納如下:

步驟5若達到收斂條件,則結(jié)束;否則返回步驟2。

3算法在飛機參數(shù)辨識中的應用

3.1ATTAS試驗飛機橫向線性模型

飛機的輸入數(shù)據(jù)取自ATTAS試驗飛機的飛行數(shù)據(jù),其試驗數(shù)據(jù)及簡化線性模型來自參考文獻[1]。ATTAS試驗飛機的橫向線性模型如下所示

(29)

(30)

式中,t為時間/s,p為滾轉(zhuǎn)角速率/(rad/s),r為偏航角速率/(rad/s),β為側(cè)滑角/(rad),δα為副翼偏轉(zhuǎn)角/(rad),δr為方向舵偏轉(zhuǎn)角/(rad)。在上述模型中,將β作為輸入變量,并假定其能夠無誤差測量。ωk(·)為過程噪聲。在飛行試驗中滾轉(zhuǎn)角速率和偏航角速率作為觀測值,其觀測方程為

(31)

(32)

3.2實驗仿真設計

仿真輸入量(副翼偏轉(zhuǎn)角δσ/(rad)、方向舵偏轉(zhuǎn)角δr/(rad)和側(cè)滑角β/(rad))輸入數(shù)據(jù)的采樣時間為15s,共有750個采樣點,采樣間隔為0.02s,其隨時間的變化如圖1所示。未知向量Θ的真值如表1所示。仿真狀態(tài)向量初值取為量測值的第一組值,即p0=9.662×10-3,r0=-9.856×10-4;初始狀態(tài)向量協(xié)方差為零,即P0=014×14;過程噪聲方差和量測噪聲方差分別取為Q=10-8I2和R=10-6I2。在上述輸入向量及參數(shù)設定的基礎上,可得到輸出(滾轉(zhuǎn)角速率p和偏航角速率r)如圖2所示。

將本文提出的算法(EM-CKS)與文獻[20]中的算法(EM-URTSS)用于ATTAS試驗飛機的橫向線性氣動模型進行辨識。算法在MATLAB環(huán)境下運行,其最大迭代次數(shù)設為i=5 000,并且將達到規(guī)定迭代次數(shù)作為唯一終止條件。

圖1　輸入量隨時間的變化

3.3實驗結(jié)果及分析

通過仿真實驗可知,EM-CKS方法比EM-URTSS算法的計算時長略短,其單步迭代計算機時分別為6.691有s和7.057s。原因在于EM-CKS方法的Sigma點個數(shù)少于EM-URTSS,且在計算Sigma點之前,能夠提前對采樣點進行計算,其權(quán)值為固定值。EM-URTSS方法在計算Sigma點時運算相對于EM-CKS稍顯復雜,均值和方差的權(quán)值則需分別計算,且每一組點的權(quán)值都不相同。

圖2　輸出量隨時間的變化

EM-CKS和EM-URTSS參數(shù)估計值如表1所示。

表1　未知參數(shù)的真值與估計值

由表中可知EM-CKS與EM-URTSS的估計精度基本一致。2種方法的待估計參數(shù)(僅選擇其中2個)隨迭代次數(shù)的收斂曲線如圖3所示。從圖3中能夠明顯看出,在參數(shù)估計結(jié)果和估計精度基本一致的情況下,EM-CKS方法的收斂速度顯著快于EM-URTSS方法,只需很少的迭代次數(shù)就可以達到收斂狀態(tài)。在2種算法均使用EM進行噪聲估計的情況下,這種現(xiàn)象說明本文算法所選用的CKS在保證估計精度的同時,還具有較好的收斂特性。

圖3　參數(shù)隨迭代次數(shù)的收斂曲線

EM-CKS和EM-URTSS方法對噪聲方差有著基本類似的估計值,與仿真真值接近,其均為

噪聲方差Q和R的跡隨迭代次數(shù)的收斂曲線如圖4所示。圖像表明EM-CKS和EM-URTSS方法都能對未知參數(shù)進行很好的估計,但是EM-CKS的收斂速度明顯高于EM-URTSS。

圖4　噪聲方差的跡隨迭代次數(shù)的收斂曲線

4結(jié)論

本文提出了基于期望最大化(EM)算法與求容積卡爾曼平滑器(CKS)的聯(lián)合估計參數(shù)辨識算法,對算法進行了詳細的推導和設計,給出了算法的計算步驟。算例驗證表明,本文提出的算法能夠同時對噪聲、狀態(tài)量及未知參數(shù)進行估計,具有較好的辨識效果;算法能夠有效地應用于飛行器氣動參數(shù)辨識問題中,通過飛機的輸入及響應對飛行器氣動模型中的未知參數(shù)進行辨識。通過與其他算法在計算性能方面進行對比,發(fā)現(xiàn)本文提出的算法辨識精度高、收斂速度快,具有比較明顯的優(yōu)勢。

參考文獻：

[1]JategaonkarRV.FlightVehicleSystemIdentification:ATimeDomainMethodology[M].Reston,VA:AIAA, 2006

[2]HamelPG,JategaonkarRV.EvolutionofFlightVehicleSystemIdentification[J].JournalofAircraft, 1996, 33(1): 9-28

[3]OwensB,BrandonJay,CroomMark,etal.OverviewofDynamicTestTechniquesforFlightDynamicsResearchatNASALaRC[C]∥25thAIAAAerodynamicMeasurementTechnologyandGroundTestingConference,FluidDynamicsandCo-LocatedConferences, 2006

[4]IliffKW.ParameterEstimationforFlightVehicles[J].JournalofGuidance,Control,andDynamics, 1989, 12(5): 609-622

[5]KutluayU.AnApplicationofEquationErrorMethodtoAerodynamicModelIdentificationAndParameterEstimationofaGlidingFlightVehicle[C]∥AIAAAtmosphericFlightMechanicsConference,Chicato,Illinois, 2009

[6]RohlfsM.IdentificationofNon-LinearDerivativeModelsFromBo105FlightTestData[J].AeronauticalJournal, 1998,102(1011): 1-8

[7]LiCW,ZouXH.MaximumLikelihoodMethodBasedonInteriorPointAlgorithmforAircraftParameterIdentification[J].JounalofAircraft, 2005, 42(5): 1355-1358

[8]ParisAC,AlaverdiO.NonlinearAerodynamicModelExtractionfromFlight-TestDatafortheS-3BViking[J].JournalofAircraft, 2005, 42(1): 26-32

[9]KalmanRE.ANewApproachtoLinearFilteringandPredictionTheory[J].TransonASMEJofBasicEngineering, 1960, 82(D):35-46

[10]JazwinskiAH.StochasticProcessesandFilteringTheory[M].NewYork:Academic, 1970: 235-237

[11]JulierSJ,UhlmannJK.UnscentedFilteringandNonlinearEstimation[J].ProceedingsoftheIEEE, 2004, 92(3): 401-422

[12]ItoK,XiongK.GaussianFiltersforNonlinearFilteringProblems[J].IEEETransonAutomaticControl, 2000,45(5): 910-927

[13]ArasaratnamI,HaykinS.CubatureKalmanFilters[J].IEEETransonAutomaticControl, 2009, 54(6): 1254-1269

[14]SarkkaS.UnscentedRauch-Tung-StriebelSmoother[J].IEEETransonAutomaticControl, 2008, 53(3): 845-849

[15]ArasaratnamI,HaykinS.CubatureKalmanSmoothers[J].Automatica, 2011, 47: 2245-2250

[16]DempsterA,LairdN,RubinD.MaximumLikelihoodfromIncompleteDataViaTheEMAlgorithm[J].JournaloftheRoyalStatisticalSociety, 1977, 39(B):1-38

[17]LangeK.AGradientAlgorithmLocallyEquivalenttotheEMAlgorithm[J].JournaloftheRoyalStatisticalSociety, 1995, 57(2), 425-437

[18]BavdekarVA,DeshpandeAP,PatwardhanSC.IdentificationofProcessandMeasurementNoiseCovarianceforStateandParameterEstimationUsingExtendedKalmanFilter[J].JournalofProcessControl, 2011,21(4): 585-601

[19]Sch?nTB,WillsA,NinnessB.SystemIdentificationofNonlinearState-SpaceModels[J].Automatica, 2011, 47(1): 39-49

[20]YokoyamaN.ParameterEstimationofAircraftDynamicsviaUnscentedSmootherwithExpectation-MaximizationAlgorithm[J].JournalofGuidance,Control,andDynamics, 2011, 34(2): 426-436

Aircraft System Identification via Cubature Kalman Smoother with Expectation Maximization Algorithm

Liu Yilei1,2, Zhang Wei1,3, Liu Minhang1,2, Wang Hongwei1,3, Zun Junyi1,3

Abstract:This paper developed a novel system identification algorithm to estimate parameter of aircraft dynamics modeled in state space. The developed method utilizes the cubature Kalman smoother to estimate the state and unknown parameters, combined with expectation-maximization algorithm, which estimates the statistics-unknown parameters, i.e., the mean and covariance of an initial state, and the covariance of both process noise and measurement noise. To reduce the computational cost with considerable accuracy decline, the cubature Kalman smoother is employed to approximate the expectation values in the expectation maximization. Further, the analytical forms of unknown statistics parameters are given in the maximization step, which makes the nonconvex numerical optimization unnecessary. Its effectiveness is demonstrated through one problem of estimating aircraft aerodynamic parameters. The result shows that the proposed algorithm is of high accuracy as well as converge faster compared with other algorithms.

Keywords:angular velocity, cost reductior, MATLAB, parameter extraction, cubature Kalman smoother, cubature Kalman filter, cubature rule, parameter identification, expectation Maximization, aircraft system identification

收稿日期：2016-04-01基金項目：國家自然科學基金(11472222、61473227)、陜西省自然科學基金(2015JM6304)、航空科學基金(20151353018)及航天技術(shù)支撐基金(2014-HT-XGD)資助

作者簡介：劉藝蕾(1990—)，女，中航工業(yè)第一飛機設計研究院助理工程師，主要從事適航技術(shù)研究。

中圖分類號：V217

文獻標志碼：A

文章編號:1000-2758(2016)04-0587-06