径向共线多星库仑编队飞行构型保持研究*

孙云龙,袁长清,李政广

0 引 言

近年来，利用库仑力对卫星进行编队成为一种新兴的编队技术.KING等[1]首次提出卫星库仑力编队的概念.通过卫星上的离子喷射器喷射离子，仅消耗数瓦特的功率使卫星表面带电，利用卫星间的静电引力或斥力形成控制力来对编队实现控制，以实现预期目标.与传统推进器的基于动量守恒依靠喷射气体或离子对卫星实现控制相比，库仑力控制不会产生腐蚀性的羽状物，不会影响临近卫星上的敏感元器件，直接延长了卫星的在轨运行寿命.另外对于近距离编队，库仑控制力的精度可以达到毫牛的数量级，恰满足编队要求.

文献[2]研究了二星库仑编队在地月共线平动点和三角平动点处轨道径向动力学与轨道径向线性动态分析.文献[3]采用Lyapunov反馈控制法探究了考虑太阳摄动时二星库仑编队圆轨道和地月共线平动点处轨道径向稳定性问题.文献[4]研究了编队反馈控制策略以实现虚拟结构控制.此控制策略基于推进器能力来实现对卫星三维运动的控制.文献[5]分析了旋转二星编队的稳定性，结果表明：徳拜长度大于分离距离时，非线性旋转运动局部稳定，否则运动不稳定，在忽略徳拜效应时，受扰动的面外运动总是稳定的.文献[6]研究了针对平动点处二体库仑绳系虚拟结构重构的非线性动力学和闭环控制，该控制策略对针对执行定位任务和鲁棒重构任务的卫星编队也是有效的.文献[7]设计了针对一维受限库仑编队的二阶电荷反馈控制策略，并通过探究系统总能量分析了库仑编队相对对称运动条件.文献[8]通过研究得出共线三星旋转的解族，设计出来基于线性模型的反馈控制以实现共线虚拟绳系系统的稳定.若系统角动量等于预估的角动量，则系统区域渐进稳定，据此可计算出名义(所需)电荷.文献[9]探究了共线三星编队的不变构型解并证明出：对于任何共线编队，存在无限组开环电荷解使编队可保持平衡.文献[10]还探究了共线三星库仑编队的解族，分析了不变构型解存在的区域，以确定轨迹可能的范围.考虑到编队构型维持和机动时固有的不稳定性的挑战，文献[11]证明了对于径向相连构型，面内扰动渐进稳定.文献[12]等针对面外运动进行了推广分析，并设计出简单控制律，成功消除了面外扰动的变化范围.在文献[13]中，将传统电推进技术与库仑力结合设计出了混合推进策略以适用于近距离编队飞行.文献[14]给出了针对静态四星编队电荷解的分析工具.Lee等[15]设计了基于Lyapunov法的电荷反馈率来控制二星编队的分离距离和分离率.研究表明径向三星编队也有这些性质，电荷反馈控制率可维持编队面内运动.

本文以地球同步轨道为卫星编队的运行环境，建模时考虑德拜效应，首先建立了共线N星库仑力编队沿径向方向的相对运动动力学模型.然后以四星编队为验证算例，针对变量对非线性方程进行线性化，根据线性化方程设计出LQR控制律并进行数值仿真.考虑到建模误差，添加积分项设计出改进型的LQR控制律，并进行了数值仿真验证.

1 动力学建模

1.1 坐标系的建立

如图所示建立Hill(希尔)轨道坐标系H:{Or,Oθ,Oh}，Or指向背离地心方向，Oθ指向卫星编队质心的速度方向，Oh指向卫星轨道角速度方向，Oh与Or,Oθ构成右手系.则第i颗卫星在Hill轨道坐标系中相对于编队质心的位置矢量为ρi，也可表示为ρi=[xiyizi].

编队中第i颗卫星在地心惯性坐标系中的位置矢量表示为ri，质心在惯性坐标中的位置矢量为rc， 第i颗卫星在惯性坐标系下的位置矢量可表示为：

(1)

1.2 动能与势能

根据速度合成定律，惯性系下第i颗卫星的绝对速度表示为：

vi=vρ+vo

(2)

带入式(2)得第i颗卫星绝对速度为：

(3)

设第i颗卫星的质量为mi，则编队系统动能为：

(4)

本文中，径向方向上编队构型欲保持稳定，因此沿径向方向编队可视为刚体结构.如图2、图3所示建立本体坐标系B:{b1,b2,b3}.其中b1方向沿ρn1方向，b1与b2b3构成右手系.当编队与理想构型无速度和位置偏差时，B系与Hill轨道坐标系H完全重合.

编队质心条件为：

m1ρ1+m2ρ2+…+mnρn=0

(5)

由图2可知：

(6)

联立式(5)、(6)可得，在坐标系B下，ρ1、ρ2、…、ρn可表示为：

(7)

欧拉角(θ,ψ,φ)分别为俯仰角、偏航角和滚转角，分别表示坐标系B与希尔坐标系H之间的角度偏差.本文中不考虑卫星绕径向轴的旋转，此时滚转角φ的值为零.采用3-2-1(ψ-θ-φ)的欧拉角旋转(如图3所示).此时，坐标系B与Hill坐标系H之间的坐标转换矩阵为：

(8)

Hill坐标系H中，卫星的位置向量Hρi可表示成：

(9)

将式(7)、(8)带入(9)可以得到，在希尔系H中各卫星的位置矢量为：

(10)

从上式可以看出，各位置矢量是关于各卫星间间距li和俯仰角θ、偏航角ψ的函数.对各个位置矢量求导可以得到各卫星的速度矢量分别为：

(11)

其中：

将式(10)、(11)代入到式(4)中，可以得到径向共线N星库仑编队的系统动能T为：

(12)

其中：

对于N星编队系统，其重力势能可以表示为：

(13)

(14)

其中,u为rc的单位矢量，ti为ρi的单位矢量.

将式(10)代入式(14)中，得到重力势能Vg表达式为：

(15)

卫星在太空中运行时，在受到太阳光压、地球非球形摄动等因素的影响时，还受到太空中离子浓度的影响.处于太空的带电体受其周围带相反电荷的离子的屏蔽的影响，其电势随距离长度的增加呈现出指数衰减，此长度用徳拜长度λd表示，这种现象称作徳拜效应.徳拜效应对电势的影响可表示为：

式中q为带电量，L为距离带电体的距离，λd为徳拜长度.

太空中等离子密度在不同轨道高度处不同，徳拜长度也因而不同.近地轨道处，徳拜长度范围为2～40 cm，仅适用于卫星机构的对接操作.同步轨道处徳拜长度范围为140～1 400 m，恰适用于星间间距在几十到数百米的卫星编队.深空中其范围为8～24 m，也可运用于更近距离的编队飞行.

对于N星库仑力编队系统，每颗卫星与其他不相邻的卫星间都具有库仑力相互作用.故其库仑力势能可表示为：

(16)

其中kc为静电引力常数，qi和qj分别为第i颗卫星与第j颗卫星的电荷量，λd为德拜长度.

1.3 N星编队相对运动方程

设系统拉格朗日函数为[5]：

L=T-(Vg+Ve)

(17)

拉格朗日方程为：

(18)

q=(θ,ψ,l1,l2…ln-1)

其中Qi为第i颗星所受到不包括地球重力和库仑力的广义力.

将式(12)、(15)、(16)代入式(18)得共线N星库仑编队相对运动方程为：

Ω2(3cos2θcos2ψ-1)]-

Ω2(3cos2θcos2ψ-1)]-

……

(19)

此处Qij为卫星间电荷乘积 ，即Qij=qiqj.上述的N+1式具有很强的非线性，且卫星之间的距离li(i=1,2…n-1)与俯仰角θ和偏航角ψ高度耦合，这对于控制器的设计及仿真校验都是不小的考验.

本文选取地球同步轨道为库伦编队运行环境.在这一轨道高度处，徳拜效应对于编队中各卫星间库仑力的作用比较有限.故下文在进行建模时将略去徳拜效应的影响.对于多星编队的建模及仿真，仍不失一般性.

1.4 四星编队相对运动方程

上文求得了共线N星库仑编队的相对运动方程组.由于其中各方程具有很强的非线性，各方程间相应变量又高度耦合，因此完全对N星库仑编队进行仿真校验是不实际的.不失一般性地，为了验证动力学方程的正确，进行有效的数值模拟，下文将进行针对四星编队的动力学模型建立及数值仿真.

令a1=m1(m2+m3+m4)，a2=(m1+m2)(m3+m4)，a3=(m1+m2+m3)m4，a4=m1(m3+m4)，a5=m1m4，a6=(m1+m2)m4.

(20)

共线四星库仑编队在径向方向相连，当该编队不受到其他摄动因素的影响而达到理想构型时(编队相对于Hill轨道坐标系保持静态稳定)，此时相对运动方程满足以下平衡条件：

l1=l1f,l2=l2f,l3=l3f

(21)

式中l1f为星1与星2之间的期望距离，l2f为星2与卫星3之间的期望距离，l3f为星3与星4之间的期望距离.将上式分别带入到运动方程，可得到保持共线四星库仑静态编队平衡所需的参考电荷积为：

(22)

Q12f、Q13f、Q14f、Q23f、Q24f、Q34f为期望距离所对应的参考电荷积.对于所期望的共线四星库仑编队来说，为使库仑力作为控制力，以保持卫星编队的静态构型达到期望，各卫星的带电量需要满足以上各式.

对于俯仰和偏航运动方程分别以小俯仰角和小偏航角进行线性化处理；式中的卫星间距离l1、l2、l3以间距误差δl1、δl2、δl3进行线性化处理，电荷积Q12、Q13、Q14、Q23、Q24、Q34以电荷积误差δQ12、δQ13、δQ14、δQ23、δQ24、δQ34进行线性化处理，如下所示：

l1=l1f+δl1,l2=l2f+δl2,l3=l3f+δl3,

Q12=Q12f+δQ12,Q13=Q13f+δQ13,

Q14=Q14f+δQ14,Q23=Q23f+δQ23,

Q24=Q24f+δQ24,Q34=Q34f+δQ34

(23)

利用式(23)可得到线性化动力学方程：

3Ω2(a1δl1+a4δl2+a5δl3)-

3Ω2(a2δl2+a4δl1+a6δl3)-

3Ω2(a3δl3+a5δl1+a6δl2)-

(24)

(25)

式(25)中前两式为线性姿态动力学方程，δL与ψ相互耦合.后3个式子为编队关于径向理想构型间距误差的线性微分方程.可以通过控制电荷使得δL、ψ方程渐近稳定.

2 控制器设计

本节将基于线性化相对运动动力学模型，仅利用库仑力作为控制力，设计LQR控制律使径向共线四星库仑编队保持稳定构型.考虑到模型误差，设计了改进型的LQR控制律，以提高传统LQR控制器在控制编队卫星构型时的鲁棒性.

2.1 传统LQR控制器

卫星编队的俯仰角θ不与偏航角ψ和间距误差δl1、δl2、δl3耦合，不受库仑力控制，所以不考虑俯仰运动.下面仅利用卫星间相互作用的库仑力对编队的面内运动(即间距误差δl1、δl2、δl3和偏航角ψ)进行控制[17].首先，定义状态向量x为：

控制器期望实现的目标为x→0.

将线性化的共线四星库仑编队构型相对运动动力学模型表示为状态空间形式为：

(26)

其中矩阵A1为8×8的状态矩阵，矩阵B1为8×6控制矩阵，矩阵C1为4×8矩阵.

定义输入量：u=[δQ12δQ23δQ34δQ13δQ24δQ14]T,通过设计输入量u,使系统满足性能指标：

(27)

其中Q1和R1分别为状态变量和输入变量的加权矩阵.

当控制输入为：

(28)

性能指标泛函数最小.式中矩阵P满足Riccati微分方程，即：PA+ATP-PBR-1BTP+Q=0.

2.2 改进LQR控制器

针对地球同步轨道处四星库仑力卫星编队运行时力学环境复杂、非线性强、控制变量耦合性 强、传统LQR控制器鲁棒性不足等现象，在传统LQR控制器的基础上，考虑到未建模摄动力的影响及线性化产生的误差，设计出一种改进型的LQR控制器，以提高控制的鲁棒性.

考虑如下式所示的带有不确定项的动力学系统的状态空间表达式：

(29)

式中ω为未建模摄动力.令ω=Dv，v为解决不确定项的辅助扩展控制向量，取D=I-BB+.

为使系统在不确定扰动ω(x)的作用下渐进稳定，可通过寻找反馈控制律[uv]T,使系统满足下列目标泛函数.

uTHu+ρ2‖v‖2)dt

(30)

式中F和H均为对角矩阵；α≥0，β≥0，ρ≥0;且2ρ2‖v‖2≤β2‖x‖.

当α=ρ=0，β=1，且F=Q1，H=R1时，即转换为式(27)所示的传统LQR控制.

为了使性能指标泛函数J1最小，取控制输入为：

(31)

式中矩阵P满足Riccati微分方程，且

3 数值仿真

3.1 初始条件

本节针对径向共线四星库仑卫星编队飞行队形保持问题作为仿真算例.在Matlab/Simulink环境下进行数值仿真，对共线四星库仑力相对运动方程和控制算法的有效性进行检验.航天器参数取值[16]为：m1=m2=m3=m4=50 kg,l1f=l2f=l3f=20 m,仿真参数的取值为：kc=8.99×109N·m2/C2,Ω=2.66×10-6rad/s,δψ(0)=0.1 rad，δl1(0)=δl2(0)=δl3(0)=0.6 m，q1=1.869 7 mC,q2=42.445 mC，q3=-46.569 mC，q4=-195.103 mC.

传统LQR控制中，状态变量的加权矩阵Q1为：

Q1=diag{0.001,0.1,0.001,0.01,0.01,0.001,0.1,0.001}

输入变量的加权矩阵R1为：

R1=diag{108,108,108,107,107,106}

在改进LQR的鲁棒控制中，取α=ρ=β=1.状态变量的加权矩阵Q2为：

Q2=diag{1,10,30,0.5,0.01,0.01,0.1,0.01,0.01,0.1,0.001,0.1}

输入变量的加权矩阵R2为：

R2=diag{107,108,107,107,107,107}

3.2 仿真结果及分析

在MATLAB中进行数值模拟.仿真结果如图4所示.

图4为LQR控制策略下有无积分项姿态角ψ随轨道时间的变化情况.可以看出，与无积分项相比，偏航角误差的仿真趋势与其基本相同，收敛略早.图5、图6及图7分别表示LQR控制策略下有无积分项卫星间距l1、l2和l3的变化情况.可以看出，加入积分项的控制有效的减少了间距收敛至稳定的时间，减小了间距变化的振荡幅度.在相同轨道时间时，带积分项控制策略下间隔距离误差优于无积分项.由于积分项的补偿作用，在控制稳定阶段的稳定误差大大减少.

4 结 论

本文研究了径向共线多体库仑卫星编队在地球同步轨道处轨道动力学与控制问题.通过拉格朗日函数推导出了共线N星库仑编队在同步轨道处的一般动力学方程.以四星库仑编队作为验证算例，针对共线四星卫星编队设计了LQR控制律和加入积分项的LQR控制律.通过Matlab/Simulink数值仿真分析(仿真过程中忽略德拜效应)，验证了改进LQR控制律的有效性.仿真结果表明，改进型LQR控制律能够更快达到期望构型，该控制可行、有效.基于本文研究，可进一步研究多星库仑编队在同步轨道处的切向或法向轨道动力学.