陳強強，陳志平，李 芳

(1. 杭州電子科技大學(xué)機械工程學(xué)院，浙江 杭州 310018；2. 中國科學(xué)院國家天文臺，北京 100012)

國際GNSS服務(wù)中心(International GNSS Service, IGS)提供的精密星歷通常用于全球定位系統(tǒng)精密單點定位等數(shù)據(jù)處理中，但國際GNSS服務(wù)中心發(fā)布的精密星歷采樣間隔是15 min，而在全球定位系統(tǒng)精密定位中接收機的采樣率一般遠小于這個值，因此，需要對精密星歷進行高精度的插值。另一方面，由于國際GNSS服務(wù)中心精密星歷僅提供當天00：00：00～23：45：00時間段的星歷數(shù)據(jù)，若想僅用當天的數(shù)據(jù)得到23：45：00～24：00：00時間段的星歷數(shù)據(jù)，還需要進行外推以獲得任意時刻的衛(wèi)星坐標。

要使插值定位結(jié)果平滑、穩(wěn)定，要求插值多項式及其導(dǎo)數(shù)連續(xù)平滑，較為常見的插值方法有切比雪夫多項式插值法、牛頓多項式插值法和拉格朗日插值法等[1-3]。在目前的數(shù)據(jù)處理中，拉格朗日插值應(yīng)用較廣，文[4]指出拉格朗日插值法在本質(zhì)上與內(nèi)維爾插值法相同；文[5]將拉格朗日插值法改進，利用滑動式拉格朗日插值法獲得全球定位系統(tǒng)精密星歷；文[6]對比了拉格朗日和牛頓插值法的內(nèi)插效果，對插值結(jié)果進行了總結(jié)分析。上述研究中，多種插值方法在一定時段內(nèi)需要進行高低階組合才能達到最優(yōu)效果，且外推效果普遍不很理想。由于廣義延拓逼近法構(gòu)造的插值函數(shù)在所有同階插值函數(shù)中具有最小平方逼近誤差，為此，本文基于廣義延拓原理對全球定位系統(tǒng)精密星歷進行插值外推，利用求得的逼近函數(shù)解算任意時刻的全球定位系統(tǒng)衛(wèi)星的精密坐標。

1 廣義延拓內(nèi)插及外推模型

廣義延拓逼近法在分片邊界點上滿足插值條件，使分片間變化協(xié)調(diào)，并充分利用分片插值區(qū)域的周圍結(jié)點信息，實現(xiàn)分片區(qū)域內(nèi)部的最佳擬合，結(jié)合插值法和擬合法的優(yōu)點，在數(shù)據(jù)處理時實現(xiàn)了高精度的逼近[7-9]。

首先在整域內(nèi)進行剖分處理，以便在單元域內(nèi)尋找逼近函數(shù)及基函數(shù)，將定義域W剖分成n個互不重疊的子域Wi(i=1, 2, ...,n)：

(1)

單元域Wi有r個結(jié)點，對應(yīng)結(jié)點坐標為xe(e=1, 2, ...,r)，結(jié)點坐標xe(e=1, 2, ...,r)滿足定義域

xe-1,xe，

(2)

圖1 延拓域及相應(yīng)函數(shù)值

Fig.1 Extension domains and their corresponding function values

(3)

廣義延拓逼近法將單元域嵌套在延拓域中，吸收單元域外延拓域中的結(jié)點信息構(gòu)造單元域內(nèi)的擬合逼近函數(shù)，同時約束單元域邊界結(jié)點，實現(xiàn)相鄰單元逼近函數(shù)的協(xié)調(diào)和連續(xù)，從而達到單元域內(nèi)逼近函數(shù)的最佳擬合效果。利用廣義延拓構(gòu)造的逼近函數(shù)即可對目標變量進行插值處理。

廣義延拓插值模型不但可以用作內(nèi)插模型，也可以作為外推模型用于外延應(yīng)用的場合。對一組不斷增長的數(shù)據(jù)序列(x1,t1), (x2,t2), ..., (xi,ti), ..., (xn,tn)，i=1, 2, ...,n，已知tn以前的數(shù)據(jù)值，根據(jù)先驗數(shù)據(jù)的變化規(guī)律和趨勢，欲求tn+1時刻的xn+1值，數(shù)據(jù)外推示意圖如圖2。

圖2 數(shù)據(jù)外推示意圖
Fig.2 Data extrapolation intent

按照廣義延拓插值外推的設(shè)計理念，令tn為最新時刻，采用外推算法可得下一時刻的值xn+1。建立廣義延拓外推模型：

(4)

在上述模型中，由于插值點為最新采樣點的值，為了克服這一點數(shù)值突變帶來的誤差，考慮建立改進的廣義延拓模型：

(5)

即選擇最新的p個采樣點的平均值作為約束，改進模型的解法與原模型的解法基本一致，但數(shù)據(jù)的平穩(wěn)性更好，插值擬合曲線更平滑，外推精度更高。

2 全球定位系統(tǒng)精密星歷的內(nèi)插和外推方法

2.1 全球定位系統(tǒng)精密星歷的解算方法

由文[10]提供的解算算法，可將衛(wèi)星定位計算過程簡要描述如下：

(6)

其中，Xi，Yi，Zi為每個歷元待求的衛(wèi)星坐標。經(jīng)過坐標變換和攝動校正，最終得到t時刻該衛(wèi)星在WGS-84地心地固坐標系中的坐標計算矩陣：

(7)

由以上求解過程可知，利用星歷參數(shù)計算全球定位系統(tǒng)衛(wèi)星在某一時刻的空間位置必須提供準確的星歷數(shù)據(jù)，每定位一次需要的計算量也挺大。

2.2 全球定位系統(tǒng)精密星歷的插值可行性分析

圖3中的3條曲線分別是某顆衛(wèi)星的位置在WGS-84地心地固坐標系中的X,Y和Z分量隨時間的變化情況，歷元間隔為15 min，可見衛(wèi)星的軌道位置呈周期性變化。如圖4，衛(wèi)星位置的各分量在短時間內(nèi)變化平滑，幾乎呈線性變化。

圖3 長時間內(nèi)衛(wèi)星的空間位置
Fig.3 Space position of satellite over a long period of time

圖4 短時間內(nèi)衛(wèi)星的空間位置
Fig.4 Space position of satellite in short time

因此，衛(wèi)星的一段軌道可以用一個以時間域的插值多項式表達[11]，如果此逼近函數(shù)構(gòu)造得當，那么這種插值方法不會引入很大的衛(wèi)星位置誤差，可避免2.1節(jié)所述的復(fù)雜計算方法，使計算量大幅度減少。

2.3 基于廣義延拓的全球定位系統(tǒng)精密星歷內(nèi)插和外推方法

利用廣義延拓插值原理，對衛(wèi)星位置與速度在WGS-84地心地固坐標系中的X，Y，Z分量分別分段建立廣義延拓模型，在tk到tn這段時間內(nèi)，以衛(wèi)星位置在X方向上的分量為例，若在區(qū)間內(nèi)取兩端點作為先驗點進行約束，則廣義延拓插值模型為

(8)

其中，xk表示在插值區(qū)間左端點tk時刻衛(wèi)星位置在WGS-84地心地固坐標系中X方向上的分量；xn表示在插值區(qū)間右端點tn時刻衛(wèi)星位置在WGS-84地心地固坐標系中X方向上的分量；a1,a2,a3為待求系數(shù)；tk,tn是精密星歷更新的時間點。

對上述模型展開求解，即由

(9)

得到法方程CY=F，(9)式中，λ1,λ2為拉格朗日乘子，其中

(10)

接著，由(5)式構(gòu)造廣義延拓外推模型：

(11)

其中，變量含義見(8)式，p表示選用了最新p個衛(wèi)星位置的X分量，取平均值作為外推模型約束。對此模型的解法如下：

(12)

可得系數(shù)求解矩陣：

(13)

其中：

(14)

解得待求系數(shù)后，便可得到外推公式：

(15)

同理可得衛(wèi)星位置在其余各分量上的外推公式。

3 算例分析

現(xiàn)驗證廣義延拓插值外推模型的推算效果，采用國際GNSS服務(wù)中心提供的2017年6月6日(即GPS 1952周)采樣時間間隔為15 min的精密星歷數(shù)據(jù)，選用的衛(wèi)星為PRN15號。選取前4個采樣數(shù)據(jù)作為先驗采樣點，以拉格朗日插值法作為參照，并與精密星歷數(shù)據(jù)擬合的標準曲線對比，得到兩種精密星歷插值法的逼近效果如圖5。

圖5 兩種精密星歷插值法的插值外推效果對比
Fig.5 Comparison of interpolation and extrapolation effects of two precise ephemeris interpolation methods

由圖6可知，前4個比對點擬合的曲線為廣義延拓內(nèi)插效果，其插值的衛(wèi)星位置誤差小于5 cm，自第5個比對點開始為廣義延拓的外推效果，可見廣義延拓外推1小時左右的衛(wèi)星位置誤差小于10 cm，外推2小時的衛(wèi)星位置誤差小于20 cm，而在外推2小時以后，衛(wèi)星位置誤差開始迅速變大。

由于國際GNSS服務(wù)中心提供的精密星歷誤差精度為5 cm，因此，利用廣義延拓法對全球定位系統(tǒng)精密星歷的內(nèi)插能夠滿足精度要求，而在利用廣義延拓外推時，在30 min內(nèi)能維持精度，在較長時間內(nèi)也不會引入很大的位置誤差。

4 結(jié) 論

圖6 衛(wèi)星位置隨時間變化圖
Fig.6 Chart of satellite position versus time

采用廣義延拓法對全球定位系統(tǒng)精密星歷進行衛(wèi)星位置內(nèi)插時效果較好，可以將衛(wèi)星的一段軌道用一個以時間為坐標的模型表達，使得計算衛(wèi)星位置時在不損失精度的情況下，不再進行復(fù)雜冗長的步驟，極大地減少了計算量。此外，通過實例分析，廣義延拓外推法能夠保證外推1小時的精度滿足要求，廣義延拓內(nèi)插法可獲得全球定位系統(tǒng)任一時刻的精密星歷數(shù)據(jù)。

[1] 李征航, 黃勁松. GPS測量與數(shù)據(jù)處理[M]. 武漢: 武漢大學(xué)出版社, 2005.

[2] 孫騰科. 基于拉格朗日與切比雪夫的精密星歷插值研究[J]. 測繪與空間地理信息, 2014, 37(2): 33-37.

Sun Tengke. The research of GPS precise interpolation methods based on Lagrange and Chebyshev[J]. Geomatics & Spatial Information Technology, 2014, 37(2): 33-37.

[3] 汪威, 陳明劍, 閆建巧, 等. 北斗三類衛(wèi)星精密星歷內(nèi)插方法比較分析[J]. 全球定位系統(tǒng), 2016, 41(2): 60-65.

Wang Wei, Chen Mingjian, Yan Jianqiao, et al. Three kinds of compass satellite precise ephemeris interpolation method analysis comparative[J]. GNSS World of China, 2016, 41(2): 60-65.

[4] 王超, 郭際明, 周命端, 等. 高精度GPS數(shù)據(jù)處理中GAMIT批處理方法與實現(xiàn)[J]. 測繪信息與工程, 2012, 37(2): 10-12.

Wang Chao, Guo Jiming, Zhou mingduan, et al. A method of GAMIT batch processing and its implementation in high precise GPS data processing[J]. Journal of Geomatics, 2012, 37(2): 10-12.

[5] 雷雨, 趙丹寧, 高玉萍. 基于滑動式Lagrange插值方法的GPS精密星歷內(nèi)插分析[J]. 測繪工程, 2013, 22(2): 34-36.

Lei Yu, Zhao Danning, Gao Yuping. Analysis of interpolation for GPS precise ephemeris using sleek Lagrange interpolation[J]. Engineering of Surveying and Mapping, 2013, 22(2): 34-36.

[6] 柳笛, 逢淑濤, 董緒榮. IGS精密星歷文件的讀取及內(nèi)插方法研究[J]. 全球定位系統(tǒng), 2011, 36(5): 46-48+64.

Liu Di, Pang Shutao, Dong Xurong. Read of precise ephemeris file and research on methods of interpolation[J]. GNSS World of China, 2011, 36(5): 46-48+64.

[7] 施滸立, 顏毅華, 徐國華, 等. 工程科學(xué)中的廣義延拓逼近發(fā)[M]. 北京: 科學(xué)出版社, 2005.

[8] 魏彥飛, 耿建平, 施滸立, 等. 一種新的數(shù)據(jù)融合方法—廣義融合法[J]. 天文研究與技術(shù), 2016, 13(3): 318-325.

Wei Yanfei, Geng Jianping, Shi Huli, et al. A new method of data fusion—the generalized fusion method[J]. Astronomical Research & Technology, 2016, 13(3): 318-325.

[9] 耿建平, 衣偉, 劉成, 等. 基于廣義延拓外推的單頻周跳檢測與修復(fù)方法[J]. 天文研究與技術(shù), 2015, 12(2): 174-182.

Geng Jianping, Yi Wei, Liu Cheng, et al. A new method for detection and correction of single-frequency cycle slips using data extrapolation based on the method of generalized extended interpolation[J]. Astronomical Research & Technology, 2015, 12(2): 174-182.

[10]Jiang R B, Liu X Q. Using fourier series to fit the GPS precise ephemeris[C]// 2010 International Conference on Computer Application and System Modeling (ICCASM 2010). 2010: 96-99.

[11]Liu W P, Hao J M. A new interpolation method based on satellite physical character in using IGS precise ephemeris[J].Geodesy and Geodynamics, 2014, 5(3): 29-33.