DOI: 10.12677/pm.2025.153102, PDF, HTML, XML,  被引量    国家自然科学基金支持
作者: 杜 雨：上海理工大学理学院，上海；张志祥, 丁德锐：上海理工大学光电信息与计算机工程学院，上海
关键词: 分布式滤波；最大相关熵；状态分解；状态饱和；S-t核函数；隐私保护；Distributed Filtering； Maximum Correntropy； State Decomposition； State Saturation； S-t Kernel Function； Privacy Protection

文章引用：杜雨, 张志祥, 丁德锐. 一类状态饱和非线性系统的隐私保护分布式最大相关熵滤波：基于状态分解的方法[J]. 理论数学, 2025, 15(3): 273-287. https://doi.org/10.12677/pm.2025.153102

1. 引言

状态估计是导航[1]、目标追踪[2]、系统识别、最优控制等[3]研究领域的重要关注之一，是获取系统动态行为的重要理论与关键技术。对于线性系统，卡尔曼滤波器(KF)是利用最小二乘思想提供状态估计问题最优递归解的线性滤波器。但现实中大多数系统均为非线性系统，从而使得KF的应用具有较大的局限性。为了解决非线性条件下的状态估计问题，许多改进的卡尔曼滤波器相继涌现，例扩展卡尔曼滤波器(EKF) [4]、无迹卡尔曼滤波器(UKF) [5]、正交卡尔曼滤波器(QKF) [6]等。

上述滤波器大都依据最小均方误差(MMSE)原则所设计的，当系统遭受的噪声符合高斯分布时，它们的估计精度较高。然而，由于MMSE原则只捕捉到了噪声的二阶统计特征，当系统遭受的噪声为非高斯噪声时，此类算法对异常值表现较为敏感从而使得算法的估计性能大大降低[7]。为了解决这一问题，学者们提出了一些基于非MMSE原则的滤波器。例如，基于Huber的卡尔曼滤波器(HKF) [8]通过采用Huber函数将残差有界化从而抑制了异常值对于状态估计的影响。最大相关熵卡尔曼滤波器(MCKF) [9]是基于信息理论学习(ITL) [10]中的最大相关熵准则(MCC)提出的。实验结果表明，其在复杂非高斯噪声干扰下展示了良好的性能，优于传统的滤波方法。

需要指出的是基于ITL的滤波器性能在很大程度上依赖于核函数及其选择的参数。目前，大多数此类滤波器仍采用高斯核函数。然而，在一些实际应用场景例如敏捷目标跟踪[11]、多路径估计等，系统通常遭受到重尾非高斯噪声的干扰。对于重尾噪声，基于高斯核的滤波器难以选择适当的核带宽而往往无法充分捕捉其重尾特性；此外，在执行过程中奇异矩阵的出现通常会导致算法的崩溃，从而可能导致估计精度的显著下降。因此，高斯核函数作为代价函数并非最优选择，探索其他类型的核函数以克服高斯核函数固有的不足具有重要意义。文献[12]中提出了一种基于Cauchy核函数的最大相关熵卡尔曼滤波器(CKKF)。与使用高斯核函数的MCKF相比，CKKF在处理多维非高斯噪声时表现更优。同时，当选择不同的Cauchy核带宽时，CKKF的性能始终保持稳定。最近，基于S-t分布提出了一种新的核函数：S-t核函数[13]。与高斯核函数不同，S-t核函数具有两个参数，可以分别控制核函数的形状和带宽，从而其更有效地捕捉噪声的重尾特性。在适当选择核函数参数的情况下，当噪声为重尾非高斯分布时，STKKF的状态估计精度优于传统算法。

另一方面，分布式滤波可以在多个节点上并行处理信息，其相较于集中式滤波具有显著的优点，包括提高系统的鲁棒性、降低通信负担、提升可扩展性等。但目前部署的网络通常是开放的。尽管已经采取了大量保护措施，但数据泄露和网络攻击等问题仍然无法避免[14]。这种情况往往会导致巨大的经济损失、负面的社会影响，甚至对国家安全构成威胁。近年来，研究者提出了几种隐私保护方法，如同态加密[15]、水印[16]和注入扰动策略等。状态饱和[17]是实际系统中普遍存在的一种现象，通常由系统状态的物理限制所引起。如果未能妥善处理状态饱和，滤波器的性能可能会迅速下降。综上所述，本文针对一类状态和非线性系统，旨在设计一种基于状态分解和S-t核函数的隐私保护分布式最大相关熵滤波算法，以增强其在非高斯环境下的鲁棒性以及数据分享的隐私性。研究面临的主要挑战包括：1) 如何选择适当的隐私保护策略设计状态饱和下的熵滤波器；2) 如何分析滤波算法的隐私性。本文系统地解决了上述挑战，给出了变量分解下的方差上界的迭代计算规则，并结合传统的定点迭代技术设计了一种新型的分布式熵滤波算法，主要贡献如下：

1) 基于分解的状态，设计了一种针对状态饱和非线性系统的分布式滤波器；

2) 基于S-t核函数的损失函数和定点迭代规则，给出了私有估计与公共估计的误差协方差矩阵的上界及其分布式最大相关熵滤波器增益。

2. 问题描述

本文考虑的传感器网络由n个传感器节点组成，节点之间的通信描述为有向图 $\mathcal{G}=\left(\mathcal{N},ℰ,\mathscr{H}\right)$ ，其中 $\mathcal{N}=\left\{1,2,\cdots ,n\right\}$ 为节点集， $ℰ\subseteq \mathcal{N}\times \mathcal{N}$ 代表边连接集， $\mathscr{H}={\left[a_{ij}\right]}_{n\times n}$ 表示传感器节点间的加权邻接矩阵。对于所有节点 $i,j\in \mathcal{N}$ ，如果 $\left(i,j\right)\in ℰ$ ，则有 $a_{ij}>0$ ，即节点i可以接收到来自节点j的信息；否则 $a_{ij}=0$ ，即节点i不能接收到来自节点j的信息。显而易见，矩阵元素 $a_{ij}$ 反应了节点i与节点j之间连接的紧密程度。最后，记第i个传感器的邻居传感器集为 ${\mathcal{N}}_i=\left\{j:\left(i,j\right)\in ℰ\right\}$ ，邻居的个数为 $d_i$ 。

考虑一类带有状态饱和的离散随机线性系统，其动力学行为可描述为

$\{\begin{array}{l}{x}_{k+1}=\sigma \left(A{x}_k+{\zeta }_{k}G{x}_k\right)+{\omega }_k\\ y_{i,k+1}=C_i{x}_{k+1}+{\nu }_{i,k+1}\end{array}$ (1)

式中 $x_{k+1}\in {ℝ}^{n_x}$ 和 $y_{i,k+1}\in {ℝ}^{n_y}$ 分别表示待估计的状态变量和传感器i的测量输出； ${\omega }_k\in {ℝ}^{n_{\omega }}$ 和 ${\nu }_{i,k}\in {ℝ}^{n_{\nu }}$ 为互不相关的非高斯噪声，它们的均值均为0且对应的协方差矩阵分别为 $\mathbb{E}\left[{\omega }_k{\omega }_k^{\text{T}}\right]=W_k$ 和 $\mathbb{E}\left[{\nu }_{i,k}{\nu }_{i,k}^{\text{T}}\right]=V_{i,k}$ ；随机变量 ${\zeta }_k$ 是均值为0方差为1的噪声序列，描述了系统遭受的乘性噪声。进而，考虑到实际工程系统通常运行在一个有限的范围内，即系统遭受状态保护，故引入非线性映射 $\sigma (\cdot ):{ℝ}^n\to {ℝ}^n$ 是一描述受状态饱和影响。具体地，对于向量 $x\triangleq {\left[x_1{x}_2\cdots x_n\right]}^{\text{T}}\in {ℝ}^n$ ，定义 $\sigma \left(x\right)={\left[{\sigma }_1\left(x_1\right){\sigma }_2\left(x_2\right)\cdots {\sigma }_n\left(x_n\right)\right]}^{\text{T}}$ ，其中 ${\sigma }_i\left(x_i\right)=\mathrm{sign}\left(x_i\right)\min \left\{{\gamma }_i,\left|x_i\right|\right\}$ 。这里， ${\gamma }_i>0$ 代表饱和水平。A，G和 $C_i$ 分别代表系统矩阵和测量矩阵。

为了监控上述目标的动态行为，理想的滤波器可设计为：

$\{\begin{array}{l}{\hat{x}}_{i,k+1}^{-}=\sigma \left(A{\hat{x}}_{i,k}^{+}+\epsilon {\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}\left(A{\hat{x}}_{i,k}^{+}-A{\hat{x}}_{j,k}^{+}\right)}\right)\\ {\hat{x}}_{i,k+1}^{+}={\hat{x}}_{i,k+1}^{-}+K_{i,k+1}\left(y_{i,k+1}-C_i{\hat{x}}_{i,k+1}^{-}\right)\end{array}$ (2)

其中 ${\hat{x}}_{i,k+1}^{-}$ 和 ${\hat{x}}_{i,k+1}^{+}$ 分别代表 $k+1$ 时刻第i个传感器的一步状态预测和一步估计， ${\hat{x}}_{j,k}^{+}$ 代表来自于邻居的理想一步估计。正标量 $\epsilon \in \left(0,1/\beta \right)$ 代表一致性增益，其中 $\beta =\max \left\{{\lambda }_i\right\}$ ， ${\lambda }_i$ 是通信拓扑图的拉普拉斯矩阵的特征值。矩阵 $K_{i,k+1}$ 为待设计的滤波器增益。不失一般性，初始状态假设满足 ${\hat{x}}_{i,k+1}^{-}={\hat{x}}_{i,0}^{+}=\mathbb{E}\left\{x_0\right\}={\mu }_0\left(i\in \mathcal{N}\right)$ 。

由(2)不难发现，分布式滤波器的实现依赖于邻居信息的交互。然而，开放的网络可能导致信息的泄露，从而被外部恶意攻击者利用，造成不必要的损失。为了保证数据传输的隐私性，对传输的估计值 ${\hat{x}}_{i,k}^{+}$ 进行状态分解：

${\hat{x}}_{i,k}^{+}={\alpha }_{i,k}^{}{\hat{x}}_{i,k}^{+}+{\beta }_{i,k}{\hat{x}}_{i,k}^{+}={\hat{x}}_{i,k}^{\alpha +}+{\hat{x}}_{i,k}^{\beta +}$ (3)

其中 ${\alpha }_{i,k}=1-{\beta }_{i,k}$ ， ${\beta }_{i,k}$ 是滤波器i的动态私钥，满足提前确定的变化规则，该规则与其邻居共享，并不被窃听所获取。需要指出的是，滤波器i的邻居j将根据该动态私钥传输信息对应的估计值 ${\hat{x}}_{ij,k}^{\beta +}={\beta }_{ik}{\hat{x}}_{j,k}^{+}$ 。根据这一规则，滤波器i仅能够收到邻居分解后的信息。为此，实际设计的滤波器为

$\{\begin{array}{l}{\hat{x}}_{i,k+1}^{-}=\sigma \left(A{\hat{x}}_{i,k}^{+}+\epsilon {\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}\left({\beta }_{ik}A{\hat{x}}_{i,k}^{+}-A{\hat{x}}_{ij,k}^{\beta +}\right)}\right)\\ {\hat{x}}_{i,k+1}^{+}={\hat{x}}_{i,k+1}^{-}+K_{i,k+1}\left(y_{i,k+1}-C_i{\hat{x}}_{i,k+1}^{-}\right)\end{array}$ (4)

其中 ${\hat{x}}_{ij,k}^{\beta +}$ 为外部窃听者可以听到的、邻居之间实际交换的数据。

对于外部窃听者来说，它可以使用网络中传输的部分数据来构建自己的状态估计，其迭代规则如下所示：

$\{\begin{array}{l}{\hat{x}}_{i,k+1,e}^{-}=\sigma \left(A{\hat{x}}_{i,k,e}^{+}+\epsilon {\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}}\left(A{\hat{x}}_{i,k,e}^{+}-A{\hat{x}}_{ij,k}^{\beta +}\right)\right)\\ {\hat{x}}_{i,k+1,e}^{+}={\hat{x}}_{i,k+1,e}^{-}+K_{i,k+1}\left(y_{i,k+1}-C_i{\hat{x}}_{i,k+1,e}^{-}\right)\end{array}$ (5)

其中 ${\hat{x}}_{i,k+1,e}^{+}$ 代表外部窃听者的估计值。

注1：在上式中， $\epsilon {\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}}\left(A{\hat{x}}_{i,k,e}^{+}-A{\hat{x}}_{ij,k}^{\beta +}\right)$ 表明了窃听者是利用从其他节点获取的信息来调整其预测值。

针对系统的饱和约束，有如下引理。

引理1 [13]：假设函数 ${\sigma }_i(\cdot )$ 如下式所定义，则有：

${\sigma }_i\left(u\right)-{\sigma }_i\left(v\right)={\lambda }_i\left(u-v\right),\forall u,v\in ℝ$

其中 ${\lambda }_i\in \left[0,1\right]$ 。

记 $e_{i,k}^{+}=x_k-{\hat{x}}_{i,k}^{+}$ ， $e_{i,k}^{-}=x_k-{\hat{x}}_{i,k}^{-}$ 。考虑到系统的饱和特性，由引理1一定存在一个满足 $0\le {\Lambda }_k\le I$ 矩阵 ${\Lambda }_k$ 使得

$\begin{array}{c}{e}_{i,k+1}^{-}=x_{k+1}-{\hat{x}}_{i,k+1}^{-}\\ =\sigma \left(A{x}_k+{\zeta }_{k}G{x}_k\right)+{\omega }_k-\sigma \left(A{\hat{x}}_{i,k}^{+}+\epsilon {\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}\left({\beta }_{ik}A{\hat{x}}_{i,k}^{+}-A{\hat{x}}_{ij,k}^{\beta +}\right)}\right)\\ ={\Lambda }_k\left(A{e}_{i,k}^{+}-\epsilon {\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}}{\beta }_{ik}A\left(e_{j,k}^{+}-e_{i,k}^{+}\right)\right)+{\gamma }_k{\Lambda }_{k}G{x}_k+{\omega }_k\end{array}$

和

$\begin{array}{c}{e}_{i,k+1}^{+}=x_{k+1}-{\hat{x}}_{i,k+1}^{+}\\ =x_{k+1}-{\hat{x}}_{i,k+1}^{-}-K_{i,k+1}\left(y_{i,k+1}-C_i{\hat{x}}_{i,k+1}^{-}\right)\\ =\left(I_n-K_{i,k+1}{C}_i\right)e_{i,k+1}^{-}-K_{i,k+1}{v}_{i,k+1}.\end{array}$

在高斯噪声下，卡尔曼滤波算法实现了最小均方误差准下的最优估计。然而，当系统中的 ${\omega }_{i,k}$ 和 ${\nu }_{i,k}$ 均为服从非高斯分布时，传统的卡尔曼滤波算法并不能得到最优估计。特别地，当噪声呈现重尾分布式或存在野值的时候，可能造成滤波性能的衰减。针对这一问题，本文采用了基于高斯核函数的最大熵准则提高滤波器的滤波性能。为此，系统(1)中的测量可写为

$\left[\begin{array}{c}{\hat{x}}_{i,k}^{-}\\ y_{i,k}\end{array}\right]=\left[\begin{array}{c}{I}_{n_x}\\ C_i\end{array}\right]x_k+{\xi }_{i,k}$ (6)

其中

${\xi }_{i,k}=\left[\begin{array}{c}{\hat{x}}_{i,k}^{-}-x_k\\ {\nu }_{i,k}\end{array}\right]$

对该向量求方差，可得

$\mathbb{E}\left[{\xi }_{i,k}{\xi }_{i,k}^{\text{T}}\right]=\left[\begin{array}{cc}{P}_{i,k}^{-}& 0\\ 0& V_{i,k}\end{array}\right]=\left[\begin{array}{cc}{\Omega }_{P_{i,k}}{\Omega }_{P_{i,k}}^{\text{T}}& 0\\ 0& {\Omega }_{V_{i,k}}{\Omega }_{V_{i,k}}^{\text{T}}\end{array}\right] ={\Omega }_{i,k}{\Omega }_{i,k}^{\text{T}}$ (7)

其中 $P_{i,k}^{-}$ 代表一步预测协方差矩阵； ${\Omega }_{i,k}$ 为 $\mathbb{E}\left[{\xi }_{i,k}{\xi }_{i,k}^{\text{T}}\right]$ 的Cholesky分解矩阵。为了对向量 ${\xi }_{i,k}$ 进行归一化处理，将式(6)两端同时乘以 ${\Omega }_{i,k}^{-1}$ ，得到

${\bar{y}}_{i,k}={\mathfrak{M}}_{i,k}{x}_k+{\eta }_{i,k}$ (8)

其中 ${\eta }_{i,k}={\Omega }_{i,k}^{-1}{\xi }_{i,k}$ ，以及

${\mathfrak{M}}_{i,k}={\Omega }_{i,k}^{-1}\left[\begin{array}{c}{I}_{n_x}\\ C_i\end{array}\right].$

对于这一分解，不难发现 $\mathbb{E}\left[{\eta }_{i,k}{\eta }_{i,k}^{\text{T}}\right]=I$ ，即将 ${\eta }_{i,k}$ 化为白噪声。此外，由于系统的分布式特性，精确的方差矩阵 $P_{i,k}^{-}$ 是难以获得的。一个可行的方法是寻找一个较小的上界加以替代。根据最大熵原则，同时为了更好地捕获噪声中的重尾特征，本文使用S-t核函数来代替高斯核函数[13]，采用的损失函数如下：

$J_{n_x+n_y}\left({\hat{x}}_{i,k}^{+}\right)=\frac{1}{n_x+n_y}\left({\displaystyle \underset{\ell =1}{\overset{n_x}{\sum }}{{\rm Y}}_{v,{\lambda }_k}}{\left[{\eta }_{x_{i,k}}\right]}_{\ell }+{\displaystyle \underset{\ell =1}{\overset{n_y}{\sum }}{{\rm Y}}_{v,{\lambda }_k}{\left[{\eta }_{y_{i,k}}\right]}_{\ell }}\right)$

其中 ${\lambda }_k$ 为带宽或者核宽度，v用于控制Student’s t核函数的形状，其他向量为

${\eta }_{x_{i,k}}={\Omega }_{P_{i,k}}^{-1}\left({\hat{x}}_{i,k}^{+}-{\hat{x}}_{i,k}^{-}\right),{\eta }_{y_{i,k}}={\Omega }_{V_{i,k}}^{-1}\left(y_{i,k}-C_i{\hat{x}}_{i,k}^{+}\right),{{\rm Y}}_{v,{\lambda }_k}\left(\zeta \right)={\left(1+\frac{\zeta }{v{\lambda }_k^2}\right)}^{-\frac{v+2}{2}}$

在该指标中， ${\Omega }_{P_{i,k}}^{}$ 和 ${\Omega }_{V_{i,k}}^{}$ 是一组权重矩阵，在实现白化处理的同时，实现了状态估计稳定和状态估计准确的权衡。

综述所述，本文的设计目标如下

1) 借助于估计误差方差，获得理想的权重矩阵 ${\Omega }_{i,k}^{}$ ；

2) 设计合适的滤波增益矩阵 $K_{i,k+1}$ 获得最优估计值 ${\hat{x}}_{i,k}^{+}$ ，使得

${\hat{x}}_{i,k}^{+}=\arg {\max }_{{\hat{x}}_{i,k}^{+}}{J}_{n_x+n_y}\left({\hat{x}}_{i,k}^{+}\right).$ (9)

基于MCC下的滤波器增益计算

本节致力于设计基于MCC原则的分布式隐私保护滤波器。为此，记

$P_{ii,k}^{-}= \mathbb{E}\left\{{\hat{e}}_{i,k}^{-}{\left({\hat{e}}_{i,k}^{-}\right)}^{\text{T}}\right\},P_{ii,k}^{+}=\mathbb{E}\left\{{\hat{e}}_{i,k}^{+}{\left({\hat{e}}_{i,k}^{+}\right)}^{\text{T}}\right\}$

代表先验协方差矩阵和后验协方差矩阵。进而， ${\hat{e}}_{i,k}^{-}$ 和 ${\hat{e}}_{i,k}^{+}$ 相应的分解可以表示为：

${\hat{e}}_{i,k}^{-}={\alpha }_{i,k}{\hat{e}}_{i,k}^{-}+{\beta }_{i,k}{\hat{e}}_{i,k}^{-}={\hat{e}}_{i,k}^{\alpha -}+{\hat{e}}_{i,k}^{\beta -},\begin{array}{c}\end{array}{\hat{e}}_{i,k}^{+}={\alpha }_{i,k}{\hat{e}}_{i,k}^{+}+{\beta }_{i,k}{\hat{e}}_{i,k}^{+}={\hat{e}}_{i,k}^{\alpha +}+{\hat{e}}_{i,k}^{\beta +}.$

结合引入的状态分解规则和引理1，两组误差动态可表示为：

$\{\begin{array}{l}{\hat{e}}_{i,k+1}^{\alpha -}={\Lambda }_{k}A{\hat{e}}_{i,k}^{\alpha +}+{\gamma }_k{\alpha }_{ik}{\Lambda }_{k}G{x}_k+{\alpha }_{ik}{w}_k\\ {\hat{e}}_{i,k+1}^{\alpha +}=\left(I_n-K_{i,k+1}{C}_i\right){\hat{e}}_{i,k+1}^{\alpha -}-{\alpha }_{ik}{K}_{i,k+1}{v}_{i,k+1}\end{array}$ (10)

和

$\{\begin{array}{l}{\hat{e}}_{i,k+1}^{\beta -}={\Lambda }_k\left(A{e}_{i,k}^{\beta +}-\epsilon {\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}}A\left(e_{j,k}^{\beta +}-e_{i,k}^{\beta +}\right)\right)+{\gamma }_k{\beta }_{ik}{\Lambda }_{k}G{x}_k+{\beta }_{ik}{\omega }_k\\ {\hat{e}}_{i,k+1}^{\beta +}=\left(I_n-K_{i,k+1}{C}_i\right){\hat{e}}_{i,k+1}^{\beta -}-{\beta }_{ik}{K}_{i,k+1}{v}_{i,k+1}\end{array}$ (11)

明显地，误差动态(10)是由私钥分解部分驱动产生的，可称之为私有误差动态；误差动态(11)是由公钥分解部分驱动生成的，可称之为公共误差动态。对比窃听者的估计策略，不难发现其窃听不准确的原因正是由于私钥分解部分的缺失。

对于状态分解后的一步预测协方差矩阵和误差协方差矩阵满足如下关系：

$P_{i,k}^{+}=P_{i,k}^{\alpha +}+P_{i,k}^{\beta +}+2{\alpha }_{i,k}{\beta }_{i,k}{P}_{i,k}^{+}={\left(1-2{\alpha }_{i,k}{\beta }_{i,k}\right)}^{-1}\left(P_{i,k}^{\alpha +}+P_{i,k}^{\beta +}\right)$ (12)

和

$P_{i,k}^{-}=P_{i,k}^{\alpha -}+P_{i,k}^{\beta -}+2{\alpha }_{i,k}{\beta }_{i,k}{P}_{i,k}^{-}={\left(1-2{\alpha }_{i,k}{\beta }_{i,k}\right)}^{-1}\left(P_{i,k}^{\alpha -}+P_{i,k}^{\beta -}\right)$ (13)

定理1 假设分解误差动态(即私有误差动态(10)和公共误差动态(11))的私有和公共估计误差协方差 $P_{i,k}^{\alpha +}$ 和 $P_{i,k}^{\beta +}$ 的初值满足 $P_{i,0}^{\alpha +}={\Xi }_{i,0}^{\alpha +}$ 和 $P_{i,0}^{\beta +}={\Xi }_{i,0}^{\beta +}$ ，先验误差协方差 $P_{i,k}^{\alpha -}$ 和 $P_{i,k}^{\beta -}$ 的初值满足 $P_{i,0}^{\alpha -}={\Xi }_{i,0}^{\alpha -}$ 和 $P_{i,0}^{\beta -}={\Xi }_{i,0}^{\beta -}$ 。分布式隐私保护滤波器(4)的先验协方差矩阵和后验协方差矩阵的上界分别为

$\{\begin{array}{l}{P}_{ii,k}^{+}\le {\Xi }_{i,k}^{+}={\left(1-2{\alpha }_{i,k}{\beta }_{i,k}\right)}^{-1}\left({\Xi }_{i,k}^{\alpha +}+{\Xi }_{i,k}^{\beta +}\right),\\ P_{ii,k}^{-}\le {\Xi }_{i,k}^{-}={\left(1-2{\alpha }_{i,k}{\beta }_{i,k}\right)}^{-1}\left({\Xi }_{i,k}^{\alpha -}+{\Xi }_{i,k}^{\beta -}\right),\end{array}$ (14)

其中 ${\Xi }_{i,k}^{\alpha +}$ 、 ${\Xi }_{i,k}^{\beta +}$ 、 ${\Xi }_{i,k}^{\alpha -}$ 和 ${\Xi }_{i,k}^{\beta -}$ 满足如下递推矩阵约束

$\{\begin{array}{l}{\Xi }_{ii,k+1}^{\alpha +}=\left(I-K_{i,k+1}{C}_{i,k+1}\right){\Xi }_{ii,k+1}^{\alpha -}{\left(I-K_{i,k+1}{C}_{i,k+1}\right)}^{\text{T}}+{\alpha }_{i,k+1}^2{K}_{i,k+1}{V}_{i,k+1}{K}_{i,k+1}^{\text{T}}\\ {\Xi }_{ii,k+1}^{\alpha -}=\text{tr}\left\{A{\Xi }_{ii,k}^{\alpha +}{A}^{\text{T}}\right\}I+\text{tr}\left\{{\alpha }_{ik}^{2}G{\Psi }_k{G}^{\text{T}}\right\}I+{\alpha }_{ik}^2{W}_k\end{array}$ (15)

和

$\{\begin{array}{l}{\Xi }_{ii,k+1}^{\beta +}=\left(I-K_{i,k+1}{C}_{i,k+1}\right){\Xi }_{ii,k+1}^{\beta -}{\left(I-K_{i,k+1}{C}_{i,k+1}\right)}^{\text{T}}+{\beta }_{i,k+1}^2{K}_{i,k+1}{V}_{i,k+1}{K}_{i,k+1}^{\text{T}}\\ {\Xi }_{ii,k+1}^{\beta -}={\beta }_{ik}^2{W}_k+\text{tr}\left\{\left(1+\pi \right)A{\Xi }_{ii,k}^{\beta +}{A}^{\text{T}}+2{\epsilon }^2{d}_i\left(1+{\pi }^{-1}\right){\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}^{2}A\left({\Xi }_{ii,k}^{\beta +}+{\Xi }_{jj,k}^{\beta +}\right)A^{\text{T}}}\right\}I+\text{tr}\left\{{\beta }_{ik}^{2}G{\Psi }_k{G}^{\text{T}}\right\}I\end{array}$ (16)

这里， ${\Psi }_k=\min \left\{\text{tr}\left\{\left(1+\rho \right){\Xi }_{jj,k}^{\beta +}+\left(1+{\rho }^{-1}\right){\hat{x}}_{j,k}^{\beta +}{\hat{x}}_{j,k}^{\beta +\text{T}}\right\},\text{tr}\left\{{\displaystyle \underset{i=1}{\overset{n_X}{\sum }}{\gamma }_i^2}I+W_k\right\}\right\}$。

证明：根据方差的定义，计算可得

(17)

同理可得

$\begin{array}{c}{P}_{ii,k+1}^{\alpha +}=\mathbb{E}\left\{{\hat{e}}_{i,k+1}^{\alpha +}{\left({\hat{e}}_{i,k+1}^{\alpha +}\right)}^{\text{T}}\right\}\\ =\left(I_n-K_{i,k+1}{C}_i\right)P_{ii,k+1}^{\alpha -}{\left(I_n-K_{i,k+1}{C}_i\right)}^{\text{T}}+{\alpha }_{ik}^2{K}_{i,k+1}{V}_{i,k+1}{K}_{i,k+1}^{\text{T}}\\ P_{ii,k+1}^{\beta +}=\mathbb{E}\left\{e_{i,k+1}^{\beta +}{\left(e_{i,k+1}^{\beta +}\right)}^{\text{T}}\right\}\\ = \left(I_n-K_{i,k+1}{C}_i\right)P_{ii,k+1}^{\beta -}{\left(I_n-K_{i,k+1}{C}_i\right)}^{\text{T}}+{\beta }_{ik}^2{K}_{i,k+1}{V}_{i,k+1}{K}_{i,k+1}^{\text{T}}\end{array}$ (18)

另一方面，让我们分析系统的二阶矩，可知

$\begin{array}{c}{X}_k=\mathbb{E}\left\{x_k{x}_k^{\text{T}}\right\}\\ =\mathbb{E}\left\{\left(\sigma \left(A{x}_k+{\zeta }_{k}G{x}_k\right)+{\omega }_k\right){\left(\sigma \left(A{x}_k+{\zeta }_{k}G{x}_k\right)+{\omega }_k\right)}^{\text{T}}\right\}\\ =\mathbb{E}\left\{\sigma \left(A{x}_k+{\zeta }_{k}G{x}_k\right){\sigma }^{\text{T}}\left(A{x}_k+{\zeta }_{k}G{x}_k\right)\right\}+W_k\\ \le \text{tr}\mathbb{E}\left\{\sigma \left(A{x}_k+{\zeta }_{k}G{x}_k\right){\sigma }^{\text{T}}\left(A{x}_k+{\zeta }_{k}G{x}_k\right)\right\}I+W_k\\ \le {\displaystyle \underset{i=1}{\overset{n_X}{\sum }}{\gamma }_i^2}I+W_k\end{array}$

又根据[18]得

$\begin{array}{c}{X}_k=\mathbb{E}\left\{x_k{x}_k^{\text{T}}\right\}\\ =\mathbb{E}\left\{{\hat{x}}_{j,k}^{\beta +}{\hat{x}}_{j,k}^{\beta +\text{T}}+{\hat{x}}_{j,k}^{\beta +\text{T}}{e}_{j,k}^{\beta +\text{T}}+e_{j,k}^{\beta +\text{T}}{\hat{x}}_{j,k}^{\beta +\text{T}}+e_{j,k}^{\beta +}{e}_{j,k}^{\beta +\text{T}}\right\}\\ \le \left(1+\rho \right)P_{jj,k}^{\beta +}+\left(1+{\rho }^{-1}\right){\hat{x}}_{j,k}^{\beta +}{\hat{x}}_{j,k}^{\beta +\text{T}}\end{array}$

因此通过结合两式我们得到：

$tr\left(X_k\right)\le \min \left\{\text{tr}\left\{\left(1+\rho \right)P_{jj,k}^{\beta +}+\left(1+{\rho }^{-1}\right){\hat{x}}_{j,k}^{\beta +}{\hat{x}}_{j,k}^{\beta +\text{T}}\right\},\text{tr}\left\{{\displaystyle \underset{i=1}{\overset{n_X}{\sum }}{\gamma }_i^2}I+W_k\right\}\right\}\triangleq {\Psi }_k$

接下来，让我们处理(17)中的交叉项，不难有

(19)

综上，可得 $D_{ii,k+1}^{\beta -}$ 的上界为：

$\begin{array}{c}{P}_{ii,k+1}^{\beta -}\le {\beta }_{ik}^2{W}_k+\text{tr}\left\{\left(1+\pi \right)A{P}_{ii,k}^{\beta +}{A}^{\text{T}}+2{\epsilon }^2{d}_i\left(1+{\pi }^{-1}\right){\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}^{2}A\left(P_{ii,k}^{\beta +}+P_{jj,k}^{\beta +}\right)A^{\text{T}}}\right\}I\\ +\text{tr}\left\{{\beta }_{ik}^{2}G{\Psi }_k{G}^{\text{T}}\right\}I\end{array}$ (20)

最后根据[19]可得到分解后的一步预测协方差矩阵和误差协方差矩阵。 证毕■

定理2 假设分解误差动态(即私有误差动态(10)和公共误差动态(11))的私有和公共估计误差协方差 $P_{i,k}^{\alpha +}$ 和 $P_{i,k}^{\beta +}$ 的初值满足 $P_{i,0}^{\alpha +}={\Xi }_{i,0}^{\alpha +}$ 和 $P_{i,0}^{\beta +}={\Xi }_{i,0}^{\beta +}$ ，先验误差协方差 $P_{i,k}^{\alpha -}$ 和 $P_{i,k}^{\beta -}$ 的初值满足 $P_{i,0}^{\alpha -}={\Xi }_{i,0}^{\alpha -}$ 和 $P_{i,0}^{\beta -}={\Xi }_{i,0}^{\beta -}$ 。如果滤波增益取

$K_{i,k}={\left({\Omega }_{{\Xi }_{i,k}}{\kappa }_{x_{i,k}}^{-1}{\Omega }_{{\Xi }_{i,k}}^{\text{T}}+C_i^{\text{T}}{\Omega }_{V_{i,k}}{C}_{P_{i,k}}{\kappa }_{y_{i,k}}^{-1}{\Omega }_{V_{i,k}}^{\text{T}}{C}_i\right)}^{-1}{C}_i^{\text{T}}{\left({\Omega }_{V_{i,k}}{\kappa }_{y_{i,k}}^{-1}{\Omega }_{V_{i,k}}^{\text{T}}\right)}^{-1}$

那么，分布式隐私保护滤波器(4)的先验协方差矩阵和后验协方差矩阵的上界满足(14)，其中

$\begin{array}{c}{\kappa }_{x_{i,k}}=\mathrm{diag}\left\{\frac{v{\lambda }_k^2{{\rm Y}}_{v,{\lambda }_k}\left({\left[{\eta }_{x_{i,k}}\right]}_1\right)}{v{\lambda }_k^2+{\left[{\eta }_{x_{i,k}}\right]}_1^2},\cdots ,\frac{v{\lambda }_k^2{{\rm Y}}_{v,{\lambda }_k}\left({\left[{\eta }_{x_{i,k}}\right]}_{n_x}\right)}{v{\lambda }_k^2+{\left[{\eta }_{x_{i,k}}\right]}_{n_x}^2}\right\},\\ {\kappa }_{y_{i,k}}=\mathrm{diag}\left\{\frac{v{\lambda }_k^2{{\rm Y}}_{v,{\lambda }_k}\left({\left[{\eta }_{y_{i,k}}\right]}_1\right)}{v{\lambda }_k^2+{\left[{\eta }_{y_{i,k}}\right]}_1^2},\cdots ,\frac{v{\lambda }_k^2{{\rm Y}}_{v,{\lambda }_k}\left({\left[{\eta }_{y_{i,k}}\right]}_{n_y}\right)}{v{\lambda }_k^2+{\left[{\eta }_{y_{i,k}}\right]}_{n_y}^2}\right\},\end{array}$

证明：令

$\frac{\partial {\mathcal{J}}_{n_x+n_y}\left({\hat{x}}_{i,k}^{+}\right)}{\partial {\hat{x}}_{i,k}^{+}}=0$

可得

$\begin{array}{c}\frac{\partial {\mathcal{J}}_{n_x+n_y}\left({\hat{x}}_{i,k}^{+}\right)}{\partial {\hat{x}}_{i,k}^{+}}={\displaystyle \underset{l=1}{\overset{n_x}{\sum }}\frac{v{\lambda }_k^2}{v{\lambda }_k^2+{\eta }_{x_{i,k}}^2}{{\rm Y}}_{v,\sigma }}\left({\left[{\eta }_{x_{i,k}}\right]}_l\right) {\left[{\Omega }_{P_{i,k}}^{-1}\right]}_l^{\text{T}}{\left[{\Omega }_{P_{i,k}}^{-1}\right]}_l\left({\hat{x}}_{i,k}^{+}-{\hat{x}}_{i,k}^{-}\right)\\ -{\displaystyle \underset{l=1}{\overset{n_y}{\sum }}\frac{v{\lambda }_k^2}{v{\lambda }_k^2+{\eta }_{y_{i,k}}^2}{{\rm Y}}_{v,\sigma }}\left({\left[{\eta }_{y_{i,k}}\right]}_l\right){\left[{\Omega }_{V_{i,k}}^{-1}\right]}_l^{\text{T}}{\left[{\Omega }_{V_{i,k}}^{-1}\right]}_l\left(y_{i,k}-C_i{\hat{x}}_{i,k}^{+}\right)\\ =0.\end{array}$

上式可以化简为：

${\Omega }_{P_{i,k}}^{-\text{T}}{\kappa }_{x_{i,k}}{\Omega }_{P_{i,k}}^{-1}\left({\hat{x}}_{i,k}^{+}-{\hat{x}}_{i,k}^{-}\right)-{\Omega }_{V_{i,k}}^{-\text{T}}{\kappa }_{y_{i,k}}{\Omega }_{V_{i,k}}^{-}\left(y_{i,k}-C_i{\hat{x}}_{i,k}^{+}\right)=0$

移项可得 ${\hat{x}}_{i,k}^{+}={\hat{x}}_{i,k}^{-}+K_{i,k}\left(y_{i,k}-C_i{\hat{x}}_{i,k}^{+}\right)$ ，其中

$K_{i,k}={\left({\Omega }_{P_{i,k}}{\kappa }_{x_{i,k}}^{-1}{\Omega }_{P_{i,k}}^{\text{T}}+C_i^{\text{T}}{\Omega }_{V_{i,k}}{C}_{P_{i,k}}{\kappa }_{y_{i,k}}^{-1}{\Omega }_{V_{i,k}}^{\text{T}}{C}_i\right)}^{-1}{C}_i^{\text{T}}{\left({\Omega }_{V_{i,k}}{\kappa }_{y_{i,k}}^{-1}{\Omega }_{V_{i,k}}^{\text{T}}\right)}^{-1}$

即对应的权重矩阵取作 ${\Omega }_{Pi,k}^{}={\Omega }_{\Xi i,k}^{}$ 。因此，最优估计值 ${\hat{x}}_{i,k}^{+}$ 可以直接由下式得到

${\hat{x}}_{i,k}^{+}={\hat{x}}_{i,k}^{-}+K_{i,k}\left(y_{i,k}-B_i{\hat{x}}_{i,k}^{-}\right).$

由此可得该定理成立。

最后，本文设计的算法的具体实现步骤如下表1所示。

Table 1. Algorithm pseudocode

表1. 算法伪代码

步骤1	参数初始化。给定初始值 $x_0={\hat{x}}_{i,0}^{+},{\Xi }_{i,0}^{+}={\Xi }_{i,0}^{-}$ 令 $k=1$
步骤2	收集 $y_{i,k}$ ；通过(5)和(14)计算 ${\hat{x}}_{i,k}^{-}$ 和 ${\Xi }_{ii,k}^{\alpha -}$ ， ${\Xi }_{ii,k}^{\beta -}$ ， ${\Xi }_{ii,k}^{-}$ ；通过对 ${\Xi }_{ii,k}^{-}$ 和 $V_{ii,k}$ 的Cholesky分解计算 $C_{{\Xi }_{ii,k}}$ 和 $C_{V_{i,k}}$ ；
步骤3	记 ${\left({\hat{x}}_{i,k}^{+}\right)}_t$ 为固定点迭代步骤t处的估计值，并令 ${\left({\hat{x}}_{i,k}^{+}\right)}_0={\hat{x}}_{i,k}^{-}$ 作为初始迭代值；
步骤4	令 $t=0$ ；当 $\varpi >{\varsigma }^{*}$ 或 $t<t^{*}$ 时；更新迭代步骤 $t=t+1$ ； 1) 计算向量 ${\left({\tilde{\eta }}_{x_{i,k}}\right)}_t={\Omega }_{{\Xi }_{ii,k}}^{-1}\left({\left({\hat{x}}_{i,k}^{+}\right)}_{t-1}-{\hat{x}}_{i,k}^{-}\right)$ ，以及 ${\left({\tilde{\eta }}_{y_{i,k}}\right)}_t={\Omega }_{V_{i,k}}\left(y_{i,k}-C_i{\left({\hat{x}}_{i,k}^{+}\right)}_{t-1}\right)$ ； 2) 计算矩阵 ${\left({\tilde{\kappa }}_{x_{i,k}}\right)}_t$ 和 ${\left({\tilde{\kappa }}_{y_{i,k}}\right)}_t$ ； 3) 计算矩阵 ${\left({\tilde{\Xi }}_{ii,k}^{-}\right)}_t={\Omega }_{{\Xi }_{i,k}}{\left({\tilde{\kappa }}_{x_{i,k}}\right)}_t^{-1}{\Omega }_{{\Xi }_{i,k}}^{\text{T}}$， ${\left({\tilde{V}}_{ii,k}\right)}_t={\Omega }_{V_{i,k}}{\left({\tilde{\kappa }}_{y_{i,k}}\right)}_t^{-1}{\Omega }_{V_{i,k}}^{\text{T}}$； 4) 计算滤波增益 $K_{i,k}$ ； 5) 使用以下公式更新估计值 ${\left({\hat{x}}_{i,k}^{+}\right)}_t={\hat{x}}_{i,k}^{-}+K_{i,k}\left(y_{i,k}-B_i{\hat{x}}_{i,k}^{-}\right)$ ； 6) 计算迭代精度 $\varpi =\frac{\Vert {\left({\hat{x}}_{i,k}^{+}\right)}_t-{\left({\hat{x}}_{i,k}^{+}\right)}_{t-1}\Vert }{\Vert {\left({\hat{x}}_{i,k}^{+}\right)}_{t-1}\Vert }$ ；
步骤5	记 ${\hat{x}}_{i,k}^{+}={\left({\hat{x}}_{i,k}^{+}\right)}_t$ ；通过(15)更新 ${\Xi }_{ii,k}^{\alpha +}$ ， ${\Xi }_{ii,k}^{\beta +}$ 和 ${\Xi }_{ii,k}^{+}$ ，结束算法。

3. 隐私性分析

由于状态分解策略对于外部窃听者来说是未知的，因此利用类似于定理1的方法，可得窃听端的先验误差协方差矩阵和后验协方差矩阵为：

(21)

对应的上界为：

$\{\begin{array}{l}{\Xi }_{ii,k+1,e}^{-}=W_k+\text{tr}\left\{\left(1+\pi \right)A{\Xi }_{ii,k,e}^{+}{A}^{\text{T}}+2{\epsilon }^2{d}_i\left(1+{\pi }^{-1}\right){\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}^{2}A\left({\Xi }_{ii,k,e}^{+}+{\Xi }_{jj,k,e}^{+}\right)}{A}^{\text{T}}\right\}I+\text{tr}\left\{G{\Psi }_{k,e}{G}^{\text{T}}\right\}I\hfill \\ {\Xi }_{ii,k+1,e}^{+}=\left(I-K_{i,k+1,e}{C}_{i,k+1}\right){\Xi }_{ii,k+1,e}^{-}{\left(I-K_{i,k+1}{C}_{i,k+1}\right)}^{\text{T}}+K_{i,k+1}{V}_{i,k+1}{K}_{i,k+1}^{\text{T}}.\hfill \end{array}$ (22)

接下来，让我们比较窃听者的估计误差协方差与正常滤波器的误差协方差之间的关系。具体地，由 (14)、(15)和(16)可得

$\begin{array}{c}\frac{{\Xi }_{i,k+1}^{\alpha +}+{\Xi }_{i,k+1}^{\beta +}}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}=\frac{1-2{\alpha }_{i,k}{\beta }_{i,k}}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}\left(I-K_{i,k+1}{C}_{i,k+1}\right)\frac{{\Xi }_{ii,k+1}^{\alpha -}+{\Xi }_{ii,k+1}^{\beta -}}{1-2{\alpha }_{i,k}{\beta }_{i,k}}{\left(I-K_{i,k+1}{C}_{i,k+1}\right)}^{\text{T}}\\ +\frac{{\alpha }_{ik}^2+{\beta }_{ik}^2}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}{K}_{i,k+1}{V}_{i,k+1}{K}_{i,k+1}^{\text{T}},\end{array}$

以及

$\begin{array}{c}\frac{{\Xi }_{i,k+1}^{\alpha -}+{\Xi }_{i,k+1}^{\beta -}}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}=\frac{{\alpha }_{ik}^2+{\beta }_{ik}^2}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}{W}_k+\frac{{\alpha }_{ik}^2+{\beta }_{ik}^2}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}\text{tr}\left\{G{\Psi }_k{G}^{\text{T}}\right\}I\\ +\frac{1-2{\alpha }_{i,k}{\beta }_{i,k}}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}\text{tr}\{\left(1+\pi \right)A\frac{{\Xi }_{ii,k}^{\alpha +}+{\Xi }_{ii,k}^{\beta +}}{1-2{\alpha }_{i,k}{\beta }_{i,k}}{A}^{\text{T}}\\ +2{\epsilon }^2{d}_i\left(1+{\pi }^{-1}\right){\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}^{2}A\frac{{\Xi }_{ii,k}^{\alpha +}+{\Xi }_{jj,k}^{\alpha +}+{\Xi }_{ii,k}^{\beta +}+{\Xi }_{jj,k}^{\beta +}}{1-2{\alpha }_{i,k}{\beta }_{i,k}}{A}^{\text{T}}}\}I\\ -\frac{1-2{\alpha }_{i,k}{\beta }_{i,k}}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}\text{tr}\{\frac{\pi A{\Xi }_{ii,k}^{\alpha +}{A}^{\text{T}}}{1-2{\alpha }_{i,k}{\beta }_{i,k}}+2{\epsilon }^2{d}_i{\pi }^{-1}{\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}^{2}A\frac{{\Xi }_{ii,k}^{\alpha +}+{\Xi }_{jj,k}^{\alpha +}}{1-2{\alpha }_{i,k}{\beta }_{i,k}}{A}^{\text{T}}}\}I.\end{array}$

鉴于状态分解的复杂性，文本只讨论一类特殊的分解规则的隐私性，即 ${\alpha }_{i,k+1}^{}={\alpha }_{ik}^{}$ 且取

${\Psi }_k={\Psi }_{k,e}=\text{tr}\left\{{\displaystyle \underset{i=1}{\overset{n_X}{\sum }}{\gamma }_i^2}I+W_k\right\}$ 。在这种情况下，我们不难有

$\{\begin{array}{c}{\Xi }_{i,k+1}^{-}=W_k+\text{tr}\left\{G{\Psi }_k{G}^{\text{T}}\right\}I+\text{tr}\left\{\left(1+\pi \right)A{\Xi }_{ii,k}^{+}{A}^{\text{T}}+2{\epsilon }^2{d}_i\left(1+{\pi }^{-1}\right){\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}^{2}A\left({\Xi }_{ii,k}^{+}+{\Xi }_{jj,k}^{+}\right)A^{\text{T}}}\right\}I\\ -\text{tr}\{\frac{\pi A{\Xi }_{ii,k}^{\alpha +}{A}^{\text{T}}}{1-2{\alpha }_{i,k}{\beta }_{i,k}}+2{\epsilon }^2{d}_i{\pi }^{-1}{\displaystyle \underset{j\in {\mathcal{N}}_i}{\sum }{a}_{ij}^{2}A\frac{{\Xi }_{ii,k}^{\alpha +}+{\Xi }_{jj,k}^{\alpha +}}{1-2{\alpha }_{i,k}{\beta }_{i,k}}}{A}^T\}I,\\ {\Xi }_{i,k+1}^{+}=\left(I-K_{i,k+1}{C}_{i,k+1}\right){\Xi }_{ii,k+1}^{-}{\left(I-K_{i,k+1}{C}_{i,k+1}\right)}^{\text{T}}+K_{i,k+1}{V}_{i,k+1}{K}_{i,k+1}^{\text{T}}.\end{array}$ (23)

如果 ${\Xi }_{ii,0}^{+}={\Xi }_{ii,0,e}^{+}$ ，对比(22)和(23)，不难有 ${\Xi }_{i,k+1}^{-}<{\Xi }_{ii,k+1,e}^{-}$ 。因此，窃听者估计的状态的有效性低于正常的滤波器。因此，系统的隐私性是可以满足的。对于更一般的状态分解情形，即 ${\alpha }_{i,k+1}^{}\ne {\alpha }_{ik}^{}$ 且取

${\Psi }_k={\Psi }_{k,e}=\text{tr}\left\{{\displaystyle \underset{i=1}{\overset{n_X}{\sum }}{\gamma }_i^2}I+W_k\right\}$ 时，如果分解策略使得

${\theta }_{ik}=\frac{{\alpha }_{ik}^2+{\beta }_{ik}^2}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}=\frac{1-2{\alpha }_{i,k}{\beta }_{i,k}}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}\ge 1$ ,

那么

$\begin{array}{c}\frac{{\alpha }_{ik}^2+{\beta }_{ik}^2}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}{W}_k+\frac{{\alpha }_{ik}^2+{\beta }_{ik}^2}{1-2{\alpha }_{i,k+1}{\beta }_{i,k+1}}\text{tr}\left\{G{\Psi }_k{G}^{\text{T}}\right\}I\\ =W_k+\text{tr}\left\{G{\Psi }_k{G}^{\text{T}}\right\}I+\left({\theta }_{ik}-1\right)\left(W_k+\text{tr}\left\{G{\Psi }_k{G}^{\text{T}}\right\}I\right)\end{array}$

所采用的状态分解算法可能会降低系统的隐私性。反之，系统的隐私性的隐私性将被增强。因此，理想的分解策略应当满足 ${\alpha }_{i,k}{\beta }_{i,k}$ 是不增的。最后，不难发现矩阵 ${\Xi }_{ii,k+1}^{+}$ 反映了系统的估计性能，矩阵 ${\Xi }_{ii,k+1,e}^{+}$ 反映了系统的安全性能。因此，可以引入调节参数 $\epsilon \in \left(\text{0,1}\right)$ 构建了一个综合估计和安全性能的复合指标 ${\Xi }_{ii,k+1}^{+}+\epsilon {\Xi }_{ii,k+1,e}^{+}$ 。这时， $\epsilon$ 可以被称作是风险水平。当 $\epsilon$ 接近于零时表明系统对隐私性的关注较小，这将会造成较大的安全风险。

4. 仿真实验

在本节中，通过仿真算例来验证所提出的具有状态分解的安全分布式估计算法的有效性。所采用的微电网包含两个分布式发电机组，其物理参数与文献[19]中的一致。此外，系统中部署了4个互联传感器，其中加权邻接矩阵为

$\mathscr{H}=\left[\begin{array}{cccc}0& 0.2& 0& 0.8\\ 0.2& 0& 0& 0.3\\ 0& 0& 0& 0.5\\ 0.8& 0.3& 0.5& 0\end{array}\right]$

在不失一般性的情况下，所考虑的电力系统的初始状态平均值和设计的估计器分别选为 $x_0={\hat{x}}_{i,0}^{+}={\left[1;3;1.5;1.6;2.2;1.3;1;1;1;1;1;1\right]}^{\text{T}}$，并且取 $P_{i,0}^{\alpha +}=P_{i,0}^{\beta +}=P_{i,0}^{+}=0.1I_{12}$ ， $G=0.1\ast randn\left(48,12\right)$ 为Matlab随机生成的矩阵，另外共识系数 $\epsilon =0.1$ ，S-t核函数中两个参数分别取做 ${\lambda }_k=0.5$ ， $v=4$ ，高斯核函数的参数取为 ${\lambda }_k=0.5$ ，其它参数设置为 $\pi =1.2$ ， $\rho =0.35$ ， ${\varsigma }^{*}={10}^{-5}$ ， ${\gamma }_i=0.4$ 和 $t^{*}={10}^3$ 。通过使用Matlab软件生成变量分解参数 ${\alpha }_{ik}$ 服从均匀分布 $U\left(0.35,0.55\right)$ ，四对随机序列 $\left({\alpha }_{ik},{\beta }_{ik}\right)$ ( $i=1,2,3,4$ )的变化见图1，并且该序列无法被窃听者识别。

本文涉及的过程噪声和测量噪声是与[20]类似的高斯混合噪声：

$\begin{array}{l}{\omega }_k=p_1{\text{Rand}}_{1\omega ,k}+\left(1-p_1\right){\text{Rand}}_{2\omega ,k},\\ {\nu }_{1,k} =p_1{\text{Rand}}_{1\nu ,1k}+\left(1-p_1\right){\text{Rand}}_{2\nu ,1k},{\nu }_{2,k}=p_2{\text{Rand}}_{1\nu ,2k}+\left(1-p_2\right){\text{Rand}}_{2\nu ,2k},\\ {\nu }_{3,k}=p_3{\text{Rand}}_{1\nu ,3k}+\left(1-p_3\right){\text{Rand}}_{2\nu ,3k},{\nu }_{4,k}=p_4{\text{Rand}}_{1\nu ,4k}+\left(1-p_4\right){\text{Rand}}_{2\nu ,4k}.\end{array}$

其中 $p_1=0.3,p_2=0.1,p_3=0.8,p_3=0.2$ ，并且上述随机变量满足以下分布：

$\begin{array}{l}{\text{Rand}}_{1\omega ,k}~\mathcal{N}\left(0.1,0.05\right),{\text{Rand}}_{2\omega ,k}~\mathcal{N}\left(0.4,0.09\right),\\ {\text{Rand}}_{1\nu ,ik}~\mathcal{N}\left(0.9,0.03\right),{\text{Rand}}_{1\nu ,ik}~\mathcal{N}\left(0.6,0.01\right).\end{array}$

Figure 1. Decomposition coefficient

图1. 分解系数

图2展示了系统在状态饱和时算法使用高斯核函数(红线)和S-t核函数(蓝线)估计系统状态的仿真结果，其中上标i ( $i\in 1,2,3,4$ )表示状态向量的第i个分量。图3展示了基于两种核函数下的系统状态估计的均方误差(MSE)。在该图中，MSE的定义为[21]：

$MS{E}_k=\frac{1}{T}{\displaystyle \underset{t=1}{\overset{T}{\sum }}{\displaystyle \underset{i=1}{\overset{4}{\sum }}{\left(x_k^i-{\hat{x}}_{k|k}^i\right)}^2}}$

其中，该式中T = 100代表独立实验的数量。

图2和图3展示了使用S-t核函数的算法的估计精度要高于使用高斯核函数的算法的估计精度。另外，除了算法的估计精度外，我们还要对系统的安全性进行分析。为此，假设外部窃听者使用估计方案(5)，并且假设窃听者除了无法获得时变参数的分解参数外，其他所有参数(包括设计增益)都能够轻易获得。实验结果如图4所示，该图展示了窃听者在估计器1和4的状态估计过程中，与真实系统状态之间的差异。

Figure 2. The estimation of the 1st and 4th components of the state variable by estimators 1 and 4

图2. 估计器1、4对状态变量的第1个和第4个分量的估计

Figure 3. The MSE of the estimator 1 for the first component of the state variable using the Gaussian kernel function and the S-t kernel function

图3. 估计器1对状态变量的第1个分量分别使用高斯核函数与S-t核函数的MSE

Figure 4. The estimation of the 1st and 4th components of the state variable by an external eavesdropper

图4. 外部窃听者对状态变量的第1个和第4个分量的估计

5. 结语

基于状态饱和与状态变量分解，本文提出了一种新的以S-t核函数的分布式安全状态估计算法。本文根据状态分解规则系统地分析了私有估计误差协方差矩阵与公共估计误差协方差矩阵之间的关系。同时，利用新的S-t核函数构造区别于高斯核函数的代价函数进而构建了估计器的增益，以优化整个估计过程。状态变量分解的引入有效保障了数据传输的安全性，防止外部窃听者从中提取敏感信息。最终，结合系统状态饱和的特征利用Matlab进行仿真实验，验证了所提出算法在安全性和有效性方面的优越性能。

基金项目

国家自然科学基金(No. 62303322)；上海市教委创新基金(No. XJ2024111)资助。