康士朋，江 涛，耿盛韦，王金童



基于蒙特卡洛法的卫星分离姿态计算方法

康士朋1，江 涛1，耿盛韦2，王金童1

（1. 上海宇航系统工程研究所，上海，201109；2. 南京航空航天大学航空宇航学院，南京，201106）

为了计算以弹簧作为分离动力源的卫星在星箭分离时刻的分离姿态，建立了六自由度动力学方程，并采用龙格库塔法求解了该方程，算例结果表明：通过该方法计算得到的分离角速度与ADAMS仿真计算结果相比，最大偏差小于0.77%。基于以上研究，针对卫星、弹簧装置多个设计偏差对卫星分离姿态的影响，提出基于蒙特卡洛法的卫星分离姿态分析方法。飞行试验结果表明，由该方法得到的分离角速度与常规保守算法相比，计算精度提高35%左右，解决了工程中由于采用将各个偏差取最大值进行保守计算而导致计算结果偏差较大的问题，可为卫星设计、弹簧装置设计提供参考。

卫星；分离姿态；弹簧装置；蒙特卡洛法；龙格库塔法

0 引 言

随着计算机技术、微纳技术及微机电系统技术的迅速发展，质量在10~100 kg范围的微纳卫星以其研制周期短、低成本、高性能等优势，在通信、军事、遥感与对地测绘、空间科学试验等领域表现出较好的发展前景[1]。

微纳卫星通常以搭载或一箭多星的方式进行发射，以弹簧装置作为分离动力源，以一定的分离速度、分离角速度与运载火箭进行分离[2]。卫星分离过程中，除了不能与火箭发生碰撞，还必须将分离角速度控制在允许范围内，过大的分离角速度甚至会导致卫星无法控制姿态，直接影响飞行任务成败。

卫星分离姿态计算方法主要包括两种：数值计算方法和基于ADAMS软件的仿真计算方法。

数值计算方法包括两种：a）建立常质量航天器动力学方程，并采用龙格库塔方法进行数值求解，如：付碧红[3]、卢丽颖[4]分别以搭载星、子星为研究对象，建立并求解了三自由度动力学方程，Jeyakumar，BN Rao[5]也通过龙格库塔法得到了分离体的非线性常微分方程数值解；b）建立拉格朗日动力学方程，如：王鑫等[6]采用拉格朗日动力学原理建立了分离体姿态动力学模型，探索了分离体姿态动力学方程组数值求解方法。而随着计算机仿真技术的发展，众多学者[7~11]提出了采用ADAMS软件进行了卫星分离仿真分析，ADAMS软件基于多刚体动力学方程，采用自动变阶、变步长方法进行计算[12]。

但相比上述计算模型，实际情况中存在两类偏差：一类是卫星重量偏差、质心偏差及转动惯量偏差；另一类是弹簧装置推力偏斜、推力大小偏差及安装误差。这些偏差因素都将影响卫星的分离姿态。

在工程应用中，为考虑以上偏差因素对分离姿态的影响，通常采用的方法是将每一个偏差量都取到最大，计算极限的分离姿态。然而，用所有偏差因素的极限值计算得到的极限分离姿态产生概率非常小，远远超出了航天可靠性要求。因此，使用以上方法计算结果偏于保守。

使用以上保守的计算方法时，为了满足卫星分离指标要求，通常采取两种措施：一种是减小以上因素的设计偏差值；另一种是提高弹簧装置生产精度，并大量投产，选配出弹簧推力偏差较小的弹簧装置。以上方法在限制了卫星与弹簧装置结构设计的同时也提高了产品生产成本。

因此，本文以采用弹簧装置进行分离的微纳卫星为例，从概率分布的角度，提出基于蒙特卡洛法的考虑设计偏差影响的卫星分离姿态计算方法，并通过飞行试验验证了该方法的有效性。

1 动力学方程建立与求解

1.1 建立动力学方程

以CZ-4B运载火箭搭载发射微纳卫星（X卫星）为例，建立六自由度常质量航天器动力学方程。该卫星在4个弹簧装置作用下与运载火箭进行分离，分离方案如图1所示。

图1 卫星分离方案

由于卫星质量、转动惯量远远小于运载火箭末子级质量、转动惯量，因此分离过程中，不考虑运载火箭对卫星分离姿态的影响。

将坐标系设定于星箭分离面中心，固定在运载火箭上，以图1模型分离系统为研究对象，建立4个弹簧装置作用下的分离姿态计算分析模型，如图2所示。

图2 分离姿态计算模型

Δ，Δ，Δ—卫星质心在，，轴的偏移量；1，2，3，4—考虑安装偏差后的4个弹簧的安装距离，其中3，4为负值；F第个弹簧装置推力在向分力（=1,2,3,4）；F第个弹簧装置推力在向分力；F第个弹簧装置推力在向分力

考虑弹簧推力偏斜影响，以单个弹簧为研究对象，根据总体坐标系，将每一个弹簧装置的推力进行分解，如图3所示。

图3 弹簧推力分解示意

F—第个弹簧装置推力；—第个弹簧推力与轴夹角；—第个弹簧力在平面投影力与轴夹角

根据图3可以得到每个弹簧推力在坐标方向的分解力为

根据图2和式（1）可以得到使卫星绕三轴旋转力矩M，M，M的计算公式为

根据牛顿第二定律和动量矩定理，建立卫星与运载火箭末子级两个刚体之间的六自由度动力学方程：

1.2 动力学方程求解

上述方程维数高、变量多、且互相耦合，无法求解其解析解，可采用龙格库塔方法进行数值计算求解。

在计算过程中，弹簧推力、力矩为变化值，为了解决这一问题，本文提出将弹簧工作时间划分为多个连续子区间，在每一个子区间内，将该区间中间时刻的弹簧力作为该子区间内的弹簧推力。

首先计算弹簧装置工作时间，以4个弹簧组成的分离系统为研究对象，由于弹簧质量远小于卫星质量，卫星运动为弹簧振子简谐运动。

在初始时刻0=0，初始速度为0，弹簧变形量达到最大值max，max即为弹簧振子振幅；在弹簧工作结束时刻s，弹簧变形量降低至最小值min：

式中为分离系统简谐运动角频率。

根据以上公式计算弹簧工作时间为

式中k为第个弹簧的刚度。

式中为第个弹簧力角频率；为弹簧效率，根据工程经验取值。

本文基于MATLAB软件，采用龙格库塔数值计算方法来求解公式，求解过程如图4所示。

图4 动力学方程求解流程

1.3 计算方法验证

为了验证本文计算方法的正确性，针对同一组计算参数，采用该方法和ADAMS仿真计算两种方法进行计算，并将计算结果进行比对。

参照X卫星给出一组设计参数，如表1所示，进行卫星分离姿态计算。

表1 计算参数汇总

续表1

4#弹簧力/NFmax=250；Fmin=80 1#弹簧位置/mm107.05 2#弹簧位置/mm107.08 3#弹簧位置/mm106.99 4#弹簧位置/mm107.0 1#弹簧推力角度/(°)θ1=0.762；φ1=30 2#弹簧推力角度/(°)θ2=0.742；φ2=45 3#弹簧推力角度/(°)θ3=0.758；φ3=60 4#弹簧推力角度/(°)θ4=0.750；φ4=90 弹簧效率η0.85

注：max—弹簧装置最大推力；min—弹簧装置最小推力

根据式（1）计算得到弹簧工作时间为0.074 2 s，再将该时间划分为50份，根据1.2节公式，采用MATLAB软件对卫星分离姿态进行计算，计算结果以分离角速度为例，如图5和表2所示。

图5 数值计算结果

表2 数值计算结果汇总

在ADAMS软件中，建立仿真分析模型，卫星的质量、惯量和质心位置同表1。在地面建立总坐标系，在弹簧推力作用点施加弹簧推力，该采用行程与刚度函数表示。

设置计算时间为0.1 s、载荷步数量为50，进行计算，结果如图6、表3所示。

图6 ADAMS姿态角速度结果

表3 仿真计算结果汇总

将由两种方法得到的计算结果进行对比，得到两种方法相对差如表4所示。

表4 计算偏差汇总

经过表4对比可得，两种计算方法得到的计算结果最大相对差为0.77%，这是由于将子区间中间时刻弹簧推力简化为该子区间弹簧推力而造成的。本算例证明了1.2节计算方法的正确性。

2 考虑设计偏差的分离姿态计算方法

2.1 计算方法

蒙特卡洛方法不同于确定性的数值方法，它是通过模拟随机变量来解决数学问题的非确定性的（概率统计的或随机的）数值方法[13]。

因此，本文提出基于蒙特卡洛方法的卫星分离姿态评估方法：首先根据各个设计变量的偏差范围随机产生大量的设计参数；再根据第1节方法进行卫星分离姿态计算，将得到大量卫星分离姿态计算结果；再对这些计算结果进行统计，得到其分布规律；最后根据概率分布范围，找出合理的极限分离姿态。分析方法流程如图7所示。

图7 卫星分离姿态分析方法

2.2 计算方法验证

以CZ-4B运载火箭搭载发射X卫星为例，验证该方法有效性。

卫星、弹簧装置完成生产后，对卫星的质量、质心位置、转动惯量进行实测；对弹簧装置的推力大小、安装位置进行实测，如表5所示。

表5 产品实测参数

由于弹簧装置在设计过程中为防止弹簧推杆与弹簧壳体发生卡滞，两者之间为间隙配合，因此弹簧推力是在一定角度范围偏斜。根据弹簧推杆、弹簧壳体的加工参数，计算得到弹簧最大推力偏斜为0.622°。

针对这一项设计偏差，将采用传统保守算法以及第2.1节提出的算法对X卫星的分离姿态进行计算，并将两种方法的计算结果与飞行试验结果进行对比。

按照保守算法，将推力偏斜设置为最大值，确定卫星分离最严酷工况，如表6所示。

表6 保守算法计算参数

根据表5、表6的设计参数，在ADAMS软件中进行卫星分离姿态仿真分析，建模方法同1.3节，以分离角速度为例，计算结果为：2.892 (°)/s，-3.0 (°)/s，0.970 (°)/s。

根据2.1节提出的评估方法，设置弹簧推力偏斜角度为：∈[0,0.622°]，∈[0,360°]，在MATLAB计算程序中产生100 000个随机变量，通过计算得到100 000组计算结果，以分离角速度为例，对这些计算结果进行处理，通过程序判断计算结果服从正态分布规律，如图8至图10所示，并得到绕三轴分离角速度均值与方差，如表7所示。

图8 绕X轴分离角速度分析结果

图9 绕Y轴分离角速度分析结果

图10 绕Z轴分离角速度分析结果

表7 分离角速度均值与方差

据表7结果，按照3原则，得到置信度为99.73%的卫星绕三轴分离角速度：1.822 (°)/s，-1.940 (°)/s，0.594 (°)/s。

2014年，该卫星通过CZ-4B运载火箭搭载发射成功，根据卫星回传数据，星箭分离时刻，卫星绕三轴分离角速度分别为：0.8 (°)/s，-0.6 (°)/s，0.17 (°)/s。

对以上数据进行汇总，如表8所示。

表8 计算、试验结果数据汇总

对比表8数据可得，由第2.1节方法得到的计算结果与常规保守法计算结果相比，计算精度分别提高37%，35%，39%，该方法可以代替传统的保守算法。

综合表7和飞行结果，如果按照1σ原则，得到置信度为68.27%的分析结果，0.7 (°)/s、-0.806 (°)/s、0.198 (°)/s，不能包络飞行试验结果；如果按照2σ原则，得到置信度为95.45%的分析结果，1.261 (°)/s、-1.373 (°)/s、0.396 (°)/s，可以包络飞行试验结果；说明至少按照2σ原则给出的计算结果才可以满足使用要求，因此本文按照3σ原则给出的计算结果。

3 结 论

本文建立了考虑推力偏斜的六自由度动力学方程，采用龙格库塔方法求解该方程，算例结果表明，由该方法得到的分离角速度结果与仿真计算结果比最大相差为0.77%，验证了该方法正确性。

基于以上研究，提出基于蒙特卡洛法的考虑设计偏差影响的卫星分离姿态计算方法。算例结果表明，卫星分离角速度计算结果服从正态分布规律。飞行试验结果表明，要至少按照2σ原则给出的计算结果才可以满足使用要求，按照本文提出的3σ原则给出的计算结果与常规保守算法相比其计算精度提高35%左右。

该方法解决了工程中由于采用将各个偏差取最大值进行计算而导致计算结果偏差较大的问题,该方法可为卫星设计、弹簧装置设计提供参考。

[1] 石荣, 李潇, 邓科. 微纳卫星发展现状及在光学成像侦察中的应用[J]. 航天电子对抗, 2016: 32(1): 8-13.

[2] Liu L K, Zheng G T. Parameter analysis of PAF for whole-spacecraft vibration isolation[J]. Aerospace Science and Technology, 2007(11): 464-472.

[3] 付碧红, 杜光华. 搭载星与运载火箭分离的动力学研究[J]. 飞行力学, 2006, 24(1): 55-58.

[4] 卢丽颖, 孟宪红, 邢依琳. 卫星空间分离动力学研究[J]. 动力学与控制学报，2014, 12(2): 165-169.

[5] Jeyakumar D, Rao B N. Dynamics of satellite separation system[J]. Journal of Sound and Vibration, 2006, 297(1-2): 444-455.

[6] 王鑫, 袁晓光, 杨星. 基于拉格朗日方法的飞行器多体分离姿态动力学分析研究[J]. 西北工业大学学报: 2014, 32(1): 18-22.

[7] 王秋梅, 孟宪红, 杨庆成. 卫星二次分离方案仿真研究[J]. 系统仿真学报, 2010, 22(9): 2217-2222.

[8] 沈晓凤, 肖余之, 杜三虎, 等. 基于蒙特卡罗方法的小卫星偏心分离动力学分析[J]. 上海航天, 2014, 31(1): 12-17.

[9] 沈晓凤, 肖余之, 康志宇. 小卫星偏心分离动力学仿真模型的建立与验证[J]. 飞行力学, 2012, 30(3): 258-262.

[10] 曾文花, 吴琳娜. 多星与上面级非对称分离技术研究[J]. 上海航天, 2014, 31(2): 30-36.

[11] 张兵, 岑拯. 多星分离的ADAMS仿真[J]. 导弹与航天运载技术, 2004(2): 1-6.

[12] 胡星志. 小卫星星箭分离系统设计、分析与优化研究[D]. 长沙: 国防科技大学, 2012.

[13] 金畅. 蒙特卡洛方法中随机数发生器和随机抽样方法的研究[D]. 大连: 大连理工大学, 2005.

Calculation Method of Satellite Separation Attitude Based on Monte Carlo Method

Kang Shi-peng1, Jiang Tao1, Geng Sheng-wei2, Wang Jin-tong1

(1. Aerospace System Engineering Shanghai, Shanghai, 201109; 2. College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, 201106）

In order to calculate the separation attitude of the satellite which is separated by spring device, the 6 degrees of freedom dynamic equation is established and solved by the Runge-Kutta Method. The numerical examples indicate that the maximum deviation of separation angular velocity obtained by this method and the ADAMS simulation is less than 0.77%. Based on the above studies, a analysis method based on Monte Carlo Method is proposed to calculate satellite separation attitude considering the satellite and spring device design deviations. The results of flight test indicate that the proposed method can improve the accuracy by 35% compared with the conventional conservative algorithm. The method solves the problem that the deviation of the calculation result is large due to the conservative calculation of the maximum value of each deviation. This method can provide a reference for satellite design and spring device design.

Satellite; Separation attitude; Spring device; Monte Carlo method; Runge-Kutta method

1004-7182(2017)06-0036-06

10.7654/j.issn.1004-7182.20170609

V412.4+2

A

2017-06-04；

2017-07-11

康士朋（1983-），男，博士研究生，工程师，主要研究方向为箭体结构和星箭连接分离装置设计