基于PHD滤波的临近空间高超声速目标跟踪方法

本发明涉及信号处理技术领域，更具体的，涉及一种基于phd滤波的临近空间高超声速目标跟踪方法。

背景技术：

近年来，随着航天技术的发展，基于助推-滑翔弹道概念的高超声速飞行器逐渐成为国内外研究的热点。这类高超声速飞行器可在大气中滑翔飞行，进行大范围快速机动，具有射程远，速度快，机动性强等优点。对于单目标高超声速飞行器跟踪，以扩展卡尔曼滤波(extendedkalmanfilter,ekf)为代表的非线性滤波算法发挥了重要的作用对高机动目标取得了良好的跟踪效果，实现了对弹道式再入目标的高精度跟踪，如文献[jawahara,koteswararaos.modifiedpolarextendedkalmanfilter(mp-ekf)forbearings-onlytargettracking[j].indianjournalofenceandtechnology,2016,9(26):1-5]、[belikbv,belovsg.usingofextendedkalmanfilterformobiletargettrackinginthepassiveairbasedradarsystem☆[j].procediacomputerscience,2017,103:280-286]。

在强杂波环境中，杂波量测与目标量测难以区分，初始目标数目未知，同时，针对高超声速飞行器携带多个滑翔弹头，利用诱饵弹增强其突防能力的问题，有必要讨论高超声速多目标跟踪。传统的多目标跟踪以数据关联为核心，代表算法有多级假设检验方法(multiplehypothesistracking,mht)和联合概率数据关联(jointprobabilisticdataassociation,jpda)，二者在目标数目和杂波数增加时，算法复杂度呈指数增长，即出现“组合爆炸”。基于随机有限集(randomfinitesets,rfs)理论的多目标跟踪算法避免了数据关联，其中最重要的工作是mahler等人提出的概率假设密度(probabilityhypothesisdensity,phd)滤波器[mahlerr.p.s..multitargetbayesfilteringviafirst-ordermultitargetmoments[j].ieeetransactionsonaerospaceandelectronicsystems,2003,39(4):1152-1178]。其主要实现方法有高斯混合phd(gaussianmixture-phd,gm-phd)滤波器和序贯蒙特卡洛phd滤波器(sequentialmontecarlo-phd,smc-phd)。

目前，使用phd滤波器进行高超声速多目标跟踪存在以下两个问题：

(1)传统的phd滤波器通常假定已知初始目标的起始位置信息和新生目标强度函数，在实际应用中，不能直接得到上述信息。

(2)目标量测附近的杂波干扰往往会导致滤波器性能下降，出现错跟，目标数目估计过大的情况。

技术实现要素：

本发明为了解决传统gmphd算法需要知道初始目标的起始位置信息和新生目标强度函数，且杂波干扰导致滤波器性能下降的问题，提供了一种基于phd滤波的临近空间高超声速目标跟踪方法，其实现了强杂波下对未知数目的临近空间飞行器进行精确跟踪。

为实现上述本发明目的，采用的技术方案如下：一种基于phd滤波的临近空间高超声速目标跟踪方法，所述的方法包括步骤如下：

s1：使用初始时刻的量测进行对滤波器初始化，将前四个时刻的雷达量测值从雷达站球坐标系转到雷达站enu直角坐标系，采用两点差分法，得到初始目标强度函数；

s2：将初始目标输入初始化后的滤波器中，预测步骤按照gm-phd滤波器的过程进行计算，得到预测目标集；

s3：根据当前时刻量测与由上一时刻预测得到的目标位置之间的马氏距离，将量测划分为存活目标量测与杂波量测；

s4：利用当前时刻的量测值得到滤波器当前的更新值，完成对存活目标的更新；

s5：设置修剪门限、合并门限对步骤s4的更新公式中的高斯项进行剪枝和合并；

s6：计算修剪合并后的强度函数、提取目标状态，完成对目标状态估计。

优选地，步骤s1，初始目标强度函数的表达式如下：

其中，np(x-m)表示自变量为x，均值为m，协方差为p的高斯概率密度函数，为第j个高斯项的初始权值，为第j个高斯项的初始协方差。

表示初始目标状态，其表达式如下：

式中，x1j(k)，x2j(k)，x3j(k)为k时刻第j个量测在雷达站enu直角坐标系下x、y、z方向的坐标；为量测方程在x、y、z维度的反函数，为k时刻的第j个量测；

设所有初始高斯项权值为1。

进一步地，将初始目标输入初始化后的滤波器进行预测，其中初始目标的状态转移的方程表示如下：

x(k+1)＝f(k)x(k)+w(k)

其中，f(k)为状态转移矩阵，有

其中，fij为零矩阵，且i≠j，

式中：

p1＝(2-2αt+α²t²-2e^-αt)/(2α³)

q1＝(αt-1+e^-αt)/α²

r1＝(1-e^-αt)/α

s1＝e^-αt

w(k)表示均值为0，协方差为q(k)的高斯白噪声序列；

预测步骤具体计算公式如下：

vk|k-1(xk|z1:k-1)＝vs,k|k-1(xk|z1:k-1)+yk(xk)

其中

其中，yk(xk)为新生目标强度函数。

再进一步地，步骤s3，划分量测具体如下：

从时刻2开始，对k时刻的量测进行划分，记：

其中，sk为k时刻高斯项个数，zk为所有量测集合；

其中，rk为量测噪声协方差；

则落在门限值内的量测集合作为存活目标量测，

落在门限值外的量测集合作为杂波量测，

其中t表示门限值。

再进一步地，所述的门限值由下式确定，若pg为正确量测落入确认区域内的概率，则有

t＝-2ln(1-pg)。

再进一步地，步骤s4，使用量测集合对存活目标进行更新，更新式如下：

其中，

其中，hk为量测方程h(·)在k时刻的雅可比矩阵；kk(zk)表示杂波强度，pd表示目标检测概率、rk为量测噪声协方差矩阵；为高斯函数；hk^t为hk的转置再进一步地，需要重新计算杂波强度kk(zk)，

设v为整个观测区域面积，λ为杂波平均数，vk为新的观测区域，不考虑重叠的情况，将vk视为所有量测对应的门限区域面积之和，即

则杂波强度

再进一步地，步骤s5，所述的剪枝和合并，具体如下：

令l＝0，t_prun为修剪门限；

当集合i非空，重复以下过程：

l＝l+1，

i＝i\l，直到

其中，u_merg表示合并门限。

再进一步地，步骤s5，修剪合并后的强度函数为：

再进一步地，步骤s5，根据下式提取目标状态：

本发明的有益效果如下：

本发明能够实现对多个高机动目标的跟踪，尤其在目标飞行轨迹高度非线性的情况表现良好，同时还可以克服传统gmphd算法需要初始目标信息的缺陷，并且能够在强杂波环境下对未知数目临近空间飞行器进行跟踪。

附图说明

图1是本实施例所述的临近空间高超声速目标跟踪方法的步骤流程图。

图2是本实施例所述的助推-滑翔弹道的示意图。

图3是通过本实施例所述的方法的跟踪结果的示意图。

图4是通过本实施例所述的方法得到的前100秒的目标数目估计的效果图，图中，a表示传统ek-gmphd方法，b表示本实施例所述的目标跟踪方法。

图5是采用wasserstein距离评价目标状态估计的效果图。

具体实施方式

下面结合附图和具体实施方式对本发明做详细描述。

实施例1

如图1所示，一种基于phd滤波的临近空间高超声速目标跟踪方法，所述的方法包括步骤如下：

步骤s1：使用初始时刻的量测进行对滤波器初始化，将前四个时刻的雷达量测值从雷达站球坐标系转到雷达站enu直角坐标系，采用两点差分法，得到初始目标强度函数；

所述的初始目标强度函数的表达式如下：

其中，np(x-m)表示自变量为x，均值为m，协方差为p的高斯概率密度函数；为第j个高斯项的初始权值；为第j个高斯项的初始协方差。

表示初始目标状态，其表达式如下：

式中，x1j(k)，x2j(k)，x3j(k)为k时刻第j个量测在雷达站enu直角坐标系下x、y、z方向的坐标；为量测方程在x、y、z维度的反函数，为k时刻的第j个量测；

设所有初始高斯项权值为1。

步骤s2：将初始目标输入初始化后的滤波器中，预测步骤按照gm-phd滤波器的过程进行计算，得到预测目标集；

其中初始目标的状态转移的方程表示如下：

x(k+1)＝f(k)x(k)+w(k)

其中，f(k)为状态转移矩阵，有

其中，fij为零矩阵，且i≠j，

式中：

p1＝(2-2αt+α²t²-2e^-αt)/(2α³)

q1＝(αt-1+e^-αt)/α²

r1＝(1-e^-αt)/α

s1＝e^-αt

w(k)表示均值为0，协方差为q(k)的高斯白噪声序列。

预测步骤具体计算公式如下：

vk|k-1(xk|z1:k-1)＝vs,k|k-1(xk|z1:k-1)+yk(xk)

其中

其中，yk(xk)为新生目标强度函数。

步骤s3：根据当前时刻量测与由上一时刻预测得到的目标位置之间的马氏距离，将量测划分为存活目标量测与杂波量测；

从时刻2开始，对k时刻的量测进行划分，记：

其中，sk为k时刻高斯项个数，zk为所有量测集合；

其中，rk为量测噪声协方差；

则落在门限值内的量测集合作为存活目标量测，

落在门限值外的量测集合作为杂波量测，

其中t表示门限值。

所述的门限值由下式确定，若pg为正确量测落入确认区域内的概率，则有

t＝-2ln(1-pg)。

s4：利用当前时刻的量测值得到滤波器当前的更新值，完成对存活目标的更新；

使用量测集合对存活目标进行更新，更新式如下：

其中，

其中，hk为量测方程h(·)在k时刻的雅可比矩阵；kk(zk)表示杂波强度，pd表示目标检测概率、rk为量测噪声协方差矩阵，hk^t为hk的转置。

上述方程组中杂波强度kk(zk)需要重新计算，具体如下：

设v为整个观测区域面积，λ为杂波平均数，vk为新的观测区域，不考虑重叠的情况，将vk视为所有量测对应的门限区域面积之和，即

则杂波强度

步骤s5：设置修剪门限、合并门限对步骤s4的更新公式中的高斯项进行剪枝和合并；

令l＝0，t_prun为修剪门限；

当集合i非空，重复以下过程：

l＝l+1，

i＝i\l，直到

其中，u_merg表示合并门限。

s6：计算修剪合并后的强度函数、提取目标状态，完成对目标状态估计。

修剪合并后的强度函数为：

根据下式提取目标状态：

为了进一步对本实施例所述的临近空间高超声速目标跟踪方法的技术效果，具体举例说明如下：

在雷达站enu三维直角坐标系下跟踪两个高机动目标，扫描周期t＝1s，目标初始状态如表1：

表1目标初始状态

临近空间超高声速飞行器的空气动力模型是以下一组非线性常微分方程组，

其中，u为地球重力常数，u＝3.986005×10¹⁴m³/s²，re为地球半径，ad为阻力参数，at为转弯力参数，ac为爬升力参数。ρ为空气密度，为目标高度h的指数函数，

ρ(h)＝ρ0e^-κh

其中ρ0＝1.2250kg/m³为海平面大气密度，κ＝1.186×10^-4。

使用四阶龙哥库塔法求解微分方程组的数值解，得到仿真的弹道，微分方程组初值如表1，得到助推-滑翔弹道如图2所示。

设置机动频率α＝0.5，目标在x，y，z方向上的加加速度方差分别为7，2，9，量测噪声标准差分别为目标存活概率ps＝0.99，目标检测概率pd＝0.98，量测落入确认区域概率pg＝0.99，高斯项修剪门限t_prun＝10^-5，合并门限u_merg＝4，杂波率为2。具体的仿真实验的步骤如下：

(1)按照上述的仿真参数设置对变量进行初始化，以及按照本实施例所述的方法的步骤s1对滤波器初始化；

(2)按照本实施例所述的方法的步骤s2得到预测目标集；

(3)按照本实施例所述的方法的步骤s3对得到的量测进行划分；

(4)按照本实施例所述的方法的步骤s4利用当前时刻的量测值得到滤波器当前的更新值，其中量测方程h(·)为

其中，zk＝[rae]^t为k时刻的量测值，r为目标到雷达站的距离，a为方位角，e为俯仰角。量测噪声vk^r，vk^a，vk^e相互独立，服从零均值高斯分布。量测方程h(·)的雅可比矩阵hk为

其中，

(5)按照本实施例所述的方法的步骤s5对高斯项进行剪枝和合并；

(6)按照本实施例所述的方法的步骤s6提取目标状态；

(7)循环步骤s2-s6直到结束。

最终的跟踪结果见附图3，最终得到的前100秒的目标数目估计如附图4所示，其中a代表传统ek-gmphd方法，b代表本实施例所述的目标跟踪方法。

可以看出，传统ek-gmphd方法会造成目标数目过估，采用目标估计绝对误差的时间平均作为目标数目估计的评价指标，即：

其中，为第k秒目标数目估计值，xk为第k秒真实目标数目，传统ek-gmphd方法的目标数目估计误差为0.1875，本实施例所述的目标跟踪方法的目标数目估计误差为0.0643，仅为传统方法的34.27％。可见使用本实施例所述的目标跟踪方法效果明显优于传统方法，目标数目估计更准确。

采用wasserstein距离评价目标状态估计精度，wasserstein距离是度量两组点集相似程度的指标，其数值越小，目标估计精度越高，实验结果见附图5，其中a代表传统ek-gmphd方法，b代表本实施例所述的目标跟踪方法。

可以看出，使用本实施例所述的目标跟踪方法目标状态估计精度要优于传统方法。传统ek-gmphd方法的wasserstein距离的时间平均值为37704，本实施例所述的目标跟踪方法的wasserstein距离的时间平均值为19340，仅为传统ek-gmphd方法的51.29％。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。