EM算法在廣播星歷計算衛(wèi)星位置中的應(yīng)用

肖琴琴，宋迎春，杜 琨

（中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083）

在衛(wèi)星導(dǎo)航定位中，為了確定用戶的位置，必須首先通過衛(wèi)星星歷計算衛(wèi)星的位置，因此，獲取GPS衛(wèi)星的位置是求得接收機(jī)坐標(biāo)的必要環(huán)節(jié)。目前，獲得衛(wèi)星位置的途徑主要有兩種：一種是通過廣播星歷計算得到，一種是通過事后的精密星歷直接獲取。廣播星歷是在接收機(jī)觀測時接收衛(wèi)星發(fā)播的導(dǎo)航電文經(jīng)譯碼取得的實時星歷，而精密星歷是由有關(guān)單位提供的事后處理星歷。其中廣播星歷的精度偏低［1］，約為2m，而由IGS提供的精密星歷精度較高，優(yōu)于5cm，但是精密星歷需觀測工作結(jié)束11d后才可用［2］。GPS廣播星歷雖然精度較精密星歷低，但它具有實時、易獲取的特點，常在動態(tài)實時導(dǎo)航定位中得到應(yīng)用。因此，對實時性要求較高的定位工作而言，只能通過廣播星歷獲得衛(wèi)星位置。

使用廣播星歷直接計算衛(wèi)星位置可以通過文獻(xiàn)［3］中的公式進(jìn)行求解，然而當(dāng)采樣頻率較高需多次計算衛(wèi)星的位置時，直接計算的方法就會占用很多的內(nèi)存并耗費很長時間。考慮到衛(wèi)星位置是一個連續(xù)變化的過程，在一定時間段內(nèi)可將衛(wèi)星坐標(biāo)表示為時間多項式，在內(nèi)存中僅保存多項式的系數(shù)，那么在計算衛(wèi)星坐標(biāo)時，只要調(diào)出多項式系數(shù)即可。這與按固定公式計算衛(wèi)星坐標(biāo)的方法相比，大大提高了數(shù)據(jù)處理的效率。

由于衛(wèi)星遮擋、衛(wèi)星信號傳播限制等原因，可能造成部分有用信息丟失，導(dǎo)致數(shù)據(jù)不完全，或者計算出來的星歷參數(shù)值具有一定的誤差，從而使廣播星歷計算出的衛(wèi)星坐標(biāo)與其真值之間有較大的偏差，另外每計算一個觀測歷元的衛(wèi)星坐標(biāo)也需要較大的計算量，因此如何提高衛(wèi)星坐標(biāo)計算的效率及精度是值得探討的問題。EM（expectation maximization）算法是近年發(fā)展很快且應(yīng)用很廣的一種算法，它能夠在觀測數(shù)據(jù)的基礎(chǔ)上增加一些“潛在數(shù)據(jù)”，從而簡化計算并完成一系列簡單的極大化或模擬，因此可以考慮將EM算法應(yīng)用到衛(wèi)星的軌道標(biāo)準(zhǔn)化。由于不同衛(wèi)星的精度等級不一樣，因此在選擇擬合點的個數(shù)以及擬合階數(shù)時也會有所差別。當(dāng)發(fā)現(xiàn)某些擬合點的誤差太大時，如果直接刪除進(jìn)行降階處理必然導(dǎo)致誤差增大，這時如利用EM算法在觀測數(shù)據(jù)的基礎(chǔ)上加一些“潛在數(shù)據(jù)”，可使得擬合的精度和可靠性仍能滿足要求。

1 常用的插值法和擬合法

利用多項式進(jìn)行插值擬合通常有很多種方法，如拉格朗日多項式插值、埃爾米特插值、普通的多項式擬合、勒讓德多項式擬合和切比雪夫多項式擬合。但目前最常用的多項式是切比雪夫多項式，關(guān)于切比雪夫擬合的原理已在文獻(xiàn)［4］詳細(xì)闡述。

用n階切比雪夫多項式來逼近時間段［t0，t0＋Δt］中的衛(wèi)星星歷時，先要將變量t∈［t0，t0＋Δt］變換為變量τ∈［－1，1］：

于是，衛(wèi)星坐標(biāo)可表示為

式中：n為多項式的階數(shù)，Cxi，Cyi，Czi為切比雪夫多項式的系數(shù)。

切比雪夫多項式Ti的遞推公式為

根據(jù)已知的衛(wèi)星坐標(biāo)，用最小二乘法擬合出多項式系數(shù)Cxi，Cyi，Czi后，就可以計算該時段內(nèi)任意時刻的衛(wèi)星位置。

2 EM算法及其在衛(wèi)星軌道計算中的應(yīng)用

廣播星歷文件經(jīng)常會由于衛(wèi)星信號限制等原因，導(dǎo)致部分?jǐn)?shù)據(jù)失真或缺失，在數(shù)據(jù)缺失或數(shù)據(jù)量較少的情況下，會導(dǎo)致坐標(biāo)的計算精度降低。采用EM算法，添加與衛(wèi)星軌道信息相關(guān)的“潛在數(shù)據(jù)”，可以有效解決這一問題，大大提高衛(wèi)星坐標(biāo)擬合的精度。

2.1 EM算法的原理

EM算法［5］是一種迭代算法，它的每一次迭代由兩步組成：E步（求期望）和M步（極大化）［6］。一般地，以表示θ的基于觀測數(shù)據(jù)的后驗分布密度，稱為觀測數(shù)據(jù)后驗分布；以P（θ/Y，Z）表示添加數(shù)據(jù)（缺失數(shù)據(jù)）Z后得到的關(guān)于θ的后驗分布密度函數(shù)，稱為完全數(shù)據(jù)后驗分布。P（Z/θ，Y）表示在給定θ和觀測數(shù)據(jù)Y下潛在數(shù)據(jù)（缺失數(shù)據(jù)）Z的條件分布密度函數(shù)，其目的是計算觀測后驗分布的參數(shù)。

EM算法的進(jìn)行如式（4）、式（5）所示。記θi為第i＋1次迭代開始時后驗參數(shù)的估計值，則第i＋1次迭代的兩步為：

E步：將P（θ/Y，Z）或logP（θ/Y，Z）關(guān)于Z的條件分布求期望，從而把Z積掉，即

M步：將Q（θ/θi，Y）極 大 化，即 找 一 個 點θi＋1，使

如此形成一次迭代θi→→θi＋1，θi＋1∈M（θi），這里M（θi）是在整個參數(shù)空間Ω內(nèi)使得Q（θ/θi，Y）最大的θ的值所組成的集合。將上述E步和M步進(jìn)行 迭 代 直 至 ‖θi＋1－θi‖ 或 者充分小時為止［7－9］。

2.2 EM算法在廣播星歷計算衛(wèi)星軌道中的應(yīng)用

這里以8階切比雪夫多項式擬合為例，衛(wèi)星的坐標(biāo)可以表示為8階切比雪夫多項式：

則誤差方程為

廣播星歷文件是每隔2h進(jìn)行一次更新，在用多項式擬合其函數(shù)模型時，可先對某個時刻的廣播星歷進(jìn)行外推1h，把某些加密時刻的數(shù)據(jù)當(dāng)成是缺失數(shù)據(jù)，應(yīng)用EM算法進(jìn)行處理，可提高擬合的精度。假設(shè)在式（8）中l(wèi)k為缺失數(shù)據(jù)（加密時刻的衛(wèi)星坐標(biāo)數(shù)據(jù)），已知誤差向量V服從高斯正態(tài)分布，利用EM算法建立似然函數(shù)方程：

缺失數(shù)據(jù)lk的條件分布概率密度函數(shù)為

式中：θi為經(jīng)過i次迭代后的參數(shù)值。

則得到EM算法的期望步：

EM算法的極大化步：將Q（θ/θi，Y）極大化，積分后對各參數(shù)求偏導(dǎo)數(shù)，計算得到參數(shù)θi＋1，將θi＋1代入E步進(jìn)行迭代，循環(huán)直至‖θi＋1－θi‖或者充分小時為止。

3 算例分析

由于GPS廣播星歷的外推時間一般為1h［10］，因此，只對廣播星歷中的某一歷元前、后各1h內(nèi)的衛(wèi)星位置進(jìn)行擬合。本文采用從Internet網(wǎng)上下載的GPS衛(wèi)星2010－09－05的廣播星歷abpo2480.10n數(shù)據(jù)（RINEX格式），對19號衛(wèi)星在1：00－3：00時段內(nèi)X方向上的坐標(biāo)值進(jìn)行擬合，并將擬合出的結(jié)果分別與直接由廣播星歷和IGS精密星歷計算出的軌道進(jìn)行對比。由于內(nèi)插點較多，本文僅選取靠近缺失值附近的8個點值來進(jìn)行對比。計算結(jié)果見表1、表2和圖1、圖2，其中表1與圖1是擬合值與廣播星歷計算值之間的比較結(jié)果，表2與圖2是擬合值與IGS精密星歷的比較結(jié)果。

圖2 與精密星歷計算值的絕對誤差曲線

表2 與IGS擬合結(jié)果精度比較 m

由計算結(jié)果可以看出：

1）當(dāng)用于切比雪夫多項式擬合的數(shù)據(jù)精度較好且數(shù)據(jù)完整的情況下，使用切比雪夫多項式內(nèi)插獲得的衛(wèi)星位置和直接使用廣播星歷用戶算法計算出的結(jié)果非常接近，且其誤差在毫米級以內(nèi)，這充分說明使用切比雪夫多項式來擬合衛(wèi)星軌道是合理的，如表1和圖1所示。

圖1 與廣播星歷計算值的絕對誤差曲線

表1 與廣播星歷直接擬合結(jié)果精度比較 m

2）在缺失數(shù)據(jù)的情況下，使用降階的切比雪夫多項式計算出的衛(wèi)星位置將會產(chǎn)生明顯的偏差，且在缺失點附近尤為顯著，偏差大小達(dá)到數(shù)據(jù)完整時的1000多倍，如表1所示。這也說明，在數(shù)據(jù)缺失時有必要采取一定的措施來提高切比雪夫多項式擬合的精度和可靠性。

3）通過引入EM算法來生成缺失數(shù)據(jù)的潛在值，解決了因數(shù)據(jù)不足需要對切比雪夫多項式進(jìn)行降階處理的問題，且使用生成的潛在數(shù)據(jù)與其它數(shù)據(jù)一起擬合出的衛(wèi)星位置的精度和可靠性與數(shù)據(jù)完整時基本相當(dāng)，如表1和圖1所示。同時也表明，在數(shù)據(jù)缺失的情況下，使用EM算法來生成缺失數(shù)據(jù)的潛在值是可行的，且其精度和可靠性能夠滿足需求。

4）將表2、圖2與表1、圖1進(jìn)行比較可知，無論是將直接使用GPS廣播星歷用戶算法計算出的衛(wèi)星位置作為參考值，還是將IGS精密星歷內(nèi)插出的衛(wèi)星位置作為參考值，其結(jié)論都是相同的，只不過前者的效果尤為顯著，其原因是廣播星歷本身的精度要明顯地低于IGS精密星歷。

4 結(jié)束語

在動態(tài)導(dǎo)航定位中，數(shù)據(jù)的采樣頻率越來越高，如果直接使用廣播星歷用戶算法來計算衛(wèi)星的位置，必然會大大增加導(dǎo)航定位計算的內(nèi)存或負(fù)荷，從而導(dǎo)致實時性難以滿足。此時，選擇較少的采樣點并使用廣播星歷用戶算法計算相應(yīng)時刻的衛(wèi)星位置，再用切比雪夫多項式擬合，將會節(jié)約導(dǎo)航定位計算占用的內(nèi)存并減少計算負(fù)荷。但是，由于各方面因素的影響，直接使用廣播星歷計算出的衛(wèi)星位置個數(shù)及精度會受到一定的影響，導(dǎo)致部分所需節(jié)點的觀測數(shù)據(jù)缺失或不可用，此時必須使用降階的切比雪夫多項式進(jìn)行擬合，其精度和可靠性將會大大降低。本文引入的EM算法能夠生成缺失或不可用節(jié)點位置的潛在值，彌補了因數(shù)據(jù)缺少切比雪夫多項式需要降階處理的不足，且使用EM算法后的切比雪夫多項式不需降階，其擬合結(jié)果的精度和可靠性與數(shù)據(jù)完整時基本相同。因此，EM算法確實是一種解決部分衛(wèi)星軌道數(shù)據(jù)缺失或不可用時獲取任意時刻的衛(wèi)星位置的較好方法。

［1］邱蕾，廖遠(yuǎn)琴，花向紅.基于IGS精密星歷的衛(wèi)星坐標(biāo)插值［J］.測繪工程，2008，17（4）：15－18.

［2］樓益棟，劉萬科，張小紅.GPS衛(wèi)星星歷精度分析［J］.測繪信息與工程，2003，28（6）：4－6.

［3］徐紹銓，張華海，楊志強.GPS測量原理及應(yīng)用［M］.武漢：武漢大學(xué)出版社，2002.

［4］余鵬，孫學(xué)金，趙世軍.GPS定位中衛(wèi)星坐標(biāo)計算的切比雪夫多項式擬合法［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］林東方，宋迎春，肖琴琴.缺失數(shù)據(jù)下基于EM算法的GPS高程擬合［J］.大地測量與地球動力學(xué)，2012，32（2）：1－4.

［8］楊基棟.EM算法理論及其應(yīng)用［J］.安慶師范學(xué)院學(xué)報：自然科學(xué)版，2009，15（4）：30－35

［9］曾傳璜，張鑫，張晶晶，等.EM算法的研究［J］.軟件導(dǎo)刊，2008，7（9）：97－98

［10］羅力，黃聲享.單組廣播星歷精度分析及其衛(wèi)星軌道擬合研究［J］.測繪信息與工程，2010，35（1）：21－22.