[发明专利]一种基于DPM模型的RB粒子滤波时间同步方法

权利要求书

1.一种基于DPM模型的RB粒子滤波时间同步方法，其特征在于该方法包括以下步骤：

(1)将双向时间信息交换模型等效成动态马尔可夫模型；把双向时间信息交换过程中的不确定时延看作是马尔可夫过程的观测噪声；把时钟偏差看成是马尔可夫模型的状态变量；

(2)马尔可夫动态模型中观测噪声未知；为了提高噪声估计准确度，利用DPM模型描述观测噪声分布，确定待估计的未知参数Θ=(θ₁,…,θ_L)；

(3)利用RB粒子滤波估计新状态变量ψ＝[Θ,x]，利用粒子滤波估计Θ，利用卡尔曼滤波估计x，得到x的估计最终得到时钟偏差估计

(4)调整节点A的时钟，消除时钟偏差实现A、B节点时间同步。

2.如权利要求1所述的基于DPM模型的RB粒子滤波时间同步方法，其特征在于步骤（2）中的DPM模型描述观测噪声分布为

n_k～f(n_k|θ_k)

θ_k～G      (1)

G～DP(G₀,α)

其中，θ_k是未知参数，G₀是基础分布，α是集中度参数；每次进行双向时间信息交换时，观测噪声n_k可能服从不同分布，L次信息交换后，观测噪声服从混合分布，未知参数Θ=(θ₁,…,θ_L)；

利用DPM模型描述观测噪声分布，需要估计未知参数Θ=(θ₁,…,θ_L)；要实现时间同步，需要估计状态变量x；定义新状态变量ψ＝[Θ,x]表示所有需要估计的变量。

3.如权利要求1所述的基于DPM模型的RB粒子滤波时间同步方法，其特征在于步骤(3)中估计新状态变量ψ的具体方法：利用RB粒子滤波估计新状态变量ψ；将Θ看成是ψ的非线性部分，将x看成是ψ的线性部分；利用粒子滤波估计Θ，利用卡尔曼滤波估计x，得到x的估计

利用RB粒子滤波估计ψ，具体过程是：

(a)依据贝叶斯理论，状态变量ψ的后验概率密度函数可分解为：

(b)利用粒子滤波估计Θ；

(b.1)根据θ_k的全条件分布，采样N个粒子：

θ₁|α,G₀～G₀

θ2|θ1,α,G0~1αδθ1+αα+1G0]]>

       （3）

…

θN|θ1,...,θN-1,α,G0~Σi=1N-11α+i-1δθi+αα+n-1G0]]>

得到粒子集{θ₁,…,θ_N}，初始化粒子权重w_i,i＝1,…,N为1/N；可表示为：

(b.2)利用观测数据重新计算粒子权重，并归一化粒子权重；计算粒子θ_i权重的过程为：依次计算粒子θ_i=[μ_i,Σ_i]表示的高斯分布N(μ_i,Σ_i)产生噪声n_k＝z_k-2x_k,k＝1,…,L的概率：

p(nk|θi)=12πΣiexp(-(nk-μi)22Σi)k=1,...,L---(5)]]>

选取p(n_k|θ_i),k＝1,…,L中最大值maxp(n_k|θ_i)，计算粒子权重将使p(n_k|θ_i),k＝1,…,L取得最大值的噪声n_k看作由高斯分布N(μ_i,Σ_i)产生，可理解为噪声n_k为高斯分布N(μ_i,Σ_i)所拥有；按照上述步骤，依次计算完粒子权重后，一个噪声可能被多个高斯分布共有，或者，有的噪声还没有确定分布；因为一个噪声值只能来源于一个高斯分布，对于第一种情况，比较多个高斯分布产生该同一噪声的概率，确定概率最大的高斯分布作为该噪声的分布，并将其它高斯分布所对应的粒子从粒子集中移除；对于没有确定分布的噪声值n_i，将高斯分布N(μ_i,Σ_i)作为其分布，Σ_i由G₀产生，并将粒子θ_i＝[n_i,Σ_i]加到粒子集中，计算粒子权重经过以上处理后，每个噪声都确定了自己服从的高斯分布，得到了新粒子集{θ₁,…,θ_m}，1≤m≤L，相应的粒子权重为{w₁,…,w_m}；依据粒子权重，从粒子集{θ₁,…,θ_m}，1≤m≤L中重新采样N个粒子{θ′₁,…,θ′_N}，重采样粒子的权重为{w′₁,…,w′_N}，计算重采样粒子的归一化权重：

wi=wi′Σi=1Nwi′---(6)]]>

后验可重新表示为：

(b.3)Θ的估计值为：

Θ＝(θ₁,…,θ_m)   （8）

利用粒子滤波估计出参数Θ后，现在利用卡尔曼滤波估计状态变量x；因为观测噪声n_k服从高斯分布，且动态马尔可夫模型是线性高斯的，可表示为：

其中，为高斯概率密度函数N(x_i(θ′_i),P(θ′_i))，x_i(θ′_i)，Pi(θ′_i)由卡尔曼滤波得到；因为不同粒子θ′_i,i＝1,…N代表着不同高斯分布的均值和方差，不同粒子θ′_i对应不同的x_i(θ′_i)，P_i(θ′_i)，需要对每个粒子进行一次卡尔曼滤波：

x^i(θi′)=x0]]>

P^i(θi′)=P0+Σi′]]>

Bi(θi′)=2P^i(θi′)[4P^i(θi′)+Q]-1---(10)]]>

xi(θi′)=x^i(θi′)+Bi(θi′)[zi-2x^i(θi′)]]]>

Pi(θi′)=P^i(θi′)-2Bi(θi′)P^i(θi′)]]>

x₀是状态变量初始值，P₀是初始协方差；

x^=Σi=1Nwixi(θi′)]]>

    （11）

P^=Σi=1Nwi{Pi(θi′)+[xi(θi′)-x^][xi(θi′)-x^]T}]]>

是状态变量x的估计。