适合关键业务型工业物联网(Industrial internet of things, IIoT) 应用的工业无线传感器网络(Industrial wireless sensor networks, IWSN)需要部署在典型严酷工业环境中, 具有关键要求^[1]. 在IIoT中, 数据融合、能量管理、传输调度、定位和追踪协议要求所有的节点运行在公共的时间基准^[2]. 通过使用统计信号处理观点, 很多时钟同步协议归类为参数估计问题, 已经在不同的网络随机延迟分布下, 评估了一些经典协议的性能^[3]. 由于IWSN中存在许多不可靠的因素, 例如包延迟、丢失、节点失效等, 导致这些同步协议算法的性能大大地降低, 甚至算法失败. 特别指出, 国内学者针对时钟同步问题的前沿理论研究, 可参见文献[4]. 国际前沿研究中, 一致性算法^[5-6]在分析中没有考虑包延迟, 因此, 在收敛的时钟参数中产生了大的均方误差. 在文献[7]中, 鲁棒平均时间同步(Robust average time synchronization, RoATS)提出来相应于边界通信延迟, 以鲁棒的方法调节节点的时钟频率和时钟偏移. 在文献[8]中, 提出一组基于2阶一致性的随机时钟同步协议. 已经证明, 在这样的方法下, 可以总是调节算法参数, 且在概率均方意义下达到时钟同步. 然而, 一致性算法^[7-8]不能在分析中量化消息延迟的影响, 特别是包延迟. 作为耦合状态模型的相对时钟模型依赖于交互平均以达到同步补偿或协作控制过程. 因此, 已经存在的协议和一致性平均方法受限于网络规模, 随着网络规模的扩展, 时钟同步收敛性能将下降.

IIoT在规模网络下分布式精确时钟同步问题仍然是处于目前的研究热点. 作为耦合状态模型的相对时钟状态模型依赖于协作控制过程. 受限于网络尺寸的一致性时钟同步方法的收敛率研究集中在收敛性(稳定性), 难以量化和分析收敛率, 特别是对于不可靠网络.

当存在消息延迟时(固定时延), 用于文献[2]的卡尔曼滤波(Kalman filtering)导出了用于时钟偏移追踪的消息传递方法(双向信息交换). 事实上, 文献[2]的时钟模型具有新的不同本质(节点的时钟状态依赖于状态空间模型度量下自身节点状态的噪声扰动). 文献[2]的绝对时钟模型建模为解耦状态模型, 但量测模型仍然是绝对状态的辅助量测过程. 事实上, 存在状态耦合关系. 根据现代控制理论观点, 通过分析文献[2], 发现无法满足系统的能观测性(更坚实的理论依据参见附录A). 因此本文首先针对绝对时钟模型的能观测性的问题.

考虑包丢失(过度的延迟被认为是丢失)和节点失效作为节点量测信息丢失, 为了量化丢失的量测信息的影响, 已经广泛地采用两类通道模型: 1)独立同分布(Independent and identically distributed, i.i.d.)模型, 建模包丢失过程为一个独立同分布过程^[9]; 2)在马尔科夫模型下, 包处理过程描述为二元马尔科夫链^[10].

在独立同分布(i.i.d.)模型下, 文献[2]和文献[11$ - $13]集中在间隙卡尔曼滤波的稳定性, 此时仅一个传感器传输它的原始量测, 且存在一个临界丢包率. 在临界丢包率以上, 状态估计误差协方差矩阵的均值将发散为无穷^[9]. 在文献[9]中, 也给出了临界丢包率的上界和下界. 对于一般的向量系统, 众所周知难以显式地表达临界丢包率. 也是源于处理文献[9]的局限性, 文献[13]中对于一步可观测系统下界可显示为"紧"的, 包括文献[12]中所示的非退化系统. 然而, 文献[14]中的一个反例显示临界丢包率严格地位于下界和上界之间.

针对含多传感器的网络化状态估计问题, 许多研究者广泛研究了针对不同传感器, 约束于可能不同的包丢失情形^[15-20]. 在文献[21]中, 开发了一种新的回归矩阵技术, 用于研究针对于MMSE (Minimum mean square error)估计器的必要与充分稳定条件, 它适用于单一传感器和多传感器情形.

在"工业物联网"的前期研究中, 文献[4]针对不可靠量测的"相对时钟模型", 基于适合在线运算的Tubes-MPC (Model predictive control)控制方法完成"指数收敛性"量化分析. 文献[1]已经开展的工业物联网TDMA (Time division multiple access)确定性调度算法中, 维护高可靠性的时间精度具有重要意义. 文献[22]作为本文重要研究背景的"工业物联网"精确时钟同步的紧时隙需求, 本文提出的MMSE等价变换能观测性观测模型符合噪声量测的实际应用模型, 对文献[22]的理想观测器进行了深度分析扩展.

本文的建模与适用性: 1) 建模特点: 文献[2]的网络拓扑来源于实际的网络应用拓扑结构(如工业无线ISA100、WirelessHart、WIA-PA中的Mesh网络拓扑), 具体的网络协议可参考前期文献[22-24]和文献[25]. 2) 协议特点: 对于本文及文献[2]的分布式网络拓扑(全分布、短时延Cyber-Physical环境要求), 具体应用于无线Mesh网络(Wireless mesh network, WMN)^[23-24], 将在工业物联网TDMA独立时隙间实现无冲突的数据传输, 以有效提高资源利用率^[25].

为了不依赖于特定的相对参考节点或网络结构, 文献[2]中绝对状态模型建模为绝对状态空间模型, 状态相对于本地时钟的理想状态(状态空间"0"点). 本文通过对文献[2]的"能观测性‘的分析, 建立了MMSE等价变换的能观测性观测模型, 观测模型的能观测性是系统状态确定性估计的必要条件. 由提出的新的解耦观测模型, 建立了基于文献[2, 15]的规模网络的修正的卡尔曼滤波, 建立了时钟同步的量化分析的收敛条件.

本文的主要贡献是:

1) 基于前期的研究^{[4, 22, 26-27]}, 本文基于卡尔曼滤波的精确时钟同步研究, 从"网络化控制理论观点"准确定位文献[2]的研究实质(参见命题1), 分析了双向信息交换的物理量测模型和在BMU (Basic measurement unit)下建立的适用于大规模网络的"能观测性"精确量测模型.

2) 通过使用MMSE等价能观测性量测模型, 能够得到确定性的绝对状态的估计值(基于理想状态(1, 0)). 能够解除节点间的状态耦合与控制耦合. 时钟同步问题被转化为对于本地时钟理想状态(1, 0)的分布式卡尔曼滤波精确设定点状态追踪问题.

3) 通过MMSE等价能观测性的状态解耦的量测模型, 实现了对工业物联网的时钟同步精度的量化分析.

注1.   本文使用有$ N $个传感器节点$ \{S_1, S_2, $ $ \cdots, S_N\} $的网络为研究对象, 对于任意节点$ S_i $, $ i\in $ $ \{1, 2, \cdots, N\} $. $ \mathcal{N}_i $表示$ S_i $的邻居节点集合, $ \mathcal{N}_i = $ $ \{\mathcal{N}_i(j) = m_j| m_1, \cdots, m_j, \cdots, m_{|\mathcal{N}_i|}\} $ $ (m_j\in \mathcal{N}_i $, $ j $ $ \in \{1, \cdots, |\mathcal{N}_i|\}) $; $ \mathcal{N}_j $表示$ S_j $的邻居节点集合, $ \mathcal{N}_j $ $ = \{\mathcal{N}_j(o) = m_o| m_1, \cdots, m_o, \cdots, m_{|\mathcal{N}_j|}\} $ ($ m_o $ $ \in $ $ \mathcal{N}_j $, $ o\in\{1, \cdots, |\mathcal{N}_j|\} $). $ \mathcal{N}_{[i]} $表示包含$ {S_i} $及其邻居节点的局部子系统, $ \mathcal{N}_{[i]} = \{\mathcal{N}_i, i\} $, 元素个数$ |\mathcal{N}_{[i]}| $ $ = $ $ |\mathcal{N}_i|+1 $. $ (\cdot)^{M\times N} $表示任意$ M $行$ N $列的矩阵.

注2.  不可靠量测丢包变量定义, 参见表 1.

表 1  不可靠量测丢包变量定义

Table 1  Definition of unreliable packet loss variables

名称	定义
$\pi_{i, k}^{\{i, j\}}$	描述$S_i$和$S_j$间第$k$次信息交换. 交换成功, 则$\pi_{i, k}^{\{i, j\}}=1$; 否则$\pi_{i, k}^{\{i, j\}}=0$. 多链路下, $S_i$与任意邻居$S_{m_j}$间则为$\pi_{i, k}^{\{i, m_j\}}$
$\phi_{i, j}$	描述$S_i$和$S_j$间的接收概率, $P(\pi_{i, k}^{\{i, j\}}=1)=\phi_{i, j}$, $P(\pi_{i, k}^{\{i, j\}}$ $=0)=1-\phi_{i, j}$. 多链路下, $S_i$与任意邻居$S_{m_j}$间则为$\phi_{i, m_j}$
$\boldsymbol{\chi}$	表示一个同步周期, $S_i$与$|\mathcal{N}_i|$个邻居的$|\mathcal{N}_i|$次信息交换中, 丢包组合$\boldsymbol{\chi}=\{\chi_1, \chi_2, \cdots, \chi_g, \cdots, \chi_{2^{|\mathcal{N}_i|}}|g\in \{1$, $2$, $\cdots$, $2^{|\mathcal{N}_i|}\}\}$, $\chi_g$为第$g$种丢包情况, 共$2^{|\mathcal{N}_i|}$种
$\pi_{i, g, k}^{\{i, m_j\}}$	描述第$k$个采样周期第$g$种丢包情况时($\chi_g$), $S_i$与$S_{m_j}$信息交换是否成功地取值. 若成功, $\pi_{i, g, k}^{\{i, m_j\}}=1$, 否则$\pi_{i, g, k}^{\{i, m_j\}}$ $=$ $0$
$\boldsymbol{L}_{i, g, k}$	如$\pi_{i, g, k}^{\{i, m_j\}}$, 则$\boldsymbol{L}_{i, g, k}={\rm diag}\{\pi_{i, g, k}^{\{i, m_1\}}, \cdots, \pi_{i, g, k}^{\{i, m_j\}}$, $\cdots$, $\pi_{i, g, k}^{\{i, m_{|\mathcal{N}_i|}\}}\}$, 表示$S_i$与所有邻居在$\chi_g$情况时的丢包系数矩阵
$\eta_{i, g, k}$	与所有邻居在$\chi_g$时的概率$\eta_{i, g, k}=\prod_{j=1}^{|\mathcal{N}_i|} \alpha_{i, m_j, k}$
$\alpha_{i, m_j, k}$	表示$S_i$与$S_{m_j}$信息交换概率
	$\alpha_{i, m_j, k}=\begin{cases} \phi_{i, m_j}, &\mbox{若}~ \pi_{i, k}^{\{i, m_j\}}=1\\ 1-\phi_{i, m_j}, &\mbox{若}~ \pi_{i, k}^{\{i, m_j\}}=0 \end{cases}$
$\boldsymbol{\gamma}_{i, g, k}$	第$k$步绝对量测信息$\boldsymbol{\gamma}_{i, k}$在$\chi_g$情况时的具体表达式, $\boldsymbol{\gamma}_{i, g, k}$为构造的值, 实际量测中不可获得
$\boldsymbol{z}_{i, g}(k)$	第$k$步相对量测信息$\boldsymbol{z}_{i}(k)$在$\chi_g$情况时的具体表达式, $\boldsymbol{z}_{i, g}(k)$为一个实际中能够量测的值

下载: 导出CSV 

| 显示表格

1.   系统解耦状态模型

文献[2]中, 考虑了传感器时钟偏斜和累积时钟偏移的演化模型

时钟状态$ \mathit{\boldsymbol{x}}_i(k) = \begin{bmatrix} \beta_i(k)\\ \vartheta_i(k) \end{bmatrix} $, $ \beta_i(k) $和$ \vartheta_i(k) $分别是在第$ k\Delta $次采样中的瞬时偏斜和累积偏移. 状态转移矩阵$ \mathit{\boldsymbol{A}}_i = \begin{bmatrix} 1&0\\ \Delta \tau_{0}&1 \end{bmatrix} $, $ \mathit{\boldsymbol{w}}_i(k) $是高斯噪声, 且均值为0和协方差矩阵$ \mathit{\boldsymbol{Q}}_i = \sigma_{q}^{2} $ $ \times $ $ \begin{bmatrix}\Delta&0\\0&\Delta(\Delta+1)(2\Delta+1)\tau_{0}^{2}/6\end{bmatrix} $, 并且状态转移常数$ \mathit{\boldsymbol{b}}_i $ $ = $ $ \begin{bmatrix}0\\ -\Delta \tau_0 \end{bmatrix} $. $ \tau_0 $是采样周期, $ \Delta \tau_0 $ $ (\Delta\gg 1) $称为同步周期.

为了达到全网的时钟同步, 所有节点的时钟读数$ c_i(k) $必须被调节为公共的值. 假定UTC (Universal time coordinated)作为准确的参考时间, 也就是$ \beta(k) = 1 $和$ \vartheta(k) = 0 $, 并且基于离散时钟读数模型$ c_i(k) = k\Delta \tau_0+\vartheta(k-1)+[\beta(k-1)-1]\Delta \tau_0 $, 全网时钟同步任务将对应于UTC追踪时变的时钟偏斜$ \{\beta_i(k)\}_{i = 1}^{N} $和累积偏移$ \{\vartheta_i(k)\}_{i = 1}^{N} $. 对比于文献[2], 假定节点$ S_1 $为具有准确时钟的参考节点实现绝对时钟量测模型是不合理的(更坚实的理论依据参见附录A). 方程(1)中的节点时钟状态建模源于自身节点扰动的物理模型(而不依赖于其他节点), 称方程(1)为绝对时钟状态模型.

2.   本地能观测性解耦观测器

2.1   大规模系统本地观测器

现有研究均采用相对时钟状态模型, 及相对时钟量测模型以获取时钟状态的观测值. 因采用相对时钟模型, 则收敛性能依赖于网络规模. 本节将讨论我们提出的从文献[2]的绝对状态空间模型出发, 构建能观测性系统量测方程的新方法.

节点$ {S_i} $和它的邻居节点$ {S_j} $之间的双向时间戳交换如图 1. 不同于文献[2], 我们在量测关系上平等地对待$ {S_i} $和$ {S_j} $, 例如式(10), 同时存在$ z_{i, k}^{\{ i, j\} } $和$ z_{j, k}^{\{ j, i\} } $.

图 1  一组时钟信息交换过程($z_{i, k}^{\{i, j\}}$和$z_{j, k}^{\{j, i\}}$互为对称量测信息)

Fig. 1  A set of clock information exchange processes ($z_{i, k}^{\{i, j\}}$ and $z_{j, k}^{\{j, i\}}$ are symmetrical measurement information each other)

下载: 全尺寸图片 幻灯片

事实上, IWSN有一系列的不可靠的因素, 会对时间戳丢失产生影响. 因此, 我们必须考虑时间戳丢失的影响(全面的丢包假设参见文献[21, 28]; 可扩展说明参见第3节). 在一轮双向信息交换过程中(如图 1所示), 无论哪一部分的时间戳丢失, 它将影响时间戳交换轮次的成功, 只要两个相邻节点在有限的时间内没有接收到相应的时间戳, 就认为这一轮时钟信息交换失败, 也就是说, 包发生了丢失.

对于时间戳交换事件仅有两种情形: 交换成功或者失败^[15]. 使用伯努利随机变量$ \pi_{i, k}^{\{i, j\}} $ ($ \pi_{j, k}^{\{j, i\}} = $ $ \pi_{i, k}^{\{i, j\}} $)描述节点$ S_i $和节点$ S_j $间第$ k $次信息交换的丢包事件. 如果成功, 则$ \pi_{i, k}^{\{i, j\}} = 1 $; 否则$ \pi_{i, k}^{\{i, j\}} = 0 $ (见表 1). 假定随机变量$ \pi_{i, k}^{\{i, j\}} $ $ (k\in {\bf N}) $遵守独立同分布^[15], 则$ P(\pi_{i, k}^{\{i, j\}} = 1) = \phi_{i, j} $, $ P(\pi_{i, k}^{\{i, j\}} = 0) = 1 $ $ - $ $ \phi_{i, j} $, $ \phi_{i, j} $为节点$ S_i $和节点$ S_j $间的接收概率; 假定随机变量$ \pi_{i, k}^{\{i, j\}} $和$ \pi_{l, k}^{\{l, n\}} $彼此独立$ (l\neq i $且$ n\neq j $或$ l\neq j $且$ n\neq i) $, 则$ P(\pi_{i, k}^{\{i, j\}}\pi_{l, k}^{\{l, n\}}) = P(\pi_{i, k}^{\{i, j\}}) $ $ \times $ $ P(\pi_{l, k}^{\{l, n\}}) $.

2.1.1   问题的提出

根据附录A的详细分析, 文献[2]利用双向信息交换机制建立的系统量测模型, 结合完全解耦的节点时钟状态转移方程(1), 由此得到的系统时钟状态空间模型, 无法满足系统状态的能观测性.

由附录A, 从现代控制理论视点, 对于子系统的量测模型, 不管是一对邻居节点的一组量测, 还是多对邻居节点的多组量测, 子系统的量测矩阵$ \mathit{\boldsymbol{C}} $都不可逆, 即根据当前步的观测值不能马上确定系统当前步的状态, 进一步得知子系统的能观测性矩阵$ \mathit{\boldsymbol{W}}_o $的秩不为$ 2|\mathcal{N}_{[i]}| $ $ (Rank(\mathit{\boldsymbol{W}}_o) = 2|\mathcal{N}_{[i]}|-2) $, 即无论多少组观测值都不能确定子系统当前步的状态, 换句话说, 子系统是不可观测的.

2.1.2   物理过程分析、策略、MMSE等价变换

根据现代控制理论, 通过量测值确定系统内部的状态, 正确地选择合适的量测模型是保证系统能观测性的关键. 为此, 我们引入BMU, 深入研究双向信息交换的物理过程本质, 希望能够建立一个辅助的绝对时钟量测模型.

结合辅助量测模型, 通过线性变换, 将文献[2]的相对量测模型转换为一个具有MMSE性能等价的解耦量测模型, 从而保证系统的能观测性. 参考这个等价变换过程, 构造一个修改的卡尔曼滤波观测器. 将其描述为下面的定理1.

定理1.  考虑式(8), (10), (11)和式(15), BMU的能观测性已构建MMSE等价性能量测的必要条件: 建立方程(8)$ \sim $(16)为系统能观测性量测模型的演化线索. 在BMU分析模型下, 构造符合物理过程的辅助量测模型, 通过MMSE量测性能等价变换, 建立实现子系统能观测性的量测模型, 实现子系统的能观测性.

已经证明(参见定理2), 通过将方程(11)中的独立量测噪声线性转换为相关量测噪声(由文献[29, 6.6节]信号处理实例, 易获得量测"相关噪声"的理解). 运用分布式卡尔曼滤波能够实施MMSE量测.

证明.

1) 能观测性策略. BMUs能观测性的清晰证明可以参考附录B. 这一部分是方程(8)$ \sim $(16)实现能观测性证明的精炼版本, 提示我们理解关键点.

2) 能观测性的量测模型–物理过程分析. 这一部分将分析双向信息状态量测的物理模型与MMSE等价线性变换.

对于时钟读数模型, 可以表达为标准时间和累积时钟偏移形式

根据如图 1所示的信息交换的过程, 仅考虑信息交换各自时间之间的关系, 也就是$ T_{2, k}^{\{ i, j\} } $和$ T_{3, k}^{\{ i, j\} } $间, $ T_{1, k}^{\{ i, j\} } $和$ T_{4, k}^{\{ i, j\} } $间的关系, 能够得到下面的方程

当节点$ {S_i} $在第$ k $次同步周期发送时钟信息时, 这里$ k\Delta {\tau _0} $是参考时间, 而$ d_{ij} $是节点$ S_i $和$ S_j $间的固定时延(且$ {d_{ij}} $是对称的^{[2, 4]}), $ D $是量测节点$ S_j $发送应答信息的相应时间; $ X_{k}^{\{i, j\}} $是在节点$ S_i $和节点$ S_j $传输过程的随机时延; $ Y_{k}^{\{i, j\}} $是在节点$ S_j $和节点$ S_i $传输过程的随机时延.

注意到, 由统计信号处理的建模观点, 由文献[29]中的6.6节, 文献[2]的建模参数能够通过技术方法获得(例如, 量测方程(文献[2]中的式(17))的噪声量测的统计方差). 而观测方程(7)和(8) 则无法直接获得. 因此, 文献[2]中的式(17)做为建模的假设前提是合理的.

1) 辅助参考绝对量测状态是真实存在的, 但不可直接获得

考虑随机时延$ X_{k}^{\{i, j\}} $和$ {\rm Y}_{k}^{\{i, j\}} $, 包含几个相互独立的随机过程, 按照中心极限定理, $ X_{k}^{\{i, j\}} $和$ Y_{k}^{\{i, j\}} $符合高斯延迟模型. 文献[30]证实了将随机时延建为高斯分布模型的置信度为99.8 %. 因此, 假设随机时延$ X_{k}^{\{i, j\}} $和$ Y_{k}^{\{i, j\}} $服从均值为0, 方差为$ \sigma _v^2 $的高斯分布. 根据式(3)和式(6), 节点$ S_i $的观测方程为

如果每一个节点$ S_j $从节点$ S_i $处接收到发送的信息, 响应立即给出, 我们能够使$ D = 0 $. 关于节点$ S_i $邻居的所有的观测方程(7)组集在一起, 能够得到

式中, $ \mathit{\boldsymbol{x}}(k) = [\mathit{\boldsymbol{x}}_{1}^{\rm T}(k), \cdots, \mathit{\boldsymbol{x}}_{i}^{\rm T}(k), \cdots, \mathit{\boldsymbol{x}}_{N}^{\rm T}(k) ]^{\rm T} $ $ \in $ $ {\bf R}^{2N\times 1} $是全局状态向量; $ \mathit{\boldsymbol{\pi}}_{i, k} = {\rm diag}\{\pi_{i, k}^{\{i, m_1\}}, \cdots, $ $ \pi_{i, k}^{\{i, m_j\}}, \cdots , \pi_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}}\}\in {\bf R}^{|\mathcal{N}_i|\times |\mathcal{N}_i|} $是丢包矩阵; $ \mathit{\boldsymbol{\gamma}}_{i, k} = [\mathit{\boldsymbol{\gamma}}_{i, k}^{\{i, m_1\}}, \cdots, \mathit{\boldsymbol{\gamma}}_{i, k}^{\{i, m_j\}}, \cdots, \mathit{\boldsymbol{\gamma}}_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}} ]^{\rm T} $ $ \in $ $ {\bf R}^{|\mathcal{N}_i|\times 1} $是系统的状态量测; $ \mathit{\boldsymbol{M}}_i\in {\bf R}^{|\mathcal{N}_i|\times 2N} $是量测矩阵, 满足$ \mathit{\boldsymbol{M}}_i(j, m) = \begin{cases} 2, &{\rm{若}}\; m = 2i\\ 0, &{\rm{否则}} \end{cases} $; $ \mathit{\boldsymbol{\zeta}}_i(k) $ $ = $ $ [\mathit{\boldsymbol{\zeta}}_{i, m_1}(k), \cdots, \mathit{\boldsymbol{\zeta}}_{i, m_j}(k), \cdots, \mathit{\boldsymbol{\zeta}}_{i, m_|\mathcal{N}_i|}(k)]^{\rm T}\in {\bf R}^{|\mathcal{N}_i|} $是量测噪声, 假定一次量测中节点$ S_i $与所有邻居节点进行信息交换的量测噪声相同, 则可改为$ \mathit{\boldsymbol{\zeta}}_{i, m_1}(k) $ $ = $ $ \cdots = \mathit{\boldsymbol{\zeta}}_{i, m_j}(k) = \cdots = \mathit{\boldsymbol{\zeta}}_{i, m_{|\mathcal{N}_i|}}(k) $, 符合0均值和协方差$ { R}_{\zeta_i} $.

针对节点$ S_i $的局部子系统的绝对时钟状态量测方程可表达为

式中, $ \mathit{\boldsymbol{x}}_{\mathcal{N}_{[i]}}(k) = \mathit{\boldsymbol{\Lambda}}_i\pmb{x}(k)\in {\bf R}^{2|\mathcal{N}_{[i]}|\times 1} $表示局部子系统时钟状态向量, 包含节点$ S_i $及其所有邻居节点$ S_{m_j} $的状态向量, 且$ \mathit{\boldsymbol{x}}_{\mathcal{N}_{[i]}}(k) = \mathit{\boldsymbol{\Lambda}}_i\pmb{x}(k) $; 量测系数向量$ \tilde{\mathit{\boldsymbol{M}}}_i = \mathit{\boldsymbol{M}}_i\pmb{\Lambda}_{i}^{\rm T} = [ \tilde{\mathit{\boldsymbol{M}}}_{i, m_1}^{\rm T}, \cdots, \tilde{\mathit{\boldsymbol{M}}}_{i, m_j}^{\rm T}, \cdots $, $ \tilde{\mathit{\boldsymbol{M}}}_{i, m_{|\mathcal{N}_i|}}^{\rm T} ]^{\rm T} $是矩阵$ \mathit{\boldsymbol{M}}_i $排除所有节点$ {S}_i $非邻居节点对应列后的约化矩阵. 且$ \mathit{\boldsymbol{\Lambda}}_{i} \in {\bf R}^{2|\mathcal{N}_{[i]}|\times 2N} $是对应于节点$ S_i $的提取矩阵, 满足

式中, $ n \in \{ 1, 2, $ $ \cdots, 2|{{{N}}_{[i]}}| \} $, 且$ m \in \{ 1, 2, \cdots $, $ 2N \} $.

因为存在标准时间$ k\Delta {\tau _0} $和固定时延, 不能直接获得$ \mathit{\boldsymbol{\gamma}}_{i, k} $, 因此, 需要另一个量测模型. 根据式(3)$ \sim $ (6), 我们能够得到节点$ S_i $的观测模型(对称量测信息$ \mathit{\boldsymbol{z}}_{j, k}^{\{j, i\}} $可以通过另一组时间戳获得, 存在于辅助量测节点中)

2) 双向信息交换、采集辅助相对状态量测

以节点$ S_i $为中心的所有观测方程(10) (即方程组(10)的第1个方程)进行组集:

式中, $ \mathit{\boldsymbol{z}}_i(k) = [\mathit{\boldsymbol{z}}_{i, k}^{\{i, m_1\}}, \cdots, \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_j\}}, \cdots, $ $ \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $ $ \in $ $ {\bf R}^{|\mathcal{N}_i|} $是系统的状态量测. 在实际中, $ \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_j\}} $和对应的对称量测信息$ \mathit{\boldsymbol{z}}_{m_j, k}^{\{m_j, i\}} $, 用于更新绝对状态辅助量测, 并通过辅助量测的更新消息彼此互传(参见附录B). $ \mathit{\boldsymbol{H}}_i\in {\bf R}^{|\mathcal{N}_i|\times 2N} $是量测矩阵, 满足

$ {{\pmb\upsilon} _i}(k) \in {{\bf{R}}^{{{|}}{{{N}}_{{i}}}{{|}}}} $是量测噪声, 符合0均值和协方差矩阵为$ {{\pmb R}_{{\upsilon _i}}} $.

注意到矩阵$ {\pmb H}_i $的特殊结构, 能够从观测模型(11)中提取它自身节点和所有邻居节点的状态

式中, $ \tilde{\mathit{\boldsymbol{H}}}_i\in {\bf R}^{|\mathcal{N}_i|\times 2|\mathcal{N}_{[i]}|} $是局部子系统的观测矩阵, 且$ \tilde{\mathit{\boldsymbol{H}}}_i = \mathit{\boldsymbol{H}}_i\pmb{\Lambda}_{i}^{\rm T} $.

3) 组集参考绝对量和相对量的量测

对于所有$ S_i $, 式(8)和式(11)作为行交替地组集在一起, 我们能够得到全局量测模型

式中, $ \mathit{\boldsymbol{\lambda}}(k) = {\rm blkdiag}\{(\mathit{\boldsymbol{\lambda}}_1(k) $, $ \cdots $, $ \mathit{\boldsymbol{\lambda}}_i(k) $, $ \cdots, $ $ \mathit{\boldsymbol{\lambda}}_N(k)\} $是全局丢包矩阵, $ \mathit{\boldsymbol{\lambda}}_i(k) = {\rm diag}\{\pi_{i, k}^{\{i, m_1\}}, $ $ \pi_{i, k}^{\{i, m_1\}} $, $ \cdots $, $ \pi_{i, k}^{\{i, m_j\}} $, $ \pi_{i, k}^{\{i, m_j\}} $, $ \cdots $, $ \pi_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}} $, $ \pi_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}})\} $; $ \mathit{\boldsymbol{Z}}(k) = [\mathit{\boldsymbol{Z}}_{1}^{\rm T}(k), \cdots, \mathit{\boldsymbol{Z}}_{i}^{\rm T}(k), \cdots $, $ \mathit{\boldsymbol{Z}}_{N}^{\rm T}(k)]^{\rm T} $是状态量测, $ \mathit{\boldsymbol{Z}}_{i}(k) = [\mathit{\boldsymbol{\gamma}}_{i, k}^{\{i, m_1\}}, \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_1\}}, \cdots $, $ \mathit{\boldsymbol{\gamma}}_{i, k}^{\{i, m_j\}}, $ $ \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_j\}}, \cdots, \mathit{\boldsymbol{\gamma}}_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}}, \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $ $ \in $ $ {\bf R}^{2|\mathcal{N}_i|} $; $ \mathit{\boldsymbol{G}}(k) $ $ = $ $ [\mathit{\boldsymbol{G}}_{1}^{\rm T}(k), \cdots, \mathit{\boldsymbol{G}}_{i}^{\rm T}(k), \cdots, \mathit{\boldsymbol{G}}_{N}^{\rm T}(k)]^{\rm T} $是量测矩阵.

另外, $ \mathit{\boldsymbol{\varpi}}(k) = [\mathit{\boldsymbol{\varpi}}_{1}^{\rm T}(k), \cdots, \mathit{\boldsymbol{\varpi}}_{i}^{\rm T}(k), \cdots, $ $ \mathit{\boldsymbol{\varpi}}_{N}^{\rm T}(k)]^{\rm T} $是0均值, 方差为$ \mathit{\boldsymbol{R}}_{\varpi} $的高斯噪声, $ \mathit{\boldsymbol{R}}_{\varpi} = {\rm blkdiag}\{\mathit{\boldsymbol{R}}_{\varpi_{1}}, \mathit{\boldsymbol{R}}_{\varpi_{2}}, $ $ \cdots, \mathit{\boldsymbol{R}}_{\varpi_{N}}\} $, $ \mathit{\boldsymbol{\varpi}}_{i}(k)\in {\bf R}^{2|\mathcal{N}_i|} $是0均值, 且方差为$ \mathit{\boldsymbol{R}}_{\varpi_{i}}\in {\bf R}^{2|\mathcal{N}_i|\times 2|\mathcal{N}_i|} $的高斯噪声.

4) 局部子系统量测的MMSE等价线性变换

在双向信息交换和对称量测信息条件下, 因为$ \pi_{j, k}^{\{j, i\}} = \pi_{i, k}^{\{i, j\}} $, 对量测模型(13)实施线性行变换, 能够得到特定结构的新的量测模型

式中, $ \mathit{\boldsymbol{Y}}(k) = [\mathit{\boldsymbol{Y}}_{1}^{\rm T}(k), \cdots, \mathit{\boldsymbol{Y}}_{i}^{\rm T}(k), \cdots, $ $ \mathit{\boldsymbol{Y}}_{N}^{\rm T}(k)]^{\rm T} $是状态量测, $ \mathit{\boldsymbol{Y}}_{i}(k) = \mathit{\boldsymbol{P}}_i\pmb{Z}(k)\in {\bf R}^{2|\mathcal{N}_i|} $; $ \mathit{\boldsymbol{\Gamma}} = [\mathit{\boldsymbol{\Gamma}}_{1}^{\rm T} $, $ \cdots, \mathit{\boldsymbol{\Gamma}}_{i}^{\rm T}, \cdots, \mathit{\boldsymbol{\Gamma}}_{N}^{\rm T}]^{\rm T} $是量测矩阵, $ \mathit{\boldsymbol{\Gamma}}_i = \mathit{\boldsymbol{P}}_i\pmb{G}\in {\bf R}^{2|\mathcal{N}_i|\times 2N} $, 满足$ \mathit{\boldsymbol{\Gamma}}_i(n, m) = \begin{cases} 2&{\rm{若}}\ m = 2i\\ 0&{\rm{否则}} \end{cases} $.

另外, $ \mathit{\boldsymbol{\psi}}(k) $ $ = $ $ [\mathit{\boldsymbol{\psi}}_{1}^{\rm T}(k) $, $ \cdots $, $ \mathit{\boldsymbol{\psi}}_{i}^{\rm T}(k) $, $ \cdots, $ $ \mathit{\boldsymbol{\psi}}_{N}^{\rm T}(k)]^{\rm T} $是0均值, 方差为$ \mathit{\boldsymbol{R}}_{\psi} $的高斯噪声, $ \mathit{\boldsymbol{\psi}}_{i}(k) $ $ = $ $ \mathit{\boldsymbol{P}}_i\pmb{\varpi}(k)\in {\bf R}^{2|\mathcal{N}_i|\times 1} $且方差为$ \mathit{\boldsymbol{R}}_{\psi} $ $ = $ $ \mathit{\boldsymbol{P}}\mathit{\boldsymbol{R}}_{\varpi}\mathit{\boldsymbol{P}}^{\rm T} $, 其中$ \mathit{\boldsymbol{P}}(k) = [\mathit{\boldsymbol{P}}_{1}^{\rm T}(k), \cdots $, $ \mathit{\boldsymbol{P}}_{i}^{\rm T}(k), \cdots, \mathit{\boldsymbol{P}}_{N}^{\rm T}(k)]^{\rm T} $, 且$ \mathit{\boldsymbol{P}} $是可逆矩阵, $ \mathit{\boldsymbol{P}}_i\in {\bf R}^{2|\mathcal{N}_i|\times 2\sum_{i = 1}^{N}|\mathcal{N}_i|} $.

式中, $ n $ $ \in $ $ \{1 $, $ 2 $, $ \cdots $, $ 2|\mathcal{N}_{i}|\} $, $ m $ $ \in $ $ \{1 $, $ 2 $, $ \cdots $, $ \sum_{k = 1}^{N}2|\mathcal{N}_k|\} $. $ N_j(o) $参见注1. $ \mathit{\boldsymbol{P}}_i(n, m) $直观的理解参见附录D-2中的"物理变换矩阵".

5) 能观测性系统量测方程

在实际的大规模网络中, 通过使用方程(13)的信息交换, 能够实现如方程(14)所示的MMSE等价变换(附录B). 为了证明量测方程(11)在辅助量测方程(8)的协助下存在MMSE量测性能解耦, 对应于式(11), 在式(14)中抽取偶数行, 我们得到一个解耦的量测方程, 在MMSE量测的等价变换上等价于能观测性方程(13)$ \sim $(15).

式中, $ \mathit{\boldsymbol{y}}_{i, k} = \mathit{\boldsymbol{V}}_i\pmb{Y}_i\in {\bf R}^{|\mathcal{N}_i|} $是系统状态量测; $ \mathit{\boldsymbol{C}}_i $ $ = $ $ \mathit{\boldsymbol{V}}_i\pmb{\Gamma}_i $ $ \in $ $ {\bf R}^{|\mathcal{N}_i|\times 2N} $是量测矩阵, 满足

$ \mathit{\boldsymbol{\nu}}_i(k) = V_i\pmb{\psi}_{i} \in {\bf R}^{|\mathcal{N}_i|} $是量测噪声, 符合0均值和协方差矩阵$ \mathit{\boldsymbol{R}}_{\nu_i} $ $ = $ $ E[\mathit{\boldsymbol{\nu}}_i(k) \mathit{\boldsymbol{\nu}}_i^{\rm T}(k)] = \mathit{\boldsymbol{V}}_i\pmb{P}_i\pmb{R}_{\varpi} \mathit{\boldsymbol{P}}_i^{\rm T}\pmb{V}_i^{\rm T} $, 并且$ \mathit{\boldsymbol{V}}_i $ $ \in $ $ {\bf R}^{|\mathcal{N}_i|\times 2|\mathcal{N}_i|} $是对应于节点$ S_i $的提取矩阵, 满足

注意到矩阵$ \mathit{\boldsymbol{C}}_i $的特殊结构, $ \mathit{\boldsymbol{C}}_i $中对应于其他节点的值(除了$ S_i $)是0. 因此, 对于节点$ S_i $的解耦量测模型为

式中, $ \bar{\mathit{\boldsymbol{C}}}_i $是解耦量测矩阵, $ \bar{\mathit{\boldsymbol{C}}}_i = \begin{bmatrix} 0&2\\ 0&2\\ \vdots&\vdots \end{bmatrix}\in {\bf R}^{|\mathcal{N}_i|\times 2} $.

为了区分观测方程(11)和(15), 我们称观测方程(11)是耦合观测方程, 且观测方程(15)是解耦观测方程. 在下面的设计过程中, 将解耦观测方程视为理论分析工具.

定理2.  在BMU能观测性量测模型的必要条件下, 系统存在MMSE性能等价变换.

证明.

1) 基于BMU的MMSE等价变换(拆分观点).

证明参见附录C.

2) 多量测BMU的MMSE等价变换(含丢包矩阵的拆分观点).

证明参见附录D.

2.2   含量测丢包的卡尔曼滤波

对于本文的多链路观测模型, 任一节点$ S_i $共有$ |\mathcal{N}_i| $个邻居节点, 则在对称双向信息交换时钟同步机制下的一次时钟同步信息交换中, $ S_i $共进行了$ |\mathcal{N}_i| $次信息交换. 其中每一次信息交换都有可能会发生丢包事件, $ \chi_g $表示第$ g $种丢包情况(参见表 1).

由$ \mathit{\boldsymbol{L}}_{i, g, k} $ $ = $ $ {\rm diag}\{\pi_{i, g, k}^{\{i, m_1\}} $, $ \cdots $, $ \pi_{i, g, k}^{\{i, m_j\}} $, $ \cdots, $ $ \pi_{i, g, k}^{\{i, m_{|\mathcal{N}_i|}\}}\} $ ($ \mathit{\boldsymbol{L}}_{i, g, k} $参见附录D-1), 将式(9)和式(12)扩展为

式中, $ \tilde{\mathit{\boldsymbol{M}}}_{i, g, k} $ $ = $ $ \mathit{\boldsymbol{L}}_{i, g, k} [\tilde{\mathit{\boldsymbol{M}}}_{i}^{\{i, m_1\}}, \cdots, $ $ \tilde{\mathit{\boldsymbol{M}}}_{i}^{\{i, m_j\}}, \cdots, $ $ \tilde{\mathit{\boldsymbol{M}}}_{i}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $和$ \tilde{\mathit{\boldsymbol{H}}}_{i, g, k} $ $ = $ $ \mathit{\boldsymbol{L}}_{i, g, k} [\tilde{\mathit{\boldsymbol{H}}}_{i}^{\{i, m_1\}} $, $ \cdots, $ $ \tilde{\mathit{\boldsymbol{H}}}_{i}^{\{i, m_j\}}, \cdots, \tilde{\mathit{\boldsymbol{H}}}_{i}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $表示绝对和相对量测矩阵$ \tilde{\mathit{\boldsymbol{M}}}_i $和$ \tilde{\mathit{\boldsymbol{H}}}_i $中含丢包信息交换下的量测矩阵, 其中$ \mathit{\boldsymbol{\zeta}}_{i, g}(k) = \mathit{\boldsymbol{L}}_{i, g, k}[\mathit{\boldsymbol{\zeta}}_i(k), \mathit{\boldsymbol{\zeta}}_i(k), \cdots]^{\rm T} $和$ \mathit{\boldsymbol{\upsilon}}_{i, g}(k) = $ $ \mathit{\boldsymbol{L}}_{i, g, k}[\mathit{\boldsymbol{\upsilon}}_i(k), \mathit{\boldsymbol{\upsilon}}_i(k), \cdots]^{\rm T} $表示量测噪声.

根据式(16) $ \sim $ (18), 可得$ \chi_g $下的能观性量测方程

式中, $ \bar{\mathit{\boldsymbol{C}}}_{i, g, k} $ $ = $ $ \mathit{\boldsymbol{L}}_{i, g, k}[\bar{\mathit{\boldsymbol{C}}}_{i}^{\{i, m_1\}} $, $ \cdots $, $ \bar{\mathit{\boldsymbol{C}}}_{i}^{\{i, m_j\}} $, $ \cdots, $ $ \bar{\mathit{\boldsymbol{C}}}_{i}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $表示解耦量测方程的量测矩阵$ \bar{\mathit{\boldsymbol{C}}}_i $中含丢包信息交换下的量测矩阵, $ \mathit{\boldsymbol{\nu}}_{i, g}(k) = $ $ \mathit{\boldsymbol{L}}_{i, g, k} $ $ \times $ $ [\mathit{\boldsymbol{\nu}}_i(k), \mathit{\boldsymbol{\nu}}_i(k), \cdots]^{\rm T} $表示量测噪声, $ \mathit{\boldsymbol{R}}_{\nu_{i, g, k}} $表示量测噪声协方差, $ \mathit{\boldsymbol{R}}_{\nu_{i, g, k}} = \mathit{\boldsymbol{L}}_{i, g, k}\times \mathit{\boldsymbol{R}}_{\nu_{i, k}} $.

使用解耦观测模型(15)和耦合观测模型(11), 能够建立解耦观测器和耦合观测器. 并针对观测器量测的MMSE性能等价变换进行了深入分析. 不同于文献[2], 已严格证明存在卡尔曼滤波的MMSE等价变换(参见定理2).

考虑系统的状态空间方程(1)和解耦观测模型(16), 沿着文献[9, 15]中的分析, 我们组合了所有的丢包情形, 并创建了修正的卡尔曼滤波:

式中, $ \hat{\mathit{\boldsymbol{x}}}_i(k|k) $和$ \hat{\mathit{\boldsymbol{x}}}_i(k|k-1) $表示$ k $时刻的估计状态和预报状态, 且$ \mathit{\boldsymbol{P}}_i(k|k) $和$ \mathit{\boldsymbol{P}}_i(k|k-1) $是后向和前向误差协方差矩阵. 集散式观测器(20) $ \sim $ (24)用于理论分析.

考虑系统状态方程和第3.1.2节中耦合量测方程(12), 以及辅助绝对量测方程(9), 能够建立修正的分布式卡尔曼滤波. 耦合量测方程(12)仅处理它自身节点和所有邻居节点的状态. 修正卡尔曼滤波的量测更新将写为下面的分布式形式

式中, $ \mathit{\boldsymbol{\gamma}}_{\mathcal{N}_i, g, k} = \mathit{\boldsymbol{L}}_{i, g, k}\times \mathit{\boldsymbol{\gamma}}_{\mathcal{N}_i, k} = \mathit{\boldsymbol{L}}_{i, g, k} [\mathit{\boldsymbol{\gamma}}_{m_1, k}^{\{m_1, i\}} $, $ \cdots $, $ \mathit{\boldsymbol{\gamma}}_{m_j, k}^{\{m_j, i\}}, \cdots, \mathit{\boldsymbol{\gamma}}_{m_{|\mathcal{N}_i|}, k}^{\{m_{|\mathcal{N}_i|}, i\}}]^{\rm T} $, $ \mathit{\boldsymbol{\gamma}}_{m_j, k}^{\{m_j, i\}} $是节点$ S_{m_j} $的绝对状态量测信息, 且$ \tilde{\mathit{\boldsymbol{M}}}_{\mathcal{N}_i, g, k} = \mathit{\boldsymbol{L}}_{i, g, k} $ $ \times $ $ \tilde{\mathit{\boldsymbol{M}}}_{\mathcal{N}_i} = \mathit{\boldsymbol{L}}_{i, g, k}\times [\mathit{\boldsymbol{M}}_{m_1}^{\{m_1, i\}}, \cdots, \mathit{\boldsymbol{M}}_{m_j}^{\{m_j, i\}} $, $ \cdots $, $ \mathit{\boldsymbol{M}}_{m_{|\mathcal{N}_i|}}^{\{m_{|\mathcal{N}_i|}, i\}}]^{\rm T}\pmb{\Lambda}^{\rm T} $, $ \mathit{\boldsymbol{M}}_{m_j}^{\{m_j, i\}} $是量测矩阵$ \mathit{\boldsymbol{M}}_{m_j} $对应于$ \mathit{\boldsymbol{\gamma}}_{m_j, k}^{\{m_j, i\}} $的行.

注1.  在实际算法中, 实际将用$ \bar{\mathit{\boldsymbol{\gamma}}}_{\mathcal{N}_i, k} $替换$ \mathit{\boldsymbol{\gamma}}_{\mathcal{N}_i, k} $, 满足$ \bar{\mathit{\boldsymbol{\gamma}}}_{\mathcal{N}_i, k+1} = \bar{\mathit{\boldsymbol{\gamma}}}_{\mathcal{N}_i, k} +\mathit{\boldsymbol{z}}_{\mathcal{N}_i, k}-\mathit{\boldsymbol{z}}_{\mathcal{N}_i, k+1} $, $ \bar{\mathit{\boldsymbol{\gamma}}}_{\mathcal{N}_i, 0} $ $ = $ $ [(T_{1, 0}^{\{m_1, i\}}+T_{4, 0}^{\{m_1, i\}}-2d_{m_1, i}), \cdots $, $ (T_{1, 0}^{\{m_j, i\}} $ $ + $ $ T_{4, 0}^{\{m_j, i\}}-2d_{m_j, i}) $, $ \cdots, (T_{1, 0}^{\{m_{|\mathcal{N}_i}|, i\}}+T_{4, 0}^{\{m_{|\mathcal{N}_i|}, i\}} $ $ - $ $ 2d_{m_{|\mathcal{N}_i|}, i})]^{\rm T}-D $和$ \mathit{\boldsymbol{Z}}_{\mathcal{N}_i, k} = [\mathit{\boldsymbol{z}}_{m_1, k}^{\{m_1, i\}} $, $ \cdots, \mathit{\boldsymbol{z}}_{m_j, k}^{\{m_j, i\}} $, $ \cdots $, $ \mathit{\boldsymbol{z}}_{m_{|\mathcal{N}_i|}, k}^{\{m_{|\mathcal{N}_i|}, i\}}]^{\rm T} $.

$ \eta_{i, g, k} $ (参见表 1)满足$ \eta_{i, g, k} = \prod_{m_j = 1}^{|\mathcal{N}_i|}\alpha_{i, m_j, k} $, 其中, $ \alpha_{i, m_j, k} = \begin{cases} \phi_{i, m_j}, &\pi_{i, k}^{\{i, m_j\}} = 1\\ 1-\phi_{i, m_j}, &\pi_{i, k}^{\{i, m_j\}} = 0 \end{cases} $, 则根据式(24)和式(25), 可以分别求出节点状态估计值$ \hat{\mathit{\boldsymbol{x}}}_i(k|k) $和估计误差协方差值E$ [\hat{\mathit{\boldsymbol{x}}}(k|k)] $的期望$ \mathit{\boldsymbol{P}}_i(k|k) $和E$ [\hat{\mathit{\boldsymbol{x}}}_i(k|k)] $. 用更新的期望值作为下一次更新的初始值, 可以得到

代入耦合量测方程(12), 以及辅助绝对量测方程(9)

3.   时钟同步长期运行稳定性统计边界

在绝对时钟建模(通过必要的基础性噪声建模, 以量化时钟同步精度, 方便用于大数据分析)的基础上, 能够得到基本包丢失假设模型下的时钟同步精度分析的保守边界. 更多引入Markov等包丢失模型^[21], 将进一步细化时钟同步精度的边界分析.

我们已经提出与模型具有MMSE量测性能等价的解耦模型(16). 子系统中心节点的MMSE优化性能(通过邻居辅助测量节点)与系统中其他子系统中心节点的MMSE优化性能一致(相邻子系统有重叠). 因此, 本节针对本地子系统中节点的"可靠性边界"分析能够应用于全局系统中的所有节点.

卡尔曼滤波稳定性依赖于误差协方差矩阵$ \mathit{\boldsymbol{P}}_i(k) $ $ = \mathit{\boldsymbol{P}}_i(k|k-1) $收敛于某一有限矩阵. 从前面的讨论, 我们知道量测更新式(27)和式(28)在量测丢失的情形下, 理论上是MMSE性能等价的, 且将式(23)和式(26)代入式(21)

我们知道, 对于没有量测丢失的理想通信网络($ \mathit{\boldsymbol{\pi}}_{i, k} $ $ = \mathit{\boldsymbol{I}} $), 偶对$ (\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{Q}}_i) $的稳定性和偶对$ (\mathit{\boldsymbol{A}}, \bar{\mathit{\boldsymbol{C}}}_i) $的可检测性能够确保$ \mathit{\boldsymbol{P}}_i(k) $从任一初始条件收敛为唯一值^[31]. 在量测丢失的情形, 对于给定的初始条件$ \mathit{\boldsymbol{P}}_i(0) $, 沿着时间传播的误差协方差矩阵$ \{\mathit{\boldsymbol{P}}_i(k)\}_{0}^{\infty} $ (按时间排列的序列)是一个随机过程. 考虑平均意义上的收敛, 也就是, 当$ k\to \infty $时, $ {\rm E}[\mathit{\boldsymbol{P}}_i(k+1)] = {\rm E}[{\rm E}[\mathit{\boldsymbol{P}}_i(k+1)|\mathit{\boldsymbol{P}}_i(k)]]\leq \infty $, 且如文献[15]中那样, 建立相应的收敛条件和稳定域边界. 如文献[15]那样, 我们以如下定理的形式给出收敛条件:

定理3(时钟同步的收敛条件).  对于给定的稳定系统(1)和(16), 也就是矩阵偶对$ (\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{Q}}_i) $和$ (\mathit{\boldsymbol{A}} $, $ \bar{\mathit{\boldsymbol{C}}}_i) $是可控和可观测的, $ \exists 0\leq \{\phi_{c, 1}, \cdots, \phi_{c, |\mathcal{N}_i|}\}\leq $ $ 1 $, 使得$ \lim\nolimits_{k\to \infty}{\rm E}[\mathit{\boldsymbol{P}}_i(k)] = \infty $, 如果所有量测的接收率都在临界值, 除了邻居节点$ j $, 即$ 0\leq \phi_{j}\leq \phi_{c, j} $.

根据定理3, 状态估计过程的稳定性(时钟同步过程)依赖于$ {\rm E}[\mathit{\boldsymbol{P}}_i(k)] $的收敛性, 依次依赖于相对独立的邻居节点的量测的接收率. 由于在多通信链路中随机的量测丢失, 求解误差协方差的稳态边界非常困难. 文献[15]通过对比文献[9]中研究的聚集量测方案, 计算出收敛解的上界和下界, 和对应的临界接收率的上界和下界, 并且得到一个连结的$ k $维域, 这里: 1)在内域的内部, 系统被定义为不稳定; 2)在外域的外部, 系统被定义为稳定; 3)当处于两者之间, 系统是不确定的. 考虑系统(1), 矩阵$ \mathit{\boldsymbol{A}} $的最大特征值是1, 过程噪声$ \mathit{\boldsymbol{Q}}_i $很小, 产生量测的临界接收率下界$ \phi_{c, lb}^{i} = 0 $, 并且在汇聚情景下收敛解$ \mathit{\boldsymbol{P}}_{c, lb}^{i} $相应的下界接近于0. 此外, 对于时钟同步系统, 应该更加关注稳定性, 因此, 我们焦点集中在外域的外部.

量测的临界接收率的上界能够通过以下的优化问题进行求解^[15]:

式中, $ \mathit{\boldsymbol{\Theta}}(\bar{\mathit{\boldsymbol{C}}}_i, \mathit{\boldsymbol{P}}_{c, ub}^{i}, \mathit{\boldsymbol{R}}_{\nu_i}) = \bar{\mathit{\boldsymbol{C}}}_{i}^{\rm T}(\bar{\mathit{\boldsymbol{C}}}_i \pmb{P}_{c, ub}^{i}\bar{\mathit{\boldsymbol{C}}}_{i}^{\rm T}+ \mathit{\boldsymbol{R}}_{\nu_i})^{-1} $ $ \times $ $ \bar{\mathit{\boldsymbol{C}}}_i $; $ \mathit{\boldsymbol{u}}_i(\cdot) $是矩阵$ \mathit{\boldsymbol{P}}_{c, ub}^{i} $的特征向量; $ \mathit{\boldsymbol{S}}_i = [s_i(1) $ $ 0 $; $ 0 $ $ s_i(2)] $是松弛矩阵, 满足$ \mathit{\boldsymbol{S}}_i\geq 0 $; 对于$ g = 1, 2 $, $ \cdots $, $ 2^{|\mathcal{N}_i|} $的概率变量$ \alpha_{i, g} = {\rm E}[\eta_{i, g, k}] $, $ \phi_{c, ub}^i $通过以下线性矩阵不等式问题求解

其中,

误差协方差$ P_{c, ub}^{i} $的上界可以通过求解下面的半正定问题得到

其中, $ \phi = \phi_{c, ub}^i $.

4.   系统状态可靠性与性能评估

4.1   时钟状态估计精度性能对比

1) 事实上, 考虑运行在独立子系统上的并行分布式滤波算法, 从第4.1.1节可知不完全能观的量测模型将导致卡尔曼滤波不稳定.

2) 由附录B的"BMU能观测性证明"、附录C的"MMSE等价变换证明", 能够移除作用于文献[2]中的式(17)不能观量测模型下的强制"K对角化"条件(试图解决卡尔曼滤波不稳定的状态估计).

4.1.1   问题的进一步阐释

1) 本地状态估计的必要条件

本文状态观测器设计目标之一是使观测误差$ \mathit{\boldsymbol{e}}(k) $的极限趋于有限值式(33), 与之等价的命题是:

命题1.  若为"时钟状态追踪" (Clock state tracking)系统, 即通过引入合适的输出反馈, 使"时钟状态追踪"误差系统(34)渐近稳定(即讨论误差系统(34)的镇定问题).

其中, $ \mathit{\boldsymbol{e}}(k|k) = \mathit{\boldsymbol{x}}(k)-\hat{\mathit{\boldsymbol{x}}}(k|k) $. 对于卡尔曼滤波, 使式(34)渐近稳定, 即求解合适的卡尔曼增益$ \mathit{\boldsymbol{K}} $使特征值$ |\lambda(\mathit{\boldsymbol{A}}-\mathit{\boldsymbol{K}}\mathit{\boldsymbol{C}})|<1 $.

在经典控制理论中, 即是输出反馈系统的镇定问题. 本文以式(38)为例, 若给定一组期望的特征值$ \{\lambda_1(\mathit{\boldsymbol{A}}-\mathit{\boldsymbol{KC}}), \cdots, \lambda_4(\mathit{\boldsymbol{A}}-\mathit{\boldsymbol{KC}})\} $($ |\lambda_i(\mathit{\boldsymbol{A}}-\mathit{\boldsymbol{KC}})|<1 $, 则存在对应的特征方程

其中, $ a_0, \cdots, a_3 $是期望特征方程对应的系数. 令$ \bar{\mathit{\boldsymbol{A}}} $ $ = $ $ \mathit{\boldsymbol{A}}-\mathit{\boldsymbol{KC}} $, 根据Cayley-Hamilton定理, 有矩阵方程$ d(\bar{\mathit{\boldsymbol{A}}}) = 0 $.

得到Ackermann公式

$ d(\mathit{\boldsymbol{A}}) $是一个给定常数, 若想通过式(37) 求解得到$ \mathit{\boldsymbol{K}} $, 需保证$ [\mathit{\boldsymbol{CA}}^3\; \; \mathit{\boldsymbol{CA}}^2\; \; \mathit{\boldsymbol{CA}}\; \; \mathit{\boldsymbol{C}}]^{\rm T} $存在Moore-Penrose逆. 因此, $ [ \mathit{\boldsymbol{CA}}^3\; \; \mathit{\boldsymbol{CA}}^2\; \; \mathit{\boldsymbol{CA}}\; \; \mathit{\boldsymbol{C}}]^{\rm T} $一定是列满秩矩阵. 对于系统(39), 通过"观测器状态空间极点配置"的讨论^[32-33], 仍然等价于须式(39)在BMU下式(38) $ (\mathit{\boldsymbol{A}}, \bar{\mathit{\boldsymbol{C}}}_i^{\{i, j\}}) $能观. 对于卡尔曼滤波在能观性条件下的唯一稳态解的讨论, 可参考文献[34-35] (已超出本文的讨论范围.)

所以, 时钟参数确定性估计的必要条件是保证系统的能观测性(参见附录B).

4.1.2   性能评估模型

1) BMU的能观测性

a) 事实上, $ S_i $和$ S_j $间的相对状态能够由$ \mathit{\boldsymbol{z}}_{i}^{\{i, j\}}(k) $ (或$ \mathit{\boldsymbol{z}}_j^{\{j, i\}}(k) $)测量(文献[2]中的式(17)). 但$ S_j $和$ S_i $的状态不可观(证明见附录A), 仅能确定$ \mathit{\boldsymbol{x}}_i(k) $和$ \mathit{\boldsymbol{x}}_j(k) $相对量测值.

b) 结合式(1)与式(16), 当系统在完全量测下时(不发生丢包), 能够得到关于节点$ S_i $的BMU能观测性量测模型(参见附录B)为

BMU的能观测性证明详见附录B.

2) 子系统的能观测性(BMU的线性组合)

证明详见附录B.1. 性能评估请参见第4.1.3节.

4.1.3   仿真与性能评估

实验参数设置$ \mathit{\boldsymbol{A}}_i = \begin{bmatrix} 1&0\\ 0.1&1 \end{bmatrix} $, $ \mathit{\boldsymbol{A}} = {\rm diag}\{\mathit{\boldsymbol{A}}_i $, $ \mathit{\boldsymbol{A}}_j\} $. $ \mathit{\boldsymbol{A}}_j $与$ \mathit{\boldsymbol{A}}_i $相同. $ \mathit{\boldsymbol{w}}_i(k) = \begin{bmatrix} 1\\ 0.1 \end{bmatrix}\mathit{\boldsymbol{u}}_i(k) $, $ \sigma_{R_P}^{2} = 0.1 $. $ \mathit{\boldsymbol{Q}}_i = \mathit{\boldsymbol{Q}}_j = \begin{bmatrix} 2.7\times 10^{-15}&2.7\times 10^{-16}\\ 2.7\times 10^{-16}&2.7\times 10^{-17} \end{bmatrix} $, $ \mathit{\boldsymbol{Q}} = {\rm diag}\{\mathit{\boldsymbol{Q}}_i, \mathit{\boldsymbol{Q}}_j\} $. $ \mathit{\boldsymbol{P}}(0|0) = {\rm diag}(\mathit{\boldsymbol{P}}_i(0|0), \mathit{\boldsymbol{P}}_j(0|0)) $, $ \mathit{\boldsymbol{P}}_i(0|0) = 100\pmb{I}_{2\times 2} $, $ \hat{\mathit{\boldsymbol{x}}}_i(0|0) = \hat{\mathit{\boldsymbol{x}}}_j(0|0) = \begin{bmatrix} 10\\ 10 \end{bmatrix} $, $ \hat{\mathit{\boldsymbol{x}}}(0|0) $ $ = \begin{bmatrix} \hat{\mathit{\boldsymbol{x}}}_i(0|0)\\ \hat{\mathit{\boldsymbol{x}}}_j(0|0) \end{bmatrix} $.

定义根平方均方误差(Root average mean square error, RAMSE)为

其中, $ \mathit{\boldsymbol{\varsigma}}\in\{\beta, \vartheta\} $, $ N $是节点数量, $ M $为试验次数.

定义根平方误差(Root mean square error, RMSE)为

其中, $ \mathit{\boldsymbol{\varsigma}}\in\{\beta, \vartheta\} $.

1) 能观测性对卡尔曼滤波的影响

参考方程(1)与方程(12), 系统的相对量测模型为

本实例BMU的方程定义详见附录A中式(A3).

根据附录B, 我们认识到若想确定性估计可控的时钟参数, 必须保证系统的能观测性. 式(39)中, 在BMU下, 能观测性矩阵的秩是2, 意味着式(39)是一个不完全能观测的系统(附录A). 本实例比较基于式(39)和基于式(38)的RMSE和RAMSE, 以此评估系统的能观测性对卡尔曼滤波稳定性的影响.

由图 2和图 3可以看出, 系统能观测性的缺失造成了其对应的卡尔曼滤波不稳定, 即$ \beta(k) $和$ \hat{\beta}(k) $之间的偏差, 即$ RAMSE(\beta(k)) $, 随着步数的增加而增大. 并且$ RMSE(\beta(k)) $ (或$ RMSE(\vartheta(k)) $)具有同样的变化趋势. 而图中基于式(38)的卡尔曼滤波则可以正常完成滤波.

图 2  可观性对BMU RAMSE系统稳定性的影响

Fig. 2  Effect of observability for system stability of a BMU-RAMSE

下载: 全尺寸图片 幻灯片

图 3  可观性对BMU-RMSE系统稳定性的影响

Fig. 3  Effect of observability for system stability of a BMU-RMSE

下载: 全尺寸图片 幻灯片

2) 解耦量测模型的可扩展性

解耦量测模型具有MMSE性能的可扩展性. 本实例将比较不同规模的子系统的滤波性能(本质为以BMU为单元的多量测链路子系统). 根据式(38), 以$ S_i $为中心, 有$ \mathcal{N}_i $个邻居节点的子系统作如下扩展

其中, $ \mathit{\boldsymbol{y}}_{i, k} = [\mathit{\boldsymbol{\gamma}}_{m_1, k}^{\{m_1, i\}} - $ $ \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_1\}} $, $ \cdots, \mathit{\boldsymbol{\gamma}}_{m_j, k} ^{\{m_j, i\}}-\mathit{\boldsymbol{z}}_{i, k}^{\{i, m_j\}}, \cdots, $ $ \mathit{\boldsymbol{\gamma}}_{m_{|N_i|}, k}^{\{m_{|N_i|}, i\}}- \mathit{\boldsymbol{z}}_{i, k}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $, $ \mathit{\boldsymbol{y}}_{i, k}\in {\bf R}^{|N_i|} $, $ \bar{\mathit{\boldsymbol{C}}}_i $ $ = $ $ \begin{bmatrix} 0&2\\ 0&2\\ \vdots&\vdots\\ \end{bmatrix}\in {\bf R}^{|\mathcal{N}_i|\times 2} $. 本实例中试验次数$ M $ $ = $ $ 80 $, 仿真结果如图 4和图 5所示.

图 4  BMU系统RAMSE估计性能

Fig. 4  Estimation performance with BMU systems-RAMSE

下载: 全尺寸图片 幻灯片

图 5  BMU系统RMSE估计性能

Fig. 5  Estimation performance with BMU systems-RMSE

下载: 全尺寸图片 幻灯片

通过观察图 4和图 5, 我们能够发现, 随着邻居节点个数的增多, 子系统在$ k $时刻可以获得更多的量测值. 卡尔曼滤波的估计性能也会因此由子系统多链路量测"过采样"而得到改善. 从图中可以看出, 有10个邻居节点的子系统的RAMSE下降最快并且也最先达到稳态.

4.2   不可靠工业物联网状态追踪性能评估

通过定理1、定理2及MMSE等价变换分析, 子系统中心节点的MMSE优化性能与其他子系统中心节点的MMSE性能一致. 不失一般性, 本节中本地子系统时钟同步状态的可靠性性能评估是规模可扩展的.

在仿真中, 选择4个节点, 如图 6所示. 采样周期$ \tau_0 = 0.1s $且$ \Delta = 1 $, 且$ \sigma_{q}^2 = 2.7\times 10^{-10} $. 以节点$ S_1 $为例, 邻居节点是$ S_2 $, $ S_3 $和$ S_4 $, 相应的量测方差为: $ \mathit{\boldsymbol{R}}_{\zeta_1} = \mathit{\boldsymbol{R}}_{\upsilon_i} = {\rm diag}\{0.25, 0.125, 0.5\} $, $ \mathit{\boldsymbol{R}}_{\nu_1} = $ $ {\rm diag}\{0.5, 0.25, 1\} $.

图 6  测量的临界接受率—上限

Fig. 6  Critical acceptance rate of measurement — upper bound

下载: 全尺寸图片 幻灯片

考虑汇聚的情景, 通过求解能够获得关于$ \phi_c^1 $的上界, 且计算结果为$ \phi_{c, ub}^1 = 0.01 $. 通过求解, 能够计算相应的收敛的误差协方差的上界, 且通过求解, 能够计算每个邻居节点的临界接收率的边界.

对于$ \phi_{c, ub}^1>0.01\to 1 $, 重复以上的步骤, 对于节点$ S_1 $, 我们能够计算在不同的潜在通信链路质量下对应的时钟同步准确度的均衡值. 同样, 节点$ S_2 $, $ S_3 $和$ S_4 $也能够通过求解和获得均衡值.

图 6显示从每一个邻居节点的量测临界接收率与收敛的误差协方差矩阵的迹对应的上界. 对于时钟同步系统中, 依赖于估计的状态所需要的准确度, 容易通过图 6推断出相应的网络设施质量的需求.

例如, 要求tr($ \mathit{\boldsymbol{P}}_1 $)的稳态值为0.001, 则对于$ \phi_{1, c, ub} >0.0855 $, $ \phi_{2, c, ub}>0.041 $, $ \phi_{3, c, ub}>0.2155 $能够定义为稳定. 我们做1 000次Monte Carlo实验(如图 7所示). 图中显示在丢包情形下(虚线)和没有丢包的情形下(实线), 迭代次数与误差协方差矩阵的迹的平均的关系, 且经过验证, 经过有限时间后, 在统计意义下, 误差协方差矩阵的迹(如图 7所示)没有超过相应的上界.

图 7  Monte Carlo实验

Fig. 7  Monte Carlo experiment

下载: 全尺寸图片 幻灯片

针对本文的实验验证与进一步研究, 相关的研究中(参见一致性算法^{[5, 7-8]}和文献[36]) 已经存在全面的阐述.对于复杂的大数据工程的适应性验证非常重要, 仍需结合具体的工程实现细节进行分析.

5.   结论

时钟同步系统的规划和设计中, 重要的方面是同步性能的权衡. 我们已经提出能观测性MMSE等价变换的状态解耦量测方程及全分布式修改卡尔曼滤波(适用于网络规模化扩展). 通过理论分析时钟同步的条件和计算统计同步误差相应的上界, 当前论文在时钟同步准确性与潜在的通信链路质量间已经作了量化均衡. 针对规模化工业物联网的硬实时调度任务的"紧时隙", 如"工业物联网确定性调度中TDMA紧时隙时间精度边界可靠性分析"具有重要的研究意义.

本文已经提供了包丢失假设下的可扩展研究的基本框架与线索. 进一步, 我们能够追踪工业无线网络Cyber-Physical条件下状态估计的前沿研究, 以全面考察规模化精确时钟同步的状态估计精度. 针对时钟同步具体应用领域, 本文已经提供了必需的(理论分析的) 研究参考价值; 进一步扩展研究的框架线索(如Markov模型^[21]等); 基于综合优化角度, 在更全面完整包丢失信道模型假设下(如文献[28]等), 则能在本文卡尔曼滤波建模基础上扩展时钟同步无线频谱感知的状态估计策略.

附录A.   问题的提出

文献[2]从一对节点出发, 利用双向信息交换机制建立系统量测模型, 结合完全解耦的节点的时钟状态转移方程(1), 得到一个子系统的时钟状态空间模型(文献[2]中的式(17). 我们以一对节点$ S_i $和$ S_j $ (BMU)为例, 分析基本量测单元的能观测性. 对于$ S_i $和$ S_j $节点分别有下面的系统

其中, 对于$ m = i, j $, $ \mathit{\boldsymbol{A}}_m = \begin{bmatrix} 1&0\\ \Delta \tau_0&1 \end{bmatrix} $, 时钟状态$ \mathit{\boldsymbol{x}}_m(l_k) $ $ = \begin{bmatrix} \beta_m(l_k)\\ \theta_m(l_k) \end{bmatrix} $, $ \mathit{\boldsymbol{w}}_m(l_k) $和$ \mathit{\boldsymbol{\nu}}_m(l_k) $分别为过程噪声和量测噪声, $ \mathit{\boldsymbol{b}}_m = \begin{bmatrix} 0\\ -\Delta \tau_0 \end{bmatrix} $. 结合式(A1)和(A2), 可以得到关于BMU系统的状态空间模型

其中, 子系统状态$ \mathit{\boldsymbol{x}}(k) = \begin{bmatrix} \mathit{\boldsymbol{x}}_i(l_k)\\ \mathit{\boldsymbol{x}}_j(l_k) \end{bmatrix} $, 状态转移矩阵$ \mathit{\boldsymbol{A}} = \begin{bmatrix} \mathit{\boldsymbol{A}}_i&0\\ 0&\mathit{\boldsymbol{A}}_j \end{bmatrix} $, 常数项$ \mathit{\boldsymbol{b}} = \begin{bmatrix} \mathit{\boldsymbol{b}}_i\\ \mathit{\boldsymbol{b}}_j \end{bmatrix} $, 过程噪声$ \mathit{\boldsymbol{w}}(k) $ $ = $ $ \begin{bmatrix} \mathit{\boldsymbol{w}}_i(l_k)\\ \mathit{\boldsymbol{w}}_j(l_k) \end{bmatrix} $, 观测值$ \mathit{\boldsymbol{Z}}(k) = \begin{bmatrix} \mathit{\boldsymbol{z}}_{i, l_k}\\ \mathit{\boldsymbol{z}}_{j, l_k} \end{bmatrix} $, 观测矩阵$ \mathit{\boldsymbol{C}} = \left[{\begin{array}{rrrr} 0&-2&0&2\\ 0&2&0&-2 \end{array}}\right] $, 观测噪声$ \mathit{\boldsymbol{\nu}}(k) = $ $ \begin{bmatrix} \mathit{\boldsymbol{\nu}}_i(l_k)\\ \mathit{\boldsymbol{\nu_j}}(l_k) \end{bmatrix} $. 定义: 对于系统(A3), 若能根据第$ k $步及以后若干步的观测$ \mathit{\boldsymbol{Z}}(k) $, $ \mathit{\boldsymbol{Z}}(k+1), \cdots $, $ \mathit{\boldsymbol{Z}}(k+N- $ $ 1) $, 就能唯一地确定出第$ k $步的每一状态$ \mathit{\boldsymbol{x}}(k) $, 则称系统在第$ k $步是能观测的. 若系统在任何一步上都能观测, 则称系统是完全能观测的, 或简称能观测的.

从系统(A3)的观测方程可以看出, 第$ k $步的观测矩阵的秩为1($ \mathit{\boldsymbol{x}}(k) $的维数为4), 那么由$ \mathit{\boldsymbol{Z}}(k) $不能把$ \mathit{\boldsymbol{x}}(k) $唯一地确定出来. 考虑再加3步观测$ \mathit{\boldsymbol{Z}}(k+1) $, $ \mathit{\boldsymbol{Z}}(k+2) $, $ \mathit{\boldsymbol{Z}}(k+3) $, 根据系统的观测方程和状态方程不难得出

$ \mathit{\boldsymbol{W}}_o(N) = \begin{bmatrix} \mathit{\boldsymbol{C}}&\mathit{\boldsymbol{CA}}&\cdots&\mathit{\boldsymbol{CA}}^{N-1} \end{bmatrix}^{\rm T} $, 则

可以很容易知道, 矩阵$ \mathit{\boldsymbol{W}}_o(4) $秩为2, 由方程(A4)知, 通过$ \mathit{\boldsymbol{Z}}(k) $, $ \mathit{\boldsymbol{Z}}(k+1) $, $ \mathit{\boldsymbol{Z}}(k+2) $, $ \mathit{\boldsymbol{Z}}(k+3) $不能把$ \mathit{\boldsymbol{x}}(k) $唯一地确定出来(不存在最小二乘解), 即使再增加更多的观测值, 矩阵$ \mathit{\boldsymbol{W}}_o(N)\; (N>4) $的秩仍然为2, 即还是不能把$ \mathit{\boldsymbol{x}}(k) $唯一地确定出来(只能唯一确定偏斜和偏移的相对量(最小二乘解)). 因此基本量测单元是不可观的.

对于由多组基本量测单元构成的子系统, 由于基本量测单元之间的相互独立, 类似一组BMU的情况, 得子系统的能观测性矩阵$ \mathit{\boldsymbol{W}}_o(N) $ $ (N>2|\mathcal{N}_{[i]}|) $的秩为$ 2|\mathcal{N}_{[i]}|-2 $ (($ \mathit{\boldsymbol{x}}(k) $的维数为$ 2|\mathcal{N}_{[i]}| $), 能够唯一确定所有邻居节点$ S_j $ $ (j\in \mathcal{N}_i) $与中心节点$ S_i $的相对偏移量$ \vartheta_j(k)-\vartheta_i(k) $和相对偏斜量$ \beta_j(k)- \beta_i(k) $. 不能把$ \mathit{\boldsymbol{x}}(k) $唯一地确定出来(不存在最小二乘解), 所以子系统是不可观的.

附录B.   BMU能观测性证明

通过上述研究得知系统(A1)是不可观的, 使用多组量测值仅能确定中心节点$ S_i $与邻居节点$ S_{j} $之间的时钟状态相对量, 为能确定$ S_i $的状态, 通过分析双向信息交换机制, 对每一个节点建立一个自身状态量测模型

其中, $ \mathit{\boldsymbol{x}}(k) = \begin{bmatrix} \mathit{\boldsymbol{x}}_i(k)\\ \mathit{\boldsymbol{x}}_j(k) \end{bmatrix} $, $ \mathit{\boldsymbol{M}}_i^{\{i, j\}} = [ 0\; \; 2\; \; 0\; \; 0 ] $和$ \mathit{\boldsymbol{M}}_j^{\{j, i\}} $ $ = [0\; \; 0\; \; 0\; \; 2] $; $ \mathit{\boldsymbol{\zeta}}_i^{\{i, j\}}(k) $和$ \mathit{\boldsymbol{\zeta}}_j^{\{j, i\}}(k) $为相互独立的量测噪声. 以$ S_i $为例, 借助辅助量测(B2), 对系统(A1)中的量测进行的MMSE等价变换可以描述为

即, 令$ \mathit{\boldsymbol{y}}_{i, k}^{\{i, j\}} = \mathit{\boldsymbol{\gamma}}_{j, k}^{\{j, i\}}-z_i^{\{i, j\}}(k) $, $ \mathit{\boldsymbol{C}}_i^{\{i, j\}} = \mathit{\boldsymbol{M}}_j^{\{j, i\}}- $ $ \mathit{\boldsymbol{H}}_i^{\{i, j\}} = [ 0\; \; 2\; \; 0\; \; 0 ] $, $ \mathit{\boldsymbol{\nu}}_i^{\{i, j\}}(k) = \mathit{\boldsymbol{\zeta}}_j^{\{j, i\}}(k)-\mathit{\boldsymbol{\upsilon}}_i^{\{i, j\}}(k) $, 其中, $ \mathit{\boldsymbol{H}}_i^{\{i, j\}} = [0\; \; -2\; \; 0\; \; 2 ] $, 则得$ \mathit{\boldsymbol{y}}_{i, k}^{\{i, j\}} = \mathit{\boldsymbol{C}}_{i}^{\{i, j\}}\mathit{\boldsymbol{x}}(k) $ $ + $ $ \mathit{\boldsymbol{\nu}}_i^{\{i, j\}}(k) $. 方程(B3)可以描述成$ S_i $依靠量测节点$ S_j $达到观测自身状态. 下面通过能观测加以证明(说明), 系统可以简化成

其中, $ \bar{\mathit{\boldsymbol{C}}}_i^{\{i, j\}} = [0\; \; 2] $.

从系统(B4)的观测方程可以看出, 第$ k $步的观测矩阵的秩为1 ($ \mathit{\boldsymbol{x}}_i(k) $中包含$ \beta_i(k) $和$ \vartheta_i(k) $), 那么由$ \mathit{\boldsymbol{y}}_{i, k}^{\{i, j\}} $不能把$ \mathit{\boldsymbol{x}}_i(k) $唯一地确定出来. 考虑再加1步的观测$ \mathit{\boldsymbol{y}}_{i, k+1}^{\{i, j\}} $, 根据系统(B4)的观测方程和状态方程不难得出

令$ \bar{\mathit{\boldsymbol{W}}}_{i, o}^{\{i, j\}}(N) = \begin{bmatrix} \bar{\mathit{\boldsymbol{C}}}_i^{\{i, j\}}\\ \bar{\mathit{\boldsymbol{C}}}_i^{\{i, j\}}\mathit{\boldsymbol{A}}\\ \vdots\\ \bar{\mathit{\boldsymbol{C}}}_i^{\{i, j\}}\mathit{\boldsymbol{A}}^{N-1} \end{bmatrix} $, 计算$ \bar{\mathit{\boldsymbol{W}}}_{i, o}^{\{i, j\}}(2) $ $ = \begin{bmatrix} 0&2\\ 2\Delta \tau_0\times 1&2 \end{bmatrix} $, 可以很容易知道矩阵$ \bar{\mathit{\boldsymbol{W}}}_{i, o}^{\{i, j\}}(2) $秩为2, 由方程(B5)知, 通过$ \mathit{\boldsymbol{y}}_{i, k}^{\{i, j\}} $, $ \mathit{\boldsymbol{y}}_{i, k+1}^{\{i, j\}} $能把$ \mathit{\boldsymbol{x}}_i(k) $唯一地确定出来, 若再增加更多的观测值, 矩阵$ \bar{\mathit{\boldsymbol{W}}}_{i, o}^{\{i, j\}}(N) $ $ (N>2) $的秩仍然为2, 仍然能把$ \mathit{\boldsymbol{x}}_i(k) $唯一地确定出来(此时为最小二乘解). 类似地, 通过$ \mathit{\boldsymbol{y}}_{j, k}^{\{j, i\}} = \mathit{\boldsymbol{\gamma}}_{i, k}^{\{i, j\}}- \mathit{\boldsymbol{z}}_j^{\{j, i\}}(k) $和$ \mathit{\boldsymbol{y}}_{j, k+1}^{\{j, i\}} $就能把$ \mathit{\boldsymbol{x}}_j(k) $唯一地确定出来. 因此在该量测方程下BMU是能观测性的.

B.1   由BMU组合扩展的系统层面的能观测性证明

对于无线传感器网络(WSN) $ G(V, E) $, 可以分成$ |E| $组独立的BMU, 其中$ |E| $代表无向图$ G $中所有边数. $ S_i $ $ (i = \{1, 2, \cdots, N\}) $为中心及其所有周围邻居节点组成的子系统, 分成$ |\mathcal{N}_i| $组BMU, $ |\mathcal{N}_i| $为$ S_i $的邻居节点数. $ BMU(S_i) = \{BMU(S_i, S_j)| $ $ (S_i, S_j)\in E, $且$ S_i, S_j\in V\} $, 其中, $ BMU(S_i, S_j) $表示由$ S_i $和$ S_j $构成的BMU, 借助辅助量测$ \mathit{\boldsymbol{\gamma}}_{i, k} $ $ = $ $ \mathit{\boldsymbol{M}}_i\pmb{x}(k)+\mathit{\boldsymbol{\zeta}}_i(k) $, 进行MMSE等价变换(如式(13)$ \sim $(15)所示), 可以得到解耦的量测方程$ \mathit{\boldsymbol{y}}_{i, k} $ $ = $ $ \bar{\mathit{\boldsymbol{C}}}_i \pmb{x}_i(k)+\mathit{\boldsymbol{\nu}}_i(k) $ (如式(16), 不含丢包时), 类似BMU的能观测证明过程, 可以得到能观测性矩阵

易知rank$ (\bar{\mathit{\boldsymbol{W}}}_{i, o}(2)) = 2 $, 说明子系统也是能观测性的.

B.2   详细的算法过程(量测消息传输、量测更新)

正是因为在实际中时钟状态的直接量测难以获得, 所以有必要通过一个可用量测模型和一个时钟状态转移模型确定时钟的状态.

针对BMU, 通过式(B3), 我们构造出一个具有能观测性的量测方程, 并利用作为时钟状态的直接量测$ \bar{\mathit{\boldsymbol{\gamma}}}_{j, k}^{\{j, i\}} $客观存在的事实, 从理论上说明(证明) MMSE等价变换的可行性.

为在实际应用中实现式(B3)的量测步骤(特指MMSE等价量测), 由替换观测量$ \bar{\mathit{\boldsymbol{\gamma}}}_{j, k}^{\{j, i\}} $来代替$ \mathit{\boldsymbol{\gamma}}_{j, k}^{\{j, i\}} $, $ \bar{\mathit{\boldsymbol{\gamma}}}_{j, k}^{\{j, i\}} = \bar{\mathit{\boldsymbol{\gamma}}}_{j, k-1}^{\{j, i\}}+\mathit{\boldsymbol{z}}_{j, k-1}^{\{j, i\}}-\mathit{\boldsymbol{z}}_{j, k}^{\{j, i\}} $和$ \bar{\mathit{\boldsymbol{\gamma}}}_{j, 0}^{\{j, i\}} = $ $ \mathit{\boldsymbol{\gamma}}_{j, 0}^{\{j, i\}} $ $ = T_{1, 0}^{\{j, i\}}+T_{4, 0}^{\{j, i\}}-2d_{ij}-D $. 在实际双向信息交换过程中, 对于$ S_i $而言, 我们知道$ \mathit{\boldsymbol{z}}_i^{\{i, j\}}(k) $在本地, 而$ \mathit{\boldsymbol{z}}_{j, k}^{\{j, i\}} $和$ \mathit{\boldsymbol{z}}_{j, k-1}^{\{j.i\}} $位于$ S_j $中, 所以在$ S_j $上更新$ \bar{\mathit{\boldsymbol{\gamma}}}_{j, k}^{\{j, i\}} $, 然后向$ S_i $发送量测信息$ \bar{\mathit{\boldsymbol{\gamma}}}_{j, k}^{\{j, i\}} $.

附录C.   可靠量测下MMSE等价变换

整个网络由多个相互独立的BMU组成, 我们以一个BMU来阐述这一过程. 以BMU ($ S_i-S_j $)为例, 通过双向信息交换机制, 可以得到下面的量测模型

其中, $ \mathit{\boldsymbol{H}}_i = [0\; \; -2\; \; 0\; \; 2 ] $和$ \mathit{\boldsymbol{H}}_j = [ 0\; \; 2\; \; 0\; \; -2 ] $; $ \mathit{\boldsymbol{x}}(k) $ $ = \begin{bmatrix} \mathit{\boldsymbol{x}}_i(k)\\ \mathit{\boldsymbol{x}}_j(k) \end{bmatrix} $, 对于$ m = i, j $, 时钟状态$ \mathit{\boldsymbol{x}}_m(k) = \begin{bmatrix} \beta_m(k)\\ \vartheta_m(k) \end{bmatrix} $; 量测噪声$ \mathit{\boldsymbol{\nu}}_i(k) $和$ \mathit{\boldsymbol{\nu}}_j(k) $是相互独立的高斯噪声.

通过能观测性分析, 知道在量测模型(C1)和(C2)下的BMU系统是不可观的. 为解决BMU系统的不可观问题, 以达到通过观测值确定系统状态, 引入下面的辅助量测模型(绝对量测模型)

其中, $ \mathit{\boldsymbol{M}}_i = [0\; \; 2\; \; 0\; \; 0 ] $和$ \mathit{\boldsymbol{M}}_j = [ 0\; \; 0\; \; 0\; \; 2 ] $; 量测噪声$ \mathit{\boldsymbol{\zeta}}_i(k) $和$ \mathit{\boldsymbol{\zeta}}_j(k) $为相互独立的高斯噪声.

借助辅助量测(C3)和(C4), 对量测模型(C1)和(C2)进行的MMSE等价变换可以描述为

其中, 令$ \mathit{\boldsymbol{y}}_{i, k} = \mathit{\boldsymbol{\gamma}}_{j, k}-\mathit{\boldsymbol{z}}_{i, k} $和$ \mathit{\boldsymbol{C}}_i = \mathit{\boldsymbol{M}}_j -\mathit{\boldsymbol{H}}_i = [0 $ 2 0 $ 0 ] $; $ \mathit{\boldsymbol{y}}_{j, k} = \mathit{\boldsymbol{\gamma}}_{i, k} -\mathit{\boldsymbol{z}}_{j, k} $和$ \mathit{\boldsymbol{C}}_j = \mathit{\boldsymbol{M}}_i-\mathit{\boldsymbol{H}}_j = [ 0\; 0\; 0\; $$ 2 ] $; $ \mathit{\boldsymbol{\nu}}_i(k) = \mathit{\boldsymbol{\zeta}}_j(k)\mathit{\boldsymbol{-\upsilon}}_i(k) $和$ \mathit{\boldsymbol{\nu}}_j(k) = \mathit{\boldsymbol{\zeta}}_i(k)-\mathit{\boldsymbol{\upsilon}}_j(k) $. 方程(C5) (方程(C6))可以描述成$ S_i $ ($ S_j $)依靠量测节点$ S_j $ ($ S_i $)达到观测自身状态.

C.1   证明过程(拆分的观点)

因为变换(C5)和(C6)具有对称性, 只需证明其中一个即可. 以$ S_i $为例, 为证明式(C5)的变换过程是MMSE等价变换, 我们以具有MMSE量测性能的卡尔曼滤波为观测器, 证明时钟状态$ \mathit{\boldsymbol{x}}_i(k) $在变换之前的量测模型(C1)和(C4)下的估计误差协方差矩阵与在变换之后的量测模型(C5)下的估计误差协方差矩阵趋向一致(具有相同的收敛值).

BMU的时钟状态转移方程为

其中, 状态转移矩阵$ \mathit{\boldsymbol{A}} = \begin{bmatrix} \mathit{\boldsymbol{A}}_i&0\\ 0&\mathit{\boldsymbol{A}}_j \end{bmatrix} $, 常数项$ \mathit{\boldsymbol{b}} $ $ = $ $ \begin{bmatrix} \mathit{\boldsymbol{b}}_i\\ \mathit{\boldsymbol{b}}_j \end{bmatrix} $, 过程噪声$ \mathit{\boldsymbol{w}}(k) = \begin{bmatrix} \mathit{\boldsymbol{w}}_i(k)\\ \mathit{\boldsymbol{w}}_j(k) \end{bmatrix} $是均值为0和协方差为$ \mathit{\boldsymbol{Q}} = \begin{bmatrix} \mathit{\boldsymbol{Q}}_i&0\\ 0&\mathit{\boldsymbol{Q}}_j \end{bmatrix} $的高斯噪声.

组合量测方程(C1)和(C4), 得变换之前的量测模型为

其中, 量测值$ \mathit{\boldsymbol{Z}}_{i, j}(k) = \begin{bmatrix} \mathit{\boldsymbol{\gamma}}_{j, k}\\ \mathit{\boldsymbol{z}}_{i, k} \end{bmatrix} $, 量测矩阵$ \mathit{\boldsymbol{G}}_{i, j} = \begin{bmatrix} \mathit{\boldsymbol{M}}_j\\ \mathit{\boldsymbol{H}}_i \end{bmatrix} = \begin{bmatrix} 0&0&0&2\\ 0&-2&0&2 \end{bmatrix} $, 量测噪声$ \mathit{\boldsymbol{\varpi}}_{i, j}(k) = \begin{bmatrix} \mathit{\boldsymbol{\zeta}}_j(k)\\ \mathit{\boldsymbol{\upsilon}}_i(k) \end{bmatrix} $是均值为0和协方差为$ \mathit{\boldsymbol{R}}_{\varpi_{i, j}}\in {\bf R}^{2\times 2} $的高斯噪声.

由方程(C7)和(C8)组成的变换前BMU系统在标准卡尔曼滤波下的验前估计误差协方差矩阵(State error covariance prediction) $ \mathit{\boldsymbol{S}}_{i, k}^{\rm pre} $满足

因为变换之后的量测模型(C5)是状态解耦的, 对于$ S_j $的状态是不可观的, 而为了确保BMU系统的能观测性以及不影响$ S_i $节点状态的估计误差协方差, 需要(再次)借助辅助量测(C4) (算法过程解释为需要由辅助量测输出(C4), 见附录B.2), 通过消息传递参与$ S_i $节点的状态估计. 组合量测方程(C4)和(C5)得变换之后的量测模型为

其中, $ \mathit{\boldsymbol{Y}}_{i, j}(k) = \begin{bmatrix} \mathit{\boldsymbol{\gamma}}_{j, k}\\ \mathit{\boldsymbol{\gamma}}_{j, k}-\mathit{\boldsymbol{z}}_{i, k} \end{bmatrix} $是量测值, 量测矩阵$ \mathit{\boldsymbol{\Gamma}}_{i, j} $ $ = $ $ \begin{bmatrix} \mathit{\boldsymbol{M}}_j\\ \mathit{\boldsymbol{M}}_j-\mathit{\boldsymbol{H}}_i \end{bmatrix} = \begin{bmatrix} 0&0&0&2\\ 0&2&0&0 \end{bmatrix} $, $ \mathit{\boldsymbol{\psi}}_{i, j}(k) = \begin{bmatrix} \mathit{\boldsymbol{\zeta}}_j(k)\\ \mathit{\boldsymbol{\zeta}}_j(k)-\mathit{\boldsymbol{\nu}}_i(k) \end{bmatrix} $是量测噪声, 均值为0和协方差为$ \mathit{\boldsymbol{R}}_{\psi_{i, j}} $ $ \in $ $ {\bf R}^{2\times 2} $的高斯噪声. 很容易得到一个可逆矩阵$ \mathit{\boldsymbol{P}}_{i, j} = \begin{bmatrix} 1&0\\ 1&-1 \end{bmatrix} $, 使得$ \mathit{\boldsymbol{\Gamma}}_{i, j} = \mathit{\boldsymbol{P}}_{i, j}\mathit{\boldsymbol{G}}_{i, j} $和$ \mathit{\boldsymbol{\psi}}_{i, j}(k) $ $ = $ $ \mathit{\boldsymbol{P}}_{i, j}\mathit{\boldsymbol{\varpi}}_{i, j}(k) $以及$ \mathit{\boldsymbol{R}}_{\psi_{i, j}} = \mathit{\boldsymbol{P}}_{i, j}\mathit{\boldsymbol{R}}_{\varpi_{i, j}}\mathit{\boldsymbol{P}}_{i, j}^{{\rm T}} $成立.

由方程(C7)和(C10)组成的变换前BMU系统在标准卡尔曼滤波下的验前估计误差协方差矩阵(State error covariance prediction)满足

将$ \mathit{\boldsymbol{\Gamma}}_{i, j} = \mathit{\boldsymbol{P}}_{i, j}\mathit{\boldsymbol{G}}_{i, j} $和$ \mathit{\boldsymbol{R}}_{\psi_{i, j}}{ = \mathit{\boldsymbol{P}}_{i, j}\mathit{\boldsymbol{R}}_{\mathit{\boldsymbol{\varpi}}_{i, j}}}\mathit{\boldsymbol{P}}_{i, j}^{\rm T} $代入方程(C11), 并利用可逆矩阵的性质$ \mathit{\boldsymbol{P}}_{i, j}\mathit{\boldsymbol{P}}_{i, j}^{-1} = \mathit{\boldsymbol{I}} $, 方程(C11)可以重写为

对比迭代式(C9)和(C12), 知(C9)和(C12)有相同的收敛解(卡尔曼滤波的协方差矩阵收敛值与初始值无关, 它是唯一的), 甚至, 当设置相同的初始值时, 即$ \mathit{\boldsymbol{S}}_{i, 0}^{\rm post} = \mathit{\boldsymbol{S}}_{i, 0}^{\rm pre} $, 可以得到$ \mathit{\boldsymbol{S}}_{i, k}^{\rm post} = \mathit{\boldsymbol{S}}_{i, k}^{\rm pre} $, 对于$ \forall k\in {\bf N} $. 所以量测方程(C10)和(C8)在MMSE量测性能上是等价的, 也显示了量测方程(C1)在辅助量测方程(C4)的帮助下存在MMSE量测性能解耦. 由此容易得到时钟状态$ \mathit{\boldsymbol{x}}_i(k) $在变换之前的量测模型(C1)和(C4)下的估计误差协方差矩阵与在变换之后的量测模型(C5)下的估计误差协方差矩阵趋向一致(具有相同的收敛值).

附录D.   多链路量测的MMSE等价变换证明(拆分)

D.1   绝对量测与相对量测中丢包系数矩阵定义

将不可靠网络下的绝对量测方程(9)和相对量测方程(12)展开可表示为

式中, $ \tilde{\mathit{\boldsymbol{M}}}_i^{\{i, m_j\}} $和$ \tilde{\mathit{\boldsymbol{H}}}_i^{\{i, m_j\}} $, $ m_j\in \mathcal{N}_i $, 表示$ \tilde{\mathit{\boldsymbol{M}}}_i $和$ \tilde{\mathit{\boldsymbol{H}}}_i $中与邻居节点$ S_{m_j} $对应的行.

本节将以绝对量测方程和相对量测方程(D1)和(D2)为例, 分析不可靠状态观测模型的具体表现形式. 为了更清晰讨论一般情况, 做以下分析:

丢包系数矩阵$ \mathit{\boldsymbol{L}}_{i, g, k} = {\rm diag} \{\pi_{i, g, k}^{\{i, m_1\}}, \cdots $, $ \pi_{i, g, k}^{\{i, m_j\}}, \cdots, \pi_{i, g, k}^{\{i, m_{|\mathcal{N}_i|}\}}\}^{\rm T} $ (参见表 1). 假定当$ g = 1 $时, 为第1种丢包情况发生, 此时$ \pi_{i, g, k}^{\{i, m_1\}} = \pi_{i, g, k}^{\{i, m_2\}} $ $ = \cdots = \pi_{i, g, k}^{\{i, m_{|\mathcal{N}_i|}\}} = 0 $, 表示节点$ S_i $与所有的邻居节点信息交换失败, $ \mathit{\boldsymbol{\gamma}}_{i, g, k} = 0 $, $ \mathit{\boldsymbol{z}}_{i, g}(k) = 0 $; 当$ g = 2^{|\mathcal{N}_i|} $时, 表示第$ 2^{|\mathcal{N}_i|} $种丢包情况发生, 此时$ \pi_{i, g, k}^{\{i, m_1\}} $ $ = \pi_{i, g, k}^{\{i, m_2\}} $ $ = $ $ \cdots $ $ = $ $ \pi_{i, g, k}^{\{i, m_{|\mathcal{N}_i|}\}} = 1 $, 表示节点$ S_i $与所有的邻居节点信息交换均成功, 丢包情况未发生, 此时$ \mathit{\boldsymbol{\gamma}}_{i, g, k} = \tilde{\mathit{\boldsymbol{M}}}_{i, k}\mathit{\boldsymbol{x}}_{\mathcal{N}_{[i]}}(k)+\mathit{\boldsymbol{\zeta}}_{i}(k) $, $ \mathit{\boldsymbol{z}}_{i, g}(k) = $ $ \tilde{\mathit{\boldsymbol{H}}}_{i, k}\mathit{\boldsymbol{x}}_{\mathcal{N}_{[i]}}(k) $ $ + $ $ \mathit{\boldsymbol{\upsilon}}_{i}(k) $. 因此, 绝对量测方程和相对量测方程展开表达为

根据丢包系数矩阵$ \mathit{\boldsymbol{L}}_{i, g, k} $的定义, 将式(D3)和式(D4)重写为

其中, $ \tilde{\mathit{\boldsymbol{M}}}_{i, g, k} = \mathit{\boldsymbol{L}}_{i, g, k} [\tilde{\mathit{\boldsymbol{M}}}_{i}^ {\{i, m_1\}}, \cdots, \tilde{\mathit{\boldsymbol{M}}}_{i}^{\{i, m_j\}}, \cdots $, $ \tilde{\mathit{\boldsymbol{M}}}_{i}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $和$ \tilde{\mathit{\boldsymbol{H}}}_{i, g, k} = \mathit{\boldsymbol{L}}_{i, g, k} [\tilde{\mathit{\boldsymbol{H}}}_{i}^{\{i, m_1\}}, \cdots $, $ \tilde{\mathit{\boldsymbol{H}}}_{i}^{\{i, m_j\}}, \cdots, \tilde{\mathit{\boldsymbol{H}}}_{i}^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $表示绝对量测方程和相对量测方程的量测矩阵$ \tilde{\mathit{\boldsymbol{M}}}_i $和$ \tilde{\mathit{\boldsymbol{H}}}_i $中含丢包信息交换下的邻居节点的量测矩阵, $ \mathit{\boldsymbol{\zeta}}_{i, g}(k) = \mathit{\boldsymbol{L}}_{i, g, k}[\mathit{\boldsymbol{\zeta}}_i(k) $, $ \mathit{\boldsymbol{\zeta}}_i(k), \cdots]^{\rm T} $和$ \mathit{\boldsymbol{\upsilon}}_{i, g}(k) = \mathit{\boldsymbol{L}}_{i, g, k}[\mathit{\boldsymbol{\upsilon}}_i(k), \mathit{\boldsymbol{\upsilon}}_i(k), \cdots]^{\rm T} $表示相应的量测噪声.

D.2   多量测的MMSE等价变换证明(拆分, 可靠量测)

有$ |\mathcal{N}_i| $个邻居节点的$ S_i $, 在不丢包情况下, 系统的时钟状态方程为

其中, 状态转移矩阵$ \mathit{\boldsymbol{A}} = {\rm blkdiag}\{\mathit{\boldsymbol{A}}_i, \mathit{\boldsymbol{A}}_{m_1}, \cdots, $ $ \mathit{\boldsymbol{A}}_{m_j}, \cdots, \mathit{\boldsymbol{A}}_{m_{|\mathcal{N}_i|}}\} $, 常数项$ \mathit{\boldsymbol{b}} = {\rm diag}\{\mathit{\boldsymbol{b}}_i, \mathit{\boldsymbol{b}}_{m_1}, \cdots, $ $ \mathit{\boldsymbol{b}}_{m_j} $, $ \cdots, $ $ \mathit{\boldsymbol{b}}_{m_{|\mathcal{N}_i|}}\} $, 过程噪声$ \mathit{\boldsymbol{w}}(k) $ $ = $ $ [\mathit{\boldsymbol{w}}_i(k) $, $ \mathit{\boldsymbol{w}}_{m_1}(k), \cdots, \mathit{\boldsymbol{w}}_{m_j}(k), \cdots, \mathit{\boldsymbol{w}}_{m_{|\mathcal{N}_i|}}(k)]^{\rm T} $是均值为0和协方差为$ \mathit{\boldsymbol{Q}} = {\rm diag}\{\mathit{\boldsymbol{Q}}_i, \mathit{\boldsymbol{Q}}_{m_1}, \cdots, \mathit{\boldsymbol{Q}}_{m_j}, \cdots $, $ \mathit{\boldsymbol{Q}}_{m_{|\mathcal{N}_i|}}\} $的高斯噪声.

1) 物理层面(PHY-level). $ S_i $相对其任意邻居节点$ S_{m_j} $的绝对量测方程和相对量测方程可分别表示为

2) 算法层面(ALG-level). 任意邻居节点$ S_{m_j} $相对$ S_i $的绝对量测方程为

算法层面将方程(D10)与方程(D9)交替组合, 得到以$ S_i $为中心的组合量测方程, 即对应于以下形式

将方程(D11)进行能观性变换(见附录C), 可以得到新的变换后的量测方程

组合所有邻居相对$ S_i $的变换, 则可分别得到式(D11)和(D12)的组集表达式为

其中, $ \bar{\mathit{\boldsymbol{Z}}}_i(k) $和$ \bar{\mathit{\boldsymbol{Y}}}_i(k) $分别是在算法层面上进行能观性行变换前与变换后节点$ S_i{ } $相对于所有的邻居节点的组合量测方程的状态量测值, 满足$ \bar{\mathit{\boldsymbol{Z}}}_i(k) = $ $ [\bar{\mathit{\boldsymbol{Z}}}_i^{\{i, m_1\}}, \cdots, \bar{\mathit{\boldsymbol{Z}}}_i^{\{i, m_j\}}, \cdots, \bar{\mathit{\boldsymbol{Z}}}_i^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $, $ \bar{\mathit{\boldsymbol{Y}}}_i(k) = $ $ [\bar{\mathit{\boldsymbol{Y}}}_i^{\{i, m_1\}} $, $ \cdots, \bar{\mathit{\boldsymbol{Y}}}_i^{\{i, m_j\}}, \cdots, \bar{\mathit{\boldsymbol{Y}}}_i^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $; $ \bar{\mathit{\boldsymbol{G}}}_i $和$ \bar{\mathit{\boldsymbol{\Gamma}}}_{i} $是相应的量测矩阵, 满足$ \bar{\mathit{\boldsymbol{G}}}_i = $ $ [\bar{\mathit{\boldsymbol{G}}}_i^{\{i, m_1\}}, $ $ \cdots, \bar{\mathit{\boldsymbol{G}}}_i^{\{i, m_j\}}, \cdots, \bar{\mathit{\boldsymbol{G}}}_i^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $, $ \bar{\mathit{\boldsymbol{\Gamma}}}_i = $ $ [\bar{\mathit{\boldsymbol{\Gamma}}}_i^{\{i, m_1\}} $, $ \cdots, $ $ \bar{\mathit{\boldsymbol{\Gamma}}}_i^{\{i, m_j\}}, \cdots, \bar{\mathit{\boldsymbol{\Gamma}}}_i^{\{i, m_{|\mathcal{N}_i|}\}}]^{\rm T} $; $ \bar{\mathit{\boldsymbol{\varpi}}}_i(k) $和$ \bar{\mathit{\boldsymbol{\psi}}}_i(k) $是均值为0, 方差分别为$ \mathit{\boldsymbol{R}}_{i, \bar{\varpi}_i} $和$ \mathit{\boldsymbol{R}}_{i, \bar{\psi}_i} $的高斯噪声, 满足$ \bar{\mathit{\boldsymbol{\varpi}}}_i(k) = [\bar{\mathit{\boldsymbol{\varpi}}}_i^{\{i, m_1\}}(k) $, $ \cdots $, $ \bar{\mathit{\boldsymbol{\varpi}}}_i^{\{i, m_j\}}(k), \cdots $, $ \bar{\mathit{\boldsymbol{\varpi}}}_i^{\{i, m_{|\mathcal{N}_i|}\}}(k)]^{\rm T} $, $ \bar{\mathit{\boldsymbol{\psi}}}_i(k) = $ $ [\bar{\mathit{\boldsymbol{\psi}}}_i^{\{i, m_1\}}(k) $, $ \cdots, $ $ \bar{\mathit{\boldsymbol{\psi}}}_i^{\{i, m_j\}}(k) $, $ \cdots, $ $ \bar{\mathit{\boldsymbol{\psi}}}_i^{\{i, m_{|\mathcal{N}_i|}\}}(k)]^{\rm T} $和$ \mathit{\boldsymbol{R}}_{i, \bar{\varpi}} = {\rm blkdiag}\{\mathit{\boldsymbol{R}}_{i, \bar{\varpi}_1}, \cdots, \mathit{\boldsymbol{R}}_{i, \bar{\varpi}_j}, \cdots, \mathit{\boldsymbol{R}}_{i, \bar{\varpi}_{|\mathcal{N}_i|}}\} $, $ \mathit{\boldsymbol{R}}_{i, \bar{\psi}} $ $ = {\rm blkdiag}\{\mathit{\boldsymbol{R}}_{i, \bar{\psi}_1}, \cdots, \mathit{\boldsymbol{R}}_{i, \bar{\psi}_j}, \cdots, \mathit{\boldsymbol{R}}_{i, \bar{\psi}_{|\mathcal{N}_i|}}\} $. 要想证明关于节点$ S_i $的方程(D13)与(D14)是否存在MMSE等价变换, 就要证明估计误差协方差矩阵的递推方程:

当设定相同的初始误差协方差矩阵时, 即$ \mathit{\boldsymbol{S}}_{i, 0}^{\rm pre} = \mathit{\boldsymbol{S}}_{i, 0}^{\rm post} $时, $ \mathit{\boldsymbol{S}}_{i, k}^{\rm pre} = \mathit{\boldsymbol{S}}_{i, k}^{\rm post} $, 那么方程(D13)与(D14)存在MMSE等价变换. 要使$ \mathit{\boldsymbol{S}}_{i, k}^{\rm pre} = \mathit{\boldsymbol{S}}_{i, k}^{\rm post} $, 根据附录C所述, 需要证明是否存在可逆矩阵$ \mathit{\boldsymbol{P}}_i $, 使得$ \bar{\mathit{\boldsymbol{\Gamma}}}_i = \mathit{\boldsymbol{P}}_i\bar{\mathit{\boldsymbol{G}}_i} $, $ \mathit{\boldsymbol{R}}_{i, \bar{\psi}} = \mathit{\boldsymbol{P}}_i\pmb{R}_{i, \bar{\varpi}} \mathit{\boldsymbol{P}}_{i}^{\rm T} $和$ \mathit{\boldsymbol{P}}_i\pmb{P}_{i}^{-1} = \mathit{\boldsymbol{P}}_i^{\rm T}(\mathit{\boldsymbol{P}}_i^{\rm T})^{-1} $.

3) 算法变换矩阵. 计算证明, 对于方程(D13), 存在这样一个可逆矩阵$ \mathit{\boldsymbol{P}}_i^I $满足方程(D13)到(D14)的变换, $ \mathit{\boldsymbol{P}}_i^I\in {\bf R}^{2|\mathcal{N}_i|\times 2|\mathcal{N}_i|} $满足: $ \mathit{\boldsymbol{P}}_i^I = {\rm blkdiag}\{\mathit{\boldsymbol{\rho}}, \mathit{\boldsymbol{\rho}}, \cdots, \mathit{\boldsymbol{\rho}}\} $, 其中$ \mathit{\boldsymbol{\rho}} = \begin{bmatrix} 1&0\\ 1&-1 \end{bmatrix} $. 通过$ \mathit{\boldsymbol{P}}_i^I $, 完成如式(D12)能观性线性变换. 而且, 可靠网络下节点多链路量测算法层面存在MMSE等价变换.

4) 物理变换矩阵. 参考量测方程(D8), 令量测模型(13)和(14)的物理变换矩阵$ \mathit{\boldsymbol{P}}_i = \mathit{\boldsymbol{P}}_i^{II} $. 同理, 则存在物理层面的MMSE等价变换(证明略).

D.3   多量测的MMSE等价变换证明(拆分, 不可靠量测)

由附录D.1可知, 在不可靠网络下, 节点$ S_i $的绝对量测方程和相对量测方程可分别表示为

考虑丢包以及$ \pi_{i, k}^{\{i, j\}} = \pi_{j, k}^{\{j, i\}} $, 参考式(D10), 根据式(D17)将节点$ S_i $的邻居节点关于节点$ S_i $的绝对量测方程改写为

本节将在算法层面对不可靠网络下单节点多链路的MMSE等价变换进行证明.

将式(D17)与式(D18)交替组合, 得到不可靠网络下节点$ S_i $的组合量测模型:

其中, $ \mathit{\boldsymbol{W}}_{i, g, k} = {\rm diag}\{\pi_{i, g, k}^{\{i, m_1\}}, \pi_{i, g, k}^{\{i, m_1\}}, \cdots, \pi_{i, g, k}^{\{i, m_j\}}, $ $ \pi_{i, g, k}^{\{i, m_j\}} $, $ \cdots, \pi_{i, g, k}^{\{i, m_{|\mathcal{N}_i|}\}}, \pi_{i, g, k}^{\{i, m_{|\mathcal{N}_i|}\}}\}\in {\bf R}^{2|\mathcal{N}_i|\times 2|\mathcal{N}_i|} $.

将式(D19)与式(D18)交替组合, 进行特定行变换, 得到不可靠网络下变换后的节点$ S_i $的组合量测模型:

通过对比上一节可靠网络下的MMSE等价变换证明, 我们可以知道要想证明在丢包情况下的MMSE等价, 同样只需证明是否存在可逆矩阵$ \mathit{\boldsymbol{P}}_{i, g} $, 使得$ \bar{\mathit{\boldsymbol{\Gamma}}}_{i, g} = \mathit{\boldsymbol{P}}_{i, g}\bar{\mathit{\boldsymbol{G}}}_{i, g} $, $ \mathit{\boldsymbol{R}}_{i, g, \bar{\psi}} = \mathit{\boldsymbol{P}}_{i, g}\mathit{\boldsymbol{R}}_{i, g, \bar{\varpi}}\mathit{\boldsymbol{P}}_{i, g}^{\rm T} $和$ \mathit{\boldsymbol{P}}_{i, g}\mathit{\boldsymbol{P}}_{i, g}^{-1} $ $ = \mathit{\boldsymbol{P}}_{i, g}^{\rm T}(\mathit{\boldsymbol{P}}_{i, g}^{\rm T})^{-1} $.

容易得到, 存在这样一个可逆矩阵$ \mathit{\boldsymbol{P}}_{i, g}^{I} $, 且$ \mathit{\boldsymbol{P}}_{i, g}^{I} $ $ = $ $ \mathit{\boldsymbol{P}}_i^{I} $. 因此不可靠网络下单节点多链路量测在算法层面存在MMSE等价变换.

同理, 在物理层面上, 可以证明存在可逆矩阵$ \mathit{\boldsymbol{P}}_{i, g}^{II} $, 且$ \mathit{\boldsymbol{P}}_{i, g}^{II} = \mathit{\boldsymbol{P}}_{i}^{II} $. 因此不可靠网络下单节点多链路量测在物理层面同样存在MMSE等价变换.