新生未知自适应噪声协方差TPHD滤波的目标跟踪方法

新生未知自适应噪声协方差tphd滤波的目标跟踪方法
技术领域
1.本发明属于变分贝叶斯滤波技术领域和雷达信号处理技术领域，具体涉及一种新生未知自适应噪声协方差tphd滤波的目标跟踪方法，可用于空中预警、导弹防御和海洋监测等雷达系统中。


背景技术：

2.多目标跟踪技术正在受到越来越多的关注，而现代作战环境中的多目标跟踪受许多不确定因素的影响，因此存在这如下的难点：(1)新生目标出现的位置信息未知，有遗漏新生目标的情况；(2)由于不同时刻跟踪的环境不同，导致量测噪声的协方差矩阵未知或随着时间变化；(3)在对检测到的目标进行跟踪滤波的同时，无法精确完成对目标轨迹的关联，造成轨迹的错误和遗漏。
3.传统的多目标跟踪技术采用的是数据关联技术，但在目标个数较多时会产生计算爆炸的情况。随机有限集理论的提出为跟踪多目标和确认轨迹提供了新的思路，在此之后出现的概率假设密度(phd)滤波器将多目标应用于单目标状态空间进行贝叶斯递推，并且实现了线性高斯环境下高斯混合的闭合解。考虑到传统概率假设密度滤波不能区分出各目标的轨迹，angel等提出了轨迹概率假设密度滤波(tphd)，既可以不评估量测到目标的关联假设，又使用轨迹作为研究对象而不是一组目标以实现不同目标轨迹的区分。tphd滤波在轨迹轨迹精度和滤波时间成本上要明显优于标记phd滤波。


技术实现要素：

4.解决的技术问题：为解决在传感器跟踪场景中目标新生强度未知和量测噪声协方差未知的条件下，传统目标跟踪技术不能准确跟踪导致滤波精度下降的问题，本发明提供一种新生未知自适应噪声协方差tphd滤波的目标跟踪方法，其核心技术在于，以高斯混合轨迹概率假设密度滤波为基础，使用当前时刻获取的量测驱动目标新生强度，解决任意时刻目标新生强度未知的问题；将量测噪声协方差建模为服从逆伽马分布的随机矩阵，采用变分贝叶斯作为框架进行递推滤波，实现同时自适应估计量测噪声协方差和多轨迹联合后验概率密度。
5.技术方案：
6.一种新生未知自适应噪声协方差tphd滤波的目标跟踪方法，所述目标跟踪方法包括以下步骤：
7.s1，初始时刻k＝0时，进行初始化；初始化过程包括以下子步骤：
8.s11，初始化轨迹x＝(t，x
1：i
)的基本状态(wk，mk，pk)，其中x泛指单个目标，t是轨迹启时时刻，i是轨迹的持续时间，mk和pk分别是持续时间内轨迹状态均值和误差协方差矩阵，wk是高斯分量权值；设置轨迹窗口长度为l，单位时间的滤波是针对不多于l个步长下的轨迹进行的预测更新；
9.s12，初始化场景的基本参数：设定量测噪声协方差矩阵r＝diag(β1/α1，
…
βd/αd)，
其中α1，
…
，αd和β1，
…
，βd分别是逆伽马分布的自由度和尺度参数，d是量测噪声协方差矩阵的维数；设定衰减因子ρ；初始化状态转移矩阵f、过程噪声q、量测转移矩阵h、幸存概率ps、检测概率pd，以及场景中杂波率λc和区域大小v；
10.s2，使用k时刻获得的量测集驱动产生目标新生强度，j是量测数目的序号指引，jz是量测总数目，γ作为新生的标记，将k时刻的量测位置信息设定为当前的新生目标均值新生目标权值设定为新生目标误差协方差表示为
11.s3，在变分贝叶斯框架下进行轨迹概率假设密度预测，具体包括以下子步骤：
12.s31，对k-1时刻幸存的轨迹x进行预测，将预测结果和步骤s2的新生目标合并作为所有的预测后的高斯分量；
13.s32，将k-1时刻各幸存轨迹的噪声协方差的参数分别乘以衰减因子ρ；
14.s4，在变分贝叶斯框架下进行轨迹概率假设密度迭代更新，具体包括以下子步骤：
15.s41，根据k时刻量测集zk对预测后的轨迹集进行迭代更新，记录循环次数n，采用变分贝叶斯近似对轨迹状态和噪声协方差矩阵进行解耦合，当轨迹状态趋于定值时停止迭代，确定最终轨迹状态和噪声协方差矩阵；
16.s42，依据检测概率对轨迹集进行漏检更新；
17.s43，将步骤s42和步骤s43的更新结果合并为最终的更新结果；
18.s5，根据轨迹权重在变分贝叶斯框架下进行轨迹概率假设密度剪枝吸收，设定权重阈值，将权值小于权重阈值的分量直接进行删除；设定距离阈值，将权值小于距离阈值的分量直接进行吸收合并；
19.s6，在变分贝叶斯框架下进行轨迹概率假设密度状态提取，提取出权重大于0.5对应的高斯分量和噪声协方差矩阵，如果在滤波时间内则转入步骤s2，将预设的初始噪声协方差直接作为k+1时刻初始量测噪声协方差的设定值，超出滤波时间则目标跟踪结束。
20.进一步地，步骤s2中，使用k时刻获得的量测集驱动产生的目标新生强度d
γ，k|k-1
(x，r)为：
[0021][0022][0023][0024][0025]
式中n(x)表示轨迹x服从正态分布，ig(r)表示量测协方差r服从逆伽马分布，表示逆伽马参数的维度d的序号，α和β表示逆伽马分布自由度和尺度参数，h表示量测转移矩阵。
[0026]
进一步地，步骤s3中，采用下述公式对各轨迹高斯分量与逆伽马的混合乘积d
k|k-1
(x，r)进行预测：
[0027]dk|k-1
(x，r)＝d
s，k|k-1
(x，r)+d
γ，k|k-1
(x，r)
[0028][0029][0030][0031][0032][0033][0034]
式中，为轨迹j的状态转移矩阵，f为单目标的状态转移矩阵，q是过程噪声协方差，ρ是衰减因子，和是k-1时刻的均值和协方差，和是预测的均值和协方差，和分别为k-1时刻的逆伽马分布的自由度和尺度参数，和分别为逆伽马分布的自由度和尺度参数。
[0035]
进一步地，步骤s41中，采用下述公式，根据量测值计算新息调整均值和误差协方差并不断迭代更新，得到最优的轨迹状态和量测噪声协方差：
[0036][0037][0038][0039][0040]
其中，h是量测转移密度，是轨迹量测转移密度，是预测的位置信息，是更新后的自由度；
[0041]
为求取最优参数，调整更新部件，记录循环次数n，开始循环：
[0042][0043][0044][0045][0046][0047][0048]
其中，表示更新后的量测协方差，表示n次迭代后逆伽马分布的尺度参数，是新息的协方差，是n次迭代后的轨迹均值，是n次迭代后的轨迹协方差。
[0049]
当基本处于定值时，停止循环，输出以及以此计算更新后的各高
斯分量权值
[0050][0051]
其中，表示新息的协方差，表示预测权值，pd表示检测概率，λcv是杂波强度。
[0052]
进一步地，步骤s5中，设定距离阈值u，如果分量满足以下条件：则直接进行吸收合并，其中表示不同的轨迹均值，表示轨迹协方差。
[0053]
本发明以轨迹概率假设密度滤波为基础，区分各个目标轨迹的同时，重点解决目标新生强度未知和量测噪声的协方差矩阵未知的问题。
[0054]
有益效果：
[0055]
(1)本发明的新生未知自适应噪声协方差tphd滤波的目标跟踪方法，利用轨迹概率假设密度滤波(tphd)为基础，跟踪各目标的同时，区分出不同的轨迹。
[0056]
(2)本发明的新生未知自适应噪声协方差tphd滤波的目标跟踪方法，利用当前时刻获取的量测驱动目标新生强度，解决任意时刻目标新生强度未知的问题。
[0057]
(3)本发明的新生未知自适应噪声协方差tphd滤波的目标跟踪方法，利用采用变分贝叶斯作为框架，将量测噪声协方差建模为服从逆伽马分布的随机矩阵，提高轨迹概率假设密度滤波算法的鲁棒性和跟踪精度，实现同时自适应估计量测噪声协方差和多轨迹联合后验概率密度。
附图说明
[0058]
图1为本发明实施例的新生未知自适应噪声协方差tphd滤波的目标跟踪方法流程图。
[0059]
图2是本发明实施例中仿真1的真实轨迹图。
[0060]
图3是本发明实施例中仿真1的轨迹估计图。
[0061]
图4是本发明实施例中仿真1经100次蒙特卡洛实验的平均噪声标准差估计图。
[0062]
图5是本发明实施例中仿真1经100次蒙特卡洛实验的平均gospa误差对比图。
[0063]
图6是本发明实施例中仿真1经100次蒙特卡洛实验的平均漏检目标代价对比图。
[0064]
图7是本发明实施例中仿真2的真实轨迹图。
[0065]
图8是本发明实施例中仿真2的轨迹估计图。
[0066]
图9是本发明实施例中仿真2经100次蒙特卡洛实验的平均噪声标准差估计图。
[0067]
图10是本发明实施例中仿真3经100次蒙特卡洛实验的平均噪声标准差估计图。
[0068]
图11是本发明实施例中仿真4经100次蒙特卡洛实验的平均噪声标准差估计图。
具体实施方式
[0069]
下面的实施例可使本专业技术人员更全面地理解本发明，但不以任何方式限制本发明。
[0070]
本实施例采用量测驱动实现当前时间内目标新生强度的估计和采用变分贝叶斯作为框架进行递推迭代出量测噪声协方差的基于轨迹概率假设密度滤波的多目标跟踪方法。
[0071]
新生未知自适应噪声协方差tphd滤波的目标跟踪方法为：以高斯混合tphd滤波为基础，使用当前时刻获取的量测驱动目标新生强度，解决任意时刻目标新生强度未知的问题；将量测噪声协方差建模为服从逆伽马分布的随机矩阵，采用变分贝叶斯作为框架进行递推滤波，根据量测值计算并调整新息不断迭代更新，实现同时自适应估计量测噪声协方差和多轨迹联合后验概率密度。
[0072]
参见图1，本实施例公开了一种新生未知自适应噪声协方差tphd滤波的目标跟踪方法，所述目标跟踪方法包括以下步骤：
[0073]
步骤1、初始时刻k＝0时，初始化具体内容：
[0074]
(1a)初始化轨迹x＝(t，x
1：i
)的基本状态(wk，mk，pk)，其中t是轨迹启时时刻，i是轨迹的持续时间，mk和pk分别是持续时间内轨迹状态均值和误差协方差矩阵，高斯分量权值为wk，设置轨迹窗口长度为l即单位时间的滤波是针对不多于l个步长下的轨迹进行的预测更新；
[0075]
(1b)初始化场景的基本参数。设定量测噪声协方差矩阵r＝diag(β1/α1，
…
βd/αd)，其中α和β分别是逆伽马分布的自由度和尺度参数，d是量测噪声协方差矩阵的维数，ρ是设定的衰减因子。初始化状态转移矩阵f，过程噪声q，量测转移矩阵h，幸存概率ps，检测概率pd，场景中杂波率和区域大小分别为λc和v。
[0076]
步骤2、使用k时刻获得的量测集驱动产生目标新生强度d
γ，k|k-1
(x，r)，将k时刻的量测位置信息设定为当前的新生目标均值新生目标权值设定为新生目标误差协方差可以表示为
[0077]
目标新生强度的设定：
[0078][0079][0080][0081][0082]
式中n(x)表示轨迹x服从正态分布，ig(r)表示量测协方差r服从逆伽马分布，表示逆伽马参数的维度d的序号，α和β表示逆伽马分布自由度和尺度参数，h表示量测转移矩阵。
[0083]
步骤3、在变分贝叶斯框架下进行轨迹概率假设密度预测，具体步骤：
[0084]
(3a)对k-1时刻幸存的轨迹x进行预测，和步骤2的新生目标合并作为所有的预测后的高斯分量；
[0085]
(3b)对k-1时刻各幸存轨迹的噪声协方差的参数分别乘以衰减因子ρ。值得注意的是，逆伽马分布的自由度和尺度参数原基础上增加了一个因子ρ，从而通过vb方法通过迭代提高了稳定性。
[0086]
各轨迹高斯分量与逆伽马的混合乘积进行的预测，具体为：
[0087]dk|k-1
(x，r)＝d
s，k|k-1
(x，r)+d
γ，k|k-1
(x，r)
[0088][0089][0090][0091][0092][0093][0094]
式中，为轨迹j的状态转移矩阵，f为单目标的状态转移矩阵，q是过程噪声协方差，ρ是衰减因子，和是k-1时刻的均值和协方差，和是预测的均值和协方差，和分别为k-1时刻的逆伽马分布的自由度和尺度参数，和分别为预测的逆伽马分布的自由度和尺度参数。
[0095]
步骤4、在变分贝叶斯框架下进行轨迹概率假设密度迭代更新，具体步骤：
[0096]
(4a)根据k时刻量测集zk对预测后的轨迹集进行迭代更新，记录循环次数n，采用变分贝叶斯近似对轨迹状态和噪声协方差矩阵进行解耦合，当轨迹状态趋于定值时停止迭代，确定最终轨迹状态和噪声协方差矩阵；
[0097]
根据量测值计算新息调整均值和误差协方差并不断迭代更新，得到最优的轨迹状态和量测噪声协方差，具体为：
[0098][0099][0100][0101][0102]
其中，h是量测转移密度，是轨迹量测转移密度，是预测的位置信息，是更新后的自由度。
[0103]
为求取最优参数，调整更新部件，记录循环次数n，开始循环：
[0104][0105][0106][0107]
[0108][0109][0110]
其中，表示更新后的量测协方差，表示n次迭代后逆伽马分布的尺度参数，是新息的协方差，是n次迭代后的轨迹均值，是n次迭代后的轨迹协方差。
[0111]
当基本处于定值时，停止循环，输出以及以此计算更新后的各高斯分量权值
[0112][0113]
其中，表示新息的协方差，表示预测权值，pd表示检测概率，λcv是杂波强度。
[0114]
(4b)依据检测概率对轨迹集进行漏检更新，(4a)和(4b)结果合并为更新结果。
[0115]
步骤5、根据轨迹权重在变分贝叶斯框架下进行轨迹概率假设密度剪枝吸收，设定权重阈值，对权值小于阈值的分量直接进行删除；设定距离阈值u，将权值小于阈值的分量
[0116][0117]
其中表示不同的轨迹均值，表示轨迹协方差。
[0118]
直接进行吸收合并，这样减少了滤波中的高斯分量，提高运算效率。
[0119]
步骤6、在变分贝叶斯框架下进行轨迹概率假设密度状态提取，提取出权重大于0.5对应的高斯分量和噪声协方差矩阵，如果在滤波时间内则继续进行步骤2，将预设的初始噪声协方差直接作为k+1时刻初始量测噪声协方差的设定值，超出滤波时间则目标跟踪结束。
[0120]
本发明的算法和处理方法已通过验证，取得了满意的应用效果：
[0121]
1.实验条件：假设存在多个目标在水平直角坐标系x-y平面上作匀速直线运动，目标的运动状态可以表示为{x，v
x
，y，vy}，其中{x，y}分别是目标在水平直角坐标系上x方向和y方向的位置信息，{v
x
，vy}分别是目标在水平直角坐标系上x方向和y方向的速度信息。初始化状态转移矩阵f，过程噪声q，量测转移矩阵h，幸存概率ps，检测概率pd，场景中杂波率和区域大小分别为λc和v，单位时间t取1。
[0122]
[0123][0124][0125]
ps＝0.99。
[0126]
pd＝0.98。
[0127]
λc＝5 v＝10002＝106。
[0128]
评价参数gospa＝10。
[0129]
衰减因子ρ为0.95.
[0130]
场景1中设定噪声协方差矩阵r参数未知且固定，r＝diag(4，4)。
[0131]
场景2中设定噪声协方差矩阵r参数未知且时变，当时刻t小于90时r＝diag(4，4)，当时刻t超过90时r＝diag(16，16)。
[0132]
场景3中设定噪声协方差矩阵r参数未知且时变，当时刻t小于50时r＝diag(4，4)，当时刻t超过50时且小于100时r＝diag(16，16)，当时刻t超过100时且小于140时r＝diag(25，25)，当时刻t超过140时r＝diag(16，16)。
[0133]
场景4中设定噪声协方差矩阵r参数未知且时变，伴随时间步长成sin函数变化。
[0134]
四个场景r参数初始值设定为α1＝α2＝1 β1＝β2＝1。
[0135]
2.仿真内容：
[0136]
仿真1：设计4个目标均作匀速直线运动，各目标出生位置速度信息分别为{100，50/7，400，20/7}、{450，-20/7，300，25/7}、{100，80/13，150，10/13}、{100，7.5/2，400，-15/4}，出生死亡时间分别为{2，75}、{2，75}、{15，85}、{20，100}，跟踪时长设定为100，图2是本场景的真实轨迹图，圆圈表示起始三角表示消亡，图3是本场景的轨迹估计图，图4是本场景的噪声标准差估计图。图5是本场景的gospa误差对比图，将错误偏大噪声协方差估计和错误偏小噪声协方差估计得到的gospa误差与本发明gospa误差及真实噪声协方差估计进行对比，图6是本场景的上述四种估计漏检目标代价对比图。
[0137]
仿真2：设计4个目标均作匀速直线运动，各目标出生位置速度信息分别为{100，50/14，400，20/14}、{450，-20/14，300，25/14}、{100，40/13，150，5/13}、{100，7.5/4，400，-15/8}，出生死亡时间分别为{2，149}、{2，149}、{15，156}、{20，181}，跟踪时长设定为180，作100次蒙特卡洛实验，图7是本场景的真实轨迹图，圆圈表示起始三角表示消亡，图8是本场景的轨迹估计图，图9是本场景的噪声标准差估计图。
[0138]
仿真3：设计4个目标均作匀速直线运动，跟踪参数与仿真2相同，作100次蒙特卡洛实验，图10是本场景的噪声标准差估计图。
[0139]
仿真4：设计4个目标均作匀速直线运动，各目标出生位置速度信息分别为{100，
50/56，400，20/56}、{450，-20/56，300，25/56}、{100，40/52，150，5/52}、{100，7.5/32，400，-15/24}，出生死亡时间分别为{2，593}、{2，593}、{15，582}、{20，667}，跟踪时长设定为666，作100次蒙特卡洛实验，图11是本场景的噪声标准差估计图。
[0140]
3.仿真结果分析：
[0141]
从图3和图8可以看出，本专利所提方法在没有设置新生目标强度的条件下仍然可以很好地实现对多目标的跟踪。
[0142]
从图4、图9、图10、图11可以看出，本专利所提方法既可以对未知不变的噪声协方差进行估计，也可以实时地估计出未知时变的噪声协方差并取得不错的效果。
[0143]
从图5和图6可以看出，本专利所提方法gospa误差要小于噪声协方差错误估计的tphd滤波，基本近似于真实估计噪声协方差的轨迹概率假设密度滤波，滤波精度高于错误估计的滤波算法；本专利所提方法漏检目标也要小于噪声协方差错误估计的轨迹概率假设密度滤波，基本近似于真实估计噪声协方差的轨迹概率假设密度滤波。
[0144]
仿真结果表明了本专利所提新生未知自适应噪声协方差tphd滤波的目标跟踪方法的有效性。
[0145]
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。