航天器的高精度姿态控制离不开精确的姿态确定，精确的姿态确定离不开精确的安装误差标定，如通过姿态敏感器(星敏、地敏、陀螺、磁强计等)进行姿态确定[7]。对于复杂航天器，航天器上往往安装可以两轴旋转的伺服机构，而同时为了扩大视域，把成像器安装在伺服机构上，通过伺服机构的指向控制来增大成像器的成像范围；在一些需要伺服机构进行高精度跟踪的任务中，决定跟踪精度的因素有跟踪控制算法、伺服机构硬件性能、伺服机构的安装误差、成像器的安装误差等。由于运载发射过程对结构安装的影响，伺服机构的安装误差、成像器的安装误差的确定需要运用在轨标定技术来精确计算。很多情况下，安装器件的安装误差在不同的敏感轴之间会产生耦合效应[8-10]。标定安装误差在一定程度上与姿态确定十分类似，相当于将未知的姿态等效为未知的安装。因此，姿态确定方法与安装误差标定方法互相借鉴。具有明显几何意义的姿态确定算法主要有TRIAD法、QUEST法等[11]，这些方法不需要复杂庞大的数据支撑，仅仅依赖几个矢量轴的确定即可完全确定姿态[12-13]，由于这些方法计算简单，相对而言应用性很高[14-16]。然而，对于超过六维的耦合安装误差标定，目前的安装误差计算方法仅仅能够确定一个误差安装矩阵，并不能完成多维安装误差标定。
多个串联机械机构安装标定的难点在于建立标定方程组时会发现不同的安装之间互相影响，其内部的耦合效应导致无法采用消元的方法将其中某个安装解耦出来，若采用联立方程组的方法进行解算，由于每个安装有3个未知参数，当安装维数为2时，联立方程组将包括多个含6个未知数的非线性方程，只能通过复杂的数值方法解算，并且难以保证解的精度为零。
因此，针对多维安装难以标定的问题，本文从以下几个方面展开：首先，将多维安装标定问题进行具体的描述，通过建立数学方程来分析标定计算面临的问题；然后，介绍本文的多维安装解耦标定方法，包含矩阵方程建立、串联安装解耦、安装参数确定、一般情形推广；最后，通过仿真分析，对无噪声、有噪声两种情况下标定方法的准确性进行验证。
1 问题的提出 航天器上的光学载荷安装标定对提高指向控制精度、获取高精度目标信息具有重要的作用。在航天器本体上安装一个可转动机构，该可转动机构上固连一个光学观测设备。该可转动机构可以按照预期的姿态进行机动来控制光学观测设备达到一定的指向。一般而言，可转动机构在其零位时与航天器本体之间的坐标变换矩阵、可转动机构与光学探测设备之间的坐标变换矩阵需要进行高精度标定，以保证光学观测设备达到更加精确的指向。航天器标定机构示意图如图 1所示。
 |
| 图 1 航天器标定机构示意图 Fig. 1 Schematic diagram of spacecraft calibration mechanism |
| 图选项 |
图 1中，航天器本体与可转动机构零位之间的坐标变换矩阵记作A，为未知量；可转动机构零位与可转动机构之间的坐标变换矩阵记作B，为已知量；可转动机构与光学观测器之间的坐标变换矩阵记作C，为未知量。
在已知航天器本体在赤道惯性坐标系下的姿态(记作矩阵M)、某恒星目标的赤道惯性系下的方向信息(记作矢量Tar)、该恒星目标在光学观测器中的方向信息(记作矢量View)后，按照三维空间的矢量转换可得：
 | (1) |
令：r=M* Tar，式(1)可以简化为
 | (2) |
此时，式(2)中的C、A为求解对象，方程中其他变量通过观测恒星目标、获取光学探测器中的该星点方位得到。
研究式(2)，发现在C、A中间仍然存在一个可变矩阵B。虽然矩阵B为已知量(可以通过已知可转动平台运动角度得出)，但是要准确解出C、A，需要解决它们之间耦合的问题。当展开式(2)会得到一个方程组，在方程组中C、A的总共6个未知数无法彼此独立，同时方程为非线性方程，若试图通过代数直接解方程的方式，求解过程会相当复杂，通过繁琐的数值解算方法或许能够得到接近准确解的结果。
通常，多次观测能够积累充分的观测数据，这些数据在理论上满足方程C*Bi*A*ri=View i，从最小二乘意义上，可以将标定问题做如下描述：
已知n个坐标变换矩阵Bi、n个矢量ri、n个矢量View i, 求出坐标变换矩阵C、A的值，使得
 | (3) |
取极小值。式中：λi为加权系数。
解方程式(3)使得J取极小值是困难的，主要难度在于J表达式中存在众多的非线性项，一般情况采用数值计算的方式进行，而数值计算方式是繁琐的、有时也难以收敛，因此并不适用于在轨标定。
本文试图绕过复杂的、不适用在轨标定的数值计算方式，建立一种简化的、准确的计算方式，该计算方式基于观测路径的规划设计，结合数学公式推导，灵活运用特征向量的性质，并将2维安装标定推广到n维安装标定。
2 多维安装解耦标定方法 2.1 矩阵方程建立 第1节得到矢量方程C*B*A*r= View后，进一步通过另一组观测数据B、r′、View ′得到另外一组矢量方程，将两组矢量方程一并写出：
 | (4) |
令：
 | (5) |
有以下方程成立：
 | (6) |
即
 | (7) |
注意到[a b c]为单位正交阵，它的逆为它的转置，所以可以得出矩阵方程：
 | (8) |
用D来表示等式右边的矩阵，可以得出以下矩阵方程：
 | (9) |
于是问题转换为如何通过上述矩阵方程得出C、A的值。
2.2 串联安装解耦 在得到式(9)的矩阵方程后，需进行方程解耦。在经过多次观测后可得出3个矩阵方程：
 | (10) |
对式(10)进行矩阵运算得到
 | (11) |
将方程组(11)化为含2个矩阵方程的方程组，将A、C解耦：
 | (12) |
进一步整理后有
 | (13) |
令：
 | (14) |
所以解耦后仅含A的方程组可以写为
 | (15) |
此时，经过解耦运算，得到了仅含A的方程组，问题变为如何根据解耦后仅含A的方程组确定A的值。
2.3 安装参数确定 2.2节得到了仅含A的方程组，该方程组同样属于非线性方程组，耦合性由A与C耦合转变成了A与AT的耦合。本节灵活运用特征向量的性质，通过一定的数学运算，对解耦A、C后的方程组进一步解耦，以得到A的解析表达式。
由于M1、M2、N1和N2均为单位正交阵，所以其一定存在特征值为1的特征向量，该特征向量考虑方向重合与方向相反时有2个，将这些特征向量写为±pM1、±pM2、±pN1和±pN2。
此时要求pM1、pM2线性无关，pN1、pN2线性无关。
考虑特征向量的性质，有式(16)成立：
 | (16) |
对解耦后仅含A的方程组右乘特征向量：
 | (17) |
结合式(16)、式(17)得到
 | (18) |
 | (19) |
所以AT*pM1、AT*pM2分别是N1、N2的特征值为1的特征向量。
此时，在N1、N2的特征向量中，各选择其中一个特征向量pNα、pNβ，如果方程有解，则pM1、pM2的夹角和pNα、pNβ的夹角相同，满足：
 | (20) |
则满足以下方程：
 | (21) |
或
 | (22) |
将pM1、pM2，pNα、pNβ或pM1、pM2，－pNα、－pNβ看成2个不同姿态的刚体上的矢量，同2.1节的处理方式类似，得出A*、A* *。
分别令A=A*和A=A* *，并代入到式(15)中来确定A=A*还是A=A* *。
进而得出：
 | (23) |
此时运用“特征向量”的性质，将耦合的非线性方程组进行了解耦，即得到了解出A、C的方法和公式。该处理方式具备一定的启发性，可以将该思路推广到更多维度的未知矩阵解耦标定中去。
2.4 一般情形推广 针对一般情形，问题可以描述为：
定理??已知以下等式：
 | (24) |
则可通过有限次观测，来确定A1，A2，…，An。式中:Ai未知且固定，Bi已知且可变，r、View已知且可变。
证明??同第2.1节的处理类似，通过两次观测将式(24)转变为矩阵方程：
 | (25) |
令：Pn=An*Bn－1*…*B2*A2*B1*A1，将矩阵方程式(25)等价写为
 | (26) |
按照第2.2节、2.3节的处理方式，可以得出An和Pn-1的值。因此，可以得出结论：
已知Pn的值，可以按照一定的数学运算得出An和Pn-1的值，进一步通过Pn-1得到An-1和Pn-2的值, ……, 最终得出A1的值。于是将n维安装解耦，标定出所有维度的安装矩阵。??证毕
3 仿真分析 仿真分析基于2维安装解耦标定，给定某航天器的初始轨道及飞行姿态、可转动平台零位在航天器上的安装矩阵、可转动平台相对于其零位的坐标转换矩阵、部分恒星J2000系下的方位数据、光学观测器中的星点方位数据作为仿真算法的数据输入，另外考虑光学观测器的噪声影响，对观测噪声对最终标定结果的影响定量分析。
仿真流程如图 2所示。
 |
| 图 2 仿真流程 Fig. 2 Simulation process |
| 图选项 |
航天器的初始轨道在J2000系下的位置坐标为(10，20，-6 778) km, 速度矢量为(7，1，0.5) km/s。航天器姿态保持与轨道坐标系重合，航天器运动轨迹如图 3所示。
 |
| 图 3 航天器三维轨迹 Fig. 3 Three-dimensional trajectories of spacecraft |
| 图选项 |
可转动平台零位在航天器上的安装矩阵对应的四元数为
 |
观测器相对于可转动平台的安装矩阵对应的四元数为
 |
可转动平台相对于其零位的坐标转换矩阵在6次观测中分别为
第1次、第2次：
 |
第3次、第4次：
 |
第5次、第6次：
 |
部分恒星的方位数据如表 1所示，光学观测器中的星点方位数据如表 2所示。
表 1 恒星方位数据(J2000系下单位矢量) Table 1 Stellar orientation data (unit vector in J2000 coordinate system)
| 坐标 | 第1次 | 第2次 | 第3次 | 第4次 | 第5次 | 第6次 |
| x | 0.775 810 6 | 0.416 382 6 | -0.440 612 21 | -0.946 873 74 | -0.384 996 277 0 | 0.372 642 24 |
| y | 0.602 428 5 | 0.170 308 8 | -0.728 514 67 | -0.174 203 86 | 0.768 706 644 40 | 0.885 212 69 |
| z | 0.187 610 4 | 0.893 095 9 | 0.524 525 74 | -0.270 338 9 | -0.510 752 348 5 | -0.278 453 29 |
表选项
表 2 星点方位数据(观测器系下单位矢量) Table 2 Star orientation data (unit vector in observer coordinate system)
| 坐标 | 第1次 | 第2次 | 第3次 | 第4次 | 第5次 | 第6次 |
| x | 0.984 807 75 | 0.766 044 443 | 0.5 | 0.939 692 62 | 0.642 787 60 | 0.866 025 4 |
| y | 0.172 987 39 | 0.604 022 77 | 0.75 | 0.336 824 1 | 0.694 272 044 | 0.482 962 91 |
| z | 0.015 134 43 | 0.219 846 31 | 0.433 012 70 | 0.059 391 174 | 0.323 744 370 | 0.129 409 52 |
表选项
3.1 观测无噪声仿真 在没有任何误差时，按照第2节的标定方法得出的标定计算结果为
A对应四元数为
 |
C对应的四元数为
 |
第1节中的J值为：3.166 183 689 589 54×10-30，几乎为零。结果表明A、C的解算结果qA、qC与真实值qA0、qC0一致，计算精度几乎为零。因此，解耦标定方法是正确的，在没有噪声干扰的情况下，该方法能够完全将未知的安装参数解出，与理论推导结论一致。
3.2 观测有噪声仿真 当存在观测误差，对光学观测器的噪声从0°增加到0.5°，步长取15″, 标定误差随着噪声变化的曲线如图 4所示。
 |
| 图 4 A的标定误差随观测噪声变化 Fig. 4 Variation of calibration error of A with observation noise |
| 图选项 |
图 4表明在观测噪声从0°增加到0.5°时，标定方法解出的未知矩阵A对应的四元数与真实矩阵A对应的四元数之间作差时，其描述2个四元数是否接近的角度差呈现增加的趋势，基本上标定误差的量级随噪声水平的提高同比增加。
图 5表明在观测噪声从0°增加到0.5°时，未知矩阵C的标定误差呈现增加的趋势，基本上标定误差的量级随噪声水平的提高同比增加。
 |
| 图 5 C的标定误差随观测噪声变化 Fig. 5 Variation of calibration error of C with observation noise |
| 图选项 |
图 6给出了最小二乘意义上描述标定准确性的J值变化情况，随着噪声水平的提高，其值也相应的增加，证明噪声水平越大，对算法的精度影响越大。
 |
| 图 6 J的标定误差随观测噪声变化 Fig. 6 Variation of calibration error of J with observation noise |
| 图选项 |
仿真给出了无噪声影响和有噪声影响下标定方法的准确性。仿真结果得出了无噪声影响时标定精度几乎为零的结论，证明了标定方法的正确性。在观测噪声水平提高时，仿真结果表明标定方法的精度会受到观测噪声的影响，基本上影响程度与噪声水平呈现等量级关系。
4 结论 1) 提出了多维安装解耦标定方法，有助于提高指向控制精度、获取高精度目标信息。
2) 针对多个串联机械机构安装标定的问题，建立了多维安装耦合标定的矩阵方程，分析了矩阵方程成立应当满足的必要条件，通过灵活运用“特征向量”的数学特性，经过严谨的数学推导，实现了未知矩阵参数之间的解耦，得出了以成像器观测为基础的在轨标定数学方法，并推广到n维串联安装标定，克服了标定方程组内部的耦合效应。
3) 对解耦标定方法进行了数学仿真分析，建立了仿真计算流程。从无噪声仿真中验证的标定方法的正确性；从有观测噪声仿真中定量分析了噪声对标定方法的影响，得出了影响大小与噪声水平同量级的结论。仿真分析证明了多维安装解耦标定方法的有效性、准确性。
本文的局限性在于由于巧妙使用了特征向量的性质，使用的观测数据需满足本文规定的一些条件如相同B矩阵将使用两次，因此并没有将所有的观测数据引入计算以更大限度地提高计算精度。未来将在引入任意大量数据时最大化提高计算精度方面进行研究。
参考文献
| [1] | LIKINS P W. Spacecraft attitude dynamics and control-a personal perspective on early developments[J]. Journal of Guidance, Control, and Dynamics, 1986, 9(2): 129-134. |
| [2] | 李青.充液挠性系统动力学分析及在航天工程中的应用研究[D].北京: 清华大学, 2010: 50-75. LI Q.Dynamic analysis of liquid-filled flexible systems and its application studies on aerospace engineering[D].Beijing: Tsinghua University, 2010: 50-75(in Chinese). |
| [3] | ZOU A M, KUMAR K D. Adaptive attitude control of spacecraft without velocity measurements using Chebyshev neural network[J]. Acta Astronautica, 2010, 66(5-6): 769-779. DOI:10.1016/j.actaastro.2009.08.020 |
| [4] | PECHEV C, ALEXANDRE N. Feedback-based steering law for control moment gyros[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(3): 848-855. DOI:10.2514/1.27351 |
| [5] | ARDUINI C, BAIOCCO P. Active magnetic damping attitude control for gravity gradient stabilized spacecraft[J]. Journal of Guidance, Control, and Dynamics, 1997, 20(1): 117-122. |
| [6] | 周军. 航天器控制原理[M]. 西安: 西北工业大学出版社, 2000: 62-68. ZHOU J. Principle of spacecraft control[M]. Xi'an: Press of Northwestern Polytechnical University, 2000: 62-68. (in Chinese) |
| [7] | 黄圳圭. 航天器姿态动力学[M]. 长沙: 国防科技大学出版社, 1997: 7-23. HUANG Z G. Spacecraft attitude dynamics[M]. Changsha: Press of National University of Defense Technology, 1997: 7-23. (in Chinese) |
| [8] | 刘俊, 石云波, 李杰. 微惯性技术[M]. 北京: 电子工业出版社, 2005: 15-33. LIU J, SHI Y B, LI J. Micro inertia technology[M]. Beijing: Publishing House of Electronics Industry, 2005: 15-33. (in Chinese) |
| [9] | ZHU R, ZHOU Z Y. Calibration of three-dimensional intergrated sensors for improved system accuracy[J]. Sensors and Actuators, 2006, 127: 340-344. DOI:10.1016/j.sna.2005.12.001 |
| [10] | 刘秉, 闫建国. 一种小型惯性测量单元的精确标定技术[J]. 计算机测量与控制, 2009, 17(11): 2170-2172. LIU B, YAN J G. Precision-calibration technique for inertial measurement units[J]. Computer Measurement and Control, 2009, 17(11): 2170-2172. (in Chinese) |
| [11] | SHUSTER M D, OH S D. Three-axis attitude determination from vector observations[J]. Journal of Guidance, Control, and Dynamics, 1981, 4(1): 70-77. |
| [12] | MARKLEY F L, MORTARI D. Quaternion attitude estimation using vector observations[J]. Journal of the Astronautical Sciences, 2000, 48(2): 359-380. |
| [13] | LERNER G M. Three-axis attitude determination[J]. Spacecraft Attitude Determination and Control, 1978, 73: 420-428. |
| [14] | RENATO B, FABIO C. A robust optimization approach for magnetic spacecraft attitude stabilization[J]. Journal of Optimization Theory and Applications, 2017, 173(3): 994-1012. |
| [15] | CHAK Y C, VARATHARAJOO R, RAZOUMNY Y. Disturbance observer-based fuzzy control for flexible spacecraft combined attitude & sun tracking system[J]. Acta Astronautica, 2017, 133: 302-310. DOI:10.1016/j.actaastro.2016.12.028 |
| [16] | STEVENSON D, SCHAUB H. Electrostatic spacecraft rate and attitude control-experimental results and performance considerations[J]. Acta Astronautica, 2016, 119: 22-33. DOI:10.1016/j.actaastro.2015.10.023 |
