文章快速检索  
  高级检索
基于决策不确定性的多目标跟踪传感器管理
田晨1, 裴扬1,2, 侯鹏1, 赵倩1     
1. 西北工业大学 航空学院, 西安 710072;
2. 光电控制技术重点实验室, 洛阳 471023
摘要: 针对高杂波、电子干扰环境, 在量测驱动的多目标滤波框架下提出了一种基于决策不确定性的传感器管理方法。首先, 根据部分可观测马尔科夫决策过程的理论, 给出了基于Rényi信息增量的传感器管理一般方法。其次, 综合考虑决策过程的信息完整性、信息质量、信息的内涵等因素, 在量测驱动的自适应滤波框架下, 基于目标运动态势评估多目标决策不确定性水平, 并选取最大决策不确定性目标。最后, 以最大决策不确定性目标的信息增量最大化为准则进行传感器分配方案的求解。仿真实验表明所提方法能够有效抑制电子干扰、杂波对多目标跟踪及传感器分配的影响, 与基于威胁的传感器管理方法相比, 所提方法的平均最优子模式分配(OSPA)距离及平均计算时长均显著降低, 且在高杂波、电子干扰情形下具有较高的可靠性。
关键词: 传感器管理    多目标跟踪    战术重要性标绘    量测驱动    部分可观马尔科夫决策过程    
Decision uncertainty based sensor management for multi-target tracking
TIAN Chen1, PEI Yang1,2, HOU Peng1, ZHAO Qian1     
1. School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China;
2. Science and Technology on Electro-Optic Control Laboratory, Luoyang 471023, China
Abstract: For electronic countermeasures and dense clutter environments, a sensor management algorithm based on decision uncertainty using the measurement-driven multi-target filter is proposed. First, according to the theory of partially observable Markov decision process, a general sensor management approach based on Rényi divergence is presented. Meanwhile, taking into account the information integrity, information quality and information connotation in the decision-making process, we evaluate the multi-target decision uncertainty level based on the target motion situation in the measurement-driven adaptive filtering framework, subsequently selecting the maximum decision uncertainty target. Finally, the sensor allocation scheme is solved with the maximum information gain of the maximum decision uncertainty target as the criterion. The simulation results show that the proposed algorithm can effectively suppress the influence of electronic countermeasures and dense clutter on multi-target tracking and sensor management. Compared with the threat-based sensor management algorithm, the average Optimal Sub-Pattern Assignment (OSPA) distance and the average calculation time are significantly reduced. In cases of dense clutter and electronic countermeasures, the proposed algorithm has high reliability.
Keywords: sensor management    multi-target tracking    tactical significance map    measurement-driven    partially observable Markov decision process    

随着军事技术的发展, 信息化作战逐渐成为现代战争的主要形式, 这也为传感器系统的有效使用提出了更高的要求[1]。特别是在区域防空中, 要实现对空中来犯目标的有效拦截, 及早发现并稳定跟踪是关键。在综合利用多传感器对多目标进行跟踪的过程中, 为实现传感器资源的高效利用, 需要对传感器接收到的多目标信息做实时评估, 保证有限的传感器资源优先分配给更感兴趣的目标[2]。从本质上讲, 依据多目标信息进行传感器管理属于多属性决策问题[3], 需要对传感器接收到的信息进行分析、排序、评价和择优, 最终决策出更切合实际的传感器管理方案。

目前为止, 共有3类基于贝叶斯理论的传感器管理方法, 即基于任务的管理方法、基于信息论的管理方法和基于风险的管理方法[1, 4]。其中, 基于风险的管理方法重点关注由作战决策所造成的潜在损失及其发生的概率, 具有良好的实际应用价值, 已成为传感器管理领域的研究热点[5-7]。威胁等级作为重要的目标状态信息, 经常被用来量化决策风险。文献[1]为了降低空中目标威胁评估结果的不准确性和传感器辐射所带来的潜在损失, 提出了一种基于风险的多传感器管理方法。文献[8]提出了一种基于威胁度的传感器管理方法, 将目标的威胁度视为与目标状态相关的函数, 并对目标威胁度的不确定性进行最小化处理, 可应用于目标区域监视与空中交通管制。在此基础上, 文献[3]基于随机有限集的多目标滤波器提出一种基于目标威胁度评估的传感器控制策略。然而, 上述文献在进行传感器管理决策时均利用目标可观测的运动特性或目标状态滤波信息, 决策结果严重依赖于传感器对目标的量测信息, 但是仅采用目标可观测的信息却无法完全表征用于决策的主要信息属性, 例如传感器对目标量测信息的完整性和信息的质量等。在高杂波、电子对抗等实际战场环境下, 一方面, 高功率的电子干扰使得传感器对目标的量测质量变差, 量测误差变大, 甚至会出现目标量测不完整等目标“暂消”现象[9];另一方面, 大量的杂波虚警参与到滤波过程与传感器管理过程, 不仅会严重影响滤波器对目标数及多目标状态的估计, 还会造成计算效率的急剧下降, 严重影响传感器资源分配的合理性和实时性[10]。如何在高杂波、电子对抗等复杂环境下, 从获得的众多传感器信息中快速地决策出合理的传感器管理方案, 对多目标跟踪系统和传感器决策系统都提出了非常高的要求[11]

近些年来, 基于随机有限集(Random Finite Set, RFS)的理论被广泛应用到多目标跟踪领域[12], 其中具有代表性的是概率假设密度(Probability Hypothesis Density, PHD)等滤波算法, 该算法可以避免多目标跟踪中的数据关联问题, 同时能够实现对目标的检测与估计[13]。目前, 已有许多基于RFS框架的传感器管理策略[3, 14-15], 但是这些策略均是在无干扰、低杂波等较为理想的环境下展开的, 在高杂波环境下, 大量的杂波虚警参与到滤波与传感器管理过程, 导致算法效率严重退化, 从而影响传感器分配的实时性。并且在PHD滤波器的应用过程中, 通常假设新生目标的PHD是先验已知的[16]。而在实际中, 新生目标的先验信息一般并不能准确获取。虽然已有文献[17]利用传感器量测信息自适应构建目标出生PHD, 但截止到目前, 高杂波、电子对抗环境下, 结合RFS多目标滤波器进行目标状态评估的传感器决策策略, 仍然没有得到系统性的方法研究。

针对上述问题, 本文在多目标跟踪背景下提出了一种基于决策不确定性的传感器管理方法。首先, 建立了基于部分可观测马尔科夫决策过程(Partially Observable Markov Decision Process, POMDP)的传感器管理模型;其次, 利用量测驱动的自适应序贯蒙特卡罗PHD(Sequential Monte Carlo PHD, SMC-PHD)滤波器, 通过门控方法对量测的属性进行区分, 从而抑制杂波对已有目标跟踪的影响, 获得对已有目标的高精度滤波;然后综合考虑用于决策过程的信息完整性(传感器对目标的探测)、信息的质量(传感器对目标的量测精度)、信息的内涵(目标的威胁度)等因素, 建立了目标决策不确定性的评估方法, 并选取决策不确定性最大的目标为优先跟踪目标;在RFS框架下, 基于Rényi信息增量作为传感器分配的评价函数, 在信息增量最大化的准则下实现基于目标决策不确定性的传感器管理;最后通过仿真实验从多目标状态估计精度、计算效率、可靠性3个方面验证了本文算法的有效性。

值得指出的是, 考虑到JDL(Joint Directors of Labs)模型[8], 基于威胁的传感器管理方法构成了从基于JDL1级的传感器管理迈向JDL2/3级传感器管理的第一步, 而基于决策不确定性的传感器管理是JDL2/3级传感器管理的进一步深入。这一转变使得传感器管理更贴近战场的实际需求。

1 基于POMDP的传感器管理模型

传感器以雷达为例, 典型的任务场景如图 1所示。敌方多个空中目标在远距离支援干扰飞机的掩护下穿越我方防御阵地。我方部署N部传感器对敌空中M个目标进行评估, 并将获取的量测信息发送到控制中心, 控制中心根据获得的信息制定相应的传感器管理方案, 并以此控制各传感器工作[1]。由于传感器量测的不确定性、杂波、远距离支援干扰等, 控制中心接收到的信息也是不确定的。因此, 本文研究的传感器管理问题本质上是一个不确定信息下的决策问题, 可以用POMDP对该问题进行建模[15]

图 1 典型任务场景 Fig. 1 Typical task scenario
1.1 传感器量测模型

考虑目标在二维笛卡尔坐标系下的运动, 单目标状态向量可表示为$\boldsymbol{x}_{k}=\left[x_{k}, y_{k}, \dot{x}_{k}, \dot{y}_{k}\right]^{\mathrm{T}}$, 分别表示xy方向的位置、速度。在大功率压制干扰下, 由于干信比增大, 造成被干扰传感器对目标的量测误差增大、探测概率下降, 从而导致较大的跟踪误差, 甚至会出现目标暂消现象。针对此问题, 本文采用一种新的量测模型, 表示为

$ {\mathit{\boldsymbol{z}}_k} = {[{r_k},{\beta _k}]^{\rm{T}}} = \xi \cdot \mathit{\boldsymbol{h}}({\mathit{\boldsymbol{x}}_k}) + {\mathit{\boldsymbol{V}}_k} $ (1)
 

式中:zk为传感器的量测向量;rkβk分别为目标相对于传感器的距离和方位角;h(·)为量测矩阵;Vk为零均值高斯噪声, 其协方差阵为P=diag(σr2, σβ2), σr2σβ2分别为距离、方位角测量误差的方差, 与传感器输出端目标的信噪比(S/N)j有关, 可表示为[18]

$ {{\sigma _r} = \frac{{\sqrt 3 {V_{\rm{c}}}}}{{2\pi B\sqrt {2{{(S/N)}_j}} }}} $ (2)
 
$ {{\sigma _\beta } = \frac{{\sqrt 3 {\lambda _{\rm{s}}}}}{{\pi D\sqrt {2{{(S/N)}_j}} }}} $ (3)
 

式中:Vc为光速;B为传感器波形带宽;λs为传感器波长;D为传感器孔径直径。电子干扰条件下, 传感器输出端信噪比(S/N)j可表示为

$ {(S/N)_j} = \frac{1}{{J/S + 1/{R_{{\rm{SN}}}}}} $ (4)
 

式中:J/S为干扰条件下传感器输入端干信比;RSN为无干扰条件下目标的信噪比。

ξ为传感器的量测系数, 是一个随机数:

$ \xi = \{ 0,1\} $ (5)
 

P{ξ=1}=Pd, P{ξ=0}=1-Pd, Pd为传感器对目标的探测概率。通常, 探测概率Pd是虚警概率Pfa和传感器输出端目标信噪比(S/N)j的函数, 详细计算可见文献[18-19]。

1.2 传感器管理

对目标跟踪而言, 通常参与跟踪的传感器越多, 跟踪性能越好。但是, 在许多特殊的应用场合, 由于传感器数量、感知范围、使用寿命、带宽、性能等的限制, 以及多传感器在时间和空间上的配准问题, 每个时刻仅允许使用一个传感器进行跟踪, 然后将跟踪结果发送到控制中心[11]。对于控制中心而言, 其传感器选择的决策结果是不确定的, 不仅与目标的跟踪性能有关, 同时也与控制中心获得的目标信息的完整性、信息的质量以及目标信息的内涵(目标的威胁程度)等有关。在这种情况下, 当目标飞抵不同区域时, 控制中心需要根据系统对目标的需求实时地选择量测信息完整的、信息质量高的传感器对更感兴趣的目标进行优先跟踪, 即最优跟踪传感器选择问题。为此, 对每一个传感器分配方案νUk给定一个相应的评价函数$\mathcal{R}(\nu)$, 最优的控制准则是使k时刻向后H(H≥1)步的评价函数$\mathcal{R} k+H(\nu)$最大化对应的最优分配序列uk, 表示为

$ {u_k} = {\rm{arg}}\mathop {{\rm{max}}}\limits_{\nu \in {U_{k + 1:k + H}}} E[{\mathcal{R}_{k + H}}(\nu )] $ (6)
 

式中:Uk+1:k+H表示向后H步总的分配方案集合。随着传感器、目标数量及H的增大, 控制集合Uk+1:k+H的势会呈指数增长, 因此为了便于计算, 本文基于“近视”(“Myopic”)方案进行传感器选择的研究, 且假定每个时刻仅使用一个传感器对所有目标进行跟踪。

为了从整体上评估跟踪系统的性能, 本文采用基于多目标整体信息增益的传感器选择方案, 即以多目标信息增益最大化作为评价函数。Rényi信息增量通过比较概率密度函数的近似程度来表示当前状态下信息的差异, 可用于强调某个局部信息, 且对多目标先验、后验概率密度函数的分布没有高斯限制, 更具灵活性和普适性[20]。因此, 本文采用Rényi信息增量作为评价函数。在RFS框架下, 假设k-1时刻多目标的先验概率密度函数和后验概率密度函数分别为pk|k-1(Xk|Z1:k-1)和pk|k(Xk|Z1:k-1, Zk(ν)), 则Rényi信息增量可表示为

$ \begin{array}{l} \mathcal{R}(\nu ) = \\ {\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} 1/\alpha - \lg \int {{{\left[ {{p_{k\mid k}}\left( {{X_k}\mid {Z_{1:k - 1}},{Z_k}(\nu )} \right)} \right]}^\alpha } \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} {[{p_{k|k - 1}}({X_k}|{Z_{1;k - 1}})]^{1 - \alpha }}{\rm{d}}{X_k} \end{array} $ (7)
 

式中:α为调整2个概率密度函数尾部重合程度的参数, 0 < α < 1。α=0.5强调概率密度函数的尾部, 可以对2个相似的概率分布达到最佳辨别, 尤其是闪烁噪声对概率密度函数造成的拖尾效应[20]。因此在本文中, 选用α=0.5。

2 自适应SMC-PHD滤波器

典型的PHD滤波器预测方程和更新方程分别为[13]

$ {D_{k|k - 1}}(\mathit{\boldsymbol{x}}) = \langle {P_{\rm{s}}}{f_{k|k - 1}}(\mathit{\boldsymbol{x}}| \cdot ),{D_{k - 1|k - 1}}\rangle + {\gamma _k}(\mathit{\boldsymbol{x}}) $ (8)
 
$ \begin{array}{l} {D_{k|k}}(\mathit{\boldsymbol{x}}) = (1 - {P_{{\rm{d}},k}}(\mathit{\boldsymbol{x}})){D_{k|k - 1}}(\mathit{\boldsymbol{x}}) + \\ {\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} \sum\limits_{z \in {\mathit{\boldsymbol{Z}}_k}} {\frac{{{P_{{\rm{d}},k}}(\mathit{\boldsymbol{x}}){g_k}(\mathit{\boldsymbol{z}}|\mathit{\boldsymbol{x}}){D_{k|k - 1}}(\mathit{\boldsymbol{x}})}}{{\langle {P_{{\rm{d}},k}}(\mathit{\boldsymbol{x}}){g_k}(\mathit{\boldsymbol{z}}| \cdot ),{D_{k|k - 1}}\rangle + {\kappa _k}(\mathit{\boldsymbol{z}})}}} \end{array} $ (9)
 

式中:Ps为单目标存活概率;fk|k-1(x|·)为单目标状态转移函数;γk(x)为k时刻新生目标RFS的PHD;gk(z|·)为单目标观测似然函数, 〈m, n〉=∫n(x)m(x)dxPd, k(x)为单目标探测概率;κk(z)为k时刻杂波RFS的PHD。

在上述PHD滤波器中, 量测集Zk是唯一的输入变量, 新生目标PHD通常假设是已知的, 而在大多数实际问题中该先验信息并不能获取。此时, 量测集Zk中与新生目标相关的量测对于已有目标的更新是冗余的[16]。特别是在高杂波、电子对抗环境下, 新生目标量测中含有大量的杂波虚警, 会严重影响滤波器对目标数及多目标状态的估计, 不仅会造成传感器资源的浪费, 更会导致滤波算法性能及效率的下降, 甚至失效。针对高杂波、电子对抗环境, 本文采用自适应的PHD滤波器对不同属性的量测进行区分并独立滤波, 以降低新生目标量测及杂波虚警对于已有目标滤波的影响, 同时显著提高算法的效率。一般而言, 杂波虚警与新生目标量测并不能直接区分, 因此本文在最优传感器选择时, 仅考虑已有目标, 对于新生目标, 需要在量测被确认后才考虑最优传感器的分配, 从而抑制杂波对传感器管理决策的影响, 避免将有限的传感器资源浪费在大量的杂波虚警上。

假设已有目标和新生目标生成的量测是可分的, 则Zk可表示为

$ {\mathit{\boldsymbol{Z}}_k} = {\mathit{\boldsymbol{Z}}_{k,{\rm{b}}}} \cup {\mathit{\boldsymbol{Z}}_{k,{\rm{e}}}} \cup {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k} $ (10)
 

式中:Zk, bZk, eΓk分别表示由新生目标、已有目标和杂波虚警生成的量测RFS。同时, 用变量ε来标记已有目标和新生目标的PHD[16], 即

$ {D_{k|k}}( \cdot ,\varepsilon ) = \left\{ {\begin{array}{*{20}{l}} {{D_{k|k}}( \cdot ,0)}&{\mathit{\boldsymbol{z}} \in {\mathit{\boldsymbol{Z}}_{k,{\rm{e}}}}}\\ {{D_{k|k}}( \cdot ,1)}&{\mathit{\boldsymbol{z}} \in {\mathit{\boldsymbol{Z}}_{k,{\rm{b}}}}} \end{array}} \right. $ (11)
 

为了解决PHD滤波方程中的积分运算问题, 本文利用SMC-PHD滤波器处理密集杂波环境下的多目标跟踪问题, 并基于门控方法实现已有目标和新生目标量测的区分。

2.1 已有目标滤波

在PHD滤波器中, k-1时刻的新生目标PHD和已有目标的PHD组成该时刻全部目标的PHD, 假设可用粒子集{wk-1, e(n), xk-1, e(n)}n=1Nk-1表示, Nk-1表示k-1时刻近似目标PHD的总粒子数。则在滤波器预测步骤, 目标的预测PHD可表示为

$ {D_{k|k - 1}}(\mathit{\boldsymbol{x}},0) \approx \sum\limits_{n = 1}^{{N_{k - 1}}} {w_{k|k - 1,{\rm{e}}}^{(n)}} {\delta _{\mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}}}(\mathit{\boldsymbol{x}}) $ (12)
 

式中:δxk, e(n)(x)表示Dirac delta函数, 对应的粒子和权值为

$ \mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)} \backsim {q_k}( \cdot |\mathit{\boldsymbol{x}}_{k - 1,{\rm{e}}}^{(n)},{\mathit{\boldsymbol{Z}}_k})\quad n = 1,2, \cdots ,{N_{k - 1}} $ (13)
 
$ \begin{array}{l} w_{k|k - 1,{\rm{e}}}^{(n)} = \frac{{{P_{\rm{s}}}(\mathit{\boldsymbol{x}}_{k - 1,{\rm{e}}}^{(n)}){f_{k|k - 1}}(\mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}|\mathit{\boldsymbol{x}}_{k - 1,{\rm{e}}}^{(n)})}}{{{q_k}(\mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}|\mathit{\boldsymbol{x}}_{k - 1,{\rm{e}}}^{(n)},{\mathit{\boldsymbol{Z}}_k})}}w_{k - 1,{\rm{e}}}^{(n)}\\ {\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} n = 1,2, \cdots ,{N_{k - 1}} \end{array} $ (14)
 

式中: $q_{k}\left(\cdot \mid \boldsymbol{x}_{k-1, \mathrm{e}}^{(n)}, \boldsymbol{Z}_{k}\right)$为重要性密度函数, 且一般取$q_{k}\left(\cdot \mid \boldsymbol{x}_{k-1, \mathrm{e}}^{(n)}, \boldsymbol{Z}_{k}\right)=f_{k \mid k-1}\left(\cdot \mid \boldsymbol{x}_{k-1, \mathrm{e}}^{(n)}\right)$

对于已有目标, 其更新的PHD可表示为

$ {D_{k|k}}(\mathit{\boldsymbol{x}},0) \approx \sum\limits_{n = 1}^{{N_{k - 1}}} {\mathit{\boldsymbol{w}}_{k|k,{\rm{e}}}^{(n)}} {\delta _{\mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}}}(\mathit{\boldsymbol{x}}) $ (15)
 

式中:

$ \begin{array}{l} w_{k|k,{\rm{e}}}^{(n)} = (1 - {P_{{\rm{d}},k}}(\mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}))w_{k|k - 1,{\rm{e}}}^{(n)} + \\ {\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} \sum\limits_{z \in {\mathit{\boldsymbol{Z}}_{k,{\rm{e}}}}} {\frac{{{P_{{\rm{d}},k}}\left( {\mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}} \right){g_k}\left( {\mathit{\boldsymbol{z}}\mid \mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}} \right)w_{k\mid k - 1,{\rm{e}}}^{(n)}}}{{{\kappa _{k,{\rm{e}}}}(\mathit{\boldsymbol{z}}) + \sum\limits_{n = 1}^{N_{k - 1}^{\rm{e}}} {{P_{{\rm{d}},k}}} \left( {\mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}} \right){g_k}\left( {\mathit{\boldsymbol{z}}\mid \mathit{\boldsymbol{x}}_{k,{\rm{e}}}^{(n)}} \right)w_{k\mid k - 1,{\rm{e}}}^{(n)}}}} \end{array} $ (16)
 

得到加权粒子集$\left\{{w}_{k \mid k, \mathrm{e}}^{(n)} \boldsymbol{, } \boldsymbol{x}_{k, \mathrm{e}}^{(n)}\right\}_{n=1}^{N_{k-1}}$后, 已有目标数可表示为

$ {\hat N_{k,{\rm{e}}}} = \sum\limits_{n = 1}^{{N_{k - 1}}} {w_{k|k,{\rm{e}}}^{(n)}} $ (17)
 

对于SMC-PHD滤波器而言, 为了缓解粒子退化问题, 通常需要重采样步骤消除具有较低更新权值的粒子。重采样的粒子数Nk依据当前估计的目标数和分配给每个目标的粒子数η来确定[16]

$ {N_k} = [\eta \cdot {\dot N_{k,{\rm{e}}}}] $ (18)
 

式中:[·]表示取整运算。经过重采样后, 可用新的加权粒子集$\left\{w_{k}^{(n)}, \boldsymbol{x}_{k}^{(n)}\right\}_{n=1}^{N_{k}}$近似表示Dk|k(x, 0)。

一般地, 粒子权值越低, 其所代表的真实目标信息就越少, 特别是在高杂波、电子对抗等复杂环境下。重采样时消除较低权值的粒子, 可以避免大量与目标信息无关的粒子直接参与滤波过程, 提高了滤波器的有效性。尽管可能会损失极少的目标信息, 但是在本文采用的自适应PHD滤波框架下, 这样的目标信息会在下一次迭代中以新生目标的属性继续出现。因此重采样步骤中, 消除较低权值的粒子并不会直接影响目标有效信息的完整性。

2.2 新生目标滤波

对于新生目标PHD, 本文采用量测驱动的自适应方法, 这意味着目标已经被传感器探测到, 即Pd, k=1在此处滤波框架下一直适用于新生目标。假设量测值zk, biZk, b(i=1, 2, …, mk, b)已知, mk, bZk, b中量测的数目, 则新生目标可以通过在似然函数gk(zk, bi|x(i))非零区域进行状态粒子采样得到[17]

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{x}}_{k,{\rm{b}}}^{(i,n)} \backsim N(\mathit{\boldsymbol{x}};{\mathit{\boldsymbol{h}}^{ - 1}}(\mathit{\boldsymbol{z}}_{k,{\rm{b}}}^i),{\mathit{\boldsymbol{H}}_*}\mathit{\boldsymbol{RH}}_*^{\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} {\kern 1pt} n = 1,2, \cdots ,{M_{\rm{b}}}} \end{array} $ (19)
 

式中:h-1(·)为h(·)的逆运算;H*h-1(·)的雅克比矩阵;Mb表示从每个量测采样的粒子数。每一个zk, bi都被视为一个可能源自新生目标的量测。显然, 此时的新生目标中含有大量的杂波虚警。多个潜在目标的PHD在单目标状态空间是可加的, 即

$ {D_{k|k}}(\mathit{\boldsymbol{x}},1) = \sum\limits_{i = 1}^{{m_{k,{\rm{b}}}}} {D_{k|k}^i} (\mathit{\boldsymbol{x}},1) $ (20)
 

式中:Dk|ki(x, 1)为单目标点集{xk, i}的PHD, 可近似表示为

$ D_{k|k}^i(\mathit{\boldsymbol{x}},1) \approx \sum\limits_{n = 1}^{{M_{\rm{b}}}} {w_{k|k,{\rm{b}}}^{(i,n)}} {\delta _{\mathit{\boldsymbol{x}}_{k,{\rm{b}}}^{(i,n)}}}(\mathit{\boldsymbol{x}}) $ (21)
 

其中, 每个粒子的初始化权值为

$ w_{k|k,{\rm{b}}}^{(i,n)} = \frac{1}{{{M_{\rm{b}}}}} $ (22)
 

为了较好地初始化新生目标, Mb的值需要依据实际情况选取。当杂波信息未知时, 引入平衡因子μ=1/mk, b, (mk, b>0)来平衡杂波量测的影响[15], 其含义是依据可能的新生目标量测数目自适应调整每个量测都源于新生目标的信任度, 可以随着mk, b的变化而变化。此时, wk|k, b(i, n)=μ/Mb

当采样和权值分配对Zk, b中所有的量测都处理后, 则Dk|k(x, 1)的加权粒子集近似{wk, b(n), xk, b(n)}Nkbn=1可表示为

$ \begin{array}{l} \{ w_{k,{\rm{b}}}^{(n)},\mathit{\boldsymbol{x}}_{k,{\rm{b}}}^{(n)}\} _{n = 1}^{N_k^{\rm{b}}} = \{ \omega _{k|k,{\rm{b}}}^{(1,n)},\mathit{\boldsymbol{x}}_{k,{\rm{b}}}^{(1,n)}\} _{n = 1}^{{M_{\rm{b}}}} \cup \\ {\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} \{ w_{k|k,{\rm{b}}}^{(2,n)},x_{k,{\rm{b}}}^{(2,n)}\} _{n = 1}^{{M_{\rm{b}}}} \cup \cdots \cup \{ w_{k|k,{\rm{b}}}^{({m_{k,{\rm{b}}}},n)},x_{k,{\rm{b}}}^{({m_{k,{\rm{b}}}},n)}\} _{n = 1}^{{M_{\rm{b}}}} \end{array} $ (23)
 

式中:Nkb=Mbmk, b表示当前与新生目标对应的粒子总数。

2.3 门控方法

与文献[3]所采用的的标准PHD滤波器相比, 上述量测驱动自适应滤波器对已有目标和新生目标独立滤波, 由于排除了杂波等无关量测对已有目标更新的影响, 可以显著改善已有目标的滤波准确性。同时, 独立滤波便于对不同属性的目标量测分配传感器, 在最优传感器选择时, 仅考虑已有目标, 对于含有大量杂波虚警的新生目标, 并不在第一时间考虑, 从而降低了杂波等对目标跟踪及传感器管理决策的影响。在滤波过程中, 由于含有大量杂波虚警的新生目标粒子并未参与到SMC-PHD滤波器的权值更新过程中, 而且基于量测的新生目标采样处理具有较低的计算复杂度, 从而可以显著提高滤波算法的效率, 保证了传感器管理决策的实时性。文献[10]提出了一种基于统计距离的验证窗口以区分新生目标量测与已有目标量测, 在此基础上, 本文根据新息加权平方和在传感器量测空间Es上定义一个验证区域Ω(z), 其半径为

$ {r_{\rm{d}}} = {\rm{max}}\{ \sqrt {{d^2}(\mathit{\boldsymbol{z}},\mathit{\boldsymbol{h}}(\mathit{\boldsymbol{x}}))} :\mathit{\boldsymbol{x}} \in {E_{\rm{s}}}(\mathit{\boldsymbol{x}})\} $ (24)
 

式中:d2(z, h(x))=(z-h(x))TSk-1(z-h(x))为新息加权平方和。假设目标的预测粒子集{xk, e(n), wk|k-1, e(n)}n=1Nk-1, 中间变量集XCWC分别用于存储和单个候选量测zk, j关联的粒子和对应的权值, 用来综合评估每个量测所关联的粒子数目和粒子权值和[16]。定义阈值参数τγ, 从而确定当前量测的属性。该门控方法的伪代码如表 1所示。

表 1 门控方法伪代码 Table 1 Pseudo code of gating method
门控方法伪代码
步骤1 初始化Zk, b=0, Zk, e=0
步骤2 对zk, jZk, 确定与其关联粒子和相应的权值, 计算验证区域半径rd, 初始化XC=Ø, WC
  for each xk, e(n), n=1, 2, …, Nk-1 do
    if d(zk, j, h(xk, e(n)))≤rd do
      XC=XC∪{xk, e(n)}, WC=WC∪{wk|k-1, e(n)}
    end if
  end for
步骤3 计算XC中粒子数和WC中的权值和
  nc=|XC|
  wc=∑wk|k-1, e(i), wk|k-1, e(i)WC and i=1, 2, …, nc
步骤4 确定量测属性
  if ncτ and wcγ do
    Zk, e=Zk, e∪{zk, j}
  else
    Zk, b=Zk, b∪{zk, j}
  end if
步骤5 转至步骤2, 直到处理完所有的量测

上述过程利用了粒子在空间中的分布信息和相应的权值信息, 可以显著提高量测区分的准确性[16]。在滤波初始时刻, 所有的量测都被认为源于新生目标。阈值参数τγ会对量测属性有一定的影响。较小的τγ, 可能会将新生目标量测误认为已有目标量测, 在这些新生目标量测中可能存在杂波虚警等, 从而造成滤波精度与滤波效率下降, 甚至传感器资源的严重浪费。较大的τγ, 可能会导致将已有目标量测误认为新生目标量测, 在该时刻最大决策不确定性目标选择时将不会再考虑该目标, 进而可能影响最优传感器的方案选择, 但是并不会影响算法的有效性。因此, 在实际中可以保守地将阈值τγ设置为一个相对较大的值。

在上述SMC-PHD滤波过程中, 新生目标和已有目标的滤波相互独立, 通过门控方法可以尽可能地排除杂波等无关量测参与到已有目标更新过程中, 从而显著提高杂波环境下已有目标的滤波准确性及算法的效率。虽然杂波量测也可能会被误用于生成新生目标粒子, 但是这样的粒子绝大多数会在下一次的滤波迭代过程中被自然消除, 而真实的新生目标会结合最新的量测被验证。因此, 本文在最优传感器选择时, 仅考虑已有目标, 对于新生目标, 需要下一次迭代确认后才考虑。这样既提高了对已有目标的跟踪精度, 又提高了传感器管理算法的效率, 同时可以避免将有限的传感器资源浪费在大量的杂波虚警上。

3 基于目标决策不确定性的传感器管理

对于控制中心而言, 传感器分配的决策结果是不确定的, 受很多因素的影响, 特别是在高杂波、电子对抗等复杂环境下, 高功率的电子干扰使得传感器对目标的量测质量变差, 量测误差变大, 甚至会出现目标量测不完整等目标暂消现象。因此, 需要对控制中心接收到的目标信息进行分析、评价和择优, 以最大程度的降低监视区域目标信息的不确定性, 最终决策出更切合实际的传感器分配方案。本文仅从传感器对目标的探测和跟踪角度出发, 选取易于量化且对决策结果最具影响的因素:信息的完整性、信息的质量、信息的内涵等, 对决策的不确定性进行描述。信息的完整性表征传感器对目标状态的量测是否完整, 杂波、电子干扰等会导致传感器对目标的探测概率下降, 甚至出现目标暂消现象, 加大了决策的不确定性。信息的质量表示传感器对目标量测的精度, 量测精度越高, 滤波结果越可信, 信息的质量就越高, 决策的不确定性就越低。信息的内涵表示传感器对目标的滤波信息中所蕴含的目标的威胁度信息, 目标威胁越高, 决策不确定性就越高。在本文中, 通过将探测概率引入传感器量测方程, 用以模拟传感器对目标量测的完整性。同时, 选取3种常见的目标运动特性(目标速度、航向及距离)对目标威胁度进行建模。并基于自适应的SMC-PHD滤波器估计多目标状态, 分别考虑目标信息的完整性、质量、内涵等因素, 依据多目标运动态势快速评估多目标决策不确定性水平, 并从中选择出当前时刻最大决策不确定性目标, 最终决策得到评价准则最优的传感器分配方案。

3.1 已有目标预测PHD及状态提取

在自适应SMC-PHD滤波器中, k-1时刻新生目标和已有目标的PHD组成该时刻全部目标的PHD。然而此时的新生目标中含有大量的杂波虚警, 需要结合k时刻最新的量测信息进行确认和估计。因此, 为了抑制杂波虚警对传感器分配决策的影响, 在最优传感器选择时, 仅考虑k-1时刻的已有目标, 暂时不考虑新生目标。

假设k-1时刻的已有目标的后验PHD可由一组带权值的随机样本粒子集$\left\{ {\mathit{w}_{k - 1, {\rm{e}}}^{(n)}, \mathit{\boldsymbol{x}}_{k - 1, {\rm{e}}}^{(n)}} \right\}_{n = 1}^{N_{k - 1}^{\rm{e}}}$表示, Nk-1e表示k-1时刻已有目标的粒子总数。依据SMC-PHD滤波算法, 用一组带权值的随机样本粒子集$\left\{ {w_{k\mid k - 1}^{(n)}, \mathit{\boldsymbol{\tilde x}}_k^{(n)}} \right\}_{n = 1}^{N_{k - 1}^{\rm{e}}}$来表示已有目标的预测PHD函数, 即

$ {D_{k|k - 1}}(\mathit{\boldsymbol{x}}) \approx \sum\limits_{n = 1}^{N_{k - 1}^{\rm{e}}} {w_{k|k - 1,{\rm{e}}}^{(n)}} {\delta _{\mathit{\boldsymbol{\tilde x}}_k^{(n)}}}(\mathit{\boldsymbol{x}}) $ (25)
 

通常, 目标状态集合可以通过加权粒子集利用峰值提取技术得到。常见的方法包括k-means[21]、有限混合模型[22]、CLEAN[23]等。但上述方法大多需要复杂的迭代计算, 影响算法的实时性、可靠性。因此本文选择文献[24]中的MEAP(Multiple Expected A Posteriori)方法解决多目标状态提取问题, 得到目标状态集合${{\mathrm{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\boldsymbol{X}}}}_{k\mid k-1}}$, 利用该集合生成分配方案所需的量测集合。

3.2 决策不确定性的确立及分配方案的决策

本文分别从信息的完整性(探测概率)、信息的质量(测量精度)、信息的内涵(目标威胁度)3个方面描述控制中心决策的不确定性。不同影响因素的综合评估通常采用线性加权方法[8], 需要根据不同的应用场景对权值进行赋值, 环境适应性较差。为了解决线性加权法不能体现各因素间的非线性关系问题, 本文基于战术重要性标绘(Tactical Significance Map, TSM)方法[25], 提出了一种改进的TSM(Improved TSM, ITSM)方法, 依据提取的各目标预测状态${{\mathrm{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\boldsymbol{x}}}}_{k\mid k-1}}\in {{\mathrm{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\boldsymbol{X}}}}_{k\mid k-1}}$来确定当前时刻各目标的决策不确定性。

假设当前时刻传感器对目标i的探测概率为Pd(i), 测量误差方差分别为σr, i2σβ, i2。目标的状态矢量可表示为$\boldsymbol{x}_{i}=\left[\boldsymbol{p}_{i}, \dot{\boldsymbol{p}}_{i}\right]^{\mathrm{T}}=\left[x_{i}, y_{i}, \dot{x}_{i}, \right.\left.\dot{y}_{i}\right]^{\mathrm{T}} $, 其中pi$\dot{\boldsymbol{p}}_{i}$分别表示目标i的位置矢量和速度矢量。假设我方重要战略防御区域中心的位置为x0=[px, 0, py, 0]T, 则目标与战略防御区域中心的距离为$d_{i}=\sqrt{\left(x_{i}-p_{x, 0}\right)^{2}+\left(y_{i}-p_{y, 0}\right)^{2}}$。若仅考虑目标的战术重要性(威胁度)时的TSM函数表达式为[25]

$ {f_{{\rm{TSM}}}}({\mathit{\boldsymbol{x}}_i}) \backsim {\rm{exp}}\left( { - \frac{{d_i^2}}{{2\sigma _{{\rm{TSM}}}^2}}} \right) $ (26)
 

式(26)表明, 目标距离防御区域中心越近, 则其威胁度越高, 反之亦然。其中, σTSM取决于目标的速度和航向。若目标以较高的速度朝向防御区域中心运动, 则TSM函数值应该很高;若该目标朝防御区域中心反方向运动, 则较高的速度应该对应更低的TSM函数值。因此可将目标速度和航向相互影响的σTSM描述为[25]

$ {\sigma _{{\rm{TSM}}}}({\mathit{\boldsymbol{p}}_i},{\mathit{\boldsymbol{\dot p}}_i}) = \left( {1 - \frac{{\theta ({\mathit{\boldsymbol{p}}_i},{{\mathit{\boldsymbol{\dot p}}}_i})}}{\pi }} \right)({k_0}\left\| {{{\mathit{\boldsymbol{\dot p}}}_i}} \right\| + {m_0}) $ (27)
 

式中:k0m0为正常数; θ(pi, $\dot{\boldsymbol{p}}_{i}$)为目标i的航向角, 即目标位置矢量和速度矢量的夹角,

$ \theta ({\mathit{\boldsymbol{p}}_i},{\mathit{\boldsymbol{\dot p}}_i}) = {\rm{arccos}}\frac{{({\mathit{\boldsymbol{p}}_i},{{\mathit{\boldsymbol{\dot p}}}_i})}}{{\left\| {{\mathit{\boldsymbol{p}}_i}} \right\|\left\| {{{\mathit{\boldsymbol{\dot p}}}_i}} \right\|}} $ (28)
 

式中:($\boldsymbol{p}_{i}, \dot{\boldsymbol{p}}_{i}$)表示标量积。

综上可得仅考虑目标的战术重要性时TSM函数的最终表达式为[25]

$ \begin{array}{l} {f_{{\rm{ TSM }}}}({x_i}) = \\ {\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} {\rm{exp}}\left( { - \frac{{d_i^2}}{{2{{\left( {1 - \frac{{\theta ({\mathit{\boldsymbol{p}}_i},{{\mathit{\boldsymbol{\dot p}}}_i})}}{\pi }} \right)}^2}{{({k_0}\left\| {{{\mathit{\boldsymbol{\dot p}}}_i}} \right\| + {m_0})}^2}}}} \right) \end{array} $ (29)
 

对于控制中心而言, 仅考虑目标的战术重要性显然无法完全表征用于决策的信息属性。通常, 目标可观测的运动特性依赖于传感器对目标的探测概率和测量精度。高杂波、电子对抗等会显著影响传感器对目标量测信息的完整性及信息的质量, 这不仅会影响滤波器的精度及效率, 也会严重影响传感器管理的决策结果。因此, 本文提出了一种综合考虑信息完整性、信息质量以及信息内涵的ITSM方法, 表达式为

$ \begin{array}{l} {f_{{\rm{ ITSM }}}}({\mathit{\boldsymbol{x}}_i}) = \\ {\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} {\rm{exp}}\left( { - \frac{{d_i^2}}{{2{{(1 - \frac{{\theta ({\mathit{\boldsymbol{p}}_i},{{\mathit{\boldsymbol{\dot p}}}_i})}}{\pi })}^2}{{({k_0}\left\| {{{\mathit{\boldsymbol{\dot p}}}_i}} \right\| + {m_0})}^2}}}} \right)\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} {\rm{exp}}\left( { - \frac{{P_{\rm{d}}^{(i)}}}{{{k_1}\sigma _r^2 + {k_2}\sigma _\beta ^2}}} \right) \end{array} $ (30)
 

式中:k1k2为正常数, 且k1+k2=1。显然, 传感器对目标的探测概率越大, 表明信息越完整, 目标的决策不确定性就越小;传感器对目标的测量误差方差越大, 表明信息质量越差, 目标的决策不确定性就越大。

依据式(30), 可在预测多目标状态集合${{\mathrm{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\boldsymbol{X}}}}_{k\mid k-1}}$中确定当前时刻最大决策不确定性目标${{\mathrm{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\boldsymbol{x}}}}_{\text{de}, k\left| k-1 \right.}}$, 利用目标状态提取过程, 提取当前时刻最大决策不确定性目标所对应的粒子集合。在局部状态空间上, 目标的粒子集可加权近似最大决策不确定性目标的分布特性[3]。假设最大决策不确定性目标的粒子集为$\left\{ {w_{{\rm{de}}, k\mid k - 1}^{(i)}, \mathit{\boldsymbol{x}}_{{\rm{de}}, k\mid k - 1}^{(i)}} \right\}_{i = 1}^{{N_{\rm{T}}}}$, NT表示最大决策不确定性目标对应的粒子总数。则最大决策不确定性目标的预测PHD可近似表示为

$ {D_{{\rm{de}},k|k - 1}}(\mathit{\boldsymbol{x}}) \approx \sum\limits_{i = 1}^{{N_{\rm{T}}}} {w_{{\rm{de}},k|k - 1}^{(i)}} {\delta _{\mathit{\boldsymbol{x}}_{{\rm{de}},k|k - 1}^{(i)}}}(\mathit{\boldsymbol{x}}) $ (31)
 

本文以最大决策不确定性目标的Rényi信息增量作为评价指标, 目的是使得最大决策不确定性目标的跟踪性能达到最优, 从而降低监视区域内目标的决策不确定性。在本文ITSM函数的计算中, 已充分考虑了量测不准确、密集杂波、电子干扰等对决策结果的影响, 因此, 此处可在检测概率Pd, k(x)=1以及不考虑量测噪声、杂波和电子干扰等理想情况下, 对每个目标仅生成一个相应的量测。最大决策不确定性目标的理想量测集可表示为

$ {\mathit{\boldsymbol{z}}_k} = \mathit{\boldsymbol{h}}({\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over x} }}_{{\rm{de}},k|k - 1}}) $ (32)
 

在执行传感器分配方案时, 按照式(32)所给定的最大决策不确定性目标的理想量测集去对已有目标的预测PHD的粒子集进行更新, 则最大决策不确定性目标的更新PHD可近似表示为

$ {D_{{\rm{de}},k|k}}(\mathit{\boldsymbol{x}}) \approx \sum\limits_{i = 1}^{{N_{\rm{T}}}} {w_{{\rm{de}},k|k}^{(i)}} {\delta _{\mathit{\boldsymbol{x}}_{{\rm{de}},k|k - 1}^{(i)}}}(\mathit{\boldsymbol{x}}) $ (33)
 

式中:

$ \tau _{{\rm{de}},k|k}^{(i)} \approx w_{{\rm{de}},k|k - 1}^{(i)} \cdot \sum\limits_{z \in {\mathit{\boldsymbol{Z}}_k}} {\frac{{{g_k}(z|x_{{\rm{de}},k|k - 1}^{(i)})}}{{\sum\limits_{i = 1}^{{N_{\rm{T}}}} {w_{{\rm{de}},k|k - 1}^{(i)}} {g_k}(\mathit{\boldsymbol{z}}|\mathit{\boldsymbol{x}}_{{\rm{de}},k|k - 1}^{(i)})}}} $ (34)
 

则以PHD形式表示的Rényi信息增量评价函数可表示为[20]

$ \begin{array}{*{20}{l}} {\mathcal{R}(\nu ) = \int {{D_{k|k - 1}}} (\mathit{\boldsymbol{x}}){\rm{d}}\mathit{\boldsymbol{x}} + \frac{\alpha }{{1 - \alpha }}\int {{D_{k|k}}} (\mathit{\boldsymbol{x}}){\rm{d}}\mathit{\boldsymbol{x}} - }\\ {{\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} \frac{1}{{1 - \alpha }}\int {{D_{k|k - 1}}} {{(\mathit{\boldsymbol{x}})}^{1 - \alpha }}{D_{k|k}}{{(\mathit{\boldsymbol{x}})}^\alpha }{\rm{d}}\mathit{\boldsymbol{x}}} \end{array} $ (35)
 

进而得到反映传感器分配前后最大决策不确定性目标信息增量的评价函数为[20]

$ \begin{array}{*{20}{c}} {{\mathcal{R}_{{\rm{de}}}}(\nu ) \approx \sum\limits_{i = 1}^{{N_T}} {w_{{\rm{de}},k|k - 1}^{(i)}} + \frac{\alpha }{{1 - \alpha }}\sum\limits_{i = 1}^{{N_{\rm{T}}}} {w_{{\rm{de}},k|k}^{(i)}} - }\\ {\frac{1}{{1 - \alpha }}\sum\limits_{i = 1}^{{N_T}} {{{(w_{{\rm{de}},k|k - 1}^{(i)})}^{1 - \alpha }}} {{(w_{{\rm{de}},k|k}^{(i)})}^\alpha }} \end{array} $ (36)
 

如前文所述, 在本文中选用α=0.5, 因此评价函数可改写为

$ \begin{array}{*{20}{c}} {{\mathcal{R}_{{\rm{de}}}}(\nu ) \approx \sum\limits_{i = 1}^{{N_{\rm{T}}}} {w_{{\rm{de}},k|k - 1}^{(i)}} + \sum\limits_{i = 1}^{{N_{\rm{T}}}} {w_{{\rm{de}},k|k}^{(i)}} - }\\ {2\sum\limits_{i = 1}^{{N_{\rm{T}}}} {\sqrt {w_{{\rm{de}},k|k - 1}^{(i)}w_{{\rm{de}},k|k}^{(i)}} } } \end{array} $ (37)
 

对每一个传感器分配方案νUk, 利用式(37)计算评价函数Rde(ν), 最优的控制准则是使评价函数最大化对应的最优分配方案uk

3.3 算法流程

在SMC-PHD滤波框架下, 本文算法的流程如图 2所示。

图 2 本文所提算法流程 Fig. 2 Algorithm framework proposed in this paper

针对高杂波、电子对抗等实战环境, 本文算法兼顾了目标跟踪与传感器管理过程, 可以较好地抑制高杂波、电子对抗对多目标跟踪及传感器管理决策的影响。首先, 本文采用的自适应SMC-PHD滤波器可以基本排除杂波等无关量测对已有目标滤波的影响, 从而显著提高已有目标的滤波准确性, 为传感器管理提供准确的目标信息;其次, 考虑到电子干扰对目标量测完整性、精度等的影响, 本文综合考虑了目标信息的完整性、信息的质量及信息的内涵等因素, 最大程度的消除了电子干扰对多目标跟踪及传感器管理决策的影响, 保证了算法在电子对抗环境下的可靠性;再次, 本文算法在滤波器的更新过程中, 由于大量的杂波虚警等不参与滤波器的粒子权值计算过程, 因此可以显著提高计算效率, 算法具有更好的实时性。最后, 本文在最优分配传感器选择时, 仅考虑已有目标, 对于潜在的新生目标, 在下一次迭代中结合最新的量测被验证后才开始考虑, 从而在提高计算效率的同时降低了杂波虚警对传感器分配决策的影响。

4 仿真 4.1 仿真场景

为了评估所提出传感器管理算法的性能, 本文采用文献[24]中的仿真场景, 设定监控区域为半径为2 000 m的圆, 防御中心设定为圆心, 共设有3部同类型传感器, 位置分别为(0, 0) m、(-1 000, -500) m、(1 000, -500) m, 参数如表 2所示。

表 2 传感器参数 Table 2 Sensor parameters
传感器参数 数值
特征探测概率 0.5
特征探测距离/km 5
特征RCS/m2 10
特征虚警概率 10-6
恒虚警处理器参考单元 24
天线增益/dB 40
发射功率/W 2 000

T=1 s为传感器采样周期, 仿真时长100 s, 整个仿真过程中共有10个相同类型的目标在不同的时间分别产生于4个不同区域。目标动态过程采用协调转弯模型, 可将单目标状态向量增广为$\boldsymbol{x}_{k}=\left[x_{k}, y_{k}, \dot{x}_{k}, \dot{y}_{k}, \omega_{k}\right]^{\mathrm{T}}$, ω为目标转弯速率, 则目标状态方程可表示为

$ {\mathit{\boldsymbol{x}}_k} = f({\omega _{k - 1}}){\mathit{\boldsymbol{x}}_{k - 1}} + \mathit{\boldsymbol{G}}{\mathit{\boldsymbol{W}}_{k - 1}} $

式中:

$ \begin{array}{l} f\left( \omega \right) = \\ \left[ {\begin{array}{*{20}{c}} 1&0&{{\omega ^{ - 1}}{\rm{sin}}(\omega T)}&{{\omega ^{ - 1}}({\rm{cos}}(\omega T) - 1)}&0\\ 0&1&{{\omega ^{ - 1}}({\rm{cos}}(\omega T) - 1)}&{{\omega ^{ - 1}}{\rm{sin}}(\omega T)}&0\\ 0&0&{{\rm{cos}}(\omega T)}&{ - {\rm{sin}}(\omega T)}&0\\ 0&0&{{\rm{sin}}(\omega T)}&{{\rm{cos}}(\omega T)}&0\\ 0&0&0&0&1 \end{array}} \right], \end{array} $

$\boldsymbol{G}=\left[\begin{array}{ccc} T^{2} / 2 & 0 & 0 \\ T^{2} / 2 & 0 & 0 \\ 0 & T & 0 \\ 0 & T & 0 \\ 0 & 0 & T \end{array}\right], \boldsymbol{W}_{k-1}=\left[W_{x}, W_{y}, W_{\omega}\right]^{\mathrm{T}}$为零均值高斯噪声, 且cov(Wx)=cov(Wy)=15 m/s, cov(Wω)=π/180 rad/s。真实的目标轨迹如图 3所示。

图 3 目标运动真实轨迹 Fig. 3 Trajectories of true targets

支援干扰机位于(0, 2 000) m, 参数如表 3所示。本文所有仿真均采用MATLAB R2012a, 硬件配置为Intel(R) Xeon(R) CPU E5620 @2.4 GHz 2.4 GHz。蒙特卡罗仿真次数为50。在SMC-PHD滤波器的实现过程中, 单目标粒子采样数为1 000, 并规定最小采样数为600。门控方法中, 选择阈值参数τ=100和γ=10-4[16]。目标存活概率Ps=0.99, 杂波是一个泊松RFS, 且在监控区域内服从均匀分布。设置ITSM函数参数为k0=500, m0=1 250[3], k1=0.4, k2=0.6。另外, 本文采用OSPA(Optimal Sub-Pattern Assignment)距离[24]评估多目标跟踪的性能。其定义为:设多目标状态集合为X={x1, x2, …, xn}, 相应的状态估计集合为Y={y1, y2, …, ym}, 若nm, 则OSPA距离定义为

表 3 干扰参数 Table 3 Jamming parameters
干扰参数 数值
干扰功率/W 50
增益/dB 5
损耗/dB 8.5
工作带宽比 5
干扰起始时间/s 10
干扰终止时间/s 80
干扰角度范围/(°) (60, 120)
$ \begin{array}{*{20}{l}} {\bar d_p^{(c)}(\mathit{\boldsymbol{X}},\mathit{\boldsymbol{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} {{\left( {\frac{1}{m}(\mathop {{\rm{min}}}\limits_{\pi \in {\mathit{\Pi} _m}} \sum\limits_{i = 1}^n {{d^{(c)}}} {{({\mathit{\boldsymbol{x}}_i},{\mathit{\boldsymbol{y}}_{\pi (i)}})}^p} + {c^p}(m - n))} \right)}^{1/p}}} \end{array} $ (38)
 

式中:d(c)(x, y)=min(c, ‖x-y‖), Πm为所有{1, 2, …, m}的排列构成的集合。如果n>m, 则dp(c)(Y, X)=dp(c)(X, Y)。如果n=m=0, 则dp(c)(X, Y)=0。距离阶次p≥1, 截断系数c>0。仿真中选取p=1, c=100。

4.2 仿真实验

仿真分别从多目标状态估计精度、计算效率、可靠性3个方面验证本文算法的有效性, 同时仿真中设置了几组不同的对比方法。方法0是基于决策不确定性的传感器管理方法, 即本文方法。方法1是在本文滤波框架下, 仅使用传感器1进行目标跟踪的无传感器管理方法。方法2是在文献[3]方法框架下, 基于多目标整体的信息增量作为评价函数的传感器管理方法。方法3是文献[3]提出的基于最大威胁度目标的传感器管理方法。

设置杂波密度λ=0.000 5 m-2, 不同方案的仿真结果对比如下。图 4为4种方法的多目标状态估计OSPA距离统计对比。在整个仿真过程中, 无论干扰是否存在, 本文方法的平均OSPA距离都是所有方案中最小的。这是由于本文提出的基于决策不确定性的传感器管理方法采用自适应的SMC-PHD滤波器, 降低了杂波对已有目标滤波的影响;同时, 本文方法兼顾传感器量测信息的完整性、质量及目标的威胁度, 可以对目标信息进行择优, 在传感器管理决策时, 排除了杂波虚警的干扰, 从而有效抑制杂波、电子对抗对多目标跟踪及传感器管理决策的影响。方法1采用无传感器管理的方法, 在电子干扰情况下(10~80 s)近乎失效, 其平均OSPA距离最大, 这也说明了电子干扰对传感器管理决策的严重影响。方法2的跟踪效果整体上要略优于方法3, 这与文献[3]中的研究结果也是完全吻合的。同时可以看出, 在监视区域仅存在杂波虚警的情形下(85~100 s), 本文方法可以基本消除杂波对多目标跟踪的影响, 其他方法由于并未对高杂波进行有效的抑制, 对目标数的估计存在较大的误差, 导致其平均OSPA远高于本文方法。

图 4 不同方法OSPA对比 Fig. 4 Comparison of OSPA in different methods

图 5为本文方法与方法3方法对目标数估计的对比。在整个仿真过程中, 方法3方法总体上对目标数过估计, 这是由于方法3方法无法有效抑制杂波虚警, 将虚警误判为真实目标, 不仅影响真实目标的跟踪精度, 而且严重浪费传感器资源。与之相比, 本文方法总体上可以做到对目标数的正确估计。一方面是由于本文提出的传感器管理方法兼顾目标的威胁度和传感器信息的完整性、质量等, 可以对目标信息进行择优, 并优先跟踪决策不确定性最大的目标, 避免了方法3中目标信息考虑不全面而造成的目标状态估计性能退化的现象;另一方面, 本文采用的自适应SMC-PHD滤波器不仅可以避免传统PHD方法对于新生目标先验信息的要求, 而且可以较好的抑制杂波对已有目标的影响, 显著提高已有目标的跟踪精度。

图 5 不同方法目标数评估对比 Fig. 5 Comparison of target number estimation in different methods

单次仿真实验时, 方法0与方法3的传感器选择分别如图 6图 7所示。在电子干扰情形下, 由于传感器1处于干扰机的正前方, 干扰机对传感器1形成主瓣干扰, 因此与方法3相比, 本文方法更多选择干扰较弱的传感器2或传感器3。特别是在监视区域内无真实目标的情况下(85~100 s), 本文方法可以较好地识别杂波虚警, 避免将有限的传感器资源浪费在杂波虚警上, 同时可以在整体上获得较低的平均OSPA距离。

图 6 方法0传感器选择 Fig. 6 Sensor selection result of method 0
图 7 方法3传感器选择 Fig. 7 Sensor selection result of method 3

图 8为本文方法与方法1、方法3在计算效率上的对比。由于方法1并未采用传感器管理, 因此其计算效率最高。在无干扰情形下, 本文方法的计算效率与方法1基本相当, 这说明本文的传感器管理算法并未造成过多的计算负担。特别地, 在整个仿真过程中本文方法的计算效率都明显高于方法3, 其平均计算时长降低了74.3%, 这一方面是由于本文采用的自适应SMC-PHD滤波器仅对已有目标进行更新, 而且基于量测的新生目标采样处理具有较低的计算复杂度, 可以显著提高滤波算法的效率;另一方面, 本文方法在最优传感器的选择时仅考虑已有目标, 并未在第一时间考虑含有大量杂波虚警的新生目标, 从而显著提高传感器管理的效率。

图 8 不同方法计算时长对比 Fig. 8 Comparison of calculation time in different methods

算法的可靠性主要通过整个仿真过程中杂波密度和干扰功率对平均OSPA的影响来验证。本文借鉴文献[3]和文献[26]中对于杂波密度的处理方式, 分别取杂波密度值为0.000 002 2[3]、0.000 01、0.000 1、0.000 5、0.001、0.005 m-2, 并设置2组对比情形。情形1为整个仿真周期0~100 s, 情形2为仅考虑电子对抗存在时的仿真时段10~80 s。本文方法与方法3在不同杂波密度不同情形下的平均OSPA对比如图 9所示。可以看出, 在不同的杂波密度下, 本文方法的平均OSPA均低于方法3, 在杂波密度较低(0.000 002 2)时, 本文方法的平均OSPA略低于方法3, 此时影响OSPA的主要因素是电子干扰, 而本文方法不仅可以抑制杂波, 还可以抑制电子干扰对多目标跟踪及传感器管理决策的影响。在杂波密度为0.000 5时, 本文方法效果最为明显, 平均OSPA相较于方法3显著降低。在高杂波密度(0.005)下, 本文方法的平均OSPA甚至低于杂波密度为0.000 1时的方法3方法, 从而验证本文算法在高杂波环境下的可靠性。在2组不同的情形下, 当杂波密度小于0.000 1时, 方法3在情形1下的平均OSPA低于情形2, 而当杂波密度大于0.000 1时, 正好相反。这是因为对于方法3而言, 当杂波密度较低时, 对平均OSPA影响较为显著的是电子干扰, 当杂波密度较高时, 杂波密度对平均OSPA的影响更为显著。而本文方法在整个仿真过程中, 2种情形下平均OSPA随杂波密度的增长趋势基本相同, 这也证明了本文方法对杂波具有很好的抑制作用, 在高杂波环境下具有较高的可靠性。

图 9 不同情形下平均OSPA随杂波密度变化的对比 Fig. 9 Comparison of average OSPA variation with clutter density in different situations

图 10为本文方法与方法3的平均OSPA随干扰功率的变化对比。可以看出在整个干扰功率变化过程中, 本文方法的平均OSPA明显小于方法3。在干扰功率小于100 W时, 本文方法的平均OSPA基本保持不变, 这是由于本文提出的基于决策不确定性的传感器管理方法可以兼顾目标的威胁度和传感器信息的完整性、质量等, 可以在电子干扰环境下合理的选择量测信息完整的、信息质量高的传感器对目标进行跟踪, 从而最大程度地消除了干扰对目标状态估计的影响。当干扰功率较大时(>100 W), 2种方法的平均OSPA均随着干扰功率的增大而增大, 但是与方法3相比, 本文方法的增加趋势相对缓慢, 平均增长速率约为方法3的74%, 这是因为电子干扰具有方向性, 同样的干扰功率下, 干扰机对不同方向的传感器的干扰效果不同。而本文方法在传感器管理时, 需要在多个不同方位的传感器中综合选择信息更完整的、信息质量更高的传感器对目标进行跟踪, 可以较为有效的抑制电子干扰对目标跟踪的影响, 验证了本文方法在电子干扰环境下的可靠性。

图 10 平均OSPA随干扰功率变化的对比 Fig. 10 Comparison of average OSPA variation with jamming power
5 结论

本文针对高杂波、电子干扰环境, 在自适应SMC-PHD滤波框架下提出了一种基于决策不确定性的传感器管理方法。与基于威胁的传感器管理方法相比, 结论如下:

1) 本文方法可以较好地抑制电子干扰、杂波对多目标估计、传感器管理决策的影响, 相比于对比方法, 其计算效率和跟踪精度明显提升。

2) 在不同杂波密度、干扰功率下, 本文方法性能均优于对比方法, 从而验证本文算法在高杂波、电子对抗环境下的可靠性。

本文仅从传感器对目标的探测和跟踪角度出发, 对传感器决策的不确定性进行描述。在今后的研究中, 可以在传感器管理中考虑更多的影响因素, 如射频隐身、数据链时延等。

参考文献
[1] 张昀普, 单甘霖. 面向空中目标威胁评估的多传感器管理方法[J]. 航空学报, 2019, 40(11): 323218.
ZHANG Y P, SHAN G L. Multi-sensor management approach for aerial target threat assessment[J]. Acta Aeronautica et Astronautica Sinica, 2019, 40(11): 323218. (in Chinese)
Cited By in Cnki | Click to display the text
[2] 高晓光, 李飞, 万开方. 数据丢包环境下的多传感器协同跟踪策略研究[J]. 系统工程与电子技术, 2018, 40(11): 2450-2458.
GAO X G, LI F, WAN K F. Research on multi-sensor cooperative tracking strategy in data packet loss environment[J]. Systems Engineering and Electronics, 2018, 40(11): 2450-2458. (in Chinese)
Cited By in Cnki (3) | Click to display the text
[3] 陈辉, 贺忠良, 连峰, 等. 多目标跟踪中基于目标威胁度评估的传感器控制方法研究[J]. 电子与信息学报, 2018, 40(12): 2861-2867.
CHEN H, HE Z L, LIAN F, et al. Threat assessment based sensor control for multi-target tracking[J]. Journal of Electronics and Information Technology, 2018, 40(12): 2861-2867. (in Chinese)
Cited By in Cnki (5) | Click to display the text
[4] PANG C, SHAN G L, DUAN X S, et al. A multi-mode sensor management approach in the missions of target detecting and tracking[J]. Electronics, 2019, 8(1): 1-18.
Click to display the text
[5] 闫涛, 韩崇昭, 张光华. 空中目标传感器管理方法综述[J]. 航空学报, 2018, 39(10): 022209.
YAN T, HAN C Z, ZHANG G H. An overview of sensor management approaches for aerial target[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(10): 022209. (in Chinese)
Cited By in Cnki (7) | Click to display the text
[6] MARTIN S. Risk-based sensor resource management for field of view constrained sensors[C]//Proceedings of IEEE International Conference on Information Fusion. Piscataway: IEEE Press, 2015: 2041-2048.
[7] MARCOS E G, DOMINIQUE M, PHILIIPPE V, et al. A risk-based sensor management using random finite sets and POMDP[C]//Proceedings of IEEE International Conference on Information Fusion. Piscataway: IEEE Press, 2017: 1588-1596.
[8] KATSILIERIS F, DRIESSEN H, YAROVOY A. Threat-based sensor management for target tracking[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(4): 2772-2785.
Click to display the text
[9] HOU J, JING Z R, YANG Y. Target tracking in standoff jammer using unscented kalman filter and particle filter with negative information[J]. Journal of Shanghai Jiaotong University (Science), 2014, 19(2): 181-189.
Click to display the text
[10] ZHENG Y M, SHI Z G, LU R X, et al. An efficient data-driven particle PHD filter for multitarget tracking[J]. IEEE Transactions on Industrial Informatics, 2013, 9(4): 2318-2326.
Click to display the text
[11] GOSTAR A K, HOSEINNEZHAD R, BAB H A. Multi-bernoulli sensor-selection for multi-target tracking with unknown clutter and detection profiles[J]. Signal Processing, 2016, 119: 28-42.
Click to display the text
[12] TIAN M C, BO Y M, CHEN Z M, et al. A new improved firefly clustering algorithm for SMC-PHD filter[J]. Applied Soft Computing Journal, 2019, 85: 105840.
Click to display the text
[13] MAHLER R P S. Advances in statistical multisource multitarget information fusion[M]. Norwood: Artech House, 2014: 825-860.
[14] WANG X Y, HOSEINNEZHAD R, GOSTAR A K, et al. Multi-sensor control for multi-object bayes filters[J]. Signal Processing, 2018, 142: 260-270.
Click to display the text
[15] HOANG H G, VO B T. Sensor management for multi-target tracking via multi-bernoulli filtering[J]. Automatica, 2014, 50(4): 1135-1142.
Click to display the text
[16] SI W J, WANG L W, QU Z Y. A Measurement-driven adaptive probability hypothesis density filter for multitarget tracking[J]. Chinese Journal of Aeronautics, 2015, 28(6): 1689-1698.
Click to display the text
[17] 董鹏, 敬忠良, 雷明, 等. 基于关联的自适应新生目标强度CPHD滤波[J]. 系统工程与电子技术, 2016, 38(4): 725-731.
DONG P, JING Z L, LEI M, et al. Association based adaptive target birth intensity CPHD filter[J]. Systems Engineering and Electronics, 2016, 38(4): 725-731. (in Chinese)
Cited By in Cnki (2) | Click to display the text
[18] BARTON D K. Radar system analysis and modeling[M]. Norwood: Artech House, 2005: 387-388.
[19] 张睿文, 宋笔锋, 裴扬, 等. 基于ABMS的飞机拦截作战效能评估方法[J]. 系统工程与电子技术, 2018, 40(2): 322-329.
ZHANG R W, SONG B F, PEI Y, et al. Evaluation method for operational effectiveness of aircraft interception based on ABMS[J]. Systems Engineering and Electronics, 2018, 40(2): 322-329. (in Chinese)
Cited By in Cnki (3) | Click to display the text
[20] RISTIC B, VO B N, CLARK D. A note on the reward function for PHD filters with sensor control[J]. IEEE Transactions on Aerospace and Electronic Systems, 2011, 47(2): 1521-1529.
Click to display the text
[21] CLARK D E, BELL J. Multi-target state estimation and track continuity for the particle PHD filter[J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(4): 1441-1453.
Click to display the text
[22] LIU W F, HAN C Z, LIAN F, et al. Multitarget state extraction for the PHD filter using MCMC approach[J]. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(2): 864-883.
Click to display the text
[23] TOBIAS M, LANTERMAN A D. Techniques for birth-particle placement in the probability hypothesis density particle filter applied to passive radar[J]. IET Radar, Sonar and Navigation, 2008, 2(5): 351-365.
Click to display the text
[24] LI T C, JUAN C M, SUN S D, et al. Multi-EAP:Extended EAP for multi-estimate extraction for SMC-PHD filter[J]. Chinese Journal of Aeronautics, 2017, 30(1): 368-379.
Click to display the text
[25] CAI S, RAN X, WANG C. A targets prioritizing method based on clustering coefficient TSM[C]//International Conference on Advances in Materials, Machinery, Electrical Engineering, 2017: 864-869.
[26] YANG X J, XING K Y, FENG X L. Maneuvering target tracking in dense clutter based on particle filtering[J]. Chinese Journal of Aeronautics, 2011, 24(2): 171-180.
Click to display the text
http://dx.doi.org/10.7527/S1000-6893.2020.23781
中国航空学会和北京航空航天大学主办。
0

文章信息

田晨, 裴扬, 侯鹏, 赵倩
TIAN Chen, PEI Yang, HOU Peng, ZHAO Qian
基于决策不确定性的多目标跟踪传感器管理
Decision uncertainty based sensor management for multi-target tracking
航空学报, 2020, 41(10): 323781.
Acta Aeronautica et Astronautica Sinica, 2020, 41(10): 323781.
http://dx.doi.org/10.7527/S1000-6893.2020.23781

文章历史

收稿日期: 2019-12-30
退修日期: 2020-04-07
录用日期: 2020-07-03
网络出版时间: 2020-07-08 17:15

相关文章

工作空间