杜 琨，宋迎春，肖琴琴

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

0 引 言

利用GPS信號進(jìn)行導(dǎo)航定位時(shí)，需要已知GPS衛(wèi)星在空間的瞬時(shí)位置，衛(wèi)星星歷的精度將直接影響定位的精度和結(jié)果，因此，精確計(jì)算GPS衛(wèi)星在空間的瞬時(shí)位置是精密定位的前提和基礎(chǔ)［1－2]。隨著現(xiàn)代測量、導(dǎo)航、氣象等領(lǐng)域?qū)πl(wèi)星星歷的精度和實(shí)時(shí)性的要求越來越高，對利用GPS廣播星歷和精密星歷來準(zhǔn)確計(jì)算衛(wèi)星在空間的瞬時(shí)位置，變得越來越重要。目前，已有很多文獻(xiàn)研究了利用精密星歷文件計(jì)算任意時(shí)刻的衛(wèi)星空間位置，且取得了很好的效果，然而這些研究工作都是基于數(shù)據(jù)完全，即星歷文件數(shù)據(jù)準(zhǔn)確、信息量較大無缺失數(shù)據(jù)時(shí)進(jìn)行研究的［3－4]。在實(shí)際應(yīng)用中，精密星歷文件并不完全準(zhǔn)確，由于衛(wèi)星遮擋、衛(wèi)星信號傳播限制等多方面的原因，可能造成部分有用信息丟失，導(dǎo)致數(shù)據(jù)不完全。在這種情況下，如果仍采用常規(guī)的計(jì)算方法，則可能受到數(shù)據(jù)不完全的影響。對于信息量不足，通常只能降階對其擬合或插值，這兩種情況都會造成很大的誤差，影響擬合或插值的效果。已有研究表明，對于15min歷元間隔的精密星歷而言，如果要精確至10－8，用8階朗格朗日內(nèi)插已足夠保證精度，而切比雪夫多項(xiàng)式擬合至少需要12階才能達(dá)到厘米級精度［3]。利用切比雪夫多項(xiàng)式擬合要達(dá)到厘米級精度，至少需要12組數(shù)據(jù)才能對其進(jìn)行擬合。但由于精密星歷是每15min一組，數(shù)據(jù)量相對較少，所以對于短弧段用切比雪夫多項(xiàng)式擬合不能滿足要求。EM算法是一種常用的缺失數(shù)據(jù)處理方法，目前，已被廣泛應(yīng)用于測量領(lǐng)域。2010年，惠沈盈等把EM算法應(yīng)用于數(shù)據(jù)缺失時(shí)的動態(tài)定位解算［5]，2011年，林東方等又把EM算法應(yīng)用于數(shù)據(jù)缺失時(shí)的GPS高程擬合［6]，都取得了很好的效果。通過對EM算法的研究，對不完全數(shù)據(jù)添加與衛(wèi)星軌道信息相關(guān)的“潛在數(shù)據(jù)”，能大大提高衛(wèi)星擬合的精度。

1 常用的插值法和擬合法

精密星歷文件是以離散的形式、按一定的時(shí)間間隔（通常為15min）給出衛(wèi)星在空間的三維坐標(biāo)、三維改正速度以及衛(wèi)星鐘改正數(shù)等信息。而在GPS事后數(shù)據(jù)處理中，經(jīng)常要用到任意時(shí)刻的衛(wèi)星坐標(biāo)，因此需要對精密星歷進(jìn)行插值或者擬合以獲得采樣間隔更高的衛(wèi)星軌道信息。利用精密星歷文件計(jì)算衛(wèi)星位置以及運(yùn)動速度通常采用拉格朗日插值和切比雪夫擬合。

關(guān)于拉格朗日插值和切比雪夫擬合的原理已經(jīng)在大量文獻(xiàn)中有過說明，僅以切比雪夫多項(xiàng)式擬合為例。用n階切比雪夫多項(xiàng)式來逼近時(shí)間段［t0，t0＋Δt]中的衛(wèi)星星歷時(shí)，先要將變量t∈［t0，t0＋Δt]變換為變量τ∈［－1，1]：

式中：n為多項(xiàng)式的階數(shù)；Cxi，Cyi，Czi為切比雪夫多項(xiàng)式Ti的系數(shù)。切比雪夫多項(xiàng)式Ti的遞推公式如下：

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

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

IGS精密星歷文件經(jīng)常因?yàn)橛?jì)算失誤、文件傳播等原因?qū)е虏糠謹(jǐn)?shù)據(jù)失真或缺失。在數(shù)據(jù)缺失或數(shù)據(jù)量較少的情況下，無法采用高階的插值法或擬合法。這時(shí)，只能降階對其插值或擬合，從而導(dǎo)致計(jì)算精度降低。采用EM算法，添加與衛(wèi)星軌道信息相關(guān)的“潛在數(shù)據(jù)”，可以有效解決這一問題，大大提高衛(wèi)星擬合的精度。

2.1 EM算法的原理

EM算法是一種迭代算法。它的每一次迭代由兩步組成：E步（求期望）和 M步（極大化）。一般，以P（θ／Y）表示θ的基于觀測數(shù)據(jù)的后驗(yàn)分布密度，稱為觀測數(shù)據(jù)后驗(yàn)分布。以P（θ／Y，Z）表示添加數(shù)據(jù)（缺失數(shù)據(jù)）Z后得到的關(guān)于θ的后驗(yàn)分布密度函數(shù)，稱為完全數(shù)據(jù)后驗(yàn)分布。P（Z／θ，Y）表示在給定θ和觀測數(shù)據(jù)Y下潛在數(shù)據(jù)（缺失數(shù)據(jù)）Z的條件分布密度函數(shù)。我們的目的是計(jì)算觀測后驗(yàn)分布P（θ／Y）的參數(shù)，EM 算法如下進(jìn)行，記θi為第i＋1次迭代開始時(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‖或者‖（Q（θi＋1／θi，Y）－Ω（θ／Qi，Y）‖充分小時(shí)為止［5－7]。

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

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

則誤差方程為

式中：V為誤差向量；C為擬合系數(shù)；T表示時(shí)間參數(shù)切比雪夫遞推公式；l為觀測向量（衛(wèi)星坐標(biāo)）。可以表示為

精密星歷文件是以離散的形式、按一定的時(shí)間間隔（通常為15min）給出衛(wèi)星在空間的三維坐標(biāo)、三維改正速度以及衛(wèi)星鐘改正數(shù)等信息。在用多項(xiàng)式擬合其函數(shù)模型時(shí)，可以先對精密星歷進(jìn)行加密，而把加密時(shí)刻的數(shù)據(jù)當(dāng)成是缺失數(shù)據(jù)，應(yīng)用EM算法進(jìn)行處理，從而提高擬合的精度。假設(shè)在上式中l(wèi)n為缺失數(shù)據(jù)（加密時(shí)刻的衛(wèi)星坐標(biāo)數(shù)據(jù)），已知誤差向量V服從高斯正態(tài)分布，利用EM算法建立似然函數(shù)方程P（θ／Y，Z），

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

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

則得到EM算法的期望步：

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

3 算例分析

本例選取的數(shù)據(jù)是從IGS數(shù)據(jù)中心下載的2012年1月29日的精密星歷。為了使擬合結(jié)果達(dá)到厘米級精度，這里選取1號衛(wèi)星0時(shí)0分0秒每隔30min共12個(gè)歷元的X方向的坐標(biāo)值作為擬合節(jié)點(diǎn)的數(shù)據(jù)，而將15min間隔的數(shù)據(jù)作為已知數(shù)據(jù)進(jìn)行對比。

分別采用12階切比雪夫多項(xiàng)式擬合、同時(shí)假設(shè)擬合節(jié)點(diǎn)中第五個(gè)歷元的數(shù)據(jù)為粗差或者丟失，在這種情況下降階采用11階切比雪夫多項(xiàng)式擬合以及采用EM算法處理后，再采用12階切比雪夫多項(xiàng)式擬合，計(jì)算擬合弧段內(nèi)每15min歷元節(jié)點(diǎn)的結(jié)果如表1所示。

表1 擬合結(jié)果精度比較（單位：m）

在缺失數(shù)據(jù)下，降階應(yīng)用11階擬合與采用EM算法處理后，再采用12階切比雪夫擬合，其絕對誤差曲線如圖1所示。

圖1 絕對誤差曲線

從圖1可以看出，由于4號和5號節(jié)點(diǎn)間的歷元數(shù)據(jù)缺失，從而導(dǎo)致4號節(jié)點(diǎn)的擬合結(jié)果產(chǎn)生明顯的振蕩。而5號節(jié)點(diǎn)處于整個(gè)擬合弧段的中間，因此擬合結(jié)果較好。

4 結(jié) 論

由于精密星歷是每15min一組，數(shù)據(jù)量相對較少，所以對于短弧段用切比雪夫多項(xiàng)式擬合不能滿足要求。當(dāng)不能采用12階切比雪夫多項(xiàng)式擬合時(shí)，只能采用降階擬合的方法，即采用11階或者更低階數(shù)的切比雪夫多項(xiàng)式來求取任意時(shí)刻的衛(wèi)星空間坐標(biāo)，此時(shí)，可能導(dǎo)致擬合精度無法滿足要求，采用EM算法，通過對其添加與計(jì)算衛(wèi)星位置有關(guān)的“潛在數(shù)據(jù)”，能有效解決這一問題，大大提高擬合的精度。由實(shí)驗(yàn)結(jié)果可以看出，在用切比雪夫多項(xiàng)式擬合內(nèi)插擬合弧段內(nèi)任意時(shí)刻的衛(wèi)星位置時(shí)，擬合弧段兩端節(jié)點(diǎn)的擬合效果較差，使用EM算法也出現(xiàn)了同樣的問題。在數(shù)據(jù)缺失的情況下，與缺失數(shù)據(jù)相鄰的節(jié)點(diǎn)擬合精度要相對較差。

［1]徐紹銓，張華海，楊志強(qiáng)，等.GPS測量原理及應(yīng)用［M].武漢：武漢大學(xué)出版社，2007：11－86.

［2]李征航，張小紅.衛(wèi)星導(dǎo)航定位新技術(shù)及高精度數(shù)據(jù)處理方法［M].武漢：武漢大學(xué)出版社，2009：1－24.

［3]茍長龍.廣播星歷插值和精密星歷外推方法研究［D].長沙：中南大學(xué)，2009：18－22.

［4]萬亞豪，張書畢，侯東陽.GPS精密星歷插值法與擬合法的精度比較［J].全球定位系統(tǒng)，2011，36（2）：1－4.

［5]惠沈盈.觀測數(shù)據(jù)不完全的動態(tài)定位算法研究［D].長沙：中南大學(xué)，2010：12－14.

［6]林東方，宋迎春，肖琴琴.缺失數(shù)據(jù)下基于EM算法的GPS高程擬合［J].大地測量與地球動力學(xué)，2012，32（2）：1－4.

［7]楊基棟.EM算法理論及其應(yīng)用［J].安慶師范學(xué)院學(xué)報(bào)·自然科學(xué)版，2009：1－5.