航天器相对运动建模及周期性相对运动求解

刘柔妮，陈 杰，孔祥龙，周世宏

(1.上海卫星工程研究所，上海 201109； 2.上海航天技术研究院，上海 201109)

0 引言

航天器相对运动建模是编队飞行等航天领域的研究热点，研究空间自然条件下周期性相对绕飞轨道对长期编队飞行航天器有十分重要的意义。通过找到不消耗能量的自然绕飞轨道，可节省编队航天器系统的能量损耗，延长系统寿命。文献[1]利用Hill方程和Clohessy-Willshire方程(简称HCW方程)进行编队卫星构型设计，在圆参考轨道、无摄动和线性化的假设条件下，找到了周期性相对运动轨道存在的初始条件。LAWDEN[2]研究了椭圆参考轨道的相对运动。TCHAUNER等[3]以真近点角为独立变量，通过限制参考基准轨道偏心率避免了Lawden方程存在的奇点问题。TILLERSON等[4]研究了椭圆轨道中编队飞行卫星的控制问题。ALFRIEND等[5]率先采用相对轨道要素描述相对运动。肖业伦等[6]也采用相对轨道要素描述相对运动。之后，韩潮等[7-8]在该概念的基础上进行了一定的研究工作。但以上研究成果均存在不足：Hill方程和Clohessy-Willshire方程使用条件苛刻，需要圆参考轨道、线性化的假设前提，工程应用受到限制；真近点角、相对轨道要素等描述方式较难考虑周期性相对运动的初始条件，求解难度较大。

本文采用哈密尔顿-雅可比(HJ)方程和正则摄动理论，引入适当的正则变量描述相对轨道运动，建立适用于椭圆参考轨道，便于考虑各种摄动的相对运动模型。经推导发现，在各种摄动下相对运动的哈密尔顿方程可分解为线性项和摄动项，而相对运动在任意摄动情况下的解可通过求解包含正则变量的微分方程和修正线性项得到，且模型满足摄动叠加原理，即如已得到一定摄动下的相对运动状态，则其他摄动下的相对运动状态只需通过摄动项的叠加即可得到。得出周期性相对运动条件后，用时域配点法和列文伯格-马夸尔特法(LM)对该条件的初值求解，可得到不消耗任何燃料的周期性绕飞轨道。

1 椭圆参考轨道相对运动建模

图1 旋转Euler-Hill参考坐标系示意图Fig.1 Schematic diagram of rotating Euler-Hill reference coordinate system

建立拉格朗日函数[9]，设主航天器在任意椭圆轨道上运动，考虑活动坐标系之间的相对导数关系，可得从航天器的速度

(1)

(2)

动能表示为

(3)

势能Ep(先只考虑正常引力势)为

Ep= -μ/‖r2‖=-μ/‖r1+ρ‖=

(4)

式中：μ为地球引力势；Pk(cosα)为勒让德多项式,其中,

(5)

从而得拉格朗日函数

(6)

只取正常引力势中的第1，2阶可得

(7)

哈密尔顿系统的正则动量(px,py,pz)为

(8)

(9)

由式(9)可见，哈密尔顿函数不显含时间t，则根据哈密尔顿动力学原理可知，HJ方程有如下形式，即

(10)

考虑由(q,p)→(β,α)的正则变换，其中，q,p,β,α为正则变量。S为雅可比完全积分，其形式一般为

S=W(q1,q2,…,qf;α1,α2,…,αf)-

E(α1,α2,…,αf)t

(11)

式中：W(q1,q2,…,qf;α1,α2,…,αf)为哈密顿特征函数，(q1,q2,…,qf)为广义坐标；E(α1,α2,…,αf)为不含时间t的函数，(α1,α2,…,αf)为积分常数；t为时间；f为系统的自由度。则正则变换表示为

(12)

系统自由度f现为3，令式(11)中E(α1,α2,…,αf)=-α′1，则可得雅可比完全积分S为

(13)

从而可将HJ方程表示为W(x,y,z)的偏微分方程，即

(14)

将z方向单独分离出来，W(x,y,z)=W′(x,y)+W3(z)，则HJ方程变为

(15)

(16)

将式(15)积分得

(17)

再令

(18)

同时令

(19)

由式(19)可得

(20)

(21)

积分得

W1(x)=

(22)

因此，可将雅可比完全积分S表示为

(23)

(24)

(25)

(26)

再根据正则变换式(12)可得βi，经过一系列化简和小量忽略可得表达式

(27)

(28)

再联立式(24)，(26)，(27)，(28)，解出x(t)，y(t)关于α1，β1，α3，β3的表达式，经化简后可得

(29)

(30)

由此可得x，y方向上的运动形式与正则变量α1，β1，α3，β3的关系。如已知α1，β1，α3，β3的函数表达式，就可得到x(t)和y(t)的表达式，从而可知从航天器相对于主航天器在x-y平面内的运动模式。

由上文推导过程可知，z方向的运动和x，y方向是完全解耦的。可先研究x-y平面内的运动特性，再通过分别控制x-y平面的运动和z方向的运动得到三维空间的运动模式。因此，本文主要针对x-y平面内的运动特性进行分析，未专门对z方向上的运动进行研究。

以上推导都是基于正常引力势的第1，2阶项得出，将此时得到的哈密尔顿函数(9)记为H(0)，当考虑引力势的更高阶项或J2项时，将哈密尔顿函数记作H=H(0)+H(1)+H(2)+…+H(n)。其中，H(1)，H(2)，…，H(n)是考虑摄动或非线性引力势附加的哈密尔顿函数ΔH，其被称为摄动哈密尔顿函数，即ΔH=H(1)+H(2)+…+H(n)。

当考虑摄动时，由哈密尔顿理论可能较难得出相对运动的显示解析表达式，但无摄情况下的运动模型已精确得到，即式(29)，(30)。因此，根据正则摄动理论，当ΔH很小时，可对无摄条件下得到的动力学模型作适当修正，从而得到摄动下动力学模型的近似解。

将不考虑摄动情况下精确的哈密尔顿函数记作H0(p,q,t)，则可将雅可比完全积分S(q,α,t)作为广义函数，得到从(q,p)→(β,α)的正则变换，且变换后的新的哈密尔顿函数为零，则(β,α)是一组与时间无关的正则变量，由此可以得出相应的在无摄情况下关于模型的精确解。当考虑摄动时，将哈密尔顿函数写成

H(p,q,t)=H0(p,q,t)+ΔH(p,q,t)

(31)

考虑到对于给定坐标变换的正则性与哈密尔顿函数的具体形式无关，因此，由S(q,α,t)作为广义函数而产生的正则变换(q,p)→(β,α)在摄动条件下还是一组正则变换，只是在摄动下得到的新哈密尔顿函数将不再为零，新的正则变量(β,α)也不再是与时间无关的变量。新哈密尔顿函数的形式为

(32)

因此，正则变量满足关系式

(33)

(34)

求解上述2f个方程(f为自由度，本文有4个自由度)，得到αi，βi关于时间的函数表达式，则由(q,p)→(β,α)的变换，可得相应的qi，pi关于时间的表达式，从而得到动力学模型。

将上述理论运用到考虑非线性引力势第3，4阶，以及J2项的建模中。非线性引力势第3阶的摄动哈密尔顿函数可表示为

(35)

非线性引力势第4阶的摄动哈密尔顿函数可表示为

24μx2y2+24μx2z2-6μy2z2)

(36)

考虑J2摄动项的摄动哈密尔顿函数为

(37)

式中：Z=(r+x)sinisinθ+ysinicosθ+zcosi；引力常数μ=3.960 05×105km3/s2；地球赤道平均半径Re=6 378.140 km；带谐系数J2=0.001 082 63；i为参考轨道倾角；θ为真近点角。

将上述摄动哈密尔顿函数代入求解后可得αi(t)，βi(t)关于时间t的表达式，再代入(q,p)→(β,α)变换关系式中，即式(29)，(30)中，可得x(t)和y(t)在摄动情况下的表达式，即摄动情况下的模型为

(38)

y(t)=-3α3(t)(t+β1(t))+β3(t)+

(39)

(40)

综上，求解摄动情况下的动力学模型只需求解一组非线性方程组，即

(41)

(42)

则可得动力学方程(38)，(39)。

2 周期性相对运动条件

椭圆基准轨道的动力学模型为式(38)，(39)，经分析不难发现，获得周期为T的周期性相对运动的条件为：∀t满足关系式

α1(t)=α1(t+T)

(43)

α3(t)=α3(t+T)

(44)

f1(t)=β3(t)-3α3(t)(β1(t)+t),

f1(t)=f1(t+T)

(45)

即

α1(t)=α1(t+T)

(46)

α3(t)=α3(t+T)

(47)

β1(t)=β1(t+T)

(48)

(49)

3 周期性相对运动初始化条件求解

将周期条件截断为含有N项的傅里叶级数的形式，即

(50)

(51)

(52)

(53)

(54)

则周期性相对运动的初值条件为

(55)

(56)

(57)

(58)

时域配点法是一种非线性系统的求解方法。它将常微分方程组转化为非线性代数方程组，再用其他求解非线性代数方程组的方法对得到的方程组进行求解，避免了直接求解常微分方程组，对于周期解的处理十分有效。该方法在文献[10-11]中被提出，本文不对此作详细说明。

(59)

利用时域配点法进行相应变换，式(59)左端有

(60)

因此，微分方程组转化后的非线性代数方程组表示为

(61)

求解该非线性代数方程组，可得4×(2N+1)个配点的值，有

(62)

(63)

由此可求得4×(2N+1)个傅里叶系数，将其代入式(55)，(56)，(57)，(58)可得周期运动初值。

采用LM方法对式(61)的非线性代数方程组

(64)

(65)

进行求解，可避免采用牛顿-拉斐森迭代法(NR)进行求解时由于方程形式过于复杂，雅可比矩阵出现奇异或病态严重，无法继续迭代的情况。

图2 参考椭圆轨道偏心率为0.02仿真结果Fig.2 Simulation results of elliptical orbit with eccentricity of 0.02

采用一种全局收敛的LM算法[12-15]，考虑矢量形式的非线性方程组F(x)=0。利用线性方程组的解dk作为在点xk上的一个搜索方向：((Jk)T·Jk+μkI)dk=-(Jk)TFk。其中，I为单位矩阵，Jk为雅可比矩阵。通过引进非负参数μk可保证矩阵((Jk)TJk+μkI)非奇异，同时也能避免出现过大的‖dk‖。同时，对于迭代步((Jk)TJk+μkI)dk=-(Jk)TFk，选取参数μk=αk(θ‖Fk‖+(1-θ)‖(Jk)TFk‖),θ∈[0,1]。其中，对于参数αk的选取采用信赖域技巧来修正。定义价值函数φ(x)=‖F(x)‖2，第k步迭代的实际下降量Ak=‖Fk‖2-‖F(xk+dk)‖2,预估下降量Pk=‖Fk‖2-‖Fk+Jkdk‖2，实际下降量与预估下降量的比值rk=Ak/Pk，由此来决定是否接受试探步dk和调整迭代过程中αk因子的大小。rk越大，表明目标函数φ(x)=‖F(x)‖2下降越多，则接受dk，同时希望下一试探步dk+1更长，故减小αk；反之，rk越小，则拒绝dk，同时增大αk。

4 仿真算例

对考虑正常引力势的非线性3，4阶项进行仿真，取椭圆参考轨道长半轴a=8 000 km，过近地点时刻tp=0，偏心率e=0.02，其余3个轨道根数均取为0，将傅氏级数截断为N=3，在一个周期上均匀地取7个配点，仿真6个轨道周期，结果如图2所示。

为定量描述每个周期的轨道漂移量的大小，以衡量周期运动的好坏，引入如下指标来粗略描述经过一个周期T的漂移量，即

(66)

经计算可得ΔT=0.315 9 m。由此可见：每过一圈的漂移量是很小的。这充分说明周期条件初值在用于轨道计算时能保持较好的周期性，在轨道保持和控制等方面具有工程应用价值。仿真结果验证了模型和求解方法的有效性。

5 结束语

本文主要基于哈密尔顿力学建立椭圆参考轨道的相对运动模型，利用正则摄动理论将哈密尔顿函数分解为线性项和摄动项。线性项的精确解很容易得到，摄动项则可通过线性项的修正获得，且各摄动项满足线性叠加的特点，各种轨道摄动和人为控制易于加到运动模型中。同时，针对模型特点推导了周期性相对运动条件，利用时域配点法和LM方法对周期性相对运动的初始值进行求解，得到了不消耗任何燃料的周期性绕飞轨道，为空间交会、编队飞行、机动伴飞等航天领域的研究提供了有益参考。