CN101795123A - 分布式非线性滤波方法 - Google Patents

Abstract

一种信号处理技术领域的分布式非线性滤波方法，包括：网络中的每一个传感器节点对目标状态方程和测量方程进行加权统计线性化处理得到加统计线性化的系统矩阵、测量矩阵以及线性化后的系统噪声与测量噪声；传感器节点与其相邻节点进行通信以交换局部信息贡献，并根据动态协同滤波器获得全局信息贡献；传感器节点根据加权统计线性化后的目标状态方程对目标当前时刻的状态进行预测，并计算对应的方差阵；然后根据步骤二中的由协同滤波器获得的全局信息贡献对目标状态进行更新，获得目标在当前时刻的状态估值。本发明减小了滤波处理的通信复杂度。当应用于分散无迹卡尔曼滤波器的通信量为O(2450)，然而本方法所需通信量减少至O(230)。

Description

分布式非线性滤波方法

技术领域

本发明涉及的是一种信号处理技术领域的方法，具体是一种可扩展规模的分布式非线性滤波方法。

背景技术

随着计算机网络、无线通信和微小型系统的发展，融合以上三种技术的无线传感器网络应运而生。无线传感器网络带来了一种全新的信息获取与处理模式，将深刻影响着信息技术的未来发展。目标跟踪是无线传感器网络最具代表性和挑战性的应用。目标跟踪可分为集中式和分布式两种方式。集中式跟踪可以实现高精度的数据处理，但由于所需数据传输量大、对融合中心要求苛刻、处理时间长，容易产生丢包和延时，从而使得跟踪精度降低。分布式是无线传感器网络的本质特性，也是解决跟踪复杂性与节点能力不足之间矛盾的有力武器。

经过对现有技术的检索发现，传统的分散(decentralized)滤波与融合方法，如分散卡尔曼滤波器和分散无迹卡尔曼滤波器(T.Vercauteren，and X.Wang，“Decentralized sigma-pointinformation filters for target tracking in collaborative sensor networks，”IEEE T.SignalProcessing，vol.53，no.8，pp.2997-3009，Aug.2005)，都是全局到全局的，即每个节点都需要与网络中的所有其它节点或者融合中心进行通信，因此，以上两种滤波或融合方法所需的通信复杂度是O(N＊N)，其中N是网络中的传感器节点或智能体的数目)。很明显，由于通信复杂度高以上方法是不能扩展规模的，尤其对于大规模传感器网络是不适用的。

发明内容

本发明针对现有技术存在的上述不足，提供一种分布式非线性滤波方法，网络中的每一个节点只需与其相邻节点进行信息交换，进而基于动态协同滤波处理，减小了滤波处理的通信复杂度。当应用于分散无迹卡尔曼滤波器的通信量为O(2450)，然而本方法所需通信量减少至O(230)。

本发明通过如下技术方案实现的，本发明包括以下步骤：

步骤一、网络中的每一个传感器节点对目标状态方程和测量方程进行加权统计线性化处理得到加统计线性化的系统矩阵、测量矩阵以及线性化后的系统噪声与测量噪声；

所述的加权统计线性化处理包括目标状态方程和测量方程，具体如下：

目标状态方程x(k+1)＝F(k)x(k)+b^x(k)+ω(k)

测量方程 $z_i\left(k\right)={\stackrel{\‾}{H}}_i\left(k\right)x\left(k\right)+b^{z_i}\left(k\right)+{\stackrel{\‾}{\mathrm{\υ}}}_i\left(k\right)$

其中：F(k)和H_i(k)分别为加权统计线性化以后的目标状态矩阵和测量矩阵；而b^x(k)和是由线性化引起的误差项，ω(k)和υ_i(k)线性化以后的过程噪声和测量噪声。

步骤二、传感器节点与其相邻节点进行通信以交换局部信息贡献，并根据动态协同滤波器获得全局信息贡献。

所述的局部信息贡献是指：

$U_i\left(k\right)={\stackrel{\‾}{H}}_i^T\left(k\right){\stackrel{\‾}{R}}_i^{-1}\left(k\right){\stackrel{\‾}{H}}_i\left(k\right),$ 以及 $u_i\left(k\right)={\stackrel{\‾}{H}}_i^T\left(k\right){\stackrel{\‾}{R}}_i^{-1}\left(k\right)[z_i\left(k\right)-b_i^z\left(k\right)];$

其中：R_i ^-1(k)为加权统计线性化后的测量噪声的方差阵，H_i(k)为加权统计线性化以后的测量矩阵，

是由线性化引起的误差项。

所述的动态协同滤波器是指：

其中：k是采样步数；δ是协同滤波器的更新步长，增益β＞0且β～O(1/λ₂)，λ₂是网络图模型的Laplacian矩阵L的第二小特征值，协同滤波器的输出

和

分别为对其输入U_i(k)或u_i(k)的估计。

步骤三、传感器节点根据加权统计线性化后的目标状态方程对目标当前时刻的状态进行预测，并计算对应的方差阵；然后根据步骤二中的由协同滤波器获得的全局信息贡献对目标状态进行更新，获得目标在当前时刻的状态估值。

与现有技术相比，本发明具有如下有益效果：

1、网络中的每一个节点只需与其相邻节点(而不是与网络中所有其它节点)进行信息交换，基于动态协同滤波处理，所有节点都能对目标的状态估计达成一致。这样，大大减小了滤波处理的通信复杂度；

2、在降低通信量的同时，本方法具有与集中式融合方法和分散融合方法可比的跟踪性能；

3、此外只要整个网络是连通的，整个网络中所有节点都能达到一致性，从而增强了处理在通信失败和节点失效情况下的鲁棒性；

因此，本方法在大规模传感器网络和多智能体系统等军用和民用领域有广泛的应用前景。

附图说明

图1为本发明流程图。

图2为50×50方形区域传感器部署及通信链接示意图。

图3为实施例x方向的均方根比较图。

图4为实施例y方向的均方根比较图。

图5为移动目标的真实航迹及估计效果示意图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的实施过程，但本发明的保护范围不限于下述的实施例。

如图1所示，本实施例包括以下步骤：

步骤1、网络中的每一个传感器节点对目标状态方程和测量方程进行加权统计线性化。具体来说，首先根据当前时刻的目标状态估值

和方差P_xx生成(2n+1)个sigma点{χ_j(k|k)，ω_j}：

${\mathrm{\χ}}_0\left(k\right|k)=\hat{x}\left(k\right|k),{\mathrm{\ω}}_0=\frac{\mathrm{\κ}}{n+\mathrm{\κ}}$ (公式一)

${\mathrm{\χ}}_j\left(k\right|k)=\hat{x}\left(k\right|k)+{\left(\sqrt{(n+\mathrm{\κ})P_{\mathrm{xx}}}\right)}_j,$ ${\mathrm{\ω}}_j=\frac{1}{2(n+\mathrm{\κ})},j=\mathrm{1,2},..,n$ (公式二)

${\mathrm{\χ}}_{j+n}\left(k\right|k)=\hat{x}\left(k\right|k)-{\left(\sqrt{(n+\mathrm{\κ})P_{\mathrm{xx}}}\right)}_j,$ ${\mathrm{\ω}}_{j+n}=\frac{1}{2(n+\mathrm{\κ})},j=\mathrm{1,2},..,n$ (公式三)

其中：尺度参数κ通常取0或3-n(n为目标状态变量的维数)，

表示P矩阵Cholesky分解的第j行；对sigma点的要求在于：χ_j(k|k)的均值与方差与当前时刻目标状态估计的先验信息一致，即 $\stackrel{\‾}{x}=\hat{x},$ P_xx＝P_xx，其中 $\stackrel{\‾}{x}={\mathrm{\Σ}}_{j=1}^r{\mathrm{\ω}}_j{\mathrm{\χ}}_j,$ ${\stackrel{\‾}{P}}_{\mathrm{xx}}={\mathrm{\Σ}}_{j=1}^r{\mathrm{\ω}}_j({\mathrm{\χ}}_j-\stackrel{\‾}{x}){({\mathrm{\χ}}_j-\stackrel{\‾}{x})}^T,$ 且 ${\mathrm{\Σ}}_{j=1}^r{\mathrm{\ω}}_j=1.$

然后，每个sigma点都通过目标状态非线性函数进行传播；再根据(公式一)-(公式三)产生(2n+1)个sigma点{χ_j(k+1|k)，ω_j}并且通过传感器节点的非线性测量方程进行传播，得到：

x(k+1)＝F(k)x(k)+b^x(k)+ω(k)(公式四)

$z_i\left(k\right)={\stackrel{\‾}{H}}_i\left(k\right)x\left(k\right)+b^{z_i}\left(k\right)+{\stackrel{\‾}{\mathrm{\υ}}}_i\left(k\right)$ (公式五)

其中：F(k)和H_i(k)分别为加权统计线性化以后的目标状态矩阵和测量矩阵；而b^x(k)和

是由线性化引起的误差项，ω(k)和υ_i(k)线性化以后的过程噪声和测量噪声。从而根据加权统计线性化的思想，(公式四)-(公式五)中的相关矩阵可以计算如下：

$\stackrel{\‾}{F}\left(k\right)={\stackrel{\‾}{P}}_{\mathrm{xx}}^T(k+1,k+1|k){\stackrel{\‾}{P}}_{\mathrm{xx}}^{-1}\left(k\right|k)$

b^x(k)＝x(k+1|k)-F(k)x(k|k)(公式六)

${\stackrel{\‾}{H}}_i\left(k\right)={\stackrel{\‾}{P}}_{{\mathrm{xz}}_i}^T(k+1|k){\stackrel{\‾}{P}}_{\mathrm{xx}}^{-1}(k+1|k)$

$b^{z_i}\left(k\right)={\stackrel{\‾}{z}}_i\left(k\right|k-1)-{\stackrel{\‾}{H}}_i\left(k\right)\stackrel{\‾}{x}\left(k\right|k-1)$ (公式七)

加权统计线性化后的过程噪声和测量噪声分别为

ω(k)＝ω(k)+ε^x(k)(公式八)

${\stackrel{\‾}{\mathrm{\υ}}}_i\left(k\right)={\mathrm{\υ}}_i\left(k\right)+{\mathrm{\ϵ}}_i^{z_i}\left(k\right)$ (公式九)

其均值为零，方差分别为

$\mathrm{Var}\left(\stackrel{\‾}{\mathrm{\ω}}\left(k\right)\right)=\stackrel{\‾}{Q}\left(k\right)=Q\left(k\right)+{\stackrel{\‾}{P}}_{\mathrm{\ϵ\ϵ}}^x\left(k\right)$

(公式八)

$=Q\left(k\right)+{\stackrel{\‾}{P}}_{\mathrm{xx}}(k+1|k)-\stackrel{\‾}{F}\left(k\right){\stackrel{\‾}{P}}_{\mathrm{xx}}\left(k\right|k){\stackrel{\‾}{F}}^T\left(k\right)$

$\mathrm{Var}\left({\stackrel{\‾}{\mathrm{\υ}}}_i\left(k\right)\right)={\stackrel{\‾}{R}}_i\left(k\right)=R_i\left(k\right)+{\stackrel{\‾}{P}}_{\mathrm{\ϵ\ϵ}}^{z_i}\left(k\right)$

(公式九)

$=R_i\left(k\right)+{\stackrel{\‾}{P}}_{z_i{z}_i}\left(k\right|k-1)-{\stackrel{\‾}{H}}_i\left(k\right){\stackrel{\‾}{P}}_{\mathrm{xx}}\left(k\right|k-1){\stackrel{\‾}{H}}_i^T\left(k\right)$

步骤2、网络中的每一个传感器节点根据(公式十)和(公式十一)并行地更新局部信息贡献，其中的相关变量由步骤1得到；

$U_i\left(k\right)={\stackrel{\‾}{H}}_i^T\left(k\right){\stackrel{\‾}{R}}_i^{-1}\left(k\right){\stackrel{\‾}{H}}_i\left(k\right)$ (公式十)

$u_i\left(k\right)={\stackrel{\‾}{H}}_i^T\left(k\right){\stackrel{\‾}{R}}_i^{-1}\left(k\right)[z_i\left(k\right)-b_i^z\left(k\right).]$ (公式十一)

步骤3、网络中的每一个传感器节点与其相邻节点通信以交换u_i(k)和U_i(k)；

步骤4、网络中的每一个传感器节点根据以下协同滤波器估计全局信息贡献：

(公式十二)

(公式十三)

步骤5、网络中的每一个传感器节点根据(公式十四)和(公式十五)进行状态预测；即得到目标状态的预测值

和对应的方差阵P_i(k|k-1)：

${\hat{x}}_i\left(k\right|k-1)=\stackrel{\‾}{F}\left(k\right){\hat{x}}_i(k-1|k-1)+b^x\left(k\right)$ (公式十四)

P_i(k|k-1)＝F(k)P_i(k-1|k-1)F^T(k)+Q^μ(k)(公式十五)

其中Q^μ(k)＝NQ(k)。

步骤6、网络中的每一个传感器节点根据(公式十六)和(公式十七)进行状态更新，即

$P_i\left(k\right|k)={[{\left(P_i\left(k\right|k-1)\right)}^{-1}+{\hat{U}}_i\left(k\right)]}^{-1}$ (公式十六)

${\hat{x}}_i\left(k\right|k)={\hat{x}}_i\left(k\right|k-1)+P_i\left(k\right|k)[{\hat{u}}_i\left(k\right)-{\hat{U}}_i\left(k\right){\hat{x}}_i\left(k\right|k-1)]$ (公式十七)

步骤7、下一个采样周期内，网络中的每一个传感器循环步骤1到步骤6。

如图2所示，本实施例仿真条件如下：考虑N＝50个传感器节点随机分布在50米×50米的区域的传感器网络。选择一个圆周运动目标跟踪的例子，假设目标运动的状态方程如下：

x(k+1)＝Fx(k)+Gω(k)

其中x(k)＝[x_k，y_k]^T是目标的状态向量，表示目标在第k个采样点上的位置；本实例取F＝[1-0.05；0.051]，G＝0.025*I₂，采样步长T＝0.025。目标的运动轨迹如图2中的大圆圈所示，并假设第i个传感器节点具有如下测量方程

y_i(k)＝a/‖(x(k)，y(k))-(x_s(i)，y_s(i))‖+υ_i(k)

其中a＝40，‖(x(k)，y(k))-(x_s(i)，y_s(i))‖目标与第i个传感器的距离，测量噪声υ_i(k)的方差为 $R_i\left(k\right)=\sqrt{i}$ (i＝1，2，...，50)。

本实施例仿真内容及结果：经过100次蒙特卡罗仿真实验，将本实施例的可扩展规模的分布式非线性滤波方法与集中式融合方法(采用无迹卡尔曼滤波器)以及分布式融合方法(假设每个传感器都根据自己的测量进行状态估计，然后将估计值及方差送至融合中心，融合中心采用加权平均法进行航迹融合)进行比较。

如图3和图4所示，比较结果中分别列出了第26，11，44，16，45个传感器节点的估计结果，它们的连接度分别为2，4，6，7和8。很明显，本实施例可扩展规模的分布式非线性滤波方法获得与分布式融合方法等同(如11)甚至更好(如44，16和45)的跟踪精度。更重要的是本实施例方法获得的跟踪精度非常接近集中式融合方法。由于没有信息丢失，众所周知，后者最优融合估计。三者的区别在于集中式融合方法和分布式融合方法都需要各局部传感器节点与融合中心进行通信，存在带融合中心传感器网络的固有缺陷：第一，这两种融合方法都是不可扩展规模的；第二，如果出现融合中心故障，那么整个系统就会瘫痪。相反，本实施例方法中没有特定的融合中心，每个传感器节点都对目标状态进行估计并达到一致，同时每个传感器节点都能达到接近最优的估计精度。每个节点的估计一致性允许用户通过查询网络中的任何一个节点获得目标的状态估计。

如图5所示，在不同采样时刻上所有传感器节点基于本实施例可扩展规模的分布式非线性滤波方法对目标状态的估值。可见估计的一致性越来越高，且各个传感器节点对目标状态的估计如同粒子围绕在目标真实轨迹的周围。另外，与分散无迹卡尔曼滤波方法相比，本实施例方法所需通信量由O(2450)减少为O(230)，这样在很大程度上减少了网络拥塞的可能性。

Claims (4)

1.一种分布式非线性滤波方法，其特征在于，包括以下步骤：

步骤一、网络中的每一个传感器节点对目标状态方程和测量方程进行加权统计线性化处理得到加统计线性化的系统矩阵、测量矩阵以及线性化后的系统噪声与测量噪声；

步骤二、传感器节点与其相邻节点进行通信以交换局部信息贡献，并根据动态协同滤波器获得全局信息贡献；

步骤三、传感器节点根据加权统计线性化后的目标状态方程对目标当前时刻的状态进行预测，并计算对应的方差阵；然后根据步骤二中的由协同滤波器获得的全局信息贡献对目标状态进行更新，获得目标在当前时刻的状态估值。

2.根据权利要求1所述的分布式非线性滤波方法，其特征是，步骤一中所述的加权统计线性化处理包括目标状态方程和测量方程，具体如下：

目标状态方程x(k+1)＝F(k)x(k)+b^x(k)+ω(k)

测量方程 $z_i\left(k\right)={\stackrel{\‾}{H}}_i\left(k\right)x\left(k\right)+b^{z_i}\left(k\right)+{\stackrel{\‾}{\mathrm{\υ}}}_i\left(k\right)$

其中：F(k)和H_i(k)分别为加权统计线性化以后的目标状态矩阵和测量矩阵，b^x(k)和

是由线性化引起的误差项，ω(k)和υ_i(k)为线性化以后的过程噪声和测量噪声。

3.根据权利要求1所述的分布式非线性滤波方法，其特征是，步骤二中所述的局部信息贡献是指：

$U_i\left(k\right)={\stackrel{\‾}{H}}_i^T{\left(k\right)\stackrel{\‾}{R}}_i^{-1}\left(k\right){\stackrel{\‾}{H}}_i\left(k\right),$ 以及

$u_i\left(k\right)={\stackrel{\‾}{H}}_i^T{\left(k\right)\stackrel{\‾}{R}}_i^{-1}\left(k\right)[z_i\left(k\right)-b_i^z\left(k\right)];$

其中：R_i ^-1(k)为加权统计线性化后的测量噪声的方差阵，H_i(k)为加权统计线性化以后的测量矩阵，是由线性化引起的误差项。

4.根据权利要求1所述的分布式非线性滤波方法，其特征是，步骤二中所述的动态协同滤波器是指：

${\hat{u}}_i\left(k\right)={\hat{u}}_i(k-1)+\mathrm{\δ\β}\underset{j\∈{\ℵ}_i}{\mathrm{\Σ}}({\hat{u}}_j(k-1)-{\hat{u}}_i(k-1))+u_i\left(k\right)-u_i(k-1);$

${\hat{U}}_i\left(k\right)={\hat{U}}_i(k-1)+\mathrm{\δ\β}\underset{j\∈{\ℵ}_i}{\mathrm{\Σ}}({\hat{U}}_j(k-1)-{\hat{U}}_i(k-1))+U_i\left(k\right)-U_i(k-1);$

其中：k是采样步数，δ是协同滤波器的更新步长，增益β＞0且β～O(1/λ₂)，λ₂是网络图模型的Laplacian矩阵L的特征值，协同滤波器的输出

和

分别为对其输入U_i(k)或u_i(k)的估计。

Images