CN109738919B - 一种用于gps接收机自主预测星历的方法 - Google Patents

Abstract

本发明公开了一种用于GPS接收机自主预测星历的方法，包括：步骤1：根据输入的广播星历计算星历有效时间内卫星的位置速度矢量，作为精密定轨的观测值；步骤2：根据观测值进行精密定轨，利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量，以及太阳光压模型参数中的卫星面质比；步骤3：进行轨道预报，以定轨弧段初始时刻的的精确位置速度为起始，采用卫星面质比作为卫星等效面质比参数进行轨道外推，外推到开机时刻；步骤4：将开机时刻的位置速度量拟合成广播星历，给出预报星历，将预报星历直接应用于接收机位置的定位解算。极大地缩小了非热启动时的接收机首次定位时长，降低对芯片存储空间的要求，有利于芯片的成本控制。

Description

一种用于GPS接收机自主预测星历的方法

技术领域

本发明属于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(t₀,X₀；t) (4)

其中：X₀为待估状态量；

将式(2)在参考状态量

并舍去二阶和二阶以上的小量，得到：

其中：参考状态量

为待估状态量X₀的近似值；

Y(X*,t)是观测量的理论计算值，是测量矩阵，记为

是状态转移矩阵，记为Φ(t₀,t)；

令：

得到精密定轨基本方程：

y＝Bx₀ (8)

其中：y为残差；

通过迭代方式由观测采样数据求解精密定轨基本方程(8)，得到待估状态量X₀的改正值

给出历元t₀时刻达到精度要求的状态量

则：

其中：下标k表示引用了k次采样数据；重复迭代i次直至满足精度要求，得到：

步骤2中：

精密定轨基本方程(8)的求解采用最小二乘估计，

通过式(11)计算得到：

其中：

表示从t₀到t_l时刻的状态转移矩阵，H_l表示t_l时刻的测量矩阵，y_l表示t_l时刻的残差向量；

矩阵B通过式(12)计算：

测量矩阵H_X通过式(13)计算：

其中：x,y,z表示三个位置分量；

状态转移矩阵Φ(t₀,t)通过状态方程(4)得到，状态转移矩阵Φ(t₀,t)通过式(14)计算：

步骤2中得到状态转移矩阵Φ(t₀,t)的具体方法为：

令：

X(t)＝Φ(t_,t₀)X₀ (15)

则：

状态转移矩阵Φ(t₀,t)所满足的条件通过对状态微分方程(3)的线性化得到，有：

其中：状态微分方程F为：

记：

则方程(17)的系数矩阵A：

将方程(17)转换为一阶方程组：

其中：

C₁(t),C₂(t),C₃(t)均用参考值X^*进行计算；

对式(21)进行积分得到

通过

得到状态转移矩阵Φ(t₀,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表示地球半径，C_nm和S_nm表示未归一化的重力势系数，根据式(27)把归一化后的重力势系数

和

转化成未归一化前的C_nm和S_nm：

V_nm和W_nm的计算公式如下：

已知：

利用式(29)置m＝0获得带谐项V_n0，相应的值均等于0，递归使用式(31)：

式(31)中

过程采用式(28)，↓过程采用式(29)，计算得到所有的V_nm和W_nm，

加速度

等于势能u的梯度，直接由V_nm和W_nm计算得出：

其中：

地球带谐摄动加速度、田谐摄动加速度对位置的偏导数根据：

通过式(40)计算得到：

地球带谐摄动加速度、田谐摄动加速度对速度的偏导数：

步骤3中日月引力引起的卫星加速度的具体确定方法为：

在以地球为中心的参考系下，太阳和月球的摄动表示为：

其中：r表示卫星位置，s表示太阳或月亮的位置；

低精度日月坐标具体计算方法如下：

T＝(MJD-51544.5)/36525.0 (43)

摄动的计算基于月球平黄经L₀、月球平近点角l、太阳平近点角l′、太阳平黄经和月球平黄经之间的差D以及月球平升交距角F，5个基础角；

L₀＝(0.606433+1336.851344T-floor(0.606433+1336.851344T)) (44)

记相对于2000年黄道和春分点的月球黄经为L_M，利用式(45)计算值，有：

L_M＝2π(L₀+dL/1296.0e3-floor(L₀+dL/1296.0e3)) (47)

月球黄纬B_M表示为：

月球的地心距r_M表示为：

月球坐标r_M为：

其中：ε为黄赤交角，ε＝23°.43929111；

假设地球围绕太阳作无摄运动，通过式(51)描述太阳相对于地球和黄道在2000年附近几十年的椭圆轨道：

位置坐标利用式(51)由开普勒轨道公式得到；

太阳的黄经L_⊙和距离r_sun如下：

通过旋转转换成赤道坐标系的笛卡尔坐标：

由于黄纬B_⊙在1′精度下为0，将式(53)简化为：

其中：r_sun表示太阳位置；

日月摄动加速度对位置的偏导数：

其中：r表示卫星位置，s表示太阳或月亮的位置；

日月摄动加速度对速度的偏导数：

步骤3中太阳光压引起的卫星加速度的具体确定方法为：

通过式(57)计算太阳光压加速度：

其中：ν为蚀因子；P_⊙为地球附近太阳光压强参数，其值约为4.56×10-6N/m2；光压系数C_R＝1+ε，卫星制造中典型材料的反射系数ε为0.2～0.9；A/m为卫星的面质比；r_⊙和r_⊙分别是卫星到太阳的向量及其模；AU为天文单位表示距离，1AU＝149597870.691km；

太阳光压加速度对卫星位置的偏导数：

其中：r表示卫星位置，r_sun表示太阳位置；

太阳光压加速度对卫星速度的偏导数：

太阳光压对卫星等效面质比的偏导数：

同时在两个坐标系中表述的光压摄动力，一个是基于地心的BYD坐标系，另一个是基于卫星本体的XYZ坐标系；

其中：r_sun、r为惯性系下太阳和卫星位置，e_D定义为太阳至卫星方向的单位矢量，e_x,e_y,e_z为星固坐标系坐标轴的单位矢量，其中：e_z指向地球中心，e_y＝e_z×e_D，e_x＝e_y×e_z，e_B＝e_D×e_y；

u定义为卫星在轨道平面上距升交点的角度，u₀代表的是太阳的升交角距；

D(β)＝λ(SRP(1)·D₀+D_C2cos(2β)+D_C4cos(4β)) (62)

其中：D_C2＝-0.813×10^-9m/s²，D_C4＝0.517×10^-9m/s²，β定义为太阳相对于卫星轨道平面的高度角D₀为ROCK模型计算出来的太阳辐射压产生的加速度的理论值，SRP为待估参数；

Y(β)＝SRP(2)·D₀+Y_Ccos(2β) (63)

其中：Y_C＝0.067×10^-9m/s²；

B(β)＝SRP(3)·D₀+B_Ccos(2β) (64)

其中：B_C＝-0.385×10^-9m/s²；

X₁(β)＝SRP(4)·D₀+X₁₀+X_1Ccos(2β)+X_1Ssin(2β) (65)

其中：X₁₀＝-0.015×10^-9m/s²，X_1C＝-0.018×10^-9m/s²，X_1S＝-0.033×10^-9m/s²

X₃(β)＝SRP(5)·D₀+X₃₀+X_3Ccos(2β)+X_3Ssin(2β) (66)

其中：X₃₀＝0.004×10^-9m/s²，X_3C＝-0.046×10^-9m/s²，X_3S＝-0.398×10^-9m/s²；

z(β)＝SRP(6)·D₀+z₀+z_C2cos(2β)+z_S2sin(2β)+z_c4cos(4β)+z_S4cos(4β) (67)

其中：Z₀(BlockII)＝1.024×10^-9m/s²，Z₀(BlockIIA)＝0.979×10^-9m/s²,Z_C2＝0.519×10^-9m/s²，Z_S2＝0.125×10^-9m/s²，

Z_C4＝0.047×10^-9m/s²，Z_S4＝-0.045×10^-9m/s²。

步骤3采用Adams-Cowell积分法求解式(22)进行轨道预报的预报校正方法为：

9-1：采用Adams显示公式计算速度，Cowell显示公式计算位置；

9-2：预报：

9-3：根据预报的t_n+1时刻的y_n+1，

求解相应的t_n+1处右函数值f_n+1；

9-4：校正：

9-5：回到9-2，重复9-2～9-4，至求得所有时刻的卫星位置和速度。

步骤4的具体方法为：

10-1：轨道根数拟合；具体为：

计算卫星位置需用到16个星历参数，其中星历表参考历元t_oe是已知的，需要拟合15个参数；待估状态参数向量和观测方程为：

γ＝γ(X_,t) (69)

其中：X为参考历元时刻的星历参数，Y为一个含m(m≥15)个观测量的观测列向量；

设X_i为参数X在第i次迭代的初值，将观测方程在所给初值处展开，并舍去二阶和二阶以上的小量后得：

其中：Y(X_i,t)为用参考历元t_oe时刻广播星历参数初值计算的卫星位置，

分别为相应星历参数的改正值，

为观测量对星历参数的偏导数；令：

L＝Y-Y(X_i,t) (72)

得误差方程：

V＝HδX_i-L (73)

由最小二乘原理有：

δX_i＝(H^TH)^-1H^TL (74)

则第i次迭代后的星历参数估值为：

X＝X_i+δX_i (75)

在程序中，所选用的迭代结束条件为：

D_i-D_i-1＜1e-8

其中：k是定轨所用的历元数，n是待估状态参数X的维数，下标i表第i次迭代；

10-2：求解卫星位置对星历参数的偏导数：

根据卫星位置对星历参数的偏导数计算公式：

给出用来推导GPS16参数偏导数的用户位置的最终计算表达式：

其中：

卫星位置对轨道长半径的平方根的偏导数，

卫星位置对轨道偏心率的偏导数，

卫星位置对平近点角的偏导数，

卫星位置对平均角速度的偏导数，

卫星位置对升交点赤经的偏导数，

卫星位置对升交点赤经变化率的偏导数，

卫星位置对轨道倾角的偏导数，

卫星位置对轨道倾角变化率的偏导数，

卫星位置对近地点角距的偏导数，

卫星位置对升交点余弦调和项校正的偏导数，

卫星位置对升交点正弦调和项校正的偏导数，

卫星位置对轨道半径余弦调和项校正的偏导数，

卫星位置对轨道半径正弦调和项校正的偏导数，

卫星位置对轨道倾角余弦调和项校正的偏导数，

卫星位置对轨道倾角正弦调和项校正的偏导数；

10-3：选取轨道根数拟合法初值；

以二体问题的方法求解初始6个密切开普勒轨道根数，其余9个摄动调和参数初值设为基于两组位置矢量的位置确定：

r_a,r_b取某一固定惯性坐标系下的位置，在程序中取周内秒为0时刻的ECEF坐标系的三轴方向为固定坐标系作为自定义惯性坐标系，计算卫星在坐标系下t_a,t_b时刻的位置；

其中：θ₁等于地球球自转角速度与t_a时刻对应的周内秒的成乘积，θ₂等于地球自转角速度与t_b时刻对应的周内秒的乘积；

高斯矢量：

W＝e_a×e₀

其中：

e_a为r_a的单位矢量，e₀为r₀的单位矢量，e_a与e₀相互垂直；

倾角i和升交点赤经Ω如下给出：

纬度角：

将根据两组位置矢量和时间间隔来确定轨道的问题转化为同求解轨道和向径组成的扇形和三角形的比例问题；

令r_a和r_b表示卫星最在时间t_a和t_b的地心距；由矢量r_a和r_b围成的三角形的面积Δ依赖于边长r_a和r_b以及两矢量之间的夹角v_b-v_a，假定该角度小于180°；三角形面积为：

其中：v_a和v_b为不同时刻的真近点角；

扇形面积S是指矢量r_a和r_b以及轨道弧线围城的区域，根据开普勒第二定律此面积大小与t_a和t_b之间的时间差成正比：

其中：a和e为轨道的长半轴和轨道的偏心率；代入半通径p＝a(1-e²)，得：

其中，η是扇形与三角形面积之比，归一化的时间间隔τ定义如下：

采用二体问题方程式的已知量来替代半通径，将η转换成包含两个方程的系统：

其中：正、辅助变量为：

由式(82)得到l的最小值为

根据g确定η，g的值等于t_b时刻偏近点角与t_a时刻偏近点角的差值的1/2；消掉g得到超越方程：

其中：W函数如下定义：

变量ω对椭圆轨道始终为正，并且小于1位迭代确定η，使用割线法：

求解式(86)的根：

选取初始值：

η₁＝η₀+0.1 η₂＝η₀ (87)

由Hansen近似法给出：

解得η后，根据式(79)，半通径用时间间隔τ表示如下：

其中：由r_a和r_b确定的三角形面积如下：

根据圆锥曲线方程得出轨道偏心率，求解e·cosν时，有：

由于：

得到两个方程：

求解t_a时刻的偏心率和真近点角：

近地点辐角是纬度角和真近点角之差：

ω＝u_a-v_a (95)

根据半通径和偏心率，得到长半轴：

最后，由开普勒方程计算平近点角M_a：

M_a＝E_a-e·sinE_a (rad) (97)

偏近点角E_a满足方程：

10-4：将预测出的卫星地心地固坐标系进行拟合，得到与GPS接口控制文件结构一致的预测广播星历。

与现有技术相比，本发明具有以下有益效果：

通过精密定轨，确定定轨时刻的精确位置速度，并以定轨时刻的精确位置速度为起始，进行轨道外推，外推到开机时刻，将得到的开机时刻附近的位置速度量拟合成广播星历，得到与GPS接口控制文件结构一致的预测广播星历，此预测广播星历可直接应用于接收机位置的解算，无需等待广播星历的收集，大幅缩短非热启动时的首次定位时长；在关机超过2小时再次开机时，首次定位时长从30秒减少至在10秒以内。同时，预测广播星历有效期为4小时，占用空间远小于存储4小时的卫星地心地固坐标，降低对芯片存储空间的要求，有利于芯片的成本控制。

附图说明

图1为本发明的方法流程框图。

具体实施方式

为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

下面结合附图对本发明做进一步详细描述：

参见图1，本发明一种用于GPS接收机自主预测星历的方法，包括以下步骤：

步骤1：根据输入的广播星历计算星历有效时间内卫星的位置速度矢量，作为精密定轨的观测值；其中，广播星历的输入可以通过以下方式：

对于每个卫星，从导航电文中获得广播星历至GPS接收机，但是不以此为限。

步骤2：根据观测值进行精密定轨，利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量，以及太阳光压模型参数中的卫星面质比。

定轨计算实际是利用含有误差的测量数据，结合航天器受到的各种外力作用，获取航天器状态的最佳估计值。由于在航天定轨的同时确定一些轨道相关的动力学或测量上的物理参数和几何参数，又称为精密定轨——轨道确定与参数估计。

本定轨程序中，首先，根据输入的一组或多组数据星历，计算得到一段时间(t₀,t₁…t_k)的卫星位置/速度观测量。然后，利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量，以及太阳光压模型参数中的卫星面质比Am2。该面质比值将作为常数用于后面的轨道预报中。

具体方法为：建立航天器的状态量模型：

其中：r,v表示t时刻卫星的位置速度矢量，β表示待估参数卫星等效面质比；

观测量Y与航天器的状态量X的量测方程：

Y＝Y(X_,t) (2)

其中：Y表示观测量的测量值；X(t)满足轨道动力学方程：

求解式(3)得到状态方程：

X(t)＝X(t₀,X₀；t) (4)

其中：X₀为待估状态量；

将式(2)在参考状态量

并舍去二阶和二阶以上的小量，得到：

其中：参考状态量

为待估状态量X₀的近似值；

Y(X^*,t)是观测量的理论计算值，是测量矩阵，记为H_X；

是状态转移矩阵，记为Φ(t₀,t)；

令：

得到精密定轨基本方程：

y＝Bx₀ (8)

其中：y为残差；

通过迭代方式根据观测采样数据Y_j(j＝1,…k)求解精密定轨基本方程(8)，关机前保存的一组有效期星历，通过按照GPS接口规范文档可以计算出星历有效期期间任意时刻的卫星位置作为观测采样数据；得到待估状态量X₀的改正值

给出历元t₀时刻达到精度要求的状态量

则：

(9)

其中：下标k表示引用了k次采样数据；重复迭代i次直至满足精度要求，得到：

其中：精密定轨基本方程(8)的求解采用最小二乘估计，

通过式(11)计算得到：

B_l＝H_lΦ_l,0

其中：

表示从t₀到t_l时刻的状态转移矩阵，H_l表示t_l时刻的测量矩阵，y_l表示t_l时刻的残差向量；

矩阵B通过式(12)计算：

测量矩阵H_X通过式(13)计算：

其中：x,y,z表示三个位置分量；

状态转移矩阵Φ(t₀,t)通过状态方程(4)得到，状态转移矩阵Φ(t₀,t)通过式(14)计算：

得到状态转移矩阵Φ(t₀,t)的具体方法为：令：X(t)＝Φ(t_,t₀)X₀ (15)

则：

状态转移矩阵Φ(t₀,t)所满足的条件通过对状态微分方程(3)的线性化得到，有：

其中：状态微分方程F为：

记：

则方程(17)的系数矩阵A：

将方程(17)转换为一阶方程组：

其中：

C₁(t),C₂(t),C₃(t)均用参考值X^*进行计算；对式(21)进行积分得到

通过

得到状态转移矩阵Φ(t₀,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表示地球半径，C_nm和S_nm表示未归一化的重力势系数，根据式(27)把归一化后的重力势系数

和

转化成未归一化前的C_nm和S_nm：

V_nm和W_nm的计算公式如下：

已知：

利用式(29)置m＝0获得带谐项V_n0，相应的值W_n0均等于0，递归使用式(31)：

(31)中

过程采用式(28)，↓过程采用式(29)，计算得到所有的V_nm和W_nm，

加速度

等于势能u的梯度，直接由V_nm和W_nm计算得出：

其中：

地球带谐摄动加速度、田谐摄动加速度对位置的偏导数根据：

通过式(40)计算得到：

地球带谐摄动加速度、田谐摄动加速度对速度的偏导数：

其中：日月引力引起的卫星加速度的具体确定方法为：

在以地球为中心的参考系下，太阳和月球的摄动表示为：

其中：r表示卫星位置，s表示太阳或月亮的位置；低精度日月坐标具体计算方法如下：

T＝(MJD-51544.5)/36525.0 (43)

摄动的计算基于月球平黄经L₀、月球平近点角l、太阳平近点角l′、太阳平黄经和月球平黄经之间的差D以及月球平升交距角F，5个基础角；

L₀＝(0.606433+1336.851344T-floor(0.606433+1336.851344T)) (44)

记相对于2000年黄道和春分点的月球黄经为L_M，利用式(45)计算值，有：

L_M＝2π(L₀+dL/1296.0e3-floor(L₀+dL/1296.0e3)) (47)

月球黄纬B_M表示为：

月球的地心距r_M表示为：

月球坐标r_M为：

其中：ε为黄赤交角，ε＝23°.43929111；

假设地球围绕太阳作无摄运动，通过式(51)描述太阳相对于地球和黄道在2000年附近几十年的椭圆轨道：

位置坐标利用式(51)由开普勒轨道公式得到；

太阳的黄经L_⊙和距离r_sun如下：

通过旋转转换成赤道坐标系的笛卡尔坐标：

由于黄纬B_⊙在1′精度下为0，将式(53)简化为：

其中：r_sun表示太阳位置；

日月摄动加速度对位置的偏导数：

其中：r表示卫星位置，s表示太阳或月亮的位置；

日月摄动加速度对速度的偏导数：

其中：太阳光压引起的卫星加速度的具体确定方法为：

通过式(57)计算太阳光压加速度：

其中：ν为蚀因子；P_⊙为地球附近太阳光压强参数，其值约为4.56×10-6N/m2；光压系数C_R＝1+ε，卫星制造中典型材料的反射系数ε为0.2～0.9；A/m为卫星的面质比；r_⊙和r_⊙分别是卫星到太阳的向量及其模；AU为天文单位表示距离，1AU＝149597870.691km；

太阳光压加速度对卫星位置的偏导数：

其中：r表示卫星位置，r_sun表示太阳位置；

太阳光压加速度对卫星速度的偏导数：

太阳光压对卫星等效面质比的偏导数：

同时在两个坐标系中表述的光压摄动力，一个是基于地心的BYD坐标系，另一个是基于卫星本体的XYZ坐标系；

其中：r_sun、r为惯性系下太阳和卫星位置，e_D定义为太阳至卫星方向的单位矢量，e_x,e_y,e_z为星固坐标系坐标轴的单位矢量，其中：e_z指向地球中心，e_y＝e_z×e_D，e_x＝e_y×e_z，e_B＝e_D×e_y；

u定义为卫星在轨道平面上距升交点的角度，u₀代表的是太阳的升交角距；

D(β)＝λ(SRP(1)·D₀+D_C2cos(2β)+D_C4cos(4β)) (62)

其中：D_C2＝-0.813×10^-9m/s²，D_C4＝0.517×10^-9m/s²，β定义为太阳相对于卫星轨道平面的高度角，D₀为ROCK模型计算出来的太阳辐射压产生的加速度的理论值，SRP为待估参数；

Y(β)＝SRP(2)·D₀+γ_Ccos(2β) (63)

其中：γ_C＝0.067×10^-9m/s²；

B(β)＝SRP(3)·D₀+B_Ccos(2β) (64)

其中：B_C＝-0.385×10^-9m/s²；

X₁(β)＝SRP(4)·D₀+X₁₀+X_1Ccos(2β)+X_1Ssin(2β) (65)

其中：X₁₀＝-0.015×10^-9m/s²，X_1C＝-0.018×10^-9m/s²，X_1S＝-0.033×10^-9m/s²

X₃(β)＝SRP(5)·D₀+X₃₀+X_3Ccos(2β)+X_3Ssin(2β) (66)

其中：X₃₀＝0.004×10^-9m/s²，X_3C＝-0.046×10^-9m/s²，X_3S＝-0.398×10^-9m/s²；

z(β)＝SRP(6)·D₀+z₀+z_C2cos(2β)+z_S2sin(2β)+z_c4cos(4β)+z_S4cos(4β) (67)

其中：Z₀(BlockII)＝1.024×10^-9m/s²，Z₀(BlockIIA)＝0.979×10^-9m/s²,Z_C2＝0.519×10^-9m/s²，Z_S2＝0.125×10^-9m/s²，Z_C4＝0.047×10^-9m/s²，Z_S4＝-0.045×10^-9m/s²。

其中采用Adams-Cowell积分法求解式(22)进行轨道预报的具体内容如下：

Adams-Cowell是多步法，计算速度较快，但多步法需要用到前两个以上时刻的计算结果，因此它不是自起步的。在实际应用中，多步法都是需要用单步法来进行起步的，本发明采用单步的RK8方法来起步。A-C方法是Adams方法和Cowell方法的结合，因此在介绍AC法前，首先给出Adams和Cowell这两种算法的计算过程。

1、Adams算法：适用于一阶常微分方程形如：

显示公式：

其中：

表示二项式系数，即

隐式公式：

2、Cowell算法

Cowell方法可对二阶微分方程直接积分。二阶常微分方程如：

可以看出，该方程的右函数中不含

(对卫星轨道运动方程来说，就是右函数不含速度项，仅考虑引力作用时的运动方程就是如此)，因此在每一步计算中，只需直接给出y_n即可，而不用去计算

这样采用Cowell方法较为方便。

Cowell显示公式：

其中：

Cowell隐式公式：

且

3、Adams-Cowell算法

一般来说，Cowell公式比同阶的Adams公式精度要高。对于k阶Adams公式，截断误差为O(h^k+1)，而对于k阶Cowell公式截断误差为O(h^k+2)。但是在卫星运动方程中，常常是要出现

也就是速度项的。对于初值问题：

可以将Adams方法和Cowell方法联合起来使用，由Adams公式提供

(速度矢量)，而用Cowell提供y(位置矢量)，这将比单用Adams公式更加有效。另外，隐式方法比显示方法精度高，且数值稳定性好，但是需要一些迭代过程来求解，计算比较复杂。因此在求解卫星运动方程中，常常使用预估—校正(PECE)算法，将显示公式和隐式公式结合起来使用。P即预报，用显示公式将解从t_n预报到t_n+1点处；E即估计，用t_n+1点处预报值计算相应的右函数值；C即改正，将计算出的右函数值通过隐式公式得到改正的预报值；最后一个E仍是估计，但是是用修正后的t_n+1点处预报值再去计算右函数，并以此作为下一步积分的起始值。

预报校正方法为：

9-1：采用Adams显示公式计算速度，Cowell显示公式计算位置；

9-2：预报：

9-3：根据预报的t_n+1时刻的y_n+1，

求解相应的t_n+1处右函数值f_n+1；

9-4：校正：

9-5：回到9-2，重复9-2～9-4，至求得所有时刻的卫星位置和速度。

步骤4：将开机时刻的位置速度量拟合成广播星历，给出预报星历，将预报星历直接应用于接收机位置的定位解算，具体方法为：

10-1：轨道根数拟合；具体为：

计算卫星位置需用到16个星历参数，其中星历表参考历元t_oe是已知的，需要拟合15个参数；待估状态参数向量和观测方程为：

Y＝Y(X_,t) (69)

其中：X为参考历元时刻的星历参数，Y为一个含m(m≥15)个观测量的观测列向量；

设X_i为参数X在第i次迭代的初值，将观测方程在所给初值处展开，并舍去二阶和二阶以上的小量后得：

其中：Y(X_i,t)为用参考历元t_oe时刻广播星历参数初值计算的卫星位置，

分别为相应星历参数的改正值，

为观测量对星历参数的偏导数；令：

L＝Y-γ(X_i,t) (72)

可得误差方程：V＝HδX_i-L；

由最小二乘原理有：

δX_i＝(H^TH)^-1H^TL (74)

则第i次迭代后的星历参数估值为：

X＝X_i+δX_i (75)

在程序中，所选用的迭代结束条件为：

其中：k是定轨所用的历元数，n是待估状态参数X的维数，下标i表第i次迭代；

10-2：求解卫星位置对星历参数的偏导数：

根据卫星位置对星历参数的偏导数计算公式：

给出用来推导GPS16参数偏导数的用户位置的最终计算表达式：

其中：

卫星位置对轨道长半径的平方根的偏导数，

卫星位置对轨道偏心率的偏导数，

卫星位置对平近点角的偏导数，

卫星位置对平均角速度的偏导数，

卫星位置对升交点赤经的偏导数，

卫星位置对升交点赤经变化率的偏导数，

卫星位置对轨道倾角的偏导数，

卫星位置对轨道倾角变化率的偏导数，

卫星位置对近地点角距的偏导数，

卫星位置对升交点余弦调和项校正的偏导数，

卫星位置对升交点正弦调和项校正的偏导数，

卫星位置对轨道半径余弦调和项校正的偏导数，

卫星位置对轨道半径正弦调和项校正的偏导数，

卫星位置对轨道倾角余弦调和项校正的偏导数，

卫星位置对轨道倾角正弦调和项校正的偏导数；

10-3：选取轨道根数拟合法初值；

以二体问题的方法求解初始6个密切开普勒轨道根数，其余9个摄动调和参数初值设为基于两组位置矢量的位置确定：

r_a,r_b取某一固定惯性坐标系下的位置，在程序中取周内秒为0时刻的ECEF坐标系的三轴方向为固定坐标系作为自定义惯性坐标系，计算卫星在坐标系下t_a,t_b时刻的位置；

其中：θ₁等于地球球自转角速度与t_a时刻对应的周内秒的成乘积，θ₂等于地球自转角速度与t_b时刻对应的周内秒的乘积；

高斯矢量：

W＝e_a×e₀

其中：

e_a为r_a的单位矢量，e₀为r₀的单位矢量，e_a与e₀相互垂直；

倾角i和升交点赤经Ω如下给出：

纬度角：

将根据两组位置矢量和时间间隔来确定轨道的问题转化为同求解轨道和向径组成的扇形和三角形的比例问题；

令r_a和r_b表示卫星最在时间t_a和t_b的地心距；由矢量r_a和r_b围成的三角形的面积Δ依赖于边长r_a和r_b以及两矢量之间的夹角v_b-v_a，假定该角度小于180°；三角形面积为：

其中：v_a和v_b为不同时刻的真近点角；

扇形面积S是指矢量r_a和r_b以及轨道弧线围城的区域，根据开普勒第二定律此面积大小与t_a和t_b之间的时间差成正比：

其中：a和e为轨道的长半轴和轨道的偏心率；代入半通径p＝a(1-e²)，得：

其中，η是扇形与三角形面积之比，归一化的时间间隔τ定义为：

采用二体问题方程式的已知量来替代半通径，将η转换成包含两个方程的系统：

其中：正、辅助变量为：

由式(82)得到l的最小值为

根据g确定η，g的值等于t_b时刻偏近点角与t_a时刻偏近点角的差值的1/2；消掉g得到超越方程：

其中：W函数如下定义：

变量ω对椭圆轨道始终为正，并且小于1位迭代确定η，使用割线法：

求解式(86)的根：

选取初始值：

η₁＝η₀+0.1 η₂＝η₀ (87)

由Hansen近似法给出：

解得η后，根据式(79)，半通径用时间间隔τ表示如下：

其中：由r_a和r_b确定的三角形面积为：

根据圆锥曲线方程得出轨道偏心率，求解e·cosν时，有：

由于：

得到两个方程：

求解t_a时刻的偏心率和真近点角：

近地点辐角是纬度角和真近点角之差：ω＝u_a-v_a(95)

根据半通径和偏心率，得到长半轴：

最后，由开普勒方程计算平近点角M_a：M_a＝E_a-e·sinE_a (rad) (97)

偏近点角E_a满足方程：

10-4：将预测出的卫星地心地固坐标系进行拟合，得到与GPS接口控制文件结构一致的预测广播星历。

仿真实验表明，在关机超过2小时再次开机时，没有星历预测首次定位时长为30秒以上，采用本方法首次定位时长在10秒以内。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (10)

1.一种用于GPS接收机自主预测星历的方法，其特征在于，包括以下步骤：

步骤1：根据输入的广播星历计算星历有效时间内卫星的位置速度矢量，作为精密定轨的观测值；

步骤2：根据所述观测值进行精密定轨，利用定轨程序计算得到定轨弧段初始时刻的精确位置速度矢量，以及太阳光压模型参数中的卫星面质比；

步骤3：进行轨道预报，以所述定轨弧段初始时刻的精确位置速度为起始，采用所述卫星面质比作为卫星等效面质比参数进行轨道外推，外推到开机时刻，得到开机时刻的位置速度矢量；

步骤4：将所述开机时刻的位置速度矢量拟合成广播星历，直接应用于接收机的定位解算。

2.根据权利要求1所述的用于GPS接收机自主预测星历的方法，其特征在于，所述步骤2中进行精密定轨的具体方法为：

建立航天器的状态量模型：

其中：r,v分别表示t时刻卫星的位置矢量和速度矢量，β表示待估参数卫星等效面质比；

观测量与航天器的状态量X的量测方程：

Y＝Y(X,t) (2)

其中：Y表示观测量的测量值且只包含位置矢量；X(t)满足轨道动力学方程：

求解式(3)得到状态方程：

X(t)＝X(t₀,X₀；t) (4)

其中：X₀为待估状态量；

将式(2)在参考状态量

并舍去二阶和二阶以上的小量，得到：

其中：参考状态量

为待估状态量X₀的近似值；

Y(X^*,t)是观测量的理论计算值，是测量矩阵，记为H_X；

是状态转移矩阵，记为Φ(t₀,t)；

令：

得到精密定轨基本方程：

y＝Bx₀ (8)

其中：y为残差；

通过迭代方式由观测采样数据求解精密定轨基本方程(8)，得到待估状态量X₀的改正值

给出历元t₀时刻达到精度要求的状态量

则：

其中：下标k表示引用了k次采样数据；重复迭代_i次直至满足精度要求，得到：

3.根据权利要求2所述的用于GPS接收机自主预测星历的方法，其特征在于，所述步骤2中：

精密定轨基本方程(8)的求解采用最小二乘估计，

通过式(11)计算得到：

其中：

表示从t₀到t_l时刻的状态转移矩阵，H_l表示t_l时刻的测量矩阵，y_l表示t_l时刻的残差向量；

矩阵B通过式(12)计算：

测量矩阵H_X通过式(13)计算：

其中：x,y,z表示三个位置分量；

状态转移矩阵Φ(t₀,t)通过状态方程(4)得到，状态转移矩阵Φ(t₀,t)通过式(14)计算：

4.根据权利要求3所述的用于GPS接收机自主预测星历的方法，其特征在于，所述步骤2中得到状态转移矩阵Φ(t₀,t)的具体方法为：

令：

X(t)＝Φ(t_,t₀)X₀ (15)

则：

状态转移矩阵Φ(t₀,t)所满足的条件通过对状态微分方程(3)的线性化得到，有：

其中：状态微分方程F为：

记：

则方程(17)的系数矩阵A：

将方程(17)转换为一阶方程组：

其中：

C₁(t),C₂(t),C₃(t)均用参考值X^*进行计算；

对式(21)进行积分得到z,

通过z,

得到状态转移矩阵Φ(t₀,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表示地球半径，C_nm和S_nm表示未归一化的重力势系数，根据式(27)把归一化后的重力势系数

和

转化成未归一化前的C_nm和S_nm：

V_nm和W_nm的计算公式如下：

已知：

利用式(29)置m＝0获得带谐项V_n0，相应的值均等于0，递归使用式(31)：

式(31)中

过程采用式(28)，↓过程采用式(29)，计算得到所有的V_nm和W_nm，

加速度

等于势能u的梯度，直接由V_nm和W_nm计算得出：

其中：

地球带谐摄动加速度、田谐摄动加速度对位置的偏导数根据：

通过式(40)计算得到：

地球带谐摄动加速度、田谐摄动加速度对速度的偏导数：

7.根据权利要求5所述的用于GPS接收机自主预测星历的方法，其特征在于，所述步骤3中日月引力引起的卫星加速度的具体确定方法为：

在以地球为中心的参考系下，太阳和月球的摄动表示为：

其中：r表示卫星位置，s表示太阳或月亮位置，当s为太阳位置r_sun时，得到太阳的摄动，当s为月球位置r_moon时，得到月球的摄动；

低精度日月坐标具体计算方法如下：

T＝(MJD-51544.5)/36525.0 (43)

摄动的计算基于月球平黄经L₀、月球平近点角l、太阳平近点角l′、太阳平黄经和月球平黄经之间的差D以及月球平升交距角F，5个基础角；

L₀＝(0.606433+1336.851344T-floor(0.606433+1336.851344T)) (44)

记相对于2000年黄道和春分点的月球黄经为L_M，利用式(45)计算值，有：

L_M＝2π(L₀+dL/1296.0e3-floor(L₀+dL/1296.0e3)) (47)

月球黄纬B_M表示为：

月球的地心距r_M表示为：

月球坐标r_M为：

其中：ε为黄赤交角，ε＝23.43929111°；

假设地球围绕太阳作无摄运动，通过式(51)描述太阳相对于地球和黄道在2000年附近几十年的椭圆轨道：

位置坐标利用式(51)由开普勒轨道公式得到；

太阳的黄经L_⊙和位置r_sun如下：

通过旋转转换成赤道坐标系的笛卡尔坐标：

由于黄纬B_⊙在1′精度下为0，将式(53)简化为：

其中：r_sun表示太阳位置；

日月摄动加速度对位置的偏导数：

其中：r表示卫星位置，s表示太阳或月亮的位置；

日月摄动加速度对速度的偏导数：

8.根据权利要求5所述的用于GPS接收机自主预测星历的方法，其特征在于，所述步骤3中太阳光压引起的卫星加速度的具体确定方法为：

通过式(57)计算太阳光压加速度：

其中：ν为蚀因子；P_⊙为地球附近太阳光压强参数，其值约为4.56×10-6N/m2；光压系数C_R＝1+ε_reflect，卫星制造中典型材料的反射系数ε_reflect为0.2～0.9；A/m为卫星的面质比；r_⊙和r_⊙分别是卫星到太阳的向量及其模；AU为天文单位表示距离，1AU＝149597870.691km；

太阳光压加速度对卫星位置的偏导数：

其中：r表示卫星位置，r_sun表示太阳位置；

太阳光压加速度对卫星速度的偏导数：

太阳光压对卫星等效面质比的偏导数：

同时在两个坐标系中表述的光压摄动力，一个是基于地心的BYD坐标系，另一个是基于卫星本体的XYZ坐标系；

其中：r_sun、r为太阳和卫星位置，e_D定义为太阳至卫星方向的单位矢量，e_x,e_y,e_z为星固坐标系坐标轴的单位矢量，其中：e_z指向地球中心，e_y＝e_z×e_D，e_x＝e_y×e_z，e_B＝e_D×e_y；

u定义为卫星在轨道平面上距升交点的角度，u₀代表的是太阳的升交角距；

D(ξ)＝λ(SRP(1)·D₀+D_C2cos(2ξ)+D_C4cos(4ξ)) (62)

其中：D_C2＝-0.813×10^-9m/s²，D_C4＝0.517×10^-9m/s²，βξ定义为太阳相对于卫星轨道平面的高度角，D₀为ROCK模型计算出来的太阳辐射压产生的加速度的理论值，SRP为待估参数；

Y(ξ)＝SRP(2)·D₀+Y_Ccos(2ξ) (63)

其中：Y_C＝0.067×10^-9m/s²；

B(ξ)＝SRP(3)·D₀+B_Ccos(2ξ) (64)

其中：B_C＝-0.385×10^-9m/s²；

X₁(ξ)＝SRP(4)·D₀+X₁₀+X_1Ccos(2ξ)+X_1Ssin(2ξ) (65)

其中：X₁₀＝-0.015×10^-9m/s²，X_1C＝-0.018×10^-9m/s²，X_1S＝-0.033×10^-9m/s²

X₃(ξ)＝SRP(5)·D₀+X₃₀+X_3Ccos(2ξ)+X_3Ssin(2ξ) (66)

其中：X₃₀＝0.004×10^-9m/s²，X_3C＝-0.046×10^-9m/s²，X_3S＝-0.398×10^-9m/s²；

z(ξ)＝SRP(6)·D₀+z₀+z_C2cos(2ξ)+z_S2sin(2ξ)+z_c4cos(4ξ)+z_S4cos(4ξ) (67)

其中：Z₀(BlockII)＝1.024×10^-9m/s²，Z₀(BlockIIA)＝0.979×10^-9m/s²,Z_C2＝0.519×10^-9m/s²，Z_S2＝0.125×10^-9m/s²，Z_C4＝0.047×10^-9m/s²，Z_S4＝-0.045×10^-9m/s²。

9.根据权利要求5所述的用于GPS接收机自主预测星历的方法，其特征在于，所述步骤3采用Adams-Cowell积分法求解式(22)进行轨道预报的预报校正方法为：

9-1：采用Adams显示公式计算速度，Cowell显示公式计算位置；

9-2：预报：

9-3：根据预报的t_n+1时刻的y_n+1，

求解相应的t_n+1处右函数值f_n+1；

9-4：校正：

9-5：回到9-2，重复9-2～9-4，至求得所有时刻的卫星位置和速度。

10.根据权利要求1所述的用于GPS接收机自主预测星历的方法，其特征在于，所述步骤4的具体方法为：

10-1：轨道根数拟合；具体为：

计算卫星位置需用到16个星历参数，其中星历表参考历元t_oe是已知的，需要拟合15个参数；待估状态参数向量和观测方程为：

Y＝Y(X_,t) (69)

其中：X为参考历元时刻的星历参数，Y为一个含m，m≥15个观测量的观测列向量；

设X_i为参数X在第i次迭代的初值，将观测方程在所给初值处展开，并舍去二阶和二阶以上的小量后得：

其中：Y(X_i,t)为用参考历元t_oe时刻广播星历参数初值计算的卫星位置，

分别为相应星历参数的改正值，

为观测量对星历参数的偏导数；令：

L＝Y-Y(X_i,t) (72)

得误差方程：

V＝HδX_i-L (73)

由最小二乘原理有：

δX_i＝(H^TH)^-1H^TL (74)

则第i次迭代后的星历参数估值为：

X＝X_i+δX_i (75)

在程序中，所选用的迭代结束条件为：

其中：k是定轨所用的历元数，n是待估状态参数X的维数，下标i表第i次迭代；

10-2：求解卫星位置对星历参数的偏导数：

根据卫星位置对星历参数的偏导数计算公式：

给出用来推导GPS16参数偏导数的用户位置的最终计算表达式：

其中：

卫星位置对轨道长半径的平方根的偏导数，

卫星位置对轨道偏心率的偏导数，

卫星位置对平近点角的偏导数，

卫星位置对平均角速度的偏导数，

卫星位置对升交点赤经的偏导数，

卫星位置对升交点赤经变化率的偏导数，

卫星位置对轨道倾角的偏导数，

卫星位置对轨道倾角变化率的偏导数，

卫星位置对近地点角距的偏导数，

卫星位置对升交点余弦调和项校正的偏导数，

卫星位置对升交点正弦调和项校正的偏导数，

卫星位置对轨道半径余弦调和项校正的偏导数，

卫星位置对轨道半径正弦调和项校正的偏导数，

卫星位置对轨道倾角余弦调和项校正的偏导数，

卫星位置对轨道倾角正弦调和项校正的偏导数；

10-3：选取轨道根数拟合法初值；

以二体问题的方法求解初始6个密切开普勒轨道根数，其余9个摄动调和参数初值设为基于两组位置矢量的位置确定：

r_a,r_b取某一固定惯性坐标系下的位置，在程序中取周内秒为0时刻的ECEF坐标系的三轴方向为固定坐标系作为自定义惯性坐标系，计算卫星在坐标系下t_a,t_b时刻的位置；

其中：θ₁等于地球球自转角速度与t_a时刻对应的周内秒的成乘积，θ₂等于地球自转角速度与t_b时刻对应的周内秒的乘积；

高斯矢量：

W＝e_a×e₀

其中：

r₀＝r_b-(r_b·e_a)e_a，

e_a为r_a的单位矢量，e₀为r₀的单位矢量，e_a与e₀相互垂直；

倾角i和升交点赤经Ω如下给出：

纬度角：

将根据两组位置矢量和时间间隔来确定轨道的问题转化为同求解轨道和向径组成的扇形和三角形的比例问题；

令r_a和r_b表示卫星最在时间t_a和t_b的地心距；由矢量r_a和r_b围成的三角形的面积Δ依赖于边长r_a和r_b以及两矢量之间的夹角v_b-v_a，假定该角度小于180°；三角形面积为：

其中：v_a和v_b为不同时刻的真近点角；

扇形面积S是指矢量r_a和r_b以及轨道弧线围城的区域，根据开普勒第二定律此面积大小与t_a和t_b之间的时间差成正比：

其中：a和e为轨道的长半轴和轨道的偏心率；代入半通径p＝a(1-e²)，得：

其中，η是扇形与三角形面积之比，归一化的时间间隔τ定义如下：

采用二体问题方程式的已知量来替代半通径，将η转换成包含两个方程的系统：

其中：正、辅助变量为：

由式(82)得到l的最小值为

根据g确定η，g的值等于t_b时刻偏近点角与t_a时刻偏近点角的差值的1/2；消掉g得到超越方程：

其中：W函数如下定义：

变量ω对椭圆轨道始终为正，并且小于1位迭代确定η，使用割线法：

求解式(86)的根：

选取初始值：

η₁＝η₀+0.1 η₂＝η₀ (87)

由Hansen近似法给出：

解得_η后，根据式(79)，半通径用时间间隔τ表示如下：

其中：由r_a和r_b确定的三角形面积如下：

根据圆锥曲线方程得出轨道偏心率，求解e·cosν时，有：

由于：

得到两个方程：

求解t_a时刻的偏心率和真近点角：

近地点辐角是纬度角和真近点角之差：

ω＝u_a-v_a (95)

根据半通径和偏心率，得到长半轴：

最后，由开普勒方程计算平近点角M_a：

M_a＝E_a-e·sinE_a，rad (97)

偏近点角E_a满足方程：

10-4：将预测出的卫星地心地固坐标系进行拟合，得到与GPS接口控制文件结构一致的预测广播星历。

Images