基于期望最大化算法和求容积卡尔曼平滑器的气动参数辨识算法

刘艺蕾, 张炜, 刘明航, 王宏伟, 左军毅



基于期望最大化算法和求容积卡尔曼平滑器的气动参数辨识算法

刘艺蕾1,2, 张炜1,3, 刘明航1,2, 王宏伟1,3, 左军毅1,3

摘要：针对初值和噪声统计特性未知情形下的飞行器系统辨识的问题，提出了基于期望最大化(expectation maximization, EM)和求容积卡尔曼平滑器(cubature Kalman smoother, CKS)的辨识算法。该算法用期望最大化算法对初值和噪声的统计特性进行估计；用求容积卡尔曼平滑器估计状态向量和未知参数。在期望最大化算法的求期望步骤中，所求的期望值通过求容积规则获得，用较少的采样点保证了估计精度；在期望最大化算法的最大化步骤中，未知量的最优值以解析解形式给出，减小了计算量。仿真结果说明，该算法在飞行器气动参数辨识问题中，能给出较好的辨识结果。与其他方法的对比验证说明新算法具有辨识精度高、收敛速度快等优点。

关键词：参数辨识；期望最大化；求容积规则；求容积平滑器；飞机参数辨识

飞行器气动力模型对于飞行器设计的很多任务环节如飞行控制系统的设计等都是必不可少的，其处理结果常被作为新研飞机和改型飞机最终定型的依据[1-3]。

目前，飞行器气动力参数辨识领域已经存在一些较为经典的算法[2-4]。传统的参数辨识算法[5-8]往往假设飞行试验中描述飞行器运动的系统方程不受系统噪声的干扰，然而在实际的飞行试验中系统噪声由于各种原因是不可避免的，强行忽略会产生较大估计误差。针对这种情况，可以采用基于状态增广的滤波方法，如卡尔曼滤波(KF)[9]，扩展卡尔曼滤波(EKF)[10]，无迹卡尔曼滤波(UKF)[11]和求容积卡尔曼滤波(CKF)[13]等。相对于滤波，平滑方法因为利用更多的量测信息而获得更高精度的估计，例如扩展卡尔曼平滑器(EKS)，正反向无迹平滑器[14](URTSS)和求容积卡尔曼平滑器(CKS)[15]等。其中，求容积卡尔曼平滑器估计精度高、计算量相对较小且适用于高维状态量的估计问题，具有一定的优势。

然而，上述滤波器及平滑器在进行参数估计时需要噪声的方差统计特性，而不准确的甚至错误的噪声假设往往对估计结果的精度影响很大，进而影响飞行器气动力参数的辨识。因此，需要另外的方法对噪声协方差及状态的初值等统计量进行估计。期望最大化(EM)算法是由Dempster等[16]提出的极大似然估计的一种经典算法，其数值稳定性好[17]并且应用广泛。因此本文提出一种基于EM和CKS的联合估计算法，其中CKS用于估计状态向量和未知参数，EM算法估计状态均值初值、状态协方差初值和噪声协方差等统计量。近期，运用EM算法结合估计算法的参数辨识技术也取得了一些进展。Bavdekar等[18]将EM算法与EKF相结合，Schön等[19]将粒子滤波与EM算法相结合，但仍存在精度与计算量等的问题。Yokoyama[20]在2011年提出了将URTSS与EM算法相结合的参数估计技术，并且推导了EM在此类问题中的解析形式，提高了计算精度并减小了计算量；然而URTSS解决高维问题时效果不佳且响应速度慢[20]。因此本文将CKS与EM算法相结合，有望使算法具有性能稳定、精度较高、收敛速度快以及计算量相对较小等特点。

1问题描述

考虑以下非线性离散系统状态空间模型中，对未知参数θ的估计问题：

(1)

(2)

为了对未知参数进行估计,需要对状态向量进行扩维。状态扩维后的状态向量表示为xk]T,原系统方程可重新构造为

(3)

2基于EM的辨识算法

将EM算法应用于上述非线性离散系统状态空间模型中,将扩维状态向量XN作为EM中的不可观测数据,量测向量ZN作为EM中的可观测数据,λ为未知统计量。由EM算法的原理可知,其代价函数为

(4)

在给定参数λ的条件下,系统状态和观测值得对数联合概率分布为

(5)

(6)

式中

内蒙古锡盟～江苏泰州±800kV特高压直流输电工程天津段路径总长196.153km，线路经过天津市蓟县、宝坻区、武清区、西青区、静海区、河北省廊坊市6个地市级行政区。全线以平地为主，沿线树木较多。线路天津段跨越河道18条，蓄滞洪区5个。依据《中华人民共和国防洪法》规定：在洪泛区、蓄滞洪区内建设非防洪建设项目，应当就洪水对建设项目可能产生的影响和建设项目对防洪可能产生的影响作出评价，编制洪水影响评价报告，提出防御措施。

(7)

(8)

(9)

(10)

(11)

(12)

(13)

式中

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

为了更加简洁明了,本文提出的算法步骤归纳如下:

步骤5若达到收敛条件,则结束;否则返回步骤2。

3算法在飞机参数辨识中的应用

3.1ATTAS试验飞机横向线性模型

飞机的输入数据取自ATTAS试验飞机的飞行数据,其试验数据及简化线性模型来自参考文献[1]。ATTAS试验飞机的横向线性模型如下所示

(29)

(30)

式中,t为时间/s,p为滚转角速率/(rad/s),r为偏航角速率/(rad/s),β为侧滑角/(rad),δα为副翼偏转角/(rad),δr为方向舵偏转角/(rad)。在上述模型中,将β作为输入变量,并假定其能够无误差测量。ωk(·)为过程噪声。在飞行试验中滚转角速率和偏航角速率作为观测值,其观测方程为

(31)

(32)

3.2实验仿真设计

仿真输入量(副翼偏转角δσ/(rad)、方向舵偏转角δr/(rad)和侧滑角β/(rad))输入数据的采样时间为15s,共有750个采样点,采样间隔为0.02s,其随时间的变化如图1所示。未知向量Θ的真值如表1所示。仿真状态向量初值取为量测值的第一组值,即p0=9.662×10-3,r0=-9.856×10-4;初始状态向量协方差为零,即P0=014×14;过程噪声方差和量测噪声方差分别取为Q=10-8I2和R=10-6I2。在上述输入向量及参数设定的基础上,可得到输出(滚转角速率p和偏航角速率r)如图2所示。

将本文提出的算法(EM-CKS)与文献[20]中的算法(EM-URTSS)用于ATTAS试验飞机的横向线性气动模型进行辨识。算法在MATLAB环境下运行,其最大迭代次数设为i=5 000,并且将达到规定迭代次数作为唯一终止条件。

图1　输入量随时间的变化

3.3实验结果及分析

通过仿真实验可知,EM-CKS方法比EM-URTSS算法的计算时长略短,其单步迭代计算机时分别为6.691有s和7.057s。原因在于EM-CKS方法的Sigma点个数少于EM-URTSS,且在计算Sigma点之前,能够提前对采样点进行计算,其权值为固定值。EM-URTSS方法在计算Sigma点时运算相对于EM-CKS稍显复杂,均值和方差的权值则需分别计算,且每一组点的权值都不相同。

图2　输出量随时间的变化

EM-CKS和EM-URTSS参数估计值如表1所示。

表1　未知参数的真值与估计值

由表中可知EM-CKS与EM-URTSS的估计精度基本一致。2种方法的待估计参数(仅选择其中2个)随迭代次数的收敛曲线如图3所示。从图3中能够明显看出,在参数估计结果和估计精度基本一致的情况下,EM-CKS方法的收敛速度显著快于EM-URTSS方法,只需很少的迭代次数就可以达到收敛状态。在2种算法均使用EM进行噪声估计的情况下,这种现象说明本文算法所选用的CKS在保证估计精度的同时,还具有较好的收敛特性。

图3　参数随迭代次数的收敛曲线

EM-CKS和EM-URTSS方法对噪声方差有着基本类似的估计值,与仿真真值接近,其均为

噪声方差Q和R的迹随迭代次数的收敛曲线如图4所示。图像表明EM-CKS和EM-URTSS方法都能对未知参数进行很好的估计,但是EM-CKS的收敛速度明显高于EM-URTSS。

图4　噪声方差的迹随迭代次数的收敛曲线

4结论

本文提出了基于期望最大化(EM)算法与求容积卡尔曼平滑器(CKS)的联合估计参数辨识算法,对算法进行了详细的推导和设计,给出了算法的计算步骤。算例验证表明,本文提出的算法能够同时对噪声、状态量及未知参数进行估计,具有较好的辨识效果;算法能够有效地应用于飞行器气动参数辨识问题中,通过飞机的输入及响应对飞行器气动模型中的未知参数进行辨识。通过与其他算法在计算性能方面进行对比,发现本文提出的算法辨识精度高、收敛速度快,具有比较明显的优势。

参考文献：

[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)及航天技术支撑基金(2014-HT-XGD)资助

作者简介：刘艺蕾(1990—)，女，中航工业第一飞机设计研究院助理工程师，主要从事适航技术研究。

中图分类号：V217

文献标志码：A

文章编号:1000-2758(2016)04-0587-06