CN109738919A - 一种用于gps接收机自主预测星历的方法 - Google Patents
一种用于gps接收机自主预测星历的方法 Download PDFInfo
- Publication number
- CN109738919A CN109738919A CN201910150844.7A CN201910150844A CN109738919A CN 109738919 A CN109738919 A CN 109738919A CN 201910150844 A CN201910150844 A CN 201910150844A CN 109738919 A CN109738919 A CN 109738919A
- Authority
- CN
- China
- Prior art keywords
- satellite
- formula
- partial derivative
- ephemeris
- orbit
- 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.)
- Granted
Links
Landscapes
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开了一种用于GPS接收机自主预测星历的方法,包括:步骤1:根据输入的广播星历计算星历有效时间内卫星的位置速度矢量,作为精密定轨的观测值;步骤2:根据观测值进行精密定轨,利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量,以及太阳光压模型参数中的卫星面质比;步骤3:进行轨道预报,以定轨弧段初始时刻的的精确位置速度为起始,采用卫星面质比作为卫星等效面质比参数进行轨道外推,外推到开机时刻;步骤4:将开机时刻的位置速度量拟合成广播星历,给出预报星历,将预报星历直接应用于接收机位置的定位解算。极大地缩小了非热启动时的接收机首次定位时长,降低对芯片存储空间的要求,有利于芯片的成本控制。
Description
技术领域
本发明属于GPS接收机领域,涉及一种用于GPS接收机自主预测星历的方法。
背景技术
随着汽车及电子计算机技术的发展和人们生活需求的提高,全球卫星定位系统(GPS)得到了广泛的应用,越来越多的人开始使用GPS接收机和导航系统。GPS接收机在车载导航和移动消费电子等领域迅速发展。随着应用领域的不断扩大和使用群体的迅速增加,使用者对GPS接收机的功能和性能提出了越来越多的要求,如灵敏度(Sensitivity)和定位精度(Accuracy)等。此外,开机后首次定位的时间(TTFF,Time To First Fix)也是其中十分重要的一项性能指标。使用者希望在GPS接收机开机后很短的时间便可以定位,目前的GPS接收机在进行热启动时,首次定位时长为1秒钟。但是,当关机时长超过热启动有效期2小时时,首次定位时长约30秒。目前采用的方法为,接收机利用关机前的广播星历信息自动预测GNSS星历,在接收机重新开机时,利用预测GNSS星历加速GNSS信号捕获。
但是,目前的现有方法对首次定位时长改进甚微。接收机从开机后到完成首次定位需经历捕获、跟踪、收集广播星历几个阶段,其中收集广播星历阶段耗时最长约30秒,而在开阔天空下捕获耗时仅需不到1秒。目前虽能加快GNSS信号捕获,但并不能显著改善首次定位时长。占用大量存储器空间,将预测轨道以地心地固坐标系ECEF的形式存储在非易失性存储其中,随着预测时长的增长,所需空间线性增长。
发明内容
本发明的目的在于克服上述现有技术的缺点,提供一种用于GPS接收机自主预测星历的方法。
为达到上述目的,本发明采用以下技术方案予以实现:
一种用于GPS接收机自主预测星历的方法,包括以下步骤:
步骤1:根据输入的广播星历计算星历有效时间内卫星的位置速度矢量,作为精密定轨的观测值;
步骤2:根据所述观测值进行精密定轨,利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量,以及太阳光压模型参数中的卫星面质比;
步骤3:进行轨道预报,以所述定轨弧段初始时刻的精确位置速度为起始,采用所述卫星面质比作为卫星等效面质比参数进行轨道外推,外推到开机时刻,得到开机时刻的位置速度矢量;
步骤4:将所述开机时刻的位置速度矢量拟合成广播星历,直接应用于接收机的定位解算。
本发明进一步的改进在于:
步骤2中进行精密定轨的具体方法为:
建立航天器的状态量模型:
其中:r,v分别表示t时刻卫星的位置矢量和速度矢量,β表示待估参数卫星等效面质比;
观测量与航天器的状态量X的量测方程:
Y=Y(X,t) (2)
其中:γ表示观测量的测量值且只包含位置矢量;X(t)满足轨道动力学方程:
求解式(3)得到状态方程:
X(t)=X(t0,X0;t) (4)
其中:X0为待估状态量;
将式(2)在参考状态量并舍去二阶和二阶以上的小量,得到:
其中:参考状态量为待估状态量X0的近似值;Y(X*,t)是观测量的理论计算值,是测量矩阵,记为是状态转移矩阵,记为Φ(t0,t);
令:
得到精密定轨基本方程:
y=Bx0 (8)
其中:y为残差;
通过迭代方式由观测采样数据求解精密定轨基本方程(8),得到待估状态量X0的改正值给出历元t0时刻达到精度要求的状态量则:
其中:下标k表示引用了k次采样数据;重复迭代i次直至满足精度要求,得到:
步骤2中:
精密定轨基本方程(8)的求解采用最小二乘估计,通过式(11)计算得到:
其中:表示从t0到tl时刻的状态转移矩阵,Hl表示tl时刻的测量矩阵,yl表示tl时刻的残差向量;
矩阵B通过式(12)计算:
测量矩阵HX通过式(13)计算:
其中:x,y,z表示三个位置分量;
状态转移矩阵Φ(t0,t)通过状态方程(4)得到,状态转移矩阵Φ(t0,t)通过式(14)计算:
步骤2中得到状态转移矩阵Φ(t0,t)的具体方法为:
令:
X(t)=Φ(t,t0)X0 (15)
则:
状态转移矩阵Φ(t0,t)所满足的条件通过对状态微分方程(3)的线性化得到,有:
其中:状态微分方程F为:
记:
则方程(17)的系数矩阵A:
将方程(17)转换为一阶方程组:
其中:C1(t),C2(t),C3(t)均用参考值X*进行计算;
对式(21)进行积分得到通过得到状态转移矩阵Φ(t0,t)。
步骤3的具体方法为:
将轨道预报问题描述为一阶常微分方程初值问题:
其中:卫星的加速度且满足r为卫星位置,为卫星速度;f(t,y)为重力场、日月引力和太阳光压引起的卫星加速度;
采用Adams-Cowell积分法求解式(22)进行轨道预报,并采用8阶龙哥库塔方法作为Adams-Cowell积分法的起步单步法。
步骤3中重力场引起的卫星加速度采用引力加速度估计函数法计算,重力场的偏导数计算方法:
6-1:将重力场分为地球质心引力部分和地球带谐摄动加速度、田谐摄动部分两部分;
6-2:地球质心引力加速度偏导数计算:
地球中心引力加速度表示如下:
其中,r=[x y z]T表示卫星在地固系下的位置矢量,r=|r|;
对式(23)求偏导数得到地球质心引力加速度对位置的偏导数:
地球质心引力加速度对速度的偏导数:
6-3:地球带谐摄动加速度、田谐摄动加速度的偏导数计算:
地球重力势:
其中:R表示地球半径,Cnm和Snm表示未归一化的重力势系数,根据式(27)把归一化后的重力势系数和转化成未归一化前的Cnm和Snm:
Vnm和Wnm的计算公式如下:
已知:
利用式(29)置m=0获得带谐项Vn0,相应的值均等于0,递归使用式(31):
式(31)中过程采用式(28),↓过程采用式(29),计算得到所有的Vnm和Wnm,
加速度等于势能u的梯度,直接由Vnm和Wnm计算得出:
其中:
地球带谐摄动加速度、田谐摄动加速度对位置的偏导数根据:
通过式(40)计算得到:
地球带谐摄动加速度、田谐摄动加速度对速度的偏导数:
步骤3中日月引力引起的卫星加速度的具体确定方法为:
在以地球为中心的参考系下,太阳和月球的摄动表示为:
其中:r表示卫星位置,s表示太阳或月亮的位置;
低精度日月坐标具体计算方法如下:
T=(MJD-51544.5)/36525.0 (43)
摄动的计算基于月球平黄经L0、月球平近点角l、太阳平近点角l′、太阳平黄经和月球平黄经之间的差D以及月球平升交距角F,5个基础角;
L0=(0.606433+1336.851344T-floor(0.606433+1336.851344T)) (44)
记相对于2000年黄道和春分点的月球黄经为LM,利用式(45)计算值,有:
LM=2π(L0+dL/1296.0e3-floor(L0+dL/1296.0e3)) (47)
月球黄纬BM表示为:
月球的地心距rM表示为:
月球坐标rM为:
其中:ε为黄赤交角,ε=23°.43929111;
假设地球围绕太阳作无摄运动,通过式(51)描述太阳相对于地球和黄道在2000年附近几十年的椭圆轨道:
位置坐标利用式(51)由开普勒轨道公式得到;
太阳的黄经L⊙和距离rsun如下:
通过旋转转换成赤道坐标系的笛卡尔坐标:
由于黄纬B⊙在1′精度下为0,将式(53)简化为:
其中:rsun表示太阳位置;
日月摄动加速度对位置的偏导数:
其中:r表示卫星位置,s表示太阳或月亮的位置;
日月摄动加速度对速度的偏导数:
步骤3中太阳光压引起的卫星加速度的具体确定方法为:
通过式(57)计算太阳光压加速度:
其中:ν为蚀因子;P⊙为地球附近太阳光压强参数,其值约为4.56×10-6N/m2;光压系数CR=1+ε,卫星制造中典型材料的反射系数ε为0.2~0.9;A/m为卫星的面质比;r⊙和r⊙分别是卫星到太阳的向量及其模;AU为天文单位表示距离,1AU=149597870.691km;
太阳光压加速度对卫星位置的偏导数:
其中:r表示卫星位置,rsun表示太阳位置;
太阳光压加速度对卫星速度的偏导数:
太阳光压对卫星等效面质比的偏导数:
同时在两个坐标系中表述的光压摄动力,一个是基于地心的BYD坐标系,另一个是基于卫星本体的XYZ坐标系;
其中:rsun、r为惯性系下太阳和卫星位置,eD定义为太阳至卫星方向的单位矢量,ex,ey,ez为星固坐标系坐标轴的单位矢量,其中:ez指向地球中心,ey=ez×eD,ex=ey×ez,eB=eD×ey;
u定义为卫星在轨道平面上距升交点的角度,u0代表的是太阳的升交角距;
D(β)=λ(SRP(1)·D0+DC2cos(2β)+DC4cos(4β)) (62)
其中:DC2=-0.813×10-9m/s2,DC4=0.517×10-9m/s2,β定义为太阳相对于卫星轨道平面的高度角D0为ROCK模型计算出来的太阳辐射压产生的加速度的理论值,SRP为待估参数;
Y(β)=SRP(2)·D0+YCcos(2β) (63)
其中:YC=0.067×10-9m/s2;
B(β)=SRP(3)·D0+BCcos(2β) (64)
其中:BC=-0.385×10-9m/s2;
X1(β)=SRP(4)·D0+X10+X1Ccos(2β)+X1Ssin(2β) (65)
其中:X10=-0.015×10-9m/s2,X1C=-0.018×10-9m/s2,X1S=-0.033×10-9m/s2
X3(β)=SRP(5)·D0+X30+X3Ccos(2β)+X3Ssin(2β) (66)
其中:X30=0.004×10-9m/s2,X3C=-0.046×10-9m/s2,X3S=-0.398×10-9m/s2;
z(β)=SRP(6)·D0+z0+zC2cos(2β)+zS2sin(2β)+zc4cos(4β)+zS4cos(4β) (67)
其中:Z0(BlockII)=1.024×10-9m/s2,Z0(BlockIIA)=0.979×10-9m/s2,ZC2=0.519×10-9m/s2,ZS2=0.125×10-9m/s2,
ZC4=0.047×10-9m/s2,ZS4=-0.045×10-9m/s2。
步骤3采用Adams-Cowell积分法求解式(22)进行轨道预报的预报校正方法为:
9-1:采用Adams显示公式计算速度,Cowell显示公式计算位置;
9-2:预报:
9-3:根据预报的tn+1时刻的yn+1,求解相应的tn+1处右函数值fn+1;
9-4:校正:
9-5:回到9-2,重复9-2~9-4,至求得所有时刻的卫星位置和速度。
步骤4的具体方法为:
10-1:轨道根数拟合;具体为:
计算卫星位置需用到16个星历参数,其中星历表参考历元toe是已知的,需要拟合15个参数;待估状态参数向量和观测方程为:
γ=γ(X,t) (69)
其中:X为参考历元时刻的星历参数,Y为一个含m(m≥15)个观测量的观测列向量;
设Xi为参数X在第i次迭代的初值,将观测方程在所给初值处展开,并舍去二阶和二阶以上的小量后得:
其中:Y(Xi,t)为用参考历元toe时刻广播星历参数初值计算的卫星位置,分别为相应星历参数的改正值,为观测量对星历参数的偏导数;令:
L=Y-Y(Xi,t) (72)
得误差方程:
V=HδXi-L (73)
由最小二乘原理有:
δXi=(HTH)-1HTL (74)
则第i次迭代后的星历参数估值为:
X=Xi+δXi (75)
在程序中,所选用的迭代结束条件为:
Di-Di-1<1e-8
其中:k是定轨所用的历元数,n是待估状态参数X的维数,下标i表第i次迭代;
10-2:求解卫星位置对星历参数的偏导数:
根据卫星位置对星历参数的偏导数计算公式:
给出用来推导GPS16参数偏导数的用户位置的最终计算表达式:
其中:卫星位置对轨道长半径的平方根的偏导数,卫星位置对轨道偏心率的偏导数,卫星位置对平近点角的偏导数,卫星位置对平均角速度的偏导数,卫星位置对升交点赤经的偏导数,卫星位置对升交点赤经变化率的偏导数,卫星位置对轨道倾角的偏导数,卫星位置对轨道倾角变化率的偏导数,卫星位置对近地点角距的偏导数,卫星位置对升交点余弦调和项校正的偏导数,卫星位置对升交点正弦调和项校正的偏导数,卫星位置对轨道半径余弦调和项校正的偏导数,卫星位置对轨道半径正弦调和项校正的偏导数,卫星位置对轨道倾角余弦调和项校正的偏导数,卫星位置对轨道倾角正弦调和项校正的偏导数;
10-3:选取轨道根数拟合法初值;
以二体问题的方法求解初始6个密切开普勒轨道根数,其余9个摄动调和参数初值设为基于两组位置矢量的位置确定:
ra,rb取某一固定惯性坐标系下的位置,在程序中取周内秒为0时刻的ECEF坐标系的三轴方向为固定坐标系作为自定义惯性坐标系,计算卫星在坐标系下ta,tb时刻的位置;
其中:θ1等于地球球自转角速度与ta时刻对应的周内秒的成乘积,θ2等于地球自转角速度与tb时刻对应的周内秒的乘积;
高斯矢量:
W=ea×e0
其中:ea为ra的单位矢量,e0为r0的单位矢量,ea与e0相互垂直;
倾角i和升交点赤经Ω如下给出:
纬度角:
将根据两组位置矢量和时间间隔来确定轨道的问题转化为同求解轨道和向径组成的扇形和三角形的比例问题;
令ra和rb表示卫星最在时间ta和tb的地心距;由矢量ra和rb围成的三角形的面积Δ依赖于边长ra和rb以及两矢量之间的夹角vb-va,假定该角度小于180°;三角形面积为:
其中:va和vb为不同时刻的真近点角;
扇形面积S是指矢量ra和rb以及轨道弧线围城的区域,根据开普勒第二定律此面积大小与ta和tb之间的时间差成正比:
其中:a和e为轨道的长半轴和轨道的偏心率;代入半通径p=a(1-e2),得:
其中,η是扇形与三角形面积之比,归一化的时间间隔τ定义如下:
采用二体问题方程式的已知量来替代半通径,将η转换成包含两个方程的系统:
其中:正、辅助变量为:
由式(82)得到l的最小值为
根据g确定η,g的值等于tb时刻偏近点角与ta时刻偏近点角的差值的1/2;消掉g得到超越方程:
其中:W函数如下定义:
变量ω对椭圆轨道始终为正,并且小于1位迭代确定η,使用割线法:
求解式(86)的根:
选取初始值:
η1=η0+0.1 η2=η0 (87)
由Hansen近似法给出:
解得η后,根据式(79),半通径用时间间隔τ表示如下:
其中:由ra和rb确定的三角形面积如下:
根据圆锥曲线方程得出轨道偏心率,求解e·cosν时,有:
由于:
得到两个方程:
求解ta时刻的偏心率和真近点角:
近地点辐角是纬度角和真近点角之差:
ω=ua-va (95)
根据半通径和偏心率,得到长半轴:
最后,由开普勒方程计算平近点角Ma:
Ma=Ea-e·sinEa (rad) (97)
偏近点角Ea满足方程:
10-4:将预测出的卫星地心地固坐标系进行拟合,得到与GPS接口控制文件结构一致的预测广播星历。
与现有技术相比,本发明具有以下有益效果:
通过精密定轨,确定定轨时刻的精确位置速度,并以定轨时刻的精确位置速度为起始,进行轨道外推,外推到开机时刻,将得到的开机时刻附近的位置速度量拟合成广播星历,得到与GPS接口控制文件结构一致的预测广播星历,此预测广播星历可直接应用于接收机位置的解算,无需等待广播星历的收集,大幅缩短非热启动时的首次定位时长;在关机超过2小时再次开机时,首次定位时长从30秒减少至在10秒以内。同时,预测广播星历有效期为4小时,占用空间远小于存储4小时的卫星地心地固坐标,降低对芯片存储空间的要求,有利于芯片的成本控制。
附图说明
图1为本发明的方法流程框图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分的实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都应当属于本发明保护的范围。术语“包括”和“具有”以及他们的任何变形,意图在于覆盖不排他的包含,例如,包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元,而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
下面结合附图对本发明做进一步详细描述:
参见图1,本发明一种用于GPS接收机自主预测星历的方法,包括以下步骤:
步骤1:根据输入的广播星历计算星历有效时间内卫星的位置速度矢量,作为精密定轨的观测值;其中,广播星历的输入可以通过以下方式:
对于每个卫星,从导航电文中获得广播星历至GPS接收机,但是不以此为限。
步骤2:根据观测值进行精密定轨,利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量,以及太阳光压模型参数中的卫星面质比。
定轨计算实际是利用含有误差的测量数据,结合航天器受到的各种外力作用,获取航天器状态的最佳估计值。由于在航天定轨的同时确定一些轨道相关的动力学或测量上的物理参数和几何参数,又称为精密定轨——轨道确定与参数估计。
本定轨程序中,首先,根据输入的一组或多组数据星历,计算得到一段时间(t0,t1…tk)的卫星位置/速度观测量。然后,利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量,以及太阳光压模型参数中的卫星面质比Am2。该面质比值将作为常数用于后面的轨道预报中。
具体方法为:建立航天器的状态量模型:
其中:r,v表示t时刻卫星的位置速度矢量,β表示待估参数卫星等效面质比;
观测量Y与航天器的状态量X的量测方程:
Y=Y(X,t) (2)
其中:Y表示观测量的测量值;X(t)满足轨道动力学方程:
求解式(3)得到状态方程:
X(t)=X(t0,X0;t) (4)
其中:X0为待估状态量;
将式(2)在参考状态量并舍去二阶和二阶以上的小量,得到:
其中:参考状态量为待估状态量X0的近似值;Y(X*,t)是观测量的理论计算值,是测量矩阵,记为HX;是状态转移矩阵,记为Φ(t0,t);
令:
得到精密定轨基本方程:
y=Bx0 (8)
其中:y为残差;
通过迭代方式根据观测采样数据Yj(j=1,…k)求解精密定轨基本方程(8),关机前保存的一组有效期星历,通过按照GPS接口规范文档可以计算出星历有效期期间任意时刻的卫星位置作为观测采样数据;得到待估状态量X0的改正值给出历元t0时刻达到精度要求的状态量则:(9)
其中:下标k表示引用了k次采样数据;重复迭代i次直至满足精度要求,得到:
其中:精密定轨基本方程(8)的求解采用最小二乘估计,通过式(11)计算得到:
Bl=HlΦl,0
其中:表示从t0到tl时刻的状态转移矩阵,Hl表示tl时刻的测量矩阵,yl表示tl时刻的残差向量;
矩阵B通过式(12)计算:
测量矩阵HX通过式(13)计算:
其中:x,y,z表示三个位置分量;
状态转移矩阵Φ(t0,t)通过状态方程(4)得到,状态转移矩阵Φ(t0,t)通过式(14)计算:
得到状态转移矩阵Φ(t0,t)的具体方法为:令:X(t)=Φ(t,t0)X0 (15)
则:
状态转移矩阵Φ(t0,t)所满足的条件通过对状态微分方程(3)的线性化得到,有:
其中:状态微分方程F为:
记:
则方程(17)的系数矩阵A:
将方程(17)转换为一阶方程组:
其中:C1(t),C2(t),C3(t)均用参考值X*进行计算;对式(21)进行积分得到通过得到状态转移矩阵Φ(t0,t)。
步骤3:进行轨道预报,以定轨弧段初始时刻的的精确位置速度为起始,采用卫星面质比作为卫星等效面质比参数进行轨道外推,外推到开机时刻;将轨道预报问题描述为一阶常微分方程初值问题:
其中:卫星的加速度且满足r为卫星位置,为卫星速度;f(t,y)为重力场、日月引力和太阳光压引起的卫星加速度;
采用Adams-Cowell积分法求解式(22)进行轨道预报,并采用8阶龙哥库塔方法作为Adams-Cowell积分法的起步单步法。
其中:重力场引起的卫星加速度采用引力加速度估计函数法计算,重力场的偏导数计算方法为:
6-1:将重力场分为地球质心引力部分和地球带谐摄动加速度、田谐摄动部分两部分;
6-2:地球质心引力加速度偏导数计算:
地球中心引力加速度表示如下:
其中,r=[x y z]T表示卫星在地固系下的位置矢量,r=|r|;
对式(23)求偏导数得到地球质心引力加速度对位置的偏导数:
地球质心引力加速度对速度的偏导数:
6-3:地球带谐摄动加速度、田谐摄动加速度的偏导数计算:
地球重力势:
其中:R表示地球半径,Cnm和Snm表示未归一化的重力势系数,根据式(27)把归一化后的重力势系数和转化成未归一化前的Cnm和Snm:
Vnm和Wnm的计算公式如下:
已知:
利用式(29)置m=0获得带谐项Vn0,相应的值Wn0均等于0,递归使用式(31):
(31)中过程采用式(28),↓过程采用式(29),计算得到所有的Vnm和Wnm,
加速度等于势能u的梯度,直接由Vnm和Wnm计算得出:
其中:
地球带谐摄动加速度、田谐摄动加速度对位置的偏导数根据:
通过式(40)计算得到:
地球带谐摄动加速度、田谐摄动加速度对速度的偏导数:
其中:日月引力引起的卫星加速度的具体确定方法为:
在以地球为中心的参考系下,太阳和月球的摄动表示为:
其中:r表示卫星位置,s表示太阳或月亮的位置;低精度日月坐标具体计算方法如下:
T=(MJD-51544.5)/36525.0 (43)
摄动的计算基于月球平黄经L0、月球平近点角l、太阳平近点角l′、太阳平黄经和月球平黄经之间的差D以及月球平升交距角F,5个基础角;
L0=(0.606433+1336.851344T-floor(0.606433+1336.851344T)) (44)
记相对于2000年黄道和春分点的月球黄经为LM,利用式(45)计算值,有:
LM=2π(L0+dL/1296.0e3-floor(L0+dL/1296.0e3)) (47)
月球黄纬BM表示为:
月球的地心距rM表示为:
月球坐标rM为:
其中:ε为黄赤交角,ε=23°.43929111;
假设地球围绕太阳作无摄运动,通过式(51)描述太阳相对于地球和黄道在2000年附近几十年的椭圆轨道:
位置坐标利用式(51)由开普勒轨道公式得到;
太阳的黄经L⊙和距离rsun如下:
通过旋转转换成赤道坐标系的笛卡尔坐标:
由于黄纬B⊙在1′精度下为0,将式(53)简化为:
其中:rsun表示太阳位置;
日月摄动加速度对位置的偏导数:
其中:r表示卫星位置,s表示太阳或月亮的位置;
日月摄动加速度对速度的偏导数:
其中:太阳光压引起的卫星加速度的具体确定方法为:
通过式(57)计算太阳光压加速度:
其中:ν为蚀因子;P⊙为地球附近太阳光压强参数,其值约为4.56×10-6N/m2;光压系数CR=1+ε,卫星制造中典型材料的反射系数ε为0.2~0.9;A/m为卫星的面质比;r⊙和r⊙分别是卫星到太阳的向量及其模;AU为天文单位表示距离,1AU=149597870.691km;
太阳光压加速度对卫星位置的偏导数:
其中:r表示卫星位置,rsun表示太阳位置;
太阳光压加速度对卫星速度的偏导数:
太阳光压对卫星等效面质比的偏导数:
同时在两个坐标系中表述的光压摄动力,一个是基于地心的BYD坐标系,另一个是基于卫星本体的XYZ坐标系;
其中:rsun、r为惯性系下太阳和卫星位置,eD定义为太阳至卫星方向的单位矢量,ex,ey,ez为星固坐标系坐标轴的单位矢量,其中:ez指向地球中心,ey=ez×eD,ex=ey×ez,eB=eD×ey;
u定义为卫星在轨道平面上距升交点的角度,u0代表的是太阳的升交角距;
D(β)=λ(SRP(1)·D0+DC2cos(2β)+DC4cos(4β)) (62)
其中:DC2=-0.813×10-9m/s2,DC4=0.517×10-9m/s2,β定义为太阳相对于卫星轨道平面的高度角,D0为ROCK模型计算出来的太阳辐射压产生的加速度的理论值,SRP为待估参数;
Y(β)=SRP(2)·D0+γCcos(2β) (63)
其中:γC=0.067×10-9m/s2;
B(β)=SRP(3)·D0+BCcos(2β) (64)
其中:BC=-0.385×10-9m/s2;
X1(β)=SRP(4)·D0+X10+X1Ccos(2β)+X1Ssin(2β) (65)
其中:X10=-0.015×10-9m/s2,X1C=-0.018×10-9m/s2,X1S=-0.033×10-9m/s2
X3(β)=SRP(5)·D0+X30+X3Ccos(2β)+X3Ssin(2β) (66)
其中:X30=0.004×10-9m/s2,X3C=-0.046×10-9m/s2,X3S=-0.398×10-9m/s2;
z(β)=SRP(6)·D0+z0+zC2cos(2β)+zS2sin(2β)+zc4cos(4β)+zS4cos(4β) (67)
其中:Z0(BlockII)=1.024×10-9m/s2,Z0(BlockIIA)=0.979×10-9m/s2,ZC2=0.519×10-9m/s2,ZS2=0.125×10-9m/s2,ZC4=0.047×10-9m/s2,ZS4=-0.045×10-9m/s2。
其中采用Adams-Cowell积分法求解式(22)进行轨道预报的具体内容如下:
Adams-Cowell是多步法,计算速度较快,但多步法需要用到前两个以上时刻的计算结果,因此它不是自起步的。在实际应用中,多步法都是需要用单步法来进行起步的,本发明采用单步的RK8方法来起步。A-C方法是Adams方法和Cowell方法的结合,因此在介绍AC法前,首先给出Adams和Cowell这两种算法的计算过程。
1、Adams算法:适用于一阶常微分方程形如:
显示公式:
其中:表示二项式系数,即
隐式公式:
2、Cowell算法
Cowell方法可对二阶微分方程直接积分。二阶常微分方程如:
可以看出,该方程的右函数中不含(对卫星轨道运动方程来说,就是右函数不含速度项,仅考虑引力作用时的运动方程就是如此),因此在每一步计算中,只需直接给出yn即可,而不用去计算这样采用Cowell方法较为方便。
Cowell显示公式:
其中:
Cowell隐式公式:
且
3、Adams-Cowell算法
一般来说,Cowell公式比同阶的Adams公式精度要高。对于k阶Adams公式,截断误差为O(hk+1),而对于k阶Cowell公式截断误差为O(hk+2)。但是在卫星运动方程中,常常是要出现也就是速度项的。对于初值问题:可以将Adams方法和Cowell方法联合起来使用,由Adams公式提供(速度矢量),而用Cowell提供y(位置矢量),这将比单用Adams公式更加有效。另外,隐式方法比显示方法精度高,且数值稳定性好,但是需要一些迭代过程来求解,计算比较复杂。因此在求解卫星运动方程中,常常使用预估—校正(PECE)算法,将显示公式和隐式公式结合起来使用。P即预报,用显示公式将解从tn预报到tn+1点处;E即估计,用tn+1点处预报值计算相应的右函数值;C即改正,将计算出的右函数值通过隐式公式得到改正的预报值;最后一个E仍是估计,但是是用修正后的tn+1点处预报值再去计算右函数,并以此作为下一步积分的起始值。
预报校正方法为:
9-1:采用Adams显示公式计算速度,Cowell显示公式计算位置;
9-2:预报:
9-3:根据预报的tn+1时刻的yn+1,求解相应的tn+1处右函数值fn+1;
9-4:校正:
9-5:回到9-2,重复9-2~9-4,至求得所有时刻的卫星位置和速度。
步骤4:将开机时刻的位置速度量拟合成广播星历,给出预报星历,将预报星历直接应用于接收机位置的定位解算,具体方法为:
10-1:轨道根数拟合;具体为:
计算卫星位置需用到16个星历参数,其中星历表参考历元toe是已知的,需要拟合15个参数;待估状态参数向量和观测方程为:
Y=Y(X,t) (69)
其中:X为参考历元时刻的星历参数,Y为一个含m(m≥15)个观测量的观测列向量;
设Xi为参数X在第i次迭代的初值,将观测方程在所给初值处展开,并舍去二阶和二阶以上的小量后得:
其中:Y(Xi,t)为用参考历元toe时刻广播星历参数初值计算的卫星位置,分别为相应星历参数的改正值,为观测量对星历参数的偏导数;令:
L=Y-γ(Xi,t) (72)
可得误差方程:V=HδXi-L;
由最小二乘原理有:
δXi=(HTH)-1HTL (74)
则第i次迭代后的星历参数估值为:
X=Xi+δXi (75)
在程序中,所选用的迭代结束条件为:
其中:k是定轨所用的历元数,n是待估状态参数X的维数,下标i表第i次迭代;
10-2:求解卫星位置对星历参数的偏导数:
根据卫星位置对星历参数的偏导数计算公式:
给出用来推导GPS16参数偏导数的用户位置的最终计算表达式:
其中:卫星位置对轨道长半径的平方根的偏导数,卫星位置对轨道偏心率的偏导数,卫星位置对平近点角的偏导数,卫星位置对平均角速度的偏导数,卫星位置对升交点赤经的偏导数,卫星位置对升交点赤经变化率的偏导数,卫星位置对轨道倾角的偏导数,卫星位置对轨道倾角变化率的偏导数,卫星位置对近地点角距的偏导数,卫星位置对升交点余弦调和项校正的偏导数,卫星位置对升交点正弦调和项校正的偏导数,卫星位置对轨道半径余弦调和项校正的偏导数,卫星位置对轨道半径正弦调和项校正的偏导数,卫星位置对轨道倾角余弦调和项校正的偏导数,卫星位置对轨道倾角正弦调和项校正的偏导数;
10-3:选取轨道根数拟合法初值;
以二体问题的方法求解初始6个密切开普勒轨道根数,其余9个摄动调和参数初值设为基于两组位置矢量的位置确定:
ra,rb取某一固定惯性坐标系下的位置,在程序中取周内秒为0时刻的ECEF坐标系的三轴方向为固定坐标系作为自定义惯性坐标系,计算卫星在坐标系下ta,tb时刻的位置;
其中:θ1等于地球球自转角速度与ta时刻对应的周内秒的成乘积,θ2等于地球自转角速度与tb时刻对应的周内秒的乘积;
高斯矢量:
W=ea×e0
其中:ea为ra的单位矢量,e0为r0的单位矢量,ea与e0相互垂直;
倾角i和升交点赤经Ω如下给出:
纬度角:
将根据两组位置矢量和时间间隔来确定轨道的问题转化为同求解轨道和向径组成的扇形和三角形的比例问题;
令ra和rb表示卫星最在时间ta和tb的地心距;由矢量ra和rb围成的三角形的面积Δ依赖于边长ra和rb以及两矢量之间的夹角vb-va,假定该角度小于180°;三角形面积为:
其中:va和vb为不同时刻的真近点角;
扇形面积S是指矢量ra和rb以及轨道弧线围城的区域,根据开普勒第二定律此面积大小与ta和tb之间的时间差成正比:
其中:a和e为轨道的长半轴和轨道的偏心率;代入半通径p=a(1-e2),得:
其中,η是扇形与三角形面积之比,归一化的时间间隔τ定义为:
采用二体问题方程式的已知量来替代半通径,将η转换成包含两个方程的系统:
其中:正、辅助变量为:
由式(82)得到l的最小值为
根据g确定η,g的值等于tb时刻偏近点角与ta时刻偏近点角的差值的1/2;消掉g得到超越方程:
其中:W函数如下定义:
变量ω对椭圆轨道始终为正,并且小于1位迭代确定η,使用割线法:
求解式(86)的根:
选取初始值:
η1=η0+0.1 η2=η0 (87)
由Hansen近似法给出:
解得η后,根据式(79),半通径用时间间隔τ表示如下:
其中:由ra和rb确定的三角形面积为:
根据圆锥曲线方程得出轨道偏心率,求解e·cosν时,有:
由于:
得到两个方程:
求解ta时刻的偏心率和真近点角:
近地点辐角是纬度角和真近点角之差:ω=ua-va(95)
根据半通径和偏心率,得到长半轴:
最后,由开普勒方程计算平近点角Ma:Ma=Ea-e·sinEa (rad) (97)
偏近点角Ea满足方程:
10-4:将预测出的卫星地心地固坐标系进行拟合,得到与GPS接口控制文件结构一致的预测广播星历。
仿真实验表明,在关机超过2小时再次开机时,没有星历预测首次定位时长为30秒以上,采用本方法首次定位时长在10秒以内。
以上内容仅为说明本发明的技术思想,不能以此限定本发明的保护范围,凡是按照本发明提出的技术思想,在技术方案基础上所做的任何改动,均落入本发明权利要求书的保护范围之内。
Claims (10)
1.一种用于GPS接收机自主预测星历的方法,其特征在于,包括以下步骤:
步骤1:根据输入的广播星历计算星历有效时间内卫星的位置速度矢量,作为精密定轨的观测值;
步骤2:根据所述观测值进行精密定轨,利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量,以及太阳光压模型参数中的卫星面质比;
步骤3:进行轨道预报,以所述定轨弧段初始时刻的精确位置速度为起始,采用所述卫星面质比作为卫星等效面质比参数进行轨道外推,外推到开机时刻,得到开机时刻的位置速度矢量;
步骤4:将所述开机时刻的位置速度矢量拟合成广播星历,直接应用于接收机的定位解算。
2.根据权利要求1所述的用于GPS接收机自主预测GPS星历的方法,其特征在于,所述步骤2中进行精密定轨的具体方法为:
建立航天器的状态量模型:
其中:r,v分别表示t时刻卫星的位置矢量和速度矢量,β表示待估参数卫星等效面质比;
观测量与航天器的状态量X的量测方程:
Y=Y(X,t) (2)
其中:Y表示观测量的测量值且只包含位置矢量;X(t)满足轨道动力学方程:
求解式(3)得到状态方程:
X(t)=X(t0,X0;t) (4)
其中:X0为待估状态量;
将式(2)在参考状态量并舍去二阶和二阶以上的小量,得到:
其中:参考状态量为待估状态量X0的近似值;Y(X*,t)是观测量的理论计算值,是测量矩阵,记为HX;是状态转移矩阵,记为Φ(t0,t);
令:
得到精密定轨基本方程:
y=Bx0 (8)
其中:y为残差;
通过迭代方式由观测采样数据求解精密定轨基本方程(8),得到待估状态量X0的改正值给出历元t0时刻达到精度要求的状态量则:
其中:下标k表示引用了k次采样数据;重复迭代i次直至满足精度要求,得到:
3.根据权利要求2所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤2中:
精密定轨基本方程(8)的求解采用最小二乘估计,通过式(11)计算得到:
其中:表示从t0到tl时刻的状态转移矩阵,Hl表示tl时刻的测量矩阵,yl表示tl时刻的残差向量;
矩阵B通过式(12)计算:
测量矩阵HX通过式(13)计算:
其中:x,y,z表示三个位置分量;
状态转移矩阵Φ(t0,t)通过状态方程(4)得到,状态转移矩阵Φ(t0,t)通过式(14)计算:
4.根据权利要求3所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤2中得到状态转移矩阵Φ(t0,t)的具体方法为:
令:
X(t)=Φ(t,t0)X0 (15)
则:
状态转移矩阵Φ(t0,t)所满足的条件通过对状态微分方程(3)的线性化得到,有:
其中:状态微分方程F为:
记:
则方程(17)的系数矩阵A:
将方程(17)转换为一阶方程组:
其中:C1(t),C2(t),C3(t)均用参考值X*进行计算;
对式(21)进行积分得到通过得到状态转移矩阵Φ(t0,t)。
5.根据权利要求1所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤3的具体方法为:
将轨道预报问题描述为一阶常微分方程初值问题:
其中:卫星的加速度且满足r为卫星位置,为卫星速度;f(t,y)为重力场、日月引力和太阳光压引起的卫星加速度;
采用Adams-Cowell积分法求解式(22)进行轨道预报,并采用8阶龙哥库塔方法作为Adams-Cowell积分法的起步单步法。
6.根据权利要求5所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤3中重力场引起的卫星加速度采用引力加速度估计函数法计算,重力场的偏导数计算方法:
6-1:将重力场分为地球质心引力部分和地球带谐摄动加速度、田谐摄动部分两部分;
6-2:地球质心引力加速度偏导数计算:
地球中心引力加速度表示如下:
其中,r=[x y z]T表示卫星在地固系下的位置矢量,r=|r|;
对式(23)求偏导数得到地球质心引力加速度对位置的偏导数:
地球质心引力加速度对速度的偏导数:
6-3:地球带谐摄动加速度、田谐摄动加速度的偏导数计算:
地球重力势:
其中:R表示地球半径,Cnm和Snm表示未归一化的重力势系数,根据式(27)把归一化后的重力势系数和转化成未归一化前的Cnm和Snm:
Vnm和Wnm的计算公式如下:
已知:
利用式(29)置m=0获得带谐项Vn0,相应的值均等于0,递归使用式(31):
式(31)中过程采用式(28),↓过程采用式(29),计算得到所有的Vnm和Wnm,加速度等于势能u的梯度,直接由Vnm和Wnm计算得出:
其中:
地球带谐摄动加速度、田谐摄动加速度对位置的偏导数根据:
通过式(40)计算得到:
地球带谐摄动加速度、田谐摄动加速度对速度的偏导数:
7.根据权利要求5所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤3中日月引力引起的卫星加速度的具体确定方法为:
在以地球为中心的参考系下,太阳和月球的摄动表示为:
其中:r表示卫星位置,s表示太阳或月亮的位置;
低精度日月坐标具体计算方法如下:
T=(MJD-51544.5)/36525.0 (43)
摄动的计算基于月球平黄经L0、月球平近点角l、太阳平近点角l′、太阳平黄经和月球平黄经之间的差D以及月球平升交距角F,5个基础角;
L0=(0.606433+1336.851344T-floor(0.606433+1336.851344T)) (44)
记相对于2000年黄道和春分点的月球黄经为LM,利用式(45)计算值,有:
LM=2π(L0+dL/1296.0e3-floor(L0+dL/1296.0e3)) (47)
月球黄纬BM表示为:
月球的地心距rM表示为:
月球坐标rM为:
其中:ε为黄赤交角,ε=23°.43929111;
假设地球围绕太阳作无摄运动,通过式(51)描述太阳相对于地球和黄道在2000年附近几十年的椭圆轨道:
位置坐标利用式(51)由开普勒轨道公式得到;
太阳的黄经L⊙和距离rsun如下:
通过旋转转换成赤道坐标系的笛卡尔坐标:
由于黄纬B⊙在1′精度下为0,将式(53)简化为:
其中:rsun表示太阳位置;
日月摄动加速度对位置的偏导数:
其中:r表示卫星位置,s表示太阳或月亮的位置;
日月摄动加速度对速度的偏导数:
8.根据权利要求5所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤3中太阳光压引起的卫星加速度的具体确定方法为:
通过式(57)计算太阳光压加速度:
其中:ν为蚀因子;P⊙为地球附近太阳光压强参数,其值约为4.56×10-6N/m2;光压系数CR=1+ε,卫星制造中典型材料的反射系数ε为0.2~0.9;A/m为卫星的面质比;r⊙和r⊙分别是卫星到太阳的向量及其模;AU为天文单位表示距离,1AU=149597870.691km;
太阳光压加速度对卫星位置的偏导数:
其中:r表示卫星位置,rsun表示太阳位置;
太阳光压加速度对卫星速度的偏导数:
太阳光压对卫星等效面质比的偏导数:
同时在两个坐标系中表述的光压摄动力,一个是基于地心的BYD坐标系,另一个是基于卫星本体的XYZ坐标系;
其中:rsun、r为惯性系下太阳和卫星位置,eD定义为太阳至卫星方向的单位矢量,ex,ey,ez为星固坐标系坐标轴的单位矢量,其中:ez指向地球中心,ey=ez×eD,ex=ey×ez,eB=eD×ey;
u定义为卫星在轨道平面上距升交点的角度,u0代表的是太阳的升交角距;
D(β)=λ(SRP(1)·D0+DC2cos(2β)+DC4cos(4β)) (62)
其中:DC2=-0.813×10-9m/s2,DC4=0.517×10-9m/s2,β定义为太阳相对于卫星轨道平面的高度角D0为ROCK模型计算出来的太阳辐射压产生的加速度的理论值,SRP为待估参数;
Y(β)=SRP(2)·D0+YCcos(2β) (63)
其中:γC=0.067×10-9m/s2;
B(β)=SRP(3)·D0+BCcos(2β) (64)
其中:BC=-0.385×10-9m/s2;
X1(β)=SRP(4)·D0+X10+X1Ccos(2β)+X1Ssin(2β) (65)
其中:X10=-0.015×10-9m/s2,X1C=-0.018×10-9m/s2,X1S=-0.033×10-9m/s2
X3(β)=SRP(5)·D0+X30+X3Ccos(2β)+X3Ssin(2β) (66)
其中:X30=0.004×10-9m/s2,X3C=-0.046×10-9m/s2,X3S=-0.398×10-9m/s2;
z(β)=SRP(6)·D0+z0+zC2cos(2β)+zS2sin(2β)+zc4cos(4β)+zS4cos(4β) (67)
其中:Z0(BlockII)=1.024×10-9m/s2,Z0(BlockIIA)=0.979×10-9m/s2,ZC2=0.519×10-9m/s2,ZS2=0.125×10-9m/s2,ZC4=0.047×10-9m/s2,ZS4=-0.045×10-9m/s2。
9.根据权利要求5所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤3采用Adams-Cowell积分法求解式(22)进行轨道预报的预报校正方法为:
9-1:采用Adams显示公式计算速度,Cowell显示公式计算位置;
9-2:预报:
9-3:根据预报的tn+1时刻的yn+1,求解相应的tn+1处右函数值fn+1;
9-4:校正:
9-5:回到9-2,重复9-2~9-4,至求得所有时刻的卫星位置和速度。
10.根据权利要求1所述的用于GPS接收机自主预测星历的方法,其特征在于,所述步骤4的具体方法为:
10-1:轨道根数拟合;具体为:
计算卫星位置需用到16个星历参数,其中星历表参考历元toe是已知的,需要拟合15个参数;待估状态参数向量和观测方程为:
Y=Y(X,t) (69)
其中:X为参考历元时刻的星历参数,Y为一个含m(m≥15)个观测量的观测列向量;
设Xi为参数X在第i次迭代的初值,将观测方程在所给初值处展开,并舍去二阶和二阶以上的小量后得:
其中:Y(Xi,t)为用参考历元toe时刻广播星历参数初值计算的卫星位置,分别为相应星历参数的改正值,为观测量对星历参数的偏导数;令:
L=Y-Y(Xi,t) (72)
得误差方程:
V=HδXi-L (73)
由最小二乘原理有:
δXi=(HTH)-1HTL (74)
则第i次迭代后的星历参数估值为:
X=Xi+δXi (75)
在程序中,所选用的迭代结束条件为:
其中:k是定轨所用的历元数,n是待估状态参数X的维数,下标i表第i次迭代;
10-2:求解卫星位置对星历参数的偏导数:
根据卫星位置对星历参数的偏导数计算公式:
给出用来推导GPS16参数偏导数的用户位置的最终计算表达式:
其中:卫星位置对轨道长半径的平方根的偏导数,卫星位置对轨道偏心率的偏导数,卫星位置对平近点角的偏导数,卫星位置对平均角速度的偏导数,卫星位置对升交点赤经的偏导数,卫星位置对升交点赤经变化率的偏导数,卫星位置对轨道倾角的偏导数,卫星位置对轨道倾角变化率的偏导数,卫星位置对近地点角距的偏导数,卫星位置对升交点余弦调和项校正的偏导数,卫星位置对升交点正弦调和项校正的偏导数,卫星位置对轨道半径余弦调和项校正的偏导数,卫星位置对轨道半径正弦调和项校正的偏导数,卫星位置对轨道倾角余弦调和项校正的偏导数,卫星位置对轨道倾角正弦调和项校正的偏导数;
10-3:选取轨道根数拟合法初值;
以二体问题的方法求解初始6个密切开普勒轨道根数,其余9个摄动调和参数初值设为基于两组位置矢量的位置确定:
ra,rb取某一固定惯性坐标系下的位置,在程序中取周内秒为0时刻的ECEF坐标系的三轴方向为固定坐标系作为自定义惯性坐标系,计算卫星在坐标系下ta,tb时刻的位置;
其中:θ1等于地球球自转角速度与ta时刻对应的周内秒的成乘积,θ2等于地球自转角速度与tb时刻对应的周内秒的乘积;
高斯矢量:
W=ea×e0
其中:r0=rb-(rb·ea)ea,ea为ra的单位矢量,e0为r0的单位矢量,ea与e0相互垂直;
倾角i和升交点赤经Ω如下给出:
纬度角:
将根据两组位置矢量和时间间隔来确定轨道的问题转化为同求解轨道和向径组成的扇形和三角形的比例问题;
令ra和rb表示卫星最在时间ta和tb的地心距;由矢量ra和rb围成的三角形的面积Δ依赖于边长ra和rb以及两矢量之间的夹角vb-va,假定该角度小于180°;三角形面积为:
其中:va和vb为不同时刻的真近点角;
扇形面积S是指矢量ra和rb以及轨道弧线围城的区域,根据开普勒第二定律此面积大小与ta和tb之间的时间差成正比:
其中:a和e为轨道的长半轴和轨道的偏心率;代入半通径p=a(1-e2),得:
其中,η是扇形与三角形面积之比,归一化的时间间隔τ定义如下:
采用二体问题方程式的已知量来替代半通径,将η转换成包含两个方程的系统:
其中:正、辅助变量为:
由式(82)得到l的最小值为
根据g确定η,g的值等于tb时刻偏近点角与ta时刻偏近点角的差值的1/2;消掉g得到超越方程:
其中:W函数如下定义:
变量ω对椭圆轨道始终为正,并且小于1位迭代确定η,使用割线法:
求解式(86)的根:
选取初始值:
η1=η0+0.1 η2=η0 (87)
由Hansen近似法给出:
解得η后,根据式(79),半通径用时间间隔τ表示如下:
其中:由ra和rb确定的三角形面积如下:
根据圆锥曲线方程得出轨道偏心率,求解e·cosν时,有:
由于:
得到两个方程:
求解ta时刻的偏心率和真近点角:
近地点辐角是纬度角和真近点角之差:
ω=ua-va (95)
根据半通径和偏心率,得到长半轴:
最后,由开普勒方程计算平近点角Ma:
Ma=Ea-e·sinEa (rad) (97)
偏近点角Ea满足方程:
10-4:将预测出的卫星地心地固坐标系进行拟合,得到与GPS接口控制文件结构一致的预测广播星历。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910150844.7A CN109738919B (zh) | 2019-02-28 | 2019-02-28 | 一种用于gps接收机自主预测星历的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910150844.7A CN109738919B (zh) | 2019-02-28 | 2019-02-28 | 一种用于gps接收机自主预测星历的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109738919A true CN109738919A (zh) | 2019-05-10 |
CN109738919B CN109738919B (zh) | 2020-12-15 |
Family
ID=66368772
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910150844.7A Active CN109738919B (zh) | 2019-02-28 | 2019-02-28 | 一种用于gps接收机自主预测星历的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109738919B (zh) |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110646819A (zh) * | 2019-10-09 | 2020-01-03 | 四川灵通电讯有限公司 | 低轨卫星星历预报装置及应用方法 |
CN111123980A (zh) * | 2019-12-31 | 2020-05-08 | 贵阳欧比特宇航科技有限公司 | 一种卫星飞临时刻与拍摄范围的计算方法 |
CN111209523A (zh) * | 2020-01-06 | 2020-05-29 | 中国科学院紫金山天文台 | 适用于大偏心率轨道密集星历精密计算的快速处理方法 |
CN111650613A (zh) * | 2020-05-30 | 2020-09-11 | 重庆金美通信有限责任公司 | 一种分布式星历计算方法 |
CN112161632A (zh) * | 2020-09-23 | 2021-01-01 | 北京航空航天大学 | 一种基于相对位置矢量测量的卫星编队初始定位算法 |
CN112666583A (zh) * | 2020-12-15 | 2021-04-16 | 上海卫星工程研究所 | 适应gnss接收机输出状态的单拍轨道递推方法及系统 |
CN113359510A (zh) * | 2021-06-04 | 2021-09-07 | 北京理工大学 | 北斗卫星导航系统信号模拟器数据实时仿真系统及方法 |
CN114440886A (zh) * | 2021-12-30 | 2022-05-06 | 上海航天控制技术研究所 | 一种大偏心率轨道高精度轨道计算方法 |
CN115308779A (zh) * | 2021-05-07 | 2022-11-08 | 华为技术有限公司 | 星历预报方法和星历预报装置 |
CN115792982A (zh) * | 2022-11-07 | 2023-03-14 | 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) | 北斗导航卫星广播星历参数拟合方法及存储介质 |
CN117194869A (zh) * | 2023-11-07 | 2023-12-08 | 中国科学院国家授时中心 | 顾及姿态的低轨卫星天线相位中心预报及拟合方法 |
CN117544214A (zh) * | 2023-05-23 | 2024-02-09 | 国家卫星海洋应用中心 | 卫星对地数据传输方法、装置、设备及可读存储介质 |
Citations (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20080111738A1 (en) * | 2006-11-10 | 2008-05-15 | Shaowei Han | Method and apparatus in standalone positioning without broadcast ephemeris |
CN101435863A (zh) * | 2008-12-25 | 2009-05-20 | 武汉大学 | 一种导航卫星实时精密定轨的方法 |
WO2012125293A2 (en) * | 2011-03-11 | 2012-09-20 | Sorce4 Llc. | Offline ephemeris prediction |
US8538682B1 (en) * | 2009-04-24 | 2013-09-17 | Qualcomm Incorporated | Systems and methods for satellite navigation using locally generated ephemeris data |
US20130257648A1 (en) * | 2007-02-05 | 2013-10-03 | Qualcomm Incorporated | Prediction refresh method for ephemeris extensions |
US20150070213A1 (en) * | 2013-09-12 | 2015-03-12 | Furuno Electric Company, Ltd. | Location determination in multi-system gnns environment using conversion of data into a unified format |
CN104459746A (zh) * | 2013-09-12 | 2015-03-25 | 马维尔国际贸易有限公司 | 用于利用预测星历的位置确定的方法、系统和设备 |
US20150153455A1 (en) * | 2000-11-17 | 2015-06-04 | Broadcom Corporation | Method and apparatus for generating and distributing satellite tracking information |
CN106569241A (zh) * | 2016-09-27 | 2017-04-19 | 北京航空航天大学 | 一种基于gnss的单频高精度定位方法 |
CN107065025A (zh) * | 2017-01-13 | 2017-08-18 | 北京航空航天大学 | 一种基于重力梯度不变量的轨道要素估计方法 |
CN107153209A (zh) * | 2017-07-06 | 2017-09-12 | 武汉大学 | 一种短弧段低轨导航卫星实时精密定轨方法 |
CN108761507A (zh) * | 2018-05-21 | 2018-11-06 | 中国人民解放军战略支援部队信息工程大学 | 基于短弧定轨和预报的导航卫星轨道快速恢复方法 |
-
2019
- 2019-02-28 CN CN201910150844.7A patent/CN109738919B/zh active Active
Patent Citations (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20150153455A1 (en) * | 2000-11-17 | 2015-06-04 | Broadcom Corporation | Method and apparatus for generating and distributing satellite tracking information |
US20080111738A1 (en) * | 2006-11-10 | 2008-05-15 | Shaowei Han | Method and apparatus in standalone positioning without broadcast ephemeris |
US20130257648A1 (en) * | 2007-02-05 | 2013-10-03 | Qualcomm Incorporated | Prediction refresh method for ephemeris extensions |
CN101435863A (zh) * | 2008-12-25 | 2009-05-20 | 武汉大学 | 一种导航卫星实时精密定轨的方法 |
US8538682B1 (en) * | 2009-04-24 | 2013-09-17 | Qualcomm Incorporated | Systems and methods for satellite navigation using locally generated ephemeris data |
WO2012125293A2 (en) * | 2011-03-11 | 2012-09-20 | Sorce4 Llc. | Offline ephemeris prediction |
US20150070213A1 (en) * | 2013-09-12 | 2015-03-12 | Furuno Electric Company, Ltd. | Location determination in multi-system gnns environment using conversion of data into a unified format |
CN104459746A (zh) * | 2013-09-12 | 2015-03-25 | 马维尔国际贸易有限公司 | 用于利用预测星历的位置确定的方法、系统和设备 |
CN106569241A (zh) * | 2016-09-27 | 2017-04-19 | 北京航空航天大学 | 一种基于gnss的单频高精度定位方法 |
CN107065025A (zh) * | 2017-01-13 | 2017-08-18 | 北京航空航天大学 | 一种基于重力梯度不变量的轨道要素估计方法 |
CN107153209A (zh) * | 2017-07-06 | 2017-09-12 | 武汉大学 | 一种短弧段低轨导航卫星实时精密定轨方法 |
CN108761507A (zh) * | 2018-05-21 | 2018-11-06 | 中国人民解放军战略支援部队信息工程大学 | 基于短弧定轨和预报的导航卫星轨道快速恢复方法 |
Non-Patent Citations (4)
Title |
---|
HUI HU ET AL.: "Extrapolation and Fitting Algorithms for GLONASS Satellite Orbit", 《2009 THIRD INTERNATIONAL SYMPOSIUM ON INTELLIGENT INFORMATION TECHNOLOGY APPLICATION》 * |
MARCO ANGHILERI ET AL.: "Reduced Navigation Data for a Fast First Fix", 《 2012 6TH ESA WORKSHOP ON SATELLITE NAVIGATION TECHNOLOGIES (NAVITEC 2012) & EUROPEAN WORKSHOP ON GNSS SIGNALS AND SIGNAL PROCESSING》 * |
楼益栋 等: "基于移动终端广播星历的自主轨道预报及精度分析", 《测绘通报》 * |
闫付龙 等: "GPS信号模拟器轨道外推及星历拟合算法设计与验证", 《第四届中国卫星导航学术年会电子文集》 * |
Cited By (22)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110646819A (zh) * | 2019-10-09 | 2020-01-03 | 四川灵通电讯有限公司 | 低轨卫星星历预报装置及应用方法 |
CN111123980A (zh) * | 2019-12-31 | 2020-05-08 | 贵阳欧比特宇航科技有限公司 | 一种卫星飞临时刻与拍摄范围的计算方法 |
WO2021139129A1 (zh) * | 2020-01-06 | 2021-07-15 | 中国科学院紫金山天文台 | 适用于大偏心率轨道密集星历精密计算的快速处理方法 |
CN111209523A (zh) * | 2020-01-06 | 2020-05-29 | 中国科学院紫金山天文台 | 适用于大偏心率轨道密集星历精密计算的快速处理方法 |
CN111209523B (zh) * | 2020-01-06 | 2020-12-29 | 中国科学院紫金山天文台 | 适用于大偏心率轨道密集星历精密计算的快速处理方法 |
US11319094B2 (en) * | 2020-01-06 | 2022-05-03 | Purple Mountain Observatory, Chinese Academy Of Sciences | Method for accurately and efficiently calculating dense ephemeris of high-eccentricity orbit |
CN111650613A (zh) * | 2020-05-30 | 2020-09-11 | 重庆金美通信有限责任公司 | 一种分布式星历计算方法 |
CN112161632A (zh) * | 2020-09-23 | 2021-01-01 | 北京航空航天大学 | 一种基于相对位置矢量测量的卫星编队初始定位算法 |
CN112666583A (zh) * | 2020-12-15 | 2021-04-16 | 上海卫星工程研究所 | 适应gnss接收机输出状态的单拍轨道递推方法及系统 |
CN115308779B (zh) * | 2021-05-07 | 2023-11-03 | 华为技术有限公司 | 星历预报方法和星历预报装置 |
CN115308779A (zh) * | 2021-05-07 | 2022-11-08 | 华为技术有限公司 | 星历预报方法和星历预报装置 |
WO2022233211A1 (zh) * | 2021-05-07 | 2022-11-10 | 华为技术有限公司 | 星历预报方法和星历预报装置 |
CN113359510A (zh) * | 2021-06-04 | 2021-09-07 | 北京理工大学 | 北斗卫星导航系统信号模拟器数据实时仿真系统及方法 |
CN113359510B (zh) * | 2021-06-04 | 2023-01-31 | 北京理工大学 | 北斗卫星导航系统信号模拟器数据实时仿真系统 |
CN114440886A (zh) * | 2021-12-30 | 2022-05-06 | 上海航天控制技术研究所 | 一种大偏心率轨道高精度轨道计算方法 |
CN114440886B (zh) * | 2021-12-30 | 2023-09-05 | 上海航天控制技术研究所 | 一种大偏心率轨道高精度轨道计算方法 |
CN115792982B (zh) * | 2022-11-07 | 2023-10-20 | 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) | 北斗导航卫星广播星历参数拟合方法及存储介质 |
CN115792982A (zh) * | 2022-11-07 | 2023-03-14 | 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) | 北斗导航卫星广播星历参数拟合方法及存储介质 |
CN117544214A (zh) * | 2023-05-23 | 2024-02-09 | 国家卫星海洋应用中心 | 卫星对地数据传输方法、装置、设备及可读存储介质 |
CN117544214B (zh) * | 2023-05-23 | 2024-05-24 | 国家卫星海洋应用中心 | 卫星对地数据传输方法、装置、设备及可读存储介质 |
CN117194869A (zh) * | 2023-11-07 | 2023-12-08 | 中国科学院国家授时中心 | 顾及姿态的低轨卫星天线相位中心预报及拟合方法 |
CN117194869B (zh) * | 2023-11-07 | 2024-03-19 | 中国科学院国家授时中心 | 顾及姿态的低轨卫星天线相位中心预报及拟合方法 |
Also Published As
Publication number | Publication date |
---|---|
CN109738919B (zh) | 2020-12-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109738919A (zh) | 一种用于gps接收机自主预测星历的方法 | |
Soffel et al. | Relativity in astrometry, celestial mechanics and geodesy | |
Green | Spherical astronomy | |
Ning et al. | Autonomous satellite navigation using starlight refraction angle measurements | |
CN102288201B (zh) | 用于星敏感器的精度测量方法 | |
WO2013070518A1 (en) | Solar timer using gps technology | |
CN102706363B (zh) | 一种高精度星敏感器的精度测量方法 | |
CN102243311B (zh) | 一种x射线脉冲星导航使用的选星方法 | |
Macdonald et al. | Extension of the sun-synchronous orbit | |
CN102288200A (zh) | 用于星敏感器的精度测量系统 | |
Hugentobler et al. | Satellite orbits and attitude | |
Peale | Determination of parameters related to the interior of Mercury | |
Zuber et al. | The Psyche gravity investigation | |
Standish | The JPL planetary ephemerides | |
Kucharski et al. | Confirmation of gravitationally induced attitude drift of spinning satellite Ajisai with Graz high repetition rate SLR data | |
CN102607597B (zh) | 星敏感器的三轴精度表述与测量方法 | |
CN104567868A (zh) | 基于ins修正的机载长航时天文导航系统的方法 | |
Batista et al. | Constellation design of a lunar global positioning system using cubesats and chip-scale atomic clocks | |
Jenness et al. | Pal: a positional astronomy library | |
Standish et al. | New accuracy levels for Solar System ephemerides | |
Jacobs et al. | The extragalactic and solar system celestial frames: Accuracy, stability, and interconnection | |
Barderas et al. | Observations of Phobos shadow: Analysis of parameters connecting Earth–Mars reference frames | |
De Lafontaine et al. | A review of satellite lifetime and orbit decay prediction | |
Pitjeva | VLBI data are the basis for orientation of planetary ephemerides with respect to ICRF2 and improvement of other ephemeris parameters | |
Hsu et al. | Long-term prediction of GPS satellite orbit |
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 |