閆志闖，徐君毅，李 巖，余春平

1. 信息工程大學(xué)，河南 鄭州，450001；2.測繪信息技術(shù)總站，陜西 西安，710054；3.第一測繪導(dǎo)航基地， 遼寧 大連，116031

?

基于配置積分器的GPS衛(wèi)星精密星歷軌道擬合

閆志闖1,2，徐君毅2，李 巖3，余春平2

1. 信息工程大學(xué)，河南 鄭州，450001；2.測繪信息技術(shù)總站，陜西 西安，710054；3.第一測繪導(dǎo)航基地， 遼寧 大連，116031

本文詳細(xì)推導(dǎo)了配置積分器的基本原理，并針對衛(wèi)星進(jìn)入地影或月影時(shí)配置積分器需要變步長的問題，提出了將變步長配置積分器轉(zhuǎn)換為定步長的改進(jìn)方法；同時(shí)，利用CODE精密星歷作為虛擬觀測量，進(jìn)行軌道擬合實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明：與GPS精密星歷相比，擬合后的軌道在R、T和N三個(gè)方向RMS均在0.01m以內(nèi)。

配置積分器；軌道擬合；軌道確定；全球定位系統(tǒng)

1 引 言

衛(wèi)星定軌中需要求解三個(gè)微分方程，即：運(yùn)動微分方程、狀態(tài)轉(zhuǎn)移矩陣微分方程和敏感矩陣微分方程[1]。若通過解析法直接求得三個(gè)常微分方程解的顯式表達(dá)式，通常很困難甚至不可能。因此多采用數(shù)值解法，采用RK積分器、RKF積分器等單步積分器，也有Adams、Cowell等多步積分器。此外，還有配置積分法，該方法已經(jīng)成功應(yīng)用于國外Bernese軟件中[2]。文獻(xiàn)[3]分析了不同步長和積分器階數(shù)對積分的影響，得出采用積分步長為1小時(shí)、10階多項(xiàng)式的配置積分器進(jìn)行GPS精密星歷軌道擬合是合適的；文獻(xiàn)[4]比較了其他積分器和配置積分器在軌道擬合中的差別，指出了配置積分器自身的優(yōu)勢。此外，配置積分器還應(yīng)用于低軌衛(wèi)星定軌和重力場恢復(fù)、導(dǎo)航衛(wèi)星精密定軌中[5-9]。本文基于自主開發(fā)的衛(wèi)星精密定軌定位軟件對CODE提供的SP3精密星歷進(jìn)行軌道擬合，得到了較好的結(jié)果；同時(shí)，提出了將變步長配置積分器轉(zhuǎn)換為定步長的改進(jìn)算法。

2 配置積分器

2.1 基本原理

配置積分器的核心思想是：在一個(gè)不是很長的積分區(qū)間[ak,ak+1)，k=0,1,2…表示積分區(qū)間編號，將衛(wèi)星的位置、狀態(tài)轉(zhuǎn)移矩陣和敏感矩陣在邊界ak處用q次多項(xiàng)式泰勒展開；核心問題是：要已知衛(wèi)星位置速度、狀態(tài)轉(zhuǎn)移矩陣和敏感矩陣在ak點(diǎn)處的直到q次導(dǎo)數(shù)(q從0開始計(jì)數(shù))。此時(shí)，衛(wèi)星在區(qū)間[ak,ak+1)內(nèi)任意時(shí)刻的位置速度、狀態(tài)轉(zhuǎn)移矩陣和敏感矩陣可以由該泰勒多項(xiàng)式表示。

假設(shè)時(shí)間區(qū)間[ak,ak+1)內(nèi)時(shí)刻t處的狀態(tài)量z(t)在初始a時(shí)刻q次泰勒展開為：

(1)

(2)

其中，r表示衛(wèi)星的位置，φ和s為壓縮成列向量后的狀態(tài)轉(zhuǎn)移矩陣和敏感矩陣。為了形式需要，下文對其按照矩陣形式，根據(jù)文獻(xiàn)[9]以及文獻(xiàn)[1]中對三個(gè)微分方程的定義， 直接給出三個(gè)二

階微分方程：

(3)

令ak=a，ak+1=b，將該區(qū)間平均分為q-2段共q-1個(gè)節(jié)點(diǎn)，區(qū)間中的節(jié)點(diǎn)分別為ti，則有

(4)

對式(1)方程兩邊求2次導(dǎo)數(shù)，則得到:

(5)

動態(tài)生成的落腳點(diǎn)在于生成性。但它不是自然而然就能實(shí)現(xiàn)的，而是要通過一定的方法和實(shí)踐過程。過程哲學(xué)觀代表人物懷特海認(rèn)為，“世界事物分為事件的世界和永恒客體的世界。事件世界中的一切都處于變化的過程之中，過程就是機(jī)體各個(gè)因子之間有內(nèi)在聯(lián)系的、持續(xù)的創(chuàng)造活動，它表現(xiàn)為機(jī)體可以轉(zhuǎn)化為另一機(jī)體，因而整個(gè)世界就表現(xiàn)為一種活動的過程” [2]。但這個(gè)生成的過程和效果是由教師和學(xué)生的互動來決定的，而且要采取多元化的策略。如果互動形式過于單一，譬如缺乏層次性和深度、互動內(nèi)容過少等，則難以收到良好的效果。

(6)

將第0到第q-2共q-1個(gè)節(jié)點(diǎn)按照式(6)寫成矩陣形式為：

(7)

(8)

則可以得到初始時(shí)刻狀態(tài)量z的從2階到q階導(dǎo)數(shù)構(gòu)成的矩陣為：

(9)

式(9)中的任意第i列對應(yīng)初始時(shí)刻狀態(tài)量z的i階導(dǎo)數(shù)，加上式(3)中三個(gè)微分方程的初值作為狀態(tài)量z的零階導(dǎo)和一階導(dǎo)數(shù)，最終可以得到式(1)中的初始時(shí)刻狀態(tài)量z從0階直到q階導(dǎo)數(shù)；然后通過式(1)得到第0到q-2節(jié)點(diǎn)的狀態(tài)量，進(jìn)行迭代計(jì)算。通常需要迭代多次才能達(dá)到理想效果，迭代終止條件為tq-2節(jié)點(diǎn)處的衛(wèi)星位置前后兩次差值小于某一指定數(shù)值，該值通常取0.01m即可。

2.2 配置積分器的起步

2.3 算法改進(jìn)與實(shí)現(xiàn)

(10)

實(shí)際相當(dāng)于對式(1)中的系數(shù)采用以步長b-a作為單位。因此，采用歸一化時(shí)間作為單位和非歸一化時(shí)間作為單位時(shí)，式(1)中的泰勒系數(shù)之間的關(guān)系為：

(11)

3 數(shù)據(jù)實(shí)驗(yàn)與結(jié)果分析

為了驗(yàn)證配置積分器的積分效果，采用2010年年積日為4～13共10天的CODE精密星歷作為虛擬觀察量，動力學(xué)模型如表1所示。其中，積分器步長設(shè)置為1小時(shí)，多項(xiàng)式階數(shù)設(shè)定為10階，最大迭代次數(shù)設(shè)置為3次。選取GPS02和GPS03星擬合結(jié)果，在RTN三方向與CODE精密星歷比較，RTN三方向殘差如圖1和圖2所示，RMS三方向殘差如圖3和圖4所示。

表1 GPS軌道擬合動力學(xué)模型

地球引力場模型潮汐模型光壓模型相對論效應(yīng)攝動12×12EGM2008模型固體潮，海潮，極潮，固體極潮，海洋極潮，采用IERS2010推薦模型9參數(shù)太陽光壓模型IERS2010推薦模型

圖1 GPS02星RTN三方向殘差

圖2 GPS03星RTN三方向殘差

圖3 GPS02星RTN三方向RMS

圖4 GPS03星RTN三方向RMS

由以上結(jié)果可知：

(1)由圖1、圖2可以看出， GPS02星在R、T、N三方向殘差都在厘米級，除了T和N方向極個(gè)別歷元外，均不超過0.01m；GPS03星在R和N方向殘差不超過0.01m，在N方向殘差不超過0.02m；整體上看，兩顆星三個(gè)方向殘差均不超過0.02m。

(2)由圖3、圖4可以看出，兩顆星10天的擬合精度均比較穩(wěn)定，GPS02星T方向RMS最大為0.007m，R和N方向均在0.005m以內(nèi)；GPS03星N方向RMS稍差，在第5天和第13天接近0.01m，R和T方向較好，均不超過0.05m；整體上看，兩顆星三個(gè)方向RMS均不超過0.01m。

(3)綜合以上分析可以看出，配置積分器用于GPS衛(wèi)星精密星歷擬合，可以實(shí)現(xiàn)相對原精密星歷毫米級精度的連續(xù)軌道，并且該軌道可以用于事后GPS精密定位和星載GPS低軌衛(wèi)星精密定軌中。

4 結(jié) 論

與傳統(tǒng)的積分器相比，利用配置積分法對GPS軌道進(jìn)行數(shù)值積分時(shí)由于積分步長較大，因此所需的計(jì)算時(shí)間較少。然而，當(dāng)衛(wèi)星進(jìn)入地影或月影時(shí)，衛(wèi)星受力不連續(xù)，需要精確計(jì)算衛(wèi)星進(jìn)出地影的時(shí)間，最終導(dǎo)致上述積分區(qū)間不是等步長，產(chǎn)生變步長問題。采用本文提出的改進(jìn)方法，可以將變步長問題轉(zhuǎn)換為定步長問題，在簡化編程實(shí)現(xiàn)、提高計(jì)算效率的同時(shí)，保證了軌道擬合的精度。因此，本文提出的配置積分器改進(jìn)方法在衛(wèi)星定軌中具有一定的應(yīng)用價(jià)值。

[1]閆志闖, 張傳定, 杜蘭等. 基于面向?qū)ο蟮男l(wèi)星軌道確定優(yōu)化方法研究 [J]. 大地測量與地球動力學(xué), 2013, 33(2): 96-100.

[2]DACH R, HUGENTOBLER U, FRIDEZ P, et al. Bernese GPS Software Version 5.0 [M]. Bern: Astronomical Institute, University of Bern, 2007.

[3]李征航, 龔曉穎, 劉萬科. 配置法在衛(wèi)星軌道積分中的應(yīng)用 [J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2010, 35(3): 253-256.

[4]常志巧, 郝金明, 張成軍. 配置積分方法在GPS衛(wèi)星軌道確定中的應(yīng)用 [J]. 海洋測繪, 2006, 26(2): 31-33.

[5]BEUTLER G, J GGI A, MERVART L, et al. The celestial mechanics approach: application to data of the GRACE mission [J]. Journal of Geodesy, 2010, 84(11): 661-81.

[6]BEUTLER G, BROCKMANN E, HUGENTOBLER U, et al. Combining consecutive short arcs into long arcs for precise and efficient GPS Orbit Determination [J]. Journal of Geodesy, 1996, 70(5): 287-99.

[7]BEUTLER G, J GGI A, HUGENTOBLER U, et al. Efficient satellite orbit modelling using pseudo-stochastic parameters [J]. Journal of Geodesy, 2006, 80(7): 353-72.

[8]BEUTLER G, J GGI A, MERVART L, et al. The celestial mechanics approach: theoretical foundations [J]. Journal of Geodesy, 2010, 84(10): 605-24.

[9]BEUTLER G. Methods of Celestial Mechanics Volume 1: Physical, Mathematical, and Numerical Principles [M]. New York: Springer-Verlag, 2005.

GPS Satellite Precise Ephemeris Orbit Fitting Based on Collocation Integration

Yan Zhichuang1,2，Xu Junyi2, Li Yan3, Yu Chunping2

1. Information Engineering University, Zhengzhou 450001, China 2.Technical Division of Surveying and Mapping, Xi’an 710054, China 3.The First Surveying, Mapping & Navigation Base Dalian 116031, China

This paper deduces the principle of collocation integration in detail. In view of the variable step problem caused by earth or moon shadow, the paper proposes an optimum method based on matrix operation to improve the efficiency. The precise ephemeris of CODE is used as the virtual measurement to fit the orbit. The experiments results show that the RMS in the R, T and N direction is less than 0.01m compared to the GPS precise ephemeris.

collocation integration; orbit fitting; orbit determination; GPS

2015-03-30。

國家自然科學(xué)基金資助項(xiàng)目(41104047；41174026)。

閆志闖(1982—)，男,博士研究生，主要從事衛(wèi)星大地測量與導(dǎo)航方面的研究。

P237；TP751.1

A