2. 中国运载火箭技术研究院, 北京 100076;
3. 中国航天科技集团有限公司, 北京 100048
2. China Academy of Launch Vehicle Technology, Beijing 100076, China;
3. China Aerospace Science and Technology Corporation, Beijing 100048, China
与传统弹道式目标不同，高超声速滑翔目标(Hypersonic Glide Target, HGT)具有大升阻比的气动外形，在临近空间以大于5的马赫数飞行，依靠气动力实现远距离滑翔和大范围机动[1]。随着HGT相关技术日趋成熟，对该类目标的跟踪引起了广泛的关注，主要的难点在于较低的飞行高度降低了雷达可探测的时间、所建立的跟踪模型难以囊括其复杂多变的机动模式[2-3]。
与跟踪其他机动目标相比，跟踪HGT存在一定的特殊性但也存在普遍性，相同的条件下量测传感器精度完全由所建立的目标运动模型和采用的滤波算法决定，因此如何实现高精度的HGT跟踪也是主要涉及这两方面。目标运动模型用于表征目标的真实运动，特别是机动目标，滤波算法基于目标运动模型从噪声量测中提取出目标的运动状态。常用的运动模型包括常速度(Constant Velovity, CV)模型(属于非机动模型)常加速度(Constant Acceleration, CA)模型、Singer模型、当前统计(Current Statistics, CS)模型、Jerk模型等，除了CV其他均属于机动模型[4]。由于HGT依靠气动力实施机动，因此气动力加速度对于跟踪系统来说是主要的未知输入项，对其建模的优劣一定程度上决定了目标运动模型的匹配性[5]。文献[2]将HGT气动力加速度中的未知气动系数项建模成维纳随机过程，并扩展至系统状态对其进行滤波估计，实现了稳定跟踪。文献[6]中采用与文献[2]中类似的建模方法，分别对HGT非跳跃式和跳跃式2种飞行轨迹进行了跟踪滤波，发现当非跳跃式轨迹中的倾侧角发生符号翻转时由于目标运动模型匹配不及时会引起较大的估计误差。文献[3, 7]进一步考虑了各气动系数之间的耦合，给出了更加精确的气动力加速度模型，但是需要一定先验控制量的假设。上述模型均是从动力学的角度出发，对气动力加速度进行精细建模，称之为“动力学”模型。当然，从运动学的角度也可直接将HGT的加速度或气动力加速度建模成已有的机动模型，但是与“动力学”模型相比通常缺乏理论根据性和清晰的物理意义，跟踪精度对模型的参数较为敏感，需设置合适[6, 8]。
由于HGT机动模式多样，利用单一固定模型的适应性较差，难以获得高精度的跟踪，甚至针对HGT的某种机动形式现有的机动模型均不能较好地匹配，例如HGT为了调整侧向运动会发生倾侧角符号的切换[9]。文献[10]基于非零均值正弦波相关模型对目标加速度进行均值补偿和方差自适应，提高了跟踪高超声速跳跃滑翔式目标时的模型适应性。文献[11]在认为控制变量服从一阶时滞过程的前提下，实现了气动参数模型与飞行状态的自适应匹配，仿真结果表明当HGT目标发生机动时所提模型的跟踪性能明显优于传统模型。文献[6]利用具有不同机动频率的CS模型构成的交互式多模型(Interacting Multiple Model, IMM)对HGT进行跟踪，并对目标气动力加速度方差进行自适应调整，估计精度优于传统的CS模型。文献[12]设计了一种自适应变结构的IMM算法用于跟踪HGT，提高了跟踪精度和效率。上述基于自适应模型的跟踪方法通常需要一定的假设条件，其普适性和准确性有待商榷，难以实用；而多模型方法中由于缺乏可靠先验信息使模型集个数和模型转移概率难以确定，模型集过少不能全面地表征目标的运动，模型集冗余会增加计算量、引起模型间的竞争，跟踪精度甚至得不到提升。此外，在跟踪HGT的仿真验证时，被跟踪轨迹的设计大都未结合HGT实际可能采用的轨迹设计方法且未考虑完成制导任务的需求，致使跟踪结果的可信度有所欠缺。
如前所述，目标跟踪精度和滤波算法也紧密相关，相同目标运动模型采用不同的滤波算法通常会产生不同的估计精度，在目标跟踪中自适应或鲁棒的滤波方法已经广泛地用于处理模型误差。针对系统模型中未知的噪声统计特性，提出了多种噪声估计器，其中基于极大后验准则的Sage-Husa估计器由于结构简单、原理清晰得到了大量关注[13-15]。但是，文献[16]指出Sage-Husa只有当过程噪声或量测噪声已知其中之一才是有效的。而且，引入估计器后易引起滤波的不收敛[17]。对于状态模型误差，基于正交性原理的强跟踪滤波(Strong Tracking Filtering, STF)表现出较强的鲁棒性，其通过渐消因子放大预测状态误差的协方差，强迫模型误差出现时正交性原理依然能够尽可能满足，进而调整增益矩阵来增加新观测数据的比重，以降低状态估计误差[18]。文献[19]将STF用于跟踪脉冲机动卫星，可快速地降低脉冲机动引起的估计误差，实现了高精度稳定跟踪。文献[20]将STF应用至高阶容积卡尔曼滤波(Kalman Filtering, KF)，增强了其跟踪目标时的鲁棒性。文献[21]提出了一种归一化的STF用于跟踪小脉冲机动的航天器。由此，针对HGT运动模型误差几乎无法避免的情况，也可配合鲁棒的滤波算法来提高跟踪精度。
众所周知，系统噪声的统计特性对卡尔曼滤波结果的影响较大。上述所提跟踪方法大多通过不同方式对系统过程噪声的协方差实施自适应调整，最终提高跟踪滤波精度。文献[22]在跟踪机动飞机时通过CS模型中的目标加速度递归方程和均值输入估计方法提出了一种自适应CS模型，实现了对过程噪声协方差的自适应调整。为解决跟踪机动目标时Singer模型参数的失配问题，文献[23]采用多模型方法中的似然函数并根据其变化自适应地调整最大加速度参数。对机动模型中的关键参数实施自适应调整的方法结构简单、原理清晰，对跟踪HGT具有一定借鉴意义。
为了提高HGT跟踪方法的适应性和精度，本文在现有研究的基础上提出了一种机动频率自适应的跟踪方法。首先将气动力加速度建模成形式简洁的Singer模型，基于正交性原理和无迹卡尔曼滤波算法利用滤波新息计算得到可反映状态模型误差大小的调整因子，然后用于放大Singer模型中的机动频率参数，实现对当前模型误差的自适应。仿真结果表明所提方法的适应性好、计算量小，能够对HGT 2种典型的机动飞行轨迹实现高精度跟踪。
1 目标跟踪模型目标跟踪模型通常由状态方程和量测方程组成，所建模型与目标真实运动模式的匹配程度直接影响着跟踪精度。
1.1 状态方程图 1给出了HGT相对于地基雷达再入运动的示意图，oxyz表示雷达坐标系，也即东北天坐标系。考虑地球旋转，则HGT的相对运动方程可表示为
|
| 图 1 目标相对地基雷达的再入运动 Fig. 1 Target reentry motion relative to ground-based radar |
| $ \frac{{{\delta ^2}\mathit{\boldsymbol{r}}}}{{\delta {t^2}}} = \mathit{\boldsymbol{g}} + \mathit{\boldsymbol{a}} - 2\mathit{\boldsymbol{\omega }} \times \frac{{\delta \mathit{\boldsymbol{r}}}}{{\delta t}} - \mathit{\boldsymbol{\omega }} \times (\mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{r}}) $ | （1） |
式中：δr/δt和δ2r/δt2分别为相对速度v和相对加速度；g和a分别为目标所受的引力加速度和气动力加速度；ω和r分别为地球自转角速度和目标地心矢径；2ω×(δr/δt)和ω×(ω×r)分别为科式加速度项和牵连加速度项。
为了获得系统的状态方程，需将式(1)投影至雷达坐标系下。设v=[vx, vy, vz]T，则相对加速度的各分量为
| $ \frac{{{\delta ^2}\mathit{\boldsymbol{r}}}}{{\delta {t^2}}} = \left[ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}{v_x}}}{{{\rm{d}}t}}}\\ {\frac{{{\rm{d}}{v_y}}}{{{\rm{d}}t}}}\\ {\frac{{{\rm{d}}{v_z}}}{{{\rm{d}}t}}} \end{array}} \right] $ | （2） |
考虑椭圆地球，则引力加速度的形式为[24]
| $ \mathit{\boldsymbol{g}} = {g_r}{\mathit{\boldsymbol{r}}^0} + {g_\omega }{\mathit{\boldsymbol{\omega }}^0} $ | （3） |
式中：r0和ω0分别为r和ω方向上的单位矢量，其各分量分别为
| $ {\mathit{\boldsymbol{r}}^0} = \frac{1}{r}\left[ {\begin{array}{*{20}{c}} {{p_x}}\\ {{p_y} - {R_{\rm{o}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}}\\ {{p_z} + {R_{\rm{o}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}} \end{array}} \right], {\mathit{\boldsymbol{\omega }}^0} = \left[ {\begin{array}{*{20}{c}} 0\\ {{\rm{cos}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}}}\\ {{\rm{sin}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}}} \end{array}} \right] $ | （4） |
其中：r为目标地心矢径r的模; px、py和pz为p的3个分量; Ro为Ro的模；Bo为地理纬度; μo=Bo－ϕo，ϕo为地心纬度。gr和gω的表达式分别为
| $ \left\{ {\begin{array}{*{20}{l}} {{g_r} = - \frac{\mu }{{{r^2}}}\left[ {1 + \frac{3}{2}{J_2}{{\left( {\frac{{{a_{\rm{E}}}}}{r}} \right)}^2}(1 - 5{\rm{si}}{{\rm{n}}^2}{\phi _{\rm{o}}})} \right]}\\ {{g_\omega } = - 3\frac{\mu }{{{r^2}}}{J_2}{{\left( {\frac{{{a_{\rm{E}}}}}{r}} \right)}^2}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\phi _{\rm{o}}}} \end{array}} \right. $ | （5） |
式中：μ为地球引力常数; J2为带谐修正系数; aE为地球的长半轴。式(5)中，为了计算简便近似地用地基雷达纬度ϕo替代目标纬度ϕ。
科式和牵连加速度项的表达式分别为
| $ 2\mathit{\boldsymbol{\omega }} \times \frac{{\delta \mathit{\boldsymbol{r}}}}{{\delta t}} = 2\omega \left[ {\begin{array}{*{20}{c}} { - {v_y}{\rm{sin}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}} + {v_z}{\rm{cos}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}}}\\ {{v_x}{\rm{sin}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}}}\\ { - {v_x}{\rm{cos}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}}} \end{array}} \right] $ | （6） |
| $ \mathit{\boldsymbol{\omega }} \times (\mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{r}}) = {\omega ^2}\left[ {\begin{array}{*{20}{c}} { - {r_x}}\\ {{r_z}{\rm{sin}}{\kern 1pt} {B_{\rm{o}}}{\rm{cos}}{\kern 1pt} {B_{\rm{o}}} - {r_y}{\rm{si}}{{\rm{n}}^2}{B_{\rm{o}}}}\\ {{r_y}{\rm{sin}}{\kern 1pt} {B_{\rm{o}}}{\rm{cos}}{\kern 1pt} {B_{\rm{o}}} - {r_z}{\rm{co}}{{\rm{s}}^2}{B_{\rm{o}}}} \end{array}} \right] $ | （7） |
式中：r=[rx, ry, rz]T; ω为ω的模。
气动力加速度作为HGT跟踪系统中主要的未知输入项，需进行重点考虑。这里直接从运动学的角度，将a建模成现有的机动模型。考虑HGT一直受到气动力的作用，且机动强度有所起伏，选用介于CV模型和CA模型之间的Singer模型来匹配a的真实变化，因此有[4, 25]
| $ \mathit{\boldsymbol{a}} = \left[ {\begin{array}{*{20}{l}} {{a_x}}\\ {{a_y}}\\ {{a_z}} \end{array}} \right], \left\{ {\begin{array}{*{20}{l}} {{{\dot a}_x} = - {\lambda _x}{a_x} + {w_x}}\\ {{{\dot a}_y} = - {\lambda _y}{a_y} + {w_y}}\\ {{{\dot a}_z} = - {\lambda _z}{a_z} + {w_z}} \end{array}} \right. $ | （8） |
式中：λx、λy和λz分别为3个方向上的机动频率；wx、wy和wz为互不相关的零均值高斯白噪声，功率谱密度分别为2λxσ2、2λyσ2和2λzσ2，σ2为当前时刻目标加速度的方差，表达式为
| $ {\sigma ^2} = \frac{{a_{{\rm{max}}}^2}}{3}(1 + 4{P_{{\rm{max}}}} - {P_0}) $ | （9） |
式中：amax为目标加速度的最大幅值，目标以amax或－amax实施机动的概率均为Pmax，P0为目标加速度为零的概率。
针对上述运动学建模，系统状态方程可写
| $ \mathit{\boldsymbol{\dot X}} = \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{X}}) + \mathit{\boldsymbol{w}} $ | （10） |
式中：X为系统状态向量，具体表达式为
| $ \mathit{\boldsymbol{X}} = {[{p_x}, {p_y}, {p_z}, {v_x}, {v_y}, {v_z}, {a_x}, {a_y}, {a_z}]^{\rm{T}}} $ | （11） |
F(X)的表达式见附录A; w表示过程噪声，为互不相关的零均值白噪声向量，其自相关函数为
| $ E[\mathit{\boldsymbol{w}}({t_1}){\mathit{\boldsymbol{w}}^{\rm{T}}}({t_2})] = \mathit{\boldsymbol{W}}\delta ({t_1} - {t_2}) $ | （12） |
式中：W为对角矩阵; δ(·)为狄拉克δ函数。
为了后续的滤波估计，需将式(10)进行离散化。设采样周期为T，则式(10)转化为
| $ {\mathit{\boldsymbol{X}}_k} = {\mathit{\boldsymbol{X}}_{k - 1}} + \int_{(k - 1)T}^{kT} {\mathit{\boldsymbol{F}}(\mathit{\boldsymbol{X}})} {\rm{d}}t + {{\mathit{\boldsymbol{\varpi }}} _{k - 1}} $ | （13） |
式中：下标k表示变量在离散kT时刻的值；
| $ \begin{array}{*{20}{l}} {E[{\mathit{\boldsymbol{\varpi }}} (k){{\mathit{\boldsymbol{\varpi }}} ^{\rm{T}}}(k)] = }\\ {\int_{kT}^{(k + 1)T} {\int_{kT}^{(k + 1)T} E } [\mathit{\boldsymbol{w}}(t){\mathit{\boldsymbol{w}}^{\rm{T}}}(\tau )]{\rm{d}}t{\rm{d}}\tau = }\\ {\quad \int_{kT}^{(k + 1)T} {\int_{kT}^{(k + 1)T} \mathit{\boldsymbol{W}} } \delta (t - \tau ){\rm{d}}t{\rm{d}}\tau = \mathit{\boldsymbol{W}}T} \end{array} $ | （14） |
地基雷达提供的量测量包括相对距离R、方位角A和高低角E，如图 2所示。
|
| 图 2 地基雷达量测量 Fig. 2 Ground-based rabar measurements |
由图 2，易得量测方程为
| $ \mathit{\boldsymbol{Z}}_k^s = \mathit{\boldsymbol{G}}({\mathit{\boldsymbol{p}}_k}) + \mathit{\boldsymbol{\varepsilon }}_k^s $ | （15） |
式中：Zks为含有噪声的量测量；εsk为零均值高斯白噪声向量序列，并与
| $ \mathit{\boldsymbol{G}}(\mathit{\boldsymbol{p}}) = \left[ {\begin{array}{*{20}{c}} {\sqrt {p_x^2 + p_y^2 + p_z^2} }\\ {{\rm{arctan}}({p_y}/{p_x})}\\ {{\rm{arctan}}({p_z}/\sqrt {p_x^2 + p_y^2} )} \end{array}} \right] $ | （16） |
式中：p=[px，py，pz]T。
由于相对距离和角度之间的量级相差较大，为了降低滤波过程中出现矩阵条件数病态的可能性，将量测量转换至直角坐标系下，可表示为
| $ \mathit{\boldsymbol{Z}}_k^{\rm{r}} = {\mathit{\boldsymbol{p}}_k} + \mathit{\boldsymbol{\varepsilon }}_k^{\rm{r}} = {\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{X}}_k} + \mathit{\boldsymbol{\varepsilon }}_k^{\rm{r}} $ | （17） |
式中:Hk为系数矩阵。
对pk在Zks处进行泰勒展开有
| $ {\mathit{\boldsymbol{p}}_k} = {\mathit{\boldsymbol{G}}^{ - 1}}(\mathit{\boldsymbol{Z}}_k^{\rm{s}} - \mathit{\boldsymbol{\varepsilon }}_k^{\rm{s}}) = \mathit{\boldsymbol{Z}}_k^{\rm{r}} - \mathit{\boldsymbol{J}}(\mathit{\boldsymbol{Z}}_k^{\rm{s}})\mathit{\boldsymbol{\varepsilon }}_k^{\rm{s}} + \mathit{\boldsymbol{O}}(\mathit{\boldsymbol{\varepsilon }}_k^{\rm{s}}) $ | （18） |
式中：G－1(·)为G(·)的反函数；O(εks)为高阶项；J(Zks)为雅克比矩阵
| $ \mathit{\boldsymbol{\varepsilon }}_k^{\rm{r}} \approx \mathit{\boldsymbol{J}}(\mathit{\boldsymbol{Z}}_k^{\rm{s}})\mathit{\boldsymbol{\varepsilon }}_k^{\rm{s}} $ | （19） |
因此，可近似地认为εkr的均值为零，协方差为
| $ \mathit{\boldsymbol{R}}_k^{\rm{r}} = {\rm{cov}} (\mathit{\boldsymbol{\varepsilon }}_k^{\rm{r}}) \approx \mathit{\boldsymbol{J}}(\mathit{\boldsymbol{Z}}_k^{\rm{s}})\mathit{\boldsymbol{R}}_k^{\rm{s}}\mathit{\boldsymbol{J}}{(\mathit{\boldsymbol{Z}}_k^{\rm{s}})^{\rm{T}}} $ | （20） |
式中：
鉴于HGT机动模式多样，仅用单一固定的Singer模型难以表征气动力加速度的变化。在卡尔曼滤波过程中，新息反映着系统模型所存在的误差。假设量测方程准确的前提下，新息代表着状态模型误差，下面将其用于自适应Singer模型中的机动频率参数。起始将式(8)中的机动频率取至很小，此时Singer模型无限接近于CA模型，认为机动加速度变化平缓。当HGT机动加速度变化幅度增加时，放大机动频率参数调整过程噪声协方差。
由于状态方程是非线性的，因此只能选用非线性的滤波算法。相比于扩展卡尔曼滤波(Extended KF, EKF)，基于无迹变换的UKF省去了复杂雅克比矩阵的运算，且具有更高的理论精度[26]。因此，选用UKF作为滤波估计的框架。
2.1 机动频率自适应对于线性KF，若系统模型能够完全匹配真实模型，则滤波器输出的新息序列满足正交性原理，即[26]
| $ E({\mathit{\boldsymbol{\xi }}_{k + j}}\mathit{\boldsymbol{\xi }}_k^{\rm{T}}) = {\bf{0}}, k = 1, 2, \cdots , j = 1, 2, \cdots $ | （21） |
式中：ξk为新息向量，即量测量与其预测值的差。当式(21)满足时，意味着量测量中的有用信息被完全提取出来，因此新息序列的正交性可作为衡量滤波性能的标志。
对于次优的非线性滤波器，新息序列不可能是完全正交的。此时若是弱自相关的，即可认为其具有较好的滤波性能。以EKF为例，考虑如下形式的非线性系统：
| $ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{x}}_k} = {\mathit{\boldsymbol{f}}_{k - 1}}({\mathit{\boldsymbol{x}}_{k - 1}}) + {\mathit{\boldsymbol{q}}_{k - 1}}}\\ {{\mathit{\boldsymbol{y}}_k} = {\mathit{\boldsymbol{g}}_k}({\mathit{\boldsymbol{x}}_k}) + {\mathit{\boldsymbol{r}}_k}} \end{array}} \right. $ | （22） |
式中：xk∈Rn和yk∈Rm分别为系统状态和量测向量；fk(·)和gk(·)分别为状态和量测方程的非线性函数；qk和rk为互不相关的零均值白噪声，其协方差分别为Qk和Rk。设
| $ \begin{array}{l} E({\mathit{\boldsymbol{\xi }}_{k + j}}\mathit{\boldsymbol{\xi }}_k^{\rm{T}}) \approx {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{k + j}}[\prod\limits_{l = k + 1}^{l = k + j - 1} {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_l}} (\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_l}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_l})] \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k}E[({\mathit{\boldsymbol{x}}_k} - {{\mathit{\boldsymbol{\hat x}}}_k})\mathit{\boldsymbol{\xi }}_k^{\rm{T}}] \end{array} $ | （23） |
式中：Kk为滤波增益；Φk－1和Γk的表达式为
| $ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k - 1}} = {\left. {\frac{{\partial {\mathit{\boldsymbol{f}}_{k - 1}}}}{{\partial \mathit{\boldsymbol{x}}}}} \right|_{\mathit{\boldsymbol{x}} = {{\mathit{\boldsymbol{\hat x}}}_{k - 1}}}}, {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k} = {\left. {\frac{{\partial {\mathit{\boldsymbol{g}}_k}}}{{\partial \mathit{\boldsymbol{x}}}}} \right|_{\mathit{\boldsymbol{x}} = {{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}}}} $ | （24） |
因此，对于EKF使新息正交性满足的充分条件为
| $ E[({\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{\hat x}}_k})\mathit{\boldsymbol{\xi }}_k^{\rm{T}}] = {\bf{0}} $ | （25） |
由于
| $ {\mathit{\boldsymbol{\tilde x}}_k} = {\mathit{\boldsymbol{\tilde x}}_{k|k - 1}} - {\mathit{\boldsymbol{K}}_k}{\mathit{\boldsymbol{\xi }}_k} $ | （26） |
式中：
| $ {\mathit{\boldsymbol{y}}_k} = {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{\mathit{\boldsymbol{x}}_k} + {\mathit{\boldsymbol{g}}_k}({\mathit{\boldsymbol{\hat x}}_{k|k - 1}}) - {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{\mathit{\boldsymbol{\hat x}}_{k|k - 1}} + {\mathit{\boldsymbol{r}}_k} $ | （27） |
那么，量测量的预测值为
| $ {\mathit{\boldsymbol{\hat y}}_{k|k - 1}} = {\mathit{\boldsymbol{g}}_k}({\mathit{\boldsymbol{\hat x}}_{k|k - 1}}) $ | （28） |
因此，ξk表示为
| $ {\mathit{\boldsymbol{\xi }}_k} = {\mathit{\boldsymbol{y}}_k} - {\mathit{\boldsymbol{\hat y}}_{k|k - 1}} = {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{\mathit{\boldsymbol{\tilde x}}_{k|k - 1}} + {\mathit{\boldsymbol{r}}_k} $ | （29） |
将式(29)和式(26)代入式(25)得
| $ {\mathit{\boldsymbol{P}}_{k|k - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k^{\rm{T}} + E({\mathit{\boldsymbol{\tilde x}}_{k|k - 1}}\mathit{\boldsymbol{r}}_k^{\rm{T}}) - {\mathit{\boldsymbol{K}}_k}{\mathit{\boldsymbol{V}}_k} = {\bf{0}} $ | （30） |
式中：
| $ {\mathit{\boldsymbol{P}}_{k|k - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k^{\rm{T}} - {\mathit{\boldsymbol{K}}_k}{\mathit{\boldsymbol{V}}_k} = {\bf{0}} $ | （31） |
由于
| $ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{\mathit{\boldsymbol{P}}_{k|k - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k^{\rm{T}} = {\mathit{\boldsymbol{V}}_k} - {\mathit{\boldsymbol{R}}_k} $ | （32） |
当状态模型存在误差时，会导致式(32)不再满足，左右两端不相等的程度反映着误差的大小，并采用左右两端主对线上元素的比值对其进行量化，具体形式为
| $ {\zeta _d} = \frac{{{{({\mathit{\boldsymbol{V}}_k} - {\mathit{\boldsymbol{R}}_k})}_{dd}}}}{{{{({\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{\mathit{\boldsymbol{P}}_{k|k - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k^{\rm{T}})}_{dd}}}}\quad d = 1, 2, \cdots , m $ | （33） |
式中：ζd称为调整因子，反映了第d个量测通道上误差大小。当目标进行大幅度机动引起状态模型误差时，可利用ζd对机动模型中的机动频率实施相应调整，以降低模型误差，提高滤波精度。当存在模型误差时，Vk的真实值在实际计算中是未知的，可以通过式(34)对其进行估算：
| $ {\mathit{\boldsymbol{V}}_k} = \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{\xi }}_1}\mathit{\boldsymbol{\xi }}_1^{\rm{T}}}&{k = 1}\\ {\frac{{\kappa {\mathit{\boldsymbol{V}}_{k - 1}} + {\mathit{\boldsymbol{\xi }}_k}\mathit{\boldsymbol{\xi }}_k^{\rm{T}}}}{{1 + \kappa }}}&{k > 1} \end{array}} \right. $ | （34） |
式中：κ为遗忘因子，通常设置为0.95。
针对上述所建立的HGT跟踪模型，量测量转换后三通道上的调整因子可分别用于调整x、y和z方向上的机动频率，x方向的调整策略为
| $ {\hat \lambda _x} = \left\{ {\begin{array}{*{20}{l}} {{\zeta _x}{\lambda _x}}&{{\zeta _x} \ge {\zeta _{{\rm{max}}}}}\\ {{\lambda _x}}&{{\zeta _x} < {\zeta _{{\rm{max}}}}} \end{array}} \right. $ | （35） |
式中：ζx为雷达坐标系x方向上的调整因子；ζmax>0为调整阀值；
由于式(33)是基于EKF推导得出的，不能直接嵌入至UKF框架下，因此采用如式(36)所示的等效计算式：
| $ {\zeta _d} = \frac{{{{({\mathit{\boldsymbol{V}}_k} - {\mathit{\boldsymbol{R}}_k})}_{dd}}}}{{{{({\mathit{\boldsymbol{P}}_y} - {\mathit{\boldsymbol{R}}_k})}_{dd}}}} $ | （36） |
式中：Py为依赖式(22)所描述的模型计算得到的新息协方差矩阵的标称值，表达式为
| $ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{P}}_y} = E[({\mathit{\boldsymbol{y}}_k} - {{\mathit{\boldsymbol{\hat y}}}_{k|k - 1}}){{({\mathit{\boldsymbol{y}}_k} - {{\mathit{\boldsymbol{\hat y}}}_{k|k - 1}})}^{\rm{T}}}] = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} E[({\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{{\mathit{\boldsymbol{\tilde x}}}_{k|k - 1}} + {\mathit{\boldsymbol{r}}_k}){{({\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{{\mathit{\boldsymbol{\tilde x}}}_{k|k - 1}} + {\mathit{\boldsymbol{r}}_k})}^{\rm{T}}}] = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{\mathit{\boldsymbol{P}}_{k|k - 1}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k^{\rm{T}} + {\mathit{\boldsymbol{R}}_k}} \end{array} $ | （37） |
对于式(22)所描述的非线性系统，首先对滤波器进行如式(38)的初始化：
| $ {\mathit{\boldsymbol{\hat x}}_0} = E({\mathit{\boldsymbol{x}}_0}), {\mathit{\boldsymbol{P}}_0} = E[({\mathit{\boldsymbol{x}}_0} - {\mathit{\boldsymbol{\hat x}}_0}){({\mathit{\boldsymbol{x}}_0} - {\mathit{\boldsymbol{\hat x}}_0})^{\rm{T}}}] $ | （38） |
自适应UKF算法可分为如下几个步骤[26]：
1) 计算simga点和权重系数
| $ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\hat x}}_{k - 1}^{(0)} = {{\mathit{\boldsymbol{\hat x}}}_{k - 1}}}&{}\\ {\mathit{\boldsymbol{\hat x}}_{k - 1}^{(i)} = {{\mathit{\boldsymbol{\hat x}}}_{k - 1}} + \Delta \mathit{\boldsymbol{\hat x}}_{k - 1}^{(i)}}&{i = 1, 2, \cdots , 2n}\\ {\Delta \mathit{\boldsymbol{\hat x}}_{k - 1}^{(i)} = \sqrt {n + \chi } (\sqrt {{\mathit{\boldsymbol{P}}_{k - 1}}} )_i^{\rm{T}}}&{i = 1, 2, \cdots , n}\\ {\Delta \mathit{\boldsymbol{\hat x}}_{k - 1}^{(i + n)} = - \sqrt {n + \chi } (\sqrt {{\mathit{\boldsymbol{P}}_{k - 1}}} )_i^{\rm{T}}}&{i = 1, 2, \cdots , n} \end{array}} \right. $ | （39） |
式中：
| $ \left\{ {\begin{array}{*{20}{l}} {w_{\rm{m}}^{(0)} = \frac{\chi }{{n + \chi }}}\\ {w_{\rm{c}}^{(0)} = \frac{\chi }{{n + \chi }} + (1 - {\varphi ^2} + \beta )}\\ {w_{\rm{m}}^{(i)} = w_{\rm{c}}^{(i)} = \frac{1}{{2(n + \chi )}}\quad i = 1, 2, \cdots , 2n} \end{array}} \right. $ | （40） |
式中：β为与系统状态先验分布有关的常数，通常设为2。
2) 时间更新
| $ {\mathit{\boldsymbol{\hat x}}_{i, k|k - 1}} = {\mathit{\boldsymbol{f}}_{k - 1}}(\mathit{\boldsymbol{\hat x}}_{k - 1}^{(i)}), {\mathit{\boldsymbol{\hat x}}_{k|k - 1}} = \sum\limits_{i = 0}^{2n} {w_{\rm{m}}^{(i)}} {\mathit{\boldsymbol{\hat x}}_{i, k|k - 1}} $ | （41） |
| $ \begin{array}{l} {\mathit{\boldsymbol{P}}_{k|k - 1}} = \sum\limits_{i = 0}^{2n} {w_{\rm{c}}^{(i)}} ({{\mathit{\boldsymbol{\hat x}}}_{i, k|k - 1}} - {{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}})({{\mathit{\boldsymbol{\hat x}}}_{i, k|k - 1}} - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}}{)^{\rm{T}}} + {\mathit{\boldsymbol{Q}}_{k - 1}} \end{array} $ | （42） |
3) 机动频率自适应
| $ {\mathit{\boldsymbol{\hat y}}_{i, k|k - 1}} = {\mathit{\boldsymbol{g}}_k}(\mathit{\boldsymbol{\hat x}}_{k|k - 1}^{(i)}), {\mathit{\boldsymbol{\hat y}}_{k|k - 1}} = \sum\limits_{i = 0}^{2n} {w_{\rm{m}}^{(i)}} {\mathit{\boldsymbol{\hat y}}_{i, k|k - 1}} $ | （43） |
式中：
| $ {\mathit{\boldsymbol{\xi }}_k} = {\mathit{\boldsymbol{y}}_k} - {\mathit{\boldsymbol{\hat y}}_{k|k - 1}}, {\mathit{\boldsymbol{V}}_k} = \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{\xi }}_1}\mathit{\boldsymbol{\xi }}_1^{\rm{T}}}&{k = 1}\\ {\frac{{\kappa {\mathit{\boldsymbol{V}}_{k - 1}} + {\mathit{\boldsymbol{\xi }}_k}\mathit{\boldsymbol{\xi }}_k^{\rm{T}}}}{{1 + \kappa }}}&{k > 1} \end{array}} \right. $ | （44） |
| $ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{P}}_y} = \sum\limits_{i = 0}^{2n} {w_{\rm{c}}^{(i)}} ({{\mathit{\boldsymbol{\hat y}}}_{i, k|k - 1}} - {{\mathit{\boldsymbol{\hat y}}}_{k|k - 1}})({{\mathit{\boldsymbol{\hat y}}}_{i, k|k - 1}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\mathit{\boldsymbol{\hat y}}}_{k|k - 1}}{)^{\rm{T}}} + {\mathit{\boldsymbol{R}}_k}} \end{array} $ | （45） |
| $ {{\zeta _d} = \frac{{{{({\mathit{\boldsymbol{V}}_k} - {\mathit{\boldsymbol{R}}_k})}_{dd}}}}{{{{({\mathit{\boldsymbol{P}}_y} - {\mathit{\boldsymbol{R}}_k})}_{dd}}}}, d = 1, 2, \cdots , m} $ | （46） |
| $ {{{\mathit{\boldsymbol{\bar P}}}_{k|k - 1}} = {\mathit{\boldsymbol{P}}_{k|k - 1}} - {\mathit{\boldsymbol{Q}}_{k - 1}} + {{\mathit{\boldsymbol{\bar Q}}}_{k - 1}}} $ | （47） |
式中：
4) 量测更新
| $ \mathit{\boldsymbol{\hat y}}_{i, k|k - 1}^\prime = {\mathit{\boldsymbol{g}}_k}(\mathit{\boldsymbol{\hat x}}_{k|k - 1}^i), \mathit{\boldsymbol{\hat y}}_{k|k - 1}^\prime = \sum\limits_{i = 0}^{2n} {w_{\rm{m}}^{(i)}} \mathit{\boldsymbol{\hat y}}_{i, k|k - 1}^\prime $ | （48） |
| $ \begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\bar P}}}_y} = \sum\limits_{i = 0}^{2n} {w_{\rm{c}}^{(i)}} (\mathit{\boldsymbol{\hat y}}_{i, k|k - 1}^\prime - \mathit{\boldsymbol{\hat y}}_{k|k - 1}^\prime )(\mathit{\boldsymbol{\hat y}}_{i, k|k - 1}^\prime - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{\hat y}}_{k|k - 1}^\prime {)^{\rm{T}}} + {\mathit{\boldsymbol{R}}_k}} \end{array} $ | （49） |
| $ {\mathit{\boldsymbol{\bar P}}_{xy}} = \sum\limits_{i = 0}^{2n} {w_c^{(i)}} (\mathit{\boldsymbol{\hat x}}_{k|k - 1}^i - {\mathit{\boldsymbol{\hat x}}_{k|k - 1}}){(\mathit{\boldsymbol{\hat y}}_{i, k|k - 1}^\prime - \mathit{\boldsymbol{\hat y}}_{k|k - 1}^\prime )^{\rm{T}}} $ | （50） |
| $ {{{\mathit{\boldsymbol{\hat x}}}_k} = {{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}} + {\mathit{\boldsymbol{K}}_k}({\mathit{\boldsymbol{y}}_k} - {{\mathit{\boldsymbol{\hat y}}}_{k|k - 1}}), {\mathit{\boldsymbol{K}}_k} = {{\mathit{\boldsymbol{\bar P}}}_{xy}}\mathit{\boldsymbol{\bar P}}_\mathit{\boldsymbol{y}}^{ - 1}} $ | （51） |
| $ {{\mathit{\boldsymbol{P}}_k} = {{\mathit{\boldsymbol{\bar P}}}_{k|k - 1}} - {\mathit{\boldsymbol{K}}_k}{{\mathit{\boldsymbol{\bar P}}}_y}\mathit{\boldsymbol{K}}_k^{\rm{T}}} $ | （52） |
式中：
为了验证所提方法跟踪滤波的性能，首先设计2条典型HGT的飞行轨迹用于跟踪，并与其他几种跟踪模型进行对比。
3.1 HGT轨迹设计以美国洛马公司发布的通用飞行器CAV-H模型为研究对象，分别基于阻力加速度-速度剖面跟踪和高斯伪谱优化方法生成2条再入飞行轨迹。考虑飞行过程中受多种约束的影响，主要设置最大过载为3g、最大动压为100 kPa和最大驻点热流密度为1 800 kW/m2，控制量攻角和倾侧角的最大变化率分别为20 (°)/s和85 (°)/s，初终端的状态约束如表 1所示[12]。
跟踪阻力加速度-速度剖面是HGT再入段轨迹规划与制导中常用的一种方法，其将再入段分为初始下降段和滑翔段，以是否满足平衡滑翔条件为交班点，攻角剖面事先给定，初始下降段采用零倾侧角飞行，在滑翔段跟踪设计的阻力加速度-速度剖面用于控制纵向运动，由此确定倾侧角的幅值，侧向运动由视线方位角误差走廊进行控制，用于确定倾侧角符号的翻转时机。以纵向跳跃幅度最小为性能指标，基于高斯伪谱方法进行轨迹优化。具体两轨迹的形态和状态量的变化如图 3~图 5所示，“轨迹1”和“轨迹2”分别为基于跟踪剖面和高斯伪谱方法获得的轨迹。可看出，“轨迹1”存在倾侧角符号的翻转，会引起HGT的阶跃机动，当倾侧角不翻转时机动较平缓；“轨迹2”的倾侧角连续变化，机动幅度强于“轨迹1”中倾侧角符号不翻转时的幅度。上述2种典型轨迹可充分验证跟踪模型对不同机动模式的适应性。
|
| 图 3 HGT跟踪轨迹 Fig. 3 HGT trajectories to be tracked |
|
| 图 4 控制量随时间变化 Fig. 4 Control variables versus time |
|
| 图 5 升阻加速度随时间变化 Fig. 5 Lift and drag accelerations versus time |
为了充分地验证所提机动频率自适应方法的有效性，将其分别与基于Singer模型的方法、文献[6]中提出的基于CS模型的方法、文献[7]中提出的基于动力学模型的方法和文献[6]中提出的利用IMM实现机动频率自适应的方法进行性能对比。基于Singer模型的方法将气动力加速度建模成Singer模型，如式(8)所示，基于动力学模型的方法见附录B，基于CS模型的方法将气动力加速度建模成CS模型，形式为[6]
| $ \left\{ {\begin{array}{*{20}{l}} {{{\dot a}_x} = - {\lambda _x}{a_x} + {\lambda _x}{{\bar a}_x} + {w_x}}\\ {{{\dot a}_y} = - {\lambda _y}{a_y} + {\lambda _y}{{\bar a}_y} + {w_y}}\\ {{{\dot a}_z} = - {\lambda _z}{a_z} + {\lambda _z}{{\bar a}_z} + {w_z}} \end{array}} \right. $ | （53） |
式中：ax、ay和a z为当前时刻目标加速度的均值。与Singer模型相比，CS模型认为加速度均值不为零，方差由修正的瑞利分布确定。由于实际中目标当前加速度的均值无法获知，通常将上一步加速度的滤波估计值当作当前的均值，通过均值实现噪声方差的自适应，详细细节可见文献[27]。文献[6]中提出的基于IMM方法[28]的框图如图 6所示，其给出了三模型交互式滤波器的结构，主要包括输入交互、滤波计算、模型概率更新和输出交互，三模型采用不同机动频率的Singer模型，对机动频率进行自适应。
|
| 图 6 交互式多模型算法结构 Fig. 6 Structure of interacting multiple model algorithm |
为了后续书写方便，将基于Singer模型的方法简记为“Singer”，在“Singer”基础上引入自适应机动频率的方法记为“Adaptive Singer”，基于动力学模型的方法记为“Dynamics”，基于CS模型的方法记为“Current Statistics”，图 6所示的基于IMM的方法记为“IMM-Singer”。为了保持仿真验证的一致性，除“Adaptive Singer”中基于UKF算法对机动频率进行自适应调整外，其他方法也均采用UKF算法。每种方法进行300次的蒙特卡罗仿真，采用均方根误差(Root Mean Square Error, RMSE)来评估滤波结果的估计精度，位置RMSE的计算公式为
| $ {\rm{RMS}}{{\rm{E}}_{\rm{p}}}(k) = \sqrt {\frac{1}{M}\sum\limits_{h = 1}^M {\left\| {{\mathit{\boldsymbol{p}}_k} - {{\mathit{\boldsymbol{\hat p}}}_k}} \right\|_2^h} } $ | （54） |
式中：M为蒙特卡罗仿真次数；
将雷达部署在经纬高为(35°, 0.5°, 0)的位置，设其最大量测距离为800 km。考虑地球曲率影响，易得雷达可量测的轨迹弧段及相应状态量的变化，如图 7~图 9所示。由于图 4中t≈600 s时，倾侧角发生了符号翻转，相应地在图 9中引起了雷达坐标系下气动力加速度的突变。
|
| 图 7 可观测弧段(轨迹1) Fig. 7 Observable segment (trajectory 1) |
|
| 图 8 量测量真值(轨迹1) Fig. 8 Noise-free measurements (trajectory 1) |
|
| 图 9 Fig. 9 True aerodynamic accelerations in radar frame (trajectory 1) |
目标初始位置和速度的真实值分别为(-789.329 0, 57.556 4, 0.002 9) km和(5 430.414 0, -365.939 3, 668.137 8) m/s。“Dynamics”方法中，根据各状态估计可能的误差量级将初始状态估计误差的协方差设为P0=diag{1002, 1002, 1002, 502, 502, 502, 10-8, 10-8, 10-8}，其他模型均属于运动学模型，P0=diag{1002, 1002, 1002, 502, 502, 502, 102, 102, 102}。每次仿真初始位置和速度的估计值为真实值加随机误差，位置分量和速度分量的随机误差分别从正态分布Ν(0, σp2)和Ν(0, σv2)中抽取，气动力加速度和广义气动系数的初始估计值均设为零。根据P0，设σp100 m，σv=50 m/s。雷达量测噪声的协方差为Rsk=diag{(10 m)2, (10-3 rad)2, (10-3 rad)2}，采样周期T=0.05 s。UKF算法中，设φ=0.1，ψ=－6，β=2。
“Dynamics”方法中，与过程噪声协方差矩阵相关的对角矩阵W设为diag{10-10, 10-10, 10-10, 10-10, 10-10, 10-10, 10-12, 10-12, 10-12}。W对角线上前6个元素对应的状态方程部分建模较准确，可认为几乎不存在过程噪声，因此相对于位置和速度的量级取值很小；后3个元素对应着未知广义气动系数模型中的噪声强度，取值大小影响着模型精确程度，这里为了与“Adaptive Singer”作比较通过多次试探仿真选取了10-12，此时估计结果相对较好。由于广义气动系数本身量级较小，所以后3个元素取值也随即较小。其他方法的对角矩阵中，对角线上的前6个元素与“Dynamics”方法取值一致，后3个元素对应着气动力加速度模型中的噪声强度，可由所采用的机动模型参数确定。“Singer”方法中，amax=2g0，g0为海平面的重力加速度，取9.806 65 m/s2，Pmax=0.05，P0=0，3个方向上的机动频率取为λx=λy=λz=1/100 000，接近于CA模型。“Adaptive Singer”方法在“Singer”的基础上利用所提方法对机动频率进行自适应，调整阀值ζmax=100。“Current Statistics”方法中，机动频率和最大加速度与Singer模型中一致，当前时刻目标加速度的均值采用上一步滤波的估计值。“IMM-Singer”中，模型1、2和3的机动频率分别为
| $ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} = \{ {\Lambda _{ij}}|i, j = 1, 2, 3\} = \left[ {\begin{array}{*{20}{l}} {0.8}&{0.1}&{0.1}\\ {0.1}&{0.8}&{0.1}\\ {0.1}&{0.1}&{0.8} \end{array}} \right] $ | （55） |
式中：Λij为从模型i转移到模型j的概率。
基于上述方法和仿真设置，分别进行300次的蒙特卡罗仿真，图 10~图 12给出了位置、速度和气动力加速度的均方根误差随时间的变化。可看出，当气动力加速度突变引起状态模型误差时，“Adaptive Singer”能够快速地降低估计误差，显著优于“Singer”，说明对自适应机动频率的有效性。在700 s以后，目标机动幅值有变化但较平缓，该方法表现出轻微的自适应过度，导致精度略低于其他几种方法，此时“Singer”和“Current Statistics”的精度稍高。“IMM-Singer”中由于采用IMM方法对机动频率进行自适应，因此其估计精度与“Adaptive Singer”最为接近。“Current Statistics”理论上也具有噪声方差的自适应能力，但是当加速度均值估计不准确时，可能适得其反，其结果劣于没有自适应能力的“Singer”。“Dynamics”对HGT阶跃机动的适应性最差。
|
| 图 10 位置RMSE(轨迹1) Fig. 10 RMSE of position (trajectory 1) |
|
| 图 11 速度RMSE(轨迹1) Fig. 11 RMSE of velocity (trajectory 1) |
|
| 图 12 气动力加速度RMSE(轨迹1) Fig. 12 RMSE of aerodynamic acceleration (trajectory 1) |
将雷达部署在经纬高为(23°, 5°, 0)的位置，与“轨迹1”相同，易得雷达可量测的轨迹弧段及相应状态量的变化，如图 13~图 15所示。相比于轨迹1，轨迹2的机动是连续的，但是机动幅度的变化强于轨迹1无倾侧角翻转的时刻。
|
| 图 13 可观测弧段(轨迹2) Fig. 13 Observable segment (trajectory 2) |
|
| 图 14 量测量真值(轨迹2) Fig. 14 Noise-free measurements (trajectory 2) |
|
| 图 15 雷达坐标系下气动力加速度真值(轨迹2) Fig. 15 True aerodynamic accelerations in radar frame (trajectory 2) |
为验证各方法的适应性，仿真设置与轨迹1中保持一致，图 16~图 18给出了300次蒙特卡罗仿真后位置、速度和气动力加速度的均方根误差随时间的变化。可看出，“Adaptive Singer”方法的滤波精度最高，对HGT机动幅度变化的适应性最好，特别是在目标机动幅度变化较大时，而“IMM-Singer”次之。由于机动幅度变化要强于轨迹一非阶跃机动的时刻，导致其他几种方法的估计精度较差，难以适应HGT多样的机动模式。
|
| 图 16 位置RMSE(轨迹2) Fig. 16 RMSE of position (trajectory 2) |
|
| 图 17 速度RMSE(轨迹2) Fig. 17 RMSE of velocity (trajectory 2) |
|
| 图 18 气动力加速度RMSE(轨迹2) Fig. 18 RMSE of aerodynamic acceleration (trajectory 2) |
通过对不同轨迹的跟踪滤波表明，所设计的“Adaptive Singer”对HGT的阶跃机动和连续幅值变化的机动均具有较好的适应性。图 19给出了各方法单次跟踪滤波的计算时间，可看出“IMM-Singer”约为其他方法的2倍，计算量大。相比于其他几种方法，“Adaptive Singer”不仅滤波精度高，而且计算量几乎没有增加，实时性好。
|
| 图 19 单次跟踪滤波计算时间 Fig. 19 Computational time for single tracking filtering |
为了精确跟踪机动模式多样的高超声速滑翔目标，基于正交性原理和UKF算法提出了一种机动频率自适应的跟踪方法。通过仿真验证，得到如下3点结论：
1) 由于HGT机动形式多样，几乎无法用固定的跟踪模型将其囊括，因此状态模型误差固然会存在。
2) 介于CV和CA模型之间的Singer模型可作为表征HGT机动的基础模型，再结合其他方法合理地调整其机动频率，降低模型误差。
3) 所提的“Adaptive Singer”方法能够较好地适应HGT的多种机动形式，与利用交互式多模型进行机动频率自适应的方法相比，跟踪精度高，计算量小。
5 附录 AF(X)的表达式为
| $ \begin{array}{l} \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{X}}) = [{F_1}, {F_2}, {F_3}, {F_4}, {F_5}, {F_6}, \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {F_7}, {F_8}, {F_9}{]^{\rm{T}}} \end{array} $ | （A1） |
式中:F1=vx，F2=vy，F3=vz，F4、F5和F6的表达式分别为
| $ {F_4} = \frac{{{g_{\rm{r}}}}}{r}{p_x} + {\omega ^2}{p_x} + 2\omega ({v_y}{\rm{sin}}{B_{\rm{o}}} - {v_z}{\rm{cos}}{B_{\rm{o}}}) + {a_x} $ | （A2） |
| $ \begin{array}{l} {F_5} = \frac{{{g_{\rm{r}}}}}{r}({p_y} - {R_{\rm{o}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}) + {g_\omega }{\rm{cos}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}} + {\omega ^2}[({p_y} - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{R_{\rm{o}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}){\rm{si}}{{\rm{n}}^2}{B_{\rm{o}}} - ({p_z} + {R_{\rm{o}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}) \cdot }\\ {{\rm{sin}}{B_{\rm{o}}}{\rm{cos}}{B_{\rm{o}}}] - 2\omega {v_x}{\rm{sin}}{B_{\rm{o}}} + {a_y}} \end{array} \end{array} $ | （A3） |
| $ \begin{array}{l} {F_6} = \frac{{{g_{\rm{r}}}}}{r}({p_z} + {R_{\rm{o}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}) + {g_\omega }{\rm{cos}}{\kern 1pt} {\kern 1pt} {B_{\rm{o}}} + {\omega ^2}[({p_z} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R_{\rm{o}}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}){\rm{co}}{{\rm{s}}^2}{B_{\rm{o}}} - ({p_y} - {R_{\rm{o}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\mu _{\rm{o}}}) \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{sin}}{\kern 1pt} {B_{\rm{o}}}{\rm{cos}}{\kern 1pt} {B_{\rm{o}}}] + 2\omega {v_x}{\rm{cos}}{\kern 1pt} {B_{\rm{o}}} + {a_z} \end{array} $ | （A4） |
F7、F8和F9的表达式为
| $ \left\{ {\begin{array}{*{20}{l}} {{F_7} = - {\lambda _x}{a_x}}\\ {{F_8} = - {\lambda _y}{a_y}}\\ {{F_9} = - {\lambda _z}{a_z}} \end{array}} \right. $ | （A5） |
若从动力学角度对HGT气动加速度进行建模，由于其采用倾侧转弯控制方式，则气动力加速度的形式为[1]
| $ \mathit{\boldsymbol{a}} = \left[ {\begin{array}{*{20}{l}} {{a_x}}\\ {{a_y}}\\ {{a_z}} \end{array}} \right] = \frac{1}{2}\rho {v^2}{\mathit{\boldsymbol{C}}_{{\rm{VtoR}}}}\left[ {\begin{array}{*{20}{c}} { - \eta {C_D}}\\ {\eta {C_L}}\\ 0 \end{array}} \right] $ | （B1） |
式中：ρ为大气密度，近似为高度的指数函数; v为相对速度v的模; CVtoR为目标速度坐标系至雷达坐标系的转换矩阵; η为目标的面质比。将式(B1)中的右边部分分为可用目标位置和速度表示的量和不可表示的量，具体形式为[7, 24]
| $ \mathit{\boldsymbol{a}} = \frac{1}{2}\rho {v^2}\mathit{\boldsymbol{M}}\left[ {\begin{array}{*{20}{c}} {{\eta _D}}\\ {{\eta _{LS}}}\\ {{\eta _{LC}}} \end{array}} \right] $ | （B2） |
式中：M的表达式为
| $ \mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{l}} { - \frac{{{v_x}}}{v}}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{v_{\rm{y}}}}}{{{v_{\rm{g}}}}}}&{ - \frac{{{v_x}{v_z}}}{{{v_{\rm{g}}}v}}}\\ { - \frac{{{v_{\rm{y}}}}}{v}}&{ - \frac{{{v_x}}}{{{v_{\rm{g}}}}}}&{ - \frac{{{v_{\rm{y}}}{v_z}}}{{{v_{\rm{g}}}v}}}\\ { - \frac{{{v_z}}}{v}}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{v_{\rm{g}}}}}{v}} \end{array}} \right] $ | （B3） |
式中：
| $ \left\{ {\begin{array}{*{20}{l}} {{\eta _D} = \eta {C_D}}\\ {{\eta _{LS}} = \eta {C_L}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {{\eta _{LC}} = \eta {C_L}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma } \end{array}} \right. $ | （B4） |
式中：γ为雷达坐标系至目标速度坐标系转换过程中的滚转角。通常情况下，HGV的面质比、气动系数和制导控制规律等先验信息均难以获知，因此广义气动系数是待估计的未知量。由于HGT在飞行过程中受到热流密度、动压、过载和控制量等因素的制约，其控制变量—攻角和倾侧角的变化应相对平缓，可用维纳随机过程来表征广义气动系数的变化，形式为[6]
| $ \left\{ {\begin{array}{*{20}{l}} {{{\dot \eta }_D} = {w_D}}\\ {{{\dot \eta }_{LS}} = {w_{LS}}}\\ {{{\dot \eta }_{LC}} = {w_{LC}}} \end{array}} \right. $ | （B5） |
式中：wD、wLS和wLC为互不相关的零均值高斯白噪声。将式(11)中ax、ay和az分别替换为ηD、ηLS和ηLC即为动力学模型的状态向量，其状态方程的形式与式(10)一致。
| [1] |
张远龙.基于三维剖面的滑翔飞行器弹道规划与制导方法研究[D].长沙: 国防科技大学, 2018: 1-11. ZHANG Y L. Research on entry trajectory generation for hypersonic glide vehicles based on three-dimensional profile[D]. Changsha: National University of Defense Technology, 2018: 1-11(in Chinese). |
| [2] |
雍恩米, 钱炜褀, 何开锋. 基于雷达跟踪仿真的滑翔式再入弹道突防性能分析[J]. 宇航学报, 2012, 33(10): 1370-1376. YONG E M, QIAN W Q, HE K F. Penetration ability analysis for glide reentry trajectory based on radar tracking[J]. Journal of Astronautics, 2012, 33(10): 1370-1376. (in Chinese) |
| Cited By in Cnki (20) | Click to display the text | |
| [3] | ZHANG K, XIONG J J, FU T T. Coupled dynamic model of state estimation for hypersonic glide vehicle[J]. Journal of System Engineering and Electronics, 2018, 29(6): 1284-1292. |
| Click to display the text | |
| [4] | LI X R, JILKOV V P. Survey of maneuvering target tracking. Part I:Dynamic models[J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(4): 1333-1364. |
| Click to display the text | |
| [5] |
吴楠, 陈磊. 高超声速滑翔再入飞行器弹道估计的自适应卡尔曼滤波[J]. 航空学报, 2013, 34(8): 1960-1971. WU N, CHEN L. Adaptive Kalman filtering for trajectory estimation of hypersonic glide reentry vehicles[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(8): 1960-1971. (in Chinese) |
| Cited By in Cnki (24) | Click to display the text | |
| [6] |
李广华.高超声速滑翔飞行器运动特性分析及弹道跟踪预报方法研究[D].长沙: 国防科学技术大学, 2016: 71-99. LI G H. Motion characteristics analysis and trajectory prediction for hypersonic glide vehicles[D]. Changsha: National University of Defense Technology, 2016: 71-99(in Chinese). |
| Cited By in Cnki (8) | Click to display the text | |
| [7] |
张凯, 熊家军, 韩春耀, 等. 一种基于气动力模型的高超声速滑翔目标跟踪算法[J]. 宇航学报, 2017, 38(2): 123-130. ZHANG K, XIONG J J, HAN C Y, et al. A tracking algorithm of hypersonic glide reentry vehicle via aerodynamic model[J]. Journal of Astronautics, 2017, 38(2): 123-130. (in Chinese) |
| Cited By in Cnki (15) | Click to display the text | |
| [8] |
王国宏, 李俊杰, 张翔宇, 等. 临近空间高超声速滑跃式机动目标的跟踪模型[J]. 航空学报, 2015, 36(7): 2400-2410. WANG G H, LI J J, ZHANG X Y, et al. A tracking model for near space hypersonic slippage leap maneuvering target[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(7): 2400-2410. (in Chinese) |
| Cited By in Cnki (40) | Click to display the text | |
| [9] |
李凡, 熊家军. 临近空间高超声速跳跃滑翔式目标自适应跟踪模型[J]. 航空学报, 2018, 39(12): 322355. LI F, XIONG J J. Adaptive tracking model for near space hypersonic jumping gliding target[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(12): 322355. (in Chinese) |
| Cited By in Cnki (2) | Click to display the text | |
| [10] |
张凯, 熊家军, 付婷婷, 等. 高超声速滑翔导弹气动参数自适应跟踪建模[J]. 国防科技大学学报, 2019, 41(1): 101-107. ZHANG K, XIONG J J, FU T T, et al. Aerodynamic parametric modeling of hypersonic gliding missile for adaptive tracking[J]. Journal of National University of Defense Technology, 2019, 41(1): 101-107. (in Chinese) |
| Cited By in Cnki (1) | Click to display the text | |
| [11] |
肖楚晗, 李炯, 雷虎民, 等. 基于AVSIMM算法的高超声速再入滑翔目标跟踪[J]. 北京航空航天大学学报, 2019, 45(2): 413-421. XIAO C H, LI J, LEI H M, et al. Hypersonic non-powered reentry gliding target tracking based on AVSIMM algorithm[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(2): 413-421. (in Chinese) |
| Cited By in Cnki (125) | Click to display the text | |
| [12] |
何睿智.高超声速助推滑翔飞行器全程弹道规划方法研究[D].长沙: 国防科技大学, 2017: 35-49. HE R Z. Study of all-course trajectory planning approach for hypersonic boost-glide vehicles[D]. Changsha: National University of Defense Technology, 2017: 35-49(in Chinese). |
| Cited By in Cnki (1) | Click to display the text | |
| [13] | MEHRA R K. Approaches to adaptive filtering[J]. IEEE Transactions on Automatic Control, 1972, 17(5): 693-698. |
| Click to display the text | |
| [14] |
刘畅, 杨锁昌, 汪连栋, 等. 基于自适应强跟踪CQKF的目标跟踪算法[J]. 北京航空航天大学学报, 2018, 44(5): 982-990. LIU C, YANG S C, WANG L D, et al. Target tracking algorithm based on adaptive strong tracking CQKF[J]. Journal of Beijing University of Aeronautics and Astronautics, 2018, 44(5): 982-990. (in Chinese) |
| Cited By in Cnki (3) | Click to display the text | |
| [15] | WANG N, LI L Y, WANG Q. Adaptive UKF-based parameter estimation for Bouc-Wen model of magnetorheological and elastomer materials[J]. Journal of Aerospace Engineering, 2019, 32(1): 0401830. |
| Click to display the text | |
| [16] | YUAN Y X, GAO W G. An optimal adaptive Kalman filter[J]. Journal of Geodesy, 2006, 80(4): 177-183. |
| Click to display the text | |
| [17] | HUANG Y L, ZHANG Y G, WU Z M, et al. A novel adaptive Kalman filter with inaccurate process and measurement noise covariance matrices[J]. IEEE Transactions on Automatic Control, 2018, 63(2): 594-601. |
| Click to display the text | |
| [18] | ZHOU D H, FRANK P M. Strong tracking filtering of nonlinear time-varying stochastic systems with coloured noise:Application to parameter estimation and empirical robustness analysis[J]. International Journal of Control, 1996, 65(2): 295-307. |
| Click to display the text | |
| [19] | WANG Y D, SUN S M, LI L. Adaptively robust unscented Kalman filter for tracking a maneuvering vehicle[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(5): 1696-1701. |
| Click to display the text | |
| [20] |
崔乃刚, 张龙, 王小刚, 等. 自适应高阶容积卡尔曼滤波在目标跟踪中的应用[J]. 航空学报, 2015, 36(12): 3885-3895. CUI N G, ZHANG L, WANG X G, et al. Application of adaptive high-degree cubature Kalman filter in target tracking[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(12): 3885-3895. (in Chinese) |
| Cited By in Cnki (52) | Click to display the text | |
| [21] | JIANG Y Z, MA P B, BAOYIN H X. Residual-normalized strong tracking filter for tracking a noncooperative maneuvering spacecraft[J]. Journal of Guidance, Control, and Dynamics, 2019, 42(10): 2304-2309. |
| Click to display the text | |
| [22] | ZHANG H W, XIE J W, GE J A, et al. Strong tracking SCKF based on adaptive CS model for manoeuvring aircraft tracking[J]. IET Radar, Sonar & Navigation, 2018, 12(7): 742-749. |
| Click to display the text | |
| [23] |
蒋冬婷, 宁静, 万洪容. 基于似然函数的自适应Singer模型滤波算法[J]. 西南师范大学学报(自然科学版), 2019, 44(1): 89-94. JIANG D T, NING J, WAN H R. An adaptive Singer model filter based on likelihood function[J]. Journal of Southwest China Normal University (Natural Science Edition), 2019, 44(1): 89-94. (in Chinese) |
| Cited By in Cnki (132) | Click to display the text | |
| [24] | LI X R, JILKOV V P. Survey of maneuvering target tracking. part Ⅱ:Motion models of ballistic and space targets[J]. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(1): 96-119. |
| Click to display the text | |
| [25] | SINGER R A. Estimating optimal tracking filter performance for manned maneuvering targets[J]. IEEE Transactions on Aerospace and Electronic Systems, 1970, 6(4): 473-483. |
| Click to display the text | |
| [26] |
赵琳. 非线性系统滤波理论[M]. 北京: 国防工业出版社, 2012: 52-134. ZHAO L. Nonlinear system filtering theory[M]. Beijing: National Defense Industry Press, 2012: 52-134. (in Chinese) |
| [27] | ZHOU H R, KUMAR K S P. A "current" statistical model and adaptive algorithm for estimating maneuvering targets[J]. Journal of Guidance, Control, and Dynamics, 1984, 7(5): 596-602. |
| Click to display the text | |
| [28] | LI X R, JILKOV V P. Survey of maneuvering target tracking. part V:Multiple-model methods[J]. IEEE Transactions on Aerospace and Electronic Systems, 2005, 41(4): 1255-1321. |
| Click to display the text |