CN104331087A - 一种鲁棒的水下传感器网络目标跟踪方法 - Google Patents

Abstract

本发明公开了一种鲁棒的水下传感器网络目标跟踪方法，包括以下步骤：(1)利用正态逆伽玛分布模型对水下传感器节点的测量偏置以及非高斯随机测量噪声进行建模；(2)利用变分贝叶斯近似方法对正态逆伽玛分布的模型参数的更新过程进行求解；(3)结合扩展卡尔曼滤波算法实现对目标状态以及传感器节点测量偏置的估计，并对非高斯测量噪声具有一定的鲁棒性。本发明提供的鲁棒的水下传感器网络目标跟踪方法，充分考虑了传感器节点本身存在测量偏置以及复杂的水下环境导致的非高斯随机测量噪声对目标状态估计的影响，利用正态逆伽玛分布模型对测量偏置和非高斯随机测量噪声进行建模，实现了鲁棒地同时估计目标状态和测量偏置。

Description

一种鲁棒的水下传感器网络目标跟踪方法

技术领域

本发明涉及水下传感器网络技术，具体涉及一种鲁棒的水下传感器网络目标跟踪方法。

背景技术

水下传感器网络(Underwater Sensor Networks，UWSN)在海洋环境监测、海底资源探测、灾难预警、辅助导航、闯入目标的监测与跟踪等领域具有广阔的应用前景。近年来，随着海洋经济的兴起以及各国对海洋权益的日益重视，水下传感器网络成为各个国家的科研机构和学者们共同关注的海洋热点课题之一，目标跟踪是水下传感器网络的一项重要应用。

在传感器节点的测量系统中主要存在两种形式的误差：随机误差和系统误差。在传统的无线传感器网络目标跟踪算法研究中，通常假设传感器节点的随机测量误差被表征为服从高斯分布的加性随机测量噪声，并假设传感器节点的系统误差为零(即假设传感器节点都经过校准)。但是，在实际工程应用中往往会有不符合高斯分布的测量野值(Outlier)存在，对目标状态估计结果造成不可忽略的影响。特别是考虑到水下传感器网络所处的复杂的水下环境，服从零均值高斯分布的随机测量噪声不足于表征实际复杂水下环境下的测量噪声。传感器节点的系统误差一般表现为一种确定性的测量偏置，并且这种测量偏置无法像随机测量噪声一样通过滤波算法进行处理，需要对其进行估计并对测量系统进行系统误差补偿。

学生t分布是一个重尾分布，比高斯噪声更能表征实际随机测量噪声的概率分布。但是由于服从学生t分布的测量噪声无法直接利用卡尔曼滤波器或粒子滤波器进行处理，因此一般把学生t分布近似成方差矩阵变化或增加了方差矩阵权重系数的高斯分布，然后利用伽玛分布、逆伽玛分布、威沙特分布以及逆威沙特分布对方差矩阵或方差矩阵的权重系数进行拟合，并利用变分贝叶斯方法对模型参数更新过程进行求解。基于极大似然的测量偏置估计法和基于伪测量的测量偏置估计法被用于测量偏置的估计与补偿。但是，能够同时处理测量系统的随机测量噪声和固定测量偏置的算法则很少。尽管Ozkan等在2013年的《Marginalized adaptive particle filtering fornonlinear models with unknown time-varying noise parameters》一文中提出了基于粒子滤波器和正态-逆威沙特分布的滤波算法来同时估计高斯测量噪声的均值和方差，但是粒子的快速退化现象大大限制了该算法的应用场景。

发明内容

本发明提供了一种鲁棒的水下传感器网络目标跟踪方法，能够充分处理水下传感器节点的测量偏置以及复杂水下环境导致的非高斯随机测量噪声，在进行目标状态估计的同时，估计出了传感器节点的测量偏置，并对非高斯随机测量噪声具有一定的鲁棒性，满足了复杂水下环境下水下传感器网络目标跟踪的需求。

一种鲁棒的水下传感器网络目标跟踪方法，包括以下步骤：

(1)建立水下传感器节点的测量偏置以及非高斯随机测量噪声的正态逆伽玛分布模型；

(2)利用变分贝叶斯近似方法求解所述的正态逆伽玛分布模型得到的模型参数的更新公式；

(3)利用所述的更新公式，采用扩展卡尔曼滤波算法对目标状态和传感器节点测量偏置进行迭代估计，直至目标状态的估计结果和测量偏置的估计结果同时收敛，并以最后一次迭代估计得到的目标状态和传感器节点测量偏置作为最终估计结果。

本发明充分考虑了未经校准或校准不完全的传感器节点本身可能存在确定性的测量偏置，对测量偏置进行在线估计和补偿，可以避免测量偏置对目标跟踪结果造成的不可忽略的影响。

在实际工程环境中可能会存在不符合高斯分布的测量野值(Outlier)存在，对目标跟踪结果造成不可忽略的影响，充分考虑了复杂水下环境下非高斯随机测量噪声对目标跟踪的影响，对非高斯随机测量噪声具有一定的鲁棒性。

本发明中水下传感器节点进行纯距离测量，因此测量值的维度为1。目标状态包括目标的位置和速度。

所述步骤(1)具体包括如下步骤：

(1-1)利用非零均值且方差变化的高斯分布表示传感器节点的测量残差，即得到的测量残差高斯模型如下：

$p\left(e_n\right|{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)=N\left(e_n\right|{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2);$

其中，μ_n、分别为测量残差高斯分布模型的均值和方差；

e_n为测量残差：

e_n＝z_n-h_n(x_n)，

其中，z_n为传感器节点的纯距离测量值；h_n(x_n)为纯距离测量方程，即运动目标与传感器之间的欧式距离；为运动目标状态，x_1，n、x_2，n和x_3，n为运动目标在三维笛卡尔坐标系内的[x,y,z]三个方位上的位置坐标，和为运动目标在三维笛卡尔坐标系内的[x,y,z]三个方位上的速度。

(1-2)利用逆伽玛分布建立测量残差高斯模型的方差的模型，得到的方差分布模型如下：

$p\left({\mathrm{\σ}}_n^2\right)={\mathrm{\Γ}}^{-1}({\mathrm{\α}}_n,{\mathrm{\β}}_n),$

其中，α_n和β_n分别为逆伽玛分布的模型参数；

(1-3)利用正态分布对测量残差高斯模型的均值μ_n进行建模，得到的均值分布模型如下：

$p\left({\mathrm{\μ}}_n\right|{\mathrm{\σ}}_n^2)=N(m_n,\frac{{\mathrm{\σ}}_n^2}{{\mathrm{\λ}}_n});$

其中，m_n、λ_n分别为正态分布的模型参数；

(1-4)以所述的方差分布模型和均值分布模型作为模型参数，利用正态逆伽玛分布同时表征测量残差高斯模型的均值和方差，得到的方差-均值的正态逆伽玛分布模型具体如下：

$p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)=N-{\mathrm{\Γ}}^{-1}(m_n,{\mathrm{\λ}}_n,{\mathrm{\α}}_n,{\mathrm{\β}}_n),$

其中，N-Γ^-1表示正态逆伽玛分布。

所述步骤(2)具体包括如下步骤：

(2-1)利用变分贝叶斯近似方法把测量残差高斯模型的均值和方差的后验分布近似成两个分布的乘积，具体如下：

$p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n})\≈Q_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)Q_{\mathrm{\σ}}\left({\mathrm{\σ}}_n^2\right).;$

(2-2)利用变分贝叶斯近似方法求解Q_μ(μ_n)，得到模型参数m_n和λ_n的更新过程；

(2-3)利用变分贝叶斯近似方法求解得到模型参数α_n和β_n的更新过程。

所述步骤(2-2)得到的模型参数m_n和λ_n的更新过程如下：

$\begin{array}{c}{\widehat{\mathrm{\λ}}}_n={\stackrel{\‾}{\mathrm{\λ}}}_n+1\\ {\hat{m}}_n=\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}{\stackrel{\‾}{m}}_n+\frac{1}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}{e}_n\end{array},$

其中，和分别表示n时刻的λ和m估计值，和分布表示对n时刻λ和m的估计值更新后得到的修正值。

所述步骤(2-2)得到的模型参数α_n和β_n的更新过程如下：

$\begin{array}{c}{\widehat{\mathrm{\α}}}_n={\stackrel{\‾}{\mathrm{\α}}}_n+\frac{1}{2}\\ {\widehat{\mathrm{\β}}}_n={\stackrel{\‾}{\mathrm{\β}}}_n+\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}\frac{{(e_n-{\stackrel{\‾}{m}}_n)}^2}{2}\end{array},$

其中，和分别表示n时刻的α和β估计值，和分布表示对n时刻α和β的估计值更新后得到的修正值。

所述步骤(3)包括如下步骤：

(3-1)根据目标的运动模型对目标进行目标状态预测和正态逆伽玛分布模型的模型参数预测，其中，目标的运动模型与实际环境相关；

(3-2)利用模型参数预测结果设置参数后验估计迭代过程的初值；

(3-3)利用估计的当前时刻n的模型参数μ_n、以及当前时刻n的测量值z_n估计最新的目标状态；

(3-4)基于最新的目标状态估计，按照步骤(2-2)和(2-3)中的推导结果对参数m_n、λ_n、α_n和β_n进行更新；

(3-5)重复步骤(3-3)和(3-4)，直至目标状态估计结果和测量偏置估计结果同时收敛。

目标状态估计结果和测量偏置估计结果同时收敛的评判方式为：在连续两次执行步骤(3-3)和(3-4)过程中，如果后一次次执行过程得到的目标状态估计结果和测量偏置估计结果，相比前一次执行得到的相应估计结果的变化率小于5％，则认为收敛。

本发明鲁棒的水下传感器网络目标跟踪方法，充分考虑了水下传感器节点的测量偏置以及复杂水下环境导致的非高斯随机测量噪声，通过对测量偏置进行在线估计和补偿，降低了测量偏置对目标跟踪结果的影响；通过利用非高斯模型对随机测量噪声进行建模，降低了测量野值对目标跟踪结果的影响。本发明方法在进行目标状态估计的同时，估计出了传感器节点的测量偏置，并对非高斯随机测量噪声具有一定的鲁棒性，满足了复杂水下环境下水下传感器网络目标跟踪的需求。

附图说明

图1为本发明方法的流程图。

具体实施方式

为了更为具体地描述本发明，下面结合附图和具体实施方式对本发明的鲁棒的水下传感器网络目标跟踪方法做详细描述。

如图1所示的一种鲁棒的水下传感器网络目标跟踪方法，包括以下步骤：

(1)利用正态逆伽玛分布模型对水下传感器节点的测量偏置以及非高斯随机测量噪声进行建模；

(2)利用变分贝叶斯近似方法对正态逆伽玛分布的模型参数的更新过程进行求解；

(3)利用所述的更新公式，采用扩展卡尔曼滤波算法对目标状态和传感器节点测量偏置进行迭代估计，直至目标状态的估计结果和测量偏置估计结果同时收敛，并以最后一次迭代估计得到的目标状态和传感器节点测量偏置作为最终估计结果。

其中，水下传感器节点进行纯距离测量，因此传感器节点的测量值的维度为1。

本实施例中步骤(1)具体包括如下步骤：

(1-1)利用非零均值且方差变化的高斯分布对传感器节点的测量残差进行建模；

水下运动目标的运动模型可以用如下式子表示：

$x_n=f_n\left(x_{n-1}\right)+u_n,u_n\stackrel{i.i.d.}{~}N(0,Q)$

其中，为运动目标在时刻的状态，x_1，n、x_2，n和x_3，n为运动目标在三维笛卡尔坐标系内的[x,y,z]三个方位上的位置坐标，和为运动目标在三维笛卡尔坐标系内的[x,y,z]三个方位上的速度。u_n为服从高斯分布的过程噪声。f_n()为运动目标的过程转移方程，具体形式与实际情况下目标的运动模型相关，本发明中可以认为是固定速度运动模型。考虑到传感器节点的测量偏置和非高斯随机测量噪声，传感器节点的测量模型可以表示为：

$z_n=h_n\left(x_n\right)+{\mathrm{\υ}}_n,{\mathrm{\υ}}_n\stackrel{i.i.d.}{~}N({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)$

其中，z_n为传感器节点的纯距离测量值，υ_n为包含了测量偏置和随机测量噪声的测量噪声，μ_n表示传感器节点的测量偏置，为高斯噪声模型的方差并且该方差值未知。令测量残差e_n＝z_n-h_n(x_n)，则可以得出：

$p\left(e_n\right|{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)=N\left(e_n\right|{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp (-\frac{{(e_n-{\mathrm{\μ}}_n)}^2}{2{\mathrm{\σ}}_n^2})$

(1-2)利用逆伽玛分布对测量残差高斯模型的方差进行建模：

$p\left({\mathrm{\σ}}_n^2\right)={\mathrm{\Γ}}^{-1}({\mathrm{\α}}_n,{\mathrm{\β}}_n)=\frac{{\mathrm{\β}}_n^{{\mathrm{\α}}_n}}{\mathrm{\Γ}\left({\mathrm{\α}}_n\right)}{\left(\frac{1}{{\mathrm{\σ}}_n^2}\right)}^{{\mathrm{\α}}_n+1}\exp (-\frac{{\mathrm{\β}}_n}{{\mathrm{\σ}}_n^2})$

(1-3)利用正态分布对测量残差高斯模型的均值进行建模：

$p\left({\mathrm{\μ}}_n\right|{\mathrm{\σ}}_n^2)=N(m_n,\frac{{\mathrm{\σ}}_n^2}{{\mathrm{\λ}}_n})=\frac{\sqrt{{\mathrm{\λ}}_n}}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\exp (-\frac{{\mathrm{\λ}}_n{({\mathrm{\μ}}_n-m_n)}^2}{2{\mathrm{\σ}}_n^2})$

(1-4)结合(1-2)中的逆伽玛分布和(1-3)中的正态分布，利用正态逆伽玛分布同时表征测量残差高斯模型的均值和方差：

$\begin{array}{c}p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)\\ =p\left({\mathrm{\μ}}_n\right|{\mathrm{\σ}}_n^2)p\left({\mathrm{\σ}}_n^2\right)\\ =N-{\mathrm{\Γ}}^{-1}(m_n,{\mathrm{\λ}}_n,{\mathrm{\α}}_n,{\mathrm{\β}}_n)\\ =\frac{\sqrt{{\mathrm{\λ}}_n}}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_n}\frac{{\mathrm{\β}}_n^{{\mathrm{\α}}_n}}{\mathrm{\Γ}\left({\mathrm{\α}}_n\right)}{\left(\frac{1}{{\mathrm{\σ}}_n^2}\right)}^{{\mathrm{\α}}_n+1}\exp (-\frac{2{\mathrm{\β}}_n+{\mathrm{\λ}}_n{({\mathrm{\μ}}_n-m_n)}^2}{2{\mathrm{\σ}}_n^2})\end{array}$

本实施例中步骤(2)具体包括如下步骤：

(2-1)测量残差高斯模型的均值和方差的后验分布表现出非高斯特性，利用变分贝叶斯近似方法把该非高斯后验分布近似成两个分布的乘积，即 $p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n})\≈Q_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)Q_{\mathrm{\σ}}\left({\mathrm{\σ}}_n^2\right);$

令 $p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n})=N-{\mathrm{\Γ}}^{-1}({\hat{m}}_n,{\widehat{\mathrm{\λ}}}_n,{\widehat{\mathrm{\α}}}_n,{\widehat{\mathrm{\β}}}_n)$ 表示均值和方差参数基于测量残差e_1：n的后验分布，令 $p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n-1})=N-{\mathrm{\Γ}}^{-1}({\stackrel{\‾}{m}}_n,{\stackrel{\‾}{\mathrm{\λ}}}_n,{\stackrel{\‾}{\mathrm{\α}}}_n,{\stackrel{\‾}{\mathrm{\β}}}_n)$ 表示均值和方差参数基于测量残差e_1：n-1的先验分布。

在变分贝叶斯近似方法中，通过近似分布和真实后验分布的KL-散度的最小化来对近似分布进行求解。该近似分布和真实分布的KL-散度定义为：

$\begin{array}{c}\mathrm{KL}\left(Q_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)Q_{\mathrm{\σ}}\left({\mathrm{\σ}}_n^2\right)\right|\left|p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n})\right)\\ ={\∫Q}_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)Q_{\mathrm{\σ}}\left({\mathrm{\σ}}_n^2\right)\log \frac{Q_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)Q_{\mathrm{\σ}}\left({\mathrm{\σ}}_n^2\right)}{p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n})}d{\mathrm{\μ}}_{n}d{\mathrm{\σ}}_n^2\end{array}.$

根据变分贝叶斯近似方法基于最小化KL-散度的求解过程，近似分布Q_μ(μ_n)(或)相对于(或Q_μ(μ_n))的解的形式如下：

$Q_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)\∝\exp (\∫\log p(e_n,{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n-1})Q_{\mathrm{\σ}}\left({\mathrm{\σ}}_n\right)d{\mathrm{\σ}}_n^2)$

$Q_{\mathrm{\Σ\μ}}\left({\mathrm{\σ}}_n^2\right)\∝\exp (\∫\log p(e_n,{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n-1})Q_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)d{\mathrm{\μ}}_n)$

(2-2)利用变分贝叶斯近似方法求解Q_μ(μ_n)，进而得到模型参数m_n和λ_n的更新过程；

由于正态逆伽玛分布是单变量正态分布的均值和方差的共轭分布，因此Q_μ(μ_n)仍然服从正态分布，则通过化简可以得出参数更新过程如下：

$\begin{array}{c}{\widehat{\mathrm{\λ}}}_n={\stackrel{\‾}{\mathrm{\λ}}}_n+1\\ {\hat{m}}_n=\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}{\stackrel{\‾}{m}}_n+\frac{1}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}{e}_n\end{array}$

(2-3)利用变分贝叶斯近似方法求解进而得到模型参数α_n和β_n的更新过程。

由于正态逆伽玛分布是单变量正态分布的均值和方差的共轭分布，因此仍然服从逆伽玛分布，则通过化简可以得出参数更新过程如下：

$\begin{array}{c}{\widehat{\mathrm{\α}}}_n={\stackrel{\‾}{\mathrm{\α}}}_n+\frac{1}{2}\\ {\widehat{\mathrm{\β}}}_n={\stackrel{\‾}{\mathrm{\β}}}_n+\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}\frac{{(e_n-{\stackrel{\‾}{m}}_n)}^2}{2}\end{array}.$

本实施例中步骤(3)具体包括如下步骤：

(3-1)进行目标状态预测和正态逆伽玛分布参数预测；

首先，根据目标运动模型进行目标状态预测：

${\stackrel{\‾}{x}}_n=F_n{\hat{x}}_{n-1}$

${\stackrel{\‾}{P}}_n=F_n{\hat{P}}_{n-1}{F}_n^T+Q_n$

其中，和为上一时刻的目标状态和目标状态协方差的后验估计，和为当前时刻的目标状态和目标状态协方差的先验预测，F_n为(1-1)中过程转移方程的雅克比矩阵，Q_n为步骤(1-1)中过程噪声的协方差矩阵。

其次，对正态逆伽玛分布参数进行线性预测：

${\stackrel{\‾}{\mathrm{\λ}}}_n=\mathrm{\ρ}{\widehat{\mathrm{\λ}}}_{n-1},$

${\stackrel{\‾}{m}}_n={\hat{m}}_{n-1},$

${\stackrel{\‾}{\mathrm{\α}}}_n=\mathrm{\ρ}{\widehat{\mathrm{\α}}}_{n-1},$

${\stackrel{\‾}{\mathrm{\β}}}_n=\mathrm{\ρ}{\widehat{\mathrm{\β}}}_{n-1}.$

其中，ρ为线性预测系数，取值根据实际应用场景设定，本实施例中ρ＝1-e^-4。

(3-2)初始化正态逆伽玛分布参数后验估计过程迭代初值，即利用正态逆伽玛分布参数预测结果设置参数后验估计迭代过程的初值；

设置参数后验估计过程中参数迭代初值：

${\widehat{\mathrm{\λ}}}_n={\stackrel{\‾}{\mathrm{\λ}}}_n+1$

${\widehat{\mathrm{\α}}}_n={\stackrel{\‾}{\mathrm{\α}}}_n+\frac{1}{2}.$

${\hat{m}}_n={\stackrel{\‾}{m}}_n$

${\widehat{\mathrm{\β}}}_n={\stackrel{\‾}{\mathrm{\β}}}_n$

(3-3)利用基于模型参数后验估计结果计算出的μ_n、以及最新的测量值z_n估计最新的目标状态；

首先，根据正态逆伽玛分布参数后验估计计算测量偏置以及测量噪声方差

其次，根据测量值以及估计出的测量偏置和测量噪声方差计算目标状态后验估计：

$\begin{array}{c}{K}_n=\frac{{\stackrel{\‾}{P}}_n{H}_n^T}{H_n{\stackrel{\‾}{P}}_n{H}_n^T+{\mathrm{\σ}}_n^2}\\ {\hat{x}}_n={\stackrel{\‾}{x}}_n+K_n(z_n+h_n\left({\stackrel{\‾}{x}}_n\right)-{\mathrm{\μ}}_n)\\ {\hat{P}}_n=(I-K_n{H}_n){\stackrel{\‾}{P}}_n\end{array},$

其中，K_n为计算得到的最佳卡尔曼增益；z_n为传感器节点的测量值，h_n()为纯距离观测方程，即运动目标与传感器之间的欧式距离；H_n为h_n()的雅克比矩阵；I为单位矩阵；和为n时刻的目标状态和目标状态协方差的后验估计，和为n时刻的目标状态和目标状态协方差的先验预测。

(3-4)基于最新的n时刻的目标状态估计，按照步骤(2-2)和(2-3)中的推导结果对参数m_n和β_n进行更新；

${\hat{m}}_n=\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}{\stackrel{\‾}{m}}_n+\frac{1}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}(z_n-h_n\left({\hat{x}}_n\right))$

${\widehat{\mathrm{\β}}}_n={\stackrel{\‾}{\mathrm{\β}}}_n+\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}\frac{{(z_n-h_n\left({\hat{x}}_n\right)-{\stackrel{\‾}{m}}_n)}^2}{2}$

(3-5)重复执行步骤(3-3)和(3-4)，直至目标状态估计结果和测量偏置估计结果同时收敛，输出最后一次迭代得到的目标状态后验估计和正态逆伽玛分布参数(即模型参数)后验估计

目标状态估计结果和测量偏置估计结果同时收敛的评判方式为：在连续两次执行步骤(3-3)和(3-4)过程中，如果后一次执行过程输出的目标状态估计结果和测量偏置估计结果，相比前一次执行结果的变化率小于5％，则认为收敛。

每次重复步骤(3-3)和(3-4)时，使用的和的初值不变，和需要采用上一次执行过程更新后(3-4)中输出的结果作为初始值。

需要特别指出的时，上述过程只是描述了本发明算法的执行过程，需要用到的参数的初值与具体环境有关：运动目标在0时刻的初始状态和与具体环境有关，可以被设置为0，可以被设置为5，可以被设置为2，和具体环境下的噪声水平有关，ρ可以被设置为1-e^-4。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (7)

1.一种鲁棒的水下传感器网络目标跟踪方法，其特征在于，包括以下步骤：

(1)建立水下传感器节点的测量偏置以及非高斯随机测量噪声的正态逆伽玛分布模型；

(2)利用变分贝叶斯近似方法求解所述的正态逆伽玛分布模型得到的模型参数的更新公式；

(3)利用所述的更新公式，采用扩展卡尔曼滤波算法对目标状态和传感器节点的测量偏置进行迭代估计，直至目标状态的估计结果和测量偏置的估计结果同时收敛，并以最后一次迭代估计得到的目标状态和传感器节点测量偏置作为最终估计结果。

2.如权利要求1所述的鲁棒的水下传感器网络目标跟踪方法，其特征在于，所述步骤(1)具体包括如下步骤：

(1-1)利用非零均值且方差变化的高斯分布表示传感器节点的测量残差，即得到的测量残差高斯模型如下：

$p\left(e_n\right|{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)=N\left(e_n\right|{\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2);$

其中，μ_n、分别为测量残差高斯分布模型的均值和方差；

e_n为测量残差：

e_n=z_n-h_n(x_n)，

其中，z_n为传感器节点的纯距离测量值；h_n(x_n)为纯距离测量方程，即运动目标与传感器之间的欧式距离；为运动目标状态，x_1，n、x_2，n和x_3，n为运动目标在三维笛卡尔坐标系内的[x,y,z]三个方位上的位置坐标，和为运动目标在三维笛卡尔坐标系内的[x,y,z]三个方位上的速度。

(1-2)利用逆伽玛分布建立测量残差高斯模型的方差的模型，得到的方差分布模型如下：

$p\left({\mathrm{\σ}}_n^2\right)={\mathrm{\Γ}}^{-1}({\mathrm{\α}}_n,{\mathrm{\β}}_n),$

其中，α_n和β_n分别为逆伽玛分布的模型参数；

(1-3)利用正态分布对测量残差高斯模型的均值μ_n进行建模，得到的均值分布模型如下：

$p\left({\mathrm{\μ}}_n\right|{\mathrm{\σ}}_n^2)=N(m_n,\frac{{\mathrm{\σ}}_n^2}{{\mathrm{\λ}}_n});$

其中，m_n、λ_n分别为正态分布的模型参数；

(1-4)以所述的方差分布模型和均值分布模型作为模型参数，利用正态逆伽玛分布同时表征测量残差高斯模型的均值和方差，得到的方差-均值的正态逆伽玛分布模型具体如下：

$p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2)=N-{\mathrm{\Γ}}^{-1}(m_n,{\mathrm{\λ}}_n,{\mathrm{\α}}_n,{\mathrm{\β}}_n),$

其中，N-Γ^-1表示正态逆伽玛分布。

3.根据权利要求2所述的鲁棒的水下传感器网络目标跟踪方法，其特征在于，所述步骤(2)具体包括如下步骤：

(2-1)利用变分贝叶斯近似方法把测量残差高斯模型的均值和方差的后验分布近似成两个分布的乘积，具体如下：

$p({\mathrm{\μ}}_n,{\mathrm{\σ}}_n^2|e_{1:n})\≈Q_{\mathrm{\μ}}\left({\mathrm{\μ}}_n\right)Q_{\mathrm{\σ}}\left({\mathrm{\σ}}_n^2\right).;$

(2-2)利用变分贝叶斯近似方法求解Q_μ(μ_n)，得到模型参数m_n和λ_n的更新过程；

(2-3)利用变分贝叶斯近似方法求解得到模型参数α_n和β_n的更新过程。

4.根据权利要求3所述的鲁棒的水下传感器网络目标跟踪方法，其特征在于，所述步骤(2-2)得到的模型参数m_n和λ_n的更新过程如下：

$\begin{array}{c}{\widehat{\mathrm{\λ}}}_n={\stackrel{\‾}{\mathrm{\λ}}}_n+1\\ {\hat{m}}_n=\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}{\stackrel{\‾}{m}}_n+\frac{1}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}{e}_n\end{array},$

其中，和分别表示n时刻的λ和m估计值，和分布表示对n时刻λ和m的估计值更新后得到的修正值。

5.根据权利要求3所述的鲁棒的水下传感器网络目标跟踪方法，其特征在于，所述步骤(2-2)得到的模型参数α_n和β_n的更新过程如下：

$\begin{array}{c}{\widehat{\mathrm{\α}}}_n={\stackrel{\‾}{\mathrm{\α}}}_n+\frac{1}{2}\\ {\widehat{\mathrm{\β}}}_n={\stackrel{\‾}{\mathrm{\β}}}_n+\frac{{\stackrel{\‾}{\mathrm{\λ}}}_n}{{\stackrel{\‾}{\mathrm{\λ}}}_n+1}\frac{{(e_n-{\stackrel{\‾}{m}}_n)}^2}{2}\end{array},$

其中，和分别表示n时刻的α和β估计值，和分布表示对n时刻α和β的估计值更新后得到的修正值。

6.如权利要求4或5所述的鲁棒的水下传感器网络目标跟踪方法，其特征在于，所述步骤(3)包括如下步骤：

(3-1)根据目标的运动模型对目标进行目标状态预测和正态逆伽玛分布模型的模型参数预测；

(3-2)利用模型参数预测结果设置参数后验估计迭代过程的初值；

(3-3)利用估计的当前时刻n的模型参数μ_n、以及当前时刻n的测量值z_n估计最新的目标状态；

(3-4)基于最新的目标状态估计，按照步骤(2-2)和(2-3)中的推导得到的更新过程对参数m_n、λ_n、α_n和β_n进行更新；

(3-5)重复步骤(3-3)和(3-4)，直至目标状态估计结果和参数估计结果同时收敛。

7.如权利要求6所述的鲁棒的水下传感器网络目标跟踪方法，其特征在于，目标状态估计结果和测量偏置估计结果同时收敛的评判方式为：在连续两次执行步骤(3-3)和(3-4)过程中，如果后一次次执行过程得到的目标状态估计结果和测量偏置估计结果，相比前一次执行得到的相应估计结果的变化率小于5％，则认为收敛。