一种飞行器最大似然参数辨识方法与流程

本发明涉及飞行器参数辨识，具体涉及一种飞行器最大似然参数辨识方法。

背景技术：

1、本节中的陈述仅提供与本公开相关的背景信息，并且可能不构成现有技术。

2、利用飞行试验数据辨识稳定和控制导数，修正地面预测结果，是飞行器设计与评估的一个重要环节。似然函数因其卓越的理论性质，已成为飞行器参数辨识最广泛使用的估计准则。从本质上来说，飞行器参数辨识就是在给定模型结构的情况下，采用统计估计方法确定模型中的未知参数，以实现飞行动力学系统输出与飞行试验实测数据的最优匹配。

3、对于同时包含过程噪声和测量噪声的动力学系统，有两类不同的参数估计方法：一类通过优化某一代价函数来估计未知的系统参数和噪声统计量；另一类通过人为地将系统参数定义为附加状态变量，从而将参数估计问题转化为滤波问题。

4、滤波器误差法(fem)属于第一类参数估计方法，采用高斯-牛顿法优化关于系统参数和噪声统计量的似然函数。与最常用的仅考虑测量噪声的输出误差法(oem)相比，在由湍流产生过程噪声的情况下，应用fem方法通常能够获得更准确的气动参数估计结果。但fem方法要求待辨识参数的初始估计在最优估计附近，并且算法较复杂，可能存在收敛性和数值奇异性问题。

5、在第二类参数估计方法中，将未知的系统参数作为附加状态变量，从而将参数估计问题转化为增广状态估计问题。kalman滤波器为线性动力学模型提供了状态变量的最佳估计。然而，即使系统状态方程是线性的，增广后的状态方程也总是非线性的，需要采用扩展kalman滤波器(ekf)或无迹kalman滤波器(ukf)。这些传统的滤波器需要先验指定过程噪声和测量噪声协方差。滤波方法都是基于最优滤波，在确定时间点使用测量数据更新未知量估计，从初始时间递推到最终时间。相比之下，平滑方法基于整个测量历史计算未知量的最优估计。因此，它有可能提供比滤波方法更准确的估计。一类候选的平滑方法是使用非线性规划，它们包括非凸数值优化算法，如序列二次规划和内点法，其计算昂贵且不容易实现。相比之下，无迹rauch–tung–striebel平滑器(urtss)的计算成本低且相对容易实现。urtss使用最优后向平滑，递归计算对ukf前向估计的校正。通过使用无迹变换，urtss比基于一阶近似的扩展kalman平滑器更精确地逼近概率密度函数(pdf)。然而，滤波方法和平滑方法的一个主要局限性是它需要另一种方法来估计问题中的统计量，即初始状态(或初始增广状态)的均值和协方差以及噪声协方差。yokoyama在文献(yokoyama n.parameterestimation of aircraft dynamics via unscented smoother with expectation-maximization algorithm.j guidance,2011,34(2):426–436)中，针对非线性动力学系统，发展了将平滑器与期望最大化(em)算法相结合的参数估计方法，但没有给出参数估计准度的评价方法。

技术实现思路

1、本发明的目的在于：针对同时包含过程噪声和测量噪声的飞行器动力学系统，提供一种具有鲁棒收敛性的最大似然参数辨识方法，以解决滤波和平滑类方法中的未知噪声统计特性估计问题和参数估计准度评价问题，克服fem方法对待辨识参数初始估计的依赖性和可能存在的收敛性和数值奇异性问题。

2、本发明的技术方案如下：

3、一种飞行器最大似然参数辨识方法，包括：

4、步骤s100：将未知的系统参数β视为附加状态变量，从而将飞行器动力学系统表示为离散形式的增广非线性动力学系统，待辨识参数为系统参数β、初始增广状态均值和协方差p0、过程噪声协方差q、测量噪声协方差r；

5、步骤s200：对于给定的p0、q、r初始估计，采用平方根无迹卡尔曼滤波器进行前向状态估计；

6、步骤s300：采用无迹rauch–tung–striebel平滑器进行后向状态平滑；

7、步骤s400：更新p0、q、r的估计；

8、步骤s500：重复执行步骤s200～步骤s400，直至收敛或达到指定迭代步数，计算参数估计值

9、进一步地，所述步骤s100，包括：

10、步骤s101：将飞行器动力学系统的数学模型表示为下列离散形式的非线性动力学系统：

11、xk+1＝f(xk,uk,β)+wk

12、yk＝h(xk,uk,β)+vk

13、式中，

14、xk61为k+1时刻的nx维状态向量，k为采样时刻标识，k＝0,1,…,n，nx表示状态量x的维度；

15、f为非线性状态函数向量；

16、xk为k时刻的nx维状态向量；

17、uk为k时刻的nu维输入向量，nu表示输入量u的维度；

18、wk为k时刻的nx维过程噪声向量；

19、yk为k时刻的ny维观测向量，ny为观测量y的维度；

20、h为非线性观测函数向量；

21、vk为ny维测量噪声向量；

22、步骤s102：定义na维增广状态向量从而构建增广状态方程，表示为：

23、

24、其中，na表示增广状态向量xa的维度，na＝nx+nβ，nβ为系统参数β的维度；

25、上标t表示转置；

26、为k+1时刻的增广状态向量；

27、βk+1为k+1时刻的nβ维系统参数向量；

28、βk为k时刻的nβ维系统参数向量。

29、进一步地，所述步骤s200，包括：

30、步骤s201：初始化：计算平方根设置k＝0；设置表示初始增广状态的期望，表示初始增广状态协方差的平方根；

31、步骤s202：定义(2na+1)个sigma点及其关联权重如下：

32、

33、其中，表示k时刻的第j个增广状态sigma点；

34、表示的均值，表示在条件yk下增广状态向量的概率密度函数，

35、κ为二次缩放参数；

36、表示的协方差平方根；

37、表示的第j列；

38、调节参数α决定了sigma点在周围的分布；对于高斯分布，调节参数β＝2；

39、表示第j个sigma点均值的权值；

40、表示第j个sigma点协方差的权值；

41、步骤s203：进行时间更新；

42、步骤s204：进行测量更新；

43、步骤s205：将k增加1；如果k＝n，退出；否则，返回步骤s202。

44、进一步地，所述步骤s203，包括：

45、计算sigma点和

46、

47、其中，为在条件yk下k+1时刻第j个增广状态sigma点的预测值；为在条件yk下k+1时刻第j个观测量sigma点的预测值；

48、和分别表示向量中与x和β相对应的部分；

49、和分别表示向量中与x和β相对应的部分；

50、uk+1为k+1时刻的nu维输入向量；

51、计算增广状态和观测量的均值：

52、

53、其中，表示在条件yk下k+1时刻增广状态向量的预测值；

54、表示在条件yk下k+1时刻观测向量的预测值；

55、计算增广状态协方差平方根：

56、

57、其中，表示为计算而定义的中间变量；

58、qr(·)和cholupdate(·)为算子；

59、表示第1～2na个sigma点协方差的权值；

60、表示k+1时刻第1～2na个增广状态sigma点的预测值；

61、表示k+1时刻增广状态协方差平方根的预测值；

62、表示k+1时刻第0个增广状态sigma点的预测值；

63、表示第0个sigma点协方差的权值；

64、计算观测协方差平方根：

65、

66、其中，表示为计算而定义的中间变量；

67、表示k+1时刻第1～2na个观测量sigma点的预测值；

68、表示k+1时刻观测协方差平方根的预测值；

69、表示k+1时刻第0个观测量sigma点的预测值；

70、计算互相关协方差：

71、

72、其中，表示k+1时刻增广状态xa与观测量y互相关协方差的预测值。进一步地，所述步骤s204，包括：

73、计算最优增益：

74、

75、kk+1表示k+1时刻的增益矩阵；

76、计算增广状态均值：

77、

78、表示k+1时刻增广状态向量的校正值；

79、yk+1表示k+1时刻观测向量测量值；

80、计算增广状态协方差平方根：

81、

82、其中，u表示为计算而定义的中间变量；

83、表示k+1时刻增广状态协方差的校正值。

84、进一步地，所述步骤s300，包括：

85、步骤s301：将均值和协方差设为平方根无迹卡尔曼滤波器得到的值；设置k＝n-1；

86、步骤s302：计算ck+1、dk、

87、

88、其中，ck+1表示与的互相关协方差；

89、dk表示为计算和计算而定义的中间变量；

90、表示在条件yn下k+1时刻增广状态向量的期望(平滑结果)；

91、表示ukf中k时刻增广状态协方差校正值；

92、表示在条件yn下k+1时刻增广状态向量的协方差(平滑结果)；

93、表示ukf中k+1时刻增广状态协方差预测值；

94、步骤s303：计算

95、

96、其中，表示在条件yn下k时刻观测向量的期望(平滑结果)；

97、和分别为增广状态向量中与x和β相关的部分；

98、步骤s304：如果k＝0，则退出；否则，将k减1并返回步骤s302。

99、进一步地，所述步骤s400中的p0、q、r采用下式进行更新：

100、

101、其中，表示在条件yn下初始增广状态向量的期望；

102、表示在条件yn下初始增广状态向量的协方差；

103、表示中与x相关的部分，表示在条件yn下k时刻增广状态向量的协方差。

104、进一步地，所述步骤s500中的参数估计值采用下式计算：

105、

106、进一步地，还包括：

107、步骤s600：计算参数估计的标准差

108、进一步地，所述步骤s600，包括：

109、步骤s601：采用中心差分来计算观测灵敏度sk，即对于中的每一个元素上的小扰动+δβi和-δβi，分别进行执行步骤s200和步骤s300，计算扰动响应和相应的灵敏度为：

110、

111、其中，表示小扰动+δβi下的k时刻观测量平滑值；

112、表示小扰动-δβi下的k时刻观测量平滑值；

113、为偏导符号；

114、观测灵敏度为：

115、

116、sk表示k时刻的观测灵敏度；

117、步骤s602：计算信息矩阵m和输出残差自相关函数

118、

119、υk表示k时刻的残差；

120、步骤s603：计算参数估计的协方差矩阵：

121、

122、其中，表示参数估计的协方差矩阵；

123、si和sj分别表示i和j时刻的观测灵敏度；

124、步骤s604：计算参数估计的标准差：

125、

126、其中，表示参数估计的标准差；

127、diag(·)表示由矩阵对角线元素构成的向量。

128、与现有的技术相比本发明的有益效果是：

129、1.本发明给出的飞行器参数最大似然估计具有优良的统计特性，是渐进一致和渐进有效的无偏估计；

130、2.本发明采用无迹变换，适用于一般非线性动力学系统；

131、3.本发明采用srukf代替标准ukf，有效避免由于计算截断误差累积可能导致状态协方差矩阵失去正定性甚至对称性，进而导致滤波算法中断的问题，并在一定程度上减少计算量；

132、4.本发明具有很好的鲁棒收敛性，不依赖于系统参数和噪声统计特性的初始估计；

133、5.本发明在获得系统参数最大似然估计的同时，给出比常用的cramer-rao界更客观的参数估计标准差。

技术特征：

1.一种飞行器最大似然参数辨识方法，其特征在于，包括：

2.根据权利要求1所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s100，包括：

3.根据权利要求2所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s200，包括：

4.根据权利要求3所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s203，包括：

5.根据权利要求4所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s204，包括：

6.根据权利要求5所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s300，包括：

7.根据权利要求1所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s400中的p0、q、r采用下式进行更新：

8.根据权利要求1所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s500中的参数估计值采用下式计算：

9.根据权利要求1所述的一种飞行器最大似然参数辨识方法，其特征在于，还包括：

10.根据权利要求9所述的一种飞行器最大似然参数辨识方法，其特征在于，所述步骤s600，包括：

技术总结
本发明公开了一种飞行器最大似然参数辨识方法，涉及飞行器参数辨识技术领域，包括：首先将飞行器动力学系统表示为离散形式的增广非线性动力学系统，待辨识参数为系统参数β、初始增广状态均值和协方差P<subgt;0</subgt;、过程噪声协方差Q、测量噪声协方差R；再对于给定的P<subgt;0</subgt;、Q、R初始估计，采用平方根无迹卡尔曼滤波器进行前向状态估计；然后采用无迹Rauch–Tung–Striebel平滑器进行后向状态平滑；再更新P<subgt;0</subgt;、Q、R的估计；最后重复执行上述步骤，直至收敛或达到指定迭代步数，计算参数估计值本发明，具有很好的鲁棒收敛性。

技术研发人员：汪清,郑凤麒,丁娣,岳茜,赵暾,程艳青,刘滔,周宇
受保护的技术使用者：中国空气动力研究与发展中心计算空气动力研究所
技术研发日：
技术公布日：2024/9/23