一种北斗非机动型广播星历故障实时监测方法

1.本发明属于导航技术领域，具体涉及一种北斗非机动型广播星历故障实时监测方法，特别是一种基于星历参数混合外推的北斗非机动型广播星历故障实时监测方法。


背景技术：

2.将卫星导航系统应用于飞机精密进近，可有效提高飞行安全和效率。飞机精密进近导航的完好性和连续性是航空用户最关注的性能需求，其中卫星广播星历故障是影响飞机精密进近的主要风险源之一。为保障与生命安全相关用户的定位完好性，对广播星历故障的监测是必不可少的。广播星历故障可分为机动性广播星历故障和非机动型广播星历故障。由于卫星机动发生的概率很小，所以非机动型广播星历故障是影响定位完好性的主要因素。
3.全球定位系统监测非机动型广播星历故障的方法是基于卫星运动的轨道重合期，采用零阶保持的方式将上一轨道周期的广播星历参数外推到当前时刻，计算卫星位置偏差，构建监测模型。由于每个广播星历参数对卫星位置的影响大小不同，对所有参数采用同一外推方式是不合理的。除此之外，北斗混合星座的构型及卫星运动的轨道重合周期和全球定位系统存在显著差异，所以传统非机动型广播星历故障监测方法无法直接应用于北斗系统。综上所述，设计一种适用于北斗多类型卫星的非机动型广播星历故障实时监测方法具有重要意义。


技术实现要素：

4.针对上述现有技术，本发明要解决的技术问题是提供一种基于星历参数混合外推的北斗非机动型广播星历故障实时监测方法，在保证较低计算量的同时得到更高精度的广播星历参数预测值，进一步获得更高精度的卫星位置偏差。
5.为解决上述技术问题，本发明的一种北斗非机动型广播星历故障实时监测方法，包括以下步骤：
6.步骤1：基于卫星的轨道运动重合周期及位置偏差选取广播星历参数外推时间；
7.步骤2：计算广播星历参数外推时间间隔的星历参数差值，分析其统计分布特性；
8.步骤3：卫星位置偏差对广播星历参数变化的灵敏度分析；
9.步骤4：根据灵敏度分析中广播星历参数对卫星位置的影响大小，对广播星历参数进行影响等级分类；
10.步骤5：对不同影响等级中的广播星历参数选择不同阶数的多项式外推；
11.步骤6：基于卫星位置偏差的统计分布特性，构建检验统计量和最小可检测误差；
12.步骤7：根据系统的连续性要求和完好性风险要求获取检验统计量阈值，通过阈值和检验统计量比较的方式，实现北斗非机动型广播星历故障实时监测，具体为：若计算得到的检验统计量小于检验统计量阈值，则满足系统的连续性要求，否则发生误警；在满足连续性要求的前提下，若计算的最小可检测误差小于其给定检测门限，则满足系统的完好性风
险要求，否则发生漏检。
13.本发明还包括：
14.1.步骤1中广播星历参数外推时间为1天。
15.2.步骤2中计算广播星历参数外推时间间隔的星历参数差值，分析其统计分布特性具体为：
16.设当前k时刻，相对于星历参考时间t
oe
，广播星历参数为p
k
，相邻日时刻的历史广播星历参数为p
k
‑1，计算广播星历参数偏差可表示为：
17.δp＝p
k
‑
p
k
‑118.其中，p代表由15个广播星历参数组成的矩阵；
19.通过计算相邻两天的广播星历参数差值得到其统计分布特性。
20.3.步骤3中卫星位置偏差对广播星历参数变化的灵敏度分析具体为：
21.通过使用北斗接口控制文件
‑
2013中定义的卫星位置非线性函数计算出卫星位置，卫星位置计算模型具体为：
[0022][0023]
其中，x、y、z为卫星三轴位置的计算函数，t为当前时刻，相对于t
oe
，p
i
表示不同广播星历参数，i＝1，2，
…
，15；
[0024]
卫星位置对广播星历参数变化的灵敏度系数通过卫星位置函数分别对每个广播星历参数的微分计算，具体为：
[0025][0026]
其中，δx(t)、δy(t)、δz(t)分别为卫星三轴位置误差，δp表示星历参数的变化量；
[0027]
简写为：
[0028]
δr(t)＝aδp
[0029]
其中，δr(t)表示t时刻卫星位置误差；
[0030]
则广播星历参数变化对卫星位置误差的影响为：
[0031][0032]
其中，σ(δr
i
)是卫星位置误差标准差，σ(δp
i
)是相邻日广播星历参数差值的标准差；
[0033]
根据上述公式可以得到每个广播星历参数变化分别对北斗地球静止轨道卫星、倾斜地球同步轨道卫星和中圆地球轨道卫星三种轨道卫星位置的影响大小。
[0034]
4.步骤4中根据灵敏度分析中广播星历参数对卫星位置的影响大小，对广播星历参数进行影响等级分类具体为：根据步骤3中的灵敏度分析结果，将对北斗地球静止轨道卫星、倾斜地球同步轨道卫星和中圆地球轨道卫星三种轨道卫星位置影响程度不同的广播星历参数分为三类：将对卫星位置影响误差超过1000m的广播星历参数归类为影响等级一，将影响卫星位置误差低于100m的广播星历参数归类为影响等级三，将其余广播星历参数归类为影响等级二。
[0035]
5.步骤5中对不同影响等级中的广播星历参数选择不同阶数的多项式外推具体为：
[0036]
根据步骤4中的广播星历参数影响等级分类，对影响等级三范围内的星历参数采用一阶的外推方式，具体为：
[0037]
p
k
＝p
k
‑1+(p
k
‑1‑
p
k
‑2)
[0038]
其中，p
k
表示当前广播星历的k时刻，p
k
‑
i
表示第i天前广播星历的k时刻；
[0039]
对影响等级一中内的参数采用三阶的外推方式，具体为：
[0040]
p
k
＝p
k
‑1+(p
k
‑1‑
p
k
‑2)+[(p
k
‑1‑
p
k
‑2)
‑
(p
k
‑2‑
p
k
‑3)]+{[(p
k
‑1‑
p
k
‑2)
‑
(p
k
‑2‑
p
k
‑3)]
‑
[(p
k
‑2‑
p
k
‑3)
‑
(p
k
‑3‑
p
k
‑4)]}
[0041]
对影响等级二中的参数采用二阶的外推方式，具体为：
[0042]
p
k
＝p
k
‑1+(p
k
‑1‑
p
k
‑2)+[(p
k
‑1‑
p
k
‑2)
‑
(p
k
‑2‑
p
k
‑3)]
[0043]
6.步骤6中基于卫星位置偏差的统计分布特性，构建检验统计量和最小可检测误差具体为：
[0044]
无星历异常条件下，计算当前时刻的卫星位置，采用步骤5中的外推方法将前一天的广播星历外推到当前时刻并计算预测卫星位置，得到卫星位置偏差δr；
[0045]
根据相邻日广播星历参数差值的统计分布特性，假设无星历异常情况下的位置偏差同样服从该统计分布；
[0046]
则计算检验统计量s满足：
[0047][0048]
其中，∑
δr
是卫星位置偏差的协方差矩阵；
[0049]
星历异常条件下的最小可检测误差为，
[0050][0051]
其中，λ是与系统完好性风险要求有关的参数。
[0052]
本发明的有益效果：本发明涉及广播星历参数自身的统计分布特性及卫星位置对广播星历参数变化的灵敏度分析，通过对影响卫星位置大小不同的广播星历参数采用高低阶多项式混合外推方式，构建故障监测模型，进而实现北斗非机动型广播星历故障实时监测。本发明基于广播星历参数变化的统计分布特性以及卫星位置对广播星历参数变化的灵敏度分析，对影响卫星位置大小不同的广播星历参数采用高低阶多项式混合外推方式，即对影响卫星位置较大的参数采用高阶外推方式，对影响卫星位置较小的参数采用低阶外推方式，通过当前广播星历参数计算卫星位置偏差和检验统计量，然后结合检验统计量阈值实现北斗非机动型广播星历故障实时监测。本发明基于星历参数混合的外推方式，可以在保证较低计算量的同时得到更高精度的广播星历参数预测值和卫星位置偏差，有助于构建非机动型广播星历故障的实时监测模型。本发明充分利用卫星位置对广播星历参数的灵敏度分析，对广播星历参数采用高低阶多项式混合外推的方式，即对影响卫星位置较大的参数采用高阶外推方式，对影响卫星位置较小的参数采用低阶外推方式，在保证较低计算量的同时得到更高精度的广播星历参数预测值，进一步获得更高精度的卫星位置偏差，有助于构建监测模型，可实现北斗非机动型广播星历故障实时监测。
附图说明
[0053]
图1是运用本发明的一种北斗非机动型广播星历故障实时监测方法的实施例示意图。
具体实施方式
[0054]
下面结合说明书附图和具体实施方式对本发明做进一步说明。
[0055]
本发明具体包括：步骤1，基于卫星的轨道运动重合周期及位置偏差选取广播星历参数外推时间；步骤2，计算广播星历参数外推时间间隔的星历参数差值，分析其统计分布特性；步骤3，卫星位置偏差对广播星历参数变化的灵敏度分析；步骤4，根据灵敏度分析中广播星历参数对卫星位置的影响大小，对广播星历参数进行影响等级分类；步骤5，对不同影响等级中的广播星历参数选择不同阶数的多项式外推；步骤6，基于卫星位置偏差的统计分布特性，构建检验统计量和最小可检测误差；步骤7，根据系统的连续性要求和完好性风险要求获取检验统计量阈值，通过阈值和检验统计量比较的方式，实现北斗非机动型广播星历故障实时监测。步骤2中基于相邻日广播星历参数差值的统计分布特性得到卫星位置偏差服从的统计分布。步骤3中分析卫星位置对广播星历参数变化的灵敏度，进而得到每个广播星历参数对卫星位置的影响大小。步骤4中基于步骤3中的灵敏度分析，对影响卫星位置不同的广播星历参数进行影响等级分类，有助于采用不同的方式外推广播星历参数，在
保证较小计算量的同时获取高精度预测值。步骤5中基于步骤4中的广播星历参数分类，对卫星位置影响较大的广播星历参数可采用高阶的外推方式，对卫星位置影响较小的广播星历参数可采用低阶的外推方式，将前一天的广播星历参数外推到当前时刻。采用这种星历参数混合外推的方式可以得到更高精度的广播星历参数预测值，进一步得到更高精度的卫星位置偏差，有助于构建星历故障监测模型。
[0056]
结合图1，本发明基于星历参数混合外推的北斗非机动型广播星历故障实时监测方法，具体步骤包括：
[0057]
步骤1，基于卫星的轨道运动重合周期及位置偏差选取广播星历参数外推时间。
[0058]
北斗地球静止轨道卫星、倾斜地球同步轨道卫星和中圆地球轨道卫星三种轨道卫星运动的轨道周期分别约为1天、1天、7天。由于7天间隔的卫星位置偏差过大，不利于构建广播星历故障监测模型，故选取1天间隔作为广播星历参数的外推时间。
[0059]
步骤2，计算广播星历参数外推时间间隔的星历参数差值，分析其统计分布特性。
[0060]
基于步骤1中选取的广播星历参数外推时间，计算广播星历参数差值。设当前k时刻(相对于星历参考时间t
oe
)广播星历参数为p
k
，相邻日时刻的历史广播星历参数为p
k
‑1，计算广播星历参数偏差可表示为，
[0061]
δp＝p
k
‑
p
k
‑1ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0062]
其中，p代表由15个广播星历参数组成的矩阵。
[0063]
通过计算相邻两天的广播星历参数差值可知其统计分布特性，该特性反映了卫星位置偏差的统计分布特性，在计算检验统计量时可假设卫星位置偏差与广播星历参数差值服从相同的统计分布特性。
[0064]
步骤3，卫星位置偏差对广播星历参数变化的灵敏度分析。
[0065]
通过使用北斗接口控制文件
‑
2013中定义的卫星位置非线性函数可以计算出卫星位置。卫星位置计算模型可简写为，
[0066][0067]
其中，x、y、z为卫星三轴位置的计算函数，t为当前时刻(相对于t
oe
)，p
i
表示不同广播星历参数，i＝1，2，
…
，15。
[0068]
卫星位置对广播星历参数变化的灵敏度系数可以通过卫星位置函数分别对每个广播星历参数的微分计算，公式如下，
[0069]
[0070]
其中，δx(t)、δy(t)、δz(t)分别为卫星三轴位置误差，δp表示星历参数的变化量。
[0071]
式(3)右边第一个矩阵为广播星历参数对卫星位置误差的灵敏度矩阵，记为矩阵a，矩阵a内元素是某一方向的卫星位置对某个星历参数的微分，将式(3)简写为，
[0072]
δr(t)＝aδp
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0073]
其中，δr(t)表示t时刻卫星位置误差。
[0074]
则广播星历参数变化对卫星位置误差的影响可以表示为，
[0075][0076]
其中，σ(δr
i
)是卫星位置误差标准差，σ(δp
i
)是相邻日广播星历参数差值的标准差。
[0077]
根据上述公式可以得到每个广播星历参数变化分别对北斗地球静止轨道卫星、倾斜地球同步轨道卫星和中圆地球轨道卫星三种轨道卫星位置的影响大小。
[0078]
步骤4，根据灵敏度分析中广播星历参数对卫星位置的影响大小对广播星历参数进行影响等级分类。
[0079]
基于步骤3中的灵敏度分析结果，将对北斗地球静止轨道卫星、倾斜地球同步轨道卫星和中圆地球轨道卫星三种轨道卫星位置影响程度不同的广播星历参数分为三类。将对卫星位置影响误差超过1000m的广播星历参数归类为影响等级一，将影响卫星位置误差低于100m的广播星历参数归类为影响等级三，将其余广播星历参数归类为影响等级二。
[0080]
步骤5，对不同影响等级中的广播星历参数选择不同阶数的多项式外推。
[0081]
基于步骤4中的广播星历参数影响等级分类，对影响等级一范围内的广播星历参数可采高阶的外推方式，对影响等级三范围内的广播星历参数可采用低阶的外推方式。
[0082]
对影响等级三范围内的星历参数采用一阶的外推方式，公式如下，
[0083]
p
k
＝p
k
‑1+(p
k
‑1‑
p
k
‑2)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0084]
其中，p
k
表示当前广播星历的k时刻，p
k
‑
i
表示第i天前广播星历的k时刻。
[0085]
对影响等级一中内的参数采用三阶的外推方式，
[0086][0087]
对影响等级二中的参数采用二阶的外推方式，
[0088]
p
k
＝p
k
‑1+(p
k
‑1‑
p
k
‑2)+[(p
k
‑1‑
p
k
‑2)
‑
(p
k
‑2‑
p
k
‑3)]
ꢀꢀꢀꢀꢀꢀ
(8)
[0089]
采用这种高低阶多项式混合外推先前的广播星历参数，在降低算法复杂度的同时，还可获得更高精度的广播星历参数预测值，进一步获得更高精度的卫星位置偏差，有助于构建监测模型。
[0090]
步骤6，基于卫星位置偏差的统计分布特性，构建检验统计量和最小可检测误差。
[0091]
无星历异常条件下，计算当前时刻的卫星位置，采用步骤5中的外推方法将前一天的广播星历外推到当前时刻并计算预测卫星位置，得到卫星位置偏差δr。
[0092]
根据相邻日广播星历参数差值的统计分布特性，假设无星历异常情况下的位置偏差同样服从该统计分布。
[0093]
基于上述卫星位置偏差及其分布特性，计算检验统计量s，
[0094][0095]
其中，∑
δr
是卫星位置偏差的协方差矩阵。
[0096]
星历异常条件下的最小可检测误差为，
[0097][0098]
其中，λ是与系统完好性风险要求有关的参数。
[0099]
步骤7，根据系统的连续性要求和完好性风险要求获取检验统计量阈值，通过阈值和检验统计量比较的方式，实现北斗非机动型广播星历故障实时监测。
[0100]
若计算得到的检验统计量小于检验统计量阈值，则满足系统的连续性要求，否则发生误警。在满足连续性要求的前提下，若计算的最小可检测误差小于其检测门限，则满足系统的完好性风险要求，否则发生漏检。基于以上步骤完成对北斗非机动型广播星历故障的监测。
[0101]
本发明还可有其他多种实施例，在不偏离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明做出各种相应的调整，但这些相应的调整都应属于本发明所附的权利要求的保护范围。