一种改进粒子滤波算法实现的多径参数估计

国 强，刘雪萌，周 凯

(哈尔滨工程大学 信息与通信工程学院，黑龙江 哈尔滨 150001)

随着全球导航卫星系统(Global Navigation Satellite System，GNSS)的应用领域日益广泛，人们对其定位精度的要求越来越高。但由于卫星信号在传播过程中会经障碍物产生反射[1]，致使接收到的信号往往由直达信号与多径信号两部分组成，进而引起多径误差。多径误差具有位置上不相关及时间上不确定特性，使其很难通过差分技术消除，因此它是影响全球导航卫星系统定位精度的主要误差源之一[2]。

目前，对于全球导航卫星系统中抗多径干扰研究的方法大致可分为3类：一是基于接收机天线的设计[3-5]；二是基于接收机跟踪环路的改进[6-7]；三是基于数据处理的多径抑制[8-9]。前两者在实现过程中均面临一定困难：基于天线端的设计，会增加天线成本，并不是一个经济的选择；而基于接收机跟踪环路的改进则需要消耗较多的硬件资源。因此，基于数据处理的多径抑制技术得到了广泛应用。数据处理技术的核心是信号参数估计，也正是笔者所研究的内容。

文献[10]采用一种基于高斯和近似的扩展切片高斯混合滤波器实现直接信号时延和多径相关输出的参数估计，但其算法复杂度较高。随着压缩感知技术的兴起，一种基于压缩感知技术的GPS多径估计方法被提出[11]，而最新研究显示传统的压缩感知以稀疏结构为先验信息来进行信号恢复，使得压缩感知的性能并不高。文献[12]针对雷达系统参考信道中存在多径干扰的问题，提出一种基于改进恒模算法的参考通道多径干扰抑制算法，该算法与传统的盲均衡方法相比，具有较好的收敛性，多径抑制性能有所提高。在信号特征提取领域，针对多径信号产生的影响，一种基于时变AR模型的多径信号抑制算法被提出[13]，该算法最终可精确地估计出多径时延及幅度因子，进而实现多径信号的逐级抑制。由于大多数多径抑制算法仅限于高斯噪声，导致在非高斯环境下这些算法性能降低，为了解决这一问题有学者采用基于相关熵准则的自适应滤波多径抑制算法[14]，提高多径抑制性能。近年来，基于粒子滤波(Particle Filter，PF)算法的参数估计越来越普遍，王宪等[15]提出一种基于粒子滤波的GPS多径估计，虽然解决了算法运算量过大的问题，但是随着时间的增加粒子退化、粒子多样性降低现象仍然存在。因此基于PF的改进算法大量涌现：在状态和参数联合估计过程中采用双重粒子滤波算法[16]可避免粒子退化现象，但其计算量相对较大；文献[17]中提出基于遗传粒子滤波(Genetic Algorithm-Particle Filter，GA-PF)算法实现的多径跟踪，可以缓解粒子滤波算法中粒子多样性降低的现象，但由于遗传算法自身存在着过早收敛的问题，导致多径跟踪的效果并不理想；利用粒子群优化算法(Particle Swarm Optimization，PSO)使得采样粒子集向似然函数较大的区域移动[18]，该算法虽然减小了参数估计所需的粒子数目但却没有解决粒子退化问题；基于扩展卡尔曼滤波(Extended Kalman Filter，EKF)改进的粒子滤波[19]，利用EKF产生粒子滤波的重要性密度函数，然而，由于实际上很多系统的状态方程是非线性的，使得EKF在线性化程中忽略了泰勒展开式的高阶项[20]，产生截断误差；基于差分进化(Differential Evolution，DE)改进粒子滤波的多径估计算法[21]，在差分进化算法中缩放因子F与交叉因子CR的选取严重影响算法性能，使得差分进化算法的寻优能力和收敛速度相互矛盾，参数估计存在一定的误差。

为克服上述问题，笔者提出一种无迹差分粒子滤波算法(Unscented Improved Differential Evolution-Particle Filter，UIDE-PF)。该算法首先将无迹卡尔曼波(Unscented Kalman Filter，UKF)算法“嵌入”到粒子滤波算法重要性采样阶段，得到合理的采样分布，避免了粒子退化问题。然后，采用改进差分进化(Improved Differential Evolution，IDE)算法对重要性采样后的粒子进行自适应变异、自适应交叉以及选择操作更新粒子，从而提高粒子多样性，进而减小参数估计结果的波动幅度及均方根误差。

1 系统模型建立

全球导航卫星系统在进行定位时，其接收机接收到的信号通常可表示为

r0(t)=A0c(t-τ0)cos(ω0t+φ) ，

(1)

其中，A0表示直达信号幅度；φ表示直达信号相位延迟；c(t)表示扩频码序列；τ0表示直达信号时延；ω0表示载波相位；式中忽略导航电文。然而，实际接收信号还包含多径及噪声干扰，此时信号模型表示为

r1(kTs)=A0c(kTs-k0Ts)cos(ω0kTs+φ0)+

(2)

式(2)表示数字化后中频信号，M表示多径信号的路数；Ai、kiTs、φi分别表示第i路多径信号的幅值、相对时延、相对相位延迟，其中，k表示第k个采样点，Ts表示采样周期；k0=τ0/Ts，表示τ0对应采样点数；ki=τi/Ts，表示τi对应采样点数；n(kTs)表示高斯白噪声。对于静态环境而言，一般是指接收机周围环境变化相对较小，如平原、村落等地。此时接收机接收信号以双径模型为例[18]，具体形式表示如下：

r1(kTs)=A0c(kTs-k0Ts)cos(ω0kTs+φ0)+

A1c(kTs-k0Ts-k1Ts)cos(ω0kTs+φ0+φ1)+n(kTs) 。

(3)

(4)

Y=[R1(A0，Ai，η，ki，d1)，R2(A0，Ai，η，ki，d2)，…，RL(A0，Ai，η，ki，dL)] ，

(5)

其中，L表示相关器个数。

为了将无迹差分粒子滤波算法用于多径参数估计，需要利用p(xk|xk-1)描述参数的演化特性[15]：

p(xi|xi-1，xi-2，…，x1)=p(xi|xi-1) 。

(6)

根据式(6)得到系统具体模型，其状态方程与观测方程如下：

(7)

其中，H表示状态转移矩阵；ω表示过程噪声；R表示观测矩阵；v表示观测噪声；Xn表示待估计参数。笔者的主要目的是用观测向量Y估计出状态向量Xn。

2 基于改进算法的多径参数估计

2.1 改进差分进化算法

差分进化(Differential Evolution，DE)算法是一种智能优化算法，标准差分进化算法主要包括以下步骤[22]：

(1) 初始化种群。其表达式为

(8)

(2) 变异过程。在差分进化算法中变异策略表达式为

(9)

(3) 交叉过程。该过程的主要目的是产生一个试验向量，其交叉过程表达式为

(10)

其中，CR一般在0～1之间。CR选取决定着目标向量xj，i和变异向量vj，i在试验向量中所占比例的大小。

(11)

其中，f(·)表示适应度函数。由于笔者的主要目的在于实现信号参数估计，因此可将参数估计过程中直达信号幅度、多径信号幅度、本地估计偏差、多径相对时延4种参数估计结果的绝对误差之和作为适应度函数，并保留适应度函数较小的个体。

(5) 边界条件处理。对于超出范围的个体按如下表达式处理：

(12)

由于差分进化算法中的F和CR往往是固定值，不能满足算法各个阶段的需求。因此笔者提出一种自适应策略，使F和CR动态变化，改进后的F和CR表达式如下：

(13)

其中，G表示最大进化次数，g表示当前进化次数。由式(13)可以看出，在进化前期g取值较小，此时获得较大的F和较小的CR，有利于提高算法的全局搜索能力；在进化后期g取值较大，此时获得较小的F和较大的CR，有利于局部搜索并提高收敛速度。

2.2 多径参数估计的实现过程

尽管粒子滤波算法可以应用在任何非高斯非线性系统中，但其存在的粒子退化、多样性降低问题仍不可避免。因此，针对该问题采用无迹差分粒子滤波算法实现信号参数估计，基本原理如下：

(1) 初始化。由于多径信号的幅度小于直达信号的幅度，且归一化幅度均在0～1之间；多径信号延时一般在0～1.5码片之间，即各参数的取值范围如下：

(14)

其中，A0、A1分别表示直达信号、多径信号幅度；τ表示多径相对时延。按式(14)产生NP个粒子。

(2) 重要性采样。利用UKF算法为每个粒子计算其均值及协方差，并不断更新粒子。

① 利用无迹变换(Unscented Transformation，UT)变化计算2n+1个Sigma点集：

(15)

其中，P表示初始状态协方差矩阵；i表示对应的采样点；n表示状态参数的维数，文中n=4。

② 计算采样点的相应权值：

(16)

其中，m表示均值；c表示协方差；上标表示对应的采样点；λ表示缩放比例；α和β表示常数。

③ 对2n+1个Sigma点集进行下一步预测：

X(i)(k+1|k)=f[k，X(i)(k|k)]，i=1，…，2n+1 ，

(17)

其中，f(·)表示状态函数。

④ 对步骤③中的Sigma点集加权求均值进而得到状态向量的一步预测结果，并计算协方差矩阵：

(18)

⑤ 根据式(18)再次使用无迹变换，得到新的Sigma点集如下：

(19)

⑥ 将式(19)代入式(7)中的观测方程，得到一步预测观测值：

Z(i)(k+1|k)=h[X(i)(k+1|k)] ，

(20)

其中，h表示观测函数。

⑦ 将步骤⑥得到的观测值通过加权平均，求得系统预测均值及协方差：

(21)

⑧ 计算卡尔曼增益：

(22)

⑨ 更新状态向量及协方差：

(23)

(3)重采样。采用IDE算法代替传统重采样过程，得到最优粒子并计算其归一化权重：

(24)

其中，i表示样本中第i个粒子；k表示迭代次数。

(4)根据步骤(1)～(3)可以得到采用无迹差分粒子滤波算法实现多径参数估计的结果，其表达式如下：

(25)

根据以上对无迹差分粒子滤波算法原理的分析，给出其流程框图，如图1所示。

图1 无迹差分粒子滤波算法流程框图

图2描述了在第500次迭代时采用无迹差分粒子滤波算法(UIDE-PF)与标准粒子滤波算法(PF)在参数估计过程中粒子权重的大小。从图2(a)可以看出，采用标准粒子滤波算法在参数估计过程中，大权重的粒子不断被复制，小权重的粒子逐渐被淘汰，导致粒子多样性逐渐降低；而从图2(b)可以看出，采用无迹差分粒子滤波算法在该过程中几乎不存在粒子权重相等的现象。由此可见，无迹差分粒子滤波算法相比于标准粒子滤波算法可有效避免粒子多样性降低问题。

(a) 采用PF算法得到的粒子权重

3 仿真结果及分析

为了验证无迹差分粒子滤波算法在参数估计过程中的有效性，将采用该算法得到的参数估计结果与无迹卡尔曼粒子滤波(Unscented Kalman Particle Filter，UPF)算法、PF算法及文献[23]中的算法进行对比分析。为了实验的严谨性，4种算法参数保持一致且与文献[23]中算法参数设置相同：粒子数M=40；过程噪声协方差Q=10-4eye(4)；各支路与即时码之间的相关间距为±0.5 Tc，±0.3 Tc，±0.1 Tc；状态转移矩阵F=eye(4)；各参数初始状态X0=[0.5 0.7 0.2 0.4]T。其中，UKF算法中状态维数L=4；缩放比例参数λ=α2L-L，采样点分布状态由α来确定，经过多次实验得出当α=0.001时算法性能较好；而权系数β通常是非负数，其取值通常为β=2。IDE算法中Fmax=0.6，Fmin=0.2，CRmax=0.8，CRmin=0.2。仿真过程进行500次迭代，仿真结果如图3所示。

图3描述了各参数的估计结果及其对应的均方根误差，其中扩展卡尔曼粒子滤波(Extended Kalman Particle Filter，EPF)算法所代表的为文献[23]中所提及的算法。从图中可以看出文献[23]中所采用的算法虽然解决了EKF算法对初值敏感的问题，但其线性化过程引起的误差，导致估计结果波动幅度及均方根误差相对较大；而UIDE-PF算法则有效减小了参数估计结果的波动幅度及均方根误差。

(a) 直达信号幅度估计

图4描述了采用上述估计结果实现多径信号重构，并去除多径干扰后的相关函数值。从图中可以看出采用UIDE-PF算法进行多径抑制后得到的相关函数曲线更平滑，且最接近真实值。同时，从图的形状上来看，采用UIDE-PF算法得到的相关函数曲线图在码延时为0时，对应的幅值恰好为曲线的峰值点且附近没有发生畸变。

(a) 采用PF算法实现多径抑制

表1给出了4种不同算法每次迭代所消耗的时间。由于UPF算法与UIDE-PF算法在参数估计过程中需要对每个粒子进行均值及协方差的计算，导致两种算法时间消耗较长。但是对于以时速为30 m/s行驶的汽车而言，要使接收机定位精度达到2 m左右，即算法完成一次运行时间不应超过66.67 ms。而笔者所提出的算法运行一次大约需15 ms，满足接收机在高速环境下的定位需求。

表1 4种不同算法每次迭代所消耗的时间

4 结束语

笔者提出了一种无迹差分粒子滤波算法。该算法可有效地解决粒子滤波算法在参数估计过程中存在的粒子退化及多样性降低的问题。相对于文中提到的粒子滤波算法、无迹粒子滤波算法及对比文献中的算法，笔者所提算法在满足了接收机对定位实时性要求的情况下，提高了参数估计精度，减小了估计结果均方根误差。接下来，将针对如何减小算法的运算量，进一步提高参数估计精度展开研究。