CN108566178B - 一种非稳态随机机会网络特征值滤波方法 - Google Patents

Abstract

本发明公开了一种非稳态随机机会网络特征值滤波方法，首先获取不稳定状态下的特征值信号，判断获取的不稳定状态下特征值信号是否为混沌信号；然后针对特征值信号的特性，采用扩展卡尔曼滤波算法进行特征值滤波。本发明在基于扩展卡尔曼滤波进行特征值滤波过程中，依据扩展卡尔曼滤波算法特性，通过对滤波结果的分析，验证滤波方式符合滤波精度要求，解决了非稳态下随机机会网络特征值滤波方式不足的问题。

Description

一种非稳态随机机会网络特征值滤波方法

技术领域

本发明属于移动通信技术领域，涉及一种特征值滤波方法，具体涉及一种采用扩展卡尔曼滤波进行非稳态随机机会网络特征值滤波的方法。

背景技术

随机机会网络(Opportunistic Networks,ON)的核心是节点在进行网络通信时，通信的目的节点与源节点可以不存在完整的路径，取而代之的是利用节点的移动时节点相遇的情况，随机建立相遇节点的通讯连接。而节点之间采取的消息传递方法时通过节点移动随机建立的通讯连接进行逐跳信息传递，其路由模式可以概括为“信息存储—信息携带—信息转发”。随着随机机会网络的研究，其在许多领域，例如传感器设备，车载系统等有着广泛的应用。

机会网络相关技术虽然有广阔的发展前景，但目前仍然面临诸多挑战：①现阶段的机会网络的研究工作缺乏对于通讯节点移动初期，即不稳定状态下相关网络特征的描述，因此研究不可谓之完整。②网络拓扑的相关参数作为信号的研究过程中，由于其包含噪声的因素，对于网络拓扑的刻画会出现偏差。

由此可见，研究随机机会网络必须解决的问题是：对于网络节点初期的不稳定情况应该深入刻画研究。某些拓扑结构相关的全局状态参量完全可以揭示随机机会网络自身的结构特征，比如，通过Fielder向量可以判断出每个网络节点的连接性。由于目前广泛采用一些分布式信号处理算法来实现无基础设施机会网络的数据融合、估算、路由、压缩与检测等任务，这些分布式信号处理算法主要包括分布式估算(包括一致性(consensus)计算、扩散(diffusion)和流言(gossip)算法)，以及最近出现的分布式滤波。而网络拓扑和拓扑结构相关的全局状态参量可以对所运行的分布式信号的相关处理产生重大影响。因此，这就要求本实施例需要选择合适的滤波算法，以解决随机机会网络所面临的关键性问题。

混沌信号(Chaotic Signal,CS)是指确定性动力学系统表现出的不可预测的、类似随机性的信号。混沌现象存在于诸多研究领域，如基础物理学，宇宙天体运动，天气状况预测等。混沌现象研究之初，人们缺乏对于混沌理论的相关了解，错把其当作噪声处理，随着混沌理论的研究，人们开始认识到混沌存在的研究价值。混沌在通讯技术等领域中也有关键技术的不断发展，混沌时间序列分析作为混沌理论的一个重要组成部分在电路系统负载短期预测、网络关键参数滤波、DNA结构序列分析等领域有着广泛的应用。

扩展卡尔曼滤波(Extended Kalman Filter，EKF)是一种常用的滤波工具，其基本原理是通过系统的输入和输出来进行数据观测，对系统的状态进行最优估计的滤波算法。扩展卡尔曼滤波算法自提出以来，它在信息通信，系统的调控与信号滤波等领域发挥了非常重要的作用，成为相关领域最有效、最常见的计算和建模工具之一。扩展卡尔曼滤波的本质就是结合上一状态系统测量值与预测值来估计当前系统的状态参数。经典的卡尔曼滤波只能解决线性系统的滤波问题，对于机会网络中网络节点组成的非线性系统则不再适用。与传统的卡尔曼滤波只能在现行系统发挥作用不同的是，扩展卡尔曼滤波可以处理非线性系统的滤波问题。这主要因为扩展卡尔曼滤波对非线性函数的泰勒展开式进行一阶线性化截取处理，从而将原本的非线性问题转化为线性问题进行处理。这样一来就解决了传统卡尔曼滤波不能应用于非线性系统的问题。扩展卡尔曼滤波应用于非线性系统的状态估计已经的到了相关领域学术界的认可，并广泛使用。

发明内容

本发明针对非稳态下特征值信号的混沌特性，提供了一种特征值滤波方法，采用扩展卡尔曼滤波解决了非稳态下特征值信号滤波方式不足的问题。

本发明所采用的技术方案是：一种非稳态随机机会网络特征值滤波方法，其特征在于，包括以下步骤：

步骤1：获取不稳定状态下的特征值信号，判断获取的不稳定状态下特征值信号是否为混沌信号；

步骤2：针对特征值信号的特性，采用扩展卡尔曼滤波算法进行特征值滤波。

本发明在基于扩展卡尔曼滤波进行特征值滤波过程中，依据扩展卡尔曼滤波算法特性，通过对滤波结果的分析，验证滤波方式符合滤波精度要求，解决了非稳态下随机机会网络特征值滤波方式不足的问题。

附图说明

图1为本发明实施例的圆内单个节点RWP模型示意图；

图2为本发明实施例的扩展卡尔曼滤波递推流程图；

图3为本发明实施例的状态测量值1示意图；

图4为本发明实施例的状态测量值2示意图；

图5为本发明实施例的非稳态特征值信号1示意图；

图6为本发明实施例的非稳态特征值信号2示意图；

图7为本发明实施例的最小嵌入维数示意图；

图8为本发明实施例的关联维数示意图；

图9为本发明实施例的特征值λ₂滤波示意图；

图10为本发明实施例的特征值λ₁₀滤波示意图；

图11为本发明实施例的扩展卡尔曼滤波收敛性分析实例示意图；

图12为本发明实施例的MSE算法判断流程图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明开展的不稳定状态下随机机会网络的特征值滤波，是对于估算网络拓扑结构相关的全局状态信息工作的补充，鉴定随机机会网络在不稳定情况下的特征值信号符合混沌特征，分析系统在不稳定状态下网络拓扑信息中的特征值信号的特性的过程采用计算最大李雅普诺夫指数，鉴定其为混沌信号，并以扩展卡尔曼滤波解决特征值信号滤波方式不足的问题。

本发明的研究对象为随机机会网络的节点移动模型。由于拓扑结构相关的全局状态参数包括特征值、特征值、网络直径等，因此研究中需要选取具有重要价值的参数。本发明中，随机机会网络的“特征值”指的是网络邻接矩阵的特征值，网络邻居矩阵的特征值是一种能够用于标识网络关键信息，帮助实现网重构的关键状态参数。

在实际情况下，机会网络设备之间存在结构的差异，不同通信设备之间的通信方式和通讯距离存在差异。而且由于设备数量有限，通常不可能模拟大型机会网络。因此，本发明利用网络仿真建立节点移动模型。在此基础上，研究了非稳态随机机会网络的特征。

本发明中，针对随机机会网络非稳态下特征值滤波的研究方式采用的是采用网络仿真的方式进行，并不是在真实网络中进行实验。因此，本实施例需要将真实网络抽象为随机机会网络相关运动节点，而节点具有相应的运动模型的特点。根据研究现状和实际存在的问题，本发明研究过程选择了具有“通信”功能的节点移动模型。以RWP模型为例，“随机机会网络”在本发明中具有以下特点：

(1)节点的运动是在一个平面的区域，区域为圆形且有界。

(2)通信节点数目选择大于400个，目的是为了减少随机移动模型本身所具有的移动模式的偶然性。

(3)节点有独立运动的特点。节点之间是相互独立存在的，节点之间的通信方式是双向通信。为了确保机会网络本身的稀疏性特点，网络节点的通信范围与整个运动区域的比值应该小于0.3。

(4)网络节点的运动速度和节点停滞时间在一定范围内是可变化的。

(5)节点的移动准则按照标准RWP模型中节点的运动准则来设置。

采取以上节点的模拟规则的原因，是为了尽量使移动模型贴近于真实的随机网络的节点运动情况，使得后续实验结论能有实际的参考价值。

实体移动模型之间有着各自的特点，因此移动模型的选取，要根据研究的实际，结合研究所需要的实体移动模型的相关性能来进行。表1给出了几种主要实体移动模型的相关特性分析，通过分析移动模型的特性，为之后的滤波工作选取合适的实验模型。

表1实体移动模型分析

SRM模型与GM模型构建相对难度较大，并且具有一定的群组移动模型的特性，因此不适合作为机会网络的节点移动模型。而RWP模型、RW模型与RD模型都是易于构建，且节点的速率与方向与上一时刻相互独立的实体移动模型。理论上，这三种模型都可以作为机会网络的节点移动模型，德国学者Bettstetter在针对RWP模型的研究中，给出了网络特征值的计算公式：

其中n为节点数，r₀为通信半径r与节点活动半径R的比值，p为节点平均移动速度。该公式的计算结果精度较高，但其局限性在于只能计算特征值的均值而非某一时刻的准确值。

基于这一研究现状的不足，为了实现特征值的实时估算，本发明拟选取RWP模型作为机会网络的节点移动模型。通过研究RWP模型特征值的实时变化特性，进而以滤波的方式实现特征值的实时预测。目前针对圆形区域的RWP的研究相对匮乏，因此本发明的网络模型主要以圆形区域的RWP模型为主。图1给出了圆形区域内RWP模型单个节点的运动轨迹图。

若要研究不稳定状态下的机会网络，先必须将不稳定的网络状态与稳定的网络状态进行区别，这是进行不稳定状态下网络网络拓扑关键信息，即特征值信号获取的前提。

在圆形区域内，随机机会网络(RWP模型)在节点分布达到平衡前处于不稳定状态。对于具有N个节点的RWP模型中处于静止状态的节点，其分布特性满足随机分布，设为f(x)；处于运动状态的节点，也满足随机的分布特性，设为g(x)。

现假设t₀时刻，RWP模型中有k个节点处于静止状态，因此有N-k个节点处于运动状态；此时F(x)满足：

其中，N为总节点个数，k为静止的节点个数，f(x)为静止节点的分布函数，g(x)表示运动状态节点的分布函数。

定义处于静止状态与所有节点比值为p，即：

理论情况下，p的值必为定值。现在系统中加入m个新节点(可以看作给系统加入扰动)，加入的m个节点满足均匀分布，因此t₀时刻系统的分布特性满足：

由于加入了新的节点，RWP模型的平衡性已经被打破。在网络区域内，如果只考虑分布特性，可以把原有的N个节点和新加入的m个节点看做是两个独立运行的RWP网络。

定义新加入的m个节点的分布特性为y(x)，当t→∞时，对于这m个节点必有：

y(x)＝m·p·f(x)+m·(1-p)·g(x) (5)；

而原有的N个节点的分布特性保持不变。综上所得，当t→∞时，满足：

将上式带入得：

F_t→∞(x)＝(k+m·p)·f(x)+(N-k+m·(1-p))·g(x) (7)；

将式(3)代入得：

当N→∞且m为定值时，必有：

N+m＝N (9)

因此公式(8)可变形为：

F_t→∞(x)＝k·f(x)+(N-k)·g(x) (10)

由于公式(2)与式(10)具有相同的表达形式，因此当RWP模型中的节点，系统节点的初始移动过程中且达到分布平衡前，系统处于不稳定状态，既基于RWP模型的机会网络仿真过程中，网络节点的分布特性不是固定的。节点在初始时刻满足随机分布，但随着模型的不断运行，节点的分布特性不再满足随机分布，转而趋于满足中心概率密度大，边缘概率密度小的分布。

最后，节点在中心区域最大，而在边界区域的为0，这种分布代表网络达到稳态。故在其达到稳定状态之前，都是不稳定的状态。

由于特征值产生于网络节点的连通度矩阵，因此其代表了节点的代数联通信息，令随机变量X＝(X,Y)，表示某时刻网络节点在系统区域A内的位置，因此节点的空间分布概率密度函数可表示为：

节点X位于A的子区域A¹的概率为在这个区域对f_X(x)进行积分，即满足：

上式可以看做是多个移动节点(n>>1)的整体平均，因此nP可以表示为A¹内任意时刻的节点数。由上式可知，理论情况下，仿真区域内的节点在某个子区域的出现的概率是不固定的。因此，对于仿真区域内的所有节点，其出现在整个仿真区域满足不稳定的分布特性。

需要注意的是，上述证明针对的是理论情况，也就是说理论情况下RWP模型的特征值曲线保持不变。因此如果要获取不稳定情况下的特征值曲线，取在特征值达到常数之前的信号图像即可。

1963年，Lorenz发表的“决定性的非周期流”指出，气候条件变化不能重复和无法长期进行气候预测之间存在相关联系，这就是系统非周期性与系统不可预测性之间的联系。他还发现了混沌现象有对系统初始条件的极端敏感的特性。60年代后期，随着混沌理论的发展，混沌理论研究逐渐深入。1971年，研究者在开展对于耗散系统的研究中正式的引入了混沌理论中奇异吸引子的概念。1975年，李天岩和J.A.Yorke等人提出了对于混沌的科学界定。李天岩和J.A.Yorke指出：在确定性一维系统中，如果出现周期3，则该系统也能对应地出现其他周期，因此能够出现混沌现象。在确定性系统中能够发现混沌现象改变了人们研究宇宙时认为宇宙是可预测系统的观点。利用确定性方程，无法找到稳定的模型，但得到了随机的系统参数，这完全否定了拉普拉斯决定论的可预见性这一观点。但是人们也发现许多过去被定性为噪声的信号实际上可以使用混沌理论进行处理。包含固有规则的“噪声”不同于实际噪声，是可以应用混沌理论进行解释和处理的。

1983年，格里波基通过相关实验证明了只要最大李雅普诺夫(Lyapunovexponent)指数大于零，就可以肯定混沌存在于相关系统中，因此判断一个非线性系统产生的信号是否为混沌信号时，需要测量它的李雅普诺夫指数λ是否为正值，如果测量值为正，则证明信号为混沌信号。李雅普诺夫指数作为一个参数代表了混沌系统的基本特征，表示了系统对初值的极端敏感性这一特性。由两个初始值生成的轨迹与时间成指数关系。李雅普诺夫指数是对这种现象的定量描述。它表示相空间中相邻轨道之间的平均指数收敛率或散度。通过最大李雅普诺夫指数的计算，可以证明混沌存在于动力系统中，如果最大李雅普诺夫指数大于0，则系统存在混沌。一个正的李雅普诺夫指数，表明运动系统的相空间中的两条轨迹无论初始距离等特性关联性如何，差异会随着时间参数变化呈现指数率的增加，如果是一个简单的m维流形(m＝1或m＝2是曲线或面)，则前m李雅普诺夫指数为零，剩下的李雅普诺夫指数为负。不管系统是否具有消散性，只要最大李雅普诺夫指数λ₁＞0，就会出现混沌。

李雅普诺夫指数的计算有着非常多的方式，本发明采用其中的基于相空间重构的计算方法。首先对于前面获取的特征值信号x(n),选取一定的延迟时间τ,(τ＝kΔt,k是正数),构造N_m个m维矢量X_i

X_i＝(x_i,x_i+τ,…,x_i+(m-1)τ) (13)

式中Δt是特征值时间序列采样的时间间隔，i＝1,2,…,N_m,N_m＝N-(m-1)τ.m称为嵌入维数，τ为延迟时间；

重构相空间的目的时正确反映原系统的相关特征，因此需要选取合适的延迟时间τ和相关嵌入维数m，通过计算互信息函数的第一极小值的方式确定延迟时间τ，互信息函数表示的是两个随机变量间普遍性随机关联的数值关系。特征值时间序列x(n)的互信息定义为：

式中,P是时间序列x(n)中互信息所对应的概率分布，Pi是xn落入第i个划分的区间里的概率，P_ij(τ)是x_n落入第i个划分的区间的同时x_(n+τ)落入第j个区间的几率的表示，当I(τ)数值为最小时，即为理论上合适的时间序列延迟τ，随后，从相空间N_m个点里面任意选择一个参考点X_j,计算剩余部分(N_m-1)个点到X_i的相对距离：

其中，i＝1,2,…,N_m；I表示式(14)获取的互信息函数第一极小值；t为特征值时间序列采样的时间；对所有Xi(i＝1,2,…,N_m)采用上述计算方式，可以算取关联积分

式中，θ()表示的是Heaviside函数，且上式满足于以下公式：

θ(x)＝1(x＞0) (7)

用logC_m(r)对log(r)做图形变换，当r的数值趋近于0时，满足C_m(r)∝r^D,因此可以推算得到：

当D的数值不因为相空间维数m的升高而改变时，可以判断其即为关联维数Dc。

在进行实际的实验时,通过画出logC_m(r)-log(r)曲线，不考虑r极小时产生的噪声区域和r极大时C_m(r)中对应的的饱和区域，可以算得中间直线部分的斜率，它表示是本实施例要求的关联维数D。

通过对logC_m(r)-log(r)所对应曲线进行相应线性拟合，通过上述方式计算获得的直线斜率即为关联维数。随后对同一特征值时间序列用数值不同的嵌入维m重构相空间，当D的数值逐渐趋于定值时，即可获到关联维数Dc。与此对应获得的m值便是合适的嵌入维数。

本实施例需要计算的李雅普诺夫指数是表示在相空间区域内相互靠近轨迹的平均算数发散率的相关数值特征。前面通过计算得到了m维重构的相空间。即：

X＝{X_i},i＝1,2,…,N_m (19)

两条起始点相邻较近的曲线L₁和L₂,起始点分别表示为x₀和x₀+Δx₀。以初始值所对应曲线为下一步计算的基准曲线，并以x₀+Δx₀为初始值的曲线为邻近曲线，t时刻基准曲线和邻近曲线上的点为x(x₀,t)和x(x₀+Δx₀,t),记w(x₀,t)＝x(x₀+Δx₀,t)-x(x₀,t)。则两条邻近曲线沿w这一方向的平均算数发散率为：

式中，w0＝w(x₀,0)，在m维相空间中w的全体组成相应跟随相曲线运动的m维切空间。选择此切空间一组算数基底{ei,i＝1,2,…,n}，并将获得的数据由大到小进行排列：λ1≥λ2≥…≥λn，这便是非稳态下特征值时间序列所对应的李雅普诺夫指数。

自从卡尔曼滤波算法被提出以来，它就成为信息通信、控制系统和信号处理中最基本、最常用的计算和建模工具之一。扩展卡尔曼滤波的实质是通过预测数据与实际测量数相结合来估计当前系统的状态参数。

传统卡尔曼滤波只能解决线性系统的滤波问题，由于机会网络节点组成的是非线性系统，所以传统卡尔曼滤波方式无法直接应用于非稳态下的特征值滤波。

扩展卡尔曼滤波(Extended Kalman Filter)与传统卡尔曼滤波只能应用于线性系统的滤波工作不同，它可以用于非线性系统的滤波处理中，图2展示了扩展卡尔曼滤波的相关递推流程。

扩展卡尔曼滤波的基本原理是首先通过上一时刻系统最优预测值和将要施加的控制量来预测当前状态。即：

其中，A_t为状态转移矩阵，B_t为输入控制矩阵，μ_t-1表示上一时刻最优预测值，u_t为本时刻施加的控制量。

计算了当前时刻的预测值之后，接下来需要进行卡尔曼增益的计算，因此我们需要获得预测值所对应的协方差Σ_t：

其中，Σ_t-1表示t-1时刻预测均值，R_t表示预测方差。

上述相关参数计算结束之后，需要计算原始输入信号的相关比例系数：

其中，K_t为卡尔曼增益，C_t为状态转移矩阵，z_t为系统测量值，

为系统状态参数。

通过计算当前时刻的预测值和信号比例系数之后，本实施例获取计算状态估计值所需的条件参数完毕。这时就可以采用下式状态方程算出估计值：

其中，Q_t表示系统预测方差，

表示系统预测均值。

如图3和图4表示了实验中在两组不同特征值信号的状态测量值。

每次预测当前时段的状态信息需要使用上一时段的协方差来进行计算，所以还需要计算当前状态的协方差用于下一时刻的迭代。它的计算过程采用之前计算的卡尔曼增益进行函数构造：

在算法中，从初始状态开始，不断计算当前状态的均值和方差，进行迭代计算。到此为止，上述滤波算法只适用于线性的、满足高斯分布的信号系统，但非稳态下随机机会网络特征值信号经鉴定为混沌信号，系统并不是线性，这时就需要对上述计算过程进行相应的处理：当前状态的概率分布是关于上一状态和将要执行的控制量的二元函数，具体表达式为：

x_t＝g(u_t,x_t-1)+ε_t

z_t＝h(x_t)+δ_t (26)

其中，x_t为t时刻系统的状态，g(u_t,x_t-1)为状态转移矩阵，z_t为系统的观测向量，h(x_t)为观测矩阵，ε_t和δ_t分别为过程噪声和观测噪声。

为了令卡尔曼滤波思想能够解决非线性系统中的状态估计问题，首先要做的就是将g(u_t,x_t-1)和h(x_t)用泰勒级数展开，将其线性化，只取一次项为一阶EKF滤波。具体如式(27)和(28)所示：

g(u_t,x_t-1)≈g(u_t,μ_t-1)+g′(u_t,μ_t-1)(x_t-1-μ_t-1) (27)

g(u_t,x_t-1)在上一状态估计的最优值处取一阶导数，在h(x_t)当前时刻预测值处取一阶导数。这样采用扩展卡尔曼滤波方式进行的非稳态下随机机会网络特征值滤波的滤波器理论构造完毕。

本实施例的滤波实验平台采用了MATLAB软件，MATLAB作为一款被广泛应用的商业数学问题处理软件，在数据分析、矩阵计算和算法模拟等方面有着极佳的表现。使用MATLAB作为随机机会网络相关实验数据的分析与算法编程等方面的研究平台至少具有以下优点：

(1)相比于一般的数值计算工具，MATLAB在数值计算方面有着极高的效率，通过优化数值计算方法可以节省数值计算的时间。

(2)在研究机会网络拓扑结构的全局状态参数时，获得网络邻接矩阵的特征值是一个重要的步骤。MATLAB擅长矩阵计算方面的工作应用，其相关函数能快速准确地求出矩阵的所有特征值，具有非常高的实验效率。

因此，本发明选取了MATLAB工具作为非稳态下随机机会网络特征值滤波实验的平台。

特征值滤波中，测量函数的选取与特征值不同。对于节点数为N的网络邻接矩阵，理论上存在N个特征值，每个特征值的大小虽然不一定相同，但相互之间具有恒定的大小关系。对于一段时间内的不同特征值的变化曲线，总体上具有很大的相关性。基于这一特性，本发明中主要将与被预测值相关的其他多个特征值作为测量值。

特征值滤波中，对于节点数为N的网络邻接矩阵，理论上存在N个特征值，对于一段时间内的不同特征值的变化曲线，总体上具有很大的相关性。如图5至6是便是两组截取的不稳定状态下特征值信号图形。该RWP模型的参数为：300个通信半径为4m的节点，在半径为20m的圆形区域内，移动速度为[5,15](m/s)，停留时间为[8,16](s)。

根据如上所述方法对图7和图8所示的两组典型非稳定状态下特征值信号进行分析，首先采用延时坐标法对特征值进行相空间重构，采用互信息法确定其延迟时间,得到正常情况下的特征值延迟时间为24,18。通过上述方法计算得到两组特征值信号的最小嵌入维数分别是3、4，见图7。两组特征值信号对应计算得到的关联维数分别是2.7561,2.8515,见图8。

继而对于图5和图6所示的两组不稳态下特征值信号进行最大李雅普诺夫指数估算,得到第一组信号为λ₁＝1.4998,而第二组的最大李雅普诺夫指数为λ₂＝1.4088.同样为混沌信号。至此验证不稳定状态下的特征值信号为混沌信号。

扩展卡尔曼滤波算法的一般形式为：假设现在需要预测λ_k，可以将已知的λ_k-1和λ_k+1作为测量值，按一定的比例系数代入滤波过程，滤波过程中的观测噪声的的值同样可以是两个测量值的方差的均值。以300个通信半径为4m的节点，在半径为20m的圆形区域内，速度为[5,15](m/s)，停留时间为[10,20](s)的RWP模型为例，图9至图10给出了扩展卡尔曼波算法对不同特征值预测结果。

卡尔曼滤波本身具有多次迭代之后滤波发散的特性，如前所述，根据扩展卡尔曼滤波算法的滤波特征值提出的，因为预测值与测量范围相关的测量值，可以有效地限制的预测结果，所以需要分析扩展卡尔曼滤波算法的收敛性。

图11给出了半径为20m的圆形区域内节点数为400，节点速度为[5,10](m/s)，节点停滞为[5,15](s)，节点通讯半径为2m的RWP多个时刻特征值的预测结果。

在拟合值的值之前的4次迭代映射，不包含在预测结果中。从图中可以清楚地看出，预测精度可以从第5次迭代到第11次迭代仍然保持，随后的预测结果逐渐开始发散。

为了量化预测结果，表2给出了第五次迭代开始时的特征值测量(拟合曲线)、预测值和实际值，并给出了预测值与实际值之间的绝对误差。滤波算法的收敛性是由绝对误差的变化来判断的。

表2扩展卡尔曼滤波收敛情况

从表中可以明显看出，预测预测值与实际值之间的绝对误差，从时刻开始逐渐变大，后续的绝对误差已远大于0.1，因此可以认为从时刻开始，滤波算法开始发散。传统卡尔曼滤波算法无法有效地解决这一问题，根据大量实验统计，扩展卡尔曼滤波算法一般只能在10次内精确预测，然后逐渐偏离。

可以清楚地看到，预测值和实际值之间的绝对误差增加是从t₁₁时刻进行，随后的绝对误差远大于0.1。因此，可以认为从t₁₂的瞬间，滤波算法开始发散。根据大量实验数据，扩展卡尔曼滤波算法只能准确预测10次，然后逐渐偏离。

为了检验扩展卡尔曼滤波过程的预测精度，本发明选取均方误差(单位：dB)进行误差分析。数理统计领域中均方误差是指:”参数估计值与参数真值之差平方的期望值”，记为MSE。MSE是衡量“平均误差”的一高效且便捷的方法，MSE可以判断数据的变化情况，MSE计算所得的值越小，说明实验相关结果精确度越高，反之则越低。

图12为均方误差法对于实验结果的判断流程。

设n个测量值的误差为ε₁、ε₂……ε_n，则这组测量值的标准误差等于：

对于特征值滤波，需要平衡预测的效率与精度，即在保证一定的预测精度的情况下，尽量减小系统的建模难度和时间复杂度。为了使滤波器具有更好的实用性，本发明尽量简化了滤波过程，因此规定当预测的均方误差不大于-20dB时(不带单位的计算结果为0.01)，预测精度符合要求。

以下将通过数据分析，对比扩展卡尔曼滤波算法与传统滤波算法对同一特征值的预测误差分析。表3表示三个不同特征值滤波实验进行了对应的六组实验进行分析，误差采用上述MSE的计算方法进行计算。

表3特征值预测均方误差对比

实验数据表明，基于非稳态圆形的RWP模型的特征值滤波，扩展卡尔曼滤波算法的预测误差基本满足MSE误差小于-20dB这一滤波精度要求，滤波精度的均值在-28dB左右，通过以上实验可以证明，根据三组特征值滤波实验的精度结果，扩展卡尔曼滤波能够满足滤波实验的精度要求。

综上所得，基于扩展卡尔曼滤波算法的非稳态下随机机会网络特征值滤波，能够在保证对特征值能够实时获取的情况下，滤波结果有较高的预测精度，使得非稳态下随机机会网络特征值滤波方式不足的问题能够得到解决。

特征值滤波的目的之一是实现网络重构，因此基于扩展卡尔曼滤波算法的预测误差分析，主要是对特征值预测误差能否满足矩阵重构的分析。

虽然扩展卡尔曼滤波滤波算法具有较高的预测精度，但是由于预测值与实际值必然存在一定的误差，并且每个特征值的预测误差也存在一定的差异，因此这里需要通过误差分析来判断改进算法是否满足网络精确重构的条件。公式定义如下：

其中Δλ_i表示第i个特征值的扰动值。根据上式可知，单个特征值预测误差均值

应满足：

由式(31)可得：

对于随机机会网络的特征值实时预估，特征值的波动性越大，特征值滤波的预测精度也会随之降低，因此一般情况下无法直接分析。前面已经讨论了所有特征值的方差变化特性，可以看出，靠近两边的特征值方差较大。根据计算，扩展卡尔曼滤波对于单个特征值预测的绝对误差为3×10^-3，因此10个特征值的误差之和为0.03。

接下来计算其他特征值的预测误差，对于其他的特征值，由于方差较小，因此采用扩展卡尔曼滤波算法依然可以保证较高的预测精度。根据对网络中所有特征值的方差分析研究，结果表明，方差较大的一部分特征值的方差趋于满足幂律分布，假设这部分特征值的方差为σ，特征值序号为x，方差与特征值序号之间的关系可表示为：

lg(σ)＝alg(x)+b (33)

其中ɑ和b为系数。而其他方差较小特征值的方差趋于满足线性分布，方差与特征值序号之间的关系可表示为：

σ＝cx+d (34)

其中c和d为系数。由于特征值的预测误差与其方差总体上成正比，对于扩展卡尔曼滤波算法，针对方差最大的特征值预测误差约为8×10^-3，方差第100大的特征值预测误差约为8×10^-4，将两组参数代入式(33)求取系数，再反过来计算这50个特征值的预测误差之和，计算结果约为0.114。最后根据式(34)计算可得，当节点数不大于1000时，预测误差可以满足矩阵重构的要求。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (4)

1.一种非稳态随机机会网络特征值滤波方法，应用于移动通信技术领域；其特征在于，包括以下步骤：

步骤1：获取不稳定状态下的特征值信号，判断获取的不稳定状态下特征值信号是否为混沌信号；

所述特征值是网络邻接矩阵的特征值；

步骤1的具体实现包括以下子步骤：

步骤1.1：搭建具有“通信”功能的节点移动模型，对随机机会网络不稳定状态提出具体界定，获取不稳定状态下的特征值信号；

步骤1.2：通过最大李雅普诺夫指数的计算，判断获取的不稳定状态下特征值信号是否为混沌信号；

其中，判断李雅普诺夫指数λ是否为正值，如果测量值为正，则信号为混沌信号；所述李雅普诺夫指数λ，其计算过程为：

首先对于获取的特征值信号x(n)，选取一定的延迟时间τ，构造N_m个m维矢量X_i：

X_i＝(x_i,x_i+τ,…,x_i+(m-1)τ) (1)

式中，τ＝kΔt，k是正数，Δt是特征值时间序列采样的时间间隔，i＝1,2,…,N_m,N_m＝N-(m-1)τ，m为嵌入维数，τ为延迟时间；

选取合适的延迟时间τ和相关嵌入维数m，通过计算互信息函数的第一极小值的方式确定延迟时间τ，互信息函数表示的是两个随机变量间普遍性随机关联的数值关系；特征值时间序列x(n)的互信息定义为：

式中,P是时间序列x(n)中互信息所对应的概率分布，Pi是x_n落入第i个划分的区间里的概率，P_ij(τ)是x_n落入第i个划分的区间的同时x_(n+τ)落入第j个区间的几率的表示；

当I(τ)数值为最小时，为理论上合适的时间序列延迟τ；随后从相空间N_m个点里面任意选择一个参考点X_j，计算剩余部分(N_m-1)个点到X_i的相对距离：

其中，i＝1,2,…,N_m；I表示式(2)获取的互信息函数第一极小值；t为特征值时间序列采样的时间；对所有X_i采用上述计算方式，算取关联积分：

式中，θ()表示的是Heaviside函数，且当x＞0时：

θ(x)＝1 (5)

用logC_m(r)对log(r)做图形变换，当r的数值趋近于0时，满足C_m(r)∝r^D；因此推算得到：

当D的数值不因为相空间维数m的升高而改变时，则判断其为关联维数Dc；

通过画出logC_m(r)-log(r)曲线，不考虑r极小时产生的噪声区域和r极大时C_m(r)中对应的的饱和区域，算得中间直线部分的斜率，为关联维数D；

通过对logC_m(r)-log(r)所对应曲线进行相应线性拟合，同时计算获得的直线斜率即为关联维数；随后对同一特征值时间序列用数值不同的嵌入维m重构相空间，当D的数值逐渐趋于定值时，即获到关联维数Dc；与此对应获得的m值便是合适的嵌入维数；

前面计算得到了m维重构的相空间，即：

X＝{X_i},i＝1,2,…,N_m (7)

取两条起始点相邻较近的曲线L₁和L₂，起始点分别表示为x₀和x₀+Δx₀，以初始值所对应曲线为下一步计算的基准曲线，并以x₀+Δx₀为初始值的曲线为邻近曲线，t时刻基准曲线和邻近曲线上的点为x(x₀,t)和x(x₀+Δx₀,t)，记w(x₀,t)＝x(x₀+Δx₀,t)-x(x₀,t)；则两条邻近曲线沿w这一方向的平均算数发散率为：

式中，w0＝w(x₀,0)，在m维相空间中w的全体组成相应跟随相曲线运动的m维切空间；选择此切空间一组算数基底{ei,i＝1,2,…,n}，并将获得的数据由大到小进行排列：λ₁≥λ₂≥…≥λ_n，这便是非稳态下特征值时间序列所对应的李雅普诺夫指数；

步骤2：针对特征值信号的特性，采用扩展卡尔曼滤波算法进行特征值滤波。

2.根据权利要求1所述的非稳态随机机会网络特征值滤波方法，其特征在于：步骤1.1中所述具有“通信”功能的节点移动模型为RWP模型，所述RWP模型具有以下特点：

特点1：节点的运动是在一个平面的区域，区域为圆形且有界；

特点2：通信节点数大于400个；

特点3：节点有独立运动的特点；节点之间是相互独立存在的，节点之间的通信方式是双向通信；网络节点的通信范围与整个运动区域的比值小于0.3；

特点4：网络节点的运动速度和节点停滞时间在一定范围内是可变化的；

特点5：节点的移动准则按照标准RWP模型中节点的运动准则来设置。

3.根据权利要求1所述的非稳态随机机会网络特征值滤波方法，其特征在于：步骤1.1中所述获取不稳定状态下的特征值信号，是在特征值曲线中，取在特征值达到常数之前的信号。

4.根据权利要求1-3任意一项所述的非稳态随机机会网络特征值滤波方法，其特征在于，步骤2中所述扩展卡尔曼滤波算法，具体实现包括以下子步骤：

步骤2.1：通过上一时刻最优预测值和将要施加的控制量来预测当前状态，即：

其中，A_t为状态转移矩阵，B_t为输入控制矩阵，μ_t-1表示上一时刻最优预测值，u_t为本时刻施加的控制量；

步骤2.2：获得预测值所对应的协方差

其中，Σ_t-1表示t-1时刻预测均值，R_t表示预测方差；

步骤2.3：计算原始输入信号的相关比例系数：

其中，K_t为卡尔曼增益，C_t为状态转移矩阵，z_t为系统测量值，

为系统状态参数；

步骤2.4：采用下式状态方程计算估计值：

其中Q_t表示系统预测方差，

表示系统预测均值；

步骤2.5：计算当前状态的协方差用于下一时刻的迭代；采用卡尔曼增益进行函数构造：

步骤2.6：令当前状态的概率分布是关于上一状态和将要执行的控制量的二元函数：

x_t＝g(u_t,x_t-1)+ε_t

z_t＝h(x_t)+δ_t (14)

其中，x_t为t时刻系统的状态，g(u_t,x_t-1)为状态转移矩阵，z_t为系统的观测向量，h(x_t)为观测矩阵，ε_t和δ_t分别为过程噪声和观测噪声；

将g(u_t,x_t-1)和h(x_t)用泰勒级数展开，将其线性化，只取一次项为一阶EKF滤波，具体如式(15)和(16)所示：

g(u_t,x_t-1)≈g(u_t,μ_t-1)+g′(u_t,μ_t-1)(x_t-1-μ_t-1) (15)

其中，g(u_t,x_t-1)在上一状态估计的最优值处取一阶导数，h(x_t)为当前时刻预测值处取一阶导数；这样采用扩展卡尔曼滤波方式进行的非稳态下随机机会网络特征值滤波的滤波器构造完毕。

Images