噪声方差阵R和初始估计误差协方 差阵P;
[0081] 步骤2、状态量变化率解算：根据真空速投影到机体系的映射关系，计算状态量在 当前k时刻的变化率，具体步骤包括：
[0082] 步骤201、载入k时刻捷联惯导系统输出的机体系加速度，包括机头方向加速度 A、机身右翼加速度匕和机身垂向加速度虼，加速度单位均为ft/s2 (英尺每平方秒）；
[0083] 步骤202、基于真空速投影到机体系的映射关系，进行微分，计算状态量的变化 率：
[0087] 式中，冷、和f为k时刻状态量的变化率；ak、0 #V〇为k时刻状态量的估 计值，初始时刻状态量的估计值取初值；
[0088] 步骤3、状态量一步预测：根据状态量变化率，结合采样步长T，对状态量进行一步 预测：
[0092]式中，.、.爲+:1和[.1^.+1为k时刻状态量的一步预测；
[0093] 步骤4、一步预测均方误差解算：计算状态一步转移矩阵F并离散化，结合k时刻 的估计误差协方差阵Pk，解算一步预测均方误差Pk+lik，具体步骤包括：
[0094] 步骤401、根据k时刻状态量的估计值和机体系加速度，计算状态一步转移矩阵F， F中的元素为：
[0103]F3j3= 0 ；
[0104] 步骤402、对F离散化，得到〇k+lik=I+FT，I为单位矩阵；
[0105] 步骤403、结合k时刻的估计误差协方差阵Pk，解算一步预测均方误差Pk+1,k:
[0106]
[0107] 式中，PkSk时刻的估计误差协方差阵，初始时刻的估计误差协方差阵为P;
[0108] 步骤5、滤波增益解算：根据捷联惯导系统和飞行控制系统提供的信息，建立气动 参数与大气参数的函数关系，解算量测矩阵，具体步骤包括：
[0109] 步骤501、载入下一时刻k+1时刻的捷联惯导系统输出的机体系角速度、飞行控制 系统的飞行控制量和空气密度P，空气密度的单位为slug/ft3 (斯勒格每立方英尺）；
[0110] 步骤502、根据气动模型，结合捷联惯导系统和飞行控制系统提供的信息，将气动 系数改写为大气参数的函数，解算量测矩阵H，H中元素为：
[0129] 式中，山为函数中的常数系数，由气动系数的函数系数、载体的速度、角速度、姿态 以及转动惯量解算获得，i= 1~43 ;
[0130] 步骤503、解算滤波增益
[0131] 步骤6、估计均方误差解算：根据滤波增益Kk+1、量测矩阵H、一步预测均方误差 Pk+lk和测量噪声方差阵R，解算估计均方误差Pk+1:
[0132]
[0133] 步骤7、状态估计：结合状态量的一步预测，根据载体的气动模型，解算k+1时刻的 观测量估计值，从而实现状态量的估计，即大气参数的估计，具体步骤包括：
[0134] 步骤701、载入k+1时刻捷联惯导系统输出的机体系角加速度，包括机头方向角加 速度為+1、机身右翼角加速度金+1和机身垂向角加速度鳥，:角加速度单位均为rad/s2;载 入k+1时刻的飞行控制量、机体系加速度'以+1、Fi+1、化+1;
[0135] 步骤702、将状态量在k时刻的一步预测、k+1时刻的飞行控制量、k+1时刻的机体 系加速度和k+1时刻的机体系角加速度代入到气动模型中，解算观测量估计值：
[0142] 式中，d44为函数中的常数，由气动系数的函数系数、载体的速度、角速度、姿态以及 转动惯量解算获得；
[0143]步骤703、记k+1时刻的观测量为
k+1时刻的观测量估计值为
，状态量在 k时刻的一步预测戈
，状态量在k+1时刻的一步预测为
则有状态估计值
在获得状态估 计值后，返回步骤2，继续进行后续的大气参数解算。
[0144] 为了评价本发明提出的真空速解算方法的性能，设计了仿真程序，在matlab平台 上对算法进行验证，结构如图2所示，该仿真程序包括以下步骤：
[0145] (1)设置航迹，根据载体气动模型，进行飞行器配平；
[0146] (2)进行飞行仿真，生成与航迹匹配的飞行控制量、惯导数据、大气参数；
[0147] (3)对飞行控制量、气动模型、惯导数据注入误差；
[0148] (4)使用注入误差后的飞行控制量、气动模型、惯导数据进行大气参数解算，与飞 行仿真得到的真实值进行对比分析；
[0149] 实验结果如图3所示，表明本发明的大气参数解算结果与真实值基本重合，攻角 误差均值为〇. 16268°，侧滑角误差均值为0.05505°，真空速误差均值为-4. 82748m/s，证 明了本发明解算大气参数方法的正确性和有效性。
[0150] 以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人 员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应 视为本发明的保护范围。
【主权项】
1. 一种基于捷联惯导与飞行控制系统的大气参数解算方法，其特征在于：包括w下步 骤： 步骤1、初始信息设置：状态量为攻角a、侧滑角P和真空速Vt，设置攻角a的初值为 载体俯仰角0，侧滑角P的初值为零，攻角、侧滑角单位均为rad,真空速Vt的初值为载体 速度Vb，单位为ft/s;设置系统噪声方差阵Q、量测噪声方差阵R和初始估计误差协方差阵 P; 步骤2、状态量变化率解算：根据真空速投影到机体系的映射关系，计算状态量在当前k时刻的变化率，具体步骤包括： 步骤201、载入k时刻捷联惯导系统输出的机体系加速度，包括机头方向加速度捉、机 身右翼加速度&和机身垂向加速度听，加速度单位均为ft/s2; 步骤202、基于真空速投影到机体系的映射关系，进行微分，计算状态量的变化率：式中，&、#和跨为k时刻状态量的变化率；ak、Pk和Vt,k为k时刻状态量的估计值， 初始时刻状态量的估计值取初值； 步骤3、状态量一步预测：根据状态量变化率，结合采样步长T，对状态量进行一步预 测：式中，馬+1、么1和巧W为k时刻状态量的一步预测； 步骤4、一步预测均方误差解算：计算状态一步转移矩阵F并离散化，结合k时刻的估 计误差协方差阵Pk，解算一步预测均方误差Pw,k，具体步骤包括： 步骤401、根据k时刻状态量的估计值和机体系加速度，计算状态一步转移矩阵F，F中 的元素为：?3.3= 0 ； 步骤402、对F离散化，得到〇w,k=I+FT，I为单位矩阵； 步骤403、结合k时刻的估计误差协方差阵Pk，解算一步预测均方误差Pw,k:式中，Pk为k时刻的估计误差协方差阵，初始时刻的估计误差协方差阵为P; 步骤5、滤波增益解算：根据捷联惯导系统和飞行控制系统提供的信息，建立气动参数 与大气参数的函数关系，解算量测矩阵，具体步骤包括： 步骤501、载入下一时刻k+1时刻的捷联惯导系统输出的机体系角速度、飞行控制系统 的飞行控制量和空气密度P，空气密度的单位为slug/ft3; 步骤502、根据气动模型，结合捷联惯导系统和飞行控制系统提供的信息，将气动系数 改写为大气参数的函数，解算量测矩阵H，H中元素为：式中，di为函数中的常数系数，由气动系数的函数系数、载体的速度、角速度、姿态W及 转动惯量解算获得，i= 1~43 ; 步骤503、解算滤波增益步骤6、估计均方误差解算：根据滤波增益Kw、量测矩阵H、一步预测均方误差和 测量噪声方差阵R，解算估计均方误差Pw:步骤7、状态估计：结合状态量的一步预测，根据载体的气动模型，解算k+1时刻的观测 量估计值，从而实现状态量的估计，即大气参数的估计，具体步骤包括： 步骤701、载入k+1时刻捷联惯导系统输出的机体系角加速度，包括机头方向角加速度 每,1、机身右翼角加速度这+,和机身垂向角加速度為+1，角加速度单位均为rad/s2;载入k+1 时刻的飞行控制量、机体系加速度化+1、n+,、馬 步骤702、将状态量在k时刻的一步预测、k+1时刻的飞行控制量、k+1时刻的机体系加 速度和k+1时刻的机体系角加速度代入到气动模型中，解算观测量估计值：式中，d44为函数中的常数，由气动系数的函数系数、载体的速度、角速度、姿态W及转动 惯量解算获得； 步骤703、记k+1时刻的观测量为k+1时刻的观测量估计值为状态量在 k时刻的一步预测为氧w= 反+1巧.w了，状态量在k+1时刻的一步预测为 氣+1 =[解+1A+1K.W了，则有状态估计值义。1 =式+1 +K;+,(Z;_+i-马.+1);在获得状态估 计值后，返回步骤2，继续进行后续的大气参数解算。
【专利摘要】本发明公开了一种基于捷联惯导与飞行控制系统的大气参数解算方法，基于攻角、侧滑角的定义，建立真空速投影到机体系的映射关系，结合高空大气流动速度稳定的特点，利用惯导信息得到大气参数的一步预测模型，根据飞行器气动模型，建立气动参数与大气参数的函数关系，构建起惯导量测值与大气参数的数学关系，通过扩展卡尔曼滤波的方法对大气参数进行求解。该方法为大气参数的测量提供了冗余手段，克服了高超声速、大机动情况下大气参数的测量难点，能够在飞行包线内实时提供大气参数。
【IPC分类】G06F17/16, G06F17/15, G01W1/02
【公开号】CN105005099
【申请号】CN201510388978
【发明人】李荣冰, 陆辰, 刘建业, 雷廷万, 郭毅, 廖自威
【申请人】南京航空航天大学
【公开日】2015年10月28日
【申请日】2015年7月3日