WO2021082189A1 - 基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法 - Google Patents

基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法 Download PDF

Info

Publication number
WO2021082189A1
WO2021082189A1 PCT/CN2019/123837 CN2019123837W WO2021082189A1 WO 2021082189 A1 WO2021082189 A1 WO 2021082189A1 CN 2019123837 W CN2019123837 W CN 2019123837W WO 2021082189 A1 WO2021082189 A1 WO 2021082189A1
Authority
WO
WIPO (PCT)
Prior art keywords
satellite
maneuvering
orbit
beidou satellite
parameters
Prior art date
Application number
PCT/CN2019/123837
Other languages
English (en)
French (fr)
Inventor
许小龙
赵齐乐
周泉
Original Assignee
中海北斗深圳导航技术有限公司
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 中海北斗深圳导航技术有限公司 filed Critical 中海北斗深圳导航技术有限公司
Publication of WO2021082189A1 publication Critical patent/WO2021082189A1/zh

Links

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/23—Testing, monitoring, correcting or calibrating of receiver elements

Definitions

  • the invention relates to the field of GNSS satellite-based augmentation systems, in particular to the field of Beidou satellite precision orbit calculation and in-orbit abnormal maneuver detection, and specifically refers to a Beidou satellite maneuver and anomaly detection method based on precision orbit and ground station data.
  • GNSS satellite orbit resolution accuracy is the core link of the entire high-precision positioning, and its resolution accuracy is greatly affected by satellite orbit maneuvering, but from the perspective of satellite systems, orbit maintenance is an indispensable link. Therefore, a more reliable mobile detection method has become a problem that needs to be solved for GNSS precise orbit determination.
  • the BeiDou Navigation Satellite System is a satellite navigation system independently developed in China.
  • the system construction is divided into three stages.
  • the first stage system consists of 2 Geostationary Earth Orbit satellites (GEO);
  • the second stage system consists of 5 GEO satellites, 5 Inclined GeoSynchronous Orbit satellites (IGSO) and 4 medium-orbiting earths Orbiting satellite (Medium Earth Orbit satellite, MEO) composition.
  • the third stage system consists of 3 GEO, 3 IGSO and 24 MEO satellites. Due to the complex composition of the Beidou system constellation and frequent orbit maintenance, the active injection of fuel by the satellite during satellite maneuvering will destroy the force properties of the satellite in space and will seriously affect the precision of precise orbit determination.
  • the purpose of the present invention is to detect the maneuver time of the Beidou satellite system and provide a method that can accurately determine the satellite maneuver time.
  • the method uses the Beidou satellite system's post-precision orbit and ground station observation data to detect the start and end of the Beidou satellite maneuver At the moment, the maneuvering magnitude can be judged by passing.
  • This method is also suitable for the detection of abnormal maneuvers in other satellite navigation systems, and has good generalization.
  • the Beidou satellite maneuvering and anomaly detection method based on the precision orbit and ground station data of the present invention has the following composition: the Beidou satellite maneuvering and anomaly detection method based on the precision orbit and ground station data, its main features are Combining the predicted precision orbit of the Beidou satellite with the satellite position calculated by the ground station, it is determined whether the Beidou satellite is maneuvering or abnormal in its orbit.
  • the method includes the following steps:
  • step S4 Perform long-term orbit integration on the state parameters obtained in step S3 to obtain the Beidou satellite space position parameters for a period of time in the integration direction;
  • step S6 Compare the satellite space coordinate parameters in step S4 and step S5, convert the coordinate difference to the orbital coordinate system, and determine whether the Beidou satellite has abnormal maneuvers based on a priori parameters. If a maneuver occurs, give the time when the maneuver occurred and the size of the maneuver. ;
  • the completion time of the Beidou satellite maneuvering is calculated by integrating the precise orbit forward after the maneuvering, and repeating steps S4 to S6.
  • the orbit integration method according to step S1 is characterized in that the 6/7th order of Runge-Kutta integration method is used to calculate the initial 10 epoch positions, and then the Adams integration method is used for subsequent epoch continuous integration;
  • the sixth-order and seventh-order formulas of Runge-Kutta integration method can be expressed as:
  • the explicit formula and the implicit formula are used at the same time, and the explicit formula is first calculated The approximate value on the step point, and then the approximate value is corrected by the implicit formula to get the required .
  • step S5 the method for calculating the satellite spatial position based on the globally distributed Beidou ground continuous tracking stations is characterized in that the Beidou satellite is solved by the fixed ground station ground coordinates and the receiver clock error parameters
  • the spatial position, using the Beidou satellite pseudorange and phase dual-frequency observation equations are as follows:
  • the observation method to obtain the partial derivative of the satellite position parameter can be expressed as:
  • M is the tropospheric projection function, with They are the approximate value of the parameter to be estimated and the random error calculated from the original observation; after linearization, the parameter to be estimated takes the distance length as the unit, taking into account the initial state parameters of the satellite, as follows:
  • step S5 the coordinate parameters in the ground-fixed coordinate system are converted to the J2000.0 inertial coordinate system, and the implementation process is as follows:
  • Is the precession matrix Is the nutation matrix
  • Is the earth rotation matrix Is the pole shift matrix.
  • the matrix can be accurately determined according to the time of the epoch;
  • the matrix can be calculated according to 1980IAU nutation theory;
  • the matrix is obtained by interpolating IAT-UT1 from the time bulletin of IERS;
  • the extreme value of the matrix needs to be obtained from the difference of extreme values published by IERS.
  • step S6 the conversion of the coordinate difference into the orbital coordinate system is implemented as follows:
  • the unit vector of the coordinate axis of the orbital coordinate system in the inertial system can use the satellite position in the inertial system With velocity vector To mean:
  • step S6 the given prior parameters are used to determine whether the Beidou satellite has abnormal maneuvers, and the implementation process is as follows:
  • the maneuvering direction and size of the coordinates in the orbital coordinate system can be visually displayed.
  • the criterion for judging the occurrence of satellite maneuvering is that the R direction of the difference between the epochs before and after the coordinate difference is greater than 1 meter. Or the N direction is greater than 5 meters, or the T direction is greater than 10 meters, the formula is expressed as:
  • s represents the state of maneuvering, where 1 means maneuvering occurs, and 0 means no maneuvering occurs; It is the coordinate difference calculated between the previous epoch and the previous epoch, which can be expressed as:
  • R , N and T represent the radial, normal and tangential directions of the orbital coordinate system, respectively.
  • the beneficial effect of the present invention is to detect the maneuvering time of the Beidou satellite system and provide a method that can accurately determine the satellite maneuvering time.
  • the method uses the Beidou satellite system's post-precision orbit and ground station observation data to detect the start and At the end time, the maneuvering magnitude can be judged by the pass. This method is also suitable for the detection of abnormal maneuvers in other satellite navigation systems, and has good generalization.
  • Figure 1 is a step diagram of Beidou satellite maneuvering and anomaly detection method based on precise orbit and ground station data
  • Figure 2 is a schematic diagram of the orbit integration method
  • Figure 3 is a schematic diagram of the calculation method of Beidou satellite space position.
  • the present invention creates the use of Beidou precision orbit afterwards and ground station tracking data to detect the abnormal maneuvering of Beidou satellites in orbit, mainly through the following technical solutions:
  • the numerical integration method uses Runge-Kutta 6/7 order single-step integration method to start, and the integral formula is:
  • step S1 the 6/7th order of Runge-Kutta integration method is used to calculate the initial 10 epoch positions, and then the Adams integration method is used to perform continuous integration of subsequent epochs;
  • step S3 Use the satellite position parameters in step S2 and the state transition matrix to fit the state parameters in (1) to obtain more accurate state parameters, where the modified orbit adopts post-precision orbits, and more accurate orbit state parameters can be obtained through this step;
  • Orbital integration method adopts Runge-Kutta Start with 6/7 single-step method, and then integrate by Adams multi-step method to calculate the orbit of the continuous interval in the future.
  • the orbit integration process is shown in Figure 2;
  • M is the tropospheric projection function, with They are the approximate value of the parameter to be estimated and the random error calculated from the original observation. After linearization, the parameters to be estimated are in the unit of distance length, taking into account the initial state parameters of the satellite, as follows:
  • step S5 the coordinate parameters in the ground-fixed coordinate system are converted to the J2000.0 inertial coordinate system, and the realization process is as follows: r represents the coordinate of a point in the protocol geocentric inertial system corresponding to epoch J2000.0 , With R as the coordinates in the fixed coordinate system of the protocol earth, there are:
  • Is the precession matrix Is the nutation matrix
  • Is the earth rotation matrix Is the pole shift matrix
  • the matrix can be accurately determined according to the time of the epoch;
  • the matrix can be calculated according to 1980IAU nutation theory;
  • the matrix is obtained by interpolating IAT-UT1 from the time bulletin of IERS;
  • the extreme value of the matrix needs to be obtained from the difference of extreme values published by IERS.
  • step S6 the coordinate difference is converted into the orbital coordinate system, and the implementation process is as follows:
  • the unit vector of the coordinate axis of the orbital coordinate system in the inertial system can use the satellite position in the inertial system With velocity vector To mean:
  • the given a priori parameter determines whether the Beidou satellite has an abnormal maneuver, and if a maneuver occurs, the time and the size of the maneuver are given; the implementation process of judging whether the Beidou satellite has an abnormal maneuver is as follows:
  • the maneuvering direction and size of the coordinates in the orbital coordinate system can be visually displayed.
  • the criterion for judging the occurrence of satellite maneuvering is that the R direction of the difference between the epochs before and after the coordinate difference is greater than 1 meter. Or the N direction is greater than 5 meters, or the T direction is greater than 10 meters, the formula is expressed as:
  • s represents the state of maneuvering, where 1 means maneuvering occurs, and 0 means no maneuvering occurs; It is the coordinate difference calculated between the previous epoch and the previous epoch, which can be expressed as:
  • R , N and T represent the radial, normal and tangential directions of the orbital coordinate system, respectively.
  • steps S4 to S6 are repeated to calculate the end time of the Beidou satellite maneuver by integrating the precise orbit forward after the maneuver is completed.

Landscapes

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

Abstract

涉及GNSS星基增强系统领域，尤其涉及北斗卫星精密轨道解算及在轨运行异常机动探测领域，具体是指一种基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法，其主要特点是结合预报的北斗卫星精密轨道与地面台站计算得到的卫星位置做比对，判定北斗卫星在轨运行是否发生机动及异常，由此能够准确确定卫星机动时间。利用北斗卫星系统事后精密轨道和地面台站观测数据，探测北斗卫星机动开始及结束时刻，能够对机动量级进行判别；同样适用于其他卫星导航系统异常机动探测，具有较好的推广性。

Description

基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法 技术领域
本发明涉及GNSS星基增强系统领域，尤其涉及北斗卫星精密轨道解算及在轨运行异常机动探测领域，具体是指一种基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法。
 背景技术
GNSS技术发展至今，一直致力于为用户提供实时、高精度、高可靠性的导航授时服务。传统实时定位采用广播星历解算卫星轨道、钟差参数，实时解算用户空间位置参数，但是受限于广播星历精度较差，用户仅能达到米级定位精度。为了满足高精度定位用户需求，对GNSS高精度卫星轨道求解提出了较高的要求，广域高精度定位用户通过引入GNSS精密轨道作为空间基准，可消除卫星轨道部分引入的误差，提高用户定位解算精度，GNSS卫星轨道解算精度是整个高精度定位的核心环节，其解算精度受卫星轨道机动影响较大，但从卫星系统角度来看，轨道维持是必不可少的环节。因此，较为可靠的机动探测方法成为GNSS精密定轨需要解决问题。
北斗卫星导航系统（BeiDou Navigation Satellite System，BDS）是中国独立发展的卫星导航系统，系统建设分为了三个阶段。第一阶段系统由2颗地球同步卫星（Geostationary Earth Orbit satellites，GEO）组成；第二阶段系统由5颗GEO卫星、5颗倾斜地球同步卫星（Inclined GeoSynchronous Orbit satellite，IGSO）和4颗中轨地球轨道卫星（Medium Earth Orbit satellite，MEO）组成。；第三阶段系统由3颗GEO、3颗IGSO和24颗MEO卫星组成。由于北斗系统星座构成复杂，轨道维持频繁，卫星机动时卫星主动喷射燃料将破坏卫星空间受力属性，将严重影响精密定轨精度。
准确探测北斗卫星在轨运行期间的机动事件，对于机动期间卫星精度定轨精度提高有一定的作用。
 技术问题
由于北斗系统星座构成复杂，轨道维持频繁，卫星机动时卫星主动喷射燃料将破坏卫星空间受力属性，将严重影响精密定轨精度。
 技术解决方案
本发明的目的是针对北斗卫星系统的机动时间进行探测，提供一种能够准确确定卫星机动时间的方法，该方法利用北斗卫星系统事后精密轨道和地面台站观测数据，探测北斗卫星机动开始及结束时刻，通过能够对机动量级进行判别。该方法同样适用于其他卫星导航系统异常机动探测，具有较好的推广性。
为了实现上述目的，本发明的基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法具有如下构成：该基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法，其主要特点是结合预报的北斗卫星精密轨道与地面台站计算得到的卫星位置做比对，判定北斗卫星在轨运行是否发生机动及异常。
所述的方法包括以下步骤：
S1、从事后解算的北斗卫星精密轨道中获取参考时刻卫星状态参数，包括：卫星在J2000.0惯性坐标系下空间位置参数
Figure 771777dest_path_image001
，空间速度参数
Figure 105807dest_path_image002
；
S2、对参考时刻北斗卫星状态参数轨道积分，获取连续弧段的北斗卫星轨道参数及各时刻卫星位置与步骤（1）中状态参数的状态转移矩阵；
S3、利用步骤S2中的卫星位置参数和状态转移矩阵拟合S1中的状态参数，获取较为较为准确的状态参数；
S4、对步骤S3中获的状态参数进行长时间轨道积分，求得积分方向一段时间北斗卫星空间位置参数；
S5、收集能够跟踪到北斗卫星的地面台站观测数据，利用其精确已知的坐标参数和地面钟差参数，计算北斗卫星空间位置，并将地固坐标系中坐标参数转换到惯性坐标系统中，其中地固系定义如下，坐标系原点与包含海洋和大气层在内的整个地球质心重合，坐标轴指向按照BIH1984.0确定的指向，Z轴与BIH1984.0瞬时自转轴重合，X轴指向BIH1984.0时刻本初子午线与地球赤道交点，Y轴与Z轴和X轴组成右手坐标系；
S6、对比步骤S4和步骤S5中卫星空间坐标参数，将坐标差异转换到轨道坐标系中，给定先验参数，判别北斗卫星是否发生异常机动，如果发生机动则给出机动发生时刻、机动大小；
S7、当检测出机动时间发生时，通过对机动结束后的精密轨道向前积分，重复步骤S4至S6，计算出北斗卫星机动结束时刻。
优选的，根据步骤S1所述的轨道积分方法，其特征在于采用Runge-Kutta积分方法6/7阶计算初始10个历元位置，随后采用Adams积分方法进行后续历元连续积分；
Runge-Kutta积分方法6阶和7阶公式可表示为：
Figure 838139dest_path_image003
和
Figure 890409dest_path_image004
相应的系数如表所示：
Figure 749912dest_path_image005
Adams积分显式公式可表示为：
 
Figure 938448dest_path_image006
Adams积分隐式公式可表示为：
 
Figure 841682dest_path_image007
在轨道数值积分计算时，显式公式与隐式公式同时采用，先由显式公式计算出
Figure 381247dest_path_image008
步点上近似值，再由隐式公式校正该近似值，得出需要的
Figure 775932dest_path_image009
。
优选的，步骤S5中，所述的根据全球分布的北斗地面连续跟踪台站计算卫星空间位置的方法，其特征在于，通过固定地面台站地固系坐标和接收机钟差参数，求解北斗卫星空间位置，利用北斗卫星伪距和相位双频观测值方程如下：
Figure 84554dest_path_image010
其中：
Figure 158689dest_path_image011
和
Figure 185551dest_path_image012
分别表示载波和伪距无电离组合观测值，
Figure 121277dest_path_image013
为卫星至测站的几何距离，
Figure 284405dest_path_image014
为真空中光速，
Figure 795021dest_path_image015
和
Figure 43599dest_path_image016
分别表示吸收了接收机端或卫星端硬件延迟的钟差参数，
Figure 783016dest_path_image017
为对流层延迟参数，
Figure 925285dest_path_image018
为其它需要顾及的改正，
Figure 216589dest_path_image019
为无电离层组合相位观测值波长，
Figure 218043dest_path_image020
为无电离层组合相位观测值模糊度，可表示为：
 
Figure 761151dest_path_image021
Figure 898871dest_path_image022
为载波相位和伪距测量值的误差项；观测方法对卫星位置参数求偏导数可表示为：
 
Figure 485710dest_path_image023
其中：
Figure 708881dest_path_image024
，
Figure 623109dest_path_image025
， M为对流层投影函数，
Figure 474390dest_path_image026
和
Figure 638655dest_path_image027
分别为待估参数近似值和原始观测计算得到的随机误差；经线性化后，待估参数以距离长度为单位，顾及卫星初始状态参数，有：
Figure 349122dest_path_image028
 
其中设计矩阵
Figure 968453dest_path_image029
和待估参数可表示为：
Figure 80766dest_path_image030
和
Figure 9408dest_path_image031
和
Figure 207171dest_path_image032
至此，观测方程线性化完成，通过线性估计方法进行待估参数的求解。
进一步，步骤S5中，所述的并将地固坐标系中坐标参数转换到J2000.0惯性坐标系统中，其实现过程如下：
以 r表示某点在历元J2000.0对应的协议地心惯性系中的坐标，以 R表示协议地球固定坐标系中的坐标，则有：
Figure 895772dest_path_image033
 其中：
Figure 862591dest_path_image034
为岁差矩阵；
Figure 227713dest_path_image035
为章动矩阵；
Figure 647193dest_path_image036
为地球自转矩阵；
Figure 873907dest_path_image037
为极移矩阵。上述转换矩阵中，
Figure 695232dest_path_image038
矩阵根据历元的时刻可以精确确定；
Figure 231256dest_path_image039
矩阵可根据1980IAU章动理论来计算；
Figure 138032dest_path_image040
矩阵通过从IERS的时间公报中内插求出IAT-UT1进而求得；
Figure 293070dest_path_image037
矩阵需要的极移值需从由IERS公布的极移值差值得到。
优选的，步骤S6中，所述将坐标差异转换到轨道坐标系中，其实现过程如下：
轨道坐标系的坐标轴在惯性系中的单位矢量可以借助于惯性系中的卫星位置
Figure 106917dest_path_image041
与速度向量
Figure 689208dest_path_image042
来表示，即：
Figure 207914dest_path_image043
 若轨道坐标系中某一点坐标为
Figure 510851dest_path_image044
，则该点在惯性系中坐标
Figure 41189dest_path_image045
为：
 
Figure 653436dest_path_image046
。
进一步，步骤S6中，所述给定先验参数，判别北斗卫星是否发生异常机动，其实现过程如下：
将计算得到的卫星坐标差异转换到轨道坐标系后，即可直观显示坐标在轨道坐标系下的机动方向和大小，判断卫星机动发生的标准为坐标差异前后历元间差分 R方向大于1米，或 N方向大于5米，或 T方向大于10米，公式表示为：
Figure 534805dest_path_image047
 
式中， s表示机动发生状态，其中1表示发生机动，0表示未发生机动；
Figure 907011dest_path_image048
为前后历元间计算的坐标差异，可表示为：
 
Figure 557435dest_path_image049
R 、 N 、 T分别代表轨道坐标系的径向、法向和切向。
有益效果
本发明的有益效果是针对北斗卫星系统的机动时间进行探测，提供一种能够准确确定卫星机动时间的方法，该方法利用北斗卫星系统事后精密轨道和地面台站观测数据，探测北斗卫星机动开始及结束时刻，通过能够对机动量级进行判别。该方法同样适用于其他卫星导航系统异常机动探测，具有较好的推广性。
 附图说明
图1是基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法步骤图；
图2是轨道积分方法示意图；
图3是北斗卫星空间位置计算方法示意图。
 本发明的最佳实施方式
为了能够清楚地描述本发明的技术内容，下面结合具体实施案例来进行进一步描述。
如图1所示，本发明创造了使用北斗事后精密轨道和地面台站跟踪数据进行北斗卫星在轨运行异常机动探测，主要通过以下技术方案来实现：
S1、从事后解算的北斗卫星精密轨道中获取参考时刻卫星状态参数，包括：卫星在J2000.0惯性坐标系下空间位置参数
Figure 340584dest_path_image001
，空间速度参数
Figure 709248dest_path_image050
；
S2、对参考时刻北斗卫星状态参数轨道积分，获取连续弧段的北斗卫星轨道参数及各时刻卫星位置关于初始状态参数的转移矩阵；
其中，数值积分方法采用Runge-Kutta 6/7阶单步积分法起步，积分公式为：
 
Figure 619566dest_path_image051
和
Figure 390076dest_path_image052
 相应的系数如表所示：
Figure 344126dest_path_image053
。
优选的，根据步骤S1中，采用Runge-Kutta积分方法6/7阶计算初始10个历元位置，随后采用Adams积分方法进行后续历元连续积分；
S3、利用步骤S2中的卫星位置参数和状态转移矩阵拟合（1）中的状态参数，获取较为准确的状态参数，其中修正轨道采用事后精密轨道，通过本步骤可以获得较为准确轨道状态参数；
S4、对步骤S3中获的精确轨道状态参数进行长时间轨道积分，求得积分方向一段时间内北斗卫星空间位置参数，实现未来一段时间的轨道预报。轨道积分方法采用Runge-Kutta 6/7单步法起步，随后用Adams多步法积分，计算未来连续区间轨道，轨道积分过程如图2所示；
S5，收集跟踪到北斗卫星的地面台站观测数据，利用其精确已知的坐标参数和地面钟差参数，计算北斗卫星空间位置，过程如图3所示，并将地固坐标系中坐标参数转换到J2000.0惯性坐标系统中。利用北斗卫星伪距和相位双频观测值方程如下：
 
Figure 200086dest_path_image054
其中：
Figure 911166dest_path_image055
和
Figure 270603dest_path_image056
分别表示载波和伪距无电离组合观测值，
Figure 661133dest_path_image057
为卫星至测站的几何距离，
Figure 348598dest_path_image058
为真空中光速，
Figure 256511dest_path_image059
和
Figure 736034dest_path_image060
分别表示吸收了接收机端或卫星端硬件延迟的钟差参数，
Figure 31886dest_path_image061
为对流层延迟参数，
Figure 472225dest_path_image062
为其它需要顾及的改正，
Figure 652671dest_path_image063
为无电离层组合相位观测值波长，
Figure 111334dest_path_image064
为无电离层组合相位观测值模糊度，可表示为：
 
Figure 719033dest_path_image065
Figure 646669dest_path_image066
为载波相位和伪距测量值的误差项。观测方法对卫星位置参数求偏导数可表示为：
 
Figure 365226dest_path_image067
其中：
Figure 678396dest_path_image068
，
Figure 456996dest_path_image069
，
M为对流层投影函数，
Figure 868998dest_path_image070
和
Figure 391247dest_path_image071
分别为待估参数近似值和原始观测计算得到的随机误差。经线性化后，待估参数以距离长度为单位，顾及卫星初始状态参数，有：
 
Figure 824502dest_path_image072
其中设计矩阵
Figure 774003dest_path_image073
和待估参数可表示为：
Figure 410652dest_path_image074
Figure 736591dest_path_image075
Figure 24353dest_path_image076
 
 至此，观测方程线性化完成，通过线性估计方法进行待估参数的求解。
进一步，步骤S5中，将地固坐标系中坐标参数转换到J2000.0惯性坐标系统中，其实现过程如下：以 r表示某点在历元J2000.0对应的协议地心惯性系中的坐标，以 R表示协议地球固定坐标系中的坐标，则有： 
Figure 144756dest_path_image077
其中：
Figure 534280dest_path_image078
为岁差矩阵；
Figure 398331dest_path_image079
为章动矩阵；
Figure 540599dest_path_image080
为地球自转矩阵；
Figure 97483dest_path_image081
为极移矩阵；上述转换矩阵中，
Figure 708724dest_path_image078
矩阵根据历元的时刻可以精确确定；
Figure 376465dest_path_image079
矩阵可根据1980IAU章动理论来计算；
Figure 373240dest_path_image080
矩阵通过从IERS的时间公报中内插求出IAT-UT1进而求得；
Figure 101025dest_path_image081
矩阵需要的极移值需从由IERS公布的极移值差值得到。
 优选的，步骤S6中,将坐标差异转换到轨道坐标系中，其实现过程如下：
轨道坐标系的坐标轴在惯性系中的单位矢量可以借助于惯性系中的卫星位置
Figure 196632dest_path_image082
与速度向量
Figure 402486dest_path_image083
来表示，即：
Figure 519346dest_path_image084
 
若轨道坐标系中某一点坐标为
Figure 418032dest_path_image085
，则该点在惯性系中坐标
Figure 738286dest_path_image086
为：
 
Figure 747830dest_path_image087
。
优选的，步骤S6中，所述给定先验参数，判别北斗卫星是否发生异常机动，如果发生机动则给出机动发生时刻、机动大小；判别北斗卫星是否发生异常机动，其实现过程如下：
将计算得到的卫星坐标差异转换到轨道坐标系后，即可直观显示坐标在轨道坐标系下的机动方向和大小，判断卫星机动发生的标准为坐标差异前后历元间差分 R方向大于1米，或 N方向大于5米，或 T方向大于10米，公式表示为：
 
Figure 984777dest_path_image088
式中， s表示机动发生状态，其中1表示发生机动，0表示未发生机动；
Figure 788785dest_path_image089
为前后历元间计算的坐标差异，可表示为： 
Figure 861914dest_path_image090
R 、 N 、 T分别代表轨道坐标系的径向、法向和切向。
优选的，S7中，当检测出机动事件发生时，通过对机动结束后的精密轨道向前积分，重复步骤S4至S6，计算出北斗卫星机动结束时刻。
需要说明的是，上述各技术特征继续相互组合，形成未在上面列举的各种实施例，均视为本发明说明书记载的范围；并且，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (6)

  1. 一种基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法，其特征在于，结合预报的北斗卫星精密轨道与地面台站计算得到的卫星位置做比对，判定北斗卫星在轨运行是否发生机动及异常，并包括以下步骤：
    S1、从事后解算的北斗卫星精密轨道中获取参考时刻卫星状态参数，包括：卫星在J2000.0惯性坐标系下空间位置参数
    Figure 210361dest_path_image001
    ，空间速度参数
    Figure 626299dest_path_image002
    ，其中J2000.0惯性坐标系定义如下，坐标轴原点与地球质心重合，Z轴与J2000.0时刻地球瞬时自转轴重合，X轴指向J2000.0时刻瞬时真春分点，Y轴与Z轴和X轴组成右手坐标系；
    S2、对参考时刻北斗卫星状态参数轨道积分，获取连续弧段的北斗卫星轨道参数及各时刻卫星位置关于初始状态参数的转移矩阵；
    S3、利用步骤S2中的卫星位置参数和状态转移矩阵拟合S1中的状态参数，获取较为准确的状态参数；
    S4、对步骤S3中获的状态参数进行长时间轨道积分，求得积分方向一段时间北斗卫星空间位置参数；
    S5、收集能够跟踪到北斗卫星的地面台站观测数据，利用其精确已知的坐标参数和地面钟差参数，计算北斗卫星空间位置，并将地固坐标系中坐标参数转换到惯性坐标系统中，其中地固系定义如下，坐标系原点与包含海洋和大气层在内的整个地球质心重合，坐标轴指向按照BIH1984.0确定的指向，Z轴与BIH1984.0瞬时自转轴重合，X轴指向BIH1984.0时刻本初子午线与地球赤道交点，Y轴与Z轴和X轴组成右手坐标系；
    S6、对比步骤S4和步骤S5中卫星空间坐标参数，将坐标差异转换到轨道坐标系中，给定先验参数，判别北斗卫星是否发生异常机动，如果发生机动则给出机动发生时刻、机动大小；
    S7、当检测出机动事件发生时，通过对机动结束后的精密轨道向前积分，重复步骤S4至S6，计算出北斗卫星机动结束时刻。
  2. 根据权利要求1所述的基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法，其特征在于，采用Runge-Kutta积分方法6/7阶计算初始10个历元位置，随后采用Adams积分方法进行后续历元连续积分；Runge-Kutta积分方法6阶和7阶公式可表示为：
    Figure 362174dest_path_image003
     
    相应的系数如表所示：
    Figure 905282dest_path_image004
    Adams积分显式公式可表示为：
    Figure 43002dest_path_image005
    Adams积分隐式公式可表示为：
    Figure 629841dest_path_image006
    在轨道数值积分计算时，显式公式与隐式公式同时采用，先由显式公式计算出
    Figure 853012dest_path_image007
    步点上近似值，再由隐式公式校正该近似值，得出需要的
    Figure 199811dest_path_image008
    。
  3. 根据权利要求1中基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法，其特征在于，第S5步骤中，根据全球分布的北斗地面连续跟踪台站计算卫星空间位置的方法，通过固定地面台站地固系坐标和接收机钟差参数，求解北斗卫星空间位置，具体的实现如下：
    利用北斗卫星伪距和相位双频观测值方程如下：
    Figure 457617dest_path_image009
    其中：
    Figure 949778dest_path_image010
    和
    Figure 925824dest_path_image011
    分别表示载波和伪距无电离组合观测值，
    Figure 935369dest_path_image012
    为卫星至测站的几何距离，
    Figure 677976dest_path_image013
    为真空中光速，
    Figure 341038dest_path_image014
    和
    Figure 273222dest_path_image015
    分别表示吸收了接收机端或卫星端硬件延迟的钟差参数，
    Figure 961824dest_path_image016
    为对流层延迟参数，
    Figure 928643dest_path_image017
    为其它需要顾及的改正，
    Figure 28186dest_path_image018
    为无电离层组合相位观测值波长，
    Figure 713245dest_path_image019
    为无电离层组合相位观测值模糊度，可表示为：
    Figure 939958dest_path_image020
    Figure 495704dest_path_image021
    为载波相位和伪距测量值的误差项；观测方法对卫星位置参数求偏导数可表示为：
    Figure 31728dest_path_image022
    其中：
    Figure 938504dest_path_image023
    ，
    Figure 968908dest_path_image024
    ， M为对流层投影函数，
    Figure 910319dest_path_image025
    和
    Figure 617244dest_path_image026
    分别为待估参数近似值和原始观测计算得到的随机误差；经线性化后，待估参数以距离长度为单位，顾及卫星初始状态参数，有：
    Figure 11316dest_path_image027
    其中设计矩阵
    Figure 970045dest_path_image028
    和待估参数可表示为：
    Figure 638399dest_path_image029
    至此，观测方程线性化完成，通过线性估计方法进行待估参数的求解。
  4. 根据权利要求1中基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法，其特征在于，步骤S5中，所述的并将地固坐标系中坐标参数转换到J2000.0惯性坐标系统中，其实现过程如下：
    以 r表示某点在历元J2000.0对应的协议地心惯性系（CIS）中的坐标，以 R表示协议地球固定坐标系（CTS）中的坐标，则有：
    Figure 126012dest_path_image030
    其中：
    Figure 132015dest_path_image031
    为岁差矩阵；
    Figure 628855dest_path_image032
    为章动矩阵；
    Figure 154645dest_path_image033
    为地球自转矩阵；
    Figure 78739dest_path_image034
    为极移矩阵；上述转换矩阵中，
    Figure 40879dest_path_image035
    矩阵根据历元的时刻可以精确确定；
    Figure 951197dest_path_image036
    矩阵可根据1980IAU章动理论来计算；
    Figure 721707dest_path_image037
    矩阵通过从IERS的时间公报中内插求出IAT-UT1进而求得；
    Figure 941336dest_path_image038
    矩阵需要的极移值需从由IERS公布的极移值差值得到。
  5. 根据权利要求1中基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法其特征在于，步骤S6中,所述将坐标差异转换到轨道坐标系中，其实现过程如下：轨道坐标系的坐标轴在惯性系中的单位矢量可以借助于惯性系中的卫星位置
    Figure 531717dest_path_image039
    与速度向量
    Figure 245726dest_path_image040
    来表示，即：
    Figure 870743dest_path_image041
    若轨道坐标系中某一点坐标为
    Figure 261273dest_path_image042
    ，则该点在惯性系中坐标
    Figure 338950dest_path_image043
    为：
    Figure 588141dest_path_image044
    。
  6. 根据权利要求1中基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法，其特征在于，步骤S6中，所述给定先验参数，判别北斗卫星是否发生异常机动，其实现过程如下：
    将计算得到的卫星坐标差异转换到轨道坐标系后，即可直观显示坐标在轨道坐标系下的机动方向和大小，判断卫星机动发生的标准为坐标差异前后历元间差分 R方向大于1米，或 N方向大于5米，或 T方向大于10米，公式表示为：
    Figure 67664dest_path_image045
    式中， s表示机动发生状态，其中1表示发生机动，0表示未发生机动；
    Figure 629096dest_path_image046
    为前后历元间计算的坐标差异，可表示为：
    Figure 459648dest_path_image047
    R 、N 、T分别代表轨道坐标系的径向、法向和切向。
PCT/CN2019/123837 2019-10-30 2019-12-07 基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法 WO2021082189A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201911046908.5 2019-10-30
CN201911046908.5A CN111308515A (zh) 2019-10-30 2019-10-30 基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法

Publications (1)

Publication Number Publication Date
WO2021082189A1 true WO2021082189A1 (zh) 2021-05-06

Family

ID=71145004

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/123837 WO2021082189A1 (zh) 2019-10-30 2019-12-07 基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法

Country Status (2)

Country Link
CN (1) CN111308515A (zh)
WO (1) WO2021082189A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115617023A (zh) * 2022-12-05 2023-01-17 中国西安卫星测控中心 航天器姿控系统异常定位方法和装置

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114114345B (zh) * 2022-01-27 2022-05-13 航天宏图信息技术股份有限公司 一种基于大数据的北斗rnss智能运维方法
CN115883612A (zh) * 2023-02-22 2023-03-31 深圳市北斗云信息技术有限公司 低功耗高精度多参数自适应的北斗物联网模组
CN116698049B (zh) * 2023-08-08 2023-11-07 南京航空航天大学 一种空间非合作目标无源探测初始定轨段的机动检测方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101853027A (zh) * 2010-05-21 2010-10-06 武汉大学 实时精密定轨中轨道的星载快速多步积分方法
US20120166084A1 (en) * 2010-12-23 2012-06-28 Electronics And Telecommunicatons Research Institute Method for fast and precise orbit propogation including maneuver
CN109001776A (zh) * 2018-06-04 2018-12-14 北京未来导航科技有限公司 一种基于云计算的导航数据处理方法及系统
CN109387859A (zh) * 2017-08-14 2019-02-26 千寻位置网络有限公司 基于地面跟踪站产生长期卫星轨道和钟差的方法和装置

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8120529B2 (en) * 2008-09-11 2012-02-21 California Institute Of Technology Method and apparatus for autonomous, in-receiver prediction of GNSS ephemerides
CN106556851B (zh) * 2016-11-25 2017-09-15 中国测绘科学研究院 一种船载gnss辅助北斗导航卫星定轨方法
CN107390233B (zh) * 2017-07-18 2020-04-17 武汉大学 一种低轨卫星导航增强电离层延迟改正参数方法
CN109001972B (zh) * 2018-08-13 2020-06-12 中国科学院国家授时中心 一种北斗广域授时系统与方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101853027A (zh) * 2010-05-21 2010-10-06 武汉大学 实时精密定轨中轨道的星载快速多步积分方法
US20120166084A1 (en) * 2010-12-23 2012-06-28 Electronics And Telecommunicatons Research Institute Method for fast and precise orbit propogation including maneuver
CN109387859A (zh) * 2017-08-14 2019-02-26 千寻位置网络有限公司 基于地面跟踪站产生长期卫星轨道和钟差的方法和装置
CN109001776A (zh) * 2018-06-04 2018-12-14 北京未来导航科技有限公司 一种基于云计算的导航数据处理方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
YAN XING-YUAN , HUANG GUAN-WEN , ZHANG RUI , ZHANG QIN: "A Method Based on Broadcast Ephemeris to Detect BDS Satellite Orbital Maneuver", JOURNAL OF NAVIGATION AND POSITIONING, vol. 3, no. 3, 20 September 2015 (2015-09-20), pages 35 - 38, XP055808550, ISSN: 2095-4999, DOI: 10.16547/j.cnki.10-1096.20150307 *
ZHANG RUI: "Research on Key Technologies of BDS Precise Orbit Determination", CHINESE DOCTORAL DISSERTATIONS FULL-TEXT DATABASE, 27 April 2016 (2016-04-27), pages 1 - 158, XP055808556, ISSN: 1674-022X *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115617023A (zh) * 2022-12-05 2023-01-17 中国西安卫星测控中心 航天器姿控系统异常定位方法和装置

Also Published As

Publication number Publication date
CN111308515A (zh) 2020-06-19

Similar Documents

Publication Publication Date Title
WO2021082189A1 (zh) 基于精密轨道和地面台站数据的北斗卫星机动及异常探测方法
US11726213B2 (en) Fast and precise positioning method and system
CN110275186B (zh) Leo卫星增强的gnss电离层归一化与融合建模方法
CN110275192B (zh) 一种基于智能手机的高精度单点定位方法与装置
CN103033188B (zh) 基于综合孔径观测的导航卫星自主时间同步方法
CN105699999B (zh) 一种固定北斗地基增强系统基准站窄巷模糊度的方法
CN111596321B (zh) 利用非差改正数的多gnss多路径误差恒星日滤波方法及系统
CN111505679B (zh) 一种基于星载gnss的leo初轨确定方法
CN109901206B (zh) 一种基于低轨卫星无线电测距信号的单星定位与授时方法
CN106324629A (zh) 一种bds_gps_glonass融合精密单点定位方法
CN111856534B (zh) 智能终端的双模gnss载波精密单点定位方法及系统
CN110007326B (zh) 一种用于星基增强系统的双频测距误差参数生成方法
CN104865587A (zh) 基于geo星座粗定位的北斗接收机快速定位方法
CN111399015B (zh) 一种适用于船舶交通管理系统的双模卫星融合定位方法
Wu et al. Algorithm theoretical basis document for GRACE level-1B data processing V1. 2
CN113687402B (zh) 一种顾及卫星轨道误差的低轨导航增强实时定位方法
CN114236580A (zh) 一种基于B2b信号的低轨卫星实时定轨与时频同步方法
CN113581501B (zh) 一种适用于组网低轨卫星联合定轨的系统及方法
CN111323796A (zh) 一种gnss接收机高采样钟差解算方法
CN111551975A (zh) Bds/gps参考站低高度角卫星整周模糊度确定方法
CN116299594A (zh) 一种基于PPP-B2b信号的LEO实时定轨与时频同步方法
CN112731504B (zh) 对月球探测器自主定轨的方法及装置
CN115236706B (zh) 一种星间链路单向测距观测数据处理方法及系统
Li et al. Doppler location algorithm and performance analysis of low orbit constellation
CN117782080B (zh) 基于PPP-B2b/INS的实时空基导航系统及方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19950439

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205A DATED 11/03/2022)