CN104898686A - 一种基于反作用轮摩擦特性的抗干扰姿态控制方法 - Google Patents

Abstract

一种基于反作用轮摩擦特性的抗干扰姿态控制方法，针对反作用轮过零摩擦以及姿控精度较低的问题，首先构建含有反作用轮摩擦特性的航天器耦合动力学方程，其次针对由于反作用轮带来的摩擦干扰设计摩擦干扰估计器，对反作用轮摩擦进行实时估计；再次，针对航天器系统中的摩擦干扰估计误差及环境干扰，设计标称控制器进行抑制；最后通过联合标称控制器与摩擦干扰估计器，设计复合分层抗干扰控制器，完成多源干扰影响下的航天器抗干扰姿态控制；本方法在传统的航天器动力学方程中加入执行机构动力学，能够更好的分析由于执行机构特性带来的控制精度下降问题，可用于航空航天领域高精度对地观测卫星，太空望远镜等航天器的高精度姿态控制。

Description

一种基于反作用轮摩擦特性的抗干扰姿态控制方法

技术领域

本发明涉及一种基于反作用轮摩擦特性的抗干扰姿态控制方法，特别是针对包含反作用轮摩擦、环境力矩等多源干扰下的航天器抗干扰姿态控制方法，该方法可用于航空航天领域高精度观测卫星，太空望远镜等航天器的高精度姿态控制。

背景技术

航天器姿态控制系统的可靠性和长期工作稳定性一直是航天器研制中的关键课题。反作用轮是航天器姿态控制系统中最重要的惯性执行部件之一，近几年发射的长寿命、高精度、多功能三轴稳定卫星，几乎毫不例外地利用反作用轮作为主要执行部件。在航天器姿控系统中，反作用轮按照姿控系统指令，提供合适的控制力矩，校正航天器的姿态偏差，完成某种预定的姿态调整。反作用轮作为航天器三轴稳定控制系统主要的执行部件，其性能对姿态控制系统具有决定性的影响，提高反作用轮控制系统的性能对卫星姿态控制系统具有重要意义。反作用轮系统内部存在的多种干扰和噪声，如摩擦力矩干扰和速率测量噪声等，会使反作用轮执行指令产生误差，给航天器姿态控制系统带来干扰。此外航天器在轨飞行过程中，除了受到反作用轮带来的摩擦力矩干扰等，还将受到各种外部太空环境干扰力矩的作用，如气动力矩、太阳辐射光压、重力梯度力矩、空间碎片撞击力矩以及磁力矩等。对于航天器系统而言，来自上述外部环境干扰以及反作用轮的干扰，都将严重影响到姿态控制系统的精度和稳定度，在航天器设计的过程中，必须考虑上述干扰的影响。

航天器姿态控制器目前设计过程一般都是基于模型进行设计，准确的模型参数能够提高航天器的姿态控制精度，以往的研究过程中姿态控制系统的设计与执行机构的特性联系并不紧密，忽略了执行机构动态特性，此外对于航天器所遭受的太空环境干扰，也忽略部分干扰特性，不够精细，使得航天器动力学模型在精确性上有所欠缺，导致所设计的控制方法不够精确，影响航天器系统的姿态精度和稳定度，在面向实际工程对象时，会存在问题。

发明内容

本发明的技术解决问题是：针对使用反作用轮的航天器系统，克服现有技术的不足，并充分利用反作用轮的动态特性，建立包含反作用轮摩擦特性的抗干扰姿态控制方法，并同时针对航天系统存在的外部环境干扰力矩以及反作用摩擦等多源干扰，提供一种基于反作用动态特性的抗干扰姿态控制方法，解决多源干扰下航天器的姿态控制问题，提高航天器姿态系统的控制精度和稳定度。

本发明的技术解决方案为：一种基于反作用轮摩擦特性的抗干扰姿态控制方法，其特征在于包括以下步骤：针对反作用轮摩擦以及姿控精度较低的问题，首先构建含有反作用轮摩擦特性的航天器耦合动力学方程，其次针对由于反作用轮带来的摩擦干扰设计摩擦干扰估计器，对反作用轮摩擦进行实时估计；再次，针对航天器系统中存在的环境干扰力矩以及摩擦干扰估计误差，设计标称控制器进行抑制；最后通过联合标称控制器与摩擦干扰估计器，设计复合分层抗干扰控制器，完成多源干扰影响下的航天器抗干扰姿态控制；

第一步：构建含有反作用轮摩擦特性的航天器耦合动力学模型

当航天器本体坐标系和轨道坐标系之间的欧拉角很小时,姿态运动学可以小角度线性化,得到如下的姿态动力学模型：

${\mathrm{\Σ}}_1:\left\{\begin{array}{c}{I}_1\stackrel{\·\·}{\mathrm{\φ}}\left(t\right)-n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\ψ}}\left(t\right)+4n^2(I_2-I_3)\mathrm{\φ}\left(t\right)=J_1{\stackrel{\·}{\mathrm{\Ω}}}_1\left(t\right)+T_{d1}\left(t\right)\\ I_2\stackrel{\·\·}{\mathrm{\θ}}\left(t\right)+3n^2(I_1-I_3)\mathrm{\θ}\left(t\right)=J_2{\stackrel{\·}{\mathrm{\Ω}}}_2\left(t\right)+T_{d2}\left(t\right)\\ I_3\stackrel{\·\·}{\mathrm{\ψ}}\left(t\right)+n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\φ}}\left(t\right)+n^2(I_2-I_1)\mathrm{\ψ}\left(t\right)=J_3{\stackrel{\·}{\mathrm{\Ω}}}_3\left(t\right)+T_{d3}\left(t\right)\end{array}\right.$

其中，I₁,I₂,I₃分别为航天器的三轴转动惯量；φ(t),θ(t),ψ(t)分别为三轴姿态角，分别为三轴姿态角速度，分别为三轴姿态角加速度；J₁,J₂,J₃分别为三轴上安装的反作用轮的转动惯量；分别表示航天器三轴上的反作用轮中电机的转动角加速度，T_d1(t),T_d2(t),T_d3(t)分别为航天器三轴受到的环境干扰；n为航天器轨道角速度；

用于航天器姿态控制的反作用轮系统一般是由微型直流无刷电机、反作用轮本体、以及控制器线路构成，是一个相对独立的机电执行部件,反作用轮的转动方程为：

$\left\{\begin{array}{c}{J}_i{\stackrel{\·}{\mathrm{\Ω}}}_i\left(t\right)=k_{T_i}*i_i\left(t\right)+M_{\mathrm{fi}}\left(t\right)\\ E_i\left(t\right)=k_{\mathrm{Ei}}*{\mathrm{\Ω}}_i\left(t\right)+R_i*i_i\left(t\right)+L_i*{\stackrel{\·}{i}}_i\left(t\right)\end{array}\right.$

其中，J_i为反作用轮的转动惯量,Ω_i(t)为反作用轮中电机的转动角速度，为反作用轮中电机的转动角加速度，k_Ti为反作用轮中电机的电流力矩系数，i_i(t)为反作用轮中电机的电流，k_Ti*i_i(t)为反作用轮的驱动力矩，M_fi(t)为反作用轮摩擦力矩；对于反作用轮使用的无刷直流电机，电流力矩系数k_Ti在数值上与反作用轮中的反电动势系数相等,k_Ei为反作用轮使用电机的反电动势系数，E_i(t)为加在反作用轮两端的控制电压，R_i为反作用轮使用电机的电阻值，L_i为反作用轮使用电机的电感，在正常使用过程中，在数值上大小接近0，在实际处理过程中，这一项可以忽略，因此反作用轮转动方程进一步可以描述为；

$\left\{\begin{array}{c}{J}_i{\stackrel{\·}{\mathrm{\Ω}}}_i\left(t\right)=k_{\mathrm{Ti}}*i_i\left(t\right)+M_{\mathrm{fi}}\left(t\right)\\ E_i\left(t\right)=k_{\mathrm{Ei}}*{\mathrm{\Ω}}_i\left(t\right)+R_i*i_i\left(t\right)\end{array}\right.$

对于分装在航天器三轴上的反作用轮，航天器的控制力矩指令可以表示为：

${\mathrm{\Σ}}_2:\left\{\begin{array}{c}{J}_1{\stackrel{\·}{\mathrm{\Ω}}}_1\left(t\right)=k_{T1}{i}_1\left(t\right)+M_{f1}\left(t\right)\\ J_2{\stackrel{\·}{\mathrm{\Ω}}}_2\left(t\right)=K_{T2}{i}_2\left(t\right)+M_{f2}\left(t\right)\\ J_3{\stackrel{\·}{\mathrm{\Ω}}}_3\left(t\right)=K_{T3}{i}_3\left(t\right)+M_{f3}\left(t\right)\end{array}\right.$

其中，k_T1,k_T2,k_T3为安装在航天器三轴上的反作用轮中电机的电流力矩系数，i₁(t),i₂(t),i₃(t)为安装在航天器三轴上的反作用轮中电机的电流，M_f1(t),M_f2(t),M_f3(t)为安装在航天器三轴上的反作用轮的摩擦力矩；进一步，三轴上安装的反作用轮中所用电机的电流方程为：

${\mathrm{\Σ}}_3:\left\{\begin{array}{c}{i}_1\left(t\right)=\frac{1}{R_1}[E_1\left(t\right)-k_{E1}({\mathrm{\Ω}}_1\left(t\right)-\stackrel{\·}{\mathrm{\φ}}\left(t\right))]\\ i_2\left(t\right)=\frac{1}{R_2}[E_2\left(t\right)-k_{E2}({\mathrm{\Ω}}_2\left(t\right)-\stackrel{\·}{\mathrm{\θ}}\left(t\right))]\\ i_3\left(t\right)=\frac{1}{R_3}[E_3\left(t\right)-k_{E3}({\mathrm{\Ω}}_3\left(t\right)-\stackrel{\·}{\mathrm{\ψ}}\left(t\right))]\end{array}\right.$

其中，R₁,R₂,R₃为安装在航天器三轴上的反作用轮电机电阻，E₁(t),E₂(t),E₃(t)为安装在航天器三轴上的三个反作用轮的控制电压；k_E1,k_E2,k_E3是安装在航天器三轴上的三个反作用轮中电机的反电动势系数；

进一步，建立包含反作用轮摩擦特性的航天器姿态系统动力学模型为：

${\mathrm{\Σ}}_4:\left\{\begin{array}{c}{I}_1\stackrel{\·\·}{\mathrm{\φ}}\left(t\right)-n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\ψ}}\left(t\right)-\frac{k_{T1}*k_{E1}}{R_1}\stackrel{\·}{\mathrm{\φ}}\left(t\right)+4n^2(I_2-I_3)\mathrm{\φ}\left(t\right)\\ =\frac{k_{T1}}{R_1}{E}_1\left(t\right)-\frac{k_{T1}*k_{E1}}{R_1}*{\mathrm{\Ω}}_1\left(t\right)+M_{f1}\left(t\right)+T_{d1}\left(t\right)\\ I_2\stackrel{\·\·}{\mathrm{\θ}}\left(t\right)+3n^2(I_1-I_3)\mathrm{\θ}\left(t\right)-\frac{k_{T2}*k_{E2}}{R_2}\stackrel{\·}{\mathrm{\θ}}\left(t\right)\\ =\frac{k_{T2}}{R_2}{E}_2\left(t\right)-\frac{k_{T2}*k_{E2}}{R_2}*{\mathrm{\Ω}}_2\left(t\right)+M_{f2}\left(t\right)+T_{d2}\left(t\right)\\ I_3\stackrel{\·\·}{\mathrm{\ψ}}\left(t\right)+n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\φ}}\left(t\right)+n^2(I_2-I_1)\mathrm{\ψ}\left(t\right)-\frac{k_{T3}*k_{E3}}{R_3}\stackrel{\·}{\mathrm{\ψ}}\left(t\right)\\ =\frac{k_{T3}}{R_3}{E}_3\left(t\right)-\frac{k_{T3}*k_{E3}}{R_3}*{\mathrm{\Ω}}_3\left(t\right)+M_{f3}\left(t\right)+T_{d3}\left(t\right)\end{array}\right.$

将式∑₄写成矩阵形式如下：

$P\stackrel{\·\·}{p}\left(t\right)+C\stackrel{\·}{p}\left(t\right)+\mathrm{Qp}\left(t\right)=B_{E}E\left(t\right)+B_{\mathrm{\Ω}}\mathrm{\Ω}\left(t\right)+B_m{M}_f\left(t\right)+B_w{T}_d\left(t\right)$

其中，状态p(t)＝[φ(t),θ(t),ψ(t)]^T，T_d(t)＝[T_d1(t),T_d2(t),T_d3(t)]^T为航天器三轴受到的环境干扰力矩，E(t)＝[E₁(t),E₂(t),E₃(t)]^T为安装在航天器三轴上的三个反作用轮的控制电压，Ω(t)＝[Ω₁(t),Ω₂(t),Ω₃(t)]^T为安装在航天器三轴上的反作用轮中电机的角速度，M_f(t)＝[M_f1(t),M_f2(t),M_f3(t)]^T为安装在航天器三轴上的反作用轮的摩擦力矩干扰；式中P＝diag{I₁,I₂,I₃}，B_m＝I,B_w＝I，Q＝diag{4n²(I₂-I₃),3n²(I₁-I₃),n²(I₂-I₁)}, $B_{\mathrm{\Ω}}=\mathrm{diag}\{-\frac{k_{T1}*k_{E1}}{R_1},-\frac{k_{T2}*k_{E2}}{R_2},-\frac{k_{T3}*k_{E3}}{R_3}\},$ $C=\left[\begin{array}{ccc}-\frac{k_{T1}*k_{E1}}{R_1}& 0& -n(I_1-I_2+I_3)\\ 0& -\frac{k_{T2}*k_{E2}}{R_2}& 0\\ n(I_1-I_2+I_3)& 0& -\frac{k_{T3}*k_{E3}}{R_3}\end{array}\right];$

进一步可以将包含反作用轮摩擦特性的航天器姿态系统动力学模型写成状态空间表达式形式如下：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+B_{0}E\left(t\right)+B_1{M}_f\left(t\right)+B_2{T}_d\left(t\right)$

其中，x(t)、E(t)分别为系统的状态变量和控制输入，E(t)＝[E₁(t),E₂(t),E₃(t)]^T, $x\left(t\right)={\left[\begin{array}{cccc}{e}_p\left(t\right)& {\∫}_{t-1}^t\mathrm{\Ω}\left(\mathrm{\τ}\right)d\left(t\right)& {\stackrel{\·}{e}}_p\left(t\right)& \mathrm{\Ω}\left(t\right)\end{array}\right]}^T,$ 其中e_p(t)＝p(t)-p_d(t)，p_d(t)为航天器系统的期望姿态信息，p(t)为航天器系统当前最新的姿态信息，e_p(t)为航天器当前姿态与期望姿态的偏差，A为系统阵，B₀为控制输入矩阵，B₁为摩擦力矩干扰输入矩阵，B₂为环境干扰力矩的输入矩阵，其对应的表达式如下：

$A=\left[\begin{array}{cccc}0& 0& I& 0\\ 0& 0& 0& I\\ -P^{-1}Q& 0& -P^{-1}C& P^{-1}{B}_{\mathrm{\Ω}}\\ 0& 0& \frac{K_T*K_E}{J*R}& -\frac{K_T*K_E}{J*R}\end{array}\right],B_0=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_E\\ \frac{K_T}{J*R}\end{array}\right],B_1=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_m\\ \frac{1}{J}\end{array}\right],B_2=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_w\\ 0\end{array}\right];$

其中，K_T＝diag{k_T1,k_T2,k_T3}，J＝diag{J₁,J₂,J₃}，K_E＝diag{k_E1,k_E2,k_E3}，R＝diag{R₁,R₂,R₃}；

第二步：设计摩擦干扰估计器

对于上述反作用轮摩擦干扰，设计摩擦干扰估计器，对摩擦干扰的估计表征为：

$\left\{\begin{array}{c}{\hat{M}}_f\left(t\right)=r\left(t\right)+\mathrm{Nx}\left(t\right)\\ \stackrel{\·}{r}\left(t\right)=-N{B}_1(r\left(t\right)+\mathrm{Nx}\left(t\right))-N(\mathrm{Ax}\left(t\right)+B_0\left(t\right))\end{array}\right.$

其中，N为摩擦干扰估计器增益，为反作用轮摩擦干扰的估计值，r(t)为摩擦干扰估计器中的辅助变量；

第三步：设计标称控制器

在通过摩擦干扰估计器完成对反作用轮摩擦进行估计后，对于航天器系统中还存在的环境干扰力矩以及摩擦干扰估计误差干扰，进一步设计具有干扰抑制能力的标称控制器，其对应控制律为：

$E_K\left(t\right)=K{\left[\begin{array}{cccc}{e}_p\left(t\right)& {\∫}_{t-1}^t\mathrm{\Ω}\left(\mathrm{\τ}\right)d\left(t\right)& {\stackrel{\·}{e}}_p\left(t\right)& \mathrm{\Ω}\left(t\right)\end{array}\right]}^T$

其中，E_K(t)为当前时刻标称控制器的控制量，K为标称控制器的控制增益；

第四步：设计航天器复合分层抗干扰姿态控制器

在分别完成反作用轮摩擦观测器与标称控制器基础上，进行复合，得到多源干扰下航天器复合分层抗干扰姿态控制器，其对应的控制律为：

$E\left(t\right)=E_K\left(t\right)-\frac{R}{K_T}{\hat{M}}_f\left(t\right)$

其中，E_K(t)为复合分层抗干扰姿态控制器的输出，用于控制安装在航天器三轴上的反作用轮完成姿态控制；

本发明与现有技术相比的优点在于：本发明在航天器动力学建模过程中，加入了执行机构反作用轮的动态特性，相对于目前存在的姿态控制方法，充分利用已有的信息，更贴合工程实际，能更好的解决工程中出现的问题；此外，针对航天器系统还存在的外部环境干扰力矩，通过引入摩擦干扰估计器与标称控制器构成复合控制器，对摩擦干扰的估计和补偿与对航天器系统受到的外部环境干扰的抑制相结合，能够有效的提高航天器的姿态控制精度和稳定度。

附图说明

图1为本发明一种基于反作用轮动态特性的抗干扰姿态控制方法的设计流程图。

具体实施方式

如图1所示，本发明具体实现步骤如下(卫星作为一种典型的航天器，本实施例中以一类微纳卫星为例来说明方法的具体实现)：

1、构建含有反作用轮摩擦特性的卫星耦合动力学模型

当卫星本体坐标系和轨道坐标系之间的欧拉角很小时,姿态运动学可以小角度线性化,得到如下的姿态动力学模型：

${\mathrm{\Σ}}_1:\left\{\begin{array}{c}{I}_1\stackrel{\·\·}{\mathrm{\φ}}\left(t\right)-n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\ψ}}\left(t\right)+4n^2(I_2-I_3)\mathrm{\φ}\left(t\right)=J_1{\stackrel{\·}{\mathrm{\Ω}}}_1\left(t\right)+T_{d1}\left(t\right)\\ I_2\stackrel{\·\·}{\mathrm{\θ}}\left(t\right)+3n^2(I_1-I_3)\mathrm{\θ}\left(t\right)=J_2{\stackrel{\·}{\mathrm{\Ω}}}_2\left(t\right)+T_{d2}\left(t\right)\\ I_3\stackrel{\·\·}{\mathrm{\ψ}}\left(t\right)+n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\φ}}\left(t\right)+n^2(I_2-I_1)\mathrm{\ψ}\left(t\right)=J_3{\stackrel{\·}{\mathrm{\Ω}}}_3\left(t\right)+T_{d3}\left(t\right)\end{array}\right.$

其中，I₁,I₂,I₃分别为本实施例中微纳卫星的三轴转动惯量，其对应数值为I₁＝5.5，I₂＝6.2，I₃＝2.2；φ(t),θ(t),ψ(t)分别为三轴姿态角，分别为三轴姿态角速度，分别为三轴姿态角加速度；J₁,J₂,J₃分别为三轴上安装的反作用轮的转动惯量，在此实施例中，采用同一类型反作用轮，其转动惯量对应数值为J₁＝0.00273，J₂＝0.00273，J₃＝0.00273；分别表示卫星三轴上的反作用轮的转速，T_d1(t),T_d2(t),T_d3(t)分别为航天器三轴受到的环境干扰力矩，其对应数值为T_d1＝2.4×10^-4(3cos10t+1)，T_d2＝2.4×10^-4(3cos10t++1.6sin10t)，T_d3＝2.3×10^-4(3sin10t+1)，n为航天器轨道角速度，此实施例中n＝0.0011rad/s；

用于卫星姿态控制的反作用轮系统一般是由微型直流无刷电机、反作用轮本体、以及控制器线路构成，是一个相对独立的机电执行部件,反作用轮的转动方程为：

$\left\{\begin{array}{c}{J}_i{\stackrel{\·}{\mathrm{\Ω}}}_i\left(t\right)=k_{T_i}*i_i\left(t\right)+M_{\mathrm{fi}}\left(t\right)\\ E_i\left(t\right)=k_{\mathrm{Ei}}*{\mathrm{\Ω}}_i\left(t\right)+R_i*i_i\left(t\right)+L_i*{\stackrel{\·}{i}}_i\left(t\right)\end{array}\right.$

其中，J_i为反作用轮的转动惯量，Ω_i(t)为反作用轮中电机的转动角速度，为反作用轮中电机的转动角加速度，k_Ti为反作用轮中电机的电流力矩系数，i_i(t)为反作用轮中电机的电流，k_Ti*i_i(t)为反作用轮的驱动力矩，M_fi(t)为反作用轮摩擦力矩；对于反作用轮使用的无刷直流电机，电流力矩系数k_Ti在数值上与反作用轮中的反电动势系数相等,k_Ei为反作用轮使用电机的反电动势系数，E_i(t)为加在反作用轮两端的控制电压，R_i为反作用轮使用电机的电阻值，L_i为反作用轮使用电机的电感，在此实施例中，三轴安装的反作用轮均为同一类型电机驱动，反作用轮的转动惯量数值大小为0.00273，对应的反电动势系数均为0.02，电流力矩系数均为0.02，电机中电阻值均为1Ω，电感大小均为0.0002H，在正常使用过程中，在数值上大小接近0，在实际处理过程中，这一项可以忽略，因此反作用轮转动方程进一步可以描述为；

$\left\{\begin{array}{c}{J}_i{\stackrel{\·}{\mathrm{\Ω}}}_i\left(t\right)=k_{\mathrm{Ti}}*i_i\left(t\right)+M_{\mathrm{fi}}\left(t\right)\\ E_i\left(t\right)=k_{\mathrm{Ei}}*{\mathrm{\Ω}}_i\left(t\right)+R_i*i_i\left(t\right)\end{array}\right.$

对于分装在卫星三轴上的反作用轮，卫星的控制力矩指令可以表示为：

${\mathrm{\Σ}}_2:\left\{\begin{array}{c}{J}_1{\stackrel{\·}{\mathrm{\Ω}}}_1\left(t\right)=k_{T1}{i}_1\left(t\right)+M_{f1}\left(t\right)\\ J_2{\stackrel{\·}{\mathrm{\Ω}}}_2\left(t\right)=K_{T2}{i}_2\left(t\right)+M_{f2}\left(t\right)\\ J_3{\stackrel{\·}{\mathrm{\Ω}}}_3\left(t\right)=K_{T3}{i}_3\left(t\right)+M_{f3}\left(t\right)\end{array}\right.$

其中，k_T1,k_T2,k_T3为安装在卫星三轴上的反作用轮中电机的电流力矩系数，其数值大小均为0.02，i₁(t),i₂(t),i₃(t)为安装在卫星三轴上的反作用轮中电机的电流，M_f1(t),M_f2(t),M_f3(t)为安装在卫星三轴上的反作用轮的摩擦力矩；进一步，三轴上安装的反作用轮中所用电机的电流方程为：

${\mathrm{\Σ}}_3:\left\{\begin{array}{c}{i}_1\left(t\right)=\frac{1}{R_1}[E_1\left(t\right)-k_{E1}({\mathrm{\Ω}}_1\left(t\right)-\stackrel{\·}{\mathrm{\φ}}\left(t\right))]\\ i_2\left(t\right)=\frac{1}{R_2}[E_2\left(t\right)-k_{E2}({\mathrm{\Ω}}_2\left(t\right)-\stackrel{\·}{\mathrm{\θ}}\left(t\right))]\\ i_3\left(t\right)=\frac{1}{R_3}[E_3\left(t\right)-k_{E3}({\mathrm{\Ω}}_3\left(t\right)-\stackrel{\·}{\mathrm{\ψ}}\left(t\right))]\end{array}\right.$

其中，R₁,R₂,R₃为安装在卫星三轴上的反作用轮电机电阻，其数值大小均为1，E₁(t),E₂(t),E₃(t)为安装在卫星三轴上的三个反作用轮的控制电压；k_E1,k_E2,k_E3是安装在卫星三轴上的三个反作用轮中电机的反电动势系数，其数值大小均为0.02；

进一步，建立包含反作用轮摩擦特性的卫星姿态系统动力学模型为：

${\mathrm{\Σ}}_4:\left\{\begin{array}{c}{I}_1\stackrel{\·\·}{\mathrm{\φ}}\left(t\right)-n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\ψ}}\left(t\right)-\frac{k_{T1}*k_{E1}}{R_1}\stackrel{\·}{\mathrm{\φ}}\left(t\right)+4n^2(I_2-I_3)\mathrm{\φ}\left(t\right)\\ =\frac{k_{T1}}{R_1}{E}_1\left(t\right)-\frac{k_{T1}*k_{E1}}{R_1}*{\mathrm{\Ω}}_1\left(t\right)+M_{f1}\left(t\right)+T_{d1}\left(t\right)\\ I_2\stackrel{\·\·}{\mathrm{\θ}}\left(t\right)+3n^2(I_1-I_3)\mathrm{\θ}\left(t\right)-\frac{k_{T2}*k_{E2}}{R_2}\stackrel{\·}{\mathrm{\θ}}\left(t\right)\\ =\frac{k_{T2}}{R_2}{E}_2\left(t\right)-\frac{k_{T2}*k_{E2}}{R_2}*{\mathrm{\Ω}}_2\left(t\right)+M_{f2}\left(t\right)+T_{d2}\left(t\right)\\ I_3\stackrel{\·\·}{\mathrm{\ψ}}\left(t\right)+n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\φ}}\left(t\right)+n^2(I_2-I_1)\mathrm{\ψ}\left(t\right)-\frac{k_{T3}*k_{E3}}{R_3}\stackrel{\·}{\mathrm{\ψ}}\left(t\right)\\ =\frac{k_{T3}}{R_3}{E}_3\left(t\right)-\frac{k_{T3}*k_{E3}}{R_3}*{\mathrm{\Ω}}_3\left(t\right)+M_{f3}\left(t\right)+T_{d3}\left(t\right)\end{array}\right.$

将式∑₄写成矩阵形式如下：

$P\stackrel{\·\·}{p}\left(t\right)+C\stackrel{\·}{p}\left(t\right)+\mathrm{Qp}\left(t\right)=B_{E}E\left(t\right)+B_{\mathrm{\Ω}}\mathrm{\Ω}\left(t\right)+B_m{M}_f\left(t\right)+B_w{T}_d\left(t\right)$

其中，状态p(t)＝[φ(t),θ(t),ψ(t)]^T，T_d(t)＝[T_d1(t),T_d2(t),T_d3(t)]^T为卫星三轴受到的环境干扰力矩，E(t)＝[E₁(t),E₂(t),E₃(t)]^T为安装在卫星三轴上的三个反作用轮的控制电压，由控制器解算所得，Ω(t)＝[Ω₁(t),Ω₂(t),Ω₃(t)]^T为安装在卫星三轴上的反作用轮中电机的角速度，M_f(t)＝[M_f1,M_f2,M_f3]^T为安装在卫星三轴上的反作用轮的摩擦力矩干扰；式中P＝diag{I₁,I₂,I₃},B_m＝I,B_w＝I，Q＝diag{4n²(I₂-I₃),3n²(I₁-I₃),n²(I₂-I₁)}, $B_{\mathrm{\Ω}}=\mathrm{diag}\{-\frac{k_{T1}*k_{E1}}{R_1},-\frac{k_{T2}*k_{E2}}{R_2},-\frac{k_{T3}*k_{E3}}{R_3}\},$ $C=\left[\begin{array}{ccc}-\frac{k_{T1}*k_{E1}}{R_1}& 0& -n(I_1-I_2+I_3)\\ 0& -\frac{k_{T2}*k_{E2}}{R_2}& 0\\ n(I_1-I_2+I_3)& 0& -\frac{k_{T3}*k_{E3}}{R_3}\end{array}\right];$

进一步可以将包含反作用轮摩擦特性的卫星姿态系统动力学模型写成状态空间表达式形式如下：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+B_{0}E\left(t\right)+B_1{M}_f\left(t\right)+B_2{T}_d\left(t\right)$

其中，x(t)、E(t)分别为系统的状态变量和控制输入，E(t)＝[E₁(t),E₁(t),E₁(t)]^T $x\left(t\right)={\left[\begin{array}{cccc}{e}_p\left(t\right)& {\stackrel{\·}{e}}_p\left(t\right)& {\∫}_{t-1}^t\mathrm{\Ω}\left(\mathrm{\τ}\right)d\left(t\right)& \mathrm{\Ω}\left(t\right)\end{array}\right]}^T,$ 其中e_p(t)＝p(t)-p_d(t)，p_d(t)为航天器系统的期望姿态信息，p(t)为航天器系统当前最新的姿态信息，e_p(t)为航天器当前姿态与期望姿态的偏差；A为系统阵，B₀为控制输入矩阵，B₁为摩擦力矩干扰输入矩阵，B₂为环境干扰力矩的输入矩阵，其对应的表达式如下：

$A=\left[\begin{array}{cccc}0& 0& I& 0\\ 0& 0& 0& I\\ -P^{-1}Q& 0& -P^{-1}C& P^{-1}{B}_{\mathrm{\Ω}}\\ 0& 0& \frac{K_T*K_E}{J*R}& -\frac{K_T*K_E}{J*R}\end{array}\right],B_0=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_E\\ \frac{K_T}{J*R}\end{array}\right],B_1=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_m\\ \frac{1}{J}\end{array}\right],B_2=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_w\\ 0\end{array}\right];$

其中，K_T＝diag{k_T1,k_T2,k_T3}，J＝diag{J₁,J₂,J₃}，K_E＝diag{k_e1,k_e2,k_e3}，R＝diag{R₁,R₂,R₃}；

2、设计摩擦干扰估计器

对于上述反作用轮摩擦干扰，设计摩擦干扰估计器，对摩擦干扰的估计表征为：

$\left\{\begin{array}{c}{\hat{M}}_f\left(t\right)=r\left(t\right)+\mathrm{Nx}\left(t\right)\\ \stackrel{\·}{r}\left(t\right)=-N{B}_1(r\left(t\right)+\mathrm{Nx}\left(t\right))-N(\mathrm{Ax}\left(t\right)+B_0\left(t\right))\end{array}\right.$

其中，N为摩擦干扰估计器的增益，在此实施例中，以卫星俯仰轴姿态控制为例，N的取值大小为[0 0.0008 0 0.1367]，为反作用轮摩擦干扰的估计值，r(t)为摩擦干扰估计器中的辅助变量；

3、设计标称控制器

在通过摩擦干扰估计器完成对反作用轮摩擦进行估计后，对于卫星系统中还存在的环境干扰力矩以及摩擦干扰估计误差干扰，进一步设计具有干扰抑制能力的标称控制器，在此实施例中，选用的标称控制律为PD控制律，其具体表现形式为：

$E_K\left(t\right)=\mathrm{Kx}\left(t\right)=K_P\left[\begin{array}{c}{e}_p\left(t\right)\\ {\∫}_{t-1}^t\mathrm{\Ω}\left(\mathrm{\τ}\right)d\left(t\right)\end{array}\right]+K_D\left[\begin{array}{c}{\stackrel{\·}{e}}_p\left(t\right)\\ \mathrm{\Ω}\left(t\right)\end{array}\right]$

其中，E_K(t)为当前时刻PD控制器的控制量，K为PD控制器的控制增益,以卫星俯仰轴姿态控制为例，PD控制器的控制增益为K＝[0.2160 2.3540 -0.0438 -0.2467]，K_P为PD控制器比例增益，K_P取值为[0.2160 -0.0438]，K_D为PD控制器微分增益，K_D取值为[2.3540 -0.2467]；

4、设计卫星复合分层抗干扰姿态控制器

在分别完成反作用轮摩擦观测器与标称控制器基础上，对两种控制器进行复合，得到多源干扰下卫星复合分层抗干扰姿态控制器，其对应的控制律为：

$E\left(t\right)=E_K\left(t\right)-\frac{R}{K_T}{\hat{M}}_f\left(t\right)=E_K\left(t\right)-50{\hat{M}}_f\left(t\right)$

其中，E_K(t)为复合分层抗干扰姿态控制器的输出，用于控制安装在微纳卫星三轴上的反作用轮完成姿态控制。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (1)

1.一种基于反作用轮摩擦特性的抗干扰姿态控制方法，其特征在于，包括以下步骤：针对反作用轮摩擦以及姿控精度较低的问题，首先，构建含有反作用轮摩擦特性的航天器耦合动力学模型，其次，针对由于反作用轮带来的摩擦干扰设计摩擦干扰估计器，对反作用轮摩擦进行实时估计；再次，针对航天器系统中存在的摩擦干扰估计误差及环境干扰力矩，设计标称控制器进行抑制；最后，通过联合标称控制器与摩擦干扰估计器，设计复合分层抗干扰控制器，完成多源干扰影响下的航天器抗干扰姿态控制；

第一步：构建含有反作用轮摩擦特性的航天器耦合动力学模型

当航天器本体坐标系和轨道坐标系之间的欧拉角很小时,姿态运动学可以小角度线性化,得到如下的姿态动力学模型：

${\mathrm{\Σ}}_1:\left\{\begin{array}{c}{I}_1\stackrel{\·\·}{\mathrm{\φ}}\left(t\right)-n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\ψ}}\left(t\right)+4n^2(I_2-I_3)\mathrm{\φ}\left(t\right)=J_1{\stackrel{\·}{\mathrm{\Ω}}}_1\left(t\right)+T_{d1}\left(t\right)\\ I_2\stackrel{\·\·}{\mathrm{\θ}}\left(t\right)+3n^2(I_1-I_3)\mathrm{\θ}\left(t\right)=J_2{\stackrel{\·}{\mathrm{\Ω}}}_2\left(t\right)+T_{d2}\left(t\right)\\ I_3\stackrel{\·\·}{\mathrm{\ψ}}\left(t\right)+n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\φ}}\left(t\right)+n^2(I_2-I_1)\mathrm{\ψ}\left(t\right)=J_3{\stackrel{\·}{\mathrm{\Ω}}}_3\left(t\right)+T_{d3}\left(t\right)\end{array}\right.$

其中，I₁,I₂,I₃分别为航天器的三轴转动惯量；φ(t),θ(t),ψ(t)分别为三轴姿态角，分别为三轴姿态角速度，分别为三轴姿态角加速度；J₁,J₂,J₃分别为三轴上安装的反作用轮的转动惯量；分别表示航天器三轴上的反作用轮中电机的转动角加速度，T_d1(t),T_d2(t),T_d3(t)分别为航天器三轴受到的环境干扰力矩；n为航天器轨道角速度；

用于航天器姿态控制的反作用轮系统一般是由微型直流无刷电机、反作用轮本体、以及控制器线路构成，是一个相对独立的机电执行部件,反作用轮的转动方程为：

$\left\{\begin{array}{c}{J}_i{\stackrel{\·}{\mathrm{\Ω}}}_i\left(t\right)=k_{T_i}*i_i\left(t\right)+M_{\mathrm{fi}}\left(t\right)\\ E_i\left(t\right)=k_{\mathrm{Ei}}*{\mathrm{\Ω}}_i\left(t\right)+R_i*i_i\left(t\right)+L_i*{\stackrel{\·}{i}}_i\left(t\right)\end{array}\right.$

其中，J_i为反作用轮的转动惯量,Ω_i(t)为反作用轮中电机的转动角速度，为反作用轮中电机的转动角加速度，k_Ti为反作用轮中电机的电流力矩系数，i_i(t)为反作用轮中电机的电流，k_Ti*i_i(t)为反作用轮的驱动力矩，M_fi(t)为反作用轮摩擦力矩；对于反作用轮使用的无刷直流电机，电流力矩系数k_Ti在数值上与反作用轮中的反电动势系数相等,k_Ei为反作用轮使用电机的反电动势系数，E_i(t)为加在反作用轮两端的控制电压，R_i为反作用轮使用电机的电阻值，L_i为反作用轮使用电机的电感，在正常使用过程中，在数值上大小接近0，在实际处理过程中，这一项可以忽略，因此反作用轮转动方程进一步可以描述为；

$\left\{\begin{array}{c}{J}_i{\stackrel{\·}{\mathrm{\Ω}}}_i\left(t\right)=k_{\mathrm{Ti}}*i_i\left(t\right)+M_{\mathrm{fi}}\left(t\right)\\ E_i\left(t\right)=k_{\mathrm{Ei}}*{\mathrm{\Ω}}_i\left(t\right)+R_i*i_i\left(t\right)\end{array}\right.$

对于分装在航天器三轴上的反作用轮，航天器的控制力矩指令可以表示为：

${\mathrm{\Σ}}_2:\left\{\begin{array}{c}{J}_1{\stackrel{\·}{\mathrm{\Ω}}}_1\left(t\right)=k_{T1}{i}_1\left(t\right)+M_{f1}\left(t\right)\\ J_2{\stackrel{\·}{\mathrm{\Ω}}}_2\left(t\right)=k_{T2}{i}_2\left(t\right)+M_{f2}\left(t\right)\\ J_3{\stackrel{\·}{\mathrm{\Ω}}}_3\left(t\right)=K_{T3}{i}_3\left(t\right)+M_{f3}\left(t\right)\end{array}\right.$

其中，k_T1,k_T2,k_T3为安装在航天器三轴上的反作用轮中电机的电流力矩系数，i₁(t),i₂(t),i₃(t)为安装在航天器三轴上的反作用轮中电机的电流，M_f1(t),M_f2(t),M_f3(t)为安装在航天器三轴上的反作用轮的摩擦力矩；进一步，三轴上安装的反作用轮中所用电机的电流方程为：

${\mathrm{\Σ}}_3:\left\{\begin{array}{c}{i}_1\left(t\right)=\frac{1}{R_1}[E_1\left(t\right)-k_{E1}({\mathrm{\Ω}}_1\left(t\right)-\stackrel{\·}{\mathrm{\φ}}\left(t\right))]\\ i_2\left(t\right)=\frac{1}{R_2}[E_2\left(t\right)-k_{E2}({\mathrm{\Ω}}_2\left(t\right)-\stackrel{\·}{\mathrm{\θ}}\left(t\right))]\\ i_3\left(t\right)=\frac{1}{R_3}[E_3\left(t\right)-k_{E3}({\mathrm{\Ω}}_3\left(t\right)-\stackrel{\·}{\mathrm{\ψ}}\left(t\right))]\end{array}\right.$

其中，R₁,R₂,R₃为安装在航天器三轴上的反作用轮电机电阻，E₁(t),E₂(t),E₃(t)为安装在航天器三轴上的三个反作用轮的控制电压；k_E1,k_E2,k_E3是安装在航天器三轴上的三个反作用轮中电机的反电动势系数；

进一步，建立包含反作用轮摩擦特性的航天器姿态系统动力学模型为：

${\mathrm{\Σ}}_4:\left\{\begin{array}{c}{I}_1\stackrel{\·\·}{\mathrm{\φ}}\left(t\right)-n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\ψ}}\left(t\right)-\frac{k_{T1}*k_{E1}}{R_1}\stackrel{\·}{\mathrm{\φ}}\left(t\right)+{4n}^2(I_2-I_3)\mathrm{\φ}\left(t\right)\\ =\frac{k_{T1}}{R_1}{E}_1\left(t\right)-\frac{k_{T1}*k_{E1}}{R_1}*{\mathrm{\Ω}}_1\left(t\right)+M_{f1}\left(t\right)+T_{d1}\left(t\right)\\ I_2\stackrel{\·\·}{\mathrm{\θ}}\left(t\right)+{3n}^2(I_1-I_3)\mathrm{\θ}\left(t\right)-\frac{k_{T2}*k_{E2}}{R_2}\stackrel{\·}{\mathrm{\θ}}\left(t\right)\\ =\frac{k_{T2}}{R_2}{E}_2\left(t\right)-\frac{k_{T2}*k_{E2}}{R_2}*{\mathrm{\Ω}}_2\left(t\right)+M_{f2}\left(t\right)+T_{d2}\left(t\right)\\ I_3\stackrel{\·\·}{\mathrm{\ψ}}\left(t\right)+n(I_1-I_2+I_3)\stackrel{\·}{\mathrm{\φ}}\left(t\right)+n^2(I_2-I_1)\mathrm{\ψ}\left(t\right)-\frac{k_{T3}*k_{E3}}{R_3}\stackrel{\·}{\mathrm{\ψ}}\left(t\right)\\ =\frac{k_{T3}}{R_3}{E}_3\left(t\right)-\frac{k_{T3}*k_{E3}}{R_3}*{\mathrm{\Ω}}_3\left(t\right)+M_{f3}+T_{d3}\left(t\right)\end{array}\right.$

将式Σ₄写成矩阵形式如下：

$P\stackrel{\·\·}{p}\left(t\right)+C\stackrel{\·}{p}\left(t\right)+\mathrm{Qp}\left(t\right)=B_{E}E\left(t\right)+B_{\mathrm{\Ω}}\mathrm{\Ω}\left(t\right)+B_m{M}_f\left(t\right)+B_w{T}_d\left(t\right)$

其中，状态p(t)＝[φ(t),θ(t),ψ(t)]^T，T_d(t)＝[T_d1(t),T_d2(t),T_d3(t)]^T为航天器三轴受到的环境的干扰力矩，E(t)＝[E₁(t),E₂(t),E₃(t)]^T为安装在航天器三轴上的三个反作用轮的控制电压，Ω(t)＝[Ω₁(t),Ω₂(t),Ω₃(t)]^T为安装在航天器三轴上的反作用轮中电机的角速度，M_f(t)＝[M_f1(t),M_f2(t),M_f3(t)]^T为安装在航天器三轴上的反作用轮的摩擦力矩干扰；式中P＝diag{I₁,I₂,I₃}，B_m＝I,B_w＝I， $B_E=\mathrm{diag}\{\frac{k_{T1}}{R_1},\frac{k_{T2}}{R_2},\frac{k_{T3}}{R_3}\},$ Q＝diag{4n²(I₂-I₃),3n²(I₁-I₃),n²(I₂-I₁)}, $B_{\mathrm{\Ω}}=\mathrm{diag}\{-\frac{k_{T1}*k_{E1}}{R_1},-\frac{k_{T2}*k_{E2}}{R_2},-\frac{k_{T3}*k_{E3}}{R_3}\},$ $C=\left[\begin{array}{ccc}-\frac{k_{T1}*k_{E1}}{R_1}& 0& -n(I_1-I_2+I_3)\\ 0& -\frac{k_{T2}*k_{E2}}{R_2}& 0\\ n(I_1-I_2+I_3)& 0& -\frac{k_{T3}*k_{E3}}{R_3}\end{array}\right];$

进一步可以将包含反作用轮摩擦特性的航天器姿态系统动力学模型写成状态空间表达式形式如下：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+B_{0}E\left(t\right)+B_1{M}_f\left(t\right)+B_2{T}_d\left(t\right)$

其中，x(t)、E(t)分别为系统的状态变量和控制输入，E(t)＝[E₁(t),E₂(t),E₃(t)]^T， $x\left(t\right)={\left[\begin{array}{cccc}{e}_p\left(t\right)& {\∫}_{t-1}^t\mathrm{\Ω}\left(\mathrm{\τ}\right)d\left(t\right)& {\stackrel{.}{e}}_p\left(t\right)& \mathrm{\Ω}\left(t\right)\end{array}\right]}^T,$ 其中e_p(t)＝p(t)-p_d(t)，p_d(t)为航天器系统的期望姿态信息，p(t)为航天器系统当前最新的姿态信息，e_p(t)为航天器当前姿态与期望姿态的偏差，A为系统阵，B₀为控制输入矩阵，B₁为摩擦力矩干扰输入矩阵，B₂为环境干扰力矩的输入矩阵，其对应的表达式如下：

$A=\left[\begin{array}{cccc}0& 0& I& 0\\ 0& 0& 0& I\\ -P^{-1}Q& 0& -P^{-1}C& P^{-1}{B}_{\mathrm{\Ω}}\\ 0& 0& \frac{K_T*K_E}{J*R}& -\frac{K_T*K_E}{J*R}\end{array}\right],$ $B_0=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_E\\ \frac{K_T}{J*R}\end{array}\right],$ $B_1=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_m\\ \frac{1}{J}\end{array}\right],$ $B_2=\left[\begin{array}{c}0\\ 0\\ P^{-1}{B}_w\\ 0\end{array}\right];$

其中，K_T＝diag{k_T1,k_T2,k_T3}，J＝diag{J₁,J₂,J₃}，K_E＝diag{k_E1,k_E2,k_E3}，R＝diag{R₁,R₂,R₃}；

第二步：设计摩擦干扰估计器

对于上述反作用轮摩擦干扰，设计摩擦干扰估计器，对摩擦干扰的估计表征为：

$\left\{\begin{array}{c}{\hat{M}}_f\left(t\right)=r\left(t\right)+\mathrm{Nx}\left(t\right)\\ \stackrel{\·}{r}\left(t\right)=-{\mathrm{NB}}_1(r\left(t\right)+\mathrm{Nx}\left(t\right))-N(\mathrm{Ax}\left(t\right)+B_{0}u\left(t\right))\end{array}\right.$

其中，N为摩擦干扰估计器增益,为反作用轮摩擦干扰的估计值,r(t)为摩擦干扰估计器中的辅助变量；

第三步：设计标称控制器

通过摩擦干扰估计器完成对反作用轮摩擦进行估计后，对于航天器系统中还存在的环境干扰力矩以及摩擦干扰估计误差，进一步设计具有干扰抑制能力的标称控制器，其对应控制律为：

$E_K\left(t\right)=K{\left[\begin{array}{cccc}{e}_p\left(t\right)& {\∫}_{t-1}^t\mathrm{\Ω}\left(\mathrm{\τ}\right)d\left(t\right)& {\stackrel{\·}{e}}_p\left(t\right)& \mathrm{\Ω}\left(t\right)\end{array}\right]}^T$

其中，E_K(t)为当前时刻标称控制器的控制量，K为标称控制器的控制增益；

第四步：设计航天器复合分层抗干扰姿态控制器

在分别完成反作用轮摩擦观测器与标称控制器基础上，进行复合，得到多源干扰下航天器复合分层抗干扰姿态控制器，其对应的控制律为：

$E\left(t\right)=E_K\left(t\right)-\frac{R}{K_T}{\hat{M}}_f\left(t\right)$

其中，E(t)为复合分层抗干扰姿态控制器的输出，用于控制安装在航天器三轴上的反作用轮完成姿态控制。