作为一种兼具战略威慑和战术打击能力的新型武器, 临近空间高超声速飞行器(near space hypersonic vehicle, NSHV)具有飞行高度高、飞行速度快、突防能力强、打击威胁大等特点, 通常飞行于20~100 km的临近空域^[1-2], 飞行速度大于5 Ma。近年来, 以美俄为首的军事大国加大开展高超声速飞行器核心技术和攻防体系的研究, 逐步迈入工程研制阶段, 并不断有印、法、日、英等国家加入NSHV的研制行列^[3-4]。可以预见, 未来几年内，NSHV技术领域的竞争将越发激烈, 各国空天安全将面临新的挑战, 对NSHV的早期预警、稳定跟踪和防御拦截成为迫切需求。因此, 为解决NSHV难以进行三维稳定跟踪的问题, 进一步对其跟踪技术及理论进行研究尤为必要。

对NSHV跟踪的问题实质上是利用最新量测信息不断对系统预测状态进行修正的过程, 主要包括两个核心内容: 机动模型和滤波算法, 其中机动模型决定了系统预测状态的发展趋势, 滤波算法用来修正系统的预测误差^[5-7]。由于NSHV高机动的特性, 采用常规卡尔曼滤波在跟踪时，易发生精度下降甚至滤波发散的现象。因此, 建立符合NSHV运动特性的机动模型和设计更为鲁棒的滤波算法是提高NSHV跟踪精度、实现稳定跟踪的可行路径。在机动模型方面, 单一的经典机动模型如匀速(constant velocity, CV)模型、匀加速(constant acceleration, CA)模型、当前统计(current statistical, CS)模型、Singer模型、Jerk模型等难以进行稳定有效的跟踪^[8-10], 系统预测误差较大。文献[11]以CS模型为中心, 结合CV模型, 设计了一种修正变结构交互多模型(interactive multiple model, IMM)算法, 使模型切换更为灵活。文献[12]结合NSHV目标的周期性滑跃特点, 对加速度的时间相关性进行建模, 设计了一种针对NSHV特性的跟踪算法。文献[13]分析了NSHV目标的机动特性, 将加速度建模为衰减震荡自相关的随机过程, 并分析了模型的适应性。在滤波算法方面, 以常规卡尔曼滤波算法为基础, 引入渐消因子，从而实时调整过程噪声协方差，克服滤波发散^[14-16]。文献[17]分析了造成Kalman滤波发散的原因, 并对算法的一步预测值进行了修正, 提高了滤波稳定性。文献[18]对传统的衰减记忆算法进行了改进, 提出了一种可以智能衰减陈旧信息, 合理补偿当前信息的强跟踪算法, 改善了跟踪效果。考虑NSHV目标的机动特性, 单一的从机动模型或滤波算法的角度出发，实现可靠跟踪较为困难, 不能适应其复杂机动特点。

基于此, 本文首先分析了NSHV目标的机动特性, 将其三维运动轨迹分解到纵向弹道和横向弹道, 在两个弹道上分别建立适合其运动特征的跟踪模型。其次，基于量测噪声协方差和残差设定了修正门限, 引入多重渐消因子，实时调节各通道的权重, 并设计了一种限定记忆长度的指数加权新息协方差估计方法，对陈旧量测信息进行限定记忆, 提高新近量测数据的影响。最后通过仿真实验对算法进行验证。

NSHV目标的典型机动模式主要有平衡滑翔和跳跃滑翔两种^[19-20], 其中平衡滑翔是目标飞行过程中在纵向平面内受到的升力与重力和离心力达到平衡时所处的状态, 当其纵向平面受力不满足平衡条件时, 目标会偏离平衡滑翔的加速度, 形成跳跃滑翔弹道, 其纵向平面的轨迹如图 1所示。在平衡滑翔时目标状态较为稳定, 而在跳跃滑翔时目标的跟踪相对更为困难, 因此, 本文主要研究对象为跳跃滑翔NSHV目标。

攻角和倾侧角是NSHV目标的主要控制参量。其中, 攻角主要影响目标纵向弹道的机动形式, 其设定受飞行走廊的约束, 通常采用线性、常值或分段函数进行描述：

(1) $\frac{L \cos \upsilon}{M}+\left(\frac{v^{2}}{r}-g\right) \cos \theta=0$

式中: L为目标所受升力; M为目标质量; υ为倾侧角; θ为速度倾角; r为地心距; g为重力加速度。

攻角取值在满足式(1)的条件下, 目标处于平衡滑翔状态, 在其他取值情况下, 目标处于跳跃滑翔状态^[21]。而倾侧角对横纵向弹道都有影响, 决定了目标横向弹道的机动形式, 其设定同样受飞行走廊的约束, 主要是为了规避探测、机动突防, 可以进行类似大C或S型横向机动。倾侧角的设定主要是有利于进攻方进行打击摧毁, 受敌方实际战场环境的影响较大, 但复杂的横向机动会增大被探测的概率, 降低攻击时效。

基于上述分析, 将NSHV目标空间机动弹道分为纵向机动弹道和横向机动弹道。在纵向弹道上, 根据加速度衰减震荡特性, 将其建模为二阶时间自相关的零均值随机过程^[13]; 在横向弹道上, 其加速度变化相对简单, 采用Singer模型和CA模型进行IMM滤波。最后通过并行滤波和合并处理形成完整的NSHV目标三维跟踪。

根据NSHV目标纵向方向上加速度的衰减震荡特性, 将加速度建模为二阶时间相关函数^[13], 即

(2) $R_{a}(\tau)=\mathrm{E}[a(t+\tau) a(t)]=\sigma_{a}^{2} \underbrace{\mathrm{e}^{-\alpha|\tau|}}_{\mathrm{I}} \underbrace{\cos (\beta \tau)}_{\mathrm{II}}$

式中：第Ⅰ部分体现了加速度的指数衰减特性；第Ⅱ部分体现了加速度的周期特性；τ为时间差；σ_a²为加速度方差；α为最大相关衰减量；β=2π/ΔT_c为机动震荡频率；ΔT_c为震荡周期。

该相关函数的功率谱密度为

(3) $\begin{gathered}S(w)=\frac{\sigma_{a}^{2}}{2}\left\{\left[\frac{1}{\alpha-(\beta-w) \mathrm{j}}+\frac{1}{\alpha+(\beta+w) \mathrm{j}}\right]+\right. \\\left.\left[\frac{1}{\alpha+(\beta-w) \mathrm{j}}+\frac{1}{\alpha-(\beta+w) \mathrm{j}}\right]\right\}\end{gathered}$

式中: w为角速率。

令G(w)=[α-(β-w)j][α+(β+w)j], 则有

(4) $\begin{gathered}S(w)=2 \alpha \sigma_{a}^{2} \frac{\left(\sqrt{\alpha^{2}+\beta^{2}}+w \mathrm{j}\right)}{G(w)} \frac{\left(\sqrt{\alpha^{2}+\beta^{2}}-w \mathrm{j}\right)}{G(-w)}= \\W(w \mathrm{j}) H(w \mathrm{j}) H(-w \mathrm{j})\end{gathered}$

式中：W(wj)为白噪声的傅里叶变换。可得到加速度a(t)的微分方程为:

(5) $\ddot{a}(t)+2 \dot{\alpha}(t)+\left(\alpha^{2}+\beta^{2}\right) a(t)=\dot{w}(t)+\sqrt{\alpha^{2}+\beta^{2}} w(t)$

式中: w(t)是方差为2ασ_a²的零均值高斯白噪声, 其状态方程为

(6) $\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d} t}\left[\begin{array}{l}x(t) \\\dot{x}(t) \\a(t) \\\dot{a}(t)\end{array}\right]=\left[\begin{array}{cccc}0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 \\0 & 0 & -\left(\alpha^{2}+\beta^{2}\right) & -2 \alpha\end{array}\right]\left[\begin{array}{l}x(t) \\\dot{x}(t) \\a(t) \\\dot{a}(t)\end{array}\right]+\\\left[\begin{array}{cc}0 & 0 \\0 & 0 \\0 & 0 \\\sqrt{\alpha^{2}+\beta^{2}} & 1\end{array}\right]\left[\begin{array}{c}w(t) \\\dot{w}(t)\end{array}\right]\end{array}$

对式(6)进行离散化, 则其状态转移矩阵为

(7) $\boldsymbol{F}_{D O}(T, \alpha, \beta)=\left[\begin{array}{llll}1 & T & f_{1}(T) & f_{2}(T) \\0 & 1 & f_{3}(T) & f_{4}(T) \\0 & 0 & f_{5}(T) & f_{6}(T) \\0 & 0 & f_{7}(T) & f_{8}(T)\end{array}\right]$

式中：

(8) $\left\{\begin{aligned}f_{1}(T)=& a_{2}(T) \\f_{2}(T)=& a_{3}(T) \\f_{3}(T)=& T-\left(\alpha^{2}+\beta^{2}\right) a_{3}(T) \\f_{4}(T)=& a_{2}(T)-2 \alpha a_{3}(T) \\f_{5}(T)=& 1-\left(\alpha^{2}+\beta^{2}\right) a_{2}(T)+2 \alpha\left(\alpha^{2}+\beta^{2}\right) a_{3}(T) \\f_{6}(T)=& T-2 \alpha a_{2}(T)+\left(3 \alpha^{2}-\beta^{2}\right) a_{3}(T) \\f_{7}(T)=&\left(\alpha^{2}+\beta^{2}\right) T+2 \alpha\left(\alpha^{2}+\beta^{2}\right) a_{2}(T)-\\&\left(\alpha^{2}+\beta^{2}\right)\left(3 \alpha^{2}-\beta^{2}\right) a_{3}(T) \\f_{8}(T)=& 1-2 \alpha T+\left(3 \alpha^{2}-\beta^{2}\right) a_{2}(T)-\\& 4 \alpha\left(\alpha^{2}-\beta^{2}\right) a_{3}(T)\end{aligned}\right.$

(9) $\left\{\begin{array}{l}a_{0}(t)=1 \\a_{1}(t)=t \\a_{2}(t)=\frac{\left(\alpha^{3}-3 \alpha \beta^{2}\right) \sin \beta t+\left(3 \alpha^{2} \beta-\beta^{3}\right) \cos \beta t}{\beta\left(\alpha^{2}+\beta^{2}\right)^{2}} \mathrm{e}^{-\alpha t}+ \\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{2 \alpha}{\alpha^{2}+\beta^{2}} a_{1}(t)-\frac{3 \alpha^{2}-\beta^{2}}{\left(\alpha^{2}+\beta^{2}\right)^{2}} a_{0}(t) \\a_{3}(t)=\frac{\left(\alpha^{2}-\beta^{2}\right) \sin \beta t+2 \alpha \beta \cos \beta t}{\beta\left(\alpha^{2}+\beta^{2}\right)^{2}} \mathrm{e}^{-\alpha t}+ \\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{1}{\alpha^{2}+\beta^{2}} a_{1}(t)-\frac{2 \alpha}{\left(\alpha^{2}+\beta^{2}\right)^{2}} a_{0}(t)\end{array}\right.$

其过程噪声为

(10) $\begin{gathered}\boldsymbol{W}(k)=\int_{k T}^{(k+1) T} \mathrm{e}^{A((k+1) T-\tau)} \boldsymbol{B} \boldsymbol{w}(\tau) \mathrm{d} \tau= \\\int_{k T}^{(k+1) T}\left[\begin{array}{ll}\sqrt{\alpha^{2}+\beta^{2}} f_{2}(\tau) & f_{2}(\tau) \\\sqrt{\alpha^{2}+\beta^{2}} f_{4}(\tau) & f_{4}(\tau) \\\sqrt{\alpha^{2}+\beta^{2}} f_{6}(\tau) & f_{6}(\tau) \\\sqrt{\alpha^{2}+\beta^{2}} f_{8}(\tau) & f_{8}(\tau)\end{array}\right] \boldsymbol{w}(\tau) \mathrm{d} \tau\end{gathered}$

过程噪声协方差Q_DO(k)=E(W(k)W^T(k)), 具体详见文献[13]。

在横向弹道上, 其加速度变化相对简单, 滤波器采用CA模型和Singer模型。CA模型将加速度建模为白噪声, 其离散化时间状态^[5]方程为

(11) $\boldsymbol{x}(k+1)=\boldsymbol{F}_{\mathrm{CA}} \boldsymbol{x}(k)+\boldsymbol{w}_{\mathrm{CA}}(k)$

式中：状态转移矩阵为

(12) $\boldsymbol{F}_{\mathrm{CA}}=\left[\begin{array}{ccc}1 & T & T^{2} / 2 \\0 & 1 & T \\0 & 0 & 1\end{array}\right]$

协方差矩阵Q_CT为

(13) $\boldsymbol{Q}_{\mathrm{CT}}=S_{w}\left[\begin{array}{ccc}T^{5} / 20 & T^{4} / 8 & T^{3} / 6 \\T^{4} / 8 & T^{3} / 3 & T^{2} / 2 \\T^{3} / 6 & T^{2} / 2 & T\end{array}\right]$

式中：S_w为白噪声的功率谱密度。

Singer模型由Singer在1970年提出, 将加速度建模为具有指数自相关的零均值随机过程^[22], 自相关量为

(14) $R_{a}(\tau)=\mathrm{E}(a(t+\tau) a(t))=\sigma_{a}^{2} \mathrm{e}^{-\alpha|\tau|}$

加速度微分方程为

(15) $\dot{a}(t)=-\alpha a(t)+w(t)$

此处α为的机动时间常数的倒数, 即机动频率, 与纵向弹道运动模型中的α含义类似。其离散化时间状态方程为

(16) $\boldsymbol{x}(k+1)=\boldsymbol{F}_{\text {Singer }} \boldsymbol{x}(k)+\boldsymbol{w}_{\text {Singer }}(k)$

式中：

(17) $\boldsymbol{F}_{\text {Singer }}=\left[\begin{array}{ccc}1 & T & \frac{\alpha T-1+\mathrm{e}^{-\alpha T}}{\alpha^{2}} \\0 & 1 & \frac{1-\mathrm{e}^{-\alpha T}}{\alpha} \\0 & 1 & \mathrm{e}^{-\alpha T}\end{array}\right]$

其过程噪声协方差矩阵参见文献[5]。

根据第1节分析可知, 跳跃滑翔NSHV目标的加速度具有震荡特性, 其机动样式相对灵活多变, 单纯采用常规卡尔曼滤波算法进行滤波易发生跟踪精度下降。因此, 本文设计了一种基于多重渐消因子的自适应卡尔曼滤波方法, 增强滤波算法的自适应跟踪能力。

设采样周期为T, 离散系统的时间状态方程和量测方程可表示为

(18) $\boldsymbol{X}(k)=\boldsymbol{F}(k-1) \boldsymbol{X}(k-1)+\boldsymbol{\varGamma}(k-1) \boldsymbol{W}(k-1)$

(19) $\boldsymbol{Z}(k)=\boldsymbol{H}(k) \boldsymbol{X}(k)+\boldsymbol{V}(k)$

式中: X(k)和Z(k)分别表示目标k时刻的状态向量和量测向量; F(k)表示k时刻的状态转移矩阵; H(k)表示k时刻的量测矩阵; Γ(k-1)为k时刻的系统噪声矩阵; W(k)和V(k)分别表示k时刻的状态噪声和量测噪声。

典型的离散时间卡尔曼滤波过程可表示为

(20) $\hat{\boldsymbol{X}}(k+1 / k)=\boldsymbol{F}(k) \hat{\boldsymbol{X}}(k / k)$

(21) $\boldsymbol{P}(k+1 / k)=\boldsymbol{F}(k) \boldsymbol{P}(k / k) \boldsymbol{F}^{\mathrm{T}}(k)+\boldsymbol{\varGamma}(k) \boldsymbol{Q}(k) \boldsymbol{\varGamma}^{\mathrm{T}}(k)$

(22) $\boldsymbol{r}(k+1)=\boldsymbol{Z}(k+1)-\boldsymbol{H}(k+1) \hat{\boldsymbol{X}}(k+1 / k)$

(23) $\boldsymbol{S}(k+1)=\boldsymbol{H}(k+1) \boldsymbol{P}(k+1 / k) \boldsymbol{H}^{\mathrm{T}}(k+1)+\boldsymbol{R}(k+1)$

(24) $\boldsymbol{K}(k+1)=\boldsymbol{P}(k+1 / k) \boldsymbol{H}^{\mathrm{T}}(k+1) \boldsymbol{S}^{-1}(k+1)$

(25) $\hat{\boldsymbol{X}}(k+1 / k+1)=\hat{\boldsymbol{X}}(k+1 / k)+\boldsymbol{K}(k+1) \boldsymbol{r}(k+1)$

(26) $\boldsymbol{P}(k+1 / k+1)=\boldsymbol{P}(k+1 / k)-\boldsymbol{K}(k+1) \boldsymbol{S}(k+1) \boldsymbol{K}^{\mathrm{T}}(k+1)$

在给定目标状态向量初始值和估计状态向量的协方差初始矩阵后, 滤波便可以根据最新量测信息进行持续递推。

一般而言, 造成卡尔曼滤波精度下降或发散的原因主要有3种: 一是建立的目标机动模型不准确; 二是目标运动状态发生突变; 三是滤波迭代过程中估计误差的传递^[17]。对卡尔曼滤波进行修正的本质就是利用最新量测值对预测值进行修正, 减小估计误差。而当预测值偏离实际值较大时, 增大修正的比重, 防止滤波发散, 这就是传统渐消卡尔曼滤波的基本思路。但是这种修正实际上依赖量测值的准确性, 也就是当量测值本身误差非常大时, 这种修正的意义非常有限。同时, 传统渐消卡尔曼滤波的单一渐消因子在处理多通道滤波时, 对各通道的调节能力相同, 导致滤波精度降低。因此, 本文设计了一种多重渐消因子的自适应卡尔曼滤波方法, 能自适应调节各通道的过程噪声协方差权重, 并设定了发散程度判断门限, 降低量测误差的影响, 提高滤波器整体性能。

由于量测值本身是存在误差的, 为降低量测误差的影响, 保证滤波器的优良性能, 在对预测值进行修正之前设定修正门限。只有当预测值与量测值的误差大于修正门限时, 才启动多重渐消自适应滤波, 否则继续执行常规卡尔曼滤波。这里, 用r(k)衡量预测值与量测值的误差, 用量测误差协方差的均方值表示量测值最大误差, 门限设定为$ {\rm{gate}} = \chi \sqrt {R{_{i, i}}\left( k \right)} $, 其中χ为门限系数, 取值范围一般为0~1, R_{i, i}(k)表示k时刻量测向量第i行第i列的量测噪声。本文将门限系数χ设定为0.5。

式(21)中如果引入标量渐消因子ρ, 则可表示为

(27) $\boldsymbol{P}(k+1 / k)=\rho(k) \boldsymbol{F}(k) \boldsymbol{P}(k / k) \boldsymbol{F}^{\mathrm{T}}(k)+\boldsymbol{\varGamma}(k) \boldsymbol{Q}(k) \boldsymbol{\varGamma}^{\mathrm{T}}(k)$

式中：ρ(k)≥1, 则标准卡尔曼滤波将变为单渐消因子卡尔曼滤波^[23]。

标量渐消因子ρ的确定方法为

(28) $\left\{\begin{array}{l}\rho(k)=\max \{1, \operatorname{tr}(\boldsymbol{N}(k)) / \operatorname{tr}(\boldsymbol{M}(k))\} \\\boldsymbol{N}(k)=\hat{\boldsymbol{C}}(k)-\boldsymbol{H}(k) \boldsymbol{\varGamma}(k-1) \boldsymbol{Q}(k-1) \boldsymbol{\varGamma}^{\mathrm{T}}(k-1) \boldsymbol{H}^{\mathrm{T}}(k)-\boldsymbol{R}(k) \\\boldsymbol{M}(k)=\boldsymbol{H}(k) \boldsymbol{\varGamma}(k-1) \boldsymbol{P}(k-1) \boldsymbol{\varGamma}^{\mathrm{T}}(k-1) \boldsymbol{H}^{\mathrm{T}}(k)\end{array}\right.$

式中: $ \mathit{\boldsymbol{\hat C}}(k)$为新息协方差估计值, 由开窗法得到其估计值^[24]为

(29) $\hat{\boldsymbol{C}}(k)=\frac{1}{N} \sum\limits_{j=0}^{N-1}\left[\boldsymbol{r}(k-j) \boldsymbol{r}^{\mathrm{T}}(k-j)\right]$

式中: N为滑动窗口长度。

从式(27)~式(29)可以看出, 当处理多通道滤波时, 单渐消因子ρ(k)对各通道滤波的调节能力相同。为避免这种缺陷, 本文引入多重渐消因子ρ(k), 能根据各通道的实际滤波情况对各通道进行实时调节。

(30) $\boldsymbol{\rho}(k)=\operatorname{diag}\left\{\rho_{1}, \rho_{2}, \cdots, \rho_{i}, \cdots, \rho_{n}\right\}$

式中: n为状态向量维数; ρ(i)表示第i维状态向量对应的渐消因子。

(31) $\boldsymbol{P}(k+1 / k)=\boldsymbol{\rho}(k) \boldsymbol{F}(k) \boldsymbol{P}(k / k) \boldsymbol{F}^{\mathrm{T}}(k)+\boldsymbol{\varGamma}(k) \boldsymbol{Q}(k) \boldsymbol{\varGamma}^{\mathrm{T}}(k)$

定义矩阵J(k)= Γ(k-1)P(k-1)Γ^T(k-1), 则有M(k)=H(k)J(k)H^T(k)。而量测矩阵H(k)中的元素通常为0或1两种取值, 故由ρ(k)M(k)=N(k), 可得出多重渐消因子的计算方式为：

(32) $\rho_{i}(k)=\left\{\begin{array}{l}\max \left\{1, \frac{N_{i, i}(k)}{J_{i, i}(k)}\right\}, H_{i}(k)=1 \\1, \text { 其他 }\end{array}\right.$

式中: H_i(k)表示k时刻量测矩阵中的第i个元素。

当ρ(k)=I时, 多重渐消滤波退化为常规卡尔曼滤波。从式(32)可以分析, ρ_i(k)的取值受新息协方差估计值$ \mathit{\boldsymbol{\hat C}}(k)$影响较大, 因此$ \mathit{\boldsymbol{\hat C}}(k)$估计值直接影响ρ_i(k)的调节效果。而$ \mathit{\boldsymbol{\hat C}}(k)$又受量测值的影响。为尽可能消除量测误差给滤波带来的不稳定性, 可以根据应用场景的实际情况, 给ρ_i(k)设定取值上限, 避免滤波过于震荡。

由式(29)可以看出, 经开窗法计算得到的新息协方差估计值$ \mathit{\boldsymbol{\hat C}}(k)$是对历史数据的算术平均, 不能增加新近量测信息的权重值。为增加最新量测数据的权值, 使陈旧数据随时间的拉长而所占比重越来越低, 设计了一种限定记忆长度的指数加权新息协方差估计方法。

设限定新息协方差序列的记忆长度为N(N≤k), 加权系数为χ (i), 在时刻k应满足

(33) $\left\{\begin{array}{l}\chi(i) b=\chi(i-1) \\\sum\limits_{i=k-N+1}^{k} \chi(i)=1\end{array}\right.$

式中: b为遗忘因子, 通常取值范围为0.7~0.95。确定遗忘因子b后, 根据数列性质, 有:

(34) $\left\{\begin{array}{l}b^{0}+b^{1}+b^{2}+\cdots+b^{N-1}=\frac{1-b^{N}}{1-b} \\\left(b^{0}+b^{1}+b^{2}+\cdots+b^{N-1}\right) \frac{1-b}{1-b^{N}}=1\end{array}\right.$

因此, $ \chi {\rm{ }}\left( i \right) = \frac{{1{\rm{ - }}b}}{{1{\rm{ - }}{b^N}}}{b^{k - i}}$。

下面推导新息协方差的递推公式。k时刻的新息协方差估计$ \mathit{\boldsymbol{\hat C}}(k)$为

(35) $\begin{aligned}\hat{\boldsymbol{C}}(k)=& \sum\limits_{i=k-N+1}^{k} \chi(i) \boldsymbol{r}(i) \boldsymbol{r}^{\mathrm{T}}(i)={\frac{1-b}{1-b^{N}} \sum\limits_{i=k-N+1}^{k} b^{k-i} \boldsymbol{r}(i) \boldsymbol{r}^{\mathrm{T}}(i)=} \\&{\frac{1-b}{1-b^{N}}[\boldsymbol{r}(k) \boldsymbol{r}^{\mathrm{T}}(k)+\underbrace{\sum\limits_{i=k-N+1}^{k-1} b^{k-i} \boldsymbol{r}(i) \boldsymbol{r}^{\mathrm{T}}(i)}_{I}]}\end{aligned}$

k-1时刻的新息协方差估计$ \mathit{\boldsymbol{\hat C}}(k-1)$为

(36) $\begin{gathered}\hat{\boldsymbol{C}}(k-1)=\sum\limits_{i=k-N}^{k-1} \chi(i) \boldsymbol{r}(i) \boldsymbol{r}^{\mathrm{T}}(i)= \\\frac{1-b}{1-b^{N}} \sum\limits_{i=k-N}^{k-1} b^{k-i-1} \boldsymbol{r}(i) \boldsymbol{r}^{\mathrm{T}}(i)= \\\frac{1-b}{1-b^{N}} \frac{1}{b}\left[b^{N} \boldsymbol{r}(k-N) \boldsymbol{r}^{\mathrm{T}}(k-N)+\sum\limits_{i=k-N+1}^{k-1} b^{k-i} \boldsymbol{r}(i) \boldsymbol{r}^{\mathrm{T}}(i)\right]\end{gathered}$

从式(35)和式(36)可以看出, 式(35)的Ⅰ项可以用$ \mathit{\boldsymbol{\hat C}}(k-1)$来表示, 由此可得到新息协方差估计值的递推公式:

(37) $\hat{\boldsymbol{C}}(k)=b \hat{\boldsymbol{C}}(k-1)+\frac{1-b}{1-b^{N}}\left[\boldsymbol{r}(k) \boldsymbol{r}^{\mathrm{T}}(k)-b^{N} \boldsymbol{r}(k-N) \boldsymbol{r}^{\mathrm{T}}(k-N)\right]$

由式(37)可得出, 新息协方差估计值的大小依赖于记忆长度N和遗忘因子b的取值, 其中N和b取值的选取主要取决于系统状态的变化剧烈程度, 对于变化越平缓的系统, 记忆长度N的取值应越大, 遗忘因子b的取值同样也应越大, 以此得到更准确的新息协方差估计值。

在跟踪过程中, 雷达量测数据通常在球坐标系下建立, 而卡尔曼滤波器在笛卡尔坐标系中工作。为解决这一问题, 本文采用修正无偏量测转换^[25-26]将球坐标系下的量测数据转换到笛卡尔坐标系中, 得到其转换后的量测值为

(38) $\left\{\begin{array}{l}x=\lambda_{E}^{-1} \lambda_{A}^{-1} r_{m} \cos A \cos E \\y=\lambda_{E}^{-1} \lambda_{A}^{-1} r_{m} \sin A \cos E \\z=\lambda_{E}^{-1} r_{m} \sin E\end{array}\right.$

式中: A为方位角; E为俯仰角; r_m为目标距离; λ_E=e^-σ_E2/2, σ_E²为俯仰角误差平方, λ_A=e^-σ_A2/2, σ_A²为方位角误差平方。量测噪声协方差R为

(39) $\boldsymbol{R}=\left[\begin{array}{ccc}r_{x x} & r_{x y} & r_{x z} \\r_{y x} & r_{y y} & r_{y z} \\r_{z x} & r_{z y} & r_{z z}\end{array}\right]$

式中: R为对称矩阵, 修正无偏量测转换后R内各元素为

(40) $\left\{\begin{aligned}r_{x x}=& {\left[\left(\lambda_{A} \lambda_{E}\right)^{-2}-2\right] r_{m}^{2} \cos ^{2} A \cos ^{2} E+\frac{1}{4}\left(r_{m}^{2}+\sigma_{r}^{2}\right) \cdot } \\& {\left[1+\lambda^{\prime}{}_{A} \cos (2 A)\right]\left[1+\lambda^{\prime}{}_{E} \cos (2 E)\right] } \\r_{x y}=& {\left[\left(\lambda_{A} \lambda_{E}\right)^{-2}-2\right] r_{m}^{2} \sin A \cos A \cos ^{2} E+\frac{1}{4}\left(r_{m}^{2}+\sigma_{r}^{2}\right) \cdot } \\& \lambda^{\prime}{}_{A} \sin (2 A)\left(1+\lambda^{\prime}{}_{E} \cos 2 E\right)\\r_{x z}=&\left(\lambda_{E}^{-2} \lambda_{A}^{-1}-\lambda_{A}^{-1}-\lambda_{A}\right) r_{m}^{2} \cos A \sin E \cos E+\frac{1}{2}\left(r_{m}^{2}+\sigma_{r}^{2}\right)\cdot \\& \lambda_{A} \lambda^{\prime}{}_{E} \sin A \cos (2 E) \\r_{y y}=& {\left[\left(\lambda_{A} \lambda_{E}\right)^{-2}-2\right] r_{m}^{2} \sin ^{2} A \cos ^{2} E+\frac{1}{4}\left(r_{m}^{2}+\sigma_{r}^{2}\right) \cdot } \\& {\left[1-\lambda^{\prime}{}_{A} \cos (2 A)\right]\left[1+\lambda^{\prime}{}_{E} \cos (2 E)\right] }\\r_{y z}=&\left(\lambda_{E}^{-2} \lambda_{A}^{-1}-\lambda_{A}^{-1}-\lambda_{A}\right) r_{m}^{2} \sin A \sin E \cos E+\frac{1}{2}\left(r_{m}^{2}+\sigma_{r}^{2}\right) \cdot \\& \lambda_{A} \lambda^{\prime}{ }_{E} \sin A \sin (2 E) \\r_{z z}=&\left(\lambda_{E}^{-2}-2\right) r_{m}^{2} \sin ^{2} E+\frac{1}{2}\left(r_{m}^{2}+\sigma_{r}^{2}\right)\left[1-\lambda^{\prime}{}_{A} \cos (2 A)\right]\end{aligned}\right.$

式中: λ′_A=e^-2σ_A2=λ_A⁴; λ′_E=e^-2σ_E2=λ_E⁴; σ_r²为测距误差平方。

参考美国官方公布的CAV-H飞行试验, 设置如下仿真场景: 高超声速滑翔飞行器目标由发射器助推到50 km高度并释放后进行滑翔, 其再入速度为15 Ma, 攻角和倾侧角模式如图 3所示。

攻角模式表达式为

(41) $\alpha_{c}=\left\{\begin{array}{l}\alpha_{\max }, V>V_{1} \\\alpha_{\operatorname{mid}}+\alpha_{\text {bal }} \sin \left[\frac{{\rm{ \mathsf{ π} }}}{V_{1}-V_{2}}\left(V-V_{\text {mid }}\right)\right], V_{2}<V \leqslant V_{1} \\\alpha_{\max -K}, V \leqslant V_{2}\end{array}\right.$

式中: v_mid=(v₁+v₂)/2;α_mid=(α_max+α_max-K)/2;α_bal=(α_max-α_max-K)/2;α_max=25°; α_{max_K}=11°; v₁=4 000 m/s; v₂=2 000 m/s。

假定观测雷达地理坐标为[12°, 1.5°, 1 km], 采样间隔为0.1 s, 距离误差为100 m, 俯仰角和方位角误差均为0.1°, 雷达观测过程中无地面遮挡, 目标与雷达距离约为300 km。

为验证本文算法性能, 设置两个仿真实验, 实验1主要验证本文算法相对其他算法跟踪NSHV目标时的精度, 实验2主要验证本文算法在不同量测参数下的跟踪稳定性。

选取CA+CV+CS的IMM滤波算法、文献[12]所提的跟踪算法以及基于本文运动模型的常规卡尔曼滤波算法作为对比, 分析各模型的跟踪效果。其中, 目标在各坐标轴对应的加速度方差σ_a²=10², 加加速度方差$ \sigma _{\dot a}^2 = {1^2}$, 各模型的具体参数设置如下:

(1) CA+CV+CS的IMM滤波算法: 机动频率为1/20, 最大加速度a_max=10g, 模型先验概率u₀=[0.1，0.8，0.1];

(2) 文献[12]所提跟踪算法: 角速率w₀=0.05;

(3) 本文运动模型: 机动震荡频率β=0.06 rad/s, 最大相关衰减量α=1/300。

蒙特卡罗次数为100, 得到各模型对应的位置均方根误差(root mean squre errer, RMSE)如图 4所示, 速度估计和速度RMSE分别如图 5和图 6所示。

从位置跟踪误差可以看出, 由于CA+CV+CS的IMM滤波算法对目标运动模式建模不够准确, 在50~150 s的时间区间, 其纵向和横向弹道跟踪误差都存在明显起伏, 且收敛速度较慢。文献[12]提出的跟踪算法在跟踪类周期运动目标时能保持较好的性能, 能快速收敛并保持较为稳定的跟踪效果, 没有明显剧烈的误差起伏, 但由于NSHV机动并不是严格周期性运动, 所以其总体跟踪误差较本文运动模型稍差。本文提出的运动模型能较准确地描述NSHV的运动模式, 所以其跟踪效果较好, 但由于常规卡尔曼滤波器对误差的修正有限, 所以在0~80 s的跟踪阶段, 其收敛速度要逊于本文算法。

从速度误差来看, 本文算法收敛速度最快, 且整体较为平稳, 精度最高。由于NSHV目标的机动特性, CA+CV+CS的IMM滤波算法速度跟踪误差起伏最大, 其他2种算法整体差别不明显, 本文运动模型的常规卡尔曼滤波算法估计误差精度略优。

这里需要说明的是, 从图 4和图 6中可以看出, 在进行到280 s以后, 所有算法的误差都有所增加, 导致这一现象的主要有两个可能原因: 第一, 从图 5可以看出, 280 s以后速度下降至4 000 m/s以下。根据图 3(b)飞行器控制模式中攻角模式的设定可知, 速度下降至4 000 m/s以下会引起攻角的变化, 攻角的变化带来飞行器的机动变化, 由此会造成跟踪模型跟踪误差增大。因此, 280 s以后所有算法的误差都有所增加可能是由于飞行器机动造成的；第二, 随着目标的机动, 目标与雷达距离越来越远, 量测值误差有所增加, 对跟踪模型的修正效果减弱, 引起跟踪误差上升。

位置与速度误差统计平均如表 1所示。

为检验本文算法在不同量测参数下的跟踪性能, 在保持其他各种仿真条件不变的情况下, 设置了如表 2所示的几种不同量测参数。

经过仿真实验分析, 可得到不同滤波器参数下, 本文算法的位置和速度跟踪误差如表 3所示。

从表 3可以看出, 量测误差协方差增大1倍后, 算法的位置估计误差和速度估计误差分别增大61%、48%, 采样周期增大1倍后, 算法的位置估计误差和速度估计误差分别增大83%、66%, 量测误差协方差和采样周期同时增大1倍后, 算法的位置估计误差和速度估计误差分别增大119%、86%。这是因为采样周期增大后, 跟踪会变得“迟钝”, 收敛速度变慢, 导致误差上升; 量测误差协方差增大后, 导致量测值对预测值修正变弱, 收敛速度变慢, 跟踪误差上升。但总体上, 本文算法在不同量测参数下保持了较为稳定的跟踪性能, 能较好适应滤波器的参数变化。

第3.2节所述仿真实验可以表明:

(1) 对于NSHV目标的跟踪, 不同跟踪模型在其机动样式稳定时, 都能保持较小的误差并进行稳定跟踪, 但当NSHV目标机动样式发生改变时, 由于跟踪模型失配, 不能反映目标真实运动状态, CA+CV+CS的IMM滤波算法、文献[12]所提的跟踪算法、基于本文运动模型的常规卡尔曼滤波算法都会发生明显的误差上升, 而本文算法则可以根据滤波发散情况对渐消因子进行实时调整, 利用最新量测值对预测值进行修正, 所以跟踪效果要优于常规卡尔曼滤波, 且滤波收敛后依然能保持稳定的滤波性能。

(2) 由于跟踪性能的优劣与量测参数直接相关, 从仿真实验2可以看出, 本文算法在不同的量测参数下, 滤波性能有所差异, 但总体还是保持了较为稳定的跟踪, 没有发生明显的滤波发散, 说明本文算法具有一定的“抗干扰能力”, 稳定性较好。

针对NSHV目标跳跃滑翔的机动特性, 本文将其三维运动轨迹分解到纵向弹道和横向弹道, 在两个弹道分别建立适合其运动特征的跟踪模型, 并针对常规卡尔曼滤波在跟踪误差较大时易发生精度下降的现象, 构建了多重渐消因子的自适应卡尔曼滤波, 提高了收敛速度和跟踪精度, 并与常用的几类跟踪模型进行了仿真对比, 结果表明本文算法在跟踪NSHV目标时精度和稳定性均较高。