CN101799524A - 全球导航卫星系统接收机自主完备性监测的方法 - Google Patents

Abstract

本发明公开了一种全球导航卫星系统接收机自主完备性监测的方法，包括：接收卫星星历及卫星观测数据，作为全球导航卫星系统GNSS定位的原始数据；获取卫星伪距测量误差方程的系数阵G阵和常数项，根据卫星高度角获取权阵；根据误差方程和卫星权阵构建法方程；根据最小二乘法进行定位求解；根据卫星数和相应的弃真概率和纳伪概率进行粗差检验；若存在粗差，进行粗差识别：将误差方程的系数阵G进行QR分解，得到相应的奇偶检校阵P和p；对P和p进行矩阵转换，得到新的矩阵A；计算矩阵A各向量之间的距离，组成相关距离矩阵B；逐一对矩阵B中各个向量距离进行聚类分析，判别向量之间的亲疏关系，从而识别粗差；将探测出的粗差剔除，回代误差方程，重新进行定位，供用户使用。

Description

全球导航卫星系统接收机自主完备性监测的方法

技术领域

本发明涉及卫星导航领域，特别是在全球导航卫星系统中对用户接收机进行自主完备性监测的方法。

背景技术

导航是利用某种方法或手段，引导交通工具或其他运动物体从一个位置移动到另一个位置的过程。当前常用的导航方法主要包括惯性导航、无线电导航和卫星导航等方法。卫星导航是通过接收导航卫星发送的导航定位信号，并将导航卫星作为动态已知点，为运动载体实时提供全球、全天候、高精度的位置、速度和时间信息，进而完成各种导航任务。卫星导航定位系统从最初的多普勒卫星导航系统发展到了全球导航卫星系统GNSS(GlobalNavigation Satellite System，简称GNSS)，其中美国的GPS(Global PositioningSystem)和前苏联的GLONASS(GLObal NAvigation Satellite System)是第二代卫星导航定位系统的代表。COMPASS北斗导航系统是中国自主开发的导航定位系统。全球导航卫星系统GNSS可以在全世界的陆海空范围内实现实时连续高精度的位置速度以及时间等导航信息，在各类用户中得到了广泛应用。

GNSS的指标可用四个参数表示：精度、持续性、有效性和完备性。精度是指任意时刻估计位置与真实位置的符合程度，即系统的真实精度或绝对精度。连续性是指导航系统在给定的使用条件下，在规定的时间内，以规定的性能完成其成功率的概率。有效性是指系统能为运动载体提供可用的导航服务的时间的概率。完备性是指导航系统发生任何故障或者误差超限时，无法用于导航和定位时，系统向用户及时发出报警的能力。

系统的精度可以根据需要和可能进行调整与控制，但系统的完备性信息在任何时候都是不可少的，它是用户对系统所提供信息的可信程度的一种度量，包括系统给用户提供及时有效的警告信息的能力，是保证用户安全性的重要参数。

用户自主式完备性监测(Receiver Autonomous Integrity Monitoring，缩写为RAIM)，是指根据用户接收机的多余观测值监测用户定位结果的完备性。最典型的情况是发现一颗有问题的卫星对定位结果的影响。卫星的问题可能是钟差改正信息不正确或者卫星不能正常工作等。

在现有的RAIM算法方面，比较普遍的是利用距多余观测量实现粗差卫星的探测和识别。由于用户站没有已知的位置坐标，只能通过其单位权中误差来进行粗差的探测，并对残差进行分析，判断粗差存在于哪颗卫星的观测值上。目前常用的方法有如下几种：

一、最小二乘残差法

用户站的基本观测方程为：

L＝Gδx+εW  (W是权)

则可得到单位权方差为：

${\widehat{\mathrm{\σ}}}_0^2=\frac{V^T\mathrm{WV}}{n-4}$

当在给定的置信度水平α下大于

的限值时，则说明观测值中含有粗差。

粗差的探测原则如下：

对GPS定位而言，残差V服从正态分布，V～N(0，σ²)，则可有V^TPV服从自由度为(n-4)的χ²分布，即

${\mathrm{\χ}}^2=\frac{V^T\mathrm{PV}}{{\mathrm{\σ}}_0^2}=\frac{(n-4){\mathrm{\σ}}^2}{{\mathrm{\σ}}_0^2}~{\mathrm{\χ}}^2(n-4)$

其中

为已知的先验值。

对统计量σ²进行χ²检验，

设原假设

备选假设

因

$P\{\frac{V^T\mathrm{PV}}{{\mathrm{\σ}}_0^2}>{\mathrm{\χ}}_{\mathrm{\α}}^2(n-4)\}=\mathrm{\α},$

或

P{σ²＜k}＝α，

式中：

$k=\frac{{\mathrm{\χ}}_{\mathrm{\α}}^2(n-4){\mathrm{\σ}}_0^2}{n-4},$

当σ²＞k时，则拒绝H₀，接受H₁；否则，接受H₀。

二、Parity方法

在Parity方法中，由差分单点定位基本观测方程和误差方程：

L＝Gδx+ε

V＝Gδx-L    W(权)

可得：

δx＝(G^TWG)^-1 G^TWL

对上式中的L向量做一线性转换。

$\left[\begin{array}{c}\mathrm{\δx}\\ ...\\ p\end{array}\right]=\left[\begin{array}{c}{\left(G^T\mathrm{WG}\right)}^{-1}{G}^{T}W\\ ............\\ P\end{array}\right]L$

则可得Parity向量定义式：

p＝PL

Parity方法的粗差分离如下：

由Parity向量定义式p＝PL，并把L用它的等量(GδX+ε)来代替，由于正交性质PG＝0，故有：

p＝Pε

假设忽略观测值噪声，在某颗卫星上考虑存在一偏差b，若偏差b存在某颗卫星上，则其在Parity空间的投影结果是：

$p(\mathrm{for\; bias\; on\; sat}.4)=\left[\begin{array}{c}{p}_{14}\\ p_{24}\end{array}\right]b$

由作用在某颗卫星上的偏差引起的Parity向量，肯定位于斜率是p₂₄/p₁₄的一条直线(特征偏差线)上。图1是Parity空间中显示的由作用在某颗卫星上的偏差引起的Parity向量p。

故此，每一颗卫星都有自己的特征偏差线，其斜率由向量P的各列的元素决定，即：

卫星的特征偏差线斜率＝p_2i/p_1i，i＝1，2，…，6

则可得法则：含有粗差的卫星就是那颗特征偏差线与观测的Parity向量p重合的卫星。

三、Parity相关系数方法

由Parity方法， $\underset{(n-4)\×1}{p}=\underset{(n-4)\×n}{p}\underset{n\×1}{\mathrm{\ϵ}}$

矩阵展开式：

$\mathrm{P\ϵ}=\left[\begin{array}{cccc}{P}_{11}& P_{12}& ...& P_{1n}\\ P_{21}& P_{22}& ...& P_{2n}\\ ...& ...& ...& ...\\ P_{m1}& P_{m2}& ...& P_{\mathrm{mn}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ ...\\ {\mathrm{\ϵ}}_n\end{array}\right]=\left[\begin{array}{c}{P}_{11}\\ P_{21}\\ ..\\ P_{m1}\end{array}\right]{\mathrm{\ϵ}}_1+\left[\begin{array}{c}{P}_{21}\\ P_{22}\\ ...\\ P_{m2}\end{array}\right]{\mathrm{\ϵ}}_2+...+\left[\begin{array}{c}{P}_{1n}\\ P_{2n}\\ ...\\ P_{\mathrm{mn}}\end{array}\right]{\mathrm{\ϵ}}_n=\left[\begin{array}{c}{p}_1\\ p_2\\ ...\\ p_m\end{array}\right]$

其中m＝n-4，n为卫星数；ε_i(i＝1，2，...，n)是观测值i的观测误差，是数值变量。

设：P_i＝[P_1i P_2i...P_mi]^T(i＝1，2，...n)

则有：

ε₁P₁+ε₂P₂+...+ε_nP_n＝p

可称P_i为观测量i的误差对Parity向量p的影响向量。

P_i由卫星位置的几何图形矩阵决定，ε_i由观测量的函数特性所决定，因此ε_iP_i是由观测量的误差和卫星几何图形矩阵共同确定的。当某观测量出现粗差时，都将表现为该观测量的||ε_iP_i||值较大，即向量ε_iP_i的模较大，并在左半部分各项中占有优势，即所占份额比较大，从而右半部分主要受到该有粗差的ε_iP_i的影响，所以，Parity向量p与误差观测量的影响向量P_i具有相当强的相关性。

现表征P_i与p相关的数值指标为

${\mathrm{\ρ}}_{P_i,p}=\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}(P_{\mathrm{ji}}-{\stackrel{\‾}{P}}_i)(p_j-\stackrel{\‾}{p})}{[{\left(\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{(P_{\mathrm{ji}}-{\stackrel{\‾}{P}}_i)}^2\right)\left(\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{(p_j-\stackrel{\‾}{p})}^2\right)]}^{\frac{1}{2}}}$

称

为P_i与p的相关系数。

是一数值量，无量纲，它反映了P_i与p的相关程度。

当观测量ε中含有一个粗差时，

越接近于1，相关性越强，说明Parity向量p受到观测量i的误差ε_i的影响越大，则其出现粗差的可能性就越大；而

越接近于0，观测量i的误差ε_i对Parity向量p的影响小，则其出现粗差的可能性就越小。当观测量ε中含有多个粗差时，Parity向量p的变化规律，将表现为来自这些粗差对其影响的叠加。

以上这些方法都是基于伪距观测测量的测量误差，通过概率和矩阵变换的方法来进行粗差的识别和探测。发明人在研究和应用的过程中发现：这些RAIM算法中主要涉及到粗差检验以及单个粗差的识别和探测，虽然相关系数法可以发现多个粗差，但相关系数本身具有相关性，还无法真正做到多个粗差的识别和探测，本发明创造寻求采用相关距离和聚类分析的方法实现多个粗差的识别和探测。

发明内容

有鉴于此，本发明的目的在于提供一种全球导航卫星系统接收机自主完备性监测的方法，用于实现多个粗差的识别和探测，提高用户使用的可靠性和可用性。

本发明的实施例提供了一种全球导航卫星系统接收机自主完备性监测的方法，包括：

接收卫星星历及卫星观测数据，作为全球导航卫星系统GNSS定位的原始数据；

获取卫星伪距测量误差方程的系数阵G阵和常数项，根据卫星高度角获取权阵；

根据误差方程和卫星权阵构建法方程；

根据最小二乘法进行定位求解；

根据卫星数和相应的弃真概率和纳伪概率进行粗差检验；

若存在粗差，进行粗差识别：将误差方程的系数阵G进行QR分解，得到相应的奇偶检校阵P和p；对P和p进行矩阵转换，得到新的矩阵A；

计算矩阵A各向量之间的距离，组成相关距离矩阵B；

逐一对矩阵B中各个向量距离进行聚类分析，判别向量之间的亲疏关系，从而识别粗差；

将探测出的粗差剔除，回代误差方程，重新进行定位，供用户使用。

本发明将误差方程系数阵进行QR分解，形成奇偶检矩阵，并在矩阵变换的基础上进行相关距离的计算，利用聚类分析的方法识别多个粗差，克服了现有技术中RAIM算法由于相关系数自身相关的缺陷而使得粗差探测失真的技术问题，该方法具有如下优点：

1、利用矩阵变换计算相似距离，方法更加简便；

2、利用聚类分析计算矩阵各向量之间的亲疏关系，使得粗差能够更加清晰的识别出来；

3、本方法可以应用在GNSS技术上，既可以利用在多种卫星导航系统的组合系统上面，同时可以将GNSS与其他导航辅助设施共同应用，如组合惯性导航系统，组合气压计等，更加便于导航用户使用。

附图说明

图1是现有技术中Parity空间中显示的由作用在某颗卫星上的偏差引起的Parity向量p的示意图；

图2为本发明实施例中GNSS系统接收机自主完备性监测的方法的流程图。

具体实施方式

本发明实施例采用相关距离和聚类分析的方法实现多个粗差的识别和探测。聚类分析是一种研究事物分类的多元统计方法，它可避免以往对事物分类主要靠经验和定性知识进行所带来的主观性、任意性等缺点，正确地反映事物的内在联系，并能定量地表达出事物之间的相似性和差异性，是一种定量的、客观的、科学的分类方法。

假设若干样品(某种事物)，其性质用选定的某些指标(属性、特征)加以刻划，聚类分析根据分类目的的不同分为两种：一种是对样品分类，称为Q型聚类分析；另一种对指标分类，称为R型聚类分析。

以Q型聚类分析为例，就是根据样品之间的相似性或亲疏程度，用数学方法把它们逐步地分型划类，最后得到一个能够反映样品之间亲疏关系的客观的分类系统。聚类分析最后能够得到一张反映样品之间亲疏关系的聚类图，比较自然、客观地描述了分类对象的各个样品之间的差异和联系。

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明作进一步的详细描述。

图2是本实施例提供的GNSS系统接收机自主完备性监测的方法的流程图，该方法具体包括：

步骤201、接收卫星星历，接收卫星观测数据，作为GNSS定位的原始数据。GNSS导航用户利用接收机接收卫星星历，即可计算卫星位置，由于卫星星历是每两个小时发播一次，可以采用内插的方法得到瞬时的卫星位置；用户同时接收卫星观测数据，卫星的观测数据一般包括伪距观测值、精码观测值以及载波相位观测值等，可以根据定位要求的不同使用不同的数据，从而得到不同的定位精度。卫星星历和卫星观测数据都作为GNSS定位的原始数据。

步骤202、获取卫星伪距测量误差方程的系数阵G阵和常数项，根据卫星高度角计算权阵。

根据卫星星历计算的卫星和用户接收机的各种观测数据，计算卫星伪距测量误差方程的系数阵G阵和常数项：

GPS单站单星的C/A码伪距观测方程为：

${\mathrm{\ρ}}_k^j\left(i\right)=|{\stackrel{\→}{X}}_s\left(t^j\left(\mathrm{GPS}\left(i\right)\right)\right)-{\stackrel{\→}{X}}_r|+C(\mathrm{\Δ}{t}_k\left(i\right)-\mathrm{\Δ}{t}^j\left(i\right))$

$+\mathrm{C\Δ}{t}_w+\mathrm{\Δ}{t}_{\mathrm{ion}}+\mathrm{\Δ}{t}_{\mathrm{tro}}+\mathrm{\Δ}{t}_{\mathrm{rel}}$

其中：

为测站k在时刻i对卫星j的伪距观测值；

为信号发射时刻的卫星坐标，

为测站坐标；Δt_k(i)为接收机钟差，它是接收机钟面时t_k(i)与GPS标准时的差值，Δt_k(i)＝t_k(i)-t_k(GPS(i))；Δt^j(i)是卫星钟差，它为接收信号离开卫星时，卫星钟面时t^j(i)与GPS标准时的差值，包括卫星钟的物理钟差和SA频率抖动产生的钟差，Δt^j(i)＝t^j(i)-t^j(GPS(i))；C为光速；CΔt_w、Δt_ion、Δt_tro、Δt_rel分别为地球自转改正、电离层改正、对流层改正和相对论延迟。

根据伪距观测方程式，伪距误差方程为：

$v^j=|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_r|-C\·\mathrm{\Δt}+\underset{k=1}{\overset{5}{\mathrm{\Σ}}}\mathrm{\Δ}{\mathrm{\ρ}}_k-{\mathrm{\ρ}}^j$

线性化后并由一个历元上的n个观测值组成的误差方程为：

V＝Gδx-L

其中： $G=\left[\begin{array}{cccc}{l}^1\left(t\right)& m^1\left(t\right)& n^1\left(t\right)& -1\\ l^2\left(t\right)& m^2\left(t\right)& n^{2\left(t\right)}& -1\\ ......& ......& & \\ l^n\left(t\right)& m^n\left(t\right)& n^n\left(t\right)& -1\end{array}\right]$

$l\left(t\right)=\frac{X_s-X_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$ $m\left(t\right)=\frac{Y_s-Y_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$ $n\left(t\right)=\frac{Z_s-Z_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$

δx＝[δx_r1，δx_r2，δx_r3，CΔt]^T，CΔt为以米为单位的接收机钟差。

从上述GPS的C/A码误差方程推广到GNSS导航系统，若GNSS包括GPS/GLONASS/GALILEO/COMPASS，每个系统在地面接收机各自包含一个接收机钟差，则GNSS一个历元上的n个观测值组成的误差方程为：

V＝Gδx-L

其中： $G=\left[\begin{array}{ccccccc}{l}^1\left(t\right)& m^1\left(t\right)& n^1\left(t\right)& -E& 0& 0& 0\\ l^2\left(t\right)& m^2\left(t\right)& n^2\left(t\right)& 0& -E& 0& 0\\ l^3\left(t\right)& m^3\left(t\right)& n^3\left(t\right)& 0& 0& -E& 0\\ l^4\left(t\right)& m^4\left(t\right)& n^4\left(t\right)& 0& 0& 0& -E\end{array}\right]$

$l\left(t\right)=\frac{X_s-X_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$ $m\left(t\right)=\frac{Y_s-Y_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$ $n\left(t\right)=\frac{Z_s-Z_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$

δx＝[δx_r1，δx_r2，δx_r3，CΔt₁，CΔt₂，CΔt₃，CΔt₄]^T

其中lⁱ(t)，mⁱ(t)，nⁱ(t)分别代表不同的导航系统的站星卫星结构系数，如l¹(t)，m¹(t)，n¹(t)代表GPS系统，l²(t)，m²(t)，n²(t)代表GLONASS系统，l³(t)，m³(t)，n³(t)代表GALILEO系统，l⁴(t)，m⁴(t)，n⁴(t)代表COMPASS系统等。

根据卫星高度角计算权阵：

权W的定权公式：根据卫星高度角定权，设高度角为90°时权最大，为W0；高度角为5°时权最小，为1。为了简化计算，假定权随高度角线性变化，则高度角为E的观测值权为：

$W=\frac{(W_0-1)E+90-5W_0}{85}$

W₀的值可以取3或4。

步骤203、根据误差方程和卫星权阵构建法方程。根据误差方程以及相应的卫星权阵，可以构建法方程，得到法方程如下：

(G^TWG)δx＝G^TWL

步骤204、根据最小二乘法进行定位求解。根据法方程，可以按照最小二乘法求解用户站位置，其中最小二乘解及其精度为：

δx＝(G^TWG)^-1 G^TWL

$D_x={\widehat{\mathrm{\σ}}}_0^2{\left(A^T\mathrm{WG}\right)}^{-1}$

${\widehat{\mathrm{\σ}}}_0^2=\frac{L^T\mathrm{WL}+L^T\mathrm{WG\δx}}{n-7}$

步骤205、根据卫星数和相应的弃真概率和纳伪概率进行粗差检验。根据卫星数和相应的弃真概率α和纳伪概率β进行粗差检验：

设原假设

备选假设

因

$P\{\frac{V^T\mathrm{PV}}{{\mathrm{\σ}}_0^2}>{\mathrm{\χ}}_{\mathrm{\α}}^2(n-4)\}=\mathrm{\δ}(\mathrm{\α},\mathrm{\β}),$

或

P{σ²＜k}＝δ(α，β)，

式中：

$k=\frac{{\mathrm{\χ}}_{\mathrm{\α},\mathrm{\β}}^2(n-4){\mathrm{\σ}}_0^2}{n-4},$

当σ²＞k时，则拒绝H₀，接受H₁；否则，接受H₀。

步骤206、若存在粗差，进行粗差识别：将误差方程的系数阵G进行QR分解，得到相应的奇偶检校阵P和p；对P和p进行矩阵转换，得到新的矩阵A。

根据步骤205，如果判断存在粗差的话，就需要将粗差识别并予以处理。此时可以将误差方程的系数阵G进行QR分解，得到相应的奇偶检校阵P和p。

在最小二乘方法中，由上述单点定位基本观测方程和误差方程：

L＝Gδx+ε

V＝Gδx-L    W(权)

则用户定位未知参数方程：

δx＝(G^TWG)^-1G^TWL

此时可对L向量做一线性转换，

$\left[\begin{array}{c}\mathrm{\δx}\\ ...\\ p\end{array}\right]=\left[\begin{array}{c}{\left(G^T\mathrm{WG}\right)}^{-1}{W}^{T}P\\ ............\\ P\end{array}\right]L$

则可得QR奇偶检效法向重定义式：

p＝PL

QR奇偶检效法向量p_(n-4)×1的性质：

E[p]＝0

E[pp^T]＝Cov[p]＝σ²I

QR奇偶检效法向量转换矩阵P是一特殊转换，它把n维观测空间转化为(n-4)维的Parity空间。P有以下特殊性质：(1)P的行与G的列正交；(2)P的行相互正交；(3)P的行都进行了标准化，每一行的大小都是单位1。

在QR奇偶检效法方法中，经证明可得：

p^Tp＝V^TV

此时对P和p进行矩阵转换，得到新的矩阵A；

由p＝Pε可得其矩阵展开式：

$\mathrm{P\ϵ}=\left[\begin{array}{cccc}{P}_{11}& P_{12}& ...& P_{1n}\\ P_{21}& P_{22}& ...& P_{2n}\\ ...& ...& ...& ...\\ P_{m1}& P_{m2}& ...& P_{\mathrm{mn}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ ...\\ {\mathrm{\ϵ}}_{\mathrm{nm}}\end{array}\right]=\left[\begin{array}{c}{P}_{11}\\ P_{21}\\ ...\\ P_{m1}\end{array}\right]{\mathrm{\ϵ}}_1+\left[\begin{array}{c}{P}_{21}\\ P_{22}\\ ...\\ P_{m2}\end{array}\right]{\mathrm{\ϵ}}_2+...+\left[\begin{array}{c}{P}_{1n}\\ P_{2n}\\ ...\\ P_{\mathrm{mn}}\end{array}\right]\mathrm{\ϵ}=\left[\begin{array}{c}{p}_1\\ p_2\\ ...\\ p_m\end{array}\right]$

其中m＝n-7，n为卫星数；ε_i(i＝1，2，...，n)是观测值i的观测误差，是数值变量。

设：P_i＝[P_1i P_2i...P_mi]^T(i＝1，2，...n)

则有：

ε₁P₁+ε₂P₂+...+ε_nP_n＝p

进行矩阵形式的转换，可以得到：

$\left[\begin{array}{cccc}{P}_{11}{\mathrm{\ϵ}}_1& P_{21}{\mathrm{\ϵ}}_2& P_{1n}{\mathrm{\ϵ}}_n& p_1\\ P_{21}{\mathrm{\ϵ}}_1& P_{22}{\mathrm{\ϵ}}_2& P_{2n}{\mathrm{\ϵ}}_n& p_2\\ ...& ...& ...& ...\\ P_{m1}{\mathrm{\ϵ}}_1& P_{m2}{\mathrm{\ϵ}}_2& P_{\mathrm{mn}}{\mathrm{\ϵ}}_n& p_m\end{array}\right]=\left[\begin{array}{ccccc}{a}_1& a_2& ...& a_n& a_{n+1}\end{array}\right]=A$

即将左边部分与右边部分合并处理。

步骤207、计算矩阵A各向量之间的距离，组成相关距离矩阵B。根据步骤206得到的矩阵A，计算各个向量a间相关距离。

相关距离是衡量样品间亲疏关系或相似程度的统计量。将样品看做m维空间的点，用点间的距离来表示样品之间的亲疏关系。距离越小，表明样品之间的关系越密切；反之，距离越大，相似程度越差。

相关距离可以采用如下公式：

$d_{\mathrm{ij}}^{\′}\left(2\right)={\left[\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(x_{\mathrm{ik}}-x_{\mathrm{jk}})}^2\right]}^{\frac{1}{2}}$

为a_i与a_j的相关距离，反映了各个向量之间的相关程度。可以采用

数值指标，判别粗差存在于哪个观测值上。

假设当观测量ε_i中含有一个粗差时，则该列向量a_i与最后一列向量a_n+1的相关距离

会明显变小，说明相关性越强，说明最后一列a_n+1受到观测量i的误差ε_i的影响越大，则其出现粗差的可能性就越大；而其他列向量a_i与最后一列向量a_n+1的相关距离

会变大，且具有同向性。当观测量ε中含有多个粗差时，相关距离

的变化规律，将表现为来自这些粗差对其影响的叠加。

计算矩阵A各向量之间的距离，组成相关距离矩阵B；

$\left[\begin{array}{cccc}{b}_{11}& b_{21}& b_{m1}& b_{(n+1)1}\\ b_{12}& b_{22}& b_{m2}& b_{(n+1)2}\\ ...& ...& ...& ...\\ b_{1(n+1)}& b_{2(n+1)}& b_{m(n+1)}& b_{(n+1)(n+1)}\end{array}\right]=B$

其中 $b_{\mathrm{ij}}={\left[\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(a_{\mathrm{ik}}-a_{\mathrm{jk}})}^2\right]}^{\frac{1}{2}}.$

步骤208、根据聚类分析的方法，逐一对各个向量距离进行归类，判别向量之间的亲疏关系，从而识别粗差。

聚类方法的基本思想是：先将n个样品各自看成一类，然后规定样品之间的距离和类与类之间的距离；开始，因每个样品自成一类，类与类之间的距离与样品之间的距离是相等的，选择距离最小的一对并成一个新类，计算新类和其它类的距离，再将距离最近的两类合并，这样每次减少一类，直至所有的样品都成一类为止。本实施例采用最小距离法进行粗差与偶然误差的分类。

令d_ij表示样品i和样品j的距离，G1、G2...表示类，D_pq表示类G_p与类G_q之间的距离。

类与类之间的距离定义为两类之间最近样品的距离：

$D_{\mathrm{pq}}=\underset{i\∈G_p,j\∈G_q}{\min }{d}_{\mathrm{ij}}$

用最短距离法聚类的步骤如下：

1)规定样品之间的距离，计算样品两两之间的距离，得到一个对称矩阵，记作D(0)，开始每个样品自成一类，这时显然D_pq＝d_pq。

2)定义类与类之间的距离为两类之间最近样品的距离。

3)选择D(0)的最小元素，设为D_pq，则将G_p和G_q合并成一新类，记为G_r，G_r＝{G_p G_q}。

4)计算新类与其它类的距离：

$D_{\mathrm{rk}}=\underset{i\∈G_r,j\∈G_k}{\min }{d}_{\mathrm{ij}}=\min \{\underset{i\∈G_p,j\∈G_k}{\min }{d}_{\mathrm{ij}},\underset{i\∈G_q,j\∈G_k}{\min }{d}_{\mathrm{ij}}\}$

$=\min \{D_{\mathrm{pk}},D_{\mathrm{qk}}\}$

将D(0)中p、q行，p、q列并成一个新行新列，新行新列对应G_r，所得到的矩阵记作D(1)。

5)对D(1)重复上述对D(0)所采用的两个步骤3)和4)，得到D(2)，如此执行下去直到所有的样品都成为一类为止。

根据上述步骤，对计算的各类相关距离矩阵B进行聚类分析，可以得到最后一列向量与其他向量的亲疏关系，从而识别粗差。

步骤209、将探测出的粗差剔除，回代误差方程，重新进行定位，供用户使用。

为了防止发生弃真错误和纳伪错误，将探测出的粗差集剔除，然后将不含粗差的卫星回代误差方程，重新进行定位，然后继续进行粗差检验，若检验通过，即认为已经不含有粗差，定位结果可供用户放心使用。

本实施例将误差方程系数阵进行QR分解，形成奇偶检矩阵，并在矩阵变换的基础上进行相关距离的计算，利用聚类分析的方法识别多个粗差，克服了现有技术中RAIM算法由于相关系数自身相关的缺陷而使得粗差探测失真的技术问题。

总之，以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (7)

1.一种全球导航卫星系统接收机自主完备性监测的方法，其特征在于，包括：

接收卫星星历及卫星观测数据，作为全球导航卫星系统GNSS定位的原始数据；

获取卫星伪距测量误差方程的系数阵G阵和常数项，根据卫星高度角获取权阵；

根据误差方程和卫星权阵构建法方程；

根据最小二乘法进行定位求解；

根据卫星数和相应的弃真概率和纳伪概率进行粗差检验；

若存在粗差，进行粗差识别：将误差方程的系数阵G进行QR分解，得到相应的奇偶检校阵P和p；对P和p进行矩阵转换，得到新的矩阵A；

计算矩阵A各向量之间的距离，组成相关距离矩阵B；

逐一对矩阵B中各个向量距离进行聚类分析，判别向量之间的亲疏关系，从而识别粗差；

将探测出的粗差剔除，回代误差方程，重新进行定位，供用户使用。

2.根据权利要求1所述的方法，其特征在于，所述接收卫星星历及卫星观测数据具体包括：

GNSS用户利用接收机接收卫星星历，采用内插的方法得到瞬时的卫星位置；用户同时接收卫星观测数据，包括：伪距观测值、精码观测值以及载波相位观测值。

3.根据权利要求1所述的方法，其特征在于，所述GNSS包括全球定位系统GPS、全球导航卫星系统GLONASS、伽俐略系统GALILEO及北斗系统COMPASS。

4.根据权利要求1、2或3所述的方法，其特征在于，所述获取卫星伪距测量误差方程的系数阵G阵和常数项具体包括：

GNSS一个历元上的n个观测值组成的误差方程为：

V＝Gδx-L

其中： $G=\left[\begin{array}{ccccccc}{l}^1\left(t\right)& m^1\left(t\right)& n^1\left(t\right)& -E& 0& 0& 0\\ l^2\left(t\right)& m^2\left(t\right)& n^2\left(t\right)& 0& -E& 0& 0\\ l^3\left(t\right)& m^3\left(t\right)& n^3\left(t\right)& 0& 0& -E& 0\\ l^4\left(t\right)& m^4\left(t\right)& n^4\left(t\right)& 0& 0& 0& -E\end{array}\right]$

$l\left(t\right)=\frac{X_s-X_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$ $m\left(t\right)=\frac{Y_s-Y_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$ $n\left(t\right)=\frac{Z_s-Z_{r0}}{|{\stackrel{\→}{x}}_s-{\stackrel{\→}{x}}_{r0}|}$

δx＝[δx_r1，δx_r2，δx_r3，CΔt₁，CΔt₂，CΔt₃，CΔt₄]^T。

5.根据权利要求1、2或3所述的方法，其特征在于，所述对P和p进行矩阵转换得到矩阵A具体包括：

ε₁P₁+ε₂P₂+...+ε_nP_n＝p进行矩阵形式的转换，则

$\left[\begin{array}{cccc}{P}_{11}{\mathrm{\ϵ}}_1& P_{21}{\mathrm{\ϵ}}_2& P_{1n}{\mathrm{\ϵ}}_n& p_1\\ P_{21}{\mathrm{\ϵ}}_1& P_{22}{\mathrm{\ϵ}}_2& P_{2n}{\mathrm{\ϵ}}_n& p_2\\ ...& ...& ...& ...\\ P_{m1}{\mathrm{\ϵ}}_1& P_{m2}{\mathrm{\ϵ}}_2& P_{\mathrm{mn}}{\mathrm{\ϵ}}_n& p_m\end{array}\right]=\left[\begin{array}{ccccc}{a}_1& a_2& ...& a_n& a_{n+1}\end{array}\right]=A.$

6.根据权利要求1、2或3所述的方法，其特征在于，所述组成相关距离矩阵B具体包括：

获取矩阵A各向量之间的距离，组成相关距离矩阵B：

$\left[\begin{array}{cccc}{b}_{11}& b_{21}& b_{m1}& b_{(n+1)1}\\ b_{12}& b_{22}& b_{m2}& b_{(n+1)2}\\ ...& ...& ...& ...\\ b_{1(n+1)}& b_{2(n+1)}& b_{m(n+1)}& b_{(n+1)(n+1)}\end{array}\right]=B$

其中 $b_{\mathrm{ij}}={\left[\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(a_{\mathrm{ik}}-a_{\mathrm{jk}})}^2\right]}^{\frac{1}{2}}.$

b_ij为a_i与a_j的相关距离，反映了各个向量之间的相关程度，采用b_ij数值指标，判别粗差存在于哪个观测值上。

7.根据权利要求1、2或3所述的方法，其特征在于，所述对矩阵B中各个向量距离进行聚类分析具体包括：

采用最小距离法进行粗差与偶然误差的分类，对计算的各类相关距离矩阵B进行聚类分析，得到最后一列向量与其他向量的亲疏关系，从而识别粗差。

Images