基于MRF最大后验估计和Gardner环的WSN时间同步算法matlab仿真-育师

1.本算法理论知识

1.1 马尔科夫随机场的最大后验估计

1.2 Gardner环

3.部分源码

4.仿真结论

无线传感器网络（WSN）由大量低成本、低功耗传感节点组成，节点间协同采集环境数据，精准时间同步是数据融合、时隙调度、目标定位的基础。传统时间同步算法如TPSN、FTSP大多仅校正节点时钟偏移，忽略信道噪声、晶振频率漂移带来的相位抖动与频率偏差；同时多数算法采用开环校正机制，无法持续跟踪时变时钟参数。Gardner位同步环是经典的闭环定时同步结构，能够实时跟踪信号相位偏差，适合动态修正节点本地时钟相位，但独立Gardner环缺乏全网协同约束，仅能实现点对点同步。马尔科夫随机场（MRF）具备空间邻域关联建模能力，结合最大后验概率估计（MAP）可以融合全网相邻节点时钟观测信息，构建全局时钟约束模型，抑制测量噪声带来的同步误差。

1.本算法理论知识

通过一种基于马尔科夫随机场的最大后验估计方法对无线传感器网络中不相邻的两个接收节点在多个参考广播消息条件下的相位偏差进行估计，然后通过Gardner环对时钟漂移进行实时的估计和跟踪，实现时间同步过程。

1.1 马尔科夫随机场的最大后验估计

MRF核心性质为马尔科夫局部性：任意节点随机变量的概率分布仅依赖其邻域节点，与网络其余节点条件独立。将WSN每个传感器节点映射为MRF图模型顶点，存在通信链路的相邻节点构建边。定义随机场：

其中，N为网络节点总数。

根据贝叶斯理论，最大后验估计求解最优时钟参数组合，最大化后验概率：

1.2Gardner环

Gardner环是无数据辅助定时同步算法，不需要预先获知发送码元信息，非常适合低带宽WSN同步报文。其结构如下图所示：

Gardner环主要由四部分组成：内插器、时钟误差提取模块、环路滤波器以及控制器模块。

Gardner环利用每个码元中间采样点与码元边界采样点计算定时误差函数。设接收离散采样序列r(kT_s)，T_s为采样间隔，码元周期T。定时误差检测函数：

误差信号e(n)反映本地采样时钟与发送端时钟相位偏差。误差经过环路滤波后输出控制量调整本地采样时钟频率与相位。环路滤波器采用一阶低通结构：

k_p比例系数，k_i积分系数，输出u(n)调整节点本地时钟计数器，持续修正相位漂移。在本文算法中，Gardner环接收MRF-MAP输出的时钟粗参数作为环路初始值，完成高精度闭环动态跟踪。

3.部分源码

clc; clear; close all; warning off; rng('default') rng(3); %网络区域范围 %25000mm x 25000mm RK = 25000; CK = 25000; %网络节点个数 %100 Nnote = 100; %节点通信半径 %40m Rad = 4000; %延时初始值 %0 delay0 = 0; %误差初始值 %0 error0 = 0; %判决次数 %250 Ncheck = 250; %两次判决时延 %6.5s Dtime = 6.5; %仿真时间 %35min Stime = 35*60; %数据包传输速率 %50kb/s Rate = 50e3; Ve = 3e8; %能耗系数 alpha = 1e-3; %测量值个数 N = 50; W = 0.1; %点随机分布 X = RK*rand(1,Nnote); Y = CK*rand(1,Nnote); figure; plot(X,Y,'b.'); for i = 1:Nnote %RBS算法要求每个节点与其他N个节点交换该消息的本地时间，计算两者之间的本地时间差 Rnot = i; indx = 0; dist = []; De = []; Vc = []; %时钟漂移估计 for j = 1:Nnote if (i == j) == 0 dist = sqrt((X(i) - X(j))^2 + (Y(i) - Y(j))^2); %time indx = indx + 1; De(indx) = 1e6*dist/Ve; end end %原来的算法 Dey = De + randn(size(De)); delta = (sum(Dey.*De)-length(De)*mean(De)*mean(Dey))/(sum(De.*De)-length(De)*mean(De)*mean(De)); v = mean(Dey)-delta*mean(De); TDe = []; for ii = 1:length(De) TDe(ii) = v*De(ii) + delta; end %增加Gardner环 TDe = func_gardner(TDe); %计算不相邻之间的相位偏差 indxx = 0 ; fai = 0; for k1 = 1:Nnote-1 for k2 = k1:Nnote-1 %不相邻编号节点 if ((k1 == k2) == 0) & (k1 > k2+1 | k1 < k2+1) %随机大小的参考消息数量 m = round(20*rand)+1; for k3 = 1:m deltaT(k3) = abs(TDe(k1) - TDe(k2))*(1 + rand/10); end indxx = indxx + 1; fai(indxx) = sum(deltaT)/m; end end end %势 Vc = deltaT.^2/length(deltaT); %N个测量值 tmps = randperm(length(fai)); x = fai(1:N); U = sum(Vc); %fai估计值 deltan = 11.357; delta = 11.357; %近似计算 fai_(i)= abs((1-W)*sum(x)/deltan^2 * (deltan^2*delta^2)/(N*delta^2 + deltan^2) - W*U); end A12-26