CN108181916A - 小卫星相对姿态的控制方法及装置 - Google Patents
小卫星相对姿态的控制方法及装置 Download PDFInfo
- Publication number
- CN108181916A CN108181916A CN201711478457.3A CN201711478457A CN108181916A CN 108181916 A CN108181916 A CN 108181916A CN 201711478457 A CN201711478457 A CN 201711478457A CN 108181916 A CN108181916 A CN 108181916A
- Authority
- CN
- China
- Prior art keywords
- satellite
- control
- attitude
- moment
- matrix
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 32
- 238000005259 measurement Methods 0.000 claims abstract description 70
- 239000011159 matrix material Substances 0.000 claims description 134
- 238000006243 chemical reaction Methods 0.000 claims description 60
- 230000006641 stabilisation Effects 0.000 claims description 20
- 238000011105 stabilization Methods 0.000 claims description 20
- 230000006698 induction Effects 0.000 claims description 13
- 238000005070 sampling Methods 0.000 claims description 13
- 230000007704 transition Effects 0.000 claims description 12
- 230000017105 transposition Effects 0.000 claims description 12
- 101100001674 Emericella variicolor andI gene Proteins 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 6
- 238000005295 random walk Methods 0.000 claims description 6
- 230000009467 reduction Effects 0.000 claims description 6
- 230000006340 racemization Effects 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 4
- 238000004891 communication Methods 0.000 abstract description 7
- 238000005516 engineering process Methods 0.000 description 9
- 238000010586 diagram Methods 0.000 description 8
- 230000015572 biosynthetic process Effects 0.000 description 6
- 238000005755 formation reaction Methods 0.000 description 6
- 108091092878 Microsatellite Proteins 0.000 description 4
- 230000008901 benefit Effects 0.000 description 4
- 239000006185 dispersion Substances 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 230000005855 radiation Effects 0.000 description 3
- 241000208340 Araliaceae Species 0.000 description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 235000008434 ginseng Nutrition 0.000 description 2
- 238000011161 development Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
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的控制权，导致我国军事上的应用受到影响；而视觉导航技术有效测量范围有限，无法进行中远程的卫星相对姿态确定；并且激光发生器的体积和功耗较大，不适宜在微小卫星上进行使用。同时，运用利用绝对测量的结果进行差分而获得星间相对状态的方法的确定精度无法保证，普遍精度不高。而对于存在星间相对测量的相对姿态确定技术，会进一步增加卫星的部件和功能，使得卫星的任务变得更多更繁琐，***变得更为复杂。因此，提高卫星相对姿态确定的精度，是编队卫星发展的重要方向。
相关技术中，卫星间的相对姿态控制方式目前主要包括主从控制方式和行为控制方式两种。主从控制方式为编队中有一颗主星，其余卫星为从星，对主星实施单独的姿态控制之后，使得从星跟踪主星的姿态运动，最终使得从星和主星都达到期望的控制目标，但是一旦主星姿态失控，就会使得整个编队姿态失控；而行为控制方式并无主星、从星之分，而是由编队的整体期望状态、自身状态和相邻卫星的姿态来确定的，并且可以达到高精度的姿态控制，但是控制方法较为复杂，目前还停留于理论研究方面。
发明内容
本发明旨在至少在一定程度上解决相关技术中的技术问题之一。
为此，本发明的一个目的在于提出一种小卫星相对姿态的控制方法，该方法可以实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
本发明的另一个目的在于提出一种小卫星相对姿态的控制装置。
为达到上述目的，本发明一方面实施例提出了一种小卫星相对姿态的控制方法，包括以下步骤：对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；根据所述共享的姿态信息得到所述双卫星之间的相对姿态；通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制。
本发明实施例的小卫星相对姿态的控制方法，可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，并利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
另外，根据本发明上述实施例的小卫星相对姿态的控制方法还可以具有以下附加的技术特征：
进一步地，在本发明的一个实施例中，所述获取所述单颗卫星的卫星姿态的估计值，进一步包括：
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，为***的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为***的状态变量，B(t)为输入矩阵，u(t)为***输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为当t＝(k-1)T时，为初始值，当t＝kT时，进一步得到预测值其中，T为采样周期值，t为时间变量。
将所述微分方程线性化之后得到状态误差向量利用数学关系得到线性干扰方程即其中，Δx为状态误差向量,为误差四元数的矢量部分，b为陀螺常值漂移，为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，为误差四元数的矢量部分随时间的导数，为陀螺常值漂移随时间的导数，为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵，并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵。
对线性干扰方程和量测方程进行基本解阵离散化，即将离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且I为单位矩阵，T为采样周期，为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且Vk为量测噪声矩阵；
预测方差，所述预测方差的公式为：其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，为量测矩阵的转置；
计算方差，所述计算方差的公式为：其中，I为单位矩阵，Pk为k时刻的方差，为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数通过方程求得那么所述估计值为：其中，δqk为误差姿态四元数，为误差姿态四元数的矢量部分，为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，为k时刻状态误差向量的估计值。
进一步地，在本发明的一个实施例中，所述根据所述共享的姿态信息得到所述双卫星之间的相对姿态，进一步包括：
根据两颗卫星的所述相对姿态四元数求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
其中，为q的矢量部分，A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
进一步地，在本发明的一个实施例中，所述相对姿态控制目标计算控制力矩控制率，进一步包括：
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，为第二颗卫星的角速度随时间的导数，为相对角速度的叉乘矩阵，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，所述控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制，所述控制过程分为三个阶段，进一步包括：
第一阶段，进入消旋模式。磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；
第二阶段，进入单星姿态稳定模式。启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量；
第三阶段，进入高精度相对姿态控制模式。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
为达到上述目的，本发明另一方面实施例提出了一种小卫星相对姿态的控制装置，包括：采集模块，用于对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；第一获取模块，用于获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；第二获取模块，用于获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；计算模块，用于根据所述共享的姿态信息得到所述双卫星之间的相对姿态；控制模块，用于通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制。
本发明实施例的小卫星相对姿态的控制装置，可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
另外，根据本发明上述实施例的小卫星相对姿态的控制装置还可以具有以下附加的技术特征：
进一步地，在本发明的一个实施例中，所述第二获取模块还用于对于具有如下形式的***：
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，为***的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为***的状态变量，B(t)为输入矩阵，u(t)为***输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为当t＝(k-1)T时，为初始值，当t＝kT时，进一步得到预测值其中，T为采样周期值，t为时间变量；
将所述微分方程线性化之后得到状态误差向量利用数学关系得到线性干扰方程即其中，Δx为状态误差向量,为误差四元数的矢量部分，b为陀螺常值漂移，为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，为误差四元数的矢量部分随时间的导数，为陀螺常值漂移随时间的导数，为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵。并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵；
对线性干扰方程和量测方程进行基本解阵离散化，即将离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且I为单位矩阵，T为采样周期，为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且Vk为量测噪声矩阵；
预测方差，所述预测方差的公式为：其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，为量测矩阵的转置；
计算方差，所述计算方差的公式为：其中，I为单位矩阵，Pk为k时刻的方差，为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数通过方程求得那么所述估计值为：其中，δqk为误差姿态四元数，为误差姿态四元数的矢量部分，为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，为k时刻状态误差向量的估计值。
进一步地，在本发明的一个实施例中，所述计算模块还用于根据两颗卫星的所述相对姿态四元数求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
其中，为q的矢量部分，A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
进一步地，在本发明的一个实施例中，所述控制模块在采用PD控制可以得到所需控制力矩的控制率为：
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，为第二颗卫星的角速度随时间的导数，为相对角速度的叉乘矩阵，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，所述控制模块进一步用于控制磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；以及启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量，以及改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω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中，获取单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩。
本发明实施例利用四元数表示的卫星姿态运动学方程为对于使用反作用轮作为执行器的动力学方程为其中，I为卫星的转动惯量，h为反作用轮相对卫星的角动量，T为卫星所受的总力矩。单颗卫星所受的干扰力矩包括四种，重力梯度力矩利用求得，其中，μ为地球引力常数，r为轨道半径，已在初始条件设置中给出，A为轨道系和本体系之间的姿态矩阵；地磁力矩利用Tm＝Mr×B，其中Mr为卫星磁矩，B为卫星所在处地磁场的磁感应强度，并且气动力矩和太阳辐射力矩均采用最大值代替。
在步骤S103中，获取单颗卫星的卫星姿态的估计值，以将单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息。
本发明实施例可以利用卫星敏感器获取姿态信息，并通过卡尔曼滤波的方式进行处理，设计控制器，通过执行器控制卫星姿态。其中，如图2所示，卫星敏感器可以包括高精度星敏感器、陀螺仪和三轴磁强计等，需要说明的是，卫星敏感器也可以称作姿态敏感器。例如，如图3所示，本发明实施例可以将单颗卫星控制至三轴稳定状态，并通过星间链路使两颗卫星(如卫星1和卫星2)共享另一颗卫星的姿态信息。
具体而言，由星敏感器的量测模型得到目前的姿态四元数，其中qv为误差四元数；由陀螺仪的量测模型u＝ω+b+η1得到目前的角速度，其中，输出u为卫星的真实角速度，ω为陀螺的测量输出，b为陀螺的常值漂移，η1为陀螺的量测噪声，而陀螺的量测噪声可以认为是均值为零的白噪声。
进一步地，在本发明的一个实施例中，获取单颗卫星的卫星姿态的估计值，进一步包括：
本发明实施例所研究的***具有如下形式：
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，为***的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为***的状态变量，B(t)为输入矩阵，u(t)为***输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为当t＝(k-1)T时，为初始值，当t＝kT时，进一步得到预测值其中，T为采样周期值，t为时间变量。
将所述微分方程线性化之后得到状态误差向量利用数学关系得到线性干扰方程即其中，Δx为状态误差向量,为误差四元数的矢量部分，b为陀螺常值漂移，为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，为误差四元数的矢量部分随时间的导数，为陀螺常值漂移随时间的导数，为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵。并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵。
对线性干扰方程和量测方程进行基本解阵离散化，即将离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且I为单位矩阵，T为采样周期，为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且Vk为量测噪声矩阵。
预测方差，所述预测方差的公式为：其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，为量测矩阵的转置；
计算方差，所述计算方差的公式为：其中，I为单位矩阵，Pk为k时刻的方差，为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数通过方程求得那么所述估计值为：其中，δqk为误差姿态四元数，为误差姿态四元数的矢量部分，为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，为k时刻状态误差向量的估计值。
具体而言，根据卡尔曼滤波的方式，对量测得到的结果进行滤波，包括以下几个步骤：
S301：选定状态变量其中，q为卫星姿态四元数，b为陀螺常值漂移。根据原始微分方程以为初始值(t＝(k-1)T)，求出t＝kT(T为采样周期值)的值，得到“一步预测值”
S302：将微分方程线性化之后得到状态误差向量
S303：利用数学关系得到线性干扰方程也就是其中，
S304：对方差进行预测，即为进行的运算，其中，Qk-1为噪声，
S305：利用求得k时刻的增益；
S306：计算方差，由得到方差量；
S307：利用星敏的量测方程求得误差姿态四元数之后即可利用方程求得
S308：由求得估计状态。
在步骤S104中，根据共享的姿态信息得到双卫星之间的相对姿态。
进一步地，在本发明的一个实施例中，根据共享的姿态信息得到双卫星之间的相对姿态，进一步包括：根据两颗卫星的所述相对姿态四元数求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
其中，为q的矢量部分，A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
可以理解的是，由两颗卫星(如卫星1和卫星2)的本体系相对参考系的姿态四元数，可以知道姿态矩阵为：
其中，因此可以得到卫星2的本体系到卫星1的本体系的坐标转换矩阵为：根据ω12＝ω1-A21ω2即可得到两星之间的相对角速度；利用四元数的性质，由可以得到两颗卫星的相对姿态。
在步骤S105中，通过相对姿态动力学模型和相对姿态控制目标计算控制力矩控制率，并控制任一卫星的执行器完成任一卫星对于另一卫星的相对姿态机动控制。
可以理解的是，如图2和图3所示，执行器可以为反作用轮和磁力矩器，需要说明的是，执行器也可以称为执行机构。
进一步地，在本发明的一个实施例中，相对姿态控制目标计算控制力矩控制率，进一步包括：
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，为第二颗卫星的角速度随时间的导数，为相对角速度的叉乘矩阵，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，所述控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制，所述控制过程分为三个阶段，进一步包括：
第一阶段，进入消旋模式。磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；
第二阶段，进入单星姿态稳定模式。启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量。
第三阶段，进入高精度相对姿态控制模式。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
可以理解的是，姿态控制过程可以分为三个阶段，第一阶段为对星箭分离之后的卫星进行初始消旋，此阶段下只有磁力矩器进行工作，可以得到期望力矩的控制率为Mc＝-Kdω，并合理选取参数，根据数学关系得到磁力矩器产生磁矩公式为如图4所示；当卫星角速度下降至0.5deg/s以下之后，进入第二阶段，即三轴稳定控制阶段，启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为并且合理选取相应参数kp和kd；但是因为反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此，在动量轮开启后，同样调整磁力矩器，按照的规律产生磁矩，其中K为参数，并对反作用轮产生的多余动量进行卸载，使得反作用轮保持在较低转速下运行。最终控制结果如图5-8所示。第三阶段为改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差，ω12为相对角速度，Kd，Kp为控制参数。具体而言，本发明实施例对控制参数合理选取，以期望相对姿态q12＝[30 0 0]T为控制目标进行控制，最终控制结果如图9、10、11所示。
在本发明的一个具体实施例中，在高精度单颗卫星姿态确定和控制的研究基础上，利用MATLABSimulink仿真平台对两颗卫星的相对姿态确定和控制问题进行了数学仿真，卫星质量M＝23.3kg，在单颗卫星的姿态确定和控制精度满足如表1所示的条件时，两颗卫星的相对姿态确定和控制精度可以达到如表2所示的程度。其中，表1为两颗卫星的相对姿态确定和控制精度参数表，表2为两颗卫星的相对姿态确定和控制精度结果表。
表1
表2
需要说明的是，本发明实施例的单颗卫星的姿态稳定度要优于±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还用于对于我们所研究的具有如下形式的***：
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，为***的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为***的状态变量，B(t)为输入矩阵，u(t)为***输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为当t＝(k-1)T时，为初始值，当t＝kT时，进一步得到预测值其中，T为采样周期值，t为时间变量。
将所述微分方程线性化之后得到状态误差向量利用数学关系得到线性干扰方程即其中，Δx为状态误差向量,为误差四元数的矢量部分，b为陀螺常值漂移，为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，为误差四元数的矢量部分随时间的导数，为陀螺常值漂移随时间的导数，为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵。并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵。
对线性干扰方程和量测方程进行基本解阵离散化，即将离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且I为单位矩阵，T为采样周期，为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且Vk为量测噪声矩阵。
预测方差，所述预测方差的公式为：其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，为量测矩阵的转置；
计算方差，所述计算方差的公式为：其中，I为单位矩阵，Pk为k时刻的方差，为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数通过方程求得那么所述估计值为：其中，δqk为误差姿态四元数，为误差姿态四元数的矢量部分，为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，为k时刻状态误差向量的估计值。
进一步地，在本发明的一个实施例中，计算模块400还用于根据两颗卫星的相对姿态四元数求得A21；相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
其中，为q的矢量部分，A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
进一步地，在本发明的一个实施例中，控制模块500在采用PD控制可以得到所需控制力矩的控制率为：
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，为第二颗卫星的角速度随时间的导数，为相对角速度的叉乘矩阵，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
进一步地，在本发明的一个实施例中，控制模块500进一步用于控制磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；以及启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
需要说明的是，前述对小卫星相对姿态的控制方法实施例的解释说明也适用于该实施例的小卫星相对姿态的控制装置，此处不再赘述。
根据本发明实施例提出的小卫星相对姿态的控制装置，可以通过两颗卫星星间链路的通讯使其中一个卫星获得另一颗卫星的量测信息，确定出两颗卫星之间的相对姿态，利用控制器和执行器完成相对姿态控制，从而实现对两颗卫星的相对姿态的高精度控制，进而提高控制的准确性和可靠性，简单易实现。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
Claims (10)
1.一种小卫星相对姿态的控制方法，其特征在于，包括以下步骤：
对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；
获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；
获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；
根据所述共享的姿态信息得到所述双卫星之间的相对姿态；以及
通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制。
2.根据权利要求1所述的小卫星相对姿态的控制方法，其特征在于，所述获取所述单颗卫星的卫星姿态的估计值，进一步包括：
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，为***的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为***的状态变量，B(t)为输入矩阵，u(t)为***输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为当t＝(k-1)T时，为初始值，当t＝kT时，进一步得到预测值其中，T为采样周期值，t为时间变量。
将所述微分方程线性化之后得到状态误差向量利用数学关系得到线性干扰方程即其中，Δx为状态误差向量，为误差四元数的矢量部分，b为陀螺常值漂移，为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，为误差四元数的矢量部分随时间的导数，为陀螺常值漂移随时间的导数，为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵，并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵；
对线性干扰方程和量测方程进行基本解阵离散化，即将离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且I为单位矩阵，T为采样周期，为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且Vk为量测噪声矩阵；
预测方差，所述预测方差的公式为：其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，为量测矩阵的转置；
计算方差，所述计算方差的公式为：其中，I为单位矩阵，Pk为k时刻的方差，为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数通过方程求得那么所述估计值为：其中，δqk为误差姿态四元数，为误差姿态四元数的矢量部分，为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，为k时刻状态误差向量的估计值。
3.根据权利要求1所述的小卫星相对姿态的控制方法，其特征在于，根据所述共享的姿态信息得到所述双卫星之间的相对姿态，进一步包括：
根据两颗卫星的所述相对姿态四元数求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
其中，为q的矢量部分，A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
4.根据权利要求3所述的小卫星相对姿态的控制方法，其特征在于，所述相对姿态控制目标计算控制力矩控制率，进一步包括：
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，为第二颗卫星的角速度随时间的导数，为相对角速度的叉乘矩阵，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
5.根据权利要求4所述的小卫星相对姿态的控制方法，其特征在于，所述控制所述任一卫星的执行器完成其对于所述另一卫星的相对姿态机动控制，所述控制过程分为三个阶段，进一步包括：
第一阶段，进入消旋模式。磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；
第二阶段，进入单星姿态稳定模式。启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和，因此调整磁力矩器，按照的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量。
第三阶段，进入高精度相对姿态控制模式。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照的控制率产生控制力矩，达到相对控制目标。其中，T为相对姿态控制力矩，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
6.一种小卫星相对姿态的控制装置，其特征在于，包括：
采集模块，用于对单颗卫星进行姿态确定与控制，得到轨道模型、卫星参数、卫星敏感器参数、卫星执行器参数、干扰力矩模型参数和卫星初始状态参数；
第一获取模块，用于获取所述单颗卫星的姿态运动学方程和动力学方程，得到干扰力矩；
第二获取模块，用于获取所述单颗卫星的卫星姿态的估计值，以将所述单颗卫星控制至三轴稳定状态，并通过星间链路使双卫星中任一卫星共享另一卫星的姿态信息；
计算模块，用于根据所述共享的姿态信息得到所述双卫星之间的相对姿态；以及
控制模块，用于通过相对姿态动力学模型和所述相对姿态控制目标计算控制力矩控制率，并控制所述任一卫星的执行器完成所述任一卫星对于所述另一卫星的相对姿态机动控制。
7.根据权利要求6所述的小卫星相对姿态的控制装置，其特征在于，所述第二获取模块还用于对于具有如下形式的***：
Z(t)＝h[x(t),t]+ν(t)
其中，t为时间变量，为***的状态变量关于时间的导数，f[x(t),t]和h[x(t),t]为状态变量与时间的向量函数，x(t)为***的状态变量，B(t)为输入矩阵，u(t)为***输入，Z(t)为量测，ν(t)为量测噪声；
选定状态变量其中，x(t)为状态变量，q为卫星的姿态四元数，b为陀螺常值漂移；
微分方程为当t＝(k-1)T时，为初始值，当t＝kT时，进一步得到预测值其中，T为采样周期值，t为时间变量；
将所述微分方程线性化之后得到状态误差向量利用数学关系得到线性干扰方程即其中，Δx为状态误差向量，为误差四元数的矢量部分，b为陀螺常值漂移，为陀螺常值漂移的估计值，Δb为陀螺常值漂移的误差，为误差四元数的矢量部分随时间的导数，为陀螺常值漂移随时间的导数，为角速度的估计值的叉乘矩阵，η2(t)为陀螺的随机游走噪声，η1(t)为陀螺量测噪声，I3为三阶单位矩阵，03×3为三阶零矩阵。并且量测方程也需要进行线性化，结果为ΔZ(t)＝H(t)Δx(t)+ν(t)，其中，ΔZ(t)量测误差向量，H(t)为量测矩阵；
对线性干扰方程和量测方程进行基本解阵离散化，即将离散化，得到离散化后的方程为：
Δxk＝Φk,k-1Δxk-1+Wk-1
ΔZk＝HkδXk+Vk
其中，Δxk为k时刻的状态误差向量，Φk,k-1为状态转移矩阵，且I为单位矩阵，T为采样周期，为k-1时刻的状态量的估计值，Wk-1为状态噪声矩阵，ΔZk为k时刻的量测误差向量，Hk为量测矩阵，且Vk为量测噪声矩阵。
预测方差，所述预测方差的公式为：其中，Pk|k-1为根据k-1时刻预测k时刻的方差，Pk-1为k-1时刻的方差，为Φk,k-1的转置，Qk-1为噪声矩阵；
k时刻的增益为：其中，Kk为k时刻的增益矩阵，Hk为量测矩阵，Rk为噪声矩阵，为量测矩阵的转置；
计算方差，所述计算方差的公式为：其中，I为单位矩阵，Pk为k时刻的方差，为增益矩阵的转置；
利用星敏的量测方程求得误差姿态四元数通过方程求得那么所述估计值为：其中，δqk为误差姿态四元数，为误差姿态四元数的矢量部分，为k时刻估计姿态四元数的逆，qm,k为k时刻星敏输出的姿态四元数，为k时刻状态误差向量的估计值。
8.根据权利要求6所述的小卫星相对姿态的控制装置，其特征在于，所述计算模块还用于根据两颗卫星的所述相对姿态四元数求得A21；
相对角速度为：
ω12＝ω1-A21ω2，
其中，A21为所述两颗卫星本体系之间的坐标转换矩阵，q1为第一颗卫星的姿态四元数，为第二颗卫星的姿态四元数的逆，q12为第二颗卫星相对于第一颗卫星的相对姿态四元数，ω12为第二颗卫星相对于第一颗卫星的角速度，ω1为第一颗卫星的角速度，ω2为第二颗卫星的角速度；
转换矩阵和姿态四元数之间的关系可由下式获取：
其中，为q的矢量部分，A为转换矩阵，q1,q2,q3,q4为四元数q的四个分量。
9.根据权利要8所述的小卫星相对姿态的控制装置，其特征在于，所述控制模块在采用PD控制可以得到所需控制力矩的控制率为：
当目标卫星三轴稳定时，可以认为ω2＝0；相应得到简化的控制力矩为：
其中，T为相对姿态控制力矩，A21为所述两颗卫星本体系之间的坐标转换矩阵，ω12为相对角速度，ω2为第二颗卫星的角速度，为第二颗卫星的角速度随时间的导数，为相对角速度的叉乘矩阵，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，I1为第一颗卫星的转动惯量，Kd，Kp为控制参数。
10.根据权利要9所述的小卫星相对姿态的控制装置，其特征在于，所述控制模块进一步用于控制磁力矩器进行工作，得到期望力矩的控制率为Mc＝-Kmω，根据数学关系得到磁力矩器产生磁矩公式为其中，Km为控制参数，ω为卫星角速度，B为卫星位置的磁感应强度，Mc为期望控制力矩，M为磁矩；以及启用反作用轮，采用转速模式对反作用轮进行控制，由PD控制产生控制率为在动量轮开启后，反作用轮因为吸收干扰力矩而使得角动量偏离标称值，随着贮存的动量增加，轮子达到额定转速而饱和。因此调整磁力矩器，按照的规律产生磁矩，对反作用轮产生的多余动量进行卸载，其中，Tc为控制力矩，为卫星姿态四元数的矢量部分，kp，kd为控制参数，ω为卫星角速度，M为磁矩，K为增益系数，B为卫星位置的磁感应强度，ΔHW为反作用轮产生的多余的角动量。改变其中一颗卫星的控制率，根据期望控制状态控制其反作用轮按照的控制率产生控制力矩，达到相对控制目标，其中，T为相对姿态控制力矩，为第二颗卫星相对于第一颗卫星的相对姿态四元数和期望相对姿态四元数之间的差的矢量部分，ω12为相对角速度，Kd，Kp为控制参数。
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 true CN108181916A (zh) | 2018-06-19 |
| CN108181916B 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) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109143280A (zh) * | 2018-10-10 | 2019-01-04 | 上海微小卫星工程中心 | 一种卫星集成状态闭环测试***以及相应的测试方法 |
| CN110109470A (zh) * | 2019-04-09 | 2019-08-09 | 西安电子科技大学 | 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制*** |
| CN110146082A (zh) * | 2019-05-05 | 2019-08-20 | 中国人民解放军63921部队 | 利用测速数据实时估计航天器异常姿态的方法和设备 |
| CN112461511A (zh) * | 2020-11-10 | 2021-03-09 | 中国科学院长春光学精密机械与物理研究所 | 浮空平台望远镜指向获取方法、装置、设备及存储介质 |
| CN112611373A (zh) * | 2020-12-30 | 2021-04-06 | 哈尔滨工业大学 | 一种近地空间短波红外星敏感器流场气动热效应分析方法 |
| CN113022894A (zh) * | 2021-03-08 | 2021-06-25 | 航天科工空间工程发展有限公司 | 一种用于微小卫星的相对姿态确定方法 |
| CN113204250A (zh) * | 2021-04-29 | 2021-08-03 | 西安电子科技大学 | 强动态条件下的卫星编队相对位置鲁棒高精度估计方法 |
| CN115771624A (zh) * | 2023-02-13 | 2023-03-10 | 北京航空航天大学 | 一种基于强化学习的自适应卫星姿轨控制方法 |
Citations (4)
| 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 | 西北工业大学 | 一种利用双阵列天线进行微小卫星相对姿态确定的方法 |
-
2017
- 2017-12-29 CN CN201711478457.3A patent/CN108181916B/zh active Active
Patent Citations (4)
| 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)
| Title |
|---|
| FANG YUANKUN 等: "Natural formation designof nano/micro spacecraft cluster subject to typical tasks", 《IEEE》 * |
| 吴云华: "编队飞行卫星相对姿态确定与控制方法研究", 《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》 * |
| 岳洋: "基于粒子滤波的三轴稳定卫星姿态确定算法的研究", 《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》 * |
| 方元坤 等: "一种星间相对状态估计方法", 《宇航学报》 * |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113495497A (zh) * | 2018-10-10 | 2021-10-12 | 上海微小卫星工程中心 | 卫星模拟在轨工况闭环测试*** |
| CN113495497B (zh) * | 2018-10-10 | 2024-05-14 | 上海微小卫星工程中心 | 卫星模拟在轨工况闭环测试*** |
| CN109143280A (zh) * | 2018-10-10 | 2019-01-04 | 上海微小卫星工程中心 | 一种卫星集成状态闭环测试***以及相应的测试方法 |
| CN110109470A (zh) * | 2019-04-09 | 2019-08-09 | 西安电子科技大学 | 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制*** |
| CN110109470B (zh) * | 2019-04-09 | 2021-10-29 | 西安电子科技大学 | 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制*** |
| CN110146082A (zh) * | 2019-05-05 | 2019-08-20 | 中国人民解放军63921部队 | 利用测速数据实时估计航天器异常姿态的方法和设备 |
| CN112461511B (zh) * | 2020-11-10 | 2022-03-25 | 中国科学院长春光学精密机械与物理研究所 | 浮空平台望远镜指向获取方法、装置、设备及存储介质 |
| CN112461511A (zh) * | 2020-11-10 | 2021-03-09 | 中国科学院长春光学精密机械与物理研究所 | 浮空平台望远镜指向获取方法、装置、设备及存储介质 |
| CN112611373A (zh) * | 2020-12-30 | 2021-04-06 | 哈尔滨工业大学 | 一种近地空间短波红外星敏感器流场气动热效应分析方法 |
| CN113022894A (zh) * | 2021-03-08 | 2021-06-25 | 航天科工空间工程发展有限公司 | 一种用于微小卫星的相对姿态确定方法 |
| CN113204250A (zh) * | 2021-04-29 | 2021-08-03 | 西安电子科技大学 | 强动态条件下的卫星编队相对位置鲁棒高精度估计方法 |
| CN115771624A (zh) * | 2023-02-13 | 2023-03-10 | 北京航空航天大学 | 一种基于强化学习的自适应卫星姿轨控制方法 |
| CN115771624B (zh) * | 2023-02-13 | 2023-05-26 | 北京航空航天大学 | 一种基于强化学习的自适应卫星姿轨控制方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108181916B (zh) | 2020-04-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108181916A (zh) | 小卫星相对姿态的控制方法及装置 | |
| CN102878995B (zh) | 一种静止轨道卫星自主导航方法 | |
| CN104848860B (zh) | 一种敏捷卫星成像过程姿态机动规划方法 | |
| CN104527994B (zh) | 异面交叉快变轨道固定时间稳定姿态指向跟踪控制方法 | |
| CN104898642B (zh) | 一种用于航天器姿态控制算法的集成测试仿真*** | |
| Hong | Fuzzy logic based closed-loop strapdown attitude system for unmanned aerial vehicle (UAV) | |
| CN109573105B (zh) | 末子级留轨应用子***姿态控制方法 | |
| US8131409B2 (en) | Gyroless transfer orbit sun acquisition using only wing current measurement feedback | |
| CN103090870B (zh) | 一种基于mems传感器的航天器姿态测量方法 | |
| CN101246012B (zh) | 一种基于鲁棒耗散滤波的组合导航方法 | |
| CN104374388B (zh) | 一种基于偏振光传感器的航姿测定方法 | |
| US20140231589A1 (en) | Gyroless Three-Axis Sun Acquisition Using Sun Sensor and Unscented Kalman Filter | |
| CN109823571A (zh) | 一种遥感微纳卫星的多阶段姿态控制方法 | |
| CN103112603B (zh) | 欠驱动高速自旋卫星建立正常姿态的方法 | |
| CN107402516B (zh) | 基于联合执行机构的递阶饱和模糊pd姿态控制方法 | |
| CN102809377A (zh) | 飞行器惯性/气动模型组合导航方法 | |
| CN104570742A (zh) | 基于前馈pid控制的异面交叉快变轨道快速高精度相对指向控制方法 | |
| CN104698485A (zh) | 基于bd、gps及mems的组合导航***及导航方法 | |
| CN110304279A (zh) | 一种电推进卫星的质心在轨标定补偿方法 | |
| CN106494646A (zh) | 一种磁阻尼重捕地球控制方法及*** | |
| CN110285815A (zh) | 一种可在轨全程应用的微纳卫星多源信息姿态确定方法 | |
| CN103019247A (zh) | 一种火星探测器无陀螺自主空间姿态机动控制方法 | |
| CN102425980A (zh) | 利用加速度计实现过载驾驶仪的控制方法 | |
| CN106767925A (zh) | 带双轴转位机构的惯导***三位置参数辨识对准方法 | |
| CN102114918B (zh) | 一种基于多速率敏感器组合定姿的姿控反馈回路 |
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 |