CN116381738A - 一种低轨卫星星载时钟钟差高精度预报方法 - Google Patents

一种低轨卫星星载时钟钟差高精度预报方法 Download PDF

Info

Publication number
CN116381738A
CN116381738A CN202310216962.XA CN202310216962A CN116381738A CN 116381738 A CN116381738 A CN 116381738A CN 202310216962 A CN202310216962 A CN 202310216962A CN 116381738 A CN116381738 A CN 116381738A
Authority
CN
China
Prior art keywords
model
clock
polynomial
data
forecasting
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.)
Pending
Application number
CN202310216962.XA
Other languages
English (en)
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.)
Tongji University
Original Assignee
Tongji University
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 Tongji University filed Critical Tongji University
Priority to CN202310216962.XA priority Critical patent/CN116381738A/zh
Publication of CN116381738A publication Critical patent/CN116381738A/zh
Pending legal-status Critical Current

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/27—Acquisition or tracking or demodulation of signals transmitted by the system creating, predicting or correcting ephemeris or almanac data within the receiver
    • 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)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Signal Processing (AREA)
  • Power Engineering (AREA)
  • Radio Relay Systems (AREA)

Abstract

本发明提出一种低轨卫星星载时钟钟差高精度预报方法，步骤包括：对获取的原始数据进行粗差及跳变探测，剔除探测到的异常数据得到干净的钟差；使用拓展最小二乘谐波估计的方法确定低轨卫星钟差的函数模型，并分析模型中谐波项的累积能量占比，确定钟差的主要周期项；附加模型的主要周期项，以预报的精度为指标分析多项式和周期参数独立估计时，不同预报弧长下的最佳拟合弧长和预报精度；根据最佳拟合弧长对低轨卫星的钟差进行高进度的预报该方法率先提出了对特性更加复杂的低轨卫星时钟预报方法，能够准确获取钟差函数模型并全面考虑到模型预报效果与预报模型、拟合方法和拟合弧长的关系，从而有效提升低轨卫星时钟钟差预报精度。

Description

一种低轨卫星星载时钟钟差高精度预报方法
技术领域
本发明涉及到低轨卫星增强GNSS实时高精度应用领域，尤其涉及到一种高精度低轨卫星时钟钟差预报方法。
背景技术
需要时空信息赋能的产业在快速地发展，智慧城市、空间科学研究、国家安全等诸多领域都需要高精度的位置服务，市场和用户给导航定位技术提出了越来越高的需求。全球导航卫星***(Global Navigation Satellite System,GNSS)是常用的导航定位技术，它可以通过实时动态RTK技术和精密单点定位PPP技术两种方式为用户提供高精度的位置服务。RTK技术可以实现快速、高精度的定位，但是需要参考站的支持。PPP技术可以实现全球范围内的高精度定位，但是需要一定时间收敛到较高精度。因此，现有的GNSS技术都无法满足在全球范围内快速高精度位置服务的需求，瞬时厘米级的定位仍然需要参考站的支持。为提升***的服务能力，各大GNSS***都在进行发展，比如GPS的Block III卫星支持3频信号，GLONASS第三代卫星开始支持码分多址等，但是它们基本定位模式没有发生本质的变化，并没有突破全球快速高精度定位的技术瓶颈。
近年来，星链、鸿雁等低轨卫星星座在快速的发展，目前近地轨道上已经有数以千计的低轨卫星为全球提供互联网接入服务。与GNSS的中高轨卫星相比，低轨卫星轨道高度更低，具有快速变化的几何构型以及更强的对地信号强度。一方面可以接收导航信号，成为天基的动态监测监测站。另一方面，也可以向下播发信号，成为近地轨道上的导航卫星。当作为导航卫星时，低轨卫星快速变化的几何构型可以大幅度缩短PPP的收敛时间，同时大量的低轨卫星也能为城市峡谷、城市高架等观测受限区域提供更多的观测值支撑。
对于实时高精度低轨增强应用，特别是实时PPP而言，高精度实时钟差产品至关重要。为了满足实时用户的需求，IGS从13年就开始通过RTS服务(Real-time Service)为用户提供实时精密轨道钟差改正产品，但是在RTS服务中断或是不可用的情况下，就只能利用预报模型对钟差进行预报。对于GNSS的钟差预报已经有大量研究，其中附有周期的多项式模型能够很好地描述钟差的物理和周期特性，具有较好的预报效果并被广泛应用；而对于低轨卫星钟差而言，相关研究比较匮乏。此外，低轨卫星受限于成本，所搭载的卫星钟与GNSS的原子钟相比稳定性更低，钟差特性也更为复杂。
发明内容
本发明的目的在于，针对目前低轨卫星增强GNSS实时应用时RTS服务中断的痛点，提供一种高精度的低轨卫星钟差预报方法，用于提高低轨增强GNSS实时应用在RTS服务缺乏情况下的可用性。
为了达到上述目的，本发明技术方案如下：
一种低轨卫星钟差高精度预报方法，包括如下步骤：
步骤S1：对获取的原始数据进行粗差及跳变探测，剔除探测到的异常数据得到干净的钟差；
步骤S2：使用拓展最小二乘谐波估计的方法确定低轨卫星钟差的函数模型，并分析模型中谐波项的累积能量占比，确定钟差的主要周期项；
步骤S3：附加模型的主要周期项，以预报的精度为指标分析多项式和周期参数独立估计时，不同预报弧长下的最佳拟合弧长和预报精度；
步骤S4：根据S3确定的最佳拟合弧长对低轨卫星的钟差进行高进度的预报。
综上所述，本发明一种低轨卫星星载时钟钟差高精度预报方法，能够有效避免数据异常或是缺失对数据分析的影响，获取精确的低轨卫星钟差函数模型，在附加主周期的钟差模型基础上充分考虑不同参数对于拟合弧长要求的差异，实现高精度的低轨卫星钟差预报。
具体的，与现有技术相比，本发明具有以下优点：
现有的钟差预报方法只针对于GNSS的卫星，对于星载钟稳定性更差、特性更加复杂的低轨卫星钟差预报缺少相应的钟差预报方法。现有的钟差周期多项式函数模型的确定只考虑到了模型的周期项部分，缺乏对模型多项式部分的确定，从到导致周期探测不准确，影响模型的预报效果。现用周期探测方法对于低轨卫星钟差中的异常数据敏感，并且分辨率较低，对低轨卫星钟差的周期分析效果较差。本发明借助了最小二乘谐波估计进行函数模型确定的优越性，使用周期累计能量分析，精确确定低轨卫星钟差的函数模型，同时充分考虑了周期项和多项式参数拟合对拟合弧长要求的差异，能够实现高精度的低轨卫星钟差预报，填补了低轨钟差预报的空白。
附图说明
图1为低轨卫星星载时钟钟差高精度预报方法的流程示意图。
具体实施方式
下面将结合示意图对本发明的具体实施方式进行更详细的描述。根据下列描述本发明的优点和特征将更清楚。
参考图1，本发明一优选实施例中，一种低轨卫星星载时钟钟差高精度预报方法包括：
步骤S1：对获取的原始数据进行粗差及跳变探测，剔除探测到的异常数据得到干净的钟差；
步骤S2：使用拓展最小二乘谐波估计的方法确定低轨卫星钟差的函数模型，并分析模型中谐波项的累积能量占比，确定钟差的主要周期项；
步骤S3：附加模型的主要周期项，以预报的精度为指标分析多项式和周期参数独立估计时，不同预报弧长下的最佳拟合弧长和预报精度；
步骤S4：根据S3确定的最佳拟合弧长对低轨卫星的钟差进行高进度的预报。
所述步骤S1包括：
步骤S11：将原始钟差数据由相位数据转换为频率数据；
原始钟差序列中的异常数据，如粗差及跳变等在序列中的特征不明显，需要转换为频率数据才能进行识别和定位。对于原始钟差序列x(t)，其相位数据定义为：
φ(t)＝2πf0x(t) (1)
其中φ为相位数据，t为时间，π为圆周率，f0为时钟的标称频率。其对应的频率数据可以由相位数据进行转换，公式为：
Figure BDA0004115261020000031
其中
Figure BDA0004115261020000032
为频率数据的均值，τ为时间间隔。使用公式(2)将原始的钟差数据转换为频率数据，提供给S12进行异常数据的探测。
步骤S12：使用绝对中位差(MAD)的方法对频率数据中的粗差和跳变进行识别；
MAD法相较于中误差法避免了反复求取中误差的复杂计算，具有较好的异常数据识别能力的同时效率也更高，MAD的定义如下：
MAD＝Median{|yi-med|/0.6745} (3)
公式中，Median代表取中位数的操作，i为标注的自然数，yi为第i个历元的频率数据，med为频率数据序列中的中位数，即med＝Meidan{yi},系数0.6745能使得MAD的值约等于正态分布数据序列中的标准差，异常数据的判断条件为：
|yi-med|>n·MAD (4)
对于频率序列中的每一个频率，如果满足公式(4)，则认为是异常数据进行标记，公式中的n为自然数，一般取3到5间的整数
步骤S13：对识别出的异常数据进行剔除后，将频率数据还原为相位数据。
所述步骤S2包括：
步骤S21：构建最小二乘谐波估计的原始假设和备选假设；
周期多项式模型由多项式部分和周期项部分组成，对一个由周期多项式模型描述的时间序列，其函数表达式如下：
Figure BDA0004115261020000041
其中E(*)代表数学期望，y代表钟差序列，D(*)代表数据的方差，Q代表方差矩阵，钟差y的方差为Qy。其中A为模型的多项式设计矩阵，Ak为周期项的设计矩阵，k,q为标识的自然数，x,xk为对应的参数向量，周期项由三角函数表示，则
Figure BDA0004115261020000042
式中，ωk为Ak对应的周期项的角速度，ak,bk为对应的参数。为了找到公式(5)中最显著的多项式阶数或是最显著的周期项，基于公式(5)提出原始建设和备选假设：
H0:E(y)＝A0x0,D(y)＝Qy (7)
与
Ha:E(y)＝A0x0+Aqxq,D(y)＝Qy (8)
其中A0＝[AA1…Aq-1]，为模型目前的设计矩阵，Aq为在A0模型的基础上附加的部分，在周期多项式中为某阶多项式或是某个频率的设计矩阵。
Figure BDA0004115261020000043
为估值，/>
Figure BDA0004115261020000044
和/>
Figure BDA0004115261020000045
分别为原始假设和备选假设的残差向量，对于备选假设Ha而言，最显著的附加模型满足/>
Figure BDA0004115261020000046
在Qy下的模最小，以频率ωq为例，则满足
Figure BDA0004115261020000047
其中，
Figure BDA0004115261020000048
代表残差向量在/>
Figure BDA0004115261020000049
空间下的模，为对(9)进行简化，对于(8)而言，其法方程可以表示为：
Figure BDA00041152610200000410
对(10)进行分块求解，可以得到
Figure BDA00041152610200000411
和/>
Figure BDA00041152610200000412
的分块解析式
Figure BDA00041152610200000413
其中
Figure BDA0004115261020000051
为矩阵操作，具体的为/>
Figure BDA0004115261020000052
Figure BDA0004115261020000053
其中I为单位矩阵。基于(11)可以得到钟差估值/>
Figure BDA0004115261020000054
Figure BDA0004115261020000055
基于公式(12)可以得到残差向量e关于设计矩阵A0和Aq的表达式，这可以简化附加模型的探测定位
Figure BDA0004115261020000056
公式(13)给出了一个重要的结论，当模型附加额外参数时，模型对于原始数据拟合的残差会减小，减小的数值为
Figure BDA0004115261020000057
与原有设计矩阵和附加模型的设计矩阵有关。考虑到钟差数据间的协方差为0(即Qy＝I)，基于(13)，公式(9)的简化表达为
Figure BDA0004115261020000058
公式(14)为周期项的探测公式，对于多项式部分，则需要将周期项的设计矩阵Aj替换为k阶多项式的设计矩阵Bk，则对应的探测公式为：
Figure BDA0004115261020000059
其中
Figure BDA00041152610200000510
与上面的进行的/>
Figure BDA00041152610200000511
矩阵操作一致，但/>
Figure BDA00041152610200000512
l这里为多项式的阶数，一般小于4的自然数。
步骤S22：构造假设检验的统计量；
由(13)可知，附加额外的参数都可以使得模型的拟合效果更好，但随着模型参数的增加，模型会过拟合，导致模型冗余甚至预报效果变差。因此需要设置截至条件防止模型过拟合。根据广义似然比检验(GLR-test)的原理，当如下公式成立时，会对原假设H0进行拒绝：
Figure BDA00041152610200000513
其中，σ2为数据的方差，a为0到1间的实数，用于给定检验的严苛程度，左式分号上下为似然函数，原始假设和备选假设的似然函数为
Figure BDA00041152610200000514
Figure BDA00041152610200000515
e为自然常数，det()为矩阵的取行列式操作，Ry为y的权逆阵，即Qy＝σ2Ry，m为自然数，这里表示数据的数量。通过最大化似然函数，(16)能够推导为如下的形式
Figure BDA0004115261020000061
根据之前的定义，m,n,q分别为数据个数、多项式阶数与周期项个数。当只考虑白噪声时，原始假设和备选假设的统计量可以基于(11)和(19)构造为
Figure BDA0004115261020000062
其中，
Figure BDA0004115261020000063
是备选假设下的方差估值。由于/>
Figure BDA0004115261020000064
Figure BDA0004115261020000065
统计量T2在原始假设下满足中心Fisher分布，即
T2～F(2,m-n-2q) (21)
对于每次探测到的附加模型，需要构建统计量进行(21)所给出的假设检验，如果通过检验则认为备选假设下附加了额外的参数后模型已经过拟合。
S23：进行模型的多项式识别；
基于一次多项式模型，使用S21与S22所构建的拓展最小二乘谐波估计进行模型多项式阶数的识别。每轮探测中分别探测满足(14)和(15)的周期项和多项式阶数，其中数值更大的认为是原始假设下最为显著的附加模型；进一步进行假设检验，不通过则将附加的参数整合到模型中进行下一轮的探测，直到检验通过为止。完成全部的探测后，保留模型的多项式部分。
S24：进行模型的周期项部分识别；
基于S23确定的多项式模型，重复S23中的步骤，不同的是本次探测中认为多项式模型已经正确确定，只进行周期项的探测。当检验通过时，认为确定了钟差的函数模型。
S25：确定模型的主要周期项；
对于钟差而言，往往几个主要的周期项就包括了周期信号中的大部分能量，能够有效拟合钟差数据，同时也能避免振幅较小的周期估计不准对预报的影响，进而实现高精度的钟差预报。对确定的钟差函数模型，进一步以能量分析的方式来确定模型中的主要周期项。定义周期部分的能量为：
Figure BDA0004115261020000066
(22)式表示频率为ωk的周期项在整个序列上的能量。计算了所有周期项的能量后，按能量大小计算周期项的累计能量，当累计能量占比大于整个周期信号总能量的80％以上时，则认为确定了钟差的主周期。
所述步骤S3包括：
步骤S31：根据S2确定的函数模型，构建钟差预报模型；
由S2中确定的多项式阶数和主周期，钟差的函数模型可以表示为
Figure BDA0004115261020000071
其中n为S2确定的周期项阶数，ci为多项式系数，q为确定的主周期个数，j为标注的自然数，ωj,aj,bj为周期参数，ε为误差。为了减小多项式常数项c0估计不准对预报的影响，对(23)进行一阶差分来消除常数项，则钟差的一阶差可以表示为
Figure BDA0004115261020000072
当完成了参数aj,bj,cj的估计后，则可以构建如下的钟差预报模型
Figure BDA0004115261020000073
其中tlast为进行参数拟合的最后一个历元。
步骤S32：使用不同拟合弧度进行参数拟合。
当多项式参数和周期参数分别估计时，模型具有更好的预报效果。首先对周期参数进行估计，使用24小时的弧长拟合模型，获得aj,bj参数的估值。对于多项式参数需要确定最佳的拟合弧长，在固定周期参数的前提下从30分钟开始按30秒的步长增加拟合弧长到24小时，估计每一个弧长下的多项式参数。
步骤S33：使用S32得到的模型参数进行不同预报弧长的钟差预报。
步骤S34：分析不同预报弧长下最佳拟合弧长。
计算S25所确定的预报模型的预报精度，预报的精度定义为预报钟差与真实钟差之间的均方根误差RMSE
Figure BDA0004115261020000074
其中yp为模型的预报值，y为真实的钟差，N为预报的钟差个数。为了使得结果更加准确，分析的结果为一段时间内(如几天)的统计平均值。对于某一个预报弧长下，其所有拟合弧长的预报结果中预报精度最高的拟合弧长则为该预报弧长的最佳拟合弧长。确定所有预报弧长的最佳拟合弧长。
所述步骤S4包括：
步骤S41：对每一个需要预报的弧长，使用S3确定的最佳拟合弧长拟合预报模型。
步骤S42：使用S41拟合的预报模型对钟差进行预报。
综上所述，本发明一种低轨卫星星载时钟钟差高精度预报方法，能够有效避免数据异常或是缺失对数据分析的影响，获取精确的低轨卫星钟差函数模型，在附加主周期的钟差模型基础上充分考虑不同参数对于拟合弧长要求的差异，实现高精度的低轨卫星钟差预报。
具体的，与现有技术相比，本发明具有以下优点：
现有的钟差预报方法只针对于GNSS的卫星，对于星载钟稳定性更差、特性更加复杂的低轨卫星钟差预报缺少相应的钟差预报方法。现有的钟差周期多项式函数模型的确定只考虑到了模型的周期项部分，缺乏对模型多项式部分的确定，从到导致周期探测不准确，影响模型的预报效果。现用周期探测方法对于低轨卫星钟差中的异常数据敏感，并且分辨率较低，对低轨卫星钟差的周期分析效果较差。本发明借助了最小二乘谐波估计进行函数模型确定的优越性，使用周期累计能量分析，精确确定低轨卫星钟差的函数模型，同时充分考虑了周期项和多项式参数拟合对拟合弧长要求的差异，能够实现高精度的低轨卫星钟差预报，填补了低轨钟差预报的空白。
上述仅为本发明的优选实施例而已，并不对本发明起到任何限制作用。任何所属技术领域的技术人员，在不脱离本发明的技术方案的范围内，对本发明揭露的技术方案和技术内容做任何形式的等同替换或修改等变动，均属未脱离本发明的技术方案的内容，仍属于本发明的保护范围之内。

Claims (5)

1.一种低轨卫星星载时钟钟差高精度预报方法，其特征在于，包括如下步骤：
步骤S1：对获取的原始数据进行粗差及跳变探测，剔除探测到的异常数据得到干净的钟差；
步骤S2：使用拓展最小二乘谐波估计的方法确定低轨卫星钟差的函数模型，并分析模型中谐波项的累积能量占比，确定钟差的主要周期项；
步骤S3：附加模型的主要周期项，以预报的精度为指标分析多项式和周期参数独立估计时，不同预报弧长下的最佳拟合弧长和预报精度；
步骤S4：根据S3确定的最佳拟合弧长对低轨卫星的钟差进行高进度的预报。
2.如权利要求1所述的低轨卫星星载时钟钟差高精度预报方法，其特征在于，所述步骤S1包括：
步骤S11：将原始钟差数据由相位数据转换为频率数据；
原始钟差序列中的异常数据，如粗差及跳变等在序列中的特征不明显，需要转换为频率数据才能进行识别和定位；对于原始钟差序列x(t)，其相位数据定义为：
φ(t)＝2πf0x(t) (1)
其中φ为相位数据，t为时间，π为圆周率，f0为时钟的标称频率；其对应的频率数据可以由相位数据进行转换，公式为：
Figure FDA0004115260990000011
其中
Figure FDA0004115260990000012
为频率数据均值，τ为时间间隔；使用公式(2)将原始的钟差数据转换为频率数据，提供给S12进行异常数据的探测；
步骤S12：使用绝对中位差MAD的方法对频率数据中的粗差和跳变进行识别；
MAD法相较于中误差法避免了反复求取中误差的复杂计算，MAD的定义如下：
MAD＝Median{|yi-med|/0.6745} (3)
公式中，Median代表取中位数的操作，i为标注的自然数，yi为第i个历元的频率数据，med为频率数据序列中的中位数，即med＝Median{yi}，系数0.6745能使得MAD的值约等于正态分布数据序列中的标准差，异常数据的判断条件为：
|yi-med|＞n·MAD (4)
对于频率序列中的每一个频率，如果满足公式(4)，则认为是异常数据进行标记，公式中的n为自然数，一般取3到5间的整数；
步骤S13：对识别出的异常数据进行剔除后，将频率数据还原为相位数据。
3.如权利要求1所述的低轨卫星星载时钟钟差高精度预报方法，其特征在于，所述步骤S2包括：
步骤S21：构建最小二乘谐波估计的原始假设和备选假设；
周期多项式模型由多项式部分和周期项部分组成，对一个由周期多项式模型描述的时间序列，其函数表达式如下：
Figure FDA0004115260990000021
其中E(*)代表数学期望，y代表钟差序列，D(*)代表数据的方差，Q代表方差矩阵，钟差y的方差为Qy；其中A为模型的多项式设计矩阵，Ak为周期项的设计矩阵，k,q为标识的自然数，x,xk为对应的参数向量，周期项由三角函数表示，则
Figure FDA0004115260990000022
式中，ωk为Ak对应的周期项的角速度，ak,bk为对应的参数；为了找到公式(5)中最显著的多项式阶数或是最显著的周期项，基于公式(5)提出原始建设和备选假设：
H0:E(y)＝A0x0,D(y)＝Qy (7)
与
Ha:E(y)＝A0x0+Aqxq,D(y)＝Qy (8)
其中A0＝[AA1…Aq-1]，为模型目前的设计矩阵，Aq为在A0模型的基础上附加的部分，在周期多项式中为某阶多项式或是某个频率的设计矩阵；
Figure FDA0004115260990000023
为估值，/>
Figure FDA0004115260990000024
和/>
Figure FDA0004115260990000025
分别为原始假设和备选假设的残差向量，对于备选假设Ha而言，最显著的附加模型满足/>
Figure FDA0004115260990000026
在Qy下的模最小，以频率ωq为例，则满足
Figure FDA0004115260990000027
其中，
Figure FDA0004115260990000028
代表残差向量在/>
Figure FDA0004115260990000029
空间下的模，为对(9)进行简化，对于(8)而言，其法方程表示为：
Figure FDA00041152609900000210
对(10)进行分块求解，得到
Figure FDA00041152609900000211
和/>
Figure FDA00041152609900000212
的分块解析式
Figure FDA00041152609900000213
其中
Figure FDA0004115260990000031
为矩阵操作，具体的为/>
Figure FDA0004115260990000032
Figure FDA0004115260990000033
其中I为单位矩阵；基于(11)得到钟差估值/>
Figure FDA0004115260990000034
Figure FDA0004115260990000035
基于公式(12)得到残差向量e关于设计矩阵A0和Aq的表达式，这简化附加模型的探测定位
Figure FDA0004115260990000036
公式(13)给出了一个重要的结论，当模型附加额外参数时，模型对于原始数据拟合的残差会减小，减小的数值为
Figure FDA0004115260990000037
与原有设计矩阵和附加模型的设计矩阵有关；考虑到钟差数据间的协方差为0(即Qy＝I)，基于(13)，公式(9)的简化表达为
Figure FDA0004115260990000038
公式(14)为周期项的探测公式，对于多项式部分，则需要将周期项的设计矩阵Aj替换为k阶多项式的设计矩阵Bk，则对应的探测公式为：
Figure FDA0004115260990000039
其中
Figure FDA00041152609900000310
与上面的进行的/>
Figure FDA00041152609900000311
矩阵操作一致，但/>
Figure FDA00041152609900000312
l为多项式的阶数，为小于4的自然数；
步骤S22：构造假设检验的统计量；
设置截至条件防止模型过拟合；根据广义似然比检验(GLR-test)的原理，当下式成立时，会对原假设H0进行拒绝：
Figure FDA00041152609900000313
其中，σ2为数据的方差，a为0到1间的实数，用于给定检验的严苛程度，左式分号上下为似然函数，原始假设和备选假设的似然函数为
Figure FDA00041152609900000314
Figure FDA00041152609900000315
e为自然常数，det()为矩阵的取行列式操作，Ry为y的权逆阵，即Qy＝σ2Ry，m为自然数，这里表示数据的数量；通过最大化似然函数，(16)能够推导为如下的形式
Figure FDA0004115260990000041
根据之前的定义，m,n,q分别为数据个数、多项式阶数与周期项个数；当只考虑白噪声时，原始假设和备选假设的统计量T2基于(11)和(19)构造为
Figure FDA0004115260990000042
其中，
Figure FDA0004115260990000046
是备选假设下的方差估值；由于/>
Figure FDA0004115260990000043
Figure FDA0004115260990000044
统计量T2在原始假设下满足中心Fisher分布，即
T2～F(2,m-n-2q) (21)
对于每次探测到的附加模型，需要构建统计量进行(21)所给出的假设检验，如果通过检验则认为备选假设下附加了额外的参数后模型已经过拟合；
S23：进行模型的多项式识别；
基于一次多项式模型，使用S21与S22所构建的拓展最小二乘谐波估计进行模型多项式阶数的识别；每轮探测中分别探测满足(14)和(15)的周期项和多项式阶数，其中数值更大的认为是原始假设下最为显著的附加模型；进一步进行假设检验，不通过则将附加的参数整合到模型中进行下一轮的探测，直到检验通过为止；完成全部的探测后，保留模型的多项式部分；
S24：进行模型的周期项部分识别；
基于S23确定的多项式模型，重复S23中的步骤，不同的是本次探测中认为多项式模型已经正确确定，只进行周期项的探测；当检验通过时，认为确定了钟差的函数模型；
S25：确定模型的主要周期项；
对于钟差而言，往往几个主要的周期项就包括了周期信号中的大部分能量，能够有效拟合钟差数据，同时也能避免振幅较小的周期估计不准对预报的影响，进而实现高精度的钟差预报；对确定的钟差函数模型，进一步以能量分析的方式来确定模型中的主要周期项；定义周期部分的能量为：
Figure FDA0004115260990000045
(22)式表示频率为ωk的周期项在整个序列上的能量；计算了所有周期项的能量后，按能量大小计算周期项的累计能量，当累计能量占比大于整个周期信号总能量的80％以上时，则认为确定了钟差的主周期。
4.如权利要求1所述的低轨卫星星载时钟钟差高精度预报方法，其特征在于，所述步骤S3包括：
步骤S31：根据S2确定的函数模型，构建钟差预报模型；
由S2中确定的多项式阶数和主周期，钟差的函数模型可以表示为
Figure FDA0004115260990000051
其中n为S2确定的周期项阶数，ci为多项式系数，q为确定的主周期个数，j为标注的自然数，ωj,aj,bj为周期参数，ε为误差；为了减小多项式常数项c0估计不准对预报的影响，对(23)进行一阶差分来消除常数项，则钟差的一阶差可以表示为
Figure FDA0004115260990000052
当完成了参数aj,bj,cj的估计后，则可以构建如下的钟差预报模型
Figure FDA0004115260990000053
其中tlast为进行参数拟合的最后一个历元；
步骤S32：使用不同拟合弧度进行参数拟合；
当多项式参数和周期参数分别估计时，模型具有更好的预报效果；首先对周期参数进行估计，使用24小时的弧长拟合模型，获得aj,bj参数的估值；对于多项式参数需要确定最佳的拟合弧长，在固定周期参数的前提下从30分钟开始按30秒的步长增加拟合弧长到24小时，估计每一个弧长下的多项式参数；
步骤S33：使用S32得到的模型参数进行不同预报弧长的钟差预报；
步骤S34：分析不同预报弧长下最佳拟合弧长；
计算S25所确定的预报模型的预报精度，预报的精度定义为预报钟差与真实钟差之间的均方根误差RMSE
Figure FDA0004115260990000054
其中yp为模型的预报值，y为真实的钟差，N为预报的钟差个数；为了使得结果更加准确，分析的结果为一段时间内的统计平均值；对于某一个预报弧长下，其所有拟合弧长的预报结果中预报精度最高的拟合弧长则为该预报弧长的最佳拟合弧长；确定所有预报弧长的最佳拟合弧长。
5.如权利要求1所述的低轨卫星星载时钟钟差高精度预报方法，其特征在于，所述步骤S4包括：
步骤S41：对每一个需要预报的弧长，使用S3确定的最佳拟合弧长拟合预报模型；
步骤S42：使用S41拟合的预报模型对钟差进行预报。
CN202310216962.XA 2023-03-07 2023-03-07 一种低轨卫星星载时钟钟差高精度预报方法 Pending CN116381738A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310216962.XA CN116381738A (zh) 2023-03-07 2023-03-07 一种低轨卫星星载时钟钟差高精度预报方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310216962.XA CN116381738A (zh) 2023-03-07 2023-03-07 一种低轨卫星星载时钟钟差高精度预报方法

Publications (1)

Publication Number Publication Date
CN116381738A true CN116381738A (zh) 2023-07-04

Family

ID=86970280

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310216962.XA Pending CN116381738A (zh) 2023-03-07 2023-03-07 一种低轨卫星星载时钟钟差高精度预报方法

Country Status (1)

Country Link
CN (1) CN116381738A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117111110A (zh) * 2023-07-11 2023-11-24 武汉纺织大学 一种卫星钟差数据短期预报方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117111110A (zh) * 2023-07-11 2023-11-24 武汉纺织大学 一种卫星钟差数据短期预报方法
CN117111110B (zh) * 2023-07-11 2024-03-08 武汉纺织大学 一种卫星钟差数据短期预报方法

Similar Documents

Publication Publication Date Title
CN107356947B (zh) 基于单频导航卫星数据确定卫星差分伪距偏差的方法
CN108508461B (zh) 基于gnss载波相位高精度定位完好性监测方法
Farrell et al. Differential GPS reference station algorithm-design and analysis
RU2591953C2 (ru) Навигационная система и способ разрешения целочисленных неоднозначностей с использованием ограничения неоднозначности двойной разности
EP2746811B1 (en) Methods for generating accuracy information on an ionosphere model for satellite navigation applications
Ge et al. A new data processing strategy for huge GNSS global networks
CN106249256B (zh) 基于粒子群优化算法的实时glonass相位偏差估计方法
CN109709579B (zh) 一种基于用户测距误差实时估计的gnss卫星星历故障检测方法
CN108076662A (zh) 具有使用未组合公式来解算模糊度的能力的gnss接收机
CN109085628A (zh) 一种整周模糊度的固定方法及***
Krypiak-Gregorczyk et al. Carrier phase bias estimation of geometry-free linear combination of GNSS signals for ionospheric TEC modeling
CN113138402B (zh) 基于rtk的模糊度固定方法及装置、存储介质
CN116381738A (zh) 一种低轨卫星星载时钟钟差高精度预报方法
CN107505642A (zh) 一种ins辅助的实时bds单频周跳探测方法
CN115792980A (zh) 一种模型和数据双重驱动的gnss rtk定位选星方法及***
CN114895336B (zh) 一种gnss实时动态定位***中参考站观测值的预报方法
CN106371092A (zh) 一种基于gps与强震仪观测自适应组合的形变监测方法
CN110068848B (zh) 一种高性能rtk处理技术方法
CN111522032A (zh) 一种北斗三号***用户完好性处理的优化方法及优化装置
CN111194001A (zh) Lte指纹定位校正的方法、装置及***
CN113031036B (zh) 基于GNSS 30s采样频率数据的电离层相位闪烁因子构建方法
El-Rabbany et al. Effect of temporal physical correlation on accuracy estimation in GPS relative positioning
CN114355410B (zh) 基于并行计算的卫星导航实时精密单点定位***及方法
CN116010770A (zh) 一种基于网络应用的时钟钟差估计预测方法
CN115755131A (zh) 一种卫星定位的方法、装置及介质

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