EM算法在广播星历计算卫星位置中的应用

肖琴琴，宋迎春，杜 琨

（中南大学 地球科学与信息物理学院，湖南 长沙 410083）

在卫星导航定位中，为了确定用户的位置，必须首先通过卫星星历计算卫星的位置，因此，获取GPS卫星的位置是求得接收机坐标的必要环节。目前，获得卫星位置的途径主要有两种：一种是通过广播星历计算得到，一种是通过事后的精密星历直接获取。广播星历是在接收机观测时接收卫星发播的导航电文经译码取得的实时星历，而精密星历是由有关单位提供的事后处理星历。其中广播星历的精度偏低［1］，约为2m，而由IGS提供的精密星历精度较高，优于5cm，但是精密星历需观测工作结束11d后才可用［2］。GPS广播星历虽然精度较精密星历低，但它具有实时、易获取的特点，常在动态实时导航定位中得到应用。因此，对实时性要求较高的定位工作而言，只能通过广播星历获得卫星位置。

使用广播星历直接计算卫星位置可以通过文献［3］中的公式进行求解，然而当采样频率较高需多次计算卫星的位置时，直接计算的方法就会占用很多的内存并耗费很长时间。考虑到卫星位置是一个连续变化的过程，在一定时间段内可将卫星坐标表示为时间多项式，在内存中仅保存多项式的系数，那么在计算卫星坐标时，只要调出多项式系数即可。这与按固定公式计算卫星坐标的方法相比，大大提高了数据处理的效率。

由于卫星遮挡、卫星信号传播限制等原因，可能造成部分有用信息丢失，导致数据不完全，或者计算出来的星历参数值具有一定的误差，从而使广播星历计算出的卫星坐标与其真值之间有较大的偏差，另外每计算一个观测历元的卫星坐标也需要较大的计算量，因此如何提高卫星坐标计算的效率及精度是值得探讨的问题。EM（expectation maximization）算法是近年发展很快且应用很广的一种算法，它能够在观测数据的基础上增加一些“潜在数据”，从而简化计算并完成一系列简单的极大化或模拟，因此可以考虑将EM算法应用到卫星的轨道标准化。由于不同卫星的精度等级不一样，因此在选择拟合点的个数以及拟合阶数时也会有所差别。当发现某些拟合点的误差太大时，如果直接删除进行降阶处理必然导致误差增大，这时如利用EM算法在观测数据的基础上加一些“潜在数据”，可使得拟合的精度和可靠性仍能满足要求。

1 常用的插值法和拟合法

利用多项式进行插值拟合通常有很多种方法，如拉格朗日多项式插值、埃尔米特插值、普通的多项式拟合、勒让德多项式拟合和切比雪夫多项式拟合。但目前最常用的多项式是切比雪夫多项式，关于切比雪夫拟合的原理已在文献［4］详细阐述。

用n阶切比雪夫多项式来逼近时间段［t0，t0＋Δt］中的卫星星历时，先要将变量t∈［t0，t0＋Δt］变换为变量τ∈［－1，1］：

于是，卫星坐标可表示为

式中：n为多项式的阶数，Cxi，Cyi，Czi为切比雪夫多项式的系数。

切比雪夫多项式Ti的递推公式为

根据已知的卫星坐标，用最小二乘法拟合出多项式系数Cxi，Cyi，Czi后，就可以计算该时段内任意时刻的卫星位置。

2 EM算法及其在卫星轨道计算中的应用

广播星历文件经常会由于卫星信号限制等原因，导致部分数据失真或缺失，在数据缺失或数据量较少的情况下，会导致坐标的计算精度降低。采用EM算法，添加与卫星轨道信息相关的“潜在数据”，可以有效解决这一问题，大大提高卫星坐标拟合的精度。

2.1 EM算法的原理

EM算法［5］是一种迭代算法，它的每一次迭代由两步组成：E步（求期望）和M步（极大化）［6］。一般地，以表示θ的基于观测数据的后验分布密度，称为观测数据后验分布；以P（θ/Y，Z）表示添加数据（缺失数据）Z后得到的关于θ的后验分布密度函数，称为完全数据后验分布。P（Z/θ，Y）表示在给定θ和观测数据Y下潜在数据（缺失数据）Z的条件分布密度函数，其目的是计算观测后验分布的参数。

EM算法的进行如式（4）、式（5）所示。记θi为第i＋1次迭代开始时后验参数的估计值，则第i＋1次迭代的两步为：

E步：将P（θ/Y，Z）或logP（θ/Y，Z）关于Z的条件分布求期望，从而把Z积掉，即

M步：将Q（θ/θi，Y）极 大 化，即 找 一 个 点θi＋1，使

如此形成一次迭代θi→→θi＋1，θi＋1∈M（θi），这里M（θi）是在整个参数空间Ω内使得Q（θ/θi，Y）最大的θ的值所组成的集合。将上述E步和M步进行 迭 代 直 至 ‖θi＋1－θi‖ 或 者充分小时为止［7－9］。

2.2 EM算法在广播星历计算卫星轨道中的应用

这里以8阶切比雪夫多项式拟合为例，卫星的坐标可以表示为8阶切比雪夫多项式：

则误差方程为

广播星历文件是每隔2h进行一次更新，在用多项式拟合其函数模型时，可先对某个时刻的广播星历进行外推1h，把某些加密时刻的数据当成是缺失数据，应用EM算法进行处理，可提高拟合的精度。假设在式（8）中lk为缺失数据（加密时刻的卫星坐标数据），已知误差向量V服从高斯正态分布，利用EM算法建立似然函数方程：

缺失数据lk的条件分布概率密度函数为

式中：θi为经过i次迭代后的参数值。

则得到EM算法的期望步：

EM算法的极大化步：将Q（θ/θi，Y）极大化，积分后对各参数求偏导数，计算得到参数θi＋1，将θi＋1代入E步进行迭代，循环直至‖θi＋1－θi‖或者充分小时为止。

3 算例分析

由于GPS广播星历的外推时间一般为1h［10］，因此，只对广播星历中的某一历元前、后各1h内的卫星位置进行拟合。本文采用从Internet网上下载的GPS卫星2010－09－05的广播星历abpo2480.10n数据（RINEX格式），对19号卫星在1：00－3：00时段内X方向上的坐标值进行拟合，并将拟合出的结果分别与直接由广播星历和IGS精密星历计算出的轨道进行对比。由于内插点较多，本文仅选取靠近缺失值附近的8个点值来进行对比。计算结果见表1、表2和图1、图2，其中表1与图1是拟合值与广播星历计算值之间的比较结果，表2与图2是拟合值与IGS精密星历的比较结果。

图2 与精密星历计算值的绝对误差曲线

表2 与IGS拟合结果精度比较 m

由计算结果可以看出：

1）当用于切比雪夫多项式拟合的数据精度较好且数据完整的情况下，使用切比雪夫多项式内插获得的卫星位置和直接使用广播星历用户算法计算出的结果非常接近，且其误差在毫米级以内，这充分说明使用切比雪夫多项式来拟合卫星轨道是合理的，如表1和图1所示。

图1 与广播星历计算值的绝对误差曲线

表1 与广播星历直接拟合结果精度比较 m

2）在缺失数据的情况下，使用降阶的切比雪夫多项式计算出的卫星位置将会产生明显的偏差，且在缺失点附近尤为显著，偏差大小达到数据完整时的1000多倍，如表1所示。这也说明，在数据缺失时有必要采取一定的措施来提高切比雪夫多项式拟合的精度和可靠性。

3）通过引入EM算法来生成缺失数据的潜在值，解决了因数据不足需要对切比雪夫多项式进行降阶处理的问题，且使用生成的潜在数据与其它数据一起拟合出的卫星位置的精度和可靠性与数据完整时基本相当，如表1和图1所示。同时也表明，在数据缺失的情况下，使用EM算法来生成缺失数据的潜在值是可行的，且其精度和可靠性能够满足需求。

4）将表2、图2与表1、图1进行比较可知，无论是将直接使用GPS广播星历用户算法计算出的卫星位置作为参考值，还是将IGS精密星历内插出的卫星位置作为参考值，其结论都是相同的，只不过前者的效果尤为显著，其原因是广播星历本身的精度要明显地低于IGS精密星历。

4 结束语

在动态导航定位中，数据的采样频率越来越高，如果直接使用广播星历用户算法来计算卫星的位置，必然会大大增加导航定位计算的内存或负荷，从而导致实时性难以满足。此时，选择较少的采样点并使用广播星历用户算法计算相应时刻的卫星位置，再用切比雪夫多项式拟合，将会节约导航定位计算占用的内存并减少计算负荷。但是，由于各方面因素的影响，直接使用广播星历计算出的卫星位置个数及精度会受到一定的影响，导致部分所需节点的观测数据缺失或不可用，此时必须使用降阶的切比雪夫多项式进行拟合，其精度和可靠性将会大大降低。本文引入的EM算法能够生成缺失或不可用节点位置的潜在值，弥补了因数据缺少切比雪夫多项式需要降阶处理的不足，且使用EM算法后的切比雪夫多项式不需降阶，其拟合结果的精度和可靠性与数据完整时基本相同。因此，EM算法确实是一种解决部分卫星轨道数据缺失或不可用时获取任意时刻的卫星位置的较好方法。

［1］邱蕾，廖远琴，花向红.基于IGS精密星历的卫星坐标插值［J］.测绘工程，2008，17（4）：15－18.

［2］楼益栋，刘万科，张小红.GPS卫星星历精度分析［J］.测绘信息与工程，2003，28（6）：4－6.

［3］徐绍铨，张华海，杨志强.GPS测量原理及应用［M］.武汉：武汉大学出版社，2002.

［4］余鹏，孙学金，赵世军.GPS定位中卫星坐标计算的切比雪夫多项式拟合法［J］.气象科技，2004，32（3）：198－201.

［5］DEMPSTER，A.P，LAIRD，N.M.，and Rubin，D.B.Maximum Likelihood From Incomplete Data Via the EM Algorithm［J］.Journal of the Royal Statistifcal Society B，1977，39（1）：1－38.

［6］GRAHAM C.G，JUAN C.A，Approximate EM Algorithms for Parameter and State Estimation in Nonlinear Stochastic Models［J］.Proceedings of the 44th IEEE Conference on Decision and Control，and the European Control Conference，2005：368－373

［7］林东方，宋迎春，肖琴琴.缺失数据下基于EM算法的GPS高程拟合［J］.大地测量与地球动力学，2012，32（2）：1－4.

［8］杨基栋.EM算法理论及其应用［J］.安庆师范学院学报：自然科学版，2009，15（4）：30－35

［9］曾传璜，张鑫，张晶晶，等.EM算法的研究［J］.软件导刊，2008，7（9）：97－98

［10］罗力，黄声享.单组广播星历精度分析及其卫星轨道拟合研究［J］.测绘信息与工程，2010，35（1）：21－22.