遍歷理論在航天動力學領(lǐng)域的應(yīng)用*

白 雪，徐 明，姚 闖

0 引 言

近年來，微分幾何和動力系統(tǒng)理論的發(fā)展為非線性動力學提供有力的分析工具；同時也促使航天動力學等傳統(tǒng)領(lǐng)域基于現(xiàn)代數(shù)學觀點反思其理論框架，并得出一些超出經(jīng)典理論的創(chuàng)新性結(jié)論[1-2].此外，航天動力學以工程應(yīng)用為背景，需要借助動力系統(tǒng)等相關(guān)數(shù)學理論進行嚴謹化處理，以驗證或改進相關(guān)經(jīng)驗公式和方法.

軌道設(shè)計的一個重要任務(wù)是使得衛(wèi)星滿足一定的覆蓋特性和日照時間等約束；而除特定軌道外，星下點軌跡一般不會重復出現(xiàn)，特別在考慮攝動影響后，解析分析覆蓋頻率和受曬因子等基本不可能.計算平均化參數(shù)往往需要統(tǒng)計長時間內(nèi)(5-10年)衛(wèi)星的累積量，這將耗費巨大的計算量，甚至不可實現(xiàn).為解決上述問題，本文將研究遍歷理論在航天動力學領(lǐng)域的應(yīng)用，將無限時間的累積問題成功轉(zhuǎn)換為有限空間的積分運算.

LO[3]最早提出應(yīng)用遍歷理論計算長期覆蓋，將復雜的微分方程求解轉(zhuǎn)化為簡單的定積分問題.對于圓軌道，所給出的預(yù)報誤差不超出0.2%，但對于橢圓軌道預(yù)報誤差超出16%.而且，LO在其一系列文章僅描繪了算法構(gòu)造的思路，但缺乏嚴格的數(shù)學證明.LO所發(fā)展的算法無法解決橢圓軌道不同于圓軌道的根本原因，即問題的升維[4].LO針對橢圓軌道覆蓋構(gòu)造的不變測度僅限于抽象動力系統(tǒng)，且需要計算向量場的散度，進而數(shù)值求解一組微分方程，計算量將會很大，很難應(yīng)用于本問題.XU等[5]對橢圓軌道長期覆蓋問題提出了有效的遍歷算法，但未見到其應(yīng)用遍歷理論在其他方面的計算研究.洪元等[6]利用太陽相對于軌道面入射角計算受曬因子，以判斷軌道極端工況出現(xiàn)的日期，但僅限于一年內(nèi)的短期受曬因子計算，對于長期受曬因子并沒有合適的計算方法.郗曉寧等[7]利用球面三角形推導了圓軌道受曬因子的計算公式，但并未對橢圓軌道情況進行討論，存在一定局限性.因此，本文應(yīng)用Birkhoff-Khinchin定理將長期覆蓋和受曬因子等無限時間的累積問題成功轉(zhuǎn)換為有限空間的二重積分運算，以此來解決上述研究存在的局限性問題.

本文將系統(tǒng)地研究遍歷理論在航天動力學領(lǐng)域的應(yīng)用，包括橢圓軌道長期覆蓋和受曬因子等長期參數(shù)的平均化計算；基于動力系統(tǒng)理論構(gòu)造上述問題的不變測度、相空間和相流、實值函數(shù)等遍歷等價條件，從而完成遍歷算法的數(shù)學嚴謹化證明.研究結(jié)果表明，遍歷理論是解決航天動力學領(lǐng)域長期參數(shù)計算的有力工具.

1 遍歷理論

遍歷理論源于統(tǒng)計力學中的Boltzmann遍歷假設(shè)，主張用統(tǒng)計的觀點來研究確定性動力系統(tǒng).遍歷理論研究群對測度空間的保測作用[5-6].遍歷理論中最基本且最重要的是Birkhoff-Khinchin定理[8-9].

設(shè)動力系統(tǒng)(M,μ,φt),其相空間定義為微分流形M；保測單參數(shù)變換群φt,稱為相流；不變測度μ，使得μ(φt(x))=μ,?x∈M,且μ(M)=1.φt的任一不變集，若其測度只可能為1或0，則稱該保測變換φt遍歷.

設(shè)M上μ可加的實(或復)值函數(shù)f∈L1(M,μ),定義f時間平均為

(1)

空間平均為

(2)

定理1(Birkhoff-Khinchin)[3-4].設(shè)動力系統(tǒng)(M,μ,φt),f∈L1(M,μ),則有

(3)

即，f時間平均等于其空間平均.

Birkhoff-Khinchin定理的詳細討論及證明參見文獻[3]和[4].應(yīng)用遍歷理論，可以將特定的時間積分轉(zhuǎn)化為空間積分.對于本文中涉及到的長期覆蓋和受曬因子計算，采用空間積分將很大程度減少計算量，使計算更加容易.

2 基于遍歷理論的長期覆蓋預(yù)報算法

2.1 橢圓軌道

航天器運行于Kepler軌道上，因攝動(本文只考慮J2項)影響，其軌道根數(shù)均不是常數(shù).考慮覆蓋的長期計算，可以用平均軌道根數(shù)來描述其軌道特征.在平均軌道根數(shù)描述下，半長軸a、偏心率e、軌道傾角i均為常數(shù)；升交點赤經(jīng)Ω、近地點幅角ω、平近點角Ma變化如下：

(4)

式中，CΩ,Cω,CMα為一周期內(nèi)的平均變化率，顯然是J2的函數(shù)[10].因此，平均軌道根數(shù)Ω,ω,Mα將唯一確定軌道的狀態(tài).記為(Ω,ω,Mα)∈T3,其中T3=S1×S1×S1為3維環(huán)面.

若CΩ,Cω,CMα有理無關(guān)，則式(4)描述的軌線在T3上致密分布且可遍歷；若CΩ,Cω,CMα有理相關(guān)，則軌線在T3上存在共振層且不再遍歷[8].

記集合

S1={(CΩ,Cω,CM)∈R3|CΩ,Cω,CM有理無關(guān)}

S2={(CΩ,Cω,CM)∈R3|CΩ,Cω,CM有理相關(guān)}

顯然，S1為全測度(Lebesgue測度)集,S2為零測度集.在概率意義上，?(CΩ,Cω,CMα)∈R3,有如下關(guān)系：

(5)

其中，P[A]為事件A發(fā)生的概率.

故本文不變測度的構(gòu)造，可以僅考慮CΩ,Cω,CMα有理無關(guān)的情況.

星下點的緯度φ和經(jīng)度θ，極徑r及軌道根數(shù)有如下關(guān)系[10]：

(6)

(7)

其中，?為真近點角，ωe為地球自轉(zhuǎn)角速度.

2.2 相流與不變測度

與圓軌道不同，橢圓軌道半徑的變化r∈[rmin,rmax]=[a(1-e),a(1+e)]將會使覆蓋描述升維[4]，圓軌道的相流(星下點)及相空間(地球面)已不再適用于橢圓軌道；橢圓軌道的相空間M定義為半徑分別為rmin和rmax的球面所圍的球面環(huán)；相流φt為航天器在M中的真實運動軌跡，是關(guān)于時間t的單參數(shù)變換群.

構(gòu)造M上的非負函數(shù)μ的微分形式為

dμ=κdΩ∧dω∧dMα

(8)

其中，κ為待定常數(shù)，∧為外積.

引理1.取Σ為M的冪集，則κ若取合適值時，{M,Σ,μ}構(gòu)成概率空間.

證明.Σ為M的冪集，即Σ是M上的σ代數(shù)，故對NM，有NΣ成立.

根據(jù)式(8)的定義，在N上取Lebesgue積分，可得

(9)

(10)

式(10)表明，μ(N)的值可由N上的Lebesgue積分得到.由Lebesgue積分的可數(shù)可加性，可以得到μ具有可數(shù)可加性；由式(10)同樣可以得到μ(φ)=0.故μ構(gòu)成Σ上的測度.

定理2(保測變換)[11].令S為有限測度空間(X,Π,v)上的一個可測變換，則為S-不變的當且僅當對所有的h=L∞(X),有

(11)

該定理的詳細討論及證明參見文獻[11].

引理2.單參數(shù)變換群φt:M→M關(guān)于μ是保測的.

證明.由于Σ為M的冪集，故對?A∈Σ,有φt-1(A)∈Σ恒成立，即φt是可測的.

對?g(Ω,ω,Mα)∈L∞(X),有

(12)

考慮φt的作用，即式(4)，可得

d(Ω+CΩt)d(ω+Cωt)d(Ma+CMat)

(13)

故由定理2可知，φt是保測的.

2.3 實值函數(shù)

航天器覆蓋角λ與極徑r的關(guān)系[10]為：

(14)

其中，γ為允許最小視角，Re為地球半徑，見圖1所示.

不同的軌道半徑對應(yīng)不同的覆蓋角，r在[rmin,rmax]內(nèi)變化所圍成的區(qū)域(station mask，SM)為地面站的可見區(qū)域，見圖2所示.

地面站對航天器的可見頻率ρ,定義為

(15)

其中，P(T)為0到T時間內(nèi)航天器處于地面站可見區(qū)域的總時間.

根據(jù)式(1)及ρ的定義，定義實值函數(shù)f為SM的特征函數(shù)，即

(16)

則可得

(17)

由于SM?M,容易驗證f為M上μ可加的實值函數(shù)且f∈L1(M,μ),即f滿足Birkhoff-Khinchin定理條件.

2.4 覆蓋的遍歷算法

2.4.1κ值的確定

顯然，對于本文的M及N(=SM)，其上的Lebesgue積分與Riemannian積分在數(shù)值上是相等的.

微分F1，可得

(18)

(19)

在M上展開式(19)，可得

(20)

進而可得

κ=1/2π3

(21)

2.4.2 長期覆蓋的遍歷算法

上面定義的動力系統(tǒng)(M,μ,φt)和實值函數(shù)f滿足Birkhoff-Khinchin定理的條件，故可得地面站對航天器的可見頻率ρ為

(22)

為了便于計算，應(yīng)用變換F1和F2，得到三重積分式：

(23)

展開上式，可得

(24)

其中，φ1、φ2、g(φ)的表達式為[3]：

φ1=max(φ0-λ，-Li)

(25)

φ2=min(φ0+λ,Li)

(26)

(27)

而φ0為地面站緯度，覆蓋角λ由式(14)給出，r由式(7)給出.

進一步化簡三重積分式(24)，可得二重積分式：

(28)

由此可得如下定理：

定理3.地面站對橢圓軌道航天器的可見頻率ρ，可以轉(zhuǎn)化為二重積分式(28)計算，且這個性質(zhì)幾乎處處成立.

由于共振層的出現(xiàn)，式(28)不適用于零測度集S2.故針對集合S2的遍歷算法有待于進一步的研究.

2.4.3 退化情況

當橢圓軌道退化為圓軌道時，覆蓋頻率ρC為

(29)

(30)

顯然，式(30)的表達式與文獻[3]的計算公式相同.但這里相空間為半徑分別為軌道半徑的球面，相流為航天器在該球面中的真實運動軌跡；而文獻[3]定義的相空間為地球面，相流為星下點軌跡.兩者不同，對于圓軌道，可以證明兩者是拓撲等價的.

2.4.4 數(shù)值驗證

由于圓軌道是橢圓軌道的特殊情況，為了驗證所得結(jié)果的正確性，僅需驗證文獻[3-4]覆蓋預(yù)報誤差較大的情況.橢圓軌道的偏心率均取為0.05，允許最小視角γ取為0°，理論結(jié)果是應(yīng)用復合Simpson公式通過計算式(22)得到的，而數(shù)值結(jié)果則是通過軌道積分得到的(仿真時間為5年)，比較結(jié)果見表1.顯然，利用本文推導的公式得到橢圓軌道的覆蓋預(yù)報誤差不超出0.2%，即明顯改進了LO的結(jié)果[3-4].受曬因子數(shù)值仿真結(jié)果與長期覆蓋結(jié)果類似，五年內(nèi)的誤差均不超過0.2%.

表1 橢圓軌道覆蓋的數(shù)值與理論結(jié)果比較(偏心率為0.05)Tab.1 The numerical and theoretical result comparison of the elliptical orbitcoverage(the eccentricity is 0.05)

3 遍歷理論在受曬因子計算中的應(yīng)用

太陽光能是星上能源的重要來源.圓軌道日光受曬因子用ks表示，其定義為航天器繞地球一周中受太陽照射時掃過的地心角與2π的比值.正常工作的衛(wèi)星，需要在其壽命內(nèi)滿足一定的日光受曬因子.因此，長期受曬因子也是軌道設(shè)計的重要指標.

通常，航天動力學教科書只給出了短期受曬因子的計算，即圓軌道的解析算法和橢圓軌道的數(shù)值算法[10].然而，隨著太陽方位角的變化，短期受曬因子不再是常值；而且，衛(wèi)星軌道在各種攝動(主要是J2攝動)下的變形，使得受曬因子的計算極其復雜；而長期受曬因子的計算，需要進行長時間(5～10年)的軌道積分，這將耗費巨大的計算資源.

幸運的是，這類時間平均的問題，在某些情況下可以轉(zhuǎn)化為對空間的平均問題.LO[3-4]應(yīng)用遍歷理論研究了圓軌道的長期覆蓋，將微分方程的求解轉(zhuǎn)化為一重定積分；第2節(jié)填補了LO的理論缺陷，并將覆蓋的遍歷算法成功擴展到橢圓軌道.本節(jié)將把遍歷理論應(yīng)用到長期受曬因子的計算，并針對圓和橢圓軌道分別給出具體算法.

3.1 短期受曬因子(The short-term non-eclipse factor，SNF)

日光照射地球時，在地球背向日光的一面將產(chǎn)生地影，當衛(wèi)星飛進地影時，將不受日光照射，此時稱為星蝕；反之，衛(wèi)星將受日光照射，稱為受曬.由于太陽尺寸遠大于地球，地球在日光照射下的地影由本影和半影區(qū)組成.地球本影區(qū)為一圓錐體，對于大部分軌道高度的衛(wèi)星，該錐體的錐角很小，可近似將日光看作平行光.

對于短期受曬因子的計算，可以忽略地球繞太陽的公轉(zhuǎn)，僅需考慮衛(wèi)星繞地球的運動.因此，給定太陽的方位，便可以判斷出相應(yīng)的陰影區(qū)(shadow region，SR)和受曬區(qū)(lighting region，LR).

短期受曬因子為動力系統(tǒng)的瞬時值，可有瞬時軌道根數(shù)獲得；因此，其算法對無攝和有攝軌道均適用.

圓軌道存在SNF的解析算法，其短期受曬因子kS的計算公式如下：

(31)

式中α和η定義為

(32)

其中，h為軌道高度，Re為地球半徑，i為軌道傾角，而iS=23.439°為黃赤交角，φ為太陽赤經(jīng)，Ω為升交點赤經(jīng).

顯然，由物理直觀和式(31)都可得到SNF恒大于0.5的事實.

對于無攝的圓軌道來說，判斷衛(wèi)星是否位于SR區(qū)域，僅需已知太陽赤經(jīng)φ和軌道幅角u的數(shù)值.這里取軌道高度為200 km、軌道傾角為61°、升交點赤徑為0°，短期受曬因子的φ-u分布如圖3所示.圖中黑色區(qū)域為SR，白色區(qū)域為LR.方位角為φ0的kS值為φ=φ0所在列中位于LR內(nèi)u值總和與360°的比值.

對于給定軌道，短期受曬因子僅為φ的函數(shù)，圖4給出了SNF與φ的關(guān)系.對于圓軌道，SNF值ks(φ)是φ的周期函數(shù).

橢圓軌道不存在受曬因子的解析表達式，一般需要應(yīng)用數(shù)值算法求解.通過坐標變換，可以得到衛(wèi)星在太陽觀測坐標系下的表示，即(xv,yv,zv)，其中太陽觀測坐標系(Sun-observer frame，SOF)定義如下：原點取為地心、x軸由地球指向太陽、z軸垂直于黃道面、y軸由右手法則確定.

當xv<0且yv2+zv2

3.2 長期受曬因子(The Long-term Non-eclipse Factor，LNF)

長期受曬因子LNF，定義為

(33)

其中P(T)為衛(wèi)星處于太陽照射下的總時間.

顯然，LNF為SNF的統(tǒng)計性結(jié)果，并可應(yīng)用一些平均化方法以簡化算法.對于非攝動圓軌道，ks(φ)的周期性和圓軌道的均一性表明：LNF可以理解為SNF在[0,TOrbit]范圍內(nèi)的時間平均，即

(34)

但攝動和偏心率的存在將使得ρ十分復雜，一般需要長時間(5～10年或更長)數(shù)值積分.

3.3 長期受曬因子的遍歷算法

根據(jù)第2.3節(jié)的分析，在J2攝動的影響下，不變測度可構(gòu)造為

dμ=κdΩ∧dω∧dMα∧dφ

(35)

為了得到LNF，定義實值函數(shù)f為LR和SR上的特征函數(shù)，即

(36)

相流定義為衛(wèi)星在太陽觀察坐標系下的運動軌跡，相空間取為R3.

根據(jù)Birkhoff-Khinchin定理，長期受曬因子ρ可表示為

(37)

常數(shù)κ基于μ(M)=1的事實獲得，即

(38)

3.3.1 無攝圓軌道受曬因子的遍歷算法

對于無攝圓軌道，僅需φ和u就可以判斷衛(wèi)星是否處于SR，因此不變測度將由4維退化為2維，相應(yīng)的不變測度的微元為

dμ=κdu∧dφ

(39)

由式(38)可以得到κ的值為

κ=1/(2π)2

(40)

則長期受曬因子ρ為

(41)

顯然，該結(jié)果與式(34)一致，這也驗證了前節(jié)時間平均方法的有效性.

3.3.2 受攝圓軌道受曬因子的遍歷算法

對于受攝圓軌道，升交點赤徑Ω不再是常值；因此，判斷衛(wèi)星是否位于LR區(qū)域不僅需要φ和u的值還需要Ω的值.公式(31)仍然有效，不過需要考慮獨立變量Ω的影響，即改進為ks(φ,Ω).

根據(jù)式(38)，可得κ的值為

κ=1/(2π)3

(42)

長期受曬因子ρ為

(43)

3.3.3 受攝橢圓軌道受曬因子的遍歷算法

對于受攝橢圓軌道，軌道半徑不再是常值，判斷衛(wèi)星是否位于LR區(qū)域?qū)⑿枰?、u、Ω和r.

前節(jié)已經(jīng)詳細敘述橢圓軌道衛(wèi)星是否處于LR區(qū)域的判斷方法.這里記

(44)

變量u和r有內(nèi)在的關(guān)聯(lián)性，故φ、u、Ω和r不具有遍歷性；而變量φ、Ma、Ω和近地點幅角ω具有遍歷性，幸運的是，兩組變量存在如下變換：

(45)

(46)

3.3.4 無攝橢圓軌道受曬因子的遍歷算法

無攝條件下Ω和保持不變，僅需Ma和φ即可判斷衛(wèi)星的受曬狀態(tài)，則不變測度μ將退化為

dμ=κ·dMα∧dφ

(47)

長期受曬因子ρ為

(48)

4 結(jié) 論

本文系統(tǒng)地研究遍歷理論在航天動力學領(lǐng)域的應(yīng)用，包括橢圓軌道長期覆蓋和受曬因子等長期參數(shù)的平均化計算.從無攝和J2項攝動下平均軌道根數(shù)的長期和長周期項出發(fā)，構(gòu)造非共振條件下覆蓋和受曬因子的測度描述，并證明該測度的不變性；給出相關(guān)相空間和相流定義，并證明該相流的保測度性；應(yīng)用Birkhoff-Khinchin定理將長期覆蓋和受曬因子等無限時間的累積問題成功轉(zhuǎn)換為有限空間的二重積分運算.

研究結(jié)果表明，遍歷理論是解決航天動力學領(lǐng)域長期參數(shù)計算的有力工具；借助動力系統(tǒng)等相關(guān)數(shù)學理論進行嚴謹化處理，可驗證或改進相關(guān)經(jīng)驗公式和方法，并提高經(jīng)典軌道理論中覆蓋率與受曬因子的計算效率.