一种基于拒绝采样的粒子滤波算法的非线性动态系统信号 处理方法 技术领域 本发明涉及基于粒子滤波算法的信号处理方法， 要求保护的技术方案属于信号处 理、 人工智能和计算机视觉领域。
背景技术 动态系统的状态估计问题涉及很多领域， 尤其是信号处理、 人工智能和计算机视 觉领域。传统的卡尔曼滤波只适用于线性高斯系统， 而扩展卡尔曼滤波也只能应对系统的 弱非线性。因此， 适用于非线性、 非高斯系统的粒子滤波备受关注。
粒子滤波是一种基于蒙特卡洛模拟和递推贝叶斯估计的滤波方法。它采用粒子 描述状态空间， 使用一组带权重的粒子近似表示系统的后验概率密度， 并通过模型方程和 观测信息实现递推的估计过程。常见的粒子滤波算法包括 SIR 粒子滤波、 辅助粒子滤波 (APF)、 正则粒子滤波 (RPF)、 高斯粒子滤波 (GPF) 和无迹粒子滤波 (UPF)。
粒子滤波算法中的两个关键技术是建议分布的选取和重采样算法。 虽然好的建议 分布可以延缓粒子权重的退化问题， 但最优的建议分布在一般情况下是无法获取的。 同时， 虽然重采样算法解决了粒子权重的退化问题， 但常见的重采样算法只是对粒子的简单复制 和剔除， 从而又导致了粒子匮乏问题 ( 粒子多样性的丧失 )。
发明内容 为了克服现有基于粒子滤波算法的动态系统信号处理方法中的粒子权重退化和 粒子匮乏、 无法有效处理非线性和非高斯信号的不足， 本发明提出一种能有效减小粒子权 重退化、 粒子多样性好， 有效处理非线性和非高斯的信号的基于拒绝采样的粒子滤波算法 的非线性动态系统信号处理方法。
为了解决上述技术问题提出的技术方案为 ：
一种基于拒绝采样的粒子滤波算法的非线性动态系统信号处理方法， 采用粒子描 述动态系统的状态空间， 设非线性动态系统的状态空间模型为 ：
xk ＝ f(xk-1)+vk-1
zk ＝ h(xk)+nk
其中， xk 和 zk 分别表示系统在 k 时刻的状态和观测值， f(xk-1) 和 h(xk) 分别表示 系统的状态转移方程和观测方程， vk-1 和 nk 分别表示系统噪声和观测噪声 ；
对于每一个新采样的粒子， 首先计算其被接受的概率， 然后判断其是否被接受， 递 推过程包括以下步骤 ：
第一步， 根据 k-1 时刻的 N 个粒子 xk-1i， i ＝ 1， 2，…， N， N 为自然数， 估计观测状 态 zk 的发生概率 p(zk) ：
其中， μki 为给定 xk-1i 时 xk 的某个特征值， 取为均值 E(xk|xk-1i) 或一个采样值 xki ： p(xk|xk-1i)。
第二步， 初始化 k 时刻接受的粒子数 R 为 0 ；
第三步， 若 R ＜ N， 则重复下述过程 ： 对于每一个粒子 xk-1i， 经采样得到 xkij， 计算 i ij i i ij ij i xk 的接受概率为 p · r ， 其中系数 p 与粒子 xk-1 的采样次数 m 成反比， 且 r 的计算方法 如下 ：
第四步， 从接受的 R 个粒子中选取 N 个代表系统状态的后验概率密度。 第五步， 输出系统状态的估计值本发明的技术构思为 ： 该算法的一个递推过程包括以下基本步骤 ：
1)、 根据 k-1 时刻的 N 个粒子估计观测状态 zk 的发生概率。
2)、 初始化 k 时刻接受的粒子数 R 为 0。
3)、 若 R ＜ N， 则重复下述过程 ： 对于每一个粒子 xk-1i， 经采样得到 xkij， 计算 xkij 的 接受概率为 pi·rij， 其中系数 pi 与粒子 xk-1i 的采样次数 mi 成反比。
4)、 从接受的 R 个粒子中选取 N 个代表系统状态的后验概率密度。
5)、 输出。
本发明具有以下优点 ：
1、 采用拒绝采样算法使接受的粒子具有相同的权重， 省去了重采样算法， 避免了 传统粒子滤波算法中的粒子权重退化和粒子匮乏问题。
2、 通过选择不同的采样接受策略可以减少算法的运行时间。
附图说明 图 1 为基于拒绝采样的粒子滤波算法流程图。
图 2 为基于拒绝采样的粒子滤波算法与其他粒子滤波算法在不同粒子数下的均 方根误差均值对比曲线示意图。
图 3 为基于拒绝采样的粒子滤波算法与其他粒子滤波算法在不同粒子数下的运 行时间对比曲线示意图。
图 4 为基于拒绝采样的粒子滤波算法在不同 mi-pi 下的均方根误差均值对比曲线
示意图。
图 5 为基于拒绝采样的粒子滤波算法在不同 mi-pi 下的运行时间对比曲线示意图。具体实施方式
下面结合附图和实施例对本发明作进一步说明。
参照图 1 ～图 5， 一种基于拒绝采样的粒子滤波算法的非线性动态系统信号处理 方法， 采用粒子描述动态系统的状态空间， 设非线性动态系统的状态空间模型为 ：
xk ＝ f(xk-1)+vk-1
zk ＝ h(xk)+nk
其中， xk 和 zk 分别表示系统在 k 时刻的状态和观测值， f(xk-1) 和 h(xk) 分别表示 系统的状态转移方程和观测方程， vk-1 和 nk 分别表示系统噪声和观测噪声 ；
对于每一个新采样的粒子， 首先计算其被接受的概率， 然后判断其是否被接受， 递 推过程包括以下步骤 ：
第一步， 根据 k-1 时刻的 N 个粒子 xk-1i(i ＝ 1， 2，…， N) 估计观测状态 zk 的发生 概率 p(zk) ：
其中， μki 为给定 xk-1i 时 xk 的某个特征值， 取为均值 E(xk|xk-1i) 或一个采样值 xki ： p(xk|xk-1i)。
第二步， 初始化 k 时刻接受的粒子数 R 为 0 ；
第三步， 若 R ＜ N， 则重复下述过程 ： 对于每一个粒子 xk-1i， 经采样得到 xkij， 计算 i ij i i ij ij i xk 的接受概率为 p · r ， 其中系数 p 与粒子 xk-1 的采样次数 m 成反比， 且 r 的计算方法 如下 ：
第四步， 从接受的 R 个粒子中选取 N 个代表系统状态的后验概率密度。 第五步， 输出系统状态的估计值本实施例通过一个非线性动态系统的状态估计对本发明和其它几种粒子滤波算 法进行比较。系统的状态空间模型如下 ：
其中， 系统噪声方差和观测噪声方差分别取为 10 和 1。p(zk) 的计算采用如下方法：
粒子的采样接受使用如下分治策略 ： 对于 p(zk|μki) ＜ p(zk) 的粒子， mi ＝ m1 ＝ 1， pi ＝ p1 ＝ 5 ； 对于其他粒子， mi ＝ m2 ＝ 5， pi ＝ p2 ＝ 1。设定观测时间为 100， 运行次数 为 100， 在粒子数 N 分别取 50、 100、 150、 200、 250、 300 时， 本发明提出的算法与其他粒子滤波 算法所产生的均方根误差 (RMSE) 均值和所运行的时间分别如图 2 和图 3 所示。
图 2 中， 本发明算法的 RMSE 均值在粒子数目小于 150 时明显优于其它算法， 在粒 子数目大于 150 时和 UPF 算法接近。图 3 中， 本发明算法在不同粒子数目下的运行时间略 多于 APF 算法， 约为 UPF 算法的 55％。
图 4、 图 5 分别为本发明算法在上述分治策略中 m1、 p2 取 1， p 1、 m2 取不同值时的 RMSE 均值和运行时间对比曲线。在不同的 p1、 m2 取值下， RMSE 均值基本不变， 而运行时间 可以部分减少。