基于双频非差无电离层延迟组合观测量的导航卫星(GPS)精密星历计算算法与实践

苏贝

(哈尔滨市勘察测绘研究院，黑龙江 哈尔滨 150010)

1 引 言

精密定轨理论的研究主要开始于20 世纪50年代末，经过半个多世纪的发展已日趋完善，定轨精度随着卫星跟踪技术的发展，测量精度的提高以及轨道动力学模型的演变和精细化得到了明显的提升［1］。

在这之中，IGS 扮演了重要的角色，它汇集了世界知名的数据分析中心，包括JPL(Jet Propulsion Laboratory，美国喷气推进实验室)、CODE(Center for Orbit Determination in Europe，欧洲定轨中心)、GFZ(Geo Forschungs Zentrum，德国地学研究中心)、MIT(Massachusetts Institute of Technology，美国麻省理工学院)、SIO(Scripps Institution of Oceanography，美国斯克里普斯海洋学研究所)、USNO(U．S．Naval Observatory，美国海军天文台)、NOAA(National Oceanic and Atmospheric Administration，美国国家海洋和大气管理局)等，计算并发布精密星历，为地球科学研究、多学科应用、商业开发及教育提供支撑。在此基础上，各家数据分析中心形成了业界知名的卫星导航定位数据分析软件，如JPL 的GIPSY，CODE 的BERNESE，MIT 的GAMIT 等。

2 基本算法

从解微分方程的角度来看，精密定轨就是将一个常微分方程初值问题转化为边值问题，由边值条件来确定初值［2］。这些微分方程通常也需要应用数值方法来解算。

2.1 观测模型

导航卫星精密定轨的基本观测量是L1 和L2 载波相位非差无电离层延迟线性组合观测量电离层延迟可抵消掉。

基本观测方程如下:

式中，λ1、λ2为L1/L2 载 波 波 长(m);f1、f2为L1/L2载波频率(Hz);φL1、φL2为L1/L2 载波相位观测值(周);ρ 为站星几何距离(m);dt、dT为接收机和卫星的时钟钟差(s);Vtron为信号对流层路径延迟(m);NLC为无电离层延迟线性组合的载波相位偏差(m);△pcvs、△pcvr为卫星和接收机的天线相位中心变化(m);△rel为相对论效应改正(m);△phw为相位缠绕效应改正(m);εLC为无电离层延迟线性组合测量噪声(m)。

式(1)里的站星几何距离ρ 通过解光时方程得到的，该方程考虑了卫星运动和接收机钟差导致的误差，为:

式中，rs(t)为卫星的空间位置(m);r'r为观测站的空间位置(m);U(t)为惯性系到地固系的转换矩阵;为卫星和接收机的天线位置改正(m);为测站位移(m)。

地固坐标通过转换矩阵可以转换成惯性系坐标，考虑岁差章动及地球自转，转换矩阵为:

式中，N(t)、P(t)为岁差章动矩阵;Rx、Ry、Rz分别为绕x、y、z 轴的旋转矩阵;GAST 为格林尼治真恒星时;GMST 为格林尼治平恒星时;xp、yp为地极偏移;UT1－UTC 为地球自转角偏;△ψ、ε 为章动经线分量和倾角;r 为世界时和恒星时之比。

卫星位置rs须使用统一的参考时间t 来表达，通常是采用GPS 时。因信号传播时间与接收机钟差之和△t 通常不超过0.1 s，时间同步方程可简化为式(8)，该方程仅考虑卫星速度和地球点质量效应。同时，式(9)为φLC相对于时刻t 卫星位置rs的偏导数。

式中，GME为地球引力系数;vs(t)为t 时刻卫星s的速度(m/s);Φrr(t，t0)为t0时刻到t 时刻的卫星状态转移矩阵。

应用映射函数NMF，对流层延迟Vtron表示如下:

式中，ZTD 为天顶总对流层延迟(m);ZHD 为天顶对流层静力学延迟(m);GE、GN为对流层大气水平梯度参数;Mdry、Mwet为映射函数的干、湿分量;Az、El 为卫星的方位角、高度角;P0为平均海水面处的大气压(hPa);φ、H 为测站的纬度、大地高(m)。

卫星天线相位中线变化表示如下:

式中，Esat→eci为星固坐标系到惯性系的转换;△satao为以星固坐标表示的天线偏移;θ 为相对于地面站的天底角。

接收机天线相位中心变化表示如下:

式中，Elocal→ecef为测站坐标系到地固坐标系的转换;△ecc为测站标志到天线参考点的距离(m);△apcr1、△apcr2为天线参考点到L1、L2 载波相位中心的偏差(m);△pcvr1、△pcvr2为L1、L2 载波相位中心偏差(m)。由地球潮汐引起的测站位移表示如下:

式中，△solid为地球固体潮引起的位移(m);△ocean为海洋负荷引起的位移(m);△polar为极潮引起的位移(m)。

相对论效应和相位缠绕效应表示如下:

式中，D、D'为卫星、接收机的有效偶极向量;x，y，z为本地接收机单位向量(东、北、高度);x'，y'，z'为卫星体坐标单位向量;k 为卫星至接收机的单位向量。

表1归纳了观测模型中精密测量改正的相关信息。

表1 精密测量改正项

2.2 卫星轨道模型

一般而言，卫星轨道的高度决定其所受扰动的程度。就精密定轨而言，任何大于10－10m/s2量级的加速度都需要考虑进来。因此，计算GPS 卫星轨道时，模型必须考虑到至少8 阶重力位影响及相关潮汐效应，还有日月引力、太阳辐射压和相对论效应。

卫星轨道方程如下式所示，为卫星位置和速度的常微分方程。其中，状态转移矩阵采用方程数值积分方式解算。表2列出了计算时所采用的一些精密卫星轨道模型。

式中，ageop为重力位加速度改正(m/s2);a3rdbody为三体重力加速度(m/s2);asrp为太阳辐射压加速度(m/s2);arel为由相对论效应产生的加速度(m/s2)。

表2 精密卫星轨道模型

2.3 其他状态转移模型

其他的状态转移模型如下:卫星钟差采用一阶高斯马尔科夫;接收机钟为白噪声;对流层参数为随机游走;地球自转参数为随机游走。在一个弧段内，载波相位偏差认为是一个固定值。如果在该弧段探测到了周跳，则将其重新初始化。

2.4 参数估计

参数估计采用扩展卡尔曼滤波的方法。下列方程组分别表示观测更新和时序更新。除了前向滤波外，计算时还可采用后向滤波及其更平滑的方法。

对于GPS 卫星定轨而言，卡尔曼滤波器所估计的状态向量，包括卫星位置和速度、卫星时钟、接收机时钟、对流层参数。此外，如必要，卫星的太阳辐射压和地球自转参数也能纳入进来。在估计中，卫星和接收机的时钟当作相对于参考时刻而变化的时钟来处理，并固定参考时刻为0。

输入观测向量zk后，tk时刻的状态向量^xk通过卡尔曼滤波估计过程逐步得到。

3 算例

3.1 HRBCORS

哈尔滨双星导航服务系统(GNSS)由哈尔滨市勘察测绘研究院于2007年建成，一期包含5 个连续运行参考站点，分别是阿城、哈西、对青、方台、双城。每站配备一台双频72 通道接收机，型号为NetR5，可同时接收GPS 和GLONASS 卫星信号，是本次计算的数据来源，如图1所示。

图1 哈尔滨市连续运行参考站分布图

3.2 轨道精度

在上面叙述的情况下，以开源的GAMIT 软件代码为基础，添加了非差相位观测方程解算模块，修改了部分参数估计算法，得到了GPS 卫星的定轨结果，然后与IGS 综合精密星历比较，得出结论。后者由IGS 组织应用全球范围内的GPS 站网，对所有分析中心给出的结果加权平均得到，其中每家分析中心应用不同的精密分析软件包进行数据分析，由此得到的轨道精度一般能够达到5 cm以内。

表3列出了卫星轨道的估计精度，以IGS 综合精密星历为标准，给出了所有卫星位置均方跟误差的平均值。图2比较了所有卫星的轨道估计误差。

表3 卫星轨道精度

图2 各卫星的轨道估计误差

4 结 语

本文叙述了一种GPS 卫星精密定轨算法，结合实测数据计算了卫星星历，与IGS 综合精密星历相比较，轨道外符合精度略大于5 cm，即卫星位置的三维均方差约为5.2 cm。作者下一步可继续进行精密钟差估计的工作，为区域非差分用户提供精密定位和授时服务。

［1］刘林．航天器轨道理论［M］．北京:国防工业出版社，2000．

［2］施闯，赵齐乐，李敏等．北斗卫星导航系统的精密定轨与定位研究［J］．中国科学:地球科学，2012，42:845～861．

［3］郑作亚，黄城，卢秀山．星载GPS 精密定轨进展及其数学模型［J］．大地测量与地球动力学，2007，27:112～118．

［4］周建华，杨龙，徐波等．一种导航卫星中长期轨道预报方法［J］．测绘学报，2011，40(S):39～45．

［5］李敏．多模GNSS 融合精密定轨理论及其应用研究［D］．武汉:武汉大学，2011．

［6］韩保民．基于星载GPS 的低轨卫星几何法定轨理论研究［D］．武汉:中科院测量与地球物理研究所，2003．

［7］叶世榕．GPS 非差相位精密单点定位理论与实现［D］．武汉:武汉大学，2002．

［8］Oliver M．etc．．Satellite Orbits:Models，Methods and Applications［M］．Springer－Verlag，Heidelberg，2000．