md_files/科研/陈茂森论文.md

陈茂森论文

随机移动网络系统的稳定性

马尔科夫链与网络平均度推导

1.马尔科夫链的基本概念

马尔科夫链描述的是这样一种随机过程：系统在若干个可能的状态中变化，下一时刻所处状态只依赖于当前状态，而与过去的状态无关，这就是所谓的“无记忆性”或马尔科夫性。

无记忆性意味着，对于任何 $s, t \ge 0$，

P(T > s+t \mid T > s) = P(T > t).

假设你已经等待了 s 分钟，那么再等待至少 t 分钟的概率，和你一开始就等待至少 t 分钟的概率完全相同。

在所有概率分布里，只有指数分布

P(T>t) = e^{-\lambda t}

具有这种“无记忆性”特征：

P(T>s+t \mid T>s) = \frac{P(T>s+t)}{P(T>s)} = \frac{e^{-\lambda(s+t)}}{e^{-\lambda s}} = e^{-\lambda t} = P(T>t).

链路状态的马尔科夫模型

考虑网络中每条链路的动态行为，其状态空间为：

• 状态0：链路断开
• 状态1：链路连通

定义概率函数：

• $p_1(t)$：时刻 t 处于连通状态的概率
• $p_0(t) = 1 - p_1(t)$：断开概率

同时，我们假设链路从一个状态转移到另一个状态需要等待一段时间，这段等待时间通常服从指数分布（论文中通过 KS 检验确认）：

• 从断开（0）到连通（1）的等待时间 T_{01} \sim \text{Exp}(\lambda_{01})
• 从连通（1）到断开（0）的等待时间 T_{10} \sim \text{Exp}(\lambda_{10})

其中，\lambda_{01} 和 \lambda_{10} 为转移速率，表示单位时间内事件（转移）发生的平均次数

2.推导单条链路的连通概率

根据连续时间马尔科夫链的理论，我们可以写出状态转移的微分方程。对于状态1（连通状态），概率 p_1(t) 的变化率由两个部分组成：

1. 当链路处于状态0时，以速率 \lambda_{01} 变为状态1。这部分概率增加的速率为

   
   \lambda_{01} \, p_0(t)=\lambda_{01} (1-p_1(t)).
   

2. 当链路处于状态1时，以速率 \lambda_{10} 转换为状态0。这部分使 p_1(t) 减少，其速率为

   
   \lambda_{10} \, p_1(t).
   

所以，p_1(t) 的微分方程写成：

\frac{d p_1(t)}{dt} = \lambda_{01} \, (1-p_1(t)) - \lambda_{10} \, p_1(t).

这个方程可以整理为：

\frac{d p_1(t)}{dt} + (\lambda_{01}+\lambda_{10}) \, p_1(t) = \lambda_{01}.

这其实是一个一阶线性微分方程，其标准求解方法是求解其齐次解与非齐次解。

3. 求解微分方程

整个微分方程的通解为：

p_1(t)= \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} + C\, e^{-(\lambda_{01}+\lambda_{10})t}.

利用初始条件 $p_1(0)=p_1^0$（初始时刻链路连通的概率），我们可以求出 $C$：

即

C = p_1^0 -\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}}.

所以，链路在任意时刻 t 连通的概率为：

p_1(t)= \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} + \left( p_1^0 -\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} \right)e^{-(\lambda_{01}+\lambda_{10})t}.

这就是单条链路的连通概率函数，描述了从任意初始条件出发，经过一段时间后，链路达到平衡状态的过程。

4.推导网络平均度的变化函数

在一个由 N 个节点构成的网络中，每个节点都与其它节点进行通信（不考虑自环），因此每个节点最多有 N-1 个邻居。对于任意一对节点 i 和 $j$，它们之间链路连通的概率 $p_1(t)$（假设所有链路独立且同分布）。

• 某个节点 i 在时刻 t 的度 $d_i(t)$可以写作：

  
  d_i(t)= \sum_{\substack{j=1 \\ j\neq i}}^N p_{ij}(t),
  

  其中 $p_{ij}(t)=p_1(t)$。

• 因此，每个节点的期望度为：

  
  E[d_i(t)]=(N-1)p_1(t).
  

• 网络平均度就是对所有节点的期望度取平均，由于网络中每个节点都遵循相同统计规律，所以网络平均度可表示为：

  
  \bar{d}(t) = \frac{1}{N}\sum_{i=1}^{N} E[d_i(t)] = \frac{N \cdot (N-1)p_1(t)}{N} = (N-1)p_1(t)
  

将我们前面得到的 p_1(t) 表达式代入，就得到网络平均度随时间变化的表达式：

\bar{d}(t)= (N-1)\left[ \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} + \left( p_1^0 -\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}} \right)e^{-(\lambda_{01}+\lambda_{10})t}\right].

这就是网络平均度的变化函数：

• 网络开始时每条链路的连通概率为 p_1^0
• 这种调整过程符合指数衰减规律，即偏离平衡值的部分按 e^{-(\lambda_{01}+\lambda_{10})t} 衰减
• 当 t 趋向无穷大时，指数项 e^{-(\lambda_{01}+\lambda_{10})t} 衰减为0，网络平均度趋向于 $(N-1)\frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}}$，这就是网络达到平衡状态后的理论平均度。

特征信号参数的平稳性

证明系统在平衡态下具有统计上的稳定性。

从节点空间分布证明平稳性。

设节点在模型区域的坐标为 $(X,Y)$，其分布概率密度函数写为

f(x,y).

那么节点在模型子区域 R_1 中出现的概率为

P_{R_1}=\int_{R_1} f(x,y) \,dx\,dy.

在平衡状态下，理论上节点的位置分布 f(x,y) 保持不变，即每个区域内节点出现的概率 P_{R_1} 是常数，不随时间变化。证明在平衡状态下节点分布稳定。

扰动后的恢复能力

• 静止节点分布特性满足均匀分布,概率密度函数为 $g(x,y)$；
• 运动节点的概率密度函数为 $h(x,y)$；
• 在时刻 t_0 时，网络中共有 N 个节点，其中有 s 个静止，故静止节点的比例为 $p=\frac{s}{N}$。

节点整体的分布概率密度函数可写为

f(x,y)=p\, g(x,y)+(1-p)\, h(x,y).

在平衡状态下，p 的理论值为一个常数，所以 f(x,y) 不随时间变化，从而网络连通度稳定。

接下来，考虑外界扰动的影响：假设在时刻 $t_1$ 新加入 m 个符合均匀分布的节点，

扰动后的总分布（$t_1$时刻后）

• 新加入的 m 个节点是静止的，其分布为 g(x,y)
• 此时网络的总节点数 $N+m$：
  • 静止节点总数：s+m
  • 运动节点总数：$N-s$（原有运动节点数不变）

因此，扰动后的分布为：

f(x,y,t_1) = \frac{s + m}{N + m} g(x,y) + \frac{N - s}{N + m} h(x,y)

f(x,y,t_1) = p' \cdot g(x,y) + (1-p') \cdot h(x,y)

其中 $p' = \frac{s + m}{N + m}$。

近似处理（当 N, s \gg m 时）

p' = \frac{s + m}{N + m} \approx \frac{s}{N} = p

因此，扰动后的分布近似为：

f(x,y,t_1) \approx p \cdot g(x,y) + (1-p) \cdot h(x,y)

这与初始平衡态的分布相同，说明网络在扰动后恢复了平衡态。

系统稳定性分析

平衡点及误差坐标

论文第 2.1 节推导出，单条链路连通概率 p_1(t) 满足

\dot p_1(t)
  = -(\lambda_{01}+\lambda_{10})\,p_1(t) \;+\;\lambda_{01}.
  \tag{2‑18}

网络有 N 个节点，平均度

d(t) = (N-1)\,p_1(t).

设平衡连通概率 p_1^* 满足 $\dot p_1=0$，解得 平衡点

p_1^* = \frac{\lambda_{01}}{\lambda_{01}+\lambda_{10}},
  \quad
  d^* \;=\;(N-1)\,p_1^*.

定义误差（偏离平衡的量）

e(t)=d(t)-d^*.

其中

d(t)=(N-1)\,p_1(t),
\qquad
d^*=(N-1)\,p_1^*.

将 d 和 d^* 代入

e(t)
=d(t)-d^*
=(N-1)\,p_1(t)\;-\;(N-1)\,p_1^*
=(N-1)\,\bigl[p_1(t)-p_1^*\bigr].

解得

p_1(t)-p_1^* \;=\;\frac{e(t)}{\,N-1\,}
\quad\Longrightarrow\quad
p_1(t)
=\frac{e(t)}{\,N-1\,}+p_1^*.

误差求导

\dot e(t) = \frac{d}{dt}\bigl[(N-1)(p_1-p_1^*)\bigr]
            = (N-1)\,\dot p_1(t),

得到

\dot e
  =(N-1)\Bigl[-(\lambda_{01}+\lambda_{10})\Bigl(\tfrac{e}{N-1}+p_1^*\Bigr)
              +\lambda_{01}\Bigr].\\\dot e = -(\lambda_{01}+\lambda_{10})\,e.

这就是把原来以 p_1 为自变量的微分方程，转写成以 "偏离平衡量" e 为自变量的形式

记常数

  c = \lambda_{01}+\lambda_{10} >0,

则误差模型就是一维线性常微分方程：

  \dot e = -\,c\,e.

构造李雅普诺夫函数

对一维系统 $\dot e=-ce$（$c>0$），自然选取

V(e)=e^2

作为李雅普诺夫函数，理由是：

• V(e)>0 当且仅当 $e\neq0$；
• 平衡点 e=0 时，$V(0)=0$。

计算 V 的时间导数

对 V 关于时间求导：

\dot V(e)
= \frac{d}{dt}\bigl(e^2\bigr)
= 2\,e\,\dot e
= 2\,e\,\bigl(-c\,e\bigr)
= -2c\,e^2.

因为 c>0 且 $e^2\ge0$，所以

\boxed{\dot V(e)\;=\;-2c\,e^2\;\le\;0.}

• 当 e\neq0 时，$\dot V<0$；
• 当 e=0 时，$\dot V=0$。

这正是“半负定”（negative semi-definite）的定义。

结论

李雅普诺夫第二类定理告诉我们：

若存在一个函数 V(e) 在平衡点处为 0、在邻域内正定，且其导数 \dot V(e) 在该邻域内为半负定，则平衡点 $e=0$（即 $d=d^*$）是稳定的。

由于我们已经构造了满足上述条件的 $V(e)=e^2$，并验证了 $\dot V(e)\le0$，故平衡态 d=d^* 是 李雅普诺夫意义下稳定 的。

网络特征谱参数的估算

由于邻接矩阵不能保证半正定性，因此会产生幂迭代估算过程不能收敛的问题。需构造A^T A

基于奇异值分解改进幂迭代估算（集中式）

输入：矩阵 $B = A^T A$，目标特征值数量 $k$，收敛阈值 \delta
输出：前 k 个特征值 \lambda_1' \geq \lambda_2' \geq \dots \geq \lambda_k' 及对应特征向量 u_1', u_2', \dots, u_k'

1. 初始化

1. 随机生成初始非零向量 $v^{(0)}$，归一化：

   
   v^{(0)} \gets \frac{v^{(0)}}{\|v^{(0)}\|_2}
   

2. 设置已求得的特征值数量 $n \gets 0$，剩余矩阵 B_{\text{res}} \gets B

2. 迭代求前k个特征值与特征向量

While n < k:

1. 幂迭代求当前最大特征值与特征向量

   • 初始化向量 $v^{(0)}$（若 $n=0$，用随机向量；否则用与已求特征向量正交的向量）

   • Repeat: a. 计算 $v^{(t+1)} \gets B_{\text{res}} v^{(t)}$​ b. 归一化：

     
     v^{(t+1)}\gets \frac{v^{(t+1)}}{\|v^{(t+1)}\|_2}
     

     c. 计算 Rayleigh 商：

     
     y^{(t)} = \frac{(v^{(t)})^T B_{\text{res}} v^{(t)}}{(v^{(t)})^T v^{(t)}}
     

     d. Until $|y^{(t)} - y^{(t-1)}| < \delta$（收敛）

   • 记录当前特征值与特征向量：

     
     \lambda_{n+1}' \gets y^{(t)}, \quad u_{n+1}' \gets v^{(t)}
     

2. 收缩矩阵以移除已求特征分量 每次收缩操作将已求得的特征值从矩阵中“移除”，使得剩余矩阵的谱（特征值集合）中次大特征值“升级”为最大特征值。

   • 更新剩余矩阵：

     
     B_{\text{res}} \gets B_{\text{res}} - \lambda_{n+1}' u_{n+1}' (u_{n+1}')^T
     

   • 确保 B_{\text{res}} 的对称性（数值修正）

3. 增量计数

   • n \gets n + 1

瑞利商公式

1. 集中式：

   
   y(k)= \frac{x(k)^T A x(k)}{x(k)^T x(k)}
   

2. 分布式一致性计算：

   
   y(k) = \frac{\sum_{i=1}^N x_i(k) b_i(k)}{\sum_{i=1}^N x_i^2(k)}
   

   其中

   
   b_i(k) = \sum_j a_{ij} x_j(k)
   

两者是等价的：

考虑一个简单的2×2矩阵

A = \begin{pmatrix}2 & 1\\1 & 2\end{pmatrix}, \quad x = \begin{pmatrix}2\\1\end{pmatrix}.

1. 集中式计算

y= \frac{x^T A x}{x^T x} 
= \frac{\begin{pmatrix}2 & 1\end{pmatrix} \begin{pmatrix}5\\4\end{pmatrix}}{2^2 + 1^2} 
= \frac{10 + 4}{5} 
= \frac{14}{5} = 2.8.

2. 分布式计算
   各节点分别计算本地观测值

​ 节点1的计算：

b_1 = a_{11}x_1 + a_{12}x_2 = 2 \cdot 2 + 1 \cdot 1 = 5.

​ 节点2的计算：

b_2 = a_{21}x_1 + a_{22}x_2 = 1 \cdot 2 + 2 \cdot 1 = 4.

​ 然后通过全网共识计算

y = \frac{2 \cdot 5 + 1 \cdot 4}{2^2 + 1^2} 
= \frac{14}{5} = 2.8.

主要符号表

符号	类型	含义	存储/计算位置
n	下标	当前计算的奇异值序号（从0开始）	全局共识
K	常量	需要计算的前$K$大奇异值总数	预设参数
j,k	下标	节点编号（$j$表示当前节点）	本地存储
𝒩_j	集合	节点$j$的邻居节点集合	本地拓扑信息
a_{jk}	矩阵元素	邻接矩阵$A$中节点$j$与$k$的连接权值	节点$j$本地存储
v_{n,j}^{(t)}	向量分量	第$n$个右奇异向量在节点$j$的分量（第$t$次迭代）	节点$j$存储
u_{n,j}	向量分量	第$n$个左奇异向量在节点$j$的分量	节点$j$计算存储
\sigma_n	标量	第$n$个奇异值	全局共识存储
\delta	标量	收敛阈值	预设参数

分布式幂迭代求前$K$大奇异对

While n < K:

1. 初始化：

   • 若 $n = 0$：

     • 各节点$j$随机初始化 v_{0,j}^{(0)} \sim \mathcal{N}(0,1)

   • 若 $n > 0$：

     • 分布式Gram-Schmidt正交化：

       
       v_{n,j}^{(0)} \gets v_{n,j}^{(0)} - \sum_{m=0}^{n-1} \underbrace{\text{Consensus}\left(\sum_k v_{m,k} v_{n,k}^{(0)}\right)}_{\text{全局内积}\langle v_m, v_n^{(0)} \rangle} v_{m,j}
       

     • 分布式归一化：

       
       v_{n,j}^{(0)} \gets \frac{v_{n,j}^{(0)}}{\sqrt{\text{Consensus}\left(\sum_k (v_{n,k}^{(0)})^2\right)}}
       

2. 迭代计算：

   • Repeat： a. 第一轮通信（计算$z=Av$）：

     
     z_j^{(t)} = \sum_{k \in 𝒩_j} a_{jk} v_{n,k}^{(t)} \quad \text{(邻居交换$v_{n,k}^{(t)}$)}
     

     b. 第二轮通信（计算$y=A^T z$）：

     
        y_j^{(t+1)} = \sum_{k \in 𝒩_j} a_{kj} z_k^{(t)} \quad \text{(邻居交换$z_k^{(t)}$)}
     

     c. 隐式收缩（$n>0$时）：

     
        y_j^{(t+1)} \gets y_j^{(t+1)} - \sum_{m=0}^{n-1} \sigma_m^2 v_{m,j} \cdot \underbrace{\text{Consensus}\left(\sum_k v_{m,k} y_k^{(t+1)}\right)}_{\text{投影系数计算}}
     

     d. 归一化：

     
        v_{n,j}^{(t+1)} = \frac{y_j^{(t+1)}}{\sqrt{\text{Consensus}\left(\sum_k (y_k^{(t+1)})^2\right)}}
     

     e. 计算Rayleigh商：

     
     \lambda^{(t)} = \text{Consensus}\left(\sum_k v_{n,k}^{(t)} y_k^{(t+1)}\right)
     

     f. 终止条件：

     
     \text{If } \frac{|\lambda^{(t)} - \lambda^{(t-1)}|}{|\lambda^{(t)}|} < \delta \text{ then break}
     

3. 保存结果：

   
   \sigma_n = \sqrt{\lambda^{(\text{final})}}, \quad v_{n,j} = v_{n,j}^{(\text{final})}
   

   • 所有节点同步 n \gets n + 1

分布式计算左奇异向量u_{n,j}

对于邻接矩阵 $A \in \mathbb{R}^{N \times N}$，其奇异值分解为：

A = U \Sigma V^T

其中：

• U 的列向量 \{u_n\} 是左奇异向量
• V 的列向量 \{v_n\} 是右奇异向量
• \Sigma 是对角矩阵，元素 \sigma_n 为奇异值

左奇异向量的定义关系：

A v_n = \sigma_n u_n \quad \Rightarrow \quad u_n = \frac{1}{\sigma_n} A v_n

展开为分量形式（对第 j 个分量）：

u_{n,j} = \frac{1}{\sigma_n} \sum_{k=1}^N a_{jk} v_{n,k}

输入：\sigma_n, $v_{n,j}$（来自幂迭代最终结果）

For n = 0 to $K-1$：

1. 本地计算：

   
   u_{n,j} = \frac{1}{\sigma_n} \sum_{k \in 𝒩_j} a_{jk} v_{n,k} \quad \text{(需邻居节点发送$v_{n,k}$)}
   

2. 正交归一化：

   • For m = 0 to $n-1$：

     
     u_{n,j} \gets u_{n,j} - \text{Consensus}\left(\sum_k u_{m,k} u_{n,k}\right) \cdot u_{m,j}
     

   • 归一化：

     
     u_{n,j} \gets \frac{u_{n,j}}{\sqrt{\text{Consensus}\left(\sum_k u_{n,k}^2\right)}}
     

分布式重构邻接矩阵A

输入：\sigma_n, u_{n,j}, v_{n,k}

For 每个节点$j$并行执行：

1. 对每个邻居$k \in 𝒩_j$：

   • 请求节点$k$发送$v_{n,k}$（$n=0,...,K-1$）

   • 计算：

     
     a_{jk} = \sum_{n=0}^{K-1} \sigma_n u_{n,j} v_{n,k}
     

2. 非邻居元素：

   
   a_{jk} = 0 \quad \text{for} \quad k \notin 𝒩_j
   

非稳态下动态特征参数的估算

一致性控制策略

1. 异步更新模型

   • 节点仅在离散时刻 t_k^i 接收邻居信息，更新自身状态 $x_i(t)$。
   • 各节点的状态更新时刻是独立的

2. 延时处理

   • 若检测到延时，节点选择最新收到的邻居状态替代旧值（避免使用过期数据）。

3. 一致性协议设计

   • 无时延系统

     
     \dot{x}_i = \sum_{j \in N(t_k^i, i)} a_{ij}(t_k^i) \left( x_j(t_k^i) - x_i(t) \right)
     

     参数	含义
     \dot{x}_i	节点 i 的状态变化率（导数），表示 x_i 随时间的变化速度。
     x_i(t)	节点 i 在时刻 t 的本地状态值（如特征估计、传感器数据等）。
     x_j(t_k^i)	节点 i 在 t_k^i 时刻收到的邻居节点 j 的状态值。
     N(t_k^i, i)	节点 i 在 t_k^i 时刻的邻居集合（可直接通信的节点）。
     a_{ij}(t_k^i)	权重因子，控制邻居 j 对节点 i 的影响权重，满足 $\sum_j a_{ij} = 1$。

   • 有时延系统

     
     \dot{x}_i = \sum_{j \in N(t_k^i, i)} a_{ij}(t_k^i) \left( x_j(t_k^i - \tau_{ij}^k) - x_i(t) \right)
     

     参数	含义
     \tau_{ij}^k	节点 j 到 i 在时刻 t_k^i 的信息传输延时。
     t_k^i - \tau_{ij}^k	节点 i 实际使用的邻居状态 x_j 的有效时刻（扣除延时）。

   • 权重\alpha_{ij}

     
     \text{有有效邻居时 } a_{ij}(t) = \begin{cases} 
     \frac{\alpha_{ij}(t_k^i)}{\sum_{s \in N(t_k^i, i)} \alpha_{is}(t_k^i)}, & \text{若 } j \in N(t_k^i, i) \\
     0, & \text{若 } j \notin N(t_k^i, i)
     \end{cases}
     

     
     \text{无有效邻居时 } a_{ij}(t) = \begin{cases} 
     1, & \text{若 } j = i \\
     0, & \text{若 } j \neq i
     \end{cases}
     

4. 通信拓扑定义

   • 引入 $G^0(t)$：实际成功通信的瞬时拓扑（非理想链路 $G(t)$），强调有效信息传递而非物理连通性。

收敛性分析

动态网络的收敛条件：

• 在节点移动导致的异步通信和随机延时下，只要网络拓扑满足有限时间内的联合连通性（即时间窗口内信息能传递到全网），所有节点的状态 $x_i$ 最终会收敛到同一全局值。 （平均代数连通度 > 0 == 动态网络拓扑的平均拉普拉斯矩阵的第二小特征值>0 ）

• 无需时刻连通：允许瞬时断连，但长期需保证信息能通过动态链路传递。

基于 UKF 的滤波估算

	KF	EKF	UKF
线性要求	严格线性	弱非线性	强非线性
可微要求	-	必须可微	不要求
计算复杂度	低	中	中
适用场景	线性系统	平滑非线性	剧烈非线性

本文基于UKF：

• 采用**确定性采样（Sigma点）**直接近似非线性分布
• 完全规避对 f(x) 和 h(x) 的求导需求
• 保持高斯系统假设
• 允许函数不连续/不可微
• 适应拓扑突变等非线性情况

UKF 具体步骤

符号说明

• $i$: 节点索引，N 为总节点数
• $x_i(k)$: 节点 i 在时刻 k 的状态分量 (x)
• $b_i(k)$: 节点 i 的本地状态估计值 （相当于$Ax$）
• $a_{ij}$: 邻接矩阵元素（链路权重）
• $Q_k, R_k$: 过程噪声与观测噪声协方差
• $\mathcal{X}_{i,j}$: 节点 i 的第 j 个 Sigma 点
• $W_j^{(m)}, W_j^{(c)}$: Sigma 点权重（均值和协方差）

Step 1: 分布式初始化

1. 节点状态初始化：
   • 每个节点 i 随机生成初始状态分量 $x_i(0)$。
   • 本地状态估计 b_i(0) 初始化为 $x_i(0)$。

---

Step 2: 生成 Sigma 点（确定性采样）

在每个节点本地执行：

1. 计算 Sigma 点：

   
   \begin{aligned}
   \mathcal{X}_{i,0} &= \hat{b}_{i,k-1} \\
   \mathcal{X}_{i,j} &= \hat{b}_{i,k-1} + \left( \sqrt{(n+\lambda) P_{i,k-1}} \right)_j \quad (j=1,\dots,n) \\
   \mathcal{X}_{i,j+n} &= \hat{b}_{i,k-1} - \left( \sqrt{(n+\lambda) P_{i,k-1}} \right)_j \quad (j=1,\dots,n)
   \end{aligned}
   

   • $\lambda = \alpha^2 (n + \kappa) - n$（缩放因子，\alpha 控制分布范围，\kappa 通常取 0）
   • $\sqrt{(n+\lambda) P}$ 为协方差矩阵的平方根（如 Cholesky 分解）

2. 计算 Sigma 点权重：

   
   \begin{aligned}
   W_0^{(m)} &= \frac{\lambda}{n + \lambda} \quad &\text{(中心点均值权重)} \\
   W_0^{(c)} &= \frac{\lambda}{n + \lambda} + (1 - \alpha^2 + \beta) \quad &\text{(中心点协方差权重)} \\
   W_j^{(m)} = W_j^{(c)} &= \frac{1}{2(n + \lambda)} \quad (j=1,\dots,2n) \quad &\text{(对称点权重)}
   \end{aligned}
   

   • $\beta$ 为高阶矩调节参数（高斯分布时取 2 最优）

---

Step 3: 预测步骤（时间更新）

1. 传播 Sigma 点：

   
   \mathcal{X}_{i,j,k|k-1}^* = f(\mathcal{X}_{i,j,k-1}) + q_k \quad (j=0,\dots,2n)
   

   • $f(\cdot)$ 为非线性状态转移函数
   • $q_k$ 为过程噪声 ，反映网络拓扑动态变化（如节点移动导致的链路扰动）。

2. 计算预测均值和协方差：

   
   \hat{b}_{i,k|k-1} = \sum_{j=0}^{2n} W_j^{(m)} \mathcal{X}_{i,j,k|k-1}^*
   

   
   P_{i,k|k-1} = \sum_{j=0}^{2n} W_j^{(c)} \left( \mathcal{X}_{i,j,k|k-1}^* - \hat{b}_{i,k|k-1} \right) \left( \mathcal{X}_{i,j,k|k-1}^* - \hat{b}_{i,k|k-1} \right)^T + Q_k
   

   • $Q_k$ 为过程噪声协方差

---

Step 4: 分布式观测生成

1. 邻居状态融合：

   • 节点 i 从邻居 j 获取其本地观测值 $ b_{j,\text{local}}(k) $：

   
   b_{j,\text{local}}(k) = \sum_{l=1}^N a_{jl} x_l(k) \quad \text{(节点 $ j $ 对邻居状态的加权融合)}
   

   • 节点 i 综合邻居信息生成自身观测：

     
     b_i^H(k) = \sum_{j=1}^N a_{ji} b_{j,\text{local}}(k) + r_k
     

   • 注：r_k 为通信噪声，反映信息传输误差（如延时、丢包）。

---

Step 5: 观测更新（测量更新）

1. 观测 Sigma 点：

   
   \mathcal{Z}_{i,j,k|k-1} = h(\mathcal{X}_{i,j,k|k-1}^*) + r_k \quad (j=0,\dots,2n)
   

   • $h(\cdot)$ 为非线性观测函数

2. 计算观测统计量：

   
   \hat{z}_{i,k|k-1} = \sum_{j=0}^{2n} W_j^{(m)} \mathcal{Z}_{i,j,k|k-1}
   

   
   P_{i,zz} = \sum_{j=0}^{2n} W_j^{(c)} \left( \mathcal{Z}_{i,j,k|k-1} - \hat{z}_{i,k|k-1} \right) \left( \mathcal{Z}_{i,j,k|k-1} - \hat{z}_{i,k|k-1} \right)^T + R_k
   

   
   P_{i,xz} = \sum_{j=0}^{2n} W_j^{(c)} \left( \mathcal{X}_{i,j,k|k-1}^* - \hat{b}_{i,k|k-1} \right) \left( \mathcal{Z}_{i,j,k|k-1} - \hat{z}_{i,k|k-1} \right)^T
   

   • $R_k$ 为观测噪声协方差

3. 计算卡尔曼增益并更新状态：

   
   K_{i,k} = P_{i,xz} P_{i,zz}^{-1}
   

   
   \hat{b}_{i,k|k} = \hat{b}_{i,k|k-1} + K_{i,k} \left( b_i^H(k) - \hat{z}_{i,k|k-1} \right)
   

   
   P_{i,k|k} = P_{i,k|k-1} - K_{i,k} P_{i,zz} K_{i,k}^T
   

---

Step 6: 全局一致性计算

1. 瑞利商计算：

   • 所有节点通过一致性协议交换 $\hat{b}_{i,k|k}$，计算全局状态：

     
     y(k) = \frac{\sum_{i=1}^N x_i(k) \hat{b}_{i,k|k}}{\sum_{i=1}^N x_i^2(k)}
     

2. 正交化：

   • 更新本地状态分量(相当于幂迭代$x=Ax$再归一化)：

     
     x_i(k+1) = \frac{\hat{b}_{i,k|k}}{\|\hat{b}(k)\|_2}
     

---

Step 7: 收敛判断

• 若 y(k) 收敛，输出 $\sigma = \sqrt{y(k)}$；否则返回 Step 2。

稳态下动态特征参数的估算

稳态下，网络拓扑变化趋于平稳，奇异值的理论曲线不再随时间变化（实际值因噪声围绕理论值波动）。此时采用集中式多观测值卡尔曼滤波

多观测值滤波算法

• 核心思想：利用相邻奇异值的有序性约束（$\sigma_{n-1} \leq \sigma_n \leq \sigma_{n+1}$），构造双观测值作为上下界，限制估计范围。

• 观测值生成：
  对第$n$大奇异值$\sigma_n$，其观测值$y_n$由相邻奇异值线性组合：

  
  y_n = C_1 \sigma_{n-1} + C_2 \sigma_{n+1}
  

  • 系数$C_1, C_2$：根据$\sigma_{n-1}$和$\sigma_{n+1}$的权重动态调整（如距离比例）。
  • 物理意义：将$\sigma_{n-1}$和$\sigma_{n+1}$作为$\sigma_n$的下界和上界，避免单观测值因噪声导致的估计偏离。

疑问：

第三章的目的是什么？先分解再重构的意义在？

状态转移函数和观测函数怎么来？UKF每次预测单奇异值，如何同时预测K个呢？

卡尔曼滤波 观测值怎么来？是否需要拟合历史数据生成观测值？还是根据第三章分布式幂迭代求真实的特征值？