遍历理论在航天动力学领域的应用*

白 雪，徐 明，姚 闯

0 引 言

近年来，微分几何和动力系统理论的发展为非线性动力学提供有力的分析工具；同时也促使航天动力学等传统领域基于现代数学观点反思其理论框架，并得出一些超出经典理论的创新性结论[1-2].此外，航天动力学以工程应用为背景，需要借助动力系统等相关数学理论进行严谨化处理，以验证或改进相关经验公式和方法.

轨道设计的一个重要任务是使得卫星满足一定的覆盖特性和日照时间等约束；而除特定轨道外，星下点轨迹一般不会重复出现，特别在考虑摄动影响后，解析分析覆盖频率和受晒因子等基本不可能.计算平均化参数往往需要统计长时间内(5-10年)卫星的累积量，这将耗费巨大的计算量，甚至不可实现.为解决上述问题，本文将研究遍历理论在航天动力学领域的应用，将无限时间的累积问题成功转换为有限空间的积分运算.

LO[3]最早提出应用遍历理论计算长期覆盖，将复杂的微分方程求解转化为简单的定积分问题.对于圆轨道，所给出的预报误差不超出0.2%，但对于椭圆轨道预报误差超出16%.而且，LO在其一系列文章仅描绘了算法构造的思路，但缺乏严格的数学证明.LO所发展的算法无法解决椭圆轨道不同于圆轨道的根本原因，即问题的升维[4].LO针对椭圆轨道覆盖构造的不变测度仅限于抽象动力系统，且需要计算向量场的散度，进而数值求解一组微分方程，计算量将会很大，很难应用于本问题.XU等[5]对椭圆轨道长期覆盖问题提出了有效的遍历算法，但未见到其应用遍历理论在其他方面的计算研究.洪元等[6]利用太阳相对于轨道面入射角计算受晒因子，以判断轨道极端工况出现的日期，但仅限于一年内的短期受晒因子计算，对于长期受晒因子并没有合适的计算方法.郗晓宁等[7]利用球面三角形推导了圆轨道受晒因子的计算公式，但并未对椭圆轨道情况进行讨论，存在一定局限性.因此，本文应用Birkhoff-Khinchin定理将长期覆盖和受晒因子等无限时间的累积问题成功转换为有限空间的二重积分运算，以此来解决上述研究存在的局限性问题.

本文将系统地研究遍历理论在航天动力学领域的应用，包括椭圆轨道长期覆盖和受晒因子等长期参数的平均化计算；基于动力系统理论构造上述问题的不变测度、相空间和相流、实值函数等遍历等价条件，从而完成遍历算法的数学严谨化证明.研究结果表明，遍历理论是解决航天动力学领域长期参数计算的有力工具.

1 遍历理论

遍历理论源于统计力学中的Boltzmann遍历假设，主张用统计的观点来研究确定性动力系统.遍历理论研究群对测度空间的保测作用[5-6].遍历理论中最基本且最重要的是Birkhoff-Khinchin定理[8-9].

设动力系统(M,μ,φt),其相空间定义为微分流形M；保测单参数变换群φt,称为相流；不变测度μ，使得μ(φt(x))=μ,∀x∈M,且μ(M)=1.φt的任一不变集，若其测度只可能为1或0，则称该保测变换φt遍历.

设M上μ可加的实(或复)值函数f∈L1(M,μ),定义f时间平均为

(1)

空间平均为

(2)

定理1(Birkhoff-Khinchin)[3-4].设动力系统(M,μ,φt),f∈L1(M,μ),则有

(3)

即，f时间平均等于其空间平均.

Birkhoff-Khinchin定理的详细讨论及证明参见文献[3]和[4].应用遍历理论，可以将特定的时间积分转化为空间积分.对于本文中涉及到的长期覆盖和受晒因子计算，采用空间积分将很大程度减少计算量，使计算更加容易.

2 基于遍历理论的长期覆盖预报算法

2.1 椭圆轨道

航天器运行于Kepler轨道上，因摄动(本文只考虑J2项)影响，其轨道根数均不是常数.考虑覆盖的长期计算，可以用平均轨道根数来描述其轨道特征.在平均轨道根数描述下，半长轴a、偏心率e、轨道倾角i均为常数；升交点赤经Ω、近地点幅角ω、平近点角Ma变化如下：

(4)

式中，CΩ,Cω,CMα为一周期内的平均变化率，显然是J2的函数[10].因此，平均轨道根数Ω,ω,Mα将唯一确定轨道的状态.记为(Ω,ω,Mα)∈T3,其中T3=S1×S1×S1为3维环面.

若CΩ,Cω,CMα有理无关，则式(4)描述的轨线在T3上致密分布且可遍历；若CΩ,Cω,CMα有理相关，则轨线在T3上存在共振层且不再遍历[8].

记集合

S1={(CΩ,Cω,CM)∈R3|CΩ,Cω,CM有理无关}

S2={(CΩ,Cω,CM)∈R3|CΩ,Cω,CM有理相关}

显然，S1为全测度(Lebesgue测度)集,S2为零测度集.在概率意义上，∀(CΩ,Cω,CMα)∈R3,有如下关系：

(5)

其中，P[A]为事件A发生的概率.

故本文不变测度的构造，可以仅考虑CΩ,Cω,CMα有理无关的情况.

星下点的纬度φ和经度θ，极径r及轨道根数有如下关系[10]：

(6)

(7)

其中，ϑ为真近点角，ωe为地球自转角速度.

2.2 相流与不变测度

与圆轨道不同，椭圆轨道半径的变化r∈[rmin,rmax]=[a(1-e),a(1+e)]将会使覆盖描述升维[4]，圆轨道的相流(星下点)及相空间(地球面)已不再适用于椭圆轨道；椭圆轨道的相空间M定义为半径分别为rmin和rmax的球面所围的球面环；相流φt为航天器在M中的真实运动轨迹，是关于时间t的单参数变换群.

构造M上的非负函数μ的微分形式为

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

(8)

其中，κ为待定常数，∧为外积.

引理1.取Σ为M的幂集，则κ若取合适值时，{M,Σ,μ}构成概率空间.

证明.Σ为M的幂集，即Σ是M上的σ代数，故对NM，有NΣ成立.

根据式(8)的定义，在N上取Lebesgue积分，可得

(9)

(10)

式(10)表明，μ(N)的值可由N上的Lebesgue积分得到.由Lebesgue积分的可数可加性，可以得到μ具有可数可加性；由式(10)同样可以得到μ(φ)=0.故μ构成Σ上的测度.

定理2(保测变换)[11].令S为有限测度空间(X,Π,v)上的一个可测变换，则为S-不变的当且仅当对所有的h=L∞(X),有

(11)

该定理的详细讨论及证明参见文献[11].

引理2.单参数变换群φt:M→M关于μ是保测的.

证明.由于Σ为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 实值函数

航天器覆盖角λ与极径r的关系[10]为：

(14)

其中，γ为允许最小视角，Re为地球半径，见图1所示.

不同的轨道半径对应不同的覆盖角，r在[rmin,rmax]内变化所围成的区域(station mask，SM)为地面站的可见区域，见图2所示.

地面站对航天器的可见频率ρ,定义为

(15)

其中，P(T)为0到T时间内航天器处于地面站可见区域的总时间.

根据式(1)及ρ的定义，定义实值函数f为SM的特征函数，即

(16)

则可得

(17)

由于SM⊂M,容易验证f为M上μ可加的实值函数且f∈L1(M,μ),即f满足Birkhoff-Khinchin定理条件.

2.4 覆盖的遍历算法

2.4.1κ值的确定

显然，对于本文的M及N(=SM)，其上的Lebesgue积分与Riemannian积分在数值上是相等的.

微分F1，可得

(18)

(19)

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

(20)

进而可得

κ=1/2π3

(21)

2.4.2 长期覆盖的遍历算法

上面定义的动力系统(M,μ,φt)和实值函数f满足Birkhoff-Khinchin定理的条件，故可得地面站对航天器的可见频率ρ为

(22)

为了便于计算，应用变换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.地面站对椭圆轨道航天器的可见频率ρ，可以转化为二重积分式(28)计算，且这个性质几乎处处成立.

由于共振层的出现，式(28)不适用于零测度集S2.故针对集合S2的遍历算法有待于进一步的研究.

2.4.3 退化情况

当椭圆轨道退化为圆轨道时，覆盖频率ρC为

(29)

(30)

显然，式(30)的表达式与文献[3]的计算公式相同.但这里相空间为半径分别为轨道半径的球面，相流为航天器在该球面中的真实运动轨迹；而文献[3]定义的相空间为地球面，相流为星下点轨迹.两者不同，对于圆轨道，可以证明两者是拓扑等价的.

2.4.4 数值验证

由于圆轨道是椭圆轨道的特殊情况，为了验证所得结果的正确性，仅需验证文献[3-4]覆盖预报误差较大的情况.椭圆轨道的偏心率均取为0.05，允许最小视角γ取为0°，理论结果是应用复合Simpson公式通过计算式(22)得到的，而数值结果则是通过轨道积分得到的(仿真时间为5年)，比较结果见表1.显然，利用本文推导的公式得到椭圆轨道的覆盖预报误差不超出0.2%，即明显改进了LO的结果[3-4].受晒因子数值仿真结果与长期覆盖结果类似，五年内的误差均不超过0.2%.

表1 椭圆轨道覆盖的数值与理论结果比较(偏心率为0.05)Tab.1 The numerical and theoretical result comparison of the elliptical orbitcoverage(the eccentricity is 0.05)

3 遍历理论在受晒因子计算中的应用

太阳光能是星上能源的重要来源.圆轨道日光受晒因子用ks表示，其定义为航天器绕地球一周中受太阳照射时扫过的地心角与2π的比值.正常工作的卫星，需要在其寿命内满足一定的日光受晒因子.因此，长期受晒因子也是轨道设计的重要指标.

通常，航天动力学教科书只给出了短期受晒因子的计算，即圆轨道的解析算法和椭圆轨道的数值算法[10].然而，随着太阳方位角的变化，短期受晒因子不再是常值；而且，卫星轨道在各种摄动(主要是J2摄动)下的变形，使得受晒因子的计算极其复杂；而长期受晒因子的计算，需要进行长时间(5～10年)的轨道积分，这将耗费巨大的计算资源.

幸运的是，这类时间平均的问题，在某些情况下可以转化为对空间的平均问题.LO[3-4]应用遍历理论研究了圆轨道的长期覆盖，将微分方程的求解转化为一重定积分；第2节填补了LO的理论缺陷，并将覆盖的遍历算法成功扩展到椭圆轨道.本节将把遍历理论应用到长期受晒因子的计算，并针对圆和椭圆轨道分别给出具体算法.

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

日光照射地球时，在地球背向日光的一面将产生地影，当卫星飞进地影时，将不受日光照射，此时称为星蚀；反之，卫星将受日光照射，称为受晒.由于太阳尺寸远大于地球，地球在日光照射下的地影由本影和半影区组成.地球本影区为一圆锥体，对于大部分轨道高度的卫星，该锥体的锥角很小，可近似将日光看作平行光.

对于短期受晒因子的计算，可以忽略地球绕太阳的公转，仅需考虑卫星绕地球的运动.因此，给定太阳的方位，便可以判断出相应的阴影区(shadow region，SR)和受晒区(lighting region，LR).

短期受晒因子为动力系统的瞬时值，可有瞬时轨道根数获得；因此，其算法对无摄和有摄轨道均适用.

圆轨道存在SNF的解析算法，其短期受晒因子kS的计算公式如下：

(31)

式中α和η定义为

(32)

其中，h为轨道高度，Re为地球半径，i为轨道倾角，而iS=23.439°为黄赤交角，φ为太阳赤经，Ω为升交点赤经.

显然，由物理直观和式(31)都可得到SNF恒大于0.5的事实.

对于无摄的圆轨道来说，判断卫星是否位于SR区域，仅需已知太阳赤经φ和轨道幅角u的数值.这里取轨道高度为200 km、轨道倾角为61°、升交点赤径为0°，短期受晒因子的φ-u分布如图3所示.图中黑色区域为SR，白色区域为LR.方位角为φ0的kS值为φ=φ0所在列中位于LR内u值总和与360°的比值.

对于给定轨道，短期受晒因子仅为φ的函数，图4给出了SNF与φ的关系.对于圆轨道，SNF值ks(φ)是φ的周期函数.

椭圆轨道不存在受晒因子的解析表达式，一般需要应用数值算法求解.通过坐标变换，可以得到卫星在太阳观测坐标系下的表示，即(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)为卫星处于太阳照射下的总时间.

显然，LNF为SNF的统计性结果，并可应用一些平均化方法以简化算法.对于非摄动圆轨道，ks(φ)的周期性和圆轨道的均一性表明：LNF可以理解为SNF在[0,TOrbit]范围内的时间平均，即

(34)

但摄动和偏心率的存在将使得ρ十分复杂，一般需要长时间(5～10年或更长)数值积分.

3.3 长期受晒因子的遍历算法

根据第2.3节的分析，在J2摄动的影响下，不变测度可构造为

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

(35)

为了得到LNF，定义实值函数f为LR和SR上的特征函数，即

(36)

相流定义为卫星在太阳观察坐标系下的运动轨迹，相空间取为R3.

根据Birkhoff-Khinchin定理，长期受晒因子ρ可表示为

(37)

常数κ基于μ(M)=1的事实获得，即

(38)

3.3.1 无摄圆轨道受晒因子的遍历算法

对于无摄圆轨道，仅需φ和u就可以判断卫星是否处于SR，因此不变测度将由4维退化为2维，相应的不变测度的微元为

dμ=κdu∧dφ

(39)

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

κ=1/(2π)2

(40)

则长期受晒因子ρ为

(41)

显然，该结果与式(34)一致，这也验证了前节时间平均方法的有效性.

3.3.2 受摄圆轨道受晒因子的遍历算法

对于受摄圆轨道，升交点赤径Ω不再是常值；因此，判断卫星是否位于LR区域不仅需要φ和u的值还需要Ω的值.公式(31)仍然有效，不过需要考虑独立变量Ω的影响，即改进为ks(φ,Ω).

根据式(38)，可得κ的值为

κ=1/(2π)3

(42)

长期受晒因子ρ为

(43)

3.3.3 受摄椭圆轨道受晒因子的遍历算法

对于受摄椭圆轨道，轨道半径不再是常值，判断卫星是否位于LR区域将需要φ、u、Ω和r.

前节已经详细叙述椭圆轨道卫星是否处于LR区域的判断方法.这里记

(44)

变量u和r有内在的关联性，故φ、u、Ω和r不具有遍历性；而变量φ、Ma、Ω和近地点幅角ω具有遍历性，幸运的是，两组变量存在如下变换：

(45)

(46)

3.3.4 无摄椭圆轨道受晒因子的遍历算法

无摄条件下Ω和保持不变，仅需Ma和φ即可判断卫星的受晒状态，则不变测度μ将退化为

dμ=κ·dMα∧dφ

(47)

长期受晒因子ρ为

(48)

4 结 论

本文系统地研究遍历理论在航天动力学领域的应用，包括椭圆轨道长期覆盖和受晒因子等长期参数的平均化计算.从无摄和J2项摄动下平均轨道根数的长期和长周期项出发，构造非共振条件下覆盖和受晒因子的测度描述，并证明该测度的不变性；给出相关相空间和相流定义，并证明该相流的保测度性；应用Birkhoff-Khinchin定理将长期覆盖和受晒因子等无限时间的累积问题成功转换为有限空间的二重积分运算.

研究结果表明，遍历理论是解决航天动力学领域长期参数计算的有力工具；借助动力系统等相关数学理论进行严谨化处理，可验证或改进相关经验公式和方法，并提高经典轨道理论中覆盖率与受晒因子的计算效率.