引言

在全球导航卫星系统逐步发展的过程中，对其性能的需求在不断提升。然而，由于卫星导航信号的脆弱性，导致其易受到干扰影响。其中多径干扰是指来自卫星的信号, 被接收机天线周围物体反射后，与直达卫星信号一起进入天线，从而对接收机形成干扰。多径干扰导致的直接后果是产生跟踪和测量误差, 降低定位精确度。相关研究指出，多径干扰导致的伪距误差可以达到米级甚至百米级，从而严重危及到系统的可靠性和定位精度。因此, 多径干扰抑制技术一直是卫星导航接收机设计领域的研究热点^[1-5]。

窄带相关技术通过减小早晚相关器间隔来降低多径对码跟踪环路的影响，是目前较为常用的一种多径干扰抑制技术。Michael S B对窄带相关器模型进行了分析^[6]，Cannon M E等分析了窄带相关器技术的性能^[7]，Novatel公司的Dierendonck等对其理论和性能进行了较为详细的阐述分析^[8]。理论上，相关器间隔越小，对多径干扰环境下的误差改善越好。然而，这是基于接收机信号通道为无限带宽，而这个假设条件在实际中无法得到满足。当相关器间隔小于通道等效带宽的倒数时，码相位跟踪误差将趋于一个常数，因此不能通过无限降低相关器间隔的方法来减小码相位跟踪误差。

多径估计延迟锁定环技术(Multipath estimating delay lock loop, MEDLL)是NovAtel公司1992年提出的多相关器接收机技术^[9]。MEDLL的核心思想是其接收机中的多个相关器获得相关函数的多个采样值，进而根据最大似然准则进行迭代优化，这使得较之于窄带相关技术，MEDLL的复杂度更高。本团队提出的基于WRELAX算法是通过最小化非线性最小二乘(None-linear least square, NLS)代价函数的估计方法，它逐个估计出直达和多径信号的时延信息，然而该方法同样存在计算量大的问题^[10]。为了进一步减小算法运算量，本文在WRELAX算法的基础上，结合卫星导航信号相关函数的特点，对接收数据与本地参考信号的互相关函数、本地参考信号的自相关函数进行加窗截取，取出相关函数中信息量较大的部分，即相关函数主峰值附近的范围，从而大大减小了WRELAX参数估计中迭代处理的数据长度，也就大大减少了算法的运算量，便于工程化。

1 数据模型

以GPS卫星信号为例，其结论可直接推广到其他GNSS系统。卫星发射信号可表示为

$ x\left( t \right) = d\left( t \right)c\left( t \right)\cos \left( {{\omega _c}t} \right) $

(1)

式中：d(t)为导航电文，c(t)为C/A码。考虑接收机可收到P条反射的多径干扰，此时接收到的信号模型可以表示为

$ \bar y\left( t \right) = \sum\limits_{p = 0}^P {{{\alpha '}_p}d\left( {t - {\tau _p}} \right)c\left( {t - {\tau _p}} \right)\cos \left( {{\omega _c}\left( {t - {\tau _p}} \right) + {{\varphi '}_p}} \right) + \bar e\left( t \right)} $

(2)

式中：α′ _p, τ_p, φ′ _p分别为第p条路径信号的幅度、时延和初相位，其中p=0表示直达卫星信号，$ \bar e\left( t \right)$表示接收机热噪声。由于导航电文的重复周期远大于C/A码的码重复周期，而在时延估计过程中需要的数据块长度一般都为一个C/A码周期，因此在信号重构信号过程中可忽略导航电文跳变产生的影响，此时的数据可表示为

$ \bar y\left( t \right) = \sum\limits_{p = 0}^P {{{\alpha '}_p}c\left( {t - {\tau _p}} \right)\cos \left( {{\omega _c}\left( {t - {\tau _p}} \right) + {{\varphi '}_p}} \right) + \bar e\left( t \right)} $

(3)

由于卫星距地面较远，可认为多径与直达信号具有相同的多普勒频率ω_d0^[10]，进一步化简式(3)得

$ \bar y\left( t \right) = \sum\limits_{p = 0}^P {{{\alpha '}_p}c\left( {t - {\tau _p}} \right)\cos \left( {{\omega _{d0}}t + {\varphi _p}} \right) + \bar e\left( t \right)} $

(4)

式中：φ₀=φ′ ₀+2πR₀/λ为直达信号相位，φ_p=φ₀+Δφ_p为多径干扰相位，Δφ_p为第p条多径相对于直达信号的额外相位差。假设载波恢复完全准确，那么经过载波环路解调后的信号可表示为

$ \bar {\bar y}\left( t \right) = \sum\limits_{p = 0}^P {{{\alpha '}_p}c\left( {t - {\tau _p}} \right)\cos \left( {{\varphi _p}} \right) + \bar {\bar e}\left( t \right)} $

(5)

因此，本文处理的目的是从接收到的混合数据中得到直达卫星信号的准确参数估计。

2 本文算法

基于非线性最小二乘准则的WRELAX算法^[11-14]能够将多维非线性寻优搜索问题转化为一系列一维搜索，从而大大减少了算法的运算量。与现有的时延估计算法相比，WRELAX更为系统和有效，且对信号波形没有限制。下面介绍基于WRELAX的GNSS多径干扰抑制算法^[10]。若对载波解调后的式(5)进一步作希尔伯特变换，得到对应的复数据可表示为

$ y\left( t \right) = \sum\limits_{p = 0}^P {{{\alpha '}_p}c\left( {t - {\tau _p}} \right){{\rm{e}}^{{\rm{j}}{\varphi _p}}} + e\left( t \right)} $

(6)

进一步对式(6)数据以T_s为间隔进行采样，可得

$ y\left( {n{T_{\rm{s}}}} \right) = \sum\limits_{p = 0}^P {{{\alpha '}_p}c\left( {n{T_{\rm{s}}} - {\tau _p}} \right)\exp \left( {{\rm{j}}{\varphi _p}} \right) + e\left( {n{T_{\rm{s}}}} \right)} \;\;\;\;\;n = 1,2, \cdots ,N $

(7)

其中N为总的采样点数。若将式(7)中相位exp(jφ_p)与幅度α′ _p合计为一个新量，并用α_p来表示该量，则式(7)可表示为

$ y\left( {n{T_{\rm{s}}}} \right) = \sum\limits_{p = 0}^P {{\alpha _p}c\left( {n{T_{\rm{s}}} - {\tau _p}} \right) + e\left( {n{T_{\rm{s}}}} \right)} $

(8)

通过准确估计上述模型中的时延信息来抑制多径干扰。

为了得到直达信号和多径信号的初始码相位，引入相关函数

$ r\left( {n{T_{\rm{s}}}} \right) = {\rm{corr}}\left( {c\left( {n{T_{\rm{s}}}} \right),y\left( {n{T_{\rm{s}}}} \right)} \right) $

(9)

式中corr(·)代表相关运算。因为信号和噪声之间是不相关的，所以式(9)可以进一步表示为

$ r\left( {n{T_{\rm{s}}}} \right) = \sum\limits_{p = 0}^P {{\alpha _p}{r^c}\left( {n{T_{\rm{s}}} - {\tau _p}} \right) + w\left( {n{T_{\rm{s}}}} \right)} $

(10)

其中r^c(nT_s)为参考信号c(nT_s)的自相关函数，w(nT_s)=corr (c(nT_s), e(nT_s))为噪声与参考信号互相关的结果。

对于GPS卫星信号来说，C/A码的波形如图 1所示，在整个C/A码的持续时间内，信号幅值不断变化，也就是C /A码的信息等概率分布在整个持续时间内。C/A码的相关函数如图 2所示，可以发现其相关函数的信息量集中在±1个码片T_c范围内。因此，考虑对接收信号与本地参考信号的相关函数以及参考相关函数进行加窗截取处理，即只截取相关函数信息量较大的部分。由于当多径干扰与直达信号延迟超过1.5个码片时将不对码跟踪环路产生影响，故在截取相关函数时截取范围取±2个码片。

对式(10)加窗处理后可进一步表示为

$ g\left( {n{T_{\rm{s}}}} \right) = \left[ {\sum\limits_{p = 0}^P {{\alpha _p}{{\bar r}^c}\left( {n{T_{\rm{s}}} - {\tau _p}} \right)} } \right] + \bar w\left( {n{T_{\rm{s}}}} \right) $

(11)

其中${\bar r^c}\left( {n{T_{\rm{s}}}} \right) $对应于式(10)中r^c(nT_s－τ_p)加窗后的结果，$\bar w\left( {n{T_{\rm{s}}}} \right) $窗函数加权后的噪声与参考信号互相关的结果。对式(11)做DFT得

$ {g_f}\left( k \right) = \left[ {\bar r_f^c\left( k \right)\sum\limits_{p = 0}^P {{\alpha _p}{{\rm{e}}^{{\rm{j}}{\omega _p}k}}} } \right] + \bar w\left( k \right) $

(12)

其中g_f(k)，$ \bar r_f^c\left( k \right)$和$\bar w\left( k \right) $分别为$g\left( {n{T_{\rm{s}}}} \right) $,${{\bar r}^c}\left( {n{T_{\rm{s}}}} \right) $和$\bar w\left( {n{T_{\rm{s}}}} \right) $对应的离散傅里叶变换。

定义NLS代价函数为

$ {C_1}\left( {\left\{ {{\omega _p},{\alpha _p}} \right\}_{p = 0}^P} \right) = \sum\limits_{k = - N/2}^{N/2} {{{\left\| {{g_f}\left( k \right) - \left[ {\bar r_f^c\left( k \right)\sum\limits_{p = 0}^P {{\alpha _p}{{\rm{e}}^{{\rm{j}}{\omega _p}k}}} } \right]} \right\|}^2}} $

(13)

通过最小化代价函数$ {C_1}\left( {\left\{ {{\omega _p}, {\alpha _p}} \right\}_{p = 0}^P} \right)$可得到未知参数的估计值。这里采用WRELAX算法来求解这个NLS代价函数，假设截取后的数据长度为N′，并令

$ \mathit{\boldsymbol{a}}\left( {{\omega _p}} \right) = {\left[ {{{\rm{e}}^{{\rm{j}}{\omega _p}\left( { - N'/2} \right)}},{{\rm{e}}^{{\rm{j}}{\omega _p}\left( { - N'/2 + 1} \right)}}, \cdots ,{{\rm{e}}^{{\rm{j}}{\omega _p}\left( {N'/2 - 1} \right)}}} \right]^{\rm{T}}} $

(14)

$ {\mathit{\boldsymbol{g}}_f} = {\left[ {{g_f}\left( { - N'/2} \right),{g_f}\left( { - N'/2 + 1} \right), \cdots ,{g_f}\left( { - N'/2 - 1} \right)} \right]^{\rm{T}}} $

(15)

$ {{\mathit{\boldsymbol{\bar R}}}_f} = {\rm{diag}}\left[ {\bar r_f^c\left( { - N'/2} \right),\bar r_f^c\left( { - N'/2 + 1} \right), \cdots ,\bar r_f^c\left( {N'/2 - 1} \right)} \right] $

(16)

那么式(13)的代价函数可以重新表示为

$ {C_2}\left[ {\left\{ {{\omega _p},{\alpha _p}} \right\}_{p = 0}^P} \right] = {\left\| {{\mathit{\boldsymbol{g}}_f} - \sum\limits_{p = 0}^P {{\alpha _p}{{\mathit{\boldsymbol{\bar R}}}_f}\mathit{\boldsymbol{a}}\left( {{\omega _p}} \right)} } \right\|^2} $

(17)

进一步假设$\left\{ {{{\hat w}_q}, {{\hat \alpha }_q}} \right\}_{q = 1, q \ne p}^P $已经被估计出来，则有

$ {\mathit{\boldsymbol{g}}_{fp}} = {\mathit{\boldsymbol{g}}_f} - \sum\limits_{\begin{array}{*{20}{c}} {q = 0}\\ {q \ne p} \end{array}}^P {{{\hat \alpha }_q}\left[ {{{\mathit{\boldsymbol{\bar R}}}_f}\mathit{\boldsymbol{a}}\left( {{{\hat \omega }_p}} \right)} \right]} $

(18)

将式(18)代入上面的代价函数，可以得到

$ {C_3}\left[ {\left( {{\omega _p},{\alpha _p}} \right)_{p = 0}^P} \right] = {\left\| {{\mathit{\boldsymbol{g}}_{fp}} - {\alpha _p}{{\mathit{\boldsymbol{\bar R}}}_f}\mathit{\boldsymbol{a}}\left( {{\omega _p}} \right)} \right\|^2} $

(19)

最小化代价函数C₃，分别得到$ {{{\hat w}_p}}$和$ {{{\hat \alpha }_p}}$的估计值为

$ {{\hat \omega }_p} = \arg \mathop {\max }\limits_{{\omega _p}} {\left| {{\mathit{\boldsymbol{a}}^{\rm{H}}}\left( {{\omega _p}} \right)\left( {\mathit{\boldsymbol{\bar R}}_f^ * {\mathit{\boldsymbol{g}}_{fp}}} \right)} \right|^2} $

(20)

且

$ {{\hat \alpha }_p} = \frac{{{\mathit{\boldsymbol{a}}^{\rm{H}}}\left( {{\omega _p}} \right)\left( {\mathit{\boldsymbol{\bar R}}_f^ * {\mathit{\boldsymbol{g}}_{fp}}} \right)}}{{\left\| {{{\mathit{\boldsymbol{\bar R}}}_f}} \right\|_{\rm{F}}^2}}\left| {_{{\omega _p} = {{\hat \omega }_p}}} \right. $

(21)

因此，${{{\hat w}_p}} $可以通过$|\mathit{\boldsymbol{a}}{^{\rm H}}\left( {{{\hat w}_p}} \right)\left( {\mathit{\boldsymbol{\bar R}}_f^ * \mathit{\boldsymbol{g}}{_{fp}}} \right){|^2} $的最大值所对应的位置求得，而$ |\mathit{\boldsymbol{a}}{^{\rm H}}\left( {{{\hat w}_p}} \right)\left( {\mathit{\boldsymbol{\bar R}}_f^ * \mathit{\boldsymbol{g}}{_{fp}}} \right){|^2}$实际上可以通过对$ {\mathit{\boldsymbol{\bar R}}_f^ * \mathit{\boldsymbol{g}}{_{fp}}}$进行补零的FFT来计算。而$ {{{\hat \alpha }_p}}$可通过计算$ \mathit{\boldsymbol{a}}{^{\rm H}}\left( {{{\hat w}_p}} \right)\left( {\mathit{\boldsymbol{\bar R}}_f^ * \mathit{\boldsymbol{g}}{_{fp}}} \right)/\parallel {{\mathit{\boldsymbol{\bar R}}}_f}\parallel _{\rm{F}}^2$的复幅度得到。至此，可以发现WRELAX将最小化NLS代价函数对应的多维搜索问题转换为一系列一维搜索，并进一步通过一维FFT实现参数的估计。整个算法基本流程如图 3所示。

综上所述，卫星信号和多径的参数$\left\{ {{{\hat w}_p}, {{\hat \alpha }_p}} \right\}_{p = 0}^1 $估计步骤如下：

(1) 假设一路信号，即认为不含多径干扰，只有直达卫星信号的情况，此时P=0，按式(20, 21)估计出$\left\{ {{{\hat w}_0}, {{\hat \alpha }_0}} \right\} $。

(2) 假设有两路信号P=1(直达信号和一路多径信号)，用步骤(1)得到的$\left\{ {{{\hat w}_0}, {{\hat \alpha }_0}} \right\} $和${{\mathit{\boldsymbol{\bar R}}}_f} $根据式(18)得到g_f1，接着利用式(20, 21)估计出$\left\{ {{{\hat w}_1}, {{\hat \alpha }_1}} \right\} $，即为第2条信号的估计值。然后再利用$\left\{ {{{\hat w}_1}, {{\hat \alpha }_1}} \right\} $和式(18)减去第2条信号，重新得到g_f0，再用g_f0重新确定$\left\{ {{{\hat w}_0}, {{\hat \alpha }_0}} \right\} $。

重复以上过程，直至收敛，这样就可以求出估计值$ \left\{ {{{\hat w}_p}, {{\hat \alpha }_p}} \right\}_{p = 0}^1$。进一步利用关系式ω_p=－2πτ_p/NT_s，可以求出信号的时延估计值$ {{\hat \tau }_p}$。

由于GPS C/A码的总长度为1 023个码片，而对相关函数加窗后其数据长度与窗函数的长度有关，若选取以相关函数峰值为中心的±2个码片范围，那么后续参数估计过程中的数据长度将会变成原相关域WRELAX算法的4/1 023，运算量也必然减少为原来的4/1 023，因此本节介绍的加窗后处理的WRELAX算法运算量大大减小，从而为工程实现提供了有利条件。

3 仿真实验 3.1 加窗的影响分析实验

仿真产生码延迟为0的卫星信号，采样频率为5 MHz，计算该卫星信号与本地参考信号之间的互相关函数如图 4所示。可以看出，经过加窗后，相关函数的长度大大减小，但是在窗函数范围内的相关函数主峰值形状保持不变，窗函数之外的相关函数其余部分会有信息损失，但是这些对于参数估计不起决定作用。

加窗后前后WRELAX算法的代价函数如图 5所示，可以看出，加窗前后代价函数的主瓣相当，即加窗本身对参数估计不会产生严重影响，同上面的结论一致，虽然加窗使数据长度缩短，但是其含有信息的部分却完全保留。

3.2 各种算法性能对比实验

用卫星导航软件接收机处理仿真生成卫星信号和一路多径信号。设置直达信号与多径干扰的幅值比(Direct to multipath ratio, DMR)分别取0.1, 0.3, 0.5, 0.7，相对延迟τ₂－τ₁从0~1.5个码片，无噪声、无限带宽条件下，分别使用传统相关器、窄相关器、MEDLL和本文所提两种基于信号分离的时延估计算法得到的多径误差包络如图 6所示。

从图 6的实验结果可以看出，相关器间隔d=1时传统相关器存在较大的时延估计误差，随着相关器间隔的缩小，时延估计误差也逐渐减小。MEDLL技术、相关模型WRELAX及加窗后处理的算法作为高分辨率的时延估计方法能提供比窄带相关器更好的性能。为了进一步观察这3种基于时延估计方法的性能，对上述实验结果进行放大得到图 7。从图 7中可以看出，当DMR=0.1~0.3时，这些方法的性能相当，也就是说，多径干扰比较弱时(实际情况中经过反射后的多径干扰都会产生严重的衰减)，加窗后处理的方法在保证能以较低运算量实现的前提下也能保证其性能不受较大损失；当多径干扰逐渐增强，如DMR=0.5时，加窗后处理的方法性能会逐渐劣于相关域WRELAX算法和MEDLL，但是其性能仍优于窄带相关技术。加窗后处理的方法由于没有考虑主峰值之外的部分，因此参数估计的迭代过程会导致这部分的信息有所损失，最终导致其估计性能较之原始的WRELAX算法有所下降。

本文所提算法是从减少运算量的角度，通过缩减参数搜索范围对文献[15]所提WRELAX算法的改进，文献[15]给出了噪声情况下WRELAX的性能分析，证明了相对于窄相关和MEDLL等常用的多径抑制技术，WRELAX算法具有良好的性能。通过理论分析和仿真实验证明了所提降低运算量的算法与原来的WRELAX算法性能相当，由于新算法只是巧妙地利用卫星信号自身良好的自相关特性来降低参数估计的搜索区间，因此不会影响算法的抗噪声性能。

4 结束语

本文根据GNSS信号相关函数的特点，提出了一种加窗处理的相关域WRELAX多径干扰抑制算法。该算法通过对接收数据与本地参考信号的相关函数以及本地参考信号的自相关函数加窗处理，使得参与参数估计的相关函数长度大大减小，而由于相关函数的信息量主要集中在相关峰值附近的范围内，因此加窗后相关函数的信息量不会有明显损失。最后，通过仿真实验进一步验证了本文算法能够在以较低运算量实现的前提下，同时保证性能不会受到太大损失，尤其是多径干扰较弱时。