CN108535746B - 一种探测gnss卫星轨道机动的方法 - Google Patents

一种探测gnss卫星轨道机动的方法 Download PDF

Info

Publication number
CN108535746B
CN108535746B CN201810161688.XA CN201810161688A CN108535746B CN 108535746 B CN108535746 B CN 108535746B CN 201810161688 A CN201810161688 A CN 201810161688A CN 108535746 B CN108535746 B CN 108535746B
Authority
CN
China
Prior art keywords
satellite
calculation
formula
calculating
observation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810161688.XA
Other languages
English (en)
Other versions
CN108535746A (zh
Inventor
袁运斌
叶飞
张宝成
潭冰峰
欧吉坤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Geodesy and Geophysics of CAS
Original Assignee
Institute of Geodesy and Geophysics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Geodesy and Geophysics of CAS filed Critical Institute of Geodesy and Geophysics of CAS
Priority to CN201810161688.XA priority Critical patent/CN108535746B/zh
Publication of CN108535746A publication Critical patent/CN108535746A/zh
Application granted granted Critical
Publication of CN108535746B publication Critical patent/CN108535746B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/25—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS
    • G01S19/258—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS relating to the satellite constellation, e.g. almanac, ephemeris data, lists of satellites in view
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Signal Processing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

一种探测GNSS卫星轨道机动的方法，该方法包括以下步骤：先通过文件或者网络获取全球测站GNSS观测值和广播星历；再将全球GNSS观测网分成多个子网，然后为每个子网选择一个装备有原子钟的测站作为基准站，其他测站作为参考站，再为子网内每个测站构建线性化后的消电离层组合平滑伪距观测方程，然后通过两个线程进行卫星机动判断。本设计不仅可以全天候对所有卫星进行可靠的单历元机动探测，而且能为将来实时卫星机动提供一个可行的思路。

Description

一种探测GNSS卫星轨道机动的方法
技术领域
本发明属于“测绘科学与技术”学科中的“大地测量”技术领域，尤其涉及一种探测GNSS卫星轨道机动的方法，主要适用于提高机动卫星定轨精度以及机动卫星的有效利用等方面。
背景技术
卫星定轨作为全球导航卫星系统(GNSS)的一项关键技术，被广泛研究以消除或削弱影响其准确性和稳定性的因素。此外，全球导航卫星系统用户也期望高精度和可靠的轨道产品。但是，无论是进行后处理的精确轨道确定，还是使用预测的轨道，确定卫星机动都很重要。在机动的情况下，需要在多天解的精密轨道确定中为机动卫星设置单独的子弧。在实时应用中，我们通常使用预测轨道，称为广播星历和超快速产品，而使用未被发现的机动卫星是有风险的。因此，实时检测卫星机动以保持导航系统服务性能至关重要。
卫星激光测距可以监测地球静止轨道卫星机动。然而，需要一个相应的GNSS卫星激光测量系统。类似地，使用双向传输测距站点观测到的传输测距数据可用于检测卫星机动，然而，对于中地球轨道和倾斜地球同步轨道卫星，必须要建立大量的双向传输测距站点，以实现全天的监测。此外，还可以使用雷达数据来检测卫星机动，雷达数据可以为空间目标提供运动状态监测。不幸的是，对于GNSS卫星，也需要额外的设备。使用广播星历，可以直接简单地在后期模式下检测卫星机动，但该方法不能满足实时性要求，时间分辨率也不理想。
全球GNSS观测数据和广播星历可以由IGS、MGEX和iGMAS免费提供，这些数据可以充分用于全天候监测卫星机动。此外，IGS在最近几年已经开发了从地面站实时传输数据的能力。因而使用这些实时数据实时检测全天候卫星机动的条件正逐渐可用。CODE(http：//www.aiub.unibe.ch/download/BSWUSER52/GEN/SAT_yyyy.CRX，其中yyyy是四位数年份)和美国安全国土部的导航中心(https://www.navcen.uscg.gov/？pageName＝selectNanuByNumber)虽然公布这类卫星机动的信息，但是该信息不是实时提供的。
如何利用GNSS测站跟踪网探测GNSS卫星机动一直是GNSS领域的研究热点，其中斯坦福大学和欧洲定轨中心团队对此进行了深入研究，也得到了广泛应用。但是斯坦福大学所提出的方法只适用于区域观测网络。此外，当消除接收机时钟钟差时，该方法易受电离层延迟，卫星机动和异常值的影响。而欧洲定轨中心提出的算法仅适用于后处理模式。
因此，现有的方法并不能完全满足实时性、全天候和高可靠性的要求。为了给未来实现实时可靠的机动探测提供另一个可行的思路，需要基于GNSS测站跟踪网来研究一种可靠的单历元解算方法，从而实现卫星机动探测。
发明内容
本发明的目的是克服现有技术中存在的只适用于区域观测网络、误判卫星机动，且当消除接收机时钟钟差时，易受电离层延迟、卫星机动和异常值的影响，或者只适用于事后处理模式的缺陷与问题，提供一种基于广播星历和全球参考网实现全天候探测卫星机动，且能为将来实时卫星机动提供一个可行思路的探测GNSS卫星轨道机动的方法。
为实现以上目的，本发明的技术解决方案是：一种探测GNSS卫星轨道机动的方法，该方法包括以下步骤：
第一步：通过文件或者网络获取全球测站GNSS观测值和广播星历；
第二步：先将全球GNSS观测网分成多个子网，再为每个子网选择一个装备有原子钟的测站作为基准站，其他测站作为参考站，然后将每个子网分线程执行第三步和第四步；
第三步：为子网内每个测站构建线性化后的消电离层组合平滑伪距观测方程：
Figure GDA0002439491930000021
式(1)中，
Figure GDA0002439491930000022
表示伪距残差，
Figure GDA0002439491930000023
表示视线方向的单位向量，Δs表示广播星历在径向dr、切向da、法向dc三个方向的误差，dtr表示接收机钟差，dts表示卫星钟差，c表示光速，ε表示噪声；
第四步、卫星机动探测：
A、根据广播星历提供的轨道和钟差简化式(1)得到式(2)，并根据式(2)计算基准站钟差；
Figure GDA0002439491930000031
式(2)中，m表示基准站接收机，
Figure GDA0002439491930000032
表示方差；则基准站钟差计算步骤如下：
A1、计算伪距残差
Figure GDA0002439491930000033
的中位数：
Figure GDA0002439491930000034
式(3)中，s＝1，2，...，N，N表示基准站观测到的卫星数，median{.}是计算序列
Figure GDA0002439491930000035
的中位数；
A2、计算伪距残差
Figure GDA0002439491930000036
的初始残差：
Figure GDA0002439491930000037
A3、计算伪距残差
Figure GDA0002439491930000038
的抗差方差因子：
Figure GDA0002439491930000039
A4、计算伪距残差
Figure GDA00024394919300000310
的抗差等价权：
Figure GDA00024394919300000311
式(6)中，c0＝2，3，4或者c0＝5；
A5、计算基准站钟差：
Figure GDA00024394919300000312
B、当已知基准站接收机钟差后，得到如下式(8)，然后消除其他参考站的接收机钟差；
Figure GDA0002439491930000041
式(8)中，
Figure GDA0002439491930000042
表示基准站m与其他参考站R之间的伪距残差之差，R＝1，2，...J，J表示子网内其他参考站的数量；则消除其他参考站钟差步骤如下：
B1、计算
Figure GDA0002439491930000043
的中位数：
Figure GDA0002439491930000044
B2、计算
Figure GDA0002439491930000045
的初始残差：
Figure GDA0002439491930000046
B3、计算
Figure GDA0002439491930000047
的抗差方差因子：
Figure GDA0002439491930000048
B4、计算
Figure GDA0002439491930000049
的抗差等价权：
Figure GDA00024394919300000410
B5、计算参考站接收机钟差：
Figure GDA00024394919300000411
B6、根据计算的基准站钟差和其他参考站接收机钟差得到如下方程：
Figure GDA0002439491930000051
Figure GDA0002439491930000052
B7、计算
Figure GDA0002439491930000053
和
Figure GDA0002439491930000054
C、假定s1为将要探测的卫星号，分别通过两个线程进行卫星机动判断：
C1、根据线程一得到卫星机动的第一个判断条件，线程一包括以下步骤：
C11、根据下式获取探测的变量
Figure GDA0002439491930000055
Figure GDA0002439491930000056
C12、计算
Figure GDA0002439491930000057
的中位数：
Figure GDA0002439491930000058
C13、计算
Figure GDA0002439491930000059
的初始残差：
Figure GDA00024394919300000510
C14、计算
Figure GDA00024394919300000511
的初始抗差方差因子：
Figure GDA00024394919300000512
C15、计算
Figure GDA00024394919300000513
的抗差因子和等价权：
Figure GDA00024394919300000514
Figure GDA0002439491930000061
C16、计算
Figure GDA0002439491930000062
的抗差方差因子：
Figure GDA0002439491930000063
C17、第一个判断条件为：
Figure GDA0002439491930000064
C2、根据线程二得到卫星机动的第二个判断条件，线程二包括以下步骤：
C21、根据广播星历的精度设置阈值Threshold1等于10米；
C22、根据式(15)计算Δs1(dr1 da1 dc1)；
C23、计算三维误差
Figure GDA0002439491930000065
Figure GDA0002439491930000066
C24、若卫星在连续十个历元没有发生机动或者异常情况，则得到一个序列如式(25)，而阈值Threshold2的计算如式(26)：
Figure GDA0002439491930000067
Threshold2＝c0(median{Ds1}/0.6745) (26)；
C25、第二个判断条件为：
Figure GDA0002439491930000068
C3、卫星机动探测的综合评估：
若在c0个连续历元同时满足方程(23.a)和(27.a)，则认为卫星s1在历元t发生了机动；
若在历元t只有(23.a)或(27.a)满足，则认为卫星s1在历元t不可用；反之，则认为卫星s1在历元t可用。
步骤C3中，评估一整天的卫星机动时，若子网中有任何卫星发生机动，则认为该卫星发生了机动。
与现有技术相比，本发明的有益效果为：
本发明一种探测GNSS卫星轨道机动的方法中充分利用GNSS全球观测网数据和广播星历，采用抗差方法，提供了一种更可靠的基于广播星历和全球参考网实现全天候探测卫星机动的单历元解算方法。因此，本发明可以全天候对所有卫星进行可靠的单历元机动探测，并且为将来实时卫星机动提供了一个可行的思路。
附图说明
图1是本发明的实施方案流程图。
图2是本发明实施例中探测卫星机动的测站分布图。
图3是本发明实施例中在2017/8/25用本发明三步法计算出的G09同步伪距残差。
图4是本发明实施例中在2017/8/25用本发明三步法计算出的G09三维误差。
图5是本发明实施例中在2017/8/25用传统方法计算出的G09同步伪距残差。
图6是本发明实施例中在2017/9/12用本发明三步法计算出的G07同步伪距残差。
图7是本发明实施例中在2017/9/12用本发明三步法计算出的G07三维误差。
图8是本发明实施例中在2017/9/12用传统方法计算出的G07同步伪距残差。
图9是本发明实施例中精密定轨所用的测站分布图。
图10是本发明实施例中在2017/8/25用本发明三步法得到的机动信息在精密定轨后与CODE轨道比较的RMS。
图11是本发明实施例中在2017/8/25用本发明三步法得到的机动信息在精密定轨后与CODE轨道比较的STD。
图12是本发明实施例中在2017/8/25用本发明传统法得到的机动信息在精密定轨后与CODE轨道在三个方向的比较结果。
图13是本发明实施例中在2017/9/12用本发明三步法得到的机动信息在精密定轨后与CODE轨道比较的RMS。
图14是本发明实施例中在2017/9/12用本发明三步法得到的机动信息在精密定轨后与CODE轨道比较的STD。
图2中，1、2、3、4、5、6分别表示一个子网。
图10、图11、图13、图14中，每三个条状为一组，每组中条状从左到右依次表示dx、dy、dz。
具体实施方式
以下结合附图说明和具体实施方式对本发明作进一步详细的说明。
一种探测GNSS卫星轨道机动的方法，该方法包括以下步骤：
第一步：通过文件或者网络获取全球测站GNSS观测值和广播星历；
第二步：先将全球GNSS观测网分成多个子网，再为每个子网选择一个装备有原子钟的测站作为基准站，其他测站作为参考站，然后将每个子网分线程执行第三步和第四步；
第三步：为子网内每个测站构建线性化后的消电离层组合平滑伪距观测方程：
Figure GDA0002439491930000081
式(1)中，
Figure GDA0002439491930000082
表示伪距残差，
Figure GDA0002439491930000083
表示视线方向的单位向量，Δs表示广播星历在径向dr、切向da、法向dc三个方向的误差，dtr表示接收机钟差，dts表示卫星钟差，c表示光速，ε表示噪声；
第四步、卫星机动探测：
A、根据广播星历提供的轨道和钟差简化式(1)得到式(2)，并根据式(2)计算基准站钟差；
Figure GDA0002439491930000084
式(2)中，m表示基准站接收机，
Figure GDA0002439491930000085
表示方差；则基准站钟差计算步骤如下：
A1、计算伪距残差
Figure GDA0002439491930000086
的中位数：
Figure GDA0002439491930000091
式(3)中，s＝1，2，...，N，N表示基准站观测到的卫星数，median{.}是计算序列
Figure GDA0002439491930000092
的中位数；
A2、计算伪距残差
Figure GDA0002439491930000093
的初始残差：
Figure GDA0002439491930000094
A3、计算伪距残差
Figure GDA0002439491930000095
的抗差方差因子：
σ0＝median{|vs}/0.6745 (5)；
A4、计算伪距残差
Figure GDA0002439491930000096
的抗差等价权：
Figure GDA0002439491930000097
式(6)中，c0＝2，3，4或者c0＝5；
A5、计算基准站钟差：
Figure GDA0002439491930000098
B、当已知基准站接收机钟差后，得到如下式(8)，然后消除其他参考站的接收机钟差；
Figure GDA0002439491930000099
式(8)中，
Figure GDA00024394919300000910
表示基准站m与其他参考站R之间的伪距残差之差，R＝1，2，...J，J表示子网内其他参考站的数量；则消除其他参考站钟差步骤如下：
B1、计算
Figure GDA0002439491930000101
的中位数：
Figure GDA0002439491930000102
B2、计算
Figure GDA0002439491930000103
的初始残差：
Figure GDA0002439491930000104
B3、计算
Figure GDA0002439491930000105
的抗差方差因子：
Figure GDA0002439491930000106
B4、计算
Figure GDA0002439491930000107
的抗差等价权：
Figure GDA0002439491930000108
B5、计算参考站接收机钟差：
Figure GDA0002439491930000109
B6、根据计算的基准站钟差和其他参考站接收机钟差得到如下方程：
Figure GDA00024394919300001010
Figure GDA00024394919300001011
B7、计算
Figure GDA00024394919300001012
和
Figure GDA00024394919300001013
C、假定s1为将要探测的卫星号，分别通过两个线程进行卫星机动判断：
C1、根据线程一得到卫星机动的第一个判断条件，线程一包括以下步骤：
C11、根据下式获取探测的变量
Figure GDA0002439491930000111
Figure GDA0002439491930000112
C12、计算
Figure GDA0002439491930000113
的中位数：
Figure GDA0002439491930000114
C13、计算
Figure GDA0002439491930000115
的初始残差：
Figure GDA0002439491930000116
C14、计算
Figure GDA0002439491930000117
的初始抗差方差因子：
Figure GDA0002439491930000118
C15、计算
Figure GDA0002439491930000119
的抗差因子和等价权：
Figure GDA00024394919300001110
Figure GDA00024394919300001111
C16、计算
Figure GDA00024394919300001112
的抗差方差因子：
Figure GDA00024394919300001113
C17、第一个判断条件为：
Figure GDA0002439491930000121
C2、根据线程二得到卫星机动的第二个判断条件，线程二包括以下步骤：
C21、根据广播星历的精度设置阈值Threshold1等于10米；
C22、根据式(15)计算Δs1(dr1 da1 dc1)；
C23、计算三维误差
Figure GDA0002439491930000122
Figure GDA0002439491930000123
C24、若卫星在连续十个历元没有发生机动或者异常情况，则得到一个序列如式(25)，而阈值Threshold2的计算如式(26)：
Figure GDA0002439491930000124
Threshold2＝c0(median{Ds1}/0.6745) (26)；
C25、第二个判断条件为：
Figure GDA0002439491930000125
C3、卫星机动探测的综合评估：
若在c0个连续历元同时满足方程(23.a)和(27.a)，则认为卫星s1在历元t发生了机动；
若在历元t只有(23.a)或(27.a)满足，则认为卫星s1在历元t不可用；反之，则认为卫星s1在历元t可用。
步骤C3中，评估一整天的卫星机动时，若子网中有任何卫星发生机动，则认为该卫星发生了机动。
本发明的原理说明如下：
本设计涉及探测GNSS卫星机动以及精密定轨的解算方法，主要适用于提高机动卫星定轨精度以及机动卫星的有效利用等方面。本设计目的是克服现有技术中存在的只适用于区域观测网络、误判卫星机动，而且当消除接收机时钟钟差时，该技术易受电离层延迟、卫星机动和异常值的影响；或者现有技术中只适用于事后处理模式的缺陷与问题，构建了一种探测卫星机动的可靠的单历元解算方法(三步法)，提供了一种更可靠的基于广播星历和全球参考网实现全天候探测卫星机动的单历元解算方法，并为将来实时卫星机动提供一个可行的思路。
实施例：
参见图1，一种探测GNSS卫星轨道机动的方法，该方法包括以下步骤：
第一步：通过文件或者网络获取全球测站GNSS观测值和广播星历；
第二步：先将全球GNSS观测网分成六个子网(可以根据实际需要调整)，由于时间是相对的，再为每个子网选择一个装备有原子钟的测站作为基准站，其他测站作为参考站，然后将每个子网分线程执行第三步和第四步；
第三步：为子网内每个测站构建线性化后的消电离层组合平滑伪距观测方程：
Figure GDA0002439491930000131
式(1)中，
Figure GDA0002439491930000132
表示伪距残差，
Figure GDA0002439491930000133
表示视线方向的单位向量，Δs表示广播星历在径向dr、切向da、法向dc三个方向的误差，dtr表示接收机钟差，dts表示卫星钟差，c表示光速，ε表示噪声；
第四步、卫星机动探测：
A、根据广播星历提供的轨道和钟差简化式(1)得到式(2)，并根据式(2)计算基准站钟差；
Figure GDA0002439491930000134
式(2)中，m表示基准站接收机，
Figure GDA0002439491930000135
表示方差；则基准站钟差计算步骤如下：
A1、计算伪距残差
Figure GDA0002439491930000136
的中位数：
Figure GDA0002439491930000141
式(3)中，s＝1，2，...，N，N表示基准站观测到的卫星数，median{.}是计算序列
Figure GDA0002439491930000142
的中位数；
A2、计算伪距残差
Figure GDA0002439491930000143
的初始残差：
Figure GDA0002439491930000144
A3、计算伪距残差
Figure GDA0002439491930000145
的抗差方差因子：
σ0＝median{|vs|}/0.6745 (5)；
A4、计算伪距残差
Figure GDA0002439491930000146
的抗差等价权：
Figure GDA0002439491930000147
式(6)中，c0＝2，3，4或者c0＝5；
A5、计算基准站钟差：
Figure GDA0002439491930000148
在计算基准站钟差后，通过步骤B可以消除其他参考站的接收机钟差；
B、当已知基准站接收机钟差后，得到如下式(8)，然后消除其他参考站的接收机钟差；
Figure GDA0002439491930000149
式(8)中，
Figure GDA00024394919300001410
表示基准站m与其他参考站R之间的伪距残差之差，R＝1，2，...J，J表示子网内其他参考站的数量；则消除其他参考站钟差步骤如下：
B1、计算
Figure GDA0002439491930000151
的中位数：
Figure GDA0002439491930000152
B2、计算
Figure GDA0002439491930000153
的初始残差：
Figure GDA0002439491930000154
B3、计算
Figure GDA0002439491930000155
的抗差方差因子：
Figure GDA0002439491930000156
B4、计算
Figure GDA0002439491930000157
的抗差等价权：
Figure GDA0002439491930000158
B5、计算参考站接收机钟差：
Figure GDA0002439491930000159
B6、根据计算的基准站钟差和其他参考站接收机钟差得到如下方程：
Figure GDA00024394919300001510
Figure GDA00024394919300001511
B7、计算
Figure GDA0002439491930000161
和
Figure GDA0002439491930000162
经过上述步骤后，由于
Figure GDA0002439491930000163
和
Figure GDA0002439491930000164
都含有星历误差，因此可以利用这两组数据在步骤C中分别通过两个线程进行卫星机动判断；
C、假定s1为将要探测的卫星号，分别通过两个线程进行卫星机动判断：
C1、根据线程一得到卫星机动的第一个判断条件，线程一包括以下步骤：
C11、根据下式获取探测的变量
Figure GDA0002439491930000165
Figure GDA0002439491930000166
C12、计算
Figure GDA0002439491930000167
的中位数：
Figure GDA0002439491930000168
C13、计算
Figure GDA0002439491930000169
的初始残差：
Figure GDA00024394919300001610
C14、计算
Figure GDA00024394919300001611
的初始抗差方差因子：
Figure GDA00024394919300001612
C15、计算
Figure GDA00024394919300001613
的抗差因子和等价权：
Figure GDA00024394919300001614
Figure GDA00024394919300001615
C16、计算
Figure GDA00024394919300001616
的抗差方差因子：
Figure GDA0002439491930000171
C17、第一个判断条件为：
Figure GDA0002439491930000172
在得到第一个判断条件后，可以通过步骤C2得到第二个判断条件；
C2、根据线程二得到卫星机动的第二个判断条件，线程二包括以下步骤：
C21、根据广播星历的精度设置阈值Threshold1等于10米；
C22、根据式(15)计算Δs1(dr1 da1 dc1)；
C23、计算三维误差
Figure GDA0002439491930000173
Figure GDA0002439491930000174
C24、若卫星在连续十个历元没有发生机动或者异常情况，则得到一个序列如式(25)，而阈值Threshold2的计算如式(26)：
Figure GDA0002439491930000175
Threshold2＝c0(median{Ds1}/0.6745) (26)；
C25、第二个判断条件为：
Figure GDA0002439491930000176
得到两个判断条件后，在步骤C3中进行卫星机动探测的综合评估；
C3、卫星机动探测的综合评估：
若在c0个连续历元同时满足方程(23.a)和(27.a)，则认为卫星s1在历元t发生了机动；
若在历元t只有(23.a)或(27.a)满足，则认为卫星s1在历元t不可用；反之，则认为卫星s1在历元t可用；
通过以上的综合评估，就可以得到卫星机动信息；评估一整天的卫星机动时，若子网中有任何卫星发生机动，则认为该卫星发生了机动。
为了测试本设计三步法的性能，我们处理了来自图2所示的111个全球分布测站的观测值，如图2所示，每个圈代表一个子网，五角星代表每个子网中的基准站，GPS数据由国际GNSS服务(IGS)的地壳动力学数据信息系统(CDDIS)提供。此外，为了对本设计方法的性能进行综合评估，使用交叉验证方式，将本设计方法的结果与以CRX(由CODE提供)结尾的文件中的信息，来自导航中心的信息以及来自广播星历的健康标志进行比较；同时，也与斯坦福大学提出的传统方法进行了比较。为了进行说明，选取有可能发生卫星机动的两天的数据：2017年8月25日和2017年9月12日。此外，这两天的数据也用来被分析使用本设计方法和传统方法检测到的机动信息对精确定轨(POD)的影响，从而也可以从侧面验证本设计方法的优势，详情如下：
1、卫星G09测试(2017年8月25日)。图3和图4分别显示了本设计三步法的同步伪距残差和三维误差；图5显示了使用传统方法的同步伪距残差；表1总结了从四个验证方式和传统方法收集的G09状态。
表1 2017年8月25日卫星G09的验证结果
Figure GDA0002439491930000181
如图3-图5和表1所示，四种验证和传统方法都显示G09在2017年8月25日出现异常，表明G09在当天操纵。通过比较五种验证方法显示，广播星历提供了比其他验证方法更大的范围，这将导致G09可用的时间减少2小时；此外，由CODE提供的以CRX结尾的文件没有提供机动结束时间；然而，使用本设计三步法的机动起始时间与以CRX终止的文件所提供的时间类似，并且本设计三步法的机动终止时间接近导航中心提供的时间；与传统方法相比，使用本设计三步法的机动起始时间更接近CODE提供的机动起始时间。
2、卫星G07测试(2017年9月12日)。图6和图7分别显示本设计三步法的同步伪距残差和三维误差；图8显示了使用传统方法的同步伪距残差；表2总结了从四个验证方式和传统方法收集的G07状态。
表2 2017年9月12日卫星G07的验证结果
Figure GDA0002439491930000191
如图6-图8和表2所示，本设计三步法和广播星历的结果表明，G07只有两个不能使用的时间段。以CRX结尾的文件信息也显示G07没有机动，而导航中心和传统方法的信息显示有两个机动；另外，由于数据丢失，GPS信号可能没有被接收机接收。因此，卫星G07在2017年9月12日没有发生机动，这可能只是由于卫星时钟偏移的调整，而且在当天调整了两次。
3、卫星机动信息对精确定轨(POD)影响的全面测试。在这个测试中，将计算的轨道(在设计中称为IGG/IGG01轨道)与由CODE提供的最终轨道(在设计中称为COD轨道)进行比较，两种情况：1)使用的是本设计三步法所探测的机动信息；2)使用的是传统方法所探测的机动信息。其中，处理了2017年8月25日和2017年9月12日的IGS观测数据，图9是本设计实施例中精密定轨所用的测站分布图。
3.1、基于2017年8月25日数据的POD测试，将使用本设计三步法检测到的机动信息的计算轨道(IGG)与X、Y和Z方向上的CODE轨道进行比较；图10和图11显示了与差异相对应的均方根(RMS)和标准偏差(STD)值；图12显示使用传统方法检测到的机动信息的计算轨道IGG01与COD之间的轨道差异；G09的统计结果列于表3中。
表3卫星G09的轨道差异统计结果
Figure GDA0002439491930000201
图10-图11显示，使用本设计三步法检测到的机动信息导致在计算的轨道(IGG)与CODE的轨道之间的差异的RMS和STD值在X、Y和Z方向小于3cm。相反，图12和表3表明，使用传统的方法检测到的机动信息导致在计算的轨道(IGG01)和CODE的轨道差异的RMS和STD值在X、Y和Z方向大于1m。
3.2、基于2017年9月12日数据的POD测试，将计算轨道(IGG)与CODE在X、Y和Z方向的轨道进行比较。图13和图14显示出了相应的RMS和STD值，由于这里主要验证G07号卫星在2017年9月12日没有发生机动，因此，在精密定轨时不需要分段处理，若此精密定轨结果正常，则该颗卫星确实没有发生机动。
如图13-图14所示，2017年9月12日的G07卫星的确不需要分成两部分来确定轨道，而且解算精度也是和其他卫星一样的；相反，说明传统的方法误判了机动。
测试结果表明利用本设计三步法可以全天候提供具有单历元解特征的可靠的卫星机动信息，并且可以为未来实时可靠的机动检测提供了另一个可行的思路。

Claims (1)

1.一种探测GNSS卫星轨道机动的方法，其特征在于，该方法包括以下步骤：
第一步：通过文件或者网络获取全球测站GNSS观测值和广播星历；
第二步：先将全球GNSS观测网分成多个子网，再为每个子网选择一个装备有原子钟的测站作为基准站，其他测站作为参考站，然后将每个子网分线程执行第三步和第四步；
第三步：为子网内每个测站构建线性化后的消电离层组合平滑伪距观测方程：
Figure FDA0002439491920000011
式(1)中，
Figure FDA0002439491920000012
表示伪距残差，
Figure FDA0002439491920000013
表示视线方向的单位向量，Δs表示广播星历在径向dr、切向da、法向dc三个方向的误差，dtr表示接收机钟差，dts表示卫星钟差，c表示光速，ε表示噪声；
第四步、卫星机动探测：
A、根据广播星历提供的轨道和钟差简化式(1)得到式(2)，并根据式(2)计算基准站钟差；
Figure FDA0002439491920000014
式(2)中，m表示基准站接收机，
Figure FDA0002439491920000015
表示方差；则基准站钟差计算步骤如下：
A1、计算伪距残差
Figure FDA0002439491920000016
的中位数：
Figure FDA0002439491920000017
式(3)中，s＝1，2，...，N，N表示基准站观测到的卫星数，median{.}是计算序列
Figure FDA0002439491920000018
的中位数；
A2、计算伪距残差
Figure FDA0002439491920000019
的初始残差：
Figure FDA0002439491920000021
A3、计算伪距残差
Figure FDA0002439491920000022
的抗差方差因子：
σ0＝median{|vs|}/0.6745 (5)；
A4、计算伪距残差
Figure FDA0002439491920000023
的抗差等价权：
Figure FDA0002439491920000024
式(6)中，c0＝2，3，4或者c0＝5；
A5、计算基准站钟差：
B、当已知基准站接收机钟差后，得到如下式(8)，然后消除其他参考站的接收机钟差；
Figure FDA0002439491920000026
式(8)中，
Figure FDA0002439491920000027
表示基准站m与其他参考站R之间的伪距残差之差，R＝1，2，...J，J表示子网内其他参考站的数量；则消除其他参考站钟差步骤如下：
B1、计算
Figure FDA0002439491920000028
的中位数：
Figure FDA0002439491920000029
B2、计算
Figure FDA00024394919200000210
的初始残差：
Figure FDA0002439491920000031
B3、计算
Figure FDA0002439491920000032
的抗差方差因子：
Figure FDA0002439491920000033
B4、计算
Figure FDA0002439491920000034
的抗差等价权：
Figure FDA0002439491920000035
B5、计算参考站接收机钟差：
Figure FDA0002439491920000036
B6、根据计算的基准站钟差和其他参考站接收机钟差得到如下方程：
Figure FDA0002439491920000037
Figure FDA0002439491920000038
B7、计算
Figure FDA0002439491920000039
和
Figure FDA00024394919200000310
C、假定s1为将要探测的卫星号，分别通过两个线程进行卫星机动判断：
C1、根据线程一得到卫星机动的第一个判断条件，线程一包括以下步骤：
C11、根据下式获取探测的变量
Figure FDA00024394919200000311
Figure FDA0002439491920000041
C12、计算
Figure FDA0002439491920000042
的中位数：
Figure FDA0002439491920000043
C13、计算
Figure FDA0002439491920000044
的初始残差：
Figure FDA0002439491920000045
C14、计算
Figure FDA0002439491920000046
的初始抗差方差因子：
Figure FDA0002439491920000047
C15、计算
Figure FDA0002439491920000048
的抗差因子和等价权：
Figure FDA0002439491920000049
Figure FDA00024394919200000410
C16、计算
Figure FDA00024394919200000411
的抗差方差因子：
Figure FDA00024394919200000412
C17、第一个判断条件为：
Figure FDA00024394919200000413
C2、根据线程二得到卫星机动的第二个判断条件，线程二包括以下步骤：
C21、根据广播星历的精度设置阈值Threshold1等于10米；
C22、根据式(15)计算Δs1(dr1 da1 dc1)；
C23、计算三维误差
Figure FDA0002439491920000051
Figure FDA0002439491920000052
C24、若卫星在连续十个历元没有发生机动或者异常情况，则得到一个序列如式(25)，而阈值Threshold2的计算如式(26)：
Figure FDA0002439491920000053
Threshold2＝c0(median{Ds1}/0.6745) (26)；
C25、第二个判断条件为：
Figure FDA0002439491920000054
C3、卫星机动探测的综合评估：
若在c0个连续历元同时满足方程(23.a)和(27.a)，则认为卫星s1在历元t发生了机动；
若在历元t只有(23.a)或(27.a)满足，则认为卫星s1在历元t不可用；反之，则认为卫星s1在历元t可用；
评估一整天的卫星机动时，若子网中有任何卫星发生机动，则认为该卫星发生了机动。
CN201810161688.XA 2018-02-27 2018-02-27 一种探测gnss卫星轨道机动的方法 Active CN108535746B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810161688.XA CN108535746B (zh) 2018-02-27 2018-02-27 一种探测gnss卫星轨道机动的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810161688.XA CN108535746B (zh) 2018-02-27 2018-02-27 一种探测gnss卫星轨道机动的方法

Publications (2)

Publication Number Publication Date
CN108535746A CN108535746A (zh) 2018-09-14
CN108535746B true CN108535746B (zh) 2020-07-21

Family

ID=63485835

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810161688.XA Active CN108535746B (zh) 2018-02-27 2018-02-27 一种探测gnss卫星轨道机动的方法

Country Status (1)

Country Link
CN (1) CN108535746B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111308504B (zh) * 2019-11-26 2023-10-31 中国科学院国家授时中心 一种基于相位观测值实时探测北斗卫星轨道机动的方法
CN111505677A (zh) * 2020-04-15 2020-08-07 中国科学院国家授时中心 一种基于地面参考站观测的geo卫星轨道机动修复方法
CN111959828B (zh) * 2020-10-21 2020-12-29 中国人民解放军国防科技大学 基于非线性偏差演化的航天器轨道机动检测方法和装置
CN113536197A (zh) * 2021-07-13 2021-10-22 中国科学院国家授时中心 一种卫星轨道机动时段的探测方法及系统
CN114690210A (zh) * 2022-04-18 2022-07-01 山东大学 一种基于多普勒观测值的北斗卫星机动探测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2088448A2 (en) * 2008-02-06 2009-08-12 Broadcom Corporation Method and apparatus for improving accuracy and/or integrity of long-term-orbit information for a global-navigation-satellite system
CN103869344A (zh) * 2012-12-13 2014-06-18 东莞市泰斗微电子科技有限公司 一种抗差估计方法
CN105629278A (zh) * 2014-11-21 2016-06-01 桂林电子科技大学 一种高精度gnss伪距单点定位的互差中值加权定位方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106885569A (zh) * 2017-02-24 2017-06-23 南京理工大学 一种强机动条件下的弹载深组合arckf滤波方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2088448A2 (en) * 2008-02-06 2009-08-12 Broadcom Corporation Method and apparatus for improving accuracy and/or integrity of long-term-orbit information for a global-navigation-satellite system
CN103869344A (zh) * 2012-12-13 2014-06-18 东莞市泰斗微电子科技有限公司 一种抗差估计方法
CN105629278A (zh) * 2014-11-21 2016-06-01 桂林电子科技大学 一种高精度gnss伪距单点定位的互差中值加权定位方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A Real-Time Robust Method to Detect BeiDou GEO/IGSO Orbital Maneuvers;Guanwen Huang;《sensors》;20171129;第1-17页 *
A Robust Method to Detect BeiDou Navigation Satellite System Orbit Maneuvering/Anomalies and Its Applications to Precise Orbit Determination;Fei Ye;《sensors》;20170516;第1-17页 *
BDS广播星历的实时异常探测;刘朝英;《第五届中国卫星导航学术年会电子文集》;20140306;第1-6页 *
GNSS卫星运行异常状态监测及实现;范顺西;《第七届中国卫星导航学术年会论文集》;20160531;第1-8页 *

Also Published As

Publication number Publication date
CN108535746A (zh) 2018-09-14

Similar Documents

Publication Publication Date Title
CN108535746B (zh) 一种探测gnss卫星轨道机动的方法
El‐Mowafy et al. Integrity monitoring for positioning of intelligent transport systems using integrated RTK‐GNSS, IMU and vehicle odometer
CN110187364B (zh) 一种低轨导航增强精密改正数据生成、上注系统及方法
CN108508461B (zh) 基于gnss载波相位高精度定位完好性监测方法
US5917445A (en) GPS multipath detection method and system
Ochieng et al. Urban road transport navigation: performance of the global positioning system after selective availability
Heng Safe satellite navigation with multiple constellations: global monitoring of GPS and GLONASS signal-in-space anomalies
CN101419275B (zh) 基于多接收机的局域机场监视方法和系统
Hassan et al. A review of system integration and current integrity monitoring methods for positioning in intelligent transport systems
CN106813664A (zh) 一种车载导航方法和装置
Li et al. Review of PPP–RTK: Achievements, challenges, and opportunities
CN108387169B (zh) 一种基于实时大气产品的gnss形变监测系统
CN110221320B (zh) 一种基于抛物面天线观测的北斗频间偏差测定方法
CN108490459A (zh) 精度与风险均衡应用于gnss位置服务的方法及系统
CN112114341B (zh) 低轨卫星协同测频无源定位方法
CN115767430A (zh) 一种基于北斗的石化领域精准时空信息处理与服务系统
Melgard et al. G2-the first real-time GPS and GLONASS precise orbit and clock service
CN110458089B (zh) 一种基于高低轨光学卫星观测的海上目标关联系统及方法
CN115826016A (zh) 一种北斗双频星基增强改正数及完好性参数解算的方法
Paziewski et al. Enhanced wide-area multi-GNSS RTK and rapid static positioning in the presence of ionospheric disturbances
Elsayed et al. Bounding of correlated double-differenced GNSS observation errors using NRTK for precise positioning of autonomous vehicles
Gaglione et al. Benefit of GNSS multiconstellation in position and velocity domain
Seepersad Improving reliability and assessing performance of global navigation satellite system precise point positioning ambiguity resolution
Lee et al. Gnss fault monitoring using android devices
Adjrad et al. 3D-mapping-aided GNSS exploiting Galileo for better accuracy in dense urban environments

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant