航天器捕获制动过程的质量特性在轨辨识方法研究

韩 宇，曹 涛

（1.上海航天控制技术研究所，上海 201400；2.上海市空间智能控制技术重点实验室，上海 201400）

航天器捕获制动过程的质量特性在轨辨识方法研究

韩 宇1，2，曹 涛1，2

（1.上海航天控制技术研究所，上海 201400；2.上海市空间智能控制技术重点实验室，上海 201400）

针对航天器系统质量、质心位置和惯性矩阵的在轨辨识问题提出一种解决方法：将系统所有未知参数以组的形式进行划分，每组未知参数都可以转化为线性表示形式，从而将一个非线性系统的参数辨识问题转化为若干个线性的子参数辨识问题；用递推式最小二乘法对每个子参数辨识问题进行求解，在对某组参数求解时所需的其他未知参数则用其估计值代替。通过该方法可将复杂非线性系统转化为若干线性系统实现系统参数在轨辨识。通过数值仿真对采用推进器激励的航天器的总质量、质心位置和惯性矩阵进行辨识，验证了方法的有效性。

捕获制动；最小二乘；质量特性辨识；推进器激励

1 引言

地外天体探测变轨捕获制动过程中，推力施加的效果和姿态的变化会使得航天器的质量特性发生变化［1］。进行地外天体探测的航天器往往本身质量较大，而捕获制动过程由于近心点的位置而时间较短，变轨速度增量大，所以都采用具有推力大、比冲小、推进剂消耗多的大推力发动机进行大推力变轨，以保证任务成功率和对近心点高度范围的控制［2］。在发射前利用地面手段计算出航天器的质量特性仅为近似准确的标准数据［3］，当航天器运行时，燃料消耗、结构变形（如天线等）、以及对航天器进行在轨维修导致的潜在载荷消耗的影响，均会导致航天器的质量特性发生改变，无法获得高精度航天器质量特性［4］。因此，对捕获制动过程的航天器进行高效可靠的在轨辨识是获取质量特性的重要途径。

目前为止，有不少学者就航天器质量特性在轨辨识这一问题开展研究。早在上世纪80年代，Bergmann等［5］就基于高斯二阶滤波对航天器质量特性进行了辨识，然而该算法过于复杂不适合在轨应用，且辨识模型中忽略了ω×Iω项，影响了模型准确性。Wilson等［6］将航天器运动学近似线性化，利用速率陀螺的测量数据，实现系统的惯量和质心位置的在线辨识，但不能保证模型的精确性。国内这一领域的研究不多，且大部分仅能实现对单一质量特性的辨识。刘伟霞等［7］通过星上陀螺测量角速率信息，基于扩展卡尔曼滤波和最小二乘方法，用两步对航天器的转动惯量完成了辨识；黄河等［8］通过闭环控制来辨识航天器转动惯量；朱东方等［9］基于扩展卡尔曼滤波，考虑柔性附件对航天器姿态的影响，辨识了航天器的转动惯量。目前航天器质量参数在轨辨识仍存在两个重要问题：1）航天器建模的准确性，辨识中需考虑系统中存在的耦合情况，以提高辨识精度；2）需要同时对航天器质量、质心位置、系统惯量进行辨识。

针对以上问题，本文提出一种对一般非线性系统均有效的在轨辨识算法：将未知参数按组划分，从而将对未知参数的在轨辨识问题转化为对每组参数的辨识问题。在对每组参数的辨识过程中应用其他参数的估计值进行计算。将该方法应用于捕获制动过程中的质量参数在轨辨识问题可以航天器的质量参数划分为质心位置、转动惯量矩阵以及转动惯量矩阵的逆四组子问题，然后通过推进器激励，利用角速度和线加速度的采样信息来进行在轨辨识。

2 辨识算法

本节提出一种对最小二乘估计的改进算法，使其可有效地在轨辨识一般形式的非线性系统，并通过一个例子来说明算法特点及普遍性。

线性最小二乘估计问题的标准形式可以记作式（1）［8］：

其中，b为无噪声测量矢量，ε为测量噪声矢量，x为待识别的参数，矩阵A包含了系统已知的变量和参数，处于无噪声状态。“≅”表示无噪声情况下相等。最小二乘估计的解＾x可以使得误差A＾x-b的平方和最小［10］，求解如式（2）：

大多数问题并不能刚好转化为标准形式Ax＝b+ε，如矩阵A中存在噪声，或待辨识参数x并不能如标准形式般在已知参数A和测量矢量b间形成线性关系。以往多为直接将所求问题模型转化为最小二乘标准形式，如去掉未知参数的耦合项等［5］。但对航天器去掉系统非线性部分会导致模型不精确。为了更直观地描述这一情况，下面举例说明。考虑某系统如式（3）：

式中，b为测量值，c1、c2和c12均为已知数值，x1和x2为待辨识的未知参数。转化为辨识方程令求解。考虑到系统噪声的存在，元素三会与元素一和元素二的乘积产生偏差，从而导致系统变形。若忽略元素三的存在，又与原系统不符。

针对这类问题，本文提出一种方案。对式（3）所示系统，可将其视作两个方程，第一个方程中将参数x2视作已知变量，仅需识别x1；第二个方程则正相反。则系统方程转化为式（4）：

其中，第一个方程待辨识参数为x1，将x2视作已知参数，取其最优化估计值＾x2，第二个方程则正相反，则可分别求解＾x1、＾x2。换言之，若对系统的两个参数＾x1、＾x2分别求解，考虑每个参数的估计值是充分准确的，则可在对某个参数辨识的过程中共享其余参数的估计值。式（4）的辨识结果可以通过递推最小二乘估计得出。

为验证该方法可行性，本文进行如下仿真：令c1＝2、c2＝3和c12＝5，考虑系统输入白噪声的功率为0.1，按照前文所提出的两种情况进行仿真：第一种情况A＝［c1，c2，c12］，x＝第二种情况A＝仿真结果如表1所示。

表1 数学仿真结果Table 1 Results of Mathematical Simulation

由仿真结果可以看出，待辨识参数存在耦合时，采用其他待辨识参数的估计值的算法准确度更高。求解实际问题时，所建系统模型往往不能完全描述系统真实工作状态，常存在一些系统模型中未包含的附加影响。因而若上述方法在计算过程中所产生的误差远小于这些附加影响产生的误差，该算法就值得进一步研究。

这种系统参数的辨识方法可以拓展应用于任意未知参数x的情况。其中，x对其中包含的数组量不限制，而每个数组中参数个数也不限制。该辨识方法对其中使用的估计算法不加以限制，上文采用的是递推最小二乘法，但若采用其它估计方法也不会对系统的结构造成影响。

x中存在多个数组的情况下的辨识方法如式（5）所示，其中下标表明该参数所处的组别，而每组的参数个数不定。

其中，第一行为对未知参数x1进行辨识的递推最小二乘估计（RLSID）的回归方程，其余方程均类似。

采用RLSID作为估计算法时，对某参数辨识时需要应用到其余参数的估计值，这就使得结果的精确度与初值的选择有着直接的关系。

考虑实际应用的情况，通常一个系统待辨识参数中每个参数的不确定性不同。以航天器为例，其质量、质心位置参数、转动惯量矩阵对角线元素等，均会由于地外天体捕获过程中燃料的消耗、推力施加的效果而变化，因而待辨识。而这几个待辨识参数中如质量主要与燃料消耗相关，可通过喷气时长等估算，相对而言不确定性较低；而质心位置参数与燃料消耗、结构变形等情况均相关，相对更难估算，不确定性较高。由此易知，以质量的估计值来求解质心位置，其结果的精确性可得到提高，反之则精确性会降低。由此，对高不确定性参数辨识的过程中应用低不确定性的参数的估计值可以提高结果的精确性。

考虑到测量值往往与一些参数更加直接相关，所以本章所提出的这种参数辨识方法对系统更具合理性。例如，在航天器的推进器点火时仅产生力矩，陀螺测量的转动力矩相对于航天器质心是相对独立的。在测量质心位置时，这组数据更多的是基于测量噪声而非物理现象。因此，在对质心位置辨识时采用的测量数据中应尽量避免对陀螺测量数据的使用。由于本章所提出的方法会对几个待辨识数据分别进行处理，可以更加合理的选择和应用测量数据。

如果待辨识参数的初值在设定时不确定性过高，当测量噪声偏高时，可能导致估计值中的某个值或某些值的偏离。此时就彰显了估计值独立性的重要性，因为具备独立性的估计值可以帮助系统进行异常值检测，并从开始进行预防。如果估计误差的协方差设定合理，且测量值不受异常值影响，则毋需考虑这一情况。

3 航天器质量参数在轨辨识

3.1 航天器系统建模

航天器捕获制动过程包含姿态运动和平动运动，针对这两种情况分别介绍如下。

考虑航天器配备飞轮和推进器作为执行机构，并装配了三轴速率陀螺和加速度计测量提供航天器角速度和线速度。以轨道坐标系为参考系，记本体坐标系相对于轨道坐标系的轨道角速度为ω，且ω∈，则航天器的姿态动力学方程可写为式（6）：

其中，I∈ RR3×3表示航天器的转动惯量矩阵，需要根据实际质心位置进行测量；τ∈ RR3表示航天器所受到的力矩。

当推进器作为航天器的执行机构时，航天器所受到的总力矩可以记为式（7）：

其中，n航天器搭载的推进器数量，Li∈ RR3表示第i个推进器在体坐标系下的x-y-z的位置信息；Di∈ RR3表示体坐标系下第i个推进器的推力方向的单位矢量；Ski为一个标量，表示第i个推进器的幅值比例参数，包括排气的影响以及多个推进器同时点火时推力下降的影响；Fnom，i为常数，其值对应了相应推进器名义上提供的推力；Fbias，i表示第i个推进器的恒定非名义推力；Frandom，i表示第i个推进器的脉冲稳定非名义推力；Tki值为0或1，表示在k时刻第i个推进器是否点火的有效值；τdisturb∈ RR3表示由外界干扰源（如拖拽、重力梯度、CMG、RWA等）作用于航天器的所有力矩的总和。

考虑航天器配备加速度计，可以测得航天器在本体坐标系下的平移加速度为x¨body，且x¨body∈，则航天器的平动方程可写为式（8）：

综上，由推进器作为执行器的航天器转动方程和平动方程可记为式（10）、（11）：

3.2 质量参数辨识

航天器转动方程中待辨识的参数包括质量m、质心位置（包含L）、转动惯量矩阵以及其逆。但是这些参数之间相互存在耦合运算，且无法直接简化为线性形式，即Ax≅b。通过前一节所提出的参数辨识算法，将航天器质量辨识问题转化为几个闭环的子参数辨识问题，易保证每个子问题都是可以简化为线性形式并通过RLSID进行求解。

根据我国国内官方报道以及巴基斯坦国内的报道，以及出台的相关安全报告中，很少有中国人员（投资人员、项目工人）在巴基斯坦受到恐怖组织或者恐怖分子的袭击，仅有的专门针对中国在巴投资的人员的袭击更是为数不多。从2004年以来，我在巴基斯坦利益受到20多起恐怖袭击，仅仅是个位到十位的变化（发生次数）。

参数辨识过程中的初始值需选择最优估计值（如标称值），对每个参数的估计误差的协方差矩阵的设定会依据初始标称值的可信度。RLSID在更新的过程中会考量敏感器误差协方差矩阵，每个子辨识问题都需要应用到其他参数的最新估计值，例如在对质心位置进行估计时，需要应用到转动惯量的最新估计值。质量参数辨识的优化通过对系统方程中其他对估计值具备显著影响的参数的拓展来实现。

影响航天器质量特性的一些参数可以准确得到，如燃料的损耗变化，该参数可以由点火时间（Burn-Time-Integration，BTI）计算得到。根据每个推进器的BTI数据计算其耗损的燃料，综合所有推进器的燃料损耗可以推算航天器质量、质心位置、转动惯量等质量特性。然而航天器质量特性估计也需要考虑BTI存在偏差的情况，BTI的偏差在质量参数辨识过程中可以拓展为一个或几个未知参数的变化。

仅受力矩作用时，航天器仅存在转动运动而不存在平动运动，此时仅需测得的航天器角速度便可辨识得到航天器的转动惯量。此时转动惯量的辨识结果是相对独立的，由第2节可知，这种相对独立的辨识结果更加精确，因此转动惯量的估计值在航天器仅受力矩作用时进行更新。而若要辨识航天器的质心位置，航天器必受到平移的力因此质心位置的估计值仅在受到平移力时更新。

下文提出一种估计算法，通过陀螺提供的角速度测量信息，实现对航天器转动惯量矩阵、转动惯量的逆矩阵以及质心位置的在轨辨识。质心位置C决定了航天器体坐标系的原点，也决定了各个推进器在体坐标系的位置，即L的值。同样的，质心位置的偏差ΔC，即实际质心位置C与标称质心位置Cnom之间的差也会对L产生影响。Cnom的取值可以确定，包括由综合点火时间计算而得的相应变化值；Lnom的值可以由此计算而得；C、L和ΔC并不能准确得知；可以由参数辨识算法得到估计值，从而计算得到和。对转动惯量矩阵和转动惯量的逆矩阵的求解与其类似。如式（12）～（15）所示。

应用前文给出的参数辨识算法，将转动方程（10）转化为三个子方程，每个子方程中包含一个待辨识参数（C、I、I-1），并化简为可以进行最小二乘估计进行求解的标准形式。三个子方程中相应的变量A、x和b分别如式（16）～（18）所列：

仅用陀螺测量值对质心位置估计的标准形式的方程变量为式（19）：

若将矢积量作为干扰项处理，得到的转动惯量逆矩阵的估计的标准形式方程为式（21），转动惯量矩阵的估计值标准形式为式（22）：

但若将矢积量作为干扰处理，则需要尽可能的剥离矢积量与待辨识参数之间的联系，而不能对矢积量加以应用。而考虑航天器模型中矢积量的存在是有重要意义的，因此，下面给出将矢积量作为相关项处理，所得到的转动惯量矩阵的估计的标准形式方程如式（23）：

考虑航天器惯量特性直接与回转量相关，且可依据航天器转动更新测量值，因此辨识航天器的惯量特性时，首选式（23）。然而，如果航天器的对称性或角速率特性使回转量可忽略，式（23）就不具备特别的运算优势了。

上文所述的方程中均包含其他变量的估计值。因此，当k取值发生变化时，每个估计值均需要更新，并相应更新标准形式方程中矩阵Ak、向量bk，从而依据RLSID实时求解。

4 仿真验证

为对所提出的方法进行验证，在Matlab／Simulink中建立了航天器的动力学模型。设定航天器的待辨识质量特性为：

航天器装配的12个推进器布局如图1，图中的箭头方向为喷气矢量方向，与产生的推力方向相反。各推进器的幅值均为2 N，安装位置及推力矢量方向如表1。同时配备了陀螺仪和加速度计测量系统的角加速度和线加速度。

模拟航天器地外天体捕获制动过程开展仿真，仿真时间为100 s，初始姿态机动，机动过程航天器的姿态由＝[000]T转为＝[6°-7°8°]T，待姿态稳定后开始参数辨识阶段，在(20，50)s以推进器激励，推进器的作用顺序为：f1，f2→f7，f8→f3，f9→f5，f11→f4，f12→f6，f10。各推进器均以最大推力工作，每对推进器作用时间为5 s，总机动时间为30 s。由表2可知，该过程中的航天器仅受力矩作用，航天器的惯性矩阵辨识在这段时间进行，在（50，60）s时间段推力器作用顺序为f1→f2，每个推进器作用时间为5 s。这段时间内航天器受力作用，可以对航天器质量及质心位置进行辨识。

表2 推进器安装位置及推力矢量Table 2 The positions and vectors of thrusters

仿真过程中的航天器角速度和欧拉角变化情况分别如图2、图3所示，可见航天器的姿态和角速度在机动过程后可达到稳定状态，在参数辨识阶段，航天器的角速度量级较小，有利于提高辨识精度。

各质量特性的辨识误差如表2所示，质量和质心位置的辨识精度可达10-3，惯性矩阵的辨识精度可达10-4。

表3 质量特性辨识误差Table 3 Error of mass characteristics identification

5 结论

本文基于最小二乘估计算法提出的在轨辨识方法，通过递归的方法减小各个待辨识参数间的耦合问题，可辨识所有质量特性参数，在验证算例中质量和质心位置的辨识精度可达10-3量级，惯性矩阵的辨识精度可达10-4量级。

本文的仿真过程主要目的是对所提出的辨识方法进行原理性验证，因此没有引入测量误差。未来将进一步考虑工程性应用问题，根据实际航天器研制任务中的具体情况，加入测量误差，以对算法的性能进行更充分的验证，并加以工程实现。

（References）

［1］ 麻永平.绕月探测飞行控制［M］.北京：国防工业出版社，2010.Ma Yongping.Spacecraft Control of Circumlunar Exploration［M］.Beijing：National Defense Industry Press，2010.（in Chinese）

［2］ 刘玥，荆武兴.利用虚拟卫星法求解火星探测器近火点制动策略［J］.哈尔滨工业大学学报，2013，45（1）：14-18.Liu Yue，Jing Wuxing.Mars probe near-center braking strategy using virtual satellite method［J］.Journal of Harbin Institute of Technology，2013，45（1）：14-18.（in Chinese）

［3］ 王洪鑫，徐在峰，赵科，等.航天器质量特性测试技术新进展［J］.航天器环境工程，2011，28（2）：171-174.Wang Hongxin，Xu Zaifeng，Zhao Ke，et.Recent advances of mass property measuring technology for spacecraft［J］.Spacecraft Environment Engineering，2011，28（2）：171-174.（in Chinese）

［4］ Mohan S，Miller D W.Operational impact of mass property update for on-orbit assembly［C］／／AIAA SpaceOps 2006 Conference，Rome，Italy，AIAA 2006-5658，June 19-23，2006.

［5］ Bergmann E V，Walker B K，Levy D R.Mass property estimation for control of asymmetrical satellites［J］.Journal of Guidance，Control，and Dynamics，1987，10（2）：483-492.

［6］ Wilson E，Lages C，Mah R.On-line gyro-based，mass-property identification for thruster-controlled spacecraft using recursive least squares［C］／／The 45th Midwest Symposium on Circuits and Systems，Tulsa，Oklahoma，Aug.4-7，2002.

［7］ 刘伟霞，熊智，郁丰，等.组合航天器转动惯量在轨两步辨识标定［J］.中国空间科学技术，2013（2）：32-39.Liu Weixia，Xiong Zhi，Yu Feng，et al.On-orbit calibration technique based on the two-step moment of inertia identification of the combination spacecraft［J］.Chinese Space Science and Technology，2013（2）：32-39.（in Chinese）

［8］ 黄河，周军，刘莹莹.航天器转动惯量在线辨识［J］.系统仿真学报，2010，22（5）：1117-1120.Huang He，Zhou Jun，Liu Yingying.On-line identification of spacecraft moment of inertia［J］.Journal of System Simulation，2010，22（5）：1117-1120.（in Chinese）

［9］ 朱东方，王卫华，宋婷，等.复杂挠性航天器转动惯量在线辨识算法研究［J］.上海航天，2015，32（5）：1-8.Zhu Dongfang，Wang Weihua，Song Ting，et al.On-line identification of flexible spacecraft moment of inertia［J］.Aerospace Shanghai，2015，32（5）：1-8.（in Chinese）

［10］ 石贤良，吴成富.基于MATLAB的最小二乘法参数辨识与仿真［J］.微处理机，2005，12（6）：44-46.Shi Xianliang，Wu Chengfu.Rls parameter identification and emulate based on Matlab／Simulink［J］.Microprocessors，2005，12（6）：44-46.（in Chinese）

［11］ Xiong L，Ma H D，Fang H Z，et al.Anomaly detection of spacecraft based on least squares support vector machine［C］／／Prognostics and System Health Management Conference（PHM-Shenzhen），2011.IEEE，2011：1-6.

［12］ Milhano T，Sequeira J，Di Sotto E.Spacecraft parameter identification using S-estimators［C］／／9th International ESA Conference on Guidance，Navigation＆Control Systems，O-porto，Portugal.June 2014.

［13］ Gunnarsson S.Combining tracking and regularization in recursive least squares identification［C］／／Decision and Control，1996.，Proceedings of the 35th IEEE Conference on.IEEE，1996，3：2551-2552.

Research on Identification of Mass Characteristics of Spacecraft during Capture Brake

HAN Yu1，2，CAO Tao1，2
（1.Shanghai Aerospace Control Technology Institute，Shanghai 201400，China；2.Shanghai Key Laboratory of Aerospace Intelligent Control Technology，Shanghai 201400，China）

Targeting the on-orbit identification of the system mass，the position of the center of mass and the inertia matrix of a spacecraft，a method was proposed in this paper for identifying the unknown parameters in a system having a set of governing equations describing its behavior that cannot be put into regression form with the unknown parameters linearly represented.In this method，the vector of the unknown parameters was segmented into a plurality of groups where each individual group of the unknown parameters could be isolated linearly by manipulation of the said equations.Multiple concurrent and independent recursive least squares identification of each said group run，

treating other unknown parameters appearing in their regression equation as if they were known perfectly，with said values provided by recursive least squares estimation from the other groups，thereby enabling the use of fast，compact，efficient linear algorithms to solve problems that would otherwise require nonlinear solution approaches.The validity of this method was verified by numerical simulation in identification of the total mass，the position of the center of mass and inertia matrix for a spacecraft with thruster actuation.

capture brake；least square；mass characteristics identification；thruster actuation

V412.4；V448.25

A

1674-5825（2017）06-0724-07

2017-02-09；

2017-09-18

载人航天预先研究项目（060101）

韩宇，女，博士，工程师，研究方向为航天器姿态控制、容错控制。E-mail：tabubu＠126.com

（责任编辑：龙晋伟）