CN104749587A - 接收机伪距故障监测方法和接收机 - Google Patents

Abstract

本申请提出了一种接收机伪距故障监测方法和接收机。该接收机伪距故障监测方法包括：获得接收机与多颗卫星之间的伪距；根据所获得的伪距，计算接收机的位置和钟差；根据所获得的伪距以及所计算的接收机的位置和钟差，获得卫星几何结构投影矩阵S和残差向量b；设定伪距误差容忍度δ，根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少；以及根据所计算的伪距误差向量ε，监测接收机伪距故障。

Description

接收机伪距故障监测方法和接收机

技术领域

本申请涉及接收机伪距故障监测方法和接收机。

背景技术

全球导航卫星系统（GNSS）能够在全球范围内全天候地为具有接收设备的用户提供实时、准确的位置、速度和时间信息。目前，导航产业正在迅速发展，并逐渐渗透到各个行业当中。在军事需求和经济利益的推动下，许多国家和地区都在积极地升级或建设自己的卫星导航系统。随着多个导航系统的建成并开放运行，未来的导航终端将进入多模多频的时代。

图1所示显示了现有的卫星导航系统，车载或手持接收机接收来自卫星的广播信号（电磁波），当卫星数量大于等于4颗时，构成最基本的定位系统。定位系统可以是如图1所示的包括GPS，GLONASS，BeiDou等导航系统的多模多频定位系统，也可以是一种导航系统，例如GPS导航系统或者BeiDou导航系统。

每颗卫星上都具有高精度的时钟，地面站计算出各个时钟相对于基准时间系统的修正值，并通过上行链路注入到卫星的星历存储单元。卫星在广播星历中播发时钟修正信息，用户接收后可以对卫星星钟进行修正，从而可以近似认为各个卫星的时钟相互同步。用户从接收的卫星导航信号中提取信号发射时刻，和用户的当前时刻做差后，得到信号的传输时间。将此时间乘以光速后，得到卫星与用户之间的伪距。由于用户的时间信息并不精确，以及提取信号发射时刻时会受到噪声等因素的影响，因而伪距并非卫星到用户之间的准确距离。各个卫星的位置不同，因此各个卫星到用户之间的伪距并不相同，伪距中包含的误差大小也并不相同。

但是，对于每一次定位解算，每一颗卫星的伪距误差是确定的，误差的大小反映了该颗卫星的伪距观测量是否存在异常。当某颗卫星的伪距观测量不存在异常时，它的伪距误差应服从零均值的高斯分布，其期望值为零。当某颗卫星的伪距观测值发生故障时，它的伪距误差将不再服从零均值的高斯分布，其期望值为故障偏差。

可能导致伪距观测故障的原因有很多。比如：卫星星钟发生故障，电离层、对流层出现异常，用户接收位置处存在严重多径等。

接收机自主完好性监测（RAIM）方法是接收机内部的、用于检测并排除异常观测量并给出定位结果及其可靠性的算法。RAIM方法的性能决定了接收机能否在高可靠性要求的场景下得以应用，从一定程度上制约着卫星定位应用的延伸与拓展。

传统RAIM方法仅适用于单星故障的检测与排除。随着多个导航系统的建设，多模多频接收机逐渐引入，RAIM方法亟需进行改进以满足多星故障场景下的应用需求。目前的多星故障RAIM方法主要分为两种。一种是在单历元时刻通过分组方法进行多星故障检测；另一种是联合利用多个历元时刻的信息进行多星故障检测。

Georg Schroth等人在2008年提出了RANCO算法。该方法需要按照一定规则将可见星分成若干个子集，然后对每个子集进行解算，查找存在异常的伪距观测量。但是，当可见卫星较多时，子集的划分比较困难，遍历所有子集进行计算也会使算法过于复杂。

Ilaria Martini等人在2006年提出了联合利用多个历元时刻的信息进行多星故障检测的方法。其基本思想是，在某个历元时刻，卫星的几何结构构成一个投影矩阵，误差向量经过投影变换得到残差向量。在下一个历元时刻，卫星的几何结构发生变化，得到新的投影矩阵和残差向量。因此，在这两个历元时刻，得到了误差向量在两个投影坐标系下的不同表示，通过重构可以得到误差向量。这种方法能恢复出所有卫星的测距误差。但是，进行重构的条件是投影矩阵之间要线性无关，误差向量在历元时刻之间保持不变。对于地面低速运动的用户，卫星的几何结构变化缓慢，因而相邻历元时刻投影矩阵之间有较强的相关性。这种相关性在进行测距误差重构时会引入较强的病态性，导致计算误差很大。在测距误差较小时，重构结果不准确。而且，误差向量在历元时刻之间保持不变在实际中也较难保证。

发明内容

本申请的目的是提供一种至少能够降低算法复杂度的接收机伪距故障监测方法和接收机。

根据本申请的一个方面，公开了一种接收机伪距故障监测方法，包括：获得接收机与多颗卫星之间的伪距；根据所获得的伪距，计算接收机的位置和钟差；根据所获得的伪距以及所计算的接收机的位置和钟差，获得卫星几何结构投影矩阵S和残差向量b；设定伪距误差容忍度δ，根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少；以及根据所计算的伪距误差向量ε，监测接收机伪距故障。

根据本申请的另一个方面，公开了一种接收机，能够监测伪距故障，所述接收机包括：接收模块，获得接收机与多颗卫星之间的伪距，并根据所获得的伪距，计算接收机的位置和钟差；处理模块，根据所获得的伪距以及所计算的接收机的位置和钟差，获得卫星几何结构投影矩阵S和残差向量b；伪距误差向量计算模块，设定伪距误差容忍度δ，根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少；以及监测模块，根据所计算的伪距误差向量ε，监测接收机伪距故障。

附图说明

图1所示显示了现有的卫星导航系统。

图2显示了根据本申请的一种实施方式的接收机伪距故障监测方法流程图。

图3显示了根据本申请的一种实施方式的接收机的示意方框图。

图4显示了根据本申请的一种实施方式的接收机的示意方框图

具体实施方式

下面参照附图对本申请公开的接收机伪距故障监测方法和接收机进行详细说明。为简明起见，本申请各实施例的说明中，相同或类似的装置使用了相同或相似的附图标记。

图2显示了根据本申请的一种实施方式的接收机伪距故障监测方法流程图。如图所示，在步骤110中，获得接收机与多颗卫星之间的伪距。例如，对于n颗卫星组成的定位系统，可以获得的伪距为 $\mathrm{\ρ}=\left[\begin{array}{c}{\mathrm{\ρ}}_1\\ {\mathrm{\ρ}}_2\\ .\\ .\\ .\\ {\mathrm{\ρ}}_n\end{array}\right].$

在步骤120中，根据所获得的伪距，计算接收机的位置(x_u,y_u,z_u)和钟差δ_u。

在步骤130中，根据所获得的伪距ρ以及所计算的接收机的位置(x_u,y_u,z_u)和钟差δ_u，获得卫星几何结构投影矩阵S和残差向量b。

在步骤140中，设定伪距误差容忍度δ；根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少。

在步骤150中，根据所计算的伪距误差向量ε，监测接收机伪距故障。

在卫星定位接收机的定位解算阶段，可以利用卫星位置和观测伪距构造如下方程组，求解用户的三维位置和时间偏差，如式（1）所示。

$\left\{\begin{array}{c}\sqrt{{(x_{s,1}-x_u)}^2+{(y_{s,1}-y_u)}^2+{(z_{s,1}-z_u)}^2}+{\mathrm{\δ}}_u+{\mathrm{\ϵ}}_1={\mathrm{\ρ}}_1\\ \sqrt{{(x_{s,2}-x_u)}^2+{(y_{s,2}-y_u)}^2+{(z_{s,2}-z_u)}^2}+{\mathrm{\δ}}_u+{\mathrm{\ϵ}}_2={\mathrm{\ρ}}_2\\ \sqrt{{(x_{s,n}-x_u)}^2+{(y_{s,n}-y_u)}^2+{(z_{s,n}-z_u)}^2}+{\mathrm{\δ}}_u+{\mathrm{\ϵ}}_n={\mathrm{\ρ}}_n\end{array}\right.---\left(1\right)$

其中：(x_u,y_u,z_u)为接收机的三维位置，δ_u为接收机钟差，(x_s,i,y_s,i,z_s,i)表示第i颗卫星的位置，ρ_i为用户到第i颗卫星的伪距，ε_i为伪距测量误差。i=1,2,…,n，n为可见卫星总数。

使用牛顿迭代法求解上述非线性方程组，首先需要对方程左边进行线性化。在近似解处进行一阶泰勒展开。第i个方程可写为：

$d{\mathrm{\ρ}}_i=\frac{-(x_{s,i}-{\hat{x}}_u)}{{\hat{r}}_i}\mathrm{dx}+\frac{-(y_{s,i}-{\hat{y}}_u)}{{\hat{r}}_i}\mathrm{dy}+\frac{-(z_{s,i}-{\hat{z}}_u)}{{\hat{r}}_i}\mathrm{dz}+d{\mathrm{\δ}}_u---\left(2\right)$

其中：为伪距增量， ${\hat{r}}_i=\sqrt{{(x_{s,i}-{\hat{x}}_u)}^2+{(y_{s,i}-{\hat{y}}_u)}^2+{(z_{s,i}-{\hat{z}}_u)}^2}$ 为第i颗卫星到用户近似位置的标量距离。 $\mathrm{dx}=x_u-{\hat{x}}_u,$ $\mathrm{dy}=y_u-{\hat{y}}_u,$ $\mathrm{dz}=z_u-{\hat{z}}_u,$ $d{\mathrm{\δ}}_u={\mathrm{\δ}}_u-{\widehat{\mathrm{\δ}}}_u.$

线性化后的矩阵方程式为：

Gx+ε=b    （3）

其中：

$G=\left[\begin{array}{cccc}\frac{-(x_{s,1}-{\hat{x}}_u)}{{\hat{r}}_1}& \frac{-(y_{s,1}-{\hat{y}}_u)}{{\hat{r}}_1}& \frac{-(z_{s,1}-{\hat{z}}_u)}{{\hat{r}}_1}& 1\\ \frac{-(x_{s,2}-{\hat{x}}_u)}{{\hat{r}}_2}& \frac{-(y_{s,2}-{\hat{y}}_u)}{{\hat{r}}_2}& \frac{-(z_{s,2}-{\hat{z}}_u)}{{\hat{r}}_2}& 1\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \frac{-(x_{s,n}-{\hat{x}}_u)}{{\hat{r}}_n}& \frac{-(y_{s,n}-{\hat{y}}_n)}{{\hat{r}}_n}& \frac{-(z_{s,n}-{\hat{z}}_u)}{{\hat{r}}_n}& 1\end{array}\right]---\left(4\right)$

$x=\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{dy}\\ \mathrm{dz}\\ d{\mathrm{\δ}}_u\end{array}\right],\mathrm{\ϵ}=\left[\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_n\end{array}\right],b=\left[\begin{array}{c}d{\mathrm{\ρ}}_1\\ d{\mathrm{\ρ}}_2\\ .\\ .\\ .\\ d{\mathrm{\ρ}}_n\end{array}\right]---\left(5\right)$

线性化矩阵G由当前时刻各颗可见卫星的相对位置决定，称之为几何矩阵。可以使用最小二乘法进行求解x为：

$\tilde{x}={\left(G^{T}G\right)}^{-1}{G}^{T}b---\left(6\right)$

在卫星信号正常接收情况下，伪距测量误差项ε应服从零均值的高斯分布。考虑到各个卫星到达用户的距离和空间传输环境不同，可以通过加权系数矩阵W为接收机所获得的与多颗卫星之间的伪距进行加权。经过加权后，伪距测量误差向量ε中各项方差不尽相同。不失一般性，我们考虑加权最小二乘解。定义加权系数矩阵W为：

$W=\mathrm{Cov}{(\mathrm{\ϵ},\mathrm{\ϵ})}^{-1}\stackrel{\mathrm{\Δ}}{=}{D}^{T}D---\left(7\right)$

则（3）式的加权最小二乘解为：

$\tilde{x}={\left(G^T\mathrm{WG}\right)}^{-1}{G}^T\mathrm{Wb}---\left(8\right)$

利用（8）式，根据定义计算伪距残差向量

$\tilde{b}:=b-G\tilde{x}=(I-G{\left(G^T\mathrm{WG}\right)}^{-1}{G}^{T}W)b---\left(9\right)$

将（3）式代入到（9）式中，得到：

$\tilde{b}=(I-G{\left(G^T\mathrm{WG}\right)}^{-1}{G}^{T}W)(\mathrm{Gx}+\mathrm{\ϵ})=(I-G{\left(G^T\mathrm{WG}\right)}^{-1}{G}^{T}W)\mathrm{\ϵ}$ $---\left(10\right)$

其中，

$S\stackrel{\mathrm{\Δ}}{=}I-G{\left(G^T\mathrm{WG}\right)}^{-1}{G}^{T}W---\left(11\right)$

则有：

$\mathrm{S\ϵ}=\tilde{b}---\left(12\right)$

式给出了矩阵S，误差向量ε和残差向量的关系。称矩阵S为卫星几何结构投影矩阵，通过投影变换，将误差向量ε映射为残差向量b。

这样，就获得了卫星几何结构投影矩阵S和残差向量b。从而，如上述步骤140所述，可以根据伪距误差方程Sε=b，计算伪距误差向量ε，使得伪距误差向量ε中的各伪距误差高于伪距误差容忍度δ的个数最少。

根据本申请的一种实施方式，伪距误差向量中应有尽可能多的元素服从零均值高斯分布。从而可以在该条件下，计算误差向量ε。也就是说，伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数尽可能少。其中，伪距误差容忍度是用于判断向量ε中的测距误差是否算作故障的门限，若向量ε中的某个元素的绝对值大于伪距误差容忍度δ，则认为该测距值存在故障。根据一种实施例，容忍度δ可以根据伪距测量的标准差σ来选取。根据一种实施例，δ可以取值为2σ到8σ之间的数值。例如，δ可以取值为3σ或者6σ。

为了计算误差向量ε，首先给出两个定义。

定义1：N(δ,x)，δ＞0表示向量x中的元素的绝对值大于δ的个数。

定义2：表示向量x在容忍度为δ下的l₀范数。

根据上述定义，为了使得伪距误差向量ε中的各伪距误差高于伪距误差容忍度δ的个数最少，可以求解使得

根据本申请的一种实施方式，为了求解可以首先构造第一目标函数：

$f\left(\mathrm{\ϵ}\right)={\left|\right|\mathrm{\ϵ}\left|\right|}_{l_0,\mathrm{\δ}}+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2---\left(13\right)$

由于l₀范数的存在，目标函数1并不是凸函数。求解使f(ε)达到最小值的向量ε较为困难。为此，对进行放大，构造新的目标函数，使其为凸函数。

将近似为‖ε‖₁，构造第二目标函数

$g\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\δ}}{\left|\right|\mathrm{\ϵ}\left|\right|}_1+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2---\left(14\right)$

这样，可以利用稀疏求解方法通过求解第二目标函数g(ε)的极小值对应的向量ε₀，该向量ε₀即为计算所得到的伪距误差向量ε。其中，稀疏求解算法例如可以包括LARS（Least Angle Regression）算法、BP（BasisPursuit）算法、IRWLS（Iterative Reweighted Least Squares）算法、IST（Iterative Soft Thresholding）算法、MP（Matching Pursuit）算法、OMP（Orthogonal Matching Pursuit）算法、stOMP（Stagewise OrthogonalMatching Pursuit）算法、Stepwise算法、PFP（Polytope Faces Pursuit）算法、Lasso（Least Absolute Shrinkage and Selection Operator）算法等。

实际上，接收机获得与多颗卫星之间的伪距时，可能导致伪距观测故障的原因有很多。比如：卫星星钟发生故障，电离层、对流层出现异常，用户接收位置处存在严重多径等。由于卫星数目的增多，比如单一导航系统中搜到了多颗卫星，或者联合使用多个卫星系统进行定位时，包含测距故障的卫星数目占卫星总数的比例较小。这样，采用上述稀疏求解算法求解测距故障，可以降低计算量。

下面以LARS算法为例，说明如何根据第二目标函数g(ε)求解伪距误差向量ε。

对第二目标函数g(ε)求梯度，得到：

$\∂g\left(\mathrm{\ϵ}\right)=S^T(\mathrm{S\ϵ}-\tilde{b})+\mathrm{\λz},\left\{\begin{array}{c}z\left[i\right]=\mathrm{sgn}\left(\mathrm{\ϵ}\right[i\left]\right),\mathrm{\ϵ}\left[i\right]\≠0\\ z\left[i\right]\∈[-\mathrm{1,1}],\mathrm{\ϵ}\left[i\right]=0\end{array}\right.---\left(15\right)$

在（15）式中，当λ→∞，显然ε=0为待求解。因此，求解第二目标函数的问题可以转化为：在λ逐渐减小并趋近的过程中，如何选取ε和z，使得

LARS方法提供了一种构造向量ε和z的方法，描述如下：

1）初始向量ε₀=0，定义解集P={j:ε[j]≠0}。则初始解集P₀为空集，即ε₀中不含非零元素。为利用S中各列通过线性组合恢复向量所得到的残差，c₀=S^Ty₀表示S中各列与向量y₀的相关值；

2）C=max_j(|c₀[j]|)，构造解集P₁={j:|c₀[j]|=C}。则 即在解集P₁下，c₀中相应位置的元素的绝对值均为C。求解向量（Δε中的其他元素为0），使得 令λ=C，z=S^TSΔε。此时，位于解集P₁中的向量有

$|\∂g\left(x_{P_1}\right)|=c_{0,P_1}+{\mathrm{CS}}_{P_1}^T{S}_{P_1}\mathrm{\Δ}{\mathrm{\ϵ}}_{P_1}=c_{0,P_1}-\mathrm{sgn}\left(c_{0,P_1}\right)\left|c_{0,P_1}\right|=0.$

3） $\∀k\∈{P_1}^C,$ 令ε₁=ε₀+γ_kΔε， $y_1={\mathrm{S\ϵ}}_1-\tilde{b},$ $c_1=S^T{y}_1=$ $S^T(S({\mathrm{\ϵ}}_0+{\mathrm{\γ}}_k\mathrm{\Δ\ϵ})-\tilde{b})=c_0+{\mathrm{\γ}}_k{S}^T\mathrm{S\Δ\ϵ}.$ γ_k的选取应满足 $\left|c_{1,P_1}\right|=\left|c_{1,k}\right|,$ 则有

$|c_{0,P_1}+{\mathrm{\γ}}_k{S}_{P_1}^T{S}_{P_1}{\mathrm{\Δ\ϵ}}_{P_1}|=\left(c_0\right[k]+{\mathrm{\γ}}_k{S}_k^T\mathrm{S\Δ\ϵ})1^{\left|P_1\right|},$

即 $|c_{0,P_1}-{\mathrm{\γ}}_k\mathrm{sgn}\left(c_{0,P_1}\right)|=\left|c_0\right[k]+{\mathrm{\γ}}_{k}z[k\left]\right|1^{\left|P_1\right|},$

进而|λ-γ_k|=|c₀[k]+γ_kz[k]|。

4）遍历所有k∈P₁ ^C，按照3）中的公式求解γ_k，选取γ=min⁺ _kγ_k。即从所有大于零的γ_k中选取幅度最小的，并将k添加到解集中，P₂=P₁∪k。此时求解向量（Δε中的其他元素为0），使得令λ=λ-γ，z=S^TSΔε。此时，位于解集P₂中的向量有

从4）中看出，每次向解集中添加元素后，λ将按照步长γ逐渐减小。按照3）、4）中所示的方法，逐步构造向量ε，直到λ减小到

通过上述稀疏求解LARS算法，就可以求解第二目标函数g(ε)，获得伪距误差向量ε。

由于限定了条件伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少，采用了稀疏求解算法，根据本申请的接收机伪距故障监测方法的复杂度降低。

根据本申请的一种实施方式，在利用稀疏求解方法计算伪距误差向量ε前，先对投影矩阵S和和残差向量b进行预处理，并通过预处理后的投影矩阵S和残差向量b重构伪距误差向量ε。

根据一种实施例，对投影矩阵S和残差向量b进行预处理的方法可以包括：计算S₁=WSW^T，W为定位方程中使用加权最小二乘法时的权系数矩阵；计算Z=diag(z)，其中j=1,2,…,n；获得预处理后的投影矩阵算以及预处理后的残差向量

将预处理后的投影矩阵和残差向量带入第二目标函数得到预处理后的第二目标函数：

使用稀疏求解算法，例如LARS算法法对预处理后的第二目标函数进行求解，获得中间向量ε_in，重构伪距误差向量ε=W^TZε_in。

根据伪距误差方程Sε=b，利用稀疏求解方法（例如LARS）计算伪距误差向量ε时，需要满足如下条件：

${\mathrm{\Σ}}_{i=1}^n{S}_{\mathrm{ij}}=0,{\mathrm{\Σ}}_{i=1}^n{S}_{\mathrm{ij}}^2=1,{\mathrm{\Σ}}_{i=1}^n{\tilde{b}}_i=0,j=\mathrm{1,2},...,n---\left(16\right)$

对于（11）式所定义的投影矩阵S，其列向量的和可以表示为：

c=1_1*nS=1_1*n(I-G(G^TWG)^-1G^TW)    （17）

即向量c中的第j个元素表示矩阵S的第j列的和。令：

S₁=WSW^T    （18）

按照（17）式求解S₁的列和，得到：

c₁=1_1*nS₁=1_1*n(WW^T-WG(G^TWG)^-1G^TWW^T)

=1_1*n(I-WG(G^TWG)-¹G^T)WW^T

=(1_1*n-((G(G^TWG)^-1G^TW)1_n*1)^T)WW^T    （19）

上式中最后一个等号利用了W^T=W，可由（7）式中对权矩阵的定义直接得到。从（4）式中给出的几何矩阵G的定义看到，向量1_n*1与矩阵G的最后一列相等，因而由最小二乘解的意义可以得到：

(G(G^TWG)^-1G^TW)1_n*1=1_n*1    （20）

将（20）式代入（19）式，可得c₁=0，即S₁矩阵的列和均为0。

计算j=1,2,…,n。Z=diag(z)为对角阵。记

S₂=S₁Z=WSW^TZ    （21）

则S₂满足各列的列和为零且各列的元素的平方和为1。

考虑如下方程：

$S_2\mathrm{\ϵ}=W\tilde{b}---\left(22\right)$

按照（17）式所示方法计算向量的列和，将（10）式代入，可得

$\mathrm{sum}\left(W\tilde{b}\right)=1_{1*n}(W-\mathrm{WG}{\left(G^T\mathrm{WG}\right)}^{-1}{G}^{T}W)\mathrm{\ϵ}$

$=1_{1*n}(I-\mathrm{WG}{\left(G^T\mathrm{WG}\right)}^{-1}{G}^T)\mathrm{W\ϵ}$

$=(1_{1*n}-{\left(\left(G{\left(G^T\mathrm{WG}\right)}^{-1}{G}^{T}W\right)1_{n*1}\right)}^T)\mathrm{W\ϵ}---\left(28\right)$

由（20）式，可知

综上，方程（22）满足（16）式所述的LARS方法使用条件。下面考察方程（22）与方程（12）的解之间的关系。

将方程（22）展开，得到：

${\mathrm{WSW}}^T\mathrm{Zx}=W\tilde{b}---\left(23\right)$

对照（12）式，可得：

ε=W^TZx    （24）。

根据本申请的一种实施方式，可以进一步根据接收机伪距故障监测结果优化所获得的伪距误差向量ε。

设定优化门限t，选取|ε[i]|＞t[i]的元素ε[i]构成集合C。根据一种实施例，优化门限可以根据伪距测量的标准差σ来选取。例如，δ可以取值为3σ或者6σ。

若集合C为空集，则优化后的伪距误差向量ε等于优化前的伪距误差向量。

否则若集合C为非空集，选取集合C中绝对值最大的k个元素，其中表示下取整，n表示可见卫星数目，记录其在向量ε中的行号，记为集合E。

根据集合E筛选投影矩阵S中的列S_E，计算e=S_E\b。

使用LARS方法对方程求解，得到ε=W^TZx。

根据集合E对向量ε的相应行进行修正，其中ε_newE=ε_E+e，得到优化后的误差向量ε_new。其中，下标E表示根据集合E选取向量中相应的元素，表示E的补集，下标表示根据集合选取出向量中的相应元素。例如向量a=[a₁,a₂,a₃,a₄,a₅]，集合E=[1,4,5]，则a_E=[a₁,a₄,a₅]，

这样，就可以根据优化后的误差向量ε=ε_new，监测接收机伪距故障。

为了实现接收机伪距故障监测，例如，可以预定故障门限；将所计算得到的伪距误差向量ε中的各伪距误差与预定的故障门限进行比较，伪距误差向量ε中的伪距误差高于预定的故障门限所对应的伪距即为故障伪距。

故障门限可以根据实际应用需求进行设定也可根据测距标准差确定。此外，对不同的卫星可以设定不同的门限值。根据一种实施例，当根据实际需求设定时，可以根据需要的定位精度设定故障门限。例如，假定定位偏差不能超过10米，若第k颗卫星发生15米测距偏差会导致定位偏差达到10米，则可以设定该颗卫星的故障门限为15米。根据另一种实施例，故障门限可以根据伪距测量的标准差σ来选取。根据一种实施例，故障门限可以取值为2σ到8σ之间的数值。例如，故障门限可以取值为3σ或者6σ。

通过伪距故障监测，就可以发现例如卫星星钟故障，电离层、对流层异常，用户接收位置处存在严重多径等故障。

根据本申请的一种实施方式，还可以根据计算所得到的伪距误差向量ε对定位结果进行修正，从而在不舍弃故障卫星观测量的情况下，得到修正的定位结果。

初始定位结果为 $x_{\mathrm{ini}}=\left[\begin{array}{c}x\\ y\\ z\\ {\mathrm{\δ}}_u\end{array}\right],$ 重构的伪距误差向量为ε，则修正后的定位结果为x_opt：x_opt=x_ini-Pε。其中，其中，(x,y,z)为接收机的三维位置，δ_u为接收机钟差，矩阵P=(G^TWG)^-1G^TW，G为当前时刻各颗可见卫星的相对位置所决定的几何矩阵，W为各个伪距的加权系数矩阵。

图3显示了根据本申请的一种实施方式的接收机200，能够监测伪距故障。如图3所示，接收机200包括接收模块210，处理模块220，伪距误差向量计算模块230，监测模块240。

接收模块210，获得接收机与多颗卫星之间的伪距，并根据所获得的伪距，计算接收机的位置和钟差；

处理模块220，根据所获得的伪距以及所计算的接收机的位置和钟差，获得卫星几何结构投影矩阵S和残差向量b；

伪距误差向量计算模块230，设定伪距误差容忍度δ，根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少；以及

监测模块240，根据所计算的伪距误差向量ε，监测接收机伪距故障。

根据一种实施方式，伪距误差向量计算模块230可以构建第一目标函数

$f\left(\mathrm{\ϵ}\right)={\left|\right|\mathrm{\ϵ}\left|\right|}_{l_0,\mathrm{\δ}}+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2.$

伪距误差向量计算模块230还可以将近似为‖ε‖₁，并构造第二目标函数

$g\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\δ}}{\left|\right|\mathrm{\ϵ}\left|\right|}_1+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2,$

其中表示向量中元素的绝对值小于δ的个数，‖·‖₂表示向量的l₂范数，‖·‖₁表示向量的l₁范数。‖b‖₁=|b₁|+|b₂|+…|b_n|。

获得第二目标函数后，伪距误差向量计算模块230利用稀疏求解方法通过求解第二目标函数g(ε)的极小值，获得对应的向量ε₀。这样，希望求解的伪距误差向量ε=ε₀。

根据一种实施方式，如图4所示，接收机200可以进一步包括预计算模块250。

预计算模块250将处理模块220获得卫星几何结构投影矩阵S和残差向量b进行预处理，计算S₁=WSW^T，W为各个伪距的加权系数矩阵，计算Z=diag(z)，其中获得预处理后的投影矩阵算并获得预处理后的残差向量

伪距误差向量计算模块230将预计算模块250预处理后的投影矩阵和残差向量带入第二目标函数得到预处理后的第二目标函数： $\tilde{g}\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\δ}}{\left|\right|\mathrm{\ϵ}\left|\right|}_1+\frac{1}{2}{\left|\right|\tilde{b}-\tilde{S}\mathrm{\ϵ}\left|\right|}_2^2.$

伪距误差向量计算模块230利用稀疏求解方法通过求解预处理后的第二目标函数的极小值，获得对应的中间向量ε_in，所述伪距误差向量计算模块230重构伪距误差向量ε=W^TZε_in。

根据一种实施方式，如图4所示，接收机200可以进一步包括伪距误差向量优化模块260，根据接收机伪距故障监测结果，优化所获得的伪距误差向量ε。

伪距误差向量优化模块260设定优化门限t，选取|ε[i]|＞t[i]的元素ε[i]构成集合C。

若集合C为空集，则伪距误差向量优化模块260将优化前的伪距误差向量直接作为优化后的伪距误差向量ε；否则伪距误差向量优化模块260选取集合C中绝对值最大的k个元素，其中表示下取整，n表示可见卫星数目，记录其在向量ε中的行号，记为集合E。

伪距误差向量优化模块260根据集合E筛选投影矩阵S中的列S_E，计算e=S_E\b。

伪距误差向量优化模块260使用稀疏求解LARS算法对方程求解，得到优化中间向量ε_in2=W^TZx。

伪距误差向量优化模块260根据集合E对中间向量ε_in2、的相应行进行修正，其中ε_newE=ε_in2E+e，将ε_newE与组合得到优化后的误差向量ε_new，伪距误差向量优化模块获得优化后的误差向量ε=ε_new。其中，下标E表示根据集合E选取向量中相应的元素，表示E的补集，下标表示根据集合选取出向量中的相应元素。

根据一种实施方式，接收机200的监测模块240可以预定故障门限，将所计算得到的伪距误差向量ε中的各伪距误差与预定的故障门限进行比较，并将伪距误差向量ε中的伪距误差高于预定的故障门限所对应的伪距作为故障伪距。

根据一种实施方式，如图4所示，接收机200的监测模块240可以进一步包括定位修正模块241。定位修正模块241根据所计算的伪距误差向量，修正所计算的接收机的位置和钟差，

x_opt=x_ini-Pε

其中，x_opt表示修正后定位结果，x_ini表示初始定位结果， $X=\left[\begin{array}{c}x\\ y\\ z\\ {\mathrm{\δ}}_u\end{array}\right],$ 其中，(x,y,z)为接收机的三维位置，δ_u为接收机钟差，矩阵P=(G^TWG)^-1G^TW，G为当前时刻各颗可见卫星的相对位置所决定的几何矩阵，W为各个伪距的加权系数矩阵。

以上参照附图对本申请的示例性的实施方案进行了描述。本领域技术人员应该理解，上述实施方案仅仅是为了说明的目的而所举的示例，而不是用来进行限制，凡在本申请的教导和权利要求保护范围下所作的任何修改、等同替换等，均应包含在本申请要求保护的范围内。

Claims (12)

1.一种接收机伪距故障监测方法，包括：

获得接收机与多颗卫星之间的伪距；

根据所获得的伪距，计算接收机的位置和钟差；

根据所获得的伪距以及所计算的接收机的位置和钟差，获得卫星几何结构投影矩阵S和残差向量b；

设定伪距误差容忍度δ，根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少；以及

根据所计算的伪距误差向量ε，监测接收机伪距故障。

2.如权利要求1所述的方法，其中，根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少的步骤包括：

构建第一目标函数

$f\left(\mathrm{\ϵ}\right)={\left|\right|\mathrm{\ϵ}\left|\right|}_{l_0,\mathrm{\δ}}+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2;$

将近似为‖ε‖₁，构造第二目标函数

$g\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\δ}}{\left|\right|\mathrm{\ϵ}\left|\right|}_1+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2;$

利用稀疏求解方法通过求解第二目标函数g(ε)的极小值，获得对应的向量ε₀；

计算伪距误差向量ε=ε₀,

其中表示向量中元素的绝对值小于δ的个数，‖·‖₂表示向量的l₂范数，‖·‖₁表示向量的l₁范数。‖b‖₁=|b₁|+|b₂|+…|b_n|。

3.如权利要求2所述的方法，其中，所述方法进一步包括：

在利用稀疏求解方法计算伪距误差向量ε前，对投影矩阵S和残差向量b进行预处理：计算S₁=WSW^T，W为各个伪距的加权系数矩阵，计算Z=diag(z)，其中j=1,2,…,n，获得预处理后的投影矩阵算并获得预处理后的残差向量

将预处理后的投影矩阵和残差向量带入第二目标函数得到预处理后的第二目标函数： $\tilde{g}\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\δ}}{\left|\right|\mathrm{\ϵ}\left|\right|}_1+\frac{1}{2}{\left|\right|\tilde{b}-\tilde{S}\mathrm{\ϵ}\left|\right|}_2^2;$

利用稀疏求解方法通过求解预处理后的第二目标函数的极小值，获得对应的中间向量ε_in；

重构伪距误差向量ε=W^TZε_in。

4.如权利要求3所述的方法，其中，所述方法进一步包括根据接收机伪距故障监测结果，优化所获得的伪距误差向量ε：

设定优化门限t，选取|ε[i]|＞t[i]的元素ε[i]构成集合C；

若集合C为空集，则优化后的伪距误差向量ε等于优化前的伪距误差向量；否则选取集合C中绝对值最大的k个元素，其中表示下取整，n表示可见卫星数目，记录其在向量ε中的行号，记为集合E；

根据集合E筛选投影矩阵S中的列S_E，计算e=S_E\b；

使用稀疏求解LARS算法对方程求解，得到优化中间向量ε_in2=W^TZx；

根据集合E对中间向量ε_in2的相应行进行修正，其中ε_newE=ε_in2E+e，将ε_newE与组合得到优化后的误差向量ε_new，则优化后的误差向量ε=ε_new，

其中，下标E表示根据集合E选取向量中相应的元素，表示E的补集，下标表示根据集合选取出向量中的相应元素。

5.如权利要求1-4中任一项所述的方法，其中，所述根据所计算的伪距误差向量，监测接收机伪距故障的步骤包括：

预定故障门限，将所计算得到的伪距误差向量ε中的各伪距误差与预定的故障门限进行比较，伪距误差向量ε中的伪距误差高于预定的故障门限所对应的伪距即为故障伪距。

6.如权利要求1-4中任一项所述的方法，其中，所述方法进一步包括：根据所计算的伪距误差向量，修正所计算的接收机的位置和钟差，

x_opt=x_ini-Pε

其中，x_opt表示修正后定位结果，x_ini表示初始定位结果， $X=\left[\begin{array}{c}x\\ y\\ z\\ {\mathrm{\δ}}_u\end{array}\right],$ 其中，(x,y,z)为接收机的三维位置，δ_u为接收机钟差，矩阵P=(G^TWG)^-1G^TW，G为当前时刻各颗可见卫星的相对位置所决定的几何矩阵，W为各个伪距的加权系数矩阵。

7.一种接收机，能够监测伪距故障，所述接收机包括：

接收模块，获得接收机与多颗卫星之间的伪距，并根据所获得的伪距，计算接收机的位置和钟差；

处理模块，根据所获得的伪距以及所计算的接收机的位置和钟差，获得卫星几何结构投影矩阵S和残差向量b；

伪距误差向量计算模块，设定伪距误差容忍度δ，根据伪距误差方程Sε=b，计算伪距误差向量ε，使得所述伪距误差向量ε中的各伪距误差高于所述伪距误差容忍度δ的个数最少；以及

监测模块，根据所计算的伪距误差向量ε，监测接收机伪距故障。

8.如权利要求7所述的接收机，其中，

所述伪距误差向量计算模块构建第一目标函数

$f\left(\mathrm{\ϵ}\right)={\left|\right|\mathrm{\ϵ}\left|\right|}_{l_0,\mathrm{\δ}}+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2;$

伪距误差向量计算模块将近似为‖ε‖₁，并构造第二目标函数

$g\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\δ}}{\left|\right|\mathrm{\ϵ}\left|\right|}_1+\frac{1}{2}{\left|\right|b-\mathrm{S\ϵ}\left|\right|}_2^2;$

伪距误差向量计算模块利用稀疏求解方法通过求解第二目标函数g(ε)的极小值，获得对应的向量ε₀；以及

伪距误差向量计算模块计算伪距误差向量ε=ε₀,

其中表示向量中元素的绝对值小于δ的个数，‖·‖₂表示向量的l₂范数，‖·‖₁表示向量的l₁范数。‖b‖₁=|b₁|+|b₂|+…|b_n|。

9.如权利要求8所述的接收机，其中，所述接收机进一步包括预计算模块，其中，

所述预计算模块将所述处理模块获得卫星几何结构投影矩阵S和残差向量b进行预处理，计算S₁=WSW^T，W为各个伪距的加权系数矩阵，计算Z=diag(z)，其中j=1,2,…,n，获得预处理后的投影矩阵算并获得预处理后的残差向量

所述伪距误差向量计算模块将所述预计算模块预处理后的投影矩阵和残差向量带入第二目标函数得到预处理后的第二目标函数： $\tilde{g}\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\δ}}{\left|\right|\mathrm{\ϵ}\left|\right|}_1+\frac{1}{2}{\left|\right|\tilde{b}-\tilde{S}\mathrm{\ϵ}\left|\right|}_2^2;$

所述伪距误差向量计算模块利用稀疏求解方法通过求解预处理后的第二目标函数的极小值，获得对应的中间向量ε_in；以及

所述伪距误差向量计算模块重构伪距误差向量ε=W^TZε_in。

10.如权利要求9所述的接收机，其中，所述接收机进一步包括伪距误差向量优化模块，根据接收机伪距故障监测结果，优化所获得的伪距误差向量ε，其中，

所述伪距误差向量优化模块设定优化门限t，选取|ε[i]|＞t[i]的元素ε[i]构成集合C；

若集合C为空集，则所述伪距误差向量优化模块将优化前的伪距误差向量直接作为优化后的伪距误差向量ε；否则所述伪距误差向量优化模块选取集合C中绝对值最大的k个元素，其中表示下取整，n表示可见卫星数目，记录其在向量ε中的行号，记为集合E；

所述伪距误差向量优化模块根据集合E筛选投影矩阵S中的列S_E，计算e=S_E\b；

所述伪距误差向量优化模块使用稀疏求解LARS算法对方程求解，得到优化中间向量ε_in2=W^TZx；

所述伪距误差向量优化模块根据集合E对中间向量ε_in2、的相应行进行修正，其中ε_newE=ε_in2E+e，将ε_newE与组合得到优化后的误差向量ε_new，所述伪距误差向量优化模块获得优化后的误差向量ε=ε_new，

其中，下标E表示根据集合E选取向量中相应的元素，表示E的补集，下标表示根据集合选取出向量中的相应元素。

11.如权利要求7-10中任一项所述的接收机，其中，所述监测模块预定故障门限，将所计算得到的伪距误差向量ε中的各伪距误差与预定的故障门限进行比较，并将伪距误差向量ε中的伪距误差高于预定的故障门限所对应的伪距作为故障伪距。

12.如权利要求7-10中任一项所述的接收机，其中，所述接收机的监测模块进一步包括定位修正模块：

所述定位修正模块根据所计算的伪距误差向量，修正所计算的接收机的位置和钟差，

x_opt=x_ini-Pε

其中，x_opt表示修正后定位结果，x_ini表示初始定位结果， $X=\left[\begin{array}{c}x\\ y\\ z\\ {\mathrm{\δ}}_u\end{array}\right],$ 其中，(x,y,z)为接收机的三维位置，δ_u为接收机钟差，矩阵P=(G^TWG)^-1G^TW，G为当前时刻各颗可见卫星的相对位置所决定的几何矩阵，W为各个伪距的加权系数矩阵。