蛇形机动是战斗机飞机飞行员在隐蔽接敌、机动规避及协同探测等战术动作中经常采用的战术机动类型^[1-2].在一定的空战态势下，飞行员进行蛇形机动反映了其战术意图，因此在空战中如何快速而又准确地识别出目标的蛇形机动模式对于判明其战术意图以及评估当前的战场态势具有十分重要的意义^[2-3].转弯角速度是识别蛇形机动模式最为关键的运动参量，通过精确辨识出机动目标的转弯角速度，可以进一步提高蛇形机动模式的识别率，然而，在实际空战过程中，由于机动目标的非合作特性，导致蛇形机动目标的转弯角速度无法被探测方直接获取，这也就为精确辨识蛇形机动目标的转弯角速度提出了迫切要求^[2].

目前，机动目标转弯角速度辨识的问题已经得到了许多学者的关注.文献[4-7]将转弯角速度作为额外的状态变量，从而对状态进行扩维处理，然后使用非线性滤波方法辨识角速度和估计状态变量.缺点是角速度估计的精度依赖于目标初始状态及协方差、过程噪声和量测噪声等因素^[2]；文献[8-10]分别根据转弯角速度与目标速度方向角、速度、加速度之间的物理关系，首先采用非线性滤波算法估计出目标速度方向角、速度、加速度等状态量，而后间接辨识出转弯角速度.Yuan等^[11]通过引入距离变化率量测，采用Cmin方法^[12]并结合航向角的变化估计出两个可能的转弯角速度的值，基于此设计包含匀速模型和两个匀速转弯模型的IMM算法进行状态估计.Frencl等^[13]通过引入距离变化率量测估计转弯角速度，并根据速度矢量在距离方向上的投影与距离变化率估计量之间的差异，对转弯角速度进行进一步修正以提高估计精度，然后依据最优线性无偏估计准则，利用粒子滤波进行状态估计.缺点是没有考虑角速度与目标状态之间相互耦合关系，增大了辨识风险^[14].

上述文献在估计目标状态和辨识转弯角速度时，都假设量测数据能够实时到达数据处理中心，但是，在实际的目标跟踪过程中，由于通信信道带宽的限制导致量测抵达数据处理中心时不可避免地存在随机延迟现象，从而造成用于估计与辨识的量测无法及时更新^[15-16]，并且在经典的估计理论当中，一般假设量测噪声与系统噪声为相互独立或相关的高斯白噪声，但是，在实际的空战过程中，由于复杂的作战环境以及目标本身所存在的散射现象，使得雷达的量测噪声不再服从高斯分布，而是尾部较长的“闪烁”噪声^[17-18].在实际空战中，由于蛇形机动目标的非合作特性，导致了蛇形机动目标的转弯角速度是未知且时变的，并且与目标状态相互耦合、彼此影响，一方面目标状态估计性能的提升有利于转弯角速度的辨识，另一方面转弯角速度辨识性能的提升有利于目标状态估计精度的提高，但是上述文献中所提到的转弯角速度辨识方法均没有考虑转弯角速度与目标状态之间的耦合关系以及这种耦合关系所带来的影响，并且均采取的是开环序贯处理方式，没有形成闭环反馈^[2].

针对上述问题，本文通过将目标状态估计与转弯角速度辨识联合考虑，基于EM算法框架，通过充分考虑量测随机延迟特性和“闪烁”噪声条件，重构了量测随机延迟下的粒子滤波器的似然函数，同时，为了避免采样粒子多样性丧失现象的发生，采用粒子群算法进一步优化采样过程，并结合后向模拟粒子平滑相关理论，设计了一步量测随机延迟条件下具有提前终止条件的后向模拟粒子平滑器用以计算后验平滑概率密度和联合分布密度，并采用牛顿迭代法的优化思想，不断迭代求得目标状态与转弯角速度闭环形式的联合优化解，并将滑窗思想引入到EM算法框架当中以降低算法的整体运行时间.该算法主要包括：E-step基于当前的角速度并采用粒子群优化的粒子滤波器与后向模拟粒子平滑器获得的后验平滑概率密度和联合分布密度，近似计算完整的对数函数似然函数的期望；M-step通过极大化似然函数进而更新获得下一次迭代的角速度估计量.

1 问题描述

蛇形机动可以分解成若干具有不同转弯角速度的圆弧转弯机动，因此蛇形机动目标转弯角速度的辨识问题相应地转化为圆弧转弯机动的角速度辨识问题.在量测一步随机延迟以及非高斯噪声背景下，假设目标在二维平面中运动，转弯机动模型状态方程为：

$ {\mathit{\boldsymbol{x}}_k} = {\mathit{\boldsymbol{f}}_{k - 1}}\left( {{\mathit{\boldsymbol{x}}_{k - 1}},{\mathit{\Omega }_{k - 1}}} \right) + {\mathit{\boldsymbol{w}}_{k - 1}}, $

$ \begin{array}{l} {\mathit{\boldsymbol{f}}_{k - 1}}\left( {{\mathit{\boldsymbol{x}}_{k - 1}},{\mathit{\Omega }_{k - 1}}} \right) = \\ \;\;\;\left[ {\begin{array}{*{20}{c}} 1&{\frac{{\sin \left( {{\mathit{\Omega }_{k - 1}}T} \right)}}{{{\mathit{\Omega }_{k - 1}}}}}&0&{\frac{{\cos \left( {{\mathit{\Omega }_{k - 1}}T} \right) - 1}}{{{\mathit{\Omega }_{k - 1}}}}}\\ 0&{\cos \left( {{\mathit{\Omega }_{k - 1}}T} \right)}&0&{ - \sin \left( {{\mathit{\Omega }_{k - 1}}T} \right)}\\ 0&{\frac{{1 - \cos \left( {{\mathit{\Omega }_{k - 1}}T} \right)}}{{{\mathit{\Omega }_{k - 1}}}}}&1&{\frac{{\sin \left( {{\mathit{\Omega }_{k - 1}}T} \right)}}{{{\mathit{\Omega }_{k - 1}}}}}\\ 0&{\sin \left( {{\mathit{\Omega }_{k - 1}}T} \right)}&0&{\cos \left( {{\mathit{\Omega }_{k - 1}}T} \right)} \end{array}} \right] \cdot {\mathit{\boldsymbol{x}}_{k - 1}}. \end{array} $

式中:x_k∈Rⁿ为状态变量; Ω_k-1为待辨识的转弯角速度; x_k=${\left[ {{x_k},{{\dot x}_k},{y_k},{{\dot y}_k}} \right]^{\rm{T}}}$;w_k∈Rⁿ为零均值高斯白噪声，满足E[w_kw_l^T]=Q_kδ_kl; δ_kl为克罗尼克δ函数.初始状态x₀满足x₀~N(x₀; ${{\mathit{\boldsymbol{\hat x}}}_0}$, P₀).

一步随机延迟量测模型^[19]为

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{z}}_k} = {\mathit{\boldsymbol{h}}_k}\left( {{\mathit{\boldsymbol{x}}_k}} \right) + {\mathit{\boldsymbol{v}}_k},\\ {\mathit{\boldsymbol{y}}_k} = \left( {1 - {\gamma _k}} \right){\mathit{\boldsymbol{z}}_k} + {\gamma _k}{\mathit{\boldsymbol{z}}_{k - 1}},k > 1,{\mathit{\boldsymbol{y}}_1} = {\mathit{\boldsymbol{z}}_1}. \end{array} \right. $

式中:k∈{1, 2, …}，z_k∈R^m为理想的量测向量; y_k∈R^m为实际的量测向量; h_k(·)为非线性量测函数; v_k∈R^m为闪烁噪声并与过程噪声w_k相互独立，即E[w_kv_k^T]=0;{γ_k∈R, k > 1}为已知的服从Bernoulli分布的随机变量，取值范围在0和1之间且满足:

$ \left\{ \begin{array}{l} p\left( {{\gamma _k} = 1} \right) = E\left[ {{\gamma _k}} \right] = \lambda ,\\ p\left( {{\gamma _k} = 0} \right) = 1 - E\left[ {{\gamma _k}} \right] = 1 - \lambda ,\\ E\left[ {\left( {{\gamma _k} - {p_k}} \right)2} \right] = \left( {1 - \lambda } \right)\lambda , \end{array} \right. $

初始状态x₀与过程噪声{w_k, k≥0}，{v_k, k≥1}以及{γ_k, k≥2}相互独立.λ为量测随机延迟概率，一般假设已知.

与高斯噪声相比，闪烁噪声具有“厚尾”特征，噪声分布尾部较长.雷达跟踪的闪烁噪声统计模型可由一个高斯分布噪声和一个附加残余噪声加权和表示，附加的残余噪声可由拉普拉斯分布、t分布、均匀分布、大方差的高斯分布等替代.

在本文当中采用两个高斯分布噪声的加权和表示雷达的闪烁噪声，其概率密度函数可表示为^[18]

$ p\left( v \right) = \left( {1 - \varepsilon } \right)N\left( {v;{\mu _1},{\mathit{\boldsymbol{P}}_1}} \right) + \varepsilon N\left( {v;{\mu _2},{\mathit{\boldsymbol{P}}_2}} \right). $

式中：ε∈[0, 1]为闪烁效应的强弱; N(ω; μ_t, P_t)为均值为μ_t、方差为P_t的高斯分布在v处的概率密度.则闪烁噪声的一、二阶距为：

$ \begin{array}{*{20}{c}} {\mu = E\left( v \right) = \left( {1 - \varepsilon } \right){\mu _1} + \varepsilon {\mu _2},}\\ {P = E\left( {\left( {v - \mu } \right){{\left( {v - \mu } \right)}^{\rm{T}}}} \right) = \left( {1 - \varepsilon } \right){P_1} + \varepsilon {P_2} + \tilde P,} \end{array} $

式中${\tilde P}$=(1-ε)μ₁μ₁^T+εμ₂μ₂^T-μμ^T.

2 基于PSOPF-EM的联合估计与辨识算法

为了解决量测一步随机延迟和非高斯条件下蛇形机动目标的未知转弯角速度辨识和目标状态估计问题，本文通过将角速度辨识与目标状态估计联合考虑，提出了基于PSOPF-EM的联合估计与辨识算法.为了降低算法的整体运行时间成本，本文通过对后向模拟粒子平滑算法设置终止条件来进行算法优化，从而降低了粒子平滑算法的时间消耗，同时又引入了滑窗思想，窗口长度为l，其长度大小可依据辨识与估计精度的实际需求自行设定，用于均衡运算时间成本与精度之间的“博弈”.具体的算法框架如图 1所示.

EM算法通过极大化对数似然函数L_Ω(Y_k-l^k|Y₁^k-l-1)，求得角速度的估计值${\mathit{\hat \Omega }}$，即

$ \begin{array}{*{20}{c}} {\mathit{\hat \Omega } = \arg \max {L_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) \buildrel \Delta \over = }\\ {\arg \max \log {p_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right),} \end{array} $

式中, Y₁^k={y₁, y₂, …, y_k}为实际的延迟量测.

由于一步随机延迟量测取决于未知的状态向量，因此直接极大化不完备数据的对数似然函数L_Ω(Y_k-l^k|Y₁^k-l-1)是不可能实现的.为了解决这个问题，EM算法将缺失状态向量X_k-l-1^k引入到对数似然函数`中，构成完备数据的对数似然函数L_Ω(X_k-l-1^k, Y_k-l^k|Y₁^k-l-1)，其中X_k-l-1^k={x_k-l-1, x_k-l, …, x_k}.

依据贝叶斯准则可知

$ \begin{array}{*{20}{c}} {{p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k,\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) = }\\ {{p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right){p_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right),} \end{array} $

则完备数据的对数似然函数可表示为^[20]

$ \begin{array}{l} {L_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k,\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) = \\ \log \left( {{p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k,\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)} \right) = \\ \log \left( {{p_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right)} \right) + \log \left( {{p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)} \right) = \\ {L_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) + \log \left( {{p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)} \right). \end{array} $

给定第t次迭代时角速度的参数估计值${\mathit{\hat \Omega }_t}$，则完备数据对数似然函数的条件期望函数为

$ \begin{array}{l} Q\left( {\mathit{\Omega },{{\mathit{\hat \Omega }}_t}} \right) \buildrel \Delta \over = {E_{{{\hat \Omega }_t}}}\left\{ {{L_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k,\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right)\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right\} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{L_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) + \int {\log \left( {{p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)} \right)} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{p_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right){\rm{d}}\mathit{\boldsymbol{X}}_{k - l - 1}^k, \end{array} $

则

$ \begin{array}{*{20}{c}} {{L_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) - {L_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) = }\\ {Q\left( {\mathit{\Omega },{{\mathit{\hat \Omega }}_t}} \right) - Q\left( {{{\mathit{\hat \Omega }}_t},{{\mathit{\hat \Omega }}_t}} \right) + }\\ {\int {\log \frac{{{p_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)}}{{{p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)}}{p_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right){\rm{d}}\mathit{\boldsymbol{X}}_{k - l - 1}^k,} } \end{array} $

上式最右边的积分被称为Kullback-Leibler散度^[21]，其值非负.

则可知

$ \begin{array}{*{20}{c}} {{L_\Omega }\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) - {L_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) \ge }\\ {Q\left( {\mathit{\Omega },{{\mathit{\hat \Omega }}_t}} \right) - Q\left( {{{\mathit{\hat \Omega }}_t},{{\mathit{\hat \Omega }}_t}} \right).} \end{array} $

(1)

从式(1)可以看出如果选择Ω的一个合适值${\mathit{\hat \Omega }_{t+1}}$使得Q(${\mathit{\hat \Omega }_{t+1}}$, ${\mathit{\hat \Omega }_t}$ > Q(${\mathit{\hat \Omega }_t}$, ${\mathit{\hat \Omega }_t}$)，则相应的对数似然函数满足$L_{{\hat \Omega }_{t+1}}$(Y_k-l^k|Y₁^k-l-1) > $L_{{\hat \Omega }_t}$(Y_k-l^k|Y₁^k-l-1).也就是说EM算法每次迭代时通过获得一个角速度值${\mathit{\hat \Omega }_{t+1}}$使得函数Q(Ω, ${\mathit{\hat \Omega }_t}$)极大化，同时也能保证对数似然函数L_Ω(Y_k-l^k|Y₁^k-l-1)在每次迭代中也是增加的.并且随着迭代次数的增加，${\mathit{\hat \Omega }_t}$也将逼近角速度的真实值，即$\mathop {\lim }\limits_{t \to \infty } {{\mathit{\hat \Omega }}_t} = \mathit{\hat \Omega }$.

根据贝叶斯准则可知

$ \begin{array}{l} {p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k,\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) = \\ \left( {\prod\limits_{i = 0}^l {{p_\Omega }\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)p\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} } \right).\\ {p_\Omega }\left( {{\mathit{\boldsymbol{x}}_{k - l - 1}}\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) = \left( {\prod\limits_{i = 0}^l {\left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.} \right)} + } \right.\\ \left. {\lambda {p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} \right) \times \left( {\prod\limits_{i = 0}^l {p\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} } \right).\\ {p_\Omega }\left( {{\mathit{\boldsymbol{x}}_{k - l - 1}}\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right), \end{array} $

(2)

其中

$ \begin{array}{l} {p_\Omega }\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}},{x_{k - i - 1}}} \right.} \right) = \\ {p_\Omega }\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\gamma _{k - i}} = 0,{\mathit{\boldsymbol{x}}_{k - i}},{x_{k - i - 1}}} \right.} \right)p\left( {{\gamma _{k - i}} = 0} \right) + \\ {p_\Omega }\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\gamma _{k - i}} = 1,{\mathit{\boldsymbol{x}}_{k - i}},{x_{k - i - 1}}} \right.} \right)p\left( {{\gamma _{k - i}} = 1} \right) = \\ \left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{y}}_{k - i}} - {\mathit{\boldsymbol{h}}_{k - i}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}} \right)} \right) + \\ \lambda {p_\Omega }\left( {{\mathit{\boldsymbol{y}}_{k - i}} - {\mathit{\boldsymbol{h}}_{k - i - 1}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right)} \right) = \\ \left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.} \right) + \lambda {p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right), \end{array} $

则完备数据对数似然函数可以分解成如下

$ \begin{array}{l} {L_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k,\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) = \\ \log {p_\Omega }\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^k,\mathit{\boldsymbol{Y}}_{k - l}^k\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) = \\ \log {p_\Omega }\left( {{\mathit{\boldsymbol{x}}_{k - l - 1}}\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right) + \sum\limits_{i = 0}^l {\log {p_\Omega }\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} + \\ \sum\limits_{i = 0}^l {\log \left[ {\left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.} \right) + \lambda {p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} \right]} , \end{array} $

则Q(Ω, ${\mathit{\hat \Omega }_t}$)可表示成如下形式：

$ Q\left( {\mathit{\Omega },{{\mathit{\hat \Omega }}_t}} \right) = {I_1} + {I_2} + {I_3}, $

其中：

$ {I_1} = \int {\log {p_\Omega }\left( {{\mathit{\boldsymbol{x}}_{k - l - 1}}\left| {\mathit{\boldsymbol{Y}}_1^{k - l - 1}} \right.} \right){p_{{{\mathit{\hat \Omega }}_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - l - 1}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right){\rm{d}}{x_{k - l - 1}}} , $

(3)

$ \begin{gathered} {I_2} = \sum\limits_{i = 0}^I {\iint {\log {p_\Omega }\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right){p_{{{\mathit{\hat \Omega }}_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)}} \times \hfill \\ \;\;\;\;\;\;{\text{d}}{\mathit{\boldsymbol{x}}_{k - i}}{\text{d}}{\mathit{\boldsymbol{x}}_{k - i - 1}}, \hfill \\ \end{gathered} $

(4)

$ \begin{gathered} {I_3} = \sum\limits_{i = 0}^I {\iint {\log \left[ {\left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.} \right) + } \right.}} \hfill \\ \;\;\;\;\;\;\left. {\lambda {p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} \right] \times \hfill \\ \;\;\;\;\;\;{p_{{{\mathit{\hat \Omega }}_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right){\text{d}}{\mathit{\boldsymbol{x}}_{k - i}}{\text{d}}{\mathit{\boldsymbol{x}}_{k - i - 1}}. \hfill \\ \end{gathered} $

(5)

式中：当l=k-1时，I₁为一个常数，与待辨识量θ无关；当0≤l＜k-1时，I₁与θ相关性较弱，而且，为了满足实时性的要求，辨识算法应尽量关注当前时刻的数据或者邻近时刻的数据^[2]，与此同时，由于蛇形机动目标的转弯角速度耦合在系统状态转移矩阵当中，即转弯角速度只与先验概率分布p_Ω(x_k-i|x_k-i-1)有关，故可知在辨识转弯角速度时，只需要考虑条件期望函数的I₂部分即可.

EM算法通过不断迭代求得Q(Ω, ${\mathit{\hat \Omega }_t}$)的极大似然估计，每次迭代过程包含E-step和M-step.

2.1 E-step

E-step的主要任务是在给定第t次迭代时的角速度估计值${\mathit{\hat \Omega }_t}$，计算条件期望函数Q(Ω, ${\mathit{\hat \Omega }_t}$).计算Q(Ω, ${\mathit{\hat \Omega }_t}$)的关键在于获取平滑概率分布$p_{{\mathit{\hat \Omega }_t}}$(x_k-i|Y₁^k)和联合平滑概率分布$p_{{\mathit{\hat \Omega }_t}}$(x_k-i, x_k-i-1|Y₁^k)，其中平滑概率分布$p_{{\mathit{\hat \Omega }_t}}$(x_k-i|Y₁^k)可以采用量测一步随机延迟下基于粒子群改进的粒子滤波器与后向模拟粒子平滑器求得，而联合平滑概率分布$p_{{\mathit{\hat \Omega }_t}}$(x_k-i, x_k-i-1|Y₁^k)可通过如下方式求得.

依据贝叶斯准则，联合平滑密度分布可表示成如下形式：

$ \begin{array}{l} {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right) = {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.,\mathit{\boldsymbol{Y}}_1^k} \right) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right). \end{array} $

(6)

根据文献[21]可知

$ {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.,\mathit{\boldsymbol{Y}}_1^k} \right) = {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.,\mathit{\boldsymbol{Y}}_1^{k - i}} \right), $

(7)

其中

$ \begin{array}{l} {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.,\mathit{\boldsymbol{Y}}_1^{k - i}} \right) = \\ \frac{{{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}},{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{y}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}}{{{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{y}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}} = \\ \frac{{{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.,{\mathit{\boldsymbol{x}}_{k - i}}} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}}{{{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{y}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}}. \end{array} $

(8)

将式(7)~(8)带入式(6)中可得

$ \begin{array}{l} {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right) = {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) \times \\ \frac{{{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.,{\mathit{\boldsymbol{x}}_{k - i}}} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right)}}{{{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{y}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}}, \end{array} $

其中

$ \begin{array}{l} {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{y}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) = \\ \int {{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}},{\mathit{\boldsymbol{x}}_{k - i}},{\mathit{\boldsymbol{y}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{k - i - 1}}} = \\ \int {{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.,{\mathit{\boldsymbol{x}}_{k - i}}} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} \cdot \\ {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{k - i - 1}}. \end{array} $

(9)

式中：${p_{{{\hat \Omega }_t}}}$(x_k-i-1|Y₁^k-i-1)、${p_{{{\hat \Omega }_t}}}$(x_k-i|Y₁^k)分别为滤波概率密度和平滑概率密度.

滤波概率密度${p_{{{\hat \Omega }_t}}}$(x_k-i-1|Y₁^k-i-1)可采用基于序贯重要性采样的粒子滤波器^[22]求得.从重要性概率密度q(X_k-l-1^k-i-1|Y₁^k-i-1)进行采样，进而近似逼近滤波概率密度${p_{{{\hat \Omega }_t}}}$(X_k-l-1^k-i-1|Y₁^k-i-1)为

$ \begin{array}{*{20}{c}} {{p_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^{k - i - 1}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) \approx {{\hat p}_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^{k - i - 1}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) = }\\ {\sum\limits_{j = 1}^N {\bar w_{k - l - 1}^{k - i - 1,j}\delta \left[ {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 1} - \mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 1,j}} \right]} .} \end{array} $

(10)

式中：N为采样粒子数量; {x_k-l-1^{k-i-1, j}}_j=1^N为采样粒子集合, i∈[0, l]; w_k-l-1^{k-i-1, j}为采样粒子对应的归一化权重.w_k-l-1^{k-i-1, j}为重要性权重，并且满足:

$ w_{k - l - 1}^{k - i - 1,j} = \frac{{{p_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{X}}_{k - l - 1}^{k - i - 1,j}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}}{{q\left( {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 1,j}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}},\tilde w_{k - l - 1}^{k - i - 1,j} = \frac{{w_{k - l - 1}^{k - i - 1,j}}}{{\sum\limits_{j = 1}^N {w_{k - l - 1}^{k - i - 1,j}} }}. $

(11)

依据贝叶斯准则可知

$ \begin{array}{*{20}{c}} {q\left( {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 1,j}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) = }\\ {q\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 2},\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) \cdot q\left( {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 2}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right).} \end{array} $

(12)

依据贝叶斯准则和Markov性可知

$ \begin{array}{l} {p_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 1}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) = {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}},{\mathit{\boldsymbol{x}}_{k - i - 2}}} \right.} \right) \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{p_{{{\mathit{\hat \Omega }}_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 2}}} \right.} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{p_{{{\mathit{\hat \Omega }}_t}}}\left( {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 2}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 2}} \right.} \right), \end{array} $

(13)

则将式(12)~(13)带入到式(11)中可得

$ \begin{array}{l} {w_{k - i - 1}} = \\ \frac{{{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.,{x_{k - i - 2}}} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 2}}} \right.} \right)}}{{\mathit{\boldsymbol{q}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 2},\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right)}}{w_{k - i - 2}}. \end{array} $

(14)

为了便于重要性权重的计算，在此选择先验转移密度作为建议密度，即

$ q\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{x}}_{k - l - 1}^{k - i - 2},\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) = {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 2}}} \right.} \right), $

(15)

将式(15)带入式(14)中可得

$ {w_{k - i - 1}} = {w_{k - i - 2}}{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.,{x_{k - i - 2}}} \right). $

(16)

将式(2)带入式(16)中可知粒子x_k-i-1^j对应的权重为

$ \begin{array}{l} w_{k - i - 1}^j = {w_{k - i - 2}}\left[ {\left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {\lambda {p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i - 1}}} \right.} \right)} \right]. \end{array} $

(17)

将采样粒子x_k-i-1^j以及相应的归一化权重w_k-i-1^j带入到式(10)中可得滤波概率密度${p_{{{\hat \Omega }_t}}}$(x_k-i-1|Y₁^k-i-1)，即

$ \begin{array}{l} {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {\mathit{\boldsymbol{Y}}_1^{k - i - 1}} \right.} \right) \approx {{\hat p}_{{{\mathit{\hat \Omega }}_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {Y_1^{k - i - 1}} \right.} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{j = 1}^N {\bar w_{k - l - 1}^j\delta \left[ {{\mathit{\boldsymbol{x}}_{k - i - 1}} - \mathit{\boldsymbol{x}}_{k - l - 1}^j} \right]} . \end{array} $

(18)

由于传统的基于序贯重要性采样的粒子滤波算法采用了次优的重要性函数^[23]，所以，存在以下两个缺陷：1)粒子贫乏问题；2)为了确保状态估计精度需要大量采样粒子，造成算法的运行时间成本较高.为解决上述问题，本文将粒子群优化算法引入到粒子滤波算法当中，使得粒子群在权重值更新之前更加趋向高似然区域，从而缓解粒子匮乏现象，可以减少滤波所需的粒子数量.粒子群算法与粒子滤波算法的结合主要是利用粒子所经历过的最优状态值pbest和粒子群中的最优粒子的状态值gbest，通过下式进行更新^[23]：

$ V_k^i = \left| {{\rm{randn}}} \right|\left( {{\rm{pbest}} - \mathit{\boldsymbol{x}}_k^i} \right) + \left| {{\rm{Randn}}} \right|\left( {{\rm{gbest}} - \mathit{\boldsymbol{x}}_k^i} \right), $

(19)

$ \mathit{\boldsymbol{x}}_k^i = \mathit{\boldsymbol{x}}_k^i + V_k^i. $

(20)

式中：|randn|、|Randn|分别为正的高斯分布的随机数，可由abs[N(0, 1)]产生；pbest为当前粒子所经历的状态值中适应度函数最大的粒子的状态值；gbest为粒子群中适应度函数值最大的粒子的状态值^[24].当粒子群的最优值符合某阈值ε时，则表明粒子群中的粒子已经分布在真实状态附近，那么粒子群算法将停止优化.此时再对粒子集依据式(17)进行权重更新并依据式(11)进行归一化处理.具体的算法过程，见图 2.

考虑到后向模拟粒子平滑器的在滤波效果不好的条件下仍能收敛的良好的性质^[25]，故本文拟采用该平滑器近似计算平滑概率密度${p_{{{\hat \Omega }_t}}}$(x_k-i|Y₁^k).

将式(18)带入到式(8)、(9)可得：

$ \begin{array}{*{20}{c}} {{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i - 1}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.,\mathit{\boldsymbol{Y}}_1^{k - i}} \right) \approx }\\ {\sum\limits_{j = 1}^N {\bar w_{k - i - 1\left| {k - i} \right.}^j\delta \left[ {{\mathit{\boldsymbol{x}}_{k - i - 1}} - \mathit{\boldsymbol{x}}_{k - l - 1}^j} \right]} .} \end{array} $

(21)

式中

$ \begin{array}{l} \bar w_{k - i - 1\left| {k - i} \right.}^j = \\ \frac{{\bar w_{k - i - 1}^j{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.,\mathit{\boldsymbol{x}}_{k - i - 1}^j} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {\mathit{\boldsymbol{x}}_{k - i - 1}^j} \right.} \right)}}{{\sum\limits_{r = 1}^N {\bar w_{k - i - 1}^r{p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{y}}_{k - i}}\left| {{\mathit{\boldsymbol{x}}_{k - i}}} \right.,\mathit{\boldsymbol{x}}_{k - i - 1}^r} \right){p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {\mathit{\boldsymbol{x}}_{k - i - 1}^r} \right.} \right)} }}. \end{array} $

(22)

假设已经从平滑分布${p_{{{\hat \Omega }_t}}}$(x_k-i^k|Y₁^k)中采样得到一组粒子集合$\tilde x_{k - i}^k$，则k-i-1时刻的粒子${{\tilde x}_{k - i - 1}}$可根据式(21)、(22)计算得到的${p_{{{\hat \Omega }_t}}}$(x_k-i-1|${{\tilde x}_{k - i}}$, Y₁^k)分布采样获得.更新后的粒子集合$\left( {{{\tilde x}_{k - i - 1}}, \tilde x_{k - i}^k} \right)$被用于近似逼近平滑密度${p_{{{\hat \Omega }_t}}}$(x_k-i-1:k|Y₁^k).采用此方法，完成[k-l-1, k]时间段之间的粒子采样.这个过程重复M次，即对m=1, 2, …, M，得到粒子集合$\left\{ {\mathit{\boldsymbol{\tilde x}}_{k - i - 1}^{k, m}} \right\}_{m = 1}^M$，则平滑分布的逼近式为^[26-27]

$ {p_{{{\hat \Omega }_t}}}\left( {\mathit{\boldsymbol{x}}_{k - l - 1}^k\left| {\mathit{\boldsymbol{Z}}_1^k} \right.} \right) \approx \frac{1}{M}\sum\limits_{m = 1}^M {\delta \left( {\mathit{\boldsymbol{x}}_{k - l - 1}^k - \mathit{\boldsymbol{\tilde x}}_{k - l - 1}^{k,m}} \right)} . $

则平滑概率密度${p_{{{\hat \Omega }_t}}}$(x_k-i|Y₁^k)可表示成如下形式：

$ {p_{{{\hat \Omega }_t}}}\left( {{\mathit{\boldsymbol{x}}_{k - i}}\left| {\mathit{\boldsymbol{Y}}_1^k} \right.} \right) \approx \frac{1}{M}\sum\limits_{m = 1}^M {\delta \left( {{\mathit{\boldsymbol{x}}_{k - i}} - \mathit{\boldsymbol{\tilde x}}_{k - i}^m} \right)} . $

标准的后向模拟粒子平滑算法过程见图 3.

从图 3中可以明显看出，在进行多项分布采样时，每采样得到一个粒子索引值，就需要繁琐地计算采样粒子的平滑权重，导致算法的计算量呈二次方式增加.为了避免粒子平滑权重的复杂计算，降低运算时间成本，Douc等^[28]通过将拒绝采样(RS)方法引入到标准的后向模拟平滑算法中来，进行粒子索引值的采样，降低了时间复杂度，提升了算法运行效率.

假设系统状态转移密度函数存在上界，即p_Ω(x_k-i|x_k-i-1)≤ρ.基于这个假设，直接将粒子滤波得到的粒子权重作为多项分布采样的采样概率，即依据该概率$\left\{ {\bar w_{k - i - 1}^j} \right\}_{j = 1}^N$进行多项式采样获得了粒子索引值I(m)~Cat($\left\{ {\bar w_{k - i - 1}^j} \right\}_{j = 1}^N$)，依据粒子索引值集合可以获得采样粒子集合，相应的权重为$\left\{ {\tilde w_{k - i - 1|k - i}^{j, m}} \right\}_{j = 1}^N $.则基于拒绝采样思想采样时粒子被接受的概率可由下式获得：

$ \frac{{\tilde w_{k - i - 1\left| {k - i} \right.}^{I\left( m \right),m}}}{{w_{k - i - 1}^{I\left( m \right)}}} \propto {p_\Omega }\left( {\mathit{\boldsymbol{\tilde x}}_{k - i}^m\left| {\mathit{\boldsymbol{x}}_{k - i - 1}^{I\left( m \right)}} \right.} \right) \le \rho , $

也就是说，通过建议分布采样得到的粒子如果满足p_Ω($\mathit{\boldsymbol{\tilde x}}_{k - i}^m$|x_k-i-1^I(m))≤ρ则会被接受；否则被拒绝，然后重新按照上述过程进行采样，直至满足接受条件为止，然后再进行下一个平滑粒子的采样.采用这种方式虽然从理论上来讲避免了粒子平滑权重的反复的计算，降低了标准后向模拟粒子平滑算法的运行时间成本，但是在实际使用拒绝采样法获取平滑粒子的过程中，有可能存在一些采样粒子无法满足接受条件，导致采样过程持续进行，甚至会导致死循环，这些低接受概率的粒子对于算法精度的贡献率不仅不高，反而造成算法运行时间成本的大幅增加.为了改善这种情况，本文依据拒绝采样法的后向模拟粒子平滑器以及标准后向模拟粒子平滑算法两者的运算特点，提出一种将二者混合的平滑算法，即为拒绝采样设定最大允许迭代次数K(K值的设定涉及到运算成本与精度的“博弈”，可以是提前设定的常数或者一个时间序列集合，也可以是依据采样过程的实际情况在线设定的自适应值)，达到最大迭代次数后，将被拒绝的采样粒子输入到图 3算法2中再进行计算采样，算法过程见图 3.

将通过图 3计算得到的联合平滑概率分布函数${p_{{{\hat \Omega }_t}}}$(X_k-l-1^k|Y₁^k)，带入到式(3)~(5)中，得到I₁、I₂和I₃的近似表达式为：

$ {{\hat I}_1} \approx \frac{1}{M}\sum\limits_{j = 1}^M {\log {p_\Omega }\left( {\mathit{\boldsymbol{\tilde x}}_{k - l - 1}^j} \right)} , $

(23)

$ {{\hat I}_2} \approx \sum\limits_{i = 0}^l {\frac{1}{M}\sum\limits_{j = 1}^M {\log {p_\Omega }\left( {\mathit{\boldsymbol{\tilde x}}_{k - i}^j\left| {\mathit{\boldsymbol{\tilde x}}_{k - i - 1}^j} \right.} \right)} } , $

(24)

$ \begin{array}{l} {{\hat I}_3} \approx \sum\limits_{i = 0}^l {\frac{1}{M}\sum\limits_{j = 1}^M {\log \left[ {\left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i}}\left| {\mathit{\boldsymbol{\tilde x}}_{k - i}^j} \right.} \right) + } \right.} } \\ \;\;\;\;\;\;\left. {\lambda {p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i - 1}}\left| {\mathit{\boldsymbol{\tilde x}}_{k - i - 1}^j} \right.} \right)} \right]. \end{array} $

(25)

将式(23)~(25)代入式(2)即可得到${\mathit{\hat \Omega }}$(Ω, ${\mathit{\hat \Omega }_t}$)的近似表达式.

2.2 M-step

M-step的主要任务是解决Q(Ω, ${\mathit{\hat \Omega }_t}$)的极大化问题，即求解使Q(Ω, ${\mathit{\hat \Omega }_t}$)满足极大值时所对应的${\mathit{\hat \Omega }_{t+1}}$值，用于EM算法的下一次迭代更新，在此采用${\mathit{\hat Q }}$(Ω, ${\mathit{\hat \Omega }_t}$)代替Q(Ω, ${\mathit{\hat \Omega }_t}$)为

$ {{\mathit{\hat \Omega }}_{t + 1}} = \arg \mathop {\max }\limits_\mathit{\Omega } Q\left( {\Omega ,{{\mathit{\hat \Omega }}_t}} \right). $

本文采用梯度下降法求取${\mathit{\hat Q}}$(Ω, ${\mathit{\hat \Omega }_t}$)的极大值，即：

$ \frac{\partial }{{\partial \mathit{\Omega }}}\hat Q\left( {\mathit{\Omega },{{\mathit{\hat \Omega }}_t}} \right) = \frac{\partial }{{\partial \mathit{\Omega }}}{{\hat I}_1} + \frac{\partial }{{\partial \mathit{\Omega }}}{{\hat I}_2} + \frac{\partial }{{\partial \mathit{\Omega }}}{{\hat I}_3}, $

(26)

$ \frac{\partial }{{\partial \mathit{\Omega }}}{{\hat I}_1} = \frac{1}{M}\sum\limits_{j = 1}^M {\frac{{\partial \log {p_\Omega }\left( {\mathit{\boldsymbol{x}}_{k - l - 1}^j} \right)}}{{\partial \mathit{\Omega }}}} , $

(27)

$ \frac{\partial }{{\partial \mathit{\Omega }}}{{\hat I}_2} = \frac{1}{M}\sum\limits_{i = 0}^l {\sum\limits_{j = 1}^M {\frac{{\partial \log {p_\Omega }\left( {\mathit{\boldsymbol{x}}_{k - i}^j\left| {\mathit{\boldsymbol{x}}_{k - i - 1}^j} \right.} \right)}}{{\partial \mathit{\Omega }}}} } , $

(28)

$ \begin{array}{*{20}{c}} {\frac{\partial }{{\partial \mathit{\Omega }}}{{\hat I}_3} = \frac{1}{M}\sum\limits_{i = 0}^l {\sum\limits_{j = 1}^M {} } }\\ {\frac{{\partial \log \left[ {\left( {1 - \lambda } \right){p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i}}\left| {\mathit{\boldsymbol{\tilde x}}_{k - i}^j} \right.} \right) + \lambda {p_\Omega }\left( {{\mathit{\boldsymbol{z}}_{k - i - 1}}\left| {\mathit{\boldsymbol{\tilde x}}_{k - i - 1}^j} \right.} \right)} \right]}}{{\partial \mathit{\Omega }}}.} \end{array} $

(29)

本文采用的是滞后窗口，即辨识k-l时刻的转弯角速度Ω的值需要使用[k-l, k]时刻的量测.

可将本文算法总结如下，见图 5.

3 仿真分析

本文利用水平方向上的转弯机动非线性动态模型，仿真出一条蛇形机动轨迹^[3].假设机动目标在1~80 s以Ω₁=-1.5°/s作转弯运动，在k=81 s时转弯角速度突变为Ω₂=2°/s作转弯运动，并持续到200 s.设置转弯角速度初始值为Ω₀=-0.5°/s，采样周期T=1，q₁=0.1 m²/s³，w_k是零均值高斯白噪声，其协方差为Q.

$ \mathit{\boldsymbol{Q}} = {\rm{diag}}\left( {{q_1}\mathit{\boldsymbol{M}},{q_1}\mathit{\boldsymbol{M}}} \right),\mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {{T^3}/3}&{{T^2}/2}\\ {{T^2}/2}&T \end{array}} \right], $

载机雷达的非线性量测方程为

$ {\mathit{\boldsymbol{z}}_k} = \mathit{\boldsymbol{h}}\left( {{\mathit{\boldsymbol{x}}_k}} \right) + {\mathit{\boldsymbol{v}}_k} = \left[ \begin{array}{l} \sqrt {x_k^2 + y_k^2} \\ \arctan \left( {{y_k}/{x_k}} \right) \end{array} \right] + {\mathit{\boldsymbol{v}}_k}. $

量测噪声v_k是采用两个高斯分布的加权和合成的闪烁噪声，具体参数设置如下：两个高斯分布均值均为0, 协方差分别为R_{1, k}=[σ_{1, r}²σ_{1, φ}²]，R_{2, k}=[σ_{2, r}²σ_{2, φ}²]σ_{1, r}=10 m，σ_{1, φ}=$ \sqrt 10 $ mrad，σ_{2, r}=100 m，σ_{2, φ}=10 mrad，ε=0.3.量测一步任意延迟概率一般假设为已知量，在此设为λ=0.3.

本文所提出的基于EM框架的联合估计与辨识算法和扩维法的初始状态以及协方差分别设置为：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_0} = {{\left[ {\begin{array}{*{20}{c}} {10\;{\rm{km}}}&{0.3\;{\rm{km/s}}}&{10\;{\rm{km}}}&{0.3\;{\rm{km/s}}} \end{array}} \right]}^{\rm{T}}},}\\ {{\mathit{\boldsymbol{P}}_0} = {\rm{diag}}\left[ {\begin{array}{*{20}{c}} {100\;{{\rm{m}}^2}}&{10\;{{\rm{m}}^2}/{{\rm{s}}^2}}&{100\;{{\rm{m}}^2}}&{10\;{{\rm{m}}^2}/{{\rm{s}}^2}} \end{array}} \right].} \end{array} $

扩维法初始状态和协方差的设置如下：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{x}}_0^a = {{\left[ {10\;{\rm{km}}\;{\rm{0}}{\rm{.3}}\;{\rm{km/s}}\;{\rm{10}}\;{\rm{km}}\;{\rm{0}}{\rm{.3}}\;{\rm{km/s}}\; - {\rm{0}}{\rm{.}}{{\rm{5}}^ \circ }{\rm{/s}}} \right]}^{\rm{T}}},}\\ {P_0^a = {\rm{diag}}\left[ {100\;{{\rm{m}}^2}\;10\;{{\rm{m}}^2}/{{\rm{s}}^2}\;100\;{{\rm{m}}^2}\;10\;{{\rm{m}}^2}/{{\rm{s}}^2}100\left( {{\rm{mrad/s}}} \right)} \right],} \end{array} $

为了评估分析本文提出的算法性能，首先将本文提出的算法与传统的扩维法进行对比分析.本文提出的联合估计与辨识算法采用采取滑动滞后窗口策略，即使用未来几拍的量测估计与辨识当前目标状态和转弯角速度，出于目标跟踪实时性的要求，在保证辨识和估计精度的基础上应选择较小的滑动滞后窗口长度，一般滑动滞后窗口选择为5，其最大迭代次数为5次，粒子数为100，标准的粒子滤波算法、带一步任意量测的PF算法以及带一步任意量测的PSO-PF算法粒子数均为1 000，量测任意延迟概率均为0.3，各执行50次蒙特卡洛仿真.图 6、7分别显示了上述两类方法对角速度辨识和目标状态估计的效果，从图 6、7中可以看出，本文提出的联合估计与辨识算法在角速度辨识和目标状态估计上比传统的扩维法误差小、精度高，也就是说本文提出的算法在使用了很少的粒子情况下仍然能够保证目标状态估计与角速度辨识的收敛性，确保估计与辨识的精度，究其原因主要有:1)本文提出的一步量测随机延迟条件下基于粒子群改进的粒子滤波与后向模拟粒子平滑算法的良好的收敛性使然.本文首先通过采用粒子群算法改进粒子采样过程，充分保持了采样粒子的多样性，改善了粒子滤波的估计效果，其次，后向模拟粒子平滑器本身具有优良的特性，即使在滤波效果较差的时候依然能够保持一定的收敛性和估计精度，这两点也就确保了在使用低粒子数条件下仍能维持一定的收敛性.从图 5可以看出，本文提出的PSOPF-EM算法主要包括粒子平滑和极大似然估计两个方面，极大似然估计的数据来源于粒子滤波与平滑的输出，稳定收敛的粒子滤波与平滑输出对与极大似然估计值的收敛性影响很大，在此次仿真当中，可以明显看出，本文提出的算法在粒子数量很少的情况下依然能够保持收敛且精度很高，反之，粒子平滑在使用了很少的粒子条件下仍然是收敛的，这间接验证了本文提出的粒子滤波与后向模拟粒子平滑算法的优越性；2)采用了闭环反馈的处理方式，通过E-step和M-step的反复迭代不断修正角速度辨识与目标状态估计的误差，从而提高了辨识与估计的精度.从图 6、7中还可以看出，当量测延迟概率为0.3时，基于粒子滤波的算法整体辨识与估计精度整体较带一步任意延迟量测的CKF算法高，但是，在闪烁噪声条件下，带一步任意延迟量测的CKF算法并没有发散，只是估计效果变差了，这主要是因为，本文采取的闪烁噪声是由两个零均值的高斯白噪声加权组成的，依然符合高斯分布，所以该算法才不至于发散，依然保持一定的收敛性.从图 6、7中依然可以看出本文提出的PSO-PF算法估计精度明显高于标准的粒子滤波算法和带一步任意量测延迟的粒子滤波算法，并且，带一步任意量测延迟的粒子滤波算法的估计精度远高于标准的粒子滤波，并且这种精度优势会随着延迟概率的增加越发明显，而标准的粒子滤波在量测延迟条件下会存在发散现象，并不适合该条件下的目标状态估计与角速度的辨识.

其次，从该算法的本身结构着手，对该算法进行评估分析.该算法主要包括：粒子滤波与粒子平滑算法以及极大似然估计的迭代寻优两个方面，也就是说本文提出的算法服从双收敛性，因此本文主要从粒子数量和EM算法最大迭代次数两个方面对该算法的性能进行评估分析.首先分析粒子数对算法性能的影响，粒子数分别为50、100、150、200，采取控制变量法的思想，最大迭代次数均为5次，滑动窗口均为5，各执行50次蒙特卡洛仿真如图 8所示.图 9为一次蒙特卡洛仿真转弯角速度的辨识结果，从图 9中可以看出，在蛇形机动第1阶段，当粒子数量100以上时，粒子数量对于转弯角度的辨识的收敛性以及辨识的精度影响不大，当角速度发生突变后，粒子数100以上的辨识结果也较为平稳，整体看来，当最大迭代次数为5，滞后窗口长度为5时，粒子数量维持在100以上，可以确保角速度辨识的精度以及稳定性.从图 9、10中可以发现当粒子数量低于100时，转弯角速度的辨识结果不稳定，存在发散的趋势，转弯角速度辨识与目标状态估计风险较大，图 11关于状态变量估计的均方根误差更能说明这一点.从图 11中可更加直观的看出粒子数位50时目标状态估计发散，当粒子数大于100时状态估计效果收敛且在蛇形机动第1阶段，粒子数量在目标状态估计精度方面的优势并没有体现出来，精度相差不大，当角速度发生突变后，发现粒子数量越大，状态估计的波动就越大，尤其在位置量y方向的均方根误差较为明显突出，但从整体看来，当粒子数量大于100时，发现粒子数量对于目标状态估计的精度影响不大，这是因为当粒子数量大于100时，本文提出的算法就已经能够保证辨识与估计收敛了，所以粒子数量对于精度的整体影响不大，但是从表 1还可以看出，粒子数越多，消耗的时间越长，从算法的时间成本以及估计与辨识的精度两方面综合考虑粒子数为100效果较好.

对于迭代次数而言，将窗口长度设置为5，最大迭代次数分别设置为3、5、8、10，各执行50次蒙特卡洛仿真.图 12为一次蒙特卡洛仿真转弯角速度的辨识结果，从图 12中可以看出，当迭代次数低于5时，该算法并不能保证转弯角速度辨识收敛，当迭代次数大于5时，随着迭代次数的增加，转弯角速度辨识的精度较高并且稳定，这一点在图 13中也得到了证实.从图 14还可以看出，目标状态4个分量的均方根误差随着迭代次数增加，逐渐减少，即该算法估计的目标状态整体精度就越高，尤其是当角速度发生突变之后，效果更加明显.但是从表 2显示的不同迭代次数下计算k=30 s时的角速度和状态所耗费的时间来看，随着迭代次数的增加，该算法的运算时间相应地就越长，结合图 13、14转弯角速度以及目标状态的均方根误差来看，当迭代次数大于5时，通过增加迭代次数带来的精度收益与算法运算时间成本来讲性价比不高，即当迭代次数大于5时，由迭代次数带来的精度效益远远比不上时间带来的损耗，降低该算法的实时性.

4 结论

1) 在E-step，通过充分考虑量测一步随机延迟特性及非高斯量测噪声，重构了粒子滤波算法中的粒子权重更新公式，并将拒绝采样思想引入到后向模拟粒子平滑器当中，并相应地设置拒绝采样终止条件，进一步优化了后向模拟粒子平滑器算法，提高了平滑算法的执行效率，然后采用改进的粒子滤波器与后向模拟粒子平滑器进一步获取目标状态的平滑估计.

2) 在M-step，通过采用数值优化算法极大化条件似然函数，从而获得转弯角速度的估计量，用于下一次算法迭代.通过E-step和M-step的不断迭代，进而获得转弯角速度的闭环形式的优化解，降低了辨识风险，提高了状态估计精度.

3) 仿真实验表明，在量测随机延迟及非高斯量测噪声条件下，本文提出的联合估计与辨识算法性能优于传统的扩维法.