CN112595319B - 一种模型自适应补偿的返回弹道估计方法 - Google Patents

Abstract

本发明公开了一种模型自适应补偿的返回弹道估计算法，首先使用UKF滤波算法框架进行计算，在J2000.0地心惯性系下建立十维的系统状态模型；然后对滤波中的状态噪声协方差矩阵进行模型自适应补偿处理，实现滤波对弹道机动过程的适应。本发明解决了现有技术中存在的弹道估计算法抗野值能力较差、弹道曲线震荡幅度较大、光滑性稍差的问题。

Description

一种模型自适应补偿的返回弹道估计方法

技术领域

本发明属于航天导航技术领域，具体涉及一种模型自适应补偿的返回弹道估计方法。

背景技术

航天器返回弹道实时估计计算是航天器返回过程中一项重要的计算内容，但由于一些航天器返回过程制导方法复杂，或出于其它考虑，工业部门常常不愿意提供返回过程的制导率数据，导致返回弹道滤波计算时常常不能建立准确的预测模型，这就给外测飞行弹道的计算带来了一定难度。

且外测观测数据一般质量较难保证(相对于GPS、惯性导航等，外测数据质量较差)，当滤波算法对观测值的震荡幅度敏感时，弹道估计结果势必受到影响。

针对这种不能准确建立运动模型、且使用外测观测数据的返回段弹道估计问题，常采用仅需建立数值运算模型的当前统计模型、多项式模型等方法来建立扩展卡尔曼滤波计算过程，但这两类方法抗野值能力较差，在外测跟踪数据质量较差时，计算出的弹道曲线震荡幅度较大、光滑性稍差。

发明内容

本发明的目的是提供一种模型自适应补偿的返回弹道估计方法，解决了现有技术中存在的弹道估计算法抗野值能力较差、弹道曲线震荡幅度较大、光滑性稍差的问题。

本发明所采用的技术方案是，一种模型自适应补偿的返回弹道估计方法，具体按照以下步骤实施：

步骤1、使用UKF滤波算法框架进行计算，在J2000.0地心惯性系下建立十维的系统状态模型；

步骤2、对滤波中的状态噪声协方差矩阵进行模型自适应补偿处理，实现滤波对弹道机动过程的适应。

本发明的特点还在于，

步骤1具体如下：

步骤1.1、在J2000.0地心惯性系下建立系统状态向量为：

其中，x,y,z分别表示为航天器位置矢量的x分量、y分量、z分量；分别为速度矢量的x分量、y分量、z分量；a_x,a_y,a_z分别为未建模力即机动推力产生的加速度的x分量、y分量、z分量；β为单位加速度的加速度变化率；

步骤1.2、建立十维的系统状态外推模型为：

其中，g_m为返回段已建模力产生的加速度矢量，加速度矢量为x、y、z各方向的三维分量分别标记为g_mx、g_my、g_mz，已建模力包括地球引力、太阳引力、月亮引力、大气阻力，；

步骤1.3、对所述步骤1.2中的公式(1)进行数值积分，得到滤波计算的状态估计模型；

步骤1.4、观测模型建立在测站东北天坐标系下：

设某个测站在J2000.0地心惯性系下的位置、速度矢量分别记为由测站坐标可算出；航天器位置、速度矢量分别为r＝(x,y,z)、/>地心惯性系到测站东北天坐标系的转换矩阵为M_I；

则航天器在测站东北天坐标系下的位置矢量为M_I(r-R)，将该矢量的三个分量分别标记为ρ_x、ρ_y、ρ_y，

则可计算出测距值为：

方位角为：

A＝arctan(ρ_x/ρ_y) (3)

仰角为：

测距变化率为:

上述公式(2)～(5)即为建立的滤波计算的观测模型。

步骤2具体如下：

对UKF滤波过程，预测均值为将预测均值/>的速度项x、y、z三个方向上的分量标记为x_vi其中，i＝x,y,z，位置项的三分量标记为x_pi、状态更新值为/>将速度项的x、y、z三个方向上的分量标记为x_vi'，位置项三分量标记为x_pi'；

则可计算出速度模型误差量的近似值为：

σ_vi＝((x'_pi-x_pi)/Δt)

加速度模型误差近似值为：

σ_ai＝((x'_vi-x_vi)/Δt)

进而，Q矩阵自适应更新为：

λ为弹道机动频率系数；

至此，模型自适应补偿完成。

步骤2中弹道机动频率系数λ为0.001。

步骤2中取值均为0。

本发明的有益效果是，一种模型自适应补偿的返回弹道估计方法，建立了包含位置矢量、速度矢量的各分量在内的十维状态模型，能有效实现对返回段飞行弹道各分量的估计，且该模型不需要对返回段受力情况进行精确建模；算法通过对加速度及速度项的模型补偿，实现了对返回段机动过程的自适应辨识；算法采用UKF框架，不需要进行矩阵求导计算，算法框架明晰，易于实现；算法状态模型、观测模型、基于加速度模型补偿的自适应处理方法同样适用于容积卡尔曼滤波(CKF)、中心差分卡尔曼滤波(CDKF)等滤波框架下的使用。

附图说明

图1是本发明的计算流程图；

图2是本发明算法的实施例效果图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明一种模型自适应补偿的返回弹道估计方法，流程图如图1所示，具体按照以下步骤实施：

步骤1、使用UKF滤波算法框架进行计算，在J2000.0地心惯性系下建立十维的系统状态模型；

步骤1具体如下：

步骤1.1、在J2000.0地心惯性系下建立系统状态向量为：

其中，x,y,z分别表示为航天器位置矢量的x分量、y分量、z分量；分别为速度矢量的x分量、y分量、z分量；a_x,a_y,a_z分别为未建模力即机动推力产生的加速度的x分量、y分量、z分量；β为单位加速度的加速度变化率；

步骤1.2、建立十维的系统状态外推模型为：

其中，g_m为返回段已建模力产生的加速度矢量，加速度矢量为x、y、z各方向的三维分量分别标记为g_mx、g_my、g_mz，已建模力包括地球引力、太阳引力、月亮引力、大气阻力，；

步骤1.3、对所述步骤1.2中的公式(1)进行数值积分，得到滤波计算的状态估计模型；

由于返回段的外测观测数据建立在测站东北天坐标系下，且测距(标记为ρ)、方位角(标记为A)、仰角(标记为E)、测距变化率(标记为)的测量精度不同，一般来说，测距数据精度较高，而测角数据精度较低，若无差别地使用这几类数据，将可能使弹道确定精度降低。

步骤1.4、观测模型建立在测站东北天坐标系下：

设某个测站在J2000.0地心惯性系下的位置、速度矢量分别记为R,由测站坐标可算出；航天器位置、速度矢量分别为r＝(x,y,z)、/>地心惯性系到测站东北天坐标系的转换矩阵为M_I；

则航天器在测站东北天坐标系下的位置矢量为M_I(r-R)，将该矢量的三个分量分别标记为ρ_x、ρ_y、ρ_y，

则可计算出测距值为：

方位角为：

A＝arctan(ρ_x/ρ_y) (3)

仰角为

测距变化率为：

上述公式(2)～(5)即为建立的滤波计算的观测模型。

机动发生时UKF滤波的状态模型将与实际情况不匹配，则在滤波计算中，若将Q矩阵设置为固定值，则可能与航天器实际飞行过程速度的变化出现较大的偏差，导致滤波不能适应航天器的机动变化，其可能需要较长的时间才能够收敛，甚至有时会出现滤波发散情况。为此本发明对滤波过程的状态噪声协方差矩阵Q进行实时补偿，来自适应未知机动的发生。

步骤2、对滤波中的状态噪声协方差矩阵进行模型自适应补偿处理，实现滤波对弹道机动过程的适应。

步骤2具体如下：

对UKF滤波过程，预测均值为将预测均值/>的速度项x、y、z三个方向上的分量标记为x_vi其中，i＝x,y,z，位置项的三分量标记为x_pi、状态更新值为/>将速度项的x、y、z三个方向上的分量标记为x_vi'，位置项三分量标记为x_pi'；

则可计算出速度模型误差量的近似值为：

σ_vi＝((x'_pi-x_pi)/Δt)

加速度模型误差近似值为：

σ_ai＝((x'_vi-x_vi)/Δt)

进而，Q矩阵自适应更新为：

λ为弹道机动频率系数；

至此，模型自适应补偿完成。

步骤2中弹道机动频率系数λ为0.001。

步骤2中取值均为0。

实施例

1、建立系统状态方程：

设探测器t时刻的位置矢量为r＝(x,y,z)，速度矢量为在J2000.0地心惯性系下建立十维的系统状态向量为：

滤波的系统状态外推模型为：

其中，g_m为返回段已建模力(主要为地球引力、太阳引力、月亮引力、大气阻力等)产生的加速度矢量(其x、y、z各方向的三维分量分别标记为g_mx、g_my、g_mz)。

2、滤波起步计算：

采用多点最小二乘拟合计算滤波起步值：

具体地，取多个观测点的测量值(设某个点的测量值时间、测距、方位角、仰角分别为t_i、ρ_i、A_i、E_i。由测站坐标可计算出测站东北天坐标系到J2000.0地心惯性系转换矩阵为M_ti，则可计算出各点在测站东北天坐标系下的位置矢量为

R_di＝[ρ_icosE_isinA_i ρ_icosE_icosA_i ρ_isinE_i] (i＝0,1,2…)

则各点在J2000.0地心惯性坐标系下位置矢量为R_Ei＝M_tiR_di ^T。

进而对得到的位置数列，使用最小二乘曲线拟合的方法可计算出位置初值、速度初值。而加速度等其它量的初值可设置为0。

3、进行模型补偿的返回弹道UKF滤波计算

该模型补偿的返回弹道UKF滤波计算包括以下步骤：

(a)UKF采样

基于系统状态变量与协方差矩阵P对系统状态量进行UKF采样，此处可取α＝0.5，β＝2，κ＝0超维球比例对称采样；

对k-1时刻状态量的采样点计算公式为

式中：λ＝α²(n+κ)-n，n为状态向量维数；

(b)时间更新计算：

状态预测均值为：

状态协方差为：

其中，W_X为状态均值变换权系数，W_P为协方差变换权系数，Q为状态协方差矩阵。

此步骤中，对Q矩阵进行自适应处理，处理方法为：

对于UKF滤波过程，预测均值的速度项x、y、z三个方向上的分量简记为x_vi，位置项三分量为x_pi、状态更新值/>的速度项的x、y、z三个方向上的分量简记为x_vi’，位置项记为x_pi’，

则速度模型误差量近似值为：

σ_vi＝((x'_pi-x_pi)/Δt) (i＝x,y,z)

加速度模型误差近似值为：

σ_ai＝((x'_vi-x_vi)/Δt) (i＝x,y,z)

进而，Q矩阵自适应更新为：

式中，λ为弹道机动频率系数，0＜λ＜1，当机动加速度小时，或需要提高抗差效果时，λ应取小量，例如取0.001；可给定为较小的固定值(例如取为0)。

这样就实现了Q的自适应更新。

(c)计算观测量预测值(i＝0，1，…，2n)：

计算方法为：

设某个测站在J2000.0地心惯性系下的位置、速度矢量分别记为R,(由测站坐标可算出)；航天器位置、速度矢量分别为r＝(x,y,z)、/>地心惯性系到测站东北天坐标系的转换矩阵为M_I。

则航天器在测站东北天坐标系下的位置矢量为M_I(r-R)，将该矢量的三个分量分别标记为ρ_x、ρ_y、ρ_y。

则可计算出测距值为：

方位角为：

A＝arctan(ρ_x/ρ_y)

仰角为：

测距变化率为：这样就建立了滤波计算的观测模型。

(d)计算观测预测均值与方差计算：

观测预测均值为：

观测预测协方差为：

其中，R为事先给定的观测噪声协方差矩阵。

(e)计算最优滤波增益K_k，并更新系统状态和协方差矩阵/>

滤波增益为：

系统状态更新公式为：

协方差更新公式为：

以某任务返回段理论弹道为基础，仿真外测观测数据，并为其加行高斯白噪声偏差，使用本文算法进行了计算，计算结果见图2。

图2中浅色线为理论弹道计算的高度曲线，深色线为对仿真外测数据，采用本发明方法计算的高度曲线。从图中可见，本方法计算结果与理论值基本一致，只有个别点偏差较大。由此可见，本发明的算法是有效的。

Claims (4)

1.一种模型自适应补偿的返回弹道估计方法，其特征在于，具体按照以下步骤实施：

步骤1、使用UKF滤波算法框架进行计算，在J2000.0地心惯性系下建立十维的系统状态模型；

所述步骤1具体如下：

步骤1.1、在J2000.0地心惯性系下建立系统状态向量为：

其中，x,y,z分别表示为航天器位置矢量的x分量、y分量、z分量；分别为速度矢量的x分量、y分量、z分量；a_x,a_y,a_z分别为未建模力即机动推力产生的加速度的x分量、y分量、z分量；β为单位加速度的加速度变化率；

步骤1.2、建立十维的系统状态外推模型为：

其中，g_m为返回段已建模力产生的加速度矢量，加速度矢量为x、y、z各方向的三维分量分别标记为g_mx、g_my、g_mz，已建模力包括地球引力、太阳引力、月亮引力、大气阻力；

步骤1.3、对所述步骤1.2中的公式(1)进行数值积分，得到滤波计算的状态估计模型；

步骤1.4、观测模型建立在测站东北天坐标系下：

设某个测站在J2000.0地心惯性系下的位置、速度矢量分别记为R,由测站坐标可算出；航天器位置、速度矢量分别为r＝(x,y,z)、/>地心惯性系到测站东北天坐标系的转换矩阵为M_I；

则航天器在测站东北天坐标系下的位置矢量为M_I(r-R)，将该矢量的三个分量分别标记为ρ_x、ρ_y、ρ_z，

则可计算出测距值为：

方位角为：

A＝arctan(ρ_x/ρ_y) (3)

仰角为：

测距变化率为：

上述公式(2)～(5)即为建立的滤波计算的观测模型；

步骤2、对滤波中的状态噪声协方差矩阵进行模型自适应补偿处理，实现滤波对弹道机动过程的适应。

2.根据权利要求1所述的一种模型自适应补偿的返回弹道估计方法，其特征在于，所述步骤2具体如下：

对UKF滤波过程，预测均值为将预测均值/>的速度项x、y、z三个方向上的分量标记为x_vi其中，i＝x,y,z，位置项的三分量标记为x_pi、状态更新值为/>将速度项的x、y、z三个方向上的分量标记为x_vi'，位置项三分量标记为x_pi'；

则可计算出速度模型误差量的近似值为：

σ_vi＝((x'_pi-x_pi)/Δt)

加速度模型误差近似值为：

σ_ai＝((x'_vi-x_vi)/Δt)

进而，Q矩阵自适应更新为：

λ为弹道机动频率系数；

至此，模型自适应补偿完成。

3.根据权利要求2所述的一种模型自适应补偿的返回弹道估计方法，其特征在于，所述步骤2中弹道机动频率系数λ为0.001。

4.根据权利要求3所述的一种模型自适应补偿的返回弹道估计方法，其特征在于，所述步骤2中取值均为0。