CN102213766B - 于卫星导航接收器中抑制多重路径误差的方法与装置 - Google Patents

Abstract

本发明公开了于卫星导航接收器中抑制多重路径误差的方法与装置，包含搜寻导航卫星信号、由卫星接收并处理复变信号、以一追随开回路追踪信号并基于直接与反射信号的延迟与多普勒频率测量值决定接收器坐标、速度与准确时间等步骤，而该方法包含：根据导航参数形成一二维的累积功率网格、利用中心化的累积功率网格计算单一路径信号修正值，并决定多重路径是否发生，若是则利用权重改进修正值。利用权重改进修正值时，基于一中心化的累积功率向量形成一概度函数，并根据该概度函数的全域极大值改进修正值。

Description

于卫星导航接收器中抑制多重路径误差的方法与装置

技术领域

本发明关于一种卫星无线导航，特别是关于一种在全球导航卫星系统(GlobalNavigating Satellite System，GNSS)的导航接收器中求得多重路径情况下的准距离(pseudo-range)与准速度(pseudo-velocity)的方法与装置。

背景技术

目前，世界上已经存在数种全球导航卫星系统。已上线运作的系统有美国的全球定位系统(Global Positioning System，GPS)，俄罗斯的全球导航卫星系统(Global Navigating Satellite System，GLONASS)，欧洲太空总署建构中的伽利略(Galileo)系统，以及中国的北斗系统。本发明可应用于上述或其他全球导航卫星系统的导航接收器。

在都市高楼形成的峡谷中，会导致导航接收器频繁地失去与卫星间的视线(line-of-sight)路径，接收信号的功率也因此而变弱或消失。在此类情况下，传统中利用延迟锁定回路与频率锁定回路的闭回路(close-loop)自动信号追踪常会受到干扰。

以开回路(open-loop)追踪卫星信号的方法，如美国专利US6633255「Open-Loop Tracking」，则可描述如下：针对每一卫星，根据导航解计算延迟预测与多普勒频率测量；形成一二维(两维度分别为延迟与多普勒频率)功率网格(power grid)，使其中心由预测值定义；针对前述预测求得修正值，并以修正值与预测值的加总作为延迟与多普勒频率的测量值。此种追踪技术可确保低讯杂比下的运作。

当导航解存在时，在开回路追踪技术中失去部份卫星的信号并不会造成重大影响；也就是说，测量仍需建立于至少三到四个卫星。

来自直接路径与反射路径的信号会导致延迟与多普勒频率的误差。图1示意多重路径信号传播影响测量的例子：如图1a所示，针对一个静止的接收器，或是一个沿着与反射面平行方向移动的接收器，只会发生延迟测量误差，因为由直接路径与反射路径至同一卫星的距离变化率是相同的。

如图1b所示，若接收器沿指向反射面的方向移动，延迟测量的误差会被加总至多普勒频率测量的误差，因为由直接路径与反射路径至同一卫星的距离变化率会是相异的。

已有一些方法用以在反射路径存在的情况下降低准距离测量的误差，如：运用窄相关(Theory and Performance of Narrow Correlator Spacing in a GPSreceiver.A.J.van Dierendonck，Pat Fenton，Tom Ford.Journal of TheInstitute of Navigation，Vol.39，No.3，1992.)、根据准乱数序列(pseudo-random sequence，PRN)的信号副本计算自相关函数(autocorrelationfunction，ACF)的选通方法(美国专利US7436356，申请日2006年3月24日)，基于自相关函数的前缘计算跨越前缘与后缘片段并对自相关函数取样加以权重的方法(美国专利US7436356，申请日2006年3月24日)。这些方法可实现于低成本、大量生产的接收器中以降低准距离的误差，但精确度却不够高。

此外，也有些减轻多重路径的方法是针对多重路径的信号进行最佳化运算，以运用最大概度(maximum likelihood)测量直接与反射路径的信号参数。估计最佳值并将反射信号由接收器排除即可实现反射信号的抑制。这些方法应用在高成本的大地接收器，需要昂贵的硬体才能实现。

以适应性天线阵列抑制反射路径信号的技术需要多个接收天线，也需要对这些天线接收的信号进行额外的处理，增加消费设备的复杂程度。

与本发明相近的解决方案见于美国专利US6031881，公开于2000年2月2日。此技术的缺点包括：必须储存输入信号的取样值，必须以高频执行量化(针对GPS，其取样率需达20MHz)，也必须对取样值进行即时(real-time)处理。因此，必须具备高容量存储器与高速信号处理才能实现此种技术。

此习知技术的另一缺点是受限的累积时间。累积时间受限的原因是接收信号会被发送信号中导航讯息的未知数据位元所调变(在GPS与GLONASS中是每间隔20ms发生)，使接收的载波信号功率C与杂讯功率N0间的比率会有下限的限制；在1Hz频带中，此比率C/N0需大于30dB Hz。

发明内容

本发明的目的之一是消除多重路径对延迟与多普勒频率测量的影响，并提出一种能有效应用于弱信号且在实现时不需庞大运算能量的方法。重要技术特征有：提供位元同步下形成的累积功率网格，建立概度函数排除干扰因素以抑制多重路径误差。

本发明解决方案如下。

本发明可为一种于卫星导航接收器(简称接收器)中抑制多重路径误差的方法，卫星导航接收器包含搜寻多颗导航卫星的信号、由各卫星接收并处理复变信号(complex signal)、以一追随开回路追踪信号并基于直接信号与反射信号的延迟与多普勒频率测量值决定该卫星导航接收器的坐标、速度与准确时间的步骤，依据下列描述：

根据多个导航参数计算延迟与多普勒频率的预测值以根据延迟与多普勒频率定向一二维的累积功率网格，此累积功率网格系依据对复变信号的接收与处理所形成，而该些导航参数包括接收器的坐标、速度、准确时间，以及根据测量到的时间与接收到的卫星星历所求得的卫星运动信息。在形成累积功率网格的各元素时，本发明会在一同调累积时段中进行复变信号加总，此一同调累积时段同步于一给定卫星的数据位元，同调累积时段长短符合一位元长度。基于对累积功率网格的处理，包括对累积功率网格的中心化，即可为延迟与多普勒频率的预测值计算修正值，并计算预测值与修正值的加总以作为延迟与多普勒频率的测量值。然后，便可对测量值外插一预设时间并计算出准距离与准速度的测量值，并根据对所有卫星所求得的测量值决定接收器的坐标与速度。

在计算延迟与多普勒频率的修正值时，包含：利用中心化的累积功率网格计算单一路径的修正值；判断多重路径是否发生，若为肯定的，则利用多个权重改进该修正值。利用该些权重改进修正值时，系基于一中心化的累积功率向量形成一概度函数，并根据概度函数的全域极大值改进修正值。

在第一实施例中，在基于中心化的累积功率向量形成概度函数时，包含下列子步骤：利用多个权重向量而为每一累积功率向量计算多个加权总和，使这些权重向量对应直接信号与反射信号的延迟预测值；针对该些权重向量的两两纯量内积所形成的矩阵求出反矩阵；并根据前述两子步骤的结果计算出二次形式的概度函数。

在第二实施例中，当基于中心化的累积功率向量形成概度函数时，包含下列子步骤：针对每一累积功率向量，利用多个取决于直接信号与反射信号延迟的标准正交向量计算多个加权总和，并计算这多个加权总和的平方总和以求出概度函数。

为降低多普勒频率测量误差所形成的概度函数可同理类推。

在第一例中，形成概度函数包含下列子步骤：利用多个权重向量而为每一累积功率向量计算多个加权总和，使这些权重向量对应直接信号与反射信号的多普勒频率预测值；针对该些权重向量的两两纯量内积所形成的矩阵求出反矩阵；并根据前述两子步骤的结果计算出二次形式的概度函数。

在第二例中，形成概度函数包含下列子步骤：针对每一累积功率向量，利用多个取决于直接信号与反射信号多普勒频率的标准正交向量计算多个加权总和，并计算这多个加权总和的平方总和以求出概度函数。

为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合附图，作详细说明如下。

附图说明

图1示意的是多重路径传播中的各种干扰，图1a中只有延迟干扰，图1b中则有延迟干扰与多普勒频率漂移。

图2绘出的方块图示意的是以一开回路系统进行追踪的卫星导航接收器。

图3示意的是在累积功率网格中形成一解析单位(bin)的流程。

图4举例说明在一多重路径信号接收下所形成的二维累积功率网格。

图5示意中心化前后的累积功率网格，对应多普勒频率下最大信号功率。

图6示意累积功率网格的中心化流程。

图7示意测量值计算的流程以消除多重路径影响。

图8示意修正值计算流程的一种实施例，其是以信号向量S^ij＝Q^-1/2P^ij及反矩阵(S^ijTS^ij)^-1消除多重路径的影响。

图9示意修正值计算流程的一种实施例，其是利用标准正交向量V^ij削减多重路径的影响。

图10举例说明一多重路径信号的概度函数y(τ₀，τ₁)。

图11示意概度函数y(τ₀，τ₁)的投影以及直接与反射路径的延迟估计。

图12示意的是本发明针对直接路径信号延迟计算所得的估计误差。

图13示意的是依据窄相关而对直接路径信号延迟计算所得的估计误差。

图14示意的是依据选通相关而对直接路径信号延迟计算所得的估计误差。

图15比较本发明与窄相关技术的准距离估计误差。

图16示意窄相关技术下的坐标测量误差。

图17示意本发明的坐标测量误差。

【主要元件符号说明】

1  模拟部份

2  数字VLSI电路

3  程式化处理器

2.1  天线

2.2-2.9、3.1-3.5、6.1-6.2、7.1-7.6、8.1-8.7、9.1-9.6  步骤

具体实施方式

以下揭露本发明的一种实施例。

图2的方块图示意一种以开回路系统进行追踪的全球导航卫星系统接收器(亦称为导航接收器或接收器)。

在此导航接收器中，信号处理分成模拟(单元1)与数字(单元2与3)两部份。在数字信号处理中，有一部分需要高速计算，但不需处理复杂演算法，可由一数字极大尺度(digital very large scale)集成电路(下称数字VLSI电路)2来进行。另一部份需动用复杂演算法但只需低速计算，则可由一程式化处理器3来执行。

射频无线电电磁波信号由一天线(2.1)接收、转换为电子信号并放大。然后，步骤2.2进行模拟射频转换，也就是将射频转换为中频，并在时间轴上进行模拟信号的离散取样，亦对信号电平进行量化。步骤2.3用一振荡器取得一参考频率，用以进行射频转换与离散取样。然后，在步骤2.4，数字信号会以频率降转(frequencydownconversion)的方式进行数字转移。为达此目的，数字信号被转移至一预设频带频谱界定下最接近频率0的位置，并依据奈奎斯频率(Nyquist frequency)降转。此信号X(t_i)是一数字复变信号，其数字内容量有限，对存储器资源的角度来说十分经济。此信号中叠加性的混合了从各卫星接收到的信号与杂讯。接下来，此复变信号X(t_i)会被用于步骤2.5，以形成并中心化一累积功率网格Z(亦称为网格或功率网格)。此外，延迟与多普勒频率预测值P＝(τ_p，f_p)^T(P为预测值向量，τ_p为延迟的预测值，f_p则是多普勒频率的预测值)、网格中各元素相对于预测值的延迟偏移量{Δτ_n}与多普勒频率偏移量{Δf_l}都会被用于网格的形成。

然后，在步骤2.6中，基于中心化的累积功率网格Z，便可计算延迟与多普勒频率(相对于预测值)的修正值ΔΦ＝(Δτ_p，Δf_p)^T，其中，Δτ_p为延迟的修正值，Δf_p则为多普勒频率的修正值。

在步骤2.7中，前述求出的修正值会和延迟与多普勒频率的预测值加总，以求出延迟与多普勒频率的测量值：Φ＝P+ΔΦ；其中，Φ＝(τ，f)^T即是由延迟τ与多普勒频率f的测量值所形成的测量值向量。

将上述求出的测量值Φ与卫星轨迹及速度(卫星星历)综合在一起，便能在步骤2.8中决定接收器的位置、速度与时间(譬如说是根据延迟与多普勒频率测量值而以标准的方法计算出准距离与准速度)。根据求出的坐标、速度与时间，以及卫星轨迹及速度的数据，就可在步骤2.9中计算延迟与多普勒频率的预测值P。

步骤2.3与2.4可于模拟部份1中进行，步骤2.4与2.5在数字VLSI电路2中进行，步骤2.6至2.9则由程式化处理器3进行。

当1Hz频带中的载波功率与杂讯功率C/N₀＞10dB Hz(C为载波功率，N₀为杂讯功率)，本发明开回路追踪技术可确保弱信号下的安全运作。

在累积功率网格中为足标n与k形成一解析单位(bin)的流程图示于图3(n代表延迟的解析单位足标，k代表多普勒频率的解析单位足标)。

在步骤3.1中，一字码混合器(code mixer)将复变信号取样值乘上一字码的实数副本值：

$X_n^{\left(1\right)}\left(t_i\right)=X\left(t_i\right)*C_p\left(\frac{t_i-{\mathrm{\τ}}_n^{\mathrm{mix}}}{1-\frac{f_p}{F}}\right),---\left(1\right)$

其中：

X_n ⁽¹⁾(t_i)为第n个字码混合器输出的复变信号；

X(t_i)为步骤2.4输出的复变信号的取样值；

i为复变信号取样的足标；

t_i＝δ*i为第i个信号取样的时间；

δ为信号取样的时间间隔；

C_p(t)是时间t时的信号字码副本值；

τ_n ^mix＝τ_p+Δτ_n为第n个字码混合器的副本信号延迟；

τ_p为延迟预测值；

Δτ_n为第n个字码混合器的延迟偏移量；

f_p为多普勒频率的预测值；

F为载波的射频频率。

一载波混合器(步骤3.2)将字码混合器3.1输出的复变信号乘以一复变指数：

$X_{\mathrm{nl}}^{\left(2\right)}\left(t_i\right)=X_n^{\left(1\right)}\left(t_i\right)*e^{-j*2\mathrm{\π}*f_l^{\mathrm{mix}}*t_i},---\left(2\right)$

其中：

X_nl ⁽²⁾(t_i)为第n个字码混合器、第l个载波混合器输出的复变信号；

j为单位虚数；

f_l ^mix＝f_p+Δf_l为第l个载波混合器的多普勒频率值；

f_p为多普勒频率的预测值；

Δf_l为第l个载波混合器的频率偏移量。

在一段与一给定卫星的数据位元同步的同调累积时段T_coh＝N_coh*δ中，其中同调累积的时间与位元长度相匹配(对GPS与GLONASS而言，位元长度例如是20ms)，同调累积器(亦称同调积分器，步骤3.3)将复变信号加总起来：

$X_{\mathrm{nl}}^{\left(3\right)}\left(t_p\right)=\underset{i=(p-1)*\mathrm{Ncoh}}{\overset{p*N_{\mathrm{coh}}}{\mathrm{\Σ}}}{X}_{\mathrm{nl}}^{\left(2\right)}\left(t_i\right),---\left(3\right)$

其中：

X⁽³⁾ _nl(t_p)为第p个同调累积时段中同调积分器的输出；

t_p＝δ*N_coh*p为第p个同调累积时段的时间；

N_coh是在同调累积时段中复变信号的取样总数。

针对同调积分器输出的复变信号，一功率计算单元(步骤3.4)计算复变信号的模数(绝对值)的平方：

$Z_{\mathrm{nl}}^{\left(4\right)}\left(t_p\right)={\mathrm{Re}}^2\left[X_{\mathrm{nl}}^{\left(3\right)}\left(t_p\right)\right]+{\mathrm{Im}}^2\left[X_{\mathrm{nl}}^{\left(3\right)}\left(t_p\right)\right],---\left(4\right)$

其中：

Z⁽⁴⁾ _nl(t_p)为第p个同调累积时段对应的实数信号功率；

Re()、Im()为分别为一复变信号的实数部份与虚数部份。

由于同调积分器输出的复变信号X⁽³⁾ _nl(t_p)与数据位元同步，此信号的功率不会受信号中数据正负号的影响。

一非同调累积器(亦称非同步积分器，步骤3.5)在一非同调累积时段T_inc＝T_coh*N_inc中将功率累加在一起：

$Z_{\mathrm{nl}}=\underset{p=0}{\overset{N_{\mathrm{inc}}}{\mathrm{\Σ}}}{Z}_{\mathrm{nl}}^{\left(4\right)}\left(t_p\right),---\left(5\right)$

其中：

Z_nl为非同步积分器针对一非同调累积时段所输出的信号功率；

N_inc则是非同调累积时段中实数信号功率内的取样总数。

针对其他相对于预测值的延迟偏移量Δτ_n与多普勒频率偏移量Δf_l，累积功率网格中的每一元素都可用相似的步骤计算。图4即是举例说明多重路径信号接收下所形成的累积功率网格；若接收的信号是只影响延迟的反射信号(见图1a)，其所对应的累积功率网格会有所不同。

前述计算出来的累积功率(累积功率网格)中包括了信号功率与杂讯功率的加总。未知的杂讯功率成为一项干扰的因素。为排除累积杂讯，累积功率网格中不包含信号的部份会用来补偿杂讯功率，即为中心化。

对应同一多普勒频率下的一个信号峰值，图5示意的是中心化前后的累积功率网格。在累积功率网格中，延迟值小于直接信号延迟的对应元素会被用来进行中心化。

中心化的流程示于图6。中心化是针对累积功率网格中不包含信号的部份将各解析单位加以平均，即步骤6.1；再由累积功率网格的其他(或所有)元素中将平均所得的杂讯功率减去，即步骤6.2。

一开始，观察一峰值信号频率下累积功率向量随信号副本延迟的变化情形，可考虑一统计模型。

对应峰值信号频率f_l ^max，累积功率网格的中心化向量Z的模型可表述为：

其中，

Z＝(z_l1，z_f2...，z_lN)^T为中心化后累积功率向量，对应于同一峰值信号频率下的所有延迟量；

N为同一频率下功率网格元素的总数；τ＝(τ₀，τ₁...，τ_K)^T为直接信号与反射信号延迟所形成的向量，足标k＝0对应于直接信号；

K+1为路径的总数；

a＝(a₀，a₁，...，a_K)^T是由来自直接与反射路径的信号的振幅所形成的向量；

是来自直接与反射路径信号的信号的相位差所形成的向量，其中

P(τ)＝(R₀₀，R₀₁，...，R_0K，R₁₁，R₁₂，...，R_1K.，R_KK)是一个维度为(N*((K+1)*(K+2)/2))的矩阵，取决于直接路径信号的延迟τ₀与各反射路径信号的延迟τ₁，τ₂，...，τ_K；

R_km是由累积功率网格元素与第k条路径及第m条路径的信号自相关函数乘积所形成的向量，其定义为：

R_km＝(r(τ₁ ^mix-τ_k)r(τ₁ ^mix-τ_m)，r(τ₂ ^mix-τ_k)r(τ₂ ^mix-τ_m)，...，r(τ_N ^mix-τ_k)r(τ_N ^mix-τ_m))^T

τ_k，τ_m分别为第k条路径与第m条路径上的传播延迟；

r(τ)为自相关函数(ACF)；

τ_n ^mix是第n个字码混合器与累积功率网格中第n个元素的信号延迟(n＝1，...N)；

是取决于直接与反射路径上振幅a与信号相位差

的向量，其定义为：

ξ＝(ξ₁...，ξ_N)^T是中心化累积功率网格中对应杂讯波动的乱数向量，其数学期望值为0，共变异数矩阵(covariance matrix)Q＝<ξξ^T>。

由于非同调累积时段通常大于或等于1秒，累积功率网格中非同调累积的总数N_inc等于50或更多。在此情形下，杂讯分布趋近于数学期望值为0的常态分布，可由已知的共变异矩阵Q来完全定义。

在真实城市峡谷的情形下，反射路径的数目会非常多。然而，对测量误差最主要的影响通常肇因于单一最强的反射路径。如此，在两路径(一直接与一反射，K＝1)的情形下，会有5个未知的参数，其为：

直接与反射路径的振幅平方a₀ ²与a₁ ²；

直接与反射路径的相位差的余弦值

直接与反射路径的信号延迟τ＝(τ₀，τ₁)^T。

在一概度函数系统中有5个非线性参数。求解此系统会是一个复杂的计算问题。

在两个路径(直接与反射)的情形下，矩阵P的维度为N*3，向量θ的维度为3*1，N为同一频率下累积功率网格的元素数量。

在两路径情形下，矩阵P与向量θ可写为：

P(τ)＝{P₀(τ₀)，P₁(τ₁)，P₂(τ₀，τ₁)}，            (7)

其中：

P₀(τ₀)，P₁(τ₁)，P₂(τ₀，τ₁)为(N*1)向量，取决于直接与反射路径的延迟，其元素为：

P_0n(τ₀)＝r²(τ_n ^mix-τ₀)，

P_1n(τ₁)＝r²(τ_n ^mix-τ₁)，

P_2n(τ₀，τ₁)＝r(τ_n ^mix-τ₀)r(τ_n ^mix-τ₁)，

其中n＝1，...，N，

r(τ)为自相关函数(ACF)，

a₀，a₁为直接与反射路径的信号振幅，

为直接与反射路径的信号相位差。

若信号的振幅与相位差在累积功率网格的累积时段中改变，就要将式(8)中的线性参数向量替换为在同一累积时段中平均后的向量θ。在累积功率网格的累积时段中，矩阵P的变化则可被忽略。

准乱数的导航信号会在射频与数字转换步骤2.2与2.4中改变，改变的本质为确定性的，可由使用的滤波器转换函数完全定义。这些信号在计算后会被用来预先计算自相关函数r(τ)，通过信号与其累积的副本的乘算。自相关函数r(τ)取决于信号副本的种类、接收器频宽与接收器中模拟与数字滤波器的振幅-频率特性。

自相关函数可预先计算并制作为表格，以简化运算。如此，自相关函数r(τ)可被视为一已知的函数，且与接收信号与其副本间的延迟相关。自相关函数的形式亦会和GPS中准乱数序列(pseudo-random sequence，PRS)的数目与GLONASS中频率频道的数目有关。

由于非同调累积时段相较于同调累积时段来说已经足够长(如1秒对20ms)，杂讯分布会被常态化。因此，高斯分布模型会被用于本方法的合成。

将高斯杂讯分布套用至等式(6)的统计模型时，最大概度法(maximumlikelihood method，MLM)是联合测量直接与反射路径参数时趋近最佳化的方法。对应等式(6)的模型，与中心化累积功率网格Z的向量相关的概度函数(likelihoodfunction，LF)可写为：

y(τ₀，τ₁，θ)＝ln(p(θ，τ|Z)＝const-0.5*log(det(2πQ))-0.5*(Z-Pθ)^TQ^-1(Z-Pθ)-＞max  (9)

其中，

p(Z)为中心化累积功率Z向量的概率分布密度函数，其可定义为：

p(Z)＝{det(2πQ)}^-1/2exp{-0.5*(Z-P(τ)θ)^TQ^-1(Z-P(τ)θ)}。

此概度函数为非线性函数，取决于5个非线性参数。等式(9)表示需要针对其中的所有参数找出概度函数的极大值。在5个非线性概度函数形成的系统中求解会是一个复杂的计算问题。

为了针对延迟τ₀及τ₁找出等式(9)中概度函数的极大值，首先可利用最大概度θ_ml的估计值，其与参数θ的线性部份有关，可由下列外显型式的等式求得：

θ_ml＝[P(τ)^TQ^-1P(τ)]^-1P(τ)^TQ^-1Z            (10)

然后，求得的估计值θ_ml可被代入至等式(9)，以决定概度相对于参数τ₀、τ₁的极大值y(τ₀，τ₁)：

y(τ₀，τ₁)＝-0.5*(Z^TQ^-1Z-Z^TQ^-1P(τ)(P(τ)^TQ^-1P(τ))^-1P(τ)^TQ^-1Z)-＞max(11)

再引用下列的额外设定：

S＝Q^-1/2P(τ)为杂讯白化(noise whitening)之后的信号向量，取决于直接与反射路径的延迟；

X＝Q^-1/2Z为杂讯白化后累积功率的中心化向量。

将S与X代入至等式(11)后，概度函数可被重新表示为延迟τ₀，、τ₁的函数：

y(τ₀，τ₁)＝-0.5*X^T[E-S(S^TS)^-1S^T]X＝-0.5*X^TX+0.5*X^TS

(S^TS)^-1S^TX-＞max      (12)

等式(12)的第一个被加数与参数τ₀、τ₁无关，故只需计算第二个被加数的最大值。因此，通过计算等式(12)的第二个被加数可为延迟参数Δτ₀、Δτ₁定义出一个修正测量的方法。

在开回路追踪中计算测量值以消除多重路径影响的流程示于图7。在步骤7.1中，二维中心化累积功率网格Z被用于二维内插、最小平方法或其他已知方法以计算单一路径信号的延迟与多普勒频率修正值ΔΦ＝(Δτ，Δf)^T。然后，在步骤7.2中，则根据网格形状与一参考形状(对应于一单一路径信号)间的差异以检测多重路径。若差异未超过一多重路径检测的临界值，则在步骤7.6中根据修正值计算延迟与多普勒频率的估计值(测量值)Φ＝(τ，f)^T＝P+ΔΦ。

若检测到多重路径，则在步骤7.3中改进修正值ΔΦ(详述于后)，并在步骤7.6中根据修正值计算延迟与多普勒频率的测量值。针对反射路径下对单一路径信号进行修正的方法进行误差分析，便可为直接路径信号定义出可能的延迟偏移量Δτ_i＝τ-τ₀的组合。在反射信号对直接路径延迟的测量精确度有重大影响的区带(band)中，可选出反射路径信号对应的延迟Δτ_j。

在改善步骤7.1的修正值时，可求出多个权重，譬如说是使用步骤7.5的权重库。为此，直接与反射信号的延迟区带可被分为多个时段；为各个时段，步骤7.4会对计算概度函数所需的权重先行计算(可以不必动用即时系统)，并将这些先行计算好的系数写入至权重库中。

根据等式(12)测量各相关集合的流程示于图8。

步骤8.1进行的是杂讯白化(将中心化累积功率网格Z的向量乘以矩阵Q^-1/2)：

$\left\{X_n\right\}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{Q}_{\mathrm{nk}}^{-1/2}{Z}_k---\left(13\right)$

举例来说，可根据丘列斯基(Cholesky)法分解出两个三角矩阵并进行三角矩阵逆运算以计算矩阵Q^-1/2[Trefethen L.N.，Bau D.Numerical Linear Algebra，SIAM，1997，p.172]。要注意的是，杂讯白化只能降低测量杂讯，对多重路径导致的偏移误差没有效果。经由模型显示，即使不经白化(使Q＝E，E为单位矩阵)，在多重路径下对直接信号的延迟测量仍会保持无偏移的状态(图12)。因此，杂讯白化在本发明方法中并不是强制一定要进行的。

在步骤8.2中，会针对每一组直接与反射信号延迟τⁱ ₀，、τ^j ₁计算出三组权重，分别对应向量S₁ ^ij，S₂ ^ij，S₃ ^ij；在步骤8.3中，则利用这些权重对向量X的各元素计算加权总和：

$h_q^{\mathrm{ij}}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_n{S}_{\mathrm{qn}}^{\mathrm{ij}}---\left(14\right),$

其中，

q＝1，2，3代表矩阵S^ij的直列向量(column vector)足标，对应于一组直接与反射信号延迟；

i，j为概度函数格点上分别与直接及反射信号延迟相对应的足标。

步骤8.3之后，在步骤8.5中会将加权总和h₁ ^ij、h₂ ^ij与h₃ ^ij配对相乘，并配合步骤8.4中求得的3*3的反矩阵R^ij＝(S^ijTS^ij)^-1进行加权相加以计算概度函数：

${\mathrm{\&epsiv;}}^{\mathrm{ij}}=\underset{g=1}{\overset{3}{\mathrm{\&Sigma;}}}\underset{q=1}{\overset{3}{\mathrm{\&Sigma;}}}{h}_g^{\mathrm{ij}}{h}_q^{\mathrm{ij}}{R}_{\mathrm{gq}}^{\mathrm{ij}}---\left(15\right)$

此运算代表一二次形式的计算。

计算结果为LF{ε^ij}，于步骤8.6中储存于存储器。

此程序会对每一对直接与反射路径的延迟重复进行(每一对延迟关联于对应的权重向量S^ij与反矩阵R^ij的值)。之后，在步骤8.7，就可根据所有直接与反射信号的延迟参数τⁱ ₀，、τ^j ₁所对应的全部概度函数值{ε^ij}内插决定全域极大值。

根据窄相关下针对准距离误差所进行的分析(图13)，可为直接与反射路径决定一延迟搜寻区带。一开始可先针对累积功率网格的极大值位置决定其所对应的延迟τ_max。然后，便可在τ_max相关的区带中选出直接路径对应的延迟τⁱ ₀：

τ_max-15m/C＜τⁱ ₀＜τ_max+15m/C，

其中C为光速。

在反射路径会影响直接路径延迟估计准确度的区带中，可选出反射路径延迟τ^j ₁，也就是：

τ_max+15m/C＜τ^j ₁＜τ_max+300m/C。

由模型分析可知，搜寻直接路径延迟的步级应在1到2m/C，搜寻反射路径延迟的步级则需要10到20m/C。

对应两(一直接与一反射)路径信号的概度函数举例示意于图10与11。

二维概度函数最大值的抛物线近似值被用于求出延迟的估计值

分别对应直接与反射路径。

图12至14的比较性模型分析将相同情形下运用本发明、窄相关与选通相关(strobe correlator)所分别获得的精确度加以比较(接收信号频率的频宽ΔF＝6MHz，反射路径信号衰减等于6dB)。

如图12所示，针对直接路径信号的延迟测量，本发明的偏移误差非常小(小于0.12m/C)，相对优于窄相关技术(图13，误差偏移超过15m/C)与选通相关技术(图14，偏移超过7m/C)。

在以上描述(图8)中，改进修正值的方法会用到向量S＝Q^-1/2P、X＝Q^-1/2Z与反矩阵(S^ijTS^ij)^-1。此方案的缺点在于要和向量S一起储存矩阵R＝(S^TS)^-1，以计算步骤8.5中的二次形式。以下将近一步考虑改进修正值的另一种方法，其系运用和直接与反射路径延迟相依的标准正交(orthonormal)向量。

矩阵(S^TS)为3*3的对称非负定义(non-negatively defined)矩阵，此矩阵可表为以下形式

(S^TS)＝UΛU^T，        (16)

其中，

U为一个维度为3*3的方阵，由标准正交的固有向量(eigenvector)形成；

Λ＝diag(λ₁，λ₂，λ₃)则是一个具有非负固有值(eigenvalue)的对角矩阵。

据此，反矩阵(S^TS)^-1就可表示为

(S^TS)^-1＝UΛ^-1U^T＝(UΛ^-1/2)(UΛ^-1/2)^T    (17)

而矩阵S(S^TS)^-1S^T则可写为：

S(S^TS)^-1S^T＝(SUΛ^-1/2)(SUΛ^-1/2)^T＝VV^T    (18)

其中V＝(V₁，V₂，V₃)＝SUΛ^-1/2是一个N*3矩阵，其直列向量V₁、V₂、V₃为矩阵S(S^TS)^-1S^T的三个标准正交固有向量。此矩阵仅取决于直接与反射路径的延迟(τ₀，τ₁)。将等式(18)代入等式(12)后，概度函数的计算可写成：

y(τ₀，τ₁)＝const+0.5*X^TV(τ₀，τ₁)V(τ₀，τ₁)T^X-＞max(19)

利用标准正交向量可简化本发明方法的数字实施。在此情形下，不需计算反矩阵，二次形式为对角线类型，可以使用整数运算。

另一种相似的计算简化方法是根据丘列斯基法将矩阵S^TS分解为两个三角矩阵(如前述参考文献[Trefethen L.N.，Bau D.Numerical Linear Algebra，SIAM，1997，p.172])：

(S^TS)＝Δ^T Δ，                (20)

其中：

Δ为一下三角(lower triangular)矩阵，

Δ^T为一上三角(upper triangular)矩阵。

在此情形下，反矩阵等于

(S^TS)^-1＝(Δ^TΔ)^-1＝Δ^-1(Δ^-1)^T。        (21)

转换为三角矩阵后，运算就更为简化。在此情形下，标准正交向量V可根据下列等式预先计算：

V＝(V₁，V₂，V₃)＝SΔ^-1。            (22)

其他更进一步的处理则与等式(19)固有向量相似。

标准正交向量V亦可根据格兰-史密特(Gram-Schmidt)标准正交化程序而由向量S求得([Trefethen L.N.，Bau D.Numerical Linear Algebra，SIAM，1997，p.56])。

图9示意的是利用信号的标准正交向量消除多重路径影响。在此实施例与前述讨论实施例中，可包括(但非必要)杂讯白化(步骤9.1)，并包括将概度函数的值储存于存储器(步骤9.5)以及内插求得全域极大值(步骤9.6)。

与先前实施例不同的是，当在步骤9.3中为每一向量计算加权总和时，是使用标准正交向量V^ij＝(V^ij ₁，V^ij ₂，V^ij ₃)来进行：

$h_g^{\mathrm{ij}}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_n{V}_{\mathrm{gn}}^{\mathrm{ij}}---\left(23\right)$

其中g＝1，2，3为固有向量的足标。

此外，在步骤9.4中，加权总和h^ij _g的平方和计算为：

${\mathrm{\&epsiv;}}^{\mathrm{ij}}=\underset{g=1}{\overset{3}{\mathrm{\&Sigma;}}}{\left(h_g^{\mathrm{ij}}\right)}^2---\left(24\right)$

此为概度函数之值，对应于直接与反射路径的信号延迟τⁱ ₀，τ^j ₁。

在同频率反射路径存在的情形下，减少直接路径延迟测量误差的方法就如前面讨论过的。

对一个移动的消费者来说，其所观察到的多重路径信号会因多普勒效应而有频域中的位移(offset)。在此情形下，直接路径信号的多普勒位移不会偏移，反射路径的多普勒位移则取决于反射面的位置(见图1b)。针对此种情形，本发明方法可用来降低多普勒频率测量的误差。为达此目的，可在累积功率网格中针对各频率与最大功率对应的延迟τ_n ^max所对应的元素形成一向量并加以运用。在依据反射路径信号的频率差异进行多重路径检测时，可比较信号累积功率的频率分布。

在此情形下，矩阵P会与等式(7)不同，其可表示为下列形式：

P＝P(f)＝(P₀，P₁，P₂)，            (25)

其中：

f＝(f₀，f₁)^T是未知的多普勒位移f₀与厶形成的向量，分别对应直接路径与反射路径；

P₀，P₁，P₂分别为L*1维度的信号向量，取决于直接与反射路径的多普勒频率，其元素为：

P_0l＝sinc²[π(f_l ^mix-f₀)T_coh]，

P_1l＝sinc[π(f_l ^mix-f₀)T_coh]sinc[π(f_l ^mix-f₁)T_coh]，

P_2l＝sinc²[π(f_l ^mix-f₁)T_coh]，

l＝1，...，L代表累积功率网格中对应频率的足标；

f_l ^mix为第l个载波混合器的频率(亦为累积功率网格的第l个元素)，

T_coh为同调累积时段。

函数sinc(x)的定义为：

sinc(x)＝sin(x)/[N_coh*sin(x/N_coh)]。

中心化累积功率网格的解析单位(bin)可用来形成向量Z＝(Z_1n，Z_2n，...，Z_Ln)^T，其格点对应最大功率延迟τ_n ^max下的各频率。明确的修正值Δf₀，Δf₁可由概度函数的最大值求得，而概度函数可利用信号向量得到，如等式(26)：

y(f₀，f₁)＝-0.5*X^T[E-S(S^TS)^-1S^T]X＝-0,5*X^TX+0，5*X^TS(S^TS)^-1S^TX-＞max，

(26)

其中：

S＝Q^-1/2P(f₀，f₁)为(杂讯白化后)信号向量形成的矩阵，取决于直接与反射路径的多普勒位移；

X＝Q^-1/2Z为(杂讯白化后)的中心化累积功率向量，与峰值信号延迟及该延迟下的所有频率相对应。

等式(26)亦可用信号的标准正交向量来简化：

y(f₀，f₁)＝const+0,5*X^TV(f₀，f₁)V(f₀，f₁)^TX-＞max，    (27)

其中：

V(f₀，f₁)为标准正交向量形成的矩阵，取决于多普勒频率。

此方法的实现可由图8、9中与延迟相关的实施例类推而得。

为确认本发明的效果，图15、16与17呈现比较性测试的结果。

图15显示的是多重路径信号下准距离的误差，分别对应本发明(曲线1，实线)与窄相关延迟估计法(曲线2，虚线)。在此测试中，累积功率网格的解析单位等于0.125μs；多重路径信号则由一多重频道模拟器(multichannel simulator)求出，使用的反射参数为：衰减-6dB与延迟0.5μs；直接与反射路径的多普勒位移差异等于0.012Hz(对应静态接收器)，接收器频宽为6MHz。

运用本发明方法后，准距离的误差降低了4到5倍。

图16示意的是以窄相关技术求得的点分布图谱(fuzzball)，图17则是以本发明方法求得的点分布图谱。由两图比较可知，在决定接受器的坐标时，本发明可将水平面的坐标误差ΔX，ΔY降低4到5倍。

由此可知，本发明的延迟与多普勒频率测量方法可运用功率网格数据而在弱信号情形下抑制多重路径误差。

复变信号的同调累积是在同调累积时段中以位元同步方式进行；在GPS与GLONASS系统中，同调累积时段符合的数据位元长度等于20ms。这可避免传输信号中位元长度未知的导航讯息干扰，并使可接受的讯杂比低于数据解码临界值(data decoding threshold)。

若数据位元已知或可被解码，对应这些位元的复变信号载波相位反转就可被补偿。在此情况下，同调累积时段可根据参考频率振荡器的特性来决定，通常等于十分之几秒。针对此种情况，根据累积功率网格抑制多重路径误差亦可在低讯杂比下进行。

根据累积功率网格数据检测多重路径后，本发明可有效抑制多重路径在延迟(准距离)与多普勒频率(准速度)上导致的测量误差。

本发明让抑制多重路径误差所需的运算可每秒执行一次，并可实现于程式化的处理器。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明，本领域技术人员在不脱离本发明的精神和范围内，当可作些许的更动与润饰，因此本发明的保护范围当以权利要求所界定的为准。

Claims (10)

1.一种于卫星导航接收器中抑制多重路径误差的方法，包含搜寻多个导航卫星的信号、由各该卫星接收并处理复变信号、以一追随开回路追踪信号并基于直接信号与反射信号的延迟与多普勒频率测量值决定该卫星导航接收器的一坐标、一速度与一准确时间的步骤，该方法包含：

根据多个导航参数计算该延迟与该多普勒频率的预测值以根据该延迟与该多普勒频率定向一二维的累积功率网格，该累积功率网格是依据该接收与处理的该复变信号所形成，而该些导航参数包括该坐标、该速度、该准确时间，以及根据测量到的时间与接收到的卫星星历求得的卫星运动信息；

其中，在形成该累积功率网格的各元素时，在一同调累积时段中进行复变信号加总，该同调累积时段同步于一给定卫星的数据位元，该同调累积时段长短符合一位元长度；以及

基于该累积功率网格的处理，包括该累积功率网格的中心化，为该延迟与该多普勒频率的该预测值计算修正值，并计算该预测值与该修正值的加总以作为该延迟与该多普勒频率的测量值，再对该测量值外插一预设时间并计算准距离与准速度的测量值，并根据对所有卫星所求得的该测量值决定该卫星导航接收器的该坐标与该速度；

其中，计算该延迟与该多普勒频率的该修正值时，包含：

利用该中心化的该累积功率网格计算单一路径的修正值；

决定多重路径是否发生；以及

若为肯定的，则利用多个权重改进该修正值；

其中，当利用该些权重改进该修正值时，基于一中心化的累积功率向量形成一概度函数，并根据该概度函数的全域极大值改进该修正值。

2.如权利要求1所述的方法，其特征在于，基于该中心化的该累积功率向量形成该概度函数的步骤包含下列子步骤：

利用多个权重向量而为每一该累积功率向量计算多个加权总和，该些权重向量对应该直接信号与该反射信号的该延迟的该预测值；

针对该些权重向量的两两纯量内积所形成的矩阵求出反矩阵；以及

根据前述两子步骤的结果计算该概度函数，使该概度函数呈一二次形式。

3.如权利要求1所述的方法，其特征在于，基于该中心化的该累积功率向量形成该概度函数的步骤包含下列子步骤：

针对每一累积功率向量，利用多个取决于该直接信号与该反射信号的该延迟的标准正交向量计算多个加权总和；以及

计算该多个加权总和的平方总和以求出该概度函数。

4.如权利要求1所述的方法，其特征在于，基于该中心化的该累积功率向量形成该概度函数时的步骤包含下列子步骤：

利用多个权重向量而为每一该累积功率向量计算多个加权总和，该些权重向量对应该直接信号与该反射信号的该多普勒频率的该预测值；

针对该些权重向量的两两纯量内积所形成的矩阵求出反矩阵；以及

根据前述两子步骤的结果计算该概度函数，使该概度函数呈一二次形式。

5.如权利要求1所述的方法，其特征在于，基于该中心化的该累积功率向量形成该概度函数的步骤包含下列子步骤：

针对每一该累积功率向量，利用多个取决于该直接信号与该反射信号的该多普勒频率的标准正交向量计算多个加权总和；以及

计算该多个加权总和的平方总和以求出该概度函数。

6.一种于卫星导航接收器中抑制多重路径误差的装置，包含：

搜寻多导航卫星的信号的装置；

由每一该卫星接收并处理复变信号的装置；

以一追随开回路追踪信号并基于直接信号与反射信号的延迟与多普勒频率测量值决定该卫星导航接收器的一坐标、一速度与一准确时间的装置；

根据多个导航参数计算该延迟与该多普勒频率的预测值的装置，该装置根据该延迟与该多普勒频率定向一二维的累积功率网格，该累积功率网格是依据该接收与处理的该复变信号所形成，而该些导航参数包括该坐标、该速度、该准确时间，以及根据测量到的时间与接收到的卫星星历求得的卫星运动信息；

其中，在形成该累积功率网格的各元素时，在一同调累积时段中进行复变信号加总，该同调累积时段同步于一给定卫星的数据位元，该同调累积时段长短符合一位元长度；以及

基于该累积功率网格的处理，对该累积功率网格中心化，为该延迟与该多普勒频率的该预测值计算修正值，并计算该预测值与该修正值的加总以作为该延迟与该多普勒频率的测量值，再对该测量值外插一预设时间并计算准距离与准速度的测量值，并根据对所有卫星所求得的该测量值决定该卫星导航接收器的该坐标与该速度的装置；

其中，计算该延迟与该多普勒频率的该修正值的该装置包含：

利用该中心化的该累积功率网格计算单一路径的修正值的装置；

决定多重路径是否发生的装置；以及

若为肯定的，则利用多个权重改进该修正值的装置，其中，当利用该权重改进该修正值时，系基于一中心化的累积功率向量形成一概度函数，并根据该概度函数的全域极大值改进该修正值。

7.如权利要求6所述的装置，其特征在于，利用该些权重改进该修正值的装置包含下列装置以基于该中心化的该累积功率向量形成该概度函数：

利用多个权重向量而为每一该累积功率向量计算多个加权总和的装置，该些权重向量对应该直接信号与该反射信号的该延迟的该预测值；

针对该些权重向量的两两纯量内积所形成的矩阵求出反矩阵的装置；以及

利用该些加权总和与该反矩阵计算该概度函数的装置，使该概度函数呈一二次形式。

8.如权利要求6所述的装置，其特征在于，利用该些权重改进该修正值的装置包含下列装置以基于该中心化的该累积功率向量形成该概度函数：

针对每一累积功率向量，利用多个取决于该直接信号与该反射信号的该延迟的标准正交向量计算多个加权总和的装置；以及

计算该多个加权总和的平方总和以求出该概度函数的装置。

9.如权利要求6所述的装置，其特征在于，利用该些权重改进该修正值的装置包含下列装置以基于该中心化的该累积功率向量形成该概度函数：

利用多个权重向量而为每一该累积功率向量计算多个加权总和的装置，该些权重向量对应该直接信号与该反射信号的该多普勒频率的该预测值；

针对该些权重向量的两两纯量内积所形成的矩阵求出反矩阵的装置；以及

利用该些加权总和与该反矩阵计算该概度函数的装置，使该概度函数呈一二次形式。

10.如权利要求6所述的装置，其特征在于，利用该些权重改进该修正值的装置包含下列装置以基于该中心化的该累积功率向量形成该概度函数：

针对每一累积功率向量，利用多个取决于该直接信号与该反射信号的该多普勒频率的标准正交向量计算多个加权总和的装置；以及

计算该多个加权总和的平方总和以求出该概度函数的装置。

Images