CN115358325A - 未知概率Skew及重尾噪声下的目标跟踪方法 - Google Patents

Abstract

本发明公开了一种未知概率Skew及重尾噪声下的目标跟踪方法，包含以下步骤：首先，在分布式目标跟踪框架下，根据未知参数所服从的特性进行先验建模；其次通过标准变分贝叶斯方法定点迭代联合推断出目标状态、系统偏差、噪声协方差等后验分布参数；最后依据协方差交叉融合策略实现对局部状态估计值的融合与修正。所提方法在目标跟踪过程中综合考虑了未知概率随机出现的过程噪声、异常重尾量测噪声和未知且时变的系统偏差的多重影响，能够有效地估计出目标状态、噪声协方差、系统偏差等未知参数，进而提高了目标的跟踪精度，同时具有较好的自适应性和鲁棒性。

Description

未知概率Skew及重尾噪声下的目标跟踪方法

技术领域

本发明涉及多传感器目标跟踪技术领域，尤其涉及一种过程噪声为未知概率随机出现Skew噪声或重尾噪声、未知且时变系统偏差及重尾量测噪声的目标跟踪方法。

背景技术

目前，当今战场环境复杂多变，电磁环境愈发复杂，采用单一传感器信息对导弹进行制导必然面临探测不精确、易受敌欺骗干扰、目标信息易中断等问题。因此，在传感器网络中，从各种传感器收集的信息被协同融合，以提高系统整体性能。然而融合过程依赖于传感器的精确配准，此外实际工程中目标往往受外界环境的影响，严重的过程噪声可能使目标产生离群值，故在目标跟踪过程中过程噪声的影响就无法被忽略，此外传感器在对目标观测时同样会遭受量测噪声的影响，使用错误/不准确的过程噪声和量测噪声可能会导致严重的估计误差，甚至使滤波发散，故进行传感器偏差配准的同时对噪声信息进行估计，进而提高目标的跟踪精度。

一些传统的配准算法如最小二乘(LS)法和最大似然估计(MLE)法，在不考虑噪声或噪声很小的情况下能够有效地估计出系统偏差，然而这些算法属于离线类算法，缺乏实时性。近些年来为了减轻传感器的通信负担，提高传感器的工作效率，分布式传感器网络应运而生，Okello等人提出了一种用于分布式航迹级配准的等效测量方法，即将传感器偏差扩维至状态向量。该方法设计了两个并行扩展卡尔曼滤波器(EKF)，用于估计量测航迹关联后同一目标的状态和传感器偏差。虽然该算法易于实现并且在线处理，实时性得到了提高，并且考虑了噪声影响，但噪声仅为高斯白噪声且时不变的情况下，收敛不合理且状态扩维增加了通信负担。Y.Huang等人在目标跟踪中应用了变分贝叶斯(VB)方法，VB是一种在线处理算法，综合考虑了过程噪声和量测噪声的影响，首先对噪声信息协方差矩阵进行先验建模，其次采用变分迭代求解,由于考虑了噪声影响提高了目标跟踪精度，但未考虑系统偏差的影响。Y.Huang等人另外运用变分迭代的思想，除了考虑量测噪声的影响外，还将系统偏差以高斯分布建模，然后利用VB定点迭代估计出系统偏差和量测噪声协方差矩阵等后验参数，但未考虑目标本身可能遭受外界过程噪声的影响。此外，目前状态过程噪声为未知概率随机出现重尾噪声或Skew噪声的研究还未有相关学者提出。

综上，现有的算法较少能兼顾处理目标跟踪过程中出现的各种问题，具体如下：

1.在进行目标跟踪过程中，目标本身可能随机遭受未知概率随机出现的重尾过程噪声或Skew过程噪声的影响，进而使目标自身产生离群值，对目标跟踪造成困难；

2.传感器量测值可能遭受异常重尾量测噪声和未知且时变的系统偏差的双重影响，进而使传感器量测信息不准确，甚至丢失目标；

3.实际工程中工作环境复杂多变，电磁环境愈发复杂，采用单一传感器信息、必然面临探测不精确、易受敌欺骗干扰、目标信息易中断等问题。

发明内容

本发明的目的是提供一种未知概率Skew及重尾噪声下的目标跟踪方法，能够处理多传感器跟踪目标过程中涉及过程噪声为未知概率随机出现的Skew噪声或重尾噪声、未知且时变系统偏差以及重尾量测噪声，进而提高了目标跟踪精确度和鲁棒性。

本发明采用的技术方案为：

未知概率Skew及重尾噪声下的目标跟踪方法，包括

包括以下步骤：

S1：多传感器共同跟踪同一目标的物理场景，在目标跟踪过程中综合考虑目标状态过程噪声、传感器量测噪声和系统偏差的影响，故建立目标的状态模型、量测模型和系统偏差模型；

S2：设置k-1时刻变分迭代求解参数的初值和初始化期望，k＝1,2,3...,ts，ts为仿真总时间；

S3：建立k时刻传感器s的目标状态一步预测概率密度P(x_k,s∣z_1:k-1)模型、量测似然概率密度P(z_k,s∣x_k,s)模型和对数联合概率密度

模型；

S4：变分贝叶斯近似后验：

步骤S3构建的模型参数之间存在相互耦合，无法解析求解未知参数的后验，引入变分

贝叶斯方法近似求解S3中的未知参数的后验，其求解公式如下：

其中，

表示关于

求期望操作，

表示取

集的某一个元素，

是

集中除

外剩余元素，

相对未知参数

是不相关的常数，

表示以自然常数

为底的对数近似后验概率密度；

S5：设置变分迭代总次数为

i表示变分第i次迭代，

S6：更新求解状态预测协方差矩阵P_k|k-1,s和量测噪声协方差矩阵R_k,s；

S7：更新求解形状参数β_k、系统偏差η_k,s和目标状态x_k,s；

S8：更新求解辅助变量λ_k、γ_k及其对应自由度参数ω_k、ν_k；

S9：更新求解随机变量ξ_k和Skew噪声发生概率π_k；

S10：重复步骤S6-S9，直到达到预先设定变分迭代总次数或满足迭代终止条件:

ε一般取1e-15，则迭代结束；

S11：输出变分迭代结果:

输出迭代结果：

和Σ_k|k,s，并赋值作为k+1时刻变分迭代初始值；

S12：多传感器分布式融合反馈:

在分布式网络中，在其局部滤波器完成偏差配准后仅依赖单个传感器进行目标跟踪往往很难达到预期的精度，故需将步骤S11变分迭代求解得到的各局部滤波器状态值

以及其协方差矩阵P_k|k,s信息进行融合并将融合后的最优目标状态

及其协方差矩阵

反馈给各局部滤波器。

所述步骤S1中离散目标状态模型、量测模型和系统偏差模型的贝叶斯理论统计公式具体如下所示：

x_k,s＝F_k-1x_k-1,s+w_k-1

z_k,s＝H_kx_k,s+η_k,s+v_k,s

η_k,s＝η_k-1,s+n_k-1,s s＝1,2,3,...,S

其中，x_k,s是k时刻传感器s的目标状态,k＝1,2,3...,ts，

(x_k,y_k)为位置坐标，

表示对应的速度，F_k-1为k-1时刻到k时刻状态转移矩阵，w_k-1为k-1时刻过程噪声,其噪声协方差矩阵为Q_k-1；z_k,s为k时刻传感器s的量测值，

H_k为k时刻量测转移矩阵；v_k,s为k时刻传感器s的量测噪声，其噪声协方差矩阵为R_k,s；η_k,s为k时刻传感器s的系统偏差，n_k-1,s为k-1时刻传感器s的系统偏差噪声，其噪声协方差矩阵为Q_η,k-1。

所述步骤S2变分参数设置及初始化期望，具体如下所示：

S2.1：设置k-1时刻变分参数初值如下：

初始估计状态及其误差协方差矩阵

P_k-1|k-1,s；量测噪声协方差矩阵R_k-1,s的自由度参数

和尺度矩阵为

自由度参数ω_k-1,ν_k-1的形状参数

和速率参数

和

Skew噪声发生概率π_k-1的形状参数

系统偏差η_k-1,s的均值

及其协方差矩阵Σ_k-1|k-1,s；标称即非真实过程噪声协方差矩阵

和标称量测噪声协方差矩阵

形状参数β_k-1的均值及其协方差矩阵

P_g,k-1|k-1；k时刻传感器s的量测值z_k,s；

S2.2：k时刻参数时间更新：

参数时间更新如下：状态一步预测：

及其预测状态误差协方差矩阵：

P_k|k-1,s的自由度参数

和尺度矩阵T_k|k-1,s＝τ_pP_k|k-1,s；量测噪声协方差矩阵R_k,s的自由度参数

和尺度矩阵U_k|k-1,s＝ρU_k-1|k-1,s；自由度ω_k、ν_k的形状参数

和速率参数

Skew噪声发生概率π_k的形状参数

系统偏差的均值

及其协方差矩阵Σ_k|k-1,s＝(Σ_k-1|k-1,s+Q_η,k-1)/ρ；

其中，(·)^T表示求转置操作，ρ表示变分遗忘因子，起自适应调节参数的作用，取1-exp(-4)，τ_p表示预测误差协方差的调谐因子；

S2.3：初始化k时刻期望值：

初始期望状态及其协方差矩阵：

初始化如下期望值：辅助变量λ_k、γ_k期望E⁰[λ_k]＝E⁰[γ_k]＝1；自由度参数ω_k、ν_k的期望E⁰[ω_k]、E⁰[ν_k]；伯努利变量ξ_k的期望E⁰[ξ_k]＝1；系统偏差及其协方差的期望：E⁰[η_k]、E⁰[P_η,k]；上标E⁰[·]表示求第0次变分迭代的期望。

所述步骤S3中建立P(x_k,s∣z_1:k-1)和P(z_k,s∣x_k,s)先验模型，具体包含以下步骤：

S3.1：建立状态一步预测概率密度P(x_k,s∣z_1:k-1)先验模型：

在实际工程应用中，由脉冲干扰、异常值和建模伪影引起的多种非高斯噪声通常具有Student-t分布或偏斜(Skew)分布，若状态过程噪声w_k-1遵循Skew分布，其概率密度函数P(w_k-1)被表示为高斯-伽马混合分布的形式，具体如下所示：

其中，N(·；μ,Σ)表示高斯分布，其均值μ，协方差矩阵为Σ，G(·；a,b)表示伽马分布，其形状参数a，速率参数b，s(·)和δ(·)分别是正偏斜和尺度函数，λ_k＞0为混合参数，ω_k为自由度参数，形状参数β_k控制Skew分布的对称性和偏度；β_k≠0是非对称的，若s(λ_k)＝δ(λ_k)＝λ_k为典型Skew分布；β_k＝0是对称的，若δ(λ_k)＝λ_k，Skew分布退化为典型Student-t分布；

为了判别k时刻目标状态x_k服从Student-t分布还是Skew分布，引入二元伯努利随机变量ξ_k，根据未知参数所服从的特性和步骤S1贝叶斯统计理论公式，对状态预测协方差矩阵P_k|k-1,s、随机变量ξ_k、Skew噪声发生概率π_k以及自由度参数ω_k建立先验模型以及结合过程噪声概率密度函数P(w_k-1)，则P(x_k,s∣z_1:k-1)被表示为：

其中，IW(·；t,T)表示逆威沙特分布，其自由度为t，尺度矩阵为T；Bn(·,π_k)表示伯努利分布，其取值为1的概率为π_k；Be(·；a,b)表示贝塔分布，其形状参数为a和b，常作为二元伯努利分布的共轭先验；ξ_k取值为0和1，ξ_k＝1表示状态x_k服从Skew分布，否则ξ_k＝0表示状态x_k服从Student-t分布；

S3.2：建立量测似然概率密度P(z_k,s∣x_k,s)先验模型：

考虑量测噪声为异常重尾噪声，以Student-t分布建模，引入辅助变量γ_k写成高斯-伽马混合分布形式，根据未知参数所服从的特性和步骤S1贝叶斯统计理论公式，对量测噪声协方差矩阵R_k,s、系统偏差η_k,s和自由度参数ν_k建立先验模型，则P(z_k,s∣x_k,s)被表示为：

S3.3：构建对数联合概率密度模型

根据步骤S3.1建立的状态一步概率密度P(x_k,s∣z_1:k-1)先验模型和步骤S3.2建立量测似然概率密度P(z_k,s∣x_k,s)先验模型，重写成分层高斯的形式，则对数联合概率密度

被表示为：

其中，

表示未知参数集合，Z表示传感器s的量测数据集，

表示以自然常数为底的对数联合概率密度。

所述步骤S6中更新求解P_k|k-1,s和R_k,s，具体如下所示：

分别令

应用步骤S4变分迭代公式，由指数分布共轭性，则近似后验q(P_k|k-1,s)和q(R_k,s)分别被更新为IW分布：

分别求解近似后验qⁱ⁺¹(P_k|k-1,s)和qⁱ⁺¹(R_k,s)的自由度参数：

和尺度矩阵参数：

由IW分布性质计算期望Eⁱ⁺¹[P_k|k-1,s]，Eⁱ⁺¹[R_k,s]，计算修正的状态预测协方差矩阵

和量测噪声协方差矩阵

其中，Eⁱ[·]表示求第i次变分迭代的期望，上标i表示第i次变分迭代，

所述步骤S7中更新求解β_k、η_k,s和x_k,s，具体如下所示：

分别令

和

应用步骤S4变分迭代公式，由指数分布共轭性，近似后验q(β_k)、q(η_k,s)和q(x_k,s)分别被更新为高斯分布：

类似卡尔曼滤波量测更新步分别求解近似后验qⁱ⁺¹(β_k)和qⁱ⁺¹(η_k,s)均值：

及其对应的协方差矩阵

由高斯分布性质，分别计算期望：Eⁱ⁺¹[β_k]、E^{i +1}[P_β,k]；

计算期望Eⁱ⁺¹[η_k,s]和修正系统偏差协方差矩阵期望

由卡尔曼滤波量测更新步计算

和

所述步骤S8中更新求解λ_k、γ_k、ω_k、ν_k，具体如下所示：

S8.1：更新辅助变量λ_k和γ_k：

令

应用步骤S4变分迭代公式，近似后验被建模成广义逆高斯分布GIG,即：

求解近似后验qⁱ⁺¹(λ_k)的形状参数

和

由现存标准GIG分布求解公式计算期望Eⁱ⁺¹[λ_k]、Eⁱ⁺¹[logλ_k]；

令

应用步骤S4变分迭代公式，由指数分布共轭性，近似后验q(γ_k)被更新为伽马分布：

求解近似后验qⁱ⁺¹(γ_k)的形状参数和速率参数：

由伽马分布性质计算期望Eⁱ⁺¹[γ_k]、Eⁱ⁺¹[logγ_k]；

S8.2：更新自由度参数ω_k和ν_k：

分别令

和

应用步骤S4变分迭代公式，由指数分布共轭性，近似后验q(ω_k)和q(ν_k)被更新为伽马分布：

求解近似后验qⁱ⁺¹(ω_k)和qⁱ⁺¹(ν_k)的形状参数：

和速率参数：

由伽马分布性质计算期望Eⁱ⁺¹[ω_k]和Eⁱ⁺¹[v_k]。

所述步骤S9中更新求解ξ_k和π_k，具体如下所示：

S9.1：更新随机变量ξ_k：

令

应用步骤S4变分迭代公式，则近似后验q(ξ_k)被更新为伯努利分布，分别计算ξ_k取0的概率Prⁱ⁺¹(ξ_k＝0)和ξ_k取1的Prⁱ⁺¹(ξ_k＝1)，由Bn分布性质，计算期望Eⁱ⁺¹[ξ_k]；

S9.2：更新Skew噪声发生概率π_k：

令

应用步骤S4变分迭代公式，由指数分布共轭性，q(π_k)被更新为贝塔分布：

求解近似后验qⁱ⁺¹(π_k)的形状参数：

由Be分布性质，计算自然对数期望Eⁱ⁺¹[logπ_k]、Eⁱ⁺¹[log(1-π_k)]。

所述步骤S12传感器分布式融合反馈具体包含：

首先，融合中心接收到各局部滤波器的状态估计，按照协方差交叉融合准则，其融合计算方法如下所示：

其中，

表示融合后的最优状态估计值及其对应的协方差矩阵；

其次，将状态及其协方差估计值反馈给各个局部滤波器，其反馈计算如下：

其中，κ_k,s为反馈权重系数，随局部状态估计协方差矩阵自适应调节，满足：

其中，||·||_F表示Frobenius范数；即对于任意矩阵D，

本发明的有益效果为：通过上述技术方案，本发明针对现有的问题在分布式融合框架下提出了一种新的未知概率Skew及重尾噪声下的目标跟踪方法：

首先，引入变分贝叶斯方法，考虑目标自身受过程噪声为未知概率随机出现的重尾噪声或Skew噪声的影响，更加适应实际工程应用中外界环境的不确定性；

其次，考虑了异常重尾量测噪声和未知且时变系统偏差对目标跟踪过程中的双重影响，即修正了传感器的量测信息值，间接提高了目标跟踪的精度；

最后，引入多传感器分布式融合反馈算法，在线融合各局部滤波器的目标状态及其对应的协方差矩阵并反馈给各局部滤波器，实时修正了目标状态及其对应的协方差矩阵，弥补了单一传感器的不足，使系统具有一定鲁棒性，提高了目标跟踪的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明包括以下步骤：

S1：考虑多传感器共同跟踪同一目标的物理场景，在目标跟踪过程中综合考虑目标状态过程噪声、传感器量测噪声和系统偏差的影响，故建立目标的状态模型、量测模型和系统偏差模型；

S2：设置k-1时刻变分迭代求解参数的初值和初始化期望，k＝1,2,3...,ts；ts为仿真总时间；

S3：建立k时刻传感器s的目标状态一步预测概率密度P(x_k,s∣z_1:k-1)模型和量测似然概率密度P(z_k,s∣x_k,s)模型和对数联合概率密度

模型；

S4：变分贝叶斯近似后验：

步骤S3构建的模型参数之间存在相互耦合，无法解析求解未知参数的后验，引入变分贝叶斯方法近似求解S3中的未知参数的后验，其求解公式如下：

其中，

表示关于

求期望操作，

表示取

集的某一个元素，

是

集中除

外剩余元素，

相对未知参数

是不相关的常数，

表示以自然常数为底的对数近似后验概率密度。

S5：设置变分迭代总次数为

i表示变分第i次迭代，

S6：更新求解状态预测协方差矩阵P_k|k-1,s和量测噪声协方差矩阵R_k,s；

S7：更新求解形状参数β_k、系统偏差η_k,s和目标状态x_k,s；

S8：更新求解辅助变量λ_k、γ_k及其对应自由度参数ω_k、ν_k；

S9：更新求解随机变量ξ_k和Skew噪声发生概率π_k；

S10：重复步骤S6-S9，直到达到预先设定变分迭代总次数或满足迭代终止条件:

ε一般取1e-15，则迭代结束；

S11：输出变分迭代结果:

输出迭代结果：

和Σ_k|k,s，并赋值作为k+1时刻变分迭代初始值；

S12：多传感器分布式融合反馈:

在分布式网络中，在其局部滤波器完成偏差配准后仅依赖单个传感器进行目标跟踪往往很难达到预期的精度，故需将步骤S11变分迭代求解得到的各局部滤波器状态值

以及其协方差矩阵P_k|k,s信息进行融合并将融合后的最优目标状态

及其协方差矩阵

反馈给各局部滤波器。

所述步骤S1中离散目标状态模型、量测模型和系统偏差模型的贝叶斯理论统计公式具体如下所示：

x_k,s＝F_k-1x_k-1,s+w_k-1

z_k,s＝H_kx_k,s+η_k,s+v_k,s

η_k,s＝η_k-1,s+n_k-1,s s＝1,2,3,...,S

其中，x_k,s是k时刻传感器s的目标状态,k＝1,2,3...,ts，

(x_k,y_k)为位置坐标，

表示对应的速度，F_k-1为k-1时刻到k时刻状态转移矩阵，w_k-1为k-1时刻过程噪声,其噪声协方差矩阵为Q_k-1；z_k,s为k时刻传感器s的量测值，

H_k为k时刻量测转移矩阵；v_k,s为k时刻传感器s的量测噪声，其噪声协方差矩阵为R_k,s；η_k,s为k时刻传感器s的系统偏差，n_k-1,s为k-1时刻传感器s的系统偏差噪声，其噪声协方差矩阵为Q_η,k-1。

所述步骤S2变分参数设置及初始化期望，具体如下所示：

S2.1：设置k-1时刻变分参数初值如下：

初始估计状态及其误差协方差矩阵

P_k-1|k-1,s；量测噪声协方差矩阵R_k-1,s的自由度参数

和尺度矩阵为

自由度参数ω_k-1,ν_k-1的形状参数

和速率参数

和

Skew噪声发生概率π_k-1的形状参数

系统偏差η_k-1,s的均值

及其协方差矩阵Σ_k-1|k-1,s；标称(非真实)过程噪声协方差矩阵

和标称量测噪声协方差矩阵

形状参数β_k-1的均值及其协方差矩阵

P_g,k-1|k-1。k时刻传感器s的量测值z_k,s等相关参数；

S2.2：k时刻参数时间更新：

参数时间更新如下：状态一步预测：

及其预测状态误差协方差矩阵：

P_k|k-1,s的自由度参数

和尺度矩阵T_k|k-1,s＝τ_pP_k|k-1,s；量测噪声协方差矩阵R_k,s的自由度参数

和尺度矩阵U_k|k-1,s＝ρU_k-1|k-1,s；自由度ω_k、ν_k的形状参数

和速率参数

Skew噪声发生概率π_k的形状参数

系统偏差的均值

及其协方差矩阵Σ_k|k-1,s＝(Σ_k-1|k-1,s+Q_η,k-1)/ρ。

其中，(·)^T表示求转置操作，ρ表示变分遗忘因子，起自适应调节参数的作用，通常取1-exp(-4)，τ_p表示预测误差协方差的调谐因子。

S2.3：初始化k时刻期望值：

初始期望状态及其协方差矩阵：

初始化如下期望值：辅助变量λ_k、γ_k期望E⁰[λ_k]＝E⁰[γ_k]＝1；自由度参数ω_k、ν_k的期望E⁰[ω_k]、E⁰[ν_k]；伯努利变量ξ^k的期望E⁰[ξ_k]＝1；系统偏差及其协方差的期望：E⁰[η_k]、E⁰[P_η,k]；上标E⁰[·]表示求第0次变分迭代的期望。

所述步骤S3中建立P(x_k,s∣z_1:k-1)和P(z_k,s∣x_k,s)先验模型，具体包含以下步骤：

S3.1：建立状态一步预测概率密度P(x_k,s∣z_1:k-1)先验模型：

在实际工程应用中，由脉冲干扰、异常值和建模伪影引起的多种非高斯噪声通常具有Student-t分布或偏斜(Skew)分布，若状态过程噪声w_k-1遵循Skew分布，其概率密度函数P(w_k-1)可以被表示为高斯-伽马混合分布的形式，具体如下所示：

其中，N(·；μ,Σ)表示高斯分布，其均值μ，协方差矩阵为Σ，G(·；a,b)表示伽马分布，其形状参数a，速率参数b，s(·)和δ(·)分别是正偏斜和尺度函数，λ_k＞0为混合参数，ω_k为自由度参数，形状参数β_k控制Skew分布的对称性和偏度。β_k≠0是非对称的，若s(λ_k)＝δ(λ_k)＝λ_k为典型Skew分布；β_k＝0是对称的，若δ(λ_k)＝λ_k，Skew分布退化为典型Student-t分布。

为了判别k时刻目标状态x_k服从Student-t分布还是Skew分布，引入二元伯努利随机变量ξ_k，根据未知参数所服从的特性和步骤S1贝叶斯统计理论公式，选择合适的分布分别对状态预测协方差矩阵P_k|k-1,s、随机变量ξ_k、Skew噪声发生概率π_k以及自由度参数ω_k建立先验模型以及结合过程噪声概率密度函数P(w_k-1)，则P(x_k,s∣z_1:k-1)可以被表示为：

其中，IW(·；t,T)表示逆威沙特分布，其自由度为^t，尺度矩阵为T；Bn(·,π_k)表示伯努利分布，其取值为1的概率为π_k；Be(·；a,b)表示贝塔分布，其形状参数为a和b，常作为二元伯努利分布的共轭先验；ξ_k取值为0和1，ξ_k＝1表示状态x_k服从Skew分布，否则ξ_k＝0表示状态x_k服从Student-t分布。

S3.2：建立量测似然概率密度P(z_k,s∣x_k,s)先验模型：

考虑量测噪声为异常重尾噪声，以Student-t分布建模，引入辅助变量γ_k写成高斯-伽马混合分布形式，根据未知参数所服从的特性和步骤S1贝叶斯统计理论公式，选择合适的分布对量测噪声协方差矩阵R_k,s、系统偏差η_k,s和自由度参数ν_k建立先验模型，则P(z_k,s∣x_k,s)可以被表示为：

S3.3：构建对数联合概率密度模型

根据步骤S3.1建立的状态一步概率密度P(x_k,s∣z_1:k-1)先验模型和步骤S3.2建立量测似然概率密度P(z_k,s∣x_k,s)先验模型，重写成分层高斯的形式，则对数联合概率密度

可以被表示为：

其中，

表示未知参数集合，Z表示传感器s的量测数据集，

表示以自然常数为底的对数联合概率密度。

所述步骤S6中更新求解P_k|k-1,s和R_k,s，具体如下所示：所述步骤S6中更新求解P_k|k-1,s和R_k,s，具体如下所示：

分别令

应用步骤S4变分迭代公式，由指数分布共轭性，则近似后验q(P_k|k-1,s)和q(R_k,s)分别被更新为IW分布：

S6.1：更新状态预测协方差矩阵P_k|k-1,s：

由相同项对应系数相等，则近似后验q(P_k|k-1)的自由度及其尺度矩阵

可以被计算为：

其中，i表示变分第i次变分迭代，辅助参数

辅助参数

则期望Eⁱ⁺¹[P_k|k-1]、Eⁱ[A_k]、Eⁱ[B_k]被计算为：

其中，Eⁱ[·]表示求第i次变分迭代的期望，

定义Eⁱ⁺¹[q_k]＝Eⁱ[β_k]/Eⁱ[λ_k]，修正的状态预测协方差矩阵

可以被计算为：

S6.2：更新量测噪声协方差矩阵R_k,s：

由相同项对应系数相等，则近似后验q(R_k,s)的自由度及其尺度矩阵

可以被计算为：

其中，辅助参数C_k＝(z_k,s-H_kx_k,s-η_k,s)(z_k,s-H_kx_k,s-η_k,s)^T，则期望Eⁱ⁺¹[R_k,s]和Eⁱ[C_k]被计算为：

其中，定义

修正量测噪声协方差矩阵

所述步骤S7中更新求解β_k、η_k,s和x_k,s，具体如下所示：

分别令

和

应用步骤S4变分迭代公式，由指数分布共轭性，近似后验q(β_k)、q(η_k,s)和q(x_k,s)分别被更新为高斯分布：

S7.1：更新形状参数β_k：

类似卡尔曼滤波的量测更新步，则近似后验qⁱ⁺¹(β_k)的均值及其协方差矩阵

可以被计算为：

期望Eⁱ⁺¹[β_k]、Eⁱ⁺¹[P_β,k]和

可以分别被计算为：

S7.2：更新系统偏差η_k,s：

类似卡尔曼滤波的量测更新步，则近似后验qⁱ⁺¹(η_k,s)的均值及其协方差矩阵

可以被计算为：

期望Eⁱ⁺¹[η_k,s]、Eⁱ⁺¹[P_η,k]可以分别被计算为：

S7.3：更新目标状态x_k,s：

根据卡尔曼滤波的量测更新步，则近似后验qⁱ⁺¹(x_k,s)的均值及其协方差矩阵

可以被计算为：

所述步骤S8中更新求解λ_k、γ_k、ω_k、ν_k，具体如下所示：

S8.1：更新辅助变量λ_k和γ_k：

令

应用步骤S4变分迭代公式，近似后验被建模成广义逆高斯分布(GIG),即：

其中，h_k,g_k,p_k是广义逆高斯的形状参数。则公式被更新为：

则Eⁱ⁺¹[λ_k]、Eⁱ⁺¹[logλ_k]的期望由下式近似计算：

其中，K_p(·)表示二阶(second kind)修正的贝塞尔函数，下标p表示索引。

令

应用步骤S4变分迭代公式，由指数分布共轭性，近似后验q(γ_k)被更新为伽马分布：

近似后验qⁱ⁺¹(γ_k)的形状参数和速率参数：

可以被计算为：

期望Eⁱ⁺¹[γ_k]、Eⁱ⁺¹[logγ_k]可以分别被计算为：

其中，ψ(·)表示digamma函数。

S8.2：更新自由度参数ω_k和ν_k：

分别令

和

应用步骤S4变分迭代公式，由指数分布共轭性，则近似后验q(ω_k)和q(ν_k)被更新为伽马分布：

由相同项对应系数相等，则近似后验qⁱ⁺¹(ω_k)和qⁱ⁺¹(ν_k)的形状参数：

和速率参数：

可以被计算为：

期望Eⁱ⁺¹[ω_k]、Eⁱ⁺¹[ν_k]可以被计算为：

所述步骤S9中更新求解ξ_k和π_k，具体如下所示：

S9.1：更新随机变量ξ_k：

令

应用步骤S4变分迭代公式，由指数族分布的共轭性，则近似后验q(ξ_k)被更新为伯努利分布，ξ_k取0和1的概率可以被计算为：

Prⁱ⁺¹(ξ_k＝0)＝Λⁱ⁺¹exp{Eⁱ[log(1-π_k)]

+0.5nEⁱ⁺¹[logλ_k]-0.5Eⁱ⁺¹[λ_k]Eⁱ[A_k](Eⁱ[P_k|k-1,s])^-1}

Prⁱ⁺¹(ξ_k＝1)＝Λⁱ⁺¹exp{Eⁱ[logπ_k]

+0.5nEⁱ⁺¹[logλ_k]-0.5Eⁱ⁺¹[λ_k]Eⁱ[B_k](Eⁱ[P_k|k-1,s])^-1}

其中，exp(·)表示以自然常数e为底的指数函数，辅助参数期望Eⁱ[A_k]、Eⁱ[B_k]在步骤S6.1被给出，

是归一化常数，期望Eⁱ⁺¹[ξ_k]可以计算为：

S9.2：更新Skew噪声发生概率π_k：

令

应用步骤S4变分迭代公式，由指数分布共轭性，q(π_k)被更新为贝塔分布：

由相同项对应系数相等，则形状参数

为：

期望Eⁱ⁺¹[logπ_k]、Eⁱ⁺¹[log(1-π_k)]可以被计算为：

所述步骤S12传感器分布式融合反馈具体包含：

首先，融合中心接收到各局部滤波器的状态估计，按照协方差交叉融合准则，其融合计算方法如下所示：

其中，

表示融合后的最优状态估计值及其对应的协方差矩阵。

其次，将状态及其协方差估计值反馈给各个局部滤波器，其反馈计算如下：

其中，κ_k,s为反馈权重系数，随局部状态估计协方差矩阵自适应调节，满足：

其中，||·||_F表示Frobenius范数。即对于任意矩阵D，

本发明通过上述技术方案，首先，通过引入变分贝叶斯方法，考虑目标自身受过程噪声为未知概率随机出现的重尾噪声或Skew噪声的影响，更加适应实际工程应用中外界环境的不确定性；其次，考虑了异常重尾量测噪声和未知且时变系统偏差对目标跟踪过程中的双重影响，即修正了传感器的量测信息值，间接提高了目标跟踪的精度；最后，引入多传感器分布式融合反馈算法，在线融合各局部滤波器的目标状态及其对应的协方差矩阵并反馈给各局部滤波器，实时修正了目标状态及其对应的协方差矩阵，弥补了单一传感器的不足，使系统具有一定鲁棒性，提高了目标跟踪的精度。

本发明在分布式融合框架下提出了一种新的未知概率Skew及重尾噪声下的目标跟踪方法，在步骤S2.2引入遗忘因子ρ，起自适应调节参数的作用，考虑多个传感器共同跟踪同一个目标的物理场景，不仅在步骤S6、S7和S9辨识并解决了过程噪声为未知概率随机出现的重尾噪声或Skew噪声；还在步骤S6、S7解决了量测噪声为异常重尾噪声以及系统偏差未知且时变的问题；以及在步骤S12提出分布式传感器融合反馈算法，不仅减轻了传感器的通信负担，还提高了目标跟踪精度。

在本发明的描述中，需要说明的是，对于方位词，如有术语“中心”，“横向”、“纵向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示方位和位置关系为基于附图所示的方位或位置关系，仅是为了便于叙述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定方位构造和操作，不能理解为限制本发明的具体保护范围。

需要说明的是，本申请的说明书和权利要求书中的术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

注意，上述仅为本发明的较佳实施例及运用技术原理。本领域技术人员会理解，本发明不限于这里所述的特定实施例，对本领域技术人员来说能够进行各种明显的变化、重新调整和替代而不会脱离本发明的保护范围。因此，虽然通过以上实施例对本发明进行较详细的说明，但本发明不限于这里所述的特定实施例，在不脱离本发明构思的情况下，还可以包括更多其他等有效实施例，而本发明的范围由所附的权利要求范围决定。

Claims (9)

1.未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：包括

包括以下步骤：

S1：多传感器共同跟踪同一目标的物理场景，在目标跟踪过程中综合考虑目标状态过程噪声、传感器量测噪声和系统偏差的影响，故建立目标的状态模型、量测模型和系统偏差模型；

S2：设置k-1时刻变分迭代求解参数的初值和初始化期望，k＝1,2,3...,ts，ts为仿真总时间；

S3：建立k时刻传感器s的目标状态一步预测概率密度P(x_k,s∣z_1:k-1)模型、量测似然概率密度P(z_k,s∣x_k,s)模型和对数联合概率密度

模型；

S4：变分贝叶斯近似后验：

步骤S3构建的模型参数之间存在相互耦合，无法解析求解未知参数的后验，引入变分贝叶斯方法近似求解S3中的未知参数的后验，其求解公式如下：

其中，

表示关于

求期望操作，

表示取

集的某一个元素，

是

集中除

外剩余元素，

相对未知参数

是不相关的常数，

表示以自然常数为底的对数近似后验概率密度；

S5：设置变分迭代总次数为

i表示变分第i次迭代，i＝0:

S6：更新求解状态预测协方差矩阵P_k|k-1,s和量测噪声协方差矩阵R_k,s；

S7：更新求解形状参数β_k、系统偏差η_k,s和目标状态x_k,s；

S8：更新求解辅助变量λ_k、γ_k及其对应自由度参数ω_k、ν_k；

S9：更新求解随机变量ξ_k和Skew噪声发生概率π_k；

S10：重复步骤S6-S9，直到达到预先设定变分迭代总次数或满足迭代终止条件:

ε一般取1e-15，则迭代结束；

S11：输出变分迭代结果:

输出迭代结果：

和Σ_k|k,s，并赋值作为k+1时刻变分迭代初始值；

S12：多传感器分布式融合反馈:

在分布式网络中，在其局部滤波器完成偏差配准后仅依赖单个传感器进行目标跟踪往往很难达到预期的精度，故需将步骤S11变分迭代求解得到的各局部滤波器状态值

以及其协方差矩阵P_k|k,s信息进行融合并将融合后的最优目标状态

及其协方差矩阵

反馈给各局部滤波器。

2.根据权利要求1所述的未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：所述步骤S1中离散目标状态模型、量测模型和系统偏差模型的贝叶斯理论统计公式具体如下所示：

x_k,s＝F_k-1x_k-1,s+w_k-1

z_k,s＝H_kx_k,s+η_k,s+v_k,s

η_k,s＝η_k-1,s+n_k-1,s s＝1,2,3,...,S

其中，x_k,s是k时刻传感器s的目标状态,k＝1,2,3...,ts，

(x_k,y_k)

为位置坐标，

表示对应的速度，F_k-1为k-1时刻到k时刻状态转移矩阵，w_k-1为k-1时刻过程噪声,其噪声协方差矩阵为Q_k-1；z_k,s为k时刻传感器s的量测值，

H_k为k时刻量测转移矩阵；v_k,s为k时刻传感器s的量测噪声，其噪声协方差矩阵为R_k,s；η_k,s为k时刻传感器s的系统偏差，n_k-1,s为k-1时刻传感器s的系统偏差噪声，其噪声协方差矩阵为Q_η,k-1。

3.根据权利要求1所述的未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：所述步骤S2变分参数设置及初始化期望，具体如下所示：

S2.1：设置k-1时刻变分参数初值如下：

初始估计状态及其误差协方差矩阵

P_k-1|k-1,s；量测噪声协方差矩阵R_k-1,s的自由度参数

和尺度矩阵为

自由度参数ω_k-1,ν_k-1的形状参数

和速率参数

和

Skew噪声发生概率π_k-1的形状参数

系统偏差η_k-1,s的均值

及其协方差矩阵Σ_k-1|k-1,s；标称即非真实过程噪声协方差矩阵

和标称量测噪声协方差矩阵

形状参数β_k-1的均值及其协方差矩阵

P_g,k-1|k-1；k时刻传感器s的量测值z_k,s；

S2.2：k时刻参数时间更新：

参数时间更新如下：状态一步预测：

及其预测状态误差协方差矩阵：

P_k|k-1,s的自由度参数

和尺度矩阵T_k|k-1,s＝τ_pP_k|k-1,s；量测噪声协方差矩阵R_k,s的自由度参数

和尺度矩阵U_k|k-1,s＝ρU_k-1|k-1,s；自由度ω_k、ν_k的形状参数

和速率参数

Skew噪声发生概率π_k的形状参数

系统偏差的均值

及其协方差矩阵Σ_k|k-1,s＝(Σ_k-1|k-1,s+Q_η,k-1)/ρ；

其中，(·)^T表示求转置操作，ρ表示变分遗忘因子，起自适应调节参数的作用，取1-exp(-4)，τ_p表示预测误差协方差的调谐因子；

S2.3：初始化k时刻期望值：

初始期望状态及其协方差矩阵：

初始化如下期望值：辅助变量λ_k、γ_k期望E⁰[λ_k]＝E⁰[γ_k]＝1；自由度参数ω_k、ν_k的期望E⁰[ω_k]、E⁰[ν_k]；伯努利变量ξ_k的期望E⁰[ξ_k]＝1；系统偏差及其协方差的期望：E⁰[η_k]、E⁰[P_η,k]；上标E⁰[·]表示求第0次变分迭代的期望。

4.根据权利要求1所述的未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：所述步骤S3中建立P(x_k,s∣z_1:k-1)和P(z_k,s∣x_k,s)先验模型，具体包含以下步骤：

S3.1：建立状态一步预测概率密度P(x_k,s∣z_1:k-1)先验模型：

在实际工程应用中，由脉冲干扰、异常值和建模伪影引起的多种非高斯噪声通常具有Student-t分布或偏斜(Skew)分布，若状态过程噪声w_k-1遵循Skew分布，

其概率密度函数P(w_k-1)被表示为高斯-伽马混合分布的形式，具体如下所示：

其中，N(·；μ,Σ)表示高斯分布，其均值μ，协方差矩阵为Σ，G(·；a,b)表示伽马分布，其形状参数a，速率参数b，s(·)和δ(·)分别是正偏斜和尺度函数，λ_k＞0为混合参数，ω_k为自由度参数，形状参数β_k控制Skew分布的对称性和偏度；β_k≠0是非对称的，若s(λ_k)＝δ(λ_k)＝λ_k为典型Skew分布；β_k＝0是对称的，若δ(λ_k)＝λ_k，Skew分布退化为典型Student-t分布；

为了判别k时刻目标状态x_k服从Student-t分布还是Skew分布，引入二元伯努利随机变量ξ_k，根据未知参数所服从的特性和步骤S1贝叶斯统计理论公式，对状态预测协方差矩阵P_k|k-1,s、随机变量ξ_k、Skew噪声发生概率π_k以及自由度参数ω_k建立先验模型以及结合过程噪声概率密度函数P(w_k-1)，则P(x_k,s∣z_1:k-1)被表示为：

其中，IW(·；t,T)表示逆威沙特分布，其自由度为t，尺度矩阵为T；Bn(·,π_k)表示伯努利分布，其取值为1的概率为π_k；Be(·；a,b)表示贝塔分布，其形状参数为a和b，常作为二元伯努利分布的共轭先验；ξ_k取值为0和1，ξ_k＝1表示状态x_k服从Skew分布，否则ξ_k＝0表示状态x_k服从Student-t分布；

S3.2：建立量测似然概率密度P(z_k,s∣x_k,s)先验模型：

考虑量测噪声为异常重尾噪声，以Student-t分布建模，引入辅助变量γ_k写成高斯-伽马混合分布形式，根据未知参数所服从的特性和步骤S1贝叶斯统计理论公式，对量测噪声协方差矩阵R_k,s、系统偏差η_k,s和自由度参数ν_k建立先验模型，则P(z_k,s∣x_k,s)被表示为：

S3.3：构建对数联合概率密度模型

根据步骤S3.1建立的状态一步概率密度P(x_k,s∣z_1:k-1)先验模型和步骤S3.2建立量测似然概率密度P(z_k,s∣x_k,s)先验模型，重写成分层高斯的形式，则对数联合概率密度

被表示为：

其中，

表示未知参数集合，Z表示传感器s的量测数据集，

表示以自然常数为底的对数联合概率密度。

5.根据权利要求1所述的未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：所述步骤S6中更新求解P_k|k-1,s和R_k,s，具体如下所示：

分别令

应用步骤S4变分迭代公式，由指数分布共轭性，则近似后验q(P_k|k-1,s)和q(R_k,s)分别被更新为IW分布：

分别求解近似后验qⁱ⁺¹(P_k|k-1,s)和qⁱ⁺¹(R_k,s)的自由度参数：

和尺度矩阵参数：

由IW分布性质计算期望Eⁱ⁺¹[P_k|k-1,s]，Eⁱ⁺¹[R_k,s]，计算修正的状态预测协方差矩阵

和量测噪声协方差矩阵

其中，Eⁱ[·]表示求第i次变分迭代的期望，上标i表示第i次变分迭代，i＝0:

6.根据权利要求1所述的未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：所述步骤S7中更新求解β_k、η_k,s和x_k,s，具体如下所示：

分别令

和

应用步骤S4变分迭代公式，由指数分布共轭性，

近似后验q(β_k)、q(η_k,s)和q(x_k,s)分别被更新为高斯分布：

类似卡尔曼滤波量测更新步分别求解近似后验qⁱ⁺¹(β_k)和qⁱ⁺¹(η_k,s)均值：

及其对应的协方差矩阵

由高斯分布性质，分别计算期望：Eⁱ⁺¹[β_k]、Eⁱ⁺¹[P_β,k]；

计算期望Eⁱ⁺¹[η_k,s]和修正系统偏差协方差矩阵期望

由卡尔曼滤波量测更新步计算

和

7.根据权利要求1所述的未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：所述步骤S8中更新求解λ_k、γ_k、ω_k、ν_k，具体如下所示：

S8.1：更新辅助变量λ_k和γ_k：

令

应用步骤S4变分迭代公式，近似后验被建模成广义逆高斯分布GIG,即：

求解近似后验qⁱ⁺¹(λ_k)的形状参数

和

由现存标准GIG分布求解公式计算期望Eⁱ⁺¹[λ_k]、Eⁱ⁺¹[logλ_k]；

令

应用步骤S4变分迭代公式，由指数分布共轭性，近似后验q(γ_k)被更新为伽马分布：

求解近似后验qⁱ⁺¹(γ_k)的形状参数和速率参数：

由伽马分布性质计算期望Eⁱ⁺¹[γ_k]、Eⁱ⁺¹[logγ_k]；

S8.2：更新自由度参数ω_k和ν_k：

分别令

和

应用步骤S4变分迭代公式，由指数分布共轭性，近似后验q(ω_k)

和q(ν_k)被更新为伽马分布：

求解近似后验qⁱ⁺¹(ω_k)和qⁱ⁺¹(ν_k)的形状参数：

和速率参数：

由伽马分布性质计算期望Eⁱ⁺¹[ω_k]和Eⁱ⁺¹[v_k]。

8.根据权利要求1所述的未知概率Skew和重尾噪声下基于变分贝叶斯的目标跟踪方法，其特征在于：所述步骤S9中更新求解ξ_k和π_k，具体如下所示：

S9.1：更新随机变量ξ_k：

令

应用步骤S4变分迭代公式，则近似后验q(ξ_k)被更新为伯努利分布，分别计算ξ_k取0的概率Prⁱ⁺¹(ξ_k＝0)和ξ_k取1的Prⁱ⁺¹(ξ_k＝1)，由Bn分布性质，计算期望Eⁱ⁺¹[ξ_k]；

S9.2：更新Skew噪声发生概率π_k：

令

应用步骤S4变分迭代公式，由指数分布共轭性，q(π_k)被更新为贝塔分布：

求解近似后验qⁱ⁺¹(π_k)的形状参数：

由Be分布性质，计算自然对数期望Eⁱ⁺¹[logπ_k]、Eⁱ⁺¹[log(1-π_k)]。

9.根据权利要求1所述的未知概率Skew及重尾噪声下的目标跟踪方法，其特征在于：所述步骤S12传感器分布式融合反馈具体包含：

首先，融合中心接收到各局部滤波器的状态估计，按照协方差交叉融合准则，其融合计算方法如下所示：

其中，

表示融合后的最优状态估计值及其对应的协方差矩阵；

其次，将状态及其协方差估计值反馈给各个局部滤波器，其反馈计算如下：

其中，κ_k,s为反馈权重系数，随局部状态估计协方差矩阵自适应调节，满足：

其中，||·||_F表示Frobenius范数；即对于任意矩阵D，

Images