CN108959734B - 一种基于实时递推太阳光压力矩辨识方法及系统 - Google Patents

Abstract

本发明公开了一种基于实时递推太阳光压力矩辨识方法及系统，其中，该方法包括如下步骤：使得通过姿态反馈控制维持卫星的三轴稳定姿态；建立太阳光压力矩估计模型；利用轨道系角动量方程建立太阳光压力矩和角动量的动力学方程，并将动力学方程离散化；按照固定的时间间隔采集卫星在轨的角动量,根据角动量变化计算在轨实际太阳光压力矩；然后计算三轴修正增益系数；对太阳光压估计模型系数进行修正；利用卫星地方时角级数系列矩阵更新预测误差方差；根据太阳光压系数系列得到为太阳光压系数。本发明解决了卫星在轨的太阳光压力矩计算,并提高了太阳光压力矩计算的精度，在进行太阳光压计算时不需要考虑卫星各个表面的几何形状、光学特征等因素。

Description

一种基于实时递推太阳光压力矩辨识方法及系统

技术领域

本发明属于空间飞行器动力学建模与计算技术领域，尤其涉及一种基于实时递推太阳光压力矩辨识方法及系统。

背景技术

本方法为卫星在轨的太阳光压力矩的在轨计算方法，精确获取太阳光压的大小是提高卫星控制精度和性能的有效手段。专利CN201310525249.X于2014年公开了一种高轨光学成像卫星的太阳光压计算方法,专利所述太阳光压反射模型包括用于表征卫星本体的立方体和用于表征卫星帆板的方板,通过计算立方体和帆板与太阳的关系，利用反射模型计算太阳光压。该方法在实施时需要对卫星立方体和帆板的各个表面光学反射系数进行精确试验测定。该专利所述方法依赖于卫星的光学几何形状的精度，和卫星表面光学反射率精度。然后，由于卫星结构复杂、表面并不平整等因素会导致试验结果与在轨的真实光学反射结果有较大误差，从而导致太阳光压计算误差较大。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供一种基于实时递推太阳光压力矩辨识方法及系统，解决了卫星在轨的太阳光压力矩计算,并提高了太阳光压力矩计算的精度，在进行太阳光压计算时不需要考虑卫星各个表面的几何形状、光学特征等因素。本发明利用卫星用于进行姿态控制的动量装置对太阳光压进行实时的在线估计，从而达到太阳光压计算的目的。

本发明目的通过以下技术方案予以实现：根据本发明的一个方面，提供了一种基于实时递推太阳光压力矩辨识方法，所述方法包括如下步骤：(1)在卫星上布置三轴动量装置，使得通过姿态反馈控制维持卫星的三轴稳定姿态；(2)建立太阳光压力矩估计模型，太阳光压力矩估计模型中太阳光压力矩是卫星地方时角的级数表示；(3)利用轨道系角动量方程建立太阳光压力矩和角动量的动力学方程，并将动力学方程离散化；(4)按照固定的时间间隔采集卫星在轨的角动量,根据角动量变化计算在轨实际太阳光压力矩；(5)根据步骤(2)中的太阳光压估计模型得到卫星地方时角和卫星地方时角级数系列矩阵；然后计算三轴修正增益系数；(6)利用步骤(4)中的在轨实际太阳光压力矩和太阳光压估计值对太阳光压估计模型系数进行修正；(7)利用卫星地方时角级数系列矩阵φ_k更新预测误差方差；(8)每间隔Δt时间,重复步骤(4)～步骤(7),得到时间系列下的太阳光压系数系列，根据太阳光压系数系列得到为太阳光压系数。

上述基于实时递推太阳光压力矩辨识方法中，在步骤(1)中，在卫星保持三轴姿态稳定的条件下卫星的姿态变化率

上述基于实时递推太阳光压力矩辨识方法中，在步骤(2)中，太阳光压力矩估计模型为:

φ＝[1 cos(θ_s) cos(2θ_s) cos(3θ_s) cos(4θ_s) sin(θ_s) sin(2θ_s) sin(3θ_s) sin(4θ_s)]

α^x＝[α^x0 α^x1 a^x2 α^x3 a^x4 b^x1 b^x2 b^x3 b^x4]^T

α^y＝[α^y0 α^y1 a^y2 α^y3 a^y4 b^y1 b^y2 b^y3 b^y4]^T

α^z＝[α^z0 α^z1 a^z2 α^z3 a^z4 b^z1 b^z2 b^z3 b^z4]^T

其中,表示太阳光压力矩估计值，θ_s表示卫星地方时角，α^x,α^y,α^z为需要辨识的光压力矩模型系数矩阵,φ表示卫星地方时角级数系列矩阵。

上述基于实时递推太阳光压力矩辨识方法中，在步骤(3)中，太阳光压力矩和角动量的动力学方程为：

其中，H^x,H^y,H^z为步骤(1)所述的动量装置合成的角动量在轨道系中的描述,ω_o为轨道角速度，T^x为x轴太阳光压力矩，T^y为y轴太阳光压力矩，T^z为z轴太阳光压力矩。

上述基于实时递推太阳光压力矩辨识方法中，在步骤(3)中，将动力学方程离散化后整理得：

其中，Δt为角动量采集时间间隔Δt＝t_k-t_k-1,ω_o为轨道角速度，表示t_k时刻的三轴太阳光压力矩,/>表示t_k-1时刻动量装置合成的三轴角动量，/>表示t_k时刻动量装置合成的三轴角动量,i＝x,y,z。

上述基于实时递推太阳光压力矩辨识方法中，在步骤(4)中，利用t_k-1，t_k时刻的角动量 变化计算在轨实际太阳光压力矩/>

上述基于实时递推太阳光压力矩辨识方法中，在步骤(5)中，三轴修正增益系数的公式如下：

其中，i＝x,y,z，表示t_k时刻的最小二乘法的预测误差方差，/>表示t_k时刻的测量噪声方差，Kⁱ表示三轴修正增益系数。

上述基于实时递推太阳光压力矩辨识方法中，在步骤(6)中，利用步骤(4)中的在轨实际太阳光压力矩和太阳光压估计值对太阳光压估计模型系数进行修正包括：利用t_k-1时刻系数和t_k时刻的卫星地方时角级数系列矩阵φ_k计算t_k时刻太阳光压估计值/>利用实际值与估计值误差计算t_k时刻的系数/>

其中，t_k时刻太阳光压估计值的公式如下：

t_k时刻的系数的公式如下：

其中，表示t_k时刻的a^x,a^y,a^z；/>表示t_k时刻的太阳光压力矩估计值

上述基于实时递推太阳光压力矩辨识方法中，在步骤(7)中，利用卫星地方时角级数系列矩阵φ_k更新预测误差方差通过如下公式实现：

上述基于实时递推太阳光压力矩辨识方法中，在步骤(8)中，每间隔Δt时间,重复步骤(4)～步骤(7),得到时间系列下的太阳光压系数系列，根据太阳光压系数系列得到为太阳光压系数包括：每间隔Δt时间,重复步骤(4)～步骤(7),得到时间系列t₀,t₁...t_k,t_k+1...t_n下的n+1个太阳光压系数系列当n→∞时，/>的稳态数值即为太阳光压系数。

根据本发明的另一方面，还提供了一种基于实时递推太阳光压力矩辨识系统，包括：第一模块，用于在卫星上布置三轴动量装置，使得通过姿态反馈控制维持卫星的三轴稳定姿态；第二模块，用于建立太阳光压力矩估计模型，太阳光压力矩估计模型中太阳光压力矩是卫星地方时角的级数表示；第三模块，用于利用轨道系角动量方程建立太阳光压力矩和角动量的动力学方程，并将动力学方程离散化；第四模块，用于按照固定的时间间隔采集卫星在轨的角动量,根据角动量变化计算在轨实际太阳光压力矩；第五模块，用于根据第二模块中的太阳光压估计模型得到卫星地方时角和卫星地方时角级数系列矩阵；然后计算三轴修正增益系数；第六模块，用于利用第上模块中的在轨实际太阳光压力矩和太阳光压估计值对太阳光压估计模型系数进行修正；第七模块，用于利用卫星地方时角级数系列矩阵φ_k更新预测误差方差；第八模块，用于得到时间系列下的太阳光压系数系列，根据太阳光压系数系列得到为太阳光压系数。

根据本发明的又一个方面，提供了一个或多个机器可读介质，其上存储有指令，当由一个或多个处理器执行时，使得设备执行本发明的一个方面的一个或多个的方法。

本发明与现有技术相比具有如下有益效果：与现有的技术相比，本发明不需要具体知道卫星的外形尺寸、表面光学特性等参数。本方法只需要利用卫星用于进行姿态控制的三轴动量交换装置通过稳定的姿态控制，让动量装置吸收太阳光压力矩，然后通过动量装置的变化规律实时辨识计算太阳光压。本方法是一种递推的算法，当递推稳定后即可获得太阳光压的大小，并且随着太阳季节的变化、太阳照射方向变化等因素本算法会依据在轨进行连续实时调整，这是基于几何特征的方法不能做到的，因而能够获得连续的较高的太阳光压力矩精度。

附图说明

通过阅读下文优选实施方式的详细描述，各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。而且在整个附图中，用相同的参考符号表示相同的部件。在附图中：

图1为本发明基于实时递推太阳光压力矩辨识方法的流程框图；

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。

本实施例提供了一种基于实时递推太阳光压力矩辨识方法，包括如下步骤：

(1)在卫星上布置三轴动量装置，动量装置的作用是通过姿态反馈控制维持卫星的三轴稳定姿态。在卫星保持三轴姿态稳定的条件下卫星的姿态变化率因此卫星自生的动量变化可以忽略。在闭环控制作用下太阳光压力矩干扰被动量装置完全吸收，动量装置的动量变化反映了太阳光压力矩的大小。动量装置可以是动量轮群、控制力矩陀螺群或者动量轮和控制力矩陀螺混合群构成。

(2)建立太阳光压力矩估计模型，太阳光压力矩估计模型中太阳光压力矩是卫星地方时角的级数表示。太阳光压估计模型为:

φ＝[1 cos(θ_s) cos(2θ_s) cos(3θ_s) cos(4θ_s) sin(θ_s) sin(2θ_s) sin(3θ_s) sin(4θ_s)] (2)

α^x＝[α^x0 α^x1 a^x2 α^x3 a^x4 b^x1 b^x2 b^x3 b^x4]^T

α^y＝[α^y0 α^y1 a^y2 α^y3 a^y4 b^y1 b^y2 b^y3 b^y4]^T

α^z＝[α^z0 α^z1 a^z2 α^z3 a^z4 b^z1 b^z2 b^z3 b^z4]^T

其中,表示太阳光压力矩估计值，θ_s表示卫星地方时角，α^x,α^y,α^z为需要辨识的光压力矩模型系数矩阵,φ表示卫星地方时角级数系列矩阵。

卫星在轨飞行过程中，太阳光与卫星各个表面发生作用从而产生光压力矩，随着卫星的运动太阳光方向与卫星各个面夹角变化，从而导致力矩变化，即卫星各个表面法线矢量与太阳光矢量之间的夹角具有明显的周期性，光压力矩也相应的表现为周期性。因此，将太阳光压力矩的大小采用为多阶傅立叶级数进行表示。本方法中取4阶，实际使用可不限于4阶。

(3)利用轨道系角动量方程建立太阳光压力矩和角动量的动力学关系，并将动力学方程离散化：

将上述方程离散化后整理得：

其中,H^x,H^y,H^z为步骤(1)所述的动量装置合成的角动量在轨道系中的描述,Δt为角动量采集时间间隔Δt＝t_k-t_k-1,ω_o为轨道角速度。表示t_k时刻的三轴太阳光压力矩,表示t_k-1时刻动量装置合成的三轴角动量，/>表示t_k时刻动量装置合成的三轴角动量,i＝x,y,z。

(4)按照固定的时间间隔采集卫星在轨的角动量,利用t_k-1，t_k时刻的角动量 变化计算在轨实际太阳光压力矩/>具体计算利用步骤(3)之公式(4)。

(5)计算卫星地方时角和卫星地方时角级数系列矩阵,然后计算三轴修正增益系数。具体为利用步骤(2)之公式(2)计算卫星地方时角和卫星地方时角级数系列矩阵φ_k，然后按照如下公式计算三轴修正增益系数。具体如下：

这里表示t_k时刻的最小二乘法的预测误差方差，/>表示t_k时刻的测量噪声方差，Kⁱ表示三轴修正增益系数.

(6)利用在轨的实际太阳光压力矩和太阳光压估计值对太阳光压模型系数进行修正。具体为利用t_k-1时刻系数和t_k时刻的卫星地方时角级数系列矩阵φ_k计算t_k时刻太阳光压估计值/>利用实际值与估计值误差计算t_k时刻的系数/>

这里，表示t_k时刻的a^x,a^y,a^z；/>表示t_k时刻的太阳光压力矩估计值在步骤(6)的计算中，需要对推力器的喷气力矩进行排除。例如，当|T_k|>0.005时认为是推力器产生的干扰，则当前数据不引入计算/>维持/>的数值，保证递推过程稳定连续。

(7)利用卫星地方时角级数系列矩阵φ_k更新预测误差方差。具体如下：

(8)每间隔Δt时间,重复步骤(4)～步骤(7),这样计算得到时间系列t₀,t₁...t_k,t_k+1...t_n下的n+1个太阳光压系数系列当n足够大时，/>的稳态数值即为太阳光压系数，从而达到辨识的目的。这里第一步/>的数值一般取为0，n的数值一般需要满足n·Δt大于一个轨道周期。例如对于同步轨道卫星一个轨道周期为86400秒，那么应满足n·Δt>86400。

本实施例还提供了一种基于实时递推太阳光压力矩辨识系统，该系统包括：第一模块，用于在卫星上布置三轴动量装置，使得通过姿态反馈控制维持卫星的三轴稳定姿态；第二模块，用于建立太阳光压力矩估计模型，太阳光压力矩估计模型中太阳光压力矩是卫星地方时角的级数表示；第三模块，用于利用轨道系角动量方程建立太阳光压力矩和角动量的动力学方程，并将动力学方程离散化；第四模块，用于按照固定的时间间隔采集卫星在轨的角动量,根据角动量变化计算在轨实际太阳光压力矩；第五模块，用于根据第二模块中的太阳光压估计模型得到卫星地方时角和卫星地方时角级数系列矩阵；然后计算三轴修正增益系数；第六模块，用于利用第上模块中的在轨实际太阳光压力矩和太阳光压估计值对太阳光压估计模型系数进行修正；第七模块，用于利用卫星地方时角级数系列矩阵φ_k更新预测误差方差；第八模块，用于得到时间系列下的太阳光压系数系列，根据太阳光压系数系列得到为太阳光压系数。

本实施例还提供了一个或多个机器可读介质，其上存储有指令，当由一个或多个处理器执行时，使得设备执行上述的一个或多个的方法。

本实施例不需要具体知道卫星的外形尺寸、表面光学特性等参数。本方法只需要利用卫星用于进行姿态控制的三轴动量交换装置通过稳定的姿态控制，让动量装置吸收太阳光压力矩，然后通过动量装置的变化规律实时辨识计算太阳光压。本方法是一种递推的算法，当递推稳定后即可获得太阳光压的大小，并且随着太阳季节的变化、太阳照射方向变化等因素本算法会依据在轨进行连续实时调整，这是基于几何特征的方法不能做到的，因而能够获得连续的较高的太阳光压力矩精度。

以上所述的实施例只是本发明较优选的具体实施方式，本领域的技术人员在本发明技术方案范围内进行的通常变化和替换都应包含在本发明的保护范围内。

Claims (10)

1.一种基于实时递推太阳光压力矩辨识方法，其特征在于，所述方法包括如下步骤：

(1)在卫星上布置三轴动量装置，使得通过姿态反馈控制维持卫星的三轴稳定姿态；

(2)建立太阳光压力矩估计模型，太阳光压力矩估计模型中太阳光压力矩是卫星地方时角的级数表示；

(3)利用轨道系角动量方程建立太阳光压力矩和角动量的动力学方程，并将动力学方程离散化；

(4)按照固定的时间间隔采集卫星在轨的角动量,根据角动量变化计算在轨实际太阳光压力矩；

(5)根据步骤(2)中的太阳光压力矩估计模型得到卫星地方时角和卫星地方时角级数系列矩阵；然后计算三轴修正增益系数；

(6)利用步骤(4)中的在轨实际太阳光压力矩和太阳光压估计值对太阳光压力矩估计模型系数进行修正；

(7)利用卫星地方时角级数系列矩阵更新预测误差方差；

(8)每间隔Δt时间,重复步骤(4)～步骤(7),得到时间系列下的太阳光压系数系列，根据太阳光压系数系列得到为太阳光压系数。

2.根据权利要求1所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(1)中，在卫星保持三轴姿态稳定的条件下卫星的姿态变化率

3.根据权利要求1所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(2)中，太阳光压力矩估计模型为:

φ＝[1 cos(θ_s) cos(2θ_s) cos(3θ_s) cos(4θ_s) sin(θ_s) sin(2θ_s) sin(3θ_s) sin(4θ_s)]

α^x＝[α^x0 α^x1 a^x2 α^x3 a^x4 b^x1 b^x2 b^x3 b^x4]^T

α^y＝[α^y0 α^y1 a^y2 α^y3 a^y4 b^y1 b^y2 b^y3 b^y4]^T

α^z＝[α^z0 α^z1 a^z2 α^z3 a^z4 b^z1 b^z2 b^z3 b^z4]^T

其中,表示太阳光压力矩估计值，θ_s表示卫星地方时角，α^x,α^y,α^z为需要辨识的太阳光压力矩估计模型系数矩阵,φ表示卫星地方时角级数系列矩阵。

4.根据权利要求3所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(3)中，太阳光压力矩和角动量的动力学方程为：

其中，H^x,H^y,H^z为步骤(1)所述的动量装置合成的角动量在轨道系中的描述,ω_o为轨道角速度，T^x为x轴太阳光压力矩，T^y为y轴太阳光压力矩，T^z为z轴太阳光压力矩。

5.根据权利要求4所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(3)中，将动力学方程离散化后整理得：

其中，Δt为角动量采集时间间隔Δt＝t_k-t_k-1,ω_o为轨道角速度，表示t_k时刻的三轴太阳光压力矩,/>表示t_k-1时刻动量装置合成的三轴角动量，/>表示t_k时刻动量装置合成的三轴角动量,i＝x,y,z。

6.根据权利要求5所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(4)中，利用t_k-1，t_k时刻的角动量 变化计算在轨实际太阳光压力矩/>

7.根据权利要求6所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(5)中，三轴修正增益系数的公式如下：

其中，i＝x,y,z，表示t_k时刻的最小二乘法的预测误差方差，/>表示t_k时刻的测量噪声方差，Kⁱ表示三轴修正增益系数，φ_k为t_k时刻的卫星地方时角级数系列矩阵。

8.根据权利要求7所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(6)中，利用步骤(4)中的在轨实际太阳光压力矩和太阳光压力矩估计值对太阳光压力矩估计模型系数进行修正包括：利用t_k-1时刻系数和t_k时刻的卫星地方时角级数系列矩阵φ_k计算t_k时刻太阳光压力矩估计值/>利用实际值与估计值误差计算t_k时刻的系数/>

其中，t_k时刻太阳光压力矩估计值的公式如下：

t_k时刻的系数的公式如下：

其中，i＝x,y,z表示t_k时刻的/>表示t_k时刻的太阳光压力矩估计值

9.根据权利要求8所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(7)中，利用卫星地方时角级数系列矩阵φ_k更新预测误差方差通过如下公式实现：

10.根据权利要求9所述的基于实时递推太阳光压力矩辨识方法，其特征在于：在步骤(8)中，每间隔Δt时间,重复步骤(4)～步骤(7),得到时间系列下的太阳光压系数系列，根据太阳光压系数系列得到为太阳光压系数包括：每间隔Δt时间,重复步骤(4)～步骤(7),得到时间系列t₀,t₁...t_k,t_k+1...t_n下的n+1个太阳光压系数系列当n→∞时，/>的稳态数值即为太阳光压系数。