1 引言

多普勒信息对雷达运动目标跟踪性能的改善一直备受关注，通过目标多普勒观测信息，跟踪器能够更为精确地估计航迹参数。然而雷达的多普勒盲区(Doppler Blind Zone, DBZ)使得目标可能出现连续漏检问题，致使多目标跟踪性能恶化，因而提高雷达在多普勒盲区下的多目标跟踪性能对提高多普勒雷达探测性能具有重要意义。

近年来，基于随机有限集(Random Finite Set, RFS)的多目标跟踪算法^[1,2]为多目标跟踪问题的研究提供了新的思路，因可避免复杂的数据关联，备受国内外学者推崇。该类算法已扩展到更一般的航迹形成^[3]、自适应出生强度^[4]和杂波密度^[5]、传感器空间配准^[6]和分布式融合^[7,8]等一系列问题。本文基于RFS框架下的高斯混合概率假设密度(Gaussian Mixture Probability Hypothesis Density, GM-PHD)滤波器^[9,10]对多普勒盲区下的多目标跟踪问题开展研究。

围绕多普勒盲区下的多目标跟踪问题，文献[11–13]通过将多普勒盲区建模成依赖目标状态的检测概率，将GM-PHD滤波器应用到多普勒盲区中的目标跟踪，估计的目标数量更稳定。然而，它们仅利用了与盲区有关的最小可检测速度(Minimum Detectable Velocity, MDV)信息，未充分利用多普勒量测，且未提供严格的推导和详细的实现步骤。一般地，多普勒盲区的宽度取决于MDV，因而，MDV为重要跟踪参数。为此，本文通过将考虑了MDV的检测概率模型引入GM-PHD更新式中，提出了多普勒盲区下带MDV的GM-PHD跟踪算法。

2 引入MDV信息的检测概率模型

假设第k时刻目标的状态为 ${x}_k = \left[ {{x_k}} \ \right.{y_k} \ {z_k}$ ${\left. {{{\dot x}_k} \ {{\dot y}_k} \ {{\dot z}_k}} \right]^{\rm{T}}}$ ，( $x_k,\;y_k,\;z_k$ )为位置分量，( $\dot x_k,\;\dot y_k,\;\dot z_k$ )为速度分量， ${\left[ \cdot \right]^{\rm T}}$ 表示矩阵转置操作。传感器状态为 ${x}_k^s = {\left[ {\begin{array}{*{20}{c}} {x_k^s} & {y_k^s} & {z_k^s} & {\dot x_k^s} & {\dot y_k^s} & {\dot z_k^s} \end{array}} \right]^{\rm{T}}}$ 。分别令 $\dot r_k$ 和 $\dot r_{{\rm{c}},k}$ 为多普勒量测以及该状态对应位置处的杂波多普勒^[14]，则

$\dot r_k^ \,\,{\rm{ = }}\frac{{\left( {{x_k} - x_k^s} \right)\left( {{{\dot x}_k} - \dot x_k^s} \right) + \left( {{y_k} - y_k^s} \right)\left( {{{\dot y}_k} - \dot y_k^s} \right) + \left( {{z_k} - z_k^s} \right)\left( {{{\dot y}_k} - \dot z_k^s} \right)}}{{\sqrt {{{\left( {{x_k} - x_k^s} \right)}^2} + {{\left( {{y_k} - y_k^s} \right)}^2} + {{\left( {{z_k} - z_k^s} \right)}^2}} }}$

(1)

$\dot r_{{\rm c},k} = \! - \frac{{\dot x_k^s\left( {x_k - x_k^s} \right) + \dot y_k^s\left( {y_k - y_k^s} \right) + \dot z_k^s\left( {z_k - z_k^s} \right)}}{{\sqrt {{{\left( {x_k - x_k^s} \right)}^2} + {{\left( {y_k - y_k^s} \right)}^2} + {{\left( {z_k - z_k^s} \right)}^2}} }} \ \ \ $

(2)

杂波凹口函数为目标多普勒与杂波多普勒之差n_c ，可表示为：

$n_{\rm{c}} \!=\! \dot r_k \!-\! \dot r_{{\rm{c}},k} \!=\! \frac{{\dot x_k\left( {x_k \!-\! x_k^s} \right) \!+\! \dot y_k\left( {y_k \!-\! y_k^s} \right) \!+\! \dot z_k\left( {z_k \!-\! z_k^s} \right)}}{{\sqrt {{{\left( {x_k \!-\! x_k^s} \right)}^2} \!\!+\!\! {{\left( {y_k \!-\! y_k^s} \right)}^2} \!\!+\!\! {{\left( {z_k \!-\! z_k^s} \right)}^2}} }}$

(3)

众所周知，凹口技术在抑制杂波的同时也会对低多普勒频率的运动目标的检测造成影响^[15]。检测概率 $p_{{\rm{D}},k}\left( \boldsymbol{x} \right)$ 通常是目标状态的函数，特别对地面运动目标指示雷达而言，多普勒盲区会严重影响雷达检测概率。当 $n_{\rm{c}} < {\rm{MDV}}$ 时， $p_{{\rm{D}},k}\left( \boldsymbol{x} \right)$ 趋于0；相反当 $n_{\rm{c}} \gg {\rm{MDV}}$ 时， $p_{{\rm{D}},k}\left( \boldsymbol{x} \right)$ 趋近饱和值p_D。因而，检测概率可以表示为：

$p_{{\rm D},k}\left( \boldsymbol{x} \right) \approx p_{\rm D}\left[ {1 - {{\rm e}^{ - {{\left( {n_{\rm c}\left( x \right)/{\rm{MDV}}} \right)}^2}\log 2}}} \right]$

(4)

其中，p_D为未考虑DBZ的常规检测概率。对n_c在 ${{\hat {\boldsymbol{x}}}_{k|k - 1}} = {\left[ {{{\hat x}_{k|k - 1}} \ {{\hat y}_{k|k - 1}} \ {{\hat z}_{k|k - 1}} \ {{\hat {\dot {x}}}_{k|k - 1}} \ {{\hat {\dot {y}}}_{k|k - 1}} \ {{\hat {\dot {z}}}_{k|k - 1}}} \right]^{\rm{T}}}$ 处进行1阶泰勒近似，有

$\begin{array}{l}n_{\rm{c}}\left( {{\boldsymbol{x}}_k} \right) \approx n_{\rm{c}}\left( {{\hat {\boldsymbol{x}}}_{k|k - 1}} \right){\rm{ + }}\frac{{\partial n_{\rm{c}}}}{{\partial {\boldsymbol{x}}_k}}| {_{{{x}}_k = \hat {\boldsymbol{x}}_{k|k - 1}}} \left( {{\boldsymbol{x}}_k - {\hat {\boldsymbol{x}}}_{k|k - 1}} \right)\\\;\;\;\;\;\;\;\;\;\;\; = n_{\rm{c}}\left( {{\hat {\boldsymbol{x}}}_{k|k - 1}} \right) \\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \frac{{\partial n_{\rm{c}}}}{{\partial {\boldsymbol{x}}_k}}| {_{{\boldsymbol{x}}_k = \hat {{x}}_{k|k - 1}}} {\hat {{x}}}_{k|k - 1}{\rm{ + }}\frac{{\partial n_{\rm{c}}}}{{\partial {\boldsymbol{x}}_k}}| {_{{{x}}_k = \hat {{x}}_{k|k - 1}}} {{x}}_k\\\;\;\;\;\;\;\;\;\;\;\; = {y_{{\rm{f}}}}\left( {{\hat {\boldsymbol{x}}}_{k|k - 1}} \right) - {{\boldsymbol{H}}}_{\rm{f}}\left( {{\hat {\boldsymbol{x}}}_{k|k - 1}} \right){\boldsymbol{x}}_k\end{array}$

(5)

其中， $y_{\rm{f}}\left( {{{\hat {\boldsymbol{x}}}}_{k|k - 1}} \right)$ 表示伪造量测函数，其表达式和对应的伪造量测矩阵分别为：

$y_{\rm{f}} \!=\! y_{\rm{f}}\left( {{{\hat {\boldsymbol{x}}}}_{k|k - 1}} \right) \!=\! n_{\rm{c}}\left( {{{\hat {\boldsymbol{x}}}}_{k|k - 1}} \right) \!+\! {{\boldsymbol{H}}}_{\rm{f}}\left( {{{\hat {\boldsymbol{x}}}}_{k|k - 1}} \right){{\hat {\boldsymbol{x}}}}_{k|k - 1}$

(6)

$ {\boldsymbol{H}}_{\rm{f}}\!\left( {{{\hat {\boldsymbol{x}}}}_{k|k - 1}} \right) \!\!=\!\! - \frac{{\partial n_{\rm{c}}}}{{\partial {\boldsymbol{x}}_k}}| {_{{{x}_k = {\hat {{x}}_{k|k - 1}}}} \!\!=\!\! {\left[\!\!\! {\begin{array}{*{20}{c}}{{n_1}}\!\! &\!\! {{n_2}}\!\! &\!\! {{n_3}} \!\!&\!\! {{n_4}} \!\!&\!\! {{n_5}} \!\!&\!\! {{n_6}}\end{array}}\!\!\! \right]^{\rm{T}}}} $

(7)

其中

${{n}_{1}}=-\left( {{{\hat{\dot{x}}}}_{k|k-1}}-{{{\hat{n}}}_{\text{c}}}\cos {{{\hat{a}}}_{k}}\cos {{{\hat{e}}}_{k}} \right)/{{\hat{r}}_{k}}$ , ${{n}_{2}}=-\left( {{{\hat{\dot{y}}}}_{k|k-1}}-{{{\hat{n}}}_{\text{c}}}\sin {{{\hat{a}}}_{k}}\cos {{{\hat{e}}}_{k}} \right)/{{\hat{r}}_{k}}$ , ${n_3} = - $ $\left( {{{\hat{\dot{z}}}}_{k|k-1}}-{{{\hat{n}}}_{c}}\sin {{{\hat{e}}}_{k}} \right)/{{\hat{r}}_{k}}$ , ${n_4} = - \cos {\hat a_k}\cos {\hat e_k}$ , ${n_5} = - \sin {\hat a_k}\cos {\hat e_k}$ , ${n_6} = - \sin {\hat e_k}$ , ${\hat r_k} \mathop = \limits^\Delta \sqrt {{{\left( {\hat x_{k|k - 1} - x_k^s} \right)}^2} + {{\left( {\hat y_{k|k - 1} - y_k^s} \right)}^2} + {{\left( {\hat z_{k|k - 1} - z_k^s} \right)}^2}} $ , ${\hat a_k} \mathop = \limits^\Delta {\rm{atan2}}\left[ {\left( {\hat y_{k|k - 1} - y_k^s} \right),\left( {\hat x_{k|k - 1} - x_k^s} \right)} \right]$ , ${\hat e_k} \mathop = \limits^\Delta \arctan \left[ {\left( {\hat z_{k|k - 1} - z_k^s} \right)\Big/\sqrt {{{\left( {\hat x_{k|k - 1} - x_k^s} \right)}^2} + {{\left( {\hat y_{k|k - 1} - y_k^s} \right)}^2}} } \right]$ , $ {{\hat{n}}_{c}}={{n}_{\rm{c}}}\left( {{\widehat{\mathit{\boldsymbol{x}}}}_{k|k-1}} \right)=\frac{{{{\hat{\dot{x}}}}_{k|k-1}}\left( {{{\hat{x}}}_{k|k-1}}-x_{k}^{s} \right)+{{{\hat{\dot{y}}}}_{k|k-1}}\left( {{{\hat{y}}}_{k|k-1}}-y_{k}^{s} \right)+{{{\hat{\dot{z}}}}_{k|k-1}}\left( {{{\hat{z}}}_{k|k-1}}-z_{k}^{s} \right)}{\sqrt{{{\left( {{{\hat{x}}}_{k|k-1}}-x_{k}^{s} \right)}^{2}}+{{\left( {{{\hat{y}}}_{k|k-1}}-y_{k}^{s} \right)}^{2}}+{{\left( {{{\hat{z}}}_{k|k-1}}-z_{k}^{s} \right)}^{2}}}} $ 。

将近似 $n_{\rm{c}}\left( {\boldsymbol{x}} \right)$ 代入到式(4)中可得

$p_{{\rm{D}},k}\!\left( \!{\boldsymbol{x}} \!\right) \!=\! p_{\rm{D}}\! \!\left[ {1 \!-\! {c_{\rm{f}}}{\cal N}\!\left( \!{y_{\rm{f}}\left( {{{{\hat {\boldsymbol{x}}}}_{k|k - 1}}} \!\right);\!{\boldsymbol{H}}_{\rm{f}}\!\left(\! {{{{\hat {\boldsymbol{x}}}}_{k|k - 1}}} \!\right){\boldsymbol{x}},{R_{\rm{f}}}} \right)} \right]$

(8)

其中， ${c_{\rm{f}}} = {\rm{MDV}}\sqrt {{{π}} /\log 2} $ 表示归一化因子， ${R_{\rm{f}}} = {\rm{MD}}{{\rm{V}}^2}/\left( {2\log 2} \right)$ 表示伪造量测在多普勒域中的方差。

3 多普勒盲区下的GM-PHD滤波器设计

为实现落入多普勒盲区情况下的目标稳定跟踪，本节提出了一种考虑MDV和多普勒信息的GM-PHD滤波器。下面给出该滤波器的设计思路和具体实现过程。假设PHD滤波器的后验强度 ${\upsilon _k}\left( {\boldsymbol{x}} \right)$ 传递通过如下预测和更新两步进行^[16,17]

${\upsilon _{k|k - 1}}\left( {\boldsymbol{x}} \right) = \int {{p_{{{\rm{D}},k}}}\left( {{ζ}} \right)} {f_{k|k - 1}}\left( {{\boldsymbol{x}}|{{ζ}}} \right){\upsilon _{k - 1}}\left( {{ζ}} \right){\rm d}{{ζ}} \\ \quad \quad \quad \quad \quad +\!\! \int {{\beta _{k|k - 1}}\left( {{\boldsymbol{x}}|{{ζ}}} \right)} {\upsilon _{k - 1}}\left( {{ζ}} \right){\rm d}{{ζ}} + {\gamma _k}\left( {\boldsymbol{x}} \right)$

(9)

$ {\upsilon _k}\left( {\boldsymbol{x}} \right) \!=\! \left[ {1 \!-\! {p_{{{\rm{D}},k}}}\left( {\boldsymbol{x}} \right)} \right]\!{\upsilon _{k|k - 1}}\left( {\boldsymbol{x}} \right) \\ \quad \quad \quad \ +\!\! \sum\limits_{{{z}} \in {Z}_k} {\frac{{{p_{{{\rm{D}},k}}}\left( {\boldsymbol{x}} \right)g_k\left( {{{z}}|{{x}}} \right){\upsilon _{k|k - 1}}\left( {\boldsymbol{x}} \right)}}{{{κ} _k\left( {\boldsymbol{z}} \right) \!+\!\!\displaystyle \!\int\! {{p_{{{\rm{D}},k}}}\left( {\boldsymbol{x}} \right)g_k\left( {{{z}}|{\boldsymbol{x}}} \right){\upsilon _{k|k - 1}}\left( {\boldsymbol{x}} \right){\rm d}{\boldsymbol{x}}} }}} $

(10)

其中， ${\gamma _k}\left( \cdot \right)$ 表示第k时刻出生目标RFS ${{{Γ}} _k}$ 的强度， ${\beta _{k|k - 1}}\left( {{\boldsymbol{x}}|{\boldsymbol{ζ}}} \right)$ 表示由前一状态 ${\boldsymbol{ζ}}$ 得到的衍生目标RFS ${B_{k|k - 1}}\left( {\boldsymbol{ζ}} \right)$ 的强度， $\kappa _k\left( {\boldsymbol{z}} \right)$ 表示杂波RFS $\kappa_k $ 的强度。假设雷达的位置量测和多普勒量测不相关，且位置量测函数为线性函数，则其观测似然函数 $g_k\left( {{\boldsymbol{z}}|{\boldsymbol{x}}} \right)$ 为：

$g_k\left( {{\boldsymbol{z}}|{\boldsymbol{x}}} \right) = {\cal N}\left( {{{\boldsymbol{y}}_{\rm{c}}};{\boldsymbol{H}}_{{\rm{c}},k}{\boldsymbol{x}},{\boldsymbol{R}}_{{\rm{c}},k}} \right){\cal N}\left( {{y_{\rm{d}}};h_{\rm{d}}({\boldsymbol{x}}),\sigma _{\rm{d}}^2} \right)$

(11)

其中， z 为传统量测，包括位置分量 ${\boldsymbol{y}_{\rm{c}}}$ 和多普勒分量 ${\boldsymbol{y}_{\rm{d}}}$ , ${\boldsymbol{H}}_{{\rm{d}},k}$ 为对应位置分量的观测矩阵， ${\boldsymbol{R}}_{{\rm{c}},k}$ 为对应位置分量的观测噪声协方差矩阵， $h_{\rm{d}}( \cdot )$ 是非线性的多普勒观测函数， $\sigma _{\rm{d}}$ 是多普勒量测噪声的标准差。

由式(9)和式(10)可知，检测概率对更新强度有一定程度的影响，但不影响预测强度。因此仅需关注推导考虑MDV信息后的更新强度公式。下面引述一个定理^[18]。

定理  已知合适维度的 H , R , m , P ，且假定矩阵 R 和 P 为半正定矩阵，则

${\cal N}\left( {{{z}};{\boldsymbol{Hx}},{\boldsymbol{R}}} \right){\cal N}\left( {{\boldsymbol{x}};{\boldsymbol{m}},{\boldsymbol{P}}} \right)\\ \quad \ =\! {\cal N}\!\left( \!{{{z}};{\boldsymbol{Hm}},{\boldsymbol{S}}} \right)\!{\cal N}\!\left( \!{\boldsymbol{x};{\boldsymbol{m}} \!\!+\!\! {\boldsymbol{G}}\left(\! {{{z}} \!\!-\!\! {\boldsymbol{Hm}}} \!\right),{\boldsymbol{P}}\! \!-\!\! {\boldsymbol{GS}}{{\boldsymbol{G}}^{\rm{T}}}} \right)$

(12)

式中， ${\boldsymbol{S}} = {\boldsymbol{HP}}{{\boldsymbol{H}}^{\rm{T}}} + {\boldsymbol{R}}$ , ${\boldsymbol{G}} = {\boldsymbol{P}}{{\boldsymbol{H}}^{\rm{T}}}{{\boldsymbol{S}}^{ - 1}}$ 。

假定第k 时刻的预测强度具有以下GM形式：

${\upsilon _{k|k - 1}}\left( {\boldsymbol{x}} \right) = \sum\limits_{j = 1}^{{J_{k|k - 1}}} {w_{k|k - 1}^{\left( j \right)}{\cal N}\left( {{\boldsymbol{x}};{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)},{\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}} \right)} $

(13)

其中， $w_{k|k - 1}^{\left( j \right)}$ 是预测分量的权重、 ${\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}$ 是预测分量的状态， ${\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}$ 是预测分量的协方差。将式(8)的检测概率以及式(11)的似然函数代入式(10)，运用上述定理整理后最终可得以下更新强度公式：

$\begin{array}{l}{\upsilon _k}\left( {\boldsymbol{x}} \right) = {\sum\limits_{j = 1}^{{J_{k|k}}} {w_{k|k}^{\left( j \right)}N} \left( {{\boldsymbol{x}};{\boldsymbol{m}}_{k|k}^{\left( j \right)},{\boldsymbol{P}}_{k|k}^{\left( j \right)}} \right)} \\\;\;\;\;\;\;\;\;\; \; = \sum\limits_{j = 1}^{{J_{k|k - 1}}} {w_{k|k,0}^{\left( j \right)}{\cal N}\left( {{\boldsymbol{x}};{\boldsymbol{m}}_{k|k,0}^{\left( j \right)},{\boldsymbol{P}}_{k|k,0}^{\left( j \right)}} \right)}\\ \quad \quad \quad \quad +\! \sum\limits_{{{z}} \in {{Z}}_k} \!\!{\sum\limits_{j = 1}^{{J_{k|k - 1}}}\!\! {w_{k|k}^{\left( j \right)}\left( \!{{z}} \!\right)\!{\cal N}\!\left(\! {{\boldsymbol{x}};{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( \!{{z}} \!\right),{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left(\! {{z}} \right)} \!\right)} } \\ \quad \quad \quad \quad +\sum\limits_{j = 1}^{{J_{k|k - 1}}} {w_{k|k,{\rm f}}^{\left( j \right)}{\cal N}\left( {{\boldsymbol{x}};{\boldsymbol{m}}_{k|k,{\rm f}}^{\left( j \right)},{\boldsymbol{P}}_{k|k,{\rm f}}^{\left( j \right)}} \right)} \\ \quad \quad \quad \quad + \!\!\sum\limits_{{{z}} \in {{Z}}_k}\!\!{\sum\limits_{j = 1}^{{J_{k|k - 1}}} \!\!{w_{k|k}^{\left( j \right)}\left(\! {{{{z}}_{\rm f}}}\! \right)\!{\cal N}\!\left( \!{{\boldsymbol{x}};{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{{z}}_{\rm f}}} \right),{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left(\! {{{{z}}_{\rm f}}} \!\right)} \!\right)} } \end{array}$

(14)

式中， z _f 表示增强量测。分量 $\left\{w_{k|k,0}^{\left( j \right)},\;{\boldsymbol{m}}_{k|k,0}^{\left( j \right)},\;{\boldsymbol{P}}_{k|k,0}^{\left( j \right)}\right\}$ 由常规丢失量测更新得到， $\left\{\!w_{k|k}^{\left( j \right)}(\!{{z}}\!),\;{\boldsymbol{m}}_{k|k}^{\left( j \right)}\!({{z}}\!),\;{\boldsymbol{P}}_{k|k}^{\left( j \right)}(\!{{z}}\!)\!\right\}$ 由常规量测更新得到， $\left\{w_{k|k,{\rm f}}^{\left( j \right)},\;{\boldsymbol{m}}_{k|k,{\rm f}}^{\left( j \right)},\;{\boldsymbol{P}}_{k|k,{\rm f}}^{\left( j \right)}\right\}$ 由伪量测更新得到， $\left\{w_{k|k}^{\left( j \right)}(\!{{{z}}_{\rm f}}\!),\;{\boldsymbol{m}}_{k|k}^{\left( j \right)}(\!{{{z}}_{\rm f}}\!),\;{\boldsymbol{P}}_{k|k}^{\left( j \right)}(\!{{{z}}_{\rm f}}\!)\right\}$ 由增强量测更新得到，即

$\!\!\!\!\!\!\!\!\!\!\! w_{k|k,0}^{\left( j \right)} = \left( {1 - p_{\rm{D}}} \right)w_{k|k - 1}^{\left( j \right)} \quad \quad \quad \quad $

(15)

$ {\boldsymbol{m}}_{k|k,0}^{\left( j \right)} = {\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}, {\boldsymbol{P}}_{k|k,0}^{\left( j \right)} = {\boldsymbol{P}}_{k|k - 1}^{\left( j \right)} $

(16)

通过对位置量测和多普勒量测序贯处理可得分量 $\left\{ {w_{k|k}^{\left( j \right)}\left( {{z}} \right),{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{z}} \right),{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{z}} \right)} \right\}$ ，即

$ w_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) = \frac{{p_{\rm{D}}w_{k|k - 1}^{\left( j \right)}q_k^{\left( j \right)}\left( {\boldsymbol{z}} \right)}}{{\kappa _k\left( {\boldsymbol{z}} \right) + {w_{{\rm{sum}}}}}} $

(17)

$ {\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) = {\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right) + {\boldsymbol{G}}_{d,k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)\left( {{{\boldsymbol{y}}_{\rm{d}}} - h_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right)} \right) $

(18)

$ {\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) = \left[ {{\boldsymbol{I}} - {\boldsymbol{G}}_{{\rm{d}},k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right){\boldsymbol{H}}_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right)} \right]{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right) $

(19)

式中，均值 ${\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)$ 和协方差 ${\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)$ 分别用量测 z (其包括位置 y _c和多普勒观测 y _d)更新，均值 ${\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)$ 和协方差 ${\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)$ 分别由位置量测 y _c更新，即

$ {\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right) = {\boldsymbol{m}}_{k|k - 1}^{\left( j \right)} + {\boldsymbol{G}}_{{\rm{c}},k}^{\left( j \right)}\left( {{\boldsymbol{y}}_{\rm{c}} - {\boldsymbol{H}}_{{\rm{c}},k}{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right) $

(20)

$ {\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right) = \left[ {{\boldsymbol{I}} - {\boldsymbol{G}}_{{\rm{c}},k}^{\left( j \right)}{\boldsymbol{H}}_{{\rm{c}},k}} \right]{\boldsymbol{P}}_{k|k - 1}^{\left( j \right)} $

(21)

对应位置分量的增益和新息协方差分别为：

$ {\boldsymbol{G}}_{{\rm{c}},k}^{\left( j \right)} = {\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}{\boldsymbol{H}}_{{\rm{c}},k}{\left[ {{{\mathit{\pmb{Ξ}}}}_{{\rm{c}},k}^{\left( j \right)}} \right]^{ - 1}} $

(22)

$ {{\mathit{\pmb{Ξ}}}}_{{\rm{c}},k}^{\left( j \right)} = {\boldsymbol{H}}_{{\rm{c}},k}{\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}{\left[ {{\boldsymbol{H}}_{{\rm{c}},k}} \right]^{\rm{T}}} + {\boldsymbol{R}}_{{\rm{c}},k} $

(23)

在式(18)和式(19)中的对应多普勒分量的增益 ${\boldsymbol{G}}_{{\rm{d}},k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)$ 则为：

$ {\boldsymbol{G}}_{{\rm{d}},k}^{\left( j \right)}\left( {{y_{\rm{c}}}} \right) = {\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right){\boldsymbol{H}}_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right){\left[ {\Xi _{{\rm{d}},k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right]^{ - 1}} $

(24)

其中，相应的新息协方差为：

$ {{\mathit{\pmb{Ξ}}}}_{{\rm{d}},k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_c}} \right) = {\boldsymbol{H}}_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_c}} \right)} \right){\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right){\left[ {{\boldsymbol{H}}_d\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right)} \right]^{\rm{T}}} + \sigma _{\rm{d}}^2 $

(25)

在式(19)、式(24)和式(25)中， ${\boldsymbol{H}}_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right)$ 为 $\dot r_k$ 对 ${\boldsymbol{x}}_k$ 在 ${\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)$ 处的雅克比矩阵，为：

$ {\boldsymbol{H}}_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right) = \frac{{\partial \dot r_k}}{{\partial {\boldsymbol{x}}_k}}\left| {_{x_k = m_{k|k}^{\left( j \right)}\left( {{y_{\rm{c}}}} \right)}} \right. = \left[ {\begin{array}{*{20}{c}} {h_1^{(j)}} & {h_2^{(j)}} & {h_3^{(j)}} & {h_4^{(j)}} & {h_5^{(j)}} & {h_6^{(j)}} \end{array}} \right] $

(26)

其中，

$\mathit{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{\mathit{\boldsymbol{y}}}_{\rm{c}}} \right)={{\left[ \begin{matrix} \hat{x}_{k|k}^{(j)} & \hat{y}_{k|k}^{(j)} & \hat{z}_{k|k}^{(j)} & \hat{\dot{x}}_{k|k}^{(j)} & \hat{\dot{y}}_{k|k}^{(j)} & \hat{\dot{z}}_{k|k}^{(j)} \\ \end{matrix} \right]}^{\rm{T}}}$ $h_1^{(j)} \!=\!\! \left[ {\left( {\hat{\dot{x}}_{k|k}^{(j)} \!\!-\!\! \dot x_k^s} \right) \!\!-\!\! h_{\rm{d}}\!\left(\! {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)}\! \right)\!\cos \hat a_k^{(j)}\cos \hat e_k^{(j)}} \right]\!\Big/\hat r_k^{(j)}$ $h_2^{(j)} \!= \!\!\left[ {\left( {\hat{\dot{y}}_{k|k}^{(j)}\!\! -\!\! \dot y_k^s} \right) \!\!-\!\! h_{\rm{d}}\left( \!{{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( \!{{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right)\sin \hat a_k^{(j)}\cos \hat e_k^{(j)}} \!\right]\!\Big/\hat r_k^{(j)}$ $h_3^{(j)} = \left[ {\left( {\hat{\dot{z}}_{k|k}^{(j)} \!\!-\!\! \dot z_k^s} \right) \!\!- h_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right)\sin \hat e_k^{(j)}} \right]\Big/\hat r_k^{(j)}$ $h_4^{(j)} \!\!=\!\cos \hat a_k^{(j)}\!\cos \hat e_k^{(j)}$ , $h_5^{(j)} \!\!=\! \sin \hat a_k^{(j)}\!\cos \hat e_k^{(j)}, $ $h_6^{(j)} \!\!=\! \sin \hat e_k^{(j)},$ $\hat r_k^{(j)} \!= \! \sqrt {{{\left( {\hat x_{k|k}^{(j)} \! - \!x_k^s} \right)}^2} \! \!+ \! \! {{\left( {\hat y_{k|k}^{(j)} \! - \! y_k^s} \right)}^2} \! \!+ \! \! {{\left( {\hat z_{k|k}^{(j)} - z_k^s} \right)}^2}} $ $\hat a_k^{(j)} = {\rm{atan2}}\left[ {\left( {\hat y_{k|k}^{(j)} - y_k^s} \right),\left( {\hat x_{k|k}^{(j)} - x_k^s} \right)} \right]$

$\hat e_k^{(j)} \!\!= \! \arctan \! \! \left[ \! {\left( {\hat z_{k|k}^{(j)} \! \! - \! \! z_k^s} \right) \!\Big/ \!\sqrt {{{\left( \! {\hat x_{k|k}^{(j)} \! \!- \!\! x_k^s} \!\right)}^2} \!\! \!+ \! \!{{\left( \! {\hat y_{k|k}^{(j)} \! \!- \! \! y_k^s} \!\right)}^2}} } \right]$

分量 $\left\{ {w_{k|k,{\rm{f}}}^{\left( j \right)},{\boldsymbol{m}}_{k|k,{\rm{f}}}^{\left( j \right)},{\boldsymbol{P}}_{k|k,{\rm{f}}}^{\left( j \right)}} \right\}$ 的计算过程如下：

$ w_{k|k,{\rm{f}}}^{\left( j \right)} = \frac{{p_{\rm{D}}}}{{1 - p_{\rm{D}}}}{c_f}{\cal N}\left( {y_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right);{\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right){\boldsymbol{m}}_{k|k - 1}^{\left( j \right)},\Xi _{{\rm{f}},k}^{\left( j \right)}} \right)w_{k|k,0}^{\left( j \right)} $

(27)

$ {\boldsymbol{m}}_{k|k,{\rm{f}}}^{\left( j \right)} = {\boldsymbol{m}}_{k|k - 1}^{\left( j \right)} + {\boldsymbol{K}}_{{\rm{f}},k}^{\left( j \right)}\left( {y_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right) - {\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right){\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right) $

(28)

$ {\boldsymbol{P}}_{k|k,{\rm{f}}}^{\left( j \right)} = \left[ {{\boldsymbol{I}} - {\boldsymbol{K}}_{{\rm{f}},k}^{\left( j \right)}{\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right)} \right]{\boldsymbol{P}}_{k|k - 1}^{\left( j \right)} $

(29)

其中，相应的增益 ${{K}}_{{\rm{f}},k}^{\left( j \right)}$ 为：

$ {\boldsymbol{K}}_{{\rm{f}},k}^{\left( j \right)} = {\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}{\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right){\left[ {\Xi _{{\rm{f}},k}^{\left( j \right)}} \right]^{ - 1}} $

(30)

式中，相应的新息协方差为：

$ \Xi _{{\rm{f}},k}^{\left( j \right)} = {\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right){\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}{\left[ {{\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right)} \right]^{\rm{T}}} + {R_{\rm{f}}} $

(31)

在分量 $\left\{ {w_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right\}$ 基础上，进一步对伪量测更新可得分量$\left\{ {w_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right),{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right),{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right)} \right\}$，其中

$ w_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right) = - {c_{\rm{f}}}{\cal N}\left( {y_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right);{\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right){\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),\Xi _{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right)w_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) $

(32)

$ {\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right) = {\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) + {\boldsymbol{G}}_{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)\left( {y_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right) - {\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right){\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right)$

(33)

${\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right) = \left[ {{\boldsymbol{I}} - {\boldsymbol{G}}_{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right){\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right)} \right]{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)$

(34)

其中，相应的增益 ${\boldsymbol{G}}_{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)$ 为：

$ {\boldsymbol{G}}_{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) = {\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right){\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right){\left[ {\Xi _{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right]^{ - 1}} $

(35)

式中

$ \Xi _{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) = {\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right){\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right){\left[ {{\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right)} \right]^{\rm{T}}} + {R_{\rm{f}}} $

(36)

在式(17)分母中，将杂波强度表示为：

$ \kappa _k\left( {\boldsymbol{z}} \right) = \kappa _{{\rm{c}},k}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)\kappa _{{\rm{d}},k}\left( {{{\boldsymbol{y}}_{\rm{d}}}} \right) $

(37)

其中， $\kappa _{{\rm{c}},k}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)$ 为位置分量的杂波强度， $\kappa _{{\rm{d}},k}\left( {{{\boldsymbol{y}}_{\rm{d}}}} \right)$ 为多普勒分量对应的杂波强度， ${w_{{\rm{sum}}}}$ 为：

$ {w_{{\rm{sum}}}} = p_{\rm{D}}\sum\limits_{j = 1}^{{J_{k|k - 1}}} {w_{k|k - 1}^{\left( j \right)}q_k^{\left( j \right)}\left( {\boldsymbol{z}} \right)} - {c_{\rm{f}}}p_{\rm{D}}\sum\limits_{j = 1}^{{J_{k|k - 1}}} {w_{k|k - 1}^{\left( j \right)}q_k^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right)} $

(38)

式中，

$ q_k^{\left( j \right)}\left( {\boldsymbol{z}} \right) = {\cal N}\left( {{\boldsymbol{y}}_{\rm{c}};{\boldsymbol{H}}_{{\rm{c}},k}{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)},{{\mathit{\pmb{Ξ}}}}_{{\rm{c}},k}^{\left( j \right)}} \right){\cal N}\left( {y_{\rm{d}};h_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right),{{\mathit{\pmb{Ξ}}}}_{{\rm{d}},k}^{\left( j \right)}\left( {{{\boldsymbol{y}}_{\rm{c}}}} \right)} \right) $

(39)

$ q_k^{\left( j \right)}\left( {{{\boldsymbol{z}}_f}} \right) = q_k^{\left( j \right)}\left( {\boldsymbol{z}} \right){\cal N}\left( {y_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right);{\boldsymbol{H}}_{\rm{f}}\left( {{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right){\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),\Xi _{{\rm{f}},k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right) $

(40)

值得指出的是，在量测更新后由于高斯分量数目随时间不断增长，还需要优化。为了便于描述，本节将以上所提出的考虑MDV和多普勒信息的GM-PHD滤波器简记为GM-PHD-D&MDV。当MDV=0时，c_f=0，因此，所提的多普勒盲区下GM-PHD滤波器退化为带多普勒信息的GM-PHD(GM-PHD-D)滤波器^[19]。代入检测概率到更新强度式(28)中，最终产生了数目成倍增加的高斯分量。换言之，除 ${J_{k|k - 1}}$ 个分量 $\left\{ {w_{k|k,0}^{\left( j \right)},{\boldsymbol{m}}_{k|k,0}^{\left( j \right)},{\boldsymbol{P}}_{k|k,0}^{\left( j \right)}} \right\}$ 和 $\left| {{\boldsymbol{Z}}_k} \right|{J_{k|k - 1}}$ 个分量 $\left\{ {w_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right\}$ 外，还得到了额外的相同数目的分量，包括 ${J_{k|k - 1}}$ 个分量 $\left\{ {w_{k|k,{\rm{f}}}^{\left( j \right)},\;{\boldsymbol{m}}_{k|k,{\rm{f}}}^{\left( j \right)},\;{\boldsymbol{P}}_{k|k,{\rm{f}}}^{\left( j \right)}} \right\}$ 和 $\left| {{\boldsymbol{Z}}_k} \right|{J_{k|k - 1}}$ 个分量 $\left\{ {w_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right),\;{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right),\;{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right)} \right\}$ 。因为在更新后得到了更多数目的分量，提出的滤波器设计思路能更好地解释滤波器行为。分量 $\left\{ {w_{k|k,0}^{\left( j \right)},{\boldsymbol{m}}_{k|k,0}^{\left( j \right)},{\boldsymbol{P}}_{k|k,0}^{\left( j \right)}} \right\}$ 处理传统漏检情况、 $\left\{ {w_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right),{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right)} \right\}$ 处理传统观测情况，而 $\left\{ {w_{k|k,{\rm{f}}}^{\left( j \right)},{\boldsymbol{m}}_{k|k,{\rm{f}}}^{\left( j \right)},{\boldsymbol{P}}_{k|k,{\rm{f}}}^{\left( j \right)}} \right\}$ 处理因多普勒盲区造成的目标漏检情况，$\left\{ {w_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right),\;{\boldsymbol{m}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right),\;{\boldsymbol{P}}_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right)} \right\}$处理增强量测。对于第3种情况，当目标落入多普勒盲区内，将导致相应分量权值较大，因此避免了因剪枝门限导致删除相应航迹的问题，从这个意义上讲，漏检其实具有一定的价值。对于最后的情况，相应分量的权重 $w_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right)$ 为负数，见式(32)。在实际执行中，由于负权重小于剪枝门限，这些分量将被删掉，不过鉴于目标势的均值由总权值之和得到，将这些负权值赋予相应的正权重分量，即 $w_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) = w_{k|k}^{\left( j \right)}\left( {\boldsymbol{z}} \right) + w_{k|k}^{\left( j \right)}\left( {{{\boldsymbol{z}}_{\rm{f}}}} \right)$ 。此外，式(32)表示无论状态是否落入DBZ内，每一个分量将分裂为两个。为保存DBZ内对应目标的分量，仅在目标状态落入杂波凹口时才开始分裂，即状态 $\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)},{\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}} \right)$ 满足 $n_{\rm{c}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right) \le {\rm{MDV}} + \sqrt {\Xi _{{\rm{d}},k}^{\left( j \right)}} $(其中，$\Xi _{{\rm{d}},k}^{\left( j \right)} = {\boldsymbol{H}}_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right){\boldsymbol{P}}_{k|k - 1}^{\left( j \right)}{\left[ {{\boldsymbol{H}}_{\rm{d}}\left( {{\boldsymbol{m}}_{k|k - 1}^{\left( j \right)}} \right)} \right]^{\rm{T}}} + \sigma _{\rm{d}}^{\rm{2}}$)的分量才执行分裂，反之不产生分裂。为简便，该近似滤波器被记为GM-PHD-D&MDV1。

4 仿真实验与结果分析

仿真的仿真平台是MATLAB 8.0，本节通过与传统GM-PHD-D滤波器^[11,12]和GM-PHD滤波器^[13](未带多普勒观测)和GM-PHD-D&MDV滤波器的比较来验证所提出的GM-PHD-D&MDV1滤波器的有效性。本节假设在x-y平面内的各目标均为线性高斯动态方程

$ {\boldsymbol{x}}_k = {{\boldsymbol{F}}_{k - 1}}{\boldsymbol{x}}_{k - 1} + {{\boldsymbol{v}}_{k - 1}} $

(41)

式中， ${\boldsymbol{x}}_k = {\left[ {\begin{array}{*{20}{c}} {x_k} & {y_k} & {\dot x_k} & {\dot y_k} \end{array}} \right]^{\rm{T}}}$ 为目标k时刻的状态，对应的状态转移矩阵为：

$ {{\boldsymbol{F}}_{k - 1}} = \left[ {\begin{array}{*{20}{c}} 1 & {{\tau _k}}\\ 0 & 1 \end{array}} \right] \otimes {{\boldsymbol{I}}_2} $

(42)

其中， I_n 为单位矩阵， ${\tau _k}$ =1 s为采样间隔， $ \otimes $ 为内积， ${{\boldsymbol{v}}_{k - 1}}$ 为均值为零的高斯白噪声，对应的协方差矩阵为：

$ {{\boldsymbol{Q}}_{k - 1}} = \left[ {\begin{array}{*{20}{c}} {\tau _k^4/4} & {\tau _k^3/2}\\ {\tau _k^3/2} & {\tau _k^2} \end{array}} \right] \otimes {\rm{blkdiag}}\left( {\sigma _x^2,\sigma _y^2} \right) $

(43)

式中， $\sigma _x^2$ 和 $\sigma _y^2$ 分别表示在x轴和y轴方向的加速度过程噪声方差，blkdiag表示块对角矩阵。

假设每个目标的存活概率为 ${p_{{\rm{S}},k}} = 0.99$ ，检测概率为 $p_{\rm{D}} = 0.98$ ，观测遵循式(29)的量测模型，其中， H _{c, k} = [ I ₂ 　 0₂]， ${{R}}_{{\rm c},k} = \sigma _{\rm c}^2{{{I}}_2}$ , $\sigma _{\rm{d}} = 0.5\;{\rm{m/s}}$ ，位置分量的量测噪声标准差 $\sigma _{\rm{c}} = 10 \;{\rm{ m}}$ 。在监视区域[–1000, 1000](m)×[–1000, 1000](m)内假定杂波均匀分布，即 $\kappa _{{\rm c},k}( \cdot ) = \lambda _{\rm c}V \, u( \cdot )$ , $\kappa _{{\rm d},k} = 1/2{v_{\max }}$ , $\lambda _{\rm c}$ 表示单位体积的平均杂波数目，V表示监视区域的体积， $u( \cdot )$ 为均匀密度， ${v_{\max }} = 35\;{\rm{m/s}}$ 表示传感器可检测的最大速度。

为方便比较，假定所有滤波器起始位置固定。剪枝参数为T=10^–5，合并门限为U=4，以及最大高斯分量数目为J_max=100。提取多目标的门限设置为0.5。初始目标状态分别设为 ${{x}}_{0,1} = \left[ { - 500\;{\rm{m}}} \quad \right.200\;{\rm{m}}$ ${\left. {10\;{\rm{m/s}} \ \ 0\;{\rm{m/s}}} \right]^{\rm{T}}}$ 和 ${{x}}_{0,2} = \left[ { - 500 \ {\rm{m \; - 200\;m \; 10 \ m/s}}} \right.$ ${\left. {0\;{\rm{m/s}}} \right]^{\rm{T}}}$ , ${\sigma _x} = {\sigma _y} = 5\;{\rm{m/}}{{\rm{s}}^{\rm{2}}}$ 。假定初始目标RFS服从泊松分布，相应的强度为：

$ {\gamma _k}\left( x \right) = 0.1\sum\limits_{j = 1}^2 {{\cal N}\left( {{\boldsymbol{x}};{\boldsymbol{m}}_{\gamma ,k}^{\left( j \right)},{\boldsymbol{P}}_{\gamma ,k}^{\left( j \right)}} \right)} $

(44)

式中， ${\boldsymbol{m}}_{\gamma ,k}^{\left( 1 \right)} = {\left[ { - 500\;{\rm{m}},200\;{\rm{m}},0\;{\rm{m/s}},0\;{\rm{m/s}}} \right]^{\rm{T}}}$ ${\boldsymbol{m}}_{\gamma ,k}^{\left( 2 \right)} = {\left[ { - 500\;{\rm{m}}, - 200\;{\rm{m}},0\;{\rm{m/s}},0\;{\rm{m/s}}} \right]^{\rm{T}}}$ , ${\boldsymbol{P}}_{\gamma ,k}^{\left( j \right)} = {\rm{blkdiag}}{\left( {100\;{\rm{m}},100\;{\rm{m}},25\;{\rm{m/s}},25\;{\rm{m/s}}} \right)^2}$,j=1, 2。

杂波、目标的真实航迹以及传感器位置以及两个目标的真实多普勒以及MDVs随时间的变化关系分别如图1和图2所示。可以看到，当k=50时，两目标相对传感器切向飞行，因而，此时的多普勒接近0。此外，随着MDV的增大，目标通过多普勒盲区所需时间变长，导致漏检数目增多。

当MDV=1时，不同算法的估计结果如图3所示。当目标未处于多普勒盲区时，所有滤波器都能成功跟踪两个目标，其中GM-PHD滤波器性能最差，而其他三者有着相似性能。当目标落入多普勒盲区时，由MDV造成的漏检使得所有滤波器均无法跟踪目标。然而，当目标飞出多普勒盲区后，GM-PHD-D&MDV和GM-PHD-D&MDV1都可再次对目标进行跟踪，而传统GM-PHD和GM-PHD-D滤波器则无法对目标进行再次稳定跟踪。因此，考虑MDV和多普勒信息后的新滤波器能有效解决漏检问题，保持生成目标稳定航迹。

下面利用圆位置误差概率(Circular Position Error Probability, CPEP)对滤波器的航迹丢失性能进行评估，CPEP的数学表达式为：

${\rm{CPE}}{{\rm{P}}_k}\left( r \right) = \frac{1}{{\left| {{{\boldsymbol{X}}_k}} \right|}}\sum\limits_{x \in {X_k}} {{\rho _k}\left( {{\boldsymbol{x}},r} \right)} $

(45)

其中， X_k 是真实目标状态的集合， H = [ I ₂ 　 0₂], ${\!\left\| {\cdot} \right\|_2}$ 是2范数，对于位置误差半径r, ${\rho _k}\left( {{\boldsymbol{x}},r} \right) = {\rm{Pr}}\left\{ {{{\left\| {{\boldsymbol{H\hat x}} - {\boldsymbol{Hx}}} \right\|}_2} > r{\rm{ for all }}{\boldsymbol{\hat x}} \in {{{\boldsymbol{\hat X}}}_k}} \right\}$， ${\rm{Pr}}( \cdot )$ 代表事件概率。另外，也计算最优子模式分配(Optimal Sub-Pattern Assignment, OSPA)性能。具体而言，CPEP用以评价航迹丢失性能；OSPA Card是指OSPA中势(Cardinality)估计误差。OSPA Loc是指OSPA的定位(Localization)误差。

假设CPEP半径r=20 m, OSPA中参数的阶p=2和门限c=20 m，图4–图5分别给出了利用Monte Carlo仿真得到的不同MDV值下的算法跟踪性能。从图中可以看出，引入了多普勒量测后的GM-PHD-D&MDV和GM-PHD-D&MDV1滤波器具有相似的CPEP和OSPA性能，且均优于未带多普勒量测的GM-PHD滤波器。在第2阶段，所有滤波器都不能跟踪多普勒盲区下的目标，这是因为高斯分量的权值小于状态提取门限，使得目标状态未被提取。原因是因为GM-PHD-D&MDV和GM-PHD-D&MDV1滤波器可根据式(27)获得高于剪枝门限的权重，从而避免航迹被门限删除。当多普勒盲区遮蔽目标时，传统GM-PHD-D和GM-PHD滤波器不具有有效保存对应目标分量的机制。因此当运动目标重新飞出多普勒盲区时，所提出的GM-PHD-D&MDV和GM-PHD-D&MDV1滤波器可重新跟踪上目标，而传统GM-PHD-D和GM-PHD滤波器却无法实现对目标的再跟踪。

对提出的GM-PHD-D&MDV和GM-PHD-D&MDV1滤波器，需要注意以下3点：(1)GM-PHD-D&MDV1与GM-PHD-D&MDV的性能非常类似，不过，从图6可以看出前者速度快于后者，图6中，绝对时间为每次实验中100步执行的平均时间。比如，当MDV=1时，GM-PHD-D&MDV1仅需11.05 s，而GM-PHD-D&MDV却需22.87 s。相对而言，前者耗时仅是后者的48.29%。因此，GM-PHD-D&MDV1方法可有效降低计算量，性能损失并不明显。(2)随着MDV变大，多普勒盲区遮蔽后的CPEP并未降到第1阶段的最初水平，且第1阶段与第3阶段的CPEP差距增加，这是因为，漏检越多，分量保存难度越大，从而使得盲区遮蔽后目标失跟的概率升高。(3)对于较大的MDV，伪量测更新的分量权重更易高于提取门限，于是，虚假航迹数目将增加，从而，目标进入盲区前，OSPA性能中势分量将增大。换句话说，在MDV较大时，所提滤波器的性能改善的代价是增加了虚假航迹的数目。然而，在MDV较小时，改善跟踪性能的同时并未明显增多虚假航迹。此外，在第3阶段，GM-PHD和GM-PHD-D的OSPA中的定位性能似乎优于GM-PHD-D&MDV和GM-PHD-D&MDV1。实际上，造成这一现象的原因是对应的OSPA势性能接近门限值。从而，GM-PHD-D和GM-PHD滤波器的总体OSPA性能仍然较GM-PHD-D&MDV和GM-PHD-D&MDV1的性能差。

5 结论

本文针对多普勒盲区中的多目标跟踪问题，将考虑了MDV的检测概率模型代入到标准GM-PHD更新式中，提出了目标落入多普勒盲区情况下的GM-PHD运动目标跟踪算法，在改善跟踪性能的同时显著减少了虚假航迹数量，且跟踪时间较快。