CN111474960A - 基于控制量特征辅助的临空高超声速目标跟踪方法 - Google Patents

Abstract

本发明提供基于控制量特征辅助的临空高超声速目标跟踪方法，包括基于控制量建模的方法和基于姿态角估计控制量的两种方法；方法一包括(1)控制量建模；(2)控制量离散化；(3)控制量加入状态方程；(4)目标跟踪模型；(5)系统量测方程；(6)基于平方根容积卡尔曼的跟踪算法。方法二包括(1)在弹体坐标系中计算目标的攻角和倾侧角；(2)根据光学传感器获取的目标姿态角信息估计攻角和倾侧角；(3)对攻角和倾侧角的递推方程离散化；(4)在目标的状态变量中增加攻角和倾侧角；(5)在量测方程中增加目标的姿态角观测；(6)通过平方根容积卡尔曼滤波对包含目标位置、速度、气动加速度和控制量的状态进行递推滤波。本发明的方法通过对高超声速目标控制量的建模和基于姿态角的估计方法，可为机动目标临空高超声速飞行器的跟踪提供了一种新的解决思路，通过案例实验表明本文方法的可行性和有效性。

Description

基于控制量特征辅助的临空高超声速目标跟踪方法

技术领域

本发明涉及高超声速目标探测跟踪领域技术，具体涉及一种基于控制量特 征辅助的临空高超声速目标跟踪方法。

背景技术

临近空间高超声速目标(Near Space Hypersonic Target，NSHT)，具有飞行跨 域大、速度快、气动力参数变化复杂等特点，给空天防御系统带来严峻挑战。 这类目标在飞行过程中具有纵向大空域非线性不规则机动及横向机动的能力， 且由于其高超声速飞行，RCS弱小，急剧增大了跟踪的不稳定性。实现对这类 目标的跟踪是对其拦截的基础。目前主要是从目标运动模型、滤波跟踪算法等 方面开展了对NSHT的跟踪研究。根据Li X.R.的《Survey of maneuvring target tracking.Part II：Motion models of ballistic andspace targets》(IEEE Transactions on Aerospace and Electronic Systems，2010，46(1)：96-119.)研究结论，针对空间目标， 使用简单的运动模型会损失很多重要信息，而动力学模型利用目标的运动特性 信息，跟踪效果较好。因此，为了实现对这类目标的跟踪，需要充分利用关于 目标运动特性的先验信息，利用尽可能多的特性信息。

李海宁等在《面向跟踪的吸气式高超声速飞行器动力学建模》(航空学报 2014，35(6)：1651-1664.)中提出，根据吸气式高超声速目标的动力学特性，在重 力转弯模型框架中，面向目标跟踪的动力学混合模型集在跟踪NSHT时具有一 定的优势，但作者给出的非解析解结构难以应用于线性化方法的滤波器设计。 张凯等在《一种基于气动力模型的高超声速滑翔目标跟踪算法》(宇航学报， 2017，38(2)：123-130.)中针对高超声速滑翔目标的跟踪问题，提出了一种基于气 动力模型的目标跟踪算法，利用转弯力参数与爬升力参数之间的先验知识，对 传统气动力模型进行了合理的改进，但仿真中没有考虑控制变量变化时的适用 性。翟岱亮等在《基于自适应IMM的高超声速飞行器轨迹预测》(航空学 报，2016，37(11)3466-3475.)中，使用基于高超声速目标的气动力模型和自适应 IMM，对固定攻角固定而和倾侧角不同变化条件下的高超声速目标进行了跟踪， 但该文只是对气动参数进行建模，没有考虑到控制量因素。Fan Y.等在《A cost-effective tracking algorithm forhypersonic glide vehicle maneuver based on modified aerodynamic model》(Applied Sciences，2016，6(312)：1-12)中用匀速高斯马 尔科夫过程对攻角和侧滑角进行建模，建立了修改的气动力模型，仿真表明了 算法的有效性，但该建模方法难以进一步提高跟踪精度。

目前，高超声速目标的跟踪方法在跟踪模型方面主要从运动特性、统计特 性和动力学模型建模三个角度进行研究，在跟踪算法方面主要从改进非线性跟 踪算法角度进行研究，这些方法虽然在提高跟踪精度方面有一定的提高，但由 于该方法缺乏对高超声速目标特征信息的深入分析和利用，很难进一步提高跟 踪精度。因此，迫切需要针对高超声速目标的特点，深入分析其机动特性的本 质，提出有创新性的跟踪方法，较好地解决该类目标的跟踪问题。

根据高超声速目标的运动特性分析可知，在滑翔段，高超声速目标的攻角 和倾侧角分别决定目标的纵向运动特性和横向机动特性。攻角作为辅助控制参 数，通常是固定的，不能频繁跳变，一般在标称攻角基础上进行小范围微调。 同时，作为主控制量的倾侧角也只能在小范围内调整。随着探测跟踪传感器性 能的提升，探测跟踪系统不但可以输出目标的运动信息，也可以提取目标的姿 态角等特征信息。与运动目标其它特征信息相比，姿态角信息更能反映机动目 标的机动特性和运动方式的不确定性。因此，在高超声速目标属性约束和传感 器性能约束的基础上，针对高超声速目标的特性进行深入分析，可以考虑在传 统目标跟踪的基础上，采用控制量建模的方法和基于姿态角估计的方法，进一 步提高高超声速目标的跟踪精度，研究思路示意图如图1所示。

发明内容

针对临近空间高超声速目标运动特性复杂难以跟踪的问题，考虑采用在机 动目标跟踪方面具有潜力的特征辅助跟踪技术，以期进行一步提高NSHT的跟 踪精度。具体来说，就是分别采用控制量建模以及基于姿态角估计控制量的方 法，建立跟踪模型，并采用动力学模型完成对NSHT的跟踪。

首先，对目标控制量特征与临近空间高超声速目标机动特性的关系进行分 析，然后在动力学模型的基础上，给出控制量特征的建模方法，提出基于控制 量建模的跟踪方法；然后在光学传感器获取到的目标姿态信息的基础上，提出 利用姿态角估计控制量的跟踪方法；最后，基于目标跟踪的动力学模型，提出 基于高超声速目标两种特征辅助跟踪方法的高超声速目标跟踪方法。

如上所述，本发明提供一种基于控制量特征辅助的临空高超声速目标跟踪 方法，采用对控制量特征直接建模，具体包括下列步骤：

Step1.控制量建模

把攻角和倾侧角这两个控制变量看作匀速高斯马尔科夫过程，用加性白噪 声来对攻角和倾侧角的变化率进行建模，表达式如下：

式中，α表示攻角，v表示倾侧角，t表示时间变量，ω_t，α和ω_t，v是未知控制输入 引起的在攻角和倾侧角上一定程度的白噪声；

Step2.控制量离散化

式(1)可以进一步离散化为：

式中，T是采样周期的时间间隔，α_k+1、v_k+1分别表示攻角和倾侧角的一步预测 值；认为攻角和倾侧角在离散时刻k+1和k的变化分别是时间间隔T乘以白噪声 ω_α，k和ω_v，k；

Step3.控制量加入状态模型

在完成姿态角建模后，结合机动目标跟踪中状态方程扩维的思路，将攻角 和倾侧角加入状态方程中进行更新，对目标状态向量进行扩维，则目标的状态 变量变为11维的向量：

x(k)＝[r(k)，v(k)，a_VTC(k)，α(k)，v(k)]^T (3)

式中，k表示时间序列变量；r(k)＝[x(k)，y(k)，z(k)]^T和

分别为目标在东北天坐标系ENU中第k时刻的位置和速度向量，a_VTC(k)＝[α_V(k)，α_T(k)，α_C(k)]^T为目标在半速度VTC坐标系下的加速度，α(k)和 v(k)分别是攻角和倾侧角；

Step4.目标跟踪模型

建立用于目标跟踪的离散模型为：

x_k+1＝Φ_kx_k+B_ku_k+Γ_kω_k (4)

式中，x_k为式(3)中x(k)的简化表示，x_k、x_k+1分别为第k+1和k时刻目标的 状态向量，

为控制输入向量，g_xk，g_yk，g_zk为重力 加速度在ENU坐标系三个方向上的分量，能够通过平面、球形或椭球三种地球 模型给出，

分别为机动加速度在VTC坐标系下加速度的均值， Φ_k＝[I_3×3 TI_3×3P₀vT₁ 0_3×2；0_3×3 I_3×3 P₀vT₂ 0_3×2；0_3×3 0_3×3 T₃ 0_3×2；0_3×3 0_3×3 0_3×3 T₄]^T为状态转移 矩阵，

为输入矩阵， Γ_k＝[P₀vΓ₁；P₀vΓ₂；Γ₃；Γ₄]^T为噪声输入矩阵，ω_k为协方差矩阵为Q_k的VTC坐标 系下的加速度过程噪声，状态转移矩阵Φ_k和噪声输入矩阵Γ_k按式(2)进行修正， 能够建立控制量姿态角辅助跟踪的离散状态方程，如式(4)所示；其中，

I_3×3＝blkdiag(1，1，1)，0_3×3＝blkdiag(0，0，0)，0_3×2＝[0 0；0 0；0 0]，

T₁＝[T₁₁η₁₁ T₁₂η₂₁ T₁₃η₃₁；T₂₁η₁₁T₂₂η₂₁ T₂₃η₃₁；T₃₁η₁₁ T₃₂η₂₁ T₃₃η₃₁]，

T₂＝[T₁₁η₁₂ T₁₂η₂₂ T₁₃η₃₂；T₂₁η₁₂ T₂₂η₂₂ T₂₃η₃₂；T₃₁η₁₂ T₃₂η₂₂ T₃₃η₃₂]

T₃＝[η₁₃ 0 0；0 η₂₃ 0；0 0 η₃₃]，T₄＝[1 0；0 1]，

B₁＝[υ_VT₁₁ξ₁₁ υ_TT₁₂ξ₂₁ υ_CT₁₃ξ₃₁；υ_VT₂₁ξ₁₁ υ_TT₂₂ξ₂₁ υ_CT₂₃ξ₃₁；υ_VT₃₁ξ₁₁ υ_TT₃₂ξ₂₁ υ_CT₃₃ξ₃₁]，

B₂＝[υ_VT₁₁ξ₁₂ υ_TT₁₂ξ₂₂ υ_CT₁₃ξ₃₂；υ_VT₂₁ξ₁₂ υ_TT₂₂ξ₂₂ υ_CT₂₃ξ₃₂；υ_VT₃₁ξ₁₂ υ_TT₃₂ξ₂₂ υ_CT₃₃ξ₃₂]，

B₃＝[υ_Vξ₁₃ 0 0；0 υ_Tξ₂₃ 0；0 0 υ_Cξ₃₃]，

Γ₁＝[T₁₁ξ₁₁ T₁₂ξ₂₁ T₁₃ξ₃₁；T₂₁ξ₁₁ T₂₂ξ₂₁ T₂₃ξ₃₁；T₃₁ξ₁₁ T₃₂ξ₂₁ T₃₃ξ₃₁]，

Γ₂＝[T₁₁ξ₁₂ T₁₂ξ₂₂ T₁₃ξ₃₂；T₂₁ξ₁₂ T₂₂ξ₂₂ T₂₃ξ₃₂；T₃₁ξ₁₂ T₃₂ξ₂₂ T₃₃ξ₃₂]，

Γ₃＝[ξ₁₃ 0 0；0 ξ₂₃ 0；0 0 ξ₃₃]，Γ₄＝[T 0；0 T]，P₀为自由流压强，

为目标平均速度，

为目标在水平面的 平均速度，

T₃₂＝0，

υ_V、υ_T、υ_C分别为目标在VTC坐标系三坐标轴向上的机动频率；

Step5.系统量测方程

假设系统探测中，只有雷达传感器，则系统的量测z_k为：

式中，R_k、A_k、E_k、σ分别表示测量的距离、方位角、俯仰角、量测噪声，、 x、y、z为ENU坐标系中的位置分量，h(x_k)为系统非线性量测方程，ΔR、ΔA 和ΔE分别为雷达的距离量测误差、方位量测误差和俯仰量测误差；

Step6.基于平方根容积卡尔曼滤波的跟踪算法

针对目标状态模型和观测模型的非线性特点，采用平方根容积卡尔曼滤波 SRCKF进行目标的状态滤波；SRCKF直接对协方差矩阵的平方根进行迭代计 算，避免了每步迭代的矩阵开方运算，能够克服CKF算法中滤波发散的问题， 从而提高滤波的精度；

合并状态方程和量测方程为：

式中，

f(·)表示非线性状态方程，W_k-1是均值为0，方差为Q_k-1的高斯过程噪声向量，h(X_k)是非线性的量测方程，量测噪声V_k是方差为R_k的 零均值高斯白噪声矢量。

本发明还提供另一种基于控制量特征辅助的临空高超声速目标跟踪方法， 采用基于姿态角的控制量特征估计，具体包括下列步骤：

利用光学传感器获取的目标姿态信息，对高超声速飞行器的攻角和倾侧角 进行估计，以期进行一步提高目标的跟踪精度；具体处理流程如下：

STEP1.在弹体坐标系中计算目标的攻角和倾侧角

在弹体坐标系与速度坐标系中计算攻角α和倾侧角v的方法如下：

式中，v_dx、v_dy、v_dz是速度矢量在弹体坐标系上的分量，

为目标的飞行速度；

STEP2.根据光学传感器获取的目标姿态角信息估计攻角和倾侧角

假设光学传感器可以测量地理坐标系与弹体坐标系的三个欧拉角，分别为 偏航角ψ、俯仰角ε和滚转角γ，则由地理坐标系到弹体坐标系的变换矩阵M_D为：

则地理坐标系的速度矢量与弹体坐标系的速度矢量的关系为：

式中，v_E、v_N、v_U分别表示目标在地理坐标系中速度分量；

由此，根据式(19)、(20)和式(21)建立由光学传感器观测的目标姿态角信息 估计的目标攻角和倾侧角；

STEP3.对攻角和倾侧角的递推方程离散化

由于NSHT的动力学模型为时域系统，对目标进行跟踪时，需要将其离散 化，离散化后，将控制变量的递推方程离散化为：

式中，α_k+1、v_k+1分别表示攻角和倾侧角的一步预测值，T为采样周期；认为攻角 和倾侧角在离散时刻k和k-1的变化分别是时间间隔T乘以零均值的白噪声ω_α，k和ω_v，k；

STEP4.在目标的状态变量中增加攻角和倾侧角

建立两个控制量的递推方程后，把通过姿态角估计的两个控制量加入目标 的状态方程，则状态方程修正为：

x(k)＝[r(k)，v(k)，a_VTC(k)，α(k)，v(k)]^T (11)

式中各变量的含义与前面相同，只是攻角α(k)和倾侧角v(k)是基于姿态角的估计值；

STEP5.系统的量测模型

系统量测模型中除了距离、方位和俯仰测量信息，还有基于光学传感器获 取的目标的姿态信息的量测模型为：

式中，在观测向量z_k中，前三项分别为距离R_k、方位A_k和俯仰角E_k，ψ_k、θ_k和 γ_k分别为光学传感器观测的偏航角、俯仰角和滚转角；各量测噪声均是已知的 独立白噪声，量测噪声方差可以表示为

式中，

和

分别是雷达量测的距离、方位、俯仰和光学传感器观测的航角、俯仰角和滚转角方差；

STEP6.基于平方根容积卡尔曼滤波的跟踪算法

根据上面的非线性量测方程和非线性状态方程，采用SRCKF方法对目标位 置、速度、气动加速度和控制量的状态进行递推滤波。

在本发明的一个实施例中，SRCKF算法步骤如下：

Step1初始化

设置状态初值x_0|0，协方差矩阵初值P_0|0，协方差矩阵平方根因子的初始值 S_0|0，其中

Step2时间更新：k＝1，2，…

(1)计算当前状态的2n个容积点X_i，k-1|k-1，i＝1，2，…，2n

其中，S_k-1|k-1为k-1时刻协方差矩阵P_k-1的平方根因子，

为k-1时刻的状态 估计值，

为基本容积点，m＝2n为容积点的总数，其中n为状态维数， [1]_i为点集[1]的第i个点；

(2)计算容积点经过非线性状态转移函数的预测值，i＝1，2，…，2n：

(3)结合权值和容积点预测值，SRCKF采用相等权值，预测状态可计算为 m＝2n：

(4)计算预测误差协方差矩阵的平方根因子：

S_k|k-1＝Tria([χ_k|k-1 S_Q，k-1]) (17)

其中，S_Q，k-1是过程噪声Q_k-1的平方根且有中间变量χ_k|k-1，

其中，算法S＝Tria(A)意为先对矩阵A进行QR分解，得到一个正规正交矩 阵B和一个上三角矩阵C，令S＝C^T，得到的S即为上三角矩阵；

Step3量测更新：k＝1，2，…

(1)计算更新的状态容积点，i＝1，2，…，m：

(2)计算容积点经量测非线性函数传播后的值：

Z_i，k|k-1＝h(X_i，k|k-1) (19)

(3)估计预测量测：

(4)估计新息的协方差矩阵：

S_zz，k|k-1＝Tria([Z_k|k-1 S_R，k]) (21)

其中，

(5)估计互协方差矩阵：

其中，

(6)估计SRCKF的滤波增益：

(7)基于k时刻新的量测z_k，更新系统状态

(8)更新误差协方差矩阵的平方根因子

S_k|k＝Tria([χ_k|k-1-W_kZ_k|k-1 W_kS_R，k]) (25)

本发明针对临空高超声速目标的跟踪问题，提出基于控制量特征进行辅助 跟踪的思路，并给出控制量特征辅助的跟踪方法，为临空高超声速目标的跟踪 提供了新的解决思路。

附图说明

图1基于控制量建模和姿态角估计的特征辅助跟踪思路示意图；

图2示出NSHT机动方式1条件下的跟踪结果，其中图2(a)示出NSHT运 动轨迹，图2(b)示出NSHT位置跟踪误差，图2(c)示出NSHT速度跟踪误差， 图2(d)示出NSHT加速度跟踪误差；

图3示出NSHT机动方式2条件下的跟踪结果，其中图3(a)示出NSHT运 动轨迹，图3(b)示出NSHT位置跟踪误差，图3(c)示出NSHT速度跟踪误差， 图2(d)示出NSHT加速度跟踪误差；

图4示出NSHT机动方式3条件下的跟踪结果，其中图4(a)示出目标运动 轨迹，图4(b)示出NSHT位置跟踪误差，图4(c)示出NSHT速度跟踪误差，图 4(d)示出NSHT加速度跟踪误差。

具体实施方式

下面结合附图详细介绍本发明。

1.对控制量特征直接建模

Step1.控制量建模

攻角和倾侧角作为作战对手的一个重要秘密，防御方甚至无法知道这两个 控制变量的变化规律。因此，需要对这两个控制变量进行有效的建模。NSHT 在临近空间以超过5马赫的速度机动，当攻角和倾侧角变化剧烈时，NSHT需要 承受结构上的破坏。在NSHT的实际控制中，为了确保稳定性，这两个控制变 量常常在某一度数范围内连续变化。通过以上分析可知，这两个控制变量的变 化率在短时间内变化不剧烈。

因此，可以把攻角和倾侧角这两个控制变量看作匀速高斯马尔科夫过程， 可以用加性白噪声来对攻角和倾侧角的变化率进行建模，表达式如下：

式中，α表示攻角，v表示倾侧角，t表示时间变量，ω_t，α和ω_t，v是未知控制输入 引起的在攻角和倾侧角上一定程度的白噪声。

Step2.控制量离散化

式(1)可以进一步离散化为：

式中，T是采样周期的时间间隔，α_k+1、v_k+1分别表示攻角和倾侧角的一步预测 值。可以认为攻角和倾侧角在离散时刻k+1和k的变化分别是时间间隔T乘以白 噪声ω_α，k和ω_v，k。

Step3.控制量加入状态模型

在完成姿态角建模后，结合机动目标跟踪中状态方程扩维的思路，将攻角 和倾侧角加入状态方程中进行更新，对目标状态向量进行扩维，则目标的状态 变量变为11维的向量：

x(k)＝[r(k)，v(k)，a_VTC(k)，α(k)，v(k)]^T (28)

式中，k表示时间序列变量；r(k)＝[x(k)，y(k)，z(k)]^T和

分别为目标在东北天坐标系(ENU)中第k时刻的位置和速度向量， a_VTC(k)＝[α_V(k)，α_T(k)，α_C(k)]^T为目标在半速度(Velocity-Turn-Climb，VTC)坐标 系下的加速度，α(k)和v(k)分别是攻角和倾侧角。

Step4.目标跟踪模型

建立用于目标跟踪的离散模型为：

x_k+1＝Φ_kx_k+B_ku_k+Γ_kω_k (29)

式中，x_k为式(3)中x(k)的简化表示，x_k、x_k+1分别为第k+1和k时刻目标的 状态向量，

为控制输入向量，g_xk，g_yk，g_zk为重力 加速度在ENU坐标系三个方向上的分量，可通过平面、球形或椭球三种地球模 型给出，

分别为机动加速度在VTC坐标系下加速度的均值， Φ_k＝[I_3×3 TI_3×3P₀vT₁ 0_3×2；0_3×3 I_3×3 P₀vT₂ 0_3×2；0_3×3 0_3×3 T₃ 0_3×2；0_3×3 0_3×3 0_3×3 T₄]^T为状态转移 矩阵，

为输入矩阵， Γ_k＝[P₀vΓ₁；P₀vΓ₂；Γ₃；Γ₄]^T为噪声输入矩阵，ω_k为协方差矩阵为Q_k的VTC坐标 系下的加速度过程噪声，状态转移矩阵Φ_k和噪声输入矩阵Γ_k按式(2)进行修正， 即可建立控制量姿态角辅助跟踪的离散状态方程，如式(4)所示。其中，

I_3×3＝blkdiag(1，1，1)，0_3×3＝blkdiag(0，0，0)，0_3×2＝[0 0；0 0；0 0]，

T₁＝[T₁₁η₁₁ T₁₂η₂₁ T₁₃η₃₁；T₂₁η₁₁ T₂₂η₂₁ T₂₃η₃₁；T₃₁η₁₁ T₃₂η₂₁ T₃₃η₃₁]，

T₂＝[T₁₁η₁₂ T₁₂η₂₂ T₁₃η₃₂；T₂₁η₁₂ T₂₂η₂₂ T₂₃η₃₂；T₃₁η₁₂ T₃₂η₂₂ T₃₃η₃₂]

T₃＝[η₁₃ 0 0；0 η₂₃ 0；0 0 η₂₃]，T₄＝[1 0；0 1]，

B₁＝[υ_VT₁₁ξ₁₁ υ_TT₁₂ξ₂₁ υ_CT₁₃ξ₃₁；υ_VT₂₁ξ₁₁ υ_TT₂₂ξ₂₁ υ_CT₂₃ξ₃₁；υ_VT₃₁ξ₁₁ υ_TT₂₂ξ₂₁ υ_CT₂₃ξ₃₁]，

B₂＝[υ_VT₁₁ξ₁₂ υ_TT₁₂ξ₂₂ υ_CT₁₃ξ₃₂；υ_VT₂₁ξ₁₂ υ_TT₂₂ξ₂₂ υ_CT₂₃ξ₃₂；υ_VT₃₁ξ₁₂ υ_TT₃₂ξ₂₂ υ_CT₃₃ξ₃₂]，

B₃＝[υ_Vξ₁₃ 0 0；0 υ_Tξ₂₃ 0；0 0 υ_Cξ₃₃]，

Γ₁＝[T₁₁ξ₁₁ T₁₂ξ₂₁ T₁₃ξ₃₁；T₂₁ξ₁₁ T₂₂ξ₂₁ T₂₃ξ₃₁；T₃₁ξ₁₁ T₃₂ξ₂₁ T₃₃ξ₃₁]，

Γ₂＝[T₁₁ξ₁₂ T₁₂ξ₂₂ T₁₃ξ₃₂；T₂₁ξ₁₂ T₂₂ξ₂₂ T₂₃ξ₃₂；T₃₁ξ₁₂ T₃₂ξ₂₂ T₃₃ξ₃₂]，

Γ₃＝[ξ₁₃ 0 0；0 ξ₂₃ 0；0 0 ξ₂₃]，Γ₄＝[T 0；0 T]，P₀为自由流压强，

为目标平均速度，

为目标在水平面的 平均速度，

T₃₂＝0，

υ_V、υ_T、υ_C分别为目标在VTC坐标系三坐标轴向上的机动频率。

Step5.系统量测方程

假设系统探测中，只有雷达传感器，则系统的量测z_k为：

式中，R_k、A_k、E_k、σ分别表示测量的距离、方位角、俯仰角、量测噪声，、 x、y、z为ENU坐标系中的位置分量，h(x_k)为系统非线性量测方程，ΔR、ΔA 和ΔE分别为雷达的距离量测误差、方位量测误差和俯仰量测误差，更详细的描 述见Li X.R.的《Survey ofmaneuvering target tracking-part III：measurement models》(Signal and DataProcessing of Small Targets，Bellingham，WA，USA：SPIE Press，2001：423-446)。

Step6.基于平方根容积卡尔曼滤波的跟踪算法

针对目标状态模型和观测模型的非线性特点，考虑采用相比其它非线性滤 波算法，计算量小、实时性好且滤波稳定性更好的平方根容积卡尔曼滤波 (Square-rootCubature Kalman Filter，SRCKF)进行目标的状态滤波。SRCKF直接对协 方差矩阵的平方根进行迭代计算，避免了每步迭代的矩阵开方运算，可克服CKF 算法中滤波发散的问题，从而提高滤波的精度。

合并状态方程和量测方程为：

式中，

f(·)表示非线性状态方程，W_k-1是均值为0，方差为Q_k-1的高斯过程噪声向量，h(X_k)是非线性的量测方程，量测噪声V_k是方差为R_k的 零均值高斯白噪声矢量。

SRCKF算法步骤如下：

Step1初始化

设置状态初值x_0|0，协方差矩阵初值P_0|0，协方差矩阵平方根因子的初始值 S_0|0，其中

Step2时间更新(k＝1，2，…)

(1)计算当前状态的2n个容积点X_i，k-1|k-1(i＝1，2，…，2n)

其中，S_k-1|k-1为k-1时刻协方差矩阵P_k-1的平方根因子，

为k-1时刻的状态 估计值，

为基本容积点，m＝2n为容积点的总数，其中n为状态维数， [1]_i为点集[1]的第i个点。

(2)计算容积点经过非线性状态转移函数的预测值(i＝1，2，…，2n)：

(3)结合权值和容积点预测值，SRCKF采用相等权值，预测状态可计算为 (m＝2n)：

(4)计算预测误差协方差矩阵的平方根因子：

S_k|k-1＝Tria([χ_k|k-1 S_Q，k-1]) (35)

其中，S_Q，k-1是过程噪声Q_k-1的平方根且有中间变量χ_k|k-1，

其中，算法S＝Tria(A)意为先对矩阵A进行QR分解，得到一个正规正交矩 阵B和一个上三角矩阵C，令S＝C^T，得到的S即为上三角矩阵。

Step3量测更新(k＝1，2，…)

(1)计算更新的状态容积点(i＝1，2，…，m)：

(2)计算容积点经量测非线性函数传播后的值：

Z_i，k|k-1＝h(X_i，k|k-1) (37)

(3)估计预测量测：

(4)估计新息的协方差矩阵：

S_zz，k|k-1＝Tria([Z_k|k-1 S_R，k]) (39)

其中，

(5)估计互协方差矩阵：

其中，

(6)估计SRCKF的滤波增益：

(7)基于k时刻新的量测z_k，更新系统状态

(8)更新误差协方差矩阵的平方根因子

S_k|k＝Tria([χ_k|k-1-W_kZ_k|k-1 W_kS_R，k]) (43)

2.基于姿态角的控制量特征估计

由于跟踪模型中很难对这类目标的气动力进行准确建模，导致跟踪精度难 以进一步提高。随着探测跟踪传感器性能的提升，探测跟踪系统不但可以输出 目标的运动信息，也可以提取目标的姿态角等特征信息。与运动目标其它特征 信息相比，姿态角信息更能反映机动目标的机动特性和运动方式的不确定性。 因此，使用目标姿态角进行辅助跟踪，可望进一步提高机动目标跟踪性能。在 对攻角和倾侧角进行建模的基础上，本发明利用光学传感器获取的目标姿态信 息，对高超声速飞行器的攻角和倾侧角进行估计，以期进行一步提高目标的跟 踪精度。具体处理流程如下：

Step1.在弹体坐标系中计算目标的攻角和倾侧角

在弹体坐标系与速度坐标系中计算攻角α和倾侧角v的方法如下：

式中，v_dx、v_dy、v_dz是速度矢量在弹体坐标系上的分量，

为目标的飞行速度。

Step2.根据光学传感器获取的目标姿态角信息估计攻角和倾侧角

假设光学传感器可以测量地理坐标系与弹体坐标系的三个欧拉角，分别为 偏航角ψ、俯仰角ε和滚转角γ，则由地理坐标系到弹体坐标系的变换矩阵M_D为：

则地理坐标系的速度矢量与弹体坐标系的速度矢量的关系为：

式中，v_E、v_N、v_U分别表示目标在地理坐标系中速度分量。

由此，根据式(19)、(20)和式(21)可以建立由光学传感器观测的目标姿态角 信息估计的目标攻角和倾侧角。

Step3.对攻角和倾侧角的递推方程离散化

由于NSHT的动力学模型为时域系统，对目标进行跟踪时，需要将其离散 化，离散化后，控制变量的递推方程可以离散化为：

式中，α_k+1、v_k+1分别表示攻角和倾侧角的一步预测值，T为采样周期。可以认为 攻角和倾侧角在离散时刻k和k-1的变化分别是时间间隔T乘以零均值的白噪 声ω_α，k和ω_v，k。

Step4.在目标的状态变量中增加攻角和倾侧角

建立两个控制量的递推方程后，可把通过姿态角估计的两个控制量加入目 标的状态方程，则状态方程可修正为：

x(k)＝[r(k)，v(k)，a_VTC(k)，α(k)，v(k)]^T (48)

式中各变量的含义与前面相同，只是攻角α(k)和倾侧角v(k)是基于姿态角的估计值。

Step5.系统的量测模型

系统量测模型中除了距离、方位和俯仰测量信息，还有基于光学传感器获 取的目标的姿态信息的量测模型为：

式中，在观测向量z_k中，前三项分别为距离R_k、方位A_k和俯仰角E_k，ψ_k、θ_k和 γ_k分别为光学传感器观测的偏航角、俯仰角和滚转角。各量测噪声均是已知的 独立白噪声，量测噪声方差可以表示为

式中，

和

分别是雷达量测的距离、方位、俯仰和光学传感器观测的航角、俯仰角和滚转角方差。

Step6.基于平方根容积卡尔曼滤波的跟踪算法

根据上面的非线性量测方程和非线性状态方程，同样采用SRCKF方法对目 标位置、速度、气动加速度和控制量的状态进行递推滤波。

上面基于姿态角的控制量估计方法，完成了对高机动目标关键控制参数的 估计，是对相应运动模式高超声速目标跟踪难题提供一种解决思路。

具体实施例

为验证本发明算法的可行性和有效性，以助推滑翔式高超声速飞行器作为 仿真分析对象，针对临空高超声速目标的机动特性，根据翟岱亮在《面向轨迹 预测的高超声速飞行器气动性能分析》(固体火箭技术，2017，40(1)：115-120)中 设定的仿真场景，设定目标仿真初始状态为 s₁＝[0km 0km 40km 2.4km/s 0km/s 0km/s]^T，雷达部署位置为[500000m 0m 0m]^T，雷 达的距离、俯仰和方位跟踪精度为[30m 0.05° 0.05°]^T。

为了充分验证本文算法的有效性，设计了高超声速目标的三种典型运动方 式：

(1)目标保持固定攻角为10°，固定倾侧角为20°飞行。

(2)目标攻角按幅度分别为6.5°和10.5°、周期为100s矩形波变化，倾 侧角始终保持0°。

(3)目标攻角固定为10°，倾侧角按幅度为20°、周期为100s的正弦波 变化。

根据前面的仿真场景，分别对三种典型机动模式下，使用本发明提出的控 制量建模方法(简称特征辅助方法1)和基于姿态角的控制量估计方法(简称特 征辅助方法2)对高超声速目标进行跟踪，并分别与翟岱亮在《面向轨迹预测的 高超声速飞行器气动性能分析》使用气动参数作为控制量的运动模型(简称为 原气动模型)的目标跟踪仿真结果进行比较。实验中，进行50次Monte Carlo 仿真实验，以验证本发明所提算法的有效性。

由图1可知，虽然在机动方式1的条件下，NSHT既存在纵向机动也存在横 向机动，但三种跟踪方法均可实现对机动目标的跟踪。通过位置、速度以及加 速度跟踪的误差曲线可知，对控制量建模方法和基于姿态角估计控制量的方法 取得了更好的跟踪效果。这充分说明了使用控制量特征进行机动目标辅助跟踪 的可行性。

由图2可知，目标在机动方式2的条件下，目标在纵向上的位置、速度和 加速度都发生了变化，但三种跟踪方法均可实现对目标的跟踪。通过位置、速 度以及加速度跟踪的误差曲线可知，对控制量建模和基于姿态角估计控制量的 方法取得了更好的跟踪效果。同时，两种特征辅助的跟踪方法可以取得更快的 收敛速度，由于特征辅助跟踪2可以精确估计控制量，取得了更好的跟踪效果。

由图3可知，目标在机动方式3下的跟踪效果几乎可以取得和前面两种机 动方式一致的跟踪结论，充分表明特征辅助跟踪方法在各种机动方式下的鲁棒 性。

通过分析图1～图3可知：基于特征辅助的两种跟踪方法都比原来的气动模 型取得更好的跟踪效果，由于控制量特征与NSHT机动的约束关系密切，通过 对控制量进行建模可以获取更好的跟踪效果。基于姿态角估计攻角和倾侧角控 制量的特征辅助跟踪方法的收敛速度和跟踪精度要优于对控制量特征直接建模 的跟踪方法，这是因为基于姿态角估计控制量的方法可以比直接对控制量进行 匀速高斯马尔科夫建模的方法获得更好的精度。

通过比较三种机动模式的目标轨迹的跟踪结果可以发现，在对控制量进行 建模后，虽然目标以机动方式1和机动方式3飞行时存在横向机动，但两种方 式下的跟踪精度与目标以机动方式2飞行时的跟踪精度相当，并且加速度跟踪 效果更好。这是因为目标以固定倾侧角飞行或当倾侧角变化时的飞行过程中， 虽然横向机动范围较大，导致雷达的观测误差增大，但是由气动加速度与飞行 器的倾侧角存在相互约束的关系，通过对倾侧角建模，相当于建立了气动力加 速度的描述关系，提高了目标横向机动时的跟踪精度。相关仿真研究进一步说 明对控制量特征进行建模对高超声速目标辅助跟踪的有效性。

Claims (4)

1.一种基于控制量特征辅助的临空高超声速目标跟踪方法，其特征在于，采用对控制量特征直接建模，具体包括下列步骤：

Step1.控制量建模

把攻角和倾侧角这两个控制变量看作匀速高斯马尔科夫过程，用加性白噪声来对攻角和倾侧角的变化率进行建模，表达式如下：

式中，α表示攻角，v表示倾侧角，t表示时间变量，ω_t，α和ω_t，v是未知控制输入引起的在攻角和倾侧角上一定程度的白噪声；

Step2.控制量离散化

式(1)可以进一步离散化为：

式中，T是采样周期的时间间隔，α_k+1、v_k+1分别表示攻角和倾侧角的一步预测值；认为攻角和倾侧角在离散时刻k+1和k的变化分别是时间间隔T乘以白噪声ω_α，k和ω_v，k；

Step3.控制量加入状态模型

在完成姿态角建模后，结合机动目标跟踪中状态方程扩维的思路，将攻角和倾侧角加入状态方程中进行更新，对目标状态向量进行扩维，则目标的状态变量变为11维的向量：

x(k)＝[r(k)，v(k)，a_VTC(k)，α(k)，v(k)]^T (53)

式中，k表示时间序列变量；r(k)＝[x(k)，y(k)，z(k)]^T和

分别为目标在东北天坐标系ENU中第k时刻的位置和速度向量，a_VTC(k)＝[α_V(k)，α_T(k)，α_C(k)]^T为目标在半速度VTC坐标系下的加速度，α(k)和v(k)分别是攻角和倾侧角；

Step4.目标跟踪模型

建立用于目标跟踪的离散模型为：

x_k+1＝Φ_kx_k+B_ku_k+Γ_kω_k (54)

式中，x_k为式(3)中x(k)的简化表示，x_k、x_k+1分别为第k+1和k时刻目标的状态向量，

为控制输入向量，g_xk，g_yk，g_zk为重力加速度在ENU坐标系三个方向上的分量，能够通过平面、球形或椭球三种地球模型给出，

分别为机动加速度在VTC坐标系下加速度的均值，Φ_k＝[I_3×3 TI_3×3 P₀vT₁ 0_3×2；0_3×3 I_3×3P₀vT₂ 0_3×2；0_3×3 0_3×3 T₃ 0_3×2；0_3×3 0_3×3 0_3×3 T₄]^T为状态转移矩阵，

为输入矩阵，Γ_k＝[P₀vΓ₁；P₀vΓ₂；Γ₃；Γ₄]^T为噪声输入矩阵，ω_k为协方差矩阵为Q_k的VTC坐标系下的加速度过程噪声，状态转移矩阵Φ_k和噪声输入矩阵Γ_k按式(2)进行修正，能够建立控制量姿态角辅助跟踪的离散状态方程，如式(4)所示；其中，

I_3×3＝blkdiag(1，1，1)，0_3×3＝blkdiag(0，0，0)，0_3×2＝[0 0；0 0；0 0]，

T₁＝[T₁₁η₁₁ T₁₂η₂₁ T₁₃η₃₁；T₂₁η₁₁ T₂₂η₂₁ T₂₃η₃₁；T₃₁η₁₁ T₃₂η₂₁ T₃₃η₃₁]，

T₂＝[T₁₁η₁₂ T₁₂η₂₂ T₁₃η₃₂；T₂₁η₁₂ T₂₂η₂₂ T₂₃η₃₂；T₃₁η₁₂ T₃₂η₂₂ T₃₃η₃₂]

T₃＝[η₁₃ 0 0；0 η₂₃0；0 0 η₃₃]，T₄＝[1 0；0 1]，

B₁＝[υ_VT₁₁ξ₁₁ υ_TT₁₂ξ₂₁ υ_CT₁₃ξ₃₁；υ_VT₂₁ξ₁₁ υ_TT₂₂ξ₂₁ υ_CT₂₃ξ₃₁；υ_VT₃₁ξ₁₁ υ_TT₃₂ξ₂₁ υ_CT₃₃ξ₃₁]，

B₂＝[υ_VT₁₁ξ₁₂ υ_TT₁₂ξ₂₂ υ_CT₁₃ξ₃₂；υ_VT₂₁ξ₁₂ υ_TT₂₂ξ₂₂ υ_CT₂₃ξ₃₂；υ_VT₃₁ξ₁₂ υ_TT₃₂ξ₂₂ υ_CT₃₃ξ₃₂]，

B₃＝[υ_Vξ₁₃ 0 0；0 υ_Tξ₂₃ 0；0 0 υ_Cξ₃₃]，

Γ₁＝[T₁₁ξ₁₁ T₁₂ξ₂₁ T₁₃ξ₃₁；T₂₁ξ₁₁ T₂₂ξ₂₁ T₂₃ξ₃₁；T₃₁ξ₁₁ T₃₂ξ₂₁ T₃₃ξ₃₁]，

Γ₂＝[T₁₁ξ₁₂ T₁₂ξ₂₂ T₁₃ξ₃₂；T₂₁ξ₁₂ T₂₂ξ₂₂ T₂₃ξ₃₂；T₃₁ξ₁₂ T₃₂ξ₂₂ T₃₃ξ₃₂]，

Γ₃＝[ξ₁₃ 0 0；0 ξ₂₃ 0；0 0 ξ₃₃]，Γ₄＝[T 0；0 T]，P₀为自由流压强，

为目标平均速度，

为目标在水平面的平均速度，

T₃₂＝0，

υ_V、υ_T、υ_C分别为目标在VTC坐标系三坐标轴向上的机动频率；

Step5.系统量测方程

假设系统探测中，只有雷达传感器，则系统的量测z_k为：

式中，P_k、A_k、E_k、σ分别表示测量的距离、方位角、俯仰角、量测噪声，、x、y、z为ENU坐标系中的位置分量，h(x_k)为系统非线性量测方程，ΔR、ΔA和ΔE分别为雷达的距离量测误差、方位量测误差和俯仰量测误差；

Step6.基于平方根容积卡尔曼滤波的跟踪算法

针对目标状态模型和观测模型的非线性特点，采用平方根容积卡尔曼滤波SRCKF进行目标的状态滤波；SRCKF直接对协方差矩阵的平方根进行迭代计算，避免了每步迭代的矩阵开方运算，能够克服CKF算法中滤波发散的问题，从而提高滤波的精度；

合并状态方程和量测方程为：

式中，

f(·)表示非线性状态方程，W_k-1是均值为0，方差为Q_k-1的高斯过程噪声向量，h(X_k)是非线性的量测方程，量测噪声V_k是方差为R_k的零均值高斯白噪声矢量。

2.如权利要求1所述的临空高超声速目标跟踪方法，其特征在于，SRCKF算法步骤如下：

Step1)初始化

设置状态初值x_0|0，协方差矩阵初值P_0|0，协方差矩阵平方根因子的初始值S_0|0，其中

Step2)时间更新：k＝1，2，…

(1)计算当前状态的2n个容积点X_i，k-1|k-1，i＝1，2，…，2n

其中，S_k-1|k-1为k-1时刻协方差矩阵P_k-1的平方根因子，

为k-1时刻的状态估计值，

为基本容积点，m＝2n为容积点的总数，其中n为状态维数，[1]_i为点集[1]的第i个点；

(2)计算容积点经过非线性状态转移函数的预测值，i＝1，2，…，2n：

(3)结合权值和容积点预测值，SRCKF采用相等权值，预测状态可计算为m＝2n：

(4)计算预测误差协方差矩阵的平方根因子：

S_k|k-1＝Tria([χ_k|k-1 S_Q，k-1]) (60)

其中，S_Q，k-1是过程噪声Q_k-1的平方根且有中间变量χ_k|k-1，

其中，算法S＝Tria(A)意为先对矩阵A进行QR分解，得到一个正规正交矩阵B和一个上三角矩阵C，令S＝C^T，得到的S即为上三角矩阵；

Step3)量测更新：k＝1，2，…

(1)计算更新的状态容积点，i＝1，2，…，m：

(2)计算容积点经量测非线性函数传播后的值：

Z_i，k|k-1＝h(X_i，k|k-1) (62)

(3)估计预测量测：

(4)估计新息的协方差矩阵：

S_zz，k|k-1＝Tria([Z_k|k-1S_R，k]) (64)

其中，

(5)估计互协方差矩阵：

其中，

(6)估计SRCKF的滤波增益：

(7)基于k时刻新的量测z_k，更新系统状态

(8)更新误差协方差矩阵的平方根因子

S_k|k＝Tria([χ_k|k-1-W_kZ_k|k-1 W_kS_R，k]) (68)。

3.一种基于控制量特征辅助的临空高超声速目标跟踪方法，其特征在于，采用基于姿态角的控制量特征估计，具体包括下列步骤：

利用光学传感器获取的目标姿态信息，对高超声速飞行器的攻角和倾侧角进行估计，以期进行一步提高目标的跟踪精度；具体处理流程如下：

STEP1.在弹体坐标系中计算目标的攻角和倾侧角

在弹体坐标系与速度坐标系中计算攻角α和倾侧角v的方法如下：

式中，v_dx、v_dy、v_dz是速度矢量在弹体坐标系上的分量，

为目标的飞行速度；

STEP2.根据光学传感器获取的目标姿态角信息估计攻角和倾侧角

假设光学传感器可以测量地理坐标系与弹体坐标系的三个欧拉角，分别为偏航角ψ、俯仰角ε和滚转角γ，则由地理坐标系到弹体坐标系的变换矩阵M_D为：

则地理坐标系的速度矢量与弹体坐标系的速度矢量的关系为：

式中，v_E、v_N、v_U分别表示目标在地理坐标系中速度分量；

由此，根据式(19)、(20)和式(21)建立由光学传感器观测的目标姿态角信息估计的目标攻角和倾侧角；

STEP3.对攻角和倾侧角的递推方程离散化

由于NSHT的动力学模型为时域系统，对目标进行跟踪时，需要将其离散化，离散化后，将控制变量的递推方程离散化为：

式中，α_k+1、v_k+1分别表示攻角和倾侧角的一步预测值，T为采样周期；认为攻角和倾侧角在离散时刻k和k-1的变化分别是时间间隔T乘以零均值的白噪声ω_α，k和ω_v，k；

STEP4.在目标的状态变量中增加攻角和倾侧角

建立两个控制量的递推方程后，把通过姿态角估计的两个控制量加入目标的状态方程，则状态方程修正为：

x(k)＝[r(k)，v，(k)，a_VTC(k)，α(k)，v(k)]^T (73)

式中各变量的含义与前面相同，只是攻角α(k)和倾侧角v(k)是基于姿态角的估计值；

STEP5.系统的量测模型

系统量测模型中除了距离、方位和俯仰测量信息，还有基于光学传感器获取的目标的姿态信息的量测模型为：

式中，在观测向量z_k中，前三项分别为距离R_k、方位A_k和俯仰角E_k，ψ_k、θ_k和γ_k分别为光学传感器观测的偏航角、俯仰角和滚转角；各量测噪声均是已知的独立白噪声，量测噪声方差可以表示为

式中，

和

分别是雷达量测的距离、方位、俯仰和光学传感器观测的航角、俯仰角和滚转角方差；

STEP6.基于平方根容积卡尔曼滤波的跟踪算法

根据上面的非线性量测方程和非线性状态方程，采用SRCKF方法对目标位置、速度、气动加速度和控制量的状态进行递推滤波。

4.如权利要求3所述的临空高超声速目标跟踪方法，其特征在于，SRCKF算法步骤如下：

Step1初始化

设置状态初值x_0|0，协方差矩阵初值P_0|0，协方差矩阵平方根因子的初始值S_0|0，其中

Step2时间更新：k＝1，2，…

(1)计算当前状态的2n个容积点X_i，k-1|k-1，i＝1，2，…，2n

其中，S_k-1|k-1为k-1时刻协方差矩阵P_k-1的平方根因子，

为k-1时刻的状态估计值，

为基本容积点，m＝2n为容积点的总数，其中n为状态维数，[1]_i为点集[1]的第i个点；

(2)计算容积点经过非线性状态转移函数的预测值，i＝1，2，…，2n：

(3)结合权值和容积点预测值，SRCKF采用相等权值，预测状态可计算为m＝2n：

(4)计算预测误差协方差矩阵的平方根因子：

S_k|k-1＝Tria([χ_k|k-1 S_Q，k-1]) (79)

其中，S_Q，k-1是过程噪声Q_k-1的平方根且有中间变量χ_k|k-1，

其中，算法S＝Tria(A)意为先对矩阵A进行QR分解，得到一个正规正交矩阵B和一个上三角矩阵C，令S＝C^T，得到的S即为上三角矩阵；

Step3量测更新：k＝1，2，…

(1)计算更新的状态容积点，i＝1，2，…，m：

(2)计算容积点经量测非线性函数传播后的值：

Z_i，k|k-1＝h(X_i，k|k-1) (81)

(3)估计预测量测：

(4)估计新息的协方差矩阵：

S_zz，k|k-1＝Tria([Z_k|k-1S_R，k]) (83)

其中，

(5)估计互协方差矩阵：

其中，

(6)估计SRCKF的滤波增益：

(7)基于k时刻新的量测z_k，更新系统状态

(8)更新误差协方差矩阵的平方根因子

S_k|k＝Tria([χ_k|k-1-W_kZ_k|k-1W_kS_R，k]) (87)。

Images