一种自适应容积卡尔曼滤波飞行器气动参数辨识方法与流程

本发明属于参数辨识，特别是涉及一种自适应容积卡尔曼滤波飞行器气动参数辨识方法。

背景技术：

1、参数辨识方法包括离线和在线两种，离线方法如最小二乘法、递推最小二乘法及带遗忘因子的最小二乘法等。在线方法则通常采用卡尔曼滤波、扩展卡尔曼滤波和无迹卡尔曼滤波等。相较于离线方法，在线方法能够实时估计模型参数，因此在控制系统模型中广泛应用，可随时修订参数以提高系统精度，因此更为普遍。

2、在飞行器参数辨识领域，扩展卡尔曼滤波(ekf)是最常见的方法。例如，有研究者利用ekf对战术导弹的参数进行辨识。然而，对于复杂的非线性系统，由于难以求解雅可比矩阵和海森矩阵，ekf存在计算量大和估计精度低的问题。为了解决这些问题，有人提出了一种改进的扩展卡尔曼滤波方法，基于迭代滤波理论，通过在ekf估计点上重新进行泰勒级数展开并迭代，有效减少了线性化误差，提高了辨识精度。然而，这并未从根本上解决计算复杂性的问题。

3、有学者提出一种基于无迹卡尔曼滤波(ukf)的巡飞弹气动参数在线辨识方法，相较于ekf，ukf具有更高的精度，并避免了求解雅可比矩阵和海森矩阵的问题。但是，ukf在处理高维强非线性系统时存在性能不稳定的情况。而且，在实际应用中，这些方法往往会面临协方差矩阵负定的问题。

4、基于此，本发明提出一种自适应容积卡尔曼滤波飞行器气动参数辨识方法。

技术实现思路

1、针对以上技术问题，本发明提供一种自适应容积卡尔曼滤波飞行器气动参数辨识方法。

2、本发明解决其技术问题采用的技术方案是：

3、一种自适应容积卡尔曼滤波飞行器气动参数辨识方法，所述方法包括以下步骤：

4、s100：获取飞行器的飞行参数，建立战术导弹的状态方程与量测方程；

5、s200：将待辨识参数加入状态方程中，建立增广之后的状态方程；

6、s300：对滤波器进行初始化，确定初始时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值；

7、s400：时间更新：根据增广之后的状态方程、上一时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值和过程噪声协方差矩阵来预测当前时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值；

8、s500：量测更新：根据量测方程、当前时刻的增广状态变量估计值、增广状态变量误差协方差矩阵估计值和测量噪声协方差矩阵计算得到更新后的状态估计和更新后的协方差矩阵；

9、s600：递推估计噪声方差，返回s400进行迭代，直至达到预设的迭代次数，得到最终的待辨识的参数估计。

10、优选地，s100中状态方程具体为：

11、

12、其中，vx、vy、vz分别是弹体坐标系轴向、法向和侧向三个分量的线速度；ωx、ωy、ωz分别是三个方向的角速度；ix、iy、iz分别是滚转、偏航和俯仰通道的转动惯量，x,y,z分别为轴向力、法向力、横向力，γ,ψ为滚转，偏航，俯仰角，mx,my,mz为各个方向上的力矩，h为飞行高度。

13、优选地，s100中的量测方程具体为：

14、ωx＝ωx+v1

15、ωy＝ωy+v2

16、ωz＝ωz+v3

17、

18、ψ＝ψ+v8

19、γ＝γ+v9

20、h＝h+v10

21、α＝arctg(-vy/vx)+v11

22、β＝arcsin(vz/v)+v12

23、其中，x、y、z分别为轴向力、法向力、横向力；nx、ny、nz分别为轴向力、法向力、横向力的过载；g0为重力加速度；v为速度；v1～v12分别为各观测量噪声，为高斯噪声服从高斯分布；α,β分别为攻角和侧滑角。

24、优选地，s200中待辨识参数具体为

25、

26、其中，cx,cy,cz分别为轴向力、法向力和横向力的气动力，cx0、为轴向力系数；cy0、为法向力系数；cz0、为横向力系数；

27、轴向力、法向力、横向力具体为：

28、

29、其中，q代表动压，表示气动流动时候的能量密度，计算公式为sm是导弹的一个特征面积；

30、s200具体为：

31、

32、z(k)＝h(x(k),u(k))+v(k)

33、其中，为增广状态矢量，其中，增广状态矢量为z(k)为量测矢量，量测矢量为u(k)为控制输入，控制输入参数为δx，δy，δz，w(k)，v(k)为系统噪声。

34、优选地，s300具体为：

35、

36、其中，表示初始状态的估计值，p0表示误差协方差矩阵。

37、优选地，s400包括：

38、s410：使用svd对协方差矩阵pk-1进行分解，进一步获得容积点，具体为：

39、

40、其中，s为对角矩阵，s＝diag{s1,s2,…,sn}经过协方差矩阵的svd分解，u为左奇异矩阵，v为右奇异矩阵；

41、pk-1＝uk-1sk-1vtk-1

42、

43、其中，xi,k-1表示在时刻k-1的第i个容积点，ui,k-1表示协方差矩阵pk-1分解得到的左奇异向量矩阵中的第i列，表示协方差矩阵pk-1分解得到的奇异值矩阵中的第i个奇异值的平方根，表示在时刻k-1的状态预测的均值；i＝1,2,m表示容积采样点个数，使用3阶容积规则时容积采样点总数是状态维数的2倍，即m＝2n；[1]i表示点集[1]中的第i个点，其中符号[1]表示完整全对称点集；

44、[1]代表n维空间的点集，即：

45、

46、s420：计算通过非线性状态方程传播容积点：

47、x*i,k＝f(xi,k-1,uk-1)

48、其中，xi,k-1表示在k-1时刻的第i个容积点，f代表状态转移，x*i,k代表通过状态转移函数f变换后的第i个状态预测点；uk-1是k-1时刻的控制输入或外部输入；

49、s430：计算k-1时刻的状态预测值和误差协方差矩阵的值：

50、

51、其中，qk-1是过程噪声协方差矩阵。

52、优选地，s500包括：

53、s510：对pk|k-1进行svd分解，具体为：

54、pk|k-1＝uk|k-1sk|k-1vtk|k-1

55、

56、s520：通过非线性量测方程计算传播的容积点，具体为：

57、zi,k|k-1＝h(xi,k|k-1)

58、其中，h是量测方程，zi,k|k-1是通过量测方程将状态映射到观测空间中得到的观测的预测值；

59、s530：估计k时刻的观测预测值，并计算测量协方差矩阵和状态与测量的交叉协方差矩阵，具体为：

60、

61、其中，是通过容积点数值积分得到的预测测量均值，pzz,k|k-1是测量协方差矩阵，rk是测量噪声协方差矩阵，描述了测量噪声的统计特性，pxz,k|k-1是状态与测量的交叉协方差矩阵，描述了状态与测量之间的协方差；

62、s540：估计滤波增益和k时刻状态估计值与协方差矩阵，具体为：

63、wk＝pxz,k|k-1pxz,k|k-1-1

64、

65、其中，wk是卡尔曼增益矩阵，用于将预测状态与测量信息结合；是更新后的状态估计，zk是实际测量值，是预测的测量均值，pk是更新后的协方差矩阵。

66、优选地，s600中递推估计噪声方差，具体为：

67、

68、其中，为加强新观测信息的作用，引入渐消记忆指数，将a＝(1-bk)/(1-b)代入即可，b为遗忘因子，0<b<1；

69、在实际应用中，式存在减号，使得噪声方差阵可能失去正定性而导致滤波发散，使用噪声有偏估计器可阻止滤波发散：

70、

71、因此提出如下改进的噪声估计器：

72、

73、其中，z是任何非零向量，通常选择标准基向量或随机向量作为测试向量，用来判断矩阵和是否正定；

74、在估计噪声方差的过程中，当为非半正定或者为非正定时，使用第二个公式来对其进行修正；

75、如果不满足为非半正定或者为非正定，采用第一个公式来进行估计；

76、如果考虑到时间效率的重要性，并愿意在一定程度上牺牲滤波的精确度，采用第三个公式来进行噪声的有偏估计。

77、上述一种自适应容积卡尔曼滤波飞行器气动参数辨识方法，结合了基于奇异值分解的卡尔曼滤波算法和sage-husa自适应滤波，实现了对协方差矩阵的奇异值分解，并根据实时测量信息修正噪声方差，从而提高了参数估计的稳定性、精度和适应性。这一创新技术能够在飞行器参数辨识等领域中广泛应用，为实时参数估计提供了更加鲁棒和精确的解决方案，为系统性能的提升提供了重要保障。