适用于卫星姿态控制的干扰力矩辨识方法

摘要

适用于卫星姿态控制的干扰力矩辨识方法，属于航天器姿态控制技术领域。为了解决传统滤波算法无法在高精度姿态控制任务中分离测量噪声和干扰力矩的问题。所述方法包括如下步骤：步骤一、根据待辨识的卫星姿态控制系统，建立带有未知干扰力矩的小量化的卫星姿态控制系统模型；步骤二、根据步骤一建立的卫星姿态控制系统模型获得未知输入观测器，采用获得的未知输入观测器估计包含噪声的干扰力矩；步骤三、采用离散傅里叶变换和离散傅里叶反变换来离线处理步骤二估计的包含噪声的干扰力矩，获得去除噪声后的干扰力矩的估计结果；步骤四、对步骤三中的估计结果采用傅里叶级数拟合得到干扰力矩的数学表达式。它用于在轨卫星进行姿态控制。

说明书

技术领域

本发明属于航天器姿态控制技术领域。

背景技术

三轴稳定卫星在轨运行期间会受到空间诸多干扰力矩的影响，例如重力梯度力矩、太 阳光压力矩、气动力矩以及卫星有效载荷转动部件产生的干扰力矩等。特别是对姿态控制 精度要求较高的卫星而言，干扰力矩会严重影响卫星的姿态控制精度，并最终导致卫星在 轨飞行任务无法顺利完成，例如高分辨率对地观测卫星。在姿态控制算法中对干扰力矩进 行补偿，能够大幅提高卫星的姿态控制精度，而对干扰力矩补偿则必须首先准确辨识出卫 星受到的干扰力矩。

目前，在轨卫星进行姿态控制时，往往利用卫星轨道参数，预先估算出干扰力矩的大 小，在轨直接在姿态控制算法中进行补偿，而不是进行实时估计和补偿。姿态控制算法对 预先估算出的干扰力矩常值进行补偿时，受轨道高度变化、转动部件运动等复杂因素影响， 无法对卫星实际受到的干扰力矩进行准确补偿，降低了卫星姿态控制的精度。干扰力矩在 轨实时估计的难点是很难将系统的噪声和干扰力矩分离，特别是当卫星处于稳态时，其受 卫星角度和角速度测量精度的影响更为明显，且执行机构输出力矩中噪声幅值可能比干扰 力矩幅值还要大。通过传统的滤波算法去噪声效果一般，有一定的时间延迟，不能够将干 扰力矩从噪声中分离出来。

发明内容

本发明的目的是为了解决传统滤波算法无法在高精度姿态控制任务中分离测量噪声和 干扰力矩的问题，本发明提供一种适用于卫星姿态控制的干扰力矩辨识方法。

本发明的适用于卫星姿态控制的干扰力矩辨识方法，所述方法包括如下步骤：

步骤一、根据待辨识的卫星姿态控制系统，建立带有未知干扰力矩的小量化的卫星姿 态控制系统模型；

步骤二、根据步骤一建立的卫星姿态控制系统模型获得未知输入观测器，采用获得的 未知输入观测器估计包含噪声的干扰力矩；

步骤三、采用离散傅里叶变换和离散傅里叶反变换来离线处理步骤二估计的包含噪声 的干扰力矩，获得去除噪声后的干扰力矩的估计结果；

步骤四、对步骤三中的估计结果采用傅里叶级数拟合得到干扰力矩的数学表达式。

所述步骤一中，根据待辨识的卫星姿态控制系统，建立带有未知干扰力矩的小量化的 卫星姿态控制系统模型：

$\left(\begin{array}{c}\stackrel{·}{x}=\mathrm{Ax}+B_{u}u+B_{d}d+w\\ y=\mathrm{Cx}+v\end{array}\right)---\left(1\right)$

其中，卫星姿态控制系统状态矩阵$A=\left(\begin{array}{cc}-A_2^{-1}{A}_0& -A_2^{-1}{A}_1\\ 0_{3\times 3}& I_{3\times 3}\end{array}\right),$输入矩阵$B_u\left(\begin{array}{c}{A}_2^{-1}{G}_u\\ 0_{3\times 3}\end{array}\right),$输入矩阵$B_d\left(\begin{array}{c}{A}_2^{-1}{G}_d\\ 0_{3\times 3}\end{array}\right),$观测矩阵C＝I_6×6，所述卫星姿态控制系统扰动w和测量噪声v为 彼此不相关的零均值白噪声；卫星姿态控制系统输入量u是3维矢量，卫星姿态控制系统 干扰量d是3维矢量；系统状态量x是6维矢量，其中 是卫星欧拉角，[ω_x ω_y ω_z]^T是卫星姿态相对参考坐标系的转动角速度在 xyz星体坐标系中的表示，为x的导数；系统输出量y是6为矢量，是系统状态量x的 测量量；

上式中，

$A_0={\omega }_0^2\left(\begin{array}{ccc}4(I_2-I_3)& & \\ & 3(I_1-I_3)& \\ & & I_2-I_1\end{array}\right)$

$A_1={\omega }_0\left(\begin{array}{ccc}& & -I_1+I_2-I_3\\ & 0& \\ I_1-I_2+I_3& & \end{array}\right)$

A₂＝I_b，G_d＝G_u＝I_3×3

其中，[0 -ω₀ 0]^T为卫星运动的牵连角速度，I_3×3为3×3阶的单位阵；

I_b是卫星转动惯量矩阵，定义为

$I_b=\left(\begin{array}{ccc}{I}_1& -I_{12}& -I_{13}\\ -I_{12}& I_2& -I_{23}\\ -I_{13}& -I_{23}& I_3\end{array}\right)$

所述步骤二中，根据步骤一建立的卫星姿态控制系统模型获得未知输入观测器，采用 获得的未知输入观测器估计包含噪声的干扰力矩：

由于式(1)中B_u＝B_d，令B_u＝B_d＝B，故式(1)离散化后为：

$\left(\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Bd}\left(k\right)+w_k\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+v_k\end{array}\right)---\left(2\right)$

根据式(2)获得未知输入观测器，采用获得的未知输入观测器具体估计为：

$\hat{d}(k-1)=M_{\mathrm{cb}}(y\left(k\right)-\mathrm{CAx}(k-1)-\mathrm{CBu}(k-1))---\left(3\right)$

式中，表示在第k步估计的第(k-1)步的干扰力矩；

矩阵M_cb满足条件：

$M_{\mathrm{cb}}\mathrm{CB}=I_{k_d\times k_d}---\left(4\right)$

其中，rank(CB)＝rank(B)＝k_d。

所述步骤二中，根据式(2)获得未知输入观测器为：

y(k+1)-CAx(k)-CBu(k)＝CBd(k)+Cw_k+v_k+1   (5)；

将式(5)两边同时左乘矩阵M_cb，得到：

M_cb(y(k+1)-CAx(k)-CBu(k))＝d(k)+M_cb(Cw_k+v_k+1)   (6)

从而获得卫星此时位置的干扰力矩的输入偏差d(k)的估计为：

$\hat{d}\left(k\right)=M_{\mathrm{cb}}(y(k+1)-\mathrm{CAx}\left(k\right)-\mathrm{CBu}\left(k\right))---\left(7\right)$

所以在第k步获得第(k-1)步的干扰力矩的估计为：

$\hat{d}(k-1)=M_{\mathrm{cb}}(y\left(k\right)-\mathrm{CAx}(k-1)-\mathrm{CBu}(k-1)).$

所述步骤一中忽略了转动惯量矩阵中的惯量积，只考虑主轴转动惯量元素：

$I_b=\left(\begin{array}{ccc}{I}_1& & \\ & I_2& \\ & & I_3\end{array}\right)---\left(8\right)$

进而得到式(1)。

所述步骤三中，采用离散傅里叶变换和离散傅里叶反变换来离线处理步骤二估计的包 含噪声的干扰力矩，获得去除噪声后的干扰力矩的估计结果包括：

采用离散傅里叶变换将步骤二获得的干扰力矩的数据序列从时域变换到频域，

选择截断频率，将高于该频率的谱值全部置零，获得去除噪声后的时间序列，所述截 断频率大于干扰力矩的变化频率且与所述变化频率同一量级；

再采用傅里叶反变换重建去除噪声后的时间序列，获得去除噪声后的干扰力矩的估计 结果。

所述步骤四中，对步骤三获得的估计结果进行傅里叶级数拟合，获得干扰力矩的数学 表达式。

本发明的有益效果在于，充分利用卫星环境力矩频率特别小的特性，提出一种基于未 知输入观测器与傅里叶变换去除噪声的干扰力矩估计算法，这种方法与传统的时延低通滤 波器相比去噪声效果更好，而且没有时间延迟。

附图说明

图1为具体实施方式一所述的适用于卫星姿态控制的干扰力矩辨识方法的原理示意 图。

图2为未知输入观测器计算的包含噪声的干扰力矩估计值随时间变化的曲线示意图。 其中，上图为x轴的干扰力矩估计值随时间变化的曲线示意图，中图为y轴的干扰力矩估 计值随时间变化的曲线示意图，下图为z轴的干扰力矩估计值随时间变化的曲线示意图。

图3为针对x轴去除噪声后的干扰力矩估计曲线。

图4为针对y轴去除噪声后的干扰力矩估计曲线。

图5为针对z轴去除噪声后的干扰力矩估计曲线。

具体实施方式

具体实施方式一：结合图1说明本实施方式，本实施方式所述的适用于卫星姿态控制 的干扰力矩辨识方法，所述方法包括如下步骤：

步骤一、根据待辨识的卫星姿态控制系统，建立带有未知干扰力矩的小量化的卫星姿 态控制系统模型；

步骤二、根据步骤一建立的卫星姿态控制系统模型获得未知输入观测器，采用获得的 未知输入观测器估计包含噪声的干扰力矩；

步骤三、采用离散傅里叶变换和离散傅里叶反变换来离线处理步骤二估计的包含噪声 的干扰力矩，获得去除噪声后的干扰力矩的估计结果；

步骤四、对步骤三中的估计结果采用傅里叶级数拟合得到干扰力矩的数学表达式。

具体实施方式二：本实施方式是对具体实施方式一所述的适用于卫星姿态控制的干扰 力矩辨识方法的进一步限定，所述步骤一的具体过程：

选择卫星欧拉角作为姿态参数，其中，是滚动角，θ是俯仰角，ψ是 偏航角，按照3-1-2的转序，即先绕z轴转ψ角，再绕x轴转角，最后绕y轴转θ角， 得到卫星的运动学方程：

式中，转动角速度ω＝[ω_x ω_y ω_z]^T；

卫星的姿态动力学方程为：

$I_b\stackrel{·}{\omega }+\omega \times \left(I_b\omega \right)=u_c+d---\left(6\right)$

式中，u_c＝[u₁ u₂ u₃]^T为星体所受控制力矩，d＝[d₁ d₂ d₃]^T表示干扰力矩。

所述姿态动力学方程可以改写成如下形式：

$\stackrel{·}{\omega }={I_b}^{-1}(u-\omega \times I_b\omega +d)---\left(7\right)$

I_b是卫星的转动惯量矩阵：

$I_b=\left(\begin{array}{ccc}{I}_1& -I_{12}& -I_{13}\\ -I_{12}& I_2& -I_{23}\\ -I_{13}& -I_{23}& I_3\end{array}\right)---\left(8\right)$

其中，对角线元素为卫星绕星体坐标轴的转动惯量，其他元素为惯量积。忽略惯量 积后，式(7)可展开为：

将上式整理成矩阵的形式为

$A_2\stackrel{··}{\alpha }+A_1\stackrel{·}{\alpha }+A_0\alpha =G_{d}d+G_{u}u---\left(10\right)$

其中，

$A_0={\omega }_0^2\left(\begin{array}{ccc}4(I_2-I_3)& & \\ & 3(I_1-I_3)& \\ & & I_2-I_1\end{array}\right)---\left(11\right)$

$A_1={\omega }_0\left(\begin{array}{ccc}& & -I_1+I_2-I_3\\ & 0& \\ I_1-I_2+I_3& & \end{array}\right)---\left(12\right)$

A₂＝I_b   (13)

G_d＝G_u＝I_3×3   (14)

令状态变量为

$x=\left(\begin{array}{c}\omega \\ \alpha \end{array}\right)---\left(15\right)$

定义卫星姿态控制系统输出变量y＝x，再考虑系统扰动和测量噪声，则建立系统的 状态方程为：

$\left(\begin{array}{c}\stackrel{·}{x}=\mathrm{Ax}+B_{u}u+B_{d}d+w\\ y=\mathrm{Cx}+v\end{array}\right)---\left(16\right)$

其中，$A=\left(\begin{array}{cc}-A_2^{-1}{A}_0& -A_2^{-1}{A}_1\\ 0_{3\times 3}& I_{3\times 3}\end{array}\right),B_u=\left(\begin{array}{c}{A}_2^{-1}{G}_u\\ 0_{3\times 3}\end{array}\right),B_d=\left(\begin{array}{c}{A}_2^{-1}{G}_d\\ 0_{3\times 3}\end{array}\right),$C＝I_6×6，I_6×6为6×6阶 的单位阵，I_3×3为3×3阶的单位阵，w和v分别表示系统扰动和测量噪声，是彼此不相关 的零均值白噪声。

具体实施方式三：本实施方式是对具体实施方式一所述的适用于卫星姿态控制的干扰 力矩辨识方法的进一步限定，所述步骤二包括如下步骤：

步骤二一、将步骤一中的模型离散化。由式(14)可知B_u＝B_d，于是可设式(16)中 B_u＝B_d＝B，取采样时间△t，得到带有未知干扰的系统离散状态空间表达式为：

$\left(\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Bd}\left(k\right)+w_k\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+v_k\end{array}\right)---\left(17\right)$

步骤二二、为了计算未知干扰输入项，必须要满足条件：

rank(CB)＝rank(B)＝k_d   (18)

而在式(17)中，由C是单位阵可知CB＝B，满足式(18)。

步骤二三、获得未知输入状态观测器。由式(17)可得：

y(k+1)-CAx(k)-CBu(k)＝CBd(k)+Cw_k+v_k+1   (19)

在步骤二二中的条件下，存在矩阵M_cb满足

$M_{\mathrm{cb}}\mathrm{CB}=I_{k_d\times k_d}---\left(20\right)$

即M_cb是矩阵CB的广义逆矩阵。将式(19)同时左乘矩阵M_cb，得到

M_cb(y(k+1)-CAx(k)-CBu(k))＝d(k)+M_cb(Cw_k+v_k+1)   (21)

从而获得卫星当前位置的干扰力矩输入偏差d(k)的估计：

$\hat{d}\left(k\right)=M_{\mathrm{cb}}(y(k+1)-\mathrm{CAx}\left(k\right)-\mathrm{CBu}\left(k\right))---\left(22\right)$

于是在第k步可以获得第(k-1)步的干扰力矩的估计公式

$\hat{d}(k-1)=M_{\mathrm{cb}}(y\left(k\right)-\mathrm{CAx}(k-1)-\mathrm{CBu}(k-1))---\left(23\right)$

卫星姿态控制系统中，测量噪声造成执行机构的控制指令力矩中包含的噪声幅值和环 境干扰力矩幅值近似或比它更大，故步骤二中估计出的结果包含环境干扰力矩和噪声两部 分。

具体实施方式四：本实施方式是对具体实施方式三所述的适用于卫星姿态控制的干扰 力矩辨识方法的进一步限定，步骤三包括如下步骤：

步骤三一、使用离散傅里叶变换将获得的干扰力矩的数据序列从时域变换到频域。由 式(23)得到的数据序列表示为{x_n}_N＝{x₀,x₁,x₂,…,x_n,…x_N-1}，其中，序列用花括号表示， xn是序列中的第n个元素(n＝0,1,2,…,N-1)。那么离散傅里叶变换得到的频域序列可以 表示为{X_k}_N＝{X₀,X₁,X₂,…X_k,…X_N-1}，其中有

$X_k=\underset{n=0}{\overset{N-1}{\Sigma }}{x}_n{w}^{-\mathrm{kn}},w=e^{j\frac{2\pi }{N}}---\left(24\right)$

式中，w满足以下性质：

w^*＝w^-1 w^l+mN＝w^l   (25)

将上述过程写成矩阵形式：

步骤三二、选择截断频率，将高于截断频率的谱值置零。由于环境干扰力矩的变化角 速度约等于轨道角速度ω_f＝10^-3rad/s，变化频率比噪声小得多，故而截断频率ω_c只需选 择和环境干扰力矩的变化频率同一量级稍大一些。

干扰力矩中太阳光压力矩和重力梯度力矩是用三角函数表示的周期函数，变化的角速 度是轨道角速度，约为ω_f＝10^-3rad/s，而当卫星姿态稳定的时候，卫星的控制器输出指 令力矩接近于高斯白噪声。所以干扰力矩的频率比白噪声的频率低的多，因而将高于轨道 频率的噪声频率基本都去掉后，剩下的就是干扰力矩的信息。

步骤三三、再将新的频域序列使用离散傅里叶反变换到时域，即可得到去除噪声后的 干扰力矩估计结果。针对频域序列{X_k}_N，有

$x_n=\frac{1}{N}\underset{k=0}{\overset{N-1}{\Sigma }}{X}_k{w}^{\mathrm{kn}}---\left(27\right)$

写成矩阵形式，可得

本发明提供一个具体实施例：

下面结合图2-5说明所述实施例，所述实施例为对以上实施方式的进一步说明，以仅 使用星敏感器为测量元件、PD控制器为控制律的卫星姿态控制系统为例，说明所设计的 干扰力矩辨识方法的合理性，具体过程为：

星敏感器的测量精度为(6,6,15)〞，干扰力矩为

$T_d=\left(\begin{array}{c}2.7\times {10}^{-4}\cos \left(0.001t\right)+4\times {10}^{-5}\\ 3.6\times {10}^{-4}\cos \left(0.001t\right)+1.8{\times 10}^{-4}\sin \left(0.001t\right)+4\times {10}^{-5}\\ -1.8\times {10}^{-4}\sin \left(0.001t\right)-4\times {10}^{-5}\end{array}\right)N·m$

卫星的转动惯量矩阵为

$I_b=\left(\begin{array}{ccc}4132.8& 0& 0\\ 0& 4322.9& 0\\ 0& 0& 7321.2\end{array}\right)\mathrm{kg}·m^2$

PD控制器参数为

Kp＝[23 24 35.8]

Kd＝[453 485 701]

本实施例采用Matlab/Simulink软件编制程序，采用ODE4算法，仿真时长设置为两 个轨道周期12566s，未知输入观测器输出值选取姿态稳定阶段的100s到仿真结束的区段， 截断频率ω_c选择0.005Hz。最后获得干扰力矩的数学表达式为

$T_{d_e}=\left(\begin{array}{c}3.977784521114195\times {10}^{-5}+2.696439242633374\times {10}^{-4}\cos \left(0.001000109588861t\right)\\ +2.069270485807362\times {10}^{-7}\sin \left(0.001000109588861t\right)\\ 4.002908671791025\times {10}^{-5}+3.600827718633706\times {10}^{-4}\cos \left(0.000999997407920t\right)\\ +1.801322753154364\times {10}^{-4}\sin \left(0.000999997407920t\right)\\ -3.952621062305346\times {10}^{-5}+1.970607900737777\times {10}^{-6}\cos \left(0.001000880148122t\right)\\ -1.799016066139477\times {10}^{-4}\sin \left(0.001000880148122t\right)\end{array}\right).$