一种抑制多径干扰的码跟踪环路及方法与流程

本发明涉及卫星导航技术领域，具体涉及一种抑制多径干扰的码跟踪环路及方法。

背景技术：

多径干扰在森林、城市峡谷、室内等弱信号环境中广泛存在。传统卫星导航接收机由于经过穿透损耗，信号的幅度和相位有所改变，使得信号严重失真，进而导致接收机对卫星导航信号的码跟踪将难以实现。鉴于此，针对微弱卫星信号需研究一种具有高精度且能够抑制多径干扰的码跟踪方法。

国内外相关学者对如何减小一码片内多径干扰带来的伪距测量和定位误差问题做了大量的研究。目前，利用信号估计理论实现多径抑制的方法主要包括：最大似然估计，medll(multipletipatheliminationdll)算法，卡尔曼滤波(kalmanfilter，kf)算法，扩展卡尔曼滤波(extendedkalmanfilter，ekf)算法，无迹卡尔曼滤波算(unscentedkalmanfilter，ukf)法，粒子滤波(particlefilter，pf)算法等。

最大似然估计算法是在解扩前采用最大似然估计原理对中频数据进行多径估计，恢复出接近真实值的直达信号，并用于后续的解扩处理。medll算法通过将信号估计与环路结合到一起去减小多径误差，该技术须通过多个相关器得到的大量量测后再做最大似然估计，从而得到直射信号与多径信号的位置和幅度估计。kf算法处理的是线性问题，对于非线性滤波问题处理效果较差；解决非线性问题最经典的方法是ekf算法，但在用ekf算法进行码跟踪时，只能对弱非线性系统得到理想的滤波结果，而对于强非线性系统，ekf算法在线性化过程中会引入更多的截断误差，使得码跟踪精度降低；ukf算法是基于无迹变换和ekf的算法框架，对于任意的非线性系统，使用ukf都能够获得精确到三阶矩的系统后验均值和协方差估计，但与ekf一样对非线性系统的后验概率密度需做高斯假设，对于一般的非高斯分布模型该方法则不适用。解决较强的非线性问题常用方法是粒子滤波，目前已有文献提出把粒子滤波应用于卫星导航接收机码跟踪中，以抑制码跟踪环中的多径干扰，但是在这类滤波估计中，观测值是唯一可获得的信息源，然而由于噪声误差、系统误差、观测误差和设备误差等误差的存在，观测值必将带有一定的误差，这些误差的存在会严重影响滤波估计的精度，降低环路跟踪的性能；同时，粒子滤波算法中将仅把转移密度作为重要性密度函数，采样粒子时并没有把当前时刻量测信息考虑进去，使得可获得量测信息未得到充分的利用，大大影响滤波的估计精度。

技术实现要素：

本发明针对现有把粒子滤波应用于卫星导航接收机码跟踪中所存在的滤波估计精度不足的问题，提供一种抑制多径干扰的码跟踪环路及方法。

为解决上述问题，本发明是通过以下技术方案实现的：

一种抑制多径干扰的码跟踪方法，包括步骤如下：

步骤1、对接收得到中频信号与本地载波信号进行混频后，得到同相和正交信号；

步骤2、将同相和正交信号均分为n条同相和正交支路信号，每条同相和正交支路信号分别与伪码发生器所产生的本地伪码进行相干积分计算后，得到n条相干积分结果并缓存；

步骤3、利用最小二乘原理对n条相干积分结果进行修正处理后，得到一组无偏、且均方误差和最小的修正观测量，并将其作为后续粒子滤波的量测信息；

步骤4、把无迹变换的思想引入粒子滤波算法中，并结合利用最小二乘原理修正得到的量测信息对伴随有多径干扰的码环进行跟踪，得到直达信号与本地伪码的码相位差；

步骤5，将得到的码相位差进行滤波后反馈至本地伪码数控振荡器中，伪码数控振荡器根据码相位差相应的调节输出频率，从而调节伪码发生器的本地伪码的生成，直到直达信号与本地伪码的码相位差收敛于零值附近。

上述步骤3的具体子步骤为：

步骤3.1、将各条支路的相干积分结果写成矩阵形式并建立测量方程；

f＝m·g+p(t)

步骤3.2、对测量方程中的状态估计矩阵g进行最小二乘估计，得到估计矩阵

步骤3.3、根据估计矩阵得到一组无偏、且均方误差和最小的修正观测量，并将其作为后续粒子滤波的量测信息；

其中，f为缓存同一条相关支路的相干积分值的矩阵；m是通过载波环构成的量测矩阵；p(t)是噪声矢量；g为状态估计矩阵，g＝[g1,g2]^t；是由最小二乘方法估计得到的矩阵。

上述步骤4的具体子步骤为：

步骤4.1、粒子初始化，即根据初始状态x0的先验概率密度p(x0)生成粒子的初始值并将其作为当前粒子集合；

步骤4.2、对当前粒子集合采用无迹卡尔曼滤波算法进行重要性权重的计算，得到重要性权值

步骤4.3、对步骤4.2所得到的重要性权值进行归一化，得到归一化权重

步骤4.4、计算衡量退化程度的指标并将衡量退化程度的指标与预设的门限值nthr进行比较：当时，将步骤4.3得到的归一化权值及所对应的粒子集合映射为一个等权重的粒子集合且把该等权重的粒子集合作为步骤4.1中的当前粒子集合，并转至步骤4.2；否则，转至步骤4.5；

步骤4.5、利用步骤4.3中得到的归一化权重估计直达信号与本地伪码的码相位差

其中，表示通过引入无迹变换的粒子滤波计算得到的归一化权重；表示直达信号与本地伪码的码相位差的估计值；i＝1,2,…,n，n为相干积分修正支路的条数。

实现上述方法的一种抑制多径干扰的码跟踪环路，包括本地载波模块、混频器、伪码发生器、伪码数控振荡器、环路滤波器、无迹粒子滤波模块、相关器、积分累加器、缓存模块和最小二乘模块；

上述相关器、积分累加器、缓存模块和最小二乘模块的数量相同，且均为n个，n为大于1的正整数；1个相关器、1个积分累加器、1个缓存模块和1个最小二乘模块组成1条相干积分修正支路，其中相关器的2个输入端分别形成该条相干积分修正支路的2个输入端，相关器的输出端连接积分累加器的输入端，积分累加器的输出端连接缓存模块的输入端，缓存模块的输出端连接最小二乘模块的输入端，最小二乘模块的输出端形成该条相干积分修正支路的输出端；

外部输入的中频数据连接混频器的一个输入端，混频器的另一个输入端连接本地载波模块的输出端；混频器的输出端分为n路，每一路分别连接1条相干积分修正支路的一个输入端；n条相干积分修正支路的输出端分别与无迹粒子滤波模块的n个输入端相连，无迹粒子滤波模块的输出端连接环路滤波器的输入端，环路滤波器的输出端连接伪码数控振荡器的输入端，伪码数控振荡器的输出端连接伪码发生器的输入端，伪码发生器的输出端分为n路，每一路分别连接1条相干积分修正支路的另一个输入端。

与现有技术相比，本发明通过最小二乘算法对i/q支路相干积分值进行修正并作为粒子滤波的观测量，同时为提高滤波的估计精度，采用无迹转换的思想取代粒子滤波算法中仅用转移密度函数产生粒子的方式，粒子的产生过程中合理引入了由最小二乘原理修正得到的当前时刻的量测信息，从而引导粒子向高似然区域移动，最小二乘原理和无迹变换思想二者的结合能够有效减小有偏相干积分值所造成的码相位跟踪误差，提高弱信号环境下环路的跟踪性能。

附图说明

图1是一种抑制多径干扰的码跟踪环路的架构示意图。

图2是无迹粒子滤波的原理简述框图。

图3是一种抑制多径干扰的码跟踪方法的流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实例，并参照附图，对本发明进一步详细说明。

针对弱信号环境下，传统接收机所采用的粒子滤波码环跟踪算法，把有偏的相干积分值直接作为滤波所需要的观测值，但在这种滤波估计中，作为唯一信息源的观测值本身带有噪声误差、系统误差、观测误差、设备误差等各种误差，这些误差势必会影响后续滤波的精度的问题，进而导致传统接收机跟踪精度低，跟踪门限高的问题。本发明提出一种基于改进的无迹粒子滤波的码跟踪环路及方法，利用最小二乘方法修正每一条i/q支路的相干积分值，得到一组无偏，均方误差和最小的修正观测量，减小了由于有偏相干积分值给粒子滤波器带来的误差，再应用无迹粒子滤波原理构建数学模型对码环的码相位进行最优估计，即应用无迹转换的思想取代粒子滤波中仅利用重要性采样密度采样粒子的方式，把使用最小二乘原理修正得到的观测信息引入到重要性密度函数的生成里，粒子不再仅依赖转移密度函数生成，而它的产生更多的考虑进去了通过最小二乘原理修正得到的相对优质的量测信息，最小二乘原理与无迹变换思想二者的结合进而实现更高精度的码相位跟踪。

本发明在传统的粒子滤波算法前引入最小二乘算法来修正有偏的相干积分值，修正得到一组无偏且均方误差和最小的值作为后续粒子滤波的观测信息，同时为减少粒子退化的影响，本发明采用ukf(无迹卡尔曼滤波)来计算粒子滤波中的重要性权重，改进后的码跟踪环架构，即一种抑制多径干扰的码跟踪环路，如图1所示，包括本地载波模块、混频器、伪码发生器、伪码nco(伪码数控振荡器)、环路滤波器、无迹粒子滤波模块、相关器、积分累加器、缓存模块和最小二乘模块。上述相关器、积分累加器、缓存模块和最小二乘模块的数量相同，且均为n个，n为大于1的正整数。1个相关器、1个积分累加器、1个缓存模块和1个最小二乘模块组成1条相干积分修正支路，其中相关器的2个输入端分别形成该条相干积分修正支路的2个输入端，相关器的输出端连接积分累加器的输入端，积分累加器的输出端连接缓存模块的输入端，缓存模块的输出端连接最小二乘模块的输入端，最小二乘模块的输出端形成该条相干积分修正支路的输出端。外部输入的中频数据连接混频器的一个输入端，混频器的另一个输入端连接本地载波模块的输出端；混频器的输出端分为n路，每一路分别连接1条相干积分修正支路的一个输入端；n条相干积分修正支路的输出端分别与无迹粒子滤波模块的n个输入端相连，无迹粒子滤波模块的输出端连接环路滤波器的输入端，环路滤波器的输出端连接伪码nco的输入端，伪码nco的输出端连接伪码发生器的输入端，伪码发生器的输出端分为n路，每一路分别连接1条相干积分修正支路的另一个输入端。接收中频信号与本地载波模块混频后得到同相、正交两路(即i、q两路)信号，然后与各相关支路的本地伪码做相关累加运算得到相干积分值，相干积分值缓存于缓存模块中，然后输入最小二乘模块做相干积分值修正，修正后的相干积分值作为粒子滤波器的量测信息，并采用ukf计算重要性权重，不断滤除权重低的粒子，最后估计得到直达信号与本地伪码的码相位差。

为了能够准确地跟踪输入信号的码相位差和幅值，粒子滤波器的状态向量表示为：

x＝[τ0τ1α1]^t(1)

式(1)中，τ0为接收卫星信号中直达信号与本地伪码的码相位差，τ1为多径信号与本地伪码的码相位差，α1为多径信号相对直达信号的幅度比例。

将i、q两路经过最小二乘算法处理过的无偏的相干积分值作为滤波器的观测信息量zk，并把它表示成一个与上述状态向量有关的方程为：

式(2)中rk为量测噪声。由上式可以得知，把通过最小二乘方法修正过的观测量与待估计的状态向量建立了一种非线性关系，下面依据修正过的观测信息zk，采用粒子滤波算法对码环环路进行滤波处理，估计得出直达信号、多径信号的码相位差τ0、τ1，以及多径信号的幅度比例α1。

改进后的多径抑制码跟踪环在很大程度上克服了传统粒子滤波对码相位跟踪时存在的跟踪门限低，跟踪精度低的缺陷，同时成功地解决了低信噪比环境下，跟踪不稳定的问题。

上述码跟踪环路所实现的一种抑制多径干扰的码跟踪环路及方法，如图3所示，包括如下几个步骤：

步骤1，接收机天线接收得到的信号经过下变频和采样步骤的处理后，信号被转化成数字中频信号，再将数字中频信号信号分别与正弦、余弦复制载波相乘得到i、q两路，最后将其送入到伪码相关器中与本地伪码进行相干积分计算，得到各支路相干积分的结果并缓存。

步骤1-1，接收机天线接收得到的信号经过下变频和采样步骤的处理后，信号被转化成数字中频信号。

步骤1-2，将上一步骤得到的数字中频信号分别与正弦、余弦复制载波相乘得到i、q两路，然后将其送入到伪码相关器中与本地伪码进行相干积分计算。

步骤1-3，在缓存模块中缓存各支路的相干积分值。

步骤2，利用最小二乘原理对有偏的相干积分值进行修正处理，把修正后的各条相关支路的相干积分值组合起来构成后续粒子滤波的量测信息。

步骤2-1，接收机在对数字中频信号去载波和伪码相关处理后得到的第j条相关支路的相干积分结果为：

其中，τ0和τ1分别是直达信号和多径信号的码相位时延；α1是多径信号相对于直达信号的幅值比例，a是本地载波幅度；a是直达信号载波幅度；fe是载波频率差；是直达信号与本地载波的相位差；是多径信号与本地载波的相位差；p(t)是服从n(0,σ²)的高斯白噪声；rc(·)是伪码的自相关函数。这里的fe、均是由载波跟踪环得到。

步骤2-2，把上述相干积分结果写成矩阵形式并建立测量方程。上述相干积分结果表示成矩阵形式如下：

把上式写成矢量形式，可建立如下测量方程：

f＝m·g+p(t)(5)

其中，f为缓存同一条相关支路的相干积分值的矩阵；m是通过载波环得到的量测矩阵；p(t)是噪声矢量；g视为状态估计矩阵，g＝[g1,g2]^t。

g1＝aarc(δj+τ0)

g2＝α1aarc(δj+τ1)

其中，α1是多径信号相对于直达信号的幅值比例，a是本地载波幅度；a是直达信号载波幅度；rc(·)是伪码的自相关函数；τ0为接收卫星信号中直达信号与本地伪码的码相位差，τ1为多径信号与本地伪码的码相位差。

步骤2-3，根据最小二乘估计原理，得到状态矩阵g的估计值公式如下：

步骤2-4，根据每一条相关支路的状态估计矩阵得到一组修正的观测值，并作为粒子滤波的观测矢量。

最小二乘估计的目的是从中获得无偏、均方误差和最小的估计值，也即得到直达信号和多径信号的相干积分的无偏、均方误差和最小的估计值，并用该估计值来修正观测值f(t)，这里把通过最小二乘修正过的观测值记为

上式中的是以i路的第j条超前支路的修正相干积分值为例说明的，同法可求i路的第j条滞后支路的修正相干积分值q路的第j条超前支路的修正相干积分值q路的第j条滞后支路的修正相干积分值把每一相关支路修正后的相干积分值集合起来组成一组无偏、均方误差和最小的修正观测矢量zk，用于后续粒子滤波的估计。

上式中的n表示某种相干支路的总条数。

步骤3，把无迹变换的思想引入粒子滤波算法中，应用无迹粒子滤波算法并结合修正后的量测信息对伴随有多径干扰的码环进行跟踪，得到直达信号与本地伪码的高精度码相位差。无迹粒子滤波算法如图2所示。

粒子滤波算法是应用序贯蒙特卡罗方法用粒子实现后验概率密度计算的过程，也即利用一组带权重的随机样本去构建后验概率密度函数(pdf)，以样本均值代替积分运算，进而获得状态的最小方差估计的算法。为了更清楚地应用数学语言表示该算法的实现过程，粒子滤波器的非线性系统模型定义为：

xk＝fk-1(xk-1)+wk(9)

zk＝hk(xk)+vk(10)

式(9)表示状态方程，xk的定义如式(1)所示，fk-1(·)是状态xk-1的非线性函数，wk是过程噪声且状态转移符合一阶markov过程。式(10)表示量测方程，hk(·)是状态xk的非线性函数，vk是量测噪声序列。粒子滤波的目的就是利用含有噪声的量测信息估计目标状态xk，也即式(10)描述的方程结合式(9)递推地估计系统状态xk。

粒子滤波中，后验概率密度函数p(xk/zk)计算的准确程度直接关系到整个算法的精确度。k时刻的后验概率密度函数p(xk/zk)的粒子表示如下：

式(11)中，表示k时刻第个粒子的权重，δ(·)是狄利克雷函数，xk是k时刻的目标状态，是k时刻的第i个粒子，n表示用于估计目标的粒子的总个数。已知k-1时刻后验概率密度p(xk-1/zk-1)的条件下，k时刻的后验密度函数的求解一般需要经历预测和更新两个步骤，下面给出预测、更新方程：

p(xk/zk-1)＝∫p(xk/xk-1)·p(xk-1/zk-1)dxk-1(12)

其中：

p(zk/zk-1)＝∫p(zk/xk)·p(xk/zk-1)dxk(14)

式(12)、(13)中，p(xk/zk-1)是k时刻预测的后验概率密度函数，p(xk/xk-1)是转移密度，p(zk/xk)表示似然函数，p(xk/zk)是更新的k时刻的后验概率密度函数。如何把(12)、(13)、(14)式的数学理论公式求解后验概率密度函数p(xk/zk)用粒子的形式表示成式(11)，下面给出粒子滤波的具体实现步骤：

步骤3-1，初始化。

根据初始状态x0的先验概率密度p(x0)生成粒子的初始值

步骤3-2，粒子采集。

粒子集合是从重要性密度中采样获取，为了方便一般把转移密度作为重要性密度进行样本的采集。

步骤3-3，采用ukf算法进行重要性权重的计算。

步骤3-4，权重归一化。

步骤3-5，重采样。

由于重要性函数的方差是随着时间的增加而增加的，而方差的增加会给精确估计带来影响，导致样本枯竭问题，称之为粒子退化。衡量退化程度的指标定义为

式(18)中，值越小，粒子退化情况越严重，为了粒子防止退化现象，须进行重采样步骤。重采样也即舍弃权重小的粒子，复制权重大的粒子。为了说明重采样步骤，这里设置一个门限值nthr。若时，进行重采样，将原来的带权粒子集合映射为等权重的粒子集合

步骤3-6，输出码相位差等状态估计值。

步骤4，将得到的直达信号与本地伪码的码相位差滤波后反馈输回本地伪码数控振荡器中，伪码数控振荡器根据码相位差相应的调节输出频率，从而调节伪码发生器的本地伪码的生成。

重复步骤1-4，直到直达信号与本地伪码的码相位差收敛于零值附近。

由上述实施例可以看出，本发明提出的一种抑制多径干扰的改进码跟踪环路，把经过最小二乘原理修正过后的无偏量测值应用于粒子滤波算法中，并通过采用ukf来计算粒子滤波中的重要性采样密度，采用一系列确定的粒子集合来逼近状态的后验概率密度，提高了粒子滤波的精度和稳定性，实现了在低信噪比条件下的多径抑制，从而提高了码相位的跟踪精度和稳定性。

需要说明的是，尽管以上本发明所述的实施例是说明性的，但这并非是对本发明的限制，因此本发明并不局限于上述具体实施方式中。在不脱离本发明原理的情况下，凡是本领域技术人员在本发明的启示下获得的其它实施方式，均视为在本发明的保护之内。