CN109781374A - 一种实时在线快速估计飞行器推力的方法 - Google Patents

Abstract

本发明涉及一种实时在线快速估计飞行器推力的方法，不依靠飞行器弹载传感器设备信息，仅依靠外测数据处理系统，通过构建飞行器的混合动力学模型，用物理模型描述飞行器的地球引力和大气阻力等能精确建模的部分，用数学模型描述不能精确建模的推力部分，并通过无迹滤波（UKF）的估计算法，实现了飞行器推力的在线实时估计，为飞行器试验鉴定、精度评估提供非常重要的支撑。

Description

一种实时在线快速估计飞行器推力的方法

技术领域

本发明适用于航空航天测控领域，尤其适用于导弹、卫星的测控系统。

背景技术

飞行器推力是保持飞行器按照设计路线前进的主要动力，如何在飞行器飞行过程中实时估计飞行器推力值是试验鉴定、飞行安控等飞行状态评估的重要手段。尤其是对于非合作飞行器而言，比如多级导弹，实时估计其推力值可以快速获得加速度曲线，从而根据加速度模板确定导弹型号，对于导弹预警、防空反导具有极为重要的意义。在线快速估计飞行器推力是不依靠飞行器弹载陀螺仪等遥测传感器，仅根据地面雷达站的测控数据构建观测方程，并根据目标运动特性构建状态方程，建立实时滤波系统，在线输出目标的轨迹和飞行器推力的过程。由于飞行需求的多样性，很多飞行器需要完成复杂的机动特性，推力控制往往很难构建确定的动力学模型，不仅造成弹道确定的复杂，更造成了推力计算的高难度，特别对于非合作飞行器目标，无法通过遥测信息获取推力曲线，只能通过地面的雷达、光学设备等外测设备获取目标的相对位置、角度、速度等外测信息。如何通过外测信息在线估计飞行器的推力特性对靶场测控、预警安全具有极高的应用价值。

通常的做法是构建带有参数的推力动力学模型，利用外测数据的测量信息估计该参数，并预测下一时刻的参数，然后再根据下一时刻的外测数据修正参数，最后达到收敛，并输出推力数据。该方法适合于推力特性特别简单的形式，且已经知道目标的飞行特性。对于非合作飞行器而言，无法构建推力模型。另外及时对于合作飞行器，当面临高机动特性时，状态方程建模很难构建推力动力学模型，通常采用描述目标运动趋势的运动学模型作为状态方程，无法实时输出飞行器的推力信息，主要通过差分的方式获取目标的加速度信息，推力估计精度收到很大影响。

因此，目前针对在线实时估计飞行器推力的问题，尚未见到合适的综合处理方法。本发明将提供一种不依赖弹载遥测传感器仅利用外测信息实时在线估计飞行器推力的方法，对于飞行器的试验鉴定、飞行评估、安全监控等都具有非常实用的价值。

发明内容

在飞行器进行高速机动时，无法准确构建飞行器的动力学方程，造成飞行器在线估计困难，并无法输出推力曲线。本发明可有效解决该类问题，不依赖飞行器自身的传感器设备，仅依靠地面测控数据，就可以在线输出飞行器的推力曲线，并且可以提高在线估计飞行器弹道的精度。针对现有技术的不足，本发明提供一种能根据外测信息在线估计飞行器推力的方法，使得不依赖飞行器自身传感器信息，仅依靠外测数据就能估计飞行器推力。

本发明的具体技术方案如下：

一种实时在线快速估计飞行器推力的方法，包括以下步骤：

S1、建立外测数据的测量方程；

本发明利用的外测数据处理设备为连续波雷达，对于其它类型的雷达设备或者光学设备同样适应。相关模型如下。

雷达设备测量数据，与地球固连，测量数据包含斜距、斜距变化率、方位角和俯仰角，模型描述为，

其中，X_b(t)＝(x_b(t),y_b(t),z_b(t))^T,V_xb＝(v_xb(t),v_yb(t),v_zb(t))^T为飞行器在测站坐标系下的位置和速度分量，ε_R(t),ε_α(t),ε_β(t),ε_v(t)分别为斜距、方位角、俯仰角、和斜距变化率的测量模型随机误差。

综上，将测量模型记为：

H_i(t)＝h_i(X(t),V_x(t))+ε,i＝1,2,3…,n (1)

其中H为t时刻第i个雷达设备获取的测量数据，数据格式为X(t)为惯性坐标系下的飞行器位置矢量，V_x(t)为惯性坐标系下的飞行器速度矢量；惯性坐标系和测站坐标系可以由坐标系转换矩阵完成转换。

S2、构建带推力模型的飞行器混合动力学模型：

混合动力学模型主要由两部分构成，一部分是飞行器地球引力、大气阻力的动力学模型；另一部分是构建的飞行器推力运动学模型；

S201：构建不含推力的飞行器动力学方程

根据牛顿第二定律飞行器在惯性坐标系中的运动方程为：

式中，为地球非球形引力；

为大气阻力；

地球潮汐(包括固体潮、海潮和大气潮汐)使地球对卫星引力的变化部分；

相对论效应对卫星的影响；

太阳辐射压力；

μ地球引力常数。

对于一般精度的飞行器而言，摄动力可以只考虑到地球非球形引力和大气阻力；

1.地球非球形引力。

上式中，分别为卫星的地心矢径、地心经度、地心纬度，P_n(x)为勒让德多项式，缔合勒让德多项式，J_n为带谐系数，A_nm,B_nm为田谐系数。

具体地球非球形引力摄动对卫星轨道的影响以及上式中各个量的计算见文献。现在比较常用的引力场模型有：JGM-2、JGM-3、GRIM4-S4、EGM96、GRIM5-S1、GEM-T3等。

2.大气阻力。

由于地球大气的存在，大气阻力摄动也就成为了影响飞行器运动的重要因素。但是由于地球大气组成的多样性，再加上它们的分布随着高度、纬度和太阳照射条件的不同而使得大气密度分布出现了相当的不规则性和随机性，从而使得飞行器在大气中飞行的受力情况也相当复杂。

为简单起见，建立大气静止情况下的大气对卫星运动的阻力，即

上式中，S为飞行器的横截面积，c_d为阻力系数，m为飞行器质量，ρ为飞行器所在位置的大气密度，本发明选取最常用的常标高球面大气模式，ρ＝ρ₀exp(r₀-r/H)式中，ρ₀为r₀处的大气密度，r₀为轨道近地点的矢径，H为密度标高，假定H为常数，r为飞行器矢径。为飞行器相对于大气的速度，

其中，ω_d为大气旋转速度，一般认为低层大气旋转速度为地球自转速度ω_e。

现在比较常见的大气模型有：CIRA86、MSIS90、JACCHIA系列、DTM大气模型等。

根据飞行器地球非球形引力和大气阻力模型，构建飞行器动力学状态方程

设表示为大气阻力在x,y,z三方向引起的摄动加速度，则有：

简记状态方程为：

其中，状态向量X＝(x，y，z，v_x，v_y，v_z)^T代表J2000.0坐标下飞行器的飞行状态，μ为引力常数，R_e为地球半径。

S202：构建的飞行器推力运动学模型

飞行器运动的加速度变化具有明显的时序相关特性，因此根据飞行器加速度的变化特点建立反映这种相关性的模型，设推力a_u(t)满足微分方程的一阶Gauss-Markov过程，

其中A(t)和B(t)为系数矩阵，u(t)是高斯噪声向量，其分量满足，

E(u(t))＝0,E(u(t)u^T(t))＝Iδ(t-τ)

其中，T₁,T₂,T₃为相关时间。

其中，b_j(j＝1,2,3)作为指定常数看待。

S203：混合状态模型

飞行器动力学模型描述了飞行器在各种复杂力的作用下的客观现实，推力运动学模型则描述了推力综合作用下的状态表征，将两种模型结合起来，可以组成混合状态模型。

飞行器的运动状态方程为：

具体可描述为，

对状态向量进行扩维，定义新的状态变量，可得

由此获得新的状态方程如下：

根据测量方程(1)，构建在线滤波系统如下：

上述(9)式即为本发明中建立的带有飞行器推力估计的跟踪测量系统，待估计的状态变量包括惯性坐标系下的飞行器位置和速度向量X(t)、推力估计值a_u(t)；

S3、基于无迹滤波(Unscented Kalman Filter,UKF)的在线估计算法：

根据构建的系统状态方程和测量方程(1)可见，测量模型和状态方程都是复杂的非线性模型，因此，本发明中，在线估计算法采用适用于非线性系统模型的无迹滤波(UKF)；算法的步骤如下：

S301、初始点设置及相应的采样点选取；

初始状态的均值x₀和协方差P_xx，对称采样点为：

其中，n_D表示待估计的状态变量的维数，λ＝α²(n+κ)-n_D，κ是一个比例参数，通常设置为0或者3-n_D，α设为一个较小的正数(1≥α≥0.0001)，对应的权值为：

其中W_i ^m为均值加权所用的权值，W_i ^c为协方差加权所用的权值，通常β＝1+α²，对于高斯分布，通常取β＝2；

S302、预测下一时刻点；

将公式(9)中状态方程和测量模型的非线性变换应用于k时刻采样的每个采样点，得到非线性转换后的点集χ_i(k+1|k)，进而利用转换后的点集计算其均值和协方差用于估计：

状态采样点下一时刻的状态预测:

χ_i(k+1|k)＝f[χ_i(k|k)] (13)

相应的状态预测协方差矩阵计算:

状态采样点下一时刻测量预测:

H_mi(k+1|k)＝h[χ_i(k+1|k)] (16)

相应的测量预测协方差矩阵以及测量与状态之间的协方差矩阵计算:

S303、更新预测值；

利用测量更新状态预测(14)的增益矩阵计算：

状态更新：

状态更新之后相应的状态协方差矩阵计算：

P_xx(k+1|k+1)＝P(k+1|k)-K(k+1)P_yy(k+1|k)K^T(k+1) (22)

即为初始时刻之后下一时刻最终的估计值；然后重复S301-S303步骤，即得到各个时刻状态变量的估计值。

本发明的有益效果是：

1)考虑到飞行器推力动力学建模的复杂性，从推力加速度的运动特性着手，构建了推力的数学模型和能够建模的动力学模型融合在一起构建了混合状态方程，提高了S302中状态预测的精度；

2)利用混合状态方程和外部测控数据的测量方程，可以实现推力的在线估计，并且可以提高S303中飞行器弹道的估计值精度。

附图说明

图1为实施例中在仿真条件下估计推力和实际推力误差差图，

图2为实施例中在仿真条件下推力估计值图，

图3为实施例中在仿真条件下估计的飞行器位置误差图，

图4为实施例中在仿真条件下估计的飞行器速度误差图，

图5为实施例中仿真的飞行器速度曲线图，

图6为实施例中仿真的飞行器位置曲线图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

给出仿真参数如下：

1)地面雷达站四个，测量随机误差标准差为:斜距1m，方位角20角秒，俯仰角20角秒，斜距变化率1cm/s；随机误差服从正态分布。

2)飞行器为单级火箭加速，在200秒时刻推力结束，推力仿真考虑质量变化因素；

3)飞行器弹道仿真考虑地球非球形引力J2阶次，大气阻力摄动；

4)初始位置误差为：[1km,1km,1km]，初始速度误差为：[5m/s,5m/s,5m/s]；

5)外测采样频率1Hz；

S1、建立外测雷达的测量模型：

首先仿真飞行器在惯性坐标系下的状态数据，并根据坐标系转换公式，求得飞行器在每个测站坐标系下的坐标。根据测量模型(1),获得每个雷达的测量数据；

S2、飞行器混合动力学模型构建：

飞行器的动力学模型按照公式(2)计算，其中地球非球形引力仅考虑J2项，其余摄动仅考虑大气阻力摄动，其中的μ＝398600441500000；

S3、基于无迹滤波(Unscented Kalman Filter,UKF)的在线估计算法：

推力初始估计值取为：[0，0，0]；

初始位置误差为：[1km,1km,1km]，初始速度误差为：[5m/s,5m/s,5m/s]；

待估计的状态变量的维数n_D＝10，κ＝3-n_D，α＝1，β＝2，λ＝-9.97，γ＝0.173，进而按照UKF滤波计算流程进行各个时刻状态估计；

基于建立的系统模型(7)和(9)，当第1至200秒飞行器在推力作用下飞行时，系统可以实时估计出飞行器的飞行推力，推力的实时估计精度和推力结果如图1和图2，飞行器轨迹结果如图3和图4所示；飞行器位置和速度曲线图如图5和图6所示。

根据上述实验结果可见，对于机动飞行器目标而言，推力变化无法用确定模型构建状态方程，本发明采用的混合动力学模型，用物理模型描述地球引力和大气阻力的运动趋势，用数学模型描述推力的运动趋势，由此构成的混合动力学模型可以很好的描述目标的运动特性，因此在估计过程中，可以非常准确的估计出飞行器的推力，绘制推力曲线，为飞行器飞行鉴定和试验评估提供非常有价值的数据支撑。

Claims (4)

1.一种实时在线快速估计飞行器推力的方法，不依赖飞行器自身传感器信息，仅依靠外测数据就能估计飞行器推力，其特征在于，包括以下步骤：

S1、建立外测数据的测量方程模型；

S2、构建带推力模型的的飞行器混合动力学模型；

混合动力学模型包括两部分，一部分是飞行器地球引力、大气阻力的动力学模型；另一部分是构建的飞行器的推力估计模型；

S201：构建不含推力的飞行器动力学方程；

S202：构建的飞行器的推力估计模型

S203：混合状态模型；

S3、基于无迹滤波在线估计算法，

S301、初始点设置及相应的采样点选取；

S302、预测下一时刻点；

S303、更新预测值。

2.根据权利要求1所述的一种实时在线快速估计飞行器推力的方法，其特征在于，所述步骤S1建立外测数据的测量方程模型，具体步骤如下：

雷达设备测量数据，与地球固连，测量数据包含斜距、斜距变化率、方位角和俯仰角，模型描述为，

其中，X_b(t)＝(x_b(t),y_b(t),z_b(t))^T,为飞行器在测站坐标系下的位置和速度分量，ε_R(t),ε_α(t),ε_β(t),ε_v(t)分别为斜距、方位角、俯仰角、和斜距变化率的测量模型随机误差，

综上，将测量模型记为：

H_i(t)＝h_i(X(t),V_x(t))+ε,i＝1,2,3…,n (1)

其中H为t时刻第i个雷达设备获取的测量数据，数据格式为为惯性坐标系下的飞行器位置矢量，V_x(t)惯性坐标系下的飞行器速度矢量。

3.根据权利要求1所述的一种实时在线快速估计飞行器推力的方法，其特征在于，所述步骤S2构建带推力模型的的飞行器混合动力学模型，具体步骤如下：

混合动力学模型包括两部分，一部分是飞行器地球引力、大气阻力的动力学模型；另一部分是构建的飞行器的推力估计模型，

S201：构建不含推力的飞行器动力学方程

根据牛顿第二定律飞行器在惯性坐标系中的运动方程为：

式中，为地球非球形引力，

为大气阻力，

地球潮汐使地球对卫星引力的变化部分，

相对论效应对卫星的影响，

太阳辐射压力，

μ地球引力常数，

对于一般精度的飞行器而言，摄动力只考虑到地球非球形引力和大气阻力，

(1)、地球非球形引力，

上式中，分别为卫星的地心矢径、地心经度、地心纬度，P_n(x)为勒让德多项式，缔合勒让德多项式，J_n为带谐系数，A_nm,B_nm为田谐系数，

(2)、大气阻力

简单起见，建立大气静止情况下的大气对卫星运动的阻力，即

上式中，S为飞行器的横截面积，c_d为阻力系数，m为飞行器质量，ρ为飞行器所在位置的大气密度，ρ＝ρ₀exp(r₀-r/H)式中，ρ₀为r₀处的大气密度，r₀为轨道近地点的矢径，H为密度标高，假定H为常数，r为飞行器矢径，为飞行器相对于大气的速度，

其中，ω_d为大气旋转速度，一般认为低层大气旋转速度为地球自转速度ω_e，

根据飞行器地球非球形引力和大气阻力模型，构建飞行器动力学状态方程

设表示为大气阻力在x,y,z三方向引起的摄动加速度，则有：

记状态方程为：

其中，状态向量X＝(x,y,z,v_x,v_y,v_z)^T代表J2000.0坐标下飞行器的飞行状态，μ为引力常数，R_e为地球半径，

S202：构建的飞行器推力运动学模型，

根据飞行器加速度的变化特点建立反映这种相关性的模型，设推力a_u(t)满足微分方程的一阶Gauss-Markov过程，

其中A(t)和B(t)为系数矩阵，u(t)是高斯噪声向量，其分量满足，

E(u(t))＝0,E(u(t)u^T(t))＝Iδ(t-τ)

其中，T₁,T₂,T₃为相关时间，

其中，b_j(j＝1,2,3)作为指定常数，

S203：混合状态模型

飞行器动力学模型描述飞行器在各种复杂力的作用下的客观现实，推力运动学模型则描述了推力综合作用下的状态表征，将两种模型结合起来，组成混合状态模型，

飞行器的运动状态方程为：

具体描述为，

对状态向量进行扩维，定义新的状态变量，可得

由此获得新的状态方程如下：

根据测量方程(1)，构建在线滤波系统如下：

上述(9)式即为本发明中建立的带有飞行器推力估计的跟踪测量系统，待估计的状态变量包括惯性坐标系下的飞行器位置和速度向量X(t)、推力估计值a_u(t)。

4.根据权利要求1所述的一种实时在线快速估计飞行器推力的方法，其特征在于，所述步骤S3基于无迹滤波在线估计算法，具体步骤如下：

S301、初始点设置及相应的采样点选取：

初始状态的均值x₀和协方差P_xx，对称采样点为：

其中，n_D表示待估计的状态变量的维数，λ＝α²(n+κ)-n_D，κ是一个比例参数，通常设置为0或者3-n_D，α设为一个较小的正数，1≥α≥0.0001，对应的权值为：

其中W_i ^m为均值加权所用的权值，W_i ^c为协方差加权所用的权值，通常β＝1+α²，对于高斯分布，取β＝2；

S302、预测下一时刻点；

将公式(9)中状态方程和测量模型的非线性变换应用于k时刻采样的每个采样点，得到非线性转换后的点集χ_i(k+1|k)，进而利用转换后的点集计算其均值和协方差用于估计：

状态采样点下一时刻的状态预测:

χ_i(k+1|k)＝f[χ_i(k|k)] (13)

相应的状态预测协方差矩阵计算:

状态采样点下一时刻测量预测:

H_mi(k+1|k)＝h[χ_i(k+1|k)] (16)

相应的测量预测协方差矩阵以及测量与状态之间的协方差矩阵计算:

S303、更新预测值：

利用测量更新状态预测(14)的增益矩阵计算：

状态更新：

状态更新之后相应的状态协方差矩阵计算：

P_xx(k+1|k+1)＝P(k+1|k)-K(k+1)P_yy(k+1|k)K^T(k+1) (22)

即为初始时刻之后下一时刻最终的估计值；然后重复以上步骤，即得到各个时刻状态变量的估计值。