任家棟， 張大力， 曾慶雙

(哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱 150001)

在目前的航天器相對(duì)運(yùn)動(dòng)框架中，動(dòng)力學(xué)模型呈多樣化趨勢(shì)，具體取決于兩星的距離、運(yùn)行軌道的扁率、運(yùn)動(dòng)遞推的運(yùn)算復(fù)雜度等因素，難以用單一標(biāo)準(zhǔn)來(lái)明確某個(gè)動(dòng)力學(xué)模型更加適合或更加有效[1]。相對(duì)運(yùn)動(dòng)描述方法分為相對(duì)軌道要素描述法和直角坐標(biāo)描述法2類[2]。本文主要針對(duì)相對(duì)導(dǎo)航制導(dǎo)設(shè)計(jì)應(yīng)用場(chǎng)景，如二階段濾波、非線性預(yù)測(cè)估計(jì)等濾波方法需要明確的干擾誤差特性，基于此，對(duì)直角坐標(biāo)系描述方法進(jìn)行分析。

多年來(lái)，研究人員在CW方程[3]的基礎(chǔ)上，從近圓和近距離等理想化假設(shè)條件下出發(fā)，向保留高階項(xiàng)、考慮軌道攝動(dòng)、橢圓軌道等幾方面進(jìn)行了擴(kuò)展研究。文獻(xiàn)[4-5]分別研究了圓軌道假設(shè)下二階、三階項(xiàng)對(duì)相對(duì)運(yùn)動(dòng)的影響。Lawden[6]以真近點(diǎn)角為參變量，給出了橢圓軌道下的相對(duì)動(dòng)力學(xué)解析描述，稱為L(zhǎng)awden方程，已經(jīng)廣泛應(yīng)用于橢圓軌道的相對(duì)運(yùn)動(dòng)控制領(lǐng)域。文獻(xiàn)[7-8]采用逐次逼近的QV(Quadratic-Volterra solution)方法研究了橢圓軌道下考慮高階項(xiàng)的分析模型。

軌道環(huán)境攝動(dòng)，如地球非球形引力攝動(dòng)、大氣阻力攝動(dòng)、太陽(yáng)光壓攝動(dòng)等均會(huì)導(dǎo)致期望軌跡產(chǎn)生偏差，是影響編隊(duì)飛行的關(guān)鍵因素?；贑W方程，文獻(xiàn)[9-10]提出了一種地球非球形引力攝動(dòng)的補(bǔ)償模型，提升了相對(duì)運(yùn)動(dòng)精度。文獻(xiàn)[11]將變化幅度較小的地心距及軌道傾角平均化，推導(dǎo)了考慮J2項(xiàng)攝動(dòng)干擾的近地軌道相對(duì)運(yùn)動(dòng)解析方程。文獻(xiàn)[12]采用分析法研究了三體引力對(duì)衛(wèi)星相對(duì)運(yùn)動(dòng)的影響，三體引力主要對(duì)對(duì)軌道面內(nèi)運(yùn)動(dòng)產(chǎn)生長(zhǎng)周期作用。太陽(yáng)光壓是高軌衛(wèi)星的主要攝動(dòng)元素之一，對(duì)于高面質(zhì)比(high area-to-mass ratio,HAMR)目標(biāo)其影響甚至占主導(dǎo)地位。太陽(yáng)光壓攝動(dòng)一般考慮為非保守力，但在全軌道考慮時(shí)會(huì)呈現(xiàn)一些近似保守力的影響特性，如平均軌道半長(zhǎng)軸不變、軌道偏心率或傾角以年或年的倍數(shù)周期進(jìn)行波動(dòng)[13-14]。

目前對(duì)CW方程的各誤差要素的影響分析成果較多，但缺少地球靜止軌道小傾角、小偏心率等特性下針對(duì)相對(duì)導(dǎo)航制導(dǎo)應(yīng)用場(chǎng)景的分析，本文結(jié)合地球靜止軌道的特點(diǎn)，推導(dǎo)CW方程的建模誤差，并對(duì)其主要誤差項(xiàng)進(jìn)行分析，以得到誤差的規(guī)律和量化特性。

1 高軌環(huán)境分析與建模

衛(wèi)星在軌道上始終受到各種攝動(dòng)力的作用，主要包括：地球非球形攝動(dòng)，高層大氣的氣動(dòng)力攝動(dòng)，太陽(yáng)、月亮的引力攝動(dòng)，太陽(yáng)光輻射壓力攝動(dòng)等。在攝動(dòng)力的作用下，衛(wèi)星軌道偏心率e、周期、升交點(diǎn)赤經(jīng)Ω和軌道傾角i等不斷演化，呈現(xiàn)周期或非周期特性。雖然攝動(dòng)力量級(jí)通常較小，但其長(zhǎng)期影響不容忽視，需對(duì)各攝動(dòng)的特性進(jìn)行詳細(xì)數(shù)值分析。

由于大氣阻力隨軌道高度的增加而急劇減小，地球靜止軌道攝動(dòng)主要考慮地球非球形引力攝動(dòng)、日月引力攝動(dòng)(也稱三體引力攝動(dòng))、太陽(yáng)光壓攝動(dòng)，其量級(jí)特性見表1[15]。

表1 地球靜止軌道各攝動(dòng)力與理想地球引力比值

下面以柯氏定理為基礎(chǔ)對(duì)相對(duì)運(yùn)動(dòng)模型進(jìn)行推導(dǎo)。慣性系下的矢量在轉(zhuǎn)動(dòng)坐標(biāo)系下求導(dǎo)為:

(1)

目標(biāo)星(target)和追蹤星(chaser)的絕對(duì)動(dòng)力學(xué)方程分別為：

(2)

式中：rt、rc分別為目標(biāo)星和追蹤星在慣性坐標(biāo)系下的位置矢量；at、ac分別為目標(biāo)星和追蹤星的環(huán)境攝動(dòng)加速度。

圖1 地心軌道坐標(biāo)系Fig.1 Geocentric coordinate system

相對(duì)運(yùn)動(dòng)通常在衛(wèi)星軌道坐標(biāo)系下研究，軌道坐標(biāo)系采用VVLH坐標(biāo)系，定義為OoXoYoZo。

軌道坐標(biāo)系OoXoYoZo不是慣性系，因此在軌道坐標(biāo)系下推導(dǎo)相對(duì)運(yùn)動(dòng)方程需考慮OoXoYoZo的轉(zhuǎn)動(dòng)信息。為了便于研究其轉(zhuǎn)動(dòng)特性，建立地心軌道系OXYZ，原點(diǎn)在地心，相當(dāng)于軌道坐標(biāo)系平移到地球球心。因此地心軌道系OXYZ和星體軌道系OoXoYoZo有同樣的轉(zhuǎn)動(dòng)特性。

(3)

顯然，ωx=0，即軌道坐標(biāo)系不存在X軸方向的轉(zhuǎn)動(dòng)分量。

根據(jù)柯氏定理求二階導(dǎo)：

(4)

式中(Fx,Fy,Fz)為軌道攝動(dòng)加速度。解式(4)，并結(jié)合兩體運(yùn)動(dòng)軌道動(dòng)力學(xué)可得瞬時(shí)軌道根數(shù)表示的瞬時(shí)軌道角速度和角加速度：

(5)

式中：a為軌道半長(zhǎng)軸；ec為偏心率；f為緯度幅角；Fx/r是軌道攝動(dòng)力所帶來(lái)的。

根據(jù)苛氏定理，在追蹤星軌道坐標(biāo)系下的兩星相對(duì)運(yùn)動(dòng)方程表示為：

(6)

(7)

式中：O(n2)為地球球形引力場(chǎng)的2階以上高階項(xiàng)引力差；Δatc為軌道環(huán)境攝動(dòng)力差。式(7)可改寫為：

(8)

若忽略ΔF的影響，上述方程是線性方程，即著名的Hill方程或CW方程，該方程被廣泛用于交會(huì)對(duì)接動(dòng)力學(xué)分析與控制問(wèn)題中。為了得到CW方程誤差的量化特性，有必要對(duì)ΔF的各組成部分進(jìn)行詳細(xì)研究。

2 相對(duì)運(yùn)動(dòng)模型誤差推導(dǎo)

分析可知，CW方程是近似方程，其誤差體現(xiàn)在加速度項(xiàng)ΔF。將ΔF分解為2部分：一部分是地球球形引力場(chǎng)的高階項(xiàng)和環(huán)境攝動(dòng)力差等未建模作用力ΔF1；另一部分是圓軌道假設(shè)下，由于實(shí)際偏心率的存在等帶來(lái)CW方程的未建模作用力ΔF2。

2.1 圓軌道假設(shè)下的建模誤差

首先，對(duì)地球球形引力場(chǎng)高階項(xiàng)和環(huán)境攝動(dòng)力差等未建模作用力ΔF1進(jìn)行量化分析，該部分誤差可改寫為：

ΔF1=O(n2)+Δatc

(9)

式中：Δatc表示空間環(huán)境中各攝動(dòng)力的綜合作用；O(n2)表示地球球形引力場(chǎng)兩星作用力差的二階及其高階項(xiàng)，若只考慮二階項(xiàng)，有：

(10)

其次，對(duì)于日月引力對(duì)相對(duì)運(yùn)動(dòng)的作用力參見式(2)。日月引力對(duì)兩星作用力的差消去共同部分，作用形式與地球球形引力場(chǎng)對(duì)兩星作用力的形式相同。同相對(duì)動(dòng)力學(xué)中地球球形引力場(chǎng)的一階近似過(guò)程(7)，在太陽(yáng)或月球中心慣性坐標(biāo)系下，對(duì)兩星的日月引力差取一階形式為：

(11)

式中：mt、mc表示在月球中心慣性坐標(biāo)系下目標(biāo)星和追蹤星的位置矢量；st、sc表示在太陽(yáng)中心慣性坐標(biāo)系下目標(biāo)星和追蹤星位置矢量；sx、sy、sz表示衛(wèi)星在太陽(yáng)軌道系下的相對(duì)位置坐標(biāo)；mx、my、mz表示衛(wèi)星在月球軌道系下的相對(duì)位置坐標(biāo)。

進(jìn)一步的，由表1可見，地球非球形引力攝動(dòng)中，各攝動(dòng)項(xiàng)形式相似，帶諧項(xiàng)J2項(xiàng)引力勢(shì)能作用力最大，較其他項(xiàng)高約3個(gè)數(shù)量級(jí)。在地心慣性系下，J2攝動(dòng)作用力表示為：

(12)

略去r-7以上高階項(xiàng)，式(12)可改寫為：

(13)

(14)

傾角i對(duì)于地球靜止軌道衛(wèi)星正常情況下是小量，小于0.1°，進(jìn)行一階展開得：

y=-sinisinΩ(xt-xc)-sinicosΩ(yt-yc)

-cosi(zt-zc)≈-(zt-zc)

(15)

將式(15)代入式(14)，J2項(xiàng)的攝動(dòng)差在追蹤星軌道坐標(biāo)系可表示為：

(16)

綜上所述，地球球形引力場(chǎng)的高階項(xiàng)和環(huán)境攝動(dòng)力差等未建模作用力ΔF1可統(tǒng)一表示為：

(17)

2.2 非圓軌道的相關(guān)誤差

由圓軌道假設(shè)等帶來(lái)CW方程的未建模誤差ΔF2表達(dá)式：

(18)

其中：

(19)

(20)

ΔF2表達(dá)式中有相對(duì)速度項(xiàng)，為了獲得相對(duì)位置表示的統(tǒng)一表達(dá)式。將速度項(xiàng)用相對(duì)位置表示。

在追蹤星軌道系下，目標(biāo)星相對(duì)運(yùn)動(dòng)的一階閉環(huán)形式，那么垂直于軌道平面運(yùn)動(dòng)，即軌道法向運(yùn)動(dòng)可以表示為：

(21)

易見，該方向運(yùn)動(dòng)為正弦周期振蕩，振蕩角速度為追蹤星瞬時(shí)軌道角速度。軌道面內(nèi)是橢圓運(yùn)動(dòng)，可表示為：

(22)

令Δl為一軌時(shí)間內(nèi)兩星距離變化量，zc為橢圓運(yùn)動(dòng)Z軸的中心坐標(biāo)，表示為：

(23)

因此，X軸相對(duì)速度可表示為：

(24)

將式(5)、(21)、(24)代入ΔF2中，忽略兩階小量得軌道面法向力攝動(dòng)為：

(25)

2.3 誤差綜合

統(tǒng)計(jì)各誤差的作用效果如表2，包含CW方程推導(dǎo)過(guò)程中引入的主要誤差源。地球球形引力場(chǎng)線性化中二次項(xiàng)誤差與衛(wèi)星高度的四次方成反比，與兩星的距離平方成正比；地球田諧項(xiàng)攝動(dòng)可以忽略；帶諧項(xiàng)攝動(dòng)中J2項(xiàng)誤差與衛(wèi)星高度五次方成反比，與兩星的距離成正比；日月引力的攝動(dòng)作用與衛(wèi)星高度相關(guān)性弱，與兩星的距離成正比；圓軌道假設(shè)中軌道偏心率相關(guān)誤差與衛(wèi)星高度的三次方成反比，與兩星的距離成正比；軌道面法向力攝動(dòng)通過(guò)影響軌道坐標(biāo)系的角速度耦合相對(duì)運(yùn)動(dòng)引入誤差，與衛(wèi)星高度的一次方成反比，與兩星的距離和一軌漂移量成正比。

表2 軌道攝動(dòng)力帶來(lái)的誤差Table 2 The error of the orbital perturbation force

綜上所述，CW方程的誤差一方面與衛(wèi)星高度有關(guān)，會(huì)隨衛(wèi)星高度的增加急速變??；另一方面與兩星的距離有關(guān)，包含線性相關(guān)和平方相關(guān)項(xiàng)。

3 地球靜止軌道相對(duì)運(yùn)動(dòng)模型誤差特性分析

本節(jié)對(duì)地球靜止軌道相對(duì)運(yùn)動(dòng)模型誤差進(jìn)行特性分析，得到其誤差的量級(jí)以及與運(yùn)動(dòng)構(gòu)型的關(guān)系，為后面章節(jié)的相對(duì)導(dǎo)航設(shè)計(jì)提供依據(jù)。

前節(jié)分析可得：日月引力攝動(dòng)、軌道面法向力攝動(dòng)、J2攝動(dòng)量級(jí)相同，且都與兩星的距離成正比；地球球形引力場(chǎng)二次項(xiàng)在兩星距離小于1 km時(shí)與J2等攝動(dòng)同階，在兩星距離大于1 km時(shí)攝動(dòng)力隨距離平方增加；偏心率相關(guān)攝動(dòng)大小與軌道偏心率成正比，在衛(wèi)星軌道偏心率小于1×10-5時(shí)攝動(dòng)力與J2等攝動(dòng)同階。CW方程在高軌誤差的作用形式變得復(fù)雜，相關(guān)因素較多，沒(méi)有明確的主要誤差，不易補(bǔ)償。

但對(duì)于典型的在軌飛行軌道，通常地球靜止軌道衛(wèi)星軌道偏心率小于0.001，兩星相對(duì)距離大于1 km，因此地球球形引力場(chǎng)的二次項(xiàng)和偏心率相關(guān)攝動(dòng)可以作為CW方程的主要誤差源，用以分析CW方程誤差的量級(jí)和作用形式。

利用軌道仿真軟件STK的HPOP軌道輸出，對(duì)主要攝動(dòng)和總攝動(dòng)的特性進(jìn)行數(shù)值分析。主要攝動(dòng)項(xiàng)包括地球引力二次項(xiàng)和偏心率相關(guān)攝動(dòng)，總攝動(dòng)在動(dòng)力學(xué)模型中按式(8)右端引出。

分析條件為：追蹤星偏心率0.001，兩星相距10 km，伴飛橢圓1 km。主要攝動(dòng)項(xiàng)的量化特性見圖2，總攝動(dòng)項(xiàng)的量化特性見圖3。

圖2 地球靜止軌道主要攝動(dòng)累加Fig.2 The main perturbation of the geostationary orbit

圖3 地球靜止軌道CW方程各軸總攝動(dòng)Fig.3 The total perturbation of each axis of the CW equation of geostationary orbit

仿真可見，上述主要攝動(dòng)項(xiàng)在各類誤差中占主要部分，與總誤差量級(jí)相同，且振蕩和常偏特性類似，能夠反映地球靜止軌道相對(duì)運(yùn)動(dòng)模型CW方程的誤差特性。圖3中的波動(dòng)部分是目標(biāo)進(jìn)出陰影過(guò)程中的太陽(yáng)光壓攝動(dòng)差。

下面對(duì)地球引力二次項(xiàng)和偏心率相關(guān)攝動(dòng)項(xiàng)的量化特性進(jìn)行獨(dú)立分析。圖4是對(duì)表2中，地球球形引力二階項(xiàng)在兩星10 km穩(wěn)定伴飛工況下各軸作用力的量化特性。相對(duì)運(yùn)動(dòng)橢圓和軌道面外振蕩控制的比較小，此時(shí)地球球形引力場(chǎng)二階項(xiàng)在軌道面法向Y軸作用力很小，該項(xiàng)與兩星的相對(duì)距離沒(méi)有關(guān)系；地球球形引力場(chǎng)的二階項(xiàng)對(duì)X軸和Z軸的攝動(dòng)作用與兩星相對(duì)距離和伴飛橢圓尺寸均相關(guān)，由于相對(duì)運(yùn)動(dòng)在Z軸為周期運(yùn)動(dòng)，X軸攝動(dòng)項(xiàng)表現(xiàn)出周期振蕩；Z軸的攝動(dòng)作用主要特點(diǎn)是常偏，在穩(wěn)定伴飛期間，由于Y軸振蕩小，Z項(xiàng)的攝動(dòng)表現(xiàn)為跟兩星距離平方成正比的常偏和一個(gè)與伴飛橢圓相關(guān)的小幅振蕩的綜合。

圖4 地球球形引力場(chǎng)二階項(xiàng)攝動(dòng)各軸分量Fig.4 Second order perturbation of the earth′s spherical gravitational field perturbs each axial component

圖5是對(duì)表2中偏心率相關(guān)項(xiàng)在兩星10 km穩(wěn)定伴飛工況下各軸作用力的仿真輸出，呈現(xiàn)周期性的振蕩特性。從相對(duì)運(yùn)動(dòng)方程CW方程的穩(wěn)態(tài)解析式(22)可得，各軸相對(duì)位置呈現(xiàn)周期性變化，偏心率相關(guān)項(xiàng)各軸攝動(dòng)耦合因素多，偏心率和相對(duì)運(yùn)動(dòng)橢圓尺寸聯(lián)合影響振蕩的幅值，相對(duì)運(yùn)動(dòng)的當(dāng)前位置和當(dāng)前軌道幅角等綜合影響振蕩的相位。

圖5 偏心率相關(guān)攝動(dòng)各軸分量Fig.5 Eccentricity correlative perturbation of each axis component

綜合分析表明，地球球形引力二次項(xiàng)和偏心率相關(guān)攝動(dòng)項(xiàng)反映了CW方程在地球靜止軌道的主要誤差特性，可以在相對(duì)導(dǎo)航濾波設(shè)計(jì)中用以定量分析狀態(tài)方程的誤差情況，或設(shè)計(jì)干擾估計(jì)模型。地球球形引力場(chǎng)二次項(xiàng)的誤差在各軸呈現(xiàn)周期振蕩特性，并隨兩星距離變大，在Z軸產(chǎn)生常值偏差分量；偏心率相關(guān)攝動(dòng)與兩星相對(duì)構(gòu)型線性相關(guān)，各軸均呈現(xiàn)為周期振蕩特性。

5 結(jié)論

1)考慮地球球形引力場(chǎng)高階項(xiàng)、地球非球形引力攝動(dòng)、三體引力攝動(dòng)和太陽(yáng)光壓攝動(dòng)，分析了地球靜止軌道環(huán)境中各攝動(dòng)力的對(duì)相對(duì)運(yùn)動(dòng)的作用量級(jí)，對(duì)兩星相對(duì)運(yùn)動(dòng)方程進(jìn)行了建模分析，得到了完整的誤差要素。

2)為研究CW方程的誤差特性，對(duì)加速度項(xiàng)進(jìn)行深入分析，給出了各誤差源的近似表達(dá)式、作用量級(jí)和影響因素，分析表明，CW方程的誤差隨衛(wèi)星高度的增加急速變小，并且與相對(duì)距離存在線性和平方復(fù)合的相關(guān)特性。

3)對(duì)典型的相對(duì)運(yùn)動(dòng)型態(tài)進(jìn)行了數(shù)值仿真，結(jié)果表明地球球形引力二次項(xiàng)、偏心率相關(guān)項(xiàng)以及J2項(xiàng)相關(guān)誤差為CW方程的主要誤差來(lái)源，在軌道系坐標(biāo)系各軸向具體呈現(xiàn)為常值偏差和軌道周期振蕩組合的分布特性。