一种基于变分贝叶斯学习的井下传感器网络目标跟踪方法

1.本发明涉及煤矿井下安全领域，具体涉及一种基于变分贝叶斯学习的井下传感器网络目标跟踪方法。


背景技术：

2.由于井下无线信号传播是在封闭、不规律的空间内进行，这将造成无线信号传播过程中发生反射、折射和绕射，因此最后接收到的信号通常为多路径的合成波。此外，井下工作环境复杂，各种生产设备会对传感器的无线探测信号造成干扰，这些因素导致井下传感器网络进行目标跟踪模型中过程噪声和量测噪声服从厚尾非高斯分布，且噪声参数未知，这些都会影响井下目标跟踪定位结果的精确度。
3.现有针对非高斯噪声下的目标跟踪方法有高斯求和滤波、粒子滤波、基于m估计的滤波方法和极大熵滤波等方法。高斯求和滤波和粒子滤波算法由于计算复杂度高，在低成本的传感器网络中难以直接应用。基于m估计的滤波方法和极大熵滤波方法本质是一种野值剔除的方法，存在自由参数选择缺乏理论依据且估计精度低的问题。
4.公开号cn103778320a的中国专利公开了一种基于变分贝叶斯多传感器量化融合目标跟踪方法，该方法使用变分贝叶斯方法进行自适应强跟踪信息滤波，对量测噪声的未知方差进行估计，但是该方法仅考虑过程噪声与量测噪声近似为零均值高斯白噪声的情形。公开号cn105682220a的中国专利公开了一种基于无线传感器网络的矿山井下定位系统，该系统使用最小二乘法对rssi测距环境参数进行修正，消除各种干扰对测量数据的影响，但是该方法仅在rssi信号强度值和移动节点到接收端节点的距离值同时存在均值为零、方差相同的噪声误差时，才能得到最好的参数拟合结果。两者均没有考虑到井下受到多路径效应以及电磁干扰等引起过程噪声与量测噪声不再服从高斯分布的问题，建模不准确，因此在井下目标跟踪中会存在精度下降甚至算法发散的问题。


技术实现要素：

5.针对上述存在的技术不足，本发明的目的是提供一种基于变分贝叶斯学习的井下传感器网络目标跟踪方法，其能够在矿井下复杂的工作环境中对工作人员和设备进行精准定位，保障井下安全生产，减少人员伤亡，助力智慧矿山建设。
6.为解决上述技术问题，本发明采用如下技术方案：
7.本发明提供一种基于变分贝叶斯学习的井下传感器网络目标跟踪方法，包括以下步骤：
8.步骤1、建立模型，构建井下目标跟踪问题的状态空间模型；
9.步骤2、时间更新，得到系统状态当前时刻的预测值；
10.步骤3、量测更新，基于变分贝叶斯学习方法迭代计算系统状态变量与噪声参数的后验估计；
11.步骤4、输出结果，结束循环迭代，利用状态估计结果得到目标的位置等信息。
12.优选地，步骤1中，针对在井下封闭环境下，目标跟踪问题中过程噪声和量测噪声服从厚尾非高斯分布，采用椭球分布对过程噪声和量测噪声进行建模，利用inverse-wishart分布建模过程噪声和量测噪声散度矩阵的参数不确定性，具体方法如下：
13.步骤1-1：给出线性系统的标准状态空间模型：
14.xi＝fix
i-1
+ωi15.z
i,k
＝h
i,k
xi+ν
i,k
16.xi表示在i时刻的系统状态，包含目标的位置、速度、加速度信息，z
i,k
表示在i时刻的传感器1
…
k的测量值，fi和h
i,k
表示已知的状态转移矩阵和量测矩阵，过程噪声ωi与量测噪声υ
i,k
相互独立，不同传感器之间的量测噪声相互独立；
17.步骤1-2：为构建厚尾非高斯过程和量测噪声模型，假设ωi和υ
i,k
服从椭球分布：
18.p(ωi)＝ε(ωi；αi,qi,ψi)
19.p(υ
i,k
)＝ε(υ
i,k
；β
i,k
,r
i,k
,ζ
i,k
)
20.式中αi与β
i,k
分别为过程噪声和量测噪声分布均值，qi和r
i,k
分别为过程噪声和量测噪声散度矩阵，ψi与ζ
i,k
分别为过程噪声和量测噪声密度发生函数；
21.步骤1-3：利用椭球分布的条件分布为高斯分布进行简化：
22.p(ωi|λi)＝n(0,g(λi)-1
qi)
23.p(υ
i,k
|γ
i,k
)＝n(0,g(λi)-1ri,k
)
24.式中λi与γ
i,k
分别为过程噪声和量测噪声的尺度变量，g(λi)-1
和g(γ
i,k
)-1
分别是λi和γ
i,k
的函数；
25.步骤1-4：定义尺度变量λi和γ
i,k
的先验分布：
[0026][0027][0028]
式中和分别是λi和γ
i,k
的参数集合；
[0029]
步骤1-5：利用inverse-wishart分布建模过程噪声散度矩阵qi和量测噪声散度矩阵r
i,k
：
[0030]
p(qi)＝iw(qi；t0,t0)
[0031]
p(r
i,k
)＝iw(r
i,k
；v0,v0)
[0032]
式中t0，v0，t0，v0是p(qi)和p(r
i,k
)的自由度参数和逆尺度矩阵。
[0033]
优选地，步骤2中，根据i-1时刻的状态估计预测i时刻状态估计，具体方法如下：
[0034]
步骤2-1：系统状态xi，其期望值为相应的估计协方差矩阵为∑i，输入初值∑
i-1|i-1
、q0、t0、t0、v0、v0；
[0035]
步骤2-2：时间更新
[0036][0037]
∑
i|i-1
＝fi∑
i-1|i-1fit
+q0。
[0038]
式中为在i时刻的一步预测状态期望值，∑
i|i-1
为相应的的一步预测估计误差协方差矩阵，为在i-1时刻的估计状态期望值，∑
i-1|i-1
为相应的估计误差协方差矩
阵。
[0039]
优选地，步骤3中，通过变分贝叶斯方法得到系统状态变量与噪声参数的后验估计，通过循环迭代进行状态更新，具体方法如下：
[0040]
步骤3-1：对模型参数进行初始化：
[0041][0042]
步骤3-2：进行n-1次循环迭代；
[0043]
步骤3-3：计算期望值：
[0044]
参数的期望值为：
[0045][0046]
式中与分别为q
(l+1)
(qi)的自由度参数和逆尺度矩阵，n
x
为系统变量xi的维数；
[0047]
参数的期望值为：
[0048][0049]
参数λi的期望值e
(l)
(λi)和参数γ
i,k
的期望值e
(l)
(γ
i,k
)根据其具体分布确定；
[0050]
步骤3-4：更新其中和分别由以下两式求得：
[0051][0052][0053]
式中l为固定点迭代次数，为估计状态的误差协方差矩阵，e
(l)
[g(γ
i,k
)]、e
(l)
[μi]、e
(l)
[g(λi)]分别为g(γ
i,k
)、g(λi)、μi的期望值；
[0054]
步骤3-5：更新其中和分别由以下两式求得：
[0055][0056][0057]
式中为μi的期望值，为的误差协方差矩阵，e
(l+1)
[xi]为系统状态xi的期望值；
[0058]
步骤3-6：更新q
(l+1)
(λi)和q
(l+1)
(γ
i,k
)，其具体分布根据选定的先验分布借助变分贝叶斯学习确定；
[0059]
步骤3-7：更新其中和分别由以下两式确定：
[0060][0061]
[0062]
式中为的辅助矩阵；
[0063]
步骤3-8：更新其中和分别由以下两式确定：
[0064][0065][0066]
式中为的辅助矩阵。
[0067]
优选地，步骤4中，在结束循环迭代后，输出最优估计状态和协方差矩阵
[0068]
本发明的有益效果在于：本发明提出的基于变分贝叶斯学习的井下传感器网络目标跟踪方法，能够在过程噪声和量测噪声为非高斯噪声并且噪声参数不精确的情况下实现精确跟踪定位，并且考虑了每个传感器噪声分布的差异，与现有技术相比，本发明具有建模准确、计算简单、且估计精度高的特点，能够实现井下等复杂环境下的传感器网络目标跟踪。
附图说明
[0069]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0070]
图1为本发明实施例提供的一种基于变分贝叶斯学习的井下传感器网络目标跟踪方法的步骤流程图；
[0071]
图2为本发明与现有变分贝叶斯方法、现有基于t分布的变分贝叶斯方法及经典卡尔曼滤波方法的位置平均均方根误差对比图；
[0072]
图3为本发明提出的方法与现有变分贝叶斯方法、现有基于t分布的变分贝叶斯方法及经典卡尔曼滤波方法的速度平均均方根误差对比图；
[0073]
图4为本发明的实际应用场景图。
具体实施方式
[0074]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0075]
如图1至图4所示，一种基于变分贝叶斯学习的井下传感器网络目标跟踪方法，具体步骤如下：
[0076]
步骤1、建立模型，构建井下目标跟踪问题的状态空间模型；
[0077]
针对在井下等封闭环境下目标跟踪问题中过程噪声和量测噪声服从厚尾非高斯分布，采用椭球分布对过程噪声和量测噪声进行建模，利用inverse-wishart分布建模过程噪声和量测噪声散度矩阵的参数不确定性，具体方法如下：
[0078]
步骤1-1：给出线性系统的标准状态空间模型：
[0079]
xi＝fix
i-1
+ωi[0080]zi,k
＝h
i,k
xi+ν
i,k
[0081]
式中xi表示在i时刻的系统状态，通常包含目标的位置、速度、加速度等信息，z
i,k
表示在i时刻的传感器1
…
k的测量值，fi和h
i,k
表示已知的状态转移矩阵和量测矩阵，过程噪声ωi与量测噪声υ
i,k
相互独立，不同传感器之间的量测噪声相互独立；
[0082]
步骤1-2：为构建厚尾非高斯过程和量测噪声模型，假设ωi和υ
i,k
服从椭球分布：
[0083]
p(ωi)＝ε(ωi；αi,qi,ψi)
[0084]
p(υ
i,k
)＝ε(υ
i,k
；β
i,k
,r
i,k
,ζ
i,k
)
[0085]
式中αi与β
i,k
分别为过程噪声和量测噪声分布均值，qi和r
i,k
分别为过程噪声和量测噪声散度矩阵，ψi与ζ
i,k
分别为过程噪声和量测噪声密度发生函数；
[0086]
步骤1-3：利用椭球分布的条件分布为高斯分布进行简化：
[0087]
p(ωi|λi)＝n(0,g(λi)-1
qi)
[0088]
p(υ
i,k
|γ
i,k
)＝n(0,g(λi)-1ri,k
)
[0089]
式中λi与γ
i,k
分别为过程噪声和量测噪声的尺度变量，g(λi)-1
和g(γ
i,k
)-1
分别是λi和γ
i,k
的函数；
[0090]
步骤1-4：定义尺度变量λi和γ
i,k
的先验分布：
[0091][0092][0093]
式中和分别是λi和γ
i,k
的参数集合；
[0094]
步骤1-5：利用inverse-wishart分布构建过程噪声协方差矩阵qi和量测噪声协方差矩阵r
i,k
：
[0095]
p(qi)＝iw(qi；t0,t0)
[0096]
p(r
i,k
)＝iw(r
i,k
；v0,v0)
[0097]
式中t0，v0，t0，v0是p(qi)和p(r
i,k
)的自由度参数和逆尺度矩阵。
[0098]
步骤2、时间更新，得到系统状态当前时刻的预测值；
[0099]
根据i-1时刻的状态估计预测i时刻状态估计，具体方法如下：
[0100]
步骤2-1：系统状态xi，其期望值为相应的估计协方差矩阵为∑i，输入初值∑
i-1|i-1
、q0、t0、t0、v0、v0；
[0101]
步骤2-2：时间更新
[0102][0103]
∑
i|i-1
＝fi∑
i-1|i-1fit
+q0[0104]
式中为在i时刻的一步预测状态期望值，∑
i|i-1
为相应的的一步预测估计误差协方差矩阵，为在i-1时刻的估计状态期望值，∑
i-1|i-1
为相应的估计误差协方差矩
阵。
[0105]
步骤3、量测更新，计算系统状态变量与其他参数的后验估计，进行状态更新；
[0106]
通过变分贝叶斯学习的方法得到系统状态变量与噪声参数的后验估计，通过循环迭代进行状态更新，具体方法如下：
[0107]
步骤3-1：对模型参数进行初始化：
[0108][0109]
步骤3-2：进行n-1次循环迭代；
[0110]
步骤3-3：计算期望值：
[0111]
参数的期望值为：
[0112][0113]
式中与分别为q
(l+1)
(qi)的自由度参数和逆尺度矩阵，n
x
为系统变量xi的维数；
[0114]
参数的期望值为：
[0115][0116]
式中与分别为q
(l+1)
(r
i,k
)的自由度参数和逆尺度矩阵，n
zk
为传感器测量值z
i,k
的维数；
[0117]
参数λi的期望值e
(l)
(λi)和参数γ
i,k
的期望值e
(l)
(γ
i,k
)根据其具体分布确定；
[0118]
步骤3-4：更新其中和分别由以下两式求得：
[0119][0120][0121]
式中l为固定点迭代次数，为估计状态的误差协方差矩阵，e
(l)
[g(γ
i,k
)]、e
(l)
[μi]、e
(l)
[g(λi)]分别为g(γ
i,k
)、g(λi)、μi的期望值；
[0122]
步骤3-5：更新其中和分别由以下两式求得：
[0123][0124][0125]
式中为μi的期望值，为的误差协方差矩阵，e
(l+1)
[xi]为系统状态xi的期望值；步骤3-6：更新q
(l+1)
(λi)和q
(l+1)
(γ
i,k
)，其具体分布根据选定的先验分布借助变分贝叶斯学习确定；
[0126]
步骤3-7：更新其中和分别由以下两式确定：
[0127]
[0128][0129]
式中为的辅助矩阵；
[0130]
步骤3-8：更新其中和分别由以下两式确定：
[0131][0132][0133]
式中为的辅助矩阵。
[0134]
步骤4、输出结果，结束循环迭代，利用状态估计结果得到目标的位置等信息；
[0135]
在结束循环迭代后，输出最优估计状态和协方差矩阵
[0136]
以下对本发明的效果做进一步的验证：
[0137]
选用二维平面基于多个测距、测速传感器构成的传感器网络进行人员跟踪的例子来验证本发明提出的方法的效果。
[0138]
选择状态转移矩阵为：
[0139][0140]
式中t＝1s是采样时间，i2是2
×
2单位矩阵，02是2
×
2零矩阵；
[0141]
传感器网络包含k＝1,2,...,10个传感器，当k＝1、3、5、7、9时，量测矩阵h
i,k
＝[i
2 02]；
[0142]
当k＝2、4、6、8、10时，量测矩阵h
i,k
＝[0
2 i2]。
[0143]
选择系统状态向量为：
[0144]
xi＝[p
x,i p
y,i v
x,i v
y,i
]
t
[0145]
式中(p
x,i
,p
y,i
)和(v
x,i
,v
y,i
)分别是在x和y方向上的位置和速度；
[0146]
选择初始状态为：
[0147][0148]
选择协方差矩阵为：
[0149]
p0＝diag[1m2,1m2,1m2s-2
,1m2s-2
]
t
[0150]
选用平均均方根误差(armse)计算本发明提出的方法、现有变分贝叶斯方法的滤波方法、现有基于t分布的变分贝叶斯滤波方法及经典卡尔曼滤波方法的估计误差。位置的平均均方根误差的计算公式为：
[0151][0152]
式中rmse
pos
(k)的计算公式为:
[0153]
[0154]
式中下标s是蒙特卡洛模拟实验的索引，是传感器k的检测到的位置的估计，蒙特卡洛模拟实验次数mc＝1000，速度的平均均方根误差armse
vel
的计算方式与上述方法类似。过程噪声和量测噪声均利用高斯混合分布来模拟典型的厚尾非高斯噪声，用来模拟多路径、野值、电磁干扰、目标机动等引起的噪声变化。
[0155]
图2、图3分别是采用本发明提出的方法、现有变分贝叶斯方法、现有基于t分布的变分贝叶斯方法及经典卡尔曼滤波对上述模型进行估计的位置平均均方根误差和速度平均均方根误差，通过比较显示，本发明提出的方法对状态估计的精确度明显优于现有变分贝斯方法、现有基于t分布的变分贝叶斯方法及经典卡尔曼滤波方法。
[0156]
本发明的实际应用场景如图4所示，在封闭的井下巷道内，通过传感器网络测得传感与井下人员或者设备上标签之间的距离、速度、角度等原始信息，然后根据本发明提出的系统模型和滤波算法对数据进行处理，实现对井下人员和设备的跟踪，实时获得目标的位置、速度等信息。
[0157]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。