CN108181916B - 小卫星相对姿态的控制方法及装置 - Google Patents

小卫星相对姿态的控制方法及装置 Download PDF

Info

Publication number
CN108181916B
CN108181916B CN201711478457.3A CN201711478457A CN108181916B CN 108181916 B CN108181916 B CN 108181916B CN 201711478457 A CN201711478457 A CN 201711478457A CN 108181916 B CN108181916 B CN 108181916B
Authority
CN
China
Prior art keywords
satellite
control
attitude
matrix
relative
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201711478457.3A
Other languages
English (en)
Other versions
CN108181916A (zh
Inventor
孟子阳
郝运
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Original Assignee
Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tsinghua University filed Critical Tsinghua University
Priority to CN201711478457.3A priority Critical patent/CN108181916B/zh
Publication of CN108181916A publication Critical patent/CN108181916A/zh
Application granted granted Critical
Publication of CN108181916B publication Critical patent/CN108181916B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/244—Spacecraft control systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明公开了一种小卫星相对姿态的控制方法及装置。其中，方法包括：对单颗卫星进行姿态确定与控制；得到干扰力矩；获取单颗卫星的姿态估计值，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；根据共享的姿态信息得到双卫星之间的相对姿态；通过相对姿态动力学模型和相对姿态控制目标计算控制力矩控制率，并控制任一卫星的执行器完成任一卫星对于另一卫星的相对姿态机动控制。该方法可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。

Description

小卫星相对姿态的控制方法及装置
技术领域
本发明涉及航天器姿态测量和控制技术领域，特别涉及一种小卫星相对姿态的控制方法及装置。
背景技术
传统大卫星重量较大，研制周期较长，随着时间的推移，各种缺点已经导致大卫星的研究进程被制约发展。微小卫星以重量轻、低成本、周期短等优点，逐渐进入研究者的视野并大力发展。但是，微小卫星的体积小导致其功能较为单一，因此可以采用多个卫星编队来实现更多的功能和任务，从而弥补微小卫星的缺点，并且实现比大卫星更丰富的功能。其中，对两颗卫星进行相对姿态确定和控制，是卫星编队的基础，也是让编队卫星能够达到预定状态正常工作的重要保证。
目前，若要充分发挥编队卫星的优势，主要依赖于卫星间的相对姿态的确定和控制精度，而国内外对于卫星之间相对姿态确定主要利用绝对测量的结果进行差分而获得星间相对状态，包括基于GPS(Global Positioning System，全球定位系统)的相对姿态确定技术、基于视觉的相对姿态确定技术、基于激光测量的相对姿态确定技术等。但是，美国对GPS的控制权，导致我国军事上的应用受到影响；而视觉导航技术有效测量范围有限，无法进行中远程的卫星相对姿态确定；并且激光发生器的体积和功耗较大，不适宜在微小卫星上进行使用。同时，运用利用绝对测量的结果进行差分而获得星间相对状态的方法的确定精度无法保证，普遍精度不高。而对于存在星间相对测量的相对姿态确定技术，会进一步增加卫星的部件和功能，使得卫星的任务变得更多更繁琐，系统变得更为复杂。因此，提高卫星相对姿态确定的精度，是编队卫星发展的重要方向。
相关技术中，卫星间的相对姿态控制方式目前主要包括主从控制方式和行为控制方式两种。主从控制方式为编队中有一颗主星，其余卫星为从星，对主星实施单独的姿态控制之后，使得从星跟踪主星的姿态运动，最终使得从星和主星都达到期望的控制目标，但是一旦主星姿态失控，就会使得整个编队姿态失控；而行为控制方式并无主星、从星之分，而是由编队的整体期望状态、自身状态和相邻卫星的姿态来确定的，并且可以达到高精度的姿态控制，但是控制方法较为复杂，目前还停留于理论研究方面。
发明内容
本发明旨在至少在一定程度上解决相关技术中的技术问题之一。
为此，本发明的一个目的在于提出一种小卫星相对姿态的控制方法，该方法可以实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
本发明的另一个目的在于提出一种小卫星相对姿态的控制装置。
为达到上述目的，本发明一方面实施例提出了一种小卫星相对姿态的控制方法，包括以下步骤：对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；根据所述共享的姿态信息得到所述双卫星之间的相对姿态；通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制。
本发明实施例的小卫星相对姿态的控制方法，可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，并利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
另外，根据本发明上述实施例的小卫星相对姿态的控制方法还可以具有以下附加的技术特征：
进一步地，在本发明的一个实施例中，所述获取所述单颗卫星的卫星姿态的估计值，进一步包括：
Figure BDA0001533306750000021
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，
Figure BDA0001533306750000022
为系统的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为系统的状态变量，B(t)为输入矩阵，u(t)为系统输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量
Figure BDA0001533306750000023
其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为
Figure BDA0001533306750000024
当t＝(k-1)T时，
Figure BDA0001533306750000025
为初始值，当t＝kT时，进一步得到预测值
Figure BDA0001533306750000031
其中，T为采样周期值，t为时间变量。
将所述微分方程线性化之后得到状态误差向量
Figure BDA0001533306750000032
利用数学关系得到线性干扰方程
Figure BDA0001533306750000033
即
Figure BDA0001533306750000034
其中，
Figure BDA0001533306750000035
Δx为状态误差向量,
Figure BDA0001533306750000036
为误差四元数的矢量部分，b为陀螺常值漂移，
Figure BDA0001533306750000037
为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，
Figure BDA0001533306750000038
为误差四元数的矢量部分随时间的导数，
Figure BDA0001533306750000039
为陀螺常值漂移随时间的导数，
Figure BDA00015333067500000310
为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵，并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵。
对线性干扰方程和量测方程进行基本解阵离散化，即将
Figure BDA00015333067500000311
离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且
Figure BDA00015333067500000312
I为单位矩阵，T为采样周期，
Figure BDA00015333067500000313
为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且
Figure BDA00015333067500000314
Vk为量测噪声矩阵；
预测方差，所述预测方差的公式为：
Figure BDA00015333067500000315
其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，
Figure BDA00015333067500000316
为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：
Figure BDA0001533306750000041
其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，
Figure BDA0001533306750000042
为量测矩阵的转置；
计算方差，所述计算方差的公式为：
Figure BDA0001533306750000043
其中，I为单位矩阵，Pk为k时刻的方差，
Figure BDA0001533306750000044
为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数
Figure BDA0001533306750000045
通过方程
Figure BDA0001533306750000046
求得
Figure BDA0001533306750000047
那么所述估计值为：
Figure BDA0001533306750000048
其中，δqk为误差姿态四元数，
Figure BDA0001533306750000049
为误差姿态四元数的矢量部分，
Figure BDA00015333067500000410
为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，
Figure BDA00015333067500000411
为k时刻状态误差向量的估计值。
进一步地，在本发明的一个实施例中，所述根据所述共享的姿态信息得到所述双卫星之间的相对姿态，进一步包括：
根据两颗卫星的所述相对姿态四元数
Figure BDA00015333067500000412
求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，
Figure BDA00015333067500000413
为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
Figure BDA00015333067500000414
其中，
Figure BDA00015333067500000415
为q的矢量部分，
Figure BDA00015333067500000416
A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
进一步地，在本发明的一个实施例中，所述相对姿态控制目标计算控制力矩控制率，进一步包括：
Figure BDA00015333067500000417
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
Figure BDA0001533306750000051
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，
Figure BDA0001533306750000052
为第二颗卫星的角速度随时间的导数，
Figure BDA0001533306750000053
为相对角速度的叉乘矩阵，
Figure BDA0001533306750000054
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，所述控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制，所述控制过程分为三个阶段，进一步包括：
第一阶段，进入消旋模式。磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为
Figure BDA0001533306750000055
其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；
第二阶段，进入单星姿态稳定模式。启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为
Figure BDA0001533306750000056
在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照
Figure BDA0001533306750000057
的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，
Figure BDA0001533306750000058
为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量；
第三阶段，进入高精度相对姿态控制模式。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照
Figure BDA0001533306750000059
的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，
Figure BDA00015333067500000510
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
为达到上述目的，本发明另一方面实施例提出了一种小卫星相对姿态的控制装置，包括：采集模块，用于对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；第一获取模块，用于获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；第二获取模块，用于获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；计算模块，用于根据所述共享的姿态信息得到所述双卫星之间的相对姿态；控制模块，用于通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制。
本发明实施例的小卫星相对姿态的控制装置，可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
另外，根据本发明上述实施例的小卫星相对姿态的控制装置还可以具有以下附加的技术特征：
进一步地，在本发明的一个实施例中，所述第二获取模块还用于对于具有如下形式的系统：
Figure BDA0001533306750000061
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，
Figure BDA0001533306750000062
为系统的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为系统的状态变量，B(t)为输入矩阵，u(t)为系统输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量
Figure BDA0001533306750000063
其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为
Figure BDA0001533306750000064
当t＝(k-1)T时，
Figure BDA0001533306750000065
为初始值，当t＝kT时，进一步得到预测值
Figure BDA0001533306750000066
其中，T为采样周期值，t为时间变量；
将所述微分方程线性化之后得到状态误差向量
Figure BDA0001533306750000067
利用数学关系得到线性干扰方程
Figure BDA0001533306750000068
即
Figure BDA0001533306750000069
其中，
Figure BDA0001533306750000071
Δx为状态误差向量,
Figure BDA0001533306750000072
为误差四元数的矢量部分，b为陀螺常值漂移，
Figure BDA0001533306750000073
为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，
Figure BDA0001533306750000074
为误差四元数的矢量部分随时间的导数，
Figure BDA0001533306750000075
为陀螺常值漂移随时间的导数，
Figure BDA0001533306750000076
为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵。并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵；
对线性干扰方程和量测方程进行基本解阵离散化，即将
Figure BDA0001533306750000077
离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且
Figure BDA0001533306750000078
I为单位矩阵，T为采样周期，
Figure BDA0001533306750000079
为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且
Figure BDA00015333067500000710
Vk为量测噪声矩阵；
预测方差，所述预测方差的公式为：
Figure BDA00015333067500000711
其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，
Figure BDA00015333067500000712
为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：
Figure BDA00015333067500000713
其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，
Figure BDA00015333067500000714
为量测矩阵的转置；
计算方差，所述计算方差的公式为：
Figure BDA00015333067500000715
其中，I为单位矩阵，Pk为k时刻的方差，
Figure BDA00015333067500000716
为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数
Figure BDA0001533306750000081
通过方程
Figure BDA0001533306750000082
求得
Figure BDA0001533306750000083
那么所述估计值为：
Figure BDA0001533306750000084
其中，δqk为误差姿态四元数，
Figure BDA0001533306750000085
为误差姿态四元数的矢量部分，
Figure BDA0001533306750000086
为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，
Figure BDA0001533306750000087
为k时刻状态误差向量的估计值。
进一步地，在本发明的一个实施例中，所述计算模块还用于根据两颗卫星的所述相对姿态四元数
Figure BDA0001533306750000088
求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，
Figure BDA0001533306750000089
为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
Figure BDA00015333067500000810
其中，
Figure BDA00015333067500000811
为q的矢量部分，
Figure BDA00015333067500000812
A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
进一步地，在本发明的一个实施例中，所述控制模块在采用PD控制可以得到所需控制力矩的控制率为：
Figure BDA00015333067500000813
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
Figure BDA00015333067500000814
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，
Figure BDA00015333067500000815
为第二颗卫星的角速度随时间的导数，
Figure BDA00015333067500000816
为相对角速度的叉乘矩阵，
Figure BDA00015333067500000817
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，所述控制模块进一步用于控制磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为
Figure BDA0001533306750000091
其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；以及启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为
Figure BDA0001533306750000092
在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照
Figure BDA0001533306750000093
的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，
Figure BDA0001533306750000094
为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量，以及改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照
Figure BDA0001533306750000095
的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，
Figure BDA0001533306750000096
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解，其中：
图1为根据本发明一个实施例的小卫星相对姿态的控制方法的流程图；
图2为根据本发明一个实施例的单颗卫星姿态确定和控制系统的总体方案的示意图；
图3为根据本发明一个实施例的两颗卫星相对姿态确定和控制系统的总体方案的示意图；
图4为根据本发明一个实施例的单颗卫星消旋阶段控制模式的功能示意图；
图5为根据本发明一个实施例的单颗卫星三轴稳定后的角速度波动曲线的示意图；
图6为根据本发明一个实施例的单颗卫星三轴稳定后角速度的估计值和真值的差的示意图；
图7为根据本发明一个实施例的单颗卫星三轴稳定后的欧拉角波动曲线的示意图；
图8为根据本发明一个实施例的单颗卫星三轴稳定后欧拉角的估计值和真值的差的示意图；
图9为根据本发明一个实施例的相对姿态控制后，两颗卫星相对角速度的波动曲线的示意图；
图10为根据本发明一个实施例的相对姿态控制后，两颗卫星相对姿态角的波动曲线的示意图；
图11为根据本发明一个实施例的相对姿态控制后，两颗卫星相对姿态角和目标相对姿态角之间的误差波动曲线的示意图；
图12为根据本发明一个实施例的小卫星相对姿态的控制装置的结构示意图。
具体实施方式
下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
下面参照附图描述根据本发明实施例提出的小卫星相对姿态的控制方法及装置，首先将参照附图描述根据本发明实施例提出的小卫星相对姿态的控制方法。
图1是本发明一个实施例的小卫星相对姿态的控制方法的流程图。
如图1所示，该小卫星相对姿态的控制方法包括以下步骤：
在步骤S101中，对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数。
其中，可以设置轨道模型为近地卫星轨道的SGP4轨道模型；卫星参数包括卫星质量M、卫星体积V、卫星初始姿态四元数q0、卫星初始角速度ω0；卫星敏感器参数包括陀螺常值漂移b、陀螺量测噪声η1；卫星执行器参数包括磁力矩器单轴最大磁矩Mm、反作用轮转动惯量Jw、反作用轮最大输出角动量h；干扰力矩模型参数包括地球引力常数μ、卫星磁矩Mr、气动力矩Ta、太阳辐射力矩Ts。
举例而言，本发明实施例可以设置卫星质量M为23.3kg，体积为324mm×360mm×449mm，设置陀螺常值漂移为3°/h、陀螺随机游走噪声标准差为4.2543×10-9rad/s2，星敏感器指向精度为3″，设置反作用轮转动惯量矩阵为diag[10-4 10-4 10-4]、反作用轮最大输出角动量为0.1Nms，设置地球引力常数为398600.5km3/s2、卫星磁矩为0.15A·m2、气动力矩为1.514×10-7N·m、太阳辐射力矩为2.746×10-8N·m，设置两颗卫星初态三轴角速度为[0.0698rad/s 0.0698rad/s 0.0698rad/s]T，初态姿态四元数为[-0.13261 0.388210.18852 0.8923]T(即欧拉角为[26.1793° -5.18° 48.2297°]T)。
在步骤S102中，获取单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩。
本发明实施例利用四元数表示的卫星姿态运动学方程为
Figure BDA0001533306750000111
对于使用反作用轮作为执行器的动力学方程为
Figure BDA0001533306750000112
其中，
Figure BDA0001533306750000113
I为卫星的转动惯量，h为反作用轮相对卫星的角动量，T为卫星所受的总力矩。单颗卫星所受的干扰力矩包括四种，重力梯度力矩利用
Figure BDA0001533306750000114
求得，其中，μ为地球引力常数，r为轨道半径，已在初始条件设置中给出，A为轨道系和本体系之间的姿态矩阵；地磁力矩利用Tm＝Mr×B，其中Mr为卫星磁矩，B为卫星所在处地磁场的磁感应强度，并且气动力矩和太阳辐射力矩均采用最大值代替。
在步骤S103中，获取单颗卫星的卫星姿态的估计值，以将单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息。
本发明实施例可以利用卫星敏感器获取姿态信息，并通过卡尔曼滤波的方式进行处理，设计控制器，通过执行器控制卫星姿态。其中，如图2所示，卫星敏感器可以包括高精度星敏感器、陀螺仪和三轴磁强计等，需要说明的是，卫星敏感器也可以称作姿态敏感器。例如，如图3所示，本发明实施例可以将单颗卫星控制至三轴稳定状态，并通过星间链路使两颗卫星(如卫星1和卫星2)共享另一颗卫星的姿态信息。
具体而言，由星敏感器的量测模型
Figure BDA0001533306750000115
得到目前的姿态四元数，其中qv为误差四元数；由陀螺仪的量测模型u＝ω+b+η1得到目前的角速度，其中，输出u为卫星的真实角速度，ω为陀螺的测量输出，b为陀螺的常值漂移，η1为陀螺的量测噪声，而陀螺的量测噪声可以认为是均值为零的白噪声。
进一步地，在本发明的一个实施例中，获取单颗卫星的卫星姿态的估计值，进一步包括：
本发明实施例所研究的系统具有如下形式：
Figure BDA0001533306750000121
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，
Figure BDA00015333067500001215
为系统的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为系统的状态变量，B(t)为输入矩阵，u(t)为系统输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量
Figure BDA0001533306750000122
其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为
Figure BDA0001533306750000123
当t＝(k-1)T时，
Figure BDA0001533306750000124
为初始值，当t＝kT时，进一步得到预测值
Figure BDA0001533306750000125
其中，T为采样周期值，t为时间变量。
将所述微分方程线性化之后得到状态误差向量
Figure BDA0001533306750000126
利用数学关系得到线性干扰方程
Figure BDA0001533306750000127
即
Figure BDA0001533306750000128
其中，
Figure BDA0001533306750000129
Δx为状态误差向量,
Figure BDA00015333067500001210
为误差四元数的矢量部分，b为陀螺常值漂移，
Figure BDA00015333067500001211
为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，
Figure BDA00015333067500001212
为误差四元数的矢量部分随时间的导数，
Figure BDA00015333067500001213
为陀螺常值漂移随时间的导数，
Figure BDA00015333067500001214
为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵。并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵。
对线性干扰方程和量测方程进行基本解阵离散化，即将
Figure BDA0001533306750000131
离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且
Figure BDA0001533306750000132
I为单位矩阵，T为采样周期，
Figure BDA0001533306750000133
为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且
Figure BDA0001533306750000134
Vk为量测噪声矩阵。
预测方差，所述预测方差的公式为：
Figure BDA0001533306750000135
其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，
Figure BDA0001533306750000136
为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：
Figure BDA0001533306750000137
其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，
Figure BDA0001533306750000138
为量测矩阵的转置；
计算方差，所述计算方差的公式为：
Figure BDA0001533306750000139
其中，I为单位矩阵，Pk为k时刻的方差，
Figure BDA00015333067500001310
为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数
Figure BDA00015333067500001311
通过方程
Figure BDA00015333067500001312
求得
Figure BDA00015333067500001313
那么所述估计值为：
Figure BDA00015333067500001314
其中，δqk为误差姿态四元数，
Figure BDA00015333067500001315
为误差姿态四元数的矢量部分，
Figure BDA00015333067500001316
为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，
Figure BDA00015333067500001317
为k时刻状态误差向量的估计值。
具体而言，根据卡尔曼滤波的方式，对量测得到的结果进行滤波，包括以下几个步骤：
S301：选定状态变量
Figure BDA00015333067500001318
其中，q为卫星姿态四元数，b为陀螺常值漂移。根据原始微分方程
Figure BDA0001533306750000141
以
Figure BDA0001533306750000142
为初始值(t＝(k-1)T)，求出t＝kT(T为采样周期值)的值，得到“一步预测值”
Figure BDA0001533306750000143
S302：将微分方程线性化之后得到状态误差向量
Figure BDA0001533306750000144
S303：利用数学关系得到线性干扰方程
Figure BDA0001533306750000145
也就是
Figure BDA0001533306750000146
其中，
Figure BDA0001533306750000147
S304：对方差进行预测，即为进行
Figure BDA0001533306750000148
的运算，其中，Qk-1为噪声，
Figure BDA0001533306750000149
S305：利用
Figure BDA00015333067500001410
求得k时刻的增益；
S306：计算方差，由
Figure BDA00015333067500001411
得到方差量；
S307：利用星敏的量测方程求得误差姿态四元数
Figure BDA00015333067500001412
之后即可利用方程
Figure BDA00015333067500001413
求得
Figure BDA00015333067500001414
S308：由
Figure BDA00015333067500001415
求得估计状态。
在步骤S104中，根据共享的姿态信息得到双卫星之间的相对姿态。
进一步地，在本发明的一个实施例中，根据共享的姿态信息得到双卫星之间的相对姿态，进一步包括：根据两颗卫星的所述相对姿态四元数
Figure BDA00015333067500001416
求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，
Figure BDA00015333067500001417
为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
Figure BDA0001533306750000151
其中，
Figure BDA0001533306750000152
为q的矢量部分，
Figure BDA0001533306750000153
A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
可以理解的是，由两颗卫星(如卫星1和卫星2)的本体系相对参考系的姿态四元数，可以知道姿态矩阵为：
Figure BDA0001533306750000154
其中，
Figure BDA0001533306750000155
因此可以得到卫星2的本体系到卫星1的本体系的坐标转换矩阵为：
Figure BDA0001533306750000156
根据ω12＝ω1-A21ω2即可得到两星之间的相对角速度；利用四元数的性质，由
Figure BDA0001533306750000157
可以得到两颗卫星的相对姿态。
在步骤S105中，通过相对姿态动力学模型和相对姿态控制目标计算控制力矩控制率，并控制任一卫星的执行器完成任一卫星对于另一卫星的相对姿态机动控制。
可以理解的是，如图2和图3所示，执行器可以为反作用轮和磁力矩器，需要说明的是，执行器也可以称为执行机构。
进一步地，在本发明的一个实施例中，相对姿态控制目标计算控制力矩控制率，进一步包括：
Figure BDA0001533306750000158
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
Figure BDA0001533306750000159
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，
Figure BDA00015333067500001510
为第二颗卫星的角速度随时间的导数，
Figure BDA00015333067500001511
为相对角速度的叉乘矩阵，
Figure BDA00015333067500001512
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，所述控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制，所述控制过程分为三个阶段，进一步包括：
第一阶段，进入消旋模式。磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为
Figure BDA0001533306750000161
其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；
第二阶段，进入单星姿态稳定模式。启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为
Figure BDA0001533306750000162
在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照
Figure BDA0001533306750000163
的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，
Figure BDA0001533306750000164
为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量。
第三阶段，进入高精度相对姿态控制模式。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照
Figure BDA0001533306750000165
的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，
Figure BDA0001533306750000166
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
可以理解的是，姿态控制过程可以分为三个阶段，第一阶段为对星箭分离之后的卫星进行初始消旋，此阶段下只有磁力矩器进行工作，可以得到期望力矩的控制率为Mc＝-Kdω，并合理选取参数，根据数学关系得到磁力矩器产生磁矩公式为
Figure BDA0001533306750000167
如图4所示；当卫星角速度下降至0.5deg/s以下之后，进入第二阶段，即三轴稳定控制阶段，启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为
Figure BDA0001533306750000168
并且合理选取相应参数kp和kd；但是因为反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此，在动量轮开启后，同样调整磁力矩器，按照
Figure BDA0001533306750000169
的规律产生磁矩，其中K为参数，并对反作用轮产生的多余动量进行卸载，使得反作用轮保持在较低转速下运行。最终控制结果如图5-8所示。第三阶段为改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照
Figure BDA0001533306750000171
的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，
Figure BDA0001533306750000172
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差，ω12为相对角速度，Kd，Kp为控制参数。具体而言，本发明实施例对控制参数合理选取，以期望相对姿态q12＝[300 0]T为控制目标进行控制，最终控制结果如图9、10、11所示。
在本发明的一个具体实施例中，在高精度单颗卫星姿态确定和控制的研究基础上，利用MATLABSimulink仿真平台对两颗卫星的相对姿态确定和控制问题进行了数学仿真，卫星质量M＝23.3kg，在单颗卫星的姿态确定和控制精度满足如表1所示的条件时，两颗卫星的相对姿态确定和控制精度可以达到如表2所示的程度。其中，表1为两颗卫星的相对姿态确定和控制精度参数表，表2为两颗卫星的相对姿态确定和控制精度结果表。
表1
Figure BDA0001533306750000173
表2
Figure BDA0001533306750000174
需要说明的是，本发明实施例的单颗卫星的姿态稳定度要优于±6×10-4deg/s，角速度的确定精度要优于0.0003deg/s，同时，姿态角的确定精度优于0.002deg，控制精度优于0.02deg并在上述的基础上，两颗卫星的相对姿态控制精度优于0.03deg。
综上，本发明实施例在高精度单颗卫星姿态确定和控制的研究基础上，建立了两颗卫星的相对姿态确定和控制的闭环系统，从而达到较高的相对姿态确定和控制精度。在本发明实施例的方法中以高精度星敏感器，陀螺仪，三轴磁强计作为姿态敏感器，其中，星敏感器和陀螺仪作为主要姿态测量部件，执行器可以为动量轮和磁力矩器。本发明实施例利用各敏感器和执行器对单星姿态进行控制，同时在两星之间建立信息通道，将两颗卫星的姿态信息通信传输，并利用相对姿态动力学方程设计控制器，最终利用执行器件完成两星之间的相对姿态控制。
根据本发明实施例提出的小卫星相对姿态的控制方法，可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
其次参照附图描述根据本发明实施例提出的小卫星相对姿态的控制装置。
图12是本发明一个实施例的小卫星相对姿态的控制装置的结构示意图。
如图12所示，该小卫星相对姿态的控制装置10包括：采集模块100、第一获取模块200、第二获取模块300、计算模块400和控制模块500。
其中，采集模块100用于对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数。第一获取模块200用于获取单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩。第二获取模块300用于获取单颗卫星的卫星姿态的估计值，以将单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息。计算模块400用于根据共享的姿态信息得到双卫星之间的相对姿态。控制模块500用于通过相对姿态动力学模型和相对姿态控制目标计算控制力矩控制率，并控制任一卫星的执行器完成任一卫星对于另一卫星的相对姿态机动控制。本发明实施例的装置10可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
进一步地，在本发明的一个实施例中，第二获取模块300还用于对于我们所研究的具有如下形式的系统：
Figure BDA0001533306750000181
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，
Figure BDA0001533306750000182
为系统的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为系统的状态变量，B(t)为输入矩阵，u(t)为系统输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量
Figure BDA0001533306750000183
其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为
Figure BDA0001533306750000184
当t＝(k-1)T时，
Figure BDA0001533306750000185
为初始值，当t＝kT时，进一步得到预测值
Figure BDA0001533306750000191
其中，T为采样周期值，t为时间变量。
将所述微分方程线性化之后得到状态误差向量
Figure BDA0001533306750000192
利用数学关系得到线性干扰方程
Figure BDA0001533306750000193
即
Figure BDA0001533306750000194
其中，
Figure BDA0001533306750000195
Δx为状态误差向量,
Figure BDA0001533306750000196
为误差四元数的矢量部分，b为陀螺常值漂移，
Figure BDA0001533306750000197
为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，
Figure BDA0001533306750000198
为误差四元数的矢量部分随时间的导数，
Figure BDA0001533306750000199
为陀螺常值漂移随时间的导数，
Figure BDA00015333067500001910
为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵。并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵。
对线性干扰方程和量测方程进行基本解阵离散化，即将
Figure BDA00015333067500001911
离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且
Figure BDA00015333067500001912
I为单位矩阵，T为采样周期，
Figure BDA00015333067500001913
为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且
Figure BDA00015333067500001914
Vk为量测噪声矩阵。
预测方差，所述预测方差的公式为：
Figure BDA00015333067500001915
其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，
Figure BDA00015333067500001916
为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：
Figure BDA0001533306750000201
其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，
Figure BDA0001533306750000202
为量测矩阵的转置；
计算方差，所述计算方差的公式为：
Figure BDA0001533306750000203
其中，I为单位矩阵，Pk为k时刻的方差，
Figure BDA0001533306750000204
为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数
Figure BDA0001533306750000205
通过方程
Figure BDA0001533306750000206
求得
Figure BDA0001533306750000207
那么所述估计值为：
Figure BDA0001533306750000208
其中，δqk为误差姿态四元数，
Figure BDA0001533306750000209
为误差姿态四元数的矢量部分，
Figure BDA00015333067500002010
为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，
Figure BDA00015333067500002011
为k时刻状态误差向量的估计值。
进一步地，在本发明的一个实施例中，计算模块400还用于根据两颗卫星的相对姿态四元数
Figure BDA00015333067500002012
求得A21；相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，
Figure BDA00015333067500002013
为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
Figure BDA00015333067500002014
其中，
Figure BDA00015333067500002015
为q的矢量部分，
Figure BDA00015333067500002016
A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
进一步地，在本发明的一个实施例中，控制模块500在采用PD控制可以得到所需控制力矩的控制率为：
Figure BDA00015333067500002017
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
Figure BDA0001533306750000211
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，
Figure BDA0001533306750000212
为第二颗卫星的角速度随时间的导数，
Figure BDA0001533306750000213
为相对角速度的叉乘矩阵，
Figure BDA0001533306750000214
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，控制模块500进一步用于控制磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为
Figure BDA0001533306750000215
其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；以及启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为
Figure BDA0001533306750000216
在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照
Figure BDA0001533306750000217
的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，
Figure BDA0001533306750000218
为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照
Figure BDA0001533306750000219
的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，
Figure BDA00015333067500002110
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
需要说明的是，前述对小卫星相对姿态的控制方法实施例的解释说明也适用于该实施例的小卫星相对姿态的控制装置，此处不再赘述。
根据本发明实施例提出的小卫星相对姿态的控制装置，可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (8)

1.一种小卫星相对姿态的控制方法，其特征在于，包括以下步骤：
对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；
获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；
获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；
根据所述共享的姿态信息得到所述双卫星之间的相对姿态；以及
通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制，其中，所述控制所述任一卫星的执行器完成其对于所述另一卫星的相对姿态机动控制，所述控制过程分为三个阶段，进一步包括：
第一阶段，进入消旋模式，磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为
Figure FDA0002221046360000011
其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；
第二阶段，进入单星姿态稳定模式，启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为
Figure FDA0002221046360000012
在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和，因此调整磁力矩器，按照
Figure FDA0002221046360000013
的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，
Figure FDA0002221046360000014
为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量；
第三阶段，进入高精度相对姿态控制模式，改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照
Figure FDA0002221046360000015
的控制率产生控制力矩，达到相对控制目标；其中，T为相对姿态控制力矩，
Figure FDA0002221046360000016
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
2.根据权利要求1所述的小卫星相对姿态的控制方法，其特征在于，所述获取所述单颗卫星的卫星姿态的估计值，进一步包括：
Figure FDA0002221046360000021
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，
Figure FDA0002221046360000022
为系统的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为系统的状态变量，B(t)为输入矩阵，u(t)为系统输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量
Figure FDA0002221046360000023
其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为
Figure FDA0002221046360000024
当t＝(k-1)Ta时，
Figure FDA0002221046360000025
为初始值，当t＝kTa时，进一步得到预测值
Figure FDA0002221046360000026
其中，Ta为采样周期值，t为时间变量；
将所述微分方程线性化之后得到状态误差向量
Figure FDA0002221046360000027
利用数学关系得到线性干扰方程
Figure FDA0002221046360000028
即
Figure FDA0002221046360000029
其中，
Figure FDA00022210463600000210
Δx为状态误差向量，
Figure FDA00022210463600000211
为误差四元数的矢量部分，b为陀螺常值漂移，
Figure FDA00022210463600000212
为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，
Figure FDA00022210463600000213
为误差四元数的矢量部分随时间的导数，
Figure FDA00022210463600000214
为陀螺常值漂移随时间的导数，
Figure FDA00022210463600000215
为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵，并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵；
对线性干扰方程和量测方程进行基本解阵离散化，即将
Figure FDA00022210463600000216
离散化，得到离散化后的方程为：
Figure FDA0002221046360000031
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且
Figure FDA0002221046360000032
I为单位矩阵，Ta为采样周期，
Figure FDA0002221046360000033
为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且
Figure FDA0002221046360000034
Vk为量测噪声矩阵；
预测方差，所述预测方差的公式为：
Figure FDA0002221046360000035
其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，
Figure FDA0002221046360000036
为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：
Figure FDA0002221046360000037
其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，
Figure FDA0002221046360000038
为量测矩阵的转置；
计算方差，所述计算方差的公式为：
Figure FDA0002221046360000039
其中，I为单位矩阵，Pk为k时刻的方差，
Figure FDA00022210463600000310
为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数
Figure FDA00022210463600000311
通过方程
Figure FDA00022210463600000312
求得
Figure FDA00022210463600000313
那么所述估计值为：
Figure FDA00022210463600000314
其中，δqk为误差姿态四元数，
Figure FDA00022210463600000315
为误差姿态四元数的矢量部分，
Figure FDA00022210463600000316
为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，
Figure FDA00022210463600000317
为k时刻状态误差向量的估计值。
3.根据权利要求2所述的小卫星相对姿态的控制方法，其特征在于，根据所述共享的姿态信息得到所述双卫星之间的相对姿态，进一步包括：
根据两颗卫星的所述相对姿态四元数
Figure FDA00022210463600000318
求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，
Figure FDA00022210463600000319
为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
Figure FDA0002221046360000041
其中，
Figure FDA0002221046360000042
为q的矢量部分，
Figure FDA0002221046360000043
A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
4.根据权利要求3所述的小卫星相对姿态的控制方法，其特征在于，所述相对姿态控制目标计算控制力矩控制率，进一步包括：
Figure FDA0002221046360000044
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
Figure FDA0002221046360000045
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，
Figure FDA0002221046360000046
为第二颗卫星的角速度随时间的导数，
Figure FDA0002221046360000047
为相对角速度的叉乘矩阵，
Figure FDA0002221046360000048
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
5.一种小卫星相对姿态的控制装置，其特征在于，包括：
采集模块，用于对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；
第一获取模块，用于获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；
第二获取模块，用于获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；
计算模块，用于根据所述共享的姿态信息得到所述双卫星之间的相对姿态；以及
控制模块，用于通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制，其中，所述控制模块进一步用于控制磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为
Figure FDA0002221046360000051
其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；以及启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为
Figure FDA0002221046360000052
在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和；因此调整磁力矩器，按照
Figure FDA0002221046360000053
的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，
Figure FDA0002221046360000054
为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量；改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照
Figure FDA0002221046360000055
的控制率产生控制力矩，达到相对控制目标，其中，T为相对姿态控制力矩，
Figure FDA0002221046360000056
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
6.根据权利要求5所述的小卫星相对姿态的控制装置，其特征在于，所述第二获取模块还用于对于具有如下形式的系统：
Figure FDA0002221046360000057
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，
Figure FDA0002221046360000058
为系统的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为系统的状态变量，B(t)为输入矩阵，u(t)为系统输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量
Figure FDA0002221046360000059
其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为
Figure FDA00022210463600000510
当t＝(k-1)Ta时，
Figure FDA00022210463600000511
为初始值，当t＝kTa时，进一步得到预测值
Figure FDA00022210463600000512
其中，Ta为采样周期值，t为时间变量；
将所述微分方程线性化之后得到状态误差向量
Figure FDA0002221046360000061
利用数学关系得到线性干扰方程
Figure FDA0002221046360000062
即
Figure FDA0002221046360000063
其中，
Figure FDA0002221046360000064
Δx为状态误差向量，
Figure FDA0002221046360000065
为误差四元数的矢量部分，b为陀螺常值漂移，
Figure FDA0002221046360000066
为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，
Figure FDA0002221046360000067
为误差四元数的矢量部分随时间的导数，
Figure FDA0002221046360000068
为陀螺常值漂移随时间的导数，
Figure FDA0002221046360000069
为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵；并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵；
对线性干扰方程和量测方程进行基本解阵离散化，即将
Figure FDA00022210463600000610
离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且
Figure FDA00022210463600000611
I为单位矩阵，Ta为采样周期，
Figure FDA00022210463600000612
为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且
Figure FDA00022210463600000613
Vk为量测噪声矩阵；
预测方差，所述预测方差的公式为：
Figure FDA00022210463600000614
其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，
Figure FDA00022210463600000615
为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：
Figure FDA00022210463600000616
其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，
Figure FDA0002221046360000071
为量测矩阵的转置；
计算方差，所述计算方差的公式为：
Figure FDA0002221046360000072
其中，I为单位矩阵，Pk为k时刻的方差，
Figure FDA0002221046360000073
为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数
Figure FDA0002221046360000074
通过方程
Figure FDA0002221046360000075
求得
Figure FDA0002221046360000076
那么所述估计值为：
Figure FDA0002221046360000077
其中，δqk为误差姿态四元数，
Figure FDA0002221046360000078
为误差姿态四元数的矢量部分，
Figure FDA0002221046360000079
为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，
Figure FDA00022210463600000710
为k时刻状态误差向量的估计值。
7.根据权利要求6所述的小卫星相对姿态的控制装置，其特征在于，所述计算模块还用于根据两颗卫星的所述相对姿态四元数
Figure FDA00022210463600000711
求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，
Figure FDA00022210463600000712
为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
Figure FDA00022210463600000713
其中，
Figure FDA00022210463600000714
为q的矢量部分，
Figure FDA00022210463600000715
A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
8.根据权利要7所述的小卫星相对姿态的控制装置，其特征在于，所述控制模块在采用PD控制可以得到所需控制力矩的控制率为：
Figure FDA00022210463600000716
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
Figure FDA00022210463600000717
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，
Figure FDA0002221046360000081
为第二颗卫星的角速度随时间的导数，
Figure FDA0002221046360000082
为相对角速度的叉乘矩阵，
Figure FDA0002221046360000083
为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
CN201711478457.3A 2017-12-29 2017-12-29 小卫星相对姿态的控制方法及装置 Active CN108181916B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711478457.3A CN108181916B (zh) 2017-12-29 2017-12-29 小卫星相对姿态的控制方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711478457.3A CN108181916B (zh) 2017-12-29 2017-12-29 小卫星相对姿态的控制方法及装置

Publications (2)

Publication Number Publication Date
CN108181916A CN108181916A (zh) 2018-06-19
CN108181916B true CN108181916B (zh) 2020-04-24

Family

ID=62549033

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711478457.3A Active CN108181916B (zh) 2017-12-29 2017-12-29 小卫星相对姿态的控制方法及装置

Country Status (1)

Country Link
CN (1) CN108181916B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109143280B (zh) * 2018-10-10 2021-05-07 上海微小卫星工程中心 一种卫星集成状态闭环测试系统以及相应的测试方法
CN110109470B (zh) * 2019-04-09 2021-10-29 西安电子科技大学 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制系统
CN110146082B (zh) * 2019-05-05 2021-03-19 中国人民解放军63921部队 利用测速数据实时估计航天器异常姿态的方法和设备
CN112461511B (zh) * 2020-11-10 2022-03-25 中国科学院长春光学精密机械与物理研究所 浮空平台望远镜指向获取方法、装置、设备及存储介质
CN112611373B (zh) * 2020-12-30 2022-05-20 哈尔滨工业大学 一种近地空间短波红外星敏感器流场气动热效应分析方法
CN113022894B (zh) * 2021-03-08 2022-04-19 航天科工空间工程发展有限公司 一种用于微小卫星的相对姿态确定方法
CN113204250B (zh) * 2021-04-29 2022-03-08 西安电子科技大学 强动态条件下的卫星编队相对位置鲁棒高精度估计方法
CN115771624B (zh) * 2023-02-13 2023-05-26 北京航空航天大学 一种基于强化学习的自适应卫星姿轨控制方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6341249B1 (en) * 1999-02-11 2002-01-22 Guang Qian Xing Autonomous unified on-board orbit and attitude control system for satellites
CN102419597A (zh) * 2011-12-05 2012-04-18 哈尔滨工业大学 一种限定相对姿态的大规模编队航天器姿态一致控制方法
CN103455035A (zh) * 2013-08-26 2013-12-18 北京理工大学 基于反步设计和非线性反馈的pd+姿态控制律设计方法
CN107270909A (zh) * 2017-06-09 2017-10-20 西北工业大学 一种利用双阵列天线进行微小卫星相对姿态确定的方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6341249B1 (en) * 1999-02-11 2002-01-22 Guang Qian Xing Autonomous unified on-board orbit and attitude control system for satellites
CN102419597A (zh) * 2011-12-05 2012-04-18 哈尔滨工业大学 一种限定相对姿态的大规模编队航天器姿态一致控制方法
CN103455035A (zh) * 2013-08-26 2013-12-18 北京理工大学 基于反步设计和非线性反馈的pd+姿态控制律设计方法
CN107270909A (zh) * 2017-06-09 2017-10-20 西北工业大学 一种利用双阵列天线进行微小卫星相对姿态确定的方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Natural formation designof nano/micro spacecraft cluster subject to typical tasks;FANG Yuankun 等;《IEEE》;20171218;全文 *
一种星间相对状态估计方法;方元坤 等;《宇航学报》;20181231;全文 *
基于粒子滤波的三轴稳定卫星姿态确定算法的研究;岳洋;《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》;20111215;全文 *
编队飞行卫星相对姿态确定与控制方法研究;吴云华;《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》;20061215;第I、11-27、45-47页 *

Also Published As

Publication number Publication date
CN108181916A (zh) 2018-06-19

Similar Documents

Publication Publication Date Title
CN108181916B (zh) 小卫星相对姿态的控制方法及装置
US20140231589A1 (en) Gyroless Three-Axis Sun Acquisition Using Sun Sensor and Unscented Kalman Filter
JP3027734B2 (ja) 衛星の位置を機内で自主的に求める方法と装置
CN103017760B (zh) 一种大椭圆轨道火星探测器自主对火定向方法
CN108680186B (zh) 基于重力仪平台的捷联式惯导系统非线性初始对准方法
CN106767797B (zh) 一种基于对偶四元数的惯性/gps组合导航方法
CN110304279B (zh) 一种电推进卫星的质心在轨标定补偿方法
Abdelrahman et al. Sigma-point Kalman filtering for spacecraft attitude and rate estimation using magnetometer measurements
CN108959734B (zh) 一种基于实时递推太阳光压力矩辨识方法及系统
Sugimura et al. Attitude determination and control system for nadir pointing using magnetorquer and magnetometer
Srivastava et al. Attitude determination and control system for a leo debris chaser small satellite
Somov et al. Health checking autonomous attitude control system of Earth-observing miniature satellite in initial orientation modes
Al-Jlailaty et al. Efficient attitude estimators: A tutorial and survey
Stearns et al. Multiple model adaptive estimation of satellite attitude using MEMS gyros
CN110667892B (zh) 基于地磁测量的卫星消旋控制方法
CN107132850B (zh) 基于角速度跟踪的变轨姿态保持控制方法
Okasha et al. Relative motion guidance, navigation and control for autonomous orbital rendezvous
Reijneveld et al. Attitude control system of the Delfi-n3Xt satellite
Post et al. Nanosatellite air bearing tests of fault-tolerant sliding-mode attitude control with Unscented Kalman Filter
Fan et al. Attitude optimization control method of agile optical small satellite for nonparallel ground track imaging
CN112393835B (zh) 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法
CN108871312A (zh) 一种重力梯度仪及星敏感器的联合定姿方法
CN114802818A (zh) 晨昏轨道卫星及其对日姿态计算方法、导引方法
CN108089434A (zh) 一种基于磁强计的皮纳卫星姿态捕获方法
Hajiyev et al. Integration of algebraic method and EKF for attitude determination of small information satellites

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant