CN115342815B - 反大气层内或临近空间机动目标视线角速率估计方法 - Google Patents

Abstract

反大气层内或临近空间机动目标视线角速率估计方法，它属于导弹制导控制技术领域。本发明解决了由于目标加速度的估计结果不精确以及拦截导弹弹体的冲击振动，导致对目标与导弹视线角速率估计的精度低的问题。本发明基于机动目标跟踪滤波器跟踪得到目标加速度信息，再结合目标加速度信息和冲击振动带来的扰动信息，采用无迹卡尔曼粒子滤波算法对视线角速率进行滤波估计，克服了由于现有的目标与导弹角速度估计方法中目标加速度的不精确以及弹体的冲击振动所带来的视线角速率估计精度低的问题，进而根据本发明的视线角速率的估计结果提高制导精度。本发明方法可以应用于导弹制导控制技术领域。

Description

反大气层内或临近空间机动目标视线角速率估计方法

技术领域

本发明属于导弹制导控制技术领域，具体涉及反大气层内或临近空间机动目标视线角速率估计方法。

背景技术

在制导技术领域，目标的持续机动在很大程度上增加了视线角速率高精度提取的难度，现有的视线角速率提取方法不能满足拦截导弹或拦截器直接碰撞目标的精确制导要求，必须研究与目标机动条件相适应的新的高精度视线角速率提取方法。

近些年来，为了降低其成本，提高其可靠性，精确制导武器越来越多地采用小型捷联或半捷联导引头，这种导引头的特点是可以测量目标—导弹之间的相对视线角，如果是主动式导引头还可以提供距离信息，但捷联或半捷联导引头不能直接输出实际制导律需要的视线角速率信息。

针对安装捷联或半捷联导引头的导弹拦截处于大气层内或临近空间的机动目标这一应用背景，由于目标加速度的估计结果不精确以及拦截导弹弹体的冲击振动，导致对目标-导弹视线角速率估计的精度也较低，目前还没有一种行之有效的视线角速率高精度提取方法。

发明内容

本发明的目的是为解决由于目标加速度的估计结果不精确以及拦截导弹弹体的冲击振动，导致对目标与导弹视线角速率估计的精度低的问题，而提出的反大气层内或临近空间机动目标视线角速率估计方法。

本发明为解决上述技术问题所采取的技术方案是：

反大气层内或临近空间机动目标视线角速率估计方法，所述方法具体包括以下步骤：

步骤一、测算导弹加速度信息；

步骤二、基于十二阶机动目标模型建立跟踪滤波器，利用建立的跟踪滤波器跟踪目标加速度信息；

步骤三、基于导弹加速度信息和目标加速度信息，并考虑弹体冲击振动的影响，对目标与导弹视线角速率进行估计。

进一步地，所述步骤一的具体过程为：

建立视线坐标系o′x₄y₄z₄：原点o′位于导引头回转中心；o′x₄轴与目标-导弹视线一致，由导引头回转中心指向目标方向为正；o′y₄轴位于包含o′x₄轴的铅垂面内，与o′x₄轴垂直且以指向上方为正；o′z₄轴按右手定则确定；

将导弹上加速度计输出的加速度在视线坐标系的o′y₄轴方向上的投影记为a_ε，将导弹上加速度计输出的加速度/>在视线坐标系的o′z₄轴方向上的投影记为a_β，a_ε和a_β即为导弹法向加速度。

进一步地，所述步骤二的具体过程为：

建立速度坐标系O₃x₅y₅z₅：将原点O₃定义在目标飞行器的质心上，O₃x₅轴沿着目标飞行器速度方向，O₃y₅轴定义在包含目标飞行器纵轴的纵向对称平面内，并与O₃x₅轴垂直，O₃z₅轴由右手定则确定；

建立弹道坐标系O₃x₆y₆z₆：O₃x₆轴沿着目标飞行器速度方向，O₃y₆轴定义在包含速度矢量的铅垂面内，并与O₃x₆轴垂直，O₃z₆轴由右手定则确定；

将作用在目标飞行器上的总的空气动力沿着速度坐标系分为阻力、升力和侧向力：

其中，X为阻力，Y为升力，Z为侧向力，q为动压头，c_x是无量纲的阻力系数，c_y是无量纲的升力系数，c_z为无量纲的侧向力系数，s为目标飞行器的特征面积；

q＝ρv²/2 (2)

其中，ρ代表空气密度，v代表目标飞行器的飞行速度；

由空气动力提供的速度坐标系下的O₃x₅轴、O₃y₅轴、O₃z₅轴的加速度为：

其中，a_X为空气动力提供的速度坐标系下的O₃x₅轴方向加速度，a_Y为空气动力提供的速度坐标系下的O₃y₅轴方向加速度，a_Z为空气动力提供的速度坐标系下的O₃z₅轴方向加速度，m是目标飞行器的质量，v_x是目标飞行器在速度坐标系O₃x₅轴方向的飞行速度，v_y是目标飞行器在速度坐标系O₃y₅轴方向的飞行速度，v_z是目标飞行器在速度坐标系O₃z₅轴方向的飞行速度；

利用弹道坐标系来代替速度坐标系，将速度坐标系下的由空气动力引起的加速度转入发射点惯性坐标系下：

其中，a_tx、a_ty、a_tz为发射点惯性坐标系下的加速度，C_1→4是发射点惯性坐标系与弹道坐标系之间的转换矩阵，θ代表弹道倾角(即O₃x₆轴与水平面O_Fx_Fz_F之间的夹角)，当拦截导弹的速度矢量位于O_Fx_Fz_F面上方时，θ为正；ψ_v代表弹道偏角(即O₃x₆轴映在O_Fx_Fz_F面上的投影与O_Fx_F轴之间的夹角)，迎着O_Fy_F轴方向从上往下看，若投影转到O_Fx_F轴为顺时针方向时，ψ_V为正；

根据万有引力定律得：

其中，g₀代表地球对目标飞行器的引力，r为目标飞行器与地心的距离，μ＝GM，G为引力常量，M为地球质量，R_e为地球半径，x、y和z是目标飞行器在发射点惯性坐标系下的坐标；

联合公式(4)和公式(5)，求得作用在目标飞行器上的发射点惯性坐标系下的总的加速度为：

其中，a_x、a_y、a_z分别代表作用在目标飞行器上的发射点惯性坐标系下x轴、y轴、z轴方向的总的加速度；

引入与气动参数有关的状态量Z_x、Z_y、Z_z：

则将目标飞行器的总加速度表示为：

再取状态量k_x、k_y、k_z分别作为状态量Z_x、Z_y、Z_z的导数，选取目标飞行器的位置信息和速度信息作为状态量，则得到的十二阶机动目标模型为：

其中，为x的一阶导数，/>为v_x的一阶导数，/>为Z_x的一阶导数，/>为k_x的一阶导数，/>为y的一阶导数，/>为v_y的一阶导数，/>为Z_y的一阶导数，/>为k_y的一阶导数，/>为z的一阶导数，/>为v_z的一阶导数，/>为Z_z的一阶导数，/>为k_z的一阶导数，w_n是零均值的高斯白噪声；

跟踪滤波器的测量信息由基站(地面的目标跟踪装置)测得的位置信息提供，则跟踪滤波器的测量方程为：

其中，z₁、z₂、z₃是基站测得的目标飞行器在发射点惯性系下的三个轴上的位置信息；

选取跟踪滤波器的状态向量X为：

X＝[x v_x Z_x k_x y v_y Z_y k_y z v_z Z_z k_z]^T (11)

将公式(11)代入公式(9)中，得到跟踪滤波器的状态方程为：

其中，x′₁＝x，x′₂＝v_x，x′₃＝Z_x，x′₄＝k_x，x′₅＝y，x′₆＝v_y，x′₇＝Z_y，x′₈＝k_y，x′₉＝z，x′₁₀＝v_z，x′₁₁＝Z_z，x′₁₂＝k_z；

将公式(12)的连续状态方程转变为离散的状态方程，再对十二阶机动目标模型采用扩展卡尔曼滤波算法设计跟踪滤波器；

根据设计的跟踪滤波器对目标加速度a_tx、a_ty和a_tz进行估计，将目标加速度a_tx、a_ty和a_tz的估计值从发射点惯性坐标系投影到视线坐标系后，得到目标加速度a_tε和a_tβ。

进一步地，所述目标飞行器的特征面积s为：

对于有翼式目标飞行器，目标飞行器的特征面积s为弹翼面积；对于无翼式目标飞行器，目标飞行器的特征面积s为弹身最大横截面的面积。

进一步地，所述将公式(12)的连续状态方程转变为离散的状态方程采用的是欧拉法。

进一步地，所述步骤三的具体过程为：

目标与导弹之间的相对距离矢量利用目标的距离矢量减去导弹的距离矢量得到：

其中，为目标与导弹之间的相对距离矢量，/>为发射点惯性坐标系下目标与原点之间的距离矢量，/>为发射点惯性坐标系下导弹与原点之间的距离矢量；

则目标与导弹之间的相对距离矢量的绝对倒数与相对导数之间的关系为：

其中，为目标与导弹之间的相对距离矢量在发射点惯性坐标系下对时间的绝对导数，/>为目标与导弹之间的相对距离矢量在视线坐标系下对时间的相对导数，/>为视线坐标系相对于发射点惯性坐标系的旋转角速度；

将公式(14)表达为公式(15)的矩阵形式：

其中，r表示目标与导弹距离，表示目标与导弹距离的导数；

其中，和/>分别为目标和导弹在视线坐标系下的速度矢量，且

其中，q_ε为视线俯仰角，q_β为视线偏航角，为q_ε的一阶导数，/>为q_β的一阶导数；

将公式(16)代入公式(15)得：

令矢量为：

对公式(18)两边相对时间求导得：

其中，a_t和a分别为目标和导弹在视线坐标系下的加速度；

根据矢量的绝对导数与相对导数之间的关系，将矢量代入公式(19)得：

其中，通过结合公式(17)和公式(18)求得：

其中，为/>一阶导数，/>为/>一阶导数，/>为/>一阶导数；

结合公式(16)、公式(20)和公式(21)得：

将公式(15)代入公式(22)中，得到目标与导弹之间相对运动的动力学方程为：

其中，a_r、a_ε、a_β分别是拦截导弹在视线坐标系下o′x₄、o′y₄、o′z₄三个轴方向上的加速度信息，a_tr、a_tε、a_tβ分别是机动目标在视线坐标系下o′x₄、o′y₄、o′z₄三个轴方向上的加速度信息；

根据公式(23)的动力学方程选取状态向量X′为：

将公式(24)代入公式(23)中，得到制导滤波器的状态方程为：

其中，x″₁＝q_ε，x″₃＝q_β，/>

选取导引头测得的视线俯仰角和视线偏航角作为测量值，得到制导滤波器的测量方程为：

其中，y₁为导引头测得的视线俯仰角的测量值，y₂为导引头测得的视线偏航角的测量值；

对公式(25)的状态方程进行离散化，得：

其中，T为采样间隔，x″₁(k)为x″₁在k时刻的值，x″₂(k)为x″₂在k时刻的值，x″₃(k)为x″₃在k时刻的值，x″₄(k)为x″₄在k时刻的值；

考虑模型误差和弹体冲击振动影响的制导滤波器的测量方程为：

z(k)＝y(k)+e(k)+v(k) (28)

其中，z(k)为考虑模型误差和弹体冲击振动影响的k时刻的量测值，y(k)为不考虑模型误差和弹体冲击振动影响的k时刻的量测值，v(k)为k时刻的模型误差，e(k)为k时刻由弹体冲击引起的在导引头处的挠性变形信息；

结合公式(27)和公式(28)得到制导滤波器的模型：

其中，W代表制导滤波器的过程噪声，V为量测的系统性偏差与测量噪声之和，e为由弹体冲击引起的在导引头处的挠性变形信息，f(X′)为状态转移矩阵，H为量测转移矩阵，Z为目标量测；

基于公式(29)的制导滤波器模型和无迹卡尔曼粒子滤波算法对目标与导弹视线俯仰角速率和视线偏航角速率进行估计。

本发明的有益效果是：

本发明基于机动目标跟踪滤波器跟踪得到目标加速度信息，再结合目标加速度信息和冲击振动带来的扰动信息，采用无迹卡尔曼粒子滤波算法对视线角速率进行滤波估计，采用本发明方法可以精确估计出视线角速率，克服了由于现有的目标与导弹角速度估计方法中目标加速度的不精确以及弹体的冲击振动所带来的视线角速率估计精度低的问题，进而根据本发明的视线角速率的估计结果提高制导精度。

附图说明

图1是发射点惯性坐标系与地心惯性坐标系之间关系的示意图；

图2是发射点惯性坐标系和视线坐标系之间关系的示意图；

图3是导弹弹体坐标系示意图；

图4a是x轴加速度跟踪结果图；

图4b是y轴加速度跟踪结果图；

图4c是z轴加速度跟踪结果图；

图5a是九阶机动目标模型与十二阶目标机动目标模型的位置误差对比图；

图中，ekf9_x代表九阶机动目标模型，ekf_x代表十二阶目标机动目标模型；

图5b是九阶机动目标模型与十二阶目标机动目标模型的速度误差对比图；

图5c是九阶机动目标模型与十二阶目标机动目标模型的加速度误差对比图；

图6是UPF视线俯仰角速率的滤波估计结果图；

图7是UPF视线偏航角速率的滤波估计结果图；

图8是UPF视线俯仰角速率滤波误差图；

图9是UPF视线偏航角速率滤波误差图；

图10是四种制导滤波器的视线俯仰角速率的滤波估计误差对比图；

图11是四种制导滤波器的视线偏航角速率的滤波估计误差对比图；

图12是目标与导弹相对飞行轨迹图；

图13是目标与导弹的高度变化图；

图14a是目标在飞行过程中的速度曲线图；

图14b是目标在飞行过程中的弹道倾角曲线图；

图14c是目标在飞行过程中的弹道偏角曲线图；

图15a是目标在飞行过程中的速度曲线图；

图15b是目标在飞行过程中的弹道倾角曲线图；

图15c是目标在飞行过程中的弹道偏角曲线图；

图16是导弹飞行角速度的示意图；

图17a是导弹飞行过程中末制导环节的俯仰角变化曲线图；

图17b是导弹飞行过程中末制导环节的偏航角变化曲线图；

图17c是导弹飞行过程中末制导环节的滚转角变化曲线图；

图18是末制导拦截的脱靶量变化情况图；

图19是末制导拦截的燃料消耗情况图。

具体实施方式

具体实施方式一、本实施方式提供了一种针对气动力精确模型的反大气层内或临近空间机动目标视线角速率估计方法，以克服已有的目标与导弹角速度估计方法中目标加速度的不精确以及弹体的冲击振动带来的视线角速率估计精度低的缺陷，从而提高了制导精度。

本发明通过下述步骤实现：

1、测算目标与导弹之间的相对距离和相对速度，具体过程如下：

地心惯性坐标系O_Ix_Iy_Iz_I：设原点O_I位于地球地心处，O_Iy_I轴为指向发射点的方向，O_Ix_I轴位于由地心、发射点、目标点构成的弹道面与水平面的交线处，指向目标的方向为正，O_Iz_I轴方向通过右手定则来确定。

发射点惯性坐标系O_Fx_Fy_Fz_F：定义原点O_F位于发射点处，O_Fy_F轴位于铅垂线上，方向向上，O_Fx_F轴位于弹道面内垂直于O_Fy_F轴的方向上，正方向规定为指向目标的方向，O_Fz_F轴方向用右手定则来确定，发射点惯性坐标系与地心惯性坐标系的关系如图1所示。

导弹在发射点惯性坐标系的位置为(x、y、z)，其速度在该坐标系中的投影为(V_x、V_y、V_z)，这些信息由弹上导航系统提供；目标在发射点惯性坐标系的位置为(x_t、y_t、z_t)，其速度在该坐标系中的投影为(V_tx、V_ty、V_tz)，这些信息由弹上的目标跟踪滤波器提供；记Δx＝x-x_t，Δy＝y-y_t，Δz＝z-z_t，ΔV_x＝V_x-V_tx，ΔV_y＝V_y-V_ty，ΔV_z＝V_z-V_tz，则目标与导弹相对距离为相对速度为/>在末制导过程中，始终有R＞0，/>

2、测算导弹加速度

视线坐标系定义：视线坐标系o′x₄y₄z₄原点o′位于导引头回转中心,o′x₄轴与目标—导弹视线一致,由导引头回转中心指向目标为正，o′y₄轴位于包含o′x₄轴的铅锤面内，与o′x₄轴垂直，指向上方为正，o′z₄轴按右手定则确定。如图2所示。由发射点惯性坐标系转换到视线坐标系定义了视线俯仰角q_ε和视线偏航角q_β。

视线俯仰角q_ε：发射点惯性坐标系原点与目标飞行器质心的连线即视线轴与水平面之间的夹角，如果视线在水平面之上，则q_ε为正。

视线偏航角q_β：发射点惯性坐标系原点与目标飞行器质心的连线及视线轴在水平面Ox₀z₀面上的投影与Ox₀轴之间的夹角，迎着Oy₀轴方向从上往下看，若投影转到Ox₀轴为顺时针方向时，q_β为正。

弹体坐标系定义：由图3所示，导弹弹体坐标系ox₁y₁z₁的原点o位于导弹质心上。ox₁轴与弹体纵轴重合，指向头部为正，oy₁轴在弹体纵向对称平面内，垂直ox₁轴，指向上方为正，oz₁轴方向按右手定则确定。

测算导弹加速度在视线坐标系o′x₄y₄z₄的o′y₄和o′z₄轴方向上的分量a_ε和a_β。导弹上加速度计输出的加速度是沿着弹体坐标系ox₁y₁z₁投影的，把它投影到视线坐标系即可求得导弹法向加速度a_ε和a_β；

3、测算目标加速度

3.1基于九阶机动目标模型的跟踪滤波器

习惯上，一般将作用在导弹上的总的空气动力沿着速度坐标系可以分为阻力X、升力Y、侧向力Z，其计算公式为：

其中，q为动压头，其计算方式为：

q＝ρv²/2 (3-2)

ρ代表空气密度，v代表目标飞行器的飞行速度，c_x是没有量纲的阻力系数，c_y是没有量纲的升力系数，c_z为无量纲的侧向力系数，S为导弹的特征面积，对于有翼式导弹来说，特征面积为弹翼面积，对于无翼式导弹，特征面积为弹身最大横截面的面积。在方向的规定上，一般将速度的负方向定义为阻力的方向，升力的方向与O₃y₅轴的方向相同，侧向力的方向与O₃z₅轴的方向相同。

由空气动力提供的速度坐标系下的x轴、y轴、z轴的加速度为：

因为目标飞行器为轴对称式的，其飞行过程中的侧向力由侧滑角产生，且实际飞行过程中的滚转角较小，可忽略不计，因此可由弹道坐标系来代替速度坐标系。代入发射点惯性系与弹道系之间的转换矩阵C_1→4，将速度坐标系下的由空气动力引起的加速度转入发射点惯性系下，可得：

假设地球为均匀的球体，忽略地球形状摄动问题，地球对飞行器的引力只与其距地心的距离的平方成反比，根据万有引力定律可得：

其中，μ＝GM，G为引力常量，M为地球质量。联合公式(3-4)和公式(3-5)，根据飞行器的飞行动力学方程，可求得作用在目标飞行器上的发射点惯性系下的总的加速度为：

为了简化最终的机动目标模型，此处引入三个与气动参数有关的状态量Z_x、Z_y、Z_z：

将公式(3-7)代入公式(3-6)中，得：

结合公式(3-8)，根据飞行器的飞行动力学方程，可得：

选取滤波器的状态量为：

x＝[x v_x Z_x y v_y Z_y z v_z Z_z]^T (3-10)

滤波器的测量信息由地面的目标跟踪装置测得的位置信息提供，因此其测量方程为：

将公式(3-10)代入公式(3-9)中，得到滤波器的状态方程为：

3.2基于十二阶机动目标模型的跟踪滤波器

为了简化最终的机动目标模型，此处引入三个与气动参数有关的状态量Z_x、Z_y、Z_z：

则目标飞行器的加速度可表示为：

为了能够更加直接的表现出目标的加速度信息，可以再取三个状态量k_x、k_y、k_z分别作为Z_x、Z_y、Z_z这三个状态量的导数，结合飞行动力学公式，选取目标的位置信息和速度信息作为状态量，可得到机动目标模型为：

滤波器的测量信息由地面的目标跟踪装置测得的位置信息提供，因此其测量方程为：

选取滤波器的状态量为：

x＝[x v_x Z_x k_x y v_y Z_y k_y z v_z Z_z k_z]^T (3-17)

将公式(3-17)代入公式(3-15)中，得到滤波器的状态方程为：

机动目标跟踪滤波器的模型是非线性的，因此可以采用扩展卡尔曼滤波算法来对其进行加速度估计。在设计扩展卡尔曼滤波器前，应先将公式(3-18)的连续状态通过欧拉法转变为离散的状态方程。欧拉法的原理即为近似迭代，其公式为：

将公式(3-19)代入公式(3-18)中，得到滤波器的离散化状态方程与测量方程为：

结合公式(3-20)和公式(3-21)可得到滤波器的模型：

其中W代表滤波器的过程噪声，其方差为Q，方差越大，代表加速度受到的扰动影响越大，V为地面基站测量时带来的系统性偏差与测量噪声的和。分别对九阶机动目标模型和十二阶机动目标模型采用扩展卡尔曼滤波算法设计跟踪滤波器，其步骤如下所示：

第一步，对非线性模型进行线性化处理，得到雅可比矩阵φ_k：

H_k本身为线性的，不需要再进行泰勒展开。

第二步，预测状态先验估计值x_k+1|k：

x_k+1|k＝f(x_k|k) (3-24)

第三步，预测协方差矩阵的先验估计值P_k+1|k：

P_k+1|k＝φ_kP_k|kφ_k ^T+Q (3-25)

第四步，根据得到的先验估计值计算状态增益矩阵K_k+1：

K_k+1＝P_k+1|kH_k+1 ^T(H_k+1P_k+1|kH_k+1 ^T+R)^-1 (3-26)

第五步，根据状态增益矩阵修正得到状态后验估计值x_k+1|k+1：

x_k+1|k+1＝x_k+1|k+K_k+1(z_k+1-Hx_k+1|k) (3-27)

第六步，根据状态增益矩阵修正得到协方差后验估计值P_k+1|k+1：

P_k+1|k+1＝(I-K_k+1H)P_k+1|k (3-28)

在实际中，将非线性系统近似线性化处理这种方法很常见，因此扩展卡尔曼滤波算法的使用也更为常见。

把目标加速度a_tx、a_ty和a_tz的估计值从发射点惯性系投影到视线坐标系后，即得到a_tε和a_tβ。而目标位置x_t、y_t和z_t，以及目标速度V_tx、V_ty和V_tz的估计值则代入到下列方程中：Δx＝x-x_t，Δy＝y-y_t，Δz＝z-z_t，ΔV_x＝V_x-V_tx，ΔV_y＝V_y-V_ty，ΔV_z＝V_z-V_tz，然后就可以计算出目标—导弹相对距离和相对速度/>

4、视线俯仰角和视线偏航角的提取方法

4.1制导滤波器模型

导引头输出的只是目标—导弹视线相对导引头测量坐标系的视线角Δq_ε和Δq_β。导引头测量坐标系的圆点在导引头回转中心o′，在不考虑导引头安装误差及导引头安装位置挠性变形的情况下，其三个轴分别平行于弹体坐标系的ox₁轴、oy₁轴和oz₁轴。根据发射点惯性系与视线坐标系之间的转换关系来建立目标-导弹相对运动方程。

目标与导弹之间的相对距离矢量通过目标的距离矢量减去导弹的距离矢量来得到：

其中为发射点惯性系下目标与原点之间的距离矢量，/>为发射点惯性系下导弹与原点之间的距离矢量。对矢量求导可以分为绝对导数和相对导数，它们之间有如下关系：将矢量在发射点惯性系中对时间的导数减去同一矢量在动坐标系中对时间的导数求其差，等于该矢量与动坐标系的转动角速度的乘积。将公式(4-1)代入上述定理中，可得到目标-导弹相对距离矢量的绝对倒数与相对导数之间的关系：

其中，为相对距离矢量在发射点惯性坐标系下对时间的导数，/>即为视线坐标系下相对距离矢量对时间的导数，/>即为视线坐标系相对于发射点惯性系的旋转角速度，将公式(4-2)用矩阵的形式表达，可得：

其中，和/>分别为目标和导弹在视线系下的速度矢量，且/>

将公式(4-4)代入公式(4-3)可得：

令

对公式(4-6)两边相对时间求导，可得：

根据矢量的绝对导数与相对导数之间的关系，代入矢量得：

其中可通过结合公式(4-5)、公式(4-6)求得：

结合公式(4-4)、公式(4-8)、公式(4-9)，得：

将公式(4-3)代入公式(4-10)中，得到目标-导弹相对运动的动力学方程为：

其中a_r、a_ε、a_β分别是拦截导弹在视线坐标系下三个轴x、y、z方向上的加速度信息，a_tr、a_tε、a_tβ分别是机动目标在视线坐标系下三个轴x、y、z方向上的加速度信息。

考虑到弹体冲击振动对视线角速率的测量带来的影响，需要先建立弹体冲击振动模型。当纵向通道的开关脉冲发动机开关机时会引起弹体俯仰方向上的弹性变形，当侧向通道的开关脉冲发动机开关机时会引起弹体偏航方向上的弹性变形，因此在分析时，可以将冲击振动假设为单向通道来处理。只考虑单向通道的情况下，根据实际实验过程中产生的振动波形，采用模态分析的方法，来得到弹体的冲击振动模型。根据实际的波形的振动频率、振动幅度，采用连续时间下的传递函数为：

考虑到弹体的冲击振动可以假设为单向通道的来处理，因此可将冲击振动模型离散化，当作导引头测量过程中的扰动来处理。根据导引头的采样时间，选取Δt为0.0005s。一般传递函数离散化的方法有输入在采样时间间隔中保持不变的零阶保持法、输入在采样时间间隔内分段线性的一阶保持器法、脉冲响应不变离散化法和双线性变化法，此处采用的是双线性变换法进行离散化，双线性变换法的原理即为用两点之间的梯形面积来代替这两点之间的数值积分值，在使用过程中，只需要将连续传递函数中的s替换成公式(4-13)所示形式即可：

将公式(4-13)代入到公式(4-12)中，得

将其转为差分方程的形式，可得：

e_k＝a₁e_k-1-a₀e_k-2+b*(u_k-u_k-2) (4-15)

式中的u_k代表轨控发动机在导引头处的弹体坐标系下产生的冲击，e_k即为由弹体冲击引起的在导引头处的挠性变形信息。再根据公式(4-11)中得到的目标-导弹相对运动的动力学方程，选取状态量为：

将公式(4-16)代入公式(4-11)中，得到制导滤波器的状态方程为：

其中a_tε和a_tβ为视线坐标系下目标的加速度信息，a_ε和a_β为视线系下拦截导弹的加速度信息，选取导引头测得的视线俯仰角、视线偏航角作为测量值，得到制导滤波器的测量方程为：

再对上述状态方程和测量方程设计滤波器之前，需要先将其离散化，状态空间方程离散化一般有两个常用的方法，分别是欧拉法和零阶保持法，由于此模型比较复杂，因此此处采用的是欧拉法，其原理类似于近似微分，将欧拉公式应用于公式(4-17)中，得：

再将测量方程离散化，并考虑到模型误差，得到：

因为导引头测量中，弹体的冲击振动会影响测量结果，因此可以将公式(4-15)中得到的弹体的冲击振动引起的在导引头处的挠性变形作为测量过程中的扰动信息，添加到公式(4-20)中，因此制导滤波器的测量方程为：

z(k)＝y(k)+e(k)+v(k) (4-21)

结合公式(4-20)和公式(4-21)可得滤波器的模型为：

其中W代表滤波器的过程噪声，其方差为Q，方差越大，代表受到的扰动影响越大，V为量测的系统性偏差与测量噪声的和。

4.2扩展卡尔曼滤波器

采用扩展卡尔曼滤波算法这一非线性滤波算法来解决这一问题，有关扩展卡尔曼滤波算法的原理在前面已经做了详细的介绍，此处不再赘述，滤波器的设计步骤如下：

第一步，设定初始状态X₀、P₀；

第二步，利用公式(4-22)中的状态方程与测量方程做第一步预测：

x(k+1|k)＝f(x(k|k)) (4-23)

z(k+1|k)＝Hx(k+1|k) (4-24)

第三步，利用泰勒展开的方式求雅可比矩阵φ_k：

其中，

R₂₁＝(x₄(k)²(sin²(x₁(k))-cos²(x₁(k))))*Δt (4-26)

R₂₄＝-2x₄(k)sin(x₁(k))cos(x₁(k))*Δt (4-28)

R₄₂＝2x₄(k)tan(x₁(k))*Δt (4-30)

第四步，得到协方差矩阵的预测值：

P(k+1|k)＝φ_k+1|kP(k|k)φ_k+1|k ^T+Q (4-32)

第五步，计算卡尔曼增益矩阵：

K_k+1＝P(k+1|k)H^T(HP(k+1|k)H^T+R)^-1 (4-33)

第六步，更新状态与协方差矩阵：

x(k+1|k+1)＝x(k+1|k)+K_k+1(z(k+1|k+1)-z(k+1|k)) (4-34)

P(k+1|k+1)＝(I-K_k+1HP(k+1|k)) (4-35)

4.3粒子滤波器

粒子滤波算法最初是由贝叶斯滤波演化而来的，结合滤波器模型，粒子滤波算法的步骤为：

第一步，初始化；

第二步，循环

FOR i＝1:N

(1)通过先验概率p(x_k|x_k-1)采样，而p(x_k|x_k-1)可由状态方程得到；

(2)选取重要性概率密度函数，并计算权重，权重的计算方式为：

END FOR

第三步，权重归一化处理，得到

第四步，重采样，粒子集变为/>

第五步，计算k时刻的状态估计值

4.4扩展卡尔曼粒子滤波器

结合扩展卡尔曼滤波算法与粒子滤波算法的优点，可以设计出扩展卡尔曼粒子滤波(EPF)算法，此算法的原理即通过扩展卡尔曼滤波算法根据采样粒子进行状态预测与协方差预测，得到加权粒子与方差矩阵，再依据此加权粒子与方差矩阵计算权值并进行重采样，滤波器的性能得到了进一步的提高。EPF滤波算法的滤波流程如下。

k＝0时：初始化粒子，通过先验概率p(x_k|x_k-1)生成采样粒子p(x_k|x_k-1)依据状态方程得到。

k＝1,2,…：(1)加权粒子的预测，计算此非线性系统的雅可比矩阵，并将其分别记为

采用扩展卡尔曼滤波算法进行重采样得到加权粒子与方差矩阵：

(2)计算粒子的权重，并进行归一化处理；

(3)重采样处理得到处理后的粒子，重采样后的粒子权重均变为以此得到新的粒子集/>

(4)计算状态估计值以及每个粒子对应的方差阵。

4.5无迹卡尔曼粒子滤波器

采用无迹卡尔曼粒子滤波算法来设计制导滤波器，其计算步骤如下：

(1)初始化k＝0；当i＝1,2,…,N，通过计算先验概率p(x_k|x_k-1)生成采样粒子；

(2)重要性概率密度采样

将每个采样点通过无迹卡尔曼滤波算法得到新的加权粒子与加权矩阵/>并据此算法设定分布函数/>因此可以得到重要性概率密度函数为：

(3)计算权重，并做归一化处理；

(4)重采样处理，得到处理之后的粒子集

实验部分

(1)目标加速度估计

仿真基于X-51A飞行器对临近空间机动目标进行拦截的末制导过程，在末制导开始之前，设定初始时刻目标在发射点惯性系下的位置信息与速度信息为：

发射点惯性坐标系下目标初始位置：

[x_t0 y_t0 z_t0]＝[204.44 26.688 -25.362]km

目标初始速度：

[v_tx0 v_ty0 v_tz0]＝[-1513.749 46.097 -561.077]m/s

基于十二阶机动目标模型的跟踪滤波器的状态估计值的初始值设为：

P₀＝0.01*eye(12,12)

根据以上初始条件进行仿真，得到在发射点惯性系坐标下x轴、y轴、z轴方向上的加速度真实值与估计值的对比如图4a至图4c所示。

为了能更加清楚的对比九阶机动目标模型与十二阶目标机动目标模型，图5a至图5c展示了针对这两种模型的滤波结果的详细对比。

通过图5a至图5c可以看出，本发明中设计的改进的基于新的十二阶机动目标模型的滤波器相比于九阶机动目标模型而言具有更加良好的跟踪效果，其滤波精度也有了显著的提高，以上的误差对比图表明了此十二阶的模型在应用于针对临近空间机动目标的跟踪方面具有一定的理论工程意义。

(2)视线角速率提取

制导滤波器的仿真是基于采用动能拦截器对临近空间机动目标进行拦截的场景。导弹初始时刻的位置坐标和速度坐标为：

[x_m y_m z_m]＝[134.523 23.852 -25.464]km

[v_mx v_my v_mz]＝[1872.897 158.960 -511.337]m/s

目标初始时刻的位置坐标和速度坐标为：

[x_t y_t z_t]＝[204.44 26.688 -25.362]km

[v_tx v_ty v_tz]＝[-1513.749 46.097 -561.077]m/s

视线俯仰角速率、视线偏航角速率的真实值、测量值、估计值对比如图6和图7所示，视线俯仰角速率、视线偏航角速率的滤波估计值与真实值之间的误差如图8和图9所示。

为了能够更加清晰的对比这四种滤波算法的精度，将四种制导滤波器的视线俯仰角速率的滤波估计误差与视线偏航角的滤波估计误差分别放到同一张图中来进行对比分析，如图10和图11所示。

计算得到的各个滤波器的均方误差整理，如表1所示。

表1均方误差对比

从图10、图11和表1可看出，无迹卡尔曼粒子滤波算法设计的制导滤波器能更加贴近于真实值，从而在降低脱靶量方面有很大的帮助。

(3)末制导拦截仿真

搭建初始时刻的仿真环境设置如下：

发射点惯性坐标系下目标初始位置：[x_t0 y_t0 z_t0]＝[204.44 26.688 -25.362]km

目标初始速度：[v_tx0 v_ty0 v_tz0]＝[-1513.749 46.097 -561.077]m/s

导弹初始位置：[x_m0 y_m0 z_m0]＝[134.523 23.852 -25.464]km/s

导弹初始速度：[v_mx0 v_my0 v_mz0]＝[1872.897 158.960 -511.337]m/s

导弹初始姿态角：

仿真最终脱靶量为0.4116m，末制导环节时，目标与导弹的相对飞行过程如图12所示，目标与导弹的高度变化信息如图13所示。

目标在飞行过程中的速度、弹道倾角、弹道偏角如图14a、图14b、图14c所示，导弹在飞行过程中的速度、弹道倾角、弹道偏角如图15a、图15b、图15c所示。

导弹飞行过程中，末制导环节的x轴、y轴、z轴方向上的角速度如图16所示。

导弹飞行过程中末制导环节的俯仰角变化如图17a所示，偏航角变化如图17b所示，滚转角变化如图17c所示。

通过上述仿真实验表明，在给定的初始条件下，所设计的目标跟踪滤波器与制导滤波器都能实现导弹对临近空间机动目标的拦截，为了能更加充分的表现出目标跟踪滤波器和制导滤波器这两个滤波器对导弹末制导拦截的精度带来的影响，采用蒙特卡洛仿真的方法进行一百次随机仿真实验，末制导拦截的脱靶量变化情况如图18所示，燃料消耗情况如图19所示。

根据图18计算脱靶量在0.6m以下的概率为94％，脱靶量在1m以下的概率为99％，且经过100次蒙特卡洛仿真之后的平均脱靶量为0.3720m。导弹在末制导初始时刻所携带燃料为28kg，根据图19可知，末制导过程中燃料始终保持充足，证明了上述机动目标模型与制导滤波器的有效性与实用性。

Claims (4)

1.反大气层内或临近空间机动目标视线角速率估计方法，其特征在于，所述方法具体包括以下步骤：

步骤一、测算导弹加速度信息；

步骤二、基于十二阶机动目标模型建立跟踪滤波器，利用建立的跟踪滤波器跟踪目标加速度信息；

所述步骤二的具体过程为：

建立速度坐标系O₃x₅y₅z₅：将原点O₃定义在目标飞行器的质心上，O₃x₅轴沿着目标飞行器速度方向，O₃y₅轴定义在包含目标飞行器纵轴的纵向对称平面内，并与O₃x₅轴垂直，O₃z₅轴由右手定则确定；

建立弹道坐标系O₃x₆y₆z₆：O₃x₆轴沿着目标飞行器速度方向，O₃y₆轴定义在包含速度矢量的铅垂面内，并与O₃x₆轴垂直，O₃z₆轴由右手定则确定；

将作用在目标飞行器上的总的空气动力沿着速度坐标系分为阻力、升力和侧向力：

其中，X为阻力，Y为升力，Z为侧向力，q为动压头，c_x是无量纲的阻力系数，c_y是无量纲的升力系数，c_z为无量纲的侧向力系数，s为目标飞行器的特征面积；

q＝ρv²/2 (2)

其中，ρ代表空气密度，v代表目标飞行器的飞行速度；

由空气动力提供的速度坐标系下的O₃x₅轴、O₃y₅轴、O₃z₅轴的加速度为：

其中，a_X为空气动力提供的速度坐标系下的O₃x₅轴方向加速度，a_Y为空气动力提供的速度坐标系下的O₃y₅轴方向加速度，a_Z为空气动力提供的速度坐标系下的O₃z₅轴方向加速度，m是目标飞行器的质量，v_x是目标飞行器在速度坐标系O₃x₅轴方向的飞行速度，v_y是目标飞行器在速度坐标系O₃y₅轴方向的飞行速度，v_z是目标飞行器在速度坐标系O₃z₅轴方向的飞行速度；

利用弹道坐标系来代替速度坐标系，将速度坐标系下的由空气动力引起的加速度转入发射点惯性坐标系下：

其中，a_tx、a_ty、a_tz为发射点惯性坐标系下的加速度，C_1→4是发射点惯性坐标系与弹道坐标系之间的转换矩阵，θ代表弹道倾角，ψ_v代表弹道偏角；

根据万有引力定律得：

其中，g₀代表地球对目标飞行器的引力，r为目标飞行器与地心的距离，μ＝GM，G为引力常量，M为地球质量，R_e为地球半径，x、y和z是目标飞行器在发射点惯性坐标系下的坐标；

联合公式(4)和公式(5)，求得作用在目标飞行器上的发射点惯性坐标系下的总的加速度为：

其中，a_x、a_y、a_z分别代表作用在目标飞行器上的发射点惯性坐标系下x轴、y轴、z轴方向的总的加速度；

引入状态量Z_x、Z_y、Z_z：

则将目标飞行器的总加速度表示为：

再取状态量k_x、k_y、k_z分别作为状态量Z_x、Z_y、Z_z的导数，选取目标飞行器的位置信息和速度信息作为状态量，则得到的十二阶机动目标模型为：

其中，为x的一阶导数，/>为v_x的一阶导数，/>为Z_x的一阶导数，/>为k_x的一阶导数，/>为y的一阶导数，/>为v_y的一阶导数，/>为Z_y的一阶导数，/>为k_y的一阶导数，/>为z的一阶导数，/>为v_z的一阶导数，/>为Z_z的一阶导数，/>为k_z的一阶导数，w_n是零均值的高斯白噪声；

跟踪滤波器的测量信息由基站测得的位置信息提供，则跟踪滤波器的测量方程为：

其中，z₁、z₂、z₃是基站测得的目标飞行器在发射点惯性系下的三个轴上的位置信息；

选取跟踪滤波器的状态向量X为：

X＝[x v_x Z_x k_x y v_y Z_y k_y z v_z Z_z k_z]^T (11)

将公式(11)代入公式(9)中，得到跟踪滤波器的状态方程为：

其中，x′₁＝x，x′₂＝v_x，x′₃＝Z_x，x′₄＝k_x，x′₅＝y，x′₆＝v_y，x′₇＝Z_y，x′₈＝k_y，x′₉＝z，x′₁₀＝v_z，x′₁₁＝Z_z，x′₁₂＝k_z；

将公式(12)的连续状态方程转变为离散的状态方程，再对十二阶机动目标模型采用扩展卡尔曼滤波算法设计跟踪滤波器；

根据设计的跟踪滤波器对目标加速度a_tx、a_ty和a_tz进行估计，将目标加速度a_tx、a_ty和a_tz的估计值从发射点惯性坐标系投影到视线坐标系后，得到目标加速度a_tε和a_tβ；

步骤三、基于导弹加速度信息和目标加速度信息，并考虑弹体冲击振动的影响，对目标与导弹视线角速率进行估计；

所述步骤三的具体过程为：

目标与导弹之间的相对距离矢量利用目标的距离矢量减去导弹的距离矢量得到：

其中，为目标与导弹之间的相对距离矢量，/>为发射点惯性坐标系下目标与原点之间的距离矢量，/>为发射点惯性坐标系下导弹与原点之间的距离矢量；

则目标与导弹之间的相对距离矢量的绝对倒数与相对导数之间的关系为：

其中，为目标与导弹之间的相对距离矢量在发射点惯性坐标系下对时间的绝对导数，/>为目标与导弹之间的相对距离矢量在视线坐标系下对时间的相对导数，/>为视线坐标系相对于发射点惯性坐标系的旋转角速度；

将公式(14)表达为公式(15)的矩阵形式：

其中，r表示目标与导弹距离，表示目标与导弹距离的导数；

其中，和/>分别为目标和导弹在视线坐标系下的速度矢量，且

其中，q_ε为视线俯仰角，q_β为视线偏航角，为q_ε的一阶导数，/>为q_β的一阶导数；

将公式(16)代入公式(15)得：

令矢量为：

对公式(18)两边相对时间求导得：

其中，a_t和a分别为目标和导弹在视线坐标系下的加速度；

根据矢量的绝对导数与相对导数之间的关系，将矢量代入公式(19)得：

其中，通过结合公式(17)和公式(18)求得：

其中，为/>一阶导数，/>为/>一阶导数，/>为/>一阶导数；

结合公式(16)、公式(20)和公式(21)得：

将公式(15)代入公式(22)中，得到目标与导弹之间相对运动的动力学方程为：

其中，a_r、a_ε、a_β分别是拦截导弹在视线坐标系下o′x₄、o′y₄、o′z₄三个轴方向上的加速度信息，a_tr、a_tε、a_tβ分别是机动目标在视线坐标系下o′x₄、o′y₄、o′z₄三个轴方向上的加速度信息；

根据公式(23)的动力学方程选取状态向量X′为：

将公式(24)代入公式(23)中，得到制导滤波器的状态方程为：

其中，x″₁＝q_ε，x″₃＝q_β，/>

选取导引头测得的视线俯仰角和视线偏航角作为测量值，得到制导滤波器的测量方程为：

其中，y₁为导引头测得的视线俯仰角的测量值，y₂为导引头测得的视线偏航角的测量值；

对公式(25)的状态方程进行离散化，得：

其中，T为采样间隔，x″₁(k)为x″₁在k时刻的值，x″₂(k)为x″₂在k时刻的值，x″₃(k)为x″₃在k时刻的值，x″₄(k)为x″₄在k时刻的值；

考虑模型误差和弹体冲击振动影响的制导滤波器的测量方程为：

z(k)＝y(k)+e(k)+v(k) (28)

其中，z(k)为考虑模型误差和弹体冲击振动影响的k时刻的量测值，y(k)为不考虑模型误差和弹体冲击振动影响的k时刻的量测值，v(k)为k时刻的模型误差，e(k)为k时刻由弹体冲击引起的在导引头处的挠性变形信息；

结合公式(27)和公式(28)得到制导滤波器的模型：

其中，W代表制导滤波器的过程噪声，V为量测的系统性偏差与测量噪声之和，e为由弹体冲击引起的在导引头处的挠性变形信息，f(X′)为状态转移矩阵，H为量测转移矩阵，Z为目标量测；

基于公式(29)的制导滤波器模型和无迹卡尔曼粒子滤波算法对目标与导弹视线俯仰角速率和视线偏航角速率进行估计。

2.根据权利要求1所述的反大气层内或临近空间机动目标视线角速率估计方法，其特征在于，所述步骤一的具体过程为：

建立视线坐标系o′x₄y₄z₄：原点o′位于导引头回转中心；o′x₄轴与目标-导弹视线一致，由导引头回转中心指向目标方向为正；o′y₄轴位于包含o′x₄轴的铅垂面内，与o′x₄轴垂直且以指向上方为正；o′z₄轴按右手定则确定；

将导弹上加速度计输出的加速度在视线坐标系的o′y₄轴方向上的投影记为a_ε，将导弹上加速度计输出的加速度/>在视线坐标系的o′z₄轴方向上的投影记为a_β，a_ε和a_β即为导弹法向加速度。

3.根据权利要求2所述的反大气层内或临近空间机动目标视线角速率估计方法，其特征在于，所述目标飞行器的特征面积s为：

对于有翼式目标飞行器，目标飞行器的特征面积s为弹翼面积；对于无翼式目标飞行器，目标飞行器的特征面积s为弹身最大横截面的面积。

4.根据权利要求3所述的反大气层内或临近空间机动目标视线角速率估计方法，其特征在于，所述将公式(12)的连续状态方程转变为离散的状态方程采用的是欧拉法。