CN106697333A - 一种航天器轨道控制策略的鲁棒性分析方法 - Google Patents

Abstract

本发明涉及一种轨道控制策略的鲁棒性分析方法，属于航天器轨道动力学与控制领域。本发明通过高斯轨道要素摄动方程对航天器的轨道运动进行建模；对低轨卫星的非球形摄动以及大气阻力摄动进行分析；设计轨道保持策略对航天器轨道要素进行保持控制；采用微分修正算法提高轨道保持的精度；建立航天器运行过程中的位置误差，速度误差，发动机推力误差模型；设计了航天器控制误差均值，方差，以及误差分布比例的计算模型；利用蒙特卡洛仿真方法对带有误差下的航天器轨道控制策略进行仿真分析，建立轨道控制策略的鲁棒性评价体系。在控制策略设计中，充分考虑星上轨道控制的实际情况，以简便可行为前提，贴合实际情况，保证该方法在实际工程中的可行性。

Description

一种航天器轨道控制策略的鲁棒性分析方法

技术领域

本发明涉及一种航天器轨道控制策略的鲁棒性分析方法，尤其涉及航天器轨道控制方法和随机误差存在下的控制误差评价体系，属于航天器轨道动力学与控制领域。

背景技术

航天器在理想的无摄动环境下能够遵守Kepler轨道的运动规律，进行长期稳定的在轨运行。然而现实环境中不可避免的存在各种轨道摄动因素，例如地球扁率摄动、大气摄动、太阳光压摄动、日月摄动等，导致航天器偏离预定轨道，从而对航天器执行相关任务，尤其是精密任务，带来偏差和不便。因此对长期在轨的航天器进行轨道保持控制研究是一项必要的工作，同时也是航天器执行其他任务的基础。

轨道维持策略有大推力脉冲式控制方案，小推力连续控制方案等。相对于小推力连续控制方案，脉冲控制方案具有很多优点，如工程上更容易实现，可靠性高，计算效率也更高。因此本发明选择脉冲式控制方法进行轨道控制。

在理想轨控过程中，首先需要对当前时刻航天器所处的位置进行测量，并与标称轨道进行对比，从而确定轨道参数机动量的大小，然后启动轨控发动机，产生脉冲对轨道进行修正，从而达到轨道保持的目的。但是在轨道确定和发动机推力执行上都不可避免的存在误差，这直接影响了轨道控制的性能，甚至决定着任务的成败。

鲁棒性分析是设计轨道维持策略的一个重要环节。通过鲁棒性分析，轨道控制策略的实用性和可靠性得到验证，为控制策略的在轨使用提供了依据和保障。而目前的鲁棒性分析则仅仅只针对一种控制方法或是某个特定的误差。文献(Howell,K.C.,and Barden,B.T.,“Trajectory design and station-keeping for multiple spacecraft information near the Sun-Earth L1point,”IAF 50th International AstronauticalCongress,1999)研究了一种平动点附近航天器编队飞行的轨道保持策略，并采用蒙特卡洛仿真验证了其具有良好的鲁棒性。文献(Rozanov,M.,and Guelman,M.,“Aeroassistedorbital maneuvering with variable structure control,”Acta Astronautica,Vol.62,No.1,2008,pp.9-17)研究了一种变结构闭环控制方法用于应对大气捕获时大气密度的不确定性问题，并进行了蒙特卡洛仿真以评价该方法的鲁棒性。目前的轨道控制策略鲁棒性分析具有一定的针对性，并没有从系统的角度建立评价卫星轨道控制策略的鲁棒性评价体系，因此本发明旨在提出一种航天器轨道控制策略的鲁棒性的分析方法，建立一种航天器轨控策略的鲁棒性评价体系，从而为轨控策略的实际应用提供依据和保障。

发明内容

本发明的目的是提供一种航天器轨道控制策略的鲁棒性分析方法，该方法普遍适用于任何轨道控制策略，所考虑的误差种类也不受限制，能对任一轨控策略的鲁棒性进行分析，并得出可靠结论，为该策略的实际应用提供依据和保障。

本发明的目的通过以下技术方案实现。

本发明公开的一种航天器轨道控制策略的鲁棒性分析方法。本发明通过Gauss轨道要素摄动方程以及平均轨道要素对航天器的轨道运动进行建模；对低轨卫星的非球形摄动以及大气阻力摄动进行分析；轨道保持策略对航天器轨道要素进行保持控制；采用微分修正算法提高轨道保持的精度；建立航天器运行过程中的随机误差模型；利用蒙特卡洛仿真方法对带有误差的航天器轨道控制策略进行仿真；设计航天器控制误差均值，标准差，以及误差分布比例的计算模型，并对数据结构进行统计分析，建立轨道控制策略的鲁棒性评价体系，用于分析轨道控制策略的鲁棒性。

一种航天器轨道控制策略的鲁棒性分析方法，包括如下步骤：

步骤一：通过Gauss轨道要素摄动方程和平均轨道要素建立航天器轨道运动模型。

采用Gauss轨道要素摄动法建立的航天器六个轨道要素半场轴a，偏心率e，倾角i，近地点角距ω，升交点赤经Ω和真近点角θ的轨道运动模型如下所示：

其中，

f_r，f_u，f_h是定义在轨道坐标系S_o的卫星摄动力在径向，横向(即在轨道平面内垂直于矢经，沿纬度幅角u增大的方向)，和副法向分量(即垂直于轨道平面，与动量矩H同向)。p是半通径，μ是地球引力常数，r是地心与航天器之间的距离。

航天器轨道要素平均化方法如下：

公式(1)的向量表达形式如下：

其中

σ＝(a,e,i,Ω,ω,M)^T (3)

这里ε＝O(J₂)即通常所取的小参数。等式右端的f₁和f₂为非球形引力加速度的部分，它们分别对应J₂项和其他非球形引力项，并有

f₁/f₀＝O(ε),f₂/f₀＝O(ε²) (5)

公式(2)是一个非常典型的小参数方程。本步骤中采用古在由秀提出的根据非线性力学中平均法的思想，针对主要带谐项摄动的平均根数法。平均根数法将摄动变化项Δσ⁽¹⁾(t),Δσ⁽²⁾(t)，……中不同性质的项(长期项，长周期项和短周期项)区分开，更重要的是在此前提下将参考解取为只包含长期变化项的平均根数。

对Gauss轨道要素摄动方程中积分得出的瞬时轨道要素，利用古在优秀方法计算出其平均轨道要素，从而用于后续的控制和鲁棒性分析。

步骤二：对低轨卫星的非球形摄动以及大气阻力摄动进行分析

步骤2.1：地球非球形摄动；

地球非球形摄动中，由于地球的旋转，在运行时间超过一天时，扇形谐函数和田形谐函数的影响或多或少地被“平均掉”，因此对于大多数应用情况，只需考虑带形谐函数的影响。其形式如下：

其中，J_n是n阶带形谐函数的系数，P_n和P′_n是n阶勒让德多项式及其导数，x，y，z是惯性坐标系S_i下的航天器三维坐标，且R_E是地球的平均赤道半径。因此由式(6)得到的惯性系下的非球形摄动需要经过坐标转换，以符合式(1)中对于摄动力的定义环境。

由惯性系到轨道系的坐标转换分为三步：

其中，u是纬度幅角，u＝ω+θ。

转换矩阵具体形式如下：

因此，满足公式(1)中定义的摄动力的地球非球形摄动加速度在轨道系的分量形式如下：

步骤3.2：大气阻力摄动；

由空气动力学可知，在大气中运动的航天器受到空气阻力为

其中，C_D为阻力系数，ρ为大气密度，S为横截面积，v_a是航天器质心相对于大气的速度。为简化计算，假设大气是球对称的；大气不随地球旋转；航天器的迎风面积不变。于是大气阻力只沿引起沿轨道切线方向的负加速度，该加速度写在速度坐标系中如下：

f_n，f_t，f_h分别是大气阻力摄动加速度在速度坐标系的三个分量，其中f_t沿速度方向，f_h垂直于轨道平面指向角动量方向，f_n与其余两个向量垂直，满足右手定则。将速度坐标系S_v中的大气摄动加速度转换到轨道坐标系S_o中,即该坐标系的转换只需要绕速度坐标系的Z轴进行一次欧拉转换，转换矩阵为：

因此，满足公式(1)中定义的摄动力的大气阻力摄动加速度在轨道系的分量形式如下：

将地球非球形摄动力和大气阻力摄动力带入公式(1)中的摄动力中，即

则得到了带有摄动的航天器轨道运动模型。

步骤三：使用轨道保持策略对航天器轨道进行控制。

名义轨道的轨道要素记为σ^*，由步骤二中的轨道动力学模型数值积分所得任意时刻航天器的密切轨道要素记为σ(t)，则轨道要素偏移量记为δσ(t)＝σ(t)-σ^*，记可容忍的最大轨道要素偏差区间为Δσ，轨道保持采取偏差-修正的控制策略，即当则施加脉冲进行轨道修正，其余轨道要素的控制同理可得。

步骤四：对步骤三中计算所需速度脉冲采用微分修正算法进行修正，从而完成精确的轨道机动。

步骤五：建立误差模型

本控制系统中，所考虑的误差是任意实际情况中会出现在航天器轨道控制过程中的随机误差，误差的数量不受限制。一般地，为分析航天器轨控策略的鲁棒性，卫星定轨误差，发动机执行误差是主要考虑的误差类型。定轨误差分别加在卫星的三个方向的位置和速度上，发动机执行误差则加在三个方向的发动机脉冲分量上。假设各误差均服从正态分布。

步骤六：蒙特卡洛仿真分析

采用带有微分修正的轨控策略对带有摄动的航天器轨道进行控制，并加入步骤五中的误差模型，采用MATLAB中normrnd函数对具有同一控制误差标准差的航天器轨道重复100次仿真，消除随机性，模拟体现真实情况。

步骤七、计算模型

本发明中以所有未落在控制范围中的点的均值和分布情况作为控制误差的评价方法。首先计算点的均值和标准差，因此控制误差的计算方法如下所示：

其中，

σ'＝0,Δσ∈[c,d] (15)

Δσ是轨道要素相对于其初始值的变化量。最大容许误差范围是[c,d]，其中c和d分别是控制范围的下边界和上边界。N为所有数据点的总个数。在式(15)中，当轨道要素的变化量在最大容许误差范围内时，认为控制误差为零。当轨道要素的变化量超出了控制范围，σ'为超出控制范围的控制结果点到最大容许误差边界的距离。第二个参数是控制误差标准差，它表示了控制误差的分布情况，具体计算公式如下：

第三个参数是控制误差分布比例，表示超出控制范围的点数占总数据点的比例。具体计算公式如下：

这三个参数可充分反映控制误差的各项分布。因此，可通过计算ξ，k的值，对在某种误差存在的情况下的轨道控制性能进行评价，并给出该控制策略鲁棒性的分析结果，为轨控策略的实际应用提供依据和保障。

步骤五所述的微分修正方法为：

航天器的运动方程(公式(1))可简化为

其中x是描述航天器运动的状态量。

假设航天器的初始位置状态向量为x₀，在施加轨道机动ΔV的情况下到达位置为x^t，而其目标到达位置为x^d。

约束函数F(X)表示航天器目标位置与施加机动后的位置之差，如下所示：

F(X)＝x^t-x^d (20)

当F(X)＝0，则表示当前施加的轨道机动ΔV能精确地将航天器机动到目标位置。而为了寻找合适的ΔV，采用牛顿迭代法进行实现。根据Taylor多项式的展开法则，省略后面的高阶项，可以得到如下表达式：

F(X^j)+DF(X^j)(X^j+1-X^j)＝0 (21)

其中DF(X)表示约束函数对状态变量的偏导数，即状态变量的雅克比矩阵，形式如下：

该矩阵表示了轨道要素从初态转移到末态的状态转移矩阵，具体的数值大小可通过给定初值的函数积分得到。

给定一个合理的初值，将式(22)带入式(21)中，则可通过迭代计算得到新的X^j向量，理论上

||F(X^j+1)||＜||F(X^j)|| (23)

于是，通过数次迭代计算，可使F(X)函数的值趋向于0。当

||F(X^j+1)||＜ε (24)

其中||F(X^j+1)||是F(X^j+1)的模。ε是一给定的极小量。通过上述方法，最终找到合适的机动向量，从而完成精确轨道机动。

有益效果

1、本发明公开的一种轨道控制策略的鲁棒性分析方法，建立了一套评价轨道策略鲁棒性的标准体系，该体系方法适用所有轨道策略的鲁棒性分析，其中误差的种类和数量亦不受限制。通过该方法对轨控策略进行鲁棒性分析，为航天器轨道控制的实用性和可靠性提供了理论依据。

2、本发明公开的一种轨道控制策略的鲁棒性分析方法，考虑了任意随机误差对轨道控制策略的影响。借助蒙特卡洛仿真，并设计了三个控制误差参数对庞大的轨控数据结果进行分析，得到航天器轨道控制策略在随机误差存在的情况下的轨控性能分析结果。

3、本发明公开的一种轨道控制策略的鲁棒性分析方法。可通过改变控制范围，分析对不同最大允许控制范围下的控制策略的鲁棒性，从而得到该控制策略在不同任务中的随机误差大小与控制性能之间的联系，以便同一控制策略在不同轨道控制任务中的应用。

附图说明

图1为本发明实施例1轨道控制结果，其中：图1(a)为半长轴的变化量、图1(b)为偏心率的变化量、图1(c)为近地点角距的变化量；

图2为位置误差存在的情况下该脉冲式轨道维持策略在给定仿真时长下的脉冲数；

图3为位置误差下半场轴和近地点角距的控制误差，其中图3(a)为半长轴的控制误差、图3(b)为近地点角距的控制误差；

图4为速度误差存在的情况下该脉冲式轨道维持策略在给定仿真时长下的脉冲数；

图5为速度误差下半场轴和近地点角距的控制误差，其中图5(a)为半长轴的控制误差、图5(b)为近地点角距的控制误差；

图6为发动机推力误差存在的情况下该脉冲式轨道维持策略在给定仿真时长下的脉冲数；

图7为发动机推力误差下半场轴和近地点角距的控制误差，其中图7(a)为半长轴的控制误差、图7(b)为近地点角距的控制误差。

具体实施方式

为了更好地说明本发明的目的和优点，下面结合附图对本发明的具体实施方式作进一步详细的说明。

实施例1：

本实施例公开的一种轨道控制策略的鲁棒性分析方法，为验证该方法，首先，选取高度500公里的运行在近圆太阳同步轨道上的卫星作为主要研究对象。其卫星参数如下表所示。

表1 卫星仿真参数

其中，T表示一个轨道周期。

步骤一：利用Gauss轨道要素摄动方程及平均轨道要素建立其轨道动力学的模型。

Gauss轨道要素摄动方程如公式(1)所示。所研究卫星运行在近地轨道，因此航天器所受摄动力中包含了地球非球形摄动。由于J₂项摄动的影响，航天器轨道要素会产生短周期的振荡，这对精确的轨道控制产生干扰和影响，因此为滤去短周期项，以分析航天器轨道要素的长周期变化情况，并采用轨控策略对其进行控制。采用古在由秀提出的平均根数法，对航天器轨道要素进行平均化处理。

由于本实施例背景为近圆的太阳同步轨道卫星，因此这里采用同时消除e＝0和通约奇点的摄动计算方法。为解决小e问题，引用无奇点变量，即a，i，Ω，ξ＝ecosω，η＝-esinω，λ＝M+ω。并仍然统一记为向量形式σ。根据构造摄动解的方法，利用原平均根数法的结果，很容易给出同时消除e＝0和通约奇点的摄动计算方法，对一阶解其短周期项表示为

这里是完整的一阶短周期项。其中

J₂，J₃，J₄分别是地球非球形摄动中的二阶，三阶和四阶带形谐函数系数。

采用这种方法对所求的瞬时轨道要素进行平均化处理，略去了一阶短周期项。这样做的好处主要是通过滤掉短周期项的振荡，从而更清楚地看出轨道的长期变化，更有利于准确判断是否需要施加轨道机动，避免发动机多次开机导致不必要的燃料消耗和错误的轨道机动。

步骤二：对低轨卫星的非球形摄动以及大气阻力摄动进行分析

对低轨近圆轨道，非球形摄动中的J₂，J₃，J₄项和大气阻力摄动项对航天器的影响显著，因此该卫星轨道动力学模型中主要考虑以上摄动影响。

具体的非球形摄动项则是公式(6)中求和上界变为4的表达式，如下所示：

并将公式(33)进行坐标转换，具体坐标转换矩阵为L_oi，因此则得到了非球形摄动的表达形式。大气阻力摄动如公式(12)所示，本实施例中使用的大气密度模型是美国国家标准大气模型(SA76)。

步骤三：使用轨道保持策略使航天器轨道保持在名义轨道附近的给定范围内。半长轴在名义轨道附近的范围是[-500,0]m，近地点角距是[-3,0]deg。

采用双脉冲闭环控制策略，具体如下：

航天器运行在高度为500km的太阳同步轨道，受地球非球形摄动的影响，近地点角距ω和升交点赤经Ω发生长期漂移，其中Ω的长期变化率恰好满足所研究对象的近圆太阳同步轨道的要求，故轨道保持中只需修正ω的长期变化。大气阻力摄动引起a，e的长期变化，故要进行轨道保持。倾角i没有长期变化，只有短周期变化，故不作修正。因此控制指标主要有半长轴a，偏心率e和近地点角距ω。

名义轨道的轨道要素记为半长轴a^*，偏心率e^*，近地点角距ω^*，任意时刻航天器的密切轨道要素记为a(t)，e(t)，i(t)，则轨道要素偏移量记为δa(t)＝a(t)-a^*，δe(t)＝e(t)-e^*，δi(t)＝i(t)-i^*，记可容忍的最大轨道要素偏差区间为Δa轨道保持采取偏差-修正的控制策略，即当则施加脉冲进行轨道修正，其余轨道要素的控制同理可得。

其中，Δv_r，Δv_u，Δv_h是速度脉冲在轨道系下的三个分量，其下标的意义与公式(1)中摄动力下相同，这里不再赘述。

由式(34)可知，当真近点角θ为0°或180°，摄动力在轨道坐标系的横向分量f_u仅引起a，e的变化，径向分量f_r仅引起ω的变化。这说明了在近地点和远地点施加脉冲的控制方法使轨道保持策略解耦，避免了复杂的计算，更有利于星上具体实现。因此对半长轴、偏心率采用经典的霍曼轨道转移策略，对近地点角距则在近地点的时候采用偏差修正的控制方式。将控制变量分为两组：半长轴、偏心率和近地点角距，采用分开控制的方式，从而便于程序后续的扩展，并且可将不同的控制指标分别加在不同的控制变量上。因此，在轨道保持的过程中，将出现三种工作模式：

(1)只控制半长轴和偏心率

(2)只控制近地点角距

(3)半长轴、偏心率、近地点角距全控。

对卫星轨道进行保持，采用近地点-远地点的双脉冲闭环控制策略，将卫星轨道保持在一定范围内，从而确保卫星完成相应观测任务。

步骤四：采用微分修正算法提高轨道保持的精度

由于本文中控制变量仅仅是a，e，ω。因此状态向量表示为：

约束函数F(X)表示航天器目标位置与施加机动后的位置之差，如下所示：

状态变量的雅克比矩阵，形式如下：

其余变量及步骤如发明内容中步骤四所示，因此不再赘述。

通过上述方法，最终找到合适的机动向量，从而完成精确的轨道机动。采用上面叙述的脉冲式轨道控制方法结合微分修正算法对航天器的半长轴，偏心率及近地点角距进行轨道保持控制。图1为轨道控制的结果。

步骤五：建立误差模型

本控制系统中，考虑了卫星定轨误差，发动机执行误差。定轨误差分别加在卫星的三个方向的位置和速度上，发动机执行误差则加在三个方向的发动机脉冲分量上。假设各误差均服从正态分布，均值为零，标准差取值见后续仿真分析。

步骤六：蒙特卡洛仿真

步骤6.1：带有位置误差的蒙特卡洛仿真

首先选取位置误差标准差在[0,50]m之间，每隔10m选取一个点，整体分析控制误差的走势和控制策略的能力。数据分析结果表明当标准差为7m时，轨道控制策略已基本失去了将半长轴控制在既定范围内的能力，因此没有必要分析位置误差标准差大于7m的情况，因此在后续的分析中，选取位置误差标准差在[0,7]m之间，每隔1m选取一个点，进行分析。

步骤6.2：带有速度误差的蒙特卡洛仿真

首先选取速度误差标准差在[10^-7,10^-1]m之间，整体分析控制误差的走势和控制策略的能力。数据分析结果表明速度误差标准差在10^-3m/s和10^-2m/s之间，控制误差的标准差增长迅速，当标准差为10^-2m/s时，该控制策略已基本失去了将轨道要素控制在既定范围内的能力，因此在后续的分析中，选取速度误差标准差在[10^-3,10^-2]m之间，每隔0.002m/s选取一个点，进行分析。

步骤6.3：带有发动机推力误差的蒙特卡洛仿真

发动机推力误差标准差为当前所需脉冲量乘以某个系数，该系数取为从0.1到0.6的6个参数，对每个参数进行100次蒙特卡洛仿真，并对不同误差下的轨控情况进行比较。之所以不考虑发动机推力误差标准差系数大于0.6的情况是由于该系数表征了发动机推力的误差，系数为0.6表明该随机误差的标准差是当前所需脉冲量的0.6，该数据已经很大，此时发动机的设计是不达标的，因此不予考虑系数大于0.6的情况。

步骤七：计算控制误差的均值，方差和控制误差比例三个参数，并利用这三个参数对控制结果进行分析，评估该轨道控制策略的鲁棒性。

由于不同任务对轨控策略的精度要求不同，因此本实施例中仅提出一种一般性的区间定义方法。当误差均值和标准差均小于最大容许控制范围的20％，控制误差比例在50％以内时，认为该控制策略的效果可以接受。当误差的分布区间等于或者大于最大容许控制范围，控制误差比例大于50％时，认为该控制策略失去了将航天器控制在既定范围内的能力。在这两者之间的区间认为轨道控制策略的性能较差但仍能实现对航天器轨道的粗略控制。

步骤7.1：位置误差下的鲁棒性分析；

图二表示了在不同位置误差标准差下，该轨道控制策略用于控制半长轴和近地点角距的脉冲次数。从图中可以看到用于控制半长轴的脉冲次数远大于用于控制近地点角距的脉冲次数。

图三表示了随位置误差标准差的增大，半长轴和近地点角距的控制误差的变化情况。该曲线是基于控制误差的数据点，采用MATLAB中的polyfit函数进行二阶拟合后的曲线。从图中可以看出，随着位置误差标准差的增大，半长轴和近地点角距的控制误差均值和方差均呈增大趋势。当位置误差标准差小于1.52m时，半长轴的控制误差均值和标准差分别为66m和126m。此时，控制误差主要分布在[-60,192]m之间，考虑到初始的半长轴控制范围为500m，因此得出了大多数的轨控数据点主要分布在692m的范围以内。此时，控制误差比例为36.23％。从以上三个数据中，可以得出将近40％的轨控数据点已超出了给定的轨道控制范围，但其偏差并不大，因此可以得出在位置误差均值小于1.52m时，该轨道控制方法的控制效果仍可以接受。当位置误差为4m时，误差数据的变化范围已接近轨道控制范围。控制误差的均值和标准差分别为196m和292m，另外，控制误差比例达到了58.62％。这三个数据表明该轨道保持策略已失去了对航天器的轨道进行保持控制的能力，因此当位置误差大于4m时，认为该轨控方法已失效。在位置误差标准差介于1.52m和4m之间时，该轨控策略的性能较差，但仍能实现对航天器半长轴的粗略的轨道保持作用。

近地点角距的控制误差均值和方差在位置误差较小时，增长缓慢，当位置误差标准差大于4m后增长速度有所增加。当位置误差标准差为7m时，控制误差的均值和标准差为0.39deg和0.60deg。这两个数据相对于6deg的控制范围而言，仍然很小，并且此时控制误差比例仅为23.48％，这说明仅有不到四分之一的轨控数据点超出了轨控范围。

因此，综上可知，该轨控策略在位置误差小于1.52m时能较好地将半场轴控制在500m范围内，在位置误差达到7m时，仍能精确地将近地点角距控制在6deg范围内，表现出较好的鲁棒性。

步骤7.2：速度误差下的鲁棒性分析；

图4表明了随着速度误差标准差的增大，用于控制半长轴和近地点角距的速度脉冲不断增大，其中，用于控制半长轴的速度脉冲增幅是近地点角距的8倍。

从图5可以看出，随着速度误差标准差的增大，半长轴和近地点角距的控制误差均值和标准差均不断增大。当速度误差标准差为0.00153m/s时，半长轴的控制误差均值和标准差分别为66m和125m，控制误差的比例为36.56％。根据这三个数据可以得出，该控制策略能以较小的误差将半长轴控制在给定的范围内。当速度误差标准差为0.00314m/s时，控制误差的均值和标准差分别为167m和250m,控制误差的分布范围与半长轴的控制区间近似相等。此时，控制误差的比例达到了53.87％，因此控制策略的效果很差，当速度误差标准差大于0.00314m/s时，该控制策略逐渐失去了将半长轴控制在给定范围内的能力。

近地点角距的在速度误差下的表现与在位置误差下相似。当速度误差标准差小于0.004m/s时，控制误差的均值和标准差增长缓慢，因此当速度误差标准差在该区间时，轨控数据点相对于近地点角距的控制范围偏离不大。当速度误差标准差为0.006m/s时，控制误差的均值和标准差分别为0.64deg和1.17deg，因此近地点角距的控制误差范围被扩大到3.34deg，控制误差的比例为29.39％。以上数据表明了控制策略仍然具有将近地点角距控制在给定范围附近的能力，只是控制精度较低。当速度误差标准差大于0.006m/s，控制误差的均值和标准差增长加快，当速度误差标准差是0.01m/s时，控制误差数据的分布范围和给定的控制范围相当，此时控制误差的均值和标准差分别为2.24deg和4.35deg，另外41.09％的轨控数据点已超出了给定的控制范围，因此此时该轨道控制策略已不再具有将近地点角距控制在给定范围内的能力。

因此，综上可知，该轨控策略在速度误差小于0.00153m/s时能较好地将半场轴控制在500m范围内，在速度误差小于0.004m/s m时，仍能将近地点角距控制在6deg范围内。

步骤7.3：发动机推力误差下的鲁棒性分析；

图6表明随着发动机推力误差系数的不断增大，轨道控制误差不断增大。图7中，随着发动机推力误差标准差的系数不断增大，半长轴和近地点角距的控制误差不断增大。当发动机推力误差标准差系数小于0.3时，半长轴控制误差的均值和标准差均在零值附近很小的范围内变化，轨道控制策略的效果与无随机误差下的控制效果差别不大。当发动机推力误差标准差系数是0.4时，控制误差的均值和标准差分别为71m和182m。控制误差的比例为26.09％。这说明此时轨道控制策略的精度仍较好，轨控结果可以接受。当发动机推力误差标准差系数达到50％，控制误差的均值和标准差分别为168m和340m，这说明大部分的控制误差分布在700m的范围内。此时超出控制范围的数据点占总轨控数据点的29.76％。这说明轨控策略的精度已较差，当标准差系数超过60％，轨控数据点的分布将会超过1500m。此时该轨道控制策略已失效。

对近地点角距而言，当误差系数为60％，控制误差的均值和标准差分别为0.46deg和0.99deg。并且此时仅有19.92％的数据超出了控制范围，因此可以得知此时轨控策略仍能够将近地点角距控制在给定的范围附近，轨控策略的效果仍较好。

因此，综上可知，该轨控策略在速度误差小于0.3时能较好地将半场轴控制在500m范围内，在速度误差达到0.6时，仍能将近地点角距控制在6deg范围内，表现出良好的鲁棒性。

步骤7.4：不同控制范围下的鲁棒性分析；

在不同的航天器任务中，对轨道的控制精度要求也会不同，因此轨道要素的控制范围并不仅仅局限于半长轴500m和近地点角距6deg。为不失一般性，这里同样进行了半长轴为250m和1000m，以及近地点角距为3deg和12deg的误差存在情况下的轨道控制情况仿真。表2展示了不同控制范围下的半长轴控制误差。第一行表示了三个控制范围250m，500m，和1000m，表的左侧一列是随机误差的种类和大小，其余部分是在给定控制范围和随机误差种类大小的情况下的控制误差均值和标准差。表3是表2对应控制范围和随机误差下的控制误差比例。

表2 半长轴控制误差.

表3 半长轴控制误差比例.

表三中黑体数值表示了控制误差的变化范围与该控制策略下的轨道要素控制范围近似。从表中可以看出，随着控制范围的增大，半长轴对控制误差的容忍能力也在增加。当控制范围为250m时，在位置误差标准差为2m时，控制误差的变化范围与给定的控制范围近似相等，而对于控制范围500m的情况下，则该位置误差标准差扩大到4m。速度误差的情况也类似。而对于发动机推力误差而言，当发动机推力误差标准差系数达到60％，轨道控制的性能才表现较差，当该系数小于40％，控制误差的值均比位置误差和速度误差情况下的误差均值和标准差小。由此可推断出，该轨道控制策略在发动机推力误差方面有较好的鲁棒性。当随机误差相同，控制误差的均值和方差随控制范围的增大而增大，这主要是由于轨道机动的频率下随着控制范围的增大而下降，对半长轴的控制越来越松散，因此控制误差才会不断增大。

表4 近地点角距的控制误差.

表5 近地点角距的控制误差比例.

表4表示了不同近地点角距的控制范围下的控制误差均值和标准差。表头的含义与表2类似，因此不做赘述。表5表示了表4对应情况下的控制误差比例。

与半长轴控制误差不同的是，当位置误差标准差达到6m或者发动机推力误差标准差系数达到60％，近地点角距的控制误差仍然很小，其占对应控制范围的比例仍然较小。这说明该控制策略在控制近地点角距是对位置误差和发动机推力误差有很好的鲁棒性。当速度误差标准差在0.008m/s和0.01m/s之间时，控制误差的分布范围与对应的控制范围数值相当。具体来说，对于控制范围为3deg和6deg的情况下，当速度误差标准差在0.008m/s及大于0.008m/s时，控制误差的变化范围也分别达到了3deg和6deg。而对于控制范围为12deg的情况下，该速度误差标准差则大约为0.01m/s。这说明随着控制范围的增大，该控制策略容许的随机误差也在不断增大。

根据以上给定的半长轴和近地点角距控制范围下的轨控误差结果可知，轨道控制策略对近地点角距的保持的鲁棒性优于半长轴。联系实际情况中，目前的轨道技术下位置误差和速度误差可分别达到厘米级和毫米级，因此该轨道控制策略在对低轨近圆太阳同步轨道的轨道保持中是可行并且可靠的。为保证轨控策略不失效，在半长轴和近地点角距控制范围分别为250m和3deg时，位置误差应小于2m，速度误差应小于0.002m/s，发动机推力误差标准差系数应小于60％。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种航天器轨道控制策略的鲁棒性分析方法，其特征在于：包括如下步骤：

步骤一：通过Gauss轨道要素摄动方程和平均轨道要素建立航天器轨道运动模型；

采用Gauss轨道要素摄动法建立的航天器六个轨道要素半场轴a，偏心率e，倾角i，近地点角距ω，升交点赤经Ω和真近点角θ的轨道运动模型如下所示：

$\left\{\begin{array}{c}\frac{da}{dt}=\frac{2a^2}{\sqrt{\mathrm{\μ}p}}(e{\mathrm{sin\θf}}_r+(1+e\cos \mathrm{\θ}\left)f_u\right)\\ \frac{de}{dt}=\sqrt{\frac{p}{\mathrm{\μ}}}({\mathrm{sin\θf}}_r+(\left(1+\frac{r}{p}\right)\cos \mathrm{\θ}+\frac{er}{p}\left)f_u\right)\\ \frac{di}{dt}=\frac{r\cos (\mathrm{\ω}+\mathrm{\θ})}{\sqrt{\mathrm{\μ}p}}{f}_h\\ \frac{d\mathrm{\Ω}}{dt}=\frac{r\sin (\mathrm{\ω}+\mathrm{\θ})}{\sqrt{\mathrm{\μ}p}\sin i}{f}_h\\ \frac{d\mathrm{\ω}}{dt}=\frac{1}{e}\sqrt{\frac{p}{\mathrm{\μ}}}(-{\mathrm{cos\θf}}_r+(1+\frac{r}{p}){\mathrm{sin\θf}}_u-\frac{er}{p}\sin (\mathrm{\ω}+\mathrm{\θ}\left){\mathrm{ctgif}}_h\right)\\ \frac{d\mathrm{\θ}}{dt}=\frac{\sqrt{\mathrm{\μ}p}}{r^2}+\frac{\cos \mathrm{\θ}}{e}\sqrt{\frac{p}{\mathrm{\μ}}}{f}_r-\frac{\sin \mathrm{\θ}}{e}(1+\frac{r}{p})\sqrt{\frac{p}{\mathrm{\μ}}}{f}_u\end{array}\right.---\left(1\right)$

其中，

$r=\frac{p}{1+e\cos \mathrm{\θ}}=\frac{a(1-e^2)}{1+e\cos \mathrm{\θ}}$

f_r，f_u，f_h是定义在轨道坐标系S_o的卫星摄动力在径向，横向(即在轨道平面内垂直于矢经，沿纬度幅角u增大的方向)，和副法向分量(即垂直于轨道平面，与动量矩H同向)；p是半通径，μ是地球引力常数，r是地心与航天器之间的距离；

航天器轨道要素平均化方法如下：

公式(1)的向量表达形式如下：

$\frac{d\mathrm{\σ}}{dt}=f_0\left(a\right)+f_1(\mathrm{\σ},\mathrm{\ϵ})+f_2(\mathrm{\σ},t,{\mathrm{\ϵ}}^2)---\left(2\right)$

其中

σ＝(a,e,i,Ω,ω,M)^T (3)

$\left\{\begin{array}{c}{f}_0\left(a\right)=\mathrm{\δ}n,n=\sqrt{\mathrm{\μ}}{a}^{-3/2}\\ \mathrm{\δ}={\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\end{array}\right)}^T\end{array}\right.---\left(4\right)$

这里ε＝O(J₂)即通常所取的小参数；等式右端的f₁和f₂为非球形引力加速度的部分，它们分别对应J₂项和其他非球形引力项，并有

f₁/f₀＝O(ε),f₂/f₀＝O(ε²) (5)

公式(2)是一个非常典型的小参数方程；本步骤中采用古在由秀提出的根据非线性力学中平均法的思想，针对主要带谐项摄动的平均根数法；平均根数法将摄动变化项△σ⁽¹⁾(t),△σ⁽²⁾(t)，……中不同性质的项(长期项，长周期项和短周期项)区分开，更重要的是在此前提下将参考解取为只包含长期变化项的平均根数；

对Gauss轨道要素摄动方程中积分得出的瞬时轨道要素，利用古在优秀方法计算出其平均轨道要素，从而用于后续的控制和鲁棒性分析；

步骤二：对低轨卫星的非球形摄动以及大气阻力摄动进行分析

步骤2.1：地球非球形摄动；

地球非球形摄动中，由于地球的旋转，在运行时间超过一天时，扇形谐函数和田形谐函数的影响或多或少地被“平均掉”，因此对于大多数应用情况，只需考虑带形谐函数的影响；其形式如下：

其中，J_n是n阶带形谐函数的系数，P_n和P′_n是n阶勒让德多项式及其导数，x，y，z是惯性坐标系S_i下的航天器三维坐标，且R_E是地球的平均赤道半径；因此由式(6)得到的惯性系下的非球形摄动需要经过坐标转换，以符合式(1)中对于摄动力的定义环境；

由惯性系到轨道系的坐标转换分为三步：

其中，u是纬度幅角，u＝ω+θ；

转换矩阵具体形式如下：

$L_{oi}=L_{R_{z\left(u\right)}}{L}_{R_{x\left(i\right)}}{L}_{R_{z\left(\mathrm{\Ω}\right)}}=\left[\begin{array}{ccc}\cos u& \sin u& 0\\ -\sin u& \cos u& 0\\ 0& 0& 0\end{array}\right]*\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos i& \sin i\\ 0& -\sin i& \cos i\end{array}\right]*\left[\begin{array}{ccc}\cos \mathrm{\Ω}& \sin \mathrm{\Ω}& 0\\ -\sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right]$

因此，满足公式(1)中定义的摄动力的地球非球形摄动加速度在轨道系的分量形式如下：

$f_o^E=L_{oi}\left[\begin{array}{c}{\mathrm{\Δg}}_x\\ {\mathrm{\Δg}}_y\\ {\mathrm{\Δg}}_z\end{array}\right]---\left(7\right)$

步骤3.2：大气阻力摄动；

由空气动力学可知，在大气中运动的航天器受到空气阻力为

$D=C_D\frac{1}{2}{\mathrm{\ρv}}_a^{2}S---\left(8\right)$

其中，C_D为阻力系数，ρ为大气密度，S为横截面积，v_a是航天器质心相对于大气的速度；为简化计算，假设大气是球对称的；大气不随地球旋转；航天器的迎风面积不变；于是大气阻力只沿引起沿轨道切线方向的负加速度，该加速度写在速度坐标系中如下：

$\left\{\begin{array}{c}{f}_n=0\\ f_t=-\frac{1}{2m}{C}_D{\mathrm{S\ρv}}^2\\ f_h=0\end{array}\right.---\left(9\right)$

f_n，f_t，f_h分别是大气阻力摄动加速度在速度坐标系的三个分量，其中f_t沿速度方向，f_h垂直于轨道平面指向角动量方向，f_n与其余两个向量垂直，满足右手定则；将速度坐标系S_v中的大气摄动加速度转换到轨道坐标系S_o中,即该坐标系的转换只需要绕速度坐标系的Z轴进行一次欧拉转换，转换矩阵为：

$R_z\left(\mathrm{\γ}\right)=\left[\begin{array}{ccc}\cos \left(\mathrm{\γ}\right)& \sin \left(\mathrm{\γ}\right)& 0\\ -\sin \left(\mathrm{\γ}\right)& \cos \left(\mathrm{\γ}\right)& 0\\ 0& 0& 1\end{array}\right]---\left(10\right)$

$\sin \mathrm{\γ}=e\sqrt{\frac{\mathrm{\μ}}{p}}\frac{\sin \mathrm{\θ}}{v}---\left(11\right)$

因此，满足公式(1)中定义的摄动力的大气阻力摄动加速度在轨道系的分量形式如下：

$f_o^D=R_z\left(\mathrm{\γ}\right)\left[\begin{array}{c}{f}_n\\ f_t\\ f_h\end{array}\right]---\left(12\right)$

将地球非球形摄动力和大气阻力摄动力带入公式(1)中的摄动力中，即

$\left[\begin{array}{c}{f}_n\\ f_t\\ f_h\end{array}\right]=f_o^E+f_o^D---\left(13\right)$

则得到了带有摄动的航天器轨道运动模型；

步骤三：使用轨道保持策略对航天器轨道进行控制；

名义轨道的轨道要素记为σ^*，由步骤二中的轨道动力学模型数值积分所得任意时刻航天器的密切轨道要素记为σ(t)，则轨道要素偏移量记为δσ(t)＝σ(t)-σ^*，记可容忍的最大轨道要素偏差区间为△σ，轨道保持采取偏差修正的控制策略，即当则施加脉冲进行轨道修正，其余轨道要素的控制同理可得；

步骤四：对步骤三中计算所需速度脉冲采用微分修正算法进行修正，从而完成精确的轨道机动；

步骤五：建立误差模型

定轨误差分别加在卫星的三个方向的位置和速度上，发动机执行误差则加在三个方向的发动机脉冲分量上；随机误差均服从正态分布。

步骤六：蒙特卡洛仿真分析

采用带有微分修正的轨控策略对带有摄动的航天器轨道进行控制，并加入步骤五中的误差模型，采用MATLAB中normrnd函数对具有同一控制误差标准差的航天器轨道重复多次仿真，消除随机性，模拟体现真实情况；

步骤七、计算模型

本发明中对控制误差的分析主要集中于所有未落在控制范围中的点的均值和分布情况作为控制误差的评价方法；首先计算点的均值和标准差，因此控制误差的计算方法如下所示：

$\stackrel{\‾}{\mathrm{\σ}}=\frac{\Σ{\mathrm{\σ}}^{\′}}{N}---\left(14\right)$

其中，

σ'＝0,△σ∈[c,d] (15)

$\left\{\begin{array}{cc}{\mathrm{\&sigma;}}^{\&prime;}=\left|\mathrm{\&Delta;}\mathrm{\&sigma;}-c\right|,& \mathrm{\&Delta;}\mathrm{\&sigma;}<c\\ {\mathrm{\&sigma;}}^{\&prime;}=\left|\mathrm{\&Delta;}\mathrm{\&sigma;}-d\right|,& \mathrm{\&Delta;}\mathrm{\&sigma;}>d\end{array}\right.---\left(16\right)$

△σ是轨道要素相对于其初始值的变化量；最大容许误差范围是[c,d]，其中c和d分别是控制范围的下边界和上边界；N为所有数据点的总个数；在式(15)中，当轨道要素的变化量在最大容许误差范围内时，认为控制误差为零；当轨道要素的变化量超出了控制范围，σ'为超出控制范围的控制结果点到最大容许误差边界的距离；第二个参数是控制误差标准差，它表示了控制误差的分布情况，具体计算公式如下：

$\mathrm{\&xi;}=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\&Sigma;}}{({\mathrm{\&sigma;}}^{\&prime;}-\stackrel{\&OverBar;}{\mathrm{\&sigma;}})}^2}---\left(17\right)$

第三个参数是控制误差分布比例，表示超出控制范围的点占总数据点的比例；具体计算公式如下：

$k=\frac{n}{N}---\left(18\right)$

这三个参数可充分反映控制误差的各项分布；因此，可通过计算ξ，k的值对其在误差存在的情况下的轨道控制性能进行评价，并给出该控制策略鲁棒性的分析结果，为轨控策略的实际应用提供依据和保障。

2.如权利要求1所述的一种航天器轨道控制策略的鲁棒性分析方法，其特征在于：步骤五所述的微分修正方法为：

航天器的运动方程(公式(1))可简化为

$\stackrel{\&CenterDot;}{x}=f(t,x)---\left(19\right)$

其中x是描述航天器运动的状态量；

假设航天器的初始位置状态向量为x₀，在施加轨道机动△V的情况下到达位置为x^t，而其目标到达位置为x^d；

约束函数F(X)表示航天器目标位置与施加机动后的位置之差，如下所示：

F(X)＝x^t-x^d (20)

当F(X)＝0，则表示当前施加的轨道机动△V能精确地将航天器机动到目标位置；而为了寻找合适的△V，采用牛顿迭代法进行实现；根据Taylor多项式的展开法则，省略后面的高阶项，可以得到如下表达式：

F(X^j)+DF(X^j)(X^j+1-X^j)＝0 (21)

其中DF(X)表示约束函数对状态变量的偏导数，即状态变量的雅克比矩阵，形式如下：

$\mathrm{\&Phi;}(t^t,t_0)=DF\left(X\right)=\frac{\&part;F\left(X\right)}{\&part;X}---\left(22\right)$

该矩阵表示了轨道要素从初态转移到末态的状态转移矩阵，具体的数值大小可通过给定初值的函数积分得到；

给定一个合理的初值，将式(22)带入式(21)中，则可通过迭代计算得到新的X^j向量，理论上

||F(X^j+1)||<||F(X^j)|| (23)

于是，通过数次迭代计算，可使F(X)函数的值趋向于0；当

||F(X^j+1)||<ε (24)

其中||F(X^j+1)||是F(X^j+1)的模；ε是一给定的极小量；通过上述方法，最终找到合适的机动向量，从而完成精确轨道机动。