CN102073054B

CN102073054B - 地基增强系统的完好性监测方法

Info

Publication number: CN102073054B
Application number: CN2010105274071A
Authority: CN
Inventors: 张军; 朱衍波; 刘霞; 陈海龙
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2010-10-26
Filing date: 2010-10-26
Publication date: 2012-11-14
Anticipated expiration: 2030-10-26
Also published as: CN102073054A

Abstract

本发明提供一种地基增强系统的完好性监测方法，其中包括：对接收到的卫星导航信号的载噪比进行实时归一化处理，得到归一化载噪比；对所述归一化的载噪比进行实时的时间累积，得到当前时刻载噪比累积和；根据所述当前时刻载噪比累积和通过马尔可夫链计算得到监测门限阈值的当前值；将所述当前时刻载噪比累积和与所述当前值进行比较，若超过所述当前值，则发出告警信息。本发明在GBAS载噪比完好性监测的数据处理过程中，采用了一种误差偏差累积的数学统计方法对导航卫星信号的信噪比进行了实时监测，简化了计算量，并且也实现了阈值的实时更新，从而提高了定位结果的精度。

Description

地基增强系统的完好性监测方法

技术领域

本发明涉及一种地基增强系统的完好性监测方法，属于卫星导航技术领域。

背景技术

全球导航卫星系统(Global Navigation Satellite System，简称：GNSS)的定位原理是：用户根据接收到的导航卫星发出的卫星导航信号计算用户接收机到达卫星的距离和卫星的实时位置，当用户同时计算出到达四颗卫星的距离时，就可以解算出该用户自己的实时三维位置。然而，用户测得的卫星到用户接收机的距离中实际上还包含卫星的星历、星钟、电离层和对流层延迟以及接收机热噪声和多径效应等引起的误差，因此测得的这个距离实际上为伪距。例如，对于全球定位系统(Global Positioning System，简称：GPS)用户来说，直接通过伪距计算的定位结果大约有12.5米的误差。因此，为了获得更高精度的定位结果，出现了一些辅助系统，如地基增强系统(GroundBased Augmentation System，简称：GBAS)。

GBAS的基本原理是利用定位误差空间和时间的相关性(即认为机载用户和地面监测站在100km内的定位误差是相同的，即公共误差)，机载用户可以通过接收地面站播发的伪距差分校正值和系统完好性信息来提高导航的精度和完好性。但采用现有GBAS进行完好性监测时通过需要根据当前的地形、接收机的位置和卫星的位置分别计算不同的阈值，并且需要进行大量的数据统计，计算非常复杂；并且计算出的阈值不能实时更新，因此不能及时反映最新的环境变化，从而影响了定位结果的精度。

发明内容

本发明提供一种地基增强系统的完好性监测方法，用以实现阈值的实时更新，以提高定位结果的精度。

本发明提供一种地基增强系统的完好性监测方法，其中包括：

对接收到的卫星导航信号的载噪比进行实时归一化处理，得到归一化载噪比；

对所述归一化的载噪比进行实时的时间累积，得到当前时刻载噪比累积和；

根据所述当前时刻载噪比累积和通过马尔可夫链计算得到监测门限阈值的当前值；

将所述当前时刻载噪比累积和与所述当前值进行比较，若超过所述当前值，则发出告警信息。

本发明在GBAS载噪比完好性监测的数据处理过程中，采用了一种误差偏差累积(CUSUM)的数学统计方法对导航卫星信号的信噪比进行了实时监测，简化了计算量，并且也实现了阈值的实时更新，从而提高了定位结果的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明所述地基增强系统的完好性监测方法实施例的流程图；

图2为图1中步骤300的可选计算方法的流程图；

图3为图2中步骤310的可选计算方法的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明所述地基增强系统的完好性监测方法实施例的流程图，如图所示，该方法包括如下步骤：

步骤100，对接收到的卫星导航信号的载噪比进行实时归一化处理，得到归一化载噪比。

具体地，可以根据公式

进行所述归一化处理，得到所述归一化载噪比X_i，其中：i＝1，2，…n，n表示实时监测的次数，C/N₀(i)表示第i次监测时得到的载噪比，表示第i次监测时得到的载噪比的均值，

表示第i次监测时得到的载噪比的方差，X_i服从N(0，1)的标准正态分布。其中，所述卫星导航信号的载噪比作在本实施例中作为监测量，有关监测量的详细介绍将在后续内容进行说明。

步骤200，对所述归一化的载噪比进行实时的时间累积，得到当前时刻载噪比累积和。

其中，所述当前时刻载噪比累积和是指对前n个时刻的归一化载噪比数据进行累加后得到的值。具体地，可以根据公式

对所述归一化的载噪比进行实时的时间累积，得到当前时刻载噪比累积和S_n。

步骤300，根据所述当前时刻载噪比累积和通过马尔可夫链计算得到监测门限阈值的当前值。

具体计算过程将在后续内容进行详细说明。

步骤400，将所述当前时刻载噪比累积和与所述当前值进行比较，若超过所述当前值则发出告警信息。

具体地，若

则发出告警信息，否则表明GNSS系统当前处于受控状态，而无需发出告警信息。

本实施例所述方法在GBAS载噪比完好性监测的数据处理过程中，采用了一种误差偏差累积(CUSUM)的数学统计方法对导航卫星信号的信噪比进行了实时监测，简化了计算量，并且也实现了阈值的实时更新，从而提高了定位结果的精度。

以下对步骤100中提到的监测量进行说明：

对被监测对象的卫星导航信号进行连续采样得到一组的载噪比x₁，x₂，...，x_v，...作为监测量，下标表示采样时刻，在v时刻之前被监测对象处于受控状态，即未发生故障，而在v时刻以后则处于失控状态，即发生了故障。本实施例所述方法就是要及时有效地找到开始出现问题的v时刻。

在受控状态下，监测量x_i，i＝1，2，…，v-1服从正态分布且方差相同(独立同分布)，即：x_i～N(u₀，σ²)，其中，u₀表示x_i的均值，σ²表示x_i的方差，其概率密度函数为：

f (x_{i}) = \frac{1}{\sqrt{2 π} σ} e^{\frac{- {(x_{i} - μ_{0})}^{2}}{2 σ^{2}}} - - - (1)

在失控状态下，监测量x_i，i＝v，…服从正态分布且方差相同(独立同分布)即：x_i～N(u₁，σ²)，其中，均值u₁与受控状态下的x_i的均值u₀不同，但方差相同。其概率密度函数为，

f (x_{i}) = \frac{1}{\sqrt{2 π} σ} e^{\frac{- {(x_{i} - μ_{1})}^{2}}{2 σ^{2}}} - - - (2)

对随机变量x_i进行归一化处理，

并代入式(1)、(2)则两种情况下的随机变量分别服从正态分布X_i～N(0，1)，i＝1，2，…，`’，X_i～N(u₁-u₀，1)，i＝v，…。

进行归一化处理后，可以归结为随机变量X_i的变点v的假设检验问题：设一系列随机变量X₁，X₂，...

H₀：X_i～N(u₀，σ²)，i＝1，2…

H₁：存在v(v＜∞)，X_i～N(u₀，σ²)，i＝1，2，…，v

X_i～N(u₁，σ²)，i＝v+1；

归一化处理后，即令：

则变为

H₀：Y_i～N(0，1)，i＝1，2…

H₁：存在v(v＜∞)，Y_i～N(0，1)，i＝1，2，…，v；

Y_i～N(u₁-μ₀，1)，i＝v+1；

用极大似然估计进行参数估计，极大似然比为：

L_{n, v} = \frac{Π_{i = 1}^{v} \frac{1}{\sqrt{2 π}} e^{\frac{- {y_{i}}^{2}}{2}} Π_{i = v + 1}^{n} \frac{1}{\sqrt{2 π}} e^{\frac{- {(y_{i} - δ)}^{2}}{2}}}{Π_{i = 1}^{n} \frac{1}{\sqrt{2 π}} e^{\frac{- {y_{i}}^{2}}{2}}} = Π_{i = v + 1}^{n} e^{y_{i} δ - \frac{δ^{2}}{2}}

其中，δ＝u₁-μ₀，是故障的平均跳变幅度，两边同时取对数得到：

\ln L_{n, v} = δ Σ_{i = v + 1}^{n} (y_{i} - \frac{δ}{2})

则对数似然比的统计量为：

Z_{n} = \max_{1 \leq v \leq n} {δ Σ_{i = v + 1}^{n} (y_{i} - \frac{δ}{2})} = \max_{1 \leq v \leq n} {Σ_{i = v + 1}^{n} (y_{i} - \frac{δ}{2})}

这就是最简单的CUSUM监测统计量，若事先选定了一个门限h，当Z_i≤h，i＝1，2，3，…，n说明到n时刻为止，监测量仍然处于受控状态，当Z_n＞h时，则表明，n时刻，监测量处于失控状态。

另外，为了改善CUSUM的监测性能也可以将Z_n＝max{0，Z_n-1+y_n}作为监测统计量，具体说明如下：

令，

则：其中，

由于，

因此：

Z_{n} - Z_{n - 1} = {\tilde{y}}_{n} - \min {0, {\tilde{S}}_{n} - \min_{1 \leq v \leq n - 1} {\tilde{S}}_{v}} = \max {{\tilde{y}}_{n}, {\tilde{y}}_{n} - {\tilde{S}}_{n} + \min_{1 \leq v \leq n - 1} {\tilde{S}}_{v}}

= \max {{\tilde{y}}_{n}, \min_{1 \leq v \leq n - 1} {\tilde{S}}_{v} - {\tilde{S}}_{n - 1}} = \max {{\tilde{y}}_{n}, - Z_{n - 1}}

所以，

Z_{n} = \max {0, Z_{n - 1} + {\tilde{y}}_{n}},

Z_{n} = \max {0, Z_{n - 1} + y_{n} - \frac{δ}{2}}

令，k＝δ/2，则Z_n＝max{0，Z_n-1+y_n-k}

这样，CUSUM由两个参数所决定(k，h)，k表示故障空间与正常空间的距离。

下面，分析一下CUSUM的不足：

Z_n可改写为

Z_{n} = Σ_{i = 1}^{n} y_{i} - nk

两边同时除以n，得

则Z_n＞h变为

其中，k为一个常数，当n很大时，在受控状态下，

趋近于零，所以对于小于k的故障监测不出来。

其实，从上式可以看出，其判决器的实际物理意义是：看监测量的均值更接近于受控状态的均值还是失控状态的均值，而做出其是处于哪个状态的判决，例如当其

大于三倍的标准差时，就认为其处于失控状态了。所以，可以看出，其监测门限与监测量的标准差有关。

又因为，Y₁，Y₂，…是独立同分布的Y_i～N(u₁-μ₀，1)，

所以例如当

超出其一倍的标准差

就不再处于受控状态了，这样可以根据监测量的数量动态更新门限阈值，就可以对幅度较小、较大的故障都可以有效地监测。

改进后的监测统计量为Z_n＝max{0，Z_n-1+y_n}

当

报警，m初值为零，当Z_n＝0，时，更新m，m＝n；其门限根据监测量的数量实时更新。

以下说明上述步骤300的计算过程，如图2所示，可以包括如下步骤：

步骤310，通过马尔可夫链计算转移概率

p_{n} (h) = \Pr (S_{n} < 0 | 0 \leq S_{i} \leq \sqrt{n} h, i = 1,2, . . ., n - 1) = \Pr (S_{n} < 0 | 0 \leq S_{n - 1} \leq \sqrt{n} h)

和

q_{n} (h) = \Pr (S_{n} > h | 0 \leq S_{i} \leq \sqrt{n} h, i = 1,2, . . ., n - 1) = \Pr (S_{n} > h | 0 \leq S_{n - 1} \leq \sqrt{n} h)

其中：P_n(h)表示所述当前时刻载噪比累积和从受控状态转移到大于所述前时刻监测门限阈值h的状态的转移概率，q_n(h)表示当前时刻载噪比累积和从受控状态转移到小于零状态的转移概率；

步骤320，估计监测门限阈值h的初值。

例如，h＝1。

步骤330，根据公式一

在所述初值的条件下计算受控状态下的平均运行时间ARL₀，其中：

p_{0} (h) = Σ_{i = 1}^{\infty} {p_{i} (h) Π_{j = 1}^{i - 1} (1 - p_{j} (h) - q_{j} (h))}

q_{0} (h) = Σ_{i = 1}^{\infty} {i (p_{i} (h) + q_{i} (h)) Π_{j = 1}^{i - 1} (1 - p_{j} (h) - q_{j} (h))}

p_{n} (h) = \Pr (S_{n} < 0 | 0 \leq S_{i} \leq \sqrt{n} h, i = 1,2, . . ., n - 1)

q_{n} (h) = \Pr (S_{n} > h | 0 \leq S_{i} \leq \sqrt{n} h, i = 1,2, . . ., n - 1)

n＝1，2，3，…，S_n表示其第n个时刻所处的状态；

由Z_n的定义知，若

则报警，其中，m＝max{j：j＜n，S_j＝0}，在一次检验中，Z_n是带有两个吸收壁的随机游动，它要么被0吸收，要么被h吸收，p₀(h)、q₀(h)分别表示被0和h吸收的概率；

ARL₀是指监测系统从被监测对象处于受控状态开始运行，直到处于失控状态的平均运行时间，理想情况下，ARL₀是无穷大，即永远都没有故障；ARL₁是指监测系统从被监测对象处于失控状态开始，直到处于失控状态的平均运行时间，因此，ARL₁越小越好。

其中，在受控状态下，GNSS系统的完好性风险不应该超过2*10^-7，而上述ARL₀是完好性风险的倒数，因此，在受控状态下，要求ARL₀＞5×10⁶。在上算上述监测门限阈值h时，可以首先给出一个h的初值，然后用本步骤进行反复迭代更新h，直到该值满足ARL₀的上述要求为止。

步骤340，将所述受控状态下的平均运行时间ARL₀与误警概率α_f的倒数进行比较，当满足不等式

时，将所述初值作为所述当前值。

步骤350，当不满足所述不等式时，则增加所述初值，并根据所述公式一在增加后的初值条件下计算所述受控状态下的平均运行时间ARL₀，直到满足所述不等式时，将能够满足该不等式的初值作为所述当前值。

以下详细介绍上述步骤310中，如图3所示，通过马尔可夫链计算转移概率的过程如下：

步骤311，将所述当前时刻载噪比累积和S_n离散化为如下三种状态：

状态0 S₀(d)＝0

状态i S_i(d)∈[(i-1)Δ，iΔ]i＝1，2，...，M

状态M+1S_M+1＞h

其中：M表示离散化间距，

步骤312，根据如下公式计算马尔可夫状态转移概率矩阵Pr：

\Pr (a < S_{n} < b | c < S_{n - 1} < d)

= \frac{[F (b - c) + 4 F (b - m) + F (b - f)] - [F (a - c) + 4 F (a - m) + F (a - f)}{6}

其中：F(·)表示X_n的分布函数，m表示[c，f]的区间中间点，即m＝(c+d)/2。所述马尔可夫状态转移概率矩阵Pr中的元素表示为：

P_i，j＝Pr(S_n+1(h)∈Statej|S_n(h)∈State i)i，j＝0，1，...，M+1

马尔可夫状态转移概率矩阵Pr如下式所示。

步骤313，根据马尔可夫状态转移概率矩阵Pr计算得到所述马尔可夫链计算转移概率p_n(h)和q_n(h)。

其中，p_n(h)＝P_i0，S_n-1∈State i；q_n(h)＝P_iM+1，S_n-1∈State M+1。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims

1.一种地基增强系统的完好性监测方法，其特征在于，包括：

所述根据所述当前时刻载噪比累积和通过马尔可夫链计算得到监测门限阈值的当前值包括：

通过马尔可夫链计算转移概率p_n(h)＝Pr(S_n＞h|0≤S_i≤h，i＝1，2，...，n-1)和q_n(h)＝Pr(S_n＜0|0≤S_i≤h，i＝1，2，...，n-1)，其中：n表示实时监测的次数，S_n表示当前时刻载噪比累积和，P_n(h)表示所述当前时刻载噪比累积和从受控状态转移到大于监测门限阈值h的状态的转移概率，q_n(h)表示当前时刻载噪比累积和从受控状态转移到小于零状态的转移概率；

估计监测门限阈值h的初值；

根据公式一