劉隆迪，杜蘭，劉澤軍，張中凱，周佩元，黃俊迦

(信息工程大學(xué) 地理空間信息學(xué)院,鄭州 450001)

0 引言

依據(jù)機(jī)構(gòu)間空間碎片協(xié)調(diào)委員會(huì)(IADC)的劃定原則，地球同步軌道(GEO)衛(wèi)星帶定義為衛(wèi)星緯度約在-15°～15°、軌道高度介于35 586～35 986 km 的空間范圍[1].為挖掘GEO 有限的軌位資源，采用多星共位策略[2].因此，一方面GEO 帶內(nèi)的空間目標(biāo)密集(截至2022年4月，僅北美防空司令部公開發(fā)布的數(shù)量多達(dá)1 026 顆)，戰(zhàn)略價(jià)值高且運(yùn)行穩(wěn)定.另一方面，GEO 帶內(nèi)的廢棄飛行器通常不能進(jìn)入大氣層自行墜毀，越來越密集的GEO 帶合作/非合作目標(biāo)及碎片等給空間目標(biāo)庫(kù)的維護(hù)與管理帶來挑戰(zhàn).

國(guó)際上，空間編目目標(biāo)采用初軌+分析解的慣性系下軌道外推的通用方法[3-4].初軌采用雙行軌道根數(shù)(TLE)，能夠表征一定精度意義下的軌道初值和少量攝動(dòng)參數(shù)；外推軌道則是基于簡(jiǎn)化動(dòng)力學(xué)模型下的一套復(fù)雜的分析解公式，如簡(jiǎn)化常規(guī)/深空攝動(dòng)4 模型(SGP4/SDP4)等.這種TLE-SDP4 模型具有通用和高效的優(yōu)點(diǎn)，但是對(duì)于GEO 帶目標(biāo)，大多數(shù)地面用戶更關(guān)注其對(duì)地的運(yùn)動(dòng)特性.此外，為滿足公開程序的調(diào)用條件或精度需求等，用戶須頻繁地進(jìn)行目標(biāo)平位置與真位置的轉(zhuǎn)換計(jì)算，且軌道精度受限(外推1 天在千米級(jí))[5]，公式繁瑣，不便于用戶理解.

針對(duì)上述問題，采用了數(shù)值積分等方式獲得高精度外推軌道、再對(duì)外推軌道進(jìn)行數(shù)學(xué)函數(shù)擬合表達(dá)的兩步法.文獻(xiàn)[6]和文獻(xiàn)[7]根據(jù)切比雪夫多項(xiàng)式、傅里葉級(jí)數(shù)和正弦函數(shù)提出分點(diǎn)根數(shù)擬合法，用戶使用方便，但其表征5 天軌道的參數(shù)多達(dá)343 個(gè)，且軌道的擬合精度仍為百米級(jí).文獻(xiàn)[8]和文獻(xiàn)[9]均利用切比雪夫多項(xiàng)式擬合GEO 軌道，擬合精度達(dá)到厘米級(jí)，但擬合時(shí)長(zhǎng)僅為1 h和2 h，軌道更新過于頻繁.此外，上述方法均沒有考慮軌道的對(duì)地特性.

本文直接針對(duì)12～24 h 對(duì)地運(yùn)動(dòng)，對(duì)GEO 帶目標(biāo)軌跡進(jìn)行參數(shù)優(yōu)化表征，設(shè)計(jì)了一種結(jié)合經(jīng)典開普勒根數(shù)和多項(xiàng)式-傅里葉級(jí)數(shù)的混合函數(shù)型參數(shù)星歷模型.長(zhǎng)期擬合實(shí)驗(yàn)表明，1 d 擬合精度優(yōu)于20 m，12 h 擬合精度優(yōu)于0.15 m，適用于目前的空間編目軌道精度水平和更新周期.此外，若要求發(fā)布的目標(biāo)軌道有效期較長(zhǎng)，可采取同時(shí)發(fā)布M(M＞1)天的2M或M組星歷的靈活方法，且保證每組的擬合精度一致.

1 定點(diǎn)附近的運(yùn)動(dòng)特征

1.1 定點(diǎn)R、T、N 坐標(biāo)系

在地心地固坐標(biāo)系(ECEF)下，理想GEO 軌道的星下點(diǎn)為赤道上的固定點(diǎn)，唯一自由度是定點(diǎn)經(jīng)度.不失一般性，GEO 帶目標(biāo)的星下點(diǎn)軌跡，表現(xiàn)為關(guān)于定點(diǎn)經(jīng)度處的小8 字，8 字的上下幅度范圍由軌道傾角決定.

為分析GEO 的空間小8 字的短期局域位置變化，新建一個(gè)定點(diǎn)徑向(R)、沿跡方向(T)、外法向(N)坐標(biāo)系.原點(diǎn)位于地心，令參考時(shí)刻的GEO 位置為“定點(diǎn)”，則坐標(biāo)系的三軸分別指向R、T、N.特別地，若GEO 的傾角為0，則定點(diǎn)落在赤道上，且T和N分別為沿赤道方向和垂直赤道面方向.

將其他時(shí)刻的ECEF 位置轉(zhuǎn)換至定點(diǎn)R、T、N坐標(biāo)系，即可分析GEO 在“定點(diǎn)”附近的運(yùn)動(dòng)特征.轉(zhuǎn)換過程如下：

1)計(jì)算參考時(shí)刻t0的開普勒根數(shù)

為簡(jiǎn)便起見，不必采用地心慣性坐標(biāo)系(ECI)，而是定義一個(gè)準(zhǔn)ECI，即令ECI繞Z軸旋轉(zhuǎn)，使其X軸與t0時(shí)刻ECEF的X軸指向保持一致.此時(shí)，ECEF和準(zhǔn)ECI 下的位置向量保持一致，速度轉(zhuǎn)換至準(zhǔn)ECI下要考慮地球自轉(zhuǎn)角速度.利用準(zhǔn)ECI 下的位置和速度向量計(jì)算相應(yīng)的開普勒根數(shù).

2)將t時(shí)刻準(zhǔn)ECI 下的位置向量旋轉(zhuǎn)至t0時(shí)刻的定點(diǎn)R、T、N坐標(biāo)系

式中：Ri()表 示繞i軸(i=x,y,z)的旋轉(zhuǎn)矩陣，其中的旋轉(zhuǎn)角來自t0的開普勒根數(shù)，即u0(u0=ω0+f0)、i0和Ω0分別為t0時(shí)刻的衛(wèi)星緯度幅角、軌道傾角和升交點(diǎn)赤經(jīng).

需要說明的是，這里的旋轉(zhuǎn)參數(shù)采用的是開普勒根數(shù)，也可以采用適用于GEO 的無奇點(diǎn)根數(shù)，不會(huì)影響坐標(biāo)旋轉(zhuǎn)后的指向精度.此外，這里可以統(tǒng)一進(jìn)行位置時(shí)序的整批轉(zhuǎn)換.

3)將t時(shí)刻的位置向量從準(zhǔn)ECI(小寫表示)轉(zhuǎn)換至ECEF(大寫表示)

若僅考慮地球均速自轉(zhuǎn)運(yùn)動(dòng)，有

式中：S=ωe(t-t0)，是t與t0時(shí)刻的格林尼 治時(shí)角之差；ωe是地球自轉(zhuǎn)平均角速率.這里慣性系向非慣性系轉(zhuǎn)換，只能逐點(diǎn)計(jì)算.

4)任意t時(shí)刻的位置向量從ECEF 轉(zhuǎn)換至定點(diǎn)R、T、N坐標(biāo)系

將式(2)代入式(1)，有

特別地，若軌道傾角無限趨于0，式(3)可進(jìn)一步簡(jiǎn)化為

此時(shí)，ECEF 轉(zhuǎn)換至定點(diǎn)R、T、N坐標(biāo)系，僅需要沿赤道旋轉(zhuǎn)至t0時(shí)刻的定點(diǎn)經(jīng)度 λe.

1.2 R、T、N 位置分量的運(yùn)動(dòng)分析

GEO 的主要攝動(dòng)力包括地球非球形引力攝動(dòng)、日月引力攝動(dòng)和光壓等.在自然攝動(dòng)作用下，軌道根數(shù)和位置分量表現(xiàn)為長(zhǎng)期、長(zhǎng)周期和短周期攝動(dòng)運(yùn)動(dòng)特性的疊加[10-11].

在定點(diǎn)R、T、N坐標(biāo)系下，可以方便地觀察GEO目標(biāo)的對(duì)地局域運(yùn)動(dòng)特性.以北斗衛(wèi)星導(dǎo)航系統(tǒng)(BDS)的BDS-G8 衛(wèi)星(軌道傾角約為1°)為例，圖1 給出了27 d 的三分量位置變化，其中R方向扣除了靜止軌道半長(zhǎng)軸的參考值aref.

圖1 BDS-G8 衛(wèi)星在定點(diǎn)R、T、N 坐標(biāo)系下的位置分量時(shí)序

由圖1 可知，在各種攝動(dòng)力作用下，位置三分量均呈現(xiàn)出趨勢(shì)項(xiàng)和軌道周期的諧振疊加變化.R和T方向有明顯的趨勢(shì)項(xiàng)，T方向長(zhǎng)期累積最為顯著；三分量均有衛(wèi)星軌道周期的1 ×T運(yùn)動(dòng)，N方向上還有約半月的長(zhǎng)周期項(xiàng)，其他短周期項(xiàng)不顯著.

若僅考慮1 天的尺度，可以將長(zhǎng)周期攝動(dòng)視為長(zhǎng)期攝動(dòng).為此，首先采用二階多項(xiàng)式進(jìn)行趨勢(shì)項(xiàng)擬合，并對(duì)去除趨勢(shì)項(xiàng)的1 天數(shù)據(jù)進(jìn)行時(shí)頻變換，提取軌道周期的諧振項(xiàng)振幅.

表1 中給出了前8 階諧振項(xiàng)的幅度(以km為單位)及累計(jì)占比.可以看出，R、T、N三方向振幅能量中，衛(wèi)星軌道運(yùn)動(dòng)周期(1 ×T)項(xiàng)最大，其次依序?yàn)橹芷谥C振項(xiàng)，且前三項(xiàng)諧振的幅度能量之和分別占各自分量總量的98.2%、97.7%和76.2%.雖然其中N方向能量分布較為分散，但是其短周期攝動(dòng)的總體幅度遠(yuǎn)小于其他兩個(gè)方向.

表1 8 階傅里葉級(jí)數(shù)的幅度值與累計(jì)占比

2 星歷參數(shù)設(shè)計(jì)和用戶算法

2.1 星歷參數(shù)集

定點(diǎn)R、T、N坐標(biāo)系下，三軸方向能表征出GEO 目標(biāo)的局域運(yùn)動(dòng)特性，其中R方向表征GEO 目標(biāo)徑向運(yùn)動(dòng)狀態(tài)，T方向表征GEO 目標(biāo)關(guān)于定點(diǎn)的東西運(yùn)動(dòng)，N方向表征GEO 目標(biāo)受第三體引力等攝動(dòng)力的影響所引起的周期振蕩且振幅不受GEO 軌道傾角影響.因此，定點(diǎn)R、T、N坐標(biāo)系下的位置三分量采取統(tǒng)一的標(biāo)準(zhǔn)進(jìn)行擬合表征.

在1 天尺度下，GEO 目標(biāo)在在定點(diǎn)R、T、N坐標(biāo)系下的位置三分量可以表征為趨勢(shì)項(xiàng)和短周期項(xiàng)的疊加.此外，由圖1 可知，在1 天尺度下，該三分量在趨勢(shì)項(xiàng)上更多呈現(xiàn)出二階多項(xiàng)式的特征.1.2 節(jié)中的8 階傅里葉級(jí)擬合結(jié)果表明對(duì)濾去趨勢(shì)項(xiàng)的位置三分量可采用3 階傅里葉級(jí)數(shù)擬合.為此，GEO 目標(biāo)在定點(diǎn)R、T、N坐標(biāo)系下的位置三分量可以統(tǒng)一優(yōu)化表達(dá)為二次多項(xiàng)式和3 階軌道周期諧振項(xiàng)的主要運(yùn)動(dòng)疊加.

為此，結(jié)合開普勒根數(shù)和多項(xiàng)式-傅里葉級(jí)數(shù)設(shè)計(jì)了ECEF 下混合函數(shù)型星歷.星歷參數(shù)集分為兩大類，一類是定點(diǎn)R、T、N坐標(biāo)系轉(zhuǎn)換至ECEF 的軌道根數(shù)子集，一類是定點(diǎn)R、T、N坐標(biāo)系下的數(shù)學(xué)函數(shù)系數(shù)集.共計(jì)31 個(gè)參數(shù)的星歷參數(shù)集為：

式中：第一組的4 個(gè)參數(shù)分別為星歷參考時(shí)刻toe及其在準(zhǔn)ECI 系下的緯度幅角、軌道傾角和升交點(diǎn)赤經(jīng)；后三組分別為定點(diǎn)R、T、N坐標(biāo)系下每個(gè)方向的9 個(gè)系數(shù)，即二次多項(xiàng)式系數(shù)k=(k0,k1,k2)T以及一至三階諧振短周期項(xiàng)的正弦系數(shù)a=(a1,a2,a3)T和余弦系數(shù)b=(b1,b2,b3)T.注意后面三組系數(shù)需要擬合計(jì)算給出.

已知ECEF 下GEO 目標(biāo)1 天弧長(zhǎng)軌道位置，其星歷參考時(shí)刻為toe，擬合計(jì)算混合函數(shù)型星歷步驟如下：

1)通過內(nèi)插計(jì)算得到星歷參考時(shí)刻toe的速度矢量；

2)將ECEF 下的位置矢量和步驟1)中得到的速度矢量轉(zhuǎn)換到準(zhǔn)ECI 下；

3)計(jì)算星歷參考時(shí)刻toe下，GEO 衛(wèi)星在準(zhǔn)ECI 系下的緯度幅角u0、軌道傾角i0和升交點(diǎn)赤經(jīng) Ω0；

4)根據(jù)式(1)將步驟2)中得到的準(zhǔn)ECI 下的位置矢量轉(zhuǎn)換至toe時(shí)刻的定點(diǎn)R、T、N坐標(biāo)系；

5)對(duì)步驟4)中得到的定點(diǎn)R、T、N坐標(biāo)系下的三個(gè)分量分別進(jìn)行最小二乘擬合.這里每個(gè)分量擬合系數(shù)包括二次多項(xiàng)式系數(shù)、正弦系數(shù)和余弦系數(shù)共9 個(gè)系數(shù).

2.2 GEO 位置的用戶算法

已知參考時(shí)刻toe的星歷參數(shù)，計(jì)算星歷有效時(shí)段內(nèi)任意時(shí)刻t的位置算法如下：

1)計(jì)算定點(diǎn)R、T、N坐標(biāo)下的位置向量

式中，aref為地球同步軌道半長(zhǎng)軸的參考值，取為42 165.76 km.

此外，為顧及模型的通用性，利用了GEO 的靜地特性，將短周期項(xiàng)的本征角頻率統(tǒng)一取為地球平均自轉(zhuǎn)角速率 ωe.

2)轉(zhuǎn)換至ECEF 下的位置向量

與式(3)相一致，有

其中，Ωe和ue均包含了t與toe時(shí)刻間的地球自轉(zhuǎn)效應(yīng)，即

3 星歷擬合實(shí)驗(yàn)

顧及定點(diǎn)經(jīng)度和軌道傾角分布，表2 中選用了兩類GEO 目標(biāo)：

表2 GEO 目標(biāo)信息

1)BDS-G8 衛(wèi)星(傾角1°)：從武漢大學(xué)IGS 數(shù)據(jù)中心下載了1年的精密星歷，數(shù)據(jù)間隔1 min；其中，有約每月一次的東西機(jī)動(dòng)和每半年一次的南北機(jī)動(dòng)，機(jī)動(dòng)期間無精密星歷數(shù)據(jù).

2)兩個(gè)退役GEO 目標(biāo)：美國(guó)國(guó)防衛(wèi)星通信系統(tǒng)衛(wèi)星DSCS 2-2(傾角6.6°)和俄羅斯彩虹衛(wèi)星RADUGA 1-8(傾角8.2°)，初軌采用TLE，外推1 a 軌道的攝動(dòng)力包括10×10 的地球非球形引力攝動(dòng)、日月引力攝動(dòng)和太陽(yáng)光壓，退役后不再施加軌道機(jī)動(dòng)，數(shù)據(jù)間隔1 min.外推軌道位置已轉(zhuǎn)換至ECEF.

將1 a 的實(shí)驗(yàn)數(shù)據(jù)分別按1 天和12 h 分組，其中1 天分組的組間有12 h 重疊，兩個(gè)分組組數(shù)均為726 組.采用最小二乘擬合解算星歷模型中的混合函數(shù)系數(shù)項(xiàng)，參考?xì)v元均取各分組的中間時(shí)刻.

3.1 1 天星歷擬合精度

為了比對(duì)混合函數(shù)型星歷，這里也利用10 階切比雪夫多項(xiàng)式分別對(duì)ECEF 下的x、y、z位置分量時(shí)序進(jìn)行建模，兩者的擬合參數(shù)個(gè)數(shù)一致.參數(shù)個(gè)數(shù)相同的條件下，混合函數(shù)的擬合精度優(yōu)于切比雪夫多項(xiàng)式的參數(shù)擬合.圖2 分別給出了兩者對(duì)BDS-G8 任意1 天數(shù)據(jù)的擬合殘差分布.從最大振幅來看，混合函數(shù)約為切比雪夫多項(xiàng)式的一半，龍格現(xiàn)象也更小.從R、T、N三方向的擬合RMS 來看，切比雪夫多項(xiàng)式為4.8 m、15.2 m和7.8 m，混合函數(shù)為4.9 m、3.0 m和5.6 m，顯然后者在T方向表現(xiàn)最佳，其次是N方向.此外，混合函數(shù)顧及到前三階的周期諧振項(xiàng)，因此周期震蕩成為擬合殘差的主周期項(xiàng).

圖2 兩種星歷參數(shù)模型擬合BDS-G8 1 d 軌道的殘差

3.2 星歷擬合時(shí)序穩(wěn)定性

為驗(yàn)證混合函數(shù)型星歷模型的長(zhǎng)期穩(wěn)定性，圖3給出了三個(gè)GEO 實(shí)驗(yàn)?zāi)繕?biāo)的1 a 星歷參數(shù)擬合的殘差RMS 時(shí)序.注意其中的兩段柱狀陰影分別對(duì)應(yīng)春分和秋分附近的兩次地影季.可以看出：1)由于軌道面定向參數(shù)略有不同[12]，GEO 目標(biāo)的地影季并不完全一致；2)地影季基本上不會(huì)影響到1 天星歷的擬合精度，這是因?yàn)椋刻爝M(jìn)出地影的最大時(shí)長(zhǎng)不超過72 min[13]，造成的光壓間斷與1 天擬合時(shí)長(zhǎng)相比可以忽略.

圖3 三個(gè)GEO 目標(biāo)1年的1 天軌道擬合殘差RMS 時(shí)序

表3為一年內(nèi)星歷擬合殘差RMS 的統(tǒng)計(jì)結(jié)果.可以看出，混合函數(shù)型星歷參數(shù)模型對(duì)三個(gè)GEO 實(shí)驗(yàn)?zāi)繕?biāo)的1 天軌道表征能力基本相當(dāng)，不受GEO 定點(diǎn)經(jīng)度與軌道傾角的影響.R、T、N位置分量的平均RMS 均優(yōu)于10 m，最大RMS 不超過20 m.因此，混合函數(shù)型星歷模型具有長(zhǎng)期的擬合穩(wěn)定性，可用于GEO 帶多目標(biāo)的軌道表征.

表3 混合函數(shù)型參數(shù)擬合1 天軌道殘差RMS 統(tǒng)計(jì) m

3.3 不同弧長(zhǎng)的星歷擬合精度

GEO 帶目標(biāo)除了日常軌道監(jiān)測(cè)，還可能有更高軌道精度的任務(wù)需求.例如空間態(tài)勢(shì)感知(SSA)，它對(duì)重點(diǎn)監(jiān)測(cè)目標(biāo)或碰撞預(yù)警分析要求徑向優(yōu)于50 m[14].

不改變星歷模型的條件下，可以通過縮短擬合時(shí)段、提高星歷更新頻率的方式給予精度改善.仍采用BDS-G8 衛(wèi)星的數(shù)據(jù)，圖4 給出了混合函數(shù)對(duì)同一參考?xì)v元12～28 h 弧段的擬合精度變化.可以發(fā)現(xiàn)，擬合時(shí)段在12～20 h，表征精度能控制在1 m 內(nèi)；當(dāng)擬合時(shí)長(zhǎng)超過24 h，T和N方向表征精度迅速下降.因此，混合函數(shù)型星歷模型擬合精度與擬合弧長(zhǎng)呈負(fù)相關(guān)，擬合弧長(zhǎng)越長(zhǎng)，擬合精度越低.傅里葉級(jí)數(shù)中的角頻率為地球自轉(zhuǎn)角頻率，為此不推薦混合函數(shù)型星歷模型用于弧長(zhǎng)超過24 h 的軌道擬合.

圖4 混合函數(shù)對(duì)12～28 h 弧段的擬合精度

以12 h 擬合時(shí)段為例，圖5 給出了三個(gè)GEO 實(shí)驗(yàn)?zāi)繕?biāo)的1 a 星歷參數(shù)擬合的殘差RMS 時(shí)序.與圖3的24 h 擬合相比：1)12 h 的擬合精度提高了近兩個(gè)數(shù)量級(jí)，位置三分量的RMS 均優(yōu)于0.15 m，且N分量最優(yōu)；2)12 h 擬合受地影季的影響，體現(xiàn)在R分量上，但對(duì)總體影響較小.因此，混合函數(shù)型星歷模型可根據(jù)SSA 的精度需求靈活調(diào)整12～24 h 的擬合時(shí)段.

圖5 三個(gè)GEO 目標(biāo)1 a 的12 h 軌道的擬合殘差RMS 時(shí)序

4 總結(jié)

空間高軌GEO 編目目標(biāo)采用通用的TLE-SDP4模型計(jì)算軌道，這種分析解軌道表達(dá)精度較低，使用復(fù)雜，不易理解.純數(shù)學(xué)函數(shù)如切比雪夫多項(xiàng)式和傅里葉級(jí)數(shù)等用于表征軌道，雖然計(jì)算簡(jiǎn)單，但表征時(shí)長(zhǎng)與表征精度的矛盾突出，且參數(shù)無物理意義.

本文針對(duì)GEO 對(duì)地運(yùn)動(dòng)特性，提出一種結(jié)合經(jīng)典開普勒根數(shù)和多項(xiàng)式-傅里葉級(jí)數(shù)的混合函數(shù)的星歷模型.該模型可靈活適應(yīng)用戶對(duì)軌道擬合精度的需求，多目標(biāo)長(zhǎng)期擬合實(shí)驗(yàn)表明，1 天擬合RMS 均值優(yōu)于20 m，12 h 軌道擬合精度優(yōu)于0.15 m.

星歷模型也能夠適應(yīng)用戶對(duì)軌道外推時(shí)長(zhǎng)的需求.例如同時(shí)解算M天的2M或M組星歷，適用于多天軌道發(fā)布.