技术领域
[0001] 本
发明涉及航空航天捷联惯性导航技术领域，具体是涉及一种适用于机翼挠曲变形的机载自适应传递对准算法。
背景技术
[0002] 捷联惯性
导航系统具有反应时间短、可靠性高、体积小、重量轻等优点，广泛应用于飞机、舰船、导弹等军用和民用导航领域，具有重要的国防意义和巨大的经济效益。
[0003] 机载武器上的惯导系统传递对准是指在载机运动条件下，机载武器上的惯导系统利用飞机上处于导航状态的高
精度主惯导系统信息进行初始对准的方法。能否快速而准确地对机载武器上的惯导系统进行初始对准，在很大程度上决定着飞机等作战平台上武器系统的快速反应和精确打击能
力。因此，动
基座传递对准技术是飞机、舰船、潜艇等运动平台发射武器的一项关键技术。
[0004] 传统传递对准算法一般采用“速度加
姿态匹配”，该匹配方法虽然有很多优点，但是受机翼形变影响很大，尤其是飞机地面
加速爬升段以及飞机空中大机动时影响非常明显，从而影响对准精度。
发明内容
[0005] 本发明的目的是为了克服上述背景技术的不足，提供一种适用于机翼挠曲变形的机载自适应传递对准算法，在机翼大形变时能够实现精确对准，降低对准时对飞机机动的要求。
[0006] 本发明提供一种适用于机翼挠曲变形的机载自适应传递对准算法，包括以下步骤：
[0007] 根据机载主惯导信息计算子惯导初始
位置、速度、姿态信息，并进行惯导捷联解算；
[0008] 根据飞行高度，对观测噪声Rk进行调整，并进行子惯导的对准解算。
[0009] 在上述方案的
基础上，当滑跑
起飞机翼形变大时，使
滤波器工作在速度匹配模式下，将观测噪声Rk中姿态对应项调大；飞机升空后机翼形变稳定时，使滤波器工作在速度加姿态匹配模式下，将观测噪声Rk中姿态对应项调小。
[0010] 在上述方案的基础上，所述根据机载主惯导信息计算子惯导初始位置、速度、姿态信息，具体包括以下步骤：
[0011] 根据机载主惯导下发的数据，对弹载子惯导的初值进行初始化；
[0012] 根据初始化时刻机载的姿态数据计算四元数初值，定义姿态四元数如下：
[0013] Q(q0,q1,q2,q3)＝q0+q1i+q2j+q3k
[0014] 计算公式如下：
[0015]
[0016]
[0017]
[0018]
[0019] 其中：q0(0),q1(0),q2(0),q3(0)为弹载子惯导姿态四元数初值，无量纲； γ,ψ为
俯仰角、
滚转角和
偏航角初值，均取自主惯导对应时刻值，单位：弧度。
[0020] 在上述方案的基础上，根据机载主惯导下发的数据，对弹载子惯导的初值进行初始化，具体包括以下步骤：
[0021] 位置和速度均以机载主惯导的值进行赋值：
[0022] Ls0＝Lm
[0023] λs0＝λm
[0024] hs0＝hm
[0025]
[0026]
[0027]
[0028] 其中：Ls0,λs0,hs0, 为弹载子惯导地理纬度、经度、高度、东向速度、北向速度与天向速度初值，单位分别为：弧度、弧度、米、米/秒、米/秒、米/秒；Lm,λm,hm,为机载主惯导下发的地理纬度、经度、高度、东向速度、北向速度与天向速度，单位分别为：弧度、弧度、米、米/秒、米/秒、米/秒。
[0029] 在上述方案的基础上，所述进行惯导捷联解算，具体包括以下步骤：
[0030] 进行姿态更新计算：
[0031]
[0032]
[0033] 其中：Δθx,Δθy,Δθz为本体系三轴角增量，单位：弧度；
[0034] q0(tk+1),q1(tk+1),q2(tk+1),q3(tk+1)为递推更新后的姿态四元数；
[0035] q0(tk),q1(tk),q2(tk),q3(tk)为递推更新前的姿态四元数；
[0036] 进行速度更新计算：
[0037]
[0038]
[0039]
[0040]
[0041] 其中：ΔVx,ΔVy,ΔVz为弹体
坐标系三个方向的视速度增量，单位：米/秒；
[0042] 为递推下一步的北向、天向与东向速度，单位均为：米/秒；
[0043] 为北向、天向与东向的牵连加速度，单位：米/秒2；
[0044] 为北向、天向与东向的哥氏加速度，单位：米/秒2；
[0046] 进行位置更新计算：
[0047]
[0048]
[0049]
[0050] 其中：Lk,λk,hk为当前步的地理纬度、经度与高度值，单位分别为：弧度、弧度与米；
[0051] Lk+1,λk+1,hk+1为递推下一步的地理纬度、经度与高度值，单位分别为：弧度、弧度与米；
[0052] 为当前步的北向、东向与天向速度，单位均为：米/秒。
[0053] 在上述方案的基础上，所述进行子惯导的对准解算，具体包括以下步骤：
[0054] 步骤201，使用卡尔曼滤波进行动基座传递对准，将子惯导姿态失准角估计出来；选取子惯导三轴方向的姿态失准角、东向与北向速度误差、3轴陀螺漂移与安装偏差角作为状态向量：
[0055] X＝[φx,φy,φz,δVE,δVN,εx,εy,εz,ψx,ψy,ψz]T
[0056] 其中：X为系统状态；
[0057] φx,φy,φz为子惯导姿态失准角，单位：弧度；
[0058] δVE,δVN为子惯导东向、北向速度误差，单位：米/秒；
[0059] εx,εy,εz为子惯导三轴陀螺漂移，单位：弧度/秒；
[0060] ψx,ψy,ψz为子惯导三轴安装偏差角，单位：弧度；
[0061] 步骤202，根据子惯导脉冲数和子惯导工具误差补偿系数矩阵，计算一个
采样周期th内子惯导的本体系三轴速度增量，分别记为：ΔVx,ΔVy,ΔVz，则导航系视速度增量计算公式为：
[0062]
[0063] 其中： 为导航系三轴速度增量；
[0064] 步骤203，进行状态转移矩阵更新：状态转移矩阵实时计算的相关项如下：
[0065]
[0066]
[0067]
[0068]
[0069]
[0070]
[0071] 其中： 为附加
角速度，采用主惯导信息计算的值；
[0072] Ak(1:3,1:3)为矩阵Ak的1至3行与1至3列构成的3×3矩阵；
[0073] Ak(1:3,6:8)为矩阵Ak的1至3行与6至8列构成的3×3矩阵；
[0074] 步骤204，进行滤波时间更新计算：
[0075] Xk/k-1＝Ak/k-1Xk-1
[0076]
[0077] 其中：Xk/k-1为下一步预测状态；
[0078] Pk/k-1为下一步预测均方差；
[0079] Qk为过程噪声方差矩阵；
[0080] 步骤205，进行观测向量更新：其观测值计算如下：
[0081]
[0082]
[0083]
[0084]
[0085]
[0086] 其中： 为弹载子惯导东向、北向速度，单位：米/秒；
[0087] 为机载主惯导东向、北向速度，单位：米/秒；
[0089] 观测关系矩阵Hk初始化如下：
[0090]
[0091] 其需要实时计算的相关项如下：
[0092]
[0093]
[0094]
[0095] 步骤207，进行滤波测量更新计算：
[0096]
[0097] Pk＝(I-KkHk)Pk/k-1
[0098] Xk＝Xk/k-1+Kk(Zk-HkXk/k-1)
[0099] 其中：Kk为滤波增益；
[0100] Pk为估计均方差矩阵；
[0101] Xk为滤波估计状态向量；
[0102] Rk为量测噪声，通过判断
巡航高度，调节Rk中姿态对应项的大小，使对准在速度匹配及速度加姿态匹配之间切换。
[0103] 在上述方案的基础上，根据飞行高度，对观测噪声Rk进行调整，并进行子惯导的对准解算，还包括以下步骤：量测更新完成解算后，将滤波估计状态向量Xk+1和估计均方差矩阵Pk+1作为初值赋给下一周期的更新计算。
[0104] 在上述方案的基础上，使用卡尔曼滤波进行动基座传递对准，将子惯导姿态失准角估计出来，
加速度计漂移、加速度计和
陀螺仪零偏均为随机常值。
[0106] 本发明根据机载主惯导信息计算子惯导初始位置、速度、姿态信息，并进行惯导捷联解算；根据飞行高度，对观测噪声Rk进行调整，并进行子惯导的对准解算，实现滤波算法的自适应，有效地提高了速度加姿态匹配算法的环境适应能力及对准精度；算法不仅实现简单，运行可靠，而且快速性好、自适应性强，具有较佳的工程应用价值。
附图说明
[0107] 图1是本发明
实施例的适用于机翼挠曲变形的机载自适应传递对准算法的流程示意图；
[0108] 图2是本发明实施例的适用于机翼挠曲变形的机载自适应传递对准算法的原理示意图。
具体实施方式
[0109] 下面结合附图及具体实施例对本发明作进一步的详细描述。
[0110] 参见图1所示，本发明实施例提供一种适用于机翼挠曲变形的机载自适应传递对准算法，包括以下步骤：
[0111] 根据机载主惯导信息计算子惯导初始位置、速度、姿态信息，并进行惯导捷联解算；
[0112] 根据飞行高度，对观测噪声Rk进行调整，并进行子惯导的对准解算。
[0113] 作为优选的实施方式，根据飞行高度，对观测噪声Rk进行调整，具体包括以下步骤：当滑跑起飞机翼形变大时，使滤波器工作在速度匹配模式下，将观测噪声Rk中姿态对应项调大；飞机升空后机翼形变稳定时，使滤波器工作在速度加姿态匹配模式下，将观测噪声Rk中姿态对应项调小。
[0114] 本发明实施例根据机载主惯导信息计算子惯导初始位置、速度、姿态信息，并进行惯导捷联解算；根据飞行高度，对观测噪声进行调整，并进行子惯导的对准解算，实现滤波算法的自适应，有效地提高了速度加姿态匹配算法的环境适应能力及对准精度；算法不仅实现简单，运行可靠，而且快速性好、自适应性强，具有较佳的工程应用价值。
[0115] 作为优选的实施方式，所述根据机载主惯导信息计算子惯导初始位置、速度、姿态信息，具体包括以下步骤：
[0116] 根据机载主惯导下发的数据，对弹载子惯导的初值进行初始化；位置和速度均以机载主惯导的值进行赋值：
[0117] Ls0＝Lm
[0118] λs0＝λm
[0119] hs0＝hm
[0120]
[0121]
[0122]
[0123] 其中：Ls0,λs0,hs0, 为弹载子惯导地理纬度、经度、高度、东向速度、北向速度与天向速度初值，单位分别为：弧度、弧度、米、米/秒、米/秒、米/秒；Lm,λm,hm,为机载主惯导下发的地理纬度、经度、高度、东向速度、北向速度与天向速度，单位分别为：弧度、弧度、米、米/秒、米/秒、米/秒。
[0124] 根据初始化时刻机载的姿态数据计算四元数初值，定义姿态四元数如下：
[0125] Q(q0,q1,q2,q3)＝q0+q1i+q2j+q3k
[0126] 计算公式如下：
[0127]
[0128]
[0129]
[0130]
[0131] 其中：q0(0),q1(0),q2(0),q3(0)为弹载子惯导姿态四元数初值，无量纲； γ,ψ为俯仰角、滚转角和偏航角初值，均取自主惯导对应时刻值，单位：弧度。
[0132] 作为优选的实施方式，所述进行惯导捷联解算，具体包括以下步骤：
[0133] 进行姿态更新计算：
[0134]
[0135]
[0136] 其中：Δθx,Δθy,Δθz为本体系三轴角增量，单位：弧度；
[0137] q0(tk+1),q1(tk+1),q2(tk+1),q3(tk+1)为递推更新后的姿态四元数；
[0138] q0(tk),q1(tk),q2(tk),q3(tk)为递推更新前的姿态四元数；
[0139] 进行速度更新计算：
[0140]
[0141]
[0142]
[0143]
[0144] 其中：ΔVx,ΔVy,ΔVz为弹体坐标系三个方向的视速度增量，单位：米/秒；
[0145] 为递推下一步的北向、天向与东向速度，单位均为：米/秒；
[0146] 为北向、天向与东向的牵连加速度，单位：米/秒2；
[0147] 为北向、天向与东向的哥氏加速度，单位：米/秒2；
[0148] gk为当前重力加速度，米/秒2；
[0149] 进行位置更新计算：
[0150]
[0151]
[0152]
[0153] 其中：Lk,λk,hk为当前步的地理纬度、经度与高度值，单位分别为：弧度、弧度与米；
[0154] Lk+1,λk+1,hk+1为递推下一步的地理纬度、经度与高度值，单位分别为：弧度、弧度与米；
[0155] 为当前步的北向、东向与天向速度，单位均为：米/秒。
[0156] 作为优选的实施方式，参见图2所示，所述进行子惯导的对准解算，具体包括以下步骤：
[0157] 步骤201，使用卡尔曼滤波进行动基座传递对准，将子惯导姿态失准角估计出来，加速度计漂移、加速度计和陀螺仪零偏均为随机常值。
[0158] 选取子惯导三轴方向的姿态失准角、东向与北向速度误差、3轴陀螺漂移与安装偏差角作为状态向量：
[0159] X＝[φx,φy,φz,δVE,δVN,εx,εy,εz,ψx,ψy,ψz]T
[0160] 其中：X为系统状态；
[0161] φx,φy,φz为子惯导姿态失准角，单位：弧度；
[0162] δVE,δVN为子惯导东向、北向速度误差，单位：米/秒；
[0163] εx,εy,εz为子惯导三轴陀螺漂移，单位：弧度/秒；
[0164] ψx,ψy,ψz为子惯导三轴安装偏差角，单位：弧度；
[0165] 步骤202，根据子惯导脉冲数和子惯导工具误差补偿系数矩阵，计算一个采样周期th内子惯导的本体系三轴速度增量，分别记为：ΔVx,ΔVy,ΔVz，则导航系视速度增量计算公式为：
[0166]
[0167] 其中： 为导航系三轴速度增量；
[0168] 步骤203，进行状态转移矩阵更新：状态转移矩阵实时计算的相关项如下：
[0169]
[0170]
[0171]
[0172]
[0173]
[0174]
[0175] 其中： 为附加角速度，采用主惯导信息计算的值；
[0176] Ak(1:3,1:3)为矩阵Ak的1至3行与1至3列构成的3×3矩阵；
[0177] Ak(1:3,6:8)为矩阵Ak的1至3行与6至8列构成的3×3矩阵；
[0178] 步骤204，进行滤波时间更新计算：
[0179] Xk/k-1＝Ak/k-1Xk-1
[0180]
[0181] 其中：Xk/k-1为下一步预测状态；
[0182] Pk/k-1为下一步预测均方差；
[0183] Qk为过程噪声方差矩阵；
[0184] 步骤205，进行观测向量更新：其观测值计算如下：
[0185]
[0186]
[0187]
[0188]
[0189]
[0190] 其中： 为弹载子惯导东向、北向速度，单位：米/秒；
[0191] 为机载主惯导东向、北向速度，单位：米/秒；
[0192] 步骤206，进行观测矩阵更新：
[0193] 观测关系矩阵Hk初始化如下：
[0194]
[0195] 其需要实时计算的相关项如下：
[0196]
[0197]
[0198]
[0199] 步骤207，进行滤波测量更新计算：
[0200]
[0201] Pk＝(I-KkHk)Pk/k-1
[0202] Xk＝Xk/k-1+Kk(Zk-HkXk/k-1)
[0203] 其中：Kk为滤波增益；
[0204] Pk为估计均方差矩阵；
[0205] Xk为滤波估计状态向量；
[0206] Rk为量测噪声，通过判断巡航高度，调节Rk中姿态对应项的大小，使对准模式在速度匹配及速度加姿态匹配之间进行切换。
[0207] 步骤208，根据飞行高度，对观测噪声Rk进行调整，并进行子惯导的对准解算，还包括以下步骤：量测更新完成解算后，将滤波估计状态向量Xk+1和估计均方差矩阵Pk+1作为初值赋给下一周期的更新计算。
[0208] 卡尔曼滤波器是线性最小方差估计器，在设计系统时为了获得“最优”的滤波性能，需要根据工程经验选择合适的滤波参数。
[0209] 本领域的技术人员可以对本发明实施例进行各种
修改和变型，倘若这些修改和变型在本发明
权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围之内。
说明书中未详细描述的内容为本领域技术人员公知的现有技术。
