肖琴琴，宋迎春，杜 琨

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

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

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

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

1 常用的插值法和擬合法

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

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

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

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

切比雪夫多項(xiàng)式Ti的遞推公式為

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

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

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

2.1 EM算法的原理

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

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

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

M步：將Q（θ/θi，Y）極 大 化，即 找 一 個(gè) 點(diǎn)θi＋1，使

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

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

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

則誤差方程為

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

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

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

則得到EM算法的期望步：

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

3 算例分析

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

圖2 與精密星歷計(jì)算值的絕對(duì)誤差曲線(xiàn)

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

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

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

圖1 與廣播星歷計(jì)算值的絕對(duì)誤差曲線(xiàn)

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

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

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

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

4 結(jié)束語(yǔ)

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

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

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

［3］徐紹銓?zhuān)瑥埲A海，楊志強(qiáng).GPS測(cè)量原理及應(yīng)用［M］.武漢：武漢大學(xué)出版社，2002.

［4］余鵬，孫學(xué)金，趙世軍.GPS定位中衛(wèi)星坐標(biāo)計(jì)算的切比雪夫多項(xiàng)式擬合法［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］.大地測(cè)量與地球動(dòng)力學(xué)，2012，32（2）：1－4.

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

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

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